Dynamic Phasor Finite Element Modeling of Grid-Connected DFIG Considering Winding Space Harmonics

This paper presents a new technique for quantifying the winding MMF space harmonics of wind-driven Doubly-Fed Induction Generator (DFIG) for power quality assessments. Instead of following the commercial electromagnetic transient programs in modeling electric machines using lumped-parameters models, the paper uses a combined Dynamic Phase modeling and Finite Element Method (FEM) to model the generator accurately while reducing the computational burden and simulation time associated with the conventional time-stepped finite element analysis. The technique has the ability to efficiently model core nonlinearity and stator/rotor MMF space harmonics. The results of a C++ code written by the authors of the new dynamic phasors FEM technique are compared to the conventional time-stepping FEM to show the capability of the method to produce comparable results in shorter simulation time.


I. INTRODUCTION
W ITH the increased attention towards power quality in electric generators and the power systems receiving their generated power, the need for more research about the sources of harmonic pollution is on the rise. Harmonics are sinusoidal voltages or currents with frequencies other than the nominal 50 or 60Hz superimposed on the original waveform. This harmonic pollution leads to harmful consequences such as: decreased life-time of equipment, increased power losses in machines and capacitor banks, malfunction of control systems and interference with data lines [1]- [4]. Today's power system is heavily polluted by voltage and current harmonics due to nonlinear characteristics of a broad range of equipment and systems especially electric machines [1].
Unlike the widely used electric generator's modeling technique of ignoring the emf and current harmonics in most of the commercial power system simulators, the actual construction of electric machines imposes time harmonics in the electric current waveform. The discrete distribution of winding around the circumference of the machine's stator and rotor, nonuniform stator and rotor inner circumference facing the airgap due to slotting effect, and rotor eccentricity caused by manufacturing errors, shaft bending or bearing faults, cause space harmonics in the airgap flux density waveform. These space harmonics are translated into time harmonics and interharmonics in the machine's back-emf and current. This paper aims to present an accurate and efficient way of modeling time harmonics and interharmonics in the stator current waveform of a grid-connected Doubly Fed Induction Generator caused by winding space harmonics in order to be used as an accurate way of quantifying its impact on the generator's power quality when connected to an external power system. In order to model the impact of DFIG's space-harmonic content on its connected power system, the machine has to be co-simulated with the system in a real-time environment. The co-simulation of electric machines and power system has always been the target of extensive research to tackle its famous problems of long simulation times and high computational burden. Because of the slow electromechanical transients in electric motors and generators, and the faster electromagnetic transients in system's inductive and capacitive elements, co-simulation studies have to be performed at very small time-steps to account for the higher stator current harmonic frequencies associated with winding space-harmonic.
For the traditional Electro-Magnetic Transient Programs (EMTPs) dealing with power system simulation, electric machines are usually modeled using lumped-parameters models in natural abc or dqo frames [5]. In spite of the simplicity and small time consumption associated with lumped-parameter models, they are originally derived assuming perfect winding distribution, ignoring slotting effect and ignoring core nonlinearity [6]. Therefore, several research attempts tried modifying lumped-parameters models to accommodate winding space harmonics. In [7], the authors used a coupled-circuit approach to model space harmonics in Induction Machines. Every set of harmonics with possible interaction are grouped together and transformed to a specific axes plane. However, it cannot be transformed to dq space unless every single harmonic is dealt with individually which reduces the model's benefit. In [8], the authors discussed a space-vector model for electrical machines. Space harmonics are included by adding corresponding harmonic impedances to the dynamic lumped-parameters model. However, the authors did not state the calculation method for the inductances representing space harmonics and ignored several nonlinearities for the sake of simplicity. In [9], the authors modified the traditional equivalent circuit of squirrel cage Induction Machine to accommodate space and time harmonics by adding corresponding impedances and branches allowing for interaction between different harmonic components. The model lacks the ability to model many other magnetic aspects like slotting effect, saturation and eddy currents. In [10], the authors presented a dynamic model representing space harmonics in the airgap of a five-phase Induction Machine. They used simple analytical model to evaluate self and mutual inductances associated with each harmonic. For saturation and slotting effects, the authors used analytical factors modifying the inductance calculation. However, more accurate representation is needed to represent nonlinearities. In [11], the authors developed an analytical model for a radial-flux external-rotor permanentmagnet synchronous machine to be used as a part of an energy storage flywheel. They modeled the machine analytically using the vector form of Poisson's equation. Winding space harmonics are represented by a space-varying function used for calculating the current density vector. However, the general analytic model lacks the ability to model local phenomenon especially local and global core saturation. In [12], the authors presented a lumped-parameters model to represent low order space harmonics in six-phase Induction Machine based on Vector Space Decomposition. The method decomposes the machine's variables into a torque producing subspace and other multiple orthogonal subspaces representing various harmonics. However, the method uses high number of equations especially for considering high number of harmonics which increases the complexity and reduces the method flexibility.
In order to better represent the physical winding distribution and the resultant space harmonics, Winding Function Theory has been introduced in [13] and modified in [14] to accommodate varying airgap length. In [15], the authors used a traditional abc lumped-parameters model to dynamically model the machine while using Winding Function theory to represent space harmonics in healthy and faulty conditions. However, it lacks the accuracy of modeling saturation and eddy currents. In [16], the author used Winding Function theory to represent the space harmonics in the winding inductance representation to be fed to a traditional lumpedparameters model. The ideal stepped pattern of winding distribution utilized in the study reduces the accuracy along with neglecting saturation, slotting and eddy currents. In [17], the authors presented an optimized selection of currents to reduce the torque ripples in Synchronous Reluctance machine using vector control. They used the traditional lumpedparameters model with a winding-Function determined inductances. Most of the nonlinearities are not included in the analysis. In [18], the authors used a Fourier expansion to represent space harmonics represented by Winding Function approach. In addition, they used coordinate transformation to form the required inductance matrices of the traditional lumped-parameters model. Although the simulated results matches the experimental results, no mention of representing core saturation, slotting or eddy current calculation. In [19], this paper presents a transient model for the simulation of cage induction motors accounting for the actual topology of the machine based on the space vector theory. The model adopts an analytical approach to evaluate the Fourier components of the airgap Magneto-Motive Force. The authors tried to reduce the computational burden by sacrificing the accuracy in terms of assuming perfect rectangular MMF wave (as the simplest winding function) and ignoring saturation.
For better representation of the machine's physical topology and local magnetic phenomena, Magnetic Equivalent Circuit method has its popularity. In [20], the authors used Magnetic Equivalent Circuit method to represent the actual topology of Induction Machine. Winding space harmonics, saturation and other nonlinearities have been represented in a fairly accurate way compared to lumped-parameters model and less complicated compared to FEM. However, in order to get the same accuracy level of FEM, a very complicated MEC has to be used sacrificing the intended reduction in simulation time. In [21], the authors used Fourier analysis of the MMF waveform to model the magnetic behavior and space harmonics effects. In addition, they used Magnetic Equivalent Circuit to model teeth topology and magnetic saturation. However, a considerable mismatch has been noticed against FEM in terms of local saturation and tooth ripples due to considering a few number of magnetic paths in their Magnetic Circuit.
Finite Element Method is the choice whenever the accuracy of machine's modeling is the target. By discretizing the solution space into small triangular mesh elements, any complex geometry can be modeled in a fairely accurate way [22]. Therefore, modeling the actual winding distribution, slotting effect and core saturation can be performed efficiently by FEM in order to represent the winding space harmonics' impact on stator current's harmonic content. In [23], the authors used 2D FEM to construct flux maps linking winding flux 2 VOLUME 4, 2016 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. linkage with dq axes currents and rotor position. Utilizing FEM allowed the authors to use their flux maps to represent various nonlinearities like saturation, cross saturation and space harmonics. In [24], the authors used 2D FEM to model the various harmonics impact on power losses. They adjusted the traditional T-type equivalent circuit of Induction Motor to accommodate the harmonic power losses.
In order to simulate the generator and its connected external system dynamically, the whole finite element model of the machine has to be solved at each time step which called time-stepped finite element method (TS-FEM). In order to account for the fast system dynamics, the simulation time-step has to be significantly small which requires impractically long simulation time and high computational power. Some research efforts have been dedicated to shorten the simulation time associated with TS-FEM to facilitate co-simulation studies of machines and power system. In [25], the authors solved the machine's steady-state per-phase equivalent circuit to provide the FEM solver with an initial condition to bypass the starting period while in [26], [27] the authors used frequency-domain FEM to do the same job. In [23] and [28], the authors used FEM to prepare offline flux-current data to be used for online simulation. This paper uses a combination of Dynamic Phase modeling (DPM) technique and FEM to accelerate the precise modeling of grid-connected DFIG. In [29], the authors proposed a Model Reduction Technique to reduce the computational burden of circuit connected FEM. The paper investigates selecting the proper set of winding harmonics to compromise the accuracy and speed of simulation. In [30], the authors proposed a method of representing space harmonics in the fast complexdomain FEM by resolving the magnetic vector potential into its Fourier components. In order to represent rotor movement, separated stator and rotor domains have been adopted.
Dynamic Phase modeling techniques approximates the time-domain signal over a sliding window by a summation of complex-valued Fourier coefficients called dynamic phasors. Compared to time-domain, DPM is a much faster modeling technique. This is due to the fact that the oscillating waveforms of AC circuits become constant or slowly-varying if represented by DPM allowing to use large step sizes in numerical simulations. DPM has been used in [31] to model the dynamic behaviour of DFIG using lumped-parameters model. This paper uses a reformulation of the state variables of finite element analysis into dynamic phasors to shorten the simulation time required for co-simulation studies of a DFIG and an external power system. The research work focuses on representing the MMF space harmonics' reflection in the harmonic content of stator current using the new DPM-FEM technique. DPM-FEM technique has been entirely coded in C++ programming language along with the alternative timedomain FEM in order to get a clear comparison of their results and execution time. Section-2 reviews the basic formulation of the conventional time-domain FEM considering core nonlinearity, the basic principles of Dynamic Phasors and the formulation of the new DPM-FEM method. Section-3 compares the simulation results of the time-domain FEM written code with the results of a commercial FEM software to validate it accuracy. Section-4 presents the topology of the system used as a simulation case-study. Section-5 presents the proposed procedure for calculating the reflection of winding space harmonics into the stator current's time harmonics. Section-6 presents the comparison of simulation results of time-domain FEM and DPM-FEM techniques. Section-7 presents the research conclusions.

