蘇波,韓向科

(1.江蘇大學(xué)土木工程與力學(xué)學(xué)院, 212013, 江蘇鎮(zhèn)江; 2.中冶天工上海十三冶建設(shè)有限公司, 200092, 上海)

用于三維非均質(zhì)流場(chǎng)計(jì)算的改進(jìn)流線迎風(fēng)Petrov-Galerkin(SUPG)方法

蘇波1,韓向科2

(1.江蘇大學(xué)土木工程與力學(xué)學(xué)院, 212013, 江蘇鎮(zhèn)江; 2.中冶天工上海十三冶建設(shè)有限公司, 200092, 上海)

在流線迎風(fēng)Petrov-Galerkin(SUPG)穩(wěn)定有限元方法基礎(chǔ)上,通過(guò)引入雙時(shí)間步法和變量分裂算法,發(fā)展了一種可用于三維非均質(zhì)流場(chǎng)計(jì)算的改進(jìn)SUPG方法。該方法摒棄了傳統(tǒng)不可壓縮流動(dòng)問(wèn)題中密度為常數(shù)的假定,采用包含密度輸運(yùn)方程的流體運(yùn)動(dòng)控制方程;基于變量分裂算法,速度、壓強(qiáng)場(chǎng)采用同階插值函數(shù)進(jìn)行空間離散,使改進(jìn)SUPG方法具有簡(jiǎn)明的有限元格式,同時(shí)對(duì)速度場(chǎng)、壓強(qiáng)場(chǎng)進(jìn)行迭代求解,以降低線性代數(shù)方程組的階數(shù)。雙時(shí)間步法的引入,有利于提高SUPG方法對(duì)復(fù)雜非定常問(wèn)題的求解效率。采用該方法對(duì)非均質(zhì)、非定常三維矩形管道重力作用下的自由流動(dòng)問(wèn)題進(jìn)行了分析,研究了重力作用下兩種不同密度液體的相對(duì)運(yùn)動(dòng)過(guò)程。計(jì)算分析表明:在采用較大時(shí)間步的情況下,速度場(chǎng)和壓強(qiáng)場(chǎng)在整個(gè)流動(dòng)過(guò)程中隨時(shí)間平穩(wěn)過(guò)渡且分布光滑,沒(méi)有出現(xiàn)數(shù)值波動(dòng)現(xiàn)象;旋渦位置及其隨時(shí)間變化的規(guī)律與經(jīng)典文獻(xiàn)結(jié)果相符,沒(méi)有出現(xiàn)跳躍和不連續(xù)現(xiàn)象。算例分析表明,改進(jìn)SUPG方法具有良好的計(jì)算精度及數(shù)值穩(wěn)定性,可用于三維非均質(zhì)流動(dòng)類似問(wèn)題的研究。

不可壓流動(dòng);有限元法;分裂算法;流線迎風(fēng)

有限元方法由于可采用非結(jié)構(gòu)網(wǎng)格、易處理復(fù)雜邊界條件等優(yōu)點(diǎn),近些年來(lái)在工程問(wèn)題中的流體計(jì)算領(lǐng)域得到了廣泛研究和重視[1-5]。但是,當(dāng)對(duì)流項(xiàng)占優(yōu)時(shí),采用經(jīng)典Garlerkin法求解流體N-S方程會(huì)引起數(shù)值波動(dòng)。Brooks和Hughes針對(duì)定常對(duì)流擴(kuò)散問(wèn)題,提出了流線迎風(fēng)Petrov-Galerkin(SUPG)方法,有效地解決了對(duì)流項(xiàng)引起的數(shù)值波動(dòng)問(wèn)題[6]。Jameson提出了雙時(shí)間步法,在計(jì)算精度和計(jì)算效率上具有較高的優(yōu)越性,適合定常、非定常問(wèn)題的求解[7]。此外,對(duì)于黏性不可壓流動(dòng)問(wèn)題,通過(guò)變量分裂法[8]將速度場(chǎng)和壓強(qiáng)場(chǎng)解耦,直接避免了LBB(Ladyzhenskaya-Babǔska-Brezzi)條件的限制,可以使單元速度、壓強(qiáng)采用同階插值函數(shù),易于對(duì)混合場(chǎng)進(jìn)行有限元空間離散。

