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

内部存在钝体的平板拖曳流数值模拟分析一、问题介绍　　如图有一矩形区域，宽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算法是一种经典的基于压力的半隐式数值模拟方法，它的计算过程大致可分为以下几步：（1）假定一速度分布，记为,以此计算动量离散方程中的系数及常数项；（2）假设一个压力场；（3）依次求解动量方程，得；（4）求解压力校正方程，得；（5）根据校正压力改进速度值；（6）利用改进后的速度场求解那些通过源项物性等与速度场耦合的变量，如果 Φ变量并不影响流场，则应在速度场收敛后再求解；（7）利用改进后的速度场重新计算动量离散方程的系数，并利用改进后的压力场作为下一层次迭代计算的初值。重复上述步骤，直到获得收敛的解。三、程序设计简要介绍（一）主程序模块LOGICAL LSTOP ! 逻辑变量 LSTOP，为真时计算停止COMMON/CNTL/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 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)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,130DO 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 J=2,M2U(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 GAM1=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 时调用扩散系数为 kE