固体力学学报

求解三维固体结构裂纹扩展的有限质点法

PDF下载

杨超 ¹^,²^,^** , 张鹏飞 ³ , 汪伟 ¹ , 罗尧治 ¹

固体力学学报 | 研究论文 2025,46(1): 39-53

收起

固体力学学报 | 研究论文 2025, 46(1): 39-53

求解三维固体结构裂纹扩展的有限质点法

全屏

杨超¹^,²^,^**, 张鹏飞³, 汪伟¹, 罗尧治¹

作者信息

¹浙江大学建筑工程学院，杭州，310058

²浙江大学长三角智慧绿洲创新中心，嘉兴，314100

³中冶华天南京工程技术有限公司，南京，210019

通讯作者:

**E-mail：zj_yangchao@zju.edu.cn.

The Finite Particle Method for Solving Crack Propagation of Three-Dimensional Solids

Chao Yang¹^,²^,^**, Pengfei Zhang³, Wei Wang¹, Yaozhi Luo¹

Affiliations

¹College of Civil Engineering and Architecture, Zhejiang University, Hangzhou, 310058

²Innovation Center of Yangtze River Delta, Zhejiang University, Jiaxing, 314100

³MCC Huatian Engineering Technology Co. Ltd., Nanjing, 210019

出版时间: 2025-02-28 doi: 10.19636/j.cnki.cjsm42-1250/o3.2024.046

文章导航

摘要

收起

对三维固体裂纹扩展的数值模拟一直是较为复杂的过程，裂纹路径的不可预知性使其计算求解相当繁琐且耗时. 有限质点法是以向量力学为基础的一种结构复杂行为分析的新方法，本文将它与粘聚力模型相结合用于三维固体动态断裂问题的分析. 该方法将固体域离散为质点群，允许通过自由增减质点实现从连续到非连续变形过程的模拟，处理裂纹扩展问题有独特优势. 本文基于外加粘聚力模型构建有限质点法三维固体断裂的判别准则和粘聚元自动嵌入的断裂分析模式，提出一种基于遍历搜索的质点拓扑更新的通用处理方法，利用CUDA开发工具包编制三维固体断裂计算的GPU并行求解程序，并通过若干算例探讨方法和程序的有效性和正确性.

关键词

三维固体 / 有限质点法 / 粘聚力模型 / 裂纹扩展 / 并行计算

Abstract

收起

Simulating three-dimensional (3D) crack propagation in solid structures poses significant challenges due to the unpredictability of crack paths, complicating both computation and solution strategies. Traditional methods often face difficulties in accurately capturing arbitrary crack propagation during large deformations. The finite particle method (FPM), based on vector mechanics, offers a novel numerical approach for analyzing complex behaviors in solid mechanics. Different from conventional continuum-based methods, FPM discretizes the solid domain into a collection of finite particles, each governed by Newton's second law of motion. This particle-based formulation enables seamless transitions between continuum and non-continuum behaviors by dynamically adding or removing particles, providing significant advantages for crack propagation analysis in both static and dynamic scenarios. In this study, the FPM is extended to address the dynamic fracture in 3D solids, focusing on the challenges related to crack initiation, propagation, and branching. The FPM is combined with an extrinsic cohesive zone model (CZM) to capture the complex behaviors of fractures, avoiding the need to pre-define crack paths and effectively managing discontinuities caused by crack propagation. A discriminant criterion is developed to identify the onset of crack initiation, and an automated embedding process for cohesive elements is implemented to enable real-time simulation of fracture surfaces. To manage the evolving topologies that arise from crack propagation, we propose a general strategy based on an ergodic search algorithm, which updates the connectivity of the discretized solid model dynamically as cracks evolve. In addition, we develop a GPU-based parallel solver using the CUDA toolkit to significantly accelerate fracture computations. The accuracy and applicability of the proposed method are validated through several numerical examples, including fracture simulations of plates and beams subjected to dynamic loading. The results demonstrate the capability of the method to accurately capture the intricate details of crack initiation, growth, and interaction in 3D solids. This extended FPM framework serves as a robust tool for analyzing dynamic fractures in engineering applications, providing a versatile framework for studying delamination, material failure, and structural collapse in both research and practical settings.

Key words

three-dimensional solid / finite particle method / cohesive zone model / crack propagation / parallel calculation

引用本文

杨超, 张鹏飞, 汪伟, 罗尧治. 求解三维固体结构裂纹扩展的有限质点法. 固体力学学报, 2025 , 46 (1) : 39 -53 . DOI: 10.19636/j.cnki.cjsm42-1250/o3.2024.046

Chao Yang, Pengfei Zhang, Wei Wang, Yaozhi Luo. The Finite Particle Method for Solving Crack Propagation of Three-Dimensional Solids[J]. Chinese Journal of Solid Mechanics, 2025 , 46 (1) : 39 -53 . DOI: 10.19636/j.cnki.cjsm42-1250/o3.2024.046

正文

收起

0　引言

收起

