Topology optimization of DC-biased pot-type reactor using design sensitivity in steady state of electromagnetic field with magnetic nonlinearity

Abstract

Because reactors applied to electric vehicles are typically driven under a DC-biased current, a constant inductance property is required. When a high-performance DC-biased reactor is designed, the topology optimization (TO) is an effective design method owing to the high degree of structural change in the magnetic circuit. However, there exists no report regarding the TO method in the steady-state time domain with magnetic nonlinearity. In this paper, based on the time-domain adjoint variable method, a novel sensitivity analysis approach for the steady-state time domain with thorough consideration of the magnetic saturation is proposed. Through the proposed method, the TO of the iron core and winding structures in a pot-type reactor was carried out to enhance the DC-biased characteristics. Furthermore, the eddy current loss occurring on the aluminum plate installed outside the reactor was suppressed by the multi-objective TO. The interesting structure of a pot-type reactor that enhanced the DC-biased characteristics and reduced the eddy current loss on the outer conductor plate was also illustrated.

Keywords

DC-biased inductance design sensitivity in steady-state pot-type reactor time domain adjoint variable method topology optimization

1. Introduction

Reactors utilized in electrical circuits installed in electrical vehicles have to be downsized and made adaptable to large currents [1]. In DC-biased reactors for such currents, inductance deterioration owing to the magnetic saturation in the iron core is a serious problem because the behavior of the electrical circuit is unstable. Furthermore, the leakage magnetic field from the reactor causes malfunction of the electric circuits and the electromagnetic environment inside the xEV. Therefore, the design of a high-performance reactor with a stable inductance property and a lower leakage magnetic flux to the surrounding area is required.

To satisfy these design requirements, the topology optimization (TO) method [2], which has a high degree of structural change in the magnetic circuit, is applied to the design method. A sensitivity-based TO [3,4], which combines the electromagnetic field analysis based on the finite element method (FEM) with design sensitivity analysis utilizing the adjoint variable method (AVM) [5], is widely utilized from the viewpoints of fast convergence and derivation of a reasonable structure.

When the DC-biased characteristics of the reactor are evaluated under actual driving conditions, it is necessary to consider the magnetic saturation of the iron core and to evaluate some physical quantities in the steady state. Although frequency-domain analysis [6] can directly introduce a steady-state solution, it cannot realize an electromagnetic field that includes higher harmonic components derived from magnetic nonlinearity. Therefore, a time-domain analysis that considers the magnetic saturation of the iron core is essential. When the design sensitivity with magnetic nonlinearity is carried out, the time-domain AVM [7,8] must be applied to TO. Subsequently, the inductance of the DC-biased reactor is evaluated through the steady state of the electromagnetic field. However, a time-domain AVM for steady-state quantity has not been proposed.

In this paper, a novel analysis approach of design sensitivity analysis for the steady state in the time domain with magnetic nonlinearity is proposed. The accuracy of the design sensitivity of the objective function derived from the proposed method is verified, and the TO method using the steady-state sensitivity is applied to pot-type reactors [1]. The design goal of a pot-type reactor is to minimize the eddy current loss on the outer shielding plate induced by the leakage magnetic flux from the reactor and maximize the self-inductance with DC-biased AC current. Hybrid objectives are achieved by optimizing the topology of the magnetic core and the winding distribution.

2. Topology optimization method for steady-state electromagnetic field

2.1. Topology modeling

