廖佳文， 張力天， 周 璇， 李水鄉(xiāng)*

(1.北京大學 工學院，北京 100871; 2.北京應用物理與計算數(shù)學研究所,北京 100094)

1 引 言

在流固耦合和多介質(zhì)流體[1]等實際工程問題中經(jīng)常需要處理邊界移動問題，此時需要根據(jù)邊界的運動或變形來改變計算網(wǎng)格[2]。網(wǎng)格重構是解決此類問題的一種思路，但計算代價較大[3]。另一種經(jīng)常采用的處理方法是動網(wǎng)格法[4]。在動網(wǎng)格方法中，計算網(wǎng)格隨著邊界移動發(fā)生變形以適應新的邊界[5]。常見的動網(wǎng)格方法[6]有彈簧近似法[7]、徑向基函數(shù)插值法[8]、距離函數(shù)插值法、偏微分方程法和彈性體法等。

從顆粒填充領域的松弛算法受到啟發(fā)，周璇等[9]提出了一種新穎的基于球松弛算法的動網(wǎng)格方法，稱之為動網(wǎng)格松弛法。與彈簧法等虛擬結構類網(wǎng)格變形方法相比，動網(wǎng)格松弛法的網(wǎng)格變形效率更高，網(wǎng)格質(zhì)量更好。然而，動網(wǎng)格松弛法與徑向基函數(shù)插值法等代數(shù)插值類方法相比，在大規(guī)模網(wǎng)格變形中計算效率偏低。

已有的研究表明，采用多套不同節(jié)點數(shù)網(wǎng)格的多重網(wǎng)格策略能夠有效地提高動網(wǎng)格方法的求解效率，是動網(wǎng)格方法改進的一種重要手段。胡會朋等[10]提出了基于背景網(wǎng)格簇和Delaunay圖映射的動網(wǎng)格生成方法，該方法變形能力高于彈簧法，且能夠直接推廣到三維。孫巖等[11]提出了基于徑向基函數(shù)與混合背景網(wǎng)格的動網(wǎng)格方法，可有效減少徑向基函數(shù)基點數(shù)，有助于提高徑向基函數(shù)插值的效率。

本文在原有動網(wǎng)格松弛法的基礎上，提出了一種基于二重網(wǎng)格的動網(wǎng)格松弛法改進方法。

2 基于二重網(wǎng)格的動網(wǎng)格松弛法改進

本文引入二重網(wǎng)格策略改進了動網(wǎng)格松弛法。通過二重網(wǎng)格映射，將邊界位移傳遞到整個網(wǎng)格計算域，達到減少動網(wǎng)格松弛法迭代次數(shù)的目的，從而有效地提高了計算效率。改進后的動網(wǎng)格松弛法包含二重網(wǎng)格映射、球松弛和后優(yōu)化三個步驟。

2.1 二重網(wǎng)格映射

在原有的非結構計算網(wǎng)格之外，重新在計算區(qū)域內(nèi)生成一套更加稀疏的非結構網(wǎng)格，即粗網(wǎng)格。原有的計算網(wǎng)格稱之為細網(wǎng)格。粗細兩套網(wǎng)格形成二重網(wǎng)格。

利用面積法，建立細網(wǎng)格節(jié)點與粗網(wǎng)格三角形單元之間的對應關系[12]。

αi=si/s

(i=1,2,3)(1)

式中αi為細網(wǎng)格節(jié)點在所屬粗網(wǎng)格單元中的面積坐標，si(i=1,2,3)分別為細網(wǎng)格節(jié)點與所屬粗網(wǎng)格單元的節(jié)點連接成的三個三角形的面積，s為該粗網(wǎng)格三角形單元的面積。利用松弛算法進行粗網(wǎng)格變形后，通過映射得到細網(wǎng)格節(jié)點的新位置(x,y)為

(i=1,2,3)(2)

式中xi和yi(i=1,2,3)為變形后粗網(wǎng)格單元節(jié)點的坐標。

2.2 球松弛算法

在球松弛算法中，利用節(jié)點分布構造球分布進行松弛移動，得到新的節(jié)點分布，其主要步驟如下。

(1) 構造球分布。將節(jié)點位置作為球心，連接該節(jié)點的邊長平均值的一半作為球半徑(目的是最小化球之間的重疊量并盡可能相切)，如式(3)。采用此方法對計算域進行布球，將計算網(wǎng)格的節(jié)點分布轉(zhuǎn)化為球分布。