断裂是固体结构一种常见的失效模式，在各工程领域中经常发生源于断裂的灾难性破坏事故，如地震引起的地质构造开裂和建筑工程垮塌、交通运载工具损坏、压力管道开裂和机械构件破断等. 固体动态断裂理论虽已有长期发展，但在实际应用中仍然存在一定的局限性，尤其是三维大变形破坏问题，裂纹具有复杂的形状和任意扩展的路径，传统断裂力学中对裂纹是平直的假设不再成立. 对于裂纹在三维实体和曲面中任意扩展这类复杂断裂力学的研究往往需要借助数值方法. 目前常用的断裂力学计算方法有传统有限元与自适应网格^[1]、节点力释放^[2]、单元间内聚力模型^[3]、嵌入非连续模型^[4]等联合求解，以及离散元法、无网格法、扩展有限元法等. 有限元法是以函数连续体假设为理论基础，处理断裂不连续问题会遇到本质上的困难，需要借助特殊的模拟技术进行改进和修正^[5]；离散元法受模型几何的限制，适用于砂石等颗粒状介质的分析^[6]，且由于颗粒间隙易导致裂缝形态表征效果不佳；扩展有限元法和无网格法在一定程度上摆脱了网格依赖，但计算代价大，三维多裂纹分叉与相交求解也有待完善^[7,8]. 正是由于断裂面附近位移场、应力场的不连续性以及裂纹扩展导致的能量剧烈变化，给三维计算断裂力学带来了巨大挑战.

粘聚力模型在弹塑性断裂和界面力学领域应用非常广泛^[3]，它避免了线弹性断裂力学理论中的裂尖奇异性，可以防止断裂导致的应力与能量突变，相较于常用的J积分法不需要在材料或结构中预置裂纹，计算模型相对简单. 该模型最早由Barenblatt等提出^[9]，认为在裂纹尖端存在微小的粘聚区，该区域内应力为开裂位移的函数. 常用的粘聚力模型分为两种，即固有型粘聚力模型和外加型粘聚力模型^[10]. 对于前者，粘聚单元由于自建模初始就存在于结构分析域内，导致计算自由度大幅增加，其初始刚度罚值对整体模型刚度影响很大，尤其是对于显式计算会导致临界时间步长显著减小，计算量大幅增加；后者只在断裂条件激活的单元界面插入粘聚单元，不存在上述问题，但在断裂发生后需要动态更新单元和节点的拓扑关系，这是数值计算中非常棘手的问题. Pandolfi和Ortiz^[11]、Gürses和Miehe^[12]、Paulino等^[13]从点、线、面等角度出发提出了不同的计算模型数据存储结构，用来处理断裂过程中新生节点的归属以及单元拓扑连接关系改变问题，但对二维、三维以及高次单元分裂分离的处理方式不同，影响算法的通用性.

有别于传统分析力学架构下的计算方法，有限质点法是基于向量式力学理论提出的面向工程结构复杂行为的新型数值分析方法^[14]. 它不再以函数连续体的概念来描述分析对象，而是在空间上和时间上将求解域分别离散为相互联系的质点系统和途径单元的集合，利用牛顿定律取代连续体的偏微分方程描述质点运动，通过追踪质点系统的运动轨迹来表征结构体的运动和变形情况. 该方法中每一质点的控制方程相互独立，不涉及刚度矩阵的集成和求逆，也不需要迭代求解非线性方程组，同时摆脱了网格连续性条件限制，可以自由增减质点，方便实现连续至非连续变形过程的模拟，在结构强非线性问题分析中具有优势和灵活性. 目前该方法已应用于弹塑性大变形、结构屈曲、接触碰撞、倒塌破坏以及张力找形、折纸开合、流固耦合等结构复杂行为的分析^[15-19]，均取得了较好的效果.

本文通过理论推导和数值模拟，将有限质点法与外加型粘聚力模型相结合用于三维固体动态断裂问题的分析. 文中首先介绍三维固体的有限质点法基本理论和控制方程，并以缩减积分六面体实体元为例给出质点间作用力的推导；然后建立三棱柱形和六面体形粘聚元用于固体断裂面的描述，提出一种基于遍历搜索的质点间拓扑关系更新的通用处理方法；在此基础上，给出三维固体断裂计算的GPU并行算法实现流程；最后通过算例的模拟和分析，检验方法的正确性和有效性.

1　三维固体的有限质点法基本理论

收起

有限质点法采用标准化向量式力学分析架构，可同时求解结构体系受力后的位置变化与几何变形，其基本概念包括：1）点值描述；2）途径单元；3）刚体运动和纯变形分离的物理模式处理；4）基于牛顿第二定律的强形式广义质点运动控制方程. 关于有限质点法的基本思想在文献[14]和[18]中已有详细介绍，这里结合三维固体仅简单介绍其要点.

1.1　离散质点模型

