计算水力学代码

//Arthor：吴昱驹//e-mail：[e-mail]//date：2013/01/06//四点线性隐格式计算法进行洪水演算//由河海大学计算水力学王船海老师的Fortran程序翻译成C语言//#includeWindows.h#includestdio.h#includestdlib.h#includemath.h#defineNSR30#definePI3.14159intN,NT,IC,ID,IT;intL1,L2,IB;doubleZU,ZD,X;doubleDX,DT,CT,CN0;doubleZ1,Z2,Q1,Q2,A1,A2,B1,B2,R1,R2,U1,U2;doubleP,V,DTT;doubleRDX,SA,SB,BC,C,D,E,F,G,O;doubleY1,Y2,Y3,Y4,W,S,T;//水位流量河段长PVST糙率doubleZZ[NSR],QQ[NSR],DXX[NSR],PP[NSR],VV[NSR],SS[NSR],TT[NSR],CNO[NSR];doubleMT;voidinput();voidoutput();voidEXTERN();voidBAEXT();voidABR(inti,doubleZ,doubleQ,double*B,double*A,double*R,double*U);intmain(intargc,char*argv[]){inti;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;}elseif(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;}elseif(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);return0;}voidinput(){N=21;NT=48;IC=1;DX=1.0;DT=10.0;CT=0.75;CN0=0.02;}voidEXTERN(){intj;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;}}voidBAEXT(){intj;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];}}}voidoutput(){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]);}voidABR(inti,doubleZ,doubleQ,double*B,double*A,double*R,double*U){doubleDH;//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);}