VoltageBehindReactanceEMT.cpp 11.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
/** Voltage behind reactance (EMT)
 *
 * @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/>.
 *********************************************************************************/

23
#include "VoltageBehindReactanceEMT.h"
Markus Mirz's avatar
Markus Mirz committed
24
#include "../IntegrationMethod.h"
25
26
27

using namespace DPsim;

28
VoltageBehindReactanceEMT::VoltageBehindReactanceEMT(String name, Int node1, Int node2, Int node3,
Markus Mirz's avatar
Markus Mirz committed
29
	Real nomPower, Real nomVolt, Real nomFreq, Int poleNumber, Real nomFieldCur,
30
	Real Rs, Real Ll, Real Lmd, Real Lmd0, Real Lmq, Real Lmq0,
31
32
	Real Rfd, Real Llfd, Real Rkd, Real Llkd,
	Real Rkq1, Real Llkq1, Real Rkq2, Real Llkq2,
33
	Real inertia, bool logActive)
34
35
36
	: SynchGenBase(name, node1, node2, node3, nomPower, nomVolt, nomFreq, poleNumber, nomFieldCur,
		Rs, Ll, Lmd, Lmd0, Lmq, Lmq0, Rfd, Llfd, Rkd, Llkd, Rkq1, Llkq1, Rkq2, Llkq2,
		inertia, logActive) {
Steffen Vogel's avatar
Steffen Vogel committed
37

38
}
39
40
41
42
43
44
45

VoltageBehindReactanceEMT::~VoltageBehindReactanceEMT() {
	if (mLogActive) {
		delete mLog;
	}
}

46

47

48
49
void VoltageBehindReactanceEMT::init(Real om, Real dt,
	Real initActivePower, Real initReactivePower, Real initTerminalVolt, Real initVoltAngle, Real initFieldVoltage, Real initMechPower) {
50

51
52
53
54
55
	mResistanceMat = Matrix::Zero(3, 3);
	mResistanceMat <<
		mRs, 0, 0,
		0, mRs, 0,
		0, 0, mRs;
56

57
58
59
60
	R_load <<
		1037.8378 / mBase_Z, 0, 0,
		0, 1037.8378 / mBase_Z, 0,
		0, 0, 1037.8378 / mBase_Z;
Viviane Sapucaia's avatar
Viviane Sapucaia committed
61

62
	//Dynamic mutual inductances
63

Viviane Sapucaia's avatar
Viviane Sapucaia committed
64
	mDLmd = 1. / (1. / mLmd + 1. / mLlfd + 1. / mLlkd);
65
	if (mNumDampingWindings == 2)
Viviane Sapucaia's avatar
Viviane Sapucaia committed
66
	{
Viviane Sapucaia's avatar
Viviane Sapucaia committed
67
		mDLmq = 1. / (1. / mLmq + 1. / mLlkq1 + 1. / mLlkq2);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
68
69
70
71
72
73
74
75
76
77
78
79

		A_flux <<
			-mRkq1 / mLlkq1*(1 - mDLmq / mLlkq1), mRkq1 / mLlkq1*(mDLmq / mLlkq2), 0, 0,
			mRkq2 / mLlkq2*(mDLmq / mLlkq1), -mRkq2 / mLlkq2*(1 - mDLmq / mLlkq2), 0, 0,
			0, 0, -mRfd / mLlfd*(1 - mDLmd / mLlfd), mRfd / mLlfd*(mDLmd / mLlkd),
			0, 0, mRkd / mLlkd*(mDLmd / mLlfd), -mRkd / mLlkd*(1 - mDLmd / mLlkd);
		B_flux <<
			mRkq1*mDLmq / mLlkq1, 0,
			mRkq2*mDLmq / mLlkq2, 0,
			0, mRfd*mDLmd / mLlfd,
			0, mRkd*mDLmd / mLlkd;
	}
80

Viviane Sapucaia's avatar
Viviane Sapucaia committed
81
	else
82
	{
Viviane Sapucaia's avatar
Viviane Sapucaia committed
83
		mDLmq = 1. / (1. / mLmq + 1. / mLlkq1);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
84

Markus Mirz's avatar
Markus Mirz committed
85
86
87
88
89
90
		A_flux = Matrix::Zero(3, 3);
		B_flux = Matrix::Zero(3, 2);
		C_flux = Matrix::Zero(3, 1);
		mRotorFlux = Matrix::Zero(3, 1);
		mDqStatorCurrents = Matrix::Zero(2, 1);
		mDqStatorCurrents_hist = Matrix::Zero(2, 1);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
91
92
93
94
95
96
97
98

		A_flux <<
			-mRkq1 / mLlkq1*(1 - mDLmq / mLlkq1), 0, 0,
			0, -mRfd / mLlfd*(1 - mDLmd / mLlfd), mRfd / mLlfd*(mDLmd / mLlkd),
			0, mRkd / mLlkd*(mDLmd / mLlfd), -mRkd / mLlkd*(1 - mDLmd / mLlkd);
		B_flux <<
			mRkq1*mDLmq / mLlkq1, 0,
			0, mRfd*mDLmd / mLlfd,
99
			0, mRkd*mDLmd / mLlkd;
Viviane Sapucaia's avatar
Viviane Sapucaia committed
100
101
	}

Viviane Sapucaia's avatar
Viviane Sapucaia committed
102
103
104
	mLa = (mDLmq + mDLmd) / 3.;
	mLb = (mDLmd - mDLmq) / 3.;

105
	// steady state per unit initial value
106
	initStatesInPerUnit(initActivePower, initReactivePower, initTerminalVolt, initVoltAngle, initFieldVoltage, initMechPower);
107

108
109
110
111
112
	// Correcting variables
	mThetaMech = mThetaMech + PI/2;
	mMechTorque = -mMechTorque;
	mIq = -mIq;
	mId = -mId;
113
114

	// #### VBR Model Dynamic variables #######################################
115
	
116
117
	mDPsid = mDLmd*(mPsifd / mLlfd) + mDLmd*(mPsikd / mLlkd);

118
	if (mNumDampingWindings == 2) {
119
		mDPsiq = mDLmq*(mPsikq1 / mLlkq1) + mDLmq*(mPsikq2 / mLlkq2);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
120
121
		mDVq = mOmMech*mDPsid + mDLmq*mRkq1*(mDPsiq - mPsikq1) / (mLlkq1*mLlkq1) +
			mDLmq*mRkq2*(mDPsiq - mPsikq2) / (mLlkq2*mLlkq2) + (mRkq1 / (mLlkq1*mLlkq1) + mRkq2 / (mLlkq2*mLlkq2))*mDLmq*mDLmq*mIq;
122
		mPsimq = mDLmq*(mPsikq1 / mLlkq1 + mPsikq2 / mLlkq2 + mIq);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
123
124
	}
	else {
125
		mDPsiq = mDLmq*(mPsikq1 / mLlkq1);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
126
		mDVq = mOmMech*mDPsid + mDLmq*mRkq1*(mDPsiq - mPsikq1) / (mLlkq1*mLlkq1) + (mRkq1 / (mLlkq1*mLlkq1))*mDLmq*mDLmq*mIq;
127
		mPsimq = mDLmq*(mPsikq1 / mLlkq1 + mIq);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
128
	}
129

130
131
	mDVd = -mOmMech*mDPsiq + mDLmd*mRkd*(mDPsid - mPsikd) / (mLlkd*mLlkd) + (mDLmd / mLlfd)*mVfd +
		mDLmd*mRfd*(mDPsid - mPsifd) / (mLlfd*mLlfd) + (mRfd / (mLlfd*mLlfd) + mRkd / (mLlkd*mLlkd))*mDLmd*mDLmd*mId;
132
	mPsimd = mDLmd*(mPsifd / mLlfd + mPsikd / mLlkd + mId);
133
134
135
136

	mDVa = inverseParkTransform(mThetaMech, mDVq, mDVd, 0)(0);
	mDVb = inverseParkTransform(mThetaMech, mDVq, mDVd, 0)(1);
	mDVc = inverseParkTransform(mThetaMech, mDVq, mDVd, 0)(2);
137
	
138
139
140
141
142
	mDVabc <<
		mDVa,
		mDVb,
		mDVc;
	mDVabc_hist = mDVabc;
143

144
145
146
	mVa = inverseParkTransform(mThetaMech, mVq, mVd, mV0)(0);
	mVb = inverseParkTransform(mThetaMech, mVq, mVd, mV0)(1);
	mVc = inverseParkTransform(mThetaMech, mVq, mVd, mV0)(2);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
147

148
149
150
	mIa = inverseParkTransform(mThetaMech, mIq, mId, mI0)(0);
	mIb = inverseParkTransform(mThetaMech, mIq, mId, mI0)(1);
	mIc = inverseParkTransform(mThetaMech, mIq, mId, mI0)(2);
151
152
}

