用于三維非均質(zhì)流場計算的改進流線迎風Petrov-Galerkin(SUPG)方法

蘇波,韓向科

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

用于三維非均質(zhì)流場計算的改進流線迎風Petrov-Galerkin(SUPG)方法

蘇波1,韓向科2

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

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

不可壓流動;有限元法;分裂算法;流線迎風

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

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

1 不可壓流體運動控制方程

一般流體運動的連續(xù)性方程為

(1)

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

(2)

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

(3)

(4)

(5)

式中:ρ為流體密度;Vi為流體速度分量;xi為直角坐標分量;t為時間。式(3)、式(5)可簡寫為

(6)

(7)

式中

2 控制方程的雙時間步離散

對式(6)、式(7)采用雙時間步法[7]進行求解

(8)

(9)

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

(10)

(11)

式中

(12)

(13)

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

3 求解流體運動控制方程的變量分裂法

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

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

第2步,不考慮壓強項及ΔVi,r項,計算中間速度增量

(14)

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

(15)

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

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

(16)

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

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

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

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

(17)

式中

(18)

對式(14)采用SUPG有限元格式進行空間離散,得到輔助動量方程有限元格式

(19)

式中

(20)

方程組(15)可簡寫為泊松型壓強方程[11],采用標準Galerkin有限元格式離散,得到計算壓強的有限元格式的弱形式

Δtphtp+ΔtpΔpe=rp

(21)

式中

(22)

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

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

(23)

其中

對速度修正方程,式(16)采用標準Galerkin有限元格式離散,有限元格式為

(24)

(25)

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

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

5.1 計算模型

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

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

(26)

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

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

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

表1 雙層液體流動模型計算條件

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

5.2 結(jié)果分析

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

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

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

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

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

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

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

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

圖6 y=0.05m中面流線圖

圖7 旋渦位置隨時間的變化

圖8 分界面位置隨時間的變化對比

6 結(jié) 論

本文所提出的改進SUPG方法摒棄了密度為常數(shù)的假定,采用包含密度輸運方程的黏性不可壓流體運動控制方程組,使之可用于非均質(zhì)流場的分析。由于采用變量分裂算法,速度場和壓強場可采用同階插值函數(shù)進行空間離散,使改進SUPG方法具有簡明的有限元格式,并且速度場與壓強場可依次求解,降低了線性代數(shù)方程組的階數(shù),有利于大規(guī)模問題的求解。改進SUPG方法引入雙時間步法進行時間離散,有利于提高非定常問題的求解精度。算例分析表明,改進SUPG方法可對重力作用下三維矩形管道內(nèi)的自由流動問題進行有效分析,得到穩(wěn)定的速度場、壓強場以及渦流變化規(guī)律,從而實現(xià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] 錢若軍, 董石麟, 袁行飛.流固耦合理論研究進展 [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].力學與實踐, 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)文獻鏈接]

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

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

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

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

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

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

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

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

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

蘇波,錢若軍,韓向科.一種用于流固耦合分析的有限元網(wǎng)格簡捷更新方法.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-),男,博士,講師。

國家自然科學基金青年基金資助項目

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

時間: 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