有限体积法介绍

有限体积法1有限体积法基本原理上一章讲到的有限差分法将数值网格的节点上定义为计算节点，并在网格节点上对微分形式的流体基本方程进行离散，用网格节点上的物理量的代数方程作为原PDE的近似。在本章所要学习的有限体积法则采用了不同的离散形式。首先，有限体积法离散的是积分形式的流体力学基本方程：jp伽-nds=j「n-V^ds+jq^dQ*计算域用数值网格划分成若干小控制体。和有限差分法不同的是，有限体积法的网格定义了控制体的边界，而不是计算节点。有限体积法的计算节点定义在小控制体内部。一般有限体积法的计算节点有两种定义方法，一种是将网格节点定义在控制体的中心，另一种方法中，相邻两个控制体的计算节点到公共边界的距离相等。第一种方法的优点在于用计算节点的值作为控制体上物理量的平均值具有二阶的精度;第二种方法的好处是在控制体边界上的中心差分格式具有较高的精度。2面积分的近似采用结构化网格，在二维情况下，每一个控制体有4个面，二维情况，每一个控制体有6个表面。计算节点用大写字母表示，控制体边界和节点用小写字母表示。为了保证守恒性，控制体不能重叠，每一个面都是相邻两个控制体的唯一公共边界。控制体边界上的积分等于控制体个表面的积分的和：(2)jfdS=Ejfds伽上式中，f可以表示pnu或r—-ndn(2)0-NN00ONWnwN°nneNE口邙wWWwp0e一seeqEE。0sWS。SE°0&SS°0■0显然，为了获得边界上的积分，必须知道f在边界上的详细分布情况，这是不可能实现的，由于只是计算节点上的函数值，因此必须采用近似的方法来计算积分。整个近似过程分成两步第一步：用边界上几个点的近似积分公式第二步：边界点上的函数值用计算节点函数值的插值函数近似面积分可采用以下不同精度的积分公式：二阶精度积分：F=jfds=fS注fSeeeeeSe上式中fe为边界中点出的函数值。近似为方格中心点的值乘以方格的面积。三阶精度积分：F=jfdsrf+fSeS2ee四阶精度积分：F=jfdsr九+4f+九SeS6ee应该注意的是，采用不同精度的积分公式，在相应的边界点的插值时也应采用相应精度的插值函数。积分公式的精度越高，近似公式就越复杂。3体积分的近似和面积分相似，体积分也有不同精度的近似公式二阶精度积分公式(6)Q=jqds=qS牝qAQSe采用双二次样条函数q(x,y)=a+ax+ay+ax2+ay2+axy+ax2y+axy2+ax2y2012345678可以得到四阶精度的积分公式：(6)Q=jqdsR竺^(16q+4q+4q+4q+4q+4q+4q+4q+4q)(8)S36Pswnsseswnen^w、e4函数的插值在上节讲到的积分的近似公式中用到了非计算节点上的函数值，被积函数f中包含了多个物理量及其偏微分，如对流项fc=pnV-n，扩散项fd=「n・V4，在源项中也有类似情况，这里假定流场和流体的物性参数是已知的，物理量4及其偏导数在控制面上的值需要通过计算节点上物理量的插值得到。下面已e面为例进行讨论。4.1迎风插值(UDS)七用上游计算节点的函数值近似相当于对一阶偏导数采用迎风格式，因此用UDS来表示这种近似方法，在UDS中：「8if(v・n)>0f：if(v•n)：<0⑼UDS是唯一无条件满足有界性要求的近似格式，在数值过程中不会产生数值振荡。UDS存在数值粘性。根据Taylor公式，该格式具有一阶精度，并具有数值粘性：Vnum=(pU)AX/2(10)在多维问题中，如果流动方向和网格是斜交的，截断误差会在垂直于流动方向以及流线方向产生扩散，这是一种非常严重的误差，函数的峰值或函数值的快速变化会被抹平，为了得到高精度结果需要采用非常精细的网格。4.2线性插值(CDS)(11)(12)线性插值具有二阶精度，线性插值相当于FDM中的CDS格式，因此用CDS表示。CDS格式会产生数值振荡。对于扩散项(13)

