一種基于三階WENO重構(gòu)的無(wú)網(wǎng)格算法研究

蒲賽虎，陳紅全

（南京航空航天大學(xué) 航空宇航學(xué)院，江蘇 南京 210016）

0 引 言

無(wú)網(wǎng)格算法是繼網(wǎng)格算法之后出現(xiàn)的一種新型數(shù)值計(jì)算方法。在剖分計(jì)算區(qū)域時(shí)，由于只涉及布點(diǎn)填充，既可以利用傳統(tǒng)網(wǎng)格算法采用的網(wǎng)格點(diǎn)，也可以在需要的前提下直接布點(diǎn)，在處理復(fù)雜外形方面，更具有靈活性［1－3］。

Batina［1］自20世紀(jì)90年代初就開始無(wú)網(wǎng)格算法的應(yīng)用研究，提出用點(diǎn)云離散計(jì)算區(qū)域，代替通常的網(wǎng)格劃分，并在當(dāng)?shù)攸c(diǎn)云上，利用最小二乘法逼近計(jì)算空間導(dǎo)數(shù)，由此發(fā)展出基于Jameson中心格式的無(wú)網(wǎng)格算法。但中心格式一般數(shù)值耗散較大，而且為了抑制高頻振蕩等還必須加入人工粘性項(xiàng)。為了克服上述不足，Morinishi通過(guò)在點(diǎn)云的衛(wèi)星點(diǎn)連線中點(diǎn)處引入交界面，將一種逆風(fēng)型的Roe格式引入到無(wú)網(wǎng)格算法中，減小了數(shù)值耗散，同時(shí)通過(guò)線性逼近重構(gòu)交界面左右狀態(tài)值，使得算法能達(dá)到二階精度［2］。之后，許多學(xué)者針對(duì)具體的求解問(wèn)題，分別將不同的逆風(fēng)型格式引入到無(wú)網(wǎng)格算法中［4－6］，但整體上對(duì)于算法涉及的交界面左右狀態(tài)值的確定，大多是基于線性逼近重構(gòu)，因此其計(jì)算精度至多能達(dá)到二階。用高階的界面逼近重構(gòu)替代傳統(tǒng)的無(wú)網(wǎng)格線性逼近重構(gòu)，應(yīng)有助于提高無(wú)網(wǎng)格算法的求解精度。

本文考慮在無(wú)網(wǎng)格算法中，引入發(fā)展相對(duì)成熟的WENO（Weighted Essentially Non－Oscillatory）重構(gòu)。該重構(gòu)法最早是由Liu等人［7］在ENO重構(gòu)［8］的基礎(chǔ)上發(fā)展而來(lái)，它將ENO重構(gòu)選擇最光滑模板進(jìn)行數(shù)值逼近的處理方法，改進(jìn)為對(duì)所有模板的數(shù)值逼近進(jìn)行加權(quán)求和，通過(guò)構(gòu)造合適的權(quán)系數(shù)，在光滑區(qū)域可以達(dá)到比ENO重構(gòu)更高的精度，且具有更好的收斂性和更好的穩(wěn)健性，而在間斷附近，卻保持有ENO重構(gòu)所具有的基本無(wú)振蕩的特性。早在1996年，Jiang G.S.和Shu C.W在Liu等人研究的基礎(chǔ)上提出了三階和五階有限差分 WENO重構(gòu)［9］。本文將致力于把該三階有限差分WENO重構(gòu)引入到無(wú)網(wǎng)格算法中，發(fā)展出一種基于三階WENO重構(gòu)的無(wú)網(wǎng)格算法?；跓o(wú)網(wǎng)格點(diǎn)云結(jié)構(gòu)，先通過(guò)在點(diǎn)云的每個(gè)衛(wèi)星點(diǎn)方向上引入局部一維坐標(biāo)系，并結(jié)合虛擬點(diǎn)的設(shè)置，將有限差分WENO重構(gòu)中基于一維坐標(biāo)系的模板構(gòu)造方法，發(fā)展用于無(wú)網(wǎng)格點(diǎn)云中的模板構(gòu)造；對(duì)于模板上虛擬點(diǎn)的流場(chǎng)值，則采用一種基于最近節(jié)點(diǎn)流場(chǎng)值的插值方法確定。該虛擬點(diǎn)插值方法可利用已有的點(diǎn)云信息，因此避免了直接尋找插值點(diǎn)等耗時(shí)的步驟，簡(jiǎn)化了插值操作；在構(gòu)造的模板上，采用三階WENO重構(gòu)確定算法所涉及的交界面左右狀態(tài)值；基于上述WENO重構(gòu)方法，并結(jié)合Roe的近似Riemann求解器求得數(shù)值通量，對(duì)Euler方程進(jìn)行了求解，編程計(jì)算了典型流動(dòng)算例，驗(yàn)證了所發(fā)展的算法獲得的數(shù)值解能逼近三階精度。在此基礎(chǔ)上，給出了若干繞流算例，展示了所提算法捕捉激波和非定常激波演化等復(fù)雜流動(dòng)問(wèn)題的能力。

1 無(wú)網(wǎng)格空間導(dǎo)數(shù)逼近簡(jiǎn)介

采用無(wú)網(wǎng)格算法，流場(chǎng)中不需要網(wǎng)格，只需要有限數(shù)量的節(jié)點(diǎn)。每一個(gè)節(jié)點(diǎn)與它周圍的節(jié)點(diǎn)形成點(diǎn)云。圖1給出了一個(gè)七點(diǎn)無(wú)網(wǎng)格點(diǎn)云。其中i點(diǎn)是中心點(diǎn)，點(diǎn)1～點(diǎn)6是其衛(wèi)星點(diǎn)。

以函數(shù)f為例，在i點(diǎn)的每一個(gè)衛(wèi)星點(diǎn)上，其函數(shù)值fk可以用i點(diǎn)處的函數(shù)值通過(guò)泰勒級(jí)數(shù)展開得到：

其 中 Δxk＝xk－xi，Δyk＝y(tǒng)k－yi，aj（j＝1，…，5）表示i點(diǎn)處函數(shù)f的各階偏導(dǎo)數(shù)，即：

圖1 無(wú)網(wǎng)格點(diǎn)云Fig.1 Cloud of points for gridless method

將式（1）保留到二階項(xiàng)，則得到fk的近似值為：

則aj（j＝1，…，5） 的值可以通過(guò)求解如下的最小二乘問(wèn)題確定：

其中M表示i點(diǎn)的衛(wèi)星點(diǎn)數(shù)。通過(guò)式（4）確定aj（j＝1，…，5）的具體過(guò)程建議閱讀文獻(xiàn)［10］的詳細(xì)描述，這里給出最終的求解方程組：

其中∑是 的簡(jiǎn)寫形式。方程組（5）的解可以寫成點(diǎn)云中各點(diǎn)函數(shù)值的線性組合，即：

系數(shù)αk，βk，γk，λk，ωk由方程組（5）確定。上述各階偏導(dǎo)數(shù)也可以采用中心點(diǎn)與衛(wèi)星點(diǎn)連線中點(diǎn)處的函數(shù)值fik類似求得［10］：

式（7）中的系數(shù)αik，βik，γik，λik，ωik與式（6）中的系數(shù)αk，βk，γk，λk，ωk有如下關(guān)系：

可以看到，這些系數(shù)只與點(diǎn)云中各點(diǎn)的幾何位置有關(guān)，因此可在流場(chǎng)迭代計(jì)算前計(jì)算儲(chǔ)存好。

2 控制方程

本文研究的三階WENO重構(gòu)及無(wú)網(wǎng)格算法將基于Euler方程展開。在直角坐標(biāo)系下，守恒形式的Euler方程可寫為：

式中，W是守恒矢變量，E和F是對(duì)流矢通量，可分別寫為：

其中，ρ為密度，u、v為沿坐標(biāo)x、y上的速度分量，p為氣體的壓強(qiáng)，e是單位體積內(nèi)總能。對(duì)完全氣體滿足狀態(tài)方程

式中，空氣的比熱比γ＝1.4。

3 基于三階WENO重構(gòu)的無(wú)網(wǎng)格算法

應(yīng)用式（7），Euler方程的通量項(xiàng)可近似寫為：

Gik是i點(diǎn)和k點(diǎn)連線中點(diǎn)處的數(shù)值通量，本文按Roe的近似 Riemann解確定［11］：

其中A＝是Jacobian矩陣，和是i點(diǎn)和k點(diǎn)連線中點(diǎn)處左右兩側(cè)的守恒變量值（見圖1中所示）。若采用如下傳統(tǒng)的線性逼近重構(gòu)確定和，則空間精度能達(dá)到二階［2，4］：

其中rik是從i點(diǎn)指向k點(diǎn)的矢量，φ－和φ＋為通量限制器［2］，守恒變量的梯度 ?Wi在每個(gè)點(diǎn)上用式（6）計(jì)算給出。

如引言中所述，我們希望空間精度能進(jìn)一步提高。為此，本文將三階WENO重構(gòu)引入到無(wú)網(wǎng)格算法中，以替代式（14）的線性逼近重構(gòu)。下面介紹無(wú)網(wǎng)格點(diǎn)云上三階WENO重構(gòu)的具體實(shí)現(xiàn)方法。

3.1 無(wú)網(wǎng)格點(diǎn)云上三階WENO重構(gòu)的實(shí)現(xiàn)

本文無(wú)網(wǎng)格點(diǎn)云上的三階WENO重構(gòu)，是基于有限差分 WENO重構(gòu)發(fā)展而來(lái)，因此為了敘述方便，我們首先對(duì)三階有限差分WENO重構(gòu)做一簡(jiǎn)單的介紹，然后給出在無(wú)網(wǎng)格點(diǎn)云上實(shí)施三階WENO重構(gòu)的過(guò)程。

三階有限差分WENO重構(gòu)是基于一維坐標(biāo)系進(jìn)行的［9］，如圖2（a）所示。以函數(shù)f為例，每個(gè)交界面上左右狀態(tài)值的重構(gòu)分別基于由3個(gè)節(jié)點(diǎn)構(gòu)成的模板。以i＋1／2處左狀態(tài)值的重構(gòu)為例，其重構(gòu)模板為S＝｛i－1，i，i＋1｝（圖2（a）中的紅色點(diǎn)），將該模板分成兩個(gè)子模板S1＝｛i－1，i｝、S2＝｛i，i＋1｝。在這兩個(gè)子模板上，根據(jù)線性分布規(guī)律可分別構(gòu)造i＋1／2處的狀態(tài)值為：

根據(jù) WENO重構(gòu)思想［9］，i＋1／2處的左狀態(tài)值為和的加權(quán)和，即：

非線性權(quán)系數(shù)w1、w2定義為：

其中＝、c＝為優(yōu)化權(quán)系數(shù)，ε是為了防止在2光滑流動(dòng)區(qū)域分母為零而加入的一個(gè)很小的數(shù)，本文取為10－5。Si是模板的光滑度量系數(shù)，定義為：

圖2 無(wú)網(wǎng)格點(diǎn)云上三階WENO重構(gòu)模板的構(gòu)造Fig.2 The stencil of third order WENO reconstruction in a cloud of points

為了將上述有限差分WENO重構(gòu)發(fā)展用于無(wú)網(wǎng)格算法，首先需要確定無(wú)網(wǎng)格點(diǎn)云上的重構(gòu)模板。為此，本文基于無(wú)網(wǎng)格點(diǎn)云衛(wèi)星點(diǎn)分布特征，沿點(diǎn)云中心點(diǎn)與衛(wèi)星點(diǎn)連線方向引入了局部一維坐標(biāo)系，圖2（b）給出了沿k點(diǎn)方向引入的局部一維坐標(biāo)系。在此局部一維坐標(biāo)系中，我們?cè)趇點(diǎn)的上游方向設(shè)置一虛擬點(diǎn)j，并且j點(diǎn)和k點(diǎn)關(guān)于i點(diǎn)對(duì)稱（關(guān)于虛擬點(diǎn)上流場(chǎng)值的確定方法將在下一節(jié)中介紹）?；诖司植恳痪S坐標(biāo)系，我們參照有限差分WENO重構(gòu)方法，采用j點(diǎn)、i點(diǎn)和k點(diǎn)上的值，通過(guò)式（16）構(gòu)造i點(diǎn)和k點(diǎn)連線中點(diǎn)處的左狀態(tài)值，即：

