# A Synchronous Machine Transient Model Based upon an Algebraic Loop Accounting for Nonlinearity and Cross-Magnetization

Authors: R. Felicetti; U. Lundin
Journal: International Journal of Rotating Machinery
 (2023)

---

## Abstract

The purpose of this paper is to carry out an alternative to the present transient models for field wound synchronous machines, which is able to take into account the nonlinearity of the magnetic materials as well as the cross-magnetization. After presenting the principal model structures according to the state variables, a model based on two lookup tables for the magnetizing flux linkages is introduced and built step by step. The resulting signal flowchart shows an algebraic loop within the model, where the main flux linkage rapidly converges to its instantaneous value by simple iteration. The proof of this convergence is given for both saturated and unsaturated machine. Even though the proposed model uses the total linkage flux as state variable, as many alternative models do, it does not require the inversion of the current to flux linkage function (i.e., of lookup tables). This can spare a heavy computational task, especially with very large lookup tables. In the proposed model, the computational effort in the worst case scenario is reduced to few iterations (<10). Finally, the nonlinear behavior of the model is verified in four different transient scenarios by comparing its outcomes with those of a linear model for the same test machine.


---

## Body

## 1. Introduction

The study of synchronous machines (SM) in nonstationary behavior has taken advantage of the Park theory since many decades. The separation of the combined actions of all machine magnetomotive forces (MMFs) along the two electrically orthogonal rotor magnetic axes makes the corresponding magnetic flux linkages independent from the relative position between rotor and stator. This simplifies in turn the solution of the electromechanical differential equations (DEs) which describe the machine behavior, since the machine self-inductance and mutual inductance become time invariant. Nevertheless, the nonlinearity of the magnetic properties prevents the problem to be simply solved by the theory of linear DE or the Laplace transform. Therefore, the numerical integration of DE by means of computational resources has become the main way to simulate SM-related transient phenomena making use of ad-q axis model. Transient related to the electric system dynamic stability [1, 2], symmetrical and asymmetrical faults [3, 4], load rejection [5], and asynchronous behaviors [6–8] have done large use of nonlinear SM models based on the d-q axis theory. In more recent time, the interest for the transient model of SMs has been awaken by the applications of embedded systems (ES) to the control of SM for industrial and traction applications [9–13]. In that case, the model is used in the observer of a driven motor, in order to predict its state variables. In this specific task, the model must fulfil not only the usual requirement to be accurate in the predictions but also to deliver them quickly. This fact well represents the observer dilemma between the complexity of the model on one side and the need for its computational lightness on the other side. For example, taking into account the nonlinearity and the cross-magnetization introduced by the ferromagnetic materials in a SM model [14, 15], it improves the outcome accuracy of the model itself, but it requires more computational resources. In order to tackle this computational problem, lookup tables (LUTs) [14, 16] are increasingly preferred to analytic formulas [17, 18] for implementing the nonlinearity and the cross-magnetization in SM models. In fact, LUTs can provide the armature flux linkage for a given magnetomotive force (MMF) just by entering a multidimensional numerical table and picking up its value. All that without getting the CPU involved in any sort of time-consuming calculation. Nonetheless, for the experimental and the computational determinations of the relationship between the flux leakages and the currents, see the currents as sources and the linked fluxes as outcomes [16, 17]. For this reason, the flux-linkage-to-current LUT, which is needed in the model using the flux linkage as state variable, cannot be measured or calculated straightforwardly. It must rather be obtained through the inversion of the current-to-flux-linkage function Ψi. In an ingenious closed loop [14, 16], control strategies have been invented, which are nested in a numeric iterative procedure, in order to obtain the flux-linkage-to-current LUT moving from an available current-to-flux-linkage one. Besides, the recourse to a SM model having the current as variable of state, which makes knowingly use of the straight functionΨi, does not lighten the calculation burden anyway, since a marginal inductance matrix li of the machine and its inversion are also needed [17]. This paper offers a solution for using a SM machine transient model having all machine flux linkages as state variables, without any need for the inversion of the current-to-flux-linkage LUTs. The key idea for bypassing the issue of the Ψi inversion is the recourse to LUTs which are related to the magnetizing flux linkages along the d-axis and the q-axis rather than to the whole flux linkages. Then, the machine currents can be alternatively obtained from the distance vector between the flux linkages and the main fluxes via a constant leakage matrix, which is related to the machine stray inductances. Therefore, particular attention is paid in this paper to the construction of the leakage matrix starting from the most suitable and general SM model. Moreover, provided an input for the machine state of magnetization together with the external electrical and mechanical constraints at a given time, it is proved that the proposed model converges necessarily to the instantaneous values of the state variables. In the part dedicated to the method, after having revised the possible structures for a SM model and introduced a complete SM d-q-axis model, the paper shows the linear proportionality between the machine currents and the distance vector between main and flux linkages. Hence, the alternative flux linkage-based model is built, and its convergence to the timely state variables is demonstrated. The model is tested on a 60 kVA, 400 V, 50 Hz salient pole synchronous machine (SPSM). The method for obtaining the related LUTs by means of FEM calculations is described in the appendix, whereas the implications of the nonlinearity and the cross-magnetization on the machine behavior are discussed in the paper body. The section dedicated to the results shows the simulations obtained by means of the novel model. Four different transitory behaviors have been simulated: (1) rapid transition from nominal excitation to underexcitation and back; (2) 10% sudden increase of the grid voltage; (3) rapid drop of the driving torque (-75%) under constant power factor; and (4) severe voltage sag (-60%) lasting several cycles. The results achieved through the proposed model are compared with the ones achieved for the same SPSM by a linear model (i.e., having constant d-axis and q-axis magnetizing inductances) in order to highlight the errors introduced in the results by neglecting nonlinearity- and cross-magnetization-related phenomena in the machine model.


## 2. Materials and Methods

### 2.1. A Novel Structure for the SM Model

In [19], it has been proved that by assuming the SM flux linkages or alternatively the currents as state variables, or even considering a combination of them, several different SM equivalent models can be obtained. Nevertheless, in the practical applications, the most recurrent models [17] are those represented in Figure 1. In particular, the model in Figure 1(a) assumes the armature flux linkage vector Ψ as a state variable whereas the model in Figure 1(b) uses the armature currents vector i to that end. The last model is straightforwardly apt to host a LUT that relates the flux linkage vector to the armature currents. In fact, as soon as the machine is considered anhysteretic, the function Ψi becomes single-valued. Hence, it can be univocally approximated, in a domain D of the armature currents, by interpolating a finite number of its samples collected in a LUT. Since the anhysteretic assumption has the collateral drawback of disregarding part of the magnetic losses in the machine, they need to be reintroduced into the model later on. Fortunately, there are many available methods [9] for doing that, and this specific topic will therefore not be dealt in the present paper.Figure 1
Typical structures of a SM model having (a) the flux linkage or (b) the current as a variable of state.
(a)(b)However, choosing the current vector as a state variable procures additional complexity to the model: the construction of the incremental inductance matrix(1)l=∂Ψd∂id∂Ψd∂iq∂Ψq∂id∂Ψq∂iqon one side and the calculation of its inverse l−1 on the other side. Both [16, 17] agree in preferring the model of Figure 1(a) to that of Figure 1(b) because it leads to more simple state equations. Moreover, in ES making use of a full order observer, the flux linkage vector is always available. Conversely, if the observer does not have full rank, the model of Figure 2(b) must be forcibly used according to the state of the art, due to the lack of estimations for the flux linkages.Figure 2
The proposed model structure.The novel model proposed by the present work is represented in Figure2.The following can be observed: (I) it keeps all flux linkages, not only those of the armature, as state variables; (II) it implements the functioniΨ of Figure 1(a) by combining a constant matrix Lσ−1, to be yet defined, and a nonlinear single-valued function Ψmi, which represents the armature main flux linkages. Feature (I) spares the use of (1) in the model as well as the issue of its inversion, whereas feature (II) introduces the novelty of the suggested model. In particular, it does not require the construction of the inverse function iΨ starting from Ψi as shown in [14, 16], and it is applicable to reduced order observers as well, since its saturation model Ψmi works upon available current measurement/estimation rather than on flux linkage ones.However, the exploitation of this model requires two aspects to be clarified in advance. First of all, the machine currents are a linear combination of flux linkages and magnetizing flux linkages according to(2)i=Lσ−1Ψ−CTΨm,with Lσ−1 and C being the constant matrices to be defined. Second, the block diagram within the shaded box in Figure 2 forms an algebraic loop, which sets a convergence problem in a discrete time integration of the related DEs. At the generic integration instant tk, the highlighted block gets the estimated vector of the flux linkages Ψk as an input. Two possible outcomes can occur within the next time instant tk+1: either the generic main fluxes Ψm,j converge to their actual values at tk—providing through that an estimation for the currents ik at the same time—or they diverge making the simulation progress impossible. From a system perspective, the proof of the convergence within the loop is equivalent to demonstrate the existence of a solution i in
(3)i=Lσ−1Ψ−CTLUTi,for all possible Ψ, once the main flux Ψmi is expressed by means of LUTs
(4)Ψmi=LUTi,where
(5)LUTi=LUTdiLUTqi.


