计算传热学课程设计报告

1、计算传热学课程设计（报告）题 目:充满多孔介质的长方形截面通道内充分发展对流换热问题的数值研究学生姓名：朱鹏 齐尚超 杨鹏来 芦旭红学 号：10123106 10123107 10123108 10123103专业班级：热能与动力工程101班指导教师：黄善波 巩亮 2013年 7 月 5 日项目每人任务分配编程芦旭红，齐尚超，杨鹏来，朱鹏问题分析朱鹏（图像），齐尚超（程序）图像杨鹏来PPT制作朱鹏推导芦旭红，齐尚超，杨鹏来，朱鹏报告芦旭红热工一班组长：朱鹏 组员：芦旭红，齐尚超，杨鹏来目录1. 设计题目.31.1设计题目.31.2已知参数.42. 物理与数学模型. . .52.1物理模型.52

2、.2数学模型.53. 数值处理与程序设计.6 3.1数学模型无量纲化.6 3.2数值求解.8 3.3程序编写. .114. 程序的验证.125. 计算结果与分析.146. 结论.217. 参考文献.218. 附录.221 设计题目（多孔介质，矩形a/b,单方程）水在一长方形截面的通道中进行充分发展的层流流动，该通道内充满多孔介质。多孔介质具有良好的强化换热能力，孔隙率是其基本结构参数，据此可以计算渗透率K，惯性系数CF，有效导热系数ke,具体表达式见5。其内部充满流体时的流动和换热通常采用体积平均法进行建模，即不考虑区域内孔的微结构而假定区域内任意一点处既有流体相又有固体相。由于金属泡沫的固体

3、骨架导热系数较高，因此对于其内部的对流换热，通常采用局部非平衡模型，即考虑区域内流体温度和固体温度的差异。填充孔隙率为=0.6的多孔介质，渗透率表示为：惯性系数表示为：有效导热系数ke表示为：沿流动方向的速度方程可以简化为 （1）截面上的平均流速为wm=0.1m/s,dp/dz的值是恒定的，可以通过下式得到： （2）其中，w为沿流动方向的速度。换热方程为： （3）其中z为轴向。假设流动和换热都达到充分发展，外壁面为恒热流边界条件（qw=1000W/）,请基于局部热平衡模型，选取a=0.04m,b=0.02m,以20作为水物性的参考温度，参考数值传热学（陶文铨著）4.8节的内容1.1 设计题目1

4、. 分析孔隙率对渗透率K，惯性系数CF，有效导热系数ke的影响规律，计算随的变化关系；2. 通过能量守恒将方程（3）化为更简单的形式进而消去z（仅对恒热流条件下实施），对流动和换热方程进行无量纲化处理；3. 计算轴向无量纲速度的二维数值解，并计算摩擦系数f, 数，分析对f数的影响；4. 计算无量纲温度的的二维数值解，并计算截面流体平均温度和数；1.2 已知参数为了求得数值结果和利用结果进行分析，现给定题目相关已知量，长方形截面:a=0.04m,b=0.02m.填充孔隙率=0.6，热流密度qw=1000W/,平均流速，孔隙直径,20水物性参数（查表知）:=998.2kg/,=。渗透率K（在一定压

5、差下，孔隙结构允许流体通过的能力）惯性系数（用于研究渗流，流体力学中常用的无量纲量）压降（沿轴向方向由于增速降压时引发的压头损失）计算2 物理与数学模型2.1 物理模型图1 长方形截面通道物理模型图 图2 流体坐标2.2 数学模型渗透率表示为： （2-1）惯性系数表示为： （2-2）有效导热系数ke表示为： （2-3） 沿流动方向的速度方程可以简化为 （2-4）相应的边界条件： 截面上的平均流速为=0.1m/s,的值是恒定的，可以通过下式得到： （2-5）其中，w为沿流动方向的速度。换热方程为： （2-6）其中z为轴向。相应的边界条件： 3.数值处理与程序设计3.1数学模型无量纲化(1).流动

