結構網格上多介質可壓縮流模擬的RKDG方法楊廣輝;歐陽潔;劉帥強;楊斌鑫【摘要】針對一般方法模擬具有運動界面的多介質可壓縮流動問題計算量大、實施復雜的缺點，本文發(fā)展了一種基于非結構網格的數(shù)值模擬方法。該方法采用RKDG(Runge-KuttaDiscontinuousGalerkin)方法的弱形式求解Euler方程，用強形式求解可壓縮流場模擬中的LevelSet方程，并用SimpleFix方法耦合兩套方程的數(shù)值求解。二維多介質可壓縮流的模擬表明：該方法成功地抑制了界面附近的非物理振蕩，計算量小、實施簡單，并可有效求解具有運動界面的多介質可壓縮流動問題。％Sincethecostofgeneralmethodisexpensiveandtheimplementsiscomplicatedtosimulatethemultimediacompressibleflowwithmovinginterface,weproposeaRKDG(Runge-KuttaDiscontinuousGalerkin)finiteelementmethodbasedontheunstructuredgridsinthispaper.TheweakformofRKDGmethodisusedtosolvetheEulerequation,whilethestrongformofthatisappliedtosolvethelevelsetequationinthesimulationofcompressibleflows.ThenumericalsolutionsoftwoequationsarecoupledbytheSimpleFixmethod.Thesimulationoftwodimensionalcompressiblemultimediaflowshowsthatthenumericalmethodcaneffectivelyrestrainthenon-physicaloscillation.Moreover,themethodhastwoadvantages:oneislowcomputingcostandtheotheriseasytoimplement,whichmakeitiseffectivetosimulatethecompressiblemultimediafluidswithmovinginterfaces.期刊名稱】《工程數(shù)學學報》年(卷),期】2014(000)005【總頁數(shù)】9頁(P719-727)【關鍵詞】間斷有限元；LevelSet;非結構;界面;可壓縮流【作者】楊廣輝;歐陽潔;劉帥強;楊斌鑫【作者單位】西北工業(yè)大學應用數(shù)學系，西安710129;西北工業(yè)大學應用數(shù)學系，西安710129;西北工業(yè)大學應用數(shù)學系，西安710129;太原科技大學應用科學學院，太原030024【正文語種】中文1引言運動界面追蹤問題的研究涉及到科學工程計算、計算物理、計算化學和生物工程等學科．鑒于實際問題的求解區(qū)域往往并不規(guī)則，因此，發(fā)展適應復雜區(qū)域的高精度算法以精確捕捉運動界面非常必要．非結構網格是解決復雜區(qū)域問題最有效的技術手段?模擬運動界面的方法按參考系選取不同可分為Lagrange方法、Euler方法和Lagrange-Euler耦合方法?基于Euler網格的LevelSet方法因能提供較多的運動界面信息最為流行.LevelSet方法能精確地描述界面，較易處理復雜的物質界面；而間斷有限元(DiscontinuousGalerkin,簡稱DG)方法作為一種新型數(shù)值方法，精度高，適合復雜求解區(qū)域，因此，耦合LevelSet方程的DG方法研究是一個十分新穎的研究方向?文獻［1-3］在結構網格上采用DG方法和LevelSet方法求解了多介質可壓縮流的界面運動問題，但結構網格難以在復雜的求解區(qū)域上實施?文獻［4-8］采用DG方法和LevelSet方法在非結構網格上求解了不可壓縮流的界面運動問題.本文選擇無粘可壓縮Euler方程作為物理問題的控制方程,LevelSet方程作為運動界面追蹤的控制方程，并采用Runge-KuttaDG(簡稱RKDG)方法求解上述兩套方程、由SimpleFix[1]方法實現(xiàn)兩套方程數(shù)值求解的耦合，進而在非結構網格上模擬具有運動界面的多介質可壓縮流問題，以驗證該方法的有效性．控制方程Euler方程二維非穩(wěn)態(tài)、無粘、可壓縮氣體方程-Euler方程包括質量、動量和能量三個守恒方程，其向量形式為其中P是密度，u和v分別是x和y方向的速度，p是壓力，E是能量，且Y代表比熱常數(shù).LevelSet方程LevelSet方法最初由Osher和Sethian[9]提出，用于求解運動界面控制方程.其中u表示流場中的速度，申(x,t)表示點x到界面的符號距離函數(shù)，根據(jù)甲(x,t)零等值面的變化即可描述前沿界面的演變?求解LevelSet方程時，普通的數(shù)值方法由于內在耗散效應，即使只進行了幾個時間步甚至一個時間步的求解，申(x,t)將不再滿足是點x到界面的符號距離這一性質?為了使申(x,t)能夠繼續(xù)保持符號距離性質，通常需要采用“重新初始化”的方法?本文采用的DG方法求解LevelSet方程，可以避免重新初始化或者較長時間計算不用重新初始化.數(shù)值方法3.1RKDG方法DG方法最初由Reed和Hill［10］在1973年提出?他們采用完全間斷的分片多項式作為近似解和試驗函數(shù)空間，求解中子輸運方程(這個方程是與時間無關的雙曲型偏微分方程)?后來，Cockburn和Shu基于時間離散的Runge-Kutta方法與空間離散的DG方法發(fā)展了RKDG［11］，以求解含時間項的雙曲守恒律方程.守恒律方程的一般形式為其中U代表守恒的物理量，可以是標量，也可以是向量?采用RKDG求解方程(4),需要先對求解區(qū)域0構造剖分Th，剖分單元用Qk表示.其次定義有限元空間Vh二{vhwL2,vh|kwVh(Qk),0kwTh},Vh(0k)是單元0k上的局部有限元空間，(j=12…,np)eVh(Qk)為單元Qk上的基函數(shù)，其中np為基函數(shù)的個數(shù).不妨設u(x,t)只含一個分量，此時它為標量，記為u(x,t)?顯然，在任意單元Qk上，式(4)同樣成立，即式(5)左右兩邊乘以試驗函數(shù)lki(x)，并在單元Qk上積分后，進行一次分部積分，可得DG方法的弱形式為Euler方程的RKDG方法含物理間斷的Euler方程，如果直接采用RKDG方法離散格式求解，會在間斷處出現(xiàn)非物理振蕩，但引入限制器［11］可以抑制解的振蕩.針對非結構網格節(jié)點DG的特點，本文采用文獻［12］針對Euler方程提出的限制器(ShuangzhangTu和ShahrouzAllabadi提出，簡稱為T&A限制器)?其主要思想是限制單元上解的梯度，利用限制后的梯度，重構出單元節(jié)點上的值.本文選擇HLL［13］(Harten-Lax-vanLeer)型數(shù)值通量計算Euler方程?它具有低耗散性和復雜間斷易處理的特點，其表達式如下其中SL,SR分別為單元界面處的最小和最大速度，即q和c分別代表垂直于單元界面的速度和聲速，記標“人”的量表示Roe平均值［12］．LevelSet方程的RKDG方法不可壓縮流場中，流體的速度場散度為零，則LevelSet方程可化為類似于式(4)的守恒形式，所以可以方便地采用3.1中非結構網格上RKDG方法求解?但在可壓縮流場中，如果采用DG方法弱形式求解LevelSet方程會涉及到速度場散度的計算，而速度場散度難以直接計算.本文采用文獻［14］提出的DG方法強形式解決此問題，即式(6)左右兩邊同時加上再分部積分一次，得到強形式強形式和弱形式在數(shù)學上等價，但強形式的左端可以保留LevelSet方程原來的形式，因此，可推導出可壓縮流場中LevelSet方程的DG方法強形式為式(8)左端只需要用到流場中節(jié)點上的速度值，不需要計算速度場的散度.LevelSet方程是對流方程，采用DG方法計算時，選擇迎風型數(shù)值通量可以保證離散格式的穩(wěn)定性?其中申LgR分別表示圖1中申在單元Qk的邊界上來自單元內部的插值和來自單元外部的插值.圖1:二維情形迎風型數(shù)值通量4SimpleFix方法多介質流場模擬方法通過LevelSet方法得到運動界面的位置等幾何特征后，在流場控制方程求解中怎樣耦合這些信息十分關鍵?如果流場中只含一種流體，DG方法求解流體控制方程本身就能夠捕捉到流場中的間斷，再加上LevelSet方法輔助，將會得到更尖銳的間斷界面.但如果流場中含有多種流體，則必須利用LevelSet方法確定每種流體的求解區(qū)域，并利用界面的幾何信息結合相應的邊界處理方法，才能精確地捕捉到物質的界面?針對流場中含多種流體的情形，由于界面兩側流體不同，氣體狀態(tài)方程也不同?這時，界面附近通常有兩種處理方法：一種是在界面附近處理計算區(qū)域，如GhostFluid[2,3,15]方法；一種是在界面附近處理物理方程，如SimpleFix[1]方法?使用GhostFluid方法，需要計算虛擬和真實兩種流場，并且每個時間步上還需使用Isobaric技術在界面處延拓間斷的物理量，總的計算量非常大?一般來說，GhostFluid方法的計算量是SimpleFix方法的十幾倍?為減少計算量，本文采用SimpleFix方法實現(xiàn)LevelSet方程和Euler方程數(shù)值求解的耦合.SimpleFix方法實施過程SimpleFix方法的思想是在界面處利用Heaviside函數(shù)光滑化比熱常數(shù)Y(在界面穿過的網格單元中，“好像”存在一個介于兩種流體間的虛擬流體)，以抑制物理量求解時界面處出現(xiàn)的振蕩?因為SimpleFix方法是對物理量的控制方程做了改動，所以流場的整個求解過程都要改變，包括數(shù)值流通量和限制器等?圖1中以數(shù)值流通量為例，設界面穿過單元Qk,在單元Qk上計算邊界上的數(shù)值通量時，需要鄰接單元上物理量的信息，此時比熱常數(shù)Y要用單元Qk中的y，即單元Qk的三條邊每時每刻只能“看到”一種流體?限制器也用同樣的方法處理?不過這種改動只出現(xiàn)在界面附近的單元，其他單元不受影響?在計算中，需要通過檢驗各單元上距離函數(shù)值申(x,t)的大小來判定哪些單元是屬于界面附近的單元.SimpleFix算法的流程如下：步驟1初始化守恒變量Qh和距離函數(shù)?0,密度在界面處光滑化處理；步驟2令i=0,1N（t0—tN）;步驟2.1利用RKDG更新守恒量Qh為Q*;步驟2.2根據(jù)Q*和第i時刻的比熱常數(shù)Yi（用單元中心的比熱常數(shù)），計算出原始變量（密度P*，速度場u*，壓力p*）;步驟2.3在速度場u*下，計算LevelSet方程，更新符號距離函數(shù)申i—?+1；步驟2.4根據(jù)更新的?i+1以及相鄰單元的比熱常數(shù)Yl，Y2計算新時刻比熱常數(shù)Yi+1,其中Yl，Y2為氣態(tài)方程中比熱常數(shù)，計算比熱常數(shù)的具體公式如下步驟2.5利用Yi+1和原始變量更新守恒量步驟3輸出數(shù)據(jù)．由SimpleFix的流程知，所謂的“Fix"是對比熱常數(shù)Y而言的，即計算時需要在上下兩個時間步上固定Y和在每個單元上固定Y-5有效性驗證本節(jié)以二維情形為例，模擬Mach數(shù)為1.22的激波通過氣泡的過程，以驗證算法的有效性．當激波作用于氣泡時，氣泡將發(fā)生復雜形變，需要用高精度的計算格式以精確地追蹤運動界面．同時由于界面附近物理量（密度、比熱常數(shù)等）不連續(xù)，因而計算過程中界面附近易產生非物理振蕩，這是該問題模擬的另一個難點．圖2給出初始時刻激波打氣泡問題的示意圖，其無量綱化的求解區(qū)域為[0,325]x[-44.5,44.5],氣泡中心位于（175,0），半徑r=25，激波位于x=225處，并向左傳播．求解區(qū)域采用三角網格剖分，計算單元數(shù)為13523．求解時，區(qū)域左側和右側分別采用出流和入流邊界條件，時間方向采用二階Runge-Kutta格式離散．5.1氦氣氣泡對于氦氣泡情況，氣泡內密度小于氣泡外理想氣體的密度，比熱常數(shù)大于理想氣體的比熱常數(shù)?初始條件為-s圖3分別給出了氦氣氣泡受到激波作用后，在T=30,50,90,180時界面的位置和形狀?從圖3可以看出，在激波與球形界面作用的初期，氦氣泡的迎風面變形較大，呈橢圓形狀，而背風面仍然保持球面形狀?隨著界面的進一步演化，迎風面開始向內凹陷，背風面形狀依然基本保持球面形狀，最終迎風面隨著凹陷的加深，開始出現(xiàn)“尖釘”結構，并呈現(xiàn)類似蘑菇頭狀的形態(tài)?圖3中界面計算結果和文獻［16］中的實驗結果比較吻合.圖2:激波打氣泡問題示意圖圖3:氦氣氣泡在不同時刻時界面位置和形狀5.2R22氣泡對于R22氣泡情況，氣泡內密度大于氣泡外理想氣體的密度，比熱常數(shù)小于理想氣體的比熱常數(shù)?初始條件為圖4分別給出了T=80,150,200,400時，R22氣泡界面的位置和形狀?由于聲速在R22氣泡中傳播較小，而在氦氣泡中的傳播速度較大，因此與氦氣泡的演化形態(tài)相反，激波與R22氣泡相互作用后，氣泡從左端向里凹陷?圖4中計算結果與文獻［16沖的實驗結果比較吻合.圖4:R22氣泡在不同時刻時界面位置和形狀上述激波氣泡問題的數(shù)值模擬表明：基于RKDG方法與SimpleFix方法，可以較高分辨率準確捕捉到運動界面，實現(xiàn)可壓縮介質流動問題的數(shù)值模擬，并有效抑制界面附近的非物理振蕩.6結論本文建立了基于DG方法求解多介質可壓縮流的流動問題的統(tǒng)一框架?用SimpleFix方法實現(xiàn)了LevelSet方程和Euler方程的耦合求解，在非結構網格上模擬了含有運動界面的多介質可壓縮流流動問題，所得計算結果與實驗結果吻合較好?研究結果表明：統(tǒng)一采用DG方法求解流場控制方程及界面運動方程是模擬多介質流動的一種有效方法.該方法使得兩套方程中的數(shù)據(jù)信息傳遞簡單而且快捷；與GhostFluid方法相比，采用SimpleFix方法實現(xiàn)流場控制方程與界面運動方程的耦合求解，計算量小，實施簡單；該方法適合復雜求解區(qū)域，并且不依賴于問題的維數(shù)，所以可方便地擴展到高維情形；由于DG方法具有精度高、質量守恒的優(yōu)點，本文采用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