用于动态Biot模型的时连续强守恒时空有限元方法

《Mathematics and Computers in Simulation》：Time-continuous strongly conservative space–time finite element methods for the dynamic Biot model

　　摘要
我们研究了动态Biot模型（参见Biot，1962年），该模型描述了流体流动与固体变形之间的相互作用，包括饱和多孔介质中液体相和固体相中的波传播现象。该模型将动量平衡的双曲方程、二阶时间动态达西定律以及质量平衡的抛物线方程相结合，此处采用三场形式进行研究，关注的物理场为弹性矩阵的位移、流体速度和流体压力。本文提出了一类变分时空有限元方法，该方法结合了任意多项式阶数的连续时间Galerkin近似、位移场及其时间导数和通量场的H(div)一致逼近——对于位移则采用不连续Galerkin方法——以及分段多项式压力近似，从而在每种情况下都得到了一个inf-sup稳定的强守恒混合方法。我们证明了在空间上的组合能量范数下，时间上最大范数的误差估计。理论结果通过针对不同空间和时间多项式阶数的数值实验得到了验证。

引言
多孔弹性理论可追溯到上个世纪中叶，由Maurice Anthony Biot在一系列开创性论文中进行了深入发展，参见[1]、[2]、[3]、[4]、[5]。该理论始于固结的准静态模型[1]，随后从各向同性介质扩展到各向异性介质[2]、粘弹性材料[3]，并进一步发展到动态系统[4]、[5]。文献[4]、[5]、[6]、[7]中的模型源自非平衡热力学的普遍原理，描述了流体流动与（线性、各向异性）多孔介质变形之间的相互作用，同时考虑了波传播、能量耗散以及其他松弛效应。关于Biot所提模型的适定性、离散化及稳定性的首批基本数学结果可见于[8]、[9]、[10]、[11]、[12]、[13]、[14]，更全面的概述可参见[15]、[16]。

（热）多孔弹性的经典应用领域是岩土力学，典型问题包括评估石油开采引起的地面沉降、地下储层模拟、环境规划，以及放射性废物安全处置屏障的设计，参见例如[17]、[18]、[19]、[20]、[21]。近年来，人们发现生物组织可视为被血液和间质液等流体渗透的多孔、透水且可变形的介质，参见例如[22]、[23]，这一发现促使人们成功地将多孔弹性模型应用于研究人体各种器官的生物力学行为，如大脑[24]、[25]、[26]、肝脏[27]、肾脏[28]，这些研究有助于更好地理解从肾衰竭、创伤性脑损伤到癌症等多种疾病状态。

在所研究的模型中，尤其是生物医学应用中，物理参数的多样性使得设计和分析合适的数值方法成为一项挑战。正因如此，过去几十年来人们一直在努力开发稳定且对参数具有鲁棒性的方法，用于数值求解这些相关的耦合偏微分代数方程。关于Biot固结模型两场形式的有限元近似的严格稳定性与收敛性分析，首次见于[29]、[30]。所得到的先验误差估计既适用于半离散格式也适用于全离散格式，在这类格式中，空间离散化采用inf-sup稳定的有限元，时间推进则使用向后欧拉法。此外，还引入了一种后处理技术以提高孔隙压力的精度。此后，有许多研究探讨了这类问题的无锁定及（局部）守恒离散化的各个方面，最近还有研究采用了广义和扩展有限元方法，参见例如[31]，以及富集Galerkin方法，参见例如[32]、[33]，还有研究致力于分析裂隙异质多孔介质中的流速与波传播现象，参见[34]。直到最近才得到了准静态Biot问题的最优L2误差估计，参见[35]。

对于经典三场形式的准静态Biot问题，最初采用标准连续Galerkin近似处理位移场，同时采用混合近似处理通量-压力对——其中通量变量通过标准达西定律引入系统——这类离散化方法最早见于[36]、[37]。这些工作中的误差估计既适用于连续时间模型也适用于离散时间模型。这种方法后来也被扩展到位移场的不连续Galerkin近似，参见[38]，以及其他非一致近似方法，例如使用改进的旋转双线性元素[39]，或用于位移的Crouzeix–Raviart元素[40]。最近在[41]中，证明了一类基于H(div)一致不连续Galerkin近似位移场的强质量守恒离散化方法，在适当选取的参数相关范数下具有均匀稳定性，能够得到近乎最优的近似结果，且对这些模型参数和离散化参数（拉梅参数、渗透率与Biot-Willis参数、储存系数、时间步长和网格尺寸）都具有鲁棒性。相关连续时间与离散时间方案的误差估计见[42]，混合化技术则见[43]。在[44]中提出了一种针对三场形式Biot问题的稳定一致化方法。

文献[45]提出了准静态Biot问题的四场形式，该形式将两个混合问题耦合在一起，对于力学子问题采用基于Hellinger–Reissner原理的混合有限元方法，而对于流动子问题则采用标准的混合有限元方法。分别利用Arnold-Winther空间和Raviart–Thomas空间，为这两种耦合的混合有限元方法证明了半离散和全离散方案的最优先验误差估计。同一公式下的相同离散化方法在[46]中得到了进一步分析，给出了时间上的L∞范数误差和空间上的L2范数误差估计，这些估计对拉梅参数具有鲁棒性，并且在受限储存系数趋于零的极限情况下依然有效。

这些技术中有不少也被应用到多网络多孔弹性理论（MPET）中，研究重点包括连续模型与离散模型的稳定性、质量守恒、混合化、近乎最优近似结果、统一预处理以及分裂方案，参见例如[43]、[47]、[48]。其他三场、四场和五场模型还涉及总压力及其离散化、先验误差分析以及鲁棒预处理器，参见例如[49]、[50]、[51]、[52]。

迄今为止提到的大多数研究——前两段中的部分参考文献除外——仅提出并分析了准静态模型的数值方法。尽管这些方法多年来已被深入研究且已相当成熟，但关于动态Biot模型的稳定数值解方法的了解却少得多。不过，已有少量早期开创性的理论成果，涉及这类问题有限元近似的构建与误差分析，参见例如[53]、[54]、[55]。

出于对固体相和流体相动力学建模的兴趣，我们采用了在[4]、[5]、[7]中提出的经典三场形式动态模型，也可参见[14]。与[56]中依赖标准达西定律的模型不同，我们所研究的模型考虑了流体惯性效应，这就要求质量平衡方程（2.1c）基于通过动态达西定律（2.1b）引入的通量变量。本工作的重点与主要贡献在于，为这种更为一般的模型提出并分析一类新的守恒变分时空有限元方法。该离散化方法采用H(div)一致近似，对于通量采用一致近似，而对于位移场则采用不连续Galerkin类型的非一致近似，从而延续了[41]的方法。空间离散化的有限元空间选取方式与[41]、[57]一致，能够确保质量在点态上的守恒。我们用于时间离散化的连续时间Galerkin-Petrov方案也在[56]中有所研究，但如前所述，那是针对两场形式下的不同双曲-抛物线问题。收敛性分析中的证明思路借鉴了[56]、[58]中的观点。

本文的其余结构如下：第2节阐述了所研究的动态Biot问题的三场形式。第3节总结了本文后续会用到的部分符号，回顾了求积公式，介绍了插值算子以及一些相关的辅助结果。第4节给出了变分公式，这些公式是离散模型及后续分析的起点，并在第4.4定理中阐述了本文的主要结果。第5节详细推导了这一先验误差估计，首先回顾一些已知的辅助结果，接着证明了一些与离散化误差估计相关的界，最后给出第4.4定理的证明。第6节通过数值实验展示了我们的理论成果。最后，以简短的结论结束本文。

问题表述
在本文中，我们研究了如[1]、[7]中所提出的、适用于线性弹性、均匀、各向同性、多孔固体且充满微不可压缩流体的动态Biot模型，即以下耦合的偏微分方程组：
$$
\rho?\frac{\partial^2 u}{\partial t^2} - 2\mu \text{div}(\varepsilon(u)) - \lambda \nabla \text{div}(u) + \alpha \nabla p + \rho f \frac{\partial u}{\partial t} = f \in \Omega \times [0,T],
$$
$$
\rho f \frac{\partial^2 u}{\partial t^2} + \rho w \frac{\partial u}{\partial t} + K^{-1}w + \nabla p = 0 \in \Omega \times [0,T],
$$
$$
s_0 \frac{\partial p}{\partial t} + \alpha \text{div}(\frac{\partial u}{\partial t}) + \text{div}(w) = g \in \Omega \times [0,T],
$$
其中研究域为有界Lipschitz域$\Omega \subset \mathbb{R}^d$，$d \in \{2,3\}$，$T > 0$。模型中的物理参数为总密度$\rho? = (1 - \phi_0)\rho_s + \phi_0\rho_f > 0$。

函数空间
在本文中，我们采用Sobolev空间的标准记号。L2内积用$\cdot,\cdot$表示。在我们的分析中，需要用到Bochner空间$L^2(J;B)$、$C(J;B)$、$C_m(J;B)$，其中$m \in \mathbb{N}$，$B$为Banach空间，$J \subseteq [0,T]$。在后续内容中，向量值和矩阵值量将用粗体标出。

鉴于所提出的时间步进方法，我们将时间区间$[0,T]$划分为$N$个子区间$[t_{n-1}, t_n]$，长度为$\tau_n = t_n - t_{n-1}$，其中$0 = t_0 < t_1 < \cdots < t_N = T$。此外，我们还定义时间网格为$\{I_1, \cdots, I_N\}$。

变分公式
首先，我们引入变量$\frac{\partial u}{\partial t}$，并将方程(2.1)改写为：
$$
\rho?\frac{\partial u}{\partial t} - 2\mu \text{div}(\varepsilon(u)) - \lambda \nabla \text{div}(u) + \alpha \nabla p + \rho f \frac{\partial u}{\partial t} = f \in \Omega \times [0,T],
$$
$$
\frac{\partial u}{\partial t} - v = 0 \in \Omega \times [0,T],
$$
$$
\rho f \frac{\partial u}{\partial t} + \rho w \frac{\partial u}{\partial t} + K^{-1}w + \nabla p = 0 \in \Omega \times [0,T],
$$
$$
s_0 \frac{\partial p}{\partial t} + \alpha \text{div}\left(\frac{\partial u}{\partial t}\right) + \text{div}(w) = g \in \Omega \times [0,T].
$$
接下来，我们将给出系统(4.1)的连续时间变分公式和离散时间变分公式。后者是通过在时间上应用连续Galerkin-Petrov方法得到的。我们通过结合混合方法和H(div)一致近似来实现强守恒的空间离散化。

误差分析
在本节中，我们将为问题4.2中提出的时空变分方法的离散化误差推导出估计值。本分析以[56]、[58]为依据，但同时详细介绍了所研究的动态Biot模型及其强守恒空间离散化方法，该方法针对流动子问题采用混合方法，而针对力学子问题则采用H(div)一致不连续Galerkin方法。

数值实验
在本节中，我们研究了第4.2节中提出的连续Galerkin-Petrov（cGP）方法在时间上的数值收敛性。我们测试的时空离散化方法族为$cGP(k)$，其中$k \in \{1,2\}$，且$\ell = k$，即时间离散化采用k阶分段多项式，而在空间上则用$\ell = k$阶多项式来确定分段多项式的局部空间，例如空间$P_\ell(K)$，详见第4.2节。接下来，我们将研究离散化误差$e_u = u - u_\tau$。

结论
在本文中，我们提出了一类用于动态Biot模型的连续时间强守恒时空有限元方法。为此，我们采用了Galerkin-Petrov方案进行时间离散化，从而得到连续时间解。我们还将该方案与用于位移的不连续Galerkin方法，以及用于流体通量、流体压力和位移时间导数的混合空间离散化方法相结合。对于这种时空离散化方案，……