### 2.2. Stray Flux Linkages vs. Current Proportionality

A quite general Park model for a wound field salient pole synchronous machine (SPSM) is represented in Figure3, where two armature circuits (outside the shaded circle) and three rotor circuits (inside the shaded circle) are considered. In particular, with reference to the rotor, there are the excitation circuit (f) and a single damper circuit (D) on the d-axis and a single damper circuit (Q) on the q-axis.Figure 3
The coordinated rotor and stator circuits for a wound field SM having single damper circuits on thed- and q-axes, respectively.A more precise consideration of a round rotor SM (RRSM) would have required two more rotor circuits on both axes indeed: one for taking into account the presence of the massive and conductive rotor body and the other one for the eventual conductive slot wedges. In the same way, the consideration of a permanent magnet excited SM (PMSM) should replace the excitation winding in Figure3 with a suitable and permanent source of MMF, in order to fit in the representation. Nevertheless, the circuits of Figure 3 are general enough for presenting the proposed model in its entirety.The more simplistic assumption often done about the Park model of Figure3 is that the circuits along the d-axis share all the same magnetizing flux. In other words, their mutual inductances referred to the armature side are considered all equal to the main armature inductance on the d-axis, Lmd(6)LfD′=Lfa′=LDa′=Lmd.Canay [20], Hiramatsu et al. [21], and Guorui et al. [22] have proved that (6) does not always hold. Not only the main fluxes but also the stray fluxes contribute to the mutual linkage between the circuits laying on the same magnetic axis. Figure 4(a) shows, e.g., that the magnetic coupling between the excitation winding and the D-damper bars in a RRSM takes place through the stray fluxes too. In fact, all rotor circuits, excitation winding and damper bars, are usually laid down in the same slots.Figure 4
Different flux linkages between threed-axis circuits in a wound field synchronous machine having (a) round rotor or (b) salient pole rotor.In this way, the main fluxΦmDf, shared between excitation winding and D-damper bars, becomes larger than Φmd, i.e., the flux they have in common with the armature. In a SPSM having grid-shaped damper bars on each pole, a part Φσaf of the magnetic flux Φmaf—produced by the armature and linked to the field winding—does not link the D-damper bars. Therefore, Figure 4(b) shows that the flux Φmd linking the D-damper circuit can be smaller than the flux Φmaf shared between the armature and the excitation winding. The practical consequence of these remarks is the introduction of an additional stray inductance LfD′−Lmd in the usual Park model as shown in Figure 5. It takes into account the algebraic difference between the total main flux, which links the armature on the d-axis, and the total flux shared between field winding and D-damper circuit.Figure 5
Direct axis equivalent inductance model referred to the armature.According to the observations presented above, it results thatLfD′−Lmd can be negative for SPSM while it is always positive for RRSM.Since the circuits along theq-axis in Figure 3 are only two, there is just one quadrature main flux to be considered. Therefore, the equivalent magnetic circuit on the q-axis looks like the one represented in Figure 6.Figure 6
Quadrature axis equivalent inductance model referred to the armature.Having the circuits of Figures5 and 6 under consideration at once, the following relationship can be set between the currents and the flux linkages at play.
(7)idiqif′iD′iQ′=Lσ−1Ψd−ΨmdΨq−ΨmqΨf′−ΨmdΨD′−ΨmdΨQ′−Ψmq,all quantities used in (7) being referred to the armature side. Hence, the inverse of the 5×5 matrix Lσ, called the matrix of the stray inductance, takes the form
(8)Lσ−1=1Lσa000001Lσa00000LσD′+LmD′−LmdK−LfD′−LmdK000−LfD′−LmdKLσf′+Lmf′−LmdK000001LσQ′,where
(9)K=Lσf′+Lmf′−LmdLσD′+LmD′−Lmd−LfD′−Lmd2.The Blondel coefficient was introduced:(i)
Stray coefficient for the armature(10)σa=1−LmdLd.Laible’s coefficients are as follows [23]:
(i)
Stray coefficient between armature and field winding(11)σaf=1−Laf′2LdLf′(ii)
Stray coefficient between armature andD-damper circuit
(12)σaD=1−LaD′2LdLD′(iii)
Stray coefficient between armature and field winding(13)σfD=1−LfD′2Lf′LD′(iv)
Stray coefficient between armature andQ-damper circuit
(14)σaQ=1−LaQ′2LqLQ′(v)
Armature shielding coefficient(15)μa=1−Laf′LaD′LfD′Ld(vi)
Field winding shielding coefficient(16)μf=1−LfD′Laf′LaD′Lf′(vii)
D-damper circuit shielding coefficient
(17)μD=1−LfD′LaD′Laf′LD′,as well as the machine anisotropy factor at nominal conditions
(18)AFN=LdNLqN−1.The inverse matrix of the stray inductances ((8)) turns into
(19)Lσ−1=1LdN1σa000001σa00000μD/1−σaD+σa+μa−σaμa/1−μaσK2−σa+μa−σaμa/1−μaσK2000−σa+μa−σaμa/1−μaσK2μf/1−σaf+σa+μa−σaμa/1−μaσK2000001+AFNσQ.In (19), LdN represents the armature direct axis inductance at nominal conditions and
(20)σK2=KLd2=μD1−σaD+σa+μa−σaμa1−μaμf1−σaf+σa+μa−σaμa1−μa−σa+μa−σaμa1−μa2.When studying (10)–(17) closer, it can be recognized that the stray and the shielding coefficients can be regarded as constants. In fact, they represent ratios where numerators and denominators vary according to the same underlying nonlinear magnetization characteristics, irrespective to d-axis or q-axis consideration. Adding that LdN and AFN are also constants, it can be concluded that (19) is a constant matrix and that (7) sets a linear proportionality between the vector of the stray fluxes and the vector of the machine currents.


### 2.3. Variable Convergence within the Algebraic Loop

The armature main flux linkages on thed- and q-axes depend on the superposed action of all MMFs at work on the d- and q-axes, respectively. Once the currents of the rotor circuits have been referred to the armature side, they are armature currents at all effects, so that the sources of the main fluxes can be summarized as follows:
(21)iΘdiΘq=∑id′∑iq′=Ci,with
(22)C=1011001001.A more concise expression for (7) can be written by using (22).
(23)i=Lσ−1Ψ−Lσ−1CTΨm,which, once substituted in (21), takes the form
(24)iΘdiΘq=CLσ−1Ψ−CLσ−1CTΨm.For later considerations, it is here anticipated that(25)CLσ−1CT=1Lσaα−100γ−1,where
(26)α=KK+LσaLσD′+LmD′+Lσf′+Lmf′−2Lmd<1,(27)γ=LσQ′LσQ′+Lσa<1.When considered in relative terms, (26) and (27) become
(28)α=σK2σK2+σaμD/1−σaD+μf/1−σaf+2σa+μa−σaμa/1−μa<1,(29)γ=σQ1−1+AFNσa2σQ1−1+AFNσa2+σa1+AFN1−σaQ<1.By considering the discrete values of (5) as one possible choice among the ∞2 sets of values generated by the magnetizing flux functions ΨmdiΘd,iΘq and ΨmqiΘd,iΘq on iΘd,iΘq∈D, it is possible to generalize (5) as
(30)LUTi=ΨmdiΘd,iΘqΨmqiΘd,iΘq,so that the inverse of (30)
(31)LUT−1Ψmd,Ψmq=iΘdΨmd,ΨmqiΘqΨmd,Ψmqexists locally in D−1 as soon as
(32)det∂Ψmd∂iΘd∂Ψmq∂iΘd∂Ψmd∂iΘq∂Ψmq∂iΘq=detlm≠0.By considering (21), (24), and (31) at the same time, it is possible to obtain
(33)Ψm=−CLσ−1CT−1LUT−1Ψm−CLσ−1Ψ.Equation (33) sets a Banach fixed point problem, where the numerical sequence
(34)Ψmn+1=ϕΨmnconverges to a finite value
(35)Ψm=limn⟶∞Ψmn,if and only if, ϕ represents a contracting map of D−1 towards D or simply when
(36)M=supΨm∈Dϕ′Ψm<1.The functionϕ′ stands for the Jacobian of ϕ. Since Ψ is kept constant during the eventual convergence in (33), the Jacobian of ϕ has the following expression:
(37)ϕ′Ψm=−CLσ−1CT−1∂LUT−1∂Ψm.Due to (32) the function Ψm can be inverted, so that (37) turns into
(38)ϕ′Ψm=−CLσ−1CT−1∂LUT∂iΘ−1.Finally, taking into account (25), (30), and (38), the inequality (36) changes into
(39)supΨm∈D−αγLσa2lm,ddlm,qq−lm,dq2<1,where lm,dd, lm,dq, and lm,qq are the elements of the symmetric matrix lm accounting for the incremental magnetizing inductances of the SM.By making use of (39), it is possible to give a proof of the studied convergence (35) within the algebraic loop of the proposed SM model, for both the nonsaturated and saturated machines. In fact, when the SM is not saturated, the incremental inductances (which stand for the slopes of the tangents to the magnetizing curves on the d- and q-axes, respectively) are comparable with the magnetizing inductances (which stand for the slopes of the chords from the origin to the magnetizing curves along the d- and q-axes, respectively). In other terms,
(40)lm,ddlm,qq≅11+AFNLdN2.Moreover, the lack of saturation, i.e., the availability of plenty of nonoriented Weiss domains in the SM magnetic circuit, makes the competition for the magnetic polarization between the coordinated MMFs practically inexistent. That means(41)lm,dq≅0.Taking into account (10) and substituting (40) and (41) in (39) lead to the inequality
(42)Mnosat=−αγσa21+AFN<1.Expression (42) fulfils condition (36) for the convergence. In fact, even assuming an unrealistically large armature stray inductance of about 30%, the d-axis inductance should be at least 10 times larger than the q-axis one in order to have σa21+AFN>1.If the SM is saturated instead, any further differential increase of the coordinated MMFs will not gain any new Weiss domain to their respective directions, so that the incremental magnetizations contributed by the MMFs can be regarded as essentially performed in air. For the same reason, the magnetic anisotropy of the machine fades away makingAF≅1. Therefore,
(43)lm,ddlm,qq≅σa2LdN2.At the same time, the lack of available domains to compete for inhibits the mutual coupling between the coordinated circuits through the cross-magnetization. This results in(44)lm,dq≅0.By substituting (43) and (44) in (39), the fulfilment of (32) for the case of the machine saturation is also achieved, since
(45)Msat=−αγ<1.Due to the local monotonic nature of the magnetization functions (30), it must be expected the transition from (42) to (45) to be also monotonic, to ensure the convergence (35) for every simulated machine behavior between saturation and nonsaturation.


### 2.4. Model Implementation

The model presented in Figure2 has been implemented by means of three functional blocks shown in Figure 7(a). The block performing the integration is the same used in the state of the art and presented in Figures 1(a) and 1(b). It has essentially the goal to gain the state variable Ψ by integrating the unbalances between the impressed voltages, the resistive voltage drops, and the transformational electromotive forces, where present. The nonlinearity manager (NLM) block, in Figure 7(b), receives the current vector as an input and refers the rotor currents to the armature side by making use of the following transformation ratios:
(46)kf=LmfLaf,kD=LmDLaD,kQ=LmQLaQ.Figure 7
Model implementation: (a) fundamental blocks; (b) nonlinearity manager; (c) flux-to-current converter.The transformed currents, which provide their MMFs on the samed- or q-axes, are summed up in order to determine the armature equivalent magnetizing currents iΘd and iΘq, respectively. The LUTd reproduces the function ΨmdiΘd,iΘq shown in Figure 8(a), whereas LUTq reproduces the ΨmqiΘd,iΘq shown in Figure 8(b). One possible way to build the lookup tables for a given SM by using a FEM analysis program is duly described in the appendix to the present work.Figure 8
Functions of the magnetizing flux linkages according to (a) thed-axis and (b) the q-axis.
(a)(b)As soon as the coordinated magnetizing current is introduced in the LUTs, two values for the armature main flux linkages are returned as outputs,Ψmd and Ψmq, respectively.These quantities provide the inputs for the flux-to-current (FTC) block in Figure7(c) together with the vector of the flux linkages Ψ. The rotor-related flux linkages must be referred to the armature side by means of the ratios (46). Then, as already shown in (7), the block performs a linear combination of the linked stray fluxes through the matrix Lσ−1, in order to get an estimate of the machine currents. The currents obtained by that way are all referred to the armature side. Therefore, the rotor-related currents need to be referred back to the rotor by using the usual ratios (46).


### 2.5. Performed Test Simulations

The proposed nonlinear model for SM has been tested on the wound field SPSM (WFSPSM) shown in Figure9, the parameters of which are given in Table 1.Figure 9
A cross-section of the test SM with magnetic field lines produced by (a) the armature direct current and (b) the armature quadrature current.
(a)(b)Table 1
Main parameters of the test SM.
SymbolDescriptionValueANRated power60 kVAVNRated voltage400 VnNRated speed1500 rpmfNRated frequency50 HzcosφNRated power factor0.9δmin⋯δmaxAirgap range at the pole shoe3.2÷5.7 mmαpPole coverage0.76liActive axial length184 mmrδAirgap average radius161 mmqSlots per pole per phase6nsConductors per slot8rArmature phase resistance (p.u.)6.2·10−2rfField winding resistance (p.u.)5.8·10−3rDd-axis damp. circuit resistance (p.u.)7.7·10−3rQq-axis damp. circuit resistance (p.u.)15.3·10−3xdNNominald-axis reactance (p.u.)1.32xqNNominalq-axis reactance (p.u.)0.72AFNNominal anisotropy factor0.83σaArmature stray factor5.2%σafField winding to armature stray factor9.7%σaDArmature tod-axis damper circuit stray factor11.8%σQq-axis damper circuit stray factor13.3%μaArmature shielding coefficient1.7%μfField winding shielding coefficient9.5%μDq-axis damper circuit shielding coefficient11.7%In order to prove the capability of the suggested model in handling the nonlinearity and the cross-magnetization in the test machine, the model itself has been compared with a linear and anisotropic one, where the LUTs of Figure7(b) have been replaced by constant magnetizing inductances Lmd and Lmq. The latter has been chosen, for each led test, so as to match the initial conditions of the nonlinear model, before an external cause intervenes to push the machine away from its initial steady state. In this way, the two models are aligned at the beginning of the test, and all subsequent divergences between them can be recognized more clearly. In particular, four tests have been envisaged, that, for the nature of the chosen perturbation, end up changing the magnetizing flux and highlighting the magnetic nonlinearity of the machine.The first test foresees the sudden drop of the excitation voltage from its nominal value down to 40% of it, while the nominal driving torque of the generator is kept constant. Therefore, the generated nominal active power remains unchanged. Due to the lack of exciting MMF, the machine must require the magnetizing current from the grid, moving by that from the generation to the consumption of reactive power.The second test is based on the observation that the voltage level at the generator bus bar is responsible for the incoming/outgoing reactive power. Therefore, a grid voltage increase, when the excitation current is kept constant, results in a higher magnetizing flux in the SM. This should in turn shift the generator towards the saturation and possibly highlight its underlying nonlinear behavior.The third test reproduces the adjustment of the excitation current when the generated power is suddenly reduced down to 25% of its nominal value. In order to keep the power factor constant, the excitation current is decreased accordingly, so that the exciting MMF is drastically reduced.The fourth and final simulations deal with a low-voltage fault ride through (LVFRT), where the grid voltage experiences a severe sag, like the one (class 3) foreseen for the industrial environment, in the reference standard for EMC tests EN61000-4-34. By increasing the duration of the sag, the stability limit for the SM is detected. The comparison between the limits—achieved through the linear and nonlinear model, respectively—should show to what extent the nonlinear magnetization matters in modelling stability problems.


## 2.1. A Novel Structure for the SM Model

In [19], it has been proved that by assuming the SM flux linkages or alternatively the currents as state variables, or even considering a combination of them, several different SM equivalent models can be obtained. Nevertheless, in the practical applications, the most recurrent models [17] are those represented in Figure 1. In particular, the model in Figure 1(a) assumes the armature flux linkage vector Ψ as a state variable whereas the model in Figure 1(b) uses the armature currents vector i to that end. The last model is straightforwardly apt to host a LUT that relates the flux linkage vector to the armature currents. In fact, as soon as the machine is considered anhysteretic, the function Ψi becomes single-valued. Hence, it can be univocally approximated, in a domain D of the armature currents, by interpolating a finite number of its samples collected in a LUT. Since the anhysteretic assumption has the collateral drawback of disregarding part of the magnetic losses in the machine, they need to be reintroduced into the model later on. Fortunately, there are many available methods [9] for doing that, and this specific topic will therefore not be dealt in the present paper.Figure 1
Typical structures of a SM model having (a) the flux linkage or (b) the current as a variable of state.
(a)(b)However, choosing the current vector as a state variable procures additional complexity to the model: the construction of the incremental inductance matrix(1)l=∂Ψd∂id∂Ψd∂iq∂Ψq∂id∂Ψq∂iqon one side and the calculation of its inverse l−1 on the other side. Both [16, 17] agree in preferring the model of Figure 1(a) to that of Figure 1(b) because it leads to more simple state equations. Moreover, in ES making use of a full order observer, the flux linkage vector is always available. Conversely, if the observer does not have full rank, the model of Figure 2(b) must be forcibly used according to the state of the art, due to the lack of estimations for the flux linkages.Figure 2
The proposed model structure.The novel model proposed by the present work is represented in Figure2.The following can be observed: (I) it keeps all flux linkages, not only those of the armature, as state variables; (II) it implements the functioniΨ of Figure 1(a) by combining a constant matrix Lσ−1, to be yet defined, and a nonlinear single-valued function Ψmi, which represents the armature main flux linkages. Feature (I) spares the use of (1) in the model as well as the issue of its inversion, whereas feature (II) introduces the novelty of the suggested model. In particular, it does not require the construction of the inverse function iΨ starting from Ψi as shown in [14, 16], and it is applicable to reduced order observers as well, since its saturation model Ψmi works upon available current measurement/estimation rather than on flux linkage ones.However, the exploitation of this model requires two aspects to be clarified in advance. First of all, the machine currents are a linear combination of flux linkages and magnetizing flux linkages according to(2)i=Lσ−1Ψ−CTΨm,with Lσ−1 and C being the constant matrices to be defined. Second, the block diagram within the shaded box in Figure 2 forms an algebraic loop, which sets a convergence problem in a discrete time integration of the related DEs. At the generic integration instant tk, the highlighted block gets the estimated vector of the flux linkages Ψk as an input. Two possible outcomes can occur within the next time instant tk+1: either the generic main fluxes Ψm,j converge to their actual values at tk—providing through that an estimation for the currents ik at the same time—or they diverge making the simulation progress impossible. From a system perspective, the proof of the convergence within the loop is equivalent to demonstrate the existence of a solution i in
(3)i=Lσ−1Ψ−CTLUTi,for all possible Ψ, once the main flux Ψmi is expressed by means of LUTs
(4)Ψmi=LUTi,where
(5)LUTi=LUTdiLUTqi.


