Finite Element Computation of Transient Parameters of a Salient-Pole Synchronous Machine

: This paper presents a novel method for calculating the transient parameters of a saturated salient-pole synchronous machine (SPSM) with damper cage using ﬁnite element analysis. All detailed leakage inductances in a modiﬁed d-q model are discussed and separately extracted. In addition; the frozen permeability method is used in a 2-D ﬁnite element analysis to consider saturation and skin effects for determining the inductances. The terminal reactance parameters are obtained from all elements of the d- and q-axis equivalent circuits in all the chosen time constants during the transient process. The variation of leakage inductances of the SPSM in transient processes is also investigated and discussed. To determine transient time constant parameters; the Prony algorithm is applied in the presented method. A program developed to automatically solve the simulation and computation with the proposed method is described. This method is applied to a prototype and validated by experimental results. Some discussion about the relationship between rotor geometric design and the transient parameters of SPSMs is provided. The variation of leakage magnetic ﬁeld of SPSMs in transient processes is also investigated. This method can be applied in transient parameter estimation and optimization of SPSMs in the design stage.


Introduction
Transient parameters of electrical machines are generally applied to estimate the dynamic behavior of machines in power system analysis [1,2].The salient-pole synchronous machine (SPSM) is one of the most-used types of machine in industrial application.In hydroelectric plants, SPSMs have been applied as generators since the nineteenth century.The SPSM is also preferred as the type of generator used in ships [3].In traction drive, mine hoist, and marine propulsion, SPSMs are commonly employed as drive motors [4].Besides, SPSMs have a long history of utilization as motor-generators in pumped storage plants and as synchronous compensators in some substations.As most of the SPSMs are connected with the power grid, accurate calculation of the transient parameters is critical in determining the dynamic performance of a SPSM.
In the literature, there are many reports of various tests and methods to estimate the transient parameters of a synchronous machine, which are represented as terminal operational impedances and time constants [5,6].However, these studies mainly focus on a machine that has already been designed or produced.In the design stage, precise determination of the transient parameters of an SPSM is also important for refinement.The designers not only concentrate on the terminal parameters, but also prefer to accurately calculate each element of the magnetic equivalent circuit for design optimization.
The analytical solution is traditionally applied in machine design to calculate the transient parameters from the equivalent circuit and lumped parameter model [7,8].However, the accuracy of analytical formulas is inadequate due to the simplified assumptions and many correction factors made in the calculation procedure.The finite element (FE) method has been recently applied in various SPSM analyses [9,10], which makes the consideration of saturation and skin effects come true.By employing FE analysis, the accuracy of computation, compared with analytical solution, is significantly improved.For example, the standard three-phase sudden short-circuit test is simulated in [9], and the standstill frequency response test is simulated in [10] by using a time-stepping FE analysis to obtain terminal transient parameters.
Though several studies about transient parameter calculation by the FE method have been reported, the existing methods do not fully satisfy the demands of designers.In the traditional method given in [11] and applied in [9], transient parameters are obtained from short-circuit current waveform by wave analysis.This method is relatively subjective, and only suitable for terminal parameters' computation.For FE simulations of standstill frequency response tests [10], it is hard to follow the same magnetization saturation effects in the dynamic process, which would influence the accuracy of parameters obtained by these methods.In [12], an improved curve-fitting method based on the Levenberg-Marquardt algorithm is proposed.This method is more suitable for small machines and overcomes some drawbacks of the traditional method like subjectivity.However, the detailed magnetic equivalent circuit parameters needed by designers can't be extracted in this method.Besides, the influence of inductance variation during dynamic processes on transient parameters is not considered and represented in the afore-mentioned methods.
For designers, a novel computational method that can produce more detailed information is needed for the computation of the transient parameters of an SPSM.Several demands of the required method are as follows.The inductance parameters and the time constants are expected to be calculated directly and automatically by FE analysis.The detailed inductance parameters of the magnetic equivalent circuit should be extracted separately for design optimization.And the inductance variation during dynamic processes ought to be considered in the computation.
Accurate computation of magnetic equivalent circuit parameters using FE analysis is an important part of inductance calculation in the novel method.A lot of related work has been previously reported in the literature [13][14][15][16][17]. Magnetic vector potential (MVP) methods have been used to calculate the inductance parameters [13].In addition, the frozen permeability method was applied in [14] to consider the core saturation and skin effect of the solid conductor bar.A detailed SPSM leakage inductance model has been generated that considers the leakage flux distribution flowing across the air gap [15][16][17].In the literature [16,17], the influence of various operational statuses on leakage inductance variation is discussed, and the leakage inductances are obtained at the steady state FE simulation.But, the magnetic field distribution of an SPSM in transient process is time varying and different from steady state as investigated in this paper.This phenomenon should also be considered in the computational analysis of leakage inductance and finally in the determination of transient parameters.FE computation of leakage inductances has been previously reported and discussed in [14][15][16][17][18][19].For time constant determination, a curve fitting technique might be helpful when applied to calculate the time constants [12,20].
Basing on the previous works mentioned above, a novel method for transient parameters computation of an SPSM in the design stage is proposed in this paper.In this method, time-stepping field-circuit FE simulation and static magnetic FE simulation are employed for accurate calculation.All elements of equivalent circuit parameters can be separately extracted by this method.In particular, the gap leakage inductances [15] are considered and obtained.The inductance variation, influenced by saturation variation during the whole transient process, are investigated and considered in the determination of transient parameters.To estimate time constants, the Prony algorithm is applied.The proposed method satisfies the afore-mentioned demand of SPSM designers.Besides, a program Energies 2017, 10, 1015 3 of 18 is developed to control the FE simulation and data processing [21].By employing this program, the transient parameters can be automatically calculated in the proposed method.
The rest of this paper is organized as follows.Section 2 discusses the modification of the d-q model of an SPSM, which considers the leakage inductances separately and presents the equivalent circuits for calculating the transient parameters.The calculation procedure using a 2-D FE analysis of a three-phase sudden short-circuit test simulation is described in Section 3; the procedure includes the method to obtain the magnetic equivalent circuit parameters separately considering the saturation variation and skin effects, and the method to compute the transient parameters.Sections 4 and 5 provide the results, discussion, method validation, and conclusion.

Machine Model
The d-q model describes an SPSM with damper cage whose stator core and rotor pole are both laminated.A cross-section schematic of the SPSM is shown in Figure 1.
Energies 2017, 10, 1015 3 of 18 Besides, a program is developed to control the FE simulation and data processing [21].By employing this program, the transient parameters can be automatically calculated in the proposed method.The rest of this paper is organized as follows.Section 2 discusses the modification of the d-q model of an SPSM, which considers the leakage inductances separately and presents the equivalent circuits for calculating the transient parameters.The calculation procedure using a 2-D FE analysis of a three-phase sudden short-circuit test simulation is described in Section 3; the procedure includes the method to obtain the magnetic equivalent circuit parameters separately considering the saturation variation and skin effects, and the method to compute the transient parameters.Sections 4 and 5 provide the results, discussion, method validation, and conclusion.

Machine Model
The d-q model describes an SPSM with damper cage whose stator core and rotor pole are both laminated.A cross-section schematic of the SPSM is shown in Figure 1.

Damper Cage Model
A damper cage with inter-pole connections is modeled by the multi-loop method in [15] as a multi-loop model having Nd and Nq windings along the d-and q-axes, respectively.However, this method is too complex to obtain the integral characteristics of the damper cage in practical applications.In this paper, the damper cage is modeled as two concentric windings, kD and kQ, along the d-and q-axes, respectively.kD and kQ consist of Nd and Nq concentrated coils, respectively.The turn numbers (TN) of each coil are adjusted to keep the magnetic motive force (MMF) of kD and kQ equal to the sum of the MMF of the multi-loop model windings along the d-and q-axes, as shown in Equations ( 1) and (2).

TN
( 2, ), TN ), TN 1 , where id and jq are the windings of the multi-loop model, and iD and jQ are the coils of the kD and kQ windings, respectively.TN is the turn number of the concentrated coil, and F is the MMF of the windings. . . .N q Q are the concentrated coils of the kD and kQ windings of one pole, respectively.f is the field winding.

Damper Cage Model
A damper cage with inter-pole connections is modeled by the multi-loop method in [15] as a multi-loop model having N d and N q windings along the dand q-axes, respectively.However, this method is too complex to obtain the integral characteristics of the damper cage in practical applications.In this paper, the damper cage is modeled as two concentric windings, kD and kQ, along the dand q-axes, respectively.kD and kQ consist of N d and N q concentrated coils, respectively.The turn numbers (TN) of each coil are adjusted to keep the magnetic motive force (MMF) of kD and kQ equal to the sum of the MMF of the multi-loop model windings along the d-and q-axes, as shown in Equations ( 1) and (2).
where id and jq are the windings of the multi-loop model, and iD and jQ are the coils of the kD and kQ windings, respectively.TN is the turn number of the concentrated coil, and F is the MMF of the windings.

d-q Model Flux Linkage Equation
The flux linkage equations of the d-q model of the SPSM is given in Equation (3).
where the symbols ψ, L, M, and I denote the magnetic flux linkage, self-inductance, mutual inductance, and current of the windings, respectively.Subscripts d, q, f, kD, and kQ denote the stator dand q-axis windings, rotor field winding, and dand q-axis damper windings

Leakage Inductance of the Model
To accurately calculate the transient reactance parameters, the leakage inductances of each equivalent winding are divided into the self-leakage, gap leakage, and end leakage inductances [13].
The self-leakage inductance consists of the slot leakage and belt leakage inductances.Part of the tooth-top leakage inductance, corresponding to the tooth-top leakage flux not flowing across the air gap, is included in the slot leakage inductance.When equipped with a damper cage, the belt leakage inductance is small and considered negligible in the calculation [7].
The gap leakage inductance represents the corresponding leakage flux linkage flowing across the air gap.The gap leakage inductance consists of the following two parts: one is the zig-zag leakage inductance, which is caused by the difference between the winding MMF distributions and the slots on the other side of the air gap [7], and the other is part of the tooth-top leakage inductance, corresponding to the tooth-top leakage flux flowing across the air gap [17].It is hard to separate the zig-zag leakage and tooth-top leakage inductance; therefore, the gap leakage inductance is calculated as a whole.
A schematic of the self-leakage and gap leakage flux paths is shown in

d-q Model Flux Linkage Equation
The flux linkage equations of the d-q model of the SPSM is given in Equation (3). , where the symbols ψ, L, M, and I denote the magnetic flux linkage, self-inductance, mutual inductance, and current of the windings, respectively.Subscripts d, q, f, kD, and kQ denote the stator d-and q-axis windings, rotor field winding, and d-and q-axis damper windings

