Journal of Modern Power Systems and Clean Energy

ISSN 2196-5625 CN 32-1884/TK

网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

Integrated Equivalent Model of Permanent Magnet Synchronous Generator Based Wind Turbine for Large-scale Offshore Wind Farm Simulation  PDF

  • Ming Zou
  • Yan Wang
  • Chengyong Zhao
  • Jianzhong Xu
  • Xiaojiang Guo
  • Xu Sun
State Key Laboratory of Alternate Electrical Power System with Renewable Energy Sources, North China Electric Power University (NCEPU), Beijing China;; China Huaneng Group Clean Energy Technology Research Institute Co., Ltd., Beijing, China

Updated:2023-09-20

DOI:10.35833/MPCE.2022.000495

  • Full Text
  • Figs & Tabs
  • References
  • Authors
  • About
CITE
OUTLINE

Abstract

The high-speed simulation of large-scale offshore wind farms (OWFs) preserving the internal machine information has become a huge challenge due to the large wind turbine (WT) count and microsecond-range time step. Hence, it is undoable to investigate the internal node information of the OWF in the electro-magnetic transient (EMT) programs. To fill this gap, this paper presents an equivalent modeling method for large-scale OWF, whose accuracy and efficiency are guaranteed by integrating the individual devices of permanent magnet synchronous generator (PMSG) based WT. The node-elimination algorithm is used while the internal machine information is recursively updated. Unlike the existing aggregation methods, the developed EMT model can reflect the characteristics of each WT under different wind speeds and WT parameters without modifying the codes. The access to each WT controller is preserved so that the time-varying dynamics of all the WTs could be simulated. Comparisons of the proposed model with the detailed model in PSCAD/EMTDC have shown very high precision and high efficiency. The proposed modeling procedures can be used as reference for other types of WTs once the structures and parameters are given.

I. Introduction

DEALING with the climate change issue under “carbon emission peak and carbon neutrality” goals, the offshore wind power will get rapid development. Since the large-scale offshore wind farm (OWF) is an islanding power generation system, the characteristics of the AC system will undergo profound changes considering the high permeability of OWFs. Also, the AC system may encounter safety and operation difficulties, e.g., instability and faults [

1]-[4].

In order to accurately evaluate the impact of the integration of large-scale OWF into the AC system, the high-speed and accurate electro-magnetic transient (EMT) simulations and the equivalent modeling methods are necessary [

5], [6]. However, there are hundreds of wind turbines (WTs) in an OWF, and each WT contains dozens of electrical nodes. If a detailed EMT model is used to simulate the OWF, the ultra-high-order matrix inversion in the EMT program, e.g., PSCAD/EMTDC in this paper, would be a great challenge [7], [8]. Meanwhile, a large number of power electronic switches are included into each WT and thus very small simulation time step is required to reproduce the dynamic switching behaviors. This would in turn make the EMT simulation more time-consuming [9], [10]. Therefore, the researcher often ignores the interaction between the internal machine dynamics of the OWF and the AC system.

To address this issue, the shifted frequency phasor (SFP) modeling [

11], [12], full aggregation [13], “multi-machine” [14], and operation modeling [15] methods are proposed. The SFP modeling method works well under large time steps and has been applied to simulate wind farms [12]. However, the error increases significantly when the signal deviates from rated frequency. The full aggregation method [13] is simple and efficient, but it is not sufficiently accurate because it ignores the wake effects. The error becomes larger when most of the WTs in an OWF operate below the rated power. To surmount these drawbacks, a “multi-machine” method is proposed in [14]. The eigenvalue-based method is accurate by clustering the WTs. Nevertheless, the eigenvalue is supposed to change with the wind speeds and make the aggregated model mismatch with the original model. An adaptive robust optimization model via operation aspects is formulated in [15]. It is suitable for the optimal scheduling by relevant operators, but the acquisition of operation information is difficult for other researchers.

The existing methods have attempted to find an optimal trade-off between the accuracy of internal node information and the efficiency of the EMT simulations. Except these OWF model-based methods, some general-purpose techniques for accelerating the EMT simulations are also applicable. For example, [

16] proposes a fine-grained electrical circuit partitioning method, in which the graphics processing unit (GPU) based parallel architecture is utilized. A distributed hybrid interface of EMT and transient stability model is proposed in [17]. In [18], the system is divided into several subsystems with different time steps, and the subsystems are simulated separately with different time steps. Although these methods avoid the ultra-high-order matrix inversion of the original circuit, the WT simulation speed is still slow, which needs acceleration. The real-time all-digital simulation is also a useful tool for OWF simulation [19]. However, it requires a large number of processor cores and may not be affordable for users. In [20], a classical Thévenin’s equivalent EMT simulation model for modular multilevel converter (MMC) in high-voltage direct current (HVDC) transmission system is proposed. The employed nested fast and simultaneous solution (NFSS) algorithm can speed up the EMT simulation of electrical circuits with modular structures by eliminating internal nodes.

This paper aims to provide an electrical node-elimination based time-efficient EMT model for large-scale OWF, and an algorithm similar to NFSS used in [

20] will be applied. However, the NFSS algorithm will be further specified in this paper by considering the complicated structure of the individual WT, which includes several types of AC and DC equipment in cascade. The essence of the proposed modeling method is that it partitions the admittance matrix of WT and derives an efficient time-varying Thévenin’s equivalent model for the WT. The new method is exactly equivalent to the detailed EMT model of the WT with much reduced computational burden.

The performance of the proposed method is compared with the modeling methods in [

12]-[15] in terms of system-level accuracy, internal accuracy, speed-up capability and applicability, as shown in Table I.

