A Light Darwin Implementation of Maxwell's Equations to Quantify Resistive, Inductive and Capacitive Couplings in Windings

High operating frequency is an enabler of high key performance indicators, such as increased power density, in electrical machines. The latter enhances the cross-coupling of resistive–inductive–capacitive phenomena in windings, which may lead to significant loss in performance and reliability. The full-wave Maxwell’s equations can be employed to characterize this coupling. To address the frequency instability that arises as a result, a simplification known as the Darwin formulation can be employed, where the wave propagation effects are neglected. Still, this modification is prone to ill-conditioned systems that necessitate intricate pre-conditioning and gauging steps. To overcome these limitations, a fast 2D formulation is derived, which preserves the current continuity conservation along the model depth. This implementation is validated experimentally on a laboratory-scale medium-frequency transformer. The computed impedances for the open-and short-circuit modes of the transformer are validated using measurements and compared with the multi-conductor transmission line model that is widely adopted for the analysis mentioned above. The developed formulation demonstrates a high accuracy and outstanding frequency stability in a wide frequency range, becoming an efficient and computationally light method to investigate the interconnected resistive, inductive, and capacitive effects in windings.


I. INTRODUCTION
To enhance the power density and improve the cost-efficiency of electrical machines, the windings are exposed to a significantly increased operating frequency.Only recently, the higher operating frequency has been enabled at high powers by the semiconductor switches relying on the technology based on silicon carbides (SiC) or gallium nitrides (GaN).The latter is utilized in power converters and other drive systems.A higher fundamental frequency and the large spectrum of frequencies injected in the winding as a result lead to a significant coupling of resistive, inductive, and capacitive (RLC) effects between the turns.Accounting for these effects during the design stage is critical to provide valuable insights into resistive losses, accelerated aging of insulators, common-mode currents of high-power, high dv/dt applications, and control in high-precision systems with aggressive motions, to name only a few.
To investigate the electromagnetic behavior of windings at sufficiently low frequencies, where the RLC interactions are negligible, the electric and magnetic fields can be modeled independently.The electrostatic (ES) and electroquasistatic (EQS) models are widely used, e.g., in a conventional low-frequency power transformer. 1he magnetoquasistatic (MQS) formulation is used to model the magnetization and inductive-resistive effects. 2 However, in a wider frequency range, the coupled RLC effects can be partially accounted for by applying the multi-conductor transmission line (MTL) model, where each winding turn is treated as an independent lump and a corresponding impedance matrix is obtained enabling the quantification of the coupling effect at higher frequencies. 3,4Still, accounting for the interaction between the induced and displacement currents in the winding is limited and could be improved by employing a full-wave Maxwell formulation, for example.However, due to high differences in the material properties, e.g., electric permittivity, conductivity, and magnetic permeability, and high aspect ratio of geometrical dimensions, such as small insulation thickness, conductor cross section, and the large size of winding's and clearance distances, the full-wave formulation faces instability problems at relatively low frequencies and lead to large and ill-conditioned systems.To overcome this issue, Darwin's formulation is derived from fullwave Maxwell's equations. 5The latter neglects the wave propagation effect and considers a negligible displacement current in conducting domains.However, current developments show that it still suffers from the aforementioned problems. 6,7To overcome these issues, various techniques, such as the use of iterative solvers and Lagrange multipliers to restore symmetry in the originally unsymmetrical system, 8 two-step approach, 6,9 model order reduction approaches, 10 the use of alternative gaugings, 11,12 and two-domain continuity gauging strategy, 13 have been proposed.In the latest publications, 14,15 the challenge is addressed by using redundant variables or subspaces as preconditioners for the iterative solvers, resulting in a further increase in computations.Accordingly, the computational intricacy and the inherent three-dimensional (3D) nature are limiting Darwin's formulation to becoming a practical tool in design and validation. 16o reduce the complexity of the 3D formulation and address computational issues, this study adopts an alternative 2D derivation of Maxwell's equations, preserving the dominant effects in the coupling of magnetic and electric potentials.This approach employs geometrical symmetries and invariances of adopted winding designs to simplify the 3D field distribution within 2D cross-sectional slices.This paper is structured as follows: Sec.II provides a brief overview of the quasistatic approximations of Maxwell's equations and the available approaches in investigating RLC components of electrical machines.Section III delves into the physical assumptions and mathematical derivations of the developed formulation.Section IV details the simulation and measurements of the implemented formulation on a laboratory-scale transformer.The final discussion and conclusions are addressed in Sec.V.

