An anisotropic damage model combined with a tracking algorithm for modelling crack propagation

Abstract

The development of a simple and efficient methodologies for numerically analyzing the material fracture process is very important in the field of computational mechanics. Damage mechanics approaches are still applied to fracture numerical analyses of many engineering practice problems. This paper focuses on the numerical prediction of crack propagation and fracture behavior by the combination of anisotropic damage model and tracking algorithm. In general, anisotropic damage models may be misunderstood to be used only in the simulations of anisotropic materials. However, it can be used for the anisotropic stiffness matrix induced by the crack plane in damaged isotropic materials. Although it is well known that the anisotropic damage model is superior to the isotropic damage model in fracture simulations, most of studies have combined the isotropic damage model and tracking algorithm, and few studies combine the anisotropic damage model and tracking algorithm. The issues of successfully combining the anisotropic damage model and crack tracking algorithm are addressed in this study. The anisotropic damage model is improved and a local tracking algorithm based on crack surface discretization is also modified. Various crack propagation problems are analyzed numerically to demonstrate the superior performance of the proposed approach.

Keywords

Crack propagation damage mechanics model tracking algorithm crack surface discretization fracture

1. Introduction

Studies on the prediction of discrete crack propagation and failure behavior in a quasi-brittle material have been the focus of computational mechanics over the past three decades. Since theoretical results can be obtained only under special loading and boundary conditions, finite element method (FEM) in conjunction with a reasonable cracking model is widely used to simulate the crack problems. Virtual crack closure technique (VCCT) [31], cohesive zone method (CZM) [27,57], and smeared crack band concept (SCA) [25] based on damage mechanics has been traditionally studied for many years in fields of numerical failure analysis of various types of materials. The damage mechanics model can capture multi-failure modes in the model without the requirement of re-meshing and the presence of initial cracks, and so it is still used to predict fracture behavior of isotropic and anisotropic materials.

The isotropic damage theory was first developed by Kachanov [26] and Rabotnov [44], in which they introduced damage variable that reflected the average of micro cracks distributed within the material. Numerous studies were subsequently carried out to improve the damage model such as effective stress or strain, criteria of damage initiation and evolution, damage tensor of damaged structures, etc. The damage models based on the smeared crack band concept can simulate the propagation of discrete crack within the framework of continuum models. Many anisotropic damage models are also widely investigated to reflect degradation of material properties of damaged elements in the direction perpendicular to the assumed crack plane [29,46]. In order to remove the mesh-induced directional bias of the crack path in the numerical results of finite element analyses, the concept of characteristic length controlled by fracture energy (like CZM) has been introduced [36,59]. When predicting the discrete crack path by the anisotropic damage model, several difficulties such as stress locking or path continuity have been occurred. To solve these issues, studies using the anisotropic–isotropic transient model [25] and element-killing method [22,43,53] have also been conducted.

Many studies have been conducted over the past 20 years on discrete cracking methods to numerically analyze discrete cracks without re-meshing [64,65]. Embedded finite element method (E-FEM) or extended finite element method (X-FEM) have developed with adding discontinuity to elements or nodes [2,5,6,20,21,28,30,34,51]. The tracking algorithm (TA) is necessarily used in analyzing the crack propagation by discrete cracking methods for continuity and robustness of the crack path. Some scholars have pointed out that TA is the most important part of discrete cracking models and that discrete cracking models without TA are not able to well predict the crack path even if the degrees of freedom are added to the elements or nodes [37,55].

Typical TAs are the local tracking algorithm (LTA), scalar field method (also known as global tracking algorithm), level set method (LSM), and explicit tracking algorithm. The LTA is very simple and computationally cheap, but its use is generally limited to 2D simulation due to topological difficulties (topological difficulties are discussed in detail by Jäger et al. [23,24]). The scalar field method developed by Oliver [38,39] requires a high additional computational cost to track the crack path due to the requirement of solution of the thermal-like problem at each loading step, and its results are sensitive to user-defined perturbation parameter [45] and Dirichlet boundary condition [3]. LSM, which uses two orthogonal level set functions to determine the crack surface, also requires the high computational cost to solve the Hamilton-Jacobi type equations [35] and is somewhat cumbersome to implement [19,54,55]. It is not easy to extend LSM from 2D to 3D, and there may be additional inaccuracies when updating the level set function especially in the case of a small increment of crack [16]. In explicit TAs, crack surface mesh is constructed independently of the volume mesh and the prediction of the crack surface does not require additional computational cost. However, updating a dense crack mesh at the crack tip requires a series of techniques that have not yet been given a standard method, making it difficult for general users to use it in commercial FE software [17,42]. Jäger et al. [24] concluded that the current 3D LTAs, such as Belytschko’s local TA [4] and Holzapfel’s two-step TA [18,19], have some drawbacks in stability and robustness and that efficient 3D LTA has not been developed after comparing the performance of several TAs.

Phase field method (PFM) [32,61] is a good alternative method and widely used as the mechanics model to analyze the crack propagation and fracture. It is similar to gradient-enhanced damage mechanics method [10,15,40] and the comparison between PFM and gradient-enhanced damage mechanics model was discussed in detail by De Borst and Verhoosel [12]. Although PFM does not require any TA and naturally describe the crack patterns, it requires the obtaining of auxiliary phase field parameters at all nodes and therefore increases the computational cost. In addition, the regularization parameter is introduced to regularize the discrete crack and stabilize the equilibrium solution. However, this parameter increases the width of crack pattern to 5 ∼ 10 times of the average element size (the other discrete methods generally require crack width of the degree of average element width), thus requiring very fine mesh and increasing the number of iterations for convergence.