4.3三阶迎风格式（QUICK）和UDS类似，QUICK格式也和流动方向有关g8—g8+(1-g+g)8if(v・n)>01E2W12P/eg8—g8+(1—g+g)8if(v-n)<0*34Ee3P4EE其中:G—xlu1+Xe,pL(+xX—xe,W<1+Xe,E-"e,PP——g，g45",P(14)(15a)(15b)4.4高阶格式（4阶精度CDS）采用三次曲线可拟合出四阶精度的中心插值公式，在均匀网格中，四阶公式为:278〃+27虬-3虹-38.―P"48~~WEE(16)%8)=278-278p+8W-8侦xJ24Axe(17)5边界的处理对于对流项，在入口处一般给出了流量或函数值，在边界和对称面上流量为零，在出口处假设和出口的法向坐标无关，因此可采用迎风格式。对于扩散项则可能需要采用偏心格式。6有限体积法应用举例例：考虑一标量在已知流场中的输运过程（如图4.4所示），输运方程为:jp8v-ndS=jrv8-ndSSS(18)边界条件:8=0；北部入口边界8=1—y；西部壁面边界对称条件；南部边界梯度为0；东部出口条件u=-y，流线方程xy=c对流项:Fc=jp^v-ndS浇成8=jpv-ndS=(pu)Ay为质量通量。Se(19)forUDSforCDS(20)Ac-=min(m,0)；Ac=min(m,0)EeforUDSforCDS(20)Ac-=min(m,0)；Ac=min(m,0)EeWwAc-=min(m,0)；Ac=min(m,0)NnSsAc-p=—(Ac+Ac+AcN+Ac)S若采用UDS格式，代数方程组中各项系数为:若采用CDS格式，代数方程组中各项系数为(21)e[mt(1—人)8+mx8AcE=z&X；AcW=mxAcN=戒x；AcS=戒xssAcP=—(AcE+AcW+ANc+ASc)(22)根据连续性方程:(22)(23)相邻CV之间的关系:(24)生,p"畔,"七,P=1-Xe,W其余相邻CV有类似关系(24)扩散项采用CDS格式Fc=jrV^-ndS〜eSe'「纠AyIox)~rAy〜xE(25)代数方程组中扩散项系数为:AdE「Ayx—x；AdWAdN「Ax；AdWrAx(26)AdP=—(Ad+Ad—WE对于任意控制体+AdN+Ad)S=Q(27)A=Ac+Ad边界条件的处理：对于西部和北部边界，由于给定了函数值，对流项可直接代入函数值而无需插值，扩散项则采用一侧差分，l为任意指标P，E，W，S，N。(28)©项—PW-xP-xWIEW这里，W点和P的w边中点重合。(29)南边和西边的梯度为零，以南边为例，由于梯度为零，，代数方程变为:A。+(A+A)。+A。+A。=QWWSPPNNEEP(30)6SIMPLE方法考虑定常不可压流动问题，控制方程为：连续性方程：jpv-ndS=0SV动量方程：jpvv-ndS=jpn-VvdS—jpndS+jpbdQSVSVSVCV(31)(32)不可压缩问题求解的困难在于压力场的求解。主要原因在于压力p没有独立的方程组。先考虑一维问题：对于动量方程：(puu)-(puu)(8")I8x)w(33)若采用CDS格式(P"")+(pu")(P"")+(P"")^2^u-uu-up(P"")+(pu")(P"")+(P"")^2^AxAx22简化后得:(puu)(puu)—(puu)^2Wu+u—2u

f*AXP—pe+pw-E2W(34)根据连续性方程，u=u=u=c，则有p=p，由于相邻节点之间的压力没有联i+1ii—1i—1i+1系方程，容易造成压力交错现象。为了解决这一问题，可采用交错网格技术，即速度场和压力场采用不同的网格。以二维问题为例，交错网格的布置如下图所示：1~r~_1_■1.._.1_.1IN11fv„1iTn1%i1IW1i1lviTs11i1_1_.i..—.1—.1IS11i1主控制体为压力控制体（黑色实线网格），u的控制体（红色虚线网格）的计算节点在主控制体的e边，控制体的e，w边界通过主控制体的计算节点，v控制体（蓝色双点划线网格）的计算节点在主控制体的n边，该控制体的n，s面经过主控制体的计算节点。在u的控制体中，采用有限体积法离散可得u的代数方程:£aunb£aunbnbL*+Q+(pp-pe)Ae+Q+(pp—pN)A(35)压力场的求解采用压力校正方法。即采用预估的压力场求速度，再用连续性方程校正压力场。当连续性方程得到满足时，压力场就是真实的压力场。具体步骤如下1.2.预测压力场p将预测压力场代入动量方程，分别求解速度场矶vaueeavnn=£=£au+Q+(p—p)AanbVnb+Q+(Pp-pN)A(36)3.用连续性方程校正压力设方程的精确解为u,v,pu=u+u'；v=v+vf；p—p+p'(37)其中u',v',p'为校正量。则校正量满足方程:au'avf=£a/+Q+(pf—p‘)AEnbnbPEeavf+Q+(p'—pf)A可得：nbnbaM-「略去相邻节点速度校正量的影响，(38)生(p，aPe-p')=d(p，-p?(39a)生(p，

aPn代入连续性方程:-p，)=d-p'N)(39b)(u+uf)A—(u+u')A+(v+v')A—(v+v')A=0(40)整理得：de,pKPUp=—(uAP—dA+(d，P，=—u了ELdeWAeW—Ww,P