三維非結(jié)構(gòu)自適應多重網(wǎng)格技術(shù)

許和勇，葉正寅

(西北工業(yè)大學翼型葉柵空氣動力學國防科技重點實驗室，陜西 西安 710072)

0 引言

非結(jié)構(gòu)網(wǎng)格舍去了網(wǎng)格節(jié)點的結(jié)構(gòu)性限制，易于控制網(wǎng)格單元的大小、形狀及網(wǎng)格點位置，對復雜外形具有很強的適應能力，這使得非結(jié)構(gòu)網(wǎng)格技術(shù)在20世紀80年代末和90年代初得到了迅速的發(fā)展［1］。非結(jié)構(gòu)網(wǎng)格中的一個點周圍的點數(shù)和單元數(shù)是不固定的，所以可以方便地進行自適應加密，以提高計算的精度，文獻［2-6］對非結(jié)構(gòu)網(wǎng)格的自適應算法進行了較深入的研究。為了使得加密后的網(wǎng)格能夠反映真實的邊界形狀，文獻［2-3，6］采用的方法一般都是輸入邊界的表面函數(shù)方法，如果邊界形狀比較復雜，則形成一個表面函數(shù)來表示邊界是比較麻煩的，所以本文發(fā)展了一種較為簡易、容易操作的基于邊界外形的自適應方法，通過生成一套初始極密表面網(wǎng)格來使得加密后的網(wǎng)格保持邊界的原始幾何形狀。

另外，非結(jié)構(gòu)網(wǎng)格節(jié)點和單元的無序性也導致了網(wǎng)格參數(shù)的存儲信息量增加，計算的速度下降，所以提高非結(jié)構(gòu)網(wǎng)格的計算效率也成了近些年來的研究熱點。采用多重網(wǎng)格方法是CFD計算中提高計算效率的重要措施之一，它可以大大減少計算機時。基于結(jié)構(gòu)網(wǎng)格的多重網(wǎng)格技術(shù)目前發(fā)展得已經(jīng)非常成熟［7-8］，而基于非結(jié)構(gòu)網(wǎng)格的多重網(wǎng)格方法則僅有少量的國內(nèi)外文獻［9-11］發(fā)表。本文將自適應后的網(wǎng)格作為多重網(wǎng)格，然后采用多重網(wǎng)格算法，以提高計算速度。

1 控制方程

采用三維Euler方程作為控制方程，其守恒型積分形式為:

對于理想氣體有:

2 基于幾何外形的自適應網(wǎng)格技術(shù)

自適應網(wǎng)格加密的一般步驟為:

(1)計算得到粗網(wǎng)格上的流場解;

(2)根據(jù)選定的自適應判據(jù)，對滿足細分條件的網(wǎng)格邊作細分標志，即在該邊中點增加一個網(wǎng)格點;

(3)檢查所有的三角形面，要求每一個三角形的三條邊只允許有一條邊加點，或一個面的三條邊都加點，如圖1所示;

(4)檢查所有的四面體單元，要求每一個四面體單元的六條邊只允許有一條邊或同一個三角面的三條邊或六條邊都加點，如圖2所示;

(5)對于邊界上的加密邊，要保證其加密點落在邊界上;

(6)計算細網(wǎng)格上的流場解;

(7)返回步驟2。

圖1 三角形面剖分的兩種情況Fig.1 Two cases of face division

圖2 四面體剖分的三種形式Fig.2 Three cases of tetrahedron division

流場內(nèi)部的網(wǎng)格單元邊可以按照圖1和圖2的剖分方法直接進行加密點的添加，但是處于邊界上的單元邊在添加加密點時，如果直接取邊的中點，則會出現(xiàn)如圖3(a)所示的情況，無法使得邊界的原始形狀得到保持。所以，網(wǎng)格自適應應該是基于初始幾何外形的自適應，所要添加的點要添加到初始幾何外形上而不是添加到已有的網(wǎng)格邊上，如圖3(b)所示。

圖3 邊界點加密示意圖Fig.3 Sketch map of boundary refinement

為了使得加密后的邊界點能夠落在邊界上，本文提出了一種簡單易行的方法，通過預先生成的一套初始極密表面網(wǎng)格來將邊界網(wǎng)格加密點近似投影到邊界上，具體步驟如下:

(a)在進行網(wǎng)格自適應求解之前，對整個邊界生成一套極密的表面網(wǎng)格，作為原始邊界的近似描述，如圖4，為了觀察方便，某機翼端面的網(wǎng)格已隱去，只保留翼型邊界上的點，即沒有編號的較小點;

(b)對初始的粗網(wǎng)格進行自適應剖分，初步確定所有的加密點，如圖4中的點A、B、C、D、E等即為端面粗網(wǎng)格的邊界節(jié)點;

(c)找出粗網(wǎng)格上所有邊界加密點(例如圖4中的點P、Q)，然后將點投影到密網(wǎng)格上，從密網(wǎng)格的三角形單元中插值確定修正后的加密點坐標，作為新的加密點坐標。例如，將點P沿著線段AB的中垂線方向投影，與密網(wǎng)格的邊界邊P1P2相交于點P'，則P'即為點P的新位置，對于點Q同理，如圖5所示。

這樣就將邊界加密點從網(wǎng)格邊的中點挪動到了邊界上了，使得自適應加密后的邊界保真性較好，從而提高了最終的流場計算精度。

圖4 初始極密表面網(wǎng)格與粗網(wǎng)格邊界示意圖Fig.4 Sketch map of the initial boundaries

圖5 邊界加密點向邊界上投影后的效果圖Fig.5 Sketch map of the boundaries after refinement

3 多重網(wǎng)格技術(shù)

多重網(wǎng)格算法需要有不同單元數(shù)目的多套網(wǎng)格，形成這些網(wǎng)格的方法主要有單獨生成［12］以及從最密網(wǎng)格聚合生成［13］兩種方法。單獨生成的多重網(wǎng)格在計算時的插值關(guān)系需要查找確定，具有一定的計算量，而且存在數(shù)據(jù)傳遞的精度問題;聚合法需要預先從初始密網(wǎng)格不斷聚合生成，聚合后的粗網(wǎng)格單元不規(guī)則性太大，對計算收斂性有一定影響。而本文利用自適應后的網(wǎng)格作為多重網(wǎng)格，不但粗細網(wǎng)格之間的對應關(guān)系很容易確定，而且每一層網(wǎng)格的質(zhì)量較好。

3.1 多重網(wǎng)格的確立

本文的多重網(wǎng)格建立的步驟為:在進行自適應加密過程中，對每一層的網(wǎng)格作標記，例如記初始的粗網(wǎng)格為第一層，第一次加密后形成的網(wǎng)格為第二層，直至得到第n層網(wǎng)格，同時記錄第1層與第2層，…，第n-1層與第n層的網(wǎng)格單元之間的對應和包含關(guān)系。這樣多重網(wǎng)格計算所需要的網(wǎng)格單元對應和包含關(guān)系就確定了，示意圖如圖6所示。

圖6 不同層的多重網(wǎng)格示意圖Fig.6 Sketch map of three multigrid levels

3.2 多重網(wǎng)格的算法

時間推進中，在細網(wǎng)格上的計算格式為:

其中αi為常系數(shù)，Rf為殘值量，V為單元體積，Wf為守恒變量。

在較粗網(wǎng)格上的計算格式為:

守恒變量Wf和殘值Rf由細網(wǎng)格插值到較粗網(wǎng)格上的公式為:

其中下標c表示粗網(wǎng)格量，f表示細網(wǎng)格量，Vf和Vc分別表示細網(wǎng)格和粗網(wǎng)格的單元體積。求和是對粗網(wǎng)格中包含的所有細網(wǎng)格進行。

粗網(wǎng)格的修正量向細網(wǎng)格傳輸時，細網(wǎng)格的修正量近似等于細網(wǎng)格所在粗網(wǎng)格的修正量。

4 算例

以M6機翼的跨聲速繞流計算作為算例，計算的來流馬赫數(shù)為0.8395，迎角為3.06°。圖7(a)給出的是初始粗網(wǎng)格的機翼表面網(wǎng)格，圖7(b)為邊界加密時要用到的表面極密網(wǎng)格放大圖。網(wǎng)格邊的加密條件是兩個節(jié)點的密度值的差大于計算域平均值的1.3倍。第一次自適應加密后得到第二套較密網(wǎng)格，然后應用多重網(wǎng)格算法得到第二套網(wǎng)格上的流場解，再進行自適應加密得到第三套網(wǎng)格，如此執(zhí)行，直至得到第四套網(wǎng)格，并用多重網(wǎng)格方法計算得到流場解。各套網(wǎng)格的數(shù)據(jù)如表1所示。

圖7 M6機翼的表面粗網(wǎng)格和極密網(wǎng)格圖Fig.7 Coarse grid and fine grid of the M6 wing surface

表1 M6機翼自適應多重網(wǎng)格情況Table1 The adaptive multigrids for M6 wing calculation

圖8給出的是三次自適應后第四套網(wǎng)格的表面網(wǎng)格圖，可以看出表面的“λ”形激波被捕捉得非常清晰。圖9中虛線表示在計算最后一次流場解時不使用多重網(wǎng)格的殘值收斂情況，而實線表示采用四重網(wǎng)格求解時的殘值收斂情況，橫坐標表示CPU時間，可見對于本算例來說，多重網(wǎng)格的加速比達到5.3左右，大大縮短了計算時間。

圖8 自適應后的機翼表面網(wǎng)格Fig.8 Surface grid after refinement

圖9 多重網(wǎng)格殘值收斂比較圖Fig.9 The comparison of the calculation residual

圖10給出的是機翼展向65%和90%處的表面壓力系數(shù)分布曲線圖。圖中虛線表示原始粗網(wǎng)格計算所得到的結(jié)果，而實線表示采用本文所發(fā)展的自適應方法后得到的結(jié)果。通過和實驗數(shù)據(jù)的比較可以看出，虛線與實驗結(jié)果相差較大，主要是因為機翼前緣的加密點沒有投影到機翼表面上，使得機翼前緣形狀“失真”，所以導致壓力極值點后移，且峰值偏低;而采用了本文的基于幾何外形的自適應技術(shù)后，機翼表面上網(wǎng)格邊加密點落在機翼的表面上，保持了機翼的原始形狀，計算結(jié)果與實驗值更加吻合。

5 結(jié)論

本文發(fā)展了一種基于幾何外形的非結(jié)構(gòu)自適應多重網(wǎng)格求解技術(shù)，對初始的粗網(wǎng)格根據(jù)加密判據(jù)進行剖分，得到與流場結(jié)構(gòu)相一致的較密網(wǎng)格，從而可以用較少的網(wǎng)格數(shù)來達到較高的求解精度。在進行網(wǎng)格自適應時，通過預先生成的初始極密表面網(wǎng)格來將表面網(wǎng)格邊上的加密點投影到邊界上，使得加密后的網(wǎng)格保持了邊界的初始外形。將自適應加密后得到的各層網(wǎng)格作為多重網(wǎng)格，采用多重網(wǎng)格算法進行流場求解，提高了計算效率。M6機翼的算例表明，本文所發(fā)展的方法比較準確地捕捉到流場中的激波，而且機翼表面的壓力系數(shù)分布與實驗值更加吻合。

圖10 機翼展向位置的表面壓力系數(shù)分布曲線Fig.10 Pressure coefficient distributions at different locations on the wing surface

［1］朱自強.應用計算流體力學［M］.北京:北京航空航天大學出版社，1998.

［2］CONNELL S D，HOLMES D G.A 3D unstructured adaptive multi-grid scheme for the Euler equations［R］.AIAA 1993-3339-CP，1993.

［3］HOLMES D G，CONNELL S D.Solution of the 2D Navier-Stokes equations on unstructured adaptive grids［R］.AIAA 89-1932，1989.

［4］王平，朱自強，拓雙芬，等.三維自適應非結(jié)構(gòu)網(wǎng)格的Euler方程解［J］. 航空學報，2001，22(6):495-499.

［5］朱培燁.三維Euler方程的自適應非結(jié)構(gòu)網(wǎng)格計算［J］.航空計算技術(shù)，2002，32(4):23-26.

［6］邱征，周磊，朱培燁.三維Euler方程非結(jié)構(gòu)自適應網(wǎng)格投影和光順技術(shù)研究［J］.航空計算技術(shù)，2006，36(1):79-85.

［7］SWANSON R C，TURKEL E.Multistage schemes with multi-grid for Euler and Navier-Stokes equations［R］.NASA TP-3631，1997.

［8］楊愛明，翁培奮，喬志德.用多重網(wǎng)格方法計算旋翼跨聲速無粘流場［J］.空氣動力學學報，2004，22(3):313-318.

［9］呂宏強，伍貽兆，夏健.非結(jié)構(gòu)聚合多重網(wǎng)格法流場數(shù)值模擬研究［J］.航空學報，2002，23(1):74-79.

［10］朱培燁.三維非結(jié)構(gòu)網(wǎng)格的歐拉方程聚合多重網(wǎng)格法［J］. 航空計算技術(shù)，2004，34(3):6-9.

［11］MOHAGNA J P，NEAL T F.Agglomeration multi-grid for unstructured grid flow solver［R］.AIAA 2004-759，2004.

［12］MAVRIPLIS D J.Three dimensional unstructured multigrid for the Euler equations［J］.AIAA Journal，1992，30(7):1753-1761.

［13］VENKATAKRISHNAN V，MAVRIPLIS D J.Agglomeration multigrid for the three dimensional Euler equations［R］.ICASE Report 94-5.Hampton，Virginia:NASA Langley Research Center，1994.