Engineering Fracture Mechanics 153 (2016) 151–170 Contents lists available at ScienceDirect Engineering Fracture Mechanics journal homepage: www.elsevier .com/locate /engfracmech On the use of finite elements with a high aspect ratio for modeling cracks in quasi-brittle materials http://dx.doi.org/10.1016/j.engfracmech.2015.12.026 0013-7944/� 2015 Elsevier Ltd. All rights reserved. ⇑ Corresponding author. Tel.: +55 14 31036000. E-mail address: omanzoli@feb.unesp.br (O.L. Manzoli). 1 Tel.: +55 11 31036000. 2 Tel.: +55 11 30915677. Osvaldo L. Manzoli a,⇑, Michael A. Maedo a,1, Luís A.G. Bitencourt Jr. b,2, Eduardo A. Rodrigues b,2 a São Paulo State University – UNESP/Bauru, Av. Eng. Luiz Edmundo C. Coube 14-01, CEP 17033-360 Bauru, SP, Brazil bUniversity of São Paulo, Av. Prof. Luciano Gualberto, Trav. 3 n. 380, CEP 05508-010 São Paulo, SP, Brazil a r t i c l e i n f o Article history: Received 30 July 2015 Received in revised form 23 November 2015 Accepted 1 December 2015 Available online 29 December 2015 Keywords: Mesh fragmentation technique Interface solid finite element Arbitrary cracks Quasi-brittle materials Tension damage model a b s t r a c t A new technique for modeling cracks in quasi-brittle materials based on the use of inter- face solid finite elements is presented. This strategy named mesh fragmentation technique consists in introducing sets of standard low-order solid finite elements with a high aspect ratio in between regular (or bulk) elements of the mesh to fill the very thin gaps left by the mesh fragmentation procedure. The conception of this strategy is supported by the fact that, as the aspect ratio of a standard low-order solid finite element increases, the element strains also increase, approaching the same kinematics as the Continuum Strong Discontinuity Approach. As a consequence, the analyses can be performed integrally in the context of the continuum mechanics, and complex crack patterns can be simulated without the need of tracking algorithms. A tension damage constitutive relation between stresses and strains is proposed to describe crack formation and propagation. This constitutive model is integrated using an implicit–explicit integration scheme to avoid convergence drawbacks, commonly found in problems involving discontinuities. 2D and 3D numerical analyses are performed to show the applicability of the presented technique. Relevant aspects such as the influence of the thickness of the interface elements and mesh objectivity are investigated. The results show that the technique is able to predict satisfac- torily the behavior of structural members involving different crack patterns, including multiple cracks, without significant mesh dependency provided that unstructured meshes are used. � 2015 Elsevier Ltd. All rights reserved. 1. Introduction Although a large number of numerical models are available in literature for modeling cracks in quasi-brittle materials (e.g. concrete and rocks), each one with their own advantages and limitations, this topic is still a challenge, mainly when numerical simulations of multiple cracks in 3D analysis is desirable. In the context of the Finite Element Method (FEM), the available models can be usually classified into the ‘‘smeared crack model” or the ‘‘discrete crack model”. Both approaches were introduced for modeling cracks in concrete in the late 1960s, the http://crossmark.crossref.org/dialog/?doi=10.1016/j.engfracmech.2015.12.026&domain=pdf http://dx.doi.org/10.1016/j.engfracmech.2015.12.026 mailto:omanzoli@feb.unesp.br http://dx.doi.org/10.1016/j.engfracmech.2015.12.026 http://www.sciencedirect.com/science/journal/00137944 http://www.elsevier.com/locate/engfracmech Nomenclature n; s; t components of the local Cartesian cordinate system n unit vector – normal to the interface FE h height (thickness) of the interface FE B strain–displacement matrix t interface stresses d nodal displacement vector d; ~d scalar damage variables f t material tensile strength A area of the base of the tetrahedron FE q; r stress and strain-like internal variables Gf fracture energy xðiÞ vector that collects the coordinates of the node (i) xðiÞn ; xðiÞs ; xðiÞt coordinates of the node (i) uðiÞ n ; uðiÞ s ; uðiÞ t displacement components of the node (i) sut relative displacement between the node (1) and its projection on the element base (10) sutn; suts; sutt components of the relative displacement C fourth order elastic tensor ~Ctan algorithmic tangent operator E Young’s modulus m Poisson’s ratio e strain tensor ~e; ê components of the strain tensor r; �r nominal and effective stress tensors rnn; �rnn nominal and effective stress components normal to the base of the interface FE /; �/ damage criterion functions �s equivalent effective stress A softening parameter 152 O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 former by Rashid [1] and the latter by Ngo and Scordelis [2]. Basic characteristics of both models and a comparison between them can be found in [3–5]. In the smeared crack model, a fixed Finite Element (FE) mesh is used and the effects of the cracks are incorporated into the constitutive model (stress–strain relation) adopted to describe the material, which is usually nonlinear with a softening regime after the failure (maximum) criterion is reached. This approach is also known as a continuum type representation and the crack zone is considered to be distributed in a certain strain localization (continuum) region of the solid. Models of this class can be formulated using fixed or rotating crack approaches [6–8,4,9]. It is known that the main drawback of stan- dard smeared crack formulations is the mesh sensitivity. To remedy problems of mesh size dependency presented by the models of this class, several regularization procedures have been proposed [10]. Some research papers [11,12] have proposed fracture mechanics concepts that lead to fracture energy release regularization. However, the loss of objectivity associated with the deformation pattern is still observed. Other strategies for dealing with the mesh dependency problem have been addressed using non-local models [13,14] and gradient-enhanced approaches [15]. Discrete models, also referred to as discontinuous models by some authors, are characterized by the introduction of dis- placement or strain discontinuities into standard finite elements to represent cracks. In general, these methods can be divided into two groups, defined by the location of the displacement discontinuities. Hence, one group gathers the approaches where cracks are embedded into existing bulk finite elements [16–20], and in the other group are the approaches with interface elements in which the cracks are expected between standard bulk finite elements [21–26]. A combination of the two groups is proposed by Radulovic et al. [27]. Models of the former group are usually derived using a nodal-based for- mulation (e.g. eXtended Finite Element Method – X-FEM) [20,28–30] or element-based formulation (e.g. Strong Discontinu- ity Approach – SDA) [31,17,32,18,33,34]. A comparative study of these approaches is presented by Oliver et al. [35]. Both methods require techniques to track the crack path during the analysis in order to provide information about the position of the crack surfaces for the discontinuous kinematic enrichment. These techniques are relatively simple to represent few cracks in 2D analyses, but can be very complex and even unsuitable for problems involving multiple crack surfaces in 3D analyses [36–38]. In the X-FEM, the level set method introduced by Osher and Sethian [39] has been used to track cracks. The necessity of crack tracking schemes can be avoided by using discrete approaches based on the use of interface ele- ments (with a zero thickness, e.g. [21–23,40–44]). In this approach the kinematics of a crack is described by a displacement discontinuity (crack opening width) and a cohesive law is used to describe the material degradation. The use of this kind of approach allows complex crack patterns (crack branching, intersecting cracks and fragmentation) to be simulated automatically. O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 153 In this paper an alternative discrete approach based on the use of interface solid finite elements for modeling the forma- tion and propagation of cracks is presented. These elements with a high aspect ratio combined with a continuum damage model to describe their behavior was introduced by Manzoli et al. [45]. In this work the authors show that, as the aspect ratio of a standard three-node triangular element increases (ratio of the largest to the smallest dimension), the element’s strain also increases, approaching the same kinematics as the Continuum Strong Discontinuity Approach (CSDA) [46,47]. As a consequence, a standard continuum constitutive model, which tends toward a discrete constitutive relation as the aspect ratio increases, can be applied to describe the interface behavior. Details about the equivalence between the strain field of the interface elements and the strong discontinuity kinematics are also demonstrated in Manzoli et al. [45]. The important concept that is used from this equivalence is that bounded stress can be obtained from unbounded strains. Recently, this approach was successfully used in different applications, as for example: modeling of formation and propaga- tion of 3D cracks networks induced by shrinkage in drying soils [48] and 2D modeling of bond slip behavior of reinforcing bars embedded in concrete [49]. Hence, the present approach consists in introducing standard low-order finite elements with high aspect ratio in between all the regular (bulk) finite elements (fragmentation procedure). These interface elements define all potential crack paths, and their behavior is described by a continuum damage model. Compared with available models based on the introduction of special interface elements (e.g. cohesive zone models and zero-thickness interface elements), the proposed approach presents some advantages, such as: � It is not necessary to develop a discrete constitutive relation (stress–displacement relation), nor special integration rules to obtain the internal forces, as required by zero-thickness interface elements. � The interfaces can be modeled entirely in the continuum framework, composed of standard finite elements usually avail- able in FEM codes and continuum constitutive models. � The use of very high stiffness to reduce the artificial compliance during the elastic regime (normally required in zero- thickness interface elements) is not necessary, avoiding numerical instabilities. � The proposed discontinuous approach can become available in existing nonlinear FEM codes, practically just by perform- ing a few changes in the pre-processing phase. This paper is organized as follows. In Section 2 the idea behind the proposed mesh fragmentation technique is presented. The formulation of the interface solid finite elements for 2D and 3D cases is detailed in Section 3. In Section 4 the tension damage model used to describe the behavior of the interface elements is described. In Section 5 the proposed mesh fragmen- tation technique is assessed in 2D and 3D case studies. Finally, some concluding remarks are given in Section 6. 2. Overview of the proposed mesh fragmentation technique The proposed technique, hereafter called mesh fragmentation technique, is based on the use of interface solid finite ele- ments with a high aspect ratio (proposed by Manzoli et al. [45]) which are inserted in between standard (bulk) finite ele- ments of a finite element mesh. In this work, these interface elements are responsible for describing the crack formation and propagation via an appropriate continuum damage model. Fig. 1 illustrates the main steps of the proposed mesh fragmentation technique for 2D and 3D problems, which can be summarized as follows: 1. Generation of the standard FE mesh to be fragmented (Fig. 1(a)). 2. Separation of the finite elements by introducing gaps in between them (Fig. 1(b)). 3. Insertion of interface elements with a high aspect ratio in between the bulk finite elements (Fig. 1(c)). It is important to note in Fig. 1(b) that the gaps are in exaggerated scale for illustration. These gaps are usually very small, with a thickness around 1% of the typical size of the regular elements. The choice of the GAPs thickness is still not a closed issue. However, based on the authors’ experience, the assumption of 1% of the typical size of the regular elements seems to be a very reasonable recommendation, provided that the size of the regular elements has been chosen to accurately capture the (elastic) stress field prior to the crack formation. Of course, smaller regular elements produce more accurate responses. The thickness is limited only by the computer precision. In addition to that, triangular/tetrahedron finite elements are used for both the standard FE mesh and the interface elements introduced during the fragmentation process (see Fig. 1). With the proposed strategy, cracks can only propagate along the interface elements. A reasonable mesh refinement is nec- essary to capture the failure process. For problems in which the region where cracks are expected is known a priori, the mesh fragmentation technique may be applied only in this region of interest. This approach is very attractive since the steps men- tioned above are straightforward in implementation, giving place to an additional pre-processing stage. The mesh fragmentation approach is completed by a continuum tension damage model formulated to describe the for- mation and propagation of cracks along the interfaces. As is well known, the constitutive tangent operator may become sin- gular in the strain softening regime, and, as a consequence, the solution of the resulting systems of nonlinear equations using a fully implicit schememay not be achieved. According to Oliver et al. [50] this problemmay also be present in the numerical Fig. 1. 2D and 3D mesh fragmentation process: (a) generation of the standard FE mesh to be fragmented; (b) separation of the finite elements (with an exaggerated scale factor for clarity); (c) insertion of interface elements (depicted in gray); and (d) detail of interface elements between regular elements. 154 O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 modeling of material failure even when a powerful continuation method to pass structural points is used (e.g. arc length methods to transverse limit and turning points). To address this problem, the implicit–explicit (IMPL–EX) integration scheme proposed by Oliver et al. [50,51] is used in this work for the integration of the stress–strain relation of the damage model. 3. Interface solid finite element To describe the main features of the interface solid finite elements in 2D and 3Dmodeling, standard three-node triangular finite element and four-node tetrahedron ones, as illustrated in Fig. 2, are considered. The geometry of these elements can be characterized by the position of their nodes according to a local Cartesian coordinate system ðn; s; tÞ, defining the unit vector, n, normal to the element base, and the height, h, given by the distance between the node ð1Þ and its projection on the ele- ment base (10). Following the standard finite element approximations, the strain field of the solid finite element can be expressed by: e ¼ Bd; ð1Þ where B is the strain–displacement matrix and d is the nodal displacement vector of the element. 3.1. Decomposition of the strain tensor To show the similarity between the kinematics provided by these elements, when h tends to zero, and that associated to the CSDA, the strain tensor is divided into two parts [45]: e ¼ ~eþ ê; ð2Þ where ê collects all the components of the strain tensor that depends on h, and ~e contains the rest of the components. Thus, the components that depend on h can be written as: (a) (b) Fig. 2. Interface solid finite elements: (a) three-node triangular element and (b) four-node tetrahedron element. O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 155 ê ¼ 1 h ðn� sutÞs; ð3Þ where ð�Þs corresponds to the symmetric part of ð�Þ; � denotes a dyadic product, and sut is a vector that collects the com- ponents of the relative displacement between the node ð1Þ and its projection on the element base (10). The total strain tensor, given by Eq. (2), can then be rewritten as: e ¼ ~eþ 1 h ðn� sutÞs|fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} ê : ð4Þ As can be noted from this decomposition, when h tends to zero, the component ~e remains bounded, while the ê compo- nent is no longer bounded. As a consequence, in this situation, the element strains are related almost exclusively to the rel- ative displacement between the node ð1Þ and its projection on the element base (10). As described by Manzoli et al. [45], in the limit situation ðh ! 0Þ, the relative displacement sut becomes the measure of a displacement discontinuity (strong dis- continuity), and the structure of the strain field in Eq. (4) corresponds to the typical kinematics of the CSDA. Therefore, based on the same concepts of the CSDA, it can be stated that bounded stress can be obtained from unbounded strains by means of a continuum constitutive relation. The equivalence between the strain field of the interface finite elements (when h ! 0) and the strong discontinuity kinematics is detailed by Manzoli et al. [45]. In this work a continuum tension damage model is used to describe the dissipation effects associated to crack opening (see Section 4). 3.1.1. Components of the strain tensor for the 3-node triangular finite element According to the local Cartesian coordinate system ðn; sÞ, depicted in Fig. 2(a), the nodal coordinates of the 3-node trian- gular finite element are: xð1Þ ¼ ðh; b2Þ; xð2Þ ¼ ð0;0Þ; xð3Þ ¼ ð0; bÞ: ð5Þ Thus, the components of the strain tensor of the element, given by Eq. (4), can be expressed as: ~e ¼ 1 b 0 1 2 ðuð3Þ n � uð2Þ n Þ 1 2 ðuð3Þ n � uð2Þ n Þ ðuð3Þ s � uð2Þ s Þ " # ; ð6Þ and ê ¼ 1 h sutn 1 2 suts 1 2 suts 0 " # ; ð7Þ where uðiÞ n and uðiÞ s are the components of the displacement of node i according to the ðn; sÞ system; while sutn ¼ uð1Þ n � uð10 Þ n and suts ¼ uð1Þ s � uð10 Þ s are the components of the relative displacement sut. 3.1.2. Components of the strain tensor for the 4-node tetrahedral finite element In the same way, the nodal coordinates of the 4-node tetrahedron finite element (see Fig. 2(b)) according to a local Carte- sian coordinate system ðn; s; tÞ are: xð1Þ ¼ ðh; xð1Þs ; xð1Þt Þ; xð2Þ ¼ ð0; xð2Þs ; xð2Þt Þ; xð3Þ ¼ ð0; xð3Þs ; xð3Þt Þ; xð4Þ ¼ ð0; xð4Þs ; xð4Þt Þ: Thus, the corresponding parts of the strain tensor, given by Eq. (4), can be expressed as: ~e ¼ 1 A eenn eens eenteens eess eesteent eest eett 264 375; ð8Þ and ê ¼ 1 h sutn 1 2 suts 1 2 sutt 1 2 suts 0 0 1 2 sutt 0 0 264 375; ð9Þ 156 O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 where eenn ¼ 0; eess ¼ 1 2 ðxð3Þt � xð2Þt Þuð4Þ s þ ðxð2Þt � xð4Þt Þuð3Þ s þ ðxð4Þt � xð3Þt Þuð2Þ s h i ; eett ¼ �1 2 ðxð3Þs � xð2Þs Þuð4Þ t þ ðxð2Þs � xð4Þs Þuð3Þ t þ ðxð4Þs � xð3Þs Þuð2Þ t h i ; eens ¼ 1 4 ðxð3Þt � xð2Þt Þuð4Þ n þ ðxð2Þt � xð4Þt Þuð3Þ n þ ðxð4Þt � xð3Þt Þuð2Þ n h i ; eent ¼ �1 4 ðxð3Þs � xð2Þs Þuð4Þ n þ ðxð2Þs � xð4Þs Þuð3Þ n þ ðxð4Þs � xð3Þs Þuð2Þ n � � ; eest ¼ 1 4 ðxð3Þt � xð2Þt Þuð4Þ t þ ðxð2Þt � xð4Þt Þuð3Þ t þ ðxð4Þt � xð3Þt Þuð2Þ t h �ðxð3Þs � xð2Þs Þuð4Þ s � ðxð2Þs � xð4Þs Þuð3Þ s � ðxð4Þs � xð3Þs Þuð2Þ s � ; and A is the area of the element base with unit vector n; uðiÞ n ; uðiÞ s and uðiÞ t are the components of the displacement of node i according to the ðn; s; tÞ system; sutn ¼ uð1Þ n � uð10 Þ n ; suts ¼ uð1Þ s � uð10 Þ s and sutt ¼ uð1Þ t � uð10 Þ t are the components of the relative displacement sut between node (1) and the point corresponding to its projection on the element base (10) xð10 Þ ¼ ð0; xð1Þs ; xð1Þt Þ � � . 4. Tension damage model It is well know that constitutive models based on the Continuum Damage Mechanics Theory (CDMT) are appropriate to describe the material degradation process due to crack propagation [52–55]. For this class of models, the mechanical behav- ior of a damaged material is usually described by using the notion of the effective stress, together with the hypothesis of mechanical equivalence between the damage and the undamaged material. In this work the phenomenon of crack initiation and propagation through the interface solid finite elements is described by a tension damage model. In the following, the main ingredients of this model and the integration of the stress–strain rela- tion using an implicit–explicit integration scheme are detailed. The resulting discrete relation obtained when the interface element height, h, tends to zero is also presented. 4.1. Continuum model The tension damage model is defined by the following constitutive relation: r ¼ ð1� dÞ C : e|ffl{zffl} �r ; ð10Þ where r is the nominal stress; d 2 ½0;1� is the scalar damage variable; C is the fourth order elastic tensor; e is the strain ten- sor; and the product C : e defines the effective stress tensor �r. The damage criterion defines the elastic domain and is given by: / ¼ rnn � qðrÞ 6 0; ð11Þ where rnn is the component of the stress normal to the base of the element ðrnn ¼ n � r � nÞ, q and r are the stress and strain- like internal variables, respectively, and the function qðrÞ defines the softening law. Dividing Eq. (11) by ð1� dÞ, the damage criterion can be expressed in terms of the effective stresses as: �/ ¼ �rnn � r 6 0; ð12Þ with r ¼ q=ð1� dÞ, which controls the size of the elastic domain in the space of the effective stress. Thus, a damage evolution rule can be defined in terms of the internal variable r: dðrÞ ¼ 1� qðrÞ r : ð13Þ The constitutive model is completed by the loading–unloading conditions, given by the Kuhn–Tucker relations: �/ 6 0; _r P 0; _r�/ ¼ 0; ð14Þ and the consistency condition _r _�/ ¼ 0 if �/ ¼ 0: ð15Þ Eqs. (15) and (12) lead to the following explicit evolution law for the strain-like internal variable along pseudo-time t associated with the loading process: O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 157 rt ¼ max½�rnnðsÞ; r0� s2½0;t� : ð16Þ According to Eq. (16), r takes the maximum value that the effective tensile stress �rnn reaches during the loading process, starting from the initial value, r0, equal to the material tensile strength f t . 4.1.1. Implicit–explicit integration scheme for the tension damage model In order to avoid numerical convergence problems related to the strain softening regime, an implicit–explicit (IMPL–EX) integration scheme proposed by Oliver et al. [50] is used for the integration of the stress–strain relation. This method has been widely tested for the integration of constitutive models based on damage mechanics [51] and plasticity theory [51,56]. The IMPL–EX scheme is characterized by the guarantee of convergence (robustness), regardless of the length of the load step, and low computational costs, when compared with fully implicit schemes. For the IMPL–EX scheme only one iteration per load-step is sufficient to establish equilibrium (given by the difference between external and internal forces). As a con- sequence, there is an error associated with the use of this alternative integration scheme, which can be reduced (or con- trolled) by decreasing the load step length, as can be seen in the numerous examples performed by Oliver et al. [51]. Table 1 illustrates the steps of the IMPL–EX integration scheme applied to the tension damage model at a pseudo-time step tnþ1. The integration procedure is performed in a closed form. For a given strain tensor enþ1, a corresponding stress ten- sor ~rnþ1 is evaluated explicitly and used to fulfill the balance equation and to compute the effective algorithmic tangent operator. Note that the strain-like internal variable increment evaluated at the previous step tn is used to explicitly evaluate the approximation of the current strain-like internal variable, ~rðnþ1Þ (see step (iv) in Table 1), which in turn is used to compute the current damage variable and then the stress tensor. Finally, it is also important to note in Table 1 that the corresponding algorithmic tangent operator is always positive definite, which makes the IMPL–EX methodology very appealing for this type of application. 4.2. Discrete relation of the interface solid finite element The tension damage constitutive (stress–strain) relation is used to describe the behavior of the interface elements (see previous section). Hence, using Eqs. (10) and (4), the following relation can be written: r ¼ ð1� dÞC : ~eþ 1 h ðn� sutÞs|fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} ê 2664 3775 ¼ ð1� dÞ h C : h~eþ ðn� sutÞs� � : ð17Þ However, in the limit state of the strong discontinuity, the continuum constitutive relation (Eq. (17)) becomes the follow- ing discrete relation between the interface stresses and the components of the displacement jump: t ¼ n � r ¼ ð1� dÞK̂sut; ð18Þ Table 1 IMPL–EX integration scheme for the tension damage model. INPUT: eðnþ1Þ; �rnnðnÞ ; rðnÞ; DrðnÞ (i) Compute the effective stress tensor �rðnþ1Þ ¼ C : eðnþ1Þ (ii) Compute loading–unloading conditions if �rnnðnþ1Þ 6 rðnÞ , then keep the damage threshold: rðnþ1Þ ¼ rðnÞ else update the damage threshold: rðnþ1Þ ¼ �rnnðnþ1Þ (iii) Compute the strain-like internal variable increment Drðnþ1Þ ¼ rðnþ1Þ � rðnÞ (iv) Compute the explicit linear extrapolation of the strain-like internal variable ~rðnþ1Þ ¼ rðnÞ þ DrðnÞ DtðnÞ Dtðnþ1Þ; Dtðnþ1Þ ¼ tðnþ1Þ � tðnÞ and DtðnÞ ¼ tðnÞ � tðn�1Þ (v) Update the damage variable ~dðnþ1Þ ¼ 1� qð~rðnþ1Þ Þ ~rðnþ1Þ (vi) Compute the stress tensor ~rðnþ1Þ ¼ ð1� ~dðnþ1ÞÞ�rðnþ1Þ if �rnnðnÞ > 0 �rðnþ1Þ if �rnnðnÞ < 0 ( OUTPUT: ~rðnþ1Þ; �rnnðnþ1Þ ; rðnþ1Þ; Drðnþ1Þ Compute the algorithmic tangent operator ~Ctan ðnþ1Þ ¼ @~rðnþ1Þ @eðnþ1Þ ¼ ð1� ~dðnþ1ÞÞC if �rnnðnÞ > 0 C if �rnnðnÞ < 0 ( 158 O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 where each component of this relation for 2D and 3D interface elements is provided in Table 2. To maintain the stresses bounded when the height tends to zero ðh ! 0Þ, all components of the displacement jump must tend to zero if d ¼ 0 (elastic regime). On the other hand, the damage variable must tend to 1 in the inelastic regime with sut– 0. Therefore, in the limit situation of h ! 0, the interface element presents a rigid-damage behavior. Note that, before the stresses reach the damage criterion all components of the displacement jump are precluded via penalization of 1=h. This penalization can be seen in the equation of the discrete constitutive model that emerges from the continuum constitutive model. Note that the elastic stiffness of the discrete model (Eq. (18)) is given by the continuum elastic constants (K and G) divided by h, which corresponds to a very high stiffness for very small values of h. This provides a rigid behavior, preclud- ing the evolution of the displacement jump, before the stresses reach the damage criterion. If the Poisson’s ratio is null ðm ¼ 0Þ, the discrete relation can be expressed in terms of the Young’s modulus, E, since Table 2 Compon Dim 2D 3D 4G=3þ K ¼ E and G ¼ E=2 : ð19Þ For a monotonic increase (in mode I) of the normal displacement (s _utn > 0, sutn �� t¼0 ¼ 0, suts ¼ sutt ¼ 0), the evolution of the normal stress becomes: rnnðsutnÞ ¼ ð1� dÞ1 h Esutn ¼ 1 h Esutn if sutn 6 su0tn ð1� dÞ 1 h Esutn ¼ qðrÞ if sutn > su0tn ( ; ð20Þ with su0tn ¼ hq0=E. Assuming an exponential softening law of the form qðrÞ ¼ q0e Ahð1�r=q0Þ; ð21Þ with q0 ¼ f t , where f t is the tensile strength of the material and A is the softening parameter, the mode I fracture energy, Gf , i.e., the energy dissipated in a complete degradation of the interface element in mode I, is given as Gf ¼ Z 1 0 rnnðsutnÞdsutn ¼ Z su0tn 0 rnnðsutnÞdsutn ð22Þ þ Z 1 su0tn rnnðsutnÞdsutn ð23Þ Gf ¼ ðf 2t h=2þ f 2t =AÞ=E : ð24Þ When h tends to zero, the fracture energy tends to a non-null value given as Gf ¼ f 2t AE ; ð25Þ from which one can define the softening parameter in terms of the material properties as A ¼ f 2t Gf E : ð26Þ Fig. 3 illustrates the typical curve tension versus displacement jump (opening) normal to the interface element, which is subjected to a uniaxial tension. Note that this discrete response corresponds to a rigid-damage behavior. Finally, Table 3 summarizes the equations describing of the continuum tension damage model and its resulting discrete relation that emerges when the height of the interface solid finite element tends to zero ðh ! 0Þ. From the discrete form it is possible to note that this model is able to deal with cracks evolving in mixed mode, since the damage variable degrades all components of the effective stress. For numerical applications, a very small (but non-null) value of h can be used, corresponding to the regularized kinematics of strong discontinuity [46]. The value of h must be small enough to provide a good representation of the resulting discrete constitutive model that emerges from the continuum model. ents of the discrete relation of the interface elements for 2D and 3D cases. ension of the problem Interface stresses Displacement jump Matrix K̂ t ¼ rnn rns � sut ¼ sutn suts � K̂ ¼ 1 h 4 3Gþ K 0 0 G � t ¼ rnn rns rnt 8<: 9=; sut ¼ sutn suts sutt 8<: 9=; K̂ ¼ 1 h 4 3Gþ K 0 0 0 G 0 0 0 G 24 35 Fig. 3. Resulting discrete relation of the tension damage model. Table 3 Continuum and discrete constitutive equations for the tension damage model. Continuum model Resulting discrete model Constitutive relation r ¼ ð1� dÞ�r �r ¼ C : e t ¼ ð1� dÞK̂ysut Equivalent effective stress �s ¼ �rnn �s ¼ �rnn ¼ 1 h Esutn Damage criterion �/ ¼ �s� r 6 0 �/ ¼ �s� r 6 0 Evolution law for the strain-type variable r ¼ maxs2½0;t� ½�sðsÞ; r0� r ¼ maxs2½0;t� ½�sðsÞ; r0� Damage evolution dðrÞ ¼ 1� qðrÞ r ; d 2 ½0;1� dðrÞ ¼ 1� qðrÞ r ; d 2 ½0;1� Hardening/softening law qðrÞ ¼ f te Ahð1�r=f t Þ qðrÞ ¼ f te Ahð1�r=f t Þ y 2D? K̂ ¼ E h 1 0 0 1=2 � and 3D? K̂ ¼ E h 1 0 0 0 1=2 0 0 0 1=2 24 35. O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 159 5. Numerical examples In this section, the ability of the proposed fragmentation approach in predicting different crack patterns is assessed. The interface elements are introduced in between all regular (bulk) elements of the entire mesh or, to minimize computational effort and time process, just in part of the mesh corresponding to the region of study. 2D and 3D examples are carried out using three-node triangular and four-node tetrahedral finite elements, respectively. An element thickness of h ¼ 0:01 mm is adopted for all the cases, except when the influence of this variable is being studied. The behavior of the interface elements is described using the tension damage model while the bulk elements are assumed to be linearly elastic. The softening param- eter A of the damage model is evaluated by the Eq. (26). 5.1. 2D case studies 5.1.1. Plate with hole In this first example, a plane strain analysis of a rectangular plate with a center hole loaded vertically is carried out. For symmetry reasons, only one-fourth of the plate is analyzed, as shown in Fig. 4. This elastic benchmark problem aims to study the response obtained using two different type of meshes: a standard FE mesh and other fragmented with three different thickness h ¼ 0:005; 0:01 and 0:1 mm. For the standard elements, were adopted Young’s modulus of E ¼ 30 GPa and Pois- son’s ratio m ¼ 0:2, while for the interface elements, Young’s modulus of E ¼ 30 GPa and Poisson’s ratio m ¼ 0. Fig. 5 shows the vertical normal stress distribution for all the cases studied. As can be seen, the results present good agreement and a fragmented mesh can reproduce the response obtained using a standard FE mesh in elastic regime in a sat- isfactory manner. The continuity of the displacement field is shown in Fig. 6. 5.1.2. Three-point bending beams Herein, a series of experimental tests carried out by Le Bellégo et al. [57] on three-point bending beams are numerically analyzed in order to assess the capability of the proposed mesh fragmentation technique to predict crack propagation in mode I. Fig. 7 shows the test setup of the notched mortar beams of different heights D ¼ 80, 160 and 320 mm, of length L ¼ 4D, and thickness b ¼ 40 mm. The numerical analyses are carried out in plane stress condition. The behavior of the bulk elements is assumed to be lin- early elastic with Young’s modulus E ¼ 38:5 GPa and Poisson’s ratio m ¼ 0:24. The parameters of the interface model are: Young’s modulus E ¼ 38:5 GPa, Poisson’s ratio m ¼ 0, tensile strength f t ¼ 3:6 MPa and fracture energy Gf ¼ 0:05 N=mm. Fig. 4. Plate with hole. Geometry, boundary conditions and load. Fig. 5. Vertical normal stress: (a) standard mesh, (b) fragmented mesh with h ¼ 0:005 mm, (c) fragmented mesh with h ¼ 0:01 mm, and (d) fragmented mesh with h ¼ 0:1 mm. 160 O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 The analyses are performed using three different unstructured meshes, varying their densities at the central region, as shown in Fig. 7(b). The main characteristics of these meshes are listed in Table 4, in which Nnodes is the number of nodes, NelemðrÞ is the number of regular elements, NelemðiÞ is the number of interface elements and NelemðtÞ is the total number of elements. Fig. 8 shows the deformed meshes for the case with D ¼ 160 mm at the end of the analysis, where the crack pattern can be seen. Strains are mostly concentrated in those interface elements located near the central cross section, following a dom- inant vertical crack propagation (in mode I) from the notch, with a zigzag path imposed by the mesh geometry, which is not aligned with the expected vertical path. To emphasize the crack pattern, Fig. 8(d)–(f) also show the deformed mesh without the interface elements. The experimental vertical crack could be well predicted for all the cases. The structural responses (curves of load P versus displacement d) obtained for the different sizes and meshes are shown in Fig. 9. As can be noted, experimental and numerical responses are in good agreement, and the numerical results do not pre- sent significant mesh dependency. Fig. 6. Vertical displacement field: (a) standard mesh, (b) fragmented mesh with h ¼ 0:005 mm, (c) fragmented mesh with h ¼ 0:01 mm, and (d) fragmented mesh with h ¼ 0:1 mm. (a) (b) Fig. 7. Three-point bending beam: (a) geometry and boundary conditions and (b) the three mesh refinements used. Table 4 Characteristics of the FE meshes of the three-point bending beam. FE mesh Nnodes NelemðrÞ NelemðiÞ NelemðtÞ mesh A 6861 4747 5264 10,011 mesh B 14,335 7335 12,672 20,007 mesh C 36,748 14,983 34,962 49,945 O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 161 Fig. 8. Crack growth trajectories for different mesh refinements (D ¼ 160 mm). Fig. 9. Structural responses obtained with different specimen sizes and meshes. 162 O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 In this example the influence of the thickness h adopted for the interface solid finite elements is also investigated. This study was carried out for the beam with D ¼ 160 mm and mesh B. As can be seen in Fig. 10, the results obtained using h ¼ 0:1, 0.01, and 0:005 mm are very similar and could predict well the experimental results. Fig. 10. Structural responses for interface elements with different thickness (h in mm). Fig. 11. Geometry and boundary conditions of the four-point bending beam experimentally tested by Arrea and Ingraffea [58]. O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 163 5.1.3. Four-point bending beam In this example the four-point bending beam tested by Arrea and Ingraffea [58] is numerically analyzed. The test setup is illustrated in Fig. 11. A fine fragmented mesh is considered only in the region where the crack is expected to propagate (at central region). The linear elastic material properties adopted for the standard (bulk) finite elements were: Young’s modulus of E ¼ 28:8 GPa and Poisson’s ratio of m ¼ 0:18; and for the interface elements: Young’s modulus of E ¼ 28:8 GPa, Poisson’s ratio of m ¼ 0, fracture energy of Gf ¼ 0:1 N=mm and tensile strength of f t ¼ 2:8 MPa. Plane stress condition is considered with an out-of-plane thickness of 150 mm. As illustrated in Fig. 12, three unstructured mesh refinements have been considered in order to study the mesh depen- dency, ranging from a coarse mesh A to a fine mesh C. The main characteristics of these meshes are listed in Table 5. The structural responses in terms of load P versus d (vertical displacement between the points A and B) are depicted in Fig. 13(a), and those ones obtained in terms of load P versus d (vertical displacement at point B) are depicted in Fig. 13 (b). As can be seen in these figures, the numerical responses are in good agreement with the experimental ones. As in the previous example, no significant mesh dependency is observed. 5.1.4. Notched specimen Fig. 14(a) shows the FE model used to analyze the notched concrete specimen experimentally tested by Kobayashi et al. [59]. The specimen has a thickness of 50:8 mm and two couples of forces F1 (wedge load) and F2 (diagonal compression load) are applied increasingly, as illustrated in Fig. 14(b). For the first loading stage, F1 and F2 are applied simultaneously, while in the second one, F2 is kept constant and F1 is increased. Fig. 12. FE meshes and crack patterns. Table 5 Characteristics of the FE meshes of the four-point bending beam. FE mesh Nnodes NelemðrÞ NelemðiÞ NelemðtÞ mesh A 16,718 7972 15,204 23,176 mesh B 50,615 20,800 48,188 68,988 mesh C 115,772 44,491 112,170 156,661 Fig. 13. Structural responses of the four-point bending beam. Fig. 14. Notched concrete specimen: (a) geometry and loading conditions and (b) loading process. 164 O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 165 The behavior of the concrete (bulk) finite elements is described by a linear elastic material model with Young’s modulus of E ¼ 31 GPa and Poisson’s ratio of m ¼ 0:20. For the interface elements, the tension damage model was adopted with Young’s modulus of E ¼ 31 GPa, Poisson’s ratio of m ¼ 0, fracture energy of Gf ¼ 0:1 N=mm and tensile strength of f t ¼ 2:7 MPa. In this example a mesh dependency is also studied. Fig. 15 shows the three mesh refinements used, whose characteristics are listed in Table 6. A fine fragmented mesh is adopted only in the region where the occurrence of cracks are expected. Comparisons between experimental and numerical results in terms of crack trajectories and applied load F1 versus CMOD (Crack Mouth Opening Displacement) are shown in Fig. 16. Note that these results are in good agreement and do not present significant mesh dependency, as can be seen in Fig. 15. 5.1.5. Double-edge-notched specimen To demonstrate the ability of the technique to describe the propagation of simultaneous arbitrary cracks evolving in mixed mode, the double-edge-notched specimen submitted to a non-proportional loading, experimentally tested by Nooru-Mohamed [60], has been chosen. Fig. 17 shows the test setup discretized with tree-node triangular finite elements in plane stress condition, with an out-of-plane thickness of 50 mm. In this example the mesh of the whole specimen is frag- mented, with exception of the mesh of the rigid steel plate used to impose the loading. Fig. 15. FE meshes and crack trajectories obtained. Table 6 Characteristics of the FE meshes of the notched specimen. FE mesh Nnodes NelemðrÞ NelemðiÞ NelemðtÞ mesh A 2798 1921 2152 4073 mesh B 11,076 5725 9812 15,537 mesh C 43,827 18,671 41,364 60,035 0 Fig. 16. Results of the notched concrete specimen: (a) crack trajectories and (b) curves of F1ðNÞ versus CMOD. 166 O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 The numerical analyses are performed in two loading intervals. First, a lateral force Ps ¼ 10 kN is imposed and then, main- taining this force constant, an increasing vertical displacement Un is imposed to the plate. A linear elastic model is assumed to describe the behavior of the standard triangular finite elements with Young’s mod- ulus of E ¼ 32 GPa and Poisson’s ratio of m ¼ 0:18. The behavior of the interface solid finite elements is described by the ten- sion damage model with Young’s modulus of E ¼ 32 GPa, Poisson’s ratio of m ¼ 0, fracture energy of Gf ¼ 0:025 N=mm and tensile strength of f t ¼ 2:5 MPa. In order to study the mesh dependency, the analyses are performed using different unstructured meshes: coarse (mesh A), medium (mesh B) and very fine mesh (mesh C), as shown in Fig. 18. The main characteristics of these meshes are listed in Table 7. Fig. 19 shows the crack patterns and the corresponding structural responses. As can be observed, the numerical results predict well the experimental one, with the exception of the peak load, whose predictions are higher than that obtained experimentally for all the meshes. As can be seen, the results do not present significant mesh dependency. 5.2. 3D case study 5.2.1. Double-edge-notched specimen In order to verify the capability of the proposed mesh fragmentation technique to simulate crack propagation in 3D anal- ysis, the double-edge-notched specimen presented in the previous example (see Section 5.1.5) are numerically simulated. The same loading process, boundary conditions and material properties of the 2D study are assumed. Fig. 20 shows the three mesh refinements adopted, and Table 8 lists the characteristics of these meshes. Fig. 20 illustrates the crack patterns obtained using the three meshes considered. As can be seen, the results do not pre- sent significant mesh dependency. Moreover, as obtained in the previous example, the structural responses illustrated in steel Fig. 17. Double-edge-notched specimen: geometry and boundary conditions. Fig. 18. FE meshes and crack trajectories obtained. Table 7 Characteristics of the FE meshes of the double-edge notched. FE mesh Nnodes NelemðrÞ NelemðiÞ NelemðtÞ mesh A 2959 1794 2418 4212 mesh B 8470 3632 2280 5912 mesh C 14,826 14,826 20,324 35,150 Fig. 19. Results of the mixed mode test: (a) crack paths and (b) structural responses. Fig. 20. Crack growth trajectories for the different 3D mesh refinements adopted. Table 8 Characteristics of the 3D FE meshes of the double-edge-notched specimen. FE mesh Nnodes NelemðrÞ NelemðiÞ NelemðtÞ mesh D 7160 2537 9678 12,215 mesh E 12,974 4517 18,204 22,721 mesh F 43,398 9633 33,765 43,398 O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 167 Fig. 21. Structural responses of the double-edge-notched specimen: (a) 3D numerical analyses versus experimental response and (b) comparison between 2D and 3D responses. 168 O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 Fig. 21 show that the numerical results in 2D and 3D analyses are in good agreement, and predict well the experimental result, with the exception of the peak load. 6. Concluding remarks In this paper a discrete crack model using special interface elements is presented for modeling cracks in quasi-brittle materials. The proposed mesh fragmentation technique is based on the use of interface solid finite elements with a high aspect ratio as proposed by Manzoli et al. [45]. This strategy is supported by the fact that these elements are able to describe the kinematics of strong discontinuities in the context of the CSDA when their thickness tends to zero. Hence, standard low- order solid elements with a high aspect ratio can be used to describe the behavior of interfaces representing potential crack surfaces using a continuum constitutive model. Therefore, it is not necessary to formulate a new (solid) contact element for this mesh fragmentation technique. A simple continuum damage constitutive model based on the tensile stresses normal to the base of these elements is pro- posed to describe the complete crack formation and propagation process controlled by a prevalent mode I of fracture. It is shown that in the limit situation of strong discontinuity (when the element aspect ratio tends to infinity) the continuum damage model is equivalent to a discrete (cohesive) rigid-damage model, whose softening law can be related to the material fracture energy. Since the damage variable degrades all components of the effective stress, this model is able to deal with cracks evolving in mixed mode. An implicit–explicit integration scheme is used to integrate this constitutive model. With this strategy all the analyses presented in this paper could be performed without problems of convergence, commonly found in problems involving discontinuities. The numerical analyses carried out showed that the proposed strategy is able to predict the propagation of cracks without the need of a special scheme to track the crack paths during the analysis. For all the examples, the numerical results are in good agreement with the experimental ones presented by various researchers. No significant mesh dependency is observed, provided that unstructured meshes with a reasonable refinement are used. The misalignment of the mesh geometry only generates a zigzag crack path very close to the expected (or theoretical) crack surface, without prejudice to the structural responses and crack pattern predictions. Irregularity of the crack surface paths is inherent in this approach and may be incompatible with the theoretical smooth surfaces considered for the evaluation of the fracture energy from experimental tests. This difference may be more pronounced in 3D simulations. Therefore, this aspect should be better studied in order to clarify if the assumed fracture energy parameter of the model must reflect the degree of misalignment of the mesh in order to provide the correct energy dissipation. Finally, with the proposed technique, different crack patterns can be modeled entirely in a continuum framework involv- ing standard solid finite elements and continuum constitutive models compatible with the (regularized) strong discontinuity regime. Therefore, almost only modifications in the pre-processing stage (to introduce finite elements with high aspect ratio in between regular elements of a mesh generated with a standard mesh generator) and in the constitutive model are nec- essary to make existing nonlinear FEM codes able to use the proposed technique to simulate the effects of cracks. Acknowledgments The authors wish to acknowledge the financial support of the National Council for Scientific and Technological Develop- ment (CNPq), Coordination for the Improvement of Higher Education Personnel (CAPES) and Sao Paulo Research Foundation (FAPESP). O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 169 References [1] Rashid Y. Analysis of reinforced concrete pressure vessels. Nucl Engng Des 1968;7:334–44. [2] Ngo D, Scordelis AC. Finite element analysis of reinforced concrete beams. ACI J 1967;64(3):152–63. [3] Borst Rd, Remmers JJC, Needleman A, Abellan M-A. Discrete vs smeared crack models for concrete fracture: bridging the gap. Int J Numer Anal Meth Geomech 2004;28(7–8):583–607. http://dx.doi.org/10.1002/nag.374. [4] Rots J, Blaauwendraad J. Crack models for concrete, discrete or smeared? Fixed, multi-directional or rotating? Heron 34(1). . [5] Tejchman J, Bobiński J. Continuous and discontinuous modelling of fracture in concrete using FEM. Springer series in geomechanics and geoengineering. Berlin: Springer; 2013. http://dx.doi.org/10.1007/978-3-642-28463-2. [6] Gupta A, Akbar H. Cracking in reinforced concrete analysis. J Struct Engng 1984;110(8):1735–46. http://dx.doi.org/10.1061/(ASCE)0733-9445(1984) 110:8(1735). [7] Jirásek M, Zimmermann T. Analysis of rotating crack model. J Engng Mech 1998;124(8):842–51. http://dx.doi.org/10.1061/(ASCE)0733-9399(1998) 124:8(842). [8] Willam K, Pramono E, Sture S. Fundamental issues of smeared crack models. In: Shah S, Swartz S, editors. SEM-RILEM int. conf. on fracture of concrete and rock. Connecticut: Society of Engineering Mechanics; 1987. p. 192–207. [9] Petrangeli M, Ožbolt J. Smeared crack approaches – material modeling. J Engng Mech 1996;122(6):545–54. http://dx.doi.org/10.1061/(ASCE)0733- 9399(1996)122:6(545). [10] Rots JG, Nauta P, Kusters GMA, Blaauwendraadm J. Smeared crack approach an fracture localization in concrete. HERON 1985;30(1). [11] Baza�nt ZP, Oh BH. Crack band theory for fracture of concrete. Matér Constr 1983;93(16):155–77. [12] Carpinteri A, Chiaia B, Nemati KM. Complex fracture energy dissipation in concrete under different loading conditions. Mech Mater 1997;26 (2):93–108. http://dx.doi.org/10.1016/S0167-6636(97)00022-7. . [13] Pijaudier Cabot G, Bazant ZP. Nonlocal damage theory. J Engng Mech, ASCE 1987;113:1512–33. [14] Bažant Z, Jirásek M. Nonlocal integral formulations of plasticity and damage: survey of progress. J Engng Mech 2002;128(11):1119–49. http://dx.doi. org/10.1061/(ASCE)0733-9399(2002)128:11(1119). [15] Peerlings RHJ, de Borst R, Brekelmans WAM, de Vree JHP. Gradient enhanced damage for quasi-brittle materials. Int J Numer Meth Engng 1996;39 (19):3391–403. http://dx.doi.org/10.1002/(SICI)1097-0207(19961015)39:19<3391::AID-NME7>3.0.CO;2-D. [16] Oliver J, Huespe AE. Theoretical and computational issues in modelling material failure in strong discontinuity scenarios. Comput Methods Appl Mech Engng 2004;193(27–29):2987–3014. [17] Oliver J, Linero D, Huespe A, Manzoli O. Two-dimensional modeling of material failure in reinforced concrete by means of a continuum strong discontinuity approach. Comput Methods Appl Mech Engng 2008;197(5):332–48. [18] Manzoli O, Shing P. A general technique to embed non-uniform discontinuities into standard solid finite elements. Comput Struct 2006;84(10– 11):742–57. [19] Simo J, Oliver J, Armero F. An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. Comput Mech 1993;12 (5):277–96. http://dx.doi.org/10.1007/BF00372173. [20] Mariani S, Perego U. Extended finite element method for quasi-brittle fracture. Int J Numer Meth Engng 2003;58(1):103–26. http://dx.doi.org/10.1002/ nme.761. [21] López CM, Carol I, Aguado A. Meso-structural study of concrete fracture using interface elements. I: Numerical model and tensile behavior. Mater Struct 2008;41(3):583–99. http://dx.doi.org/10.1617/s11527-007-9314-1. [22] López CM, Carol I, Aguado A. Meso-structural study of concrete fracture using interface elements. II: Compression, biaxial and brazilian test. Mater Struct 2008;41(3):601–20. http://dx.doi.org/10.1617/s11527-007-9312-3. [23] Carol I, Prat PC, López CM. Normal/shear cracking model: application to discrete crack analysis. J Engng Mech 1997;123(8):765–73. http://dx.doi.org/ 10.1061/(ASCE)0733-9399(1997)123:8(765). [24] Needleman A. An analysis of decohesion along an imperfect interface. Int J Fract 1990;42(1):21–40. http://dx.doi.org/10.1007/BF00018611. [25] Ortiz M, Pandolfi A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. Int J Numer Meth Engng 1999;44(9):1267–82. http://dx.doi.org/10.1002/(SICI)1097-0207(19990330)44:9<1267::AID-NME486>3.0.CO;2-7. [26] Xu X-P, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 1994;42(9):1397–434. http://dx.doi.org/10.1016/ 0022-5096(94)90003-5. . [27] Radulovic R, Bruhns OT, Mosler J. Effective 3D failure simulations by combining the advantages of embedded strong discontinuity approaches and classical interface elements. Engng Fract Mech 2011;78(12):2470–85. http://dx.doi.org/10.1016/j.engfracmech.2011.06.007. . [28] Asferg JL. Modeling of concrete fracture applying the extended finite element method. Ph.D. Thesis, Technical University of Denmark; 2006. [29] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Engng Fract Mech 2002;69(7):813–33. http://dx.doi.org/10.1016/ S0013-7944(01)00128-X. . [30] Gupta P, Pereira J, Kim D-J, Duarte C, Eason T. Analysis of three-dimensional fracture mechanics problems: a non-intrusive approach using a generalized finite element method. Engng Fract Mech 2012;90:41–64. http://dx.doi.org/10.1016/j.engfracmech.2012.04.014. . [31] Feist C, Hofstetter G. Three-dimensional fracture simulations based on the SDA. Int J Numer Anal Meth Geomech 2007;31(2):189–212. http://dx.doi. org/10.1002/nag.542. [32] Wells G, Sluys L. Three-dimensional embedded discontinuity model for brittle fracture. Int J Solids Struct 2001;38(5):897–913. http://dx.doi.org/ 10.1016/S0020-7683(00)00029-9. . [33] Oliver J, Huespe A, Pulido M, Chaves E. From continuum mechanics to fracture mechanics: the strong discontinuity approach. Engng Fract Mech 2002;69(2):113–36. http://dx.doi.org/10.1016/S0013-7944(01)00060-1. . [34] Armero F, Garikipati K. An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids. Int J Solids Struct 1996;33(20-22):2863–85. http://dx.doi.org/10.1016/0020-7683(95)00257-X. . [35] Oliver J, Huespe A, Sánchez P. A comparative study on finite elements for capturing strong discontinuities: E-FEM vs X-FEM. Comput Methods Appl Mech Engng 2006;195(37–40):4732–52. http://dx.doi.org/10.1016/j.cma.2005.09.020 [John H. Argyris Memorial Issue. Part I.] . [36] Jäager P, Steinmann P, Kuhl E. Modeling three-dimensional crack propagation – a comparison of crack path tracking strategies. Int J Numer Meth Engng 2008;76(9):1328–52. http://dx.doi.org/10.1002/nme.2353. [37] Oliver J, Huespe AE, Samaniego E, Chaves EWV. Continuum approach to the numerical simulation of material failure in concrete. Int J Numer Anal Meth Geomech 2004;28:609–32. [38] Manzoli O, Claro G, Rodrigues E, Lopes J. A local-global scheme for tracking crack path in three-dimensional solids. Comput Concr 2013;12(3):261–83. http://dx.doi.org/10.12989/cac.2013.12.3.261. [39] Osher S, Sethian J. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton–Jacobi formulations. J Comput Phys 1988;79 (1):12–49. http://dx.doi.org/10.1016/0021-9991(88)90002-2. . [40] Pandolfi A, Ortiz M. An efficient adaptive procedure for three-dimensional fragmentation simulations. Engng Comput 2002;18(2):148–59. http://dx. doi.org/10.1007/s003660200013. http://refhub.elsevier.com/S0013-7944(15)00702-X/h0005 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0010 http://dx.doi.org/10.1002/nag.374 http://repository.tudelft.nl/view/ir/uuid:0a401939-2938-4f9d-a395-b6a6652b2cd9/ http://repository.tudelft.nl/view/ir/uuid:0a401939-2938-4f9d-a395-b6a6652b2cd9/ http://dx.doi.org/10.1007/978-3-642-28463-2 http://dx.doi.org/10.1061/(ASCE)0733-9445(1984)110:8(1735) http://dx.doi.org/10.1061/(ASCE)0733-9445(1984)110:8(1735) http://dx.doi.org/10.1061/(ASCE)0733-9399(1998)124:8(842) http://dx.doi.org/10.1061/(ASCE)0733-9399(1998)124:8(842) http://refhub.elsevier.com/S0013-7944(15)00702-X/h0040 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0040 http://dx.doi.org/10.1061/(ASCE)0733-9399(1996)122:6(545) http://dx.doi.org/10.1061/(ASCE)0733-9399(1996)122:6(545) http://refhub.elsevier.com/S0013-7944(15)00702-X/h9000 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0055 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0055 http://dx.doi.org/10.1016/S0167-6636(97)00022-7 http://www.sciencedirect.com/science/article/pii/S0167663697000227 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0065 http://dx.doi.org/10.1061/(ASCE)0733-9399(2002)128:11(1119) http://dx.doi.org/10.1061/(ASCE)0733-9399(2002)128:11(1119) http://dx.doi.org/10.1002/(SICI)1097-0207(19961015)39:19<3391::AID-NME7>3.0.CO;2-D http://refhub.elsevier.com/S0013-7944(15)00702-X/h0080 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0080 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0085 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0085 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0090 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0090 http://dx.doi.org/10.1007/BF00372173 http://dx.doi.org/10.1002/nme.761 http://dx.doi.org/10.1002/nme.761 http://dx.doi.org/10.1617/s11527-007-9314-1 http://dx.doi.org/10.1617/s11527-007-9312-3 http://dx.doi.org/10.1061/(ASCE)0733-9399(1997)123:8(765) http://dx.doi.org/10.1061/(ASCE)0733-9399(1997)123:8(765) http://dx.doi.org/10.1007/BF00018611 http://dx.doi.org/10.1002/(SICI)1097-0207(19990330)44:9<1267::AID-NME486>3.0.CO;2-7 http://dx.doi.org/10.1016/0022-5096(94)90003-5 http://dx.doi.org/10.1016/0022-5096(94)90003-5 http://www.sciencedirect.com/science/article/pii/0022509694900035 http://dx.doi.org/10.1016/j.engfracmech.2011.06.007 http://www.sciencedirect.com/science/article/pii/S0013794411002311 http://www.sciencedirect.com/science/article/pii/S0013794411002311 http://dx.doi.org/10.1016/S0013-7944(01)00128-X http://dx.doi.org/10.1016/S0013-7944(01)00128-X http://www.sciencedirect.com/science/article/pii/S001379440100128X http://dx.doi.org/10.1016/j.engfracmech.2012.04.014 http://www.sciencedirect.com/science/article/pii/S0013794412001476 http://www.sciencedirect.com/science/article/pii/S0013794412001476 http://dx.doi.org/10.1002/nag.542 http://dx.doi.org/10.1002/nag.542 http://dx.doi.org/10.1016/S0020-7683(00)00029-9 http://dx.doi.org/10.1016/S0020-7683(00)00029-9 http://www.sciencedirect.com/science/article/pii/S0020768300000299 http://dx.doi.org/10.1016/S0013-7944(01)00060-1 http://www.sciencedirect.com/science/article/pii/S0013794401000601 http://dx.doi.org/10.1016/0020-7683(95)00257-X http://www.sciencedirect.com/science/article/pii/002076839500257X http://www.sciencedirect.com/science/article/pii/002076839500257X http://dx.doi.org/10.1016/j.cma.2005.09.020 http://www.sciencedirect.com/science/article/pii/S0045782505005049 http://www.sciencedirect.com/science/article/pii/S0045782505005049 http://dx.doi.org/10.1002/nme.2353 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0185 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0185 http://dx.doi.org/10.12989/cac.2013.12.3.261 http://dx.doi.org/10.1016/0021-9991(88)90002-2 http://www.sciencedirect.com/science/article/pii/0021999188900022 http://dx.doi.org/10.1007/s003660200013 http://dx.doi.org/10.1007/s003660200013 170 O.L. Manzoli et al. / Engineering Fracture Mechanics 153 (2016) 151–170 [41] Molinari JF, Gazonas G, Raghupathy R, Rusinek A, Zhou F. The cohesive element approach to dynamic fragmentation: the question of energy convergence. Int J Numer Meth Engng 2007;69(3):484–503. http://dx.doi.org/10.1002/nme.1777. [42] Caballero A, Willam K, Carol I. Consistent tangent formulation for 3D interface modeling of cracking/fracture in quasi-brittle materials. Comput Methods Appl Mech Engng 2008;197(33–40):2804–22. http://dx.doi.org/10.1016/j.cma.2008.01.011. . [43] Etse G, Caggiano A, Vrech S. Multiscale failure analysis of fiber reinforced concrete based on a discrete crack model. Int J Fract 2012;178(1–2):131–46. http://dx.doi.org/10.1007/s10704-012-9733-z. [44] Caggiano A, Etse G, Martinelli E. Zero-thickness interface model formulation for failure behavior of fiber-reinforced cementitious composites. Comput Struct 2012;98–99:23–32. http://dx.doi.org/10.1016/j.compstruc.2012.01.013. . [45] Manzoli O, Gamino A, Rodrigues E, Claro G. Modeling of interfaces in two-dimensional problems using solid finite elements with high aspect ratio. Comput Struct 2012;94–95(0):70–82. http://dx.doi.org/10.1016/j.compstruc.2011.12.001. . [46] Oliver J, Cervera M, Manzoli OL. Strong discontinuities and continuum plasticity models: the strong discontinuity approach. Int J Plast 1999;15 (3):319–51. [47] Oliver J. On the discrete constitutive models induced by strong discontinuity kinematics and continuum constitutive equations. Int J Solids Struct 2000;37(48–50):7207–29. http://dx.doi.org/10.1016/S0020-7683(00)00196-7. . [48] Sánchez M, Manzoli O, Guimarães L. Modeling 3-D desiccation soil crack networks using a mesh fragmentation technique. Comput Geotech 2014;62 (0):27–39. http://dx.doi.org/10.1016/j.compgeo.2014.06.009. . [49] Rodrigues E, Manzoli O, Bitencourt Jr L, Prazeres P, Bittencourt T. Failure behavior modeling of slender reinforced concrete columns subjected to eccentric load. Latin Am J Solids Struct 2015;12(3):520–41. http://dx.doi.org/10.1590/1679-78251224. [50] Oliver J, Huespe A, Blanco S, Linero D. Stability and robustness issues in numerical modeling of material failure with the strong discontinuity approach. Comput Methods Appl Mech Engng 2006;195(52):7093–114 [computational Modelling of Concrete]. [51] Oliver J, Huespe A, Cante J. An implicit/explicit integration scheme to increase computability of non-linear material and contact/friction problems. Comput Methods Appl Mech Engng 2008;197(21–24):1865–89. http://dx.doi.org/10.1016/j.cma.2007.11.027. . [52] Kachanov LM. Introduction to continuum damage mechanics. Netherlands: Springer; 1986. [53] Lemaitre J. A course on damage mechanics. Springer; 1992. [54] Murakami S. Continuum damage mechanics – a continuum mechanics approach to the analysis of damage and fracture. Springer; 2012. [55] Cervera M, Oliver J, Manzoli O. A rate-dependent isotropic damage model for the seismic analysis of concrete dams. Earthq Engng Struct Dyn 1996;25 (9):987–1010. http://dx.doi.org/10.1002/(SICI)1096-9845(199609)25:9<987::AID-EQE599>3.0.CO;2-X. [56] Prazeres PG, Bitencourt Jr LA, Bittencourt TN, Manzoli OL. A modified implicit–explicit integration scheme: an application to elastoplasticity problems. J Braz Soc Mech Sci Engng 2016;38(1):151–61. http://dx.doi.org/10.1007/s40430-015-0343-3. [57] Le Bellégo C, Gérard B, Pijaudier-Cabot G. Chemo-mechanical effects in mortar beams subjected to water hydrolysis. J Engng Mech 2000;126 (3):266–72. http://dx.doi.org/10.1061/(ASCE)0733-9399(2000)126:3(266). [58] Arrea M, Ingraffea AR. Mixed-mode crack propagation in mortar and concrete. Technical report, 81-13, Deptartment of Structural Engineering, Cornell University, New York. [59] Kobayashi A, Hawkins N, Barker D, Liaw B. Fracture process zone of concrete. In: Shah S, editor. Application of fracture mechanics to cementitious composites. NATO ASI series, vol. 94. Netherlands: Springer; 1985. p. 25–50. [60] Nooru-Mohamed MB. Mixed-mode fracture of concrete: an experimental approach. Ph.D. Thesis, Delft University of Technology, Delft; 1992. http://dx.doi.org/10.1002/nme.1777 http://dx.doi.org/10.1016/j.cma.2008.01.011 http://www.sciencedirect.com/science/article/pii/S0045782508000352 http://www.sciencedirect.com/science/article/pii/S0045782508000352 http://dx.doi.org/10.1007/s10704-012-9733-z http://dx.doi.org/10.1016/j.compstruc.2012.01.013 http://www.sciencedirect.com/science/article/pii/S0045794912000284 http://www.sciencedirect.com/science/article/pii/S0045794912000284 http://dx.doi.org/10.1016/j.compstruc.2011.12.001 http://www.sciencedirect.com/science/article/pii/S0045794911002938 http://www.sciencedirect.com/science/article/pii/S0045794911002938 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0230 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0230 http://dx.doi.org/10.1016/S0020-7683(00)00196-7 http://www.sciencedirect.com/science/article/pii/S0020768300001967 http://www.sciencedirect.com/science/article/pii/S0020768300001967 http://dx.doi.org/10.1016/j.compgeo.2014.06.009 http://www.sciencedirect.com/science/article/pii/S0266352X14001177 http://dx.doi.org/10.1590/1679-78251224 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0250 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0250 http://dx.doi.org/10.1016/j.cma.2007.11.027 http://www.sciencedirect.com/science/article/pii/S0045782507004756 http://www.sciencedirect.com/science/article/pii/S0045782507004756 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0260 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0265 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0270 http://dx.doi.org/10.1002/(SICI)1096-9845(199609)25:9<987::AID-EQE599>3.0.CO;2-X http://dx.doi.org/10.1007/s40430-015-0343-3 http://dx.doi.org/10.1061/(ASCE)0733-9399(2000)126:3(266) http://refhub.elsevier.com/S0013-7944(15)00702-X/h0295 http://refhub.elsevier.com/S0013-7944(15)00702-X/h0295 On the use of finite elements with a high aspect ratio �for modeling cracks in quasi-brittle materials 1 Introduction 2 Overview of the proposed mesh fragmentation technique 3 Interface solid finite element 3.1 Decomposition of the strain tensor 3.1.1 Components of the strain tensor for the 3-node triangular finite element 3.1.2 Components of the strain tensor for the 4-node tetrahedral finite element 4 Tension damage model 4.1 Continuum model 4.1.1 Implicit–explicit integration scheme for the tension damage model 4.2 Discrete relation of the interface solid finite element 5 Numerical examples 5.1 2D case studies 5.1.1 Plate with hole 5.1.2 Three-point bending beams 5.1.3 Four-point bending beam 5.1.4 Notched specimen 5.1.5 Double-edge-notched specimen 5.2 3D case study 5.2.1 Double-edge-notched specimen 6 Concluding remarks Acknowledgments References