两种水面边界条件下两层流体中内波传播的数值模拟

张周昊，张洪生，王宇鑫，马婷婷

(上海海事大学 海洋科学与工程学院，上海 201306)

海洋内波是最大振幅位于海洋内部的波动，海水稳定层化和外力扰动是内波发生的必要条件。内波作为能量和动量垂向运输的重要载体，有利于海洋生态的保护，但会对海洋工程和水下舰只的安全造成很大的威胁，故深入研究海洋内波的传播特性具有重要意义。研究海洋内波主要有理论研究、物理模型试验、卫星遥感和数值模拟方法。其中数值模拟方法不仅操作相对简便而且经济高效，还可以与理论研究等方法互为对照。随着计算机技术的发展，数值模拟方法的应用越来越广泛[1-2]。

Zhang 等[3]基于Euler 方程，采用改进的SIMPLE 算法，模拟了密度为连续分层流体中内波的传播和变形，并将等水深水域的数值解和理论解进行了对比，定性分析了存在潜堤的变水深水域中不同时刻的密度分布，研究了变水深水域中内波传播的特性。李景远等[4]基于OpenFOAM 平台，引入密度输运方程，基于雷诺时均的N-S 方程，模拟了密度连续分层流体中内孤立波的传播。

虽然密度连续分层更贴近海洋的真实情况，但把实际海洋概化为两层密度不同的流体不失为一种有效的简化手段。高原雪等[5]利用FLUENT 软件，以Euler 方程为控制方程，模拟了密度为两层的水域中内孤立波的生成与传播。陈钰等[6]基于FLUENT 软件，通过给定水位和入口速度的造波方法，数值模拟了弱非线性内孤立波的传播。Terez 等[7]基于二维不可压的Boussinesq 近似的N-S 方程，用重力塌陷的造波方法数值模拟了大波幅内孤立波，通过追踪分布在混合区和密跃层的拉格朗日粒子观察运动。韩鹏等[8]对长非线性内孤立波在斜坡地形上的传播和演化进行了理论和数值模拟研究。Baie 等[9]应用k-ε模型，实现了潜艇生成内波的数值模拟。在利用FLUENT 软件进行计算时，不仅在内波的数值模拟中经常使用紊流模型，而且在表面波的数值模拟中也经常使用紊流模型[10]。

海洋表面边界条件可分为刚盖假定和自由表面假定两种边界条件[11-12]，虽然将刚盖假定视为海面近似边界条件具有一定的合理性，但实际海洋中海面并非如同刚性表面一般静止不动，因此模拟在自由表面边界条件下内波的传播特性具有重要的研究意义。本文利用FLUENT 软件，对两层密度或深度不同的流体，上边界条件分别作刚盖假定和自由表面假定，对微幅内波的传播进行数值模拟，研究在上表面两种假定下，密度场和速度场的异同及其产生原因。

1 两层流体中小振幅内波传播的理论解

小振幅内波问题可归结为二维的线性问题。假定海水中没有其他流动，x轴取在静止的两层流体交界面上，z轴取向上为正，如图1 所示，则不可压缩流体的运动方程和连续性方程可表示为[11]：

图1 稳定层化两层流体中内波传播示意Fig.1 Sketch of propagation of internal waves in density stratified two-layer fluid

式中：u和w分别为流体质点水平和垂向速度；t为时间；p为压强；ρ0=ρ0(z)为静止时分层流体稳定的密度分布，假定其只随深度变化；ρ′为流体的扰动密度；g为重力加速度。

假定上、下层流体交界面处内波频率为 ω、水平波数为k，文献[11]推导了刚盖假定下内波的频散关系式和速度场理论解；文献[12-14]推导了自由表面边界条件下内波的频散关系式和速度场理论解。刚盖假定下的频散关系式为：

式中：ρ1、ρ2分别为上、下层流体密度；b1、b2分别为上、下层流体深度。

自由表面假定下的频散关系式为：

式中：k0=ω2/g；γ=ρ1/ρ2；t1=th(kb1)，t2=th(kb2)。

2 两层流体中内波传播的数值模型

2.1 模型的建立及边界条件的设置

数值水槽如图2 所示，原点取在模型左侧的两层流体交界处。数值水槽中内波传播的有效区域为HJEF，其长度取各工况下7 倍波长L，消波区域为FEKI，消波区长度取为1.5L，设定总水深为100 m。给定振幅A为1.0 m，周期T为400 s。

图2 数值水槽示意Fig.2 Sketch of numerical wave flume

本文采用FLUENT 软件中的k-ε模型进行数值模拟，并利用VOF(Volume Of Fluid)方法追踪自由液面位置，模型的控制方程和VOF 法的介绍可参见文献[6,15]，方程的离散采用了有限体积法。通过改变HI处的边界条件来实现两种边界条件下内波传播的数值模拟。在模拟刚盖假定下内波的传播时，HI处为对称（symmetry）边界条件；模拟自由表面假定下内波的传播时，增加高度为50 m 的空气层，空气层上表面HI处为压力进口（Pressure inlet）边界条件。设HJ、JK和KI均为固壁边界条件。

模型采用结构化网格，传播区x方向网格大小为各工况波长L的2%；在两层流体交界面以及上层流体的上表面及其邻近处，z方向采用加密网格，其余位置z方向网格均取1 m。在消波区域的x方向采用1.05 倍网格大小递增的渐变网格。控制流场全局库朗数（Courant number）为0.02，让FLUENT 自动选择合适的时间步长，总模拟时间为6 000 s。

2.2 数值造波与消波

采用平板拍击的方法来实现造波[16]。造波板上下拍动，带动流体运动以达到造波的目的，其平衡位置在两层密度不同的流体交界面处。根据体积守恒定律，造波板上下拍击时所排开的流体体积，与波面从造波板边缘沿内波传播方向至无穷远处的积分相等，即：

式中：D为造波板长；z(t)为某时刻造波板偏离平衡位置的距离；η(x,t)为两层流体交界面处波面位移函数。

由于内波沿x轴正向传播，因此：

式中：c为内波相速度。将式（5）两侧对时间求导，再结合式（6），设η(x,t)=Acos(kx-ωt)，推导可得：

式中：um(t)为造波板垂向运动速度。由式（7）可见运动函数由需要生成的内波波速、振幅、波数、频率和造波板长度确定。

采用FLUENT 用户自定义函数（UDF）中的刚体运动DEFINE_CG_MOTION 宏，结合推导的运动方程，控制造波板上下平移。

在水槽的末端，即在FEKI区域设置消波层进行消波。消波的方法是在消波区内，在每个迭代步内速度u和w都要乘以消波系数µ(x)。µ(x)有线性和指数等多种形式。现选取线性形式：

式中：a为 常数，本文取为0.02；x0为消波区起始位置水平坐标；x为消波区内任一点的水平坐标；B为消波区长度，常取为1～2 倍波长L。

2.3 敏感性分析

为了验证计算结果的正确性，需对采用不同Δz大小网格的计算结果进行敏感性分析。采用ρ1=1 020 kg/m3、ρ2=1 028 kg/m3、b1=10 m、b2=90 m 的刚盖假定数值水槽作为对比分析算例，该算例比其他算例（下文论及）的数值解与相应理论解之间的差别要大。Δz分别取0.5 和1.0 m。图3 是x=1 200m、z=0 m 处流体质点垂向速度随时间的变化曲线。由图3 可见，两种网格的计算结果差别很小，与理论解的变化趋势一致。鉴于采用Δz=1.0 m时能减少计算时间，因此本文采用Δz=1.0 m的网格进行计算。

图3 (1 200,0)处垂向速度随时间的变化过程线Fig.3 Time series of vertical velocity at point (1 200,0)

3 不同情况下速度场数值解与相应理论解的比较

3.1 上下层流体不同密度差的影响

工况设置如表1所示。假定ρ1、ρ2、b1、b2四个参数及周期T，用两种假定下的色散关系求得各工况下波数k（同种工况两种假定下波数k均取3 位有效数字后近似相等）。由表1 可知，在上下层水深及波周期不变的情况下，上下层流体的密度差越大，波数k就越小，内波传播速度越快。

