Analytical Model for HF-Losses Caused by 2D Magnetic Fields in Litz Wire

The 2D magnetic field effect arises when the winding height is lower than the window height due to isolation requirements and mechanical constraints. In this situation, most analytical models for calculating HF-losses in Litz wire fail to predict the winding losses accurately since they are normally based on a 1D magnetic field assumption. Although numerical or semi-numerical methods can accurately calculate the winding losses, they are too time-consuming to be integrated into the converter optimization routines. This article provides a new loss model, which is validated to have a much lower error (<10%) than 1D field loss models (up to −45%) by numerical calculations and experiments. The model is presented for two winding transformers with layered windings and non-gapped cores with high permeability in the frequency range, where the strand diameter to skin depth ratio is equal to 1.2. Moreover, the model can also be extended for interleaved windings by separating the windings to multiple blocks. Furthermore, the computation time is less than 300µs using the algorithm provided in this article. This new model achieves a good balance between accuracy and computational time and has great potential to improve the transformer as well as the converter design.


I. INTRODUCTION
Magnetic components such as medium frequency transformers play a key role in many power electronic systems. Because of the high-frequency (HF) effects, the winding losses increase rapidly with frequency and limit the efficiency and the power density of the magnetic components as well as the power electronic systems. Litz wire (LW) consisting of numerous twisted strands is often used to minimize the winding losses. With perfect twisting in theory, each strand in the LW conducts the same current and the bundle level proximity effect is canceled [1]. Although perfect twisting is assumed, calculating the HF-losses in LW windings by numerical methods, such as FEM, results in high computational effort. By contrast, analytical models are typically much faster and are preferred in optimization routines where the losses must be evaluated several thousand times as explained in [2], [3].
In [4] four different such analytical models presented in [5], [6], [7], [8], [9], [10], have been evaluated and compared with respect to accuracy and computational time. In all of these models, a 1D external magnetic field is assumed, i.e. the fields are considered to be in parallel to the winding layers and have a constant amplitude along each field line. However, this assumption/simplification is not fulfilled, if the winding height is not the same as the (core) window height. Due to isolation requirements and mechanical constraints, the winding height must typically be lower than the window height, e.g. in medium voltage (MV) where the nominal voltage is usually between 1 kV and 35 kV, the required isolation distance between the winding and the core is larger than it is in the low voltage, where the nominal voltage is below 1 kV. In the case of a large distance between the core yoke and the winding, the magnetic field disperses out of the winding in different directions and has a more 2D distribution. In some publications (e.g. [11] and [12]), this 2D field effect is investigated as 'edge effect'. The results presented in [4] demonstrate that considerable error occurs, if the mentioned four models are used to calculate the winding losses in presence of such a 2D field. This error increases rapidly with the contribution of losses caused by the proximity effect and the distance between the winding and the yoke, which can be measured in the average window to winding height ratio (h w /h avg ). Compared to measurement results, the error of 1D field loss model [5] could be up to −45% for multi-layered windings with h w /h avg > 1. 5.
To calculate the losses caused by such a 2D magnetic field more accurately, FEM is used in [11], [13], [14], [15]. In [16], [17], [18], [19], the losses are calculated by a semi-numerical method, which calculates the magnetic field by FEM and uses this magnetic field in the analytical model to estimate the losses. However, the computational time of both methods is too long for converter optimization routines, since the loss calculation needs to be iterated numerous times. To reduce the computational time of FEM simulation, in [20] the PEEC (Partial element equivalent circuit) method is used to calculate the HF-losses based on the magnetic field calculated by 2D FEM. Moreover, in [21] and [22], the 2D LW geometry with numerous realistic strands is simplified to a homogeneous winding region with an equivalent complex permeability, of which the number of mesh elements can be substantially reduced. However, both methods for reducing the calculation time are still too computationally expensive to be integrated into optimization routines.
Compared to numerical methods, analytical models are faster and preferred in converter optimization routines. In [23], the 2D magnetic field is calculated by using the method given in [24]. However, only planar windings located in the air and above a magnetic substrate are considered. In [25], an analytical loss model considering the effects of twisting is presented, but no model for calculating the 2D magnetic field is given. In [26] and [27], the winding losses are derived based on iteratively calculated magnetic field and current density for round conductors, as the external magnetic field and the induced eddy current interact with each other. This model is accurate and suitable for any winding arrangements and core types with/without air gap(s). However, in LW the contribution of each strand to the magnetic field must be considered, and the losses converge after at least 3 or 4 iterations. Depending on the number of strands and total number of turns in litz wire windings, the total number of round conductors considered in the loss calculation can easily exceed several thousands. Hence, this model needs serval minutes to calculate the winding losses for litz wire windings. Moreover, if the windings are enclosed in the core window, the windings must be mirrored to each core boundary multiple times, as the contribution of the strand images to the magnetic field and losses must also be considered. In this case, the total number of conductors as well as the computational effort increase multiple times. To calculate the losses of litz wire windings enclosed in a core window could take several hours, which is unacceptable for converter optimization routines.
In [28], the "direct integration method" for calculating the losses caused by the proximity effect in a round conductor is proposed. In this method, the flux induced by the eddy current is integrated over the conductor area to calculate the losses in that conductor, based on the superposition of magnetic fields caused by all other conductors. To avoid this integral of flux and reduce the calculation time, the "three orthogonal fields method" is presented. Using this simplified method, the losses caused by the proximity effect can be calculated for 1000 round conductors within around 75 seconds. However, as mentioned before, litz wire windings in a transformer usually contain several thousand of strands. The calculation time is up to several minutes or several hours for litz wires enclosed in a core window. Therefore, this method is also very computationally expensive for converter optimization routines. Besides the loss model for LW, the 2D magnetic field model is also an essential part of calculating the loss of LW. In [29], [30], [31], [32], the 2D magnetic field is calculated for rectangular winding layers. As the winding is enclosed in the core, the results contain an infinite Fourier series. Although the fourier series converge relatively fast, the spatial average of the squared field in the winding region is normally required for the loss calculation, which is difficult or computationally expensive to calculate based on this Fourier series. In [5], [33], [34] the method of mirror images is applied for calculating the 2D magnetic field for solid wire and LW. To calculate the magnetic field accurately, the complete windings must be mirrored on the 4 window edges for a closed core window [3]. Thus, the method of mirror images is also relatively time consuming.
In [35], a closed-form model to estimate the HF losses caused by 2D magnetic fields is presented. This model considers only transformer windings with the same height and no isolation layer between the primary and the secondary winding. However, in real transformer the winding height of the primary and the secondary or even the winding height of each layer in the same winding could be different, as shown in [3], [34], [36], [37]. Moreover, for some assumptions made in [35], a detailed physical explanation is missing, and the model is not validated by the experiment. To overcome these limitations, the model presented in [35] is adopted in this article, and more detailed physical explanation is added. Note that this article focuses on two winding transformers with layered windings and ungapped cores with high permeability. The model can also be extended for transformers with interleaved windings. A simple example for such an extension is given in the conclusion section. The frequency range of interest is where the HF-losses in the transformer with litz wire are lower than the one with an equivalent round wire, which has the same DC resistance as the litz wire. According to [4], litz wires result only in lower HF-losses than equivalent round wires in the frequency range, where the strand diameter is smaller than the skin depth, i.e. d s /δ < 1. As the litz wire should be avoided in transformer design beyond that frequency range, this article focuses only on the frequency range (d s /δ < 1.2), where litz wires are beneficial. Furthermore, this article also presents an algorithm that is 30 times faster than the one presented in [35]. An experimental validation for the proposed model based on multiple transformer designs is also conducted in this article.
This article is arranged as follows. Section II presents the 2D magnetic field model. In Section III, the modified loss model is derived. In order to implement the new model computational efficiently, an algorithm is derived in Section IV. In Section V, the new model is numerically and experimentally evaluated and compared to the 1D-field loss model presented in [5] with respect to accuracy and computational time. The equations for calculating the magnetic field are summarized in Appendix A.

