Skip to content

Commit

Permalink
EMTP-RV Models
Browse files Browse the repository at this point in the history
  • Loading branch information
@lvanfretti
lvanfretti committed Jul 25, 2017
1 parent 8d871e8 commit 0b86730
Show file tree
Hide file tree
Showing 7 changed files with 1,007 additions and 0 deletions.
201 changes: 201 additions & 0 deletions EMTP-RV/ADN EPEC2017/DFIG_mask.dwj
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,201 @@
// Last change: HANI SAAD - JEAN MAHSEREDJIAN 04/11/2013 03:36:19 AM
//*This script is used for the DFIG wind generator model.
// This scripts allows to modify parameters that are not accessible on the main mask.

with (Math) {

//////////////////////////////////////////////////_BASIC CALCULATIONS //////////////////////////////////////////////////
//*Fundamental Period
Period = 1 / Freq;
period = Period;
w = 2 * PI * Freq;
//*Total Wind park S and P
Pgen_W = Pgen_MW * 1e6; // Rated active power of wind generator (W)
Sgen_VA = Pgen_W * 1.11; // Rated apparent power of wind generator (VA). Note that 11% of reserve power is added.
Spark_VA = Ngen * Sgen_VA; // Total aggregated wind park S (VA)
Spark_MVA = Ngen * Sgen_VA / 1e6; // Total aggregated wind park S (MVA), ASM machine uses MVA selection
Ppark_W = Ngen * Pgen_W; // Total aggregated wind park P (W) (used for the load-flow device)
Qpark_Var = Qref_pu * Spark_VA; // Total aggregated wind park Q (Var) (used for the load-flow device when Q or Q+Vac control is selected)
Vdc_V = Vdc_kV * 1e3; // DC voltage (V)
//*base calculations
Vbase_VRMSLL = Vgen_kVRMSLL * 1e3; // Rated generator voltage (VRMSLL)
Vbase_V = Vbase_VRMSLL * sqrt(2 / 3); // Peak generator voltage (V)
Ibase_A = Spark_VA / Vbase_VRMSLL * sqrt(2 / 3); // Peak generator current (A)
Zbase_Ohm = (Vgen_kVRMSLL * Vgen_kVRMSLL) / Spark_MVA; // Base impedance (Ohm)

//////////////////////////////////////////////////_ASYNCHRONOUS GENERATOR & AC/DC/AC CONVERTER PARAMETERS //////////////////////////////////////////////////
//*ASM parameters (see ASM_DFIG/ASM1)
npoles = 3; // KTH, original value is 6 Number of poles of the asynchronous generator
Rs_pu = 0.00706 // KTH, original value is 0.007 in pu
Lls_pu = 0.171; // in pu
Lmd_pu = 2.9 // in pu
Lmq_pu = 2.9 // in pu
Rr_pu = 0.005 // in pu
Llr_pu = 0.156 // in pu
Hgen_pu = 5.06 // KTH=5.04 in sec
D_pu = 0.0 // in pu
Llr_plus_Lmd_pu = Llr_pu + Lmd_pu; // in pu (used in DFIG_CONTROL/ROTOR side CONTROL)
//*Choke impedance (see ASM_DFIG/Choke)
Rchoke_pu = 0.003; // in pu
Rchoke_Ohm = 1.8034e-5; // KTH, original value is "Rchoke_pu * Zbase_Ohm " in Ohm
Lchoke_pu = 0.3 // in pu
Lchoke_H = 5.7404e-6; // KTH, original value is "Lchoke_pu * Zbase_Ohm / (2 * PI * Freq)" in H
//*AC/DC/AC Converter parameters (see ASM_DFIG/BtoB_converter)
E_Cdc = 5; // in KJ/MVA
Cdc_F = 0.01; // KTH=0.01, original value is "E_Cdc * 1e-3 / (1 / 2 * pow(Vdc_V, 2) / Sgen_VA); " DC capacitor (F)
CdcPark_F = Cdc_F * Ngen; // Aggregated DC capacitor (F)
CarrierSignal_Freq = 3e3; // PWM Carrier signal frequency (Hz) (used only for Detailed Converter Model)
CarrierSignal_ratio = CarrierSignal_Freq / Freq; // PWM carrier signal ratio (used only for Detailed Converter Model)
Rsnubber_Ohm = 3000; // IGBT/diode resistance snubber (Ohm) (used only for Detailed Converter Model)
Csnubber_F = 1e-6; // IGBT/diode capacitor snubber (F) (used only for Detailed Converter Model)
Rchopper_Ohm = 0.2; // Chopper resistance in Ohm
RchopperPark_Ohm = Rchopper_Ohm / Ngen; // Aggregation of the chopper resistance
//*High pass filter parameters (see ASM_DFIG/filter_stator_ASM )
Qfilter_Var = 0.473e6 // KTH setting, original value is "150e3 * Ngen" Reactive power generated by each filter
n1 = CarrierSignal_ratio // First tuning harmonic
n2 = CarrierSignal_ratio * 2 // Second tuning harmonic
QF = 1e3 // Quality factor
Cfilter1_F = Qfilter_Var / (Vbase_VRMSLL * Vbase_VRMSLL * 2 * PI * Freq);
Lfilter1_H = 1e6; // KTH, original value is " 1 / (Cfilter1_F * pow((n1 * Freq * 2 * PI), 2));"
Rfilter1_Ohm = 1.7 // KTH, original value is "QF * 2 * PI * n1 * Freq * Lfilter1_H;"
Cfilter2_F = Cfilter1_F;
Lfilter2_H = 1 / (Cfilter2_F * pow((n2 * Freq * 2 * PI), 2));
Rfilter2_Ohm = QF * 2 * PI * n2 * Freq * Lfilter2_H;

//////////////////////////////////////////////////_WIND TURBINE PARAMETERS //////////////////////////////////////////////////
//* Blade and Gearbox parameters (see WIND_TURBINE/Wind_Power)
blade_radius = ceil((sqrt(pow(35.25, 2) / 1.5e6 * Pgen_W)) * 100) / 100; // Blade radius (m) is in function of the rated power
gearbox_ratio = 72 * (Freq / 60); // Gearbox ratio. Note that the ratio is adjusted depending on the frequency used in order to preserve the same tip-speed value (Lambda).
ro = 1.225; // air density (kg/m^3)
K_speed = 2 * PI * Freq / (npoles / 2) * blade_radius / gearbox_ratio; // speed coefficient defined as: Lambda = K_speed / wind_speed * w_rotor;
K_air = (1 / 2) * PI * blade_radius * blade_radius * ro; // air coefficient defined as: Pair = K_air * Cp * (wind_speed)^3;
Pitch_max_deg = 45; // KTH, original value is 20 Maximum pitch angle in degree
Pitch_min_deg = 0; // KTH, original value is -1.5 Minimum pitch angle in degree

//*Shaft and train drive parameters (see WIND_TURBINE/Two_mass_model)
Ksh = 1.2; //Speed deviation damping
Dm = 1.5; //Mutual damping
Ht = 4 // Inertia constant (s)

//////////////////////////////////////////////////_TRANSFORMER PARAMETERS //////////////////////////////////////////////////
//* Wind generator transformer parameters (see Transformer)
Rtransfo_sec_Ohm = (Rtransfo_pu * 0.1 * Vgrid_kVRMSLL * Vgrid_kVRMSLL / Spark_MVA); // Secondary winding resistance, Ohm
Ltransfo_sec_H = (Xtransfo_pu * 0.1 * Vgrid_kVRMSLL * Vgrid_kVRMSLL / Spark_MVA) / (2 * PI * Freq); // Secondary winding inductance, H
Rtransfo_prim_Ohm = (Rtransfo_pu * 0.9 * Vgen_kVRMSLL * Vgen_kVRMSLL / Spark_MVA); // Primary winding resistance, Ohm
Ltransfo_prim_H = (Xtransfo_pu * 0.9 * Vgen_kVRMSLL * Vgen_kVRMSLL / Spark_MVA) / (2 * PI * Freq); // Primary winding inductance, H
Transfo_ratio = Vgen_kVRMSLL / Vgrid_kVRMSLL / sqrt(3); // Turns ratio

//////////////////////////////////////////////////_CONTROL SYSTEM PARAMETERS //////////////////////////////////////////////////
//* Pitch control (see CONTROL/Pitch_Control )
// NB: Ptich Control Gains are tuned by trial and error
PitchCtrl_ki = 20; // Integral Gain of the Pitch control
PitchCtrl_kp = 500; // KTH, original value is 50 Proportional Gain of the Pitch control
PitchCtrl_Td_d = 10; // Derivative Gain of the Pitch control
PitchCtrl_tau_d = 5e-3; // derivative filter in s
w_rotor_max_pu = 1.25; // maximum rotor speed in pu
Pitch_rate = 2 // KTH, original value is 10 Slop rate (see Pitch_Control/Rate_Limiter)

//*Rotor control (see CONTROL/Rotor_Control )
//+Rotor Inner current control
eps = 0.7; // Damping ratio
w_i = 3 / (eps * TimeConstant_i_rotor); // Impulse of the closed loop control
RotorCtrl_ki = 8 // KTH, 0.2 * ((Llr_plus_Lmd_pu) / w) * w_i * w_i; Integral Gain of the Inner Current control
RotorCtrl_kp = 0.3; // KTH, original value is 0.2 * 2 * eps * w_i * ((Llr_plus_Lmd_pu) / w); Proportional gain of the Inner Current control
//+Q control and Q+Vac control
Qctrl_kp = 0.05; //KTH, original value is 0 Proportional Gain of Q control (Only Integral part of the PI control is used)
Qctrl_ki = 5; // KTH, original value is "3 / TimeConstant_Q" Integral Gain of Q control
//+Vac control
X_tot = X_PCC_pu + Xtransfo_pu; // Total impedance in pu at wind generator connection point
VacCtrl_ki = 3 / (X_tot * TimeConstant_Vac); // Integral Gain of Vac control
VacCtrl_kp = 1; // Proportional Gain of Vac control (Only Integral part of the PI control is used)
//+P control
Pctrl_kp = 1; // Proportional Gain of P control (Only Integral part of the PI control is used)
Pctrl_ki = 3 / TimeConstant_P; // Integral Gain of P control
//+Maximum Power Point Tracking (MPPT)
Lambda_optimal = 6.635; // Optimal Lambda in order to have the rated active power for a wind speed = 11.24 m/s
Cp_max = 0.46; // Related Cp_max in order to have the rated active power for a wind speed = 11.24 m/s
K_optimal_pu = K_air * Cp_max * pow(K_speed / Lambda_optimal, 3) / Sgen_VA; // Constant coefficient defined as: P = K_optimal_pu * w^3
//+ Idq limit (see CONTROL/Rotor_Control/Idq_ref_limiter1 )
Q_priority = 1; // 1 = priority is given for reactive power (should be set for FRT capability). 0 = priority is given to active power
I_lim_Rotor_pu = 1.1; // Maximum rotor Idq limit
Id_lim_Rotor_pu = 1; // Maximum rotor Id limit, if Q_priority = 1
Iq_lim_Rotor_pu = 1.1; // Maximum rotor Iq limit, if Q_priority = 0

//*Grid control (see CONTROL/Grid_Control )
//+Grid Inner current control
eps = 0.7; // Damping ratio
w_i = 3 / (eps * TimeConstant_i_grid); // Impulse of the closed loop control
GridCtrl_ki = 500; // KTH, original value is "(Lchoke_pu / w) * w_i * w_i" Integral Gain of the Inner Current control
GridCtrl_kp = 2.5; // KTH, original value is "2 * eps * w_i * (Lchoke_pu / w)" Proportional gain of the Inner Current control
//+Vdc control
E_Cdc_J = 0.5 * CdcPark_F * Vdc_V * Vdc_V; // Energy stored in the DC capacitor (Joule)
H_Cdc = E_Cdc_J / (0.3 * Spark_VA); // Static moment of inertia (S). Note: Power converter is rated for 30% of the nominal power
eps = 1; // Damping ratio
w_Vdc = 3 / (eps * TimeConstant_Vdc); // Time response for 5% (sec)
VdcCtrl_ki = H_Cdc * w_Vdc * w_Vdc; // KTH=0.05, Proportional Gain of Vdc control
VdcCtrl_kp = 2 * eps * w_Vdc * H_Cdc; // KTH=2e-3 Integral Gain of Vdc control

//*PLL (see CONTROL/Computing_Variables/PLL)
PLL_lim_low = 12
PLL_lim_up = 12
PLL_max_delay = 0.0223
eps = 1; // damping ratio
w_i = 3 / (eps * TimeConstant_PLL);
PLL_ki = 1 * w_i * w_i; // Int. Gain integral (rad/V)
PLL_kp = 2 * eps * w_i * 1; // Prop. Gain (rad/s/V)

//*Protection system (see CONTROL/Protection )
//+Deep Voltage Sag (DVS) protection
DVS_level = 0.001; // Threshold Voltage value in pu to activate the DVS protection
DVS_hysteresis = 0.1; // Threshold Voltage value in pu to deactivate the DVS protection
//+Crowbar protection
Rcrowbar_Ohm = 0.09; // Crowbar resistance (Ohm)
Lcrowbar_H = 4.0e-05; // Crowbar inductance (H)
RcrowbarPark_Ohm = Rcrowbar_Ohm / Ngen; // Crowbar equivalent resistance (Ohm)
LcrowbarPark_H = Lcrowbar_H / Ngen; // Crowbar equivalent inductance (H)
//+Converter Overcurrent limit
I_RotorSideConv_max = 4; // Threshold Current value in pu to activate RotorSide Overcurrent protection
I_GridSideConv_max = 4; // Threshold Current value in pu to activate GridSide Overcurrent protection

//////////////////////////////////////////////////_INITIALIZATION //////////////////////////////////////////////////
//*Initialization of the Protection system (see CONTROL/Protection)
init_Prot_delay = 0.3; // Delay time initialization of the Protection devices (s) (All protections, expect the Chopper, are deactivated for the specified time )

//*Initialization of the WIND TURBINE model (see WIND_TURBINE/Two_mass_model)
Lambda_ini_guess = Lambda_optimal; // Initial Lambda guess value
Cp_init_guess = Cp_max; // Initial Cp guess value
max_wind_speed = K_speed * w_rotor_max_pu / Lambda_ini_guess; // maximum wind speed before the Pitch control is activated
if (Mean_wind_speed > max_wind_speed) { // if the wind speed value specified by the user is higher than the maximum wind speed before Pitch control is activated
Mean_wind_speed_init = max_wind_speed; // than the wind speed value is fixed at his maximum value for initialization.
}
else {
Mean_wind_speed_init = Mean_wind_speed;
}

init_w_turbine_pu = Mean_wind_speed_init * Lambda_ini_guess / K_speed * 0.989; // Initial Turbine speed in pu. The value is multiplied by a factor of losses.
init_Pturbine_pu = K_air * Cp_init_guess * pow(Mean_wind_speed_init, 3) / Sgen_VA;
init_Tturbine_pu = init_Pturbine_pu / init_w_turbine_pu; // Initial Turbine Torque in pu.

//*Initialization of _ASYNCHRONOUS GENERATOR model (See ASM_DFIG/ASM1)
init_slip_pct = (1 - init_w_turbine_pu / 1) * 100; // Initial slip (%)

/////////////////////////////////////////////////////////////_LOAD-FLOW CALCULATIONS //////////////////////////////////////////////////

//* Load-Flow of the wind generators (see LF_init1)
LF_V = Vgen_kVRMSLL * 1000 * sqrt(2 / 3); // Voltage for the LF device at the bus of the generator

if (ControlMode == 1 || ControlMode == 3) { // if Q control or Q+Vac Control is selected, than PQ constraint is activated
LF_t1 = -1; // switch trick to active PQ constraint
LF_t2 = 1e15; // switch trick to active PQ constraint
LF_Ppark1_W = init_Pturbine_pu * Spark_VA * 0.967; // Active power in W of PQ constraint. The value is multiplied by a factor of losses.
LF_Ppark2_W = 0; // Active power in W of PV constraint. The value is set to zero
} else { // if Vac control is selected, than PV constraint is activated
LF_t1 = 1e15; // switch trick to active PV constraint
LF_t2 = -1; // switch trick to active PV constraint
LF_Ppark1_W = 0; // Active power in W of PQ constraint. The value is set to zero
LF_Ppark2_W = init_Pturbine_pu * Spark_VA * 0.967; // Active power in W of PV constraint. The value is multiplied by a factor of losses.
}


}

Loading

0 comments on commit 0b86730

Please sign in to comment.