TABLE I  Comparisons of Different Modeling Methods
ReferenceSystem-level accuracyInternal accuracySpeed-up capabilityApplicability
Proposed ★★★★★ ★★★★★ ★★★★ ★★★★
[12] ★★★ ★★★★ ★★★ ★★★★
[13] ★★★★★ ★★
[14] ★★★★ ★★★ ★★★★ ★★★
[15] ★★★★ ★★★ ★★ ★★★

Note:   more stars ★ indicate better performance.

II. Model Description

The structure of a typical OWF is shown in Fig. 1. The OWF consists of four series-parallel collecting strings. Each string includes five WTs. A 35 kV AC collecting network is used to collect the power of all WTs. Then, a three-phase double-winding 35 kV/230 kV transformer is used to boost the AC voltage to a transmission level. Finally, the collected wind power is transmitted to the AC grid via a 30 km 220 kV AC cable.

Fig. 1  System structure of a typical OWF.

The permanent magnet synchronous generator (PMSG) based WT shown in Fig. 1 includes three parts: a PMSG, a 0.69 kV/35 kV transformer, and a converter system. The converter system includes back-to-back fully-rated converter (FRC) with DC chopper branch, LC filter, and smoothing reactors. The generator-side filter is used to reduce the rate of change of the voltage and the peak voltage of the generator [

21], [22].

III. Equivalent Model of PMSG

A. Mathematical Model of PMSG

The mathematical equations of PMSG in abc frame is the background knowledge of this subsection [

23]. By transforming these equations with time-varying coefficients into d-q frame and assuming that the stator winding resistance is symmetric, the mathematical equations of PMSG are written as follows.

The voltage equations are:

UdUqU000=Rs00000Rs00000Rs00000Rkd00000RkqIdIqI0IkdIkq+ddtΨdΨqΨ0ΨkdΨkq+-ωrΨqωrΨd000 (1)

Equation (1) can be rewritten in the matrix form:

U=RI+dΨ/dt+e (2)

And the flux-linkage equations are:

ΨdΨqΨ0ΨkdΨkq=Ld00Lakd00Lq00Lakq00L000Lakd00Lkd00Lakq00LkqIdIqI0IkdIkq+Ψf00Ψf0LI+Ψm (3)

where ωr is the rotor angular velocity; U is the column vector of stator and rotor voltages; I is the column vector of current; Ψ is the column vector of flux-linkage; e is the column vector of speed voltage; R and L are the constant matrices of resistance and inductance, respectively; Ψm is the column vector of permanent magnet flux; the subscripts s and r represent stator and rotor, respectively; the subscripts d, q, and 0 represent the windings of d-, q-, and 0-axis, respectively; the subscripts kd and kq represent the damping windings of d- and q-axis, respectively; and the subscripts akd and akq represent the relationship between a- and d-axis and the relationship between a- and q-axis, respectively.

The rotor mechanical equations include the torque equation given in (4) and the speed equation given in (5).

Jdωrdt=Tm-Te-KDωr (4)
ωr=dθdt (5)
Te=ΨdIq-ΨqId (6)

where J is the rotational inertia; KD is the mechanical damping coefficient; Tm is the mechanical torque; Te is the electromagnetic torque; and θ is the rotor angle between a- and d-axis.

B. Equivalent Modeling of PMSG

The common flux-linkage Ψ in (2) and (3) is eliminated to:

dIdt=-L-1(R+X)I+L-1U-L-1F (7)
L-1=-Lkd/kd00Lakd/kd00-Lkq/kq00Lakq/kq001/L000Lakd/kd00-Ld/kd00Lakq/kq00-Lq/kq          Ld'00Lakd'00Lq'00Lakq'00L0'00Lakd'00Lkd'00Lakq'00Lkq'kd=Lakd2-LdLkdkq=Lakq2-LqLkqX=0-ωrLq00-ωrLakqωrLd00ωrLakd0000000000000000F=0ωrΨf000T (8)

where X is the asymmetric reactance matrix introduced from speed voltage column vector e; and F is the column vector related to the rotor angular velocity.

Discretize (7) with the trapezoidal rule (TR) integration method [

23] and obtain (9), and then partition (9) by separating the stator and rotor sides and obtain (10) in the matrix form.

1+ΔtRsLd'/2-ΔtωrLd'Lq/20ΔtRkdLakd'/2-ΔtωrLd'Lakq/2ΔtωrLq'Ld/21+ΔtRsLq'/20ΔtωrLq'Lakd/2ΔtRkqLakq'/2001+ΔtRsL0'/200ΔtRsLakd'/2-ΔtωrLakd'Lq/201+ΔtRkdLkd'/2-ΔtωrLakd'Lakq/2ΔtωrLakq'Ld/2ΔtRsLakq'/20ΔtωrLakq'Lakd/21+ΔtRkqLkq'/2IdtIqtI0tIkdtIkqt=ΔtLd'/200000ΔtLq'/200000ΔtL0'/200ΔtLakd'/200000ΔtLakq'/2000UdtUqtU0t00+Jdt-ΔtJqt-ΔtJ0t-ΔtJkdt-ΔtJkqt-ΔtJdq0(t-Δt)Jkdq(t-Δt) (9)
A(t)B(t)C(t)D(t)Idq0(t)Ikdq(t)=X110X210Udq0(t)0+Jdq0(t-Δt)Jkdq(t-Δt) (10)

where Jdq0(t-Δt) and Jkdq(t-Δt) represent the history current sources, which are given as follows.