Leakage Inductance of the Model
To accurately calculate the transient reactance parameters, the leakage inductances of each equivalent winding are divided into the self-leakage, gap leakage, and end leakage inductances [13].
The self-leakage inductance consists of the slot leakage and belt leakage inductances.Part of the tooth-top leakage inductance, corresponding to the tooth-top leakage flux not flowing across the air gap, is included in the slot leakage inductance.When equipped with a damper cage, the belt leakage inductance is small and considered negligible in the calculation [7].
The gap leakage inductance represents the corresponding leakage flux linkage flowing across the air gap.The gap leakage inductance consists of the following two parts: one is the zig-zag leakage inductance, which is caused by the difference between the winding MMF distributions and the slots on the other side of the air gap [7], and the other is part of the tooth-top leakage inductance, corresponding to the tooth-top leakage flux flowing across the air gap [17].It is hard to separate the zig-zag leakage and tooth-top leakage inductance; therefore, the gap leakage inductance is calculated as a whole.
A schematic of the self-leakage and gap leakage flux paths is shown in Figure 2. Here, Φsl and Φslg are the self-leakage and gap leakage flux paths of a conductor in a stator slot, respectively.Φfl and Φflg are the self-leakage and gap leakage flux paths of the field winding, respectively.Φlcb and Φlcbg are the self-leakage and gap leakage flux paths of a conductor bar in a damper cage slot, respectively.For the damper cage, the leakage flux of all the damper bars are considered as d-and q-axes leakage flux of the damper cage in the model.The gap leakage inductances of the dand q-axes are different, influenced by the tapered rotor pole face; therefore, the stator leakage inductances of the dand q-axes are not equal.Additionally, the different saturation levels along the dand q-axes cause differences between the dand q-axis slot leakage inductances of the stator winding.To consider the afore-mentioned effects, the self-leakage and gap leakage inductances of the dand q-axis windings are calculated separately in this paper.
The end leakage inductances of the d, q, field, and damper windings are assumed to be constant and unaffected by saturation [14][15][16][17].
Therefore, the leakage inductances of the SPSM in this paper are defined in Equation ( 4) as follows: where L sσd , L sσq , L fσ , L Dσ , and L Qσ denote the total leakage inductances of the stator dand q-axis windings, field winding, and kD and kQ damper winding.L sld , L slq , L fl , L Dl , and L Ql are the self-leakage inductances of the d, q, field, kD and kQ windings, respectively.L slgd , L slgq , L flg , L Dlg , and L Qlg are the gap leakage inductances and L se , L se , L fe , L De , and L Qe are the end leakage inductances of the d, q, field, kD and kQ windings, respectively.
The flux distribution in the transient process is complex, time-varying, and noticeably different from the steady state.The corresponding saturation condition should be considered in the leakage inductance computation.Furthermore, the variation in saturation also leads to the variation in leakage inductance during the transient process.These effects are considered in the calculation procedure in Section 4.

Equivalent Circuit
In the transient equivalent circuit, the inductances of the rotor side, such as inductances of field and damper windings, must be transferred to the stator side.All the inductances in Equation ( 4) should be transferred to reactances conventionally.Thus, the transient equivalent circuit of the dand q-axes without resistance is shown in Figure 3, where reactances have been expressed per unit (P.U.).The same footnotes are applied between inductances in Equation ( 4) and reactances in Figure 3, transferred from the inductances for simplification.The subtransient reactances, x d " and x q ", transient reactance, x d ', and steady state reactances, x d and x q , are calculated by the equivalent circuits, as described in Equation ( 5).
x d = x sσd + x f σ x ad x Dσ x ad x f σ +x f σ x Dσ +x Dσ x ad x Qσ x aq x aq +x Qσ x d = x sσd + x ad , x q = x sσq + x aq (5) where x denotes the P.U.equivalent reactance, corresponding to the inductance transferred into the X ad base P.U.expression.

Calculation Method
The modified model and methods are applied for an SPSM designed as a dynamic simulation generator, installed at the Power System Dynamic Simulation Laboratory.The main parameters of the prototype are given in Table 1.

Calculation Method
The modified model and methods are applied for an SPSM designed as a dynamic simulation generator, installed at the Power System Dynamic Simulation Laboratory.The main parameters of the prototype are given in Table 1.

Formulation of FE Application
Two FE models of a prototype are established and applied using commercial FE analysis software: a 2-D field-circuit coupled time-stepping FE model [6] and a 2-D static magnetic field FE model [12].The same 2-D FE model of the presented prototype is used in both applications.Maxwell's equation using MVP is shown in Equation ( 6): ) + ∂ ∂y where A represents MVP, µ e and σ are effective permeability and conductivity, respectively.J is current density.The boundary condition of Figure 4 is described in Equation ( 7): Energies 2017, 10, 1015 7 of 18 The 2-D field-circuit coupled time-stepping FE prototype model is used to simulate the three-phase sudden short-circuit test [6].By making a three-phase short-circuit fault at terminals of the stator windings, the test is employed (as a standard procedure) to obtain the transient and subtransient parameters in a machine experiment [9].It should be noted that the negative resistor device, usually connected in field circuits of dynamic simulation machines in series, is not considered in this paper.
The external circuits are set and linked to the 2-D FE model as shown in Figure 4 to perform the simulation of the test, as shown in Figure 5.The resistances of the stator and field windings can be accurately calculated by using an analytical solution from design data and are set in the winding  The 2-D field-circuit coupled time-stepping FE prototype model is used to simulate the three-phase sudden short-circuit test [6].By making a three-phase short-circuit fault at terminals of the stator windings, the test is employed (as a standard procedure) to obtain the transient and subtransient parameters in a machine experiment [9].It should be noted that the negative resistor device, usually connected in field circuits of dynamic simulation machines in series, is not considered in this paper.
The external circuits are set and linked to the 2-D FE model as shown in Figure 4 to perform the simulation of the test, as shown in Figure 5.The resistances of the stator and field windings can be accurately calculated by using an analytical solution from design data and are set in the winding components in the external circuit.The end leakage inductance of the stator and field windings (L se and L fe in Figure 5), damper cage end ring equivalent resistance, and end leakage inductance of the adjacent bars (R c , L c , R ic , L ic ) are also obtained analytically [12,14].The external circuit parameters are summarized in Table 2.In the three-phase sudden short-circuit test, the SPSM is driven by the prime mover, typically a motor in the Power System Dynamic Simulation Laboratory; hence, the fluctuation of the rotation speed can be neglected.The rotor rotation speed is therefore constant at the rated synchronous speed in the mechanical setting of the time-stepping FE model.The 2-D field-circuit coupled time-stepping FE prototype model is used to simulate the three-phase sudden short-circuit test [6].By making a three-phase short-circuit fault at terminals of the stator windings, the test is employed (as a standard procedure) to obtain the transient and subtransient parameters in a machine experiment [9].It should be noted that the negative resistor device, usually connected in field circuits of dynamic simulation machines in series, is not considered in this paper.
The external circuits are set and linked to the 2-D FE model as shown in Figure 4 to perform the simulation of the test, as shown in Figure 5.The resistances of the stator and field windings can be accurately calculated by using an analytical solution from design data and are set in the winding components in the external circuit.The end leakage inductance of the stator and field windings (Lse and Lfe in Figure 5), damper cage end ring equivalent resistance, and end leakage inductance of the adjacent bars (Rc, Lc, Ric, Lic) are also obtained analytically [12,14].The external circuit parameters are summarized in Table 2.In the three-phase sudden short-circuit test, the SPSM is driven by the prime mover, typically a motor in the Power System Dynamic Simulation Laboratory; hence, the fluctuation of the rotation speed can be neglected.The rotor rotation speed is therefore constant at the rated synchronous speed in the mechanical setting of the time-stepping FE model.To consider the skin effects of the damper bars, each damper bar is divided into several segments in the 2-D FE model.The current information recorded in the field-circuit FE simulation is imported into the static magnetic FE simulation to produce a magnetic flux distribution similar to the flux distribution in time-stepping FE solution for accurately computing the inductances.

Calculation Procedure
The procedure to calculate the transient d-q axis reactance parameters is described as follow.
Step 1: The three-phase sudden short-circuit test is simulated using the 2-D field-circuit coupled time-stepping FE model.The total currents of the stator windings (I A , I B , I C ), field winding (I F ) and each segment of each damper bar (I b11 , I b12 , . . ., I b18 , I b21 , I b22 , . . . ) are obtained and exported in the simulation results at each time step in the whole simulation time range.The rotor position at each time step is also recorded.All currents and rotor position data are related to the time instants.
Each damper bar is divided into eight geometric segments in the 2-D FE model to consider the skin effects mentioned in Section 3.1 As the current density of each segment and the damper bar can be obtained from MVP calculated in the FE solution, the total current of each segment in each damper bar can be obtained by using Equation (8).The current density distributions of the damper bar segments are shown in Figure 6.
where A bki is the MVP in the damper bar region, U bk is the voltage drop of the damper bar, σ is the conductivity of the conductor bar, l ef is the effective length, k represents the serial number of the damper bar and i is the serial number of the segments in a bar region.I bk represents the total current of a damper bar and I bki is the current of a segment in a damper bar.S ki is the integration region of each segment.
Energies 2017, 10, 1015 9 of 18 frozen and exported, corresponding to the time instants, as the permeability distribution can represent the iron saturation condition.Step 2: The recorded currents and rotor positions at each calculation time instant are imported into the 2-D static magnetic FE model to obtain the same magnetic field distribution as the time-stepping transient FE application in step 1.As shown in Figure 7, the air-gap flux density distributions are the same for the field-circuit coupled time-stepping FE application in step 1 and the static magnetic FE simulation in step 2. Therefore, the magnetic field distribution and iron saturation condition in both simulations are nearly uniform.It is difficult to obtain the inductance directly; hence, the inductance computed by using the iron saturation condition in the static magnetic FE application, carried out in step 2, can accurately represent the flux distribution of the field-circuit FE application obtained in step 1.The frozen permeability method is applied to record the core saturation information of the results in step 2 for inductance computation in step 3.In this method, the permeability of all the mesh elements Energies 2017, 10, 1015 9 of 18 in the stator and rotor core regions of the results in step 2 is frozen and exported, corresponding to the time instants, as the permeability distribution can represent the iron saturation condition.Step 3: After freezing the permeability, the kD or kQ model is implemented in the circuit set of damper bar components in the static FE circuit model, corresponding to the time instants.The mesh elements in the model must remain unchanged.By importing the frozen permeability of each mesh element obtained in step 2 into the same mesh elements in the same core region, the iron saturation condition in step 2 is reproduced in the FE model of step 3 by the same permeability distribution in core regions.The inductances are calculated with the method presented in Section 3. 3 The flux distributions at the no-load condition before the short-circuit in Figure 7a and at the transient condition after the short-circuit in Figure 7b are vastly different.The variation in the flux distribution can lead to variation in the inductances of the SPSM, including leakage inductances; therefore, the variation should be investigated and considered in the computation.To calculate the inductance variation, a recursive inductance calculation is applied to each chosen time instant during the entire period, from the short-circuit to the steady short circuit state-referred to as the transient process in this paper.For simplification, only 2p (p represents pole pairs) of one rotation cycle is needed to calculate the d-and q-axis inductances to represent the inductance variation.Thus, in the period, the d-axis inductances are calculated at the time instants when the d-axis meets the stator A phase winding axis; the q-axis inductances are calculated similarly.Steps 2 and 3 are carried out and repeated at each chosen time instant of the transient process.Step 3: After freezing the permeability, the kD or kQ model is implemented in the circuit set of damper bar components in the static FE circuit model, corresponding to the time instants.The mesh elements in the model must remain unchanged.By importing the frozen permeability of each mesh element obtained in step 2 into the same mesh elements in the same core region, the iron saturation condition in step 2 is reproduced in the FE model of step 3 by the same permeability distribution in core regions.The inductances are calculated with the method presented in Section 3. 3 The flux distributions at the no-load condition before the short-circuit in Figure 7a and at the transient condition after the short-circuit in Figure 7b are vastly different.The variation in the flux distribution can lead to variation in the inductances of the SPSM, including leakage inductances; therefore, the variation should be investigated and considered in the computation.To calculate the inductance variation, a recursive inductance calculation is applied to each chosen time instant during the entire period, from the short-circuit to the steady short circuit state-referred to as the transient process in this paper.For simplification, only 2p (p represents pole pairs) of one rotation cycle is needed to calculate the dand q-axis inductances to represent the inductance variation.Thus, in the period, the d-axis inductances are calculated at the time instants when the d-axis meets the stator A phase winding axis; the q-axis inductances are calculated similarly.Steps 2 and 3 are carried out and repeated at each chosen time instant of the transient process.