表1 上下层不同密度工况设置Tab.1 Cases with different densities for upper and lower layers

在x方向传播区长度的一半处，设置上下两个测点。下测点取在两层流体交界面处（z=0 m），监测该处流体质点的垂向速度w；上测点取在上层流体的上表面（z=50m），监测该处流体质点的垂向速度w和水平速度u。计算的水气交界面处水平速度结果见图4～5。由于篇幅所限，只给出Case 1 及Case 4 的部分结果。

图4 Case 1(750,50)位置水平速度历时曲线Fig.4 Time series of horizontal velocity at point (750,50) for Case 1

图5 Case 4(1 750,50)位置水平速度历时曲线Fig.5 Time series of horizontal velocity at point (1 750,50) for Case 4

由图4～5 可知，各工况中的速度在两种边界条件下，数值解与各自的理论解吻合良好；刚盖与自由表面两种边界条件下的数值解基本无差别。对比图4 和5 可知，如其他条件不变，上下层流体之间密度差越大，上层流体表面的水平速度振幅就越大，水平运动越剧烈。鉴于上层流体表面的垂向速度，在刚盖条件下为0，在自由表面条件下尽管不等于0，但幅值很小，因此没有给出其计算结果图。

计算了波峰和波谷时垂向速度和水平速度沿水深的分布，波峰和波谷均取在波型充分稳定的时间段。限于篇幅，文中仅给出了水平速度沿水深分布的计算结果，见图6～7。在两种边界条件下，数值解与理论解基本一致。在波峰时刻，下层流体的数值解和理论解吻合得更好；在波谷时刻，上层流体的数值解和理论解吻合得更好。

图6 自由表面假定下波峰时水平速度沿水深的分布Fig.6 Distribution of horizontal velocity along water depth at the moment for wave crest under the assumption of free surface

图7 刚盖假定下波谷时水平速度沿水深的分布Fig.7 Distribution of horizontal velocity along water depth at the moment for wave trough under the assumption of rigid lid

3.2 上下层流体不同水深比的影响

工况设置如表2 所示。控制上下层流体的密度不变，改变深度比。由于篇幅所限，只给出工况较为极端的Case 5 和Case 8 的计算结果，不同位置处流体质点的速度历时曲线见图8～12。

图8 Case 5(1 200,0)位置垂向速度历时曲线Fig.8 Time series of vertical velocity at point (1 200,0) for Case 5

表2 上下层不同水深工况设置Tab.2 Cases with different water depths of upper and lower layers

由表2 可以明显看出，当上下层流体深度b1、b2互换时，波数k近似相等，这体现了色散方程的“对偶性”[17]。

图9 Case 5(1 200,10)位置垂向速度数值解历时曲线Fig.9 Time series of numerical vertical velocity at point (1 200,10) for Case 5

图10 Case 5(1 200,10)位置水平速度数值解历时曲线Fig.10 Time series of numerical horizontal velocity at point (1 200,10) for Case 5

图11 Case 8(1 200,0)位置垂向速度历时曲线Fig.11 Time series of vertical velocity at point (1 200,0) for Case 8

由图8 和11 可见，两层流体交界面处的垂向速度，无论在哪种边界条件下，数值解与理论解相比均已有较大区别，体现出一定的非线性影响；当b1小于b2时，数值解的波峰处幅值小于波谷处幅值；当b1大于b2时，数值解的波峰处幅值大于波谷处幅值。当b1小于b2时，刚盖假定下所计算的垂向速度幅值大于自由表面假定下所计算的垂向速度幅值；当b1大于b2时则相反。

图12 Case 8(1 200,90)位置水平速度数值解历时曲线Fig.12 Time series of numerical horizontal velocity at point (1 200,90) for Case 8

由图9 可见，当b1很小时，自由表面假定下的数值解有明显的波动。而当b1大于b2时，波动几乎为0，因此文中没有给出其计算结果图。图10 和12 是上层流体表面的水平速度随时间变化曲线。可见，当b1小于b2时，数值解的波峰尖陡、波谷平坦，呈现出非线性影响；而当b1大于b2时，则非线性影响微弱。在入射条件相同的情况下，b1越小上层流体的表面速度越大。

鉴于波谷时刻与波峰时刻的速度沿水深分布图相类似及两种假定下数值计算结果类似，文中仅给出了波峰时的垂向和水平速度沿水深分布，见图13～14，R 代表刚盖，F 代表自由表面。由图13 可见数值解垂向速度的极大值不再位于两层流体交界面处。b1小于b2时，垂向速度极大值的位置向下偏移，导致监测点记录的数据更小；b1大于b2时，极大值的位置向上偏移，监测点记录到的数据也更小。在图8 和11 中也体现出波峰时数值解的幅值相对于理论解有所偏小。图14 是波峰时刻水平速度沿水深的分布。由图14 可见，密跃层越靠近上层表面，上层流体内数值解和理论解差别越大，下层流体内差别越小，反之亦然。同种工况内，两种边界条件下的流体质点水平速度沿水深的分布，没有明显的区别。

图13 波峰时垂向速度沿水深的分布图Fig.13 Distribution of vertical velocity along water depth at the moment for wave crest

图14 波峰时水平速度沿水深的分布图Fig.14 Distribution of horizontal velocity along water depth at the moment for wave crest

各工况在波峰波谷时刻，水平速度与垂向速度的数值解和理论解之间的一致程度可用指标d来衡量[18-19]：

式中：x(j) 为理论解；为 理论解均值；y(j)为数值解。d=0时表示数值解与理论解完全不一致，而d=1表示数值解与理论解完美匹配。

由表3 及图4～7 可见，Case 1～Case 4 尽管下层流体密度不同，但上下两层深度相等，密跃层刚好位于水深的1/2 处，此时总体说来，数值解与理论解相差较小。Case 5～Case 8 上下两层流体密度固定，但上下两层深度差别较大，由表3 及图8～14 可见，此时数值解与理论解的差别较大。特别是当上下两层水深差值较大时，数值解与理论解的差别更大，其中又以密跃层接近静水面时表现得更为明显，这体现了非线性的影响。这是因为数值解是非线性解，而理论解是线性解。在上层水深很小时，在自由表面假定下水气交界面处出现了较为明显的垂向速度。鉴于在实际海洋中上层水深远小于下层水深，由此可以推断，对于实际的海洋分层流体，采用更为真实的自由表面假定应更为合理。

表3 各工况水平和垂向速度的数值解和理论解一致性分析Tab.3 Consistency analysis between the numerical and analytical results for the horizontal and vertical velocities of different cases

4 结语

本文基于k-ε模型和VOF 方法，以FLUENT 软件为计算平台，采用平板造波的方法，实现了刚盖假定和自由表面假定下小振幅内波传播的数值模拟，通过比较z方向不同网格大小的计算结果，检验了计算结果的合理性。

设置上下两层流体不同的密度差和深度，数值模拟了刚盖假定和自由表面假定下多种情况内波的传播；通过与各自理论解的比较，阐述了数值解与理论解之间的异同；通过对两种假定下计算结果的比较，分析了两种假定下数值计算结果之间的差别及其产生的原因。在上下层水深相同的前提下，尽管上下层流体密度差不同，但此时总体上数值解与理论解相吻合；尽管在自由表面假定下，静水面处存在垂向速度，但其数值很小。在上下两层流体密度差固定的前提下，上下两层流体深度差值的改变会较为明显地影响数值计算结果，以密跃层接近海洋上表面时表现得更为明显，此时数值解与理论解差别较大，这体现了非线性的影响；在自由表面假定下出现了较为明显的垂向速度，这体现了海洋上表面边界条件的影响；在给定入射条件下，上层流体水深越小，上表面处速度越大。

鉴于在实际海洋中，一般密跃层以上水深远小于下层水深。当计算运动幅值更大的内孤立波时，可以设想此时边界条件的影响会更加明显，这便是本研究意义所在，也是下一步将要研究的内容。本文工作也为内波与表面波之间相互作用的研究提供了参考。