改進的基于特征線的N-S方程算子分裂有限元法

水慶象，王大國，沈連山

（1西南科技大學環(huán)境與資源學院，四川綿陽621010；2大連大學信息工程學院，遼寧大連116622）

改進的基于特征線的N-S方程算子分裂有限元法

水慶象1，王大國1，沈連山2

（1西南科技大學環(huán)境與資源學院，四川綿陽621010；2大連大學信息工程學院，遼寧大連116622）

基于特征線的Navier-Stokes（N-S）方程算子分裂有限元法（CBOS有限元法）的核心是在每一個時間層上，采用算子分裂法將N-S方程的對流項和擴散項分開求解；對流項的求解過程借鑒了簡單顯式特征線時間離散，顯式求解。該文在原CBOS有限元法基礎(chǔ)上推導(dǎo)了一種更加精確的對流項顯式離散方法。通過自編程序?qū)Ψ角涣鬟M行數(shù)值模擬，表明該算法具有更高的計算精度。低雷諾數(shù)下圓柱繞流計算所得的阻力系數(shù)、升力系數(shù)、斯特勞哈數(shù)等與已有數(shù)據(jù)較為接近，表明該文中算法能較準確地模擬圓柱繞流的流場特性；最后，文中分析了Re= 200時圓柱繞流在一個周期內(nèi)所受升力變化與對應(yīng)流場中壓力和流線演化的關(guān)系。

N-S方程；CBOS有限元法；方腔流；圓柱繞流

0 引 言

N-S方程的動量方程中存在對流項導(dǎo)致方程不具備自伴隨性，尤其當雷諾數(shù)較大使得對流占優(yōu)時，方程會呈現(xiàn)出強非線性特性，采用標準Galerkin有限元法求解將導(dǎo)致數(shù)值解的振蕩。針對對流項占優(yōu)時導(dǎo)致數(shù)值解振蕩的問題，自20世紀80年代以來，已有多種有別于標準Galerkin有限元的新算法問世。目前應(yīng)用較為廣泛的有Petrov-Galerkin法[1]、Taylor-Galerkin法[2]、Galerkin最小二乘法[3]和特征線Galerkin法[4]等。其中，特征線Galerkin法的基本思想是將N-S方程的動量方程沿特征線進行時間離散，以理性形式引入穩(wěn)定項，具有比較明確的物理意義。1995年，Zienkiewicz和Codina[5]將Taylor展開法引入到特征線Galerkin法中并將其與分離算法相結(jié)合，提出了基于特征線的分離算法（CBS），該方法可以直接由N-S方程推導(dǎo)出合理的平衡耗散項，從而可以獲得穩(wěn)定的數(shù)值解。王大國等人[6]結(jié)合了CBS算法和算子分裂法[7]的優(yōu)點，提出了基于特征線的算子分裂有限元法（CBOS有限元法）求解二維非定常不可壓N-S方程，但沒有對圓柱繞流進行數(shù)值模擬。

CBOS有限元法根據(jù)N-S方程的特性，采用算子分裂法將N-S方程分裂成對流項和擴散項，這種分裂方法既能考慮方程的擴散性質(zhì)，又能突出對流占優(yōu)的特性。對流項的求解借鑒了CBS法的簡單顯式特征線時間離散，它的思想包括局部Taylor展開。然而原CBOS有限元法中為了獲得對流項的完全顯示的計算格式做了一定的近似而導(dǎo)致算法的精度降低，本文在原CBOS有限元法基礎(chǔ)上推導(dǎo)了一種更加精確的對流項顯示離散方法。通過自編程序?qū)Ψ角涣鬟M行數(shù)值模擬，表明該算法具有更高的計算精度。低雷諾數(shù)下圓柱繞流計算所得的阻力系數(shù)、升力系數(shù)、斯特勞哈數(shù)等與已有數(shù)據(jù)較為接近，表明該算法能較準確地模擬圓柱繞流的流場特性；最后本文分析了Re=200時圓柱繞流在一個周期內(nèi)所受升力變化與對應(yīng)流場中壓力和流線演化的關(guān)系。

1 控制方程

二維非定常粘性不可壓縮流動（忽略能量損失）可由連續(xù)方程和N-S方程（統(tǒng)稱N-S方程）控制，其無量綱形式表示為

式中： （ui,ui）=（u,v），u為水平方向速度，v為垂直方向速度，p為壓強，t為時間，Re為雷諾數(shù)，f1為水平方向外力，f2為垂直方向外力，（xi,xi）=（x,y）。

2 計算方法

2.1 算子分裂法

本文采用算子分裂法將控制方程（1）和（2）分成擴散項

和對流項

2.2 對流項的特征線法顯式時間離散

中的沿流線的穩(wěn)定擴散項是不同的。由于穩(wěn)定擴散項在高雷諾數(shù)或可壓縮流動計算中是抑制振蕩尤其是壓力振蕩的基礎(chǔ)[8]，所以采用改進的CBOS有限元法對流動問題進行數(shù)值模擬將會在計算域內(nèi)得到更為穩(wěn)定的壓力分布。

3 有限元求解

3.1 擴散項有限元求解

擴散項中出現(xiàn)了速度和壓力變量，采用混合有限元求解，為滿足LBB條件取單元為九節(jié)點四邊形單元，速度插值采用九節(jié)點，壓力插值采用四節(jié)點（單元的角節(jié)點），應(yīng)用標準Galerkin法對上式進行有限元空間離散，可得

式中：Nα和Nβ為二次插值函數(shù)，Ml為線性插值函數(shù)，α=1,2,…,m，β=1,2,…,m，l=1,2,…,h，m=9為插值單元內(nèi)速度節(jié)點個數(shù)，h=4為插值單元內(nèi)壓力節(jié)點個數(shù)，δij為置換算子，δijuiβ=ujβ。

3.2 對流項有限元求解

同擴散項處理方法一樣，由標準Galerkin法建立（20）式的弱形式（省略上標）

對上式最后一項分部積分（忽略邊界項的影響）得

式中：Nγ和Nη為二次插值函數(shù)，γ，η=1,2,…,m。

4 邊界條件處理及求解過程

本文采用改進的CBOS有限元法求解N-S方程，求解過程中只需給出擴散項的邊界條件，擴散項是一個拋物型方程，邊界條件見后面具體的數(shù)值模型。本文算法的求解過程如下：

（3）轉(zhuǎn)到下一時刻，重復(fù)步驟（1）、（2）。

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é)點四邊形單元，共有3 721個節(jié)點。

圖1 方腔流幾何構(gòu)型和邊界條件Fig.1 Cavity flow configuration and boundary conditions

5.2 不同雷諾數(shù)下的計算情況

圖2給出了不同雷諾數(shù)下水平速度u沿垂直中線、垂向速度v沿水平中線的分布結(jié)果。圖中實線表示本文計算結(jié)果，虛線表示原CBOS有限元法[6]的計算結(jié)果，點表示Ghia等人[9]的計算結(jié)果。圖2（a）給出Re=1 000、Δt=0.005、t=28時的結(jié)果；圖2（b）給出Re=5 000、Δt=0.001、t=110時的結(jié)果；圖2（c）給出Re=5 000、Δt=0.01、t=110時的結(jié)果。從圖2（a）和圖2（b）可以看出，改進前后的CBOS有限元法計算結(jié)果比較接近，且與Ghia等人[9]的計算結(jié)果符合較好。但在較大時間步長情況下，本文算法計算結(jié)果比原CBOS有限元法計算結(jié)果具有更高的計算精度，這從圖2（b）和圖2（c）可以看出。

圖2 速度沿中線分布（實線：本文計算結(jié)果；虛線：原CBOS有限元法[6]的計算結(jié)果；點：Ghia等人[9]的計算結(jié)果）Fig.2 Velocity along lines through geometric center

5.3 改進的CBOS法對壓力的影響

圖3 壓力沿水平中線分布（實線：本文計算結(jié)果；虛線：原CBOS有限元法[6]的計算結(jié)果；點：CBS法[8]的計算結(jié)果）Fig.3 Pressure along horizontal center line

圖3給出了方腔流壓力p沿垂直中面分布結(jié)果。圖中實線表示本文計算結(jié)果，虛線表示原CBOS有限元法[6]計算結(jié)果，點表示CBS算法（半隱式格式）[8]的計算結(jié)果。圖3（a）給出Re=1 000、Δt=0.005、t=28時的結(jié)果；圖3（b）給出Re=5 000、Δt=0.001、t=110時的結(jié)果。由圖3可以看出，本文算法計算的壓力更加接近CBS算法的計算結(jié)果，且雷諾數(shù)越大，壓力的改進越大。究其原因，是由于改進前后的CBOS有限元法引入了不同的沿流線的穩(wěn)定擴散項，相關(guān)研究[8]表明穩(wěn)定擴散項在高雷諾數(shù)或可壓縮問題中相當重要，它可以有效地控制壓力的振蕩。

6 圓柱繞流工程應(yīng)用

圓柱繞流是一個經(jīng)典的流體力學問題，是最具代表性的鈍體繞流問題，本文對Re=100和Re=200的圓柱繞流進行了數(shù)值研究。

6.1 計算區(qū)域、網(wǎng)格劃分和邊界條件

取圓柱直徑D=0.1，計算區(qū)域尺度在主流方向上取為30D，其中圓柱上游分配10D，橫向尺寸為18D，入口處指定無因次水平速度為u=1，橫向速度v=0；側(cè)壁采用可滑移邊界條件；圓柱表面為不可滑移邊界條件；出口處為自由出流邊界，并指定相對壓力為p=0。整個計算域劃分為2 964個九節(jié)點四邊形單元，共有12 114個節(jié)點，由于圓柱附近流場比較復(fù)雜，含有渦旋的生成與脫落等現(xiàn)象，需要布置較密的網(wǎng)格；而在遠離圓柱的區(qū)域則可布置較疏的網(wǎng)格以節(jié)省計算量，具體的網(wǎng)格劃分如圖4所示。

圖4 計算區(qū)域網(wǎng)格劃分Fig.4 Sketch of computation grid

6.2 阻力系數(shù)、升力系數(shù)以及斯特勞哈數(shù)的比較

（30）式中f0表示渦的自然脫落頻率。

表1 Re=100和Re=200時本文計算結(jié)果與其他文獻數(shù)據(jù)的比較Tab.1 Comparison between calculated results with data provided by other references at Re=100 and Re=200

續(xù)表1

從表1中可以看出，本文的計算結(jié)果與相關(guān)參考文獻的數(shù)值結(jié)果符合較好。圖5給出了圓柱繞流的升力系數(shù)和阻力系數(shù)時程曲線，從圖中可以發(fā)現(xiàn)阻力脈動頻率為升力的2倍。

圖5 升力和阻力時程曲線（實線：阻力系數(shù)，虛線：升力系數(shù)）Fig.5 Variations of lift coefficient and drag coefficient with time

6.3 升力的周期變化與流場的演化

圖6給出了Re=200時的固定圓柱繞流一個渦旋脫落周期內(nèi)升力系數(shù)的時程曲線，圖7給出了Re=200時的固定圓柱繞流在大致一個渦旋脫落周期內(nèi)的壓力和流線的演化過程。圖6中的●點與圖7中各時刻流場圖一一對應(yīng)。

圖6 一個渦旋脫落周期內(nèi)升力時程曲線Fig.6 Time variations of lift coefficient in a cycle

圖7（a）給出了當t=8.96時流場的壓力云圖和流線圖，可以看出，圓柱下方壓力偏低，上部壓力相對較高，這在圓柱的背流面表現(xiàn)得更為明顯，因此圓柱在總體上受到負向壓力作用。在同一時刻的流線圖中，圓柱背流面的下部存在一個貼體主渦旋區(qū)，該渦旋區(qū)同時對應(yīng)流場中的主低壓區(qū)；而在圓柱表面背流面上部出現(xiàn)流動分離趨勢，在分離點后方的壓力為高壓區(qū)，該高壓區(qū)與圓柱正上方表面的低壓分布正好形成逆壓梯度，逆壓梯度驅(qū)動近壁面流體倒流產(chǎn)生分離。

圖7 Re=200時一個周期內(nèi)的壓力和流線演化Fig.7 The development of pressure and streamline during a cycle at Re=200

隨著時間的發(fā)展，在t=9.02時可以發(fā)現(xiàn)，圓柱附近的貼體主渦開始脫離壁面，主低壓區(qū)也開始遠離圓柱，導(dǎo)致圓柱背流面下方壁面上的壓力升高；與此同時，由流線圖可知在原分離點位置形成較明顯的分離渦旋，壓力也相應(yīng)減小。上述壓力和流線的演化導(dǎo)致圓柱在總體上所受的負向壓力減小，這與升力時程曲線由a到b的變化情況相符。

在t=9.10時刻對應(yīng)于圓柱所受升力基本為零的狀態(tài)。在時間為9.02～9.10階段伴隨著圓柱背流面上部渦旋區(qū)擴大，臨近區(qū)域的壓力持續(xù)降低，圓柱尾跡區(qū)的主漩渦則向下游進一步推移，相應(yīng)的低壓區(qū)也更進一步地遠離壁面，圓柱背流面下部附近的壓力繼續(xù)增大。圓柱背流面上部壓力減小與下部壓力增大趨勢的持續(xù)發(fā)展將導(dǎo)致圓柱所受負向升力持續(xù)減小直至出現(xiàn)為零，如升力時程曲線上的c點。之后，渦流場的這一演化持續(xù)進行，當圓柱背流面下部壓力在總體上高于背流面上部壓力時，圓柱受到正向升力的作用，如圖7（d）所示在t=9.16時的流場圖。

當t=9.22時作用于圓柱上的正向升力達到最大。圓柱背流面下部的主渦旋已經(jīng)基本消散，此時流場中的貼體主渦旋區(qū)位于圓柱背流面上方，并于此形成低壓區(qū)；在圓柱背流面下部近壁面則出現(xiàn)了新的流動分離，該區(qū)域壓力達到半個周期來的最大，對比圖7（a）和圖7（e）可以發(fā)現(xiàn)兩幅圖的流線的空間分布恰好相反。此后圓柱受力狀態(tài)進入正向升力減小階段，流線以及壓力的演化規(guī)律類似于前半個周期的情況，分別將圖7（f）與圖7（b）、圖7（g）與圖7（c）、圖7（h）與圖7（d)對比可以發(fā)現(xiàn)以上各對流線和壓力分布都比較相似，只是空間分布相反。這種對稱演化過程也充分體現(xiàn)了升力的周期性變化特征。

當t=9.50時，圓柱所受升力又重新達到負向最大，相應(yīng)的渦旋分離、脫落與壓力分布也完成了一個周期的演化過程。

7 結(jié) 論

（1）在原CBOS有限元法[6]基礎(chǔ)上，推導(dǎo)了一種更加精確的對流項顯式離散格式。通過方腔流的數(shù)值模擬，發(fā)現(xiàn)當時間步長較大時，本文算法計算結(jié)果比原CBOS有限元法計算結(jié)果更接近Ghia等人[9]的計算結(jié)果。同時，相對原CBOS有限元法，本文算法壓力的計算結(jié)果與CBS法[8]壓力的計算結(jié)果更加吻合，尤其是在高雷諾數(shù)的情況下壓力的改變更加明顯，究其原因，主要是由于本文算法引入了不同的沿流線的穩(wěn)定擴散項，穩(wěn)定擴散項在高雷諾數(shù)或可壓縮問題中相當重要，它可以有效地控制壓力的振蕩。

（2）通過對低雷諾數(shù)的圓柱繞流進行數(shù)值模擬，計算所得的阻力系數(shù)、升力系數(shù)、斯特勞哈數(shù)等與參考文獻中的值較為接近，表明本文提出的算法能較準確地模擬圓柱繞流的流場特性。通過對Re= 200的圓柱繞流在一個周期內(nèi)所受升力的變化及其相應(yīng)的流場中的壓力和流線演化分析可知：當圓柱所受升力達到最大（負向最大或者正向最大）時，在圓柱后方尾流區(qū)內(nèi)的一側(cè)存在明顯的貼體渦旋，同時在圓柱表面尾流區(qū)內(nèi)的另一側(cè)由于逆壓梯度的驅(qū)動開始出現(xiàn)流動分離現(xiàn)象。

[1]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(6):1389-1396.

[2]Donea J.A Taylor-Galerkin method for convective transport problems[J].International Journal for Numerical Methods in Engineering,1984,20(1):101-119.

[3]Hughes T J R,Franca L P,Hulbert G M.A new finite element formulation for computational fluid dynamics:VIII.The Galerkin/least-squares method for advective diffusive equations[J].Computer Methods in Applied Mechanics and Engineering,1989,73(2):173-189.

[4]Pironneau O,Liou J,Tezduyar T.Characteristic-Galerkin and Galerkin/least-squares space-time formulations for the advection-diffusion equation with time-dependent domains[J].Computer Methods in Applied Mechanics and Engineering, 1992,100(1):117-141.

[5]Zienkiewicz O C,Codina R.A general algorithm for compressible and incompressible-flow,I:The split characteristic based scheme[J].International Journal for Numerical Methods in Fluids,1995,20(8-9):869-885.

[6]Wang Daguo,Wang Haijiao,Xiong Juhua,et al.Characteristic-based Operator-splitting Finite Element Method for Navier-Stokes equations[J].Science in China(Series E),2011,54(8):2157-2166.

[7]Glowinski R,Pironneau O.Finite element method for Navier-Stokes equations[J].Annual Review of Fluid Mechanics, 1992,24(1):167-204.

[8]Zienkiewicz O C,Taylor R L.有限元方法（第五版）第三卷流體動力學[M].符 松,劉揚揚譯.北京:清華大學出版社, 2008.

[9]Ghia U,Ghia K N,Shin C T.High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[J].Journal of Computational Physics,1982,48(3):387-411.

[10]Liu Z,Zheng X,Sung C H.Preconditioned multi-grid methods for unsteady incompressible flows[J].Journal of Computational Physics,2000,160(1):151-178.

[11]Wu G X,Hu Z Z.Numerical simulation of viscous flow around unrestrained cylinders[J].Journal of Fluids Structures, 2006,22(3):371-390.

[12]Xu S,Wang Z J.An immersed interface method for simulating the interaction of a fluid with moving boundaries[J].Journal of Computational Physics,2006,216(2):454-493.

[13]Mahfouz F M,Badr H M.Flow structure in the wake of a rotationally oscillation cylinder[J].Journal of Fluids Engineering, 2000,122(2):290-301.

[14]魏志理,孫德軍,尹協(xié)遠.圓柱尾跡流場中橫向振蕩翼型繞流的數(shù)值模擬[J].水動力學研究與進展,2006,21(3):299-308. Wei Zhili,Sun Dejun,Yin Xieyuan.A numerical simulation of flow around a transversely oscillating hydrofoil in the wake of a circular cylinder[J].Journal of Hydrodynamics,2006,21(3):229-308.

The pivotal ideas of characteristic-based operator-splitting(CBOS)finite element method for Navier-Stokes equations are that the equations are split into the diffusive part and the convective part by adopting the operator-splitting algorithm in each time step.The convective part can be discretized using the simple explicit characteristic temporal discretization and solved explicitly.On the basis of CBOS finite element method,an exact explicit discrete method is developed.The improved method has higher calculation accuracy through numerical simulation of lid-driven cavity flow.Furthermore,the method is used to simulate the incompressible viscous flow around cylinder with low Reynolds number.The numerical results agree well with other numerical results,which proves that it can exactly simulate the characteristics of flow around cylinder.Finally,the relationship between a cycle of lift coefficient and the corresponding development of pressure and streamline is analyzed at Re=200.

N-S equations;CBOS finite element method;lid-driven cavity flow;flow around cylinder

O35

：Adoi:10.3969/j.issn.1007-7294.2016.04.001

1007-7294（2016）04-0381-12

2015-11-15

國家自然科學基金（51349011，41072235）

水慶象（1988-），男，碩士，講師；王大國（1975-），男，博士，教授，通訊作者，E-mail:dan_wangguo@163.com。

Improved characteristic-based operator-splitting finite element for Navier-Stokes equations SHUI Qing-xiang1,WANG Da-guo1,SHEN Lian-shan2

(1 School of Environment and Resource,Southwest University of Science and Technology,Mianyang 621010,China; 2 Department of Information Engineering,Dalian University,Dalian 116622,China)