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

水慶象，王大國，沈連山

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

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

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

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

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

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

0 引 言

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

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

1 控制方程

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

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

2 計(jì)算方法

2.1 算子分裂法

本文采用算子分裂法將控制方程（1）和（2）分成擴(kuò)散項(xiàng)

和對(duì)流項(xiàng)

2.2 對(duì)流項(xiàng)的特征線法顯式時(shí)間離散

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

3 有限元求解

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

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

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

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

同擴(kuò)散項(xiàng)處理方法一樣，由標(biāo)準(zhǔn)Galerkin法建立（20）式的弱形式（省略上標(biāo)）

對(duì)上式最后一項(xiàng)分部積分（忽略邊界項(xiàng)的影響）得

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

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

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

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

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

5.1 方腔流驗(yàn)證

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

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

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

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

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

5.3 改進(jìn)的CBOS法對(duì)壓力的影響

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

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

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

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

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

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

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

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

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

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

續(xù)表1

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

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

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

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

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

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

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

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

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

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

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

7 結(jié) 論

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

（2）通過對(duì)低雷諾數(shù)的圓柱繞流進(jìn)行數(shù)值模擬，計(jì)算所得的阻力系數(shù)、升力系數(shù)、斯特勞哈數(shù)等與參考文獻(xiàn)中的值較為接近，表明本文提出的算法能較準(zhǔn)確地模擬圓柱繞流的流場特性。通過對(duì)Re= 200的圓柱繞流在一個(gè)周期內(nèi)所受升力的變化及其相應(yīng)的流場中的壓力和流線演化分析可知：當(dāng)圓柱所受升力達(dá)到最大（負(fù)向最大或者正向最大）時(shí)，在圓柱后方尾流區(qū)內(nèi)的一側(cè)存在明顯的貼體渦旋，同時(shí)在圓柱表面尾流區(qū)內(nèi)的另一側(cè)由于逆壓梯度的驅(qū)動(dòng)開始出現(xiàn)流動(dòng)分離現(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.有限元方法（第五版）第三卷流體動(dòng)力學(xué)[M].符 松,劉揚(yáng)揚(yáng)譯.北京:清華大學(xué)出版社, 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é)遠(yuǎn).圓柱尾跡流場中橫向振蕩翼型繞流的數(shù)值模擬[J].水動(dòng)力學(xué)研究與進(jìn)展,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

國家自然科學(xué)基金（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)