Dynamic Phasor Finite-Element Modeling of a DFIG for Grid Connection Studies

Cosimulation studies of electric power systems and electric machines have always been a challenge. In order to reduce the simulation time to a reasonable value, lumped-parameter electric machine models are commonly used in electric power system modeling software packages to avoid the heavy computational burden of more accurate modeling methods especially finite-element method (FEM) on the expense of less accuracy. The proposed technique in this article combines the dynamic phasor modeling technique for power system simulations with the FEM to accurately model the doubly fed induction generator while connected to the grid. The utilization of dynamic phasors enables adopting large simulation time steps resulting in a significant reduction in the simulation time compared to the conventional time-domain FEM modeling. The mathematical foundation of the proposed modeling method is presented including the generator's core saturation. Custom-written C++ codes have been developed by the authors to execute the new dynamic phasor FEM algorithm and the conventional time-domain FEM in order to fairly compare their accuracy and numerical performances. As the proposed method combines time and frequency domains, a unique capability of modeling the rotor movement can be achieved. The rotation can be represented by physically incrementing the rotor and airgap mesh as in regular time-domain solvers, by mathematically representing the rotation using the virtual blocked rotor method as in frequency-domain solvers, and the proposed method of combining the two aforementioned approaches. The three methods of modeling rotor rotation are discussed, and their simulation results are compared to give a guide to choose the proper method for the different modeling targets. The results show that the proposed dynamic phasor FEM is capable of producing comparable results to the traditional time-domain solver at a substantially reduced simulation time.

the stator and rotor slots, nonlinear characteristics of the stator and rotor core materials, skin effect, skewing of the rotor bars, end effect of the stator winding, and eddy currents can all be modeled accurately in the FEM. The field equations are solved over the mesh elements, and then, the solution is assembled over the whole nodes to get the global solution of the nodal magnetic vector potential [2]. Then, other parameters, such as inductances, magnetic flux density, winding flux linkage, back electromotive force, and electromagnetic torque, can be calculated in the postprocessing phase [3], [4], [5]. In order to use the accurate FEM to model an electric generator as a component of a larger electric and/or mechanical system, the magnetic field equations of the FEM must be solved along with the power system's electric equations and the prime mover's mechanical equations in a direct or an indirect approach. The direct approach means solving the FEM equations and the external system equations simultaneously [6], [7], [8]. The indirect approach uses the FEM to calculate the machine parameters to be used separately in state-space or lumpedparameter models [9], [10], [11].
In spite of all the aforementioned merits of the FEM in terms of accuracy and flexibility, it is always criticized by being computationally expensive and time consuming. Therefore, using FEM as the machine's modeling method inside a transient time-stepped simulation of an interconnected system requires impractically long time and computational power as the electromagnetic FEM equations must be solved entirely at each time step, which is in the very low range of microseconds in order to capture high-frequency components. Several research efforts have been presented in the literature to reduce the simulation time of the time-stepped FEM to facilitate the cosimulation studies of electric generators and their connected systems. Salon et al. [2] used the axisymmetric boundary condition of the FEM to model a one pole pitch of the machine. This reduces the size of the problem significantly. However, it works only for symmetric and healthy conditions of the machine. Koti et al. [12] used the frequency-domain steady-state equivalent circuit to provide an initialization of the transient FEM solution. Koti et al. [13], [14] used a blocked-rotor frequency-domain FEM to do the same initialization. However, these methods work only if the starting behavior is out of concern. Drobnič et al. [15] used precalculated FEM's flux-current data of a permanent magnet machine to use it for online simulation in order to reduce the simulation time of the connected system. Ionel and Popescu [16] used the FEM to model interior permanent magnet motors. In order to reduce the simulation time, the authors took advantage of the winding symmetry and modeled only three rotor positions in magnetostatic environment to construct the flux linkage waveform of each phase. Finite-element postprocessing along with analytical equations is used to estimate the motor torque and core losses. Miller et al. [17] used multiple offline FEM simulations to construct a flux-magnetomotive force diagram for the fast estimation of average torque and winding inductances. Di et al. [18] proposed a procedure for reducing the simulation time required by the FEM over three stages: 1) solving a locked-rotor current-source simulation under timeharmonic analysis with an estimated current value taken from an initial time-harmonic FEM solution or analytical methods; 2) voltage-source solution with time-harmonic analysis; and 3) switching to normal rotating solution to get the actual performance. The results show good agreement with the traditional FEM solution. However, using an estimated value of the starting current sacrifices the optimum FEM accuracy. Fu et al. [19] proposed a technique to reduce the simulation time taken by the Newton-Raphson iteration for nonlinear core representation. The method depends on automatically adjusting the calculation tolerance of the current time step by a prediction based on the information taken from the last one in order to reduce the required iterations as much as possible. The authors of [20], [21], [22], and [23] used model order reduction techniques to reduce the simulation time associated with the FEM. However, they lack the optimum accuracy especially for magnetic field distribution and use specific snapshots data to construct the model that risks the solver's accuracy. This article uses a combination of the dynamic phasor modeling (DPM) technique and the FEM to accelerate the precise modeling of a grid-connected DFIG.
DPM, based on a generalized averaging technique, is an attractive method of modeling systems, in which the small number of harmonics is enough to describe the system behavior accurately. DPM approximates the time-domain signal over a sliding window by a summation of complex-valued Fourier coefficients called dynamic phasors. Since the dynamic oscillation of dynamic phasors is much slower than the corresponding time-domain signals, much larger integration time steps can be utilized. Therefore, dynamic-phasor-based simulation has the advantage of shorter execution time and less overall computational burden compared to the detailed time-domain simulations even though DPM uses a higher number of equations in addition to be carried out using complex math.
The DPM technique has been used in several research papers to model the performance of DFIG [24], [25], [26], [27], [28]. These papers and much more, which used DPM to model electric machines especially grid-connected DFIG, used lumped-parameter models in dqo or abc frames. As discussed before, lumped-parameter models lack the required accuracy of physically modeling the machine topology and nonlinearities. Therefore, this article tries to fill in this gap by using dynamic phasors to represent the state variables of the electromagnetic equations of the FEM and the electrical equations of the connected power system in order to make use of the FEM's well-known accuracy while reducing the simulation time to a practical value. The rest of this article is organized as follows. Section II reviews the traditional formulation of the time-stepped FEM in transient time-domain environment considering both linear and nonlinear core representation. Section III discusses the proposed dynamic phasor FEM (DPM-FEM) formulation considering core nonlinearity. Section IV presents the simulation results of the starting behavior of a grid-connected DFIG in order to prove the validity of the proposed method in terms of accuracy and time consumption. Finally, Section V concludes this article.

