岩土结构大变形断裂分析的相场物质点方法转让专利
申请号 : CN202110727024.7
文献号 : CN113360992B
文献日 : 2022-02-15
发明人 : 胡志强 , 郑勇刚 , 刘宇 , 叶宏飞 , 张涵博 , 郑章成 , 刘振海
申请人 : 大连理工大学
摘要 :
权利要求 :
1.一种岩土结构大变形断裂分析的相场物质点方法，其特征在于，实施步骤如下：为了准确捕捉高孔隙率岩土材料在大变形情况下的真实力学特性，将基于有限应变连续介质力学框架开展采用耦合相场断裂模型/光滑双屈服面Drucker‑Prager cap塑性本构模型分析压力敏感岩土材料的脆性‑塑性转变断裂破坏问题的研究，提供一种岩土结构大变形断裂分析的相场物质点方法PF‑ICPDI，具体步骤如下：首先，基于有限变形运动学框架，考虑一个任意固体域 内部存在一个相场c(x,t)∈[0,1]近似的非连续边界Γ，其外部边界 包括位移边界 和外部力边界 满足 通过运动函数 将参考构型Ω0中点X映射到t时刻的当前构型Ωt中点x，ndim∈{1,2,3}表示结构维度；则变形梯度表示为：对于有限弹塑性变形，变形梯度F通过乘式分解为弹性和塑性部分：e p
F＝FF (2)其中，上标e和p分别表示变量的弹性和塑性部分，因此，变形梯度的Jacobian分解为弹性和塑性部分：
e p
J＝det(F＝JJ (3)e e p p
其中，det(·)表示行列式算子，J＝det(F)和J＝det(F)；
基于各向同性假设和有限变形框架描述固体的有限弹塑性变形，因此，固体的弹性变形通过对称的弹性左Cauchy‑Green应变张量描述，即：e e eT
b＝FF (4)此外，固体域内每一点的应力状态通过对称的Kirchhoff应力张量τ表示，其具体的本构表达式通过能量存储函数Ψ表示为：e
对弹性左Cauchy‑Green应变张量b和Kirchhoff应力张量τ进行谱分解得：(A)
其中，n 表示特征向量， 和τA分别表示张量b和τ的特征值；
基于框架不变性假设，则系统自由能函数ψ仅依赖物质点X，且弹性主对数应变可定义为：
相应的体量弹性对数应变被表示为 偏量弹性对数应变被表示为 其中 和δ＝(1,1,1)；
e e
因此，Ψ(X,b)＝ψ(X,ε)，则主Kirchhoff应力τA表示为：一个任意变形体受到体积力ρ0G和边界 上的外部面力T作用，其中ρ0表示初始固体密度，G表示重力加速度，Γ为固体域Ω内被相场变量c(x,t)∈[0,1]近似的初始裂隙面，c取0表示完全断裂区域，取1表示完整区域，整个系统忽略惯性力效应并考虑一个准静态条件，则基于虚功原理、微观力平衡法则和热力学第二定律推导出如下强形式控制方程：系统动量平衡方程
DivP+ρ0G＝0 (10)相场演化方程
并且有P·N＝T为边界 上外力边界条件， 为边界 上相场边界条件，其中N表示边界 的外法线向量，Div(·)表示散度算子， 表示梯度算子，l0表示控制光滑近似裂纹宽度的正则化模型参数，0
从上式(11)看出，基于微观力平衡法则推导的岩土材料弹塑性断裂分析的相场演化方程与广泛使用的基于变分法则推导的脆性断裂相场演化方程在表达形式完全一致，仅仅区别在于历史应变场函数 的组分来源不同；此外，虽然两者具有相同表达形式，但是二者热动力学含义不同，并且本方法推导的相场演化方程退化为标准的相场脆性断裂模型；
然后，对上述强形式系统动量平衡方程(10)和相场演化方程(11)在物质点上采用隐式求解格式进行离散，并对系统离散控制方程采用交错迭代求解策略；首先，对系统动量平衡T
方程(10)采用Galerkin法在初始构型下构造弱形式控制方程，然后再利用τ＝P·F 和ρ＝J‑1
ρ0将其转变到当前构型，即：其中，ω表示位移场任意权函数，然后，通过在物质点上进行离散上式(12)得：式中：Np表示物质点总数，Nn表示活动网格节点总数，Vp表示物质点p的粒子域体积，下标p、I分别表示变量与物质点和背景网格节点相关，φI、 分别表示背景网格广义插值数值基函数及其梯度， 表示可替代的背景网格数值基函数；
由于上式(13)涉及材料非线性和几何非线性，因此，需要先对其采用线性化操作，再构造成Newton‑Raphson增量迭代方程进行求解：KuΔuI＝Ru (14)其中，Ku表示位移场切线刚度矩阵，ΔuI表示背景网格节点增量位移向量，表示位移场节点残差力向量， 和 分别表示外部和内部节点力向量，它们的具体表示形式如下：
其中， 则αp表示材料刚度项， 表示几何刚度项，其中I表示二阶单位张量，下标i、j、k和l均表示张量分析中的自由指标；
然后，对相场演化方程(11)采用Galerkin法构造弱形式方程为：其中，为相场任意权函数，相应的，将上式(18)通过物质点进行离散得：因此，上式(19)构造相场增量迭代求解方程构造如下：KcΔcI＝Rc (20)其中，Kc为相场切线刚度矩阵，ΔcI为背景网格节点增量相场向量，Rc表示相场节点残差力向量，其具体表达形式如下：根据上述理论推导得出的压力敏感岩土材料有限变形弹塑性断裂分析的增量迭代离散形式控制方程，采用交错迭代求解策略，再结合对流粒子域插值技术和隐式物质点方法数值计算方法，即实现岩土结构大变形断裂分析的相场物质点方法，其具体实施步骤如下：步骤一、根据变形体域范围大小划分Euler计算背景网格、根据离散分辨率关系在变形体域内建立离散的Lagrange物质点模型，然后在物质点上定义物理材料属性，再根据给定初始条件初始化域内物质点的物理属性；
步骤二、每一加载步开始时初始化Euler计算背景网格，根据变形后离散物质点模型域大小对Euler计算背景网格划分活动网格及非活动网格、活动网格节点及非活动网格节点，并在离散的物质点和活动网格节点间通过对流粒子域插值技术建立点对点的映射关系，再在交错迭代求解之前通过插值映射函数将物质点携带的物理信息量映射到背景网格节点；
步骤三、通过实施交错迭代求解策略，并在Euler计算背景网格上以物质点作为积分点完成相场和位移场的分步解耦求解，先求解由位移场历史应变能函数 促进裂纹扩展的离散相场演化方程，再求解由相场变量c引起结构刚度退化的离散位移场控制方程，对于位移场控制方程中弹塑性本构关系更新通过引入相场有效应力完成计算；
步骤四、交错迭代求解结束后，将变形后的背景网格节点物理信息量通过插值映射函数映射回物质点，完成当前物质点的位置、变形及破损状态更新，并返回步骤二，进入下一加载步计算循环，直至数值计算完成。
2.根据权利要求1所述的岩土结构大变形断裂分析的相场物质点方法，其特征在于，在步骤二中，对Euler计算背景网格进行活动网格及非活动网格划分主要是为了避免因系统全局刚度矩阵产生奇异而造成数值计算困难，活动网格是指当前时刻下网格内存在一个或多个物质点或粒子域角点即表示网格处于激活状态，反之则为非活动背景网格，而活动网格所包含的网格节点即为活动网格节点。
3.根据权利要求1或2所述的岩土结构大变形断裂分析的相场物质点方法，其特征在于，
在步骤二中，采用对流粒子域插值技术需要先假定如下两个基本假设：每个物质点所在的粒子域始终保持平行四边形，以及变形梯度在粒子域内近似常数；则广义插值数值基函数及其梯度解析表达为：
其中， 和 分别表示当前加载步的粒子域向量 和 的分量，表示物质点p的粒子域角点i关于Euler计算背景网格的插值基函数。
4.根据权利要求1或2所述的岩土结构大变形断裂分析的相场物质点方法，其特征在于，在步骤三中，促进裂纹扩展的位移场历史应变能函数 采用如下方式更新：e
其中，q表示类似应变的塑性内变量，W+表示总应变能函数W(b ,q)中驱动裂纹扩展的能量部分；为了具体化W中的裂纹驱动能量W+，将采用能量分解方法，即先对弹塑性变形的总存储能量采用拉伸‑压缩分解，再进一步进行弹性和塑性部分分解，因此系统自由能函数ψ定义为：
e e e
ψ(b ,q,c)＝g(c)W+(b ,q)+W‑(b ,q) (26)2
其中，g(c)＝(1‑k)c+k表示退化函数，下标+和–分别表示对断裂有贡献和无贡献部分；
在本方法中，采用Hencky弹性应变能量函数描述岩土材料的弹性变形，并对其采用主对数应变拉伸‑压缩分解方法：其中，<·>±＝(·±|·|)/2表示Macaulay括号算子，μ和λ分别表示剪切模量和拉梅常数；因此，主Kirchhoff应力定义为考虑到岩土材料在塑性变形常伴有体积变形特性，则针对塑性能量函数采用体量塑性应变的膨胀‑压缩分解方法：
p,tot
其中，α表示1减Taylor‑Quinney系数，W 表示总的塑性功。
5.根据权利要求1或2所述的岩土结构大变形断裂分析的相场物质点方法，其特征在于，
在步骤三中，为了将相场断裂模型引入光滑双屈服面Drucker‑Prager cap塑性本构模型中，参考损伤力学中损伤有效应力概念，则在弹塑性本构关系更新中引入相场有效Kirchhoff应力张量 为：其中， 表示与相场变量c无关的相场有效Kirchhoff应力张量 的特征值，则其体量应力不变量P和偏量应力不变量Q分别定义为：其中， 表示相场有效Kirchhoff应力张量 的偏量部分，trace(·)表示张量的迹算子；
光滑双屈服面Drucker‑Prager cap塑性本构模型，其包括四个塑性流动区域：①锥顶返回映射区域，②线性非关联返回映射区域，③椭圆非关联返回映射区域，④椭圆关联返回映射区域；两种硬化法则：①线性屈服面的摩擦硬化，②椭圆屈服面的压缩硬化；
定义线性Drucker‑Prager线性屈服面函数为：其中，C表示线性屈服函数与竖向Q轴的截距，M表示线性屈服函数的斜率，与其相关的摩擦硬化法则为：
其中，M0和Mf分别表示斜率M的初值和终值，κ表示摩擦硬化控制参数，γ表示累积塑性乘子；
定义cap椭圆屈服面函数为：其中，Pi为椭圆屈服函数中心，A和B分别表示椭圆屈服函数的短轴和长轴，两个屈服面的相切条件为 切点为 则与椭圆屈服面相关的压缩硬化法则为：*
其中， 为塑性体量对数应变，Pi0为初始椭圆屈服面的中心，ε表示 的终态压缩状态下塑性体量对数应变，r为压缩硬化控制参数；
此外，区域②的线性屈服面和区域③的椭圆屈服面采用非关联流动法则，则定义塑性势能函数为：
其中， 表示与剪胀角相关的材料参数，满足对于区域④的椭圆屈服面采用关联流动法则，即塑性势能函数与屈服函数 一致，则两个塑性势能函数的切点为
压力敏感岩土材料有限变形弹塑性断裂分析中具体的光滑双屈服面Drucker‑Prager cap塑性本构模型的应力返回更新算法包含了四种不同应力返回映射更新过程，其具体实施过程如下：
(1)弹性变形预测，计算测试左Cauchy‑Green张量：eTr
(2)对测试左Cauchy‑Green张量b 进行谱分解：(3)计算测试弹性主对数应变 及其不变量Tr Tr
(4)根据方程(32)和(33)计算相场有效应力及其体量、偏量不变量： P 、Q ；
(5)检查cap椭圆屈服面函数是否满足：(5.1)若 再检查线性屈服面函数是否满足：(5.1.1)若 则进入线性非关联返回映射区域②；
(5.1.1.1)通过线性非关联返回映射算法计算真实主应变不变量 和增量塑性乘子Δγ，并根据方程(35)更新摩擦硬化M；
*
(5.1.1.2)检查真实应力是否在返回映射在区域②内：P≥P*
(5.1.1.2.1)若P≥P，则检查真实应力是否返回映射到互补锥区域①：Q≥0(5.1.1.2.1.1)若Q≥0，则算法正确，进入步骤(8)；
(5.1.1.2.1.2)若Q<0，则采用锥顶返回映射算法，再进入步骤(8)；
*
(5.1.1.2.2)若P
Tr *
(5.1.2)若 则检查相场有效应力体量不变量是否在弹性区域：P ≥PTr *
(5.1.2.1)若P ≥P，进入弹性变形阶段，即进入步骤(7)；
Tr *
(5.1.2.2)若P
(5.2)若 进入弹性变形阶段，即进入步骤(7)；
(6)进入椭圆屈服面，对本构积分内每一个Newton‑Raphson循环迭代步检查：(6.1)若 则真实应力返回映射到椭圆非关联区域③；
(6.1.1)通过椭圆非关联返回映射算法计算真实主应变不变量 和增量塑性乘子Δγ，并根据方程(35)更新摩擦硬化M，并检查塑性体积应变(6.1.1.1)若 则根据方程(37)更新压缩硬化Pi，再进入步骤(8)；
(6.1.1.2)若 则cap椭圆压缩塑性模型失效，进入纯剪切Drucker‑Prager塑性模型；
(6.2)若 则真实应力返回映射到椭圆关联区域④；
(6.2.1)通过椭圆关联返回映射算法计算真实主应变不变量 和增量塑性乘子Δγ，并根据方程(35)更新摩擦硬化M，并检查塑性体积应变(6.2.1.1)若 则根据方程(37)更新压缩硬化Pi，再进入步骤(8)；
(6.2.1.2)若 则cap椭圆压缩塑性模型失效，进入纯剪切Drucker‑prager塑性模型；
Tr Tr
(7)进入弹性变形阶段，设定M不变、Pi＝Pi,n、 (P,Q)＝(P ,Q )；
(8)计算真实弹性主对数应变 和主Kirchhoff应力τA；
ep
(9)计算一致切线算子α 用以更新材料刚度模量αp；
e
(10)根据方程(6)和(7)计算弹性左Cauchy‑Green张量b和Kirchhoff应力张量τ；
在上述实施过程的步骤(1)中 表示相对变形梯度，eTr
为b 的特征值，上标n和n+1分别表示变量属于上一加载步和当前加载步，I表示delta函数，Δu和 分别表示位移增量和插值函数梯度；在步骤(8)中真实弹性主对数应变 通过如下公式修正：
6.根据权利要求1或2所述的岩土结构大变形断裂分析的相场物质点方法，其特征在于，岩土结构大变形断裂分析的相场物质点方法PF‑ICPDI的具体数值实现过程将通过如下伪代码形式展示：
(1)建立离散物质点模型和Euler计算背景网格，定义材料参数，包括弹性：λ,μ；塑性：*
M0,Mf, C,κ,A,Pi0,ε,r；相场：l0, k, 其他：(2)加载步循环：n＝1,2,3,…,Nload(2.1)重划分Euler计算背景网格，并在离散的物质点和活动网格节点间通过对流粒子域插值技术建立联系：φIp,(2.2)映射物质点相场 到网格节点(2.3)迭代求解循环：k＝1,2,3,…,Niter(2.3.1)基于历史应变能函数 计算相场切线刚度矩阵 相场网格节点残差力向量 并求解
(2.3.2)更新网格节点相场 通过映射网格节点相场增量 更新物质点相场(2.3.3)基于Drucker‑Prager cap塑性本构模型和上一加载步物质点相场 计算位移场切线刚度 位移场网格节点残差力向量 并求解(2.3.4)通过映射网格节点位移增量 更新物质点位置 位移(2.3.5)更新物质点变形梯度 和体积(2.3.6)若相场和位移场均满足收敛条件或k>Niter，进入步骤(3)；否则，返回步骤(2.3)；
(3) 计算每个物质点的 裂纹驱动能 量W+用以 更新历史应变能 函数(4)更新物质点的粒子域向量 和(5)若n
说明书 :
岩土结构大变形断裂分析的相场物质点方法
技术领域
背景技术
工程结构的断裂破坏机理研究，若通过实验手段不仅实施困难、代价高昂，而且很难同时考
虑多种外部因素对其断裂破坏的影响。数值仿真方法作为新兴高性能数值计算技术为深入
探究岩土结构的复杂断裂过程提供了一种可行的方案，其在选用合理的本构模型和有效的
损伤法则后即可快速、准确地预测岩土结构的断裂失效行为。因此，发展准确、高效的岩土
结构大变形断裂分析的数值模拟方法尤为重要。
效模式中均产生由脆性断裂向塑性断裂转变的共性。先前研究工作发现，岩土材料的断裂
失效行为随着围压的增大将经历一个转变过程：脆性断裂‑变形带‑耗散塑性流动，即表现
出与压力相关的脆性‑塑性转变断裂行为。此外，岩土材料的脆性断裂行为将随应变加载率
的增大而减弱，则其断裂失效过程中的塑性变形将随应变加载率的增大而增大，表现出与
应变加载率敏感的脆性‑塑性转变断裂行为。本发明将重点考虑围压和应变加载率对岩土
材料弹塑性断裂破坏行为的影响。为了仔细研究压力敏感岩土材料在复杂应力状态和不同
应变加载率下的宏观响应，需要选择一个有效的本构模型去准确捕捉脆性‑塑性失效行为。
目前已发展的岩土材料本构模型主要包括以下三种：考虑剪胀塑性效应的线性屈服面模
型，如Mohr‑Coulomb模型和Drucker‑Prager模型；考虑压缩塑性效应的椭圆屈服面模型，如
剑桥模型和修正的剑桥模型；同时考虑剪胀‑压缩塑性效应的多屈服面模型，如盖帽塑性模
型、非光滑三屈服面塑性模型、光滑三屈服面塑性模型和光滑双屈服面Drucker‑Prager
cap塑性模型。其中，多屈服面模型相较单屈服面模型描述岩土材料的复杂变形行为更加全
面，而光滑双屈服面Drucker‑Prager cap塑性模型相较三屈服面模型在理论描述上更加简
单，数值实施较为容易。因此，本发明将采用光滑双屈服面Drucker‑Prager cap塑性模型开
展压力敏感岩土材料有限变形弹塑性断裂破坏研究。
化计算中的网格依赖性。先前研究工作表明，相场断裂模型是一种处理非局部损伤问题的
高效方法，其通过求解一个偏微分方程即可得出由标量近似表征的裂纹，该模型最大的优
势在于计算复杂断裂问题(如裂纹交叉、分岔及三维空间内裂纹自由扩展)无需额外算法追
踪裂纹结构。相场断裂模型在本世纪初首先由Aranson等人提出，随后Bourdin等人通过变
分法则将其推导并定义与材料断裂能量相互作用，且应用到准静态脆性断裂分析。相场断
裂模型因其具有很多显著优势而得到诸多学者的广泛青睐，随后在断裂力学领域被快速发
展并应用于求解各种复杂断裂问题，例如，准静态断裂问题、动态断裂问题、塑性断裂问题、
多相断裂问题及多物理场断裂问题等。值得注意的是，目前相场断裂模型主要是基于两种
方法推导：其一是基于变分法则，其二是基于微观力平衡法则。二者在求解脆性断裂问题时
无明显差异，但在塑性断裂分析时，前者的理论推导涉及最大塑性耗散法则，适合用于关联
流动金属材料的塑性断裂分析，而非关联流动岩土材料在实验中很难观测到最大塑性耗散
法则，因此，需要谨慎考虑使用基于变分法则推导的相场断裂模型对其进行塑性断裂分析；
后者的理论推导由于不涉及最大塑性耗散法则，对金属和岩土材料的塑性断裂分析均适
用。因此，本发明将考虑使用基于微观力平衡法则推导的相场断裂模型作为损伤函数耦合
光滑双屈服面Drucker‑Prager cap塑性本构模型探究压力敏感岩土材料的弹塑性断裂破
坏问题。
质点法(MPM)和混合有限单元‑有限体积法(FE‑FV)等，而物质点法(MPM)作为近年来被广泛
发展的无网格法之一，其相比有限单元法(FEM)在处理大变形问题时具有较大优势，已被用
于处理各类复杂问题，例如，爆炸、高速冲击、断裂损伤、材料切割、相变及自由液面流动等。
然而，针对岩土材料断裂损伤的数值模拟研究，物质点法(MPM)相较有限单元法(FEM)及其
他无网格法应用较少，且其亟待被发展用以研究岩土材料在复杂外部环境因素影响下的断
裂破坏问题。因此，本发明将提出一种耦合相场对流粒子域插值隐式物质点方法对准静态
压力敏感岩土材料有限变形弹塑性断裂破坏问题展开研究。
发明内容
断裂分析，创新性的提出了一种耦合相场对流粒子域插值隐式物质点方法(coupled
phase‑field implicit material point method with the convected particle domain
interpolation，简称PF‑ICPDI)，其目的在于解决现有技术存在的以下问题：采用相场断裂
模型用以克服传统损伤断裂模型难以准确捕捉裂纹扩展路径的不足；采用光滑双屈服面
Drucker‑Prager cap塑性本构模型用以克服传统单一屈服面塑性本构模型描述岩土材料
复杂力学行为不全面的缺点，以及简化三屈服面塑性本构模型的数值实施复杂度；克服基
于传统有限单元法(FEM)模拟岩土材料有限变形弹塑性断裂破坏时因发生网格畸变而造成
数值精度严重降低的缺点；采用隐式求解格式克服基于显式求解格式计算岩土材料准静态
弹塑性断裂问题出现计算量大、计算精度低的不足；采用对流粒子域插值技术克服传统物
质点法因物质点跨越网格造成严重数值误差的缺陷，并提高施加边界条件的准确度。
cap塑性本构模型分析压力敏感岩土材料的脆性‑塑性转变断裂破坏问题的研究，提供了一
种岩土结构大变形断裂分析的相场物质点方法，具体步骤如下：
边界 可包括位移边界 和外部力边界 满足 通过运动函数
将参考构型Ω0中点X映射到t时刻的当前构型Ωt中点x，则变形梯度可以表示为：
近似的初始裂隙面，c取0表示完全断裂区域，取1表示完整区域，整个系统忽略惯性力效应
并考虑一个准静态条件，则基于虚功原理、微观力平衡法则和热力学第二定律可以推导出
如下强形式控制方程：
滑近似裂纹宽度的正则化模型参数，0
致，仅仅区别在于历史应变场函数 的组分来源不同；此外，虽然两者具有相同表达形式，
但是二者热动力学含义不同，并且本发明推导的相场演化方程可以退化为标准的相场脆性
断裂模型；
T
平衡方程(10)采用Galerkin法在初始构型下构造弱形式控制方程，然后再利用τ＝P·F 和
‑1
ρ＝J ρ0将其转变到当前构型，即：
插值数值基函数及其梯度， 表示可替代的背景网格数值基函数；
它们的具体表示形式如下：
MPM数值计算方法即可实现本发明提出的岩土结构大变形断裂分析的相场物质点方法，结
合图3所示本发明提出的PF‑ICPDI方法计算循环示意图，其具体实施步骤如下：
给定初始条件初始化域内物质点的物理属性；
节点，并在离散的物质点和活动网格节点间通过CPDI插值技术建立点对点的映射关系，再
在交错迭代求解之前通过插值映射函数将物质点携带的物理信息量映射到背景网格节点；
的离散相场演化方程，再求解由相场变量c引起结构刚度退化的离散位移场控制方程，对于
位移场控制方程中弹塑性本构关系更新通过引入相场有效应力完成计算；
下一加载步计算循环，直至数值计算完成；
刻下网格内存在一个或多个物质点或粒子域角点即表示网格处于激活状态，相反则为非活
动网格，而活动网格所包含的网格节点即为活动网格节点；
数及其梯度解析表达为：
变形的总存储能量采用拉伸‑压缩分解，再进一步进行弹性和塑性部分分解，因此系统自由
能函数ψ可以定义为：
Kirchhoff应力张量 为：
域，④椭圆关联返回映射区域；两种硬化法则：①线性屈服面的摩擦硬化，②椭圆屈服面的
压缩硬化；
椭圆屈服面相关的压缩硬化法则为：
应力返回映射更新过程，具体实施过程如下：
Q )；
为b 的特征值，上标n和n+1分别表示变量属于上一加载步和当前加载步，I表示
delta函数，Δu和 分别表示位移增量和插值函数梯度；在步骤(8)中真实弹性主对数应
变 可以通过如下公式修正：
示：
M0,Mf, C,k,A,Pi0,ε,r；相场：l0, k, 其他： )；
体材料断裂领域的应用范围，通过采用相场断裂模型可以有效克服采用传统损伤断裂模型
难以准确捕捉裂纹扩展路径的不足，且采用相场断裂模型可以较为简单地处理复杂裂纹扩
展问题，例如裂纹交叉、分岔及裂纹在三维空间自由扩展等，并且本发明提供的PF‑ICPDI方
法采用了交错迭代求解策略不仅有效提升了计算效率，还降低了数值实施过程的复杂度；
描述岩土材料复杂力学行为不全面的缺陷，并且简化了三屈服面塑性本构模型数值实施的
复杂度，该双流动法则双硬化法则的塑性本构模型使得本发明可以很好地捕捉岩土材料在
复杂外部因素作用下的脆性向塑性转变的断裂失效行为，验证了岩土材料在断裂失效过程
中压力敏感性和应变加载率敏感性，其数值仿真结果与先前已发表工作和实验结果对比吻
合良好；
(FEM)模拟岩土材料有限变形弹塑性断裂破坏时因发生网格畸变而造成数值精度严重降低
的不足，并且采用隐式求解格式相较显式求解格式在计算岩土材料准静态弹塑性断裂问题
时可以显著提高计算效率和计算精度，还通过采用对流粒子域插值技术(CPDI)克服了传统
物质点法因物质点跨越网格造成严重数值误差的缺陷，并提高了施加边界条件的准确度，
此外，本发明提供的算法实施内容可通过简单的变换操作即可与显式物质点方法相适配。
附图说明
云图(e‑h)；
(d‑f)；
＝0.5mm)对弹塑性断裂结果的影响：(a)位移‑约束反力曲线和(b)位移‑名义体积应变曲
线；
断裂能量释放率gc＝0.45、2.25、4.394kN/mm对弹塑性断裂结果的影响：(a)位移‑约束反力
曲线和(b)位移‑名义体积应变曲线；
位移加载率Δu＝10 、5×10 mm、10 mm对弹塑性断裂结果的影响：(a)位移‑约束反力曲线
和(b)位移‑名义体积应变曲线；
具体实施方式
详细说明。
ICPDI方法处理岩土材料有限变形弹塑性断裂的能力。由于本发明采用隐式物质点法求解
相场断裂问题，故对以下所有数值实施例均采用双线性四边形网格单元并在其内部设置2
×2个物质点，且设置相场模型参数l0大于等于背景网格尺寸h来满足相场断裂数值模拟结
果合理性。特别说明，对于第二个验证Drucker‑Prager cap塑性本构模型准确性的数值实
施例，为了方便处理圆孔周边应力集中问题，我们将采用FEM方法进行本构模型有效性和准
确性验证。我们还针对主要参数对压力敏感岩土材料弹塑性断裂的影响做了详细讨论。此
外，对于前三个数值算例不考虑重力效应，并且对所有算例均假设平面应变条件。
×255个四节点矩形单元，且每个网格内设置2×2个物质点。在结构上边界施加增量位移Δ
‑6
u＝10 mm，材料参数如下表1所示。
力曲线与Miehe工作的对比结果，我们可以清楚的看出二者对比结果吻合良好。此外，通过
监测临界断裂状态发现临界断裂位移ucr＝0.0055mm和临界约束反力Fcr＝0.7147kN，其结
果与Miehe等人工作对比仅存在较小的误差。因此，通过此算例有效说明了本发明采用的相
场断裂模型在脆性断裂模拟时的有效应和准确性。
水平钻孔结构尺寸及边界条件示意图如图9(a)所示，其钻孔半径R＝107.95mm，竖向围压σV
＝32.1Mpa，水平围压σH＝9Mpa。为了模拟不同钻孔压力作用下围压应力状态及塑性变形行
为，我们设置四组钻孔压力P＝0、4、6、8MPa。其精细网格划分结果如图9(b)所示，总共离散
7168个双线性四结点矩形单元，且在钻孔周边划分64个单元。为了方便验证对比数值结果，
本实施例中不考虑摩擦硬化，且弹性自由能函数采用如下表达式：
量塑性应变云图。上述两图给出的数值结果与Spiezia等人工作吻合良好，很好地验证了本
发明采用的光滑双屈服面Drucker‑Prager cap塑性本构模型的准确性和有效性。此外，我
们可以看出随着钻孔压力P增加，钻孔周边塑性变形行为逐渐减小，这一特性反映了适当增
加钻孔压力可以有效提高钻孔围岩的稳定性。与此同时，随着钻孔压力增大，孔边体积塑性
应变逐渐由正变负，此现象表明较高的压力作用将促使剪胀机制向压缩机制转变。因此，本
发明采用的光滑双屈服面Drucker‑Prager cap塑性本构模型可以有效的捕捉岩土材料的
压力敏感特性，且该塑性本构模型的准确性也得到了很好的验证。
化区域，且设置圆圈内所有物质点的内聚力强度C为其他区域的98％。试件离散为120×300
个物质点，相应的背景网格离散为边长为h＝0.5mm的70×152个四结点矩形单元，且每个网
格内设置2×2个物质点。设置七组边界围压σc＝5、20、40、60、100、150、200MPa，在结构上方
‑3
边界施加数值向下的位移增量Δu＝10 mm。该实施例采用的类岩石材料参数如下表3所示。
和相场断裂云图。上述数值模拟结果表明，随着围压σc的增大，岩石试件的稳定性逐渐被增
强，从而可以承受较大的竖向位移载荷以至于试件经历较大的塑性变形，进而导致断裂失
效过程逐渐变长，其主要原因是因为较大的围压致使更多的压缩硬化行为在cap椭圆塑性
模型上被激活，并且由于剪胀塑性流动和相场变量的联合作用导致在临界应力状态之后出
现更慢的应变软化行为。此外，我们还可以清晰的看出岩石试件在低围压作用下的断裂行
为相较高围压作用下断裂行为表现的更加脆性，即随着围压逐渐增大，断裂行为由脆性向
塑性转变。然而，试件在高围压σc＝150、200MPa下仍出现应变局部化失效行为，其主要原因
是高围压致使结构内部微孔洞发生坍塌导致结构出现破坏失效现象，即岩土材料在高围压
下产生脆性破坏失效现象。因此，通过本发明提出的PF‑ICPDI方法模拟得到的压力敏感岩
土材料的脆性‑塑性转变断裂失效的数值结果与先前已发表工作的结果对比吻合良好，很
好地验证了本发明提出的PF‑ICPDI方法模拟压力敏感岩土材料弹塑性断裂问题的准确性
和有效性。
采用不同网格尺寸h＝0.25和0.5mm探究本发明提出的PF‑ICPDI方法的网格敏感性，该对比
实验固定裂纹宽度l0＝0.5mm和每个网格内4个物质点，通过图16所示对比结果可以看出浅
色虚线和黑色实线结果吻合良好，从而有效说明了采用本发明提出的PF‑ICPDI方法模拟岩
土材料应变软化行为不存在网格依赖性，该结果也与Choo等人工作结论一致。与此同时，我
们还通过控制网格尺寸h＝0.5mm探究不同裂纹宽度l0＝0.5和1mm对弹塑性断裂的影响。从
图16中可以看出深色虚线比黑色实线展示了更快的断裂过程，并且临界断裂强度也略有降
低，该特性与裂纹宽度对相场脆性断裂的影响一致。
件的断裂失效过程将出现更多的应变硬化和更慢的软化速率从而导致更慢的断裂过程，此
特性表现的趋势刚好与裂纹宽度l0对弹塑性断裂的影响相反。另一方面，我们还模拟了不
同竖向加载率对压力敏感岩土材料弹塑性断裂的影响，从图18可以看出，较大的竖向加载
率将导致塑性变形阶段出现更多的应变硬化行为，从而试件的整个断裂过程将经历较大的
应变变形现象，这一特点与先前已发表工作的结论一致。
个物质点，相应的背景网格离散为边长为h＝0.03m的205×104个四结点矩形单元，且每个
11
网格内设置2×2个物质点。结构右上方的刚性基底的弹性参数为：杨氏模量E＝5×10 Pa和
‑5
泊松比ν＝0.3，并且在其上部从左往右施加一个线性减小的位移载荷：Δu＝‑3.03×10 x+
‑5
5×10 m,(0≤x≤0.99m)。本实施例考虑重力效应，且其材料参数设置如下表4所示。
ICPDI方法模拟得出的图20(a)岩石边坡的横向位移云图和(b)岩石边坡的竖向位移云图将
出现不连续现象。图20(c)和(d)展示了由本发明提出的PF‑ICPDI方法模拟得出的岩石边坡
终态断裂时等效塑性应变云图和相场云图。在边坡失效模拟过程中，顶部与刚性基底接触
的位移载荷边界的左端首先出现塑性变形，则导致应变局部化带和相场断裂演化首先在此
处发生。从图20(a)中可以观察到一个有趣的结果，当应变局部化带从边坡结构的上端边界
产生后直至贯穿斜坡边界后，应变局部化带的右下端相较顶部与刚性基底接触的位移载荷
边界的左端出现更大的水平位移，此现象将导致损伤后的岩石边坡首先从右端斜坡发生坍
塌。
其准确性。其次，通过平面应变域内钻孔围岩应力应变状态模拟有效验证了本发明采用的
光滑双屈服面Drucker‑Prager cap塑性本构模型在捕捉岩土材料的压力敏感性行为的准
确性。最后，通过实施两个经典岩土结构弹塑性断裂失效数值模拟研究很好地验证了本发
明针对压力敏感高孔隙率岩土材料有限变形弹塑性断裂破坏分析所提出的PF‑ICPDI方法
理论及算法的有效性和准确性。虽然本发明展示的实施例是对二维岩土结构弹塑性断裂破
坏问题进行了有效的数值模拟，但其很容易被拓展至三维岩土结构弹塑性断裂破坏问题分
析。因此，本发明提出的岩土结构大变形断裂分析的相场物质点方法(PF‑ICPDI)是一种极
具发展前景的高性能数值算法。
择和描述实施例是为了更好的说明本发明的原理和实际应用，并且使本领域的普通技术人
员能够理解本发明从而设计适于特定用途的带有各种修改的各种实施例。