張來平山f 鄧小剛 M 張涵信 2

　　1空氣動力學國家重點實驗室，綿陽 621∞0

　　2 中國空氣動力研究與發展中心，綿陽 621儀)()

　　摘 要對應用于飛行器非定常運動的數值計算方法(包括動態網格技術和相應的數值離散格式)進行了綜 述.根據網格拓撲結構的不同，重點論述了基于結構網格的非定常計算方法和基于非結構/混合阿格的非定常 計算方法，比較了各種方法的優缺點.在基于結構網格的非定常計算方法中，重點介紹了剛性運動網格技術、 超限插值動態網格技術、重疊動網格技術、滑移動網格技術等動態結構網格生成方法，同時介紹了慣性系和非慣性系下的控制方程，討論了非定常時間離散方法、動網格計算的幾何守恒律等問題.在基于非結構/混合同 格的非定常計算方法中，重點介紹了重疊非結構動網格技術、重構非結構動網格技術、變形非結構動網格技術以及變形/重構稿合動態混合網格技術等方法，以及相應的計算格式，包括非定常時間離散、幾何守恒律計算 方法、可壓縮和不可壓縮非定常流動的計算方法、各種加速收斂技術等.在介紹國內外進展的同時，介紹了作 者在動態混合網格生成技術和相應的非定常方法方面的研究與應用工作.

　　關鍵詞：中級職稱論文,非定常運動,動態網格技術,非定常計算方法

　　1 引言

　　對速度、機動性與效率的追求，一直主導著軍 事飛行器的發展.無論是未來的戰略、戰術導彈， 還是新一代的戰斗機，都需要具備良好的機動性 和敏捷性，尤其要求在快速機動時，能夠利用流動 特征對飛行器實施高效的主動控制.為了實現飛 行器的可控快速機動，首先必須弄清飛行器的動 態氣動特性，以及與這些動態特性相對應的非定 常流動機理

　　從飛行器的運動方式和流動特征來看，非定 常流動問題可以分為以下 3 類: (1) 物體靜止而 流動本身為非定常的流動問題，如大攻角飛行的 細長體背風區的分離流動等; (2) 單個物體作剛性 運動的非定常流動問題，如飛行器的俯仰、搖滾及 其藕合運動等 j (3) 多體作相對運動或變形運動的 非定常問題，如子母彈分離、飛機外掛物投放、機 翼的氣動彈性振動、魚類的擺動、昆蟲和鳥類的撲 動等等.

　　收稿日期 :2∞9-9-4，修回日期:2010-4-23

　　解決以上非定常問題的有效途徑之一是計算 流體力學 (CFD). 隨著計算機技術的迅猛發展， CFD 己經成為一種有效的流動機理分析和氣動設 計手段.但是，對于上述的非定常流動問題，高效 高精度的非定常計算方法仍然是當今 CFD 研究 的熱點和難點.

　　為了有效地開展非定常計算方法的研究，本 文對當前國內外的非定常計算方法(包括網格技 術和相應的非定常計算格式)進行了綜述.根據網 格拓撲結構(結構網格、非結構網格和混合網格) 的不同，將非定常計算方法分為兩大類，即基于結 構網格的非定常計算方法和基于非結構/混合網 格的非定常計算方法.在評述過程中，同時介紹 了一些作者在非定常計算方法方面的研究及應用 工作.

　　2 非定常計算方法的分類

　　非定常計算方法可以按多種方式進行分類，

　　·國家自然科學基金 (108η023) 和國家重點基礎研究發展計劃 (2∞9CB72擁02) 資助項目

　　t E-mail: zha噸lp..c町dc@126.∞m

　　第 4 期 張米平等:動網格生成技術及非定常計算方法進展綜述 425

　　例如按時間離散方法，可以分為 :顯示 Runge-Kutta 方法、單時間步法、雙時間步方法、隱式迭代法等 ; 再比如，按照選取的坐標系的不同，可以分為慣性 系下的非定常計算方法和非慣性系下的非定常計 算方法.本文以計算方法采用的網格技術為基礎， 對非定常計算方法進行分類，因為基于不同的網 格技術，需要采用不同的數值計算方法.

　　眾所周知，數值計算的第 1 步是生成合適的 計算網格，即用合適的計算網格離散流場.就目前 國內外的網格技術發展情況來看[習，網格生成技 術大致可以分為 3 類，即結構網格、非結構網格和 混合網格.混合網格綜合了結構網格和非結構網 格的優勢，代表了當前和未來網格技術發展的趨 勢 (2) 就基于溫合網格的計算方法而言，目前大多 將混合網格數據結構轉化為非結構網格的形式進 行整體處理這樣使得整個流場的計算方法統一， 便于程序設計和調試.因此，就非定常計算方法而 言，將其細分為兩大類，即:

　　(1) 基于結構網格的非定常計算方法;

　　(2) 基于非結構/混合網格的非定常計算方法.

　　3 基于結構網格的非定常計算方法

　　3.1 適用于非定常計算的結構動網格生成技術 根據非定常運動方式的不同，可以采用不同

　　的結構網格.對于前述的第 1 類非定常問題，即物 體靜止而流動本身為非定常的流動問題，靜止的 剛性網格就可以滿足要求當然，由于所面對的 外形越來越復雜，傳統的統 一貼體結構網格(四­ fìed body-fìtted grids) 已無法生成合適的網格，為此 CFD 工作者發展了多塊對接結構網格 (multiblock structured grids) [3，4J 、多塊搭接結構網格(patched structured grids)[5，6] 和重疊結構網格技術 (overlap­ pingjchimera structured grids) [7，8) 等.

　　對于第 2 類非定常問題，仍然可以采用剛性 網格.根據坐標系選取方式的不同，可以采用慣性 系和非慣性系兩種方法進行非定常計算.對于第 3 類非定常問題，涉及到物體的相對運動或變形，因 此必須在每個時間步更新網格，由此動網格生成 技術成為非定常計算的關鍵技術之一，而且它是 實現工程應用中最主要的"瓶頸"問題.目前，剛 性運動網格技術 [9，10) 、超限插值(transfìnite inter­ polation) 動網格生成技術 [ll 14J 、重疊結構動網格

　　技術 [15 18] 、滑移結構動網格技術[19，20J 是幾種常 用的方法.

　　3.1.1 剛性運動網格技術 剛性運動網格技術的基本思想是令計算網格

　　隨物體一起作剛體運動.這一方法的優點是:在整 個非定常運動過程中，計算網格無須重新生成，可 以根據運動方式直接給出，因此其計算量小，并可 保持初始網格的質量.但是，這種方法僅適用于單 個剛性物體的非定常運動，對于變形體或者多體 的相對運動等復雜問題，這種方法已經不再適用. 另一方面，由于網格隨體運動，使遠場邊界的位移 和速度很大，導致運動遠場邊界的處理不易.因此 這種方法在簡單外形的非定常運動數值模擬中采 用得較多 [9，lJO

　　3.1.2 超限插值動網格生成技術 超限插值動網格生成方法的基本思想是令外

　　邊界保持靜止，物面邊界由物體運動規律或運動 方程得到，內場網格由超限插值的方法代數生成， 其計算量也較小，能夠生成相對復雜的動態網格， 但不易保證網格品質，尤其是不能保證物面網格 的正交性為此，文獻 [14] 提出了一種加權技術將 剛性運動網格技術和超限插值動網格技術結合起 來，實現了飛船等外形的動態結構網格生成.

　　加權超限插值動網格生成方法的基本思路是: 首先，由初始網格 xn 生成剛性運動網格 xmf1，再 由 η 時刻的開邊界和 n+1 時刻的物面邊界代數 插值得到 n+1 時刻的網格 xren ，最后由兩者加權 得到新時刻的動網格

　　xn+1 = (1-ω)xrefl +ωxref2，0 ω 三 1 (1)

　　只要適當選擇加權因子 ω，就既可保證物面 附近網格質量，義可使外邊界保持靜止而易于邊 界處理.文獻 [14] 構造了如下的加權因子，網格生 成實例表明效果較好.令 j 為法向網格指標，則有

　　3.1.3 重疊結構動網格技術

　　重疊網格技術始見于 1983 年 Steg町等 [7] 的 開創性 E作.重番網格技術的基本思想是在計算 域的各個子域采用阿域(亙疊部分)共享的方法來 實現信息交 換，而不是 采用邊界共享的方法 ，從而 大大減輕了子域網格生成的難度，而且能夠保證 子域的網格品質.重疊網格技術在復雜外形的數 值模擬中己經得到了廣泛的應用 [8]

　　由于重疊結構網格在處理簡單的單體組合構 刑方面的優勢，其被許多 CFD 工作者推廣應用 于多體分離非非定常問題的數值模擬 .例如在美 國國防部資助的外掛物投放 ACFD CHANLLENG I"，-，IV 項日中，許多參與該項日的研究人員采用了 重替 網格技術 [15，16J 在其他方面，也 有很好 的重 疊網格的應用實例 [17，叫在國內，張玉東等[21]利 用重疊結構動網格技術數值模擬了子母彈分離過

　　和;李亭鶴 r.r; [22J 利用多塊豆費結構動網格技術和 準定常計算方法，數值模擬了子母彈拋殼過秤:，楊 明智 [23] 等利用豆替結構動網格數值模擬了直升

　　機外掛物投放的過和.國 2 顯示了文獻 [23] 中采 用的直升機、懸翼和外掛物的重疊網格拓撲結構 和l直升機周 圍的結構網格.

　　盡管重疊結構動網格技術在實際工和應用中得到了廣泛應用，也取得了巨大的成功. 但是?

　　;在非定常運動的每一新的時間層上，重疊網格不

　　僅需要更新各子域的網格而且還需要對于域的 重香阪網格進行插值平 n更新，從而導致計算量增 大.對于流場中存在強間斷的問題，如果強間斷 穿越豆疊慶，則會導致子域邊界間的插值誤差， 由此影響到非定常問 題的計 算精度.這種插值誤 差 的長時間積累將導致非定常 計算誤差的進 一步 放大.

　　第 4 期 張米平等功網格生成技術及非定常計算方法進展綜述 427

　　國 2 重疊結構動網格生成豆豆例 [231

　　3.1.4 滑移結構動網格技術 由于重疊網格需要在每個時間層內的子選代

　　中進行重疊阪插值，而且在每個時間層之間，需

　　要重新搜尋插值關系，一方面會帶來插值誤 差，另 一方面由于搜尋插值關系，也會降低計算效率.為

　　此， CFD 工作基于搭接結構網格技術，發展了搭 接結構動網格技術，即所謂滑移結構動網格技術 ， 目前在一些大型 CFD 軟件中集成了這種網格技

　　術[凹，201

　　滑移結構動網格技術的基本思想是在運動部 件的運動軌跡周圍預先劃分出一個滑移子域.在 滑移子域和以外的區域分別生成多塊結構網格 j 在滑移子域與其他區域的交界面處?利用搭接邊 界條件與其他|天域對接從而實現整體流場的計 算.滑移網格技術在旋轉部件(如發動機內流、直 升機懸翼、螺旋槳等)、列車交會、列車過隧道等問 題中應用廣泛.國 3 顯示了從文獻 [19J 中摘錄的滑 移結構動網格實例.

　　國 3 消移結構動網格頭'例 [191

　　3.2 基于結構動網格的非定常計算格式

　　3.2.1 非慣性系卡的非定常計算方法

　　正如 3.1 節中所述，征剛性結構網格卡，可以 采用不同的參考坐標系描述控制流動方和. 如果 為了消除動網格計算所帶來的諸多問題(如兒何 守恒率、動邊界條件等)，可以采用非慣性系卡的 流動方和. 例如，對于俯仰運動問題，假設坐標系 繞 z 軸在俯仰平面內作單自由度旋轉運動，角速

　　度矢量為 ω(t) = -f2k ，記 V 為流場中某 一點的絕 對速度(在慣性坐標系卡觀察的速度)， Vr 為 該點 相對于動坐標系(國連在運動的物體上隨物體 一

　　起轉動)的速度 ，r 為該點到動坐標系原 點的向徑

　　(見國 4) ，應有V= Vr 十 ω X r . 經過不復雜的推 導，就可以得到俯仰運動 非慣性笛卡兒坐標系 下 以 (ρ，Ur ，Vr ，Wr， p) T 為原始變量 的三維非定常無量 綱化的 Navier-Stokes 方和

　　θQ ， θE θFθG 一一一一一一一 θt ' θX δY θz

　　1 (θEv θF v θGv \ ，

　　Re∞ \θX θν l θz ) ， - (4)

　　上式中各物理量的具體表達式可以參見文獻 [14]，

　　。

　　其中的 S 是岡為坐標變換而 引入的源項

　　ρ( f22x - 2 f2vr - f2ν)

　　8 = 1ρ。( f22y + 2 f2ur + βX) I (5)

　　ρ[ f22(XU r + YVr ) + 口(XVr -yur)J

　　利用非慣性系的流動 方程，可以降低計算對

　　動網格的依賴性，傳 統 的定常流網格生成方法可

　　428

　　力 學 進 展

　　2010 年第 40 卷

　　以直接應用于生成非慣性系的非定常計算網格. 但是，這種方法具有較大的局限性.第一，在非慣 性系下的方程非常復雜，尤其是在考慮多自由度 運動的情況下會更加復雜;第二，在非慣性系下的 邊界條件的實現，同樣必須進行嚴格的坐標系變

　　化的三維非定常 Navier-Stokes 方程可寫成如下守 恒形式

　　一-一-一一一

　　θ(J-1.Q) ， θtt-OS' ，。。­

　　，

　　θt θ'{ .旬 'θ〈

　　1 I θE制。F.， θG.，

　　換，形式復雜;第三，對于源項的處理是一個需要

　　Re \茍:..+茍:..+寫:")

　　(8)

　　研究的問題，處理得不好將影響計算的收斂效率; 第四，非慣性系下得到的計算結果需要經過繁瑣 的坐標變化轉換到慣性系 F進行流場顯示，帶來 諸多不便;第五，對于多體分離問題或變形運動問

　　題，該方法難以滿足要求.鑒于上述問題，目前國

　　值得注意的是，在動網格非定常計算中，變換中必 須保留時間導數項，即

　　p9

　　p9u +kxp

　　際上通行的非定常計算方法絕大部分均采用慣性 系下的流動控制方程進行計算.

　　F = ktQ + kxE +k1lF + kzG=

　　ρ9v +k1lP

　　ρ0ω +kzp

　　ρ9h - ktp

　　Yo

　　V∞

　　一一+

　　zO z

　　αt: 平衡攻角

　　8(t): 俯仰振蕩角

　　V=Vr +ω xr

　　α=αt + 8(t)

　　z

　　Fv = kxEv +k1lFv + kzGv =

　　o

　　kxTxx +k1lTX 1l +kzTxz kX T1IX +k1lT1I1I +kzT1Iz

　　。

　　kxTzx + k1lTz1l +kzTzz kxßx +k1lt馬 +kzßz

　　= kt + kxu +k1lv +kzω = kx (u - Xt) +

　　(9a)

　　(9b)

　　怡，y ，z): 非慣性坐標系

　　(xo，YO ，zo): 慣性坐標系

　　圖 4 俯仰運動非慣性坐標系統的定義

　　k1l (v - Yt) +kz (ω - Zt)

　　關于其中符號的含義可以參見文獻 [14).

　　3.2.3 非定常計算的時間離散方法

　　(10)

　　3.2.2 慣性系下的非定常計算方法 在慣性笛卡爾坐標系中，無量綱化的 三維非

　　定常 Navier-Stokes 方程可寫成如下守恒形式

　　- -ô-E+ F+-一ôG =

　　δQ. δ

　　ôt-. + θ'x .-ô-y . ôz

　　1 (δEv . ôF" . δGv 飛

　　Re∞飛 θx ' ôy ， θZ ) (6)

　　在貼體結構網格上進行計算，需要進行貼體坐標 變換，在動網格計算中，一般引入如下與時間相關 的變換

　　(7)

　　e = e(x，y，z; t.) l η=η(x，y ，z; t.) I ç = Ç(x，y ，z;t.) J

　　將物理平面中的運動網格轉換至計算平面中的正 交靜止網格.則在一般靜止曲線坐標系中，無量綱

　　非定常計算的空間離散方法，與定態網格情 況 F的計算方法類似，唯一的不同之處是 :在動網 格情況下，對流通量項中多出了網格運動項，即式 (10)中的問，協和 Zt. 由于網格運動項的出現，要 求在動網格情況下保持幾何守恒率，該問題將在 3.2.4節中討論.

　　非定常計算除了關注空間離散精度外，還要

　　重點關注時間項的離散方法.時間方向的精確離 散，才有可能得到精確的非定常計算結果.另一方 面，因為非定常問題的計算量要比定常問題的計 算量高 1"，，2 個量級，非定常計算的效率也是普遍 關注的問題.因此，時間項的離散就顯得尤其重要.

　　關于非定常計算的時間離散方法，大致可以分為如下幾種:顯式時間推進 [24l ，雙時間步隱式 迭代推進 [25l ，單時間步隱式法代推進[26] ，通量 分裂隱式迭代推進 [14] 以及半隱半顯時間推進方

　　法 [27] 等.

　　第 4 期 張來平等:動網格生成技術及非定常計算方法進展綜述 429

　　J-l +J-1 -'" - .''''

　　守恒方程 (8) 經空間離散后，可寫成如下半離 散形式

　　方程 (15)左端引入虛擬時間 (1') 的導數項，得到

　　θ(J-IQ)

　　，âQ ， T_，3Q_4Qn+Qn-l ，

　　述混合網格技術推廣應用于非定常流動的計算.

　　與結構動網格類似，將非結構網格和混合網 格推廣應用于運動物體非定常運動的方法主要有: 重疊非結構動網格技術 [37，381，重構非結構動網格 技術 l39，40l ，變形非結構動網格技術[41~47!，以及變 形/重構混合網格生成技術問 50J 等.當然，如果 利用非慣性系下的方程，計算網格可以采用剛性 網格，這與結構網格類似，這里不再重述.

　　4.1.1重疊非結構動網格技術

　　正如 3.1.3 節中所述，重疊網格方法最先在 結構網格中得到應用. 然而這種方法在網格間 信息傳遞方面的欠缺始終得不到有效解決 首 先 Chimera-hole 需要在每一個網格區域進行切 割，不可避免地要破壞物面邊界或者外圍不需要 切割的網格;其次插值需要建立在重疊網格洞邊界 (chimera-hole)的每一個邊界點上.這些"挖洞"過 程中存在的問題影響了結構重疊網格的適用范圍.

　　432 力 學 進 展

　　2010 年第 40 卷

　　非結構重疊網格的出現有效地解決了上述問 題.通過使用非結構網格，覆蓋流動區域所需要的 子網格數目可以大大減少，而且很容易將該方法 推廣到處理有相對運動的多體繞流問題下面介 紹 N幽hashi 等問提出來的在非結構重疊動網 格中用重疊網格邊界定義 (intergrid-boundarγdef­ inition) 技術進行數值傳遞的方法.該方法中子網 格間的數值傳遞，主要包括兩個步驟:

　　(1) 網格"挖洞"該步驟首先將每個子網格中 的所有節點分為兩類:活動節點和非活動節點，主

　　網格(背景網格)的重疊區內邊界則由最靠近非活

　　動節點的活動節點構成.所謂"活動"是指在非定 常計算中需要計算的網格節點;而"非活動"是指 與另外的網格域重疊而無需計算的網格節點，這 部分節點可以被"挖"去.而用來區分是否被挖去

　　的判據，則是"距壁面的最小距離飛具體方法是: 首先計算包涵某部件的子網格域內的網格點到該 部件表面的距離缸，和該節點到背景網格中部件

　　表面的距離 d2，比較二者的大小，將更靠近子網格 部件表面的網格點作為活動節點 (dl < d2). 如圖

　　5 所示，虛線表示包涵部件 A 的子網格域，其中點

　　t 距 A 表面的距離較 B 近，因此節點 4 為活動節 點參與計算;而節點 j 則相反，屬于非活動節點 圖 6 '"圖 8 則是根據這一判據所得到的網格. 由于分離物體在很短的時間步內只運動很小的距 離，每一個網格可能只運動到鄰近的三、四層網 格位置，因此在建立起每一時間步內關于節點的 八叉樹數據結構后，可以使用 N-T-N (neighbor­

　　to-neighbor) 查找技術快速確立網格之間的聯

　　系 [37，叫.

　　國 5 洞邊界搜尋與插值盡意圖

　　國 7 背景網格中 的活動部分

　　國 6 背景網格與子網格

　　國 8 子網格中的活動部分

　　第 4 期 張來平等:動網格生成技術及非定常計算方法進展綜述 433

　　(2) 數據插值.在第 1 步中確定了運動子網格 域的邊界點與背景網格域的聯系后，通過背景網 格重疊區域的活動網格得到子網格"外"邊界上 的物理值，反之亦然.

　　雖然重疊非結構動網格具有比重疊結構網格 更好的適應性，但是其仍然需要在每一個時間步 進行數據插值，因此不可避免地會帶來誤差 ，這對 長時間的非定常計算同樣是不利的.

　　4.1.2 重構非結構動網格技術

　　"重構" (static-regridding) 非結構動網格生成 技術，即在每個時間步重新生成一次網格.其優點 是: (1) 概念簡單 (2) 容易處理大變形或大位移

　　問題.其缺點是: (1) 每一步都需要插值，不可避

　　免地引入插值誤差; (2) 每個時間層之間的插值關 系需要重新獲得，由此增大了計算量.這種方法被 Grüber 和 Carstens[39) 用于渦輪發動機葉片的強迫 振動問題.Schulze[40) 采用這種辦法處理了機翼運 動的空氣彈性問題.

　　4.1.3 變形非結構動網格技術

　　"變形" (movi噸:-grid 或 deforming-grid) 非結構 動網格生成技術的基本思想是:在保持網格拓撲 結構不變的情況 F，重新分布網格節點其優點 是能夠保持網格關聯信息不變，避免了每一步的 插值.缺點是:在大尺度運動后網格質量會變的很 差，甚至有可能相交，即出現負體積單元，導致非 定常計算的失敗.非結構網格的變形技術主要有

　　以-1"幾種方法:

　　(1) 彈簧拉伸法

　　在文獻 [41] 中， Frederic 將彈簧拉伸法分為兩 類:點拉伸法 (vertex spring) 和片拉伸法 (segment spring). 點拉伸法最先用于優化初始非結構網格 的生成，其基本思想是將網格節點之間的聯系視 為等強度的彈簧，在給定邊界約束的前提下，彈簧 系統的平衡態即為優化后的網格分布 .在動網格 生成過程中，由于物體運動導致彈簧系統的平衡 態破壞，此時需要通過松弛法代法移動相關的網 格節點，保持彈簧系統的動態平衡.由 Hook 定律 可以得到所有與節點 t 相連的節點 j 對其所施的 力

　　式中 αω 是節點 4 和 j 之間的彈性系數 ，N 是與 4 相鄰的節點數目.為了使系統平衡，每個節點上的 合力必須等于 0，整理之后的法代方程為

　　N

　　Zα432?

　　z?+1= 主卡一 (37)

　　藝 αij

　　方程 (36) 作用于網格內的每一個節點，每法 代一次都使節點進→步趨于平衡.在實際應用過 程中，并不要求整個系統達到真正的平衡，因為彈 簧拉伸只是一種優化網格質量的方法，因此法代 的次數可以不必太多，過多的法代步數會影響計 算效率.在每一次的迭代中，新的節點位置 Xi 實 際上是周圍節點坐標的加權平均，權重是彈性系 數 αij. 一般取何= 1，即假設所有彈簧等強度;有 時根據需要，可以適當調整彈性系數 α妒

　　點拉伸法對于小變形或小位移問題具有良好

　　的計算效率;但是對于變形或位移較大的問題，將 無法得到高質量的計算網格為此， Batina[42) 等提 出了片拉伸法，其與點拉伸法的不同之處在于拉 伸時的平衡長度 (equilibrium length of spring) 不同 (點拉伸法的平衡長度設為 0). 該方法中，彈簧的 平衡長度定義為初始時刻的節點間距，而將 Hook 定律作用于節點的位移量，即

　　N

　　F i = 匯α甘 (ðj - ði) (38)

　　j=l

　　式中 ði ，ðj 分別是節點i，j 的位移.迭代方程即為

　　N

　　Zαijðj n

　　δ?+l= 弓一-UXi (39)

　　藝 αij i=l

　　這里， Diric劃.et 邊界條件被設定為己知的動邊界位 移，而彈簧的彈性系數取為節點間距的倒數

　　αij = I !=:= (甜)

　　\/(Xi -Xj)" + (執 - Yj)"

　　經過方程 (38) 的迭代，點 t 的坐標最終成為

　　z?ew=z?ld+δ?anal(41)

　　N

　　F= 藝創j( 的 -Zi) (36)

　　Batina[42) 首先將該方法用于生成機翼顫振的

　　動態非結構網格. 后來很多學者都采用該方法來

　　434 力 學 進 展

　　2010 年第 40 卷

　　處理動邊界問題:如 Slikkev'四r 等問用它來處理 自由表面問題;H兇san 等問!用它來處理外掛物分 離問題; Blom 和Leyland[45] 用它來處理強迫振動 和流場結構干擾問題; FI町hat 等 [46] 和 Piperno[47] 用它來處理氣動彈性的計算問題.值得注意的是， 彈簧拉伸法不僅適用于非結構網格，也適用于結 構網格，例如 Nakah回，bi 和 Deiwert[51] 用其處理網 格自適應問題.在國內，楊國偉等 [52] 利用該方法 進行了飛機氣動彈性的數值模擬.

　　(2) 基于Delaunay 背景網格的網格變形方法 最近 Liu 等 [53] 發展了基于 Delaunay 背景網

　　格的插值方法，較好地解決了非結構網格的動態 變形問題它首先按照 Delaunay 準則生成一套非 常稀疏的背景網格;然后建立計算網格節點與 D令 launay 背景網格單元的一對一映射關系這種映射 關系-旦建立，在運動變形過程中便不會更改.這 樣在 Delaunay 背景網格變形運動以后，計算網格 的節點坐標就可以依照先前建立的一對一映射關 系很快回插計算出新的坐標.與目前普遍使用的

　　彈簧拉伸法相比，該方法的最大優勢在于不用法 代計算，效率得到較大提高，而且對于較大尺度變 形問題，該方法能得到較好的動態網格.

　　Delaunay 背景網格的插值方法可以分為 4 個 步驟: (1) 背景網格的生成; (2) 映射，即建立計算

　　網格節點與背景網格單元的映射關系; (3) 背景網

　　格的運動和變形; (4) 回插，即依照新的背景網格

　　和在第 2 步中得到的映射關系做逆運算得到節點 新坐標. Bowyer[叫和 Watson(55] 指出:對于空間 位置確定的點集，存在一個唯一的三角形 (2D) 或 四面體 (3D) 構造方案. Delaunay 網格生成方法的最大缺陷是不能保證物形邊界的完整性，因此往

　　往需要通過各種方法刪除物體內部的網格，并通 過邊界修正保證邊界的完整性.然而，在這里卻 不必嚴格保持物形邊界，因為De launay 背景網格 只是用于記錄計算網格的相對位置.它只需保證 能夠進行正確的回插，其自身的邊界完整性與計 算網格沒有任何關系圖 9 顯示了通過 Delaunay 背景網格插值方法得到的魚體巡游的動態混合 網格.

　　〈二二::::::=:- C二二、 盧

　　(b)

　　(c)

　　困 9 通過 Delaunay 背景網格插值法得到的混合動網格 (53]

　　4.1.4 變形/重構混合網格生成技術 無論是彈簧拉伸法，還是基于 Delaunay 背景

　　網格的插值法，它們在處理更大變形的非定常問題

　　時(如長距離多體分離)，均存在較大的困難，因為

　　大變形將導致網格單元的質量急劇下降，有時甚至 導致網格的相交，進而影響非定常計算的精度，甚 至導致計算失敗張來平等問、常興華等悶， 50] 在以往發展的定態混合網格技術的基礎上，提出了

　　第 4 期

　　張宋平等:動網格生成技術及非定常計算方法進展綜述

　　435

　　解決此類問題的新策略，即將彈簧拉伸法與局部網 格重構結合起來生成動態網格.首先采用彈簧拉伸 法移動網格節點，然后進行網格質量檢測，如果網 格質量滿足要求，則繼續利用彈簧拉伸法進行下一 步的網格生成;如果變形后的網格不能通過質量檢 測，則在局部進行重構.這便是所謂的變形/重構混 合網格生成技術.與此類似郭正和劉君 [561 等利 用該技術進行了多體分離的動態非結構網格技術 的研究以變形體的動態混合網格生成為例，這里 簡要地給出了動態混合網格生成過程.具體的動 態混合網格的流程圖如圖 10 所示.

　　要對一個運動的物體生成動態網格，首先要 生成初始時刻的靜態網格，然后當物體運動時對 網格進行相應的調整.為了適應復雜外形，最好

　　物體變形

　　的方法當然是混合網格 [1，21 這里，采用了如下 的初始定態混合網格生成策略(以二維問題為例) : 混合網格由貼體的四邊形網格、外圍的矩形網格 和中間連接的三角形網格組成.貼體的大伸展比 的四邊形網格由層推進法生成，用來模擬邊界層. 然后在整個計算域內用四分樹方法生成覆蓋全場 的矩形網格，隨后刪除貼體結構網格附近和內部 的矩形網格，這樣在矩形網格和四邊形網格之間 形成了一個縫隙，之后用陣面推進法生成三角形 網格填補這個縫隙，并用彈簧拉伸法和"對角線交 換"方法對三角形網格進行優化.具體方法可以參 見文獻 [57"'59]. 為了適時控制初始陣面劃分以及 推進過程中新的點的插入，以得到光滑并且適用 于物理問題的網格分布，可以采用控制計算網格 空間步長分布的矩形背景網格技術 [57，60J

　　貼體四邊形間格變形 控制網格分旬的背景間格調整 三角形網格節點松馳

　　笛卡爾網格重新生成 局部三角形網格重構

　　下一時間步

　　圖 10 動態混合網格生成示意圖(彈簧拉伸法和局部網格重構法鍋合)

　　在靜態混合網格生成之后，則按如圖 10 所示 的方法生成動態混合網格.首先，物體周圍貼體的 四邊形網格隨物體變形而變形.具體過程如下:

　　(1) 在初始網格生成中記錄每一層的節點數 目、四邊形數日，并記錄生長方向的點與點之間的

　　聯系，以及點與點之間的線段長度.

　　(2) 當外形變化時，第 1 層的網格點坐標隨著 外形的變形而更新

　　(3) 求第 1層節點的新推進方向并進行法線光 滑.

　　436 力 學 進 展

　　2010 年第 40 卷

　　(4) 按照保存的初始網格信息(推進距離，節

　　點關聯信息等)，利用新的推進方向，更新第 2 層點 的新坐標.

　　(5) 以此類推，直到最后一層四邊形網格點的 坐標更新完畢.推進的步長和初始網格生成中步 長一致，無需重新求解.

　　當貼體四邊形網格更新完成后，外部的矩形 網格保持不變，中間的三角形網格利用彈簧拉伸

　　法進行節點松弛當三角形的質量過差(外接圓和 三角形面積比大于一定的值)或者四邊形網格和 矩形網格相交時，局部矩形網格和三角形網格都 重新生成.作為動態混合網格的生成實例，這里給 出了兩個動態混合網格生成實例，其中圖 11 為一

　　個周期內雙魚對擺的動態混合網格，圖 12 是戰斗

　　由于Delaunay 背景網格插值法在中小變形情 況下能快速生成高質量的動態網格，為此進一步 發展了基于Delaunay 背景網格插值法和局部網格 重構法相結合的動態混合網格生成方法.其網格 生成流程如圖 13 所示.其基本思想與前述稿合 方法相似，只是在變形網格生成過程中采用了 De­ launay 背景網格插值法，而 Delaunay 背景網格本 身的變形則通過彈簧拉伸法實施變形，如果彈簧 拉伸法無法得到合適的 Delaunay 背景網格，背景

　　網格本身也可以重構.作為實例，本文生成了三 段翼型的副翼折轉運動的動態混合網格(圖 14). 后緣副翼從初始狀態折起，變形位移很大，但是利 用 Delaunay 網格插值方法仍能在不改變網格拓撲 結構的情況下，生成全過程的動態混合網格(如圈 14(b)rv 圈 14(d) 所示).這里 Delaunay 背景網格僅 覆蓋副翼附近的區域，其外邊界點當副翼折轉時 適當旋轉.從圖 14 中可以看到，在整個運動過程 中，外場的網格保持靜止，局部網格變形，變形網格

　　第 4 期 張米平等:動網格生成技術及非定常計算方法進展綜述 437

　　的質量保持良好.對于三維問題，國 15 顯示了魚體

　　模型巡游時的動態混合網格，國 16 顯示了變形飛 機模型的動態非結構網格(表面網格和 空間網格)，

　　可以看到，由于棍合方法綜合了 Delaunay 背景網 格插值法和局部網格重構法的優勢，岡此能夠生 成復雜外形的高質量動態網格.

　　Delaunay背景網一幅至南

　　建立計算網格'與 Delaunay

　　背景網格的描值對應關系

　　國 13 動態混合網格生成示忘國 (Delaunay 背景網格插值法和局部網格重構法稍合)

　　國 14 多段典型的品IJ 翼折轉過程的動態混合同格

　　l:l '" i11 It( :20 10 r : 10 古

　　劇 15 三維魚體巡游的動態混合網格

　　YLCz Yt之z

　　國 16 變形飛機的 三維動態非結構網格

　　4.2 基于非結構/混合網格的非定常計算格式

　　4.2.1 流動控制方程 基于非結構/混合網格的流動控制方程實際上

　　與基于結構網格的流動控制方程一致.但是，在非 結構/混合網格情況下，基于貼體坐標變換的有限

　　差分法已經不再適用，這時一般采用積分形式的控

　　制方和.在運動網格情況下，一般可以將非定常 NS

　　方程寫成如下的積分形式(以 二維問題為例)

　　!JQdV+ f陽x +Frny - Q的叫Jds-

　　V 8

　　言1φ(Evnx + FVny)ds = 0 (42)

　　.& ..OO J

　　其中 Vg 為網格面的運動速度， n 為法線方向，V 和

　　8 分別為網格單元的體積分和面積分，其他物理量

　　的定義與 3.2 節相同.

　　4.2.2 時間離散方法 與基于結構網格的非定常計算方法類似，利

　　用非結構/混合網格求解非定常問題，最直接的時

　　間離散方法是多步 Runge-Kutta 顯式推進方法，但 是這種方法要求選取非常小的真實時間步長，否 則無法滿足時間精度的要求，而且時間步長的選 取需要滿足穩定性條件，對計算帶來很大不便.隱 式方法又大致分為兩類:一類是直接的隱式離散 方法，另一類是雙時間步方法 [61，叫.在實際工程

　　第 4 期 張來平等:動網格生成技術及非定常計算方法進展綜述 439

　　問題的非定常計算中，一般都選用雙時間步方法 . 在虛擬時間步遠代過程中，則可以采用成熟的加 速收斂技術，如隱式計算方法、多重網格方法、并 行計算方法等.這與基于結構網格的非定常計算 方法類似，只不過需要采用針對非結構網格的隱 式計算方法、多重網格方法和并行計算方法等.關 于子送代的加速收斂技術將在 4.2.5 節中介紹.

　　4.2.3 不可壓縮的非定常計算方法 在不可壓縮流中，非定?，F象也非常普遍.

　　對于不可壓縮問題，常用的計算方法有預處理方 法 [62 64] 和虛擬壓縮方法 [65 6η. 預處理方法是

　　在控制方程時間項中引入了預處理矩陣 (r)

　　r!JQ機 f [E[nx + F[ny - Q(Vg• n)1

　　V 8

　　R cof 肌x +FVny)…側

　　8

　　通過選取合適的矩陣，解決低速流計算的收

　　斂問題，具體參見文獻歸2rv64l. 虛擬壓縮的概念 則是 Chorin[65] 首先提出的，他通過在連續性方程中引入壓力對虛擬時間的偏導數從而把壓力場和 速度場藕合起來.

　　預處理方法和虛擬壓縮方法最先應用于定常 問題的計算，以加速不可壓縮流計算的收斂歷程， 后來一些學者結合雙時間步方法，將預處理方法 和虛擬壓縮方法應用于非定常計算 [62，66，67J 當然， 關于不可壓縮流的計算還包括 SIMPLE 算法[倒l、 微可壓模型 (SCM) 計算方法陽，而]和當前比較熱 門的技影法 [71，72] 等等，這些方法也相繼被推廣應 用于不可壓縮流的非定常計算.由于篇幅的限制， 這里不再詳述.

　　4.2.4 動邊界問題 非定常計算往往要碰到動邊界問題，解決這

　　一問題的方法有很多，例如非慣性系 F的剛性 運動網格方法、內置邊界方法 [73 75] (immersed boundary，IB) 、動網格計算方法等等.關于非慣 性系 F的非定常計算方法，在前節已經介紹，這里 不再重復

　　(1) IB 方法

　　對比邊界保形 (bound町 conforming) 的動網 格計算方法，田方法實際上是一種非邊界保形法

　　(bound町non-conforming method). IB 方法采用的 網格一般為笛卡兒網格或者簡單的圓柱網格，網 格不需要與物體邊界保持一致，因而簡化了網格 生成的難度.采用 IB 方法，需要在控制方程中引 入源項以保證物體邊界的無滑移條件;對于包含 物理邊界的網格單元還需要在控制方程中引入 質量源以保證質量守恒，以不可壓縮流為例，其微 分形式的控制方和變化為

　　(等)儼+V" (UrU儼) = -V'P+ 合2Ur +

　　(44)

　　j' V" Ur-q=O

　　其中 f 為動量力，q 為質量源. 這種方法不用考慮復雜的物體外形，可以采

　　用簡單的網格，因此很受重視，近來一些作者相繼

　　將其推廣應用于非定常問題的計算，如 Kim 等 [75] 結合坐標變換，將該方法應用于非慣性系下的非 定常計算;王亮 [76] 通過內置邊界和自適應笛卡兒 網格的方法，對三維魚體的自主推進進行了數值 模擬;邵雪明等 [77] 利用該方法模擬了二維多魚巡 游的問題等.盡管這種方面在一些簡單外形的數 值模擬中取得了較大的成功，但是對于復雜問題， 如何適當選取邊界處的分布函數，嚴格保證邊界 條件，依然是值得深入研究的問題.因此，邊界保 形的動網格方法仍然是主流方法.

　　(2) 兒何守恒律 與結構動網格計算一樣，在非結構動網格計

　　算中，必須考慮兒何守恒率. Lesoinne 和 Farhat[78]

　　闡述了動網格計算的幾何守恒率 (GCL) 問題.他 們認為動網格的計算必須滿足幾何守恒率，并分 析了任意拉格朗日·歐拉 (ALE) 有限體積法、有 限元方法、時空穩定的有限元方法等算法的兒何

　　守恒率問題.其中，時空穩定的有限元算法能自動 滿足幾何守恒率 j ALE 方法則需要保證每個網格 單元的體積變化等于每個面在運動過程中掃過的 體積之和，并提出了用"中間點原則"來求解網格 面運動速度的方法.他們還通過計算分析了兒何 守恒率對計算結果的影響(如國 17 所示)，在此類 流國輯合問題的計算中，如果不滿足幾何守恒率， 會導致計算結果的振蕩.

　　v.

　　θ

　　對于式 (42)，如果流場為均勻流，則要求下式 得以滿足

　　!= )vgndSf (45)

　　這便是所謂的兒何守恒律的數學表達.將式 (45)

　　改寫為差分形式，有

　　A; =L

　　Ern+l τrn

　　勺 勺 與ndSf (46)

　　事實上，每個單元的體積變化量等于每個面在運 動過程中掃過的體積之和，即

　　只叫 -η= 藝A巧 (47)

　　比較式 (46) 和式 (47)，可得

　　gn - .6. f

　　.6.Vf = .6.tvgndSf 、 Vgn = 全旦二 (48)

　　tdS

　　利用上述方法求出 Vgn ，帶入原控制方秤，即可 自動滿足兒何守恒率 [79 ，80)

　　4.2.5 加速收斂技術

　　對于非定常計算，當采用雙時間步方法，要求 每個時間層間的子選代具有較高的收斂效率，這 時可以將定常流計算中很多成熟的加速收斂技術

　　引入進來.用隱式方法求解時，控制方和最后變為

　　一個非線性的方程組.如果用全隱式的牛頓選代 法求解此非線性方程組，則每一時間步都需要對 矩陣求逆，計算量會非常大-般可以將方程組線 性化，采用法代法來求解代數方和組.

　　簡單的法代法如雅克比 (Jacobian) 法代、高 斯 - 賽德爾法代 (GS)[81 84] 已經在隱式方程的求 解中得到了成功的應用.但是，隨著網格數量的

　　增加，這種選代方法的收斂效率會降低.日前，應 用較多的是 LU-SGS (lower upper symmetric Gauss­ Seidel) 方法. LU-SGS 方法最初由 Jameson 和 Yoon(85) 應用于結構網格的計算中，通過雅克比矩 陣的岡式分解結合高斯-賽德爾地代，提高了收 斂效率，之后這種方法被成功地應用于非結構網 格 [86] 預處理的 GMRES (general minimum resid­ ual) 方法同樣可以提高收斂效率，尤其當結合 ILU 因式分解法時，收斂效率會更高 [87，88). Luo 等 [89，90] 采用 LU-SGS 方法作為 GMRES 方法的預 處理算子，發展了 GMRES-LU-SGS 方法，得到了 良好的收斂效率此外， Chen 和 Wang[91] 發展了

　　-種基于混合網格的塊 LU-SGS 方法，在內存量 增加不多的情況扎得到了與全隱式方法相當的 收斂效率，在實際應用中得到了驗證;后來，張來

　　平等 [80，92} 將該方法推廣應用到動網格非定常計

　　算，得到了較好的非定常計算效率.

　　在不可壓和亞跨聲速流的計算中，廣泛采用 多重網格法.多重網格法最早由 Brandt!93} 引入到 求解非線性偏微分方程的離散形式，后來這種方 法在 CFD 計算中的應用越來越廣泛.它通過在兒 套不同疏密程度的網格上循環地對代數方程進行 求解，使得不同頻率的誤差分量均勻地得到衰減， 從而極大地提高了收斂速率

　　對于非結構網格而言，生成兒套疏密程度不 同的網格，要比結構網格困難的多. 方法之一是 生成相互嵌套的網格，即疏網格不斷通過密網格 單元的聚合而得到 [94，95) 這種網格生成方法效率 高，但是當聚合的次數過多時，會使得網格質量急 劇下降(尤其是在邊界 層 內)， 有可 能導致計算的

　　第 4 期 張米平等:動網格生成技術及非定常計算方法進展綜述 441

　　失敗;可能的解決方案是在聚合過程中引入一些 限制條件，保證聚合后的網格質量.張來平等 [96J 在基于動態混合網格的不可壓縮流非定常計算中， 應用了基于混合網格的多重網格技術，進一步提

　　闊的應用前景.

　　2.5

　　2.0 • 實驗結果

　　_ ._._.-無貓流

　　高了子選代的效率.

　　另一種生成多重網格的方法是生成相互獨立 的網格，網格之間通過線性插值進行流場信息的 傳遞 [97，98J 這種方法能夠保證疏、密網格的質量 ， 可以生成多套網格，提高計算效率和精度，缺點是 網格生成難以自動化，需要重復地在同-15<:域內 生成網格.

　　5 結束語

　　本文從非定常計算采用的網格體系的不同 ， 分兩個部分對非定常計算方法進行了綜述，即基 于結構網格的非定常計算方法和基 于非結構/棍 合網格的非定常計算方法.其間，重點介紹了應 用于非定常計算的動態網格生成技術、時間離散 方法、兒何守恒率、加速收斂技術等.在介紹國內 外進展的過程中，也介紹了作者在非定常計算方 法方面的研究成果.

　　從正文的論述中，可以看到，各種非定常計算 方法均有各自的優缺點.對于基于結構動網格的 非定常計算方法，其網格拓撲簡單、計算效率高， 所以在簡單外形的非定常流動機理研究中經常采 用;對于復雜的工程實際問題，顯然基 于動態混合 網格的非定常計算方法具有更好的優勢.

　　就國內外日前的發展而言，非定常計算方法 己經取得了長足的進展，在實際工程應用中也發 揮了較好的作用.但是，應該清醒地 看到，仍有許 多工作需要深入，尤其是非定常計算的時間精度

　　問題、復雜構型的動態混合網格自動生成的 魯棒

　　性和實用化問題、計算效率的提高問題等等.當 前，大規模并行計算已經成為一種趨勢，對于計算 量巨大的非定常問題，采用并行計算技術是必然 的選擇.

　　最后需耍說明的是，在非定常計算方法的研 究與應用方面，國內外都有很多卓有成效的工作， 文中引述的資料肯定不夠全面;而且，對某些問題 的看法可能帶有個人的偏好岡此，希望讀者給予 批評指正

　　參考文獻

　　1 Thompson JF，Soni BK ，Weatherill NP. Handbook of Grid

　　Generation. USA: CRC Pr四民 1998

　　2 ∞

　　Baker TJ. Mesh generation: 町t or science? P.何可9. Aero.

　　Sci.，2 5，41: 29"'63

　　3 Steinthorson E，Liou MS，Povinelli LA. Development of 創1

　　explicit multiblockjmultigrid flow solver for viscous flows in complex geometries. A1AA 93盧2380，1993

　　4 Sheng C，Taylor L，Whitfield D. Multiblock multigrid sφ lution of three dimensional in∞mpr四sible turbulent flows

　　about appended subma.rine confìguration. AIAA 95-0203，

　　1995

　　5 Flor臼 J，Reznick SG，Holst TL，et.al. τ'ransonic Navier­

　　∞

　　Sto倒四lution for a figh飼r-like confìguration . AIAA 87- 32，1987

　　推薦期刊：《計算機應用與軟件》創刊于1984年，由上海市計算技術研究所和上海計算機軟件技術開發中心共同主辦，注重刊登反映計算機應用和軟件技術開發應用方面的新理論、新方法、新技術以及創新應用的文章，致力于創辦以創新、準確、實用為特色，突出綜述性、科學性、實用性，及時報道國內外計算機技術在科研、教學、應用方面的研究成果和發展動態的綜合性技術期刊，為國內計算機同行提供學術交流的平臺。