In this study, the design domains of the iron core and the winding were explicitly separated. The design variables of both domains are defined as follows: $\begin{eqnarray}\{{\psi}\}=\{{\psi}_{1},{\psi}_{2}\}^{T},\end{eqnarray}$ (1) where 𝜓₁ is the design variable defined in the design domain for the iron core, and 𝜓₂ is the design variable in the design domain for the winding. Next, the continuous Heaviside function H(𝜓) changing within the range 0 ≤ H(𝜓) ≤ 1 [9] is incorporated with the material properties. Here, H(𝜓) is given as: $\begin{eqnarray}H({\psi})=\frac{3}{16}\left(\frac{{\psi}}{h}\right)^{5}-\frac{5}{8}\left(\frac{{\psi}}{h}\right)^{3}+\frac{15}{16}\left(\frac{{\psi}}{h}\right)+\frac{1}{2}\quad (-h\leq {\psi}\leq h),\end{eqnarray}$ (2) where h is the transition width. When the continuous Heaviside function is used, it is easy to binarize the material distribution into ON or OFF state. Next, the magnetic reluctivity 𝜈(𝜓₁) in the design domain of the iron core structure is defined as follows: $\begin{eqnarray}v({\psi}_{1})=\{1-H({\psi}_{1})\}v_{0}+H({\psi}_{1})v_{e}(|\boldsymbol{B}|),\end{eqnarray}$ (3) where B is the magnetic flux density, 𝜈₀ is the magnetic reluctivity in a vacuum, and 𝜈_e(| B |) is the nonlinear reluctivity of the magnetic material. The source current density J(𝜓₂, t) at time t in the design domain of the winding structure is defined as follows: $\begin{eqnarray}J({\psi}_{2},t)=H({\psi}_{2})J_{0}(t),\end{eqnarray}$ (4) where J₀(t) is the source current density, which is the input to the winding at time t.

Fig. 1.

Example of objective function in steady state.

2.2. Sensitivity analysis in steady state of time domain

The objective function in steady-state is defined as follows: $\begin{eqnarray}W=\int _{t_{0}}^{t_{0}+T}w({\psi},\boldsymbol{A},\dot{\boldsymbol{A}})dt,\end{eqnarray}$ (5) where t₀ is the starting time of the steady state, T represents a single time period, A and $\dot{\boldsymbol{A}}$ indicate the magnetic vector potential and its time derivative ∂ A ∕∂t, respectively, while w is the integrand function, for example, eddy current loss, etc. To efficiently evaluate the design sensitivity ∂W∕∂𝜓, the time-domain AVM was applied. However, because conventional time-domain AVM uses the time integral between the initial and terminal time, the design sensitivity based on the specific integral interval cannot be calculated normally. To calculate ∂W∕∂𝜓 in a steady state, the objective function W₁ in the interval [0, t₀] and W₂ evaluated in the interval [0, t₀ + T] are defined as follows: $\begin{eqnarray}\displaystyle W_{1} & = & \displaystyle \int _{0}^{t_{0}}w({\psi},\boldsymbol{A},\dot{\boldsymbol{A}})dt,\end{eqnarray}$ (6) $\begin{eqnarray}\displaystyle W_{2} & = & \displaystyle \int _{0}^{t_{0}+T}w({\psi},\boldsymbol{A},\dot{\boldsymbol{A}})dt,\end{eqnarray}$ (7) where examples of W₁ and W₂ are shown in Fig. 1. The objective function W is evaluated by the subtracting the transient component W₁ from the total component W₂. Therefore, the design sensitivity of W is obtained as follows: $\begin{eqnarray}\frac{\partial W}{\partial {\psi}}=\frac{\partial W_{2}}{\partial {\psi}}-\frac{\partial W_{1}}{\partial {\psi}}.\end{eqnarray}$ (8) To evaluate ∂W₁∕∂𝜓 and ∂W₂∕∂𝜓, the augmented functions of W₁ and W₂ are defined as follows: $\begin{eqnarray}\displaystyle \overline{W}_{1} & = & \displaystyle \int _{0}^{t_{0}}\{w({\psi},\boldsymbol{A},\dot{\boldsymbol{A}})+\boldsymbol{\boldsymbol{\it\lambda}}_{1}(t)^{T}(K\boldsymbol{A}+C\dot{\boldsymbol{A}}-\boldsymbol{F})\}dt,\end{eqnarray}$ (9) $\begin{eqnarray}\displaystyle \overline{W}_{2} & = & \displaystyle \int _{0}^{t_{0}+T}\{w({\psi},\boldsymbol{A},\dot{\boldsymbol{A}})+\boldsymbol{\boldsymbol{\it\lambda}}_{2}(t)^{T}(K\boldsymbol{A}+C\dot{\boldsymbol{A}}-\boldsymbol{F})\}dt,\end{eqnarray}$ (10) where the vectors 𝝀₁ and 𝝀₂ are the adjoint variables for W₁ and W₂, respectively. The matrices K, C, and the vector F are the reluctivity-dependent stiffness matrix, conductivity-dependent matrix and J₀-dependent vector, respectively. The time-domain adjoint equation with magnetic nonlinearity is described as follows; $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle C^{T}\boldsymbol{\boldsymbol{\it\lambda}}_{l}(t_{l})=\left.-\frac{\partial w}{\partial \dot{\boldsymbol{A}}}\right|_{t=t_{l}},\end{eqnarray}$ (11) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left(K+\frac{\partial K}{\partial \boldsymbol{A}}\boldsymbol{A}\right)^{T}\boldsymbol{\boldsymbol{\it\lambda}}_{l}(t)-C^{T}\dot{\boldsymbol{\boldsymbol{\it\lambda}}}_{l}(t)=\frac{\partial }{\partial t}\left(\frac{\partial w}{\partial \dot{\boldsymbol{A}}}\right)-\frac{\partial w}{\partial \boldsymbol{A}},\end{eqnarray}$ (12) where 𝝀_l is the adjoint variable for the objective function W_l, in which the value 1 or 2 is substituted for l, and t_l is the terminal time. If l = 1, t₀ is set to t_l. If l = 2, (t₀ + T) is set to t_l. The terminal value of 𝝀_l(t) is primarily determined by solving (11), and 𝝀_l(t) is sequentially calculated in the direction of the initial time from the terminal time using the forward Euler method. The sensitivity ∂W_l∕∂𝜓 is derived from the substitution of 𝝀_l(t) in the following equation: $\begin{eqnarray}\frac{\partial W_{l}}{\partial {\psi}}=\int _{0}^{t_{l}}\left\{\frac{\partial w}{\partial {\psi}}+\boldsymbol{\boldsymbol{\it\lambda}}_{l}(t)^{T}\left(\frac{\partial K}{\partial {\psi}}\boldsymbol{A}+\frac{\partial C}{\partial {\psi}}\dot{\boldsymbol{A}}-\frac{\partial \boldsymbol{F}}{\partial {\psi}}\right)\right\}dt.\end{eqnarray}$ (13) Next, by substituting both derivatives ∂W₁∕∂𝜓 and ∂W₂∕∂𝜓 for (8), the steady-state sensitivity ∂W∕∂𝜓 is derived.

