改進節點積分的無單元Galerkin法及其在流動問題中的應用

王玉龍，歐陽潔，王曉東，蔣 濤

（西北工業大學 理學院應用數學系，陜西 西安710129）

0 引 言

在流動問題中，對流擴散方程具有非常重要的地位。如何消除數值方法在對流占優時出現的數值偽振蕩，一直是學者們關心的重點。目前，基于網格的數值方法是數值求解流動問題的主流方法，如有限差分法（Finite Difference Method，FDM）、有限元法（Finite Element Method，FEM）和有限體積法（Finite Volume Method，FVM）等，但此類方法得到的數值結果與網格質量有關，且復雜三維結構的有限元網格剖分也非常困難。

無網格方法是近年來發展起來的一類新的數值方法，該類方法基于點近似，不需要將節點連成單元，克服了網格類方法遇到的一些困難。其中，無單元Galerkin（Element Free Galerkin，EFG）法［1－4］由 于精度高、穩定性好等優點而備受學者青睞。標準的EFG方法需要在背景網格上計算積分，由于無網格近似函數不是多項式，積分不能精確計算，一般需要用到較高階的 Gauss積分［2，4］，從而增加了計算量。另外，運用標準的EFG方法求解對流占優的流動問題時會出現非物理的數值偽振蕩現象［5］。

針對EFG方法計算量大的問題，一些學者提出了采用節點積分代替背景網格積分的方法。其中，Beissel等人［6］提出了直接節點積分方法，但該方法精度較低，穩定性較差。Nagashima［7］提出了基于局部Taylor展開的節點積分方法，解決了在結構分析中直接節點積分不穩定的問題，但Nagashima的方法在計算積分權值時要構造類似背景網格的“Bucket”。Liu等人［8］提出了基于全局Taylor展開的節點積分方法，然而該方法需要計算形函數的三階導數。Dyka等人［9］提出了應力點積分，該方法需要添加輔助應力點，計算量是直接節點積分的兩倍左右。Chen等人［10］提出了穩定協調的節點積分方法，該方法需要在各節點子域的邊界上進行積分，增加了節點積分的計算量。此外，針對EFG方法求解對流占優問題時出現的數值偽振蕩現象，歐陽潔等人［5］建立了一系列穩定化方案，但在穩定化方案中需要確定穩定化參數。

基于目前的研究現狀，本文采用局部Taylor展開思想，建立了局部Taylor展開積分無單元Galerkin（Local Taylor Expansion Integration Element Free Galerkin，LTEI－EFG）法，并首次將該方法用于求解對流占優的流體力學問題。與采用背景網格積分的EFG方法相比，LTEI－EFG法能夠明顯提高計算效率且有效抑制由于對流占優引起的數值偽振蕩。

1 控制方程和EFG方法的基本原理

1.1 控制方程

非穩態對流擴散方程的一般形式為［11］：

其中，v（x，t）是速度場，k（x，t）是擴散系數，s（x，t）是源項，若u（x，t）與時間t無關，則方程（1）是穩態對流擴散方程。

1.2 EFG方法的基本原理

在函數近似和空間離散時分別采用移動最小二乘法（Moving Least Square，MLS）和Galerkin法，則得到EFG方法［1］。下面針對方程（1）給出EFG方法的基本原理。

1.2.1 解的近似方案

設函數u（x，t）的MLS近似表達式為：

其中a（x，t）是待定向量，p（x）為基函數向量，m為基函數個數?；瘮迪蛄縫（x）通常采用Pascal三角形所決定的單項式以確保其最小完備性［2］，本文采用線性基，即pT（x）＝（1，x）（m＝2，in1D），pT（x）＝（1，x，y）（m＝3，in2D）。

式（2）中的待定向量a（x，t）滿足離散加權L2范數達到最小，即a（x，t）由

確定。其中n是點x支持域內的節點個數，ui（t）是節點參數，w（x－xi）是緊支權函數。局部Taylor展開積分需要計算形函數的二階導數，而權函數的連續性直接影響到形函數的連續性［1－4］，本文取三次樣條函數［2，4］作為權函數。