根据三维固体的变形特点，有限质点法将求解域离散为一组具有空间平动自由度的质点集合，而质点间由一组能够发生三维变形的实体元相连. 固体域中的所有物理量均由各质点的点值描述，实体元没有质量且相互独立，仅用于表示质点间相互作用和联系，可通过内插函数来计算固体变形和内力. 实体元的界面连续性非必要条件，它们随质点运动发生变形，通过应力积分或者直接由经验公式求得内力与位移的关系，由此产生的内力反作用于相连质点，即被称之为“点值描述”.

如图1所示，假设某一固体在外力作用下经历一个连续的变形和运动过程，从初始时刻（t₀=0）构形V₀经中间时刻（t_a=t-τ）构形V_a至当前时刻t的构形V. 该时段内（t₀-t）分布于固体域内所有质点的运动轨迹被一组时间点（t₀，t₁，t₂，…，t）分割成若干独立的时间片段，每个时间段的运动轨迹即是一个途径单元，目的是简化内力计算和便于处理不连续行为^[18,19]. 假定断裂只发生在各途径单元的时间点上，途径单元内仍为连续的运动行为，在每个途径单元的起始点进行断裂判别，更新断裂之后的几何形态和相关质点属性（如质量、数量、材料、约束、拓扑等），再以新构形作为该途径单元内的基础架构，计算该时段内的质点轨迹. 断裂不会给计算分析框架带来本质上的改变，因而有限质点法处理此类强非线性问题具有独特的优势.

1.2　质点运动方程

从物理角度考虑，有限质点法认为所有质点均处于动平衡状态，在每个途径单元内（如t_a≤t≤t_b）的运动都遵循牛顿第二定律，且控制方程独立. 如图2所示，对于三维固体域内的任一质点α，运动变量可分解为沿坐标轴方向的3个平动自由度分量，对应坐标轴方向的3个质点力. 以

为参考，点α的全位移向量为

，途径单元t_a≤t≤t_b内质点运动方程式可表达为

(1)

式中，m_α为质点α的质量，包括该点的集中质量及相连实体元所提供的等效质量；

为质点α外力向量，包括质点受到的集中力以及由相连实体元上分布荷载转化来的等效外力，即

，n_c代表与质点α相连的实体元个数；

为质点α内力向量，可通过相连实体元的节点等效内力组集得到，即

；

为作用在质点α上的阻尼力向量^[20]；

为作用在质点α上的粘聚力向量，可通过相连粘聚元的节点等效作用力组集得到，即

，n_z代表与质点α相连的粘聚元个数.

有限质点法采用显式时间积分法求解运动方程，避免了复杂的迭代和收敛问题. 将质点的速度和加速度写成位移和时间的中央差分表达式，带入式（1）后即可得到途径单元内各时间步的质点位移、速度等物理量.

1.3　质点内力计算推导

以图3所示六面体实体元为例，在t_a和t时刻的构形分别用Ω_a和Ω表示. 节点（i=1，2，…，8）与质点（I=l，m，…，s，t）刚性联结，t_a和t时刻节点位置向量分别记作

和x_i=x_I. 实体元以t_a时刻实体元构形Ω_a为基础架构，t_a-t时段内节点位移向量为

固体质点内力仅与相连实体元的纯变形相关，所涉及到的应变、应力均通过节点的纯变形量计算，因此需要扣除Δu_i中的刚体运动分量，这是有限质点法求解大变形问题的关键. 为了估算实体元在t_a-t时段内的刚体转动量，将六面体顶面与底面之间的四边形中面M作为参考平面，参考平面的转动向量可假设为整个实体元刚体转动的近似值. 按照有限质点法虚拟逆向运动的基本思想，求解t_a-t时段内节点纯变形量

的主要步骤如下：令六面体元从构形Ω做逆向平移（-Δu₁）得到虚拟构形Ω′，然后依次发生面外逆向转动（-θ_op）和面内逆向转动（-θ_ip）得到虚拟构形

和Ω_r，构形Ω_r与Ω_a的相对位移即为所要估算的变形量.

与传统有限元法不同，这里求出的节点变量是变形位移而非全位移. 上述节点变形位移共有24个分量，去除对应于空间刚体运动模式的6个自由度后，应当只有18个独立变数. 在利用内插函数表示变形分布时，需要先求出这些独立量. 为此，引入变形坐标

，它与全域坐标的转换矩阵为

，可得到独立节点变形位移为

，其中

参照有限元法构建等参单元的方式，六面体实体元的内力计算基于等参元的形函数^[21]，其中任意一点的坐标以及变形量可由自然坐标（ξ，η，ζ）表示的等参形函数插值而来，即

以及

. 利用变形分布函数，t_a-t时段内的增量应变和增量应力可分别表示为

(2)

(3)

式中，

是将形函数求导后的变形矩阵考虑部分位移为零修正后得到的；

为节点非零独立变形位移向量；

为应力在

时的材料切线模量矩阵.

再由虚功原理可得t时刻的等效节点内力为

(4)

借助自然坐标下的选择数值积分^[22]，由式（4）可以得到节点非零独立变形位移所对应的等效内力分量，其余6个分量对应于被扣除的刚体自由度，可利用实体元在三维空间的静力平衡条件求得. 然后，通过正向运动与坐标变换，将等效节点内力

