Added Turbine Governor model

Examples/Cxx/SynchronGenerators/ExciterAndTurbine.cpp

+/** 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;
+}

Examples/Cxx/SynchronGenerators/SimpSynGenThreePhaseFault.cpp

@@ -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) {

Examples/Matlab/CompareSynGenResults.m

@@ -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

Examples/Matlab/compareDPresults_with_PLECS.m

@@ -2,7 +2,7 @@ clc
 clear
 %% read PLECS results
 
-Results_PLECS = csvread('../../vsa/Results/SynGenDqEmt_ABCFault_Simulink/Voltages_and_currents.csv'); 
+Results_PLECS = csvread('../../../vsa/Results/SynGenDqEmt_ABCFault_PLECS/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'); 
 %theta_PLECS = csvread('../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/theta.csv'); 
@@ -10,8 +10,8 @@ Results_PLECS = csvread('../../vsa/Results/SynGenDqEmt_ABCFault_Simulink/Voltage
 %% 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');
+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');
 compOffsetDP = (size(voltageDP,2) - 1) / 2;
 
@@ -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,7 +121,7 @@ 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')

Examples/Matlab/debug.log

+[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 

Source/Components/Exciter.cpp

@@ -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;
 

Source/Components/TurbineGovernor.cpp

+/** 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;
+
+}
+
+
+

Source/Components/TurbineGovernor.h

+/** 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.
+*
+* 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/>.
+*********************************************************************************/
+
+#pragma once
+
+#include "BaseComponent.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 TurbineGovernor {
+	protected:
+
+		// ### Steam Turbine Parameters ####
+
+		/// Time constant of main inlet volume and steam chest
+		Real mTa;
+		/// Time constant of reheater
+		Real mTb;
+		/// Time constant of cross over piping and LP inlet volumes
+		Real mTc;
+		/// Fraction of total turbine power generated by HP section
+		Real mFa;
+		/// Fraction of total turbine power generated by IP section
+		Real mFb;
+		/// Fraction of total turbine power generated by LP section
+		Real mFc;
+
+		// ### Regulator ####
+		/// Main droop
+		Real mK;
+
+		// ### Speed Relay and servo motor ####
+		/// Time constant of speed relay
+		Real mTsr;
+		/// Time constant of servo motor
+		Real mTsm;
+		/// Opening rate limit
+		Real mLc1 = 0.1;
+		/// Closing rate limit
+		Real mLc2 = -1;
+
+		// ### Vaariables ###
+		/// Mechanical Torque in pu (output of steam turbine)
+		Real mTm = 0;
+		/// Valve position
+		Real mVcv = 0;
+		/// Valve position changing rate
+		Real mpVcv = 0;
+		/// Boiler pressure
+		Real mPb = 1;
+		/// Speed reference
+		Real mOmRef = 1;
+		/// Speed
+		Real mOm = 0;
+		/// Power Reference
+		Real mPmRef;
+		/// Input of speed realy
+		Real Psr_in = 0;
+		/// Input of servor motor
+		Real Psm_in = 0;
+
+		
+		
+
+		bool mLogActive;
+		Logger* mLog;
+
+
+	public:
+		TurbineGovernor() { };
+		~TurbineGovernor() {};
+
+		/// Initializes exciter parameters
+		TurbineGovernor(Real Ta, Real Tb, Real Tc, Real Fa, Real Fb, Real Fc, Real K, Real Tsr, Real Tsm, Real Tm_init);
+
+		/// Performs an step to update field voltage value
+		Real step(Real mOm, Real mOmRef, Real PmRef, Real dt);
+
+
+
+
+	};
+}

Source/Components/VoltageBehindReactanceEMT.cpp

@@ -48,6 +48,11 @@ void VoltageBehindReactanceEMT::AddExciter(Real Ta, Real Ka, Real Te, Real Ke, R
 	WithExciter = true;
 }
 
+void VoltageBehindReactanceEMT::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 VoltageBehindReactanceEMT::init(Real om, Real dt,
@@ -110,6 +115,7 @@ void VoltageBehindReactanceEMT::init(Real om, Real dt,
 	// steady state per unit initial value
 	initStatesInPerUnit(initActivePower, initReactivePower, initTerminalVolt, initVoltAngle, initFieldVoltage, initMechPower);
 
+
 	// Correcting variables
 	mThetaMech = mThetaMech + PI/2;
 	mMechTorque = -mMechTorque;
@@ -191,6 +197,11 @@ void VoltageBehindReactanceEMT::stepInPerUnit(Real om, Real dt, Real time, Numer
 		mIc;
 
 	// 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* mBase_OmMech);

Source/Components/VoltageBehindReactanceEMT.h

@@ -25,6 +25,7 @@
 
 #include "SynchGenBase.h"
 #include "Exciter.h"
+#include "TurbineGovernor.h"
 
 namespace DPsim {
 
@@ -42,6 +43,10 @@ namespace DPsim {
 		Exciter mExciter;