II. THEORETICAL FOUNDATION OF FEM
The basic formulation of the widely-used time-stepped FEM in time-domain is reviewed in this section. Then, the new DPM-FEM method's detailed equations are derived based on the principles of complex-valued frequency-domain FEM. Representation of core's magnetic saturation in both timedomain and DPM environments is discussed to guarantee an accurate representation of the magnetic behavior of the machine.

A. TIME-DOMAIN FEM
The 2D space and time variation of the magnetic vector potential (A) may be found by solving the well-known matrix equation of the form [22]: where, [F e ] is the normalized applied current density to the local element called "the forcing function" and [A e ] is a three-element column representing the nodal magnetic vector potentials of the local element.
[S e ] and [T e ] are geometrical matrices related to the element coordinates and dimensions [22]. For the required time-stepping solution, the FEM magnetic equation is discretized to a time-step ∆t using Crank-Nicholson Method. In order to represent core nonlinearity, the magnetic reluctivity (ν) is considered as a function of the magnetic vector potential changing with it in a timely basis. Newton-Raphson method is used to solve the system equation taking the form: where n is the Newton-Raphson iteration index, and [G] is calculated as: where B is the flux density in Tesla. The value ∂ν ∂B 2 represents the rate of core-reluctivity change with the change of the flux density. That is where the actual magnetization curve of the core material is introduced to represent core VOLUME 4, 2016 nonlinearity. The discrete B-H data of the core material is recalculated to get the corresponding ν − B 2 values, interpolated using Cubic Spline to get a representing equation that can be used to find the value of ∂ν ∂B 2 corresponding to the current value of core flux density. The elemental matrices are assembled to represent the sharing of different elements connected at a certain node, this forms the global system equation in the form of Ax = b which is solved by Conjugate Gradient method to get the nodal magnetic vector potentials.

