Cosserat非常状态基于非局部动力学方法用于岩石类材料的脆性断裂研究

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：Cosserat non-ordinary state-based peridynamic method for brittle fracture in rock-like materials

　　熊志辉|周书伟|何思阳|庄小英|蒂蒙·拉布丘克（Timon Rabczuk）
中国上海200092，同济大学土木工程学院岩土工程系

**摘要**
本研究提出了一种基于Cosserat非局部状态的周动力学模型（CNOSBPD），用于模拟类岩石材料中的弹性-脆性断裂过程。通过引入材料点的独立微旋转自由度，基于应变能量密度等价原理推导出了与旋转相关的刚度系数的显式表达式。此外，还为Cosserat模型建立了一个广义的混合模式断裂准则，以同时确定拉伸和剪切断裂。通过数值示例验证了该模型，包括具有孔洞的板材的单轴加载、双缺口梁的四点弯曲以及预裂岩石试样的单轴压缩。结果表明，所提出的CNOSBPD模型能够合理再现裂纹起始、扩展和不稳定演化的特征，在裂纹偏移、裂纹合并以及承载能力退化方面与实验观察结果吻合良好。与忽略微旋转的非局部状态周动力学模型（NOSBPD）相比，该模型增强了对裂纹起始和扩展过程中剪切效应的描述，从而使得断裂形态和演化过程的预测更加符合实验观察结果。

**引言**
研究类岩石材料中的断裂机制对于评估岩土稳定性及确保资源开采的安全性至关重要。阐明岩石破坏前的前兆现象对于采矿活动和岩土工程项目的安全运行至关重要[1]。然而，岩石介质通常存在孔隙和微裂纹等内在缺陷，在外部载荷作用下裂纹可能在不同位置起始并随后扩展和合并，导致高度复杂的破坏模式和演化过程。因此，可靠预测裂纹扩展路径和断裂模式一直是岩石力学中的主要挑战[2][3]。此外，多裂纹相互作用和尺度效应的耦合进一步复杂化了这些断裂机制的定量表征[4][5]。

实验方法有助于研究岩石断裂机制。单轴拉伸、巴西圆盘分裂和三点/四点弯曲等测试能够直接观察裂纹的起始和扩展，并允许测试关键力学参数，如断裂韧性[6]。这些实验为理解岩石断裂行为提供了直观的证据。然而，实验研究也有局限性：受试样尺寸和边界条件的限制，实验室得到的断裂参数在不同尺度上可能存在显著差异，难以直接外推到实际工程尺度[7]。此外，实验通常成本高昂且耗时较长，由于测量技术的限制，难以全面捕捉加载过程中岩石内部的裂纹演化，包括多源起始、相互作用以及裂纹贯穿发展。因此，仅依靠实验观察不足以完全阐明岩石断裂机制，需要可靠的数值模型来提供补充分析和机理解释。

为了克服实验研究的局限性，已经开发了多种数值方法来模拟岩石断裂过程[8]。传统的连续介质方法，如有限元方法（FEM），使用网格来描述岩石的力学响应。然而，随着裂纹的扩展，网格可能会发生严重畸变，从而影响计算稳定性。因此，通常需要重新网格划分或裂纹跟踪方案来保持准确性[9]。改进的公式，包括扩展有限元方法（XFEM）和数值流形方法（NMM），部分解决了这些跟踪问题，但其实现仍然复杂且计算成本较高[10]。无网格方法，如平滑粒子动力学（SPH），减少了网格依赖性，在大变形问题中表现良好。然而，它们在施加边界条件和确保稳定性方面仍面临挑战[11]。类似地，尽管场方法越来越多地用于模拟复杂裂纹模式而无需显式跟踪，但它们通常需要仔细的参数校准，并且对于三维问题需要较高的计算资源[12]。另一方面，离散元方法（DEM）通过将岩石离散为颗粒或块体，天然适合表示节理和裂纹。然而，这些方法需要对接触模型进行广泛的校准，并可能难以捕捉整体的连续介质响应[13]。总体而言，大多数现有方法仍需要一定程度的用户干预，如重新网格划分或预定义裂纹路径，这激发了开发更健壮的断裂模拟方法的需求。

周动力学（PD）作为一种非局部数值方法，在模拟断裂和破坏方面受到了广泛关注[14]。通过用积分形式的运动方程替代经典力学中的偏微分方程，PD避免了连续可微位移场的要求。这一特性使其非常适合表示岩石和地质材料中的不连续性，如裂纹。基于键合的PD（BPD）模型最初由Silling和Askari提出[15]，通过材料点之间键合的断裂来表示材料破坏。然而，BPD受到固定泊松比的限制，这限制了其应用范围。最近的研究也从非局部变形场的角度重新审视了BPD，以阐明其潜在的连续介质结构[16]。为了解决这些限制，Silling等人提出了基于状态的PD框架[17]。该框架分为基于普通状态的PD（OSBPD）和基于非局部状态的PD（NOSBPD）。尽管OSBPD消除了固定的泊松比限制，但它仍然假设相互作用力与键合变形共线。此外，OSBPD通常需要根据“力状态”重新推导本构关系，这使得直接使用成熟的经典模型变得困难[18]。相比之下，NOSBPD引入了非局部变形梯度，将键合力映射到经典的柯西应力张量上，允许非共线的相互作用力并直接集成经典本构关系。这些特性使NOSBPD在处理复杂应力状态时具有明显优势[19]。然而，需要注意的是，NOSBPD模型容易受到虚假零能模式和数值不稳定性的影响，因此通常需要适当的稳定化或改进的公式来确保模拟的可靠性，尤其是在剪切敏感的混合模式断裂问题中[20][21]。除了岩石力学，PD还成功扩展到更广泛的复杂断裂问题，包括近乎不可压缩的超弹性裂纹扩展、弹塑性脆性断裂和强度驱动的脆性断裂[22][23]。这展示了其在不同材料类别和变形机制下的卓越适用性。

Cosserat理论为通过纳入内部微观结构效应来模拟岩石断裂提供了新的视角。在经典连续介质力学（CCM）中，通常假设材料点没有独立的内部旋转自由度。相反，Cosserat连续介质为材料点赋予了额外的局部旋转自由度以及平移自由度。该框架允许存在耦合应力，定义为单位面积上的力偶。因此，材料点可以同时进行平移和独立旋转[24]。通过引入这些微旋转自由度及其相关刚度，Cosserat理论可以捕捉由微观结构效应（如尺寸效应）驱动的宏观力学差异[25]。在含有孔洞或裂纹的介质中，Cosserat模型通常预测微旋转效应可以减少应力集中，从而提高表观韧性[26]。由于岩石通常具有颗粒骨架、胶结和天然裂缝等微观结构元素，将Cosserat理论纳入岩石断裂分析有助于更真实地表示局部旋转和局部变形。这种方法对于解释驱动断裂演化的微观机制至关重要。

近年来，一些研究人员将微观结构诱导的旋转效应整合到PD框架中，以更好地表示岩石断裂。Gerstle等人提出了早期的微极向PD概念，通过在基于键合的相互作用中引入力偶项来描述微旋转[27]。在此基础上，后续研究开发了一系列改进的模型[28][29][30]。例如，Chen等人和Guo等人分别在BPD框架内开发了多种Cosserat-PD（CPD）模型[31][32]。通过考虑横向相互作用和力偶效应，这些模型增强了颗粒旋转对裂纹扩展的影响描述。与此同时，Zhou等人提出了基于OSBPD框架的Cosserat基于普通状态的周动力学模型（COSBPD）[33]。通过引入Cosserat旋转刚度，该模型对预裂岩石试样的裂纹扩展预测与实验观察结果更为吻合。然而，大多数现有模型仍然基于键合或基于普通状态的框架构建。它们的断裂评估通常依赖于关键的键合拉伸准则。这种基于几何的准则无法有效区分拉伸断裂和剪切滑移等其他破坏模式。在剪切主导或压缩-剪切组合载荷下，它倾向于低估剪切裂纹的起始和扩展。为了解决这个问题，Wang等人和Zhou等人将基于应力-强度的断裂准则引入了NOSBPD框架，并考虑了拉伸-剪切区分[34][35]。这种方法在一定程度上改进了破坏模式的识别。尽管如此，仍有改进的空间，以更准确地区分剪切相关的裂纹和一致地表征混合模式断裂演化。总体而言，需要进一步努力将Cosserat微旋转效应整合到更通用的基于状态的框架中，并开发出能够同时表征拉伸和剪切断裂的断裂准则。

为了解决这些挑战，本研究基于Zhou等人的工作，提出了一种基于Cosserat非局部状态的周动力学（CNOSBPD）模型[33]。通过引入独立的微旋转自由度，该模型通过Cosserat能量分解和应变-能量密度等价原理确定旋转刚度系数。这些系数考虑了微曲率和反对称剪切贡献，使Cosserat微旋转效应能够一致地集成到NOSBPD框架中。因此，所得到的运动方程能够准确表示岩石微观尺度单元中剪切变形和旋转之间的耦合。与COSBPD模型和其他现有的CPD方法相比，CNOSBPD利用了非局部框架的本构映射能力，使模型能够更真实地捕捉局部旋转响应和应力重分布，同时允许直接使用经典的应力概念。在此基础上，提出了CNOSBPD的广义莫尔-库仑断裂准则。该准则克服了传统PD准则仅依赖关键键合拉伸的局限性，显著提高了模型预测剪切主导和混合模式断裂的能力。通过包括具有单个或多个预先存在的缺陷的试样的单轴压缩以及预裂梁的四点弯曲等代表性数值示例，验证了CNOSBPD模型。

本文的其余部分组织如下：第2节简要回顾了基于非局部状态的周动力学和Cosserat连续介质理论的基础知识。第3节详细介绍了CNOSBPD模型的理论推导和本构构建，还介绍了所提出的广义莫尔-库仑断裂准则及其数值实现。第4节通过将数值结果与代表性基准测试的实验数据进行比较来验证模型，包括单裂纹压缩、双裂纹压缩和四点弯曲测试。第5节讨论了所提出模型的适用性及其当前的限制。最后，第6节总结了主要结论并指出了未来的研究方向。

**NOSBPD理论**
NOSBPD是一种连续介质力学理论，以积分形式描述材料点之间的非局部相互作用。在这个框架中，材料点x处的运动方程可以表示为[36]：
$$ \rho u?(x,t) = \int_H^x T_{\xi,t} - T_{-\xi,t} \, dV \, x' + f(x,t) $$
其中$\rho$表示材料点的质量密度，$u?$是材料点的加速度向量；$H_x = x'|x' \leq \delta$表示以x为中心、视野半径为$\delta$的家族（非局部积分域）；$f(x,t)$是作用在材料上的体力密度。

**CNOSBPD模型**
对于二维问题，CCM的应变能量密度（SED）可以表示为：
$$ W_{CCM} = \frac{1}{2} \sigma_{xx} \epsilon_{xx} + \sigma_{xy} \epsilon_{xy} + \sigma_{yy} \epsilon_{yy} $$
相应地，Cosserat连续介质的应变能量密度（SED），用上标C表示，表示为：
$$ W_{CCM} = \frac{1}{2} \sigma_{xx}^C \epsilon_{xx}^C + \sigma_{xy}^C \epsilon_{xy}^C + \sigma_{yy}^C \epsilon_{yy}^C + \sigma_{yx}^C \epsilon_{yx}^C + \mu_{zx} \gamma_{zx} + \mu_{zy} \gamma_{zy} $$
这里，$\sigma_C$和$\epsilon_C$表示Cosserat连续介质中材料点的应力和应变张量。在实际应用中，这些量通常被近似为相应量的有效平均值。

**数值示例**
本节通过几个代表性边界值问题系统地验证了所提出的CNOSBPD模型，重点是其预测弹性响应、裂纹起始、扩展路径、裂纹相互作用和最终破坏模式的能力。除非另有说明，以下数值示例中使用的常规力学参数和断裂相关参数均来自相应的实验数据或基准参考文献。

**讨论**
数值结果表明，显著的微旋转主要集中在预先存在的缺陷尖端和施加边界条件的区域附近，这与高应力梯度和规定的位移约束密切相关。从物理上看，岩石材料固有的异质颗粒结构在加载下会引起局部旋转。然而，传统的NOSBPD方法通常基于柯西连续介质假设，忽略了某些重要的因素。