Jdq0(t-Δt)Jkdq(t-Δt)=2E3-A(t)-B(t)C(t)2E2-D(t)Idq0(t-Δt)Ikdq(t-Δt)+X110X210Udq0(t-Δt)0+Wdq0Wkdq (11)
E2=1001E3=100010001Wdq0=0-Δt2ΨfLq'(ωr(t)+ωr(t-Δt))0Wkdq=0-Δt2ΨfLakq'(ωr(t)+ωr(t-Δt)) (12)

The stator and rotor currents in (13) and (14), respectively, are obtained by solving (10).

Idq0(t)=Gdq0(t)Udq0(t)+Idq0EQ(t-Δt) (13)
Ikdq(t)=Gkdq(t)Udq0(t)+IkdqEQ(t-Δt) (14)

where Gdq0(t) and Gkdq(t) are the admittance matrices of stator and rotor, respectively; and Idq0EQ(t-Δt) and IkdqEQ(t-Δt) are the history current sources of stator and rotor, respectively. These values are given in (15).

Gdq0(t)=(A(t)-B(t)D-1(t)C(t))-1(X11-B(t)D-1(t)X21)Gkdq(t)=(D(t)-C(t)A-1(t)B(t))-1(X21-C(t)A-1(t)X11)Idq0EQ(t-Δt)=(A(t)-B(t)D-1(t)C(t))-1(Jdq0(t-Δt)-B(t)D-1(t)Jkdq(t-Δt))IkdqEQ(t-Δt)=(D(t)-C(t)A-1(t)B(t))-1(Jkdq(t-Δt)-C(t)A-1(t)Jdq0(t-Δt)) (15)

For the purpose of connecting the PMSG model with external network, (13) is converted to abc frame via dq-abc transformation, thereby obtaining the equivalent equation (16).

iabc(t)=Gabc(t)uabc(t)+iabcEQ(t-Δt) (16)

where the subscript abc respesents the corresponding values in abc frame.

The expressions of stator admittance matrix Gabc(t) and history current source in abc frame IabcEQ(t-Δt) are given as:

Gabc(t)=P-1(t)Gdq0(t)P(t)IabcEQ(t-Δt)=P-1(t)Idq0EQ(t-Δt) (17)

where P(t) is the Park transformation matrix.

Similarly, the electromechanical transient model of PMSG is obtained by discretizing (4) and (5) with TR method, and then we can obtain (18) and (19).

ωr(t)=Kωr(t-Δt)+M(ΔT(t)+ΔT(t-Δt)) (18)
θ(t)=θ(t-Δt)+Δt2(ωr(t)+ωr(t-Δt)) (19)
K=(2J-KDΔt)/(2J+KDΔt)M=Δt/(2J+KDΔt)ΔT(t-Δt)=Tm(t-Δt)-Te(t-Δt)ΔT(t)=Tm(t)-Te(t) (20)

According to (16)-(19), the equivalent model of PMSG is developed in EMTDC program, as shown in Fig. 2. It includes EMT equivalent circuit and electromechanical transient model. The stator admittance matrix Gabc(t) is asymmetric because of the speed voltage vector e and matrix X, which are introduced from the frame transformation and flux-linkage Ψ elimination. Thus, it is not practical to develop the Norton equivalent circuit from (16). Therefore, the PMSG equivalent circuit is represented by three controlled current sources. When the wind speed changes, the mechanical torque will change according to (18), and the rotor speed and angle will change as well. In addition, the PMSG terminal voltage, current, and admittance matrix will change accordingly.

Fig. 2  Equivalent model of PMSG.

The rotor angular velocity ωr(t), rotor angle θ(t), and PMSG terminal voltage uabc(t) are needed to calculate parameters of PMSG equivalent circuit and used for abc-dq0 transformation in (16). However, the model only stores history values, so the values at time t need to be predicted. The linear extrapolation method given in (21) is used to predict the rotor angular velocity because of the large inertia of generator and slow change of rotational speed.

ωr(t)=ωr(t-Δt)+(ωr(t-Δt)-ωr(t-2Δt)) (21)

Then, the rotor angle at time t can be obtained from (19). Single-step approximation is applied to Udq0 with the assumption that Udq0(t) is close to Udq0(t-Δt). The rationality of this assumption is valid considering that the capacitor voltage of the LC filter at the outlet of PMSG will not change abruptly, and the voltage change can be ignored between microsecond-range time steps for simulating the switching behaviors of the power electronic devices in the wind farm. In addition, the voltage in the dq0 frame is the DC constant value, and the range of steady-state change is less than the AC value in the abc frame, so the approximate expression of PMSG terminal voltage can be obtained as:

uabc(t)=P-1(t)Udq0(t)P-1(t)Udq0(t-Δt) (22)

IV. Equivalent Model of Three-phase Transformer

A. Equivalent Model of Single-phase Transformer

An actual single-phase transformer considering copper loss is shown in Fig. 3(a). Figure 3(b) shows its T-type model, which includes an ideal transformer with turning radio N:1. In the T-type model, the main parameters of the actual model, e.g., the primary-side winding resistance R1, leakage reactance L1, excitation reactance Lm, secondary-side winding resistance R2, and leakage reactance L2 are preserved. The reference direction of the primary-side port voltage vT1, primary-side current iT1, secondary-side port voltage vT2, and secondary-side current iT2 are marked in Fig. 3(c).

Fig. 3  An actual single-phase transformer. (a) Single-phase transformer. (b) T-type model. (c) Equivalent model.

Apply Kirchohoff voltage law (KVL) to the T-type circuit in Fig. 3(b), and then we can obtain:

vT1(t)=R1iT1(t)+L1diT1(t)dt+LmddtiT1(t)+iT2(t)NvT2(t)=R2iT2(t)+L2diT2(t)dt+LmNddtiT1(t)+iT2(t)N (23)

Re-organize (23) and write it in the matrix form as:

ddtiT1(t)iT2(t)=L1+LmLm/NLm/NLm/N2+L2-1vT1(t)vT2(t)-R100R2iT1(t)iT2(t) (24)

Define

R=R100R2YT=Δt2L1+LmLm/NLm/NLm/N2+L2-1 (25)

Although the ideal transformer equations are often ill-conditioned, the TR method is proven to be adequate to solve the ill-conditioned equations [

24], [25].

Then, using TR method, (24) becomes:

iT1(t)iT2(t)=YMATvT1(t)vT2(t)+JT1(t-Δt)JT2(t-Δt) (26)
YMAT=(E2+YTR)-1YTy11y12y12y22YOUT=(E2+YTR)-1(E2-YTR)JT1(t-Δt)JT2(t-Δt)=YOUTiT1(t-Δt)iT2(t-Δt)+YMATvT1(t-Δt)vT2(t-Δt) (27)

where y11 and y22 both represent the self-admittances; y12 represents the mutual-admittance; and YOUT is a coefficient matrix. These parameters are constants because the transformer is a static device. JT1(t-Δt) and JT2(t-Δt) represent the transformer history current sources, which are obtained from the last time-step EMT simulation results.

According to (24)-(27), the equivalent model of single-phase transformer is derived, as shown in Fig. 3(c).

B. Equivalent Model of Three-phase Transformer

The typical three-phase transformer connection of YN/d11 is shown in Fig. 4(a).

Fig. 4  Model of three-phase transformer. (a) Transformer connection of YN/d11. (b) Phasor diagram. (c) Equivalent model of three-phase transformer. (d) Schematic YN/d11 of connection.

Figure 4(b) shows its phasor diagram. The delta connection of the secondary-side windings blocks the zero-sequence current and thus avoids the impact on AC system.

The equivalent model of the three-phase transformer, as shown in Fig. 4(c) is obtained by connecting three single-phase transformers following the YN/d11 connection. The schematic of YN/d11 connection is shown in Fig. 4(d).

V. Equivalent Model of Converter System

The structure of converter system is shown in Fig. 5(a), where the back-to-back FRC plays an isolated role between PMSG and the grid. Therefore, the AC frequency of PMSG can vary with wind speed variations and the grid frequency remains unchanged. The filter and smoothing reactor can reduce hormonic content and limit short-circuit current, respectively.

Fig. 5  Converter system. (a) Structure. (b) Companion circuit.

Figure 5(b) shows the companion circuit of the converter system, which is obtained by discretizing the capacitors and inductors in the converter and representing all the power electronic switches via two-state switch model (TSSM).

The discretization procedures of the capacitors and inductors using TR integration methods are not provided here because they are basic knowledge and have been very well documented in [

20]. Note that, to speed up the simulation, the blocked states of the FRC are not considered because the use of the diode-bridge-based blocked circuit will interrupt the node-elimination process. However, one or several detailed WT models without node-elimination can be used in the OWF model to provide the function of simulating the blocked state and dead zone effect of the FRC, while the rest of the WTs keep using the proposed integrated equivalent models.

VI. Integrated Modeling Method of PMSG-based WT

A. Complete Circuit of PMSG-based WT

By connecting all the individual devices, the complete circuit of each PMSG-based WT (excluding the mechanical drivetrain) is obtained as shown in Fig. 6. The electrical nodes are numbered from ① to 25, among which ① to ④ are external nodes; and ⑤ to 25 are internal nodes.

Fig. 6  Complete circuit of PMSG based WT.

B. Integrated Equivalent Modeling of PMSG-based WT

The nodal admittance equation of Fig. 6 is shown in Fig. 7, which is expressed in the form of (28).

Gu=J+i (28)

Fig. 7  Nodal admittance equation of Fig. 6.

Using the idea of partitioned matrix, the equation in Fig. 7 is written in a much more simplified form in (29) by separating the internal and external nodes.

Y11Y12Y21Y22uENuIN=0JIN+iEN0 (29)

where the subscript EN represents the external nodes; and IN represents the internal nodes.

By solving (29), the internal nodes are eliminated and nodal admittance equation of the external nodes is obtained as:

(Y11-Y12Y22-1Y21)uEN=iEN-Y12Y22-1JIN (30)
YEQ=Y11-Y12Y22-1Y21=YEQ11YEQ12YEQ13YEQ14YEQ21YEQ22YEQ23YEQ24YEQ31YEQ32YEQ33YEQ34YEQ42YEQ42YEQ43YEQ44YEQij=YEQji    i=1,2,,4; j=1,2,,4IEQ=-Y12Y22-1JIN=[IEQAIEQBIEQCIEQG]T (31)

where YEQ is the nodal admittance matrix of the equivalent model with reduced order; and IEQ is the injection current of external nodes.

The integrated equivalent model of each WT that contains only four nodes is shown in Fig. 8 according to (30) and (31). The equivalent model of the OWF is obtained by replacing all WTs in the OWF by the integrated equivalent model of each WT, as shown in Fig. 8.

Fig. 8  Integrated equivalent model of each WT.

The parameters of branch admittance in Fig. 8 are obtained as:

YEQA=j=14YEQ1j//YEQ14YEQB=j=14YEQ2j//YEQ24YEQC=j=14YEQ3j//YEQ34YEQAB=YEQ12YEQBC=YEQ23YEQAC=YEQ13 (32)

Although the internal node is eliminated, the information of the internal node can be preserved by (33), which is obtained by solving (29).

uIN=Y22-1(JIN-Y21uEN) (33)

The back-to-back FRC controller will use some of the information of internal node, so each simulation step needs to update the information of internal node, which inevitably requires additional computation. Computations increase with the number of internal nodes. The methods to reduce the computational burden of inverse solutions are given in the next subsection.