II. 2D MAGNETIC FIELD CALCULATION
In this section, the model for calculating the 2D magnetic field is derived. As the magnetic field is 2D, it is difficult and not necessary to distinguish the internal field from the external field for the loss calculation. In this article, the losses caused by the proximity effect are calculated by the total magnetic field through the winding, thus, only two types of HF-losses (losses caused by the skin and the proximity effect) are considered. The losses caused by the external proximity effect P p increase much faster with frequency than the losses caused by the skin effect (P Skin ), and dominate the total losses in the high frequency range. Furthermore, the accuracy for calculating the losses caused by the skin effect depends on whether the current is equally shared between strands, which is controlled primarily by the twisting and wire termination, is not affected by the winding arrangement. Therefore, if perfect twisting is assumed, the accuracy of the loss model depends principally on whether P p is precisely calculated.
In the presented model it is assumed, that an external magnetic field exists only outside the core, as an infinite permeability of the core (μ r,core → ∞) is assumed and the amplitude of the magnetic field in the core becomes zero (H core ≈ 0). If the magnetic field distribution in the winding region (H WR ) is known, the 2D losses per unit length caused by the proximity effect for a single strand in a certain winding layer can be given by: Equation (1) implies that H WR,RMS is the spatial RMS value of the peak value over the sinusoidal variation rather than the time RMS value that is commonly used in electrical engineering. Furthermore, the area S in (1) refers to the rectangular area outlining a winding layer, since the losses are calculated layer by layer in this article. The winding region (WR) is composed by these areas as the region highlighted in dark blue in Fig. 1(b). The factor G( f ) in (1) is assumed to be related to the material, the conductor shape and the frequency and independent of the 2D magnetic field. Hence, only H WR and H WR,RMS need to be derived in order to include the 2D H-field effect in the loss calculation.

A. STRAIGHT FIELD LINE APPROXIMATION
The magnetic field in the window can be derived based on Ampere's law. In order to determine the integral path of the magnetic field, a straight field line approximation is proposed. In the following, a transformer with a one-layer primary and a two-layer secondary winding (see Fig. 1) is considered as an example for deriving the basic idea of the straight field line approximation. Since the core window is symmetric to the midline, only the upper half of the core window is shown. In Fig. 1(a) and (b), the magnetic field distribution calculated by FEM is shown. In many applications, coil formers/bobbins are used to simplify the manufacturing and to achieve a proper winding arrangement, as well as to provide the necessary electrical isolation. As the bobbins or isolation materials generate neither eddy-current losses nor influences the 2D H-field distribution. The bobbins or isolation materials can be treated as an air gap between the primary and the secondary winding. This is indicated with a green bar shown in Fig. 1(a). As the magnetic field does not change significantly within the gap between the primary and the secondary winding (cf. Fig. 1(a)), the gap does in a first approximation not affect the magnetic field distribution inside the windings. Therefore, the gap between the primary and the secondary winding could be eliminated for the winding loss calculation by moving the secondary winding right next to the primary winding as shown in Fig. 1(b), where the straight field line approximation is applied.
Since the winding height (h p or h s ) is smaller than the window height (h w ) as shown in Fig. 1(b), a share of magnetic field lines does not reach up to the yoke but enters directly into the center or the outer core legs. Therefore, the magnetic field is more 2D between the yokes and the winding region (WR) as shown in Fig. 1(b), i.e. it consists of both an x-and a y-component. Within the WR the magnetic field is approximately 1D, i.e. it has only a y-component in the considered case. Thus, the upper half core window can be divided into one 2D field region (2DFR) with a 2D (H-field) and one WR with a 1D H-field as shown in Fig. 1(b) by the yellow and the purple areas. In order to calculate the length of each field line, the curved field lines in the 2DFR as well as the WR in Fig. 1(b) are approximated by the straight field lines as given in Fig. 1(c). In the 2DFR, the norm of the magnetic field vector H norm = | H x | 2 + | H y | 2 is calculated, since in Ampere's law only the magnetic field amplitude along the field line is considered.
As the winding is symmetric to the midline and μ r does not change in the window, the flux density B distributed on the midline has only a y-component. In addition, LW is usually applied in a frequency range where the skin depth is bigger than the strand diameter. Otherwise, LW produces more losses than its equivalent round wire. In this frequency range, the winding can be replaced by a homogenously distributed current density region, which is abbreviated by HCR ( Fig. 1(c)). With the HCR, the flux density changes linearly along the midline as shown in Fig. 1(c).
Moreover, the window can be divided into 4 regions (R1-R4), based on the core parts that the H-fields enter as shown in Fig. 1(c). In R1 and R4, the H-fields enter the center or the outer core leg, whereas in R2 and R3 the H-fields enter the core yoke. The boundaries of these 4 regions can be determined by 3 points (x B1 , x FBL and x B2 shown in Fig. 1(c)). Since the coordinates of the 3 points are different in the different coordinate systems, to keep the derivation more compact, the coordinate system must fulfill the following 3 rules.
1) The x axis of the coordinate system must be in parallel to the midline. 2) The WR must be located on the positive x half-axis 3) The origin of the coordinate system must be located on the outer edge of the winding of interest. Based on the 3 rules, a coordinate system is chosen as shown in Fig. 1(c) for solving the H-field of the secondary winding. In the case of solving the H-field of the primary winding, the coordinate system must be reversed and shifted to the outer edge of the primary winding as shown in Fig. 1(c).
To calculate the H-field in the different regions, the position x FBL of the field boundary line (FBL), which separates the magnetic fields that bent to left and right, needs to be determined in the first step. To simplify the calculation, the primary and the secondary winding are considered as 2 winding blocks with the same height as shown in Fig. 2(a). As the core does not affect the location of FBL, the core can be neglected in the calculation of x FBL . Moreover, the MMF caused by primary and the secondary winding has the same amplitude but in the opposite direction. As shown in Fig. 2(a), a share of H-fields is bent to the left and closed in the left hand side of the FBL, the rest magnetic fields are bent to the right and closed in the right hand side of the FBL. Therefore, in the transition between the left bent and the right bent field lines, a straight field line must exist, i.e. the x-component H x of the magnetic field on the FBL is zero. As the same winding height of the primary and the secondary winding is assumed, scaling both windings with the same factor in height does not affect the location of FBL, consequently the winding block can be further simplified to 2 line current sources as shown in Fig. 2(b). The 2 line current sources with length −x a and x b are located both on the x-axis. The x-component of the H-field H sx at a random point (x, y) caused by a random point current in the secondary line current source can be given by (2), where J s is the current density of the secondary line current source.
The x-component of the H-field caused by the secondary line current source (H sx ) is equal to the integral of the current density J s over the interval [0, x b ] and is given by: Similar to the secondary line current source, the x-component of the H-field (H px ) caused by the primary line current source is given by: The variable θ 1 and θ 2 in (3) and (4) are given by: As the current (I pb and I sb ) in the primary and the secondary line current source flows in the opposite direction, the primary current density J pb is given by: As the sum of H px and H sx on the FBL is zero everywhere and is independent of variable y, where x FBL is the coordinate of FBL in the coordinate system in Fig. 2(b), the value of y is set to 1 for simplicity. Substituting (3), (4) and (5) into (6), (7) can be derived.
Unfortunately, there is no general and unique solution for x FBL in (7). Nevertheless, x FBL can be solved numerically for given winding widths (x a , x b ). However, the computational time (≈20 ms) of this numerical solution is relatively high.
To reduce the computational time, the numerical solution can be replaced by a fitted curve function. As shown in (7) the value of x FBL depends both on x a and x b . To avoid 3D curve fitting, the number of variables need to be reduced to one. As in the a small interval close to origin (|x| 1), the function arctan(x) can be considered as a linear function, i.e arctan(ax) = a arctan(x), x a can be extracted from the arctanfunction and eliminated. The number of variables in (7) can be reduced to 2 and expressed as a function f (x * FBL , x r ) as given by (8) The solution of x FBL in (7) can be calculated by solving equa- FBL , x r ) and the surface f = 0, as shown in Fig. 3(b). On the cut line for each x r , there is a value of x * FBL which satisfies the equation f (x * FBL , x r ) = 0. Based on these numerical results, a 2D curve fitting, which is much faster than the 3D curve fitting, can be used to predict the value of x * FBL as shown in Fig. 3(a). The equation of the fitted curve is given by: Finally, the value of x FBL in the coordinate system shown in Fig. 1(c) can be calculated as: To evaluate the derivation of x FBL , the x-component of the magnetic field H x is calculated for an exemplary winding as  shown in Fig. 4(a). The current in the primary and the secondary winding is 10 A and −10 A. The variable x r , x * FBL,CF , and x FBL for this winding setup is equal to -1.5, 0.1341 and 13.66 mm. The FBL calculated by FEM is given by the blue line in Fig. 4(a), which is not entirely straight as assumed. Nevertheless, the average value of the x coordinates of all the points on the FBL calculated by FEM is equal to 13.65 mm, which is very close to the value 13.66 mm calculated by the proposed model. In Fig. 4(b), H x along the FBL calculated by the proposed model is shown. As expected, H x on the FBL is almost zero, except the line sections, which are near the top and the bottom edge of the winding.
In order to determine the coordinates of points x B1 and x B2 as well as the field amplitude H WR , three assumptions need to be made (A1 and A3 are inherited from [35]): r A1:A linearly increasing flux (φ approx ) along the window edges (points P 1 -P 2 -P 3 -P 4 in Fig. 1(b) and P 1 -P 2 -P 3 -P 4 in Fig. 1(c) is assumed. A comparison to FEM results calculated for the considered case shown in Fig. 1(b) is given in Fig. 5(a). r A2: In the WR, a homogenously distributed current density region HCR as shown in Fig. 1(c) is assumed. Furthermore, the isolation thickness between two adjacent layers is neglected so that the magnetic field increases continuously from layer to layer. Thus, the current per unit width of the WR for the primary and the secondary winding is given by (11).
There, S p and S s are the rectangle areas of the primary and the secondary winding.
r A3: A piece-wise linear amplitude H PL of the H-field along each approximated straight field line is assumed as shown in Fig. 5(b). In the 2DFR, the amplitude of the H-field H 2DFR is considered to increase/decrease with a constant slope k p or k s along the length of the field line l for the primary or the secondary winding. Therefore, H 2DFR of the primary and the secondary winding are given by: H 2DFR,p = k p l + b p and H 2DFR,s = k s l + b s , where b p = b s = 0, since l starts at the window boundary and H 2DFR,p (l = 0) = H 2DFR,s (l = 0) = 0. In the WR, a constant amplitude of the H-field H WR along the field line is assumed, since the H-field in the WR is 1D as shown in Fig. 1(c). Note that in this assumption the MMF (magnetomotive force) which is the area surrounded by H PL and H FEM as shown in Fig. 5, must be equal.

B. MAGNETIC FIELD CALCULATION FOR SECONDARY WINDING
In this section, the RMS value of the magnetic field in the WR H WR,RMS is derived for the secondary winding. Since H WR has only a y component and is constant along the field line, the RMS value of H WR for the nth layer H WR,RMS is defined by (12). H WR,RMS can be derived in 5 steps as shown in Fig. 6, where the 1st step, i.e. the definition of the coordinate system has been discussed in the last section and the related parameters are shown. r Step 2: Since each winding layer may consist of a different number of turns and has therefore a different winding height, H WR,RMS needs to be determined for each layer separately. H WR,RMS varies from region to region, hence the region(s) where each winding layer is located in must be determined first. Therefore, the positions of x B1 and x B2 are calculated in this step with derived x FBL from the last section.
In region R1, the flux per unit length φ B1 through the line segment L 0,x B1 on the midline is equal to the flux per unit length φ P 3 −P 4 through the line segment L P 3 ,P 4 as shown in Fig. 1(c). Based on assumption A1 and the linearly increased flux density from x = 0 to x = x FBL shown at the bottom of Fig. 1(c), x B1 can be derived as given in (13). Note that for windings where the winding height of each layer is different, the average winding height h ap and h as is used for calculating d yp and d ys of the primary and the secondary winding (e.g. Analog to deriving the position x B1 , the flux per unit length φ B32 through line segment L x B3 ,x B2 is equal to the flux per unit length through line segment L P 1 ,P 2 as shown in Fig. 1(c). Based on assumption A1, the length d B23 of line segment L x B3 ,x B2 is given by (14). The position of x B2 is equal to (14) r Step 3: In this step, the lengths of the straight field lines as shown in Fig. 1(c) are calculated for each region.

1) REGION 1
In Fig. 7(a), a random magnetic field line (dashed line), which consists of two parts, the 1D part in the WR and the 2D part in the 2DFR, is shown. The 1D part in the WR is as long as half of the winding height, and the length of the 2D part in the 2DFR is equal to , where the unknown value of d bP 4 is given by (15).
Because the winding is symmetric to the midline, the total length of the field line in R1 is equal to l R1 = h s + 2d ab .

2) REGION 2
Similar to R1, the length of the 2D part of the field line is equal 2 where the unknown value of d bP3 must be derived (see Fig. 7(b). Based on assumption A1, the flux density φ x−x B1 through line segment L x,x B1 is equal to the flux per unit length through the line segment L b,P3 , so that d bP3 is given by (16).

3) REGION 3
As shown in Fig. 7(c), only a very small part of the secondary winding is located in R3 in the considered example. Similar to R2, the length of the 2D part of the field line can be calcu- Therefore, the only unknown variable is d bF which can be derived with (17), which is again based on assumption A1.
r Step 4: In this step, the amplitude of the magnetic field in the winding region H WR is calculated. Based on Ampere's law, the integral of the H-field along a closed loop, which includes distance d ab in both 2DFRs, and the winding height h s in the WR, is equal to the current through that loop as given by (18). As can be seen in Fig. 5(b), the magnetic field in the 2DFR H 2DFR is assumed to increase linearly with field line length l with a slope k p or k s , which can be derived by applying Ampere's law to a specific field line. Since the slope varies with field lines, in this article the field line which separates the flux in a region into 2 equal parts is considered as the mid field line, which is chosen to calculate the slope (k p or k s ) for each region. k s can be calculated for each region based on (18) and an example for region R2 is given by (19).
There, (18) and (19), a general form for the H-field H WR in the WR can be derived for the secondary winding as given in (20). By substituting the corresponding values of k s and d ab calculated in step 2, H WR (x) for different regions can be derived. The resulting equations are given in Appendix A. 20) r For a layer which is located in 2 regions (e.g. the first layer of secondary winding in Fig. 1(b) is located in both region 1 and 2), H 2 WR,RMS is given by (22).

