JFM | 法国学者开发粘性破碎水波数值研究方法，揭示粘度趋零时水波演化过程

1. 绪论    

波浪破碎是自然界中最常见也可能是最美丽的流体流动之一。然而，从建模的角度来看，对它的研究仍然极具挑战性。作为一种强非线性现象，通常的分析方法无法捕捉到其全部机制。此外，由于涉及不同的尺度以及自由表面流动，任何数值方法都会变得非常复杂。我们的目的是提供一个计算粘性水波数值的新框架，允许自由表面翻转，直到界面相交点。

研究波浪破碎现象的另一种方法是考虑具有较小但有限粘度的自由表面纳维-斯托克斯方程。因此可以描述各种不稳定性，包括断裂阶段后形成的充气涡。难点在于如何接近相关的粘度和表面张力消失极限（当包括粘度和表面张力时）。令人惊讶的是，近期包含粘度的研究都没有将其结果与无粘性情况进行比较。 

研究波浪破碎现象的另一种方法是考虑具有较小但有限粘度的自由表面N-S方程。因此可以描述各种不稳定性，包括断裂阶段后形成的充气涡旋丝。难点在于如何接近相关的粘度和表面张力消失极限（当包括粘度和表面张力时）。令人惊讶的是，近期包含粘性的研究都没有将其结果与无粘情况进行比较。

在这里，我们为自由表面的Navier-Stokes流动引入了一种有限元公式，并研究了一个冲入式破浪的情况。我们研究了冲入式破浪上尖锐尖端的形成。我们实现了对欧拉解的收敛，并研究了粘性边界层在界面正则化中的作用。

2. 数学公式    

我们考虑图1中所示的二维水域Ωt，假设它在x方向上是L周期的。在y方向上，Ωt被流体域的刚性底部Γb和水-空气界面Γs,t所包围，后者由时间依赖的参数化曲线γt表示。前者在研究过程中保持不变，而后者则随时间变化。我们使用静止流体的高度h0作为长度单位，√gh0作为速度单位（g为重力加速度），ρgh0作为压力单位（ρ为流体密度，假设为均匀的）。然后，雷诺数定义为Re = ρh0 √gh0/μ，其中μ为流体动力粘度，也假设为均匀的。一些早期的研究使用了深水尺度来定义他们的雷诺数，即Redw = ρgλ3/μ，其中λ是初始波的波长，这里设定为L，因此Redw = (L/h0)3/2Re。

图1：粘性水波问题初始（t=0）域的几何形状

域Ωt中的不可压缩纳维-斯托克斯方程的形式为

在水-空气界面Γs,t上施加无应力边界条件。在底部Γb，我们在切向方向使用无应力条件，并在法向方向施加无穿透条件。上述公式有时被称为“单流体”，因为已经忽略了空气密度。最后，S表示应力张量，定义为S(u) = 1/2 (∇u + (∇u)T)。

由于界面Γs,t是一个物质表面，因此需要解决一个二维对流问题来追踪γt随时间的变化，

使用这种方法，界面的形状不需要是图形的形式。这是描述破碎波的关键特性，并对水波问题的几种表述造成了困难（参见Lannes 2013，第一章）。

我们考虑将有限振幅a的简单波（线性化方程的解）作为初始条件，这样初始界面可以表示为：

其中k = 2π/L表示波数。在界面Γs,t=0上的初始速度u0由无粘性水波问题的一阶二维解的有限振幅扩展给出，

在全流体域中的初始速度u(0, ·)可以通过ka的级数展开来近似。但是，由于我们考虑的是大振幅波，我们更倾向于用数值方法计算它。这可以很容易地通过求解初始速度势φ0在Ωt=0中的拉普拉斯方程来实现（因此假设初始涡量为零）

具有边界条件

为了实现有限元表述，我们需要将方程（2.1）重写为弱形式。这是通过引入函数空间

来完成的，其中H1(Ωt)表示通常的Sobolev空间。然后，我们将速度方程乘以一个任意的测试函数v ∈ H1 Γb (Ωt)，并将不可压缩条件乘以另一个测试函数q ∈ L2(Ωt)。对整个域Ωt进行积分并利用我们的边界条件，得到以下变分表述：找到u ∈ C1([0, T); H1 Γb (Ωt))和p ∈ L∞([0, T); L2(Ωt))，使得对于所有测试函数v ∈ H1 Γb (Ωt)和q ∈ L2(Ωt)，有

3. 数值离散    

我们使用FreeFEM语言（Hecht 2012）对系统（2.8）的解进行数值近似。

一旦对方程（2.8）进行数值积分，就需要根据方程（2.2）来追踪界面的变化。这可以通过任意拉格朗日-欧拉方法（ALE方法）来实现。在每次迭代中，我们数值构造一个网格速度，求解

在这个方法中，下标h表示离散的数值边界。然后，我们使用网格速度w来移动网格顶点。因此，由于网格的移动（参考元素与全局单元之间的变换雅可比矩阵不是常数），这种方法需要在每个时间步重新计算有限元矩阵。w场会从流体方程的对流项中减去。

尽管本研究中忽略了表面张力的影响，但数值方案中很容易包含（2.8）中的毛细管项。扩展到两相流是可能的，但需要同时对两个区域进行网格划分。

尽管在理论上将这种方法扩展到三维情况是可能的，但在实践中将非常具有挑战性。首先，这是因为在三维空间中重新划分网格的问题，但也因为达到如此大的雷诺数收敛所需的自由度数量将变得难以处理，即使使用现代超级计算机也是如此。这是由于在每次时间迭代中都需要在移动流体域上重新计算矩阵。此外，我们方法的一个主要限制在于界面不可能相交。因此，一旦发生飞溅，模拟就必须停止。我们的方法不允许分析破碎后的行为，但请参见例如Iafrati（2009）、Deike等（2015）、Lubin和Glockner（2015）、Deike、Melville和Popinet（2016）以及Lubin等人（2019）。

4. 粘性耗散对界面演化的影响    

我们在数值上考虑了一个初始域Ωt=0，其中L = 2π，h0 = 1，a = 0.5。在后续的计算中，我们忽略了表面张力σ = 0。域Ωt在水-气界面上被离散化为最多4000个顶点。形成大约27万个三角形，因此Navier-Stokes问题约有120万个未知数。

图2：界面随时间演化的示意图，雷诺数Re=10^6。

图2所展示模拟相对应的动画

我们进行了从10^2到10^6的雷诺数范围内的模拟。图2展示了Re = 10^6时的结果，图3则展示了不同雷诺数下界面的比较。所产生的波浪表现为倾覆式破碎波。我们应该强调，由于初始条件是无旋的，除了靠近表面的一个薄边界层外，其他地方的涡量都为零（见第5节）。涡量局部化的特点有助于在雷诺数较大时进行数值解析。

我们现在的目标是表征随着雷诺数增加时的收敛性。图3比较了不同雷诺数下相同初始条件的界面时间演化。为了比较，还包含了相同问题的欧拉方程（Dormy & Lacave 2024）的解。耗散的规则化效果清晰可见。在较大的耗散（即随着雷诺数的减小）下，悬垂区域呈圆形并更快地落下。也许更令人惊讶的是，耗散的影响主要集中在倾覆射流附近。欧拉界面似乎提供了一个极限解，随着雷诺数的增加，纳维-斯托克斯解会收敛到这个极限解。对于Re = 106，欧拉解与纳维-斯托克斯解之间只有很小的差异。这种微小的差异可能是由于Re的有限性，但也可能由于一定程度的数值扩散，因为在这个极端的雷诺数情况下，我们的数值解析能力已接近极限（见下文）。

为了量化有限雷诺数流向欧拉解的收敛性，我们必须测量各种界面位置之间的差异。（给定Re下的数值收敛性是通过改变网格大小来评估的，见附录。）由于界面不是图形，一旦波浪翻转，我们就不能使用标准范数来进行测量。因此，我们依赖于（如在Dormy & Lacave 2024中所做的）曲线之间的双向Hausdorff距离，即

图4(a)展示了给定雷诺数下每条曲线与欧拉解之间距离的时间演化。由于初始条件相同，距离是时间的增函数，直到接近飞溅。随着雷诺数的增加，粘度效应变得显著的时间点也会增加。

图3：对于不同雷诺数值的界面随时间演化的示意图，强调波浪尖端。欧拉解来自Dormy & Lacave（2024）。阴影区域对应于欧拉流体域。

图4：(a) 随着雷诺数的增加，Navier-Stokes解向欧拉解（Dormy & Lacave 2024）收敛的情况，使用Hausdorff距离来衡量；

(b) 对于不同雷诺数值，最大曲率半径随时间演化的示意图。最后四条曲线在这个比例尺下无法区分。

对于这里考虑的初始条件，即使在欧拉解（Dormy & Lacave 2024）的情况下，似乎也没有出现有限时间的楔形奇点。这可以通过引入最小曲率半径Rc(t)来评估（见图4b）。界面的曲率可以通过数值计算得出，最大值对应于波峰。然后，最小曲率半径Rc(t)是最大曲率的倒数。图4(b)展示了Rc随时间的变化。当界面自相交时，每个模拟都会中断。虽然对于足够大的雷诺数，Rc趋向于零，但在我们所有的模拟中，它始终严格为正。对于这个设置，没有获得有限时间的奇点。我们的低雷诺数情况（Re ≤ 10^3）的特征是Rc(t)较大。对于Re ≥ 10^4的曲线在图中无法区分，并且与Dormy & Lacave（2024）的欧拉模拟相吻合，这表明对于此配置没有有限时间的奇点并不是粘度的结果，而且伯努利原理（在波尖附近加速流体）对于此类初始数据并不导致奇点。因此，需要进一步研究其他初始条件和域几何形状，以研究形成此类奇点的必要条件。

5. 能量耗散    

为了进一步描述欧拉解和纳维-斯托克斯解之间的差异，我们现在研究粘性耗散的空间分布。典型的全局能量平衡方程可以通过在弱形式（2.8）中设置v = u并使用不可压缩条件来计算：

对于垂直重力加速度g = -ˆy和平坦底部，这是水波（例如Lannes 2013）的通常动能和势能，并附加了一个粘性耗散项。

