24小时热门版块排行榜    

查看: 917  |  回复: 15

我原是我以为

新虫 (小有名气)


[交流] 陶老师数值传热学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客户端
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
简单回复
syhorchid2楼
2017-04-11 22:11   回复  
我原是我以为(金币+1): 谢谢参与
2017-04-11 22:18   回复  
我原是我以为(金币+1): 谢谢参与
n 发自小木虫Android客户端
2017-04-11 22:21   回复  
我原是我以为(金币+1): 谢谢参与
2017-04-12 19:43   回复  
2017-04-12 19:48   回复  
我原是我以为(金币+1): 谢谢参与
。。。 发自小木虫Android客户端
tzynew7楼
2017-04-12 19:55   回复  
我原是我以为(金币+1): 谢谢参与
2017-04-12 20:00   回复  
我原是我以为(金币+1): 谢谢参与
xhmaohan9楼
2017-04-12 20:03   回复  
我原是我以为(金币+1): 谢谢参与
发自小木虫Android客户端
2017-04-12 20:08   回复  
我原是我以为(金币+1): 谢谢参与
发自小木虫Android客户端
41588116811楼
2017-04-12 20:12   回复  
我原是我以为(金币+1): 谢谢参与
发自小木虫IOS客户端
hydzp12楼
2017-04-12 20:42   回复  
我原是我以为(金币+1): 谢谢参与
发自小木虫IOS客户端
2017-04-12 21:13   回复  
我原是我以为(金币+1): 谢谢参与
skyish14楼
2017-04-12 21:21   回复  
我原是我以为(金币+1): 谢谢参与
发自小木虫IOS客户端
guanlianwu15楼
2017-04-12 21:32   回复  
我原是我以为(金币+1): 谢谢参与
假大空16楼
2017-04-12 21:59   回复  
我原是我以为(金币+1): 谢谢参与
相关版块跳转 我要订阅楼主 我原是我以为 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 投票:  有多少人是今天查系统知道结果的? +15 爱看书的可乐 2026-08-26 17/850 2026-08-27 00:35 by xxniao123
[基金申请] 基金不中,共勉 +5 eulota 2026-08-26 5/250 2026-08-26 23:33 by young90
[基金申请] 国自然面上复盘~欢迎讨论 +4 晴天加油 2026-08-26 5/250 2026-08-26 21:27 by 木风9012
[基金申请] 2026年8月25日国自然放榜前突然收到列入评审专家邮件,有关系吗? +25 木水思豆 2026-08-25 28/1400 2026-08-26 14:53 by draco1987
[基金申请] 出来了 +9 trojank 2026-08-26 9/450 2026-08-26 14:25 by 宝贝虫子
[基金申请] 只有每年这种时候来逛逛小木虫 +26 yaoyewhu2008 2026-08-20 28/1400 2026-08-26 13:49 by zzuzxg
[基金申请] 哪位高人中了,把查询到的截图贴出来让我看看,让我长长见识 +4 yuleib84 2026-08-26 5/250 2026-08-26 13:40 by yuleib84
[基金申请] 为什么国自然不能直接公布 +4 bjdxyxy 2026-08-26 4/200 2026-08-26 13:12 by qingmu1201
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 国际合作可查了,中了面上 +18 Ldrop2023 2026-08-26 18/900 2026-08-26 11:15 by cmrandy
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 国合可查了 +3 paperzjh 2026-08-26 3/150 2026-08-26 10:41 by LemmonTr
[基金申请] 今天务委会开完了,明天出结果吗 +19 angus9576 2026-08-25 23/1150 2026-08-26 10:03 by zp519
[基金申请] 今日不放榜?网传国自然预计 8 月 27 日可查结果 +17 医学老男孩 2026-08-20 22/1100 2026-08-25 15:36 by 医学老男孩
[基金申请] 建议基金发布提前给出明确的时间点 +13 kulium 2026-08-21 16/800 2026-08-24 16:27 by superceng
[基金申请] 让我中一个面上吧! +13 大萍1987 2026-08-20 16/800 2026-08-24 10:23 by 太傻了
[教师之家] 跳槽后在研项目怎么办? +5 简单化xn 2026-08-22 10/500 2026-08-23 12:38 by 简单化xn
[基金申请] 时间戳今天,20号变了 +5 archvillain 2026-08-20 5/250 2026-08-22 06:12 by hui_daxiao
[基金申请] 应该是下周三26日公布了吧? +4 哈哈蛤? 2026-08-21 4/200 2026-08-21 10:58 by Vivilian
[基金申请] 基金啊基金 +4 longfie172 2026-08-20 4/200 2026-08-21 08:58 by mark mao
信息提示
请填处理意见