24小时热门版块排行榜

返回列表

我原是我以为

新虫 (小有名气)

应助: 0 (幼儿园)
金币: 343.8
帖子: 163
在线: 84.1小时
虫号: 4254249

[交流] 陶老师数值传热学4.8节长方形通道对流换热相关问题

fRe结果跟书上差不多，但是Nu差别好大，有大神帮我看看嘛

#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客户端