引入補償剛度的流體力學(xué)標(biāo)準(zhǔn)伽遼金有限元研究

袁行飛，方少文，錢若軍

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

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

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

1 一維對流擴散方程的有限元解波動

1.1 解析解

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

式中：ν為運動黏性系數(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)伽遼金有限元進行求解.首先對定義域L進行剖分，共劃分J個線單元，單元長度為Δx，Δx=L／J，則有K（1，2，…，J+1）個節(jié)點.采用一維線性拉格朗日插值函數(shù)，則式（1）的有限元解為

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

1.3 有限元解波動

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

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

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

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

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

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

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

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

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

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

2.2 補償有限元基本原理及格式

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖3 Pe＝5時補償線性拉格朗日插值有限元解波動Fig.3 Solution wave of CSGFE using LLBI at Pe＝5

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

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

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

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

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

圖4 Pe＝100時補償線性拉格朗日插值有限元解波動Fig.4 Solution wave of CSGFE using LLBI at Pe＝100

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

（5）本文研究結(jié)果表明采用指數(shù)型插值函數(shù)的補償標(biāo)準(zhǔn)伽遼金有限元不僅在單元內(nèi)可精確給出變量的分布，而且單元之間的連續(xù)性更好，因而能更好地控制數(shù)值波動現(xiàn)象，取得問題較好的數(shù)值解.本文工作可為進一步研究和發(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］.北京：機械工業(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ù)值波動研究［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.