## 2.2. Stray Flux Linkages vs. Current Proportionality

A quite general Park model for a wound field salient pole synchronous machine (SPSM) is represented in Figure3, where two armature circuits (outside the shaded circle) and three rotor circuits (inside the shaded circle) are considered. In particular, with reference to the rotor, there are the excitation circuit (f) and a single damper circuit (D) on the d-axis and a single damper circuit (Q) on the q-axis.Figure 3
The coordinated rotor and stator circuits for a wound field SM having single damper circuits on thed- and q-axes, respectively.A more precise consideration of a round rotor SM (RRSM) would have required two more rotor circuits on both axes indeed: one for taking into account the presence of the massive and conductive rotor body and the other one for the eventual conductive slot wedges. In the same way, the consideration of a permanent magnet excited SM (PMSM) should replace the excitation winding in Figure3 with a suitable and permanent source of MMF, in order to fit in the representation. Nevertheless, the circuits of Figure 3 are general enough for presenting the proposed model in its entirety.The more simplistic assumption often done about the Park model of Figure3 is that the circuits along the d-axis share all the same magnetizing flux. In other words, their mutual inductances referred to the armature side are considered all equal to the main armature inductance on the d-axis, Lmd(6)LfD′=Lfa′=LDa′=Lmd.Canay [20], Hiramatsu et al. [21], and Guorui et al. [22] have proved that (6) does not always hold. Not only the main fluxes but also the stray fluxes contribute to the mutual linkage between the circuits laying on the same magnetic axis. Figure 4(a) shows, e.g., that the magnetic coupling between the excitation winding and the D-damper bars in a RRSM takes place through the stray fluxes too. In fact, all rotor circuits, excitation winding and damper bars, are usually laid down in the same slots.Figure 4
Different flux linkages between threed-axis circuits in a wound field synchronous machine having (a) round rotor or (b) salient pole rotor.In this way, the main fluxΦmDf, shared between excitation winding and D-damper bars, becomes larger than Φmd, i.e., the flux they have in common with the armature. In a SPSM having grid-shaped damper bars on each pole, a part Φσaf of the magnetic flux Φmaf—produced by the armature and linked to the field winding—does not link the D-damper bars. Therefore, Figure 4(b) shows that the flux Φmd linking the D-damper circuit can be smaller than the flux Φmaf shared between the armature and the excitation winding. The practical consequence of these remarks is the introduction of an additional stray inductance LfD′−Lmd in the usual Park model as shown in Figure 5. It takes into account the algebraic difference between the total main flux, which links the armature on the d-axis, and the total flux shared between field winding and D-damper circuit.Figure 5
Direct axis equivalent inductance model referred to the armature.According to the observations presented above, it results thatLfD′−Lmd can be negative for SPSM while it is always positive for RRSM.Since the circuits along theq-axis in Figure 3 are only two, there is just one quadrature main flux to be considered. Therefore, the equivalent magnetic circuit on the q-axis looks like the one represented in Figure 6.Figure 6
Quadrature axis equivalent inductance model referred to the armature.Having the circuits of Figures5 and 6 under consideration at once, the following relationship can be set between the currents and the flux linkages at play.
(7)idiqif′iD′iQ′=Lσ−1Ψd−ΨmdΨq−ΨmqΨf′−ΨmdΨD′−ΨmdΨQ′−Ψmq,all quantities used in (7) being referred to the armature side. Hence, the inverse of the 5×5 matrix Lσ, called the matrix of the stray inductance, takes the form
(8)Lσ−1=1Lσa000001Lσa00000LσD′+LmD′−LmdK−LfD′−LmdK000−LfD′−LmdKLσf′+Lmf′−LmdK000001LσQ′,where
(9)K=Lσf′+Lmf′−LmdLσD′+LmD′−Lmd−LfD′−Lmd2.The Blondel coefficient was introduced:(i)
Stray coefficient for the armature(10)σa=1−LmdLd.Laible’s coefficients are as follows [23]:
(i)
Stray coefficient between armature and field winding(11)σaf=1−Laf′2LdLf′(ii)
Stray coefficient between armature andD-damper circuit
(12)σaD=1−LaD′2LdLD′(iii)
Stray coefficient between armature and field winding(13)σfD=1−LfD′2Lf′LD′(iv)
Stray coefficient between armature andQ-damper circuit
(14)σaQ=1−LaQ′2LqLQ′(v)
Armature shielding coefficient(15)μa=1−Laf′LaD′LfD′Ld(vi)
Field winding shielding coefficient(16)μf=1−LfD′Laf′LaD′Lf′(vii)
D-damper circuit shielding coefficient
(17)μD=1−LfD′LaD′Laf′LD′,as well as the machine anisotropy factor at nominal conditions
(18)AFN=LdNLqN−1.The inverse matrix of the stray inductances ((8)) turns into
(19)Lσ−1=1LdN1σa000001σa00000μD/1−σaD+σa+μa−σaμa/1−μaσK2−σa+μa−σaμa/1−μaσK2000−σa+μa−σaμa/1−μaσK2μf/1−σaf+σa+μa−σaμa/1−μaσK2000001+AFNσQ.In (19), LdN represents the armature direct axis inductance at nominal conditions and
(20)σK2=KLd2=μD1−σaD+σa+μa−σaμa1−μaμf1−σaf+σa+μa−σaμa1−μa−σa+μa−σaμa1−μa2.When studying (10)–(17) closer, it can be recognized that the stray and the shielding coefficients can be regarded as constants. In fact, they represent ratios where numerators and denominators vary according to the same underlying nonlinear magnetization characteristics, irrespective to d-axis or q-axis consideration. Adding that LdN and AFN are also constants, it can be concluded that (19) is a constant matrix and that (7) sets a linear proportionality between the vector of the stray fluxes and the vector of the machine currents.


## 2.3. Variable Convergence within the Algebraic Loop

