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

劉君， 劉瑜， 陳澤棟

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

?

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

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

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

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

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

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

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

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

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

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

1　非結(jié)構(gòu)動網(wǎng)格技術(shù)

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

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

1.1物理比擬方法

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

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

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

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

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

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

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

1.3插值方法

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

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

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

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

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

1.5變拓?fù)浞椒?/p>

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

連接關(guān)系確定后計(jì)算量很少，變拓?fù)浞椒ǖ挠?jì)算效率最高，但是數(shù)據(jù)結(jié)構(gòu)變化較頻繁，需要求解器進(jìn)行改造，實(shí)現(xiàn)較為復(fù)雜，目前應(yīng)用范圍較窄。

1.6混合方法

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

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

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

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

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

2) 計(jì)算網(wǎng)格內(nèi)點(diǎn)在運(yùn)動子網(wǎng)格單元中的相對面積或體積坐標(biāo)(映射關(guān)系)。

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

4) 利用徑向基函數(shù)插值得到運(yùn)動子網(wǎng)格內(nèi)點(diǎn)的坐標(biāo)。

5) 利用映射關(guān)系得到計(jì)算網(wǎng)格內(nèi)點(diǎn)坐標(biāo)。

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

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

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

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

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

α0=arctan(2kh0)-θ0

(2)

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

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

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

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

3　離散幾何守恒律

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

(3)

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

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

(4)

(5)

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

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

3.1面積修正法

3.2給定速度的面積修正法

(6)

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

3.3速度修正法

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

(7)

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

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

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

(8)

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

3.4體積修正法

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

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

(9)

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

(10)

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

(11)

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

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

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

圖5　計(jì)算誤差Fig.5　Calculation errors

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

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

(12)

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

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

4　界面耦合算法

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

(13)

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

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

(14)

式中：vs為界面上的固體運(yùn)動速度；rs為界面上的固體運(yùn)動位移；m為質(zhì)量；Fair為氣動力。

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

5　結(jié)　論

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

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

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

[1]楊永健, 張魯民. 應(yīng)用重疊網(wǎng)格技術(shù)求解復(fù)雜組合體無粘流場[J]. 空氣動力學(xué)報(bào), 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)化技術(shù)的改進(jìn)與應(yīng)用[J]. 航空學(xué)報(bào), 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]. 空氣動力學(xué)學(xué)報(bào), 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]. 航空學(xué)報(bào), 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]劉學(xué)強(qiáng), 伍貽兆, 夏健. 多重網(wǎng)格法在非結(jié)構(gòu)網(wǎng)格中的應(yīng)用[J]. 計(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]李盾, 何躍龍, 紀(jì)楚群. 多體分離數(shù)值模擬研究與應(yīng)用[C]//北京力學(xué)會第19屆學(xué)術(shù)年會論文集. 北京： 北京力學(xué)會, 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]何躍龍, 李盾, 劉曉文. 非結(jié)構(gòu)直角網(wǎng)格動網(wǎng)格方法[C]//北京力學(xué)會第18屆學(xué)術(shù)年會論文集. 北京： 北京力學(xué)會, 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]劉君, 白曉征, 郭正. 非結(jié)構(gòu)動網(wǎng)格計(jì)算方法-及其在包含運(yùn)動界面的流場模擬中的應(yīng)用[M]. 長沙: 國防科技大學(xué)出版社, 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]郭正, 劉君, 瞿章華. 非結(jié)構(gòu)動網(wǎng)格在三維可動邊界問題中的應(yīng)用[J]. 力學(xué)學(xué)報(bào), 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)格生成方法及在流固耦合振動中的應(yī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]. 計(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)格方法的旋轉(zhuǎn)變形能力[J]. 排灌機(jī)械工程學(xué)報(bào), 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]陳炎, 張勤昭, 曹樹良, 等. 基準(zhǔn)溫度分布動網(wǎng)格生成方法的研究及應(yīng)用[J]. 北京理工大學(xué)學(xué)報(bào), 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]肖天航, 昂海松, 仝超. 大幅運(yùn)動復(fù)雜構(gòu)型撲翼動態(tài)網(wǎng)格生成的一種新方法[J]. 航空學(xué)報(bào), 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), 陳斌. 非結(jié)構(gòu)動網(wǎng)格生成的彈簧-插值聯(lián)合方法[J]. 航空學(xué)報(bào), 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]安效民, 徐敏, 陳士櫓. 二階時(shí)間精度的CFD/CSD耦合算法研究[J]. 空氣動力學(xué)學(xué)報(bào), 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).

劉君男， 博士， 教授， 博士生導(dǎo)師。主要研究方向： 空氣動力學(xué)。

Tel: 0411-84707176

E-mail: liujun65@dlut.edu.cn

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

E-mail: liuyu@nudt.edu.cn

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

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；

時(shí)間： 2016-05-1015:25

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

國家自然科學(xué)基金 (91541117)

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

10.7527/S1000-6893.2016.0141

V211.3

A

1000-6893(2016)08-2395-13

引用格式： 劉君， 劉瑜， 陳澤棟． 非結(jié)構(gòu)變形網(wǎng)格和離散幾何守恒律[J]． 航空學(xué)報(bào), 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