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

王大國(guó),任燕緯,沈連山,THAM Leslie George

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

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

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

1 控制方程

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

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

2 計(jì)算方法

2.1 算子分裂法

采用LSBOS有限元法求解式(1),計(jì)算中用算子分裂法將式(1)分成擴(kuò)散項(xiàng)(式(2))和對(duì)流項(xiàng)(式(3)):擴(kuò)散項(xiàng)

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

2.2 對(duì)流項(xiàng)最小二乘法

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

式中:Δt為計(jì)算時(shí)間步長(zhǎng)。

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

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

3 有限元求解

3.1 擴(kuò)散項(xiàng)有限元求解

擴(kuò)散項(xiàng)有限元求解的詳細(xì)推導(dǎo)過(guò)程與文獻(xiàn)[14]中擴(kuò)散項(xiàng)有限元求解的過(guò)程相同,在整個(gè)計(jì)算域合成可得總體有限元方程:

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

3.2 對(duì)流項(xiàng)有限元求解

類(lèi)似擴(kuò)散項(xiàng)的單元分析及插值,對(duì)式(7)進(jìn)行分析,得到單元有限元方程:

合成總體有限元方程:

在電商購(gòu)物平臺(tái)上有很大一部分不滿意的評(píng)價(jià)是來(lái)自于物流配送過(guò)程中造成的快件丟失、快件破損以及快遞慢遞等原因，而導(dǎo)致這些原因的主要問(wèn)題在于快遞員的服務(wù)質(zhì)量不過(guò)關(guān)，快遞員的暴力丟件和私藏快件等都是導(dǎo)致快遞丟件、破損的主要原因。圓通承接了電商平臺(tái)很大一部分的訂單，提高服務(wù)質(zhì)量勢(shì)在必行。

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

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

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

4.1 邊界條件處理

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

4.2 求解過(guò)程

已知n時(shí)刻的值,求解 n+1時(shí)刻的值時(shí),N-S方程的求解過(guò)程為:①以n時(shí)刻的速度場(chǎng)ui,n和壓力場(chǎng)pn作為初值,求解擴(kuò)散項(xiàng)(式(2)),得到 n+1時(shí)刻速度場(chǎng)的過(guò)渡值 ui,n+θ和壓力場(chǎng)pn+1;②以u(píng)i,n+θ作為初值,求解對(duì)流項(xiàng)(式(3)),得到 n+1時(shí)刻的速度場(chǎng) ui,n+1,完成n+1時(shí)刻的求解;③轉(zhuǎn)到下一時(shí)刻,重復(fù)步驟①和②。

5 數(shù)值驗(yàn)證

5.1 方腔流驗(yàn)證

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

圖1 方腔流幾何構(gòu)型和邊界條件

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

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

5.1.2 空間步長(zhǎng)對(duì)計(jì)算精度的影響

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

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

5.1.3 時(shí)間步長(zhǎng)對(duì)計(jì)算精度的影響

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

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

5.2 后臺(tái)階流驗(yàn)證

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

圖5 后臺(tái)階流的幾何構(gòu)型

5.2.1 流場(chǎng)特征分析

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

5.2.2 速度曲線對(duì)比

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

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

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

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

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

6 結(jié) 語(yǔ)

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

[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]王小華,樊洪明,何鐘怡.后臺(tái)階流動(dòng)的數(shù)值模擬[J].計(jì)算力學(xué)學(xué)報(bào),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].武漢水利水電學(xué)院學(xué)報(bào),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]李慶揚(yáng),王能超,易大義.數(shù)值分析[M].4版.北京:清華大學(xué)出版社,2001.

[16]章本照.流體力學(xué)中的有限元方法[M].北京:機(jī)械工業(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]王兵,張會(huì)強(qiáng),虞建豐,等.后臺(tái)階分離流動(dòng)中大渦結(jié)構(gòu)演變的數(shù)值模擬[J].力學(xué)季刊,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.