B. FUNDAMENTALS OF DYNAMIC PHASE MODELING
Dynamic Phase modeling concept means approximating the time-domain signal x(τ ) in a time window τ ∈ (t − T, t) with a Fourier expansion as [31]: where T is the fundamental period in seconds, ω is the fundamental frequency in rad./sec. (ω = 2π T ), X k (t) is the k th Fourier coefficient and K is the chosen set of dynamic phasors to approximate the actual wave. The value X k (t) is complex-valued and time-varying coefficient, called Dynamic Phasor and written as: ⟨x⟩ k (t).
From the first sight, the formulation of DPM in complexnumbers seems more complicated compared to regular timedomain modeling as it is required to separate real and imaginary parts of the equations which means using twice the number of equations for each harmonic. But even with much higher number of system equations, DPM leads to a significant reduction in simulation time because Dynamic Phasors are much slower than the corresponding time-domain signals, allowing for much larger simulation time-steps. However, choosing the proper set of dynamic phasors is crucial for the success of modeling the targeted system using this technique. Too low a number of considered harmonics leads to poor system representation while too many harmonics leads to a complex and a computationally-expensive numerical model. Therefore, expert knowledge of the studied system is required in order to determine the best harmonic combination to model the system with a good compromise between accuracy and complexity.

C. DYNAMIC PHASOR FEM FORMULATION
Time-domain FEM has the advantage of capturing most of the machine's physical phenomenons of concern and representing all the required harmonics if small time-step is used. Smaller time-steps along with the need for changing the representation of machine's topology for each incremental position makes time-domain FEM (especially with nonlinear core representation) computationally-expensive. In this subsection, the required formulation for merging Dynamic Phase modeling with FEM is presented.
If the objective is modeling just the steady state performance, frequency-domain FEM is used, and (1) may be written as: where ⃗ A e and ⃗ F e are complex-valued vectors of A e and F e respectively, ω e is the operating frequency in rad./sec, and j is the complex operator.
Single-frequency steady state FEM is relatively fast but cannot capture the time-change of system inputs, starting behavior or transients. In order to find a mid-point between the two extremes of time-domain and single-frequency FEM, a merge between Dynamic Phase modeling technique and FEM is used in this paper to model winding space harmonics in stator current waveform. Using DPM, a much larger timestep can be used for time-stepping solution while solving the machine's electromagnetic equation and the equations of the external systems in complex-valued frequency based environment at each time-step. The elemental Poisson's equation in DPM domain takes the form: where ⟨ ⃗ A e ⟩ k and ⟨ ⃗ F e ⟩ k are the complex-valued k th dynamic phasors of magnetic vector potentials and the applied current density respectively and ω 1 is the fundamental frequency in rad./sec. This equation clearly shows that DPM-FEM is in the mid-way between the real-valued time-stepping FEM and the complex-valued frequency-domain FEM. Just like the aforementioned time-domain method, this elemental matrices are used to form the corresponding global matrices by summing up the nodal contribution from different elements. The previous equation is solved for each harmonic at each time-step making DPM-FEM more mathematically complicated compared to regular time-domain FEM. However, the nature of DPM technique enables using much larger simulation time-steps which reduces the overall mathematical burden and shortens simulation time. In addition, solving in time-frequency domain requires no change in rotor position throughout the simulation process which saves the mathematical power and time required for re-meshing the airgap for each new rotor position. In order to model core nonlinearity in DPM-FEM, the Newton-Raphson approach is again used here, transforming the time-discretized DPM-FEM magnetic equation into:

VOLUME 4, 2016
This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. where [G f ] is calculated as: where ⟨ ⃗ A e * ⟩ is the complex-conjugate of the magnetic vector potential dynamic phasor and | ⃗ B| is the absolute value of the element's magnetic flux density. Unlike time-domain FEM, the B-H curve of the core material cannot be traced point by point in frequency-domain solution. Instead, an effective AC value of the magnetic reluctivity (ν f req ) has to be chosen as an equivalent for the whole cycle of time-variation. Different approaches can be used to identify this frequency-domain value of reluctivity, the adopted method in this study is using a value of reluctivity that gives the same average energy of the actual magnetization curve [20].

III. VALIDATING THE TIME-DOMAIN FEM WRITTEN CODE AGAINST A COMMERCIAL SOFTWARE
The time-domain FEM equations have been entirely coded by the authors in order to be used for judging the performance of DPM-FEM code for modeling space harmonics. This section is dedicated to validate the accuracy of the coded timedomain FEM by comparing its results with a commercial FEM software (FEMM). The studied machine is a winddriven doubly fed induction generator whose dimensions are listed in Table 1.
The meshing process in the written code is based on a manually distributed nodes for a single stator/rotor slot which can be then copy-rotated to get the full mesh of stator and rotor core. The airgap is automatically re-meshed at each rotor position using Delaunay Triangulation technique. Fig.  1, 2 show the automatically generated slots mesh performed by the commercial software and the manual meshing of stator and rotor slots in the authors' code respectively. The meshing process in the written code sums up to 27180 triangular elements connecting 15570 nodes.
The two simulators have been run at test values for stator and rotor phase currents as: amplitude of stator current = 10A, frequency of stator = 50Hz, amplitude of rotor current = 20A and a 0.5 slip, and the results are shown in the coming figures. Fig. 3 and 5 show the commercial software's equipotential lines distribution and magnetic flux density respectively, while Fig. Fig. 4 and 6 show the corresponding results generated by the authors' code. Fig. 7 shows the space distribution of airgap flux density at the edge of stator bore and Fig. 8 shows the harmonic analysis of the airgap flux density. The results show a high degree of correlation between the written code and the commercial software which validates the accuracy of the written FEM code. The minor mismatch in the results comes from the difference between the manual meshing of the authors' code and the automatic meshing of the commercial software package.

