Modelling crack growth in quasi-brittle materials using mesh fragmentation technique

Modelling crack growth in quasi-brittle materials using mesh fragmentation technique Maedo, Michael1; Manzoli, Osvaldo2; Bitencourt Jr., Luís3; Rodrigues, Eduardo4 ABSTRACT The modelling of discontinuities with new finite element techniques (X-FEM and E-FEM) is cumbersome due to the need of algorithm tracking to simulate the crack growth. Therefore, a new technique to capture strong discontinuities in quasi-brittle materials without the need of a crack tracking scheme is presented. The method consists of inserting solid finite elements with high aspect ratio (ratio of the largest to the smallest dimension) between the bulk elements of the original mesh. Hence, when the ultimate tensile strength is reached only these interface elements are activated. Since the analyses are performed integrally in the context of the continuum mechanics, there is no need to use tracking algorithms to simulate the formation of cracks. The study is performed using different specimens in order to cover the cases that can show up in practice. Numerical results show that the technique is able to represent crack growth. Keywords: solid finite elements with high aspect ratio, tensile damage model, quasi-brittle materials, crack growth, strong discontinuities. 1. INTRODUCTION The formation and propagation of cracks involves many engineering fields. In concrete, the crack process rules the nonlinear behaviour of the material, so that significant cracks can lead the material to a structural collapse. On the other hand, in soil mechanics the formation and propagation of cracks may affect the soil, decreasing its strength and sometimes causing erosions. In recent years, there are many techniques to represent crack growth in quasi-brittle materials. In the context of the finite element method, some methods such as E-FEM[1,2,3] and X-FEM[4,5] require tracking algorithms to represent crack growth. These tracking techniques are simple to represent a few cracks in 2D analyses, but are very complex for 3D problems. Thus, in order to describe the crack process without the need of crack tracking schemes, a new technique that consists of introducing in special interface elements between the regular elements of the mesh is proposed [6,7]. Thereby, these special elements must be able to represent the discontinuities, and consequently, the cracks develop along the boundaries of the regular elements. 1Department of Mechanical Engineering. UNESP – Sao Paulo State University (BRAZIL). [email protected] (Corresponding author) 2Department of Civil Engineering. UNESP – Sao Paulo State University (BRAZIL). oman[email protected]esp.br 3Polytechnic School. USP – University of Sao Paulo (BRAZIL). luis.bitenco[email protected] 4Polytechnic School. USP – University of Sao Paulo (BRAZIL). eduardoale[email protected] 373 Modeling crack growth in quasi-brittle materials using mesh fragmentation technique Third International Conference on Mechanical Models in Structural Engineering University of Seville. 24-26 june 2015. However, instead of zero-thickness interface elements this work presents a new technique that uses solid finite elements with high aspect ratio (ratio of the larges to the smallest dimension) to represent discontinuities. Throughout this paper it is shown that, as the aspect ratio increases, the element strain also increase, approaching the kinematics associated to the discontinuity, as in the case of Continuum Strong Discontinuity Approach (CSDA) [8,9]. Thus, based on the same principles of the CSDA, bounded stresses are obtained from unbounded strains. To model the nonlinear behaviour of the interfaces, a tensile damage constitutive model that relates stress to strain is described. Since the analyses are performed integrally in the context of the continuum mechanics, complex crack patterns and 3D cracks can be properly addressed using the technique proposed. 2. INTERFACE SOLID FINITE ELEMENT Let us consider a three-node triangular finite element, with base  given by the distance between nodes 2 and 3, and height ℎ formed by the segment between node 1 and its projection on the base of the element, as shown in Fig. 1. Figure 1. Interface triangular element in the coordinate system ,[10] Assuming the coordinate system ,, where  is the unit normal to the base of the element, the strain tensor for any point can be written as follows [10]:  =   +   (1) with = 1       0 1 2       −       0 1 2     −         −     0 0 0 0      (2) and = 1 ℎ       !  1 2  !  0 1 2  !  0 0 0 0 0      (3) 374 Maedo1, Manzoli2, Bitencourt Jr.3 and Rodrigues4 where  " and " are the components of the displacement of node # in the  and  directions, respectively; and != $− $% and !=$−$% are the components of the relative displacement between node 1 and its projection on the base (1’). The part of the strain tensor that depends on the height ℎ, which is expressed b Eq. (3), can be written in a more general form as:   = 1 ℎ   ⊗ ' !  ( (4) where ∙( is the symmetric part of ∙, '! is relative displacement between node 1 and the projection (1’) and ⊗ is the dyadic product. Thus, substituting Eq. (4) into (1), we obtain:  =   * +,-./0/ + 1 ℎ   ⊗ ' !  ( 1 2 2 2 3 2 2 2 4 -.+,-./0/ (5) Taking the limit as ℎ→0, the strain component  remains bounded while  becomes unbounded. In this case, the deformations of the element are almost exclusively defined by the displacements between node 1 and its projection on the base of the element. Note that, when ℎ tends to zero, node 1 and the projection (1’) converge to the same material point. Therefore, '! becomes the measure of the strong discontinuity. The strain field given by Eq. (5) consists of the kinematics of the CSDA, and consequently, the same applications of the CSDA can be treated with finite elements with high aspect ratio [6]. 3. TENSILE DAMAGE MODEL The tensile damage model is designed to describe crack growth in quasi-brittle materials. The damage criterion is based on the component of the tensile stress normal to the base of the element with high aspect ratio, so that the model is able to represent the crack process in mode I. It is worth mentioning that this model degrades all components of the effective stress, and therefore the crack initiation can evolve in mixed mode. 3.1. Continuum constitutive model The continuum damage model proposed is defined by the following constitutive relation: 6 =  1 − 7  6 8  (6 ) where 6 is the stress tensor, 7∈:0,1; is the damage variable, and 6 8 is the effective stress tensor expressed as: 6 8 = < ∶ >  (7 ) where <=2?@+A+2/3?D⊗D is the fourth order elastic tensor, ? is the shear modulus, A is the bulk modulus, @ is the fourth order unit tensor, and D is the second order unit tensor. 375 Modeling crack growth in quasi-brittle materials using mesh fragmentation technique Third International Conference on Mechanical Models in Structural Engineering University of Seville. 24-26 june 2015. The damage criterion is based on the component of the stress normal to the base of the element, in other words: E = F  − G  H  ≤ 0  (8 ) Where G and H are the stress and strain-like internal variables, respectively Dividing Eq. (8) by 1−7, the damage criterion can be rewritten as: E J = F J  − H ≤ 0  (9 ) Since H=G/1−7, the damage evolution rule is given by 7  H  = 1 − G  H  H (10) The constitutive model is completed by the Kuhn-Tucker relations E J ≤ 0 ;  H L ≥ 0 ; H L  E J = 0 (1 1 ) and the consistency condition H L  E J L = 0  if  E J = 0 (12) Based on Eqs. (9), (11), and (12), the evolution law for the strain-like internal variable can be expressed as: H  P  = max ( ∈ : T , U ; : V W    , X U ; (13) Thus, H assumes the maximum value that FJ reaches throughout the loading process, starting from the initial tensile strength of the material XU. 3.2. Discrete constitutive relation In the limit case of strong discontinuity, the continuum constitutive relation given by Eq. (6) becomes a discrete constitutive relation that relates the interface stresses to the components of the displacement jump. Therefore, Y = lim [ → T  \ ⋅  = lim [ → T   1 − 7   ⋅ < ∶    +    = lim [ → T   1 − 7   ⋅ < ∶ ^   + 1 ℎ   ⊗ ' !  ( _ =lim [→T1−7 ℎ⋅<∶:ℎ+⊗ '!(; = lim [ → T ` 1−7 ℎ    ⋅ < ⋅   a 1 2 2 2 2 2 3 2 2 2 2 2 4 b c ⋅ ' ! = lim [ → T b c ⋅ ' ! (14) 376 Maedo1, Manzoli2, Bitencourt Jr.3 and Rodrigues4 where Y=dFFFUef, '!=d ! ! !Uef and g c =  1 − 7  ℎ h 4 / 3 ? + A 0 0 0 ? 0 0 0 ? j (15) According to CSDA, the stresses remain bounded even for unbounded strains. In the case of the element with high aspect ratio, to maintain the stresses bounded when ℎ→0, the components of the displacement jump must tend to zero in the elastic regime 7=0. On the other hand, in the inelastic regime  '!≠0, the damage variable 7 must tends to 1 to hold the stresses bounded. Assuming Poisson’s ratio null l=0, we obtain: 4 ? 3 ⁄ + A = n  and  ? = n / 2 (1 6 ) and therefore, the discrete relation can be written in terms of the Young’s modulus. Table 1 shows the continuum and discrete constitutive expressions for the tensile damage model. Table 1. Continuum and discrete tensile damage models Continuum model Discrete model Constitutive relation 6 =  1 − 7  6 8 6 8 = < : > Y = b c ⋅ r ! Equivalent stress V W = F J  V W = F J  = 1 ℎ n  !  Damage criterion E J = V W − H ≤ 0 E J = V W − H ≤ 0 Evolution law for the strain-type variable H = max  ∈ : T , U ; : V W    , X U ; H = max  ∈ : T , U ; : V W    , X U ; Damage evolution 7  H  = 1 − G  H  H 7  H  = 1 − G  H  H Hardening law G  H  G  H  The evolution of the normal stress for a normal displacement monotonically increasing in mode I (i.e. L!≥0, !│UtT=0, != !U=0) is: F    !   =  1 − 7  ℎ n  !  = u 1 ℎ v n  !            if   ! ≤  T !  G  1 ℎ v n  !      if   ! >  T !  (17) with T!=ℎXU/n. Let us assume that the hardening law is governed by the following exponential equation: G  H  = X U x y  [  $ z { / | }  (18) where XU is the tensile strength of the material and ~ is the softening parameter. Thus, the energy dissipated in a complete degradation process in mode one of the interface element is 377 Modeling crack growth in quasi-brittle materials using mesh fragmentation technique Third International Conference on Mechanical Models in Structural Engineering University of Seville. 24-26 june 2015. ? | = • F    !   € T d  !  = • F    !   • ‚ ! ƒ T d  !  + • F    !   € •‚!ƒ d  !  = • 1 ℎ v n  !  •‚!ƒ T d  !  + • G  1 ℎ v n  !   € • ‚ ! ƒ d  !  =  X U  ℎ 2 ⁄ + X U  ~ ⁄  n (19) Taking the limit as ℎ→0, the fracture energy becomes: ? | = X U  ~ n (20) From Eq. (20), we can define the softening parameter in terms of the material properties, i.e. ~ = X U  ? | n (21) 4. NUMERICAL EXAMPLES The 2D numerical analyses were performed using three-node triangular finite element, while in 3D cases, tetrahedron finite elements were used in simulations. To describe crack growth, interfaces elements are inserted between the bulk (regular) elements of the mesh. Although in some cases the entire mesh is fragmented, in general only the part of the mesh corresponding to the region of study is fragmented, reducing computational effort and consequently time process. As previously mentioned, the interfaces are represented using elements with high aspect ratio, which are modelled using the tensile damage model equipped with exponential softening. 4.1. Gravity dam In order to show the capability of the technique to predict cracks in quasi-brittle materials, the gravity dam, tested experimentally by Carpinteri et al. [11], is analysed by assuming the geometry and boundary conditions shown in Fig. 2. The thickness of the specimen is 300 mm and the depth of the notch of the dam is 150 mm. For the numerical analyses, the behaviour of the bulk elements is governed by the linear elastic model with n=35.7 GPa, and l=0.1. On the other hand, for the elements with high aspect ratio the following material properties are assumed: n=35.7 GPa, l=0, ?|=184 N/m, and XU=3.6 MPa. 378 Maedo1, Manzoli2, Bitencourt Jr.3 and Rodrigues4 Figure 2. Gravity dam: geometry and boundary conditions. The analyses were performed using three unstructured meshes with different refinements. Fig. 3 shows the crack pattern at the end of the analysis for all meshes studied. As can be seen from Fig. 4a, the structural responses, in terms of the load vs. CMOD (Crack Mouth Opening Displacement), derived from the numerical results are compared with the experimental one. Moreover, the numerical and experimental crack trajectories were contrasted as illustrated in Fig. 4b. mesh A mesh B mesh C Figure 3. Gravity dam: deformed meshes. 80 60 40 20 100 (a) (b) Figure 4. Gravity dam: (a) structural responses; (b) crack paths. 379 Modeling crack growth in quasi-brittle materials using mesh fragmentation technique Third International Conference on Mechanical Models in Structural Engineering University of Seville. 24-26 june 2015. 4.2. Four point bending test In order to study the technique with respect to the mesh dependency, the notched beam tested by Arrea and Ingraffea [12] is simulated. The geometry and the boundary conditions of the problem are sketched in Fig. 5, where the thickness is 150 mm. Figure 5. Four point bending test: geometry and boundary conditions. The considered material properties of the bulk elements are n =28.8GPa and l=0.18. For the interface elements, the following material parameters have been adopted: n=28.8GPa, l=0,?|=100N/m, andXU=2.8MPa. Six meshes, three unstructured and three structured, were used to study mesh dependency. Fig. 6 shows the crack trajectory for the six meshes, where the A, B and C meshes are the unstructured ones, and the A’, B’ and C’ meshes are structured. mesh A mesh B mesh C mesh A’ mesh B’ mesh C’ Figure 6. Four point bending test: deformed meshes. 380 Maedo1, Manzoli2, Bitencourt Jr.3 and Rodrigues4 Fig. 7a-c illustrates the numerical and experimental results in terms of load vs. • (which is the crack mouth sliding displacement). Thus, in Fig. 7a the curves of the unstructured meshes are shown to fit very well in the experimental region. On the other hand, the structural response of the structured meshes deviates from the experimental band in some cases. Fig 7c compares all obtained results. Fig. 7d-f presents the load-deflection at the point where P is applied. Thereby, Fig. 7d and 7e show the load-deflection curve of the unstructured and structured meshes, respectively. Analogously to Fig. 7c, the response of the six meshes is shown in Fig. 7f. experimental mesh A mesh B mesh C experimental mesh A’ mesh B’ mesh C’ experimental mesh A mesh B mesh C mesh A’ mesh B’ mesh C’ (a) (b) (c) mesh A mesh B mesh C mesh A’ mesh B’ mesh C’ mesh A mesh B mesh C mesh A’ mesh B’ mesh C’ (d) (e) (f) Figure 7. Four point bending test: (a) load vs. CMSD curves of the unstructured meshes; (b) load vs. CMSD curves of the structured meshes; (c) load vs. CMSD curves of the unstructured and structured meshes; (d) load-deflection curves of the unstructured meshes; (e) load-deflection curves of the structured meshes; (f) load-deflection curves of the unstructured and structured meshes. 4.3. Three point bending test This test is performed to show that the technique is also able to predict crack growth by assuming a three-dimensional analysis. Hence, the three point bending test reported by Béllego et al. [13] was studied by considering both 2D and 3D analyses. The dimensions and the boundary conditions of the problem are indicated in Fig. 8 and the thickness of the specimen is 40 mm. Figure 8. Three point bending test: geometry and boundary conditions. 381