Kitty Live Streaming - Random Video Chat | Download Mortal Kombat Unchained PSP, PPSSPP, Android | PLAY NOW

Strong discontinuity embedded approach with standard SOS formulation: Element formulation, energy-based crack-tracking strategy, and validations

Strong discontinuity embedded approach with standard SOS formulation: Element formulation, energy-based crack-tracking strategy, and validations

Available online at www.sciencedirect.com ScienceDirect Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366 www.elsevier.com/locate/cma Strong dis...

2MB Sizes 0 Downloads 0 Views

Recommend Documents

The Strong Discontinuity Approach as a limit case of strain localization in the implicit BEM formulation
The Implicit Formulation of the Boundary Element Method (BEM) is used to deal with physically non-linear 2D-problems in

A discrete strong discontinuity approach
In this paper, strong discontinuities are embedded in finite elements to describe fracture in quasi-brittle materials. A

Strong sampling surfaces formulation for layered shells
The strong SaS formulation for the three-dimensional (3D) stress analysis of layered shells is based on a new concept of

A finite element formulation for thermal stress analysis. Part II: Finite element formulation
This paper deals with the application of the finite element method to the variational principle derived in Part I of thi

Multi-materials with strong interface: Young measures formulation
We propose a model of a multi-material with strong interface, whose thickness and stiffness are of order ε and 1ε , whic

Characteristics of methodologies for manufacturing strategy formulation
Although the need for companies to develop competitive manufacturing strategies is widely accepted, the processes or met

Finite element formulation including interlaminar stress calculations
A quadrilateral plate element is developed on the basis of utilizing the compatibility equations to obtain the in-plane

Finite element formulation of the finite rotation solid element
The paper describes the development of an 8-node solid finite element capable of undergoing both large displacements and

A probabilistic reasoning model: Formulation and control strategy
It has been recognized that past experiences of a decision maker often plays a pivotal role in solving new problem insta

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366 www.elsevier.com/locate/cma

Strong discontinuity embedded approach with standard SOS formulation: Element formulation, energy-based crack-tracking strategy, and validations Yiming Zhang a,b , Roman Lackner a,∗ , Matthias Zeiml a , Herbert A. Mang b,c a Material-Technology Innsbruck (MTI), University of Innsbruck, Technikerstraße 13, 6020 Innsbruck, Austria b Institute for Mechanics of Materials and Structures (IMWS), Vienna University of Technology, Karlsplatz 13/202, 1040 Vienna, Austria c Department of Geotechnical Engineering, Tongji University, Siping Road 1239, Shanghai, PR China

Received 18 July 2014; received in revised form 16 January 2015; accepted 2 February 2015 Available online 9 February 2015

Abstract The Strong Discontinuity embedded Approach (SDA) has proved to be a robust numerical method for simulating fracture of quasi-brittle materials such as glass and concrete. Three different numerical formulations are used in the SDA: (i) Statically Optimal Symmetric (SOS), (ii) Kinematically Optimal Symmetric (KOS), (iii) Statically and Kinematically Optimal Nonsymmetric (SKON). While the SOS formulation (standard version) of the SDA is simple for coding and provides good numerical stability (considering standard Galerkin method), several researchers have pointed out that this formulation encounters serious stress-locking. In this paper, an SOS formulation of the SDA is presented, considering elements with linear and quadratic interpolation for the displacement field. The performance of the proposed model of crack simulation is investigated by re-analysis of a pulling test, a three-point bending test, and an L-shaped panel with prescribed crack paths. The obtained results show that elements with linear interpolation encounter significant stress-locking, whereas elements with quadratic interpolation give good results without locking. Based on the eight-node quadrilateral element, which showed the best performance in the aforementioned re-analysis of the tests, an energy-based crack-tracking strategy, originally used in the framework of the XFEM, is modified and implemented into the SDA model. Several numerical benchmark tests, including an L-shaped panel test, a tension-shear test, a notched panel test, four-point bending tests with single and double notches, and a pull-out test are performed, illustrating the good performance with respected to the SDA model and the robustness of the proposed crack-tracking strategy. c 2015 Elsevier B.V. All rights reserved. ⃝

Keywords: Fracture; Quasi-brittle material; Strong Discontinuity embedded Approach (SDA); Statically Optimal Symmetric (SOS) formulation; Energy-based crack-tracking

1. Introduction Models for simulating fracture of quasi-brittle materials such as concrete in the framework of the Finite Element Method (FEM) are existing for more than half a century [1,2]. Traditional approaches include smeared-crack models ∗ Corresponding author. Tel.: +43 0 512 63500.

E-mail address: [email protected] (R. Lackner). http://dx.doi.org/10.1016/j.cma.2015.02.001 c 2015 Elsevier B.V. All rights reserved. 0045-7825/⃝

336

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

and discrete-crack models. For smeared-crack models, fracture is represented by reduction of the local stiffness and strength at Gauss points. Discrete-crack models, on the other hand, employ interface elements, such as the element proposed in [3,4], allowing the cracks to propagate along the boundaries of the finite elements. Though both models have been validated and applied successfully to simulation of fracture of quasi-brittle materials [5,6], they have several shortcomings, such as stress-locking (spurious stress transfer across a widely open crack) [7] and mesh-bias. In recent years, numerical models with cracks embedded in finite elements were developed, considering either nodal enrichment, such as the eXtended Finite Element Method (XFEM) [8] and the Partition of Unity Finite Element Method (PUFEM) [9] or elemental enrichment such as the Strong Discontinuity embedded Approach (SDA), with both models of enrichment giving similar results. However, nodal enrichment requires more computing time [10–12]. Considering the SDA, elements which locally capture discontinuities were first introduced in [13] and applied in [14–17]. In [18], the notion of Enhanced Assumed Strains (EAS) was introduced, linking SDA to classical models in plasticity and continuum damage mechanics [19,20]. In [21], a finite width of the localization zone, h, was introduced, taking into account the experimentally-observed fracture process zone [22–24]. However, in the course of cracking of quasi-brittle materials, the fracture process zone (h > 0) transforms into a macro crack (h = 0) [25], making the SDA a suitable approach to describe the fracture process. The crack inside the element is locally embedded, contributing to the element stiffness matrix and deformation response. Hereby, three different modes of condensation can be distinguished, resulting in the following formulations of the SDA [26]: • Kinematically Optimal Symmetric (KOS), see [27,28], • Statically and Kinematically Optimal Nonsymmetric (SKON), see [7,14,21,29,30], • Statically Optimal Symmetric (SOS), see [17,31,32]. Although standard KOS formulation describes the kinematics satisfactorily, giving a symmetric stiffness matrix, it does not satisfy the L 2 -orthogonality condition [28] and is not consistent with the Hu–Washizu variational principle [33]. Another shortcoming of this formulation (standard version) is the assumption of the tractions along the discontinuity to be equal to the nodal force divided by the length of the band (kinematic condition), violating the natural traction continuity condition (static condition) [26]. The standard SKON formulation provides the correct transfer of the stresses between the cohesive zone, giving a nonsymmetric stiffness matrix. This formulation combines the optimal static and kinematic equations, leading to an improved numerical performance [26], which is the most successful formulation in last decades. If this formulation is based on the standard Galerkin method, the numerical stability and convergence rate of which strongly depend on the angle between the vector describing the crack opening and the vector describing strain localization [34]. Petrov–Galerkin method shall be used for assuring better numerical stability and convergence rate [35,36]. This formulation (standard version) is not consistent with the Hu–Washizu variational principle [28]. The standard SOS formulation is based on a standard Galerkin method, giving a symmetric stiffness matrix [26]. This formulation (standard version) satisfies the L 2 -orthogonality condition and is consistent with the Hu–Washizu variational principle [28,37]. Though this formulation captures the local internal equilibrium (static condition), it cannot reproduce full stress relaxation around a widely open crack (kinematic condition), resulting in serious stresslocking [26,28]. The success of this formulation highly depends on the ability to eliminate the locking. Considering standard versions of these three formulations, the complexity of coding is similar. However, in order to satisfy the L 2 -orthogonality condition and to be consistent with the Hu–Washizu variational principle, modifications of the KOS formulation such as given in [28] shall be considered. As regards the SKON formulation, Petrov–Galerkin method shall be used in order to improve the numerical stability [35,36]. These aspects result in an increased programming effort for coding the KOS and the SKON formulations compared to the standard SOS formulation. In this paper, the authors consider a standard version of SOS formulation of the SDA. Benchmark tests, such as a pulling test, a three-point bending test, and an L-shaped panel test are considered in order to investigate the performance of the finite elements with linear and quadratic interpolation for the description of the displacement field (see Fig. 1). For predicting the propagation of cracks, an energy-based crack-tracking strategy is introduced, complementing the proposed SDA formulation, which is finally applied to several benchmark tests. This paper is organized as follows: in Section 2, an SOS formulation of the SDA, the constitutive relations, and the algorithmic formulation will be presented, assuming the direction of the crack opening to be known. The performance of different discretizations considering triangular and quadrilateral elements will be investigated in Section 3. In Section 4, the crack-tracking strategy as well as the corresponding smoothing technique for the SDA will

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

