推板造波的光滑粒子流体动力学数值模拟探讨

1、推板造波的光滑粒子流体动力学数值模拟探讨    1 前言造波问题是船舶工程和海洋工程领域十分普遍的现象，是模拟船舶与海洋结构物相互作用波浪环境的基础。过去的几十年中，造波问题的数值方法不断得到更新和发展。最早对造波问题进行研究是基于势流理论的边界元方法(BEM)。势流理论假设流体没有旋度，没有粘性，不可压缩，流体的速度势 满足拉普拉斯方程。Longuet-Higgins 和Cokelet最早采用边界元方法对二维水波问题进行研究。Baudic et al对于非线性水波的产生，传播和吸收都做了详细研究。贺五州和段文洋采用完全非线性边界元方法对二维摇板造波问题

2、进行了时域计算。该类方法的优点可以自然满足质量守恒，不存在能量耗散，计算量较小，但是对粘性效应，自由表面发生翻转、破碎现象还无法模拟。随着计算机运算能力的不断提高，很多学者发展和应用基于直接求解N-S 或者雷诺平均N-S 方程的方法。该类方法最具挑战性的问题是如何处理非稳态非线性自由表面。目前，常用的方法是VOF 方法和Level Set方法，该类方法都是通过求解每个网格上的一个标量函数来描述界面的几何特征，标量函数的输运方程形式上完全相同，同属于界面扑获类算法，都有较强的拓扑表达能力。但是目前该类数值方法对于自由表面，特别是伴随有翻卷、破碎等现象的复杂自由表面问题，所取得的成就还很有限。目前

3、，一种被称为无网格的方法引起人们关注。因为在处理大变形的情况下，不需要重新划分网格。该类方法对于复杂自由表面问题具有较大的潜力。本文所采用的光滑粒子流体动力学(SPH)方法是无网格方法的一种，该方法的计算建立在一群离散的粒子上，函数以及导数的求解可以利用周围粒子的相关值积分得到。某个粒子影响区域的大小称为支持域。支持域内越靠近中心粒子影响越明显，支持域外粒子对计算结果没有影响。SPH 模拟方法最早用于天体物理学，并由Monaghan将该方法引入自由表面流动问题。由于SPH 方法中是对每个粒子的运动轨迹进行模拟，这些粒子能够自动满足自由表面的运动学条件。目前，SPH 方法在自由表面流动问题上已经

4、取得一些成就，例如Y.M.Lo& Shao用SPH 方法结合大涡模拟方法对近岸孤立波的爬升进行模拟。Colagrossi &Landrini采用SPH 方法对二维两相流问题进行计算分析。Souto et al采用SPH 方法对液舱晃荡的相位滞后问题进行研究，之后又实现了SPH 方法对液舱晃荡横摇力矩的模拟。国内对SPH 方法的研究已经开始，但是还处于学习和消化阶段。虽然SPH 方法具有巨大潜力，但是该方法同时还存在一些问题有待解决，例如计算效率较低，数值结果不很稳定，但是该方法编写程序较容易，对于大变形问题和动边界问题的处理方便。因此，本文尝试采用SPH方法对二维推板造波问题进

5、行模拟，通过与全非线性边界元和VOF 方法比较，对SPH 方法模拟造波问题的适用性进行研究。本文的主要内容包括四个方面。首先对SPH 方法的基本原理进行简单介绍，给出SPH方法模拟自由表面流动问题的控制方程，然后对求解控制方法的若干技巧进行了介绍，包括压力与密度的显式关系，人工粘性，边界条件处理方法。然后给出二维推板造波问题计算结果，通过与全非线性边界元和VOF 方法比较，得到了SPH 方法模拟该问题的精度特点，特别是当波面出现翻转现象时，SPH 方法都能够准确的反映出来。最后对SPH 方法计算造波问题的精度进行了总结，对存在的不足提出改进方法。2 基本原理2.1SPH 方法SPH 方法的基本

6、理论包括两部分：核近似和粒子近似。核近似的思想来自函数f (x)的积分表达式(1)式中为包含x的积分域。h 是核函数的光滑长度，由于核函数W 不是 函数，故(3)式只是近似式，(3)式这种表示函数的近似方法称为核近似(Kernel approximation)。核函数有多种形式，其中三次样条核函数，高斯核函数，五次样条核函数应用得最广泛。所有核函数需要满足一些特性，如正交性，紧支性，非负性，单调性等等。求解流体动力学方程的关键除了求函数本身，更重要的是如何用核近似来计算函数导数。根据流体系统是由具有独立质量，占有独立空间的有限个粒子组成的假设，可以将连续形式的函数核近似用粒子近似表达。2.2

7、SPH 形式的流体动力学方程水波问需要求解的方程包括质量守恒方程和动量守恒方程。其中i F 表示质量力，在水波问题求解时i F 表示重力项， 表示某一矢量方向。接下来对求解上述方程的一些技巧进行介绍。3 数值解法若干技巧3.1 人工可压缩性在SPH 方法中认为流体系统具有弱可压缩性。流场密度和压力的关系可以通过显式关系直接表示出来，这种压力的确定方法就称为人工可压缩性。P = C ，对于水 = 7， 0 C 为声速，可取流体区域中粒子最大速度的10 倍。这是目前SPH 方法模拟水波问题所普遍采用的形式。3.2 人工粘性在SPH 方法的计算中，由于粒子容易发生穿透，重叠现象，使得核近似计算结果不

8、稳定，为了避免以上现象发生，Monaghan在动量方程中加入人工粘性项，主要作用可以避免粒子相互穿透，使得整个区域的粒子排列更加有序，可以较好的改善计算稳定性。3.3 边界处理当粒子位于计算区域内部时称为自由粒子。当粒子的支持域没有被边界截断， (3)式和(5)式成立。当粒子位于边界区域附近时，粒子支持域被边界截断，此时 (3)式和(5)式不再成立，边界附近处就要采用特殊的方法进行处理。Monaghan3提出过固壁边界的处理办法。为了保证固壁边界不可穿透，在边界处安放一些固壁边界粒子，当自由粒子靠近固壁粒子时，固壁粒子将施加外力使自由粒子远离固壁边界来保证自由粒子不穿透固壁。但该方法一方面缺乏

9、物理意义，另一方面当自由粒子速度较高时，固壁边界往往会施加过大的外力而使整个粒子系统不稳定。本文采用一种镜像粒子方法。在边界外侧生成一些镜像粒子，镜像粒子与对应自由粒子到边界的距离相同，镜像粒子与所对应的自由粒子的速度切向相同，法向相反，如图1 所示。这样可以保证边界处的自由粒子没有密度缺损，而且可以保证自由粒子和镜像粒子的速度累加效应在固壁边界法向上为零，满足固壁边界不可穿透条件。3.4 步进方法在解问题(8)(11)时，通常采用的步进方法主要有三种，改进Euler 法，蛙跳 (leap-frog)法和四阶的龙格-库塔法。但是通过大量的计算发现，这些显式的步进格式对计算结果影响并不明显，因此

10、本文采用了较简单的改进Euler 法，其取值与粒子数和声速有关，一般粒子数目增大或声速增大，时间步长t取值需减小。4 数值结果与分析推板造波问题是研究波浪与结构物相互作用的典型例子。目前对于该方面的研究通常还只限于小波幅非线性波的模拟，但是对于大幅强非线性波模拟，特别是伴随有翻转、破碎现象的造波问题，传统的数值方法还面临极大的挑战。本文所计算的推板造波模型如图2 所示，在平坦池底的左侧安放有造波机，右侧为一垂直固壁。由于本文主要模拟SPH 方法短期造波过程，假设波浪没有到达右侧墙壁，不存在波浪反射问题。有关避免波浪反射，增加阻尼区的研究正在进行中。坐标原点建立在初始时刻造波机与静止水面的交点，

11、x 轴水平向右为正， y 轴垂直向上为正。推板运动速度u a t w = sin ，位移s a t a w = ? cos + ，a为推板运动幅值， 表示推板运动频率。水深hd 1.0 米，水池长度= 20 f l 米。4.1 弱非线性波的模拟当取推板运动频率 = 3.14秒-1，运动周期T = 2 = 2秒，推板运动幅值a = 0.1米，由有限水深的线性波色散关系可以估算所造波数k = 2 =1.01米-1，波长 = 6.25米，并由传递函数21 可得生成的波高4 sinh2 /(2 sinh 2 ) 0.0990 = + = d d d a a kh kh kh 米， 波陡0.032 1

12、20 0 a = < ， = 0.16 < 1 2 d h 。可以判断生成波属于有限水深弱非线性波。对于该问题的计算，主要是验证SPH 方法模拟造波的精度特点。整个流体区域分布400*20=8000个粒子，每个粒子空间大小V = 2.5e ? 4 i 米2，每个粒子光滑长度= 0.05 i h 米，核函数采用三次样条核函数，人工粘性系数= 0.06 ，时间步长t = 2.0e ? 4秒，步进50000 步，在Core23.0 的PC 机上计算时间5.1 小时。图3 给出SPH 方法在不同时刻推板附近区域的粒子分布结果。由图中的结果可以看出，SPH 方法计算的自由表面非常光顺，在推板

13、和底部边界附近没有出现粒子发散和穿透现象，表明镜像粒子的边界处理办法是稳定可靠的。其中i f 称为体积分数，表示i 粒子影响域内粒子空间与影响域面积的比值， i f 越大说明该粒子周围的粒子数越多， i f 越小说明该粒子周围的粒子数越少。j V 表示i 粒子支持域内第j 个粒子占据的空间， i h 表示第i 个粒子的光滑长度， 表示与核函数有关的光滑长度伸展系数，当核函数取三次样条函数时 = 2.0。通过大量的计算表明，当体积分数 0.75 i f 就可以认为，该粒子位于自由表面上。图4 给出了BEM 方法与SPH 方法结果的比较。由于自由表面某些区域的粒子存在相互挤压，根据以上判断标准对少

14、数区域不能准确分辨存在空缺，但是从总体效果来看，该方法还是能够反映出整个波面形状的发展趋势。由自由表面形状的比较结果可以发现，在模拟线性规则波时SPH 与BEM 方法的计算结果吻合得较好，体现了SPH 方法在弱非线性自由表面波的模拟上具有较高的精度。4.2 强非线性波的模拟当增加推板运动频率 = 6.28秒-1，推板运动幅值保持a = 0.1米，运动周期T =1秒，由有限水深的线性波色散关系可以估算所造波长 = 1.56米，波数k = 4.02米-1，同样由传递函数可得生成波幅0.199 0 a = 米，波陡2 0.255 1 7 0 a = > ， = 0.64 > 1 2 d

15、h 。可以初步估计生成波属于深水强非线性波，生成的波浪会发生翻转破碎现象。对于该问题的计算，主要是验证SPH 方法在模拟强非线性波破碎波的精度特点。考虑到模拟波长减小，将计算区域长度取为=10 f l 米，既保证在模拟时间内波浪不会发生反射，又可以减少粒子数目，减少计算量。当整个流体区域分布200*20=4000 个粒子，其它因素保持不变的情况下，步进25000 步计算耗时2.5 小时。图5 给出SPH 方法从初始状态到3.0 秒时，每隔0.5 秒推板附近区域粒子分布结果。由图中结果可以看出，SPH 方法计算出的波高度大约为0.4 米，波长在1.5 米至1.8 米之间。波面形状并不光顺，但是自

16、由表面发生翻转破碎的现象还不清晰。可以看出SPH 方法模拟的最大波高处往往不连续，发生了比较模糊的破碎现象，而BEM 方法自由表面始终保持光顺没有发生破碎现象。在推板附近波浪形成阶段，两种方法计算的自由表面形状符合得较好，当波浪传播大约一个波长后，SPH 模拟波面与BEM 计算结果存在相位差，而且随着波浪传播距离增大相位差也增大。这是由于SPH 方法模拟这种较陡峭的波峰会发生破碎现象，使得生成的波长增大波速增加的结果。为了能够更加清晰的观察到翻转破碎现象，将粒子总数增加到400*40=16000，在1.5 秒至1.8 秒时刻对推板附近区域进行放大。当粒子数增加后在1.5 秒时，SPH 方法计算

17、的波形开始发生翻转，在1.6 秒和1.7 秒可以清晰的观察到波面翻转的发展过程，在1.8 秒时翻转现象结束。BEM 方法在所对应时刻还无法对该现象进行清晰的捕捉。因此，为了更好的对于波浪细节进行描述，必须使该区域内的粒子数目增加到一定范围。本文推荐能够对波浪翻转进行清晰描述的粒子半径/8 0 r a i ，其中0 a 为浪高。为了对该频率和幅值下波浪发生的翻转现象进行更加充分验证，本文还采用FLUENT软件对该问题进行了计算。FLUENT 软件对微分方程的求解采用有限体积法，自由表面捕捉采用VOF 方法。左侧推板采用动网格方法，推板运动方式编写UDF 程序加以控制，运动方式与SPH 方法完全相

18、同。图8 给出了不同时刻两种方法的自由表面的比较结果。其中SPH方法的粒子数为16,000，VOF 方法的网格数为200,000。图中可以看出在各时刻，SPH 与VOF 方法的计算结果总体上都吻合的很好，波面形状保持了较好的一致性。只是在自由表面发生翻转破碎后，复杂的自由表面形状还有待实验的验证，例如射流发生的时间和形状，气穴生成的大小和位置等等。5 结论与展望本文采用SPH 方法对推板造波问题进行模拟研究。当所造波是线性波时，由于BEM 方法已经能够对该问题进行高效、准确的计算，因此SPH 与该方法完全吻合的计算结果，可以充分证明SPH 在模拟线性波时具有较高的精度。当模拟非线性波，特别是破碎波的发展过程时，BEM 方法就显得力不从心，但是SPH 却能非常方便的实现波浪翻转，破碎，重入等复杂的自由表面现象。通过与VOF 方法的比较，两种方法所计算出的波形在不同时刻都具有较好的一致性，只是波浪破碎后所产生的射流，