Fig. 2.

Flowchart of sensitivity analysis in steady-state.

The procedure for the AVM-based sensitivity analysis in the steady state is summarized in Fig. 2. The transient nonlinear eddy current analysis is carried out using a step-by-step procedure based on the backward Euler method, with both A and $\dot{\boldsymbol{A}}$ at all time steps being stored (step 1). Next, by solving the adjoint equations shown in ((11)) and ((12)) using the step-by-step procedure based on the forward Euler method, the adjoint variable at all time steps is obtained (step 2). Subsequently, ∂W₁∕∂𝜓 is evaluated with ((13)) (step 3). Because ∂W₂∕∂𝜓 is also calculated using the same procedure, ((8)) can derive the steady-state sensitivity ∂W∕∂𝜓 (step 4).

Fig. 3.

Analysis model of pot-type reactor.

3. Optimization problem

A 2-D axisymmetric pot-type reactor model is illustrated in Fig. 3, where the design domains are set to Ω_m and Ω_e. While the structure of the magnetic core is optimized in Ω_m, the winding structure is optimized in Ω_e. The number of elements and nodes were set to 5,623 and 2,895, respectively. The optimization objective is to minimize the eddy current loss W_e induced on the outer conductive plate and to maximize the square of the DC-biased self-inductance W_L. Both W_e and W_L were calculated in the steady state (t₀ ≤ t ≤ t₀ + T). The optimization problem is formulated as follows: $\begin{eqnarray}\begin{array}{@{}ll@{}}\text{min}. & W({\psi}_{1},{\psi}_{2})=c_{e}W_{e}({\psi}_{1},{\psi}_{2})-c_{L}W_{L}({\psi}_{1},{\psi}_{2})\\[3.0pt] \text{s.t.} & g_{1}({\psi}_{1})=V_{1}({\psi}_{1})-{\alpha}_{1}V_{{\rm\Omega}\text{m}}\leq 0\\[3.0pt] & g_{2}({\psi}_{2})=V_{2}({\psi}_{2})-{\alpha}_{2}V_{{\rm\Omega}\text{w}}\leq 0\\[3.0pt] & -h\leq {\psi}_{1}\leq h\\[3.0pt] & -h\leq {\psi}_{2}\leq h,\end{array}\end{eqnarray}$ (14) where c_e and c_L represent the weighted coefficients that adjust the degree of effect of W_e and W_L on W. The function V₁(𝜓₁) is the volume of the iron core in the design domain for the iron core distribution. The soft magnetic composite EU67 is set as the nonlinear magnetic material applied to Ω_m. The constant V_{Ω_m} is the whole volume of the design domain: the constant 𝛼₁ indicates the ratio of V₁(𝜓₁) to V_{Ω_m}, while the function V₂(𝜓₁) stands for the winding volume in the design domain for the winding distribution; the constant V_{Ω_w} means the entire volume of the design domain for the winding distribution, whereas constant 𝛼₂ implies the ratio of V₂(𝜓₂) to V_{Ω_w}. The functions W_e, W_L, V₁(𝜓₁), and V₂(𝜓₂) are defined as follows. $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle W_{e}=\mathop{\sum }_{i=1}^{N}\left\{\int _{t_{0}}^{t_{0}+T}\left(2{\rm\pi}\iint _{{\rm\Omega}_{\text{c}}}r\frac{{\boldsymbol{J}_{e}}^{(i)^{2}}}{{\sigma}}drdz\right)dt\right\},\end{eqnarray}$ (15) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle W_{L}=\mathop{\sum }_{i=1}^{N}L^{(i)^{2}}=\mathop{\sum }_{i=1}^{N}\left\{\left(\frac{1}{{\omega}{I_{rms}}^{(i)}T}\right)^{2}\int _{t_{0}}^{t_{0}+T}\left(\frac{d{\Phi}^{(i)}}{dt}\right)^{2}dt\right\}^{2},\end{eqnarray}$ (16) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle V_{1}({\psi}_{1})=2{\rm\pi}\iint _{{\rm\Omega}_{\text{m}}}rH({\psi}_{1})drdz,\end{eqnarray}$ (17) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle V_{2}({\psi}_{2})=2{\rm\pi}\iint _{{\rm\Omega}_{\text{w}}}rH({\psi}_{2})drdz,\end{eqnarray}$ (18) where the variable i is the index that changes the value of DC current i_dc, as shown in Fig. 4, while N is the total number of i_dc used in the forward analysis based on the FEM, the vector $\boldsymbol{J}_{\boldsymbol{e}}^{(i)}$ is the eddy current density at the i-th i_dc, ω is the angular frequency, Φ⁽ⁱ⁾ is the linkage flux, and I_rms⁽ⁱ⁾ is the effective value at the i-th i_dc.