337

Fig. 1. Elements considered in this paper: (a) 3-node triangular element (T3), (b) 6-node triangular element (T6), (c) 4-node quadrilateral element (Q4), (d) 8-node quadrilateral element (Q8).

Fig. 2. Domain Ω with a discontinuity L.

be introduced. In Section 5, the performance of the SDA and the crack-tracking strategy is assessed by the analysis of an L-shaped panel test, a tension-shear test, a notched panel test, four-point bending tests with single and double notches, and a pull-out test, considering different discretizations. The so-obtained results, including load–displacement curves and crack paths, are compared with experimental results and results obtained from different numerical methods. This paper closes with concluding remarks given in Section 6. 2. Basics of SDA and finite-element implementation 2.1. Kinematics The SDA adopted in this paper was presented in [38,39]. A 2D domain Ω is separated by a discontinuity L into two sub-domains Ω − and Ω + , see Fig. 2. Although the displacement field u of Ω is discontinuous along L, resulting in the singularity of ∇u on L, a sub-domain Ωϕ is introduced to avoid mathematical problems with the singularity. Considering the finite-element implementation, Ωϕ is assembled by the elements passed by L, On the basis of the assumptions of strong discontinuity, enhanced assumed strains only exists inside Ωϕ , the displacement field of which can be described by ¯ u(x) = u(x) + [Hs (x) − ϕ(x)][[u]],

(1)

¯ where u(x) is the regular part, Hs (x) is the Heaviside function, and ϕ(x) is a smooth derivable function, see Fig. 3 [37], with  ∀x ∈ Ω − \ Ωϕ− , 0 (2) ϕ(x) = 1 ∀x ∈ Ω + \ Ωϕ+ ,  continuous function ∀x ∈ Ωϕ .

338

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

Fig. 3. Decomposition of the displacement field u by introduction of the function ϕ.

Hence, the strain field within Ωϕ is obtained as [30] ε(x) = ∇ S u(x) =

¯ ∇ S u(x) − ([[u]](x) ⊗ ∇ϕ) S +    ε¯ (x), ∀x ∈ Ωϕ \ L

δL ([[u]](x) ⊗ n) S ,    εδ (x), ∀x ∈ L

(3)

where (·) S denotes the symmetric part of the tensor and δL stands for the Dirac-delta distribution. In Eq. (3), ∇[[u]](x) = 0 is assumed, resulting from the underlying finite element implementation [38]. Furthermore, the regular part ε¯ (x) can be decomposed into compatible strains and enhanced strains [18], giving ε¯ (x) =

¯ − ([[u]](x) ⊗ ∇ϕ) S . ∇ S u(x)       compatible strains enhanced strains

(4)

In case the crack direction is known, the normal unit vector n and the parallel unit vector t of the failure surface, with n · t = 0, are given (see Fig. 2). Accordingly, the displacement jump at the crack, [[u]](x), is obtained as [[u]] = ζn (x)n + ζt (x)t, with ζn (x) and ζt (x) as the crack opening (crack width) along n and t, respectively. Thus, Eq. (4) can be rewritten as   ¯ ε¯ (x) = ∇ S u(x) − (n ⊗ ∇ϕ) S ζn (x) + (t ⊗ ∇ϕ) S ζt (x) ,    (5)     ε εt where εt is the total-strain tensor and  ε is the enhanced strain tensor. In case the crack direction is fixed, n˙ = 0 and ˙t = 0 and, thus, the strain rate ε˙¯ (x) is determined only by ε˙ t , ζ˙n , and ζ˙t . Moreover, considering a uniform displacement jump within one finite element (for T6, Q4, and Q8 elements, this (e) (e) technique refers to reduced integration), [[u]](x) becomes [[u ]](e) = ζn n + ζt t, giving the strain field ε¯ (e) (x) as: ε¯ (e) (x) ≈

ne  

(e)

∇ Ni

⊗ ui

S

  (e) , − (n ⊗ ∇ϕ (e) ) S ζn(e) + (t ⊗ ∇ϕ (e) ) S ζt

(6)

i=1

where n e is the number of nodes of element “e”, ui is the displacement vector of node “i”. Considering T6, Q4, and Q8 elements, reduced integration will be considered for elements experiencing cracking. This technique was suggested in other SDA models such as [40,41].1 In Appendix A, the authors also give an 1 In [41], the notion “strain-injection” is proposed by using discontinuous bifurcation analysis for specifying the integration mode (reduced vs. full) in each load step. In this paper, on the other hand, once an element experiences cracking, the reduced integration will always be considered in this element.

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

339

Fig. 4. Geometric representation of ζn , ζt , and ζeq for β = 1.

example which shows that full integration for elements experiencing cracking results in stress-locking at the element level. As regards the KOS and the SKON formulation, at first, ϕ (e) is determined [42]. Then ϕ (e) is used to compute ∇ϕ (e) [38]. For the SOS formulation adopted in this paper, on the other hand, ∇ϕ (e) is considered as a vector parallel to n, and directly derived by accounting for the conservation of cohesive energy, as will be presented later. 2.2. Constitutive relation and algorithmic formulation 2.2.1. Constitutive relation Considering the relationship of the crack opening and the traction between the surface of the crack, a mixed-mode traction–separation [43,44] is adopted. The equivalent crack opening (crack width) ζeq is defined as follows:  (7) ζeq = ζn2 + β 2 ζt2 , where β governs the contribution of Mode I and Mode II crack opening, respectively, on ζeq . For β = 1, ζeq is equal to the face-to-face distance of the crack, see Fig. 4. According to [45,46], β = 1 is well-suited for quasi-brittle materials such as concrete. Therefore, it will be considered throughout this paper. The traction components along n and t are defined as Tn and Tt , obtained as follows [44]: Tn =

Teq ζn , ζeq

Tt =

Teq ζt ζeq

with

(8) 

Teq

 ft = f t exp − ζeq , Gf

where f t is the uniaxial tensile strength and G f is the fracture energy. For the initiation of cracking, the Rankine criterion is used, reading φRK (σ ) = (n ⊗ n) : σ − f t < 0, in which σ is the stress tensor. After cracking, the corresponding yield functions for mixed-mode failure is considered as φn = (n ⊗ n) : σ − Tn = 0,

(9)

φt = (t ⊗ n) : σ − Tt = 0.

2.2.2. Algorithmic formulation Based on the classical non-associated plasticity, as proposed in [38] in the context of SDA, the authors consider the stress state of the domain at time step i as   σ i = C : εit −  εi , with  εi = (n ⊗ ∇ϕ

(10) (e) S

) ζn,i + (t ⊗ ∇ϕ

where C is the elasticity tensor.

(e) S

) ζt,i ,

340

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

The corresponding yield functions φn,i = (n ⊗ n) : σ i − Tn (ζn,i , ζt,i ) = 0,

(11)

φt,i = (t ⊗ n) : σ i − Tt (ζn,i , ζt,i ) = 0,

are used for determination of (ζn,i , ζt,i ). t ,ζ For this step i, the known state variables at a specified Gauss point are (εit , εi−1 n,i−1 , ζt,i−1 ), whereas the unknown state variables are (ζn,i , ζt,i ), with ζn,i = ζn,i−1 + 1ζn,i and ζt,i = ζt,i−1 + 1ζt,i . In the general returnmapping algorithm [47], the trial state of stress σ itr is introduced as   σ itr = C : εit − (n ⊗ ∇ϕ (e) ) S ζn,i−1 − (t ⊗ ∇ϕ (e) ) S ζt,i−1 . (12) Based on the trial state, the unknown incremental state variables, 1ζn,i and 1ζt,i , are determined from the yield function, reading2  (n ⊗ n) : σ itr − g11 1ζn,i − g12 1ζt,i − Tn = 0 (13) (t ⊗ n) : σ itr − g21 1ζn,i − g22 1ζt,i − Tt = 0 with g11 = (n ⊗ n) : C : (n ⊗ ∇ϕ (e) ) S ,

g12 = (n ⊗ n) : C : (t ⊗ ∇ϕ (e) ) S ,

g21 = (t ⊗ n) : C : (n ⊗ ∇ϕ (e) ) S ,

g22 = (t ⊗ n) : C : (t ⊗ ∇ϕ (e) ) S .

(14)

The relation between the tractions differentials, dTn and dTt , and the crack-opening parameters, dζn and dζt , is given by [44]     dTn dζn =D , (15) dTt dζt with ζn2 ζt2  β2G − f ζ t eq f    G f + f t ζeq ζn ζt G f f t ζeq  D=−

β 2 f t Teq ζeq

 G f + f t ζeq ζn ζt G f f t ζeq    β 2 ζt2 ζn2  − Gf f t ζeq

for elastoplastic loading (16)

and  D = Teq

1 0

0 β2

 for elastic unloading,

leading to the elastoplastic tangent moduli Cep , given as Cep =

dσ i = C − S, dεit

(17)

where  S = C : (n ⊗ ∇ϕ (e) ) S

 (t ⊗ ∇ϕ (e) ) S : (G − D)−1 :



 n⊗n :C (t ⊗ n) S

with  G=

g11 g21

 g12 . g22

2 For loading, i.e. if the condition (n ⊗ n) : σ tr − T > 0 is fulfilled, Eq. (13) will be solved by the Newton–Raphson method. n

(18)

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

341

Fig. 5. Strain localization: domain (element) with a discrete crack and domain (element) with SDA localization.

The condition number κ(G) of the matrix G will increase if the angle enclosed by n and ∇ϕ (e) is increasing [48]. Thus, the maximum numerical stability of the SDA model (based on standard Galerkin method) will be obtained if ∇ϕ (e) remains parallel to n [34]. 2.3. SOS formulation For ∇ϕ (e) remaining parallel to n in the SOS formulation, ∇ϕ = l n,

(19)

where l is a scalar quantity. To obtain the l, it is assumed that the dissipated energy E is the same for the SDA localization (see Fig. 5).        ζn

ζt

Tn dζn +

E= A

0

 ε

!

Tt dζt d A = Ω

0

σ : d ε dΩ .

(20)

0

  Introducing ∇ϕ = l n into  ε provides  ε = l (n ⊗ n)ζn + (t ⊗ n) S ζt . Since reduced integration technique (considering only the center Gauss point) is used for elements experiencing cracking, Eq. (20) comes into  ζn   ζt E = A Tn dζn + Tt dζt 0 !

0



=Vl 0

ζn

σ : (n ⊗ n)dζn +



ζt

 σ : (t ⊗ n) S dζt ,

(21)

0

where A is the area of the crack and V is the volume of the finite element experiencing cracking. Moreover, according to Eq. (9),3 Tn = (n ⊗ n) : σ and Tt = (t ⊗ n) : σ . Thus l = A/V can be obtained directly from Eq. (21). As only the central Gauss point is considered in the cracked finite elements, the crack area A is determined considering a crack passing through the center Gauss point, as illustrated in Fig. 6.4 Considering the element with quadratic interpolation of the displacement field (T6 and Q8), crack path may appear as curved. However, only straight crack are considered in this paper for T6 and Q8 (only crack plane are considered for H20 and H27 in Appendix E), which is more convenient for accounting the dissipation of the energy when using reduced integration. The proposed SOS formulation, considering Eq. (19) and l = A/V , deduced from equivalence of dissipated energy, is very similar to the original SOS formulation proposed in [17,32] and summarized in [26], which belongs to standard SOS formulation. 3 As (σ ) S equals σ , (t ⊗ n) : σ equals (t ⊗ n) S : σ . 4 As can be seen in Fig. 6, there is a significant discrepancy between the area of the “real crack” and of its corresponding FE modeling A. The authors want to emphasize that the direct introduction of the area of the “real crack” into Eq. (21) would be illogical. Because the stress state of the “real crack” is unknown in FEM, the corresponding cohesive energy is not available.

342

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

Fig. 6. Determination of the crack area A: (a) element experiencing cracking, (b) determination of A from a parallel crack located at the center Gauss point.

3. Numerical examples for investigating different elements As crack-tracking strategies [44,49,50] employ information from the actual state of deformation, the crack paths obtained from such strategies depend also on the numerical properties of the elements, e.g. whether or not there is locking. Therefore, for the purpose of a fair comparison between different types of elements, crack-tracking strategies are not considered in this chapter, where all of the crack paths are prescribed. 3.1. Pulling test The set-up of the pulling test, the material properties and the considered discretizations of the test specimen are shown in Fig. 7. The finite element “e1” is the root element of the prescribed crack. Hence, the crack will propagate in the order e1, e2, e3, to finally e4. Fig. 8 shows force–displacement curves, considering different types of elements. It is seen that the elements with linear interpolation functions (T3 and Q4) overestimate the cohesive stress inside the crack, while T6 underestimates the cohesive stress and Q8 gives almost the same result as the analytical solution. In order to provide more insight, the development of the crack opening ζeq is given in Fig. 9. Considering two elements along the pulling direction (“e1” and “e2”; “e3” and “e4”), large values of ζeq for both elements are obtained for T3, T6, and Q4, indicating two macro cracks. For T3 and Q4, this property causes overestimation of the cohesive

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

Fig. 7. Pulling test: geometric properties, material parameters, and discretizations.

Fig. 8. Pulling test: force–displacement curves for different types of elements.

343

344

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

Fig. 9. Pulling test: values of crack opening ζeq for different types of elements.

stress, see Fig. 8 (T3 and Q4). For T6, these two macro cracks cause underestimation of the cohesive stress, see Fig. 8 (T6). For the Q8 element, which gives only one macro crack on one crack band, good agreement with the analytical solution was found, see Fig. 8 (Q8). Different property of different elements experiencing cracking can also be observed from the deformation plots, see Fig. 10. While linear interpolation of the displacement field causes stress-locking inside two elements along the pulling direction (since the common boundary of these two elements remains linear), quadratic interpolation allows non-linear deformation of the boundary of these two elements, giving better results. On the other hand, Fig. 10 shows an hourglass mode for Q8 element, in case of reduced integration, which is discussed in Appendix B. 3.2. Three-point bending test The set-up of the three-point bending test, the material properties and the corresponding discretizations are shown in Fig. 11.5 A report of the corresponding experimental investigation is contained in [51]. Force–displacement curves and experimental results from [51] are shown in Fig. 12. It is seen that elements with linear interpolation (T3 and Q4) exhibit stress-locking in case of the irregular mesh (TMesh II and QMesh II), whereas elements with quadratic interpolation give better results, although the T6 element shows mesh bias. The Q8 element provides the best numerical results, characterized by no mesh bias. 3.3. L-shaped panel test The set-up of the L-shaped panel test, the loading and the corresponding discretizations are shown in Fig. 13. A report of the corresponding experimental investigations is contained in [52]. The prescribed crack path was obtained by means of an energy-based crack-tracking strategy in the framework of the XFEM [44]. 5 In all the figures in this paper, which show the dimensions of the investigated test specimens, the term “thickness” stands for the dimensions perpendicular to x–y plane. Since the “thickness” is relatively small in comparison to the other dimensions, all investigated problems of this paper, except the pulling test, belong to the category of plane-stress problems.

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

345

Fig. 10. Pulling test: deformation plots for different types of elements (d = 0.2 mm, scale = 1:100).

Fig. 11. Three-point bending test: geometric properties, material parameters, and discretizations.

Force–displacement curves are shown in Fig. 14. They are compared with experimental results [52] and numerical results that were obtained by means of the XFEM [53]. The results indicate that elements with linear interpolation (T3 and Q4) exhibit significant stress-locking, whereas elements with quadratic interpolation do not do so, although the

346

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

Fig. 12. Three-point bending test: force–displacement curves for different types of elements and discretizations.

T6 element underestimates the cohesive stress and the peak load. The Q8 element provides the best numerical results, characterized by excellent agreement with the result via XFEM, and by no mesh bias. In addition, considering the Q8 element, coarse discretizations may commonly used in the FEM analysis. The properties of Q8 element with coarse discretizations are discussed in Appendix C. 4. Energy-based crack-tracking strategy 4.1. Introduction of crack-tracking strategies Crack-tracking strategies are used for predicting the propagation of cracks, based on the history of loading and the deformation state of the system. Except for special formulations of SDA models, based on e.g. rotating cracks [38,54] and multiple cracks [55,56], most SDA formulations require a corresponding crack-tracking strategy for determination of the proper crack path. Crack-tracking strategies commonly used within the SDA and the XFEM are listed in Table 1. Comparisons of numerical results obtained from different crack-tracking strategies can be found in [53,63,62,64]. The local crack-tracking strategy is one of the most widely used tracking strategies. However, when applying this strategy, stress-locking may be encountered. Although this drawback may be mitigated by some simple modifications [57], a zigzag crack path is almost inevitable. In the delayed tracking strategy, based on the observation that the direction of the crack band becomes smoother with further propagation of the crack [58], at first, a smeared-crack model for simulation of the fracture process is adopted. Then, a switch to the SDA model is made. Smooth crack paths can be found in several examples in [58,59]. The drawback of this method is that an extra smeared-crack model must be combined with the SDA, taking into account a transformation from the former to the latter. In the global tracking strategy, an unknown scalar quantity is introduced at each node. It is determined at each load step, with the contour lines of this unknown quantity representing the crack path. For some examples, a smooth

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

347

Fig. 13. L-shaped panel test: geometric properties, material parameters, and discretizations.

crack path is obtained even at the beginning of loading, as observed for the tension-shear test [53] and the four-point bending test with double notches [60]. The only drawback of this method are the additional unknowns, requiring more computing effort. The partial domain tracking strategy is a modified version of the global tracking strategy, which combines the advantages of the local and the global strategy [65]. However, presently this strategy is only suitable for constantstrain triangular elements (T3 element).6 While the XFEM is able to locate the exact position of the crack tips, in case of the SDA, the cracks propagate from one boundary of an element to another one. Hence, the exact position of the crack tip is not accessible by the SDA, resulting in inexact stress states around the crack tip. Accordingly, the level-set method and the maximum energy release method, which are commonly used in the XFEM and highly rely on the stress/strain states around the crack tip, cannot be applied in the framework of the SDA without considerable modifications. Early applications of energy-based strategies of crack propagation in finite-element analysis can be found in [66]. Though based on linear-elastic fracture mechanics, contrary to the maximum circumferential stress criterion or the maximum energy release-rate criterion, the energy-based criterion also allows to predict the direction of the crack propagation for cohesive cracks [62,44]. Moreover, the energy-based crack-tracking method considers the energy of 6 Like the global tracking strategy, an unknown at each node is introduced in the partial domain tracking strategy. It is determined within a specified domain but not within the whole body. For determination of this domain, a formulation is presented, wherein only the T3 element is considered [65].

348

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

Fig. 14. L-shaped panel test: numerically-obtained force–displacement curves for different types of elements, experimental results [52], and numerical results (XFEM) in [53].

Table 1 Crack-tracking strategies.

(a) (b) (c) (d) (e) (f) (g)

Criterion

Smoothing

Suitable for

References

MPS NMPS MPS NMPS MCS MER MPE

Maximum curvature of crack Delayed crack propagation Isotherm Isotherm Level set – –

SDA SDA SDA & XFEM SDA XFEM XFEM XFEM

[57] [58,59] [49,53,60] [39] [50,61] [62] [44]

(a): Local tracking. (b): Delayed tracking. (c): Global tracking. (d): Partial domain tracking. (e): Level-set method. (f): Maximum energy release rate (J-integrals). (g): Energy-based crack-tracking. MPS: Maximum principal stress/strain. NMPS: Maximum principal nonlocal strain. MCS: Maximum circumferential stress. MER: Maximum energy release rate. MPE: Minimization of total potential energy.

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

349

Fig. 15. Deformed body with a crack.

the entire system, which does not require precise information about the stress/strain state around the crack tip, making this method suitable not only for the XFEM but also for the SDA. Hence, this method is proposed in this paper and implemented into the SDA for predicting the crack propagation. 4.2. Energy-based crack-tracking strategy for SDA A deformed body Ω with a crack is considered. The body force is denoted as f, the surface loading as t, and the cohesive force at the surface of the crack as w, see Fig. 15. The mechanical energy Ψ is given as Ψ = U + E − W,

(22)

where U is the elastic strain energy, E is the dissipated energy, and W is the work done by the applied forces. Ψ is a function of the crack-propagation angle θ and of the crack-propagation area A, i.e. Ψ = Ψ (θ, A). On the basis of linear-elastic fracture mechanics, θ and A are determined by minimizing Ψ . U , E, and W in Eq. (22) are given as follows:    ε  t  σ : d ε − ε dΩ , U= Ω

0

ζn

  E=



ζt

Tn · dζn + A

0

  W = Ω

 (23)

Tt · dζt d A, 0

u

  f du dΩ +

∂Ω

0



ut



t dut d∂Ω ,

0

where u is the displacement of a point of the body, ut is the displacement of a point on the surface of the body loaded by t, and A is the surface of the crack.7 Considering Eq. (20), U + E is given by    ε U+E= σ : dε t dΩ , (24) Ω

0

which is used by the authors for obtaining Ψ . Considering the procedure of tracking, in a prescribed angle domain (in this paper, ±30◦ with respect to the last crack direction, which is enlarged if the minimum value of Ψ found at ±30◦ ), several values of θ, as θi , i = 1 . . . n (n = 15 is taken in this paper) are tried as possible angles for further propagation (with fixed number of iteration step, e.g. 10 steps is recommend in [44]). With the corresponding Ψi at hand, the discrete values (θi , Ψi ) are fitted with a continuous function,8 giving access to Ψ = Ψ (θ ). Then, the value of θ, which defines the angle for further 7 In Eq. (23), A seems to have direct influence on E then on Ψ . Nevertheless, the authors want to emphasize, considering  ε =   and the strain localization (see Fig. 5 and Eq. (20)), A is a parameter coherent inside  ε. There is no proportional relationship between A and Ψ , see Appendix D. 8 In [62], fitting with moving least squares was suggested. In [44], it was shown that a polynomial function of order 2–4 gives sufficiently accurate results. A S V (n ⊗ n)ζn + (t ⊗ n) ζt

350

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

Fig. 16. Continuous vs. non-continuous crack path.

Fig. 17. Smoothing technique for energy-based crack-tracking in the SDA.

Fig. 18. Wrong propagation and correction by means of the smoothing method.

propagation of the crack, is obtained at the minimization of Ψ (θ ), i.e. where ∂Ψ /∂θ = 0 and ∂ 2 Ψ /∂θ 2 > 0. An example is given in Appendix D to show the relationship between the direction of the crack and the energy Ψ .

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

351

Fig. 19. L-shaped panel test: crack paths and load–displacement curves obtained from numerical analysis, compared with experimental results presented in [52], and with numerical results based on the XFEM [53].

Fig. 20. Tension-shear test: geometric properties, material parameters, and support conditions.

Considering Fig. 6 and Eq. (24), it is unnecessary to assume continuous crack path since the area of the “real crack” is not used in calculation. However, when using non-continuous crack path (every crack passing through the center Gauss point) such as used in [38,58,67,68], all surrounding elements of the element with crack tips should be checked for determining the next element experiencing cracking. For example, in Fig. 16, the crack propagates from element “e0” to element “e1” (e1 is the element with the crack tip). In present step the crack will propagate to one element neighboring “e1”. When using non-continuous crack path, the total energy considering cracking of e2, e3, e4, and e5 should all be checked respectively, bringing much more computing effort for tracking, while only the total energy considering cracking of e2 needs to be checked when using continuous crack path.

352

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

Fig. 21. Tension-shear test: crack paths and load–displacement curves obtained from numerical analysis considering different values of Fs and different discretizations, compared with experimental results [69] and numerical results obtained by the XFEM [44].

Unlike the XFEM, allowing a continuous propagation of cracks, the energy-based crack-tracking in the SDA is characterized by a mesh bias, as the propagation length is always connected to the element size. As a remedy, a smoothing technique (see Fig. 17) is employed, using similar changes in the crack direction (see θ1S in Step 2 and θ2S

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

353

Fig. 22. Notched panel test: geometric properties, material parameters and loading.

Fig. 23. Notched panel test: crack path and load–displacement curve obtained from numerical analysis, compared with experimental results [70] and numerical results from a KOS [27] and an SKON [71] formulation of the SDA.

in Step 4) for the two elements closest to the crack tip [48]. Hereby, the change in the crack direction is determined from the mean value of the unsmoothed crack directions. Furthermore, the presented smoothing technique is capable of correcting a wrong propagation (sudden turning) shown in Fig. 18. 5. Numerical examples for investigating the SDA model with the crack-tracking strategy 5.1. L-shaped panel test The test set-up, the material parameters, the loading and the corresponding discretizations of the L-shaped panel were shown in Fig. 13. The crack path and the load–displacement curve, obtained from the numerical analysis are illustrated in Fig. 19, showing good agreement with the experimental results [52] and the results based on the XFEM [53]. The location of the crack paths in the upper part of the experimental range can be explained by the small shear stiffness considered in the analysis, where β was set equal to 1, see [44]. Considering the computing time of the L-shaped panel with (i) energy-based crack-tracking strategy tt , (ii) prescribed crack path in Section 3.3 t p , the efficiency of this strategy can be evaluated as (tt − t p )/tt , which is 8.2% in QMesh I and 9.6% in QMesh II.

354

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

Fig. 24. Four-point bending tests with single and double notches: geometric properties and material parameters.

5.2. Tension-shear test The tension-shear test was first investigated experimentally in [69]. The geometric dimensions, the material parameters, the loading and the support conditions of the tension-shear test are shown in Fig. 20. The “Load Path 4” in [69] is considered in this paper. First, the horizontal load Fs is applied and kept constant. Thereafter, Fn is gradually increased, leading to a vertical displacement d and to cracking of the structure. Since different values of Fs result into different crack paths, the tension-shear test is a classical benchmark test for both numerical models and crack-tracking strategies, see [39,60,44]. The crack paths and the load–displacement curves obtained from the numerical analysis are shown in Fig. 21. The comparison with experimental results [69] and numerical results obtained by the XFEM [44] shows that the presented crack-tracking strategy and the SDA model are performing reasonably well. 5.3. Notched panel test The notched panel was experimentally investigated in [70]. The geometric dimensions, the material parameters, and the loading are shown in Fig. 22. The loads F1 and F2 are increased until a given target value is reached. Then, F2 is kept constant, while F1 (or a corresponding displacement d) are further increased. Representing a classical example, the notched panel was also investigated numerically in [12,71].

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

355

Fig. 25. Four-point bending test with a single notch: crack paths, load–displacement curves, and load-CMSD (Crack Mouth Sliding Displacement) curves obtained from numerical analysis, considering two different discretizations, compared with experimental results in [72], and numerical results obtained from an SKON formulation of the SDA [71] and an interface element [74]; remark: In [74], G f was set equal to 150 N/m, resulting in higher values of the peak load and the residual load in the post peak range.

Fig. 26. Four-point bending test with double notches: crack paths and load–displacement curve obtained from numerical analysis, compared with experimental results [73] and numerical results based on interface elements [74].

The crack path and the load–displacement curve obtained from numerical analysis are shown in Fig. 23, indicating reasonably good agreement with experimental results as well as with numerical results obtained by a KOS [27] and an SKON [71] formulation of the SDA.

356

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

Fig. 27. Pull-out test: geometric properties and material parameters.

Fig. 28. Pull-out test: crack path and load–displacement curve obtained from numerical analysis, compared with numerical results obtained by the XFEM [77] and by an SKON formulation of the SDA [78], respectively.

5.4. Four-point bending tests with single and double notches Four-point bending tests with single and double notches [72,73], are considered as relatively complicated benchmark tests for validating numerical models for the simulation of fracture processes (see [71,74–77]). Similar to [74], the steel beam at which the load is applied (see Fig. 24) was not considered in the numerical model. Loading was directly applied at points A and B based on Archimedes’ law of the lever (in SEN: FA = 0.118F, FB = 0.882F; in DEN: FA = 0.167F, FB = 0.833F). Hence, the displacement of point C is determined from the displacements at points A and B. The crack paths and the load–displacement curves obtained from the numerical analysis are shown in Figs. 25 and 26. As regards the four-point bending test with a single notch, the presented results for the crack paths, the peak value of the load, and the post-peak behavior are reasonable, when compared to numerical results obtained by other researchers. Considering the four-point bending test with double notches, the numerical results for the peak value of the load, and the post-peak behavior are similar to the experimental and numerical results in the literature. Although two crack paths are initialized, unlike as in the tension-shear test, only one crack will propagate, which was reported in [60]. 5.5. Pull-out test Because of axial symmetry (see Fig. 27), the pull-out test can be simulated by means of 2D axisymmetric analysis [51]. As pointed out in [51], there are two types of failure modes: tensile failure and compressive failure. However, only tensile failure is considered in this benchmark example. The crack path and the load–displacement curve obtained from numerical analysis are shown in Fig. 28. Consistent with investigations by other authors [77,78], who considered a quarter of the structure, F/4 represents the applied

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

357

Fig. A.29. Single-element (Q8) example: geometrical dimensions, loading conditions, and material properties.

Fig. A.30. Single-element (Q8) example: force–displacement curves illustrating the effect of full and reduced integration.

load in the load–displacement curve in Fig. 28, showing good agreement between the presented results and the ones published in [77,78]. 6. Conclusion In this paper, an SOS formulation of a strong discontinuity-embedded approach was presented. In order to avoid for stress-locking in consequence of the SOS formulation, as reported in the literature, four elements with linear and quadratic interpolation for the displacement field, respectively, were considered, proving the eight-node quadrilateral element (Q8) showing the best numerical performance. Using the Q8 element, the proposed SDA was extended by an energy-based crack-tracking strategy. In addition to crack-tracking, a smoothing technique was proposed, avoiding mesh bias and/or wrong crack propagation. The performance of the proposed formulation, on the one hand, and of the energy-based crack-tracking strategy, on the other hand, was investigated by numerous experimental and numerical validation tests. From comparison of the results, the following conclusions can be drawn: • Using elements with linear interpolation of the displacement field in an SOS formulation of the SDA, results in significant stress-locking. For such elements, SKON and modified KOS formulations are recommended. • Using elements with quadratic interpolation for the displacement field in an SOS formulation of the SDA, provides good results with no stress-locking. Triangular elements with quadratic interpolation (T6) result in underestimation of the cohesive stress and in mesh-bias. Quadrilateral elements with quadratic interpolation (Q8) provide results which agree very well with the analytical solution and with experimental and numerical results reported in the literature. No stress-locking and no mesh-bias were encountered. • The proposed energy-based crack-tracking strategy, which does not require knowledge of the precise stress state around the crack tip, is suitable for both the SDA and the XFEM. • The proposed SOS formulation of the SDA and the energy-based crack-tracking strategy provide reasonable crack paths and load–displacement curves for all numerical examples investigated in this paper. This is a strong indication of the robustness of the SDA model.

358

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

Fig. B.31. Pulling test: cases for investigation of the hourglass mode in the Q8 element: (a) discretizations and prescribed cracks, (b) deformation plots (d = 0.2 mm, scale = 1:100, F: full integration, R: reduced integration), (c) force–displacement curves.

The basic advantages of the proposed SOS formulation of the SDA together with the introduced energy-based crack-tracking strategy, are (i) easy coding, (ii) good numerical stability, and (iii) good performance as regards the load–displacement curves and the crack paths within the re-analysis of the benchmark examples. Considering the extension of the proposed model towards 3D analysis, elements with quadratic interpolation e.g. 20 node or 27 node hexahedron elements shall be used within the SDA, see Appendix E. As regards the underlying cracktracking strategy, 3D benchmark problems shall be investigated (similar to the series of 2D benchmark problems presented in this paper) in order to assess the applicability and performance of the energy-based tracking strategy with respect e.g. to global tracking strategies as proposed in [49]. Acknowledgments The authors gratefully acknowledge financial support by the Austrian Ministry for Transport, Innovation and Technology (bm.vit) within the KIRAS-project (Austrian Security Research Program), project number: 836268 “ Sicherheit von Hohlraumbauten unter Feuerlast” (“Safety of underground structures under fire loading”).

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

359

Distorted case I

root element

Distorted case II

root element

Distorted case III

root element Fig. B.32. Pulling test: cases with distorted discretizations: (a) discretizations and prescribed cracks, (b) deformation plots (d = 0.2 mm, scale = 1:100), (c) force–displacement curves.

Appendix A. Adoption of reduced integration for elements exhibiting cracking As outlined in [38,79], in case of elements with quadratic interpolation of the displacement field (T6 and Q8), the crack opening has different values at different Gauss points, with some Gauss points remaining in the elastic state in the course of cracking. This causes stress-locking, which cannot be eliminated by tracking strategies, as it occurs at the element level. As a remedy, for elements exhibiting cracking, reduced integration with one Gauss point at the element center is applied. The effect of this integration technique is illustrated by means of an example, characterized by two different load cases, for both of which the Q8 element is used, see Fig. A.29. Considering symmetry of both cases, the crack direction is assumed to be parallel to the y-axis. The resulting force–displacement curves are shown in Fig. A.30. For CASE 1, all Gauss points experience the same stress/strain states. Therefore, full and reduced integration give the same results.9 9 Although full integration for CASE 1 gives the same result as reduced integration, the numerical stability of full integration is poorer, since different Gauss points may be loaded differently, which is caused by numerical errors in the course of the iteration process.

360

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

Mesh I

Mesh II

Fig. C.33. L-shaped panel test considering coarse discretizations: Mesh I and Mesh II with corresponding prescribed crack paths.

Mesh I

Mesh II

Fig. C.34. L-shaped panel test considering coarse discretizations: numerically-obtained force–displacement curves, experimental results [52], and numerical results (XFEM) in [53].

Fig. D.35. Example for showing energy-based crack-tracking in the SDA: geometrical dimensions, loading conditions, and material properties.

For CASE 2, the stress/strain states vary inside the element, leading to different crack directions at different Gauss points for the case of full integration.10 Assuming one global crack direction for a single element, full integration results in stress-locking at the element level. Considering the shift from full to reduced integration, the stress state at center Gauss point is monitored in each iteration step (a threshold value for the tolerance of the iteration residual shall be set, avoiding numerical errors in the course of the iteration process). Once the Rankine criterion (see Section 2.2, φ R K (σ ) = (n ⊗ n) : σ − f t < 0) is violated at the center Gauss point, the reduced integration will always be considered for this finite element thereafter, which is irreversible. 10 Rotating localization (directions of cracks will change in every loading step) presented in [38,55], allowing different crack directions at different Gauss points, greatly alleviates this type of stress-locking. Nevertheless, the SDA proposed in this paper is based on fixed localization.

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

361

Fig. D.36. Example for showing energy-based crack-tracking in the SDA: (a) relationship between A and the direction of the crack η, (b) relationship between total energy of the system Ψ and the direction of the crack η. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Appendix B. Discussion of the hourglass mode of the Q8 element, introduced by reduced integration This hourglass mode is illustrated in Fig. 10. The cases shown in Fig. B.31 are considered here. Since reduced integration is only used for cracked elements, the hourglass mode will only occur if the crack passes two elements lying abreast, whereas full integration continues to be used for the elements remaining in the elastic state. Hence, the resulting global deformation is reasonable. This property results the Q8 element very robust with respect to changes of the discretization, see the results for the distorted discretization shown in Fig. B.32. Appendix C. Properties of Q8 element with coarse discretizations As an element with quadratic interpolation of the displacement field, Q8 element is commonly used with relatively coarse discretizations for elastic FEM analysis. Therefore, it is necessary to discuss the properties of Q8 element in coarse discretizations considering elastoplastic FEM analysis. The L-shaped panel is re-analyzed considering coarse discretizations shown in Fig. C.33. The corresponding force–displacement curves are shown in Fig. C.34, indicating oscillating behavior. This can be explained by delayed plastification of elements when adopting reduced integration (elements get cracking only after their center Gauss point indicates the occurrence of crack). No stress-locking is observed. Appendix D. Energy-based crack-tracking in the SDA: an example An example is designed for showing energy-based crack-tracking strategy in the SDA, see Fig. D.35. The prescribed displacements are applied to the structure in i steps. From load step 0 to load step i − 1 the structure keeps as elastic when the cracking (in element “e1”) happens between load step i − 1 and i. After cracking, 10 iteration steps

362

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

Fig. D.37. Example for showing energy-based crack-tracking in the SDA: considering element with high aspect ratio.

Fig. E.38. Elements considered for pulling test in 3D (only visible nodes and surfaces are shown): (a) 8 node hexahedron elements (H8), (b) 20 node hexahedron elements (H20), (c) 27 node hexahedron elements (H27).

are conducted. By assuming different crack direction η, size of A and the corresponding total energy Ψ of the system are obtained, giving relationship shown in Fig. D.36, indicating there is no proportional relationship between A and Ψ . Considering the tracking of the crack, the value of η with the minimum value of Ψ is the direction of the crack. Although there is no proportional relationship between A and Ψ , A has influence on Ψ (see the blue arrows in Fig. D.36). Hence, energy-based tracking considering elements with high aspect ratio may bring unreliable results, which should be avoided in discretization. Fig. D.37 shows an example (material properties and loading conditions are the same as shown in Fig. D.35) in which the value of d2 max was greatly increased, but the angle with minimum value of Ψ changes only a little, indicating the “insensitive” relationship between loading and the direction of the crack considering elements with high aspect ratio. Appendix E. Pulling test in 3D Hexahedron elements with 8 node (H8), 20 node (H20), and 27 node (H27) are used in this section for analyzing the pulling test in 3D, see Fig. E.38. The model and discretization are shown in Fig. E.39, with loading condition and material properties same as shown in Fig. 7. For H20 and H27, the crack propagates in the order e1, e2, e3, to finally e4. For H8, the crack propagates from e2, then to e3, to finally e4.11 11 Making this assumption for H8 is to show that considering SOS formulation of the SDA, even when only one macro crack is formed, the element with linear interpolation for displacement field still encounters serious stress locking.

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

