首页/文章/ 详情

POF:用于完全非线性自由表面波传播的物理信息神经网络

3月前浏览2553

上海交通大学的陆昊成等学者提出了完全非线性自由表面物理信息神经网络 (FNFS-PINNs),用于解决非线性潜流中的波传播问题。通过引入准-σ变换和参数无量纲化,FNFS-PINNs 有效地模拟了孤立波、常规波及其逆问题的传播,展示了其在处理复杂水波动力学问题中的高效性和准确性,文章发表在《Physics of Fluids》上。


     

     

摘要:这项研究提出了一种完全非线性自由表面物理信息神经网络 (FNFS-PINNs),这是在PINNs框架内的一项进步,用于解决完全非线性自由表面潜流的波传播问题。利用神经网络的非线性拟合能力,FNFS-PINNs 提供了一种处理非线性自由表面流建模复杂性的方法,拓宽了将PINNs应用于各种波传播场景的范围。改进的准-σ 坐标变换和基本方程的无量纲公式化用于将时间依赖的计算域转换为静止域,并分别对不同维度的变量尺度变化进行对齐。这些创新,以及专门的网络结构和两阶段优化过程,增强了模型在非线性水波数学公式化和可解性方面的表现。FNFS-PINNs通过三个场景进行了评估:特征非线性的孤立波传播、在高色散条件下的规则波传播,以及专注于从晚期状态回溯初始状态的非线性自由表面流逆问题。这些测试表明了FNFS-PINNs在二维垂直场景中计算孤立波和规则波传播的能力。在关注二维波传播的同时,本研究为将FNFS-PINNs扩展到其他自由表面流体问题奠定了基础,并突出了它们在解决逆问题中的潜力。

SIMPOP
1. 引言    

在过去的四十年里,开发易于使用的计算模型来模拟复杂的表面波一直是海洋和沿海工程领域的一个持久挑战。海洋波模型对广泛的自然现象和工程应用至关重要,包括海洋和大气之间的相互作用、波浪的生长和崩溃、沿海洋流的生成以及冲浪区的动力学。此外,波浪的建模和预测对于工程设计特别重要,特别是评估波浪对沿海和海上基础设施影响时。

在水波动力学的研究中,特别是在波浪破碎之前,理论和计算领域最常采用的方法是势流理论。这个框架通过忽略波浪破碎或海洋结构附近的粘性效应来简化对波浪行为的理解,从而在简化模型的同时忽略了波浪行为中固有的一些复杂性。具有自由表面的势流的控制方程是速度势的拉普拉斯方程。为了求解这个方程,势流模型将其与运动和动态自由表面边界条件以及底部边界条件一起整合。在各种势流建模技术中,边界元方法 (BEM) 在建模坡度上的波浪变浅和破碎以及三维波浪在任意海底传播方面尤为突出。最近的进展进一步扩展到了高分辨率自由表面变形和速度场的水波与沉没板相互作用的鲁棒模型。尽管其计算效率很高,BEM 由于完全填充的非对称矩阵而被认为具有数学复杂性,这使得采用高阶数值方案和并行计算策略变得复杂,从而不适用于大规模工程项目。

解决速度势的另一种方法是高阶谱 (HOS) 方法或完全非线性高分散 Boussinesq 方法,其特点是预分析计算域内的拉普拉斯方程解,从而仅需自由表面边界条件的时间积分,大大降低了计算需求。快速傅里叶变换 (FFT) 技术的结合进一步提高了其在大规模波浪建模中的计算效率。然而,拉普拉斯方程解析解的简单性缺乏对水波与沉没或浮动物体相互作用的建模能力。

物理信息神经网络 (PINNs) 因其在解决偏微分方程 (PDEs) 方面的有效性而广受关注。一个通过 TensorFlow(一个领先的机器学习平台)解决非线性 PDEs 正问题和逆问题的突破性框架已经被引入,这种方法标志着科学机器学习的显著进步,通过自动微分将物理原理与数据驱动的见解独特地结合起来。

PINNs 的应用涵盖了广泛的科学领域,展示了它们的适应性和潜力。在热传递领域,PINNs 能够建模复杂的热现象。这种方法被用于改进大气和海洋相互作用的模型,并解决量子物理问题。固体力学研究利用 PINNs 理解材料行为和结构分析。此外,PINNs 在流体动力学和湍流建模中的有效性也得到了突显。总的来说,这些应用强调了 PINNs 在推动科学研究边界方面的广泛实用性和变革潜力,并为各个领域的复杂问题提供了新的解决方案。

PINNs 已证明其在解决复杂水波问题方面的实用性,展示了其在建模自由水面现象方面的灵活性和精确性。为此,早期研究中探索了用于自由边界模拟的不同深度神经网络 (DNNs),并开发了一种根据表面轮廓变化动态更新残差点的算法。这种方法已被有效地用于解决水波逆问题,利用 Serre–Green–Naghdi 方程来估计波面轮廓和深度平均水平速度。一个专门用于自由表面水问题的 PINN 模型被报道,该模型将边界和初始条件直接纳入控制方程。该模型通过其在处理线性和小振幅水波以及晃动问题方面的应用得到了验证。然而,该模型对拉格朗日框架的需求在其应用于更广泛的波传播挑战方面引入了某些限制。同步预测和采样自由表面边界以及处理波传播中的多尺度变化仍然是 PINN 方法在解决完全非线性自由表面波传播问题中的最大挑战。

在 PINNs 框架内,我们引入了一种创新发展——完全非线性自由表面物理信息神经网络 (FNFS-PINNs),旨在解决动态自由表面完全非线性潜流中的波传播挑战。通过利用神经网络的非线性拟合能力,FNFS-PINNs 提供了一种解决非线性自由表面建模复杂性的替代解决方案。这些网络利用 PINNs 在表征物理现象方面的简单性和有效性,促进其在各种复杂波传播场景中的应用。与传统方法如 BEM 相比,FNFS-PINNs 具有若干优势。首先,在训练完成后,FNFS-PINNs 提供了问题的完全连续和可微解,而不是特定物理性质的离散值。其次,FNFS-PINNs 在具有复杂边界的场景中允许更简单的计算设置,例如仅有部分物理数据可用或边界条件包含噪声时,无需额外调整。第三,FNFS-PINNs 可以解决逆问题或时间反转问题,这通常需要在传统方法中进行算法更改,而无需额外设置。与之前的 PINN 方法相比,我们的 FNFS-PINN 更好地解决了同步预测和采样自由表面边界以及波传播中的多尺度变化问题。

我们的 FNFS-PINNs 研究涉及三个精心选择的场景:孤立波传播、复杂边界条件下的规则波传播以及水盆中的时间反转波传播。在这些场景中,仅需定义初始波面条件,FNFS-PINNs 就能独立推断水波速度场,比传统数值技术提供更简化的边界条件方法。本文其余部分的结构如下:第二节概述了完全非线性自由表面物理信息神经网络 (FNFS-PINNs) 的概念和方法框架,包括两阶段训练过程。第三节展示了三个场景的训练结果,从而验证了 FNFS-PINNs 的非线性模拟能力、对建模非线性永久波的适应性以及其在回溯初始波条件方面的能力。结论和未来研究方向在第四节中简要总结,包括本研究的核心贡献并为 FNFS-PINNs 在解决高级波动问题方面的探索指明了方向。

2. 方法    
  • A. 表面波的控制方程

考虑具有非线性自由表面的势流,在笛卡尔坐标系内使用速度势来描述不可压缩的二维流动,如图1所示。自由表面的垂直位移表示为η(x,t),而常数未扰动深度表示为h(x)。流动内的速度场定义如下:

流体域内的压力可以通过以下方程确定:

图1 完全非线性自由表面波传播问题的计算域示意图

边界由Γ表示,包括顶部的自由表面z=η,固体底面z=-h,并包含沿侧面延伸的指定边界条件。主要的控制方程是应用于流体域Ω的二维拉普拉斯方程。利用伯努利方程可以确定该域内任意位置的压力。

相应地,适用于由Γ限定的流体域Ω的连续性方程表现为势流的拉普拉斯方程,

其中下标表示对指定变量的偏微分。在自由表面边界z=η,分别有动态和运动自由表面边界条件,

其中,自由表面的常压通常设为零。动态自由表面边界条件表示自由表面的压力等于大气压力。同时,运动自由表面边界条件确保最初位于自由表面的流体粒子在运动过程中保持在界面上。流体层底部的条件是零法向速度,

在被分类为“干燥”的区域,其中为替代计算高度的指定高度,

其中δ=10^-4是流体层的最小经验厚度。其他边界或初始条件应根据具体问题进行修改。

  • B. 坐标变换与量纲分析

表面波的数学公式化使得在笛卡尔坐标系中流体域的时间依赖性,自由表面被定义为时间的函数。为了克服这一挑战并实现域的静态边界描述,需要进行坐标系的转换。坐标变换成为实现这一目标的常用方法:

此变换代替了z坐标,自由表面对应于0,海底对应于-1。

在机器学习建模中 特别是在实施神经网络时,规范化输入和输出的数据尺度对于最佳性能至关重要。同样,在物理建模领域,不同物理参数在数值范围上往往有很大差异。变量和参数的无量纲化对于协调不同维度变量变化的尺度起着关键作用。考虑一个典型的小振幅余弦波在平底上传播的问题,

借鉴传统问题和坐标变换,我们引入以下坐标变换:

如图2所示,这一变换确保了新坐标系中的波问题变量尺度的一致性,将垂直方向的坐标范围调整为固定范围,底部固定在0,自由表面在α。同时,这种方法有助于波问题无量纲参数的定义,

图2 坐标变换过程示意图

通过这种变换,初始坐标框架(x, z, t)转换为新的框架,有效地将自由表面和底部的动态边界转化为静态边界。同时,在流场动力学中,拉普拉斯方程的表征不仅限于势函数,还包括无量纲势函数和无量纲自由表面高度。

通过这种变换,初始坐标框架转换为新的框架,有效地将自由表面和底部的动态边界转化为静态边界。同时,在流场动力学中,拉普拉斯方程的表征不仅限于势函数,还包括无量纲势函数和无量纲自由表面高度。

采用这种坐标变换需要对数学公式中使用的微分算子进行修改。由于和独立于坐标,通过应用链式法则可以推导出以下关系:

这一关系通过重复应用变换后的微分算子可以获得高阶导数。

显然,坐标变换将时空变化的边界转换为固定不变的边界,便于边界上的采样,并解耦自由表面预测和采样的组合。同时值得注意的是,在坐标变换之后,原本仅与计算域相关的控制方程变得与无量纲势函数和无量纲自由表面高度相关。这样,势函数和自由表面高度可以在整个计算区域内满足拉普拉斯方程。而在之前的PINN方法中使用的笛卡尔坐标系中,自由表面高度和势函数的耦合约束仅限于自由表面边界,这使得优化PINN以获得满意解变得困难。应注意到,坐标变换后,原本仅与域内相关的控制方程被修改为同时依赖和的方程。这种修改促进了自由表面高度和势函数在整个计算域内满足拉普拉斯方程。相比之下,在笛卡尔坐标系框架下,自由表面高度和势函数的耦合约束仅限于自由表面边界。这一差异强调了采用这种准-σ坐标变换的变革性影响,通过将关键物理关系的适用范围扩展到整个域,从而增强了自由表面波传播的建模。

  • C. 物理信息神经网络设置

在建立控制方程、执行坐标变换和量纲分析之后,为波传播问题构建完全非线性自由表面物理信息神经网络 (FNFS-PINNs) 奠定了基础。这些网络通过将完全非线性自由表面波传播问题的控制方程和边界条件直接纳入损失函数项,采用基于优化的方法来解决问题。FNFS-PINNs 的示意图如图3所示。这里,采用全连接神经网络 (FNN) 架构来描述无量纲坐标和相关物理量之间的关系,如下式所示:

神经网络的输入包括变换坐标系中的时间和二维空间变量,旨在通过神经网络函数预测无量纲势函数和无量纲自由表面高度,符号表示网络的可训练参数。值得注意的是独立于,通过掩膜矩阵使所有输入失效,如下式所示:

图3 用于波传播问题的完全非线性自由表面物理信息神经网络 (FNFS-PINNs) 示意图。网络获得输入变量并将作为结果输出变量处理。为了适应网络的特殊性,完全连接的神经网络 (FCNNs) 在对变量应用掩膜后被使用。该方法采用自动微分来计算任何输出相对于输入变量的导数。通过坐标变换和无量纲参数,可以将系统描述恢复到真实物理系统。因此,物理信息部分的损失函数可以通过拉普拉斯方程和其他边界或初始条件来建立。

这里,下标i表示网络的第i个样本点。这些输入随后被送入两个结构相同但不同的完全连接神经网络 (FCNNs)。每个网络由四个隐藏层组成,每层包含256个神经元,并使用双曲正切函数 (tanh) 作为非线性激活函数。每层FCNN结构相同,因此第n层隐藏层与(n-1)层隐藏层之间的关系写为

其中Wn-1是网络可训练参数中的权重,Bn-1是偏置。

组装神经网络后,通过自动微分机制可以计算所有输出相对于输入的导数。这一过程不仅有助于神经网络的学习,还通过利用前述的坐标变换和量纲分析,可以在变换坐标系中从无量纲参数中重构实际物理参数及其导数。

在神经网络的物理信息部分,控制偏微分方程和边界条件的残差在数学公式中起着重要作用。FNFS-PINNs的总损失函数由五个不同的组成部分组成,

其中表示n个样本点的均方误差,代表预测值和目标值之间的差异,这里的目标值为零。

偏微分方程(PDE)损失量化了计算域内拉普拉斯方程的残差。此外,模型还包含了典型边界条件的残差,包括动态和运动自由表面边界条件以及底部边界条件,将这些方程的残差嵌入损失函数中。损失函数的最后

一部分针对具体问题定制,包括初始条件和左边界条件。需要注意的是,所有损失函数的组成部分都经过适当的无量纲化。这确保了所有物理参数的归一化,使其在相同数量级上,促进了混合总损失的收敛。在波传播问题中,物理量的尺度往往存在显著差异。通过对损失函数进行量纲分析,FNFS-PINNs方法可以扩展到不同的尺度,从而提高优化的收敛性。

FNFS-PINNs的优化过程如图4所示。在优化过程中,优化器通过迭代调整网络的权重和偏置来最小化总损失L。首先采用广泛认可的Adam优化器来减少L。选择Adam是因为其在处理大数据集方面的有效性及其自适应学习率能力,这显著增强了训练过程的早期阶段。在KAdam周期(本研究中KA=5000)初步优化Adam之后,采用Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS)算法40,该算法在提高解的精确度方面表现出色。与Adam不同,L-BFGS特别适用于微调阶段,能够更精确地导航损失函数的复杂景观,是在需要精确的情况下实现收敛的不可或缺的工具。这种两阶段优化方法在优化网络性能方面发挥了重要作用,在训练过程中平衡了效率和准确性。

图4 使用两阶段优化方法的FNFS-PINNs优化过程示意图

首先,实施Adam优化算法以增强训练过程的早期阶段,促进坚实的基础学习。随后,为了确保快速和精确的收敛,采用Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) 优化技术。

使用两阶段优化方法的FNFS-PINNs优化过程示意图。首先,实施Adam优化算法以增强训练过程的早期阶段,促进坚实的基础学习。随后,为了确保快速和精确的收敛,采用Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) 优化技术。

为了提高训练结果的稳健性和可靠性,采用了Glorot (Xavier)均匀初始化方法来设置神经网络的初始权重和偏置。选择这一方法是因为其能够保持各层激活的方差,从而防止梯度变得过小或过大,这是确保网络学习过程稳定的关键因素。为了增强训练过程的有效性,采用了伪随机点采样器,以指定的Ks周期(本研究中Ks=1000)对域和边界点进行系统采样,这借鉴了蒙特卡罗模拟的思想,促进了计算域的更均匀和具有代表性的采样。在指定域内计算输出值,标志着模型训练的结束。FNFS-PINN模型在DeepXDE框架内开发,使用PyTorch作为后端以提高计算效率。训练在NVIDIA A100张量核心图形处理单元 (GPU) 上执行,利用其先进的计算能力加速学习过程。

3. 数值结果与讨论    

在本节中,我们展示了在三个精心策划的场景中应用完全非线性自由表面物理信息神经网络 (FNFS-PINNs),如第二节所介绍:孤立波的传播、复杂边界条件下的规则波传播以及水池环境中的时间反演波传播问题。孤立波传播场景被用作常见基准,评估FNFS-PINNs在准确建模非线性波传播方面的能力,从而强调网络在模拟详细波动力学方面的能力。第二个场景进一步探讨了FNFS-PINNs在应对非线性色散规则波传播挑战方面的有效性。最后,第三个场景旨在验证FNFS-PINNs在解决波传播中固有的时间反演问题(涉及从后续时间点的波状态回溯初始状态)方面的效果。

为了便于后续讨论并提供明确的参考框架,我们在表I中列出了这些场景中使用的特征参数。这一列表旨在为每个案例研究的参数定义提供一个综合概览。

表I. 不同验证场景中的特征参数

  • A. 孤立波

孤立波,其特征是对称的抬升自由表面轮廓,具有单个波峰并逐渐向无穷远减小,在均匀水深上传播时保持其形状。该现象代表了一个高度非线性波动力学的典型例子,可以通过拉普拉斯方程与非线性自由表面条件相结合来求解。为了进行分析,我们参考了三阶解作为基准,并采用本文中先前建立的无量纲符号。

将移动坐标引入无量纲形式,并引入:

作为简写符号,第三阶解的无量纲自由表面轮廓和速度分量总结如下:

测试中采用的特征参数列于表I中。在计算域中,我们观察到一个孤立波向右传播,如图5(a)所示。为了获得确定的解,在域的两个侧边界施加水平速度的Neumann边界条件。此外,在初始波峰处的特定点设置条件,并规定初始条件。这需要在方程框架内引入以下约束:

其中,上标“theory”代表对应于理论解的函数。值得注意的是,训练过程并未显式提供任何速度场数据,而 FNFS-PINN 仅基于初始波面条件推导出与波传播相关的速度场。

图5 孤立波传播的自由表面高度,A = 0.1米

(a) 无量纲自由表面高度的理论解:第三阶解。(b) 无量纲自由表面高度的 FNFS-PINN 解。(c) 理论解和 FNFS-PINN 解的偏差。

在计算过程中,域内及其边界上采用了最佳数量的采样点,唯一限制是可用的硬件内存容量。采样点的分布确保了内部采样点与边界点的比例为二比一。附录图16记录了训练损失的演变,最终综合损失收敛成功,表明训练过程有效。附录还提供了训练过程的其他细节供读者参考。

图5和图6(I)展示了孤立波传播的自由表面高度以及相对于理论解的偏差,A = 0.1米。从图5(a)和5(b)中可以看出,二者具有基本一致性。图5(c)中的最大偏差小于0.02,其中 PINN 结果显示在波峰之后自由表面略有下沉,在波峰之前略有升高。为了定量比较两者随时间变化的误差,提出了误差函数随时间的变化,如下所示:

图6. 孤立波在不同时间的自由表面高度及自由表面高度误差的时间历史

(a)–(c) 不同时间的无量纲自由表面高度的理论解和 FNFS-PINN 解。(d)随时间的误差函数。

在图6(I-d)中可以观察到,误差随时间逐渐减小。这种差异可能源于两个潜在原因:首先,Fenton (1972) 的三阶解可能未能完全捕捉满足拉普拉斯方程的完全非线性孤立波的特征;其次,FNFS-PINN 可能仍然存在少量残余未收敛,且选择的离散训练点集未能覆盖整个计算空间。从图6(I-a)–6(I-c)可以明显看出,产生的误差并未显著影响整体波形,波高 \(A\) 的相对误差仅为1.63%。

在图6(II)和6(III)中,我们展示了波幅为A = 0.2和A = 0.3米时的误差。可以观察到,在A = 0.2和 0.3 米时,随时间变化的总体误差与A = 0.1米时的误差量级相同。对于A = 0.2和 0.3 米,波高A的相对误差分别达到1.55%和1.62%。事实上,我们的方法直接求解拉普拉斯方程,理论上问题的非线性不应影响其有效性。本研究的主要目的是介绍这种方法并展示其在建模波传播方面的可行性。详细分析其在高度非线性情况中的性能可能需要在未来工作中进行专门测试。我们选择了非线性较低的场景进行分析,因为我们使用的理论解只是一个近似解,不适用于高度非线性条件。这个选择是为了确保传播后的真实解与理论解之间的任何偏差可以与我们的计算方法的误差区分开来。测试表明,计算误差与这种类型的理论误差紧密相关。

图7展示了孤立波的速度场,由于 PINN 速度场与理论解之间的差异很小,图中仅显示速度场的偏差。为了提高精度并进一步减少与孤立波传播相关的误差,通过调整相关波高从0.1到图中指定时间的波高值来改进理论解。尽管理论解与观察数据之间显然一致,速度场中仍存在细微的差异。发现t = 0时,孤立波的速度场非常准确,但随着时间推移,PINN 速度场与理论解之间的差异逐渐增加。这种差异的原因与自由表面高度的偏差原因相似。值得注意的是,在t = 10时,速度偏差突显了高阶非线性孤立波的成分,进一步支持了 Fenton (1972) 的三阶解不足以描述满足拉普拉斯方程的完全非线性孤立波的观点。相比于水平速度,垂直速度的相对偏差更大,因为其的物理值较小,在拉普拉斯方程中训练到相同残差时,导致相对误差更大。

图7 孤立波的速度场,波幅为A = 0.1米

(a) 理论解的水平速度场。(b) 和 (c) 不同时刻理论解和 PINN 解的水平速度场偏差。(d)–(f) 对应的垂直速度场。

  • B. 规则波

在固定区域内生成非线性规则波通常需要复杂的造波边界条件。采用方程(17)中提出的无量纲移动坐标系,并定义无量纲自由表面轮廓和对应速度分量的三阶Stokes波解如下:

波数和波角频率满足色散关系,

各种理论方法对周期性水波的能力已由Le Méhauté系统地编目。第三阶Stokes波被采用用于建模研究中的非线性规则水波。特征参数列于表I中,由波高H、水深h和波周期确定,

在计算域内,可以观察到一个向右移动的规则波,如图8(a)所示。计算区域涵盖了一个波长和两个波周期的传播,足以进行验证。

图8 规则波传播的自由表面高度

(a) 无量纲自由表面高度的理论解:二阶解。(b) 无量纲自由表面高度的 FNFS-PINN 解。(c) 理论解和 FNFS-PINN 解的偏差。

为了确保定义明确的解,在侧边界施加规则波边界条件。需要注意的是,仅规定了自由表面及其相对于x的梯度,速度剖面由 FNFS-PINN 自动推导。传统的数值造波方法需要给定速度边界条件。FNFS-PINN 与传统的造波方法不同,仅需波面条件即可计算,具有更方便的优势。此外,在初始波峰处指定一个Dirichlet条件,明确规定了和的初始条件,需要在模型中加入以下约束:

这种方法类似于孤立波分析中使用的方法。训练损失的轨迹如附录图16所示,训练以一个综合损失指标结束,表明有效的收敛和模型的可靠性。读者可参考附录了解更多训练过程的细节。

规则波传播期间自由表面高度及其与理论解的差异如图8和图9所示。对比图8(a)和8(b)可以发现,变化很小,图8(c)中最大的差异小于0.003。需要注意的是,图8(c)显示了更复杂的模式,类似于高频波,这表明当规则波传播时,高频误差更容易发生。这些模式表明,PINNs中的误差通常具有高频特征,并且带有详细特征。

图9 不同时间规则波传播的自由表面高度及随时间变化的自由表面高度误差

(a)–(c) 不同时间的无量纲自由表面高度的理论解和 FNFS-PINN 解。(d) 随时间变化的误差函数。

图9允许更好地定量比较这些误差,其中误差函数遵循方程(23)提供的定义。在FNFS-PINN模拟规则波期间,波面差异随时间变化趋于恒定,没有明显的增加或减少模式。模拟开始和结束时的误差最小,这表明规则波边界条件可能限制了随时间推移的误差增长。相对波高误差仅为0.45%。

规则波传播速度场的全面可视化如图10所示,提供了与理论预测的水平和垂直视角比较。图(a)描绘了理论模型的水平速度场。图(b)和(c)显示了初始和结束时FNFS-PINN预测与理论模型之间的细微差异,差异局限于-0.05到0.05的窄色标范围,强调了FNFS-PINN模型的准确性。同样,图(d)显示了垂直速度场的理论剖面,而图(e)和(f)突显了初始和结束时间段PINN预测与理论解之间的微小差异。相对误差大约是横向速度的一个数量级,因为垂直速度相对较小。由于相对波高为0.1,已知速度的绝对误差在同一个数量级。这些图中颜色渐变的一致性和差异尺度的有限范围表明PINN模型与理论预测之间高度一致,确认了模型在模拟规则波传播动态方面的精确性。 

图10 规则波的速度场

(a) 理论解的水平速度场。(b) 和 (c) 不同时刻理论解和 PINN 解的偏差。(d)–(f) 对应的垂直速度场。

  • C. 时间反演问题

为了突显 FNFS-PINN 方法在便捷性方面相对于传统数值技术的优势,我们选择了一个时间反演问题作为示例。该问题的配置如图11所示,涉及一个由两堵垂直墙边界的二维水池计算域。初始时,水完全静止,自由表面定义为

在这个高斯形状的波面中,超过±2.5个标准差的区域被认为几乎是平的。线性波理论提供了一个估计的角频率,利用线性波速和水深的关系。为了增加模型的复杂性,底部边界加入了一个平面斜坡,遵循以下关系,

其中 特征水深由表I指定。这种设置允许在一个简化但有效的计算框架内进行波动力学的深入分析。

图11 时间反演问题的设置

域的侧边界为不透水的壁面。初始条件包括半个高斯波面,流体粒子处于完全静止状态。域的底部边界为一个平坦的斜坡。

图12 不同时间的时间反演问题的自由表面高度及随时间变化的误差

(a)–(e) 不同时刻无量纲自由表面高度的 BEM 解,时间分别为t' = 0, 0.25, 0.5, 0.75, 1。(f) 不同时刻的自由表面高度误差函数。

应用于时间反演问题的边界元法 (BEM) 遵循 Grilli 等人的方法。为了验证 BEM 结果,测试了两种不同的网格分辨率:Delta x = 0.1和Delta x = 0.05。如图12所示,两种网格获得的结果非常一致,确认了 BEM 实现的收敛性。用于与 FNFS-PINN 比较的结果是Delta x = 0.05。BEM 计算的时间步长为0.01秒,总计算时间为2秒。

这个问题的前向演化可以通过传统方法高效解决。在本研究中,我们采用了边界元法 (BEM)。在时间反演问题的分析中,使用了网格尺寸为Delta x = 0.05的 BEM 结果。本研究重点关注时间反演问题,特别是从t = 2的波数据推导t = 0的波条件。在不包括波浪破碎和粘性的情况下,纯粹关注势流,物理系统展示了信息保守性。势流模型假设不可压缩、无粘性和无旋流,动力学完全可以通过满足拉普拉斯方程的势函数来表示。在这种条件下,势流理论描述的流体运动是可逆的,允许系统的任何初始状态通过流体动力学方程精确预测到未来状态。这意味着在我们研究的问题设置中,信息是保守的,使得时间反演问题的解决可行。此外,我们在方程(32)中详细说明了解决时间反演问题的具体边界条件,包括侧壁的边界条件,并保持初始条件不变,从而提供有关波面及其时间导数在t = 2的信息,以及指定点的势函数的常数值,

时间反演问题的无量纲自由表面高度如图13所示,对比了传统 BEM 结果与 FNFS-PINN 方法预测的结果及其偏差。图(a)显示了 BEM 提供的数值解,其中记录了一个大致的振荡运动。在图(b)中,FNFS-PINN 解与 BEM 结果相似,表明波面高度的相似准确表示。图(c)展示了 BEM 与 FNFS-PINN 结果之间的偏差,最大偏差不到0.05。结果突显了 FNFS-PINN 模型的精确性,证明了其能够紧密逼近时间反演问题的基准 BEM 结果。

图13 时间反演问题的自由表面高度

(a) 边界元法的无量纲自由表面高度数值解。(b) FNFS-PINN 的无量纲自由表面高度解。(c) 理论解和 FNFS-PINN 解的偏差。

图14详细展示了时间反演问题在连续时间帧的无量纲自由表面高度及随时间变化的误差演变。该图展示了从t = 0到t = 2的并列图,比较了边界元法 (BEM) 解决方案和 FNFS-PINN 预测的高度在特定时间间隔内的高度,BEM 解决方案以实线表示,FNFS-PINN 预测以虚线表示。误差函数图显示了随时间变化的波面高度差异,使用了方程(23)中的相同定义。该图表明,模拟和预测高度之间的差异很小,表明整个模拟过程中模型的稳定精度。在结束时,波面上的误差最小。随着时间从给定信息点向外推移,波面上的误差逐渐增加。波面的误差函数显示出两个峰和两个谷,误差函数的谷值对应于最高势能点,峰值对应于最低势能点,也是最高动能点。这可能是由于给定条件涉及波面信息,而速度信息中的误差相对较大。总的来说,FNFS-PINN 模型的精确性由狭窄的误差范围突显,展示了其在准确捕捉时间反演问题动态方面的有效性。

图14 不同时间的时间反演问题的自由表面高度及随时间变化的自由表面高度误差

(a)–(e) 不同时刻无量纲自由表面高度的 BEM 解和 FNFS-PINN 解,时间分别为t' = 0, 0.25, 0.5, 0.75, 1。(f) 随时间变化的自由表面高度误差函数。

图15展示了 FNFS-PINN 模型重建的时间反演问题的速度场,描述了选定时间点的水平和垂直速度。图(a)到(e)依次显示了水平速度,而图(f)到(j)展示了相同时间快照的对应垂直速度。这些图清楚地揭示了在和时动能高涨的时刻,其他时间间隔则表现出相对较低的动能水平。这些细节验证了之前讨论的波面误差振荡背后的推论。

图15 时间反演问题的速度场

(a)–(e) 不同时刻的水平速度,时间分别为t' = 0, 0.25, 0.5, 0.75, 1。(f)–(j) 对应的垂直速度场。

4. 结论与未来工作    

在本研究中,我们利用物理信息神经网络 (PINNs) 的框架,创新性地提出了完全非线性自由表面物理信息神经网络 (FNFS-PINNs),以解决完全非线性潜流中的波传播问题。该方法利用神经网络在非线性拟合方面的独特能力,提供了一种替代方案来建模非线性水波传播。此外,FNFS-PINNs 利用了 PINN 在描述物理现象方面的固有便利性,使其能够简便且广泛地应用于各种复杂的波传播问题。

与传统的 PINN 方法不同,本研究引入了准-σ 变换和参数无量纲化的概念。准-σ 变换有效地将时间依赖的计算域转化为静止域,将原始笛卡尔坐标系中仅在自由表面边界上操作的自由表面高度整合到具有自由表面的势流的拉普拉斯方程中。参数无量纲化确保了不同维度变量变化的统一缩放,解决了水波计算中常见的物理尺度差异显著的问题。这为跨尺度的水波传播建模提供了基本框架。设计了专门的输入掩膜网络结构以保证解的准确性和精度。此外,采用两阶段优化过程结合了 Adam 和 L-BFGS 优化器的优势,并通过伪随机点采样器利用了蒙特卡罗模拟的原理。

我们对 FNFS-PINNs 的研究涵盖了三个不同的场景:孤立波传播、常规波传播和波传播的逆问题。孤立波场景作为基准,展示了 FNFS-PINNs 在建模非线性波传播方面的精度。第二个场景评估了 FNFS-PINNs 有效模拟非线性色散波传播的能力。最后一个场景展示了 FNFS-PINNs 解决波传播逆问题的潜力,特别是从后期观察的表面高度重构初始波场。与之前的 PINN 方法相比,FNFS-PINNs 更好地解决了同时预测和采样自由表面以及水波传播中的多尺度问题的挑战。与传统方法如 BEM 相比,FNFS-PINNs 具有显著优势。完成训练后,FNFS-PINNs 提供了一个完全连续且可微的解。FNFS-PINNs 高效地解决了通常需要大量迭代的计算流体动力学中的逆问题。

FNFS-PINNs 方法在二维波传播中的适用性和准确性已得到验证。此外,纳入 FNFS-PINNs 框架的准-σ 坐标变换和无量纲公式化策略有望解决具有自由表面的潜流问题。然而,提升 FNFS-PINNs 的训练速度和增强网络的收敛能力仍然是必要的。幸运的是,FNFS-PINNs 的高度并行化架构表明,通过多 GPU 的分布式并行处理可能显著推进 FNFS-PINNs 的实际应用。尽管本研究侧重于非破碎的二维波,非线性水波传播和破碎的三维场景将在适当时候进行研究。尽管本文未直接应用于反演问题,但对反演问题的考察隐含地突显了 FNFS-PINNs 在具有自由表面的流体动力学领域的潜力。

翻译转载自:《Physics of Fluids》:"Physics-informed neural networks for fully non-linear free surface wave propagation",作者:陆昊成、王千、汤文浩、刘桦



来源:多相流在线
非线性湍流海洋UM理论材料多尺度控制
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-08-14
最近编辑:3月前
积鼎科技
联系我们13162025768
获赞 110粉丝 112文章 306课程 0
点赞
收藏
作者推荐

JFM:CFD模拟地幔对流与板壳运动,新视角解读地球物理板块构造

