高精度DG方法用于Rayleigh-Benard自然对流问题计算分析#.docx

510152025303540高精度DG方法用于Rayleigh-Benard自然对流问题计算分析谭勤学，任静，蒋洪德(清华大学热能工程系)摘要：Rayleigh-Banard自然对流由于其流动机理的复杂性及强源项的特性，给精确模拟此类问题带来定困难。本文使用预处理间断有限元方法求解封闭方腔内的Rayleigh-Banard自然对流，一阶精度和二阶精度计算结果表明：一阶精度的计算不能有效捕捉Rayleigh-Benard自然对流的非稳定性，需二阶精度及以上的计算。有限体积及间断有限元结果对比表明：预处理间断有限元使用理想气体模型能非常有效的模拟极低Ma数下的封闭方腔内自然对流,且预处理间断有限元方法有望更准确的模拟低Ma数的流动及传热，该方法对Rayleigh-Benard自然对流问题的计算有广泛的应用前景。关键词：间断有限元，预处理方法，Rayleigh-Benard自然对流中图分类号：V231.1PreconditionDiscontinuousGalerkinMethodForRayleigh-BenardFlowTanQinxue,RenJing,JiangHongde(TsinghuaUniversity,Beijing100084)Abstract:Asthecomplexityoftheflowmechanismandbouyancysourceterm,theRayleigh-BenardNaturalConvectionishardtosimulateaccurately.Inthispaper,thepreconditioningdiscontinuousgalerkinmethodisinducedtocalculatetheRayleigh-Benardnatrualconvection.First-orderandsecond-order,sresultsshowthatsecond-orderaccuracyorhigherisneededinordertocapturetheinstabilityoftheRayleigh-Benardnatrualconvection.Theresultsoffinitevolumemethodanddiscontinuousgalerkinmethodsshowthatthediscontinuousgalerkinmethodcaneffectivesimulatethenatrualconvectionusingthefullbouyancymodal.Andtheresultsalsoindicatethediscontinuousgalerkinmethodcangetamoreaccuracyresultthanfinitevolumemethod,andwithbroadprospectstosimulatethenatrualconvection.Keywords:DiscontinuousGalerkinMethod,PreconditionMethod,Rayleigh-Benardnatrualconvection0引言热浮升驱动流动是由流体热膨胀造成的密度差,进一步产生的浮力效应引起的，并广泛存在于自然界中，如大气的自然对流现象，Rayleigh-Benard对流现象等。在实际工程应用中，热浮升力驱动的流动也广泛存在：如具有对称结构的封闭正方形腔内竖直板的自然对流换热问题,是从空调工程热环境的控制及冷冻冷藏设备内的流动和换热等工程问题中抽象出来的理论模型；旋转机械内的流动，由于受离心力哥式力的影响，在有温差情况下旋转浮升力对其内的影响非常大。这类热浮升力驱动流动其显著特征是，随着Ra数的增大其流动情况非常复杂，如RB对流其随Ra数增会出现稳定状态，稳定有规律结果进而分叉出现混沌等现象。由于这类流动机理的复杂性及强源项的特性,给精确模拟此类问题带来一定困难。IlJ断右限元方法(DGM:DiscontinuousFiniteElementMethOd)最早是由Reed和Hill在基金项目：高等学校博1学科点专项科研基金资助课题(20090002120033)作者简介：谭勤学(1986-),男，博士研究生，燃气轮机一次空气系统及高精度算法通信联系人：任静,女，教授，燃气轮机透平传热.E-mail:renjnIaiI<1973年提出的，用于中子输运方程的求解。该方法结合了有限体积方法(FVM)和有限元方法(FEM)的基本思想，特别易于处理复杂边界及边值问题，同时具有灵活处理间断的能力，且可以通过适当选取基函数，提高单元插值多项式的次数来实现任意高阶精度，另外具有容易实现自适应、并行计算等优势。由于这些优异性质，在欧洲探讨用于未来工业设计高阶CFD程序的ADlGMA计划中，间断有限元方法被认为是最有前途的算法之一。现阶段Cockburn和Shu提出的龙格库塔DGM14',以及结合混合有限元思路提出的LDGl5,(LocalDiscontinuousGaIerkin)方法在DGM框架下广泛用于Euler方程及NS方程的求解。本文旨在探讨使用高精度DGM来模拟计算热浮升力驱动的流动，以望能通过提高数值计算的精度，更好的理解热浮升力驱动机理，获得实际应用中更准确的数据。1控制方程守恒形式的三维NS方程：lU+V-Fc(U)-VFv(UyU)=S(1.1)其中E(U),乙(U,VU)分别表示无粘通量和粘性通量。对于层流问题，不考虑重力，源项S为零；对于考虑重力的自然对流等问题，控制方程形式与三维NS方程一致，但使用有效压力P班代替p,且有重力浮升力，源项S为如下计算：P物=P+0og("办)(L2)S=0(夕-R)&(P-PO)g2(P-PJg3(P-A)g可(L3)其中：石表示参考点位置，0。表示参考密度，可以使用计算域内平均密度。2预处理方法参考Weiss&Smith的预处理方法向：将NS方程中的守恒变量U用原始变量Q代替,并引入预处理矩阵，控制方程如式2.1浮+Vj(QVO)=S(2.1)Ot其中：Q=p,匕卬,",®二1Pt厂PCPJ000Pt-vrP00Pvx=vv0P0Pvyv.00PPT竺-lP匕PU>小PTH+pC0计算过程中Ur的确定是保证计算精度及稳定性的关键因素，本文中使用式2.2计算:455055606570Ur=max(0.5 小XM 嗨)(2.2)在使用预处理方法求解式2.1时，需要对无粘通量格式进行相应的修正，本文中给出文献中预处理方法常用的几种无粘通量格式:7580859095根据原始Roe格式原理，推导预处理ROe格式：小严+小访此叫"4-'渔叫Q3)同样可根据LF(Lax-FriedrichsFIUX)原理可推导相应的预处理无粘格式：I(胡)if部)'F=-(Fk+Fl-2l)-FF.-,-AO(2.3)c2rlI飒I2rlmxQ)其中4rax(4)函数表示矩阵A的最大特征值，另外对HLLC格式(TheHarten,Lax,andvanLeerwithcontactrestoration(HLLC)SCheme)预处理修正的方法可参考文献，预处理SLAU格式，预处理AUSM格式可参考文献。而粘性通量，本文采用中心格式，不需要针对预处理进行修正。3预处理间断有限元方法将式2.1与测试函数0相乘，并对其在任意单元上分部积分，得到变分形式如下：J。华。dQ+J-F-hdS-Vra=jSd.(3.1)其中：。,表示单元区域，为。表示单元边界，使用有限元空间近似函数0，,中代替函数Q,，对于间断有限元取数值解空间和测试函数空间相同，且都取为P阶多项式函数空间PJiNQh=EQ由，iWPP(3.2)i=l,=N=R4,""(33)I=I在间断有限元离散过程中取：h=i3.1 时间项离散结合式3.2及式3.3：由于其为虚拟时间积分，解的精度仅依赖于空间项的离散，与时间项无关。为简化计算过程，单元内部的值可用单元内平均值替代而不会影响到计算精度，此时可得时间项离散过程：i=Nrf3d, = Mh t k h ttt,Mjj= JQ i jd. (3.4)3.2 单元体积分结合式3.2及式3.3,可得单元体积分过程:其中i=1,.,N, j = 1,Nq,此表示体数值点数量。(3.5)3.3 单元面积分IOO结合式3.2及式3.3,可得单元面积分过程：JqF(U)hdSD.,Dij=ijWj(3.6)n,nsp其中，=1,.,%,/=1,.,",汽,表示面数值积分点数量，4J表示第j个面数值积分点的法向通量。3.4 源项积分105结合式3.2及式3.3,可得单元面积分过程：JsJn¢SJ=njSdrefT.yT.Wj(3.6).SN吸其中i=1,.,NJ=1,N,M表示面数值积分点数量，Sj表示第j个体数值积分点的源项。3.5 数值离散过程小结110总结上述离散过程，间断有限元方法离散后方程可写为下述所示半离散形。式3.7可采用多种方法求解，本文采用3步龙格库塔方法求解该方程，在此不再赘述该方法。dUt Tm .H÷jI(3.7)1151201251304计算结果及分析图1描述了该算例的物理模型及无量纲变量。，U,W的定义：方腔内左右壁面存在温差，上下壁面绝热，计算网格为50*50的壁面加密结构化网格，使用理想气体模型考虑浮升力作用，并采用预处理方法保证在非常低Ma数下DGM方法的收敛性。图1封闭方腔内自然对流问题描述图2给出了下分别采用一阶精度、二阶精度的温度分布：一阶精度计算由于耗散比较大，计算域内流场比较稳定，不会出现RayIeigh-Benard现象。二阶精度计算能较好模拟Rayleigh-Benard对流。此计算结果表明：RayIeigh-Benard流动对计算精度敏感，计算精度不足甚至会得到一个非物理解，因此有必要对该问题进行高精度计算。(a)一阶精度(b)二阶精度图2计算精度对计算结果的影响:温度分布计算Ra=IO3,io4o5o6四种情况，计算得到无量纲温度及无量纲速度分布如图3,图4,图5所示。与DaViS结果比较可知,本文预处理间断有限元方法计算得到的分布情况与其一致。135140145150(a)Ra=103(b)Ra=104(c)Ra=105(d)Ra=10ft图3Ra=IO3,IO4,105,IO6,无量纲速度U分布同参考文献云图范围0,1,间隔0.1图4Ra=IO3,IO",105,IO6,无量纲速度V分布同参考文献【川，云图范围0,1,间隔0.1(a)Ra=10(b)Ra=104(c)Ra=105(d)Ra=106图5Ra=l(P,104,105,IO6,无量纲速度W分布同参考文献Ul云图范