方腔驱动流源代码(SIMPLE算法)

1、simplesimple 算法求解方腔内粘性不可压流动算法求解方腔内粘性不可压流动 一、问题描述 假设1,0yx 的方腔内充满粘性不可压缩流体，左、右、下壁固定，上壁以22116xxu 运动，试求400,200,100re 时的定常解，方腔如图 1 所示。 图 1 方腔内流动示意图 二、离散格式 本算例采用求解不可压缩流动的经典算法， 即 simple 算法， 求解方腔内粘性不可压缩流体运动的定常解。 simple 算法的全称为 semi-implicit method for pressure-linked equations，即求解压力关联方程的半隐式算法。 采用 simple 算法时，

2、为了避免中心差分格式将 “棋盘” 型参量分布误认为是均匀分布，需要用交错网格对计算域进行离散。 交错网格 交错网格如图 2 所示，压力、密度等物理量存储在控制体ji,的中心，这个控制体称为主控制体。速度分量vu,分别存储在主控制体的ji, 2/1 和2/1,ji 位置处，标记为ji, 位置，再分别以此为中心，划分速度分量 u、v 的控制体。采用空间均匀网格，等间距离散整个求解域，如图 3 所示。 图 2 交错网格示意图 图 3 求解域离散示意图 图 3 中阴影部分代表方腔内的流动区域，阴影区域的边界代表方腔的上、下、左、右壁面， 阴影区域外面的网格节点是为边界处理需要而设定的虚拟网格节点， 后

3、面介绍边界处理方法时详细论述。 方程离散 无量纲化的守恒型不可压缩sn 方程为 0re102upuutuu 其积分形式为 ssysvssxsvsvdsndspnvdsundvtvudsndspnudsundvtudsun0re10re10 图 4 主控制体 图 5 速度 u 控制体 图 6 速度 v 控制体 采用有限体积法离散sn 方程，连续性方程在主控制体上离散 011,1,1, 11,xvvyuumjimjimjimji x 方向动量方程在速度 u 控制体上离散，时间采用前差 01,1, 111,1,1, 11,1,yppxggyffuutyxmjimjijijijijimjimji y

4、方向动量方程在速度 v 控制体上离散，时间采用前差 01,11,21,2,2, 12,1,xppyggxffvvtyxmjimjijijijijimjimji 其中，数值通量 yuuvgxuufre1,re1121 yvuvgxvvfre1,re1222 通量 11,gf 分别定义在主控制体的中心和角点，如图所示，并按照如下方法离散 1,1, 11,1, 1, 111,1, 11,1, 1, 11re141re141mjimjimjimjimjimjimjimjimjimjimjimjiuuyuuvvguuxuuuuf 通量 22,gf 分别定义在主控制体的中心和角点，如图所示，并按照如下方式

5、离散 1,1, 11,1, 1, 121,1, 11,1, 1, 12re141re141mjimjimjimjimjimjimjimjimjimjimjimjivvxvvuugvvyvvvvf 通量 11,gf和 22,gf的某些项冻结于 m 时间层，使离散化之后的方程对11,mmvu 是线性的。将离散化之后的 11,gf和 22,gf代入离散后的 x 方向和 y 方向的动量方程，整理之后得离散后的动量方程如下 001,11,1,1,1,1, 1,1,1,xppbvavayppbuauamjimjivjimqpvqpmjivjimjimjiujimqpuqpmjiuji 其中 mjiujiu

6、tyxb, yvvxayvvxaxuuyaxuuyamjimjiujimjimjiujimjimjiujimjimjiujire141,re141re141,re1411, 11,1, 1,1, 1 -, 1 -, 1, 1 tyxyvvvvxxuuyamjimjimjimjimjimjiujire241re2411, 11, 1, 1, 1, mjivjivtyxb, xuuyaxuuyayvvxayvvxamjimjivjimjimjivjimjimjivjimjimjivjire141,re141re141,re141, 11, 1, 1,1, 1,1,1,1,1, tyxxuuuuyy

7、uvxamjimjimjimjimjimjivjire241re241, 11, 1,1,1,1, 以上是 simple 算法中离散化的动量方程 三、simple 算法基本思想 simple 算法是一种解决压力-速度耦合问题的“半隐式”算法。首先给定 m 时刻猜测的速度场mmvu ,，用于计算离散动量方程中的系数和常数项。给定 m+1 时刻猜测的压力场估计值*p，迭代求解离散动量方程，得到 m+1 时刻速度场的估计值*,vu，速度场的估计值*,vu满足如下离散方程。 00*,*1,*,*,*,*, 1,*,*,xppbvavayppbuauajijivjiqpvqpjivjijijiujiqp

8、uqpjiuji 一般地， 速度场*,vu不满足离散的连续性方程， 因而需要对速度场*,vu和压力场*p进行修正。m+1 时刻的修正值和估计值有如下关系 pppvvvuuummm*1*1*1 其中，vu,和p分别速度和压力的修正量，修正量亦满足离散的动量方程 00,1, 1,xppbvavayppbuauajijivjiqpvqpjivjijijiujiqpuqpjiuji 编号为（i,j）的速度修正量vu,不仅与压力修正量p有关，还与邻近点的速度修正量有关。simple 算法的重要假定：速度的改变只与压力的改变有关，忽略邻近点对速度修正的影响。因而得到如下速度修正量 jijivjijijij

9、iujijippaxvppayu,1, 1, 修正后的速度分量 jijivjijimjijijiujijimjippaxvvppayuu,1,*,1, 1,*,1, 将修正后的速度分量代入离散后的连续性方程，得到压力修正方程 pjiqppqpjipjibpapa, 其中 *1,*,*, 1*,1,2,2, 12,2,1,21,21, 12, 1,2, 1,jijijijipjivjivjiujiujipjivjipjivjipjiujipjiujipjivvxuuybaxaxayayaaxaaxaayaaya 采用迭代法求解压力修正方程，得到压力修正量p，代入修正公式得到 m+1 时刻的速度场

10、11,mmvu和压力场1mp。将 m+1 时刻的速度场11,mmvu和压力场1mp作为新的猜测的速度场和猜测的压力场估计值， 采用上述方法计算下一个时刻的速度场和压力场， 直到满足收敛条件。 收敛判据 pjibmax, 为很小的正实数，视计算的精度要求而定。本算例中取81 e。 若0,pjib，则0,jip，此时*,1,jimjiuu，从而来自于离散动量的*, jiu满足离散的连续性方程。因此pjibmax, 可以作为收敛判据。 边界条件处理 首先对计算区域离散，并流动边界之外扩充一个虚拟网格，将真实流动的离散域包围，如图 7 所示。 图 7 虚拟网格扩充示意图 压力 p 的存储位置由图中符号

11、代表，速度分量 u 的存储位置由图中符号代表，速度分量 v 的存储位置由图中符号 代表。阴影区域表示真实的流动区域，阴影区域外的网格为虚拟网格。 速度分量 u 在 x 方向虚拟网格上的存储位置分别记为jnjuu. 1, 1, ， 在 y 方向虚拟网格上的存储位置分别记为1,0 ,niiuu。速度分量 v 在 x 方向虚拟网格上的存储位置分别记为jnjvv. 1, 0, ，在 y 方向虚拟网格上的存储位置分别记为1.1,niiuv。压力和压力修正量在x 方向虚拟网格上的存储位置分别记为jnjjnjpppp, 1, 0, 1, 0,，在 y 方向虚拟网格上的存储位置分别记为1, 0 ,1,0 ,n

12、iiniipppp。 虚拟网格处理 静止壁面的虚拟网格按图 8 所示方式处理，可以满足粘性边界条件，以下边界为例说明。 图 8 静止壁面虚拟网格处理示意图 边界外的虚拟网格存储位置上的速度分量，由边界内的真实网格对应存储位置上的速度分量沿边界对称得到。 这样可以保证虚拟网格与真实网格在对应存储位置上的速度分量大小相等，方向相反，二者在边界上的线性平均值为零，满足粘性边界条件。 上壁面有速度，为运动壁面，速度分量0, 0vu。速度分量 v 采用前述对称方法处理。速度分量 u 的处理方法略微不同。上壁面速度boundaryupu不随时间变化，虚拟网格上速度分量1, ,niu由1,niboundar

13、yupuu线性插值得到，可以保证上壁面满足边界条件，即1,1, ,2niboundaryupniuuu。 存储位置位于恰好位于静止壁面上的速度分量恒为零，不必迭代计算，例如0, 0ju。 边界外虚拟网格中心的压力和压力修正量， 认为分别等于边界附近对应真实网格中心的压力和压力修正量，即假定区域边界附近沿边界法向的压力梯度为零。 因此， 只需求解真实网格上各个存储位置上变量的值， 虚拟网格上各个存储位置上变量的值可以通过前述处理虚拟网格的方式得到。 至此，包括虚拟网格在内的所有存储位置上都被赋予了确定的值。 方程求解 x 方向动量方程 迭代法求解 x 方向动量方程，若要得到 m+1 时刻*, j

14、iu的值，需要用到 m 时刻mjimjimjimjimjiuuuuu1,1, 1, 1,，mjimjimjimjivvvv11,1, 1,九个位置的速度分量值和 m+1 时刻*,*, 1,jijipp两个位置的压力值。内点上的值，可以直接迭代得到。求 m+1 时刻速度分量*1 , 1u，需要用到 m 时刻mmmmmuuuuu0 , 12, 11 , 01 , 21 , 1,，mmmmvvvv0, 20, 11 , 21 , 1, 九个位置的速度分量值和 m+1 时刻 *1 , 1*1 , 2, pp 两个位置的压力值，包括虚拟网格在内的离散域满足计算需求。同理可得，在其他边界位置，离散域满足求

15、解 x 方向动量方程的需要。 y 方向动量方程 迭代法求解 y 方向动量方程，若要得到 m+1 时刻*, jiv的值，需要用到 m 时刻mjimjimjimjimjivvvvv1,1, 1, 1,，mjimjimjimjiuuuu11,1, 1,九个位置的速度分量值和 m+1 时刻*,*1,jijipp两个位置的压力值。内点上的值，可以直接迭代得到。求 m+1 时刻速度分量*1 , 1v，需要用到 m 时刻mmmmmvvvvv0 , 12, 11 , 01 , 21 , 1,，mmmmuuuu0, 20, 11 , 21 , 1, 九个位置的速度分量值和 m+1 时刻*1 , 1*2, 1,

16、pp 两个位置的压力值，包括虚拟网格在内的离散域满足计算需求。同理可得，在其他边界位置，离散域满足求解 y 方向动量方程的需要。 压力修正方程 迭代法求解压力修正方程，若要 m+1 时刻得到jip, 的值，需要用到 m 时刻1,1, 1, 1,jijijijipppp四个位置的压力修正值和 m+1 时刻 *, 2*1, 1*1, 1*1,*1,*, 1*1*,jijijijijijijijiuuuuuuuu ， *2,*1, 1*1, 1*, 1*, 1*1,*1,*,jijijijijijijijivvvvvvvv 十六个位置的速度分量估计值。内点上的值，可以直接迭代得到。求 m+1 时刻压

17、力修正值1 , 1p，需要用到 m 时刻0 , 12, 11 , 01 , 2,pppp 四个位置的压力修正值和 m+1 时刻 *1 , 1*0 , 0*2, 0*0 , 1*2, 1*1 , 0*1 , 2*1 , 1,uuuuuuuu *1, 1*0 , 0*0 , 2*1 , 0*1 , 2*0 , 1*2, 1*1 , 1,vvvvvvvv 十六个位置的速度分量估计值，包括虚拟网格在内的离散域满足计算需求。同理可得，在其他边界位置，离散域满足求解压力修正方程的要求。 将 m+1 时刻变量的估计值和修正值代入如下公式 pppvvvuuummm*1*1*1 即可得到 m+1 时刻变量111

18、,mmmpvu的值。 输出变量处理 涡量定义 yuxv 涡量是一个矢量，在平面流动中只有一个方向，垂直于流动平面。 流函数满足 yuxv, 流函数等值线为流线 udyvdxdyydxxd 在交错网格里， 对涡量和流函数在一个网格内做中心差分， 涡量和流函数的存储位置为真实网格的节点，如图 9 所示。 图 9 涡量和流函数存储位置示意图 涡量差分格式 yuuxvvjijijiji, 1, 1 y 方向流函数差分格式 1,jijijiyu x 方向流函数差分格式 jijijixv, 1, 任取其中一个式子计算流函数即可， 或者将两个式子的计算结果取平均值。 本算例采用的是取平均值的办法。边界上的流

