SIMPLE算法求解方腔内粘性不可压流动学习资料

1、精选优质文档-倾情为你奉上SIMPLE算法求解方腔内粘性不可压流动目录一、问题描述 假设 的方腔内充满粘性不可压缩流体，左、右、下壁固定，上壁以 运动，试求 时的定常解，方腔如图1所示。图1 方腔内流动示意图二、离散格式 本算例采用求解不可压缩流动的经典算法，即SIMPLE算法，求解方腔内粘性不可压缩流体运动的定常解。SIMPLE算法的全称为Semi-Implicit Method for Pressure-Linked Equations，即求解压力关联方程的半隐式算法。采用SIMPLE算法时，为了避免中心差分格式将“棋盘”型参量分布误认为是均匀分布，需要用交错网格对计算域进行离散。交错网格

2、交错网格如图2所示，压力、密度等物理量存储在控制体的中心，这个控制体称为主控制体。速度分量分别存储在主控制体的 和 位置处，标记为 位置，再分别以此为中心，划分速度分量u、v的控制体。采用空间均匀网格，等间距离散整个求解域，如图3所示。图2 交错网格示意图图3 求解域离散示意图图3中阴影部分代表方腔内的流动区域，阴影区域的边界代表方腔的上、下、左、右壁面，阴影区域外面的网格节点是为边界处理需要而设定的虚拟网格节点，后面介绍边界处理方法时详细论述。方程离散无量纲化的守恒型不可压缩 方程为其积分形式为 图4 主控制体 图5 速度u控制体 图6 速度v控制体采用有限体积法离散 方程，连续性方程在主控

3、制体上离散X方向动量方程在速度u控制体上离散，时间采用前差Y方向动量方程在速度v控制体上离散，时间采用前差其中，数值通量通量 分别定义在主控制体的中心和角点，如图所示，并按照如下方法离散通量 分别定义在主控制体的中心和角点，如图所示，并按照如下方式离散通量和的某些项冻结于M时间层，使离散化之后的方程对 是线性的。将离散化之后的和代入离散后的x方向和y方向的动量方程，整理之后得离散后的动量方程如下 其中 以上是SIMPLE算法中离散化的动量方程三、SIMPLE算法基本思想SIMPLE算法是一种解决压力-速度耦合问题的“半隐式”算法。首先给定M时刻猜测的速度场，用于计算离散动量方程中的系数和常数项

4、。给定M+1时刻猜测的压力场估计值，迭代求解离散动量方程，得到M+1时刻速度场的估计值，速度场的估计值满足如下离散方程。 一般地，速度场不满足离散的连续性方程，因而需要对速度场和压力场进行修正。M+1时刻的修正值和估计值有如下关系 其中，和分别速度和压力的修正量，修正量亦满足离散的动量方程 编号为（i,j）的速度修正量不仅与压力修正量有关，还与邻近点的速度修正量有关。SIMPLE算法的重要假定：速度的改变只与压力的改变有关，忽略邻近点对速度修正的影响。因而得到如下速度修正量 修正后的速度分量 将修正后的速度分量代入离散后的连续性方程，得到压力修正方程 其中 采用迭代法求解压力修正方程，得到压力

5、修正量，代入修正公式得到M+1时刻的速度场和压力场。将M+1时刻的速度场和压力场作为新的猜测的速度场和猜测的压力场估计值，采用上述方法计算下一个时刻的速度场和压力场，直到满足收敛条件。收敛判据 为很小的正实数，视计算的精度要求而定。本算例中取。若，则，此时，从而来自于离散动量的满足离散的连续性方程。因此 可以作为收敛判据。边界条件处理首先对计算区域离散，并流动边界之外扩充一个虚拟网格，将真实流动的离散域包围，如图7 所示。图7 虚拟网格扩充示意图压力p的存储位置由图中符号代表，速度分量u的存储位置由图中符号代表，速度分量v的存储位置由图中符号 代表。阴影区域表示真实的流动区域，阴影区域外的网格

6、为虚拟网格。速度分量u在x方向虚拟网格上的存储位置分别记为 ，在y方向虚拟网格上的存储位置分别记为。速度分量v在x方向虚拟网格上的存储位置分别记为 ，在y方向虚拟网格上的存储位置分别记为。压力和压力修正量在x方向虚拟网格上的存储位置分别记为，在y方向虚拟网格上的存储位置分别记为。虚拟网格处理静止壁面的虚拟网格按图8所示方式处理，可以满足粘性边界条件，以下边界为例说明。图8 静止壁面虚拟网格处理示意图边界外的虚拟网格存储位置上的速度分量，由边界内的真实网格对应存储位置上的速度分量沿边界对称得到。这样可以保证虚拟网格与真实网格在对应存储位置上的速度分量大小相等，方向相反，二者在边界上的线性平均值为