6、方程无量纲化处理 （3-1）定义无量纲速度： 定义无量纲坐标：，其中D为通道截面的某一特性尺寸，即D=a，则有： ， ， 将 ， 代入式（3-1）中，得： （3-2） （3-3） （3-4）沿流动方向速度方程边界条件处理（由于紧靠壁面的流体受粘滞力作用而导致靠壁面处侧流速为0）相应的边界条件为： (2).换热方程无量纲化处理 （3-5） ， ， 当量直径 ， ：壁温 ，：截面任意位置处流体的温度 ， ，代入式（3-5）得： （3-6）换热方程的边界条件处理（由于对于研究沿轴向的一个微元小薄层的四周壁温均视为定值，即均为）相应的边界条件为： 3.2 数值求解利用有限容积法对控制方程进行数值求解首

7、先采用外节点法对求解区域进行等距离散(见图3) 图3 对求解区域进行等距离散(1) 流动方程的离散对速度方程在控制容积上积分有： （3-7）其中： ， 源项处理： ， ， 满足0故 （为上一次迭代值）将以上各项代入式（3-7）得： （3-8） （3-9）其中： ， ， ， ，(2) 换热方程的离散对换热方程在控制容积上积分有： （3-10） 源项： ， ，（满足要求）将以上各项代入式（3-10）中，整理的： （3-11） ， ， ， ， (3)计算摩擦系数f, 数增速降压 （3-12） ， （3-13）将以上各量代入（3-12），得： （3-14） （3-15） ， 将以上各量代入（3-14）

8、，得： （3-16） （3-17）(4)推到数 假设四周对长方形通道加热， （3-18） ， ， 将以上各量代入式（3-17）得： （3-19）3.3 程序编写 程序流程图如下：图4程序流程4程序的验证4.1孔隙率=1时的验证当孔隙率=1时，水在长方形截面通道中充满多孔介质的充分发展的层流流动就相应退化成水在长方形截面通道中无多孔介质的充分发展的层流流动。其沿流动方向的速度方程相应的退化成：， 边界条件： 定义无量纲量：定义无量纲速度： ，定义无量纲坐标：，化简为： ，无量纲边界条件： 离散该方程最终可得： ， ， ， ，表一 长方形截面层流充分发展换热的Nu与fRe运行值和精确值的比较（=1

9、）项目运行值精确值误差（%）fRe61.999753620.0004Nu4.1747994.121.334.2网格独立性考核表二网格密度17*921*1131*1641*2151*2661*31Nu8.0027957.9211797.8415927.8144127.8024267.79651图5分析：对于不同网格划分之下的Nu，其值开始减小，但随着网格密度的增加变化逐渐变小，最后趋于稳定。故可认为满足网格独立性考核，<5%5 计算结果与分析5.1速度与无量纲温度分布图6图7图8图9特点：紧靠长方形截面四周壁面流动的流体由于受到粘滞力的作用，速度为0；而稍离四周壁面的流体由于所受影响较小，

10、因而速度基本稳定在平均速度左右。长方形截面通道四壁均受到恒热流加热，对于沿轴向方向截面的一个微元薄层中认为四周壁温均等，总体来说，内部流体受热少，无量纲温度越往里越来越负，无量纲温度呈现出对称性。52分析孔隙率对渗透率K，惯性系数CF，有效导热系数ke的影响规律，计算随的变化关系。l 孔隙率对渗透率K的影响 () 表三 孔隙率与渗透率K数据表00.050.10.150.20.250.30.35K09.23E-138.23E-123.11E-118.33E-111.85E-103.67E-106.77E-100.40.450.50.550.60.650.70.75K1.19E-092.01E-0

11、93.33E-095.48E-099E-091.49E-082.54E-084.5E-080.80.850.90.950.970.99K8.53E-081.82E-074.86E-072.29E-066.76E-066.47E-05图10特点： 1）.当为0时，说明该通道完全被阻隔，因而渗透能力为0，渗透率为0；当为1时，说明该通道完全不被阻隔，即认为不存在任何孔隙结构，渗透率趋于无穷大，此时谈渗透率无意义。 2）.从函数关系及图像上看，上图开始时变化趋势缓慢增加，但趋近于1时，其变化趋势猛然急增并趋于无穷大。这充分说明=1时为K=K（）的一个阶跃点。l 孔隙率对渗透率的影响 ()表四 孔隙率

12、与惯性系数数据表0.010.050.10.150.20.250.30.350.40.450.5CF3.0784051.05280.6632220.5061330.4178040.3600520.3188440.2877050.26320.2433240.2268190.550.60.650.70.750.80.850.90.951CF0.2128550.2008590.1904220.1812430.1730950.1658060.1592380.1532840.1478570.142887图11特点：1）.孔隙率接近0时，惯性系数趋于无穷大。但随着孔隙率的增加，越来越小。 2）.孔隙率接近于

13、0时，图像变化趋势很大。接近于1时，图像变化趋势较慢。l 孔隙率对渗透率的影响 ()表五 孔隙率与有效导热系数数据表00.050.10.150.20.250.3ke（W/m*K)383.8364.64345.4799326.3199307.1598287.9998268.83970.350.40.450.50.550.60.65ke（W/m*K)249.6797230.5196211.3596192.1995173.0395153.8794134.71940.70.750.80.850.90.951ke（W/m*K)115.559396.3992577.239258.0791538.91911

14、9.759050.599图12特点：在其他条件不发生变化时，与大致成线性关系，并且随的增大而逐渐减小。当时，；当时，。l 计算随的变化关系 ()表六 与的数据表wm00.10.20.30.40.50.6dp/dz0-32199.5-106487-222862-381325-581875-824514wm0.70.80.911.11.21.3dp/dz-1109240-1436054-1804956-2215946-2669023-3164189-3701442wm1.41.51.61.71.81.92dp/dz-4280783-4902212-5565728-6271333-7019025-7

15、808805-8640673图13特点：1）.从函数关系上看，当其他条件不变时，与成二次函数关系，并且随的增加而单调递减，也就是说平均流速越大，其压降也就越大。 2）. 的图像开始变化缓慢，之后变化趋势增加。5.3分析对f数的影响 ()表七 与f数据表 dp/dz0-32199.5-106487-222862-381325-581875-824514fRe0456123.4150844431569625401678824259111679701dp/dz-1109240.173-1436054-1804956-2215946-2669023-3164189-3701442fRe15713008

16、.23203425132556821531390115378082124482250652432998dp/dz-4280782.909-4902212-5565728-6271333-7019025-7808805-8640673fRe60639687.15694425747884165788836939994284171.11E+081.22E+08图14特点：1）.其他条件不变时，与之间成线性关系且为正比关系，随的增大而增大。 2）.对于在其他条件不变时，同一个平均流速下，其为定值。5.4 lnNu随lnb/a变化关系表八 lnNu与lnb/a数据表ln(b/a)ln(Nu)-4.605

17、1701862.559106707-2.3025850932.371209443-1.3862943612.214954033-0.6931471812.054434711-0.2876820721.98295445701.9658375580.2231435511.9766396940.4054651081.9999553250.5596157882.0275695270.6931471812.0557723741.6094379122.2637028832.3025850932.3738002134.6051701862.508669725图15特点：1）.从图像上看，在其他条件不变时，l

18、n()随着ln(b/a)的变化呈轴对称变化。 2）.从函数关系上考虑，当即ln(b/a)=0时，最小。时可以认为此=的一个极值点。5.5 Nu随孔隙率变化规律表九Nu与孔隙率数据表孔隙率 0.60.70.80.91Nu7.8026.2595.2594.6254.174图16特点：根据图像，Nu数随孔隙率的增大逐渐减小。（说明取不同的孔隙率的时候所对应的空隙直径也不相同）6 总结 1.速度与无量纲温度分布的结论：紧靠长方形截面四周壁面流动的流体由于受到粘滞力的作用，速度为0；而稍离四周壁面的流体由于所受影响较小，因而速度基本稳定在平均速度左右。长方形截面通道四壁均受到恒热流加热，对于沿轴向方向截

19、面的一个微元薄层中认为四周壁温均等，总体来说，内部流体受热少，无量纲温度越往里越来越负，无量纲温度呈现出对称性。 2.孔隙率的影响的结论：在其他条件不变的情况下，随着孔隙率的增大，渗透率逐渐增大，惯性系数逐渐减小，有效导热系数逐渐减小。 3.摩擦系数f与雷诺数Re以及Nu数的结论：孔隙率等于1时，在恒热流条件下，f*Re和Nu均为定值；当孔隙率逐渐减小时，f*Re和Nu逐渐增大。7参考文献1 黄善波，刘中良.计算传热学基础.中国石油大学（华东）热能与动力工程系，20092 杨世铭，陶文铨.传热学（第四版）.高等教育出版社，20073 陶文铨关于充分发展的概念及其在传热学教学、科研中的应用J教材

20、通讯， 1989(1):3133 4 陶文铨数值传热学M2版西安：西安交通大学出版社，2001：1245 辛荣昌，陶文铨非稳态导热充分发展阶段的分析解J工程热物理学报，1993，14(1)：8O一836 刘伟,范爱武,黄晓明.多孔介质传热传质理论与应用. 科学出版社, 2006.12 7 张果，刘治军. 长方形截面通道内的充分发展对流换热问题的数值研究。重庆工学院学报(自然科学版) 2007.118 W. Lu a, C.Y. Zhao a,b,*, S.A. Tassou a. Thermal analysis on metal-foam filled heat exchangers. Pa

21、rt I: Metal-foam filled pipes,20058附录附录1 主要程序#include<math.h>#include<stdio.h>#define NT 70void main()/定义变量int N,M,i,j,Iter;/Iter迭代次数 double WNTNT,W0NTNT,wNTNT,sitaNTNT,sita0NTNT,sita_average,sita_average0,De,f,Re;/W无量纲速度数组 W0上一次迭代值的无量纲速度数组 w速度数组/sita无量纲温度数组 sita0上一次迭代值的无量纲温度数组 平均无量纲温度si

22、ta_average/De当量直径 f摩擦系数 Re雷诺数double dltx,dlty,dx,dy,Imda1,Imda2,E,wm,K,D,d,Cf,x,y,Nu,a,b; /E为孔隙率 wm平均速度 K为渗透率 D特征尺寸 d孔隙直径 Cf惯性系数 Nu努塞尔数 a、b为界面尺寸double Eps,DT,DTmax,u,pf;/DT、 DTmax分别为差值和最大差值 u为水的动力粘度 pf水的密度 double ap,ae,aw,an,as,b1,b2;/无量纲速度场、温度场系数，源项 FILE *fp1,*fp2,*fp3,*fp4,*fp5;/*fp1速度场输出结果文件指针 *f

23、p2无量纲速度场输出结果文件指针 *fp3无量纲温度输出结果文件指针/*fp4速度场计算结果tecplot文件指针 *fp5无量纲温度场计算结果tecplot文件指针/读入部分参数printf("请输入x方向节点数N y方向节点数Mn");scanf("%d%d",&N,&M);printf("请输入孔隙率E：n");scanf("%lf",&E);printf("请输入界面长度a 界面宽度bn");scanf("%lf%lf",&a,&

24、;b);printf("请输入流体动力粘度u 流体密度pf 流体平均流速wmn");scanf("%lf%lf%lf",&u,&pf,&wm);/已知参数 dltx=a;/a值大小dlty=b;/b值大小 De=2*dltx*dlty/(dltx+dlty); D=0.04; /特征尺寸大小为a d=1.e-3;/孔隙直径大小 N;/x方向节点数 M;/y方向节点数 Eps=1.e-6;/计算精度 /基本参数计算 K=d*d*E*E*E/150./(1-E)*(1-E); Cf=1.75/sqrt(150.)/pow(E,2./3

25、); dx=dltx/(N-1); dy=dlty/(M-1); ae=aw=dy/dx; an=as=dx/dy; Imda1=E*D*D/K; Imda2=pf*Cf*D*D*D*D*E*E*(-u/K*wm-pf*Cf*wm*wm/sqrt(K)/sqrt(K)/(u*u); /速度边界条件处理/速度与无量纲速度的边界赋0处理for(j=0;j<M;j+) w0j=0; wN-1j=0;W0j=0; WN-1j=0; for(i=1;i<N-1;i+) wi0=0; wiM-1=0;Wi0=0; WiM-1=0; /对速度与无量纲速度初始速度场的赋值for(j=1;j<

26、M-1;j+) for(i=1;i<N-1;i+)Wij=0.0000032388;W0ij=Wij;/对上一次迭代值赋值 Iter=0; /由Gauss-Seidel迭代法计算各点温度 loop1: /内部节点处理for(j=1;j<M-1;j+)for(i=1;i<N-1;i+) b1=dx*dy/(D*D); ap=ae+aw+an+as+(Imda1-Imda2*W0ij)*b1; Wij=ae/ap*Wi+1j+aw/ap*Wi-1j+an/ap*Wij+1+as/ap*Wij-1+b1/ap; wij=Wij*D*D*E*(u/K*wm+pf*Cf*wm*wm/s

27、qrt(K)/u;/计算两次迭代最大误差 Iter=Iter+1; DTmax=0.0; for(j=1;j<M-1;j+) for(i=1;i<N-1;i+) DT=fabs(Wij-W0ij); W0ij=Wij; if(DT>DTmax)DTmax=DT; printf("Iter=%5d DT=%10.6fn",Iter,DTmax);/判断是否满足计算精度 if(DTmax>Eps)goto loop1; /不满足，返回继续迭代/速度场迭代结束，分别输出计算结果/将计算结果输出到文件中 fp1=fopen("result1.dat

28、","w"); printf("速度输出结果：n");/速度输出结果 fprintf(fp1,"速度输出结果：n"); for(j=M-1;j>=0;j-) for(i=0;i<N;i+) printf("%5.6f ",wij); fprintf(fp1,"%5.6f ",wij); printf("n"); fprintf(fp1,"n"); fclose(fp1); /将计算结果输出到文件中 fp2=fopen("re

29、sult2.dat","w"); printf("无量纲速度输出结果：n");/无量纲速度输出结果 fprintf(fp2,"无量纲速度输出结果：n");/无量纲速度输出结果 for(j=M-1;j>=0;j-) for(i=0;i<N;i+) printf("%5.6f ",Wij); fprintf(fp2,"%5.6f ",Wij); printf("n"); fprintf(fp2,"n"); fclose(fp2);/计算摩

30、擦系数f printf("摩擦系数f:n"); f=2*De*(u/K*wm+pf*Cf*wm*wm/sqrt(K)/(pf*wm*wm); printf("%5.6fn",f);/计算雷诺数Re printf("雷诺数Re:n"); Re=pf*wm*De/u; printf("%5.6fn",Re);/无量纲温度边界处理for(j=0;j<M;j+)sita0j=0; sitaN-1j=0; for(i=1;i<N-1;i+)sitai0=0; sitaiM-1=0; /对无量纲温度场的初始赋值fo

31、r(j=1;j<M-1;j+) for(i=1;i<N-1;i+) sitaij=0.2;sita0ij=sitaij; Iter=0; /由Gauss-Seidel迭代法计算各点温度 loop2: /内部节点处理for(j=1;j<M-1;j+) for(i=1;i<N-1;i+) b2=-4*wij*dx*dy/wm/(De*De); ap=ae+aw+an+as; sitaij=ae/ap*sitai+1j+aw/ap*sitai-1j+an/ap*sitaij+1+as/ap*sitaij-1+b2/ap; /计算两次迭代最大误差 Iter=Iter+1; DT

32、max=0.0; for(j=1;j<M-1;j+) for(i=1;i<N-1;i+) DT=fabs(sitaij-sita0ij); sita0ij=sitaij; if(DT>DTmax)DTmax=DT; printf("Iter=%5d DT=%10.6fn",Iter,DTmax);/判断是否满足计算精度 if(DTmax>Eps) goto loop2; /不满足，返回继续迭代/迭代结束，输出计算结果 /将无量纲温度计算结果输出到文件中 fp3=fopen("result3.dat","w");

33、 printf("无量纲温度输出结果：n"); fprintf(fp3,"无量纲温度输出结果：n"); for(j=M-1;j>=0;j-) for(i=0;i<N;i+) printf("%5.6f ",sitaij); fprintf(fp3,"%5.6f ",sitaij); printf("n"); fprintf(fp3,"n"); fclose(fp3);/计算截面流体平均无量纲温度sita_average sita_average=0; for(i=

34、1;i<N-1;i+) for(j=1;j<M-1;j+) sita_average0=wij*sitaij*dx*dy; sita_average=sita_average0+sita_average; sita_average=sita_average/(dltx*dlty*wm); printf("截面流体平均无量纲温度sita_average:n"); printf("%5.6fn",sita_average);/计算努塞尔数Nu Nu=-1./sita_average; printf("努塞尔数Nu:n"); p

35、rintf("%5.6fn",Nu);/将速度场计算结果输出到文件中 fp4=fopen("isotherm1.dat","w"); fprintf(fp4,"VARIABLES="X","Y","SPEED"n"); fprintf(fp4,"ZONE I=%d,J=%d,F=POINTn",N,M); for(j=0;j<M;j+) y=j*dy; for(i=0;i<N;i+) x=i*dx; fprintf(fp4,

36、"%10.5f%10.5f%10.6fn",x,y,wij); fclose(fp4);/将无量纲温度场计算结果输出到文件中 fp5=fopen("isotherm2.dat","w"); fprintf(fp5,"VARIABLES="X","Y","TEMPERTURE"n"); fprintf(fp5,"ZONE I=%d,J=%d,F=POINTn",N,M); for(j=0;j<M;j+) y=j*dy; for(i=

37、0;i<N;i+) x=i*dx; fprintf(fp5,"%10.5f%10.5f%10.6fn",x,y,sitaij); fclose(fp5);附录2验证程序#include<math.h>#include<stdio.h>#define NT 70void main()/定义变量int N,M,i,j,Iter;/Iter迭代次数 double WNTNT,W0NTNT,wNTNT,sitaNTNT,sita0NTNT,sita_average,sita_average0,De,f,Re;/W无量纲速度数组 W0上一次迭代值的无量纲

38、速度数组 w速度数组/sita无量纲温度数组 sita0上一次迭代值的无量纲温度数组 平均无量纲温度sita_average/De当量直径 f摩擦系数 Re雷诺数double dltx,dlty,dx,dy,E=1,wm,K,Cf,D,d,x,y,Nu,a,b,fRe,Pressure_drop; /E为孔隙率 wm平均速度 K为渗透率 D特征尺寸 d孔隙直径 Cf惯性系数 Nu努塞尔数 a、b为界面尺寸double Eps,DT,DTmax,u,pf;/DT、 DTmax分别为差值和最大差值 u为水的动力粘度 pf水的密度 double ap,ae,aw,an,as,b1,b2;/无量纲速度

39、场、温度场系数，源项 FILE *fp1,*fp2,*fp3,*fp4,*fp5;/*fp1速度场输出结果文件指针 *fp2无量纲速度场输出结果文件指针 *fp3无量纲温度输出结果文件指针/*fp4速度场计算结果tecplot文件指针 *fp5无量纲温度场计算结果tecplot文件指针/读入部分参数printf("请输入fRe:n");scanf("%lf",&fRe);printf("请输入界面长度a 界面宽度bn");scanf("%lf%lf",&a,&b);printf("

40、请输入流体动力粘度u 流体密度pf 流体平均流速wmn");scanf("%lf%lf%lf",&u,&pf,&wm);/已知参数 dltx=a;/a值大小dlty=b;/b值大小 De=2*dltx*dlty/(dltx+dlty); D=0.04; /特征尺寸大小为a d=1.e-3;/孔隙直径大小N=51;/x方向节点数M=26;/y方向节点数 Eps=1.e-6;/计算精度 Pressure_drop=-fRe*wm*u/(2*De*De); /基本参数计算 K=d*d*E*E*E/150./(1-E)*(1-E); Cf=1.75

41、/sqrt(150.)/pow(E,2./3); dx=dltx/(N-1); dy=dlty/(M-1); ae=aw=dy/dx; an=as=dx/dy; /速度边界条件处理/速度与无量纲速度的边界赋0处理for(j=0;j<M;j+) w0j=0; wN-1j=0;W0j=0; WN-1j=0; for(i=1;i<N-1;i+) wi0=0; wiM-1=0;Wi0=0; WiM-1=0; /对速度与无量纲速度初始速度场的赋值for(j=1;j<M-1;j+) for(i=1;i<N-1;i+)Wij=0.0000032388;W0ij=Wij;/对上一次迭代

42、值赋值 Iter=0; /由Gauss-Seidel迭代法计算各点温度 loop1: /内部节点处理for(j=1;j<M-1;j+)for(i=1;i<N-1;i+) b1=dx*dy/(D*D); ap=ae+aw+an+as; Wij=ae/ap*Wi+1j+aw/ap*Wi-1j+an/ap*Wij+1+as/ap*Wij-1+b1/ap; wij=-Wij*D*D*Pressure_drop/u;/计算两次迭代最大误差 Iter=Iter+1; DTmax=0.0; for(j=1;j<M-1;j+) for(i=1;i<N-1;i+) DT=fabs(Wij

43、-W0ij); W0ij=Wij; if(DT>DTmax)DTmax=DT; printf("Iter=%5d DT=%10.6fn",Iter,DTmax);/判断是否满足计算精度 if(DTmax>Eps)goto loop1; /不满足，返回继续迭代/速度场迭代结束，分别输出计算结果/将计算结果输出到文件中 fp1=fopen("result1.dat","w"); printf("速度输出结果：n");/速度输出结果 fprintf(fp1,"速度输出结果：n"); for

44、(j=M-1;j>=0;j-) for(i=0;i<N;i+) printf("%5.6f ",wij); fprintf(fp1,"%5.6f ",wij); printf("n"); fprintf(fp1,"n"); fclose(fp1); /将计算结果输出到文件中 fp2=fopen("result2.dat","w"); printf("无量纲速度输出结果：n");/无量纲速度输出结果 fprintf(fp2,"无量纲速度

45、输出结果：n");/无量纲速度输出结果 for(j=M-1;j>=0;j-) for(i=0;i<N;i+) printf("%5.6f ",Wij); fprintf(fp2,"%5.6f ",Wij); printf("n"); fprintf(fp2,"n"); fclose(fp2);/计算摩擦系数f printf("摩擦系数f:n"); f=-2.*De*Pressure_drop/(pf*wm*wm); printf("%5.6fn",f);/计算雷诺数