C. MAGNETIC FIELD CALCULATION FOR PRIMARY WINDING
The magnetic field for the primary winding in the example transformer case is calculated in this section. The derivation is similar to the one in the secondary winding. As step 5 for both windings is identical, in this section only step 1 to step 4 are discussed.
r Step 1: In order to reuse some calculation results from Section II-B and simplify the derivation, the coordinate system built for calculating the H-field of the secondary winding needs to be reversed horizontally as shown in Fig. 8. The coordinate system must fulfill the 3 rules given in Section II-A, but the origin of the coordinate system aligns now with the outer edge of the primary winding.
r Step 2: In the reversed coordinate system, the boundary positions x FBLR and x B2R as shown in Fig. 8 are equal to a t − x FBL and d B23 . r Step 3: As shown in Fig. 2(c) the primary winding is located in regions R3 and R4, and the length of straight field lines is calculated for these two regions in this step.

1) REGION 4
The derivation of the field line length is the same as it is in region R1 for the secondary winding. As can be seen in Fig. 8(a), the length of the field line in the 2DFR is equal to d ab = (x + d xi ) 2 + d 2 bP 1 , where the unknown value of d bP 1 is given by: 2) REGION 3 Similar to the derivation for secondary winding in region R3, the length of the 2D part of the field line is equal to Fig. 8(b). The unknown value of d bP 2 is given by: r Step 4: Also, the mid field line in each region is chosen for calculating the gradient k p as it is done for the secondary winding. The gradient k p3 can be calculated for example for region R3 by: There, The general equation for calculating the magnetic field in the WR for primary winding can be therefore given by: The method introduced here is based on considerably strong assumptions (A1 and A3) and cannot be used to accurately estimate the magnetic field everywhere in the window. However, an accurate estimation of the magnetic field everywhere is not necessary for the loss calculation, as the accuracy of the field calculation is not equivalent to the accuracy of loss calculation. Moreover, the equations of H WR for different regions are polynomials that can be analytically integrated and be relatively faster computed than other equations obtained by solving Maxwell's equations directly.

