基于谱Galerkin方法研究Kelvin–Voigt型非均质材料中的波传播

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：Wave propagation in Kelvin–Voigt heterogeneous materials: A spectral-Galerkin approach

为模拟三维Kelvin–Voigt黏弹性非均质材料中的波传播现象，研究团队结合Cauchy方程与Kelvin–Voigt应力?应变本构关系建立了偏微分方程模型，并提出一种基于谱Galerkin的数值算法。该方法利用定义在整个计算域上的全局基函数逼近解，并通过在材料界面施加弱L2连续性来处理非均质性。通过Crank–Nicolson或强迫项精确积分的时间离散方案求解所得常微分方程组，一系列算例验证了该方法的有效性和灵活性，为复杂地质结构（如井孔）中的波传播模拟提供了高精度工具。

　　

在地球科学、工程学和数学等领域，波在材料中的传播现象一直是研究热点。特别是在地球物理应用中，如井孔成像，了解声波如何在井孔及周围岩石中传播至关重要。这些传播过程涉及弹性波和黏弹性波，而实际地质材料往往同时表现出弹性和黏性行为，这类材料被称为黏弹性材料。当波在黏弹性材料中传播时，会发生机械能损耗，导致波幅衰减和波速频散。因此，需要能够描述这种耗散和频散效应的本构模型。在众多黏弹性模型中，Kelvin–Voigt模型因其简洁性和广泛的适用性而备受关注。然而，在三维、非均质复杂材料中对该模型进行高精度、高效率的数值模拟，仍然是一个挑战。现有的方法，如有限差分法、伪谱法、有限元法和谱元法，各有优势但也存在局限。例如，有限差分法计算高效但数值频散明显，且处理复杂界面和边界条件困难；伪谱法精度高但并行性差，且对非周期边界的处理复杂。如何发展一种兼具高精度、能灵活处理非均质和复杂几何、并便于实现的方法，是当前波传播模拟中的一个重要问题。为解决这一问题，发表在《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》上，由 Nicolò Crescenzio, Giulio G. Giusteri, Antonia Larese, Federico Piazzon, Mario Putti, Axel Turolla, Massimo Zampato 共同完成的研究，提出了一种新颖的数值算法，专门用于模拟三维、非均质、圆柱层状几何的Kelvin–Voigt黏弹性材料中的波传播。

本研究采用的关键技术方法主要包括：1. 构建了将Cauchy运动方程与三维Kelvin–Voigt黏弹性本构模型耦合的数学模型。2. 采用了谱Galerkin空间离散方法，通过全局模态基函数（而非网格点）逼近解，基函数依赖于圆柱坐标(r, θ, z)，由径向和轴向坐标的一类切比雪夫多项式组合与角度坐标的三角函数张量积构成。3. 为处理材料非均质性，将计算域划分为多个子域（每个对应一种材料），为每个子域定义独立的基函数集，然后在材料界面上施加弱L²连续性约束，通过QR分解稳定地构造满足全局连续性的新基函数集，其方法类似于Mortar元法。4. 采用Crank–Nicolson时间积分格式或对强迫项进行精确积分来求解所得的常微分方程组。

研究结果

2. Waves in Kelvin–Voigt viscoelastic heterogeneous media

本节正式给出了三维Kelvin–Voigt黏弹性材料的控制方程。动量守恒（Cauchy）方程为ρü = div(σ) + f_ext。本构模型将应力张量σ分解为弹性部分σ_el和黏性部分σ_vi之和，即σ = σ_el+ σ_vi。弹性应力由胡克定律给出：σ_el= λ_eltr(ε(u))I + 2με(u)，黏性应力由牛顿定律给出：σ_vi= λ_vitr(ε(?))I + 2ηε(?)。其中，ε(u)是应变张量，λ_el和μ是Lamé弹性系数，λ_vi和η是Lamé黏性系数。组合后得到完整的位移方程。研究还定义了描述材料衰减特性的品质因子Q_P和Q_S，它们与材料参数和波频率ω相关。

3. Weak formulation

为处理非齐次、随时间变化的Dirichlet边界条件，研究引入了提升算子L，构造了一个满足边界条件的函数L_uD。通过定义? = u - L_uD，将原问题转化为关于?的齐次Dirichlet边界问题，从而可以在与时间无关的Sobolev空间子空间（如[H₀¹(Ω)]^d）中寻找解，这为后续的变分（弱）形式离散和理论分析（存在性、唯一性、正则性）奠定了基础。

4. Spectral-Galerkin discretization in space and time

详细阐述了空间和时间离散过程。对弱形式方程进行空间离散，利用一组全局基函数{φ_j}_j=1^N对解进行近似。通过Galerkin投影，将连续的弱形式问题转化为一个关于时间二阶常微分方程的系统：M_ρ? + C? + Kw = F。然后，通过引入新变量v = ?，将此二阶系统降阶为一阶系统，便于使用标准的时间积分器（如Crank–Nicolson方法）进行求解。

5. Construction of the global modal basis functions

核心部分介绍了全局模态基函数的构造。对于研究的圆柱层状几何，基函数采用张量积形式：φ_l,m,n(r,θ,z) = R_l(r) Θ_m(θ) Z_n(z)。其中，角度部分Θ_m(θ)采用三角函数（如e^imθ），径向和轴向部分R_l(r)和Z_n(z)采用一类切比雪夫多项式的特定线性组合，以满足问题的边界条件（如轴对称、齐次Neumann条件等）。对于非均质材料，首先在每个材料子域上独立构造局部基函数。然后，在所有材料交界面上施加弱L²连续性条件，即要求解在界面两侧的L²内积相等。通过求解由此约束条件导出的线性系统（并使用QR分解以提高数值稳定性），得到一组系数矩阵，将局部基函数组合成一组满足全局弱连续的新的全局基函数。这种方法允许在不同子域（材料）上使用不同的分辨率（多项式阶数）。

6. Numerical examples

通过两个数值算例验证了方法的有效性。第一个算例是经典的Lamb问题的一个变体，模拟了均匀Kelvin–Voigt半空间表面受到垂直点力载荷时的波场。数值结果清晰显示了P波、S波和Rayleigh面波的传播，以及由于材料黏性导致的波幅衰减，验证了算法在模拟均匀黏弹性波传播方面的正确性。第二个算例模拟了波在非均质、多层圆柱介质中的传播，该几何模拟了井孔环境。结果成功追踪了波在多层介质界面处的反射和折射现象，准确捕捉了由材料属性突变引起的复杂波场，充分展示了该方法处理复杂非均质几何和材料界面的能力。

研究结论与意义

本研究成功开发并验证了一种用于模拟三维、非均质、Kelvin–Voigt型黏弹性材料中波传播的新型谱Galerkin数值算法。该方法的创新性主要体现在：1. 采用全局模态基函数（基于圆柱坐标的切比雪夫多项式与三角函数的张量积）进行空间离散，无需定义计算网格，具有谱方法的高精度特性。2. 通过引入弱L²连续性条件并利用Mortar元法的思想处理材料界面，能够自然、灵活地处理非均质分层材料，并允许子域间采用不同的多项式阶数（hp自适应）。3. 所构造的基函数天然满足圆柱几何的边界条件，特别适用于井孔地球物理等轴对称或近轴对称问题。

本研究的意义在于为解决地球物理勘探、地震工程和材料科学中复杂的波传播正演模拟问题提供了一个强有力的数值工具。相比于传统的有限差分法和伪谱法，该方法在复杂界面处理、边界条件实现和数值精度方面具有潜在优势。相比于谱元法，其采用的模态基函数和弱连续性条件提供了另一种处理非均质性的灵活途径。论文中展示的算例表明，该方法能够精确模拟波在黏弹性介质中的衰减、频散以及在复杂界面处的相互作用，对于理解和量化实际地质材料（如含流体岩石）中的波传播机理，以及优化地球物理测井和地震成像技术具有重要的理论和应用价值。未来的工作可以扩展该方法以包含更复杂的黏弹性模型（如Zener模型）、更一般的几何形状以及非线性效应。