In this paper, a procedure for the implementation of the XFEM within the FE
code Cast3M for two-dimensional fracture problems is presented. The quadrangular finite
element CRCOH for cohesive crack are constructed and involved in to Cast3M for quasibrittle materials. This finite element is enriched by the cohesive stress which is related
with a separation of the crack sides through the law of Hillerborg. The obtained results
are compared with the results of the linear elastic crack showing that the existence of cohesive stress in quasi-brittle materials reduces the separation of crack and stress at crack-tip.
So, using the above mentioned model with cohesive crack and given finite element CRCOH
is necessary. In the near future, relied upon this framework, the growth of cohesive crack
will be continued to develop.

10 trang |

Chia sẻ: huongthu9 | Lượt xem: 470 | Lượt tải: 0
Bạn đang xem nội dung tài liệu **Cast3m implementation of the extended finite element method for cohesive crack**, để tải tài liệu về máy bạn click vào nút DOWNLOAD ở trên

Vietnam Journal of Mechanics, VAST, Vol. 33, No. 1 (2011), pp. 55 – 64
CAST3M IMPLEMENTATION OF THE EXTENDED
FINITE ELEMENT METHOD FOR COHESIVE CRACK
Nguyen Truong Giang, Ngo Huong Nhu
Institute of Mechanics, VAST
Abstract. In this paper, the finite element for cohesive crack for quasi-brittle materials
is constructed by XFEM, which allows incorporation of the displacement discontinuities
in the element. The algorithm of construction and procedures for involving this finite
element into code Cast3M are presented. The numerical calculations in fracture mechan-
ics are presented to demonstrate the benefits of the proposed implementation.
Key words: Fracture, cohesive crack, quasi-brittle material, XFEM, Cast3M.
1. INTRODUCTION
The extended finite element method (XFEM ) was first developed for two-dimensional
linear elastic fracture mechanics [1, 2, 3], then it was used for crack growth, intersecting
cracks and fluid mechanics. In the recent years, XFEM has been also applied to the cohe-
sive crack model [4], [6]-[9]. The model which Hillerborg introduced in [5] has been used
extensively in non-linear fracture mechanics of quasi-brittle materials and is applied in
this paper. There are other developments of XFEM to existing codes. An implementation
of the extended finite element method for linear elastic fracture mechanics problem within
the finite element software Abaqus is introduced in [11]. In the same interest, an object
enriched finite element library, named OpenXFEM++ which has been built and used to
solve various 2D linear elastic fracture problems with great success, is presented in [12].
The finite element for the brittle crack by XFEM is constructed and involved into code
Cast3M [3]. In order to produce extendable code, the objective of this paper is to present
the algorithm of construction of new extended cohesive crack element 4 nodes and the way
to add it into an existing finite element code Cast3M.
2. COHESIVE CRACK FOR QUASI-BRITTLE MATERIALS
The fracture process in quasi-brittle materials like concrete, rocks, ceramics, fiber
reinforced composites etc., is characterized by the contribution of the resistance from
aggregate interlocking, and another inelastic effect. At the crack tip region there exists a
narrow strip, where occur non-linear phenomena, such as plasticity yielding, micro- cracks,
etc., however the rest of the body presents an elastic behaviour (Fig. 1a). These inelastic
located processes can be represented by cohesive models, which simulate the damage region
close to the crack tip. These models use to represent those non-linear phenomena such
56 Nguyen Truong Giang, Ngo Huong Nhu
as the "cohesive forces", which are applied along the "cohesive crack". These forces are
defined by relations between stresses and displacements, the material resistance decreases
proportionally to the increase of the displacements [4].
Fig. 1. Sketch showing the process zone of a cohesive law with smooth crack
closure Cohesive stress t(w) used by Hillerborg [5]
The stress closing or cohesive stress t(w) is applied to the crack surface. These
cohesive stresses depend on the crack opening w, and vary from zero at a characteristic
crack opening wc, to the tensile strength of the material ft at the tip of the crack (Figure.
1b). It is assumed that the cohesive stresses close the crack smoothly, thus even when
the un-cracked material is considered as linear elastic, and often used in the modelling
of cementitious materials, stresses are finite in the un-cracked material at the crack tip.
In other words, during Mode I crack propagation, the stress intensity factor KI is zero
and the condition for crack propagation is that the stress at the crack tip has reached
the tensile strength ft. The zone in which the cohesive stresses are located was originally
called the micro-cracked zone [5], and later the process zone or fracture zone [6].
Mathematical formulation for cohesive crack problem
Consider a domain Ω as shown in Fig. 2, containing a crack Γd. The part of the crack
where a cohesive law acted is denoted by Γcoh. Prescribed tractions F are imposed on the
boundary ΓF whereas prescribed displacements u¯ (assumed to be zero for simplicity) are
imposed on Γu. The stress field inside the domain, σ is related to the external loading F
and the tractions t+, t− in the cohesive zone through the equilibrium equations.
The equilibrium conditions for the quasi-brittle materials having cohesive crack
without body force are [7]:
∇ · σ = 0,
σ · n = F on ΓF ,
σ ·m = −t+ on Γ+
coh
,σ ·m = t− on Γ−
coh
,
(1)
m,n are the normal vectors for ΓF and Γcoh, respectively.
The equilibrium condition across surface Γcoh is :
t ≡ t− = −t+
Cast3M implementation of the extended finite element method for cohesive crack 57
G d
n F
Fig. 2. Geometry of modelled domain and notation
The displacement discontinuity w across Γd can be expressed in terms of the dis-
placement vector u computed on two sides of the discontinuity:
w = u
∣∣∣Γ+
d
− u
∣∣∣Γ−
d
Let U be the space of admissible displacements u in Ω; i.e. such that u = u¯ on Γu,
u possibly discontinuous on Γd and u ∈ C0 (C0 is the space of Sobolev) in Ω \ Γd. By
introducing the test functions v ∈ U0 (the displacements are equal zero on Γu), the weak
form of the equilibrium equations reads [7]:∫
Ω
σ (u) : ε (v) dΩ =
∫
ΓF
F · vdΓ +
∫
Γ
+
d
t+ · vdΓ+
∫
Γ
−
d
t− · vdΓ ∀v ∈ U0∫
Ω
σ (u) : ε (v) dΩ =
∫
ΓF
F · vdΓ−
∫
Γd
t · wdΓ∫
Ω
σ (u) : ε (v) dΩ+
∫
Γd
t · wdΓ =
∫
ΓF
F · vdΓ
(2)
The equation (2) is weak form of (1) and it includes cohesive stresses t on the two
crack faces.
3. THE EXTENDED FINITE ELEMENT METHOD
Partition of unity
A collection of functions ϕi(x), each belonging to a node, defined over a body Ω(x ∈
Ω) forms a partition of unity [1] if:
n∑
i=1
ϕi (x) = 1 (3)
58 Nguyen Truong Giang, Ngo Huong Nhu
n is the number of nodal points. The displacement approximation using partition of unity
can be formed by:
uh (x) =
n∑
i=1
ϕi (x)
(
m∑
α=1
ψα (x) a
α
i
)
(4)
where ψα are enrichment functions, m is the number of enrichment functions.
Extended finite element method
Finite element shape functions Ni are also a partition of unity because
n∑
i=1
Ni (x) = 1
Therefore, finite element shape functions are chosen as the partition of unity func-
tions. From equation (4), we note that the finite element space (ψ1 = 1, ψα = 0(α 6= 1)) is
a subspace of the enriched space uh. The enriched displacement approximation is written
as [2]
uh (x) =
∑
I∈N
NI (x)
uI +H (x) aI︸ ︷︷ ︸
I∈NΓ
+
4∑
α=1
Fα (x) b
α
I︸ ︷︷ ︸
I∈NΛ
(5)
where uI is the nodal displacement vector associated with the continuous part of the finite
element, aI is the nodal enriched degree of freedom vector associated with the Heaviside
function, and bαI is the nodal enriched degree vector associated with the elastic asymptotic
crack tip function. N is the set of all nodes in the mesh; NΓ is the set of nodes whose
shape function support is cut by the crack interior; and NΛ is the set of nodes whose
shape function support is cut by the crack tip (NΓ ∩NΛ = ∅) (see [3]).
The interior of crack is modeled by the generalized Heaviside enrichment function
H , where H takes the value +1 above the crack and -1 below the crack [2]
H (x, y) =
{
1 if y > 0
−1 if y < 0 (6)
The functions Fα are defined as:
{Fα (r, θ)} =
{√
r sin
(
θ
2
)
,
√
r cos
(
θ
2
)
,
√
r sin
(
θ
2
)
cos θ,
√
r cos
(
θ
2
)
cos θ
}
(7)
where (r, θ) are the local polar coordinates at the crack tip. Note that the first function
in (7),
√
r sin (θ/2), is the discontinuous across the crack faces whereas the last three
functions are continuous.
4. DISCRETISED EQUATIONS AND FINITE ELEMENT FOR THE
COHESIVE CRACK
The numerical simulation of the cohesive crack is hereby performed within the ex-
tended finite element method which is presented in the previous section. The displacement
Cast3M implementation of the extended finite element method for cohesive crack 59
field u of the body can be additively decomposed into a continuous part ucont and a dis-
continuous part udisc:
u(x) = ucont(x) + udisc(x). (8)
Inserting the Heaviside function into the discontinuous part, the displacement field
in element where the extra degrees of freedom are acted can be written as
u(x) = N (x)a+HN (x)b,
vectors a has the standard degrees and b has the extra degrees.
The displacement discontinuity across Γd can be given
w = 2N (x)b (9)
The strain field and stress field in element can be expressed
ε = Ba+HBb
σ =D(Ba+HBb)
(10)
where D relates the stress and strain.
If T is called the relation between the traction and crack displacement, the traction
can be written as
t = T .w or
{
tn
ts
}
= [T ]
{
wn
ws
}
(11)
Here, we used the cohesive law T which proposed by Hillerborg [5], the normal
traction force tn transmitted across a discontinuity is made by an exponentially decaying
function of the history parameter:
tn = ft exp
(
− ft
Gf
k
)
,
where ft is the tensile strength of the material and Gf is the fracture energy and k (internal
variable) is the largest value of the normal separation wn. The crack shear stiffness is
also made by a function of the history parameter. The shear traction ts acting on the
discontinuity surface is calculated from
ts = dint exp (hsk)ws,
where dint is the initial crack shear stiffness (k= 0), ws is the crack sliding displacement
and hs is equal to:
hs = ln (dk=1.0/dint)
where dk=1.0 is crack shear stiffness when k = 1.0.
Replace the displacement field (9), (10), (11) into the equation (2), the discrete
equations can be obtained as follow:[ ∫
Ω
BTDBdΩ
∫
Ω
HBTDBdΩ∫
Ω
HBTBdΩ
∫
Ω
HBTDBdΩ
]{
a
b
}
+
[
0 0
0 4
∫
Γd
NTTNdΓ
]{
a
b
}
=
∫
ΓF
NTFdΓ
60 Nguyen Truong Giang, Ngo Huong Nhu
The discretised form of virtual work for incremental displacements (da, db) is written
as:[ ∫
Ω
BTDBdΩ
∫
Ω
HBTDBdΩ∫
Ω
HBTBdΩ
∫
Ω
HBTDBdΩ + 4
∫
Γd
NTTNdΓ
]{
da
db
}
=
{
fext
0
}
−
{
f inta
f intb
}
(12)
So, the element matrix can be written:
Ke =
[ ∫
Ω
BTDBdΩ
∫
Ω
HBTDBdΩ∫
Ω
HBTBdΩ
∫
Ω
HBTDBdΩ+ 4
∫
Γd
NTTNdΓ
]
, (13)
where the matrices of shape functions derivatives B for node i are:
B=i
Ni,x 00 Ni,y
Ni,y Ni,x
(14)
f ext =
∫
ΓF
NTFdΓ (15)
f inta =
∫
Ω
BTσdΩ,
f intb =
∫
Ω
HBTσdΩ+ 2
∫
Γd
NT tdΓ
(16)
Numerical integration for the element with the cohesive crack
As in the standard FEM, it is necessary to perform numerical integration over the
element domain for computing the element stiffness matrix (13). The elements which
contain the crack include a displacement discontinuity due to the X-FEM formulation.
These elements must be subdivided into sub-domains in which the crack is one of the
sub-domain boundaries to carry out the numerical integrations. Therefore, the elements
non-crossed by a discontinuity are integrated by standard four-point Gauss quadrature.
When an element is crossed by a discontinuity, the domains Ω+ and Ω− on either side
of the discontinuity are triangulated into sub-domains [3]. Within each triangular sub-
domain, three-point Gauss quadrature is applied. For the cohesive crack sub-domain, two
integration points are positioned on the discontinuity in order to integrate the traction
forces.
5. CAST3M IMPLEMENTATION FOR THE COHESIVE ELEMENT
In this part, a procedure for the implementation of the XFEM within the FE code
Cast3M [10] for two-dimensional fracture problems is presented. The quadrangular element
4 nodes CRCOH is constructed and involved into Cast3M. The implementation was based
on the use of this element and performed in source code of Cast3M by language FORTRAN
and enables in modeling of different crack locations and orientations. In order to develop
XFEM in numerical calculation, a number of subroutines contain XFEM formulation and
some operators represented below such as RIGI for matrix stiffness K of cohesive element,
EPSI for the displacement ε, COMP for stress σ, BSIG for nodal force,. . . are inserted
into software Cast3M .
Cast3M implementation of the extended finite element method for cohesive crack 61
5.1. Algorithm
- Step 0: Initially, a mesh for quadrilateral element, the location of the crack and
the parameters of concrete are inputted.
- Step 1: Using level set to choose crack tip element and crack element. For the
element without crack, the standard FEM is applied.
- Step 2: Calculating of the matrix stiffness K for cohesive elements by using (13)
and then assembling in the matrix stiffness standard.
- Step 3: Applying a given loading step to solve for the displacement field
{
da
db
}
using Eq. (12).
- Step 4: Calculating the incremental deformation dε, and then determinating the
stress σ. Note that the cohesive stresses are local stresses, leading to the co-ordinate system
changing.
- Step 5: Computing the nodal forces through the program BSIG using Eq. (15) and
(16).
- Step 6: Verifying the condition of convergence. If it isn’t satisfied, go to step 3.
On the contrast, if the condition is satisfied, the numerical results as the stress,
deformation, and displacement for each the loading step can be obtained.
5.2. Numerical example and discussion
Example 1: Consider a plate of width w = 100mm and height L = 200mm with an
edge crack of length a = 20mm. The tensile stress σ =10MPa. The concretes parameters
Ec = 35000MPa, ν = 0.3, ft = 3MPa, Gf = 0.1N/mm.
Fig. 3. Comparison of the stresses
62 Nguyen Truong Giang, Ngo Huong Nhu
Fig. 4. Comparison of the deformations
Remark: Separation of crack in the linear case (=3.677x10−2 mm) is larger than
that of cohesive crack (=2.73086x10−2 mm). The stress at a crack-tip (= 45MPa) in
cohesive crack is also represents a contribution of the cohesive stress in comparison to
stress intensity (= 57 MPa) of the linear elastic crack.
Example 2: We use the same dimension of specimen for the incline crack α = 45˚,
length a = 20mm.
Fig. 5. Comparison of the deformations and the σyy for incline crack
Cast3M implementation of the extended finite element method for cohesive crack 63
Remark: In the incline crack case, the influence of the quasi-brittle material on the
crack side displacements is presented (see Fig. 5). Clearly, the stress at a crack-tip of cohe-
sive crack (σmax = 26 MPa) are smaller than that of linear elastic crack (σmax = 47 MPa).
6. CONCLUSIONS
In this paper, a procedure for the implementation of the XFEM within the FE
code Cast3M for two-dimensional fracture problems is presented. The quadrangular finite
element CRCOH for cohesive crack are constructed and involved in to Cast3M for quasi-
brittle materials. This finite element is enriched by the cohesive stress which is related
with a separation of the crack sides through the law of Hillerborg. The obtained results
are compared with the results of the linear elastic crack showing that the existence of cohe-
sive stress in quasi-brittle materials reduces the separation of crack and stress at crack-tip.
So, using the above mentioned model with cohesive crack and given finite element CRCOH
is necessary. In the near future, relied upon this framework, the growth of cohesive crack
will be continued to develop.
ACKNOWLEDGMENTS
The authors are grateful to Professor Alain Millard from French Atomic energy Com-
mission (CEA) for his assistance. This paper is completed with the partly financial support
from the NAFOSTED (National Foundation for Science and TechnologyDevelopment).
REFERENCES
[1] Melenk JM, Babuska I, The partition of unity finite element method: Basis theory and appli-
cation, Computer Methods in Applied Mechanics and Engineering, 45(5) (1999), 601-620.
[2] Moes N, Dolbow J, Belytschko T, A finite element method for crack growth without remeshing,
International Journal for Numerical Methods in Engineering, 46 (1999), 131-150.
[3] Nguyen Truong Giang, Numerical simulation the crack 2D by the finite element incorporated
the discontinuity, Vietnam Journal of Mechanics, VAST, 30(2) (2008).
[4] Estela M.R. Bueno, Tulio N. Bittencourt, Luiz F. Martha, Carlos V.A. Carvalho, A graphics
and object-oriented system for modeling 2D cohesive crack problems, European Congress on
Computational Methods in applied Sciences and Engineering, (2000).
[5] Hillerborg, A., Modéer, M. and Petersson, P.E., Analysis of crack formation and crack growth
concrete by means of fracture mechanics and finite elements, Cements and Concrete Research,
6 (1976), 773-782.
[6] Henrik Stang, John Forbes Olesen, Peter Noe Poulsen, Lars Dick Nielsen, On the application
of cohesive crack modeling in cementious materials, Journal of Materials and Structures, 40
(2007), 365-374.
[7] Stefano Mariani, Umberto Perego, Extended finite element method for quasi-brittle fracture,
International Journal for Numerical Methods in Engineering, 58 (2003), 103-126.
[8] Belytschko T, Black T, Elastic crack growth in finite elements with minimal remeshing, In-
ternational Journal for Numerical Methods in Engineering, 45 (1999), 601-620.
64 Nguyen Truong Giang, Ngo Huong Nhu
[9] Sukumar N, Prevost JH, Modeling quasi-static crack growth with the extended finite element
method, Part I: Computer implementation, International Journal of Solids and Structure, 40
(2003), 7513-7537.
[10]
[11] Giner E, Sukumar N, Tarancon JE, Fuenmayor FJ, An Abaqus implementation of the extended
finite element method, Engineering fracture mechanics, 76 (2009), 347-368.
[12] Bordas S, Phu Vinh Nguyen, Dunant C, Hung Nguyen-Dang, Guidoum A, An extended finite
element library, International Journal for Numerical Methods in Engineering, 2 (2006), 1-33.
Received March 29, 2010

Các file đính kèm theo tài liệu này:

- cast3m_implementation_of_the_extended_finite_element_method.pdf