2nd Reading

December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

International Journal of Computational Methods

Vol. 9, No. 4 (2012) 1250055 (15 pages)

c World Scientific Publishing Company
DOI: 10.1142/S0219876212500557

A TWICE-INTERPOLATION FINITE ELEMENT METHOD

(TFEM) FOR CRACK PROPAGATION PROBLEMS

S. C. WU∗,†,§ , W. H. ZHANG∗ , X. PENG‡ and B. R. MIAO∗

∗StateKey Lab of Traction Power, Southwest Jiaotong University
by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.

Chengdu 230009, China

†State Key Lab of Advanced Design and Manufacturing for Vehicle Body
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

Changsha 410082, China

‡School of Engineering, Institute of Mechanics & Advanced Materials

Cardiff University, Wales, UK

§[email protected]

Received 1 June 2010

Accepted 21 May 2011
Published 31 December 2012

The twice-interpolation finite element method (TFEM) constructs the trial function for
the Galerkin weak form through two stages of sequential interpolation without additional
degrees of freedom and achieves better accuracy and convergence compared to the con-
ventional finite element method (FEM). The TFEM has been shown to be insensitive
to the quality of the elemental mesh and thus has the potential to simulate fracture
problems. Intense examples and issues with cracking problems are investigated in this
study. It is observed that the stress intense factor (SIF) of the crack tip can be evaluated
with satisfactory accuracy. Since the present TFEM can produce more continuous nodal
gradients, better stress fields can be reproduced, especially around the crack, without
requiring more nodes. It is also shown that crack propagation can be reproduced readily
and the cracking path agrees well with the reference solution and experimental results.

Keywords: Crack propagation; twice-interpolation; finite element method; ALOF

software.

1. Introduction
The numerical simulation of crack initiation and its progress is essential for predict-
ing the service life of engineering structures. In the standard finite element method
(FEM), nodal relaxation techniques are frequently used for modeling cracks, where
each node along the free surface is split into two nodes [Swenson and Ingrafea
(1996)]. The partition of unity method [Melenk and Babuška (1996)], which does
not require frequent meshing and remeshing processes, has also emerged as a pow-
erful tool for crack propagation analysis. Moreover, mesh-free methods, such as the

§ Corresponding author.

1250055-1
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

S. C. Wu et al.

element-free Galerkin method (EFG) [Belytschko et al. (1995)] and the meshless
local Petrov–Galerkin method (MLPG) [Atluri and Zhu (1998)], which are more
flexible than the standard FEM, have been developed to solve fracture problems
with arbitrary node distributions. In addition, some other mesh-free approaches
[Rabczuk and Belytschko (2004); Zheng et al. (2008); Liu and Zhang (2009);
Chen et al. (2010); Liu et al. (2010); Liu et al. (2011)] have also been proposed
for progressive cracking problems. However, these well-known reduced mesh meth-
ods often result in high computation costs, and thus, many robust solutions and
other excellent properties developed for the FEM are not utilized. As a special exam-
ple of the partition of unity method, the extended finite element method (XFEM)
by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.

[Sukumar et al. (2000); Belytschko et al. (2009)], which offers advantages over other
numerical methods for discontinuous problems of dynamic fractures, was developed
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

and has been successfully employed in both commercial ALOF [Wu and Wu (2011)]
to solve complex three-dimensional (3D) fracture problems. If only the XFEM is
used to simulate the cracking process, the numerical accuracy cannot always be
assured. The XFEM also requires a finer mesh, especially in the moving crack tip.
By integrating the XFEM and the virtual node method (VNM) [Tang et al. (2009)]
into the ALOF, satisfactory accuracy and lower calculation costs can be achieved
for progressive cracking problems.
Recently, a novel twice-interpolation finite element method (TFEM) [Zheng
et al. (2010)] was proposed for solid mechanics problems. In the TFEM, the trial
function for the Galerkin weak form is constructed through two stages of sequential
interpolation. Similar to the standard FEM, the TFEM is based on an elemental
mesh, with the unknowns being the nodal displacements. In the standard FEM,
the trial functions are obtained by interpolation according to nodal displacement
and are represented as C 0 on the element edges. The first stage of TFEM is the
same as that of the standard FEM, but the average nodal gradients must be com-
puted for the second stage of the interpolation, unlike with the second stage of
the interpolation for the Hermite FEM. Then, the approximation functions can be
further rebuilt, in which both the nodal displacement and average gradients are
selected as the interpolation conditions. Thus, the shape functions constructed by
the twice-interpolation step, exhibit more continuous nodal gradients and higher-
order polynomial contrast compared to the standard FEM when analyzing the same
mesh. Intensive numerical benchmark problems have shown that the TFEM exhibits
better accuracy and convergence properties than standard FEM.
In the present work, the TFEM is further developed for fracture problems into
ALOF. In order to maintain the standard FEM computer procedure to the extent
possible, the nodal relaxation technique is utilized to model the crack growth. In
addition, a node projection technique is formulated to improve the modeling of the
crack evolution. The outline of this paper is as follows. Section 2 briefly reviews the
twice-interpolation strategy, and Sec. 3 details the evaluation of the stress intensity
factor and simulation of the crack growth. An intensive numerical study is conducted
in Sec. 4, and the conclusions are presented in Sec. 5.

1250055-2
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

A Twice-Interpolation Finite Element Method (TFEM ) for Crack Propagation Problems

2. Twice-Interpolation Strategy
The implementation of and variable definitions for the TFEM were detailed
[Zheng et al. (2010)]. Because the numerical problems in this work are solved based
on the triangular element, the detailed expression of the two-dimensional (2D) tri-
angular element is shown briefly.

2.1. The first stage of interpolation

In the standard FEM, the approximation of interested point x can be obtained by
u(x) = Ni u[i] + Nj u[j] + Nm u[m] , (1)
by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.

Ni = L i , Nj = L j , Nm = L m , (2)
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

in which L is the area coordinate and N is the shape function. Si , Sj , and Sm are
set of element related to nodes i, j, and m, respectively, as shown in [Zheng et al.
(2010)]. The supporting nodes for point x include all the nodes in the element set
Si , Sj , and Sm .
We define qs as the nodal displacement vector of supporting node:
qs = {q1 , q2 , . . . , qns }T , (3)
in which ns is the total number of supporting nodes. For any point in Si , Sj , and
Sm , the classical interpolation can be written in the form,


ns
u(x) = Nl (x)ql . (4)
l=1

The averaged derivatives of node i can be computed as

ns
[i]
ū[i]
,x = N̄l,x ql , (5)
l=1
[i]

[i][e]
N̄l,x = (ωe · Nl,x ), (6)
e∈Si

in which ωe is the weight function of element e and can be computed as follows:



ωe = ∆e ∆e e ∈ Si , (7)
e∈Si

where ∆e is the area of element e.

2.2. The second stage of interpolation

In second stage of interpolation, the averaged derivatives, rather than the nodal
derivatives, are adopted to reproduce the interpolated value at point x,
û(x) = ϕi u[i] + ϕix ū,x
[i] [i]
+ ϕiy ū,y + ϕj u[j] + ϕjx ū[j] [j]
,x + ϕjy ū,y

+ ϕm u[m] + ϕmx ū[m] [m]

,x + ϕmy ū,y , (8)

1250055-3
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

S. C. Wu et al.

in which ϕi , ϕix are ϕiy are constructed to satisfy the relationship of ϕi (xl ) = δil ,
ϕix,x (xl ) = δil . When i is equated to l, δil = 1; otherwise, δil = 0.

ϕi = Li + L2i Lj + L2i Lm − Li L2j − Li L2m , (9)

   
1 1
ϕix = −cj Lm Li + Li Lj Lm + cm Li Lj + Li Lj Lm ,
2 2
(10)
2 2
   
1 1
ϕiy = bj Lm L2i + Li Lj Lm − bm L2i Lj + Li Lj Lm . (11)
2 2

Similarly, ϕj , ϕjx , ϕjy , ϕm , ϕmx , and ϕmy can also be obtained by a cyclic permu-
by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.

[i] [i]
tation of the indices of i, j, and m. Substituting ū,x , ū,y and ϕ into Eq. (8) provides
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

the final trial functions of the twice-interpolation process:

ns
û(x) = N̂l (x)ql , (12)
l=1
[i] [i] [i] [j] [j] [j]
N̂l = ϕi Nl + ϕix N̄l,x + ϕiy N̄l,y + ϕj Nl + ϕjx N̄l,x + ϕjy N̄l,y
[m] [m] [m]
+ ϕm Nl + ϕmx N̄l,x + ϕmy N̄l,y . (13)

2.3. Modification when C 0 node is required

It has been found that the nodal derivatives obtained by the twice-interpolation
process are continuous. However, in solving practical problems, some nodes are
required to be recovered to C 0 , that is, the nodal property of the classical FEM,
such as the nodes located on the essential boundary or material boundary lines.
Therefore, these specified nodes should be modified. Suppose that node i is a C 0
node and that the point of interest x is in element e. Then, the average nodal
derivatives ūi,x can be obtained by
[e]
ūi,x = ui,x . (14)

Note that ūi,y can be computed similar to ūi,x .

It is found that if all the nodes in the problem domain are recovered to be C 0 ,
the TFEM will degenerate to the standard FEM. Thus, by coupling novel TFEM
and standard FEM, a satisfactory ratio of accuracy-versus-cost can be obtained.

3. Stress Intensity Factor Evaluation

3.1. Stress intensity factor
In this section, we briefly discuss the criterion used to specify the direction of the
cracking. Furthermore, we also review the domain form of the interaction integral
used to determine the mixed-mode stress intensity factors (SIF).

1250055-4
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

A Twice-Interpolation Finite Element Method (TFEM ) for Crack Propagation Problems

Under general mixed-mode loadings, the asymptotic near-tip circumferential and

shear stresses take the following form:
   
σθθ KI 1 3 cos(θ/2) + cos(3θ/2)
= √
σrθ 2πr 4 sin(θ/2) + sin(3θ/2)
 
KII 1 −3 sin(θ/2) − 3 sin(3θ/2)
+√ . (15)
2πr 4 cos(θ/2) + 3 cos(3θ/2)

In this paper, we use the maximum circumferential stress criterion, which states
that the crack propagates from the tip in a direction θc so that the circumferential
by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.

stress σθθ , obtains the maximum. The circumferential stress in the direction of the
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

crack propagation is the principal stress. Therefore, the critical angle θc , which
defines the radial direction of propagation, can be determined by setting the shear
stress in Eq. (15) equal to zero. The following equation is thus obtained:
  
1 θ 1 1
√ cos KI sin(θ) + KII (3 cos(θ) − 1) = 0. (16)
2πr 2 2 2
This leads to the equation defining the angle of the crack propagation, θc , in the
tip coordinate system.

KI sin(θc ) + KII (3 cos(θc ) − 1) = 0, (17)

where KI and KII are the pure Mode I and pure Mode II SIF, respectively. Solving
the equation for the angle of crack propagation gives the equation,
  
2
1  KI KI
θc = 2 arctan ± + 8. (18)
4 KII KII

To obtain the SIF, the domain forms of the interaction integrals are used. For
completeness, these are discussed here. The coordinates are the crack tip coordinates
with the x-axis parallel to the crack faces. For general mixed-mode problems, we
have the following relationship between the value of the J integral and the SIFs:

KI2 2
KII
J= + , (19)
E∗ E∗
where E ∗ is defined in terms of Young’s modulus of E and Poisson’s ratio of ν:


E plane stress
E∗ = E . (20)

 1 − ν 2 plane strain

(1) (1)
Two states of a crack body are considered in the computation. State 1 (σij , εij ,
(1) (2) (2) (2)
ui ) corresponds to the actual state, and State 2 (σij , εij , ui ) corresponds to

1250055-5
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

S. C. Wu et al.

an auxiliary state, which will be chosen as the asymptotic fields for Modes I and II.
The J-integral for the sum of the two states is defined as:
J (1+2) = J (1) + J (2) + M (1,2) , (21)
where M (1,2) is called the interaction integral for States 1 and 2
  2 1

(1) ∂ui (2) ∂ui
M (1,2) = W (1,2) δ1j − σij − σij nj dΓ (22)
Γ ∂x1 ∂x1
and W (1,2) is the interaction strain energy
(1) (2) (2) (1)
W (1,2) = σij εij = σij εij . (23)
by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.

After rearranging terms, expressing Eq. (21) for the combined states gives
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

2 (1) (2) (1) (2)

J (1+2) = J (1) + J (2) + ∗ (KI KI + KII KII ). (24)
E
Combining Eqs. (21) and (24) yields
2 (1) (2) (1) (2)
M (1,2) = ∗ (KI KI + KII KII ). (25)
E
Making the judicious choice of State 2 as the pure Mode I asymptotic field with
(2) (2)
KI = 1 and KII = 0, give the Mode I stress intensity factor for State 1 in terms
of the interaction integral
(1) E ∗ (1,Mode I)
KI = M . (26)
2
(2) (2)
Choosing State 2 as Mode II asymptotic fields with KI = 0 and KII = 1 gives
the Mode II stress intensity factor for State 1 in terms of the interaction integral
(1) E ∗ (1,Mode II)
KII = M . (27)
2
The contour integral in Eq. (22) is not flexible for the FEM. Therefore, we must
recast the integral into an equivalent domain form by multiplying the integrand
by a bounded weighting function, (q(x) = 1), on an open set containing the crack
tip and vanishing on an outer prescribed contour C0 . Then, for each contour Γ, as
shown in Fig. 1, the crack faces are assumed to be free of stresses and straight in
the region A, which is bounded by contour C0 . The interaction integral can then
be written as:
  2 1

(1) ∂ui (2) ∂ui
M (1,2)
= W (1,2)
δ1j − σij − σij qmj dΓ, (28)
C ∂x1 ∂x1
where C = Γ + C+ + C− + C0 and m is the outward unit normal to the contour C.
Assuming that the divergence theorem holds true and passes the limit as contour Γ
is shrunk to the crack tip, which is justified by the Dominated Convergence theorem,
the following interaction integral in domain form can be obtained
  2 1

(1) ∂ui (2) ∂ui ∂q
M (1,2) = σij − σij − W (1,2) δ1j dA, (29)
A ∂x1 ∂x1 ∂xj
where mj = −nj on Γ and mj = nj on C0 , C+ and C− .

1250055-6
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

A Twice-Interpolation Finite Element Method (TFEM ) for Crack Propagation Problems

by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.

Fig. 1. Convention at the crack tip. Domain A is enclosed by Γ, C+, C− and C0 .

Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

Fig. 2. Elements selected around the crack tip for the integral.

To numerically evaluate the integral above, domain A is defined in terms of

the length of the elements in contact with the crack tip, with this quantity being
designated as h. Then, domain A is set to be all elements that have a node within
a square. Figure 2 shows a typical set of elements for domain A with rq = 4h.
The weighting function is represented as
q(x, y) = (rq − x)(rq + x)(rq − y)(rq + y). (30)
Note that the coordinates are actually crack tip coordinates with the x1 -axis being
parallel to the crack faces.

3.2. Nodal projection and relaxation

For triangular elements, the mesh irregularity has little effect on the accuracy of the
TFEM [Zheng et al. (2010)]. Therefore, if the node along the crack moves, acceptable
solutions can still be obtained. To model the accurate path of the cracks and modify
the mesh as little as possible, the nodes along the crack path are projected onto the
crack path. In nodal projection, the total number of elements and nodes does not
change, and the nodal coordinates along cracks only require a simple modification.

1250055-7
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

S. C. Wu et al.

Fig. 3. Nodal projection.

by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

As shown in Fig. 3, the nodal projection is divided into three cases. In Case 1, the
crack intersects two edges of the triangle, and the common node of these two edges
is moved onto the crack. In Case 2, the crack intersects only one edge, and thus,
the node near the crack is moved onto the crack. In Case 3, the triangle contains
the crack tip, and the nearest node is moved to the crack tip.
After nodal projection, the crack path is represented by a set of element edges.
The traditional method, known as nodal relaxation, is employed to model the dis-
continuity. In nodal relaxation, each node on the crack path is split into two nodes,
and the elements on two sides of the crack are set to two different nodes, causing
the displacement to be discontinuous.

4. Benchmarks and Experiment

The formulated TFEM is then incorporated into the ALOF and used for solv-
ing some experimental fracture problems, benchmarked in terms of accuracy and
robustness. Numerical results and stress plots are obtained with the ALOF soft-
ware. The variable units used are the international standard unit system, unless
specially denoted.
To examine the accuracy and convergence, the relative L2 error in the displace-
ment norm and energy norm, respectively, are given by

exact − unumerical )2 dΩ
Ω (u 
ed = , (31)
Ω
(uexact )2 dΩ
 
1 exact − εnumerical )T D(εexact − εnumerical )dΩ
2 Ω (ε
ee =  , (32)
1 exact )T Dεexact dΩ
2 Ω (ε

in which the superscript “exact” and “numerical” stand for the exact solution and
a numerical result, such as the present TFEM. For the element used in the TFEM,
a four-point Hammer integration is utilized to evaluate the stiffness matrix. The
nodes located on the essential boundary and material dividing lines are all set to
be C 0 .

1250055-8
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

A Twice-Interpolation Finite Element Method (TFEM ) for Crack Propagation Problems

4.1. Edge crack plate in tension

As the first example, a plate with a crack under the action of uniform tension on
the top edge is analyzed. A sketch of the problem is shown in Zheng et al. [2008],
with the following parameters: H = 1, L = 1 and σ = 1. The material used exhibits
linear elastic behavior with a Young’s modulus of E = 1 × 1010 and ν = 0.25, and
the plane strain conditions are assumed.
The value of the Mode I stress intensity KI factors, as a function of the crack
lengths, are computed as follows:

r = a/H, (33)
by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.

F = 1.120 − 0.231r + 10.550r − 21.710r + 30.382r ,

2 3 4
(34)
√
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

KI = F σ πa. (35)

The mesh used by the TFEM contains 806 uniform nodes. The SIF obtained by
numerical methods is shown in Fig. 4. The computed SIFs agree well with the
benchmark value for all values of the crack length and for similar meshes; the
solutions are better than for the standard FEM. The stress fields in the loading
direction, obtained by both the TFEM and the standard FEM, are plotted in Fig. 5.
It is observed that the stress field obtained by the TFEM is more continuous and
smooth.

4.2. Center-oblique crack in an infinite plate

In this section, a finite plate with an oblique central crack is studied under axial
loading conditions, as shown in Zheng et al. [2008]. The material used is linear
elastic in behavior, with E = 1 × 1010 and ν = 0.25. The plate is subjected to
tensile loading of σ1 and σ2 in the x- and y-directions, respectively. The exact
stress intensity factors are a function of the angle of the crack, β, with the x axial

(a) SIF for difference length of crack. (b) Compare of relative error of SIF.

Fig. 4. SIF of the edge crack problem.

1250055-9
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

S. C. Wu et al.
by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

(a) TFEM solution. (b) FEM solution.

Fig. 5. Stresses in the loading direction.

being given by
√
KI = P sin2 β πa, (36)
√
KII = P sin β cos β πa. (37)
A finite plate with a crack is modeled by a square plate with a width of L = 10, a
crack length of a = 1.414 and an axial loading of P = 1. Since the crack length is
much smaller than the specimen’s dimensions, the numerical results are compared
to the reference solution of a crack in an infinite plate.
The mesh used by the TFEM contained 2392 uniform nodes. To improve the
accuracy, the mesh is fine around the crack tip in the stress plot shown in Fig. 6.
Modes I and II SIFs as a function of the angle β are presented in Fig. 7; the
result agrees well with the reference solution result.

Fig. 6. Stress field of σx and σy .

1250055-10
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

A Twice-Interpolation Finite Element Method (TFEM ) for Crack Propagation Problems

by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

Fig. 7. SIF of the center-oblique oblique problem.

4.3. Plate with two holes

This example concerns the growth of cracks emanating from holes in a square plate,
loaded in tension as shown in Fig. 8. The material used exhibits linear elastic prop-
erties with E = 1 × 1010 and ν = 0.25, and the plane strain conditions are assumed.
The crack path, after 28 steps of crack growth, is plotted in Fig. 8, together with a
reference crack path.

Reference solu.
The TFEM solu.

Path comparison of TFEM and the Reference

Fig. 8. Computational paths using the TFEM with a fine mesh for two holes of the plate.

1250055-11
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

S. C. Wu et al.
by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

Fig. 9. Stress field of σx and σxy after the 28th growth step.

The reference solution is presented for a very fine mesh (12,174 nodes), using
the standard FEM. The numerical solution agrees well with the reference solution,
as shown in Fig. 8. Figure 9 shows the stress distribution of σx and σxy at the 28th
growth step under the ∆a = 0.025.

4.4. Beam under four-point loading

We consider the condition of precracks of a = 0.2 and c = 0.4, with the beam
being subjected to four-point shear loading. The boundary conditions are shown in
Fig. 10, with a beam length of L = 4, a width of b = 1 and a load of P = 1.
The mesh used for the TFEM contains 144 uniform nodes. Together with the
experimental crack path, the computed crack path at step 18 is shown in Fig. 11.
It is clearly observed that using the fracture theory and the TFEM, the crack
propagation path can be simulated with satisfactory accuracy. Additional research
also shows that TFEM approaches the experimental values more closely by increas-
ing the node density and element order.
Figure 12 shows the stress distribution of σx and σy at the 18th growth step
with a cracking increment of ∆a = 0.05.

Fig. 10. Stretch of the beam under four-point loading.

1250055-12
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

A Twice-Interpolation Finite Element Method (TFEM ) for Crack Propagation Problems

Fig. 11. Computational cracking paths at the 18th step and experimental solution.
by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

Fig. 12. Stress distributions of σx and σy after the 18th growth step.

5. Conclusions and Discussions

The TFEM is advantageous compared to the standard FEM, as demonstrated by
the example fracture problems in this paper. Example problems are studied to
examine the accuracy and robustness of the TFEM. The following conclusions can
be drawn:

(a) Using a similar mesh, the TFEM can achieve better accuracy than the tradi-
tional FEM for both Mode I cracks and mixed-mode cracks, in terms of the
accuracy of the SIF of the crack tip.
(b) TFEM can achieve more smoothing stress field of cracking because of its two
stages of interpolation.
(c) The TFEM can simulate crack propagation readily and accurately by incorpo-
rating the TFEM into the ALOF platform.

It must be noted that when compared to the standard FEM, TFEM does not
increase the total degrees of freedom of a model, which can be easily misunderstood
based on the above formulations.
In the second stage of interpolation, each node possesses three interpolation
conditions, including displacements and its derivatives in x- and y-directions. As a
result, linear trial functions can be constructed with only one node, and incomplete

1250055-13
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

S. C. Wu et al.

cubic trial functions can be constructed with two nodes. Thus, using the visible
criterion, discontinuity can be easily embedded in an element. However, discontinu-
ity in an element will cause problems in integration and postprocesses. The nodal
relaxation method can be easily utilized with little modification to the FEM code,
which is always preferred by engineers. Thus, the technique of embedding disconti-
nuity into elements is not investigated in this paper, even though such as technique
is useful for the TFEM.

Acknowledgments
Financial support from the National Nature Science Foundation of China (Grant
by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.

No. 51005068), the Major State Basic Research Development program of China
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

(973 Program, Grant No. 2011CB711105-1), the Independent Research Project of

the State Key Lab of Traction Power of Southwest Jiaotong University (Grant
No. 2012TPL T18) and the Open Research Fund Program of the State Key Lab
of Advanced Design and Manufacturing for Vehicle Body (Grant No. 31115030)
is gratefully acknowledged. Authors also give sincere thanks to Dr. Zheng C and
Dr. Tang XH for their help, constructive discussion and collaboration.

References
Atluri, S. N. and Zhu, T. [1998] “A new meshless local Petrov–Galerkin (MLPG) approach
in computational mechanics,” Comput. Mech. 2(2), 117–127.
Belytschko, T., Lu, Y. Y. and Gu, L. [1995] “Crack propagation by element-free Galerkin
methods,” Eng. Fract. Mech. 51(2), 295–315.
Belytschko, T., Gracie, R. and Ventura, G. [2009] “A review of extended/generalized finite
element methods for material model,” Model. Simul. Mater. Sci. Eng. 17(4), 1–24.
Chen, L., Liu, G. R., Nourbakhsh, N. and Zeng, K. Y. [2010] “A singular edge-based
smoothed finite element method (ES-FEM) for bimaterial interface cracks,” Comput.
Mech. 45(2–3), 109–125.
Liu, G. R., Jiang, Y., Chen, L., Zhang, G. Y. and Zhang, Y. W. [2011] “A singular cell-
based smoothed radial point interpolation method for fracture problems,” Comput.
Struct. 89(13–14), 1378–1396.
Liu, G. R., Nourbakhshnia, N., Chen, L. and Zhang, Y. W. [2010] “A novel general formu-
lation for singular stress field using the ES-FEM method for the analysis of mixed-mode
cracks,” Int. J. Comput. Meth. 7(1), 191–214.
Liu, G. R. and Zhang, G. Y. [2009] “A normed G-space and weakened weak (W2 ) formu-
lation of a cell-based smoothed point interpolation method,” Int. J. Comput. Meth.
6(1), 147–179.
Melenk, J. M. and Babuška, I. [1996] “The partition of unity finite element method: Basic
theory and applications,” Comput. Meth. Appl. Mech. Eng. 139(1–4), 289–314.
Rabczuk, T. and Belytschko, T. [2004] “Cracking particles: A simplified meshfree method
for arbitrary evolving cracks,” Int. J. Numer. Meth. Eng. 61(13), 2316–2343.
Sukumar, N., Moes, N., Moran, B. and Belytschko, T. [2000] “Extended finite element
method for three-dimensional crack modeling,” Int. J. Numer. Meth. Eng. 48(11),
1549–1570.
Swenson, D. and Ingrafea, A. [1996] “Modeling mixed-mode dynamic crack propagation
using finite elements: Theory and applications,” Comput. Mech. 3(6), 381–397.

1250055-14
 2nd Reading
December 28, 2012 17:33 WSPC/0219-8762 196-IJCM 1250055

A Twice-Interpolation Finite Element Method (TFEM ) for Crack Propagation Problems

Tang, X. H., Wu, S. C., Zheng, C. and Zhang, J. H. [2009] “A novel virtual node method
for polygonal elements,” Appl. Math. Mech. Engl. Ed. 30(10), 1233–1246.
Wu, S. C. and Wu, Y. C. [2011] “ALOF: New 3D fatigue crack propagation analysis
software,” Comput. Aided Eng. 20(1), 136–140 (in Chinese).
Zheng, C., Wu, S. C., Tang, X. H. and Zhang, J. H. [2008] “A mesh-free poly-cell Galerkin
(MPG) approach for problems of elasticity and fracture,” CMES 38(2), 149–178.
Zheng, C., Wu, S. C., Tang, X. H. and Zhang, J. H. [2010] “A novel twice-interpolation
finite element method for solid mechanics problems,” Acta Mech. Sin. 26(2), 265–278.
by UNIVERSITY OF NEW ENGLAND LIBRARIES on 01/12/15. For personal use only.
Int. J. Comput. Methods 2012.09. Downloaded from www.worldscientific.com

1250055-15