混合非結(jié)構(gòu)網(wǎng)格下壁面最短距離的快速計算方法

王剛, 曾錚, 葉正寅

(西北工業(yè)大學(xué) 航空學(xué)院流體力學(xué)系, 陜西 西安 710072)

在高雷諾數(shù)流動問題的數(shù)值模擬中，通常需要在雷諾平均N-S方程中引入湍流模型來刻畫湍流的脈動效應(yīng)。目前常用的湍流模型包括一方程S-A模型及二方程(k-ε模型，k-ω模型等)及多方程模型(雷諾應(yīng)力模型等)[1]。上述的湍流模型均需要利用壁面距離作為背景參數(shù)。對于結(jié)構(gòu)網(wǎng)格而言，由于網(wǎng)格拓撲結(jié)構(gòu)的規(guī)律性，其壁面距離的計算效率比較容易保證，而非結(jié)構(gòu)網(wǎng)格舍去了網(wǎng)格連接的結(jié)構(gòu)性和正交性限制，其拓撲結(jié)構(gòu)具有很大的隨意性，導(dǎo)致壁面距離參數(shù)的確定比較復(fù)雜。目前非結(jié)構(gòu)網(wǎng)格中普遍采用的壁面距離計算方法是枚舉法，即依次計算比較空間每個單元到所有表面單元的距離，這種算法精度足夠但效率很低。因此，有必要研究適合非結(jié)構(gòu)網(wǎng)格的高效、快速的壁面距離搜索算法。國內(nèi)外對點到曲面最短距離的求法有很多文獻資料，但大部分均是在曲面方程已知的情況下求解的[2-4]。如果采用這些方法，就需要重新分塊擬合曲面方程，這需要耗費額外的工作量，不適用于計算流體力學(xué)中的離散曲面問題。而N.E.Renato等人提出的快速推進算法[5]，是針對非結(jié)構(gòu)網(wǎng)格下的基于有限元方法求解距離函數(shù)的優(yōu)化算法，其提升效率不錯，但這種方法需要控制的變量太多，且其邊界處理并不理想。T.Maruto等人采用的面集法[6]也是利用局部化的思想，將計算域限制在精確值附近的一個小區(qū)域內(nèi)，提高計算的效率。但是這種算法需要單獨考慮復(fù)雜邊界的處理，在邊界處常常會出現(xiàn)意想不到的問題。為了提高非結(jié)構(gòu)網(wǎng)格中壁面距離搜索算法的效率，本文提出一種壁面逐層推進法。該算法只需計算流場中任意一個網(wǎng)格單元到目標(biāo)面元(最短距離面元)及與其相鄰的幾個面元的距離即可，不需像枚舉法一樣計算到所有面元的距離，其算法復(fù)雜度比前述各類方法要低很多，且精度完全滿足流場求解的要求，是一種理想的局部尋優(yōu)算法。

1 壁面逐層推進法

1.1 算法原理

對于一個給定的已經(jīng)離散化的表面，令S={S1,S2,…,Sm}代表該表面集合，其中Si代表壁面上的第i個面元。流場中任意一點用p=(xp,yp,zp)表示令qi(xq,yq,zq)為Si上距離p最近的一點，則最短距離函數(shù)為

(1)

當(dāng)m→∞時，d可視為連續(xù)。在實際應(yīng)用中，取qi為Si面心，p為流場某一網(wǎng)格單元的中心。當(dāng)m充分多的情況下，距離函數(shù)d可近似為連續(xù)。

圖1 逐層推進法原理示意圖

將上述連續(xù)性假設(shè)應(yīng)用到圖1a)所示情形，第一層單元是表面網(wǎng)格單元層，Si(i=1,2,3,4,5)是表面單元編號。我們?nèi)網(wǎng)格來觀察，根據(jù)距離函數(shù)的連續(xù)性，與B網(wǎng)格最近的固壁面元信息應(yīng)當(dāng)來自于A、C、D3個網(wǎng)格單元。假設(shè)A、C、D各自存儲的最近面元編號分別為S3、S2、S4，則B的最近面元應(yīng)是這3個面元中到B距離最近的一個。

逐層推進法的另一個主要思想是將網(wǎng)格分層，從壁面單元層開始，向外逐層推進搜索計算。在沒有明顯分層結(jié)構(gòu)的非結(jié)構(gòu)網(wǎng)格中，首先設(shè)置2個數(shù)組RS和RS1，前者用來存儲當(dāng)前層的網(wǎng)格信息，后者存儲下一層的網(wǎng)格信息，網(wǎng)格信息中增加一個“染色”標(biāo)記FLAG用來標(biāo)記該單元是否已經(jīng)存儲了最短距離，并以此來尋找下一層的網(wǎng)格。只要RS中的單元全部被染色，則所有RS中單元未被染色的鄰居就視作空間的下一層。如圖1b)，第1層全被染色后(灰色單元)，找到所有它們的未染色鄰居(斜線單元)，視作第2層。該層染色后，標(biāo)有數(shù)字的網(wǎng)格單元視作第3層，依此類推。 當(dāng)整個空間網(wǎng)格中每一個單元都被染色，程序結(jié)束。

逐層推進法將連續(xù)性假設(shè)與分層思想結(jié)合起來，在圖1a)中，第1層為表面網(wǎng)格單元層，很容易確定它們的格心分別到自身所含的壁面單元最近。對于第2層的B網(wǎng)格，它作為第1層A網(wǎng)格的未染色鄰居被找到，計算它到S3的距離并將該距離及S3存入B的信息體，再計算它的所有鄰居到S3的距離，分別存入各自對應(yīng)的信息體中。當(dāng)C網(wǎng)格作為第1層E網(wǎng)格的未染色鄰居被找到時，同樣進行如上操作，B網(wǎng)格作為它的鄰居計算了到S2的距離，如果該距離小于之前存儲的距離，則更新B的信息體。當(dāng)一層全部計算完畢后將它們?nèi)旧?，繼續(xù)下一層的計算，直到空間所有網(wǎng)格均被染色。

1.2 搜索算法的進一步改進

在實際應(yīng)用中如果遇到很復(fù)雜的拓撲結(jié)構(gòu)，前述算法可能會出現(xiàn)由于網(wǎng)格疏密不同而造成的反向推進?？疾烊鐖D2所示的增升裝置構(gòu)型的計算網(wǎng)格進行距離計算結(jié)果，很容易看出圖2a)的計算結(jié)果有誤。出現(xiàn)該問題的原因是網(wǎng)格的疏密不同造成了逐層推進時信息傳遞出現(xiàn)了錯誤。被深色“吃掉”的區(qū)域由于該處表面網(wǎng)格密集導(dǎo)致推進速度遠慢于兩側(cè)，當(dāng)兩側(cè)推進到較遠時，本屬于密集區(qū)域的網(wǎng)格單元會作為稀疏區(qū)域單元的未染色鄰居被找到，此時密集區(qū)域內(nèi)層的單元尚未被染色，因而被認作了是外層單元的下一層。最后，密集區(qū)域的正反2種推進方式會在某一處相遇，形成一個交匯邊界，在交匯邊界的內(nèi)側(cè)，其所存距離是正確的，而外側(cè)所存距離則是錯誤的，因為外側(cè)網(wǎng)格的最近表面被認為來自物面的稀疏區(qū)域，如圖3a)所示。

圖2 未改進的算法可能遇到的問題示意圖

圖3 問題原因及解決方式

從壁面1向外推進，淺色單元表示儲存了正確的面元，即壁面1，深色單元儲存的是錯誤的面元，即壁面2。交匯時發(fā)生圖3b)所示情形，算法停止。

本文對前述算法進行了改進。解決方法是在交匯邊界處計算鄰居單元的距離時，若發(fā)現(xiàn)所存距離需要更新但該單元已被染色，此時需要將該單元的FLAG重新置為FALSE，即去除染色，使推進過程可以重新由近場向遠場進行。如圖3c)，深色單元被依次去除染色，程序可以繼續(xù)進行。正反兩種推進方式還可能發(fā)生圖3d)所示的交匯情形，處理方法是一樣的。

2 算例分析

以類空天飛機外形[7]、三維增升裝置[8]及DLRF6翼身組合體[9]的計算網(wǎng)格為考察對象，將逐層推進法的壁面距離計算結(jié)果與枚舉法比較，測試算法的精度，并考察逐層推進法的計算結(jié)果在整個流場空間的平均相對誤差以及最大相對誤差。令

(2)

(3)

平均相對誤差

(4)

最大相對誤差

(5)

利用逐層推進法的壁面距離計算結(jié)果，對上述外形的黏性流場進行RANS方程求解計算，計算采用的求解算法參見文獻[10]。通過流場計算結(jié)果(截面Cp分布)與實驗值的比較，驗證逐層推進法計算的壁面距離在流場計算中的有效性。

2.1 類空天飛機模形

本例的類空天飛機模型采用混合非結(jié)構(gòu)網(wǎng)格，網(wǎng)格單元總數(shù)為1 845 705個，表面網(wǎng)格總數(shù)為61 442個。其流場計算狀態(tài)為馬赫數(shù)4.94，雷諾數(shù)5.26×107，迎角為0°，采用的湍流模型是S-A一方程模型。表1給出了最短物面距離的計算時間比較結(jié)果。

在本例中枚舉法采用的是八核并行的運行機制，而逐層推進法采用單核計算。由表1可知，枚舉法的計算時間為逐層推進法的58倍。如果單純考慮實際計算量，逐層推進法的計算效率約為枚舉法的464倍。

表1 2種算法計算時間比較(空天飛機)

圖4 空天飛機模型網(wǎng)格切片圖 圖5 空天飛機壁面距離等值線圖

圖4給出的是機體某鉛垂截面的網(wǎng)格切片圖，圖5給出該截面計算所得的壁面距離等值線圖。圖中枚舉法的等值線設(shè)為虛線，逐層推進法的等值線設(shè)為實線，實線和虛線吻合，說明距離計算結(jié)果精度較好。新算法在整個流場中的平均相對誤差δave=8.257×10-5，最大相對誤差δmax=0.019，相對誤差在各區(qū)間的分布如表2所示。

表2 類空天飛機模型壁面距離相對誤差分布范圍

雖然表2中有少數(shù)的網(wǎng)格其壁面距離的相對誤差在1%到10%之間，但其在總網(wǎng)格中所占的比例極小，且位置均在附面層外，因此對流場計算結(jié)果不產(chǎn)生實質(zhì)性影響。圖6給出了子午線截面處Cp的計算值與實驗值對比。圖中計算值與實驗值符合良好，進一步證明逐層擴散法計算得到的壁面距離具備有效的精度，可以適用于本算例的高超聲速流場計算。

圖6 子午線截面處Cp計算值與實驗值對比圖

2.2 增升裝置模形

該模型采用結(jié)構(gòu)化網(wǎng)格，共有5 957 632個網(wǎng)格單元，表面有65 248個單元。其流場計算狀態(tài)為馬赫數(shù)0.2，雷諾數(shù)4.3×106，迎角為13°。計算采用S-A一方程湍流模型。2種算法的最短物面距離計算時間如表3所示。

表3 2種算法計算時間比較(增升裝置)

由表3可知，枚舉法的計算時間為逐層推進法的65.7倍，實際計算量枚舉法約為逐層推進法的525.6倍。本算例由于和類空天飛機模型的表面網(wǎng)格數(shù)目大致相當(dāng)，因此效率提高倍數(shù)也基本相同。

圖7 增升裝置模型網(wǎng)格切片圖 圖8 增升裝置壁面距離等值線圖

7給出了增升裝置模型機翼某鉛垂截面的網(wǎng)格切片圖，這是一個結(jié)構(gòu)化網(wǎng)格。圖8是該截面距離等值線示意圖。由圖可見，等值線分布均勻，與枚舉法吻合良好，說明逐層擴散法在此算例中計算精度較好。以枚舉法的結(jié)果作為標(biāo)準(zhǔn)，其中整個流場平均相對誤差δave=8.25×10-7，最大相對誤差δmax=0.119。相對誤差在各區(qū)間的分布列表如下：

表4 增生裝置模型壁面距離相對誤差分布范圍

由表4可得，絕大部分網(wǎng)格的相對誤差分布在一個極小的范圍內(nèi)。將逐層擴散法計算所得的壁面距離引入流場求解中，圖9給出了翼展方向50%處的截面Cp的計算值與實驗值對比，結(jié)果顯示計算值與實驗值吻合良好。該算例的流場結(jié)果表明，逐層擴散法計算的壁面距離同樣可以良好適用于結(jié)構(gòu)化網(wǎng)格。

圖9 增升裝置翼展50%截面Cp計算值與實驗值對比

2.3 DLR-F6翼身組合體模型

該模型來自第二屆AIAA 阻力預(yù)測會議，本文測試的是非結(jié)構(gòu)混合網(wǎng)格，其網(wǎng)格單元總數(shù)為13 641 688個，表面網(wǎng)格單元有250 245個。流場計算狀態(tài)為馬赫數(shù)0.75，雷諾數(shù)3.0×106，迎角為1°，采用MSST兩方程模型。兩種算法的最短物面距離計算時間如表5所示：

表5 2種算法計算時間比較(DLR-F6)

由表5知,對于本例中表面網(wǎng)格數(shù)量達到25萬的大型網(wǎng)格，枚舉法的計算時間為逐層推進法的121倍，實際計算量為逐層推進法的968倍，可見所提算法在壁面網(wǎng)格數(shù)目很大的情況下相比枚舉法更具優(yōu)越性。

圖10 DLR-F6模型網(wǎng)格切片圖 圖11 DLR-F6模型距離等值線圖

圖10和圖11分別給出了通過發(fā)動機軸線的鉛垂截面的網(wǎng)格切片和距離等值線示意圖。逐層擴散法計算的壁面距離等值線在壁面外圍呈均勻分布，與枚舉法吻合良好，證實逐層擴散法計算的壁面距離在此算例中精度較好。以枚舉法的結(jié)果作為標(biāo)準(zhǔn)，其中整個流場平均相對誤差δave=2.129 812×10-4，最大相對誤差δmax=0.311 541。相對誤差在各區(qū)間的分布列表如下：

