基于復(fù)雜網(wǎng)格處理的高精度數(shù)值積分技術(shù)

林義淼

摘要：目前，對飛機(jī)表面進(jìn)行網(wǎng)格劃分并進(jìn)行流體計(jì)算以優(yōu)化整體結(jié)構(gòu)的手段被廣泛運(yùn)用，然而由于外形復(fù)雜，采用笛卡爾網(wǎng)格劃分時(shí)會(huì)出現(xiàn)帶曲邊或曲面的復(fù)雜網(wǎng)格，若直接進(jìn)行數(shù)值積分則精度有限。該文提出了對復(fù)雜網(wǎng)格進(jìn)行映射變換處理后再進(jìn)行數(shù)值積分，并且數(shù)值實(shí)驗(yàn)表明采用該處理方法大大提高了精度。

關(guān)鍵詞：網(wǎng)格劃分；數(shù)值積分；復(fù)雜網(wǎng)格處理

中圖分類號(hào)：TP393 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1009-3044（2015）29-0179-02

High-accuracy Numerical Integration Technology Based on Processing Complex Grids

LIN Yi-miao

（School of Mathematics and System Sciences， Beihang University， Beijing 100191， China）

Abstract： Currently， the method of aircraft surface meshing and computational fluid dynamics to optimize the overall structure is widely used， however， due to complex shape， there exist complex grids with curved edge or surface when using Descartes grid. The accuracy is limited if doing the numerical integration directly. This paper proposes that processing complex grids with mapping transformation before numerical integration， and the numerical experiment shows that using the method of processing complex grids can greatly improve the accuracy.

Key words： meshing； numerical integration； processing complex grids

隨著數(shù)字化風(fēng)洞的不斷研究，使用計(jì)算機(jī)模擬飛機(jī)在風(fēng)洞中的表現(xiàn)，大大降低實(shí)驗(yàn)成本來優(yōu)化飛機(jī)結(jié)構(gòu)，將飛機(jī)表面進(jìn)行網(wǎng)格劃分并且進(jìn)行流體力學(xué)計(jì)算是在所難免的手段。

然而，由于飛機(jī)外形結(jié)構(gòu)復(fù)雜，采用笛卡爾網(wǎng)格對飛機(jī)表面進(jìn)行網(wǎng)格劃分時(shí)，所出現(xiàn)的數(shù)值積分區(qū)域就不是標(biāo)準(zhǔn)的積分區(qū)域了。例如二維翼型網(wǎng)格劃分后會(huì)出現(xiàn)帶曲邊邊界的三角形和四邊形網(wǎng)格，三維翼型網(wǎng)格劃分后會(huì)出現(xiàn)帶曲面的四面體網(wǎng)格。

目前，一般的數(shù)值積分技術(shù)中，對于帶曲邊邊界的網(wǎng)格是直接線段化處理，然而這樣做精度有限，而如果將網(wǎng)格加密處理，精度提高有限卻會(huì)增加計(jì)算負(fù)擔(dān)和降低計(jì)算效率，所以在本文將提出一種保留曲邊邊界進(jìn)行處理并且提高了數(shù)值積分精度的方法。

1 復(fù)雜網(wǎng)格介紹

網(wǎng)格劃分屬于流體力學(xué)計(jì)算中的預(yù)處理部分，是將計(jì)算區(qū)域劃分為較小的、不重疊的子域或單元（網(wǎng)格），從而得到有限單元及節(jié)點(diǎn)。而單元的數(shù)目決定了CFD的計(jì)算精度，網(wǎng)格的細(xì)密程度決定了數(shù)值解的精度和計(jì)算時(shí)間。

在采用笛卡爾網(wǎng)格進(jìn)行劃分時(shí)，由于飛機(jī)外形復(fù)雜，例如對飛機(jī)機(jī)翼表面進(jìn)行網(wǎng)格劃分時(shí)，會(huì)出現(xiàn)一邊為曲邊的四邊形網(wǎng)格（或三角形網(wǎng)格）和兩邊為曲邊的四邊形網(wǎng)格（或三角形網(wǎng)格），而如今隨著不斷研究，飛機(jī)外形也是日趨復(fù)雜化，相應(yīng)地對飛機(jī)表面進(jìn)行網(wǎng)格劃分所得到的網(wǎng)格也是日漸復(fù)雜。

當(dāng)然，若是網(wǎng)格劃分足夠密，每個(gè)網(wǎng)格足夠小，那么曲邊在一定程度上也就相當(dāng)于直邊了，所以在目前的數(shù)值積分技術(shù)中，都是采用直接將曲邊看成直線進(jìn)行處理，但如此做法必然會(huì)影響到最后的數(shù)值積分精度，所以本文提出一種基于復(fù)雜網(wǎng)格處理的數(shù)值積分技術(shù)，來提高數(shù)值積分的精度。

2 復(fù)雜網(wǎng)格處理

現(xiàn)假設(shè)在二維平面[?2]上有一個(gè)帶曲邊邊界的四邊形網(wǎng)格，要求解函數(shù)[fx，y]在其積分區(qū)域[Ω]上的數(shù)值積分，一般情況下，直接將曲邊看成直線邊進(jìn)行處理，利用上面介紹的數(shù)值積分方法（Simpson或Gauss方法等）求數(shù)值積分得到所要的結(jié)果，下面來介紹一下我們的處理方法。

首先建立一般化的直角坐標(biāo)，坐標(biāo)平面為[xOy]平面，四邊形網(wǎng)格的四個(gè)頂點(diǎn)坐標(biāo)所對應(yīng)的向量依次為[γ1，γ2，γ3，γ4]，其中向量[γ2，γ4]所對應(yīng)的頂點(diǎn)連接的正是網(wǎng)格的曲邊?，F(xiàn)有另一坐標(biāo)系[O-ξη]，其上有一個(gè)單位為1的正方形區(qū)域[Ω]，接下來找到一個(gè)一一映射：

[φ：?2（ξ，η）??2（x，y）]

將區(qū)域[Ω]一一映射到區(qū)域[Ω]上。

很自然地，有如下映射關(guān)系：

[（0，1）?γ1，（0，0）?γ2，（1，1）?γ3，（1，0）?γ4]

那么就可以構(gòu)造一個(gè)一一映射：

[γ（x，y）=γ1（1-ξ）η+γ2（1-ξ）（1-η）+γ3ξη+γ4ξ（1-η）]

則有：

[x=γ1（x）（1-ξ）η+γ2（x）（1-ξ）（1-η）+γ3（x）ξη+γ4（x）ξ（1-η）y=γ1（y）（1-ξ）η+γ2（y）（1-ξ）（1-η）+γ3（y）ξη+γ4（y）ξ（1-η）] （8）

根據(jù)數(shù)學(xué)分析中積分變換的內(nèi)容有如下關(guān)系：

[Ωf（x，y）dxdy=Ωf（ξ，η）J（φ）dξdη]

其中[Ω]表示原來曲邊邊界網(wǎng)格的區(qū)域，[Ω]表示單位為1的標(biāo)準(zhǔn)正方形區(qū)域，而[J（φ）]為映射[φ]的[Jacobi]矩陣，

[J=?x?ξ?x?η?y?ξ?y?η]

我們還可進(jìn)行如下處理：

根據(jù)一一映射關(guān)系，可以假設(shè)[?2（x，y）]中曲邊的曲線函數(shù)為[Γ（ξ）（0≤ξ≤1）]，那么在曲邊的兩個(gè)端點(diǎn)處，就有如下關(guān)系：

[γ2=Γ（0），γ4=Γ（1）]

則在式（8）中就有

[x=γ1（x）（1-ξ）η+Γ（0）（x）（1-ξ）（1-η）+γ3（x）ξη+Γ（1）（x）ξ（1-η）]

其中

[Γ（0）（x）（1-ξ）（1-η）+Γ（1）（x）ξ（1-η）=（1-η）[Γ（0）（x）（1-ξ）+Γ（1）（x）ξ]]

通過對[ξ]的取值討論，容易得到上式右邊等于[（1-η）Γ（ξ）（x）]，于是式（8）就可以寫成：

[x=γ1（x）（1-ξ）η+γ3（x）ξη+Γ（ξ）（x）（1-η）y=γ1（y）（1-ξ）η+γ3（y）ξη+Γ（ξ）（y）（1-η）]

對帶曲邊的網(wǎng)格區(qū)域進(jìn)行如上處理后再利用已有的數(shù)值積分技術(shù)進(jìn)行數(shù)值積分即可得到精度提高的數(shù)值結(jié)果了。

3 數(shù)值實(shí)驗(yàn)和結(jié)果

為了便于進(jìn)行實(shí)驗(yàn)，我們?nèi)×艘粋€(gè)相對簡單的算例，并且能得到精確結(jié)果的例子，最后能夠進(jìn)行多方對比來驗(yàn)證我們的方法。

現(xiàn)假設(shè)有一個(gè)帶曲邊的四邊形網(wǎng)格，并且曲邊函數(shù)為[f（x）=ex]，另外三邊分別為[x=0，y=0，x=1]，區(qū)域表示為[Ω]。

現(xiàn)在要求解[ΩF（x，y）dxdy]的數(shù)值積分結(jié)果，其中[F（x，y）=x2+y2]，我們都知道該結(jié)果完全可以精確求解出來：

[I=ΩF（x，y）dxdy=010exx2+y2dydx=19e3+e-199]

數(shù)值結(jié)果為[I≈2.8388970421]。

而在現(xiàn)有的數(shù)值積分技術(shù)中，都是忽略了曲邊，直接連接了曲邊的兩點(diǎn)，而后利用現(xiàn)有的數(shù)值積分方法進(jìn)行數(shù)值計(jì)算。我們選取了比較常用的Simpson方法和Gauss方法進(jìn)行數(shù)值積分，得到結(jié)果如下：

[IS=2.7954089752，IG=2.8133023936.]

下面我們先用本文第三部分介紹的方法對曲邊網(wǎng)格進(jìn)行處理后有：

[φ：?2??2：（ξ，η）?（ξ，ηeξ）]

則[Jacobi]矩陣為

[J（φ）=?x?ξ?x?η?y?ξ?y?η=10ηeξeξ=eξ]

則所求的數(shù)值積分變成了

[I=ΩF（x，y）dxdy=Ω（ξ2+η2e2ξ）eξdξdη=Ωξ2eξ+η2e3ξdξdη]

其中[Ω]表示標(biāo)準(zhǔn)正面形區(qū)域。

之后我們依然采用現(xiàn)有的Simpson方法和Gauss方法，并且取同樣的節(jié)點(diǎn)進(jìn)行數(shù)值積分計(jì)算，得到的結(jié)果為：

[I′S=2.8384170008，I′G=2.8388960361.]

我們可以對以上數(shù)據(jù)進(jìn)行簡單的處理和誤差分析，得到結(jié)果如表1所示。

其中[α]表示對于網(wǎng)格曲邊直接連接曲邊兩點(diǎn)做線段化處理，而[β]表示對網(wǎng)格曲邊進(jìn)行映射變換處理，處理方法見本文第三部分。

從上表可以很容易看出，對曲邊邊界進(jìn)行映射變換后再進(jìn)行數(shù)值積分的話，精確度有顯著提高，盡管這只是一個(gè)簡單算例，但也說明了我們對于帶曲邊的網(wǎng)格進(jìn)行適當(dāng)?shù)奶幚碓贁?shù)值積分的話，的確可以提高精度，相信這種處理方法和思想在高精度數(shù)值積分技術(shù)中會(huì)發(fā)揮其應(yīng)有的作用。

參考文獻(xiàn)：

[1] 徐利治，周蘊(yùn)時(shí).高維數(shù)值積分[M]. 科學(xué)出版社，1980.

[2] 徐利治，周蘊(yùn)時(shí)，何天曉.高維數(shù)值積分選講[M]. 安徽教育出版社，1985.

[3] Philip J.Davis， Philip Rabinowitz著，馮振興，伍富良譯.數(shù)值積分法[M].高等教育出版社，1986.

[4] Curtis F.Gerald， Patrick O.Wheatley.應(yīng)用數(shù)值分析[M]. 呂淑娟，譯.機(jī)械工業(yè)出版，2006.