基于修正光滑粒子流体动力学算法的低能量耗散数值波浪水池开发1)

2022-07-10黄晓婷孙鹏楠吕鸿冠钟诗蕴

力学学报 2022年6期

黄晓婷孙鹏楠 ,2) 吕鸿冠钟诗蕴

* (中山大学海洋工程与技术学院,广东珠海 519082)

† (南方海洋科学与工程广东省试验室(珠海),广东珠海 519082)

引言

随着海洋开发逐渐走向深远海,海洋环境也变得越发恶劣和复杂,准确预报结构物在极端海浪环境下的水动力性能是确保海洋工程装备安全性和可靠性的重要前提.由于试验手段存在成本高、试验周期长等缺点,数值模拟方法在波浪-结构相互作用、海岸水动力学等[1]应用上越发广泛.其中,建立高精度、高效率的数值水池是预报结构物水动力性能的重要手段之一.国内外学者采用了许多数值方法建立数值水池,开展波浪传播特性的研究.基于势流理论,Grilli 等[2]、Sung[3]、Baudic 等[4]采用高阶BEM 方法(拉格朗日网格法)对控制方程进行离散,构建数值水池.但无黏无旋的理想流体假设和波浪翻卷、破碎时的网格畸变使得强非线性波浪的动态演变过程较难模拟.对于欧拉网格法,如有限体积法(FVM) 等,文献[5] 采用开源程序库OpenFOAM构建了三维数值水池,生成的波浪精度较高,但在远距离传播时,受网格尺寸的限制,易造成数值耗散问题,且为了模拟自由液面的演化过程,常常需要结合复杂的VOF[6-7]等自由液面捕捉方法.

相比于网格化算法,近年来光滑粒子流体动力学(smoothed particle hydrodynamics,SPH)方法在船海领域得到广泛应用[8-10],尤其在波浪与结构物相互作用研究上优势突出[11-12].SPH 方法的无网格特性使其易于处理大变形问题,完全避免了网格畸变;凭借其拉格朗日粒子特性,能轻松实现粒子运动路径的追踪,有利于处理自由液面,方便地实现波浪翻卷、破碎等现象的模拟.但SPH 模拟中存在一定的数值耗散[13-15],常常引起能量的非物理性衰减.例如,在波浪远距离传播模拟过程中,能量非物理性耗散表现为波高的非物理性减小.如果波浪衰减不能得到有效补偿或改善,将会限制SPH 数值水池应用的物理尺度,无法实现远距离、长时间海浪环境演化的模拟.

目前,提高SPH 方法的能量守恒性,减少数值耗散,成为SPH 数值水池开发中广受关注的问题.研究人员主要从光滑长度及核函数选择上着手解决该问题.Colagrossi 等[16]的研究表明,选择较大的光滑长度系数可以减少模拟重力波时的数值耗散量,但此举容易造成较高的计算成本和较长的模拟时间.Zhang 等[17]采用不同核函数进行规则波传播的模拟,发现在SPH 模拟中,Quintic 核函数的动能衰减率更小.Macià等[18]的研究表明采用Wendland C2 核函数进行模拟可以有效减少动能衰减.Meng等[19]提出了一种TENO-SPH 模型与算法,提升了SPH 的数值精度.此外,加密粒子也能起到减少耗散的效果.

对于具有自由液面的流体运动模拟,自由液面附近粒子的支持域被截断,导致物理量及其导数的计算出现较大误差.即使是处于流体内部的粒子,当其分布不均匀时,往往也无法保证理论上的二阶精度.这些误差是传统SPH 格式数值耗散的主要来源之一[20].因此,在SPH 算法中,施加核函数修正技术[21]也是一种减少能量耗散的重要手段.Wen 等[22]指出,采用核函数梯度修正(kernel gradient correction,KGC)可以显著降低波浪传播的衰减,而无需增大光滑长度,但缺点是粒子间的相互作用力不再具有对称性,从而无法保证动量守恒.