Inductance Calculation
In the stator self-leakage inductance calculation, the frozen permeability of the stator core is imported into the FE model [18].The rotor parts in the FE model should be removed to extract the stator self-leakage magnetic field.Technically, a boundary between the rotor region and air gap region is built in to the FE model.The boundary condition is set the same as Equation (7).Thus, the excitation of stator currents can only generate stator self-leakage field in the model.The self-leakage inductance can be obtained using the stored magnetic energy computation and circuit equations as described in [14].At the time instants when the d-axis meets the A phase winding axis, the d-axis stator self-leakage inductance is calculated, as shown in Equations ( 9)- (11).As the non-linear core regions are processed linearly by the frozen permeability method, the influence of non-linear core saturation is considered in inductance calculation.Besides, after the linear processing, both the field Equation ( 9) and circuit Equation ( 10) can be used in this case for stator self-leakage inductance computation.
where W m is the stored magnetic energy for the whole domain, S is the integration region, H is the magnetic field intensity and B is the magnetic flux density.L s ' is the self-inductance of the stator phase winding, and M s ' is the mutual inductance between two phase windings.
In step 3 of the static FE simulation, I A is set to 1 A, and I B and I C are set to −0.5 A. W m can be directly obtained from Equation ( 9) by using the FE result.As only the stator self-leakage magnetic field is generated, L s ' and M s ' are equal to the self-leakage inductance of a stator phase winding (L spσ ) and the mutual leakage inductance between two stator phase windings (M spσ ), respectively, in this simulation [14].Thus, the currents in Equation ( 10) and result of Equation ( 9) are applied in Equation ( 11) to obtain the d-axis stator self-leakage inductance.
where L sld is the d-axis stator self-leakage inductance.
The flux distribution of the stator d-axis winding's self-leakage flux is shown in Figure 8a.The same method is employed to compute the q-axis stator self-leakage inductance at the time instants when the q-axis meets the stator A phase winding axis.Additionally, the self-leakage inductances of the field and damper windings are calculated by importing the rotor frozen permeability magnetic characteristic setting in the FE model, exciting only I F , I D , or I Q .
The flux linkages of windings are obtained from MVP computed by FE solutions, as shown in Equation (12).The dand qaxis flux linkages are calculated when the d-axis and q-axis meet with the stator A phase winding axis, respectively.
where A Uk and A Lk are the MVP of the upper and lower sides of coil k in winding i, respectively, β is the total number of coils in winding i, l ef is the axial effective length of prototype.The frozen permeabilities of the stator and rotor cores are imported into the FE model to calculate the armature reaction inductances.Only I F is excited to calculate the mutual inductance between the stator d-axis winding and field winding M adf , and only I Q is excited to calculate the mutual inductance between the stator q-axis winding and kQ winding M aqQ .The armature reaction inductances are calculated as: where ψ d and ψ q are the flux linkages of the dand q-axis stator windings, respectively, i b , i fb , and i Qb are the P.U.current base values of the stator winding, field winding, and kQ damper winding, respectively.The flux distribution of only I F excitation to compute M adf is shown in Figure 8b.The gap leakage inductances are also calculated by importing both the stator and rotor permeabilities.The stator gap leakage inductance of the d-axis is calculated in Equation ( 14).The q-axis gap leakage inductances are obtained in a similar manner.The field gap leakage inductance is computed in Equation ( 15), and the gap leakage inductances of kD and kQ windings are obtained similarly.
The length of the transient process from the short-circuit to the steady short state is measured by a transient FE application.The transient process consists of a subtransient process as the period of the damper current attenuation, until the aperiodic component of the damper current decays to less than 5 percent of the initial value.The leakage inductance of the damper windings is computed for all the chosen time instants in the subtransient process; the other leakage inductances are obtained for all the chosen time instants during the whole transient process.
Finally, the subtransient and transient reactances xd", xq", and xd' for each selected time instant are calculated on a P.U.basis in Equations ( 4) and ( 5) by separately transferring all inductances to a P.U.basis.
The subtransient reactance parameters Xd" and Xq", as well as the transient reactance parameter Xd', of the SPSM are defined as the mean value of xd", xq", and xd' at all time instants from the short-circuit to the time when the transient aperiodic component of the kD, kQ, and field winding The gap leakage inductances are also calculated by importing both the stator and rotor permeabilities.The stator gap leakage inductance of the d-axis is calculated in Equation ( 14).The q-axis gap leakage inductances are obtained in a similar manner.The field gap leakage inductance is computed in Equation ( 15), and the gap leakage inductances of kD and kQ windings are obtained similarly.
The length of the transient process from the short-circuit to the steady short state is measured by a transient FE application.The transient process consists of a subtransient process as the period of the damper current attenuation, until the aperiodic component of the damper current decays to less than 5 percent of the initial value.The leakage inductance of the damper windings is computed for all the chosen time instants in the subtransient process; the other leakage inductances are obtained for all the chosen time instants during the whole transient process.
Finally, the subtransient and transient reactances x d ", x q ", and x d ' for each selected time instant are calculated on a P.U.basis in Equations ( 4) and ( 5) by separately transferring all inductances to a P.U.basis.
The subtransient reactance parameters X d " and X q ", as well as the transient reactance parameter X d ', of the SPSM are defined as the mean value of x d ", x q ", and x d ' at all time instants from the short-circuit to the time when the transient aperiodic component of the kD, kQ, and field winding currents, respectively, decay to 20% of the initial value; here, the saturation variation in the constant terminal reactance parameters of the SPSM is taken into account.

