Navier-Stokes方程在流体力学中有很重要的意义,是一个非线性偏微分方程,对于Navier-Stokes方程问题的数值求解,国内外学者都进行了大量研究.2004年,曹伟[1]应用高阶Taylor展开式给出了时间增量的Galerkin方法,并对不可压缩流体建立了流动分析模型.2005年,Nagrath等[2]通过SUPG(streamline upwind Petrov Galerkin)方法数值模拟了不可压缩Navier-Stokes方程.2007年,马飞遥等[3]基于粗网格及细网格上的函数空间,将定常的不可压缩的Navier-Stokes方程转化为一个线性问题来进行求解.2008年,骆艳、冯民富[4]对可压缩线性Navier-Stokes方程,利用压力梯度投影与间断有限元法相结合,给出了一种求解此方程的稳定格式.2011年,Cho等[5]采用Taylor-Galerkin方法离散LS函数的平流方程,采用P1P1分裂有限元方法求解Navier-Stokes方程及基于结构网格模拟了不可压缩问题.2013年,Heimann等[6]在结构网格上基于非对称内部惩罚方法,给出了求解不相容不可压缩两相流问题的间断Galerkin方法.2013年,章争荣[7]以Galerkin加权余量积分式为基础,对不可压缩黏性流动Navier-Stokes方程速度与压力耦合直接数值求解.2015年,郭虹平等[8]采用两相流问题的统一间断有限元框架对气泡上升过程进行了模拟.2016年,秦望龙等[9]基于六面体网格单元,采用插值方法将物面的四边形面网格单元构造为弯曲面网格单元,用于三维可压缩Euler方程和Navier-Stokes方程的求解,并通过对三维圆球绕流以及旋转流线体绕流进行数值求解,验证了DG方法的高精度特性.2018年,Amrouche等[10]证明了不可压缩Navier-Stokes方程关于Navier边界条件下的Stokes算子与滑移条件有关且在能量空间中收敛到无滑移边界条件的Navier-Stokes方程的解,在线性系统和非线性系统的弱解和强解相对于滑移系数都是一致有界的.2019年,Kirk等[11]提出了压阻混合非连续Galerkin法(pressure-robust hybridized discontinuous Galerkin method)证明了定常Navier-Stokes方程的速度误差估计与压力无关,此外,还表明速度和压力的估计都是最优的, 数值算例表明此方法是全局最优收敛的.
本文基于分裂算法,采用Adams-Bashforth二阶时间离散格式[12]和IPDG空间离散格式来求解不可压Navier-Stokes方程.需要指出的是,鉴于非结构网格模拟复杂区域上实际问题的最有效技术手段是建立DG求解框架,目前基于非结构网格实现Navier-Stokes方程问题的统一DG求解寥寥无几.所以本文框架将LDG通量直接推广到区域,采用IPDG方法求解不可压Navier-Stokes方程,两端强加Dirichlet边界条件似乎破坏了最优收敛阶,但在实际计算中,避免精确计算这些内积,优化了算法实施和编程实现.结果证实达到了最优阶收敛,且收敛到O(hN+1/2).
二维不可压缩流体动力学Navier-Stokes方程描述的是瞬态、不可压、等温流动问题,其形式如下:
)u=-p+ν2u, x∈Ω,
u=0,
其中u=(u,v)为速度的分量,p为标量压力场.写成守恒通量形式为
·F=-p+ν2u, u=0,
(1)
其中通量F为
考虑刚性稳定的时间步长法[12],将第n层到第n+1层每个时间步的推进分为3个阶段.第一步,用Adams-Bashforth二阶格式对守恒律部分作显式积分计算预估速度;第二步,将更新的速度分量投影到弱散度为零的空间计算压力;第三步,黏性项用隐式格式处理更新中间速度,最后求解un+1.
定义非线性函数N(u)=(Nx,Ny)=·F,则
Nx(u)=
(2)
第一步,用二阶Adams-Bashforth方法将其表示为
(3)
当系数γ0=0,α0=1,α1=0,β0=1,β1=0时,其格式降为向前Euler法,这与初始条件u0和v0保持一致;对以后时间步长取γ0=3/2,α0=2,α1=-1/2,β0=2,β1=-1.
第二步,压力投影步.为了系统完备,我们寻求使得中间速度场的散度为零,由来更新,即满足
(4)
在流出边界取Dirichlet压力边界条件;在流入边界取Neumann边界条件,即有
其中,Du/Dt=∂u/∂t+·F(u)为质量导数,wn×un是时间t=tn时的涡度.再用更新中间速度场
第三步,对速度分量方程求解,即
2un+1.
(5)
将对流方程(3)、压力投影方程(4)及黏性方程(5)相加,得
2un+1.
利用迎风DG方法来提高非线性对流方程(3)对流占优的稳定性,用内部惩罚DG方法处理压力投影方程(4)及黏性方程(5),可以避免单元间的混杂解跳跃.
第一步,对非线性对流方程(3)选Lax-Friedrichs通量,利用迎风DG方法[12]得
(φh,N)Ωk=(φh,INF(uh))∂Ωk+
uh)∂Ωk, ∀φh∈Vh,
(6)
其中被三角剖分成k个单元Ωk,PN(Ωk)为定义在Ωk上的N阶多项式空间,跳跃为u流入边界条件由u+给定,INf是f的N阶插值.
第二步,对压力投影方程(4)利用内部惩罚方法[12]处理,记ΓD和ΓN分别为Dirichlet和Neumann边界条件,且Γ为Ωk中单元边界的集合,∂Ω=ΓD∪ΓN,Γ0=Γ\∂Ω,取内部惩罚通量为uh}}-τu,对压力ph∈Vh,∀φh∈Vh得
φh·φh·phph}}·φhdΓ+
ph
φh·p0dΓ+
(7)
其中,惩罚因子[13]
第三步,对式(5)利用内部惩罚方法[12],对∀φh∈Vh,速度的x分量uh∈Vh满足如下方程:
φh·φh·uh
uh·τkdΩ-
(8)
速度的y分量vh∈Vh满足如下方程:
φh·φh·vh
vh·τkdΩ-
(9)
当非线性项与压力梯度相平衡,黏性项与动量方程的时间导数相平衡时,方程精确解为
u=-sin(2πy)e-ν·4π2t, v=sin(2πx)e-ν·4π2t, p=-cos(2πx)cos(2πy)e-ν·8π2t.
当ν很小时,时间依赖流具有解析解,可以定量检验数值算法的稳定性和精度.图1是时间依赖流的计算区域、边界条件及网格剖分[14].图2给出了ν=1/100, 1/1 000, 1/3 200, 1/5 000时时间依赖流的流线分布.观察图2可知ν取不同值时,方腔涡流左、右下角均出现一级涡,且随着ν的减小,涡的尺寸在增加,原始涡流慢慢基本趋于方腔的几何中心.当ν=1/3 200时,方腔右下角既出现一级涡还出现了一个较小的次级涡;当ν=1/5 000时,方腔右下角除了一级涡外还出现了两个很小的次级涡;由此说明模拟结果可准确地捕捉方腔右下角出现次级涡的现象.而表1表明当N=8时,由于时间误差的主导,使收敛阶为次优,当N为其他值时,几乎达到了最优收敛阶.该算例验证了本文方法具有很好的稳定性和较高的计算精度.
(a) 计算区域和边界条件(b) 网格剖分
(a) The calculation region and boundary conditions(b) Grid generation
图1 求解涡流问题示意图
Fig. 1 Schematic diagram for solving the eddy current problem
图2 不同黏性系数ν下时间依赖流的流线分布
Fig. 2 Streamline distributions of time-dependent flow under different viscosity coefficients ν
表1 涡流问题的不同逼近阶对应u与p最大节点误差及收敛阶
Table 1 Maximum u and p nodal errors of the vortex solution corresponding to different approximation orders
Nu errorhh/2h/4convergence orderp errorhh/2h/4convergence order18.4E-12.9E-15.4E-22.02.1E-05.4E-19.4E-22.227.4E-23.6E-21.0E-21.43.1E-18.8E-21.8E-82.031.2E-14.7E-35.6E-43.92.3E-11.6E-28.8E-44.041.0E-22.1E-34.5E-53.92.8E-24.1E-31.2E-43.951.8E-27.5E-53.0E-66.33.2E-23.3E-46.8E-66.166.0E-44.1E-51.7E-75.99.7E-46.7E-54.9E-75.571.4E-36.3E-74.5E-87.51.9E-33.9E-64.1E-87.781.8E-53.3E-72.7E-84.72.7E-55.4E-72.8E-85.0
考虑稍微复杂的流体问题,给出依赖时间的二维渠道流,在渠道中间有一圆形障碍物,此圆柱放置在偏离渠道中心线处,因此,一旦由零初始条件慢慢形成抛物型流时,就开始产生涡流.初始时刻T=0,流体假设为零,流入边界条件为
流出边界假设压力为零.因为此问题不存在解析解,定义圆柱的阻力系数为圆柱体上总力在水平方向分量,即
圆柱的升力系数为圆柱体上总力在垂直方向分量,即
其中为沿圆柱表面的外法向量,取压力差为Δp=p(-0.05,0)-p(0.05,0).图3是圆柱绕流的计算区域、边界条件及网格剖分.图4给出了单元K=236,每个单元用N=10阶逼近非稳态圆柱流的模拟结果,将文献[15-16]的结果作为标准的基础进行比较,表2显示了识别最大升力与阻力及何时达到这些极大值.
(a) 计算区域和边界条件
(a) The calculation region and boundary conditions
(b) 网格剖分
(b) Grid generation
图3 求解管道问题示意图
Fig. 3 Schematic diagram for solving pipeline problems
(a) x-速度 (b) y-速度
(a) The x-velocity (b) The y-velocity
(c) 压力 (d) 涡度
(c) The pressure (d) The vorticity
图4 非稳态圆柱流对应单元个数K=236, N=10阶的模拟结果
Fig. 4 Simulation results for K=236, N=10 in the unsteady pipe flow case
表2 圆柱管道流问题的最大升力、阻力系数及在t=8时的压力差
Table 2 Maximum lift, drag coefficients and final pressure drops in the pipe flow case for t=8
KNCdClΔp(t=8)11562.9530.497-0.11023682.9420.488-0.112236102.9550.479-0.112ref. [16]-2.9500.478-0.111
本文建立了一种求解二维Navier-Stokes方程黏弹流体流动的DG方法,采用非结构网格进行数值模拟,对于圆形障碍物,首先计算速度的梯度,再抽出在圆柱边上的速度,然后估计积分,最后计算出围绕圆柱的积分来估计升力及阻力.然而对大型时间问题的模拟存在很多挑战,二阶微分算子的存在可引进很小时间步长,对不可压缩Navier-Stokes情形,采用半隐式时间积分可避免小步长问题.
[1] 曹伟. 黏性不可压缩流体流动前沿的数值模拟[J]. 力学学报, 2004, 36(5): 583-588.(CAO Wei. Numerical simulation for the flow front of viscous incompressible fluid[J]. Chinese Journal of Theoretical and Applied Mechanics, 2004, 36(5): 583-588.(in Chinese))
[2] NAGRATH S, JANSEN K E, JR LAHEY R T. Computation of incompressible bubble dynamics with a stabilized finite element level set method[J]. Computer Methods in Applied Mechanics and Engineering, 2005, 194(42/44): 4565-4587.
[3] 马飞遥, 马逸尘, 沃维丰. 基于二重网格的定常Navier-Stokes方程的局部和并行有限元算法[J]. 应用数学和力学, 2007, 28(1): 25-33.(MA Feiyao, MA Yichen,WO Weifeng. Local and parallel finite element algorithms based on two-grid discretization for steady Navier-Stokes equations[J]. Applied Mathematics and Mechanics, 2007, 28(1): 25-33.(in Chinese))
[4] 骆艳, 冯民富. 可压缩Navier-Stokes方程的压力梯度局部投影间断有限元法[J]. 应用数学和力学, 2008, 29(2): 157-168.(LUO Yan, FENG Minfu. Discontinuous element pressure gradient stabilizations for the compressible Navier-Stokes equations based on local projections[J]. Applied Mathematics and Mechanics, 2008, 29(2): 157-168.(in Chinese))
[5] CHO M H, CHOI H G, YOO J Y. A direct reinitialization approach of level-set/splitting finite element method for simulating incompressible two-phase flows[J]. International Journal for Numerical Methods in Fluids, 2011, 67(11): 1637-1654.
[6] HEIMANN F, ENGWER C, IPPISCH O, et al. An unfitted interior penalty discontinuous Galerkin method for incompressible Navier-Stokes two-phase flow[J]. International Journal for Numerical Methods in Fluids, 2013, 71(3): 269-293.
[7] 章争荣. 不可压缩黏性流动N-S方程直接耦合数值求解的流形方法[C]//中国力学大会: 2013论文摘要集. 西安, 2013.(ZHANG Zhengrong. Manifold method for directly coupled numerical solution of N-S equations of incompressible viscous flow [C]//China Mechanics Conference: 2013 Abstracts. Xi’an, 2013.(in Chinese))
[8] 郭虹平, 欧阳洁. 气液两相流的间断有限元模拟[J]. 计算物理, 2015, 32(2): 160-168.(GUO Hongping, OUYANG Jie. Simulation of gas-liquid two-phase flows with discontinuous Galerkin method[J]. Chinese Journal of Computational Physics, 2015, 32(2): 160-168.(in Chinese))
[9] 秦望龙, 吕宏强, 伍贻兆, 等. 三维可压缩Navier-Stokes方程的间断Galerkin有限元方法研究[J]. 空气动力学学报, 2016, 34(5): 617-624.(QIN Wanglong, LÜ Hongqiang, WU Yizhao, et al. Discontinuous Galerkin method for 3-D compressible Navier-Stokes equations[J]. Acta Aerodynamica Sinica, 2016, 34(5): 617-624.(in Chinese))
[10] AMROUCHE C, ESCOBEDO M, GHOSH A. Semigroup theory for the Stokes operator with Navier boundary condition on spaces[J]. Archive for Rational Mechanics & Analysis, 2018, 223(2): 1-60.
[11] KIRK K, RHEBERGEN S. Analysis of a pressure-robust hybridized discontinuous Galerkin method for the stationary Navier-Stokes equations[J]. Journal of Scientific Computing, 2019, 81(2): 881-897.
[12] KARNIADAKIS G E, ISRAELI M, ORSZAG S A. High-order splitting methods for the incompressible Navier-Stokes equations[J]. Journal of Computational Physics, 1991, 97(2): 414-443.
[13] YACOUBI A E, XU S, WANG Z J. A new method for computing particle collisions in Navier-Stokes flows[J]. Journal of Computational Physics, 2019, 399: 108919.
[14] PAL S, HALOI R. On solution to the Navier-Stokes equations with Navier slip boundary condition for three dimensional incompressible fluid[J]. Acta Mathematica Scientia, 2019, 39(6): 1628-1638.
[15] SHAHBAZI K, FISCHER P F, ETHIER C R. A high-order discontinuous Galerkin method for the unsteady incompressible Navier-Stokes equations[J]. Journal of Computational Physics, 2007, 222(1): 391-407.
[16] JOHN V. Reference values for drag and lift of a two-dimensional time-dependent flow around a cylinder[J]. International Journal for Numerical Methods in Fluids, 2010, 44(7): 777-788.
陈亚飞(1992— ),女,硕士生(E-mail: 610556349@qq.com);
郑云英(1973— ),女,教授,博士(通讯作者. E-mail: zhengyunying@eyou.com).