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

SIMPLE算法求解方腔内粘性不可压流动目录一、问题描述...........................................................................2二、离散格式...........................................................................3交错网格..........................................................................3方程离散..........................................................................4三、SIMPLE算法基本思想......................................................7边界条件处理..................................................................8虚拟网格处理..................................................................9方程求解........................................................................11输出变量处理................................................................12SIMPLE算法流程图.......................................................15四、程序中主要变量的意义.................................................15五、计算结果与讨论.............................................................17函数最大值....................................................................17变量等值线图................................................................18主要结论........................................................................22六、源程序.............................................................................22一、问题描述假设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方向虚