Recently, the damage model with TA have attracted much attention in the simulation of fracture problem. Cervera and his team have found that the crack path can be well predicted when a very simple and easy 2D local crack tracking algorithm is added to the isotropic CDM model [9]. The two-dimensional C1 continuous TA, based on the postponement of the crack path fixed moment, was developed and incorporated into the isotropic damage model to reduce the mesh-induced directional bias of the simulated crack path [34,52]. A damage model combined with a 2D LTA capable of simulating the crack initiation at the internal point of the discrete domain and the crack propagation along two opposite orientations, was proposed by Saloustros [49]. In order to solve the sharp turns of the crack path, a combination of a 2D LTA with self-correction ability of the crack path and an advanced damage model has suggested by Yun [60]. Pela investigated an implicit orthotropic damage model based on the mapping relationship and its combination with 2D LTA to predict fracture behavior of common orthotropic material [41]. The Cervera’s LTA [9] was added to the multi-directional crack model based on CDM to analyze the failure process of heterogeneous materials. Since the extension of LTA from 2D to 3D is not trivial, many studies using LTA have limited to 2D. For addressing the topological difficulties, a 3D LTA based on crack surface discretization is proposed by Yun and Wang [61], where crack path is tracked by the connection of crack points located on the sides of solid mesh.

LTA can save much more computational cost compared to other TAs, but most of the studies that combine CDM and TA are limited in 2D, since the extension of TA from 2D to 3D is not trivial and the efficient 3D LTA has not been developed. A desirable 3D LTA is proposed based on crack surface discretization [61], where crack path is tracked by the connection of crack points located on the sides of solid mesh. However, when the crack point is placed near the node, many track lines are needed around this node. That is, the LTA in the study by Yun and Wang [61] was partly independent of the solid mesh. So the improvement of LTA [61] should be performed for forming the fully independent crack surface.

Most of studies have combined the isotropic CDM and LTA. In general, it is well known that anisotropic damage models have a high ability to predict crack patterns compared to isotropic damage models. It is noted that anisotropy in the anisotropic damage model means anisotropy in the damaged stiffness matrix. However, since the traditional anisotropic damage models produce stress locking behavior appeared at the end of damage evolution, the anisotropic damage model is improved to solve this problem.

In this paper, an approach that combines the improved anisotropic damage model and the simple LTA is proposed for failure analysis of quasi-brittle material to increase continuity of crack path and tracking quality. The purpose of the present study is to improve the numerical approach based on damage mechanics in two main aspects: the anisotropic damage model and the crack tracking algorithm. The present damage model is improved to be able to consider the effect of shear stress appeared due to the rotation of stress tensor during damage evolution. By adding the damage variable directly to the stiffness matrix of the second order tensor type, it eliminates an asymmetry of damaged stiffness matrix and has its own ability to track crack path. The anisotropic damage model is constructed by ignoring some Poisson’s ratios during damage evolution to prevent the spurious stress locking. The necessary derivation was performed to guarantee the correct energy dissipation to eliminate the drawback that the dissipated energy is inaccurate when using many anisotropic damage models. In addition, the LTA based on the crack surface discretization is also improved by forming the fully independent crack surface from the mesh. In fact, in the study by Yun and Wang [61] it was partly independent of the solid mesh, so crack path is not completely independent of solid mesh. The mesh dependency of crack path can be avoided by releasing the constraint that crack points should be located on the element faces. By using the maximum crack inclination angle criterion, it prevents the sharp turns of crack path that is out of true. The high accuracy and stability of the current approach in the prediction of crack propagation path are verified by numerical simulation of various benchmark examples.

2. Continuum damage model

An improved anisotropic damage model that can predict the crack path is proposed in this section. Assuming that the energy dissipated due to damage degradation of the finite element is consumed in cracking, anisotropic constitutive equations are constructed to allow degradation of the material properties only in the direction perpendicular to the crack surface. The proposed CDM model not only avoids stress locking, but also be easy to use in common FE software. It also solves the asymmetric problem of the damaged constitutive tensor that would occur when the damage tensor is substituted.

2.1. Constitutive equation

The concept of effective stress is generally used in simulating the fracture process of materials by damage models. The effective stress is expressed as follows. $\begin{eqnarray}\displaystyle \tilde{\boldsymbol{{\sigma}}}=\boldsymbol{C}\boldsymbol{{\varepsilon}},\quad \boldsymbol{{\sigma}}=\boldsymbol{C}_{d}\boldsymbol{{\varepsilon}}, & & \displaystyle\end{eqnarray}$ (1) where $\boldsymbol{{\sigma}},\tilde{\boldsymbol{{\sigma}}}$ and 𝜺 are the Cauchy stress tensor, the effective stress tensor and the strain tensor respectively, C and C _d are the undamaged and damaged elastic stiffness tensors respectively. In isotropic damage models, all components of the stiffness matrix are reduced to the same ratio and expressed by the following equation. $\begin{eqnarray}\displaystyle \boldsymbol{C}_{d}=(1-d)\boldsymbol{C}. & & \displaystyle\end{eqnarray}$ (2)

The damaged elastic matrix of the isotropic CDM has no information related to the crack plane in the model since the constitutive model is constructed assuming that micro cracks are randomly distributed. Cervera [8] pointed out that the energy amount derived from Rankine’s damage model is not the equal to the total energy dissipated in crack formation and is underestimated. In the anisotropic damage model, only some of the stiffness components associated with the crack surface are reduced after damage by accepting the damage tensor. In general, the second order damage tensor is constructed for isotropic damage model so that only stiffness components in the direction perpendicular to the crack plane can be decreased. The damage tensor is then applied to the fourth order flexibility tensor to present material damage. Mathematical notation of these models is simple but lacks physical intuition and has drawbacks of asymmetry of the flexibility matrix. To overcome these drawbacks, studies to modify the damage tensor symmetrically have also been carried out.

In mesoscopic fracture simulation of composite laminates, the method of entering the damage variables directly into the second order compliance tensor is generally applied to distinguish between fiber and the matrix damage. This method is applied to the damage analysis of isotropic materials in the present research due to their more physical institution. It is important that the damaged stiffness matrix of the element can accurately account for the presence of the crack plane and the displacement jumps in anisotropic CDM and that the energy dissipated by stiffness reduction is exactly equal to the energy consumed for the cracking in the element.

The damaged compliance tensor is defined as follows in the local crack coordinates. $\begin{eqnarray}\displaystyle S_{d}^{\prime }=\left(\begin{array}{@{}cccccc@{}}S_{11}/(1-D) & 0 & 0 & 0 & 0 & 0\\ 0 & S_{22} & S_{23} & 0 & 0 & 0\\ 0 & S_{23} & S_{33} & 0 & 0 & 0\\ 0 & 0 & 0 & S_{44}/(1-D) & 0 & 0\\ 0 & 0 & 0 & 0 & S_{55} & 0\\ 0 & 0 & 0 & 0 & 0 & S_{66}/(1-D)\end{array}\right), & & \displaystyle\end{eqnarray}$ (3) where S ₁₁, S ₂₂, S ₃₃, S ₂₃ and S _ii, i = 4,5,6 are the component of compliance matrix of undamaged material respectively, D is the damage variable. The local crack coordinates are described in detail in Section 3. The damaged stiffness tensor in the local crack coordinates is $C_{d}^{\prime }=(S_{d}^{\prime })^{-1}$ .

Unlike previous studies on the anisotropic damage model, the effect of Poisson’s ratios 𝜈₁₂ and 𝜈₁₃ in Eq. (3) is ignored during the damage evolution. Setting Poisson’s ratios 𝜈₁₂ and 𝜈₁₃ as zero will prevent that spurious stress transfers to the other direction due to the large deformation at the maximum principal stress direction. It is based on the hypothesis that stress 𝜎₁₁ only increases the displacement perpendicular to the crack surface for cracking of the element. This hypothesis is consistent with crack propagation principles.

From the compliance matrix in the local coordinate system given in Eq. (3), the stiffness matrix in the local coordinate system is obtained using its inverse and then the transformation matrix is used to obtain C _d. Since it is difficult to insert the damage variable into the stiffness matrix in the global coordinates directly, the expression in the local coordinates is used.

The stiffness matrix of a damaged element in global Cartesian coordinates is expressed as $\begin{eqnarray}\displaystyle \boldsymbol{C}_{d}=\boldsymbol{T}_{{\sigma}}\boldsymbol{C}_{d}^{\prime }\boldsymbol{T}_{{\sigma}}^{T}, & & \displaystyle\end{eqnarray}$ (4) where T _𝜎 is the transformation matrix of stress. The superscript “ ′ ” refers to the local crack coordinates, but the superscript will be omitted in the latter part because only the local coordinate system of the element is considered.

2.2. Damage evolution

The loading function 𝜙 is defined as the normalized equivalent stress. $\begin{eqnarray}\displaystyle {\phi}={\displaystyle \frac{\sqrt{\langle \tilde{{\sigma}}_{11}\rangle ^{2}+{\alpha}\left(\tilde{{\tau}}_{12}^{2}+\tilde{{\tau}}_{13}^{2}\right)}}{{\sigma}_{f}}}, & & \displaystyle\end{eqnarray}$ (5) where 〈⋅〉 is the McCauley operator, $\tilde{{\sigma}}_{1}$ , $\tilde{{\tau}}_{12}$ , $\tilde{{\tau}}_{13}$ are the components of effective stress tensor, 𝛼 is the parameter given in the middle of this subsection, and 𝜎_f is the tension strength under uniaxial stress state.

The effective stress is defined based on the hypothesis of strain equivalence. In order to take into account the effects of tensile and shear stresses on the damage variable together, the tensile and shear components of the effective stress tensor are used when constructing the loading function.

Since the direction perpendicular to the first principal stress is chosen as the x ^′-axis of local coordinates, the shear stresses $\tilde{{\tau}}_{12},\tilde{{\tau}}_{13}$ at the moment of damage initiation are equal to zero and Eq. (5) is simplified in the assessment of damage starting: ${\phi}=\frac{\langle \tilde{{\sigma}}_{11}\rangle }{{\sigma}_{f}}$ .

The damage activation function F is expressed as: $\begin{eqnarray}\displaystyle F={\phi}-r\leq 0, & & \displaystyle\end{eqnarray}$ (6) where r is the elastic domain threshold, which is expressed in time step t + Δt as: $\begin{eqnarray}\displaystyle r_{t+{\rm\Delta}t}=\max _{{\tau}=\overline{0,t}}\{1,\max \{{\phi}_{{\tau}}\}\}, & & \displaystyle\end{eqnarray}$ (7) where 𝜙_𝜏 is the value of the loading function at time 𝜏. It is assumed that the damage variable in Eq. (3) varies according to the exponential softening low. $\begin{eqnarray}\displaystyle D(r)=1-{\displaystyle \frac{1}{r}}\exp \{-2H_{S}(r-1)\}, & & \displaystyle\end{eqnarray}$ (8) where r is the elastic domain threshold, which is the internal variable, and H _S is the softening parameter. The initial value of r for the undamaged material is equal to one and increases with damage evolution.

The total energy per unit volume dissipated during the damage process is expressed as: $\begin{eqnarray}\displaystyle G_{V}=\int _{t=0}^{+\infty }dW. & & \displaystyle\end{eqnarray}$ (9)

The strain energy is $\begin{eqnarray}\displaystyle W={\displaystyle \frac{1}{2}}\boldsymbol{{\sigma}}:{\epsilon}={\displaystyle \frac{1}{2}}\mathop{\sum }_{i,j=1}^{6}{\sigma}_{i}S_{ij}{\sigma}_{j}, & & \displaystyle\end{eqnarray}$ (10) where the stress vector is (𝜎₁, 𝜎₂, 𝜎₃, 𝜎₄, 𝜎₅, 𝜎₆)^T = (𝜎₁₁, 𝜎₂₂, 𝜎₃₃, 𝜏₁₂, 𝜏₂₃, 𝜏₃₁)^T.

Substituting Eq. (10) to Eq. (9) $\begin{eqnarray}\displaystyle G_{V}=\int _{0}^{+\infty }{\dot{W}}∼dt & = & \displaystyle \int _{0}^{+\infty }{\displaystyle \frac{1}{2}}\mathop{\sum }_{i,j=1}^{6}{\sigma}_{i}{\displaystyle \frac{dS_{ij}}{dD}}{\sigma}_{j}{\dot{D}}dt={\displaystyle \frac{1}{2}}\int _{1}^{+\infty }\mathop{\sum }_{i,j=1}^{6}{\sigma}_{i}{\displaystyle \frac{dS_{ij}}{dD}}{\sigma}_{j}D^{\prime }dr\nonumber\\ \displaystyle & = & \displaystyle {\displaystyle \frac{1}{2}}\int _{1}^{+\infty }({\sigma}_{11}^{2}S_{11}+{\tau}_{12}^{2}S_{44}+{\tau}_{13}^{2}S_{66}){\displaystyle \frac{-D^{\prime }}{(1-D)^{2}}}dr.\end{eqnarray}$ (11)

From Eq. (1) and Eq. (3) $\begin{eqnarray}\displaystyle \tilde{\boldsymbol{{\sigma}}}=\boldsymbol{C}\boldsymbol{{\varepsilon}}=\boldsymbol{C}\boldsymbol{S}_{d}\boldsymbol{{\sigma}}. & & \displaystyle\end{eqnarray}$ (12)

From Eq. (12) $\begin{eqnarray}\displaystyle \tilde{{\sigma}}_{11}={\displaystyle \frac{{\sigma}_{11}}{1-D}},\quad \tilde{{\tau}}_{12}={\displaystyle \frac{{\tau}_{12}}{1-D_{s}}},\quad \tilde{{\tau}}_{13}={\displaystyle \frac{{\tau}_{13}}{1-D_{s}}}. & & \displaystyle\end{eqnarray}$ (13)

Substituting Eq. (13) to Eq. (11) $\begin{eqnarray}G_{V}={\displaystyle \frac{S_{11}}{2}}\int _{1}^{+\infty }\left(\tilde{{\sigma}}_{11}^{2}+{\displaystyle \frac{S_{44}}{S_{11}}}\tilde{{\tau}}_{12}^{2}+{\displaystyle \frac{S_{66}}{S_{11}}}\tilde{{\tau}}_{13}^{2}\right)D^{\prime }dr.\end{eqnarray}$

In isotropic material, the components of compliance matrix are as follows: S ₁₁ = 1∕E, S ₄₄ = S ₆₆ = 1∕G = 2(1 + 𝜈)∕E. At this point, the parameter 𝛼 is defined as follows: 𝛼 = 2(1 + 𝜈).

Then the above expression can be rewritten as the following simple form. $\begin{eqnarray}\displaystyle G_{V}={\displaystyle \frac{{\sigma}_{f}^{2}}{2E}}\int _{1}^{+\infty }r^{2}D^{\prime }dr. & & \displaystyle\end{eqnarray}$ (14)

Substituting Eq. (8) to Eq. (14) and integrating, the following equation is obtained. $\begin{eqnarray}\displaystyle G_{V}=\left(1+{\displaystyle \frac{1}{H_{S}}}\right){\displaystyle \frac{{\sigma}_{f}^{2}}{2E}}. & & \displaystyle\end{eqnarray}$ (15)

It is noted that if Poisson’s ratios 𝜈_ij are present in Eq. (3), the derivation maybe more difficult, and it is very difficult to build a model to guarantee the correct energy dissipation, as can be seen from the above derivation. That is, the derivation guarantees the correct energy dissipation to eliminate the drawback that the dissipated energy is inaccurate when using many anisotropic damage models.

The characteristic length l _e of element is introduced to reduce the mesh dependency of the simulation results. Then the relationship between energy dissipated per unit volume and the energy release rate is given as follows [36,48,62]. $\begin{eqnarray}\displaystyle G_{V}={\displaystyle \frac{G_{c}}{l_{e}}}, & & \displaystyle\end{eqnarray}$ (16) where G _c is the energy release rate of material, which is used to measure the material strength in linear elastic fracture mechanics.

From Eq. (15) and Eq. (16), the softening parameter H _S is expressed in terms of the characteristic length of the element. $\begin{eqnarray}\displaystyle H_{S}={\displaystyle \frac{l_{e}}{l_{mat}-l_{e}}}, & & \displaystyle\end{eqnarray}$ (17) where l _mat is the material parameter given as follows: $l_{mat}=∼2EG_{c}/{\sigma}_{f}^{2}$ . If the energy dissipated in the element is close enough to the energy release rate, the element is deleted by the killing element technique and the complete crack path is formed in the element [22,60]. It is noted that the introduction of the characteristic length of the element can significantly reduce the mesh-induced directional bias of the crack path in the FE progressive failure analysis [59].

3. Crack tracking algorithm

In the past, many scholars have found that damage models fail in many fracture problems. However, the combination of damage model and efficient TA brings many improvements in the prediction of crack path and fracture behavior of the material, and can be sufficiently overcome the problems discussed in the past. It can prevent the simultaneous damage of many elements due to large stress at the crack tip when using only SCA, and guarantee that only the elements which will belong to the crack path are damaged. And it can also eliminate mesh dependency of simulated crack path.

In this section, a local 3D crack tracking algorithm proposed by Yun and Wang [61] is updated. The present LTA has low computational cost and high stability and solves the 3D topological difficulties which are the main weakness of 3D LTA. Since the present algorithm is based on crack surface discretization, it is easily extended from 2D local TA to 3D and combined well with the damage mechanics model.

3.1. The basic concept of crack surface discretization

In crack surface discretization, the crack surface is divided into crack triangles and each crack triangle contains three crack points. The crack surface is determined by creating the crack points and forming the crack triangles, and the elements through which the crack surface passes are identified to be damaged. The initial crack points are first created and then the vectors are drawn in the crack propagation direction to generate new crack point successively as in 2D LTA. In the present LTA, the crack lines starting at each initial crack point are ploted from Runge–Kutta principle [60]. A new crack point is created at the end point of crack propagation vector starting at an old crack point. The path of the crack points (crack line) draws the shape of the flow line in the fluid flow as shown in Fig. 1a. The surface that contains all the crack points and crack lines is the crack face we are trying to get. How to draw crack propagation vectors and how to form crack triangles is explained in detail in the following subsection. In the previous LTAs, the intersection between the element faces and the plane perpendicular to the maximum principal stress is assumed to be the crack plane of the element. In the case where the crack planes are entered into an element from several element faces, the multiple facets shared in an element will cause topological problems [3]. So the constraint that one crack triangle is always placed inside an element is eliminated to solve the topological problems of the traditional 3D LTA.

Fig. 1.

Formation of crack surface: (a) crack lines and local coordinates; (b) crack triangles.

3.2. Formation of crack surface

Details on how to create crack propagation direction, local coordinate systems, crack points, and crack triangles are described here.

3.2.1. Creation of initial crack points

It is relatively easy to form the required number of initial crack points at the initial crack tip in the problem given the initial cracks. In the case where no initial crack is given, the crack root element that first meets the damage initiation condition (Eq. (6)) is searched from the geometry boundary. We assume that the crack starts only at the boundary of the solid mesh [9]. The first initial crack point is created at the face center point of the of crack root element located at the solid boundary. Two next initial crack points are generated on the boundary at a certain distance by using the intersection line between the boundary and the plane perpendicular to the maximum principal stress at the first initial crack point. Then new initial crack points are got from these two initial crack points and the initial crack points are created sequentially. Creation of the initial crack points is continued until they reach the other boundary of the solid model. For problems with multiple cracks, a crack exclusion radius is used, which is a radius that allows a new crack root at a certain distance from the current crack tip points (this radius is explained in detail by Yun et al. [60]).

3.2.2. Crack propagation direction and the local coordinate system

The crack propagation direction T at the new crack point is determined along with the following equation when the crack point is created. $\begin{eqnarray}\displaystyle \boldsymbol{S}=\boldsymbol{N}\times \boldsymbol{T}_{0},\quad \boldsymbol{T}=\boldsymbol{S}\times \boldsymbol{N}, & & \displaystyle\end{eqnarray}$ (18) where N is the unit vector of first principal stress in non-local stress calculated around a new crack point. The use of non-local stress can remedy the drawback that the crack plane very much dependent on the stress-state of a single element when using the local stress. T ₀ is the crack propagation vector at the previous crack tip point, and S is the tangent unit vector of the crack surface perpendicular to N and T ₀. The vector T ₀ at the initial crack points is assumed to be a unit vector perpendicular to the geometry boundary. The local coordinate system is generated at the crack point using three vectors N , T , S , and these are used for damaged stiffness matrix in the crack path. The crack propagation direction T is also used to determine the next crack point.

3.2.3. Generation of crack points and crack triangles

Once the element with a crack point satisfies the damage initiation criterion, the generation of a new crack point is performed. The vector is drawn along crack propagation direction T and the tip point is set to the new crack point and the element containing this new crack point is set as the crack tip point. Since the vector is independently of solid mesh, there is no topological problem at this moment. A newly created crack point is defined as belonging to the adjacent element having the common face with the current crack tip element.

Next, the crack triangles are generated to form the crack surface. When one crack point is created, two crack triangles are generated except for the crack points in the first and last crack lines (Fig. 1b). A triangle is drawn consisting of the new crack point, the previous crack point on the same crack line, and the crack tip point on the adjacent crack line. The sum of these crack triangles becomes the crack surface. This method of generating cracked triangles is very simple but very efficient for crack tracking.

3.3. Crack path management

In TA combined with SCA, the elements that intersect the crack surface should be identified and damage evolution should be performed on these elements. Elements with crack points can be included in the crack path because they are tracked along with creating a crack point. The local coordinate system for the anisotropic damage of these elements is chosen as the coordinate system of the crack points on the element. Among the elements cut by the crack triangles, there may be some elements that do not have the crack points. Elements that do not belong to the crack path and intersect the crack triangle are added to the crack path by finding among the elements around the crack triangle. The local coordinate system of these elements is that of the new crack point when the crack triangle is generated.

A set of crack tip elements, which have the elements with the tip point of the crack line, is used to evaluate the possibility of crack propagation at each crack line. If a crack tip element meets the damage initiation condition, a new crack line is drawn within this element and a crack point and two crack triangles are generated. The number of crack tip element is equal to the number of developing crack lines. The set of crack tip elements is updated by adding the element having the new crack tip point to the set and removing the previous tip element from the set.

The crack path correction is applied to the following two cases to improve the reliability of the crack path tracked in the failure analysis [60]. The maximum crack inclined angle criterion is applied to avoid sharp turns of crack propagation direction which is often met in numerical simulations of cracked structure by CDM [9,52], E-FEM [3,45] and X-FEM [13,42]. If the following inequality is not satisfied, the maximum crack inclined angle criterion is performed in the present LTA. $\begin{eqnarray}\displaystyle Arg(\boldsymbol{T},\bar{\boldsymbol{T}})\leq {\alpha}_{max}, & & \displaystyle\end{eqnarray}$ (19) where T is the propagation direction vector at the considered crack point, $\bar{\boldsymbol{T}}$ is the possible propagation vector obtained by the previous three propagation vectors at a considered crack line, and 𝛼_max is the maximum allowable inclined angle defined by the user. Once the angle between T and $\bar{\boldsymbol{T}}$ at the current crack point exceeds the maximum allowable inclined angle 𝛼_max, T is corrected to $\bar{\boldsymbol{T}}$ .

Elements with a damage variable D greater than 0.999 are considered to have fully dissipated the crack energy and are deleted by the killing element technique. The complete crack formation is accomplished by deleting the damaged element, while the damaged living elements perform energy release only in the direction perpendicular to the selected crack plane of the element.

4. Results and discussion

In this section, the effectiveness and robustness of our approach is verified through solving three 3D benchmark problems.

4.1. The planar crack propagation under bending

The prediction of crack path of a concrete beam subject to three-point bending is considered here. An initial crack of deep 50 mm and wide 5 mm is located in the centerline of the bottom plane of the concrete beam [9,47,63]. The geometry and load for three-point bending test are depicted in Fig. 2. The material properties are given as follows: E = 20 kN/mm², 𝜈 = 0.2, 𝜎_f = 2.4 N/mm², G _c = 0.113 N/mm. Calculation is performed with an imposed maximum vertical displacement of 1 mm which is incrementally applied to the beam with an incremental displacement Δu = 0.01 mm.

Fig. 2.

Three-point bending test.

Fig. 3.

Deformed shapes in different meshes (×100): (a) Mesh 1; (b) Mesh 2.

In this example, the capability of our approach for tracking the planar crack path is tested under three-point bending. In all of the examples in this study, the parameter 𝛼_max for applying the maximum crack inclined angle criteria is set to 30°.

Fig. 4.

Contour plots of damage variable (×1): (a) Mesh 1; (b) Mesh 2.

Fig. 5.

Load-displacement curves for three-point bending beam.

Fig. 6.

Three-point bending specimen with the curved crack path.

Fig. 7.

Results by isotropic damage model: (a) deformed shape (×5); (b) contour of damage variable (×1).

Fig. 8.

Results by our anisotropic damage model without TA: (a) deformed shape (×5); (b) contour of damage variable (×1).

Two tetrahedron meshes are employed to consider mesh-induced directional bias of the simulated crack path. The simulated deformed shapes in two different meshes are shown in Fig. 3. Although different meshes are used for the simulations, the predicted cracks in both cases rise from the tips of the initial crack to the loading location in a planar fashion. The contour plots of damage variable (at the final time step) are shown in Fig. 4. In the figure, the red parts of the crack path represent almost fully damaged elements of which damage variable is approximately one. The deformed shapes and contours of damage variable calculated in the case of two different meshes show that there is very little mesh induced directional bias by our approach.

A comparison of the numerical load-displacement curves and experiment [47] is shown in Fig. 5. It can be easily found from the figure that two predicted load-displacement curves are very similar and agree well with the experiment. It is clear that the present approach has a high capability of predicting the planar crack propagation and load bearing characteristics in the fracture analysis of quasi-brittle material.

4.2. Mixed mode failure

The mixed mode failure problem of a beam subject to three-point bending is solved to test the tracking ability of two damage models themselves in this example [60]. The geometry and the load conditions are presented in Fig. 6. The initial crack is located at the bottom plane of the beam offset by 0.8 mm from the centerline of the bottom plane and an imposed vertical displacement is applied at the centerline of the top plane. The following material data is used in numerical analysis: E = 0.1 kN/mm², 𝜈 = 0, f _t = 0.5 N/mm², G _f = 0.01 N/mm.

This example explains the difference in the prediction of crack propagation path of isotropic damage model and anisotropic damage model. The deformed shapes and the contours of damage variable simulated by isotropic damage model and anisotropic damage model are shown in Figs 7–8. It can be seen that isotropic damage model does not predict the curved crack path at all but the anisotropic damage model predicts the curved crack path appropriately in this example. It indicates that the anisotropic damage model has far better capability than the isotropic damage model in crack propagation prediction. It can be seen from this example that the use of the anisotropic damage model is far superior to the failure problem analysis compared to the isotropic damage model.

However, it should be noted that the simulation results by traditional anisotropic damage models without the TA depend on the mesh quality as mentioned above. Figure 9 shows the deformed shape and contour of damage variable when adding our TA. The prediction ability of the crack path becomes very high when the TA is added with anisotropic damage model.

Fig. 9.

Results by our anisotropic damage model with LTA: (a) deformed shape (×5); (b) contour of damage variable (×1).

Fig. 10.

Comparison of load-displacement curves.

The load-displacement curves by the three methods used above are compared in Fig. 10. The magnitudes of the maximum loads by isotropic and anisotropic damage models, are less predicted than the ones by the present approach with LTA and X-FEM due to the difference in the crack path. The load-displacement curves by our approach and one by X-FEM are very similar, although the last parts have slight differences in the unloading state.

4.3. Four-point shear test

In this example, the single edge notched (SEN) beam specimen subjected to shear loading, which is experimented by Schlangen [50] and has been investigated by several numerical simulations [1,7,11,52,58], is examined to test whether our approach can predict sliding crack surface well.

Fig. 11.

Four-point shear test: geometries and loading conditions.

Fig. 12.

Contours of damage variable in different element sizes (×1): (a) 3990 elements; (b) 15470 elements.

Fig. 13.

Deformed shapes in different element sizes (×30): (a) 3990 elements; (b) 15470 elements.

Fig. 14.

Crack surfaces tracked by numerical analyses.

The SEN beam which has a notch with 5 mm wide and 20 mm deep at the top is enforced by the monotonic increasing crack mouth sliding displacement (CMSD). The experimental setup of specimen is shown in Fig. 11. The material properties are as follows: E = 30 kN/mm², 𝜈 = 0.2, 𝜎_t = 2.4 N/mm², G _f = 0.08 N/mm.

Two different FE meshes with 3990 elements (element size 5 mm) and 15470 elements (element size 3 mm) are used in the calculation to account for the effect of mesh size on the crack path. Comparisons of numerical results for damage contours, deformed shapes and crack paths in different meshes are shown in Fig. 12–Fig. 14 respectively. 7 and 12 initial crack points are used respectively for tracking crack surfaces in two meshes. In two meshes, the smooth curvilinear crack from the notch tip to the bottom side are all well predicted and similar in simulations by the present approach. The curvilinear crack surfaces are all similar in two structured meshes of different sizes, indicating that there is little mesh quality sensitivity in simulation results by the present model.

A comparison between the crack surface by experiments [50] and the crack surfaces simulated by our approach is shown in Fig. 15. As shown, the numerically simulated crack surfaces in two cases are also well coincided with crack surfaces obtained by Schlangen’s experiment. With this example, the accuracy and mesh independence of the simulation results by our approach are fully demonstrated.

Fig. 15.

Comparison of numerical crack surfaces with experiments for SEN beam.

5. Conclusions

In the paper, the anisotropic damage model combined with LTA is proposed. Although damage models regularized by a non-local stress measures are known to fail in many benchmark examples in the past, we have developed the approach addressing these problems by the combination of the improved anisotropic damage model with TA. Stability and efficiency of our damage model in crack path tracking are demonstrated by solving the curved crack propagation problem. The 3D extension of the method combining the damage model and LTA has been made possible by the 3D LTA developed in this paper. The fact that the present approach is not only highly capable of fracture prediction but also is not affected by the mesh quality, is proved by evaluating the effects of mesh sizes and mesh biases on crack path.

Footnotes

Conflict of interest

None to report.

References

Alfaiate

Wells

and Sluys

, On the use of embedded discontinuity elements with crack path continuity for mode-I and mixed-mode fracture, Engineering Fracture Mechanics 69 (2002), 661–686.

Alsahly

Callari

and Meschke

, An algorithm based on incompatible modes for the global tracking of strong discontinuities in shear localization analyses, Computer Methods in Applied Mechanics and Engineering 330 (2018), 33–63.

Annavarapu

Settgast

R.R.

Vitali

, A local crack-tracking strategy to model three-dimensional crack propagation with embedded methods, Computer Methods in Applied Mechanics and Engineering 311 (2016), 815–837.

Areias

P.M.A.

and Belytschko

, Analysis of three-dimensional crack initiation and propagation using the extended finite element method, International Journal for Numerical Methods in Engineering 63 (2005), 760–788.

Azdoud

Han

and Lubineau

, The morphing method as a flexible tool for adaptive blocal/non-local simulation of static fracture, Computational Mechanics 54 (2014), 711–722.

Benvenuti

, An effective XFEM with equivalent eigenstrain for stress intensity factors of homogeneous plates, Computer Methods in Applied Mechanics and Engineering 321 (2017), 427–454.

Cervenka

and Cervenka

, On the Uniqueness of Numerical Solutions of Shear Failure of Deep Concrete Beams: Comparison of Smeared and Discrete Crack Approaches. Computational Modeling of Structures. EURO-C 2010, Taylor & Francis Group, 2010.

Cervera

and Chiumenti

, Mesh objective tensile cracking via a local continuum damage model and a crack tracking technique, Computer Methods in Applied Mechanics and Engineering 196 (2006), 304–320.

Cervera

Pela

Clemente

, A crack-tracking technique for localized damage in quasi-brittle materials, Engineering Fracture Mechanics 77 (2010), 2431–2450.

10.

De Borst

Benallal

and Heeres

O.M.

, A gradient-enhanced damage approach to fracture, Journal of Physics IV C6 (1996), 491–502.

11.

De Borst

and Pamin

, Gradient plasticity in numerical simulation of concrete cracking, European Journal of Mechanics A/Solids 15 (1996), 295–320.

12.

De Borst

and Verhoosel

C.V.

, Gradient damage vs phase-field approaches for fracture: similarities and differences, Computer Methods in Applied Mechanics and Engineering 312 (2016), 78–94.

13.

Dumstorff

and Meschke

, Crack propagation criteria in the framework of X-FEM-based structural analyses, International Journal for Numerical and Analytical Methods in Geomechanics 31 (2010), 239–259.

14.

Feist

and Hofstetter

, Three-dimensional fracture simulations based on the SDA, International Journal for Numerical and Analytical Methods in Geomechanics 31 (2010), 189–212.

15.

Frémond

and Nedjar

, Damage, gradient of damage and principle of virtual power, International Journal of Solids & Structures 33 (1996), 1083–1103.

16.

Fries

and Baydoun

, Crack propagation with the extended finite element method and a hybrid explicit–implicit crack description, International Journal for Numerical Methods in Engineering 89 (2012), 1527–1558.

17.

Garzon

O’Hara

Duarte

C.A.

, Improvements of explicit crack surface representation and update within the generalized finite element method with application to three-dimensional crack coalescence, International Journal for Numerical Methods in Engineering 97 (2014), 231–273.

18.

Gasser

T.C.

and Holzapfel

G.A.

, Modeling 3D crack propagation in unreinforced concrete using PUFEM, Computer Methods in Applied Mechanics and Engineering 194 (2005), 2859–2896.

19.

Gasser

T.C.

and Holzapfel

G.A.

, 3D crack propagation in unreinforced concrete: A two-step algorithm for tracking 3D crack paths, Computer Methods in Applied Mechanics and Engineering 195 (2006), 5198–5219.

20.

Han

Lubineau

Azdoud

, A morphing approach to couple state-based peridynamics with classical continuum mechanics, Computer Methods in Applied Mechanics and Engineering 301 (2016), 336–358.

21.

Han

Lubineau

and Azdoud

, Adaptive coupling between damage mechanics and peridynamics: A route for objective simulation of material degradation up to complete failure, Journal of the Mechanics and Physics of Solids 94 (2016), 453–472.

22.

Ismail

Zaïri

Nait-Abdelaziz

, How cracks affect the contact characteristics during impact of solid particles on glass surfaces: A computational study using anisotropic continuum damage mechanics, International Journal of Impact Engineering 40 (2012), 10–15.

23.

Jäger

Steinmann

and Kuhl

, Towards the treatment of boundary conditions for global crack path tracking in three-dimensional brittle fracture, Computational Mechanics 45 (2009), 91–107.

24.

Jäger

Steinmann

and Kuhl

, Modeling three-dimensional crack propagation. A comparison of crack path tracking strategies, International Journal for Numerical Methods in Engineering 76 (2010), 1328–1352.

25.

Jirasek

and Zimmermann

, Rotating crack model with transition to scalar damage, Journal of Engineering Mechanics 124 (1998), 277–284.

26.

Kachanov

L.M.

, On the time to failure under creep condition, Izv Akad Nauk USSR, Otd Techn Nauk 8 (1958), 26–31, [In Russian].

27.

Khisamitov

and Meschke

, Variational approach to interface element modeling of brittle fracture propagation, Computer Methods in Applied Mechanics and Engineering 328 (2018), 452–476.

28.

Kim

and Armero

, Three-dimensional finite elements with embedded strong discontinuities for the analysis of solids at failure in the finite deformation range, Computer Methods in Applied Mechanics and Engineering 317 (2017), 890–926.

29.

Kishta

Giry

Richard

, A discrete anisotropic damage constitutive law with an enhanced mixed-mode kinematics: Application to RC shear walls, Engineering Fracture Mechanics 184 (2017), 121–140.

30.

Kishta

Richard

Giry

, Strong discontinuity analysis of a class of anisotropic continuum damage constitutive models. Part II: Concrete material application, Mechanics Research Communications 86 (2017), 27–31.

31.

Krueger

, Virtual crack closure technique: History, approach, and applications, Applied Mechanics Reviews 57 (2004), 109–143.