受文献[23]启发,采用对称化的修正核函数导数格式,对动量方程进行离散,有利于减少数值耗散,但文献[23]的算法需要搜索自由液面,这将带来额外的计算负担.为此,本文基于δ-SPH 算法,开发一种新的SPH 数值波浪水池,采用另一种核函数梯度修正方式,减少SPH 数值波浪水池中能量非物理性耗散的同时,确保SPH 算法格式简单、计算稳定、动量守恒.

本文具体安排如下:数值方法中,首先介绍δ-SPH 模型的控制方程,随后论述本文修正算法的推导过程.数值结果中,首先通过振荡液滴基准算例,验证和讨论本修正算法对非物理性能量耗散现象的改善效果.随后,在数值波浪水池的构建上,对规则波与不规则波传播进行模拟.先通过数值波高与试验波高的对比,论证传统δ-SPH 数值水池中的波浪衰减现象;随后,通过传统算法结果与本文修正算法结果的对比,说明本文修正算法在低能量耗散数值波浪水池开发上的有效性.

1 数值算法

1.1 δ-SPH 模型

SPH 方法是一种具有拉格朗日特性的粒子法.在计算时,流场被离散成粒子,通过赋予粒子物理属性,追踪粒子的拉格朗日运动来实现流场物理演化过程的模拟.在δ-SPH 计算模型中,粒子密度、速度和位置变化的控制方程[24]为

式中,n=2,光滑长度h=αs∆x,αs为光滑长度系数,∆x为粒子间距.α 为黏性系数,在SPH 模拟中,通过在消波区内增大黏性系数实现数值消波,具体设置为

式中,L为数值水池总长度,L−Ldamp为消波区长度.在本文模拟中,采用两种常用的光滑长度系数αs=1.3,2.0 进行计算.其中 αs=1.3 时,由于核函数紧支域内粒子数量更少,因此计算量小,但是离散误差相比 αs=2.0 时更大,数值耗散也会更大.本文将开发改进SPH 算法,力争在αs=1.3 时,也获得较好的计算精度和较小的数值耗散.

1.2 低能量耗散的修正δ-SPHC 算法

在传统SPH 方法中,梯度算子粒子近似[25]可写为

然而,在粒子分布不均匀或核函数截断情况下,该离散形式不能较好地确保计算精度.对此,许多修正方法被提出,其中常用的有修正光滑粒子法(corrective smoothed particle method,CSPM)[26-27].梯度算子采用CSPM 修正后为

经Wen 等[22]和Zago 等[23]研究启发,采用上述CSPM 算法中修正的核函数梯度,有助于改善SPH 计算中能量的衰减问题,但是为了保证动量守恒,也就是确保粒子对(i↔j)之间相互作用力的对称性,需要对CSPM 格式进行改进.对此,Zago 等[23]提出了一种改进后的修正格式,但是需要显式搜索自由液面,这带来了额外的计算负担,而本文基于CSPM 修正形式,参考Sun 等[28]处理梯度算子的方法,采用另一种具有对称性的梯度算子离散形式,推导过程见下:

首先,根据核函数导数的特性[25]有

本文修正算法记为δ-SPHC.对比Zago 等[23]的算法,本修正算法避免了自由液面搜索的步骤;另外,相对于δ-SPH,δ-SPHC修正算法并不引入新的核函数修正矩阵,只是再次调用密度耗散项中的矩阵Mi对动量方程中压力梯度算子进行修正,因此不会造成计算成本和时间的显著增加.下文将采用振荡液滴算例进行修正算法抗能量非物理性衰减效果的测试和讨论,并应用于数值波浪水池规则波与不规则波的数值模拟和试验验证.本文的数值积分采用四阶龙格库塔法,具体可参考文献[31].

2 振荡液滴算例验证

2.1 收敛性分析与精度验证

振荡液滴算例无需施加固壁边界,初始条件简单,被广泛应用于SPH 等数值新算法的验证[32].Antuono 等[33]采用振荡液滴算例验证了δ-SPH 模型中的能量守恒问题.基于此,本文采用该算例对修正算法δ-SPHC减少能量非物理性衰减的效果进行验证.

