*3.3. Leakage Inductance of CW ELII MFT*

Figure 7 shows the simplified magnetic field distribution of CW ELII MFT. The integrated leakage inductance of CW ELII MFT consists of additional leakage inductance *L*σa, primary leakage inductance *L*σp, and secondary leakage inductance *L*σs. Since the leakage inductance of the primary winding and secondary winding are difficult to be separated physically, the two are usually combined as the winding leakage inductance *L*σw, so the total leakage inductance of the transformer is given as

$$L\_{\sigma} = L\_{\sigma \mathbf{a}} + L\_{\sigma \mathbf{p}} + L\_{\sigma \mathbf{s}}' = L\_{\sigma \mathbf{a}} + L\_{\sigma \mathbf{w}} \tag{18}$$

The Rogowski model and Dowell model are classical calculation methods for transformer leakage inductance. The Rogowski model is a 2D static model in which the internal characteristics of the winding (for example, conductor thickness and inter-turn distance) and the effect of eddy current in the conductor are ignored. The Dowell model is a 1D dynamic model, and describes the influence of the eddy current on the leakage inductance based on the assumption of parallel flux and uniform arrangement of conductors [21]. The Dowell model is suitable for cases with eddy current effects but has an obvious error in cases with low utilization of the core window height [22]. Researchers introduced Rogowski coefficient into the Dowell model, making up for the deficiency of Dowell model [22,23]. The hybrid model is given as [22]

$$L\_{\rm dW-H} = \frac{N^2 \mu\_0 l\_{\rm w}}{h\_{\rm eq}} \left[ \frac{d\_{\rm cs} m\_{\rm s}}{3} F\_{\rm 8W} + \frac{d\_{\rm cp} m\_{\rm P}}{3} F\_{\rm pw} + d\_{\rm d} + \frac{d\_{\rm is} (m\_{\rm s} - 1)(2m\_{\rm s} - 1)}{6m\_{\rm s}} + \frac{d\_{\rm ip} \left(m\_{\rm p} - 1\right)(2m\_{\rm p} - 1)}{6m\_{\rm p}} \right] \tag{19}$$

$$F\_{\rm W} = \frac{\left[ (4m^2 - 1) \frac{\sin \Lambda 2\Lambda - \sin 2\Lambda}{\cosh 2\Lambda - \cos 2\Lambda} - 2(m^2 - 1) \frac{\sinh \Lambda - \sin \Lambda}{\cosh \Lambda - \cos \Lambda} \right]}{2m^2 \Delta} \tag{20}$$

$$
\Delta = d\_{\rm c} \sqrt{\mu\_0 \pi f r \frac{nh\_{\rm c}}{h\_{\rm wd}}} \tag{21}
$$

$$h\_{\rm eq} = \frac{h\_{\rm w}}{K\_{\rm R}}\tag{22}$$

$$K\_{\rm R} = 1 - \frac{1 - e^{-\pi h\_{\rm w}/d\_{\rm w}}}{\pi h\_{\rm w}/d\_{\rm w}} \tag{23}$$

where *N* is the turns of the winding; *l*<sup>w</sup> is the average length of turns; *m*<sup>p</sup> and *m*<sup>s</sup> are respectively the numbers of primary and secondary winding layers; Δ is the ratio of conductor thickness to skin depth; *σ* is conductivity; *f* is frequency; *μ*<sup>0</sup> is the vacuum permeability; *n* is the number of conductors per layer; *K*<sup>R</sup> is Rogowski coefficient; and *h*eq is the corrected magnetic circuit length. As shown in Figure 8, *d*<sup>d</sup> is the main insulation distance; *d*wp, *d*ws, and *d*w, are respectively the width of the primary, secondary, and the entire winding; *h*<sup>w</sup> is the height of winding; *h*wd is the window height; *d*ip and *d*is are inter-layer distances; *d*cp and *d*cs are the widths of single conductors; *h*cp or *h*cs are the heights of conductors.

**Figure 7.** Simplified magnetic field distribution of CW ELII MFT.

