改進(jìn)節(jié)點(diǎn)積分的無(wú)單元Galerkin法及其在流動(dòng)問(wèn)題中的應(yīng)用

王玉龍，歐陽(yáng)潔，王曉東，蔣 濤

（西北工業(yè)大學(xué) 理學(xué)院應(yīng)用數(shù)學(xué)系，陜西 西安710129）

0 引 言

在流動(dòng)問(wèn)題中，對(duì)流擴(kuò)散方程具有非常重要的地位。如何消除數(shù)值方法在對(duì)流占優(yōu)時(shí)出現(xiàn)的數(shù)值偽振蕩，一直是學(xué)者們關(guān)心的重點(diǎn)。目前，基于網(wǎng)格的數(shù)值方法是數(shù)值求解流動(dòng)問(wèn)題的主流方法，如有限差分法（Finite Difference Method，F(xiàn)DM）、有限元法（Finite Element Method，F(xiàn)EM）和有限體積法（Finite Volume Method，F(xiàn)VM）等，但此類方法得到的數(shù)值結(jié)果與網(wǎng)格質(zhì)量有關(guān)，且復(fù)雜三維結(jié)構(gòu)的有限元網(wǎng)格剖分也非常困難。

無(wú)網(wǎng)格方法是近年來(lái)發(fā)展起來(lái)的一類新的數(shù)值方法，該類方法基于點(diǎn)近似，不需要將節(jié)點(diǎn)連成單元，克服了網(wǎng)格類方法遇到的一些困難。其中，無(wú)單元Galerkin（Element Free Galerkin，EFG）法［1－4］由 于精度高、穩(wěn)定性好等優(yōu)點(diǎn)而備受學(xué)者青睞。標(biāo)準(zhǔn)的EFG方法需要在背景網(wǎng)格上計(jì)算積分，由于無(wú)網(wǎng)格近似函數(shù)不是多項(xiàng)式，積分不能精確計(jì)算，一般需要用到較高階的 Gauss積分［2，4］，從而增加了計(jì)算量。另外，運(yùn)用標(biāo)準(zhǔn)的EFG方法求解對(duì)流占優(yōu)的流動(dòng)問(wèn)題時(shí)會(huì)出現(xiàn)非物理的數(shù)值偽振蕩現(xiàn)象［5］。

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

基于目前的研究現(xiàn)狀，本文采用局部Taylor展開(kāi)思想，建立了局部Taylor展開(kāi)積分無(wú)單元Galerkin（Local Taylor Expansion Integration Element Free Galerkin，LTEI－EFG）法，并首次將該方法用于求解對(duì)流占優(yōu)的流體力學(xué)問(wèn)題。與采用背景網(wǎng)格積分的EFG方法相比，LTEI－EFG法能夠明顯提高計(jì)算效率且有效抑制由于對(duì)流占優(yōu)引起的數(shù)值偽振蕩。

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

1.1 控制方程

非穩(wěn)態(tài)對(duì)流擴(kuò)散方程的一般形式為［11］：

其中，v（x，t）是速度場(chǎng)，k（x，t）是擴(kuò)散系數(shù)，s（x，t）是源項(xiàng)，若u（x，t）與時(shí)間t無(wú)關(guān)，則方程（1）是穩(wěn)態(tài)對(duì)流擴(kuò)散方程。

1.2 EFG方法的基本原理

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

1.2.1 解的近似方案

設(shè)函數(shù)u（x，t）的MLS近似表達(dá)式為：

其中a（x，t）是待定向量，p（x）為基函數(shù)向量，m為基函數(shù)個(gè)數(shù)。基函數(shù)向量p（x）通常采用Pascal三角形所決定的單項(xiàng)式以確保其最小完備性［2］，本文采用線性基，即pT（x）＝（1，x）（m＝2，in1D），pT（x）＝（1，x，y）（m＝3，in2D）。

式（2）中的待定向量a（x，t）滿足離散加權(quán)L2范數(shù)達(dá)到最小，即a（x，t）由

確定。其中n是點(diǎn)x支持域內(nèi)的節(jié)點(diǎn)個(gè)數(shù)，ui（t）是節(jié)點(diǎn)參數(shù)，w（x－xi）是緊支權(quán)函數(shù)。局部Taylor展開(kāi)積分需要計(jì)算形函數(shù)的二階導(dǎo)數(shù)，而權(quán)函數(shù)的連續(xù)性直接影響到形函數(shù)的連續(xù)性［1－4］，本文取三次樣條函數(shù)［2，4］作為權(quán)函數(shù)。

式（3）關(guān)于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形函數(shù)矩陣。

1.2.2 方程離散方案

為簡(jiǎn)單起見(jiàn)，本文采用θ加權(quán)法［9］進(jìn)行時(shí)間離散，其它時(shí)間離散方法可參見(jiàn)有關(guān)文獻(xiàn)［9］。

應(yīng)用θ加權(quán)法，方程（1）可寫(xiě)為以下形式：

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

運(yùn)用Galerkin原理進(jìn)行空間離散，則方程（6）的弱形式為：

其中ω∈是檢驗(yàn)函數(shù)，取為MLS形函數(shù)，

2 局部Taylor展開(kāi)的節(jié)點(diǎn)積分方案

本文基于局部Taylor展開(kāi)計(jì)算式（7）中的積分，并將采用局部Taylor展開(kāi)積分的EFG方法記為L(zhǎng)TEI－EFG。

2.1 局部Taylor展開(kāi)

2.2 局部Taylor展開(kāi)的節(jié)點(diǎn)積分

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

其中Ω＝∪，Ωi∩Ωj＝Φ，i≠j，且每個(gè)Ωi（i＝1，2，…，N）包含一個(gè)節(jié)點(diǎn)，式（9）是時(shí)間項(xiàng)積分，式（10）和式（12）是擴(kuò)散項(xiàng)積分，式（11）和式（13）是對(duì)流項(xiàng)積分。

對(duì)任意節(jié)點(diǎn) （xi，yi）∈Ωi，下面分別以y）dΩ和（x，y）dΩ為例介紹局部 Taylor展開(kāi)的節(jié)點(diǎn)積分，其計(jì)算公式如下：

定義式（14）～式（17）中的積分權(quán)系數(shù)如下：

權(quán)系數(shù)的確定將在2.3節(jié)介紹。

局部Taylor展開(kāi)積分公式中包含了形函數(shù)的二階導(dǎo)數(shù)。一方面，該二階導(dǎo)數(shù)是在對(duì)被積函數(shù)做一階近似時(shí)出現(xiàn)的，與直接節(jié)點(diǎn)積分的零階近似相比，提高了計(jì)算精度；另一方面，二階導(dǎo)數(shù)具有穩(wěn)定化的作用，能夠有效消除對(duì)流占優(yōu)導(dǎo)致的數(shù)值偽振蕩。

2.3 權(quán)系數(shù)的確定

本文采用矩形支持域，且設(shè)x方向和y方向平均節(jié)點(diǎn)間距均為Δx。對(duì)任意節(jié)點(diǎn)（xi，yi），其影響半徑為dri，支持域半徑為dri·Δx，當(dāng)（xi，yi）是內(nèi)節(jié)點(diǎn)時(shí)取對(duì)應(yīng)的子區(qū)域Ωi為圓形，當(dāng)（xi，yi）是邊界節(jié)點(diǎn)時(shí)取對(duì)應(yīng)的子區(qū)域Ωi為扇形，則節(jié)點(diǎn)（xi，yi）處權(quán)系數(shù)W1的計(jì)算公式如下［6］：

其中MΩ是區(qū)域Ω的面積，Ν′是節(jié)點(diǎn)個(gè)數(shù)，fi＝，對(duì)內(nèi)節(jié)點(diǎn)θ＝2π，θ＝0，對(duì)邊界節(jié)點(diǎn)θ和θ FBFB是扇形區(qū)域?qū)?yīng)的角度，如圖1所示［7］。

圖1 權(quán)系數(shù)示意圖Fig.1 Schematic diagram of the weight coefficient

下面采用文獻(xiàn)［7］的方法計(jì)算W2：

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

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

3 數(shù)值算例

下面通過(guò)兩個(gè)數(shù)值算例檢驗(yàn)LTEI－EFG法和采用直接節(jié)點(diǎn)積分的EFG方法的計(jì)算效果。將采用直接節(jié)點(diǎn)積分的EFG方法記為NI－EFG。本文算例均在Intel Pentium 4CPU 2.60GHz PC機(jī)上計(jì)算。

3.1 一維定常對(duì)流擴(kuò)散問(wèn)題

考慮一維定常對(duì)流擴(kuò)散方程［11］

其精確解為u（x）

圖2給出了分別取a＝1.0、v＝0.1和a＝1.0、v＝0.01時(shí)用EFG方法、NI－EFG法和LTEI－EFG法計(jì)算的數(shù)值解和精確解的比較。計(jì)算時(shí)在［0，1］區(qū)間內(nèi)均勻布置31個(gè)節(jié)點(diǎn)，EFG方法用30個(gè)積分背景網(wǎng)格，每個(gè)積分背景網(wǎng)格用5點(diǎn)Gauss積分，影響半徑為1.2，NI－EFG法和LTEI－EFG法的影響半徑均為1.5。由圖2可見(jiàn)，當(dāng)v較大時(shí)，EFG方法和LTEIEFG法的數(shù)值解與精確解吻合很好，而NI－EFG法的數(shù)值解出現(xiàn)不穩(wěn)定現(xiàn)象［6］；當(dāng)v較小時(shí)，EFG方法和NI－EFG法的數(shù)值解均出現(xiàn)了數(shù)值偽振蕩，且NI－EFG法的偽振蕩更加劇烈，而LTEI－EFG法的數(shù)值解并未出現(xiàn)偽振蕩。因此可知，（1）擴(kuò)散占優(yōu)時(shí)EFG方法有很好的計(jì)算精度，但對(duì)流占優(yōu)時(shí)EFG方法的數(shù)值解會(huì)出現(xiàn)數(shù)值偽振蕩［12］；（2）擴(kuò)散占優(yōu)時(shí) NIEFG法的數(shù)值解存在不穩(wěn)定問(wèn)題，對(duì)流占優(yōu)時(shí)NIEFG法的數(shù)值解會(huì)出現(xiàn)數(shù)值偽振蕩，且比EFG方法的數(shù)值偽振蕩更劇烈；（3）擴(kuò)散占優(yōu)或?qū)α髡純?yōu)時(shí)LTEI－EFG法的數(shù)值解均未出現(xiàn)數(shù)值偽振蕩，與精確解一致。

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

表1比較了取不同節(jié)點(diǎn)數(shù)計(jì)算時(shí)EFG方法、NIEFG法和LTEI－EFG法組裝剛度矩陣花費(fèi)的時(shí)間，其中EFG方法分別采用3點(diǎn)和5點(diǎn)Gauss積分。由表1可見(jiàn)，LTEI－EFG法和NI－EFG法組裝剛度矩陣花費(fèi)的時(shí)間基本相同，且均小于采用3點(diǎn)和5點(diǎn)Gauss積分的EFG方法。

表1 不同節(jié)點(diǎn)數(shù)時(shí)組裝剛度矩陣的CPU時(shí)間（單位：ms）Table 1 CPU time of estimating stiffness matrix for different nodes（unit：ms）

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

式中為節(jié)點(diǎn)i處的函數(shù)數(shù)值解，ui為節(jié)點(diǎn)i處的函數(shù)精確解。

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

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

3.2 二維非定常對(duì)流擴(kuò)散問(wèn)題

考慮二維Burgers方程［11］：

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

圖3～圖5分別給出了Re＝200時(shí)，用EFG方法、NI－EFG法和LTEI－EFG法求解時(shí)在不同時(shí)刻的數(shù)值解。此時(shí)Burgers方程中的對(duì)流項(xiàng)占優(yōu)，在過(guò)點(diǎn)（0，1），（1，0）的對(duì)角線上將會(huì)產(chǎn)生間斷。

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

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

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

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

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

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

由表1和表3可見(jiàn)，一維情況下LTEI－EFG法需要的計(jì)算時(shí)間大約是采用3點(diǎn)Gauss積分EFG方法的40%；二維情況下LTEI－EFG法需要的計(jì)算時(shí)間大約是采用3×3Gauss積分EFG方法的18%。因而，LTEI－EFG法大幅度提高了采用背景網(wǎng)格積分EFG方法的計(jì)算效率，且問(wèn)題維數(shù)越高，計(jì)算效率提高幅度越大。

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

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

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

4 結(jié) 論

建立了局部Taylor展開(kāi)積分無(wú)單元Galerkin法，并對(duì)一維定常對(duì)流擴(kuò)散方程和二維Burgers方程進(jìn)行了求解。所得結(jié)論如下：

（1）求解流動(dòng)問(wèn)題時(shí)，LTEI－EFG法的計(jì)算效率明顯高于采用背景網(wǎng)格積分的EFG方法。

（2）采用直接節(jié)點(diǎn)積分的EFG方法求解對(duì)流占優(yōu)問(wèn)題時(shí)會(huì)出現(xiàn)數(shù)值偽振蕩，且比用背景網(wǎng)格積分的EFG方法數(shù)值偽振蕩更加劇烈。而LTEI－EFG法可以很好地消除由于對(duì)流占優(yōu)引起的數(shù)值偽振蕩，并且該方法不需要選擇穩(wěn)定化參數(shù)。

（3）LTEI－EFG法完全不需要網(wǎng)格，是一種純無(wú)網(wǎng)格方法。

［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］李九紅，程玉民.無(wú)網(wǎng)格方法的研究進(jìn)展與展望［J］.力學(xué)季刊，2006，24：143－152.

［4］張雄，劉巖.無(wú)網(wǎng)格法［M］.清華大學(xué)出版社，Springer，2004.

［5］歐陽(yáng)潔，張林，張小華.求解非穩(wěn)態(tài)對(duì)流占優(yōu)問(wèn)題的非標(biāo)準(zhǔn)無(wú)網(wǎng)格Galerkin方法［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，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.