SIMPLE算法求解方腔内粘性不可压流动

1、SIMPLE算法求解方腔内粘性不可压流动m录一、问题描述2二、离散格式3交错网格3方程离散4三、SIMPLE算法基本思想 7边界条件处理8虚拟网格处理9方程求解11输出变虽处理12SIMPLE算法流程图15四、程序中主要变虽的意义 15五、计算结果与讨论 17函数最大值17变虽等值线图18主要结论22六、源程序22问题描述假设0 <x,y <1的方腔内充满粘性不可压缩流体，左、右、下壁固定,22ll = -1上壁以1所示。u = 16 x（1 x ）运动，试求Re =100 , 200 , 400 时的定常解，万腔如图图1方腔内流动示意图:、离散格式本算例采用求解不可压缩流动的经典

2、算法，即SIMPLE算法，求解方腔内粘性不可压缩流体运动的定常解。SIMPLE算法的全称为 Semi-Implicit Method for Pressure-Linked Equations , 即求解压力关联方程的半隐式算法。采用SIMPLE算法时，为了避免中心差分格式将“棋盘”型参量分布误认为是均匀分布， 需要用交错网格对计算域进行离散。交错网格交错网格如图2所示，压力、密度等物理量存储在控制体（i,j ）的中心，这个控制体称为主控制体。速度分量 （u , V ）分别存储在主控制体的 （i +1/ 2, j ）和（i, j +1/ 2 ）位置处， 标记为（i, j ）位置，再分别以此为中

3、心，划分速度分量 u、v的控制体。采用空间均匀网格， 等间距离散整个求解域，如图 3所示。图2交错网格示意图<AA) <AA) <:A,A) (A；“50> P也9吃野A)AO <AA) (_A_> CA）涝财A5E oAo (A) (A 一J,P：鸟 JcA)C一A,> cA>oAo <AEJ) <ALi) (Au) (A) (LJ> oAVQA 邳术知Avi.rk). <A)(Ap/ %A)oAU0T0<vt:oPi p 知An-iLJ) <AEJ> cALl <Lh)图3求解域离散示意图图3中

4、阴影部分代表方腔内的流动区域，阴影区域的边界代表方腔的上、下、左、右壁面，阴影区域外面的网格节点是为边界处理需要而设定的虚拟网格节点，后面介绍边界处理方法时详细论述。方程离散无量纲化的守恒型不可压缩N S方程为'、 *U =0；U_1 o , U U P - U = 0：tRe其积分形式为, n UdS =0S；:u1 . . 1 一dV -，in U udS i pnxdS - n udS =0V ftSSRe S31,一 dVi n UvdSi pnydS - : n '、vdS =0V .:t SSRe S图4主控制体图5 速度u控制体图6速度v控制体采用有限体积法离散N

5、 - S方程，连续性方程在主控制体上离散M 1 M 1.M 1 M 1Ui'j -uu,j ：y Vi'j -Vi'j.：x =°x方向动量方程在速度.'：x：ym 1t j'Y方向动量方程在速度u控制体上离散，时间采用前差-UiMj FI'1 FiL 祠 Gj Gi'.Lx piMv控制体上离散，时间采用前差M-pi,j -y =0L x L y M 1 MV|,jV|,j.:tF 2 -F 2 x G 2 -G 2 "v pM 1 -pM 1 x=0i,j I J, j x I ,j i,j J . v pi,j5

6、 pi,j .人其中，数值通量1 ；：u11 ；:u, G = uv Re ;xRe : y1 :v 21 ;v, G = uv Re ;xRe ;y通量F（）G（）分别定义在主控制体的中心和角点，如图所示，并按照如下方法离散F 1 =u2F 2 =v2M M -1 M 1 ui'j ui 1,j ui'j)-M 1 M 1 ui -1'j Ui'jRexG 1 = 1 ViM1,j4M MT M 1 Vi'j ui 1,j - ui'j)-M 1 M 1 ui 1'j -ui'jRey通量F f）G?）分别定义在主控制体的中心和

