基于有限点法的自由面流动的数值模拟

卢 雨，胡安康,2，刘亚冲

(1.哈尔滨工程大学 船舶工程学院，黑龙江 哈尔滨 150001；2.中集船舶海洋工程设计研究院，上海 201206)



基于有限点法的自由面流动的数值模拟

卢 雨1，胡安康1,2，刘亚冲1

(1.哈尔滨工程大学 船舶工程学院，黑龙江 哈尔滨 150001；2.中集船舶海洋工程设计研究院，上海 201206)

为了准确模拟复杂变形的自由面流动，应用空间粒子法中的有限点法(FPM)，即利用移动最小二乘(MLS)法近似流场函数，采用投影法分步求解压力泊松方程，迭代更新流场各粒子的物理信息，数值模拟大变形下的自由面流动问题.数值模拟结果表明：采用有限点法可以较准确地模拟出自由面流动过程中的复杂变形，如迸溅、翻滚、破碎以及入水等现象，与试验结果吻合较好.验证了采用的有限点法对于处理大变形自由面流动问题具有较高的可靠性及较好的灵活性.

拉格朗日观点；有限点法(FPM)；移动最小二乘法；自由面流动；溃坝流

自由面流动问题在许多工程计算中备受关注,然而流体的自由面判断是数值模拟中比较棘手的问题.众所周知,非线性Navier-Stokes方程是求解流体流动的控制方程.从流体力学的观点出发,描述流体的运动有2种方法：拉格朗日(Lagrangian)法和欧拉(Euler)法.欧拉观点借助于由固定空间位置点构建的网格来数值模拟流体的流动问题.Hirt等[1]基于欧拉观点使用VOF法对自由边界进行数值模拟.Renardy等[2]对VOF法进行改进,提出基于表面张力的抛物线型网格重建技术,使得自由边界的模拟更准确.不断变化的自由面不能与固定的网格节点位置时刻保持一致,因此自由面边界条件没有准确地应用于欧拉法中进行求解.

拉格朗日法着眼于流体质点,通过合适的边界条件,可以求解出这些控制点自身携带的物理属性(速度、压强等),进而可得这些流体质点的空间位置,故基于拉格朗日观点对自由面位置的确定便可迎刃而解.相比于传统基于欧拉观点的空间网格重构技术,拉格朗日法在处理大变形、复杂几何形状等问题上具有显著的优点,可以更精确地捕捉不断变化的边界面形状,并能够减小由大幅变形引起的数值误差.

基于拉格朗日观点求解流体流动问题主要采用空间粒子法.在问题域中及其边界布置粒子,排除了空间固定网格的拓扑关系,使用这些自由移动的粒子来表达流场信息,同时粒子法避免了因求解Navier-Stokes方程中的对流项离散而产生的数值耗散.目前,计算流体力学中的空间粒子法主要有：光滑粒子法(SPH)[3-7]、移动粒子半隐式法(MPS)[8-12]以及有限点法(FPM)[13-17]等.

本文基于FPM对自由面流动进行数值模拟研究.采用投影法对不可压缩Navier-Stokes方程进行求解,压力泊松方程中的空间导数由移动最小二乘法近似获得.通过对比FPM计算结果与试验数据,验证了拉格朗日观点下的有限点法在求解自由面流动问题中具有较高的可信度,得到了较满意的结论.

1 FPM数值方法

1.1 控制方程

在拉格朗日理论框架下,流体运动控制方程包括连续方程和动量守恒方程,对于不可压缩流体,可以分别写为

(1)

(2)

式中：ρ为流体密度；ui为直角坐标系中对应i轴的速度分量,相应地xi为对应i轴的方向矢量；p为压力；μ为流体的运动黏性系数；fi为源项(通常取重力加速度gi).由于是从拉格朗日观点出发,式(1)和(2)右端的时间导数项以物质导数的形式给出.

1.2 投影法

(3)

(4)

(5)

式(5)中增加了压力梯度项,通过引入不可压缩流体的连续方程,可以显式获得关于压力的泊松方程,即为

(6)

将式(5)沿边界Γ的外单位法向量n投影,得到关于压力p的Neumann边界条件：

(7)

由于边界面的不可穿透性,可得

(8)

1.3 移动最小二乘法函数近似

首先构造待求函数u(x)在控制点x处的支持域内(支持半径为h)的近似函数Πu(x),则可令

(9)

式中：bT(x)=[b1(x),b2(x),…,bk(x)],其中bi(x)为基函数,k为基函数的个数；aT(x)=[a1(x),a2(x),…,ak(x)],其中ai(x)为待定系数.本文计算中单项式基函数取二次基.

二维空间单项式基函数为

bT(x)=[1,x,y,x2,xy,y2];k=6.

(10)

三维空间单项式基函数为

bT(x)=[1,x,y,x2,xy,y2,yz,z2,xz];k=10.

(11)

欲求式(9)中的a(x),则转化为Πu(x)与u(x)的加权平方和J的极值问题.对J关于a(x)求偏导,则可得

A(x)a(x)=B(x)U.

(12)

式中：

(13)

B(x)=[w1(x)b(x),w2(x)b(x),…,wm(x)b(x)].

(14)

U为支持域内粒子点函数向量,即U=[u1(x),u2(x),…,um(x)].联合式(9)、(12),可得控制点x处的近似函数Πu(x)表达式：

Πu(x)=N(x)U=bT(x)A-1(x)B(x)U.

(15)

式中：N(x)为控制点x的支持域内形函数,即N(x)bT(x)A-1(x)B(x),则在该粒子点x处的待求函数是通过其支持域内各粒子的考察函数与形函数近似得到.在 FPM 计算中使用移动最小二乘法时所用到的权重函数w(x)常取距离权函数,具体表达式为

w(x)=

(16)

1.4 压力泊松方程导数模型

在利用投影法求解压力泊松方程的过程中,产生了压力函数的一阶导数及二阶导数项.根据移动最小二乘法的思想可知,流场中任意粒子点的待求函数均可由支持域内的周围粒子信息及形函数来表述.考虑通用性,给出任意函数的一阶及二阶导数模型.

考察流域(粒子总数为n)中任意粒子,将基函数及权重函数写成全局矩阵形式：

(17)

根据式(12),可令

S(x)a(x)=r(x).

(18)

式中：S(x)=VW(x)VT,r(x)=VW(x)u.对式(18)求一阶导数,可得

Sa′+S′a=r′.

(19)

于是

a′=S-1(r′-S′S-1r).

(20)

对式(18)求二阶导数,可得

a″=S-1[(r″-S″S-1r)-2S′S-1(r′-S′S-1r].

(21)根据式(20),可得移动最小二乘近似函数的一阶导数：

(Πu)′=(b′)T·a+bT·a′=

(b′)T·S-1r+bT·S-1(r′-S′S-1r).

(22)

联合式(20)、(21),可得移动最小二乘近似函数的二阶导数：

(Πu)″=(b″)T·a+2(b′)T·a′+bT·a″= (b″)T·S-1r+2(b′)T·S-1(r′-S′S-1r)+ bT·S-1((r″-S″S-1r)-2S′S-1(r′-S′S-1r).

(23)

借助式(22)、(23),可以离散压力泊松方程,最后形成关于压力的大型稀疏矩阵.采用较通用的双共轭梯度稳定法(bi-conjugate gradients stabilized method)[18]来迭代求解.为了满足Dirichlet条件,在计算开始时须指定每个粒子点的初始压力.

1.5 自由面判断

FPM 对自由液面的判断采用文献[19]的处理方式,依据粒子数密度的概念来区分自由液面与内部流体粒子,如图1所示.

图1 自由面粒子与流体内部粒子示意图Fig.1 Points clouds in free surface and flow domain

当粒子数密度〈ni〉满足：

〈ni〉<βn0.

(24)

即可判定为自由面粒子,其中n0为初始粒子数密度.在求解压力 Poisson 方程时,自由面粒子的压力赋值为 0.β是一参数,对自由面的判断有一定的影响,一般取值为β=0.80～0.99.

1.6 计算流程

图2 FPM计算流程图Fig.2 FPM computational procedure

针对上述FPM 的数值方法,给出具体实施的计算流程图,如图2所示.在数值计算过程中,由于采用移动最小二乘法近似空间场函数,在计算过程中涉及到求解逆矩阵的问题,会加大计算量,对计算效率产生了一定的影响.

2 FPM数值模拟算例

选用2个计算模型：二维溃坝模型、带障碍物三维溃坝模型.将数值计算结果与试验结果进行对比,分析FPM在模拟大变形自由面流动问题时的可靠性.

2.1 二维溃坝模型

溃坝问题作为无网格技术中的经典算例,常常被各粒子法用来模拟验证自由面流动的准确性.本文计算的二维溃坝模型的计算域为一敞篷正方形,边长为0.584 m.初始水柱高0.292 m,宽0.146 m.具体的空间位置见图3(a),实验模型如图 3(b)所示, 该实验是由 Koshizuka 等[20]完成的.计算初始时刻用细梁挡住水柱,随即迅速抽出右侧细梁,水柱在重力的作用下开始坍塌,计算中所有壁面采用不可滑移边界条件,水的密度为 1 000 kg/m3,运动黏性系数取10-6m2/s,重力加速度为 9.81 m/s2,粒子支持域半径取为0.007 m.由于流域粒子总数过少,会对计算的精度产生一定的影响,而过多则会耗费大量的计算时间.通过不断调整计算参数,综合考虑准确度与计算效率,根据相应的几何尺寸最终确定流体粒子总数为 5 846.

图3 溃坝模型几何尺寸及模型试验Fig.3 Dimension of dam break and model trial

图4给出1.0 s内不同时刻溃坝模型的数值模拟结果与试验现象.从数值计算结果来看,当t=0.2 s时,水柱已经坍塌, 且前缘沿着流域底部迅速向右侧流动,这一数值计算结果与实验现象基本吻合.当t= 0.3 s 时水柱撞击到右侧壁面,在惯性的驱使下继续沿壁面向上运动.当运动到 0.4 s 时上升的水柱产生飞溅的水花,且动量逐渐减小.直到t=0.6 s 左右,由于重力场的存在使得水流开始回落.当达到0.7 s 时,回落水柱与底部水流混合,水面发生翻卷,隆起水柱形似蘑菇状,且右侧壁面的水柱逐渐变薄.直到t=0.8 s 时,上升的水柱被全部拉回底部,并开始向左侧运动.在本文计算中自由面流动采用单相流数值模拟,不考虑空气,因此实际中的波面会将一定组分的空气裹入水中,使得水面变形更复杂,使得数值计算与试验现象存在一定的偏差,但整个运动的数值捕捉与试验现象吻合较好.总的来说,采用FPM可以较好地模拟出溃坝这种复杂大变形的自由面流动问题,证明有限点法具有较好的灵活性.

图5、6列出了不同时刻溃坝模型粒子速度(上)与压力(下)的数值计算结果.由图5、6可知,在溃坝流的初始阶段,靠近自由面的流体粒子具有较高的速度,远离自由面的粒子速度较小,在这一较大的速度梯度下,流场中静压力分布发生剧烈变化.随着流体运动趋于稳定,速度梯度逐渐减小使得静压力分布规则化.当t=0.4 s和t=0.8 s 时,均由于水柱的撞击作用,在流域右侧壁面与底部产生了较大的压力,图5、6给出这一现象的完美诠释.

为了从定量角度验证数值模拟结果的可靠性,将水柱前缘和侧壁处的高度值进行无量纲化,与试验数据进行对比,如图 7、8所示,试验数据来自 Martin 等[21]的工作.从结果来看,数值计算结果与实验基本吻合,说明FPM具有较高的可靠性.

2.2 带障碍物三维溃坝

为了进一步研究大变形溃坝流问题,计算采用三维溃坝模型且流动前方带有障碍物搁置,几何尺寸为长(h)×宽(L)×高(2h),流域尺寸:长(4L)×宽(L)×高(3L),其中L=0.146 m,h=0.024 m.初始流体模型见图9(a),试验模型[22]见图9(b).粒子支持域半径取值为0.015 m,其他计算参数的设置与算例1相同,流体粒子总数为 41 420.

图4 流场数值计算结果(t=0.1～0.8 s,△t=0.1 s)与试验结果(t=0.2～0.8 s,△t=0.2 s)的对比Fig.4 Flow numerical and experimental results of dam break problem at different time (FPM: t=0.1～0.8 s, △t=0.1 s; EXP: t=0.2～0.8 s, △t=0.2 s)

图5 溃坝模型在不同时刻(t=0.1,0.4,0.8 s)粒子速度场计算结果Fig.5 Computed velocity fields for dam break flow at different time(t=0.1,0.4,0.8 s)

图6 溃坝模型在不同时刻(t=0.1,0.4,0.8 s)粒子压力场计算结果Fig.6 Computed pressure fields for dam break flow at different time(t=0.1,0.4,0.8 s)

图7 溃坝流水柱波前的计算结果与实验结果对比Fig.7 Comparison between calculated position of leading edge with experimental data for dam break problem

图8 溃坝流水柱高度的计算结果与实验结果对比Fig.8 Comparison of computed height of collapsing water column with experimental results for dam break problem

图9 带障碍物三维溃坝模型几何尺寸及模型试验Fig.9 Three-dimension dam break with obstacle and model trial

不同时刻的数值计算结果与实验现象见图10、11.观察图10、11可知,当流体开始运动到0.1 s时,溃坝流已基本形成.当达到0.15 s左右时,流体遇到前方障碍物,与其发生强烈的砰击,随后隆起水柱呈“水舌”状并向右前方迸射,期间伴有溅落水花.当t=0.4 s时,迸射的水柱撞击到流域的右侧壁面,并开始向下坠落.在实际流动过程中,流体与壁面封闭的过程中会夹杂一定量的空气,而这一成分会参与并影响流体的运动,所以与单相流数值模拟结果存在一定的偏差,在t=0.5 s时刻数值计算与实验现象产生了较明显的差别,故采用两相流的数值模拟来改进这一缺陷.从整体的现象对比来看,采用FPM数值模拟带障碍物的三维溃坝给出了较真实的流动计算结果.

图12列出在0.6～3.0 s 之间不同时刻FPM的计算结果.当t=1.1 s左右时,从右侧翻滚而来的水流撞击到流域左侧壁面并伴有小距离的水柱升高,之后在重力作用下被卷及入水.当运动到1.8 s附近时,水流再一次砰击右壁面并被反弹,这一过程反复交替,使得流体动能逐渐减少,波面翻滚逐渐减弱,最终流动水面趋于平静.

图10 带障碍物三维溃坝模型在不同时刻(t=0.1,0.2,0.3,0.5 s)的数值模拟与试验结果对比Fig.10 Numerical and experimental results of three-dimension dam break with obstacle at different time (t=0.1,0.2,0.3,0.5 s)

图11 带障碍物三维溃坝模型在t=0.15 s和t=0.4 s时刻的数值模拟结果Fig.11 Numerical simulation of three-dimension dam break with obstacle at intermediate time(t=0.15, 0.4 s)

3 结 语

本文采用空间粒子法中的有限点法(FPM)对大变形的自由面流动进行数值模拟.流体流动的控制方程采用分步求解压力的投影法来求解,流场函数及空间导数采用移动最小二乘法来近似得到,空间离散的压力泊松方程利用较通用的双共轭梯度稳定法来求解.将数值计算结果与试验现象进行对比,发现两者吻合较好,验证了FPM的可靠性.从整个流体运动过程分析可知,采用FPM较准确地模拟出了水面的迸溅、翻滚、破碎以及入水等现象,可见拉格朗日理论框架下的有限点法在模拟大变形流动问题中体现出了极大的优势.相比于欧拉理论的网格技术,FPM具有更好的灵活性及较强的自由面模拟能力.有限点法是一种新型的数值算法,它的理论基础不如成熟的有限元法、有限体积法等完善,故在一些复杂边界条件的处理和大量工程实例面前,FPM的应用存在一定的局限性.对该算法的收敛性、稳定性、误差估计以及多相流等理论问题须进行更多、更深入的研究.

图12 带障碍物三维溃坝模型在后续时刻(t=0.6～3.0 s)的数值模拟结果Fig.12 Numerical results of three-dimension dam break with obstacle at subsequent times (t=0.6～3.0 s)

[1] HIRT C W, NICHOLS B D. Volume of fluid (VOF) method for the dynamics of free boundaries [J]. Journal of Computational Physics, 1981, 39(1): 201-225．

[2] RENARDY Y, RENARDY M. A parabolic reconstruction of surface tension for the volume-of-fluid method [J]. Journal of Computational Physics, 2002, 183(2): 400-421．

[3] CLEARY P W, PRAKASH M. Discrete-element modelling and smoothed particle hydrodynamics: potential in the environmental sciences [J]. The Royal Society A: Mathematical. Physical. Engineering Sciences, 2004, 362(1822): 2003-2030．

[4] ATA R, SOULAIMANI A. A stabilized SPH method for inviscid shallow water flows [J]. Numerical Methods in Fluids, 2005, 47(2): 139-159．

[5] FANG J N, PARRIAUX A, RENTSCHLER M, et al. Improved SPH methods for simulating free surface flows of viscous fluids [J]. Applied Numerical Mathematics, 2009, 59(2): 251-271．

[6] FANG J N, OWENS R G, TACHER L, et al. A numerical study of the SPH method for simulating transient viscoelastic free surface flows [J]. Journal of Non-Newtonian Fluid Mechanics, 2006, 139 (1/2): 68-84．

[7] ELLERO M, KROGER M, HESS S. Viscoelastic flows studied by smoothed particle dynamics [J]. Journal of Non-Newtonian Fluid Mechanics, 2002, 105 (1): 35-51．

[8] PAN X J, ZHANG H X, LU Y T. Numerical simulation of viscous liquid sloshing by moving-particle semi-implicit method [J]. Journal of Marine Science and Application, 2008, 7(3): 184-189．

[9] KHAYYER A, GOTOH H. Modified moving particle semi-implicit methods for the prediction of 2D wave impact pressure [J]. Coastal Engineering, 2009, 56(4): 419-440．

[10] PAN X J, ZHANG H X, LU Y T. Moving-particle semi-implicit method for vortex patterns and roll damping of 2D ship sections [J]. China Ocean Engineering, 2008, 22(3): 399-407．

[11] GOTOH H, SAKAI T. Key issues in the particle method for computation of wave breaking [J]. Coastal Engineering, 2006, 53(2/3): 171-179．

[12] GOTOH H, IKARI H, MEMITA T, et al. Lagrangian particle method for simulation of wave overtopping on a vertical seawall [J]. Coastal Engineering Journal, 2005, 47(2/3): 157-181．

[13] ONATE E, IDELSOHN S, ZIENKIEWICZ O, et al. A finite point method in computational mechanics: applications to convective transport and fluid flow [J]. International Journal for Numerical Methods in Engineering, 1996, 39(22): 3839-3866．

[14] ONATE E, SACCO S, IDELSOHN S. A finite point method for incompressible flow problems [J]. Computing and Visualization in Science, 2000, 3(1/2): 67-75．

[15] ONATE E, PERAZZO F, MIQUEL J. A finite point method for elasticity problems [J]. Computers and Structures, 2001, 79(22-25): 2151-2163．

[16] SHU C, DING H, CHEN H Q, et al. An upwind local RBF-DQ method for simulation of inviscid compressible flows [J]. Computer Methods in Applied Mechanics and Engineering, 2005, 194(18-20): 2001-2017．

[17] BOROOMAND B, NAJJAR M, ONATE E. The generalized finite point method [J]. Computational Mechanics, 2009, 44(2): 173-190．

[18] GU T X, ZUO X Y, LIU X P, et al. An improved parallel hybrid bi-conjugate gradient method suitable for distributed parallel computing [J]. Journal of Computational and Applied Mathematics, 2009, 226 (1): 55-65．

[19] 徐刚,段文洋. 自由面流动模拟的MPS算法研究[J]. 哈尔滨工程大学学报，2008, 29(6): 539-543. XU Gang, DUAN Wen-yang. An investigation of MPS algorithm for free surface flow simulation [J]. Journal of Harbin Engineering University, 2008, 29(6): 539-543．

[20] KOSHIZUKA S, OKA Y. Moving-particle semi-implicit method for fragmentation of incompressible fluid [J]. Nuclear Science and Engineering, 1996, 123(3): 421-434．

[21] MARTIN J C, MOYCE W J. An experimental study of the collapse of liquid columns on a rigid horizontal plane [J]. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 1952,244(882):312-324．[22] KOSHIZUKA S, OKA Y, TAMAKO H. A particle method for calculating splashing of incompressible viscous fluid [C]∥International Conference on Mathematics and Computations, Reactor Physics, and Environmental Analyses. Portland: [s.n.], 1995: 1514-1521.

Numerical simulation of free surface flow based on finite point method

LU Yu1, HU An-kang1,2, LIU Ya-chong1

(1.CollegeofShipbuildingEngineering,HarbinEngineeringUniversity,Harbin150001,China; 2.CIMCOceanEngineeringDesignandResearchInstituteLimitedCompany,Shanghai201206,China)

The finite point method (FPM) was employed to accurately simulate the free surface flow with complicated deformation. The projection method was adopted to fulfill the incompressibility of flow. The moving least squares (MLS) approach was introduced for solving the pressure Poisson equation and updating physics information of particles in flow filed over time steps. Then the flow with a large deformation of free surface was simulated. The numerical results show that the finite point method can accurately provide more detailed information and features of complicated free surface flow, such as the phenomena of water jump, rolled, merged and broken waves. The numerical results agreed well with the experimental data. The presented finite point method was proved to have excellent flexibility and reliability in dealing with complex free surface flows.

Lagrangian view; finite point method (FPM); moving least squares method; free surface flow; dam breaking flow

2014-12-30. 浙江大学学报(工学版)网址： www.journals.zju.edu.cn/eng

国家自然科学基金资助项目(51379040,51409063).

卢雨(1988-)，男，博士生，从事流体力学的研究. ORCID：0000-0001-7859-2876. E-mail：luyu90627@126.com 通信联系人：胡安康，女，教授. E-mail：ankang.hu@cimc.com

10.3785/j.issn.1008-973X.2016.01.009

U 661

A

1008-973X(2016)01-0055-08