Fig. 4.

DC-biased AC current.

Table 1

Parameters for topology optimization

h (mm)	ϵ_FEM	ϵ_AVM	ϵ_NR	ϵ_opt
1	10⁻⁶	10⁻¹²	10⁻³	10⁻⁵

4. Optimization results

Here, to consider the constraints g₁ and g₂, the objective function W is combined with the constraints by the augmented Lagrangian method [10]. The steepest descent method is adopted as the optimization method.

The parameters for TO are illustrated in Table 1. The constants ϵ_FEM and ϵ_AVM represent the convergence criterion values based on the relative norm in a linear system derived from FEM and AVM, respectively. The incomplete Cholesky conjugate gradient method is applied to the solver for linear systems derived from FEM and AVM. The constant ϵ_NR is the convergence criterion for the Newton–Raphson method, in which the iteration is terminated when the maximum change of magnetic flux becomes a value less than ϵ_NR. The constant ϵ_opt is the convergence criterion value for the optimization method in which the maximum change in W becomes less than ϵ_opt.

4.1. Results of sensitivity analysis

To verify the accuracy of the sensitivity analysis derived from the proposed method, the sensitivities derived from the time-domain AVM were compared with the finite-difference method (FDM). The formula for the sensitivity-based FDM is as follows: $\begin{eqnarray}\frac{\partial W}{\partial {\psi}}\approx \frac{W({\psi}+{\rm\Delta}{\psi})-W({\psi})}{{\rm\Delta}{\psi}},\end{eqnarray}$ (19) where Δ𝜓 is infinitesimal with respect to the design variable. The frequency was set to 50 Hz, and the number of time steps per period (0.02 s) was set to 32. The constants T and t₀ were set to 0.02 s and 0.08 s, respectively. The design variables 𝛹₁ and 𝛹₂ for all nodes were set to 0.5, respectively.

