Impedance Modeling Oriented Toward the Early Prediction of High-Frequency Response for Permanent Magnet Synchronous Machines

This article presents and validates a method for early estimation of the high-frequency impedance response for a real permanent magnet synchronous machine. The main contribution of this work lies in the simplicity of the data required to tackle the prediction, which do not require any experimental tuning or 3-D simulations. The required inputs for the method are the 2-D cross section of the machine and the material characteristics of the conductors, insulation, and iron laminations. An approach using the finite element method coupled with a circuit simulator is utilized to obtain the machine common-mode and differential-mode impedances over a broad frequency range. Experimental measurements are performed up to the resonance band to validate the methodology. In addition, the influence of the machine rotor and housing on the impedance curves is experimentally assessed. A machine with concentrated winding is used as a case study machine to test the model.

appearance of wide bandgap semiconductors, high-frequency phenomena are attracting increased research attention. Following this technology trend, the assessment of the high-frequency response throughout the drive design process becomes of capital importance.
Several predictive models for EMI and electrical machine failure were developed over the last years [3]- [5]. The rotating electrical machine model plays a key role in the full drive simulation. These models represent the machine impedance in different propagation modes over a broad frequency range, usually starting from tenths of kHz until the MHz range. The two main trends for motor modeling are the data-based models obtained from measured impedance curves and the pure numerical models using geometry, material characteristics, and machine configuration parameters.
The main types of data-based models are the behavioral or black-box models and the physics-based circuits. In the first type, fitting and synthesis techniques, such as vector fitting [6], are utilized to map the measured impedance of the machine [7]- [9]. Physics-based circuits attempt to describe the real highfrequency behavior of the machine by using passive circuit elements. The resonant frequencies of the measured impedances are used to compute the values of the different passive elements in the circuit. Some examples of physics-based circuit models can be found in [10]- [13]. Data-based models need a real machine and measurement equipment to obtain the impedance of the different propagation modes. For this reason, it is not a valid methodology for early-stage design when these resources are not yet available.
Numerical models can overcome data-based models limitations since no machine prototype and measurement equipment are needed. However, the complexity of the high-frequency phenomena makes this method challenging in terms of computational power and accuracy. The most utilized approach is to obtain capacitive, inductive, and resistive parameters by using different finite element method (FEM) solvers and to assemble them in a distributed parameter equivalent circuit of the machine winding. In [14], a 3-D model was built to take into account eddy currents within single motor laminations with arbitrary iron permeability, which was proven to affect the inductive and resistive parameters. Mihaila et al. [15] developed a full machine model to predict interturn voltage spikes with a SPICE simulated circuit, where only the capacitive couplings were precisely computed by using the FEM. In [16], a simplified coil was modeled in detail considering skin and proximity effects in the conductors where the iron core was modeled in 2-D assuming that flux do not penetrate the core at high frequencies. In [17], a 3-D model was developed focusing on the magnetodynamic solver for resistance and inductance computation, but no experimental validation was provided. In [18], both common-mode (CM) and differentialmode (DM) impedances of a full induction machine was modeled, including frequency-dependent permeability, in order to have improved accuracy. Measured data were utilized to estimate the eddy-current-associated losses, which adjusts the damping of the simulated curves. Boucenna et al. [19] developed a model for (CM) impedance prediction, where dielectric loss was accounted for by performing an experimental tuning. This model was tested and validated on a simplified stator sample. Radja et al. [20] focused on the usage of FEM to predict the transient voltage distribution in machine windings. The method considers skin and proximity effects in single conductors but no data about the influence of laminated iron cores was found. In [21], voltage stress and current distribution prediction in individual conductors by utilizing FEM methods were targeted. Randomness of the conductor distribution was studied, and the model was fully validated on a simplified stator sample. End-winding sections and laminated iron core behavior were neglected in this work. In [22], a simplified method to predict machine terminal overvoltages by utilizing a 2-D FEM approach, with no experimental validation, was described. Ryu et al. [23] developed a method for induction machine CM impedance estimation, including the study of the variation of some design parameters. Results about the extraction of the different multiconductor transmission line parameters were not presented and discussed in this work. In [24], a completely numerical methodology to predict both CM and DM impedances was introduced for a case study permanent magnet synchronous machine (PMSM). However, the impedance measurement did not cover the resonant frequency band.
This article presents and validates a completely numerical approach to predict qualitative trends of the high-frequency impedance response for a PMSM at early design stages. Simplicity is the main asset of the method, since computationally expensive 3-D simulations and experimental tuning are avoided. In this way, the parameters of the equivalent circuit are estimated only using 2-D FEM methods, including the complex electromagnetic behavior of laminated iron cores. The machine selected to implement the method is a 8/12 PMSM with concentrated winding and double layer slots. Experimental measurements are performed using machine samples with different constructional features to reduce the validation uncertainties. The rest of this article is organized as follows. Section II explains the modeling method, the initial assumptions, and the input data. Section III assess the FEM analysis to obtain the equivalent circuit parameters and the equivalent circuit assembly. Section IV analyzes the experimental measurements and discusses the main outcomes of the study. Finally, Section V concludes this article.

