| 查看: 750 | 回复: 0 | |||
[求助]
求陶老师数值传热学4.8节长方形截面对流换热C程序，或者帮忙看看程序哪里不对
|
|
求大神帮忙看看loop2，求Nu过程，哪里不对，fRe和书上差不多 #include <math.h> #include <stdio.h> #define NT 70 void main() { /*定义变量*/ int N,M,i,j,Iter; double W[NT][NT],W0[NT][NT],sita[NT][NT],sita0[NT][NT],T[NT][NT],De,Re,Nu; double dltx,dlty,dx,dy,Wm,D,x,y,a,b; double Eps,DT,DTmax,SD; double ap,aw,ae,as,an,JF,b1,b2,JA,ap1; FILE*fp; /*读入参数*/ printf("请输入x方向节点数N y方向节点数M\n"  ;scanf("%d%d",&N,&M); printf("请输入界面长度a 界面宽度b\n" ;scanf("%lf%lf",&a,&b); /*已知参数*/ SD=0; JF=0; dltx=a; dlty=b; De=2*dltx*dlty/(dltx+dlty); D=a;/*特征尺寸*/ N; M; Eps=1.e-6;/*计算精度*/ /*基本参数*/ dx=dltx/(N-1); dy=dlty/(M-1); ae=aw=dy/dx; an=as=dx/dy; /*速度边界条件处理*/ /*固体表面速度为0值*/ for(j=0;j<M;j++) { W[0][j]=0; W[N-1][j]=0; } for(i=1;i<N-1;i++) { W[0]=0; W[M-1]=0; } /*速度场赋初值*/ for(j=1;j<M-1;j++) for(i=1;i<N-1;i++) { W[j]=0.00008; W0[j]=W[j];/*对上一次迭代值赋值*/ } Iter=0; /*Gauss-Seidel迭代计算速度*/ loop1: for(j=1;j<M-1;j++) { for(i=1;i<N-1;i++) { ap=ae+aw+an+as; b1=(dx*dy)/(D*D); W[j]=ae/ap*W[i+1][j]+aw/ap*W[i-1][j]+an/ap*W[j+1]+as/ap*W[j-1]+b1/ap; } } /*计算两次迭代最大误差*/ Iter=Iter+1; DTmax=0.0; for(j=1;j<M-1;j++) for(i=1;i<N-1;i++) { DT=fabs(W[j]-W0[j]); W0[j]=W[j]; if(DT>DTmax) DTmax=DT; } if(DTmax>Eps) goto loop1; for(j=1;j<M-1;j++) for(i=1;i<N-1;i++) { SD+=W[j]*dx*dy; } Wm=SD/(a*b); printf("无量纲速度平均值:\n" ;printf("Wm=%5.6f\n",Wm); /*边界条件处理*/ /*固体表面无量纲为0,取T=JA*sita进行迭代求JA*/ for(j=0;j<M;j++) { sita[0][j]=0; sita[N-1][j]=0; } for(i=1;i<N-1;i++) { sita[0]=0; sita[M-1]=0; } /*sita场赋初值*/ for(j=1;j<M-1;j++) for(i=1;i<N-1;i++) { sita[j]=0.003; sita0[j]=sita[j];/*对上一次迭代值赋值*/ } Iter=0; /*Gauss-Seidel迭代计算*/ loop2: for(j=1;j<M-1;j++) { for(i=1;i<N-1;i++) { JF+=W[j]*sita[j]; } } JA=a*b/JF/Wm; for(j=1;j<M-1;j++) for(i=1;i<N-1;i++) { b2=dx*dy*JA*W[j]/(D*D*Wm); ap1=ap-b2; sita[j]=ae/ap1*sita[i+1][j]+aw/ap1*sita[i-1][j]+an/ap1*sita[j+1]+as/ap1*sita[j-1]; } /*计算两次迭代最大误差*/ Iter=Iter+1; DTmax=0.0; for(j=1;j<M-1;j++) for(i=1;i<N-1;i++) { DT=fabs(sita[j]-sita0[j]); sita0[j]=sita[j]; if(DT>DTmax) DTmax=DT; } if(DTmax>Eps) goto loop2; /*Re计算*/ printf("Re数为:\n" ;Re=2*(De/D)*(De/D)/Wm; printf("%5.6f\n",Re); /*Nu计算*/ printf("Nu数为:\n" ;Nu=0.25*(De/D)*(De/D)*JA; printf("%5.6f\n",Nu); for(j=1;j<M-1;j++) for(i=1;i<N-1;i++) { T[j]=JA*sita[j]; } fp=fopen("result.dat","w" ;printf("无量纲温度:\n" ;for(j=M-1;j>=0;j--) { for(i=0;i<N;i++) { fprintf(fp,"%5.6f",T[j]); } printf("\n" ;} fclose(fp); } 发自小木虫IOS客户端 |
回复此楼» 猜你喜欢
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有3人回复
-
基底STO,薄膜SRO，XRD里面的振荡，是laue震荡还是kiessig振荡? 怎么判断？
已经有2人回复
-
物理学I论文润色/翻译怎么收费?
已经有224人回复
-
B区学生调剂-兰州交通大学材料科学与工程学院
已经有12人回复
-
山西大同大学物理学专业还有调剂名额，欢迎调剂！
已经有13人回复
-
桂林理工大学物理学专业招收调剂
已经有14人回复
-
VASP 的一组 GPU / CPU 基准测试记录
已经有0人回复
-
美国顶级物理期刊《应用物理快报》（APL）的编辑欺骗和歧视作者及AIP的官僚主义傲慢
已经有7人回复
-
津理工大学晶体材料全国重点实验室刘红军教授课题组招收博士生一名
已经有0人回复
-
【原创讨论】从电子约束到物质编辑：一套可迭代的环形磁场科技树
已经有0人回复