非線性權(quán)系數(shù)w1、w2也按式（17）計(jì)算，此時(shí)模板的光滑度量系數(shù)相應(yīng)變?yōu)椋?/p>

根據(jù)對(duì)稱性，i點(diǎn)和k點(diǎn)連線中點(diǎn)處右狀態(tài)值可以類似得到：

非線性權(quán)系數(shù)w1、w2也按式（17）計(jì)算，此時(shí)模板的光滑度量系數(shù)相應(yīng)變?yōu)椋?/p>

式（21）和式（22）中的fl表示上述局部一維坐標(biāo)系中，i點(diǎn)關(guān)于k點(diǎn)的對(duì)稱點(diǎn)l上的函數(shù)值（圖2b中沒有畫出）。

將式（19）和式（21）應(yīng)用于守恒變量的重構(gòu)，就可得到中心點(diǎn)和衛(wèi)星點(diǎn)連線中點(diǎn)處守恒變量的左右狀態(tài)值。再將左右狀態(tài)值帶入式（13），就可得到每個(gè)衛(wèi)星點(diǎn)方向上的數(shù)值通量。

3.2 虛擬點(diǎn)上流場(chǎng)值的確定

在上述WENO重構(gòu)模板構(gòu)造過(guò)程中我們?cè)O(shè)置了虛擬點(diǎn)，這些虛擬點(diǎn)上的流場(chǎng)值可以通過(guò)插值確定。一種最直接的做法是搜索到虛擬點(diǎn)周圍一定范圍內(nèi)的所有節(jié)點(diǎn)，作為插值點(diǎn)，然后根據(jù)插值方法確定虛擬點(diǎn)上的值。這種做法涉及到逐一搜索插值點(diǎn)的過(guò)程，可能需要花費(fèi)較多的機(jī)時(shí)，而隨后進(jìn)行的插值系數(shù)確定過(guò)程，也可能涉及矩陣求逆等操作，同樣較為耗時(shí)?？紤]到對(duì)本文采用的無(wú)網(wǎng)格算法而言，每個(gè)節(jié)點(diǎn)的點(diǎn)云結(jié)構(gòu)及求導(dǎo)系數(shù)（即式（6）中的系數(shù)）本來(lái)就已經(jīng)確定，我們?nèi)舫浞掷眠@些已有信息，則有望簡(jiǎn)化插值處理過(guò)程?；诖?，本文提出了一種基于最近點(diǎn)流場(chǎng)值的插值處理方法。以圖2（b）中虛擬點(diǎn)j的流場(chǎng)值確定為例，該插值方法的實(shí)施過(guò)程如下：

（1）把引入虛擬點(diǎn)j的點(diǎn)云中的節(jié)點(diǎn)作為候選點(diǎn)（即圖2b中的i點(diǎn)及其衛(wèi)星點(diǎn)），通過(guò)比較與j點(diǎn)的距離，找到這些候選點(diǎn)中距離j點(diǎn)最近的點(diǎn)，為圖2（b）中的p點(diǎn)；

（2）虛擬點(diǎn)j上的流場(chǎng)值則可基于p點(diǎn)的流場(chǎng)值，按如下方式逼近：

其中Δxjp＝xj－xp，Δyjp＝y(tǒng)j－yp。將各階偏導(dǎo)數(shù)的表達(dá)式（6）帶入，則式（23）可以進(jìn)一步寫為：

將含有相同函數(shù)值的項(xiàng)進(jìn)行合并，式（24）最終可整理成如下簡(jiǎn)單形式：

其中Mp表示p點(diǎn)的衛(wèi)星點(diǎn)總數(shù)，各插值系數(shù)具體形式為：

可以看到，上述插值方法實(shí)質(zhì)上是直接取p點(diǎn)點(diǎn)云中的節(jié)點(diǎn)作為虛擬點(diǎn)j的插值點(diǎn)，因此避免了逐一搜索插值點(diǎn)的過(guò)程，節(jié)約了計(jì)算時(shí)間。另外，插值系數(shù)也是基于點(diǎn)云已有的求導(dǎo)系數(shù)，直接采用式（26）通過(guò)簡(jiǎn)單的代數(shù)運(yùn)算得到，而不需額外引入插值系數(shù)，編程實(shí)現(xiàn)簡(jiǎn)單?？梢钥吹?，上述插值系數(shù)也只與幾何量有關(guān)，因此同樣可以在流場(chǎng)迭代計(jì)算前計(jì)算存儲(chǔ)好。

3.3 時(shí)間離散格式

采用上述基于WENO重構(gòu)的無(wú)網(wǎng)格算法離散Euler方程的空間導(dǎo)數(shù)后，可以得到方程的半離散形式如下：

其中，Ri為的離散形式。

為了達(dá)到空間和時(shí)間的一致高階精度，本文采用Shu C W 和Osher提出的三階TVD Runge－Kutta時(shí)間離散格式［12］：

其中Δt為設(shè)定的時(shí)間步長(zhǎng)。具體求解時(shí)還涉及到邊界條件的處理，本文對(duì)各種邊界條件沿用了有限差分WENO重構(gòu)邊界條件的處理方法（如對(duì)周期性邊界條件通過(guò)設(shè)置輔助點(diǎn)的方式處理），限于篇幅，不再詳細(xì)描述，這里列出文獻(xiàn)［13］供參考。

4 算例與分析

本文已采用上述基于三階WENO重構(gòu)的無(wú)網(wǎng)格算法進(jìn)行了編程，并對(duì)典型流動(dòng)問(wèn)題進(jìn)行了計(jì)算。這里首先給出求解模型問(wèn)題典型流動(dòng)的計(jì)算結(jié)果，以驗(yàn)證所提算法的精度，在此基礎(chǔ)上對(duì)存在間斷的流動(dòng)問(wèn)題進(jìn)行了數(shù)值模擬。

4.1 二維線性模型方程的求解算例

二維線性模型方程的精確解很容易求得，因此常用于檢測(cè)算法的精度，方程形式及初值如下：

利用特征線法，可求得精確解為u（t，x，y） ＝sin（2π（x＋y－2t））。本文計(jì)算區(qū)域取為［0，1］×［0，1］，邊界條件為周期性邊界條件。我們分別采用規(guī)則分布的節(jié)點(diǎn)（如圖3（a））和不規(guī)則分布的節(jié)點(diǎn)（如圖3b）進(jìn)行了計(jì)算。表1為t＝0.2時(shí)刻，采用不同密度的規(guī)則分布的節(jié)點(diǎn)進(jìn)行的精度分析（表中L1表示所有節(jié)點(diǎn)上誤差的平均值，L∞表示所有節(jié)點(diǎn)上誤差的最大值），表2為同一時(shí)刻，采用不同密度的不規(guī)則分布的節(jié)點(diǎn)進(jìn)行的精度分析。從表1和表2可以看出，采用上述兩種布點(diǎn)形式，所發(fā)展的無(wú)網(wǎng)格算法總體都能逼近三階精度。圖4還給出了采用41×41個(gè)規(guī)則分布的節(jié)點(diǎn)，經(jīng)過(guò)較長(zhǎng)時(shí)間后t＝2時(shí)刻），沿y＝0.5解的分布，作為比較，還同時(shí)給出了基于線性逼近重構(gòu)的計(jì)算結(jié)果?？梢钥吹剑谌AWENO重構(gòu)的計(jì)算結(jié)果與精確解符合更好，在極值點(diǎn)附近也沒有觀察到明顯耗散。

圖3 兩種不同的布點(diǎn)形式Fig.3 Two different types of point distribution

表1 二維線性模型方程的精度分析（采用規(guī)則分布的點(diǎn)）Table 1 Accuracy of 2Dlinear equation（regular point distribution）

表2 二維線性模型方程的精度分析（采用不規(guī)則分布的點(diǎn)）Table 2 Accuracy of 2Dlinear equation（irregular point distribution）

圖4 基于線性逼近重構(gòu)與基于WENO重構(gòu)得到的y＝0.5位置解的分布（t＝2）Fig.4 Comparison of results using linear reconstruction and WENO reconstruction along y＝0.5at t＝2

4.2 二維Euler方程的求解算例

4.2.1 等熵渦問(wèn)題

