袁行飛，方少文，錢若軍

（1.浙江大學(xué) 空間結(jié)構(gòu)研究中心，浙江 杭州310058；2.同濟(jì)大學(xué) 土木工程學(xué)院，上海200092）

流體力學(xué)有限元法是流體力學(xué)數(shù)值計(jì)算中的一種重要方法，其主要是采用了標(biāo)準(zhǔn)伽遼金有限元格式，并且采用了傳統(tǒng)的多項(xiàng)式插值函數(shù)來得到形函數(shù).由于流體Navier-Stokes方程的解析解至今無法得到，不少學(xué)者對簡單的一維對流擴(kuò)散方程進(jìn)行研究，將一維對流擴(kuò)散方程蛻化為常系數(shù)微分方程，以常微分方程的有限元解與解析解對比，發(fā)現(xiàn)了有限元解的波動(dòng)［1-2］.從某種意義上說，流體力學(xué)有限元法的發(fā)展過程也是如何解決數(shù)值波動(dòng)的過程，很多方法都是圍繞著解的波動(dòng)問題提出的.國內(nèi)外學(xué)者對流體力學(xué)有限元法進(jìn)行了長期研究，并取得了豐碩成果，如美國麻省理工學(xué)院教授 Klaus-Jürgen Bathe領(lǐng)導(dǎo)的研究小組提出的基于流動(dòng)條件的插值方法（flow-condition-based interpolation，F(xiàn)CBI）已廣泛應(yīng)用于土木、機(jī)械等領(lǐng)域［3-5］.英國Swansea大學(xué)以Zienkiewicz教授為首的研究小組提出的基于特 征 線 的 分 裂 算 法 （Characteristic-Based Split，CBS）［6-7］現(xiàn)在也已得到了較為廣泛的應(yīng)用.Heinrich和Zienkiewicz等人在1976年提出的迎風(fēng)有限元法，經(jīng)過不斷的修正和發(fā)展，現(xiàn)在已應(yīng)用到三維流體力學(xué)的計(jì)算［8-10］.國內(nèi)較早對流體力學(xué)有限元法展開系統(tǒng)研究的學(xué)者有章本照［11］、劉希云［12］、楊曜根［13］等.

本文對一維定常對流擴(kuò)散方程標(biāo)準(zhǔn)伽遼金有限元解的波動(dòng)問題進(jìn)行了研究，指出除單元內(nèi)插值函數(shù)連續(xù)性差外，單元間形函數(shù)一階導(dǎo)數(shù)不連續(xù)，即屬于不同單元同一節(jié)點(diǎn)的左右一階導(dǎo)數(shù)不相等是導(dǎo)致有限元解波動(dòng)的另一主要原因.進(jìn)而提出采用補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元進(jìn)行分析，通過在單元與單元之間的節(jié)點(diǎn)上增加一個(gè)“不平衡力”，即在單元間引入補(bǔ)償項(xiàng)來提高單元間形函數(shù)一階導(dǎo)數(shù)的連續(xù)性，來減少數(shù)值波動(dòng).針對線性拉格朗日插值和指數(shù)型插值分別探討了補(bǔ)償項(xiàng)中補(bǔ)償剛度表達(dá)式，研究了不同修正系數(shù)下有限元解的波動(dòng)情況，確定了最優(yōu)修正系數(shù).以常見流體為例對比了采用指數(shù)型插值函數(shù)的補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元解與解析解計(jì)算結(jié)果，驗(yàn)證了其有效性.本文工作可為進(jìn)一步研究和發(fā)展三維流體力學(xué)有限元法提供思路.

1 一維對流擴(kuò)散方程的有限元解波動(dòng)

1.1 解析解

對于在L內(nèi)定義，待求未知變量為φ的一維定常對流擴(kuò)散方程

式中：ν為運(yùn)動(dòng)黏性系數(shù)；u為速度；φ為因變量；x為坐標(biāo)，0≤x≤L.

式（1）是變系數(shù)偏微分方程，通常難以求解.作為問題的近似，假設(shè)u為常數(shù)，則式（1）蛻化為常系數(shù)微分方程.解該方程得

式中：D1，D2均為待定常數(shù)，通過引入邊界條件φ|x=0=φ0和φ|x=L=φL可惟一確定.式（1）的解析解為

1.2 常用有限元解

現(xiàn)采用標(biāo)準(zhǔn)伽遼金有限元進(jìn)行求解.首先對定義域L進(jìn)行剖分，共劃分J個(gè)線單元，單元長度為Δx，Δx=L／J，則有K（1，2，…，J+1）個(gè)節(jié)點(diǎn).采用一維線性拉格朗日插值函數(shù)，則式（1）的有限元解為

式中：Pe稱為貝克萊（Péclet）數(shù)，Pe=uΔx／2ν；c1，c2均為待定常數(shù)，通過引入邊界條件φ|x=0=φj和φ|x=Δx=φj+1可確定.式（1）的有限元解為

1.3 有限元解波動(dòng)

取L=1，J=10，ν=0.2，對流速度u分別為0.04和20時(shí)，可得相應(yīng)單元的Pe和φ值.如將每個(gè)節(jié)點(diǎn)的坐標(biāo)值代入式（3），則可得到相應(yīng)的解析解.

圖1為該問題的解析解和常用有限元解對比.由圖1可見，當(dāng)Pe較小時(shí)，有限元解和解析解十分接近；隨著Pe的增加，有限元解與解析解差異逐漸增大.當(dāng)Pe＞1時(shí)，有限元解在解析解附近來回跳動(dòng)，出現(xiàn)波動(dòng)現(xiàn)象［2］.

圖1 一維對流擴(kuò)散方程解析解和有限元解對比Fig.1 Comparison of analytical solution and finite element solution

通常將一維對流擴(kuò)散方程的有限元解產(chǎn)生數(shù)值波動(dòng)原因歸納為：網(wǎng)格的尺寸、對流速度以及黏性系數(shù).已有的解決方法大都針對一維線性拉格朗日插值函數(shù)進(jìn)行研究，提出加密網(wǎng)格、修正運(yùn)動(dòng)黏性系數(shù)、修正對流速度等方法來改善有限元數(shù)值波動(dòng).作者認(rèn)為引起數(shù)值波動(dòng)的兩個(gè)本質(zhì)原因是線性拉格朗日插值函數(shù)連續(xù)性差以及單元之間形函數(shù)一階導(dǎo)數(shù)不連續(xù).文獻(xiàn)［14］對第一個(gè)原因進(jìn)行了探索，研究結(jié)果表明采用連續(xù)性較好的插值函數(shù)（如三次多項(xiàng)式埃爾米特插值函數(shù)、指數(shù)型插值函數(shù)等）是改善數(shù)值波動(dòng)的有效方法.但這些插值函數(shù)仍不能保證節(jié)點(diǎn)處不同單元的連續(xù)性.以一維問題為例，同一節(jié)點(diǎn)屬于前后兩個(gè)單元，前一單元右節(jié)點(diǎn)的一階導(dǎo)數(shù)不等于后一單元左節(jié)點(diǎn)的一階導(dǎo)數(shù).為解決這一問題，本文通過引入補(bǔ)償項(xiàng)來改善單元之間節(jié)點(diǎn)處的連續(xù)性，從而改善一維對流擴(kuò)散方程標(biāo)準(zhǔn)伽遼金有限元數(shù)值波動(dòng)，為研究和發(fā)展三維流體力學(xué)有限元法提供思路.

2 補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元基本原理

2.1 標(biāo)準(zhǔn)伽遼金有限元格式

方程（1）對流擴(kuò)散微分方程可寫成等價(jià)的變分公式，當(dāng)對此變分公式采用標(biāo)準(zhǔn)伽遼金有限元法求解時(shí)，需先給出如下變分公式的弱形式：

采用分部積分，并引入Gauss-Green公式，對含有函數(shù)二階導(dǎo)數(shù)的擴(kuò)散項(xiàng)降階，得

因此一維對流擴(kuò)散方程的標(biāo)準(zhǔn)伽遼金有限元格式為

式中：N為單元的形函數(shù)；φ=Nφe；φe為函數(shù)φ在單元節(jié)點(diǎn)處的值.

2.2 補(bǔ)償有限元基本原理及格式

對一維對流擴(kuò)散問題的連續(xù)性進(jìn)行分析，發(fā)現(xiàn)引起數(shù)值波動(dòng)一個(gè)重要的原因是由于單元之間函數(shù)的一階導(dǎo)數(shù)不連續(xù).如圖2所示，為了提高單元間函數(shù)一階導(dǎo)數(shù)的連續(xù)性，減小波動(dòng)，可在單元之間即節(jié)點(diǎn)處引入補(bǔ)償項(xiàng)