表6 DLR-F6模型壁面距離相對誤差分布范圍

由表6可見，壁面距離的相對誤差在1%以上相比起前2個算例增多，但網(wǎng)格單元總數(shù)是前2個算例的2倍以上，相對誤差較大的單元占網(wǎng)格單元總數(shù)的比例并未變大。圖12給出了翼展方向37.7%的截面Cp的計算值與實驗值對比，由圖可見計算值和實驗值吻合良好。該算例流場計算結(jié)果表明，對于大型復(fù)雜非結(jié)構(gòu)網(wǎng)格，逐層擴散法計算的壁面距離精度較好，具有良好的適應(yīng)性。

3 結(jié) 論

本文提出了一種計算流場網(wǎng)格單元到壁面最短距離的逐層推進算法，并給出了該算法的原理和求解過程。該算法主要針對非結(jié)構(gòu)網(wǎng)格由于拓撲結(jié)構(gòu)的隨意性導(dǎo)致壁面距離求解困難的問題。通過3個典型算例對本文提出的壁面最短距離逐層推進算法進行了考察和驗證計算。結(jié)果表明，本文算法可有效處理任意復(fù)雜的邊界，既可用于非結(jié)構(gòu)網(wǎng)格的計算，也適用于結(jié)構(gòu)網(wǎng)格，具有很好的普適性。該算法較枚舉法的計算效率高2～3個量級，精度完全滿足湍流模型的求解要求。隨著壁面復(fù)雜程度的提升，該算法相比起枚舉法的優(yōu)越性會愈發(fā)明顯。

參考文獻：

[1] David C W. Turbulence Modeling for CFD[M]. 2nd Edition. DCW Industries, 2000:30-39

[2] 徐汝鋒，陳志同，陳五一. 計算點到曲面最短距離的網(wǎng)格法[J]. 計算機集成制造系統(tǒng), 2011, 17(1): 95-100

Xu Rufeng, Chen Zhitong, Chen Wuyi. Grid Algorithm for Calculating the Shortest Distance from Spatial Point to Free-Form Surface[J]. Computer Integrated Manufacturing Systems, 2011, 17(1): 95-100 (in Chinese)

[3] 蘇智劍，吳序堂，毛世民.遺傳算法在求解空間曲線與曲面間最短距離中的應(yīng)用[J].機械設(shè)計與制造，2003(6):56-57

Su Zhijian, Wu Xutang, Mao Shimin. Genetic Algorithms for Solving the Minimum Distance between Bezier Curves and Surfaces[J]. Machinery Design and Manufacture, 2003(6): 56-57 (in Chinese)

[4] 陳麗萍，陳燕，胡德金.一種快速完備的自由曲線和曲面間最短距離求取算法[J]. 上海交通大學(xué)學(xué)報，2003,37(Suppl 2)：41-44

Chen Liping, Chen Yan, Hu Dejin. A High-Efficiency Algorithm to Calculate the Shortest Distance between Free-Form Curve and Free-Form Surface[J]. Journal of Shanghai Jiaotong University, 2003,37(Suppl 2)：41-44 (in Chinese)

[5] Renato N E, Marcos A D, Alvaro L G. Simple Finite Element-Based Computation of Distance Functions in Unstructured Grids[J]. Int J Numer Meth Engng, 2007, 72: 1095-1110

[6] Mauro T, David A S. Calculating Particle-to-Wall Distances in Unstructured Computational Fluid Dynamic Models[J]. Applied Mathematical Modelling, 2001, 25: 803-814

[7] 李素循. 典型外形高超聲速流動特性[M]. 北京：國防工業(yè)出版社， 2007

Li Suxun. Hypersonic Flow Characteristics of Typical Appearance[M]. Beijing: National Defense Industry Press, 2007 (in Chinese)

[8] Rumsey C L, Slotnick J P, Long M. Summary of the First AIAA CFD High-Lift Prediction Workshop[J]. Journal of Aircraft, 2011,48(6):2068-2079

[9] Hemsch M J, Morrison J H. Statistical Analysis of CFD Solutions from 2nd Drag Prediction Workshop[R]. AIAA-2004-556

[10] 王剛，葉正寅.三維非結(jié)構(gòu)混合網(wǎng)格生成與N-S方程求解[J]. 航空學(xué)報，2003,24(5)：385-390

Wang Gang, Ye Zhenyin. Generation of Three Dimensional Mixed and Unstructured Grids and Its Application in Solving Navier Stokes Equations[J]. Acta Aeronautica et Astronautica Sinica, 2003,24(5)：385-390 (in Chinese)