在對(duì)黏性不可壓流動(dòng)問(wèn)題進(jìn)行有限元求解時(shí),通常假定密度為常數(shù)[2,9],即采用均質(zhì)流場(chǎng)進(jìn)行求解。這樣雖然簡(jiǎn)化了黏性不可壓流動(dòng)問(wèn)題的計(jì)算過(guò)程,但同時(shí)給問(wèn)題的分析帶來(lái)局限性。為能對(duì)非均質(zhì)流場(chǎng)進(jìn)行精確分析,本文摒棄了密度為常數(shù)的假定,采用原始黏性不可壓流體運(yùn)動(dòng)控制方程,繼而綜合運(yùn)用上述各種數(shù)值方法,提出一種適用于非定常黏性不可壓流動(dòng)問(wèn)題求解的改進(jìn)SUPG有限元方法,并給出相應(yīng)的有限元計(jì)算步驟和格式。最后,對(duì)三維管道中變密度液體在重力作用下的自由流動(dòng)問(wèn)題進(jìn)行分析,以驗(yàn)證該方法的有效性。

1 不可壓流體運(yùn)動(dòng)控制方程

一般流體運(yùn)動(dòng)的連續(xù)性方程為

(1)

對(duì)于黏性不可壓流體,根據(jù)定義,流體的密度ρ在運(yùn)動(dòng)過(guò)程中保持不變,即

(2)

當(dāng)假定密度ρ為常數(shù)時(shí),式(2)自動(dòng)滿足,質(zhì)量守恒僅需滿足速度散度為0的條件。但是,對(duì)于非均勻流場(chǎng),式(2)中隱含了密度場(chǎng)運(yùn)動(dòng)條件,是不可忽略的,因此,對(duì)于非均質(zhì)黏性不可壓流體,其運(yùn)動(dòng)控制方程為

(3)

(4)

(5)

式中:ρ為流體密度;Vi為流體速度分量;xi為直角坐標(biāo)分量;t為時(shí)間。式(3)、式(5)可簡(jiǎn)寫(xiě)為

(6)

(7)

式中

2 控制方程的雙時(shí)間步離散

對(duì)式(6)、式(7)采用雙時(shí)間步法[7]進(jìn)行求解

(8)

(9)

對(duì)式(8)、式(9)中關(guān)于虛擬時(shí)間tp的導(dǎo)數(shù)項(xiàng)采用加權(quán)隱格式[10]進(jìn)行時(shí)間離散,關(guān)于物理時(shí)間的導(dǎo)數(shù)項(xiàng)采用二階精度三點(diǎn)向后差分格式,得

(10)

(11)

式中

(12)

(13)

上角標(biāo)t表示物理時(shí)刻t的已知變量;tp表示虛擬時(shí)刻tp的已知變量;t+Δt表示物理時(shí)刻t+Δt的待求變量;tp+Δtp表示虛擬時(shí)刻tp+Δtp的待求變量;Δt為物理時(shí)間步長(zhǎng);Δtp為虛擬時(shí)間步長(zhǎng);tp+ΔtpΔρ為虛擬時(shí)間的密度增量,tp+ΔtpΔρ=tp+Δtpρ-tpρ;tp+ΔtpΔV為虛擬時(shí)間的速度增量,tp+ΔtpΔV=tp+ΔtpV-tpV;θρ、θV、θp為時(shí)間離散權(quán)系數(shù)。

3 求解流體運(yùn)動(dòng)控制方程的變量分裂法

采用變量分裂算法[8],對(duì)式(10)、式(11)方程組的速度場(chǎng)、壓強(qiáng)場(chǎng)進(jìn)行迭代求解,計(jì)算步驟如下:

第1步,由式(10)求出時(shí)刻tp+Δtp的密度增量tp+ΔtpΔρ及密度全量tp+Δtpρ=tpρ+tp+ΔtpΔρ;

第2步,不考慮壓強(qiáng)項(xiàng)及ΔVi,r項(xiàng),計(jì)算中間速度增量

(14)

