劉君， 劉瑜， 陳澤棟

1.大連理工大學 航空航天學院， 大連　116024 2.裝備學院 航天裝備系， 北京　101416 3.大連理工大學 工程力學系， 大連　116024

?

非結構變形網(wǎng)格和離散幾何守恒律

劉君1,*， 劉瑜2， 陳澤棟3

1.大連理工大學 航空航天學院， 大連116024 2.裝備學院 航天裝備系， 北京101416 3.大連理工大學 工程力學系， 大連116024

數(shù)值模擬流固耦合問題或多體分離問題的非定常流動時，常采用基于任意拉格朗日-歐拉(ALE)方程的有限體積法，涉及到變形網(wǎng)格和離散幾何守恒律。在對變形網(wǎng)格算法進行綜述時，按照構造思想分為物理比擬、橢圓光順、插值、運動子網(wǎng)格(MSA)及其混合法共5類，分別介紹了基本原理、研究現(xiàn)狀和適用范圍，通過算例比較表明，徑向基函數(shù)(RBFs)和運動子網(wǎng)格相結合的混合方法既有很好的變形能力，也有較高的計算效率，值得進一步發(fā)展和推廣。在介紹了離散幾何守恒律(DGCL)概念之后，采用二維幾何模型進行分析，指出其機理是離散過程中體積增量與網(wǎng)格面元掃過的體積不相等造成的，把目前國內外應用的算法分為面積修正法、給定速度的面積修正法、速度修正法和體積修正法共4類，對其應用范圍和存在的問題進行討論，認為提出的體積修正算法既可以保證流固界面條件，也可以用于時間多層格式。

非結構網(wǎng)格； 網(wǎng)格變形； 非定常流動； 幾何守恒律； 界面耦合

有相對運動邊界的非定常流動在自然界很常見，比如飛鳥、游魚、走獸、以及人體內的血液流動等。這類流場具有相同的特點，即物體的形狀、位置隨時間變化，固體與流體相互作用，此類動力學系統(tǒng)常稱為流固耦合(Fluid-Structure)問題。在航天與航空領域還有一類是流場中有多個物體，它們之間存在相對運動，而且大部分情況下運動還受到流場的影響。例如，外掛物與載機分離、座艙蓋/座椅的彈射、多級火箭的級間分離、多彈頭再入、子母彈拋撒、爆炸彈片飛散和沖擊波驅動物體運動等，需要采用流體方程和六自由度剛體運動方程(Fluid-Rigid)耦合計算，為區(qū)別于前一類流固耦合問題，本文稱為多體分離問題。

對非定常流動進行數(shù)值模擬時，物理空間隨邊界運動而變化，導致計算中網(wǎng)格也要隨時間變化，故需要發(fā)展與流場計算過程相關的網(wǎng)格技術。根據(jù)計算中流場內網(wǎng)格數(shù)目是否變化分為兩類算法：一種是網(wǎng)格沒有變形、通過增減網(wǎng)格數(shù)目實現(xiàn)物理空間變化，目前主要有笛卡兒自適應網(wǎng)格和運動嵌套網(wǎng)格；另一種是保持計算網(wǎng)格數(shù)目不變、采用變形動網(wǎng)格技術描述流場物理空間隨時間的變化。

嵌套網(wǎng)格又稱為重疊網(wǎng)格，采用多塊貼體網(wǎng)格處理復雜外形，在物理邊界處或網(wǎng)格重疊區(qū)內建立人工邊界，稱為“洞”邊界。“洞”內單元是不參與計算的無效網(wǎng)格，“洞”外單元是參與計算的有效網(wǎng)格，“洞”邊界和相鄰子塊之間采用插值公式建立信息傳遞關系。國內采用嵌套網(wǎng)格計算定常流動的文獻較多，最早的應用可以追溯到1991年楊永健和張魯民研究帶有4個助推火箭的復雜流場模擬[1]，閻超等在自動判明邊界和快速挖洞等方面有算法的改進，在復雜外形的應用方面也取得實效[2]。近幾年國內采用子塊在靜止背景塊上運動的嵌套網(wǎng)格技術模擬包含運動邊界的流動[2-5]，由于計算過程中各子網(wǎng)格塊存在相對運動，子網(wǎng)格相對位置隨時間變化，因此每一個時間步都需要重新確定“洞”邊界，并重新判斷每一個網(wǎng)格單元是否為有效網(wǎng)格，這對于網(wǎng)格量數(shù)百萬乃至上千萬的應用問題將顯著增加計算量，在外形復雜、分塊較多情況下，描述子塊之間關系時比較繁瑣，算法的程序實現(xiàn)也較為困難。

笛卡兒網(wǎng)格(Cartesein Grids)采用四邊形(二維)或六面體(三維)形狀的網(wǎng)格，并與直角坐標系的坐標軸平行，是一種出現(xiàn)最早的用于計算流體力學(CFD)計算的網(wǎng)格劃分方法。笛卡兒網(wǎng)格具有生成簡單、流動求解器容易實現(xiàn)等優(yōu)點，但是難于描述復雜幾何外形。為了提高邊界附近流場的計算精度，通常采用局部網(wǎng)格加密技術或網(wǎng)格分區(qū)技術，在幾何形狀不規(guī)則的壁面附近進行特殊處理，這些邊界處理算法會導致數(shù)據(jù)結構變復雜，降低計算效率，因此笛卡兒網(wǎng)格一直沒有成為主流。采用自適應笛卡兒網(wǎng)格可以實現(xiàn)運動邊界，但在計算過程中需要建立網(wǎng)格之間的動態(tài)數(shù)據(jù)鏈表，且每一步計算都需要重新判斷邊界單元，在網(wǎng)格量較大的情況下，大量繁雜的搜索和插值過程將消耗巨大的計算機資源，甚至會遠大于求解流體力學方程本身。目前國內對笛卡兒網(wǎng)格的研究和應用非常少，其中李盾等采用該技術模擬多體分離問題取得很好效果[6-7]。

變形動網(wǎng)格早期主要用于激波裝配法模擬定常流場，近20年逐步回歸到邊界運動問題的計算。本文作者所在課題組從1998年開始從事變形動網(wǎng)格研究，并于2009年出版相關專著[8]，下面將對這一領域的研究歷史和最新進展進行綜述。提出以上幾種網(wǎng)格技術的年代很早，例如，變形動網(wǎng)格在20世紀70年代曾被廣泛應用，形成的相關文獻非常多，因此本文選擇文獻時對國外主要采用標志性的，對國內盡可能全面地介紹研究現(xiàn)狀。

1　非結構動網(wǎng)格技術

網(wǎng)格生成技術曾是CFD的重要內容，先后出現(xiàn)了結構網(wǎng)格、疊合網(wǎng)格、搭接網(wǎng)格、非結構網(wǎng)格、直角網(wǎng)格、混合網(wǎng)格等提法。按照數(shù)據(jù)結構是否有序，主要分為結構和非結構兩類，非結構網(wǎng)格除了相鄰單元編碼無序外，也包括控制體形狀和面元數(shù)目不受限制。嚴格的說，混合網(wǎng)格是指分塊以后有些采用有序編碼、有些采用無序編碼的算法，現(xiàn)在把非結構網(wǎng)格中出現(xiàn)六面體形狀的網(wǎng)格稱為“混合網(wǎng)格”的提法不準確。

很多文獻按照求解方程類型把網(wǎng)格變形方法分為代數(shù)方法和偏微分方程方法兩類，本文按照構造思想分為物理比擬、橢圓光順、插值、運動子網(wǎng)格及其混合這5種方法。網(wǎng)格變形有3個基本要求：變形中不出現(xiàn)失效單元，保持較好品質；能夠精確描述邊界，即保證網(wǎng)格是貼體的；不能造成CFD求解器效率大幅度降低。下面按照這3條標準來評價這5種方法。

1.1物理比擬方法

這是一類將網(wǎng)格點的運動比擬為某種物理過程的網(wǎng)格變形方法，常用的包括彈簧比擬法和固體比擬法。

1990年Batina首次采用彈簧比擬法實現(xiàn)三角形網(wǎng)格的變形[9]。它的原理是將網(wǎng)格單元的邊視為彈簧，整個計算網(wǎng)格構成一個彈簧系統(tǒng)，在彈簧系統(tǒng)內部節(jié)點(網(wǎng)格點)建立力的平衡方程組，邊界的運動或變形導致彈簧系統(tǒng)受到拉伸或壓縮，彈簧節(jié)點受力發(fā)生變化，通過求解整個彈簧系統(tǒng)節(jié)點受力的平衡方程組來確定網(wǎng)格節(jié)點新的位置。這種方法只能描述彈簧的拉壓變形，扭轉變形易導致網(wǎng)格折穿(Snap-Through)現(xiàn)象。1998年Farhat等在二維網(wǎng)格的彈簧節(jié)點處通過施加扭轉彈簧增加了網(wǎng)格單元的抗扭轉能力，并取得良好效果，但在推廣到三維網(wǎng)格時遇到了困難[10-11]。2000年Blom根據(jù)彈簧邊所對應的三角形單元的頂角對線彈簧剛度系數(shù)進行修正，提出了半扭轉彈簧方法[12]。2003年郭正等針對四面體單元，采用不包含彈簧邊的二面角作為頂角，把該方法推廣到了三維情形，同時通過求解固體導熱方程得到內部節(jié)點的溫度，并將其作為參數(shù)用來增加運動邊界附近網(wǎng)格層的彈簧剛度系數(shù)，提高了網(wǎng)格變形能力[13]。彈簧方法只適用于三角形網(wǎng)格和四面體網(wǎng)格，應用于其他類型網(wǎng)格則需要剖分重構，例如，六面體單元剖分6個四面體網(wǎng)格[8]，增加了計算量和實施難度，限制了該方法的應用。2005年Bottasso等提出球點彈簧方法，在以網(wǎng)格邊構成的彈簧系統(tǒng)內再添加新彈簧，這些新彈簧的起點位于網(wǎng)格節(jié)點，終點止于該點在所對面內的垂足，如此一來，格點便為多面球所包圍，從而避免折穿，可以用于四邊形和六面體單元的網(wǎng)格變形[14-15]。彈簧比擬法所求解的力平衡方程組對角占優(yōu)，采用Jacobi迭代或Seidel迭代較快收斂，因此得到大量應用，但其變形能力比固體比擬法差。

固體比擬法是Tezduyar等于1992年提出的，它將計算域視為彈性體，將運動邊界的位移作為施加在彈性體上的變形載荷，在每個網(wǎng)格節(jié)點上建立一組彈性力學基本方程，通過計算該方程得到在上述變形載荷下各個網(wǎng)格節(jié)點的位移，從而確定新的網(wǎng)格節(jié)點位置[16]。彈性體方法的模型完備，具有比較成熟的計算方法，適用于任何網(wǎng)格類型，通過調整彈性模量、泊松比、添加源項等改進措施，使其具有很強的網(wǎng)格變形能力，至今依然在應用[17-26]。

還有其他的物理比擬方法，如Kennon等提出的流動比擬法[27]和陳炎等提出的溫度體比擬法[28-31]，因其缺少普適性，這里不再介紹。

物理比擬方法中，彈性體法需要求解偏微分方程組，計算量大，不適用于大規(guī)模問題，彈簧法使用簡單，目前應用最為廣泛，但是變形能力較差。

1.2網(wǎng)格光順方法

基于偏微分方程的網(wǎng)格生成算法很早就被廣泛應用[32]，1996年L?hner等用于網(wǎng)格變形，通過求解Laplace方程把邊界運動逐步傳遞到計算區(qū)域內部網(wǎng)格點，由于橢圓型方程具有在求解域內部不會產生極值的數(shù)學特性，這種方法得到的網(wǎng)格變形光滑過渡，因此稱為網(wǎng)格光順方法[33]。計算均勻擴散系數(shù)的Laplace方程得到的位移主要和網(wǎng)格點所處位置相關，同樣變形相對于不同網(wǎng)格單元的效果不一樣，通過調整擴散系數(shù)讓小單元經歷較小變形、大單元經歷較大變形，提高網(wǎng)格變形能力[34-38]。近幾年探索從雙調諧方程和雙橢圓方程對網(wǎng)格實施變形[39-40]。由于彈性力學基本方程和Laplace方程的數(shù)學性質類似，實際應用也表明網(wǎng)格光順法和彈性體法相似，具有較強的網(wǎng)格變形能力，但是計算量大，實施較復雜。

1.3插值方法

采用插值函數(shù)根據(jù)已知邊界網(wǎng)格點位移插值得到內部網(wǎng)格點位移，分為顯式和隱式兩種。沿著網(wǎng)格線根據(jù)確定的邊界位置，采用代數(shù)插值、多項式插值、樣條插值等直接計算網(wǎng)格內點位置的顯式插值方法，在激波裝配法、氣動彈性計算中使用幾十年了，由于這種簡單插值僅適用于結構網(wǎng)格和簡單外形，因此很少作為變形網(wǎng)格技術討論。近幾年通過采用基于距離的插值函數(shù)建立了適用于復雜外形、任意形狀網(wǎng)格的顯式插值[41-44]，其中Shepard提出的逆距加權插值(Inverse Distance Weighting Interpolation， IDW)[45],采用關于未知點到已知點距離倒數(shù)的函數(shù)作為權重，通過已知點集的數(shù)據(jù)的加權平均得到未知點處的數(shù)據(jù)，取得很好效果[46-48]。顯式插值具有很高的計算效率，適合于大規(guī)模的并行計算，但是不考慮連接關系的“點對點”特性使其變形后網(wǎng)格質量很快下降，在外形或邊界運動等情況下，其效果遠不如考慮相鄰點影響的隱式插值。常用的隱式插值是徑向基函數(shù)(Radial Basis Functions， RBFs)法，采用點與點之間徑向距離的函數(shù)，在已知運動邊界離散點數(shù)據(jù)條件下，通過求解線性方程組得到插值函數(shù)的系數(shù)，然后得到內部網(wǎng)格點位置。該算法數(shù)據(jù)結構簡單，可以應用于任何類型網(wǎng)格，并且有完備的理論基礎[49]。2007年Boer等用于網(wǎng)格變形，發(fā)現(xiàn)計算插值系數(shù)需要的線性方程組規(guī)模是邊界點數(shù)目的平方，對于三維問題計算量很大[50]，后來對這種方法進行改進減少了計算量，但算法變得較為復雜[51-54]。

