本發(fā)明公開了一種基于埃爾米特插值基本加權無振蕩格式的全流場模擬方法，涉及計算流體力學工程技術領域。

背景技術：

在計算流體力學數(shù)值模擬中，高精度格式的構造及其應用一直是研究的重點內容，因為高精度格式能夠對全流場進行精確模擬，能夠很好的模擬出解的結構以及準確的捕捉激波位置。1983年，harten首次提出了tvd(totalvariationdiminishing)格式，并在此基礎上與osher于1987年提出了eno(essentiallynon-oscillatory)高精度格式，eno格式的主要思想是在逐次擴展的模板中選出最光滑的模板來構造多項式求出單元邊界處的值，進而達到高精度、高分辨率和在間斷處基本無振蕩的效果。但是，在方法的實現(xiàn)過程中，eno格式會造成計算結果的浪費，導致計算效率不高，因此，liu，osher和chan等提出了weno(weightedessentiallynon-oscillatory)格式，提高了計算結果的利用率并且使得r階精度的eno格式提高到r+1階精度。1996年，jiang和shu提出了一種新的weno格式，使得數(shù)值精度能夠提高到2r-1階。隨后，qiu和shu發(fā)展了此類方法，命名為hweno格式并將其作為限制器用于rkdg(runge-kuttadiscontinuousgalerkin)方法。來年，他們將這種方法推廣到二維情形。2008年，zhu和qiu構造出具有四階精度的有限體積hweno格式。這些經(jīng)典的數(shù)值計算方法對網(wǎng)格質量要求較高，不能直接應用于計算流場中包含復雜物面邊界的可壓繞流問題，另一方面，非結構網(wǎng)格上的高精度格式構造復雜，所以將高精度格式在笛卡爾網(wǎng)格中實現(xiàn)處理含復雜物面的全流場模擬問題是有必要的。對于這些含有復雜幾何外形物體的流場數(shù)值模擬，由于物面邊界與笛卡爾網(wǎng)格相交的任意性，所以在物面邊界處會產(chǎn)生解的偽振蕩。針對此問題，dadone等系統(tǒng)地研究了貼體網(wǎng)格中復雜幾何外形的物面邊界處理方法，將其命名為st(symmetrytechnique)方法和ccst(curvaturecorrectedsymmetrytechnique)方法，并將這些浸入邊界方法推廣到了非貼體的笛卡爾網(wǎng)格中，命名為gbcm(ghostbody-cellmethod)方法。另外，forrer等于1998年提出了一種浸入邊界法，命名為fgcm(forrer’sghostcellmethod)方法。這些方法均能有效處理復雜物面邊界條件，能有效降低所用計算網(wǎng)格的生成難度。

另外，一階精度的數(shù)值方法在捕捉激波時不會出現(xiàn)非物理的數(shù)值振蕩但會過度抹平強間斷，而往往強間斷對問題的后續(xù)研究有著重要意義，因此引進高精度數(shù)值計算格式進行相關數(shù)值模擬是有必要的。當計算區(qū)域中出現(xiàn)復雜物面時，許多基于結構網(wǎng)格設計出來的高精度計算格式不能直接應用于相關數(shù)值計算，需要先生成復雜的貼體網(wǎng)格，再構造出基于此貼體網(wǎng)格的高精度計算格式，步驟復雜繁瑣，推廣困難。

技術實現(xiàn)要素：

本發(fā)明的針對現(xiàn)有技術中的不足，提供一種基于埃爾米特插值基本加權無振蕩格式的全流場模擬方法，能針對各種復雜物體繞流問題，全流場統(tǒng)一采用笛卡爾網(wǎng)格和高精度計算格式進行數(shù)值模擬。本發(fā)明首先給出笛卡爾網(wǎng)格下有限體積robusthweno高精度數(shù)值計算格式的構造，此格式能更好地推廣到高維情形以及自適應和移動網(wǎng)格情形，隨后針對笛卡爾網(wǎng)格的非貼體特性，采用合適的邊界處理方法處理物面邊界條件，能將這二者有效結合并應用于具有復雜幾何外形物面的無粘可壓縮繞流問題的數(shù)值模擬。相比于weno格式，hweno格式通過hermite插值方法構造多項式，既需要函數(shù)值又需要其導數(shù)值，因此加強了各個網(wǎng)格單元之間的聯(lián)系，使得格式更加緊湊和穩(wěn)定，而robusthweno格式相對于hweno格式增強了格式的魯棒性，能夠更好地推廣。

為實現(xiàn)上述目的，本發(fā)明采用以下技術方案：

在笛卡爾坐標下，流場計算區(qū)域內含有物體，先通過虛擬單元浸入邊界法對物體內的網(wǎng)格點進行賦值形成全流場，然后利用robusthweno格式對全流場進行數(shù)值模擬，具體步驟包括：

1)利用虛擬單元浸入邊界法對物體內部的網(wǎng)格單元點進行賦值，使含有物體的流場形成單一介質的全流場；

2)建立無粘流體的控制方程，所述控制方程是一個關于時間變量t和關于空間變量x，y的方程組，對控制方程空間變量進行離散得到半離散的有限體積格式，具體步驟如下：

1.1分別對空間變量x，y求導得到方程組，增加控制方程的數(shù)量；

1.2對控制方程兩邊同時在目標單元上進行積分得到有限體積格式，所述有限體積格式包含時間導數(shù)項和空間積分項；

1.3對所述空間積分項利用gauss求積公式求解從而得到關于時間導數(shù)項的半離散有限體積格式；

3)對時間變量使用三階tvdrunge-kutta離散公式從而將半離散有限體積格式變成時空全離散有限體積格式；

4)根據(jù)時空全離散有限體積格式得到下一時間層上的流場值，依次迭代，得到全流場穩(wěn)定時的數(shù)值模擬。

為優(yōu)化上述技術方案，采取的具體措施還包括：

步驟1)中，利用st和fgcm兩種虛擬單元方法處理物面邊界條件，使得流場近似為單一介質的全流場，具體步驟包括：

a、先對計算區(qū)域內的點進行分類，分為物體內部的虛擬單元點和物體外部點；

b、找到所有虛擬單元點關于物面的對稱點，確定對稱點落入的最近網(wǎng)格單元i1，以及其上下左右的網(wǎng)格單元i2，i3，i4，i5；

c、利用網(wǎng)格單元i1，i2，i3，i4，i5上的流場值，通過插值公式，確定對稱點流場值；

d、通過st和fgcm虛擬單元方法，得到關于虛擬單元點和對稱點的物理量關系，從而得到虛擬單元點的流場值，將含有物體的流場形成成單一介質的全流場。

步驟1.1中，控制方程為：

其中，u＝(ρ，ρu，ρv，e)^t表示守恒變量，f(u)＝(ρu，ρu²+p，ρuv，u(e+p))^t，g(u)＝(ρv，ρuv，ρv²+p，v(e+p))^t，f(u)，g(u)表示通量，ut表示u對t求導，f(u)x表示f(u)對x求導，g(u)y表示g(u)對y求導，ρ，p，u，v，e分別表示流體密度、壓強、水平方向速度、豎直方向速度和能量，t表示轉置，u0表示初始狀態(tài)值；

對控制方程分別對空間變量x和y求導得到方程：

其中，表示u對x求導，表示u對y求導，h(u，v)x和q(u，w)x表示h(u，v)q(u，w)分別對x求導，r(u，v)y和s(u，w)y表示r(u，v)s(u，w)分別對y求導，h(u，v)＝f′(u)v，r(u，v)＝g′(u)v，q(u，w)＝f′(u)w，s(u，w)＝g′(u)w，f′(u)，g′(u)分別表示f(u)，g(u)對u求導。

步驟1.2中，對控制方程在目標單元iij內積分有：

其中，表示u對x求導，表示u對y求導，h(u，v)x和q(u，w)x表示h(u，v)q(u，w)分別對x求導，r(u，v)y和s(u，w)y表示r(u，v)s(u，w)分別對y求導，h(u，v)＝f′(u)v，r(u，v)＝g′(u)v，q(u，w)＝f′(u)w，s(u，w)＝g′(u)w，f′(u)，g′(u)分別表示f(u)，g(u)對u求導。

步驟1.3中，利用gauss求積公式求解，具體過程如下：

a、在目標單元邊界上取8個gauss點，分別為和并將目標單元及周圍共9個單元組成母模板；

b、在母模板中選取8個子模板，每個子模板必須包括目標單元；

c、利用hermite插值方法在8個子模板中求出目標單元每個gauss點的重構多項式pl(x，y)，l＝1，…，8，其中l(wèi)表示對應子模板序號，則有：

其中，分別代表u在網(wǎng)格單元ii-1，j-1、ii，j-1、ii+1，j-1、ii-1，j、iij、ii+1，j、ii-1，j+1、ii，j+1、ii+1，j+1的平均值，分別代表v在網(wǎng)格單元ii-1，j、ii+1，j的平均值，分別代表w在網(wǎng)格單元ii，j-1、ii，j+1的平均值，h為網(wǎng)格步長；

d、取線性權γl，計算光滑指示器βl，用于衡量重構多項式pl(x，y)，l＝1，…，8在目標單元上的光滑程度，計算公式為：

其中，α為多重指標，d是求導算子，r是多項式pl(x，y)的最高次數(shù)；

e、計算非線性權ωl，其計算公式為：

其中，為計算過程中的過渡值，βl為光滑指示器，ε＝10^-6；

f、求出u在每個gauss點gk處的近似表達式為：

其中gk表示第k個gauss點，k＝1，…，8；

g、重復步驟c到步驟f，求出v，w在每個gauss點gk處的近似表達式；

h、利用gauss求積公式求解：

其中，表示目標單元的邊界長度，ak為gauss求積公式的權重系數(shù)，為單位外法向向量，和用數(shù)值通量代替；

i、將計算結果代入含有時間導數(shù)項的半離散有限體積格式，得到關于時間導數(shù)的常微分方程。

步驟3)中，三階tvdrunge-kutta離散公式為：

其中，為中間過渡值，δt為時間步長，上標n表示第n時間層，為積分近似值。

步驟4)中，所述時空全離散有限體積格式為關于時間層的迭代公式，初始狀態(tài)值已知，通過迭代公式求出下一時間層的流場值，依次得到穩(wěn)定時的全流場數(shù)值模擬。

本發(fā)明的有益效果是：將笛卡爾網(wǎng)格下有限體積robusthweno格式與浸入邊界虛擬單元方法結合，應用于具有復雜幾何外形的無粘流動問題數(shù)值模擬中，這樣避免了復雜網(wǎng)格的生成和非結構網(wǎng)格下高精度格式的構造，而相比于weno格式和hweno格式，robusthweno格式增強了格式的緊湊性和魯棒性，并且可以很好地推廣到更高維的和移動網(wǎng)格下的數(shù)值模擬中。傳統(tǒng)的數(shù)值方法在解決具有真實斜面的雙馬赫反射問題，圓柱繞流問題和翼型繞流問題時需要將做一定的等價轉化甚至是無法解決的，但本發(fā)明能穩(wěn)健，準確，基本無振蕩地捕捉到強激波和接觸間斷并同時在解的光滑區(qū)域保持時空一致高階精度。

附圖說明

圖1是大模板t示意圖。

圖2是笛卡爾網(wǎng)格中網(wǎng)格單元與物面邊界相交示意圖。

圖3a-3b是實施例一中，前臺階問題分別采用st和fgcm兩種邊界處理方法的密度等值線圖。

圖4a-4b是實施例一中，前臺階問題分別采用st和fgcm兩種邊界處理方法的壓強等值線圖。

圖5a-5b是實施例二中，雙馬赫問題分別采用st和fgcm兩種邊界處理方法的密度等值線圖。

圖6a-6b是實施例二中，雙馬赫問題分別采用st和fgcm兩種邊界處理方法的壓強等值線圖。

圖7a-7b是實施例三中，圓柱繞流問題分別采用st和fgcm兩種邊界處理方法的密度等值線圖。

圖8a-8b是實施例三中，圓柱繞流問題分別采用st和fgcm兩種邊界處理方法的馬赫數(shù)等值線圖。

圖9a-9b是實施例四中，naca0012翼型繞流問題分別采用st和fgcm兩種邊界處理方法的壓強等值線圖。

圖10a-10b是實施例四中，naca0012翼型繞流問題分別采用st和fgcm兩種邊界處理方法的翼型表面壓力系數(shù)圖。

具體實施方式

現(xiàn)在結合附圖對本發(fā)明作進一步詳細的說明。

一、有限體積ronusthweno格式的構造

考慮二維euler方程：

其中，t表示時間變量，x，y表示空間變量，u＝(ρ，ρu，ρv，e)^t表示守恒變量，f(u)＝(ρu，ρu²+p，ρuv，u(e+p))^t，g(u)＝(ρv，ρuv，ρv²+p，v(e+p))^t，f(u)，g(u)表示通量，f(u)x表示f(u)對x求導，g(u)y表示g(u)對y求導，ρ，p，u，v，e分別表示流體密度、壓強、水平方向速度、豎直方向速度以及能量，t表示轉置，u0表示初始狀態(tài)值。

假設網(wǎng)格單元步長是固定的，xi+1/2-xi-1/2＝δx，yj+1/2-yj-1/2＝δy，網(wǎng)格中心為網(wǎng)格單元為iij＝[xi-1/2，xi+1/2]×[yj-1/2，yj+1/2]，其中下標i、j為坐標序號。

對方程(1)求導有：

其中，表示u對x求導，表示u對y求導，h(u，v)x和q(u，w)x表示h(u，v)q(u，w)分別對x求導，r(u，v)y和s(u，w)y表示r(u，v)s(u，w)分別對y求導，其中h(u，v)＝f′(u)v，r(u，v)＝g′(u)v，q(u，w)＝f′(u)w，s(u，w)＝g′(u)w，f′(u)，g′(u)分別表示f(u)，g(u)對u求導。

定義表示u在網(wǎng)格單元iij內的平均值，表示v在網(wǎng)格單元iij內的平均值，表示w在網(wǎng)格單元iij內的平均值，即

首先，求每個gauss點gk處u(gk，t)，v(gk，t)，w(gk，t)的值，具體步驟如下：

步驟1，如圖1所示，將目標單元iij以及周圍共9個網(wǎng)格單元組成一個母模板，且假設網(wǎng)格步長都為h，如圖1，其中i5為目標單元iij，并記分別代表u，v，w在網(wǎng)格單元im的平均值(m＝1，…，9)。

步驟2，在母模板中選擇8個子模板分別為：s1＝{i1，i2，i4，i5}，s2＝{i2，i3，i5，i6}，s3＝{i4，i5，i7，i8}，s4＝{i5，i6，i7，i9}，s5＝{i1，i2，i3，i4，i5，i7}，s6＝{i1，i2，i3，i5，i6，i9}，s7＝{i1，i4，i5，i7，i8，i9}，s8＝{i3，i5，i6，i7，i8，i9}，其中ii為對應序號的網(wǎng)格單元。

步驟3，在每個子模板中分別求出每個gauss點gk處u，v，w的重構多項式pn(x，y)，n＝1，…，8。

步驟3.1，u的重構多項式pn(x，y)，n＝1，…，8求解過程如下：

步驟3.1.1，在子模板s1，s2，s3，s4上構造多項式pn(x，y)，n＝1，…，4，使其滿足：

其中

n＝1，k＝1，2，4，5，kx＝4，ky＝2；n＝2，k＝2，3，5，6，kx＝6，ky＝2；

n＝3，k＝4，5，7，8，kx＝4，ky＝8；n＝4，k＝5，6，8，9，kx＝6，ky＝8；

步驟3.1.2，在子模板s5，s6，s7，s8上構造多項式pn(x，y)，n＝5，…，8，使其滿足：

其中

n＝5，k＝1，2，3，4，5，7；n＝6，k＝1，2，3，5，6，9；

n＝7，k＝1，4，5，7，8，9；n＝8，k＝3，5，6，7，8，9；

步驟3.1.3，得到每個子模板在每個gauss點gk處的插值多項式pn(x，y)，n＝1，…，8，如下：

步驟3.2，v的重構多項式pn(x，y)，n＝1，…，8求解過程如下：

步驟3.2.1，在子模板s1，s2，s3，s4上構造多項式pn(x，y)，n＝1，…，4，使其滿足：

n＝1，k＝1，2，4，5，kx＝1，4，5，ky＝1，2，5；n＝2，k＝2，3，5，6，kx＝3，5，6，ky＝2，3，5；

n＝3，k＝4，5，7，8，kx＝4，5，7，ky＝5，7，8；n＝4，k＝5，6，8，9，kx＝5，6，9，ky＝5，8，9；

步驟3.2.2，在子模板s5，s6，s7，s8上構造多項式pn(x，y)，n＝5，…，8，使其滿足：

其中

n＝5，kx＝1，2，3，4，5，7；n＝6，kx＝1，2，3，5，6，9；

n＝7，kx＝1，4，5，7，8，9；n＝8，kx＝3，5，6，7，8，9；

步驟3.2.3，得到每個子模板在每個gauss點gk處的插值多項式pn(x，y)，n＝1，…，8，如下：

其中，表示pn(x，y)對x求偏導(n＝1，…，4)，h為步長；

步驟3.3，w的重構多項式pn(x，y)，n＝1，…，8求解過程如下：

步驟3.3.1，在子模板s1，s2，s3，s4上構造多項式pn(x，y)，n＝1，…，4，使其滿足：

其中

n＝1，k＝1，2，4，5，kx＝1，4，5，ky＝1，2，5；n＝2，k＝2，3，5，6，kx＝3，5，6，ky＝2，3，5；

n＝3，k＝4，5，7，8，kx＝4，5，7，ky＝5，7，8；n＝4，k＝5，6，8，9，kx＝5，6，9，ky＝5，8，9；

步驟3.3.2，在子模板s5，s6，s7，s8上構造多項式pn(x，y)，n＝5，…，8，使其滿足：

其中

n＝5，ky＝1，2，3，4，5，7；n＝6，ky＝1，2，3，5，6，9；

n＝7，ky＝1，4，5，7，8，9；n＝8，ky＝3，5，6，7，8，9；

步驟3.3.3，得到每個子模板在每個gauss點gk處的插值多項式pn(x，y)，n＝1，…，8，如下：

其中，表示pn(x，y)對y求偏導(n＝1，…，4)，h為步長。

步驟4，取每個gauss點gk處的線性權為

步驟5，計算光滑指示器βl，l＝1，…，8，計算公式如下：

步驟5.1，重構u的光滑指示器計算公式為：

步驟5.2，重構v的光滑指示器計算公式為：

步驟5.3，重構w的光滑指示器計算公式為：

其中，a是多重指標，|i5|是目標單元的面積，d是求導算子，下標l是子模板對應的序號。

步驟6，通過線性權gl和光滑指示器βl計算非線性權ωl：

其中，為計算中的過渡值，βl為光滑指示器，ε取值10^-6，防止分母為零。

步驟7，進一步得到u，v，w在每個gauss點gk處的近似值：

其中，表示pl(gk)對x求偏導，表示pl(gk)對y求偏導，pl(gk)表示將gauss點的坐標gk代入相應多項式pl(x，y)后的取值；

其次，對方程(1)、(2)、(3)在目標單元iij內積分得到半離散的有限體積格式如下：

其中，表示u在網(wǎng)格單元iij內的平均值，表示v在網(wǎng)格單元iij內的平均值，表示w在網(wǎng)格單元iij內的平均值，即f＝(f(u)，g(u))^t，h＝(h(u，v)，r(u，v))^t，q＝(q(u，w)，s(u，w))^t，t表示轉置，為外法向向量，等式右邊的積分項可用q點gauss求積公式求得，如下：

其中，表示目標單元的邊界長度，ak為gauss求積公式的權重系數(shù)，為單位外法向向量，用數(shù)值通量代替；

最后，利用三階tvdrunge-kutta離散公式：

得到時空全離散有限體積格式，其中為中間過渡值，δt為時間步長，上標n表示第n時間層，為積分近似值。

二、虛擬單元方法

在實際計算當中，由于笛卡爾網(wǎng)格的非貼體性，有限體積robusthweno格式在處理復雜幾何外形的物面可壓繞流問題時非常困難，在網(wǎng)格和物面的相交處會產(chǎn)生非物理振蕩。為彌補這個不足，本發(fā)明采用浸入邊界法處理物面邊界。如圖2所示，實心三角形表示物面內部虛擬網(wǎng)格單元點，實心正方形表示虛擬網(wǎng)格單元點關于物面邊界的對稱點。

我們需要通過物面邊界條件以及動量關系等利用流體內b點的值確定物體內a點值，進而進行全流場的數(shù)值模擬。下面介紹兩種不同的方法確定a點流場值和b點流場值之間的關系并進行數(shù)值實驗，數(shù)值模擬出計算區(qū)域內包含復雜幾何外形物面的全流場。

2.1、st(symmetrytechnique)方法

這種方法是確定a點和b點流場值最簡單最直接的方法，計算公式如下：

上式中r為密度，p為壓強，為速度矢量，表示物面處的單位外法線向量，這樣給出的邊界條件可以滿足壁面無穿透條件。

2.2、fgcm(forrer’sghostcellmethod)方法

forrer等提出了一種虛擬單元浸入邊界法，其思想是a點的速度仍然使用一般的對稱反射，而密度和壓強通過以下公式求得：

其中，xw為過a點的法向向量與物面邊界相交的點，為物面的單位外法向量，同樣地，這樣給出的邊界條件也滿足壁面無穿透條件。

三、有限體積robusthweno格式與虛擬單元方法的結合

