Finite element implementation of the coupled criterion for numerical simulations of crack initiation and propagation in brittle or quasi-brittle materials

Jia Li^a*, Dominique Leguillon^b

^a LSPM, CNRS UPR 3407, Institut Galileacute;e, Universiteacute; Sorbonne Paris Citeacute;, Universiteacute; Paris 13, F-93430 Villetaneuse, France

^b Institut Jean Le Rond dAlembert, CNRS UMR 7190, Sorbonne Universiteacute;s, UPMC Universiteacute; Paris 06, F-75005 Paris, France.

*: Corresponding author

LSPM, CNRS UPR 3407

Institut Galileacute;e, Universiteacute; Paris 13

99 Avenue Jean-Baptiste Cleacute;ment,

93430 Villetaneuse, France

Tel : 33 1 49 40 28 89 Fax : 33 1 49 40 39 38

E-mail : jia.li@univ-paris13.fr

Abstract

In this paper, we present a numerical model for simulating two-dimensional cracking process in brittle or quasi-brittle materials. This model is based on the coupled criterion which enters in the general framework of Finite Fracture Mechanics. It states that both the energy and stress conditions should be verified for the crack creation in solids. Special approaches were developed in drawing correctly the potential crack paths and to minimize the potential energy with respect to all predicted crack paths. This numerical model can be used to predict both the crack initiation and crack growth in complicated two-dimensional structures. The multiple cracking processes can also be simulated by using this model and statistical distributions of flaws can be accounted for. The efficiency and accuracy of the proposed approach were approved by means of several numerical examples for which the solutions are well-known or well-studied in the literature. Moreover, it was also shown that the numerical results are mesh-independent when the stress field is correctly evaluated.

Keywords: Crack initiation, crack growth, multiple cracking, finite element analysis, brittle or quasi-brittle materials

1. Introduction

Numerical prediction of fracture in solids is a classical challenge for mechanics and computational sciences. Developing efficient and accurate numerical models to simulate fracture behaviour has been extensively studied in the last two decades. On the one hand, many fracture models were established on the basis of continuum damage mechanics [1-5]. On the other hand, the fracture-mechanics-based models such as the cohesive model [6-8], the extended-finite element method (XFEM) [9-11], or the non-local fracture models [12-14] have also been largely investigated. Even though huge success has been achieved with all these fracture models, numeral drawbacks or restrictions still exist in various mechanical or computational aspects. Developing new efficient numerical models is always an actual need for fracture prediction and description.

In studying the crack onset from a sharp notch point, Leguillon [15] noticed that application of a single conventional fracture criterion, namely the Griffith energy criterion or the stress criterion, is not suitable and often leads to paradox predictions. It was concluded that when fracture occurs both these criteria are fulfilled simultaneously, even if one often hides the other. Both the criteria are necessary conditions and together they form a sufficient one. Application of this concept on crack onset predictions at a V-notch point shows remarkably satisfactory agreement with respect to experimental results. Afterward, this concept gets increasing popularity in fracture prediction of brittle or quasi-brittle materials. A large number of fracture cases were studied by using this criterion [16-27 among others].

Even though the accuracy and efficiency of this criterion were approved by means of numerous real structure analyses, these studies were often performed case by case by using particular resolution techniques such as asymptotic analysis, local finite element modelling etc. Therefore, we believe that it is necessary and useful to implement this criterion into a numerical code such that more complicated fracture problems can be dealt with. This is the main motivation of the present work.

In the present work, various numerical efforts were made for developing algorithms to predict automatically the crack paths. By means of these efforts, the Leguillon criterion was successfully implemented into a finite element code with which two-dimensional fractures in brittle or quasi-brittle materials can readily be simulated. Moreover, several classical fracture problems were studied in order to assess the efficiency of the proposed numerical model. The numerical results were compared with some well-known solutions in the literature. These comparisons show that the numerical results issued from the proposed numerical method are fully satisfactory.

2. Description of the numerical fracture model

The coupled criterion is based on the Finite Fracture Mechanics (FFM) in which the Griffith criterion is written in incremental form:

where is the change in potential energy, is the newly created crack surface, G is the so-called incremental energy release rate and is its critical value representing the fracture toughness of the material. This criterion is not sufficient since it depends on the newly created surface. For example, this criterion is difficult to fulfill when the stress concentration is weak and when is small. In order to remedy this shortcoming, the stress criterions have to be added. For brittle or quasi-brittle materials, the maximal principal stress criterion is an appropriate candidate:

where and are respectively the maximal principal stress and the ultimate stress of the material. The above condition must be fulfilled on the whole surface S. Combination of these two criteria allows the determination of

剩余内容已隐藏，支付完成后下载完整资料

---

基于有限元方法实现脆性或准脆性材料中裂缝起裂和扩展的数值模拟的耦合标准

李佳^a*, Dominique Leguillon^b

^a LSPM，CNRS UPR3407，伽利略研究所，巴黎太阳城索邦大学，巴黎大学13，F-93430维尔塔纳斯，法国

^b让·勒朗·达朗贝尔学院，CNRS UMR7190，索邦大学，UPMC巴黎大学06，F-75005法国巴黎。

*：通讯作者

LSPM，CNRS UPR3407

伽利略研究所，巴黎大学13

第9大道让 - 巴蒂斯特克莱门特，

93430维尔塔纳斯，法国

联系电话：331 4940 2889传真：331 4940 3938

电子信箱：jia.li@univ-paris13.fr

摘 要

本文提出了一种用于模拟脆性或准脆性材料中裂缝形成过程的二维数值模型。该模型基于有限断裂力学一般框架下的耦合标准。本文表明，断裂能和应力两个条件都能作为固体中裂缝产生的验证依据。开发出正确绘制潜在裂缝路径并降低所有预测裂缝路径的势能的特殊方法。该数值模型不仅可用于预测复杂二维结构中裂缝的起裂和扩展，还可用于模拟多个开裂过程，并且可以考虑裂缝的统计分布。本文通过几个数值示例来验证所提出的方法的高效和准确性，所采用解决方案是公认经过充分研究的。此外，本文还表明，当正确评估应力场时，数值结果与网格无关。

关键词：裂缝起裂，裂缝扩展，多裂缝，有限元分析，脆性或准脆性材料

简介

对固体中裂缝进行数值分析，是断裂力学和计算科学的经典挑战。在过去二十年中，已经广泛研究并开发出高效和精确的数值模型来模拟断裂行为。一方面，许多断裂模型建立在连续损伤力学的基础上^[1-5]，另一方面，基于断裂力学的模型如粘性模型^[6-8]，扩展有限元法（XFEM）^[9-11]，非局部断裂模型^[12-14]，都被大量研究。尽管这些断裂力学模型已经取得了巨大的成功，但是在一些力学或计算方面仍然存在数字缺陷或限制。开发新的、高效的数值模型始终是对裂缝预测和描述的真实需求。

Leguillon在进行裂缝从尖锐切口起裂的研究中注意到，应用单一的传统断裂标准，即Griffith能量标准或应力标准是不合适的，并且经常导致悖论预测。结论是，当断裂发生时，这些标准应同时实现，即便其中一个经常被隐藏。这两个标准都是必要的条件，它们一起形成了一个充分的条件。这个概念在预测V形切口裂缝发生的应用中显示出与实验结果相当令人满意的一致性。之后，这个概念在脆性或准脆性材料的断裂预测中越来越受欢迎。大量的断裂示都例使用了Leguillon准则^[16-27]。

虽然许多真实结构的研究分析证明了Leguillon准则的准确性和高效性，但是这些研究通常通过使用了特定的分析技术，如渐近分析，局部有限元建模等。因此，我们认为将此准则落实为数值代码是必要和有用的，以便可以处理更复杂的断裂问题。这是目前本工作的主要动机。

在本工作中，我们进行了大量数值分析工作，开发出自动预测裂缝路径的算法。通过这些努力，Leguillon准则被成功地落实为有限元代码，利用该有限元代码可以容易地模拟脆性或准脆性材料中的二维断裂。此外，为了评价提出的数值模型的效率，本文研究了几个经典断裂问题，并将数值结果与文献中的一些众所周知的解决方案进行比较。对比表明，本文提出的数值方法计算出结果是完全令人满意的。

数值断裂模型的描述

耦合的标准基于有限断裂力学（FFM），其中Griffith断裂判据以增量形式写为：

其中指势能的变化，指新形成的裂缝面面积，指能量释放率，指能量释放率的临界值，表示材料的断裂韧性指标。这个标准是不够的，因为它取决于新形成的裂缝面。例如，当应力集中弱且小时，该标准难以实现。为了弥补这一缺点，我们必须增加压力标准。对于脆性或准脆性材料，最大主应力判据是适当的选择：

其中和分别是材料的最大主应力和极限应力。上述条件必须在整个表面上满足。这两个标准的组合允许确定新创建的表面和断裂时的载荷。在平面应变框架内，Leguillon准则表明，当且仅当应力和能量条件都满足时，才产生新的裂纹，换句话说，沿着新创建的裂缝路径的所有点都满足应力准则，并且相应的裂缝路径满足Griffith 断裂判据^[15]。

Leguillon准则在概念上是清晰且有说服力的。其在各种学术或工程断裂问题中已经显示出显著的效率和准确性。当断裂位置和裂纹路径是预先已知的时，Leguillon准则应用容易。 但是，对于复杂结构，这些数据并不总是适用。因此，Leguillon准则的数值实现中的第一个任务就是确定裂缝开展的轨迹。为此，我们采用了一些简单但物理上合理的假设。下面，我们提出实现Leguillon准则的数值细节。本文采用三角形单元进行模拟。应当注意，由于解缺乏平滑性，因此在弹性奇点附近增加插值度是没有意义的。

形成潜在新裂缝的位置。

根据Leguillon准则，断裂的必要条件之一是假定的裂纹路径必须满足应力判据。因此，我们应该首先找到满足的所有区域。为了便于进一步简化算法我们提出补充假设，假设裂纹从第一主应力局部最大的点开始起裂。因此，我们的第一个任务是找到满足的区域，以及对于给定的规定荷载下的局部最大值。

为此，我们取反的符号，使得局部最大值成为局部最小值。然后我们使用所谓的水线阈值算法来找到这些区域。基本思想是考虑表面的应力分布，并将水源放置在每个区域最低的地势中。接下来，整个地势从水源处逐渐地充满，并且坝被放置在不同水源相遇处。图1显示了根据已建立算法模拟的1D表面的漫射过程^[28]。

大坝	大坝
(a)	(b) 大坝
(c)	(d)

图1：水线阈值算法的过程(a):水源放置在最低点处，增加水位允许找到相邻的单元，然后创建第一集水盆；(b) 通过增加水位，发现另一个区别的局部最小值，从该最小值产生第二集水盆地；(c)第一集和第二集水盆相遇的位置放置坝；(d)继续这一过程，直到，并且区别出含有每个局部最小值的集水盆。

围绕给定最小值的区域中的所有元素构成与该最小值相关联的集水盆地。分水线是划分相邻集水盆的区域。水上升达到的第一点就是处于最高应力值的位置，然后所有元素逐渐达到最低应力水平。

集水盆用整数值标记。每个集水盆包含局部最大应力，潜在的裂缝由此产生。

通常，经典水线阈值算法将导致输入应力分布的过分割，特别是对于具有许多区域最小值干扰的数值结果。在这种情况下，建议使用去噪方法预处理地形表面。在本工作中，使用中值滤波技术。中值滤波器的主要思想是对逐个元素绘制应力分布图，用相邻元素的值的中值代替元素的值。该算法允许提取N个区域，其中每个区域包含如“山峰”的局部最大第一主应力。

绘制潜在的裂纹路径

在每个区域中，潜在裂纹都可能产生。根据最大主应力准则，新产生的裂纹的方向垂直于第一主应力。这个标准从每个“山峰”处绘制主应力轨迹。为此，我们首先提出以下准则：

1. 沿着从“山峰”开始的主应力轨迹，主应力单调减小，即：

这里指源自“山峰”点的主应力轨迹的曲线坐标，指沿应力轨迹的增量；

1. 沿着从“山峰”开始的主应力轨迹，到山峰的距离单调增加。这个说法表明了裂纹不可逆。 准则如下：

这里，和是原点与点和之间的距离;

因此，我们构建以下算法来确定裂纹路径：

1. 从对应于中的“山峰”的点，我们先绘制预测裂纹路径的右分支。为此，我们计算坐标，其中是第一主应力轨迹在点的角度。如果需要，改变，使得点位于相邻元素中；
2. 如果点不在中，则转到步骤5。否则，从点以及在此点处第一主应力轨迹的角度开始，计算

且

如果需要，改变，使得点和位于相邻元素中；

1. 比较和，并选择对应于和的最大值的坐标。
2. 令并转到步骤2，直到不再位于中。
3. 从点，通过遵循与上述相同的算法绘制预测裂纹路径的左分支。

对于每个区域，绘制潜在的裂纹路径。可以通过累加，获得裂纹路径的长度。