该算例中,初始速度V(u,v)、初始半径为R的圆形液滴在向心体积力场F中振荡,其中F(x,y)=−B2r(x,y),B=1,R=1 .图1 展示了初始时刻t0时液滴的速度场分布[33],可表示为

图1 振荡液滴初始速度分布Fig.1 Initial velocity distribution

A(t0) 为控制液滴初始速度场的参数.在本文模拟中,A(t0)=1 .初始压力场分布如图2 所示,压力场[33]表示为

图2 振荡液滴初始压力分布Fig.2 Initial pressure distribution

首先,为检验SPH 修正算法的收敛性,分别在粒子分辨率为R/∆x=50,R/∆x=100,R/∆x=200 时,开展振荡液滴运动特性研究.由于振荡液滴质心与向心体积力F的汇聚点重合,液滴不受其他外力作用,液滴在原位保持振荡.因此,理论上振荡液滴的角动量、线动量都守恒,即

图3 展示在光滑长度系数 αs=1.3 和 αs=2.0 条件下,不同粒子分辨率时,角动量MA随时间的变化情况.可见,在三种粒子分辨率下,MA的值围绕理论值MA=0 波动的量级均在10−5以内,且随着粒子分辨率的增加,角动量MA的误差范围逐渐变小.在粒子分辨率R/∆x=200 时,MA趋近于零.即在R/∆x=200时,角动量误差逐渐收敛于0,说明该修正算法在提高粒子分辨率时能实现良好的角动量守恒性.

图3 不同粒子分辨率下角动量 MA 的时历曲线Fig.3 Time evolution of the angular momentum MA with different particle resolutions

图4 展示了光滑长度系数 αs=1.3 和 αs=2.0 条件下,R/∆x=200 分辨率时,线动量ML随时间的变化情况.观察可得,线动量误差在−10−14～ 10−14范围内波动,可以认为本修正算法的线动量实现了精确守恒.

图4 线动量时历曲线Fig.4 Time evolution of the linear momentum

同时,定义液滴振荡过程的机械能衰减率为Eper(t)=(E(t)−Eint)/Eint,其中Eint为初始机械能,E(t) 为t时刻的机械能.液滴振荡过程的机械能衰减率随时间变化结果如图5 所示.可以发现,随着粒子分辨率的提高,机械能衰减率Eper越来越小.在R/∆x=100,R/∆x=200 时,两者机械能衰减率基本接近,计算时长达t=120 时,机械能衰减控制在5%以内.

图5 αs=2.0,不同粒子分辨率时机械能衰减率 Eper 时历曲线Fig.5 Time evolution of the decay rate of mechanical energy at different particle resolutions with αs=2.0

图6 展示了R/∆x=200 时,振荡液滴在不同时刻的理论形态与修正算法δ-SPHC模拟的形态.可以发现,在两种光滑长度条件下,δ-SPHC模拟的振荡形态均与理论值吻合良好,且在t=8.75T周期后,液滴仍保持较为光滑的形态.由此可见,本修正算法具有较高的精度.

图6 αs=1.3(上)和αs=2.0(下)时,振荡液滴形态:理论值(虚线)和SPH 结果(压力云图)Fig.6 With smoothing length coefficient αs=1.3 (top) and αs=2.0 (bottom),the shapes of oscillating droplet:theoretical results (dash line) and SPH results (pressure contour)

综上可见,修正SPH 算法在粒子分辨率R/∆x=200时,具有较好的数值精度和动量与能量守恒性.因此,下文将在粒子分辨率R/∆x=200 时,依托该算例进一步讨论不同算法抗能量衰减效果的差异.

2.2 传统δ-SPH 与修正δ-SPHC 算法的对比研究

基于上节收敛性分析与精度验证,本节在粒子分辨率R/∆x=200 条件下,开展传统δ-SPH 模型与改进δ-SPHC模型抗能量衰减效果的对比研究.图7展示了δ-SPH 与δ-SPHC模拟的振荡液滴动能随时间的变化曲线.可见,在传统算法δ-SPH 结果中,αs=1.3 时,动能随着时间的衰减十分明显.在αs=2.0时,衰减有所减缓.但是,取h=2∆x相比h=1.3∆x,由于相邻粒子数量的增多,会导致额外的计算成本和时间,且从图7 的局部放大图可以看出,在较长时间的计算下,仍会出现轻微动能衰减.对比修正算法δ-SPHC,在光滑长度系数αs=1.3 时,动能幅值基本保持不变,已体现出较好的抗能量衰减效果;在光滑长度系数αs=2.0 时,本修正算法实现了良好的能量守恒.此外,由局部放大图可以看出,在修正算法αs=2.0 时,动能幅值与初始幅值始终保持一致,取得了较为精确的计算结果.而在δ-SPHC结果中,αs=1.3时,动能第一个峰值较初始值高,主要原因在第一个周期的计算中,粒子会从格子分布状态(图8(a))突然过渡到一种各项同性的均匀分布状态(图8(b)).在粒子分布变化过程中,存在一些粒子的瞬时分布十分不均匀,此时,这不均匀性会造成.也就是说,部分粒子的运动并不是由于压力梯度造成,而是粒子分布的不均匀性导致的数值加速度产生,这对粒子的总动能造成了一定的非物理性影响.但是,从图7 下方的放大图可以看出,这种由于粒子分布突变引起的动能误差控制在1%以内.类似现象在Colagrossi 等[34]的研究结果也有过报道,可通过改善初始粒子分布得以解决.

图7 δ-SPH 与δ-SPHC 动能比较:总体图(上)与局部放大图(下),其中黑矩形框为放大区域Fig.7 The comparison between time evolutions of kinetic energy between δ-SPH (top) and δ-SPHC (bottom),an enlarged zone view of the results in the dark rectangle is shown in the bottom panel

图8 振荡液滴粒子分布Fig.8 Particle distribution of oscillating droplet

为进一步突显修正算法的低能量耗散特性,图9呈现了两种算法中,机械能衰减率的时历曲线.δ-SPHC算法在两种光滑长度系数下的计算结果中,机械能衰减率均较小,基本接近零,说明本文修正算法能够有效减少机械能的非物理性耗散.其中,在光滑长度系数 αs=2.0 时效果更为显著.值得一提的是,δ-SPHC算法在αs=2.0 时的机械能衰减率相比δ-SPH 算法在αs=1.3 时的机械能衰减率更小,表明本修正SPH 算法的一大优势是能在较小光滑长度系数下实现较为准确的模拟,从而显著减少计算时间.

图9 不同算法和光滑长度系数下,振荡液滴机械能衰减率时历曲线Fig.9 Time evolution of the decay rate of mechanical energy for the oscillating droplet with different SPH models and different smoothing length coefficient

以上算例在配有Intel(R) Core(TM)I9-10900 K CPU 和48 GB 内存的个人电脑上进行计算.

表1 比较了在最高粒子分辨率R/∆x=200,αs=2.0 参数下,δ-SPH 与δ-SPHC模型模拟振荡液滴的计算时长.可见,两种算法单步计算时间相当.具体而言,虽然修正算法δ-SPHC对压力梯度项进行了核函数修正,但相对δ-SPH 算法,计算时间仅增加1.09 倍,计算量增加并不显著.在下节中,该算法将被应用于数值波浪水池,验证其改善SPH 模拟远距离波浪传播时的波高抗衰减效果.

表1 δ-SPH 算法与δ-SPHC 算法计算效率比较Table 1 Comparison of computational efficiency between δ-SPH and δ-SPHC schemes

3 SPH 数值造波与试验验证

3.1 试验与数值水池设置

3.1.1 试验水池

作者在中山大学海洋工程与技术学院波浪水槽中开展造波试验.如图10 所示,该水池长度L=16 m,水深d=0.266 m,采用推板造波方式生成规则波和不规则波,推板运动曲线如图11 所示.试验时,波高仪设置在距离造波板平衡位置xprobe=2.37 m,4.37 m,5.37 m,6.37 m,15 m 处.

图10 中山大学波浪试验水槽Fig.10 Experimental wave tank in Sun Yat-sen University