式（3）關于a（x，t）最小化得如下線性方程組：

其中A＝pT（x）W（x）p（x）

如果A（x）可逆，則a（x，t）＝A－1（x）B（x）u（t），將a（x，t）代入式（2），得：

其中N（x）是MLS形函數矩陣。

1.2.2 方程離散方案

為簡單起見，本文采用θ加權法［9］進行時間離散，其它時間離散方法可參見有關文獻［9］。

應用θ加權法，方程（1）可寫為以下形式：

其中 Δu＝un＋1－un，v＝（un，vn）。當θ＝1／2時為Crank－Nicolson格式，該格式在時間域上具有二階精度。本文算例均取θ＝1／2。

運用Galerkin原理進行空間離散，則方程（6）的弱形式為：

其中ω∈是檢驗函數，取為MLS形函數，

2 局部Taylor展開的節點積分方案

本文基于局部Taylor展開計算式（7）中的積分，并將采用局部Taylor展開積分的EFG方法記為LTEI－EFG。

2.1 局部Taylor展開

2.2 局部Taylor展開的節點積分

由式（7）知，需要計算下列積分：

其中Ω＝∪，Ωi∩Ωj＝Φ，i≠j，且每個Ωi（i＝1，2，…，N）包含一個節點，式（9）是時間項積分，式（10）和式（12）是擴散項積分，式（11）和式（13）是對流項積分。

對任意節點 （xi，yi）∈Ωi，下面分別以y）dΩ和（x，y）dΩ為例介紹局部 Taylor展開的節點積分，其計算公式如下：

定義式（14）～式（17）中的積分權系數如下：

權系數的確定將在2.3節介紹。

局部Taylor展開積分公式中包含了形函數的二階導數。一方面，該二階導數是在對被積函數做一階近似時出現的，與直接節點積分的零階近似相比，提高了計算精度；另一方面，二階導數具有穩定化的作用，能夠有效消除對流占優導致的數值偽振蕩。

2.3 權系數的確定

本文采用矩形支持域，且設x方向和y方向平均節點間距均為Δx。對任意節點（xi，yi），其影響半徑為dri，支持域半徑為dri·Δx，當（xi，yi）是內節點時取對應的子區域Ωi為圓形，當（xi，yi）是邊界節點時取對應的子區域Ωi為扇形，則節點（xi，yi）處權系數W1的計算公式如下［6］：

其中MΩ是區域Ω的面積，Ν′是節點個數，fi＝，對內節點θ＝2π，θ＝0，對邊界節點θ和θ FBFB是扇形區域對應的角度，如圖1所示［7］。

圖1 權系數示意圖Fig.1 Schematic diagram of the weight coefficient

下面采用文獻［7］的方法計算W2：

同理可得W3、W4、W5和W6的計算公式如下：

在式（19）～式（23）中，R

3 數值算例

下面通過兩個數值算例檢驗LTEI－EFG法和采用直接節點積分的EFG方法的計算效果。將采用直接節點積分的EFG方法記為NI－EFG。本文算例均在Intel Pentium 4CPU 2.60GHz PC機上計算。

3.1 一維定常對流擴散問題

考慮一維定常對流擴散方程［11］

其精確解為u（x）

圖2給出了分別取a＝1.0、v＝0.1和a＝1.0、v＝0.01時用EFG方法、NI－EFG法和LTEI－EFG法計算的數值解和精確解的比較。計算時在［0，1］區間內均勻布置31個節點，EFG方法用30個積分背景網格，每個積分背景網格用5點Gauss積分，影響半徑為1.2，NI－EFG法和LTEI－EFG法的影響半徑均為1.5。由圖2可見，當v較大時，EFG方法和LTEIEFG法的數值解與精確解吻合很好，而NI－EFG法的數值解出現不穩定現象［6］；當v較小時，EFG方法和NI－EFG法的數值解均出現了數值偽振蕩，且NI－EFG法的偽振蕩更加劇烈，而LTEI－EFG法的數值解并未出現偽振蕩。因此可知，（1）擴散占優時EFG方法有很好的計算精度，但對流占優時EFG方法的數值解會出現數值偽振蕩［12］；（2）擴散占優時 NIEFG法的數值解存在不穩定問題，對流占優時NIEFG法的數值解會出現數值偽振蕩，且比EFG方法的數值偽振蕩更劇烈；（3）擴散占優或對流占優時LTEI－EFG法的數值解均未出現數值偽振蕩，與精確解一致。

圖2 EFG方法、NI－EFG法和LTEI－EFG法的數值解和精確解的比較Fig.2 Comparison of between exact solution and numerical solution for EFG、NI－EFG and LTEI－EFG

表1比較了取不同節點數計算時EFG方法、NIEFG法和LTEI－EFG法組裝剛度矩陣花費的時間，其中EFG方法分別采用3點和5點Gauss積分。由表1可見，LTEI－EFG法和NI－EFG法組裝剛度矩陣花費的時間基本相同，且均小于采用3點和5點Gauss積分的EFG方法。

表1 不同節點數時組裝剛度矩陣的CPU時間（單位：ms）Table 1 CPU time of estimating stiffness matrix for different nodes（unit：ms）

為了對LTEI－EFG法進行精度分析，采用如下的相對誤差［2］：

式中為節點i處的函數數值解，ui為節點i處的函數精確解。

取v＝0.1，在［0，1］區間內分別均勻布置11、21、41、81、161、321、641和1281個節點計算，結果如表2所示，其中收斂率計算公式見文獻［2］。由表2可見，采用線性基的LTEI－EFG法仍然具有很高的收斂率，這是由于權函數的使用使得MLS近似形函數具有高階光滑性［2，4］。

表2 LTEI－EFG法的收斂率Table 2 Rate of convergence for LTEI－EFG method

3.2 二維非定常對流擴散問題

考慮二維Burgers方程［11］：

其中u和v分別是沿x軸和y軸的速度，Re是Reynolds數，（x，y，t）∈Ω×（0，T］，Ω＝（0，1）×（0，1），該方程具有如下對稱性：u（x，y，t）＝v（y，x，t），u（x，y，t）＝－u（1－x，1－y，t），下面只給出u的計算結果。

圖3～圖5分別給出了Re＝200時，用EFG方法、NI－EFG法和LTEI－EFG法求解時在不同時刻的數值解。此時Burgers方程中的對流項占優，在過點（0，1），（1，0）的對角線上將會產生間斷。

圖3 EFG方法在不同時刻的數值解（Re＝200）Fig.3 Numerical solution of EFG method for different times（Re＝200）

計算時均勻布置21×21個節點，時間步長為0.01s，EFG方法用3×3Gauss積分，影響半徑為1.2，NI－EFG法和LTEI－EFG法的影響半徑均為1.5。由圖3和圖4可見，EFG方法和NI－EFG法的數值解均出現了非物理的數值偽振蕩，且NI－EFG的偽振蕩更加劇烈。由圖5可見，LTEI－EFG法明顯抑制了數值偽振蕩，其數值解與文獻［12］的結果吻合。

圖4 NI－EFG法在不同時刻的數值解（Re＝200）Fig.4 Numerical solution of NI－EFG method at different times（Re＝200）

圖5 LTEI－EFG法在不同時刻的數值解（Re＝200）Fig.5 Numerical solution of LTEI－EFG method at different times（Re＝200）

表3比較了EFG方法、NI－EFG法和LTEI－EFG法在不同時間步長下計算到t＝1.0s時組裝剛度矩陣花費的時間，其中均勻布置21×21個節點，EFG方法分別采用3×3和5×5Gauss積分。由表3可見：（1）隨著時間步長的減小，組裝剛度矩陣的次數增加，因而花費的時間越多；（2）LTEI－EFG法和 NI－EFG法組裝剛度矩陣花費的時間基本相同，且均遠小于采用3×3點和5×5點Gauss積分的EFG方法。

表3 不同時間步長時組裝剛度矩陣的CPU時間（單位：s）Table 3 CPU time of estimating stiffness matrix for different time step（unit：s）

由表1和表3可見，一維情況下LTEI－EFG法需要的計算時間大約是采用3點Gauss積分EFG方法的40%；二維情況下LTEI－EFG法需要的計算時間大約是采用3×3Gauss積分EFG方法的18%。因而，LTEI－EFG法大幅度提高了采用背景網格積分EFG方法的計算效率，且問題維數越高，計算效率提高幅度越大。

為了進一步檢驗LTEI－EFG法在求解強對流問題時的有效性，取Re＝1000，均勻布置41×41個節點進行計算。圖6給出了計算到t＝1.0s時的數值解，圖7給出了沿直線y＝0.0處在不同時刻的數值解。由圖6和圖7可見，對于充分大的Re數，LTEIEFG法仍然能夠消除數值偽振蕩，其數值解與文獻［11－12］的結果一致，且LTEI－EFG法消除偽振蕩的效果更好。

圖6 t＝1.0s時LTEI－EFG法的數值解Fig.6 Numerical solution of LTEI－EFG method at t＝1.0s（Re＝1000）

圖7 不同時刻時沿直線y＝0的數值解Fig.7 Numerical solution along the line y＝0at different time

4 結 論

建立了局部Taylor展開積分無單元Galerkin法，并對一維定常對流擴散方程和二維Burgers方程進行了求解。所得結論如下：

（1）求解流動問題時，LTEI－EFG法的計算效率明顯高于采用背景網格積分的EFG方法。

（2）采用直接節點積分的EFG方法求解對流占優問題時會出現數值偽振蕩，且比用背景網格積分的EFG方法數值偽振蕩更加劇烈。而LTEI－EFG法可以很好地消除由于對流占優引起的數值偽振蕩，并且該方法不需要選擇穩定化參數。

（3）LTEI－EFG法完全不需要網格，是一種純無網格方法。

［1］BELYTSCHKO T，LU Y Y，GU L.Element－free galerkin methods［J］.Int.J.Numer.MethodsEng，1994，37：229－256.

［2］LIU G R，GU Y T.An introduction to meshfree methods and their programming［M］.Springer，2005.

［3］李九紅，程玉民.無網格方法的研究進展與展望［J］.力學季刊，2006，24：143－152.

［4］張雄，劉巖.無網格法［M］.清華大學出版社，Springer，2004.

［5］歐陽潔，張林，張小華.求解非穩態對流占優問題的非標準無網格Galerkin方法［J］.空氣動力學學報，2007，25：287－293.

［6］BEISSEL S，BELYTSCHKO T.Nodal integration of the element－free Galerkin method［J］.Comput.Methods Appl.Mech.Eng.，1996，139：49－74.

［7］NAGASHIMA T.Node－by－node meshless approach and its application to structural analyses［J］.Int.J.Numer.MethodsEng.，1999，46：341－385.

［8］LIU G R，ZHANG G Y，WANG Y Y，et al.A nodal integration technique for meshfree radial point interpolation method（NI－RPIM）［J］.Int.J.SolidsandStructures，2007，44：3840－3860.

［9］DYKA C T，INGEL R P.An approach for tensile instablity in smoothed particle hydrodynamics［J］.ComputStruct.，1994，57：573－580.

［10］CHEN J S，PAN C，WU C T.A Lagrangian reproducing kernel particle method for metal forming analysis［J］.Comput.Mech.，1998，22：289－307.

［11］DONEA J，HUERTA A.Finite element methods for flow problems［M］.Wiley，2003.

［12］ZHANG X H，OUYANG J，ZHANG L.Element－Free Charateristic Galerkin method for Burgers'equation［J］.Eng.Anal.BoundaryElem.，2009，33：356－362.