Commit c0e49ff5 authored by Markus Mirz's avatar Markus Mirz
Browse files

Merge branch 'development' of git.rwth-aachen.de:acs/core/simulation/dpsim into development

merge


Former-commit-id: 4f3fdf0b
parents 27edef64 47bcce2a
2b4d04e670b37372df0d3056b8be7bb2d7fe7e96
\ No newline at end of file
79462c4b18037d1e60faeefde281bab314188107
\ No newline at end of file
/** Synchron Generator Tests
*
* @file
* @author Markus Mirz <mmirz@eonerc.rwth-aachen.de>
* @copyright 2017, Institute for Automation of Complex Power Systems, EONERC
* @license GNU General Public License (version 3)
*
* DPsim
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*********************************************************************************/
#include "Simulation.h"
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"),
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<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<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);
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<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);
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::DynPhasor, 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<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);
// 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;
}
\ No newline at end of file
......@@ -66,13 +66,13 @@ int main() {
Real Llkq2 = 0;
// Declare circuit components
ElementPtr gen = std::make_shared<VoltageBehindReactanceEMT>("gen", 1, 2, 3,
ElementPtr gen = std::make_shared<SynchronGeneratorDP>("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);
......@@ -81,25 +81,25 @@ int main() {
circElements.push_back(r3);
// Declare circuit components for resistance change
Real breakerRes = 103.78378;
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);
Real breakerRes = 0.001;
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);
circElementsBreakerOn.push_back(rBreaker2);
circElementsBreakerOn.push_back(rBreaker3);
//circElementsBreakerOn.push_back(r1);
//circElementsBreakerOn.push_back(r2);
//circElementsBreakerOn.push_back(r3);
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);
tf = 0.3; dt = 0.000001; t = 0;
Int downSampling = 50;
Simulation newSim(circElements, om, dt, tf, log, SimulationType::DynPhasor, downSampling);
newSim.setNumericalMethod(NumericalMethod::Trapezoidal_flux);
newSim.addSystemTopology(circElementsBreakerOn);
newSim.switchSystemMatrix(0);
......@@ -111,9 +111,9 @@ 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<SynchronGeneratorDP> genPtr = std::dynamic_pointer_cast<SynchronGeneratorDP>(gen);
genPtr->init(om, dt, initActivePower, initReactivePower, initTerminalVolt, initVoltAngle, fieldVoltage, mechPower);
genPtr->AddExciter(Ta, Ka, Te, Ke, Tf, Kf, Tr, Lmd, Rfd);
//genPtr->AddExciter(Ta, Ka, Te, Ke, Tf, Kf, Tr, Lmd, Rfd);
// Calculate initial values for circuit at generator connection point
Real initApparentPower = sqrt(pow(initActivePower, 2) + pow(initReactivePower, 2));
......@@ -129,8 +129,8 @@ int main() {
Real lastLogTime = 0;
Real logTimeStep = 0.0001;
newSim.setSwitchTime(2, 1);
newSim.setSwitchTime(6, 0);
newSim.setSwitchTime(0.1, 1);
newSim.setSwitchTime(0.2, 0);
// Main Simulation Loop
while (newSim.getTime() < tf) {
......
......@@ -3,14 +3,14 @@
%% read PLECS results
Results_PLECS = csvread('../../../vsa/Results/SimpSynGenDqEmtExciter_LoadChange_Simulink/Voltages_and_currents.csv');
Te_PLECS = csvread('../../../vsa/Results/SimpSynGenDqEmtExciter_LoadChange_Simulink/electrical_torque.csv');
omega_PLECS = csvread('../../../vsa/Results/SimpSynGenDqEmtExciter_LoadChange_Simulink/omega.csv');
Results_PLECS = csvread('../../../vsa/Results/TestExciterAndTurbine/SynGenDqEMTExciter_loadChange_Simulink/Voltages_and_currents.csv');
%Te_PLECS = csvread('../../../vsa/Results/SynGenDqEmt_ABCFault_PLECS/electrical_torque.csv');
omega_PLECS = csvread('../../../vsa/Results/TestExciterAndTurbine/SynGenDqEMTExciter_loadChange_Simulink/omega.csv');
%theta_PLECS = csvread('../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/theta.csv');
%% read results from c++ simulation
VoltageVector = csvread('../../../vsa/Results/Testing/data_vt.csv');
CurrentVector = csvread('../../../vsa/Results/Testing/data_j.csv');
Log_SynGen = csvread('../../../vsa/Results/Testing/SynGen_gen.csv');
VoltageVector = csvread('../../../vsa/Results/TestExciterAndTurbine/SynGenVBREmtExciter_LoadChange_DPsim/data_vt.csv');
CurrentVector = csvread('../../../vsa/Results/TestExciterAndTurbine/SynGenVBREmtExciter_LoadChange_DPsim/data_j.csv');
Log_SynGen = csvread('../../../vsa/Results/TestExciterAndTurbine/SynGenVBREmtExciter_LoadChange_DPsim/SynGen_gen.csv');
%% Plot
figure(1)
hold off
......@@ -71,21 +71,21 @@ legend('ia DPSim','ia Simulink');
figure(3)
hold off
plot(Log_SynGen(:,1),Log_SynGen(:,8));
plot(Log_SynGen(:,1),Log_SynGen(:,9));
hold on
plot(Results_PLECS(:,1),omega_PLECS);
title('Rotor speed');
legend('\omega DPSim','\omega PLECS');
legend('\omega DPSim','\omega Simulink');
% %
% figure(4)
% hold off
% plot(Log_SynGen(:,1),Log_SynGen(:,20));
% hold on
% plot(Results_PLECS(:,1),-Te_PLECS);
%
figure(4)
hold off
plot(Log_SynGen(:,1),Log_SynGen(:,7));
hold on
plot(Results_PLECS(:,1),-Te_PLECS);
title('Electrical Torque');
legend('Te DPSim','Te PLECS');
% title('Electrical Torque');
% legend('Te DPSim','Te PLECS');
%
% figure(5)
% hold off
......
......@@ -2,17 +2,17 @@ clc
clear
%% read PLECS results
Results_PLECS = csvread('../../vsa/Results/SynGenDqEmt_ABCFault_Simulink/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/SynGenVBRDynPh_ABCFault_DPsim/data_vt.csv');
currentDP = csvread('../../vsa/Results/SynGenVBRDynPh_ABCFault_DPsim/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
......@@ -128,14 +128,14 @@ 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
......
[1207/155452:WARNING:resource_bundle.cc(291)] locale_file_path.empty() for locale
[1207/155452:WARNING:resource_bundle.cc(291)] locale_file_path.empty() for locale
......@@ -38,6 +38,7 @@ Exciter::Exciter(Real Ta, Real Ka, Real Te, Real Ke, Real Tf, Real Kf, Real Tr,
mLad = Lad;
mRfd = Rfd;
mVf = Vf_init*(mLad / mRfd);
mVf = 0;
}
......@@ -47,15 +48,15 @@ Real Exciter::step(Real mVd, Real mVq, Real Vref, Real dt) {
mVh = sqrt(pow(mVd, 2.) + pow(mVq, 2.));
// Voltage Transducer equation
mVm = Trapezoidal(mVm, -1, 1, dt / mTr, mVh);
mVm = Euler(mVm, -1, 1, dt / mTr, mVh);
// Stabilizing feedback equation
mVfl = mVfl + dt*mVis / mTf;
mVis = (mKf / mTf)*mVf - mVfl;
// Amplifier equation
mVr = Trapezoidal(mVr, -1, mKa, dt / mTa, Vref - mVm - mVis);
mVr = Euler(mVr, -1, mKa, dt / mTa, Vref - mVm - mVis);
// Exciter
mVse = 0.0039*exp(mVf*1.555);
mVf = Trapezoidal(mVf, -mKe, 1, dt / mTe, mVr - mVse);
mVf = Euler(mVf, -mKe, 1, dt / mTe, mVr - mVse);
return (mRfd / mLad)*mVf;
......
/** Turbine Governor
*
* @author Markus Mirz <mmirz@eonerc.rwth-aachen.de>
* @copyright 2017, Institute for Automation of Complex Power Systems, EONERC
* @license GNU General Public License (version 3)
*
* DPsim
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*********************************************************************************/
#include "TurbineGovernor.h"
#include "../IntegrationMethod.h"
using namespace DPsim;
TurbineGovernor::TurbineGovernor(Real Ta, Real Tb, Real Tc, Real Fa, Real Fb, Real Fc, Real K, Real Tsr, Real Tsm, Real Tm_init)
{
mTa = Ta;
mTb = Tb;
mTc = Tc;
mFa = Fa;
mFb = Fb;
mFc = Fc;
mK = K;
mTsr = Tsr;
mTsm = Tsm;
mTm = Tm_init;
}
Real TurbineGovernor::step(Real Om, Real OmRef, Real PmRef, Real dt) {
// ### Governing ###
// Input of speed relay
Psr_in = PmRef + (OmRef - Om)*mK;
// Input of servor motor
Psm_in = Euler(Psm_in, -1, 1, dt / mTsr, Psr_in);
// rate of change of valve
mpVcv = (Psm_in - mVcv) / mTsm;
if (mpVcv >= 0.1)
mpVcv = 0.1;
else if (mpVcv <= -1)
mpVcv = -1;
//Valve position
mVcv = mVcv + dt*mpVcv;
if (mVcv >= 1)
mVcv = 1;
else if (mVcv <= 0)
mVcv = 0;
//### Turbine ###
// Simplified equation
mTm = Euler(mTm, -1 / mTb, 1 / mTb, mpVcv*mFa, dt, mVcv);
return mTm;
}
/** Turbine Governor
*
* @file
* @author Markus Mirz <mmirz@eonerc.rwth-aachen.de>
* @copyright 2017, Institute for Automation of Complex Power Systems, EONERC
* @license GNU General Public License (version 3)
*
* DPsim
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* any later version.