ACCURATE SIMULATION OF MIXED-MODE COHESIVE CRACK PROPAGATION IN PLAIN AND FIBRE-REINFORCED CONCRETE USING XFEM B.L. Karihaloo1∗ and Q.Z. Xiao 2 1)
School of Engineering, Cardiff University, UK
2)
LUSAS FEA Ltd, Kingston-upon-Thames, UK
ABSTRACT This paper discusses the asymptotic fields ahead of mixed mode frictional cohesive cracks in plain and fibre-reinforced concrete. These fields have been derived after reformatting the cohesive-law of plain concrete in a special but universal polynomial containing fractional or integer powers. For a high performance fibre-reinforced concrete, the cohesive law is already known as a polynomial. This polynomial form ensures that the radial and angular variations of the asymptotic fields are separable as in the Williams expansions for a traction-free crack. The coefficients of the expansions however depend nonlinearly on the softening law and the boundary conditions. The leading terms of the true displacement asymptotic field are derived explicitly. These are especially useful as the enrichment functions at the tip of a mixed mode cohesive crack for accurate simulation of crack growth in plain and fibre-reinforced concrete using the extended finite element (XFEM). Several examples of these simulations will be presented to highlight the advantages of XFEM. KEYWORDS Cohesive crack, plain concrete, fibre-reinforced concrete, asymptotic fields, XFEM
INTRODUCTION The cohesive zone (or crack) model of Hillerborg et al. [1] has been extensively used in the study of localization and failure in concrete structures. de Borst et al. [2] have given a concise overview of the various ways for the numerical implementation of the cohesive zone methodology. The knowledge of the asymptotic crack tip displacement fields is especially useful in the recently developed extended/generalized finite element methodology (XFEM) (see, e.g. [3 - 6]). The XFEM enriches the standard local FE approximations with known information about the problem, with the use of the partition of unity (PU). It avoids meshes conforming with the discontinuity and adaptive remeshing as the discontinuity grows as is the case with the FEM. Karihaloo and Xiao [6] have demonstrated that for a crack with traction-free faces, when the crack tip asymptotic field is available and used as enrichment function, the XFEM not only avoids using a mesh conforming with the crack but is also more accurate than FEM. Hence XFEM can use a much coarser mesh around the crack tip. However, when the enrichment function differs from the true asymptotic crack tip field, the mesh needs to be refined in the same manner as in the FEM. Thus it is necessary to know the true asymptotic displacement field around a cohesive crack tip in order to exploit fully the advantages of XFEM. The application of the XFEM to the simulation of the growth of cohesive cracks in quasi-brittle materials has received considerable attention. However, adjacent to the cohesive crack tip, the enrichment function is chosen as the jump function [7 - 9] or a (branch) function which does not represent the true asymptotic nature of the displacement/stress field there [10]. Since there is no ∗
Corresponding author
[email protected]
singularity at the tip of a cohesive crack, a stress criterion is often used to judge the initiation and propagation of the crack. Therefore a reliable analysis of cohesive crack propagation requires an accurate knowledge of the crack tip field. However, although no singularity exists at the tip of a cohesive crack, the stresses obtained by direct differentiation of the displacements are not accurate, and cannot be used to predict accurately the growth of the tip, exactly as in the traction-free cracks. Recently, the authors [11, 12] obtained universal asymptotic expansions at a cohesive crack tip, analogous to the Williams’ expansions at a traction-free crack tip for any normal cohesion-separation law (i.e. softening law) that can be expressed in a special polynomial with integer or fractional powers. This special form ensures that the radial and angular variations of the asymptotic fields are separable as in the Williams expansions. The coefficients of the expansions of course depend nonlinearly on the softening law and the boundary conditions. They demonstrated that many commonly-used cohesionseparation laws, e.g. rectangular, linear, bilinear and exponential, can indeed be expressed very accurately in this special form. They also obtained universal asymptotic expansions when the cohesive crack faces are subjected to Coulomb friction. Xiao et al. [13] used the true crack tip asymptotic displacement field obtained by Xiao & Karihaloo [11, 12] as crack tip enrichment function in the XFEM for the simulation of mode I cohesive crack propagation in quasi-brittle materials to improve the prediction of the stress field ahead of the cohesive crack even with a coarse mesh. In the cohesive cracks, the friction is considered for a finite opening. In this sense frictional cohesive cracks are different from the frictional contact of crack faces, where the crack faces are in contact and not open. However, in cohesive cracks, although the crack faces are not in contact because of the applied cohesive stresses, frictional forces can come into play between the faces when there is relative sliding. Many studies on the mixed mode cohesive cracks can also be found in the literature, but there is doubt about the accuracy of the cohesion-sliding relation because it is difficult to isolate it from frictional forces between the rough cohesive crack faces in quasi-brittle materials such as concrete. This paper will extend the methodology proposed in [13] to the simulation of general mixed-mode cohesive frictional crack propagation in plain and fibre-reinforced concrete. The true crack tip mixedmode asymptotic displacement field will be used as the enrichment function at the cohesive crack tip both with and without Coulomb friction. Typical specimens made of plain and fibre-reinforced concrete will be analyzed to demonstrate the usefulness and potential of this methodology.
GENERAL POLYNOMIAL REINFORCED CONCRETE
COHESION-SEPARATION
LAWS
FOR
PLAIN
AND
FIBRE-
Cornelissen et al. [14] introduced the following exponential cohesion-separation relation to fit their results from uni-axial tests on double edge notched normal and lightweight concrete panels (details of the test set up can be found in Karihaloo [15])
⎛ w = f ⎜⎜ ft ⎝ wc
σ
⎞ w ⎟⎟ − f (1) ; ⎠ wc
⎛ w f ⎜⎜ ⎝ wc
3 w ⎞ ⎡ ⎛ w ⎞ ⎤ −C2 wc ⎟⎟ = ⎢1 + ⎜⎜ C1 ⎟⎟ ⎥ e ⎠ ⎢⎣ ⎝ wc ⎠ ⎥⎦
(1)
Xiao and Karihaloo [11] showed that following polynomial L
⎛
L
⎞
i =1
⎝
i =1
⎠
σˆ y = 1 + ∑ α i wˆ (2 3 )i − ⎜1 + ∑ α i ⎟ wˆ (2 3 )( L +1)
(2)
not only fits the experimental results of Cornelissen et al. [14], but also many commonly-used cohesion-separation laws, e.g., rectangular, linear, and bilinear. In (1), σ and ft are the stress normal to the cohesive crack face and the uni-axial tensile strength, respectively; w and wc are the opening
displacement of the cohesive crack faces, and the critical opening displacement of the pre-existing ˆ = w wc . C1, C2 and αi, i = 1 ∼ L, macro-crack tip when it begins to grow. In (2), σˆ y = σ y f t and w are fitting parameters and L is an integer. The choice of the special form of the cohesion-separation relation (2) involving as it does fractional ˆ may seem strange but it ensures that the asymptotic fields at the cohesive crack tip are powers of w separable into radial and angular variations exactly as the Williams expansions at a traction-free crack tip. Karihaloo and Xiao [12] showed that separable asymptotic crack tip fields are also obtainable for the ˆ cohesion-separation relation not involving fractional powers of w L
⎛
L
⎞
i =1
⎝
i =1
⎠
σˆ y = 1 + ∑ α i wˆ 2i − ⎜1 + ∑ α i ⎟ wˆ 2( L +1)
(3)
where L is an integer. In [11,12] it is shown that Eqs (2,3) capture many commonly used cohesionseparation laws. For the high-performance fibre-reinforced concrete, CARDIFRC [16], tensile tests showed that the post-peak cohesion-separation curve can be accurately represented by the polynomial
σˆ y = 116.76 wˆ 7 - 468.08 wˆ 6 + 738.12 wˆ 5 – 612.10 wˆ 4 + 272.71 wˆ 3 – 59.33 wˆ 2 + 2.89 wˆ + 1.00 The above cohesion-separation law can be easily expressed in the general form (2) or (3). Note that for the CARDIFRC material, ft = 16 MPa and wc = 6.5 mm.
MATHEMATICAL PRELIMINARIES The mathematical formulation is described in [11,12], so that only a brief overview will suffice here. For plane problems, the stresses and displacements in the Cartesian coordinate system (see, Fig. 1) can be expressed in terms of two analytic functions φ(z) and χ(z) of the complex variable z = reiθ
σ x + σ y = 2[φ ' ( z ) + φ ' ( z )] σ y − σ x + 2iτ xy = 2[ zφ '' ( z ) + χ '' ( z )]
(4)
2μ (u + iv) = κφ ( z ) − zφ ( z ) − χ ( z ) '
'
where a prime denotes differentiation with respect to z and an overbar the complex conjugate. In (4), μ = E 2(1 + ν ) is the shear modulus; the Kolosov constant is κ = 3 − 4ν for plane strain or
[
]
κ = (3 − ν ) (1 + ν )
for plane stress; E and ν are Young’s modulus and Poisson’s ratio, respectively.
For a general plane mixed mode I + II problem, the complex functions φ(z) and χ(z) can be chosen as series of complex eigenvalue Goursat functions
φ ( z ) = ∑ An z λ = ∑ An r λ e iλ θ n
n =0
n
n =0
n
χ ( z ) = ∑ Bn z λ n =0
n +1
= ∑ Bn r λn +1e i ( λn +1)θ
(5)
n =0
where the complex coefficients are An = a1n + ia2n and Bn = b1n + ib2n. The eigenvalues λn and coefficients a1n, a2n, b1n and b2n are real.
y ft
σy(w) w
r
θ KI = 0
x
traction-free cohesive zone lp crack Fig. 1: A traction-free crack terminating in a cohesive (or fracture process) zone (FPZ) with residual stress transfer capacity σy(w) whose faces close smoothly near its tip (KI = 0). The material outside the FPZ is linear elastic, but within the FPZ is softening.
Substitution of the complex functions (5) into (4) gives the complete series expansions of the displacements and stresses near the tip of the crack, for details, see [11,12]. These solutions need to satisfy the proper symmetry conditions along the line of extension ahead of the cohesive crack, and boundary conditions on the cohesive crack faces. Considering frictional cohesive crack faces, the boundary conditions on the cohesive crack (0≤ r ≤ lp ) are
σy
θ =π
=σy
θ = −π
≠ 0 , τ xy
θ =π
= τ xy
θ = −π
= −μ f σ y
θ = ±π
≠0
(6)
where μf equals the positive or negative value of the coefficient of kinetic friction, which is assumed to be constant, depending on the relative sliding direction of the two crack faces. Specifically, μf > 0 when δ > 0 and μf < 0 when δ < 0. The length of the process (cohesive) zone lp is either prescribed (i.e. an initial cohesive zone exists before the loading is applied, and does not propagate under the ˆ = 1 in the normal cohesion-separation relation present loading) or is determined by the condition w (2) or (3) at the instant of growth of the pre-existing traction-free crack. ASYMPTOTIC FIELDS AT THE TIP OF A MIXED-MODE COULOMB FRICTIONAL COHESIVE CRACK In principle, a cohesive relationship can also be considered in the tangential direction for quasi-brittle materials. However, this is a contentious issue, since it is difficult to separate the cohesive-sliding relation from the frictional force between the rough cohesive crack faces. Hence, we consider the Coulomb friction between the crack faces instead of a tangential cohesive relationship. With the assumption that (for details, see [11, 12])
μ f a1n − a 2 n ≠ 0
(7)
the separated crack tip asymptotic solutions can be obtained after satisfying the boundary conditions (6), as shown in [11,12]. The complete asymptotic solutions are composed of two parts, corresponding to integer and non-integer eigenvalues, respectively. The expansions corresponding to integer eigenvalues -- giving non-vanishing cohesive stresses – are
λn = n + 1 ,
n = 0, 1, 2, L
(8)
{
2 μ u = ∑ r n +1 a1In ⎡⎣κ cos ( n + 1) θ − ( n + 1) cos ( n − 1) θ − ( n + 2 ) μ f sin ( n + 1) θ ⎤⎦ n =0
+ a2I n ⎡⎣ − (κ + n ) sin ( n + 1) θ + ( n + 1) sin ( n − 1) θ ⎤⎦
(9)
}
− ( n + 2 ) b1In ⎡⎣ cos ( n + 1) θ + μ f sin ( n + 1) θ ⎤⎦
{
2 μ v = ∑ r n +1 a1In ⎡⎣κ sin ( n + 1) θ + ( n + 1) sin ( n − 1) θ − ( n + 2 ) μ f cos ( n + 1) θ ⎤⎦ n=0
+ a2I n ⎡⎣(κ − n ) cos ( n + 1) θ + ( n + 1) cos ( n − 1) θ ⎤⎦
(10)
}
+ ( n + 2 ) b1In ⎡⎣sin ( n + 1) θ − μ f cos ( n + 1) θ ⎤⎦
{
σ x = ∑ r n ( n + 1) a1In ⎡⎣ 2 cos nθ − n cos ( n − 2 ) θ − ( n + 2 ) μ f sin nθ ⎤⎦ n =0
+ ( n + 1) a2I n ⎡⎣ − ( n + 2 ) sin nθ + n sin ( n − 2 ) θ ⎤⎦
(11)
}
− ( n + 1)( n + 2 ) b1In ( cos nθ + μ f sin nθ )
{
σ y = ∑ r n ( n + 1) a1In ⎡⎣ 2 cos nθ + n cos ( n − 2 ) θ + ( n + 2 ) μ f sin nθ ⎤⎦ n=0
+ ( n + 1) a2I n ⎡⎣( n − 2 ) sin nθ − n sin ( n − 2 ) θ ⎤⎦
(12)
}
+ ( n + 1)( n + 2 ) b1In ( cos nθ + μ f sin nθ )
{
τ xy = ∑ r n ( n + 1) a1In ⎡⎣ n sin ( n − 2 ) θ − ( n + 2 ) μ f cos nθ ⎤⎦ n =0
+ n ( n + 1) a2I n ⎡⎣cos ( n − 2 ) θ − cos nθ ⎤⎦
(13)
}
+ ( n + 1)( n + 2 ) b1In ( sin nθ − μ f cos nθ )
The superscript “I” distinguishes coefficients associated with integer eigenvalues. The expansions corresponding to non-integer (fractional) eigenvalues -- giving non-vanishing crack opening -- are
λn =
2n + 3 , 2
2μu = ∑ r n =0
2 n +3 2
n = 0, 1, 2, L ⎧ f ⎡⎛ 2n + 1 ⎞ 2n + 3 2n + 3 2n − 1 ⎤ θ− θ⎥ cos ⎨a1n ⎢⎜ κ + ⎟ cos 2 ⎠ 2 2 2 ⎦ ⎩ ⎣⎝
2n + 5 ⎞ 2n + 3 2n + 3 2n − 1 ⎤ ⎫ ⎡ ⎛ θ+ θ ⎥⎬ + a2fn ⎢ − ⎜ κ + sin ⎟ sin 2 ⎠ 2 2 2 ⎣ ⎝ ⎦⎭
(14)
(15)
2μ v = ∑ r
2 n +3 2
n=0
⎧ f ⎡⎛ 2n + 1 ⎞ 2n + 3 2n + 3 2n − 1 ⎤ θ+ θ⎥ sin ⎨ a1n ⎢⎜ κ − ⎟ sin 2 ⎠ 2 2 2 ⎦ ⎩ ⎣⎝
2n + 5 ⎞ 2n + 3 2n + 3 2n − 1 ⎤ ⎫ ⎡⎛ θ+ θ ⎥⎬ + a2fn ⎢⎜ κ − cos ⎟ cos 2 ⎠ 2 2 2 ⎣⎝ ⎦⎭
σx = ∑r
2 n +1 2
n =0
⎡ 2n + 3 f ⎛ 2n + 5 2n + 1 2n + 1 2n − 3 ⎞ ⎢ 2 a1n ⎜ 2 cos 2 θ − 2 cos 2 θ ⎟ ⎝ ⎠ ⎣
2n + 3 f ⎛ 2n + 9 2n + 1 2n + 1 2n − 3 ⎞ ⎤ θ+ θ ⎟⎥ a2 n ⎜ − sin sin + 2 2 2 2 2 ⎝ ⎠⎦
σy = ∑r
2 n +1 2
n=0
⎡ 2n + 3 f ⎛ −2n + 3 2n + 1 2n + 1 2n − 3 ⎞ ⎢ 2 a1n ⎜ 2 cos 2 θ + 2 cos 2 θ ⎟ ⎝ ⎠ ⎣ 2n + 1 2 n + 3 f ⎛ 2n + 1 2n − 3 ⎞ ⎤ θ − sin θ ⎟⎥ + a2 n ⎜ sin 2 2 2 2 ⎝ ⎠⎦
τ xy = ∑ r n =0
2 n +1 2
⎡ 2n + 1 2n + 3 f ⎛ 2n − 3 2n + 1 ⎞ θ − sin θ⎟ a1n ⎜ sin ⎢ 2 2 2 2 ⎝ ⎠ ⎣
(16)
(17)
(18)
(19)
2n + 3 f ⎛ 2n + 1 2n − 3 2n + 5 2n + 1 ⎞ ⎤ θ− θ ⎟⎥ + a2 n ⎜ cos cos 2 2 2 2 ⎝ 2 ⎠⎦ The superscript “f” distinguishes coefficients associated with non-integer (fractional) eigenvalues.
If
μf = 0,
(8) – (19) reduce to the asymptotic crack tip field of a frictionless cohesive crack with
normal cohesive separation. Furthermore, we can make the following remarks on the relationship between the asymptotic fields of frictionless cohesive cracks and those of the traction free cracks (Williams’ expansions): 1.
For integer eigenvalues, the terms in (9)- (13) corresponding to traction free cracks, but the terms associated with cracks,
a 2I n are identical to those for
a1In and b1In are different; in traction free
b1In are dependent on a1In ; while in frictionless cohesive cracks a1In and b1In are
independent, resulting in non-vanishing normal cohesive tractions on the crack faces. 2.
The expansions (15)-(19) corresponding to non-integer eigenvalues are the same as for traction free cracks; only the first singular terms corresponding to the eigenvalue 1/2 in traction free cracks have been excluded from frictionless cohesive cracks, since there is no singularity at their tips.
3.
For cracks opened by crack face loadings, the above asymptotic fields for frictionless cohesive cracks can be readily used with the inclusion of the singular terms corresponding to the eigenvalue 1/2.
The displacements corresponding to the1st integer eigenvalue are
{ 2 μ v = r {a
} cos θ )}
I 2μ u = r a10I ⎡⎣(κ − 1) cos θ − 2μ f sin θ ⎤⎦ − a20 (κ + 1) sin θ − 2b10I ( cos θ + μ f sin θ ) I 10
I ⎡⎣(κ − 1) sin θ − 2μ f cos θ ⎤⎦ + a20 (κ + 1) cos θ +2b10I ( sin θ − μ f
(20)
The displacements corresponding to the 2nd integer eigenvalue are
{
I 2 μ u = r 2 a11I (κ cos 2θ − 2 − 3μ f sin 2θ ) + a21 ⎡⎣ − (κ + 1) sin 2θ + 2 ⎤⎦
}
−3b11I ( cos 2θ + μ f sin 2θ )
{
I 2 μ v = r 2 a11I (κ sin 2θ + 2 − 3μ f cos 2θ ) + a21 ⎡⎣(κ − 1) cos 2θ + 2 ⎤⎦
}
+3b11I ( sin 2θ − μ f cos 2θ )
(21)
The displacements corresponding to the1st non-integer (fractional) eigenvalue are 3 ⎧ ⎡⎛ ⎡ ⎛ 1⎞ 3 3 1 ⎤ 5⎞ 3 3 1 ⎤⎫ 2 μ u = r 2 ⎨a10f ⎢⎜ κ + ⎟ cos θ − cos θ ⎥ + a20f ⎢ − ⎜ κ + ⎟ sin θ − sin θ ⎥ ⎬ 2⎠ 2 2 2 ⎦ 2⎠ 2 2 2 ⎦⎭ ⎣ ⎝ ⎩ ⎣⎝
3 ⎧ ⎡⎛ ⎡⎛ 1⎞ 3 3 1 ⎤ 5⎞ 3 3 1 ⎤⎫ 2 μ v = r 2 ⎨a10f ⎢⎜ κ − ⎟ sin θ − sin θ ⎥ + a20f ⎢⎜ κ − ⎟ cos θ + cos θ ⎥ ⎬ 2⎠ 2 2 2 ⎦ 2⎠ 2 2 2 ⎦⎭ ⎣⎝ ⎩ ⎣⎝
(22)
The displacements corresponding to the 2nd non-integer eigenvalue are
⎧ ⎡⎛ ⎡ ⎛ 3⎞ 5 5 1 ⎤ 7⎞ 5 5 1 ⎤⎫ 2 μ u = r ⎨a11f ⎢⎜ κ + ⎟ cos θ − cos θ ⎥ + a21f ⎢ − ⎜ κ + ⎟ sin θ + sin θ ⎥ ⎬ 2⎠ 2 2 2 ⎦ 2⎠ 2 2 2 ⎦⎭ ⎣ ⎝ ⎩ ⎣⎝ 5 2
⎧ ⎡⎛ ⎡⎛ 3⎞ 5 5 1 ⎤ 7⎞ 5 5 1 ⎤⎫ 2 μ v = r ⎨ a11f ⎢⎜ κ − ⎟ sin θ + sin θ ⎥ + a21f ⎢⎜ κ − ⎟ cos θ + cos θ ⎥ ⎬ 2⎠ 2 2 2 ⎦ 2⎠ 2 2 2 ⎦⎭ ⎣⎝ ⎩ ⎣⎝ 5 2
(23)
For integer eigenvalues the opening displacement (COD) behind the cohesive zone tip and the sliding displacement of the cohesive crack faces vanish
w = v θ =π − v θ = −π = 0 ,
δ = u θ =π − u θ = −π = 0 , while σy and τxy are non-zero along the cohesive crack faces with
(24)
σˆ y =
σy ft
= ∑ cn r n = 1 + ∑ cn r n n =0
θ = ± π , 0≤ r ≤l p
(25)
n =1
where
(n + 1)(n + 2)(a1In + b1In )cos nπ
cn =
(26)
ft
with
2(a10I + b10I ) c0 = =1 ft since
σ y θ = ±π = f t
(27)
at r = 0.
For non-integer eigenvalues
σ y θ = ±π = 0 , wˆ =
δ
w wc
0≤ r ≤l p
dn =
Denote
(28)
= ∑ dnr 0≤ r ≤ l p
=−
2 n +3 2
n =0
κ +1 ∑ a2fn r μ n =0
3
= d0r 2 + ∑ dnr
2 n +3 2
(29)
n =1
2 n +3 2
sin
2n + 3 π 2
2n + 3 κ +1 f a1n sin π 2 μwc
(30)
(31)
d 0 = d 0 , d n = d n d 0 (n ≥ 1) ,
d0 = −
af κ +1 f 2n + 3 a10 , d n = − 1fn sin π 2 μwc a10
(32)
By enforcing the cohesion-separation law (2) or (3), non-linear relationships between dn and cn are obtained. The general relationships can be found in [11,12]. In [17] these general relationships have been specialized to two types of cohesion-separation law, namely
σˆ y = 1
(33)
and
σˆ y = 1 − wˆ 2 ( L +1)
,
L = 1, 2, 3, L
(34)
The cohesion-separation law (33) is akin to the Dugdale [18] model of plasticity in thin metallic sheets, which is also frequently used in fibre-reinforced polymeric materials, whereas the law (34) is commonly used in studying the fracture process zone embedded within the large plastic zone in metals.
IMPLEMENTATION OF THE ASYMPTOTIC FIELDS IN XFEM To model the cohesive cracks in the XFEM [3-6], a standard local FE displacement approximation around the crack is enriched with discontinuous Heaviside functions along the crack faces behind the crack tip including the open traction-free part, and the crack tip asymptotic displacement fields at nodes surrounding the cohesive crack tip using the PU. The approximation of displacements for an element can be expressed in the following form
⎧u h ( x)⎫ ⎧u 0i ⎫ ⎨ h ⎬ = ∑ φi ( x)⎨ ⎬ + ⎩v0i ⎭ ⎩v ( x) ⎭ i∈I
⎧ b1 j ⎫ ⎧u ⎫ φ j ( x) H ( x)⎨ ⎬ + ∑ φ m ( x)⎨ m ⎬ ∑ j∈J ∩ I ⎩v m ⎭ ⎩b2 j ⎭ m∈M k ∩ I
( tip k )
(35)
Where I is the set of all nodes in the element, (u0i, v0i) are the regular degrees of freedom at node i, φi is the FE shape function associated with node i, J is the subset of nodes whose support is intersected by the crack but do not cover any cohesive crack tips, the function H(x) is the Heaviside function centred on the crack discontinuity, and (b1j, b2j) are the corresponding additional degrees of freedom. Mk is the subset of nodes that are enriched around the cohesive crack tip k with the asymptotic displacements u(tip k) and v(tip k). The general mixed-mode asymptotic field for each integer eigenvalue has three independent terms, e.g. (20)-(21), and each non-integer (fractional) eigenvalue has two independent terms, e.g. (22)-(23). The fields for integer and non-integer eigenvalues need to be used together to produce the crack opening and non-vanishing tractions. In order to improve the accuracy of stresses, we may use eight terms corresponding to one noninteger eigenvalue and two integer eigenvalues as crack tip enrichment functions in (35)
⎧u ⎫ ⎨ ⎬ ⎩v ⎭
(tip )
= Φ(r ,θ )q (tip )
(36)
where q(tip) ⎯ [ a10 , a20 , b10 , a10 , a20 , a11 , a21 , b11 ]T are additional nodal degrees of freedom at the I
I
I
f
f
I
I
I
enriched nodes to be solved together with nodal displacements of conventional FEM. Φ in (36) is the matrix formed by (r, θ) terms dependent on the additional unknown coefficients in (20)(22). For a mode I frictionless cohesive crack, the leading asymptotic displacement terms which correspond to a non-integer eigenvalue adopted in [13] that gives a normal displacement discontinuity over the cohesive-crack faces can be used as crack tip enrichment functions:
u= v=
r 3 2 ⎡⎛ 1⎞ 3 3 1 ⎤ a1 ⎢⎜ κ + ⎟ cos θ − cos θ ⎥ 2μ ⎣⎝ 2⎠ 2 2 2 ⎦
(37)
r 3 2 ⎡⎛ 1⎞ 3 3 1 ⎤ a1 ⎢⎜ κ − ⎟ sin θ − sin θ ⎥ 2 μ ⎣⎝ 2⎠ 2 2 2 ⎦
(38)
There are two possible ways to implement the above expansions in the XFEM. The first is to consider directly the non-linear relationships between dn and cn and to obtain the cohesive stresses from the expansions. This avoids iterations on the cohesion-separation law, but requires solution of a system of algebraic equations with nonlinear constraints. The second way is to treat cn and dn as independent variables, and satisfy the cohesion-separation law (2) or (3) iteratively. The cohesive stresses in this case are obtained from the cohesion-separation law but not the expansions. However, these
expansions can be used to smooth the numerically computed results. This way is generally more convenient to implement and will be used below.
EXAMPLES, RESULTS AND DISCUSSION
In this section, we will analyze a typical mode I cohesive cracking problem of plain and fibrereinforced concrete using the above asymptotic fields and XFEM. We consider three point bend beams without any initial crack (Fig. 2) made of plain concrete with linear, bilinear and constant softening laws, and of CARDIFRC. We have actually assumed a very small initial crack of length 0.1mm at the bottom midpoint of the beam. A state of plane strain is considered for all specimens. The geometrical parameters l = 4b, t = b (t is the specimen thickness in the out-of-plane direction). As in Moës and Belytschko [10], x-direction of nodes with coordinates (0, 0) and (0, l) and y-direction of the node with coordinates (b, l/2) are constrained; the central “point” load is distributed over two elements. The potential fracture locus coincides with the beam axis of symmetry (y = l/2 = 2b). In XFEM, the crack is modelled by enriching the nodes on the crack faces with jump and branch functions without the double nodes that are used in the traditional FEM [19]. The conventional 4-node bilinear isoparametric Q4 elements are used as background elements. The first layer of nodes surrounding the cohesive crack tip (the elements that include the crack tip k are defined as the first layer elements of the crack tip with enriched nodes; the nodes in the first layer elements are called the first layer enriched nodes) are enriched with asymptotic displacements. Other nodes on the crack faces are enriched with jump functions.
P
b
x y l P/2
Fig. 2: Three-point bend beam (TPB).
P/2
600
550
550
500
500
450
450
y
y
600
400
400
350
350
300
300
0
50
x
100
150
0
50
x
100
150
(a) (b) Fig. 3: Coarse (a) and fine (b) mesh for one half of the specimen to the left of mid-span (Fig. 2).
We consider first the three-point bend plain concrete beams with b = 150mm. Two meshes, as shown in Fig. 3, are used in the analysis. The coarser mesh consists of 50×100 = 5000 rectangular elements, giving a total of 5151 nodes. The finer mesh consists of 150×120 = 18000 rectangular elements, giving a total of 18271 nodes. Both meshes are uniformly divided in x-direction. For the coarser mesh, the central 50 layers of elements have an identical height (y-direction) of 3mm; the remaining elements have an identical height of 9mm. Therefore, elements in the central zone are 3×3 mm2 squares. For the finer mesh, the central 60 layers of elements have an identical height of 1mm; the remaining elements have an identical height of 9mm. Therefore elements in the central zone are 1×1 mm2 squares. The finer mesh is only used for linear and bilinear softening to investigate whether the results have any mesh sensitivity. In the simulation, the first increment of the cohesive crack is 4.4mm; thereafter the cohesive crack propagates by a segment of length 3mm after each step in the coarser mesh, and by three segments of length 1mm each in the finer mesh. Three-point bend plain concrete beam with linear softening The initiation and growth of a cohesive crack in such a beam has been studied extensively by Carpinteri and Colombo [19] using the FEM and the node release technique. The material properties with linear softening law (Fig. 4) are
f t = 3.19MPa, GF = 50 Nm −1
E = 36.5GPa, ν = 0.1,
where E is Young’s modulus, ν Poisson’s ratio, and GF the specific fracture energy, and
σ y = f t (1 −
w ) wc
1 0.8 0.6
σy /f t
Linear law Bilinear law
0.4 0.2
(w 1 /w c , f 1 /f t )
0 0
0.2
0.4 0.6 w /w c
0.8
1
Fig. 4: Linear and bilinear softening laws.
Since a low Poisson’s ratio of 0.1 is used, the results with the load being distributed over two elements are believed to be close to [19] where a plane stress condition is assumed and a concentrated load was considered. Asymptotic displacements (37) and (38) are used as crack tip enrichment functions. The nondimensional load-midspan deflection curves are shown in Fig. 5. They agree very well with the results of Carpinteri and Colombo [19]. As the specific fracture energy is rather moderate (GF = 50 Nm-1), the snap-back in the global load-deflection curve (Fig. 5) is not pronounced.
0.25
A Coarse mesh
Load / (f tbt )
0.2
B
0.15
Fine mesh Carpinteri & Colombo (1989)
0.1 0.05
C
0 0
5 10 15 4 Deflection / b × 10
20
Fig. 5: The non-dimensional load-midspan deflection curves of the three-point bend beam
Three-point bend plain concrete beam with bilinear softening We now report the results of the above beam with the properties and a bilinear softening law obtained for a real normal strength concrete, studied in [13]. The material properties are E = 36.9 GPa, ν = 0.2, and the bilinear softening law (Fig. 4) is described by
w( f t − f1 ) ⎧ ], 0 ≤ w < w1 ⎪⎪ f t [1 − w f 1 t σy =⎨ f 1 wc w ⎪ (1 − ), w1 ≤ w ≤ wc ⎪⎩ ( wc − w1 ) t wc
(39)
with
f t = 3.14MPa, f1 = 0.455MPa, wc = 0.279mm, w1 = 0.0373mm, G F = 122 N / m As for the linear softening law, asymptotic displacements (37) and (38) are used as crack tip enrichment functions. The non-dimensional load-midspan deflection curves are shown in Fig. 6. The evolution of the cohesive zone size as the cohesive tip travels through the beam is shown in Fig. 7. These global properties reveal very weak mesh size sensitivity.
0.3 A
Load/(f tbt )
0.25
Coarse mesh Fine mesh
B
0.2 0.15
C
0.1 0.05 0 0
5 10 15 4 Deflection/b ×10
20
Fig. 6: The non-dimensional load-midspan deflection curves of the three-point bend beam with a bilinear softening law.
1 Coarse mesh
FPZ / b
0.8
Fine mesh B
0.6 0.4
(1)
C
A
(2)
0.2 0 0
0.2 0.4 0.6 0.8 Total crack length / b
1
Fig. 7: Evolution of the cohesive zone size as the cohesive tip travels through the beam. (1) and (2) correspond to the first and second branches of the bilinear softening diagram (Fig. 4).
Three-point bend plain concrete beam with constant softening As the leading asymptotic displacement terms from the two cohesion-separation laws (33) and (34) are the same [17], we reanalyse the three point bend beam (TPB) analysed above under the bilinear softening law using the constant traction law (33) with
Load/(ftbt )
ft = 3.14MPa, wc = 0.03886mm, GF = 122 N / m
0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0
Bilinear Constant
0
5
10 Deflection/b ×10
15
20
4
Fig. 8: The non-dimensional load-midspan deflection curves of the TPB with constant traction and bilinear softening laws.
Furthermore, we use (20) - (22) as crack tip enrichment functions instead of the leading asymptotic displacement terms of a mode I frictionless cohesive crack (37) and (38) adopted in [13]. The nondimensional load-midspan deflection curve is shown in Fig. 8 and compared with that of the bilinear law (Fig. 6). It is clear that although the tensile strength ft and specific fracture energy GF are the same, the initial branch of the cohesion-separation law produces very different responses: the constant traction law produces a higher peak load but a more brittle overall response with a significant snapback. Three point bend fibre reinforced concrete beam We now consider the smaller three point bend beam b = 100mm made of CARDIFRC. We found it convenient to approximate the polynomial form of its cohesion separation relationship given earlier by the bilinear relationship (39) but now the material parameters are E = 48GPa, ft = 16MPa, w1 = 2.29mm, f1 = 3.9088MPa and Poisson’s ratio ν = 0.2.The constraint and loading conditions are identical to the plain concrete beam. The discretised mesh is also identical to that shown in Fig. 3a after scaling the coordinates by a factor of 2/3. Asymptotic displacements (37) and (38) are used as crack tip enrichment functions. In the simulation, the first increment of the cohesive crack is 2.9mm, thereafter the cohesive crack propagates by a segment of length 2mm after each step. When the cohesive crack is developed over the entire depth of the specimen at a deflection of about 3mm, jump functions are used to enrich all nodes on the crack faces. The simulation is continued by increasing the deflection by 0.5mm in each step.The load-midspan deflection curves are shown in Fig. 9. The overall shape of the curve, peak load and the corresponding displacement from the XFEM simulation using the bilinear law agrees very well with the experiments. The effect of large deformation has been ignored in the present simulation, which may be significant at the final stage. If this effect is taken into account, the agreement may be further improved. 70.00
60.00
Load (kN)
50.00
40.00 HM 1 Test 1 30.00
20.00
10.00
0.00 0.0000
1.0000
2.0000
3.0000
4.0000
5.0000
6.0000
7.0000
Displacement (mm)
(a) Load –deflection curves from experiments and the hinge model
8.0000
9.0000
10.0000
70 60
Load (kN)
50 40 30 20 10 0 0
1
2
3
4
5
6
7
8
9
10
Loading point deflection (mm) (b) Load –deflection curve from the XFEM simulation Fig. 9. The load-midspan deflection curves of the CARDIFRC three-point bend beam.
Conclusion The solved model I problems have revealed the efficiency of the XFEM enriched with true crack tip asymptotic displacements in simulating cracking process of plain and fibre reinforced concrete. The method may be more advantageous in simulating more complicated mix-mode cracking process.
REFERENCES [1] Hillerborg A., Modeer E., Petersson P.E.: Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research 1976; 6: 773-781. [2] de Borst R., Gutierrez M.A., Wells G.N., Remmers J.J.C., Askes H.: Cohesive-zone models, higher-order continuum theories and reliability methods for computational failure analysis. International Journal for Numerical Methods in Engineering 2004; 60: 289-315.
[3] Moës N., Dolbow J., Belytschko T.: A finite element method for crack growth without remeshing. International Journal of Numerical Methods in Engineering 46 (1999) 131-150. [4] Strouboulis T., Copps K., Babuška I.: The generalized finite element method. Computer Methods in Applied Mechanics and Engineering 190 (2001) 4081-4193. [5] Babuška I., Banerjee U., Osborn J.E.: Survey of meshless and generalized finite element methods: A unified approach. Acta Numerica 12 (2003) 1–125. [6] Karihaloo B.L., Xiao Q.Z.: Modelling of stationary and growing cracks in FE framework without remeshing: a state-of-the-art review. Comp Struct 81 (2003) 119-129. [7] Wells G.N., Sluys L.J.: A new method for modelling cohesive cracks using finite elements. International Journal for Numerical Methods in Engineering 2001; 50: 2667-2682. [8] Hansbo A., Hansbo P.: A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Computer Methods in Applied Mechanics and Engineering 2004; 193: 3523-3540. [9] Zi G., Belytschko T.: New crack-tip elements for XFEM and applications to cohesive cracks. International Journal for Numerical Methods in Engineering 2003; 57: 2221-2240. [10] Moës N., Belytschko T.: Extended finite element method for cohesive crack growth. Engineerig Fracture Mechanics 2002; 69: 813-833. [11] Xiao Q.Z., Karihaloo B.L.: Asymptotic fields at frictionless and frictional cohesive crack tips in quasi-brittle materials. Journal of Mechanics of Materials and Structures. 1, 2006, 881-910. [12] Karihaloo B.L., Xiao Q.Z.: Asymptotic fields at the tip of a cohesive crack. International Journal of Fracture, 150, 2008, 55-74. [13] Xiao Q.Z., Karihaloo B.L., Liu X.Y.: Incremental–secant modulus iteration scheme and stress recovery for simulating cracking process in quasi-brittle materials using XFEM. International Journal for Numerical Methods in Engineering. 69, 2007, 2606-2635. [14] Cornelissen H.A.W., Hordijk D.A., Reinhardt H.W.: Experimental determination of crack softening characteristics of normal and lightweight concrete. Heron 31 (1986) 45-56 [15] Karihaloo B.L.: Fracture Mechanics and Structural Concrete. Addison Wesley Longman, UK, 1995
[16] Benson S.D.P., Karihaloo B.L.: CARDIFRC – Development and Mechanical Properties. Part III: Uniaxial Tensile Response and other Mechanical Properties”, Magazine of Concrete Research, 57, 2005, 433—443 [17] Karihaloo B.L., Xiao Q.Z.: Asymptotic fields ahead of mixed mode frictional cohesive cracks, Zeitschrift für Angewandte Mathematik und Mechanik (in press) [18] Dugdale D.S.: Yielding of steel sheets containing slits. Journal of the Mechanics and Physics of Solids 8 (1960)100108 [19] Carpinteri A., Colombo G.: Numerical analysis of catastrophic softening behaviour (snap-back instability). Computers and Structures 31 (1989) 607-636.