7、角点，如图所示，并按照如下方式离散F C= 1 (v4G ()= l(u4MMi 1,j vi,jM 1Vi 1,jM .1-Vi,j-1Re _ yM 1M .1vE1, j vi, jMMM 1M 11M 1Miui 1, jui, jvi1,j-vi,j)-vi -1,j-vi,j通量F C）G（而F（2 ） G（协某些项冻结于 M时间层，使离散化之后的方程对uM *, vM +是线性的。将离散化之后的F（）,G（届F（2）,G（2 '弋入离散后的x方向和y方向的动量方程,整理之后得离散后的动量方程如下uai , jM -1 u I, jua p,quM 1p,qubi, jM

8、1 pi1,j -M !；-pi,j1yvM1vM 1vM耶M -1 .: xai, jvI, ja p,qvp,q bi, jpi,j 1 -pi,j其中上x. ：v:.tMui,juai 1,juai,j 11 M =y Ui 1,j _41 M =*x _ vi,j _4Mu u i, jM vi .1,j1 Mai,j = y ； ui 1,j,juai-1, jMui,jM'ui-1,juai ,j1 M=-.、：x |： vi ,j _1_4M, vi -1, j Ji,jRe ： x J 4, vi 1,j 一 vi,j _1 - vi 1 j q)十 2 1 十 AxA

9、yRe "ybiVjtMv.I, j1 Mai,j 1 = ."X Vi,j_4v .I 1 Mai 1,j= y - Ui,j 1_4M'Vi,j 1Re Lyvai,j 4MUi,jARexvai,j=-x - viMj J - viMj-IL41 M=-y ”,j 1Re LyM uj,j -Rexv .1 MMai,j 7'x vi,j 1 -ui,j一4Re *y J1 MM M'V - ui,j 1 ui,j u"j 1M uIj2_x：y以上是SIMPLE算法中离散化的动量方程L SIMPLE算法基本思想SIMPLE算法是一种

10、解决压力-速度耦合问题的“半隐式”算法。首先给定M时刻猜测的速度场uM ,vM ,用于计算离散动量方程中的系数和常数项。给定 M+1时刻猜测的压力场 估计值p* ,迭代求解离散动量方程，得到 M+1时刻速度场的估计值 u* ,v* ,速度场的估计*值u ,v满足如下离散万程。u *u *u*a u 一、 a u -b -ip. p ：v = 0i, j i,jp,q p,q i, Jpi 1, J pi,j l yv * v *. v*a v 一 ： a v b i p p ：x = 0i, j i, jp,q p,q i, j p i, j -1pi , jw因而需要对速度场u , v和压力

11、场p 、一 * *般地，速度场u ,v不满足离散的连续性方程,进行修正。M+1时刻的修正值和估计值有如下关系M福1uM1vM -f1p其中，u ： v '和p '分别速度和压力的修正量，修正量亦满足离散的动量方程au u 、 au u bu - j p p ： v = 0i, j i, jp ,q p ,q i, j pi 1, j p i ,j vav v - * av v bv - 1 p p x = 0i, j i, jp,q p,q i, j p i, j -1pi , j编号为（i,j ）的速度修正量u : v '不仅与压力修正量 p有关，还与邻近点的速度修正

12、量有关。SIMPLE算法的重要假定：速度的改变只与压力的改变有关，忽略邻近点对速度修 正的影响。因而得到如下速度修正量匚yui,j =pi n 一 pi,jai,jLXvi,j = 一 v Pi,j 1 一 pi,jai,jM : u i, ji*:y.=ui,j - uPi .1,jai,j- pi,jM 1 Vi,j*x- vi, j - vPi,j 1ai,j- p,j修正后的速度分量将修正后的速度分量代入离散后的连续性方程，得到压力修正方程其中ai 1, ji,ji,jaai,jPi,ji l,ja -i .1, jp, qpp,qi,jini J,ji ,ji,ji, j Jp*bi

13、,j = - ：y ui,j -山$- xvi,j I”采用迭代法求解压力修正方程，得到压力修正量 p，代入修正公式得到 M+1时刻的速度场uM*,vM*和压力场pM*。将M+1时刻的速度场uM*,vM*和压力场pM*作为新的猜测的速度场和猜测的压力场估计值，采用上述方法计算下一个时刻的速度场和压力场，直到满足收敛条件。收敛判据Max (b" )< &&为很小的正实数，视计算的精度要求而定。本算例中取8 =1e - 8。若b" = 0 ,贝U p j = 0 ,此时U " * = U * j ,从而来自于离散动量的 U j, j满足离散的连续

14、性方程。因此Max (bipj )< &可以作为收敛判据。边界条件处理首先对计算区域离散，并流动边界之外扩充一个虚拟网格，将真实流动的离散域包围, 如图7所示。7虚拟网格扩充示意图压力p的存储位置由图中符号, 代表，速度分量u的存储位置由图中符号 -代表,速度分量v的存储位置由图中符号 代表。阴影区域表示真实的流动区域，阴影区域外的网格为虚拟网格。速度分量u在x方向虚拟网格上的存储位置分别记为u,j, un.j，在y方向虚拟网格上的存储位置分别记为ui,0 , ui,n*。速度分量V在X方向虚拟网格上的存储位置分别记为Vo,j, Vn+.j ，在y方向虚拟网格上的存储位置分别记为

15、Vi, ui.n书。压力和压力修正量在X方向虚拟网格上的存储位置分别记为Po,j,Pn*j, p0,j, p：毗j，在y方向虚拟网格上的存储位置分别记为PioPi", Pi'。，Pi*。虚拟网格处理静止壁面的虚拟网格按图 8所示方式处理，可以满足粘性边界条件，以下边界为例说 明。边界外的虚拟网格存储位置上的速度分量，由边界内的真实网格对应存储位置上的速 度分量沿边界对称得到。这样可以保证虚拟网格与真实网格在对应存储位置上的速度分量大 小相等，方向相反，二者在边界上的线性平均值为零，满足粘性边界条件。上壁面有速度，为运动壁面，速度分量U #0, V = 0。速度分量v采用前述对

16、称方法处理。速度分量U的处理方法略微不同。上壁面速度Uupboundary不随时间变化，虚拟网格上速度分量U i,n噂由U up boundary , U成线性插值得到，可以保证上壁面满足边界条件，即U i,n *= 2 * u up boundary 一 Ui,n_J。存储位置位于恰好位于静止壁面上的速度分量恒为零，不必迭代计算，例如u0,j =0。边界外虚拟网格中心的压力和压力修正量，认为分别等于边界附近对应真实网格中心的压力和压力修正量，即假定区域边界附近沿边界法向的压力梯度为零。因此，只需求解真实网格上各个存储位置上变量的值，虚拟网格上各个存储位置上变量的值可以通过前述处理虚拟网格的方

17、式得到。至此，包括虚拟网格在内的所有存储位置上都被赋予了确定的值。方程求解X方向动量方程*迭代法求解x万向动量万程，若要得到M+1时刻Ui 的值，需要用到 M时刻,JM MMMMM M MMUi,J,U+，U4，Ui,j*UiH , Vi,j,Vi*,Vi,j,Vj*上九个位置的速度分重值和M+1时刻* * 一 . 、 、 . *p ,p.两个位置的压力值。内点上的值，可以直接迭代得到。求M+1时刻速度分量U1 1 ,i ! I1, Ji , J,7m kh i-TI* -h"i f M MM M MM M M Mrz.,、». , /、,、.八(=t / (.>需要

18、用到M时刻U1,1 ,U2,1,Uo,1,U1,2,U1, 0, V1, 1,v2,1,v1,0,v2,0九个位置的速度分量值和M+1时*刻P2,1 , P1,1两个位置的压力值，包括虚拟网格在内的离散域满足计算需求。同理可得，在其他边界位置，离散域满足求解X方向动量方程的需要。Y方向动量方程迭代法求解y方向动量方程，若要得到M+1时刻v"的值，需要用到M时刻M M M M MM MMMVi,j , Vf J, Vi,j,Vi,j*,Vi,j，Ui,j ,Ui+,j ,Ui,j，Ui, J 九个位置的速度 重值和 M+1 时刻* * 一 . 、 、 . *PijPij两个位置的压力值

19、。内点上的值，可以直接迭代碍到。求M+1时刻速度分重v1,1 ,M M M M Mm M M M一，需要用到M时刻V1,1,V2,1,Vo,1,V1,2,V1,o , U1,1 ,U2,1 ,U1,0 ,U2,o九个位置的速度分量值和M+1时、.*刻P1,2 , P1,1两个位置的压力值，包括虚拟网格在内的离散域满足计算需求。同理可得，在其他边界位置，离散域满足求解y方向动量方程的需要。压力修正方程迭代法求解压力修正方程，若要M+1时刻得到p：的值，需要用到M时刻Pi,j, Pi'j，PJ+，Pi：J四个位置的压力修正值和M+1时刻*Ui,j , U i + j , Uj,j, U i

20、, j , Uj,j，Ui,j+ Ui,j，Uiu,j ,*Vi, J，Vi,J1， Vi , J J， ViT, J，Vi M, J，ViT, J ，V" J W，Vi,J -2十六个位置的速度分量估计值。内点上的值，可以直接迭代得到。 求M+1时刻压力修正值P；1，需要用到M时刻P2,1, p0,1 , P；2, p1,0四个位置的压力修正值和M+1时刻* *U1,1 ,U 2,1>U0,1 ,U1,2 ,U1,0 ,U0 ,2 >U0,0 ,U _1,1* *V1,1 ,V1,2,V1,0 ,V2 ,1 ,V0,1 ,V2 ,0 ,V0, 0 ,V1, J十六个位置

21、的速度分量估计值，包括虚拟网格在内的离散域满足计算需求。同理可得, 在其他边界位置，离散域满足求解压力修正方程的要求。将M+1时刻变量的估计值和修正值代入如下公式M u1*-=uuM v1*-=vvMP 1*=pp一一-、一M -即可得到M+1时刻变量u1M -,v1M -1,p的值。输出变虽处理涡量定义v c u£0 =CL': xc y涡量是一个矢量，在平面流动中只有一个方向，垂直于流动平面。流函数满足v =,u -f xf y流函数等值线为流线兴d =dxdy=-vdxudyxy在交错网格里，对涡量和流函数在一个网格内做中心差分，涡量和流函数的存储位置为真实网格的节点，

22、如图 9所示。图9涡量和流函数存储位置示意图涡量差分格式vi 1 j 一 vi, j ui i,j 一 ui,j co =xyy方向流函数差分格式= yui, jyu i, j i, j -ix方向流函数差分格式"i,j = - XVi,j"i_i,j任取其中一个式子计算流函数即可，或者将两个式子的计算结果取平均值。本算例采用的是取平均值的办法。边界上的流函数值为零。需要输出用于后处理的变量主要有位置坐标x、y,速度分量u、v,压力p,流函数巾和涡量切。由于SIMPLE算法采用交错网格，不同性质的变量存储在网格的不同位置，在变量输出时，有必要将所有变量统一到真实网格的网格节

23、点上，便于后处理，如图10所示。图10变量统一后存储位置示意图c 1质=2 ui,j ' ui,j 1网格节点上的速度分量u?为节点上下两速度分量的平均值，即网格节点上的速度分量。为节点左右两速度分量的平均值，即1= vi,j vi 1,j2网格节点上的压力 ？为节点周围四个主控制体中心压力的平均值,1p =- Pi,j - Pi,j 1 - Pi 1,j- Pi 1,j 14边界处的压力平均采用相同的公式，需要用到虚拟网格中心的压力，角点依然。SIMPLE算法流程图统一变量存储位置，后处理J结是图11 SIMPLE算法流程图四、程序中主要变量的意义主程序中主要变量的意义变量意义cou

24、nter_max预设最大循环次数grid_number、nX方向和y方向的离散网格数目t、h时间步长和空间步长Re雷诺数d收敛判据，压力修正方程中bj绝对值的最大值i , J数组u、vX方向和y方向速度分量数组 p、p_wave压力和压力修正量数组 flow_function、 vortex_function流函数和涡量函数数组 u_scalar速度绝对值数组 d_array、g、u_0计算中需要用到的一些辅助数组子程序1,即求解x方向动量方程程序中主要变量的意义变量意义nX方向和y方向的离散网格数目t、h时间步长和空间步长Re雷诺数a、b、c、d、e迭代求解x方向动量方程所需的辅助变量 一数

25、组u、v、pX方向、y方向速度分量和压力场子程序2,即求解y方向动量方程程序中主要变量的意义变量意义nX方向和y方向的离散网格数目t、h时间步长和空间步长Re雷诺数a、b、c、d、e迭代求解x方向动量方程所需的辅助变量 一数组u、v、pX方向、y方向速度分量和压力场子程序3,即求解压力修正方程程序中主要变量的意义变量意义nX方向和y方向的离散网格数目t、h时间步长和空间步长Re雷诺数a、b、c、e、f压力修正方程中各压力修正项系数nd压力修正方程中0p 项 , Ja_u1、 a_u2辅助变量，与u方向速度修正式有关a_v1、 a_v2辅助变量，与v方向速度修正式有关error求解压力修正方程的

26、迭代收敛判据数组u、vX方向、y方向速度分量数组 p、p_wave压力和压力修正量五、计算结果与讨论本算例分别对网格数目为 20 X 20,40X 40两种情况，数值求解了 Re=100,200,400时的方 腔流动的定常解，计算结果列举如下。函数最大值离散网格为20X 20数值计算得到的流函数最大值Re100200400流函数最大值0.29269532920.27242897280.2427700376对应x坐标0.400.450.45对应y坐标0.650.600.60离散网格为20X 20数值计算得到的上壁面涡函数最大值Re100200400上壁面涡函数最大值91.722027587810

27、4.8169276505117.6126516741对应x坐标0.800.750.75对应y坐标1.001.001.00离散网格为20X 20数值计算得到的方腔中线涡函数最大值Re100200400中线上涡函数最大值25.400771350538.822847263853.8979317133对应x坐标0.500.500.50对应y坐标0.951.001.00离散网格为40X 40数值计算得到的流函数最大值Re100200400流函数最大值0.31935425740.31347035100.2972193667对应x坐标0.4250.4500.475对应y坐标0.6250.6000.575离散

28、网格为40X 40数值计算得到的上壁面涡函数最大值Re100200400上壁面涡函数最大值119.0699676898149.6299870265181.4581848045对应x坐标0.8250.8000.800对应y坐标1.001.001.00离散网格为40X 40数值计算得到的方腔中线涡函数最大值Re100200400中线上涡函数最大值32.418443541041.919028264845.2529732845对应x坐标0.500.500.50对应y坐标0.9750.9750.975从表中数据可以看到，随着雷诺数增大，流函数最大值减小，最大流函数的位置也更趋 于方腔的中心。这是因为雷诺

29、数增大， 粘性效应减弱，上壁面运动带动放腔内流体运动的能 力降低。上壁面涡量最大值位于上壁面的速度最大处，这是因为在该处速度分量 u沿y方向的的梯度最大，速度分量v沿x方向的梯度为零（壁面粘性边界条件），因而剪切效应最强。 随着雷诺数增大，上壁面涡量最大值增大。中线上涡量最大值出现在上壁面附近，但不位于上壁面。离散网格为 20X 20,雷诺数为100时，数值计算结果表明，中线上最大涡量出现 在y=0.95处；雷诺数为200和400时，数值计算结果表明，中线上最大涡量出现在y=1.00处。离散网格为40X 40,对于雷诺数为100、200和400三种情况，数值计算结果表明，中 线上最大涡量出现在

30、 y=0.975处。这说明增大离散网格数目，即加密网格，可以增加数值计 算对涡量最大值的捕捉精度。变虽等值线图速度大小等值线图12 Re=100时速度大小等值线W 6 2 o O 0 d二图15 Re=100时流函数等值线图13 Re=200时速度大小等值线图14 Re=400时速度大小等值线流函数等值线FLOW-F? K.- J-:- & -b 2，-J .宝 2 2- .1 .S D U。Q Q 00.0.14 0.12 0.1 0 08AM 0.02 0 OOH 0 0061 0G通 -UflC&4 -0.U12503 。为 02 口网 022口诵 口 16 0 14 0

31、.12 0.1 0.09 0.00 0.M 0.02 0.00770AM1 Q-0Q125J3.00M fl.0125Fl W-F网3也M.22.2怙怖14SRes200FLCWF=gQ" 0.00&10 1 0.0fi! OJM 冲0.0077 0.OT51 0图16 Re=200时流函数等值线图18 Re=10 0时涡量等值线20X20ReMOOFLg=成 M=。册 0 16 0 14B 0 1210 1 0® 006 DM.GQ2 H 血oott Ottfli . C:l -Dm .-四维 -M125PLCW FR«-4CCn761GMJ5 21-e

32、禅U121MIKM02汕网浙顷" 以o.QQGUQ.o.o.o.oa.。"#图17 Re=400时流函数等值线涡量等值线20X20£,04-4fMTTDJLm图19 Re=200时涡量等值线图20 Re=400时涡量等值线从以上各个变量的等值线图中可明显看出， 相同雷诺数不同网格密度时， 各个变量相应 的等值线形状基本相同， 但40 X 40网格的等值线明显比 20 X 20网格的等值线光滑。 这说明 加密网格可以提高数值计算的精度。从流函数等值线中可以看到，Re=100和200时，方腔左下角处流函数均匀，Re=400时,方腔左下角处流函数不均匀，出现了与方腔中心

33、流函数等值线不一致的流函数等值线。这说明在增大雷诺数，流体粘性减弱的情况下，方腔下壁面角点处可能出现回流。因此本算例对Re=1600的情形在40X 40网格上做了数值计算，结果如图所示。图中流函数等值线分布表 明，在方腔左下角附近，确实出现了回流，但是回流区的流函数值很小，说明回流强度很小。40X40FLOW-FO.O.O.O.O.C5O.O.0.O.O.0.0077O.OQ61 0 -0.0046 -0.0084-0.0125图21 Re=1600时流函数等值线主要结论加密离散网格，可以提高数值计算的精度。随着雷诺数的增加， 流体粘性效应减弱，上壁面驱动流体的能力减弱，放腔内流函数整体减小，

34、下壁面角点附近出现回流，上壁面涡量增大。壁面上的涡量可以为正， 也可以为负，涡量的最大值出现在速度分量u沿y方向梯度与速度分量v沿x方向梯度只差最大的地方，因而不一定出现在壁面上。六、源程序program mainimplicit noneinteger , parameter : counter_max=50000 !设置取大的时间推进步数为五万步 real(8) , parameter : standard=0.00000001!小数点后取八位integer : i,j,i_max,j_max,n,counter,grid_number!grid_numer 为网格数目real(8) ,

35、allocatableu(:,:),v(:,:),p(:,:),p_wave(:,:),d_array(:,:),g(:),u_0(:,:),flow_function(:,:),flow_function_u(:,:),flow_function_v(:,:),vortex_function(:,:),u_scalar(:,:)real(8) : h,t,Re !h 为空间步长，t 为时间步长，R刮 Reynolds number real(8) : d,p_wave_maxwrite (*,*)"输入离散网格数目：”read (*,*)grid_numberwrite (*,*)

36、" 输入 Reynolds number :"read (*,*)Recounter=1 !迭代步数计数器h=1.0/real (grid_number)n=grid_numbert=0.001 allocate (u(-1:n+1,0:n+1),v(0:n+1,-1:n+1),p(0:n+1,0:n+1),p_wave(0:n+1,0:n+1),d_array(n,n),g(0:n),u_0(-1:n+1,0:n+1),flow_function(0:n,0:n),vortex_function(0:n,0:n),u_scalar(0:n,0:n),flow_functi

37、on_u(0:n,0:n),flow_function_v(0:n,0:n)!声明二维数组的长度，p_wave代表修正压力u=0!赋初值v=0do i=1,n-1u(i,n+1)=-16.0*(real (i)*h)*2*(1.0-(real (i)*h)*2)*2-u(i,n)end dop=0p_wave=0end dod=1.0do while(d>=standard)if (counter>counter_max) exit!如果时间推进步数大于五万步，跳出循环u0=ucallsub1_x_momentum(u,v,p,h,t,n,Re)!求解离散的x方向动量方程calls

38、ub2_y_momentum(u_0,v,p,h,t,n,Re)!求解离散的y方向动量方程callsub3_pressure(u,v,p,p_wave,n,t,Re)!求解压力修正方程do i=1,ndo j=1,nd_array(i,j)= 10.0*abs(u(i,j)-u(i-1,j)+v(i,j)-v(i,j-1)end doend dod=maxval (d_array)!d 作为收敛判据do j=1,n d_array(i,j)=abs(p_wave(i,j) !只考虑实际的压力波动，即边界内的压力波动end do end do p_wave_max=maxval (d_array

39、)counter=counter+1 write (*,*)"d=",d,"p_wave_max=",p_wave_max,counter-1 end do write (*,*)counter-1flow_function_u = 0!求流函数flow_function_v = 0 do i=1,n-1 do j=1,n flow_function_u(i,j)= h*u(i,j)+flow_function_u(i,j-1) end do end do do j=1,n-1 do i=1,n flow_function_v(i,j)= -1.0*h*

40、v(i,j)+flow_function_v(i-1,j) end do end do do i=0,n do j=0,n flow_function(i,j)= 0.5*(flow_function_u(i,j)+flow_function_v(i,j) end do end do do i=1,n-1!求涡量函数!left walldo j=1,n-1 vortex_function(i,j)= (v(i+1,j)-v(i,j)/h - (u(i,j+1)-u(i,j)/h end do end do !求边界上的涡量函数 do j=1,n-1 vortex_function(0,j)=

41、(v(1,j)-v(0,j)/h - (u(0,j+1)-u(0,j)/hdo j=1,n-1vortex_function(n,j)= (v(n+1,j)-v(n,j)/h - (u(n,j+1)-u(n,j)/hend dodo i=1,n-1vortex_function(i,n)= (v(i+1,n)-v(i,n)/h - (u(i,n+1)-u(i,n)/hend dodo i=1,n-1vortex_function(i,0)= (v(i+1,0)-v(i,0)/h - (u(i,1)-u(i,0)/hend do!right wall!up wall!down wall!计算角点

42、涡量函数vortex_function(0,0)= 0.5*(vortex_function(0,1) + vortex_function(1,0)vortex_function(0,n)= 0.5*(vortex_function(0,n-1) + vortex_function(1,n)vortex_function(n,0)= 0.5*(vortex_function(n-1,0) + vortex_function(n,1)vortex_function(n,n)= 0.5*(vortex_function(n,n-1) + vortex_function(n-1,n)do j=0,n

43、!将u速度节点与网格节点重合do i=0,ng(i)=(u(i,j+1)+u(i,j)*0.5end dodo i=0,nu(i,j)=g(i)end doend dodo i=0,n!将v速度节点与网格节点重合do j=0,ng(j)=(v(i+1,j)+v(i,j)*0.5end dodo j=0,nv(i,j)=g(j)end doend dodo j=0,n!求节点速度绝对值do i=0,nu_scalar(i,j)=sqrt (u(i,j)*u(i,j) + v(i,j)*v(i,j)end doend dodo i=0,n!将压力节点与网格节点重合do j=0,np_wave(i,

44、j)=0.25*(p(i,j)+p(i+1,j)+p(i,j+1)+p(i+1,j+1)end do end dop=p_waveopen(unit=10,file="data.dat")write (10,”('TITLE = ""Grid=',I3,',Re=',f5.1,""")")n,Re!输出双引号时需要打出两个双引号，否则出错write (10,*)'VARIABLES = "X" "Y" "U" &qu

45、ot;V" "U-SCALAR" "P" "FLOW-F" "VORTEX-F" '!这种情况不需要打两个双引号就可以输出双引号write (10,"(/)")write (10,"('ZONE N=',I6,', E=',I6,', ZONETYPE=FEQuadrilateral,DATAPACKING=POINT')")(n+1)*(n+1),n*ndo j=0,ndo i=0,nwrite (10,&q

46、uot;(f15.10,f15.10,f15.10,f15.10,f15.10,f15.10,f15.10,f15.10)")real (i)*h, real (j)*h,u(i,j),v(i,j),u_scalar(i,j),p(i,j),flow_function(i,j),vortex_function(i,j)end doend dodo j=0,n-1!对网格节点进行编号，便于 tecplot软件后处理i= j*(n+1)counter=1do while (counter<=n)write (10,"(I5,I5,I5,I5)")i+counte

47、r, i+counter+1, i+counter+1+n+1, i+counter+n+1counter=counter+1 end doend doopen(unit=11,file="function_max.txt")d=0.0do j=0,n!流函数最大值do i=0,nif (d<=flow_function(i,j) thend=flow_function(i,j)i_max=ij_max=jend ifend doreal (i_max)*h, real (j_max)*hwrite (11,*)"flow_functioin_max:&qu

48、ot;write (11,”(f15.10,I5,I5,f15.10,f15.10)")d,i_max,j_max,d=0.0!上壁面涡函数最大值do i=0,nif (d<=vortex_function(i,n) thend=vortex_function(i,n)i_max=iend ifend dowrite (11,*)"up_wall_vortex_functioin_max:"write (11,”(f15.10,I5,f15.10)")d,i_max,real (i_max)*hd=0.0!中竖线上涡函数最大值do j=0,nif

49、(d<=vortex_function(n/2,j) thend=vortex_function(n/2,j)j_max=jend ifend dowrite (11,*)"middle_line_vortex_functioin_max:"write (11,”(f15.10,I5,f15.10)")d,j_max,real (j_max)*hstopendsubroutine sub1_x_momentum(u,v,p,h,t,n,Re) implicit none integer : i,j,n!注意保证数组传递相匹配real(8) : u(-1:n+

50、1,0:n+1),v(0:n+1,-1:n+1),p(0:n+1,0:n+1)real(8) : h,t,a,b,c,d,e,f,Redo j=1,ndo i=1,n-1a=h*(0.25*(u(i+1,j)-u(i-1,j)+2.0/(h*Re) + h*(0.25*(v(i,j)+v(i+1,j)-v(i,j-1)-v(i+1,j-1)+2.0/(h*Re) + h*h/tb=h*(0.25*(u(i,j)+u(i-1,j)+1.0/(h*Re) * u(i-1,j)c=h*(0.25*(u(i+1,j)+u(i,j)-1.0/(h*Re) * u(i+1,j)d=h*(0.25*(v(i

51、,j-1)+v(i+1,j-1)+1.0/(h*Re) * u(i,j-1)e=h*(0.25*(v(i,j)+v(i+1,j)-1.0/(h*Re) * u(i,j+1)f=h*h/t*u(i,j) - h*(p(i+1,j)-p(i,j)u(i,j)=(b-c+d-e+f)/aend dou(-1,j)=-1.0*u(1,j)u(n+1,j)=-1.0*u(n-1,j)end dodo i=1,n-1u(i,0)=-1.0*u(i,1)u(i,n+1)=-16.0*( real (i)*h)*2*(1.0-( real (i)*h)*2)*2-u(i,n) end doreturnend subroutine sub1_x_momentumsubroutine sub2_y_momentum(u,v,p,h,t,n,Re) implicit noneinteger : i,j,n!注意保证数组传递相匹配real(8) : u(-1:n+1,0:n+1),v(0:n+1,-1:n+1),p(0: