一种基于流固耦合框架的动力地形热力学模拟方法

本发明涉及数据模拟，特别涉及一种基于流固耦合框架的动力地形热力学模拟方法。

背景技术：

1、在地球动力学、地质工程和地震预测等研究领域，准确模拟地壳和地幔在数百万年时间尺度上的长期变形、热对流及断层活动，对于理解板块运动、造山作用、盆地演化和地震周期性活动等重大地质过程具有基础性科学意义。地球内部物质表现出复杂的流固耦合与热力学耦合特性：地幔对流本质上是高黏度流体在温度梯度驱动下的热对流过程，而地壳和岩石圈则在长期构造应力作用下表现出弹性或弹塑性固体行为，两者的相互作用——即地幔对流驱动岩石圈变形、岩石圈变形又反过来影响地幔对流格局和热输运过程——是驱动板块构造运动、造山运动、地震活动等核心地质现象的根本物理机制。因此，建立能够统一描述上述流固耦合、热对流和复杂材料流变行为的数值模拟框架，是计算地球动力学领域长期追求的核心目标。

2、在现有的数值模拟方法中，早期的大地构造模拟代码，如fullsack于1995年发展的sopale程序，采用了任意拉格朗日-欧拉方法，在“蠕变流”极限下求解斯托克斯方程，实现了对黏塑性材料大变形行为的模拟。然而，此类方法通常将地壳和地幔统一视为具有不同高黏度的流体介质，难以准确刻画流固界面处的强耦合效应，特别是地壳与地幔之间的莫霍面，也无法完整描述固体岩石圈的弹性或弹塑性有限变形行为。近年来，以wick于2013年在deal.ii有限元库中实现的单相流固耦合求解器为代表，基于任意拉格朗日-欧拉框架的流固耦合数值方法取得了重要进展。该方法采用双调和方程构造流体网格运动以处理大变形，并通过牛顿法求解强耦合的非线性系统，在生物医学和工程领域的流固耦合问题中表现出良好的适用性。然而，该求解框架的设计目标主要针对生物医学和工程应用，缺乏针对地球动力学问题中长时间尺度、大规模三维模拟需求的专门优化，且未集成热力学方程，无法描述温度场与流固耦合过程之间的相互反馈。

3、在热力学模拟方面，以deal.ii有限元库的step-32教程程序为代表，展示了如何基于boussinesq近似大规模并行求解描述地幔对流的热对流方程。该程序能够有效模拟地幔热对流的温度场和速度场演化，具备良好的并行可扩展性，但未与固体变形模块耦合，将岩石圈简化为对流系统的上边界条件而非可变形固体，因而无法完整描述地幔对流驱动岩石圈变形、岩石圈变形影响地幔对流这一闭环耦合过程。

4、综合来看，现有技术的主要不足在于：一方面，现有的流固耦合求解框架缺乏与热力学模块的深度集成，无法建立热-力-流完全耦合的统一数值模型；另一方面，现有的地幔热对流求解器未与固体力学模块耦合，难以模拟热膨胀、热应力以及变形对热输运的反馈等关键物理过程。上述不足导致当前缺乏一种能够在一个统一的任意拉格朗日-欧拉变分框架下，将流体力学方程、固体力学方程和热力学方程进行单块耦合求解，难以高精度、高效率地模拟地幔对流驱动岩石圈变形这一完整物理过程的动力地形热力学模拟方法。

技术实现思路

1、本发明实施例的目的是提供一种基于流固耦合框架的动力地形热力学模拟方法，通过建立热-力-流多物理场在统一任意拉格朗日-欧拉框架下的单块耦合变分模型并进行联合求解，提高了地幔对流驱动岩石圈变形模拟的物理保真度，解决了现有技术无法在同一数值框架中准确刻画流体运动、固体变形与热输运三者之间闭环耦合反馈的技术问题。

2、为解决上述技术问题，本发明实施例的第一方面提供了一种基于流固耦合框架的动力地形热力学模拟方法，包括如下步骤：

3、步骤s1：加载目标地质区域的初始温度场，在每一个时间步起始阶段，利用前两个已知时刻的温度场执行温度外推，得到当前时间步的预测温度场；所述温度外推为以已知时刻温度场为输入、通过时间外插构造当前时刻温度场近似值的操作；

4、步骤s2：以所述预测温度场为热力学初值，在同一牛顿迭代框架下依次执行求解流体域动量守恒与质量守恒方程、更新黏度参数、求解固体域有限变形方程和求解全域热力学方程，并在流固交界面施加速度连续条件，构成对所有物理场变量进行联合求解的非线性块耦合迭代；

5、步骤s3：判断当前迭代的残差是否收敛，若未收敛则调整牛顿迭代步长后返回步骤s2继续迭代，若收敛则依据收敛解更新网格与场量，自适应确定下一时间步的步长，执行速度外推后进入下一时间步循环，直至模拟结束，输出动力地形演化结果。

6、进一步地，步骤s2中所述对所有物理场变量进行联合求解包括：

7、s21：识别所述非线性块耦合迭代中代数方程组的块稀疏结构，按物理场类型将未知量划分为第一变量组和第二变量组；

8、s22：针对划分后的变量组构造schur补矩阵，形成以所述第一变量组为未知量的缩减系统；

9、s23：求解所述缩减系统，得到第一变量组的修正量；

10、s24：将所述第一变量组的修正量回代，计算第二变量组的修正量；

11、s25：以所述第一变量组的修正量和所述第二变量组的修正量联合更新所有物理场变量，完成本次联合求解。

12、进一步地，步骤s21中所述识别所述非线性块耦合迭代中代数方程组的块稀疏结构，包括：

13、s211：将所述代数方程组按流体速度、压力、固体位移和温度四个物理场变量进行分块排列，得到块稀疏矩阵；

14、s212：识别所述块稀疏矩阵中流体速度子矩阵的对称正定性和固体位移子矩阵的块对角性，所述对称正定性用于指导几何多重网格预处理器的选取，所述块对角性用于指导代数多重网格预处理器的选取。

15、进一步地，步骤s21中所述按物理场类型划分的方式为：以流固界面位移和温度作为第一变量组，以流体速度、压力和除界面位移外的固体位移作为第二变量组，使得控制流固耦合与热力耦合的关键物理量归入第一变量组，可由第一变量组通过耦合关系驱动的从属物理量归入第二变量组。

16、进一步地，步骤s221中对所述第二变量组的子矩阵构造其近似逆采用混合多重网格预处理：

17、s221a：对流体速度与压力对应的子矩阵块施加几何多重网格预处理，得到该子矩阵块的近似逆；

18、s221b：对固体位移对应的子矩阵块施加代数多重网格预处理，得到该子矩阵块的近似逆；

19、s221c：将两步得到的近似逆按分块对角结构组装，形成所述第二变量组子矩阵的整体近似逆。

20、进一步地，步骤s3中所述自适应确定下一时间步的步长，包括：

21、s31：基于数值稳定性条件计算候选步长；

22、s32：判断当前模拟时间所处的数量级区间，根据所述数量级区间对所述候选步长施加对数时间尺度调节，得到最终采用的下一时间步步长；

23、s33：以所述最终采用的下一时间步步长更新当前时间步长，作为下一时间步循环的时间推进步长。

24、进一步地，步骤s32中所述判断当前模拟时间所处的数量级区间，根据所述数量级区间对所述候选步长施加对数时间尺度调节，得到最终采用的下一时间步步长，包括：

25、s321：将模拟总时长划分为多个数量级区间；

26、s322：判断当前模拟时间所处的数量级区间，当模拟进程从前一数量级区间进入后一数量级区间时，按预设的区间步长映射关系调整所述候选步长，使得调整后的步长的对数与模拟进程呈线性关系，得到所述最终采用的下一时间步步长。

27、进一步地，步骤s2中所述流体域动量守恒与质量守恒方程采用boussinesq近似的不可压缩navier-stokes方程，所述固体域有限变形方程采用以piola-kirchhoff应力张量描述的超弹性或黏弹塑性本构方程，且步骤s3中所述更新网格与场量中的网格运动由双调和方程控制，所述双调和方程为以拉普拉斯算子的复合作为控制算子的网格运动方程。

28、进一步地，所述流体域动量守恒与质量守恒方程、固体域有限变形方程和全域热力学方程在参考构型下构造为统一变分弱形式，所述统一变分弱形式以流体速度、压力、固体位移和温度为未知量，包含流体动量项、流体对流项、流体应力项、压力散度项、固体惯性项、固体应力项、热瞬态与对流项、热扩散项以及浮力项、体积力项和内热源项，所述流固交界面处的速度连续与应力平衡条件满足在所述统一变分弱形式。

29、进一步地，在步骤s1至s3的时间推进过程中施加多机制数值稳定化处理，所述多机制数值稳定化处理包括：步骤s2中以向后二阶差分格式对时间导数项进行离散、步骤s2中在求解流体域方程时进行压力正则化处理、以及步骤s2中在求解热力学方程时施加基于熵黏度的人工黏性处理。

30、相应地，本发明实施例的第二方面提供了一种电子设备，包括：至少一个处理器；以及与所述至少一个处理器连接的存储器；其中，所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被所述至少一个处理器执行，以使所述至少一个处理器执行上述基于流固耦合框架的动力地形热力学模拟方法。

31、相应地，本发明实施例的第三方面提供了一种计算机可读存储介质，其上存储有计算机指令，该指令被处理器执行时实现上述基于流固耦合框架的动力地形热力学模拟方法。

32、本发明实施例的上述技术方案具有如下有益的技术效果：

33、1. 通过识别多物理场耦合代数方程组的块稀疏结构，按物理场类型将未知量划分为包含流固界面位移和温度的第一变量组与包含流体速度、压力和其余固体位移的第二变量组，并针对不同物理场对应的子矩阵块分别施加几何多重网格与代数多重网格预处理，显著加速了大规模非线性系统的收敛速度，实现了比通用求解器更高的计算效率，解决了现有技术在处理流-固-热三方耦合问题时因未充分利用系统矩阵结构特性而导致的计算效率低下问题；

34、2. 通过将流体力学方程、固体力学方程和热力学方程在统一的任意拉格朗日-欧拉框架下构造为以流体速度、压力、固体位移和温度为未知量的单块耦合变分模型，使流固交界面处的速度连续与应力平衡条件在该模型中自动满足，实现了对地幔对流驱动岩石圈变形过程中热驱动、力传递和变形反馈的闭环耦合模拟，从根本上克服了现有技术因将流体、固体与热力学分开处理而无法完整刻画热-力-流反馈机制的缺陷；

35、3. 通过在时间推进过程中将基于数值稳定性条件的候选步长计算与对数时间尺度调节相结合，根据模拟进程所处的数量级区间自适应调整步长，并同时施加向后二阶差分格式时间离散、压力正则化处理和基于熵黏度的人工黏性处理，保障了百万年时间跨度积分过程的数值稳定性和精度，解决了现有技术难以兼顾长时间积分效率与多物理场耦合数值稳定性的技术问题。