CN115410663B - 动态冲击/接触弹塑性大变形断裂分析显式相场物质点法 - Google Patents

动态冲击/接触弹塑性大变形断裂分析显式相场物质点法 Download PDF

Info

Publication number
CN115410663B
CN115410663B CN202210978276.1A CN202210978276A CN115410663B CN 115410663 B CN115410663 B CN 115410663B CN 202210978276 A CN202210978276 A CN 202210978276A CN 115410663 B CN115410663 B CN 115410663B
Authority
CN
China
Prior art keywords
phase field
fracture
vector
plastic
field
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
Application number
CN202210978276.1A
Other languages
English (en)
Other versions
CN115410663A (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
63921 Troops of PLA
Original Assignee
Dalian University of Technology
63921 Troops of PLA
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, 63921 Troops of PLA filed Critical Dalian University of Technology
Priority to CN202210978276.1A priority Critical patent/CN115410663B/zh
Publication of CN115410663A publication Critical patent/CN115410663A/zh
Application granted granted Critical
Publication of CN115410663B publication Critical patent/CN115410663B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C60/00—Computational materials science, i.e. ICT specially adapted for investigating the physical or chemical properties of materials or phenomena associated with their design, synthesis, processing, characterisation or utilisation
    • 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
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
    • 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
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Computing Systems (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

本发明属于岩土结构断裂分析技术领域，提出了一种动态冲击/接触弹塑性大变形断裂分析显式相场物质点法，为岩土材料动态断裂破坏研究提供了一种全新的数值计算方法。在该方法中，提出了一种基于微观力平衡法则推导的显式相场断裂模型，其既可用于分析脆性断裂问题，也可用于求解弹塑性断裂失效问题。还发展了相应的耦合显式相场‑塑性模型，可以有效预测岩土材料的复杂脆性‑塑性断裂失效行为，其相较传统耦合损伤塑性本构模型，数值实施简单、计算效率高。此外，该方法通过搭载显式物质点法，并采用相场‑位移场交错求解策略和粒子接触算法，能够稳定、高效地求解接触和大变形等强非线性大规模断裂破坏问题。

Description

动态冲击/接触弹塑性大变形断裂分析显式相场物质点法
技术领域
本发明涉及岩土结构断裂分析技术领域，尤其涉及一种动态冲击/接触弹塑性大变形断裂分析显式相场物质点法。
背景技术
岩土材料(如岩石、岩土)被广泛应用于土木工程和采矿工程等领域，其因复杂的内部微结构而具有复杂的力学特性。研究岩土材料的失效机理对诸多工程实际应用具有重要意义，例如页岩气及石油开采、隧道开挖和地下存储二氧化碳等。目前，数值仿真作为一种极具发展前景的新型技术可以辅助我们更深入的探究岩土材料的失效过程。为了全面准确地预测材料的失效行为，在断裂数值仿真中选择合适的本构模型和损伤模型至关重要。然而，对于实施稳定、准确的岩土材料动态冲击/接触弹塑性大变形断裂模拟，目前仍是一个极具挑战性的课题，主要因为其需要全面细致地考虑稳定的裂纹扩展模拟与材料、几何及接触非线性的有效融合。
对于断裂模拟，目前已经发展了诸多处理裂纹扩展问题的数值方法，主要可以分为两类：离散裂缝表征方法和弥散裂缝表征方法。首先，对于离散裂缝表征方法，主要有基于格里菲斯理论发展的多种基于有限单元网格类方法或模型，如单元删除法、虚拟裂纹闭合技术(VCCT)、内聚力模型(CZM)、扩展有限单元法(XFEM)和裂纹驱动构形力法。这些方法通常存在如下一个或多个限制缺陷：需要引入额外裂纹追踪算法捕捉裂纹扩展、分叉及交汇或需要引入额外自由度等。其次，对于弥散裂缝表征方法或模型，主要有近场动力学模型(PD)、相场断裂模型和裂纹粒子法(CPM)等。值得一提的是弥散裂纹表征方法相比离散裂纹表征方法具有一个显著优势是通过场变量表征裂缝而不需要引入额外的裂纹追踪算法，一定程度上提升了计算效率和降低了实施难度。此外，在上述弥散裂缝表征方法中，相场断裂模型可以很好的考虑断裂分析与材料、几何及接触非线性的融合，并已被广泛应用于分析诸多复杂断裂破坏问题。
在相场断裂模型中，裂纹面的演化可以通过求解偏微分方程得出标量相场进行表征。最初，针对准静态脆性断裂分析，Francfort和Marigo提出了基于变分法则推导的相场断裂模型。然后，Bourdin等人发展了正则化近似的变分相场断裂模型的数值实施流程，并且定义了相场与材料应变能相互作用。他们的工作实现了相场与位移场的耦合作用，可以通过求解耦合场控制方程自然的得出脆性断裂分析中裂纹扩展路径。随后，Miehe等人从梯度损伤模型概念中解释了相场变量含义，并区分了拉伸、压缩变形对驱动裂纹扩展的不同贡献作用。近年来，相场断裂模型被广泛应用于处理各种复杂断裂破坏问题，例如脆性断裂、塑性断裂、多相断裂以及多物理场断裂问题。需要注意的是，目前相场断裂模型的推导方式主要分为两类：基于变分法则推导和基于微观力平衡法则推导。两种推导方式的最主要区别在于，针对塑性断裂分析，基于微观力平衡法则推导的相场断裂模型不涉及最大塑性耗散法则。先前研究工作也已表明了基于微观力平衡法则推导的相场断裂模型更适合用于分析岩土材料的塑性断裂破坏问题。此外，研究学者目前普遍采用隐式相场公式来分析各类断裂问题，但是隐式相场公式在分析动态断裂问题时明显效率低下，特别是考虑大规模的复杂裂纹扩展问题。最近，Wang等人基于Miehe等人工作发展了可用于分析类岩石材料的准静态和动态断裂失效行为的显式相场方法。Ren等人提出了针对动态断裂分析的子步迭代格式的显式相场断裂模型。Rad和Shen发展了基于GPU并行计算的显式相场方法可用于分析动态脆性断裂问题。上述工作所发展的显式相场断裂模型均是基于变分法则推导，因此，针对岩土材料动态弹塑性断裂问题亟需发展一种基于微观力平衡法则推导的显式相场断裂模型。
在弹塑性断裂模拟中，需要采用合适的本构模型去准确刻画岩土材料的复杂力学行为。而对于岩土材料最常用的塑性本构模型主要有单屈服面塑性模型和多屈服面塑性模型，其中，单屈服面塑性模型可根据描述材料变形特性不同分为以下两种：一是描述剪胀塑性变形，例如Mohr-Coulomb本构模型和Drucker-Prager本构模型；二是描述压缩塑性变形，例如Cam-Clay本构模型和修正的Cam-Clay本构模型。其次，对于多屈服面塑性模型，其结合了上述两种单屈服面塑性模型的各自优势可以同时表征剪胀和压缩塑性变形特性，常用的模型主要有光滑双屈服面塑性模型和三屈服面塑性模型。相关数值和实验研究表明岩土材料在失效过程中通常表现出一种典型的脆性-塑性失效转变行为，而采用多屈服面塑性模型可以更加全面地描述岩土材料弹塑性断裂分析中复杂力学特性。因此，Choo和Sun发展了一种耦合相场-塑性模型用于分析岩土材料脆性-塑性失效转变行为。我们课题组提供了一种隐式相场物质点法用于模拟岩土材料准静态脆性-塑性失效转变过程。
岩土材料断裂破坏行为往往涉及大变形和接触等变形，因此在其高效数值分析方法构造方面需要有效考虑这些变形行为。近年来，物质点法(MPM)被广泛应用于求解各种接触和大变形问题，因其使用拉格朗日物质点结合欧拉背景计算网格的描述方式可以很好地克服大变形模拟中网格畸变问题，并得到业界一致肯定。在MPM中，物理量从物质点映射到背景计算网格用以求解控制方程，变形场在背景网格计算完成后再映射回物质点用以更新物理量。MPM因能稳定可靠的模拟大变形现象而成功吸引了诸多学者基于此开展各种强非线性问题研究，例如，断裂、高速冲击、爆炸等问题。为了降低传统MPM中因物质点跨越网格产生的数值误差，多种改进的插值方式被相继提出，如对流粒子域插值技术(CPDI)、B样条插值物质点法(BS-MPM)和广义插值物质点法(GIMP)。值得注意的是，在MPM中使用CPDI可以更加便捷、准确的施加位移、围压等边界条件。
因此，本发明将提出一种针对岩土材料动态冲击/接触弹塑性大变形断裂问题分析的显式相场物质点法。其中，主要提出了一种基于微观力平衡法则和热力学第二定律推导的显式相场断裂模型，并发展了一种耦合显式相场-塑性模型用以捕捉压力敏感岩土材料的复杂失效行为。此外，还将采用粒子接触算法将此方法拓展至可以处理多体接触弹塑性动态失效问题。
发明内容
针对岩土材料动态冲击/接触弹塑性断裂破坏高效数值分析，本发明创新性的提出了一种动态冲击/接触弹塑性大变形断裂分析显式相场物质点法(英文简称：ePF-CPDI)，其目的在于：首先，为了克服传统损伤模型难以准确捕捉裂纹扩展路径的不足，以及解决隐式相场断裂模型分析动态断裂及大规模裂纹扩展问题的低效性和基于变分法则推导的显式相场断裂模型对岩土材料失效行为分析的理论不适用性，本发明基于微观力平衡法则和热力学第二定律提出适用于岩土材料断裂破坏分析的显式率相关相场断裂模型；其次，为了克服传统单屈服面本构模型难以准确有效的描述压力敏感岩土材料的复杂失效行为和简化三屈服面本构模型的操作复杂度，本发明采用光滑双屈服面本构模型结合提出的显式相场断裂模型发展了一种耦合显式相场-塑性断裂模型；此外，本发明还采用对流粒子域插值技术(CPDI)用以消除传统物质点法中因物质点跨越网格产生的数值噪声，并采用粒子接触算法实现对冲击/接触断裂强非线性问题分析；最后，本发明旨在解决现有相场有限元法在分析大变形断裂破坏问题时因网格畸变造成的数值精度严重下降的缺点和现有隐式相场物质点法在分析动态断裂破坏和大规模裂纹扩展问题时计算效率低下的缺陷。
本发明的技术方案：一种动态冲击/接触弹塑性大变形断裂分析显式相场物质点法，具体步骤如下：
步骤一，定义欧拉背景网格，建立离散物质点模型并定义离散物质点的物理材料参数，初始化物质点变量；物理材料参数包括弹性参数、塑性参数和相场参数；弹性参数包括剪切模量λ和拉梅常数μ；塑性参数包括线性屈服函数的斜率初值M0和斜率终值Mf、塑性势函数斜率
Figure BDA0003799319660000031
线性屈服函数与竖向Q轴的截距C、摩擦硬化控制参数κ、椭圆屈服函数的短轴A、初始椭圆屈服面的中心Pi0、终态压缩状态下塑性体量对数应变ε*和压缩硬化控制参数r；相场参数包括相场模型正则化参数l0、临界断裂能量释放率/>
Figure BDA0003799319660000032
相场粘性参数η、相场模型参数k和初始历史应变场函数/>
Figure BDA0003799319660000033
初始化物质点变量包括初始物质点位置向量xp，0、初始物质点位移向量up，0、初始物质点速度向量/>
Figure BDA0003799319660000034
初始物质点相场cp，0、物质点质量mp、初始物质点体积Vp，0、初始物质点变形梯度张量Fp，0和两个初始物质点粒子域向量r1，0和r2，0；
步骤二，初始化欧拉背景网格，通过CPDI插值技术建立离散物质点与欧拉背景网格之间的映射关系，然后将离散物质点上的物理材料属性映射到欧拉背景网格节点上；
CPDI插值技术中定义每个离散物质点都存在一个平行四边形粒子域，粒子域的变形更新为r1，n+1＝Fp，n+1·r1，0和r2，n+1＝Fp，n+1·r2，0，其中(r1，0，r2，0)和(r1，n+1，r2，n+1)分别表示初始时刻和当前时刻的物质点粒子域向量，Fp，n+1表示当前时刻的变形梯度张量；当前构型下的广义插值函数及其梯度解析表达为：
Figure BDA0003799319660000041
Figure BDA0003799319660000042
其中，(r1x，n+1，r1y，n+1)和(r2x，n+1，r2y，n+1)分别表示当前时刻tn+1的物质点粒子域向量r1，n+1和r2，n+1的分量，
Figure BDA0003799319660000043
为离散物质点p的粒子域角点i关于欧拉背景网格的插值基函数，/>
Figure BDA0003799319660000044
表示在当前构型下的梯度算子，Vp表示离散物质点p的粒子域体积，下标I表示背景网格节点；
步骤三，通过粒子接触算法判别接触物质点并计算接触力，将接触力映射到位移场的节点外力向量，再通过显式时间积分方法结合交错求解策略分别计算相场和位移场的离散控制方程，其中通过位移场的历史应变场函数
Figure BDA00037993196600000416
完成计算节点相场率/>
Figure BDA0003799319660000045
通过考虑物质点相场cp的耦合显式相场-塑性模型完成计算位移场的节点加速度/>
Figure BDA0003799319660000046
和节点速度/>
Figure BDA0003799319660000047
在显式时间积分方法中，将0到当前时刻tn+1内的时间间隔[0，tn+1]离散成若干时间增量0＝t0＜t1＜…＜tn＜tn+1，并定义当前时间增量Δtn+1＝tn+1-tn，时间增量满足克朗CFL条件以确保结果的稳定性；
利用交错求解策略分别求解相场和位移场离散控制方程的具体操作方式如下：
对于位移场，通过tn时刻的物质点位移向量up，n、物质点速度向量
Figure BDA0003799319660000048
和节点速度向量/>
Figure BDA0003799319660000049
更新tn+1时刻的物质点位移向量up，n+1、物质点速度向量/>
Figure BDA00037993196600000410
节点速度向量/>
Figure BDA00037993196600000411
和节点加速度向量/>
Figure BDA00037993196600000412
具体操作如下；
Figure BDA00037993196600000413
Figure BDA00037993196600000414
Figure BDA00037993196600000415
Figure BDA0003799319660000051
其中，MI为节点集中质量，
Figure BDA0003799319660000052
和/>
Figure BDA0003799319660000053
分别表示位移场的节点外力向量和节点内力向量，Nn表示总的背景网格节点数；
Figure BDA0003799319660000054
Figure BDA0003799319660000055
Figure BDA0003799319660000056
其中，ρp表示当前离散物质点密度，Jp为变形梯度张量F的Jacobian行列式，Np表示总的离散物质点数，τp表示Kirchhoff应力张量，G表示重力加速度向量，
Figure BDA0003799319660000057
为可替代的背景网格数值形函数，T表示指定外力向量；
对于相场，通过tn时刻的离散物质点相场cp，n和tn+1时刻的背景网格节点相场率
Figure BDA00037993196600000523
更新tn+1时刻的离散物质点相场cp，n+1，具体操作如下；
Figure BDA0003799319660000058
Figure BDA0003799319660000059
其中，DI为节点集中相场粘性，YI为节点相场残余力；
Figure BDA00037993196600000510
Figure BDA00037993196600000511
在粒子接触算法中，两个离散物质点
Figure BDA00037993196600000512
和/>
Figure BDA00037993196600000513
分别来自两个不同的固体，当两个固体进入接触状态时满足如下关系；
Figure BDA00037993196600000514
其中，
Figure BDA00037993196600000515
表示穿透距离，/>
Figure BDA00037993196600000516
和/>
Figure BDA00037993196600000517
分别表示离散物质点/>
Figure BDA00037993196600000518
和/>
Figure BDA00037993196600000519
的半径，定义为
Figure BDA00037993196600000520
Figure BDA00037993196600000521
Figure BDA00037993196600000522
表示二者相对距离；
在当前时间步n+1内根据局部线动量守恒法则存在
Figure BDA0003799319660000061
其中/>
Figure BDA0003799319660000062
和/>
Figure BDA0003799319660000063
分别表示接触力作用在离散物质点/>
Figure BDA0003799319660000064
和/>
Figure BDA0003799319660000065
的中心，因此，非侵彻条件(14)重构为
Figure BDA0003799319660000066
其中
Figure BDA0003799319660000067
Figure BDA0003799319660000068
其中，
Figure BDA0003799319660000069
和/>
Figure BDA00037993196600000610
分别表示离散物质点/>
Figure BDA00037993196600000611
和/>
Figure BDA00037993196600000612
的质量，Δtn+1表示当前增量时间步长；
将公式(16)和(17)带入公式(15)，获得接触力；
Figure BDA00037993196600000613
其中，
Figure BDA00037993196600000614
和/>
Figure BDA00037993196600000615
分别表示n+1时间步的法相接触力向量和单位法向向量，当考虑库伦摩擦接触条件时，有切向接触力向量
Figure BDA00037993196600000616
其中，
Figure BDA00037993196600000617
表示摩擦系数，/>
Figure BDA00037993196600000618
表示切向单位向量，
Figure BDA00037993196600000619
Figure BDA00037993196600000620
表示切线方向的离散物质点/>
Figure BDA00037993196600000621
和/>
Figure BDA00037993196600000622
的相对速度向量，/>
Figure BDA00037993196600000623
表示离散物质点/>
Figure BDA00037993196600000624
和/>
Figure BDA00037993196600000625
的相对速度向量；
因此，离散物质点
Figure BDA00037993196600000626
和/>
Figure BDA00037993196600000627
的总的接触力向量分别为
Figure BDA00037993196600000628
将离散物质点
Figure BDA00037993196600000629
和/>
Figure BDA00037993196600000630
的接触力向量分别映射到各自对应的节点外力向量/>
Figure BDA00037993196600000631
Figure BDA00037993196600000632
其中，下标II和
Figure BDA00037993196600000637
分别表示两个不同的固体的背景计算网格上的节点；
步骤四，在完成步骤三的当前时间步的相场和位移场离散控制方程计算后，将背景网格上相关物理信息映射回离散物质点用于更新离散物质点的相关变量cp、
Figure BDA00037993196600000633
up、Fp和
Figure BDA00037993196600000634
根据
Figure BDA00037993196600000635
对历史应变场函数进行更新，具体操作如下：
Figure BDA00037993196600000636
Figure BDA0003799319660000071
其中，
Figure BDA0003799319660000072
为裂纹驱动弹性应变能，裂纹驱动塑性应变能/>
Figure BDA0003799319660000073
其中α为1减Taylor-Quinney系数，Wp，tot表示总的塑性功，Jp为塑性变形梯度张量Fp的Jacobian行列式；
步骤五，存储和输出相关变量信息，返回步骤二，进入下一时间步，直至计算完成。
所述耦合显式相场-塑性模型根据显式相场断裂模型和光滑双屈服面塑性本构模型发展，其包含四个返回映射区域：锥顶返回映射区域、线性非关联返回映射区域、椭圆非关联返回映射区域和椭圆关联返回映射区域；包括两种硬化法则：线性屈服面的摩擦硬化法则和椭圆屈服面的压缩硬化法则；
其中，所述显式相场断裂模型依据微观力平衡法则和热力学第二定律得出：考虑一个含有裂纹面Γ的任意固体域，指定外力向量T＝P·N作用于外力边界
Figure BDA0003799319660000074
其中P表示第一PK应力张量，N表示边界/>
Figure BDA0003799319660000075
的外法线方向单位向量，/>
Figure BDA0003799319660000076
表示位移边界，在耦合损伤-塑性系统中设定存在用内部微观力π和微观表面力ζ表征的微观力平衡法则，则有限变形动态断裂问题在参考构型下的动量平衡方程和相场微观力平衡方程分别为
Figure BDA0003799319660000077
Divξ+π＝0，ξ·N＝0 (26)
其中，
Figure BDA0003799319660000078
表示加速度向量，ρ0表示初始密度，u表示位移向量，/>
Figure BDA0003799319660000079
表示边界/>
Figure BDA00037993196600000719
上的指定位移，Div(·)为参考构型下散度算子，ξ表示微观力系统的外部作用力；
基于上述平衡方程，从热力学分析角度出发，耦合损伤-塑性模型的能量平衡方程为
Figure BDA00037993196600000710
其中，e为内能的单位质量密度，
Figure BDA00037993196600000711
为相场率，/>
Figure BDA00037993196600000712
为变形梯度变化率张量；
然后，通过采用功共轭关系
Figure BDA00037993196600000713
在上式(27)中引入Kirchhoff应力张量τ和变形率张量d，基于热力学第二定律的耦合系统非负耗散不等式为/>
Figure BDA00037993196600000714
其中，
Figure BDA00037993196600000715
为能量耗散，Ψ为单位体积的存储能量密度，其不包含创建裂纹面的能量，则
Figure BDA00037993196600000716
其中，W表示无损伤材料的存储能量，q为类似应变的塑性内变量向量，
Figure BDA00037993196600000717
为塑性内变量向量变化率，/>
Figure BDA00037993196600000718
为弹性左Cauchy-Green应变张量be的变化率，退化函数g(c)＝(1-k)c2+k，k满足0＜k＜＜1；
因此，公式(28)重构为
Figure BDA0003799319660000081
其中，
Figure BDA0003799319660000082
为Lie导数算子；
设定公式(31)满足任何热力学过程，并且通过π＝πen+πdis将分解内部微观力π为能量部分πen和耗散部分πdis，得如下关系式
Figure BDA0003799319660000083
Figure BDA0003799319660000084
Figure BDA0003799319660000085
Figure BDA0003799319660000086
Figure BDA0003799319660000087
Figure BDA0003799319660000088
其中，
Figure BDA0003799319660000089
为塑性势函数，则流动法则定义为/>
Figure BDA00037993196600000810
Figure BDA00037993196600000811
为非负的塑性乘子的率；
将公式(34)和(35)带入公式(36)，重构断裂耗散不等式为
Figure BDA00037993196600000812
其中，正则化裂纹表面密度函数Γc的率形式为
Figure BDA00037993196600000813
其中，/>
Figure BDA00037993196600000814
通过构造如下局部阈值函数限制局部裂纹驱动力场
Figure BDA00037993196600000815
/>
其中，
Figure BDA00037993196600000816
β是与相场变量c相互作用的裂纹驱动力场，可知当tc＜0暗示着无耗散裂纹的累积；
然后，通过借鉴相场变分框架中相场变量的演化使能量耗散最大化假设，求解断裂耗散函数的约束优化问题
Figure BDA00037993196600000817
其中，Λ为拉格朗日乘子，
Figure BDA00037993196600000818
则产生如下约束条件
Figure BDA0003799319660000091
其次，将扩展的断裂耗散函数
Figure BDA0003799319660000092
代入能量平衡方程(27)取代微观力作用项
Figure BDA0003799319660000093
则非负的耦合系统耗散不等式(28)重写为
Figure BDA0003799319660000094
因此，公式(43)重构为
Figure BDA0003799319660000095
则得出如下关系式
Figure BDA0003799319660000096
Figure BDA0003799319660000097
Figure BDA0003799319660000098
Figure BDA0003799319660000099
Figure BDA00037993196600000910
Figure BDA00037993196600000911
通过引入相场粘性参数η＞0对断裂耗散公式(40)重构为
Figure BDA00037993196600000912
其中，相场粘性参数为确保相场断裂模拟的数值稳定性须满足
Figure BDA00037993196600000913
为Macaulay括号算子，下标+和-分别表示对断裂有贡献和无贡献部分，对上式采用<·>+是为了满足相场演化的不可逆条件，因此，此约束优化问题的必要条件是/>
Figure BDA00037993196600000914
将公式(46)带入上式(52)，并考虑引入历史应变场函数
Figure BDA00037993196600000915
保证弹塑性变形中裂纹扩展的不可逆性，得参考构型下所述显式相场断裂模型的强形式控制方程为
Figure BDA00037993196600000916
其中，历史应变场函数满足
Figure BDA00037993196600000917
W+为W中驱动裂纹扩展能量；
在相场弹塑性断裂分析中，拉-压分解后的裂纹驱动弹性应变能采用如下
Figure BDA0003799319660000101
则Kirchhoff应力张量推导为
Figure BDA0003799319660000102
其中，τA为Kirchhoff应力张量τ的分量；
借鉴耦合损伤-塑性力学中有效应力概念，引入相场有效Kirchhoff应力张量如下
Figure BDA0003799319660000103
其中
Figure BDA0003799319660000104
为相场有效Kirchhoff应力张量/>
Figure BDA0003799319660000105
的分量，则其体积不变量P和偏量不变量Q分别为
Figure BDA0003799319660000106
其中，
Figure BDA0003799319660000107
为偏量应力张量，I为二阶单位张量，trace(·)表示张量的迹算子；
定义Drucker-Prager线性屈服面函数及其摩擦硬化法则分别为
Figure BDA0003799319660000108
其中，C表示线性屈服函数与竖向Q轴的截距，M表示线性屈服函数的斜率，M0和Mf分别表示斜率M的初值和终值，κ表示摩擦硬化控制参数，γ表示累积塑性乘子；
定义Cap椭圆屈服面函数及其压缩硬化法则分别为
Figure BDA0003799319660000109
其中，Pi为椭圆屈服函数中心，
Figure BDA00037993196600001010
为塑性体量对数应变；则两个屈服面满足光滑连接的相切条件为/>
Figure BDA00037993196600001011
切点为/>
Figure BDA00037993196600001012
此外，针对椭圆关联返回映射区域采用塑性关联流动法则，则采用Cap椭圆屈服面函数
Figure BDA00037993196600001013
作为其塑性势函数，而对于线性非关联返回映射区域和椭圆非关联返回映射区域采用塑性非关联流动法则，则其塑性势函数定义如下
Figure BDA00037993196600001014
其中，存在满足
Figure BDA00037993196600001015
并且两个塑性势函数的切点为/>
Figure BDA00037993196600001016
所述耦合显式相场-塑性模型的具体实施流程如下：
步骤1、弹性变形预测：
Figure BDA00037993196600001017
其中，
Figure BDA0003799319660000111
表示相对变形梯度张量，Δu表示增量位移向量；
步骤2、应变谱分解：beTr＝n·λeTr·nT；
其中，λeTr和n分别表示预测弹性左Cauchy-Green应变张量beTr的特征值和特征向量；
步骤3、计算预测应变不变量
Figure BDA0003799319660000112
和/>
Figure BDA0003799319660000113
其中，预测弹性对数应变向量
Figure BDA0003799319660000114
其体积不变量和偏量不变量分别为/>
Figure BDA0003799319660000115
和/>
Figure BDA0003799319660000116
其中预测弹性对数应变向量的偏量应变向量
Figure BDA0003799319660000117
Figure BDA0003799319660000118
和单位向量δ＝(1，1，1)；
步骤4、根据公式(56)和(57)预测弹性应力不变量PTr和QTr；
步骤5、根据公式(59)检查：
Figure BDA0003799319660000119
步骤5.1当
Figure BDA00037993196600001110
根据公式(58)检查/>
Figure BDA00037993196600001111
步骤5.1.1当
Figure BDA00037993196600001112
进入线性非关联返回映射区域；
步骤5.1.1.1采用线性非关联返回映射算法：计算
Figure BDA00037993196600001113
和Δγ，并根据公式(58)更新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 BDA00037993196600001114
检查应力状态是否处于弹性区域PTr≥P*；
步骤5.1.2.1当PTr≥P*：进入弹性变形步骤7；
步骤5.1.2.2当PTr＜P*：进入步骤6；
步骤5.2当
Figure BDA00037993196600001115
进入弹性变形步骤7；
步骤6对每一个迭代步检查应力状态是否在椭圆非关联返回映射区域或椭圆关联返回映射区域：
Figure BDA00037993196600001116
步骤6.1当
Figure BDA00037993196600001117
进入椭圆非关联返回映射区域；
步骤6.1.1采用椭圆非关联返回映射算法：计算
Figure BDA00037993196600001118
和Δγ，根据公式(58)更新M，并检查/>
Figure BDA00037993196600001119
/>
步骤6.1.1.1当
Figure BDA00037993196600001120
根据公式(59)更新Pi，并进入步骤8；
步骤6.1.1.2当
Figure BDA00037993196600001121
Cap椭圆压缩塑性模型失效，进入纯剪切Drucker-Prager塑性模型；
步骤6.2当
Figure BDA0003799319660000121
进入椭圆关联返回映射区域；
步骤6.2.1采用椭圆关联返回映射算法：计算
Figure BDA0003799319660000122
和Δγ，根据公式(58)更新M，并检查/>
Figure BDA0003799319660000123
步骤6.2.1.1当
Figure BDA0003799319660000124
根据公式(59)更新Pi，并进入步骤8；
步骤6.2.1.2当
Figure BDA0003799319660000125
Cap椭圆压缩塑性模型失效，进入纯剪切Drucker-Prager塑性模型；
步骤7弹性变形：设定M＝Mn、Pi＝Pi，n、
Figure BDA0003799319660000126
和(P，Q)＝(PTr，QTr)
步骤8计算弹性主对数应变向量εe和Kirchhoff应力张量的分量τA；
步骤9计算弹性左Cauchy-Green应变张量be＝n·exp(2εe)·nT，并根据公式(55)更新Kirchhoff应力张量τ。
基于有限应变假设，如图1所示，参考构型Ω中点X将通过变形映射函数
Figure BDA0003799319660000127
映射到t时刻的当前构型Ωt中点x，则变形梯度张量为
Figure BDA0003799319660000128
对变形梯度张量F实施乘法分解，则有
F＝FeFp (62)
其中，上标e和p分别表示弹性变形和塑性变形相关变量，后续公式中上角标相同标注代表相同含义；
然后，弹性主对数应变向量εe可以表示为
Figure BDA0003799319660000129
其中，λe和n分别表示弹性左Cauchy-Green应变张量be的特征值和特征向量，此外，弹性对数应变向量的体积不变量和偏量不变量分别为
Figure BDA00037993196600001210
和/>
Figure BDA00037993196600001211
其中弹性对数应变向量的偏量应变向量/>
Figure BDA00037993196600001212
和单位向量δ＝(1，1，1)；
在相场断裂模型中，裂纹可以通过如下正则化裂纹表面密度函数Γc表示
Figure BDA00037993196600001213
其中，相场变量c∈[0，1]，c＝0表示材料完全断裂，c＝1表示材料完整无损伤，
Figure BDA00037993196600001214
表示在参考构型下的梯度算子，l0表示控制裂纹宽度的相场模型正则化参数，当l0→0则Γc由弥散裂纹逼近离散裂纹；
则裂纹表面演化速率可以表示为
Figure BDA0003799319660000131
/>
由此可得
Figure BDA0003799319660000132
为相场演化不可逆限制条件，是相场变量c及其梯度/>
Figure BDA0003799319660000133
演化的限制性条件，其具体可以表达为
Figure BDA0003799319660000134
其中，
Figure BDA0003799319660000135
表示相场率；
所述参考构型下所述显式相场断裂模型的强形式控制方程采用Galerkin法在物质点上进行空间离散，同时考虑两个转变关系引入Kirchhoff应力张量τ＝P·FT和当前密度ρ＝J-1ρ0，J为变形梯度张量F的Jacobian行列式，则动量平衡方程在当前构型下的弱形式为
Figure BDA0003799319660000136
其中，ω为位移场任意权函数，然后上式(71)在物质点可以离散为
Figure BDA0003799319660000137
其中，
Figure BDA0003799319660000138
表示在当前构型下的梯度算子，Vp表示物质点p的粒子域体积，φI、/>
Figure BDA0003799319660000139
分别表示背景网格广义插值数值基函数及其梯度，/>
Figure BDA00037993196600001310
为可替代的背景网格数值形函数，Nn和Np分别表示总的网格节点数和总的物质点数，下标p、I分别表示与物质点和背景网格节点相关的变量；
通过对公式(72)中左端第一项采用集中质量矩阵，且考虑ω的任意性，则上式可以构造为如下矩阵形式
Figure BDA00037993196600001311
其中，
Figure BDA00037993196600001312
表示节点加速度向量，MI为节点集中质量，/>
Figure BDA00037993196600001313
和/>
Figure BDA00037993196600001314
分别表示位移场的节点外力向量和节点内力向量
同理，对显式相场控制方程构造当前构型下的弱形式方程为
Figure BDA00037993196600001315
其中，
Figure BDA00037993196600001416
为相场任意权函数，则对上式(741)在物质点上进行离散可得/>
Figure BDA0003799319660000141
相应的，对上式(75)中左端项采用集中相场粘性矩阵，且考虑
Figure BDA00037993196600001417
的任意性，则上式可以构造为如下矩阵形式
Figure BDA0003799319660000142
其中，
Figure BDA0003799319660000143
为节点相场率，DI为节点集中相场粘性，YI为节点相场残余力。
根据上述理论推导，给出了本发明提出的大变形框架下显式相场物质点法的详细求解格式，其通过结合耦合显式相场-塑性模型、对流粒子域插值技术(CPDI)、粒子接触算法和交错求解策略即可实现本发明所要发展的动态冲击/接触弹塑性大变形断裂分析显式相场物质点法。
根据图6所示本发明提出的ePF-CPDI方法的计算流程图，其具体的实施过程如下所示：
1)建立物质点离散模型和定义材料参数(弹性：λ，λ；塑性：M0，Mf，
Figure BDA0003799319660000144
C，κ，A，Pi0，ε*，r；相场：l0，/>
Figure BDA00037993196600001418
η，k，/>
Figure BDA0003799319660000145
初始变量：xp，0，up，0，/>
Figure BDA0003799319660000146
Cp，0，mp，Vp，0，Fp，0，r1，0，r2，0)；
2)时间步循环n＝1，2，…，Nload
2.1)初始化背景网格：根据式(1)和(2)计算插值函数φIp及其梯度
Figure BDA0003799319660000147
2.2)根据公式(14)搜索接触物质点，并根据公式(18)-(20)计算物质点接触力
Figure BDA0003799319660000148
2.3)映射物质点相场cp，n到网格节点cI，n，并计算其梯度
Figure BDA0003799319660000149
2.4)根据公式(12)和(13)计算节点集中相场粘性矩阵DI，n和节点相场残余力向量YI，n；
2.5)更新节点相场率
Figure BDA00037993196600001410
2.6)根据公式(11)更新物质点相场cp，n；
2.7)映射物质点动量
Figure BDA00037993196600001411
到网格节点/>
Figure BDA00037993196600001412
2.8)计算节点外力向量
Figure BDA00037993196600001413
并根据公式(21)和(22)映射物质点接触力/>
Figure BDA00037993196600001414
2.9)根据公式(7)计算节点集中质量矩阵MI，n，根据耦合显式相场-塑性模型和公式(8)计算节点内力向量
Figure BDA00037993196600001415
2.10)施加Dirichlet边界条件；
2.11)计算节点加速度向量
Figure BDA0003799319660000151
2.12)根据公式(4)-(6)更新物质点变量：
Figure BDA0003799319660000152
up，n和xp，n；
2.13)更新物质点变形梯度张量Fp，n和体积Vp，n；
2.14)更新粒子域向量r1，n和r2，n；
2.15)根据公式(23)计算W+，并更新
Figure BDA0003799319660000153
3)返回步骤2)直至计算结束；
其中，Nload表示总的时间步数。
本发明的有益效果：(1)本发明提供的一种动态冲击/接触弹塑性大变形断裂分析显式相场物质点法，为岩土材料动态断裂破坏研究提供了一种全新的数值计算方法，由于该方法搭载显式物质点法，相较传统基于网格类方法，能有效克服网格畸变问题，故其显著优势在于能够很好地处理大变形和接触等强非线性断裂破坏问题，并且能自动处理复杂的裂纹分叉、交汇及三维空间自由扩展等问题，且该方法还能通过替换塑性本构模型扩展至其他材料断裂破坏分析，并能通过嵌入多物理场耦合理论扩展至复杂多场耦合断裂破坏分析，如液固耦合断裂分析、热力耦合断裂分析和热-水-力三场耦合断裂破坏分析等；
(2)本发明提供的一种动态冲击/接触弹塑性大变形断裂分析显式相场物质点法，在该方法中，发展了一种显式率相关相场断裂模型，它的理论推导过程由于是基于微观力平衡法则和热力学第二定律，故其不涉及最大塑性耗散法则，适用针对于非关联流动特性的岩土材料的断裂破坏分析，并且均适用于预测脆性和弹塑性断裂破坏行为，此外，该显式相场断裂模型还能推广至其他材料的断裂破坏研究，如金属、凝胶软材料等；
(3)本发明提供的一种动态冲击/接触弹塑性大变形断裂分析显式相场物质点法，采用光滑双屈服面塑性本构模型和显式相场断裂模型发展了一种耦合显式相场-塑性模型，可以有效预测岩土材料的复杂脆性-塑性转变失效行为，相较传统耦合损伤塑性本构模型，其数值实施较为容易，并且该耦合显式相场-塑性模型中的光滑双屈服面塑性本构模型通过简单操作即可退化至单屈服面的线性剪切Drucker-prager塑性本构模型，可兼容单屈服面塑性本构模型和多屈服面塑性本构模型的多种用途，实现一模两用功能；
(4)本发明提供的一种动态冲击/接触弹塑性大变形断裂分析显式相场物质点法，采用相场-位移场交错求解策略结合显式时间积分方案，其相较隐式交错迭代求解方案，极大程度上降低了数值实施的难度，且提高了计算效率，为大规模复杂断裂破坏问题的并行计算研究提供了可行性方案。
附图说明
图1为本发明的含有裂纹界面Γ的任意固体域Ω的变形映射示意图；
图2为本发明提出的ePF-CPDI方法计算循环示意图；
图3为本发明的CPDI插值技术中物质点与背景网格联系示意图；
■为背景网格节点；●为初始物质点；
Figure BDA0003799319660000162
为MP与GN建立联系；/>
Figure BDA0003799319660000163
为初始粒子域；
Figure BDA0003799319660000164
为粒子域向量；●为变形后物质点；·为粒子域角点；/>
Figure BDA0003799319660000167
为变形后粒子域；
图4为本发明的粒子接触算法示意图；
图5为本发明的光滑双屈服面塑性本构模型的应力返回映射示意图；
图6为本发明的一种动态冲击/接触弹塑性大变形断裂分析显式相场物质点法(ePF-CPDI)的操作流程图；
图7为本发明的动态脆性断裂分析实施例1结构和边界条件示意图；
图8为本发明的实施例1的由三种数值方法得到的(a)弹性应变能、(b)断裂耗散能和(c)裂纹尖端速度时程曲线对比结果图；
图9为本发明的实施例1的不同时刻相场断裂演化云图：(a)t＝2.5×10-5s、(b)t＝5.0×10-5s、(c)t＝6.5×10-5s和(d)t＝8.0×10-5s，其中相场值小于c＜0.05的物质点已被去除，且结构变形被方大50倍；
图10为本发明的准静态弹塑性断裂分析实施例2结构和边界条件示意图；
图11为本发明的实施例2的由三种数值方法得到的(a)名义偏差应力和(b)名义体积应变随名义竖向应变关系曲线对比图；
图12为本发明的实施例2的两种网格宽度下的终态断裂结果与PF-ICPDI方法对比结果：(a-c)等效塑性应变云图和(d-f)相场断裂云图；
图13为本发明的多体动态接触的实施例3结构和边界条件示意图；
图14为本发明的实施例3的三种工况下断裂耗散能时程曲线；
图15为本发明的实施例3的工况1下不同时刻(a-f)等效塑性应变云图和(g-l)相场断裂云图：(a,g)t＝0s、(b,h)t＝3.5×10-5s、(c,i)t＝6.5×10-5s、(d,j)t＝9.5×10-5s、(e,k)t＝1.25×10-4s和(f,l)t＝1.5×10-4s，其中相场值小于c＜0.1的物质点已被去除；
图16为本发明的实施例3的工况2下不同时刻(a-f)等效塑性应变云图和(g-l)相场断裂云图：(a,g)t＝0s、(b,h)t＝3.5×10-5s、(c,i)t＝6.5×10-5s、(d,j)t＝9.5×10-5s、(e,k)t＝1.25×10-4s和(f,l)t＝1.5×10-4s，其中相场值小于c＜0.1的物质点已被去除。
具体实施方式
下面结合附图和实施例对本发明的性能做出进一步详细说明。以下实施例用于说明本发明，但不能用来限制本发明的适用范围。
为了使本发明的目的、技术方案和具体实施效果展示更加清晰明了，下面通过三个具体实施例结合附图7～16对本发明提出的ePF-CPDI方法的准确性、可靠性和优异性能作进一步的详细说明。
首先，我们将实施一个动态脆性断裂模拟来验证本发明所提出的显式相场断裂模型的准确性和有效性。然后，考虑了一个准静态弹塑性断裂实施例来说明所发展的耦合显式相场-塑性模型的准确性。此外，还通过开展多体动态接触弹塑性断裂实施例来进一步说明本发明提出的ePF-CPDI方法在处理岩土材料复杂断裂失效问题的可靠性和优异性能。前两个实施例中的参考解分别来自于我们课题组提出的PF-ICPDI方法、Kakouris等人和Choo等人的工作，并且在第一个实施例中关于统计裂纹尖端扩展速度采用vtip，n＝(xtip，n-xtip，n-1)/Δtn，其中xtip为相场值c＝0.25等值线所处裂纹尖端的物质点坐标。在所有实施例中，均采用双线性四边形网格单元并在其内部设置至少2×2个物质点，且设置相场模型参数l0大于等于背景网格尺寸h来满足相场断裂数值模拟结果合理性。所有实施例均不考虑重力效应。
(1)实施例1：预制裂纹矩形板动态拉伸脆性断裂实验(附图7～9)
实施例1考虑一个标准的预制裂纹矩形板动态拉伸脆性断裂算例，其几何结构和边界条件如图7所示。通过在中心裂纹线两侧各删除两排网格内的物质点以实现预制裂纹，并上下两端施加等大反向的拉伸冲击载荷σt＝1MPa。对预制裂纹矩形板考虑三种离散划分：254400、1020800和4089600个物质点划分，相对应网格尺寸h＝0.25、0.125和0.0625mm。增量时间步长设定Δt＝1×10-8s，瑞利波速为vR＝2125m/s。此算例考虑平面应变假设条件。材料参数如下表1所示。
表1.预制裂纹矩形板动态拉伸脆性断裂实验材料参数表
Figure BDA0003799319660000171
E-杨氏模量、v-泊松比
图8给出了分别采用ePF-CPDI方法、PF-ICPDI方法和Kakouris等人工作得到的弹性应变能、断裂耗散能和裂纹尖端速度时程曲线对比结果，我们可以看出三者结果吻合良好。由于ePF-CPDI方法采用与Kakouris等人工作一致的网格宽度h＝0.125mm，而PF-ICPDI方法中采用较稀疏的网格尺寸h＝0.25mm，故ePF-CPDI方法得到的结果相比PF-ICPDI方法更接近于Kakouris等人工作中的结果。此外，我们还通过计算效率对比发现，PF-ICPDI方法采用较稀疏的网格尺寸计算一个载荷步所花费的时间是采用精细网格尺寸的ePF-CPDI方法的25倍，即PF-ICPDI方法用时60s，ePF-CPDI方法用时2.4s，究其原因是ePF-CPDI方法采用显式时间积分求解策略，无需N-R迭代求解矩阵线性方程组的计算优势，这也反映出ePF-CPDI方法具有极大潜力去求解大规模复杂断裂问题。图9展示了不同时刻变形体内相场断裂演化云图，其结果与Kakouris等人工作一致。综合上述分析，证明了本发明提出的显式相场断裂模型的准确性和有效性。
(2)实施例2：平面应变压缩弹塑性断裂失效实验(附图10～12)
实施例2模拟如图10所示的岩石结构在围压作用下的平面应变压缩失效问题，具体结构和边界条件如图所示，图中圆圈为引导结构非均匀变形的弱化区域，且设置圆圈内所有物质点的内聚力强度C为其他区域的98％。试件离散为120×300个物质点，相应的背景网格离散为边长为h＝0.5mm的70×154个四结点矩形单元，且每个网格内设置2×2个物质点。在其左、右边界施加围压σc＝5Mpa，上边界施加竖直向下的位移边界条件Δu，在该方法中，竖向位移将通过常速度形式施加v0＝Δu/Δt，其中v0＝0.2m/s、Δt＝2.5×10-8s。具体材料如下表2所示。
表2.平面应变压缩弹塑性断裂失效实验材料参数表
Figure BDA0003799319660000181
图11展示了本发明的ePF-CPDI方法得出的名义偏差应力、名义体积应变随名义轴向应变变化曲线与参考解的对比图，从图中可以清晰的看出三种方法的结果均吻合良好。此外，图11和12显示了采用两种网格尺寸得到的数值结果几乎完全一致，从而说明了本发明提出的ePF-CPDI方法在分析弹塑性断裂问题时无网格敏感性特性。因此，上述分析结果验证了本发明提出的耦合显式相场-塑性模型和ePF-CPDI方法的准确性。
(3)实施例3：多体接触弹塑性压缩断裂破坏实验(附图13～16)
实施例3考虑了如图13所示多体接触弹塑性压缩断裂破坏模拟，其具体结构和边界条件如图所示，采用平面应变假设条件。背景网格尺寸设定h＝0.05mm并采用网格物质点密度2×2的均匀划分方式。整体多体结构共离散为203108个物质点，并在上端平板结构施加一个竖直向下的速度边界v＝16m/s，两侧和下边界均认为是刚体和完全固定边界。在此算例中，对两个固体颗粒考虑两种刚度工况：(i)两个颗粒杨氏模量相同E1＝E2＝10500MPa、(ii)颗粒2比颗粒1的杨氏模量小E1＝10500MPa＞E2＝3500MPa。增量时间步长设定Δt＝5×10-9s。两个颗粒的材料参数如下表3所示。
表3.多体接触弹塑性压缩断裂破坏实验材料参数表
Figure BDA0003799319660000182
Figure BDA0003799319660000191
图14给出了两种工况的耗散能量时程曲线，图15和16分别展示了两种工况的不同时刻颗粒的变形失效行为。从给出的两种工况数值结果可以看出，当上端平板下移接触到上面颗粒1后致使颗粒1先发生破损，而在终态失效状态下下面颗粒2出现变形最大且破损最严重的现象。在工况1中(如图15所示)，上面颗粒1先出现贯穿裂纹如图15(c，i)所示，当上端继续压缩，颗粒2紧跟着发生破损。在工况2中(如图16所示)，由于下面颗粒2较软，故其先出现贯穿裂纹并遭受更多的压缩变形。通过两种工况对比，我们还能看出当减小颗粒2的刚度，颗粒1的破损情况出现明显好转。上述结果分析与工程实际应用中缓冲吸能结构的工作原理是一致的，从而说明了本发明提出的ePF-CPDI方法在分析岩土材料接触弹塑性断裂失效问题的有效性和优异性能。
综上所述，上述三个实施例分别从不同层次详细验证了本发明所提出的ePF-CPDI方法的准确性和有效性，并从计算效率和计算规模方面展示了该方法的显著优势，也证明了其发展的必要性。与此同时，还很好地说明了本发明所发展的显式相场断裂模型既能退化用于分析脆性断裂问题，也能用于岩土材料的弹塑性断裂失效模拟。并通过第三个拓展性实施例有效的说明了本发明在处理岩土材料复杂弹塑性接触断裂失效问题的能力。因此，本发明提出的动态冲击/接触弹塑性大变形断裂分析显式相场物质点法(ePF-CPDI)是一种极具发展前景的高性能数值方法。
本发明的实施例是为了示例和描述起见而给出的，而并不是无遗漏的或者将本发明限于所公开的形式。很多修改和变化对于本领域的普通技术人员而言是显而易见的。选择和描述实施例是为了更好的说明本发明的原理和实际应用，并且使本领域的普通技术人员能够理解本发明从而设计适于特定用途的带有各种修改的各种实施例。

