CN113360992B - 岩土结构大变形断裂分析的相场物质点方法 - Google Patents
岩土结构大变形断裂分析的相场物质点方法 Download PDFInfo
- Publication number
- CN113360992B CN113360992B CN202110727024.7A CN202110727024A CN113360992B CN 113360992 B CN113360992 B CN 113360992B CN 202110727024 A CN202110727024 A CN 202110727024A CN 113360992 B CN113360992 B CN 113360992B
- Authority
- CN
- China
- Prior art keywords
- phase field
- plastic
- deformation
- fracture
- 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.)
- Active
Links
- 239000000463 material Substances 0.000 title claims abstract description 146
- 238000000034 method Methods 0.000 title claims abstract description 105
- 239000002689 soil Substances 0.000 title claims abstract description 58
- 238000004364 calculation method Methods 0.000 claims abstract description 28
- 230000006870 function Effects 0.000 claims description 73
- 238000006073 displacement reaction Methods 0.000 claims description 44
- 238000004458 analytical method Methods 0.000 claims description 42
- 238000013507 mapping Methods 0.000 claims description 36
- 239000002245 particle Substances 0.000 claims description 35
- 238000007906 compression Methods 0.000 claims description 25
- 239000013598 vector Substances 0.000 claims description 24
- 239000000126 substance Substances 0.000 claims description 23
- 230000006835 compression Effects 0.000 claims description 22
- 230000008569 process Effects 0.000 claims description 18
- 238000004422 calculation algorithm Methods 0.000 claims description 17
- 238000012360 testing method Methods 0.000 claims description 17
- 230000005489 elastic deformation Effects 0.000 claims description 14
- 239000007787 solid Substances 0.000 claims description 13
- 230000000694 effects Effects 0.000 claims description 12
- 238000000354 decomposition reaction Methods 0.000 claims description 11
- 238000011160 research Methods 0.000 claims description 11
- 238000005516 engineering process Methods 0.000 claims description 10
- 239000011159 matrix material Substances 0.000 claims description 10
- 238000005381 potential energy Methods 0.000 claims description 6
- 230000015556 catabolic process Effects 0.000 claims description 5
- 238000006731 degradation reaction Methods 0.000 claims description 5
- 210000000746 body region Anatomy 0.000 claims description 4
- 230000005484 gravity Effects 0.000 claims description 4
- 230000001737 promoting effect Effects 0.000 claims description 4
- 230000007704 transition Effects 0.000 claims description 4
- 230000036961 partial effect Effects 0.000 claims description 3
- 229920000034 Plastomer Polymers 0.000 claims description 2
- 230000001133 acceleration Effects 0.000 claims description 2
- 230000000295 complement effect Effects 0.000 claims description 2
- 125000004122 cyclic group Chemical group 0.000 claims description 2
- 238000004146 energy storage Methods 0.000 claims description 2
- VWDWKYIASSYTQR-UHFFFAOYSA-N sodium nitrate Chemical compound [Na+].[O-][N+]([O-])=O VWDWKYIASSYTQR-UHFFFAOYSA-N 0.000 claims description 2
- 238000001228 spectrum Methods 0.000 claims description 2
- 230000001808 coupling effect Effects 0.000 abstract 1
- 206010017076 Fracture Diseases 0.000 description 141
- 208000010392 Bone Fractures Diseases 0.000 description 140
- 239000011435 rock Substances 0.000 description 27
- 230000006399 behavior Effects 0.000 description 24
- 238000004088 simulation Methods 0.000 description 20
- 238000005553 drilling Methods 0.000 description 14
- 230000009471 action Effects 0.000 description 13
- 238000002474 experimental method Methods 0.000 description 12
- 230000007547 defect Effects 0.000 description 9
- 238000010586 diagram Methods 0.000 description 9
- 238000006243 chemical reaction Methods 0.000 description 8
- 230000009466 transformation Effects 0.000 description 6
- 230000035945 sensitivity Effects 0.000 description 5
- 230000008901 benefit Effects 0.000 description 4
- 230000002829 reductive effect Effects 0.000 description 4
- 239000000243 solution Substances 0.000 description 4
- 230000007246 mechanism Effects 0.000 description 3
- 239000000758 substrate Substances 0.000 description 3
- 206010010214 Compression fracture Diseases 0.000 description 2
- 238000000418 atomic force spectrum Methods 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 230000004807 localization Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000005482 strain hardening Methods 0.000 description 2
- 206010053206 Fracture displacement Diseases 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000001747 exhibiting effect Effects 0.000 description 1
- 238000004880 explosion Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 239000007788 liquid Substances 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 229910052751 metal Inorganic materials 0.000 description 1
- 239000007769 metal material Substances 0.000 description 1
- 150000002739 metals Chemical class 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 239000011148 porous material Substances 0.000 description 1
- 239000012088 reference solution Substances 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000000452 restraining effect Effects 0.000 description 1
- 239000011343 solid material Substances 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
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所示,考虑一个任意固体域(ndim∈{1,2,3}表示结构维度)内部存在一个相场c(x,t)∈[0,1]近似的非连续边界Γ,其外部边界可包括位移边界和外部力边界满足通过运动函数将参考构型Ω0中点X映射到t时刻的当前构型Ωt中点x,则变形梯度可以表示为:
对于有限弹塑性变形,变形梯度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应力张量τ表示,其具体的本构表达式可以通过能量存储函数Ψ表示为:
对弹性左Cauchy-Green应变张量be和Kirchhoff应力张量τ进行谱分解可得:
基于框架不变性假设,则系统自由能函数ψ仅依赖物质点X,且弹性主对数应变可定义为:
因此,Ψ(X,be)=ψ(X,εe),则主Kirchhoff应力τA可以表示为:
如图2所示,一个任意变形体受到体积力ρ0G和边界上的外部面力T作用,其中ρ0表示初始固体密度,G表示重力加速度,Γ为固体域Ω内被相场变量c(x,t)∈[0,1]近似的初始裂隙面,c取0表示完全断裂区域,取1表示完整区域,整个系统忽略惯性力效应并考虑一个准静态条件,则基于虚功原理、微观力平衡法则和热力学第二定律可以推导出如下强形式控制方程:
线动量平衡方程
DivP+ρ0G=0 (10)
相场演化方程
并且有P·N=T为边界上外力边界条件,为边界上相场边界条件,其中N表示边界的外法线向量,Div(·)表示散度算子,表示梯度算子,l0表示控制光滑近似裂纹宽度的正则化模型参数,0<k<<1表示避免刚度阵病态的模型参数,表示临界能量释放率,表示历史场应变能函数;
从上式(11)可以看出,基于微观力平衡法则推导的岩土材料弹塑性断裂分析的相场演化方程与广泛使用的基于变分法则推导的脆性断裂相场演化方程在表达形式完全一致,仅仅区别在于历史应变场函数的组分来源不同;此外,虽然两者具有相同表达形式,但是二者热动力学含义不同,并且本发明推导的相场演化方程可以退化为标准的相场脆性断裂模型;
然后,对上述强形式系统动量平衡方程(10)和相场演化方程(11)在物质点上采用隐式求解格式进行离散,并对系统离散控制方程采用交错迭代求解策略;首先,对系统动量平衡方程(10)采用Galerkin法在初始构型下构造弱形式控制方程,然后再利用τ=P·FT和ρ=J-1ρ0将其转变到当前构型,即:
其中,ω表示位移场任意权函数,然后,通过在物质点上进行离散上式(12)可得:
式中:Np表示物质点总数,Nn表示活动网格节点总数,Vp表示物质点p的粒子域体积,下标p、I分别表示变量与物质点和背景网格节点相关,φI、分别表示背景网格广义插值数值基函数及其梯度,表示可替代的背景网格数值基函数;
由于上式(13)涉及材料非线性和几何非线性,因此,需要先对其采用线性化操作,再构造成Newton-Raphson增量迭代方程进行求解:
KuΔuI=Ru (14)
然后,对相场演化方程(11)采用Galerkin法构造弱形式方程为:
因此,上式(19)构造相场增量迭代求解方程构造如下:
KcΔcI=Rc (20)
其中,Kc为相场切线刚度矩阵,ΔcI为背景网格节点增量相场向量,Rc表示相场节点残差力向量,其具体表达形式如下:
根据上述理论推导得出的压力敏感岩土材料有限变形弹塑性断裂分析的增量迭代离散形式控制方程,采用交错迭代求解策略,再结合对流粒子域插值技术(CPDI)和隐式MPM数值计算方法即可实现本发明提出的岩土结构大变形断裂分析的相场物质点方法,结合图3所示本发明提出的PF-ICPDI方法计算循环示意图,其具体实施步骤如下:
步骤一、根据变形体域范围大小划分Euler计算背景网格、根据离散分辨率关系在变形体域内建立离散的Lagrange物质点模型,然后在物质点上定义物理材料属性,再根据给定初始条件初始化域内物质点的物理属性;
步骤二、每一加载步开始时初始化Euler计算背景网格,根据变形后离散物质点模型域大小对Euler计算背景网格划分活动网格及非活动网格、活动网格节点及非活动网格节点,并在离散的物质点和活动网格节点间通过CPDI插值技术建立点对点的映射关系,再在交错迭代求解之前通过插值映射函数将物质点携带的物理信息量映射到背景网格节点;
步骤三、通过实施交错迭代求解策略,并在Euler计算背景网格上以物质点作为积分点完成相场和位移场的分步解耦求解,先求解由位移场历史应变能函数促进裂纹扩展的离散相场演化方程,再求解由相场变量c引起结构刚度退化的离散位移场控制方程,对于位移场控制方程中弹塑性本构关系更新通过引入相场有效应力完成计算;
步骤四、交错迭代求解结束后,将变形后的背景网格节点物理信息量通过插值映射函数映射回物质点,完成当前物质点的位置、变形及破损状态更新,并返回步骤二,进入下一加载步计算循环,直至数值计算完成;
在步骤二中,对Euler计算背景网格进行活动网格及非活动网格划分主要是为了避免因系统全局刚度矩阵产生奇异而造成数值计算困难,如图3所示,活动网格是指当前时刻下网格内存在一个或多个物质点或粒子域角点即表示网格处于激活状态,相反则为非活动网格,而活动网格所包含的网格节点即为活动网格节点;
在步骤二中,采用CPDI插值技术需要先假定如下两个基本假设:每个物质点所在的粒子域始终保持平行四边形,以及变形梯度在粒子域内近似常数;则广义插值数值基函数及其梯度解析表达为:
其中,q表示类似应变的塑性内变量,W+表示总应变能函数W(be,q)中驱动裂纹扩展的能量部分;为了具体化W中的裂纹驱动能量W+,我们将采用能量分解方法,即先对弹塑性变形的总存储能量采用拉伸-压缩分解,再进一步进行弹性和塑性部分分解,因此系统自由能函数ψ可以定义为:
ψ(be,q,c)=g(c)W+(be,q)+W-(be,q) (26)
其中,g(c)=(1-k)c2+k表示退化函数,下标+和–分别表示对断裂有贡献和无贡献部分;
在本发明中,采用Hencky弹性应变能量函数描述岩土材料的弹性变形,并对其采用主对数应变的拉伸-压缩分解方法:
考虑到岩土材料在塑性变形常伴有体积变形特性,则针对塑性能量函数采用体量塑性应变的膨胀-压缩分解方法:
其中,α表示1减Taylor-Quinney系数,Wp,tot表示总的塑性功;
如图4所示光滑双屈服面Drucker-Prager cap塑性本构模型,其主要包括四个塑性流动区域:①锥顶返回映射区域,②线性非关联返回映射区域,③椭圆非关联返回映射区域,④椭圆关联返回映射区域;两种硬化法则:①线性屈服面的摩擦硬化,②椭圆屈服面的压缩硬化;
定义线性Drucker-Prager线性屈服面函数为:
其中,C表示线性屈服函数与竖向Q轴的截距,M表示线性屈服函数的斜率,与其相关的摩擦硬化法则为:
其中,M0和Mf分别表示斜率M的初值和终值,κ表示摩擦硬化控制参数,γ表示累积塑性乘子;
定义cap椭圆屈服面函数为:
此外,区域②的线性屈服面和区域③的椭圆屈服面采用非关联流动法则,则定义塑性势能函数为:
如图5所示为压力敏感岩土材料有限变形弹塑性断裂分析中具体的光滑双屈服面Drucker-Prager cap塑性本构模型的应力返回更新算法流程图,其包含了图4中四种不同应力返回映射更新过程,具体实施过程如下:
(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.1)若PTr≥P*,进入弹性变形阶段,即进入步骤(7);
(5.1.2.2)若PTr<P*,进入步骤(6);
(9)计算一致切线算子αep用以更新材料刚度模量αp;
(10)根据方程(6)和(7)计算弹性左Cauchy-Green张量be和Kirchhoff应力张量τ;
在上述实施过程的步骤(1)中表示相对变形梯度,为beTr的特征值,上标n和n+1分别表示变量属于上一加载步和当前加载步,I表示delta函数,Δu和分别表示位移增量和插值函数梯度;在步骤(8)中真实弹性主对数应变可以通过如下公式修正:
结合图1所示的PF-ICPDI操作流程图和图3所示的计算循环示意图,岩土结构大变形断裂分析的相场物质点方法(PF-ICPDI)的具体数值实现过程将通过如下伪代码形式展示:
(2)加载步循环:n=1,2,3,…,Nload
(2.3)迭代求解循环:k=1,2,3,…,Niter
(2.3.6)若相场和位移场均满足收敛条件或k>Niter,进入步骤(3);否则,返回步骤(2.3);
(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作用下不同临界断裂能量释放率gc=0.45、2.25、4.394kN/mm2对弹塑性断裂结果的影响:(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单开口板拉伸破坏实验材料参数表
图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个单元。为了方便验证对比数值结果,本实施例中不考虑摩擦硬化,且弹性自由能函数采用如下表达式:
表2平面应变域内水平钻孔压缩实验材料参数表
图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平面应变压缩断裂破坏实验类岩石材料参数表
在不同围压σ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中可以看出深色虚线比黑色实线展示了更快的断裂过程,并且临界断裂强度也略有降低,该特性与裂纹宽度对相场脆性断裂的影响一致。
不同临界断裂能量释放率对压力敏感岩土材料有限变形弹塑性断裂影响的数值对比结果如图17所示。从图中我们可以清晰的看出随着临界能量释放率值得增大,试件的断裂失效过程将出现更多的应变硬化和更慢的软化速率从而导致更慢的断裂过程,此特性表现的趋势刚好与裂纹宽度l0对弹塑性断裂的影响相反。另一方面,我们还模拟了不同竖向加载率对压力敏感岩土材料弹塑性断裂的影响,从图18可以看出,较大的竖向加载率将导致塑性变形阶段出现更多的应变硬化行为,从而试件的整个断裂过程将经历较大的应变变形现象,这一特点与先前已发表工作的结论一致。
(4)实施例4:边坡失效模拟(附图19~20)
最后一个实施例将考虑岩石边坡失效模拟来验证本文方法在研究经典岩土结构的局部化变形失效的能力。边坡结构的几何尺寸和边界条件如图19所示,其被离散为60364个物质点,相应的背景网格离散为边长为h=0.03m的205×104个四结点矩形单元,且每个网格内设置2×2个物质点。结构右上方的刚性基底的弹性参数为:杨氏模量E=5×1011Pa和泊松比ν=0.3,并且在其上部从左往右施加一个线性减小的位移载荷:Δu=-3.03×10-5x+5×10-5m,(0≤x≤0.99m)。本实施例考虑重力效应,且其材料参数设置如下表4所示。
表4平面应变压缩断裂破坏实验类岩石材料参数表
因为相场断裂模型在岩土材料弹塑性断裂模拟中可以被认为是一种软化法则,则岩石边坡的失效行为将出现在塑性变形区的应变局部化带中,因此,由本发明提出的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,具体步骤如下:
首先,基于有限变形运动学框架,考虑一个任意固体域内部存在一个相场c(x,t)∈[0,1]近似的非连续边界Γ,其外部边界包括位移边界和外部力边界满足通过运动函数将参考构型Ω0中点X映射到t时刻的当前构型Ωt中点x,ndim∈{1,2,3}表示结构维度;则变形梯度表示为:
对于有限弹塑性变形,变形梯度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应力张量τ表示,其具体的本构表达式通过能量存储函数Ψ表示为:
对弹性左Cauchy-Green应变张量be和Kirchhoff应力张量τ进行谱分解得:
基于框架不变性假设,则系统自由能函数ψ仅依赖物质点X,且弹性主对数应变可定义为:
因此,Ψ(X,be)=ψ(X,εe),则主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<k<<1表示避免刚度阵病态的模型参数,表示临界能量释放率,表示历史场应变能函数;
从上式(11)看出,基于微观力平衡法则推导的岩土材料弹塑性断裂分析的相场演化方程与广泛使用的基于变分法则推导的脆性断裂相场演化方程在表达形式完全一致,仅仅区别在于历史应变场函数的组分来源不同;此外,虽然两者具有相同表达形式,但是二者热动力学含义不同,并且本方法推导的相场演化方程退化为标准的相场脆性断裂模型;
然后,对上述强形式系统动量平衡方程(10)和相场演化方程(11)在物质点上采用隐式求解格式进行离散,并对系统离散控制方程采用交错迭代求解策略;首先,对系统动量平衡方程(10)采用Galerkin法在初始构型下构造弱形式控制方程,然后再利用τ=P·FT和ρ=J-1ρ0将其转变到当前构型,即:
其中,ω表示位移场任意权函数,然后,通过在物质点上进行离散上式(12)得:
式中:Np表示物质点总数,Nn表示活动网格节点总数,Vp表示物质点p的粒子域体积,下标p、I分别表示变量与物质点和背景网格节点相关,φI、分别表示背景网格广义插值数值基函数及其梯度,表示可替代的背景网格数值基函数;
由于上式(13)涉及材料非线性和几何非线性,因此,需要先对其采用线性化操作,再构造成Newton-Raphson增量迭代方程进行求解:
KuΔuI=Ru (14)
然后,对相场演化方程(11)采用Galerkin法构造弱形式方程为:
因此,上式(19)构造相场增量迭代求解方程构造如下:
KcΔcI=Rc (20)
其中,Kc为相场切线刚度矩阵,ΔcI为背景网格节点增量相场向量,Rc表示相场节点残差力向量,其具体表达形式如下:
根据上述理论推导得出的压力敏感岩土材料有限变形弹塑性断裂分析的增量迭代离散形式控制方程,采用交错迭代求解策略,再结合对流粒子域插值技术和隐式物质点方法数值计算方法,即实现岩土结构大变形断裂分析的相场物质点方法,其具体实施步骤如下:
步骤一、根据变形体域范围大小划分Euler计算背景网格、根据离散分辨率关系在变形体域内建立离散的Lagrange物质点模型,然后在物质点上定义物理材料属性,再根据给定初始条件初始化域内物质点的物理属性;
步骤二、每一加载步开始时初始化Euler计算背景网格,根据变形后离散物质点模型域大小对Euler计算背景网格划分活动网格及非活动网格、活动网格节点及非活动网格节点,并在离散的物质点和活动网格节点间通过对流粒子域插值技术建立点对点的映射关系,再在交错迭代求解之前通过插值映射函数将物质点携带的物理信息量映射到背景网格节点;
步骤三、通过实施交错迭代求解策略,并在Euler计算背景网格上以物质点作为积分点完成相场和位移场的分步解耦求解,先求解由位移场历史应变能函数促进裂纹扩展的离散相场演化方程,再求解由相场变量c引起结构刚度退化的离散位移场控制方程,对于位移场控制方程中弹塑性本构关系更新通过引入相场有效应力完成计算;
步骤四、交错迭代求解结束后,将变形后的背景网格节点物理信息量通过插值映射函数映射回物质点,完成当前物质点的位置、变形及破损状态更新,并返回步骤二,进入下一加载步计算循环,直至数值计算完成。
2.根据权利要求1所述的岩土结构大变形断裂分析的相场物质点方法,其特征在于,
在步骤二中,对Euler计算背景网格进行活动网格及非活动网格划分主要是为了避免因系统全局刚度矩阵产生奇异而造成数值计算困难,活动网格是指当前时刻下网格内存在一个或多个物质点或粒子域角点即表示网格处于激活状态,反之则为非活动背景网格,而活动网格所包含的网格节点即为活动网格节点。
其中,q表示类似应变的塑性内变量,W+表示总应变能函数W(be,q)中驱动裂纹扩展的能量部分;为了具体化W中的裂纹驱动能量W+,将采用能量分解方法,即先对弹塑性变形的总存储能量采用拉伸-压缩分解,再进一步进行弹性和塑性部分分解,因此系统自由能函数ψ定义为:
ψ(be,q,c)=g(c)W+(be,q)+W-(be,q) (26)
其中,g(c)=(1-k)c2+k表示退化函数,下标+和–分别表示对断裂有贡献和无贡献部分;
在本方法中,采用Hencky弹性应变能量函数描述岩土材料的弹性变形,并对其采用主对数应变拉伸-压缩分解方法:
考虑到岩土材料在塑性变形常伴有体积变形特性,则针对塑性能量函数采用体量塑性应变的膨胀-压缩分解方法:
其中,α表示1减Taylor-Quinney系数,Wp,tot表示总的塑性功。
5.根据权利要求1或2所述的岩土结构大变形断裂分析的相场物质点方法,其特征在于,
光滑双屈服面Drucker-Prager cap塑性本构模型,其包括四个塑性流动区域:①锥顶返回映射区域,②线性非关联返回映射区域,③椭圆非关联返回映射区域,④椭圆关联返回映射区域;两种硬化法则:①线性屈服面的摩擦硬化,②椭圆屈服面的压缩硬化;
定义线性Drucker-Prager线性屈服面函数为:
其中,C表示线性屈服函数与竖向Q轴的截距,M表示线性屈服函数的斜率,与其相关的摩擦硬化法则为:
其中,M0和Mf分别表示斜率M的初值和终值,κ表示摩擦硬化控制参数,γ表示累积塑性乘子;
定义cap椭圆屈服面函数为:
此外,区域②的线性屈服面和区域③的椭圆屈服面采用非关联流动法则,则定义塑性势能函数为:
压力敏感岩土材料有限变形弹塑性断裂分析中具体的光滑双屈服面Drucker-Pragercap塑性本构模型的应力返回更新算法包含了四种不同应力返回映射更新过程,其具体实施过程如下:
(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.1)若PTr≥P*,进入弹性变形阶段,即进入步骤(7);
(5.1.2.2)若PTr<P*,进入步骤(6);
(9)计算一致切线算子αep用以更新材料刚度模量αp;
(10)根据方程(6)和(7)计算弹性左Cauchy-Green张量be和Kirchhoff应力张量τ;
在上述实施过程的步骤(1)中表示相对变形梯度,为beTr的特征值,上标n和n+1分别表示变量属于上一加载步和当前加载步,I表示delta函数,Δu和分别表示位移增量和插值函数梯度;在步骤(8)中真实弹性主对数应变通过如下公式修正:
6.根据权利要求1或2所述的岩土结构大变形断裂分析的相场物质点方法,其特征在于,岩土结构大变形断裂分析的相场物质点方法PF-ICPDI的具体数值实现过程将通过如下伪代码形式展示:
(2)加载步循环:n=1,2,3,…,Nload
(2.3)迭代求解循环:k=1,2,3,…,Niter
(2.3.6)若相场和位移场均满足收敛条件或k>Niter,进入步骤(3);否则,返回步骤(2.3);
(5)若n<Nload进入步骤(2);否则,计算结束。
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 CN113360992A (zh) | 2021-09-07 |
CN113360992B true 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) |
Families Citing this family (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113820474B (zh) * | 2021-11-24 | 2022-02-15 | 中南大学 | 一种模拟脆性岩石复合型裂纹扩展的相场方法 |
CN114186456B (zh) * | 2021-12-02 | 2022-06-14 | 大连理工大学 | 结构冲击弹塑性断裂分析的时间间断态基近场动力学方法 |
CN114491831B (zh) * | 2021-12-24 | 2023-07-18 | 哈尔滨工业大学 | 一种基于断裂相场法的非均匀材料弥散裂纹j积分方法 |
CN115203847B (zh) * | 2022-07-15 | 2024-05-28 | 燕山大学 | 一种基于mpm的各向异性相场断裂算法的仿真方法 |
CN115410663B (zh) * | 2022-08-16 | 2023-06-02 | 大连理工大学 | 动态冲击/接触弹塑性大变形断裂分析显式相场物质点法 |
CN115618691B (zh) * | 2022-11-10 | 2024-01-26 | 东南大学 | 基于纤维增强复合材料各向异性损伤破裂的相场分析方法 |
CN116358439B (zh) * | 2023-06-02 | 2023-08-25 | 山东省地质科学研究院 | 一种岩石有限应变测量方法、系统、电子设备及存储介质 |
CN116911147B (zh) * | 2023-07-11 | 2024-01-23 | 中国科学院计算机网络信息中心 | 一种基于三维自适应划分的物质点仿真并行方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110298105A (zh) * | 2019-06-26 | 2019-10-01 | 大连理工大学 | 饱和多孔介质大变形分析的ccpdi-impm方法 |
CN110705165A (zh) * | 2019-10-08 | 2020-01-17 | 中国石油大学(华东) | 一种构建岩石材料弹塑性-损伤耦合力学本构模型的方法 |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US8280709B2 (en) * | 2008-10-03 | 2012-10-02 | 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 | 中国石油大学(北京) | 一种微观窜流分析方法及装置 |
-
2021
- 2021-06-29 CN CN202110727024.7A patent/CN113360992B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110298105A (zh) * | 2019-06-26 | 2019-10-01 | 大连理工大学 | 饱和多孔介质大变形分析的ccpdi-impm方法 |
CN110705165A (zh) * | 2019-10-08 | 2020-01-17 | 中国石油大学(华东) | 一种构建岩石材料弹塑性-损伤耦合力学本构模型的方法 |
Non-Patent Citations (5)
Title |
---|
A CONCURRENT MULTISCALE METHOD FOR SIMULATION OF CRACK PROPAGATION;Jingkai Wu等;《Acta Mechanica Solida Sinica》;20150630;第28卷(第3期);235-251 * |
Phase-Field Material Point Method for dynamic brittle fracture with isotropic and anisotropic surface energy;Emmanouil G等;《Computer Methods in Applied Mechanics and Engineering》;20191201;第357卷;1-39 * |
岩体热-水-化-力耦合近场动力学模型及数值模拟研究;王允腾;《中国优秀博硕士学位论文全文数据库(博士)基础科学辑(月刊)》;20201015(第10期);A011-16 * |
耦合热弹塑性动力问题的广义插值物质点法研究;陶俊;《中国优秀博硕士学位论文全文数据库(博士)工程科技Ⅱ辑(月刊)》;20200115(第01期);C039-35 * |
饱和多孔介质动力学大变形分析耦合对流粒子域插值方法;张洪武;《应用力学学报》;20121031;第29卷(第5期);494-501 * |
Also Published As
Publication number | Publication date |
---|---|
CN113360992A (zh) | 2021-09-07 |
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 | |
Zhao et al. | Stabilized material point methods for coupled large deformation and fluid flow in porous materials | |
Bandara et al. | Coupling of soil deformation and pore fluid flow using material point method | |
Lai et al. | Peridynamics simulations of geomaterial fragmentation by impulse loads | |
Pasquariello et al. | A cut-cell finite volume–finite element coupling approach for fluid–structure interaction in compressible flow | |
Popov et al. | SLIM3D: A tool for three-dimensional thermomechanical modeling of lithospheric deformation with elasto-visco-plastic rheology | |
Cervera et al. | Mixed stabilized finite element methods in nonlinear solid mechanics: Part I: Formulation | |
JP5911077B2 (ja) | 空気と水と土骨格の連成計算装置および連成計算方法並びに連成計算プログラム | |
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) | 动态冲击/接触弹塑性大变形断裂分析显式相场物质点法 | |
Fan et al. | A four-way enhanced numerical manifold method for crack propagation and failure analysis of rock slopes | |
Acosta et al. | Development of an implicit contact technique for the material point method | |
Iaconeta et al. | An implicit material point method applied to granular flows | |
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 | |
Areias et al. | Two‐scale method for shear bands: thermal effects and variable bandwidth | |
Chen et al. | Non-local continuum damage model for poro-viscoelastic porous media | |
Cyril et al. | Smooth particle hydrodynamics for the analysis of stresses in soil around Jack-in Pile | |
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 | |
Chandra et al. | Soil-structure interaction simulation of landslides impacting a structure using an implicit material point method | |
Habte | Numerical and constitutive modelling of monotonic and cyclic loading in variably saturated soils |
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 |