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

Abstract

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 elemental 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 experimental data and representative numerical simulations than the linear model.

Keywords:
fracture; quasi-brittle material; elemental enrichment finite element method (E-FEM)

Resumo

Neste artigo, um modelo numérico com descontinuidades fortes é apresentado para abordar problemas de fratura em materiais quase frágeis. É implementada uma formulação não simétrica estaticamente e cinemáticamente consistente. A descontinuidade forte no campo de deslocamentos é representada usando o método do elemento finito com enriquecimento elementar (E-FEM). Em outras palavras, a descontinuidade forte é introduzida no elemento finito e os graus de liberdade adicionais são condensados em nível de elemento, permitindo a implementação em códigos computacionais existentes. Dois modelos constitutivos são utilizados para analisar o comportamento da região fissurada, um linear e outro exponencial. Os resultados do modelo exponencial estão mais próximos dos obtidos em dados experimentais e simulações numéricas representativas do que o modelo linear.

Palavras-chave:
fratura; material quase-frágil; método de elementos finitos com enriquecimento elementar (E-FEM)

1. 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 [¹].

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 [²], [³], [⁴], and [⁵]) and smeared cracks (see [⁶], [⁷], [⁸] and [⁹]). 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 ([¹⁰], [¹¹]) and the elemental enrichment (E-FEM), involving finite elements with internal degrees of freedom to represent strong or weak discontinuities ([¹²], [⁸], [¹³], [¹⁴]).

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 [¹⁵].

This paper presents a strong discontinuity model proposed by [¹²]. The model is nonsymmetrical (SKON) according to the nomenclature of [¹⁵] 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 [¹⁶]).

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.

2. 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 [⁸].

In the implemented model [¹²], 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 $\widehat{b}$ are applied. The surface is divided into two parts: S_u, where the essential boundary conditions, $u=\widehat{u}$, are applied; and S_t, in which the natural boundary conditions are applied (Figure 1a) [¹⁵].

Figure 1
Boundary condition in the domain. a) continuous, b) with inner interface

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 [⁸].

The field equations that govern the problem can be coupled on a variational principle according to equation (1):

$\underset{V}{{\displaystyle \int }}\delta {\epsilon }^T\sigma \left(\epsilon \right)dV+\delta \underset{V}{{\displaystyle \int }}{\sigma }^T\left(\nabla u-\epsilon \right)dV=\underset{S_t}{{\displaystyle \int }}\delta u^T\widehat{t}dS+\underset{V}{{\displaystyle \int }}\delta u^T\widehat{b}dV$ (1)

in which δ represents variation, σ(ε) the stress obtained from the constitutive relations, ?u the displacement gradient.

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 [¹⁷].

$\nabla u=\epsilon$ (2)

$\sigma \left(\epsilon \right)=\sigma$ (3)

$\nabla \sigma .n=\widehat{b}$ (4)

$n\sigma =\widehat{t}$ (5)

where n is the outward normal vector to the boundary.

3. 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 [⁸], [¹⁵] and [¹⁸].

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 [¹⁹]. Therefore, the field of displacements with discontinuities is:

$u\approx N{d}_N+N_C{d}_C$ (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:

$\epsilon \approx B{d}_N+Ge$ (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 [⁸].

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:

$\sigma \approx Ss$ (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_Cd_C ) = B_Cd_C, we have:

$\delta d_{\text{N}}^T\underset{V}{{\displaystyle \int }}{B}^T\sigma \left(\epsilon \right)dV+\delta {\epsilon }^T\underset{V}{{\displaystyle \int }}{G}^T\left[\sigma \left(\epsilon \right)-Ss\right]dV+\delta s^T\underset{V}{{\displaystyle \int }}{S}^T\left(B_C{d}_C-Ge\right)dV+\delta d_{\text{C}}^T\underset{V}{{\displaystyle \int }}{B}_{\text{C}}^{T}SsdV=\delta d_{\text{N}}^T{f}_{EXT}+\delta d_C^T{f}_C$ (9)

In that, f_EXT are the usual external forces and f_c are the additional external forcers defined by:

$f_{EXT}=\underset{V}{{\displaystyle \int }}{N}^T\widehat{b}dV+\underset{S_t}{{\displaystyle \int }}{N}^T\widehat{t}dS$ (10)

and

$f_c=\underset{V}{{\displaystyle \int }}{N}_{\text{C}}^T\widehat{b}dV+\underset{S_t}{{\displaystyle \int }}{N}_{\text{C}}^T\widehat{t}dS$ (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:

$\underset{V}{{\displaystyle \int }}{B}^T\sigma \left(\epsilon \right)dV=f_{EXT}$ (12)

$\underset{V}{{\displaystyle \int }}{G}^T\sigma \left(\epsilon \right)dV-\left(\underset{V}{{\displaystyle \int }}{G}^{T}SdV\right)s=0$ (13)

$\left(\underset{V}{{\displaystyle \int }}{S}^T{B}_{C}dV\right)d_C-\left(\underset{V}{{\displaystyle \int }}{S}^{T}GdV\right)e=0$ (14)

$\left(\underset{V}{{\displaystyle \int }}{B}_C^{T}SdV\right)s=0$ (15)

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:

$\sigma =D\epsilon \approx D\left(B{d}_N+Ge\right)$ (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:

$\underset{V}{{\displaystyle \int }}\left[\begin{array}{c}\begin{array}{cc}{B}^{T}DB& B^{T}DG\\ G^{T}DB& G^{T}DG\end{array}\begin{array}{cc}0& 0\\ -G^{T}S& 0\end{array}\\ \begin{array}{cc}0& -S^{T}G\\ 0& 0\end{array}\begin{array}{cc}0& S^T{B}_C\\ B_C^{T}S& 0\end{array}\end{array}\right]dV\left\{\begin{array}{c}{d}_N\\ e\\ s\\ d_C\end{array}\right\}=\left\{\begin{array}{c}{f}_{EXT}\\ 0\\ 0\\ \end{array}\right\}$ (17)

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 [⁸].

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ⁱ _EXT .

According to [¹⁵], 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 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 [¹⁹], for more information see [¹⁶].

4. Asymmetric model implemented

In this paper the asymmetric model proposed by [¹²] is implemented to represent the strong discontinuities, with the following characteristics:

• Consider the entire element as a minimum domain in the localization of strains, instead of working at the integration point level.

• The strain localization within the finite element is considered as a displacement discontinuity line incorporated in the element domain.

• Two constitutive relations are defined to represent the material behavior when the localization is started at the element. A stress-displacement relation for the discontinuity line, related to fracture energy, and a stress-strain relation for the element domain.

• The elements resulting from this formulation are non-conforming.

• 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).

Figure 2
Finite element with localization line a) one-dimensional, b) two-dimensional

Before the opening of cracks, displacements are obtained according to equation (18).

$u\approx NU$ (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 _^V1 and_^V2 . It is assumed that _^V2 suffers an increment of rigid body displacement (e) regarding _^V1 .

In order to obtain the same strain for both subdomains, the interpolations of lines 1 and 2 shown in Figure 2a) are adopted for _^V1 and _^V2 . Therefore, the displacements for each subdomain considering the line of discontinuity becomes:

$u_1=N\left(U-\varphi e\right)$ (19)

and

$u_2=N\left(U-\varphi e\right)+e$ (20)

where ϕ =(0 1)^T. Deriving the equations (19) and (20), the same incremental strain is obtained. Therefore, for any point into subdomains _^V1 and _^V2 .

$\epsilon =B\left(U-\Phi e\right)$ (21)

For the two-dimensional case (Figure 2b)) the concepts presented above can be generalized considering:

$e=d_C=Re\text{'}$ (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’):

$e^{\prime }=\left[\begin{array}{c}{u}^c\\ v^c\end{array}\right]$ (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:

$\varphi ={\left[\begin{array}{cc}{\varphi }_1& {\varphi }_2\end{array}\begin{array}{cc}{\varphi }_3& {\varphi }_4\end{array}\right]}^T$ (24)

where each submatrix ϕ_n is a matrix of 2x2 dimension, defined according to:

${\varphi }_n=\{\begin{array}{c}0to{\text{V}}_1\\ Ito{\text{V}}_1\end{array}$

where I represents the identity matrix.

Therefore, the displacements u₁ and u₂ takes the form.

$u_1=N\left(U-\varphi Re\text{'}\right)$ (25)

and

$u_2=N\left(U-\varphi Re\text{'}\right)+Re\text{'}$ (26)

where

$d_N=U-\varphi Re\text{'}$ (27)

which is the displacement that causes strains.

For two points near the line of discontinuity, one on _^V1 and another on _^V2, the displacement is given by:

$u^{+}-u^{-}=Re\text{'}$ (28)

The strain for the two-dimensional case is:

$\epsilon =B\left(U-\varphi Re\text{'}\right)=B{d}_N$ (29)

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 [¹⁶], we have:

$\underset{S_L}{{\displaystyle \int }}tdS=\underset{S_L}{{\displaystyle \int }}{P}^T\sigma dS$ (30)

The matrix P selects the stress components that will be transmitted through the crack.

The constitutive relation for the localization line is defined as:

$\begin{array}{l}t=D^{cr}e\text{'}\\ \left[\begin{array}{c}{t}_1\\ t_2\end{array}\right]=\begin{array}{cc}{E}_t\left(u^c\right)& 0\\ 0& G_t\end{array}\left[\begin{array}{c}{u}^c\\ v^c\end{array}\right]\end{array}$ (31)

where _^Et (_^uc ) is adjusted according to the fracture energy,_^Gf , and _^Gt is related to the shear modulus.

For the element domain:

${}_{}{}^{t+\text{Δ}t}\sigma ={}_{}{}^t\sigma +D^e\Delta \epsilon$ (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.

For the elements in the region where the localization does not occur, the same relation presented in (32) is used.

Considering the strain in (29), replacing (32) in (12):

$\underset{V}{{\displaystyle \int }}{B}^T{D}^{e}B\left(U-\varphi Re\text{'}\right)dV={}_{}{}^{t+\text{Δ}t}f{}_{EXT}-{}_{}{}^{t}F$ (33)

where

${}_{}{}^{t+\text{Δ}t}f{}_{EXT}=\underset{S}{{\displaystyle \int }}{N}^T{}_{}{}^{t+\text{Δ}t}pdS$ (34)

and

${}_{}{}^{t}F=\underset{V}{{\displaystyle \int }}{B}^T{}_{}{}^t\sigma dV$ (35)

Replacing (32) and (31) in the additional equilibrium condition (30), we have:

$\left(\underset{S_L}{{\displaystyle \int }}{P}^T{D}^{e}BdS\right)U=\left[\underset{S_L}{{\displaystyle \int }}\left(P^T{D}^{e}B\varphi R+D^{cr}\right)dS\right]e\text{'}$ (36)

where

$S_{uu}=\underset{S_L}{{\displaystyle \int }}{P}^T{D}^{e}BdS$ (37)

and

$S_{cc}=\underset{S_L}{{\displaystyle \int }}\left(P^T{D}^{e}B\varphi R+D^{cr}\right)dS$ (38)

Therefore, replacing (37) and (38) in (36) and isolating the additional displacement, we obtain:

$e^{\prime }=S_{cc}^{-1}{S}_{uu}U$ (39)

From (39) we have the condensation of the internal degrees of freedom of equation (33) resulting in:

$\left[\underset{V}{{\displaystyle \int }}{B}^T{D}^{e}B\left(I-\varphi R{S}_{cc}^{-1}{S}_{uu}\right)dV\right]U={}_{}{}^{t+\text{Δ}t}f{}_{EXT}-{}_{}{}^{t}F$ (40)

where

${}_{}{}^{t}K{}^{*}=\underset{V}{{\displaystyle \int }}{B}^T{D}^{e}B\left(I-\varphi R{S}_{cc}^{-1}{S}_{uu}\right)dV$ (41)

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 [¹²], 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:

${}_x{}^{t}K{}^{*}U={}_{}{}^{t+\text{Δ}t}f{}_{EXT}-{}_{}{}^{t}F$ (42)

To obtain the equilibrium at time $t+\Delta 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:

$\begin{array}{c}{}_x{}^{t}K{}^{*}{U}^{\left(i\right)}={}_{}{}^{t+\text{Δ}t}f{}_{EXT}-{}_{}{}^{t}F{}^{\left(i-1\right)}\\ U^{\left(i\right)}=U^{\left(i-1\right)}+\Delta U^{\left(i\right)}\end{array}$ (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 [⁸].

To obtain the equilibrium within the element, equation (30), presents some particularities. Reference [²⁰] 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

${}_{}{}^{t+\text{Δ}t}e{\text{'}}_{\left(0\right)}^{\left(i\right)}={}_{}{}^t{e}^{\prime }$ (44)

b) From the displacement obtained by the global equation system, determine the displacement that causes strain in the element

${}_{}{}^{t+\text{Δ}t}d{{}_{Ne}}_{\left(k\right)}^{\left(i\right)}={}_{}{}^{t+\text{Δ}t}U{}^{\left(i\right)}-{}_{}{}^{t+\text{Δ}t}\varphi Re{\text{'}}_{\left(k\right)}^{\left(i\right)}$ (45)

c) Calculate the incremental stress in the element domain, using (32)

d) Calculate the incremental stress in the discontinuity line using (31)

e) Obtain the increment of crack opening

$\Delta e{\text{'}}_{\left(k\right)}^{\left(i\right)}=S_{cc}^{-1}\left[\underset{S_L}{{\displaystyle \int }}{P}^T{\sigma }_{\left(k\right)}^{\left(i\right)}dS-\underset{S_L}{{\displaystyle \int }}{t}_{\left(k\right)}^{\left(i\right)}dS\right]$ (46)

f) Update the crack opening value

${}_{}{}^{t+\text{Δ}t}e{\text{'}}_{\left(k+1\right)}^{\left(i\right)}={}_{}{}^{t+\text{Δ}t}e{\text{'}}_{\left(k\right)}^{\left(i\right)}+\text{Δ}e{\text{'}}_{\left(k+1\right)}^{\left(i\right)}$ (47)

g) Verify convergence

$if(\left|\text{Δ}e\text{'}\right|<Tolerance)ok$ (48)

$elsek=k+1$ (49)

h) Return to item b)

5. Constitutive relations and crack propagation criterion

The constitutive relations used in this work are presented in Figure 3.

Figure 3
Constitutive relations: a) linear elastic; b) linear softening; c) exponential softening

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).

$\sigma =D^e\epsilon$ (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 [²¹]. 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),

$w_{\text{max}}=\frac{2G_f}{f_t}$ (51)

The equation that represents this constitutive relation is given by:

$\sigma =f_t\left(1-\frac{w}{w_{\text{max}}}\right)$ (52)

The exponential constitutive relation, schematized in Figure 3c), is based on the model presented by [¹²]. The equations that represent the exponential softening curve is:

$\sigma =f_t\text{exp}\left(-aw\right)$ (53)

$a=\alpha f_t/G_f$ (54)

A limit is imposed for the factor exp(-aw) = 0.05 in order to obtain a maximum crack opening value (_^wmá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 one of the principal stresses (in modulus), |σ₁| or |σ₂|, reaches the tensile strength (ft), the crack opens [²²].

6. Numerical examples

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

6.1 Four-point bending test

In this first example, a simply supported beam loaded by two symmetric disposed forces, tested experimentally by [²³], 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.

Figure 4
Beam geometry and boundary conditions - four-point bending test

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.

Table 1
Meshes for the four-point bending test

Name	Total of elements	Elements at crack line
M1	423	15
M2	548	20
M3	690	30
M4	822	30

Table 2
Input parameters for the four-point bending test

Geometric properties
Beam
Length	450 mm
Height	100 mm
Thickness	50 mm
Notch
Length	2,5 mm
Height	10 mm
Thickness	50 mm
Material properties
Tensile strength (ft)	3,0 MPa
Young modulus (E)	38000 MPa
Poisson coefficient (ν)	0,2
Fracture energy (Gf)	0,06 N/mm

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

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.

Figure 6
Constitutive models to the crack line

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 [²³] 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 [²³].

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

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.

Figure 8
Load x Deflection of the beam in the load points for the 4 meshes. (R1 and R2 results are from reference [²³])

6.2 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 [²⁴] 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.

Table 3
Meshes for the three-point bending test

Name	Total of elements	Elements at crack line
M1	381	15
M2	721	25
M3	811	25

Figure 9
Beam geometry and boundary conditions - three-point bending test

Figure 10 shows the meshes used with discretization detail at the notched region.

Figure 10
Discretization of the meshes with detail in the region of the notch - three-point bending test

Table 4 shows the input data used for this example and correspond to experiments done by [²⁴].

Table 4
Input parameters for the three-point bending test

Geometric properties
Beam
Length	2000 mm
Height	200 mm
Thickness	50 mm
Notch
Length	20 mm
Height	100 mm
Thickness	50 mm
Material properties
Tensile strength (ft)	3,33 MPa
Young modulus (E)	30000 MPa
Poisson coefficient (ν)	0,2
Fracture energy (Gf)	0,124 N/mm

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 softening 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.

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

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 present results are satisfactory. Numerical results tend to converge to a smoother response with the increase of refinement.

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

7. 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 [¹²]. 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 [¹²] appeared to be the best solution because the peak stress was correctly estimated and the softening curve presented a more realistic behavior, when compared with experimental data. In the linear model [²¹], 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 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.

8. Acknowledgements

The authors thank the Brazilian Government for support through CAPES and CNPq Fellowship and UFRGS for support.

9. References

• ^[1]
  DRIEMEIER L. Contribuição ao estudo da localização de deformações com modelos constitutivos de dano e plasticidade, São Carlos, 1999, Tese (doutorado) - Escola de Engenharia de São Carlos, 101p.
• ^[2]
  NEEDLEMAN A. A Continuum Model for Void Nucleation by Inclusion Debonding. Journal of Applied Mechanics, vol.54, n.3, 1987; p.525-531.
• ^[3]
  PLANAS J., ELICES M., GUINEA G. V., GÓMEZ F. J., CENDÓN D. A., ARBILLA I. Generalizations and specializations of cohesive crack models. Engineering Fracture Mechanics, vol.70, n.14, 2003; p.1759-1776.
• ^[4]
  LENS L. N., BITTENCOURT E., D’AVILA V. M. R. Constitutive models for cohesive zones in mixed-mode fracture of plain concrete. Engineering Fracture Mechanics, vol.76, n.14, 2009; p.2281-2297.
• ^[5]
  GEIßLER G., NETZKER C., KALISKE M. Discrete crack path prediction by an adaptive cohesive crack model. Engineering Fracture Mechanics, vol.77, n.18, 2010; p.3541-3557.
• ^[6]
  BAZANT Z., OH B. Crack band theory of concrete. Materials & Constructions, vol.16, n.93, 1983; p.155-177.
• ^[7]
  JIRÁSEK M., Nonlocal models for damage and fracture: Comparison of approaches. International Journal of Solids and Structures, vol.35, n.31-32, 1998; p.4133-4145.
• ^[8]
  D’AVILA V. M. R., BRISOTTO D. S., BITTENCOURT E . Numerical Simulation of Cracking in Reinforced Members by an Embedded Model. Engineering Computations, v. 25, 2008; p. 739-763.
• ^[9]
  HARIRI-ARDEBILI M. A., SEYED-KOLBADI S. M. Seismic cracking and instability of concrete dams: Smeared crack approach. Engineering Failure Analysis, vol.52, June 2015; p.45-60.
• ^[10]
  JIRÁSEK M., BELYTSCHKO T. Computational resolution of strong discontinuities. In: Fifth World Congress on Computational Mechanics, Viena, 2002, Anais, Viena, 2002, p.7-12.
• ^[11]
  WELLS G. N. Discontinuous modeling of strain localization and failure, Delft, 2001, Tese (doutorado), Delft University of Technology, 169p.
• ^[12]
  DVORKIN E. N., CUITIÑO A. M., GIOIA G. Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions. International Journal for Numerical Methods in Engineering, vol.30, 1990; p.541-564.
• ^[13]
  HUESPE A. E., NEEDLEMAN A., OLIVER J., SÁNCHEZ P. J. A finite thickness band method for ductile fracture analysis. International Journal of Plasticity, vol.25, n.12, 2009; p.2349-2365.
• ^[14]
  ZHANG Y., LACKNER R., ZEIML M., MANG H. A. Strong discontinuity embedded approach with standard SOS formulation: Element formulation, energy-based crack-tracking strategy, and validations. Computer Methods in Applied Mechanics and Engineering, vol.287, 2015; p.335-366.
• ^[15]
  JIRÁSEK M. Comparative study on finite elements with embedded discontinuities. Computer methods in applied mechanics and engineering, vol.188, n.1-3, 2000; p.307-330.
• ^16]
  DA SILVA C. Z. Uso de descontinuidades fortes na simulação de problemas de fratura, Porto Alegre, 2015, Dissertação (mestrado) - Programa de Pós-Graduação em Engenharia Civil, Universidade Federal do Rio Grande do Sul, 85p.
• ^[17]
  SORIANO H. L., LIMA S. S. Método de Elementos Finitos em Análise de Estruturas, São Paulo: Edusp, 2003, 608p.
• ^[18]
  SPENCER B. Embedded Crack elements for analysis of reinforced concrete, Boulder, 2000, Ph.D. Thesis Proposal, University of Colorado - Boulder, 38p.
• ^[19]
  OLIVER J., HUESPE A. E., SAMANIEGO E. A study on finite elements for capturing strong discontinuities. International Journal for Numerical Methods in Engineering, vol.56, n.14, 2003; p.2135-2161.
• ^[20]
  DVORKIN E. N., ASSANELLI A. P. 2D finite elements with displacement interpolated embedded localization lines: The analysis of fracture in frictional materials. Computer Methods in Applied Mechanics and Engineering, vol.90, n.1-3, 1991; p.829-844.
• ^[21]
  HILLERBORG A., MODÉER M., PETERSSON P.-E. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research, vol.6, n.6, 1976; p.773-781.
• ^[22]
  HIBBELER R. C. Resistência dos Materiais, São Paulo: Pearson Education, 7th ed, 2010, 688p.
• ^[23]
  LOTFI H. R., SHING P. B. Embedded Representation of Fracture in Concrete with Mixed Finite Elements. International Journal for Numerical Methods in Engineering, v.38, 1995; p.1307-1325.
• ^[24]
  PETERSSON P.-E. Crack growth and development of fracture zones in plain concrete and similar materials, Lund, 1981, Report TVBM 1006 - Division of building materials, Lund Institute of Technology, 174p.

Publication Dates

• Publication in this collection
  Apr 2018

History

• Received
  26 May 2017
• Accepted
  11 Jan 2018