王剛, 曾錚, 葉正寅

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

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

1 壁面逐層推進(jìn)法

1.1 算法原理

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

(1)

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

圖1 逐層推進(jìn)法原理示意圖

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

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

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

1.2 搜索算法的進(jìn)一步改進(jìn)

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

圖2 未改進(jìn)的算法可能遇到的問題示意圖

圖3 問題原因及解決方式

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

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

2 算例分析

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

(2)

(3)

平均相對(duì)誤差

(4)

最大相對(duì)誤差

(5)

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

2.1 類空天飛機(jī)模形

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

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

表1 2種算法計(jì)算時(shí)間比較(空天飛機(jī))

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

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

表2 類空天飛機(jī)模型壁面距離相對(duì)誤差分布范圍

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

圖6 子午線截面處Cp計(jì)算值與實(shí)驗(yàn)值對(duì)比圖

2.2 增升裝置模形

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

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

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

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

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

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

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

圖9 增升裝置翼展50%截面Cp計(jì)算值與實(shí)驗(yàn)值對(duì)比

2.3 DLR-F6翼身組合體模型

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

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

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

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

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

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

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

3 結(jié) 論

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

參考文獻(xiàn)：

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

[2] 徐汝鋒，陳志同，陳五一. 計(jì)算點(diǎn)到曲面最短距離的網(wǎng)格法[J]. 計(jì)算機(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].機(jī)械設(shè)計(jì)與制造，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é)報(bào)，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] 李素循. 典型外形高超聲速流動(dòng)特性[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é)報(bào)，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)