19、函数值为零。 需要输出用于后处理的变量主要有位置坐标 x、y，速度分量 u、v，压力 p，流函数和涡量。由于 simple 算法采用交错网格，不同性质的变量存储在网格的不同位置，在变量输出时，有必要将所有变量统一到真实网格的网格节点上，便于后处理，如图 10 所示。 图 10 变量统一后存储位置示意图 网格节点上的速度分量u 为节点上下两速度分量的平均值，即1,21jijiuuu 网格节点上的速度分量v 为节点左右两速度分量的平均值，即jijivvv, 1,21 网格节点上的压力p 为节点周围四个主控制体中心压力的平均值，即 1, 1, 11,41jijijijippppp 边界处的压力平均采

20、用相同的公式，需要用到虚拟网格中心的压力，角点依然。 simple 算法流程图 图 11 simple 算法流程图 四、计算结果与讨论 本算例分别对网格数目为 2020,4040 两种情况， 数值求解了 re=100,200,400 时的方腔流动的定常解，计算结果列举如下。 函数最大值 离散网格为 2020 数值计算得到的流函数最大值 re 100 200 400 流函数最大值 0.2926953292 0.2724289728 0.2427700376 对应 x 坐标 0.40 0.45 0.45 对应 y 坐标 0.65 0.60 0.60 离散网格为 2020 数值计算得到的上壁面涡函数

21、最大值 re 100 200 400 上壁面涡函数最大值 91.7220275878 104.8169276505 117.6126516741 对应 x 坐标 0.80 0.75 0.75 对应 y 坐标 1.00 1.00 1.00 离散网格为 2020 数值计算得到的方腔中线涡函数最大值 re 100 200 400 中线上涡函数最大值 25.4007713505 38.8228472638 53.8979317133 对应 x 坐标 0.50 0.50 0.50 对应 y 坐标 0.95 1.00 1.00 离散网格为 4040 数值计算得到的流函数最大值 re 100 200 400

22、 流函数最大值 0.3193542574 0.3134703510 0.2972193667 对应 x 坐标 0.425 0.450 0.475 对应 y 坐标 0.625 0.600 0.575 离散网格为 4040 数值计算得到的上壁面涡函数最大值 re 100 200 400 上壁面涡函数最大值 119.0699676898 149.6299870265 181.4581848045 对应 x 坐标 0.825 0.800 0.800 对应 y 坐标 1.00 1.00 1.00 离散网格为 4040 数值计算得到的方腔中线涡函数最大值 re 100 200 400 中线上涡函数最大值

23、32.4184435410 41.9190282648 45.2529732845 对应 x 坐标 0.50 0.50 0.50 对应 y 坐标 0.975 0.975 0.975 从表中数据可以看到，随着雷诺数增大，流函数最大值减小，最大流函数的位置也更趋于方腔的中心。这是因为雷诺数增大，粘性效应减弱，上壁面运动带动放腔内流体运动的能力降低。 上壁面涡量最大值位于上壁面的速度最大处， 这是因为在该处速度分量 u 沿 y 方向的的梯度最大，速度分量 v 沿 x 方向的梯度为零（壁面粘性边界条件） ，因而剪切效应最强。随着雷诺数增大，上壁面涡量最大值增大。中线上涡量最大值出现在上壁面附近，但不位

24、于上壁面。离散网格为 2020，雷诺数为 100 时，数值计算结果表明，中线上最大涡量出现在 y=0.95 处；雷诺数为 200 和 400 时，数值计算结果表明，中线上最大涡量出现在 y=1.00处。离散网格为 4040，对于雷诺数为 100、200 和 400 三种情况，数值计算结果表明，中线上最大涡量出现在 y=0.975 处。这说明增大离散网格数目，即加密网格，可以增加数值计算对涡量最大值的捕捉精度。 变量等值线图 速度大小等值线 图 12 re=100 时速度大小等值线 图 13 re=200 时速度大小等值线 图 14 re=400 时速度大小等值线 流函数等值线 图 15 re=