第3步,將中間速度增量tp+ΔtpΔVint及密度全量tp+Δtpρ代入下式,求解壓強(qiáng)tp+Δtpp

(15)

第4步,根據(jù)tp+ΔtpΔVint和tp+θpΔtpp,求解速度增量

tp+ΔtpΔV=tp+ΔtpΔVint-

(16)

式(10)、式(14)～(16)即構(gòu)成了基于變量分裂方法的非定常不可壓流體運(yùn)動(dòng)的基本方程。當(dāng)流場(chǎng)為均質(zhì)流體,即密度不隨空間改變時(shí),連續(xù)性方程(3)自動(dòng)滿足,則上述方法退化為一般不可壓流動(dòng)分裂算法。

4 基于SUPG方法的不可壓流動(dòng)穩(wěn)定有限元格式

對(duì)式(10)采用SUPG方法[1,6]進(jìn)行空間離散,得到連續(xù)性方程的有限元格式

mρtp+ΔtpΔρe=rρ

(17)

式中

(18)

對(duì)式(14)采用SUPG有限元格式進(jìn)行空間離散,得到輔助動(dòng)量方程有限元格式

(19)

式中

(20)

方程組(15)可簡(jiǎn)寫(xiě)為泊松型壓強(qiáng)方程[11],采用標(biāo)準(zhǔn)Galerkin有限元格式離散,得到計(jì)算壓強(qiáng)的有限元格式的弱形式

Δtphtp+ΔtpΔpe=rp

(21)

式中

(22)

rp=θp(tp+Δtpgp-Δtphtp+Δtppe-tp+Δtpfρ)+

(1-θp)(tpgp-Δtphtppe-tpfρ)

(23)

其中

對(duì)速度修正方程,式(16)采用標(biāo)準(zhǔn)Galerkin有限元格式離散,有限元格式為

(24)

(25)

式(17)、(19)、(21)、(24)即構(gòu)成了可用于非均質(zhì)流場(chǎng)計(jì)算的改進(jìn)SUPG穩(wěn)定有限元格式。

5 非均質(zhì)、非定常三維矩形管道重力作用下的自由流動(dòng)問(wèn)題

5.1 計(jì)算模型

在非均質(zhì)流動(dòng)算例中,雙層流體在重力作用下的自由流動(dòng)問(wèn)題[12]是一個(gè)十分典型的算例。如圖1所示,有一矩形流場(chǎng),流場(chǎng)內(nèi)盛有2種不同密度的流體,密度分別為ρ1和ρ2,且ρ2>ρ1,密度較大的流體位于流場(chǎng)上部,密度較小的流體位于流場(chǎng)下部。參考二維算例[13],初始密度分布為

ρ(x,y,z,t=0)=

(26)

式中:η(x)=-0.1dcos(2πx/d);d為指定長(zhǎng)度。

(a)流場(chǎng)幾何尺寸 (b)有限元網(wǎng)格劃分

有限元網(wǎng)格劃分形式見(jiàn)圖1b,計(jì)算雷諾數(shù)定義為Re=ρ1d3/2g1/2/μ[14],其中g(shù)為重力加速度,μ為液體的動(dòng)力黏性系數(shù)。其他計(jì)算條件見(jiàn)表1。

表1 雙層液體流動(dòng)模型計(jì)算條件

在重力作用下,流場(chǎng)中2種液體發(fā)生相對(duì)運(yùn)動(dòng)。雖然流場(chǎng)形狀簡(jiǎn)單,但隨著上部液體向下流動(dòng),2種液體逐漸融合交織,導(dǎo)致流動(dòng)狀態(tài)十分復(fù)雜,流動(dòng)對(duì)網(wǎng)格形式、邊界條件等因素十分敏感,容易引起數(shù)值求解失敗[13]。

5.2 結(jié)果分析

