《内部存在钝体的平板拖曳流数值模拟分析》

1、内部存在钝体的平板拖曳流数值模拟分析一、问题介绍如图有一矩形区域，宽W=2.5m,高H=1m。在矩形区域内，有一钝体，钝体位于下边界中央，宽0.1m，高0.5m。区域左壁面温度0K，下壁面50K；上、右壁面为100K。以水为介质，动力粘度0.001N·m/s2，密度1000kg/m3。我们分别计算了顶盖以0.001m/s，0.002m/s，0.005m/s，0.01m/四种速度向右运动时的流场，对应的雷诺数分别为1000，2000，5000，10000。 二、SIMPLE算法程序简介 本次数值模拟使用的是Fortran编写的SIMPLE算法程序，SIMPLE算法是一种经典的基于压力的

2、半隐式数值模拟方法，它的计算过程大致可分为以下几步：（1）假定一速度分布，记为,以此计算动量离散方程中的系数及常数项；（2）假设一个压力场；（3）依次求解动量方程，得；（4）求解压力校正方程，得；（5）根据校正压力改进速度值；（6）利用改进后的速度场求解那些通过源项物性等与速度场耦合的变量，如果 变量并不影响流场，则应在速度场收敛后再求解；（7）利用改进后的速度场重新计算动量离散方程的系数，并利用改进后的压力场作为下一层次迭代计算的初值。重复上述步骤，直到获得收敛的解。三、程序设计简要介绍（一）主程序模块LOGICAL LSTOP ! 逻辑变量 LSTOP，为真时计算停止COMMON/CNTL

3、/LSTOPCALL USER(1) ! 绘制网格CALL SETUP1 ! 设置交错网格CALL USER(2) ! 设置初值，包括松弛因子、迭代步数、物性参数等CALL USER(3) ! 密度的计算（可压缩流体密度计算时需要，本算例设置为常量）CALL USER(4) ! 设定边界条件CALL USER(5) ! 数据输出IF(LSTOP)CALL PRINT ! 调用打印模块IF(LSTOP) STOP ! 停止程序CALL SETUP2 ! 主要计算核心模块GO TO 10END（二）USER(1)模块（绘制网格）MODE=1 ! 坐标为直角坐标系C* USER DATA: B X

4、H X W (UNIT: M)H=1W=2.5U_REFF=1.E-8! 设置内部速度为1.E-8，为了在接下来的计算中使分母不为0C*L1=250 ! 设置 X 向网格数M1=100 ! 设置 Y 向网格数XL=WYL=HCALL SUPPLY(1) ! 调用 SUPPLY 函数，确定网格点坐标RETURN（三）USER(2)模块（设置初值）RELAX(1)=0.45RELAX(2)=0.45RELAX(11)=0.8RELAX(13)=0.5 ! 设置各参数对应的松弛因子LAST=5000 ! 设置迭代步数为 5000AMU=0.001 ! 设置动力粘度为 0.001(PA·S)

5、AMU_TEMP=1. ! 设置导热系数 K为1DO 100 J=2,M2DO 100 I=2,L2V(I,J)=0.0U(I,J)=U_REFFRHO(I,J)=1000. ! 设置密度为1000kg/m3GAM(I,J)=AMUT(I,J)=0. ! 对各变量赋初值100 CONTINUEDO 1000 I=120,130 DO 1000 J=2,50 ! 确定钝体位置V(I,J)=0.0U(I,J)=0.0GAM(I,J)=AMU*2.E6 ! 将钝体内部网格点的动力粘度设置为无穷大1000 CONTINUE（四）USER(4)模块（设置边界条件）DO 301 I=2,L2DO 301

6、J=2,M2 U(I,M1)=1.E-3！上表面x方向速度为0.001m/sV(I,M1)=0.U(I,1)=0. ! 下表面固定V(I,2)=0.U(2,J)=0.V(1,J)=0. ! 左表面固定U(L1,J)=0.V(L1,J)=0 ! 右表面固定301 CONTINUEDO 302 I=1,L1DO 302 J=1,M1T(I,1)=50 ! 下边界温度T(1,J)=0! 左边界温度T(I,M1)=100 ! 上边界温度T(L1,J)=100 ! 右边界温度302 CONTINUE（五）USER(6)模块（设置扩散系数）IF(NF.EQ.1 .OR. NF.EQ.2) THEN GAM

7、1=AMU ! 计算U 、V时调用扩散系数DO 520 I=1,L1IF(NF.NE.2)GAM(I,M1)=0.0IF(NF.NE.2)GAM(I,1) =0.0 ! 设置边界扩散系数，定义为 0520 CONTINUEEND IFIF(NF.EQ.4)THENGAM1=AMU_TEMP ! 计算 U 和 V 时调用扩散系数为 kEND IFDO 521 I=1,L1DO 521 J=1,M1IF(I>119.AND.I<131.AND.J>1.AND.J<51)THENGAM(I,J)=GAM1*2.E5 ! 钝体内部扩散系数为无穷大ELSEGAM(I,J)=GAM

8、1END IF521CONTINUE四、计算结果及分析（一）Re=10000图4.1(a) Re=10000 流线图以Re=10000作典型分析，从如上流线图发现，一共出现9个漩涡。分别是：在钝体右边出现了一个顺时针漩涡，同时出现与这个大涡旋转方向相反的两个角涡；在钝体左边出现了一个逆时针漩涡，同时出现与大涡旋转方向相反的两个角涡；再顶盖中部附近产生了一个顺时针小漩涡；在钝体上方产生了两个逆时针小漩涡。一个有趣的现象是钝体把整个流场分割成了两部分，而这两部分有一定的对称性，分别被一个尺寸较大的涡占据，而在两个大涡的底面又各有两个角涡，但是这两部分中的涡转动方向却全都是相反的。我们这样解释涡的形

9、成机理：顶盖向右拖拽后，使钝体右边流体产生了顺时针的角动量，这部分动量一部分产生了顺时针的大漩涡，一部分流向左边与顶盖的驱动共同使顶盖中部附近产生了一个顺时针漩涡及在钝体上方产生了两个逆时针小漩涡。钝体左边大漩涡的产生，一部分是源于右边大漩涡的能量，还有一部分是顶盖中部漩涡的一部分能量分离所得。四个角涡的出现则是因为大漩涡对拐角壁面附近的碰撞及离心力对流体的作用，形成了二次回流。图4.1(b) Re=10000 速度场云图从如上速度场云图发现，钝体右边大漩涡的流体速度及顶盖中间附近的流体速度明显很大，钝体左边大漩涡的流速及钝体上面的两个小漩涡流速很小，这在一定程度上可以验证对图4.1(a)流线

10、图形成的解释。图4.1(c) Re=10000 压力场云图压力场云图中可以观察到右侧区域压力梯度较大，与速度场云图有一定的相似性，同样有助于理解图4.1(a)流线图的形成机理。图4.1(d) Re=10000 温度场云图给定的温度条件是左壁面温度0K，下壁面50K；上、右壁面为100K。由于我们设置钝体导热系数无穷大，故其内部温度等于下壁面温度，可在云图中看到钝体表面区域温度均与下壁面相同，对温度场分布有明显影响。（二）Re=1000 图4.2(a) Re=1000 流线图图4.2(b) Re=1000 速度场云图图4.2(c) Re=1000 压力场云图图4.2(d) Re=1000 温度场

11、云图4.3 Re=2000图4.3(a) Re=2000 流线图图4.3(b) Re=2000 速度场云图图4.3(c) Re=2000 压力场云图图4.3(d) Re=2000 温度场云图4.4 Re=5000图4.4(a) Re=5000 流线图图4.4(b) Re=5000 速度场云图图4.4(c) Re=5000 压力场云图图4.4(d) Re=5000 温度场云图比较四种不同雷诺数下的流体流动情况：相同点：右侧大漩涡包含了整个流场的大部分能量，另一小部分集中在顶盖中部附近的涡中，而左侧大漩涡虽然尺寸较大，但是流速非常慢，压力梯度也较小。流线图、速度场云图、压力场云图都可说明钝体右边的大漩涡和顶盖中部附近的小漩涡是除四个角涡以外的其他漩涡产生的动力来源。不同点：随着雷诺数的减小，在钝体上面的两个小漩涡会消失。五、后续工作展望本文利用Fortran程序计算了方腔拖曳流在雷诺数分别为1000、2000、5000及1000