32.

Miehe

Hofacker

and Welschinger

, A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits, Computer Methods in Applied Mechanics and Engineering 199 (2010), 2765–2778.

33.

Miehe

Welschinger

and Hofacker

, Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations, International Journal for Numerical Methods in Engineering 83 (2010), 1273–1311.

34.

Mosler

and Meschke

, Embedded crack vs. smeared crack models: A comparison of elementwise discontinuous crack path approaches with emphasis on mesh bias, Computer Methods in Applied Mechanics and Engineering 193 (2004), 3351–3375.

35.

Nagashima

and Sawada

, Development of a damage propagation analysis system based on level set XFEM using the cohesive zone model, Computers and Structures 174 (2016), 42–53.

36.

Oliver

, A consistent characteristic length for smeared cracking models, International Journal for Numerical Methods in Engineering 28 (1989), 461–474.

37.

Oliver

Dias

I.F.

and Huespe

A.E.

, Crack-path field and strain-injection techniques in computational modeling of propagating material failure, Computer Methods in Applied Mechanics and Engineering 274 (2014), 289–348.

38.

Oliver

Huespe

A.E.

Samaniego

, On strategies for tracking strong discontinuities in computational failure mechanics, 2002.

39.

Oliver

Huespe

A.E.

Samaniego

, Continuum approach to the numerical simulation of material failure in concrete, International Journal for Numerical and Analytical Methods in Geomechanics 28 (2004), 609–632.

40.

Peerlings

R.H.J.

de Borst

Brekelmans

W.A.M.

and de Vree

J.H.P.

, Gradient enhanced damage for quasi-brittle materials, International Journal for Numerical Methods in Engineering 39 (1996), 3391–3403.

41.

Pelà

Cervera

Oller

, A localized mapped damage model for orthotropic materials, Engineering Fracture Mechanics 124 (2014), 196–216.

42.

Pereira

J.P.

Duarte

C.A.

and Jiao

, Three-dimensional crack growth with hp-generalized finite element and face offsetting methods, Computational Mechanics 46 (2010), 431–453.

43.

Petracca

Pelà

Rossi

, Micro-scale continuous and discrete numerical models for nonlinear analysis of masonry shear walls, Construction and Building Materials 149 (2017), 296–314.

44.

Rabotnov

Y.N.

, On the equation of state of creep, Proceedings of the Institution of Mechanical Engineers Conference 178 (1963), 117–122.

45.

Riccardi

Kishta

and Richard

, A step-by-step global crack-tracking approach in EFEM simulations of quasi-brittle materials, Engineering Fracture Mechanics 170 (2016), 44–58.

46.

Richard

Kishta

Giry

, Strong discontinuity analysis of a class of anisotropic continuum damage constitutive models. Part I: Theoretical considerations, Mechanics Research Communications 86 (2017), 32–36.

47.

Rots

J.G.

, Computational modeling of concrete fracture, PhD thesis, Delft University of Technology, 1988.

48.

Rots

J.G.

and Blaauwendraad

, Crack models for concrete, discrete or smeared? Fixed, multi-directional or rotating? HERON 34 (1989).

49.

Saloustros

Cervera

and Roca

, Finite element modelling of internal and multiple localized cracks, Computational Mechanics 59 (2016), 299–316.

50.

Schlangen

, Experimental and numerical analysis of fracture processes in concrete, PhD Thesis, Delft University of Technology, Delft, 1993.

51.

Slobbe

A.T.

, Propagation and band width of smeared cracks, 2015.

52.

Slobbe

A.T.

Hendriks

M.A.N.

and Rots

J.G.

, Smoothing the propagation of smeared cracks, Engineering Fracture Mechanics 132 (2014), 147–168.

53.

Sun

and Khaleel

M.A.

, Modeling of glass fracture damage using continuum damage mechanics-static spherical indentation, International Journal of Damage Mechanics 13 (2004), 263–285.

54.

Ventura

Budyn

and Belytschko

, Vector level sets for description of propagating cracks in finite elements, International Journal for Numerical Methods in Engineering 58 (2003), 1571–1592.

55.

Wang

Bui

T.Q.

, 3-D local mesh refinement XFEM with variable-node hexahedron elements for extraction of stress intensity factors of straight and curved planar cracks, Computer Methods in Applied Mechanics and Engineering 313 (2017), 375–405.

56.

Wells

G.N.

, Discontinuous Modelling of Strain Localisation and Failure, University of Cambridge, 2001.

57.

Wells

G.N.

and Sluys

L.J.

, A new method for modelling cohesive cracks using finite elements, International Journal for Numerical Methods in Engineering 50 (2001), 2667–2682.

58.

Wells

G.N.

and Sluys

L.J.

, Three-dimensional embedded discontinuity model for brittle fracture, International Journal of Solids & Structures 38 (2001), 897–913.

59.

Xin

S.H.

and Wen

H.M.

, A progressive damage model for fiber reinforced plastic composites subjected to impact loading, International Journal of Impact Engineering 75 (2015), 40–52.

60.

Yun

Kim

T.J.

Jang

P.S.

, An improved crack tracking algorithm with selfcorrection ability of crack path and its application in continuum damage model, International Journal for Numerical Methods in Engineering (2018).

61.

Yun

and Wang

, A computational methodology for simulating quasi-brittle fracture problems, Computers and Structures 215 (2019), 65–79.

62.

Yun

Wang

, A damage model based on the introduction of a crack direction parameter for FRP composites under quasi-static load, Composite Structures 184 (2018), 388–399.

63.

Yun

Wang

Ronald

, An advanced continuum damage mechanics model for predicting the crack progress process based on the consideration of the influence of crack direction under quasi-static load, International Journal of Mechanical Sciences 130 (2017), 487–496.

64.

Zhang

Dissanayake

Saputra

, Three-dimensional damage analysis by the scaled boundary finite element method, Computers and Structures 206 (2018), 1–17.

65.

Zhang

Z.H.

Yang

Z.J.

Liu

G.H.

, An adaptive scaled boundary finite element method by subdividing subdomains for elastodynamic problems, Science China Technological Sciences 54 (2011), 101–110.