**Figure 8.** Winding parameters.

According to the assumption of the Dowell model, the main insulation layer has the same width *d*<sup>d</sup> in the perimeter direction; conductor layers are arranged uniformly with the same inter-layer distance, i.e., *d*is and *d*ip. As shown in Figure 7, the width of the main insulation in the CW ELII MFT is not equal in the circumference direction. In addition, when there is an internal air channel in the winding, the inter-layer distance will not be equal. Therefore, the original model needs to be improved to enable it to be used in the CW ELII MFT.

Assuming that the cross-sectional area of the main insulation layer is *S*d, the inductance of the main insulation layer can be obtained based on the magnetic field energy

$$L\_{\rm d} = \frac{N^2 S\_{\rm d} \mu\_0}{h\_{\rm eq}} \tag{24}$$

Assuming that conductor layers are arranged non-uniformly with the inter-layer distances *d*i-1, *d*i-2,... , *d*i-*q*,... , *d*i-(*m*−1), then the inter-layer inductance is given as

$$L\_{\rm i} = \frac{\mu\_0 N^2 l\_{\rm w}}{h\_{\rm eq} m\_{\rm P}^2} \sum\_{q=1}^{q=m-1} q^2 d\_{\rm i-q} \tag{25}$$

Replace the static terms in (19) by (24) and (25), and the mean length of turns of the primary and secondary windings, *l*wp and *l*ws, are used to calculate the leakage inductance. The improved hybrid model is given as

$$L\_{\rm dW-H\rm I} = \frac{N^2 \mu\_0}{h\_{\rm eq}} \left[ \frac{l\_{\rm WP} d\_{\rm SP} m\_{\rm P}}{3} F\_{\rm WP} + \frac{l\_{\rm WS} d\_{\rm AG} m\_{\rm S}}{3} F\_{\rm WS} + S\_{\rm d} + \frac{l\_{\rm WP}}{m\_{\rm P}^2} \sum\_{q=1}^{q-m\_{\rm P}-1} q^2 d\_{\rm ip-q} + \frac{l\_{\rm new}}{m\_{\rm s}^2} \sum\_{r=1}^{r-m\_{\rm s}-1} r^2 d\_{\rm iso-r} \right] \tag{26}$$

where, *d*ip-*<sup>q</sup>* and *d*is-*<sup>r</sup>* refer to the width of the *q*th and *r*th inter-layer distance of the primary and secondary winding from the outside to the inside.

A winding model with a non-uniform arrangement of conductors is shown in Figure 9. The model parameters are: 18 turns of the primary winding, 11 turns of the secondary, the width of the main insulating layer is 8 mm, the thickness of the copper foil is 1 mm, the height is 220 mm, the height of the window is 245 mm. An 8 mm-wide air channel is located between the *q* layer and the *q* + 1 layer in the primary, and the other turns of the insulating layer have a width of 0.2 mm.

**Figure 9.** Winding model with non-uniform arrangement of conductor layers.

As shown in Figure 10, different models and FEM are used to calculate the leakage inductance of the example winding with different air channel locations. In the original and hybrid models, the average inter-layer width is used for the calculation. Both the FEM results and the calculation results of the proposed model increase as the position of the channel gradually approaches the main insulating layer, while the original model and the hybrid model cannot reflect the influence of the air channel on the leakage inductance. Based on the FEM results, the average error of the proposed model is 1.9%, and the maximum error is 2.4%, both smaller than the other two models; when *q* = 10–12, the errors of the three models are relatively close, all less than 5%.

**Figure 10.** Comparison of FEM and analytical calculation results. (**a**) Comparison of results, (**b**) Comparison of errors.

The additional leakage inductance is composed of multiple air-gap inductances. In [24], an empirical formula is provided for calculating the inductance of the lumped air-gap inductor. In [25], the 2D expressions of air-gap reluctance are derived based on the Schwarz– Christoffel transformation.

A simplified air gap is shown in Figure 11, and its reluctance is composed of the internal reluctance *R*in and the fringing reluctance *R*fr in parallel. The expressions are given as [26]

$$R\_{\rm in} = \frac{d\_{\rm g}}{\mu\_0 A\_{\rm ac}} \tag{27}$$

$$R\_{\rm fr} = \frac{\pi}{\mu\_0 C \ln\left(\frac{2h + d\_{\rm g}}{d\_{\rm g}}\right)}\tag{28}$$

$$R\_{\text{g}-i} = R\_{\text{fr}} |\, |R\_{\text{in}} \tag{29}$$

where *d*<sup>i</sup> is the height of the air gap, *A*ac is the cross-sectional area of the core, *μ*<sup>0</sup> is air permeability, and *C* is the perimeter of the cross-section.

**Figure 11.** Simplified model of the air gap.

Thus, the additional leakage inductance is given as

$$L\_{\text{ca}} = \frac{N^2}{\sum R\_{\text{g}} - i} \tag{30}$$

where *N* is the number of turns of the winding.

The FEM simulation is carried out to verify the calculation method for the additional leakage inductance, and the model and parameters are shown in Figure 12a. As shown in Figure 12b, as the air gap height increases, the inductance value is greatly reduced. Therefore, in the MFT design process, the height and number of the air gap are used to adjust the inductance value.

**Figure 12.** Calculation and simulation of the additional leakage inductance. (**a**) FEM model (**b**) comparison of the calculated and simulated values.

#### **4. Thermal Modeling and Analysis**

For the design of the air–water combined cooling MFT, a thermal model is necessary. The steady-state thermal model mainly includes the heat source, thermal resistance, and network topology. The heat source, namely the core and winding losses, has been discussed above. Therefore, this section mainly discusses the modeling of thermal resistance and network topology.

#### *4.1. Convective Thermal Resistance of the Air Channel*

The air channel is an important heat dissipation structure of transformer windings. The air channel can be simplified as a parallel-plate channel for thermal modeling. As shown in Figure 13, the length of the air channel is *L*, the cross-section is *W*×*S*, and *W* >> *S*. The hydraulic diameter is given as [27]

$$d\_h = \frac{4A}{P} = \frac{4WS}{2(W+S)} \approx 2S\tag{31}$$

where *A* is the cross-section area and *P* is the inner perimeter.

**Figure 13.** Parallel-plate air channel.

Since the convection model of the parallel-plate channel in literature is not complete, a combination of the laminar-flow model [28] for the parallel-plate channel and a turbulentflow model [28] for the non-circular tube is adopted in this paper. At the same time, linear interpolation is used in the transition period between laminar flow and turbulent flow. The Nusselt number calculation formula of the parallel-plate air channel is given as

