ä¸éè¦ã
#include<stdio.h>#include<math.h>
#define PI 3.14159
int main(){
double hap[42],qp[42],ha[42],q[42],qpu[42],qu[42],vcav[42],hapr[1601],icav[42],c3[42],am[42],zeta[42],dz[42],r[42],har[42],has[42];
double xl=1000.;
double d=7.;
double f[42]={0,.015,.015,.015,.415,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015,.015};
double g=32.2;
double rp=.000137;
double q0=320.;
double ck=4.77e7;
double e=4.32e9;
double w=.0833;
double rg=1716.;
double te=570.;
double rho=1.923;
double hb=33.5;
double amm=2.16e-5;
double amd=4.e-5;
double hs0=33.5;
double cmk=3.e-8;
double el[42]={0,5.,5.,5.,28.,28.,5.,5.,5.,5.,5.,5.,5.,5.,5.,5.};
double hv=2.95;
double hh[35]={0,60.694,55.75,51.08,46.96,43.39,40.37,37.9,35.97,34.5,
33.78,33.5,33.78,34.5,35.97,37.9,40.37,43.39,46.96,
51.08,55.75,60.964};
int nh=21;
double dtt=.5;
int iii=5;
double tmax=13.1;
int n=10;
int ipr=18;
int kit=5;
double tol=.05;
int n2,iter,i,k,ns,i1,j;
double x3,c2,b,po,dmdt,dt,xpp,x6,x2,qr,th,hhh,x1,ar,ff,t,th1,c1,z1,cm2,x4,z2,xp,cp,dfdh,gam,cm,x5,aa,qs,dh,s;
FILE *fp;
fp=fopen("cprintf.txt","a+");
//! namelist/din/xl,d,f,g,rp,el,q0,hv,ck,e,w,rg,te,rho,hh,amd,hb;
//! 2 amm,hs0,cmk,tol,tmax,n,ipr,kit,nh,dtt,iii;
//! 10 read(5,din,fclose(fp);return 0=99);
ns=n+1; //!模æç管éååç段æ°ï¼nsæ¯èç¹æ°
n2=n/10; //!n2çå«ä¹æªç¥
ar=.7854*d*d; //! aræ¯ç®¡éçé¢ç§¯
gam=rho*g; //!æ¯é
c1=1.+ck*d/(e*w); //! å ¬å¼8-6ä¸ç1+kl*d/ee;
c2=rg*te*ck/(c1*gam*gam); //!111页ä¸çc3=c2*m/rçc2*/r/rï¼c2ä½8-7ä¹å³æ¯è¯´å¯¹åºä¹¦ä¸çæ¯c3;
//!aa对åºçæ¯ä¹¦ä¸çå ¬å¼8-6;
aa=sqrt(ck/(rho*c1)); //! å ¬å¼8-6;
b=aa/(g*ar);
//!æ¶é´é´é
dt=xl/(aa*n); //!对åºçæ¶æ¥é¿
cm2=(hs0-hv)/amd; //!ç»å¯¹è¸æ±½ååï¼è¿ä¸ªå ¬å¼å°åºæ¯æä¹å¾å°çï¼
//!hhè§å®ä¸ºæ¶é´çå½æ°ï¼ç¨æ¥è¡¨ç¤ºç±æ³µäº§ççç»å¯¹å头
//!zæ¯åºæ°´ç®¡çé«åº¦
//!rpæ¯å¨æ³µçä¸æ¸¸ç®¡ç³»ä¸ï¼ä»å ¥å£å°æªé¢1é´å æ¬å ¨é¨æ失çä¸ä¸ªæ失系æ°ã
ha[1]=hh[1]-rp*q0*q0-el[1];
for(i=1;i<=ns;i++){ //20
dz[i]=el[i+1]-el[i];
r[i]=f[i]*xl/(2.*g*d*ar*ar*n);
if(i>1)ha[i]=ha[i-1]-r[i-1]*q0*q0-dz[i-1]; //!稳ææ¶çå个æªé¢ä¸çå头
hap[i]=ha[i];
har[i]=ha[i];
has[i]=ha[i];
q[i]=q0;
qu[i]=q0;
//!没æèªç±ç©ºæ°æ¶ï¼a'=aå³å ¬å¼8-13ä¸çc3=0;
c3[i]=0.;
//!èµ·å§æ¶å»çæå¼ç³»æ°ä¸º1ï¼å³ä¸ç¨æå¼ï¼ä½ç¨æ¯èµåå¼
zeta[i]=1.;
//!icav[i]åvcav[i]ç¨æ¥æ§å¶æ°ç©´ççï¼å¦æ没ææ°è¡ç¼ 身å置为0;
icav[i]=0;
vcav[i]=0.;
l20: am[i]=0.; //!am[i]表示çæ¯c3=c2*m/r/rä¸çmï¼èµ·å§æ¶å»ç空æ°çè´¨é为0;
}
hapr[1]=ha[iii];
th=dt*n/xl;
//!????è¿ä¸ªå ¬å¼çç©çæä¹æ¯ä»ä¹å¢ï¼xlæ¯ç®¡çº¿æ»çé¿åº¦ ï¼næ¯ç®¡çº¿åæç段æ°ï¼dtæ¯æ¶é´é´é,thæ¯æ³¢éçåæ°ï¼ä½æ¯ç©çæä¹æ¯ä»ä¹å¢ï¼
t=0.;
k=0;
fprintf(fp,"a,xl,d,f=%8.1lf%8.1lf%8.4lf%8.4lf\n",aa,xl,d,f[1]);
fprintf(fp,"h0,q0,rp=%8.2lf%8.2lf%8.5lf\n",ha[ns],q0,rp);
fprintf(fp,"ck,e,w =%11.4e%11.4e%8.3lf\n",ck,e,w);
fprintf(fp,"hv,rg,te,rho=%9.3lf%9.3lf%9.3lf%8.2lf\n",hv,rg,te,rho);
fprintf(fp,"amd,amm,cmk,hs0=%11.4e%11.4e%11.4e%8.3lf\n",amd,amm,cmk,hs0);
fprintf(fp,"g,tmax,b,dt,tol=%8.3lf%8.1lf%8.1lf%8.4lf%8.4lf\n",g,tmax,b,dt,tol);
fprintf(fp,"hb,dtt,nh=%8.1lf%8.4lf%4d\n",hb,dtt,nh);
fprintf(fp,"n,ipr,kit,iii%4d%4d%4d%4d\n",n,ipr,kit,iii);
fprintf(fp," pressure heads,discharges,air mass,zeta,and cav size\n\n");
fprintf(fp," time x/l= 0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0\n");
l30: fprintf(fp,"0%7.3lf,ha=",t);
for(i=1;i<=ns;i+=n2)fprintf(fp,"%8.2lf",ha[i]);
fprintf(fp,"\n q=");
for(i=1;i<=ns;i+=n2)fprintf(fp,"%8.3lf",q[i]);
fprintf(fp,"\n am=");
for(i=1;i<=ns;i+=n2)fprintf(fp,"%8.1e",am[i]);
fprintf(fp,"\n z=");
for(i=1;i<=ns;i+=n2)fprintf(fp,"%8.4f",zeta[i]);
fprintf(fp,"\n ,vcav=");
for(i=1;i<=ns;i+=n2)fprintf(fp,"%8.5f",vcav[i]);
fprintf(fp,"\n");
l40: k=k+1;
t=t+dt;
if(t>tmax)goto l145;
for(i=1;i<=n;i++){ //50
if(am[i]==amm)goto l50;
dmdt=cmk*(hs0-cm2*am[i]-.5*(ha[i]+ha[i+1])); //!dmdt表示çå«ä¹æªç¥
if(dmdt>0.);
am[i]=am[i]+dmdt*dt;
if(am[i]>amm) am[i]=amm;
c3[i]=am[i]*c2; //!!!!书ä¸111页çc3;
l50: continue;
}
//c upstream boundary;
i=t/dtt+1; //!dttå°åºæ¯ä»ä¹
if(i>=nh)goto l53; //!nhçå«ä¹æªç¥
th1=(t-(i-1)*dtt)/dtt; //!dttåth1çå«ä¹æªç¥
hhh=hh[i]*(1.-th1)+th1*hh[i+1]-el[1]; //!hhhçå«ä¹æªç¥
goto l54;
l53: hhh=hh[nh]-el[1];
l54: iter=0;
if(rp==0.)hap[1]=hhh;
x4=sqrt(c3[1]);
l55: iter=iter+1;
if(abs(hap[1]-has[1])>.001)goto l56;
zeta[1]=1./sqrt(1.+c3[1]/pow(hap[1],2)); //!æ¤å¤æé®é¢ï¼çç¥çä¸å¯¹ã éè¦éªè¯
goto l59;
l56: x5=sqrt(c3[1]+hap[1]*hap[1]);
xpp=x5-x4*log((x4+x5)/hap[1]);
x5=sqrt(c3[1]+has[1]*has[1]);
zeta[1]=(hap[1]-has[1])/(xpp-x5+x4*log((x4+x5)/has[1]));
l59: has[1]=ha[1]+zeta[1]*(ha[2]-ha[1]);
qs=q[1]+zeta[1]*(qu[2]-q[1]);
x5=sqrt(c3[1]+has[1]*has[1]);
cm=x5-x4*log((x4+x5)/has[1])-qs*(b-r[1]*abs(qs))+dz[1];
if(rp==0.)goto l64;
for(i=1;i<=kit;i++){ //63
s=1.;
x2=hhh-hap[1];
if(x2<0) goto l57;
else if(x2==0) goto l62;
else goto l58;
l57: s=-1.;
l58: if(abs(x2)<1.e-5)goto l62;
x1=s*b*sqrt(s*x2*rp); //!x1ï¼x3æä¹ä¸æç¡®
x3=.5*b/sqrt(s*x2*rp);
goto l61;
l62: x1=0.;
x3=0.;
l61: x5=sqrt(c3[1]+hap[1]*hap[1]);
x6=(x4+x5)/hap[1];
ff=x5-x4*log(x6)-cm-x1;
dfdh=1./x6+x4/hap[1]+x3;
dh=-ff/dfdh;
hap[1]=hap[1]+dh;
if(abs(dh)<tol)goto l64;
l63: if(hap[1]<hv) hap[1]=hv;
}
l64: if(iter<3)goto l55;
x5=sqrt(c3[1]+hap[1]*hap[1]);
qp[1]=(x5-x4*log((x4+x5)/hap[1])-cm)/b;
l65: qpu[1]=qp[1];
//c interior section;
for(i=2;i<=n;i++){ //90
i1=i-1;
iter=1;
x1=sqrt(c3[i1]);
x4=sqrt(c3[i]);
goto l67;
l66: iter=iter+1;
z1=zeta[i1];
z2=zeta[i];
if(iter>5)goto l90; //!!!!æ¯å¤§äºè¿æ¯å°äºæªå®
l67: if(abs(hap[i]-hap[i1])>0.001)goto l68;
zeta[i1]=1./sqrt(1.+c3[i1]/pow(hap[i],2));
goto l69;
l68: x2=sqrt(c3[i1]+hap[i]*hap[i]);
xp=x2-x1*log((x1+x2)/hap[i]);
x2=sqrt(c3[i1]+hap[i1]*hap[i1]);
zeta[i1]=(hap[i]-hap[i1])/(xp-x2+x1*log((x1+x2)/hap[i1]));
l69: if(abs(hap[i]-has[i])>0.001)goto l70;
zeta[i]=1./sqrt(1.+c3[i]/pow(hap[i],2));
goto l72;
l70: x5=sqrt(c3[i]+hap[i]*hap[i]);
xpp=x5-x4*log((x4+x5)/hap[i]);
x5=sqrt(c3[i]+has[i]*has[i]);
zeta[i]=(hap[i]-has[i])/(xpp-x5+x4*log((x4+x5)/has[i]));
l72: if(iter==1||icav[i]==1)goto l74;
if(abs(zeta[i1]-z1)<.001&&abs(zeta[i]-z2)<0.001)goto l90;
l74: hap[i1]=ha[i]-zeta[i1]*(ha[i]-ha[i1]);
qr=qu[i]-zeta[i1]*(qu[i]-q[i1]);
x2=sqrt(c3[i1]+hap[i1]*hap[i1]);
has[i]=ha[i]-zeta[i]*(ha[i]-ha[i+1]);
qs=q[i]-zeta[i]*(q[i]-qu[i+1]);
x5=sqrt(c3[i]+has[i]*has[i]);
cp=x2-x1*log((x1+x2)/hap[i1])+qr*(b-r[i1]*abs(qr))-dz[i1];
cm=x5-x4*log((x4+x5)/has[i])-qs*(b-r[i]*abs(qs))+dz[i];
if(icav[i]==1)goto l85;
l80: for(j=1;j<=kit;j++){ //89
x2=sqrt(c3[i1]+hap[i]*hap[i]);
x3=(x1+x2)/hap[i];
x5=sqrt(c3[i]+hap[i]*hap[i]);
x6=(x4+x5)/hap[i];
ff=x2-x1*log(x3)+x5-x4*log(x6)-cp-cm;
dfdh=1./x3+1./x6+(x1+x4)/hap[i];
dh=-ff/dfdh;
hap[i]=hap[i]+dh;
if(abs(dh)<tol)goto l82;
l89: if(hap[i]<hv) hap[i]=hv-0.0001;
}
l82: if(hap[i]<hv) goto l85;
x2=sqrt(c3[i1]+hap[i]*hap[i]);
qpu[i]=(cp-x2+x1*log((x1+x2)/hap[i]))/b;
qp[i]=qpu[i];
goto l66;
l85: hap[i]=hv;
icav[i]=1;
if(iter==4)goto l86;
if(icav[i1]==0||icav[i+1]==0)goto l66;
l86: x2=sqrt(c3[i1]+hap[i]*hap[i]);
x3=(x1+x2)/hap[i];
x5=sqrt(c3[i]+hap[i]*hap[i]);
x6=(x4+x5)/hap[i];
qpu[i]=(cp-x2+x1*log(x3))/b;
qp[i]=(x5-x4*log(x6)-cm)/b;
vcav[i]=vcav[i]+.5*dt*(qp[i]+q[i]-qpu[i]-qu[i]);
if(vcav[i]>0.)goto l90;
icav[i]=0;
vcav[i]=0.;
goto l80;
l90: continue;
}
//c downstream boundary;
iter=0;
x1=sqrt(c3[n]);
l109: iter=iter+1;
if(abs(hap[ns]-hap[n])>0.001)goto l110;
zeta[n]=1./sqrt(1.+c3[n]/pow(hap[ns],2));
goto l112;
l110: x2=sqrt(c3[n]+hap[ns]*hap[ns]);
xp=x2-x1*log((x1+x2)/hap[ns]);
x2=sqrt(c3[n]+hap[n]*hap[n]);
zeta[n]=(hap[ns]-hap[n])/(xp-x2+x1*log((x1+x2)/hap[n]));
l112: hap[n]=ha[ns]-zeta[n]*(ha[ns]-ha[n]);
if(iter<4)goto l109;
x2=sqrt(c3[n]+hap[n]*hap[n]);
qr=qu[ns]-zeta[n]*(qu[ns]-q[n]);
cp=x2-x1*log((x1+x2)/hap[n])+qr*(b-r[n]*abs(qr))-dz[n];
x2=sqrt(c3[n]+hap[n]*hap[n]);
qp[ns]=(cp-x2+x1*log((x1+x2)/hap[ns]))/b;
qpu[ns]=qp[ns];
for(i=1;i<=ns;i++){ //140
po=ha[i];
ha[i]=hap[i];
hap[i]=2.*ha[i]-po;
if(hap[i]<hv)hap[i]=hv;
q[i]=qp[i];
l140: qu[i]=qpu[i];
}
hapr[k+1]=ha[iii];
if(ha[iii]>hs0)goto l30;
if(k/ipr*ipr-k<0) goto l40;
else if(k/ipr*ipr-k==0) goto l30;
else goto l40;
l145: for(i=1;i<=k;i++){
fprintf(fp,"%6.1lf",hapr[i]);
if(i%10==0)fprintf(fp,"\n");
}
//! goto l10;
l99: fclose(fp);
return 0;
}
温馨提示：内容为网友见解，仅供参考
无其他回答
相似回答