153

154
void VoltageBehindReactanceEMT::step(SystemModel& system, Real time) {
155

156
157
158
	


159
160
	stepInPerUnit(system.getOmega(), system.getTimeStep(), time, system.getNumMethod());

161
162
	R_load = system.getCurrentSystemMatrix().inverse() / mBase_Z;

163
164
165
166
167
168
169
170
171
172
	// Update current source accordingly
	if (mNode1 >= 0) {
		system.addRealToRightSideVector(mNode1, -mIa*mBase_i);
	}
	if (mNode2 >= 0) {
		system.addRealToRightSideVector(mNode2, -mIb*mBase_i);
	}
	if (mNode3 >= 0) {
		system.addRealToRightSideVector(mNode3, -mIc*mBase_i);
	}
173

174
	if (mLogActive) {
175
176
		Matrix logValues(getRotorFluxes().rows() + getDqStatorCurrents().rows() + 3, 1);
		logValues << getRotorFluxes(), getDqStatorCurrents(), getElectricalTorque(), getRotationalSpeed(), getRotorPosition();
177
178
		mLog->LogDataLine(time, logValues);
	}
179

180

181
182
}

183

184
void VoltageBehindReactanceEMT::stepInPerUnit(Real om, Real dt, Real time, NumericalMethod numMethod) {
185
186
187
188
189
190

	mIabc <<
		mIa,
		mIb,
		mIc;

191
192
	// Calculate mechanical variables with euler
	mElecTorque = (mPsimd*mIq - mPsimq*mId);
193
	mOmMech = mOmMech + dt * (1. / (2. * mH) * (mElecTorque - mMechTorque));
Viviane Sapucaia's avatar
Viviane Sapucaia committed
194
	mThetaMech = mThetaMech + dt * (mOmMech* mBase_OmMech);
Steffen Vogel's avatar
Steffen Vogel committed
195

196

197
	// Calculate Inductance matrix and its derivative
198
	CalculateLandpL();
Viviane Sapucaia's avatar
Viviane Sapucaia committed
199

200
	// Solve circuit - calculate stator currents
201
	if (numMethod == NumericalMethod::Trapezoidal_flux) {
202
		mIabc = Trapezoidal(mIabc, -mDInductanceMat.inverse()*(mResistanceMat + R_load + pmDInductanceMat), mDInductanceMat.inverse(), dt*mBase_OmElec, -mDVabc, -mDVabc_hist);
203
204
205
		
	}
		
206
207
	else if (numMethod == NumericalMethod::Euler) 
		mIabc = Euler(mIabc, -mDInductanceMat.inverse()*(mResistanceMat + R_load + pmDInductanceMat), mDInductanceMat.inverse(), dt*mBase_OmElec, -mDVabc);
208
209
210
211
212
213
	


	mIa_hist = mIa;
	mIb_hist = mIb;
	mIc_hist = mIc;
214

Viviane Sapucaia's avatar
Viviane Sapucaia committed
215
216
217
	mIa = mIabc(0);
	mIb = mIabc(1);
	mIc = mIabc(2);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
218
219
	Real mIq_hist = mIq;
	Real mId_hist = mId;
Viviane Sapucaia's avatar
Viviane Sapucaia committed
220
221
	mIq = parkTransform(mThetaMech, mIa, mIb, mIc)(0);
	mId = parkTransform(mThetaMech, mIa, mIb, mIc)(1);
Steffen Vogel's avatar
Steffen Vogel committed
222
	mI0 = parkTransform(mThetaMech, mIa, mIb, mIc)(2);
223

Viviane Sapucaia's avatar
Viviane Sapucaia committed
224
	// Calculate rotor flux likanges
225
	if (mNumDampingWindings == 2)
Viviane Sapucaia's avatar
Viviane Sapucaia committed
226
227
228
229
230
231
232
	{
		C_flux <<
			0,
			0,
			mVfd,
			0;

233
		mDqStatorCurrents <<
Viviane Sapucaia's avatar
Viviane Sapucaia committed
234
235
			mIq,
			mId;
236
237

		mDqStatorCurrents_hist <<
Viviane Sapucaia's avatar
Viviane Sapucaia committed
238
239
240
			mIq_hist,
			mId_hist;

241
		mRotorFlux <<
Viviane Sapucaia's avatar
Viviane Sapucaia committed
242
243
244
245
246
			mPsikq1,
			mPsikq2,
			mPsifd,
			mPsikd;

247
248
249
250
		if (numMethod == NumericalMethod::Trapezoidal_flux)
			mRotorFlux = Trapezoidal(mRotorFlux, A_flux, B_flux, C_flux, dt*mBase_OmElec, mDqStatorCurrents, mDqStatorCurrents_hist);
		else if (numMethod == NumericalMethod::Euler)
			mRotorFlux = Euler(mRotorFlux, A_flux, B_flux, C_flux, dt*mBase_OmElec, mDqStatorCurrents);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
251

252
253
254
255
		mPsikq1 = mRotorFlux(0);
		mPsikq2 = mRotorFlux(1);
		mPsifd = mRotorFlux(2);
		mPsikd = mRotorFlux(3);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
256

Viviane Sapucaia's avatar
Viviane Sapucaia committed
257
	}
Viviane Sapucaia's avatar
Viviane Sapucaia committed
258
259
260
261
262
263
264
	else
	{
		C_flux <<
			0,
			mVfd,
			0;

265
		mDqStatorCurrents <<
Viviane Sapucaia's avatar
Viviane Sapucaia committed
266
267
			mIq,
			mId;
268
		mDqStatorCurrents_hist <<
Viviane Sapucaia's avatar
Viviane Sapucaia committed
269
270
271
			mIq_hist,
			mId_hist;

272
		mRotorFlux <<
Viviane Sapucaia's avatar
Viviane Sapucaia committed
273
274
275
276
			mPsikq1,
			mPsifd,
			mPsikd;

277
278
279
280
		if (numMethod == NumericalMethod::Trapezoidal_flux)
			mRotorFlux = Trapezoidal(mRotorFlux, A_flux, B_flux, C_flux, dt*mBase_OmElec, mDqStatorCurrents, mDqStatorCurrents_hist);
		else if (numMethod == NumericalMethod::Euler)
			mRotorFlux = Euler(mRotorFlux, A_flux, B_flux, C_flux, dt*mBase_OmElec, mDqStatorCurrents);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
281

282
283
284
		mPsikq1 = mRotorFlux(0);
		mPsifd = mRotorFlux(1);
		mPsikd = mRotorFlux(2);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
285
	}
286

Steffen Vogel's avatar
Steffen Vogel committed
287

288
	// Calculate dynamic flux likages
289
	if (mNumDampingWindings == 2) {
Viviane Sapucaia's avatar
Viviane Sapucaia committed
290
		mDPsiq = mDLmq*(mPsikq1 / mLlkq1) + mDLmq*(mPsikq2 / mLlkq2);
291
		mPsimq = mDLmq*(mPsikq1 / mLlkq1 + mPsikq2 / mLlkq2 + mIq);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
292
293
294
	}
	else {
		mDPsiq = mDLmq*(mPsikq1 / mLlkq1);
295
		mPsimq = mDLmq*(mPsikq1 / mLlkq1 + mIq);
Viviane Sapucaia's avatar
Viviane Sapucaia committed
296
	}
297
	mDPsid = mDLmd*(mPsifd / mLlfd) + mDLmd*(mPsikd / mLlkd);
298
	mPsimd = mDLmd*(mPsifd / mLlfd + mPsikd / mLlkd + mId);
Steffen Vogel's avatar
Steffen Vogel committed
299

300
301

	// Calculate dynamic voltages
302
	if (mNumDampingWindings == 2) {
Viviane Sapucaia's avatar
Viviane Sapucaia committed
303
304
305
306
307
308
		mDVq = mOmMech*mDPsid + mDLmq*mRkq1*(mDPsiq - mPsikq1) / (mLlkq1*mLlkq1) +
			mDLmq*mRkq2*(mDPsiq - mPsikq2) / (mLlkq2*mLlkq2) + (mRkq1 / (mLlkq1*mLlkq1) + mRkq2 / (mLlkq2*mLlkq2))*mDLmq*mDLmq*mIq;
	}
	else {
		mDVq = mOmMech*mDPsid + mDLmq*mRkq1*(mDPsiq - mPsikq1) / (mLlkq1*mLlkq1) + (mRkq1 / (mLlkq1*mLlkq1))*mDLmq*mDLmq*mIq;
	}
309
310
	mDVd = -mOmMech*mDPsiq + mDLmd*mRkd*(mDPsid - mPsikd) / (mLlkd*mLlkd) + (mDLmd / mLlfd)*mVfd +
		mDLmd*mRfd*(mDPsid - mPsifd) / (mLlfd*mLlfd) + (mRfd / (mLlfd*mLlfd) + mRkd / (mLlkd*mLlkd))*mDLmd*mDLmd*mId;
311
312
313
314

	mDVa = inverseParkTransform(mThetaMech, mDVq, mDVd, 0)(0);
	mDVb = inverseParkTransform(mThetaMech, mDVq, mDVd, 0)(1);
	mDVc = inverseParkTransform(mThetaMech, mDVq, mDVd, 0)(2);
315
316
317
318
319
320
	mDVabc_hist = mDVabc;
	mDVabc <<
		mDVa,
		mDVb,
		mDVc;

321
322
323



324

Viviane Sapucaia's avatar
Viviane Sapucaia committed
325
326
327
}


328
void VoltageBehindReactanceEMT::postStep(SystemModel& system) {
Steffen Vogel's avatar
Steffen Vogel committed
329

330
}
331

332
333
334
335
336
337
338
339
340
341
void VoltageBehindReactanceEMT::CalculateLandpL() {
	mDInductanceMat <<
		mLl + mLa - mLb*cos(2 * mThetaMech), -mLa / 2 - mLb*cos(2 * mThetaMech - 2 * PI / 3), -mLa / 2 - mLb*cos(2 * mThetaMech + 2 * PI / 3),
		-mLa / 2 - mLb*cos(2 * mThetaMech - 2 * PI / 3), mLl + mLa - mLb*cos(2 * mThetaMech - 4 * PI / 3), -mLa / 2 - mLb*cos(2 * mThetaMech),
		-mLa / 2 - mLb*cos(2 * mThetaMech + 2 * PI / 3), -mLa / 2 - mLb*cos(2 * mThetaMech), mLl + mLa - mLb*cos(2 * mThetaMech + 4 * PI / 3);
	pmDInductanceMat <<
		mLb*sin(2 * mThetaMech), mLb*sin(2 * mThetaMech - 2 * PI / 3), mLb*sin(2 * mThetaMech + 2 * PI / 3),
		mLb*sin(2 * mThetaMech - 2 * PI / 3), mLb*sin(2 * mThetaMech - 4 * PI / 3), mLb*sin(2 * mThetaMech),
		mLb*sin(2 * mThetaMech + 2 * PI / 3), mLb*sin(2 * mThetaMech), mLb*sin(2 * mThetaMech + 4 * PI / 3);
	pmDInductanceMat = pmDInductanceMat * 2 * mOmMech;
342
343
344
}


Markus Mirz's avatar
Markus Mirz committed
345
Matrix VoltageBehindReactanceEMT::parkTransform(Real theta, Real a, Real b, Real c) {
346

Markus Mirz's avatar
Markus Mirz committed
347
	Matrix dq0vector(3, 1);
348

Markus Mirz's avatar
Markus Mirz committed
349
	Real q, d;
350

Steffen Vogel's avatar
Steffen Vogel committed
351
	q = 2. / 3. * cos(theta)*a + 2. / 3. * cos(theta - 2. * M_PI / 3.)*b + 2. / 3. * cos(theta + 2. * M_PI / 3.)*c;
352
	d = 2. / 3. * sin(theta)*a + 2. / 3. * sin(theta - 2. * M_PI / 3.)*b + 2. / 3. * sin(theta + 2. * M_PI / 3.)*c;
353

354
355
356
357
358
359
360
361
362
	dq0vector << q,
		d,
		0;

	return dq0vector;
}



Markus Mirz's avatar
Markus Mirz committed
363
Matrix VoltageBehindReactanceEMT::inverseParkTransform(Real theta, Real q, Real d, Real zero) {
364

Markus Mirz's avatar
Markus Mirz committed
365
	Matrix abcVector(3, 1);
366

Markus Mirz's avatar
Markus Mirz committed
367
	Real a, b, c;
368

Viviane Sapucaia's avatar
Viviane Sapucaia committed
369
370
371
	a = cos(theta)*q + sin(theta)*d + 1.*zero;
	b = cos(theta - 2. * M_PI / 3.)*q + sin(theta - 2. * M_PI / 3.)*d + 1.*zero;
	c = cos(theta + 2. * M_PI / 3.)*q + sin(theta + 2. * M_PI / 3.)*d + 1.*zero;
372
373
374
375
376
377

	abcVector << a,
		b,
		c;

	return abcVector;
Viviane Sapucaia's avatar
Viviane Sapucaia committed
378
379
}