转到当前构形的全局坐标系下，最后反向组集到各相连质点上.

以上给出了有限质点法六面体实体元等效节点内力的主要推导过程，四面体实体元计算步骤类似，具体公式可参见文献[23]和[24].

2　三维固体的断裂计算

收起

2.1　离散质点模型

本文采用外加粘聚力模型^[25]发展有限质点法的三维固体断裂计算方法. 如图4所示，在相邻实体元的界面上建立局部坐标系（o-s，t，n），其中n轴沿着界面法向，t和n轴位于界面内. 将相关的实体元节点应力转换至局部坐标系并取平均值，得到该界面上的正应力和剪应力分别为（σ_n，τ_s，τ_t），根据二次应力判别准则，单元界面上的等效应力σ_eff可以表示为

(5)

式中，β为拉剪应力耦合系数.

当界面等效应力σ_eff达到材料破坏应力σ_c时，界面断裂条件被激活，裂纹开始萌生，六面体或三棱柱粘聚元将分别自动嵌入相应几何形状的实体元之间，以模拟两个相邻实体元界面的开裂过程，如图5所示.

2.2　质点分裂与拓扑更新

有限质点法计算过程中，对于外加粘聚力模型，在质点处于可分类状态之前时，粘聚元两侧的质点是共用的，如图5中的1_d和1_u共用质点1. 当界面断裂条件被激活且围绕质点的界面形成至少一个完整断裂面时，实体元会沿着若干嵌入粘聚元的界面发生分裂，同时裂纹前端质点将被复制，相应实体元的连接拓扑关系也将发生改变.

根据从属实体元的个数，质点间的连接面（对于平面固体，退化为连接线）可分为内部面（线）和外部面（线）两类. 前者从属于两个实体元，后者仅从属于一个实体元，分别如图6中虚线和实线所示. 其中，满足断裂条件的内部面（线）将自动变为外部面（线）. 于是，可定义如下质点分裂判定准则：当某质点所从属的内部面（线）满足断裂条件转化为外部面（线）之后，若该质点所连接的实体元集合无法通过公共面（边）两两相连形成连通域，则该质点将被断开并重新复制. 如果绕质点有N（N≥1）个独立连通域，则该质点需要被复制N-1次，新生成质点的坐标、约束、速度等参数的初值与被复制点在分裂前的数值相同，其质量、内力、外力则需要按照新的拓扑关系重新分配和集成.

为了便于描述和理解，以图6所示的线-面拓扑关系为例，内部有四个粘聚元C₁、C₂、C₃、C₄依次嵌入. 当内部线MN、IJ满足断裂条件而转为外部线后，C₁、C₂自动嵌入，此时质点M所连接的实体元连通域被打断，因此质点M将被复制并生成新的质点M′，而质点I和J虽然满足所属内部线断裂的条件，但在它们断开之后所连接的实体元仍能通过公共边两两相连形成连通域，所以质点不发生分裂. 之后，内部线NP、PQ相继达到断裂条件，质点N、P、Q先后发生复制和分离，裂缝沿粘聚元C₁、C₃、C₄向前扩展.

对于质点连通域数量，本文利用遍历搜索算法^[26]来判断质点从属实体元的连通性. 令任一质点所连接实体元的邻接矩阵为Ψ，其中Ψ（p，q）代表实体元p与实体元q是否相连，1表示相连，0表示不相连，显然邻接矩阵Ψ是对称的. 以图6所示质点M为例，其邻接矩阵中Ψ（1，2）=Ψ（2，1）=Ψ（2，3）=Ψ（3，2）=Ψ（4，5）=Ψ（5，4）=Ψ（5，6）=Ψ（6，5）=Ψ（6，7）=Ψ（7，6）=Ψ（7，8）=Ψ（8，7）=1，其余元素均为0. 调用遍历搜索算法，从任意实体元开始出发，因内部线MN达到断裂条件，实体元3和实体元4不再相连，质点M从属的实体元只有部分被访问遍历，此时将其余未被访问到的实体元相连接的质点替换为新生质点M′，则可完成新旧质点的单元归属问题. 当裂纹发生分叉、相交等复杂现象时，例如质点M相连接的内部线有多条达到断裂条件，如图7所示，质点M将被复制多次，其数学描述仅是质点M的邻接矩阵D_M发生改变，上述方法仍然适用.

需要指出的是，虽然三维固体裂纹体的描述较二维复杂，但本质上只是实体元由面元变成体元、共用线拓展成了共用面而稍有差异，对于被断裂面分割后的绕质点连通域数量判断仍可采用同样的搜索方法，因此上述点-线-面拓扑关系更新的处理方法同样适用于三维固体点-面-体关系的描述.

2.3　粘聚力计算

在开裂过程中，粘聚质点之间不会立即完全脱离相互作用，因为断裂界面上有粘聚力，它将随着预定义的张力-位移关系发生变化. 本文重点旨在发展断裂模拟算法，这里仅以常用的线性单调模型^[27]为例介绍粘聚力的计算.