II. MODELING METHOD DESCRIPTION
The high-frequency model described in this study is based on the coupling of different FEM solvers with a SPICE circuit simulator. Each single conductor inside the machine slots is modeled to take into account interturn couplings. The equivalent circuit is built following multiconductor transmission line principles, where each conductor in the winding is represented by an equivalent circuit. Fig. 1 shows a block diagram of the modeling method.
Transmission lines are defined by their per-unit length parameters, which are the resistance, inductance, capacitance, and conductance. The dielectric material conductance is neglected since it is affected by a wide variety of parameters such as temperature, moisture, ageing, or pressure. It is not possible to find tabulated data about this parameter, and it can be estimated only by using experimental methods [19]. The capacitances are computed by using an electrostatic FEM solver while the inductances and resistances use a magnetodynamic approach. This is a common procedure to obtain RLC parameters, as shown in [20], [21], and [25].
It is assumed that the conductors are electrically short. A circuit can be considered electrically small when its largest spatial dimension is smaller than one-tenth of a wavelength.  This assumption implies that each single conductor can be represented by a lumped equivalent circuit. It is also necessary to avoid restrictive computational efforts when assembling the full machine equivalent circuit. The pi-equivalent circuit type is selected to ease the circuit practical implementation. This circuit supports the single representation of the frequency-dependent resistances and inductances for each conductor. Fig. 2 shows the pi-equivalent lumped parameter circuit and its RLC parameters. Note that resistances and inductances are embedded inside the frequency-dependent impedance.
The elements not included in the 2-D cross section, such as end winding, housing, and busbars are neglected. The consideration of the end-winding sections would imply the usage of 3-D simulations as shown in [19] and [26], and thus, it is avoided in the present modeling approach at a accuracy price. In addition, the end-winding sections for concentrated winding are normally very short compared with the distributed topology. Fig. 3 shows the machine cross section, and Fig. 4 shows a single tooth geometry. The conductor distribution is determined according to the distance from the tooth symmetry axis to the position of last conductor in a half slot. These two lengths are highlighted in Fig. 4 as h 1 and h 2 . In addition, conductor diameter, slot wall insulation thickness, and fill factor are normally known or assumed, which restricts the location of conductors in the case of randomly wound coils [27]. Conductor location is precisely known for the machine under study. A star connection with one  branch per-phase and four coils connected in series for each branch is utilized. All material characteristics are assumed to be isotropic over the 2-D surfaces of the geometry. Table I summarizes the most important material and geometric characteristics of the machine under study.

A. Electrostatic FE Analysis
An electrostatic solver is used to obtain parasitic slot capacitances [28]. The governing equation for the electrostatic solver is given by [29] −ε∇ 2 where ε is the permittivity of the homogeneous region, V is the electric potential, and ρ is the charge density. The capacitance extraction is addressed by using the surface charge in the modeled conductors. The slot walls are considered as the reference ground. Only one conductor is excited at a  time with a unitary voltage, inducing a surface charge density in the nonexcited slot elements. The different capacitances are obtained by using the following expression: where q is the surface charge density on the different slot elements.
The capacitance values are organized in a matrix, where the off-diagonal terms represent the interturn capacitances, and the diagonal entries are the turn-to-ground couplings. Fig. 5 shows the capacitance matrix for a machine slot, and Fig. 6 shows the equipotential lines and a voltage density contour plot when a single conductor is excited in the slot. From the figures it is possible to observe that the strongest interturn capacitive coupling occurs between neighboring conductors. For turn-to-ground capacitances, the strongest coupling takes place in conductors directly facing the slot walls.

