非结构网格上多介质可压缩流模拟的RKDG方法

结构网格上多介质可压缩流模拟的RKDG方法杨广辉;欧阳洁;刘帅强;杨斌鑫【摘要】针对一般方法模拟具有运动界面的多介质可压缩流动问题计算量大、实施复杂的缺点，本文发展了一种基于非结构网格的数值模拟方法。该方法采用RKDG(Runge-KuttaDiscontinuousGalerkin)方法的弱形式求解Euler方程，用强形式求解可压缩流场模拟中的LevelSet方程，并用SimpleFix方法耦合两套方程的数值求解。二维多介质可压缩流的模拟表明：该方法成功地抑制了界面附近的非物理振荡，计算量小、实施简单，并可有效求解具有运动界面的多介质可压缩流动问题。％Sincethecostofgeneralmethodisexpensiveandtheimplementsiscomplicatedtosimulatethemultimediacompressibleflowwithmovinginterface,weproposeaRKDG(Runge-KuttaDiscontinuousGalerkin)finiteelementmethodbasedontheunstructuredgridsinthispaper.TheweakformofRKDGmethodisusedtosolvetheEulerequation,whilethestrongformofthatisappliedtosolvethelevelsetequationinthesimulationofcompressibleflows.ThenumericalsolutionsoftwoequationsarecoupledbytheSimpleFixmethod.Thesimulationoftwodimensionalcompressiblemultimediaflowshowsthatthenumericalmethodcaneffectivelyrestrainthenon-physicaloscillation.Moreover,themethodhastwoadvantages:oneislowcomputingcostandtheotheriseasytoimplement,whichmakeitiseffectivetosimulatethecompressiblemultimediafluidswithmovinginterfaces.期刊名称】《工程数学学报》年(卷),期】2014(000)005【总页数】9页(P719-727)【关键词】间断有限元；LevelSet;非结构;界面;可压缩流【作者】杨广辉;欧阳洁;刘帅强;杨斌鑫【作者单位】西北工业大学应用数学系，西安710129;西北工业大学应用数学系，西安710129;西北工业大学应用数学系，西安710129;太原科技大学应用科学学院，太原030024【正文语种】中文1引言运动界面追踪问题的研究涉及到科学工程计算、计算物理、计算化学和生物工程等学科．鉴于实际问题的求解区域往往并不规则，因此，发展适应复杂区域的高精度算法以精确捕捉运动界面非常必要．非结构网格是解决复杂区域问题最有效的技术手段•模拟运动界面的方法按参考系选取不同可分为Lagrange方法、Euler方法和Lagrange-Euler耦合方法•基于Euler网格的LevelSet方法因能提供较多的运动界面信息最为流行.LevelSet方法能精确地描述界面，较易处理复杂的物质界面；而间断有限元(DiscontinuousGalerkin,简称DG)方法作为一种新型数值方法，精度高，适合复杂求解区域，因此，耦合LevelSet方程的DG方法研究是一个十分新颖的研究方向•文献［1-3］在结构网格上采用DG方法和LevelSet方法求解了多介质可压缩流的界面运动问题，但结构网格难以在复杂的求解区域上实施•文献［4-8］采用DG方法和LevelSet方法在非结构网格上求解了不可压缩流的界面运动问题.本文选择无粘可压缩Euler方程作为物理问题的控制方程,LevelSet方程作为运动界面追踪的控制方程，并采用Runge-KuttaDG(简称RKDG)方法求解上述两套方程、由SimpleFix[1]方法实现两套方程数值求解的耦合，进而在非结构网格上模拟具有运动界面的多介质可压缩流问题，以验证该方法的有效性．控制方程Euler方程二维非稳态、无粘、可压缩气体方程-Euler方程包括质量、动量和能量三个守恒方程，其向量形式为其中P是密度，u和v分别是x和y方向的速度，p是压力，E是能量，且Y代表比热常数.LevelSet方程LevelSet方法最初由Osher和Sethian[9]提出，用于求解运动界面控制方程.其中u表示流场中的速度，申(x,t)表示点x到界面的符号距离函数，根据甲(x,t)零等值面的变化即可描述前沿界面的演变•求解LevelSet方程时，普通的数值方法由于内在耗散效应，即使只进行了几个时间步甚至一个时间步的求解，申(x,t)将不再满足是点x到界面的符号距离这一性质•为了使申(x,t)能够继续保持符号距离性质，通常需要采用“重新初始化”的方法•本文采用的DG方法求解LevelSet方程，可以避免重新初始化或者较长时间计算不用重新初始化.数值方法3.1RKDG方法DG方法最初由Reed和Hill［10］在1973年提出•他们采用完全间断的分片多项式作为近似解和试验函数空间，求解中子输运方程(这个方程是与时间无关的双曲型偏微分方程)•后来，Cockburn和Shu基于时间离散的Runge-Kutta方法与空间离散的DG方法发展了RKDG［11］，以求解含时间项的双曲守恒律方程.守恒律方程的一般形式为其中U代表守恒的物理量，可以是标量，也可以是向量•采用RKDG求解方程(4),需要先对求解区域0构造剖分Th，剖分单元用Qk表示.其次定义有限元空间Vh二{vhwL2,vh|kwVh(Qk),0kwTh},Vh(0k)是单元0k上的局部有限元空间，(j=12…,np)eVh(Qk)为单元Qk上的基函数，其中np为基函数的个数.不妨设u(x,t)只含一个分量，此时它为标量，记为u(x,t)•显然，在任意单元Qk上，式(4)同样成立，即式(5)左右两边乘以试验函数lki(x)，并在单元Qk上积分后，进行一次分部积分，可得DG方法的弱形式为Euler方程的RKDG方法含物理间断的Euler方程，如果直接采用RKDG方法离散格式求解，会在间断处出现非物理振荡，但引入限制器［11］可以抑制解的振荡.针对非结构网格节点DG的特点，本文采用文献［12］针对Euler方程提出的限制器(ShuangzhangTu和ShahrouzAllabadi提出，简称为T&A限制器)•其主要思想是限制单元上解的梯度，利用限制后的梯度，重构出单元节点上的值.本文选择HLL［13］(Harten-Lax-vanLeer)型数值通量计算Euler方程•它具有低耗散性和复杂间断易处理的特点，其表达式如下其中SL,SR分别为单元界面处的最小和最大速度，即q和c分别代表垂直于单元界面的速度和声速，记标“人”的量表示Roe平均值［12］．LevelSet方程的RKDG方法不可压缩流场中，流体的速度场散度为零，则LevelSet方程可化为类似于式(4)的守恒形式，所以可以方便地采用3.1中非结构网格上RKDG方法求解•但在可压缩流场中，如果采用DG方法弱形式求解LevelSet方程会涉及到速度场散度的计算，而速度场散度难以直接计算.本文采用文献［14］提出的DG方法强形式解决此问题，即式(6)左右两边同时加上再分部积分一次，得到强形式强形式和弱形式在数学上等价，但强形式的左端可以保留LevelSet方程原来的形式，因此，可推导出可压缩流场中LevelSet方程的DG方法强形式为式(8)左端只需要用到流场中节点上的速度值，不需要计算速度场的散度.LevelSet方程是对流方程，采用DG方法计算时，选择迎风型数值通量可以保证离散格式的稳定性•其中申LgR分别表示图1中申在单元Qk的边界上来自单元内部的插值和来自单元外部的插值.图1:二维情形迎风型数值通量4SimpleFix方法多介质流场模拟方法通过LevelSet方法得到运动界面的位置等几何特征后，在流场控制方程求解中怎样耦合这些信息十分关键•如果流场中只含一种流体，DG方法求解流体控制方程本身就能够捕捉到流场中的间断，再加上LevelSet方法辅助，将会得到更尖锐的间断界面.但如果流场中含有多种流体，则必须利用LevelSet方法确定每种流体的求解区域，并利用界面的几何信息结合相应的边界处理方法，才能精确地捕捉到物质的界面•针对流场中含多种流体的情形，由于界面两侧流体不同，气体状态方程也不同•这时，界面附近通常有两种处理方法：一种是在界面附近处理计算区域，如GhostFluid[2,3,15]方法；一种是在界面附近处理物理方程，如SimpleFix[1]方法•使用GhostFluid方法，需要计算虚拟和真实两种流场，并且每个时间步上还需使用Isobaric技术在界面处延拓间断的物理量，总的计算量非常大•一般来说，GhostFluid方法的计算量是SimpleFix方法的十几倍•为减少计算量，本文采用SimpleFix方法实现LevelSet方程和Euler方程数值求解的耦合.SimpleFix方法实施过程SimpleFix方法的思想是在界面处利用Heaviside函数光滑化比热常数Y(在界面穿过的网格单元中，“好像”存在一个介于两种流体间的虚拟流体)，以抑制物理量求解时界面处出现的振荡•因为SimpleFix方法是对物理量的控制方程做了改动，所以流场的整个求解过程都要改变，包括数值流通量和限制器等•图1中以数值流通量为例，设界面穿过单元Qk,在单元Qk上计算边界上的数值通量时，需要邻接单元上物理量的信息，此时比热常数Y要用单元Qk中的y，即单元Qk的三条边每时每刻只能“看到”一种流体•限制器也用同样的方法处理•不过这种改动只出现在界面附近的单元，其他单元不受影响•在计算中，需要通过检验各单元上距离函数值申(x,t)的大小来判定哪些单元是属于界面附近的单元.SimpleFix算法的流程如下：步骤1初始化守恒变量Qh和距离函数©0,密度在界面处光滑化处理；步骤2令i=0,1N（t0—tN）;步骤2.1利用RKDG更新守恒量Qh为Q*;步骤2.2根据Q*和第i时刻的比热常数Yi（用单元中心的比热常数），计算出原始变量（密度P*，速度场u*，压力p*）;步骤2.3在速度场u*下，计算LevelSet方程，更新符号距离函数申i—®+1；步骤2.4根据更新的©i+1以及相邻单元的比热常数Yl，Y2计算新时刻比热常数Yi+1,其中Yl，Y2为气态方程中比热常数，计算比热常数的具体公式如下步骤2.5利用Yi+1和原始变量更新守恒量步骤3输出数据．由SimpleFix的流程知，所谓的“Fix"是对比热常数Y而言的，即计算时需要在上下两个时间步上固定Y和在每个单元上固定Y-5有效性验证本节以二维情形为例，模拟Mach数为1.22的激波通过气泡的过程，以验证算法的有效性．当激波作用于气泡时，气泡将发生复杂形变，需要用高精度的计算格式以精确地追踪运动界面．同时由于界面附近物理量（密度、比热常数等）不连续，因而计算过程中界面附近易产生非物理振荡，这是该问题模拟的另一个难点．图2给出初始时刻激波打气泡问题的示意图，其无量纲化的求解区域为[0,325]x[-44.5,44.5],气泡中心位于（175,0），半径r=25，激波位于x=225处，并向左传播．求解区域采用三角网格剖分，计算单元数为13523．求解时，区域左侧和右侧分别采用出流和入流边界条件，时间方向采用二阶Runge-Kutta格式离散．5.1氦气气泡对于氦气泡情况，气泡内密度小于气泡外理想气体的密度，比热常数大于理想气体的比热常数•初始条件为-s图3分别给出了氦气气泡受到激波作用后，在T=30,50,90,180时界面的位置和形状•从图3可以看出，在激波与球形界面作用的初期，氦气泡的迎风面变形较大，呈椭圆形状，而背风面仍然保持球面形状•随着界面的进一步演化，迎风面开始向内凹陷，背风面形状依然基本保持球面形状，最终迎风面随着凹陷的加深，开始出现“尖钉”结构，并呈现类似蘑菇头状的形态•图3中界面计算结果和文献［16］中的实验结果比较吻合.图2:激波打气泡问题示意图图3:氦气气泡在不同时刻时界面位置和形状5.2R22气泡对于R22气泡情况，气泡内密度大于气泡外理想气体的密度，比热常数小于理想气体的比热常数•初始条件为图4分别给出了T=80,150,200,400时，R22气泡界面的位置和形状•由于声速在R22气泡中传播较小，而在氦气泡中的传播速度较大，因此与氦气泡的演化形态相反，激波与R22气泡相互作用后，气泡从左端向里凹陷•图4中计算结果与文献［16冲的实验结果比较吻合.图4:R22气泡在不同时刻时界面位置和形状上述激波气泡问题的数值模拟表明：基于RKDG方法与SimpleFix方法，可以较高分辨率准确捕捉到运动界面，实现可压缩介质流动问题的数值模拟，并有效抑制界面附近的非物理振荡.6结论本文建立了基于DG方法求解多介质可压缩流的流动问题的统一框架•用SimpleFix方法实现了LevelSet方程和Euler方程的耦合求解，在非结构网格上模拟了含有运动界面的多介质可压缩流流动问题，所得计算结果与实验结果吻合较好•研究结果表明：统一采用DG方法求解流场控制方程及界面运动方程是模拟多介质流动的一种有效方法.该方法使得两套方程中的数据信息传递简单而且快捷；与GhostFluid方法相比，采用SimpleFix方法实现流场控制方程与界面运动方程的耦合求解，计算量小，实施简单；该方法适合复杂求解区域，并且不依赖于问题的维数，所以可方便地扩展到高维情形；由于DG方法具有精度高、质量守恒的优点，本文采用DG方法求解LevelSet方程，可以避免重新初始化或者较长时间计算不采用重新初始化.参考文献：JorickN.ARunge-Kuttadiscontinuous-Galerkinlevelsetmethodforunsteadycompressibletwo-fluidflow[D].Delft:DelftUniversityofTechnology,2006QiuJX,LiuTG,KhooBC.Runge-KuttadiscontinuousGalerkinmethodsforcompressibletwo-mediumlfowsimulations:one-dimensionalcase[J].JournalofComputationalPhysics,2007,222(1):353-373QiuJX,LiuTG,KhooBC,etal.Simulationsofcompressibletwo-mediumflowbyRunge-KuttadiscontinuousGalerkinmethodswiththeghostfluidmethod[J].CommunicationsinComputationalPhysics,2008,3(2):479-504GroossJ,HesthavenJS.AlevelsetdiscontinuousGalerkinmethodforfreesurfaceflows[J].ComputerMethodsinAppliedMechanicsandEngineering,2006,195(25-28):3406-3429MarchandiseE,RemacleJF.Astabilizedfiniteelementmethodusingadiscontinuouslevelsetapproachforsolvingtwophaseincompressibleflows[J].JournalofComputationalPhysics,2006,219(2):780-800MarchandiseE,GeuzaineP.Astabilizedfiniteelementmethodusingadiscontinuouslevelsetapproachforthecomputationofbubbledynamics[J].JournalofComputationalPhysics,2007,225(1):949-974NguyenVT,PeraireJ.AdiscontinuousGalerkinfronttrackingmethodfortwo-phaseflowswithsurfacetension[J].Computers&Fluids,2009,39(1):1-14TurekS,MierkaO,HysingS,etal.Ahighorder3DFEM-levelsetapproachformultiphaseflowswithapplicationtomonodispersedropletgeneration[J].InternationalJournalforNumericalMethodinFluids,2010,00(1):1-10OsherS,SethianJ.Frontspropagatingwithcurvaturedependentspeed:algorithmsbasedonHamilton-Jacobiformulations[J].JournalofComputationalPhysics,1988,79(1):12-49ReedWH,HillTR.Triangularmeshmethodsfortheneutrontransportequation[R].TechnicalReportLA-UR-73-479,LosAlamosScientificLaboratory,1973ShuCW.Runge-KuttadiscontinuousGalerkinmethodsforconvection-dominate