王 碩， 龐宇飛， 肖素梅， 劉 楊， 齊 龍， 何雨陽， 謝冬香

(1.西南科技大學(xué) 制造科學(xué)與工程學(xué)院,綿陽 621010； 2.中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室,綿陽 621000)

1 引 言

隨著計算機(jī)科學(xué)的發(fā)展，計算力大幅度提高，使得計算流體力學(xué)(computational fluid danamics)在各個行業(yè)中變得越來越重要，而網(wǎng)格生成在CFD中扮演著重要的角色，且是數(shù)值計算的重要步驟。網(wǎng)格中附面層網(wǎng)格是為了更好地捕捉近物面薄層的流動物理特性，影響著CFD計算的精度。目前生成附面層網(wǎng)格的自動性及可靠性較差，而且需要大量的人工交互來完成，使得附面層網(wǎng)格的生成費時費力。近年來，很多學(xué)者采用混合網(wǎng)格的策略生成半結(jié)構(gòu)化的棱柱網(wǎng)格，以完成附面層網(wǎng)格的生成[1-5]。

前沿層進(jìn)法[6-10]為經(jīng)典附面層網(wǎng)格生成方法。該方法先選取物面網(wǎng)格為前沿層，然后基于該前沿層逐層沿法向推進(jìn)，生成的網(wǎng)格作為前沿層，繼續(xù)重復(fù)上述操作。文獻(xiàn)[11]基于前沿層進(jìn)法，發(fā)展了一套附面層網(wǎng)格生成方法，其中包括前沿法向和附面層生成高度的調(diào)整方法。

張來平等[12]針對附面層三棱柱網(wǎng)格進(jìn)行了研究，提出了一種聚合的方法。首先生成各向異性的四面體網(wǎng)格，然后再將四面體網(wǎng)格進(jìn)行聚合，從而生成三棱柱網(wǎng)格。但是聚合后會有少數(shù)質(zhì)量較差的網(wǎng)格單元。孫巖[13]提出一種交互式棱柱網(wǎng)格生成方法,通過人工交互生成附面層邊界點的空間推進(jìn)網(wǎng)格，然后利用徑向基函數(shù)插值方法，得到附面層每一層內(nèi)部的網(wǎng)格點。該方法利用人工交互生成邊界點作為約束條件，有效避免了附面層網(wǎng)格交叉，并可以對局部網(wǎng)格進(jìn)行調(diào)整，最終生成高質(zhì)量的棱柱網(wǎng)格。但該方法需要人工交互生成附面層網(wǎng)格的邊界，且附面層推進(jìn)過程中會存在一定的誤差，需要后期進(jìn)行優(yōu)化。孫巖等[14]提出了一種基于約束框架的棱柱網(wǎng)格生成方法。首先，利用TFI(Transfinite Interpolation)生成粗背景結(jié)構(gòu)框架網(wǎng)格，再通過背景網(wǎng)格和表面網(wǎng)格的對應(yīng)關(guān)系插值得到附面層網(wǎng)格。該方法在附面層網(wǎng)格質(zhì)量控制和局部修改方面有較高的效率和顯著的優(yōu)勢，但背景框架網(wǎng)格生成比較困難。Lu等[15]提出一種結(jié)構(gòu)多塊附面層網(wǎng)格自動生成技術(shù),是一種由外而內(nèi)的附面層網(wǎng)格生成技術(shù)。首先，從表面網(wǎng)格中提取幾何特征，再基于幾何特征構(gòu)造附面層網(wǎng)格框架，然后在框架約束內(nèi)通過TFI生成結(jié)構(gòu)附面層網(wǎng)格。該方法能夠快速生成結(jié)構(gòu)附面層網(wǎng)格，但結(jié)構(gòu)網(wǎng)格的生成由于需要進(jìn)行網(wǎng)格拓?fù)湓O(shè)計，不僅費時，且對軟件使用者的經(jīng)驗有一定要求。

本文在前文的研究基礎(chǔ)上，針對非結(jié)構(gòu)附面層網(wǎng)格生成的自動性和局部網(wǎng)格質(zhì)量控制能力進(jìn)行研究，提出了一種棱柱附面層網(wǎng)格生成方法。

2 附面層棱柱網(wǎng)格生成方法

本文發(fā)展了文獻(xiàn)[15]的以分而治之的附面層網(wǎng)格生成思想，用構(gòu)造的方法生成非結(jié)構(gòu)附面層網(wǎng)格。該方法生成的網(wǎng)格精度和全局一致性較高，能夠有效避免復(fù)雜外形附面層網(wǎng)格局部及全局交叉現(xiàn)象。

2.1 方法思路

引進(jìn)自動化的表面網(wǎng)格分片算法[16]，自動將表面網(wǎng)格分成區(qū)域可控的多個網(wǎng)格面片，并以此生成多塊網(wǎng)格框架；利用結(jié)構(gòu)網(wǎng)格附面層自動生成技術(shù)[15]的框架線算法，并加以改進(jìn)，利用徑向基函數(shù)插值[17]方法生成框架線中的輪廓線；在框架線的基礎(chǔ)上，利用徑向基函數(shù)插值方法構(gòu)造附面層頂面網(wǎng)格；通過線性插值生成附面層頂層網(wǎng)格點與表面網(wǎng)格點之間的附面層內(nèi)部網(wǎng)格點，進(jìn)而生成附面層網(wǎng)格。避免了文獻(xiàn)[13]在外形變化劇烈的復(fù)雜條件下，精細(xì)附面層網(wǎng)格生成的插值精度擾動抑制難題。圖1為附面層網(wǎng)格生成算法的流程。

圖1 附面層網(wǎng)格生成示意圖

2.2 徑向基函數(shù)插值

徑向基函數(shù)插值在動網(wǎng)格中已經(jīng)有較好應(yīng)用。該方法通過選取基點之間的距離構(gòu)建插值基，能夠適應(yīng)復(fù)雜邊界。徑向基函數(shù)插值的一般表現(xiàn)形式為[18]

(1)

根據(jù)式(1)得到關(guān)于插值系數(shù)αi的線性方程組，表示成矩陣的形式為

(2)

求解矩陣方程(2)可以得到插值系數(shù)αi。然后通過式(1)能得到任一點r的函數(shù)值及位移。

2.3 生成流程

算法的詳細(xì)說明如下。

(1) 輸入表面網(wǎng)格R。

(2) 支撐線的生成。按照一定的幾何特性選取部分邊界點作為支撐點，如圖1(a)選取A，B，C和D四個支撐點。按照網(wǎng)格點的法向方向和輸入的附面層高度h，算出支撐線AE，BF，CG和DH?？梢酝ㄟ^平均支撐點周圍網(wǎng)格單元的法向得到每個支撐點的法向方向，

(3)

式中np為支撐點的法向方向，Nf為支撐點p相鄰網(wǎng)格單元總數(shù)，ni為第i個網(wǎng)格單元的法向。

該方法在個別點法向生成存在一定問題，這是附面層生成求法向所面臨的共同難題。

(3) 輪廓線的生成。定義徑向基函數(shù)插值的被插值函數(shù)f(r)為選取基點的位移，r為表面網(wǎng)格邊界點。點A和點B的位移變形量為h，將點A和點B作為徑向基函數(shù)插值的基點，首先計算出插值系數(shù)αi，然后將線AB上的非支撐網(wǎng)格點進(jìn)行插值計算，得到線AB上非支撐網(wǎng)格點的位移變形量，從而得到EF上對應(yīng)的網(wǎng)格點。同理得到FG，F(xiàn)H和HE上對應(yīng)的網(wǎng)格點，生成輪廓線上的網(wǎng)格點Ro。

對比文獻(xiàn)[15]輪廓線的生成是線性連接兩個相鄰的支撐線端點，這需要較好的網(wǎng)格邊界。而本文方法可以處理比較復(fù)雜的網(wǎng)格邊界的輪廓線生成，對網(wǎng)格拓?fù)湟筝^低。

(4) 附面層頂層網(wǎng)格點的生成。求解插值系數(shù)αi時，徑向基函數(shù)插值中f(r)為表面邊界網(wǎng)格點的位移ΔR，r為表面網(wǎng)格邊界點，已知N個邊界點Rb為基點，邊界點的位移變形量為

ΔR=Ro-Rb

(4)

利用徑向基函數(shù)插值可以得到求解插值系數(shù)的線性方程組

ΔRx=Mαx， ΔRy=Mαy， ΔRz=Mαz

(5～7)

在求解表面網(wǎng)格內(nèi)部點的位移時，被插值函數(shù)f(r)為待求值，內(nèi)部網(wǎng)格的位移表示為ΔV。

ΔV=Dα=DM-1ΔR=HΔR

(8)

式(8)寫成矩陣的形式為

(9)

求解矩陣方程(9)計算出表面網(wǎng)格內(nèi)部點的位移值ΔV，然后得到附面層頂層網(wǎng)格點Rt。把表面網(wǎng)格點的拓?fù)潢P(guān)系映射到附面層頂層網(wǎng)格點上，得到附面層頂層網(wǎng)格的拓?fù)潢P(guān)系，生成附面層頂層網(wǎng)格。

(5 )附面層內(nèi)部網(wǎng)格點的生成。連接RtR，將其按照一定的增長率分成k份，得到內(nèi)部網(wǎng)格點Rp，通過表面網(wǎng)格拓?fù)潢P(guān)系和附面層頂層網(wǎng)格的鏈接關(guān)系，生成三棱柱附面層網(wǎng)格。

2.4 算法測試

圖2所示為法向推進(jìn)法和構(gòu)造法的二維算法。圖2(a)表明法向推進(jìn)法在局部因推進(jìn)過高會出現(xiàn)網(wǎng)格交叉的現(xiàn)象，一般通過降低附面層高度來避免，而圖2(b)通過框架線的約束，能有效解決此問題。

圖2 法向推進(jìn)法和構(gòu)造法

3 算例驗證

通過4種典型幾何外形驗證本文提出的附面層網(wǎng)格生成方法，分別是V型槽、圓球、吊艙和F6機(jī)翼。測試在3.5 GHz主頻、16 GB內(nèi)存的計算機(jī)上進(jìn)行。

3.1 V型槽

V型槽是各種幾何形體中的一種常見結(jié)構(gòu)。如飛行器的機(jī)翼與機(jī)身連接處等。圖3為V型槽的附面層網(wǎng)格，圖中V型槽槽底附近的附面層網(wǎng)格不是沿著法向生長，而是在框架線和附面層頂層網(wǎng)格的約束下，沿著與法向有一定角度的方向生長，保證了附面層網(wǎng)格生成時不會出現(xiàn)交叉現(xiàn)象。

V型槽表面網(wǎng)格共有1074個網(wǎng)格點。設(shè)置20層附面層層數(shù)、0.001 mm的首層高度和1.2的增長率。附面層的網(wǎng)格單元最小角為22.5°，最大角為135.3°，生成時間為0.4 s。

圖3 V型槽的附面層網(wǎng)格

3.2 圓球

球面是最基本的幾何曲面。圖4所示為圓球的附面層網(wǎng)格，可以看出附面層網(wǎng)格的正交性較好，整體附面層網(wǎng)格質(zhì)量較好。

圓球表面網(wǎng)格共有5024個網(wǎng)格點。設(shè)置20層附面層生長層數(shù)、0.001 mm的首層高度和1.2的增長率。附面層的網(wǎng)格單元最小角為25.4°，最大角為124.3°，生成時間為1.2 s。

圖4 圓球的附面層網(wǎng)格

3.3 吊艙

吊艙是外凸面和內(nèi)凹面都存在的幾何形體。前緣是環(huán)形的凸曲面，后緣是個環(huán)形薄片。圖5(a)為吊艙表面網(wǎng)格和輪廓線。圖5(b,c)分別為吊艙的切片圖，圖5(d)為局部切片圖。從圖5(b～d)可以看出，采用本文方法生成的附面層網(wǎng)格，在吊艙特殊位置處，有效地避免了附面層網(wǎng)格生成過程中網(wǎng)格點法向量交叉的問題，提高了附面層網(wǎng)格質(zhì)量。

圖5 吊艙的附面層網(wǎng)格

吊艙表面網(wǎng)格共有19097個網(wǎng)格點。設(shè)置20層附面層生長層數(shù)、0.001 mm的首層高度和1.2的增長率。附面層的網(wǎng)格單元最小角為6.5°，最大角為171°，生成時間為21 s。

3.4 F6機(jī)翼

圖6(a)所示為F6機(jī)翼的表面網(wǎng)格和輪廓線。圖6(b)為F6機(jī)翼附面層的切片圖。圖6(c)為F6局部切片圖。從圖6(b,c)可以看出，通過本方法進(jìn)行附面層網(wǎng)格生成，在F6機(jī)翼的前后緣以及機(jī)翼邊界處，可以有效避免利用法向推進(jìn)法帶來的法向交叉問題。因為本方法選取少數(shù)網(wǎng)格點作為基點進(jìn)行法向生成，所以避免了交叉現(xiàn)象，較好地提高了網(wǎng)格質(zhì)量。

圖6 F6機(jī)翼的附面層網(wǎng)格

機(jī)翼表面網(wǎng)格共有12970個網(wǎng)格點。設(shè)置20層附面層生長層數(shù)、0.001 mm的首層高度和1.2的增長率。附面層的網(wǎng)格單元最小角為5.2°，最大角為161°，生成時間為31.4 s。

4 結(jié) 論

本文提出了一種基于框架的附面層棱柱網(wǎng)格生成方法，并通過幾種典型外形對該方法進(jìn)行了測試，成功生成了附面層網(wǎng)格，得到以下結(jié)論。通過構(gòu)建的框架線和附面層頂層網(wǎng)格點對附面層內(nèi)部網(wǎng)格點的約束，有效避免了附面層網(wǎng)格的交叉現(xiàn)象，保證了附面層網(wǎng)格的生成；表面網(wǎng)格發(fā)生變化時，僅需要對發(fā)生變化的局部區(qū)域重新生成附面層網(wǎng)格，即可快速重新生成整體的附面層網(wǎng)格，大幅降低了修改附面層網(wǎng)格的時間；大量減少了人工交互的工作，為進(jìn)一步附面層網(wǎng)格生成自動化提供了參考。

本文在生成附面層內(nèi)部的網(wǎng)格點時，運用了簡單線性插值。表面網(wǎng)格形狀比較復(fù)雜的時候，會造成附面層網(wǎng)格正交性較差的現(xiàn)象。將復(fù)雜的表面網(wǎng)格進(jìn)行合理網(wǎng)格分區(qū),達(dá)到復(fù)雜表面網(wǎng)格簡單化的目的，或者選擇一種合適插值方法能有效地解決這個問題，下一步工作將進(jìn)一步嘗試和完善。