圖2給出了流場(chǎng)在不同時(shí)刻的密度分布規(guī)律。在重力作用下,矩形流場(chǎng)上部密度較大的液體向下流動(dòng)。當(dāng)t=1 s時(shí),2種液體分界面逐漸擴(kuò)大,但流場(chǎng)分布規(guī)律仍和初始分布規(guī)律基本相似;當(dāng)t=1.5s時(shí),上部液體繼續(xù)向下流動(dòng),且在兩液體交界面處開(kāi)始形成旋渦,但2種液體仍然具有較為明顯的分界面;當(dāng)t=1.75s時(shí),與t=1.5s時(shí)的分布規(guī)律基本一致,2種液體交界面處旋渦更加明顯;當(dāng)t=2 s時(shí),2種液體的交界面逐漸擴(kuò)大,并相互交織,旋渦繼續(xù)增強(qiáng),旋渦中心逐漸上移;當(dāng)t>2 s時(shí),流場(chǎng)的流動(dòng)特性開(kāi)始發(fā)生改變,隨著2種液體的交界面繼續(xù)擴(kuò)大,交界面逐漸模糊,流場(chǎng)中部2種液體已相互滲透,流動(dòng)情況愈加復(fù)雜,這與文獻(xiàn)[13-14]中給出的分布規(guī)律相同。

圖3～圖5給出了整個(gè)流動(dòng)過(guò)程中的壓強(qiáng)和速度場(chǎng)分布規(guī)律。如圖所示,在整個(gè)時(shí)段內(nèi),盡管流動(dòng)狀態(tài)復(fù)雜,但壓強(qiáng)和速度分布光滑,隨時(shí)間變化過(guò)渡平穩(wěn),可以得到穩(wěn)定的壓強(qiáng)場(chǎng)和速度場(chǎng),沒(méi)有出現(xiàn)數(shù)值波動(dòng)現(xiàn)象。

圖6和圖7分別給出了流動(dòng)過(guò)程中的流線圖和流場(chǎng)中旋渦位置隨時(shí)間的變化規(guī)律。在整個(gè)流動(dòng)過(guò)程中,旋渦位置在x方向變化幅度較小,在約t=1 s之前旋渦略向左側(cè)靠攏,隨后逐漸向右側(cè)邊發(fā)展;相比之下,旋渦位置在z方向變化幅度較為明顯,在t=0.5s之前z向坐標(biāo)值變化較為緩慢,隨后迅速變小,至約1.8 s左右又逐漸增大,整體呈正弦曲線規(guī)律。圖8給出了流場(chǎng)的上部和下部液體分界面位置隨時(shí)間的變化規(guī)律,可以看出,本文的計(jì)算結(jié)果與文獻(xiàn)[14]給出的計(jì)算結(jié)果十分吻合。

為使求解收斂,文獻(xiàn)[14]采用的歸一化時(shí)間步為Δtnon=5×10-4,文獻(xiàn)[15]選擇的歸一化時(shí)間步為Δtnon=0.001 25(At)1/2=8.84×10-4,而本算例采用較大的歸一化時(shí)間步:Δtnon=t(At)1/2=0.1(At)1/2=7.07×10-2。計(jì)算表明,在整個(gè)流動(dòng)過(guò)程中旋渦位置變化平穩(wěn),沒(méi)有出現(xiàn)任何跳躍和不連續(xù)現(xiàn)象,依然可以獲得良好的計(jì)算精度,可見(jiàn)本文提出的改進(jìn)SUPG方法在時(shí)間方向具有良好的穩(wěn)定性能。

圖2 密度ρ的分布規(guī)律

圖3 速度Vx的分布規(guī)律

圖4 速度Vz分布規(guī)律

圖5 壓強(qiáng)p的分布規(guī)律

圖6 y=0.05m中面流線圖

圖7 旋渦位置隨時(shí)間的變化

圖8 分界面位置隨時(shí)間的變化對(duì)比

6 結(jié) 論

本文所提出的改進(jìn)SUPG方法摒棄了密度為常數(shù)的假定,采用包含密度輸運(yùn)方程的黏性不可壓流體運(yùn)動(dòng)控制方程組,使之可用于非均質(zhì)流場(chǎng)的分析。由于采用變量分裂算法,速度場(chǎng)和壓強(qiáng)場(chǎng)可采用同階插值函數(shù)進(jìn)行空間離散,使改進(jìn)SUPG方法具有簡(jiǎn)明的有限元格式,并且速度場(chǎng)與壓強(qiáng)場(chǎng)可依次求解,降低了線性代數(shù)方程組的階數(shù),有利于大規(guī)模問(wèn)題的求解。改進(jìn)SUPG方法引入雙時(shí)間步法進(jìn)行時(shí)間離散,有利于提高非定常問(wèn)題的求解精度。算例分析表明,改進(jìn)SUPG方法可對(duì)重力作用下三維矩形管道內(nèi)的自由流動(dòng)問(wèn)題進(jìn)行有效分析,得到穩(wěn)定的速度場(chǎng)、壓強(qiáng)場(chǎng)以及渦流變化規(guī)律,從而實(shí)現(xiàn)對(duì)非均勻、非定常等復(fù)雜流動(dòng)問(wèn)題的求解。

[1] HUANG Cheng, ZHOU Dai, BAO Yan, et al.A stabilized finite element technique and its application for turbulent flow with high Reynolds number [J].Wind and Structures, 2011, 14(5): 465-480.

[2] ZIENKIEWICZ O C, TAYLOR R L.The finite element method for fluid dynamics [M].5th ed.Amsterdam, The Netherlands: Elsevier, 2000.

[3] 錢(qián)若軍, 董石麟, 袁行飛.流固耦合理論研究進(jìn)展 [J].空間結(jié)構(gòu), 2008, 4(1): 3-15.QIAN Ruo-jun, DONG Shi-lin, YUAN Xing-fei.Advances in research on fluid-structure interaction theory [J].Spacial Structures, 2008, 4(1): 3-15.

[4] 韓向科, 蘇波.黏性不可壓流體的一種FCBIS三角形單元 [J].力學(xué)與實(shí)踐, 2010, 32(3): 22-25.HAN Xiangke, SU Bo.A FCBIS triangular element for incompressible fluid flows [J].Mechanics in Engineering, 2010, 32(3): 22-25.

[5] KOHNO H, BATHE K J.A nine-node quadrilateral FCBI element for incompressible fluid flows [J].Communications in Numerical Methods in Engineering, 2006, 22(8): 917-931.

[6] BROOKS A N, HUGHES T J R.Streamline upwind/Petrov-Galerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations [J].Computer Methods in Applied Mechanics and Engineering, 1982, 32(1/2/3): 199-259.

[7] JAMESON A.Time dependent calculation using multigrid with applications to unsteady flows past airfoils and wings [C]∥Proceedings of AIAA 10th Computational Fluid Dynamics Conference.Reston, VA, USA: AIAA, 1991: 1-8.

[8] NITHIARASU P, CODINA R, ZIENKIEWICZ O C.The characteristic-based split (CBS) scheme: a unified approach to fluid dynamics [J].Numerical Methods in Engineering, 2006, 66(10): 1514-1546.

[9] ADINA R&D Inc.ADINA CFD&FSI: theory and modeling guide [M].Watertown, MA, USA: ADINA R&D Inc., 2005.

[10]DONEA J, HUERTA A.Finite element methods for flow problems [M].London, UK: John Wiley & Sons, 2003: 92-93.

[11]GRESHO P M, SANI R L.Incompressible flow and the finite element method [M].New York, USA: John Wiley & Sons, 2000: 302.

[12]TRYGGVASON G.Numerical simulation of the Rayleigh-Taylor instability [J].Journal of Computational Physics, 1988, 75(2): 253-282.

[13]CALGARO C, CREUSE E, GOUDON T.A hybrid finite volume-finite element method for variable density incompressible flows [J].Journal of Computational Physics, 2008, 227(9): 4671-4696.

[14]GUERMOND J L, QUARTAPELLE L.A projection FEM for variable density incompressible flows [J].Journal of Computational Physics, 2000, 165(1): 167-188.

[15]GUERMOND J L, SALGADO A.A splitting method for incompressible flows with variable density based on a pressure Poisson equation [J].Journal of Computational Physics, 2009, 228(8): 2834-2846.

[本刊相關(guān)文獻(xiàn)鏈接]

趙立波,徐龍起,熱合曼艾比布力,等.矩形微懸臂梁的流固耦合諧振頻率分析.2013,47(11):60-64.[doi:10.7652/xjtuxb201311011]

何聯(lián)格,左正興,向建華.氣缸蓋中兩相流沸騰換熱熱機(jī)耦合仿真分析.2013,47(7):23-28.[doi:10.7652/xjtuxb201307 005]

李連忠,牛文全,魏正英,等.微灌壓力調(diào)節(jié)器調(diào)壓特性流固耦合數(shù)值計(jì)算及分析.2013,47(5):131-136.[doi:10.7652/xjtuxb201305024]

何聯(lián)格,左正興,向建華.氣缸蓋冷卻水腔內(nèi)兩相流動(dòng)沸騰傳熱仿真研究.2013,47(1):21-26.[doi:10.7652/xjtuxb 201301005]

張穎莉,種道彤,劉繼平,等.方管內(nèi)混合對(duì)流與管壁導(dǎo)熱耦合換熱的數(shù)值模擬.2012,46(5):19-24.[doi:10.7652/xjtuxb201205004]

劉小民,王星,許運(yùn)賓.運(yùn)動(dòng)罐體內(nèi)液體晃動(dòng)的雙向流固耦合數(shù)值分析.2012,46(5):120-126.[doi:10.7652/xjtuxb201205021]

張慧賢,寇子明,吳娟,等.液壓激波作用下管道流固耦合的動(dòng)力學(xué)建模.2012,46(3):94-99.[doi:10.7652/xjtuxb201203 017]

康偉,張家忠,李凱倫.利用本征正交分解的非線性Galerkin降維方法.2011,45(11):58-62.[doi:10.7652/xjtuxb201111 011]

梅冠華,張家忠.時(shí)滯慣性流形在三維壁板顫振數(shù)值分析中的應(yīng)用.2011,45(9):40-46.[doi:10.7652/xjtuxb201109008]

蘇波,錢(qián)若軍,韓向科.一種用于流固耦合分析的有限元網(wǎng)格簡(jiǎn)捷更新方法.2011,45(3):16-24.[doi:10.7652/xjtuxb 201103003]

(編輯 葛趙青)

ModifiedStreamlineUpwindPetrov-Galerkin(SUPG)StrategyforCalculationofThree-DimensionalHeterogeneousFlowProblems

SU Bo1,HAN Xiangke2

(1.Faculty of Civil Engineering and Mechanics, Jiangsu University, Zhenjiang, Jiangsu 212013, China;2.CMTCC Shanghai Shisanye Construction Co., Ltd., Shanghai 200092, China)

By introducing the dual-time stepping method and the variable splitting algorithm, a modified SUPG strategy is developed based on the stable streamline upwind Petrov-Galerkin (SUPG) finite element method, where the traditional assumption of constant density for incompressible flow problem is abandoned and the density transportation equation is introduced into the governing equations.The velocity and pressure fields are discretized with interpolating function of the same order, thus finite element scheme of the modified SUPG strategy gets simple and clear, and the order of algebraic equations is reduced.The dual-time stepping method is also introduced to enhance calculation stability of the SUPG strategy for complex unsteady problems.A free flow with heterogeneous, unsteady field of three-dimensional rectangular pipe under gravity and the whole relative motion between two types of liquids with different density are analyzed.The calculation results indicate that the velocity and pressure fields distribute and transmit smoothly with time and no numerical wave occurs in the case of greater time steps; the vortex position and its varying regulation coincide well with the results in the classic literatures without jumps and discontinuities.Example proves the numerical stability and accuracy of the modified SUPG method.

incompressible flow; finite element method; splitting algorithm; streamline upwind Petrov-Galerkin (SUPG) method

10.7652/xjtuxb201403023

2013-08-05。

蘇波(1977-),男,博士,講師。

國(guó)家自然科學(xué)基金青年基金資助項(xiàng)目

(51108210);江蘇省博士后基金資助項(xiàng)目(1301048C);國(guó)家自然科學(xué)基金專項(xiàng)數(shù)學(xué)天元基金資助項(xiàng)目(11226308)。

時(shí)間: 2013-12-25

O357.1

:A

:0253-987X(2014)03-0128-07

網(wǎng)絡(luò)出版地址: http:∥www.cnki.net/kcms/detail/61.1069.T.20131225.1702.005.html