最大能量释放率准则

显然，即使所有预测的裂纹路径都可以通过使用上述算法来绘制，但它们可能不会全都成为真正的裂纹。其中哪些变成真正的裂缝取决于能量标准的实现。显然，必须考虑可能裂纹扩展的所有组合，以便找到最大能量释放率。如果我们预测到N个潜在的裂纹路径，并且每个路径都可能变成真实的裂纹，则真实的裂纹可能有种组合。因此，必须执行种线性弹性有限元分辨率：一种用于裂纹扩展前的结构，表示所有可能的裂纹扩展组合的结构。这些组合中，使G值最大并验证Griffith标准的组合将成为真正的裂缝。

这种详尽的方法在数学上是严格的，但是当N大时可能导致大量的计算工作。在某些情况下，不同的潜在裂缝之间相互作用很小。因此，同时产生多个裂纹的能量释放率始终小于产生单个裂纹的最大能量释放速率。例如，令和是产生潜在裂纹和的能量释放率,，令 为同时产生裂纹和的能量释放率。当这两个潜在裂缝之间的相互作用较弱时，我们有。在这种情况下，我们只需要评估每个潜在裂纹的N个能量释放速率，以确定哪个裂纹成为真实的裂纹。在这种情况下，计算任务可以大大降低。

形成新裂缝

在有限元法的框架中存在形成新裂缝的不同技术。我们认为所提出的算法允许使用以下方法：

1. 沿着预测的裂缝路径去除元素^[14]；
2. 通过创建一个新的裂缝来弥补预测裂缝路径附近的元素^[29]；
3. 引入沿着预测的裂缝路径的不连续性^[10][30]。

很明显，在所有这些方法中，第一个是最简单的。由于这个原因，它用于本文中。实际上，我们发现与预测的裂纹路径相交的所有元素，然后从有限元模型中去除它们。

整体算法

为了避免复杂的快照分析并保证解决方案的存在，外部荷载的施加，通过在结构的边界上施加位移直到前置等级来逐步完成。通过使用上述技术，我们可以轻松地构建一个整体算法如下：

1. 在每个加载步骤中，进行第一有限元计算，以确定结构中的应力分布，然后计算总势能；
2. 得到的应力分布，我们首先找到N个区域，其中满足的区域，及每个区域具有局部最大值；
3. 在每个这样的区域中，通过跟随主应力轨迹绘制潜在的裂纹路径；
4. 通过激活或不激活每个潜在的裂纹路径，对所有可能的组合执行个有限元计算，并找到相应的总势能；
5. 选择能量释放率增量最大的裂纹路径。如果满足能量准则，则通过消除预测裂纹路径中的所有元素来创建相应的新裂纹。转到步骤1；
6. 如果能量准则未被验证，则增加外部负载并转到步骤1。

验证和数值示例

为了验证，提出的数值断裂模型可以通过假定一些特殊的应力状态来完成，在这些特殊的应力状态下，断裂模式在文献中是众所周知的，且被很好的研究。此外，还通过这些示例讨论了模型的效率，特别是网格独立性和多重破解过程建模的能力。

1. 均匀单轴拉伸

第一个实例是经受均匀单轴张力的二维板。Leguillon标准在这种情况下相当于最大应力标准。从理论上讲，当拉伸应力均匀时，断裂可以在板的任何点处起始。在目前的工作中，我们倾向于断裂发生在样品的中间部分。

为了检查数值结果的网格独立性，用四个不同的网格进行计算，网格从细到粗依次变化。数值计算结果表明，所有网格都应严格遵循最大应力断裂准则。图2表示用不同网格获得的裂纹路径。黑线代表预测的裂缝路径。红色的元素代表要删除的元素。对于这种特殊的应力状态，所提出的数值模型与断裂时的临界载荷以及裂纹路径是和网格完全无关的。

(a) d=0.015mm (b) d=0.03mm (c) d=0.0625mm (d) d=0.125mm

图2：不同的尺寸网格划分的板在均匀张力下的断裂模拟

1. 型裂缝

本模型验证中的第二个例子是一块具有型裂缝中心裂纹的板的远端进行加载。在这种情况下，Leguillon准则相当于Griffith准则。

考虑一块二维方板，其中心处有一条长度为2a的裂纹。样品尺寸为2W<!--

剩余内容已隐藏，支付完成后下载完整资料

---

资料编号：[137875]，资料为PDF文档或Word文档，PDF文档可免费转换为Word