24小时热门版块排行榜

返回列表

我原是我以为

新虫 (小有名气)

应助: 0 (幼儿园)
金币: 343.8
散金: 34
红花: 2
帖子: 163
在线: 84.1小时
虫号: 4254249
注册: 2015-11-29
性别: GG
专业: 工程热物理与能源利用

[求助] 求陶老师数值传热学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客户端

回复此楼

» 猜你喜欢

求标准粉末衍射卡号 ICDD 01-076-1802 已经有0人回复
新西兰Robinson研究所招收全奖PhD 已经有0人回复
物理学I论文润色/翻译怎么收费? 已经有260人回复
石墨烯转移--二氧化硅衬底石墨烯已经有0人回复
笼目材料中量子自旋液体基态的证据已经有0人回复
数学教学论硕士可以读数学物理博士吗？已经有0人回复
德国亥姆霍兹Hereon中心汉堡分部招镁合金腐蚀裂变SCC课题方向2026公派博士生已经有4人回复
澳门大学应用物理及材料工程研究院潘晖教授课题组诚招博士后已经有11人回复
求助NH4V4O10晶体的CIF文件已经有0人回复
英国全奖博士招聘-深度学习与量子物理已经有0人回复

1楼 2017-04-12 19:48:44

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

相关版块跳转我要订阅楼主我原是我以为的主题更新

返回列表