A step-by-step global crack-tracking approach in E-FEM simulations of quasi-brittle materials

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Introduction
In the last decade, the Embedded Finite Element Method (E-FEM) has gained wide popularity for the description of cracking phenomena [1][2][3][4][5][6]. Due to the local nature of the kinematic enrichment, this approach presents some advantages concerning the computational effort with respect, for example, to the Extended Finite Element Method (X-FEM) [7][8][9][10][11]. Indeed the additional degrees of freedom can be statically condensed, since they are expressed at the element level in terms of crack opening displacement components [5], thus leaving unchanged the dimension of the global system of equations. As a counterpart, since no information is available in the vicinity of cracked elements, the continuity of the crack path is not intrinsically guaranteed and a loss in objectivity may be encountered in numerical simulations [12,13]. For this reason, the E-FEM is usually associated with the adoption of tracking algorithms leading to C 0 -continuous crack paths [14][15][16][17]. Different crack-tracking strategies can be found in the literature. In general, two families may be distinguished: local tracking and global tracking algorithms [14,5,18,19].
Local tracking relies on geometry-based or energy-based schemes, which are applied for each element able to crack. The main difference between the two strategies stands in the evaluation of the direction of the propagating discontinuity. In case ⇑ Corresponding author.
E-mail address: benjamin.richard@cea.fr (B. Richard). of a geometrical approach, the crack orientation is given by the assumed failure criterion, e.g. accordingly to Rankine it is supposed to be perpendicular to the maximum tensile principal direction [20,5,21], whereas in case of an energetic approach, it is computed from linear-elastic fracture mechanics (LEFM) by minimization of the mechanical energy [18,22]. From the knowledge of the root of each discontinuity, i.e. the material point experiencing failure and not associated to any pre-existing crack path, it is possible to make the input point of the crack inside an element match the output point of the crack present in an adjacent element. Local tracking techniques are able to reproduce continuous crack paths in a robust manner by exploiting the informations of nearby elements. Nevertheless, their implementation in the case of multiple crack problems may be cumbersome and this strategy can loose much of its robustness.
Global tracking has been introduced to overcome the limits of local strategies when dealing with multi-cracking modeling [14,5]. The main idea is to trace the envelopes of the tangent vector field to the discontinuities as the isovalues of the temperature field of a heat conduction-like problem in the case of steady state conditions and no internal sources. Dirichlet and Neumann boundary conditions must be prescribed on the respective portions of the boundary, in terms of fixed temperature values and (null) heat flux. Once the assumed failure criterion is fulfilled for the first time in a certain material point, the latter becomes the root of a new discontinuity line, which can be traced as the isovalue passing through that point. Numerical applications can be found in [23,18,24]. The main advantage of this approach is the fact that no information from the neighbourhood of the cracked elements is required to perform the analysis: indeed, since the isovalues are available at every point of the domain, only root element coordinates shall be provided in order to ensure a continuous crack path.
The finite element formulation of the heat conduction-like problem is straightforward. However, the stiffness-like matrix deduced from the anisotropic conductivity tensor reveals to be singular and a user-defined perturbation (isotropic algorithmic conductivity) is introduced in numerical simulations to circumvent this drawback. The dependence of the solution on this parameter may then represents a limitation for the application of global tracking, since its value changes for each specific structural problem. In addition, it is found out that the capability of the thermal-like isovalues to envelop the vector tangent field is reduced as soon as cracking occurs. This fact, due to the rotation of the principal stress directions outside the region crossed by the discontinuity, may lead to a loss of continuity of the crack path and, in the worst case, to a wrong evaluation of the enriched shape functions whenever the element domain is not decomposed properly [12,25]. Such a circumstance becomes increasingly critical as the evolution of the principal stress field is important, see Fig. 1.
The objective of this paper is to provide an alternative formulation of the global crack-tracking strategy able to improve the performance of this technique with respect to the aforementioned issues. This paper is organized as follows. In Section 2 a new physical interpretation of the problem is given in terms of Navier-Stokes equations, where the concept of numerical diffusion is introduced in order provide a stable and consistent solution of the initially ill-posed discrete problem. A revised algorithm for ensuring continuous crack paths in case of step-by-step analysis is presented in Section 3, considering the evolution of the root for the choice of the isovalue enveloping the propagating discontinuity. Section  application of the proposed model to the E-FEM by means of benchmark tests. Section 5 then concludes with a critical comparison between the proposed approach and the original formulation.

