蒲賽虎，陳紅全

（南京航空航天大學 航空宇航學院，江蘇 南京 210016）

0 引 言

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

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

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

1 無網(wǎng)格空間導數(shù)逼近簡介

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

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

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

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

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

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

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

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

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

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

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

2 控制方程

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

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

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

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

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

應用式（7），Euler方程的通量項可近似寫為：

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

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

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

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

3.1 無網(wǎng)格點云上三階WENO重構的實現(xiàn)

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

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

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

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

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

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

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

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

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

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

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

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

3.2 虛擬點上流場值的確定

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

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

（2）虛擬點j上的流場值則可基于p點的流場值，按如下方式逼近：

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

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

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

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

3.3 時間離散格式

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

其中，Ri為的離散形式。

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

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

4 算例與分析

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

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

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

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

圖3 兩種不同的布點形式Fig.3 Two different types of point distribution

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

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

圖4 基于線性逼近重構與基于WENO重構得到的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 等熵渦問題

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

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

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

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

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

4.2.2 激波管流動問題

激波管流動問題常用于檢測算法對激波，接觸間斷等的捕捉能力。本算例計算區(qū)域為［0 ，1］×［0 ，0 .1］的矩形區(qū)域，其中布置了76090個不規(guī)則分布的節(jié)點。初始條件設置為［14］：

本文分別采用基于線性逼近重構和三階WENO重構的無網(wǎng)格算法進行了計算，算至時刻t＝0.0039。圖6為流體密度在激波和接觸間斷附近的計算結果，圖7則為對應的在膨脹波附近的計算結果。整體上，基于三階WENO重構的無網(wǎng)格算法捕捉的激波、接觸間斷和膨脹波較線性逼近重構更靠近精確值，與文獻［14］五階有限差分 WENO重構（x方向點的間距與本文的相同）結果接近（見圖6、圖7），且在膨脹波附近也沒有明顯的振蕩，體現(xiàn)出WENO重構 “基本無振蕩”的特性（見圖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的計算區(qū)域，布置了101×61個規(guī)則分布的節(jié)點。外邊界采用自由來流條件，物面采用無穿透條件，上下出流邊界采用外推法處理［15］，圖8分別給出了無網(wǎng)格算法基于線性逼近重構以及WENO重構的壓強等值線，圖中可看出，兩種方法都捕捉到了脫體弓形激波，激波距駐點的距離都為0.68r（這與文獻［15］計算出的距離吻合），但基于WENO重構的無網(wǎng)格算法捕捉的激波更加清晰，表明后者三階WENO重構具有更高的分辨率。圖9給出了無網(wǎng)格算法基于線性逼近重構以及WENO重構的收斂曲線，整體看來，兩種方法的收斂曲線比較接近。另外，我們還對計算時間進行了統(tǒng)計，兩種方法迭代5000步所需的時間分別為197s和203s（CPU為Intel Core I7－860）。這一結果展示出基于WENO重構的無網(wǎng)格算法在計算效率上能與傳統(tǒng)基于線性逼近重構的無網(wǎng)格算法相當。

4.2.4 激波過彎道繞雙圓柱流場

在前面算例對發(fā)展的基于WENO重構的無網(wǎng)格算法進行驗證的基礎上，本節(jié)將該算法應用于求解激波過彎道繞雙圓柱流場，以展示其處理含非定常激波演化的復雜流動問題的效果。該算例幾何外形取自文［16］的實驗模型，本文計算捕捉到與實驗類似的波系演化結構，這里只給出計算結果。

圖8 計算得到的壓強等值線Fig.8 Comparison of the pressure contours obtained

圖9 收斂曲線Fig.9 Convergence history

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

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

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

5 結 論

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

［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ǎng)格算法及其應用研究［J］.空氣動力學學報，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）盛鳴劍，葉正寅，蔣超奇.附面層修正應用于無網(wǎng)格算法的研究［J］.計算力學學報，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）馬志華.自適應無網(wǎng)格及網(wǎng)格和無網(wǎng)格混合算法研究［D］.南京：南京航空航天大學，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ù)值測試與定量比較［J］.湘潭大學自然科學學報，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.