为了综合考虑粘聚元两侧界面之间的法向张开量δ_n和切向滑移量δ_s、δ_t，定义等效位移δ_eff如下：

(6)

粘聚元的引入会产生抵抗裂纹法向张开和切线滑移的粘聚力. 若当前等效位移δ_eff大于历史最大位移

，即裂纹张开，则为加载情形，内聚力会随着位移的增大而单调递减直至降为0，软化区扩展，如图8所示，直到等效位移δ_eff达到其极限位移δ_c时，断裂能完全释放，粘聚质点对则完全脱开. 当粘聚元各质点对的等效位移量δ_eff均超过了δ_c时，该粘聚元视作完全失效破坏. 于是，粘聚元应力（T_n，T_s，T_t）可以表示为

(7)

(8)

式中，δ_c为等效位移的极限值；G_c为断裂释放能，数值上等于图8中三角形围成的阴影面积.

相反，若当前等效位移δ_eff小于历史最大位移

，即裂缝闭合，则为卸载情形. 与塑性变形类似，该过程是不可逆的，相应的粘聚元应力可以表示为

(9)

特别地，从式（5）可以看出，实体元界面受压且剪应力足够大时也会导致断裂，此时粘聚元的法向张开量应满足δ_n=0. 若δ_n<0，说明实体元之间发生穿透，显然没有物理意义. 最后，将粘聚力分别视作裂纹两侧界面上的分布面力，借助有限元面荷载等效外力计算的虚功方程，通过积分容易求得粘聚元应力

在局部坐标系（o-s，t，n）下各节点的等效作用力

（六面体，i=1～4；三棱柱，i=1～3），再将其转换到全局坐标系并组集到相连接的质点上.

3　GPU并行程序编制及实现

收起

根据以上理论公式推导和数值模拟策略，本文基于自研的有限质点法GPU并行计算平台，采用C++语言和CUDA工具包进一步编制了三维固体断裂分析程序. 有限质点法并行化平台的详细信息可参考文献[23]和[28].

图9给出了基于有限质点法的三维固体断裂并行计算流程框架. 在每个时间迭代步内，依次运行断裂、质点力、运动轨迹这三类求解器，具体计算步骤如下：

（1）建立模型，读入初始条件，包括质点位置、速度、连接以及材料、几何、约束等参数；

（2）按照第2节固体断裂计算的推导，建立当前步每个实体元的界面局部坐标系，再由式（5）计算界面等效应力，并判断是否达到断裂条件，若满足则按照2.2节进行质点复制与拓扑更新，然后由式（6）-式（7）计算等效粘聚力

；

（3）按照第1.3节的推导，依次计算当前时刻每个实体元的纯变形、应变、应力和等效节点内力

；

（4）集成当前时刻的质点内力

、外力

、阻尼力

和粘聚力

，代入式（1）求解质点运动方程，得到下一时刻质点位移、速度等物理量；

（5）判断是否收敛或者最后时间步，如果满足则退出循环，否则回到第（2）步继续下一时间步的计算.

上述计算过程中，除了断裂求解器中质点复制和拓扑更新需要借助CPU端完成外，各求解器的主要计算任务均在GPU线程内并行执行. 在CUDA各线程执行过程中，对GPU内存的访问效率会直接影响整个并行求解器的计算效率，其中存储容量最大的是全局内存，它可以同时被所有线程访问，对整体计算效率的影响也最大. 程序中采用SoA数据存储方式^[29]，通过合并和对齐访问提高全局吞吐量，实现内存访问优化.

4　数值算例

收起

本节通过若干典型三维固体的裂纹扩展算例，与相关文献及软件计算结果进行对比分析，验证上述理论和所编制程序的有效性和正确性.

4.1　边缘单裂缝板片受拉

受动力拉伸作用的带缺陷矩形板片裂纹扩展是固体断裂计算的经典算例^[30,31]. 板片宽为W，高为L，厚度为h，左侧边缘中部有条长度为a的水平初始裂缝. 板片上、下端以恒定速度u均匀拉伸，由于对称性，裂缝会沿水平方向扩展直至板片完全撕裂. 几何尺寸、材料性质以及边界条件如图10所示. 为检验方法的收敛性，采用3种不同数量的质点离散形式，沿长宽厚方向质点分布分别为125×55×2、205×85×3、295×125×4、351×145×5，并分别用六面体实体元（EPM-H8）和四面体实体元（EPM-T4）进行连接，相应的计算模型记作EPM-T①～T④及EPM-H①～H④.