Claims (3)

1.一种动态冲击/接触弹塑性大变形断裂分析显式相场物质点法，其特征在于，具体步骤如下：应用于岩土材料动态冲击/接触弹塑性断裂破坏高效数值分析；
步骤一，定义欧拉背景网格，建立离散物质点模型并定义离散物质点的物理材料参数，初始化物质点变量；物理材料参数包括弹性参数、塑性参数和相场参数；弹性参数包括剪切模量λ和拉梅常数μ；塑性参数包括线性屈服函数的斜率初值Mo和斜率终值Mf、塑性势函数斜率
Figure FDA00041886990100000112
线性屈服函数与竖向Q轴的截距C、摩擦硬化控制参数κ、椭圆屈服函数的短轴A、初始椭圆屈服面的中心Pi0、终态压缩状态下塑性体量对数应变ε*和压缩硬化控制参数r；相场参数包括相场模型正则化参数l0、临界断裂能量释放率/>
Figure FDA0004188699010000011
相场粘性参数η、相场模型参数k和初始历史应变场函数/>
Figure FDA0004188699010000012
初始化物质点变量包括初始物质点位置向量xp，0、初始物质点位移向量up，0、初始物质点速度向量/>
Figure FDA0004188699010000013
初始物质点相场cp，0、物质点质量mp、初始物质点体积Vp，0、初始物质点变形梯度张量Fp，0和两个初始物质点粒子域向量r1，0和r2，0；
步骤二，初始化欧拉背景网格，通过CPDI插值技术建立离散物质点与欧拉背景网格之间的映射关系，然后将离散物质点上的物理材料属性映射到欧拉背景网格节点上；
CPDI插值技术中定义每个离散物质点都存在一个平行四边形粒子域，粒子域的变形更新为r1，n+1＝Fp，n+1·r1，0和r2，n+1＝Fp，n+1·r2，0，其中(r1，0，r2，0)和(r1，n+1，r2，n+1)分别表示初始时刻和当前时刻的物质点粒子域向量，Fp，n+1表示当前时刻的变形梯度张量；当前构型下的广义插值函数及其梯度解析表达为：
Figure FDA0004188699010000014
Figure FDA0004188699010000015
其中，(r1x，n+1，r1y，n+1)和(r2x，n+1，r2y，n+1)分别表示当前时刻tn+1的物质点粒子域向量r1，n+1和r2，n+1的分量，
Figure FDA0004188699010000016
为离散物质点p的粒子域角点i关于欧拉背景网格的插值基函数，/>
Figure FDA0004188699010000017
表示在当前构型下的梯度算子，Vp表示离散物质点p的粒子域体积，下标I表示背景网格节点；
步骤三，通过粒子接触算法判别接触物质点并计算接触力，将接触力映射到位移场的节点外力向量，再通过显式时间积分方法结合交错求解策略分别计算相场和位移场的离散控制方程，其中通过位移场的历史应变场函数
Figure FDA0004188699010000018
完成计算节点相场率/>
Figure FDA0004188699010000019
通过考虑物质点相场cp的耦合显式相场-塑性模型完成计算位移场的节点加速度/>
Figure FDA00041886990100000110
和节点速度/>
Figure FDA00041886990100000111
在显式时间积分方法中，将0到当前时刻tn+1内的时间间隔[0，tn+1]离散成若干时间增量0＝t0＜t1＜…＜tn＜tn+1，并定义当前时间增量Δtn+1＝tn+1-tn，时间增量满足克朗CFL条件以确保结果的稳定性；
利用交错求解策略分别求解相场和位移场离散控制方程的具体操作方式如下：
对于位移场，通过tn时刻的物质点位移向量up，n、物质点速度向量
Figure FDA0004188699010000021
和节点速度向量
Figure FDA0004188699010000022
更新tn+1时刻的物质点位移向量up，n+1、物质点速度向量/>
Figure FDA0004188699010000023
节点速度向量/>
Figure FDA0004188699010000024
和节点加速度向量/>
Figure FDA0004188699010000025
具体操作如下；/>
Figure FDA0004188699010000026
Figure FDA0004188699010000027
Figure FDA0004188699010000028
Figure FDA0004188699010000029
其中，MI为节点集中质量，
Figure FDA00041886990100000210
和/>
Figure FDA00041886990100000211
分别表示位移场的节点外力向量和节点内力向量，Nn表示总的背景网格节点数；
Figure FDA00041886990100000212
Figure FDA00041886990100000213
Figure FDA00041886990100000214
其中，ρp表示当前离散物质点密度，Jp为变形梯度张量F的Jacobian行列式，Np表示总的离散物质点数，τp表示Kirchhoff应力张量，G表示重力加速度向量，
Figure FDA00041886990100000215
为可替代的背景网格数值形函数，T表示指定外力向量；
对于相场，通过tn时刻的离散物质点相场cp，n和tn+1时刻的背景网格节点相场率
Figure FDA00041886990100000216
更新tn+1时刻的离散物质点相场cp，n+1，具体操作如下；
Figure FDA00041886990100000217
Figure FDA00041886990100000218
其中，DI为节点集中相场粘性，YI为节点相场残余力；
Figure FDA0004188699010000031
Figure FDA0004188699010000032
在粒子接触算法中，两个离散物质点
Figure FDA0004188699010000033
和/>
Figure FDA0004188699010000034
分别来自两个不同的固体，当两个固体进入接触状态时满足如下关系；/>
Figure FDA0004188699010000035
其中，
Figure FDA0004188699010000036
表示穿透距离，/>
Figure FDA0004188699010000037
和/>
Figure FDA0004188699010000038
分别表示离散物质点/>
Figure FDA0004188699010000039
和/>
Figure FDA00041886990100000310
的半径，定义为/>
Figure FDA00041886990100000311
Figure FDA00041886990100000312
Figure FDA00041886990100000313
表示二者相对距离；
在当前时间步n+1内根据局部线动量守恒法则存在
Figure FDA00041886990100000314
其中/>
Figure FDA00041886990100000315
和
Figure FDA00041886990100000316
分别表示接触力作用在离散物质点/>
Figure FDA00041886990100000317
和/>
Figure FDA00041886990100000318
的中心，因此，非侵彻条件(14)重构为
Figure FDA00041886990100000319
其中
Figure FDA00041886990100000320
Figure FDA00041886990100000321
其中，
Figure FDA00041886990100000322
和/>
Figure FDA00041886990100000323
分别表示离散物质点/>
Figure FDA00041886990100000324
和/>
Figure FDA00041886990100000325
的质量，Δtn+1表示当前增量时间步长；
将公式(16)和(17)带入公式(15)，获得接触力
Figure FDA00041886990100000326
其中，
Figure FDA00041886990100000327
和/>
Figure FDA00041886990100000328
分别表示n+1时间步的法向接触力向量和单位法向向量，当考虑库伦摩擦接触条件时，有切向接触力向量
Figure FDA00041886990100000329
其中，θ表示摩擦系数，
Figure FDA00041886990100000330
表示切向单位向量，
Figure FDA00041886990100000331
Figure FDA00041886990100000332
表示切线方向的离散物质点/>
Figure FDA00041886990100000333
和/>
Figure FDA00041886990100000334
的相对速度向量，/>
Figure FDA00041886990100000335
表示离散物质点/>
Figure FDA00041886990100000336
和/>
Figure FDA00041886990100000337
的相对速度向量；
因此，离散物质点
Figure FDA00041886990100000338
和/>
Figure FDA00041886990100000339
的总的接触力向量分别为
Figure FDA0004188699010000041
将离散物质点
Figure FDA0004188699010000042
和/>
Figure FDA0004188699010000043
的接触力向量分别映射到各自对应的节点外力向量
Figure FDA0004188699010000044
Figure FDA0004188699010000045
其中，下标
Figure FDA0004188699010000046
和/>
Figure FDA0004188699010000047
分别表示两个不同的固体的背景计算网格上的节点；
步骤四，在完成步骤三的当前时间步的相场和位移场离散控制方程计算后，将背景网格上相关物理信息映射回离散物质点用于更新离散物质点的相关变量cp、
Figure FDA0004188699010000048
up、Fp和/>
Figure FDA0004188699010000049
根据
Figure FDA00041886990100000410
对历史应变场函数进行更新，具体操作如下：
Figure FDA00041886990100000411
/>
Figure FDA00041886990100000412
其中，
Figure FDA00041886990100000413
为裂纹驱动弹性应变能，裂纹驱动塑性应变能/>
Figure FDA00041886990100000414
其中α为1减Taylor-Quinney系数，Wp，tot表示总的塑性功，Jp为塑性变形梯度张量Fp的Jacobian行列式；
步骤五，存储和输出相关变量信息，返回步骤二，进入下一时间步，直至计算完成。
2.根据权利要求1所述的动态冲击/接触弹塑性大变形断裂分析显式相场物质点法，其特征在于，
所述耦合显式相场-塑性模型根据显式相场断裂模型和光滑双屈服面塑性本构模型发展，其包含四个返回映射区域：锥顶返回映射区域1、线性非关联返回映射区域2、椭圆非关联返回映射区域3和椭圆关联返回映射区域4；包括两种硬化法则：线性屈服面的摩擦硬化法则和椭圆屈服面的压缩硬化法则；
其中，所述显式相场断裂模型依据微观力平衡法则和热力学第二定律得出：考虑一个含有裂纹面Γ的任意固体域，指定外力向量T＝P·N作用于外力边界
Figure FDA00041886990100000415
其中P表示第一PK应力张量，N表示边界/>
Figure FDA00041886990100000416
的外法线方向单位向量，/>
Figure FDA00041886990100000417
表示位移边界，在耦合损伤-塑性系统中设定存在用内部微观力π和微观表面力ζ表征的微观力平衡法则，则有限变形动态断裂问题在参考构型下的动量平衡方程和相场微观力平衡方程分别为
Figure FDA00041886990100000418
Divξ+π＝0，ξ·N＝0 (26)
其中，
Figure FDA0004188699010000051
表示加速度向量，ρ0表示初始密度，u表示位移向量，/>
Figure FDA0004188699010000052
表示边界/>
Figure FDA0004188699010000053
上的指定位移，Div(·)为参考构型下散度算子，ξ表示微观力系统的外部作用力；
基于上述平衡方程，从热力学分析角度出发，耦合损伤-塑性模型的能量平衡方程为
Figure FDA0004188699010000054
其中，e为内能的单位质量密度，
Figure FDA0004188699010000055
为相场率，/>
Figure FDA0004188699010000056
为变形梯度变化率张量；
然后，通过采用功共轭关系
Figure FDA0004188699010000057
在上式(27)中引入Kirchhoff应力张量τ和变形率张量d，基于热力学第二定律的耦合系统非负耗散不等式为
Figure FDA0004188699010000058
其中，
Figure FDA0004188699010000059
为能量耗散，Ψ为单位体积的存储能量密度，其不包含创建裂纹面的能量，则
Ψ(be，q，c)＝g(c)W(be，q) (29)
Figure FDA00041886990100000510
其中，W表示无损伤材料的存储能量，q为类似应变的塑性内变量向量，
Figure FDA00041886990100000511
为塑性内变量向量变化率，/>
Figure FDA00041886990100000512
为弹性左Cauchy-Green应变张量be的变化率，退化函数g(c)＝(1-k)c2+k，k满足0＜k＜＜1；
因此，公式(28)重构为
Figure FDA00041886990100000513
其中，
Figure FDA00041886990100000514
为Lie导数算子；
设定公式(31)满足任何热力学过程，并且通过π＝πen+πdis将分解内部微观力π为能量部分πen和耗散部分πdis，得如下关系式
Figure FDA00041886990100000515
Figure FDA00041886990100000516
Figure FDA00041886990100000517
Figure FDA00041886990100000518
Figure FDA00041886990100000519
Figure FDA00041886990100000520
其中，
Figure FDA00041886990100000521
为塑性势函数，则流动法则定义为/>
Figure FDA00041886990100000522
为非负的塑性乘子的率；
将公式(34)和(35)带入公式(36)，重构断裂耗散不等式为
Figure FDA0004188699010000061
其中，正则化裂纹表面密度函数Γc的率形式为
Figure FDA0004188699010000062
其中，/>
Figure FDA0004188699010000063
通过构造如下局部阈值函数限制局部裂纹驱动力场
Figure FDA0004188699010000064
其中，
Figure FDA0004188699010000065
β是与相场变量c相互作用的裂纹驱动力场；
然后，通过借鉴相场变分框架中相场变量的演化使能量耗散最大化假设，求解断裂耗散函数的约束优化问题
Figure FDA0004188699010000066
其中，Λ为拉格朗日乘子，
Figure FDA0004188699010000067
则产生如下约束条件
Figure FDA0004188699010000068
其次，将扩展的断裂耗散函数
Figure FDA0004188699010000069
代入能量平衡方程(27)取代微观力作用项
Figure FDA00041886990100000610
/>
则非负的耦合系统耗散不等式(28)重写为
Figure FDA00041886990100000611
因此，公式(43)重构为
Figure FDA00041886990100000612
则得出如下关系式
Figure FDA00041886990100000613
Figure FDA00041886990100000614
Figure FDA00041886990100000615
Figure FDA00041886990100000616
Figure FDA00041886990100000617
Figure FDA0004188699010000071
通过引入相场粘性参数η＞0对断裂耗散公式(40)重构为
Figure FDA0004188699010000072
其中，相场粘性参数为确保相场断裂模拟的数值稳定性须满足
Figure FDA0004188699010000073
为Macaulay括号算子，下标+和-分别表示对断裂有贡献和无贡献部分，对上式采用<·>+是为了满足相场演化的不可逆条件，因此，此约束优化问题的必要条件是
Figure FDA0004188699010000074
将公式(46)带入上式(52)，并考虑引入历史应变场函数
Figure FDA0004188699010000075
保证弹塑性变形中裂纹扩展的不可逆性，得参考构型下所述显式相场断裂模型的强形式控制方程为
Figure FDA0004188699010000076
其中，历史应变场函数满足
Figure FDA0004188699010000077
W+为W中驱动裂纹扩展能量；
在相场弹塑性断裂分析中，拉-压分解后的裂纹驱动弹性应变能采用如下
Figure FDA0004188699010000078
/>
则Kirchhoff应力张量推导为
Figure FDA0004188699010000079
其中，τA为Kirchhoff应力张量τ的分量；
借鉴耦合损伤-塑性力学中有效应力概念，引入相场有效Kirchhoff应力张量如下
Figure FDA00041886990100000710
其中
Figure FDA00041886990100000711
为相场有效Kirchhoff应力张量/>
Figure FDA00041886990100000712
的分量，则其体积不变量P和偏量不变量Q分别为
Figure FDA00041886990100000713
其中，
Figure FDA00041886990100000714
为偏量应力张量，I为二阶单位张量，trace(·)表示张量的迹算子；
定义Drucker-Prager线性屈服面函数及其摩擦硬化法则分别为
Figure FDA00041886990100000715
其中，γ表示累积塑性乘子；
定义Cap椭圆屈服面函数及其压缩硬化法则分别为
Figure FDA0004188699010000081
其中，Pi为椭圆屈服函数中心，
Figure FDA0004188699010000082
为塑性体量对数应变；则两个屈服面满足光滑连接的相切条件为/>
Figure FDA0004188699010000083
切点为/>
Figure FDA0004188699010000084
此外，针对椭圆关联返回映射区域采用塑性关联流动法则，则采用Cap椭圆屈服面函数
Figure FDA0004188699010000085
作为其塑性势函数，而对于线性非关联返回映射区域(2)和椭圆非关联返回映射区域(3)采用塑性非关联流动法则，则其塑性势函数定义如下
Figure FDA0004188699010000086
其中，存在满足
Figure FDA0004188699010000087
并且两个塑性势函数的切点为/>
Figure FDA0004188699010000088
3.根据权利要求2所述的动态冲击/接触弹塑性大变形断裂分析显式相场物质点法，其特征在于，所述耦合显式相场-塑性模型的具体实施流程如下：
步骤1、弹性变形预测：
Figure FDA0004188699010000089
其中，
Figure FDA00041886990100000810
表示相对变形梯度张量，Δu表示增量位移向量；
步骤2、应变谱分解：beTr＝n·λeTr·nT；
其中，λeTr和n分别表示预测弹性左Cauchy-Green应变张量beTr的特征值和特征向量；
步骤3、计算预测应变不变量
Figure FDA00041886990100000811
和/>
Figure FDA00041886990100000812
/>
其中，预测弹性对数应变向量
Figure FDA00041886990100000813
其体积不变量和偏量不变量分别为
Figure FDA00041886990100000814
和/>
Figure FDA00041886990100000815
其中预测弹性对数应变向量的偏量应变向量
Figure FDA00041886990100000816
Figure FDA00041886990100000817
和单位向量δ＝(1，1，1)；
步骤4、根据公式(56)和(57)预测弹性应力不变量PTr和QTr；
步骤5、根据公式(59)检查：
Figure FDA00041886990100000818
步骤5.1当
Figure FDA00041886990100000819
根据公式(58)检查/>
Figure FDA00041886990100000820
步骤5.1.1当
Figure FDA00041886990100000821
进入线性非关联返回映射区域；
步骤5.1.1.1采用线性非关联返回映射算法：计算
Figure FDA00041886990100000822
和Δγ，并根据公式(58)更新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 FDA0004188699010000091
检杳应力状态是否处于弹性区域PTr≥P*：
步骤5.1.2.1当PTr≥P*：进入弹性变形步骤7；
步骤5.1.2.2当pTr＜P*：进入步骤6；
步骤5.2当
Figure FDA0004188699010000092
进入弹性变形步骤7；
步骤6对每一个迭代步检查应力状态是否在椭圆非关联返回映射区域或椭圆关联返回映射区域：
Figure FDA0004188699010000093
步骤6.1当
Figure FDA0004188699010000094
进入椭圆非关联返回映射区域；
步骤6.1.1采用椭圆非关联返回映射算法：计算
Figure FDA0004188699010000095
和Δγ，根据公式(58)更新M，并检查/>
Figure FDA0004188699010000096
步骤6.1.1.1当
Figure FDA0004188699010000097
根据公式(59)更新Pi，并进入步骤8；
步骤6.1.1.2当
Figure FDA0004188699010000098
Cap椭圆压缩塑性模型失效，进入纯剪切Drucker-Prager塑性模型：
步骤6.2当
Figure FDA0004188699010000099
进入椭圆关联返回映射区域；
步骤6.2.1采用椭圆关联返回映射算法：计算
Figure FDA00041886990100000910
和Δγ，根据公式(58)更新M，并检查/>
Figure FDA00041886990100000911
步骤6.2.1.1当
Figure FDA00041886990100000912
根据公式(59)更新Pi，并进入步骤8；
步骤6.2.1.2当
Figure FDA00041886990100000913
Cap椭圆压缩塑性模型失效，进入纯剪切Drucker-Prager塑性模型：
步骤7弹性变形：设定M＝Mn、Pi＝Pi，n、
Figure FDA00041886990100000914
和(P，Q)＝(PTr，QTr)
步骤8计算弹性主对数应变向量εe和Kirchhoff应力张量的分量τA；
步骤9计算弹性左Cauchy-Green应变张量be＝n·exp(2εe)·nT，并根据公式(55)更新Kirchhoff应力张量τ。
CN202210978276.1A 2022-08-16 2022-08-16 动态冲击/接触弹塑性大变形断裂分析显式相场物质点法 Active CN115410663B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210978276.1A CN115410663B (zh) 2022-08-16 2022-08-16 动态冲击/接触弹塑性大变形断裂分析显式相场物质点法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210978276.1A CN115410663B (zh) 2022-08-16 2022-08-16 动态冲击/接触弹塑性大变形断裂分析显式相场物质点法