25、100 时流函数等值线 图 16 re=200 时流函数等值线 图 17 re=400 时流函数等值线 涡量等值线 图 18 re=100 时涡量等值线 图 19 re=200 时涡量等值线 图 20 re=400 时涡量等值线 从以上各个变量的等值线图中可明显看出， 相同雷诺数不同网格密度时， 各个变量相应的等值线形状基本相同， 但 4040 网格的等值线明显比 2020 网格的等值线光滑。 这说明加密网格可以提高数值计算的精度。 从流函数等值线中可以看到，re=100 和 200 时，方腔左下角处流函数均匀，re=400 时，方腔左下角处流函数不均匀， 出现了与方腔中心流函数等值线不一致的

26、流函数等值线。 这说明在增大雷诺数，流体粘性减弱的情况下，方腔下壁面角点处可能出现回流。因此本算例对re=1600 的情形在 4040 网格上做了数值计算，结果如图所示。图中流函数等值线分布表明，在方腔左下角附近，确实出现了回流,但是回流区的流函数值很小，说明回流强度很小。 图 21 re=1600 时流函数等值线 主要结论 加密离散网格，可以提高数值计算的精度。 随着雷诺数的增加，流体粘性效应减弱，上壁面驱动流体的能力减弱，放腔内流函数整体减小，下壁面角点附近出现回流，上壁面涡量增大。 壁面上的涡量可以为正， 也可以为负， 涡量的最大值出现在速度分量 u 沿 y 方向梯度与速度分量 v 沿

27、x 方向梯度只差最大的地方，因而不一定出现在壁面上。 五、源程序 programprogram main implicit noneimplicit none integerinteger,parameterparameter : counter_max=50000 real(8)real(8),parameterparameter : standard=0.00000001 integerinteger : i,j,i_max,j_max,n,counter,grid_number real(8)real(8),allocatableallocatable : u(:,:),v(:,:),p

28、(:,:),p_wave(:,:),d_array(:,:),g(:),u_0(:,:),flow_function(:,:),flow_function_u(:,:),flow_function_v(:,:),vortex_function(:,:),u_scalar(:,:) real(8)real(8) : h,t,re real(8)real(8) : d,p_wave_max writewrite(*,*)输入离散网格数目： readread(*,*)grid_number writewrite(*,*)输入reynolds number： readread(*,*)re count

29、er=1 h=1.0/real(grid_number) n=grid_number t=0.001 allocateallocate(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_function_u(0:n,0:n),flow_function_v(0:n,0:n) u=0 v=0 dod

30、o i=1,n-1 u(i,n+1)=-16.0*(real(i)*h)*2*(1.0-(real(i)*h)*2)*2-u(i,n) end doend do p=0 p_wave =0 d=1.0 dodo whilewhile(d=standard) ifif(countercounter_max) exitexit !如果时间推进步数大于五万步，跳出循环 u_0=u callcall sub1_x_momentum(u,v,p,h,t,n,re) callcall sub2_y_momentum(u_0,v,p,h,t,n,re) callcall sub3_pressure(u,v,

31、p,p_wave,n,t,re) dodo i=1,n dodo j=1,n d_array(i,j)= 10.0*abs(u(i,j)-u(i-1,j)+v(i,j)-v(i,j-1) end doend do end doend do d=maxval(d_array) !d作为收敛判据 dodo i=1,n dodo j=1,n d_array(i,j)=abs(p_wave(i,j) !只考虑实际的压力波动，即边界内的压力波动 end doend do end doend do p_wave_max=maxval(d_array) counter=counter+1 writewrit

32、e(*,*)d=,d,p_wave_max=,p_wave_max,counter-1 end doend do writewrite(*,*)counter-1 flow_function_u = 0 flow_function_v = 0 dodo i=1,n-1 dodo j=1,n flow_function_u(i,j)= h*u(i,j)+flow_function_u(i,j-1) end doend do end doend do dodo j=1,n-1 dodo i=1,n flow_function_v(i,j)= -1.0*h*v(i,j)+flow_function_

33、v(i-1,j) end doend do end doend do dodo i=0,n dodo j=0,n flow_function(i,j)= 0.5*(flow_function_u(i,j)+flow_function_v(i,j) end doend do end doend do dodo i=1,n-1 dodo 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 doend do end doend do dodo j=1,n-1 vortex_function(0,j)= (

34、v(1,j)-v(0,j)/h - (u(0,j+1)-u(0,j)/h end doend do dodo j=1,n-1 vortex_function(n,j)= (v(n+1,j)-v(n,j)/h - (u(n,j+1)-u(n,j)/h end doend do dodo i=1,n-1 vortex_function(i,n)= (v(i+1,n)-v(i,n)/h - (u(i,n+1)-u(i,n)/h end doend do dodo i=1,n-1 vortex_function(i,0)= (v(i+1,0)-v(i,0)/h - (u(i,1)-u(i,0)/h e

35、nd doend do 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

36、,n) dodo j=0,n dodo i=0,n g(i)=(u(i,j+1)+u(i,j)*0.5 endend dodo dodo i=0,n u(i,j)=g(i) end doend do end doend do dodo i=0,n dodo j=0,n g(j)=(v(i+1,j)+v(i,j)*0.5 end doend do dodo j=0,n v(i,j)=g(j) end doend do end doend do dodo i=0,n dodo j=0,n p_wave(i,j)=0.25*(p(i,j)+p(i+1,j)+p(i,j+1)+p(i+1,j+1) e

37、nd doend do end doend do p=p_wave openopen(unit=10,file=data.dat) writewrite(10,(title = grid=,i3,re=,f5.1,)n,re !输出双引号时需要打出两个双引号，否则出错 writewrite(10,*)variables = x y u v u-scalar p flow-f vortex-f !这种情况不需要打两个双引号就可以输出双引号 writewrite(10,(/) writewrite(10,(zone n=,i6, e=,i6, zonetype=fequadrilateral, d

38、atapacking=point)(n+1)*(n+1),n*n dodo j=0,n dodo i=0,n writewrite(10,(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 do end doend do dodo j=0,n-1 !对网格节点进行编号，便于tecplot软件后处理 i= j*(n+1) cou

39、nter=1 dodo whilewhile (counter=n) writewrite(10,(i5,i5,i5,i5)i+counter, i+counter+1, i+counter+1+n+1, i+counter+n+1 counter=counter+1 end doend do end doend do openopen(unit=11,file=function_max.txt) d=0.0 dodo j=0,n dodo i=0,n ifif (d=flow_function(i,j) thenthen d=flow_function(i,j) i_max=i j_max=

40、j end ifend if end doend do end doend do writewrite(11,*)flow_functioin_max: writewrite(11,(f15.10,i5,i5,f15.10,f15.10)d,i_max,j_max,real(i_max)*h,real(j_max)*h d=0.0 dodo i=0,n ifif (d=vortex_function(i,n) thenthen d=vortex_function(i,n) i_max=i end ifend if end doend do writewrite(11,*)up_wall_vor

41、tex_functioin_max: writewrite(11,(f15.10,i5,f15.10)d,i_max,real(i_max)*h d=0.0 dodo j=0,n ifif (d= standard) error = 0.0 dodo j=1,n !迭代法计算修正压力 dodo i=1,n a_u1= 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/t a_u2= h*(0.25*(u(i,j)-u(i-2,j)+2.0/(

42、h*re) + h*(0.25*(v(i-1,j)+v(i,j)-v(i-1,j-1)-v(i,j-1)+2.0/(h*re) + h*h/t a_v1= h*(0.25*(v(i,j+1)-v(i,j-1)+2.0/(h*re) + h*(0.25*(u(i,j+1)+u(i,j)-u(i-1,j+1)-u(i-1,j)+2.0/(h*re) + h*h/t a_v2= h*(0.25*(v(i,j)-v(i,j-2)+2.0/(h*re) + h*(0.25*(u(i,j)+u(i,j-1)-u(i-1,j)-u(i-1,j-1)+2.0/(h*re) + h*h/t a=h*h/a_u1 + h*h/a_u2 + h*h/a_v1 + h*h/a_v2 b=h*h/a_