The armature main flux linkages on thed- and q-axes depend on the superposed action of all MMFs at work on the d- and q-axes, respectively. Once the currents of the rotor circuits have been referred to the armature side, they are armature currents at all effects, so that the sources of the main fluxes can be summarized as follows:
(21)iΘdiΘq=∑id′∑iq′=Ci,with
(22)C=1011001001.A more concise expression for (7) can be written by using (22).
(23)i=Lσ−1Ψ−Lσ−1CTΨm,which, once substituted in (21), takes the form
(24)iΘdiΘq=CLσ−1Ψ−CLσ−1CTΨm.For later considerations, it is here anticipated that(25)CLσ−1CT=1Lσaα−100γ−1,where
(26)α=KK+LσaLσD′+LmD′+Lσf′+Lmf′−2Lmd<1,(27)γ=LσQ′LσQ′+Lσa<1.When considered in relative terms, (26) and (27) become
(28)α=σK2σK2+σaμD/1−σaD+μf/1−σaf+2σa+μa−σaμa/1−μa<1,(29)γ=σQ1−1+AFNσa2σQ1−1+AFNσa2+σa1+AFN1−σaQ<1.By considering the discrete values of (5) as one possible choice among the ∞2 sets of values generated by the magnetizing flux functions ΨmdiΘd,iΘq and ΨmqiΘd,iΘq on iΘd,iΘq∈D, it is possible to generalize (5) as
(30)LUTi=ΨmdiΘd,iΘqΨmqiΘd,iΘq,so that the inverse of (30)
(31)LUT−1Ψmd,Ψmq=iΘdΨmd,ΨmqiΘqΨmd,Ψmqexists locally in D−1 as soon as
(32)det∂Ψmd∂iΘd∂Ψmq∂iΘd∂Ψmd∂iΘq∂Ψmq∂iΘq=detlm≠0.By considering (21), (24), and (31) at the same time, it is possible to obtain
(33)Ψm=−CLσ−1CT−1LUT−1Ψm−CLσ−1Ψ.Equation (33) sets a Banach fixed point problem, where the numerical sequence
(34)Ψmn+1=ϕΨmnconverges to a finite value
(35)Ψm=limn⟶∞Ψmn,if and only if, ϕ represents a contracting map of D−1 towards D or simply when
(36)M=supΨm∈Dϕ′Ψm<1.The functionϕ′ stands for the Jacobian of ϕ. Since Ψ is kept constant during the eventual convergence in (33), the Jacobian of ϕ has the following expression:
(37)ϕ′Ψm=−CLσ−1CT−1∂LUT−1∂Ψm.Due to (32) the function Ψm can be inverted, so that (37) turns into
(38)ϕ′Ψm=−CLσ−1CT−1∂LUT∂iΘ−1.Finally, taking into account (25), (30), and (38), the inequality (36) changes into
(39)supΨm∈D−αγLσa2lm,ddlm,qq−lm,dq2<1,where lm,dd, lm,dq, and lm,qq are the elements of the symmetric matrix lm accounting for the incremental magnetizing inductances of the SM.By making use of (39), it is possible to give a proof of the studied convergence (35) within the algebraic loop of the proposed SM model, for both the nonsaturated and saturated machines. In fact, when the SM is not saturated, the incremental inductances (which stand for the slopes of the tangents to the magnetizing curves on the d- and q-axes, respectively) are comparable with the magnetizing inductances (which stand for the slopes of the chords from the origin to the magnetizing curves along the d- and q-axes, respectively). In other terms,
(40)lm,ddlm,qq≅11+AFNLdN2.Moreover, the lack of saturation, i.e., the availability of plenty of nonoriented Weiss domains in the SM magnetic circuit, makes the competition for the magnetic polarization between the coordinated MMFs practically inexistent. That means(41)lm,dq≅0.Taking into account (10) and substituting (40) and (41) in (39) lead to the inequality
(42)Mnosat=−αγσa21+AFN<1.Expression (42) fulfils condition (36) for the convergence. In fact, even assuming an unrealistically large armature stray inductance of about 30%, the d-axis inductance should be at least 10 times larger than the q-axis one in order to have σa21+AFN>1.If the SM is saturated instead, any further differential increase of the coordinated MMFs will not gain any new Weiss domain to their respective directions, so that the incremental magnetizations contributed by the MMFs can be regarded as essentially performed in air. For the same reason, the magnetic anisotropy of the machine fades away makingAF≅1. Therefore,
(43)lm,ddlm,qq≅σa2LdN2.At the same time, the lack of available domains to compete for inhibits the mutual coupling between the coordinated circuits through the cross-magnetization. This results in(44)lm,dq≅0.By substituting (43) and (44) in (39), the fulfilment of (32) for the case of the machine saturation is also achieved, since
(45)Msat=−αγ<1.Due to the local monotonic nature of the magnetization functions (30), it must be expected the transition from (42) to (45) to be also monotonic, to ensure the convergence (35) for every simulated machine behavior between saturation and nonsaturation.


## 2.4. Model Implementation

The model presented in Figure2 has been implemented by means of three functional blocks shown in Figure 7(a). The block performing the integration is the same used in the state of the art and presented in Figures 1(a) and 1(b). It has essentially the goal to gain the state variable Ψ by integrating the unbalances between the impressed voltages, the resistive voltage drops, and the transformational electromotive forces, where present. The nonlinearity manager (NLM) block, in Figure 7(b), receives the current vector as an input and refers the rotor currents to the armature side by making use of the following transformation ratios:
(46)kf=LmfLaf,kD=LmDLaD,kQ=LmQLaQ.Figure 7
Model implementation: (a) fundamental blocks; (b) nonlinearity manager; (c) flux-to-current converter.The transformed currents, which provide their MMFs on the samed- or q-axes, are summed up in order to determine the armature equivalent magnetizing currents iΘd and iΘq, respectively. The LUTd reproduces the function ΨmdiΘd,iΘq shown in Figure 8(a), whereas LUTq reproduces the ΨmqiΘd,iΘq shown in Figure 8(b). One possible way to build the lookup tables for a given SM by using a FEM analysis program is duly described in the appendix to the present work.Figure 8
Functions of the magnetizing flux linkages according to (a) thed-axis and (b) the q-axis.
(a)(b)As soon as the coordinated magnetizing current is introduced in the LUTs, two values for the armature main flux linkages are returned as outputs,Ψmd and Ψmq, respectively.These quantities provide the inputs for the flux-to-current (FTC) block in Figure7(c) together with the vector of the flux linkages Ψ. The rotor-related flux linkages must be referred to the armature side by means of the ratios (46). Then, as already shown in (7), the block performs a linear combination of the linked stray fluxes through the matrix Lσ−1, in order to get an estimate of the machine currents. The currents obtained by that way are all referred to the armature side. Therefore, the rotor-related currents need to be referred back to the rotor by using the usual ratios (46).


## 2.5. Performed Test Simulations

The proposed nonlinear model for SM has been tested on the wound field SPSM (WFSPSM) shown in Figure9, the parameters of which are given in Table 1.Figure 9
A cross-section of the test SM with magnetic field lines produced by (a) the armature direct current and (b) the armature quadrature current.
(a)(b)Table 1
Main parameters of the test SM.
SymbolDescriptionValueANRated power60 kVAVNRated voltage400 VnNRated speed1500 rpmfNRated frequency50 HzcosφNRated power factor0.9δmin⋯δmaxAirgap range at the pole shoe3.2÷5.7 mmαpPole coverage0.76liActive axial length184 mmrδAirgap average radius161 mmqSlots per pole per phase6nsConductors per slot8rArmature phase resistance (p.u.)6.2·10−2rfField winding resistance (p.u.)5.8·10−3rDd-axis damp. circuit resistance (p.u.)7.7·10−3rQq-axis damp. circuit resistance (p.u.)15.3·10−3xdNNominald-axis reactance (p.u.)1.32xqNNominalq-axis reactance (p.u.)0.72AFNNominal anisotropy factor0.83σaArmature stray factor5.2%σafField winding to armature stray factor9.7%σaDArmature tod-axis damper circuit stray factor11.8%σQq-axis damper circuit stray factor13.3%μaArmature shielding coefficient1.7%μfField winding shielding coefficient9.5%μDq-axis damper circuit shielding coefficient11.7%In order to prove the capability of the suggested model in handling the nonlinearity and the cross-magnetization in the test machine, the model itself has been compared with a linear and anisotropic one, where the LUTs of Figure7(b) have been replaced by constant magnetizing inductances Lmd and Lmq. The latter has been chosen, for each led test, so as to match the initial conditions of the nonlinear model, before an external cause intervenes to push the machine away from its initial steady state. In this way, the two models are aligned at the beginning of the test, and all subsequent divergences between them can be recognized more clearly. In particular, four tests have been envisaged, that, for the nature of the chosen perturbation, end up changing the magnetizing flux and highlighting the magnetic nonlinearity of the machine.The first test foresees the sudden drop of the excitation voltage from its nominal value down to 40% of it, while the nominal driving torque of the generator is kept constant. Therefore, the generated nominal active power remains unchanged. Due to the lack of exciting MMF, the machine must require the magnetizing current from the grid, moving by that from the generation to the consumption of reactive power.The second test is based on the observation that the voltage level at the generator bus bar is responsible for the incoming/outgoing reactive power. Therefore, a grid voltage increase, when the excitation current is kept constant, results in a higher magnetizing flux in the SM. This should in turn shift the generator towards the saturation and possibly highlight its underlying nonlinear behavior.The third test reproduces the adjustment of the excitation current when the generated power is suddenly reduced down to 25% of its nominal value. In order to keep the power factor constant, the excitation current is decreased accordingly, so that the exciting MMF is drastically reduced.The fourth and final simulations deal with a low-voltage fault ride through (LVFRT), where the grid voltage experiences a severe sag, like the one (class 3) foreseen for the industrial environment, in the reference standard for EMC tests EN61000-4-34. By increasing the duration of the sag, the stability limit for the SM is detected. The comparison between the limits—achieved through the linear and nonlinear model, respectively—should show to what extent the nonlinear magnetization matters in modelling stability problems.


## 3. Results and Discussion

### 3.1. Excitation Voltage Step Test

Figure10 shows the instantaneous drop of the SM excitation voltage, at the time t=3 s, in the test machine providing the nominal power at the nominal conditions. This change in the control voltage produces a subsequent and slower drop of the excitation current.Figure 10
Change in the excitation current produced by a temporary reduction of the control voltage.It is not possible to observe any substantial difference between the resulting field currents for the nonlinear and linear model, since thed-damper bars screen the field winding from the large variations of the armature MMFs occurring during the transient. Figure 11 shows that the d-axis component of the armature current drops almost to zero whereas the q-axis one increases.Figure 11
d-axis and q-axis current responses due to the adjustment of the field voltage.Since the electromagnetic torque remains unchanged, from the Park equation for the torque follows that(47)TEM=32pΨdiq−Ψqid≅32pΨdiq,where Ψd, at steady state and for id≅0, is
(48)Ψd≅Lmdif′.In Figure11, it is possible to observe that iq in the linear model is larger than that in the nonlinear one. Since the excitation current if′ is the same for the two models (Figure 10), by considering (47) and (48), it follows
(49)LmdLUTiqLUT=Lmdliniqlin,where
(50)LmdLUT>Lmdlin,because of
(51)iqLUT<iqlin.By remembering that the magnetizing inductances of the linear model are kept constant and equal to the initial magnetizing inductances of the nonlinear model, from (50) descends
(52)LmdLUT3s<t<5s>LmdLUTt<5s.Inequality (52) proves that the magnetizing inductance of the nonlinear model senses the transition from the over excitation to the under excitation, since its value increases when the magnetizing flux decreases. The transition from the over excitation to the under excitation is given also from the graph of the reactive power in Figure 12.Figure 12
Effect of the under excitation on the reactive power.The SM produces reactive power as long as the rotor MMF is strong enough. As soon as the excitation current drops, the generator starts consuming reactive power. It can also be observed that, sinceid≅0,
(53)Q≅32Lmqiq2.The larger value of theiq in the linear model observed in (51) together with (53) explains the discrepancy between the answers of the two models shown in Figure 12.


### 3.2. Grid Voltage Step Test

In this test, while the generator is providing 60% of the nominal power atcosφ=0.9 (lagging), a stepwise 10% increase of the grid voltage occurs at t=2s. In Figure 13, it can be observed that the two models give a different account of the d-axis and q-axis currents in the final steady state, as well as of the generator load angle, as it is shown in Figure 14. In particular, Figure 14 shows, firstly, that the load angle δ decreases moving from the initial steady state to the final one and, secondly, that the linear model provides a smaller load angle for the final steady state, δlin, than the one given by the nonlinear model, δLUT.Figure 13
Effect of a stepwise +10% grid voltage increase onid andiq.Figure 14
Load angle perturbation due to a +10% stepwise increase of the grid voltage at unchanged excitation and generated active power.Moreover, since the excitation of the SM and its driving torque are kept constant, the generated active power does not change when the transient response fades away, whereas the reactive power does (Figure15). The nonlinear model gives account of a drop in the reactive power provided by the generator, which is larger than the one reported by the linear model.Figure 15
Effect of a stepwise grid voltage increase on the generated active and reactive power.Figure16 helps to interpret and check the different results shown in Figures 14 and 15.Figure 16
Shift of the magnetic work point in the models during the present test.The letterA indicates the magnetic work point shared by the two models at the initial steady state, when the same q-axis voltage sets the same magnetizing flux linkage on the d-axis according to
(54)Vcosδ≅ωΨmdA.In (54), the voltage drop on the armature stray inductance is neglected, since it is irrelevant to the present line of reasoning.At the initial work pointA, the main flux linkage is set by the superposition of the MMF proportional to if′ and id, respectively. As soon as the grid voltage is incremented by ΔV>0, the magnetic flux linkage in the linear model increases up to
(55)V+ΔVcosδlin≅ωΨmdB>ωΨmdA,due to 0<δlin<δ in Figure 14. Since the model is linear, the MMF responsible for ΨmdB can be found by the vertical projection of the point B on the x-axis, which is taken on the straight s at y=ωΨmdB. It can be observed that the resulting MMF corresponding to the work point B gives a correct account of a drop of the direct axis current in the final steady state, since Δidlin<0 for graphic construction.Coming to the nonlinear model in Figure16, an increase ΔV>0 of the grid voltage produces an increase of the magnetic flux linkage up to
(56)V+ΔVcosδLUT≅ωΨmdC>ωΨmdA,due to 0<δLUT<δ in Figure 14. At the same time, it must be
(57)ωΨmdC<ωΨmdB,due to 0<δlin<δLUT in Figure 14. The MMF, which sets ΨmdC, can be determined by the vertical projection of the point C on the x-axis, which is taken on the nonlinear magnetization curve c in Figure 16 at y=ωΨmdC. It can be noticed that the resulting MMF corresponding to the work point C gives a correct account of a drop in the direct axis current for the final steady state, since ΔidLUT<0. Moreover, the graphic construction shows that ΔidLUT≪Δidlin in spite of (57), just in reason of the magnetic nonlinearity of the SM taken into account in the model containing the LUTs. The much stronger reduction of the id in the nonlinear model rather than in the linear one is the reason of the opposite tendency of the iq in the two models, shown in Figure 13. In fact, the constant machine torque during this test
(58)T≈Ψmdiq+Ψmqidsets a relationship between id and iq. The small decrease of the id in the linear model does not change the torque effectively whereas the iq is forced to decrease in order to compensate the relevant linear increase in Ψmd. On the contrary, in the nonlinear model, the drop in id is substantial, and due to the small increase of the main flux linkage Ψmd, the iq must be larger than in the linear model for keeping the torque constant.


### 3.3. Driving Torque Step Test (at Constant Power Factor)

The present test foresees a driving torque drop from 100% down to 25% of its nominal value within 0.1 s. At the same time, in order to target the same initial power factor (0.9 lagging) in the final steady state, the excitation current is decreased accordingly. Figure17(a) shows the needed adjustments of the field current whereas Figure 17(b) and the generated active power.Figure 17
The adjustment of the excitation (a) and the active power (b).Even though Figure17 does not present relevant differences between the answers of the two models, in Figure 18, they give a different account for the envelope of the armature phase current amplitude.Figure 18
Timely envelope of the phase current amplitude.The mismatch between the two answers depends mainly on the direct axis current, as Figure19 reveals by showing the different steady-state values for the reactive power.Figure 19
The different account of the reactive power generated by the SM.The difference between the values of theid in the two-models can be justified by the Park armature voltage equation along the q-axis at steady state, where the d-axis current is given by
(59)−id≅ωΨmd−VcosδωLσa.Since the load angleδ decreases due to the drastic drop of the generated active power, so that cosδ≅1, the Ψmd plays the main role in producing the difference between the two steady-state values of id provided by the models. Figure 20 shows how the magnetic work point in the linear model moves from A to B along the straight s, until the main flux ΨmdB is set through the constant magnetizing inductance LmdN and the superposed effects of the currents if′B and idlinB.Figure 20
Shift of the magnetic work point in the models during the present test.Equation (59) shows that for ωΨmd to exactly balance the quadrature axis voltage at id=0, the main flux Ψmd must be set by if′ only. This is way, by following the straight w (having slope −ωLσa) down to the point H at y=Vcosδ, it is possible to find the current if′B. Since the final value of the excitation current if′B and the quadrature axis voltage Vcosδ is shared by the two models, by prolonging w up to the nonlinear magnetization characteristic c, the final magnetic work point C for the nonlinear model is found. By comparing the points C and B, it can be recognized that
(60)LmdC>LmdA=LmdN,which gives a right account of the nonlinear property of the SM magnetic circuit. Furthermore, it can be found that
(61)−idLUT>−idlin,which explains why neglecting the nonlinearity in the SM model would lead to underestimating the value of the reactive power provided at low excitation.


### 3.4. Low-Voltage Fault Ride Through

In the following last test, two voltage sags with the same intensity (-60% of the nominal grid voltage) but different duration (14 and 18 cycles, respectively) are applied to the armature of the SM. The response of the generator is simulated by means of the two models under comparison. This test, which evaluates the capability of the SM to keep the pace when the disturbance is over, represents a LVFRT test. Figure21 shows that both machines keep the pace for a sag duration of 14 cycles. No substantial differences between the final steady-state values of the load angle can be recognized.Figure 21
Low-voltage fault ride through test with different sag durations.In the same figure, for a duration of the sag of 18 cycles, the nonlinear model proves the SM to be able to keep the pace after the sag, whereas the linear model does not.The different results provided by the two models can be interpreted in an intuitive way through the application of the well-known equivalent area criterion, as shown in Figure22. In there, the difference between the constant driving torque Tm and the electromagnetic torque T is responsible for the acceleration work on the rotor. The electromagnetic torque can be expressed by
(62)T=pω3E′VXd′sinδ+V22Xd′Xd′Xq′−1sin2δ,where Xd′ and Xq′ are the transient reactances of the SM along the d- and q-axes, respectively. Since the transient reactances are related to the stator and rotor stray inductances essentially, the nonlinearity of the magnetic circuit does not influence their values in a substantial way. Therefore, the different performances of the two models must be caused by their different transient electromotive force E′. In particular, in order to justify the results of Figure 21, it must be true that, during the voltage sag,
(63)ELUT′>Elin′.Figure 22
Equal area criterion used for explaining the transient stability. The continuous bell-shaped curve represents the synchronous torque in p.u. whereas the shaded one is the transient torque in p.u.: (a) 14-cycle sag-linear modelδmax=δL; (b) 14-cycle sag-non-linear model δ′max<δL; (c) 18-cycle sag-linear model Alin++Alin−>0 for δ=δL; (d) 18-cycle sag-non-linear model δ′max≤δL.
(a)(b)(c)(d)In studying Figure22, it emerges that the nonlinearity of the SM magnetic circuit confers to the machine itself, a given capability to withstand grid voltage sags. Neglecting the said nonlinearity in the SM model leads to underestimating the machine performances during LVFRT tests.