摘要:众所周知,在对流流体中添加移动边界会产生非同寻常的惊人动态变化,导致从千米级喀斯特地貌到行星级板块构造的壮观地质构造。一方面,运动的固体改变了周围的流场,但另一方面,流动改变了固体的运动和形状。这导致了一种双向耦合,这种耦合对流体-结构相互作用的研究和对地貌的理解具有重要意义。在这项工作中,我们研究了浮板和其下方对流流体之间的耦合。通过数值实验,我们表明这个板块的运动是由下面的流动驱动的。然而,由于板块的存在,流动结构也发生了变化,导致了“热毯”效应,即板块下方滞留的热量导致浮力和上升流,进而将板块推开。通过分析移动边界和流体之间的这种双向耦合,我们能够通过低维随机模型捕捉该板块的动力学行为。从地球物理学的角度来看,热毯效应被认为是大陆漂移的驱动力,因此理解这一机制具有超越流体动力学的意义。SIMPOP1. 绪论 地球内部让几代科学家着迷(Plummer,Carlson & Hammersley 2015)。其中,Leonardo da Vinci(1452-1519)是注意到我们地球不停的地质运动的先驱之一,因为他观察到山脉中存在海洋化石。我们现在知道地球的大陆不会停留在原地,而是经历构造运动,地球地幔中的热对流被认为是这些运动的驱动力(Kious & Tilling 1996)。当流体的温度不均匀导致密度和浮力不均匀时,就会发生热对流,因此温暖的流体上升,而寒冷的流体下沉。这里流体的定义可以非常广泛,因为现代地质学家证实,即使是地幔在大时间尺度上也像流体一样流动(Turcotte & Schubert 2002)。普朗特数(Pr)定义为地幔运动粘度和热扩散率之比,估计约为1023(Meyers等 1987)。地球的核心比表面温暖得多,不稳定的浮力足以驱动地幔对流。作为浮力和粘性效应之间相对强度的衡量标准,地幔中的瑞利数(Ra)约为106(Selley,Cocks & Plimer,2005)。在Rayleigh–bénard对流的充分研究案例中,如此高的Ra已知会导致湍流运动(Ahlers,Grossmann & Lohse 2009)。由于地幔像流体一样对流,其表面流动运送了大陆板块,导致了它们的构造运动。由于地球的空间尺度大,地幔对流的时间尺度长,板块构造的地球物理研究侧重于大陆的当前状态以及预测其重要后果,如地震(Plummer等 2015年)。另一方面,数值模拟(Howard, Malkus & Whitehead 1970; Whitehead 1972, 2023; Whitehead & Behn 2015; Mao, Zhong & Zhang 2019; Mao 2021)和实验室规模的实验(Elder 1967; Zhang & Libchaber 2000; Zhong & Zhang 2005, 2007a,b; Whitehead, Shea & Behn 2011)已被证明是理解板块构造动力学的有效手段。Elder(1967)的早期实验展示了如何在实验室中复现构造运动。在实验室中,石蜡流体层被用来模拟地幔,而漂浮在顶部的塑料薄片被用作模型大陆板块。当石蜡从下面被加热时,对流发生,由于下面对流的剪切作用,板块移动。这样一个简单的实验装置也显示了大尺度的运动;正如Zhang & Libchaber(2000)所观察到的那样,板块在流体表面的两个边界壁之间周期性地移动。通过更详细的研究(Zhong & Zhang 2005,2007a,b),浮板的尺寸对板块运动有很大影响,小板块呈现周期性运动,而大板块则被困在流体表面的中间。将移动热源放置在热对流流体顶部也会产生板块运动,Howard等(1970年)和Whitehead(1972)对此进行了实验研究。尽管这些作品中呈现的几何形状、物理参数和时间尺度与地幔对流非常不同,但它们揭示了令人惊讶的动力学,最重要的是,为大陆漂移背后的流体-结构相互作用机制提供了宝贵的见解。板块构造的数值探索在过去几十年中发展迅速。Gurnis(1988年)提供了第一个与时间相关的大陆漂移数值模拟,其中允许多个大陆合并和分离。在这项工作中,许多其他数值和理论的努力(Zhong & Gurnis 1993; Lowman & Jarvis 1993, 1995, 1999; Lowman & Gable 1999; Zhong et al. 2000; Lowman, King & Gable 2001)采用地幔的地球物理参数,使我们对地球内部的了解迅速推进。同样清楚的是,大陆板块和地幔对流之间的双向耦合导致了不同的动力学(Gurnis 1988; Zhong & Gurnis 1993; Phillips & Bunge 2005; Whitehead & Behn 2015)。最值得注意的是,观察到大板块具有更一致的运动,而小板块则倾向于偶尔运动(Gurnis 1988; Whitehead & Behn 2015)。这些观察结果最近由Mao等(2019年)和Mao(2021)通过解析数值模拟进行了检查。上述工作证实,大陆板块不仅对下方的地幔流平流是被动的,而且还通过热毯效应影响流动结构:众所周知,与大洋地壳相比,大陆地壳的热通量要低得多,因为其地壳深度较大(Mao 2021),因此大陆板块本质上是一个毯子,可以防止热量逃逸并使下方的地幔升温。温暖而轻盈的地幔往往会上升,形成向上的对流气流。随着这种流动向地球表面移动,它会分散并在大陆板块下方产生流体作用力,将板块带走。这是目前对大陆漂移的理解,热毯效应已经在数值和实验上得到验证(Gurnis 1988; Zhong & Zhang 2005, 2007a,b)。本文旨在通过低维模型为板块构造建模提供一个新的角度。在对二维周期域中的板块流相互作用进行直接数值模拟(DNS)后,我们提出了一个板块运动的随机模型,并展示了移动板块如何与其下方的对流流体流进行机械和热耦合。该模型仅保留了热对流的最基本物理性质,恢复了在完整DNS中观察到的动力学,并捕捉了Gurnis(1988)、Whitehead & Behn(2015)、Mao等(2019)和Mao(2021)看到的板块动力学转变。在下文中,我们在第2章总结了方程和数值方法,并在第3章中给出了数值结果。在第4章中系统地推导了随机模型,并在第5章中讨论了它在具有各种纵横比的对流域中的应用。最后,我们在第6章中总结并讨论了数值结果。2. 数值模型 2.1流动方程我们的数值模拟配置如图1所示,其中以x=xp位置为中心的固体板漂浮在对流流体的顶部,该流体在y方向上受到限制,在x方向上具有周期性。在整个研究中,所有长度都由流体深度H重新调整,时间由扩散时间H2/κ(其中κ是热扩散率)重新调整,温度由底部和顶部自由表面之间的温差T重新调整。流体域的x方向是周期性的,周期Γ=D/H(其中D是域宽度),因此整个计算域为x∈(0, Γ)和y∈(0, 1),如图1所示。使用Boussinesq近似,流速u=(u, v)、压力p和温度θ∈[0, 1]的所得偏微分方程为:这里,瑞利数为Ra=αgΔTH3/νκ,普朗特数为Pr=ν/κ,其中ν、α和g分别为运动粘度、流体热膨胀系数和重力加速度。可以针对地球物理地幔对流对流动解算器进行简单修改,但由于我们希望考虑更普遍的流体-结构相互作用情况,并将我们的理论应用于未来的实验室实验,因此我们在本研究中保留了流体和固体板的惯性。图1:移动板块和对流流体示意图流体域以y∈(0, 1)为界,并且在x上是周期性的,宽度为d的浮动板具有中心位置xp2.2边界条件在没有板的情况下,边界条件很简单:流速u=(u, v)在流体/固体边界处无滑移,在空气/流体界面处无剪切。温度θ在底部为1,在空气/流体边界为0。这产生了底面y=0的边界条件:在顶部表面,板有效地屏蔽了热量的逃逸,因此我们在那里取θn=0,同时将流动设置为相对于移动板无滑动。由此产生的边界条件为:或者,这些条件可以强制执行为其中1P是一个指示函数,在板下取值1,否则取值0。2.3板块动力学流体剪切力直接驱动板块运动,所以这里,up=xp(导数))是板块速度,m=ρd是具有线性密度ρ和宽度d的板块的无量纲质量,积分面积P={x|x∈(xp-d/2, xp+d/2)}是板块下的区域。2.4.参数和数值方法附录A中详细说明了用(2.4)、(2.7)和(2.8)求解(2.1)-(2.3)的数值方法。在附录A中,我们使用傅立叶-切比雪夫谱方法来获得精确的数值解。在所有模拟中,我们选择Ra=106,Pr=7.9,Γ=1–16,m=4d(其中d为板宽),与实验中的对流参数相匹配((Zhang & Libchaber 2000; Zhong & Zhang 2005)。图2:漂浮在对流流体顶部的板块运动。(a)各种大小板块的轨迹。覆盖比是Cr=d/Γ,时间t0标志着板块运动的开始。(b)板块的总位移,其中dp通过dp=|xp(导数)|定义。(c)当Cr>0.33时,平均板速度vp=<|xp(导数)|>变高,在Cr=0.56时达到最大值。在这些模拟中,Γ=4,ρ=4,Ra=106,Pr=7.9。3. 数值结果 图2(a)显示了不同大小的板块的几个轨迹,我们可以立即看到板块的大小如何影响其运动。我们定义了覆盖比Cr=d/Γ,以测量宽度为d的板覆盖了多少自由表面,其中对于本节中的所有结果,流体长宽比为Γ=4。对于小板来说,它们的净位移很小,从图2(b)所示的总位移可以更好地看出这一点。如图2(a)中绿色轨迹所示,当Cr大于0.33时,增加板尺寸会出现线性运动。这些轨迹会发生逆转,因为湍流作用力会产生有效噪声。随着Cr进一步增加,线性运动变得更加持久,因为当图2(a)中Cr →1时,板块运动的逆转变得罕见。我们注意到,在地球物理斯托克斯流模拟(Mao 2021)中已经观察到类似的动态行为,因此对于不同的流态,移动板块和下方流之间的耦合机制必须相似。从总位移dp可以看出,在Cr≈0.5时达到最大板速度,这可以通过绘制图2(c)中的时间平均板速度vp来确认。对于小板,平均速度vp仍然较低,但当Cr>0.33时,平均速度vp显著增加,并在Cr=0.5时达到最大值。为了研究动态之间的转换,流体中的典型流量和温度分布如图3所示。在图3(a)中,一个Cr=0.125的小板被放置在对流流体上,并被xm处的下行流体中心吸引(图1),在此处表面流形成一个汇。这种下沉对板块来说是一种稳定的平衡,因为任何偏离这种下沉的情况都会导致作用在板块上的恢复流体力。图3(b)进一步显示了该流汇的结构,其中y平均温度和y平均垂直流速均达到最小值。按照这种表面流动模式,板块位移xp是随机的,如图3(c)所示。由于湍流的随机作用力,板块位置受到噪声的影响,如图3(d)所示,噪声会影响板块速度,其直方图呈高斯分布。板块遇到来自流动的强大“风”的情况很少见,但并非不可能,这种风可以推动板块远离流动汇,穿过流动源,然后再次回到汇,从而导致图3(c)所示的xp跳跃。图3:(a–d)Cr=0.125和(e–h)Cr=0.5的浮板的动力学(a)Cr=0.125的板下流量和温度分布的典型快照(b)对应于(a)的y平均温度和垂直流速,阴影区域表示板块位置(c)板块位移是一个随机过程(d)板块速度是均值为0的随机变量,其概率密度函数是高斯分布,如直方图所示(e)在Cr=0.5时,典型的流量和温度分布显示板由表面流输送(f)对应于(e)的y平均温度和垂直流速(g)板块位移在时间上是线性的,表明是单向平移(h)板块速度具有非零平均值这些模拟的补充动图可在https://doi.org/10.1017/jfm.2023.1071.获得。 覆盖比为0.125的小板浮在对流流体上的动力学。上部面板显示了流场和温度场以及板的位置。下方面板显示了y轴平均温度和垂直流速。 覆盖比为0.5的大板在对流流体的顶部平移。显示了流场和温度场(上图)及其y平均值(下图)以及板的位置。图3(e)显示了Cr=0.5的板的运动情况。在这种情况下,板块运动是单向的,如图3(g)和3(h)所示,向上的速度具有非零平均值。如图3(e)和3(f)所示,移动板倾向于位于流汇和流源之间。随着表面流将板块推向其汇点,流动温度的分布也会发生变化,导致移动的板块追逐移动的表面流汇点。这是热毯效应的直接结果:当板块足够大时,其下方的温度会上升,因为热量无法从那里逃逸。这种局部升温改变了流动温度,并有效地将冷的下行流体推走,导致流动汇位置的移动。总的来说,板朝着冷流散热器移动,同时将散热器推开。因此,对于这个看似复杂的流体-结构相互作用问题,存在一个简单的动力学,我们将从这些观察中推导出一个模型。4. 随机性模型 如图3所示,y平均温度θ的变化强烈影响流型。为了捕捉θ的变化和周期性质,我们用其最低的非平凡傅立叶模式来近似它,其中xm是图1中表面流汇的位置,r=2πΓ-1是波数,常数α衡量温度变化的强度。由该温度分布引起的表面流速U(x, t)=U(x, 1, t)可近似为其中β>0为表面流动强度。事实上,这种表面流动剖面在x=xm处有一个汇点:与xm的小偏差导致x<xm时U>0,而x>xm时U<0,因此在局部,流动指向x=xm。我们注意到,该表面流动剖面与固体/流体边界处的板块速度不匹配,u和up之间的不匹配允许我们估计∂u/∂y(x, 1, t)和由此产生的剪切应力,从而导致板块加速这里,δ是动量边界层厚度(Schlichting & Gersten 2016),由Ra和Pr确定,因此我们在研究中假设其为常数。我们还包括一个标准差为σ的白噪声,代表湍流作用力。总体而言,(4.2)和(4.3)代表了同时具有大尺度环流和小尺度湍流的流动剖面,与高Ra热对流的观测结果一致(Ahlers等 2009;Huang & Zhang 2022)。为了将移动的板块模拟为热毯,我们看一下y平均热量方程,这里,我们忽略了流体平流,q(x, t)=(∂θ/∂y)(x, 1, t)(∂θ/∂y)(x, 0, t)是通过位置x的垂直热通量。假设离开流体-空气界面的热量遵循牛顿冷却定律,并且没有热量穿透板,我们可以将热量方程改写为当x∈P时指示函数1P(x)为1,否则为0,常数γ模拟冷却速率。我们现在从(4.1)中插入θ的值并在x上积分,这导致xm的常微分方程,(4.3)和(4.6)中的积分可以精确计算。定义相位角φ=r(xp−xm),我们得到(up, φ)的闭合动力系统:其中λ=Pr/(ρδ)。一旦知道(up, φ)的动态特性,就可以通过x(导数)p=up和xm=xp-r-1φ计算出(xp, xm)的动态特性。图4:随机模型的模拟轨迹(a)小盘的动态是随机游走的(b)中等大小的板块具有非零的平移速度,其方向可能发生逆转(c)大板块的平移运动更持久(d)所有模拟路径的轨迹如图2(a)该模型中有四个参数:β是表面流动强度,λ=Pr/(ρδ)是阻尼系数,γ是表面冷却速率,σ是随机流体作用力。在物理上,表面冷却速率γ受表面流动强度β的影响,因此我们在该模型中取γ=rβ,这将导致正确的动力学。其余参数可通过数值模拟进行估算。从φ(0)=0和随机值up(0)开始,图4(a)–4(c)显示了不同Cr下xp(t)的一些典型轨迹。在图4(a)中,小板(Cr=0.2)的轨迹是噪声驱动的。对于Cr=0.5的中板,图4(b)显示其轨迹由反向线性平移组成。对于Cr=0.8的大板,图4(c)表明平移是单向的。图4(c)中的板速度与图2中的完全DNS结果相当,并且该速度随着Cr的进一步增加而降低。不同尺寸板的典型位移xp如图4(d)所示,与图2(a)相似,在Cr≈0.3时,噪声驱动运动和线性运动之间存在过渡。因此,这个简单的模型包含了整个数值模拟的所有关键特征。没有噪声,动力系统(4.7)和(4.8)的临界行为可以进一步分析。对于小Cr,很容易看出(4.7)和(4.8)有up=0,φ=2πN(其中N是整数)作为平衡,它们是稳定的,反映了板块运动的被动状态。随着Cr的增加,在Cr*=1/3时出现新的平衡。当Cr>Cr*时,有可能得到非零板速度u*p=(β/π)u^*p,其中平衡φ*可以由下式确定这些状态代表板块的平移。我们注意到(4.9)和(4.10)只是Cr的函数,因此与该模型中假设的所有其他参数无关。为了恢复量纲板速度u*p,只需知道流速系数β/π。U^*p和φ*的可能值如图5所示。对于Cr<Cr*=1/3,u*p=0和φ*=2πN是唯一可能的平衡,反映了始终被表面流汇吸引的小板的被动性质。当Cr>Cr*时,新相位表现为φ*=(2N+1)πarccos([(2Cr)-1-1](cosπCr)-1),这是(4.10)的解,在较大Cr时变得稳定。它们表明较大的板块倾向于位于表面流汇(φ=2Nπ)和源[φ=(2N+1)π]之间,这证实了我们在图3中的观察。当表面流从其源头指向其汇点时(见图5b中的箭头),这些新相位表明了两种可能的板块速度,由(4.9)式给出,如图5(a)所示。此外,图5(a)与图2(c)相似,当Cr→Cr*和Cr→1时,板块速度为零,并在Cr=0.5左右达到最大值。图5:板速度和相位角的平衡值:(a)无量纲板速度up对于小Cr为0,但随着Cr增加而变成平移;(b)相位角φ。蓝线代表板块被流汇吸引的被动状态,红/橙色曲线显示平移状态的平衡阶段。表面流向标有箭头。通过这个简单的模型,我们可以清楚地看到固体板如何与下面的流体相互作用的物理过程,覆盖比Cr可以作为热毯效应强度的衡量标准。对于小Cr,热毯效应较弱,因此板块运动对流体是被动的。当Cr超过Cr*时,热毯效应会强烈到足以改变流量和温度分布,产生上升流,导致板块运动。图2(c)和图5(a)表明,板块速度在Cr≈0.5时达到峰值,并在Cr→1时消失,这也反映了热毯效应和流动对流之间的竞争。随着Cr的增加,更多的自由表面面积被板覆盖,下方的流体力在更大的范围内得到平均。由于该域可能涵盖上升流和下降流,总流体力受Cr影响。这可以在(4.3)式中看到,其中表面流速U提供了板块加速度。在(4.2)式中,U被模拟为正弦曲线,该曲线在(4.3)式中被积分,因此将(4.3)式中的积分面积增加到开放表面的一半(Cr=1/2)将覆盖U的最高贡献。进一步增加覆盖面积将因此减少U的贡献,因为积分域超过正弦函数的半个周期。在极端情况下,Cr=1表示U在整个周期内的积分,导致流体力为零,如图5(a)所示。5. 纵横比效应 到目前为止讨论的所有结果都集中在纵横比Γ=4的对流流体域,这与许多实验研究的几何形状相匹配。随着更复杂的多卷流结构的出现,改变区域纵横比肯定会影响对流流体和移动板块的动态特性(Ahlers等人,2009年)。在本节中,我们研究了Γ=2、4、8和16时的板块动力学,同时保持其他动力学参数与2–4中描述的相同。图6:对流流体在不同纵横比下的快照:(a)Γ=2,(b)Γ=8,(c)Γ=16。在这些模拟中,所有其他参数都固定在Cr=0.2、ρ=4、Ra=106和Pr=7.9。图7:长宽比(a)Γ=2,(b)Γ=8和(c)Γ=16的DNS模拟轨迹。在这些仿真中,ρ=4,Ra=106,Pr=7.9,如图2和3所示。(d)分离被动和平移状态的临界Cr*与Γ成2/3幂律比例图6显示了Γ=2–16时DNS结果的一些典型温度分布。很明显,随着Γ的增加,更多的对流单元出现,出现了更复杂的流动结构。图7显示了不同Γ下的板块轨迹;一个共同的特征再次出现,当大板块平移时,小板块移动很小。为了将我们的简单模型扩展到高展弦比的情况,我们考虑具有更复杂的空间相关性的整体温度和表面流动剖面。其中整数k是傅立叶光谱中最主要的波数。包含波数k是受对流可能具有多卷结构这一事实的启发,因此上面的温度和流量分布图表明流体中存在2k个对流卷,具有k个表面汇和k个表面源。我们跟踪一个表面下沉xm(t),这也是最低流体温度的位置。通过类似于第4章中的推导,我们可以再次获得相位φ=r(xp-xm)和板块速度up的随机动力系统:在没有噪声的情况下,很容易验证(5.3)和(5.4)具有被动平衡u*p=0和φ*=2mπ/k,其中m是任意整数。这种被动状态的雅可比矩阵是在Cr→0的极限中求J的迹和行列式因此J的两个本征值都是负的,表明小Cr的稳定被动态,并证实了我们的数值结果。当J的一个或两个本征值变为正值时,被动状态的稳定性就会丧失,因此det(J)=0提供了一个确定临界Cr*的方程。经过简化后,我们有已知参数β和Γ,(5.8)的根可以用数值确定。在实验室实验和地球物理板块构造中,可以适当估计表面流速尺度β和通气系数Γ,并用于确定Cr*,因此(5.8)提供了通过物理参数确定板块活动性的可能性。我们注意到,对于较大的Γ,Cr*通常较小,如图7所示,在此极限中(5.8)的根可以表示为这里,我们使用了关系式r~Γ-1和k~Γ,后者表明对流辊的数量与长宽比Γ成比例。等式(5.9)确实可以得到验证,如图7(d)所示,从DNS数据测得的临界Cr确实遵循2/3乘幂定律。通过(5.3)和(5.4)式的动力学可以研究更多的内容,未来的实验研究肯定可以用来进一步解决板块和下面的对流流体之间的相互作用。6. 讨论 在这项工作中,我们从数值和理论上探讨了运动平板和其下对流流体之间的机械和热耦合。受到现在和过去作品的启发(Gurnis 1988; Whitehead & Behn 2015; Mao et al. 2019; Mao 2021),我们提出了一个随机模型,表明板块大小(覆盖比)是热毯效应强度的决定因素。对于小板块,屏蔽效应较弱,板块运动对流动结构是被动的;对于大板块来说,下方的流动变得足够温暖,并形成一个上升流中心,将板块推开并导致其平移。提议的简单模型由关于流动及其与板块的机械和热耦合的最小假设组成;然而,它能够预测板块运动的动态转变。虽然直接数值模拟(DNS)是在类似于实验室实验的参数范围内进行的(Zhong & Zhang 2005,2007a,b),但这种随机模型与雷诺数无关,因此可应用于实验室和地球物理场景。实验室规模的实验通常在有界对流系统中进行,因此平板仅限于在两壁之间移动。在有界对流中,流动结构及其与板块的耦合非常不同,因为先前的工作表明,流动具有两个反向旋转的大尺度环流,其强度受板块位置的调节(Zhong & Zhang 2005,2007a,b)。该板块还受到固体边界的限制,因此一旦接触到墙壁,它就会停止移动。在这种情况下,对板块动力学进行建模并非易事,随机变分不等式(Huang et al. 2018)等先进工具可以作为分析这种板块-边界相互作用的数学手段。对流流体的周期域是实验验证我们的随机模型所必需的,这样的实验也将更接近地模拟地幔对流。由环形对流域组成的未来实验可以提供更多细节来验证随机模型。虽然这里我们只研究了单个板块的动力学,但只要能够正确模拟它们之间的相互作用,我们的数值方法就能够处理多个板块。我们目前正在研究这种相互作用,这导致了更加多样化和不可预测的动态。例如,如果两个小板块的Cr<1/3,但它们的总大小达到Cr>1/3,那么我们已经看到每个单独的板块随机移动,但它们的组合“超级大陆”可以平移。在板块构造的地球物理案例中,板块相互作用是许多火山活动和山脉形成的原因,因此了解附近板块的聚散运动可能会为这些地球物理事件背后的流体力学提供新的见解。为了简单起见,这项工作中的模拟和模型都是二维的,将我们的结果扩展到三维是当前的优先事项。使用Chebyshev–Fourier–Fourier方法,我们实现了数值求解器,用于演化位于在两个水平方向上具有周期性的三维域顶部的板。此外,通过Chebyshev–Chebyshev–Fourier方法也可以模拟球形壳上的板块构造,这是一种更接近地球物理板块构造的配置。通过分析那里的DNS结果,我们希望进一步开发我们的随机模型,并将其用于解决地球内部发生的流体-结构相互作用。最后,我们注意到地球物理板块构造比迄今为止进行的任何实验或数值模拟都要复杂得多,因为地球内部的环境如此复杂,现代科学仍在探索之中。但是尽管目前的简单模型不能完全捕捉大陆漂移的动态,我们希望它们仍然可以提供一些关于地球物理学的流体力学见解。翻译转载自:《Journal of Fluid Mechanics》:"Covering convection with a thermal blanket: numerical simulation and stochastic modelling"来源:多相流在线

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