II. FORMULATION OF TIME-STEPPED FEM IN TIME-DOMAIN ENVIRONMENT
This section presents a brief review of the basic formulation of the transient time-domain FEM. It starts from the electromagnetic diffusion equation [2] × where ν is the magnetic reluctivity, A is the magnetic vector potential, J is the external applied current density, and σ is the electric conductivity. The term (σ ∂A ∂t ) represents the eddy current induced in the winding by cutting the time-varying magnetic field. The utilized rotor and stator winding of the DFIG are of stranded wires with small cross-sectional areas, which reduces the induced eddy currents to a negligible value. In the 2-D FEM, the magnetic vector potential and the applied current density are assumed to vary only in the Cartesian x, y dimensions as The 2-D electromagnetic solution is achieved by applying the Galerkin method of weighted residuals, which matches the shape functions to the weighting functions. Consequently, the diffusion equation representing the magnetic vector potential over a single mesh element takes the matrix form of where [S e ] is the elemental stiffness matrix whose element S e,i j calculated as S e,i j = (b i b j + c i c j )/4 e , e is the element surface area, and b i, j and c i, j are the y and x coordinate difference of the element nodes, and T e,i j equals e /6 for i = j and equals e /12 for i = j. The local stiffness matrix and the local excitation vector of each element are used to construct the global system equation. The contributions of different elements that share certain node/nodes are locally summed to determine the corresponding global value for each node.
Equation (5) is valid for calculating the magnetic vector potential of current-source problems if the machine's iron core is assumed to have linear BH characteristics, i.e., the magnetic reluctivity ν is constant and independent of the magnetic vector potential throughout the solution process. In order to include core nonlinearity in the calculation, the magnetic reluctivity ν is considered a function of the local magnetic vector potential over the considered mesh element. The Newton-Raphson method is used to solve the FEM equation iteratively to determine the nonlinear value of ν using where n is the Newton-Raphson iteration index and [G e ] is calculated as where B is the magnetic flux density of the element. The value (∂ν/∂B 2 ) represents the rate of core-reluctivity change with the change of the flux density. The magnetic reluctivity is initialized by its linear value, and the global system of equations is solved for the corresponding magnetic vector potential [2]. Then, the local magnetic flux density is calculated in postprocessing and used to determine the corresponding nonlinear value of ν using the actual BH data of the core material. The cubic spline method is used to interpolate the BH data to a continuous function to be used in the calculation. The process is repeated till convergence. The system equation can be solved as Ax = b system using numerical techniques (conjugate gradient method has its potential in this matter) instead of trying to invert the very large global stiffness matrix. Fortunately, the system matrix is symmetric and sparse, which can be used to optimize the numerical solution to reduce the required memory and solution time as much as possible.
The calculated magnetic vector potential can be used in the postprocessing phase to calculate other electromagnetic quantities, such as the winding flux linkage, which can be divided by the input current to determine the winding inductance and differentiated to determine the induced back electromotive force (EMF). In addition, the magnetic vector potential can be used to determine the elemental magnetic flux density, which can then be used to determine the developed electromagnetic torque using the Maxwell stress tensor approach.

III. FORMULATION OF THE PROPOSED DYNAMIC PHASOR FEM
The time-stepped FEM performed in a time-domain environment has the superior ability of capturing all the physical phenomena such as the winding space harmonics, speed-induced slotting harmonics, and electromagnetic torque ripples. However, in order to capture all the frequencies of concern, a very small time step has to be used, which is translated into a small increment in the rotor position from a time instance to the next. Therefore, it would take a very long simulation time to pass the starting transients and reach the steady state. The purpose of this article is reducing the simulation time required for the time-stepped FEM by applying the DPM technique to the FEM's electromagnetic equations along with the equations describing the behavior of the connected power system. This section discusses the proposed dynamic phasor FEM approach and derives its equations in detail. First, the steady-state frequency-domain FEM is reviewed to be used as a base for the DPM-FEM. Second, the basic principles of the DPM technique are reviewed. Finally, the formulation of the proposed DPM-FEM is discussed.

A. STEADY-STATE FREQUENCY-DOMAIN FEM
For various modeling applications, just the steady-state fundamental component of phase current, magnetic flux density, power loss, and the average value of the electromagnetic torque is the target. For this scope of applications, the steadystate frequency-domain representation of FEM equations is the choice. In spite of being performed in a complex-math representation, which is more complicated than real-valued time-domain equations, it solves the system equations once for each harmonic frequency to get the steady-state solution. In addition, it does not require any rotor rotation; instead, the impact of the rotor speed is taken into account using what is called virtual blocked rotor (VBR) method. Just like the approach followed in the well-known induction motor's steady-state equivalent circuit, the rotor winding resistance and the rotor input voltage are divided by the rotor slip to account for the speed of rotation. The 2-D elemental diffusion equation in the frequency domain takes the form where A e and F e are complex-valued vectors of A e and F e , respectively, ω e is the system frequency in rad/s, and j is the complex operator. This equation is used for each frequency of concern one at a time.

B. PRINCIPLES OF DPM
Fourier series has been used for decades to represent periodic signals by a summation of sinusoidal waves (with constant magnitudes and phases) in order to take advantage of their ease of differentiation and integration. The DFIG's stator and rotor currents, flux linkages, and back EMFs are periodic signals once they reach the steady state. However, during transients, these signals are not perfectly periodic but can be considered almost periodic, which prevents the utilization of the regular Fourier series to decompose those signals. DPM can be used for this purpose by approximating the timedomain signal with a special Fourier expansion calculated over a time window that slides over the considered time range.
As the modeled signal is not perfectly periodic and as a sliding time window is used, the magnitudes of the Fourier coefficients are time varying. The modeled time-domain signal x(τ ) modeled over a time window τ ∈ (t − T, t ) takes the form [25] x where where T is the fundamental period in seconds, ω is the fundamental frequency in rad/s, X k (t ) is the kth Fourier coefficient, and K is the chosen set of dynamic phasors to approximate the actual wave. The value X k (t ) is a complex-valued timevarying coefficient, called the dynamic phasor and written as x k (t ). The dynamic phasor must be treated in a different way regarding the mathematical operations due to its unique nature. The derivative of x(τ ) is calculated in terms of dynamic phasors as dx dt k = dX k dt + jkωX k (11) and the product of two time-domain signals x and y is done using the convolution theorem for dynamic phasor representation as DPM offers a better tool than the regular frequency-domain analysis in terms of flexibility as it can model multiple frequencies simultaneously. In addition, it can be used to investigate the mutual effects of certain frequencies without the need for a full-detailed time-domain analysis. The most important advantage of dynamic phasors is the slower variation compared to time-domain signals, which allows the utilization of larger simulation time steps reducing the simulation time significantly. However, choosing the proper set of harmonics is vital for optimizing the performance of DPM. Considering too many harmonic components leads to a complicated solution and long simulation time, while considering fewer harmonics leads to less accuracy. Therefore, a good knowledge of the expected frequency components in the modeled system is required for compromising the accuracy and complexity. In general, DPM is more mathematically complicated than the corresponding time-domain analysis. The reason is that the original number of system equations is multiplied by the number of required harmonics. In addition, it has to be performed in complex math. However, the advantage of the slow variation and its allowance of using much larger time step outweighs the complexity, and the overall time-stepped solution in DPM is usually significantly less time consuming.

C. DYNAMIC PHASOR MODELING FEM
The aforementioned electromagnetic equations of the FEM in the time domain are rederived here for DPM. The state variables of the FEM equations, the nodal magnetic vector potentials and the phase currents, are modeled as dynamic phasors, and the time derivatives are carried out using the DPM derivative property given in (11). The elemental electromagnetic diffusion equation in DPM can be derived as where A e k and F e k are the complex-valued kth-harmonic dynamic phasors of A e and F e , respectively, and ω k is the harmonic angular frequency. This equation is used as many times as the required harmonics in order to fully represent the system. If just the fundamental frequency is considered, the FEM equations are solved after shifting the system frequency backward by the fundamental value (which becomes a dc component); this is called "shifted frequency analysis" [29].
Using the Crank-Nicolson method, the previous equation can be time-discretized in order to be used for the timestepped solution using complex-valued time-varying dynamic phasors as where ν is the linear value of the magnetic reluctivity. A complex-valued magnetic reluctivity can be used to approximately account for the magnetic hysteresis, but it is generally ignored without much loss of accuracy [2]. For the timedomain FEM, the actual BH data of the core material are followed point by point for each time step to determine the corresponding value of the magnetic reluctivity for a certain saturation level. For the DPM-FEM, the equations are solved in a time-frequency environment. Therefore, a single ac value of the magnetic reluctivity is chosen to represent a whole cycle of alternation. In the traditional frequency-domain FEM, all the variables (voltages, currents, and flux) are considered sinusoidal. However, for nonlinear core, even if the magnetic field intensity H is sinusoidal, the corresponding flux density B is not. Therefore, a special transformation is done for the actual BH data in order to be used for the frequency and time-frequency FEM. The BH data are converted so that the effective ac magnetic reluctivity gives the same average energy as the actual one [30]. First, H is assumed sinusoidal, and the corresponding B is obtained from the actual BH curve. Then, the energy E in a linear core is obtained by integration over a quarter of a cycle (T ) as For each peak value of sinusoidal H, the effective ac magnetic reluctivity can be calculated as Fig. 1 shows the dc magnetic reluctivity calculated from the actual BH data and the resulted ac effective reluctivity for the core material used in the modeled DFIG. The cubic spline method is again used to interpolate the resulted data in order to be used for the DPM-FEM solver.
The calculated dc reluctivity is used with the Newton-Raphson method to model the core nonlinearity. The matrix form of the time-discretized DPM-FEM equation is derived as where [G f ] is calculated as where A e * is the complex conjugate of the dynamic phasor of the magnetic vector potential and | B| is the absolute value of the complex-valued elemental magnetic flux density. The saturation level of the machine affects the number of iterations required for the Newton-Raphson solver to converge. Therefore, the highly saturated machines usually take more simulation time than the lightly or unsaturated machines. However, this will not affect the time-saving factor expected with DPM as the effect of the saturation level on the convergence rate is present in both time-domain and DPM-FEM solvers.
The effect of rotor speed is modeled in the conventional time-domain FEM by the actual position increment of the rotor and airgap mesh from a time step to the next, while in the steady-state frequency-domain FEM, the effect of the rotor speed is taken into account using the VBR technique, as illustrated before. Therefore, no actual change in the rotor position is modeled. For the proposed DPM-FEM, as the technique incorporates both a frequency-domain calculation and a time-stepping behavior, the effect of the rotor speed can be modeled by both the two aforementioned ways. The rotor and airgap mesh can be kept fixed, and the VBR approach is used to model speed effect by dividing the rotor resistance and input voltage by the rotor slip. This can save the time and computational power associated with changing the coordinates of rotor mesh nodes and remeshing the airgap for each time step on the expense of sacrificing the accuracy of modeling speed-induced current harmonics. On the other hand, the speed effect can be represented by incrementing the rotor position with the increment in time just like time domain in order to model the starting subharmonics and the speed induced current harmonics on the expense of a longer simulation time.

IV. SIMULATION RESULTS AND DISCUSSION
The validity of the proposed dynamic phasor FEM is tested in this section by comparing its results and computational performance against the conventional time-domain FEM. As the proposed method essentially modifies the conventional FEM's electromagnetic equations (in time and frequency domains), the authors developed a custom-written C++ code to execute the proposed DPM-FEM algorithm taking a grid-connected DFIG as the modeling target. In addition, the conventional time-domain FEM algorithm for linear and saturated core has also been coded by the authors to facilitate a fair comparison. The case study taken here is modeling the starting performance of a DFIG when connected to the public grid via a transmission line. Modeling the starting performance of electric generators using the FEM has always been a challenge due to the rapid oscillations caused by various speed-induced subharmonics in the machine's current and torque, which requires long simulation time if the conventional time-domain FEM is used. The following subsections present the simulation and computational performance results.

A. SYSTEM UNDER STUDY
The system used as a case study for testing the proposed DPM-FEM technique for simulating a grid-connected DFIG is shown in Fig. 2.
The stator of the simulated DFIG is connected to the grid via a transmission line and a transformer, while the rotor is interfaced to the grid by a back-to-back converter used to regulate the rotor power and voltage. The machine ratings are 115 kW, 3.3 kV, ten poles, and 50 Hz. The length of the  transmission line is 10 km, which is in the range of short lines. Therefore, transmission line is modeled as a series RL branch (resistance/km = 0.35 , reactance/km = 0.4 ). The transformers are also modeled as a series RL circuit (%resistance = 1.23, %reactance = 4.0). The rotor supply is modeled as a controlled voltage source to represent the back-to-back converter control action.
The authors' developed codes for the time-domain and dynamic phasor FEM considering both linear and saturated core have been executed on a PC with Intel Xeon CPU v3-3.3 GHz with 16-GB installed RAM. The machine's 2-D FEM model is discretized to 27 180 triangular mesh elements connecting 15 570 nodes. Discretizing the machine's 2-D cross-sectional area into the mesh elements can be done automatically using various mathematical methods; the most widely used one is Delaunay triangulation method [31]. For the sake of flexibility, the main mesh of the modeled DFIG is built manually in the authors' code by setting the node and element distribution for a single stator and rotor slot pitch and copy-rotating them to form the whole stator and rotor mesh. Fig. 3 shows the stator and rotor mesh adopted in the code. The airgap is remeshed automatically using Delaunay triangulation for each rotor position in the time-stepped solutions.
For grid-connection studies, the FEM electromagnetic equations are solved simultaneously with the electrical equations, which describe the interaction between the grid voltage and the machine back EMF. The time-domain grid voltage equation takes the form whereṼ ph ,Ĩ ph , andÃ are the complex-valued static vectors of phase voltages, currents, and magnetic vector potential, respectively, and ω e is the angular frequency.
In the dynamic phasor domain, the external grid voltage equation is derived as where Ã k , Ĩ ph k , Ṽ ph k , and ω k are the kth harmonic dynamic phasor of the magnetic vector potential, phase currents, input voltages, and harmonic angular frequency, respectively. The following subsections present the simulation results of the starting behavior of a grid-connected DFIG using both VBR and position incremental methods for modeling rotor speed.

B. CHOOSING THE SAMPLING RATE
The rule of thumb deduced from the literature is using 100 samples per cycle for the position-incremented stepped timedomain FEM. The custom-written DPM-FEM code of the time-domain solution is first used to test the performance of the machine under different integration samples in order to demonstrate its impact on the numerical accuracy.
The starting behavior of a grid-connected DFIG is investigated using the time-domain DPM-FEM under different integration samples starting at 200 samples per cycle down to just 4. The dynamic variation of the DFIG's electromagnetic torque is shown in Fig. 4, and a zoomed portion of the figure  for the starting and steady-state periods is shown in Figs. 5 and 6.
It can be seen from the figures that the results are almost the same for 100 samples per cycle or more. However, reducing the sampling value to a small number such as 20 samples per cycle or less leads to a significant degradation of the modeling accuracy not just in the harmonic content but also in the average steady-state value. Therefore, the least number of samples that can be adopted to keep the optimum accuracy is 100 samples per cycle if the harmonic content is of concern, and is 20 samples per cycle if just the steady-state average values are of concern.

C. SIMULATING THE STARTING BEHAVIOR OF A GRID-CONNECTED DFIG USING THE VBR METHOD
In order to test the validity of the proposed DPM-FEM algorithm for modeling the DFIG starting behavior, the modeled DFIG is simulated using both the time-domain-FEM-and the DPM-FEM-developed codes considering the core saturation effect. The time-domain FEM always models the rotor speed by incrementing the position of the rotor mesh according to the speed value and time step while automatically remeshing the airgap. In the first test case, in order to model the rotor speed in the DPM-FEM environment, the VBR method is used to represent the required effect while keeping the actual rotor and airgap mesh stationary throughout the simulation. The effect of the rotor speed is included by dividing the rotor winding resistance and voltage excitation by the rotor slip.
The starting behavior at the rated voltage of a gridconnected DFIG (which is always hard to model) is the    FEM's simulation time of 25 848.18 s. The reduction in simulation time is not at the same order of the increase in the time step because the DPM technique is more mathematically complicated than the time-domain solver.
A hundred samples have been adopted in the time-domain FEM solver to represent the full capacity of the solver's modeling accuracy, which is the reason of bearing the struggle with its associated high time consumption and computational burden, while the DPM-FEM has been carried out at the minimum allowed number of samples for representing the required harmonics (which is only the fundamental here). However, if only the steady-state fundamental values of currents and flux linkages along with the average value of the developed torque are the main concern (not the harmonic content), the minimum number of samples can be adopted for the time-domain solver, as previously discussed (20 samples per cycle). This brings down the simulation time of the time-domain FEM to 6908.54 s, which brings down the time-saving percentage of the DPM-FEM to 86.64% of the time domain after it was 96.43% with 100 samples per cycle. Even with using the least viable sampling rate for the timedomain FEM, the DPM-FEM using the VBR method is still not accurate enough in modeling the speed-induced starting and steady-state harmonics especially when only the fundamental component is considered. However, it can be a perfect fit for modeling the fundamental current, power losses, and average torque variations associated with the change in the system variables such as input mechanical torque, wind speed, and grid voltage in a significantly reduced simulation time.

D. SIMULATING THE STARTING BEHAVIOR OF A GRID-CONNECTED DFIG USING THE INCREMENTED ROTOR POSITION METHOD
In order to model the speed-induced harmonics and subharmonics causing the starting ripples, the rotor and airgap rotation is modeled physically by incrementing the rotor position from a time step to the next according to the prime mover speed. The coordinates of the rotor mesh nodes are adjusted, and the airgap is remeshed for each new position. This is the typical method used with the conventional time-domain FEM. In the DPM-FEM, the incremented-position method is carried out by modeling the machine variables by time-varying Fourier coefficients at each time step. In frequency-domain representation, the angular phase of the rotor excitation has to be adjusted according to the corresponding rotor position in order to simulate the actual distribution of current in the rotor winding. In order to match the accuracy of the traditional time-domain FEM, a relatively smaller time step has to be adopted in order to capture the speed-induced harmonics.
Figs. 12-16 show the simulation results of the DPM-FEMdeveloped code and the corresponding time-domain FEM results for the stator phase current, rotor phase current, stator winding flux-linkage, rotor winding flux linkage, and electromagnetic developed torque, respectively. The time-domain code has been executed at the least possible sampling rate of 20 samples per cycle (0.001-s time step). In order to be able to capture the same harmonics, the DPM-FEM code has also been executed at 20 samples per cycle.     As expected, if both the solvers are executed at the same sampling rate, the DPM-FEM faces a disadvantage. The numerical nature of the DPM-FEM makes it more complicated as it requires dealing with complex numbers, which calls for using the "complex" class of C++, which is much slower compared to the standard variable containers, or the need for separating the real and imaginary parts, which doubles the system size. Therefore, the DPM-FEM takes 8474.84 s, which is 22.7% more than its time-domain counterpart.
However, the DPM-FEM results almost perfectly track those of the traditional time-domain FEM, which proves the accuracy of using the regular incremented rotor position method with the DPM-FEM just as in the traditional timedomain FEM. This would be the base of using it in a mixed approach with the VBR method in the following subsection.

E. SIMULATING THE STARTING BEHAVIOR OF A GRID-CONNECTED DFIG USING A MIXED APPROACH BETWEEN THE VBR AND THE INCREMENTED ROTOR POSITION METHODS
The previous subsections discussed the utilization of the VBR and the incremented-position methods for modeling the rotor speed in the DPM-FEM. The VBR method saves the time and computational power required for remeshing the rotor and airgap while being fairly accurate for modeling the steady-state phase. However, the VBR method is not accurate enough for modeling the starting harmonics. This can be solved by using the incremented-position method on the expense of more simulation time and computational power. In order to compromise the time saving of the VBR method in the steady-state phase and the accuracy of the incremented-position method in the starting phase, a mixed approach between the two methods is proposed in this subsection.
First, the machine is modeled in the DPM-FEM environment using the incremented-position method at a relatively small time step of 0.001 s (20 samples per cycle). The rotor position is physically incremented from a time instance to the next while using the actual values of rotor resistance and applied voltage in order to represent the starting harmonics and subharmonics. Once the starting period is passed, the rotor position and the corresponding phase angles of rotor excitation are frozen at their last values to be used with the VBR method for modeling the steady-state phase in order to guarantee modeling continuity. The time step is enlarged from 0.001 s (20 samples per cycle) with the incremented-position method to the higher 0.01-s time step (two samples per cycle) with the VBR method in order to reduce the overall simulation time.
Figs. 17-21 show the simulation results of the DPM-FEM code using the proposed mixed approach between the VBR and the incremented-position methods of stator phase current, rotor phase current, stator winding flux linkage, rotor winding flux linkage, and electromagnetic torque, respectively. The switch between the incremented-position method and the VBR method is applied at 0.3-s instance. 0.3 s is an estimated value for the machine's starting period, which can be determined based on the previous experience depending on the machine design. The estimated value of the starting period affects the accuracy and computational superiority of the proposed method. Therefore, a good prior knowledge of the  studied system is a must for a proper design of the modeling approach.
The results show that the proposed method coupled with the DPM-FEM highly matches the results of the conventional time-domain FEM in both starting and steady-state phases while being executed at the total simulation time of 4027.74 s, which is only 58.3% of the corresponding time-domain FEM execution time.

V. CONCLUSION
In this article, a new modeling technique for the gridconnected DFIG was proposed by merging the dynamic phase modeling technique with the electromagnetic FEM. The formulation of the proposed method was derived after reviewing the conventional time-domain FEM. The new method took advantage of the shifted-frequency property associated with the dynamic phasors to significantly reduce the simulation time via enlarging the integration time step. The DPM-FEM modeling technique for the grid-connected DFIG was used to model the starting behavior of a grid-connected DFIG. Three different methods for modeling the rotor speed (VBR,  incremented-position method, and a mix between the two) were compared in order to characterize the optimum method for modeling the starting and steady-state behavior of the grid-connected DFIG for different modeling purposes. The simulation results proved the validity of the proposed DPM-FEM for producing comparable results to the conventional time-domain FEM while reducing the simulation time by 86.64% with the VBR method and by 41.7% with a mixed approach between the incremented-position method (for modeling the starting phase) and the VBR method (for modeling the steady-state phase). The system operator can then choose the proper method according to the level of accuracy and modeling target required by his application in order to optimize the results accuracy while saving the computational power and time. The VBR method provided the maximum time saving while poorly tracking the actual starting harmonics. The rotor incremented-position method gave the highest accuracy on the expense of the longest simulation time. By mixing the two methods, a good compromise can be achieved. The results showed that the proposed DPM-FEM method can be a game changer for the grid integration studies of electric generators.