Commit 4b7a5400 authored by Georg Martin Reinke's avatar Georg Martin Reinke
Browse files

Merge remote-tracking branch 'origin/ci' into shmem-interface

parents 01faf730 2b52e62c
stages:
- prepare
- build
docker:
stage: prepare
script:
- docker build -t dpsim-dev .
tags:
- shell
build:
stage: build
script:
- make -C Source
image: dpsim-dev
tags:
- docker
\ No newline at end of file
FROM fedora:latest
MAINTAINER Steffen Vogel <stvogel@eonerc.rwth-aachen.de>
RUN dnf -y install \
gcc-c++ \
make \
eigen3-devel \
doxygen
\ No newline at end of file
%----------------------------------------------------------
%----------------------------------------------------------
% Packages and other configurations
%----------------------------------------------------------
......
......@@ -377,6 +377,188 @@ In both cases, the function "init" is setting every initial value to zero.
The model is according to \cite{wang2010methods} and \cite{kundur1994power}.
\subsubsection{Prerequisites}
\subsubsection{Magnetic Circuits}
Considering two magnetically coupled windings with number of turns $N_1$ and $N_2$, resistance $r_1$ and $r_2$ and winding currents $i_1$ and $i_2$ respectively, the flux linkage with the respective windings, produced by the total effect of both currents are given by equations \ref{eq:FluxLinkage1} and \ref{eq:FluxLinkage2}.
\begin{equation} \label{eq:FluxLinkage1}
\Psi_1 = N_1(\Phi_{l1} + \Phi_{m1} + \Phi_{m2})
\end{equation}
\begin{equation} \label{eq:FluxLinkage2}
\Psi_2 = N_2(\Phi_{l2} + \Phi_{m2} + \Phi_{m1})
\end{equation}
\\
where
$\Phi_{m1}$ is the mutual flux linking both windings due to current $i_1$;
$\Phi_{m2}$ is the mutual flux linking both windings due to current $i_2$;
$\Phi_{l1}$ is the leakage flux linking winding 1 only;
$\Phi_{l2}$ is the leakage flux linking winding 2 only.
The flux linkage can also be expressed in terms of self and mutual inductance as in equation \ref{eq:FluxLinkageInductances}
\begin{equation} \label{eq:FluxLinkageInductances}
\begin{bmatrix}
\Psi_1\\
\Psi_2
\end{bmatrix}
=
\begin{bmatrix}
L_{11} & L_{12}\\
L_{21} & L_{22}
\end{bmatrix}
\cdot
\begin{bmatrix}
i_1\\
i_2
\end{bmatrix}
\end{equation}
\\
where $L_{11}$ and $L_{22}$ are the self inductances of windings 1 and 2 respectively and defined as the flux linkage per unit current in the same winding. $L_{12}$ and $L_{21}$ are the mutual inductance between two windings, defined as the flux linkage with one winding per unit current in the other winding.
\begin{equation}
L_{11} = \frac{N_1(\Phi_{m1} + \Phi_{l1})}{i_1}
\end{equation}
\begin{equation}
L_{22} = \frac{N_2(\Phi_{m2} + \Phi_{l2})}{i_2}
\end{equation}
\begin{equation}
L_{12} = \frac{N_1 \Phi_{m2}}{i_2}
\end{equation}
\begin{equation}
L_{21} = \frac{N_2 \Phi_{m1}}{i_1}
\end{equation}
The mutual fluxes $\Phi_{m1}$ and $\Phi_{m2}$ can be defined as in equations \ref{eq:MutualFlux1} and \ref{eq:MutualFlux2}, where $\Re_m$ is the reluctance of the path of the magnetizing fluxes.
\begin{equation} \label{eq:MutualFlux1}
\Phi_{m1} = \frac{N_1 i_1}{\Re_m}
\end{equation}
\begin{equation} \label{eq:MutualFlux2}
\Phi_{m2} = \frac{N_2 i_2}{\Re_m}
\end{equation}
From these equations, it is possible to see that
\begin{equation}
L_{12}=L_{21}= \frac{N_1 N_2}{\Re_m}
\end{equation}
The terminal voltages $e_1$ and $e_2$ are defined as following:
\begin{equation} \label{eq:VoltageEquation1}
e_1 = \frac{d \Psi_1}{d t} + r_1 i_1
\end{equation}
\begin{equation} \label{eq:VoltageEquation2}
e_2 = \frac{d \Psi_2}{d t} + r_2 i_2
\end{equation}
Equations \ref{eq:FluxLinkageInductances}, \ref{eq:VoltageEquation1} and \ref{eq:VoltageEquation2} give the performance of the magnetically coupled circuit.
\subsubsection{Synchronous Machine Basic equations}
We can derive the equations of the synchronous machine similarly to the equations of the magnetic circuit.
The flux linkage for stator phase windings a, b and c are given by
\begin{equation}
\begin{bmatrix}
\Psi_a \\
\Psi_b \\
\Psi_c
\end{bmatrix}
=
\begin{bmatrix}
-l_{aa} & -l_{ab} & -l_{ac} & l_{afd} & l_{akd} & l_{akq} \\
-l_{ab} & -l_{bb} & -l_{bc} & l_{afd} & l_{akd} & l_{akq} \\
-l_{ac} & -l_{bc} & -l_{cc} & l_{afd} & l_{akd} & l_{akq}
\end{bmatrix}
\begin{bmatrix}
i_a \\
i_b \\
i_c \\
i_{fd} \\
i_{kd} \\
i_{kq}
\end{bmatrix}
\end{equation}
where
$l_{aa}$, $l_{bb}$ and $l_{cc}$ are the self-inductances of the stator windings;
$l_{ab}$, $l_{bc}$ and $l_{ac}$ are the mutual inductances between the stator windings;
$l_{afd}$, $l_{akd}$ and $l_{akq}$ are the mutual inductances between stator and rotor windings;
$i_a$, $i_b$ and $i_c$ are the instantaneous stator currents in phases a, b and c;
$i_{fd}$, $i_{kd}$ and $i_{kd}$ are the field and amortisseur circuit currents.
The equations of the stator phase to neutral voltages are
\begin{equation}
\begin{bmatrix}
e_a \\
e_b \\
e_c
\end{bmatrix}
=
\frac{d}{dt}
\begin{bmatrix}
\Psi_a \\
\Psi_b \\
\Psi_c
\end{bmatrix}
-R_a
\begin{bmatrix}
i_a \\
i_b \\
i_c
\end{bmatrix}
\end{equation}
where $R_a$ is the armature resistance per phase.
It is important to notice that the inductances are time-varying. They vary with the rotor position, due to the variation in the reluctance of the magnetic flux path (non-uniform air gap).
The permeance of the magnetic flux path $P$ (inverse of the reluctance $\Re$) in function of the rotor position $\alpha$ can be writen as
\begin{equation}
P = P_0 + P_2 \cos 2 \alpha
\end{equation}
Inductance is directly proportional to permeance and the relative position between the phase a and the rotor is $\theta$. Therefore, we can write the self-inductance $l_{aa}$ as in equation \ref{eq:SelfInductanceA}. Knowing that the windings of phases b and c are identical to phase a and displaced from it by $120 \circ$ and $240 \circ$ respectively, we can also determine the self-inductances $l_{bb}$ and $l_{cc}$. A more detailed derivation of the self-inductances is in \cite{kundur1994power}.
\begin{equation} \label{eq:SelfInductanceA}
l_{aa} = L_{aa0} + L_{aa2} \cos 2 \theta
\end{equation}
\begin{equation} \label{eq:SelfInductanceA}
l_{bb} = L_{aa0} + L_{aa2} \cos 2 (\theta - \frac{2 \pi}{3})
\end{equation}
\begin{equation} \label{eq:SelfInductanceA}
l_{cc} = L_{aa0} + L_{aa2} \cos 2 (\theta + \frac{2 \pi}{3})
\end{equation}
The mutual inductance between two windings can be deduced evaluating the flux linking one phase due to the current in the other phase. The mutual inductance will also have a second harmonic variation due to the non-uniform air gap.
\begin{equation} \label{eq:SelfInductanceA}
l_{ab} = -L_{ab0} - L_{ab2} \cos (2 \theta + \frac{\pi}{3})
\end{equation}
\begin{equation} \label{eq:SelfInductanceA}
l_{bc} = -L_{ab0} - L_{ab2} \cos (2 \theta - \pi)
\end{equation}
\begin{equation} \label{eq:SelfInductanceA}
l_{ac} = - L_{ab0} - L_{ab2} \cos (2 \theta - \frac{\pi}{3})
\end{equation}
\subsubsection{Park's transformation}
Park's transformation is commonly used to achieve a model with static parameters:
%
\begin{equation}
......
% circuit with R, L and C
fileName1='../VisualStudio/DPsimVS2015/Logs/LeftVectorLog_SimulationExamplePiLine_0.001.csv'
%model with 3 nodes
fileName2='../VisualStudio/DPsimVS2015/Logs/LeftVectorLog_SimulationExamplePiLine2_0.001.csv';
compareDpResults(fileName1,fileName2,1,1,'Model with 3 nodes x Resistor + Inductor');
......@@ -18,5 +18,6 @@
#include "Components/VoltSourceResEMT.h"
#include "Components/IdealVoltageSource.h"
#include "Components/RxLine.h"
#include "Components/PiLine.h"
#endif // !COMPONENTS_H
\ No newline at end of file
#include "PiLine.h"
using namespace DPsim;
PiLine::PiLine(std::string name, int node1, int node2, int node3, Real resistance, Real inductance, Real capacitance) : BaseComponent(name, node1, node2, node3) {
mResistance = resistance;
mConductance = 1.0 / resistance;
mInductance = inductance;
mCapacitance = capacitance / 2;
}
void PiLine::applySystemMatrixStamp(SystemModel& system) {
Real a = system.getTimeStep() / (2. * mInductance);
Real b = system.getTimeStep() * system.getOmega() / 2.;
mGlr = a / (1 + b*b);
mGli = -a*b / (1 + b*b);
mPrevCurFacRe = (1 - b*b) / (1 + b*b);
mPrevCurFacIm = -2. * b / (1 + b*b);
// Resistive part
// Set diagonal entries
if (mNode1 >= 0) {
system.addCompToSystemMatrix(mNode1, mNode1, mConductance, 0);
}
if (mNode3 >= 0) {
system.addCompToSystemMatrix(mNode3, mNode3, mConductance, 0);
}
// Set off diagonal entries
if (mNode1 >= 0 && mNode3 >= 0) {
system.addCompToSystemMatrix(mNode1, mNode3, -mConductance, 0);
system.addCompToSystemMatrix(mNode3, mNode1, -mConductance, 0);
}
// Inductance part
// Set diagonal entries
if (mNode3 >= 0) {
system.addCompToSystemMatrix(mNode3, mNode3, mGlr, mGli);
}
if (mNode2 >= 0) {
system.addCompToSystemMatrix(mNode2, mNode2, mGlr, mGli);
}
if (mNode3 >= 0 && mNode2 >= 0) {
system.addCompToSystemMatrix(mNode3, mNode2, -mGlr, -mGli);
system.addCompToSystemMatrix(mNode2, mNode3, -mGlr, -mGli);
}
//capacitive part
mGcr = 2.0 * mCapacitance / system.getTimeStep();
mGci = system.getOmega() * mCapacitance;
if (mNode1 >= 0) {
system.addCompToSystemMatrix(mNode1, mNode1, mGcr, mGci);
}
if (mNode2 >= 0) {
system.addCompToSystemMatrix(mNode2, mNode2, mGcr, mGci);
}
}
void PiLine::init(Real om, Real dt) {
// Initialize internal state
mCurrIndRe = 0;
mCurrIndIm = 0;
mCurrCapRe1 = 0;
mCurrCapIm1 = 0;
mCurrCapRe2 = 0;
mCurrCapIm2 = 0;
mCurEqIndRe = 0;
mCurEqIndIm = 0;
mCurEqCapRe1 = 0;
mCurEqCapIm1 = 0;
mCurEqCapRe2 = 0;
mCurEqCapIm2 = 0;
mDeltaVre = 0;
mDeltaVim = 0;
mVoltageAtNode1Re = 0;
mVoltageAtNode1Im = 0;
mVoltageAtNode2Re = 0;
mVoltageAtNode2Im = 0;
}
void PiLine::step(SystemModel& system, Real time) {
// Initialize internal state inductance
mCurEqIndRe = mGlr * mDeltaVre - mGli * mDeltaVim + mPrevCurFacRe * mCurrIndRe - mPrevCurFacIm * mCurrIndIm;
mCurEqIndIm = mGli * mDeltaVre + mGlr * mDeltaVim + mPrevCurFacIm * mCurrIndRe + mPrevCurFacRe * mCurrIndIm;
if (mNode3 >= 0) {
system.addCompToRightSideVector(mNode3, -mCurEqIndRe, -mCurEqIndIm);
}
if (mNode2 >= 0) {
system.addCompToRightSideVector(mNode2, mCurEqIndRe, mCurEqIndIm);
}
// Initialize internal state capacitance 1
mCurEqCapRe1 = mCurrCapRe1 + mGcr * mVoltageAtNode1Re + mGci * mVoltageAtNode1Im;
mCurEqCapIm1 = mCurrCapIm1 + mGcr * mVoltageAtNode1Im - mGci * mVoltageAtNode1Re;
if (mNode1 >= 0) {
system.addCompToRightSideVector(mNode1, mCurEqCapRe1, mCurEqCapIm1);
}
// Initialize internal state capacitance 2
mCurEqCapRe2 = mCurrCapRe2 + mGcr * mVoltageAtNode2Re + mGci * mVoltageAtNode2Im;
mCurEqCapIm2 = mCurrCapIm2 + mGcr * mVoltageAtNode2Im - mGci * mVoltageAtNode2Re;
if (mNode2 >= 0) {
system.addCompToRightSideVector(mNode2, mCurEqCapRe2, mCurEqCapIm2);
}
}
void PiLine::postStep(SystemModel& system) {
Real vposr, vnegr, vposi, vnegi;
// extract solution
if (mNode3 >= 0) {
system.getRealFromLeftSideVector(mNode3);
vposr = system.getRealFromLeftSideVector(mNode3);
vposi = system.getImagFromLeftSideVector(mNode3);
}
if (mNode2 >= 0) {
system.getRealFromLeftSideVector(mNode2);
vnegr = system.getRealFromLeftSideVector(mNode2);
vnegi = system.getImagFromLeftSideVector(mNode2);
}
mDeltaVre = vposr - vnegr;
mDeltaVim = vposi - vnegi;
mCurrIndRe = mGlr * mDeltaVre - mGli * mDeltaVim + mCurEqIndRe;
mCurrIndIm = mGli * mDeltaVre + mGlr * mDeltaVim + mCurEqIndIm;
// extract solution
if (mNode1 >= 0) {
mVoltageAtNode1Re = system.getRealFromLeftSideVector(mNode1);
mVoltageAtNode1Im = system.getImagFromLeftSideVector(mNode1);
}
if (mNode2 >= 0) {
mVoltageAtNode2Re = system.getRealFromLeftSideVector(mNode2);
mVoltageAtNode2Im = system.getImagFromLeftSideVector(mNode2);
}
mCurrCapRe1 = mGcr * mVoltageAtNode1Re - mGci * mVoltageAtNode1Im - mCurEqCapRe1;
mCurrCapIm1 = mGci * mVoltageAtNode1Re + mGcr * mVoltageAtNode1Im - mCurEqCapIm1;
mCurrCapRe2 = mGcr * mVoltageAtNode2Re - mGci * mVoltageAtNode2Im - mCurEqCapRe2;
mCurrCapIm2 = mGci * mVoltageAtNode2Re + mGcr * mVoltageAtNode2Im - mCurEqCapIm2;
}
#pragma once
#ifndef PILINE_H
#define PILINE_H
#include "../Components.h"
namespace DPsim {
class PiLine : public BaseComponent {
protected:
///resistance in ohms
Real mResistance;
///Conductance
Real mConductance;
///Capacitance
Real mCapacitance;
///Inductance
Real mInductance;
Real mVoltageAtNode1Re;
Real mVoltageAtNode1Im;
Real mVoltageAtNode2Re;
Real mVoltageAtNode2Im;
Real mDeltaVre;
Real mDeltaVim;
Real mCurrIndRe;
Real mCurrIndIm;
Real mCurrCapRe1;
Real mCurrCapIm1;
Real mCurrCapRe2;
Real mCurrCapIm2;
Real mCurEqIndRe;
Real mCurEqIndIm;
Real mCurEqCapRe1;
Real mCurEqCapIm1;
Real mCurEqCapRe2;
Real mCurEqCapIm2;
Real mGlr;
Real mGli;
Real mGcr;
Real mGci;
Real mPrevCurFacRe;
Real mPrevCurFacIm;
public:
PiLine() { };
//PiLine(std::string nam3e, int node1, int node2, Real resistance, Real inductance);
PiLine(std::string name, int node1, int node2, int node3, Real resistance, Real inductance, Real capacitance);
void init(Real om, Real dt);
void applySystemMatrixStamp(SystemModel& system);
void applyRightSideVectorStamp(SystemModel& system) { }
void step(SystemModel& system, Real time);
void postStep(SystemModel& system);
};
}
#endif
\ No newline at end of file
......@@ -11,8 +11,6 @@ using namespace DPsim;
int main(int argc, char* argv[]) {
//villasExample();
villasDistributedExample(argc, argv);
//villasDistributedRef();
......
......@@ -395,4 +395,87 @@ void DPsim::simulationExampleRXLine3()
log.WriteLogToFile("Logs/Log_" + fileName.str() + ".log");
leftVectorLog.WriteLogToFile("Logs/LeftVectorLog_" + fileName.str() + ".csv");
rightVectorLog.WriteLogToFile("Logs/RightVectorLog_" + fileName.str() + ".csv");
}
void DPsim::simulationExamplePiLine()
{
// Define Object for saving data on a file
Logger log, leftVectorLog, rightVectorLog;
std::vector<BaseComponent*> circElements0;
circElements0.push_back(new IdealVoltageSource("v_1", 1, 0, 345, 0, 1));
circElements0.push_back(new LinearResistor("r1", 1, 2, 5));
circElements0.push_back(new PiLine("PiLine1", 2, 3, 4, 6.4, 0.186, 0.004));
circElements0.push_back(new LinearResistor("r_load", 3, 0, 150));
std::cout << "The contents of circElements0 are:";
for (std::vector<BaseComponent*>::iterator it = circElements0.begin(); it != circElements0.end(); ++it) {
std::cout << "Added " << (*it)->getName() << std::endl;
}
std::cout << '\n';
// Set up simulation
Real timeStep = 0.001;
Simulation newSim(circElements0, 2.0*M_PI*50.0, timeStep, 0.3, log);
// Main Simulation Loop
std::cout << "Start simulation." << std::endl;
while (newSim.step(log, leftVectorLog, rightVectorLog))
{
newSim.increaseByTimeStep();
updateProgressBar(newSim.getTime(), newSim.getFinalTime());
}
std::cout << "Simulation finished." << std::endl;
// Write simulation data to file
std::ostringstream fileName;
fileName << "SimulationExamplePiLine_" << timeStep;
log.WriteLogToFile("Logs/Log_" + fileName.str() + ".log");
leftVectorLog.WriteLogToFile("Logs/LeftVectorLog_" + fileName.str() + ".csv");
rightVectorLog.WriteLogToFile("Logs/RightVectorLog_" + fileName.str() + ".csv");
}
void DPsim::simulationExamplePiLine2()
{
// Define Object for saving data on a file
Logger log, leftVectorLog, rightVectorLog;
std::vector<BaseComponent*> circElements0;
circElements0.push_back(new IdealVoltageSource("v_1", 1, 0, 345, 0, 1));
circElements0.push_back(new LinearResistor("r1", 1, 2, 5));
circElements0.push_back(new Capacitor("c_1", 2, 0, 0.002));
circElements0.push_back(new LinearResistor("r_load", 2, 4, 6.4));
circElements0.push_back(new Inductor("l_1", 4, 3, 0.186));
circElements0.push_back(new Capacitor("c_2", 3, 0, 0.002));
circElements0.push_back(new LinearResistor("r_load", 3, 0, 150));
std::cout << "The contents of circElements0 are:";
for (std::vector<BaseComponent*>::iterator it = circElements0.begin(); it != circElements0.end(); ++it) {
std::cout << "Added " << (*it)->getName() << std::endl;
}
std::cout << '\n';
// Set up simulation
Real timeStep = 0.001;
Simulation newSim(circElements0, 2.0*M_PI*50.0, timeStep, 0.3, log);
// Main Simulation Loop
std::cout << "Start simulation." << std::endl;
while (newSim.step(log, leftVectorLog, rightVectorLog))
{
newSim.increaseByTimeStep();
updateProgressBar(newSim.getTime(), newSim.getFinalTime());
}
std::cout << "Simulation finished." << std::endl;
// Write simulation data to file
std::ostringstream fileName;
fileName << "SimulationExamplePiLine2_" << timeStep;
log.WriteLogToFile("Logs/Log_" + fileName.str() + ".log");
leftVectorLog.WriteLogToFile("Logs/LeftVectorLog_" + fileName.str() + ".csv");
rightVectorLog.WriteLogToFile("Logs/RightVectorLog_" + fileName.str() + ".csv");
}
\ No newline at end of file
......@@ -15,5 +15,7 @@ namespace DPsim {
void simulationExampleRXLine();
void simulationExampleRXLine2();
void simulationExampleRXLine3();
void simulationExamplePiLine();
void simulationExamplePiLine2();
}
#endif
......@@ -12,7 +12,7 @@ TEST_SOURCES = $(wildcard Tests/Test*.cpp)
TEST_OBJS = $(TEST_SOURCES:.cpp=.o)
TEST_BINS = $(TEST_SOURCES:.cpp=)
INCLUDES = -I /usr/include/eigen3
INCLUDES = -I /usr/local/include/eigen -I /usr/include/eigen3
.PHONY: tests clean
......
......@@ -47,7 +47,7 @@ Simulation::~Simulation() {
void Simulation::initialize(std::vector<BaseComponent*> newElements) {
int maxNode = 0;
Int numIdealVS = 0;
int numRxLines = 0;
int numLines = 0;