The proposed method only eliminates the circuit nodes of the WT without equivalent control. The access to each WT controller is preserved and the equivalent model only requires the controller to provide the triggering signals of the converter, which is independent of its control structure. Therefore, when the control structure needs to change, it is easy to adjust.

C. Preprocessing

To update the information of internal node, the admittance matrix inversion is needed at each time step according to (33). The matrix order increases with the number of internal nodes, which takes up a large proportion of the computational burden.

It can be observed that most of the elements in the nodal admittance matrix are constant due to the nature of static devices, and only the elements related to FRC are time-varying. Using this feature, the matrix order can be further reduced.

First, the admittance matrix of internal node is partitioned as:

Y22=A11A12A21A22 (34)

Then, the inverse matrix of (34) is:

Y22-1=A11-1(E13+A12F2A21A11-1)-A11-1A12F2-F2A21A11-1F2 (35)
F2=(A22-A21A11-1A12)-1 (36)

where E13 is a 13×13 identity matrix. The original 21×21 matrix inversion is replaced by an 8×8 matrix inversion. In addition, A11-1 is a constant matrix, which needs to be calculated only once at the beginning of simulation.

D. Modeling Procedures

Figure 9 shows the modeling procedures of the proposed method. Firstly, the EMT modeling methods of PMSG, AC transformer, back-to-back FRC, AC filter, and smoothing reactor in each WT are proposed and their companion circuits are obtained. Next, the complete original WT circuit is obtained by connecting all the above-mentioned companion circuits. Finally, the integrated equivalent modeling of each WT is carried out by eliminating all the internal nodes, thus obtaining the equivalent model which contains only four external nodes.

Fig. 9  Modelling procedures of the proposed method.

The detailed steps are as follows.

Step 1:   obtain the equivalent circuit of PMSG. First, predict the rotor angular velocity ωr(t) and rotor angle θ(t) via (19) and (21), respectively. Next, calculate the three-phase controlled current sources of PMSG iabc(t) through (16), (17), and (22).

Step 2:   obtain the equivalent circuit of three-phase double-winding transformer. First, calculate the admittance matrix YMAT and history current sources JT of single-phase transformer via (25) and (27), respectively. Next, calculate the parameters of three-phase transformer according to its connection mode.

Step 3:   obtain the companion circuit of the converter system. Discretize the capacitors and inductors in the FRC and filters of the converter system and represent the power electronic switches of the FRC through TSSM.

Step 4:   obtain the admittance matrix as shown in Fig. 7. First, connect all the circuits mentioned above and obtain the complete circuit of WT, as shown in Fig. 6. Next, write the nodal admittance equation, as shown in Fig. 7.

Step 5:   obtain the equivalent circuit of WT. First, eliminate the internal nodes of Fig. 6 through (30). Next, calculate the equivalent parameters YEQ and IEQ via (31) and (32), respectively.

Step 6:   obtain the EMTDC solution. Then, obtain the voltage uEN of the four reserved external nodes.

Step 7:   update the information of internal nodes. Use the above-mentioned external node voltages uEN and obtain the information of internal nodes through (33).

Step 8:   update the parameters related to PMSG. Using the information of internal nodes 12-14, the rotor and stator current Idq0(t) and Ikdq(t), electromagnetic torque Te(t), rotor angular velocity ωr(t) and rotor angle θ(t) are updated via (6), (13), (14), (18), and (19).

Step 9:   update the parameters related to transformer. Using the internal nodes ⑤-⑧ and external nodes ①-④, the transformer port currents iT1(t) and iT2(t) are updated through (26).

Step 10:   update all the currents of the capacitor and inductor branches in the converter system using internal nodes 6-25.

The calculation procedure at time t is completed, and by repeating Steps 1-10, computations for the next integration step at time t+Δt can be pursued.

VII. Model Validation

This section validates the accuracy and efficiency of the developed equivalent model, which is simplified as EM in this section, by comparison with the EMT time-domain simulation of a fully detailed model (DM) of the OWF using PSCAD/EMTDC.

A. Configuration and Control Strategy of Test System

The test system structure of a typical OWF is shown in Fig. 1 and the system parameters are given in Appendix A. The FRC is modulated by pulse width modulation (PWM). The grid-side control regulates the DC bus voltage and the reactive power, while the machine-side controller is used to control the real power and the AC voltage at the terminal of each PMSG. This control strategy is applied in [

26] and the user’s library of PSCAD/EMTDC [27]. The DC chopper circuit is used to protect the DC bus from overvoltage.

B. Performance of EM of Single WT Under Varying Wind Speeds

The EM of WT is used to check the accuracy by simulating conditions with varying wind speeds. Figure 10 shows the simulation results when wind speed changes from 6 m/s to 15 m/s. The system reaches the stable state at t=1 s with wind speed 10 m/s. The wind speed Vw decreases to 6 m/s at t=7 s and then ramps up from 6 m/s to 15 m/s at t=10 s. At t=17 s, the wind speed decreases from 15 m/s to 10 m/s, and then maintains at 10 m/s until the end of simulation.

Fig. 10  Simulation results for varying wind speeds. (a) Wind speed. (b) Real power at PCC. (c) AC voltage at PCC. (d) DC voltage VDC. (e) Rotor speed.

The results show that the EM is visually overlapping with DM, and the maximum relative errors are below 2.16%. The real power and the voltage at PCC PPCC and URMSPCC in Fig. 10 confirm that the individual WT model is accurate when reproducing the external behaviors. Figure 10(d) and (e) shows the features of internal nodes which confirm that the method can preserve the internal-machine information.

C. Performance of EM of Single WT in Case of Oscillation

In order to stimulate the sub-synchronous interactions between the WTs and the AC system, the regulator gain of real power outer loop steps from 0.25 to 5 at t=1 s. The legend TSSM represents the power electronic switches in DM are represented by TSSM. The simulation comparisons among DM, TSSM, and EM are shown in Fig. 11.

Fig. 11  Simulation comparisons among DM, TSSM, and EM. (a) Real power at PCC. (b) Instantaneous current at PCC.

Figure 11 shows the real power and instantaneous current at PCC, i.e., PPCC and IPCC, respectively. It can be observed that, as the regulator gain of real power outer loop steps from 0.25 to 5, the system gradually diverges and appears a sub-synchronous oscillation around 25 Hz. The oscillation frequency of DM is 25.3 Hz while the oscillation frequency of TSSM and EM is 24 Hz. The relative error of oscillation frequency is 5.14%.

The waveforms obtained from the TSSM and EM are visually overlapping, yet there is a slight phase difference between them and DM, which is because the TSSM and EM do not enjoy the interpolation feature as that of DM in the EMT program.

D. Performance of Proposed OWF Model with Different Wind Speeds Following a Three-phase-to-ground Fault

In order to further validate the proposed OWF model, i.e., EM, a typical OWF consisting of 20 WTs, as shown in Fig. 1, is simulated. Both the DM and EM are developed and compared. The WTs have different wind speeds as given in Table II.

TABLE II  Wind Speeds of WTs
WT No.

Vw

(m/s)

WT No.

Vw

(m/s)

WT No.

Vw

(m/s)

WT No.

Vw

(m/s)

1 12.6 6 13.1 11 12.8 16 12.5
2 11.3 7 11.6 12 11.7 17 11.5
3 10.5 8 10.1 13 9.8 18 10.4
4 9.5 9 9.2 14 9.0 19 9.3
5 8.2 10 8.0 15 7.4 20 7.8

Figure 12 shows the real power, voltage, and current at PCC following a temporary 50 ms three-phase-to-ground fault.

Fig. 12  Comparison results under short-circuit fault. (a) Voltage at PCC. (b) Real power at PCC. (c) Instantaneous voltage at PCC. (d) Instantaneous current at PCC.

Figure 12 indicates that the EM essentially maintains identical transient response as the DM. During the fault, the voltage at PCC uPCC drops to 0.2 p.u. of the rated voltage, and the voltage at PCC can recover to 0.9 p.u. of the rated voltage within 2 s after the voltage drop.

Figures 13-15 show the real and reactive power, voltages and AC currents, and DC voltages of the No. 1 and No. 20 WTs in the OWF. As can be observed, the simulation curves are overlapping with less than 4% maximum relative error, indicating that the information of internal node is retained. The computation time is tested in next subsection. Figure 15 illustrates that EM can properly reproduce the internal switching action procedure of the system.

Fig. 13  Comparison results of real and reactive power. (a) Real power on grid-side converter of the No. 1 WT (PG1). (b) Reactive power on grid-side converter of the No. 1 WT (QG1). (c) Real power on grid-side converter of the No. 20 WT (PG20). (d) Reactive power on grid-side converter of the No. 20 WT (QG20).

Fig. 14  Comparison results of voltages and AC currents. (a) Voltage on grid-side converter of the No. 1 WT (UGS1). (b) AC current on grid-side converter of the No. 1 WT (IGS1). (c) Voltage on grid-side converter of the No. 20 WT (UGS20). (d) AC current on grid-side converter of the No. 20 WT (IGS20).

Fig. 15  Comparison results of DC voltages. (a) DC voltage of the No. 1 WT (VDC1). (b) DC voltage of the No. 20 WT (VDC20).

E. Computation Efficiency Test

As each detailed model of WT is integrated to an equivalent model containing only four nodes, there should be a significant reduction in the simulation time (defined as CPU time in this paper). Transmission lines are not considered in all models in this subsection because the consideration of lines would result in unacceptable simulation time for DM. One or several DMs without node-elimination can be used in the OWF model to provide the function of simulating the internal fault, while the rest of the WTs keep using the EM (this model is called EDM). In order to show the impact of internal fault simulation on CPU time, the EDM is also compared in which an EM is replaced by a DM. Table III gives the CPU time for an OWF with 1, 5, 10, 20, and 50 WTs on a PC with an Intel Core i7-6700HQ CPU running at 2.60 GHz. All the simulation models are implemented on PSCAD/EMTDC Professional V4.6. The time step of the EMT solver is set to be 10 μs and the duration is 1 s. Speedup factor 1 is the ratio of DM and EM, and speedup factor 2 is the ratio of DM and EDM.

TABLE III  Comparisons of Simulation Time
Number of WTsCPU time (s)

Speedup

factor 1

Speedup

factor 2

DMEMEDM
1 18.75 4.60 18.75 4.08 1.00
5 388.33 13.60 22.91 28.56 16.95
10 1665.30 26.60 36.36 62.62 45.80
20 10950.33 54.36 73.55 201.44 148.89
50 58055.11 189.69 208.84 306.06 277.98

For an OWF with 50 WTs, the proposed EM gives a speedup factor 1 of 306.06, which is over two orders of magnitude faster than the DM, while the EM has only 4% relative error as compared to DM.

For an EDM, the CPU time will increase, but the OWF model still has a significant speedup factor. When the included number of WTs in the OWF becomes larger, the impact of internal fault simulation on the simulation speed due to using EDM will be less.

The CPU time as well as the speedup factors in Table III are also graphically shown in Fig. 16.

Fig. 16  Comparisons of CPU time and speedup factor.

An interesting observation is that the CPU time of EM is almost linear while that of the DM rises more drastically. It could be expected that with hundreds of WTs, the time savings of the proposed OWF models will be significantly larger than the DM.

VIII. Conclusion

This paper attempts to speed up the simulation of a large-scale OWF using the proposed integrated model of the WT. The model is developed by establishing the EMT circuit model of each device and eliminating all the internal nodes of the complete circuit. The features of the model are as follows.

1) The proposed model is sufficiently accurate and efficient compared with the fully detailed model of OWF. When the number of WTs increases, the acceleration is more significant. As a result, the proposed model is suitable for large-scale OWF simulation. But it still works for small-scale OWF. Besides, it provides significant speedup factors for the off-line EMT simulation, and can save hardware and computation cores to a large extent for the real-time EMT simulation.

2) The developed model preserves the characteristics of all internal nodes and the access to external circuits and controls. The control structure can be adjusted freely to meet different research demands. Although the method only considers the PMSG-based WT, the modeling procedures can be used as reference to the simulation of other types of WTs once the structures and parameters are given. The developed model cannot be directly used to simulate other types of WTs.

3) The internal fault can be achieved by replacing an EM in the OWF model as the DM. It will slightly affect the simulation time of EM, but the OWF model still has a good speedup factor.

Appendix

Appendix A

TABLE AI  System Parameters of OWF
DeviceSymbolParameter DescriptionValueDeviceSymbolParameter descriptionValue
Wind wheel Rad,turb (m) Turbine radius 50 Filter Rd (Ω) Resistance of filter 1.332
Vbwind (m/s) Rated wind speed 10 Lf (H) Inductance of filter 0.00062
Vcutin (m/s) Cut-in wind speed 3 Cf (μF) Capacitance of filter 700
Vcutout (m/s) Cut-out wind speed 25 Cd (μF) Capacitance of filter 700
PMSG SbPMSG (MVA) Rated capacity of PMSG 2 Transformer SbT (MVA) Rated capacity of transformer 2
UbPMSG (kV) Rated voltage of PMSG 0.69 vbT1 (kV) Primary-side voltage rating of transformer 35
fbPMSG (Hz) Base frequency of PMSG 25 vbT2 (kV) Secondary-side voltage rating of transformer 0.69
Rs (p.u.) Stator winding resistance 0.0017 RT (p.u.) Resistance of transformer 0.00042
Xs (p.u.) Stator leakage reactance 0.0364 LT (p.u.) Leakage inductance of transformer 0.025
Xd, Xq (p.u.) Reactances of d and q axes 0.55, 1.11

35 kV

cable

Rc1 (Ω/km) Cable resistance 0.128
Rad (p.u.) Damping winding resistance of d axis 0.055 Lc1 (H/km) Cable inductance 0.00102
Xad (p.u.) Damping winding reactance of d axis 0.62 220 kV cable Rc2 (Ω/km) Cable resistance 0.01288
Raq (p.u.) Damping winding resistance of q axis 0.183 Lc2 (H/km) Cable inductance 0.001
Xaq (p.u.) Damping winding reactance of q axis 1.175 Cc2 (μF/km) Cable capacitance 0.20559
Ψf (p.u.) Permanent magnet flux 1 Unit step-up transformer SNT (MVA) Rated capacity of transformer 200
J (s) Rotational inertia 4 vNT1 (kV) Primary-side voltage rating of transformer 220
KD (p.u.) Mechanical damping coefficient 0.01 vNT2 (kV) Secondary-side voltage rating of transformer 35
Back-to-back FRC SbFRC (MVA) Rated capacity of FRC 2

AC

system

SbAC (MVA) Rated capacity of AC system 200
UbDC (kV) Rated DC voltage 1.45 UbAC (kV) Rated voltage of AC system 220
C (μF) DC-side capacitance 15000 fbAC (Hz) Base frequency of AC system 50
fs (Hz) Switching frequency 3800 RbAC (Ω) Internal resistance of AC system 25

References

1

P. Verma, S. K. and B. Dwivedi, “A cooperative approach of frequency regulation through virtual inertia control and enhancement of low voltage ride-through in DFIG-based wind farm,” Journal of Modern Power Systems and Clean Energy, vol. 10, no. 6, pp. 1519-1530, Nov. 2022. [Baidu Scholar] 

2

M. Wang, Y. Mu, H. Jia et al., “Active power regulation for large-scale wind farms through an efficient power plant model of electric vehicles,” Applied Energy, vol. 185, no. 2, pp. 1673-1683, Jan. 2017. [Baidu Scholar] 

3

U. Karaagac, J. Mahseredjian, R. Gagnon et al., “A generic EMT-type model for wind parks with permanent magnet synchronous generator full size converter wind turbines,” IEEE Power and Energy Technology Systems Journal, vol. 6, no. 3, pp. 131-141, Sept. 2019. [Baidu Scholar] 

4

J. Ding, K. Xie, B. Hu et al., “Mixed Aleatory-epistemic Uncertainty Modeling of Wind Power Forecast Errors in Operation Reliability Evaluation of Power Systems,” Journal of Modern Power Systems and Clean Energy, vol. 10, no. 5, pp. 1174-1183, Sept. 2022. [Baidu Scholar] 

5

D. J. Trudnowski, A. Gentile, J. M. Khan et al., “Fixed-speed wind-generator and wind-park modeling for transient stability studies,” IEEE Transactions on Power Systems, vol. 19, no. 4, pp. 1911-1917, Nov. 2004. [Baidu Scholar] 

6

M. Ali, I. Ilie, J. V. Milanovic et al., “Wind farm model aggregation using probabilistic clustering,” IEEE Transactions on Power Systems, vol. 28, no. 1, pp. 309-316, Feb. 2013. [Baidu Scholar] 