7、零，满足粘性边界条件。上壁面有速度，为运动壁面，速度分量。速度分量v采用前述对称方法处理。速度分量u的处理方法略微不同。上壁面速度不随时间变化，虚拟网格上速度分量由线性插值得到，可以保证上壁面满足边界条件，即。存储位置位于恰好位于静止壁面上的速度分量恒为零，不必迭代计算，例如。边界外虚拟网格中心的压力和压力修正量，认为分别等于边界附近对应真实网格中心的压力和压力修正量，即假定区域边界附近沿边界法向的压力梯度为零。因此，只需求解真实网格上各个存储位置上变量的值，虚拟网格上各个存储位置上变量的值可以通过前述处理虚拟网格的方式得到。至此，包括虚拟网格在内的所有存储位置上都被赋予了确定的值。方程求解X

8、方向动量方程迭代法求解x方向动量方程，若要得到M+1时刻的值，需要用到M时刻，九个位置的速度分量值和M+1时刻两个位置的压力值。内点上的值，可以直接迭代得到。求M+1时刻速度分量，需要用到M时刻， 九个位置的速度分量值和M+1时刻 两个位置的压力值，包括虚拟网格在内的离散域满足计算需求。同理可得，在其他边界位置，离散域满足求解x方向动量方程的需要。Y方向动量方程迭代法求解y方向动量方程，若要得到M+1时刻的值，需要用到M时刻，九个位置的速度分量值和M+1时刻两个位置的压力值。内点上的值，可以直接迭代得到。求M+1时刻速度分量，需要用到M时刻， 九个位置的速度分量值和M+1时刻 两个位置的压力值

9、，包括虚拟网格在内的离散域满足计算需求。同理可得，在其他边界位置，离散域满足求解y方向动量方程的需要。压力修正方程迭代法求解压力修正方程，若要M+1时刻得到 的值，需要用到M时刻四个位置的压力修正值和M+1时刻 ，十六个位置的速度分量估计值。内点上的值，可以直接迭代得到。求M+1时刻压力修正值，需要用到M时刻 四个位置的压力修正值和M+1时刻十六个位置的速度分量估计值，包括虚拟网格在内的离散域满足计算需求。同理可得，在其他边界位置，离散域满足求解压力修正方程的要求。将M+1时刻变量的估计值和修正值代入如下公式即可得到M+1时刻变量的值。输出变量处理涡量定义 涡量是一个矢量，在平面流动中只有一个

10、方向，垂直于流动平面。流函数满足流函数等值线为流线在交错网格里，对涡量和流函数在一个网格内做中心差分，涡量和流函数的存储位置为真实网格的节点，如图9所示。图9 涡量和流函数存储位置示意图涡量差分格式y方向流函数差分格式 x方向流函数差分格式 任取其中一个式子计算流函数即可，或者将两个式子的计算结果取平均值。本算例采用的是取平均值的办法。边界上的流函数值为零。需要输出用于后处理的变量主要有位置坐标x、y，速度分量u、v，压力p，流函数和涡量。由于SIMPLE算法采用交错网格，不同性质的变量存储在网格的不同位置，在变量输出时，有必要将所有变量统一到真实网格的网格节点上，便于后处理，如图10所示。图

11、10 变量统一后存储位置示意图网格节点上的速度分量为节点上下两速度分量的平均值，即网格节点上的速度分量为节点左右两速度分量的平均值，即网格节点上的压力为节点周围四个主控制体中心压力的平均值，即边界处的压力平均采用相同的公式，需要用到虚拟网格中心的压力，角点依然。SIMPLE算法流程图图11 SIMPLE算法流程图四、程序中主要变量的意义主程序中主要变量的意义变量意义counter_max预设最大循环次数grid_number、nX方向和y方向的离散网格数目t、h时间步长和空间步长Re雷诺数d收敛判据，压力修正方程中绝对值的最大值数组u、vX方向和y方向速度分量数组p、p_wave压力和压力修正

12、量数组flow_function、vortex_function流函数和涡量函数数组u_scalar速度绝对值数组d_array、g、u_0计算中需要用到的一些辅助数组子程序1，即求解x方向动量方程程序中主要变量的意义变量意义nX方向和y方向的离散网格数目t、h时间步长和空间步长Re雷诺数a、b、c、d、e迭代求解x方向动量方程所需的辅助变量数组u、v、pX方向、y方向速度分量和压力场子程序2，即求解y方向动量方程程序中主要变量的意义变量意义nX方向和y方向的离散网格数目t、h时间步长和空间步长Re雷诺数a、b、c、d、e迭代求解x方向动量方程所需的辅助变量数组u、v、pX方向、y方向速度分量

13、和压力场 子程序3，即求解压力修正方程程序中主要变量的意义变量意义nX方向和y方向的离散网格数目t、h时间步长和空间步长Re雷诺数a、b、c、e、f压力修正方程中各压力修正项系数d压力修正方程中项a_u1、a_u2辅助变量，与u方向速度修正式有关a_v1、a_v2辅助变量，与v方向速度修正式有关error求解压力修正方程的迭代收敛判据数组u、vX方向、y方向速度分量数组p、p_wave压力和压力修正量五、计算结果与讨论本算例分别对网格数目为2020,4040两种情况，数值求解了Re=100,200,400时的方腔流动的定常解，计算结果列举如下。函数最大值离散网格为2020数值计算得到的流函数最

14、大值Re100200400流函数最大值0.0.0.对应x坐标0.400.450.45对应y坐标0.650.600.60离散网格为2020数值计算得到的上壁面涡函数最大值 Re100200400上壁面涡函数最大值91.104.117.对应x坐标0.800.750.75对应y坐标1.001.001.00离散网格为2020数值计算得到的方腔中线涡函数最大值 Re100200400中线上涡函数最大值25.38.53.对应x坐标0.500.500.50对应y坐标0.951.001.00离散网格为4040数值计算得到的流函数最大值Re100200400流函数最大值0.0.0.对应x坐标0.4250.450

15、0.475对应y坐标0.6250.6000.575离散网格为4040数值计算得到的上壁面涡函数最大值 Re100200400上壁面涡函数最大值119.149.181.对应x坐标0.8250.8000.800对应y坐标1.001.001.00离散网格为4040数值计算得到的方腔中线涡函数最大值 Re100200400中线上涡函数最大值32.41.45.对应x坐标0.500.500.50对应y坐标0.9750.9750.975从表中数据可以看到，随着雷诺数增大，流函数最大值减小，最大流函数的位置也更趋于方腔的中心。这是因为雷诺数增大，粘性效应减弱，上壁面运动带动放腔内流体运动的能力降低。上壁面涡量

16、最大值位于上壁面的速度最大处，这是因为在该处速度分量u沿y方向的的梯度最大，速度分量v沿x方向的梯度为零（壁面粘性边界条件），因而剪切效应最强。随着雷诺数增大，上壁面涡量最大值增大。中线上涡量最大值出现在上壁面附近，但不位于上壁面。离散网格为2020，雷诺数为100时，数值计算结果表明，中线上最大涡量出现在y=0.95处；雷诺数为200和400时，数值计算结果表明，中线上最大涡量出现在y=1.00处。离散网格为4040，对于雷诺数为100、200和400三种情况，数值计算结果表明，中线上最大涡量出现在y=0.975处。这说明增大离散网格数目，即加密网格，可以增加数值计算对涡量最大值的捕捉精度。

17、变量等值线图速度大小等值线图12 Re=100时速度大小等值线图13 Re=200时速度大小等值线图14 Re=400时速度大小等值线流函数等值线图15 Re=100时流函数等值线图16 Re=200时流函数等值线图17 Re=400时流函数等值线涡量等值线图18 Re=100时涡量等值线 图19 Re=200时涡量等值线 图20 Re=400时涡量等值线从以上各个变量的等值线图中可明显看出，相同雷诺数不同网格密度时，各个变量相应的等值线形状基本相同，但4040网格的等值线明显比2020网格的等值线光滑。这说明加密网格可以提高数值计算的精度。从流函数等值线中可以看到，Re=100和200时，方

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

19、为负，涡量的最大值出现在速度分量u沿y方向梯度与速度分量v沿x方向梯度只差最大的地方，因而不一定出现在壁面上。六、源程序program main implicit none integer,parameter : counter_max=50000 !设置最大的时间推进步数为五万步real(8),parameter : standard=0. !小数点后取八位integer : i,j,i_max,j_max,n,counter,grid_number !grid_numer为网格数目real(8),allocatable : u(:,:),v(:,:),p(:,:),p_wave(:,:),

20、d_array(:,:),g(:),u_0(:,:),flow_function(:,:),flow_function_u(:,:),flow_function_v(:,:),vortex_function(:,:),u_scalar(:,:)real(8) : h,t,Re !h为空间步长，t为时间步长，Re为Reynolds numberreal(8) : d,p_wave_maxwrite(*,*)输入离散网格数目：read(*,*)grid_numberwrite(*,*)输入Reynolds number：read(*,*)Re counter=1 !迭代步数计数器 h=1.0/re