图11 造波推板横向位移随时间变化:(a)周期为T =1 s 的规则波,(b)不规则波Fig.11 Paddle motions of the regular wave with (a) period T =1 s and(b) irregular wave

3.1.2 SPH 数值水池设置

在SPH 模拟中,固定虚粒子边界法[35,36]能实现固壁边界的精确处理,有效防止壁面穿透,且保证壁面处压力连续,因而得到了广泛应用.

本文SPH 数值模拟中,数值水池长度L=60 m,边界条件采用固定虚粒子进行施加.如图12 所示,固定虚粒子厚度不小于核函数半径,其压力从流体粒子外插得到,更多细节可参考文献[37].本文SPH模拟中,通过强制流体粒子与壁面虚粒子之间黏性力为0,实现壁面边界处的自由滑移边界条件.数值水池左侧的运动造波板也由虚粒子组成,通过控制造波板虚粒子的总体位置发生改变,实现推板造波.具体实现方式如下.

图12 SPH 数值水池示意图,红竖线表示波高仪,选取距离波高仪一个粒子间距内自由面粒子的最大高度为波面高度Fig.12 Schematic diagram of the numerical wave tank,the red vertical line represents wave gage.Maximum height of the free-surface particle within one-particle distance from the red line is measured to represent the wave elevation

强制SPH 造波板按照图11 中试验所输出的造波板位移曲线进行运动,从而确保SPH 模拟中,造波板运动路径与试验一致.同时,为克服SPH 模拟的时间步长与试验中输出的造波板运动数据时间间隔不一致,采用线性插值技术,根据物理试验中造波板运动特性,更新SPH 模拟中的造波板位移xS PH、速度vS PH及加速度aS PH,即

其中,上标SPH表示数值水池模拟中的物理量,上标k表示试验中造波板运动输出的时间步.试验中造波板的速度vk与加速度ak均采用线性差分计算获得,即

此外,值得注意的是,当波浪传播到数值水池另一端时,往往会产生波浪反射,降低造波精度,这就需要在数值水池右端构建数值消波区.类似于网格算法中增大数值黏性来实现阻尼消波,本文通过提高SPH 数值水池消波区的黏性系数,使波浪进入消波区时能够在短时间内衰减.具体设置为 :在数值水池设置消波区Ldamp,黏性系数 α 线性增大,详见公式(4).本文计算中,Ldamp=20 m.

3.2 δ-SPH 算法在高粒子分辨率下的精度验证

Colagrossi 等[16]指出δ-SPH 模型能较为准确地重现在物理黏性条件下的波浪传播和衰减过程.因此,本文基于δ-SPH 模型开发低能量耗散的数值水池.首先,在αs=2.0 条件下,验证δ-SPH 计算结果的精度和收敛性.选取d/∆x=15,30,60,即水深方向上粒子数分别为15,30,60 进行规则波模拟,在三种粒子分辨率条件下监测xprobe=2.37 处的波面高度随时间变化,如图13 所示.可以看出,数值波高结果变化趋势均与试验波高结果吻合良好,在d/∆x=60 时,数值波高收敛于试验波高.基于以上讨论,在d/∆x=60 条件下,δ-SPH 模型能较为精确实现造波及波浪传播的模拟.但是值得一提的是,δ-SPH 模型在αs<2.0 时,容易出现波高的非物理性衰减,具体论述见下一节.

图13 αs=2.0 时,不同粒子分辨率下波面高度的δ-SPH 模拟值与试验值对比Fig.13 Comparison of the wave elevations between δ-SPH results and experimental data,SPH results with different particle resolutions and αs=2.0 are provided

3.3 δ-SPH 算法中的波浪衰减现象

本节采用δ-SPH 模型,在光滑长度系数 αs=1.3 和 αs=2.0 条件下,进行规则波造波模拟.为了测定波高变化曲线,沿水池长度方向,距离造波板初始位置xprobe=2.37 m,4.37 m,5.37 m,6.37 m,15 m 处设置5 个数值波高仪.