7

W. Li, A. M. Gole, M. K. Das et al., “Research on wind farms aggregation method for electromagnetic simulation based on FDNE,” in Proceedings of 2019 IEEE PES ISGT-Europe, Bucharest, Romania, Sept. 2019, pp. 1-5. [Baidu Scholar] 

8

L. M. Fernandez, C. A. Garcia, F. Jurado et al., “Aggregation of doubly fed induction generators wind turbines under different incoming wind speeds,” in Proceedings of 2005 IEEE Russia Power Tech, St. Petersburg, Russia, Jun. 2005, pp. 1-6. [Baidu Scholar] 

9

M. A. Chowdhury, W. X. Shen, N. Hosseinzadeh et al., “A novel aggregated DFIG wind farm model using mechanical torque compensating factor,” Energy Conversion and Management, vol. 67, pp. 265-274, Mar. 2013. [Baidu Scholar] 

10

J. Zhang, F. Miao, and R. Gu, “Research on equivalent electromechanical transient modeling of PMSG-based wind farms,” in Proceedings of 2020 12th IEEE PES Asia-Pacific Power and Energy Engineering Conference (APPEEC), Nanjing, China, Sept. 2020, pp. 1-5. [Baidu Scholar] 

11

P. Zhang, J. R. Marti, and H. W. Dommel, “Shifted-frequency analysis for EMTP simulation of power-system dynamics,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 57, no. 9, pp. 2564-2574, Sept. 2010. [Baidu Scholar] 

12

Y. Li, D. Shu, F. Shi et al., “A multi-rate co-simulation of combined phasor-domain and time-domain models for large-scale wind farms,” IEEE Transactions on Energy Conversion, vol. 35, no. 1, pp. 324-335, Mar. 2020. [Baidu Scholar] 

13

N. Shabanikia, A. A. Nia, A. Tabesh et al., “Weighted dynamic aggregation modeling of induction machine-based wind farms,” IEEE Transactions on Sustainable Energy, vol. 12, no. 3, pp. 1604-1614, Jul. 2021. [Baidu Scholar] 

14

W. Teng, X. Wang, Y. Meng et al., “An improved support vector clustering approach to dynamic aggregation of large wind farms,” CSEE Journal of Power and Energy Systems, vol. 5, no. 2, pp. 215-223, Jun. 2019. [Baidu Scholar] 

15

B. Zhang and Z. Zhang, “A two-stage model for asynchronously scheduling offshore wind farm maintenance tasks and power productions,” International Journal of Electrical Power & Energy Systems, vol. 130, pp. 1-10, Sept. 2021. [Baidu Scholar] 

16

N. Lin and V. Dinavahi, “Exact nonlinear micromodeling for fine-grained parallel EMT simulation of MTDC grid interaction with wind farm,” IEEE Transactions on Industrial Electronics, vol. 66, no. 8, pp. 6427-6436, Aug. 2019. [Baidu Scholar] 

17

D. Shu, X. Xie, Q. Jiang et al., “A novel interfacing technique for distributed hybrid simulations combining EMT and transient stability models,” IEEE Transactions on Power Delivery, vol. 33, no. 1, pp. 130-140, Feb. 2018. [Baidu Scholar] 

18

T. Duan and V. Dinavahi, “Variable time-stepping parallel electromagnetic transient simulation of hybrid AC-DC grids,” IEEE Journal of Emerging and Selected Topics in Industrial Electronics, vol. 2, no. 1, pp. 90-98, Jan. 2021. [Baidu Scholar] 

19

V. Jalili-Marandi, L. Pak, and V. Dinavahi, “Real-time simulation of grid-connected wind farms using physical aggregation,” IEEE Transactions on Industrial Electronics, vol. 57, no. 9, pp. 3010-3021, Sept. 2010. [Baidu Scholar] 

20

U. N. Gnanarathna, A. M. Gole, and R. P. Jayasinghe, “Efficient modeling of modular multilevel HVDC converters (MMC) on electromagnetic transient simulation programs,” IEEE Transactions on Power Delivery, vol. 26, no. 1, pp. 316-324, Jan. 2011. [Baidu Scholar] 

21

Technical Specification for Converter Filter of Wind Turbine Generator, the People’s Republic of China Energy Industry Standard, NB/T 10437-2020. [Baidu Scholar] 

22

Wind Turbines Generator System—Full-power Converter—Part 1: Technical Condition, the People’s Republic of China National Standard, GB/T 25387.1-2021. [Baidu Scholar] 

23

P. Kundur, Power System Stability and Control. New York: McGraw-Hill, 1994, pp. 45-138. [Baidu Scholar] 

24

R. L. Burden and J. D. Faires, “Ordinary differential equations,” in Numerical Analysis, 9th ed. Boston, USA: Richard Stratton, 2011, pp. 348-352. [Baidu Scholar] 

25

H. K. Lauw and W. S. Meyer, “Universal machine modeling for the representation of rotating electric machinery in an electromagnetic transients program,” IEEE Power Engineering Review, vol. PER-2, no. 6, pp. 24-25, Jun. 1982. [Baidu Scholar] 

26

S. Dong, H. Li, and Y. Wang. “Low voltage ride through capability enhancement of PMSG-based wind turbine,” in Proceedings of International Conference on Sustainable Power Generation and Supply (SUPERGEN 2012), Hangzhou, China, Sept. 2012, pp. 1-5. [Baidu Scholar] 

27

K. Clark, N. W. Miller, and J. J. Sanchez-Gasca. “Modeling of GE wind turbine-generators for grid studies,” General Electric Co., Tech. Rep. Apr. 2010. [Baidu Scholar]