III. LOSS CALCULATION
The following calculation of winding losses is derived by modifying the two models presented in [5] and [9], which are based on Bessel and hyperbolic functions. By modifying the loss contribution P p , the accuracy of both models can be significantly improved. The modification of the model based on Bessel functions can be simply performed by substituting the external magnetic field in P p by H WR,RMS calculated for each layer from the last section. However, the model based on hyperbolic functions needs to be adapted in 2 steps. Since the calculated 2D magnetic field H WR,RMS is only applied to the losses caused by the external proximity effect, the loss equations for the skin and the proximity effect have to be separated as given in [4]. The loss equation of P p for each strand can be derived and given by: where h s = √ π 2 d s is the thickness of the transformed foil winding. Accordingly, the separated equation of losses caused by the skin effect for a single turn is given by: Note that the derived loss (27) and (30) are valid for both primary and secondary winding. Compared to the 1D field loss model in [5], the main difference is the field calculation, as the 2D field effect is considered in the proposed model. Therefore, the proposed model needs 2 more input parameters, including the distance between winding and core yoke (d y ) and the distance between winding and center core leg (d xi ). All other parameters in the proposed model can be derived based on the input parameters. Moreover, some algorithm, that detects in which region(s) that each winding layer is located, needs to be implemented in order to select the correct function(s) for calculating the magnetic field. Therefore, the losses P p must be calculated separately for each layer. Consequently, the calculation effort increases with the number of layers. To reduce calculation time due to the increasing number of layers, an algorithm which assembles the calculation of losses for different layers into a single matrix system, is introduced in the next section.

IV. ALGORITHM
To calculate the RMS value of H WR for each layer, the integral of H 2 WR needs to be calculated. In the previous work [35], the integral is calculated by using the Matlab function 'integral'. Furthermore, the integral needs to be performed at least n times for an n-layer winding. If the layers are located in more than one region, the integral must be conducted up to n + 3 times. The average calculation time for the 13 different transformer cases presented in [35] is around 8.3 ms, which can be significantly improved by using the algorithm presented in this section.
The basic idea of the algorithm is to use a single matrix system to represent all the integrals for a winding. Instead of integrating H 2 WR for each layer separately, with the matrix system the integrals of H 2 WR of all layers can be evaluated in a single step. Furthermore, this algorithm accelerates the computation further by getting rid of using the Matlab function 'integral'. The algorithm is explained for the secondary winding as an example in the following.
For the secondary winding, which is located in 3 regions (R1-R3), all the integrals of H 2 WR can be arranged in a single matrix system as given in (33). There, N IR1 -N IR3 are the where A is the coefficient matrix and X is the variable matrix. Matrix A is different for different windings and regions, e.g. the coefficient matrix for the secondary winding in region 1 is Each definite integral in (33)  where H s,Int is a 1 × 2N Is matrix (N Is = (N IR1s + N IR2s + N IR3s )) containing the integral value of H 2 WR at each node of integral intervals. For a layer which is located in a single region, such as the first layer in region R1, the integral The complete loss calculation procedure using this algorithm is shown in Fig. 9. At the beginning, the indefinite integrals of the functions H WR (x) in each region need to be calculated. The coefficient and variables are extracted from the indefinite integrals of the functions H WR (x) in each region as given in (32). Furthermore, to assemble the matrix system given in (35), the number of integrals in each region (N I ) and the operation matrix need to be further calculated based on input parameters. Finally, H 2 WR,RMS and the total winding losses can be calculated based on the matrix system. This algorithm uses a single matrix system to calculate H WR,RMS for all winding layers, which significantly reduces the calculation time of the proposed model. The detailed evaluation of this algorithm as well as the proposed model is given in the next section.

V. EVALUATION AND COMPARISON
In this section, the proposed model is evaluated by the numerical method and the experiments. 2D FEM is chosen for the numerical validation, as the proposed model is 2D, the 3D effects such as the bundle level HF-effect of litz wire and the influence of the winding part located outside the core window are not considered. In the experiment, 5 transformers based on E/ETD and pot cores and 3 different litz wires are designed and implemented. The experimental validation is done by comparing the winding resistance calculated by the proposed model to the winding resistance measured by impedance analyzer.

A. NUMERICAL VALIDATION
The field model proposed in Section II is validated first with FEM simulations for transformer designs 5 and 10 shown in Fig. 14. As in the proposed model and other loss models, the losses of each strand in the same layer are assumed to be identical and the losses are calculated layer by layer, the magnetic field is also compared layer-wise. Moreover, the proposed model uses the spatial RMS value of the magnetic field for each layer (see (1)) in the loss calculation, thus this spatial RMS value of the magnetic field is also calculated by FEM for comparison. In Fig. 10(a1) and (b1), the magnetic field calculated by the proposed model H WR,RMS and the 1D model H 1D , and the average value of the 1D field H avg are compared to the value H FEM calculated by FEM, where H 1D and H avg for the nth layer in the primary or the secondary winding are given by: In (38) and (39), N l , I and h l are the number of turns pro layer, the current in one turn, and the winding height of the considered winding layer. In (40), S n is the rectangle area including the conductor area outlining the nth winding layer, which is indicated by different colors for transformer designs 5 and 10 in Fig. 14.
As can be seen in Fig. 10(a1) and (b1), the proposed field model always underestimates the magnetic fields through the winding layers close to the center and the outer core legs. This can be seen for example for P1 and S2 in transformer design 5, and P1, S2 and S3 in transformer design 10 shown in Fig. 14. The winding layers close to the center and the outer core legs are usually located in regions R1 and R4, so that the field lines through these winding layers in the simulation are not 1D and not parallel to the winding layer as that is assumed in the proposed model. The field lines are rather 2D and even go perpendicularly through the winding layer and enter the core leg. Therefore, the resulting field line length in the simulation is smaller than it is assumed in the proposed model. This results in a higher H-field amplitude in the simulation than it in the proposed model. Moreover, the field amplitude of the winding layers located in regions R2 and R3 is mainly affected by the relative position of the primary and the secondary winding. If the height difference between the primary and the secondary winding is small, so that the field line through the winding layers located in regions R2 and R3 are relatively straight, the estimated field amplitude by the proposed model matches well to the one in the simulation, as shown in Fig. 10(a1). However, if the height difference between the primary and the secondary winding is relatively large, so that the field lines are bent towards the center of the primary or the secondary winding as shown in Fig. 11. In this case, the field lines (purple dashed line shown in Fig. 11(b)) are not straight in the winding region as assumed in the proposed model and are hence longer than the assumed field lines (red dashed line shown in Fig. 11(b)), what results in a lower field amplitude in the simulation.
Since there are 4 methods for calculating the magnetic field and 2 methods for calculating the G-factor, 4 × 2 = 8 different methods plus 1 method, which is introduced in [16] can be used for the loss calculation. In Fig. 10(a2), (a3), (b2), and (b3) seven methods are compared to FEM results. The 1D loss model presented in [5] is used for comparison, as it has a relatively smaller error than other 1D loss models. Moreover, the proposed model based on the hyperbolic functions ("New_H") is more accurate than the model based on the Bessel functions ("New_B") (see Fig. 14), thus only the method "New_H" is used. In Fig. 10, the methods "AVGB", "AVGH", "HFEMB", and "HFEMH" are the loss models based on 1D average field and Bessel/hyperbolic functions and the loss model based on H FEM and Bessel/hyperbolic functions. The last model "SFD" (squared field derivative) presented in [16] is also used for comparison. As can be seen in Fig. 10(a2) and (b2), the method with higher accuracy of field estimation exhibit also a better accuracy of loss calculation for most layers.
To compare the accuracy of different methods for loss calculation more precisely, the error for the total losses of transformer design 5 and 10 are shown in Fig. 10(a3) and (b3). For both transformer designs, the proposed model has the best accuracy (< 2%), and the method "AVGB" is the second most accurate one with an accuracy between 3% and 4%. Although the semi-numerical methods "HFEMB", "HFEMH", and "SFD" are based on the most accurate field results, the accuracy of these methods is worse than methods "New_H" and "AVGB". The possible reasons for this can be concluded in the following 2 reasons: r The losses of strands in the same layer are not identical. r The G-factor is not accurate or the (1) is not rigorous for calculating 2D losses. The first reason can be conducted from Fig. 1, as the field amplitude through each strand in the same layer is different, the losses of the strands, which are subject to different field amplitudes, are also different. Therefore, using a uniform field (spatial average field) for all strands in the same layer to calculate the losses results in inaccuracy. To verify the second reason, the G factors based on Bessel G B and hyperbolic G H functions and used in method "SFD" G SFD (see [39]), which are all derived based on a 1D field assumption and are named 1D G factors here, are compared to the G factor calculated from FEM (G FEM ). Note that the G factor G SFD is equal to the losses P p, SFD divided by H 2 FEM , which is not the same as it is defined in [39]. Moreover, the G FEM can be calculated by (41), if P Skin is assumed to be orthogonal to P p , where P t,FEM is the total losses of a strand. As LWs is usually made of the same strands and all 1D G factors depend only on the frequency, material and dimension of the strand, the 1D G factors is the same for the strands in the same winding. In Fig. 12, the G factors for the 333 strands in the primary  Fig. 1(b) are shown in (a) and (b).
winding of transformer design 5 are shown. As G FEM is shown in a grey line with 20% opacity, the darker lines or area contain multiple curves of G FEM . The G factor G FEM varies significantly for different strands, nevertheless, around 90% of the curves of G FEM are located within ±7.5% range of the average value (G FEM avg ) of G FEM , which is highlighted in yellow in Fig. 12. Moreover, the curves of the 1D G factors, which are located on the top of the profile of the G FEM curves, are very close to each other and are slightly higher than the curve of G FEM avg . The remaining 10% strands, of which the curves of G FEM are below the curve of G FEM avg , are located mainly around the left edge of the primary winding, which are highlighted in black in Fig. 1. Note that the transformer design 5 is the same as the transformer shown in Fig. 1, and strands 1 and 2, which are highlighted in magenta in Fig. 1(b), belong also to this 10% strands. These strands are close to the center core leg, thus a certain share of the magnetic field is bent and enters directly into the center core leg rather than traveling straight and entering the core yoke. In this case, the magnetic fields through these 10% strands have more 2D shape, which means the magnetic field through the strands are curved with changing amplitude rather than straight with a constant amplitude. To compare this 1D and 2D field, the current density and field line of 4 strands, which are shown in magenta in Fig. 1), are calculated by FEM, as shown in Fig. 12. The curvature of the field line decreases from strands 1 to 4, and the magnetic field tends to converge to a straight line. Correspondingly, G FEM increases from strand 1 to 4, as the smaller the curvature of the field line is, the closer is the magnetic field to 1D. The curve of G FEM of strand 4 is aligned with the 1D G factors, as the field is almost 1D. Moreover, (1) is derived based on the assumption that the magnetic field is 1D and perpendicular to the conductor. Thus, (1) is not accurate for calculating the losses P p for a round conductor subjected to a 2D magnetic field. Nevertheless, this overestimated 1D G factor used in the proposed model is kind of compensated by the underestimated magnetic field through those winding layers, which are next to the center/outer core leg. This results in accurate loss results for layer P1 and S2 in Fig. 10(a2) and for layer P1 and S3 in Fig. 10(b2). Based on the analysis above, it can be concluded that the accuracy of the loss calculation can be improved by replacing the 1D magnetic field in the 1D model [5] by other more accurate field values. However, more accurate field values do not results in more accurate losses because of the inaccurate 1D G factors and unequal losses in the different strands in the same layer. Moreover, the strands subjected to the strong 2D field are close to the center and outer core legs, thus, the error of the semi-numerical models increases with the loss contribution of the strands close to the core legs. As this article aims to provide an accurate analytical loss model, the semi-numerical models ("HFEMB", "HFEMH", and "SFD"), which uses the magnetic field calculated by FEM are not considered in the following comparison. As it is difficult to identify, which analytical model is more accurate, it is necessary to compare them for different transformer designs. In the following comparison, 12 transformer designs, which cover a wide range of core sizes (from E 20/20/5 to E 90/50/16), and winding arrangements with different window to winding height ratio are used as shown in Fig. 14. Because of the symmetry, only a single core window is simulated. Since the transformers are designed to evaluate the loss models, they are not optimized.
In order to compare different models in the same figure, an error bar figure is used. As can be seen in Fig. 13, the error bar figure is converted from a 3D figure by projecting the 3D figure on the plane C1 and extracting the peak, the median, and the minimal error for each transformer. The error bar figure shows not only the losses changing with the window to the average winding height ratio h w /h avg , but also with frequency, which is indicated by the penetration ratio = d s /δ.
In Fig. 14 the new models based on the hyperbolic ("New_H") and the Bessel functions ("New_B"), as well as models based on the 1D average field ("AVGH" and "AVGB") are compared to the 1D field loss model [5], which is widely used for calculating the HF-losses of LW. The error displayed in Fig. 14 is calculated by Error = P model −P FEM P FEM × 100%. The geometry used for the loss calculation of each transformer case is also displayed. It can be seen that the error of the proposed model based on the hyperbolic function ("New_H") is lower than all other methods for every transformer design except for transformer design 9, where the error of method "AVGB" is slightly lower than method "New_H". The error of the 1D field loss model increases with h w /h avg , and it is up to -30% when h w /h avg is bigger than 1.3, whereas the error amplitude of the proposed method "New_H" is within 8.5%. By comparing the average value of the maximum error E avg,max of each transformer design as shown in the legend of Fig. 14, the proposed method "New_H" with E avg,max = 2.69% is identified to be the most accurate one. Although the average error E avg,max of the method "AVGB" (4.32%) is bigger than the proposed method "New_H", the method of "AVGB" is easier and faster than the proposed model "New_H". Therefore, both methods are selected for the experimental validation.
Besides the accuracy, the computational time is also compared. Table 1 shows the average computational time of 3 calculation procedures, including the 1D field loss model, the proposed model with and without using the algorithm presented in this article. The computational time is the average time of calculating losses at 12 frequencies based on 12 transformer cases shown in Fig. 14 and is measured by using the Matlab function 'timeit' on a laptop with a CPU i7-8665 U. Thanks to the algorithm, the computational time of the new model given in Table 1 is around 275 µs, which is 30 times less than it without using the algorithm. Note that no parallel computing is used in all functions in the calculation, hence, only a single core is used. By integrating this new model and the algorithm into converter optimization routines, a better and faster design can be achieved.

