Commit 133ae3a5 authored by Viviane Sapucaia's avatar Viviane Sapucaia
Browse files

Created voltage behind reactance model - in construction

parent 3a3bd0f9
......@@ -20,5 +20,6 @@
#include "Components/IdealVoltageSource.h"
#include "Components/RxLine.h"
#include "Components/PiLine.h"
#include "Components/VoltageBehindReactanceEMT.h"
#endif // !COMPONENTS_H
\ No newline at end of file
#include "VoltageBehindReactanceEMT.h"
using namespace DPsim;
VoltageBehindReactanceEMT::VoltageBehindReactanceEMT(std::string name, int node1, int node2, int node3,
Real nomPower, Real nomVolt, Real nomFreq, int poleNumber, Real nomFieldCur,
Real Rs, Real Ll, Real Lmd, Real Lmq,
Real Rfd, Real Llfd, Real Rkd, Real Llkd,
Real Rkq1, Real Llkq1, Real Rkq2, Real Llkq2,
Real inertia) {
this->mNode1 = node1 - 1;
this->mNode2 = node2 - 1;
this->mNode3 = node3 - 1;
mNomPower = nomPower;
mNomVolt = nomVolt;
mNomFreq = nomFreq;
mPoleNumber = poleNumber;
mNomFieldCur = nomFieldCur;
mRs = Rs;
mLl = Ll;
mLmd = Lmd;
mLmq = Lmq;
mRfd = Rfd;
mLlfd = Llfd;
mRkd = Rkd;
mLlkd = Llkd;
mRkq1 = Rkq1;
mLlkq1 = Llkq1;
mRkq2 = Rkq2;
mLlkq2 = Llkq2;
mH = inertia;
//Dynamic mutual inductances
mDLmd = 1 / (1 / mLmd + 1 / mLlfd + 1 / mLlkd);
mDLmq = 1 / (1 / mLmq + 1 / mLlkq1 + 1 / mLlkq2);
mLa = (mDLmq + mDLmd) / 3;
mLb = (mDLmd - mDLmq) / 3;
}
void VoltageBehindReactanceEMT::init(Real om, Real dt,
Real initActivePower, Real initReactivePower, Real initTerminalVolt, Real initVoltAngle) {
mResistanceMat <<
mRs, 0, 0,
0, mRs, 0,
0, 0, mRs;
// steady state per unit initial value
initStates(initActivePower, initReactivePower, initTerminalVolt, initVoltAngle);
mVa = inverseParkTransform(mThetaMech, mVq, mVd, mV0)(0);
mVb = inverseParkTransform(mThetaMech, mVq, mVd, mV0)(1);
mVc = inverseParkTransform(mThetaMech, mVq, mVd, mV0)(2);
mIa = inverseParkTransform(mThetaMech, mIq, mId, mI0)(0);
mIb = inverseParkTransform(mThetaMech, mIq, mId, mI0)(1);
mIc = inverseParkTransform(mThetaMech, mIq, mId, mI0)(2);
}
void VoltageBehindReactanceEMT::initStates(Real initActivePower, Real initReactivePower,
Real initTerminalVolt, Real initVoltAngle) {
double init_P = initActivePower;
double init_Q = initReactivePower;
double init_S = sqrt(pow(init_P, 2.) + pow(init_Q, 2.));
double init_vt = initTerminalVolt;
double init_it = init_S / init_vt;
// power factor
double init_pf = acos(init_P / init_S);
// load angle
double init_delta = atan(((mLmq + mLl) * init_it * cos(init_pf) - mRs * init_it * sin(init_pf)) /
(init_vt + mRs * init_it * cos(init_pf) + (mLmq + mLl) * init_it * sin(init_pf)));
double init_delta_deg = init_delta / DPS_PI * 180;
// dq stator voltages and currents
double init_vd = init_vt * sin(init_delta);
double init_vq = init_vt * cos(init_delta);
double init_id = init_it * sin(init_delta + init_pf);
double init_iq = init_it * cos(init_delta + init_pf);
// rotor voltage and current
double init_ifd = (init_vq + mRs * init_iq + (mLmd + mLl) * init_id) / mLmd;
double init_vfd = mRfd * init_ifd;
// flux linkages
double init_psid = init_vq + mRs * init_iq;
double init_psiq = -init_vd - mRs * init_id;
double init_psifd = (mLmd + mLlfd) * init_ifd - mLmd * init_id;
double init_psid1 = mLmd * (init_ifd - init_id);
double init_psiq1 = -mLmq * init_iq;
double init_psiq2 = -mLmq * init_iq;
// rotor mechanical variables
double init_Te = init_P + mRs * pow(init_it, 2.);
mOmMech = 1;
mVd = init_vd;
mVq = init_vq;
mV0 = 0;
mVfd = init_vfd;
mVkd = 0;
mVkq1 = 0;
mVkq2 = 0;
mIq = init_iq;
mId = init_id;
mI0 = 0;
mIfd = init_ifd;
mIkd = 0;
mIkq1 = 0;
mIkq2 = 0;
mPsiq = init_psiq;
mPsid = init_psid;
mPsi0 = 0;
mPsifd = init_psifd;
mPsikd = init_psid1;
mPsikq1 = init_psiq1;
mPsikq2 = init_psiq2;
mDPsiq = mDLmq*(mPsikq1 / mLlkq1) + mDLmq*(mPsikq2 / mLlkq2);
mDPsid = mDLmd*(mPsifd / mLlfd) + mDLmd*(mPsikd / mLlkd);
mDVq = mOmMech*mDPsid + mDLmq*mRkq1*(mDPsiq - mPsikq1) / (mLlkq1*mLlkq1) +
mDLmq*mRkq2*(mDPsiq - mPsikq2) / (mLlkq2*mLlkq2) + (mRkq1 / (mLlkq1*mLlkq1) + mRkq2 / (mLlkq2*mLlkq2))*mDLmq*mDLmq*mIq;
mDVd = -mOmMech*mDPsiq + mDLmq*mRkd*(mDPsid - mPsikd) / (mLlkd*mLlkd) + (mDLmd / mLlfd)*mVfd +
mDLmd*mRfd*(mDPsiq - mPsifd) / (mLlfd*mLlfd) + (mRfd / (mLlfd*mLlfd) + mRkd / (mLlkd*mLlkd))*mDLmd*mDLmd*mId;
mDVa = inverseParkTransform(mThetaMech, mDVq, mDVd, 0)(0);
mDVb = inverseParkTransform(mThetaMech, mDVq, mDVd, 0)(1);
mDVc = inverseParkTransform(mThetaMech, mDVq, mDVd, 0)(2);
// Initialize mechanical angle
//mThetaMech = initVoltAngle + init_delta;
mThetaMech = initVoltAngle + init_delta - M_PI / 2;
//Initial inductance matrix
mDInductanceMat <<
mLl + mLa - mLb*cos(2 * mThetaMech), -mLa / 2 - mLb*cos(2 * mThetaMech - 2 * PI / 3), -mLa / 2 - mLb*cos(2 * mThetaMech + 2 * PI / 3),
-mLa / 2 - mLb*cos(2 * mThetaMech - 2 * PI / 3), mLl + mLa - mLb*cos(2 * mThetaMech - 4 * PI / 3), -mLa / 2 - mLb*cos(2 * mThetaMech),
-mLa / 2 - mLb*cos(2 * mThetaMech + 2 * PI / 3), -mLa / 2 - mLb*cos(2 * mThetaMech), mLl + mLa - mLb*cos(2 * mThetaMech + 4 * PI / 3);
}
void VoltageBehindReactanceEMT::step(SystemModel& system, Real fieldVoltage, Real mechPower, Real time) {
stepInPerUnit(system.getOmega(), system.getTimeStep(), fieldVoltage, mechPower, time, system.getNumMethod());
mVa = mVabc(0);
mVb = mVabc(1);
mVc = mVabc(2);
mIa = mVa / 1037.8378;
mIb = mVb / 1037.8378;
mIc = mVc / 1037.8378;
}
void VoltageBehindReactanceEMT::stepInPerUnit(Real om, Real dt, Real fieldVoltage, Real mechPower, Real time, NumericalMethod numMethod) {
mVabc <<
mVa,
mVb,
mVc;
mIabc <<
mIa,
mIb,
mIc;
mDVabc <<
mDVa,
mDVb,
mDVc;
mV_hist = (mResistanceMat - (2 / dt)*mDInductanceMat)*mIabc + mDVabc - mVabc;
mId = parkTransform(mThetaMech, mIa, mIb, mIc)(0);
mIq = parkTransform(mThetaMech, mIa, mIb, mIc)(1);
mI0 = parkTransform(mThetaMech, mIa, mIb, mIc)(2);
//usindg euler for flux
Real mPsimq = mDLmq*(mPsikq1 / mLlkq1 + mPsikq2 / mLlkq2 + mIq);
Real mPsimd = mDLmq*(mPsifd / mLlfd + mPsikd / mLlkd + mId);
mPsikq1 = mPsikq1 - dt*(mRkq1 / mLlkq1)*(mPsikq1 - mPsimq);
mPsikq2 = mPsikq2 - dt*(mRkq2 / mLlkq2)*(mPsikq2 - mPsimq);
mPsifd = mPsifd - dt*(mRfd / mLlfd)*(mPsifd - mPsimd) + mVfd;
mPsikd = mPsikd - dt*(mRkd / mLlkd)*(mPsikd - mPsimd);
// Calculate dynamic flux likages
mDPsiq = mDLmq*(mPsikq1 / mLlkq1) + mDLmq*(mPsikq2 / mLlkq2);
mDPsid = mDLmd*(mPsifd / mLlfd) + mDLmd*(mPsikd / mLlkd);
mDVq = mOmMech*mDPsid + mDLmq*mRkq1*(mDPsiq - mPsikq1) / (mLlkq1*mLlkq1) +
mDLmq*mRkq2*(mDPsiq - mPsikq2) / (mLlkq2*mLlkq2) + (mRkq1 / (mLlkq1*mLlkq1) + mRkq2 / (mLlkq2*mLlkq2))*mDLmq*mDLmq*mIq;
mDVd = -mOmMech*mDPsiq + mDLmq*mRkd*(mDPsid - mPsikd) / (mLlkd*mLlkd) + (mDLmd/mLlfd)*mVfd +
mDLmd*mRfd*(mDPsiq - mPsifd) / (mLlfd*mLlfd) + (mRfd / (mLlfd*mLlfd) + mRkd / (mLlkd*mLlkd))*mDLmd*mDLmd*mId;
mDVa = inverseParkTransform(mThetaMech, mDVq, mDVd, 0)(0);
mDVb = inverseParkTransform(mThetaMech, mDVq, mDVd, 0)(1);
mDVc = inverseParkTransform(mThetaMech, mDVq, mDVd, 0)(2);
mThetaMech_hist2 = mThetaMech_hist1;
mThetaMech_hist1 = mThetaMech;
mThetaMech = 2 * mThetaMech - mThetaMech_hist2;
mOmMech_hist2 = mOmMech_hist1;
mOmMech_hist1 = mOmMech;
mOmMech = 2 * mOmMech - mOmMech_hist2;
mVabc = (mResistanceMat + (2 / dt)*mDInductanceMat)*mIabc + mDVabc + mV_hist;
}
void VoltageBehindReactanceEMT::postStep(SystemModel& system) {
}
DPSMatrix VoltageBehindReactanceEMT::parkTransform(Real theta, double a, double b, double c) {
DPSMatrix dq0vector(3, 1);
// Park transform according to Kundur
double q, d, zero;
q = 2. / 3. * cos(theta) * a + 2. / 3. * cos(theta - 2. * M_PI / 3.)*b + 2. / 3. * cos(theta + 2. * M_PI / 3.)*c;
d = 2. / 3. * sin(theta)*a + 2. / 3. * sin(theta - 2. * M_PI / 3.)*b + 2. / 3. * sin(theta + 2. * M_PI / 3.)*c;
zero = 1. / 3. * a, 1. / 3. * b, 1. / 3. * c;
dq0vector << q,
d,
0;
return dq0vector;
}
DPSMatrix VoltageBehindReactanceEMT::inverseParkTransform(Real theta, double q, double d, double zero) {
DPSMatrix abcVector(3, 1);
double a, b, c;
// Park transform according to Kundur
a = cos(theta)*d + sin(theta)*q + 1.*zero;
b = cos(theta - 2. * M_PI / 3.)*d + sin(theta - 2. * M_PI / 3.)*q + 1.*zero;
c = cos(theta + 2. * M_PI / 3.)*d + sin(theta + 2. * M_PI / 3.)*q + 1.*zero;
abcVector << a,
b,
c;
return abcVector;
}
\ No newline at end of file
#ifndef VOLTAGEBEHINDREACTANCEEMT_H
#define VOLTAGEBEHINDREACTANCE_H
#include "BaseComponent.h"
#include "ComponentCommons.h"
namespace DPsim {
/// Synchronous generator model
/// If parInPerUnit is not set, the parameters have to be given with their respective stator or rotor
/// referred values. The calculation to per unit is performed in the initialization.
/// The case where parInPerUnit is not set will be implemented later.
/// parameter names include underscores and typical variables names found in literature instead of
/// descriptive names in order to shorten formulas and increase the readability
class VoltageBehindReactanceEMT : public BaseComponent {
protected:
// ### Machine parameters ###
/// nominal power Pn [VA]
double mNomPower;
/// nominal voltage Vn [V] (RMS)
double mNomVolt;
/// nominal frequency fn [Hz]
double mNomFreq;
/// nominal field current Ifn [A]
double mNomFieldCur;
/// stator resistance Rs[Ohm]
double mRs;
/// leakage inductance Ll [H]
double mLl;
/// d-axis mutual inductance Lmd [H]
double mLmd;
/// unsaturated d-axis mutual inductance Lmd [H]
double mLmd0;
/// q-axis mutual inductance Lmq [H]
double mLmq;
/// unsaturated q-axis mutual inductance Lmq [H]
double mLmq0;
/// field resistance Rfd [Ohm]
double mRfd;
/// field leakage inductance Llfd [H]
double mLlfd;
/// d-axis damper resistance Rkd [Ohm]
double mRkd;
/// d-axis damper leakage inductance Llkd [H]
double mLlkd;
/// q-axis damper resistance 1 Rkq1 [Ohm]
double mRkq1;
/// q-axis damper leakage inductance 1 Llkq1 [H]
double mLlkq1;
/// q-axis damper resistance 2 Rkq2 [Ohm]
double mRkq2;
/// q-axis damper leakage inductance 2 Llkq2 [H]
double mLlkq2;
/// q winding inductance
double mLaq;
/// d winding inductance
double mLad;
/// d dynamic inductance
double mDLmd;
/// q dynamic inductance
double mDLmq;
/// d dynamic flux
double mDPsid;
/// q dynamic flux
double mDPsiq;
/// Dynamic dq voltages
double mDVq;
double mDVd;
/// Dynamic abc voltages
double mDVa;
double mDVb;
double mDVc;
//Theta mech hist
double mThetaMech_hist1;
double mThetaMech_hist2;
//Theta omega hist
double mOmMech_hist1;
double mOmMech_hist2;
/// inertia J [kg*m^2]
double mJ;
/// number of poles
int mPoleNumber;
/// inertia coefficient H
double mH;
// ### Stator base values ###
/// specifies if the machine parameters are transformed to per unit
SynchGenStateType mStateType;
/// base stator voltage
double mBase_v;
/// base stator voltage RMS
double mBase_V_RMS;
/// base stator current
double mBase_i;
/// base stator current RMS
double mBase_I_RMS;
/// base stator impedance
double mBase_Z;
/// base electrical angular frequency
double mBase_OmElec;
/// base mechanical angular frequency
double mBase_OmMech;
/// base stator inductance
double mBase_L;
/// base torque
double mBase_T;
/// base flux linkage
double mBase_Psi;
/// base field current
double mBase_ifd;
/// base field voltage
double mBase_vfd;
/// base field impedance
double mBase_Zfd;
/// base field inductance
double mBase_Lfd;
// ### State variables ###
/// rotor speed omega_r
double mOmMech;
/// theta
double mThetaMech;
double mThetaMech2;
/// mechanical Power Pm [W]
double mMechPower;
/// mechanical torque
double mMechTorque;
/// electrical torque
double mElecTorque;
/// mechanical torque
double mMechTorque_past;
/// electrical torque
double mElecTorque_past;
/// rotor speed omega_r
double mOmMech_past;
/// voltage vector q d 0 kq1 kq2 df kd
DPSMatrix mVoltages2 = DPSMatrix::Zero(7, 1);
/// flux linkage vector
DPSMatrix mFluxes2 = DPSMatrix::Zero(7, 1);
/// current vector
DPSMatrix mCurrents2 = DPSMatrix::Zero(7, 1);
/// voltage vector q d 0 fd kd kq1 kq2
double mVd;
double mVq;
double mV0;
double mVfd;
double mVkd;
double mVkq1;
double mVkq2;
/// current vector q d 0 fd kd kq1 kq2
double mId;
double mIq;
double mI0;
double mIfd;
double mIkd;
double mIkq1;
double mIkq2;
double mId_past;
double mIq_past;
/// flux linkage vector q d 0 fd kd kq1 kq2
double mPsid;
double mPsiq;
double mPsi0;
double mPsifd;
double mPsikd;
double mPsikq1;
double mPsikq2;
double mPsid_past;
double mPsiq_past;
/// Interface voltage vector
double mVa;
double mVb;
double mVc;
/// Interface voltage vector
double mIa;
double mIb;
double mIc;
// ### Useful Matrices ###
/// inductance matrix
DPSMatrix mDInductanceMat = DPSMatrix::Zero(3, 3);
/// resistance matrix
DPSMatrix mResistanceMat = DPSMatrix::Zero(3, 3);
/// reactance matrix
//DPSMatrix mReactanceMat = DPSMatrix::Zero(7, 7);
/// omega - flux matrix
//DPSMatrix mOmegaFluxMat = DPSMatrix::Zero(7, 7);
/// matrix for reversing stator current directions in calculations with respect to other currents
//DPSMatrix mReverseCurrents = DPSMatrix::Zero(7, 7);
double mLa;
double mLb;
//Historical term of voltage
DPSMatrix mV_hist = DPSMatrix::Zero(3, 1);
//Historical term of voltage
DPSMatrix mVabc = DPSMatrix::Zero(3, 1);
//Historical term of current
DPSMatrix mIabc = DPSMatrix::Zero(3, 1);
//Historical term of voltage
DPSMatrix mDVabc = DPSMatrix::Zero(3, 1);
public:
VoltageBehindReactanceEMT() { };
/// Initializes the per unit or stator referred machine parameters with the machine parameters given in per unit or
/// stator referred parameters depending on the setting of parameter type.
/// The initialization mode depends on the setting of state type.
VoltageBehindReactanceEMT(std::string name, int node1, int node2, int node3,
Real nomPower, Real nomVolt, Real nomFreq, int poleNumber, Real nomFieldCur,
Real Rs, Real Ll, Real Lmd, Real Lmq,
Real Rfd, Real Llfd, Real Rkd, Real Llkd,
Real Rkq1, Real Llkq1, Real Rkq2, Real Llkq2,
Real inertia);
/// Initializes states in per unit or stator referred variables depending on the setting of the state type.
/// Function parameters have to be given in real units.
void init(Real om, Real dt,
Real initActivePower, Real initReactivePower, Real initTerminalVolt, Real initVoltAngle);
/// Initializes states in per unit. All machine parameters are assumed to be in per unit.
/// Function parameters have to be given in real units.
void initStates(Real initActivePower, Real initReactivePower,
Real initTerminalVolt, Real initVoltAngle);
/// Performs an Euler forward step with the state space model of a synchronous generator
/// to calculate the flux and current from the voltage vector.
void step(SystemModel& system, Real fieldVoltage, Real mechPower, Real time);
/// Performs an Euler forward step with the state space model of a synchronous generator
/// to calculate the flux and current from the voltage vector in per unit.
void stepInPerUnit(Real om, Real dt, Real fieldVoltage, Real mechPower, Real time, NumericalMethod numMethod);
/// Retrieves calculated voltage from simulation for next step
void postStep(SystemModel& system);
/// Park transform as described in Krause
DPSMatrix parkTransform(Real theta, double a, double b, double c);
/// Inverse Park transform as described in Krause
DPSMatrix inverseParkTransform(Real theta, double q, double d, double zero);
DPSMatrix& getVoltages() { return mVabc; }
DPSMatrix& getCurrents() { return mIabc; }
//DPSMatrix& getFluxes() { return mFluxes2; }
void init(Real om, Real dt) { }
void applySystemMatrixStamp(SystemModel& system) { }
void applyRightSideVectorStamp(SystemModel& system) { }
void step(SystemModel& system, Real time) { }
};
}
#endif
......@@ -2,6 +2,7 @@
#include <string>
#include "NetlistSim.h"
#include "Examples/SynchronGenUnitTest.h"
#include "Examples/SyncGenUnitTestVBR.h"
#include "Examples/ReferenceCircuits.h"
#include "Examples/ShmemTest.h"
#include <complex>
......@@ -44,7 +45,8 @@ int main(int argc, char* argv[]) {
//SynGenUnitTestPhaseToPhaseFault();
//SynGenUnitTestThreePhaseFault();
SynGenDPUnitTestBalancedResLoad();
//SynGenDPUnitTestBalancedResLoad();
SynGenUnitTestVBR();
//RXLineResLoad();
//VarFreqRXLineResLoad();
......
#include "SyncGenUnitTestVBR.h"
#include "../Simulation.h"
#include "../Utilities.h"
using namespace DPsim;
void DPsim::SynGenUnitTestVBR() {
// Define Object for saving data on a file
Logger log("Logs/logVBR.txt"),
synGenLogVolt("Logs/data_synGen_volt.csv"),
synGenLogCurr("Logs/data_synGen_curr.csv");