帶有障礙物潰壩流動(dòng)問題的有限元數(shù)值模擬研究

高普陽, 楊 揚(yáng)

(1.長(zhǎng)安大學(xué) 理學(xué)院，西安 710064；2.西北工業(yè)大學(xué) 航空學(xué)院，西安 710072)

1 引 言

潰壩流動(dòng)問題出現(xiàn)在海洋工程、水動(dòng)力學(xué)和泥石流滑坡等諸多領(lǐng)域，該過程會(huì)遇到運(yùn)動(dòng)界面的翻轉(zhuǎn)、飛濺及破碎等大變形情況，一直都是學(xué)者們關(guān)注的熱點(diǎn)問題[1-5]。

潰壩坍塌流動(dòng)過程會(huì)涉及到牛頓或者非牛頓流體和空氣之間的相互作用，可以將其視為兩相流動(dòng)問題。王吉飛等[6]采用Level Set方法追蹤運(yùn)動(dòng)界面，并通過SUPG格式數(shù)值求解Level Set方程。對(duì)于Navier-Stokes方程的求解，則在FEM中加入了SUPG和PSPG的穩(wěn)定化方法。Issakhov等[7]依據(jù)流體體積函數(shù)法捕捉運(yùn)動(dòng)界面，以保證流體的質(zhì)量守恒性，并考慮了湍流模型。以上學(xué)者們的研究?jī)H限于牛頓流體。

自然界中諸多流動(dòng)現(xiàn)象也常涉及到非牛頓流體，Cremonesi等[8]研究了賓漢非牛頓流體的潰壩流動(dòng)過程，并模擬了泥石流現(xiàn)象。采用顯式粒子有限元方法數(shù)值求解了拉格朗日形式的Navier-Stokes方程，但是并沒有考慮下游帶有障礙物的非牛頓流體潰壩流動(dòng)情況。目前對(duì)于非牛頓流體潰壩流動(dòng)問題，特別是下游帶有障礙物情形的研究鮮有報(bào)道。

Level Set方法可以簡(jiǎn)單和有效地追蹤運(yùn)動(dòng)界面，但其最大的缺陷就是容易引起流體的質(zhì)量誤差。Yu等[9]采用CLSVOF方法模擬了帶有障礙物的潰壩流動(dòng)，通過Level Set方法追蹤運(yùn)動(dòng)界面，同時(shí)利用VOF方法對(duì)流體質(zhì)量進(jìn)行追蹤，以保證良好的質(zhì)量守恒性，但這時(shí)需要另外求解VOF輸運(yùn)方程。Level Set方程屬于對(duì)流型方程，直接采用標(biāo)準(zhǔn)有限元方法求解時(shí)，容易引起數(shù)值振蕩現(xiàn)象。Touré等[10]將Level Set方程轉(zhuǎn)化為一個(gè)帶有參數(shù)的對(duì)流-擴(kuò)散方程，并采用穩(wěn)定化SUPG方法進(jìn)行空間離散，但這會(huì)產(chǎn)生一些額外的數(shù)值誤差。間斷有限元方法在求解雙曲型方程時(shí)具有諸多優(yōu)點(diǎn)[1,5]，因此本文采用了間斷有限元方法對(duì)Level Set 方程進(jìn)行求解，以保證計(jì)算過程的穩(wěn)定性和數(shù)值精確度。

本文主要目的是基于兩相流動(dòng)模型，在有限元算法框架下，數(shù)值模擬研究帶有障礙物的牛頓和非牛頓潰壩流動(dòng)問題。對(duì)于自由界面的追蹤，采用簡(jiǎn)單和有效的Level Set方法，并引入了一個(gè)質(zhì)量矯正技術(shù)，以保證長(zhǎng)時(shí)間計(jì)算以及界面大變形過程中的流體質(zhì)量守恒性。另外，本文在控制方程中加入了冪律型本構(gòu)模型，用以描述非牛頓流體的流動(dòng)特性。通過基于分裂格式的SUPG穩(wěn)定化有限元算法求解統(tǒng)一的Navier-Stokes方程，并采用間斷有限元方法求解雙曲型Level Set及其重新初始化方程。

2 數(shù)學(xué)模型

2.1 Level Set方程

隨著時(shí)間發(fā)展，Level Set函數(shù)值會(huì)根據(jù)速度進(jìn)行更新。對(duì)于不可壓縮流體，守恒形式的Level Set方程為

(1)

為了保證計(jì)算過程中Level Set函數(shù)的符號(hào)距離函數(shù)特性，需要求解如下的Level Set重新初始化方程[12]，即

(2)

2.2 統(tǒng)一的Navier-Stokes方程

潰壩流動(dòng)中會(huì)涉及到兩種不同的流體，為了求解方便，本文將整個(gè)流場(chǎng)控制方程寫成統(tǒng)一形式為

(3)

?u/?+(u·)u=-p+(1)·

{μ[u+(u)T]}+(1)f

(4)

(5)

2.3 冪律型本構(gòu)模型

根據(jù)流變學(xué)理論，非牛頓流體的粘度會(huì)受到剪切速率變化的影響。本文采用冪律型本構(gòu)模型[13]描述剪切應(yīng)力和形變速率之間的非線性關(guān)系，應(yīng)變率張量定義為d=[u+(u)T]/2，量級(jí)γ定義為對(duì)于冪律型本構(gòu)模型，粘度與應(yīng)變率張量的關(guān)系為

μl=μlγn - 1

(6)

式中n為冪律指數(shù)。n=1時(shí)，該模型表示牛頓流體；n<1時(shí)，該模型表示假塑性非牛頓流體；n>1時(shí)，該模型表示脹塑性非牛頓流體。

3 統(tǒng)一Navier-Stokes方程的數(shù)值求解

3.1 分裂格式及時(shí)間離散

對(duì)于統(tǒng)一Navier-Stokes，本文采用Euler向前格式進(jìn)行時(shí)間離散，對(duì)于非線性對(duì)流項(xiàng)和表面張力項(xiàng)進(jìn)行顯式處理，并隱式處理壓力項(xiàng)和粘性項(xiàng)。這樣，可以得到半離散格式

(7)

本文采用分裂格式，首先忽略方程(7)右端的壓力項(xiàng)和粘性項(xiàng)，得到雙曲型子方程

(8)

(9)

(10)

根據(jù)pn + 1，就可以計(jì)算中間速度

(11)

最后，可以得到關(guān)于速度的Helmholtz方程為

(1/ρ)·{μ[un + 1+(un + 1)T]}+(1/Δt)un + 1=

(12)

將分裂所得的子方程 (8，9，12)疊加，就可以得到半離散方程(7)。

3.2 空間離散

本文采用Ωh表示對(duì)求解區(qū)域Ω進(jìn)行正則三角剖分，Ωi(i=1,2,…,N)代表網(wǎng)格上的小三角形，有限元空間可以定義為

Ck(Ωh)={v∈C0(Ω)?Ωi∈Ωh,v|Ωi∈Pk(Ωi)}

(13)

式中Pk(Ωi)為Ωi上不低于k次的多項(xiàng)式。

對(duì)于雙曲型子方程(8)，采用SUPG方法[14]進(jìn)行空間離散，對(duì)方程兩邊同時(shí)乘以試探函數(shù)。u速度分量方程的具體形式為

(14)

另外的幾個(gè)子方程均屬于橢圓型方程，因此采用標(biāo)準(zhǔn)FEM進(jìn)行空間離散。方程(10)的空間離散形式為

(15)

類似的，可以得到方程方程(11，12)的空間離散格式。

4 Level Set方程的數(shù)值求解

在求解Navier-Stokes方程得到速度后，就可以進(jìn)一步求解Level Set方程以更新運(yùn)動(dòng)界面的位置和形態(tài)。對(duì)于Level Set方程(1)的時(shí)間離散，采用向前Euler格式

(16)

由于Level Set方程(1)屬于雙曲型方程，因而本文采用穩(wěn)定和有效的間斷有限元方法進(jìn)行空間離散。首先在之前網(wǎng)格三角剖分的基礎(chǔ)上，定義間斷有限元空間為

Dk(Ωh)={l∈L2(Ω)∶?Ωi∈Ωh,l|Ωi∈Pk(Ωi)}

(17)

另外，Lax-Friedrichs數(shù)值通量f*可以定義為

(18)

在方程(18)兩端同時(shí)乘以試探函數(shù)lh，并關(guān)于對(duì)流項(xiàng)做兩次分部積分，可得

(19)

(20)

類似的，可以得到Level Set重新初始化方程的全離散格式。計(jì)算過程中時(shí)間步長(zhǎng)的選取依據(jù)為

(21)

5 質(zhì)量矯正技術(shù)

為了保證流體的質(zhì)量守恒性，在整個(gè)求解過程中加入了一個(gè)簡(jiǎn)單的質(zhì)量矯正技術(shù)。如圖1所示，假設(shè)現(xiàn)在計(jì)算得到的運(yùn)動(dòng)界面為G，這時(shí)得到的液體質(zhì)量為S，而真實(shí)的流體質(zhì)量為Se。S與Se之間可能會(huì)存在一個(gè)比較小的誤差(假設(shè)S

(22)

(23)

圖1 質(zhì)量矯正技術(shù)

6 數(shù)值算例

根據(jù)前面發(fā)展的耦合算法，數(shù)值模擬研究牛頓、非牛頓潰壩流動(dòng)問題以及下游帶有障礙物的情形，所有的數(shù)值模擬均是基于Freefem++[15]開源計(jì)算平臺(tái)進(jìn)行。

6.1 單渦剪切流動(dòng)

首先研究單渦剪切流動(dòng)問題，用來驗(yàn)證數(shù)值方法較好的質(zhì)量守恒性。初始時(shí)刻圓心所在位置為(0.5,0.75)，圓的半徑為0.15，整個(gè)計(jì)算區(qū)域?yàn)閇0,1]×[0,1]。速度場(chǎng)設(shè)定為[16]

u=-sin2(πx)sin(2πy)，v=sin2(πy)sin(2πx)

在t=1時(shí)，將整個(gè)速度場(chǎng)進(jìn)行翻轉(zhuǎn)。因此，理論上單渦會(huì)在t=2時(shí)回到其初始位置。計(jì)算時(shí)分別采用了兩套規(guī)則的三角形網(wǎng)格1(Δx=Δy=1/100)和網(wǎng)格2(Δx=Δy=1/100)，時(shí)間步長(zhǎng)均取為Δt=0.001。

圖2中黑色實(shí)線、紅色虛線以及藍(lán)色虛線分別表示t=2時(shí)的精確解以及Mesh 1和Mesh 2上的數(shù)值解。在Mesh 1和Mesh 2上得到數(shù)值結(jié)果的相對(duì)質(zhì)量誤差分別為1.1%和0.5%，驗(yàn)證了數(shù)值算法較好的質(zhì)量守恒性。可以看出，數(shù)值解具有很好的網(wǎng)格收斂性，而且和精確解吻合較好。

圖2 t =2時(shí)刻精確解和數(shù)值解的比較

6.2 帶有障礙物的牛頓流體潰壩流動(dòng)

圖3 下游帶有障礙物潰壩流動(dòng)問題示意圖

初始時(shí)刻，液柱右側(cè)會(huì)有一個(gè)擋板，當(dāng)瞬間抽走這個(gè)擋板時(shí)，液柱會(huì)在重力的作用下發(fā)生坍塌，向右下方流動(dòng)。整個(gè)流動(dòng)過程由于中間障礙物的存在，液體前沿界面會(huì)發(fā)生碰撞和飛濺等大變形現(xiàn)象。

計(jì)算過程中用到了三套逐次加密的非結(jié)構(gòu)網(wǎng)格，網(wǎng)格1，2和3的最大網(wǎng)格尺寸分別為1/100，1/120和1/150，求解區(qū)域內(nèi)的小三角形單元數(shù)分別為8983，13386和20539，圖4(a)給出了計(jì)算過程用到的網(wǎng)格2。本文所有的數(shù)值模擬中，時(shí)間步長(zhǎng)均取為0.0002。圖4(b)展示了三套計(jì)算網(wǎng)格上，t=0.2 s時(shí)的前沿界面形態(tài)。其中，紅色、黑色和藍(lán)色線條分別表示在網(wǎng)格1，2和3上計(jì)算得到的數(shù)值結(jié)果，可以看出，耦合算法具有很好的網(wǎng)格收斂性。為了保證數(shù)值結(jié)果的精確性并減少計(jì)算量，后續(xù)的數(shù)值模擬均在網(wǎng)格2上進(jìn)行。

圖4 計(jì)算用的網(wǎng)格2及三套計(jì)算網(wǎng)格上t =0.2 s時(shí)的界面形態(tài)

圖5給出了在網(wǎng)格2上計(jì)算得到的(t=0.1 s,0.3 s和0.5 s)液體的界面位置和形態(tài)。本文數(shù)值結(jié)果和實(shí)驗(yàn)結(jié)果均吻合較好，這驗(yàn)證了耦合算法的可行性以及準(zhǔn)確性。整個(gè)流動(dòng)過程中，流體質(zhì)量的最大誤差為1.2%，證明了耦合算法較好的質(zhì)量守恒性。本文數(shù)值模擬結(jié)果的最大相對(duì)質(zhì)量誤差，要比文獻(xiàn)[10]穩(wěn)定化有限元方法計(jì)算得到的最大相對(duì)質(zhì)量誤差小一些。

圖5 不同時(shí)刻牛頓流體的前沿界面形態(tài)

圖6(a)給出了左側(cè)邊界上流體高度隨時(shí)間的變化，圖6(b)給出了障礙物左下角(點(diǎn)A)處的壓力隨時(shí)間的變化。當(dāng)液柱坍塌下來沖擊到障礙物左下點(diǎn)時(shí)，該點(diǎn)的壓力瞬間增加，峰值約為6000 Pa。可以看出，本文的數(shù)值結(jié)果和文獻(xiàn)[7]幾乎一致，進(jìn)一步驗(yàn)證了數(shù)值算法的準(zhǔn)確性。表1給出了若干時(shí)刻障礙物左下角(點(diǎn)A)處液面的垂直高度。

圖6 左邊界上液面高度及障礙物左下角(點(diǎn)A)處的壓力隨時(shí)間的變化

表1 障礙物左下角(點(diǎn)A)處若干時(shí)刻液面的垂直高度

6.3 非牛頓流體潰壩流動(dòng)

圖7給出了若干時(shí)刻液體的前沿界面形態(tài)，其中淺灰色表示本文的計(jì)算結(jié)果，深灰色表示文獻(xiàn)的計(jì)算結(jié)果。

圖7 流體前沿界面形態(tài)

另外，圖8進(jìn)一步展示了液體前沿界面位置隨著時(shí)間的變化曲線。可以看出，本文的數(shù)值結(jié)果和文獻(xiàn)[8]已有的計(jì)算結(jié)果以及實(shí)驗(yàn)結(jié)果[18]均吻合得很好，這驗(yàn)證了數(shù)值算法處理非牛頓流體潰壩流動(dòng)的準(zhǔn)確性。

圖8 x軸上液體前沿界面位置

6.4 帶有障礙物的非牛頓流體潰壩流動(dòng)

繼續(xù)研究非牛頓流體的潰壩流動(dòng)情況，整個(gè)求解區(qū)域及邊界條件和6.2節(jié)一致，所有的數(shù)值模擬均在圖4給出的網(wǎng)格2上進(jìn)行。首先考慮n=0.5時(shí)的假塑性非牛頓流體，圖9給出了若干時(shí)刻液體的自由運(yùn)動(dòng)界面形態(tài)，可以看出，這時(shí)前沿界面相比于圖5的結(jié)果更為紊亂一些，尤其是在和右壁面碰撞之后，液體的飛濺和破碎現(xiàn)象更加明顯，這些都是因?yàn)殡S著流動(dòng)速度的增加，假塑性非牛頓流體的粘度會(huì)變小，更容易發(fā)生形變。

圖9 假塑性流體(n =0.5)潰壩流動(dòng)中若干時(shí)刻的運(yùn)動(dòng)界面形態(tài)

進(jìn)一步，圖10給出了障礙物左下角(點(diǎn)A)處壓力的變化曲線，可以看到其峰值約為6600 Pa，較高于6.2節(jié)的情況。另外，流動(dòng)過程中，流體質(zhì)量的最大相對(duì)質(zhì)量誤差為1.5%。

接著，本文進(jìn)一步考慮n=1.5時(shí)的脹塑性非牛頓流體，圖11展示了若干時(shí)刻流體的運(yùn)動(dòng)界面形態(tài)，相比于圖5和圖9的結(jié)果，這時(shí)的前沿界面更加光滑。這是因?yàn)槊浰苄苑桥ｎD流體在速度增加時(shí)，粘度會(huì)變大，不容易發(fā)生飛濺和破碎等現(xiàn)象，與文獻(xiàn)[8]的描述一致。圖12展示了障礙物左下角(點(diǎn)A)處壓力隨時(shí)間的變化曲線，壓力峰值約為4200 Pa，相比于牛頓流體和假塑性非牛頓流體，該值都要低很多。此外，流動(dòng)過程中，液體的最大相對(duì)質(zhì)量誤差為0.8%。

圖10 障礙物左下角(點(diǎn)A)壓力隨時(shí)間的變化

圖11 脹塑性流體(n =1.5)潰壩流動(dòng)中若干時(shí)刻的運(yùn)動(dòng)界面形態(tài)

圖12 障礙物左下角(點(diǎn)A)壓力隨時(shí)間的變化

7 結(jié) 論

本文基于兩相流動(dòng)模型，采用Level Set前沿界面追蹤方法，在有限元框架下數(shù)值模擬研究了下游帶有障礙物的牛頓和非牛頓潰壩流動(dòng)問題，并分析了不同特性流體自由界面的發(fā)展形態(tài)。

對(duì)于牛頓流體，本文數(shù)值模擬結(jié)果和文獻(xiàn)已有報(bào)道均吻合得很好，這表明本文耦合算法具有較高的準(zhǔn)確性。對(duì)于非牛頓流體，數(shù)值結(jié)果表明，假塑性流體對(duì)于障礙物的最大沖擊力最大，而且在流動(dòng)過程中液體界面最容易出現(xiàn)破碎和飛濺等大變形現(xiàn)象；脹塑性流體對(duì)于障礙物的最大沖擊力最小，在整個(gè)流動(dòng)過程中界面最為穩(wěn)定和光滑。整個(gè)流動(dòng)過程均保證了較好的質(zhì)量守恒性，所有模擬結(jié)果中最大的質(zhì)量誤差為1.5%。剪切變稀的非牛頓流體在潰壩流動(dòng)問題中，對(duì)于下游障礙物的破壞性最大。因此，當(dāng)實(shí)際問題涉及到這類流體物質(zhì)時(shí)，需要注意對(duì)下游障礙物的加固，防止流體將障礙物沖垮。