Commit 40170493 authored by Viviane Sapucaia's avatar Viviane Sapucaia
Browse files

Added simulation of synchronous generator solving everything in dq.

parent 40510e51
......@@ -9,6 +9,7 @@
#include "Components/VoltSourceRes.h"
#include "Components/InterfacedInductor.h"
#include "Components/SynchronGeneratorEMT.h"
#include "Components/SynchronGeneratorEMTdq.h"
#include "Components/SynchronGenerator.h"
#include "Components/ComponentCommons.h"
#include "Components/VoltSourceResFreq.h"
......
......@@ -206,7 +206,7 @@ void SynchronGeneratorEMT::step(SystemModel& system, Real fieldVoltage, Real mec
}
void SynchronGeneratorEMT::stepInPerUnit(Real om, Real dt, Real fieldVoltage, Real mechPower, Real time, NumericalMethod numMethod) {
mVa = (1 / mBase_v) * mVa;
mVb = (1 / mBase_v) * mVb;
mVc = (1 / mBase_v) * mVc;
......@@ -358,24 +358,26 @@ void SynchronGeneratorEMT::stepInPerUnit(Real om, Real dt, Real fieldVoltage, Re
else if (numMethod == NumericalMethod::Trapezoidal_flux){
//mMechPower = mechPower / mNomPower;
//mMechTorque = mMechPower / mOmMech;
Real mElecTorque_hist = (mPsid*mIq - mPsiq*mId);
//mElecTorque = (mPsid*mIq - mPsiq*mId);
// Save last step values
//Real mOmMech_hist = mOmMech;
//Real mElecTorque_hist = mElecTorque;
//Real mMechTorque_hist = mMechTorque;
// Euler step forward
//mOmMech = mOmMech + dt * (1 / (2 * mH) * (mMechTorque - mElecTorque));
DPSMatrix A = mBase_OmElec*(mResistanceMat*mReactanceMat - mOmMech*mOmegaFluxMat);
DPSMatrix I(7, 7);
I <<
1, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 1;
DPSMatrix I = DPSMatrix::Identity(7, 7);
DPSMatrix Aux = I + (dt / 2) * A;
DPSMatrix Aux2 = I - (dt / 2) * A;
DPSMatrix InvAux = Aux2.inverse();
DPSMatrix Fluxes(7, 1);
Fluxes(0, 0) = mPsiq;
Fluxes(1, 0) = mPsid;
......@@ -403,7 +405,7 @@ void SynchronGeneratorEMT::stepInPerUnit(Real om, Real dt, Real fieldVoltage, Re
mPsikq2 = Fluxes(4, 0);
mPsifd = Fluxes(5, 0);
mPsikd = Fluxes(6, 0);
//Calculation of currents based on inverse of inductance matrix
mId = ((mLlfd*mLlkd + mLmd*(mLlfd + mLlkd))*mPsid - mLmd*mLlkd*mPsifd - mLlfd*mLmd*mPsikd) / detLd;
mIfd = (mLlkd*mLmd*mPsid - (mLl*mLlkd + mLmd*(mLl + mLlkd))*mPsifd + mLmd*mLl*mPsikd) / detLd;
......@@ -412,23 +414,30 @@ void SynchronGeneratorEMT::stepInPerUnit(Real om, Real dt, Real fieldVoltage, Re
mIkq1 = (mLmq*mLlkq2*mPsiq - (mLmq*(mLlkq2 + mLl) + mLl*mLlkq2)*mPsikq1 + mLmq*mLl*mPsikq2) / detLq;
mIkq2 = (mLmq*mLlkq1*mPsiq + mLmq*mLl*mPsikq1 - (mLmq*(mLlkq1 + mLl) + mLl*mLlkq1)*mPsikq2) / detLq;
mI0 = -mPsi0 / mLl;
//mMechTorque = mMechPower / mOmMech;
//mElecTorque = (mPsid*mIq - mPsiq*mId);
//Real Hist_term = (dt / (4 * mH))*(mMechTorque_hist - mElecTorque_hist) + mOmMech_hist;
//Real a = 1;
//Real b = -(Hist_term - (dt / (4 * mH))*mElecTorque);
//Real c = -(dt / (4 * mH))*mMechPower;
mMechPower = mechPower / mNomPower;
mMechTorque = mMechPower / mOmMech;
mElecTorque = (mPsid*mIq - mPsiq*mId);
//mOmMech = (-b + sqrt(b*b - 4 * a*c)) / (2.*a);
// Trapezoidal rule
Real mOmMech_hist = mOmMech;
mOmMech = mOmMech + dt / (4 * mH) * (mMechTorque - mElecTorque_hist + mMechTorque - mElecTorque);
//mOmMech = mOmMech + dt / (4 * mH) * (mMechTorque - mElecTorque_hist + mMechTorque - mElecTorque);
// Update mechanical rotor angle with respect to electrical angle
mThetaMech = mThetaMech + dt / 2 * (mOmMech * mBase_OmMech + mOmMech_hist * mBase_OmMech);
//mThetaMech = mThetaMech + (dt / 2.) * (mOmMech * mBase_OmMech + mOmMech_hist * mBase_OmMech)
mThetaMech = mThetaMech + dt * (mOmMech * mBase_OmMech);
}
else if (numMethod == NumericalMethod::Trapezoidal_current) {
Real mElecTorque_hist = (mPsid*mIq - mPsiq*mId);
//Real mElecTorque_hist = (mPsid*mIq - mPsiq*mId);
DPSMatrix A = mBase_OmElec*(mReactanceMat*mResistanceMat);
DPSMatrix B = mBase_OmElec*mReactanceMat;
......@@ -439,16 +448,6 @@ void SynchronGeneratorEMT::stepInPerUnit(Real om, Real dt, Real fieldVoltage, Re
DPSMatrix I = DPSMatrix::Identity(7, 7);
//DPSMatrix I(7, 7);
//I <<
// 1, 0, 0, 0, 0, 0, 0,
// 0, 1, 0, 0, 0, 0, 0,
// 0, 0, 1, 0, 0, 0, 0,
// 0, 0, 0, 1, 0, 0, 0,
// 0, 0, 0, 0, 1, 0, 0,
// 0, 0, 0, 0, 0, 1, 0,
// 0, 0, 0, 0, 0, 0, 1;
DPSMatrix Aux = I + (dt / 2) * A;
DPSMatrix Aux2 = I - (dt / 2) * A;
DPSMatrix InvAux = Aux2.inverse();
......@@ -491,36 +490,21 @@ void SynchronGeneratorEMT::stepInPerUnit(Real om, Real dt, Real fieldVoltage, Re
mPsikd = -mLmd*mId + mLmd*mIfd + (mLlkd + mLmd)*mIkd;
mMechPower = mechPower / mNomPower;
mMechTorque = mMechPower / mOmMech;
//mMechPower = mechPower / mNomPower;
//mMechTorque = mMechPower / mOmMech;
mElecTorque = (mPsid*mIq - mPsiq*mId);
//mElecTorque = (mPsid*mIq - mPsiq*mId);
// Trapezoidal rule
Real mOmMech_hist = mOmMech;
mOmMech = mOmMech + dt / (4 * mH) * (mMechTorque - mElecTorque_hist + mMechTorque - mElecTorque);
//Real mOmMech_hist = mOmMech;
//mOmMech = mOmMech + dt / (4 * mH) * (mMechTorque - mElecTorque_hist + mMechTorque - mElecTorque);
// Update mechanical rotor angle with respect to electrical angle
mThetaMech = mThetaMech + dt / 2 * (mOmMech * mBase_OmMech + mOmMech_hist * mBase_OmMech);
//mThetaMech = mThetaMech + dt / 2 * (mOmMech * mBase_OmMech + mOmMech_hist * mBase_OmMech);
mThetaMech = mThetaMech + dt * (mOmMech * mBase_OmMech);
}
//if (dtPsid < 0.000001)
// dtPsid = 0;
//if (dtPsiq < 0.000001)
// dtPsiq = 0;
//if (dtPsi0 < 0.000001)
// dtPsi0 = 0;
//if (dtPsifd < 0.000001)
// dtPsifd = 0;
//if (dtPsikd < 0.000001)
// dtPsikd = 0;
//if (dtPsikq1 < 0.000001)
// dtPsikq1 = 0;
//if (dtPsikq2 < 0.000001)
// dtPsikq2 = 0;
mIa = mBase_i * inverseParkTransform2(mThetaMech, mId, mIq, mI0)(0);
mIb = mBase_i * inverseParkTransform2(mThetaMech, mId, mIq, mI0)(1);
mIc = mBase_i * inverseParkTransform2(mThetaMech, mId, mIq, mI0)(2);
......@@ -548,6 +532,7 @@ void SynchronGeneratorEMT::stepInPerUnit(Real om, Real dt, Real fieldVoltage, Re
mPsikq2,
mPsifd,
mPsikd;
}
void SynchronGeneratorEMT::postStep(SystemModel& system) {
......@@ -569,6 +554,8 @@ void SynchronGeneratorEMT::postStep(SystemModel& system) {
else {
mVc = 0;
}
}
......
#include "SynchronGeneratorEMTdq.h"
using namespace DPsim;
SynchronGeneratorEMTdq::SynchronGeneratorEMTdq(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 Lmd0, Real Lmq, Real Lmq0,
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;
// base stator values
mBase_V_RMS = mNomVolt / sqrt(3);
mBase_v = mBase_V_RMS * sqrt(2);
mBase_I_RMS = mNomPower / (3 * mBase_V_RMS);
mBase_i = mBase_I_RMS * sqrt(2);
mBase_Z = mBase_v / mBase_i;
mBase_OmElec = 2 * DPS_PI * mNomFreq;
mBase_OmMech = mBase_OmElec / (mPoleNumber / 2);
mBase_L = mBase_Z / mBase_OmElec;
mBase_Psi = mBase_L * mBase_i;
mBase_T = mNomPower / mBase_OmMech;
// steady state per unit initial value
initWithPerUnitParam(Rs, Ll, Lmd, Lmd0, Lmq, Lmq0, Rfd, Llfd, Rkd, Llkd, Rkq1, Llkq1, Rkq2, Llkq2, inertia);
}
void SynchronGeneratorEMTdq::initWithPerUnitParam(
Real Rs, Real Ll, Real Lmd, Real Lmd0, Real Lmq, Real Lmq0,
Real Rfd, Real Llfd, Real Rkd, Real Llkd,
Real Rkq1, Real Llkq1, Real Rkq2, Real Llkq2,
Real H) {
// base rotor values
mBase_ifd = Lmd * mNomFieldCur;
mBase_vfd = mNomPower / mBase_ifd;
mBase_Zfd = mBase_vfd / mBase_ifd;
mBase_Lfd = mBase_Zfd / mBase_OmElec;
mRs = Rs;
mLl = Ll;
mLmd = Lmd;
mLmd0 = Lmd0;
mLmq = Lmq;
mLmq0 = Lmq0;
mRfd = Rfd;
mLlfd = Llfd;
mRkd = Rkd;
mLlkd = Llkd;
mRkq1 = Rkq1;
mLlkq1 = Llkq1;
mRkq2 = Rkq2;
mLlkq2 = Llkq2;
mH = H;
// Additional inductances according to Krause
mLaq = 1 / (1 / mLmq + 1 / mLl + 1 / mLlkq1 + 1 / mLlkq2);
mLad = 1 / (1 / mLmd + 1 / mLl + 1 / mLlkd + 1 / mLlfd);
// Determinant of Ld (inductance matrix of d axis)
detLd = (mLmd + mLl)*(-mLlfd*mLlkd - mLlfd*mLmd - mLmd*mLlkd) + mLmd*mLmd*(mLlfd + mLlkd);
// Determinant of Lq (inductance matrix of q axis)
detLq = -mLmq*mLlkq2*(mLlkq1 + mLl) - mLl*mLlkq1*(mLlkq2 + mLmq);
}
void SynchronGeneratorEMTdq::init(Real om, Real dt,
Real initActivePower, Real initReactivePower, Real initTerminalVolt, Real initVoltAngle) {
// Create matrices for state space representation
mInductanceMat <<
-(mLl + mLmq), 0, 0, mLmq, mLmq, 0, 0,
0, -(mLl + mLmd), 0, 0, 0, mLmd, mLmd,
0, 0, -mLl, 0, 0, 0, 0,
-mLmq, 0, 0, mLlkq1 + mLmq, mLmq, 0, 0,
-mLmq, 0, 0, mLmq, mLlkq2 + mLmq, 0, 0,
0, -mLmd, 0, 0, 0, mLlfd + mLmd, mLmd,
0, -mLmd, 0, 0, 0, mLmd, mLlkd + mLmd;
mResistanceMat <<
mRs, 0, 0, 0, 0, 0, 0,
0, mRs, 0, 0, 0, 0, 0,
0, 0, mRs, 0, 0, 0, 0,
0, 0, 0, -mRkq1, 0, 0, 0,
0, 0, 0, 0, -mRkq2, 0, 0,
0, 0, 0, 0, 0, -mRfd, 0,
0, 0, 0, 0, 0, 0, -mRkd;
mOmegaFluxMat <<
0, 1, 0, 0, 0, 0, 0,
-1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0;
mReactanceMat = mInductanceMat.inverse();
// steady state per unit initial value
initStatesInPerUnit(initActivePower, initReactivePower, initTerminalVolt, initVoltAngle);
/*mVa = inverseParkTransform2(mThetaMech, mVd* mBase_v, mVq* mBase_v, mV0* mBase_v)(0);
mVb = inverseParkTransform2(mThetaMech, mVd* mBase_v, mVq* mBase_v, mV0* mBase_v)(1);
mVc = inverseParkTransform2(mThetaMech, mVd* mBase_v, mVq* mBase_v, mV0* mBase_v)(2);
mIa = inverseParkTransform2(mThetaMech, mId* mBase_i, mIq* mBase_i, mI0* mBase_i)(0);
mIb = inverseParkTransform2(mThetaMech, mId* mBase_i, mIq* mBase_i, mI0* mBase_i)(1);
mIc = inverseParkTransform2(mThetaMech, mId* mBase_i, mIq* mBase_i, mI0* mBase_i)(2);*/
}
void SynchronGeneratorEMTdq::initStatesInPerUnit(Real initActivePower, Real initReactivePower,
Real initTerminalVolt, Real initVoltAngle) {
double init_P = initActivePower / mNomPower;
double init_Q = initReactivePower / mNomPower;
double init_S = sqrt(pow(init_P, 2.) + pow(init_Q, 2.));
double init_vt = initTerminalVolt / mBase_v;
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;
// Initialize mechanical angle
//mThetaMech = initVoltAngle + init_delta;
mThetaMech = initVoltAngle + init_delta - M_PI / 2;
}
void SynchronGeneratorEMTdq::step(SystemModel& system, Real fieldVoltage, Real mechPower, Real time) {
stepInPerUnit(system.getOmega(), system.getTimeStep(), fieldVoltage, mechPower, time, system.getNumMethod());
//// Update current source accordingly
//if (mNode1 >= 0) {
// system.addRealToRightSideVector(mNode1, mIa);
//}
//if (mNode2 >= 0) {
// system.addRealToRightSideVector(mNode2, mIb);
//}
//if (mNode3 >= 0) {
// system.addRealToRightSideVector(mNode3, mIc);
//}
mVq = (1037.8378 / mBase_Z)*mIq;
mVd = (1037.8378 / mBase_Z)*mId;
}
void SynchronGeneratorEMTdq::stepInPerUnit(Real om, Real dt, Real fieldVoltage, Real mechPower, Real time, NumericalMethod numMethod) {
/*mVa = (1 / mBase_v) * mVa;
mVb = (1 / mBase_v) * mVb;
mVc = (1 / mBase_v) * mVc;
mIa = (1 / mBase_i) * mIa;
mIb = (1 / mBase_i) * mIb;
mIc = (1 / mBase_i) * mIc;*/
mVfd = fieldVoltage / mBase_v;
// TODO calculate effect of changed field voltage
// dq-transform of interface voltage
//mVd = parkTransform2(mThetaMech, mVa, mVb, mVc)(0);
//mVq = parkTransform2(mThetaMech, mVa, mVb, mVc)(1);
//mV0 = parkTransform2(mThetaMech, mVa, mVb, mVc)(2);
if (numMethod == NumericalMethod::Euler) {
mMechPower = mechPower / mNomPower;
mMechTorque = mMechPower / mOmMech;
mElecTorque = (mPsid*mIq - mPsiq*mId);
// Euler step forward
mOmMech = mOmMech + dt * (1 / (2 * mH) * (mMechTorque - mElecTorque));
double dtPsid = mVd + mRs*mId + mPsiq*mOmMech;
double dtPsiq = mVq + mRs*mIq - mPsid*mOmMech;
double dtPsi0 = mV0 + mRs*mI0;
double dtPsifd = mVfd - mRfd*mIfd;
double dtPsikd = -mRkd*mIkd;
double dtPsikq1 = -mRkq1*mIkq1;
double dtPsikq2 = -mRkq2*mIkq2;
mPsid = mPsid + dt*mBase_OmElec*dtPsid;
mPsiq = mPsiq + dt*mBase_OmElec*dtPsiq;
mPsi0 = mPsi0 + dt*mBase_OmElec*dtPsi0;
mPsifd = mPsifd + dt*mBase_OmElec*dtPsifd;
mPsikd = mPsikd + dt*mBase_OmElec*dtPsikd;
mPsikq1 = mPsikq1 + dt*mBase_OmElec*dtPsikq1;
mPsikq2 = mPsikq2 + dt*mBase_OmElec*dtPsikq2;
//Calculation of currents based on inverse of inductance matrix
mId = ((mLlfd*mLlkd + mLmd*(mLlfd + mLlkd))*mPsid - mLmd*mLlkd*mPsifd - mLlfd*mLmd*mPsikd) / detLd;
mIfd = (mLlkd*mLmd*mPsid - (mLl*mLlkd + mLmd*(mLl + mLlkd))*mPsifd + mLmd*mLl*mPsikd) / detLd;
mIkd = (mLmd*mLlfd*mPsid + mLmd*mLl*mPsifd - (mLmd*(mLlfd + mLl) + mLl*mLlfd)*mPsikd) / detLd;
mIq = ((mLlkq1*mLlkq2 + mLmq*(mLlkq1 + mLlkq2))*mPsiq - mLmq*mLlkq2*mPsikq1 - mLmq*mLlkq1*mPsikq2) / detLq;
mIkq1 = (mLmq*mLlkq2*mPsiq - (mLmq*(mLlkq2 + mLl) + mLl*mLlkq2)*mPsikq1 + mLmq*mLl*mPsikq2) / detLq;
mIkq2 = (mLmq*mLlkq1*mPsiq + mLmq*mLl*mPsikq1 - (mLmq*(mLlkq1 + mLl) + mLl*mLlkq1)*mPsikq2) / detLq;
mI0 = -mPsi0 / mLl;
// Update mechanical rotor angle with respect to electrical angle
mThetaMech = mThetaMech + dt * (mOmMech * mBase_OmMech);
}
else if (numMethod == NumericalMethod::AdamBashforth) {
//Two steps Adams-Bashforth
if (time < dt) {
// calculate mechanical states
mMechPower = mechPower / mNomPower;
mMechTorque = mMechPower / mOmMech;
mElecTorque = (mPsid*mIq - mPsiq*mId);
mOmMech_past = mOmMech;
mOmMech = mOmMech + dt * (1 / (2 * mH) * (mMechTorque - mElecTorque));
double dtPsid = mVd + mRs*mId + mPsiq*mOmMech;
double dtPsiq = mVq + mRs*mIq - mPsid*mOmMech;
double dtPsi0 = mV0 + mRs*mI0;
double dtPsifd = mVfd - mRfd*mIfd;
double dtPsikd = -mRkd*mIkd;
double dtPsikq1 = -mRkq1*mIkq1;
double dtPsikq2 = -mRkq2*mIkq2;
mPsid_past = mPsid;
mPsiq_past = mPsiq;
mPsid = mPsid + dt*mBase_OmElec*dtPsid;
mPsiq = mPsiq + dt*mBase_OmElec*dtPsiq;
mPsi0 = mPsi0 + dt*mBase_OmElec*dtPsi0;
mPsifd = mPsifd + dt*mBase_OmElec*dtPsifd;
mPsikd = mPsikd + dt*mBase_OmElec*dtPsikd;
mPsikq1 = mPsikq1 + dt*mBase_OmElec*dtPsikq1;
mPsikq2 = mPsikq2 + dt*mBase_OmElec*dtPsikq2;
//Calculation of currents based on inverse of inductance matrix
mId_past = mId;
mIq_past = mIq;
mId = ((mLlfd*mLlkd + mLmd*(mLlfd + mLlkd))*mPsid - mLmd*mLlkd*mPsifd - mLlfd*mLmd*mPsikd) / detLd;
mIfd = (mLlkd*mLmd*mPsid - (mLl*mLlkd + mLmd*(mLl + mLlkd))*mPsifd + mLmd*mLl*mPsikd) / detLd;
mIkd = (mLmd*mLlfd*mPsid + mLmd*mLl*mPsifd - (mLmd*(mLlfd + mLl) + mLl*mLlfd)*mPsikd) / detLd;
mIq = ((mLlkq1*mLlkq2 + mLmq*(mLlkq1 + mLlkq2))*mPsiq - mLmq*mLlkq2*mPsikq1 - mLmq*mLlkq1*mPsikq2) / detLq;
mIkq1 = (mLmq*mLlkq2*mPsiq - (mLmq*(mLlkq2 + mLl) + mLl*mLlkq2)*mPsikq1 + mLmq*mLl*mPsikq2) / detLq;
mIkq2 = (mLmq*mLlkq1*mPsiq + mLmq*mLl*mPsikq1 - (mLmq*(mLlkq1 + mLl) + mLl*mLlkq1)*mPsikq2) / detLq;
mI0 = -mPsi0 / mLl;
// Update mechanical rotor angle with respect to electrical angle
mThetaMech = mThetaMech + dt * (mOmMech * mBase_OmMech);
}
else {
// calculate mechanical states
mMechPower = mechPower / mNomPower;
mMechTorque = mMechPower / mOmMech;
mMechTorque_past = mMechPower / mOmMech_past;
mElecTorque = (mPsid*mIq - mPsiq*mId);
mElecTorque_past = (mPsid_past*mIq - mPsiq_past*mId);
mOmMech_past = mOmMech;
mOmMech = mOmMech + (3. / 2.)*dt* (1 / (2 * mH) * (mMechTorque - mElecTorque)) - (1. / 2.)*dt* (1 / (2 * mH) * (mMechTorque_past - mElecTorque_past));
double dtPsid = mVd + mRs*mId + mPsiq*mOmMech;
double dtPsiq = mVq + mRs*mIq - mPsid*mOmMech;
double dtPsi0 = mV0 + mRs*mI0;
double dtPsifd = mVfd - mRfd*mIfd;
double dtPsikd = -mRkd*mIkd;
double dtPsikq1 = -mRkq1*mIkq1;
double dtPsikq2 = -mRkq2*mIkq2;
mPsid_past = mPsid;
mPsiq_past = mPsiq;
mPsid = mPsid + dt*mBase_OmElec*dtPsid;
mPsiq = mPsiq + dt*mBase_OmElec*dtPsiq;
mPsi0 = mPsi0 + dt*mBase_OmElec*dtPsi0;
mPsifd = mPsifd + dt*mBase_OmElec*dtPsifd;
mPsikd = mPsikd + dt*mBase_OmElec*dtPsikd;
mPsikq1 = mPsikq1 + dt*mBase_OmElec*dtPsikq1;
mPsikq2 = mPsikq2 + dt*mBase_OmElec*dtPsikq2;
//Calculation of currents based on inverse of inductance matrix
mId_past = mId;
mIq_past = mIq;
mId = ((mLlfd*mLlkd + mLmd*(mLlfd + mLlkd))*mPsid - mLmd*mLlkd*mPsifd - mLlfd*mLmd*mPsikd) / detLd;
mIfd = (mLlkd*mLmd*mPsid - (mLl*mLlkd + mLmd*(mLl + mLlkd))*mPsifd + mLmd*mLl*mPsikd) / detLd;
mIkd = (mLmd*mLlfd*mPsid + mLmd*mLl*mPsifd - (mLmd*(mLlfd + mLl) + mLl*mLlfd)*mPsikd) / detLd;
mIq = ((mLlkq1*mLlkq2 + mLmq*(mLlkq1 + mLlkq2))*mPsiq - mLmq*mLlkq2*mPsikq1 - mLmq*mLlkq1*mPsikq2) / detLq;
mIkq1 = (mLmq*mLlkq2*mPsiq - (mLmq*(mLlkq2 + mLl) + mLl*mLlkq2)*mPsikq1 + mLmq*mLl*mPsikq2) / detLq;
mIkq2 = (mLmq*mLlkq1*mPsiq + mLmq*mLl*mPsikq1 - (mLmq*(mLlkq1 + mLl) + mLl*mLlkq1)*mPsikq2) / detLq;
mI0 = -mPsi0 / mLl;
// Update mechanical rotor angle with respect to electrical angle
mThetaMech = mThetaMech + dt * (mOmMech * mBase_OmMech);
}
}
else if (numMethod == NumericalMethod::Trapezoidal_flux) {
mMechPower = mechPower / mNomPower;
mMechTorque = mMechPower / mOmMech;
mElecTorque = (mPsid*mIq - mPsiq*mId);
// Save last step values
Real mOmMech_hist = mOmMech;
Real mElecTorque_hist = mElecTorque;
Real mMechTorque_hist = mMechTorque;
// Euler step forward
mOmMech = mOmMech + dt * (1 / (2 * mH) * (mMechTorque - mElecTorque));
DPSMatrix A = mBase_OmElec*(mResistanceMat*mReactanceMat - mOmMech*mOmegaFluxMat);
DPSMatrix I = DPSMatrix::Identity(7, 7);
DPSMatrix Aux = I + (dt / 2) * A;
DPSMatrix Aux2 = I - (dt / 2) * A;
DPSMatrix InvAux = Aux2.inverse();
DPSMatrix Fluxes(7, 1);
Fluxes(0, 0) = mPsiq;
Fluxes(1, 0) = mPsid;
Fluxes(2, 0) = mPsi0;
Fluxes(3, 0) = mPsikq1;
Fluxes(4, 0) = mPsikq2;
Fluxes(5, 0) = mPsifd;
Fluxes(6, 0) = mPsikd;
DPSMatrix dqVoltages(7, 1);
dqVoltages(0, 0) = mVq;
dqVoltages(1, 0) = mVd;
dqVoltages(2, 0) = mV0;
dqVoltages(3, 0) = mVkq1;
dqVoltages(4, 0) = mVkq2;
dqVoltages(5, 0) = mVfd;
dqVoltages(6, 0) = mVkd;
Fluxes = InvAux*Aux*Fluxes + InvAux*dt*mBase_OmElec*I*dqVoltages;