Global equations
The problem of tracing the envelopes of a vector field TðxÞ in a domain X has been treated in [14,5]. For the sake of simplicity, let us focus upon the bi-dimensional case. If we indicate with NðxÞ the normal vector field to TðxÞ, we can consider a scalar function hðxÞ whose gradient is parallel to NðxÞ, i.e. such that: Hence, the following partial differential equation will hold in the domain X: Since the level contours of the function hðxÞ are orthogonal to the gradient, the envelope of the vector field TðxÞ passing through a generic point P can be defined as: The envelopes of the vector field T thus provide C 0 -continuous curves, which are well-suited to model crack-paths within the framework of the E-FEM. From now on the dependence of all the quantities on x will be omitted.

Heat conduction-like problem
Eq. (2) can be manipulated by multiplying it by the vector field T. After some analytical computations and using the same notations as in [14], the previous problem can be reformulated as the following boundary value problem for the unknown function h: with: where denotes the tensor product. The boundary value problem (4) defines a heat conduction-like problem in the domain X, where h is the temperature field and where q is the conduction flux vector. Boundary conditions are prescribed on the boundary @X ¼ @ q X [ @ h X such that @ q X \ @ h X ¼ £. More precisely, Eq. (4c) expresses the Neumann condition of a null heat flux on the set @ q X, whereas Eq. (4d) represents the Dirichlet condition of fixed temperature values on the boundary @ h X.
From expression (5), it turns out that the conductivity tensor K is singular. In order to avoid the ill-posedness of the conduction problem, the following isotropic perturbation is introduced [14,5]: where T x and T y are the Cartesian components of the vector field T and where is a user-defined numerical parameter. No rigorous criterion is available for its choice. However, it should be as small as possible in order to fulfil Eq. (4a), but sufficiently large to break down the singularity of K. The dependence of the results on this numerical parameter may then limit the applicability of global tracking, eventually leading to numerical instability issues.

Heat convective-diffusion-like problem
With the aim of overcoming the aforementioned limitation, a new interpretation of the original problem defined by Eq. (2) is here proposed. To start with, let us consider a convection-diffusion-like problem, which consists in finding a temperature field h such that: where a is a diffusion coefficient and m the unit outward normal vector to the boundary. Problem (7) can be derived from Navier-Stokes equations in case of incompressible fluids. In this context, the vector field T takes the physical meaning of a fluid velocity. Therefore Eq. (7a) describes a heat transfer process where two contributions can be distinguished: the first one of convective nature, the second one of diffusive nature. Boundary condition (7b) expresses the Neumann condition on the diffusive flux J, while Eq. (7c) assumes the same meaning as in problem (4). However, since T is now a velocity field, Dirichlet boundary conditions should be prescribed only on the portions of the boundary where the former is directed inward the domain X.
In order to better understand the preceding formulation, let us focus upon the one-dimensional case. Assuming a constant diffusion coefficient and considering, on the one hand, a first order upwind scheme for the convective term and, on the other hand, a second order centered scheme for the diffusive term, Eq. (7a) is discretized as: where Dx is the one-dimensional spatial discretization step. If we consider now a centered discretization of the convective term, it turns out that the difference between the upwind discretization and the aforementioned one is given by: Eq. (9) shows that the discretized expression of the convective term by means of the first order upwind scheme is equal to the second order centered scheme discretization of the same term plus an additional diffusive contribution. By comparing Eqs. (8) and (9), we notice that a numerical diffusive term comes out naturally and it is characterized by a diffusion coefficient TDx 2 which is function of a mesh characteristic length -Dx in the case of one-dimensional problems. This term is similar to the numerical conductivity coefficient introduced in Eq. (6) but it is no more user-defined and it tends towards zero as Dx ! 0, which means that the centered and the upwind discretization schemes are consistent. This observation constitutes the fundamentals of upwind discretization methods, classically used in fluid mechanics [26]. Higher dimensional extensions of the concept of numerical diffusion have been well-established in fluid mechanics in the case of finite element discretizations. More precisely, two formulations are considered in this study. Given the vector field TðxÞ, the following bi-dimensional extension of the concept of numerical diffusion has been proposed [26] (Streamline Upwind method -SU): where h T is a mesh characteristic length. Problem (10) describes the numerical diffusion as an isotropic process and this may lead to imprecise results if coarse meshes are adopted. As a matter of fact, since in this case a large amount of diffusion is introduced in all directions, the isovalues of the temperature-like field do not longer envelop the tangent vector field. Another formulation, which considers an anisotropic numerical diffusion only in the direction of the vector field T (Streamline Upwind Petrov Galerkin method -SUPG) has also been proposed: It can be demonstrated that the tensor TT kTk 2 has only one non-zero eigenvalue associated to T as eigenvector. The product TT kTk 2 rh then retains only the part of rh parallel to T.

Preliminary results
Concrete structures often experience mixed mode fracture mechanisms, due to the strong interaction between shear and bending. This occurrence can translate into curved crack paths developing from the zones where maximum tensile principal stresses concentrate. Therefore, global tracking must allow to reproduce such behavior. In particular, the capability of the three formulations presented in Sections 2.1 and 2.2 to envelop the principal stress field is now investigated. These are recalled in Table 1.
The well-known experiment performed by Schlangen [27] has been simulated. The specimen is a single-edge notched (SEN) plain concrete beam subjected to four-point shear (see Fig. 2).
The isovalues of the thermal-like field obtained by solving problems (4), (10) and (11)  In case of the HC formulation, a numerical perturbation ¼ 10 À5 has been considered (see [14,28] for typical values of this parameter). As it can be observed in Fig. 3b, a loss of uniqueness of the solution is encountered near the notch for the intermediate discretization.
The isovalues obtained by the HCD-SU formulation, implemented in the finite element software Cast3M-CEA [29], are depicted in Fig. 4. It can be noticed that the solution is now more diffused, which means that the prescribed temperature values are more easily transported all over the domain. This feature gives stability to the solution, but at the same time it provides less accurate results if coarse meshes are adopted. However, for finer discretizations, the results are closed to those obtained for the HC problem.
The HCD-SUPG formulation seems to be a good candidate to reduce the dependence on the mesh refinement. The results are shown in Fig. 5. In fact, the solution reveals to be unstable, even though less diffused than in the previous case. Therefore, the introduction of a numerical diffusion only in the direction of the tangent field T is not sufficient to solve out our problem. As a consequence of the above considerations, the most appropriate convection-diffusion-like formulation of the global tracking problem is the one given by Eqs. (10).

Crack-tracking algorithm
In Section 2.2 the problem of tracing the envelopes of a vector field Tðx; tÞ has been formulated in terms of Navier-Stokes equations for incompressible fluids. As for the heat-conduction-like formulation, the scalar function hðx; tÞ represents the temperature field whose isovalues describe all the possible discontinuity lines in the domain X.
The choice of the right isovalue stands on the stress distribution at time t. In particular, global tracking associates to each discontinuity line C i a root r i , i.e. the material point (or the element) at time t 0 not belonging to any crack path and satisfying for the first time the activation condition [14,5]. The discontinuity is thus represented as follows: where the activation condition is expressed in terms of a certain tensorial norm : k k and the material strength f t . The reference isovalue is considered to pass through the centroid of the root element of the discontinuity [14,5]. Table 1 Possible strategies to envelop the principal stress field.

HC
Heat conduction HCD-SU Heat convection-diffusion SU HCD-SUPG Heat convection-diffusion SUPG The previous definition implicitly assumes that the portion of isovalue associated to the active part of the crack, i.e. with points characterized by sut C i -0, does not change any more. In reality, due to the approximative nature of the finite element solution, this is not generally guaranteed: indeed, since principal stress directions are free to rotate where the material is linear elastic, the isovalues of the thermal field may evolve also inside the elements already exhibiting a crack. As a consequence, the nodes of the element domain X e traversed by the crack may be incorrectly separated, with a loss of continuity of the crack-path and stress-locking effects taking place. A possible solution would be to freeze the nodal temperature of the cracked element by imposing additional Dirichlet boundary conditions [6,28]. However, the initial boundary value problem (4) or (10) would be contradicted with respect to the initial boundary conditions and local techniques should be adopted in order to perform the analysis at each time step [25].   This drawback can be found if the root r i is not updated for t P t 0 . Thus, the tracking procedure should explicitly take into account both the active part C ia and the potential part C ip of each discontinuity C i and, at the same time, impose the match between the root r i , used to trace the prosecution of the crack, and the crack-tip. In these circumstances, the reference point for tracing the isovalue defining C ip does not coincide with the root element centroid but instead with a point belonging to one of the root element edges, i.e. the crack-tip. If we consider a time discretization ft 0 ; . . . ; t k ; . . . ; t n g, the total discontinuity C i at time t n can then be represented as: x ðkÞ r i 2 X j rðx ðkÞ From the previous definition, the potential part of the discontinuity C ip has been defined as the portion of the isovalue triggered off the root r i at time t n and characterized by linear elastic behavior. Consequently, only the portion of the isovalue that does not cross the element associated to r i should be taken into account. Since root r i divides the potential line into two parts, in order to avoid ambiguity, it seems convenient to orient the curve by setting the origin at the root itself and assume as positive the direction of the propagating discontinuity. This can be done by means of a curvilinear abscissa s i , whose origin is set to coincide with root r i (see Fig. 6).
Thus, the prosecution of the crack path C i at time t n will be associated to the positive values of s i , with origin at the root r i .
This procedure can be translated by the following steps: 1. Make the input point I C i match the crack-tip. 2. Trace of the isovalue passing trough I C i .   4. Choose the potential continuation of discontinuity C i as the part of isovalue associated to s i > 0.
The algorithms derived by problem (12) and problem (13) are now compared. The test reported in [30] has been simulated. The geometry is depicted in Fig. 7 and consists in a concrete specimen under tension, characterized by two notches with an offset in the direction of the load. An imposed displacement d is applied horizontally on the right side of the structure, while the left side is fixed in this direction. A hinged support is introduced at the upper left corner in order to forbid any rigid body motion.
In Fig. 8a the case of fixed roots r 1 and r 2 for the crack paths C 1 and C 2 is shown. It appears that the evolution of the reference isovalues due to the rotation of the principal stress directions inside linear elastic elements may lead to an imprecise evaluation of the propagating discontinuities. In particular, even if the principal stress directions are freezed inside cracked elements, the isovalues may no longer be able to envelop the tangents to the active part of the discontinuities. A change in the decomposition mode of the elemental domain X e may then occur whenever the nodes shared by adjacent elements do not belong to the same sub-domain X þ e or X À e . If the root of each discontinuity is updated, it is possible to separate its active part from its potential prosecution. This strategy allows to enforce the continuity of the crack path even if the isovalue distribution evolves during the analysis. As depicted in Fig. 8b, the requirement of a domain of unique decomposition is attained for all the cracked elements.

Numerical examples
In this section, the performances of the global crack-tracking technique applied to the E-FEM are investigated. A comparison between the version presented in [14,5] and the version proposed in this paper is made by means of bi-dimensional benchmark tests, presented in order of increasing complexity. Both the mathematical formulation and the tracking scheme are here discussed.

Double-edge notched specimen under tension
The first example is the same as the double-edge notched specimen under tension that has been preliminarily studied in Section 3. The formation of two principal crack paths has been observed experimentally in [30]. A correct evaluation of the propagating discontinuities is essential in order to properly simulate the structural response. Firstly, the fixed-root algorithm deriving from Eqs. (12) is applied in the case of the HCD-SU formulation. Three discretizations counting 496, 1691 and 6762 3-node triangular elements respectively have been considered. The corresponding load-displacement curves are plotted in Fig. 9a. The localization zones at the end of each simulation are depicted in Fig. 9b.
It can be noticed that the crack paths start developing correctly from the notches, but then, as the distance from the respective root element increases, a loss of continuity occurs in all the three simulations. As it can be observed in Fig. 9b, the localizations zones are not always defined by simple bands of elements, which means that stress-locking and spurious cracking take place. Consequently the knowledge of the crack-tip position is lost and the constitutive response is not well evaluated. In addition, numerical issues are encountered already at early stages in particular when fine meshes are adopted. These drawbacks have already been reported in [25] and may constitute a limitation to the applicability of the global tracking scheme.
If the variable-root algorithm derived from Eqs. (13) is adopted, the position of the crack-tip is always available during the analysis. This information allows to impose the continuity of each crack path by separating it into an active part and a potential part without the use of any further strategy. As shown in Fig. 10a, the numerical simulations fit pretty well with the experimental result in terms of load-displacement curves, denoting a good mesh-size independence. The localization zones, depicted in Fig. 9b, consist in a fully developed crack-band initiating at the upper notch and a partially developed crack-band starting from the lower notch.
The mathematical formulation is now discussed. A comparison between the heat conduction-like problem and the heat convection-diffusion-like problem is drawn in case of the variable-root global tracking. The load-displacement curves and The global structural response shows minor differences between the two simulations, in particular it seems that the HC approach allows a better energy dissipation in the final stage of the analysis. By observing the crack trajectories, it appears that an asymmetrical propagation takes place for the HCD-SU problem, whereas two almost identical crack paths are found for the HC problem, which is coherent with the fact that the stress distribution is symmetrical with respect to the vertical axis. This discrepancy may be due to the different physical meaning of the mathematical formulation. In particular, for the HCD-SU problem the tangent vector field represents a fluid velocity, deduced from the principal stress directions. Therefore, the solution is affected by the sense given to the velocity field. Such operation may lead to less accurate results with respect to the HC formulation if not properly done, especially in presence of singular points or high gradients in the stress distribution.

Four-point shear test on a single-edge notched beam
The second example is the numerical simulation of four-point shear of a single-edge notched beam. In literature this problem has been addressed by several authors, see e.g. [31][32][33]22]. At first, the same SEN beam configuration studied in Section 2.3 is considered. An imposed vertical displacement is applied centrally to the loading platen. A curved crack path is expected to initiate at the right lower corner of the notch.
The fixed-root algorithm is compared to the variable-root algorithm when the HCD-SU is adopted. Fig. 13a and b depict the total applied load F versus the Crack Mouth Opening Displacement (CMOD) for the two cases. In Fig. 14, the same comparison is made in terms of load-Crack Mouth Sliding Displacement (CMSD). It can be noticed that in the case of the SEN beam the global structural response is less sensitive to the evolution of the thermal-like field. This can be explained by considering the fact that the structure undergoes a sudden failure mechanism as soon as the peak-load is reached. As a consequence of this brittle behavior, a fully developed localization zone is observed once equilibrium is re-obtained. Thus, the reference isovalue at the end of the simulation does not differ substantially from the initial one. Nevertheless, in the case of fine meshes, the difference between the two approaches is more important: indeed, spurious cracking arises for the fixed-root algorithm, as it can be seen in Fig. 15a. The band of elements that undergo failure thus does not define a strictly      13 continuous crack path, although it reproduces pretty well the shape of the experimental failure surface. On the contrary, the continuity requirement is recovered by the variable-root algorithm (see Fig. 15b).
The capability to well reproduce the failure mode is now analyzed once again for the two reference formulations considered here, i.e. the HC and HCD-SU problems. This is done by means of the variable-root algorithm applied to the intermediate discretization containing 3672 elements. The global structural response is plotted in Fig. 16 in terms of load-CMOD and load-CMSD curves, while the resulting cracks paths are depicted in Fig. 17. A good agreement between the simulations can be observed, although in the case of the heat conduction-like problem a loss of uniqueness of the thermal solution has been encountered in the proximity of the peak-load. This numerical issue has obliged the authors to stop the tracking procedure in order to complete the analysis (see Fig. 18).

Conclusions and perspectives
In this paper a modified global crack-tracking strategy has been applied to the E-FEM simulation of quasi-brittle materials. Firstly, the problem of tracing the envelopes of a vector field has been addressed by comparing different mathematical formulations. The Streamline Upwind discretization of the heat convection-diffusion-like problem (HCD-SU) has been chosen as the privileged candidate for this task. Secondly, a revised algorithm for the representation of the propagating discontinuity has been presented. The main ingredient is the updating of the material point used to trace the potential part of the crack, instead of assuming a fixed-root for its representation. These issues have been discussed in detail by means of representative numerical examples and the following conclusions can be drawn: Independently from the mathematical formulation of the thermal-like problem, a fixed-root scheme does not prevent the loss of continuity of the crack path due to the rotation of the isovalues. As a consequence, spurious cracking appears and both stress-locking effects and numerical instability issues occur; The variable-root scheme is able to provide a C 0 -continuous crack path at no additional computational cost and without integrating any further technique; A good agreement with the experiments is found when coupling the variable-root algorithm to the HCD-SU formulation, although in the case of coarse meshes its precision is reduced with respect to the heat conduction-like approach (HC). However, from the authors' experience, good results have been obtained for all the spatial discretizations adopted in the analysis; Attention must be paid to the sense of the tangent vector field, especially in the case of strong gradients of the stress distribution. This consideration is important since Dirichlet boundary conditions can be applied only on the portions of the boundary where the velocity field is directed inward the domain; The HCD-SU formulation provides better stability performances with respect to the HC formulation. The latter is strongly dependent on the numerical conductivity parameter, which may be not sufficient to guarantee stable solutions and therefore to perform the global tracking procedure.
The applicability of the variable-root HCD-SU tracking strategy to three-dimensional problems has not been investigated yet, in particular the variable-root algorithm could be object of further studies. The extension of the proposed approach to branching scenarios should be also evaluated, such as the possibility to handle intersecting cracks.