//Arthor：吴昱驹
//e-mail：328522073@https://www.doczj.com/doc/ea4611166.html,
//date：2013/01/06
//四点线性隐格式计算法进行洪水演算
//由河海大学计算水力学王船海老师的Fortran程序翻译成C语言
//#include
#include
#include
#include
#define NSR 30
#define PI 3.14159
int N,NT,IC,ID,IT;
int L1,L2,IB;
double ZU,ZD,X;
double DX,DT,CT,CN0;
double Z1,Z2,Q1,Q2,A1,A2,B1,B2,R1,R2,U1,U2;
double P,V,DTT;
double RDX,SA,SB,BC,C,D,E,F,G,O;
double Y1,Y2,Y3,Y4,W,S,T;
// 水位 流量 河段长 P V S T 糙率
double ZZ[NSR],QQ[NSR],DXX[NSR],PP[NSR],VV[NSR],SS[NSR],TT[NSR],CNO[NSR];
double MT;
void input();
void output();
void EXTERN();
void BAEXT();
void ABR(int i,double Z,double Q,double *B,double *A,double *R,double *U);
int main(int argc, char *argv[])
{
int i;
input();
for(i=1;i<=N;i++)
{
DXX[i]=DX*1000;
CNO[i]=CN0;
ZZ[i]=4.0+0.025*(N+1-i);
QQ[i]=0.0;
}
MT=60/DT;
DT=DT*60;
DTT=2*DT;
L1=1;L2=N;
if(IC==1){
IB=0;ZU=4.5;
}
else if(IC==2){
IB=1;ZU=500;
}
else{
IB=1;
}
//以上是初始条件
printf("小时数 0 时段 0--> 断面1 断面6 断面11 断面16 断面21 (水位Z)\n");
for(ID=0;ID<=NT;ID++)
{
for(IT=1;IT<=MT;IT++)
{
X=ID+IT/MT;
ZD=4+1.5*sin(PI*X/12);
if(IC==1){
P=ZU;V=0.0;
}
else if(IC==2){
P=ZU;V=0.0;
}
else{
V=20E6/DT;P=V*ZZ[L1];
}
EXTERN();
ZZ[L2]=ZD;
QQ[L2]=P-V*ZD;
BAEXT();
output();
}
}
//for(i=1;i<=21;i++)
//printf("%lf\n",QQ[i]);
system("pause");
return 0;
}
void input()
{
N=21;NT=48;
IC=1;
DX=1.0;DT=10.0;CT=0.75;CN0=0.02;
}
void EXTERN()
{
int j;
Z1=ZZ[L1];
Q1=QQ[L1];
ABR(L1,Z1,Q1,&B1,&A1,&R1,&U1);
PP[L1]=P;
VV[L1]=V;
for(j=L1+1;j<=L2;j++)
{
Z2=ZZ[j];
Q2=QQ[j];
ABR(j,Z2,Q2,&B2,&A2,&R2,&U2);
DX=DXX[j-1];
RDX=DX*4.905*CNO[j]*CNO[j]/CT;
SA=RDX*(abs(U1)+0.01)/(pow(R1,1.333));
SB=RDX*(abs(U2)+0.01)/(pow(R2,1.333));
BC=(B1+B2)*0.5;
C=DX/DTT*BC/CT;
D=C*(Z1+Z2)-(1.0-CT)/CT*(Q2-Q1);
E=DX/DTT/CT-U1+SA;
G=DX/DTT/CT+U2+SB;
F=9.81*(A2+A1)*0.5;
O=DX/DTT/CT*(Q1+Q2)-(1-CT)/CT*(U2*Q2-U1*Q1);
O=O-(1-CT)/CT*F*(Z2-Z1);
if(IB==0){
Y1=D-C*P;
Y2=O+F*P;
Y3=1+C*V;
Y4=E+F*V;
W=C*Y4+F*Y3;
S=(C*Y2-F*Y1)/W;
T=(C*G-F)/W;
P=(Y1+Y3*S)/C;
V=(1+T*Y3)/C;
}
else{
Y1=C+V;
Y2=F+E*V;
Y3=O-E*P;
Y4=D+P;
W=G*Y1+Y2;
S=(G*Y4-Y3)/W;
T=(C*G-F)/W;
P=Y4-Y1*S;
V=C-Y1*T;
}
PP[j]=P;
VV[j]=V;
SS[j]=S;
TT[j]=T;
Z1=Z2;
Q1=Q2;
U1=U2;
B1=B2;
A1=A2;
R1=R2;
}
if(IB==0){
P=P/V;
V=1/V;
}
}
void BAEXT()
{
int j;
if(IB==0){
for(j=L2-1;j>=L1;j--)
{
QQ[j]=SS[j+1]-TT[j+1]*QQ[j+1];
ZZ[j]=PP[j]-VV[j]*QQ[j];
}
}
else{
for(j=L2-1;j>=L1;j--)
{
ZZ[j]=SS[j+1]-TT[j+1]*ZZ[j];
QQ[j]=PP[j]-VV[j]*ZZ[j];
}
}
}
void output()
{
printf("小
时数%2d 时段%2d--> %4.2lf %4.2lf %4.2lf %4.2lf %4.2lf\n",ID,IT,ZZ[1],ZZ[6],ZZ[11],ZZ[16],ZZ[21]);
}
void ABR(int i,double Z,double Q,double *B,double *A,double *R,double *U)
{
double DH;
//ZD=0.1*(21-i);
//DH=Z-ZD;
DH=Z-0.1*(21-i);
*B=100+6*DH;
*A=(*B+100)*DH/2;
*R=(*A)/(*B);
if(*R<=0.1)
*R=0.1;
*U=Q/(*A);
}