Publications (2)

Publication Number Publication Date
CN115410663A CN115410663A (zh) 2022-11-29
CN115410663B true CN115410663B (zh) 2023-06-02

Family

ID=84158834

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210978276.1A Active CN115410663B (zh) 2022-08-16 2022-08-16 动态冲击/接触弹塑性大变形断裂分析显式相场物质点法

Country Status (1)

Country Link
CN (1) CN115410663B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116486953B (zh) * 2023-04-23 2023-09-05 哈尔滨工业大学 一种包含微结构效应的断裂相场仿真方法
CN116358439B (zh) * 2023-06-02 2023-08-25 山东省地质科学研究院 一种岩石有限应变测量方法、系统、电子设备及存储介质
CN117034689B (zh) * 2023-08-02 2024-02-02 大连理工大学 一种基于无网格rbf映射技术的土体液化大变形分析方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112036060A (zh) * 2020-08-03 2020-12-04 武汉大学 一种用于模拟脆性材料破坏的双线性自适应相场方法
CN112051142A (zh) * 2020-08-03 2020-12-08 武汉大学 一种用于模拟脆性材料不同破坏模式的普适性相场方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110298105B (zh) * 2019-06-26 2021-04-16 大连理工大学 饱和多孔介质大变形分析的ccpdi-impm方法
CN111832102B (zh) * 2020-06-19 2024-02-06 浙江大学 一种高维随机场条件下的新型复合材料结构优化设计方法
CN112652060A (zh) * 2021-01-06 2021-04-13 上海交通大学 一种基于粒子图像测速法的多模态视触觉传感系统及方法
CN113360992B (zh) * 2021-06-29 2022-02-15 大连理工大学 岩土结构大变形断裂分析的相场物质点方法
CN114186456B (zh) * 2021-12-02 2022-06-14 大连理工大学 结构冲击弹塑性断裂分析的时间间断态基近场动力学方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112036060A (zh) * 2020-08-03 2020-12-04 武汉大学 一种用于模拟脆性材料破坏的双线性自适应相场方法
CN112051142A (zh) * 2020-08-03 2020-12-08 武汉大学 一种用于模拟脆性材料不同破坏模式的普适性相场方法