式中：λ為補(bǔ)償剛度，其物理意義類似于固體力學(xué)中發(fā)生單位轉(zhuǎn)角所需施加的彎矩；dφ／dx為單元間的轉(zhuǎn)角.針對不同的插值函數(shù)，補(bǔ)償剛度λ表達(dá)式可以不同.

上述引入補(bǔ)償項(xiàng)的方法可以不改變原來插值函數(shù)，也不會(huì)增加單元的自由度及內(nèi)插點(diǎn)數(shù)，而僅在單元與單元之間的節(jié)點(diǎn)上增加一個(gè)“不平衡力”，從而使得單元之間的連續(xù)性得到改善，較好地控制有限元解的波動(dòng)問題.

圖2 單元間的彈性約束Fig.2 Elastic constraint between different elements

引入補(bǔ)償項(xiàng)后的補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元格式為

式（10）右端項(xiàng)φe為前一步計(jì)算得到的節(jié)點(diǎn)值，通過不斷的迭代求解，φ最終收斂到較精確的數(shù)值解.

3 采用線性拉格朗日插值的補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元

3.1 線性拉格朗日插值函數(shù)

為進(jìn)行一般性研究，剖分定義域L為J個(gè)線單元，單元長度為Δx=L／J，則有K（1，2，…，J+1）個(gè)節(jié)點(diǎn).單元內(nèi)各點(diǎn)的值可用單元節(jié)點(diǎn)處值表示

3.2 補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元補(bǔ)償剛度

由于單元之間采用的是拉格朗日插值，其一階導(dǎo)數(shù)不連續(xù)，導(dǎo)致Pe＞1時(shí)有限元解存在較大波動(dòng)（圖1）.為了增強(qiáng)單元間函數(shù)一階導(dǎo)數(shù)的連續(xù)性，減小波動(dòng)，考慮在單元之間引入包含補(bǔ)償剛度λ的補(bǔ)償項(xiàng).參考迎風(fēng)有限元格式［8-10］，初步確定線性拉格朗日插值補(bǔ)償項(xiàng)中補(bǔ)償剛度λ的表達(dá)式如下：

式中：Pe=uΔx／2ν；k為修正系數(shù).

研究表明，數(shù)值波動(dòng)程度與k值有關(guān).當(dāng)k值較小時(shí)，數(shù)值解存在一定程度的波動(dòng)，隨著k值的增加數(shù)值波動(dòng)逐漸減緩；但當(dāng)k值過大時(shí)，數(shù)值解又逐漸偏離解析解.為了獲得最優(yōu)的修正系數(shù)k，令有限元值與解析解的誤差為.其中φi，φ′分別為解析解和采用補(bǔ)償線性拉格朗日插值的有限元解.

以u=8，L=1，J=10為例，當(dāng)Pe分別為5和100時(shí)，采用不同修正系數(shù)k的補(bǔ)償線性拉格朗日插值有限元解波動(dòng)情況見圖3和圖4.

圖3 Pe＝5時(shí)補(bǔ)償線性拉格朗日插值有限元解波動(dòng)Fig.3 Solution wave of CSGFE using LLBI at Pe＝5

由圖3，4可以看出，當(dāng)Pe較大時(shí)，線性拉格朗日插值出現(xiàn)了嚴(yán)重的波動(dòng)現(xiàn)象，而引入補(bǔ)償項(xiàng)以后，波動(dòng)情況得到明顯改善.對于不同的調(diào)整系數(shù)k，波動(dòng)的改善情況有所不同.當(dāng)k＜1.0時(shí)，隨著k的增加，補(bǔ)償線性拉格朗日插值的誤差越來越小；k＞1.0時(shí)，誤差又急劇上升；而當(dāng)k=1.0時(shí)，線性拉格朗日插值的誤差最小，可以獲得與解析解最為接近的數(shù)值解.因此，取k=1.0為補(bǔ)償線性拉格朗日插值的最優(yōu)修正系數(shù).此時(shí)補(bǔ)償線性拉格朗日插值能提高單元之間的連續(xù)性，改善波動(dòng)情況，取得較好的數(shù)值解.

由此確定采用線性拉格朗日插值的補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元補(bǔ)償剛度λ的表達(dá)式為

4 采用指數(shù)型插值的補(bǔ)償伽遼金有限元

4.1 指數(shù)型插值函數(shù)

