Frequency-Domain Analytical Prediction of PWM-Induced Current Harmonics in Slotless Surface-Mounted PM Synchronous Machines

Slotless surface-mounted permanent magnet (PM) synchronous machines with high pole numbers generally have relatively low inductance, which will lead to significant current harmonics under PWM-driven inverters. Current harmonics can cause additional losses and torque ripples, severely impacting temperature rises and mechanical vibrations. Accurate prediction of current harmonics is crucial for machine performance calculation and drive system design. This article proposes a frequency-domain analytical method to predict the current harmonics with high fidelity equivalent circuit models. The harmonic prediction method is validated experimentally under a wide range of operating conditions and improves the prediction accuracy by up to 75% compared to the existing analytical modeling methods. The prediction error using the proposed method is generally maintained below 10%, making it an effective tool for practical applications.

Frequency-Domain Analytical Prediction of PWM-Induced Current Harmonics in Slotless Surface-Mounted PM Synchronous Machines Xiaolong Zhang , Member, IEEE, Kiruba S. Haran , Fellow, IEEE, Philip T. Krein , Life Fellow, IEEE, and Luis J. Garcés Abstract-Slotless surface-mounted permanent magnet (PM) synchronous machines with high pole numbers generally have relatively low inductance, which will lead to significant current harmonics under PWM-driven inverters.Current harmonics can cause additional losses and torque ripples, severely impacting temperature rises and mechanical vibrations.Accurate prediction of current harmonics is crucial for machine performance calculation and drive system design.This article proposes a frequency-domain analytical method to predict the current harmonics with high fidelity equivalent circuit models.The harmonic prediction method is validated experimentally under a wide range of operating conditions and improves the prediction accuracy by up to 75% compared to the existing analytical modeling methods.The prediction error using the proposed method is generally maintained below 10%, making it an effective tool for practical applications.
Index Terms-Permanent magnet (PM) machines, pulse width modulation, carrier harmonics, total harmonic distortion, equivalent circuit.

I. INTRODUCTION
R APID development of electrified transportation has stim- ulated an increasing demand for surface-mounted permanent magnet (PM) synchronous machines (PMSMs) [1], [2].The application of high-pole-number slotless designs has facilitated their continuous development towards high efficiency and high power density.Many research works have been published in recent years on their topology comparisons, electromagnetic and mechanical design, analysis and optimization [3], [4], [5], [6], [7], [8].In particular, drive and control of slotless surfacemounted PMSMs have been a research focus due to their inherent low inductance and current harmonic issues [9], [10], [11], Xiaolong Zhang, Kiruba S. Haran, and Philip T. Krein are with the Department of Electrical and Computer Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801 USA (e-mail: xzhng157@illinois.edu;kharan@illinois.edu;krein@illinois.edu).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TEC.2023 [12], [13].The high-frequency pulse-width modulation (PWM) harmonics introduced by the inverter could impose significant losses on the motor, leading to a decline in efficiency and winding insulation life.The implementation of filters to mitigate these harmonics not only adds extra weight and loss, but also amplifies the risks of reliability.Therefore, accurate analysis and prediction of PWM current harmonics are imperative to enhancing the slotless surface-mounted PMSM system performance.
In [14], [15], [16], precalculated PWM voltage waveforms were imported into time-stepped nonlinear finite element (FE) models to generate predicted non-sinusoidal current waveforms, thus calculating the losses in interior PMSMs.In [17], [18], current control algorithm, PWM generation module, inverter circuit model, and motor finite element model were coupled together to obtain the PWM current harmonics and losses.The advantage of these FE methods is their high modeling fidelity in magnetic saturation and eddy current effects.However, building and solving FE models with solid conductor regions is quite time-consuming, therefore more computationally efficient method is desired.
To reduce computation time, the motor's electromagnetic field finite element model can be replaced with a subdomainbased electromagnetic field analytical model.[19] first used this method for calculating harmonic currents in induction motors and developed both time-domain and frequency-domain approaches.[20] applied the frequency-domain method to calculate harmonic currents and losses in a 2-pole surface-mounted PMSM, which simultaneously calculated the operational inductance.The advantage of using subdomain-based methods is a significant reduction in computation time compared to FE method.However, since the solving process involves numerical iteration to inverse a large matrix constructed from the magnetic field boundary conditions, it remains a numerical rather than an analytical method.
Equivalent circuit modeling is another method to improve computational efficiency, comprising numerical and analytical approaches.Numerical methods include those used by [13], [21], [22], [23].They numerically solve simple motor RL series equivalent circuits to obtain current waveforms and then calculate the harmonics.The advantage of this method is its simplicity and fast solving speed, but it does not consider motor saturation effects and additional losses, leading to large prediction errors.In [23], flux-linkage lookup tables with saturation information were used to replace the constant L parameter to improve accuracy.
Analytical methods of equivalent circuit modeling have been presented in [24], [25], [26], [27], [28].In [24], the motor is simplified as an equivalent L circuit, assuming that the ratio of the time-rate slope of the instantaneous current is constant.[25] used the same assumption but added a look-up table for L to include magnetic saturation and slotting effects in interior PM motors.The advantage of this method is high computational efficiency.However, the constant slope assumption of the model significantly limits its universality.The various lossy components and reduced impedance at high frequencies are absent in the model, which may cause prediction errors.Moreover, the method cannot be extended to high-order circuits, such as drive circuits containing LC filters.
Instead of indirectly calculating spectrum from time-domain circuit solution as in [24]- [25], the amplitude of harmonic currents can be directly calculated in the frequency domain [26], [27], [28].Frequency domain method is more straightforward because expressed in phasors, the harmonic current can be calculated by dividing harmonic voltage over the impedance.In [26], Park developed a frequency-domain equivalent circuit model, including LC filters, for synchronous reluctance motors in rotor reference frame to estimate rotor losses.However, the article did not provide any analytical voltage harmonic spectrum expressions.In [27]- [28], approximate analytical formulas for specific PWM sideband harmonics were derived.However, this approximate method cannot calculate all harmonic amplitudes and its equivalent circuit also neglect resistance components in the system.
Based on the above summary, we have identified two important aspects in the field of harmonic prediction that lack relevant research and need to be addressed.First, there is a lack of direct harmonic prediction formulas covering the entire frequency domain, Second, existing analytical equivalent circuit models generally used approximations that neglect resistance and rotor eddy current phenomenon including skin effects.To fill the gap, we propose an analytical method that directly calculates current spectral distribution in frequency domain.The proposed method further utilizes a high fidelity harmonic equivalent circuit that incorporates impedance-frequency characteristics in slotless surface-mounted PMSMs including rotor eddy current phenomenon and skin effects.Thus, the proposed method expands the state of the art in the following two aspects: First, the frequency-domain method provides explicit formulas for spectral distribution that covers all harmonics under different circuit topologies.This is a more direct and comprehensive computation of the harmonic distribution.Second, its equivalent circuit model considers eddy current phenomenon including skin effects, which are particularly significant factors in harmonic prediction when LC resonance is present or when voltage harmonic frequencies are high.This greatly improves prediction accuracy compared to the frequency-invariant circuit models used in existing methods The paper is structured as follows.Section II details the frequency-domain harmonics prediction method.Section III presents the equivalent circuit model, including impedancefrequency characteristics of various circuit components.Section IV validates the proposed method through the experimental results and demonstrate the improved prediction accuracy compared to conventional frequency-invariant modeling method.A wide range of operating conditions is tested and carrier sideband harmonic distortion trends with regard to some key design parameters are demonstrated.Finally, conclusions are summarized in Section V.

