王昊達(dá), 劉 南, 張 穎, 崔曉春

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

引 言

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

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

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

1 方法描述

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

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

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

圖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)格點(diǎn)變形之間的關(guān)系;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

本文通過網(wǎng)格聚合算法得到粗化后的背景網(wǎng)格點(diǎn)如圖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是各中心點(diǎn)坐標(biāo)。對(duì)于網(wǎng)格變形問題,p(x)=0。令函數(shù)值等于物面計(jì)算網(wǎng)格點(diǎn)位移, 并寫成矩陣和向量形式

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

(1)

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

因此,Cbb矩陣如下

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

φbibj=φ(‖xbi-xbj‖)

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

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

(2)

其中Avb矩陣如下

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

(3)

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

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

(a) Two-dimensional case

(a) Triangle

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

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

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

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

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

(4)

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

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

2 方法測(cè)試

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

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

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

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

對(duì)于三角形單元

對(duì)于四面體單元

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

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

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

2.2 二維矩形算例

本算例測(cè)試網(wǎng)格變形方法在大變形下的適應(yīng)能力。采用長(zhǎng)a=1 m, 寬b=0.5 m的矩形作為物面邊界, 正方形作為網(wǎng)格計(jì)算域, 網(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°。對(duì)3種方法變形后的網(wǎng)格進(jìn)行質(zhì)量評(píng)估, 分別對(duì)比最小質(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)變形而受到擠壓, 大部分計(jì)算網(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方法在大變形時(shí)對(duì)物面變形的傳遞能力較差, 因此只有少部分計(jì)算網(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)格粗化比會(huì)影響到背景網(wǎng)格的密度, 進(jìn)而影響幾何映射關(guān)系建立的精準(zhǔn)度。因此本節(jié)將研究網(wǎng)格粗化層數(shù)l、 網(wǎng)格粗化比R對(duì)網(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)先上升后下降的趨勢(shì), 當(dāng)選取第5層粗化網(wǎng)格時(shí), 最小網(wǎng)格質(zhì)量已經(jīng)接近于0, 說明此時(shí)已經(jīng)出現(xiàn)網(wǎng)格單元交叉情況。

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

但單一地通過網(wǎng)格層數(shù)選取合適的背景網(wǎng)格密度有時(shí)是很難實(shí)現(xiàn)的, 因此還可以通過調(diào)節(jié)網(wǎng)格粗化比R來進(jìn)一步選取合適的背景網(wǎng)格密度。選定第3層粗化網(wǎng)格研究R對(duì)網(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)格粗化能力越強(qiáng)、 背景網(wǎng)格的密度越小, 圖14中展示了不同R下的背景網(wǎng)格。由圖13可以看出, 隨著R的增加平均網(wǎng)格質(zhì)量先增大后減小, 最小網(wǎng)格質(zhì)量快速減小, 且當(dāng)R大于2.5時(shí)最小網(wǎng)格質(zhì)量約等于0, 出現(xiàn)網(wǎng)格單元交叉的情況。

(a) R=2

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

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

2.4 二維翼型算例

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

(a) Global view

繞翼型1/4弦線處作旋轉(zhuǎn)變形, 旋轉(zhuǎn)角度θ以Δθ=10°, 從0°順時(shí)針旋轉(zhuǎn)至90°, 并對(duì)3種方法變形后的網(wǎng)格進(jìn)行質(zhì)量評(píng)估, 分別對(duì)最小質(zhì)量系數(shù)和平均質(zhì)量系數(shù)進(jìn)行對(duì)比。翼型旋轉(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)格無法隨物面運(yùn)動(dòng), 黏性網(wǎng)格單元受擠壓影響較大, 而另外兩種方法能夠使黏性網(wǎng)格跟隨物面一起運(yùn)動(dòng)。3) 本文方法與RBF方法的網(wǎng)格質(zhì)量曲線基本相同, 說明本文方法對(duì)于黏性網(wǎng)格同樣具有較好的網(wǎng)格變形能力。

(a) Minimum quality coefficient

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

(a) Global view

(a) Global view

(a) Global view

2.5 三維HIRENASD算例

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

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

該構(gòu)型的計(jì)算域共包含網(wǎng)格節(jié)點(diǎn)個(gè)數(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所示, 彎曲變形會(huì)導(dǎo)致附面層網(wǎng)格正交性下降, 圖23、 圖24也表明著彎曲變形對(duì)網(wǎng)格正交性的影響明顯大于扭轉(zhuǎn)變形。由上述幾幅圖可知, 本文方法和RBF方法的最小質(zhì)量系數(shù)曲線和平均質(zhì)量系數(shù)曲線基本吻合, 表明本文方法在保持網(wǎng)格質(zhì)量方面與RBF方法基本一致。

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

表1 計(jì)算時(shí)間對(duì)比(單位: s)

3 結(jié)論

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

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

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

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