海塘越浪过程SPH模拟

杨洲+颜君来+曾甄

摘 要：随着计算机软硬件水平的不断发展，数值模拟技术已经越来越广泛地应用到海岸工程领域中。现今波浪爬高越浪的数值模拟技术主要有网格方法和无网格方法两种。和网格方法相比，无网格方法在计算海塘越浪的连续性和大变形问题上体现出了优势。

关键词：越浪 SPH方法 数值模拟

浙江省地理位置特殊，海岸线绵长，台风灾害频繁，一旦海塘受损，将会造成不可估量的损失。根据以往溃坝的经验，越浪是造成海塘溃损的主要原因之一。随着计算机软硬件水平的不断发展，海塘越浪的研究当中不断地引入数值模拟技术。数值模拟海塘的越浪以及计算海塘的越浪量，通过更加符合实际的边界条件等计算的越浪量，比通过规范公式计算的越浪量更加符合工程安全要求，同时模拟的越浪过程可以直观形象地了解海塘越浪的情况；数值模拟不需要专门的试验室、仪器和相关的工作人员，计算周期要短得多，方便进行多次计算，相比传统的物理模型试验，数值模拟节约了大量的人力物力，同时在试验结果的精确度上又有一定的保证。

如今海塘越浪的数值模拟计算方法主要有两种：网格方法和无网格方法。 1977年Lucy、Gingold、Monaghan]等分别提出了光滑质点流体动力学（SPH）方法，近几年来，光滑粒子流体动力学（SPH）方法因为它的拉格朗日特性，在模拟自由表面流体的应用正变得热门，但目前而言在海塘（堤）越浪的定量问题应用研究还并不广泛。

本文基于DualSPHysics开源平台，建立二维数值水槽，模拟规则波在海塘上的爬坡和越浪，并与物理模型试验结果进行了对比，验证了数值模型的可靠性。再通过数值模型计算越浪量，通过和规范公式计算得到的以及物理模型试验得到的越浪量进行对比，实现SPH方法在工程实践中的应用。

1.SPH方法介绍

SPH方法有两个主要的步骤：第一步是积分表达区域函数的函数近似；第二个步骤就是粒子近似。通过函数近似，可以将任意场函数（质量场、温度场等）进行积分表达形式，并进行粒子近似，将在粒子i处的场函数积分表达式进一步离散为迭加形式。由连续性方程、动量方程和能量方程组成流体的基本控制方程即N-S方程，并分别通过流体的质量守恒、动量守恒和能量守恒原理推导得出。

2.物理模型验证

选取某工程一个断面进行物理模型试验和SPH数值模拟试验，然后对比两个试验模型中波浪的运动过程。当试验开始波浪沿着海塘护面向上传播，在经过海塘平台处时，由于水深变浅，波浪开始变形翻卷并破碎。水体冲击防浪墙，并沿防浪墙向上爬，由于防浪墙挡浪面的弧线形式，水体冲上一定高度开始翻卷并下落。接着挡浪墙前面的水体退去，高度降低，又开始重复上面的过程。通过对比物理模型和数值模型波浪的运动过程，可以发现两个过程基本一致。

在防浪墙前面选取测量点（如图1所示）测量水体高度，由于SPH模拟的越浪过程刚开始由于水体内部还没有产生持续且稳定的紊流使得前面几个波浪一般偏大，到第五个开始水体周期内的运动过程趋于一致，故这次测量从第5个波浪开始。测量18s的结果如图2所示。

通过观察上图，在所设测量点的水体高度变化中SPH模拟的值和物模试验中的值非常吻合，最大误差范围在16%以内，主要是由于SPH方法在模拟海塘越浪过程DualSPHysics对固液边界处理不够完善以及物理模型试验中存在一个人为误差。验证了SPH方法在海塘越浪过程模拟中的可靠性。

3.工程应用

3.1越浪过程模拟

通过dualSPHysics开源平台计算的海塘越浪过程大致可以分为3个阶段，分别为：爬坡阶段、飞溅阶段、越浪流过程。首先当水流上爬至海塘中段的平台，水流及其速度都比较稳定，再往上爬，波浪翻卷变形，水舌产生较大的速度。水流遇到防浪墙并冲击防浪墙，水流改变方向沿着防浪墙往上爬，水流在防浪墙的顶端处产生局部的速度峰值，但防浪墙未能阻挡全部水流；接着进入飞溅阶段，跃入空中的水流质点速度方向均朝着斜上方，并有水滴飞溅，之后虽然水流还在不断爬坡，但是水流的速度开始变小，这是仍在不断爬坡，但水流仍保持一个较小的指向斜上方的塘顶速度，在到达海塘上空的水流下降的同时，前坡上的水流汇合，这时在前堤上的大部分速度接近于零。之后前坡上的水流开始回落，此时坡前平台附近的水流速度刚刚降落至零点，回落水流与之相冲，形成一个明显的回落波形。最后是越浪流过程，越浪流过程和前面的飞溅过程时间上有部分重合，当波高较大时，有些水质点会直接落在后坡上，而下落的水流砸落在防浪墙顶上后逐渐形成越浪流这时水流将分成两部分，靠近防浪墙的水流将冲回原水流当中，并卷入气泡，靠近后坡的水流年则顺着塘顶向后坡流去，飞溅在空中的水花也会有部分回落这部分水流当中。

在消浪平台的左端、折角处、右端以及防浪墙前端分别设置水体高度测量点，记为1#、2#、3#和4#（如图3所示）。测量点水体高度测量值见表1，测量点水体高度变化图见图4。

1#测量点初段水体高度变化曲线较平时由于前一个波浪的回落力度和下一个波浪的雍高力度相当，是此处的水体高度能够保持一段时间稳定。接着随着波浪的进一步传播由于1#测量点为一个消浪平台，波浪的传播形态能够基本保持一个完整的波形，同时由于消浪平台的高程太高，会使水体有雍高现象，故水体高度的变化曲线的峰值比谷值大；2#测量点之后就是有坡度的消浪坡面，高程逐渐抬高，波浪到2#测量点会受到阻碍，水体高程较平稳。当水体回落，由于前后高程差较大，产生较大的谷值；3#测量点由于靠近挡浪墙，水体遇到挡浪墙会沿着防浪墙往上爬，高度升高比较快，同时也是因为挡浪墙的作用，雍高了水位，是水体高度下降较缓慢；4#测量点由于在挡浪墙前，水体高度变化图跟3#测量点相差不大，但是水体高度的变化不管是升高还是降低都平缓一点。

3.2越浪量计算

和物理模型一样，越浪过程中的前几个波浪往往偏大，本文考察第5个和第6个波浪作用过程，以此计算越浪量。无风条件下三种方式得到的越浪量比较见表2。

通过对比发现通过规范公式计算的越浪量值比物理模型和SPH模型得到的值都大，这是因为物理模型和SPH模拟试验采用的是更加符合实际情况的边界条件，而计算公式中没有考虑反弧挑浪式防浪墙的挑浪作用和消浪平台的消浪作用，因此认为物理模型试验和SPH模拟试验的结果更符合真实情况。

通过SPH方法得到的越浪量几乎都略大于物理试验值，比如本文SPH方法所模拟的结果值比物理模型试验值要大大约29%，造成这样误差的原因是建立数值模型的时候简化了消浪平台，使其消浪作用减弱，另外本文采用设置一个水槽来接水最后按照所接水体体积来计算越浪量，这个方法存在人为误差，另一个计算误差是Dual SPHysics对于固液边界处理不够完善，首先，液体粒子和固壁粒子之间并未添加摩阻机制，使得波浪在斜坡堤的爬坡过程中波高偏大。其次，在使用连续密度法求解粒子密度导数时，将粒子相对速度引入了算式。由于固壁粒子和液体粒子的相对速度要明显大于液体粒子之间的相对速度，直接导致边界粒子的密度和压强大于正常值，液体粒子收排斥而在固壁粒子之间逐渐形成无粒子层。该无粒子层的厚度约在2个光滑长度左右，起到一种类似“气垫”的效果，将越浪水流托起，使其能更轻松地越过堤顶，如图4所示。

4.结论

本文基于DualSPHysics开源平台，采用光滑质点流体动力学（SPH）进行尝试，通过选取某工程一个断面的物理模型试验的结果和SPH方式模拟的海塘挡浪过程以及预设在防浪墙前的水体高度测量点位置的水体高度变化过程的对比来验证SPH方法在水动力模型试验中的可靠性，最后把SPH方法应用在漩门三期工程中，模拟海塘的越浪过程并验证越浪量的准确性。

本文SPH方法在海塘越浪中作了初步的研究，为保证该方法的稳定性和可靠性，仍需进行大量的测试和改进工作。另外，SPH方法的基础理论尚存很多待解决的问题，特别是在固壁粒子与液体粒子相互作用方面，亟需引入一种更为有效的处理方法，以便实现准确而快速的边界压强计算。

参考文献：

[1] Lucy L B. A Approach to the Testing of the Fission Hypothesis. The Astron. J. ， 1977， 8 （12）： 1013- 1024p.

[2]Gingold RA， Monaghan JJ. Smoothed Partide Hydrodynamics： Theory and Applications to Non-spherical Stars. Mon. Not. Roy. Astrou. Soc. ， 1977， 18： 375 -389P.

[3] Monaghan J J. Smoothed particle hydrodynamics[j]. Reports On Progress In Physics. 2005， 68（8）： 1703-1759.

通过对比发现通过规范公式计算的越浪量值比物理模型和SPH模型得到的值都大，这是因为物理模型和SPH模拟试验采用的是更加符合实际情况的边界条件，而计算公式中没有考虑反弧挑浪式防浪墙的挑浪作用和消浪平台的消浪作用，因此认为物理模型试验和SPH模拟试验的结果更符合真实情况。

通过SPH方法得到的越浪量几乎都略大于物理试验值，比如本文SPH方法所模拟的结果值比物理模型试验值要大大约29%，造成这样误差的原因是建立数值模型的时候简化了消浪平台，使其消浪作用减弱，另外本文采用设置一个水槽来接水最后按照所接水体体积来计算越浪量，这个方法存在人为误差，另一个计算误差是Dual SPHysics对于固液边界处理不够完善，首先，液体粒子和固壁粒子之间并未添加摩阻机制，使得波浪在斜坡堤的爬坡过程中波高偏大。其次，在使用连续密度法求解粒子密度导数时，将粒子相对速度引入了算式。由于固壁粒子和液体粒子的相对速度要明显大于液体粒子之间的相对速度，直接导致边界粒子的密度和压强大于正常值，液体粒子收排斥而在固壁粒子之间逐渐形成无粒子层。该无粒子层的厚度约在2个光滑长度左右，起到一种类似“气垫”的效果，将越浪水流托起，使其能更轻松地越过堤顶，如图4所示。

4.结论

本文基于DualSPHysics开源平台，采用光滑质点流体动力学（SPH）进行尝试，通过选取某工程一个断面的物理模型试验的结果和SPH方式模拟的海塘挡浪过程以及预设在防浪墙前的水体高度测量点位置的水体高度变化过程的对比来验证SPH方法在水动力模型试验中的可靠性，最后把SPH方法应用在漩门三期工程中，模拟海塘的越浪过程并验证越浪量的准确性。

本文SPH方法在海塘越浪中作了初步的研究，为保证该方法的稳定性和可靠性，仍需进行大量的测试和改进工作。另外，SPH方法的基础理论尚存很多待解决的问题，特别是在固壁粒子与液体粒子相互作用方面，亟需引入一种更为有效的处理方法，以便实现准确而快速的边界压强计算。

参考文献：

[1] Lucy L B. A Approach to the Testing of the Fission Hypothesis. The Astron. J. ， 1977， 8 （12）： 1013- 1024p.

[2]Gingold RA， Monaghan JJ. Smoothed Partide Hydrodynamics： Theory and Applications to Non-spherical Stars. Mon. Not. Roy. Astrou. Soc. ， 1977， 18： 375 -389P.

[3] Monaghan J J. Smoothed particle hydrodynamics[j]. Reports On Progress In Physics. 2005， 68（8）： 1703-1759.

通过对比发现通过规范公式计算的越浪量值比物理模型和SPH模型得到的值都大，这是因为物理模型和SPH模拟试验采用的是更加符合实际情况的边界条件，而计算公式中没有考虑反弧挑浪式防浪墙的挑浪作用和消浪平台的消浪作用，因此认为物理模型试验和SPH模拟试验的结果更符合真实情况。

通过SPH方法得到的越浪量几乎都略大于物理试验值，比如本文SPH方法所模拟的结果值比物理模型试验值要大大约29%，造成这样误差的原因是建立数值模型的时候简化了消浪平台，使其消浪作用减弱，另外本文采用设置一个水槽来接水最后按照所接水体体积来计算越浪量，这个方法存在人为误差，另一个计算误差是Dual SPHysics对于固液边界处理不够完善，首先，液体粒子和固壁粒子之间并未添加摩阻机制，使得波浪在斜坡堤的爬坡过程中波高偏大。其次，在使用连续密度法求解粒子密度导数时，将粒子相对速度引入了算式。由于固壁粒子和液体粒子的相对速度要明显大于液体粒子之间的相对速度，直接导致边界粒子的密度和压强大于正常值，液体粒子收排斥而在固壁粒子之间逐渐形成无粒子层。该无粒子层的厚度约在2个光滑长度左右，起到一种类似“气垫”的效果，将越浪水流托起，使其能更轻松地越过堤顶，如图4所示。

4.结论

本文基于DualSPHysics开源平台，采用光滑质点流体动力学（SPH）进行尝试，通过选取某工程一个断面的物理模型试验的结果和SPH方式模拟的海塘挡浪过程以及预设在防浪墙前的水体高度测量点位置的水体高度变化过程的对比来验证SPH方法在水动力模型试验中的可靠性，最后把SPH方法应用在漩门三期工程中，模拟海塘的越浪过程并验证越浪量的准确性。

本文SPH方法在海塘越浪中作了初步的研究，为保证该方法的稳定性和可靠性，仍需进行大量的测试和改进工作。另外，SPH方法的基础理论尚存很多待解决的问题，特别是在固壁粒子与液体粒子相互作用方面，亟需引入一种更为有效的处理方法，以便实现准确而快速的边界压强计算。

参考文献：

[1] Lucy L B. A Approach to the Testing of the Fission Hypothesis. The Astron. J. ， 1977， 8 （12）： 1013- 1024p.

[2]Gingold RA， Monaghan JJ. Smoothed Partide Hydrodynamics： Theory and Applications to Non-spherical Stars. Mon. Not. Roy. Astrou. Soc. ， 1977， 18： 375 -389P.

[3] Monaghan J J. Smoothed particle hydrodynamics[j]. Reports On Progress In Physics. 2005， 68（8）： 1703-1759.