2005simpler程序讲解陶老师解说nht第12

1、1/69 第12章 求解二维椭圆型流动与换热问题的通用程序应用举例(续）Presented by Professor W Q Tao2005-12-21, Xian2/69目 录127 例题8 旋转圆盘上的冲击流动128 块修正技术129 例题10 带中心射流的通道内的紊流换热 121 例题1 直角坐标中二维稳态无内热源的导热122 例题2 空心圆柱内的稳态热传导123 例题3 正方形管道内的充分发展对流换热124 例题4 内壁上有直肋的环形通道内充分发展对流换热125 例题6 二维突扩通道中的流动与换热126 例题7 方形通道内的复杂充分流动3/69 12-7 旋转圆盘上的冲击流动 圆柱坐标

2、系中动量方程源项的离散已知：一股流体冲击到一个旋转圆盘上，=100,Uin100， 1. 流体从右上侧流出腔体(图8）。求：腔体内的速度及压力分布 这是三维圆柱坐标系中的流动问题，根据所研究问题的特点流动系由于转盘与流体间的粘滞力而引起，三个速度分量的控制方程如下： 12-7-1 物理问题及其数学描写4/69图80.50.5Uin5/69源项源项无压力梯度，因流动系由旋转所致。6/69注意： 方向没有压力梯度，流动系由于切应力而造成。12-7-2 数值方法及程序说明（1）有三个速度分量，但是所有的 ，正如上例中所有的z 方向的导数均为零；（2）在vz,vr 的对流项中没有 ,但 出现在vr的源

3、项中，因而可以将 作为温度那样的 变量处理，与vr,vz 耦合求解，这样实际上用MODE2。 （3）对 以 作为求解变量，其控制方程为：7/69（4）源项的处理8/69迎风模式引入负的 以增加迭代计算的稳定性。9/69*- EXAMPLE EIGHT -* ENTRY GRID TITLE(1)=.VEL U. TITLE(2)=.VEL V. TITLE(3)=.STR FN. TITLE(4)=.R.VTH. TITLE(11)=PRESSURE RELAX(1)=0.8 RELAX(2)=0.8 LSOLVE(1)=.TRUE. LSOLVE(4)=.TRUE. LPRINT(1)=.T

4、RUE. LPRINT(2)=.TRUE. LPRINT(3)=.TRUE. LPRINT(4)=.TRUE. LAST=25 MODE=2 XL=0.5 YL=0.5 L1=7 M1=7 R(1)=0. CALL UGRID RETURN 将第 4个变量赋给R.VTheta10/69 ENTRY START OMEGA=100. UIN=100. DO 100 J=1,M1 DO 100 I=1,L1 U(I,J)=0. V(I,J)=0. F(I,J,4)=0. F(L1,J,4)=R(J)*2*OMEGA ！盘上线速度，引起流体运动的根本原因 100 CONTINUE U(2,2)=UI

5、N AMU=1. RETURN* ENTRY DENSE RETURN 第 4个变量赋是R.VTheta11/69 ENTRY BOUND IF(ITER.NE.0) GOTO 300 FLOWIN=RHO(1,2)*U(2,2)*YCVR(2) 300 FL=0. AR=0. DO 301 I=L3,L2 FLT=R(M1)*XCV(I)*RHO(I,M1) AR=AR+FLT 301 FL=FL+FLT*V(I,M2) ADD=(FLOWIN-FL)/AR DO 302 I=L3,L2 302 V(I,M1)=V(I,M2)+ADD RETURN 通过加一个常数获得满足总体质量守恒条件的出

6、口流速12/69 ENTRY OUTPUT IF(ITER.NE.0) GOTO 400 WRITE(8,401) 401 FORMAT(1X, ITER,7X,SMAX,11X,SSUM,10X,U(4,4), & 9X,V(4,4) 400 WRITE(8,403) ITER,SMAX,SSUM,U(4,4),V(4,4) 403 FORMAT(1X,I6,1P5E15.4) IF(ITER.EQ.LAST) CALL PRINT RETURN13/69 ENTRY GAMSOR IF(ITER.NE.0) GOTO 501 DO 500 J=1,M1 DO 500 I=1,L1 GAM(

7、I,J)=AMU 500 CONTINUE GAM(L3,M1)=0. ！出口局部单向化 GAM(L2,M1)=0. ！出口局部单向化 501 CONTINUE IF(NF.NE.2) GOTO 510 DO 502 J=3,M2 DO 502 I=2,L2 RSWM=FY(J)*F(I,J,4)+FYM(J)*F(I,J-1,4) RHOM=FY(J)*RHO(I,J)+FYM(J)*RHO(I,J-1) CON(I,J)=RHOM*RSWM*2/RMN(J)*3 ！计算Vr的源项 502 AP(I,J)=-AMU/RMN(J)*2！从界面插值得到14/69 510 IF(NF.NE.4)

8、RETURN ！在主网格上对 计算源项 DO 512 J=2,M2 DO 512 I=2,L2 AR=2.*AMU/YCVR(J) CON(I,J)=AR*F(I,J-1,4) 512 AP(I,J)=-AR RETURN END15/6912-7-3 数值计算结果分析 COMPUTATION FOR AXISYMMETRICAL SITUATION * ITER SMAX SSUM U(4,4) V(4,4) 0 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 1 3.1852E-01 0.0000E+00 3.3742E+00 4.8158E+00

9、2 3.6224E-01 1.1921E-07 2.9314E+00 7.6065E+00 3 1.1265E-01 7.4506E-09 1.8755E+00 8.5863E+00 4 6.1974E-02 -3.7253E-08 1.5199E+00 8.8029E+00 5 3.2279E-02 -3.1665E-08 1.2971E+00 8.4019E+00 6 1.7869E-02 -4.0280E-08 1.2738E+00 7.6836E+00 7 1.2370E-02 5.1223E-09 1.3363E+00 6.8852E+00 8 1.0312E-02 -1.1176E

10、-08 1.4400E+00 6.1421E+00 9 7.9294E-03 -2.9569E-08 1.5480E+00 5.5244E+00 10 5.9429E-03 4.8894E-08 1.6437E+00 5.0452E+00 11 4.6140E-03 -1.6531E-08 1.7207E+00 4.6926E+00 12 3.3741E-03 3.1199E-08 1.7787E+00 4.4432E+00 13 2.6291E-03 -5.1106E-08 1.8202E+00 4.2728E+00 16/69 14 1.9695E-03 -2.6543E-08 1.848

11、6E+00 4.1597E+00 15 1.4364E-03 6.2981E-08 1.8674E+00 4.0867E+00 16 1.0142E-03 -4.5111E-08 1.8792E+00 4.0409E+00 17 6.9815E-04 8.9640E-09 1.8864E+00 4.0129E+00 18 4.6667E-04 3.8388E-08 1.8906E+00 3.9963E+00 19 3.0389E-04 3.3469E-09 1.8929E+00 3.9868E+00 20 1.9290E-04 -1.1176E-08 1.8941E+00 3.9816E+00

12、 21 1.1830E-04 5.2169E-09 1.8946E+00 3.9790E+00 22 7.0846E-05 4.6941E-08 1.8947E+00 3.9778E+00 23 4.0823E-05 5.4388E-08 1.8947E+00 3.9773E+00 24 2.2590E-05 -8.0094E-08 1.8945E+00 3.9772E+00 25 1.1003E-05 -3.8743E-08 1.8944E+00 3.9773E+0017/69 * .VEL U. * I = 2 3 4 5 6 7 J 7 0.00E+00 0.00E+00 0.00E+0

13、0 0.00E+00 0.00E+00 0.00E+00 6 0.00E+00-1.33E+00-2.67E+00-2.12E+00-8.37E-01 0.00E+00 5 0.00E+00-1.86E+00-2.70E+00-1.86E+00-6.39E-01 0.00E+00 4 0.00E+00-2.17E-01 1.89E+00 2.90E+00 1.65E+00 0.00E+00 3 0.00E+00 1.33E+01 1.97E+01 1.92E+01 1.04E+01 0.00E+00 2 1.00E+02 8.63E+01 7.43E+01 5.99E+01 3.27E+01

14、0.00E+00 1 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 18/69 * .VEL V. * I = 1 2 3 4 5 6 7 J 7 0.00E+00 0.00E+00 0.00E+00 0.00E+00 3.99E+00 6.01E+00 0.00E+00 6 0.00E+00-1.50E+00-1.50E+00 6.18E-01 6.44E+00 8.45E+00 0.00E+00 5 0.00E+00-4.17E+00-2.98E+00 1.81E+00 1.00E+01 1.20E+01 0.00E+00 4

15、0.00E+00-6.53E+00-1.84E+00 3.98E+00 1.34E+01 1.60E+01 0.00E+00 3 0.00E+00 6.87E+00 5.96E+00 7.21E+00 1.36E+01 1.63E+01 0.00E+00 2 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 19/69 * .STR FN. * I = 2 3 4 5 6 7 J 7 5.00E-01 5.00E-01 5.00E-01 5.00E-01 3.00E-01 0.00E+00 6 5.00E-01 5.6

16、0E-01 6.20E-01 5.95E-01 3.38E-01 0.00E+00 5 5.00E-01 6.25E-01 7.15E-01 6.60E-01 3.60E-01 0.00E+00 4 5.00E-01 6.31E-01 6.67E-01 5.88E-01 3.19E-01 0.00E+00 3 5.00E-01 4.31E-01 3.72E-01 3.00E-01 1.63E-01 0.00E+00 2 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 壳体流量为常数转盘上流量为零20/69 * R. VTH * I =

17、 1 2 3 4 5 6 7 J 7 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 2.50E+01 6 0.00E+00 1.24E-01 5.24E-01 1.64E+00 5.76E+00 1.26E+01 2.02E+01 5 0.00E+00 2.02E-01 7.28E-01 1.69E+00 3.66E+00 7.75E+00 1.23E+01 4 0.00E+00 1.40E-01 4.46E-01 8.49E-01 1.53E+00 3.54E+00 6.25E+00 3 0.00E+00 5.15E-02 1.4

18、9E-01 2.47E-01 3.84E-01 1.09E+00 2.25E+00 2 0.00E+00 4.66E-03 1.84E-02 3.72E-02 5.53E-02 1.55E-01 2.50E-01 1 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 21/69 * PRESSURE * I = 1 2 3 4 5 6 7 J 7 -4.93E+02-4.81E+02-4.57E+02-3.68E+02-3.47E+02-3.61E+02-3.61E+02 6 -5.08E+02-4.96E+02-4.

19、72E+02-3.94E+02-3.61E+02-3.61E+02-3.61E+02 5 -5.38E+02-5.26E+02-5.02E+02-4.46E+02-3.89E+02-3.61E+02-3.47E+02 4 -6.85E+02-6.47E+02-5.72E+02-4.92E+02-3.60E+02-2.41E+02-1.81E+02 3 -1.15E+03-9.63E+02-5.97E+02-4.57E+02-1.85E+02 1.02E+02 2.46E+02 2 -3.01E+02-3.62E+02-4.84E+02-3.04E+02 1.83E+02 6.20E+02 8.

20、39E+02 1 0.00E+00-6.11E+01-4.27E+02-2.28E+02 3.67E+02 8.79E+02 1.10E+03 22/69 应用一般ADI方法求解代数方程时，有时收敛速度减慢的例子： 计算区域中出现第二类边界条件（尤其是绝热条件时）迭代收敛速度减慢。如长方形截面通道中的充分发展对流，如果以整个计算区域计算，迭代收敛较快，因为四个边界均为第一类边界条件，但所需节点数较多；改为1/4区域以后，节点数减少了，但迭代收敛速度放慢。 12-8 块修正技术12-8-1 块修正的需要23/6924/69 更为严重的是如右下图所示情形：从网格划分所形成的离散方程来看， ，即南北

21、节点影响系数大；但从影响的传递角度来讲，应尽量使 传入到区域内部。由于这一矛盾的因素，使这类问题求解迭代收敛速度较慢。因而需要一种能加速迭代收敛的方法。 25/6912-8-2 块修正基本思想 对二维问题，我们希望设计出一种辅助算法，这种算法能在平均意义上把边界条件的影响很快传入到区域内部，但又只须求解一个一维的代数方程组。 2. 从物理意思上来说，所谓迭代收敛的过程就是所得计算结果 逐渐满足守恒性的过程。因而所谓迭代计算一次也就是在原来不满足守恒性的解 加上一个修正值， 使 能满足守恒要求。26/69 为了只求解一个三对角阵的方程，我们不是对 每一个点的来加上各自的修正值，而是对计算区域 的

22、一个块加上一个修正值， 或 ， 使该块中修正后的值， ， 可以在整体上满足守恒性，这种方法就称为块修 正法。下面我们来讨论修正 值， 或 方程 应当怎样 来形成。27/69 以下图所示的竖直块为例来推导。12-8-3 单块修正值的代数方程系数及其边界条件 28/69 （1）设上一次迭代获得 ，第i块的修正值为 ， 则修正后的 应在块的总体上满足守恒要求，于是有： 29/69将此式展开合并写成 的代数方程后有： 其中系数为：北边界修正值为零南边界修正值为零30/69左右邻点影响边界上的 应该包括在内31/69对 其系数可仿此写出。 （2）边界条件 在本算法中认为所有变量都给出为第一类边界条件，亦

23、即 值给定的边界条件，因而边界中的修正值应为0, 无论 还是 其边界值均为零：32/6933/6912-8-4 双块、多块修正 双块修正：计算区域中有两个物性截然不同的区域时单块修正可能不适用，需要双块修正 多块修正：与多重网格类似。 34/6912-9 带中心射流的通道内的紊流换热 湍流模型及壁面函数法 已知：如图10所示，带一般中心射流的流体进入一 平行板通道。流动为紊流， 106 , Pr=0.7。 求： 采用 标准 模型及 壁面函数法， 确定 突扩区域中的流场、温度场。12-9-1 物理问题及其数学描写35/69T=400V=100 Tw=01.04.0T=100V=10图1036/6

24、9控制方程为：其中各个变量的扩散系数等如下：NF= 1 2 3 4 5 6 7 8 11 变量 U V PC T k AMUT ,GEN P 0.80.80.60.60.637/69边界条件为：（1）入口处：k 取来流平均动能的1； 按以下公式计算其中（2）壁面处按壁面函数法；（3）出口处按局部单向化处理38/69（4）对称线上，u=0, 其余法向一阶导数为零。12-9-2 数值方法及程序说明（1） 源项的处理39/69 采用局部单向化条件时，为防止迭代过程中出口法向流速出现负值对计算区域内的影响，采用提升方法避免之：outflowVmin（2）出口速度的提升40/69（3）动量离散方程中源项

25、的处理 DIMENSION T(NI,NJ),XPLUS(NI),AKE(NI,NJ),DIS(NI,NJ) EQUIVALENCE (F(1,1,4),T(1,1),(F(1,1,5),AKE(1,1),(F(1,1,6),DIS(1, &1) DIMENSION AMUT(NI,NJ),GEN(NI,NJ) EQUIVALENCE (F(1,1,7),AMUT(1,1),(F(1,1,8),GEN(1,1)* ENTRY GRID TITLE(1)=.VEL U. TITLE(2)=.VEL V. TITLE(3)=.STR FN. TITLE(4)=. TEMP . TITLE(5)=K

26、IN ENE TITLE(6)=.DISIPA. TITLE(7)=TURB VI TITLE(11)=PRESSURE TITLE(12)=DENSITY *41/69RELAX(1)=0.8RELAX(2)=0.8RELAX(5)=0.6RELAX(6)=0.6RELAX(13)=0.6LSOLVE(1)=.TRUE.LSOLVE(5)=.TRUE.LSOLVE(6)=.TRUE.LPRINT(1)=.TRUE.LPRINT(2)=.TRUE.LPRINT(3)=.TRUE.LPRINT(4)=.TRUE.LPRINT(5)=.TRUE.LPRINT(6)=.TRUE.LPRINT(7)=

27、.TRUE.LPRINT(11)=.TRUE.LAST=50XL=1.YL=4.L1=7M1=9CALL UGRIDRETURN 42/69 ENTRY START DO 100 J=1,M1 DO 100 I=1,L1 U(I,J)=0. V(I,J)=10. V(1,J)=0. V(I,2)=10. IF(I.GT.4) V(I,2)=100. T(I,J)=100. T(1,J)=0. IF(I.GT.4) T(I,1)=400. AKE(I,J)=0.005*V(I,2)*2 ！进口动能的1 DIS(I,J)=0.1*AKE(I,J)*2 100 CONTINUE 43/69 AMU=

28、1.E-6 CMU=0.09 C1=1.44 C2=1.92 PRT=0.9 PRK=1.0 PRD=1.3 PR=0.7 PRPRT=PR/PRT PFN=9.*(PRPRT-1.)/PRPRT*.25 温度壁面函数法中的P函数 CMU4=CMU*.25 RETURN ENTRY DENSE RETURN 44/69 ENTRY BOUND IF(ITER.NE.0) GOTO 300 FLOWIN=0. DO 310 I=2,L2 310 FLOWIN=FLOWIN+RHO(I,1)*V(I,2)*XCV(I) ！入口流量 300 FL=0. AFL=0. VMIN=0. DO 301 I

29、=2,L2 IF(V(I,M2).LT.0.) VMIN=AMAX1(VMIN,-V(I,M2) ！搜索Vmin AFL=AFL+RHO(I,M1)*XCV(I) 301 FL=FL+RHO(I,M1)*V(I,M2)*XCV(I) FACTOR=FLOWIN/(FL+AFL*VMIN) DO 302 I=2,L2 V(I,M1)=(V(I,M2)+VMIN)*FACTOR 302 CONTINUE DO 303 J=2,M2 AKE(L1,J)=AKE(L2,J) ！相当于充分发展 303 DIS(L1,J)=DIS(L2,J) RETURN 45/69 ENTRY OUTPUT IF(IT

30、ER.NE.0) GOTO 400 WRITE(8,401) 401 FORMAT(1X, ITER,6X,SMAX,6X,SSUM,5X,V(6,6), & 4X,T(5,6),4X,KE(5,6) 400 WRITE(8,403) ITER,SMAX,SSUM,V(6,6),T(5,6),AKE(5,6) 403 FORMAT(1X,I6,1P5E11.3) IF(ITER.EQ.LAST) CALL PRINT RETURN46/69 ENTRY GAMSOR IF(NF.EQ.3) RETURN IF(NF.NE.1) GOTO 510 REL=1.-RELAX(NGAM) DO 50

31、0 J=1,M1 DO 500 I=1,L1 AMT=CMU*RHO(I,J)*AKE(I,J)*2/(DIS(I,J)+1.E-30)！初始值 IF(ITER.EQ.0) AMUT(I,J)=AMT 500 AMUT(I,J)=RELAX(NGAM)*AMT+REL*AMUT(I,J) ！亚松弛 510 FACTOR=1. IF(NF.EQ.4) FACTOR=1./PRT IF(NF.EQ.5) FACTOR=1./PRK IF(NF.EQ.6) FACTOR=1./PRD DO 520 J=1,M1 DO 520 I=1,L1 GAM(I,J)=AMUT(I,J)*FACTOR IF(N

32、F.NE.1) GAM(L1,J)=0. ！对称线 520 GAM(I,M1)=0. ！出口局部单向化 47/69 DO 530 J=2,M2 GOTO(531,532,531,533,531,531),NF 531 GAM(1,J)=0. GOTO 530 532 XPLUS(J)=RHO(2,J)*SQRT(AKE(2,J)*CMU4*XDIF(2)/AMU GAM(1,J)=AMU IF(XPLUS(J).GT.11.5) GAM(1,J)=AMU*XPLUS(J)/(ALOG(9.*XPLUS(J)*2.5) GOTO 530 533 GAM(1,J)=AMU/PR IF(XPLUS(

33、J).GT.11.5) GAM(1,J)=AMU/PRT*XPLUS(J) & /(2.5*ALOG(9.*XPLUS(J)+PFN) 530 CONTINUE 实施壁面函数法 耗散率为对第一个内节点赋值， 因此也可以切断与边界的联系48/6949/6950/6951/69 IF(NF.NE.1) GOTO 509 DO 590 J=2,M2 DO 590 I=3,L2 CON(I,J)=(GAM(I,J)*(U(I+1,J)-U(I,J)/XCV(I) & -GAM(I-1,J)*(U(I,J)-U(I-1,J)/XCV(I-1)/XDIF(I) GAMP=GAM(I,J+1)*GAM(I-

34、1,J+1)/(GAM(I,J+1)+GAM(I-1,J+1)+1.E-30) GAMP=GAMP+GAM(I,J)*GAM(I-1,J)/(GAM(I,J)+GAM(I-1,J)+1.E-30) GAMM=GAM(I,J-1)*GAM(I-1,J-1)/(GAM(I,J-1)+GAM(I-1,J-1)+1.E-30) GAMM=GAMM+GAM(I,J)*GAM(I-1,J)/(GAM(I,J)+GAM(I-1,J)+1.E-30) CON(I,J)=CON(I,J)+(GAMP*(V(I,J+1)-V(I-1,J+1) & -GAMM*(V(I,J)-V(I-1,J)/(YCV(J)*XD

35、IF(I) AP(I,J)=0. 590 CONTINUE RETURN 计算u方程的源项52/69509 IF(NF.NE.2) GOTO 519 DO 595 J=3,M2 DO 595 I=2,L2 CON(I,J)=(GAM(I,J)*(V(I,J+1)-V(I,J)/YCV(J)- & GAM(I,J-1)*(V(I,J)-V(I,J-1)/YCV(J-1)/(YDIF(J) GAMP=GAM(I+1,J)*GAM(I+1,J-1)/(GAM(I+1,J)+GAM(I+1,J-1)+1.E-30) GAMP=GAMP+GAM(I,J)*GAM(I,J-1)/(GAM(I,J)+GAM

36、(I,J-1)+1.E-30) GAMM=GAM(I-1,J)*GAM(I-1,J-1)/(GAM(I-1,J)+GAM(I-1,J-1)+1.E-30) GAMM=GAMM+GAM(I,J)*GAM(I,J-1)/(GAM(I,J)+GAM(I,J-1)+1.E-30) CON(I,J)=CON(I,J)+(GAMP*(U(I+1,J)-U(I+1,J-1) & -GAMM*(U(I,J)-U(I,J-1)/(XCV(I)*YDIF(J) AP(I,J)=0. 595 CONTINUE RETURN 519 CONTINUE 计算v方程的源项53/69 IF(NF.NE.4) GOTO 52

37、9 DO 596 J=2,M2 DO 596 I=2,L2 CON(I,J)=0. 596 AP(I,J)=0. RETURN 计算T方程的源项54/69 529 IF(NF.NE.5) GOTO 539 DO 598 J=2,M2 DO 598 I=2,L2 DUDX=(U(I+1,J)-U(I,J)/XCV(I) DVDY=(V(I,J+1)-V(I,J)/YCV(J) IF(J.EQ.2) GOTO 581 IF(J.EQ.M2) GOTO 582 DUDY=(0.5*(U(I,J+1)-U(I,J-1)+0.5*(U(I+1,J+1)-U(I+1,J-1) & /(YDIF(J)+YD

38、IF(J+1) GOTO 583 581 DUDY=(0.5*(U(I,J+1)-U(I,J)+0.5*(U(I+1,J+1)-U(I+1,J)/YDIF(J+1) GOTO 583 582 DUDY=(0.5*(U(I,J)-U(I,J-1)+0.5*(U(I+1,J)-U(I+1,J-1) & /YDIF(J) 583 IF(I.EQ.2) GOTO 584 IF(I.EQ.L2) GOTO 585 DVDX=(0.5*(V(I+1,J)-V(I-1,J)+0.5*(V(I+1,J+1) & -V(I-1,J+1)/(XDIF(I)+XDIF(I+1)计算k方程的源项55/69 GOTO

39、586 584 DVDX=(0.5*(V(I+1,J)-V(I,J)+0.5*(V(I+1,J+1) & -V(I,J+1)/XDIF(I+1) GOTO 586 585 DVDX=(0.5*(V(I,J)-V(I-1,J)+0.5*(V(I,J+1) & -V(I-1,J+1)/XDIF(I) 586 GEN(I,J)=2.*(DUDX*2+DVDY*2)+(DUDY+DVDX)*2 CON(I,J)=GEN(I,J)*AMUT(I,J) AP(I,J)=-RHO(I,J)*DIS(I,J)/(AKE(I,J)+1.E-30) 598 CONTINUE RETURN 539 CONTINUE

40、 计算k方程的源项56/69 DO 599 J=2,M2 DO 599 I=2,L2 CON(I,J)=C1*GEN(I,J)*CMU*RHO(I,J)*AKE(I,J) AP(I,J)=-C2*RHO(I,J)*DIS(I,J)/(AKE(I,J)+1.E-30) 599 CONTINUE DO 600 J=2,M2 DISS=CMU*AKE(2,J)*1.5/(0.4*CMU4*XDIF(2) CON(2,J)=1.E30*DISS AP(2,J)=-1.E30 600 CONTINUE RETURN END计算方程的源项57/69 COMPUTATION IN CARTISIAN COO

41、RDINATES * ITER SMAX SSUM V(6,6) T(5,6) KE(5,6) 0 0.000E+00 0.000E+00 1.000E+01 1.000E+02 5.000E+01 1 9.610E+00 1.192E-05 5.691E+01 3.871E+02 5.723E+01 2 9.514E+00 -2.146E-06 5.821E+01 3.925E+02 6.918E+01 3 1.631E+01 -2.623E-06 6.358E+01 3.929E+02 8.412E+01 4 1.852E+01 -3.219E-06 6.502E+01 3.922E+02

42、 1.001E+02 5 7.633E+00 1.669E-06 6.859E+01 3.888E+02 1.158E+02 6 3.936E+00 2.921E-06 7.270E+01 3.844E+02 1.306E+02 7 2.899E+00 8.047E-07 7.596E+01 3.811E+02 1.450E+02 8 2.467E+00 -1.520E-06 7.880E+01 3.781E+02 1.588E+02 9 1.233E+00 3.767E-06 8.126E+01 3.754E+02 1.720E+02 10 5.827E-01 8.019E-07 8.348

43、E+01 3.730E+02 1.848E+02 11 4.670E-01 -3.013E-07 8.535E+01 3.708E+02 1.974E+02 12 2.132E-01 -1.404E-06 8.684E+01 3.690E+02 2.098E+02 13 1.752E-01 2.645E-06 8.799E+01 3.675E+02 2.219E+02 14 2.045E-01 7.786E-07 8.881E+01 3.663E+02 2.334E+02 15 1.997E-01 1.352E-06 8.936E+01 3.654E+02 2.443E+02 58/69 16

44、 1.952E-01 -3.356E-06 8.968E+01 3.647E+02 2.543E+02 17 1.732E-01 -5.886E-07 8.983E+01 3.642E+02 2.633E+02 18 1.515E-01 1.214E-06 8.985E+01 3.639E+02 2.711E+02 19 1.275E-01 2.459E-06 8.978E+01 3.637E+02 2.778E+02 20 1.135E-01 -1.770E-07 8.966E+01 3.636E+02 2.836E+02 21 9.660E-02 1.088E-06 8.950E+01 3

45、.635E+02 2.884E+02 22 8.860E-02 1.376E-06 8.932E+01 3.635E+02 2.924E+02 23 8.655E-02 4.222E-06 8.913E+01 3.635E+02 2.958E+02 24 8.673E-02 1.618E-06 8.894E+01 3.635E+02 2.987E+02 25 8.763E-02 6.519E-08 8.874E+01 3.635E+02 3.011E+02 26 8.823E-02 -1.248E-06 8.855E+01 3.636E+02 3.032E+02 27 8.634E-02 1.

46、515E-06 8.837E+01 3.637E+02 3.050E+02 28 8.221E-02 -9.155E-07 8.820E+01 3.637E+02 3.066E+02 29 7.629E-02 -4.168E-07 8.803E+01 3.638E+02 3.079E+02 30 6.849E-02 4.278E-06 8.789E+01 3.639E+02 3.090E+02 59/69 31 6.000E-02 -1.577E-06 8.776E+01 3.639E+02 3.100E+02 32 5.131E-02 -1.215E-06 8.764E+01 3.640E+

47、02 3.109E+02 33 4.320E-02 1.020E-06 8.753E+01 3.640E+02 3.117E+02 34 3.870E-02 -1.668E-06 8.744E+01 3.640E+02 3.125E+02 35 3.469E-02 -1.627E-06 8.736E+01 3.641E+02 3.132E+02 36 3.132E-02 2.183E-06 8.728E+01 3.641E+02 3.138E+02 37 2.813E-02 -1.673E-06 8.722E+01 3.641E+02 3.145E+02 38 2.516E-02 -2.713

48、E-06 8.715E+01 3.641E+02 3.151E+02 39 2.318E-02 7.274E-07 8.710E+01 3.641E+02 3.157E+02 40 2.092E-02 5.514E-06 8.705E+01 3.642E+02 3.163E+02 41 1.954E-02 -5.197E-07 8.700E+01 3.642E+02 3.169E+02 42 1.805E-02 -7.967E-07 8.695E+01 3.642E+02 3.174E+02 43 1.683E-02 3.801E-06 8.691E+01 3.642E+02 3.179E+0

49、2 44 1.575E-02 -3.199E-06 8.687E+01 3.642E+02 3.184E+02 45 1.476E-02 2.365E-06 8.684E+01 3.642E+02 3.188E+02 46 1.418E-02 2.495E-06 8.680E+01 3.642E+02 3.192E+02 47 1.367E-02 4.471E-06 8.678E+01 3.643E+02 3.196E+02 48 1.321E-02 5.106E-07 8.675E+01 3.643E+02 3.199E+02 49 1.282E-02 1.093E-06 8.672E+01

50、 3.643E+02 3.202E+02 50 1.222E-02 -4.708E-07 8.670E+01 3.643E+02 3.205E+02 60/69* .VEL U. * I = 2 3 4 5 6 7 J 9 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 8 0.00E+00 1.56E-02 3.75E-02 3.84E-02 2.04E-02 0.00E+00 7 0.00E+00-1.65E+00-2.68E+00-2.78E+00-1.33E+00 0.00E+00 6 0.00E+00-2.37E+00-3.

51、56E+00-3.57E+00-1.63E+00 0.00E+00 5 0.00E+00-2.38E+00-3.88E+00-3.98E+00-1.66E+00 0.00E+00 4 0.00E+00-1.39E+00-3.33E+00-3.86E+00-1.45E+00 0.00E+00 3 0.00E+00 3.74E+00-3.47E-01 -2.75E+00-8.62E-01 0.00E+00 2 0.00E+00 4.44E+00 6.55E+0 0-2.87E+00-6.77E-01 0.00E+00 1 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.

52、00E+00 0.00E+0061/69 * .VEL V. * I = 1 2 3 4 5 6 7 J 9 0.00E+00 8.87E+00 3.18E+01 4.59E+01 6.52E+01 7.82E+01 1.00E+01 8 0.00E+00 8.87E+00 3.18E+01 4.59E+01 6.52E+01 7.82E+01 1.00E+01 7 0.00E+00 4.16E+00 2.89E+01 4.56E+01 6.93E+01 8.20E+01 1.00E+01 6 0.00E+00-2.61E+00 2.55E+01 4.56E+01 7.48E+01 8.67E

53、+01 1.00E+01 5 0.00E+00-9.41E+00 2.12E+01 4.53E+01 8.15E+01 9.14E+01 1.00E+01 4 0.00E+00-1.34E+01 1.56E+01 4.38E+01 8.83E+01 9.56E+01 1.00E+01 3 0.00E+00-2.70E+00 3.98E+00 3.69E+01 9.37E+01 9.81E+01 1.00E+01 2 1.00E+01 1.00E+01 1.00E+01 1.00E+01 1.00E+02 1.00E+02 1.00E+02 62/69 * .STR FN * I = 2 3 4

54、 5 6 7 J 9 0.00E+00-1.77E+00 -8.12E+00-1.73E+01-3.03E+01-4.60E+01 8 0.00E+00-1.77E+00 -8.14E+00-1.73E+01-3.04E+01-4.60E+01 7 0.00E+00-8.31E-01 -6.61E+00-1.57E+01-2.96E+01-4.60E+01 6 0.00E+00 5.21E-01 -4.58E+00-1.37E+01 -2.87E+01-4.60E+01 5 0.00E+00 1.88E+00 -2.36E+0 0-1.14E+01-2.77E+01-4.60E+01 4 0.

55、00E+00 2.68E+00 -4.55E-01 -9.22E+00 -2.69E+01-4.60E+01 3 0.00E+00 5.39E-01 -2.57E-01 -7.64E+00 -2.64E+01-4.60E+01 2 0.00E+00-2.00E+00 -4.00E+00-6.00E+00 -2.60E+01-4.60E+01 63/69 * . TEMP . * I = 1 2 3 4 5 6 7 J 9 0.00E+00 1.00E+02 1.00E+02 1.00E+02 1.00E+02 1.00E+02 1.00E+02 8 0.00E+00 3.01E+02 3.26

56、E+02 3.39E+02 3.60E+02 3.80E+02 1.00E+02 7 0.00E+00 3.00E+02 3.21E+02 3.35E+02 3.63E+02 3.85E+02 1.00E+02 6 0.00E+00 2.93E+02 3.10E+02 3.26E+02 3.64E+02 3.89E+02 1.00E+02 5 0.00E+00 2.80E+02 2.93E+02 3.11E+02 3.67E+02 3.92E+02 1.00E+02 4 0.00E+00 2.65E+02 2.69E+02 2.88E+02 3.72E+02 3.95E+02 1.00E+02

57、 3 0.00E+00 2.52E+02 2.36E+02 2.53E+02 3.79E+02 3.97E+02 1.00E+02 2 0.00E+00 1.29E+02 1.16E+02 2.01E+02 3.90E+02 3.99E+02 1.00E+02 1 0.00E+00 1.00E+02 1.00E+02 1.00E+02 4.00E+02 4.00E+02 4.00E+02 64/69 * KIN ENE * I = 1 2 3 4 5 6 7 J 9 5.00E-01 5.00E-01 5.00E-01 5.00E-01 5.00E+01 5.00E+01 5.00E+01 8

58、 5.00E-01 1.59E+02 4.93E+02 4.65E+02 3.53E+02 2.15E+02 2.15E+02 7 5.00E-01 1.90E+02 5.34E+02 4.85E+02 3.35E+02 1.74E+02 1.74E+02 6 5.00E-01 2.20E+02 5.83E+02 5.22E+02 3.20E+02 1.37E+02 1.37E+02 5 5.00E-01 2.39E+02 6.06E+02 5.46E+02 2.94E+02 1.06E+02 1.06E+02 4 5.00E-01 2.15E+02 5.40E+02 5.31E+02 2.5

59、4E+02 8.23E+01 8.23E+01 3 5.00E-01 1.15E+02 3.30E+02 4.69E+02 2.06E+02 6.62E+01 6.62E+01 2 5.00E-01 1.88E+01 1.03E+01 3.22E+02 1.46E+02 5.55E+01 5.55E+01 1 5.00E-01 5.00E-01 5.00E-01 5.00E-01 5.00E+01 5.00E+01 5.00E+01 65/69 * .DISIPA. * I = 1 2 3 4 5 6 7 J 9 2.50E-02 2.50E-02 2.50E-02 2.50E-02 2.50E+0