363

Fig. E.39. Model of 3D analysis of pulling test.

Fig. E.40. Pulling test in 3D: force–displacement curves for different types of elements.

Fig. E.41. Pulling test in 3D: values of crack opening ζeq for different types of elements.

The force–displacement curves are shown in Fig. E.40, indicating hexahedron element with linear interpolation (H8) encountering serious stress locking while hexahedron element with quadratic interpolation (H20 and H27) giving reliable results. Furthermore, the development of the cracking opening ζeq are shown in Fig. E.41, indicating only one macro crack is formed in all types of element.12 Considering H8, even while the values of crack opening are 12 For H20 and H27, ζ of e1 element are obtained as small values, but not zero, indicating the cracking of e1. Hence the crack further propagates eq to e2 and e3 elements smoothly.

364

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

reasonable (fulfilling the static condition), cracking of elements “e2” and “e3” is unable to fully release the stress in the surrounding elements (violating the kinematic condition). While the H20 and H27 show good results comparing to the analytical solution. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

D. Ngo, A. Scordelis, Finite element analysis of reinforced concrete beams, J. Am. Concr. Inst. 64 (1967) 152–163. Y. Rashid, Ultimate strength analysis of prestressed concrete pressure vessels, Nucl. Eng. Des. 7 (1968) 334–344. R. Goodman, R. Taylor, T. Brekke, A model for the mechanics of jointed rock, J. Soil Mech. Found. Div. (ASCE) 99 (1968) 637–660. L. Herrmann, Finite element analysis of contact problems, J. Eng. Mech. (ASCE) 104 (1978) 1043–1059. P. Feenstra, R. Borst, A plasticity model and algorithm for mode-I cracking in concrete, Internat. J. Numer. Methods Engrg. 38 (1995) 2509–2529. A. Ingraffea, V. Saouma, Numerical modelling of discrete crack propagation in reinforced and plain concrete, in: G.C. Sih, A. DiTommaso (Eds.), Fracture Mechanics of Concrete, Martinus Nijhoff, Dordrecht, 1985, pp. 171–225. M. Jir´asek, T. Zimmermann, Embedded crack model: I. Basic formulation, Internat. J. Numer. Methods Engrg. 50 (2001) 1269–1290. N. Sukumar, N. Mo¨es, B. Moran, T. Belytschko, Extended finite element method for three-dimensional crack modelling, Internat. J. Numer. Methods Engrg. 48 (2000) 1549–1570. J. Melenk, I. Babuˇska, The partition of unity finite element method: basic theory and applications, Comput. Methods Appl. Mech. Engrg. 139 (1996) 289–314. M. Jir´asek, Computational resolution of strong discontinuities, in: H. Mang, F. Rammerstorfer, J. Eberhardsteiner (Eds.), Fifth World Congress on Computational Mechanics (WCCM V), Vienna University of Technology, 2002. R.I. Borja, Assumed enhanced strain and the extended finite element methods: a unification of concepts, Comput. Methods Appl. Mech. Engrg. 197 (2008) 2789–2803. J. Oliver, A. Huespe, P. Sanchez, A comparative study on finite elements for capturing strong discontinuities E-FEM vs X-FEM, Comput. Methods Appl. Mech. Engrg. 195 (2006) 4732–4752. C. Johnson, R. Scott, A finite element method for problems in perfect plasticity using discontinuous trial functions, in: W. Wunderlich, E. Stein, K.-J. Bathe (Eds.), Nonlinear Finite Element Analysis in Structural Mechanics, Springer, 1981, pp. 307–324. E.N. Dvorkin, A.M. Cuiti˜no, G. Gioia, Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions, Internat. J. Numer. Methods Engrg. 30 (1990) 541–564. E.N. Dvorkin, A.P. Assanelli, 2D finite elements with displacement interpolated embedded localization lines: the analysis of fracture in frictional materials, Comput. Methods Appl. Mech. Engrg. 90 (1991) 829–844. M. Klisinski, K. Runesson, S. Sture, Finite element with inner softening band, J. Eng. Mech. (ASCE) 117 (1991) 575–587. T. Belytschko, J. Fish, B.E. Engelmann, A finite element with embedded localization zones, Comput. Methods Appl. Mech. Engrg. 70 (1988) 59–89. J. Simo, S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Internat. J. Numer. Methods Engrg. 29 (1990) 1595–1638. J. Simo, J. Oliver, F. Armero, An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids, Comput. Mech. 12 (1993) 277–296. J. Simo, F. Armero, Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes, Internat. J. Numer. Methods Engrg. 33 (1992) 1413–1449. J. Oliver, Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part I: fundamentals, Internat. J. Numer. Methods Engrg. 39 (1996) 3575–3600. H. Horii, T. Ichinomiya, Observation of fracture process zone by laser speckle technique and governing mechanism in fracture of concrete, Int. J. Fract. 51 (1991) 19–29. F. Wittmann, X. Hu, Fracture process zone in cementitious materials, Int. J. Fract. 51 (1991) 3–18. K. Otsuka, H. Date, Fracture process zone in concrete tension specimen, Eng. Fract. Mech. 65 (2000) 111–131. J. Oliver, M. Cervera, O. Manzoli, Strong discontinuities and continuum plasticity models: the strong discontinuity approach, Int. J. Plast. 15 (1999) 319–351. M. Jir´asek, Comparative study on finite elements with embedded discontinuities, Comput. Methods Appl. Mech. Engrg. 188 (2000) 307–330. H. Lotfi, P. Shing, Embedded representation of fracture in concrete with mixed finite elements, Internat. J. Numer. Methods Engrg. 38 (1995) 1307–1325. T.C. Gasser, G.A. Holzapfel, Geometrically non-linear and consistently linearized embedded strong discontinuity models for 3D problems with an application to the dissection analysis of soft biological tissues, Comput. Methods Appl. Mech. Engrg. 192 (2003) 5059–5098. J. Oliver, Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part II: numerical simulation, Internat. J. Numer. Methods Engrg. 39 (1996) 3601–3623. F. Armero, K. Garikipati, An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids, Internat. J. Solids Structures 33 (1996) 2863–2885. R. Larsson, K. Runesson, Element-embedded localization band based on regularized displacement discontinuity, J. Eng. Mech. (ASCE) 122 (1996) 402–411. R. Larsson, K. Runesson, S. Sture, Embedded localization band in undrained soil based on regularized strong discontinuity—theory and FE-analysis, Internat. J. Solids Structures 33 (1996) 3081–3101. K. Washizu, Variational Methods in Elasticity and Plasticity, third ed., Pergamon Press, 1982.

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

365