II. FREQUENCY-DOMAIN PREDICTION METHOD
In this section, the concept of sequence harmonics is first introduced via "double decomposition analysis" process.Second, analytical expressions for positive-and negative-sequence harmonics are developed.Finally, the frequency-domain direct spectral prediction process using sequence harmonics is described.

A. Double Decomposition Analysis
To facilitate the use of equivalent circuit method, the current and voltage waveforms distorted by PWM processes should be transformed into a series of three-phase balanced quantities via "double decomposition analysis", whose processes is symbolically represented as (1) [29].
Without loss of generality, assume there is a group of threephase time series denoted as (t, y a , y b , y c ), where t is the time stamp.First, a Fourier transform (denoted as F) is performed on each phase so that the time series is decomposed into frequencydomain harmonics.The harmonics from all three phases are organized as three phasor arrays labeled with the ordinal h, i.e., (h, Y a , Y b , Y c ).Then, a phase-sequence transform (denoted as T ) is performed for each individual harmonic order h, and a new set of harmonic phasors containing the positive-, negativeand zero-sequence components are obtained.These are called sequence-harmonics, whose ordinals are re-denoted as h+, h−, and h0, respectively.The phase-sequence transform for h-th harmonics can be expressed as where α is the 120°phase-shifting unit vector, defined as The advantage of using sequence-harmonics for analysis is that each harmonic order corresponds to a set of three-phase balanced ac voltages or currents at a particular frequency.Therefore, the three-phase circuit can be simplified and treated as a single-phase equivalent circuit.Once current sequenceharmonics are obtained, abc-phase current harmonic spectra and time-domain waveforms can be obtained by consecutively performing two inverse transformations, i.e., T −1 and F −1 .For example, Y a,h and y a can be obtained by:

B. Closed-Form Analytical Voltage Spectral Distribution
This section presents the double integral harmonic coefficients for PWM voltage waveforms, and derives the analytical expressions for their sequence harmonic spectra.Some useful conclusions about PWM voltage spectral distribution characteristics are obtained.
In [30], voltage harmonic is analytically calculated with the double integral Fourier coefficients A mn , where m ∈ N is the carrier index variable, and n ∈ Z is the baseband index variable.These index variables define the frequency of each harmonic component of the switched phase leg output voltage, as (mf c + nf o ), where f c is the carrier frequency, f o is the target fundamental frequency.In other words, the nth sideband harmonic in the group of harmonics that are located around the mth carrier harmonic (i.e., within the mth carrier sideband group) has a harmonic coefficient of A mn .
For commonly used PWM algorithms, A mn is explicitly expressed with m, n, as well as modulation index M and carrier ratio p, where For example, in the case of naturally sampled sinusoidal PWM (SPWM), A mn takes a simple form when m > 0: where J n is the Bessel function of the first kind.We can define a harmonic amplitude function 2 ) term, f (n) is nonzero only when the parity of n and m is different (n is odd and m is even, or vice versa).
For the space vector PWM (SVPWM) scheme, we also have analytical expressions for A mn , as listed in the Appendix.Fig. 1(b) illustrates the f (n) function curves for SVPWM, along with the integer-sampled bars representing actual harmonics.Like SPWM curves, f (n) also decays when n increases.Generally, SVPWM has a wider spread of spectral distribution.When both n and m are odd or even, the harmonic may still exhibit small amplitude.However, both graph show that the second sideband harmonic of the first carrier group (m = 1, n = 2), the first sideband harmonic of the second carrier group (m = 2, n = 1), and the second sideband harmonic of the third carrier group (m = 3, n = 2), are generally the highest harmonic amid their respective carrier groups.
Next, let's examine the variation of |A mn | with M for these significant sideband harmonics.In Fig. 2, the amplitudes of these three most significant harmonics under SPWM and SVPWM schemes are plotted for M ranging from 0 to 1.15.It is observed that the first carrier sideband harmonic amplitude monotonically increases with increasing M , while the second and third carrier sideband harmonic amplitudes initially rise and then decline.The amplitude of the second carrier sideband harmonic is always greater than that of the third.When M is high, the first carrier sideband harmonic is greater than the second.When M is low, the second carrier sideband harmonic is greater than the first.Fig. 3 compares the analytical prediction with the experimental measurement of the line-line PWM voltage for M = 0.7, p = 20 and V dc = 60 V.It can be seen from the waveforms and spectra that the two match well and the harmonic distribution aligns with the aforementioned characteristics.Some small discrepancies can be attributed to the non-linear factors in the actual hardware, such as switch device resistance, turn on/off delay and overshoot, non-ideal dead time effect, and close-loop controller influence.Some practical factors in the inverter can be modeled with adapted expressions of A mn .For example, Dc bus voltage ripples may be considered using the analytical derivations presented in [31].Analytical harmonic coefficients accounting for dead time effects are provided in [32].
Next, we can derive the expressions for positive, negative, and zero sequence voltages.We can obtain the voltage expressions for the abc phase as follows: Fig. 4. Voltage sequence spectrum.
where V dc is the constant dc bus voltage, θ c and θ o are the initial phase angles of carrier and reference waveform signals, respectively.Angular frequencies ω c and ω o are defined by Combining ( 7) with ( 2), the phasor representation of hth positive, negative, and zero sequence voltage harmonics can be synthesized from all harmonic coefficients as where A * mn is the conjugate of A mn .Thus, F and T transforms in (1) has been done analytically for the voltage waveforms, and the tedious numerical calculation is skipped.Fig. 4 illustrates the theoretically calculated positive and negative sequence voltage harmonic spectra for a two-level inverter controlled by a SVPWM scheme with M = 0.7 and p = 20.Since the neutral point of the motor is generally not grounded and there is no path for zero-sequence current, the zero-sequence voltage harmonics are out of our interest and therefore not shown here.
Based on the derived results, following characteristics of the voltage sequence harmonic spectrum are observed: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

C. Frequency-Domain Current Harmonic Prediction
The process of frequency-domain prediction method for current harmonics is summarized in Fig. 5. First, the voltage sequence harmonic spectrum is analytically determined based on the inverter topology and PWM scheme.Second, based on the distribution of harmonic frequencies, the equivalent circuit parameters for the positive and negative h-th order harmonics are obtained, and the impedance expressions should be calculated.Third, each voltage sequence-harmonic is divided by its corresponding impedance to obtain the current sequence-harmonic, i.e., where I m denotes the machine current, Z t denotes the Thevenin equivalent impedance.Finally, we can convert the 120 sequence currents into abc phase current quantities for all harmonic orders using (11): In general, obtaining the harmonic spectrum of the threephase current is sufficient for machine performance prediction such as electromagnetic loss or torque ripple, and there is no need to continue calculating the time-domain expression.Furthermore, only harmonics within a certain bandwidth are necessarily computed.Usually, harmonics beyond 4f c contribute very little to electromagnetic losses or torque ripples.As a result, the number of harmonics (in the order of tens or hundreds) to be calculated in a frequency-domain method is much smaller than the number of time steps (typically tens of thousands) required in a time-stepping method, which leads to significantly reduced computation time.
Once the current harmonic spectrum is obtained, it is easy to calculate total harmonic distortion (THD).Here, we introduce a new set of parameters indicating the combined amplitude of the sideband harmonics near a certain carrier harmonic.By definition, m-th carrier sideband-harmonic-group distortion (CHD m ) is For completeness, the baseband-harmonic-group distortion (CHD 0 ) is also defined and expressed as The span of these harmonic groups is illustrated in the last graph of Fig. 5. THD can be found through the squared sum of CHD m , i.e.
Compared to THD, CHD m reveals the important frequency distribution of harmonic spectra and can help us better understand the trend of harmonic distortion changes with respect to drive system parameters.

III. SLOTLESS PMSM EQUIVALENT CIRCUIT MODEL
A crucial step in the frequency-domain prediction is the development of sequence harmonic equivalent circuit for the machine drive system, which includes an inverter, a dc bus, a possible filter, and a slotless PMSM, as shown in Fig. 6.This section is dedicated to such a model that includes the impedance-frequency curves of the machine and filter.Analytical expressions will be presented to illustrate their trends.We assume that the dc values of the circuit components are known through FE simulation and focus on the increment of the ac value as the frequency changes.

A. Equivalent Circuit Model for Sequence-Harmonics
The drive-system model includes two types of per-phase lumped-parameter equivalent circuits separately for fundamental harmonic and asynchronous sequence-harmonics [33], as shown in Fig. 9.Each circuit must have three parts modeled, i.e., inverter, filter and machine.
The inverter is represented as a voltage source in each circuit.The fundamental voltage harmonic, and positive-or negativesequence h-th asynchronous voltage harmonics are symbolized as V s,1 and V s,h± , respectively.In Section II-B, the determination of these voltage harmonic components have been presented.
Assume an LC filter is used between the inverter and the motor to reduce the high-frequency harmonic contents.The associated circuit parameters include the inductance (L f ) and resistance (R f ) of the inverter-side inductor, the capacitance of the shunt capacitor (C f ), and the damping resistor (R c ) in series with the capacitor.
Now our focus will be the slotless PMSM equivalent circuits.The establishment of this equivalent circuit is based on the following assumptions: 1) The armature reaction magnetic field in the airgap is a sine wave, and the harmonic leakage component in the magnetizing inductance can be neglected.
2) The saturation effect on the motor inductance can be neglected, and the synchronous inductance does not significantly change with the load current or rotor position.These assumptions are true due to the very large electromagnetic airgap in this machine class.We took a commercial off-the-shelf slotless outer-rotor PMSM shown in Fig. 7(a) as an example and established a finite element model to numerically calculate the magnetic field distribution and inductance.Its geometrical details are tabulated in Table II.Fig. 8(a) shows the   radial and tangential components of the armature reaction field in the air gap region versus rotor position.Fig. 8(b) shows the corresponding spatial harmonic spectrum.Due to the relatively large electromagnetic airgap, the armature reaction field is low, with radial and tangential fundamental components of only 12.2 mT and 8.4 mT, respectively.The spectral analysis results show that the high-order spatial harmonic components are even lower.The THD of the radial and tangential magnetic fields is only 2.2% and 3.5%, respectively.Since inductance is proportional to the square of the magnetic field magnitude, the harmonic inductance accounts for less than 0.1% of the overall magnetizing inductance, which is negligible from an impedance point of view.This means that for a three-phase balanced current excitation of any particular frequency, the armature reaction field in the air gap can be considered as a sine wave with a spatial period of p and a fixed rotating speed proportional to the excitation frequency.derive the synchronous inductance L s .First, incremental phase self-inductance and mutual inductance are calculated in the FE model.Than the synchronous inductance is obtained using the following formula: Where θ is the rotor position, L AA is the self-inductance of phase A winding, M AB is the mutual inductance between phases A and B, and M AC is the mutual inductance between phases A and C, all of which are dc incremental values.It can be seen that due to the large equivalent air gap in the slotless PMSM, the inductance barely changes with the current amplitude or rotor position.Thus, a series circuit comprising the phase resistance and a complex-valued inductance can represent the machine impedance at any particular stator frequency.
First, the fundamental-harmonic circuit, as shown in Fig. 9(a), addresses only the synchronous-frequency (ω 1 = 2πf o ) current.The machine is represented as a serial connection of the synchronous inductance L s , stator resistance R s and back electromotive force (EMF) E 1 .This is the classical synchronous machine steady-state equivalent circuit, which can be used to determine the fundamental voltage phasor V s,1 .Then, M and θ o can be determined according to The equivalent circuits for higher-order sequence-harmonics are developed considering stator asynchronous fields, as shown in Fig. 9(b).The machine impedance is determined by an equivalent circuit consisting of serially-connected stator resistance R s and complex-valued operational inductance L s (jω r ).The concept of operational inductance will be explained in more detail in Section III-C.The slip frequency ω r is needed to determine its value for each sequence-harmonic.
For a particular h, the slip frequency corresponding to the positive-and negative-sequence-harmonic currents are different.The magneto-motive force (MMF) induced by the positivesequence currents rotates in the positive direction at a speed proportional to its electrical frequency hω 1 .The MMF induced by negative-sequence currents rotates at the same speed but in the opposite direction, corresponding to an electrical frequency −hω 1 .The MMF induced by the zero sequence current, if it exists, remains stationary.Unless the winding neutral is connected to the ground, there will be no zero-sequence currents.The electrical frequency of the induced currents in the rotor circuit is the slip frequency ω r , which is determined by the relative speed between the stator MMF traveling wave and the rotor.The rotor is spinning at the synchronous frequency ω 1 at the steady state.Thus, the slip frequency for the positive-sequence h-th harmonic component is (h − 1)ω 1 ; the slip frequency for the negative-sequence h-th harmonic component is (h + 1)ω 1 .In summary, for any sequence harmonic ordinal h±, the slip frequency, or rotor circuit frequency, can be expressed as Then the operational inductance L s in the equivalent circuits is determined by evaluating the L s (jω r ) characteristic curve at ω r,h± .An equivalent-circuit model to compute L s characteristic curve is described in Section III-C.

B. Resistance of Multi-Strand/Multi-Turn Coils
In the windings of magnetic components such as inductors, transformers, and motors, the structure of multi-strand/multiturn coils is widely used.Due to the presence of alternating leakage fields, the current in each strand or turn of the coil exhibits an uneven distribution, resulting in the phenomenon that ac resistance is greater than dc resistance, known as proximity/skin effects.Generally, the higher the frequency, the greater the ac resistance.Assuming that the leakage field is distributed in a parallel manner in the winding area, one can calculate ac resistance with the following equations [34], [35]: where h c is the strand/turn height, σ c is the wire conductivity, b c is the wire width, b is the width between iron core walls, z t is the number of layers in perpendicular to the leakage flux, ξ c is called the reduced conductor height, k R is called the ac resistance factor.When the frequency is high (ξ c > 2), k R will increase proportionally with ξ c , i.e., at a rise rate of ω 1 2 .In general, R f and R s in the equivalent circuit model follows this frequency dependency trend.Fig. 10 plots the predicted and measured resistance value of a 306 μH multi-turn inductor versus excitation frequency, which consolidates the analytical expression.

C. Operational Inductance
As the eddy current density in the rotor conductive parts such as magnets changes with frequency, the motor impedance also changes with frequency.We can study the trend of impedance variation using equivalent circuits.
First, consider the simplest, fundamental-harmonic case where the stator field is synchronous with the rotor.The armature reaction field is stationary with respect to the rotor, so there is no induced eddy current in the magnets.The motor impedance consists only of the magnetizing inductance L m and the stator leakage inductance L sl , as shown in Fig. 11(a).The lumped inductance is the synchronous inductance L s .Now, let's consider a stator excitation of an arbitrary asynchronous frequency ω h ( = ω 1 ).Eddy currents are induced in the rotor magnets as the slip between stator-excited field and the spinning rotor appears.The magnetic field generated by the eddy current generally weakens the original air gap magnetic field according to Lenz's law, thereby reducing the motor inductance.In addition, eddy currents flowing in the rotor magnets with non-zero conductivity will generate Joule losses.The impedance seen from the stator winding ports not only have a reactive component but also a resistive lossy component.These effects grow stronger as the slip frequency between the stator and rotor (which is the frequency of rotor eddy currents, ω r ) increases.In summary, the motor impedance under asynchronous stator current excitation is a complex function that varies with the slip frequency.It is characterized by the so-called operational inductance, L s (jω r ).
The operational inductance can be represented by an equivalent circuit that comprises additionally a rotor branch connected in parallel to the magnetizing inductance [36], [37], [38], as shown in Fig. 11(b).The added branch includes a rotor resistance R r and a rotor leakage inductance L rl .R r is divided by the slip s to account for the frequency difference between stator and rotor currents, where s is expressed as Through simple manipulation of impedance expression from Fig. 11(b), the operational inductance is obtained as Rare earth permanent magnet materials exhibit significant skin effects under high-order harmonics due to their low resistivity and significant geometrical size.When the frequency increases, the eddy current density increases, while the penetration depth decreases, which is similar to the deep bar effects in an induction machine (IM) [39].Just like an IM, R r and L rl of a PMSM also change with rotor frequency ω r .The direct cause of the parametric variation is the circulating current induced by the alternating magnetic field that penetrates the conductor.Fig. 12(a) and (b) respectively illustrate the alternating magnetic flux lines on the rotor bar in IM and the permanent magnet in PMSM, as well as the distribution of the induced circulating currents.As the frequency increases, the circulating currents tend to concentrate on the surfaces of the conductor parallel to the magnetic flux lines with a certain skin depth.Again, we can define the reduced conductor height ξ m to characterize the skin effect: where h m is the magnet circumferential width, b m is the magnet width, σ m is the magnet conductivity, μ m is the magnet Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.permeability, g is the electromagnetic airgap length (equals to the physical airgap plus the magnet depth).
Based on the similarity between the magnetic field and circulating current distribution, we can assume that the frequencydependency functions of R r and L rl of the permanent magnet follows a pattern similar to that of the deep bar conductor in IM [40], i.e., where R r0 and L rl0 are the dc values of these parameters.Both of these dc parameters are determined by the net current distribution (not the circulating current) on the conductor's cross-section and the induced magnetic flux, which will be affected by the external magnetic circuit.Even though ( 27) and ( 28) can be used as good approximation formulas for frequency effects, analytical expression of R r0 and L rl0 for the deep bar conductor cannot be directly applied.Their values in PMSMs can be determined by FE methods.As expected, ( 27) and ( 28) suggest that R r0 increases while L rl0 decreases as ω r increases.When ω r is high (ξ m > 2), R r0 and L rl0 change approximately in the order of ω In Fig. 13, real and imaginary parts of L s (jω r ) of the surfacemounted PMSM are calculated by the analytical equivalent circuit model.The real part corresponds to the inductive component, while the imaginary part corresponds to the resistive component.The real part is always positive and monotonically decreases with frequency.The imaginary part is always negative owing to real power consumption, and its absolute value first increases and then decreases with the frequency.The analytical results are validated by motor impedance measurements with an HP 4284 A Precision LCR meter and good agreement is found.Some small discrepancies may be attributed to the distributed parasitic circuit parameters.

D. Positive and Negative Sequence Equivalent Impedance
Once we obtain the impedance-frequency characteristic curves of each component, we can then calculate the equivalent impedance required by (10).Based on the circuit model in Fig. 9(b), the Thevenin equivalent impedance Z t for the harmonic ordinal h± is derived in (29) to calculate the machine harmonic current excited by a particular voltage harmonic: where the impedance are given by and the stator's harmonic frequency ω h = hω 1 .
If we assume that all the resistance can be neglected and L s ≈ L s at the resonant frequency f res , then f res can be obtained as approximately Generally, filter parameters should be selected so that f res resides between f o and f c .If an L filter is used instead of the LC filter, the impedance expression is simply: If no filter is used, then Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Table IV provides an overview of the motor and filter system parameters used in this study.The filter inductance is set at 305 μH, whereas the filter capacitance is 60 μF, and the damping resistance is 0.2 Ω.The resonant frequency f res is 5.6 kHz.

TABLE IV CIRCUIT PARAMETERS OF THE EXPERIMENTAL SYSTEM
A review of previous literature has shown that the existing analytical models often assume frequency-invariant equivalent circuit parameters.Now, let's compare the proposed model and a frequency-invariant model and examine the change of the equivalent impedance when frequency factors are taken into account.
Firstly, consider the simplest scenario, where there is no filter and only the motor is present, represented by (33).The impedance corresponding to constant inductance, denoted as Z 0 (= ω h L s ), the positive sequence impedance Z 1 (= |Z t,h+ |), the negative sequence impedance Z 2 (= |Z t,h− |), and their percentage differences are illustrated in Fig. 14(a).It can be observed that at lower frequencies, the positive and negative sequence impedances show no significant difference compared to the impedances calculated using the frequency-invariant model.However, as the frequency increases, differences gradually become pronounced, with the positive and negative sequence impedance being lower.Moreover, the higher the frequency, the higher the percentage decrease in impedance.
Next, consider the LC circuit scenario represented by (29), and the curves are ploted in Fig. 14(b).This situation is more complex than the previous one.The impedance characteristic curve is divided by f res into two segments.Below f res , the impedance is capacitive.Above f res , the impedance becomes inductive.Around f res , due to mutual cancellation between L and C, the impedance exhibits a significant dip, constituting primarily a resistive component.
As seen from the earlier discussion, the eddy current phenomenon including skin effects manifests in the lumped parameters primarily in two ways: firstly, by increasing the resistive impedance, and secondly, by causing a decrease in the inductive impedance.At low frequencies, Z 0 , Z 1 and Z 2 show no significant difference, as the capacitance does not vary with frequency.When the frequency approaches the resonance point (5.6 kHz), the frequency effect leads to an increase in the equivalent impedance Z 1 and Z 2 , which can reach up to 10%.As the frequency surpasses the resonance point to a certain extent, it induces a decrease in Z 1 and Z 2 compared to Z 0 , and as the frequency increases, the impedance reduction becomes more pronounced.By the time we reach 40 kHz, impedance can be reduced by up to 20% compared to Z 0 .From this, it becomes evident that neglecting the impedance-frequency characteristic curve can introduce significant errors in harmonic prediction.
Once the expression for Z t is determined, we possess a set of analytical equations to determine the current harmonics.For instance, for a slotless PMSM system with an LC filter and a 2-level inverter driven by SVPWM, ( 9), (10), (11), ( 29), (30), and (36) constitute a set of close-form prediction equations.They can be implemented using the flowchart in Fig. 15.A MATLAB execution code has been uploaded to the open source code repository website GitHub [41].In the following, we will use this process to predict current harmonic spectra for an experimental setup under various settings.The effectiveness of the proposed method is validated by comparing it against actual test results.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