B. EXPERIMENTAL VALIDATION
Besides the numerical validation, experiments have been also performed to validate the proposed model. The validation is done by comparing the AC resistance R AC of different transformers calculated by the proposed model to the one measured with an impedance analyzer. A standard short circuit measurement is performed, i.e. the primary or the secondary winding is shorted during the measurement as shown in Fig. 15. Usually, the transformer designed for converters has very low winding resistance to improve the efficiency. However, such transformers are not suitable for accurate impedance measurements, since the impedance analyser can only measure the impedance accurately in a certain range, which is normally larger than the impedance of such transformers. Therefore, special measurement setup and transformer designs are required and discussed in this section.

1) MEASUREMENT SETUP
To accurately measure the AC resistance of the windings, the following 2 issues must be considered: r The impact of the core losses on the winding losses must be kept at a very low level.
r The resistance of the winding must be designed in an appropriate range, where the impedance analyser can measure the resistance accurately. To minimize the core losses, core materials with low loss density, such as N87, are used in the transformer design. Moreover, the core losses consist mainly of hysteresis and eddy current losses, which increase with frequency and the flux density in the core. To further reduce the core losses, the amplitude and frequency of the current/voltage excitation must be kept at a low level. According to the data sheet of material N87 [40], the imaginary part of the permeability, which is related to the core losses, increases rapidly after 200 kHz. Therefore, the maximum frequency of the excitation in the performed measurements is below 200 kHz, except for transformer a3). The used frequency ranges are given in Table 2.
In addition to the excitation frequency, the excitation amplitude must be kept also low. However, excessively small excitation amplitude reduces the measurement accuracy of the impedance analyser. Based on the manual of the impedance analyser (E4990 A), if other configurations such as measurement frequency, measurement time and so on, are set to fix values, the impact of the excitation/oscillator level (V osc ) on the overall measurement accuracy can be calculated. As shown in Fig. 16(a), the oscillator level with the lowest error is at 100 mV and between 200 mV and 500 mV. To keep the core loss level low, 100 mV oscillator level is used in the measurement. As the impedance analyser can only accurately measure the impedance in a certain range, the impedance needs to be in  the range of (300 m -3 M ) and the measurement frequency needs to be in the range of 20 Hz-10 MHz to achieve less than 1% accuracy, based on the accuracy map of impedance measurements with four terminal pair port given in [41]. However, this impedance range is usually hard to achieve, as at very low frequencies the impedance of the transformer is almost equal to the DC resistance of the winding, which is usually much lower than 300 m , where the measurement accuracy is not well specified. Moreover, the accuracy map given in [41] is only for the 4-terminal pair port. In the measurement, a fixture for connecting the transformer to the impedance analyser must be used. The extra measurement error caused by the fixture, which is not included in the accuracy map in [41], must be also considered. To determine the measurement accuracy for the low range of impedance values, the error caused by the 4-terminal pair port, the fixture and the combined error are calculated as a function of impedance values given in Fig. 16(b). As can be seen in Fig. 16(b), the combined error is below 5% at Z > 90 m , hence, the DC resistance of windings is designed around or above 90 m . Furthermore, the error drops rapidly with frequency, as the impedance of the transformer and the measurement accuracy of the impedance analyser increase both with the frequency. Measurement errors are visible as stochastic fluctuations in the impedance curve as shown in Fig. 18(a2). Usually, large combined error values occur only in the low frequency range ( f < 1 kHz see Fig. 18(a2)) and decrease to less than 1% at above 1 kHz.

2) TRANSFORMER DESIGN
The transformer design here is based on the constrains given by the measurement setup rather than the technical requirements of a converter design. Therefore, high power density and efficiency are not important, but the proper DC resistance of the winding and low core loss. Besides those 2 design criteria, the following points must be considered: 1) The designed transformers need to cover a wide range of the window to winding height ratio. 2) The HF-effect of the winding must be visible under 200 kHz. 3) Variable winding arrangements and LW types need to be implemented in some transformers. 4) Transformers need to be sufficiently different in terms of geometry. Based on these requirements, 5 different transformers with window to winding height ratios from 1.1 to 1.8 are designed, which satisfies requirements 1) and 4). The detailed design parameters are given in Table 2. To obtain a significant HF-effect of LW under 200 kHz, the strand diameter must be larger than the skin depth of copper at 200 kHz which is 0.15 mm. Furthermore, to achieve high DC resistance of   Tables 2 and 3. In addition, 3D printed bobbins are used as shown in Fig. 17 to achieve different winding arrangements of the primary and the secondary winding.

3) MEASUREMENT RESULTS
In total, the impedance of 5 different transformer prototypes (see Fig. 18) is measured. The AC resistance is calculated by the proposed model "New_H", "AVGB" and the 1D field loss model for comparison as shown in Fig. 18(a1)-(e1). To indicate the exact error of the analytical models compared to measured values, the error of these 3 models is also given in the Fig. 18(a2)-(e2).
The bobbin offered by the ferrite core manufacturer is used in the transformer cases (a3) and (b3) and the window is fully occupied by the winding to achieve a minimal window to winding height ratio. The secondary winding is wound directly on the primary winding. As can be seen in Fig. 18(a2), the error of all analytical models is within ±10%. The error amplitude (7.5%) of proposed model AVG_H is slightly bigger than the error amplitude of model "AVGB" but is smaller than the error amplitude (8.5%) of the 1D field model. The proposed model overestimates the total winding losses, because in the proposed model AVG_H the winding is assumed to be enclosed in the core window. However, in the transformer shown in Fig. 18(a3), E shaped cores are used and only approximately half of the winding sections is enclosed in the core. The other half of winding sections is located outside the core window, where the magnetic field and the winding losses are lower than them inside the core winding.
To evaluate the proposed model "New_H" properly, a transformer with two pot cores is built as shown in Fig. 18(b3). This ensures that the winding sections are completely enclosed in the core window. In the setups, the bobbin supplied by the core manufacturer is used. As expected, the proposed model "New_H" has the lowest error (< 1.8%) in the whole frequency range, as shown in Fig. 18(b2). The error of the model "AVGB" is slightly bigger than the proposed model "New_H". However, the error of the 1D field model is up to −13%, which is much higher than the error of the presented 2 models.
As can be seen in Fig. 18(c3), a single LW is used for both the primary and the secondary winding. The LW is first wound in two layers for the primary winding and then bent over, to be wound in the opposite direction to form the two layers of the secondary winding. This technique builds a perfect series opposing connection. As the turns ratio is one, the flux in the core generated by the primary winding is exactly canceled by the secondary winding. Therefore, the effect of core loss on the measurement is significantly reduced. The resistance curve calculated by the proposed model matches the measurement results perfectly in the complete frequency range, as shown in Fig. 18(c1). In Fig. 18(c2), the error of the proposed model is almost equal to 0, whereas the error of method "AVGB" and the 1D field loss model is up to −6.43% and −20%.
In Fig. 18(d1) and (e1) the AC resistances for transformers with a large window to winding height ratio (1.5 and 1.8) are calculated and compared. In these two transformers, two 3D printed bobbins (see Fig. 17) are used for the primary and the secondary windings with different winding heights and LW types in each transformer. In Fig. 18(d2), the proposed model predicts the AC resistance with less than 10% error in the whole frequency range, whereas the error of method "AVGB" is up to −14.7%, which is twice as much as the peak error of the proposed model "New_H" (up to −7.2%) for transformer (d3). In Fig. 18(e3), the error of method "AVGB" (up to −10.2%) is slightly bigger than the error of the proposed method "New_H" (up to −8.7%) for transformer (e3). By contrast, the error of 1D field loss model is up to −35% and −45%, which is much higher than the error of the other 2 models for transformer (d3) and (e3).
Based on both the numerical validation and the measurement results, it can be concluded that, only the proposed model "New_H" can accurately predict the winding losses with less than 10% error for transformers with a window to winding height ratio smaller than 1.5 (h w /h avg < 1.5). The   TABLE 4. Comparison of Different Models, the Semi-Numerical Models Refer to the Model "SFD", "HFEMB", and "HFEMH" model "AVGB" is comparably accurate as the proposed model "New_H" in the range h w /h avg < 1.2, however, in the range h w /h avg ≥ 1.2 the error of model "AVGB" (up to −14.7%) becomes significantly larger than the proposed model. As model "AVGB" is similarly accurate but much simpler than the proposed model "New_H" in range h w /h avg < 1.2, thus, model "AVGB" is recommended for transformers with a window to winding height ratio smaller than 1.2. For transformers with a window to winding height ratio bigger than 1.2, only the proposed model "New_H" achieves a good accuracy of less than 10%. Therefore, the proposed model is recommended for transformers, which have a window to winding height ratio bigger than 1.2. A comparison of different models is given in Table 4.

VI. CONCLUSION
Due to isolation requirements and mechanical constraints, the winding height is usually lower than the window height and therefore 2D magnetic fields occur in the window, what results in substantial error for conventional 1D field loss models [5], [6], [7], [8], [9], [10]. To improve the accuracy of the loss calculation, a new loss model, that calculates the 2D magnetic field in the winding region and modifies the 1D field loss models with this magnetic field, is proposed in this article. Furthermore, an algorithm, that converts the integral of the magnetic field for each layer into a matrix system, substantially reduces the calculation time. The model ("New_H") and the algorithm proposed in this article offer a good compromise between accuracy and computational time. By comparing to 3 semi-numerical methods and 3 other analytical models including 2 average field model and a 1D field model from [5], the proposed model is more accurate than all other analytical models for the most transformer designs, especially for the transformer designs with a large window to winding height ratio (h w /h avg > 1.2). The proposed model exhibits even a better accuracy than these 3 semi-numerical methods. The numerical and experimental validations based on multiple transformer cases indicate that the proposed model can estimate the winding losses with less than 10% error within 300 µs on a personal laptop, whereas the conventional 1D field loss model [5] has an error up to −45%. With these promising results, the proposed model will improve the optimization of the magnetic component as well as the converter designs. Furthermore, the proposed model can be easily extended to interleaved windings. As shown in Fig. 19, an exemplary interleaved winding with a PSPSP arrangement is shown. Based on the magnetic field distribution shown in Fig. 19(a), the winding window can be divided into 4 blocks (B1-B4), where the field starts always at zero and increase to a peak value and goes back to zero. Therefore, the proposed model can be applied to each winding block by dividing the repetitive winding block into 4 regions (R1-R4) as shown in Fig. 19(b). Fortunately, this interleaved winding example is symmetric to the midline (purple dashed line) in Fig. 19. Therefore, only the field value in blocks B1 and B2 needs to be calculated, as the field value in blocks B3 and B4 is equal to that in blocks B1 and B2.