[34] R.I. Borja, A finite element model for strain localization analysis of strongly discontinuous fields based on standard Galerkin approximation, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1529–1549. [35] J. Oliver, A. Huespe, S. Blanco, D. Linero, Stability and robustness issues in numerical modeling of material failure with the strong discontinuity approach, Comput. Methods Appl. Mech. Engrg. 195 (2006) 7093–7114. [36] R.I. Borja, Finite element simulation of strain localization with large deformation: capturing strong discontinuity using a Petrov–Galerkin multiscale formulation, Comput. Methods Appl. Mech. Engrg. 191 (2002) 2949–2978. [37] R. Radulovic, Numerical modeling of localized material failure by means of strong discontinuities at finite strains (Ph.D. thesis), Ruhr University, Bochum, Germany, 2010. [38] J. Mosler, G. Meschke, 3D modelling of strong discontinuities in elastoplastic solids: fixed and rotating localization formulations, Internat. J. Numer. Methods Engrg. 57 (2003) 1553–1576. [39] C. Feist, G. Hofstetter, An embedded strong discontinuity model for cracking of plain concrete, Comput. Methods Appl. Mech. Engrg. 195 (2006) 7115–7138. [40] J. Oliver, A. Huespe, Theoretical and computational issues in modelling material failure in strong discontinuity scenarios, Comput. Methods Appl. Mech. Engrg. 193 (2004) 2987–3014. [41] J. Oliver, I. Dias, A. Huespe, Crack-path field and strain-injection techniques in computational modeling of propagating material failure, Comput. Methods Appl. Mech. Engrg. 274 (2014) 289–348. [42] J. Oliver, A consistent characteristic length for smeared cracking models, Internat. J. Numer. Methods Engrg. 28 (1989) 461–474. [43] G. Camacho, M. Ortiz, Computational modelling of impact damage in brittle materials, Internat. J. Solids Structures 33 (1996) 2899–2938. [44] G. Meschke, P. Dumstorff, Energy-based modeling of cohesive and cohesionless cracks via X-FEM, Comput. Methods Appl. Mech. Engrg. 196 (2007) 2338–2357. [45] S. Mariani, U. Perego, Extended finite element method for quasi-brittle fracture, Internat. J. Numer. Methods Engrg. 58 (2003) 103–126. [46] T. Belytschko, D. Organ, C. Gerlach, Element-free Galerkin methods for dynamic fracture in concrete, Comput. Methods Appl. Mech. Engrg. 187 (2000) 385–399. [47] J. Simo, T. Hughes, Computational Inelasticity, Vol. 7, Springer, 1998. [48] Y. Zhang, Simulation methods for durability assessment of concrete structures: multifield framework and strong discontinuity embedded approach (Ph.D. thesis), Vienna University of Technology, 2013. [49] J. Oliver, A. Huespe, E. Samaniego, E. Chaves, On strategies for tracking strong discontinuities in computational failure mechanics, in: H. Mang, F. Rammerstorfer, and J. Eberhardsteiner, (Eds.) Fifth World Congress on Computational Mechanics, WCCM V, 2002. [50] M. Stolarska, D. Chopp, N. Mo¨es, T. Belytschko, Modelling crack growth by level sets in the extended finite element method, Internat. J. Numer. Methods Engrg. 51 (2001) 943–960. [51] J. Rots, Computational modeling of concrete fracture (Ph.D. thesis), Delft University of Technology, 1988. [52] B. Winkler, Traglastuntersuchungen von unbewehrten und bewehrten Betonstrukturen auf der Grundlage eines objektiven Werkstoffgesetzes f¨ur Beton (Ph.D. thesis), University of Innsbruck, 2001. [53] P. Dumstorff, G. Meschke, Crack propagation criteria in the framework of X-FEM-based structural analyses, Int. J. Numer. Anal. Methods Geomech. 31 (2007) 239–259. [54] M. Jir´asek, T. Zimmermann, Analysis of rotating crack model, J. Eng. Mech. (ASCE) 124 (1998) 842–851. [55] J. Mosler, On advanced solution strategies to overcome locking effects in strong discontinuity approaches, Internat. J. Numer. Methods Engrg. 63 (2005) 1313–1341. [56] R. Radulovic, O. Bruhns, J. Mosler, Effective 3D failure simulations by combining the advantages of embedded strong discontinuity approaches and classical interface elements, Eng. Fract. Mech. 78 (2011) 2470–2485. [57] M. Cervera, L. Pela, R. Clemente, P. Roca, A crack-tracking technique for localized damage in quasi-brittle materials, Eng. Fract. Mech. 77 (2010) 2431–2450. [58] M. Jir´asek, T. Zimmermann, Embedded crack model: II. Combination with smeared cracks, Internat. J. Numer. Methods Engrg. 50 (2001) 1291–1305. [59] Y. Theiner, G. Hofstetter, Numerical prediction of crack propagation and crack widths in concrete structures, Eng. Struct. 31 (2009) 1832–1840. [60] J. Oliver, A. Huespe, Continuum approach to material failure in strong discontinuity settings, Comput. Methods Appl. Mech. Engrg. 193 (2004) 3195–3220. [61] A. Gravouil, N. Mo¨es, T. Belytschko, Non-planar 3D crack growth by the extended finite element and level sets—part II: level set update, Internat. J. Numer. Methods Engrg. 53 (2002) 2569–2586. [62] J.F. Unger, S. Eckardt, C. K¨onke, Modelling of cohesive crack growth in concrete structures with the extended finite element method, Comput. Methods Appl. Mech. Engrg. 196 (2007) 4087–4100. [63] G. Meschke, P. Dumstorff, How does the crack know where to propagate?—A X-FEM-based study on crack propagation criteria, in: H. Yuan, F. Wittmann (Eds.), Nonlocal Modelling of Failure of Materials, 2007, pp. 201–218. Proceedings of an International Workshop (Nonlocal Modeling of Material’s Failure, NMMF 2007). [64] P. J¨ager, P. Steinmann, E. Kuhl, Modeling three-dimensional crack propagation—a comparison of crack path tracking strategies, Internat. J. Numer. Methods Engrg. 76 (2008) 1328–1352. [65] C. Feist, A numerical model for cracking of plain concrete based on the strong discontinuity approach (Ph.D. thesis), University of Innsbruck, 2004. [66] M. Xie, W. Gerstle, Energy-based cohesive crack propagation modeling, J. Eng. Mech. (ASCE) 121 (1995) 1349–1358. [67] J. Sancho, J. Planas, D. Cend´on, E. Reyes, J. G´alvez, An embedded crack model for finite element analysis of concrete fracture, Eng. Fract. Mech. 74 (2007) 75–86. [68] G. Wells, L. Sluys, Three-dimensional embedded discontinuity model for brittle fracture, Internat. J. Solids Structures 38 (2001) 897–913. [69] M. Nooru-Mohamed, Mixed-mode fracture of concrete: an experimental approach (Ph.D. thesis), Delft University of Technology, 1992.

366

Y. Zhang et al. / Comput. Methods Appl. Mech. Engrg. 287 (2015) 335–366

[70] A. Kobayashi, M. Hawkins, D. Barker, B. Liaw, Fracture process zone of concrete, in: S. Shah (Ed.), Application of Fracture Mechanics to Cementitious Composites, Martinus Nijhoff, 1985, pp. 25–50. [71] J. Oliver, A. Huespe, M. Pulido, E. Chaves, From continuum mechanics to fracture mechanics: the strong discontinuity approach, Eng. Fract. Mech. 69 (2002) 113–136. [72] M. Arrea, A. Ingraffea, Mixed-mode crack propagation in mortar and concrete, Tech. Rep., No. 81-13, Department of Structural Engineering, Cornell University, Ithaca, New York, 1982. [73] P. Bocca, A. Carpinteri, S. Valente, Mixed mode fracture of concrete, Internat. J. Solids Structures 27 (1991) 1139–1153. [74] Z. Yang, J. Chen, Fully automatic modelling of cohesive discrete crack propagation in concrete beams using local arc-length methods, Internat. J. Solids Structures 41 (2004) 801–826. [75] M. Geers, R.D. Borst, R. Peerlings, Damage and crack modeling in single-edge and double-edge notched concrete beams, Eng. Fract. Mech. 65 (2000) 247–261. [76] G. Wells, L. Sluys, A new method for modelling cohesive cracks using finite elements, Internat. J. Numer. Methods Engrg. 50 (2001) 2667–2682. [77] T.C. Gasser, G.A. Holzapfel, Modeling 3D crack propagation in unreinforced concrete using PUFEM, Comput. Methods Appl. Mech. Engrg. (2005) 2859–2896. [78] C. Feist, G. Hofstetter, Three-dimensional fracture simulations based on the SDA, Int. J. Numer. Anal. Methods Geomech. 31 (2007) 189–212. [79] J. Mosler, G. Meschke, Analysis of mode I failure in brittle materials using the strong discontinuity approach with higher order elements, in: Z. Waszczyszyn and J. Pamin, (Eds.) European Congress on Computational Mechanics, 2001.