盡管插值方法很少發(fā)生變形失敗，但是主要考慮距離信息的特性使其難以保證網(wǎng)格質量，RBFs涉及到矩陣求逆，網(wǎng)格規(guī)模較大后計算效率也值得考慮。黏性流動計算時采用長寬比較大的四邊形或六面體網(wǎng)格來分辨邊界層，常規(guī)的變形算法難以滿足要求，插值方法具有明顯的優(yōu)勢。

1.4運動子網(wǎng)格方法

運動子網(wǎng)格法早期特指Delaunay圖映射方法，是Liu等在2006年提出的一種動網(wǎng)格方法[55]，在所有動邊界網(wǎng)格點和若干控制點基礎上生成Delaunay圖，也就是子網(wǎng)格，計算網(wǎng)格包含在子網(wǎng)格內部，每個計算網(wǎng)格點必然屬于子網(wǎng)格中的一個單元，可以在子網(wǎng)格單元內部得到計算網(wǎng)格點的相對位置坐標，即得到計算網(wǎng)格節(jié)點和子網(wǎng)格單元的映射關系，邊界運動牽引子網(wǎng)格運動，利用映射關系便可以得出計算網(wǎng)格節(jié)點在變形后的坐標。Delaunay圖映射方法具有很高的計算效率，但對于復雜外形和大幅度運動，網(wǎng)格質量下降很快，重新生成需要花費時間。近年來，部分學者提出一些改進算法，肖天航等提出了雙重Delaunay圖動網(wǎng)格生成算法[56]，在一定程度上增強動網(wǎng)格的質量和變形的魯棒性，周璇等采用Delaunay圖映射方法計算運動子網(wǎng)格點相對于動邊界的相對位置，比較方便地實現(xiàn)了對彈簧近似方法的逐層邊界修正，利用彈簧比擬法計算運動子網(wǎng)格的變形[57]，可以用于黏性計算。這種背景網(wǎng)格變形牽引內部網(wǎng)格運動的思路，也不僅僅局限于Delaunay圖，例如，Lefran?ois采用的運動子網(wǎng)格包含了計算域內的點[58]，且運動子網(wǎng)格可以根據(jù)不同的問題進行調整，比Delaunay圖要靈活，因此，采用運動子網(wǎng)格方法表述這類方法比Delaunay圖映射方法更合適。

運動子網(wǎng)格方法的優(yōu)點是計算效率高，缺點是用于復雜構型生成背景網(wǎng)格需要經驗，大幅度運動特別是大角度轉動易出現(xiàn)背景網(wǎng)格交叉引起的變形失敗，重新生成背景網(wǎng)格常需要人工干預。

1.5變拓撲方法

以上列舉的網(wǎng)格變形方法，網(wǎng)格點之間的連接關系不改變，在動邊界經歷大位移和大變形時，網(wǎng)格質量必然會變差，如果放棄網(wǎng)格點之間連接關系不變的限制，則動網(wǎng)格的變形能力會進一步增強。2011年Alauzet和Olivier用彈性體比擬法模擬葉片旋轉，允許改變網(wǎng)格節(jié)點間的連接關系，極大增強了變形能力，實現(xiàn)了二維葉片的360° 旋轉[59-60]，Isola等采用對質量差的網(wǎng)格單元進行邊交換也具有同樣效果[61-62]。2013年Wang和Perssony用彈簧比擬法進行二維網(wǎng)格變形，通過改變網(wǎng)格節(jié)點間的連接關系模擬了雙翼型撲動流場[63]。

連接關系確定后計算量很少，變拓撲方法的計算效率最高，但是數(shù)據(jù)結構變化較頻繁，需要求解器進行改造，實現(xiàn)較為復雜，目前應用范圍較窄。

1.6混合方法

混合方法是指同時采用兩種以上方法的組合對網(wǎng)格進行變形。最早用于多塊網(wǎng)格求解復雜外形[64]，用彈簧比擬法計算多塊結構網(wǎng)格的塊頂點位移，然后用超限插值方法更新塊內網(wǎng)格節(jié)點坐標。Zhou和Li通過多種方法的組合，得到了一種混合網(wǎng)格變形方法[65]。

本文作者注意到這種方法的優(yōu)勢在于兼顧計算效率和網(wǎng)格形狀適應性，沿著這一思路把徑向基函數(shù)方法作為計算網(wǎng)格，比較了彈簧法和運動子網(wǎng)格方法作為變形控制背景網(wǎng)格，發(fā)現(xiàn)后者具有變形能力強、計算效率高的優(yōu)點，提出RBFs-MSA網(wǎng)格變形算法[66]，下面通過具體算例進行介紹。

2　RBFs-MSA網(wǎng)格變形算法

RBFs-MSA也屬于混合方法，在生成如圖1所示的計算網(wǎng)格時，還需要生成背景網(wǎng)格，見圖2。可以看出，這種混合方法綜合了其他算法的優(yōu)勢，計算網(wǎng)格幾何形狀任意，有效滿足黏性流動模擬需要的異向加密需求，背景網(wǎng)格為三角形(二維)或四面體(三維)，具有強的變形能力。流場計算過程中，兩套網(wǎng)格按照如下程序組合使用：

1) 在計算網(wǎng)格基礎上通過網(wǎng)格細化或者重新生成的方法得到子網(wǎng)格，子網(wǎng)格邊界點和內點根據(jù)具體問題自由控制，但是數(shù)目遠少于計算網(wǎng)格。

2) 計算網(wǎng)格內點在運動子網(wǎng)格單元中的相對面積或體積坐標(映射關系)。

3) 計算動邊界的位移，更新計算網(wǎng)格邊界點坐標，這些點中包含子網(wǎng)格的邊界點。

4) 利用徑向基函數(shù)插值得到運動子網(wǎng)格內點的坐標。

5) 利用映射關系得到計算網(wǎng)格內點坐標。

圖1　計算網(wǎng)格Fig.1　Computational mesh

圖2　背景網(wǎng)格Fig.2　Background mesh

國外常用來檢驗網(wǎng)格變形算法的經典算例是Anderson等完成的撲翼實驗[67]。弦長記為c,其上下沉浮運動和俯仰運動規(guī)律分別為

(1)式中：h和h0分別為沉浮運動的位移和振幅；ω為振動頻率；t為時間；φh和φθ分別為沉浮和俯仰運動的相位角；θ和θ0分別為俯仰角和俯仰運動的振幅。

對于周期運動，更關心的是沉浮和俯仰運動的相位角φh和φθ的差：ψ=φθ-φh=75°。在實驗中采用機械結構實現(xiàn)組合運動的振動頻率ω需要根據(jù)機構運動特性計算，Anderson引入最大名義迎角α0=15°和后緣點最大位移At定義的斯特勞哈數(shù)St=ωAt/(2πu)，u為參考速度，對應的減縮頻率k=0.5ωc/u。在給定振幅h0=0.75c和St以后，結合減縮頻率公式、斯特勞哈數(shù)定義公式以及名義迎角定義公式，可計算沉浮和俯仰組合運動過程中角振幅θ0和振動頻率ω：

α0=arctan(2kh0)-θ0

(2)

計算區(qū)域以翼型前緣點為圓心，半徑為15c的半圓后面連接長度為15c的矩形。計算網(wǎng)格采用C形貼體網(wǎng)格，總共有25 200個四邊形網(wǎng)格，翼型表面網(wǎng)格上分布有270個節(jié)點，在翼型附近采用逐層加密，距物面第一層網(wǎng)格單元的高度為0.000 02c，并在前后緣進行適當加密。背景網(wǎng)格含有408個三角形單元，224個節(jié)點，其中翼型邊界節(jié)點為21個。

本文比較了背景網(wǎng)格整體剛性運動、背景網(wǎng)格采用RBFs和兩套網(wǎng)格RBFs-MSA這3種網(wǎng)格運動方法計算十步所需的CPU時間，分別為40 s、341 s和56 s。第一種方法采用直接解析表達式，因此計算網(wǎng)格更新所用時間相比于計算流動方程的時間幾乎可忽略，由此可見RBFs-MSA方法可以比RBFs大量減少計算時間。

本文分別進行層流和湍流的模擬，湍流模式選擇Spalart-Allmaras (S-A)模型和Wilcox的k-ω模型，圖3是平均推力系數(shù)CT和平均輸入功率系數(shù)CP與實驗的比較。本文計算的推力系數(shù)與Shyy[68]和Young[69]的計算結果符合較好，所有計算值略低于實驗值。在St<0.3范圍內，本文計算的平均輸入功率系數(shù)與文獻符合較好，當St>0.4時，計算結果比文獻偏大。

圖3　平均推力系數(shù)和平均輸入功率系數(shù)比較Fig.3　Comparison between average thrust coefficients and average input power coefficients

3　離散幾何守恒律

變形網(wǎng)格技術常采用基于任意拉格朗日-歐拉(Arbitrary Lagrangian-Eulerian， ALE)方程描述的有限體積法。對于三維非定?？蓧嚎sNavier-Stokes方程為

(3)

式中：Ω為控制體；?Ω為控制邊界；Q為守恒量；Fc為對流通量；Fv為黏性通量；dV為體積微元；dS為面積微元；n為控制體邊界外法向向量；xc為網(wǎng)格運動速度。由于xc是時間和空間的函數(shù)，在離散以后可能會引入新的誤差，即使對于均勻流場也難保證式(3)嚴格成立，為此，需要考慮幾何守恒律(GeometricConservationLaw，GCL)。

幾何守恒律概念最早于1978年在AIAA會議上由Thomas和Lombard提出[70-71]，當時主要討論有限差分法，也提及ALE有限體積法，給出的幾何守恒律方程為

(4)

(5)

文獻[8]中證明以上網(wǎng)格的幾何不等引起流場誤差可以達到O(1)量級，因此，不滿足DGCL對ALE方程的有限體積法計算會產生很大影響。由于DGCL與求解流場的時空離散格式相關，國外文獻構造了很多種幾何守恒律算法，選擇的修正參數(shù)不同，本文作者把國外的算法分為3類：面積修正法(包含給定速度的面積修正法)、速度修正法和體積修正法。

圖4　離散幾何守恒律示意圖Fig.4　Schematic diagram of discrete geometricconservation law

3.1面積修正法

3.2給定速度的面積修正法

(6)

為滿足DGCL要求，Guillard和Farhat提出了兩種實現(xiàn)方法[72]：一種是在tn和tn+1時刻之間構造多套網(wǎng)格，假設網(wǎng)格按照已知規(guī)律運動和變形，在這些網(wǎng)格上分別計算通量，然后對這些通量進行加權平均；另一種是在構造的多套網(wǎng)格中得到一個平均網(wǎng)格，然后在這個平均網(wǎng)格構型上計算通量。為了保證式(4)的離散格式具有時間二階精度，需要涉及20個參數(shù)來構造多個時間層的網(wǎng)格運動速度和中間面積[73-74]。這類算法推導過程復雜，增加了計算量，另外，通量計算時所用的流場變量都是tn+1時刻的流場變量，這樣處理至少在物理意義上是不明確的，國內還沒有看到具體應用。

3.3速度修正法

Mavriplis和Yang發(fā)現(xiàn)按照Guillard和Farhat方法確定的參數(shù)不惟一，針對BDF2格式，有些情況下并不都嚴格滿足DGCL條件，于是提出采用根據(jù)DGCL條件計算出來的速度ug來代替實際的網(wǎng)格運動速度xc的修正速度算法[75]，即

(7)

Mavriplis和Yang的方法不需要引入中間網(wǎng)格，計算通量時網(wǎng)格位置與流場變量在時間上是一致的，消除了Guillard和Farhat方法中物理量與網(wǎng)格在時間上的不匹配。但是ug≠xc又會引起3.4節(jié)討論的流固耦合界面算法的精度問題。

Hyung和Yannis[76]曾提出通過在控制方程右端引入源項的方式來消除DGCL誤差，雖然出發(fā)點有所不同，但最終落腳點還是在于修正網(wǎng)格運動速度。

值得指出的是張來平等也提出過類似的思想來處理一階Euler后差格式的DGCL算法[77]，先用面積修正法計算平均面積，然后采用式(8)計算網(wǎng)格速度，即

(8)

盡管不滿足DGCL會給流場計算引入誤差，文獻[78]的數(shù)值試驗表明誤差可以達到O(1)量級，但是誤差變化規(guī)律還缺乏理論研究，有些格式表現(xiàn)為周期振蕩、有些格式表現(xiàn)為單調增加，由此引出的DGCL和ALE格式穩(wěn)定性之間的關系是目前學術界主要爭論的問題。Farhat和Geuzaine[74]認為滿足DGCL是ALE格式具有和靜止網(wǎng)格格式相同穩(wěn)定性的充分必要條件，Boffi和Gastaldi[79]給出了相反的結論，Masud[80]在滿足DGCL的前提下，對網(wǎng)格速度對穩(wěn)定性和收斂性的影響進行了研究，認為網(wǎng)格運動速度能夠影響穩(wěn)定的時間步長，且計算誤差與網(wǎng)格相對運動速度有關。

3.4體積修正法

上述3種修正算法均采用修正右端項的思路來解決問題，本文提出完全不同的新思路，即對于以上BDF2格式，按照不同的時間層，有3種修正方法：

1) 修正tn+1時間層的體積Vn+1，計算時用體積V′代替進行計算，即

(9)

2) 修正tn時間層的體積Vn，計算時用V″代替，即

(10)

3) 修正tn-1時間層的體積Vn-1，計算時用V?代替，即

(11)

在非定常計算中，采用上述不同的修正方法都能在一定程度上消除DGCL誤差，但在不同的計算問題中，實際效果存在一定差異。張來平等從截斷誤差的觀點出發(fā)，曾對均勻流和等熵渦問題采用不同的修正方法進行過測試、分析[81]。

本文構建了如下兩個模型來驗證DGCL算法。計算中離散幾何守恒律分別采用Mavriplis算法[75]、張來平算法[77]和本文提出的3種算法，除了張來平算法采用一階Euler后差格式，其余采用BDF2格式，CFL數(shù)都取為200。

第1個模型是在x∈[-2,2]和y∈[0,1]構成的矩形內充滿均勻密度ρ=1、壓力p=1/γ的靜止氣體(以氣體聲速作為速度無量綱參數(shù))，計算中采用左右兩條邊上網(wǎng)格點靜止，初始時刻上下兩條邊上中點處的網(wǎng)格點在水平方向正弦運動，從而帶動整個流場的網(wǎng)格運動和變形。計算中焓的絕對誤差隨無量綱時間t的變化如圖5(a)

圖5　計算誤差Fig.5　Calculation errors

所示，其中橫坐標為無量綱時間，縱坐標為全場焓值絕對誤差的2-范數(shù)，誤差為10-15量級，接近計算機隨機誤差。

第2個模型考核網(wǎng)格形狀的影響，取x∈[0,1]和y∈[0,1]構成的正方形計算域，結點為51×51的四邊形網(wǎng)格，保持邊界點不動，內點按照如下規(guī)律運動：

(12)

式中：(x0,y0)和Δx=Δy=0.02是網(wǎng)格的初始位置和間距。從焓的絕對誤差隨無量綱時間變化(如圖5(b)所示)看出，誤差也接近計算機隨機誤差。

需要指出的是體積修正法對非均勻流場存在守恒性問題，理論上只能用于空間二階精度的有限體積法。

4　界面耦合算法

從圖5(a)來看，本文的體積修正法和速度修正法沒有本質差別，但是考慮到流固耦合問題或多體分離問題的邊界是運動的，那么速度修正法不能保證界面耦合條件的精度。下面對此進行討論。采用時間相關法求解流體力學方程本質上是時間變量的一次積分計算，界面上流體速度記為vf，CFD在離散空間上構建，這決定了從t推進到t+Δt過程中速度保持為常數(shù)，因此，邊界位置更新為

(13)

式中：rf為界面上的流體位移。

固體運動學方程本質是時間變量的二次積分，以剛體動力學方程mrtt=Fair為例來說明(m為物體質量，rtt為物體位移對時間的二階導數(shù)，即加速度，F(xiàn)air為物體氣動力的合力)，從t推進到t+Δt過程中加速度保持為常數(shù)，采用最簡單的式(14)可以看出速度是變化的。

(14)

式中：vs為界面上的固體運動速度；rs為界面上的固體運動位移；m為質量；Fair為氣動力。

如果DGCL采用面積修正算法，與網(wǎng)格運動速度無關，以上滿足位置相等條件。如果采用假設運動規(guī)律的面積修正算法，可以在二階精度條件下實現(xiàn)位置相等[74]。如果采用速度修正法計算，網(wǎng)格運動速度ug既不能保證等于流體界面速度vf，也不等于剛體運動速度vs，更無法滿足位置相等條件，因此，這類算法在流場和結構之間信息交換過程中必然引入誤差。本文提出的體積修正算法與網(wǎng)格運動速度和面積變形均無關，使得DGCL算法與流固界面信息傳遞算法不再相關，可實現(xiàn)嚴格滿足位置或速度條件的要求。

5　結　論

1) 徑向基函數(shù)和運動子網(wǎng)格相結合的混合方法既有很好的變形能力，也有較高的計算效率，值得進一步發(fā)展和推廣。

2) 需要滿足離散幾何守恒律的本質原因是離散過程中體積增量與網(wǎng)格面元掃過的體積不相等。

3) 體積修正算法既可以保證流固界面條件，也可以用于時間多層格式，具有良好的應用前景。

[1]楊永健, 張魯民. 應用重疊網(wǎng)格技術求解復雜組合體無粘流場[J]. 空氣動力學報, 1991, 31(1): 51-57.

YANG Y J, ZHAMG L M. Numerical simulation of inviscid flow over a complicated body using an overlapping grid technique[J]. Acta Aerodynamica Sinica, 1991, 31(1): 51-57 (in Chinese).

[2]范晶晶, 閻超, 張輝. 重疊網(wǎng)格洞面優(yōu)化技術的改進與應用[J]. 航空學報, 2010, 31(6): 1127-1133.

FAN J J, YAN C, ZHANG H. Improvement of hole-surface optimization technique in overset grids and Its application[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(6): 1127-1133 (in Chinese).

[3]張培紅, 王明, 鄧有奇, 等. 武器分離及艙門開啟過程數(shù)值模擬研究[J]. 空氣動力學學報, 2013, 31(3): 277-281.

ZHANG P H, WANG M, DENG Y Q, et al. Numerical simulation of store separation and door operation[J]. Acta Aerodynamica Sinica, 2013, 31(3): 277-281 (in Chinese).

[4]肖中云, 江雄, 牟斌, 等. 并行環(huán)境下外掛物動態(tài)分離過程的數(shù)值模擬[J]. 航空學報, 2010, 31(8): 1509-1516.

XIAO Z Y, JIANG X, MOU B, et al. Numerical simulation of dynamic process of store separation in parallel environment[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(8): 1509-1516 (in Chinese).

[5]劉學強, 伍貽兆, 夏健. 多重網(wǎng)格法在非結構網(wǎng)格中的應用[J]. 計算物理, 2002, 19(4): 357-361.

LIU X Q, WU Z Y, XIA J. Application of the multigrid method in unstructured meshes[J]. Chinese Journal of Computational Phyfics, 2002, 19(4): 357-361 (in Chinese).

[6]李盾, 何躍龍, 紀楚群. 多體分離數(shù)值模擬研究與應用[C]//北京力學會第19屆學術年會論文集. 北京： 北京力學會, 2013: 121-122.

LI D, HE Y L, JI C Q. Numerical simulation research and application of multi-body separation[C]//The Beijing Society of Theoretical and Applied Mechanics 19th Annual Conference Proceedings. Beijing: The Beijing Society of Theoretical and Applied Mechanics, 2013: 121-122 (in Chinese).

[7]何躍龍, 李盾, 劉曉文. 非結構直角網(wǎng)格動網(wǎng)格方法[C]//北京力學會第18屆學術年會論文集. 北京： 北京力學會, 2012: 19-20.

HE Y L, LI D, LIU X W. Dynamic unstructured cartesian grid method[C]//The Beijing Society of Theoretical and Applied Mechanics 18th Annual Conference Proceedings. Beijing: The Beijing Society of Theoretical and Applied Mechanics, 2012: 19-20 (in Chinese).

[8]劉君, 白曉征, 郭正. 非結構動網(wǎng)格計算方法-及其在包含運動界面的流場模擬中的應用[M]. 長沙: 國防科技大學出版社, 2009.

LIU J, BAI X Z, GUO Z. Dynamic unstructured grid method and its application in simulation of unsteady flows involving moving boundaries[M]. Changsha: National University of Defense Technology Press, 2009 (in Chinese).

[9]BATINA J T. Unsteady euler airfoil solutions using unstructured dynamic meshes[J]. AIAA Journal, 1990, 28(8): 1381-1388.

[10]FARHAT C, DEGAND C, KOOBUS B, et al. Torsional springs for two-dimensional dynamic unstructured fluid meshes[J]. Computer Methods in Applied Mechanics and Engineering, 1998, 163(1): 231-245.

[11]DEGAND C, FARHAT C. A three-dimensional torsional spring analogy method for unstructured dynamic meshes[J]. Computers & Structures, 2002, 80(3): 305-316.

[12]BLOM F J. Considerations on the spring analogy[J]. International Journal for Numerical Methods in Fluids, 2000, 32(6): 647-668.

[13]郭正, 劉君, 瞿章華. 非結構動網(wǎng)格在三維可動邊界問題中的應用[J]. 力學學報, 2003, 35(2): 140-146.

GUO Z, LIU J, QU Z H. Dynamic unstructured grid method with applications to 3D unsteady flows involving moving boundaries[J]. Acta Mechanica Sinica, 2003, 35(2): 140-146 (in Chinese).

[14]BOTTASSO C L, DETOMI D, SERRA R. The ball-vertex method: A new simple spring analogy method for unstructured dynamic meshes[J]. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39): 4244-4264.

[15]ACIKGOZ N, BOTTASSO C L. A unified approach to the deformation of simplicial and non-simplicial meshes in two and three dimensions with guaranteed validity[J]. Computers & Structures, 2007, 85(11): 944-954.

[16]TEZDIUAR T E, BEHR M, MITTAL S, et al. Computation of unsteady incompressible flows with the finite element methods: Space-time formulations, iterative strategies and massively parallel implementations[J]. Asme Pressure Vessels Piping Div Publ PVP, 1992, 246(1): 7-24.

[17]JOHNSON A A, TEZDUYAR T. Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces[J]. Computer Methods in Applied Mechanics and Engineering, 1994, 119(1-2): 73-94.

[18]JOHNSON A A, TEZDUYAR T E. Advanced mesh generation and update methods for 3D flow simulations[J]. Computational Mechanics, 1999, 23(2): 130-143.

[19]STEIN K, TEZDUYAR T, BENNEY R. Mesh moving techniques for fluid-structure interactions with large displacements[J]. Journal of Applied Mechanics, 2003, 70(1): 58-63.

[20]STEIN K, TEZDUYAR T E, BENNEY R. Automatic mesh update with the solid-extension mesh moving technique[J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193(21): 2019-2032.

[21]CHIANDUSSI G, BUGEDA G, ONATE E. A simple method for automatic update of finite element meshes[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 16(1): 1-19.

[22]XU Z, ACCORSI M. Finite element mesh update methods for &uid-structure interaction simulations[J]. Finite Elements in Analysis and Design, 2004, 40(9): 1259-1269.

[23]SMITH R W. A PDE-based mesh update method for moving and deforming high reynolds number meshes: AIAA-2011-0472[R]. Reston: AIAA, 2011.

[24]SMITH R W, WRIGHT J A. A Classical elasticity-based mesh update method for moving and deforming meshes: AIAA-2010-0164[R]. Reston: AIAA, 2010.

[25]YANG Z, MAVRIPLIS D J. Unstructured dynamic meshes with higher-order time integration schemes for the unsteady Navier-Stokes equations: AIAA-2005-1222[R]. Reston: AIAA, 2005.

[26]KARMAN S L, ANDERSON W K, SAHASRABUDHE M. Mesh generation using unstructured computational meshes and elliptic partial differential equation smoothing[J]. AIAA Journal, 2006, 44(6): 1277-1286.

[27]KENNON S R, MEYERING J M, BERRY C W. Geometry based Delaunay tetrahedralization and mesh movement strategies for multi-body CFD: AIAA-1992-4575[R]. Reston: AIAA, 1992.

[28]陳炎, 曹樹良, 梁開洪, 等. 基于溫度體模型的動網(wǎng)格生成方法及在流固耦合振動中的應用[J]. 振動與沖擊, 2010, 29(4): 1-5.

CHEN Y, CAO S L, LIANG K H, et al. A new dynamic grids based on temperature analogy and its application in vibration engineering with fluid-solid interaction[J]. Journal of Vibration and Shock, 2010, 29(4): 1-5 (in Chinese).

[29]陳炎, 曹樹良, 梁開洪, 等. 溫度體動網(wǎng)格模型中控制參數(shù)的研究[J]. 計算物理, 2010, 27(3): 396-402.

CHEN Y, CAO S L, LIANG K H, et al. Parameter control in temperature analogy method[J]. Chinese Journal of Computational Phyfics, 2010, 27(3): 396-402 (in Chinese).

[30]陳炎, 張勤昭, 曹樹良. 溫度體動網(wǎng)格方法的旋轉變形能力[J]. 排灌機械工程學報, 2012, 30(4): 447-451.

CHEN Y, ZHANG Q Z, CAO S L. Rotational deformability of temperature analogy method[J]. Journal of Drainage and Irrigation Machinery Engineering, 2012, 30(4): 447-451 (in Chinese).

[31]陳炎, 張勤昭, 曹樹良, 等. 基準溫度分布動網(wǎng)格生成方法的研究及應用[J]. 北京理工大學學報, 2012, 32(9): 900-904.

CHEN Y, ZHANG Q Z, CAO S L, et al. A new method of dynamic grid generation based on reference temperature distribution[J]. Transactions of Beijing Institute of Technology, 2012, 32(9): 900-904 (in Chinese).

[32]THOMPSON J F, WARSI Z U A, MASTIN C W. Numerical grid generation: Foundations and applications[M]. Amsterdam: North-holland, 1985.

[33]L?HNER R, YANG C. Improved ALE mesh velocities for moving bodies[J]. Communications in Numerical Methods in Engineering, 1996, 12(10): 599-608.

[34]LOMTEV I, KIRBY R M, KARNIADAKIS G E. A discontinuous Galerkin ALE method for compressible viscous flows in moving domains[J]. Journal of Computational Physics, 1999, 155(1): 128-159.

[35]CALDERER R, MASUD A. A multiscale stabilized ALE formulation for incompressible flows with moving boundaries[J]. Computational Mechanics, 2010, 46(1): 185-197.

[36]HUSSAIN M, ABID M, AHMAD M, et al. A parallel implementation of ALE moving mesh technique for FSI problems using openMP[J]. International Journal of Parallel Programming, 2011, 39(6): 717-745.

[37]MASUD A, BHANABHAGVANWALA M, KHURRAM R A. An adaptive mesh rezoning scheme for moving boundary flows and fluid-structure interaction[J]. Computers & Fluids, 2007, 36(1): 77-91.

[38]MASUD A, HUGHESB T J R. A space-time Galerkin/least-squares finite element formulation of the Navier-Stokes equations for moving domain problems[J]. Computer Methods in Applied Mechanics and Engineering, 1997, 146(1): 91-126.

[39]HELENBROOK B T. Mesh deformation using the biharmonic operator[J]. International Journal for Numerical Methods in Engineering, 2003, 56(7): 1007-1021.

[40]WANG Q, HU R. Adjoint-based optimal variable stiffness mesh deformation strategy based on bi-elliptics equations[J]. International Journal for Numerical Methods in Engineering, 2012, 90(5): 659-670.

[41]YIGIT S, FER M S, HECK M. Grid movement techniques and their influence on laminar fluid-structure interaction computations[J]. Journal of Fluids and Structures, 2008, 24(6): 819-832.

[42]ALLEN C B. Parallel flow-solver and mesh motion scheme for forward flight rotor simulation: AIAA-2006-3476 [R]. Reston: AIAA, 2006.

[43]LIEFVENDAHL M, TR?ENG C. Deformation and regeneration of the computational grid for CFD with moving boundaries： AIAA-2005-1090[R]. Reston: AIAA, 2005.

[44]PERSSON P O, BONET J, PERAIRE J. Discontinuous Galerkin solution of the Navier-Stokes equations on deformable domains[J]. Computer Methods in Applied Mechanics and Engineering, 2009, 198(17): 1585-1595.

[45]SHEPARD D. A two-dimensional interpolation function for irregularly-spaced data [C]//Proceedings of the 1968 23rd ACM National Conference. New York: ACM, 1968: 517-524.

[46]WITTEVEEN J A S. Explicit and robust inverse distance weighting mesh deformation for CFD： AIAA-2010-0165 [R]. Reston: AIAA, 2010.

[47]WITTEVEEN J A S, BIJL H. Explicit mesh deformation using inverse distance weighting interpolation: AIAA-2009-3936[R]. Reston: AIAA, 2009.

[48]LUKE E, COLLINS E, BLADES E. A fast mesh deformation method using explicit interpolation[J]. Journal of Computational Physics, 2012, 231(2): 586-601.

[49]BUHMANN M. Radial basis functions[M]. Cambridge: Cambridge University Press, 2005.

[50]BOER A D, SCHOOT M S V D, BIJL H. Mesh deformation based on radial basis function interpolation[J]. Computers & Structures, 2007, 85(11): 784-795.

[51]RENDALL T C S, ALLEN C B. Efficient mesh motion using radial basis functions with data reduction algorithms[J]. Journal of Computational Physics, 2009, 228(17): 6231-6249.

[52]RENDALL T C S, ALLEN C B. Reduced surface point selection options for efficient mesh deformation using radial basis functions[J]. Journal of Computational Physics, 2010, 229(8): 2810-2820.

[53]SHENG C, ALLEN C B. Efficient mesh deformation using radial basis functions on unstructured meshes[J]. AIAA Journal, 2013, 51(3): 707-720.

[54]BOS F M, OUDHEUSDEN B W, BIJL H. Radial basis function based mesh deformation applied to simulation of flow around flapping wings[J]. Computers & Fluids, 2013, 79(6): 167-177.

[55]LIU X, QIN N, XIA H. Fast dynamic grid deformation based on Delaunay graph mapping[J]. Journal of Computational Physics, 2006, 211(2): 405-423.

[56]肖天航, 昂海松, 仝超. 大幅運動復雜構型撲翼動態(tài)網(wǎng)格生成的一種新方法[J]. 航空學報, 2008, 29(1): 41-48.

XIAO T H, ANG H S, TONG C. A new dynamic mesh generation method for large movements of flapping-wings with complex geometries[J]. Acta Aeronautica et Astronautica Sinica, 2008, 29(1): 41-48 (in Chinese).

[57]周璇, 李水鄉(xiāng), 陳斌. 非結構動網(wǎng)格生成的彈簧-插值聯(lián)合方法[J]. 航空學報, 2010, 31(7): 1389-1395.

ZHOU X, LI S X, CHEN B. Spring-interpolation approach for generating unstructured dynamic meshes[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(7): 1389-1395 (in Chinese).

[58]LEFRAN?OIS E. A simple mesh defomation technique for fluid-structure interaction based on a submesh approach[J]. International Journal for Numerical Methods in Engineering, 2008, 75(9): 1085-1101.

[59]ALAUZET F D R. A changing-topology moving mesh technique for large displacements[J]. Engineering with Computers, 2014, 30(2): 175-200.

[60]OLIVIER G, ALAUZET F. A new changing-topology ALE scheme for moving mesh unsteady simulations: AIAA-2011-0474[R]. Reston: AIAA, 2011.

[61]GUARDONE A, ISOLA D, QUARANTA G. Arbitrary lagrangian eulerian formulation for two-dimensional flows using dynamic meshes with edge swapping[J]. Journal of Computational Physics, 2011, 230(20): 7706-7722.

[62]ISOLA D. An interpolation-free two-dimensional conservative ALE scheme over adaptive unstructured grids for rotorcraft aerocraft aerodynamics[D]. Milano: Politecnico Di Milano, 2012.

[63]WANG L, PERSSONY P O. A discontinuous galerkin method for the Navier-Stokes equations on deforming domains using unstructured moving space-time meshes: AIAA-2013-2833[R]. Reston: AIAA, 2013.

[64]GOPALAKRISHNAN P, TAFTI D K. A parallel multiblock boundary fitted dynamic mesh solver for simulating flows with complex boundary movement: AIAA-2008-4142[R]. Reston: AIAA, 2008.

[65]ZHOU X, LI S. A new mesh deformation method based on disk relaxation algorithm with pre-displacement and post-smoothing[J]. Journal of Computational Physics, 2013, 235(2): 199-215.

[66]LIU Y, GUO Z, LIU J. RBFs-MSA hybrid method for mesh deformation[J]. Chinese Journal of Aeronautics, 2012, 25(4): 500-507.

[67]ANDERSON J M, STREITLIEN K, BARRETT D S, et al. Oscillating foils of high propulsive efficiency[J]. Journal of Fluid Mechanics, 1998, 360(1): 41-72.

[68]LIAN Y, SHYY W. Aerodynamics of low reynolds number plunging airfoil under gusty environment: AIAA-2007-0071[R]. Reston: AIAA, 2007.

[69]YOUNG J. Numerical simulation of the unsteady aerodynamics of flapping airfoils[D]. Canberra: University of New South Wales, 2005.

[70]THOMAS P D, LOMBARD C K. The geometric conservation law-a link beween finite-difference and finite-volume methods of flow computation on moving grids: AIAA-1978-1208[R]. Reston: AIAA, 1978.

[71]THOMAS P D, LOMBARD C K. Geometric conservation law and its application to flow computations on moving grids[J]. AIAA Journal, 1979, 17(10): 1030-1037.

[72]GUILLARD H, FARHAT C. On the significance of the geometric conservation law for flow computations on moving meshes[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 190(11): 1467-1482.

[73]GEUZAINE P, GRANDMONT C, FARHAT C. Design and analysis of ALE schemes with provable second-order time-accuracy for inviscid and viscous flow simulations[J]. Journal of Computational Physics, 2003, 191(1): 206-227.

[74]FARHAT C, GEUZAINE P. Design and analysis of robust ALE time-integrators for the solution of unsteady flow problems on moving grids[J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193(39): 4073-4095.

[75]MAVRIPLIS D J, YANG Z. Construction of the discrete geometric conservation law for high-order time-accurate simulations on dynamic meshes[J]. Journal of Computational Physics, 2006, 213(2): 557-573.

[76]HYUNG T A, YANNIS K. Strong coupled flow/structure interaction with a geometrically conservative ALE scheme on general hybrid meshes[J]. Journal of Computational Physics, 2006, 219(2): 671-696.

[77]ZHANG L P, WANG Z J. A block LU-SGS implicit dual time-stepping algorithm for hybrid dynamic meshes[J]. Computers & Fluids, 2004, 33(7): 891-916.

[78]KAMAKOTI R, SHYY W. Evaluation of geometric conservation law using pressure-based fluid solver and moving grid technique[J]. International Journal of Numerical Methods for Heat & Fluid Flow, 2004, 14(7): 851-865.

[79]BOFFI D, GASTALDI L. Stability and geometric conservation laws for ALE formulations[J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193(42): 4717-4739.

[80]MASUD A. Effects of mesh motion on the stability and convergence of ALE based formulations for moving boundary flows[J]. Computational Mechanics, 2006, 38(4-5): 430-439.

[81]CHANG X H, MA R, ZHANG L P, et al. Further study on the geometric conservation law for finite volume method on dynamic unstructured mesh[J]. Computers & Fluids, 2015, 120: 98-110.

[82]安效民, 徐敏, 陳士櫓. 二階時間精度的CFD/CSD耦合算法研究[J]. 空氣動力學學報, 2009, 27(5): 547-552.

AN X M, XU M, CHEN S L. Analysis for second order time accurate CFD/CSD coupled algorithms[J]. Acta Aerodynamica Sinica, 2009, 27(5): 547-552 (in Chinese).

劉君男， 博士， 教授， 博士生導師。主要研究方向： 空氣動力學。

Tel: 0411-84707176

E-mail: liujun65@dlut.edu.cn

劉瑜男， 博士， 講師。主要研究方向： 燃燒數(shù)值模擬。

E-mail: liuyu@nudt.edu.cn

陳澤棟男， 博士研究生。主要研究方向： 計算流體力學， 空氣動力學。

E-mail: chenzd_dut@163.com

Unstructured deforming mesh and discrete geometricconservation law

LIU Jun1,*, LIU Yu2, CHEN Zedong3

1. School of Aeronautics and Astronautics, Dalian University of Technology, Dalian116024, China 2. Department of Space Equipment, Academy of Equipment, Beijing101416, China 3. Department of Engineering Mechanics, Dalian University of Technology, Dalian116024, China

A popular method for simulating unsteady flow of fluid-structure interaction or multi-body separation is finite volume method based on arbitrary Lagrangian-Eulerian (ALE) equations, which involves mesh deformation and discrete geometric conservation law. The fundamental principle, state of the art and boundaries of validity of mesh deformation are reviewed following 5 categories of constructing ideas, e.g., physics analogy, ellipse smoothing, interpolation, moving submesh approach (MSA) and hybrid method. A hybrid method which combines the benefits of MSA and radial basis functions (RBFs) interpolation is proved to be robust and efficient via several numerical examples. After the concept of discrete geometric conservation law (DGCL) is introduced, the intrinsic mechanism of DGCL is analyzed through 2D model, which is the inequality between volume increment and the volume sweeping by the surfaces enclosing of mesh cell. The different implementations of DGCL which could be mainly categorized as area correction, area correction via assuming velocity of surface, velocity correction and volume correction are surveyed and their range of application and the existing problems are discussed. We found that the proposed volume correction method can satisfy the fluid-structure interface condition, and is also appropriate for multi-step temporal discrete schemes.

unstructured mesh; mesh deformation; unsteady flow; geometric conservation law; interface coupling

2016-01-16； Revised： 2016-02-17； Accepted： 2016-05-03； Published online： 2016-05-1015:25

National Natural Science Foundation of China (91541117)

. Tel.： 0411-84707176E-mail: liujun65@dlut.edu.cn

2016-01-16； 退修日期： 2016-02-17； 錄用日期： 2016-05-03；

時間： 2016-05-1015:25

www.cnki.net/kcms/detail/11.1929.V.20160510.1525.004.html

國家自然科學基金 (91541117)

.Tel.： 0411-84707176E-mail: liujun65@dlut.edu.cn

10.7527/S1000-6893.2016.0141

V211.3

A

1000-6893(2016)08-2395-13

引用格式： 劉君， 劉瑜， 陳澤棟． 非結構變形網(wǎng)格和離散幾何守恒律[J]． 航空學報, 2016, 37(8): 2395-2407. LIU J, LIU Y, CHEN Z D. Unstructured deforming mesh and discrete geometric conservation law[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2395-2407.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

URL： www.cnki.net/kcms/detail/11.1929.V.20160510.1525.004.html