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

蘇波,韓向科

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

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

蘇波1,韓向科2

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

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

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

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

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

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

一般流體運動的連續性方程為

(1)

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

(2)

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

(3)

(4)

(5)

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

(6)

(7)

式中

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

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

(8)

(9)

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

(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為時間離散權系數。

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步,根據tp+ΔtpΔVint和tp+θpΔtpp,求解速度增量

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

(16)

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

4 基于SUPG方法的不可壓流動穩定有限元格式

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

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)即構成了可用于非均質流場計算的改進SUPG穩定有限元格式。

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

5.1 計算模型

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

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

(26)

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

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

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

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

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

5.2 結果分析

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

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

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

為使求解收斂,文獻[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。計算表明,在整個流動過程中旋渦位置變化平穩,沒有出現任何跳躍和不連續現象,依然可以獲得良好的計算精度,可見本文提出的改進SUPG方法在時間方向具有良好的穩定性能。

圖2 密度ρ的分布規律

圖3 速度Vx的分布規律

圖4 速度Vz分布規律

圖5 壓強p的分布規律

圖6 y=0.05m中面流線圖

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

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

6 結 論

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

[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].空間結構, 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.

[本刊相關文獻鏈接]

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

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

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

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

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

劉小民,王星,許運賓.運動罐體內液體晃動的雙向流固耦合數值分析.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]

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

蘇波,錢若軍,韓向科.一種用于流固耦合分析的有限元網格簡捷更新方法.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);國家自然科學基金專項數學天元基金資助項目(11226308)。

時間: 2013-12-25

O357.1

:A

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

網絡出版地址: http:∥www.cnki.net/kcms/detail/61.1069.T.20131225.1702.005.html