(3)

式中rj為節(jié)點j處的球半徑，m為與節(jié)點j相連的節(jié)點個數(shù)，li j為節(jié)點j與節(jié)點i相連的邊長。

(2) 用球松弛算法調(diào)整球的位置，以得到區(qū)域內(nèi)球的均勻分布。如圖1所示，球松弛將相交和相離的球盡可能轉(zhuǎn)化為相切，調(diào)整后的球心位置計算公式為

Rj i=Rj+(Ri-Rj)[(ri+rj)/di j]

(4)

式中Rj i為球i與附近的球j進行松弛操作得到的新的球心位置，ri和rj分別為球i和球j的半徑，di j為球i與球j之間的球心距離。通過建立背景網(wǎng)格[13]，將鄰接球搜索局部化以減少鄰接球搜索時間。球i移動后的球心位置計算式為

(5)

圖1 球松弛[9]

在一次迭代中調(diào)整每個球的位置，經(jīng)過多次迭代直到所有球的位移量足夠小，則迭代停止，并得到新的節(jié)點分布。再利用初始網(wǎng)格的連接關系將這些節(jié)點連接起來，得到變形后的網(wǎng)格。

2.3 后優(yōu)化

為了防止非法網(wǎng)格的出現(xiàn)，本文采用文獻[9]的后優(yōu)化過程。圖2給出了后優(yōu)化過程的二維示例，對于節(jié)點j，與其連接的節(jié)點形成的網(wǎng)格區(qū)域定義為區(qū)域D,對于區(qū)域D的每條邊，向區(qū)域內(nèi)部嘗試做等邊三角形。等邊三角形的另一個頂點位置矢量為Ri j(i=1,2,…,n),n為區(qū)域D的邊長總數(shù)，節(jié)點j的新位置Rj由所有的相對位置Ri j決定，

(6)

式中hi j為節(jié)點j到區(qū)域D邊i的垂直距離。后優(yōu)化過程可以防止非法網(wǎng)格的出現(xiàn)。

圖2 后優(yōu)化算法

3 數(shù)值算例與討論

本文通過二維NACA0012翼型轉(zhuǎn)動、矩形平動和轉(zhuǎn)動以及三維梁彎曲三個典型算例驗證改進松弛法在動網(wǎng)格計算中的特性和效率。計算平臺為3 GHz主頻，16 GB內(nèi)存的臺式機。

3.1 NACA0012翼型轉(zhuǎn)動

平面NACA0012翼型在矩形區(qū)域內(nèi)逆時針轉(zhuǎn)動[9]，計算區(qū)域的大小為20C×20C，其中C為翼型的弦長。三角形計算網(wǎng)格節(jié)點數(shù)為29382。將翼型的運動過程分為若干步，每步翼型轉(zhuǎn)動角度為5°。在細網(wǎng)格節(jié)點數(shù)固定的情況下，選取多組不同疏密程度的粗網(wǎng)格進行對比。粗細網(wǎng)格節(jié)點數(shù)之比(以下稱網(wǎng)格比，在0～1之間)分別取為0.36,0.42,0.46,0.51和0.64。

圖3(a)為改進松弛法在翼型旋轉(zhuǎn)110°后的邊界局部網(wǎng)格，可以看出，網(wǎng)格變形后的物面邊界附近保持了良好的網(wǎng)格質(zhì)量[14]。圖3(b，c)分別為網(wǎng)格變形過程中的最差網(wǎng)格質(zhì)量和平均網(wǎng)格質(zhì)量對比?？梢钥闯?，改進松弛法的變形能力(極限變形量)與原動網(wǎng)格松弛法相同，但網(wǎng)格質(zhì)量稍有下降，不同網(wǎng)格比下網(wǎng)格質(zhì)量差別不大。圖3(d)為不同網(wǎng)格比下改進松弛法與原動網(wǎng)格松弛法的計算時間之比。綜合在不同網(wǎng)格比下變形后的網(wǎng)格質(zhì)量和計算時間數(shù)據(jù)可以發(fā)現(xiàn)，隨著網(wǎng)格比的增加，變形網(wǎng)格質(zhì)量是逐步提高的，但計算時間先減小后增大。NACA0012算例顯示，當網(wǎng)格比為0.51時，改進松弛法時間效率最優(yōu)，計算時間減小到原動網(wǎng)格松弛法的50%，同時變形能力保持不變，變形網(wǎng)格質(zhì)量也基本保持不變。

3.2 矩形平動和轉(zhuǎn)動算例

平面正方形區(qū)域邊長為500個長度單位，初始時矩形水平置于計算區(qū)域中心，矩形的長和寬分別為25和10個單位長度。矩形運動為逆時針旋轉(zhuǎn)同時向正方形區(qū)域的左下角平動[15]。三角形計算網(wǎng)格節(jié)點數(shù)為9012。矩形運動和網(wǎng)格變形分步進行，矩形每步旋轉(zhuǎn)30°并分別向左方和下方各移動 3個單位長度。網(wǎng)格比分別取為0.16，0.25，0.33，0.41，0.45，0.49，0.53，0.65和0.79。

圖4(a)為變形21步后的計算網(wǎng)格，圖4(b)為改進松弛法與原動網(wǎng)格松弛法的計算時間之比與網(wǎng)格比的關系。從圖4(b)可以看出，所有網(wǎng)格比下的計算時間比都小于1，表明基于二重網(wǎng)格的改進松弛法能提高動網(wǎng)格松弛法的計算效率。對不同網(wǎng)格比進行比較，計算時間隨著網(wǎng)格比的增大先減小再增大。最優(yōu)的時間比出現(xiàn)在網(wǎng)格比為0.49處，所用時間為原動網(wǎng)格松弛法計算時間的47.4%。

3.3 三維梁彎曲算例

三維長方體梁尺寸為20×2×0.5單位長度,計算區(qū)域是半徑為25個單位長度的球[13]。計算網(wǎng)格中有13712個節(jié)點，63385個四面體單元。初始時長方形梁水平置于球區(qū)域的中心，如圖5(a)所示。變形過程中梁向上彎曲，如圖5(b)所示。變形過程中梁的中心線圓弧對應的圓心角從0°變化至180°。變形過程分為5步，每步圓心角的變化為36°。

圖5(a)為梁附近的局部初始網(wǎng)格，圖5(b)為梁彎曲變形第5步后的局部網(wǎng)格(圓心角180°)，圖5(c)為改進松弛法和原動網(wǎng)格松弛法在不同網(wǎng)格比下的最差網(wǎng)格質(zhì)量對比，圖5(d)給出了改進松弛法與原動網(wǎng)格松弛法的計算時間之比與網(wǎng)格比的關系?？梢钥闯?，改進松弛法不但與原動網(wǎng)格松弛法變形能力相同，而且能夠有效地提高計算效率。網(wǎng)格比在0.52時計算效率達到最優(yōu)，這與前面的二維算例所得結論一致。但與二維算例相比，其最優(yōu)時間比僅為78%。這是因為該三維算例中球松弛迭代的次數(shù)遠小于二維算例，故二重網(wǎng)格縮減迭代次數(shù)的作用減弱，并且三維映射比二維所需計算時間更多。

4 結 論

本文提出了基于二重網(wǎng)格策略的改進動網(wǎng)格松弛法，通過二重網(wǎng)格映射來減少松弛法的迭代次數(shù)，從而提高動網(wǎng)格松弛法的計算效率。二維和三維算例的結果表明，(1)改進后的方法在保持網(wǎng)格變形能力和變形網(wǎng)格質(zhì)量基本不變的前提下，能夠有效地提高網(wǎng)格變形的計算效率； (2) 隨著粗細網(wǎng)格比的增加，平均網(wǎng)格質(zhì)量和最差網(wǎng)格質(zhì)量均單調(diào)上升，但計算時間先增大后減小。其原因在于隨著網(wǎng)格比的增加，雖然粗網(wǎng)格變形所需的計算時間不斷增加，但松弛算法的迭代次數(shù)會相應減少，因此疊加起來總的計算時間出現(xiàn)先增大后減小； (3) 粗細網(wǎng)格比在0.5左右時達到最優(yōu)計算效率。

本文研究表明，采用二重網(wǎng)格策略并選取適當?shù)拇旨毦W(wǎng)格比可以有效地提高動網(wǎng)格松弛算法的計算效率，為該方法在大規(guī)模網(wǎng)格計算中的應用提供了有力的支撐。