IV. STUDIED SYSTEM
The studied system consists of a DFIG connected to the grid via a 10Km transmission line and a transformer. The stator's three-phase winding are connected directly to the external system while the rotor winding are interfaced to the grid by a back-to-back power electronic converter used to control the rotor voltage and power. A schematic diagram of the studied  system is shown in the following figure: The DFIG's ratings are: 115 kW, 3.3 kV, 10 poles, 50 Hz. The transmission line's length is 10 km (short line). Therefore, transmission line is modelled as a series RL branch (resistance = 0.35Ω/km, reactance = 0.4Ω/km). The transformers are also modelled as a series RL circuit (%resistance = 1.23, %reactance = 4).
In time-domain, the per-phase electrical equation describing the studied system in electromagnetic transient state takes the form: This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. resistances and inductances respectively, [Q] is a (6 × no. of FE nodes) matrix representing the winding size and direction, A is the nodal magnetic vector potential, l is the machine's stack length and the term (l[Q][ ∂A ∂t ]) represents the backemf induced in the machine's winding. In dynamic phase modeling, the system electrical equation can be written as: where ⟨ ⃗ A⟩ k , ⟨ ⃗ I ph ⟩ k and ⟨ ⃗ V ph ⟩ k are the k th harmonic dynamic phasor of magnetic vector potential, phase currents and input voltages respectively. In order to co-simulate the studied power system and finite element modelled DFIG, the aforementioned system electric equation has to be solved with the FEM magnetic equation directly to get the resulted nodal magnetic vector potentials and phase currents dynamic phasors for a certain applied voltage and slip value.

V. MODELING WINDING SPACE HARMONICS USING DPM-FEM
In order to quantify the harmonic content of DFIG's stator current as an indicator for the generator's power quality in time-domain, the resulted time-stepped waveform of stator phase current simulated by the transient co-simulation of the machine's FEM model and the system electrical equation can by analysed using Discrete Fourier Transform (DFT) (or Fast Fourier Transform (FFT)) to get the magnitude and phase of the different harmonic frequencies contained in the waveform. On the other hand, Dynamic Phasor FEM calculates the equivalent steady-state frequency-domain value of current for the chosen time window. Therefore, using the regular DFT or FFT to analyze the frequency spectrum of the stator current waveform taken from DPM-FEM is not viable. In this section an alternative technique to model the space harmonics reflection into the phase current's time harmonics using DPM-FEM is proposed. The harmonic frequencies found in the stator's phase current apart from supply-frequency are caused by the induced emf's driven by rotor currents. Rotor current typically has two components: the forced current supplied by the rotor's voltage supply and the induced current driven by emf's generated in rotor coils by cutting the magnetic field created by stator current.
Firstly, the space harmonic representation of the airgap flux density waveform caused by the forced rotor current in rotor coordinate reference takes the form [32]: where B r and ⃗ B h r are the airgap flux density and its h th harmonic respectively, s is the rotor slip, θ r is the rotor mechanical position related to the speed ω r = (1 − s) ωs p where ω s is the stator operating angular frequency and p is the pole-pairs. Transferring this flux representation into stator coordinate system (θ s = θ r + ω r t), the previous equation would take the form: From the previous equations, a current of frequency (sω s ) in the rotor winding would induce space harmonics of the orders h in the airgap flux density waveform, consequently, it would be translated into speed-dependent time harmonics of the order s + h p (1 − s) in the stator's induced emf and current. On the other hand, a current of ω s frequency forced in the stator winding by its supply would produce an airgap magnetic field in the stator coordinates of the form: transferring this into rotor coordinates to get the airgap harmonics induced in rotor winding and then transferring it back to stator, the reflected-back induced flux density harmonic compilation takes the form [32]: This implies that a fundamental stator current of ω s frequency will induce space harmonics of order h in the airgap flux density which reflects back to stator (by rotor induction) time-harmonics of orders (1 − h−β p (1 − s)). Therefore, the harmonic frequencies induced in the stator current due to space harmonics of airgap field caused by the fundamentalfrequency excitation of stator and rotor takes the form: where k = 0, 1, 2, 3, ..., these harmonic orders represent both the influence of rotor forced excitation and reflection-back from stator forced excitation for healthy/symmetric machine. If higher time harmonics of stator excitation of order γ are considered, the space-harmonics induced time-harmonics in stator current takes the form [33]: Based on the previous derivation of the harmonic indices to be observed, the proposed procedure for modeling space harmonics in DPM-FEM environment is summarized in the following steps: 1) The DFIG's DPM-FEM model (considering core nonlinearity) is solved simultaneously with the external system equation as a voltage-source powered system using (8) and (11). This step calculates the fundamental stator and rotor currents along with the calculated values of FE elements' magnetic reluctivity after Newton-Raphson convergence.
2) The nonlinear values of magnetic reluctivities are frozen to be used for next steps assuming that the VOLUME 4, 2016 saturation level is determined by only the fundamental components of stator and rotor currents.
3) A current-source powered DPM-FEM solution is carried with rotor currents set to the fundamental values taken from step-1 while stator currents are set to zero to represent open-circuit. The solution is carried out considering linear core at the frozen reluctvities found in step-1 to represent saturation level set by the fundamental current.
4) The nodal magnetic vector potential calculated by step-3 is used to determine the airgap flux density waveform at the inner radius of stator using: where ⃗ B ag is the complex-valued radial flux density calculated at the discrete nodes i of the inner stator radius (r ag ), ⃗ A i is the magnetic vector potential at node i and θ i is node i mechanical position in radians. The resulted flux density waveform can be expressed as a summation of Fourier Coefficients of order h in complex numbers as: where the complex-valued Fourier Coefficients can be determined by integrating over a pole-pair as: using this equation, the real and imaginary values can be used to find the magnitude and angle of each harmonic field to represent it by a traveling sinusoidal wave that rotates around the airgap with equivalent number of pole-pairs equals p/h. 5) The space distribution of each harmonic field determined in the last step can be used to calculate the induced emf of each stator slot ( ⃗ E h slot ) using the Flux Cutting Rule (Blv) as: where ⃗ B h ag (θ slot ) is the value of space harmonic flux density at the slot position (θ slot ) and v h is the linear velocity of the harmonic field travelling wave. This EMF value can be applied to the external system impedance considering the mapping frequency calculated by (16) to determine the corresponding value of stator harmonic current. The utilization of flux density at stator bore to calculate winding EMF is viable here because the machine has an open-slot configuration.
6) The harmonic stator currents found in last step are used one at a time to determine its own induced harmonics back into stator winding after being applied to the rotor winding. A linear current-source FE solution is carried out like previous steps while forcing rotor voltage to zero (short-circuited rotor winding) at each stator harmonic field. The airgap flux density wave at stator inner diameter is determined, analyzed into Fourier coefficients and the corresponding slot EMFs are calculated using the same approach of step-4. The resultant EMFs are used as a source for the external system impedance at the mapped frequencies found by (17) to determine the harmonic currents caused by stator excitation.
7) The harmonic stator currents occurring at a matched frequency from rotor excitation and stator excitation are summed to determine the final value.

VI. SIMULATION RESULTS AND DISCUSSION
In order to prove the validity of the presented DPM-FEM method for modeling space harmonics effect on statorcurrent's frequency spectrum, the simulation results of the written code must be close enough to the corresponding timedomain FEM code while significant reduction in DPM-FEM simulation time is to be noticed. The codes have been run on a PC with Intel Xeon CPU v3 -3.3 GHz with 16.0 GB installed RAM. The system under study has been previously presented. The DFIG's stator is connected to the external system directly to be fed by three-phase balanced voltage supply at the machine's rated value. The rotor supply is modeled as a controlled voltage source to represent the backto-back converter. In order to validate the proposed procedure of considering the speed-dependent space harmonics in the DPM algorithm, the DFIG is simulated under a timechanging rotational speed. The working slip is shown in 11. The simulation results of the coded DPM and time-domain FEM of DFIG while connected to the adopted external power system are shown in Fig. 12 to 21. Fig. 10 shows the magnetic flux lines in a quarter of the machine proving that DPM-FEM gives the same flux pattern of the original time-domain one (which has been bench-marked against a commercial software earlier in Section-3). In order to test the validity of the proposed procedure of modeling the winding space harmonics influence on the harmonic content of stator current (as an indication for the generator's power quality), the proposed procedure is tested at the two different steady-state speeds corresponding to rotor slip chnaging from 0.5 to 0.25. Fig. 16 shows the airgap magnetic flux density over the span of a single pole-pair calculated by both time-domain and DPM-FEM at 0.4s time instance, a high degree of matching between the two solvers 8 VOLUME 4, 2016 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. can be noticed. In order to further prove the two solvers matching, fig. 17 shows the matched space harmonic analysis of the flux density waves resulted from the two adopted solvers. This airgap flux density spectral components are used to estimate the corresponding time-harmonics in the stator current's waveform using the proposed procedure in Section-V. Each harmonic flux represents a travelling wave with a specific magnitude, speed or rotation and direction, the flux cutting rule is used to determine the resulted induced voltage which is solved simultaneously with the external circuit to estimate the corresponding current harmonic. Fig.  18 shows the resultant frequency components of stator line current determined by FFT of the time-domain FEM resulted current and the values calculated by DPM-FEM's adopted procedure. A high degree of matching can be observed. Moreover, the same results are taken again at 1.0s time instance to prove the validity of the proposed procedure for capturing the speed-dependent harmonics. Fig. 19 to 21 show the resulted airgap flux density per pole-pair, spectral harmonic analysis of these flux waves and the resulted timeharmonics in stator line current calculated by the two solvers respectively. A high degree of matching continues to be observed with the speed-change.
The comparison of the percentage total harmonic distortion (%THD) of the stator current waveform -as a power quality indicator-between the proposed DPM-FEM and the original time-domain FEM is given in Table 2 along with the computational performance of the two solvers.
The included simulation results clearly show that the proposed DPM-FEM accurately tracks the results of the time-domain counterpart, which validates the convenience of using the proposed technique for co-simulation studies of machines and power system from the accuracy point of view. In addition, Table 2 presents the comparison of the mathematical performance of both solvers. Firstly, the proposed technique for modeling space harmonics in DPM-FEM gives a %THD value of 4.603 (at 0.4s) which is only 4.28% deviated from the time-domain counterpart, and %THD value of 4.822 (at 1.0s) which is only 3.79% deviated from the time-domain counterpart. The nature of DPM allows using much larger simulation step-size of 0.005 which is 50 times larger than the time-domain counterpart which leads to a significant reduction in total simulation time by 84.97% of the time-domain FEM solver. The high degree of accuracy paired with the superior computational performance, proves the validity of using the proposed DPF-FEM technique for co-simulation studies of grid-connected generators considering winding space harmonics representation. The minor mismatch seen in the airgap flux density calculated by both DPM-FEM and time-domain FEM shown in Fig. 16 and 19 and the resulted minor mismatch in stator current's timeharmonics, is due to the absence of speed-induced slotting harmonics in DPM-FEM as the rotor is kept stationary in DPM-FEM while being rotational in time-domain FEM. This can be clearly seen through the perfect matching of the airgap flux density of the two solvers considering a blocked-rotor condition shown in Fig. 22. This very small mismatch is acceptable for the sake of reducing time consumption and computational burden achieved by DPM-FEM.

VII. CONCLUSION
This paper investigates the accuracy and validity of using dynamic phasor finite element method for co-simulation studies of DFIG and a connected power system. More concern VOLUME 4, 2016 9 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and     has been directed to modeling the winding space harmonics impact on the harmonic content of the machine's sator current pushed to the grid as an indicator of the machine's     has been introduced along with the basic formulation of the traditional time-domain FEM. Both solvers have been coded by the authors and the accuracy of the proposed DPM-FEM method including modeling space harmonics has been proved by matching the simulation results of the traditional time-domain FEM code which has been previously benchmarked against a commercial FEM software package. The computational superiority of DPM-FEM has been noticed by the ability to be performed at 50 times larger step-size leading to a more than 6.5 times reduction in simulation time.