### 3.5. Anisotropy and Cross-Magnetization: 1 or 2 LUT?

A relevant question considers the possibility of using just one single LUT instead of two for expressing both anisotropy and cross-magnetization in the SM. In that case, theLUTd related to Ψmd could be used in first place. In order to obtain Ψmq from LUTd, its values should be scaled down, e.g., by a suitable constant which depends on the anisotropy factor AF. As long as the LUTd represents a function of id and iq, some information about the influence of the q-axis MMF on the cross-magnetization must be contained in LUTd. In fact, some energy considerations presented in [15, 24] make it possible to infer that
(64)ldq=∂Ψd∂iq=lqd=∂Ψq∂d.Therefore, as soon as(65)LUTq=11+AFLUTdT,it is evident that identity (64) is violated, because, from (65), it follows
(66)lqd=11+AFldq.Hence, it seems impossible to give a correct account of the cross-magnetization by using the information coming from just one single LUT.With reference to the anisotropy factor AF instead, in many works, the ratio betweenLd and Lq is kept constant, even though it proves to be strongly dependent on the intensities of id and iq and on their relative ratio. In Figure 23, the AFm (see the appendix for its definition) for the test SM of Figure 9 has been represented. It gives an evidence of the fact that the anisotropy of a SM is in general not constant. It tends to fade away for very deep saturation, since the magnetic circuit reluctance along the d- and q-axes becomes comparable in magnitude. Even more important is the fact that the cross-magnetization contributes to foster the anisotropy of the machine. It is possible to observe this fact in Figure 23, where the anisotropy factor shows higher values along two crossed versants, which run along the 1st to 3rd and 2nd to 4thid-iq quadrant directions. This is probably due to the competing action of the coordinated MMFs on the material polarization, which prevents one axis to reach the full saturation when the other one is excited at the same time.Figure 23
Magnetic anisotropy factor for the test SM of Figure9.


### 3.6. An Estimate for the Convergence Speed of the Algebraic Loop in the Model

Equation (33) represents a contracting map ϕD⟶D when the SM iron is unsaturated (42) as well as when it is not (43). Therefore, it is possible to calculate the number k (in excess) of the needed iterations in order to get Ψmk close enough to its asymptotic value Ψ¯m, according to a given relative accuracy ε. With M representing the upper bound of the Jacobian ϕ′ over D, it results as follows:
(67)ε=Ψ¯m−ΨmkΨ¯m−Ψm0≤Mk.After calculating the upper boundMsat=0.351 in (45) by means of (28) and (29), the maximal number of iterations (67) reducing the initial error about one thousand times ε=10−3 is k=7.


## 3.1. Excitation Voltage Step Test

Figure10 shows the instantaneous drop of the SM excitation voltage, at the time t=3 s, in the test machine providing the nominal power at the nominal conditions. This change in the control voltage produces a subsequent and slower drop of the excitation current.Figure 10
Change in the excitation current produced by a temporary reduction of the control voltage.It is not possible to observe any substantial difference between the resulting field currents for the nonlinear and linear model, since thed-damper bars screen the field winding from the large variations of the armature MMFs occurring during the transient. Figure 11 shows that the d-axis component of the armature current drops almost to zero whereas the q-axis one increases.Figure 11
d-axis and q-axis current responses due to the adjustment of the field voltage.Since the electromagnetic torque remains unchanged, from the Park equation for the torque follows that(47)TEM=32pΨdiq−Ψqid≅32pΨdiq,where Ψd, at steady state and for id≅0, is
(48)Ψd≅Lmdif′.In Figure11, it is possible to observe that iq in the linear model is larger than that in the nonlinear one. Since the excitation current if′ is the same for the two models (Figure 10), by considering (47) and (48), it follows
(49)LmdLUTiqLUT=Lmdliniqlin,where
(50)LmdLUT>Lmdlin,because of
(51)iqLUT<iqlin.By remembering that the magnetizing inductances of the linear model are kept constant and equal to the initial magnetizing inductances of the nonlinear model, from (50) descends
(52)LmdLUT3s<t<5s>LmdLUTt<5s.Inequality (52) proves that the magnetizing inductance of the nonlinear model senses the transition from the over excitation to the under excitation, since its value increases when the magnetizing flux decreases. The transition from the over excitation to the under excitation is given also from the graph of the reactive power in Figure 12.Figure 12
Effect of the under excitation on the reactive power.The SM produces reactive power as long as the rotor MMF is strong enough. As soon as the excitation current drops, the generator starts consuming reactive power. It can also be observed that, sinceid≅0,
(53)Q≅32Lmqiq2.The larger value of theiq in the linear model observed in (51) together with (53) explains the discrepancy between the answers of the two models shown in Figure 12.


## 3.2. Grid Voltage Step Test

In this test, while the generator is providing 60% of the nominal power atcosφ=0.9 (lagging), a stepwise 10% increase of the grid voltage occurs at t=2s. In Figure 13, it can be observed that the two models give a different account of the d-axis and q-axis currents in the final steady state, as well as of the generator load angle, as it is shown in Figure 14. In particular, Figure 14 shows, firstly, that the load angle δ decreases moving from the initial steady state to the final one and, secondly, that the linear model provides a smaller load angle for the final steady state, δlin, than the one given by the nonlinear model, δLUT.Figure 13
Effect of a stepwise +10% grid voltage increase onid andiq.Figure 14
Load angle perturbation due to a +10% stepwise increase of the grid voltage at unchanged excitation and generated active power.Moreover, since the excitation of the SM and its driving torque are kept constant, the generated active power does not change when the transient response fades away, whereas the reactive power does (Figure15). The nonlinear model gives account of a drop in the reactive power provided by the generator, which is larger than the one reported by the linear model.Figure 15
Effect of a stepwise grid voltage increase on the generated active and reactive power.Figure16 helps to interpret and check the different results shown in Figures 14 and 15.Figure 16
Shift of the magnetic work point in the models during the present test.The letterA indicates the magnetic work point shared by the two models at the initial steady state, when the same q-axis voltage sets the same magnetizing flux linkage on the d-axis according to
(54)Vcosδ≅ωΨmdA.In (54), the voltage drop on the armature stray inductance is neglected, since it is irrelevant to the present line of reasoning.At the initial work pointA, the main flux linkage is set by the superposition of the MMF proportional to if′ and id, respectively. As soon as the grid voltage is incremented by ΔV>0, the magnetic flux linkage in the linear model increases up to
(55)V+ΔVcosδlin≅ωΨmdB>ωΨmdA,due to 0<δlin<δ in Figure 14. Since the model is linear, the MMF responsible for ΨmdB can be found by the vertical projection of the point B on the x-axis, which is taken on the straight s at y=ωΨmdB. It can be observed that the resulting MMF corresponding to the work point B gives a correct account of a drop of the direct axis current in the final steady state, since Δidlin<0 for graphic construction.Coming to the nonlinear model in Figure16, an increase ΔV>0 of the grid voltage produces an increase of the magnetic flux linkage up to
(56)V+ΔVcosδLUT≅ωΨmdC>ωΨmdA,due to 0<δLUT<δ in Figure 14. At the same time, it must be
(57)ωΨmdC<ωΨmdB,due to 0<δlin<δLUT in Figure 14. The MMF, which sets ΨmdC, can be determined by the vertical projection of the point C on the x-axis, which is taken on the nonlinear magnetization curve c in Figure 16 at y=ωΨmdC. It can be noticed that the resulting MMF corresponding to the work point C gives a correct account of a drop in the direct axis current for the final steady state, since ΔidLUT<0. Moreover, the graphic construction shows that ΔidLUT≪Δidlin in spite of (57), just in reason of the magnetic nonlinearity of the SM taken into account in the model containing the LUTs. The much stronger reduction of the id in the nonlinear model rather than in the linear one is the reason of the opposite tendency of the iq in the two models, shown in Figure 13. In fact, the constant machine torque during this test
(58)T≈Ψmdiq+Ψmqidsets a relationship between id and iq. The small decrease of the id in the linear model does not change the torque effectively whereas the iq is forced to decrease in order to compensate the relevant linear increase in Ψmd. On the contrary, in the nonlinear model, the drop in id is substantial, and due to the small increase of the main flux linkage Ψmd, the iq must be larger than in the linear model for keeping the torque constant.


## 3.3. Driving Torque Step Test (at Constant Power Factor)

The present test foresees a driving torque drop from 100% down to 25% of its nominal value within 0.1 s. At the same time, in order to target the same initial power factor (0.9 lagging) in the final steady state, the excitation current is decreased accordingly. Figure17(a) shows the needed adjustments of the field current whereas Figure 17(b) and the generated active power.Figure 17
The adjustment of the excitation (a) and the active power (b).Even though Figure17 does not present relevant differences between the answers of the two models, in Figure 18, they give a different account for the envelope of the armature phase current amplitude.Figure 18
Timely envelope of the phase current amplitude.The mismatch between the two answers depends mainly on the direct axis current, as Figure19 reveals by showing the different steady-state values for the reactive power.Figure 19
The different account of the reactive power generated by the SM.The difference between the values of theid in the two-models can be justified by the Park armature voltage equation along the q-axis at steady state, where the d-axis current is given by
(59)−id≅ωΨmd−VcosδωLσa.Since the load angleδ decreases due to the drastic drop of the generated active power, so that cosδ≅1, the Ψmd plays the main role in producing the difference between the two steady-state values of id provided by the models. Figure 20 shows how the magnetic work point in the linear model moves from A to B along the straight s, until the main flux ΨmdB is set through the constant magnetizing inductance LmdN and the superposed effects of the currents if′B and idlinB.Figure 20
Shift of the magnetic work point in the models during the present test.Equation (59) shows that for ωΨmd to exactly balance the quadrature axis voltage at id=0, the main flux Ψmd must be set by if′ only. This is way, by following the straight w (having slope −ωLσa) down to the point H at y=Vcosδ, it is possible to find the current if′B. Since the final value of the excitation current if′B and the quadrature axis voltage Vcosδ is shared by the two models, by prolonging w up to the nonlinear magnetization characteristic c, the final magnetic work point C for the nonlinear model is found. By comparing the points C and B, it can be recognized that
(60)LmdC>LmdA=LmdN,which gives a right account of the nonlinear property of the SM magnetic circuit. Furthermore, it can be found that
(61)−idLUT>−idlin,which explains why neglecting the nonlinearity in the SM model would lead to underestimating the value of the reactive power provided at low excitation.


## 3.4. Low-Voltage Fault Ride Through

In the following last test, two voltage sags with the same intensity (-60% of the nominal grid voltage) but different duration (14 and 18 cycles, respectively) are applied to the armature of the SM. The response of the generator is simulated by means of the two models under comparison. This test, which evaluates the capability of the SM to keep the pace when the disturbance is over, represents a LVFRT test. Figure21 shows that both machines keep the pace for a sag duration of 14 cycles. No substantial differences between the final steady-state values of the load angle can be recognized.Figure 21
Low-voltage fault ride through test with different sag durations.In the same figure, for a duration of the sag of 18 cycles, the nonlinear model proves the SM to be able to keep the pace after the sag, whereas the linear model does not.The different results provided by the two models can be interpreted in an intuitive way through the application of the well-known equivalent area criterion, as shown in Figure22. In there, the difference between the constant driving torque Tm and the electromagnetic torque T is responsible for the acceleration work on the rotor. The electromagnetic torque can be expressed by
(62)T=pω3E′VXd′sinδ+V22Xd′Xd′Xq′−1sin2δ,where Xd′ and Xq′ are the transient reactances of the SM along the d- and q-axes, respectively. Since the transient reactances are related to the stator and rotor stray inductances essentially, the nonlinearity of the magnetic circuit does not influence their values in a substantial way. Therefore, the different performances of the two models must be caused by their different transient electromotive force E′. In particular, in order to justify the results of Figure 21, it must be true that, during the voltage sag,
(63)ELUT′>Elin′.Figure 22
Equal area criterion used for explaining the transient stability. The continuous bell-shaped curve represents the synchronous torque in p.u. whereas the shaded one is the transient torque in p.u.: (a) 14-cycle sag-linear modelδmax=δL; (b) 14-cycle sag-non-linear model δ′max<δL; (c) 18-cycle sag-linear model Alin++Alin−>0 for δ=δL; (d) 18-cycle sag-non-linear model δ′max≤δL.
(a)(b)(c)(d)In studying Figure22, it emerges that the nonlinearity of the SM magnetic circuit confers to the machine itself, a given capability to withstand grid voltage sags. Neglecting the said nonlinearity in the SM model leads to underestimating the machine performances during LVFRT tests.


## 3.5. Anisotropy and Cross-Magnetization: 1 or 2 LUT?

A relevant question considers the possibility of using just one single LUT instead of two for expressing both anisotropy and cross-magnetization in the SM. In that case, theLUTd related to Ψmd could be used in first place. In order to obtain Ψmq from LUTd, its values should be scaled down, e.g., by a suitable constant which depends on the anisotropy factor AF. As long as the LUTd represents a function of id and iq, some information about the influence of the q-axis MMF on the cross-magnetization must be contained in LUTd. In fact, some energy considerations presented in [15, 24] make it possible to infer that
(64)ldq=∂Ψd∂iq=lqd=∂Ψq∂d.Therefore, as soon as(65)LUTq=11+AFLUTdT,it is evident that identity (64) is violated, because, from (65), it follows
(66)lqd=11+AFldq.Hence, it seems impossible to give a correct account of the cross-magnetization by using the information coming from just one single LUT.With reference to the anisotropy factor AF instead, in many works, the ratio betweenLd and Lq is kept constant, even though it proves to be strongly dependent on the intensities of id and iq and on their relative ratio. In Figure 23, the AFm (see the appendix for its definition) for the test SM of Figure 9 has been represented. It gives an evidence of the fact that the anisotropy of a SM is in general not constant. It tends to fade away for very deep saturation, since the magnetic circuit reluctance along the d- and q-axes becomes comparable in magnitude. Even more important is the fact that the cross-magnetization contributes to foster the anisotropy of the machine. It is possible to observe this fact in Figure 23, where the anisotropy factor shows higher values along two crossed versants, which run along the 1st to 3rd and 2nd to 4thid-iq quadrant directions. This is probably due to the competing action of the coordinated MMFs on the material polarization, which prevents one axis to reach the full saturation when the other one is excited at the same time.Figure 23
Magnetic anisotropy factor for the test SM of Figure9.


## 3.6. An Estimate for the Convergence Speed of the Algebraic Loop in the Model

Equation (33) represents a contracting map ϕD⟶D when the SM iron is unsaturated (42) as well as when it is not (43). Therefore, it is possible to calculate the number k (in excess) of the needed iterations in order to get Ψmk close enough to its asymptotic value Ψ¯m, according to a given relative accuracy ε. With M representing the upper bound of the Jacobian ϕ′ over D, it results as follows:
(67)ε=Ψ¯m−ΨmkΨ¯m−Ψm0≤Mk.After calculating the upper boundMsat=0.351 in (45) by means of (28) and (29), the maximal number of iterations (67) reducing the initial error about one thousand times ε=10−3 is k=7.


## 4. Conclusions

A two-axis theory-based SM model has been presented, which takes into account the nonlinearity of the SM magnetic circuit as well as the cross-magnetization. These performances are achieved by using two LUTs that reproduce the magnetizing flux linkages of the machine armature. Following this authors’ choice, the outlined model does not need any inversion of the current-to-flux-linkage LUT, even if the model maintains the flux linkages as state variables. This is possible thanks to a LUT-centered algebraic loop, which is intrinsic into the model and whose iterations always converge to the instantaneous machine main flux. The proof for the convergence of the introduced algebraic loop has been given for both cases of nonsaturated and saturated machines. The novel SM model has been successfully tested for the simulation of a wholly specified salient pole synchronous generator. Four different transient scenarios have been considered, where as many external disturbances provoke a substantial perturbation on the machine main fluxes, so to highlight the magnetic nonlinearity at play in the model. The outcomes have been compared with those of a traditional linear Park model for the same SM, giving a punctual interpretation and explanation of the occurred discrepancies.All tests performed have shown that neglecting the nonlinearity of the machine magnetization characteristic leads to remarkable divergences between the two model outcomes. Finally, the construction of the main flux linkages LUTs for the test SM has highlighted the importance of a correct account of the machine magnetic anisotropy. In particular, the assumption of a constant ratio between thed-axis and q-axis armature inductances cannot be representative of all machine behaviors. In fact, the anisotropy of the SM depends not only on the absolute values of the direct and quadrature armature currents but also on their ratio. Moreover, it has been found out that the cross-magnetization leaves a recognizable signature on the machine magnetic anisotropy profile and that this fact cannot be reproduced by using one single LUT in the model. For all exposed reasons, the authors believe that the proposed model can be of a certain interest for those who perform simulations of SM in classic power systems applications, as well as in the control tasks of the modern electric drives. In this last field, this model can be used in reduced order observers, since its underlying mechanism for reproducing the magnetic nonlinearity relies on the use of the measured currents rather than on the flux linkages.



---

*Source: 4547086-2023-02-18.xml*