The sensitivity distributions of ∂W_e∕∂𝜓 and ∂W_L∕∂𝜓 are illustrated in Fig. 5. It can be observed that the sensitivity in the design domain of the winding is larger than that of the magnetic core. The sensitivity distribution on the line at z = 7 mm is demonstrated in Fig. 6. The sensitivity distributions are quite similar between the AVM and FDM in Figs 6(a) and (b). The maximum sensitivity values (∂W_e∕∂𝜓)_max and (∂W_L∕∂𝜓)_max are also quite similar between AVM and FDM, as shown in Table 2. The elapsed time of the AVM shown in Table 2 was much less than that of the FDM.

Fig. 5.

Magnitudes of design sensitivity.

Fig. 6.

Sensitivity on the line at z = 7 mm.

4.2. Results of topology optimization

In this section, single-region TO for Ω_m and multi-region TO for both design domains Ω_m and Ω_e were carried out. The initial structures of a single-region and multi-region TO are shown in Fig. 7(a) and Fig. 7(b), respectively. Both initial structures converged to the best objective function among several trials of topology optimization. Table 3 shows the optimization conditions. Here, the DC-biased characteristics were evaluated by setting N to 6 in Eqs ((15)) and ((16)) with totally six different DC current values (i_dc = 50, 150, 250, 350, 450, and 550 A) are used.

Table 2
Comparison of maximum sensitivity

FDM AVM Relative error (%)

(∂W_e∕∂𝜓)_max × 10⁻³ 3.855 3.887 0.81

(∂W_L∕∂𝜓)_max × 10⁻¹¹ 4.453 4.534 1.8

Elapsed time (h) 27 (Estimated time) 0.014

	FDM	AVM	Relative error (%)
(∂W_e∕∂𝜓)_max × 10⁻³	3.855	3.887	0.81
(∂W_L∕∂𝜓)_max × 10⁻¹¹	4.453	4.534	1.8
Elapsed time (h)	27 (Estimated time)	0.014

Figure 8 shows the results derived from the single-region TO. Figure 8(a) shows the optimized magnetic core derived from the minimization of W_e without the maximizing the W_L. It can be seen that a one-layer shielding suppressing the leakage magnetic flux to the outside was generated. Figure 8(b) shows the optimized magnetic core derived from the maximization of W_L without the minimization of W_e. In this structure, the material was concentrated around the air gap to reduce its reluctivity. Figure 8(c) shows the magnetic core with both the minimization of W_e and the maximization of W_L. The resultant structure resembles the average shape shown in Figs 8(a) and (b). Although the gray scale area is negated, a reasonable core is obtained to satisfy the hybrid objectives.

Figure 9 shows the structures derived from the multi-material TO. The optimized structures of the iron core shown in Fig. 9 share similar portions to those shown in Fig. 8. The winding shown in Fig. 9(a) is located at a distance from the outer aluminum conductor. The winding shown in Fig. 9(b), which is optimized to increase the W_L, is located on the z-axis side. The winding shown in Fig. 9(c) considers both the minimization of W_e and the maximization of the W_L is generated at the bottom of Ω_W.

Figure 10(a) shows the DC-biased characteristics in the structure derived from the condition c_e = 0 or 1 with c_L = 0. When i_dc is increased, the inductances of all the structures deteriorate. The inductances of the optimized structures exhibit a higher performance than that of the initial model. In Fig. 10(b), when c_e is set to 1, the structure to reduce W_e is successfully obtained. Although the inductance properties of the optimized structures are deteriorated in comparison to those of the initial structure, the eddy current loss is reduced by TO. From this result, it can be deduced that the multi-region TO is more effective than the single-region TO.

Fig. 7.

Initial topologies.

Table 3

Optimization condition

f (Hz)	Number of times steps per 1 period	t_max (s)	t₀ (s)	i_ac (A)	Number of turns	𝛼₁	𝛼₂
50	16	0.08 (4T)	0.06 (3T)	30	10	0.6	0.4

Fig. 8.

Optimized structures of single-region TO (Ω_m).

Fig. 9.

Optimized structure derived from multi-region TO (Ω_m and Ω_w).

Fig. 10.

DC-biased characteristics.

As shown in Table 4, the elapsed time for all optimizations was several hours. The converged values of V₁(𝜓₁)∕V_{Ω_m} and V₂(𝜓₂)∕V_{Ω_w} largely satisfy the constraint conditions.

Table 4

Results of topology optimization

Design domain	c _e	c _L	W_e (nJ)	W_L × 10⁻³ (H² )	Elapsed time (h)	k _opt	V₁(𝜓₁)∕V_{Ω_m}	V₂(𝜓₂)∕V_{Ω_w}
	1	0	63.1	1.35	3.2 (A)	415	0.60
Ω_m	0	10	456	2.73	3.1 (B)	360	0.61	0.40 (const.)
	1	10	193	2.16	3.8 (A)	490	0.61
	1	0	50.5	0.992	1.3 (A)	229	0.57	0.40
Ω_m and Ω_w	0	10	496	2.82	4.3 (B)	500	0.57	0.40
	1	10	175	2.34	3.1 (A)	442	0.60	0.40
Initial model	-	-	319	2.58	-	-	0.60 (const.)	0.40 (const.)

A: Intel Core i7 4790 K at 4.0 GHz with 16.0 GB RAM. B: Intel Core i7 6850 K at 3.6 GHz with 128.0 GB RAM.

5. Conclusion

To realize TO in the steady state of the time domain with magnetic nonlinearity, a design sensitivity analysis procedure based on the time-domain adjoint variable method was introduced. The proposed approach can accurately derive the design sensitivity in the steady state of the time domain. Using the proposed method for its sensitivity, a reasonable structure for a high-performance pot-type reactor was obtained with the higher DC-biased inductance and the lower leakage flux.

Footnotes

Acknowledgement

All results were derived from Mr. Masuda who was a master’s student in my laboratory. For some reasons, he could not appear as a co-author of this manuscript. I would like to express sincere gratitude to him for his significant contribution.

References

Hashimoto

Zaitsu

Hayashi

Mitani

and Inoue

, Low loss reactor composed of dust core and a copperbelt coil, Kobe Steel Engineering Reports65(2) (2015).

Bendsøe

M.P.

, Optimal shape design as a material distribution problem, Struct. Optim.1(4) (1989), 193–202.

Dyck

D.N.

and Lowther

D.A.

, Automated design of magnetic devices by optimizing material distribution, IEEE Trans. Magn.32(3) (1996), 1188–1193.

Okamoto

Hoshino

Wakao

and Tsuburaya

, Improvement of torque characteristics for a synchronous reluctance motor using MMA-based topology optimization method, IEEE Trans. Magn.54(3) (2018), 7203104.

Gitosusatro

Coulomb

J.L.

and Sabonnadiere

J.C.

, Performance derivative calculation and optimization process, IEEE Trans. Magn.25(4) (1983), 2834–2839.

Shim

Wang

and Hameyer

, Topology optimization of magnetothermal systems considering eddy current as joule heat, IEEE Trans. Magn.43(4) (2007), 1617–1620.

Park

I.-H.

Kwak

I.-G.

Lee

H.-B.

Hahn

S.-Y.

and Lee

K.-S.

, Design sensitivity analysis for transient eddy current problems using finite element discretization and adjoint variable method, IEEE Trans. Magn.32(3) (1996), 1242–1245.

Masuda

Okamoto

and Wakao

, Sensitivity analysis based on time domain adjoint variable method in electromagnetic field problem, in: The Papers of Joint Technical Meeting on Static Apparatus, SA-18-20, RM-18-20, 2018, pp. 23–28.

Sethian

J.A.

and Wiegmann

, Structural boundary design via level set and immersed interface methods, J. Comp. Phys.163(2) (2000), 489–528.

10.

Suzuki

Wakabayashi

Kuwamura

Hirono

and Okuyama

, Application of constrained optimization using augmented Lagrange functions to photo attenuation correction in SPECT, Med. J. Imag. Technol.17(3) (1999), 245–260.