﻿ 砾石单元法与SPH耦合模型在数值波浪水槽中的应用研究
 海洋科学  2020, Vol. 44 Issue (9): 100-111 PDF
http://dx.doi.org/10.11759/hykx20191018001

#### 文章信息

LI Yan-qing, BIE She-an, LI Da-ming, WANG Xin. 2020.

Application research in channel waves with SPH and GEM coupling model

Marina Sciences, 44(9): 100-111.
http://dx.doi.org/10.11759/hykx20191018001

### 文章历史

1. 天津大学水利工程仿真与安全国家重点实验室, 天津 300072;
2. 国家海洋技术中心, 天津 300112

Application research in channel waves with SPH and GEM coupling model
LI Yan-qing1, BIE She-an1, LI Da-ming1, WANG Xin2
1. State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin University, Tianjin 300072, China;
2. National Ocean Technology Center, Tianjin 300112, China
Abstract: To study the large deformation of the biphasic medium, whose fluid particles are used to describe the movement of the waves and Gravel is described by a solid element, a mathematical model of SPH-GEM coupling was established. SPH (Hydrodynamic Smooth Particles Method) is applied to establish a numerical wave flume, and the Gravel Element Method (GEM) is used to simulate the roll deformation and collapse of the Gravel under the action of waves. The GEM proposed in this article has the function of simulating the movement of single-phase gravel piling and coupling the movement of gravel and fluid particles in multiphase media. It is a complement of the Discrete Element Method (DEM). The article presents the basic principle of SPH and GEM, the similarities and differences between GEM and DEM. The SPH coupling wave flume with the GEM coupling method was explained and the mathematical model was established to simulate the water flume under the action of the waves. The wave generation process and the gravel fall, the result of the collapse deformation are compared with the results of the physical experiment. The results are basically consistent. The principle, method and application of the mathematical model coupled to SPH-GEM provide a new idea for researching a similar mathematical model.
Key words: SPH method    GEM    DEM    wave numerical flume

1 堆积砾石的受力分析 1.1 单砾石几何形态

 $R(\theta , \varphi ) = \frac{{{a_0}}}{2} + \sum\limits_{n = 1}^\infty {{{[{a_n}\cos (n\varphi ) + {b_n}\sin (n\varphi )]}_\theta }} ,$ (1)

 $\left\{ {\begin{array}{*{20}{c}} {{a_n} = \frac{2}{T}\int_{ - \frac{T}{2}}^{\frac{T}{2}} {R(\theta , \varphi )\cos (n\varphi ){\rm{d}}\varphi } }\\ {{b_n} = \frac{2}{T}\int_{ - \frac{T}{2}}^{\frac{T}{2}} {R(\theta , \varphi )\sin (n\varphi ){\rm{d}}\varphi } } \end{array}} \right.,$ (2)

anbn均为零时, 砾石颗粒半径为, a0为圆形砾石。低频谐波(n较小)反映砾石颗粒表面较大尺度的形态变化, 高频谐波(n较大)反映砾石颗粒表面较小尺度的形态变化。谐波越多, 傅里叶级数生成的轮廓线越接近于实际砾石颗粒的轮廓线。图 1a为砾石颗粒的图像采集, 图 1b为砾石颗粒的轮廓模拟。

 图 1 砾石颗粒图像和轮廓模拟 Fig. 1 Image of gravel particles and contour simulation 注: a:原始颗粒; b:颗粒轮廓模拟

 $\psi = \sqrt[3]{{{{\left( {\frac{b}{a}} \right)}^2}\frac{c}{b}}},$ (3)

 $\phi = \sqrt {\frac{b}{a}{{\left( {\frac{c}{b}} \right)}^2}} .$ (4)

 图 2 非圆颗粒的组合形状 Fig. 2 Combined shape of non-circular particles

 图 3 砾石的二维圆周矢量半径分布 Fig. 3 Distribution of the 2D circumferential vector radius of the gravel
1.2 单砾石力学球概念

 砾石形状 球度系数 沙粒种类 球度系数 球 1.000 河沙 0.860 六方八面体 0.906 海沙 0.860 正八面体 0.846 煤粉 0.700 立方体 0.806 石英沙 0.670~0.600 正四面体 0.670 无烟煤 0.610~0.400 圆柱体微粒 0.860 焦炭 0.350 页岩屑 0.320~0.290

 $f(\theta ) = {f_{r0}} + {f_r}F(\theta ),$ (5)

 $J(\theta ) = \frac{2}{5}m{r^2} + {J_s}{F_s}(\theta ),$ (6)

1.3 堆积砾石的受力状态

 图 4 砾石受力示意图 Fig. 4 Gravel strength diagram

 图 5 力学球砾石堆 Fig. 5 Mechanical ball gravel pile

 图 6 按求解层次编排的级别号 Fig. 6 A level number organized by the solution level

 图 7 按级别求出的压力分布 Fig. 7 Pressure distribution by grade

2 砾石单元法(GEM) 2.1 GEM模式

GEM是以独立的砾石整体作为一个刚性单元, 砾石表面的轮廓线是单元边界, 砾石以整体运动, 具有平动和转动的形式。砾石间可以通过接触传递力和力矩, 具有运动速度的砾石可以传递动量。在GEM模式中引入力学球的概念, 不考虑砾石的形状变化; 力学性质可以分布在力学球表面, 使力学球具有静、动摩擦系数, 不同表面位置具有不同的转动惯量。力学球之间不能压入和穿越, 运动受拟序排列控制, 将砾石刚体碰撞效果代替微弹性假定, 力用接触压力和反力模拟, 力矩用接触摩擦和重力对支撑点之矩模拟[15-16]。运动受力球动力传递模式如图 8所示。

 图 8 砾石力学球动力传递示意图 Fig. 8 Schematic of the dynamic ball transmission of gravel mechanics 注: a:接触模式; b:转动模式

2.2 GEM主要力学公式

 $\left\{ \begin{array}{l} {S_{ijx}} = {m_i}{v_x} - {m_j}{v_{x0}}\\ {S_{ijy}} = {m_i}{v_y} - {m_j}{v_{y0}} \end{array} \right.,$ (7)

 ${F_i} = \sum\limits_{j = 1}^n {\left( {f{n_{ij}} + f{\tau _{ij}}} \right)} ,$ (8)
 ${T_i} = \sum\limits_{j = 1}^n {{r_i}\left( {{f_{{\rm{s}}ij}} \cdot {n_i}} \right)} ,$ (9)

 $\left\{ \begin{array}{l} {f_1} = - {f_2}\\ {M_1} = - {M_2} \end{array} \right.,$ (10)

 ${R_{ij}} \ge \left( {{d_i} + {d_j}} \right)/2,$ (11)

2.3 GEM与DEM比较

GEM与DEM主要有以下几点不同:

(1) 砾石单元法以砾石为整体刚性单元, 主要以硬球的接触模式为主, 忽略砾石表面的变形细节; 离散单元法以粒子形式分布单元, 主要考虑软球模式的接触方式, 把颗粒间的法向力简化为弹簧和阻尼器。

(2) 砾石单元法以堆积砾石为研究对象; 离散单元法可以模拟固体颗粒群的运动, 一般不考虑砾石承压和遮蔽引起的运动次序的差别。

(3) 砾石单元法是在堆积砾石稳定性研究基础上提出来的, 在力学模式上具有外部失稳但内部稳定的物理意义; 离散单元法以模拟固体颗粒群体运动的效果为基础, 具有整体失稳的物理意义。

3 SPH与GEM耦合波浪砾石模型 3.1 SPH波浪水槽模型及其控制方程

SPH控制方程严格遵循质量守恒、动量守恒及能量守恒这三条物理守恒定律, 将描述流体运动的Navier-Stokes方程, 通过SPH法基本理论的思路对N-S方程进行空间离散化, 得到了在笛卡尔坐标系下并且适用于广义流体动力学的SPH控制方程。SPH中最常用的连续方程粒子近似式为:

 $\frac{{{\mathop{\rm d}\nolimits} {\rho _i}}}{{{\mathop{\rm d}\nolimits} t}} = \sum\limits_{j = 1}^{{N_j}} {{m_j}{v_{ij}} \cdot {\nabla _i}{W_{ij}}} ,$ (12)

 $\begin{array}{c} \frac{{{\mathop{\rm d}\nolimits} {v_i}}}{{{\mathop{\rm d}\nolimits} t}} = - \sum\limits_{j = 1} {{m_j}\left( {\frac{{{p_j}}}{{\rho _j^2}} + \frac{{{p_i}}}{{\rho _i^2}} + {\prod _{ij}}} \right){\nabla _i}{W_{ij}}} + \\ \;\;\sum\limits_{j = 1} {{m_j}\left( {\frac{{4{\upsilon _0}{x_{ij}} \cdot {\nabla _i}{W_{ij}}}}{{({\rho _i} + {\rho _j})r_{ij}^2}}} \right){v_{ij}}} , \end{array}$ (13)

 $p = B\left[ {{{\left( {\frac{\rho }{{{\rho _0}}}} \right)}^\gamma } - 1} \right],$ (14)

 图 9 数值波浪水槽模拟的造波过程和压力分布 Fig. 9 Wave-making process and pressure distribution simulated using the wave numeric channel 注:子图a—h为模型历时0.0~4.0 s变化过程

 图 10 波浪水槽数值模拟与实测数据对比 Fig. 10 Numerical simulation of the wave channel compared with the measured data
3.2 SPH与GEM耦合模型数学基础

SPH与GEM耦合模型由SPH粒子和堆积砾石构成, SPH粒子运动由方程(12)控制, 堆积砾石处为SPH粒子的固体边界, SPH粒子运动速度为

 $\left\{ \begin{array}{l} {v_n}(x,y) = 0{\rm{ }}(x,y) \in \mathit{\Gamma }\\ v(x,y) = 0{\rm{ (}}x,y{\rm{)}} \in \mathit{\Omega } \end{array} \right.,$ (15)

SPH粒子对砾石的作用力为沿砾石表面的法线方向压力的积分, 利用势流的伯努利方程为

 $\vec F{\rm{ = }}\int\limits_\Gamma { - \vec p} \cdot \vec n{\mathop{\rm d}\nolimits} s = \int_0^{2\pi } { - \left[ {{p_\infty } + \frac{1}{2}\rho \left( {v_\infty ^2 - v_\theta ^2} \right)} \right]{\mathop{\rm d}\nolimits} s} ,$ (16)

 $\left\{ \begin{array}{l} \frac{{{\mathop{\rm d}\nolimits} x}}{{{\mathop{\rm d}\nolimits} t}}{\rm{ = }}{v_{x0}} + \frac{{\sum {{F_x}} }}{m}(t - {t_0})\\ \frac{{{\mathop{\rm d}\nolimits} y}}{{{\mathop{\rm d}\nolimits} t}}{\rm{ = }}{v_y}_0 + \frac{{\sum {{F_y}} }}{m}(t - {t_0})\\ \frac{{{\mathop{\rm d}\nolimits} \theta }}{{{\mathop{\rm d}\nolimits} t}}{\rm{ = }}{\omega _0} + \frac{{\sum M }}{J}(t - {t_0}) \end{array} \right.,$ (17)

3.3 SPH与GEM耦合模型模拟过程

SPH与GEM的耦合过程主要有以下几点:

(1) 应用坝体轮廓线控制将具有级配的堆积砾石随机放置在SPH数值水槽中, 确定每一块砾石轮廓线内受控制的SPH粒子, 这些粒子称为受控粒子, 不包含SPH粒子的砾石是水面以上砾石或尺度较小的砾石, 不包含SPH粒子的砾石不会影响砾石运动计算模拟。受控粒子与非受控粒子分布见图 11

 图 11 受控粒子与非受控粒子分布 Fig. 11 Distribution of controlled and uncontrolled particles

(2) 在SPH模型边界上使用推板造波, 砾石中的受控SPH粒子参与核函数计算, 但运动速度为零, 相当于将砾石对周边SPH粒子的影响分解为受控粒子对周边SPH粒子的影响, 体现了砾石边界轮廓对SPH模式的耦合效果。

(3) 同步统计砾石的受力和周边的流速场影响(见图 12), 砾石受到的总压力可以用砾石中受控粒子的总压力代替, 受控粒子(图 12黑轮廓线内粒子)的压力分布决定了砾石运动形态。砾石周边流场取砾石轮廓以外、砾石光滑长度以内的粒子做统计分析, 图 12中黑色砾石轮廓外与红色光滑长度圆内的粒子, 黑色矢量表示砾石光滑长度, 本文取2倍的力学球半径, 这一区域的流速、流向可以处理为波浪对砾石的冲击效果, 砾石是否进入运动状态, 这也是重要判别条件之一。

 图 12 砾石光滑长度和控制区域 Fig. 12 Gravel length smooth and control area

(4) 在堆积砾石的受力分析中, 用拟序排列、虚位移原理确定砾石的承压力、反力、摩擦力和挤压力, 确定砾石进入运动状态的级别, 通常级别的确定会在每次计算前重新划分, 因为每一次砾石运动都会改变砾石堆的排列形式, 一般一次计算中只需考虑1级砾石的运动, 其他级别得不到运动的机会。

(5) 波浪力不足以撼动砾石时, 砾石保持原有状态; 波浪力撼动砾石时, 要考虑砾石的反力分布确定砾石运动形式, 一般情况为没有反力作用时为塌落, 有支撑点时要找出支撑点砾石。

(6) 砾石运动和就位应在瞬时完成, 所以在一个时间步长内就应当确定出砾石运动后的位置, 这一位置应避免与其它砾石压入, 并位于砾石运动的路径上; 砾石进入新的位置后将会放弃以前的受制粒子, 同时产生出新的受制粒子, 完成了一个计算的循环过程。

3.4 SPH与GEM耦合模型模拟结果

 图 13 防波堤模型坡面颗粒变化过程 Fig. 13 Process of changing the grain slope in the breakwater model 注:子图a—j为模型历时0.5~5.0 s变化过程

 图 14 防波堤模型初始状态 Fig. 14 Initial state of the breakwater model

 图 15 防波堤模型稳定状态 Fig. 15 Steady state of the breakwater modelf

 图 16 防波堤模型断面变化对比 Fig. 16 Comparison of the cross section of the breakwater model
3.5 耦合模型的可靠性验证

 图 17 文献实验布局示意图(单位: cm) Fig. 17 Literature experiment layout

 图 18 数值模型计算示意图 Fig. 18 Schematic of the numerical model calculation

 图 19 测点波高对比图 Fig. 19 Contrast diagram of the measured point height wave 注:子图a: 1号测点; b: 2号测点; c: 3号测点; d: 4号测点
4 波浪砾石物理模型试验 4.1 物理模型试验

 图 20 波浪试验水槽 Fig. 20 Experimental wave tank 注:子图a:试验示意图; b:试验布置图

 图 21 试验设计断面和试验设备在水槽中的分布 Fig. 21 Test design and distribution section of the test equipment in the water tank
4.2 物理模型试验结果

 图 22 波浪观察与波浪模拟的过程比较 Fig. 22 Comparison between wave observation and wave simulation

 图 23 防波堤断面1模型试验工况1动力平衡断面 Fig. 23 Section 1 of the breakwater model test Condition 1 dynamic balance section
5 结论

(1) 采用傅氏级数的方法对块石轮廓进行了曲线展开, 说明块石轮廓沿曲线的矢径分布情况, 通过分布砾石轮廓上力学性质, 可以对砾石轮廓赋予力学上的摩擦系数、转动惯量等物理量, 提出了力学球体的概念。

(2) 分析了堆积砾石的受力特点, 提出堆积砾石分级方法和拟序排列求解方法, 对静态砾石堆应用拟序排列求得的解, 与理论值误差很小, 说明拟序排列求解方法可行。

(3) 提出了GEM模型方法, 建立了SPH波浪水槽与GEM堆积砾石耦合数学模型, 模拟了砾石在波浪作用下的堤坝变形过程, 模拟结果与物理模型试验结果基本一致。

 [1] Monaghan J J. An introduction to SPH[J]. Comput Phys Commun, 1988, 48(1): 89-96. DOI:10.1016/0010-4655(88)90026-4 [2] Monaghan J J. Smoothed particle hydrodynamics[J]. Annual Rev Astron Appl, 1992, 30: 543-574. DOI:10.1146/annurev.aa.30.090192.002551 [3] Bai Ling, Li Daming, Li Yanqing, et al. Study on the droplet impact on hydrophobic surface in terms of van der Waals surface tension model[J]. Acta Physica Sinica, 2015, 64(11): 114701-114707. [4] Li Daming, Wang Zhichao, Bai Ling, et al. Investigations on the process of droplet impact on an oriﬁce plate[J]. Acta Phys Sin, 2013, 62(19): 194704. [5] Li Daming, Li Xiaoyu, Lin Yi. Numerical simulation of droplet impacting surface by SPH[J]. Science China Technological Sciences, 2011, 54(7): 1873-1880. DOI:10.1007/s11431-011-4422-0 [6] Swegle J W, Hicks D L, Attaways S W. Smoothed particle hydrodynamics stability analysis[J]. Journal of Computational Physics, 1995, 116: 123-134. DOI:10.1006/jcph.1995.1010 [7] Cundall P A. A computer model for simulating progressive large scale movements in blocky rock systems[J]. In Proc Syrup Int Rock Mech, 1971, 1: 2-8. [8] Cundall P A, Strack O D L. A discrete numerical model for granular assemblies[J]. Geotechnique, 1979, 29(1): 47-65. [9] 任冰, 金钊, 高睿, 等. 波浪与斜坡堤护面块体相互作用SPH-DEM数值模拟[J]. 大连理工大学学报, 2013, 53(2): 241-248. Ren Bing, Jin Zhao, Gao Rui, et al. SPH-DEM numerical simulation of the interaction between wave and slope embankment protection block[J]. Journal of Dalian University of Technology, 2013, 53(2): 241-248. [10] 王志超.基于SPH-DEM耦合方法的液滴冲击散粒体运动机理研究[D].天津: 天津大学, 2015. Wang Zhichao. SPH-DEM coupling method to study the motion mechanism of droplet impact granule[D]. Tianjin: Tianjin University, 2015. [11] 别社安, 姜长杰. 离散单元法与港口和海岸工程中的非连续介质问题[J]. 海洋通报, 1999, 18(1): 80-87. Bie Shean, Jiang Changjie. Discrete element method and non-continuum problems in port and coastal engineering[J]. Marine Science Bulletin, 1999, 18(1): 80-87. [12] Li Yanqing, Li Daming, Bie shean, et al. Numerical simulation for fluid droplet impact on discrete particles with coupled SPH-DEM method[J]. International Journal of Numerical Methods for Heat and Fluid Flow, 2018, 28(11): 2581-2605. DOI:10.1108/HFF-11-2017-0464 [13] Li Daming, Li Yangyang, Wang Zhichao, et al. Quantitative SEM-based shape analysis of sediment particles in the Yellow River[J]. International Journal of Sediment Research, 2016, 31(4): 341-350. DOI:10.1016/j.ijsrc.2016.05.006 [14] 钱宁. 泥沙运动力学[M]. 北京: 科学出版社, 1983: 14-21. Qian Ning. Mechanics of Sediment Transport[M]. Beijing: Science Press, 1983: 14-21. [15] Cundall P A. A computer model for simulating progressive large scale movements in blocky rock systems[J]. In Proc Syrup Int Rock Mech, 1971, 1: 2-8. [16] Cundall P A, Strack O D L. A discrete numerical model for granular assemblies[J]. Geotechnique, 1979, 29(1): 47-65. [17] 俞聿修. 随机波浪及其工程应用[M]. 大连: 大连理工大学出版社, 2000. Yu Yuxiu. Random Wave and Its Applications for Engineering[M]. Dalian: Dalian University of Technology Press, 2000.