設(shè)有限元解為φ=c1+c2eax，單元形函數(shù)為N=.其中a的取值與流場屬性有關(guān)［14］.取不同a值時(shí)指數(shù)型插值函數(shù)有限元解波動(dòng)見圖5.由圖5可見，線性拉格朗日插值有較大波動(dòng)，指數(shù)型插值的波動(dòng)情況與a值有關(guān).隨著a增大，有限元解波動(dòng)減小.當(dāng)a=u／v=80時(shí)，有限元計(jì)算結(jié)果與解析解基本一致.

圖4 Pe＝100時(shí)補(bǔ)償線性拉格朗日插值有限元解波動(dòng)Fig.4 Solution wave of CSGFE using LLBI at Pe＝100

圖5 不同a值時(shí)指數(shù)型插值函數(shù)有限元解Fig.5 Finite element solution using EFBI with different a

4.2 補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元補(bǔ)償剛度

指數(shù)插值函數(shù)在單元內(nèi)更光滑，并且多階連續(xù)，因而可以獲得問題更好的數(shù)值解.但采用不同a值時(shí)，指數(shù)插值函數(shù)仍會(huì)出現(xiàn)不同幅度的波動(dòng).其原因仍在于單元之間采用的是拉格朗日插值，其一階導(dǎo)數(shù)不連續(xù).同理，為增強(qiáng)單元間函數(shù)一階導(dǎo)數(shù)的連續(xù)性，考慮在單元之間引入包含補(bǔ)償剛度的補(bǔ)償項(xiàng)，建議如下補(bǔ)償剛度λ的表達(dá)式：

式中：Re=φΔx／2ν；k為修正系數(shù).數(shù)值波動(dòng)程度與k取值有關(guān).

以u=5，L=1，J=10為例，當(dāng)Pe分別為25，2 500時(shí)，采用不同修正系數(shù)k的補(bǔ)償指數(shù)型插值有限元解波動(dòng)情況見圖6和圖7.

圖6 Pe＝25時(shí)補(bǔ)償指數(shù)型插值有限元解波動(dòng)Fig.6 Solution wave of CSGFE using FEBI at Pe＝25

由圖6，7可知，與線性拉格朗日插值相比，補(bǔ)償指數(shù)型插值在單元內(nèi)和單元之間都有較好的連續(xù)性，能夠較顯著改善波動(dòng)情況，取得較好的數(shù)值結(jié)果.而對于不同的修正系數(shù)，當(dāng)k＜0.5時(shí)，隨著k的增加，補(bǔ)償指數(shù)型插值的誤差越來越小；當(dāng)k＞0.5時(shí)，誤差不斷增加；當(dāng)k=0.5時(shí)，補(bǔ)償指數(shù)型插值的誤差最小，可以獲得與解析解最為接近的數(shù)值解.因此取k=0.5為補(bǔ)償指數(shù)型插值的最優(yōu)修正系數(shù).此時(shí)補(bǔ)償指數(shù)型插值能提高單元之間的連續(xù)性，改善波動(dòng)情況，取得較好的數(shù)值解.

由此確定采用指數(shù)型插值的補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元補(bǔ)償剛度λ的表達(dá)式為

圖7 Pe＝2 500×103時(shí)補(bǔ)償指數(shù)型插值有限元解波動(dòng)Fig.7 Solution wave of CSGFE using FEBI at Pe＝2.5×103

圖8 不同流體的補(bǔ)償指數(shù)型插值有限元解波動(dòng)Fig.8 Solution wave of CSGFE using EBI of different fluids

采用公式（14），當(dāng)u=10，L=1，J=10時(shí)，對常溫（27℃）下空氣（νa=1.6×10-5m2·s-1）和水（νw=8.6×10-7m2·s-1）兩種常見流體的補(bǔ)償指數(shù)型插值有限元解波動(dòng)進(jìn)行了研究，結(jié)果見圖8.由圖可知，對于常見流體空氣和水，補(bǔ)償指數(shù)型插值均能取得較為理想的結(jié)果，與精確解非常接近.

5 結(jié)論

（1）本文對一維定常對流擴(kuò)散方程標(biāo)準(zhǔn)伽遼金有限元解的波動(dòng)問題進(jìn)行了研究，指出提高單元內(nèi)插值函數(shù)連續(xù)性和單元之間形函數(shù)連續(xù)性是改善一維對流擴(kuò)散方程有限元解數(shù)值波動(dòng)問題的有效途徑.

（2）本文提出的補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元在原來插值函數(shù)的基礎(chǔ)上，不增加單元的自由度及內(nèi)插點(diǎn)數(shù)，只是增加了補(bǔ)償項(xiàng)，即在單元與單元之間的節(jié)點(diǎn)上增加一個(gè)“不平衡力”，就可使得單元之間的連續(xù)性得到改善，從而更好地控制數(shù)值波動(dòng)現(xiàn)象，取得問題較好的數(shù)值解.

（3）相比于常規(guī)多項(xiàng)式插值需增加插值函數(shù)的階次來提高計(jì)算精度，本文提出的補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元格式不增加單元的自由度及內(nèi)插點(diǎn)數(shù)，即在不明顯增加計(jì)算量的前提下使得計(jì)算精度得到大幅提高.

（4）本文探討了采用線性拉格朗日插值和指數(shù)型插值的補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元補(bǔ)償剛度表達(dá)式.相比標(biāo)準(zhǔn)伽遼金有限元，補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元均能提高單元之間的連續(xù)性，顯著改善波動(dòng)情況.

（5）本文研究結(jié)果表明采用指數(shù)型插值函數(shù)的補(bǔ)償標(biāo)準(zhǔn)伽遼金有限元不僅在單元內(nèi)可精確給出變量的分布，而且單元之間的連續(xù)性更好，因而能更好地控制數(shù)值波動(dòng)現(xiàn)象，取得問題較好的數(shù)值解.本文工作可為進(jìn)一步研究和發(fā)展三維流體力學(xué)有限元法提供思路.

［1］ Zienkiewicz O C，Taylor R L.The finite element method for fluid dynamics［M］.5th ed.Elsevier：Amsterdam，2000.

［2］ Belytschko T，Liu W K，Moran B.Nonlinear finite elements for continua and structures［M］.New York：John Wiley &Sons Limited，2000.

［3］ Bathe K J.Finite element procedures［M］.New Jersey：Prentice Hall，1996：1-25.

［4］ Bathe K J，ZHANG Hou.A flow-condition-based interpolation finite element procedure for incompressible fluid flows［J］.Computers and Structures，2002，80（14）：1267.

［5］ Kohno H，Bathe K J.A nine-node quadrilateral FCBI element for incompressible fluid flows［J］.Communications in Numerical Methods in Engineering，2006，22（8）：917.

［6］ Zienkiewicz O C， Wu J.A general explicit or semi-explicit algorithm for compressible and incompressible flows ［J］.International Journal for Numerical Methods in Engineering，1992，35（3）：457.

［7］ Nithiarasu P，Codina R，Zienkiewicz O C.The characteristicbased split（CBS）scheme—a unified approach to fluid dynamics［J］.Numerical Methods in Engineering，2006，66（10）：1514.

［8］ Hughes T J R，Brooks A.A theoretical framework for Petrov-Galerkin methods with discontinuous weighting functions：Application to the streamline-upwind procedure［J］.Finite Elements in Fluids，1982，4（2）：47.

［9］ Burman E，Smith G.Analysis of the space semi-discretized SUPG method for transient convection—diffusion equations［J］.Mathematical Models and Methods in Applied Sciences，2011，21（10）：2049.

［10］ Takase S，Kashiyama K，Tanaka S，etal.Space-time SUPG finite element computation of shallow-water flows with moving shorelines［J］.Computational Mechanics，2011，48（3）：293.

［11］ 章本照.流體力學(xué)中的有限元方法［M］.北京：機(jī)械工業(yè)出版社，1986.ZHANG Benzhao.The finite element method in fluid mechanics［M］.Beijing：Mechanical Industry Press，1986.

［12］ 劉希云，趙潤祥.流體力學(xué)中的有限元與邊界元方法［M］.上海：上海交通大學(xué)出版社，1993.LIU Xiyun，ZHAO Runxiang.Finite element and boundary element methods in fluid mechanics［M］.Shanghai：Shanghai Jiaotong University Press，1993.

［13］ 楊曜根.流體力學(xué)有限元［M］.哈爾濱：哈爾濱工程大學(xué)出版社，1995.YANG Yaogen.Fluid dynamics finite element［M］.Harbin ：Harbin Engineering University Press，1995.

［14］ 方少文，袁行飛，錢若軍，等.采用不同插值函數(shù)的流體力學(xué)有限元數(shù)值波動(dòng)研究［J］.工程力學(xué)，2013，30（11）：266.FANG Shaowen，YUAN Xingfei，QIAN Ruojun，etal.Numerical wave of finite element solution in fluid mechanics using different interpolation function［J］.Engineering Mechanics，2013，30（11）：266.