动能的局部方程也可以通过将纳维-斯托克斯方程（2.1）来计算，其中u⊥表示逆时针旋转π/2，ω是二维涡量，则有：

因此，粘性耗散发生在涡量不消失的区域。

图5：(a-e) 在时间t = 2.9时，不同雷诺数值下波尖附近的涡量ω。(f) Re = 10^5情况下波尖的局部放大图

涡量在图5(a-e)中进行了表示。注意，随着雷诺数的增加，界面附近形成了一个涡量片。Lundgren & Koumoutsakos (1999) 定理指出，边界层中平均涡量的来源实际上是这种表面涡量片。有趣的是，在波尖处出现了一个小幅度正涡量的边界层（见图5f）。随着时间的推移，当界面的曲率变得重要时，涡量会增长——关于稳态情况，参见Longuet-Higgins (1992)。

图6：沿垂直于边界的直线方向的涡量截面，如图5(a-e)所示

边界层预期会按Re^(−0.5)的比例缩放——例如，参见 Landau & Lifshitz (1987) 的一般理论，Liu & Davie (1977) 关于粘性水波的特定情况，以及 Masmoudi & Rousset (2017) 对此极限的数学描述。我们在图6中展示了从界面 y = 0.925（在图5a-e中标出）开始，沿法线方向的涡量截面。Re = 10^6 情况的边界层最多分布在三到四个节点上，因此指数行为并不明显。然而，图6仍然清楚地说明了Re = 10^3 到 Re = 10^5 的边界层按 Re^(−0.5) 的比例缩放，并且与 Re = 10^6 的缩放比例相兼容。如前所述，Re = 10^2 的情况非常粘稠，因此界面不表现出与其他情况相同的特性。在 y = 0.925 处，向内指向的法向量是上升的，而其他的是下降的。这解释了图6中不同的行为。 

有趣的是，在 Re ≈ 10^4 时，涡量片的大小变得与最小曲率半径 Rc 相当，即当曲率半径作为雷诺数的函数达到最小值时（图4）。图6的另一个显著之处在于边界处的涡量值，它似乎与雷诺数关系不大。这表明了点对点收敛到欧拉问题的界面涡量，即 Baker 等人（1982）中的涡旋强度 γ。

6. 粘性的正则化效应    

界面附近涡量层的形成与边界的粘性正则化有关。实际上，我们可以按照Longuet-Higgins（1953）的方法，定义一个随时间沿界面变化的曲线坐标系（s，n）。这个坐标系有时被称为Frenet框架。界面曲率κ（向内凸时为正值）的时间演化可以表示为：

其中us和un在n = 0处求值。完整的论证可以在Longuet-Higgins（1953，第6节）中找到，其中曲率的符号约定有所不同。

图7：坐标系（s，n）的定义，以及h和κ的几何解释

定义度量因子h（见图7），连续性方程变为：

Longuet-Higgins（1953）中阐述了在上述坐标系中的完整N-S方程。我们引入流函数ψ，使得：

涡量定义为：

在没有涡量的地方，我们可以定义一个速度势φ，使得：

将φ和ψ插入连续性方程（6.2）和涡量方程（6.4）中，得到：

我们现在将流动分解为全局势流和局部粘性分量，粘性效应，局限于大小为δ = Re^(−0.5)的边界层内，在粘度消失时消失（见图5）。我们引入展开式：

此外，由于边界层的结构，我们可以将表达式（6.7）插入涡量方程（6.4）中，会得到O(δ^(−1))阶的项。然而，图5和图6突出显示涡量仍然保持O(1)阶，从而得出结论：ψ = O(δ^2)，因此ψ1 = 0。

用φ和ΨRe重写曲率演化方程（6.1），我们得到：

因此，当表面曲率达到O(δ^(−1))阶（即曲率半径为O(δ)阶）时，粘性耗散对界面的影响就会进入主导阶的平衡。相反，当表面曲率达到O(δ)阶或更小时，粘性效应会在O(δ^(−2))时间内出现。在实践中，对于适用于水波的大雷诺数，表面张力也会在类似的尺度上变得显著。

7. 结论    

本研究表明，即使在波浪破碎的情况下，粘性水波问题也会收敛到无粘解。我们强调了有限粘度的正则化效应，并量化了粘性效应变得显著时的曲率。关于在破碎波尖端可能形成有限时间奇点的问题，需要进一步的研究，涉及不同的初始条件。

附录：数值收敛性

所有模拟都是在不同大小的网格上进行的，网格大小用N来衡量，即自由表面上的点数。（有限元网格的总自由度远大于N。例如，对于Re = 10^6且N = 4000的情况，最终自由度数量为2.2 × 10^6。）根据雷诺数的不同，使用了不同的N值。图8通过使用每个案例中自由表面上具有N/2、N/4和N/8个点的网格，强调了数值收敛性。

图8：每个雷诺数的数值收敛性；左侧指示的矩形区域在右侧放大。请注意，对于Re = 10^6且N = 1000的情况，在大于2.2的时间点上是不稳定的，因此在最后一个图中仅比较了两个网格。