由以上格式的構造步驟可知，有限體積robusthweno格式是一種有限體積高精度數(shù)值計算格式。有限體積方法從積分守恒形式的流體方程出發(fā)，通過對網(wǎng)格單元上的方程進行離散來構造積分型離散方程。相比于有限差分方法，有限體積方法在網(wǎng)格處理上更加靈活，守恒性好。但是，對于在笛卡爾結構網(wǎng)格上構造的有限體積robusthweno格式，在求解雙馬赫反射、圓柱繞流等復雜界面繞流問題時是困難的。因為有限體積加權基本無振蕩格式在對網(wǎng)格單元流場值進行多項式重構時需要用到周圍(包括自身)一共9個網(wǎng)格單元，尤其在網(wǎng)格與物面邊界相交的地方，網(wǎng)格與物面以任意形式相交，問題就顯得尤為突出。面對這些問題，在雙馬赫反射數(shù)值模擬中，傳統(tǒng)的解決辦法是使入射激波傾斜，斜面便與網(wǎng)格線重合，數(shù)值模擬就變得清晰而簡單了。這樣雖然能較好地捕捉到激波的反射情況，但和實驗模擬無論是在形式上還是本質上都是有不同的。而在圓柱繞流問題或者是翼型繞流問題中，傳統(tǒng)的解決辦法是采用混合網(wǎng)格或者非結構網(wǎng)格，這都無疑增大了網(wǎng)格生成的難度以及對應網(wǎng)格上格式構造的難度。本發(fā)明將結構網(wǎng)格有限體積robusthweno格式同浸入邊界虛擬單元方法結合，有效模擬了具有真實斜面的雙馬赫反射問題，精確地捕捉了激波與斜面作用后反射的情況。同樣地，在高馬赫數(shù)圓柱繞流和翼型中也得到了較好的數(shù)值模擬結果。

下面給出四個算例作為本發(fā)明所公開方法的具體實施例。

實施例一、臺階問題。該問題是emery于1968年提出的一個用于檢驗非線性雙曲型守恒律格式的經(jīng)典算例。初始數(shù)據(jù)為水平來流馬赫數(shù)為3，密度為1.4，壓強為1，管道區(qū)域為[0，3]′[0，1]，在距離左邊界0.6處有一高度為0.2的臺階，且臺階延伸到管道的盡頭。上下邊界為反射邊界，左邊界為來流邊界，右邊界為出流邊界。圖3a-3b和圖4a-4b分別給出了t＝4時刻的密度和壓強的等值線圖。

實施例二、斜面雙馬赫反射問題。該問題是描述一個強激波入射在與平面成30°角的斜坡上后發(fā)生的變化，來流是馬赫數(shù)為10的強激波。計算區(qū)域取[0，3]′[0，2]，始于激波垂直于x軸，初始數(shù)據(jù)為：

其中的左右狀態(tài)為ul＝(8，66.009，0，563.544)^t，ur＝(1.4，0，0，2.5)^t，邊界條件處理為左右邊界分別取左右狀態(tài)值。下邊界：當時，設為反射邊界；當時，設為左狀態(tài)。上邊界：當x＞g(t)時，取右狀態(tài)；當時，取左狀態(tài)，其中圖5a-5b和圖6a-6b分別給出了t＝0.2時刻密度和壓強的等值線圖。

實施例三、圓柱繞流問題?？紤]超音速無黏流體流向直徑為1的圓柱，初始數(shù)據(jù)為水平來流馬赫數(shù)為3，密度為1，壓強為1，計算區(qū)域為[-10，10]？[10，10]，左邊界為來流邊界，右邊界為超音速出流邊界。圖7a-7b和圖8a-8b分別給出了在計算區(qū)域內兩種方法得到的密度和馬赫數(shù)等值線。

實施例四、naca0012翼型跨聲速繞流問題。來流馬赫數(shù)為0.8，攻角a＝1.25度。圖9a-9b和圖10a-10b分別給出了st方法和fgcm方法壓強等值線圖和翼型表面壓力系數(shù)圖以及殘差圖，從圖中可以看出有限體積robusthweno格式結合兩種浸入邊界方法能夠很好的捕捉到激波的位置以及激波的強弱，上表面形成一道強激波，下表面形成一道弱激波。

以上僅是本發(fā)明的優(yōu)選實施方式，本發(fā)明的保護范圍并不僅局限于上述實施例，凡屬于本發(fā)明思路下的技術方案均屬于本發(fā)明的保護范圍。應當指出，對于本技術領域的普通技術人員來說，在不脫離本發(fā)明原理前提下的若干改進和潤飾，應視為本發(fā)明的保護范圍。