B. Magnetodynamic FE Analysis
In this section, the computation of frequency-dependent resistances and inductances is presented. By using FEM, it is possible to consider skin and proximity effects in the conductors inside the slots. In addition, the flux penetration in the iron core is taken into account. The governing equation of the magnetodynamic solver is given by [29] with A the complex magnetic vector potential, σ is the electrical conductivity,Ĵ src the phasor transform of the applied source current, B the flux density, and μ the magnetic permeability. The problem is solved in the frequency domain, where all ac quantities and A are expressed in phasorial notation. 1) Iron Core Model: An iron core model taking into account eddy currents flowing within single laminations is required. These currents will shield the main magnetic flux, and at a certain frequency, this flux will be mainly confined inside the slot. This effect will influence the value of the conductor impedances over a broad frequency range.
The iron core model is based on the correlation of the laminated anisotropic block to a solid isotropic nonconducting material with a frequency-dependent magnetic permeability [30]. This frequency-dependent complex permeability is derived from the 1-D magnetic field solution along the axial direction of a lamination. Fig. 7 shows the geometry of a single lamination and the magnetic field and current density theoretical distribution. The expression for the magnetic field along the Cartesian z-axis is given by [31] where H x (z) is the magnetic field amplitude in x-direction along the z-axis, H 0 is the surface magnetic field amplitude in the edges of the lamination when z = b, k is the complex term defined as (1 − j)/δ, and δ is the penetration depth δ = 1 πf σμ . The effective complex permeability is defined as the average magnetic flux density in the lamination normalized to the surface magnetic field. By combining this definition with (4), the complex frequency-dependent effective permeability is found as follows: The relative permeability in low frequency (μ r,dc ) depends on the hysteresis loop characteristic of the material. In PMSMs, the magnetic field created by the rotor magnets will excite the magnetic domains of the soft-magnetic material. The lowamplitude high-frequency signals will follow minor hysteresis loops around the instantaneous working point, which moves on the main loop at the basic frequency [32]. It is assumed that the iron core behaves linearly in the initial part of the BH curve instead of hysteretic. In this way, it is possible to use the low-frequency initial magnetization BH curve provided in the material data sheet. Thus, the low-frequency permeability is estimated by using the average magnetic flux density in the stator at a fixed rotor position and the nonlinear BH curve where H| B avg the magnetic field strength for the obtained averaged magnetic flux density. The term B avg is obtained by averaging both B-field components (i.e., x-and y-components) over the stator mesh nodes during a magnetostatic simulation at a fixed rotor position. For other types of machines without a dc field in the rotor, the low-frequency magnetic permeability equals to the initial permeability of the material since the input signals of the standstill frequency sweep have very small amplitude.

2) Iron Core Model Validation:
A simple electromagnetic device is built to validate the iron core model. The model depth is equal to the lamination thickness. Fig. 8 shows the geometry for the simplified model, and Table II tabulates the model input  data. Three modeling approaches are implemented. The first one utilizes an extruded 3-D mesh to properly map the magnetic field and current density distribution over a single lamination (3-D). The second method is a 2-D approach with an adjusted iron core mesh to account for the iron core penetration depth (2-D   μ dc constant). The third approach uses the homogenized iron core model with complex frequency-dependent permeability (2-D μ eff (f )). Fig. 9 shows the extruded 3-D mesh.
The results from the model validation are presented in Fig. 10. The frequencies for the study are set between 1 and 100 kHz, and only ten frequency points are defined. The iron core modeled with μ eff (f ) closely follows the trend from the 3-D approach, while the detailed 2-D mesh model with constant lowfrequency permeability overestimates the shielding effect in the slot.

3) Slot Impedances Calculation:
The impedance for each conductor inside the slot is computed by separately exciting them with a real-valued unitary current. The voltage drop in the excited and nonexcited conductors represent the self-and mutual impedances, respectively. Mutual effects between conductors from different slots are not considered since the penetration depth at high frequencies is smaller than the tooth width. The mesh in the conductors is adapted to the skin depth for each simulation frequency. The contour plot in Fig. 11 shows the absolute value of the current density with real flux density lines at 10 kHz. The real and imaginary terms of the impedance are organized in matrices with one matrix for each simulation frequency. Diagonal terms represent the self-resistance and inductance while off-diagonal entries are the mutual couplings. Fig. 12 shows the resistance and inductance matrices corresponding to 10 kHz, where the strongest couplings are easily observable.

C. Equivalent Circuit Assembly
The frequency of the self-and mutual resistances and inductances is implemented by using the Laplace transform together with dependent voltage and current sources in the circuit simulator. Each impedance can be represented by using a function of the shape y(x) = ax b + c. This expression is selected as the outcome from a tradeoff between accuracy, simplicity, and method compliance. Thus, we can describe the frequency-dependent behavior by using where p n are the parameters used to build each impedance curve over frequency. These parameters are obtained by using a curve fitting algorithm for each impedance response. Self-impedances are represented by voltage-controlled current sources. The source current is defined by a transfer function and the voltage drop in the source according to ⎧ ⎪ ⎨ ⎪ ⎩ (8) where s is the Laplace variable, and G(s) is the transfer function between the source voltage and current. The mutual-impedance terms are represented by behavioral voltage sources. These sources emulate the voltage drop caused by the mutual terms by multiplying the current flowing through one conductor and the transfer function of the mutual impedance M (s). Fig. 13 shows the SPICE representation of a n number of pi-equivalent circuits, and Fig. 14 shows the comparison between FEM-obtained data and the SPICE-implemented element for single conductor impedances. The optimized fitting algorithm show discrepancies for the mutual resistances in the 10-100 kHz frequency band. The frequency-dependent parameters and capacitances are assembled forming pi-equivalent circuits for each conductor inside the slots. Interphase terms are included by lumping the mutual impedances and capacitances. A general script is utilized to build the equivalent circuit connections and to integrate all the conductor parameters. This is achieved by automatically writing SPICE "netlists," which allow to simulate a large number of circuit elements without using the software graphic interface.

IV. EXPERIMENTAL MEASUREMENT AND RESULTS DISCUSSION
This section describes the experimental setup and the different machine samples utilized for model validation. The experimental validation is performed by small-signal frequency sweep in standstill conditions. Three motor samples are utilized to experimentally evaluate the influence of the rotor insertion and the housing in the impedance response. The three samples are a machine stator, the same stator inserted into the machine housing without rotor and shaft, and the full machine. In this way, validation uncertainties are reduced by considering different constructive features. Fig. 15 shows the machine samples utilized in the study.
A precision impedance analyzer (Agilent 4294 A) with frequency band 40 Hz-110 MHZ is utilized to obtain the motor impedances. A test fixture for leaded components compatibility is used (16047E). Auxiliary connection cables are needed to properly interface the motor with the impedance analyzer. Accuracy starts to decrease beyond 10 MHz due to the parasitics of the connection leads even if calibration under short-circuit and open-circuit conditions is performed for each impedance measurement. Due to this lack of accuracy, the frequency range for the test is defined between 10 kHz and 10 MHz.
The connections for impedance measurement are implemented according to the two main propagation modes. The DM   connection is made between two phases. For the CM impedance, the connection is made between the machine housing and the three phases in short circuit. The influence of the stator grounding point is tested by changing its connection. Four different points along the axial and circumferential axis of the stator are chosen, as shown in Fig. 16. Fig. 17 shows the results for the grounding test. From this result, it is possible to confirm that the grounding position is not relevant for the measurements in the machine under study. In addition, the lack of accuracy beyond 30 MHz is visible, when the parasitic elements of the connection cables are dominant over the main machine core. The comparison between the experimental measurement of the three machine samples and the simulation results is shown in Fig. 18. The simulation results follow the same trends as the measured impedances in both propagation modes. The impedance amplitude of the simulation results is located in the same range as the experimental measurements and the resonance and antiresonance points are closely approached. The first CM resonance point ( 1 ) shows a discrepancy of 500 kHz, and the antiresonance point ( 2 ) shows a lower difference of 120 kHz. On the DM impedance curves, the main resonance point ( 3 ) presents a close agreement with a difference of 10 kHz around 570 kHz, while the second and third resonant frequencies ( 4 and 5 ) show a difference of 300 and 230 kHz, respectively. The abovementioned values are indicative and rounded up to the closest multiple of 10 kHz. A lack of damping is observed for both modes, and amplitude differences are detected for both CM and DM in the frequency range before the resonant points.
Regarding the different machine samples, slight variations are observed. The amplitude of the CM curve before the antiresonance point is influenced by the housing. Furthermore, the resonant point is slightly shifted toward a higher frequency with the insertion of the rotor. In the DM curves, the changes are observed mainly between samples with and without housing, with a very reduced influence of the rotor insertion in the studied frequency range. The induced impedance variations due to the rotor insertion and the housing slightly influence the impedance curves but they do not change the impedance behavior in the studied frequency range. This demonstrates that the machine stator represents the most important contribution to the impedance response in both propagating modes for the studied frequency band.
A sensitivity analysis of the model parameters is presented in Fig. 19 to better understand the possible sources of inaccuracies. The lack of damping observed in both impedance curves is associated with the different resistive losses, which are linked to the joule losses in the conductors, iron losses, and dielectric losses. Joule losses are properly taken into account by the utilized FEM method, which models the real field distribution within the slot, skin and proximity effects, and the frequency dependency. Iron losses can be taken into account by utilizing a parallelconnected resistance to the series branches, as described in [18]. This parameter is estimated by using measured impedance data and cannot be properly taken into account by utilizing a 2-D homogenized iron model. Dielectric losses depend on the material characteristics of the insulation. These are also estimated experimentally, as shown in [19]. The numerical estimation of such losses would significantly improve the model accuracy, since damping represents the main drawback of the presented methodology. The amplitude discrepancies over the inductive section of the DM curve are due to an error on the reactance estimation. The causes for this error are linked with the neglected end-winding sections, the consideration of nonhysteretic isotropic BH characteristic in the iron, and the utilization of the averaged magnetic flux density when estimating the lowfrequency relative permeability, which varies for the different coils at a fixed rotor position. Amplitude discrepancies in the initial part of the CM impedance curve (i.e., before the first resonant point) are due to turn-to-ground capacitances estimation errors. The causes for these errors arise from the end-winding contribution and from inaccuracies in the relative permittivity value. The effect of interturn capacitances can be observed in resonant points 2 , 4 , and 5 . Main error sources linked to the interturn capacitances are related to the material characteristics of the conductor insulation. In addition, possible geometric discrepancies associated with the machine manufacturing process and advanced development (e.g., insulation thickness irregularities, tolerances, conductor misplacement, stator welding, etc.) can affect the high-frequency impedance response.

V. CONCLUSION
A method for high-frequency CM and DM machine impedance estimation from an early design perspective was presented. The 2-D machine cross section, machine connections, and material characteristics had been the only utilized input data. The utilized methodology was based on the coupling of a circuit simulator with FEM solvers. The circuit was assembled by using pi-equivalent circuits for each conductor in the slots and the FEM to obtain its RLC components. An iron core model taking into account eddy currents within single laminations was built and validated. Three different machine samples with different constructional features had been considered to experimentally study the effect of the housing and rotor insertion. The comparison between experimental and simulation results showed how the method predicted the main amplitude trends and resonant points for both CM and DM impedances. The experimental comparison between samples proved that the contribution of the housing and the rotor to the high-frequency impedance responses was not critical for the machine under study. Possible sources of inaccuracies were analyzed by performing a sensitivity analysis of the RLC parameters. The damping associated with the dielectric losses and iron losses was the major source of error. The error on the inductive and capacitive terms was mainly associated with model assumptions and uncertainties of the material characteristics. The main benefit of the presented methodology was the estimation of the impedance response in early design stages, avoiding the usage of computationally expensive 3-D simulations, and experimental tuning of the circuit parameters.
Future work on the topic includes the detailed sensitivity analysis of the machine design parameters. This includes the manufacturing of different machine prototypes with different design parameters to experimentally assess the curves variation. The method can be implemented for different types of machines, such as machines with different constructional features (i.e., induction machines, PM-less synchronous machines, and machines with different types of stator winding). In addition, accuracy of the method should be tested for stators with distributed winding since they present a larger end-winding section. The method accuracy can be improved by performing further experimental validation in simple coils surrounded by laminated iron cores to better track model inaccuracies and by developing a reduced end-winding model that enables the fast computation of its RLC parameters. The influence of stator faults (e.g., interturn short circuits) can also be implemented in the model with the aim of predicting and identifying this type of machine faults. Finally, the measurement of rotating machine impedance characteristics on-line instead of in standstill conditions remains an open question that could help on the better prediction of EMIs and machine faults.