Also Published As

Publication number Publication date
CN115410663A (zh) 2022-11-29

Similar Documents

Publication Publication Date Title
Ni et al. Hybrid FEM and peridynamic simulation of hydraulic fracture propagation in saturated porous media
Zhou et al. Phase-field modeling of fluid-driven dynamic cracking in porous media
CN115410663B (zh) 动态冲击/接触弹塑性大变形断裂分析显式相场物质点法
Zhang et al. A fictitious crack XFEM with two new solution algorithms for cohesive crack growth modeling in concrete structures
Silling et al. Peridynamic states and constitutive modeling
CN113360992B (zh) 岩土结构大变形断裂分析的相场物质点方法
Kafaji Formulation of a dynamic material point method (MPM) for geomechanical problems
Geiger et al. Black-oil simulations for three-component, three-phase flow in fractured porous media
Lian et al. Coupling of finite element method with material point method by local multi-mesh contact method
Foster et al. Embedded strong discontinuity finite elements for fractured geomaterials with variable friction
Rundle A physical model for earthquakes: 3. Thermodynamical approach and its relation to nonclassical theories of nucleation
Saksala et al. Combined continuum damage‐embedded discontinuity model for explicit dynamic fracture analyses of quasi‐brittle materials
De Borst Some recent developments in computational modelling of concrete fracture
Santillán et al. Phase-field model for brittle fracture. Validation with experimental results and extension to dam engineering problems
Wang et al. A hybrid local/nonlocal continuum mechanics modeling and simulation of fracture in brittle materials
Sun et al. Nonlinear dynamic response and damage analysis of hydraulic arched tunnels subjected to P waves with arbitrary incoming angles
Yan et al. A three-dimensional thermal-hydro-mechanical coupling model for simulation of fracturing driven by multiphysics
Li et al. An implicit coupling finite element and peridynamic method for dynamic problems of solid mechanics with crack propagation
Fei et al. Phase‐field modeling of rock fractures with roughness
Khan Investigation of discontinuous deformation analysis for application in jointed rock masses
Liu et al. An extended finite element framework for slow‐rate frictional faulting with bulk plasticity and variable friction
Chau et al. Numerical analysis of flyer plate experiments in granite via the combined finite–discrete element method
Armaghani et al. Investigating the effect of jointed environment on the cracked concrete arch dam in 3D conditions using FEM
Duan et al. A dynamic phase field model for predicting rock fracture diversity under impact loading
Liu et al. Long-term stability analysis for high arch dam based on time-dependent deformation reinforcement theory

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