IV. MODEL VALIDATION WITH EXPERIMENTAL DATA ANALYSIS
We will scrutinize the current harmonic spectra obtained from the model prediction and experimental results under various scenarios.The prediction errors from both the proposed method and the existing method that do not consider frequency-impedance characteristics are presented and compared.

A. Experimental Setup and Procedure
The experimental system hardware is shown in Fig. 16.It comprises the same slotless PMSM (Brand name: ThinGap) as presented in the previous analysis.The inverter is a Texas Intrument's High Voltage Digital Motor Control and Power Factor Correction kit (Part Number: TMDSHVMTRPFCKIT).We ran the PM motor in speed control mode and regulated the current by synchronous PI control, using a symmetrical SVPWM algorithm.Furthermore, the dynamometer operated in torque mode and provided a load torque that kept the motor current at 4 A rms .The instantaneous machine current signals were measured by Tektronix TCP0030 A current probes with 1.5% accuracy, which were connected to a Tektronix Model MSO-4034 Oscilloscope.The oscilloscope recorded the data and stored them on a PC.We then performed an Fast Fourier Transform (FFT) on the saved data to obtain the harmonic spectrum.
We investigated how changes in carrier frequency f c , dc bus voltage V dc , filter capacitance C f , and modulation index M affected the current harmonic spectrum, as predicted by our model, and compared the findings with the experimental results.Fig. 17 depicts the predicted and experimental spectra and waveforms of six representative testing cases.Fig. 18 illustrates the trends of CHD 0 -CHD 2 with respect to the aforementioned operating parameters.We are going to discuss our observations in detail.

B. Carrier Frequency
An LC filter with 60 μF capacitance is used for studying the influence of f c .A constant operating speed for the motor is maintained, which implies that f o and back EMF amplitude remain unchanged.Additionally, we keep V dc constant, which ensures that the modulation index M remains almost unchanged.At this point, we examine the spectral and waveform characteristics under two different f c , as shown in Fig. 17(a)-(d).Here, we maintain that f o = 400 Hz, V dc = 60 V, and M = 0.8 for both cases, while their f c are 5.6 kHz and 7.2 kHz, respectively, giving rise to p of 14 and 18 for each case.
It is observed that the higher the carrier frequency, the lower the amplitudes of the current harmonic components.This is attributed to the fact that high-amplitude voltage harmonics are always distributed near integer multiples of the carrier frequency.Therefore, the higher the carrier frequency, the higher the frequencies of high-amplitude voltage harmonics, and the larger the impedance corresponding to those harmonics.This results in lower amplitude of current harmonics.As such, we observe a clear negative correlation between the carrier ratio and harmonic amplitude, as depicted in Fig. 18(a), which plots the variations of CHD 0 to CHD 2 in the range of p, or the normalized f c with f o , from 12 to 20.
Our experimental results demonstrate that all three CHDs exhibit decreasing trend with respect to a rising f c .We note that this observation is in excellent agreement with our predicted values.Given the significant influence of f c on the current harmonic amplitudes, we maintain f c constant when analyzing the effects of other parameters in the following tests.When M = 0.55, the amplitude of the second carrier sideband-group harmonics is higher than that of the first carrier sideband-group harmonics.However, when M = 1.0, the opposite is true.This is because M has a different effect on the harmonics belonging to these two individual sideband-groups.As shown in Fig. 18(b), CHD 1 monotonically increases with M , while CHD 2 first increases and then decreases.In general, the higher the value of M , the higher the harmonic content of the first sideband.Because the impedance corresponding to lower carrier order is smaller, THD is also higher for higher M .

C. Modulation Index
It is worth noting that the switch dead time results in highamplitude components in the baseband harmonic group.The voltage distortion component of dead time is generally proportional to f c and V dc , while the corresponding impedance is proportional to f o .Therefore, as M decreases and f o becomes lower, CHD 0 becomes larger.The distorted waveforms caused by dead time are also more clearly visible in Fig. 17(f).
Furthermore, the good agreement between prediction and test results shows that the proposed analytical method is valid for motor drive systems with both LC and L filters.This shows the method's versatility even with high-order circuits.

D. Dc Bus Voltage
Keeping f c and M constant, influence of V dc is studied.By reducing V dc of the Fig. 17    almost the same.The major difference is that all harmonics in the 33-V case have about half amplitudes of those in the 60-V case.This is because the PWM voltage harmonic amplitudes are directly proportional to V dc .Under the same f c , M , and system impedance, current harmonic amplitudes should follow the same linear relationship.Fig. 18(c) shows the trend of CHD with respect to V dc .CHD 1 and CHD 2 increase almost linearly with V dc , consistent with theoretical predictions.However, CHD 0 shows the opposite trend due to lower dead-time distortion at higher f o .

E. Filter Capacitance
The influence of C f variation on the current spectrum is investigated.First, we compare the cases with and without C f .In the L-filtered case of Fig. 17  the system to an LC-filtered case and the results are shown in Fig. 17(k) and (l).The system now has a higher f res of 13.8 kHz, while f c is 12 kHz.The first carrier-sideband-group harmonics will excite the resonant mode when f c < f res .The amplitude of current harmonics is even higher than that of the case when no capacitor is added.The second carrier-group harmonics are farther from f res , so the high-frequency harmonic components in the current waveforms are bypassed through the shunt C f , resulting in a decrease in the second carrier-sideband-group harmonics.
Fig. 18(d) shows the trend of CHD m for C f ranging from 0 to 60 μF.When we increase C f further beyond 10 μF, f res decreases and moves away from f c .Therefore, CHD 1 curve decreases with an increase in C f .Meanwhile, the larger the C f , the stronger its filtering effect, resulting in a lower CHD 2 .Thus, as long as f res is substantially lower than f c , adding filter capacitance is an effective means to reduce the current harmonics.As p remains constant at 30, dead-time distortion does not change much, CHD 0 curve is almost flat.

F. Sideband Carrier Harmonic Prediction Error
Next, we will compare the prediction errors between the improved model we proposed and the original frequency-invariant model.We define individual harmonic prediction error in relative percentage as P.E.(h) can be used to demonstrate the prediction error distribution over the frequency spectrum to identify any systematic prediction errors and analyze their patterns.The root mean squared (rms) weighted prediction error in relative percentage for all significant harmonics is defined as where h shall take all corresponding ordinals to significant harmonics, i.e., h = p ± 1, p ± 2, p ± 4, 2p ± 1, 3p ± 2. P.E.rms is equivalently a weighted average gauge that provides a good measure of the overall harmonic prediction accuracy.Fig. 19 illustrates the P.E.(h) distribution of the two models under three different cases.From cases (a) and (b), the original frequency-invariant model clearly exhibits systematic prediction errors, which overestimates harmonics around the resonant frequencies and underestimates them at high frequencies, aligning with our earlier impedance observations.In contrast, the proposed model effectively reduces such prediction errors.In the high-order harmonic range, the proposed model might also exhibit relatively large errors, which can be attributed to the inherently low amplitudes in this frequency range.Lowamplitude harmonics are more susceptible to external dynamic disturbances and more difficult to measure precisely.In case (c), where no capacitor is present and the equivalent circuit only consists of a large filter inductor in series with the motor, the frequency-dependent effects on the overall impedance have been diluted.Here, harmonic prediction errors mainly arise from an imperfect modeling of the PWM voltage spectrum.
Table V presents the rms prediction error in relative percentage across various scenarios.In cases corresponding to LC filtered configuration, the original frequency-invariant model, due to the lack of frequency effect consideration, may yield prediction errors up to 36.7%.The proposed model notably enhances prediction accuracy by modeling the frequency effect.The reduction rate of relative percentage error is between 25.9% to 75.7%.In cases corresponding to L-filtered configuration, frequency effect is limited and our model provides modest improvement by 7-10%.Overall, the prediction error is generally maintained below 10% with the proposed model, indicating a satisfactory prediction accuracy for practical applications.These comparative results underscore the effectiveness of the proposed method and its significant improvement compared to the conventional model.

V. CONCLUSION
This article presents an analytical current harmonic prediction method for slotless surface-mounted PMSMs.The spectral prediction is performed directly in the frequency domain using double-integral Fourier coefficients and a sequence-harmonic equivalent circuit model.The model incorporates impedancefrequency characteristics of each circuit component, thus forming an accurate representation of the Thevenin equivalent impedance.The variation trends of the sideband harmonics under the influence of changing carrier frequency, modulation index, dc bus voltage, and filter capacitance are predicted and explained by the frequency-domain model.The predictions are aligned well with the experimental results.A comparison has been made on the prediction errors from the proposed model and the conventional frequency-invariant circuit model.The proposed model comprehensively reduces prediction errors, and in certain scenarios, prediction accuracy can be improved by 75%.

APPENDIX
The Fourier coefficients A mn for voltage harmonics with regularly-sampled, symmetrical carrier SVPWM is given by [30] as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where J n is the Bessel functions of the first kind, and q = m + n p .For asymmetrical carrier SVPWM and various other PWM schemes, analytical expressions for A mn can also be found in [30].

Manuscript received 8
April 2023; revised 13 August 2023; accepted 14 November 2023.Date of publication 28 November 2023; date of current version 23 May 2024.This work was supported in part by the National Science Foundation Engineering Research Center for Power Optimization of Electro-Thermal systems (POETS) with cooperative agreements under Grant EEC-1449548, and in part by the Grainger Center for Electric Machinery and Electromechanics.Paper no. TEC-00393-2023.(Corresponding author: Xiaolong Zhang.)
symmetric about the y-axis.This allows us to focus on the variation of f (n) for n > 0. Fig. 1(a) plots the curves of f (n) for M = 0.7, with different colors representing m = 1, 2, and 3, separately.The figure also shows the curve sampled at positive integers n, which are the harmonic amplitude values.Two features are identified.First, the envelope of f (n) rapidly decays with increasing n, owing to the Bessel function term.Second, due to the modulation effect of the sin([m + n] π

Fig. 12 .
Fig. 12. Alternating flux lines and induced circulating currents in rotor conductive parts.

Fig. 14 .
Fig. 14.Machine and system impedance under positive and negative sequence harmonics (Z 0 : Impedance from frequency-invariant model, Z 1 : Positive sequence impedance from proposed model, Z 2 : Negative sequence impedance from proposed model).

Fig. 17
Fig. 17(e) -(h) illustrate the spectra and waveforms of two different M .Here, we switched to using an L filter, raised the f c to 12 kHz and kept V dc at 50 V.In the Fig. 17(e) and (f) case, f o is 200 Hz and M = 0.55.In the Fig. 17(g) and (h) case, f o is 400 Hz and M = 1.0.When M = 0.55, the amplitude of the second carrier sideband-group harmonics is higher than that of the first carrier sideband-group harmonics.However, when M = 1.0, the opposite is true.This is because M has a different effect on the harmonics belonging to these two individual sideband-groups.As shown in Fig.18(b), CHD 1 monotonically increases with M , while CHD 2 first increases and then decreases.In general, the higher the value of M , the higher the harmonic content of the first sideband.Because the impedance corresponding to lower carrier order is smaller, THD is also higher for higher M .It is worth noting that the switch dead time results in highamplitude components in the baseband harmonic group.The voltage distortion component of dead time is generally proportional to f c and V dc , while the corresponding impedance is proportional to f o .Therefore, as M decreases and f o becomes lower, CHD 0 becomes larger.The distorted waveforms caused by dead time are also more clearly visible in Fig.17(f).Furthermore, the good agreement between prediction and test results shows that the proposed analytical method is valid for motor drive systems with both LC and L filters.This shows the method's versatility even with high-order circuits.
Fig. 17(e) -(h) illustrate the spectra and waveforms of two different M .Here, we switched to using an L filter, raised the f c to 12 kHz and kept V dc at 50 V.In the Fig. 17(e) and (f) case, f o is 200 Hz and M = 0.55.In the Fig. 17(g) and (h) case, f o is 400 Hz and M = 1.0.When M = 0.55, the amplitude of the second carrier sideband-group harmonics is higher than that of the first carrier sideband-group harmonics.However, when M = 1.0, the opposite is true.This is because M has a different effect on the harmonics belonging to these two individual sideband-groups.As shown in Fig.18(b), CHD 1 monotonically increases with M , while CHD 2 first increases and then decreases.In general, the higher the value of M , the higher the harmonic content of the first sideband.Because the impedance corresponding to lower carrier order is smaller, THD is also higher for higher M .It is worth noting that the switch dead time results in highamplitude components in the baseband harmonic group.The voltage distortion component of dead time is generally proportional to f c and V dc , while the corresponding impedance is proportional to f o .Therefore, as M decreases and f o becomes lower, CHD 0 becomes larger.The distorted waveforms caused by dead time are also more clearly visible in Fig.17(f).Furthermore, the good agreement between prediction and test results shows that the proposed analytical method is valid for motor drive systems with both LC and L filters.This shows the method's versatility even with high-order circuits.
Fig. 17(e) -(h) illustrate the spectra and waveforms of two different M .Here, we switched to using an L filter, raised the f c to 12 kHz and kept V dc at 50 V.In the Fig. 17(e) and (f) case, f o is 200 Hz and M = 0.55.In the Fig. 17(g) and (h) case, f o is 400 Hz and M = 1.0.When M = 0.55, the amplitude of the second carrier sideband-group harmonics is higher than that of the first carrier sideband-group harmonics.However, when M = 1.0, the opposite is true.This is because M has a different effect on the harmonics belonging to these two individual sideband-groups.As shown in Fig.18(b), CHD 1 monotonically increases with M , while CHD 2 first increases and then decreases.In general, the higher the value of M , the higher the harmonic content of the first sideband.Because the impedance corresponding to lower carrier order is smaller, THD is also higher for higher M .It is worth noting that the switch dead time results in highamplitude components in the baseband harmonic group.The voltage distortion component of dead time is generally proportional to f c and V dc , while the corresponding impedance is proportional to f o .Therefore, as M decreases and f o becomes lower, CHD 0 becomes larger.The distorted waveforms caused by dead time are also more clearly visible in Fig.17(f).Furthermore, the good agreement between prediction and test results shows that the proposed analytical method is valid for motor drive systems with both LC and L filters.This shows the method's versatility even with high-order circuits.
Fig. 17(e) -(h) illustrate the spectra and waveforms of two different M .Here, we switched to using an L filter, raised the f c to 12 kHz and kept V dc at 50 V.In the Fig. 17(e) and (f) case, f o is 200 Hz and M = 0.55.In the Fig. 17(g) and (h) case, f o is 400 Hz and M = 1.0.When M = 0.55, the amplitude of the second carrier sideband-group harmonics is higher than that of the first carrier sideband-group harmonics.However, when M = 1.0, the opposite is true.This is because M has a different effect on the harmonics belonging to these two individual sideband-groups.As shown in Fig.18(b), CHD 1 monotonically increases with M , while CHD 2 first increases and then decreases.In general, the higher the value of M , the higher the harmonic content of the first sideband.Because the impedance corresponding to lower carrier order is smaller, THD is also higher for higher M .It is worth noting that the switch dead time results in highamplitude components in the baseband harmonic group.The voltage distortion component of dead time is generally proportional to f c and V dc , while the corresponding impedance is proportional to f o .Therefore, as M decreases and f o becomes lower, CHD 0 becomes larger.The distorted waveforms caused by dead time are also more clearly visible in Fig.17(f).Furthermore, the good agreement between prediction and test results shows that the proposed analytical method is valid for motor drive systems with both LC and L filters.This shows the method's versatility even with high-order circuits.
(a) case from 60 V to 33 V and f o from 400 Hz to 200 Hz, we obtain the spectrum and waveforms in Fig.17 (i) and (j).f c and M are kept at 7.2 kHz and 0.8, respectively.The shapes of the spectra in these two cases are

Fig. 17
Fig. 17.Predicted and Measured Current Waveforms and Spectra.
Fig. 17.Predicted and Measured Current Waveforms and Spectra.
(g), adding a 10 μF capacitance changes Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE II SLOTLESS
PMSM MACHINE GEOMETRICAL PARAMETERS

TABLE III VARIATION
RANGE OF INCREMENTAL SYNCHRONOUS INDUCTANCE OVER DIFFERENT ROTOR POSITIONS AND CURRENT AMPLITUDES

TABLE V RMS
PREDICTION ERROR COMPARISON FOR SIGNIFICANT CARRIER HARMONICS Fig. 19.Prediction Error Comparison for Significant Carrier Harmonics.