$$N\_{\rm u} = \left\{ \begin{array}{l} \left[N\_1^3 + N\_2^3 + N\_3^3\right]^{1/3}, \left(R\_{\rm e} \le 2300\right) \\\left(\frac{\left(\frac{\mu}{8}\right)R\_{\rm e}P\_{\rm e}\left[1 + \left(\frac{d\_{\rm e}}{2}\right)^{2/3}\right] \left(\frac{T\_{\rm e}}{T\_{\rm w}}\right)^{0.45}}{1 + 12.7\sqrt{\frac{L\_{\rm e}}{8}\left(P\_{\rm r}^{2/3} - 1\right)}}, \left(R\_{\rm e} \ge 10^4\right) \\\ (1 - r)N\_{\rm u2300} + rN\_{\rm u10000}, \left(2300 < R\_{\rm e} < 10^4\right) \end{array} \right. \tag{32}$$

$$R\_{\mathbf{e}} = \frac{\rho v\_{\text{in}} d\_{\mathbf{h}}}{\eta} \tag{33}$$

$$N\_1 = 7.541\tag{34}$$

$$N\_2 = 1.841 \left( R\_\text{e} P\_\text{r} d\_\text{h} / L \right)^{1/3} \tag{35}$$

$$N\_3 = \left(\frac{2}{1 + 22P\_\mathrm{r}}\right)^{1/6} \left(R\_\mathrm{e} P\_\mathrm{r} d\_\mathrm{h} / L\right)^{1/2} \tag{36}$$

$$f\_{\mathbf{r}} = \left(1.8 \log\_{10} R\_{\mathbf{e}} - 1.5\right)^{-2} \tag{37}$$

$$P\_{\mathbf{f}} = \frac{c\_{\mathbf{p}}\eta}{\lambda\_{\mathbf{f}}} \tag{38}$$

$$T\_{\rm f} = \frac{T\_{\rm i} + T\_{\rm o}}{2} \tag{39}$$

$$r = \frac{R\_{\text{e}} - 2300}{10^4 - 2300} \tag{40}$$

where *R*<sup>e</sup> is the Reynolds number; *N*1, *N*2, and *N*<sup>3</sup> are the laminar flow model coefficients; *P*<sup>r</sup> is the Prandtl number; *T*<sup>f</sup> is the average fluid temperature; *T*<sup>i</sup> is the inlet temperature, which is equal to the ambient environment temperature *T*am; *T*o is the outlet temperature; *T*<sup>w</sup> is the wall temperature; *f*<sup>r</sup> is the friction coefficient; *r* is the interpolation coefficient in the transition period; *ρ* is the fluid density; *η* is dynamic viscosity; *c*<sup>p</sup> is fluid specific heat capacity; and *λ*<sup>f</sup> is fluid thermal conductivity.

After the Nusselt number is obtained, the surface heat transfer coefficient is given as

$$\hbar = \frac{\lambda\_{\text{f}} N\_{\text{u}}}{d\_{\text{h}}} \tag{41}$$

The FEM simulation is carried out to verify the accuracy of the above calculation methods. Figure 14 shows the temperature distribution inside the air channel under different wind velocities. The simulated and calculated results are compared in Figure 15. As shown in Figure 15a, the average heat transfer coefficient increases with the higher wind velocity. The reason is that the higher wind velocity reduces the thickness of the boundary layer. As shown in Figure 15b, the shorter channel can achieve higher average heat transfer coefficients. The reason is that the thickness of the boundary layer increases with the channel length. The error between the calculated and simulated result is less than 10%.

**Figure 14.** Temperature distribution of the air channel under different wind velocities.

**Figure 15.** Comparison of the simulated and calculated heat transfer coefficients. (**a**) Different wind velocities; (**b**) different channel lengths.

After obtaining the heat transfer coefficient, the convective thermal resistance can be obtained according to the Newton cooling formula of the internal convection.

$$Q = hA \left( T\_{\rm w} - \frac{T\_{\rm i} + T\_{\rm o}}{2} \right) = \frac{\left( T\_{\rm w} - T\_{\rm f} \right)}{R\_{\rm invency}'} \tag{42}$$

In (42), the temperature difference is between the surface temperature *T*w and the fluid temperature *T*f. As the fluid temperature *T*<sup>f</sup> is related to the internal convection process, the thermal resistance *R*'inconv defined by (42) cannot be directly connected to the external environment nodes in the thermal network. Therefore, transform (42) as

$$Q = \frac{hA\left(T\_\mathrm{W} - \frac{T\_\mathrm{i} + T\_\mathrm{o}}{2}\right)}{\left(T\_\mathrm{W} - T\_\mathrm{i}\right)}(T\_\mathrm{W} - T\_\mathrm{i}) = \frac{\left(T\_\mathrm{W} - T\_\mathrm{i}\right)}{R\_\mathrm{incconv}}\tag{43}$$

The thermal resistance *R*inconv defined by (43) can be directly connected to the external environment nodes for the temperature difference is between surface temperature *T*w and inlet temperature *T*i. The proposed thermal resistance definition for the internal convection is more convenient for thermal network modeling.

#### *4.2. Conduction and Radiation Thermal Resistance*

Conduction: for a uniform cuboid with length *l*, thermal conductivity *λ*, and crosssectional area *A*, the thermal resistance is given as [27]

$$R\_{\text{cond}} = \frac{l}{\lambda A} \tag{44}$$

Radiation: for a surface with the area *A*sur, the surface temperature *T*sur, and the ambient temperature *T*am, the thermal radiation resistance is given as [27]

$$R\_{\rm rad} = \frac{T\_{\rm sur} - T\_{\rm am}}{\varepsilon\_{\rm i} \sigma A\_{\rm sur} (T\_{\rm sur}^4 - T\_{\rm am}^4)} \tag{45}$$

where *σ* is Boltzmann constant and *ε*<sup>i</sup> is emissivity factor.

#### *4.3. Thermal Network Topology*

The structure of the air–water cooled CW ELII MFT is shown in Figure 16. As shown in Figure 16b, to reduce the eddy current loss caused by the fringing flux, the tooth-slot structure is adopted in water-cooling plates adjacent to the additional core.

**Figure 16.** The structure of CW ELII MFT. (**a**) Cores; (**b**) water-cooling system; (**c**) windings; (**d**) packaging and fans.

As shown in Figure 17, the thermal network model has ten nodes. *N*1~*N*<sup>3</sup> denote the upper yoke, the limb, and the lower yoke of the main core, respectively. *N*<sup>4</sup> represents the secondary winding. *N*<sup>5</sup> represents the primary winding. *N*6~*N*<sup>8</sup> represent the upper yoke, the limb, and the lower yoke of the additional core. *N*<sup>9</sup> represents the external air. *N*<sup>10</sup> represents the water-cooling plate. The additional nodes *n*1~*n*<sup>8</sup> represent the convection and radiation surfaces of nodes *N*1~*N*<sup>8</sup> and are used to provide the wall temperature for convection and radiant thermal resistance calculations. *Q*1–*Q*<sup>8</sup> represent the losses of each node, respectively; *T*am and *T*co represent the temperature of the environment and the water-cooling plate. In the above model, the internal forced convection process of the water-cooling plate is ignored. The implication of *Rxy*-Cond/Conv/Rad is conduction/convection/radiation thermal resistance from node *x* to node *y*. In a steady-state, the governing equation is given as [12]

$$0 = AT + BU\tag{46}$$

where *A* describes the thermal resistance network of internal nodes; *B* describes the thermal resistance network between the internal nodes and the external nodes (air and water-cooling plates); *T* is the node temperature vector; *U* is the external excitation vector.

**Figure 17.** Thermal network model of the CW ELII MFT.

#### **5. Optimal Design, Analysis, and Example**

The MFT design is a non-linear, non-convex, and non-continuous problem. The optimization design and the comparative evaluation of the inductor-integrated MFT is presented in this section.

#### *5.1. Design Method*

The critical factors of optimal design include: design inputs, constraints, and optimization objectives.

#### 5.1.1. Design Inputs

As shown in Figure 18, there are many variables in the optimal design of MFT, and it is unrealistic to optimize all variables. Therefore, seven free variables are selected in this paper, including the flux density of main core *B*m, the flux density of additional core *B*ma, turns of the primary winding *N*p, the height of foils *h*c, the thickness of primary foils *d*cp, the thickness of secondary foils *d*cs, and the width of core tape *D*cs. The free variables constitute the input vector of the optimization design as

$$\mathbf{x} = \begin{bmatrix} B\_{\text{m}\prime} \, B\_{\text{ma}\prime} \, \text{N}\_{\text{p}\prime} \, \text{h}\_{\text{c}\prime} \, d\_{\text{cp}\prime} \, d\_{\text{cs}\prime} \, \text{D}\_{\text{cs}} \end{bmatrix} \tag{47}$$

**Figure 18.** Sketch of the CW ELII MFT structure.

In addition, inputs design includes many fixed parameters about operating conditions and material properties, such as the DC voltage *U*dc1 and *U*dc2, the turns ratio of the transformer *k*PS, frequency *f*, and electrical and thermal properties of materials.