在αs=1.3 时,图14 展示了计算时间为t=27 s 的波面形态,可以较为明显地看出波浪在x=10 m 处开始出现衰减.选取xprobe=2.37 m,5.37 m,15 m 处的波面高度时历曲线进行比较,如图15 所示.在xprobe=2.37 m 处,波高与试验相差较大,且随着传播距离的增大,波高逐渐衰减,在xprobe=15 m 处,波高衰减严重.可以看出,在传统δ-SPH 中,光滑长度系数 αs=1.3 条件下无法准确模拟波浪传播.图16 展示了δ-SPH 算法在光滑长度系数 αs=2.0 时计算的波面形态.与振荡液滴类似,相比光滑长度系数 αs=1.3 时,衰减速度有所减缓,与Colagrossi 等[16]的结论相符.这表明,增大 αs可以有效减缓SPH 的数值衰减,但在大尺度、长时间、远距离的海浪模拟下,光滑长度系数 αs取2.0 时,将带来巨大的计算量,不利于工程应用.且比较不同距离处波高,如图17所示,可以发现在xprobe=15 m 时,仍存在明显波高衰减.因此,改善SPH 方法在数值波浪模拟时能量的非物理性耗散,在较低光滑长度系数下实现数值水池造波的高精度高效率模拟具有重要意义.

图14 在αs=1.3,d/∆x=60 条件下,传统δ-SPH 模拟的波面形态Fig.14 Wave surface simulated by δ-SPH method with αs=1.3 and d/∆x=60

图15 在αs=1.3 ,d/∆x=60 条件下,传统δ-SPH 模型计算的不同位置波面高度时历曲线Fig.15 Wave elevation probed at different positions simulated by δ-SPH method with αs=1.3 and d/∆x=60

图16 在αs=2.0 ,d/∆x=60 条件下,传统δ-SPH 模拟的波面形态Fig.16 Wave surface simulated by δ-SPH method with αs=2.0 and d/∆x=60

图17 在αs=2.0,d/∆x=60 条件下,传统δ-SPH 模拟的不同位置波面高度时历曲线Fig.17 Wave elevation probed at different positions simulated by δ-SPH method with αs=2.0 and d/∆x=60

4 修正算法δ-SPHC 的波浪传播模拟

4.1 规则波模拟

4.1.1 收敛性分析与精度验证

首先,基于规则波传播算例,校核本文修正的δ-SPHC算法精度与收敛性.如图18 所示,在αs=1.3时,选取xprobe=2.37 m 处数值波高结果与试验结果进行比较,可以看出,水深方向不同数量粒子的计算结果略有差别,且随着d/∆x的增大,数值波高趋近于试验波高.在d/∆x=60 时,SPH 结果与试验结果吻合良好.与振荡液滴算例类似,图18 方框中的峰值可能是由于粒子由初始的格子状分布到各向同性均匀分布时引起的能量瞬时误差.试验波高曲线与SPH 数值结果曲线在初始阶段第一个峰值存在相位差,可能是由于试验波高采集与造波板运动存在一定的时间差.在αs=2.0 时,水深方向不同粒子数量的计算结果与试验结果对比如图19 所示.同样可以看到,在d/∆x=60 时,数值波高曲线与试验结果基本吻合.另外,由于物理试验水池长度仅有16 m,消波区长度设置有限,易导致试验后期消波不彻底,造成波浪周期的变化,即随着时间增加,周期略小于1 s;而数值水池长度L=60 m,能较为准确地重现规则波传播,并且确保周期保持不变.因此,可以看到,在图18、图19 中,t=20 s 之后SPH 结果的波峰和波谷时刻与试验值存在略微误差,但是波高和波谷的幅值与试验值吻合良好.基于以上讨论可以得出,在水深方向粒子数d/∆x=60 条件下,计算结果能够达到预期精度要求.因此,下文将在d/∆x=60 粒子设置下,进行修正δ-SPHC算法的抗能量非物理性衰减效果的讨论.

图18 在αs=1.3 和 d/∆x=15,30,60 条件下,δ-SPHC 计算的波面高度时历曲线与试验值对比Fig.18 Comparison between δ-SPHC results and experimental data for wave elevation:SPH results with αs=1.3 and d/∆x=15,30,60

图19 在αs=2.0 和 d/∆x=15,30,60 条件下,δ-SPHC 计算的波面高度时历曲线与试验值对比Fig.19 Comparison between δ-SPHC results and experimental data for wave elevation:SPH results with αs=2.0 and d/∆x=15,30,60

图20 呈现了在xprobe=2.37 m 处,规则波波高数值模拟结果与试验结果的对比.采用传统δ-SPH 计算时,光滑长度系数 αs=1.3 条件下,由于核函数紧支域内粒子数量少,数值结果与试验结果相差较大,而光滑长度系数 αs=2.0 时,数值波高结果接近试验结果.但值得注意的是,采用δ-SPHC,光滑长度系数αs=1.3 时,数值结果与试验结果能较好吻合,基本达到传统δ-SPH 算法在光滑长度系数 αs=2.0 时的计算精度,与上文振荡液滴算例得出的结论一致,再次说明了本文修正δ-SPHC算法的优势,即在光滑长度系数 αs=1.3,2.0 时,计算结果与试验结果均能较好吻合.图21 展示了数值水池的压力场,可见流场压力分布均匀,在运动较为剧烈的推板附近,也能确保压力场稳定.基于以上讨论,可以认为修正的δ-SPHC算法具有较高的精度.

图20 xprobe=2.37 m 处,不同SPH 算法下的波面高度数值结果与试验值对比Fig.20 Comparison between different SPH simulations and experimental data for wave elevations at xprobe=2.37 m

图21 规则波压力云图(δ-SPHC,αs=2.0,d/∆x=60)Fig.21 Pressure contour of regular wave (δ-SPHC,αs=2.0,d/∆x=60)

4.1.2 δ-SPHC的抗能量非物理性衰减结果分析

经过前期收敛性与精度验证,本节在水深方向粒子数取d/∆x=60 时,开展修正算法δ-SPHC在两种光滑长度系数 αs=1.3,2.0 下的抗能量非物理性衰减效果讨论.在αs=1.3 时,如图22 所示,观察在修正算法模拟下的自由面形态,可以看出在波浪传播20 m 时,波面仍与推板附近波面基本持平.监测不同距离下的波面高度如图23 所示,与试验结果对比可以看出,在xprobe=15 m 处,δ-SPHC的数值波高与推板附近波高基本接近,说明能量衰减现象得到极大改善.图24 展示了在αs=2.0 时,规则波传播的波面形态.可见,在远距离传播时,波面仍保持在同一高度,说明了修正算法抗能量衰减的有效性.对比图23 和图25,可见在xprobe=15 m 位置处,αs=2.0 比 αs=1.3 计算的波高更加精确和稳定,说明光滑长度系数的增大有利于减少离散误差,其中在xprobe=15 m 处出现略微振荡,原因主要是由于数值误差在远距离、长时间上的累积造成,但是相比传统算法,修正算法的精度得到了有效改善.

图22 αs=1.3,d/∆x=60 时,δ-SPHC 算法模拟的波面形态Fig.22 The wave surface simulated by δ-SPHC with αs=1.3 and d/∆x=60

图23 αs=1.3,d/∆x=60 时,δ-SPHC 模拟的不同位置波面高度时历曲线Fig.23 Wave elevations probed at different positions simulated by δ-SPHC method with αs=1.3 and d/∆x=60

图24 αs=2.0,d/∆x=60 时,δ-SPHC 算法模拟的波面形态Fig.24 The wave surface simulated by δ-SPHC with αs=2.0 and

图25 αs=2.0,d/∆x=60 时,δ-SPHC 模拟的不同位置波面高度时历曲线Fig.25 Wave elevations probed at different positions simulated by δ-SPHC method with αs=2.0 and d/∆x=60

如图26 所示,综合对比δ-SPH,δ-SPHC与试验方法在xprobe=6.37 m 的波高曲线,可见相比于传统δ-SPH 模型,修正的δ-SPHC算法在αs=1.3 时的波高结果更加准确.值得一提的是,其计算结果与传统δ-SPH 模型在αs=2.0 时波高结果基本重合,并与试验波高吻合良好.此外,采用δ-SPHC修正算法,在光滑长度系数 αs=2.0 时,波高计算结果与试验波高高度吻合.

图26 在xprobe=6.37 m 处,规则波波面高度数值模拟结果与试验结果Fig.26 Time evolution of the regular wave elevation at xprobe=6.37 m:SPH results and experimental data

为进一步定量验证δ-SPHC的模拟精度,表2 给出了试验和SPH 模拟中,xprobe=6.37 m 测点处的平均波高值及相对误差,后者定义为 ε=|HSPH−HEXP|/HEXP,HEXP为试验测量的平均波高.可见,在两种光滑长度系数下,δ-SPHC算法在较远处的波高监测结果与试验相对误差均较小.本节讨论表明,修正算法能有效解决能量非物理性衰减问题,且在较低光滑长度系数时能得到更为准确的结果.

4.1.3 不同人工黏性系数下δ-SPHC模拟结果分析

在数值水池中,能量衰减来源主要有两部分:一是由于流体物理黏性引起的能量耗散;二是本文所聚焦的由于数值误差引起的能量非物理性耗散.本文虽然未在SPH 控制方程中添加物理黏性项,但是使用了人工黏性项πiv,这等效于增加了流体的黏性.为进一步讨论 δ -SPHC在不同人工黏性系数 α 下,对长时间、远距离波浪传播模拟的能量耗散情况,选取四个不同黏性系数,即 α=0.01,0.02,0.05,0.1 开展数值模拟.图27 给出了在xprobe=6.37 m 处监测的波高曲线.可见,在α=0.01 和 α=0.02 时,波浪的衰减较小,与实验值吻合较好.在α=0.05 时,由于采用了更大的人工黏性,等效增加了物理黏性,因此波高衰减有所增加.在α=0.1 时,由于人工黏性引起的波高衰减更为明显,与实验值相比误差较大.综合以上分析,在α=0.02 时,能够获得较为稳定和精确的波高计算结果,且与试验值吻合良好.因此,本文其他算例均基于 α=0.02 进行模拟.基于此,在下文中进一步开展不规则波生成和传播的SPH 模拟和验证.

图27 δ-SPHC 在不同黏性系数α 时,xprobe=6.37 m 处波面高度时历曲线Fig.27 Time history of wave elevation at xprobe=6.37 m simulated by δ-SPHC with different α

4.2 基于δ-SPHC 的不规则波模拟

4.2.1 收敛性分析与精度验证

经过上文分析,在光滑长度系数 αs=1.3 时,修正算法δ-SPHC能够准确模拟规则波传播,而本节算法,在αs=1.3 时,模拟和验证不规则波传播.选取d/∆x=15,30,60 进行收敛性分析,如图28 所示,在水深方向粒子数d/∆x=15,30,60 时,在xprobe=2.37处探测的波高与试验值基本一致.此外,在不规则波模拟中,δ-SPHC实现了稳定的压力场预报(见图29),进一步说明了δ-SPHC方法所构建的数值水池的精度.由于d/∆x=60 时,数值结果已基本收敛于试验值,因此下文选取d/∆x=60 进行不规则波模拟结果的讨论.

图28 αs=1.3 和不同粒子分辨率条件下,不规则波δ-SPHC 模拟结果的收敛性分析Fig.28 Convergence analysis of irregular waves simulated by δ-SPHC with αs=1.3 at different particle resolutions

图29 不规则波压力云图(δ-SPHC,αs=2.0,d/∆x=60)Fig.29 Pressure contour of irregular wave (δ-SPHC,αs=2.0,d/∆x=60)

4.2.2 δ-SPHC抗能量非物理性衰减结果分析

如图30 所示,监测xprobe=6.37 处波高时历曲线,可见在传统δ-SPH 计算结果中,与试验结果相比,波高衰减明显;而修正算法结果与试验结果基本吻合,说明不规则波模拟中波高衰减现象也得到较好改善.在δ-SPHC计算框架下,数值水池生成不规则波的能量非物理性耗散问题得到解决,并能够在光滑长度系数 αs=1.3 时,实现较高的造波精度.