- Email: [email protected]
Discrete Fracture Propagation Analysis using a Robust Combined Continuation Method
Journal Pre-proof
Discrete Fracture Propagation Analysis using a Robust Combined Continuation Method Cristian Mejia , Luis. F. Paullo , Deane Roehl PII: DOI: Reference:
S0020-7683(20)30036-6 https://doi.org/10.1016/j.ijsolstr.2020.02.002 SAS 10600
To appear in:
International Journal of Solids and Structures
Received date: Revised date: Accepted date:
8 July 2019 18 December 2019 1 February 2020
Please cite this article as: Cristian Mejia , Luis. F. Paullo , Deane Roehl , Discrete Fracture Propagation Analysis using a Robust Combined Continuation Method, International Journal of Solids and Structures (2020), doi: https://doi.org/10.1016/j.ijsolstr.2020.02.002
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.
1
Discrete Fracture Propagation Analysis using a Robust Combined Continuation Method Cristian Mejia a, Luis. F. Paullo a,b, Deane Roehl a,b a
Tecgraf Institute, Pontifical Catholic University of Rio de Janeiro, Rua Marquês de São Vicente, 225, Gávea, 22453-900, Rio de Janeiro, RJ, Brasil, www.tecgraf.puc-rio.br b
Department of Civil and Environmental Engineering, Pontifical Catholic University of Rio de Janeiro, Rua Marquês de São Vicente, 225, Gávea, 22453-900, Rio de Janeiro, RJ, Brasil, www.civ.puc-rio.br
Abstract. The prediction of equilibrium paths in the presence of brittle fractures remains a challenge, despite the advances in the fracture modelling techniques. Difficulties arise in the solution of the equilibrium equations in the presence of material and structural instabilities. This work proposes a continuation method with multiple restrictions to solve the nonlinear system of equilibrium equations and overcome global convergence problems owing to the sudden energy release in fracture propagation analyses. The method combines simultaneously the cylindrical arc-length and the dissipation-energy control methods. In addition, this method allows the use of algorithms for the solution of linear systems of equations, which do not require the assembly of the global stiffness matrix. Numerical examples of tensile-mode and mixed-mode fracture propagation are presented and compared with experimental results. A zero-thickness interface element with linear softening damage model simulates mixed-mode fracture. The results show the accuracy, robustness, and effectiveness of the proposed continuation method even for challenging fracture problems. Keywords: Rock fracture, Softening behavior, Cohesive elements, Damage law, Continuation methods.
2
1
Introduction The prediction of material failure is essential in many mechanical and civil engineering
problems. In the last decades, several computational methods have been developed to model fracture processes with efficiency and accuracy, and great advances have been made. Since the first works using the finite element method to model fracture, one of the main challenges has been how to simulate crack propagation with arbitrary paths. One of the first numerical methods for fracture propagation was formulated in the framework of linear elastic fracture mechanics. According to this approach, fracture growth is associated to a mesh update to include the discontinuity[1]. An alternative approach considers each element edge as a potential crack and allows crack opening/propagation when specific criteria are satisfied [2,3]. The same concept is implemented using zero-thickness interface elements with double nodes at each boundary between solid elements [4–7]. This technique has reported very realistic representations of fracture behavior in concrete (meso-mechanical simulation of concrete specimens [8]) and rock (evolution of cracking and fractures around oil wells [9]). However, it leads to larger global systems of equations, more critically in 3D cases. The embedded discontinuity approach is another technique to insert strong discontinuities into existing finite element models [10–12]. For such, suitable condensation techniques are required in order to describe the kinematics properly, to transfer stresses correctly between the cohesive zones, and to eliminate locking [13,14]. This approach requires crack-tracking algorithms to allow crack growth through the finite element mesh. More recently, an extended finite element method has been widely used to model strong discontinuities. Initially proposed by Moës and Belytschko [15], this method does not require the definition of potential cracks and propagation paths a priori. A fracture path-tracing algorithm identifies and inserts the crack path in the model during the analysis. The level set method and the fast marching method have been developed as efficient techniques to track the crack path [16]. The extended finite element method is an attractive technique to model discontinuity. A cohesive zone model simulates the fracture process zone ahead of the crack tip, where material softening takes place [17]. A cohesive constitutive model provides the relation between the cohesive traction on the crack surface and the separation across the crack surfaces. Good overviews and discussions related to cohesive constitutive models are given, for example, in Papoulia et al. [18], Park et al. [17] and Nguyen [7]). The cohesive interface element method, which has its origin in the concept of the cohesive zone model proposed by Barenblatt [19] and Dugdale [20], has been used with great success
3
to simulate fracture processes in concrete, metals, ceramics, composites and rocks. Several research studies such as Hillerborg et al. [21], Xu and Needleman [22] and Camacho and Ortiz [2] have focused on the development of the cohesive zone concept in computational fracture mechanics. In addition, the concept of cohesive zone model has also been applied in conjunction with the extended finite element method XFEM [23]. On the other hand, numerical modeling of materials with softening behavior is prone to convergence problems owing to sudden release of energy in the post-critical regime. In order to overcome this problem, some researchers such as Chaboche et al. [24], Hamitouche et al. [25], Kayhan [26], Jiang [27] and Lapczyk et al. [28] have adopted viscous regularization techniques in discrete fracture propagation. These techniques modify the constitutive equation of the interface element using a viscous damage variable. However, they lead to stresses that drift from the traction-separation law [29,30]. In this context, the numerical analysis of structural problems with strong material softening behavior requires the use of control methods capable of following the complex loading/unloading history of the structure introduced by strong stiffness degradation. Generally, in non-linear geometric and non linear material problems the so-called continuation methods are used. Among them, one of the most used family of continuation methods are arc-length control methods [31–34]. The main pros of the arc-length methods are versatility and easy computational implementation. Nevertheless, many researchers have reported difficulties with the use of the arc-length control methods when analyzing strong material softening problems, as pointed out in Pohl et al. [35] and Paullo and Roehl [36]. The indirect strain increment method controls the monotonic increment of strain at a certain point of interest in the structure and improves convergence of the solution for some softening problems [37]. However, fracture propagation often exhibits non-monotonic strain increase, especially in the region near the crack opening. Methods based on energy dissipation also have good performance in material softening problems [38,39]. These methods require the use of an additional control at non-dissipative steps. Thus, there is the need to define a switch parameter based on a minimum amount of released energy. The definition of an adequate parameter value requires good knowledge of the system behavior. Labanda et al. [40] introduced a restriction based on energy release criterion in an augmented Lagrangian formulation to model cohesive cracks. To overcome convergence difficulties associated with material softening behavior and contact/friction interfaces, Oliver et al. proposed an implicit-explicit integration scheme, which combines implicit integration of stresses in the constitutive model with explicit
4
extrapolation of the internal variables [41]. Recently, incremental sequentially linear analysis (ISLA) has been proposed by Yu et al. [42], which combines the classical Newton-Raphson method with Sequentially Linear Analysis. In this method, each increment starts and ends with an equilibrium state. Equilibrium is updated based on the stiffness reduction of critical elements and load increments. The ISLA method demonstrated robustness to simulate the failure of quasi-brittle material under non-proportional loading. Paullo and Roehl [36] formulated a technique to combine two or more restrictions through the least squares methods, obtaining a mixed control. The mixed control showed good performance in problems with strong geometric nonlinearities and elastoplastic behavior. In addition, the combination reduces the participation of the user in the control of the analysis. This work proposes a step forward of the continuation method with combined constraints formulated by Paullo and Roehl [36] by simplifying the control equation. This simplification allows the use of the combined approach in association with direct solution strategies. In addition, the alternative strategy reduces the computational effort in the solution of the linear system of equations. To deal with strong dissipation in fracture problems the method combines the cylindrical arc-length control and the dissipated energy control. Numerical examples of fracture problems modeled with interface cohesive elements show the robustness and effectiveness of the proposed strategy. Conventional continuation methods cannot handle these examples. 2
Cohesive Interface Element
2.1 Interface element formulation The weak form of the equilibrium equation, including the interface effects, follows from the principle of virtual work. The virtual work done by the body forces ( ) in the domain and by the external traction ( energy in the domain
) on the boundary
are equal to the sum of the virtual strain
and the cohesive fracture energy on the fracture surface ( ) (see
Figure 1). ∫ where
∫ and
∫
∫
Eq. 1
are virtual strains, virtual displacements and virtual fracture opening,
respectively. In addition,
is the stress tensor in the continuum element while
cohesive traction along the fracture surface.
is the
5
𝑇𝑒𝑥
𝛤
𝛺 𝛤𝑓
𝑇𝑐
Figure 1 Schematic representation of a domain ( ) with an internal fracture ( ) under external tractions ( and cohesive tractions ( ).
)
The domain ( ) is discretized into continuum elements and zero-thickness interface elements along the crack path. The displacement field ( ̅) over the domain ( ) is interpolated from the nodal displacement vector ( ) by using shape functions ( ). ̅
Eq. 2 The local separation along the fracture ( ) is approximated from the nodal displacements
by using the global displacement-separation matrix (
). Eq. 3
where
is given by Eq. 4
is the shape function matrix,
is the local displacement-separation matrix and
is the
rotation matrix. Following the standard finite element procedures [43,44], the internal force vector (
) of an interface element is given by
∫
Eq. 5
The gradient of the internal force vector leads to the tangent matrix of the interface element, ∫ where
Eq. 6 is the constitutive matrix relating tractions and relative displacements
on the fracture surface.
6
2.2 Cohesive zone model The cohesive constitutive model can be defined directly in terms of the traction separation law. The traction on the cohesive element in the fracture process zone consists of the normal and the shear
components, which are related to the shear
and normal
separation through the constitutive equations [44]. , Here
[
[ ]{
}
Eq. 7
] is the tangential constitutive matrix of the cohesive element,
the shear elastic stiffness and
is
is the normal elastic stiffness.
This work adopts a cohesive model, which is characterized by a linear softening law. This relationship assumes linear elastic behavior up to the critical point of damage initiation, where the traction reaches the maximum cohesive traction
matching the damage initiation
criterion. The initiation of the softening process follows a quadratic failure criterion, which considers the interaction of the traction components according to [45] ( where
)
is the normal strength,
operator, defined as 〈 〉
| |
.
〈 〉
Eq. 8
/
is the shear strength and 〈 〉 is the Macaulay bracket .
The effective separation of damage initiation cohesive traction
√
〈 〉 . Here
and
corresponds to the point of maximum represent the normal and shear stress
that meet the damage initiation criterion (Eq. 8). Thus, this maximum cohesive strength considers indirectly the effect of normal and shear stress components. Softening behavior occurs in the post-critical zone, describing progressive degradation of the material stiffness. Figure 2 shows the damage model with linear softening [46,47].
7
Point of damage initiation
Traction, 𝜏
𝜏
𝑑𝜏
𝜏𝑜 𝑑 𝜏 𝑘𝑖
Post critical regime
𝐺𝑐
𝑑 𝑘𝑖 Δ𝑓 Δ𝑜 Effective displacement jump Δ
𝑘𝑖
Δ
Figure 2 Damage model with linear softening.
The area under the traction-separation curve represents the cohesive fracture energy . In the post-critical regime, the cohesive traction
̅
̅
can be defined as,
̅ ̅
Eq. 9
{
where d is the damage variable. The damage variable for linear softening can be defined as proposed in [29,48], Eq. 10
where
is the effective separation at damage initiation, √〈
for material failure,
〉
consistent tangent modulus
is the critical effective separation
is the effective displacement. Therefore, the
is obtained by differentiating the cohesive tractions with
respect to the fracture separation [49]: ̇ ̇
{
Eq. 11 .
̅
〈
〉
/.
where ̅ is the Kronecker delta,
̅
〈
〉
/
is the elastic stiffness matrix, and
Eq. 12
is a scalar value
defined as: Eq. 13
8
The consistent tangent modulus for unloading and reloading is a function of the damage variable and is given by, [
Eq. 14
]
Finally, a crack develops when the effective separation reaches the critical value of
and
the cohesive tractions become zero. Other cohesive models can be adopted by changing the damage evolution law given by Eq. 10, such as exponential softening model. 3
Solution of the nonlinear system of equations
3.1 Classical Newton-Raphson solution In a finite element context, the incremental global equilibrium can be defined by the following system of equations Eq. 15 where
is the global nonlinear internal force vector which depends on the global
displacement vector parameter, and
.
The subscript
indicates the -th analysis step,
is the load
is the reference nodal force vector. The increment of the displacements and
of the load factor for step
can be obtained as. Eq. 16
Accordingly, the incremental equilibrium in step
is obtained as Eq. 17
The solution of Eq. 17 requires a restriction in
, which is expressed as Eq. 18
As indicated in Paullo and Roehl [36], an iterative Newton-Raphson scheme is applied to the extended system defining the following 0
where
1{
}
{
is the tangent stiffness matrix,
iteration index, and
Eq. 19
}
and
is the unbalanced force vector,
is the
are the displacement vector and load factor correctors,
respectively. Following the proposal in Paullo and Roehl [36], the following error measures
9
are adopted. |
|
|,
Eq. 20
|
3.2 Cylindrical Arc-Length control method For nonlinear structural problems the cylindrical arc-length equation [31] is frequently used as restriction (Eq. 18) in the Newton Raphson scheme. In this case, the restriction is defined by (
) |
where
Eq. 21
|
is the arc length control parameter. The corrector of the load parameter,
,
follows from the solution of the quadratic equation: Eq. 22 In this equation the coefficients are: |(
)
(
)
(
)
| (
) (
Eq. 23
)
Definition of the increment size is performed by an adaptation procedure introduced by Ramm [34] and also used by Chaisomphob [50]. This scheme defines the parameter
for a
subsequent step as √ where
is the desired number of iterations and
Eq. 24 the number of iterations required to
achieve convergence. Paullo and Roehl [36] point out the need of for higher and lower limits for
in order to achieve adequate control parameter values throughout the incremental
analysis. 3.3
Dissipated energy increment control method
For problems with energy release by damage, the equation of restriction based on the increment of dissipated energy can be obtained following the formulation of Gutierrez [38]. In the geometrically linear case, we have
10
( where
Eq. 25
)
is the desired incremental dissipated energy. We obtain the load parameter
corrector by introducing Eq. 25 in the extended equilibrium system in Eq. 19, (
)
(
)
Note that
and
Eq. 26
correspond to the converged load factor and displacement vector
from the previous step, respectively. It is important to point out that for elastic incremental steps, where no dissipation occurs, this control fails and alternative ones must be adopted. 3.4
Full combined restriction
Aiming at providing a more robust solution strategy in problems with both material softening and geometric nonlinearity, Paullo Muñoz e Roehl [36] formulate a combination of two or more restrictions working simultaneously. According to [36], the general non-linear system of equations defined in Eq. 17 to Eq. 19 can include
restrictions grouped in a
-
dimensional vector defined by: ̃
where ̃
̃
{
[ ̃ ̃
̃
Eq. 27
}
̃ ], is the vector containing the control parameters. Introducing the
vector of restrictions in the Newton Raphson scheme results in the following inconsistent system of equations. 0
1{
}
{
Eq. 28
}
Paullo and Roehl [36] applied the least squares technique to obtain an approximate solution with the participation of all restrictions. Thus, the corrector of displacements and the load factor follow directly from {
}
(0
1 0
1)
0
1 {
}
Eq. 29
Control parameter values at each iteration are recalculated using the corresponding single
11
constraint equation to guarantee convergence of the solution. The converged values of the control parameters are used to estimate the predictor of the following step. In addition, each control parameter is adapted according to Eq. 24. Convergence rate and robustness of the combined constraint technique are studied in [36]. From Eq. 29 the load factor corrector ( ((
)
((
)
follows as
)
)[
)[
]( (
)
+
] ((
)
)
(
)
(
)
)
Eq. 30
)
where *(
)
(
Eq. 31
Consequently, the displacement vector corrector is * (
)
(
)
((
)
(
)
)
+
Eq. 32
3.5 Simplified combined restriction A combined restriction technique requires a full coefficient matrix
for the solution of
the linear system of equations. Consequently, the possibility of employing algorithms that take advantage of sparsity of the stiffness matrix is lost, increasing strongly the computational effort. Moreover, some solution schemes do not assemble the global stiffness matrix. In this context, a modified constraint combination technique is proposed to allow working with general solution methods. We can define: (
)
( were
Eq. 33
) and
are the displacement solutions corresponding to the unbalanced forces and
the reference forces, respectively. The displacement corrector derives from Eq. 28 to Eq. 34 Substitution Eq. 34 in the second line of Eq. 28, establishes the following inconsistent system of equations {
}
{
Eq. 35
}
Applying the least squares method to solve Eq. 35 gives the load parameter corrector {
} {
}
{
} {
}
Eq. 36
12
It is important to point out that when none of the constraint equations depend on the load parameter, i.e.,
, Eq. 30 and Eq. 35 provide the same result. If we consider the
combination of the cylindrical arc-length control and the dissipated energy control, we add the following restrictions to the iterative procedure 3,
2
2
Eq. 37
3
Substituting Eq. (34) in Eq. (33) we obtain (
)
( (
(
) (
) )
(
) (
(
))
)
( (
)
( )
(
)
(
(
)
Eq. 38
))
Alternatively, we can write the following expression for the predictor of the load parameter, which considers the positive work in the arc-length predictor contribution: ( ((
)
) )
( (
(
)
Eq. 39
))
The proposed combination technique does not require the assembly of the full coefficient matrix, reduces the computational effort, and allows using direct solvers. In addition, the proposed method combines the arc-length and the dissipated energy controls. The values of and
in the first step and first iteration are calculated from Eq. 22 and Eq. 25,
respectively, using an initial prescribed load parameter
. It is important to point out that Eq.
38 works for steps with and without dissipation. Thus, the use of a switch procedure is not necessary. For non-dissipative steps, Eq. 38 reduces to a linearized cylindrical arc-length equation. 4
Numerical Examples The zero-thickness interface element and the solution strategies proposed are implemented
in the in-house multiphysics framework called GeMA (Geo Modelling Analysis) [51]. Three examples show the robustness and effectiveness of the proposed strategy. 4.1 L-Shaped Plate (Single Fracture Trajectory) This section assesses fracture propagation in an L-Shaped Plate. Figure 3 presents the structure dimensions and loading. This example illustrates the case of a single fracture path predefined by a line of cohesive elements. The configuration and properties of the problem are taken from Pohl et al. [35]. Plane strain four-node quadrilateral elements model the
13
structure. The solid elements have linear elastic behavior with Young’s modulus E = 100kN/mm2 and Poisson’s coefficient
. The line of cohesive elements is composed by
four-node interface elements. The cohesive elements have the following constitutive kN/mm3, shear elastic stiffness
parameters: normal elastic stiffness kN/mm3, normal strength
N/mm2, shear strength
effective total displacement at failure al. [35]. The elastic stiffness values
N/mm2 and the
mm, in accordance with Pohl et and
are computed by the following relation
. In order to reduce the effect of this artificial stiffness on the elastic response, we adopt thickness
= 1 mm [52].
5 𝑚
cohesive elements
5 𝑚
λ×
5 𝑚
𝑘𝑁 𝑚
x
5 𝑚
Figure 3 Geometry and FEM mesh of L–shaped plate.
Figure 4 shows the equilibrium path for the horizontal displacement
obtained with the
cylindrical arc-length control and with the combined control methods. The cylindrical arclength control fails in the region near the maximum load factor showing its difficulty to analyze fracture propagation problems characterized by material softening behavior. The works Paullo and Roehl [36] and Pohl et al. [35] present similar problems related to the use of the arc-length control methods in structures with material non linearity. In those cases, the arc-length presented convergence problems in analyses involving plasticity with linear softening and damage. On the other hand, the proposed combined method showed good performance and was able to carry out the analysis up to the stage of total opening of all
14
cohesive elements (total crack propagation). Figure 4 also shows a very complex equilibrium path exhibiting local loading and unloading regions in which each “finger” corresponds to fracture propagation of a single cohesive element. These fingers (local snap-backs) lead to failure of traditional solution control algorithms. However, in this example, finer meshes render smoother solutions, allowing the use of standard arc-length method. Figure 5 shows the deformed shape of the structure when 9 elements reach the limit damage factor d =1. It is noticeable that the solid elements do not exhibit important geometric distortions. In addition, the main mechanism of fracture is normal traction. Cylindrical arc-length
Figure 4 Load factor
vs. horizontal displacement Ux1 considering
Proposed method
and
5.
Figure 5 Final deformed configuration where nine cohesive elements collapsed, scale factor = 5.
4.2 Plate with eccentric hole This section presents fracture propagation through a square plate with an eccentric hole
15
under uniaxial stress.
Figure 6 presents the plate dimensions, boundary and loading
conditions. The aim of this problem is to assess the performance of the proposed algorithm in snap-back situations. The geometry and properties of the problem are from Labanda et al. [40]. Four-node quadrilateral elements discretize the structure under plane stress condition. The solid elements have linear elastic behavior with Young’s modulus E = 20 GPa and Poisson’s coefficient
. The fracture behavior is modeled using four-node zero
thickness interface elements placed horizontally, as shows Figure 6. The cohesive elements have the following constitutive parameters: normal elastic stiffness elastic stiffness
GPa/mm, normal strength
GPa/mm, shear
MPa, shear strength
MPa . The effective total displacement at failure is corresponds to a fracture energy release of
mm. This value N/mm, as adopted in Labanda et al.[40]. λ×
𝑘𝑁 𝑚
5𝑚 P Cohesive elements
𝑚
5𝑚
5𝑚
𝑚
5𝑚
Figure 6 Geometry and FEM mesh of plate with hole.
Figure 7 shows the equilibrium path in terms of the load factor versus the vertical displacement of point P. This example was simulated using 200 load steps. The proposed combined method is able to track complex equilibrium paths including snap back phenomena. Figure 7 also shows that the equilibrium path obtained matches the numerical results reported by Labanda et al. [40] that were obtained with a switched energy dissipated control. Figure 8 shows four stages, A, B, C, D, according to the damage progress. Up to state A, the peak load, the fracture propagates in both directions of the hole (see Figure 8a). Then, the fracture develops on the right side of the plate producing unloading behavior until point B,
16
when the fracture reached the plate border. Between B and C fracture aperture increases. After point C, the fracture propagates through the left side of the plate until point D, when the left fracture reached the plate border (see Figure 8d).
A
C B
D
Figure 7 Load factor vs Vertical displacement in P. A
B Initiation
D C
Figure 8 Damage evolution for the plate with hole under uniaxial stress: states A, B, C and D.
Figure 9 shows the damage evolution along the interface elements placed horizontally for the states presented in Figure 8 (A, B, C and D).
17
Left side
Hole
Right side
Left side
(a)
Hole
Right side
(b)
Figure 9 Damage evolution along the fracture path for a) state A and B, and b) state C and D.
4.3 Perforated cantilever beam: Single crack trajectory and multiple snap-back path This example presents the delamination process in a perforated cantilever beam that has been extensively studied in the literature on continuation methods for fracture and softening problems. Initially studied by Verhoosel et al. [39], this structure attracts the attention of researchers, such as Pohl et al. [35], Labanda et al. [40] and Paullo and Roehl [36], once the load-displacement path exhibits several snap-back phenomena along the delamination process. Verhoosel studied this problem using an energy release control switched with load control for non-dissipative steps. Figure 10 shows the geometry, mesh, boundaries and loading conditions. The geometry and properties of the structure are from Paullo and Roehl [36] and Verhoosel et al. [39], respectively. A line of linear cohesive elements defines the single crack trajectory, with the following properties: GPa/mm,
MPa,
GPa/mm,
MPa The effective total displacement at failure is
5 mm that corresponds to a fracture energy release of
5
N/mm. Plane strain four-node quadrilateral elements discretize the linear elastic solid domain with Young’s modulus thickness of
mm.
MPa and Poisson’s coefficient
. The model has a
18
Figure 10 Geometry and finite element model.
Figure 11 shows the vertical displacement at the top of the free end of the beam vs the load factor and the deformed shape of the structure after 150 converged steps. That stage corresponds to the total collapse of five segments. As expected, the equilibrium path exhibits several load and displacement limit points, as well as several snap-back forming a “finger” as each segment collapses. This behavior agrees with the response in Verhoosel et al. [39] and shows the robustness of the proposed continuation method to trace the equilibrium path with several snap-back situations. (a)
Figure 11 Load factor
(b)
vs. vertical displacement and b) deformed shape after 150 steps with the proposed method. and 5.
4.4 Four-Point Bending Test (Multiple fracture trajectory) In this section, crack propagation in a four-point bending beam test is assessed. The geometry and properties of the system are from Gálvez et al. [53], which shows experimental results mainly for two kinds of support condition. Figure 12 shows the load and support conditions for the two models presented in [53]. The beam models have a thickness of 50 mm. The main difference between the two configurations is the support condition. This
19
difference targets at assessing different fracture mechanisms since in configuration Type-I it is expected that fracture occurs due to normal tractions, while in configuration Type-II, fracture is due to mixed normal and shear tractions, as indicated in Gálvez et al. [53]. 𝐷
5𝑚𝑚 𝑃
𝐷
𝑇𝑌𝑃𝐸
𝐼
𝐷
𝐶𝑜 𝑒𝑠𝑖𝑣𝑒 𝑟𝑒𝑔𝑖𝑜𝑛
𝑃
𝐷
𝑇𝑌𝑃𝐸
5𝑚𝑚
𝐼𝐼
𝐶𝑜 𝑒𝑠𝑖𝑣𝑒 𝑟𝑒𝑔𝑖𝑜𝑛
5𝑚𝑚 𝐷
𝐷 𝐷
𝐷
𝐷
𝐷
𝐷
𝐷
𝐷
𝐷
𝐷
5𝐷
𝐷
𝐷
𝐷
Figure 12 Geometry and properties of beams Type - I and Type – II, taken from Galvez [46].
The material behavior of the solid elements is linear elastic while the cohesive elements present linear softening behavior. Since in the experiments fractures develop in the central region of the beam, the finite element discretization includes cohesive interface elements surrounding 3-node triangular solid elements in this region. The use of triangular elements in this region aims at allowing a more alternatives for fracture paths. Outside this region, discretization follows with 4-node quadrilateral solid elements. Figure 13 shows the finite element mesh of the beam model where the bounding box identifies the region with a dense mesh of triangles with cohesive elements at its sides. In this way, the crack trajectory is not predefined. Actually, several trajectories appear during the analysis.
𝑃
𝐶𝑜 𝑒𝑠𝑖𝑣𝑒 𝑒𝑙𝑒𝑚𝑒𝑛𝑡 𝑚𝑒𝑠
Figure 13 FEM mesh for beams Type - I and Type - II
In accordance with the properties given in Galvez [46], the material has a Young’s modulus of E=38 GPa and a Poisson’s ratio of following parameters:
GPa/m,
,. The cohesive model takes the GPa/m,
MPa,
MPa
20
and
mm. Figure 14 shows the deformed shape and the accumulated scalar
damage parameter in the cohesive elements for model Type-I after the crack propagates through the cohesive region (200 converged steps). This stage corresponds to a crack-mouth opening displacement (CMOD) equal to 0.42mm. The main crack trajectory is apparent. In addition, the distribution of accumulated damage indicates that several secondary fracture trajectories exist around the main one. Figure 15 shows that the numerical main crack trajectory (solid circled line) obtained numerically matches the experimental envelope (shaded area) presented in Galvez et al. [53], indicating good performance of the cohesive crack approach to represent real fracture processes. Figure 15 also shows the numerical damage envelope for accumulated damage greater than 0.85 (dashed line) which contains the secondary trajectories. The secondary trajectories also coincide with the experimental envelope.
𝐴𝑐𝑐𝑢𝑚𝑢𝑙𝑎𝑡𝑒𝑑 𝑑𝑎𝑚𝑎𝑔𝑒
Figure 14 Deformed shape and accumulated damage factor for beam Type – I for CMOD = 0.42mm. Scale factor = 20.
21
Figure 15 Main numerical crack trajectory and damaged region (D>=0,85) for beam Type – I for CMOD = 0.42mm.
Figure 16 shows the relation between the crack-mouth opening displacement (CMOD) and load P. In Figure 16, the maximum load level agrees with the experimental envelope. In addition, the curve P vs CMOD follows coherently the experimental results. It is important to point out that these numerical results were obtained using the proposed combined method, since the cylindrical arc-length control failed in the beginning of crack propagation.
CMOD
Figure 16 Vertical Load P vs. crack mouth opening displacement (CMOD) for beam Type – I.
For model Type-II Figure 17 presents the deformed shape and the accumulated damage at
22
the Gauss integration points of the cohesive elements after final crack propagation (380 converged steps and CMOD = 0.13mm). From beam Type-I to Type-II the change in the support conditions lead to a modification of the main crack trajectory. In beam Type II damage in the cohesive elements is generated by the action of normal and tangential tractions, while in beam Type I fracture occurs mainly due to normal tractions. In addition, beam Type II presents a second important damaged zone located at the upper-left part of the cohesive region. In that case, during the fracture process, the normal traction over this zone increases. However, damage does not reach the maximum value indicating that the interface elements are in the softening process but do not reach total degradation.
𝐴𝑐𝑐𝑢𝑚𝑢𝑙𝑎𝑡𝑒𝑑 𝑑𝑎𝑚𝑎𝑔𝑒
Figure 17 Deformed shape and accumulated damage for beam Type – II for CMOD = 0.13. Scale Factor = 20.
23
Figure 18 Main crack trajectory and damaged region for beam Type – II for CMOD = 0.13.
Figure 18 shows the agreement between the numerical and the experimental results. The numerical main trajectory and the numerical damaged region around it coincide with the experimental envelope for beam type - II. Figure 19 also reflects the good agreement between the experimental and numerical results in terms of the P vs. CMOD curves. Another important aspect refers to the computational effort to reach the maximum crack propagation. Despite the complexity in the analysis, the proposed simplified combined method showed good performance, controlling the analysis until total crack propagation without interruptions in the nonlinear solution procedure. In this way, the proposed continuation method shows robustness for this kind of problem.
24
CMOD
Figure 19 Load vs. CMOD for beam Type – II.
4.5 Highly pre-stressed beam test This problem shows the robustness of the proposed combined method to simulate a prestressed beam under non-proportional loading. Figure 20 shows the beam geometry, boundary conditions and loading conditions. The geometry and properties of the beam are from Yu et al. [42]. Four-node quadrilateral elements discretize the structure under plane stress condition. Fracture behavior is modeled using four-node zero thickness interface elements, which are placed at the edge of all solid elements. The properties are Young’s modulus E = 32 GPa, Poisson’s coefficient fracture energy release of
5𝑚
, tensile strength
MPa, and
N/mm, as adopted in Yu et al. [42]. λ
5𝑚
λ
5𝑚
A P P
5𝑚
Cohesive elements at all edges
5𝑚
Figure 20 Geometry and FEM mesh of pre-stressed beam.
5𝑚
𝑚
25
The beam is subjected to the combination of horizontal load P that is contant during the simulation and vertical load
, as illustrated in Figure 20. Two cases are studied considering
P of 1 or 5 MPa. The proposed algorithm controls the vertical load, which increases up to failure and reduces afterwards. Figure 21 indicates the final crack path at the end of the simulation for both cases (prestressed beam with P =1 and P = 5 MPa). For both cases, multiple cracks are generated at the base of the beam, and after that, two main cracks propagate to the top.
a)
Multiple cracks b)
Figure 21 Prestressed beam: damage evolution and final crack pattern for a) P=1 MPa and b) P=5 MPa prestress.
The damages regions obtained are similar to those reported by Yu et al.[42]. Figure 22 shows the reaction force and vertical displacement in A ( midspan displacement) for both cases. These results match with those obtained by Yu et al.[42]. However, for 1 MPa, the curves are different because the propagation of multiple cracks produces different damage regions.
a)
b)
Figure 22 Prestressed beam: reaction force vs midspan displacement for a) 1 MPa and b) 5 MPa prestress.
26
Figure 23 shows the total number of iterations required in each load step for both cases. For the beam under a prestress of 1 MPa, the total number of steps is 100 with an average of 4 iterations per step. For the beam under prestress of 5 MPa, the curve is simulated in 200 load steps with an average of 6 iterations per step.
a)
b) 3rd attempt
3rd
2nd attempt
2nd
1st attempt
1rd
Figure 23 Number of iterations required in each load step for a) 1 MPa and b) 5 MPa prestress.
Figure 24 presents the evolution along the analysis of the arc-length increments and of the energy increments that satisfy the combined method. The contribution of both parameters ( ,
) is present throughout the analysis.
Figure 24 Evolution of the control parameter
5
and
in the combined method,
5 and
.
Conclusions and Remarks This paper presents a numerical model for fracture propagation using zero-thickness
interface elements with a linear softening cohesive model. This approach in association with a combined restriction control method is proposed to overcome convergence problems due to strong energy dissipation by fracture propagation.
The chosen restrictions were the
cylindrical arc-length and the dissipation energy. The numerical models were carried out
27
considering intrinsic cohesive zone models with very stiff initial elastic slope in order to reduce the effect of artificial compliance. An L-shaped plate with a predefined fracture path validates the solution techniques. The proposed continuation method with the combination of restrictions overcomes the convergence problems reported in [30] owing to the very complex equilibrium path generated by the local loading/unloading at crack tip. The square plate with an eccentric hole and the delamination of a perforated cantilever beam show the capability of the proposed combined method to deal with complex equilibrium paths, including multiple snap-back situations. In a third problem, the experimental results published by Galvez et al. (1998) validate both tensile mode and mixed-mode fracture in concrete under non-proportional loading. This example is a challenging benchmark for numerical fracture models considering multiple fracture paths. For this example, a target region is modelled with a zero-thickness interface element at the border of the continuum element and several fractures appear during the analysis. Furthermore, good matching of results for load-CMOD curves, ultimate-load capacity and fracture path within the experimental envelope limits highlight the strengths of the proposed strategy. The highly prestressed beam demonstrates the robustness of the proposed combined control for non-proportional loading with multiple competing cracks. According to the numerical results, the combined control overcomes the global convergence problem associated with the damage constitutive model, which also includes instabilities related to the derivative of the discontinuous Macaulay function. Finally, the combined control shows good performance to analyze the fracture propagation process in all the examples assessed. The combined restriction method is less dependent on parameter settings. The automatic step size definition also contributes in this sense. In addition, the proposed technique allows the use of optimized solution methods for systems of linear equations. REFERENCES [1]
[2] [3]
Ingraffea AR, Saouma V. Numerical modeling of discrete crack propagation in reinforced and plain concrete. Fract. Mech. Concr. Struct. Appl. Numer. Calc., Dordrecht: Springer Netherlands; 1985, p. 171–225. doi:10.1007/978-94-009-61524_4. Camacho GT, Ortiz M. Computational modelling of impact damage in brittle materials. Int J Solids Struct 1996;33:2899–938. doi:10.1016/0020-7683(95)00255-3. Ciancio D, Carol I, Cuomo M. On inter-element forces in the FEM-displacement formulation, and implications for stress recovery. Int J Numer Methods Eng 2006;66:502–28. doi:10.1002/nme.1564.
28
[4] [5] [6]
[7]
[8]
[9] [10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18] [19] [20] [21]
Rots JG, Schellekens JC. Interface elements in concrete mechanics. Comput Anal Des Concr Struct 1990;2:909–18. López CM, Carol I, Aguado A. Microstructural analysis of concrete fracture using interface elements. 2000. Carol I, López CM, Roa O. Micromechanical analysis of quasi brittle materials using fracture based interface elements. Int J Numer Methods Eng 2001;52:193–215. doi:10.1002/nme.277. Nguyen VP. Discontinuous Galerkin/extrinsic cohesive zone modeling: Implementation caveats and applications in computational fracture mechanics. Eng Fract Mech 2014;128:37–68. doi:10.1016/j.engfracmech.2014.07.003. Ciancio D, Lopez CM, Carol I, Cuomo M. New results in meso-mechanical modelling of concrete using fracture-based zero-thickness interface elements. Comput. Model. Concr. Struct., 2003, p. 171–7. Garolera D, Carol I, Papanastasiou P. Micromechanical analysis of sand production. Int J Numer Anal Methods Geomech 2019;43:1207–29. doi:10.1002/nag.2892. Dvorkin EN, Cuitiño AM, Gioia G. Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions. Int J Numer Methods Eng 1990;30:541–64. doi:10.1002/nme.1620300311. Manzoli OL, Shing PB. A general technique to embed non-uniform discontinuities into standard solid finite elements. Comput Struct 2006;84:742–57. doi:10.1016/J.COMPSTRUC.2005.10.009. Radulovic R, Bruhns OT, Mosler J. Effective 3D failure simulations by combining the advantages of embedded Strong Discontinuity Approaches and classical interface elements. Eng Fract Mech 2011;78:2470–85. doi:10.1016/J.ENGFRACMECH.2011.06.007. Jirásek M. Comparative study on finite elements with embedded discontinuities. Comput Methods Appl Mech Eng 2000;188:307–30. doi:10.1016/S00457825(99)00154-1. Zhang Y, Lackner R, Zeiml M, Mang HA. Strong discontinuity embedded approach with standard SOS formulation: Element formulation, energy-based crack-tracking strategy, and validations. Comput Methods Appl Mech Eng 2015;287:335–66. doi:10.1016/J.CMA.2015.02.001. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Methods Eng 1999;46:131–50. doi:10.1002/(SICI)10970207(19990910)46:1<131::AID-NME726>3.0.CO;2-J. Fries T-P, Baydoun M. Crack propagation with the extended finite element method and a hybrid explicit-implicit crack description. Int J Numer Methods Eng 2012;89:1527–58. doi:10.1002/nme.3299. Park K, Paulino GH, Roesler JR. A unified potential-based cohesive model of mixedmode fracture. J Mech Phys Solids 2009;57:891–908. doi:10.1016/J.JMPS.2008.10.003. Papoulia KD, Sam C-H, Vavasis SA. Time continuity in cohesive finite element modeling. Int J Numer Methods Eng 2003;58:679–701. doi:10.1002/nme.778. Barenblatt GI. The Mathematical Theory of Equilibrium Cracks in Brittle Fracture. Adv Appl Mech 1962;7:55–129. doi:10.1016/S0065-2156(08)70121-2. Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8:100– 4. doi:10.1016/0022-5096(60)90013-2. 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. Cem Concr Res 1976;6:773–81. doi:10.1016/0008-8846(76)90007-7.
29
[22] Xu X-P, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 1994;42:1397–434. doi:10.1016/0022-5096(94)90003-5. [23] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Eng Fract Mech 2002;69:813–33. doi:10.1016/S0013-7944(01)00128-X. [24] Chaboche JL, Feyel F, Monerie Y. Interface debonding models: a viscous regularization with a limited rate dependency. Int J Solids Struct 2001;38:3127–60. doi:10.1016/S0020-7683(00)00053-6. [25] Hamitouche L, Tarfaoui M, Vautrin A. An interface debonding law subject to viscous regularization for avoiding instability: Application to the delamination problems. Eng Fract Mech 2008;75:3084–100. doi:10.1016/J.ENGFRACMECH.2007.12.014. [26] Kayhan E. Application of viscous and non-local integral type regularization schemes for softening plasticity. Atilim University, 2009. [27] Jiang H. Cohesive zone model for carbon nanotube adhesive simulation and fracture/fatigue crack growth. University of Akron, 2010. [28] Lapczyk I, Hurtado JA. Progressive damage modeling in fiber-reinforced materials. Compos Part A Appl Sci Manuf 2007;38:2333–41. doi:10.1016/J.COMPOSITESA.2007.01.017. [29] Dávila CG, Camanho PP, Turon A. Cohesive elements for shells. NASA/Technical Publ 2007. [30] Simonovski I. Cohesive zone modeling of intergranular cracking in polycrystalline aggregates. Nucl Eng Des 2015;283:139–47. doi:10.1016/J.NUCENGDES.2014.09.041. [31] Riks E. An incremental approach to the solution of snapping and buckling problems. Int J Solids Struct 1979;15:529–51. doi:10.1016/0020-7683(79)90081-7. [32] Batoz J-L, Dhatt G. Incremental displacement algorithms for nonlinear problems. Int J Numer Methods Eng 1979;14:1262–7. doi:10.1002/nme.1620140811. [33] Crisfield MA. A fast incremental/iterative solution procedure that handles “snapthrough.” Comput Struct 1981;13:55–62. doi:10.1016/0045-7949(81)90108-5. [34] Ramm E. Strategies for Tracing the Nonlinear Response Near Limit Points. Nonlinear Finite Elem. Anal. Struct. Mech., Berlin, Heidelberg: Springer Berlin Heidelberg; 1981, p. 63–89. doi:10.1007/978-3-642-81589-8_5. [35] Pohl T, Ramm E, Bischoff M. Adaptive path following schemes for problems with softening. Finite Elem Anal Des 2014;86:12–22. doi:10.1016/J.FINEL.2014.02.005. [36] Paullo Muñoz LF, Roehl D. A Continuation method with combined restrictions for nonlinear structure analysis. Finite Elem Anal Des 2017;130:53–64. doi:10.1016/J.FINEL.2017.02.003. [37] Chen Z, Schreyer HL. A numerical solution scheme for softening problems involving total strain control. Comput Struct 1990;37:1043–50. doi:10.1016/00457949(90)90016-U. [38] Gutiérrez MA. Energy release control for numerical simulations of failure in quasibrittle solids. Commun Numer Methods Eng 2003;20:19–29. doi:10.1002/cnm.649. [39] Verhoosel C V., Remmers JJC, Gutiérrez MA. A dissipation-based arc-length method for robust simulation of brittle and ductile failure. Int J Numer Methods Eng 2009;77:1290–321. doi:10.1002/nme.2447. [40] Labanda NA, Giusti SM, Luccioni BM. A path-following technique implemented in a Lagrangian formulation to model quasi-brittle fracture. Eng Fract Mech 2018;194:319–36. doi:10.1016/J.ENGFRACMECH.2018.03.004. [41] Oliver J, Huespe AE, Cante JC. An implicit/explicit integration scheme to increase computability of non-linear material and contact/friction problems. Comput Methods Appl Mech Eng 2008;197:1865–89. doi:10.1016/j.cma.2007.11.027.
30
[42] Yu C, Hoogenboom PCJ, Rots JG. Incremental sequentially linear analysis to control failure for quasi-brittle materials and structures including non-proportional loading. Eng Fract Mech 2018;202:332–49. doi:10.1016/j.engfracmech.2018.07.036. [43] Zienkiewicz OC, Taylor RL, Zhu JZ. Finite element method : its basis and fundamentals. Butterworth-Heinemann; 2013. [44] Potts DM, Zdravkovic L. Finite Element Analysis in Geotechnical Engineering: Volume One - Theory. Thomas Telford Ltd; 1999. doi:10.1680/feaiget.27534. [45] Lin Ye. Role of matrix resin in delamination onset and growth in composite laminates. Compos Sci Technol 1988;33:257–77. doi:10.1016/0266-3538(88)90043-7. [46] Turon A, Costa J, Camanho PP, Dávila CG. Simulation of delamination in composites under high-cycle fatigue. Compos Part A Appl Sci Manuf 2007;38:2270–82. doi:10.1016/J.COMPOSITESA.2006.11.009. [47] Mejia C, Paullo L., Roehl D. Continuation Methods for the Simulation of Rock Fracture with Cohesive Elements. Int. Soc. Rock Mech. Rock Eng., 2015. [48] Camanho PP, Davila CG. Mixed-Mode Decohesion Finite Elements for the Simulation of Delamination in Composite Materials. 2002. [49] Turon A, Camanho PP, Costa J, Dávila CG. A damage model for the simulation of delamination in advanced composites under variable-mode loading. Mech Mater 2006;38:1072–89. doi:10.1016/j.mechmat.2005.10.003. [50] Chaisomphob T, Kanok-Nukulchai W, Nishino F. An automatic arc length control algorithm for tracing equilibrium paths of nonlinear structures. Struct Eng Eng 1988;5:227–30. doi:10.2208/jscej.1988.392_227. [51] Mende CA., Gattass M, Roehl D. The GeMA Framework - An Innovative Framework for the Development of Multiphysics and Multiscale Simulations. Proc. VII Eur. Congr. Comput. Methods Appl. Sci. Eng. ECCOMAS Congr., 2016. [52] Manzoli OL, Maedo MA, Bitencourt LAG, Rodrigues EA. On the use of finite elements with a high aspect ratio for modeling cracks in quasi-brittle materials. Eng Fract Mech 2016;153:151–70. doi:10.1016/J.ENGFRACMECH.2015.12.026. [53] Gálvez JC, Elices M, Guinea GV, Planas J. Mixed Mode Fracture of Concrete under Proportional and Nonproportional Loading. Int J Fract 1998;94:267–84. doi:10.1023/A:1007578814070. AUTHORSHIP STATEMENT: Component of the research
Author’s number
substantial contribution to conception and design
1,2,3
substantial contribution to acquisition of data
1,2,3
substantial contribution to analysis and interpretation of data
1,2,3
drafting the article
1,2,3
critically revising the article for important intellectual content
1,2,3
final approval of the version to be published
3
Conflict of Interest All authors have participated in (a) conception and design, or analysis and interpretation of the data; (b) drafting the article or revising it critically for important intellectual content; and (c) approval of the final version. This manuscript has not been submitted to, nor is under review at, another journal or other publishing venue.
31
The authors have no affiliation with any organization with a direct or indirect financial interest in the subject matter discussed in the manuscript
Discrete Fracture Propagation Analysis using a Robust Combined Continuation Method Cristian Mejia , Luis. F. Paullo , Deane Roehl PII: DOI: Reference:
S0020-7683(20)30036-6 https://doi.org/10.1016/j.ijsolstr.2020.02.002 SAS 10600
To appear in:
International Journal of Solids and Structures
Received date: Revised date: Accepted date:
8 July 2019 18 December 2019 1 February 2020
Please cite this article as: Cristian Mejia , Luis. F. Paullo , Deane Roehl , Discrete Fracture Propagation Analysis using a Robust Combined Continuation Method, International Journal of Solids and Structures (2020), doi: https://doi.org/10.1016/j.ijsolstr.2020.02.002
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.
1
Discrete Fracture Propagation Analysis using a Robust Combined Continuation Method Cristian Mejia a, Luis. F. Paullo a,b, Deane Roehl a,b a
Tecgraf Institute, Pontifical Catholic University of Rio de Janeiro, Rua Marquês de São Vicente, 225, Gávea, 22453-900, Rio de Janeiro, RJ, Brasil, www.tecgraf.puc-rio.br b
Department of Civil and Environmental Engineering, Pontifical Catholic University of Rio de Janeiro, Rua Marquês de São Vicente, 225, Gávea, 22453-900, Rio de Janeiro, RJ, Brasil, www.civ.puc-rio.br
Abstract. The prediction of equilibrium paths in the presence of brittle fractures remains a challenge, despite the advances in the fracture modelling techniques. Difficulties arise in the solution of the equilibrium equations in the presence of material and structural instabilities. This work proposes a continuation method with multiple restrictions to solve the nonlinear system of equilibrium equations and overcome global convergence problems owing to the sudden energy release in fracture propagation analyses. The method combines simultaneously the cylindrical arc-length and the dissipation-energy control methods. In addition, this method allows the use of algorithms for the solution of linear systems of equations, which do not require the assembly of the global stiffness matrix. Numerical examples of tensile-mode and mixed-mode fracture propagation are presented and compared with experimental results. A zero-thickness interface element with linear softening damage model simulates mixed-mode fracture. The results show the accuracy, robustness, and effectiveness of the proposed continuation method even for challenging fracture problems. Keywords: Rock fracture, Softening behavior, Cohesive elements, Damage law, Continuation methods.
2
1
Introduction The prediction of material failure is essential in many mechanical and civil engineering
problems. In the last decades, several computational methods have been developed to model fracture processes with efficiency and accuracy, and great advances have been made. Since the first works using the finite element method to model fracture, one of the main challenges has been how to simulate crack propagation with arbitrary paths. One of the first numerical methods for fracture propagation was formulated in the framework of linear elastic fracture mechanics. According to this approach, fracture growth is associated to a mesh update to include the discontinuity[1]. An alternative approach considers each element edge as a potential crack and allows crack opening/propagation when specific criteria are satisfied [2,3]. The same concept is implemented using zero-thickness interface elements with double nodes at each boundary between solid elements [4–7]. This technique has reported very realistic representations of fracture behavior in concrete (meso-mechanical simulation of concrete specimens [8]) and rock (evolution of cracking and fractures around oil wells [9]). However, it leads to larger global systems of equations, more critically in 3D cases. The embedded discontinuity approach is another technique to insert strong discontinuities into existing finite element models [10–12]. For such, suitable condensation techniques are required in order to describe the kinematics properly, to transfer stresses correctly between the cohesive zones, and to eliminate locking [13,14]. This approach requires crack-tracking algorithms to allow crack growth through the finite element mesh. More recently, an extended finite element method has been widely used to model strong discontinuities. Initially proposed by Moës and Belytschko [15], this method does not require the definition of potential cracks and propagation paths a priori. A fracture path-tracing algorithm identifies and inserts the crack path in the model during the analysis. The level set method and the fast marching method have been developed as efficient techniques to track the crack path [16]. The extended finite element method is an attractive technique to model discontinuity. A cohesive zone model simulates the fracture process zone ahead of the crack tip, where material softening takes place [17]. A cohesive constitutive model provides the relation between the cohesive traction on the crack surface and the separation across the crack surfaces. Good overviews and discussions related to cohesive constitutive models are given, for example, in Papoulia et al. [18], Park et al. [17] and Nguyen [7]). The cohesive interface element method, which has its origin in the concept of the cohesive zone model proposed by Barenblatt [19] and Dugdale [20], has been used with great success
3
to simulate fracture processes in concrete, metals, ceramics, composites and rocks. Several research studies such as Hillerborg et al. [21], Xu and Needleman [22] and Camacho and Ortiz [2] have focused on the development of the cohesive zone concept in computational fracture mechanics. In addition, the concept of cohesive zone model has also been applied in conjunction with the extended finite element method XFEM [23]. On the other hand, numerical modeling of materials with softening behavior is prone to convergence problems owing to sudden release of energy in the post-critical regime. In order to overcome this problem, some researchers such as Chaboche et al. [24], Hamitouche et al. [25], Kayhan [26], Jiang [27] and Lapczyk et al. [28] have adopted viscous regularization techniques in discrete fracture propagation. These techniques modify the constitutive equation of the interface element using a viscous damage variable. However, they lead to stresses that drift from the traction-separation law [29,30]. In this context, the numerical analysis of structural problems with strong material softening behavior requires the use of control methods capable of following the complex loading/unloading history of the structure introduced by strong stiffness degradation. Generally, in non-linear geometric and non linear material problems the so-called continuation methods are used. Among them, one of the most used family of continuation methods are arc-length control methods [31–34]. The main pros of the arc-length methods are versatility and easy computational implementation. Nevertheless, many researchers have reported difficulties with the use of the arc-length control methods when analyzing strong material softening problems, as pointed out in Pohl et al. [35] and Paullo and Roehl [36]. The indirect strain increment method controls the monotonic increment of strain at a certain point of interest in the structure and improves convergence of the solution for some softening problems [37]. However, fracture propagation often exhibits non-monotonic strain increase, especially in the region near the crack opening. Methods based on energy dissipation also have good performance in material softening problems [38,39]. These methods require the use of an additional control at non-dissipative steps. Thus, there is the need to define a switch parameter based on a minimum amount of released energy. The definition of an adequate parameter value requires good knowledge of the system behavior. Labanda et al. [40] introduced a restriction based on energy release criterion in an augmented Lagrangian formulation to model cohesive cracks. To overcome convergence difficulties associated with material softening behavior and contact/friction interfaces, Oliver et al. proposed an implicit-explicit integration scheme, which combines implicit integration of stresses in the constitutive model with explicit
4
extrapolation of the internal variables [41]. Recently, incremental sequentially linear analysis (ISLA) has been proposed by Yu et al. [42], which combines the classical Newton-Raphson method with Sequentially Linear Analysis. In this method, each increment starts and ends with an equilibrium state. Equilibrium is updated based on the stiffness reduction of critical elements and load increments. The ISLA method demonstrated robustness to simulate the failure of quasi-brittle material under non-proportional loading. Paullo and Roehl [36] formulated a technique to combine two or more restrictions through the least squares methods, obtaining a mixed control. The mixed control showed good performance in problems with strong geometric nonlinearities and elastoplastic behavior. In addition, the combination reduces the participation of the user in the control of the analysis. This work proposes a step forward of the continuation method with combined constraints formulated by Paullo and Roehl [36] by simplifying the control equation. This simplification allows the use of the combined approach in association with direct solution strategies. In addition, the alternative strategy reduces the computational effort in the solution of the linear system of equations. To deal with strong dissipation in fracture problems the method combines the cylindrical arc-length control and the dissipated energy control. Numerical examples of fracture problems modeled with interface cohesive elements show the robustness and effectiveness of the proposed strategy. Conventional continuation methods cannot handle these examples. 2
Cohesive Interface Element
2.1 Interface element formulation The weak form of the equilibrium equation, including the interface effects, follows from the principle of virtual work. The virtual work done by the body forces ( ) in the domain and by the external traction ( energy in the domain
) on the boundary
are equal to the sum of the virtual strain
and the cohesive fracture energy on the fracture surface ( ) (see
Figure 1). ∫ where
∫ and
∫
∫
Eq. 1
are virtual strains, virtual displacements and virtual fracture opening,
respectively. In addition,
is the stress tensor in the continuum element while
cohesive traction along the fracture surface.
is the
5
𝑇𝑒𝑥
𝛤
𝛺 𝛤𝑓
𝑇𝑐
Figure 1 Schematic representation of a domain ( ) with an internal fracture ( ) under external tractions ( and cohesive tractions ( ).
)
The domain ( ) is discretized into continuum elements and zero-thickness interface elements along the crack path. The displacement field ( ̅) over the domain ( ) is interpolated from the nodal displacement vector ( ) by using shape functions ( ). ̅
Eq. 2 The local separation along the fracture ( ) is approximated from the nodal displacements
by using the global displacement-separation matrix (
). Eq. 3
where
is given by Eq. 4
is the shape function matrix,
is the local displacement-separation matrix and
is the
rotation matrix. Following the standard finite element procedures [43,44], the internal force vector (
) of an interface element is given by
∫
Eq. 5
The gradient of the internal force vector leads to the tangent matrix of the interface element, ∫ where
Eq. 6 is the constitutive matrix relating tractions and relative displacements
on the fracture surface.
6
2.2 Cohesive zone model The cohesive constitutive model can be defined directly in terms of the traction separation law. The traction on the cohesive element in the fracture process zone consists of the normal and the shear
components, which are related to the shear
and normal
separation through the constitutive equations [44]. , Here
[
[ ]{
}
Eq. 7
] is the tangential constitutive matrix of the cohesive element,
the shear elastic stiffness and
is
is the normal elastic stiffness.
This work adopts a cohesive model, which is characterized by a linear softening law. This relationship assumes linear elastic behavior up to the critical point of damage initiation, where the traction reaches the maximum cohesive traction
matching the damage initiation
criterion. The initiation of the softening process follows a quadratic failure criterion, which considers the interaction of the traction components according to [45] ( where
)
is the normal strength,
operator, defined as 〈 〉
| |
.
〈 〉
Eq. 8
/
is the shear strength and 〈 〉 is the Macaulay bracket .
The effective separation of damage initiation cohesive traction
√
〈 〉 . Here
and
corresponds to the point of maximum represent the normal and shear stress
that meet the damage initiation criterion (Eq. 8). Thus, this maximum cohesive strength considers indirectly the effect of normal and shear stress components. Softening behavior occurs in the post-critical zone, describing progressive degradation of the material stiffness. Figure 2 shows the damage model with linear softening [46,47].
7
Point of damage initiation
Traction, 𝜏
𝜏
𝑑𝜏
𝜏𝑜 𝑑 𝜏 𝑘𝑖
Post critical regime
𝐺𝑐
𝑑 𝑘𝑖 Δ𝑓 Δ𝑜 Effective displacement jump Δ
𝑘𝑖
Δ
Figure 2 Damage model with linear softening.
The area under the traction-separation curve represents the cohesive fracture energy . In the post-critical regime, the cohesive traction
̅
̅
can be defined as,
̅ ̅
Eq. 9
{
where d is the damage variable. The damage variable for linear softening can be defined as proposed in [29,48], Eq. 10
where
is the effective separation at damage initiation, √〈
for material failure,
〉
consistent tangent modulus
is the critical effective separation
is the effective displacement. Therefore, the
is obtained by differentiating the cohesive tractions with
respect to the fracture separation [49]: ̇ ̇
{
Eq. 11 .
̅
〈
〉
/.
where ̅ is the Kronecker delta,
̅
〈
〉
/
is the elastic stiffness matrix, and
Eq. 12
is a scalar value
defined as: Eq. 13
8
The consistent tangent modulus for unloading and reloading is a function of the damage variable and is given by, [
Eq. 14
]
Finally, a crack develops when the effective separation reaches the critical value of
and
the cohesive tractions become zero. Other cohesive models can be adopted by changing the damage evolution law given by Eq. 10, such as exponential softening model. 3
Solution of the nonlinear system of equations
3.1 Classical Newton-Raphson solution In a finite element context, the incremental global equilibrium can be defined by the following system of equations Eq. 15 where
is the global nonlinear internal force vector which depends on the global
displacement vector parameter, and
.
The subscript
indicates the -th analysis step,
is the load
is the reference nodal force vector. The increment of the displacements and
of the load factor for step
can be obtained as. Eq. 16
Accordingly, the incremental equilibrium in step
is obtained as Eq. 17
The solution of Eq. 17 requires a restriction in
, which is expressed as Eq. 18
As indicated in Paullo and Roehl [36], an iterative Newton-Raphson scheme is applied to the extended system defining the following 0
where
1{
}
{
is the tangent stiffness matrix,
iteration index, and
Eq. 19
}
and
is the unbalanced force vector,
is the
are the displacement vector and load factor correctors,
respectively. Following the proposal in Paullo and Roehl [36], the following error measures
9
are adopted. |
|
|,
Eq. 20
|
3.2 Cylindrical Arc-Length control method For nonlinear structural problems the cylindrical arc-length equation [31] is frequently used as restriction (Eq. 18) in the Newton Raphson scheme. In this case, the restriction is defined by (
) |
where
Eq. 21
|
is the arc length control parameter. The corrector of the load parameter,
,
follows from the solution of the quadratic equation: Eq. 22 In this equation the coefficients are: |(
)
(
)
(
)
| (
) (
Eq. 23
)
Definition of the increment size is performed by an adaptation procedure introduced by Ramm [34] and also used by Chaisomphob [50]. This scheme defines the parameter
for a
subsequent step as √ where
is the desired number of iterations and
Eq. 24 the number of iterations required to
achieve convergence. Paullo and Roehl [36] point out the need of for higher and lower limits for
in order to achieve adequate control parameter values throughout the incremental
analysis. 3.3
Dissipated energy increment control method
For problems with energy release by damage, the equation of restriction based on the increment of dissipated energy can be obtained following the formulation of Gutierrez [38]. In the geometrically linear case, we have
10
( where
Eq. 25
)
is the desired incremental dissipated energy. We obtain the load parameter
corrector by introducing Eq. 25 in the extended equilibrium system in Eq. 19, (
)
(
)
Note that
and
Eq. 26
correspond to the converged load factor and displacement vector
from the previous step, respectively. It is important to point out that for elastic incremental steps, where no dissipation occurs, this control fails and alternative ones must be adopted. 3.4
Full combined restriction
Aiming at providing a more robust solution strategy in problems with both material softening and geometric nonlinearity, Paullo Muñoz e Roehl [36] formulate a combination of two or more restrictions working simultaneously. According to [36], the general non-linear system of equations defined in Eq. 17 to Eq. 19 can include
restrictions grouped in a
-
dimensional vector defined by: ̃
where ̃
̃
{
[ ̃ ̃
̃
Eq. 27
}
̃ ], is the vector containing the control parameters. Introducing the
vector of restrictions in the Newton Raphson scheme results in the following inconsistent system of equations. 0
1{
}
{
Eq. 28
}
Paullo and Roehl [36] applied the least squares technique to obtain an approximate solution with the participation of all restrictions. Thus, the corrector of displacements and the load factor follow directly from {
}
(0
1 0
1)
0
1 {
}
Eq. 29
Control parameter values at each iteration are recalculated using the corresponding single
11
constraint equation to guarantee convergence of the solution. The converged values of the control parameters are used to estimate the predictor of the following step. In addition, each control parameter is adapted according to Eq. 24. Convergence rate and robustness of the combined constraint technique are studied in [36]. From Eq. 29 the load factor corrector ( ((
)
((
)
follows as
)
)[
)[
]( (
)
+
] ((
)
)
(
)
(
)
)
Eq. 30
)
where *(
)
(
Eq. 31
Consequently, the displacement vector corrector is * (
)
(
)
((
)
(
)
)
+
Eq. 32
3.5 Simplified combined restriction A combined restriction technique requires a full coefficient matrix
for the solution of
the linear system of equations. Consequently, the possibility of employing algorithms that take advantage of sparsity of the stiffness matrix is lost, increasing strongly the computational effort. Moreover, some solution schemes do not assemble the global stiffness matrix. In this context, a modified constraint combination technique is proposed to allow working with general solution methods. We can define: (
)
( were
Eq. 33
) and
are the displacement solutions corresponding to the unbalanced forces and
the reference forces, respectively. The displacement corrector derives from Eq. 28 to Eq. 34 Substitution Eq. 34 in the second line of Eq. 28, establishes the following inconsistent system of equations {
}
{
Eq. 35
}
Applying the least squares method to solve Eq. 35 gives the load parameter corrector {
} {
}
{
} {
}
Eq. 36
12
It is important to point out that when none of the constraint equations depend on the load parameter, i.e.,
, Eq. 30 and Eq. 35 provide the same result. If we consider the
combination of the cylindrical arc-length control and the dissipated energy control, we add the following restrictions to the iterative procedure 3,
2
2
Eq. 37
3
Substituting Eq. (34) in Eq. (33) we obtain (
)
( (
(
) (
) )
(
) (
(
))
)
( (
)
( )
(
)
(
(
)
Eq. 38
))
Alternatively, we can write the following expression for the predictor of the load parameter, which considers the positive work in the arc-length predictor contribution: ( ((
)
) )
( (
(
)
Eq. 39
))
The proposed combination technique does not require the assembly of the full coefficient matrix, reduces the computational effort, and allows using direct solvers. In addition, the proposed method combines the arc-length and the dissipated energy controls. The values of and
in the first step and first iteration are calculated from Eq. 22 and Eq. 25,
respectively, using an initial prescribed load parameter
. It is important to point out that Eq.
38 works for steps with and without dissipation. Thus, the use of a switch procedure is not necessary. For non-dissipative steps, Eq. 38 reduces to a linearized cylindrical arc-length equation. 4
Numerical Examples The zero-thickness interface element and the solution strategies proposed are implemented
in the in-house multiphysics framework called GeMA (Geo Modelling Analysis) [51]. Three examples show the robustness and effectiveness of the proposed strategy. 4.1 L-Shaped Plate (Single Fracture Trajectory) This section assesses fracture propagation in an L-Shaped Plate. Figure 3 presents the structure dimensions and loading. This example illustrates the case of a single fracture path predefined by a line of cohesive elements. The configuration and properties of the problem are taken from Pohl et al. [35]. Plane strain four-node quadrilateral elements model the
13
structure. The solid elements have linear elastic behavior with Young’s modulus E = 100kN/mm2 and Poisson’s coefficient
. The line of cohesive elements is composed by
four-node interface elements. The cohesive elements have the following constitutive kN/mm3, shear elastic stiffness
parameters: normal elastic stiffness kN/mm3, normal strength
N/mm2, shear strength
effective total displacement at failure al. [35]. The elastic stiffness values
N/mm2 and the
mm, in accordance with Pohl et and
are computed by the following relation
. In order to reduce the effect of this artificial stiffness on the elastic response, we adopt thickness
= 1 mm [52].
5 𝑚
cohesive elements
5 𝑚
λ×
5 𝑚
𝑘𝑁 𝑚
x
5 𝑚
Figure 3 Geometry and FEM mesh of L–shaped plate.
Figure 4 shows the equilibrium path for the horizontal displacement
obtained with the
cylindrical arc-length control and with the combined control methods. The cylindrical arclength control fails in the region near the maximum load factor showing its difficulty to analyze fracture propagation problems characterized by material softening behavior. The works Paullo and Roehl [36] and Pohl et al. [35] present similar problems related to the use of the arc-length control methods in structures with material non linearity. In those cases, the arc-length presented convergence problems in analyses involving plasticity with linear softening and damage. On the other hand, the proposed combined method showed good performance and was able to carry out the analysis up to the stage of total opening of all
14
cohesive elements (total crack propagation). Figure 4 also shows a very complex equilibrium path exhibiting local loading and unloading regions in which each “finger” corresponds to fracture propagation of a single cohesive element. These fingers (local snap-backs) lead to failure of traditional solution control algorithms. However, in this example, finer meshes render smoother solutions, allowing the use of standard arc-length method. Figure 5 shows the deformed shape of the structure when 9 elements reach the limit damage factor d =1. It is noticeable that the solid elements do not exhibit important geometric distortions. In addition, the main mechanism of fracture is normal traction. Cylindrical arc-length
Figure 4 Load factor
vs. horizontal displacement Ux1 considering
Proposed method
and
5.
Figure 5 Final deformed configuration where nine cohesive elements collapsed, scale factor = 5.
4.2 Plate with eccentric hole This section presents fracture propagation through a square plate with an eccentric hole
15
under uniaxial stress.
Figure 6 presents the plate dimensions, boundary and loading
conditions. The aim of this problem is to assess the performance of the proposed algorithm in snap-back situations. The geometry and properties of the problem are from Labanda et al. [40]. Four-node quadrilateral elements discretize the structure under plane stress condition. The solid elements have linear elastic behavior with Young’s modulus E = 20 GPa and Poisson’s coefficient
. The fracture behavior is modeled using four-node zero
thickness interface elements placed horizontally, as shows Figure 6. The cohesive elements have the following constitutive parameters: normal elastic stiffness elastic stiffness
GPa/mm, normal strength
GPa/mm, shear
MPa, shear strength
MPa . The effective total displacement at failure is corresponds to a fracture energy release of
mm. This value N/mm, as adopted in Labanda et al.[40]. λ×
𝑘𝑁 𝑚
5𝑚 P Cohesive elements
𝑚
5𝑚
5𝑚
𝑚
5𝑚
Figure 6 Geometry and FEM mesh of plate with hole.
Figure 7 shows the equilibrium path in terms of the load factor versus the vertical displacement of point P. This example was simulated using 200 load steps. The proposed combined method is able to track complex equilibrium paths including snap back phenomena. Figure 7 also shows that the equilibrium path obtained matches the numerical results reported by Labanda et al. [40] that were obtained with a switched energy dissipated control. Figure 8 shows four stages, A, B, C, D, according to the damage progress. Up to state A, the peak load, the fracture propagates in both directions of the hole (see Figure 8a). Then, the fracture develops on the right side of the plate producing unloading behavior until point B,
16
when the fracture reached the plate border. Between B and C fracture aperture increases. After point C, the fracture propagates through the left side of the plate until point D, when the left fracture reached the plate border (see Figure 8d).
A
C B
D
Figure 7 Load factor vs Vertical displacement in P. A
B Initiation
D C
Figure 8 Damage evolution for the plate with hole under uniaxial stress: states A, B, C and D.
Figure 9 shows the damage evolution along the interface elements placed horizontally for the states presented in Figure 8 (A, B, C and D).
17
Left side
Hole
Right side
Left side
(a)
Hole
Right side
(b)
Figure 9 Damage evolution along the fracture path for a) state A and B, and b) state C and D.
4.3 Perforated cantilever beam: Single crack trajectory and multiple snap-back path This example presents the delamination process in a perforated cantilever beam that has been extensively studied in the literature on continuation methods for fracture and softening problems. Initially studied by Verhoosel et al. [39], this structure attracts the attention of researchers, such as Pohl et al. [35], Labanda et al. [40] and Paullo and Roehl [36], once the load-displacement path exhibits several snap-back phenomena along the delamination process. Verhoosel studied this problem using an energy release control switched with load control for non-dissipative steps. Figure 10 shows the geometry, mesh, boundaries and loading conditions. The geometry and properties of the structure are from Paullo and Roehl [36] and Verhoosel et al. [39], respectively. A line of linear cohesive elements defines the single crack trajectory, with the following properties: GPa/mm,
MPa,
GPa/mm,
MPa The effective total displacement at failure is
5 mm that corresponds to a fracture energy release of
5
N/mm. Plane strain four-node quadrilateral elements discretize the linear elastic solid domain with Young’s modulus thickness of
mm.
MPa and Poisson’s coefficient
. The model has a
18
Figure 10 Geometry and finite element model.
Figure 11 shows the vertical displacement at the top of the free end of the beam vs the load factor and the deformed shape of the structure after 150 converged steps. That stage corresponds to the total collapse of five segments. As expected, the equilibrium path exhibits several load and displacement limit points, as well as several snap-back forming a “finger” as each segment collapses. This behavior agrees with the response in Verhoosel et al. [39] and shows the robustness of the proposed continuation method to trace the equilibrium path with several snap-back situations. (a)
Figure 11 Load factor
(b)
vs. vertical displacement and b) deformed shape after 150 steps with the proposed method. and 5.
4.4 Four-Point Bending Test (Multiple fracture trajectory) In this section, crack propagation in a four-point bending beam test is assessed. The geometry and properties of the system are from Gálvez et al. [53], which shows experimental results mainly for two kinds of support condition. Figure 12 shows the load and support conditions for the two models presented in [53]. The beam models have a thickness of 50 mm. The main difference between the two configurations is the support condition. This
19
difference targets at assessing different fracture mechanisms since in configuration Type-I it is expected that fracture occurs due to normal tractions, while in configuration Type-II, fracture is due to mixed normal and shear tractions, as indicated in Gálvez et al. [53]. 𝐷
5𝑚𝑚 𝑃
𝐷
𝑇𝑌𝑃𝐸
𝐼
𝐷
𝐶𝑜 𝑒𝑠𝑖𝑣𝑒 𝑟𝑒𝑔𝑖𝑜𝑛
𝑃
𝐷
𝑇𝑌𝑃𝐸
5𝑚𝑚
𝐼𝐼
𝐶𝑜 𝑒𝑠𝑖𝑣𝑒 𝑟𝑒𝑔𝑖𝑜𝑛
5𝑚𝑚 𝐷
𝐷 𝐷
𝐷
𝐷
𝐷
𝐷
𝐷
𝐷
𝐷
𝐷
5𝐷
𝐷
𝐷
𝐷
Figure 12 Geometry and properties of beams Type - I and Type – II, taken from Galvez [46].
The material behavior of the solid elements is linear elastic while the cohesive elements present linear softening behavior. Since in the experiments fractures develop in the central region of the beam, the finite element discretization includes cohesive interface elements surrounding 3-node triangular solid elements in this region. The use of triangular elements in this region aims at allowing a more alternatives for fracture paths. Outside this region, discretization follows with 4-node quadrilateral solid elements. Figure 13 shows the finite element mesh of the beam model where the bounding box identifies the region with a dense mesh of triangles with cohesive elements at its sides. In this way, the crack trajectory is not predefined. Actually, several trajectories appear during the analysis.
𝑃
𝐶𝑜 𝑒𝑠𝑖𝑣𝑒 𝑒𝑙𝑒𝑚𝑒𝑛𝑡 𝑚𝑒𝑠
Figure 13 FEM mesh for beams Type - I and Type - II
In accordance with the properties given in Galvez [46], the material has a Young’s modulus of E=38 GPa and a Poisson’s ratio of following parameters:
GPa/m,
,. The cohesive model takes the GPa/m,
MPa,
MPa
20
and
mm. Figure 14 shows the deformed shape and the accumulated scalar
damage parameter in the cohesive elements for model Type-I after the crack propagates through the cohesive region (200 converged steps). This stage corresponds to a crack-mouth opening displacement (CMOD) equal to 0.42mm. The main crack trajectory is apparent. In addition, the distribution of accumulated damage indicates that several secondary fracture trajectories exist around the main one. Figure 15 shows that the numerical main crack trajectory (solid circled line) obtained numerically matches the experimental envelope (shaded area) presented in Galvez et al. [53], indicating good performance of the cohesive crack approach to represent real fracture processes. Figure 15 also shows the numerical damage envelope for accumulated damage greater than 0.85 (dashed line) which contains the secondary trajectories. The secondary trajectories also coincide with the experimental envelope.
𝐴𝑐𝑐𝑢𝑚𝑢𝑙𝑎𝑡𝑒𝑑 𝑑𝑎𝑚𝑎𝑔𝑒
Figure 14 Deformed shape and accumulated damage factor for beam Type – I for CMOD = 0.42mm. Scale factor = 20.
21
Figure 15 Main numerical crack trajectory and damaged region (D>=0,85) for beam Type – I for CMOD = 0.42mm.
Figure 16 shows the relation between the crack-mouth opening displacement (CMOD) and load P. In Figure 16, the maximum load level agrees with the experimental envelope. In addition, the curve P vs CMOD follows coherently the experimental results. It is important to point out that these numerical results were obtained using the proposed combined method, since the cylindrical arc-length control failed in the beginning of crack propagation.
CMOD
Figure 16 Vertical Load P vs. crack mouth opening displacement (CMOD) for beam Type – I.
For model Type-II Figure 17 presents the deformed shape and the accumulated damage at
22
the Gauss integration points of the cohesive elements after final crack propagation (380 converged steps and CMOD = 0.13mm). From beam Type-I to Type-II the change in the support conditions lead to a modification of the main crack trajectory. In beam Type II damage in the cohesive elements is generated by the action of normal and tangential tractions, while in beam Type I fracture occurs mainly due to normal tractions. In addition, beam Type II presents a second important damaged zone located at the upper-left part of the cohesive region. In that case, during the fracture process, the normal traction over this zone increases. However, damage does not reach the maximum value indicating that the interface elements are in the softening process but do not reach total degradation.
𝐴𝑐𝑐𝑢𝑚𝑢𝑙𝑎𝑡𝑒𝑑 𝑑𝑎𝑚𝑎𝑔𝑒
Figure 17 Deformed shape and accumulated damage for beam Type – II for CMOD = 0.13. Scale Factor = 20.
23
Figure 18 Main crack trajectory and damaged region for beam Type – II for CMOD = 0.13.
Figure 18 shows the agreement between the numerical and the experimental results. The numerical main trajectory and the numerical damaged region around it coincide with the experimental envelope for beam type - II. Figure 19 also reflects the good agreement between the experimental and numerical results in terms of the P vs. CMOD curves. Another important aspect refers to the computational effort to reach the maximum crack propagation. Despite the complexity in the analysis, the proposed simplified combined method showed good performance, controlling the analysis until total crack propagation without interruptions in the nonlinear solution procedure. In this way, the proposed continuation method shows robustness for this kind of problem.
24
CMOD
Figure 19 Load vs. CMOD for beam Type – II.
4.5 Highly pre-stressed beam test This problem shows the robustness of the proposed combined method to simulate a prestressed beam under non-proportional loading. Figure 20 shows the beam geometry, boundary conditions and loading conditions. The geometry and properties of the beam are from Yu et al. [42]. Four-node quadrilateral elements discretize the structure under plane stress condition. Fracture behavior is modeled using four-node zero thickness interface elements, which are placed at the edge of all solid elements. The properties are Young’s modulus E = 32 GPa, Poisson’s coefficient fracture energy release of
5𝑚
, tensile strength
MPa, and
N/mm, as adopted in Yu et al. [42]. λ
5𝑚
λ
5𝑚
A P P
5𝑚
Cohesive elements at all edges
5𝑚
Figure 20 Geometry and FEM mesh of pre-stressed beam.
5𝑚
𝑚
25
The beam is subjected to the combination of horizontal load P that is contant during the simulation and vertical load
, as illustrated in Figure 20. Two cases are studied considering
P of 1 or 5 MPa. The proposed algorithm controls the vertical load, which increases up to failure and reduces afterwards. Figure 21 indicates the final crack path at the end of the simulation for both cases (prestressed beam with P =1 and P = 5 MPa). For both cases, multiple cracks are generated at the base of the beam, and after that, two main cracks propagate to the top.
a)
Multiple cracks b)
Figure 21 Prestressed beam: damage evolution and final crack pattern for a) P=1 MPa and b) P=5 MPa prestress.
The damages regions obtained are similar to those reported by Yu et al.[42]. Figure 22 shows the reaction force and vertical displacement in A ( midspan displacement) for both cases. These results match with those obtained by Yu et al.[42]. However, for 1 MPa, the curves are different because the propagation of multiple cracks produces different damage regions.
a)
b)
Figure 22 Prestressed beam: reaction force vs midspan displacement for a) 1 MPa and b) 5 MPa prestress.
26
Figure 23 shows the total number of iterations required in each load step for both cases. For the beam under a prestress of 1 MPa, the total number of steps is 100 with an average of 4 iterations per step. For the beam under prestress of 5 MPa, the curve is simulated in 200 load steps with an average of 6 iterations per step.
a)
b) 3rd attempt
3rd
2nd attempt
2nd
1st attempt
1rd
Figure 23 Number of iterations required in each load step for a) 1 MPa and b) 5 MPa prestress.
Figure 24 presents the evolution along the analysis of the arc-length increments and of the energy increments that satisfy the combined method. The contribution of both parameters ( ,
) is present throughout the analysis.
Figure 24 Evolution of the control parameter
5
and
in the combined method,
5 and
.
Conclusions and Remarks This paper presents a numerical model for fracture propagation using zero-thickness
interface elements with a linear softening cohesive model. This approach in association with a combined restriction control method is proposed to overcome convergence problems due to strong energy dissipation by fracture propagation.
The chosen restrictions were the
cylindrical arc-length and the dissipation energy. The numerical models were carried out
27
considering intrinsic cohesive zone models with very stiff initial elastic slope in order to reduce the effect of artificial compliance. An L-shaped plate with a predefined fracture path validates the solution techniques. The proposed continuation method with the combination of restrictions overcomes the convergence problems reported in [30] owing to the very complex equilibrium path generated by the local loading/unloading at crack tip. The square plate with an eccentric hole and the delamination of a perforated cantilever beam show the capability of the proposed combined method to deal with complex equilibrium paths, including multiple snap-back situations. In a third problem, the experimental results published by Galvez et al. (1998) validate both tensile mode and mixed-mode fracture in concrete under non-proportional loading. This example is a challenging benchmark for numerical fracture models considering multiple fracture paths. For this example, a target region is modelled with a zero-thickness interface element at the border of the continuum element and several fractures appear during the analysis. Furthermore, good matching of results for load-CMOD curves, ultimate-load capacity and fracture path within the experimental envelope limits highlight the strengths of the proposed strategy. The highly prestressed beam demonstrates the robustness of the proposed combined control for non-proportional loading with multiple competing cracks. According to the numerical results, the combined control overcomes the global convergence problem associated with the damage constitutive model, which also includes instabilities related to the derivative of the discontinuous Macaulay function. Finally, the combined control shows good performance to analyze the fracture propagation process in all the examples assessed. The combined restriction method is less dependent on parameter settings. The automatic step size definition also contributes in this sense. In addition, the proposed technique allows the use of optimized solution methods for systems of linear equations. REFERENCES [1]
[2] [3]
Ingraffea AR, Saouma V. Numerical modeling of discrete crack propagation in reinforced and plain concrete. Fract. Mech. Concr. Struct. Appl. Numer. Calc., Dordrecht: Springer Netherlands; 1985, p. 171–225. doi:10.1007/978-94-009-61524_4. Camacho GT, Ortiz M. Computational modelling of impact damage in brittle materials. Int J Solids Struct 1996;33:2899–938. doi:10.1016/0020-7683(95)00255-3. Ciancio D, Carol I, Cuomo M. On inter-element forces in the FEM-displacement formulation, and implications for stress recovery. Int J Numer Methods Eng 2006;66:502–28. doi:10.1002/nme.1564.
28
[4] [5] [6]
[7]
[8]
[9] [10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18] [19] [20] [21]
Rots JG, Schellekens JC. Interface elements in concrete mechanics. Comput Anal Des Concr Struct 1990;2:909–18. López CM, Carol I, Aguado A. Microstructural analysis of concrete fracture using interface elements. 2000. Carol I, López CM, Roa O. Micromechanical analysis of quasi brittle materials using fracture based interface elements. Int J Numer Methods Eng 2001;52:193–215. doi:10.1002/nme.277. Nguyen VP. Discontinuous Galerkin/extrinsic cohesive zone modeling: Implementation caveats and applications in computational fracture mechanics. Eng Fract Mech 2014;128:37–68. doi:10.1016/j.engfracmech.2014.07.003. Ciancio D, Lopez CM, Carol I, Cuomo M. New results in meso-mechanical modelling of concrete using fracture-based zero-thickness interface elements. Comput. Model. Concr. Struct., 2003, p. 171–7. Garolera D, Carol I, Papanastasiou P. Micromechanical analysis of sand production. Int J Numer Anal Methods Geomech 2019;43:1207–29. doi:10.1002/nag.2892. Dvorkin EN, Cuitiño AM, Gioia G. Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions. Int J Numer Methods Eng 1990;30:541–64. doi:10.1002/nme.1620300311. Manzoli OL, Shing PB. A general technique to embed non-uniform discontinuities into standard solid finite elements. Comput Struct 2006;84:742–57. doi:10.1016/J.COMPSTRUC.2005.10.009. Radulovic R, Bruhns OT, Mosler J. Effective 3D failure simulations by combining the advantages of embedded Strong Discontinuity Approaches and classical interface elements. Eng Fract Mech 2011;78:2470–85. doi:10.1016/J.ENGFRACMECH.2011.06.007. Jirásek M. Comparative study on finite elements with embedded discontinuities. Comput Methods Appl Mech Eng 2000;188:307–30. doi:10.1016/S00457825(99)00154-1. Zhang Y, Lackner R, Zeiml M, Mang HA. Strong discontinuity embedded approach with standard SOS formulation: Element formulation, energy-based crack-tracking strategy, and validations. Comput Methods Appl Mech Eng 2015;287:335–66. doi:10.1016/J.CMA.2015.02.001. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Methods Eng 1999;46:131–50. doi:10.1002/(SICI)10970207(19990910)46:1<131::AID-NME726>3.0.CO;2-J. Fries T-P, Baydoun M. Crack propagation with the extended finite element method and a hybrid explicit-implicit crack description. Int J Numer Methods Eng 2012;89:1527–58. doi:10.1002/nme.3299. Park K, Paulino GH, Roesler JR. A unified potential-based cohesive model of mixedmode fracture. J Mech Phys Solids 2009;57:891–908. doi:10.1016/J.JMPS.2008.10.003. Papoulia KD, Sam C-H, Vavasis SA. Time continuity in cohesive finite element modeling. Int J Numer Methods Eng 2003;58:679–701. doi:10.1002/nme.778. Barenblatt GI. The Mathematical Theory of Equilibrium Cracks in Brittle Fracture. Adv Appl Mech 1962;7:55–129. doi:10.1016/S0065-2156(08)70121-2. Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8:100– 4. doi:10.1016/0022-5096(60)90013-2. 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. Cem Concr Res 1976;6:773–81. doi:10.1016/0008-8846(76)90007-7.
29
[22] Xu X-P, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 1994;42:1397–434. doi:10.1016/0022-5096(94)90003-5. [23] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Eng Fract Mech 2002;69:813–33. doi:10.1016/S0013-7944(01)00128-X. [24] Chaboche JL, Feyel F, Monerie Y. Interface debonding models: a viscous regularization with a limited rate dependency. Int J Solids Struct 2001;38:3127–60. doi:10.1016/S0020-7683(00)00053-6. [25] Hamitouche L, Tarfaoui M, Vautrin A. An interface debonding law subject to viscous regularization for avoiding instability: Application to the delamination problems. Eng Fract Mech 2008;75:3084–100. doi:10.1016/J.ENGFRACMECH.2007.12.014. [26] Kayhan E. Application of viscous and non-local integral type regularization schemes for softening plasticity. Atilim University, 2009. [27] Jiang H. Cohesive zone model for carbon nanotube adhesive simulation and fracture/fatigue crack growth. University of Akron, 2010. [28] Lapczyk I, Hurtado JA. Progressive damage modeling in fiber-reinforced materials. Compos Part A Appl Sci Manuf 2007;38:2333–41. doi:10.1016/J.COMPOSITESA.2007.01.017. [29] Dávila CG, Camanho PP, Turon A. Cohesive elements for shells. NASA/Technical Publ 2007. [30] Simonovski I. Cohesive zone modeling of intergranular cracking in polycrystalline aggregates. Nucl Eng Des 2015;283:139–47. doi:10.1016/J.NUCENGDES.2014.09.041. [31] Riks E. An incremental approach to the solution of snapping and buckling problems. Int J Solids Struct 1979;15:529–51. doi:10.1016/0020-7683(79)90081-7. [32] Batoz J-L, Dhatt G. Incremental displacement algorithms for nonlinear problems. Int J Numer Methods Eng 1979;14:1262–7. doi:10.1002/nme.1620140811. [33] Crisfield MA. A fast incremental/iterative solution procedure that handles “snapthrough.” Comput Struct 1981;13:55–62. doi:10.1016/0045-7949(81)90108-5. [34] Ramm E. Strategies for Tracing the Nonlinear Response Near Limit Points. Nonlinear Finite Elem. Anal. Struct. Mech., Berlin, Heidelberg: Springer Berlin Heidelberg; 1981, p. 63–89. doi:10.1007/978-3-642-81589-8_5. [35] Pohl T, Ramm E, Bischoff M. Adaptive path following schemes for problems with softening. Finite Elem Anal Des 2014;86:12–22. doi:10.1016/J.FINEL.2014.02.005. [36] Paullo Muñoz LF, Roehl D. A Continuation method with combined restrictions for nonlinear structure analysis. Finite Elem Anal Des 2017;130:53–64. doi:10.1016/J.FINEL.2017.02.003. [37] Chen Z, Schreyer HL. A numerical solution scheme for softening problems involving total strain control. Comput Struct 1990;37:1043–50. doi:10.1016/00457949(90)90016-U. [38] Gutiérrez MA. Energy release control for numerical simulations of failure in quasibrittle solids. Commun Numer Methods Eng 2003;20:19–29. doi:10.1002/cnm.649. [39] Verhoosel C V., Remmers JJC, Gutiérrez MA. A dissipation-based arc-length method for robust simulation of brittle and ductile failure. Int J Numer Methods Eng 2009;77:1290–321. doi:10.1002/nme.2447. [40] Labanda NA, Giusti SM, Luccioni BM. A path-following technique implemented in a Lagrangian formulation to model quasi-brittle fracture. Eng Fract Mech 2018;194:319–36. doi:10.1016/J.ENGFRACMECH.2018.03.004. [41] Oliver J, Huespe AE, Cante JC. An implicit/explicit integration scheme to increase computability of non-linear material and contact/friction problems. Comput Methods Appl Mech Eng 2008;197:1865–89. doi:10.1016/j.cma.2007.11.027.
30
[42] Yu C, Hoogenboom PCJ, Rots JG. Incremental sequentially linear analysis to control failure for quasi-brittle materials and structures including non-proportional loading. Eng Fract Mech 2018;202:332–49. doi:10.1016/j.engfracmech.2018.07.036. [43] Zienkiewicz OC, Taylor RL, Zhu JZ. Finite element method : its basis and fundamentals. Butterworth-Heinemann; 2013. [44] Potts DM, Zdravkovic L. Finite Element Analysis in Geotechnical Engineering: Volume One - Theory. Thomas Telford Ltd; 1999. doi:10.1680/feaiget.27534. [45] Lin Ye. Role of matrix resin in delamination onset and growth in composite laminates. Compos Sci Technol 1988;33:257–77. doi:10.1016/0266-3538(88)90043-7. [46] Turon A, Costa J, Camanho PP, Dávila CG. Simulation of delamination in composites under high-cycle fatigue. Compos Part A Appl Sci Manuf 2007;38:2270–82. doi:10.1016/J.COMPOSITESA.2006.11.009. [47] Mejia C, Paullo L., Roehl D. Continuation Methods for the Simulation of Rock Fracture with Cohesive Elements. Int. Soc. Rock Mech. Rock Eng., 2015. [48] Camanho PP, Davila CG. Mixed-Mode Decohesion Finite Elements for the Simulation of Delamination in Composite Materials. 2002. [49] Turon A, Camanho PP, Costa J, Dávila CG. A damage model for the simulation of delamination in advanced composites under variable-mode loading. Mech Mater 2006;38:1072–89. doi:10.1016/j.mechmat.2005.10.003. [50] Chaisomphob T, Kanok-Nukulchai W, Nishino F. An automatic arc length control algorithm for tracing equilibrium paths of nonlinear structures. Struct Eng Eng 1988;5:227–30. doi:10.2208/jscej.1988.392_227. [51] Mende CA., Gattass M, Roehl D. The GeMA Framework - An Innovative Framework for the Development of Multiphysics and Multiscale Simulations. Proc. VII Eur. Congr. Comput. Methods Appl. Sci. Eng. ECCOMAS Congr., 2016. [52] Manzoli OL, Maedo MA, Bitencourt LAG, Rodrigues EA. On the use of finite elements with a high aspect ratio for modeling cracks in quasi-brittle materials. Eng Fract Mech 2016;153:151–70. doi:10.1016/J.ENGFRACMECH.2015.12.026. [53] Gálvez JC, Elices M, Guinea GV, Planas J. Mixed Mode Fracture of Concrete under Proportional and Nonproportional Loading. Int J Fract 1998;94:267–84. doi:10.1023/A:1007578814070. AUTHORSHIP STATEMENT: Component of the research
Author’s number
substantial contribution to conception and design
1,2,3
substantial contribution to acquisition of data
1,2,3
substantial contribution to analysis and interpretation of data
1,2,3
drafting the article
1,2,3
critically revising the article for important intellectual content
1,2,3
final approval of the version to be published
3
Conflict of Interest All authors have participated in (a) conception and design, or analysis and interpretation of the data; (b) drafting the article or revising it critically for important intellectual content; and (c) approval of the final version. This manuscript has not been submitted to, nor is under review at, another journal or other publishing venue.
31
The authors have no affiliation with any organization with a direct or indirect financial interest in the subject matter discussed in the manuscript