本算例是典型的板片I型断裂问题. 图11给出了8种离散模型计算得到的水平裂纹扩展长度随拉伸位移的变化情况. 可以看出，质点间距的减小导致曲线斜率增大，裂纹扩展速度加快，当质点间距小于1.02 mm时（EPM-T③～④及EPM-H③～④），裂纹扩展速度基本保持不变. 不同质点数模型的计算结果与文献[31]相比，平均误差分别为12.8%、9.8%、4.1%、1.7%（T①～T④）和11.9%、9.0%、3.5%、1.2%（H①～H④）. 随着模型细化，模拟结果趋于吻合，体现了本文方法良好的收敛性. 图12还给出了由EPM-H4模型计算得到的加载过程中板片端部约束力与位移关系曲线，与文献[30]和[31]以及ABAQUS中采用固有型粘聚力模型的分析结果基本一致. 端部约束力的波动反映出板片内应力波的传递影响，当反射在水平对称轴上时，导致了裂缝的扩展. 图13给出了拉伸过程中不同时刻的板片应力σ_yy的分布云图. 板片约在102 μs时裂纹开始向右端扩展，裂尖处出现了较明显的应力集中，并向周围逐渐扩散，直至约t=172 μs裂缝到达板片的右端. 本算例验证了本文方法在求解固体I型断裂问题方面的有效性与准确性.

4.2　边缘双缺口板片受撞击

本算例分析单边双缺口板片在局部低速冲击作用下的脆性破坏过程. Kalthoff等^[32]曾对此开展过试验研究，其结果被许多学者用来验证各类数值计算理论的可靠性^{[13,25,33,34]}. 板片的几何尺寸、材料性质及初始冲击作用如图14（a）所示，两条长度50 mm的初始裂缝对称分布在板片中间上下两侧，撞击物直径与初始裂缝间距相等，板片边缘施加初始冲击速度16.54 m/s来模拟撞击作用. 利用对称性，仅取上半部分建模，并约束水平对称面的竖向位移. 计算模型采用51×101×5个质点进行离散，再将缺口前端裂纹扩展区域（图14（b）阴影部分）加密为65×95×9个离散质点，使得裂纹区域的质点间距不大于最小粘聚区尺寸^[31]，质点间用EPM-H8实体元连接.

图15给出了板片受撞击后不同时刻的变形和裂纹扩展情况. 其中，大约在t=25 μs时，裂纹从缺口尖端开始扩展，与初始裂缝间夹角为68.1°，之后倾斜向上延伸，倾斜角度略有减小，直至约t=81 μs时到达板片上端. 图16给出了板片内的最终裂纹分布，与文献[25]、[32]和[33]的分析结果基本一致，且裂缝起始扩展角与Kalthoff试验结果（70°）误差仅为2.7%，表明本文方法能有效追踪冲击荷载下的裂纹扩展路径. 此外，本算例使用GPU并行的有限质点法计算耗时仅需7.8 s，与文献[33]和[25]分别采用自适应网格和非连续Galerkin法的计算耗时25 min和108 min相比，虽然存在计算硬件和离散形式等方面的差异，但仍在一定程度体现了本文方法和模拟技术对固体断裂问题并行求解的高效性.

4.3　矩形梁受拉扭大变形

本算例对矩形梁受扭拉断裂行为进行模拟. 如图17所示的矩形截面梁，长度为L，横截面边长为a. 为模拟其在受拉扭动力荷载作用下的破坏，梁两端分别以恒定的轴向拉伸速度u和绕轴转动速度ω进行加载. 几何尺寸、材料性质以及相关计算参数如图17所示. 计算模型中共布置31×31×151个质点对矩形梁离散，质点间距满足粘聚区尺寸要求，质点间由四面体实体元（EPM-T4）连接.

矩形梁受拉扭裂缝是Ⅰ、Ⅱ、Ⅲ型断裂模式的组合. 图18比较了本文方法与ABAQUS采用粘聚力模型计算得到的加载过程中梁端荷载时程，两种方法得到的拉力和扭矩变化趋势基本一致，本文方法能够较好地捕捉到裂缝扩展阶段因断裂能释放引起的端部扭矩值波动. 图19给出了不同加载时刻矩形梁的变形及开裂情况，追踪到裂纹萌生、扩展、相交、分叉等复杂行为以及从大变形到完全破坏的完整过程，在矩形梁中部裂缝由最初萌生的四条细小裂纹逐渐扩展、相交演变成一个环形裂尖，并最终贯穿梁体. 整个计算过程中没有出现因强非线性行为而导致的数值不稳定现象，反映出本文方法处理三维固体复杂断裂问题的良好性能.

5　结论

收起

本文将有限质点法与外加型粘聚力模型结合，应用于求解三维固体的裂纹扩展问题. 建立了三棱柱形和六面体形粘聚元模型用于固体断裂面的描述，发展了采用二次等效应力作为界面开裂和粘聚元自动插入的判别准则，推导了基于线性张拉—位移关系的加、卸载等效粘聚力计算公式，可避免裂纹尖端应力奇异性，并提出了一种基于遍历搜索的适合处理多种维度固体模型拓扑关系变化的通用断裂分析模式，而无需在材料或结构中预设裂纹路径. 结合GUDA架构实现了固体结构断裂计算的并行加速求解，拓展了有限质点法仿真平台断裂模拟功能模块. 通过一系列数值算例对本文方法进行了测试，并将计算结果和相关文献及通用软件结果进行了对比，表明该方法从理论到程序均能有效求解三维固体的裂纹扩展问题，且计算高效、精度较高.

有限质点法直接使用点值作为变量，固体发生断裂仅改变了质点上的作用力和质量，而不会给整体求解框架带来本质上的改变，质点行为仍受原运动方程式的控制，计算中始终维持质点上的能量守恒，而且每个质点的运动参数都是独立求解，不需要集成整体刚度矩阵，在满足结构体系质量守恒的条件下不受连续性条件的限制而可以自由增减质点，物理概念清晰，程序通用性较好，用于分析固体的断裂行为有其独特的优势和灵活性. 本文工作也为进一步研究固体材料脱层、碎裂及金属结构疲劳断裂等更为复杂的不连续行为打下了良好基础.

基金

收起

国家自然科学基金项目(52238001; 52008366)
国家重点研发计划资助项目(2021YEE0501000)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

Miehe

, Gürses

. A robust algorithm for configurational-force-driven brittle crack propagation with R-adaptive mesh alignment[J]. International Journal for Numerical Methods in Engineering, 2007, 72(2): 127-155.

[2]

Zhuang

, O'donoghue

P E

. The recent development of analysis methodology for rapid crack propagation and arrest in gas pipelines[J]. International Journal of Fracture, 2000, 101(3): 269-290.

[3]

Tabiei

, Meng

. Improved cohesive zone model: integrating strain rate, plastic strain, variable damping, and enhanced constitutive law for fracture propagation[J]. International Journal of Fracture, 2023, 244(1): 125-148.

[4]

Jirásek

. Comparative study on finite elements with embedded discontinuities[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 188(1): 307-330.

[5]

Ting

E C

, Shih

, Wang

Y K

. Fundamentals of a vector form intrinsic finite element: Part I. Basic procedure and a plane frame element[J]. Journal of Mechanics, 2004, 20(2): 113-122.

[6]

Tavarez

F A

, Plesha

M E

. Discrete element method for modelling solid and particulate materials[J]. International Journal for Numerical Methods in Engineering, 2007, 79(4): 379-404.

[7]

Belytschko

, Gracie

. A review of extended/generalized finite element methods for material modeling[J]. Modelling and Simulation in Materials Science and Engineering, 2009, 17(4): 43001.

[8]

杨建军, 文丕华. 无网格法：理论与算法[M]. 北京: 科学出版社, 2018.

(Yang

J J

, Wen

P H

. Meshless Method: Theory and Approches[M]. Beijing: Science Press, 2018. (in Chinese))

[9]

Barenblatt

G I

. The Mathematical Theory of Equilibrium Cracks in Brittle Fracture[M]//Dryden

H L

, von Kármán

. Advances in Applied Mechanics. New York: Academic Press, 1962, 7: 55-129.

[10]

Kubair

D V

, Geubelle

P H

. Comparative analysis of extrinsic and intrinsic cohesive models of dynamic fracture[J]. International Journal of Solids and Structures, 2003, 40(15): 3853-3868.

[11]

Pandolfi

, Ortiz

. An efficient adaptive procedure for three-dimensional fragmentation simulations[J]. Engineering with Computers, 2002, 18(2): 148-159.

[12]

Gürses

, Miehe

. A computational framework of three-dimensional configurational-force-driven brittle crack propagation[J]. Computer Methods in Applied Mechanics and Engineering, 2009, 198(15): 1413-1428.

[13]

Paulino

G H

, Celes

, Espinha

, Zhang

Z J

. A general topology-based framework for adaptive insertion of cohesive elements in finite element meshes[J]. Engineering with Computers, 2008, 24(1): 59-78.

[14]

喻莹. 基于有限质点法的空间钢结构连续倒塌破坏研究[D]. 杭州: 浙江大学, 2010.

(Yu

. Progressive Collapse of Space Steel Structures Based on the Finite Particle Method[D]. Hangzhou: Zhejiang University, 2010. (in Chinese))

[15]

罗尧治, 郑延丰, 杨超, 喻莹, 俞锋, 张鹏飞. 结构复杂行为分析的有限质点法研究综述[J]. 工程力学, 2014, 31(8): 1-7, 23.

(Luo

Y Z

, Zheng

Y F

, Yang

, Yu

, Zhang

P F

. Review of the finite particle method for complex behaviors of structures[J]. Engineering Mechanics, 2014, 31(8): 1-7, 23. (in Chinese))

[16]

Zheng

, Wu

H W

, Li

S Y

, Zhang

J Y

, Yang

, Luo

Y Z

. Analysis of dynamic process of thick-panel origami structures based on finite particle method[J]. Engineering Structures, 2023, 295: 116761.

[17]

Liu

F H

, Yu

, Wang

Q H

, Luo

Y Z

. A coupled smoothed particle hydrodynamic and finite particle method: An efficient approach for fluid-solid interaction problems involving free-surface flow and solid failure[J]. Engineering Analysis with Boundary Elements, 2020, 118: 143-155.

[18]

喻莹, 罗尧治. 基于有限质点法的结构倒塌破坏研究Ⅰ：基本方法. 建筑结构学报, 2011, 32(11): 17-26.

(Yu

, Luo

Y Z

. Structural collapse analysis based on finite particle method I: Basic approach[J]. Journal of Building Structures, 2011, 32(11): 17-26. (in Chinese))

[19]

, Luo

Y Z

, Paulino

G H

. Finite Particle Method for Progressive Failure Simulation of Truss Structures[J]. Journal of Structural Engineering, 2011, 137(10): 1168-1181.

[20]

喻莹, 刘飞鸿, 王钦华, 罗尧治, 李洋. 有限质点法阻尼构造问题的研究[J]. 工程力学, 2019, 36(11): 34-40.

(Yu

, Liu

F H

, Wang

Q H

, Luo

Y Z

, Li

. Study on damping in finite particle method[J]. Engineering Mechanics, 2019, 36(11): 34-40. (in Chinese))

[21]

朱伯芳. 有限单元法原理与应用[M]. 北京: 中国水利水电出版社, 2018.

(Zhu

B F

. The Finite Element Method Theory and Applications[M]. Beijing: China Water Power Press, 2018. (in Chinese))

[22]

Hughes

T J R

. Generalization of selective integration procedures to anisotropic and nonlinear media[J]. International Journal for Numerical Methods in Engineering, 1980, 15(9): 1413-1418.

[23]

唐敬哲. 面向结构-土协同效应分析的有限质点法并行计算平台研发[D]. 杭州: 浙江大学, 2021.

(Tang

J Z

. Parallelized Finite Particle Method Computational Platform for Soil-Structure Interaction Analysis[D]. Hangzhou: Zhejiang University, 2021. (in Chinese))

[24]

张鹏飞, 罗尧治, 杨超. 基于有限质点法的三维固体弹塑性问题求解[J]. 工程力学, 2017, 34(4): 5-12.

(Zhang

P F

, Luo

Y Z

, Yang

. Elastic-plastic analysis of 3D solids using the finite particle method[J]. Engineering Mechanics, 2017, 34(4): 5-12. (in Chinese))

[25]

Nguyen

V P

. Discontinuous Galerkin/extrinsic cohesive zone modeling: Implementation caveats and applications in computational fracture mechanics[J]. Engineering Fracture Mechanics, 2014, 128: 37-68.

[26]

Thomas

H C

, Charles

E L

, Ronald

L R

, Clifford

. Introduction to Algorithms[M]. Fourth Edition ed. Cambridge: MIT Press, 2022.

[27]

Park

, Paulino

G H

. Cohesive zone models: A critical review of traction-separation relationships across fracture surfaces[J]. Applied Mechanics Reviews, 2011, 64(6): 60802.

[28]

Tang

J Z

, Zheng

Y F

, Yang

, Wang

, Luo

Y Z

. Parallelized implementation of the finite particle method for explicit dynamics in GPU[J]. Computer Modeling in Engineering & Sciences, 2020, 122(1): 5-31.

[29]

Cheng

, Grossman

, Mckercher

. Professional CUDA C Programming[M]. Indianapolis: John Wiley & Sons, 2014.

[30]

Becker

, Geuzaine

, Noels

. A one field full discontinuous Galerkin method for Kirchhoff-Love shells applied to fracture mechanics[J]. Computer Methods in Applied Mechanics and Engineering, 2011, 200(45): 3223-3241.

[31]

Zavattieri

P D

. Modeling of crack propagation in thin-walled structures using a cohesive model for shell elements[J]. Journal of Applied Mechanics, 2006, 73(6): 948-958.

[32]

Kalthoff

J F

. Shadow optical analysis of dynamic shear fracture[J]. Optical Engineering, 1988, 27(10): 271035.

[33]

Park

, Paulino

G H

, Celes

, Espinha

. Adaptive mesh refinement and coarsening for cohesive zone modeling of dynamic fracture[J]. International Journal for Numerical Methods in Engineering, 2012, 92(1): 1-35.

[34]

Zhou

, Jia

. Dynamic propagation of moving cracks in brittle materials by field-enriched finite element method[J]. Engineering Fracture Mechanics, 2024, 305: 110177.

2025年第46卷第1期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.19636/j.cnki.cjsm42-1250/o3.2024.046

接收时间：2024-10-10
首发时间：2026-03-20
出版时间：2025-02-28

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2024-10-10

基金

国家自然科学基金项目(52238001; 52008366)

国家重点研发计划资助项目(2021YEE0501000)

作者信息

¹浙江大学建筑工程学院，杭州，310058

²浙江大学长三角智慧绿洲创新中心，嘉兴，314100

³中冶华天南京工程技术有限公司，南京，210019

通讯作者:

**E-mail：zj_yangchao@zju.edu.cn.

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/gtlxxb/CN/10.19636/j.cnki.cjsm42-1250/o3.2024.046

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT