基于最小二乘的N-S方程算子分裂有限元法

王大國,任燕緯,沈連山,THAM Leslie George

(1.大連大學材料破壞力學數(shù)值試驗研究中心,遼寧 大連 116622;2.大連大學信息工程學院,遼寧 大連 116622;3.香港大學土木工程系,香港 999077)

求解N-S方程的數(shù)值計算方法主要有有限差分法[1]、有限體積法[2]和有限元法。其中,有限元法由于適合處理復雜幾何形狀和邊界條件而成為計算力學中的優(yōu)選方法,目前已經(jīng)發(fā)展了多種有限元數(shù)值計算方法,傳統(tǒng)的Galerkin有限元法實質(zhì)上采用的是中心差分格式,隨著雷諾數(shù)的增大,N-S方程中的對流項越來越強,呈現(xiàn)強非線性特性,從而引起數(shù)值解的振蕩失真。Petrov-Galerkin法(P-G法)[3]通過將權函數(shù)取為基函數(shù)的某種修正形式使計算格式具有人工耗散能力,提高了解的精度,隨后發(fā)展了多種迎風方法,如流線迎風Petrov-Galerkin法(SUPG法)[4]等。許多其他近似方法也可以得到與P-G法相同的結果,如Galerkin最小二乘近似法(GLS法)[5]等。上述各方法都包含一個對原Galerkin法的耗散修正項,產(chǎn)生了一個包括二階空間導數(shù)的附加項。Taylor-Galerkin有限元法(T-G法)[6]時間方向上的Taylor展開先于空間上的Galerkin離散,相比P-G法不需要使用特殊的權函數(shù),也無需確定人工耗散自由參數(shù)來達到高精度。在T-G法的基礎上發(fā)展了多種有限元法,如二階全展開Euler-Taylor-Galerkin有限元法(ETG法)[7]、特征線Galerkin法[8]、基于Taylor展開的分步有限元格式[9]等?；谔卣骶€的分離算法(CBS法)[10]是結合分離算法與特征線Galerkin法的一種比較新的算法。算子分裂法是Yanenko[11]在1971年提出的經(jīng)典分裂法,該算法具有格式靈活、穩(wěn)定性好等特點。Ying[12]采用算子分裂法求解流函數(shù)渦量形式的線性N-S方程,證明了近似解的收斂性,但未給出數(shù)值解。曹志先等[13]用算子分裂法求解Burgers方程,并分析了該算法的穩(wěn)定性條件。Wang等[14]提出了CBOS方法,將N-S方程分裂成擴散項和對流項,對流項求解采用特征線Galerkin法。

本文在文獻[5,11-14]的基礎上,采用基于最小二乘的N-S方程算子分裂有限元法(least-squarebased operator-splitting,簡稱LSBOS)求解二維非定常黏性不可壓縮流體N-S方程,在每個時間層上應用算子分裂技術將N-S方程分裂成擴散項和對流項,這樣可以避免兩種不同性質(zhì)的物理過程在一起求解時計算的困難。擴散項是一個拋物線方程,時間離散采用向后差分格式,空間離散采用標準Galerkin有限元法,將其結果作為求解對流項的初值。對流項是一個雙曲線方程,時間離散仍采用向后差分格式,空間離散采用最小二乘法,并且依據(jù)最小二乘法的原理直接給出權函數(shù),避免了其他有限元法依據(jù)經(jīng)驗或人為選擇權函數(shù)的困難。分別對方腔流和后臺階流進行數(shù)值模擬,在方腔流數(shù)值模擬中,給出了速度對比曲線,分析了時間步長和空間步長對數(shù)值模擬結果的影響;在后臺階流數(shù)值模擬中,描述了不同雷諾數(shù)下的流場特征和速度對比曲線,所得計算結果與試驗結果吻合很好,說明該方法有很好的收斂性和較高的精度,認為該方法具有良好的應用前景。

1 控制方程

二維非定常黏性不可壓縮流體N-S方程的無量綱形式為

式中:(u1,u2)=(u,v),u和v分別為水平方向和垂直方向的速度;t為時間;p為壓強;f1和f2分別為水平方向和垂直方向的外力;(x1,x2)=(x,y),x,y分別為水平方向和垂直方向的坐標;Re為雷諾數(shù);ρ為流體密度;U為特征速度;L為特征長度;μ為黏性系數(shù)。

2 計算方法

2.1 算子分裂法

采用LSBOS有限元法求解式(1),計算中用算子分裂法將式(1)分成擴散項(式(2))和對流項(式(3)):擴散項

式中:ui,n+θ,uj,n+θ為擴散項(式(2))在 n+1時刻的解,同時也是對流項(式(3))在 n時刻的初值;ui,n+1,uj,n+1為對流項(式(3))在 n+1時刻的解,同時也是N-S方程(式(1))在n+1時刻的解。

2.2 對流項最小二乘法

將式(3)采用向后差分格式進行時間離散,并將ui,n+θ簡記為ui:

式中:Δt為計算時間步長。

然后用最小二乘法求解式(6),其虛功方程如下:

式中:Ω為求解區(qū)域;δui,δuj為虛位移。

3 有限元求解

3.1 擴散項有限元求解

擴散項有限元求解的詳細推導過程與文獻[14]中擴散項有限元求解的過程相同,在整個計算域合成可得總體有限元方程:

式中:Ast為速度對應的剛度矩陣;Cisr為壓力對應的剛度矩陣;Ost,Fis,Bis均為載荷矩陣,具體推導參見文獻[14];s,t為速度的總體節(jié)點號;r為壓力總體節(jié)點號。求解式(8)可得 p和ui,n+θ。

3.2 對流項有限元求解

類似擴散項的單元分析及插值,對式(7)進行分析,得到單元有限元方程:

合成總體有限元方程:

在電商購物平臺上有很大一部分不滿意的評價是來自于物流配送過程中造成的快件丟失、快件破損以及快遞慢遞等原因，而導致這些原因的主要問題在于快遞員的服務質(zhì)量不過關，快遞員的暴力丟件和私藏快件等都是導致快遞丟件、破損的主要原因。圓通承接了電商平臺很大一部分的訂單，提高服務質(zhì)量勢在必行。

式中:D st為剛度矩陣;G i s為載荷矩陣。

求解式(10)可得式(1)在 n+1時刻的速度值ui,n+1。

4 邊界條件處理及N-S方程的求解過程

4.1 邊界條件處理

采用LSBOS有限元法求解N-S方程,求解過程中只需給出擴散項的邊界條件。擴散項是一個拋物線方程,邊界條件見第5節(jié)數(shù)值驗證中的具體數(shù)值模型。

4.2 求解過程

已知n時刻的值,求解 n+1時刻的值時,N-S方程的求解過程為:①以n時刻的速度場ui,n和壓力場pn作為初值,求解擴散項(式(2)),得到 n+1時刻速度場的過渡值 ui,n+θ和壓力場pn+1;②以ui,n+θ作為初值,求解對流項(式(3)),得到 n+1時刻的速度場 ui,n+1,完成n+1時刻的求解;③轉(zhuǎn)到下一時刻,重復步驟①和②。

5 數(shù)值驗證

5.1 方腔流驗證

方腔流是流體力學中用來檢驗數(shù)值模擬可靠性的經(jīng)典算例,利用本文所建立的模型對方腔流進行數(shù)值模擬。如圖1所示,方腔無量綱尺度為1×1,內(nèi)部流體的初始速度和壓力都是零;頂邊施加無量綱速度u=1,v=0,其他三邊都是固壁,施加無滑移邊界條件u=0,v=0。坐標原點在方腔左下角,在該點置相對壓力p=0。整個計算域劃分為30×30個九節(jié)點四邊形單元,共有3721個節(jié)點。

圖1 方腔流幾何構型和邊界條件

5.1.1 不同雷諾數(shù)的影響

從圖 2可以看出,本文模型計算結果與Ghia等[17]的計算結果吻合很好。本文選擇的網(wǎng)格數(shù)為30×30,而Ghia等[17]計算的網(wǎng)格數(shù)為129×129,本文算法在采用較少單元數(shù)的情況下得到了較精確的結果,這表明基于最小二乘的N-S方程算子分裂有限元法具有較高的精度和很好的收斂性。

5.1.2 空間步長對計算精度的影響

圖2 不同雷諾數(shù)下速度沿中線的分布

圖3 Re=1000,Δt=0.0025,t=50時速度沿中線的分布

5.1.3 時間步長對計算精度的影響

圖4給出了雷諾數(shù)Re=1000、計算時間t=50、網(wǎng)格數(shù)為30×30時不同時間步長下速度沿中線分布的對比。從圖4可以看出,當 Δt=0.01時,本文計算結果與Ghia等[17]計算結果的平均相對誤差絕對值為9.152%;當Δt=0.005時誤差為5.476%;當Δt=0.0025時誤差為2.123%。

圖4 Re=1000,t=50,網(wǎng)格數(shù)為30×30時速度沿中線的分布

5.2 后臺階流驗證

后臺階流的幾何構型見圖5,圖中H為后臺階的高度,h為入口處高度(H和h均為無量綱量),臺階長度L1=10H,臺階后長度L2=30H。入流特征速度為拋物線分布,U=k(y-H)(H+h-y)(其中k為待定常數(shù),根據(jù)雷諾數(shù)的不同取值不同),固壁采用無滑移條件,出口邊界為壓力零參考面。網(wǎng)格單元為九節(jié)點四邊形,垂直方向每單位網(wǎng)格長度為H/8,水平方向每單位網(wǎng)格長度為H/2,在臺階壁面附近進行局部網(wǎng)格加密。

圖5 后臺階流的幾何構型

5.2.1 流場特征分析

計算時取 H=1,h=H,Δt=0.005。圖6給出了Re=100,500,1000下流場在 t=50時的流線。從圖6(a)可以看出,低雷諾數(shù)下流線平滑,流動穩(wěn)定,在臺階后方存在唯一回流區(qū);隨著雷諾數(shù)的增大,頂部出現(xiàn)二次回流區(qū),見圖6(b);當雷諾數(shù)繼續(xù)增大時,下壁面出現(xiàn)三次回流區(qū),見圖6(c)。文獻[18]也得到了類似的分析結果,表明本文模型能比較精確地模擬出各種雷諾數(shù)下后臺階流的流場特征。

5.2.2 速度曲線對比

為了與Armaly等[19]的試驗結果進行對比,模型的幾何參數(shù)與文獻[19]的試驗布置完全一致,取H=0.49cm,h=52H/49。圖7、圖8分別給出了雷諾數(shù)Re=100,389時水平速度u在不同水平位置處沿垂向剖面的分布。

圖6 不同雷諾數(shù)下t=50時的流線

圖7 Re=100時水平速度u在不同水平位置處沿垂向剖面的分布

圖8 Re=389時水平速度u在不同水平位置處沿垂向剖面的分布

從圖7、圖 8可以看出,當 Re=100時,本文計算結果與試驗結果吻合很好;當Re=389時,盡管個別地方計算結果與試驗結果存在一定的誤差,但整體上基本吻合。在計算過程中經(jīng)過各個斷面的流體質(zhì)量最大誤差不超過3%,基本保持不變,說明數(shù)值模型是守恒的。

6 結 語

本文提出了一種新的求解二維非定常黏性不可壓縮流體N-S方程的算法,即基于最小二乘的N-S方程算子分裂有限元法。在每一個時間層上,應用算子分裂技術將N-S方程分裂為擴散項和對流項。擴散項采用標準Galerkin有限元法,并將其結果作為求解對流項的初值;對流項采用最小二乘有限元法。將對流項與擴散項分開求解,既能考慮對流占優(yōu)特點又能兼顧方程的擴散性質(zhì),同時可以避免在整個計算域求解大規(guī)模非線性方程組。對流項采用最小二乘避免了其他有限元法選擇權函數(shù)的困難。對方腔流進行了數(shù)值模擬,計算結果與標準解吻合很好,分析了空間步長和時間步長對數(shù)值模擬結果的影響,得出空間步長對計算結果的影響較小,即在較大空間步長情況下也能取得較滿意的計算結果。同時對后臺階流也進行了數(shù)值研究,描述了不同雷諾數(shù)下的流場特征和速度對比曲線,所得計算結果與試驗結果吻合較好,表明本文提出的算法具有較好的收斂性和較高的精度。

[1]DOUGLAS J,RUSSELL T F.Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures[J].SIAM Journal on Numerical Analysis,1982,19(5):871-885.

[2]SYLVAIN B,FLORENT C,JEAN-MATE H.A finite volume method to solve the Navier-Stokes equations for in compressible flows on unstructured meshes[J].Int J Them Science,2000,39:806-825.

[3]CHRISTIE I,GRIFFITHS D F,MITCHELL A R,et al.Finite element methods for second order differential equations with significant first derivatives[J].International Journal for Numerical Methods in Engineering,1976,10:1389-1396.

[4]BROOKS A N,HUNGHES T J.Streamline upwind/preov Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[J].Computer Methods in Applied Mechanics and Engineering,1932,32(1/2/3):199-259.

[5]HUNGHES T J,FRANCA L P,HULBERT G M.A new finite element formulation for computational fluid dynamicsⅧ:the Galerkin least-squares method for advective-diffusive equations[J].Computer Methods in Applied Mechanics and Engineering,1989,73(2):173-189.

[6]DONEA J.A Taylor-Galerkin method for convective transport problems[J].International Journal for Numerical Methods in Engineering,1984,20:101-119.

[7]王小華,樊洪明,何鐘怡.后臺階流動的數(shù)值模擬[J].計算力學學報,2003,20(3):361-365.

[8]MORTON K W.Generalised Galerkin methods for hyperbolic problems[J].Computer Methods in Applied Mechanics and Engineering,1985,52(1/2/3):847-871.

[9]JIANG C B,KAWAHARA M,KASHIYMA K.A Taylor-Galerkin-based finite element method for turbulent flows[J].Fluid Dynamic Research,1992,9:165-178.

[10]ZIENKINEWICZ O C,TAYLOR R L.The finite element method:fluid dynamics[M].5th edition.Oxford:Butterworth Heinemann,2000.

[11]YANENKO N N.The method of fractional steps[M].Berlin:Springer-Verlag,1971.

[12]YING Long-an.Viscosity splitting method in bounded domains[J].Science in China:Series A,1989,32(8):908-921.

[13]曹志先,魏良琰.用算子分裂法解Burgers方程[J].武漢水利水電學院學報,1991,24(2):193-201.

[14]WANG Da-guo,WANG Hai-jiao,XIONG Ju-hua,et al.Characteristic-based operator-splitting finite element method for Navier-Stokes equations[J].Science in China:Series E,2011,54(8):2157-2166.

[15]李慶揚,王能超,易大義.數(shù)值分析[M].4版.北京:清華大學出版社,2001.

[16]章本照.流體力學中的有限元方法[M].北京:機械工業(yè)出版社,1984.

[17]GHIA U,GHIA K N,SHIN C T.High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[J].Journalof Computational Physics,1982,48:387-411.

[18]王兵,張會強,虞建豐,等.后臺階分離流動中大渦結構演變的數(shù)值模擬[J].力學季刊,2003,24(2):166-173.

[19]ARMALY B F,DURST F,PEREIRA J C F,et al.Experimentaland theoretical investigation of backward-facing step flow[J].Fluid Mech,1983,127:473-496.