回复 2# 的帖子
#define PI 3.1415926
#define G 9.81
#define E_I (9.8*pow(10,10))
#define E_S (20.6*pow(10,10))
#define E_W (2.06*pow(10,9))
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
void main()
{int i,j;
double ws[50];
/*各段波速*/
double t[50];
double tsum;
double dt;
double dx[50];
double err=0.1;
int k;
/*时间层*/
double T;
/*时间*/
double Tm=12.0;
/*过渡过程计算时间*/
int n1=2;
/*并联泵的个数*/
int n2=2;
/*管段根据几何特征分的段数*/
int num=60;
/*管道特征线计算欲分段*/
double flow0=0.593;
/*单泵额定流量*/
double head0=39.4;
/*单泵额定扬程*/
double n0=970.0;
/*额定转速*/
double ef=0.8;
/*机组效率*/
double iner=1765.8;
/*机组惯量*/
double ss;
/*比转速*/
double mn;
/*正常工作时转矩*/
double es=0.0;
/*上游水库水位*/
double el=39.4;
/*下游水库水位*/
double c1;
/*计算中涉及到的中间变量*/
double resis0;
/*阀门全开时的阻力系数*/
double a[50];
double z[50][1000];
int n[50];
struct pipe_parameter
{int n;
double length;
double diameter;
double thickness;
double friction;
double angle;
};
struct pipe_parameter pi_p[2];
double c3;
double b[50],c[50],r[50];
double flow[50][1000];
double head[50][1000];
double flowp[50][1000];
double headp[50][1000];
double vo[2000];
int n3;/*数组成员个数*/
double dtvo;/*阀门开度特性时间间隔*/
double dh;
double rf[2000];
double rs[2000];
double rh[2000];
double rm[2000];
double x;
int y;
double a0,a1,b0,b1;
/*读入wh[]wm[]数组*/
double wh[89]={0.634,0.643,0.646,0.64,0.629,0.613,0.595,0.575,0.552,0.533,0.516,0.505,0.504,0.510,0.512,0.522,0.539,0.559,0.580,0.601,0.63,0.662,0.692,0.722,0.753,0.782,0.808,0.832,0.857,0.879,0.904,0.93,0.959,0.996,1.027,1.060,1.090,1.124,1.165,1.204,1.238,1.258,1.271,1.282,1.288,1.281,1.260,1.225,1.172,1.107,1.031,0.942,0.842,0.733,0.617,0.5,0.368,0.24,0.125,0.011,-0.102,-0.168,-0.255,-0.342,-0.423,-0.494,-0.556,-0.620,-0.655,-0.67,-0.66,-0.655,-0.64,-0.6,-0.57,-0.52,-0.47,-0.43,-0.36,-0.275,-0,16,-0.04,0.13,0.295,0.43,0.55,0.62,0.634};
double wm[89]={-0.684,-0.547,-0.414,-0.292,-0.187,-0.105,-0.053,-0.012,0.042,0.097,0.156,0.227,0.3,0.371,0.444,0.522,0.596,0.672,0.738,0.763,0.797,0.837,0.865,0.883,0.886,0.877,0.859,0.838,0.804,0.758,0.703,0.645,0.583,0.52,0454,0.408,0.37,0.343,0.331,0.329,0.338,0.354,0.372,0.405,0.45,0.486,0.52,0.552,0.579,0.603,0.616,0.617,0.606,0.582,0.546,0.5,0.432,0.36,0.288,0.214,0.123,0.037,-0.053,-0.161,-0.248,-0.314,-0.372,-0.58,-0.74,-0.88,-1,-1.12,-1.25,-1.37,-1.49,-1.59,-1.66,-1.69,-1.77,-1.65,-1.59,-1.52,-1.42,-1.32,-1.23,-1.1,-0.98,-0.82,-0.684};
int n4=89;/*数组成员个数*/
double dangel;
int yp;
double wh0,wm0;
double cp,cm;
double f1,f2,f1f,f1s,f2f,f2s;
double df,ds;
int d;
double resis;
double rt;
/*读入管道参数*/
pi_p[0].n=2;
pi_p[0].length=27.0;
pi_p[0].diameter=0.7;
pi_p[0].thickness=0.014;
pi_p[0].friction=0.0;
pi_p[0].angle=PI/2;
pi_p[1].n=1;
pi_p[1].length=273.0;
pi_p[1].diameter=1.0;
pi_p[1].thickness=0.014;
pi_p[1].friction=0.0;
pi_p[1].angle=atan((39.4-27.0)/273);
printf("flow0=%7.2f,head0=%7.2f\n",flow0,head0);
/*计算常数*/
ss=3.65*n0*sqrt(flow0)/(pow(head0,0.75));
/*计算比转速*/
mn=(60*1000*G*flow0*head0)/(2*PI*n0*ef);
/*计算转矩*/
printf("ss=%7.2f,mn=%7.2f\n",ss,mn);
c1=es-el;
printf("c1=%7.2f\n",c1);
for(tsum=0,i=0;i<n2;i++)
{ a[i]=PI*pi_p[i].diameter*pi_p[i].diameter/4;
ws[i]=1435/(sqrt(1+E_W*pi_p[i].diameter/(E_S*pi_p[i].thickness)));
t[i]=pi_p[i].length/ws[i];
tsum=tsum+t[i];
}
for(i=0;i<n2;i++)
{int I;
float ws1,ws2;
do
{dt=tsum/num;
I=(int)(pi_p[i].length/dt/ws[i]);
ws1=pi_p[i].length/dt/I;
ws2=pi_p[i].length/dt/(I+1);
num=num+1;
}
while((fabs(ws1-ws[i])>(err*ws[i]))&&(fabs(ws2-ws[i])>(err*ws[i])));
if(fabs(ws1-ws[i])<fabs(ws2-ws[i]))
ws[i]=ws1;
else ws[i]=ws2;
n[i]=pi_p[i].length/ws[i]/dt;
printf("%f,%d\n",ws[i],n[i]);
}
num=num-1;
printf("num=%d\n",num);
printf("dt=%f\n",dt);
c3=iner*n0*PI/(60*G*mn*dt);
printf("c3=%f\n",c3);
z[0][0]=es;
for(i=0;i<n2;i++)
{dx[i]=pi_p[i].length/n[i];
for(j=1;j<=n[i];j++)
{z[i][j]=z[i][j-1]+dx[i]*sin(pi_p[i].angle);
printf("%f\n",z[i][j]);
}
z[i+1][0]=z[i][n[i]];
}
for(i=0;i<n2;i++)
{b[i]=ws[i]/G/a[i];
c[i]=dt/a[i]*sin(pi_p[i].angle);
r[i]=pi_p[i].friction*dx[i]/2/G/pi_p[i].diameter/a[i]/a[i];
printf("%f,%f,%f\n",b[i],c[i],r[i]);
}
k=0;
T=0;
printf("k=0,T=0\n");
/*计算管道稳态参数*/
head[n2][0]=el;
for(i=n2-1;i>=0;i--)
{for(j=n[i];j>=0;j--)
{flow[i][j]=n1*flow0/pi_p[i].n;
head[i][j]=head[i+1][0]+(n[i]+1)*r[i]*flow[i][j]*fabs(flow[i][j]);
}
}
rf[0]=1;
rs[0]=1;
rf[1]=1;
rs[1]=1;
dangel=2*PI/(n4-1);
x=PI+atan(rf[0]/rs[0]);
y=(int)(x/dangel+1);
yp=y;
/*定义函数计算无量纲参数*/
a1=(wh[y+1]-wh[y])/dangel;
a0=wh[y+1]-y*a1*dangel;
b1=(wm[y+1]-wm[y])/dangel;
b0=wm[y+1]-y*b1*dangel;
wh0=a0+a1*x;
wm0=b0+b1*x;
rm[0]=(rf[0]*rf[0]+rs[0]*rs[0])*wm0;
/*计算暂态参数*/
for(k=1,T=dt;k<1000;k++,T=T+dt)
{
for(i=0;i<n2;i++)/*定义函数：计算中间节点的瞬态参数*/
{
for(j=1;j<n[i];j++)
{cp=head[i][j-1]+(b[i]+c[i])*flow[i][j-1]-r[i]*flow[i][j-1]*fabs(flow[i][j-1]);
cm=head[i][j+1]-(b[i]-c[i])*flow[i][j+1]+r[i]*flow[i][j+1]*fabs(flow[i][j+1]);
headp[i][j]=(cp+cm)/2;
if(headp[i][j]<=z[i][j]-10+0.5)
headp[i][j]=z[i][j]-10+0.5;
flowp[i][j]=(cp-headp[i][j])/b[i];
}
}
/*计算管道最末端瞬态参数*/
j=n[n2-1];
headp[n2-1][j]=el;
cp=head[n2-1][j-1]+(b[n2-1]+c[n2-1])*flow[n2-1][j-1]-r[n2-1]*flow[n2-1][j-1]*fabs(flow[n2-1][j-1]);
flowp[n2-1][j]=(cp-headp[n2-1][j])/b[n2-1];
for(i=1;i<n2;i++)/*定义函数：计算管道相接处瞬态参数*/
{j=n[i-1];
cp=head[i-1][j-1]+(b[i-1]+c[i-1])*flow[i-1][j-1]-r[i-1]*flow[i-1][j-1]*fabs(flow[i-1][j-1]);
cm=head[i][1]-(b[i]-c[i])*flow[i][1]+r[i]*flow[i][1]*fabs(flow[i][1]);
headp[i-1][j]=(pi_p[i-1].n*cp/b[i-1]+pi_p[i].n*cm/b[i])/(pi_p[i-1].n*(1/b[i-1])+pi_p[i].n*(1/b[i]));
if(headp[i-1][j]<=z[i-1][j]-10+0.5)
headp[i-1][j]=z[i-1][j]-10+0.5;
headp[i][0]=headp[i-1][j];
flowp[i-1][j]=(cp-headp[i-1][j])/b[i-1];
flowp[i][0]=pi_p[i-1].n*flowp[i-1][j]/pi_p[i].n;
}
/*定义函数计算阀门水头损失*/
/*读入开度数组vo[100]*/
vo[k]=1.0;
resis0=0;
dtvo=1.0;
d=(int)(T/dtvo+1);
if(T<Tm)
rt=vo[d]+(vo[d+1]-vo[d])*(T/dtvo+1-d);
else
rt=0;
resis=resis0/rt/rt;
dh=resis*rf[k]*fabs(rf[k])*flow0*flow0/(2*G*a[1]*a[1]);
/*定义函数：泵处瞬态参数*/
do{y=yp;
a1=(wh[y+1]-wh[y])/dangel;
a0=wh[y+1]-y*a1*dangel;
b1=(wm[y+1]-wm[y])/dangel;
b0=wm[y+1]-y*b1*dangel;
wh0=a0+a1*x;
wm0=b0+b1*x;
for(i=0;i<1000;i++)
{
cm=head[0][1]-(b[0]-c[0])*flow[0][1]+r[0]*flow[0][1]*fabs(flow[0][1]);
f1=es-cm+head0*(rf[k]*rf[k]+rs[k]*rs[k])*(a0+a1*x)-b[0]*flow0*rf[k]-dh;
f2=(rf[k]*rf[k]+rs[k]*rs[k])*(b0+b1*x)+rm[k-1]-c3*(rs[k-1]-rs[k]);
f1f=head0*(2*rf[k]*(a0+a1*x)+a1*rs[k])-2*resis*fabs(rf[k])*flow0*flow0/(2*G*a[1]*a[1])-b[0]*flow0;
f1s=head0*(2*rs[k]*(a0+a1*x)-a1*rf[k]);
f2f=2*rf[k]*(b0+b1*x)+b1*rs[k];
f2s=2*rs[k]*(b0+b1*x)-b1*rf[k]+c3;
ds=(f2/f2f-f1/f1f)/(f1s/f1f-f2s/f2f);
df=-f1/f1f-ds*f1s/f1f;
rf[k]=rf[k]+df;
rs[k]=rs[k]+ds;
if(rs[k]>=0)
x=PI+atan(rf[k]/rs[k]);
else
if(rf[k]<=0)
x=atan(rf[k]/rs[k]);
else
x=2*PI+atan(rf[k]/rs[k]);
yp=(int)(x/dangel+1);
if((fabs(df)+fabs(ds))<0.0001)
break;
}
}
while(yp!=y);
flowp[0][0]=rf[k]*flow0;
headp[0][0]=cm+b[0]*flow[0][0];
if(headp[0][0]<=z[0][0]-10+0.5)
headp[0][0]=z[0][0]-10+0.5;
rh[k]=headp[0][0]/head0;
rm[k]=(rs[k]*rs[k]+rf[k]*rf[k])*(b0+b1*x);
printf("k=%d,T=%f\n",k,T);
printf("%f,%f,%f,%f\n",rf[k],rs[k],rh[k],rm[k]);
rf[k+1]=2*rf[k]-rf[k-1];
rs[k+1]=2*rs[k]-rs[k-1];
if(rs[k+1]>=0)
x=PI+atan(rf[k+1]/rs[k+1]);
else
if(rf[k+1]<=0)
x=atan(rf[k+1]/rs[k+1]);
else
x=2*PI+atan(rf[k+1]/rs[k+1]);
yp=(int)(x/dangel+1);
for(i=0;i<n2;i++)
for(j=0;j<n[i];j++)
{
head[i][j]=headp[i][j];
flow[i][j]=flowp[i][j];
printf("%f,%f\n",flowp[i][j],headp[i][j]);
}
}
getch();
}
2008-05-12 11:22
