暂态.cpp

电力系统分析的NR法潮流计算、暂态稳定计算VC++源码。

#include<iostream.h>
#include<stdio.h>
#include<stdlib.h>
#include "stdlib.h"
#include<iomanip.h>
#include<math.h>
#include<fstream.h>
double pi = 3.14159;



void main(){
int n,G,kind;
double dt,Vg,SB,UB,F1,Tj,Xd,X2,k1,k2,Us1,Us2,XL1,U,P0,F2,xg2,q0,X1,X3,T1,T2;
double jd[50],djd[50],dw[50],jd0[50],w0[50],djd0[50],dw0[50],djda[50],dwa[50],w[50],T[50];

G = 240;
Vg = 10.5;
X2 = 0.44;
Xd = 0.3;
T1 = 300;
k1 = 0.04339;
Us1 = 14;
T2 = 280;
k2 = 1.8182;
Us2 = 14;
XL1 = 0.42;
double Line = 230;
U = 115;
P0 = 220;
F2 = 0.98;
SB = 220;
UB = 209;
Tj = 6;
F1 = 0.8;
double Vt1 = 242;
double Vt2 = 220;
//
double ub1 = UB * k1;
double ub2 = UB / k2;
double xg1 = Xd * Vg * Vg * SB * F1 / G / ub1 / ub1;
double Xt11 = Us1 * Vt1 * Vt1 * SB / T1 / UB / UB / 100;
double Xl11 = XL1 * Line * SB / UB / UB / 2;
double Xt21 = Us2 * Vt2 * Vt2 * SB / T2 / UB / UB / 100;
xg2 = X2 * Vg * Vg * SB * F1 / G / ub1 / ub1;
Tj = Tj * G / SB / F1;
X1 = xg1 + Xt11 + Xl11 + Xt21;
double s0 = P0 / F2;
q0 = sqrt(s0 * s0 - P0 * P0);
q0 = q0 / SB;
double s = (1 + q0 * X1) *(1 + q0 * X1) + (P0 * X1 / SB)*(P0 * X1 / SB);
double e = sqrt(s);
double angle0 = atan((P0 * X1 / SB) / (1 + q0 * X1)) * 180 / pi;

//***********考虑线路始端短路
kind =4;
double X12,X10,X,p2m,p3m,pe;
if (kind==1){
	X12 = (xg2 + Xt11) * (Xl11 + Xt21) / ((xg2 + Xt11) + (Xl11 + Xt21));
	X10 = Xt11 * (4 * Xl11 + Xt21) / (Xt11 + 4 * Xl11 + Xt21);
	X = X12 + X10;
	X2 = (xg1 + Xt11) + (Xl11 + Xt21) + (xg1 + Xt11) * (Xl11 + Xt21) / X;
	p2m = e / X2;
	X3 = xg1 + Xt11 + 2 * Xl11 + Xt21;
	p3m = e / X3;
	}else if (kind==2){
	X12 = (xg2 + Xt11) * (Xl11 + Xt21) / ((xg2 + Xt11) + (Xl11 + Xt21));
	X = X12;
	X2 = (xg1 + Xt11) + (Xl11 + Xt21) + (xg1 + Xt11) * (Xl11 + Xt21) / X;
	p2m = e / X2;
	X3 = xg1 + Xt11 + 2 * Xl11 + Xt21;
	p3m = e / X3;
	}else if (kind==3){
	p2m = 0;
	X3 = xg1 + Xt11 + 2 * Xl11 + Xt21;
	p3m = e / X3;
	}else{
	X12 = (xg2 + Xt11) * (Xl11 + Xt21) / ((xg2 + Xt11) + (Xl11 + Xt21));
	X10 = Xt11 * (4 * Xl11 + Xt21) / (Xt11 + 4 * Xl11 + Xt21);
	X = X12 * X10 / (X12 + X10);
	X2 = (xg1 + Xt11) + (Xl11 + Xt21) + (xg1 + Xt11) * (Xl11 + Xt21) / X;
	p2m = e / X2;
	X3 = xg1 + Xt11 + 2 * Xl11 + Xt21;
	p3m = e / X3;
	}
double angleh = 180 - atan(1 / (p3m * p3m - 1)) * 180 / pi;
s = ((angleh - angle0) * pi / 180 + p3m * cos(angleh * pi / 180) - p2m * cos(angle0 * pi / 180)) / (p3m - p2m);
double anglec = atan(sqrt(1-s*s ) / s) * 180 / pi;
if (anglec < 0 ){
	anglec = anglec + 180;
}



T1 = 0.15 ;// 故障切除时间
T2 = 1.0;//计 算 时 间
dt = 0.05;
w[0] = 1;
jd[0] = angle0;
double pt = 1;
double f0 = 50;
n = 0;

for (n=0;n<T2/dt; n++){

	if (n > (T1 / dt) ){
		pe = p3m;
	}else{
		pe = p2m;
	}
	djd[n] = (w[n] - 1) * 360 * f0;
	dw[n] = (pt - pe * sin(jd[n] * pi / 180)) / Tj;
	jd0[n + 1] = jd[n] + djd[n] * dt;
	w0[n + 1] = w[n] + dw[n] * dt;
	djd0[n + 1] = (w0[n + 1] - 1) * 360 * f0;
	dw0[n + 1] = (pt - pe * sin(jd0[n + 1] * pi / 180)) / Tj;
	djda[n + 1]= (djd[n] + djd0[n + 1]) / 2;
	dwa[n + 1] = (dw[n] + dw0[n + 1]) / 2;
	jd[n + 1] = jd[n] + djda[n + 1] * dt;
	w[n + 1] = w[n] + dwa[n + 1] * dt;
	T[n] = (n + 1) * dt;
	//输出
	cout<<"时间："<<T[n]<<"--角度"<<jd[n + 1]<<endl;
}
}