II. QUASISTATIC MAXWELL'S EQUATIONS
We begin with Maxwell's equations to describe the behavior of resting electric charge, alternating electric and magnetic fields, and induced currents.In the frequency domain, it reads as where E is the electric field intensity (V/m), D is the electric flux density or the displacement field (C/m 2 ), H is the magnetic field intensity (A/m), B is the induced magnetic flux density (T), and ρ is the electric charge density (A/m 3 ).The constitutive equations are defined to relate the field intensities H and E to the induced flux densities B and D, and the induced current density J (A/m 2 ) using μ, the magnetic permeability, and ε and σ, the electric permittivity and conductivity, respectively, To reduce the computational complexity in solving the magnetic and electric field distributions in Maxwell's equations, B and E are defined in terms of magnetic vector potential (MVP), A, and electric scalar potential (ESP), φ, as Subsequently, by substituting ( 8) and ( 9) in ( 2), the full-wave Maxwell's equation in terms of potentials is derived as The wavelength λ of an electromagnetic wave is related to the speed of light c and frequency f : λ = c/ f .The wavelength is significantly larger than the dimensions of structures involved in the machines operated in the range of a few megahertz.Therefore, the resulting wave propagation, or the radiation term, εω 2 A, is negligible compared to the near-field effects.Considering this assumption, the interaction of electric and magnetic fields through the induced current is still accounted for.Additional approximations could reduce the computational costs when capacitive and inductive effects have negligible mutual effects.One approximation, where capacitive effects are not considered and the goal is to solve magnetic fields and induced currents, is the MQS formulation.In MQS, it is assumed that the displacement (bound) currents are negligible compared to conduction (free) currents: jωε∇φ ≃ 0, which results in Similarly, when the electric field is of interest without the effector of induction, the EQS formulation is adopted to compute the capacitive effects.The derivation of EQS involves applying the divergence operator to (2), which leads to the elimination of ∇ × H. Consequently, the EQS in terms of the ESP is expressed as a conservation law, The independent analysis of MQS and EQS formulations is primarily attributed to the systems operating at relatively low frequencies and encouraged due to constraints in computational resources.In the contemporary definition of power systems, the significance of RLC couplings is obvious.][19][20]

ARTICLE pubs.aip.org/aip/adv
The representation of RLC coupling can be achieved using the MTL method. 21In this approach, the values of R, L, and C are computed through decoupled MQS and EQS formulations, 3,4 as shown in Fig. 2(c).The MTL method lacks accuracy due to the trivial approach adopted in coupling the R, L, and C effects since the contribution of conductive and displacement currents are considered within a concentrated and discontinuous lumped circuit configuration of RLC components.In addition, the computation of these parameters involves multiple computations per frequency and the number of simulations increases exponentially with the number of turns.
The Darwin formulation addresses these disadvantages, by solving the radiation-free Maxwell's equations together with the continuity equation (12), which now includes the induced currents. 22The Darwin model in 3D couples A and φ potentials, Still, the resulting formulation does not have the unique solution for A and suffers from the ill-conditioned matrix.In addition, in a practical application, a 3D representation of such a system yields a large number of degrees of freedom (DOF), requiring additional precautions and resources to solve. 6Thus, gauging and pre-conditioning methods together with iterative solvers are employed. 22n most electrical machines, a significant portion of the winding length typically resides within the winding window or slots, facilitating the 2D geometry fidelity analysis.However, in certain magnetic components, such as transformers, the length of the overhang section of the winding is comparable to that of the windings enclosed by the core within the winding window.This is shown in Figs.1(a) and 1(b), where a medium frequency transformer (MFT) and the 2D cross sections of foil windings inside the transformer window (XY cross section) and outside the winding window (YZ cross section) are illustrated.The inner and outer slices are considered as shown in the derivation, similar to the proposed double 2D method; 23 this time, however, the slices are coupled and solved together.The method accurately captures the behavior of the transformer avoiding expensive 3D simulations.As a validation step, the developed formulation is compared with the MTL model and measurements from the laboratory-scale MFT.

III. DERIVATION OF LIGHT DARWIN FORMULATION
In this section, we derive a computationally light alternative to 3D Darwin's formulation to model the foil windings of a transformer.We divide the 3D model into several subdomains (with overlap) and apply approximations to the subdomains, allowing us to considerably reduce the dimensionality of the model.We consider the inner and outer slices, as shown in Fig. 2, denoted by V I and V O , respectively.The inner slice is a box with a surface S I , a subset of the XY plane, and depth [−r I , r I ], along the z axis.The outer slice is a box with surface S O , a subset of the ZY plane, and depth [−r O , r O ], along the x axis.Furthermore, we consider the turns from the primary and secondary windings W PN and W SM , respectively, as separate subdomains: a part of turn i from the primary winding, W Pi ∩ V I , is a box with surface S IPi ⊂ S I and depth PDEs defined on the inner and outer slices, respectively.The superscripts IPi and OPi are used for the PDE variables defined for the turn i of the primary winding intersecting the inner and outer slices, respectively.The superscripts ISi and OSi are used for the secondary winding, correspondingly.We start with the derivation of the PDE for the ESP in turn i of the primary winding intersecting with the inner slice, W Pi ∩ V I .The derivation for the secondary winding is analogous to the primary winding and is omitted for brevity.The derivation for the intersection with the outer slice is analogous to the inner slice and is omitted for brevity as well.We integrate the conservation of charge over the cross section S IPi , Applying the divergence theorem to the x and y components gives The current density in the normal direction, J IPi ⋅ n, at the boundary ∂S IPi is zero because the other side of the boundary is nonconducting.We account for the penetration of the displacement current into the winding by replacing the normal direction, jωεE IPi ⋅ n, with the displacement current density from the inner slice, jωεE I ⋅ n.Note that the individual quantities ε and E ⋅ n jump at the interface between the winding and the insulator, but the current density in the normal direction should not.We obtain E IPi in the first integral in terms of potentials using (9) and substitute A IPi with A I , the MVP from the inner slice.We apply (7), interchange integral and derivative, and substitute ε for the constant εC inside the conductor.Finally, we define the average potential φ IPi av as This gives the PDE for the average potential φ IPi av of the winding passing through W I on [−r I , r I ], with The derivation for the ESP on [−r O , r O ] in W Pi ∩ V O is similar and omitted for brevity.Combined for all turns, this gives 2N independent PDEs, where N is the number of turns.To connect them, we assume that the winding segment W Pi ∩ V I is directly connected to W Pi ∩ V O , as depicted in Fig. 2. Formally, on the boundaries of each conductor segment, We continue with the derivation of the PDE for the magnetic field on V I , and the derivation for V O is omitted for brevity.The derivation for the field on the outer slice is analogous to the inner slice and is omitted for brevity.Recall that V I is S I × [−r I , r I ].Furthermore, we assume that the z component of magnetic field is zero.This implies that the x and y components of the A I are zero and the z component of A I does not vary in z.This allows us to reduce the Darwin model to 2D.Applying the above to (13a) gives the PDE for A I z on S I , ARTICLE pubs.aip.org/aip/adv with where S I nc is the non-conductive subdomain of S I , Finally, we derive the PDE for the electric field on V I .For this case, we assume that the electric field does not vary in z.This allows us to reduce the Darwin model to two dimensions.We obtain the gradient of the φ I inside a winding from the ESP obtained for the windings at the center, i.e., r = 0.By doing so, we can omit the conducting subdomains from the equation and use the ESP from the PDE for the windings as boundary conditions.Formally, for all i ∈ {0, 1, . . ., n − 1}, Applying the above to the Darwin model, (13b), gives the PDE for the ESP φ I in the non-conducting subdomain S I nc , In summary, the reduced model consists of ( 17) and the equivalent for the outer slice and the secondary winding, (20) and (24), and their equivalents for the outer slice, with boundary conditions (19) Multiplying (20) with test function β I , integrating over S I , and applying integration by parts to the Laplace term, Multiplying (24) with test function ψ I , integrating over S I nc , and applying integration by parts to the Laplace term, For numerical results, we used the following discrete function spaces: piecewise linear functions on S I for A I z and β I and similar for S O , piecewise linear functions on S I nc for φ I and ψ I and similar for S O and quadratic functions on [−r I , r I ] for φ IPi av , ψ IPi , φ ISi , and ψ SIi and similar for S O .The developed formulation is implemented using the open source FEM library Nutils. 24

IV. EXPERIMENTAL AND ANALYTICAL VALIDATION
The developed formulation is validated using a laboratory-scale MFT prototype, shown in Fig. 1(a).The length of the inner slice V I is set equal to half the length of the winding window along the z axis, and the length of the outer slice V O is set such that the total length of two slices is equal to half of the mean turn length of the winding as indicated in Fig. 2(c).The frequency response of the MFT is analyzed using the computed terminal impedance and compared with measurements and results obtained by the MTL method.The MFT is composed of a ferrite (N87) magnetic core with a relative permeability of 2400 and primary and secondary copper foil windings, with σ = 5.998 × 10 7 (S/m), that are stacked using 3D-printed coil formers.The dimensional specifications of the MFT are provided in Table I.
The MVP Az, the norm of displacement currents jεω|E|, and the norm of magnetic flux density distribution |B| calculated using the developed formulation for the inner 2D cross-sectional geometries and the short-circuit mode of operation are illustrated in Fig. 3, at the frequency of 1 Hz and 1 MHz.
To evaluate the accuracy of the developed approach, the calculated short-and open-circuit impedances obtained by the proposed Darwin model and the MTL method are compared with the measurements through a wide range of frequencies and are presented in Fig. 4. To obtain measurements, an impedance analyzer (Agilent 4294A) is used to measure the frequency response of the MFT in open-circuit (O.C.) and short-circuit (S.C.) configurations.It can be seen that the impedance computed with the 2D Darwin method shows a better accuracy compared to the MTL method on account of the cross-coupling of the fields at high frequencies.Specifically, in the S.C. mode, shown in Figs.4(c) and 4(d), the amplitude and the phase are in better agreement with the measurements.The resonance is measured at 1.63 MHz and predicted accurately using the developed approach.Contrarily, the MTL method estimates the resonance frequency at 1.  simulated results continues up to the range of 1 MHz.However, the discrepancy happens at frequencies above 1.5 MHz, where there is a shift in the phase and an inconsistency in the amplitude of the O.C. impedance.By comparing the discrepancy of the MTL results at high frequencies, it can be concluded that the 2D Darwin method shows a better accuracy compared to MTL.In addition, the disagreement on estimating the O.C. impedance calculated using 2D Darwin method could be attributed to the open-circuit operating mode, where the apparent impedance of the primary winding is mainly determined by the magnetizing impedance of the transformer, strongly determined by the core, as illustrated.Thus, the discrepancy is due to either the aspect of field distribution on the 3D cross section of the core and the resultant flux linkage or the ferrite material properties.

V. CONCLUSION
The 2D Darwin formulation is developed based on full-wave Maxwell's equations and effectively leverages the geometric invariance of the system along one axis.Furthermore, the circuit coupling of the two geometries is embedded in the formulation, and the impedance of the MFT under study in short-circuit and open-circuit configurations is calculated and compared with measurements.The comparison across a wide frequency range from 10 kHz to 10 MHz indicates good correspondence between measurements and simulations.Following the numerical implementation and experimental validation, we observe that the proposed approach becomes significantly faster compared to the existing 3D Darwin approaches and eliminates the need for multiple simulations and circuit postprocessing procedures required by MTL.The results exhibit consistency in a wide range of frequencies.The developed formulation facilitates the application of fully coupled magneto-and electroquasistatic formulations in 2D, addressing the computational cost, complexity, and frequency instability associated with the 3D Darwin approximation of Maxwell's equations.Thus, 2D Darwin approximation proves to be an efficient choice for studying phenomena such as voltage distribution across windings and for analyzing surge or common mode currents as well as AC losses in the windings of electrical machines.

FIG. 1 .
FIG. 1.(a) Laboratory-scale MFT, (b) XY and YZ cross sections of the MFT, and (c) circuit representation of the implemented open-circuit and short-circuit tests.

FIG. 2 .
FIG. 2. (a) Illustration of the inner 2D cross sections of the MFT, (b) inner 2D cross sections of the MFT, and (c) top view of the MFT as an explainer for the developed formulation implemented on a three-turn MFT model, including the inner slice V I and the outer slice V O .
and (23) to be solved for φ IPi av , φ ISi av on [−r I , r I ]), φ OPi av , φ OSi av on [−r O , r O ], A I z on S I , A O x on S O , φ I on S I nc , and φ O on S O nc .The Galerkin projection is employed by multiplying (17) with test function ψ IPi , integrating over [−r I , r I ], and applying integration by parts to the Laplace term, − (σ + jωεC) ∬ S IPi dx dy ∫

FIG. 3 .
FIG. 3. The z axis magnetic vector potential Az with (a) f = 1 (Hz) and (b) f = 1 (MHz), the norm of displacement currents jεω|E| with (c) f = 1 (Hz) and (d) f = 1 (MHz), and the norm of magnetic flux density distribution |B| with (e) f = 1 (Hz) and (f) f = 1 (MHz) for the inner cross section area of the modeled MFT.

FIG. 4 .
FIG. 4. Simulated and measured frequency response results, (a) the magnitude of open-circuit impedance, (b) the phase of open-circuit impedance, (c) the magnitude of short-circuit impedance, and (d) the phase of short-circuit impedance.

TABLE I .
Design parameters of the laboratory scale MFT.