Time Constant Computation and Developed Program
The short-circuit time constants T d ", T q ", T d ', and T a are obtained by curve fitting and applying the Prony algorithm.The Prony algorithm exponentially provides the frequency, magnitude, and damping factor information for all modes of the imported discrete transient signal.As described in [19], the Prony algorithm and least squares algorithm are used to determine the short-circuit time constant.The mathematical model of the Prony algorithm is shown in Equations ( 16) and (17).
where ŷ(k) is the approximate representation of sampling current signal.n is the number of total harmonics of the signal, and A i , θ i , B i , α i , and f i are the peak value, phase, amplitude, damping factor, and frequency of each exponential harmonic component of the fitting current, respectively.The time constant of the harmonic component, obtained by the damping factor, is: In the presented method, the d-axis short-circuit time constant, including the d-axis subtransient short-circuit time constant T d " and the d-axis transient short-circuit time constant T d ', is estimated from the armature short-circuit current using the Prony algorithm; the armature short-circuit time constant T a is also estimated this way.The q-axis subtransient short-circuit time constant T q " is obtained from the kQ current waveform.

Developed Program
To address the complex simulation and data processing, a program is developed by employing PYTHON and visual basic to control the commercial FE software and data processing.The entire calculation is automatic after applying the program to the FE model.The structure of the program is shown in Figure 9.In this program, a loop process of inductances computation whose calculation methods are described in Section 3.3 is applied to compute the inductances in each chosen time instants.
where y ∧ (k) is the approximate representation of sampling current signal.n is the number of total harmonics of the signal, and Ai, θi, Bi, αi, and fi are the peak value, phase, amplitude, damping factor, and frequency of each exponential harmonic component of the fitting current, respectively.The time constant of the harmonic component, obtained by the damping factor, is: In the presented method, the d-axis short-circuit time constant, including the d-axis subtransient short-circuit time constant Td" and the d-axis transient short-circuit time constant Td', is estimated from the armature short-circuit current using the Prony algorithm; the armature short-circuit time constant Ta is also estimated this way.The q-axis subtransient short-circuit time constant Tq" is obtained from the kQ current waveform.

Developed Program
To address the complex simulation and data processing, a program is developed by employing PYTHON and visual basic to control the commercial FE software and data processing.The entire calculation is automatic after applying the program to the FE model.The structure of the program is shown in Figure 9.In this program, a loop process of inductances computation whose calculation methods are described in Section 3.3 is applied to compute the inductances in each chosen time instants.

Results and Comparison
The transient parameters of the prototype are determined based on the abovementioned model and methods.All inductances are transferred to reactances in X ad P.U.basis, and the time instant of three-phase sudden short circuit is set to 0 s in this section.
The leakage inductances of stator, field, and damper windings and their variations in the dand q-axes are shown in Figures 10 and 11, respectively.It is shown that the gap leakage inductances play a significant role in the leakage inductances, particularly for the field winding.If the gap leakage inductances are neglected, as in the traditional design and analysis of an SPSM, the leakage inductances may be inaccurately estimated.
The magnetic flux line of the gap leakage flux, shown in Figure 12, verifies the leakage inductance model.Arrow 1 illustrates the field winding gap leakage flux line flowing through the pole shoe, air gap, and stator slot opening; Arrow 2 indicates the stator winding gap leakage flux line mainly flowing across the air gap and damper cage tooth.Compared to the megawatt-class SPSM, such as a large hydro-generator [13][14][15], the gap leakage inductance is more evident in a kilowatt-class SPSM due to the differences in geometric design, including a semi-closed stator slot and a short air gap length.and methods.All inductances are transferred to reactances in Xad P.U.basis, and the time instant of three-phase sudden short circuit is set to 0 s in this section.
The leakage inductances of stator, field, and damper windings and their variations in the d-and q-axes are shown in Figures 10 and 11, respectively.It is shown that the gap leakage inductances play a significant role in the leakage inductances, particularly for the field winding.If the gap leakage inductances are neglected, as in the traditional design and analysis of an SPSM, the leakage inductances may be inaccurately estimated.and methods.All inductances are transferred to reactances in Xad P.U.basis, and the time instant of three-phase sudden short circuit is set to 0 s in this section.
The leakage inductances of stator, field, and damper windings and their variations in the d-and q-axes are shown in Figures 10 and 11, respectively.It is shown that the gap leakage inductances play a significant role in the leakage inductances, particularly for the field winding.If the gap leakage inductances are neglected, as in the traditional design and analysis of an SPSM, the leakage inductances may be inaccurately estimated.Figures 10 and 11 also show that the leakage inductances are not constant in the transient process, particularly for the damper and field windings.The fluctuation is caused by variation of core saturation during the transient process.However, the variations in the self-leakage inductance of the stator windings are small, which indicates that the saturation variation mainly occurs in the pole shoe and pole tip, specifically between the damper cage slot and shoe edge, as shown in Figure 12.The flow through the damper cage is avoided by induced currents for conservation of flux linkage; therefore, the main flux flows through the path between the damper slot and shoe edge, leading to saturation in that area.As seen in Figure 12, there is increased saturation in the preliminary stage of the transient process after the short-circuit.In addition, the transient currents decay to the steady short-circuit state and thereby alleviate the saturation, as illustrated in the  the stator windings are small, which indicates that the saturation variation mainly occurs in the pole shoe and pole tip, specifically between the damper cage slot and shoe edge, as shown in Figure 12.The flow through the damper cage is avoided by induced currents for conservation of flux linkage; therefore, the main flux flows through the path between the damper slot and shoe edge, leading to saturation in that area.As seen in Figure 12, there is increased saturation in the preliminary stage of the transient process after the short-circuit.In addition, the transient currents decay to the steady short-circuit state and thereby alleviate the saturation, as illustrated in the growing leakage inductances in Figures 10 and 11.Therefore, the traditional theory that ignores the saturation variation in the leakage inductance computation and dynamic analysis is not impeccable.The variation of the leakage inductances should be considered to acquire precise dynamic performance.In an effort to eliminate the influence of saturation in transient processes, designers should optimize geometry of the pole shoe and damper cage slot in the design stage.
It is difficult to precisely consider the skin effects, the saturation, and magnetic distortion in the analytical solution when attempting to estimate generator performance.It is also difficult to determine the detailed characteristics of the machine from terminal data.Thus, the FE analysis should be applied for SPSM design.
The transient aperiodic component of the kD, kQ, and field winding currents can be easily obtained by the time-stepping FE simulation and data processing.The time lengths of the transient current aperiodic component attenuation from the short circuit to 20% of the initial value of the kD, kQ, and field windings, considering the calculation time instants of each winding, are approximately 0.08 s, 0.125 s, and 0.24 s, respectively.The subtransient and transient reactances x d ", x q ", and x d ' at each calculation time instant are shown in Table 3.The mean values of the subtransient and transient reactances in each computation time instants, shown in Table 3, are the subtransient and transient reactance parameters of the prototype computed by the presented method.A comparison of the prototype experimental values and the transient parameters calculated using both the proposed and conventional short-circuit methods [3] is shown in Table 4.The experimental value of the prototype is provided by the Power System Dynamic Simulation Laboratory where the presented prototype is installed and applied for dynamic simulation experiments.The experimental methods recommended in standard [11] are applied to obtain the experimental transient parameters of the prototype, including transient and steady state reactances.
The comparison demonstrates the accuracy of the proposed method for computing the transient reactance parameters.
The waveforms of the stator A phase winding current and kQ winding current, obtained by the time-stepping FE simulation, are imported into a Prony algorithm-based data processing program to compute the time constant parameters, as described in Section 3.4.The comparison of the FE simulated current and Prony algorithm fitting current are illustrated in Figure 13, and the main modal information of the fitting current is shown in Table 5.There is little error between the simulated current and fitted current waveforms, as illustrated in Figure 13.Hence, the modal information obtained from the fitting current waveform can be  There is little error between the simulated current and fitted current waveforms, as illustrated in Figure 13.Hence, the modal information obtained from the fitting current waveform can be applied to precisely estimate time constants of each component's current obtained from the FE simulation.The components of a short-circuit current have been previously described [12].In the modal information of the stator short-circuit current, mode 1 is the aperiodic component of the short-circuit current, and mode 2/3 and mode 4/5 represent the d-axis transient and subtransient AC components, respectively.Mode 7/8 and mode 9/10 denote the d-axis steady state AC component and second harmonic frequency AC component, respectively.In the modal information of the kQ short-circuit current, mode 1 is the DC component, while mode 2/3 represents the rotating frequency AC component.Each time constant can be obtained from the damping factor of the corresponding mode.The time constants and their comparison with the conventional method [11] and experimental results are shown in Table 6.In Table 6, T d ", and T d ' are the d-axis subtransient and transient short-circuit time constant respectively, T q " is the q-axis subtransient short-circuit time constant and T a is the armature short-circuit time constant.The comparison validates that the new method for computing the time constants of transient parameters is effective and accurate.In the simulation and computation of this work, a personal computer (PC) with Intel ® Xeon CPU E3-1230 V2 @ 3.3 GHz and 16 GB random-access memory (RAM) is used.And the whole computation time interval applying the developed program is about 5 h and 30 min.All the presented results can be automatically calculated by the developed program.
Generally, as the aforementioned results and discussions shown, the results reveal that the gap leakage inductances and inductance variation during transient processes have significant influences on transient parameters and should be considered in computation to obtain more precise transient parameters.As shown in Figures 10 and 11, the proposed method could extract the detailed leakage inductances separately in each computation time instant.Thus, a clear understanding of leakage inductances saturation and variation during a transient process could be obtained.Some phenomena of core saturation in transient process are investigated.Detailed information of leakage inductances can be provided to designers by this method for geometric optimization.The accuracy of the proposed method in transient parameters computation is validated by the comparison with experimental and conventional method results.The proposed method has turned out to be an effective method for the transient parameter computation and transient performance optimization of an SPSM.

Conclusions
This paper presents a novel method for calculating the transient parameters of an SPSM with a damper cage.A modified model is proposed to accurately calculate the leakage inductances, and FE analysis is used to consider the skin effects and saturation.By using this method, the leakage inductances of the magnetic equivalent circuit are separately obtained, and the inductance variation during the transient process is also investigated and discussed.The Prony method is applied to simulate the current waveform to estimate the time constant parameters.The complex FE analysis and data processing can be automatically performed by the developed program.
In the proposed method, the gap leakage inductances and inductances variation during transient process are considered, which makes the computation results precise.The detailed inductances in each time instant during a transient process can be extracted separately, which could help designers for optimization and refinement.Besides, as the transient inductances are calculated from detailed inductances and time constants are computed by curve fitting of the current waveform, parameter determination becomes objective and automatic.
The accuracy of the present method is validated by comparing the results with the conventional method.

Figure 1 .
Figure 1.Schematic of the salient-pole synchronous machine (SPSM) model.1D…NdD and 1Q…NqQ are the concentrated coils of the kD and kQ windings of one pole, respectively.f is the field winding.

Figure 1 .
Figure 1.Schematic of the salient-pole synchronous machine (SPSM) model.1D . . .N d D and 1Q . . .N q Q are the concentrated coils of the kD and kQ windings of one pole, respectively.f is the field winding.

Figure 2 .
Here, Φ sl and Φ slg are the self-leakage and gap leakage flux paths of a conductor in a stator slot, respectively.Φ fl and Φ flg are the self-leakage and gap leakage flux paths of the field winding, respectively.Φ lcb and Φ lcbg are the self-leakage and gap leakage flux paths of a conductor bar in a damper cage slot, respectively.For the damper cage, the leakage flux of all the damper bars are considered as dand q-axes leakage flux of the damper cage in the model.Energies 2017, 10, 1015 4 of 18

Figure 4 .
Figure 4. 2-D FE model cross section of prototype.

Figure 4 .
Figure 4. 2-D FE model cross section of prototype.

Figure 4 .
Figure 4. 2-D FE model cross section of prototype.

Figure 5 .
Figure 5. External circuit of field-circuit coupled FE application model.(a) Stator winding circuit; (b) Rotor field winding circuit; (c) Rotor damper winding circuit of one pole and inter-pole connection.

Figure 5 .
Figure 5. External circuit of field-circuit coupled FE application model.(a) Stator winding circuit; (b) Rotor field winding circuit; (c) Rotor damper winding circuit of one pole and inter-pole connection.

Figure 6 .Figure 7 .
Figure 6.Damper bar segments and current distribution in the damper bar at 0.01 s after the short circuit by time-stepping FE simulation.

Figure 6 .
Figure 6.Damper bar segments and current distribution in the damper bar at 0.01 s after the short circuit by time-stepping FE simulation.

Figure 6 .Figure 7 .
Figure 6.Damper bar segments and current distribution in the damper bar at 0.01 s after the short circuit by time-stepping FE simulation.

Figure 7 .
Figure 7. Air-gap flux density comparison.(a) No-load condition (b) at 0.01 s after short-circuit.B nt the result of the field-circuit coupled time-stepping FE simulation in step 1, and B ns is the result of the static magnetic FE application in step 2.

Figure 8 .
Figure 8. Flux distribution of the SPSM.(a) d-axis self-leakage flux line, (b) only field winding excitation flux line.

Figure 8 .
Figure 8. Flux distribution of the SPSM.(a) d-axis self-leakage flux line, (b) only field winding excitation flux line.

Figure 9 .
Figure 9. Structure of the developed calculation program.Figure 9. Structure of the developed calculation program.

Figure 9 .
Figure 9. Structure of the developed calculation program.Figure 9. Structure of the developed calculation program.

Figure 10 .
Figure 10.Variation in the d-axis leakage inductances during the transient process.

Figure 11 .
Figure 11.Variation in the q-axis leakage inductances during the transient process.The magnetic flux line of the gap leakage flux, shown in Figure 12, verifies the leakage inductance model.Arrow 1 illustrates the field winding gap leakage flux line flowing through the pole shoe, air gap, and stator slot opening; Arrow 2 indicates the stator winding gap leakage flux line mainly flowing across the air gap and damper cage tooth.Compared to the megawatt-class SPSM, such as a large hydro-generator [13-15], the gap leakage inductance is more evident in a kilowatt-class SPSM due to the differences in geometric design, including a semi-closed stator slot and a short air gap length.

Figure 10 .
Figure 10.Variation in the d-axis leakage inductances during the transient process.

Figure 10 .
Figure 10.Variation in the d-axis leakage inductances during the transient process.

Figure 11 .
Figure 11.Variation in the q-axis leakage inductances during the transient process.The magnetic flux line of the gap leakage flux, shown in Figure 12, verifies the leakage inductance model.Arrow 1 illustrates the field winding gap leakage flux line flowing through the pole shoe, air gap, and stator slot opening; Arrow 2 indicates the stator winding gap leakage flux line mainly flowing across the air gap and damper cage tooth.Compared to the megawatt-class SPSM, such as a large hydro-generator [13-15], the gap leakage inductance is more evident in a kilowatt-class SPSM due to the differences in geometric design, including a semi-closed stator slot and a short air gap length.

Figure 12 .
Figure 12.Flux distribution and magnetic flux line near the air gap at 0.01 s after the short circuit. 1 is the gap leakage flux line of the field winding, and 2 is the gap leakage flux line of the stator winding.

Figure 12 .
Figure 12.Flux distribution and magnetic flux line near the air gap at 0.01 s after the short circuit. 1 is the gap leakage flux line of the field winding, and 2 is the gap leakage flux line of the stator winding.

Figures 10
Figures 10 and 11 also show that the leakage inductances are not constant in the transient process, particularly for the damper and field windings.The fluctuation is caused by variation of core saturation during the transient process.However, the variations in the self-leakage inductance of

Figure 13 .
Figure 13.Comparison of current waveforms from time-stepping FE simulation and Prony fitting.(a) stator A phase winding current; (b) kQ winding current.

Figure 13 .
Figure 13.Comparison of current waveforms from time-stepping FE simulation and Prony fitting.(a) Stator A phase winding current; (b) kQ winding current.
ic 0.12 µHThe 2-D static magnetic field FE model has the same geometry as the field-circuit coupled time-stepping FE model and is applied to calculate the inductances.

Table 3 .
Subtransient and transient reactance parameters values at calculation time instants.

Table 5 .
Modal identification of fitted current.

Table 5 .
Modal identification of fitted current.

Table 6 .
Comparison of transient time constant parameters.