本節(jié)采用所發(fā)展的無(wú)網(wǎng)格算法求解二維等熵渦問(wèn)題，以檢測(cè)算法求解Euler方程的精度。設(shè)初始的平均流為 （ρ，u，v，p）＝ （1 ，1 ，1，1），在平均流 上 加入一個(gè)等熵渦：

其中r2＝（x－5）2＋（y－5）2，ε＝5，該問(wèn)題的精確解是等熵渦隨平均流等速移動(dòng)。本文計(jì)算區(qū)域?。?，10］×［0，10］，邊界條件取周期性邊界條件。與前一算例一樣，我們分別采用規(guī)則分布的節(jié)點(diǎn)和不規(guī)則分布的節(jié)點(diǎn)進(jìn)行了計(jì)算（這兩種節(jié)點(diǎn)分布形式類似于圖3，限于篇幅，這里不再展示）。表3和表4分別給出了t＝2時(shí)刻，采用不同密度的規(guī)則分布的節(jié)點(diǎn)和不規(guī)則分布的節(jié)點(diǎn)進(jìn)行的精度分析?？梢钥闯?，無(wú)論是采用規(guī)則分布的節(jié)點(diǎn)還是不規(guī)則分布的節(jié)點(diǎn)，所發(fā)展的無(wú)網(wǎng)格算法總體都能逼近三階精度。圖5還給出了采用81×81個(gè)規(guī)則分布的節(jié)點(diǎn)計(jì)算出的t＝20時(shí)刻沿y＝5流體密度的分布，作為比較，還同時(shí)給出了基于線性逼近重構(gòu)的計(jì)算結(jié)果?？梢钥吹?，基于三階WENO重構(gòu)得到的結(jié)果與精確值符合更好，特別是在極值點(diǎn)附近也沒有觀察到明顯的耗散。

表3 二維Euler方程的精度分析（采用規(guī)則分布的點(diǎn)）Table 3 Accuracy of 2DEuler equations（regular point distribution）

表4 二維Euler方程的精度分析（采用不規(guī)則分布的點(diǎn)）Table 4 Accuracy of 2DEuler equations（irregular point distribution）

圖5 基于線性逼近重構(gòu)與基于WENO重構(gòu)得到的y＝5位置的密度分布（t＝20）Fig.5 Comparison of results using linear reconstruction and WENO reconstruction along y＝5at t＝20

4.2.2 激波管流動(dòng)問(wèn)題

激波管流動(dòng)問(wèn)題常用于檢測(cè)算法對(duì)激波，接觸間斷等的捕捉能力。本算例計(jì)算區(qū)域?yàn)椋? ，1］×［0 ，0 .1］的矩形區(qū)域，其中布置了76090個(gè)不規(guī)則分布的節(jié)點(diǎn)。初始條件設(shè)置為［14］：

本文分別采用基于線性逼近重構(gòu)和三階WENO重構(gòu)的無(wú)網(wǎng)格算法進(jìn)行了計(jì)算，算至?xí)r刻t＝0.0039。圖6為流體密度在激波和接觸間斷附近的計(jì)算結(jié)果，圖7則為對(duì)應(yīng)的在膨脹波附近的計(jì)算結(jié)果。整體上，基于三階WENO重構(gòu)的無(wú)網(wǎng)格算法捕捉的激波、接觸間斷和膨脹波較線性逼近重構(gòu)更靠近精確值，與文獻(xiàn)［14］五階有限差分 WENO重構(gòu)（x方向點(diǎn)的間距與本文的相同）結(jié)果接近（見圖6、圖7），且在膨脹波附近也沒有明顯的振蕩，體現(xiàn)出WENO重構(gòu) “基本無(wú)振蕩”的特性（見圖7）。

圖6 激波和接觸間斷附近的密度分布Fig.6 Density distribution near shock and contact discontinuity

圖7 膨脹波附近的密度分布Fig.7 Density distribution near expansion wave

4.2.3 超聲速半圓柱繞流

考慮馬赫數(shù)為3的超聲速半圓柱繞流，圓柱半徑為r＝0.5?；趫D8的計(jì)算區(qū)域，布置了101×61個(gè)規(guī)則分布的節(jié)點(diǎn)。外邊界采用自由來(lái)流條件，物面采用無(wú)穿透條件，上下出流邊界采用外推法處理［15］，圖8分別給出了無(wú)網(wǎng)格算法基于線性逼近重構(gòu)以及WENO重構(gòu)的壓強(qiáng)等值線，圖中可看出，兩種方法都捕捉到了脫體弓形激波，激波距駐點(diǎn)的距離都為0.68r（這與文獻(xiàn)［15］計(jì)算出的距離吻合），但基于WENO重構(gòu)的無(wú)網(wǎng)格算法捕捉的激波更加清晰，表明后者三階WENO重構(gòu)具有更高的分辨率。圖9給出了無(wú)網(wǎng)格算法基于線性逼近重構(gòu)以及WENO重構(gòu)的收斂曲線，整體看來(lái)，兩種方法的收斂曲線比較接近。另外，我們還對(duì)計(jì)算時(shí)間進(jìn)行了統(tǒng)計(jì)，兩種方法迭代5000步所需的時(shí)間分別為197s和203s（CPU為Intel Core I7－860）。這一結(jié)果展示出基于WENO重構(gòu)的無(wú)網(wǎng)格算法在計(jì)算效率上能與傳統(tǒng)基于線性逼近重構(gòu)的無(wú)網(wǎng)格算法相當(dāng)。

4.2.4 激波過(guò)彎道繞雙圓柱流場(chǎng)

在前面算例對(duì)發(fā)展的基于WENO重構(gòu)的無(wú)網(wǎng)格算法進(jìn)行驗(yàn)證的基礎(chǔ)上，本節(jié)將該算法應(yīng)用于求解激波過(guò)彎道繞雙圓柱流場(chǎng)，以展示其處理含非定常激波演化的復(fù)雜流動(dòng)問(wèn)題的效果。該算例幾何外形取自文［16］的實(shí)驗(yàn)?zāi)Ｐ停疚挠?jì)算捕捉到與實(shí)驗(yàn)類似的波系演化結(jié)構(gòu)，這里只給出計(jì)算結(jié)果。

圖8 計(jì)算得到的壓強(qiáng)等值線Fig.8 Comparison of the pressure contours obtained

圖9 收斂曲線Fig.9 Convergence history

圖10為本算例的計(jì)算區(qū)域及采用的節(jié)點(diǎn)分布（節(jié)點(diǎn)總數(shù)為31972個(gè)）。初始時(shí)刻，在彎道上方入口處有一馬赫數(shù)為1.8的激波，隨著流動(dòng)的發(fā)展，該激波逐漸沿彎道移動(dòng)，并與雙圓柱發(fā)生碰撞。圖11（a～f）給出了六個(gè)典型時(shí)刻的密度等值線。可以看到，當(dāng)與第一個(gè)圓柱相遇時(shí)，激波從物面上反射，進(jìn)一步會(huì)出現(xiàn)激波與激波的相互作用，出現(xiàn)馬赫桿、接觸間斷等流動(dòng)現(xiàn)象。激波繼續(xù)向下游移動(dòng)，在遇到第二個(gè)圓柱時(shí)又發(fā)生反射，同時(shí)，前后圓柱間的激波發(fā)生多次相撞，使得流場(chǎng)變得非常復(fù)雜，但本文計(jì)算出的結(jié)果給出了清晰的流場(chǎng)結(jié)構(gòu)，展示了所發(fā)展的算法處理含非定常激波演化流動(dòng)的能力，有望發(fā)展用于實(shí)際的復(fù)雜流動(dòng)數(shù)值模擬。

圖10 激波過(guò)彎道繞雙圓柱流場(chǎng)計(jì)算區(qū)域及節(jié)點(diǎn)分布Fig.10 Point distribution for simulating the flow field of shock wave through curved channel around double cylinders

圖11 六個(gè)典型時(shí)刻的密度等值線Fig.11 Density contours of the flow field at six typical times

5 結(jié) 論

本文通過(guò)在無(wú)網(wǎng)格算法中引入WENO重構(gòu)，發(fā)展了基于三階WENO重構(gòu)的無(wú)網(wǎng)格算法。算例驗(yàn)證了該算法獲得的數(shù)值解如預(yù)期能逼近三階精度；與傳統(tǒng)的基于線性逼近重構(gòu)的無(wú)網(wǎng)格算法相比，在相同的布點(diǎn)情況下，本文算法捕捉的激波等間斷問(wèn)題具有更高分辨率，能體現(xiàn)出WENO重構(gòu)“基本無(wú)振蕩”的特性。由于算法只需無(wú)網(wǎng)格布點(diǎn)，不但適合處理激波管等簡(jiǎn)單外形流動(dòng)問(wèn)題，也適合處理存在多個(gè)物體相互干擾的復(fù)雜流動(dòng)問(wèn)題。

［1］BATINA J T.A gridless Euler／Navier－Stokes solution algorithm for complex aircraft applications［R］.AIAA 93－0333.

［2］MORINISHI K.Gridless type solution for high Reynolds number multielement flow fields［R］.AIAA 95－1856.

［3］CHEN H Q.An implicit gridless method and its applications［J］.ACTAAerodynamicaSinica，2002，20（2）：133－140.（in Chinese）陳紅全.隱式無(wú)網(wǎng)格算法及其應(yīng)用研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2002，20（2）：133－140.

［4］SRIDAR D，BALAKRISHNAN N.An upwind finite difference scheme for meshless solvers［J］.JournalofComputational Physics，2003，189：1－29.

［5］SHENG M J，YE Z Y，JIANG C Q.Study of the boundary layer solution coupled with gridless method［J］.ChineseJournal ofComputationalMechanics，2011，28（6）：290－925.（in Chinese）盛鳴劍，葉正寅，蔣超奇.附面層修正應(yīng)用于無(wú)網(wǎng)格算法的研究［J］.計(jì)算力學(xué)學(xué)報(bào)，2011，28（6）：290－925.

［6］MA Z H.Research of adaptive meshfree and hybridized mesh／meshfree methods［D］.Nanjing：Nanjing University of Aeronautics and Astronautics，2008.（in Chinese）馬志華.自適應(yīng)無(wú)網(wǎng)格及網(wǎng)格和無(wú)網(wǎng)格混合算法研究［D］.南京：南京航空航天大學(xué)，2008.

［7］LU X D，OSHER S，CHAN T.Weighted essentially non－oscillatory schemes［J］.JournalofComputationalPhysics，1994，115（1）：200－212.

［8］HARTEN A，ENGQUIST B，OSHER S，et al.Uniformly high－order accurate essentially non－oscillatory shock－capturing schemes［J］.JournalofComputationalPhysics，1987，71（2）：231－303.

［9］JIANG G S，SHU C W.Efficient implementation of weighted ENO schemes［J］.JournalofComputationalPhysics，1996，126（1）：202－228.

［10］MA Z H，CHEN H Q，ZHOU C H.A study of point moving adaptivity in gridless method［J］.Computermethodsinapplied mechanicsandengineering，2008，197：1926－1937.

［11］ROE P L.Approximate Riemann solvers，parameter vectors，and difference schemes［J］.JournalofComputationalPhysics，1981，23：357－372.

［12］SHU C W，STANLEY OSHER.Efficient implementation of essentially non－oscillatory shock capturing schemes［J］.Journal ofComputationalPhysics，1988，77（2）：439－471.

［13］SHU C W.Essentially non－oscillatory and weighted essentially non－oscillatory schemes for hyperbolic conservation laws［R］.NACA／CR－97－206253，1－78.

［14］YANG S P，LI S F，QU X M，et al.Numerical comparison of fifth－order WENO schemes and second－order Godunov schemes［J］.NaturalScienceJournalofXiangtanUniversity，2006，28（2）：8－12.（in Chinese）楊水平，李壽佛，屈小妹，等.五階FD－WENO格式和二階Godunov格式MUSCL的數(shù)值測(cè)試與定量比較［J］.湘潭大學(xué)自然科學(xué)學(xué)報(bào)，2006，28（2）：8－12.

［15］SJOGREEN B.High order centered difference methods for the compressible Navier－Stokes equations［J］.JournalofComputationalPhysics，1995，117（1）：67－78.

［16］YE Xichao，LE Jialing，YANG Hui.Flow field of shock wave through curved channel around double cylinder［J］.AerodynamicExperimentAndMeasurement＆Control，1996，10（2）：80－85.