21、al(grid_number) n=grid_number t=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_function_u(0:n,0:n),flow_function_v(0:n,0:n) !声明二维数组的长度,p_wave代表修正压力 u=0 !赋初值

22、 v=0 do 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 do p=0 p_wave =0 d=1.0 do while(d=standard) if(countercounter_max) exit !如果时间推进步数大于五万步，跳出循环 u_0=u call sub1_x_momentum(u,v,p,h,t,n,Re) !求解离散的x方向动量方程 call sub2_y_momentum(u_0,v,p,h,t,n,Re) !求解离散的y方向动量方程 call sub3_pressure(u

23、,v,p,p_wave,n,t,Re) !求解压力修正方程 do i=1,n do 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 do end do d=maxval(d_array) !d作为收敛判据 do i=1,n do j=1,n d_array(i,j)=abs(p_wave(i,j) !只考虑实际的压力波动，即边界内的压力波动 end do end do p_wave_max=maxval(d_array) counter=counter+1 write(*,*)d=,d,p_wave_max=,p

24、_wave_max,counter-1 end do write(*,*)counter-1 flow_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*v(i,j)+flow_function_v(i-1,j) end do end do do i=0,n do j=0,n flow_fun

25、ction(i,j)= 0.5*(flow_function_u(i,j)+flow_function_v(i,j) end do end do do i=1,n-1 !求涡量函数 do 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)= (v(1,j)-v(0,j)/h - (u(0,j+1)-u(0,j)/h !left wall end do do j=1,n-1 vortex_func

26、tion(n,j)= (v(n+1,j)-v(n,j)/h - (u(n,j+1)-u(n,j)/h !right wall end do do i=1,n-1 vortex_function(i,n)= (v(i+1,n)-v(i,n)/h - (u(i,n+1)-u(i,n)/h !up wall end do do i=1,n-1 vortex_function(i,0)= (v(i+1,0)-v(i,0)/h - (u(i,1)-u(i,0)/h !down wall end do !计算角点涡量函数 vortex_function(0,0)= 0.5*(vortex_function

27、(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 !将u 速度节点与网格节点重合 do i=0,n g(i)=(u(i,j+1)+u

28、(i,j)*0.5 end do do i=0,n u(i,j)=g(i) end do end do do i=0,n !将v 速度节点与网格节点重合 do j=0,n g(j)=(v(i+1,j)+v(i,j)*0.5 end do do j=0,n v(i,j)=g(j) end do end do do j=0,n !求节点速度绝对值 do i=0,n u_scalar(i,j)= sqrt(u(i,j)*u(i,j) + v(i,j)*v(i,j) end do end do do i=0,n !将压力节点与网格节点重合 do j=0,n p_wave(i,j)=0.25*(p(i,

29、j)+p(i+1,j)+p(i,j+1)+p(i+1,j+1) end do end do p=p_wave open(unit=10,file=data.dat) write(10,(TITLE = Grid=,I3,Re=,f5.1,)n,Re !输出双引号时需要打出两个双引号，否则出错 write(10,*)VARIABLES = X Y U V U-SCALAR P FLOW-F VORTEX-F !这种情况不需要打两个双引号就可以输出双引号 write(10,(/) write(10,(ZONE N=,I6, E=,I6, ZONETYPE=FEQuadrilateral, DATA

30、PACKING=POINT)(n+1)*(n+1),n*n do j=0,n do i=0,n write(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 do end do do j=0,n-1 !对网格节点进行编号，便于tecplot软件后处理 i= j*(n+1) counter=1 do while (counter=n

31、) write(10,(I5,I5,I5,I5)i+counter, i+counter+1, i+counter+1+n+1, i+counter+n+1 counter=counter+1 end do end do open(unit=11,file=function_max.txt) d=0.0 do j=0,n !流函数最大值 do i=0,n if (d=flow_function(i,j) then d=flow_function(i,j) i_max=i j_max=j end if end do end do write(11,*)flow_functioin_max: wr

32、ite(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 !上壁面涡函数最大值 do i=0,n if (d=vortex_function(i,n) then d=vortex_function(i,n) i_max=i end if end do write(11,*)up_wall_vortex_functioin_max: write(11,(f15.10,I5,f15.10)d,i_max,real(i_max)*h d=0.0 !中竖线上涡函数最大值 do j=0,n if

33、(d= standard) error = 0.0 do j=1,n !迭代法计算修正压力 do 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/(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

34、,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_u2 c=h*h/a_u1 d=(u(i,j)-u(i-1,j)*h + (v(i,j)-v(i,j-1)*h e=h*h/a_v2 f=h*h/a_v1 p_wave_0 = p_wave(i,j) p_wave(i,j)=(b*p_wave(i-1,j) + c*p_wave(i+1,j) + e*