python中dtt函数,dt python

没有编程基础，学PYTHON前要不要看一下C语言

不需要。

#includestdio.h

#includemath.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(i1)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(ttmax)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(dmdt0.);

    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(x20) 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(iter3)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(iter5)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).001abs(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(iter4)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-k0) 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;

}