非結(jié)構(gòu)網(wǎng)格上多介質(zhì)可壓縮流模擬的RKDG方法

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