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

SIMPLE算法求解方腔内粘性不可压流动一、问题描述假设1,0yx的方腔内充满粘性不可压缩流体，左、右、下壁固定，上壁以22116xxu运动，试求400,200,100Re时的定常解，方腔如图1所示。图1方腔内流动示意图二、离散格式本算例采用求解不可压缩流动的经典算法，即SIMPLE算法，求解方腔内粘性不可压缩流体运动的定常解。SIMPLE算法的全称为Semi-ImplicitMethodforPressure-LinkedEquations，即求解压力关联方程的半隐式算法。采用SIMPLE算法时，为了避免中心差分格式将“棋盘”型参量分布误认为是均匀分布，需要用交错网格对计算域进行离散。交错网格交错网格如图2所示，压力、密度等物理量存储在控制体ji,的中心，这个控制体称为主控制体。速度分量vu,分别存储在主控制体的ji,2/1和2/1,ji位置处，标记为ji,位置，再分别以此为中心，划分速度分量u、v的控制体。采用空间均匀网格，等间距离散整个求解域，如图3所示。图2交错网格示意图图3求解域离散示意图图3中阴影部分代表方腔内的流动区域，阴影区域的边界代表方腔的上、下、左、右壁面，阴影区域外面的网格节点是为边界处理需要而设定的虚拟网格节点，后面介绍边界处理方法时详细论述。方程离散无量纲化的守恒型不可压缩SN方程为0Re102UPUUtUU其积分形式为SSySVSSxSVSvdSndSpnvdSUndVtvudSndSpnudSUndVtudSUn0Re10Re10图4主控制体图5速度u控制体图6速度v控制体采用有限体积法离散SN方程，连续性方程在主控制体上离散011,1,1,11,xvvyuuMjiMjiMjiMjiX方向动量方程在速度u控制体上离散，时间采用前差01,1,111,1,1,11,,1,yppxGGyFFuutyxMjiMjijijijijiMjiMjiY方向动量方程在速度v控制体上离散，时间采用前差01,11,21,2,2,12,,1,xppyGGxFFvvtyxMjiMjijijijijiMjiMji其中，数值通量yuuvGxuuFRe1,Re1121yvuvGxvvFRe1,Re1222通量11,GF分别定义在主控制体的中心和角点，如图所示，并按照如下方法离散1,1,11,1,1,,111,1,11,1,1,,11Re141Re141MjiMjiMjiMjiMjiMjiMjiMjiMjiMjiMjiMjiuuyuuvvGuuxuuuuF通量22,GF分别定义在主控制体的中心和角点，如图所示，并按照如下方式离散1,1,11,1,1,,121,1,11,1,1,,12Re141Re141MjiMjiMjiMjiMjiMjiMjiMjiMjiMjiMjiMjivvxvvuuGvvyvvvvF通量11,GF和22,GF的某些项冻结于M时间层，使离散化之后的方程对11,MMvu是线性的。将离散化之后的11,GF和22,GF代入离散后的x方向和y方向的动量方程，整理之后得离散后的动量方程如下001,11,,1,,1,,1,1,1,1,,1,,xppbvavayppbuauaMjiMjivjiMqpvqpMjivjiMjiMjiujiMqpuqpMjiuji其中Mjiujiutyxb,,yvvxayvvxaxuuyaxuuyaMjiMjiujiMjiMjiujiMjiMjiujiMjiMjiujiRe141,Re141Re141,Re1411,11,1,,1,1,,1-,,1-,,1,1tyxyvvvvxxuuyaMjiMjiMjiMjiMjiMjiujiRe241Re2411,11,,1,,1,1,Mjivjivtyxb,,xuuyaxuuyayvvxayvvxaMjiMjivjiMjiMjivjiMjiMjivjiMjiMjivjiRe141,Re141Re141,Re141,11,1,1,1,,1,1,1,1,,1,tyxxuuuuyyuvxaMjiMjiMjiMjiMjiMjivjiRe241Re241,11,1,1,1,1,,以上是SIMPLE算法中离散化的动量方程三、SIMPLE算法基本思想SIMPLE算法是一种解决压力-速度耦合问题的“半隐式”算法。首先给定M时刻猜测的速度场MMvu,，用于计算离散动量方程中的系数和常数项。给定M+1时刻猜测的压力场估计值*p，迭代求解离散动量方程，得到M+1时刻速度场的估计值**,vu，速度场的估计值**,vu满足如下离散方程。00*,*1,,*,,*,,*,*,1,*,,*,,xppbvavayppbuauajijivjiqpvqpjivjijijiujiqpuqpjiuji一般地，速度场**,vu不满足离散的连续性方程，因而需要对速度场**,vu和压力场*p进行修正。M+1时刻的修正值和估计值有如下关系pppvvvuuuMMM*1*1*1其中，vu,和p分别速度和压力的修正量，修正量亦满足离散的动量方程00,1,,,,,,,,1,,,,,xppbvavayppbuauajijivjiqpvqpjivjijijiujiqpuqpjiuji编号为（i,j）的速度修正量vu,不仅与压力修正量p有关，还与邻近点的速度修正量有关。SIMPLE算法的重要假定：速度的改变只与压力的改变有关，忽略邻近点对速度修正的影响。因而得到如下速度修正量jijivjijijijiujijippaxvppayu,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时刻的速度场11,MMvu和压力场1Mp。将M+1时刻的速度场11,MMvu和压力场1Mp作为新的猜测的速度场和猜测的压力场估计值，采用上述方法计算下一个时刻的速度场和压力场，直到满足收敛条件。收敛判据pjibMax,为很小的正实数，视计算的精度要求而定。本算例中取81e。若0,pjib，则0,jip，此时*,1,jiMjiUU，从而来自于离散动量的*,jiU满足离散的连续性方程。因此pjibMax,可以作为收敛判据。边界条件处理首先对计算区域离散，并流动边界之外扩充一个虚拟网格，将真实流动的离散域包围，如图7所示。图7虚拟网格扩充示意图压力p的存储位置由图中符号代表，速度分量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,,,niiniipppp。虚拟网格处理静止壁面的虚拟网格按图8所示方式处理，可以满足粘性边界条件，以下边界为例说明。图8静止壁面虚拟网格处理示意图边界外的虚拟网格存储位置上的速度分量，由边界内的真实网格对应存储位置上的速度分量沿边界对称得到。这样可以保证虚拟网格与真实网格在对应存储位置上的速度分量大小相等，方向相反，二者在边界上的线性平均值为零，满足粘性边界条件。上壁面有速度，为运动壁面，速度分量0,0vu。速度分量v采用前述对称方法处理。速度分量u的处理方法略微不同。上壁面速度boundaryupu不随时间变化，虚拟网格上速度分量1,,niu由1,,niboundaryupuu线性插值得到，可以保证上壁面满足边界条件，即1,1,,2niboundaryupniuuu。存储位置位于恰好位于静止壁面上的速度分量恒为零，不必迭代计算，例如0,0ju。边界外虚拟网格中心的压力和压力修正量，认为分别等于边界附近对应真实网格中心的压力和压力修正量，即假定区域边界附近沿边界法向的压力梯度为零。因此，只需求解真实网格上各个存储位置上变量的值，虚拟网格上各个存储位置上变量的值可以通过前述处理虚拟网格的方式得到。至此，包括虚拟网格在内的所有存储位置上都被赋予了确定的值。方程求解X方向动量方程迭代法求解x方向动量方程，若要得到M+1时刻*,jiu的值，需要用到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两个位置的压力值，包括虚拟网格在内的离散域满足计算需求。同理可得，在其他边界位置，离散域满足求解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,pp两个位置的压力值，包括虚拟网格在内的离散域满足计算需求