Commit 49184cc1 authored by Viviane Sapucaia's avatar Viviane Sapucaia
Browse files

clearing fault when phase currents are equal to zero

parent 61bfc985
...@@ -60,10 +60,10 @@ int main() { ...@@ -60,10 +60,10 @@ int main() {
Real Llkd = 0.1713; Real Llkd = 0.1713;
Real Rkq1 = 0.0062; Real Rkq1 = 0.0062;
Real Llkq1 = 0.7252; Real Llkq1 = 0.7252;
Real Rkq2 = 0.0237; //Real Rkq2 = 0.0237;
Real Llkq2 = 0.125; //Real Llkq2 = 0.125;
//Real Rkq2 = 0; Real Rkq2 = 0;
//Real Llkq2 = 0; Real Llkq2 = 0;
// Declare circuit components // Declare circuit components
ElementPtr gen = std::make_shared<VoltageBehindReactanceEMT>("gen", 1, 2, 3, ElementPtr gen = std::make_shared<VoltageBehindReactanceEMT>("gen", 1, 2, 3,
...@@ -97,7 +97,7 @@ int main() { ...@@ -97,7 +97,7 @@ int main() {
// Set up simulation // Set up simulation
Real tf, dt, t; Real tf, dt, t;
Real om = 2.0*M_PI*60.0; Real om = 2.0*M_PI*60.0;
tf = 3; dt = 0.0001; t = 0; tf = 0.3; dt = 0.0001; t = 0;
Int downSampling = 1; Int downSampling = 1;
Simulation newSim(circElements, om, dt, tf, log, SimulationType::EMT, downSampling); Simulation newSim(circElements, om, dt, tf, log, SimulationType::EMT, downSampling);
newSim.setNumericalMethod(NumericalMethod::Trapezoidal_flux); newSim.setNumericalMethod(NumericalMethod::Trapezoidal_flux);
...@@ -128,9 +128,9 @@ int main() { ...@@ -128,9 +128,9 @@ int main() {
std::cout << newSim.getRightSideVector() << std::endl; std::cout << newSim.getRightSideVector() << std::endl;
Real lastLogTime = 0; Real lastLogTime = 0;
Real logTimeStep = 0.00005; Real logTimeStep = 0.0001;
newSim.setSwitchTime(0.1, 1); newSim.setSwitchTime(0.1, 1);
newSim.setSwitchTime(2.1, 0); newSim.setSwitchTime(0.2, 0);
// Main Simulation Loop // Main Simulation Loop
while (newSim.getTime() < tf) { while (newSim.getTime() < tf) {
......
...@@ -3,28 +3,25 @@ ...@@ -3,28 +3,25 @@
%% read PLECS results %% read PLECS results
Results_PLECS = csvread('../../vsa/Results/SynGenDqEmt_ABCFault_PLECS/Voltages_and_currents.csv'); Results_PLECS = csvread('../../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/Voltages_and_currents.csv');
Te_PLECS = csvread('../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/electrical_torque.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/SynGenVBREmt_ABCFault_PLECS/omega.csv');
%theta_PLECS = csvread('../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/theta.csv'); %theta_PLECS = csvread('../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/theta.csv');
%% read results from c++ simulation %% read results from c++ simulation
VoltageVector = csvread('../../vsa/Results/Testing/data_vt.csv'); VoltageVector = csvread('../../../vsa/Results/Testing/data_vt.csv');
CurrentVector = csvread('../../vsa/Results/Testing/data_j.csv'); CurrentVector = csvread('../../../vsa/Results/Testing/data_j.csv');
%%Log_SynGen = csvread('../../vsa/Results/Testing/SynGen_gen.csv'); Log_SynGen = csvread('../../../vsa/Results/Testing/SynGen_gen.csv');
%% Plot %% Plot
figure(1) figure(1)
hold off hold off
plot(VoltageVector(:,1),VoltageVector(:,2)); plot(VoltageVector(:,1),VoltageVector(:,2));
hold on hold on
% plot(VoltageVector(:,1),VoltageVector(:,3)); plot(VoltageVector(:,1),VoltageVector(:,3));
% plot(VoltageVector(:,1),VoltageVector(:,4)); plot(VoltageVector(:,1),VoltageVector(:,4));
plot(tout,voltages(:,1),'--') plot(Results_PLECS(:,1),Results_PLECS(:,2),'--');
% plot(tout,voltages(:,2),'--') plot(Results_PLECS(:,1),Results_PLECS(:,3),'--');
% plot(tout,voltages(:,3),'--') plot(Results_PLECS(:,1),Results_PLECS(:,4),'--');
% plot(Results_PLECS(:,1),Results_PLECS(:,2),'--');
% plot(Results_PLECS(:,1),Results_PLECS(:,3),'--');
% plot(Results_PLECS(:,1),Results_PLECS(:,4),'--');
title('Phase Voltages'); title('Phase Voltages');
legend('va DPSim','vb DPSim', 'vc DPSim','va PLECS','vb PLECS','vc PLECS'); legend('va DPSim','vb DPSim', 'vc DPSim','va PLECS','vb PLECS','vc PLECS');
...@@ -32,19 +29,14 @@ legend('va DPSim','vb DPSim', 'vc DPSim','va PLECS','vb PLECS','vc PLECS'); ...@@ -32,19 +29,14 @@ legend('va DPSim','vb DPSim', 'vc DPSim','va PLECS','vb PLECS','vc PLECS');
figure(2) figure(2)
hold off hold off
plot(CurrentVector(:,1),CurrentVector(:,2)); plot(CurrentVector(:,1),-CurrentVector(:,2));
hold on hold on
% plot(CurrentVector(:,1),CurrentVector(:,3)); plot(CurrentVector(:,1),-CurrentVector(:,3));
% plot(CurrentVector(:,1),CurrentVector(:,4)); plot(CurrentVector(:,1),-CurrentVector(:,4));
plot(tout,currents(:,1),'--') plot(Results_PLECS(:,1),Results_PLECS(:,5),'--');
% plot(tout,currents(:,2),'--') plot(Results_PLECS(:,1),Results_PLECS(:,6),'--');
% plot(tout,currents(:,3),'--') plot(Results_PLECS(:,1),Results_PLECS(:,7),'--');
% plot(Results_PLECS(:,1),Results_PLECS(:,5),'--');
% plot(Results_PLECS(:,1),Results_PLECS(:,6),'--');
% plot(Results_PLECS(:,1),Results_PLECS(:,7),'--');
title('Phase Currents'); title('Phase Currents');
legend('ia DPSim','ib DPSim','ic DPSim','ia PLECS','ib PLECS','ic PLECS'); legend('ia DPSim','ib DPSim','ic DPSim','ia PLECS','ib PLECS','ic PLECS');
...@@ -77,14 +69,14 @@ legend('ia DPSim','ib DPSim','ic DPSim','ia PLECS','ib PLECS','ic PLECS'); ...@@ -77,14 +69,14 @@ legend('ia DPSim','ib DPSim','ic DPSim','ia PLECS','ib PLECS','ic PLECS');
% title ('dq currents'); % title ('dq currents');
% legend('q','d','fd'); % legend('q','d','fd');
% figure(3) figure(3)
% hold off hold off
% plot(Log_SynGen(:,1),Log_SynGen(:,8)); plot(Log_SynGen(:,1),Log_SynGen(:,8));
% hold on hold on
% plot(Results_PLECS(:,1),omega_PLECS); plot(Results_PLECS(:,1),omega_PLECS);
%
% title('Rotor speed'); title('Rotor speed');
% legend('\omega DPSim','\omega PLECS'); legend('\omega DPSim','\omega PLECS');
% %
% figure(4) % figure(4)
% hold off % hold off
......
% Compare voltage and current of c++ simulation with voltage and currents
% from PLECS simulation
%% read results from c++ simulation Full model
VoltageVector_full = csvread('../../vsa/Results/SynGenVbrEmt_ABCLongFault_DPsim/data_vt.csv');
CurrentVector_full = csvread('../../vsa/Results/SynGenVbrEmt_ABCLongFault_DPsim/data_j.csv');
Log_SynGen_full = csvread('../../vsa/Results/SynGenVbrEmt_ABCLongFault_DPsim/SynGen_gen.csv');
%% read results from c++ simulation
VoltageVector = csvread('../../vsa/Results/SimpSynGenVbrEmt_ABCLongFault_DPsim/data_vt.csv');
CurrentVector = csvread('../../vsa/Results/SimpSynGenVbrEmt_ABCLongFault_DPsim/data_j.csv');
Log_SynGen = csvread('../../vsa/Results/SimpSynGenVbrEmt_ABCLongFault_DPsim/SynGen_gen.csv');
%% Plot
figure(1)
hold off
plot(VoltageVector(:,1),VoltageVector(:,2));
hold on
plot(VoltageVector_full(:,1),VoltageVector_full(:,2));
title('Voltage Phase a');
legend('va DPSim simplified model','va DPSim detailed model');
figure(2)
hold off
plot(VoltageVector(:,1),VoltageVector(:,3));
hold on
plot(VoltageVector_full(:,1),VoltageVector_full(:,3));
%plot(Results_Simulink(:,1),Results_Simulink(:,3),'--');
title('Voltage Phase b');
legend('vb DPSim simplified model','vb DPSim detailed model');
figure(3)
hold off
plot(VoltageVector(:,1),VoltageVector(:,4));
hold on
plot(VoltageVector_full(:,1),VoltageVector_full(:,4));
title('Voltage Phase c');
legend('vc DPSim simplified model','vc DPSim detailed model');
figure(4)
hold off
plot(CurrentVector(:,1),CurrentVector(:,2));
hold on
plot(CurrentVector_full(:,1),CurrentVector_full(:,2));
title('Current Phase a');
legend('ia DPSim simplified model','ia DPSim detailed model');
figure(5)
hold off
plot(CurrentVector(:,1),CurrentVector(:,3));
hold on
plot(CurrentVector_full(:,1),CurrentVector_full(:,3),'--');
title('Current Phase b');
legend('ib DPSim simplified model','ib DPSim detailed model');
figure(6)
hold off
plot(CurrentVector(:,1),CurrentVector(:,4));
hold on
plot(CurrentVector_full(:,1),CurrentVector_full(:,4),'--');
title('Current Phase c');
legend('ic DPSim simplified model','ic DPSim detailed model');
figure(7)
hold off
plot(Log_SynGen(:,1),Log_SynGen(:,9));
hold on
plot(Log_SynGen_full(:,1),Log_SynGen_full(:,9));
title('Rotor speed');
legend('\omega DPSim simplified model','\omega DPSim detailed model');
figure(8)
hold off
plot(Log_SynGen(:,1),Log_SynGen(:,8));
hold on
plot(Log_SynGen_full(:,1),Log_SynGen_full(:,8));
title('Electrical Torque');
legend('Te DPSim simplified model','Te DPSim detailed model');
...@@ -164,6 +164,10 @@ void SimplifiedVBR::stepInPerUnit(Real om, Real dt, Real time, NumericalMethod n ...@@ -164,6 +164,10 @@ void SimplifiedVBR::stepInPerUnit(Real om, Real dt, Real time, NumericalMethod n
mIq = mDqStatorCurrents(0, 0); mIq = mDqStatorCurrents(0, 0);
mId = mDqStatorCurrents(1, 0); mId = mDqStatorCurrents(1, 0);
mIa_hist = mIa;
mIb_hist = mIb;
mIc_hist = mIc;
mIa = inverseParkTransform(mThetaMech, mIq, mId, 0)(0); mIa = inverseParkTransform(mThetaMech, mIq, mId, 0)(0);
mIb = inverseParkTransform(mThetaMech, mIq, mId, 0)(1); mIb = inverseParkTransform(mThetaMech, mIq, mId, 0)(1);
mIc = inverseParkTransform(mThetaMech, mIq, mId, 0)(2); mIc = inverseParkTransform(mThetaMech, mIq, mId, 0)(2);
...@@ -212,19 +216,42 @@ void SimplifiedVBR::stepInPerUnit(Real om, Real dt, Real time, NumericalMethod n ...@@ -212,19 +216,42 @@ void SimplifiedVBR::stepInPerUnit(Real om, Real dt, Real time, NumericalMethod n
// Load resistance // Load resistance
if (time < 0.1 ||time > 2.1 ) if (time < 0.1)
{ {
Ra = 1037.8378 / mBase_Z;
Rb = 1037.8378 / mBase_Z;
Rc = 1037.8378 / mBase_Z;
R_load << R_load <<
1037.8378 / mBase_Z, 0, 0, Ra, 0, 0,
0, 1037.8378 / mBase_Z, 0, 0, Rb, 0,
0, 0, 1037.8378 / mBase_Z; 0, 0, Rc;
}
else if (time > 0.2)
{
if (signbit(mIa) != signbit(mIa_hist))
Ra = 1037.8378 / mBase_Z;
if (signbit(mIb) != signbit(mIb_hist))
Rb = 1037.8378 / mBase_Z;
if (signbit(mIc) != signbit(mIc_hist))
Rc = 1037.8378 / mBase_Z;
R_load <<
Ra, 0, 0,
0, Rb, 0,
0, 0, Rc;
} }
else else
{ {
Ra = 0.001 / mBase_Z;
Rb = 0.001 / mBase_Z;
Rc = 0.001 / mBase_Z;
R_load << R_load <<
0.001 / mBase_Z, 0, 0, Ra, 0, 0,
0, 0.001 / mBase_Z, 0, 0, Rb, 0,
0, 0, 0.001 / mBase_Z; 0, 0, Rc;
} }
} }
......
...@@ -55,6 +55,15 @@ namespace DPsim { ...@@ -55,6 +55,15 @@ namespace DPsim {
/// Auxiliar inductance /// Auxiliar inductance
Real mLb; Real mLb;
///Load Resistance
Real Ra;
Real Rb;
Real Rc;
Real mIa_hist;
Real mIb_hist;
Real mIc_hist;
/// d dynamic flux /// d dynamic flux
Real mDPsid; Real mDPsid;
......
...@@ -153,8 +153,13 @@ void VoltageBehindReactanceEMT::init(Real om, Real dt, ...@@ -153,8 +153,13 @@ void VoltageBehindReactanceEMT::init(Real om, Real dt,
void VoltageBehindReactanceEMT::step(SystemModel& system, Real time) { void VoltageBehindReactanceEMT::step(SystemModel& system, Real time) {
stepInPerUnit(system.getOmega(), system.getTimeStep(), time, system.getNumMethod()); stepInPerUnit(system.getOmega(), system.getTimeStep(), time, system.getNumMethod());
R_load = system.getCurrentSystemMatrix().inverse() / mBase_Z;
// Update current source accordingly // Update current source accordingly
if (mNode1 >= 0) { if (mNode1 >= 0) {
system.addRealToRightSideVector(mNode1, -mIa*mBase_i); system.addRealToRightSideVector(mNode1, -mIa*mBase_i);
...@@ -172,6 +177,7 @@ void VoltageBehindReactanceEMT::step(SystemModel& system, Real time) { ...@@ -172,6 +177,7 @@ void VoltageBehindReactanceEMT::step(SystemModel& system, Real time) {
mLog->LogDataLine(time, logValues); mLog->LogDataLine(time, logValues);
} }
} }
...@@ -187,15 +193,24 @@ void VoltageBehindReactanceEMT::stepInPerUnit(Real om, Real dt, Real time, Numer ...@@ -187,15 +193,24 @@ void VoltageBehindReactanceEMT::stepInPerUnit(Real om, Real dt, Real time, Numer
mOmMech = mOmMech + dt * (1. / (2. * mH) * (mElecTorque - mMechTorque)); mOmMech = mOmMech + dt * (1. / (2. * mH) * (mElecTorque - mMechTorque));
mThetaMech = mThetaMech + dt * (mOmMech* mBase_OmMech); mThetaMech = mThetaMech + dt * (mOmMech* mBase_OmMech);
// Calculate Inductance matrix and its derivative // Calculate Inductance matrix and its derivative
CalculateLandpL(); CalculateLandpL();
// Solve circuit - calculate stator currents // Solve circuit - calculate stator currents
if (numMethod == NumericalMethod::Trapezoidal_flux) if (numMethod == NumericalMethod::Trapezoidal_flux) {
mIabc = Trapezoidal(mIabc, -mDInductanceMat.inverse()*(mResistanceMat + R_load + pmDInductanceMat), mDInductanceMat.inverse(), dt*mBase_OmElec, -mDVabc, -mDVabc_hist); mIabc = Trapezoidal(mIabc, -mDInductanceMat.inverse()*(mResistanceMat + R_load + pmDInductanceMat), mDInductanceMat.inverse(), dt*mBase_OmElec, -mDVabc, -mDVabc_hist);
}
else if (numMethod == NumericalMethod::Euler) else if (numMethod == NumericalMethod::Euler)
mIabc = Euler(mIabc, -mDInductanceMat.inverse()*(mResistanceMat + R_load + pmDInductanceMat), mDInductanceMat.inverse(), dt*mBase_OmElec, -mDVabc); mIabc = Euler(mIabc, -mDInductanceMat.inverse()*(mResistanceMat + R_load + pmDInductanceMat), mDInductanceMat.inverse(), dt*mBase_OmElec, -mDVabc);
mIa_hist = mIa;
mIb_hist = mIb;
mIc_hist = mIc;
mIa = mIabc(0); mIa = mIabc(0);
mIb = mIabc(1); mIb = mIabc(1);
...@@ -303,21 +318,9 @@ void VoltageBehindReactanceEMT::stepInPerUnit(Real om, Real dt, Real time, Numer ...@@ -303,21 +318,9 @@ void VoltageBehindReactanceEMT::stepInPerUnit(Real om, Real dt, Real time, Numer
mDVb, mDVb,
mDVc; mDVc;
// Load resistance
if (time < 0.1 || time > 0.2)
{
R_load <<
1037.8378 / mBase_Z, 0, 0,
0, 1037.8378 / mBase_Z, 0,
0, 0, 1037.8378 / mBase_Z;
}
else
{
R_load <<
0.001 / mBase_Z, 0, 0,
0, 0.001 / mBase_Z, 0,
0, 0, 0.001 / mBase_Z;
}
} }
......
...@@ -77,6 +77,22 @@ namespace DPsim { ...@@ -77,6 +77,22 @@ namespace DPsim {
/// Interface curent phase c /// Interface curent phase c
Real mIc; Real mIc;
/// Interface curent phase a last time step
Real mIa_hist;
/// Interface curent phase b last time step
Real mIb_hist;
/// Interface curent phase c last time step
Real mIc_hist;
///Load Resistance phase a
Real Ra;
///Load Resistance phase b
Real Rb;
///Load Resistance phase c
Real Rc;
/// Magnetizing flux linkage in q axis /// Magnetizing flux linkage in q axis
Real mPsimq; Real mPsimq;
/// Magnetizing flux linkage in d axis /// Magnetizing flux linkage in d axis
......
...@@ -197,12 +197,26 @@ Int Simulation::stepGeneratorTest(Logger& leftSideVectorLog, Logger& rightSideVe ...@@ -197,12 +197,26 @@ Int Simulation::stepGeneratorTest(Logger& leftSideVectorLog, Logger& rightSideVe
(*it)->postStep(mSystemModel); (*it)->postStep(mSystemModel);
} }
if (ClearingFault)
clearFault(1, 2, 3);
if (mCurrentSwitchTimeIndex < mSwitchEventVector.size()) { if (mCurrentSwitchTimeIndex < mSwitchEventVector.size()) {
if (mTime >= mSwitchEventVector[mCurrentSwitchTimeIndex].switchTime) { if (mTime >= mSwitchEventVector[mCurrentSwitchTimeIndex].switchTime) {
switchSystemMatrix(mSwitchEventVector[mCurrentSwitchTimeIndex].systemIndex);
if (mCurrentSwitchTimeIndex == 1) {
mCurrentSwitchTimeIndex++; clearFault(1, 2, 3);
mCurrentSwitchTimeIndex++;
}
else {
switchSystemMatrix(mSwitchEventVector[mCurrentSwitchTimeIndex].systemIndex);
mElements = mElementsVector[mSwitchEventVector[mCurrentSwitchTimeIndex++].systemIndex];
}
//mCurrentSwitchTimeIndex++;
mLogger->Log(LogLevel::INFO) << "Switched to system " << mCurrentSwitchTimeIndex << " at " << mTime << std::endl; mLogger->Log(LogLevel::INFO) << "Switched to system " << mCurrentSwitchTimeIndex << " at " << mTime << std::endl;
mLogger->Log(LogLevel::INFO) << "New matrix:" << std::endl << mSystemModel.getCurrentSystemMatrix() << std::endl;
mLogger->Log(LogLevel::INFO) << "New decomp:" << std::endl << mSystemModel.getLUdecomp() << std::endl;
} }
} }
...@@ -226,52 +240,59 @@ Int Simulation::stepGeneratorTest(Logger& leftSideVectorLog, Logger& rightSideVe ...@@ -226,52 +240,59 @@ Int Simulation::stepGeneratorTest(Logger& leftSideVectorLog, Logger& rightSideVe
} }
} }
int Simulation::stepGeneratorVBR(Logger& leftSideVectorLog, Logger& rightSideVectorLog,
ElementPtr generator, Real time) {
// Set to zero because all components will add their contribution for the current time step to the current value
mSystemModel.getRightSideVector().setZero();
// Execute step for all circuit components void Simulation::clearFault(Real Node1, Real Node2, Real Node3) {
for (ElementList::iterator it = mElements.begin(); it != mElements.end(); ++it) {
(*it)->step(mSystemModel, mTime);
}
// Solve circuit for vector j with generator output current ClearingFault = true;
mSystemModel.solve();
// Execute PostStep for all components, generator states are recalculated based on new terminal voltage mIfa = getRightSideVector()(Node1 - 1);
for (ElementList::iterator it = mElements.begin(); it != mElements.end(); ++it) { mIfb = getRightSideVector()(Node2 - 1);
(*it)->postStep(mSystemModel); mIfc = getRightSideVector()(Node3 - 1);
}
if (mCurrentSwitchTimeIndex < mSwitchEventVector.size()) {
if (mTime >= mSwitchEventVector[mCurrentSwitchTimeIndex].switchTime) {
switchSystemMatrix(mSwitchEventVector[mCurrentSwitchTimeIndex].systemIndex);
mCurrentSwitchTimeIndex++; if (FirstTime == true)
mLogger->Log(LogLevel::INFO) << "Switched to system " << mCurrentSwitchTimeIndex << " at " << mTime << std::endl; {
} mIfa_hist = mIfa;
} mIfb_hist = mIfb;
mIfc_hist = mIfc;
// Save simulation step data FirstTime = false;
if (mLastLogTimeStep == 0) {
std::cout << mTime << std::endl;
leftSideVectorLog.LogDataLine(getTime(), getLeftSideVector());
rightSideVectorLog.LogDataLine(getTime(), getRightSideVector());
} }
mLastLogTimeStep++;
if (mLastLogTimeStep == mDownSampleRate) {
mLastLogTimeStep = 0;
}
if (mTime >= mFinalTime) {
return 0; if (signbit(mIfa) != signbit(mIfa_hist) && !aCleared) {
mElements.erase(mElements.begin() + 1);
addSystemTopology(mElements);
switchSystemMatrix(mSwitchEventVector.size() + NumClearedPhases);
NumClearedPhases++;
aCleared = true;
} }
else {
return 1; if (signbit(mIfb) != signbit(mIfb_hist) && !bCleared) {
mElements.erase(mElements.begin() + 2);
addSystemTopology(mElements);
switchSystemMatrix(mSwitchEventVector.size() + NumClearedPhases);
NumClearedPhases++;
bCleared = true;
}
if (signbit(mIfc) != signbit(mIfc_hist) && !cCleared) {
mElements.erase(mElements.begin() + 1);
addSystemTopology(mElements);
switchSystemMatrix(mSwitchEventVector.size() + NumClearedPhases);
NumClearedPhases++;
cCleared = true;
} }
mIfa_hist = mIfa;
mIfb_hist = mIfb;
mIfc_hist = mIfc;
if (NumClearedPhases == 3)
ClearingFault = false;
} }
void Simulation::switchSystemMatrix(Int systemMatrixIndex) { void Simulation::switchSystemMatrix(Int systemMatrixIndex) {