CN113360992A - 岩土结构大变形断裂分析的相场物质点方法 - Google Patents

岩土结构大变形断裂分析的相场物质点方法 Download PDF

Info

Publication number
CN113360992A
CN113360992A CN202110727024.7A CN202110727024A CN113360992A CN 113360992 A CN113360992 A CN 113360992A CN 202110727024 A CN202110727024 A CN 202110727024A CN 113360992 A CN113360992 A CN 113360992A
Authority
CN
China
Prior art keywords
phase field
plastic
fracture
deformation
strain
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110727024.7A
Other languages
English (en)
Other versions
CN113360992B (zh
Inventor
胡志强
郑勇刚
刘宇
叶宏飞
张涵博
郑章成
刘振海
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Dalian University of Technology
Original Assignee
Dalian University of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Dalian University of Technology filed Critical Dalian University of Technology
Priority to CN202110727024.7A priority Critical patent/CN113360992B/zh
Publication of CN113360992A publication Critical patent/CN113360992A/zh
Application granted granted Critical
Publication of CN113360992B publication Critical patent/CN113360992B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Architecture (AREA)
  • Civil Engineering (AREA)
  • Structural Engineering (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

本发明属于计算力学领域，提供了一种岩土结构大变形断裂分析的相场物质点方法，该方法考虑了高孔隙率岩土材料在复杂外部因素作用下产生脆性向塑性转变的复杂断裂现象，拓宽了物质点方法在固体材料断裂领域的应用范围。本发明采用相场断裂模型作为损伤函数可以准确、高效的捕捉复杂裂纹扩展路径，采用光滑双屈服面Drucker‑Prager cap塑性本构模型准确全面的描述了压力敏感岩土材料复杂力学行为，并通过引入相场有效应力实现了相场断裂模型与塑性本构模型的耦合作用。此外，本发明还采用CPDI插值方法提高了数值稳定性和边界施加的准确度，并实施交错迭代求解策略提升了计算效率、降低了数值实施复杂度。

Description

岩土结构大变形断裂分析的相场物质点方法
技术领域
本发明属于计算力学领域，具体涉及一种岩土结构大变形断裂分析的相场物质点方法。
背景技术
山体滑坡、桥梁坍塌、地基沉降及大坝溃坝等诸多岩土工程结构的变形过程涉及复杂的断裂破坏问题，例如脆性/塑性断裂破坏、多物理场耦合断裂破坏等。对于大型岩土工程结构的断裂破坏机理研究，若通过实验手段不仅实施困难、代价高昂，而且很难同时考虑多种外部因素对其断裂破坏的影响。数值仿真方法作为新兴高性能数值计算技术为深入探究岩土结构的复杂断裂过程提供了一种可行的方案，其在选用合理的本构模型和有效的损伤法则后即可快速、准确地预测岩土结构的断裂失效行为。因此，发展准确、高效的岩土结构大变形断裂分析的数值模拟方法尤为重要。
研究表明，工程实际中常见的岩土材料如岩石、土壤，由于它们具有复杂的内部微结构和力学特性，将导致其在不同外部因素作用下产生多种断裂失效模式，而这些断裂失效模式中均产生由脆性断裂向塑性断裂转变的共性。先前研究工作发现，岩土材料的断裂失效行为随着围压的增大将经历一个转变过程：脆性断裂-变形带-耗散塑性流动，即表现出与压力相关的脆性-塑性转变断裂行为。此外，岩土材料的脆性断裂行为将随应变加载率的增大而减弱，则其断裂失效过程中的塑性变形将随应变加载率的增大而增大，表现出与应变加载率敏感的脆性-塑性转变断裂行为。本发明将重点考虑围压和应变加载率对岩土材料弹塑性断裂破坏行为的影响。为了仔细研究压力敏感岩土材料在复杂应力状态和不同应变加载率下的宏观响应，需要选择一个有效的本构模型去准确捕捉脆性-塑性失效行为。目前已发展的岩土材料本构模型主要包括以下三种：考虑剪胀塑性效应的线性屈服面模型，如Mohr-Coulomb模型和Drucker-Prager模型；考虑压缩塑性效应的椭圆屈服面模型，如剑桥模型和修正的剑桥模型；同时考虑剪胀-压缩塑性效应的多屈服面模型，如盖帽塑性模型、非光滑三屈服面塑性模型、光滑三屈服面塑性模型和光滑双屈服面Drucker-Pragercap塑性模型。其中，多屈服面模型相较单屈服面模型描述岩土材料的复杂变形行为更加全面，而光滑双屈服面Drucker-Prager cap塑性模型相较三屈服面模型在理论描述上更加简单，数值实施较为容易。因此，本发明将采用光滑双屈服面Drucker-Prager cap塑性模型开展压力敏感岩土材料有限变形弹塑性断裂破坏研究。
针对耦合损伤-塑性模型中损伤函数的选择需要特别谨慎，因为选择合理的非局部损伤函数不仅能在应变软化过程中实现结构刚度退化，还能因其非局部特性消除应变软化计算中的网格依赖性。先前研究工作表明，相场断裂模型是一种处理非局部损伤问题的高效方法，其通过求解一个偏微分方程即可得出由标量近似表征的裂纹，该模型最大的优势在于计算复杂断裂问题(如裂纹交叉、分岔及三维空间内裂纹自由扩展)无需额外算法追踪裂纹结构。相场断裂模型在本世纪初首先由Aranson等人提出，随后Bourdin等人通过变分法则将其推导并定义与材料断裂能量相互作用，且应用到准静态脆性断裂分析。相场断裂模型因其具有很多显著优势而得到诸多学者的广泛青睐，随后在断裂力学领域被快速发展并应用于求解各种复杂断裂问题，例如，准静态断裂问题、动态断裂问题、塑性断裂问题、多相断裂问题及多物理场断裂问题等。值得注意的是，目前相场断裂模型主要是基于两种方法推导：其一是基于变分法则，其二是基于微观力平衡法则。二者在求解脆性断裂问题时无明显差异，但在塑性断裂分析时，前者的理论推导涉及最大塑性耗散法则，适合用于关联流动金属材料的塑性断裂分析，而非关联流动岩土材料在实验中很难观测到最大塑性耗散法则，因此，需要谨慎考虑使用基于变分法则推导的相场断裂模型对其进行塑性断裂分析；后者的理论推导由于不涉及最大塑性耗散法则，对金属和岩土材料的塑性断裂分析均适用。因此，本发明将考虑使用基于微观力平衡法则推导的相场断裂模型作为损伤函数耦合光滑双屈服面Drucker-Prager cap塑性本构模型探究压力敏感岩土材料的弹塑性断裂破坏问题。
目前，相场断裂模型已经基于多种数值计算方法展开各类断裂破坏问题研究，例如，有限单元法(FEM)、扩展有限元法(XFEM)、虚拟单元法(VEM)、不连续伽辽金法(DGM)、物质点法(MPM)和混合有限单元-有限体积法(FE-FV)等，而物质点法(MPM)作为近年来被广泛发展的无网格法之一，其相比有限单元法(FEM)在处理大变形问题时具有较大优势，已被用于处理各类复杂问题，例如，爆炸、高速冲击、断裂损伤、材料切割、相变及自由液面流动等。然而，针对岩土材料断裂损伤的数值模拟研究，物质点法(MPM)相较有限单元法(FEM)及其他无网格法应用较少，且其亟待被发展用以研究岩土材料在复杂外部环境因素影响下的断裂破坏问题。因此，本发明将提出一种耦合相场对流粒子域插值隐式物质点方法对准静态压力敏感岩土材料有限变形弹塑性断裂破坏问题展开研究。
发明内容
本发明要解决的技术问题：本发明基于相场断裂模型和光滑双屈服面Drucker-Prager cap塑性本构模型，以及结合对流粒子域插值技术，针对岩土材料有限变形弹塑性断裂分析，创新性的提出了一种耦合相场对流粒子域插值隐式物质点方法(coupledphase-field implicit material point method with the convected particle domaininterpolation，简称PF-ICPDI)，其目的在于解决现有技术存在的以下问题：采用相场断裂模型用以克服传统损伤断裂模型难以准确捕捉裂纹扩展路径的不足；采用光滑双屈服面Drucker-Prager cap塑性本构模型用以克服传统单一屈服面塑性本构模型描述岩土材料复杂力学行为不全面的缺点，以及简化三屈服面塑性本构模型的数值实施复杂度；克服基于传统有限单元法(FEM)模拟岩土材料有限变形弹塑性断裂破坏时因发生网格畸变而造成数值精度严重降低的缺点；采用隐式求解格式克服基于显式求解格式计算岩土材料准静态弹塑性断裂问题出现计算量大、计算精度低的不足；采用对流粒子域插值技术克服传统物质点法因物质点跨越网格造成严重数值误差的缺陷，并提高施加边界条件的准确度。
本发明的技术方案：
岩土结构大变形断裂分析的相场物质点方法(PF-ICPDI)，
为了准确捕捉高孔隙率岩土材料在大变形情况下的真实力学特性，本发明将基于有限应变连续介质力学框架开展采用耦合相场断裂模型/光滑双屈服面Drucker-Pragercap塑性本构模型分析压力敏感岩土材料的脆性-塑性转变断裂破坏问题的研究，提供了一种岩土结构大变形断裂分析的相场物质点方法，具体步骤如下：
首先，基于有限变形运动学框架，如图2所示，考虑一个任意固体域
Figure BDA0003137883050000031
(ndim∈{1,2,3}表示结构维度)内部存在一个相场c(x,t)∈[0,1]近似的非连续边界Γ，其外部边界
Figure BDA0003137883050000032
可包括位移边界
Figure BDA0003137883050000033
和外部力边界
Figure BDA0003137883050000034
满足
Figure BDA0003137883050000035
通过运动函数
Figure BDA0003137883050000036
将参考构型Ω0中点X映射到t时刻的当前构型Ωt中点x，则变形梯度可以表示为：
Figure BDA0003137883050000037
对于有限弹塑性变形，变形梯度F可通过乘式分解为弹性和塑性部分：
F＝FeFp (2)
其中，上标e和p分别表示变量的弹性和塑性部分，后继公式中相同标记均表示此含义，因此，变形梯度的Jacobian也可分解为弹性和塑性部分：
J＝det(F)＝JeJp (3)
其中，det(·)表示行列式算子，Je＝det(Fe)和Jp＝det(Fp)；
在本发明中，我们基于各向同性假设和有限变形框架描述固体的有限弹塑性变形，因此，固体的弹性变形可以通过对称的弹性左Cauchy-Green应变张量描述，即：
be＝FeFeT (4)
此外，固体域内每一点的应力状态可以通过对称的Kirchhoff应力张量τ表示，其具体的本构表达式可以通过能量存储函数Ψ表示为：
Figure BDA0003137883050000041
对弹性左Cauchy-Green应变张量be和Kirchhoff应力张量τ进行谱分解可得：
Figure BDA0003137883050000042
Figure BDA0003137883050000043
其中，n(A)表示特征向量，
Figure BDA0003137883050000044
和τA分别表示张量b和τ的特征值；
基于框架不变性假设，则系统自由能函数ψ仅依赖物质点X，且弹性主对数应变可定义为：
Figure BDA0003137883050000045
相应的体量弹性对数应变可以被表示为
Figure BDA0003137883050000046
偏量弹性对数应变可以被表示为
Figure BDA0003137883050000047
其中
Figure BDA0003137883050000048
和δ＝(1,1,1)；
因此，Ψ(X,be)＝ψ(X,εe)，则主Kirchhoff应力τA可以表示为：
Figure BDA0003137883050000049
如图2所示，一个任意变形体受到体积力ρ0G和边界
Figure BDA00031378830500000410
上的外部面力T作用，其中ρ0表示初始固体密度，G表示重力加速度，Γ为固体域Ω内被相场变量c(x,t)∈[0,1]近似的初始裂隙面，c取0表示完全断裂区域，取1表示完整区域，整个系统忽略惯性力效应并考虑一个准静态条件，则基于虚功原理、微观力平衡法则和热力学第二定律可以推导出如下强形式控制方程：
线动量平衡方程
DivP+ρ0G＝0 (10)
相场演化方程
Figure BDA00031378830500000411
并且有P·N＝T为边界
Figure BDA00031378830500000412
上外力边界条件，
Figure BDA00031378830500000413
为边界
Figure BDA00031378830500000414
上相场边界条件，其中N表示边界
Figure BDA00031378830500000415
的外法线向量，Div(·)表示散度算子，
Figure BDA00031378830500000416
表示梯度算子，l0表示控制光滑近似裂纹宽度的正则化模型参数，0<k＜＜1表示避免刚度阵病态的模型参数，
Figure BDA00031378830500000417
表示临界能量释放率，
Figure BDA00031378830500000418
表示历史场应变能函数；
从上式(11)可以看出，基于微观力平衡法则推导的岩土材料弹塑性断裂分析的相场演化方程与广泛使用的基于变分法则推导的脆性断裂相场演化方程在表达形式完全一致，仅仅区别在于历史应变场函数
Figure BDA0003137883050000051
的组分来源不同；此外，虽然两者具有相同表达形式，但是二者热动力学含义不同，并且本发明推导的相场演化方程可以退化为标准的相场脆性断裂模型；
然后，对上述强形式系统动量平衡方程(10)和相场演化方程(11)在物质点上采用隐式求解格式进行离散，并对系统离散控制方程采用交错迭代求解策略；首先，对系统动量平衡方程(10)采用Galerkin法在初始构型下构造弱形式控制方程，然后再利用τ＝P·FT和ρ＝J-1ρ0将其转变到当前构型，即：
Figure BDA0003137883050000052
其中，ω表示位移场任意权函数，然后，通过在物质点上进行离散上式(12)可得：
Figure BDA0003137883050000053
式中：Np表示物质点总数，Nn表示活动网格节点总数，Vp表示物质点p的粒子域体积，下标p、I分别表示变量与物质点和背景网格节点相关，φI、
Figure BDA0003137883050000054
分别表示背景网格广义插值数值基函数及其梯度，
Figure BDA0003137883050000055
表示可替代的背景网格数值基函数；
由于上式(13)涉及材料非线性和几何非线性，因此，需要先对其采用线性化操作，再构造成Newton-Raphson增量迭代方程进行求解：
KuΔuI＝Ru (14)
其中，Ku表示位移场切线刚度矩阵，ΔiI表示背景网格节点增量位移向量，
Figure BDA0003137883050000056
Figure BDA0003137883050000057
表示位移场节点残差力向量，
Figure BDA0003137883050000058
和
Figure BDA0003137883050000059
分别表示外部和内部节点力向量，它们的具体表示形式如下：
Figure BDA00031378830500000510
Figure BDA00031378830500000511
Figure BDA0003137883050000061
其中，
Figure BDA0003137883050000062
则αp表示材料刚度项，
Figure BDA0003137883050000063
表示几何刚度项，其中I表示二阶单位张量，下标i、j、k和l均表示张量分析中的自由指标；
然后，对相场演化方程(11)采用Galerkin法构造弱形式方程为：
Figure BDA0003137883050000064
其中，
Figure BDA0003137883050000065
为相场任意权函数，相应的，将上式(18)通过物质点进行离散可得：
Figure BDA0003137883050000066
因此，上式(19)构造相场增量迭代求解方程构造如下：
KcΔcI＝Rc (20)
其中，Kc为相场切线刚度矩阵，ΔcI为背景网格节点增量相场向量，Rc表示相场节点残差力向量，其具体表达形式如下：
Figure BDA0003137883050000067
Figure BDA0003137883050000068
根据上述理论推导得出的压力敏感岩土材料有限变形弹塑性断裂分析的增量迭代离散形式控制方程，采用交错迭代求解策略，再结合对流粒子域插值技术(CPDI)和隐式MPM数值计算方法即可实现本发明提出的岩土结构大变形断裂分析的相场物质点方法，结合图3所示本发明提出的PF-ICPDI方法计算循环示意图，其具体实施步骤如下：
步骤一、根据变形体域范围大小划分Euler计算背景网格、根据离散分辨率关系在变形体域内建立离散的Lagrange物质点模型，然后在物质点上定义物理材料属性，再根据给定初始条件初始化域内物质点的物理属性；
步骤二、每一加载步开始时初始化Euler计算背景网格，根据变形后离散物质点模型域大小对Euler计算背景网格划分活动网格及非活动网格、活动网格节点及非活动网格节点，并在离散的物质点和活动网格节点间通过CPDI插值技术建立点对点的映射关系，再在交错迭代求解之前通过插值映射函数将物质点携带的物理信息量映射到背景网格节点；
步骤三、通过实施交错迭代求解策略，并在Euler计算背景网格上以物质点作为积分点完成相场和位移场的分步解耦求解，先求解由位移场历史应变能函数
Figure BDA00031378830500000710
促进裂纹扩展的离散相场演化方程，再求解由相场变量c引起结构刚度退化的离散位移场控制方程，对于位移场控制方程中弹塑性本构关系更新通过引入相场有效应力完成计算；
步骤四、交错迭代求解结束后，将变形后的背景网格节点物理信息量通过插值映射函数映射回物质点，完成当前物质点的位置、变形及破损状态更新，并返回步骤二，进入下一加载步计算循环，直至数值计算完成；
在步骤二中，对Euler计算背景网格进行活动网格及非活动网格划分主要是为了避免因系统全局刚度矩阵产生奇异而造成数值计算困难，如图3所示，活动网格是指当前时刻下网格内存在一个或多个物质点或粒子域角点即表示网格处于激活状态，相反则为非活动网格，而活动网格所包含的网格节点即为活动网格节点；
在步骤二中，采用CPDI插值技术需要先假定如下两个基本假设：每个物质点所在的粒子域始终保持平行四边形，以及变形梯度在粒子域内近似常数；则广义插值数值基函数及其梯度解析表达为：
Figure BDA0003137883050000071
Figure BDA0003137883050000072
其中，
Figure BDA0003137883050000073
和
Figure BDA0003137883050000074
分别表示当前加载步的粒子域向量
Figure BDA0003137883050000075
和
Figure BDA0003137883050000076
的分量，
Figure BDA0003137883050000077
表示物质点p的粒子域角点i关于Euler计算背景网格的插值基函数；
在步骤三中，促进裂纹扩展的位移场历史应变能函数
Figure BDA0003137883050000078
采用如下方式更新：
Figure BDA0003137883050000079
其中，q表示类似应变的塑性内变量，W+表示总应变能函数W(be,q)中驱动裂纹扩展的能量部分；为了具体化W中的裂纹驱动能量W+，我们将采用能量分解方法，即先对弹塑性变形的总存储能量采用拉伸-压缩分解，再进一步进行弹性和塑性部分分解，因此系统自由能函数ψ可以定义为：
ψ(be,q,c)＝g(c)W+(be,q)+W-(be,q) (26)
Figure BDA0003137883050000081
Figure BDA0003137883050000082
其中，g(c)＝(1-k)c2+k表示退化函数，下标+和-分别表示对断裂有贡献和无贡献部分；
在本发明中，采用Hencky弹性应变能量函数描述岩土材料的弹性变形，并对其采用主对数应变拉伸-压缩分解方法：
Figure BDA0003137883050000083
其中，<■>±＝(■±|■|)/2表示Macaulay括号算子，μ和λ分别表示剪切模量和拉梅常数；因此，主Kirchhoff应力可以别定义为
Figure BDA0003137883050000084
考虑到岩土材料在塑性变形常伴有体积变形特性，则针对塑性能量函数采用体量塑性应变的膨胀-压缩分解方法：
Figure BDA0003137883050000085
Figure BDA0003137883050000086
其中，α表示1减Taylor-Quinney系数，Wp,tot表示总的塑性功；
在步骤三中，为了将相场断裂模型引入光滑双屈服面Drucker-Prager cap塑性本构模型中，参考损伤力学中损伤有效应力概念，则在弹塑性本构关系更新中引入相场有效应力
Figure BDA0003137883050000087
为：
Figure BDA0003137883050000088
其中，
Figure BDA0003137883050000089
表示与相场变量c无关的相场有效Kirchhoff应力张量
Figure BDA00031378830500000810
的特征值，则其体量应力不变量P和偏量应力不变量Q分别定义为：
Figure BDA00031378830500000811
其中，
Figure BDA00031378830500000812
表示相场有效Kirchhoff应力张量
Figure BDA00031378830500000813
的偏量部分，trace(·)表示张量的迹算子；
如图4所示光滑双屈服面Drucker-Prager cap塑性本构模型，其主要包括四个塑性流动区域：①锥顶返回映射区域，②线性非关联返回映射区域，③椭圆非关联返回映射区域，④椭圆关联返回映射区域；两种硬化法则：①线性屈服面的摩擦硬化，②椭圆屈服面的压缩硬化；
定义线性Drucker-Prager线性屈服面函数为：
Figure BDA0003137883050000091
其中，C表示线性屈服函数与竖向Q轴的截距，M表示线性屈服函数的斜率，与其相关的摩擦硬化法则为：
Figure BDA0003137883050000092
其中，M0和Mf分别表示斜率M的初值和终值，κ表示摩擦硬化控制参数，γ表示累积塑性乘子；
定义cap椭圆屈服面函数为：
Figure BDA0003137883050000093
其中，Pi为椭圆屈服函数中心，A和B分别表示椭圆屈服函数的短轴和长轴，两个屈服面的相切条件为
Figure BDA0003137883050000094
切点为
Figure BDA0003137883050000095
则与椭圆屈服面相关的压缩硬化法则为：
Figure BDA0003137883050000096
其中，
Figure BDA0003137883050000097
为塑性体量对数应变，Pi0为初始椭圆屈服面的中心，ε*表示
Figure BDA0003137883050000098
的终态压缩状态下塑性体量对数应变，r为压缩硬化控制参数；
此外，区域②的线性屈服面和区域③的椭圆屈服面采用非关联流动法则，则定义塑性势能函数为：
Figure BDA0003137883050000099
其中，
Figure BDA00031378830500000910
表示与剪胀角相关的材料参数，满足
Figure BDA00031378830500000911
对于区域④的椭圆屈服面采用关联流动法则，即塑性势能函数与屈服函数
Figure BDA00031378830500000912
一致，则两个塑性势能函数的切点为
Figure BDA00031378830500000913
如图5所示为压力敏感岩土材料有限变形弹塑性断裂分析中具体的光滑双屈服面Drucker-Prager cap塑性本构模型的应力返回更新算法流程图，其包含了图4中四种不同应力返回映射更新过程，具体实施过程如下：
(1)弹性变形预测，计算测试左Cauchy-Green张量：
Figure BDA00031378830500000914
(2)对测试左Cauchy-Green张量beTr进行谱分解：
Figure BDA0003137883050000101
(3)计算测试弹性主对数应变
Figure BDA0003137883050000102
及其不变量
Figure BDA0003137883050000103
(4)根据方程(32)和(33)计算相场有效应力及其体量、偏量不变量：
Figure BDA0003137883050000104
PTr、QTr；
(5)检查cap椭圆屈服面函数是否满足：
Figure BDA0003137883050000105
(5.1)若
Figure BDA0003137883050000106
再检查线性屈服面函数是否满足：
Figure BDA0003137883050000107
(5.1.1)若
Figure BDA0003137883050000108
则进入线性非关联返回映射区域②；
(5.1.1.1)通过线性非关联返回映射算法计算真实主应变不变量
Figure BDA0003137883050000109
和增量塑性乘子Δγ，并根据方程(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<P*，则真实应力返回椭圆屈服面，即进入步骤(6)；
(5.1.2)若
Figure BDA00031378830500001010
则检查相场有效应力体量不变量是否在弹性区域：PTr≥P*？
(5.1.2.1)若PTr≥P*，进入弹性变形阶段，即进入步骤(7)；
(5.1.2.2)若PTr<P*，进入步骤(6)；
(5.2)若
Figure BDA00031378830500001011
进入弹性变形阶段，即进入步骤(7)；
(6)进入椭圆屈服面，对本构积分内每一个Newton-Raphson循环迭代步检查：
Figure BDA00031378830500001012
(6.1)若
Figure BDA00031378830500001013
则真实应力返回映射到椭圆非关联区域③；
(6.1.1)通过椭圆非关联返回映射算法计算真实主应变不变量
Figure BDA00031378830500001014
和增量塑性乘子Δγ，并根据方程(35)更新摩擦硬化M，并检查塑性体积应变
Figure BDA00031378830500001015
(6.1.1.1)若
Figure BDA00031378830500001016
则根据方程(37)更新压缩硬化Pi，再进入步骤(8)；
(6.1.1.2)若
Figure BDA00031378830500001017
则cap椭圆压缩塑性模型失效，进入纯剪切Drucker-Prager塑性模型；
(6.2)若
Figure BDA00031378830500001018
则真实应力返回映射到椭圆关联区域④；
(6.2.1)通过椭圆关联返回映射算法计算真实主应变不变量
Figure BDA00031378830500001019
和增量塑性乘子Δγ，并根据方程(35)更新摩擦硬化M，并检查塑性体积应变
Figure BDA00031378830500001020
(6.2.1.1)若
Figure BDA00031378830500001021
则根据方程(37)更新压缩硬化Pi，再进入步骤(8)；
(6.2.1.2)若
Figure BDA00031378830500001022
则cap椭圆压缩塑性模型失效，进入纯剪切Drucker-prager塑性模型；
(7)进入弹性变形阶段，设定M不变、Pi＝Pi,n、
Figure BDA00031378830500001023
(P,Q）＝(PTr,QTr)；
(8)计算真实弹性主对数应变
Figure BDA00031378830500001024
和真实均质主Kirchhoff应力τA；
(9)计算一致切线算子αep用以更新材料刚度模量αp；
(10)根据方程(6)和(7)计算弹性左Cauchy-Green张量be和Kirchhoff应力张量τ；
在上述实施过程的步骤(1)中
Figure BDA0003137883050000111
表示相对变形梯度，
Figure BDA0003137883050000112
为beTr的特征值，上标n和n+1分别表示变量属于上一加载步和当前加载步，δ表示delta函数，Δu和
Figure BDA0003137883050000113
分别表示位移增量和插值函数梯度；在步骤(8)中真实弹性主对数应变
Figure BDA0003137883050000114
可以通过如下公式修正：
Figure BDA0003137883050000115
结合图1所示的PF-ICPDI操作流程图和图3所示的计算循环示意图，本发明提出的岩土结构大变形断裂分析的相场物质点方法(PF-ICPDI)的具体数值实现过程将通过如下伪代码形式展示：
(1)建立离散物质点模型和Euler计算背景网格，定义材料参数(弹性：λ,μ；塑性：M0,Mf,
Figure BDA0003137883050000116
C,κ,A,Pi0,ε*,r；相场：l0,
Figure BDA0003137883050000117
k,
Figure BDA0003137883050000118
其他：
Figure BDA0003137883050000119
)；
(2)加载步循环：n＝1,2,3,…,Nload
(2.1)重划分Euler计算背景网格，并在离散的物质点和活动网格节点间通过CPDI插值技术建立联系：φIp,
Figure BDA00031378830500001110
(2.2)映射物质点相场
Figure BDA00031378830500001111
到网格节点
Figure BDA00031378830500001112
(2.3)迭代求解循环：k＝1,2,3,…,Niter
(2.3.1)基于历史应变能函数
Figure BDA00031378830500001113
计算相场切线刚度矩阵
Figure BDA00031378830500001114
相场网格节点残差力向量
Figure BDA00031378830500001115
并求解
Figure BDA00031378830500001116
(2.3.2)更新网格节点相场
Figure BDA00031378830500001117
通过映射网格节点相场增量
Figure BDA00031378830500001118
更新物质点相场
Figure BDA00031378830500001119
(2.3.3)基于Drucker-Prager cap塑性本构模型和上一加载步物质点相场
Figure BDA00031378830500001120
计算位移场切线刚度
Figure BDA00031378830500001121
位移场网格节点残差力向量
Figure BDA00031378830500001122
并求解
Figure BDA00031378830500001123
(2.3.4)通过映射网格节点位移增量
Figure BDA00031378830500001124
更新物质点位置
Figure BDA00031378830500001125
位移
Figure BDA00031378830500001126
(2.3.5)更新物质点变形梯度
Figure BDA00031378830500001127
和体积
Figure BDA00031378830500001128
(2.3.6)若相场和位移场均满足收敛条件或k>Niter，进入步骤(3)；否则，返回步骤(2.3)；
(3)计算每个物质点的裂纹驱动能量W+用以更新历史应变能函数
Figure BDA00031378830500001129
(4)更新物质点的粒子域向量
Figure BDA00031378830500001130
和
Figure BDA00031378830500001131
(5)若n<Nload进入步骤(2)；否则，计算结束；
本发明的有益效果：
采用本发明提供的技术方案与现有技术相比，具有如下显著效果：
(1)本发明提供的一种岩土结构大变形断裂分析的相场物质点方法，其为岩土材料脆性-塑性断裂破坏分析提供了一种全新的数值仿真方法，并且拓宽了物质点方法在固体材料断裂领域的应用范围，通过采用相场断裂模型可以有效克服采用传统损伤断裂模型难以准确捕捉裂纹扩展路径的不足，且采用相场断裂模型可以较为简单地处理复杂裂纹扩展问题，例如裂纹交叉、分岔及裂纹在三维空间自由扩展等，并且本发明提供的PF-ICPDI方法采用了交错迭代求解策略不仅有效提升了计算效率，还降低了数值实施过程的复杂度；
(2)本发明提供的一种岩土结构大变形断裂分析的相场物质点方法，通过采用光滑双屈服面Drucker-Prager cap塑性本构模型有效克服了传统单一屈服面塑性本构模型描述岩土材料复杂力学行为不全面的缺陷，并且简化了三屈服面塑性本构模型数值实施的复杂度，该双流动法则双硬化法则的塑性本构模型使得本发明可以很好地捕捉岩土材料在复杂外部因素作用下的脆性向塑性转变的断裂失效行为，验证了岩土材料在断裂失效过程中压力敏感性和应变加载率敏感性，其数值仿真结果与先前已发表工作和实验结果对比吻合良好；
(3)本发明提供的一种岩土结构大变形断裂分析的相场物质点方法，采用的隐式物质点法因其在处理大变形及极端变形时具有独特优势有效克服了基于传统有限单元法(FEM)模拟岩土材料有限变形弹塑性断裂破坏时因发生网格畸变而造成数值精度严重降低的不足，并且采用隐式求解格式相较显式求解格式在计算岩土材料准静态弹塑性断裂问题时可以显著提高计算效率和计算精度，还通过采用对流粒子域插值技术(CPDI)克服了传统物质点法因物质点跨越网格造成严重数值误差的缺陷，并提高了施加边界条件的准确度，此外，本发明提供的算法实施内容可通过简单的变换操作即可与显式物质点方法相适配。
附图说明
图1为本发明的一种岩土结构大变形断裂分析的相场物质点方法(PF-ICPDI)的操作流程图；
图2为本发明的内部包含相场近似非连续边界Γ的任意固体域Ω的初始和当前构型示意图；
图3为本发明的PF-ICPDI方法计算循环示意图；
图4为本发明的光滑双屈服面Drucker-Prager cap塑性本构模型；
图5为本发明的光滑双屈服面Drucker-Prager cap塑性本构模型的应力返回更新算法流程图；
图6为本发明脆性断裂分析实施例1单开口板拉伸破坏实验结构及边界条件示意图；
图7为本发明脆性断裂分析实施例1由PF-ICPDI得出的加载位移-约束反力曲线与Miehe工作的对比结果；
图8为本发明脆性断裂分析实施例1由PF-ICPDI得出的结构在加载位移(a)u＝0.005525mm、(b)u＝0.0058mm和(c)u＝0.006mm时的相场断裂云图；
图9为本发明塑性变形分析实施例2(a)四分之一平面应变域内水平钻孔结构及边界条件示意图与(b)FEM精细离散网格；
图10为本发明塑性变形分析实施例2不同钻孔压力作用下水平钻孔周边围岩的径向(上)和环向(下)应力(单位：MPa)云图；
图11为本发明塑性变形分析实施例2不同钻孔压力作用下水平钻孔周边围岩的体量(上)和偏量(下)塑性应变云图；
图12为本发明塑性断裂分析实施例3平面应变压缩实验结构及边界条件示意图；
图13为本发明塑性断裂分析实施例3由PF-ICPDI得出的不同围压作用下(a)位移-约束反力曲线和(b)位移-名义体积应变曲线；
图14为本发明塑性断裂分析实施例3由PF-ICPDI得出的不同围压σc＝5(a,e)、20(b,f)、40(c,g)、60MPa(d,h)作用下结构的终态断裂时等效塑性应变云图(a-d)和相场断裂云图(e-h)；
图15为本发明塑性断裂分析实施例3由PF-ICPDI得出的不同围压σc＝100(a,d)、150(b,e)、200MPa(c,f)作用下结构的终态断裂时等效塑性应变云图(a-c)和相场断裂云图(d-f)；
图16为本发明塑性断裂分析实施例3由PF-ICPDI得出的围压5MPa作用下不同网格宽度h＝0.25、0.5mm(固定每个网格内4个物质点)及不同裂纹宽度参数l0＝0.5、1mm(固定h＝0.5mm)对弹塑性断裂结果的影响：(a)位移-约束反力曲线和(b)位移-名义体积应变曲线；
图17为本发明塑性断裂分析实施例3由PF-ICPDI得出的围压5MPa作用下不同临界断裂能量释放率
Figure BDA0003137883050000131
对弹塑性断裂结果的影响：(a)位移-约束反力曲线和(b)位移-名义体积应变曲线；
图18为本发明塑性断裂分析实施例3由PF-ICPDI得出的围压5MPa作用下不同竖向位移加载率Δu＝10-3、5×10-3mm、10-2mm对弹塑性断裂结果的影响：(a)位移-约束反力曲线和(b)位移-名义体积应变曲线；
图19为本发明边坡失效分析实施例4岩石边坡结构尺寸和边界条件示意图；
图20为本发明边坡失效分析实施例4由PF-ICPDI得出的岩石边坡的终态断裂云图：(a)水平位移(单位：mm)、(b)竖向位移(单位：mm)、(c)等效塑性应变和(d)相场。
具体实施方式
下面结合附图和实施例对本发明的性能做出进一步详细说明。以下实施例用于说明本发明，但不能用来限制本发明的适用范围。
为了使本发明的目的、技术方案和具体实施效果展示更加清晰明了，下面通过四个具体实施例结合附图6～20对本发明提出的PF-ICPDI方法的准确性和有效性做进一步的详细说明。
首先，通过两个标准算例分别验证本发明采用的相场模型和Drucker-Prager cap塑性本构模型的准确性和有效性。然后，通过采用两个经典岩土算例验证本发明提出的PF-ICPDI方法处理岩土材料有限变形弹塑性断裂的能力。由于本发明采用隐式物质点法求解相场断裂问题，故对以下所有数值实施例均采用双线性四边形网格单元并在其内部设置2×2个物质点，且设置相场模型参数l0大于等于背景网格尺寸h来满足相场断裂数值模拟结果合理性。特别说明，对于第二个验证Drucker-Prager cap塑性本构模型准确性的数值实施例，为了方便处理圆孔周边应力集中问题，我们将采用FEM方法进行本构模型有效性和准确性验证。我们还针对主要参数对压力敏感岩土材料弹塑性断裂的影响做了详细讨论。此外，对于前三个数值算例不考虑重力效应，并且对所有算例均假设平面应变条件。
(1)实施例1：单开口板拉伸破坏实验(附图6～8)
实施例1为标准数值算例为单开口板拉伸破坏实验，其结构尺寸和边界条件如图6所示。方形板被离散成249000个物质点，相应的背景网格被离散为边长为h＝0.004mm的255×255个四节点矩形单元，且每个网格内设置2×2个物质点。在结构上边界施加增量位移Δu＝10-6mm，材料参数如下表1所示。
表1单开口板拉伸破坏实验材料参数表
Figure BDA0003137883050000141
图8所示为单开口板拉伸破坏实验在加载位移u＝0.005525mm、u＝0.0058mm和u＝0.006mm时的相场断裂云图，图7给出了由本发明提出的PF-ICPDI得出的加载位移-约束反力曲线与Miehe工作的对比结果，我们可以清楚的看出二者对比结果吻合良好。此外，通过监测临界断裂状态发现临界断裂位移ucr＝0.0055mm和临界约束反力Fcr＝0.7147kN，其结果与Miehe等人工作对比仅存在较小的误差。因此，通过此算例有效说明了本发明采用的相场断裂模型在脆性断裂模拟时的有效应和准确性。
(2)实施例2：平面应变域内水平钻孔周边围岩应力状态及塑性变形模拟(附图9～11)
实施例2为验证光滑双屈服面Drucker-Prager cap塑性本构模型有效性和准确性的平面应变域内水平钻孔周边围岩应力状态及塑性变形行为模拟。四分之一平面应变域内水平钻孔结构尺寸及边界条件示意图如图9(a)所示，其钻孔半径R＝107.95mm，竖向围压σV＝32.1Mpa，水平围压σH＝9Mpa。为了模拟不同钻孔压力作用下围压应力状态及塑性变形行为，我们设置四组钻孔压力P＝0、4、6、8MPa。其精细网格划分结果如图9(b)所示，总共离散7168个双线性四结点矩形单元，且在钻孔周边划分64个单元。为了方便验证对比数值结果，本实施例中不考虑摩擦硬化，且弹性自由能函数采用如下表达式：
Figure BDA0003137883050000151
其中，μ0和θ为常数，K和
Figure BDA0003137883050000152
分别表示压缩模量和平均法向应力P0作用下初始弹性体积应变。具体材料参数如下表2所示。
表2平面应变域内水平钻孔压缩实验材料参数表
Figure BDA0003137883050000153
图10给出了不同钻孔压力P＝0、4、6、8MPa作用下水平钻孔周边围岩的径向和环向应力云图，图11展示了不同钻孔压力P＝0、4、6、8MPa作用下水平钻孔周边围岩的体量和偏量塑性应变云图。上述两图给出的数值结果与Spiezia等人工作吻合良好，很好地验证了本发明采用的光滑双屈服面Drucker-Prager cap塑性本构模型的准确性和有效性。此外，我们可以看出随着钻孔压力P增加，钻孔周边塑性变形行为逐渐减小，这一特性反映了适当增加钻孔压力可以有效提高钻孔围岩的稳定性。与此同时，随着钻孔压力增大，孔边体积塑性应变逐渐由正变负，此现象表明较高的压力作用将促使剪胀机制向压缩机制转变。因此，本发明采用的光滑双屈服面Drucker-Prager cap塑性本构模型可以有效的捕捉岩土材料的压力敏感特性，且该塑性本构模型的准确性也得到了很好的验证。
(3)实施例3：平面应变压缩塑性断裂失效实验(附图12～18)
实施例3为验证岩土材料弹塑性断裂行为的不同围压作用下平面应变压缩实验。图12给出了岩石试件的结构尺寸和边界条件示意图，图中圆圈为引导结构非均匀变形的弱化区域，且设置圆圈内所有物质点的内聚力强度C为其他区域的98％。试件离散为120×300个物质点，相应的背景网格离散为边长为h＝0.5mm的70×152个四结点矩形单元，且每个网格内设置2×2个物质点。设置七组边界围压σc＝5、20、40、60、100、150、200MPa，在结构上方边界施加数值向下的位移增量Δu＝10-3mm。该实施例采用的类岩石材料参数如下表3所示。
表3平面应变压缩断裂破坏实验类岩石材料参数表
Figure BDA0003137883050000161
在不同围压σc＝5、20、40、60、100、150、200MPa作用下的位移-约束反力曲线和位移-名义体积应变曲线如图13所示，图14和15给出了相应的终态断裂时等效塑性应变云图和相场断裂云图。上述数值模拟结果表明，随着围压σc的增大，岩石试件的稳定性逐渐被增强，从而可以承受较大的竖向位移载荷以至于试件经历较大的塑性变形，进而导致断裂失效过程逐渐变长，其主要原因是因为较大的围压致使更多的压缩硬化行为在cap椭圆塑性模型上被激活，并且由于剪胀塑性流动和相场变量的联合作用导致在临界应力状态之后出现更慢的应变软化行为。此外，我们还可以清晰的看出岩石试件在低围压作用下的断裂行为相较高围压作用下断裂行为表现的更加脆性，即随着围压逐渐增大，断裂行为由脆性向塑性转变。然而，试件在高围压σc＝150、200MPa下仍出现应变局部化失效行为，其主要原因是高围压致使结构内部微孔洞发生坍塌导致结构出现破坏失效现象，即岩土材料在高围压下产生脆性破坏失效现象。因此，通过本发明提出的PF-ICPDI方法模拟得到的压力敏感岩土材料的脆性-塑性转变断裂失效的数值结果与先前已发表工作的结果对比吻合良好，很好地验证了本发明提出的PF-ICPDI方法模拟压力敏感岩土材料弹塑性断裂问题的准确性和有效性。
此外，我们还讨论了主要材料参数对岩土材料有限变形弹塑性断裂的影响。为了方便结果比对，所有对比实验模拟均以围压σc＝5MPa的数值结果为基准参考解。首先，我们采用不同网格尺寸h＝0.25和0.5mm探究本发明提出的PF-ICPDI方法的网格敏感性，该对比实验固定裂纹宽度l0＝0.5mm和每个网格内4个物质点，通过图16所示对比结果可以看出浅色虚线和黑色实线结果吻合良好，从而有效说明了采用本发明提出的PF-ICPDI方法模拟岩土材料应变软化行为不存在网格依赖性，该结果也与Choo等人工作结论一致。与此同时，我们还通过控制网格尺寸h＝0.5mm探究不同裂纹宽度l0＝0.5和1mm对弹塑性断裂的影响。从图16中可以看出深色虚线比黑色实线展示了更快的断裂过程，并且临界断裂强度也略有降低，该特性与裂纹宽度对相场脆性断裂的影响一致。
不同临界断裂能量释放率
Figure BDA0003137883050000162
对压力敏感岩土材料有限变形弹塑性断裂影响的数值对比结果如图17所示。从图中我们可以清晰的看出随着临界能量释放率
Figure BDA0003137883050000163
值得增大，试件的断裂失效过程将出现更多的应变硬化和更慢的软化速率从而导致更慢的断裂过程，此特性表现的趋势刚好与裂纹宽度l0对弹塑性断裂的影响相反。另一方面，我们还模拟了不同竖向加载率对压力敏感岩土材料弹塑性断裂的影响，从图18可以看出，较大的竖向加载率将导致塑性变形阶段出现更多的应变硬化行为，从而试件的整个断裂过程将经历较大的应变变形现象，这一特点与先前已发表工作的结论一致。
(4)实施例4：边坡失效模拟(附图19～20)
最后一个实施例将考虑岩石边坡失效模拟来验证本文方法在研究经典岩土结构的局部化变形失效的能力。边坡结构的几何尺寸和边界条件如图19所示，其被离散为60364个物质点，相应的背景网格离散为边长为h＝0.03m的205×104个四结点矩形单元，且每个网格内设置2×2个物质点。结构右上方的刚性基底的弹性参数为：杨氏模量E＝5×1011Pa和泊松比v＝0.3，并且在其上部从左往右施加一个线性减小的位移载荷：Δu＝-3.03×10-5x+5×10-5m,(0≤x≤0.99m)。本实施例考虑重力效应，且其材料参数设置如下表4所示。
表4平面应变压缩断裂破坏实验类岩石材料参数表
Figure BDA0003137883050000171
因为相场断裂模型在岩土材料弹塑性断裂模拟中可以被认为是一种软化法则，则岩石边坡的失效行为将出现在塑性变形区的应变局部化带中，因此，由本发明提出的PF-ICPDI方法模拟得出的图20(a)岩石边坡的横向位移云图和(b)岩石边坡的竖向位移云图将出现不连续现象。图20(c)和(d)展示了由本发明提出的PF-ICPDI方法模拟得出的岩石边坡终态断裂时等效塑性应变云图和相场云图。在边坡失效模拟过程中，顶部与刚性基底接触的位移载荷边界的左端首先出现塑性变形，则导致应变局部化带和相场断裂演化首先在此处发生。从图20(a)中可以观察到一个有趣的结果，当应变局部化带从边坡结构的上端边界产生后直至贯穿斜坡边界后，应变局部化带的右下端相较顶部与刚性基底接触的位移载荷边界的左端出现更大的水平位移，此现象将导致损伤后的岩石边坡首先从右端斜坡发生坍塌。
综上所述，我们首先通过实施经典脆性断裂模拟研究验证了本发明采用基于微观力平衡法则及热力学第二定律推导的相场断裂模型对脆性断裂分析的适用性，并且说明了其准确性。其次，通过平面应变域内钻孔围岩应力应变状态模拟有效验证了本发明采用的光滑双屈服面Drucker-Prager cap塑性本构模型在捕捉岩土材料的压力敏感性行为的准确性。最后，通过实施两个经典岩土结构弹塑性断裂失效数值模拟研究很好地验证了本发明针对压力敏感高孔隙率岩土材料有限变形弹塑性断裂破坏分析所提出的PF-ICPDI方法理论及算法的有效性和准确性。虽然本发明展示的实施例是对二维岩土结构弹塑性断裂破坏问题进行了有效的数值模拟，但其很容易被拓展至三维岩土结构弹塑性断裂破坏问题分析。因此，本发明提出的岩土结构大变形断裂分析的相场物质点方法(PF-ICPDI)是一种极具发展前景的高性能数值算法。
本发明的实施例是为了示例和描述起见而给出的，而并不是无遗漏的或者将本发明限于所公开的形式。很多修改和变化对于本领域的普通技术人员而言是显而易见的。选择和描述实施例是为了更好的说明本发明的原理和实际应用，并且使本领域的普通技术人员能够理解本发明从而设计适于特定用途的带有各种修改的各种实施例。

Claims (6)

1.一种岩土结构大变形断裂分析的相场物质点方法，其特征在于，实施步骤如下：
为了准确捕捉高孔隙率岩土材料在大变形情况下的真实力学特性，将基于有限应变连续介质力学框架开展采用耦合相场断裂模型/光滑双屈服面Drucker-Prager cap塑性本构模型分析压力敏感岩土材料的脆性-塑性转变断裂破坏问题的研究，提供一种岩土结构大变形断裂分析的相场物质点方法PF-ICPDI，具体步骤如下：
首先，基于有限变形运动学框架，考虑一个任意固体域
Figure FDA0003137883040000011
内部存在一个相场c(x，t)∈[0，1]近似的非连续边界Γ，其外部边界
Figure FDA0003137883040000012
包括位移边界
Figure FDA0003137883040000013
和外部力边界
Figure FDA0003137883040000014
满足
Figure FDA0003137883040000015
通过运动函数
Figure FDA0003137883040000016
将参考构型Ω0中点X映射到t时刻的当前构型Ωt中点x，ndim∈{1，2，3}表示结构维度；则变形梯度表示为：
Figure FDA0003137883040000017
对于有限弹塑性变形，变形梯度F通过乘式分解为弹性和塑性部分：
F＝FeFp (2)
其中，上标e和p分别表示变量的弹性和塑性部分，因此，变形梯度的Jacobian分解为弹性和塑性部分：
J＝det(F)＝JeJp (3)
其中，det(·)表示行列式算子，Je＝det(Fe)和Jp＝det(Fp)；
基于各向同性假设和有限变形框架描述固体的有限弹塑性变形，因此，固体的弹性变形通过对称的弹性左Cauchy-Green应变张量描述，即：
be＝FeFeT (4)
此外，固体域内每一点的应力状态通过对称的Kirchhoff应力张量τ表示，其具体的本构表达式通过能量存储函数ψ表示为：
Figure FDA0003137883040000018
对弹性左Cauchy-Green应变张量be和Kirchhoff应力张量τ进行谱分解得：
Figure FDA0003137883040000019
Figure FDA00031378830400000110
其中，n(A)表示特征向量，
Figure FDA00031378830400000111
和τA分别表示张量b和τ的特征值；
基于框架不变性假设，则系统自由能函数ψ仅依赖物质点X，且弹性主对数应变可定义为：
Figure FDA00031378830400000112
相应的体量弹性对数应变被表示为
Figure FDA0003137883040000021
偏量弹性对数应变被表示为
Figure FDA0003137883040000022
其中
Figure FDA0003137883040000023
和δ＝(1，1，1)；
因此，ψ(X，be)＝ψ(X，εe)，则主Kirchhoff应力τA表示为：
Figure FDA0003137883040000024
一个任意变形体受到体积力ρ0G和边界
Figure FDA0003137883040000025
上的外部面力T作用，其中ρ0表示初始固体密度，G表示重力加速度，Γ为固体域Ω内被相场变量c(x，t)∈[0，1]近似的初始裂隙面，c取0表示完全断裂区域，取1表示完整区域，整个系统忽略惯性力效应并考虑一个准静态条件，则基于虚功原理、微观力平衡法则和热力学第二定律推导出如下强形式控制方程：
系统动量平衡方程
DivP+ρ0G＝0 (10)
相场演化方程
Figure FDA0003137883040000026
并且有P·N＝T为边界
Figure FDA0003137883040000027
上外力边界条件，
Figure FDA0003137883040000028
为边界
Figure FDA0003137883040000029
上相场边界条件，其中N表示边界
Figure FDA00031378830400000210
的外法线向量，Div(·)表示散度算子，
Figure FDA00031378830400000211
表示梯度算子，l0表示控制光滑近似裂纹宽度的正则化模型参数，0＜k＜＜1表示避免刚度阵病态的模型参数，
Figure FDA00031378830400000212
表示临界能量释放率，
Figure FDA00031378830400000213
表示历史场应变能函数；
从上式(11)看出，基于微观力平衡法则推导的岩土材料弹塑性断裂分析的相场演化方程与广泛使用的基于变分法则推导的脆性断裂相场演化方程在表达形式完全一致，仅仅区别在于历史应变场函数
Figure FDA00031378830400000214
的组分来源不同；此外，虽然两者具有相同表达形式，但是二者热动力学含义不同，并且本方法推导的相场演化方程退化为标准的相场脆性断裂模型；
然后，对上述强形式系统动量平衡方程(10)和相场演化方程(11)在物质点上采用隐式求解格式进行离散，并对系统离散控制方程采用交错迭代求解策略；首先，对系统动量平衡方程(10)采用Galerkin法在初始构型下构造弱形式控制方程，然后再利用τ＝P·FT和ρ＝J-1ρ0将其转变到当前构型，即：
Figure FDA00031378830400000215
其中，ω表示位移场任意权函数，然后，通过在物质点上进行离散上式(12)得：
Figure FDA0003137883040000031
式中：Np表示物质点总数，Nn表示活动网格节点总数，Vp表示物质点p的粒子域体积，下标p、I分别表示变量与物质点和背景网格节点相关，φI、
Figure FDA0003137883040000032
分别表示背景网格广义插值数值基函数及其梯度，
Figure FDA0003137883040000033
表示可替代的背景网格数值基函数；
由于上式(13)涉及材料非线性和几何非线性，因此，需要先对其采用线性化操作，再构造成Newton-Raphson增量迭代方程进行求解：
KuΔuI＝Ru (14)
其中，Ku表示位移场切线刚度矩阵，ΔuI表示背景网格节点增量位移向量，
Figure FDA0003137883040000034
Figure FDA0003137883040000035
表示位移场节点残差力向量，
Figure FDA0003137883040000036
和
Figure FDA0003137883040000037
分别表示外部和内部节点力向量，它们的具体表示形式如下：
Figure FDA0003137883040000038
Figure FDA0003137883040000039
Figure FDA00031378830400000310
其中，
Figure FDA00031378830400000311
则αp表示材料刚度项，
Figure FDA00031378830400000312
表示几何刚度项，其中I表示二阶单位张量，下标i、j、k和l均表示张量分析中的自由指标；
然后，对相场演化方程(11)采用Galerkin法构造弱形式方程为：
Figure FDA00031378830400000313
其中，
Figure FDA00031378830400000314
为相场任意权函数，相应的，将上式(18)通过物质点进行离散得：
Figure FDA00031378830400000315
因此，上式(19)构造相场增量迭代求解方程构造如下：
KcΔcI＝Rc (20)
其中，Kc为相场切线刚度矩阵，ΔcI为背景网格节点增量相场向量，Rc表示相场节点残差力向量，其具体表达形式如下：
Figure FDA0003137883040000041
Figure FDA0003137883040000042
根据上述理论推导得出的压力敏感岩土材料有限变形弹塑性断裂分析的增量迭代离散形式控制方程，采用交错迭代求解策略，再结合对流粒子域插值技术和隐式物质点方法数值计算方法，即实现岩土结构大变形断裂分析的相场物质点方法，其具体实施步骤如下：
步骤一、根据变形体域范围大小划分Euler计算背景网格、根据离散分辨率关系在变形体域内建立离散的Lagrange物质点模型，然后在物质点上定义物理材料属性，再根据给定初始条件初始化域内物质点的物理属性；
步骤二、每一加载步开始时初始化Euler计算背景网格，根据变形后离散物质点模型域大小对Euler计算背景网格划分活动网格及非活动网格、活动网格节点及非活动网格节点，并在离散的物质点和活动网格节点间通过对流粒子域插值技术建立点对点的映射关系，再在交错迭代求解之前通过插值映射函数将物质点携带的物理信息量映射到背景网格节点；
步骤三、通过实施交错迭代求解策略，并在Euler计算背景网格上以物质点作为积分点完成相场和位移场的分步解耦求解，先求解由位移场历史应变能函数
Figure FDA0003137883040000043
促进裂纹扩展的离散相场演化方程，再求解由相场变量c引起结构刚度退化的离散位移场控制方程，对于位移场控制方程中弹塑性本构关系更新通过引入相场有效应力完成计算；
步骤四、交错迭代求解结束后，将变形后的背景网格节点物理信息量通过插值映射函数映射回物质点，完成当前物质点的位置、变形及破损状态更新，并返回步骤二，进入下一加载步计算循环，直至数值计算完成。
2.根据权利要求1所述的岩土结构大变形断裂分析的相场物质点方法，其特征在于，
在步骤二中，对Euler计算背景网格进行活动网格及非活动网格划分主要是为了避免因系统全局刚度矩阵产生奇异而造成数值计算困难，活动网格是指当前时刻下网格内存在一个或多个物质点或粒子域角点即表示网格处于激活状态，反之则为非活动背景网格，而活动网格所包含的网格节点即为活动网格节点。
3.根据权利要求1或2所述的岩土结构大变形断裂分析的相场物质点方法，其特征在于，
在步骤二中，采用对流粒子域插值技术需要先假定如下两个基本假设：每个物质点所在的粒子域始终保持平行四边形，以及变形梯度在粒子域内近似常数；则广义插值数值基函数及其梯度解析表达为：
Figure FDA0003137883040000051
Figure FDA0003137883040000052
其中，
Figure FDA0003137883040000053
和
Figure FDA0003137883040000054
分别表示当前加载步的粒子域向量
Figure FDA0003137883040000055
和
Figure FDA0003137883040000056
的分量，
Figure FDA0003137883040000057
表示物质点p的粒子域角点i关于Euler计算背景网格的插值基函数。
4.根据权利要求1或2所述的岩土结构大变形断裂分析的相场物质点方法，其特征在于，在步骤三中，促进裂纹扩展的位移场历史应变能函数H采用如下方式更新：
Figure FDA0003137883040000058
其中，q表示类似应变的塑性内变量，W+表示总应变能函数W(be，q)中驱动裂纹扩展的能量部分；为了具体化W中的裂纹驱动能量W+，将采用能量分解方法，即先对弹塑性变形的总存储能量采用拉伸-压缩分解，再进一步进行弹性和塑性部分分解，因此系统自由能函数ψ定义为：
ψ(be，q，c)＝g(c)W+(be，q)+W-(be，q) (26)
Figure FDA00031378830400000510
Figure FDA00031378830400000511
其中，g(c)＝(1-k)c2+k表示退化函数，下标+和-分别表示对断裂有贡献和无贡献部分；
在本方法中，采用Hencky弹性应变能量函数描述岩土材料的弹性变形，并对其采用主对数应变拉伸-压缩分解方法：
Figure FDA00031378830400000512
其中，<■>±＝(■±|■|)/2表示Macaulay括号算子，μ和λ分别表示剪切模量和拉梅常数；因此，主Kirchhoff应力别定义为
Figure FDA00031378830400000513
考虑到岩土材料在塑性变形常伴有体积变形特性，则针对塑性能量函数采用体量塑性应变的膨胀-压缩分解方法：
Figure FDA0003137883040000061
Figure FDA0003137883040000062
其中，α表示1减Taylor-Quinney系数，Wp，tot表示总的塑性功。
5.根据权利要求1或2所述的岩土结构大变形断裂分析的相场物质点方法，其特征在于，
在步骤三中，为了将相场断裂模型引入光滑双屈服面Drucker-Prager cap塑性本构模型中，参考损伤力学中损伤有效应力概念，则在弹塑性本构关系更新中引入相场有效应力
Figure FDA00031378830400000612
为：
Figure FDA0003137883040000063
其中，
Figure FDA0003137883040000064
表示与相场变量c无关的相场有效Kirchhoff应力张量
Figure FDA0003137883040000065
的特征值，则其体量应力不变量P和偏量应力不变量Q分别定义为：
Figure FDA0003137883040000066
其中，
Figure FDA0003137883040000067
表示相场有效Kirchhoff应力张量
Figure FDA0003137883040000068
的偏量部分，trace(·)表示张量的迹算子；
光滑双屈服面Drucker-Prager cap塑性本构模型，其包括四个塑性流动区域：①锥顶返回映射区域，②线性非关联返回映射区域，③椭圆非关联返回映射区域，④椭圆关联返回映射区域；两种硬化法则：①线性屈服面的摩擦硬化，②椭圆屈服面的压缩硬化；
定义线性Drucker-Prager线性屈服面函数为：
Figure FDA0003137883040000069
其中，C表示线性屈服函数与竖向Q轴的截距，M表示线性屈服函数的斜率，与其相关的摩擦硬化法则为：
Figure FDA00031378830400000610
其中，M0和Mf分别表示斜率M的初值和终值，κ表示摩擦硬化控制参数，γ表示累积塑性乘子；
定义cap椭圆屈服面函数为：
Figure FDA00031378830400000611
其中，Pi为椭圆屈服函数中心，A和B分别表示椭圆屈服函数的短轴和长轴，两个屈服面的相切条件为
Figure FDA0003137883040000071
切点为
Figure FDA0003137883040000072
则与椭圆屈服面相关的压缩硬化法则为：
Figure FDA0003137883040000073
其中，
Figure FDA0003137883040000074
为塑性体量对数应变，Pi0为初始椭圆屈服面的中心，ε*表示
Figure FDA0003137883040000075
的终态压缩状态下塑性体量对数应变，r为压缩硬化控制参数；
此外，区域②的线性屈服面和区域③的椭圆屈服面采用非关联流动法则，则定义塑性势能函数为：
Figure FDA0003137883040000076
其中，
Figure FDA0003137883040000077
表示与剪胀角相关的材料参数，满足
Figure FDA0003137883040000078
对于区域④的椭圆屈服面采用关联流动法则，即塑性势能函数与屈服函数
Figure FDA0003137883040000079
一致，则两个塑性势能函数的切点为
Figure FDA00031378830400000710
压力敏感岩土材料有限变形弹塑性断裂分析中具体的光滑双屈服面Drucker-Pragercap塑性本构模型的应力返回更新算法包含了四种不同应力返回映射更新过程，其具体实施过程如下：
(1)弹性变形预测，计算测试左Cauchy-Green张量：
Figure FDA00031378830400000711
(2)对测试左Cauchy-Green张量beTr进行谱分解：
Figure FDA00031378830400000712
(3)计算测试弹性主对数应变
Figure FDA00031378830400000713
及其不变量
Figure FDA00031378830400000714
(4)根据方程(32)和(33)计算相场有效应力及其体量、偏量不变量：
Figure FDA00031378830400000715
PTr、QTr；
(5)检查cap椭圆屈服面函数是否满足：
Figure FDA00031378830400000716
(5.1)若
Figure FDA00031378830400000717
再检查线性屈服面函数是否满足：
Figure FDA00031378830400000718
(5.1.1)若
Figure FDA00031378830400000719
则进入线性非关联返回映射区域②；
(5.1.1.1)通过线性非关联返回映射算法计算真实主应变不变量
Figure FDA00031378830400000720
和增量塑性乘子Δγ，并根据方程(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＜P*，则真实应力返回椭圆屈服面，即进入步骤(6)；
(5.1.2)若
Figure FDA0003137883040000081
则检查相场有效应力体量不变量是否在弹性区域：PTr≥P*
(5.1.2.1)若PTr≥P*，进入弹性变形阶段，即进入步骤(7)；
(5.1.2.2)若PTr＜P*，进入步骤(6)；
(5.2)若
Figure FDA0003137883040000082
进入弹性变形阶段，即进入步骤(7)；
(6)进入椭圆屈服面，对本构积分内每一个Newton-Raphson循环迭代步检查：
Figure FDA0003137883040000083
(6.1)若
Figure FDA0003137883040000084
则真实应力返回映射到椭圆非关联区域③；
(6.1.1)通过椭圆非关联返回映射算法计算真实主应变不变量
Figure FDA0003137883040000085
和增量塑性乘子Δγ，并根据方程(35)更新摩擦硬化M，并检查塑性体积应变
Figure FDA0003137883040000086
(6.1.1.1)若
Figure FDA0003137883040000087
则根据方程(37)更新压缩硬化Pi，再进入步骤(8)；
(6.1.1.2)若
Figure FDA0003137883040000088
则cap椭圆压缩塑性模型失效，进入纯剪切Drucker-Prager塑性模型；
(6.2)若
Figure FDA0003137883040000089
则真实应力返回映射到椭圆关联区域④；
(6.2.1)通过椭圆关联返回映射算法计算真实主应变不变量
Figure FDA00031378830400000810
和增量塑性乘子Δγ，并根据方程(35)更新摩擦硬化M，并检查塑性体积应变
Figure FDA00031378830400000811
(6.2.1.1)若
Figure FDA00031378830400000812
则根据方程(37)更新压缩硬化Pi，再进入步骤(8)；
(6.2.1.2)若
Figure FDA00031378830400000813
则cap椭圆压缩塑性模型失效，进入纯剪切Drucker-prager塑性模型；
(7)进入弹性变形阶段，设定M不变、Pi＝Pi，n、
Figure FDA00031378830400000814
(P，Q)＝(PTr，QTr)；
(8)计算真实弹性主对数应变
Figure FDA00031378830400000815
和真实均质主Kirchhoff应力τA；
(9)计算一致切线算子αep用以更新材料刚度模量αp；
(10)根据方程(6)和(7)计算弹性左Cauchy-Green张量be和Kirchhoff应力张量τ；
在上述实施过程的步骤(1)中
Figure FDA00031378830400000816
表示相对变形梯度，
Figure FDA00031378830400000817
为betr的特征值，上标n和n+1分别表示变量属于上一加载步和当前加载步，δ表示delta函数，Δu和
Figure FDA00031378830400000818
分别表示位移增量和插值函数梯度；在步骤(8)中真实弹性主对数应变
Figure FDA00031378830400000819
通过如下公式修正：
Figure FDA00031378830400000820
6.根据权利要求1或2所述的岩土结构大变形断裂分析的相场物质点方法，其特征在于，
本发明提出的岩土结构大变形断裂分析的相场物质点方法PF-ICPDI的具体数值实现过程将通过如下伪代码形式展示：
(1)建立离散物质点模型和Euler计算背景网格，定义材料参数，包括弹性：λ，μ；塑性：M0，Mf，
Figure FDA00031378830400000821
C，κ，A，Pi0，ε*，r；相场：l0，
Figure FDA00031378830400000822
k，
Figure FDA00031378830400000823
其他：
Figure FDA00031378830400000824
(2)加载步循环：n＝1，2，3，...，Nload
(2.1)重划分Euler计算背景网格，并在离散的物质点和活动网格节点间通过对流粒子域插值技术建立联系：φIp，
Figure FDA0003137883040000091
(2.2)映射物质点相场
Figure FDA0003137883040000092
到网格节点
Figure FDA0003137883040000093
(2.3)迭代求解循环：k＝1，2，3，...，Niter
(2.3.1)基于历史应变能函数
Figure FDA0003137883040000094
计算相场切线刚度矩阵
Figure FDA0003137883040000095
相场网格节点残差力向量
Figure FDA0003137883040000096
并求解
Figure FDA0003137883040000097
(2.3.2)更新网格节点相场
Figure FDA0003137883040000098
通过映射网格节点相场增专
Figure FDA0003137883040000099
更新物质点相场
Figure FDA00031378830400000910
(2.3.3)基于Drucker-Prager cap塑性本构模型和上一加载步物质点相场
Figure FDA00031378830400000911
计算位移场切线刚度
Figure FDA00031378830400000912
位移场网格节点残差力向量
Figure FDA00031378830400000913
并求解
Figure FDA00031378830400000914
(2.3.4)通过映射网格节点位移增量
Figure FDA00031378830400000915
更新物质点位置
Figure FDA00031378830400000916
位移
Figure FDA00031378830400000917
(2.3.5)更新物质点变形梯度
Figure FDA00031378830400000918
和体积
Figure FDA00031378830400000919
(2.3.6)若相场和位移场均满足收敛条件或k＞Niter，进入步骤(3)；否则，返回步骤(2.3)；
(3)计算每个物质点的裂纹驱动能量W+用以更新历史应变能函数
Figure FDA00031378830400000920
(4)更新物质点的粒子域向量
Figure FDA00031378830400000921
和
Figure FDA00031378830400000922
(5)若n＜Nload进入步骤(2)；否则，计算结束。
CN202110727024.7A 2021-06-29 2021-06-29 岩土结构大变形断裂分析的相场物质点方法 Active CN113360992B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110727024.7A CN113360992B (zh) 2021-06-29 2021-06-29 岩土结构大变形断裂分析的相场物质点方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110727024.7A CN113360992B (zh) 2021-06-29 2021-06-29 岩土结构大变形断裂分析的相场物质点方法

Publications (2)

Publication Number Publication Date
CN113360992A true CN113360992A (zh) 2021-09-07
CN113360992B CN113360992B (zh) 2022-02-15

Family

ID=77537130

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110727024.7A Active CN113360992B (zh) 2021-06-29 2021-06-29 岩土结构大变形断裂分析的相场物质点方法

Country Status (1)

Country Link
CN (1) CN113360992B (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113820474A (zh) * 2021-11-24 2021-12-21 中南大学 一种模拟脆性岩石复合型裂纹扩展的相场方法
CN114186456A (zh) * 2021-12-02 2022-03-15 大连理工大学 结构冲击弹塑性断裂分析的时间间断态基近场动力学方法
CN114491831A (zh) * 2021-12-24 2022-05-13 哈尔滨工业大学 一种基于断裂相场法的非均匀材料弥散裂纹j积分方法
CN115203847A (zh) * 2022-07-15 2022-10-18 燕山大学 一种基于mpm的各向异性相场断裂算法的仿真方法
CN115410663A (zh) * 2022-08-16 2022-11-29 大连理工大学 动态冲击/接触弹塑性大变形断裂分析显式相场物质点法
CN115618691A (zh) * 2022-11-10 2023-01-17 东南大学 基于纤维增强复合材料各向异性损伤破裂的相场分析方法
CN116358439A (zh) * 2023-06-02 2023-06-30 山东省地质科学研究院 一种岩石有限应变测量方法、系统、电子设备及存储介质
CN116911147A (zh) * 2023-07-11 2023-10-20 中国科学院计算机网络信息中心 一种基于三维自适应划分的物质点仿真并行方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100088076A1 (en) * 2008-10-03 2010-04-08 Schlumberger Technology Corporation Fully coupled simulation for fluid flow and geomechanical properties in oilfield simulation operations
CN107066680A (zh) * 2017-02-04 2017-08-18 中国石油大学(北京) 一种微观窜流分析方法及装置
US20170282247A1 (en) * 2016-04-01 2017-10-05 Board Of Regents, The University Of Texas System Modeling of nanoparticle agglomeration and powder bed formation in microscale selective laser sintering systems
CN110298105A (zh) * 2019-06-26 2019-10-01 大连理工大学 饱和多孔介质大变形分析的ccpdi-impm方法
CN110705165A (zh) * 2019-10-08 2020-01-17 中国石油大学(华东) 一种构建岩石材料弹塑性-损伤耦合力学本构模型的方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100088076A1 (en) * 2008-10-03 2010-04-08 Schlumberger Technology Corporation Fully coupled simulation for fluid flow and geomechanical properties in oilfield simulation operations
US20170282247A1 (en) * 2016-04-01 2017-10-05 Board Of Regents, The University Of Texas System Modeling of nanoparticle agglomeration and powder bed formation in microscale selective laser sintering systems
CN107066680A (zh) * 2017-02-04 2017-08-18 中国石油大学(北京) 一种微观窜流分析方法及装置
CN110298105A (zh) * 2019-06-26 2019-10-01 大连理工大学 饱和多孔介质大变形分析的ccpdi-impm方法
CN110705165A (zh) * 2019-10-08 2020-01-17 中国石油大学(华东) 一种构建岩石材料弹塑性-损伤耦合力学本构模型的方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
EMMANOUIL G等: "Phase-Field Material Point Method for dynamic brittle fracture with isotropic and anisotropic surface energy", 《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》 *
JINGKAI WU等: "A CONCURRENT MULTISCALE METHOD FOR SIMULATION OF CRACK PROPAGATION", 《ACTA MECHANICA SOLIDA SINICA》 *
张洪武: "饱和多孔介质动力学大变形分析耦合对流粒子域插值方法", 《应用力学学报》 *
王允腾: "岩体热-水-化-力耦合近场动力学模型及数值模拟研究", 《中国优秀博硕士学位论文全文数据库(博士)基础科学辑(月刊)》 *
陶俊: "耦合热弹塑性动力问题的广义插值物质点法研究", 《中国优秀博硕士学位论文全文数据库(博士)工程科技Ⅱ辑(月刊)》 *

Cited By (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113820474B (zh) * 2021-11-24 2022-02-15 中南大学 一种模拟脆性岩石复合型裂纹扩展的相场方法
CN113820474A (zh) * 2021-11-24 2021-12-21 中南大学 一种模拟脆性岩石复合型裂纹扩展的相场方法
CN114186456A (zh) * 2021-12-02 2022-03-15 大连理工大学 结构冲击弹塑性断裂分析的时间间断态基近场动力学方法
CN114186456B (zh) * 2021-12-02 2022-06-14 大连理工大学 结构冲击弹塑性断裂分析的时间间断态基近场动力学方法
CN114491831A (zh) * 2021-12-24 2022-05-13 哈尔滨工业大学 一种基于断裂相场法的非均匀材料弥散裂纹j积分方法
CN115203847B (zh) * 2022-07-15 2024-05-28 燕山大学 一种基于mpm的各向异性相场断裂算法的仿真方法
CN115203847A (zh) * 2022-07-15 2022-10-18 燕山大学 一种基于mpm的各向异性相场断裂算法的仿真方法
CN115410663A (zh) * 2022-08-16 2022-11-29 大连理工大学 动态冲击/接触弹塑性大变形断裂分析显式相场物质点法
CN115618691B (zh) * 2022-11-10 2024-01-26 东南大学 基于纤维增强复合材料各向异性损伤破裂的相场分析方法
CN115618691A (zh) * 2022-11-10 2023-01-17 东南大学 基于纤维增强复合材料各向异性损伤破裂的相场分析方法
CN116358439A (zh) * 2023-06-02 2023-06-30 山东省地质科学研究院 一种岩石有限应变测量方法、系统、电子设备及存储介质
CN116358439B (zh) * 2023-06-02 2023-08-25 山东省地质科学研究院 一种岩石有限应变测量方法、系统、电子设备及存储介质
CN116911147A (zh) * 2023-07-11 2023-10-20 中国科学院计算机网络信息中心 一种基于三维自适应划分的物质点仿真并行方法
CN116911147B (zh) * 2023-07-11 2024-01-23 中国科学院计算机网络信息中心 一种基于三维自适应划分的物质点仿真并行方法

Also Published As

Publication number Publication date
CN113360992B (zh) 2022-02-15

Similar Documents

Publication Publication Date Title
CN113360992B (zh) 岩土结构大变形断裂分析的相场物质点方法
Fan et al. A hybrid peridynamics–SPH simulation of soil fragmentation by blast loads of buried explosive
Lai et al. Peridynamics simulations of geomaterial fragmentation by impulse loads
Bandara et al. Coupling of soil deformation and pore fluid flow using material point method
Pasquariello et al. A cut-cell finite volume–finite element coupling approach for fluid–structure interaction in compressible flow
Farhat et al. Robust and provably second‐order explicit–explicit and implicit–explicit staggered time‐integrators for highly non‐linear compressible fluid–structure interaction problems
Nguyen et al. On a family of convected particle domain interpolations in the material point method
Cervera et al. Mixed stabilized finite element methods in nonlinear solid mechanics: Part I: Formulation
Saksala et al. Combined continuum damage‐embedded discontinuity model for explicit dynamic fracture analyses of quasi‐brittle materials
Zhao et al. A generic approach to modelling flexible confined boundary conditions in SPH and its application
Chen et al. Dynamic fracture analysis of the soil-structure interaction system using the scaled boundary finite element method
Zhang et al. A coupling peridynamic approach for the consolidation and dynamic analysis of saturated porous media
CN115410663B (zh) 动态冲击/接触弹塑性大变形断裂分析显式相场物质点法
Adam et al. Thermomechanical modeling of metals at finite strains: First and mixed order finite elements
Fan et al. A four-way enhanced numerical manifold method for crack propagation and failure analysis of rock slopes
Menon et al. Updated Lagrangian unsaturated periporomechanics for extreme large deformation in unsaturated porous media
Tian et al. A challenging dam structural analysis: large-scale implicit thermo-mechanical coupled contact simulation on Tianhe-II
CN112949065A (zh) 模拟层状岩体力学行为的双尺度方法、装置、存储介质及设备
Hashimoto et al. Improvement of discontinuous deformation analysis incorporating implicit updating scheme of friction and joint strength degradation
Cyril et al. Smooth particle hydrodynamics for the analysis of stresses in soil around Jack-in Pile
Chen et al. Non-local continuum damage model for poro-viscoelastic porous media
Sun et al. An improved quadrature scheme in B-spline material point method for large-deformation problem analysis
Liu et al. A general finite deformation hypoelastic-plasticity non-ordinary state-based peridynamics model and its applications
Junior et al. Brittle fracture and hydroelastic simulations based on moving particle simulation
Bakroon et al. Geotechnical large deformation numerical analysis using implicit and explicit integration

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant