Modelling of fracture problems in quasi-brittle materials by the E-FEM Modelagem de problemas de fratura em materiais quase-frágeis pelo E-FEM

Resumo In this paper a numerical model with strong discontinuities is presented to address fracture problems in quasi-brittle materials. A non-symmetrical statically and kinematically consistent formulation is implemented. The strong discontinuity in the displacement field is represented using the el - emental enrichment finite element method (E-FEM). In other words, the strong discontinuity is introduced into the finite element and the additional degrees of freedom are condensed at the element level, allowing the implementation into existing computational codes. Two constitutive models are used to analyze the behavior of the cracked zone, linear and exponential. The exponential model results are closer to those obtained in ex- perimental data and representative numerical simulations than the linear model


Introduction
The formation and propagation of cracks is a phenomenon observed in many materials used in engineering, such as concrete, metals, ceramics and rocks. This process happens due to the formation of zones with strain localization where the concentration of damage and other inelastic effects occurs. Crack propagation occurs in arbitrary directions that can be influenced by the geometry of the structure, boundary conditions, heterogeneity or local defects of the material [1]. Therefore, one of the major challenges in discretization of crack propagation problems is the fact that the discontinuities propagate through the structure in arbitrary directions as loading evolves. Several models to represent the crack propagation can be found in the literature, as discrete cracks (see [2], [3], [4], and [5]) and smeared cracks (see [6], [7], [8] and [9]). However, due to the complexity and limitations of those models, methods that are mesh independent and allow the propagation of cracks without remeshing are also being considered.
In the context of finite elements, two new approaches have been proposed. The extended finite elements (X-FEM), based on nodal enrichment or interpolation functions enrichment associated with existing nodes ( [10], [11]) and the elemental enrichment (E-FEM), involving finite elements with internal degrees of freedom to represent strong or weak discontinuities ( [12], [8], [13], [14]). One of the major advantages of considering elemental enrichment (E-FEM) is the local enrichment feature, i.e., additional degrees of freedom are eliminated from the global solution by static condensation. The technique allows the implementation in existing finite element codes making few modifications, besides presenting advantages in terms of computational cost and convergence when compared to extended models (X-FEM), as observed by [15]. This paper presents a strong discontinuity model proposed by [12]. The model is nonsymmetrical (SKON) according to the nomenclature of [15] and belongs to the elemental enrichment elements (E-FEM). The model was implemented in the code METAFOR (METAFOR is a commercial FE code developed in the Liège University, see [16]).
The outline of the rest of this paper is as follows. Section 2 presents the Variational Principle governing the problem. Section 3 describes the approximation of the model by the finite element method. Section 4 presents the asymmetric model implemented. Section 5 provides the crack propagation criterion and constitutive relation. Section 6 shows the results and the conclusions are discussed in chapter 7. ows the results and the conclusions are discussed in chapter 7.

Variational principle
Strong discontinuity models simulate the relation between forces through the crack and the opening of the crack (discontinuity of the internal displacement field of the element). For these models, the variational principle that represents the problem should include the relation between the transmitted stresses versus crack opening [8]. In the implemented model [12], the Hu-Washizu variational principle for incorporation of discontinuities in the displacement field is used. In this principle the displacement u, strain ε, and stress σ fields are independent of each other. These fields are defined in a V domain, where volume forces are applied. The surface is divided into two parts: S u , where the essential boundary conditions, , are applied; and S t , in which the natural boundary conditions are applied (Figure 1a) [15]. The principle can be extended to a body with an internal interface S, Figure 1b), which divides the domain and the boundary conditions into two parts. A field of surface forces, t j , appears on the inner surface. This field is a function of the discontinuity of displacements through the internal interface [8].
The field equations that govern the problem can be coupled on a variational principle according to equation (1): (1) in which δ represents variation, σ(ε) the stress obtained from the constitutive relations, ∇u the displacement gradient.

Modelling of fracture problems in quasi-brittle materials by the E-FEM
The stationary condition of this principle provides the relations of strain-displacement, the constitutive relation, the differential equation of equilibrium and the static boundary condition, according to equations to , respectively. In addition, it provides geometric and mechanical boundary conditions as natural boundary conditions [17]. (2) where n is the outward normal vector to the boundary.

Finite element method approximation
The numerical modeling of strong discontinuities on solids requires the use of a formulation that correctly represents the discontinuity in the displacement field, considering the independence between the fields of stress, strain and displacement. In this item, a general formulation will be presented within the context of finite elements, based on the works of [8], [15] and [18]. In order to represent the discontinuity of displacements in the internal interface, the displacement field is represented by the sum of the continuous and discontinuous portions representing the relative motion between the two parts of the domain separated by the discontinuity [19]. Therefore, the field of displacements with discontinuities is: (6) where N are the standard shape functions, which assume the unitary value on its respective node and zero elsewhere, d N are the nodal displacements, N C are the shape functions of the additional displacement terms and d C are the additional displacement modes. The strain field can be determinate as: (7) where B are the derivatives of the standard shape functions (N), G is the matrix that contains the additional shape functions for strain, e is the vector that contains the additional strain modes [8].
Based on the variational principle defined in , the stress field, strain field and displacement field can be defined independently. Therefore, the stresses are calculated as: (8) in which S is the matrix that contains the stress interpolation function and s is the vector that contains the stress parameters.
Replacing the equations (6) to (8) in (1) and considering ∇(Nd N ) = Bd N and ∇(N C d C ) = B C d C , we have: (9) In that, f EXT are the usual external forces and f c are the additional external forcers defined by: (10) and (11) For loads applied outside the region with additional interpolation, f c = 0.
In the equation (9) B C are the additional displacement interpolation functions, G are the additional strain term shape functions, which can be defined independently, in the case where additional displacements and strains are defined independently. Due to the independence of the variables, we can obtain the discretized equations: In order to linearize the dependence of σ in relation to d N and e, the formulation is changed to incremental form (rates). For a given state, the linearized stress-strain relation is: (16) where D ≡ δσ / δε (constitutive matrix of the material). Changing the equations (12) to (15) for the rate form and replacing (16) in these equations we obtain the following set of equations: The interpolation of stresses and strains can be discontinuous. Therefore, the stress and strain parameter can be associated with only one finite element. The same happens for the additional displacement parameters (e, s, d C ). They can be condensed at the element level and global equations contain then only the degrees of freedom relative to the standard displacement, d N [8]. Therefore, the equations (12) to (15) and (17) can be rewritten to a finite element that occupies a volume V e and the external forces f EXT are replaced by the elemental contribution of internal forces f i EXT . According to [15], the formulation presented can be particularized in three cases: Kinematically optimal symmetric elements (KOS), Statically optimal symmetric elements (SOS), and Kinematically and statically optimal non-symmetric elements (SKON). The first describes the kinematic aspects satisfactorily, but leads to an inappropriate relation of stresses in the crack, the second considers the continuity of stresses through the internal interface, but does not guarantee kinematic continuity and the latter presents C. Z. S. MARASCA | E. BITTENCOURT | V. M. R. D. BESSA a better performance by using a continuity condition of natural stresses and fairly well represent kinematic continuity. A kinematically and statically optimal non-symmetric elements (SKON) model is implemented because this formulation presents more robust and reliable results than the others [19], for more information see [16].

Asymmetric model implemented
In this paper the asymmetric model proposed by [12] is implemented to represent the strong discontinuities, with the following characteristics: n Consider the entire element as a minimum domain in the localization of strains, instead of working at the integration point level. n The strain localization within the finite element is considered as a displacement discontinuity line incorporated in the element domain. n Two constitutive relations are defined to represent the material behavior when the localization is started at the element. A stressdisplacement relation for the discontinuity line, related to fracture energy, and a stress-strain relation for the element domain. n The elements resulting from this formulation are non-conforming. n This formulation represents the global effects of locating strains on a structure. Hence, it is not possible to obtain a detailed description of the stress field near the localization zone. In order to better understand the implemented asymmetric formulation, initially a simple case of a bar finite element with two nodes is analyzed, Figure 2a). Before the opening of cracks, displacements are obtained according to equation (18). (18) where U is the vector that contains the increments of nodal displacement and N contains the interpolation functions. The dashed line 3 in the Figure 2a) represents this interpolation. When the crack opens, a discontinuity of displacements is considered in center A, which divides the element into two subdomains V 1 and V 2 . It is assumed that V 2 suffers an increment of rigid body displacement (e) regarding V 1 . In order to obtain the same strain for both subdomains, the interpolations of lines 1 and 2 shown in Figure 2a) are adopted for V 1 and V 2 . Therefore, the displacements for each subdomain considering the line of discontinuity becomes: (19) and (20) where ϕ =(0 1) T . Deriving the equations (19) and (20), the same incremental strain is obtained. Therefore, for any point into subdomains V 1 and V 2 .

(21)
For the two-dimensional case (Figure 2b)) the concepts presented above can be generalized considering: (22) where e' contains the rigid body motion components associated with the localization line, evaluated at the center of the element, in the local system (x', y'): (23) and R is the rotation matrix of the local Cartesian system of discontinuity (x', y') to the coordinate system of the element. Moreover, ϕ becomes a matrix defined by: (24) where each submatrix ϕ n is a matrix of 2x2 dimension, defined according to: where I represents the identity matrix.  which is the displacement that causes strains.
For two points near the line of discontinuity, one on V 1 and another on V 2, the displacement is given by: The strain for the two-dimensional case is: After the crack is open, it is necessary to consider an internal equilibrium condition for the discontinuity line. Starting from equation (1), see more details in [16], we have: The matrix P selects the stress components that will be transmitted through the crack. The constitutive relation for the localization line is defined as: where E t (u c ) is adjusted according to the fracture energy, G f , and G t is related to the shear modulus.
For the element domain: (32) In that t and t+Δt represents pseudo-times relatives to the load increments, where t is the previous increment and t+Δt is the current load increment (total). D e is the elastic constitutive matrix of the material. is the consistent tangent stiffness matrix. It is observed that the terms in parentheses make the matrix asymmetric, which can be a problem for implementation in pre-existing codes that use only symmetric solvers. In the implemented model, only the symmetric part of the stiffness matrix of equation (41) was used because, for the analyzed problems, the difference in disregarding the asymmetric part is irrelevant. This is supported by [12], that observed that the symmetric matrix presented satisfactory results, although it required more iterations to achieve convergence. Even using only the symmetrical part of the stiffness matrix in equation (41), the presented method produces results more consistent and robust than methods based on the symmetrical formulations -in crack propagation, the conditions of continuity of stresses at the internal interface and rigid body motion between the sides of the element separated by the crack are still guaranteed. Therefore, at the global level the system of equations to be solved is given by: To obtain the equilibrium at time t t + ∆ it is necessary to iterate at the element level the equation and, at the global level (structure), using load control for the ith iteration:

(43)
The final equilibrium of the equations given in (43) is obtained using some iterative method. In this work we used the Newton-Rhapson method. The stiffness matrix is obtained by integration with 2x2 Gauss points. After the crack opens on an element, the element is considered as the minimum domain, i.e., all properties are calculated in the center of the element, using integration with 1 Gauss point. This may lead to null strains modes (spurious modes) in the use of distorted finite elements in the cracked region. In order to avoid this problem, it is interesting to use finite elements able to overcome this problem, such as the QMITC element [8].
To obtain the equilibrium within the element, equation (30), presents some particularities. Reference [20] shows an algorithm to obtain this equilibrium when a localization line is opened inside an element: a) Consider the crack opening equal to the opening of the previous global iteration k=0 local iteration counter

Constitutive relations and crack propagation criterion
The constitutive relations used in this work are presented in Figure 3.
The linear elastic constitutive relation (Figure 3a)) is used to represent the behavior of the intact material (without cracking) and to represent the discharge behavior of the material in the non-cracked region after the cracking process has begun. This relation is presented in equation (50).

(50)
where σ represents the stress in the domain, D e the elastic constitutive matrix and ε the strain. The linear softening constitutive relation used for the cracked line, Figure 3b), is the Hillerborg model [21]. The area under the graphic in the Figure 3b) represents the fracture energy. The fracture energy (G f ) and tensile strength (f t ) are characteristic of the material, then the maximum crack opening can be obtained through equation (51), The equation that represents this constitutive relation is given by: The exponential constitutive relation, schematized in Figure 3c), is based on the model presented by [12]. The equations that represent the exponential softening curve is: (53) A limit is imposed for the factor exp(-aw) = 0.05 in order to obtain a maximum crack opening value (w máx ). α is taken equal to 1.05. In the implemented model, the failure criterion proposed by Rankine was used, then failure occurs by fracture in planes of maximum tensile stress. Therefore, according to this criterion, when

Numerical examples
In this section, we present the numerical experiments to illustrate the application of the presented methodology.

Four-point bending test
In this first example, a simply supported beam loaded by two symmetric disposed forces, tested experimentally by [23], is analyzed, as seen in Figure 4. In this case, shear effects are eliminated between the points of loading application and beam is under pure bending in the region. The beam was discretized by four different finite element meshes, according to Table 1 (elements at the crack line include enriched interpolation functions described above). In Figure 5, M1 e M4 meshes are shown. In Table 2 material and geometrical data are given. Initially, the influence of the softening law on global behavior is analyzed. In Figure 6, the softening functions used are presented. In the linear model the maximal opening is smaller than in the exponential model. As it will be seen later, this difference changes global behavior. In Figure 7, load versus displacement results for the two softening functions are shown for mesh M4. Displacement is measured at loading point. Lotfi and Shing [23] analyzed the same example using two types of fracture models: one incorporated (R1) and another distributed (R2). In Figure 7, R1 case is also shown. As seen, linear softening leads to larger peak loading. Post-peak behavior is also brittler because less energy is available in the material for larger openings. In the exponential case, post-peak descent is not so steep and closer to values obtained by [23].
In Figure 8, a study of the mesh influence is shown. Meshes M1, M2, M3 and M4 (see Table 1) are considered. It can be concluded that results tend to converge to one result with mesh refinement, showing the objectivity of the formulation is fulfilled. Results R1 and R2 of the reference are also shown.

Three-point bending test
A classic test to analyze the crack propagation in mode I is the notched beam tested under three-point bending. To validate the methodology implemented in this paper, the beam tested experimentally by [24] was used as reference. Figure 9 shows the beam geometry and the boundary conditions. The beam was discretized with three different finite element meshes, according to Table 3. Figure 10 shows the meshes used with discretization detail at the notched region.

Figure 4
Beam geometry and boundary conditions -four-point bending test   Table 4 shows the input data used for this example and correspond to experiments done by [24].
For the non-cracked zone a linear elastic constitutive relation was used. In the cracked zone, the influence of the two different softening laws, linear and exponential, are again investigated. Figure 11 presents a comparison of the results obtained for the two softening laws used, for the most refined mesh (M3), and also the experimental results. As can be observed the exponential soften-ing presents results very similar to the real behavior of the tested beam. In the linear model, the maximum stresses are overestimated and the post-peak curve is steeper than the exponential model. This behavior is similar to the previous case analyzed.
A comparison among different finite element meshes is presented in Figure 12. A finer refinement of the mesh was done in the zone near the crack. It is noted that even for a coarse mesh (M1) the

Figure 5
Discretization of the meshes M1 and M4 with detail in the region of the notch -four-point bending test

Figure 6
Constitutive models to the crack line

Conclusion
The aim of this work was to implement a crack model with strong discontinuities in code METAFOR, in order to analyze the behavior of structural elements in the post-peak stage. The implemented model was based on the study proposed by [12]. Kinematically and statically optimal non-symmetric elements (SKON) formulation was implemented and the QMIT finite element (Quadrilateral with Mixed Interpolation of Tensorial Components) was used. In this work a quadrilateral bilinear element in addition to the symmetrical part of the stiffness matrix is used.
Analyzing the results presented in section 6 we could verify the ability of the model to correctly capture the behavior of the material in the regions with and without cracking. It was possible to verify that the softening laws have big influence on the behavior of the structure. The exponential softening model [12] appeared to be the best solution because the peak stress was correctly estimated and the softening curve presented a more

Figure 9
Beam geometry and boundary conditions -three-point bending test  realistic behavior, when compared with experimental data. In the linear model [21], the peak loads were overestimated when compared to experimental data. In addition, the linear model presented a more brittle softening branch than is actually observed for quasi-brittle materials.
In general, the results obtained presented a good agreement to the experimental results as well as numerical results obtained by crack models presented in the literature. Therefore, it can be stated that the implemented model is suitable for the simulation of quasi-brittle

Figure 11
Load x Deflection of the beam with linear and exponential softening

Figure 12
Load x Deflection of the beam with different mesh sizes

Figure 10
Discretization of the meshes with detail in the region of the notch -three-point bending test Modelling of fracture problems in quasi-brittle materials by the E-FEM materials. It is worth mentioning that this model has as advantage the use of relatively coarse meshes, besides the possibility of implementation in existing finite element codes.