Commit 47bcce2a authored by Viviane Sapucaia's avatar Viviane Sapucaia
Browse files

implemented exciter and turbine to VBR DP model


Former-commit-id: 41f7539f
parent 36b60faa
......@@ -25,6 +25,139 @@
using namespace DPsim;
//int main() {
// // Define Object for saving data on a file
// Logger log("log.txt"),
// vtLog("data_vt.csv"),
// jLog("data_j.csv");
//
// // Define machine parameters in per unit
// Real nomPower = 555e6;
// Real nomPhPhVoltRMS = 24e3;
// Real nomFreq = 60;
// Real nomFieldCurr = 1300;
// Int poleNum = 2;
// Real J = 2.8898e+04;
// Real H = 3.7;
//
// //Exciter
// Real Ka = 20;
// Real Ta = 0.2;
// Real Ke = 1;
// Real Te = 0.314;
// Real Kf = 0.063;
// Real Tf = 0.35;
// Real Tr = 0.02;
//
// // Turbine
// Real Ta_t = 0.3;
// Real Fa = 0.3;
// Real Tb = 7;
// Real Fb = 0.3;
// Real Tc = 0.2;
// Real Fc = 0.4;
// Real Tsr = 0.1;
// Real Tsm = 0.3;
// Real Kg = 20;
//
//
// Real Rs = 0.003;
// Real Ll = 0.15;
// Real Lmd = 1.6599;
// Real Lmd0 = 1.6599;
// Real Lmq = 1.61;
// Real Lmq0 = 1.61;
// Real Rfd = 0.0006;
// Real Llfd = 0.1648;
// Real Rkd = 0.0284;
// Real Llkd = 0.1713;
// Real Rkq1 = 0.0062;
// Real Llkq1 = 0.7252;
// Real Rkq2 = 0.0237;
// Real Llkq2 = 0.125;
// //Real Rkq2 = 0;
// //Real Llkq2 = 0;
//
// // Declare circuit components
// ElementPtr gen = std::make_shared<VoltageBehindReactanceEMT>("gen", 1, 2, 3,
// nomPower, nomPhPhVoltRMS, nomFreq, poleNum, nomFieldCurr,
// Rs, Ll, Lmd, Lmd0, Lmq, Lmq0, Rfd, Llfd, Rkd, Llkd, Rkq1, Llkq1, Rkq2, Llkq2, H, true);
// Real loadRes = 1037.8378;
// ElementPtr r1 = std::make_shared<ResistorEMT>("r1", 1, 0, loadRes);
// ElementPtr r2 = std::make_shared<ResistorEMT>("r2", 2, 0, loadRes);
// ElementPtr r3 = std::make_shared<ResistorEMT>("r3", 3, 0, loadRes);
//
// ElementList circElements;
// circElements.push_back(gen);
// circElements.push_back(r1);
// circElements.push_back(r2);
// circElements.push_back(r3);
//
// // Declare circuit components for resistance change
// Real breakerRes = 1037.8378;
// ElementPtr rBreaker1 = std::make_shared<ResistorEMT>("rbreak1", 1, 0, breakerRes);
// ElementPtr rBreaker2 = std::make_shared<ResistorEMT>("rbreak2", 2, 0, breakerRes);
// ElementPtr rBreaker3 = std::make_shared<ResistorEMT>("rbreak3", 3, 0, breakerRes);
// ElementList circElementsBreakerOn;
// circElementsBreakerOn.push_back(gen);
// circElementsBreakerOn.push_back(rBreaker1);
// circElementsBreakerOn.push_back(rBreaker2);
// circElementsBreakerOn.push_back(rBreaker3);
// circElementsBreakerOn.push_back(r1);
// circElementsBreakerOn.push_back(r2);
// circElementsBreakerOn.push_back(r3);
//
// // Set up simulation
// Real tf, dt, t;
// Real om = 2.0*M_PI*60.0;
// tf = 10; dt = 0.0001; t = 0;
// Int downSampling = 1;
// Simulation newSim(circElements, om, dt, tf, log, SimulationType::EMT, downSampling);
// newSim.setNumericalMethod(NumericalMethod::Trapezoidal_flux);
// newSim.addSystemTopology(circElementsBreakerOn);
// newSim.switchSystemMatrix(0);
//
// // Initialize generator
// Real initActivePower = 555e3;
// Real initReactivePower = 0;
// Real initTerminalVolt = 24000 / sqrt(3) * sqrt(2);
// Real initVoltAngle = -DPS_PI / 2;
// Real fieldVoltage = 7.0821;
// Real mechPower = 5.5558e5;
// shared_ptr<VoltageBehindReactanceEMT> genPtr = std::dynamic_pointer_cast<VoltageBehindReactanceEMT>(gen);
// genPtr->init(om, dt, initActivePower, initReactivePower, initTerminalVolt, initVoltAngle, fieldVoltage, mechPower);
// genPtr->AddExciter(Ta, Ka, Te, Ke, Tf, Kf, Tr, Lmd, Rfd);
// genPtr->AddGovernor(Ta_t, Tb, Tc, Fa, Fb, Fc, Kg, Tsr, Tsm, initActivePower/ nomPower);
//
// // Calculate initial values for circuit at generator connection point
// Real initApparentPower = sqrt(pow(initActivePower, 2) + pow(initReactivePower, 2));
// Real initTerminalCurr = initApparentPower / (3 * initTerminalVolt)* sqrt(2);
// Real initPowerFactor = acos(initActivePower / initApparentPower);
//
// std::cout << "A matrix:" << std::endl;
// std::cout << newSim.getSystemMatrix() << std::endl;
// std::cout << "vt vector:" << std::endl;
// std::cout << newSim.getLeftSideVector() << std::endl;
// std::cout << "j vector:" << std::endl;
// std::cout << newSim.getRightSideVector() << std::endl;
//
// Real lastLogTime = 0;
// Real logTimeStep = 0.0001;
// newSim.setSwitchTime(1, 1);
// //newSim.setSwitchTime(6, 0);
//
// // Main Simulation Loop
// while (newSim.getTime() < tf) {
// std::cout << newSim.getTime() << std::endl;
// newSim.stepGeneratorTest(vtLog, jLog, gen, newSim.getTime());
// newSim.increaseByTimeStep();
// }
//
// std::cout << "Simulation finished." << std::endl;
//
// return 0;
//}
int main() {
// Define Object for saving data on a file
Logger log("log.txt"),
......@@ -79,13 +212,13 @@ int main() {
//Real Llkq2 = 0;
// Declare circuit components
ElementPtr gen = std::make_shared<VoltageBehindReactanceEMT>("gen", 1, 2, 3,
ElementPtr gen = std::make_shared<VoltageBehindReactanceDP>("gen", 1, 2, 3,
nomPower, nomPhPhVoltRMS, nomFreq, poleNum, nomFieldCurr,
Rs, Ll, Lmd, Lmd0, Lmq, Lmq0, Rfd, Llfd, Rkd, Llkd, Rkq1, Llkq1, Rkq2, Llkq2, H, true);
Real loadRes = 1037.8378;
ElementPtr r1 = std::make_shared<ResistorEMT>("r1", 1, 0, loadRes);
ElementPtr r2 = std::make_shared<ResistorEMT>("r2", 2, 0, loadRes);
ElementPtr r3 = std::make_shared<ResistorEMT>("r3", 3, 0, loadRes);
ElementPtr r1 = std::make_shared<ResistorDP>("r1", 1, 0, loadRes);
ElementPtr r2 = std::make_shared<ResistorDP>("r2", 2, 0, loadRes);
ElementPtr r3 = std::make_shared<ResistorDP>("r3", 3, 0, loadRes);
ElementList circElements;
circElements.push_back(gen);
......@@ -95,9 +228,9 @@ int main() {
// Declare circuit components for resistance change
Real breakerRes = 1037.8378;
ElementPtr rBreaker1 = std::make_shared<ResistorEMT>("rbreak1", 1, 0, breakerRes);
ElementPtr rBreaker2 = std::make_shared<ResistorEMT>("rbreak2", 2, 0, breakerRes);
ElementPtr rBreaker3 = std::make_shared<ResistorEMT>("rbreak3", 3, 0, breakerRes);
ElementPtr rBreaker1 = std::make_shared<ResistorDP>("rbreak1", 1, 0, breakerRes);
ElementPtr rBreaker2 = std::make_shared<ResistorDP>("rbreak2", 2, 0, breakerRes);
ElementPtr rBreaker3 = std::make_shared<ResistorDP>("rbreak3", 3, 0, breakerRes);
ElementList circElementsBreakerOn;
circElementsBreakerOn.push_back(gen);
circElementsBreakerOn.push_back(rBreaker1);
......@@ -112,7 +245,7 @@ int main() {
Real om = 2.0*M_PI*60.0;
tf = 10; dt = 0.0001; t = 0;
Int downSampling = 1;
Simulation newSim(circElements, om, dt, tf, log, SimulationType::EMT, downSampling);
Simulation newSim(circElements, om, dt, tf, log, SimulationType::DynPhasor, downSampling);
newSim.setNumericalMethod(NumericalMethod::Trapezoidal_flux);
newSim.addSystemTopology(circElementsBreakerOn);
newSim.switchSystemMatrix(0);
......@@ -124,10 +257,10 @@ int main() {
Real initVoltAngle = -DPS_PI / 2;
Real fieldVoltage = 7.0821;
Real mechPower = 5.5558e5;
shared_ptr<VoltageBehindReactanceEMT> genPtr = std::dynamic_pointer_cast<VoltageBehindReactanceEMT>(gen);
shared_ptr<VoltageBehindReactanceDP> genPtr = std::dynamic_pointer_cast<VoltageBehindReactanceDP>(gen);
genPtr->init(om, dt, initActivePower, initReactivePower, initTerminalVolt, initVoltAngle, fieldVoltage, mechPower);
genPtr->AddExciter(Ta, Ka, Te, Ke, Tf, Kf, Tr, Lmd, Rfd);
genPtr->AddGovernor(Ta_t, Tb, Tc, Fa, Fb, Fc, Kg, Tsr, Tsm, initActivePower/ nomPower);
genPtr->AddGovernor(Ta_t, Tb, Tc, Fa, Fb, Fc, Kg, Tsr, Tsm, initActivePower / nomPower);
// Calculate initial values for circuit at generator connection point
Real initApparentPower = sqrt(pow(initActivePower, 2) + pow(initReactivePower, 2));
......@@ -156,4 +289,4 @@ int main() {
std::cout << "Simulation finished." << std::endl;
return 0;
}
}
\ No newline at end of file
......@@ -2,17 +2,17 @@ clc
clear
%% read PLECS results
Results_PLECS = csvread('../../../vsa/Results/SynGenDqEmt_ABCFault_PLECS/Voltages_and_currents.csv');
Results_PLECS = csvread('../../../vsa/Results/TestExciterAndTurbine/SynGenDqEMTExciter_loadChange_Simulink/Voltages_and_currents.csv');
%Te_PLECS = csvread('../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/electrical_torque.csv');
%omega_PLECS = csvread('../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/omega.csv');
omega_PLECS = csvread('../../../vsa/Results/TestExciterAndTurbine/SynGenDqEMTExciter_loadChange_Simulink/omega.csv');
%theta_PLECS = csvread('../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/theta.csv');
%% Read data from DP simulation and calculate absolute value and phase
% Read values from CSV files
voltageDP = csvread('../../../vsa/Results/Testing/data_vt.csv');
currentDP = csvread('../../../vsa/Results/Testing/data_j.csv');
%Log_SynGen = csvread('../../vsa/Results//SynGenDqDynPh_ABCFault_DPsim/Euler/SynGen_gen.csv');
voltageDP = csvread('../../../vsa/Results/TestExciterAndTurbine/SynGenVBRDPExciter_LoadChange_DPsim/data_vt.csv');
currentDP = csvread('../../../vsa/Results/TestExciterAndTurbine/SynGenVBRDPExciter_LoadChange_DPsim/data_j.csv');
Log_SynGen = csvread('../../../vsa/Results/TestExciterAndTurbine/SynGenVBRDPExciter_LoadChange_DPsim/SynGen_gen.csv');
compOffsetDP = (size(voltageDP,2) - 1) / 2;
% Calculate Voltage DP absolute value
......@@ -97,7 +97,7 @@ figure(4)
hold off
PLECSplotc = plot(Results_PLECS(:,1), Results_PLECS(:,5), '--');
hold on
DPplotc = plot(currentShiftDP(:,1),-currentShiftDP(:,2));
DPplotc = plot(currentShiftDP(:,1),currentShiftDP(:,2));
DPabsPlotc = plot(currentAbsDP(:,1),currentAbsDP(:,2));
title('Current phase A');
legend('Current Phase a Simulink', 'DP shift a', 'DP abs a')
......@@ -109,7 +109,7 @@ figure(5)
hold off
PLECSplot2c = plot(Results_PLECS(:,1), Results_PLECS(:,6), '--');
hold on
DPplot2c = plot(currentShiftDP(:,1),-currentShiftDP(:,3));
DPplot2c = plot(currentShiftDP(:,1),currentShiftDP(:,3));
DPabsPlot2c = plot(currentAbsDP(:,1),currentAbsDP(:,3));
title('Current phase B');
legend('Current Phase b Simulink', 'DP shift b', 'DP abs b')
......@@ -121,21 +121,21 @@ figure(6)
hold off
PLECSplot3c = plot(Results_PLECS(:,1), Results_PLECS(:,7), '--');
hold on
DPplot3c = plot(currentShiftDP(:,1),-currentShiftDP(:,4));
DPplot3c = plot(currentShiftDP(:,1),currentShiftDP(:,4));
DPabsPlot3c = plot(currentAbsDP(:,1),currentAbsDP(:,4));
title('Currents phase C');
legend('Current Phase c Simulink', 'DP shift c', 'DP abs c')
xlabel('time [s]')
ylabel('current [A]')
% figure(7)
% hold off
% plot(Log_SynGen(:,1),Log_SynGen(:,21));
% hold on
% plot(Results_PLECS(:,1),omega_PLECS);
%
% title('Rotor speed');
% legend('\omega DPSim','\omega PLECS');
figure(7)
hold off
plot(Log_SynGen(:,1),Log_SynGen(:,9));
hold on
plot(Results_PLECS(:,1),omega_PLECS);
title('Rotor speed');
legend('\omega DPSim','\omega PLECS');
%
% figure(8)
% hold off
......
......@@ -43,6 +43,17 @@ VoltageBehindReactanceDP::~VoltageBehindReactanceDP() {
}
void VoltageBehindReactanceDP::AddExciter(Real Ta, Real Ka, Real Te, Real Ke, Real Tf, Real Kf, Real Tr, Real Lad, Real Rfd) {
mExciter = Exciter(Ta, Ka, Te, Ke, Tf, Kf, Tr, Lad, Rfd, mVfd);
WithExciter = true;
}
void VoltageBehindReactanceDP::AddGovernor(Real Ta, Real Tb, Real Tc, Real Fa, Real Fb, Real Fc, Real K, Real Tsr, Real Tsm, Real Tm_init) {
mTurbineGovernor = TurbineGovernor(Ta, Tb, Tc, Fa, Fb, Fc, K, Tsr, Tsm, Tm_init);
WithTurbineGovernor = true;
}
void VoltageBehindReactanceDP::init(Real om, Real dt,
Real initActivePower, Real initReactivePower, Real initTerminalVolt, Real initVoltAngle, Real initFieldVoltage, Real initMechPower) {
......@@ -169,6 +180,8 @@ void VoltageBehindReactanceDP::step(SystemModel& system, Real time) {
stepInPerUnit(system.getOmega(), system.getTimeStep(), time, system.getNumMethod());
R_load = system.getCurrentSystemMatrix().inverse() / mBase_Z;
// Update current source accordingly
if (mNode1 >= 0) {
system.addCompToRightSideVector(mNode1, -mIaRe*mBase_i, -mIaIm*mBase_i);
......@@ -199,6 +212,11 @@ void VoltageBehindReactanceDP::stepInPerUnit(Real om, Real dt, Real time, Numeri
mIcIm;
// Calculate mechanical variables with euler
if (WithTurbineGovernor == true)
{
mMechTorque = -mTurbineGovernor.step(mOmMech, 1, 0.001, dt);
}
mElecTorque = (mPsimd*mIq - mPsimq*mId);
mOmMech = mOmMech + dt * (1. / (2. * mH) * (mElecTorque - mMechTorque));
mThetaMech = mThetaMech + dt * ((mOmMech - 1) * mBase_OmMech);
......@@ -225,6 +243,14 @@ void VoltageBehindReactanceDP::stepInPerUnit(Real om, Real dt, Real time, Numeri
mId = abcToDq0Transform(mThetaMech, mIaRe, mIbRe, mIcRe, mIaIm, mIbIm, mIcIm)(1);
mI0 = abcToDq0Transform(mThetaMech, mIaRe, mIbRe, mIcRe, mIaIm, mIbIm, mIcIm)(2);
if (WithExciter == true) {
// dq-transform of interface voltage
mVd = R_load(0, 0)*mId;
mVq = R_load(0, 0)*mIq;
mVfd = mExciter.step(mVd, mVq, 1, dt);
}
// Calculate rotor flux likanges
if (mNumDampingWindings == 2)
{
......@@ -327,27 +353,27 @@ void VoltageBehindReactanceDP::stepInPerUnit(Real om, Real dt, Real time, Numeri
mDVbIm,
mDVcIm;
// Load resistance
if (time < 0.1 || time > 0.2)
{
R_load <<
1037.8378 / mBase_Z, 0, 0, 0, 0, 0,
0, 1037.8378 / mBase_Z, 0, 0, 0, 0,
0, 0, 1037.8378 / mBase_Z, 0, 0, 0,
0, 0, 0, 1037.8378 / mBase_Z, 0, 0,
0, 0, 0, 0, 1037.8378 / mBase_Z, 0,
0, 0, 0, 0, 0, 1037.8378 / mBase_Z;
}
else
{
R_load <<
0.001 / mBase_Z, 0, 0, 0, 0, 0,
0, 0.001 / mBase_Z, 0, 0, 0, 0,
0, 0, 0.001 / mBase_Z, 0, 0, 0,
0, 0, 0, 0.001 / mBase_Z, 0, 0,
0, 0, 0, 0, 0.001 / mBase_Z, 0,
0, 0, 0, 0, 0, 0.001 / mBase_Z;
}
//// Load resistance
//if (time < 0.1 || time > 0.2)
//{
// R_load <<
// 1037.8378 / mBase_Z, 0, 0, 0, 0, 0,
// 0, 1037.8378 / mBase_Z, 0, 0, 0, 0,
// 0, 0, 1037.8378 / mBase_Z, 0, 0, 0,
// 0, 0, 0, 1037.8378 / mBase_Z, 0, 0,
// 0, 0, 0, 0, 1037.8378 / mBase_Z, 0,
// 0, 0, 0, 0, 0, 1037.8378 / mBase_Z;
//}
//else
//{
// R_load <<
// 0.001 / mBase_Z, 0, 0, 0, 0, 0,
// 0, 0.001 / mBase_Z, 0, 0, 0, 0,
// 0, 0, 0.001 / mBase_Z, 0, 0, 0,
// 0, 0, 0, 0.001 / mBase_Z, 0, 0,
// 0, 0, 0, 0, 0.001 / mBase_Z, 0,
// 0, 0, 0, 0, 0, 0.001 / mBase_Z;
//}
}
......
......@@ -24,6 +24,8 @@
#pragma once
#include "SynchGenBase.h"
#include "Exciter.h"
#include "TurbineGovernor.h"
namespace DPsim {
......@@ -37,6 +39,16 @@ namespace DPsim {
class VoltageBehindReactanceDP : public SynchGenBase {
protected:
/// Exciter Model
Exciter mExciter;
/// Determine if Exciter is activated
bool WithExciter = false;
/// Governor Model
TurbineGovernor mTurbineGovernor;
/// Determine if Turbine and Governor are activated
bool WithTurbineGovernor = false;
/// d dynamic inductance
Real mDLmd;
/// q dynamic inductance
......@@ -154,6 +166,11 @@ namespace DPsim {
Real Rkq1, Real Llkq1, Real Rkq2, Real Llkq2,
Real inertia, bool logActive = false);
/// Function to initialize Exciter
void AddExciter(Real Ta, Real Ka, Real Te, Real Ke, Real Tf, Real Kf, Real Tr, Real Lad, Real Rfd);
/// Function to initialize Governor and Turbine
void AddGovernor(Real Ta, Real Tb, Real Tc, Real Fa, Real Fb, Real Fc, Real K, Real Tsr, Real Tsm, Real Tm_init);
/// 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,
......
......@@ -41,10 +41,12 @@ namespace DPsim {
/// Exciter Model
Exciter mExciter;
/// Determine if Exciter is activated
bool WithExciter = false;
/// Governor Model
TurbineGovernor mTurbineGovernor;
/// Determine if Turbine and Governor are activated
bool WithTurbineGovernor = false;
/// d dynamic inductance
......@@ -152,7 +154,9 @@ namespace DPsim {
Real Rkq1, Real Llkq1, Real Rkq2, Real Llkq2,
Real inertia, bool logActive = false);
/// Function to initialize Exciter
void AddExciter(Real Ta, Real Ka, Real Te, Real Ke, Real Tf, Real Kf, Real Tr, Real Lad, Real Rfd);
/// Function to initialize Governor and Turbine
void AddGovernor(Real Ta, Real Tb, Real Tc, Real Fa, Real Fb, Real Fc, Real K, Real Tsr, Real Tsm, Real Tm_init);
/// Initializes states in per unit or stator referred variables depending on the setting of the state type.
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment