基于徑向基函數(shù)和Delaunay圖映射的高效高魯棒性的非結(jié)構(gòu)網(wǎng)格變形方法

王昊達, 劉 南, 張 穎, 崔曉春

(中國航空工業(yè)空氣動力研究院, 遼寧沈陽 110034)

引 言

在使用數(shù)值模擬方法求解復(fù)雜流場時, 須對控制方程進行空間離散和時間離散, 然后求解離散后的線性方程組得到計算網(wǎng)格上的流場變量。在對空間進行離散時往往須將計算域切分為若干計算網(wǎng)格。當(dāng)計算域的邊界運動或者變形時, 須利用網(wǎng)格變形方法生成新的網(wǎng)格。網(wǎng)格變形效率和魯棒性是限制其應(yīng)用的關(guān)鍵因素, 尤其是在氣動外形優(yōu)化、 氣動彈性計算等涉及大規(guī)模網(wǎng)格變形的問題中[1-2]。

根據(jù)采用的計算模型的不同, 將網(wǎng)格變形方法可大致分為兩類: 物理模型法、 代數(shù)插值法[1-2]。物理模型法將網(wǎng)格中的點、 線和單元視為可變形的實體, 考慮計算網(wǎng)格節(jié)點之間的連接關(guān)系, 可以有效地避免網(wǎng)格單元交叉。物理模型法主要包含: 彈簧類方法和彈性體方法等。彈簧類方法的應(yīng)用比較廣泛, 首先由Batina[3]提出, 其基本思想是將計算域內(nèi)兩節(jié)點間的連接線看作一根拉壓彈簧, 將整個網(wǎng)格看作一個彈簧系統(tǒng), 通過求解該彈簧系統(tǒng)的靜平衡方程求解節(jié)點位移。Farhat等[4]隨后提出了一種扭轉(zhuǎn)彈簧法, 用以防止網(wǎng)格單元的交叉。彈性體法最早是由Tezduyar[5]提出的, 作為彈簧類法的延伸, 其核心是將整個計算域內(nèi)的網(wǎng)格看作一整個彈性體, 通過對線性彈性方程組的求解來獲得各網(wǎng)格節(jié)點的位移。由于彈性體法要耗費大量時間迭代求解大型的控制方程組, 所以其變形效率很低, 因此工程上不適用于網(wǎng)格數(shù)量很多的復(fù)雜外形以及三維問題。代數(shù)插值法是指用某種數(shù)學(xué)插值方法將物面邊界的變化插值到內(nèi)部節(jié)點上, 其優(yōu)勢在于不需要計算域內(nèi)部網(wǎng)格節(jié)點間的連接信息, 只需要網(wǎng)格節(jié)點坐標(biāo)。主要包含: 超限插值法(transfinite interpolation, TFI)、 反距離加權(quán)函數(shù)插值法(inverse distance weighted, IDW)、 徑向基函數(shù)插值法以及基于Delaunay圖映射方法。TFI方法最早由Gaitonde等[6]提出, 雖然其變形效率高、 網(wǎng)格質(zhì)量好但僅能應(yīng)用于結(jié)構(gòu)網(wǎng)格。Liu等[7]提出的DGM方法, 通過構(gòu)建背景網(wǎng)格與計算網(wǎng)格的映射關(guān)系求解變形后的網(wǎng)格節(jié)點坐標(biāo)。該方法效率較高但不適用于大變形。與之相反, Rendall等[8-9]提出的RBF方法以及Witteveen等[10]提出的IDW方法在大變形時能夠保證較高的網(wǎng)格質(zhì)量但計算效率很低[11]。

目前的各類網(wǎng)格變形方法很難同時滿足網(wǎng)格變形效率和網(wǎng)格變形質(zhì)量的需求[2], 因此本文通過網(wǎng)格聚合算法自動建立粗化網(wǎng)格作為背景網(wǎng)格; 通過RBF方法建立背景網(wǎng)格和物面計算網(wǎng)格之間的代數(shù)關(guān)系; 并通過DGM方法建立背景網(wǎng)格和所有計算網(wǎng)格之間的映射關(guān)系, 結(jié)合RBF方法變形能力強和DGM方法變形效率高的優(yōu)點, 建立一種高效、 高魯棒性的網(wǎng)格變形方法。最后通過不同算例驗證本方法的網(wǎng)格變形能力和效率, 以及不同網(wǎng)格粗化參數(shù)對網(wǎng)格變形能力的影響。

1 方法描述

1.1 RBF&DGM網(wǎng)格變形方法

利用網(wǎng)格聚合方法結(jié)合RBF和DGM方法的優(yōu)點, 建立RBF&DGM網(wǎng)格變形方法。方法的主要思路是: 通過網(wǎng)格聚合算法對計算網(wǎng)格進行粗化并自動生成一套背景網(wǎng)格, 然后結(jié)合物面計算網(wǎng)格點、 計算域邊界點以及4個極大點構(gòu)建Delaunay控制體, 最后建立Delaunay控制體與計算網(wǎng)格點之間的幾何映射關(guān)系。當(dāng)物面邊界發(fā)生變化時, 通過RBF插值方法更新背景網(wǎng)格, 物面計算網(wǎng)格通過對應(yīng)的變化規(guī)律進行更新, 空間計算網(wǎng)格結(jié)合更新后的背景網(wǎng)格點坐標(biāo)通過建立的幾何映射關(guān)系進行更新。由于網(wǎng)格粗化是從物面邊界出發(fā)一層層推進的, 因此可以更好地保證物面附近黏性網(wǎng)格的正交性, 并且除計算邊界網(wǎng)格點、 物面邊界網(wǎng)格點外, 背景網(wǎng)格中還存在空間計算網(wǎng)格點, 因此通過背景網(wǎng)格建立的幾何映射關(guān)系更為精確, 能更好地傳遞物面邊界的變化規(guī)律。為方便程序使用, 本文針對二維網(wǎng)格以及三維網(wǎng)格均直接構(gòu)造三維Delaunay控制體進行幾何映射關(guān)系的建立。

本文方法主要分為前處理和網(wǎng)格變形兩步, 前處理的作用是自動生成背景網(wǎng)格并建立Delaunay控制體以及Delaunay控制體與計算網(wǎng)格之間的幾何映射關(guān)系, 網(wǎng)格變形的作用是利用上述幾何映射關(guān)系和物面計算網(wǎng)格點的變化規(guī)律更新計算網(wǎng)格, 流程分別如圖1、 圖2所示, 圖中下標(biāo)f表示背景網(wǎng)絡(luò)中的物面網(wǎng)格節(jié)點, 下標(biāo)s表示背景網(wǎng)格中的空間網(wǎng)格節(jié)點。

圖1 前處理步驟Fig. 1 Pre-processing steps

圖2 網(wǎng)格變形步驟Fig. 2 Grid deformation steps

前處理步驟如下:

1) 通過網(wǎng)格聚合方法建立粗網(wǎng)格作為背景網(wǎng)格;

2) 利用RBF方法建立背景網(wǎng)格中物面和空間網(wǎng)格點變形之間的關(guān)系;

3) 利用物面網(wǎng)格點、 計算域邊界點、 RBF控制點、 4個極大點建立Delaunay控制體;

4) 構(gòu)建Delaunay控制體與所有計算網(wǎng)格點之間的映射關(guān)系, 即得到每個空間網(wǎng)格點對應(yīng)的4個背景網(wǎng)格點, 以及相應(yīng)的權(quán)重系數(shù)。

網(wǎng)格變形步驟如下:

1) 根據(jù)輸入獲取物面變形量;

2) 利用RBF方法, 根據(jù)物面網(wǎng)格點變形得到背景網(wǎng)格中空間網(wǎng)格點的變形, 更新背景網(wǎng)格;

3) 利用前處理步驟中建立的幾何映射關(guān)系, 根據(jù)物面網(wǎng)格點變形和背景網(wǎng)格點變形計算得到所有空間網(wǎng)格點變形, 更新計算網(wǎng)格。

1.2 網(wǎng)格聚合方法

混合網(wǎng)格的聚合方法較結(jié)構(gòu)網(wǎng)格復(fù)雜很多, 本文參考Berglind[12]提出的方法。在該算法中, 強制將前沿陣面從物面邊界處逐層向外擴展, 通過設(shè)置每個網(wǎng)格節(jié)點的優(yōu)先級來完成(該優(yōu)先級表征每個節(jié)點與壁面的最小連接數(shù)), 以確保每一個附面層網(wǎng)格單元每一層頂點都具有相同的優(yōu)先級。

如圖3所示, 所有節(jié)點的優(yōu)先級由以下算法來確定:

圖3 網(wǎng)格節(jié)點優(yōu)先級Fig. 3 Priority numbering of all vertices

1) 設(shè)置所有網(wǎng)格節(jié)點的優(yōu)先級為p=-1;

2) 設(shè)置所有物面邊界網(wǎng)格節(jié)點的優(yōu)先級p=0;

3) 創(chuàng)建與所有優(yōu)先級為p的網(wǎng)格節(jié)點相鄰的邊構(gòu)成的前沿陣面;

4) 將所有與非負(fù)優(yōu)先級節(jié)點相鄰節(jié)點的優(yōu)先級加1,p=p+1;

5) 返回第3)步, 直到所有網(wǎng)格節(jié)點的優(yōu)先級均非負(fù)。

優(yōu)先級用于種子點選擇的標(biāo)準(zhǔn)。以下選擇標(biāo)準(zhǔn)按給定順序, 從前沿陣面的候選節(jié)點中選擇種子點:

1) 具有最低優(yōu)先級數(shù)p的網(wǎng)格節(jié)點。

2) 與已聚合的控制體間連接數(shù)最大的網(wǎng)格節(jié)點。

為了保證聚合后的粗網(wǎng)格物面附近的正交性, 本算法中規(guī)定了不同區(qū)域的不同定向聚合策略: 在附面層網(wǎng)格區(qū)域中(二維為四邊形網(wǎng)格單元、 三維為棱柱形網(wǎng)格單元), 網(wǎng)格聚合的方向只在法向方向而切向不聚合。在附面層網(wǎng)格以外的區(qū)域, 長寬比將作為定向聚合的判據(jù)。長寬比被定義為三維V/S1.5, 二維V/L2。其中V是控制體的體積,S是表面積,L是控制體積的周向長度。

粗化比R用于定義與粗網(wǎng)格控制體聚合的細(xì)網(wǎng)格控制體的數(shù)量。Rp、Rt分別表示具有棱柱狀元素和四面體的區(qū)域的粗化比, 二者在大多數(shù)情況下會有所不同。因此本文通過全局平均粗化比來控制網(wǎng)格粗化程度。

在以下的描述中, “點”特指對偶網(wǎng)格的控制單元。具體網(wǎng)格聚合算法可概述如下:

1) 在初始前沿陣面中根據(jù)標(biāo)準(zhǔn)選擇一個“種子點”(初始聚合網(wǎng)格)。

2) 檢查種子點是否位于邊界層內(nèi)。

3) 若種子點位于邊界層內(nèi), 則只與優(yōu)先級比先前聚集的節(jié)點高一個單位的相鄰細(xì)網(wǎng)格點沿法向聚合。若種子點位于邊界層之外, 則從其相鄰的控制體中聚合使其長寬比最大化的點, 生成一個粗網(wǎng)格控制單元。

4) 更新該粗網(wǎng)格控制單元的相鄰細(xì)網(wǎng)格點列表, 重復(fù)步驟3)直至達到所需的網(wǎng)格粗化比。

5) 更新前沿陣面, 將與最后被聚合的細(xì)網(wǎng)格點相鄰的未被聚合的點加入到前沿陣面。

6) 若仍有細(xì)網(wǎng)格點未完全聚合, 則重復(fù)步驟1), 直至所有細(xì)網(wǎng)格點均被聚合。

7) 將所有的孤立點與和它相鄰的一個粗網(wǎng)格控制體進行融合(清除孤立點)。

對于同一初始網(wǎng)格, 本文方法可以進行多次網(wǎng)格聚合算法對初始網(wǎng)格進行粗化, 該過程是通過給定網(wǎng)格粗化層數(shù)l確定的。

本文通過網(wǎng)格聚合算法得到粗化后的背景網(wǎng)格點如圖4所示。

(a) Original grid

1.3 RBF網(wǎng)格變形方法

RBF插值函數(shù)形式如下[13]

式中,s(x)是x處徑向基函數(shù)值,φ(‖·‖)是基函數(shù), 一般選擇Gauss或Wendland′s C2函數(shù),i代表RBF中心,xi是各中心點坐標(biāo)。對于網(wǎng)格變形問題,p(x)=0。令函數(shù)值等于物面計算網(wǎng)格點位移, 并寫成矩陣和向量形式

Δxb=Cbbax
Δyb=Cbbay
Δzb=Cbbaz

(1)

其中x方向(y和z方向與之類似)如下

因此,Cbb矩陣如下

其中徑向基函數(shù)為

φbibj=φ(‖xbi-xbj‖)

式中, 下標(biāo)b代表流場物面計算網(wǎng)格點,Nb為網(wǎng)格點數(shù), 令函數(shù)值等于空間點位移則有

Δxv=Avbax
Δyv=Avbay
Δzv=Avbaz

(2)

其中Avb矩陣如下

式中, 下標(biāo)v代表網(wǎng)格聚合后的背景網(wǎng)格點,Nv為網(wǎng)格點數(shù)。利用式(1)和(2)得到物面計算點和背景網(wǎng)格點位移之間的關(guān)系, 如下

(3)

1.4 DGM網(wǎng)格變形方法

首先在平面或空間上給定的包括流場邊界點和網(wǎng)格聚合后的點在內(nèi)的一組點的基礎(chǔ)上, 對整個計算域進行唯一的Delaunay圖三角化, 具體算法可參考Leatham[14], 從而完成計算域Delaunay控制體的構(gòu)建。隨后建立計算網(wǎng)格節(jié)點與Delaunay控制體節(jié)點之間的幾何映射關(guān)系。幾何映射關(guān)系通過相對面積(二維)和相對體積(三維)坐標(biāo)建立, 如圖5所示。結(jié)合背景網(wǎng)格, 通過DGM方法建立的Delauny控制體如圖6所示。

(a) Two-dimensional case

(a) Triangle

對于二維情況如圖5(a)所示, 任意網(wǎng)格節(jié)點O位于三角形單元△MNQ中, 由節(jié)點O和三角形組成的3個三角形△MON, △NOQ, △QOM的面積分別為S1,S2,S3, 令

對于三維情況如圖5(b)所示, 任意網(wǎng)格節(jié)點O位于四面體單元MNPQ中, 4個四面體OMNQ,OMQP,ONPQ,OMNP的體積分別為V1,V2,V3,V4, 令

式中,Si表示各頂點對應(yīng)三角形的面積,Vi表示各頂點對應(yīng)的四面體體積,S表示三角形單元面積,V表示四面體單元體積,λ表示表征幾何映射關(guān)系的定位參數(shù)。當(dāng)背景網(wǎng)格更新后, 根據(jù)定位參數(shù)可實現(xiàn)計算網(wǎng)格節(jié)點坐標(biāo)的更新。

對于二維網(wǎng)格有

對于三維網(wǎng)格有

(4)

式中,XO為任意計算網(wǎng)格節(jié)點坐標(biāo),X′為更新后的背景網(wǎng)格點坐標(biāo)。

通過式(4)可知, 當(dāng)完成背景網(wǎng)格的更新后, 通過λ即可計算出任意計算網(wǎng)格節(jié)點坐標(biāo), 完成計算網(wǎng)格的更新。

2 方法測試

通過算例對本文所建立的基于RBF&DGM的網(wǎng)格變形方法進行測試, 并與RBF、 DGM方法進行對比, 驗證本文方法的網(wǎng)格變形能力、 網(wǎng)格變形效率。同時改變不同網(wǎng)格粗化參數(shù), 研究不同網(wǎng)格粗化參數(shù)下生成的背景網(wǎng)格對網(wǎng)格變形質(zhì)量的影響。

2.1 網(wǎng)格質(zhì)量驗證標(biāo)準(zhǔn)

本文采用尺寸斜交度量函數(shù)fsize-skew定量衡量網(wǎng)格質(zhì)量[15], 計算公式如下

fsize=min(τ, 1/τ),τ=α/w

對于三角形單元

對于四面體單元

式中,fsize為相對尺寸度量函數(shù), 用于檢驗?zāi)繕?biāo)網(wǎng)格單元相比于參考網(wǎng)格單元的尺寸是否合適, 選取變形前網(wǎng)格單元作為參考單元, 則w為參考網(wǎng)格單元面積的2倍,α為目標(biāo)單元面積的2倍。當(dāng)fsize=1時, 表示目標(biāo)單元和參考單元面積相等, 當(dāng)fsize=0時, 表示目標(biāo)單元已經(jīng)退化。

fskew為形狀斜交度量函數(shù), 用于檢測網(wǎng)格單元形狀的扭曲程度, 對于三角形單元選取正三角形作為參考網(wǎng)格單元, 對于四面體單元選取正四面體作為參考網(wǎng)格單元,λ通過構(gòu)造張量度量矩陣獲得。當(dāng)fskew=1時, 表示目標(biāo)單元為正三角形或正四面體,fskew=0時, 表示目標(biāo)單元已退化。

fsize-skew為尺寸斜交度量函數(shù), 其取值范圍在0～1之間。當(dāng)單元退化或具有負(fù)面積/負(fù)體積時,fsize-skew值為0, 當(dāng)fsize-skew為1時, 說明此時變形后的網(wǎng)格為理想的正三角形或正多面體。

2.2 二維矩形算例

本算例測試網(wǎng)格變形方法在大變形下的適應(yīng)能力。采用長a=1 m, 寬b=0.5 m的矩形作為物面邊界, 正方形作為網(wǎng)格計算域, 網(wǎng)格數(shù)量9 100, 網(wǎng)格種類為非結(jié)構(gòu)網(wǎng)格且不生成黏性網(wǎng)格, 原始網(wǎng)格如圖7所示。

(a) Global view

繞矩形中心做旋轉(zhuǎn)變形, 旋轉(zhuǎn)角度θ以Δθ=10°, 從10°增加至90°。對3種方法變形后的網(wǎng)格進行質(zhì)量評估, 分別對比最小質(zhì)量系數(shù)和平均質(zhì)量系數(shù), 結(jié)果如圖8所示。

(a) Minimum quality coefficient

由圖8可知: 1) 由于矩形網(wǎng)格沒有生成黏性網(wǎng)格, 以正三角形作為斜交度量函數(shù)的參考網(wǎng)格單元, 3種方法下變形得到的網(wǎng)格的最小質(zhì)量系數(shù)和平均質(zhì)量系數(shù)均較高; 2) 隨著旋轉(zhuǎn)角度的增加, 最小網(wǎng)格質(zhì)量系數(shù)曲線快速下降, 平均網(wǎng)格質(zhì)量曲線變化較緩, 僅有小部分物面附近的網(wǎng)格因旋轉(zhuǎn)變形而受到擠壓, 大部分計算網(wǎng)格并未受到過多影響; 3) 本文方法與RBF方法的最小、 平均質(zhì)量系數(shù)曲線幾乎重合, 說明本方法具有與RBF方法相當(dāng)?shù)木W(wǎng)格變形能力。

圖9～圖11分別給出了物面旋轉(zhuǎn)角度θ=90°大變形下, 3種方法的網(wǎng)格變形情況。從圖中可以看出, 本文方法和RBF方法可以很好地將物面的變形傳遞至空間網(wǎng)格中, 而DGM方法在大變形時對物面變形的傳遞能力較差, 因此只有少部分計算網(wǎng)格受到物面變形的影響, 其平均網(wǎng)格質(zhì)量并未與其他兩種方法相差過大。

(a) Global view

(a) Global view

(a) Global view

2.3 網(wǎng)格粗化參數(shù)影響

背景網(wǎng)格的生成是本文方法的基礎(chǔ), 不同的網(wǎng)格粗化層數(shù)、 網(wǎng)格粗化比會影響到背景網(wǎng)格的密度, 進而影響幾何映射關(guān)系建立的精準(zhǔn)度。因此本節(jié)將研究網(wǎng)格粗化層數(shù)l、 網(wǎng)格粗化比R對網(wǎng)格變形質(zhì)量的影響。

圖12給出了矩形旋轉(zhuǎn)θ=90°后網(wǎng)格質(zhì)量隨網(wǎng)格層數(shù)l的變化曲線, 網(wǎng)格變形質(zhì)量以RBF方法作為參考, 徑向基函數(shù)半徑r取7.5。

圖12 不同l下的網(wǎng)格變形質(zhì)量Fig. 12 Grid deformation quality with different l

網(wǎng)格聚合算法共設(shè)置網(wǎng)格粗化層數(shù)5層, 隨著網(wǎng)格層數(shù)的增加, 背景網(wǎng)格密度逐漸變小, 分別取每一層粗化后的網(wǎng)格作為背景網(wǎng)格并建立Delaunay控制體及幾何映射關(guān)系。由圖12可以看出: 隨著l的不斷增大, 網(wǎng)格平均質(zhì)量幾乎與RBF方法一致, 網(wǎng)格最小質(zhì)量呈現(xiàn)先上升后下降的趨勢, 當(dāng)選取第5層粗化網(wǎng)格時, 最小網(wǎng)格質(zhì)量已經(jīng)接近于0, 說明此時已經(jīng)出現(xiàn)網(wǎng)格單元交叉情況。

圖12的結(jié)果表明, 背景網(wǎng)格的層數(shù)不宜過大也不宜過小, 應(yīng)選取較為適中的背景網(wǎng)格密度所對應(yīng)的網(wǎng)格層數(shù)l。當(dāng)l為1時, 此時背景網(wǎng)格密度較高, 背景網(wǎng)格單元尺寸較小, 建立的幾何映射關(guān)系過于精細(xì), 導(dǎo)致過度傳遞物面的旋轉(zhuǎn)變形, 使得物面附近的計算網(wǎng)格受擠壓更嚴(yán)重。因此, 最小網(wǎng)格質(zhì)量小于RBF方法。當(dāng)l為5時, 此時背景網(wǎng)格密度較低, 背景網(wǎng)格單元尺寸較大, 建立的幾何映射關(guān)系較為粗糙, 導(dǎo)致無法傳遞物面的旋轉(zhuǎn)變形, 使得出現(xiàn)網(wǎng)格單元交叉的情況。因此在選取網(wǎng)格粗化層數(shù)時, 應(yīng)當(dāng)選擇背景網(wǎng)格密度適中的粗化網(wǎng)格層數(shù)建立背景網(wǎng)格。

但單一地通過網(wǎng)格層數(shù)選取合適的背景網(wǎng)格密度有時是很難實現(xiàn)的, 因此還可以通過調(diào)節(jié)網(wǎng)格粗化比R來進一步選取合適的背景網(wǎng)格密度。選定第3層粗化網(wǎng)格研究R對網(wǎng)格質(zhì)量的影響, 見圖13。

圖13 不同R下的網(wǎng)格變形質(zhì)量Fig. 13 Grid deformation quality with different R

粗化比R表示與粗網(wǎng)格控制單元聚合的細(xì)網(wǎng)格控制單元的數(shù)量,R越大說明網(wǎng)格粗化能力越強、 背景網(wǎng)格的密度越小, 圖14中展示了不同R下的背景網(wǎng)格。由圖13可以看出, 隨著R的增加平均網(wǎng)格質(zhì)量先增大后減小, 最小網(wǎng)格質(zhì)量快速減小, 且當(dāng)R大于2.5時最小網(wǎng)格質(zhì)量約等于0, 出現(xiàn)網(wǎng)格單元交叉的情況。

(a) R=2

圖13的結(jié)果表明, 在固定網(wǎng)格層數(shù)的情況下, 當(dāng)R<2.75時, 隨著R增大, 生成的背景網(wǎng)格密度逐漸變小, 導(dǎo)致對物面變形的傳遞能力逐漸減小, 最終出現(xiàn)網(wǎng)格單元交叉情況, 但空間網(wǎng)格受物面變形影響較小, 因此平均網(wǎng)格質(zhì)量逐漸增加。但當(dāng)R>2.75時, 此時背景網(wǎng)格密度過大, 導(dǎo)致物面附近第1層網(wǎng)格高度變大, 受物面變形影響導(dǎo)致更多的空間網(wǎng)格出現(xiàn)網(wǎng)格單元交叉情況, 因此平均網(wǎng)格質(zhì)量開始逐漸減小。

綜上, 在構(gòu)建背景網(wǎng)格時, 應(yīng)當(dāng)將網(wǎng)格粗化層數(shù)與網(wǎng)格粗化比結(jié)合考慮, 不同的網(wǎng)格層數(shù)對應(yīng)著不同的最佳網(wǎng)格粗化比。本文建議: 為方便調(diào)節(jié)參數(shù), 在進行網(wǎng)格粗化時可選取第3層粗化網(wǎng)格作為基礎(chǔ), 若背景網(wǎng)格密度過大則適當(dāng)減小網(wǎng)格粗化比R, 若背景網(wǎng)格密度過小, 則適當(dāng)增大網(wǎng)格粗化比R, 進而選取更為合適的背景網(wǎng)格密度。

2.4 二維翼型算例

本算例測試本文方法對黏性網(wǎng)格物面正交性的保持能力。采用NACA64A010翼型, 網(wǎng)格數(shù)量為56 000, 計算域為r=30l0(l0為翼型弦長)的圓形域, 網(wǎng)格種類為非結(jié)構(gòu)網(wǎng)格, 在物面處生成25層黏性網(wǎng)格, 原始網(wǎng)格如圖15所示。

(a) Global view

繞翼型1/4弦線處作旋轉(zhuǎn)變形, 旋轉(zhuǎn)角度θ以Δθ=10°, 從0°順時針旋轉(zhuǎn)至90°, 并對3種方法變形后的網(wǎng)格進行質(zhì)量評估, 分別對最小質(zhì)量系數(shù)和平均質(zhì)量系數(shù)進行對比。翼型旋轉(zhuǎn)變形通過式(5)定義

(5)

由圖16可知: 1) 本文所用網(wǎng)格質(zhì)量系數(shù)fsize-skew中的斜交度量函數(shù)是以正三角形作為參考網(wǎng)格單元, 黏性網(wǎng)格的存在, 使得3種變形方法變形后的最小網(wǎng)格質(zhì)量系數(shù)和平均網(wǎng)格質(zhì)量系數(shù)均較小。2) 隨著旋轉(zhuǎn)角度θ的增加, DGM方法的網(wǎng)格質(zhì)量曲線快速下降, 但RBF和本文方法的網(wǎng)格質(zhì)量曲線變化較緩, 表明DGM方法物面附近黏性網(wǎng)格無法隨物面運動, 黏性網(wǎng)格單元受擠壓影響較大, 而另外兩種方法能夠使黏性網(wǎng)格跟隨物面一起運動。3) 本文方法與RBF方法的網(wǎng)格質(zhì)量曲線基本相同, 說明本文方法對于黏性網(wǎng)格同樣具有較好的網(wǎng)格變形能力。

(a) Minimum quality coefficient

圖17～圖19分別給出了翼型順時針旋轉(zhuǎn)至45°時3種方法下的網(wǎng)格變形情況。可以看出: 1) 在物面附近, 本文方法能夠與RBF方法保持一致, 較好地保證黏性網(wǎng)格的正交性; 2) 本文方法能夠較好地將物面變形傳遞至空間網(wǎng)格, 因此遠離物面的空間計算網(wǎng)格的網(wǎng)格質(zhì)量也要好于DGM方法并與RBF方法保持一致, 進一步保證了網(wǎng)格質(zhì)量。

(a) Global view

(a) Global view

(a) Global view

2.5 三維HIRENASD算例

以HIRENASD翼身組合體的機翼上彎和扭轉(zhuǎn)變形為例驗證本方法在三維復(fù)雜外形問題中的變形效率和能力。由于DGM方法的網(wǎng)格變形能力較差, 對于HIRENASD翼身組合體在彎曲小變形情況下網(wǎng)格就已出現(xiàn)負(fù)體積, 因此本算例僅與RBF方法進行對比。機翼彎曲和扭轉(zhuǎn)測試變形分別如下

式中, Δz(y)為機翼沿z向的彎曲變形位移沿展向分布,θ(y)為機翼繞前緣的扭轉(zhuǎn)角度沿展向分布,y為機翼展向,b為機翼展長,γ為幅度系數(shù)。

該構(gòu)型的計算域共包含網(wǎng)格節(jié)點個數(shù)為246 043。圖20為原始網(wǎng)格拓?fù)浜鸵砩姨幙臻g網(wǎng)格細(xì)節(jié)。

(a) Global view

(a) Global view

(a) Global view

(a) Minimum quality coefficient

(a) Minimum quality coefficient

設(shè)置彎曲變形系數(shù)γw=1和扭轉(zhuǎn)變形系數(shù)γt=0.5, 彎曲和扭轉(zhuǎn)變形后的網(wǎng)格分別如圖21、 圖22所示, 彎曲變形會導(dǎo)致附面層網(wǎng)格正交性下降, 圖23、 圖24也表明著彎曲變形對網(wǎng)格正交性的影響明顯大于扭轉(zhuǎn)變形。由上述幾幅圖可知, 本文方法和RBF方法的最小質(zhì)量系數(shù)曲線和平均質(zhì)量系數(shù)曲線基本吻合, 表明本文方法在保持網(wǎng)格質(zhì)量方面與RBF方法基本一致。

表1為本文方法與RBF方法對HIRENASD翼身組合體的網(wǎng)格變形計算時間對比。為避免誤差, 分別采用RBF方法以及本文方法執(zhí)行50次網(wǎng)格變形程序并對計算時間取平均值, 結(jié)果如表1所示, 整體網(wǎng)格計算時間由網(wǎng)格前處理時間、 網(wǎng)格變形時間兩部分構(gòu)成, 前處理的時間花費是一次性的, 通過前處理得到的Delaunay和RBF映射關(guān)系會存儲到文件或內(nèi)存中, 在后續(xù)計算中只需反復(fù)調(diào)用即可, 而網(wǎng)格變形時間的花費是需要反復(fù)疊加的, 因此本方法所花費的網(wǎng)格變形時間要遠低于RBF方法。在單次計算中本方法的計算效率提升可能并不顯著, 但對于大多數(shù)需要多次迭代更新網(wǎng)格的計算而言, 本方法的計算效率要遠高于RBF方法。且隨著網(wǎng)格數(shù)量的不斷增加, 計算效率提升的效果會更加明顯。

表1 計算時間對比(單位: s)

3 結(jié)論

本文建立了一種以網(wǎng)格聚合算法為基礎(chǔ), 結(jié)合徑向基函數(shù)插值和Delaunay圖映射的動態(tài)網(wǎng)格變形方法。通過二維、 三維網(wǎng)格算例對該方法進行了測試, 結(jié)果表明:

1) 背景網(wǎng)格的建立對網(wǎng)格變形質(zhì)量有著明顯影響, 在建立背景網(wǎng)格時應(yīng)綜合考慮網(wǎng)格粗化層數(shù)以及網(wǎng)格粗化比對背景網(wǎng)格密度的影響。在選取合適的背景網(wǎng)格密度時, 本文方法具有與標(biāo)準(zhǔn)RBF方法接近的網(wǎng)格變形質(zhì)量。

2) 該方法針對于無黏網(wǎng)格和黏性網(wǎng)格均具有較好的網(wǎng)格變形能力, 可以較好地保證黏性網(wǎng)格物面附近的網(wǎng)格正交性。

3) 該方法具有較高的網(wǎng)格變形效率, 網(wǎng)格點數(shù)量越大, 提升越明顯, 針對于三維標(biāo)準(zhǔn)HIRENASD模型, 本文方法的變形效率較RBF方法提升90%以上。