基于双曲型热力学兼容模型的一种异质多孔介质计算均质化方法

《Journal of Computational and Applied Mathematics》：A method of computational homogenization of a heterogeneous porous medium based on a hyperbolic thermodynamically compatible model

　　摘要：本文提出了一种评估含流体饱和的可变形多孔材料有效弹性模量的方法。该方法基于计算非均匀介质的应变能，并将其与等效均匀弹性介质的能量进行比较。对于多孔介质，利用双曲热力学相容方程组来描述其非均匀性，该方程组在小变形情况下进行了简化；这种表述方式使得在存在纯弹性固体或纯流体区域时能够采用扩散界面模型。通过求解一系列简单的静态应力-应变问题，比较多孔介质及其等效弹性介质的能量，就可以确定有效柔顺张量的所有分量，进而得到有效刚度张量的所有分量以及相应的有效波速。为了解决控制方程的数值问题，采用了在交错网格上的高阶有限差分方案。为了获得样品的静态平衡解，通过结合人工速度松弛的时间推进技术来实现稳态状态。一系列测试问题表明，该方法在应用于不同类型的饱和多孔样品时具有很强的灵活性。此外，还将该方法与已知理论预测的结果进行了比较。特别是对于具有简单周期性微观结构的液体饱和多孔弹性样品，HTC模型中使用的弹性模量能够被数值计算得到的有效模量准确近似。

引言：利用数值建模来确定岩石性质是数字岩石物理学领域中快速发展的一项技术，这一领域涵盖了广泛的研究内容，从获取岩石岩芯的三维数字断层图像及其后续处理（过滤和分割），到开发用于计算岩石样品各种性质的方法和算法，如孔隙度、渗透率、弹性模量、纵波和横波速度等，可参考Andr?等人[1]、[2]，Kang等人[3]，Liu等人[4]及其中的参考文献。本文的目的是开发一种算法，用于确定二维和三维非均匀介质的有效弹性模量，尤其是饱和多孔介质。确定有效模量的问题始于20世纪70年代（可参考Christensen[5]、Glüge[6]及其中的参考文献），当时由于需要描述复合材料的性质，人们开始发展一些理论方法来研究具有规则周期性结构的介质。我们将研究饱和多孔介质的有效弹性模量问题，这类介质包括岩石等。这类介质的特点是微观非均匀性分布不规则，无法用规则的周期性或准周期性结构来描述。为了确定多孔介质的有效弹性参数，自然会采用基于变形能量计算的数值方法，这类方法近年来得到了快速发展，并已成功应用于评估微结构介质和蜂窝材料的弹性模量[7]、[8]、[9]、[10]。下面所介绍方法的本质是在简单应力状态下对非均匀样品的变形能量进行数值计算。之后，假设给定材料对应于一个等效弹性介质，便可计算出其有效柔顺模量。我们方法的关键在于选择合适的数学模型来计算所研究材料的能量。由于我们的研究对象是饱和多孔介质，因此将采用在[11]、[12]中提出的、基于双曲热力学相容模型[13]、[14]构建的、用于描述两种可压缩流体混合体饱和的可变形多孔介质中的波动过程模型。属于HTC系统类的连续介质力学方程满足热力学定律（能量守恒和熵增原理），并且具有良好的数学特性（初始边界值问题的适定性），可参考[15]、[16]及其中的参考文献。此外，为了解决双曲守恒律系统而开发的先进高精度数值方法可以直接应用于解决这些问题，而HTC方程的良好特性则保证了所得数值解的可靠性。关于HTC模型的高阶方法的发展和应用的一些成果可见于[17]、[18]、[19]、[20]、[21]、[22]、[23]。另外，还有一些近期论文探讨了连续介质力学模型的热力学相容格式的发展[24]、[25]、[26]、[27]、[28]。在这些研究中，熵不等式被作为主要的演化方程，而非通常使用的总能量守恒定律（这与大多数传统的双曲偏微分方程格式有所不同），总能量守恒定律则是热力学相容离散化的结果。在饱和多孔介质的HTC模型中，该介质被视为一个宏观多相连续体，其行为由各相的质量、动量和能量平衡定律来描述，同时还会考虑各相之间的相互作用。许多研究人员都采用基于多相介质守恒定律和热力学定律的连续介质方法，可参考Marle[29]、Wilmanski[30]、Pesavento等人[31]的研究。HTC模型的一个重要特点是，其控制方程是从第一性原理推导出来的，具有数学上的正确性。实际上，这些控制方程可以通过在拉格朗日坐标系中寻找作用量泛函的极值，从欧拉-拉格朗日变分方程中得到，而平衡定律中的通量结构也是唯一确定的[16]。此外，多孔介质的HTC模型是所谓连续统一模型的推广形式，该模型使用一套控制方程来描述介质的弹性变形、弹塑性变形以及作为流体的流动[17]、[18]。完整的非线性HTC模型可用于推导饱和多孔介质中小振幅波传播的方程，这一方程可以作为在[32]、[33]中提出的、在地震学领域被广泛接受且应用广泛的Biot模型的替代方案（可参考[34]、[35]及其中的参考文献）。HTC模型能够定性描述Biot理论所描述的波场特性，比如快压缩波和慢压缩波[13]。本文还对单一流体饱和多孔介质中波的色散曲线进行了定量比较，结果显示当孔隙度为0.2时，HTC模型和Biot模型的结果非常接近。重要的是，HTC模型允许使用统一的方程组结合扩散界面方法，来模拟具有复杂非均匀结构的区域中的波场，这些区域包含纯固体、纯液体和多孔相等子区域，而无需明确指定界面边界[36]。这一优势使得HTC模型相比Biot模型更具优越性，因为Biot模型在孔隙度极低或极高的极限状态下，会错误地描述纵波速度[13]。还需指出的是，将Biot模型扩展到具有复杂流变特性的介质，如有限变形、多相饱和流体、结构损伤等，会面临诸多困难，而本文所提出的线性波场模型则是从完整的HTC模型衍生而来的，在HTC框架内可以实现对模型复杂性的处理。在本文所提出的多孔介质波场线性模型中，通过引入扩散界面方法，可以模拟体积分数为0或1的“纯”固体相和液体相之间的界面。在常用的Baer-Nunziato型两相流模型中，从最初的开创性工作[37]到后续的研究，都假定“纯”相中包含少量另一种相，以避免在求时间导数时出现矩阵奇异性。而在我们的线性HTC模型中，完整方程组中并不包含体积分数的方程，而且在介质稳态状态下，空间导数系数中的体积分数初始值在将扩散界面方法应用于其0或1的极限值时也不会出现问题。由于使用基于应变能的方法在简单应力状态下计算弹性模量需要求解静态方程组，因此我们采用时间推进法来获得模型方程的静态解。该方法是通过在方程中加入额外的耗散性速度松弛源，然后使用交错网格上的有限差分法来数值求解所得方程组。该方法通过计算已知弹性模量的介质的柔顺张量来验证其有效性。随后通过解决一些测试问题，我们得出结论：所使用的模型能够很好地匹配所研究的合成多孔弹性介质，因为模型中采用的弹性模量与合成样品的有效模量十分接近。需要指出的是，目前已有几种成熟且完善的数值均质化方法，可用于计算非均匀介质的参数。对于具有周期性结构的介质，可以利用快速傅里叶变换方法，借助周期性特点提高计算效率，从而计算出复杂微观结构的有效宏观性质[38]、[39]。另一种用于非均匀介质的成熟方法是基于有限元方法的数值均质化方法（可参考Allaire和Brizzi[40]、Otero等人[41]、Wei等人[42]、Khoei等人[43]的研究）。这种方法是通过在包含多孔介质微观结构的代表性体积元素上求解边界值问题，来确定有效的宏观性质。这种方法需要仔细处理孔隙级别的相界面，而且计算成本较高。变分均质化方法则利用变分原理来获得材料有效性质的严格上下界（可参考[44]、[45]、[46]及其中的参考文献）。在这些方法中，均质化问题被表述为在具有周期性约束的函数空间中的最小化问题。变分边界条件能够提供可靠且精确的估计值，有助于评估数值结果的误差和质量。本文所发展的方法将应变能计算与时间推进法相结合，用于获得HTC系统的稳态平衡解，这种方法在数值计算有效弹性模量方面展现出很大的潜力。为了解决该问题，采用了高效的有限差分方案，该方案不会导致过高的计算成本，同时还便于实现高效并行化。在本文中，我们将研究范围限制在由线性HTC系统描述的小变形情况。这是因为，在计算非均匀介质的有效弹性模量时，我们的目标是得到一个封闭形式的超弹性变形过程模型。这样的模型在数学和物理上都是正确的，可应用于解决地震学中的实际问题、生物材料的开发，以及制造具有特定性能的复合材料。当然，研究有限变形条件下可变形非均匀介质的有效性质也具有重要意义。然而，这个问题极为复杂，因为要建立非均匀介质的超弹性模型，就需要确定应选择哪种作为变形非线性函数的弹性势。在这种情况下，基于变形能量方法，就可以得到取决于当前应力状态和应变状态的有效弹性模量。有效模量张量可以被定义为应变能密度对宏观应变梯度的二阶导数，但它取决于计算该张量值时的应力值和应变值（可参考[47]）。因此，在有限变形的情况下，计算有效模量的问题伴随着大量不确定性。此外，还需要指出的是，试图应用本文提出的方法，即通过时间推进法求解多孔介质的完整HTC模型，将会面临巨大困难。这是因为，当在样品边界施加载荷时，样品会发生有限变形，这会导致相界面发生位移、形状改变以及孔隙度变化。为解决这类问题，有必要运用近期开发的高精度数值方法，例如Thomann和Dumbser[22]、Busto和Dumbser[28]、Río-Martín等人[48]针对多孔介质的HTC方程所采用的方法，这也是未来研究的一个可能方向。总之，本文提出了一种基于以下内容的新型数值均质化方法：•用于小变形情况的多孔介质HTC模型；•基于应变能的方法，利用HTC模型评估非均匀多孔介质的弹性模量；•一种带有人工速度衰减的时间推进法，用于求解孔弹性能方程，从而得到平衡问题的解。本文的其余结构如下：第2节介绍了充满两种流体混合物的变形多孔介质的HTC控制方程，并分析了其特性，同时还给出了平行六面体域中最简单应力状态的边界条件。第3节探讨了基于应变能确定柔顺度张量和刚度张量的方法。第4节展示了数值模拟的结果，包括对模型和方法的验证，以及对一些合成多孔介质样本以及通过断层成像得到的三维岩心样本的测试问题求解结果。文章最后为结论部分。

模型控制方程
我们将研究在充满两种流体混合物、尤其是液体与气体混合物的可变形多孔介质中，小振幅波的传播规律。此处所考虑的模型是对多孔介质通用HTC多相模型的线性化处理[12]、[13]，它能够生成适用于小振幅波场建模的微分方程，且属于一阶线性双曲型系统，可应用高效的高精度数值方法进行求解。

弹性能量与平衡状态
我们的目标是通过對加载后的平行六面体样本进行简单的数值实验，开发出一种确定有效刚度张量的方法。假设存在一种弹性介质，其一般的应力-应变关系遵循胡克定律σij=cijklεkl，其中cijkl为刚度张量。该关系式也可表示为εij=sijklσkl，此时sijkl为柔顺度张量，它是刚度张量的逆矩阵。应变弹性能量密度可表示为U=1/2σijεij=1/2cijk。

数值模拟
在本节中，我们通过一系列数值实验来验证文中提出的理论推导。用于模拟的交错网格有限差分方案借鉴了Virieu[54]、Moszo[55]、Graves[56]的研究成果，他们提出了一种基于速度-应力形式求解线性弹性理论方程的方法。该方法非常适合处理结构类似的双曲型系统，而我们也据此发展出了四阶数值方法。

结论
确定具有复杂微观结构的非均匀介质的有效参数，对于构建精确的连续介质模型至关重要，尤其是在地震勘探领域。尽管传统的均质化理论能够为周期性弹性复合材料提供可靠的估算结果，但它们并不直接适用于真实的、充满流体的多孔材料。为解决这一问题，就需要采用计算均质化方法，如数字岩心分析技术。

作者贡献声明
Evgeniy Romenski：写作——审阅与编辑，写作——初稿撰写，方法论，研究，概念设计。Galina Reshetova：写作——审阅与编辑，写作——初稿撰写，可视化，验证，软件应用，形式化分析。

利益冲突声明
作者声明不存在任何可能影响本文研究工作的已知财务利益或个人关系。

Evgeniy Romenski|Galina Reshetova
俄罗斯科学院索伯列夫数学研究所，科普秋格大街4号，新西伯利亚，630090