DEVELOPMENT OF A NUMERICAL WAVE TANK WITH A CORRECTED SMOOTHED PARTICLE HYDRODYNAMICS SCHEME TO REDUCE NONPHYSICAL ENERGY DISSIPATION
-
摘要: 目前, 无网格光滑粒子流体动力学SPH粒子法在波浪与结构物相互作用研究方面得到广泛应用, 但该方法模拟波浪远距离传播时, 常常面临严重的能量耗散问题, 导致波高非物理性降低, 给大范围海域、长时间作用下的波-物耦合作用研究带来一定困难. 对此, 本文采用一种核函数修正算法, 在确保粒子间相互作用对称性的同时, 改进压力梯度离散项的计算精度, 设法解决SPH方法中能量非物理性耗散的难题. 相较于前人减缓能量非物理性衰减的方法, 本文的修正SPH算法避免了自由液面搜索等复杂处理过程, 并能保证动量守恒特性. 数值结果中, 采用振荡液滴、规则波、不规则波等算例, 验证本修正SPH算法的准确性和有效性. 结果表明, 该修正SPH算法能准确模拟振荡液滴形态变化, 且动能保持较好守恒性. 通过数值水池与物理水池两者规则波与不规则波结果的对比分析表明, 基于本文修正SPH算法建立的数值波浪水池具有较好的抗能量衰减效果, 能实现长时间、远距离波浪传播的准确模拟. 此外, 本算法能在低光滑长度系数条件下, 实现精确模拟, 将极大缩减三维SPH模拟的时间, 从而节约计算成本.Abstract: So far, the smoothed particle hydrodynamics (SPH) method has been widely applied in the study of the interactions between water waves and structures. However, nonphysical energy dissipation is a still problem which challenges the simulation of wave-body interactions at large-scale and long duration. For example, in the SPH simulation of wave propagation to a long distance, the wave height could gradually become much smaller than the one generated near the wave maker. To tackle this problem, in this work a kernel correction algorithm is applied to the pressure gradient term in the SPH model, aiming to prevent nonphysical energy dissipation in long time simulations. The kernel correction algorithm is able to ensure the symmetry of the interaction between particle pairs, and therefore, compared with other corrective methods, the present corrected algorithm ensures the conservation of linear momentum and also avoids the complicated treatment at the free surface. Two numerical cases, i.e., the oscillating droplet and wave propagation in a numerical wave tank, are presented to test the accuracy and validity of present corrected SPH algorithm. For the oscillating droplet case, the corrected algorithm is shown to accurately simulate the evolution of the droplet shape, and the kinetic energy is dissipated much slower than traditional SPH models. Through the simulations of regular and irregular wave propagations as well as validations with experimental data, the capability of the corrected SPH algorithm to reduce nonphysical energy attenuation is demonstrated, even for wave propagation at long-term and long-distance conditions. In addition, this algorithm will be shown to be optimal for the SPH simulation at small smooth length, which contributes to save SPH computational cost significantly at three dimensional simulations.
-
引 言
随着海洋开发逐渐走向深远海, 海洋环境也变得越发恶劣和复杂, 准确预报结构物在极端海浪环境下的水动力性能是确保海洋工程装备安全性和可靠性的重要前提. 由于试验手段存在成本高、试验周期长等缺点, 数值模拟方法在波浪-结构相互作用、海岸水动力学等[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]为
$$ \left. \begin{array}{l} \dfrac{{{\text{D}}{\rho _i}}}{{{\text{D}}t}} = - {\rho _i}\displaystyle\mathop \sum \limits_j ({{\boldsymbol{u}}_j} - {{\boldsymbol{u}}_i}){\nabla _i}{W_{{{ij}}}}{V_j} + {D_i}^\rho \hfill \\ \dfrac{{{\text{D}}{{\boldsymbol{u}}_i}}}{{{\text{Dt}}}} = - \dfrac{1}{{{\rho _i}}} {\left\langle {\nabla p} \right\rangle }^C_i + {\boldsymbol{g }}+ {{\textit{π}} _i}^v \hfill \\ \dfrac{{{\text{D}}{{\boldsymbol{r}}_i}}}{{{\text{D}}t}} = {{\boldsymbol{u}}_i}\\ {p_i} = c_0^2\left( {{\rho _i} - {\rho _0}} \right) \end{array} \right\} $$ (1) 式中,
$ {{\boldsymbol{r}}_i} $ ,${{\boldsymbol{u}}_{\boldsymbol{i}}}$ ,$ \;{p_i} $ ,${\rho _i}$ 分别代表粒子$i$ 的位移、速度、压力、密度, 下标$j$ 代表粒子$i$ 紧支域内的相邻粒子.$ g $ ,${\rho _0}$ 分别代表重力加速度和流体参考密度,$ c_0^{} $ 为人工声速. 在弱可压条件下,$ c_0^{} $ 需满足${c_0} \geqslant 10\max \left( {{U_{\max }},\sqrt {{\raise0.7 ex\hbox{${{p_{\max }}}$} \mathord{\left/ {\vphantom {{{p_{\max }}} {{\rho _0}}}}\right.} \lower0.7 ex\hbox{${{\rho _0}}$}}} } \right)$ , 即根据流场中最大速度${U_{\max }}$ 或最大期望压力${p_{\max }}$ 取值.${\left\langle {\nabla p} \right\rangle _i}^C$ 为压力梯度, 它的离散方式具体见下节.对于δ-SPH模型, 连续性方程引入密度耗散项
$ {D_i}^\rho = \delta h{c_0}{D_i} $ , 动量方程中加入人工黏性项${{\textit{π}} _i}^v$ , 以提升SPH方法的计算稳定性[24]. 本文取δ = 0.1,${D_i}$ 写为$$ \left. \begin{array}{l} {D_i}^{} = 2\displaystyle\sum\limits_j {{\psi _{ji}}\dfrac{{\left( {{{\boldsymbol{r}}_j} - {{\boldsymbol{r}}_i}} \right) \cdot {\nabla _i}{W_{ij}}}}{{{{\left\| {{{\boldsymbol{r}}_j} - {{\boldsymbol{r}}_i}} \right\|}^2}{\text{ + }}0.01{h^2}}}{V_j}} \hfill \\ {\psi _{ji}} = {\rho _j} - {\rho _i} - \dfrac{1}{2}\left( {\left\langle {\nabla \rho } \right\rangle _i^{CSPM} + \left\langle {\nabla \rho } \right\rangle _j^{CSPM}} \right) \cdot(\boldsymbol{r}_j-\boldsymbol{r}_i)\hfill \\ \end{array} \right\} $$ (2) 式中,
$ \left\langle {\nabla \rho } \right\rangle _i^{CSPM} $ [24]的计算采用修正光滑粒子近似方法, 更多细节见下一节论述.人工黏性项
${{\textit{π}} _i}^v$ 写为$$ \left. \begin{array}{l} {{\textit{π}} _i}^v = \alpha h{c_0}\dfrac{{{\rho _0}}}{{{\rho _i}}}\displaystyle\mathop \sum \limits_j {\phi _{ij}}{\nabla _i}{W_{ij}}{V_j} \hfill \\ {\phi _{ij}} = 2(n + 2)\dfrac{{({{\boldsymbol{u}}_j} - {{\boldsymbol{u}}_i}) \cdot \left( {{{\boldsymbol{r}}_j} - {{\boldsymbol{r}}_i}} \right)}}{{{{\left\| {{{\boldsymbol{r}}_j} - {{\boldsymbol{r}}_i}} \right\|}^2}{\text{ + }}0.01{h^2}}} \hfill \\ \end{array} \right\} $$ (3) 式中, n =2, 光滑长度
$h = {\alpha _s}\Delta x$ ,${\alpha _s}$ 为光滑长度系数,$\Delta x$ 为粒子间距.$\alpha $ 为黏性系数, 在SPH模拟中, 通过在消波区内增大黏性系数实现数值消波, 具体设置为$$ \alpha = \left\{ {\begin{array}{*{20}{l}} {0.02},&{0 < x \leqslant {L_{{\text{damp}}}}} \\ {0.6(x - {L_{{\text{damp}}}})/(L - {L_{{\text{damp}}}})},&{{L_{{\text{damp}}}} < x \leqslant L} \end{array}} \right. $$ (4) 式中,
$ L $ 为数值水池总长度,$ L - {L_{{\text{damp}}}} $ 为消波区长度. 在本文模拟中, 采用两种常用的光滑长度系数$ {\alpha }_{s}=1.3, 2.0 $ 进行计算. 其中${\alpha _s} = 1.3$ 时, 由于核函数紧支域内粒子数量更少, 因此计算量小, 但是离散误差相比${\alpha _s} = 2.0$ 时更大, 数值耗散也会更大. 本文将开发改进SPH算法, 力争在${\alpha _s} = 1.3$ 时, 也获得较好的计算精度和较小的数值耗散.1.2 低能量耗散的修正δ-SPHC算法
在传统SPH方法中, 梯度算子粒子近似[25]可写为
$$ \left\langle {\nabla f} \right\rangle _{i}=\sum _{j}\left({f}_{j} + {f}_{i}\right){\nabla }_{i}{W}_{ij}{V}_{j} $$ (5) 因此, 在传统δ-SPH, 压力梯度算子离散
${\left\langle {\nabla p} \right\rangle _i}$ [25]为$$ \left\langle \nabla p \right\rangle _{i}=\sum _{j}\left({p}_{j} + {p}_{i}\right){\nabla }_{i}{W}_{ij}{V}_{j} $$ (6) 然而, 在粒子分布不均匀或核函数截断情况下, 该离散形式不能较好地确保计算精度. 对此, 许多修正方法被提出, 其中常用的有修正光滑粒子法(corrective smoothed particle method, CSPM)[26-27]. 梯度算子采用CSPM修正后为
$$ \left. \begin{array}{l} \left\langle \nabla f\right\rangle{ _i}^{CSPM} = \displaystyle\mathop \sum \limits_j \;({f_j} - {f_i})\;{\nabla _i}W_{ij}^{CSPM}{V_j} \hfill \\ {\nabla _i}W_{ij}^{CSPM}: = {{\boldsymbol{M}}_i}{\nabla _i}{W_{ij}} \hfill \\ {{\boldsymbol{M}}_i}: = {\left[ {\displaystyle\mathop {\sum \;}\limits_j ({{\boldsymbol{r}}_j} - {{\boldsymbol{r}}_i}) \otimes {\nabla _i}{W_{ij}}{V_j}} \right]^{ - 1}} \end{array} \right\} $$ (7) δ-SPH算法中的
$ \left\langle {\nabla \rho } \right\rangle _{}^{CSPM} $ 项采用上式进行离散. 然而, 使用$\left\langle {\nabla f } \right\rangle _{i}^{CSPM}$ 对压力梯度进行离散时, 矩阵$ {{\boldsymbol{M}}_i} $ 的存在使得粒子对($i \leftrightarrow j$ )之间的相互作用不能保持对称性, 从而不能确保动量精确守恒性, 进而影响能量的守恒性. 因此, 在δ-SPH模型中, 动量方程的压力梯度项并不使用该方式进行离散, 只对连续方程中密度梯度项进行修正.经Wen等[22]和Zago等[23]研究启发, 采用上述CSPM算法中修正的核函数梯度
$ {\nabla _i}W_{ij}^{CSPM} $ , 有助于改善SPH计算中能量的衰减问题, 但是为了保证动量守恒, 也就是确保粒子对($i \leftrightarrow j$ )之间相互作用力的对称性, 需要对CSPM格式进行改进. 对此, Zago等[23]提出了一种改进后的修正格式, 但是需要显式搜索自由液面, 这带来了额外的计算负担, 而本文基于CSPM修正形式, 参考Sun等[28]处理梯度算子的方法, 采用另一种具有对称性的梯度算子离散形式, 推导过程见下:首先, 根据核函数导数的特性[25]有
$$ {\nabla _i}{W_{ij}} = - {\nabla _j}{W_{ij}} $$ (8) 将上式代入式(5), 有
$$ \left\langle \nabla f\right\rangle { _i} = \mathop \sum \limits_j \;\left( {{f_i}{\nabla _i}{W_{ij}} - {f_j}{\nabla _j}{W_{ij}}} \right){V_j} $$ (9) 将上式
$ {\nabla _i}{W_{ij}} $ 和$ {\nabla _j}{W_{ij}} $ 采用CSPM核函数修正, 得$$ \left\langle \nabla f\right\rangle{ _i}^C = \mathop \sum \limits_j \;\left( {{f_i}{\nabla _i}W_{ij}^{CSPM} - {f_j}{\nabla _j}W_{ij}^{CSPM}} \right){V_j} $$ (10) 上式可进一步等价写为
$$ \left\langle\nabla f\right\rangle{ _i}^C = \mathop \sum \limits_j ({f_i}{{\boldsymbol{M}}_i} + {f_j}{{\boldsymbol{M}}_j}){\nabla _i}W_{ij}^{}{V_j} $$ (11) 由上式可见,
$\left\langle \nabla f\right\rangle{ _i}^C$ 的离散形式具有对称性. 事实上, 式(11)由Oger等[29]最早得出, Ganzenmüller等[30]证明其与近场动力学的算子相一致, 但是对于该离散格式在数值水池抗能量非物理性耗散方面的应用, 现有文献中鲜有论述.在本文改进的SPH算法中, 动量方程中的压力梯度项
$ \left\langle \nabla p\right\rangle{ _i}^C $ 采用公式(11)进行离散, 可写为$$ \left\langle \nabla p\right\rangle{ _i}^C = \mathop \sum \limits_j \;\left( {{p_j}{{\boldsymbol{M}}_j} + {p_i}{{\boldsymbol{M}}_i}} \right)\;{\nabla _i}W_{ij}^{}{V_j} $$ (12) 最后控制方程写为
$$ \left.\begin{array}{l} \dfrac{{{\text{D}}{\rho _i}}}{{{\text{D}}t}} = - {\rho _i}\displaystyle\mathop \sum \limits_j \;\left( {{{\boldsymbol{u}}_j} - {{\boldsymbol{u}}_i}} \right)\;{\nabla _i}{W_{{{ij}}}}{V_j} + {D_i}^\rho \hfill \\ \dfrac{{{\text{D}}{{\boldsymbol{u}}_i}}}{{{\text{D}}t}} = - \dfrac{1}{{{\rho _i}}}\displaystyle\mathop \sum \limits_j \;\left( {{p_j}{{\boldsymbol{M}}_j} + {p_i}{{\boldsymbol{M}}_i}} \right)\;{\nabla _i}W_{ij}^{}{V_j} + {\boldsymbol{g}} + {{\textit{π}} _i}^v \hfill \\ \dfrac{{{\text{D}}{{\boldsymbol{r}}_i}}}{{{\text{D}}t}} = {{\boldsymbol{u}}_i}\hfill \\ {p_i} = c_0^2\left( {{\rho _i} - {\rho _0}} \right) \end{array} \right\} $$ (13) 本文修正算法记为δ-SPHC. 对比Zago等[23]的算法, 本修正算法避免了自由液面搜索的步骤; 另外, 相对于δ-SPH, δ-SPHC修正算法并不引入新的核函数修正矩阵, 只是再次调用密度耗散项中的矩阵
$ {{\boldsymbol{M}}_i} $ 对动量方程中压力梯度算子进行修正, 因此不会造成计算成本和时间的显著增加. 下文将采用振荡液滴算例进行修正算法抗能量非物理性衰减效果的测试和讨论, 并应用于数值波浪水池规则波与不规则波的数值模拟和试验验证. 本文的数值积分采用四阶龙格库塔法, 具体可参考文献[31].2. 振荡液滴算例验证
2.1 收敛性分析与精度验证
振荡液滴算例无需施加固壁边界, 初始条件简单, 被广泛应用于SPH等数值新算法的验证[32]. Antuono等[33]采用振荡液滴算例验证了δ-SPH模型中的能量守恒问题. 基于此, 本文采用该算例对修正算法δ-SPHC减少能量非物理性衰减的效果进行验证.
该算例中, 初始速度
${\boldsymbol{V}}(u,v)$ 、初始半径为$R$ 的圆形液滴在向心体积力场${\boldsymbol{F}}$ 中振荡, 其中${\boldsymbol{F}}(x,y) = - {B^2}{\boldsymbol{r}}(x,y)$ ,$B = 1$ ,$R = 1$ . 图1展示了初始时刻${t_0}$ 时液滴的速度场分布[33], 可表示为$$ \left. \begin{array}{l} u({t_0}) = A({t_0})x \hfill \\ v({t_0}) = - A({t_0})y \end{array} \right\} $$ (14) $ A({t_0}) $ 为控制液滴初始速度场的参数. 在本文模拟中,$ A({t_0}){\text{ = }}1 $ . 初始压力场分布如图2所示, 压力场[33]表示为$$ p\left( {x,y,{t_0}} \right) = B\left[ {1 - \left( {{x^2} + {y^2}} \right)} \right] $$ (15) 首先, 为检验SPH修正算法的收敛性, 分别在粒子分辨率为
$ R/\Delta x = 50 $ ,$R/\Delta x = 100$ ,$R/\Delta x = 200$ 时, 开展振荡液滴运动特性研究. 由于振荡液滴质心与向心体积力${\boldsymbol{F}}$ 的汇聚点重合, 液滴不受其他外力作用, 液滴在原位保持振荡. 因此, 理论上振荡液滴的角动量、线动量都守恒, 即${M_A}(t) = \left| {\displaystyle\sum\limits_i {{{\boldsymbol{r}}_i}(t) \times {m_i}{{\boldsymbol{u}}_i}} (t) } \right|= 0$ ,${M_L}(t) = \left|\displaystyle\sum\limits_i {{m_i}{{\boldsymbol{u}}_i}} \right| = 0$ .图3展示在光滑长度系数
${\alpha _s}$ =1.3和${\alpha _s}$ =2.0条件下, 不同粒子分辨率时, 角动量${M_A}$ 随时间的变化情况. 可见, 在三种粒子分辨率下,${M_A}$ 的值围绕理论值${M_A} = 0$ 波动的量级均在${10^{ - 5}}$ 以内, 且随着粒子分辨率的增加, 角动量${M_A}$ 的误差范围逐渐变小. 在粒子分辨率$R/\Delta x = 200$ 时,${M_A}$ 趋近于零. 即在$R/\Delta x = 200$ 时, 角动量误差逐渐收敛于0, 说明该修正算法在提高粒子分辨率时能实现良好的角动量守恒性.![]()
图 3 不同粒子分辨率下角动量${M_A}$ 的时历曲线Figure 3. Time evolution of the angular momentum${M_A}$ with different particle resolutions图4展示了光滑长度系数
${\alpha _s} = 1.3$ 和${\alpha _s} = 2.0$ 条件下,$R/\Delta x = 200$ 分辨率时, 线动量${M_L}$ 随时间的变化情况. 观察可得, 线动量误差在−10−14 ~ 10−14范围内波动, 可以认为本修正算法的线动量实现了精确守恒.同时, 定义液滴振荡过程的机械能衰减率为
$ {E_{per}}(t) = ({{E(t) - {E_{{int} }})} \mathord{\left/ {\vphantom {{E(t) - {E_{{int} }})} {{E_{{int} }}}}} \right. } {{E_{{int} }}}} $ , 其中$ {E_{{int} }} $ 为初始机械能,$ E(t) $ 为$t$ 时刻的机械能. 液滴振荡过程的机械能衰减率随时间变化结果如图5所示. 可以发现, 随着粒子分辨率的提高, 机械能衰减率${E_{per}}$ 越来越小. 在$R/\Delta x = 100$ ,$R/\Delta x = 200$ 时, 两者机械能衰减率基本接近, 计算时长达$t = 120$ 时, 机械能衰减控制在5%以内.![]()
图 5 αs = 2.0, 不同粒子分辨率时机械能衰减率${E_{per}}$ 时历曲线Figure 5. Time evolution of the decay rate of mechanical energy at different particle resolutions with αs = 2.0图6展示了
$R/\Delta x = 200$ 时, 振荡液滴在不同时刻的理论形态与修正算法δ-SPHC模拟的形态. 可以发现, 在两种光滑长度条件下, δ-SPHC模拟的振荡形态均与理论值吻合良好, 且在t = 8.75T周期后, 液滴仍保持较为光滑的形态. 由此可见, 本修正算法具有较高的精度.综上可见, 修正SPH算法在粒子分辨率
$R/\Delta x = 200$ 时, 具有较好的数值精度和动量与能量守恒性. 因此, 下文将在粒子分辨率$R/\Delta x = 200$ 时, 依托该算例进一步讨论不同算法抗能量衰减效果的差异.![]()
图 6 αs = 1.3(上)和αs = 2.0(下)时, 振荡液滴形态: 理论值(虚线)和SPH结果(压力云图)Figure 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)2.2 传统δ-SPH与修正δ-SPHC算法的对比研究
基于上节收敛性分析与精度验证, 本节在粒子分辨率
$R/\Delta x = 200$ 条件下, 开展传统δ-SPH模型与改进δ-SPHC模型抗能量衰减效果的对比研究. 图7展示了δ-SPH与δ-SPHC模拟的振荡液滴动能随时间的变化曲线. 可见, 在传统算法δ-SPH结果中, αs = 1.3时, 动能随着时间的衰减十分明显. 在αs = 2.0时, 衰减有所减缓. 但是, 取$h = 2\Delta x$ 相比$h = 1.3\Delta x$ , 由于相邻粒子数量的增多, 会导致额外的计算成本和时间, 且从图7的局部放大图可以看出, 在较长时间的计算下, 仍会出现轻微动能衰减. 对比修正算法δ-SPHC, 在光滑长度系数αs =1.3时, 动能幅值基本保持不变, 已体现出较好的抗能量衰减效果; 在光滑长度系数αs = 2.0时, 本修正算法实现了良好的能量守恒. 此外, 由局部放大图可以看出, 在修正算法αs = 2.0时, 动能幅值与初始幅值始终保持一致, 取得了较为精确的计算结果. 而在δ-SPHC结果中,${\alpha _s} = 1.3$ 时, 动能第一个峰值较初始值高, 主要原因在第一个周期的计算中, 粒子会从格子分布状态(图8(a))突然过渡到一种各项同性的均匀分布状态(图8(b)). 在粒子分布变化过程中, 存在一些粒子的瞬时分布十分不均匀, 此时, 这不均匀性会造成$\sum\nolimits_j { {\nabla _i}{W_{ij}}{V_j} \ne 0} \;$ . 也就是说, 部分粒子的运动并不是由于压力梯度造成, 而是粒子分布的不均匀性导致的数值加速度产生, 这对粒子的总动能造成了一定的非物理性影响. 但是, 从图7下方的放大图可以看出, 这种由于粒子分布突变引起的动能误差控制在1%以内. 类似现象在Colagrossi等[34]的研究结果也有过报道, 可通过改善初始粒子分布得以解决.![]()
图 7 δ-SPH与δ-SPHC动能比较: 总体图(上)与局部放大图(下), 其中黑矩形框为放大区域Figure 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为进一步突显修正算法的低能量耗散特性, 图9呈现了两种算法中, 机械能衰减率的时历曲线. δ-SPHC算法在两种光滑长度系数下的计算结果中, 机械能衰减率均较小, 基本接近零, 说明本文修正算法能够有效减少机械能的非物理性耗散. 其中, 在光滑长度系数
${\alpha _s} = 2.0$ 时效果更为显著. 值得一提的是, δ-SPHC算法在${\alpha _s} = 2.0$ 时的机械能衰减率相比δ-SPH算法在${\alpha _s} = 1.3$ 时的机械能衰减率更小, 表明本修正SPH算法的一大优势是能在较小光滑长度系数下实现较为准确的模拟, 从而显著减少计算时间.![]()
图 9 不同算法和光滑长度系数下, 振荡液滴机械能衰减率时历曲线Figure 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{\text{/}}\Delta x$ =200,${\alpha _s}$ = 2.0参数下, δ-SPH与δ-SPHC模型模拟振荡液滴的计算时长. 可见, 两种算法单步计算时间相当. 具体而言, 虽然修正算法δ-SPHC对压力梯度项进行了核函数修正, 但相对δ-SPH算法, 计算时间仅增加1.09倍, 计算量增加并不显著. 在下节中, 该算法将被应用于数值波浪水池, 验证其改善SPH模拟远距离波浪传播时的波高抗衰减效果.表 1 δ-SPH算法与δ-SPHC算法计算效率比较Table 1. Comparison of computational efficiency between δ-SPH and δ-SPHC schemesScheme CPU time
of single step/st2/t1 δ-SPH t1 = 0.5963 1.09 δ-SPHC t2 = 0.6576 1.09 3. SPH数值造波与试验验证
3.1 试验与数值水池设置
3.1.1 试验水池
作者在中山大学海洋工程与技术学院波浪水槽中开展造波试验. 如图10所示, 该水池长度L=16 m, 水深d = 0.266 m, 采用推板造波方式生成规则波和不规则波, 推板运动曲线如图11所示. 试验时, 波高仪设置在距离造波板平衡位置
${x_{{\rm{probe}}}}$ =2.37 m, 4.37 m, 5.37 m, 6.37 m, 15 m处.![]()
图 11 造波推板横向位移随时间变化: (a)周期为T = 1 s的规则波, (b)不规则波Figure 11. Paddle motions of the regular wave with (a) period T = 1 s and (b) irregular wave3.1.2 SPH数值水池设置
在SPH模拟中, 固定虚粒子边界法[35, 36]能实现固壁边界的精确处理, 有效防止壁面穿透, 且保证壁面处压力连续, 因而得到了广泛应用.
本文SPH数值模拟中, 数值水池长度L=60 m, 边界条件采用固定虚粒子进行施加. 如图12所示, 固定虚粒子厚度不小于核函数半径, 其压力从流体粒子外插得到, 更多细节可参考文献[37]. 本文SPH模拟中, 通过强制流体粒子与壁面虚粒子之间黏性力为0, 实现壁面边界处的自由滑移边界条件. 数值水池左侧的运动造波板也由虚粒子组成, 通过控制造波板虚粒子的总体位置发生改变, 实现推板造波. 具体实现方式如下.
![]()
图 12 SPH数值水池示意图, 红竖线表示波高仪, 选取距离波高仪一个粒子间距内自由面粒子的最大高度为波面高度Figure 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模拟中的造波板位移
$ {x^{SPH}} $ 、速度$ {v^{SPH}} $ 及加速度$ {a^{SPH}} $ , 即$$ \left.\begin{array}{l} {x^{SPH}} = {x^k} + \dfrac{{({x^{k + 1}} - {x^k})(t - {t^k})}}{{{t^{k + 1}} - {t^k}}} \hfill \\ {v^{SPH}} = {v^k} + \dfrac{{({v^{k + 1}} - {v^k})(t - {t^k})}}{{{t^{k + 1}} - {t^k}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \hfill \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {a^{SPH}} = {a^k} + \dfrac{{({a^{k + 1}} - {a^k})(t - {t^k})}}{{{t^{k + 1}} - {t^k}}} \end{array}\right\} $$ (16) 其中, 上标
$SPH$ 表示数值水池模拟中的物理量, 上标$k$ 表示试验中造波板运动输出的时间步. 试验中造波板的速度${v^k}$ 与加速度${a^k}$ 均采用线性差分计算获得, 即$$ \left.\begin{array}{l} {v^k} = \dfrac{{{x^{k + 1}} - {x^k}}}{{{t^{k + 1}} - {t^k}}} \hfill \\ {a^k} = \dfrac{{{v^{k + 1}} - {v^k}}}{{{t^{k + 1}} - {t^k}}} \end{array} \right\}$$ (17) 此外, 值得注意的是, 当波浪传播到数值水池另一端时, 往往会产生波浪反射, 降低造波精度, 这就需要在数值水池右端构建数值消波区. 类似于网格算法中增大数值黏性来实现阻尼消波, 本文通过提高SPH数值水池消波区的黏性系数, 使波浪进入消波区时能够在短时间内衰减. 具体设置为 :在数值水池设置消波区
$ {L_{{\text{damp}}}} $ , 黏性系数$\alpha $ 线性增大, 详见公式(4). 本文计算中,$ L_{\text {damp }}=20 $ m.3.2 δ-SPH算法在高粒子分辨率下的精度验证
Colagrossi等[16]指出δ-SPH模型能较为准确地重现在物理黏性条件下的波浪传播和衰减过程. 因此, 本文基于δ-SPH模型开发低能量耗散的数值水池. 首先, 在
${\alpha _s}$ = 2.0条件下, 验证δ-SPH计算结果的精度和收敛性. 选取$d/\Delta x$ =15, 30, 60, 即水深方向上粒子数分别为15, 30, 60进行规则波模拟, 在三种粒子分辨率条件下监测${x_{\rm{probe}}}$ =2.37处的波面高度随时间变化, 如图13所示. 可以看出, 数值波高结果变化趋势均与试验波高结果吻合良好, 在$d/\Delta x$ = 60时, 数值波高收敛于试验波高. 基于以上讨论, 在$d/\Delta x$ = 60条件下, δ-SPH模型能较为精确实现造波及波浪传播的模拟. 但是值得一提的是, δ-SPH模型在${\alpha _s}$ < 2.0时, 容易出现波高的非物理性衰减, 具体论述见下一节.![]()
图 13${\alpha _s}$ =2.0时, 不同粒子分辨率下波面高度的δ-SPH模拟值与试验值对比Figure 13. Comparison of the wave elevations between δ-SPH results and experimental data, SPH results with different particle resolutions and${\alpha _s} = 2.0$ are provided3.3 δ-SPH算法中的波浪衰减现象
本节采用δ-SPH模型, 在光滑长度系数
${\alpha _s}$ =1.3和${\alpha _s}$ = 2.0条件下, 进行规则波造波模拟. 为了测定波高变化曲线, 沿水池长度方向, 距离造波板初始位置${x_{\rm{probe}}}$ =2.37 m, 4.37 m, 5.37 m, 6.37 m, 15 m处设置5个数值波高仪.在
${\alpha _s}$ =1.3时, 图14展示了计算时间为t = 27 s的波面形态, 可以较为明显地看出波浪在x = 10 m处开始出现衰减. 选取${x_{\rm{probe}}}$ = 2.37 m, 5.37 m, 15 m处的波面高度时历曲线进行比较, 如图15所示. 在${x_{\rm{probe}}}$ = 2.37 m处, 波高与试验相差较大,且随着传播距离的增大, 波高逐渐衰减, 在${x_{\rm{probe}}}$ = 15 m处, 波高衰减严重. 可以看出, 在传统δ-SPH中, 光滑长度系数${\alpha _s}$ = 1.3条件下无法准确模拟波浪传播. 图16展示了δ-SPH算法在光滑长度系数${\alpha _s}$ = 2.0时计算的波面形态. 与振荡液滴类似, 相比光滑长度系数${\alpha _s}$ = 1.3时, 衰减速度有所减缓, 与Colagrossi等[16]的结论相符. 这表明, 增大${\alpha _s}$ 可以有效减缓SPH的数值衰减, 但在大尺度、长时间、远距离的海浪模拟下, 光滑长度系数${\alpha _s}$ 取2.0时, 将带来巨大的计算量, 不利于工程应用. 且比较不同距离处波高, 如图17所示, 可以发现在${x_{\rm{probe}}}$ = 15 m时, 仍存在明显波高衰减. 因此, 改善SPH方法在数值波浪模拟时能量的非物理性耗散, 在较低光滑长度系数下实现数值水池造波的高精度高效率模拟具有重要意义.![]()
图 14 在${\alpha _s}$ = 1.3,$d/\Delta x$ = 60条件下, 传统δ-SPH模拟的波面形态Figure 14. Wave surface simulated by δ-SPH method with${\alpha _s}$ = 1.3 and$d/\Delta x$ = 60![]()
图 15 在${\alpha _s} = 1.3$ ,$d/\Delta x$ = 60条件下, 传统δ-SPH模型计算的不同位置波面高度时历曲线Figure 15. Wave elevation probed at different positions simulated by δ-SPH method with${\alpha _s} = 1.3$ and$d/\Delta x$ = 60![]()
图 16 在${\alpha _s} = 2.0$ ,$d/\Delta x$ = 60条件下, 传统δ-SPH模拟的波面形态Figure 16. Wave surface simulated by δ-SPH method with${\alpha _s}$ = 2.0 and$ d/\Delta x $ = 60![]()
图 17 在${\alpha _s}$ = 2.0,$ d/\Delta x $ = 60条件下, 传统δ-SPH模拟的不同位置波面高度时历曲线Figure 17. Wave elevation probed at different positions simulated by δ-SPH method with${\alpha _s}$ = 2.0 and$d/\Delta x$ = 604. 修正算法δ-SPHC的波浪传播模拟
4.1 规则波模拟
4.1.1 收敛性分析与精度验证
首先, 基于规则波传播算例, 校核本文修正的δ-SPHC算法精度与收敛性. 如图18所示, 在
${\alpha _s} = 1.3$ 时, 选取${x_{\rm{probe}}}$ = 2.37 m处数值波高结果与试验结果进行比较, 可以看出, 水深方向不同数量粒子的计算结果略有差别, 且随着$d/\Delta x$ 的增大, 数值波高趋近于试验波高. 在$d/\Delta x = 60$ 时, SPH结果与试验结果吻合良好. 与振荡液滴算例类似, 图18方框中的峰值可能是由于粒子由初始的格子状分布到各向同性均匀分布时引起的能量瞬时误差. 试验波高曲线与SPH数值结果曲线在初始阶段第一个峰值存在相位差, 可能是由于试验波高采集与造波板运动存在一定的时间差. 在${\alpha _s}{\text{ = }}2.0$ 时, 水深方向不同粒子数量的计算结果与试验结果对比如图19所示. 同样可以看到, 在$d/\Delta x$ = 60时, 数值波高曲线与试验结果基本吻合. 另外, 由于物理试验水池长度仅有16 m, 消波区长度设置有限, 易导致试验后期消波不彻底, 造成波浪周期的变化, 即随着时间增加, 周期略小于1 s; 而数值水池长度$L$ = 60 m, 能较为准确地重现规则波传播, 并且确保周期保持不变. 因此, 可以看到, 在图18、图19中, t = 20 s之后SPH结果的波峰和波谷时刻与试验值存在略微误差, 但是波高和波谷的幅值与试验值吻合良好. 基于以上讨论可以得出, 在水深方向粒子数$d/\Delta x$ = 60条件下, 计算结果能够达到预期精度要求. 因此, 下文将在$d/\Delta x$ = 60粒子设置下, 进行修正δ-SPHC算法的抗能量非物理性衰减效果的讨论.![]()
图 18 在αs = 1.3和$d/\Delta x$ = 15, 30, 60条件下, δ-SPHC计算的波面高度时历曲线与试验值对比Figure 18. Comparison between δ-SPHC results and experimental data for wave elevation: SPH results with αs = 1.3 and$d/\Delta x$ = 15, 30, 60![]()
图 19 在αs = 2.0和$d/\Delta x$ =15, 30, 60条件下, δ-SPHC计算的波面高度时历曲线与试验值对比Figure 19. Comparison between δ-SPHC results and experimental data for wave elevation: SPH results with αs = 2.0 and$d/\Delta x$ = 15, 30, 60图20呈现了在
${x_{\rm{probe}}}$ = 2.37 m处, 规则波波高数值模拟结果与试验结果的对比. 采用传统δ-SPH计算时, 光滑长度系数${\alpha _s}$ =1.3条件下, 由于核函数紧支域内粒子数量少, 数值结果与试验结果相差较大, 而光滑长度系数${\alpha _s}$ = 2.0时, 数值波高结果接近试验结果. 但值得注意的是, 采用δ-SPHC, 光滑长度系数${\alpha _s}$ = 1.3时, 数值结果与试验结果能较好吻合, 基本达到传统δ-SPH算法在光滑长度系数${\alpha _s}$ = 2.0时的计算精度, 与上文振荡液滴算例得出的结论一致, 再次说明了本文修正δ-SPHC算法的优势, 即在光滑长度系数$ {\alpha }_{s}=1.3, 2.0 $ 时, 计算结果与试验结果均能较好吻合. 图21展示了数值水池的压力场, 可见流场压力分布均匀, 在运动较为剧烈的推板附近, 也能确保压力场稳定. 基于以上讨论, 可以认为修正的δ-SPHC算法具有较高的精度.![]()
图 20${x_{{\rm{probe}}}}$ =2.37 m处, 不同SPH算法下的波面高度数值结果与试验值对比Figure 20. Comparison between different SPH simulations and experimental data for wave elevations at${x_{{\rm{probe}}}}$ = 2.37 m![]()
图 21 规则波压力云图(δ-SPHC,${\alpha _s}$ = 2.0,$d/\Delta x$ = 60)Figure 21. Pressure contour of regular wave (δ-SPHC,${\alpha _s}$ = 2.0,$d/\Delta x$ = 60)4.1.2 δ-SPHC的抗能量非物理性衰减结果分析
经过前期收敛性与精度验证, 本节在水深方向粒子数取
$d/\Delta x$ = 60时, 开展修正算法δ-SPHC在两种光滑长度系数${\alpha _s}$ = 1.3, 2.0下的抗能量非物理性衰减效果讨论. 在${\alpha _s}$ = 1.3时, 如图22所示, 观察在修正算法模拟下的自由面形态, 可以看出在波浪传播20 m时, 波面仍与推板附近波面基本持平. 监测不同距离下的波面高度如图23所示, 与试验结果对比可以看出, 在${x_{\rm{probe}}}$ = 15 m处, δ-SPHC的数值波高与推板附近波高基本接近, 说明能量衰减现象得到极大改善. 图24展示了在${\alpha _s}$ = 2.0时, 规则波传播的波面形态. 可见, 在远距离传播时, 波面仍保持在同一高度, 说明了修正算法抗能量衰减的有效性. 对比图23和图25, 可见在${x_{\rm{probe}}}$ = 15 m位置处,${\alpha _s}$ = 2.0比${\alpha _s}$ = 1.3计算的波高更加精确和稳定, 说明光滑长度系数的增大有利于减少离散误差, 其中在${x_{\rm{probe}}}$ = 15 m处出现略微振荡, 原因主要是由于数值误差在远距离、长时间上的累积造成, 但是相比传统算法, 修正算法的精度得到了有效改善.![]()
图 22${\alpha _s}$ =1.3,$d/\Delta x$ = 60时, δ-SPHC算法模拟的波面形态Figure 22. The wave surface simulated by δ-SPHC with${\alpha _s}$ = 1.3 and$d/\Delta x$ = 60![]()
图 23${\alpha _s}$ =1.3,$d/\Delta x$ =60时, δ-SPHC模拟的不同位置波面高度时历曲线Figure 23. Wave elevations probed at different positions simulated by δ-SPHC method with${\alpha _s}$ = 1.3 and$d/\Delta x$ = 60![]()
图 24${\alpha _s}$ =2.0,$d/\Delta x$ =60时, δ-SPHC算法模拟的波面形态Figure 24. The wave surface simulated by δ-SPHC with${\alpha _s}$ = 2.0 and$d/\Delta x$ =60![]()
图 25${\alpha _s}$ = 2.0,$d/\Delta x$ = 60时, δ-SPHC模拟的不同位置波面高度时历曲线Figure 25. Wave elevations probed at different positions simulated by δ-SPHC method with${\alpha _s}$ = 2.0 and$d/\Delta x$ = 60如图26所示, 综合对比δ-SPH, δ-SPHC与试验方法在
${x_{\rm{probe}}}$ = 6.37 m的波高曲线, 可见相比于传统δ-SPH模型, 修正的δ-SPHC算法在${\alpha _s} = 1.3$ 时的波高结果更加准确. 值得一提的是, 其计算结果与传统δ-SPH模型在${\alpha _s}$ =2.0时波高结果基本重合, 并与试验波高吻合良好. 此外, 采用δ-SPHC修正算法, 在光滑长度系数${\alpha _s}$ =2.0时, 波高计算结果与试验波高高度吻合.![]()
图 26 在${x_{{\rm{probe}}}}$ = 6.37 m处, 规则波波面高度数值模拟结果与试验结果Figure 26. Time evolution of the regular wave elevation at${x_{{\rm{probe}}}}$ = 6.37 m: SPH results and experimental data为进一步定量验证δ-SPHC的模拟精度, 表2给出了试验和SPH模拟中,
${x_{\rm{probe}}}$ = 6.37 m测点处的平均波高值及相对误差, 后者定义为$\varepsilon = {{\left| {{H_{{\text{SPH}}}} - {H_{{\text{EXP}}}}} \right|} \mathord{\left/ {\vphantom {{\left| {{H_{{\text{SPH}}}} - {H_{{\text{EXP}}}}} \right|} {{H_{{\text{EXP}}}}}}} \right. } {{H_{{\text{EXP}}}}}}$ ,${H_{{\text{EXP}}}}$ 为试验测量的平均波高. 可见, 在两种光滑长度系数下, δ-SPHC算法在较远处的波高监测结果与试验相对误差均较小. 本节讨论表明, 修正算法能有效解决能量非物理性衰减问题, 且在较低光滑长度系数时能得到更为准确的结果.表 2 不同算法在${{\boldsymbol{x}}_{{\bf{probe}}}}$ = 6.37 m 处平均波高对比Table 2. Comparison of average wave heights probed at 6.37 m simulated with different methodsMethod Exp. δ-SPH, αs = 1.3 δ-SPH, αs = 2.0 δ-SPHC, αs = 1.3 δ-SPHC, αs = 2.0 wave height H/m 0.0428 0.0224 0.0397 0.0419 0.0427 relative error $\varepsilon $/% − 47.66 7.24 2.10 0.23 4.1.3 不同人工黏性系数下δ-SPHC模拟结果分析
在数值水池中, 能量衰减来源主要有两部分: 一是由于流体物理黏性引起的能量耗散; 二是本文所聚焦的由于数值误差引起的能量非物理性耗散. 本文虽然未在SPH控制方程中添加物理黏性项, 但是使用了人工黏性项
${{\textit{π}} _i}^v$ , 这等效于增加了流体的黏性. 为进一步讨论$\delta $ -SPHC 在不同人工黏性系数$\alpha $ 下, 对长时间、远距离波浪传播模拟的能量耗散情况, 选取四个不同黏性系数, 即$\alpha $ = 0.01, 0.02, 0.05, 0.1开展数值模拟. 图27给出了在${x_{\rm{probe}}}$ = 6.37 m处监测的波高曲线. 可见, 在$\alpha $ = 0.01和$\alpha $ = 0.02时, 波浪的衰减较小, 与实验值吻合较好. 在$\alpha $ = 0.05时, 由于采用了更大的人工黏性, 等效增加了物理黏性, 因此波高衰减有所增加. 在$\alpha $ = 0.1时, 由于人工黏性引起的波高衰减更为明显, 与实验值相比误差较大. 综合以上分析, 在$\alpha $ = 0.02时, 能够获得较为稳定和精确的波高计算结果, 且与试验值吻合良好. 因此, 本文其他算例均基于$\alpha $ = 0.02进行模拟. 基于此, 在下文中进一步开展不规则波生成和传播的SPH模拟和验证.![]()
图 27 δ-SPHC在不同黏性系数α时,${x_{{\rm{probe}}}}$ = 6.37 m处波面高度时历曲线Figure 27. Time history of wave elevation at${x_{{\rm{probe}}}}$ = 6.37 m simulated by δ-SPHC with different α4.2 基于δ-SPHC的不规则波模拟
4.2.1 收敛性分析与精度验证
经过上文分析, 在光滑长度系数
${\alpha _s}$ = 1.3时, 修正算法δ-SPHC能够准确模拟规则波传播, 而本节算法, 在${\alpha _s}$ = 1.3时, 模拟和验证不规则波传播. 选取$d/\Delta x$ = 15, 30, 60进行收敛性分析, 如图28所示, 在水深方向粒子数$d/\Delta x$ =15, 30, 60时, 在${x_{\rm{probe}}}$ =2.37处探测的波高与试验值基本一致. 此外, 在不规则波模拟中, δ-SPHC实现了稳定的压力场预报(见图29), 进一步说明了δ-SPHC方法所构建的数值水池的精度. 由于$d/\Delta x$ = 60时, 数值结果已基本收敛于试验值, 因此下文选取$d/\Delta x$ = 60进行不规则波模拟结果的讨论.![]()
图 28${\alpha _s}$ =1.3和不同粒子分辨率条件下, 不规则波δ-SPHC模拟结果的收敛性分析Figure 28. Convergence analysis of irregular waves simulated by δ-SPHC with αs=1.3 at different particle resolutions![]()
图 29 不规则波压力云图(δ-SPHC,${\alpha _s}$ =2.0,$d/\Delta x$ =60)Figure 29. Pressure contour of irregular wave (δ-SPHC,${\alpha _s}$ =2.0,$d/\Delta x$ =60)4.2.2 δ-SPHC抗能量非物理性衰减结果分析
如图30所示, 监测
${x_{\rm{probe}}}$ = 6.37处波高时历曲线, 可见在传统δ-SPH计算结果中, 与试验结果相比, 波高衰减明显; 而修正算法结果与试验结果基本吻合, 说明不规则波模拟中波高衰减现象也得到较好改善. 在δ-SPHC计算框架下, 数值水池生成不规则波的能量非物理性耗散问题得到解决, 并能够在光滑长度系数${\alpha _s}$ = 1.3时, 实现较高的造波精度.![]()
图 30$ x_{\text {probe }} $ = 6.37处, 不规则波的波面高度时历曲线Figure 30. Time evolution of the irregular wave elevation measured at${x_{{\rm{probe}}}}$ = 6.37 m5. 结论与展望
本文给出一种修正的δ-SPHC算法, 采用对称化的核函数导数修正格式, 对压力梯度进行离散, 确保粒子对之间作用力的对称性, 实现了精确线动量守恒, 且无需自由面粒子搜索等额外的复杂数值处理过程. 通过模拟振荡液滴算例, 说明其抗能量衰减的有效性. 随后将其应用于数值波浪水池. 研究结果表明:
(1) δ-SPH模型存在能量非物理性耗散问题, 增大光滑长度系数可以有效减缓能量衰减;
(2)修正的δ-SPHC算法具有良好的动量守恒性, 克服了传统核函数梯度修正的动量不守恒问题;
(3)修正δ-SPHC算法具有较好的低能量耗散特性. 在修正算法作用下, 机械能衰减得到有效控制;
(4)在修正δ-SPHC算法作用下, 数值波浪水池波高衰减问题大为改善. 在低光滑长度系数
${\alpha _s}$ =1.3时, 计算结果与试验结果吻合良好. 因此, 本修正算法可以有效节约计算成本, 减少计算时间. 未来在模拟大尺度、长时间、远距离波浪传播和波浪与结构物相互作用时, 有望发挥重要作用. -
![]()
图 3 不同粒子分辨率下角动量
${M_A}$ 的时历曲线Figure 3. Time evolution of the angular momentum
${M_A}$ with different particle resolutions![]()
图 5 αs = 2.0, 不同粒子分辨率时机械能衰减率
${E_{per}}$ 时历曲线Figure 5. Time evolution of the decay rate of mechanical energy at different particle resolutions with αs = 2.0
![]()
图 6 αs = 1.3(上)和αs = 2.0(下)时, 振荡液滴形态: 理论值(虚线)和SPH结果(压力云图)
Figure 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)
![]()
图 7 δ-SPH与δ-SPHC动能比较: 总体图(上)与局部放大图(下), 其中黑矩形框为放大区域
Figure 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
![]()
图 9 不同算法和光滑长度系数下, 振荡液滴机械能衰减率时历曲线
Figure 9. Time evolution of the decay rate of mechanical energy for the oscillating droplet with different SPH models and different smoothing length coefficient
![]()
图 11 造波推板横向位移随时间变化: (a)周期为T = 1 s的规则波, (b)不规则波
Figure 11. Paddle motions of the regular wave with (a) period T = 1 s and (b) irregular wave
![]()
图 12 SPH数值水池示意图, 红竖线表示波高仪, 选取距离波高仪一个粒子间距内自由面粒子的最大高度为波面高度
Figure 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
![]()
图 13
${\alpha _s}$ =2.0时, 不同粒子分辨率下波面高度的δ-SPH模拟值与试验值对比Figure 13. Comparison of the wave elevations between δ-SPH results and experimental data, SPH results with different particle resolutions and
${\alpha _s} = 2.0$ are provided![]()
图 14 在
${\alpha _s}$ = 1.3,$d/\Delta x$ = 60条件下, 传统δ-SPH模拟的波面形态Figure 14. Wave surface simulated by δ-SPH method with
${\alpha _s}$ = 1.3 and$d/\Delta x$ = 60![]()
图 15 在
${\alpha _s} = 1.3$ ,$d/\Delta x$ = 60条件下, 传统δ-SPH模型计算的不同位置波面高度时历曲线Figure 15. Wave elevation probed at different positions simulated by δ-SPH method with
${\alpha _s} = 1.3$ and$d/\Delta x$ = 60![]()
图 16 在
${\alpha _s} = 2.0$ ,$d/\Delta x$ = 60条件下, 传统δ-SPH模拟的波面形态Figure 16. Wave surface simulated by δ-SPH method with
${\alpha _s}$ = 2.0 and$ d/\Delta x $ = 60![]()
图 17 在
${\alpha _s}$ = 2.0,$ d/\Delta x $ = 60条件下, 传统δ-SPH模拟的不同位置波面高度时历曲线Figure 17. Wave elevation probed at different positions simulated by δ-SPH method with
${\alpha _s}$ = 2.0 and$d/\Delta x$ = 60![]()
图 18 在αs = 1.3和
$d/\Delta x$ = 15, 30, 60条件下, δ-SPHC计算的波面高度时历曲线与试验值对比Figure 18. Comparison between δ-SPHC results and experimental data for wave elevation: SPH results with αs = 1.3 and
$d/\Delta x$ = 15, 30, 60![]()
图 19 在αs = 2.0和
$d/\Delta x$ =15, 30, 60条件下, δ-SPHC计算的波面高度时历曲线与试验值对比Figure 19. Comparison between δ-SPHC results and experimental data for wave elevation: SPH results with αs = 2.0 and
$d/\Delta x$ = 15, 30, 60![]()
图 20
${x_{{\rm{probe}}}}$ =2.37 m处, 不同SPH算法下的波面高度数值结果与试验值对比Figure 20. Comparison between different SPH simulations and experimental data for wave elevations at
${x_{{\rm{probe}}}}$ = 2.37 m![]()
图 21 规则波压力云图(δ-SPHC,
${\alpha _s}$ = 2.0,$d/\Delta x$ = 60)Figure 21. Pressure contour of regular wave (δ-SPHC,
${\alpha _s}$ = 2.0,$d/\Delta x$ = 60)![]()
图 22
${\alpha _s}$ =1.3,$d/\Delta x$ = 60时, δ-SPHC算法模拟的波面形态Figure 22. The wave surface simulated by δ-SPHC with
${\alpha _s}$ = 1.3 and$d/\Delta x$ = 60![]()
图 23
${\alpha _s}$ =1.3,$d/\Delta x$ =60时, δ-SPHC模拟的不同位置波面高度时历曲线Figure 23. Wave elevations probed at different positions simulated by δ-SPHC method with
${\alpha _s}$ = 1.3 and$d/\Delta x$ = 60![]()
图 24
${\alpha _s}$ =2.0,$d/\Delta x$ =60时, δ-SPHC算法模拟的波面形态Figure 24. The wave surface simulated by δ-SPHC with
${\alpha _s}$ = 2.0 and$d/\Delta x$ =60![]()
图 25
${\alpha _s}$ = 2.0,$d/\Delta x$ = 60时, δ-SPHC模拟的不同位置波面高度时历曲线Figure 25. Wave elevations probed at different positions simulated by δ-SPHC method with
${\alpha _s}$ = 2.0 and$d/\Delta x$ = 60![]()
图 26 在
${x_{{\rm{probe}}}}$ = 6.37 m处, 规则波波面高度数值模拟结果与试验结果Figure 26. Time evolution of the regular wave elevation at
${x_{{\rm{probe}}}}$ = 6.37 m: SPH results and experimental data![]()
图 27 δ-SPHC在不同黏性系数α时,
${x_{{\rm{probe}}}}$ = 6.37 m处波面高度时历曲线Figure 27. Time history of wave elevation at
${x_{{\rm{probe}}}}$ = 6.37 m simulated by δ-SPHC with different α![]()
图 28
${\alpha _s}$ =1.3和不同粒子分辨率条件下, 不规则波δ-SPHC模拟结果的收敛性分析Figure 28. Convergence analysis of irregular waves simulated by δ-SPHC with αs=1.3 at different particle resolutions
![]()
图 29 不规则波压力云图(δ-SPHC,
${\alpha _s}$ =2.0,$d/\Delta x$ =60)Figure 29. Pressure contour of irregular wave (δ-SPHC,
${\alpha _s}$ =2.0,$d/\Delta x$ =60)![]()
图 30
$ x_{\text {probe }} $ = 6.37处, 不规则波的波面高度时历曲线Figure 30. Time evolution of the irregular wave elevation measured at
${x_{{\rm{probe}}}}$ = 6.37 m表 1 δ-SPH算法与δ-SPHC算法计算效率比较
Table 1 Comparison of computational efficiency between δ-SPH and δ-SPHC schemes
Scheme CPU time
of single step/st2/t1 δ-SPH t1 = 0.5963 1.09 δ-SPHC t2 = 0.6576 1.09 
下载: 导出CSV
表 2 不同算法在
${{\boldsymbol{x}}_{{\bf{probe}}}}$ = 6.37 m 处平均波高对比Table 2 Comparison of average wave heights probed at 6.37 m simulated with different methods
Method Exp. δ-SPH, αs = 1.3 δ-SPH, αs = 2.0 δ-SPHC, αs = 1.3 δ-SPHC, αs = 2.0 wave height H/m 0.0428 0.0224 0.0397 0.0419 0.0427 relative error $\varepsilon $/% − 47.66 7.24 2.10 0.23
下载: 导出CSV
-
[1] Luo M, Khayyer A, Lin P. Particle methods in ocean and coastal engineering. Applied Ocean Research, 2021, 114: 102734 doi: 10.1016/j.apor.2021.102734
[2] Grilli ST, Vogelmann S, Watts P. Development of a 3D numerical wave tank for modeling tsunami generation by underwater landslides. Engineering Analysis with Boundary Elements, 2002, 26(4): 301-313 doi: 10.1016/S0955-7997(01)00113-8
[3] Sung HG. BEM computations of 3D fully nonlinear free-surface flows caused by advancing surface disturbances. International Journal of Offshore and Polar Engineering, 2008, 18(4): 292-301
[4] Baudic SF, Williams AN, Kareem A. A two-dimensional numerical wave flume—Part 1: Nonlinear wave generation, propagation, and absorption. Journal of Offshore Mechanics and Arctic Engineering, 2001, 123(2): 70-75 doi: 10.1115/1.1365117
[5] 曹洪建, 万德成. 基于naoe-FOAM-SJTU求解器构建三维数值波浪水池. 复旦学报(自然科学版), 2013, 52(5): 627-634 (Cao Hongjian, Wan Decheng. Three-dimensional numerical wave tank based on naoe-FOAM-SJTU solver. Journal of Fudan University (Natural Science) , 2013, 52(5): 627-634 (in Chinese) Cao Hongjian, Wan decheng. Three-dimensional numerical wave tank based on naoe-FOAM-SJTU solver. Journal of Fudan University (Natural Science), 2013, 52(05): 627-634 (in Chinese)
[6] 董志, 詹杰民. 基于VOF方法的数值波浪水槽以及造波、消波方法研究. 水动力学研究与进展: A辑, 2009, 24(1): 15-21 (Dong Zhi, Zhan Jiemin. Comparison of existing methods for wave generating and absorbing in VOF-based numerical tank. Chinese Journal of Hydrodynamics, 2009, 24(1): 15-21 (in Chinese) Dong Zhi, Zhan Jiemin. Comparison of existing methods for wave generating and absorbing in VOF-based numerical tank. Chinese Journal of Hydrodynamics, 2009, 24(1): 15-21 (in Chinese)
[7] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics, 1981, 39(1): 201-225 doi: 10.1016/0021-9991(81)90145-5
[8] Huang C, Zhang DH, Si YL, et al. Coupled finite particle method for simulations of wave and structure interaction. Coastal Engineering, 2018, 140: 147-160
[9] Ye T, Pan D, Huang C, et al. Smoothed particle hydrodynamics (SPH) for complex fluid flows: recent developments in methodology and applications. Physics of Fluids, 2019, 31(1): 011301 doi: 10.1063/1.5068697
[10] Zhang A, Sun P, Ming F, et al. Smoothed particle hydrodynamics and its applications in fluid-structure interactions. Journal of Hydrodynamics, Ser. B, 2017, 29(2): 187-216 doi: 10.1016/S1001-6058(16)60730-8
[11] Liu M, Zhang ZL. Smoothed particle hydrodynamics (SPH) for modeling fluid-structure interactions. Science China: Physics, Mechanics and Astronomy, 2019, 62(8): 5-42
[12] 张华, 邵颂东. 水流结构相互作用的粒子法数值仿真. 水利水电技术, 2006, 9: 44-47 Zhang Hua, Shao Songdong. Numerical simulation for fluid-structure interaction with particle method. Water Resources and Hydropower Engineering, 2006, 9: 44-47 (in Chinese)
[13] Guilcher PM, Ducorzet G, Alessandrini B, et al. Water wave propagation using SPH models//Proceedings 2nd International Spheric Workshop, 2007: 119-122
[14] Omidvar P, Norouzi H, Zarghami A. Smoothed particle hydrodynamics for water wave propagation in a channel. International Journal of Modern Physics C, 2015, 26(8): 1550085 doi: 10.1142/S0129183115500850
[15] Antuono M, Colagrossi A, Marrone S, et al. Propagation of gravity waves through an SPH scheme with numerical diffusive terms. Computer Physics Communications, 2011, 182(4): 866-877
[16] Colagrossi A, Souto-Iglesias A, Antuono M, et al. Smoothed-particle-hydrodynamics modeling of dissipation mechanisms in gravity waves. Physical Review E, 2013, 87(2): 023302 doi: 10.1103/PhysRevE.87.023302
[17] Zhang DH, Shi YX, Huang C, et al. SPH method with applications of oscillating wave surge converter. Ocean Engineering, 2018, 152: 273-285 doi: 10.1016/j.oceaneng.2018.01.057
[18] Macià F, Colagrossi A, Antuono M, et al. Benefits of using a Wendland kernel for free-surface flows//Proceedings of 6th Ercoftac Spheric Workshop on SPH Applications, 2011: 30-37
[19] Meng Z, Zhang A, Wang P, et al. A targeted essentially non-oscillatory (TENO) SPH method and its applications in hydrodynamics. Ocean Engineering, 2022, 243: 110100 doi: 10.1016/j.oceaneng.2021.110100
[20] Liu GR, Liu MB. Smoothed Particle Hydrodynamics: A Meshfree Particle Method. World Scientific, 2003
[21] Bonet J, Lok T. Variational and momentum preservation aspects of smooth particle hydrodynamic formulations. Computer Methods in Applied Mechanics & Engineering, 1999, 180(1-2): 97-115
[22] Wen H, Ren B, Yu X. An improved SPH model for turbulent hydrodynamics of a 2D oscillating water chamber. Ocean Engineering, 2018, 150: 152-166 doi: 10.1016/j.oceaneng.2017.12.047
[23] Zago V, Schulize LJ, Bilotta G, et al. Overcoming excessive numerical dissipation in SPH modeling of water waves. Coastal Engineering, 2021, 170: 104018 doi: 10.1016/j.coastaleng.2021.104018
[24] Marrone S, Antuono M, Colagrossi A, et al. δ-SPH model for simulating violent impact flows. Computer Methods in Applied Mechanics and Engineering, 2011, 200(13-16): 1526-1542 doi: 10.1016/j.cma.2010.12.016
[25] Liu MB, Liu GR. Smoothed particle hydrodynamics (SPH): An overview and recent developments. Archives of Computational Methods in Engineering, 2010, 17(1): 25-76 doi: 10.1007/s11831-010-9040-7
[26] 顾声龙, 吴玉帅, 解宏伟等. 基于CSPM方法对二维管嘴出流的数值模拟. 水利水电技术, 2016, 47(9): 39-43 (Gu Shenglong, Wu Yushuai, Jie Hongwei, et al. CSPM method-based numerical simulation on outflow from 2-D pipe nozzle. Water Resources and Hydropower Engineering, 2016, 47(9): 39-43 (in Chinese) Gu Shenglong, Wu Yushuai, Jie Hongwei et al. . CSPM method-based numerical simulation on outflow from 2-D pipe nozzle. Water Resources and Hydropower Engineering, 2016, 47(09): 39-43 (in Chinese)
[27] 刘谋斌, 宗智, 常建忠. 光滑粒子动力学方法的发展与应用. 力学进展, 2011, 41(2): 217-234 Liu Moubin, Zong Zhi, Chang Jianzhong. Developments and applications of smoothed particle hydrodynamics, Advances in Mechanics, 2011, 41(02): 217-234 (in Chinese)
[28] Sun PN, Le Touze D, Oger G, et al. An accurate FSI-SPH modeling of challenging fluid-structure interaction problems in two and three dimensions. Ocean Engineering, 2021, 221: 108552 doi: 10.1016/j.oceaneng.2020.108552
[29] Oger G, Doring M, Alessandrini B, et al. An improved SPH method: Towards higher order convergence. Journal of Computational Physics, 2007, 225(2): 1472-1492 doi: 10.1016/j.jcp.2007.01.039
[30] Ganzenmüller GC, Hiermaier S, May M. On the similarity of meshless discretizations of peridynamics and smooth-particle hydrodynamics. Computers & Structures, 2015, 150: 71-78
[31] Lyu HG, Deng R, Sun PN, et al. Study on the wedge penetrating fluid interfaces characterized by different density-ratios: Numerical investigations with a multi-phase SPH model. Ocean Engineering, 2021, 237: 109538
[32] 杨秋足, 徐绯, 王璐等. 一种基于黎曼解处理大密度比多相流SPH的改进算法. 力学学报, 2019, 51(3): 730-742 Yang Qiuzu, Xu Fei, Wang Lu, et al. An improved SPH algorithm for large density ratios multiphase flows based on riemann solution, Chinese Journal of Theoretical and Applied 2019, 51(3): 730-742 (in Chinese)
[33] Antuono M, Marrone S, Colagrossi A, et al. Energy balance in the δ-SPH scheme. Computer Methods in Applied Mechanics and Engineering, 2015, 289: 209-226 doi: 10.1016/j.cma.2015.02.004
[34] Colagrossi A, Bouscasse B, Antuono M, et al. Particle packing algorithm for SPH schemes. Computer Physics Communications, 2012, 183(8): 1641-1653 doi: 10.1016/j.cpc.2012.02.032
[35] Antuono M, Colagrossi A, Marrone S, et al. Free-surface flows solved by means of SPH schemes with numerical diffusive terms. Computer Physics Communications, 2010, 181(3): 532-549 doi: 10.1016/j.cpc.2010.12.012
[36] Cui J, Chen X, Sun P. Numerical investigation on the hydrodynamic performance of a new designed breakwater using smoothed particle hydrodynamic method. Engineering Analysis with Boundary Elements, 2021, 130: 379-403 doi: 10.1016/j.enganabound.2021.05.007
[37] Adami S, Hu XY, Adams NA. A generalized wall boundary condition for smoothed particle hydrodynamics. Journal of Computational Physics, 2012, 231(21): 7057-7075 doi: 10.1016/j.jcp.2012.05.005
-
期刊类型引用(1)
1. 管祥善，孙鹏楠，李江昊，孙龙泉. 基于光滑粒子流体动力学的波浪中航行体入水数值模拟. 空气动力学学报. 2024(02): 85-95 . 
百度学术
其他类型引用(2)
计量
- 文章访问数: 1084
- HTML全文浏览量: 337
- PDF下载量: 274
- 被引次数: 3
