A novel X-FEM based fast computational method for crack
propagation
Zhenxing Chenga, Hu Wanga*
a

State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body,
Hunan University, Changsha, 410082, P.R. China

Abstract This study suggests a fast computational method for crack propagation, which
is based on the extended finite element method (X-FEM). It is well known that the XFEM might be the most popular numerical method for crack propagation. However,
with the increase of complexity of the given problem, the size of FE model and the
number of iterative steps are increased correspondingly. To improve the efficiency of
X-FEM, an efficient computational method termed decomposed updating reanalysis
(DUR) method is suggested. For most of X-FEM simulation procedures, the change of
each iterative step is small and it will only lead a local change of stiffness matrix.
Therefore, the DUR method is proposed to predict the modified response by only
calculating the changed part of equilibrium equations. Compared with other fast
computational methods, the distinctive characteristic of the proposed method is to
update the modified stiffness matrix with a local updating strategy, which only the
changed part of stiffness matrix needs to be updated. To verify the performance of the
DUR method, several typical numerical examples have been analyzed and the results
demonstrate that this method is a highly efficient method with high accuracy.
Keywords Crack propagation, Fast computational method, Decomposed updating
reanalysis method, Extended finite element method

*

Corresponding author Tel: +86 0731 88655012; fax: +86 0731 88822051.
E-mail address: wanghu@hnu.edu.cn (Hu Wang)
1

1 Introduction
A great amount of engineering practice indicates that the quality and stability of
engineering structures are closely related to the internal crack propagation[1, 2].
Therefore, prediction of the path of crack propagation and analysis of the stability of
crack are significant for estimating the safety and the reliability of engineering
structures. There are many numerical methods have been developed to simulate crack
propagation process, such as finite element method (FEM) [3], boundary element
method (BEM) [4], meshless method [5, 6], edge-based finite element method (ESFEM) [7, 8], numerical manifold method (NMM) [9, 10], extended finite element
method (X-FEM) [11-13] and so on[14]. Compared with above methods, the X-FEM
might be the most popular numerical method for the crack propagation simulation due
to its superiority of modeling both strong and weak discontinuities within a standard
FE framework. The X-FEM was first developed by Belytschko and Black [15]. They
analyzed the crack propagation problem with minimal re-meshing. Then, Dolbow et al
[16]. and Moës et al. [17] improved this method by using a Heaviside function to
enrichment function, and it also has been extended to 3D static crack modeling by
Sukumar et al.[18]. Sequentially, the X-FEM was significantly improved by coupling
with the level set method (LSM) which is used to track both the crack position and tips
[19]. Moreover, the X-FEM has been applied to multiple engineering fields, such as
dynamic crack propagation or branching [20, 21], crack propagation in composites [22]
or shells [23-25], multi-field problems [26], multi-material problems [27, 28],
solidification [29], shear bands [30], dislocations [31] and so on. More details of the
development of X-FEM can be found in Refs [12, 32-34].
Generally, in order to improve the accuracy of simulation, a very refined mesh with a
very small increment of crack propagation or fatigue cycles should be engaged.
Correspondingly, the computational cost is expensive. Therefore, reanalysis algorithm
is used to improve the efficiency.

2

Reanalysis, as a fast computational method, is used to predict the response of modified
structures efficiently without full analysis, and reanalysis method can be divided into
two categories: direct methods (DMs) and approximate methods. DMs can update the
inverse of modified stiffness matrix quickly by Sherman-Morrison-Woodbury lemma
[35, 36] and obtain the exact response of the modified structure, but usually it can only
solve the problems of local or low-rank modifications. In recent decades, many DMs
have been suggested. Such as, Song et al. suggested a novel direct reanalysis algorithm
based on the binary tree characteristic to update the triangular factorization in sparse
matrix solution [37]. Liu et al. applied Cholesky factorization to structural reanalysis
[38]. Huang and Wang suggested an independent coefficient (IC) method for large-scale
problems with local modification [39]. Compared with DMs, approximate methods can
solve the high-rank modifications, but the exact response usually cannot be obtained.
The approximate methods mainly include local approximations (LA), global
approximations (GA) [40], iterative approximations (IA) [41] and combined
approximations (CA) [42-46]. Moreover, many other reanalysis methods have been
proposed in recent years. Zuo et al. combined reanalysis method with genetic algorithm
(GA) [47]. Sun et al. extended the reanalysis method into a structural optimization
process [48]. To improve the efficiency of reanalysis method, He et al. developed a
multiple-GPU based parallel IC reanalysis method [49]. Materna et al. applied the
reanalysis method to nonlinear problems [50].
In this study, a fast computational method is proposed to model crack propagation under
the framework of X-FEM. It is observed that modeling crack propagation by the XFEM will bring the additional DOFs in each iterative step and it will lead to a local
change of stiffness matrix. Considering this characteristic, the decomposed updating
reanalysis (DUR) method theoretically can improve the efficiency of X-FEM
significantly. Moreover, a local updating stiffness matrix strategy is suggested to
improve the efficiency of stiffness matrix assembling. Furthermore, in order to
guarantee the accuracy and efficiency of the DUR method, a local strategy has been
introduced to update the Cholesky factorization of stiffness matrix efficiently.
3

The rest of this paper is organized as follows. The basic theories of X-FEM are briefly
introduced in Section 2. The details of DUR method are described in Section 3. Then,
several numerical examples are tested in Section 4 to investigate the performance of the
DUR method. Finally, conclusions are summarized in Section 5.

2 Basics theories of XFEM
2.1 X-FEM approximation
In the X-FEM, the standard FEM shape function should be enriched by the enrichment
function. Assume that the enriched displacement approximation of X-FEM can be
defined as:

u h (x)   N I (x)u I 
I 

us tan dard

  ( x) N

J  E

J

(x)q J ,

(1)

uenrich

where N I and u I denote the standard FEM shape function and nodal degrees of
freedom (DOF), respectively. The enriched displacement approximation should be
divided into two parts: the standard finite element approximation and partition of unity
enrichment. The （x） means enrichment function while the q J is the additional
nodal degrees of freedom.

Fig. 1 An arbitrary crack line in a structured mesh
Consider an arbitrary crack in a structured mesh as shown in Fig. 1, then Eq.(1) can be
4

rewritten as
4

u (x)   N I (x)u I   H I (x) N I (x)a I     I , (x) N I (x)bI ,
h

I 

(2)

I T  1

I  S

where  is the solution domain,  S is the domain cut by crack,  T is the domain
which crack tip located, H (x) is the shifted Heaviside enrichment and  (x) is the
shifted crack tip enrichment. The details of H (x) and  (x) are given as following:
1
H ( x)  
1

Above crack
,
Below crack




 
{ (x)}4 1  r sin , cos ,sin  sin ,sin  cos  .
2
2
2
 2

(3)

(4)

The discrete X-FEM equations are obtained by substituting Eq.(2) into the principle of
virtual work. Assume that the discrete equations can be defined as
K uu
K T
 ua
 K Tub

K ua
K aa
K Tab

K ub  u  Fu 
   
K ab  a   Fa  ,
b  F 
K bb  
  b

(5)

where K uu is the traditional finite element stiffness matrix, Kua , K aa , K ab are
components with Heaviside enrichment and Kub , K ab , Kbb are components with crack
tip enrichment.
2.2 Crack propagation model
Generally, the direction and magnitude of crack propagation at each iterative step are
used to determine how the crack will propagate. The direction of crack propagation is
found from the maximum circumferential stress criterion and the crack will propagate
in the direction where   is maximum [51]. The angle of crack propagation is
defined as
2


 KI 
1  KI
,
  2 arctan
 sign K II 

8


4  K II
K
 II 



5

(6)

where  is defined in the crack tip coordinate system, K I and K II are the mixedmode stress intensity factors. The details are given in the reference [51].
There are two main quasi-static manners when modeling crack growth. The first one
assumes a constant increment of crack growth at each cycle [16]while the other option
is to assume a constant number of cycles and apply a fatigue crack growth law to predict
the crack growth increment for the fixed number of cycles [52]. In this study, a fixed
increment of crack growth a is considered.

3 Decomposed updating reanalysis method
3.1 Framework of the DUR method
The DUR method is proposed to model quasi-static crack propagation under the
framework of X-FEM and the framework of the DUR method is presented in Fig. 2.

6

Start

Initial equilibrium
equation solution

Separate equilibrium equations

No

Local Cholesky
factorization
updating strategy

Local stiffness
matrix updating
strategy

Yes
Extract balanced
equations

Calculate fundamental
solution system B

Yes

Balanced block

Extract unbalanced
equations

No

No
Calculate the
coefficients vector y

Calculate modified displacements:

Unbalanced block

Reanalysis algorithm

Crack stop
extending?
Yes
End

Fig. 2 The framework of the DUR method
It can be found that the DUR method includes three parts: local updating stiffness
matrix strategy, local Cholesky factorization updating strategy and reanalysis algorithm.
The local stiffness matrix updating strategy is suggested to improve the efficiency of
stiffness matrix assembling according to the characteristic of X-FEM. Moreover,
considering the local change of stiffness matrix, the specific reanalysis algorithm is
used to predict the modified response and improve the efficiency of the X-FEM.
Furthermore, in order to guarantee the accuracy and the efficiency of the suggested
reanalysis method, a local strategy has been introduced to update the Cholesky
factorization of stiffness matrix efficiently. More details can be found in Section 3.2,
3.3 and 3.4.

7

3.2 Local stiffness matrix updating strategy
As mentioned above, the stiffness matrix can be given as the following form, and each
component has been marked by different colors associated with Fig. 1.

Fig. 3 Color-marked stiffness matrix
Consider a small increment of crack based on Fig. 1, and the comparison is shown in
Fig. 4.

Fig. 4 A small increment of crack propagation
It can be found that the stiffness component K uu associated with the standard finite
element approximation will be constant at each iterative step of crack propagation. It
means that the change part of stiffness matrix is the enriched part, which should be
much smaller than the un-enriched part. Furthermore, it can be found that once an
element had become a split element, it will be a split element during entire
computational procedure. It implies that once an element has been enriched by
Heaviside function, the value of stiffness matrix will be a constant in the subsequent
iterative steps. Therefore, the changed part of the stiffness matrix in each iterative step
will be a small part as shown in Fig. 5.
8

Fig. 5 Changed and unchanged part of stiffness matrix
Assume that the unchanged part is a
contains a B  B matrix with

 N  A   N  A matrix and the changed part

 N  A  B ,

B   N  A  coupling matrices. If

N  A  B , only a small part of stiffness matrix needs to be updated in each iterative
steps. Obviously, it should be high efficient and the matrix assembling cost should be
significantly reduced.
3.3 The reanalysis algorithm for X-FEM
The reanalysis algorithm is used to predict the response of the current iterative step by
using the information of the first iterative step. In this study, a reanalysis algorithm is
proposed for X-FEM according to the characteristic of X-FEM. The reanalysis
algorithm avoids the full analysis after the first iterative step, and the response of the
subsequent iterative steps can be efficiently obtained. The details of the reanalysis
algorithm are described as following.
Assume that the equilibrium equation of the i-th iterative step is

K   U   F  ,
i

where U

i 

i

i

(7)

is the displacement in the i-th iterative step, and the equilibrium equation

of the first iterative step can be given as

K   U   F   .
1

1

1

(8)

Assume that the solution of the equation in the i-th iterative step can be defined as

Ui   U1  U  i  1 ,
then substitute Eq.(9) into Eq.(7),

9

(9)





K i  U1  U  Fi   i  1 ,

(10)

K i  U  Fi   K i  U1  i  1 .

(11)

or

Define the residual value of displacement δ as

δ  Fi   K i  U1  i  1 ,

(12)

K   U  δ .

(13)

then Eq.(11) can be written as
i

Consider that only a small part of stiffness matrix will change in every iterative step,
the most part of δ should be zero. Based on this property, the U

i 

can be divided

into two blocks: unbalanced and balanced blocks, according to Eq.(14):





Δ  sum K    K    δ .
i

1

(14)

If Δ  j   0 , the j-th DOF is unbalanced, otherwise the j-th DOF is balanced.
Accord to this, the Eq.(13) can be rewritten as
i
K mm
 i 
 K nm

i
 U m   0 
K mn

   ,
i
K nn   U n  δn 

(15)

where m is the number of balanced DOFs, and n is the number of unbalanced DOFs.
Equation (15) can be rewritten as
i
i
K mm
U m  K mn
U n  0

(16)

i
K nm
U m  K nni  U n  δ n .

(17)

and

Obviously, Eq.(16) is a homogeneous equation set which has infinite solutions, and the
solution U

1

Let U n  0

is one of them.

0 1 0

0 (where 1 is the k-th member of Un , and
T

k  1,2, , n ), then the fundamental solution system of Eq. (16) can be obtained as

10

i
 K i  1K mn

B   mm
,
Enn



(18)

where Enn is a rank-n unit matrix.
Define the general solution of Eq.(16) is

U  By ,

(19)

where y is a dimension-n vector.
Then, substitute Eq.(19) into Eq.(17) to find a unique solution of Eq.(13), obtain

K   K  K 
i
nn

i
nm



K mni  y  δn .

i 1
mm

(20)

Solve Eq.(20), y can be obtained, and then U can be obtained by Eq.(19).
Sequentially, the U

i 

can be obtained by Eq.(9). Because only a small part of stiffness

matrix will change in each iterative step, the number n should not be too large.
i
Therefore, Eq.(20) is a small scale problem. The key point is how to obtain K mm
from

the K   .
1

Assume that

K

1

1
K mm
 1

 K nm

1

K mn
,
1
K nn 

(21)

and then the Cholesky factorization of K   can be defined as
1

0  LTmm

L nn   0

L
K 1  LLT   mm
 L nm
L LT
  mm Tmm
 L nm L mm

LTnm 

LTnn 


L mmLTnm
T
T 
L nm L nm  L nn L nn 

,

(22)

where Lmm and Lnn are lower triangular matrices.
Compared with Eq.(21), it can be found that
i
1
K mm
 K mm
 L mm LTmm .

(23)

Therefore, the fundamental solution system B can be calculated by Eq.(18), and the
Cholesky factorization of stiffness matrix in the first iterative step can directly re-used.
i
In describe this method clearly, the K mm
can be associated with Fig. 4. Assume that

11

the left figure of Fig. 4 shows the crack position in the first iterative step, and the right
figure of Fig. 4 shows the crack position in the second iterative step. Then the stiffness
matrices of X-FEM in the first and the second iterative steps can be written as Fig. 6
and Fig. 7 respectively.

Fig. 6 The stiffness matrix of the first iterative step

Fig. 7 The stiffness matrix of the second iterative step
It can be found that the red-marked part is constant as shown in Fig. 6 and Fig. 7.
 
Compared with Eq.(15), the K mm can be give as
2

 K 1
2
K mm
  1uuT

K ua

K ua1 
,
1
K aa 

(24)

 2 
K ub
,
0 

(25)

0
,
0

(26)

K ab2 
.
 2
K bb


(27)

 
 
 
and the K mn , K nm and K nn can be give as
2

2

2

 2

K mn

K ua2

 0

K  2T
K nm2   ua2
 T
K ub
 2

K nn

 K aa2
 2 T
 
K ab

12

3.4 Local Cholesky factorization updating strategy
Generally, the accuracy and the efficiency of the reanalysis methods are depended on
the percentage of changed part, and if the percentage of changed part is too large, the
accuracy and the efficiency should be unavailable for engineering problems. Therefore,
a suitable critical value of changed percentage should be stipulated, and usually 5% is
chosen for the critical value [53].
In the DUR method, the changed part is composed by unbalanced DOFs. Then the
changed percentage  can be defined as the following form:

  n 100% ,

(28)

N

where n is the number of unbalanced DOFs and N is the number of total DOFs. In order
to guarantee the accuracy and efficiency of the reanalysis method, a suitable critical
value of  should be set. In this study, 5% is also chose as the critical value of  , so
the initial information needs to be updating when the   5% and the key issue is how
to obtain the Cholesky factorization of the initial stiffness matrix efficiently. Therefore,
a local updating Cholesky factorization strategy has been suggested. This strategy is
based on the property that only a small part of stiffness matrix will change in each
iterative step, so that the Cholesky factorization of modified stiffness matrix can be
updating based on the Cholesky factorization of initial stiffness matrix. An example is
given to explain this strategy, which is based on Fig. 6 and Fig. 7.
Assume that Fig. 6 means the initial stiffness matrix and Fig. 7 means the modified
stiffness matrix. Recall the Cholesky factorization of initial stiffness matrix takes the
form

 K 1
K 1   11
1T
K12

1   1
K12
L11 0  L111T

  1

1
K 221  L21
L22
  0

where

13

1 T

L21
1 1T
L L ,
1 T
L22 

(29)

1
L111  chol K11
1
1
L21
 L111 / K12

.

(30)

1
1 1 T
L22
 chol K 221  L21
L 21

1
Compared with Fig. 6, K 11
can be considered equivalent to

1

K11

 K uu1
  1T

K ua

1 
K ua
,
1
K aa 

(31)

1
and K 12
, K 221 can be written as

1

K12

K ub1 
  1  , K 221  K bb1 .

K ab 

(32)
 2

Then the Cholesky factorization of modified stiffness matrix L

L112
L  2
 
L 21
 2

0 
,
2
L22


can be obtained by
(33)

where
 2
1
L112  chol K11
 chol K11
 L111
 2
L212  L112 / K12

,

(34)

K ab2 
.
 2
K bb


(35)

L222  chol K 222  L212L212T

and compared with Fig. 7,
 2

K12

K ua2

 0

 2 
 K aa2
K ub
 2
 , K 22    2T
0 
K ab

1
Then L11
can be reused at subsequent iterative steps of crack propagation, and usually
1
the size of L11
is much large than other parts. Therefore, the strategy should save

much computational cost than calculate the Cholesky factorization of the entire stiffness
matrix directly.

4 Numerical examples
In order to test the accuracy and efficiency of the DUR method, three examples are
tested by the proposed methods. These three cases involve edge and center crack
propagation, concentrated and uniformed load problems, thus the performance of the
14

DUR method could be verified thoroughly. In this study, the comparison has been made
between the DUR and full analysis, and the errors of displacement, Von Mises stress
and Von Mises strain are defined by the following formulas:

Eu 

U PFR  U FA
100%
U FA

(36)

E 

σ PFR  σ FA
100%
σ FA

(37)

where U PFR , σ PFR mean the results of DUR method, and UFA , σ FA mean the results
of full analysis.
Moreover, in order to investigate the performance of the DUR method, the CPU running
time has been recorded and all the simulations were performed on an Intel(R) Core(TM)
i7-5820K 3.30GHz CPU with 32GB of memory within MATLAB R2016b in x64
Windows 7.
4.1 Edge crack in a plate with a hole
The problem shown in Fig. 8 is an adaptation of an example presented in reference [54].
The initial crack length is a0  10mm , the force F  2 104 N , and linear elastic
material behavior is assumed. The material is aluminum 7075-T6 with

E  7.17 104 MPa ,   0.33 and a plane strain state is considered. The increment of
propagation a  1mm . For this problem, the result of experimental test with
specimens 16 mm thick was given by Giner et al. [54] as shown in Fig. 9.

15

65

Ø13

120

6.5

F

Ø20

28.5

17

51

F

13

Fig. 8 The geometry of edge crack in a plate with a hole

Fig. 9 The result of edge crack in a plate with a hole in the reference [54]
Then solve this case by the DUR method and only the first iterative step need to be
calculated by full analysis method while other iterative steps should be predicted
efficiently by the DUR method. In order to investigate the accuracy of the DUR method,
The comparisons of displacement and stress between the DUR and full analysis are
illustrated in Fig. 10 and Fig. 11, respectively. It is obviously that the results of DUR
method and full analysis are almost the same. Compared with Fig. 9, it is proved that
the DUR method is accurate. Moreover, the errors of each iterative step which is defined
by Eq.(36) and Eq.(37) are shown in Fig. 12. It can be found that the DUR method is
accurate. The computational costs of the DUR and full analysis are also listed in Tab.
16

1. It shows that the computational cost of the DUR method is much cheaper than the
full analysis method. Moreover, the computational results of crack tip coordinates are
shown in Tab. 2 and the results of the DUR and full analysis are also the same.

Fig. 10 Displacement contour of edge crack in a plate with a hole

Fig. 11 Von Mises stress contour of edge crack in a plate with a hole

17

Fig. 12 The errors of each iterative step
Tab. 1 Performance comparison of edge crack in a plate with a hole
Computational Time/s

Average Errors

DUR

Full analysis

Displacement

Von Mises Stress

0.968

18.953

4.8916e-13

1.7356e-12

Tab. 2 Crack tip coordinates of edge crack in a plate with a hole
Iterative steps

0
1
5
10
15
20
25
30

Crack tip coordinates of DUR
Method

Crack tip coordinates of full analysis

X-coordinate

Y-coordinate

X-coordinate

Y-coordinate

0.000
10.000
10.998
14.994
19.993
24.986
29.918
34.605
37.695

67.500
67.500
67.559
67.727
67.805
67.560
66.766
65.059
61.523

0.000
10.000
10.998
14.994
19.993
24.986
29.918
34.605
37.695

67.500
67.500
67.559
67.727
67.805
67.560
66.766
65.059
61.523

4.2 Edge crack in a plate with a circular inclusion
As shown in Fig. 13, a plate with a circular inclusion is considered. Assume that there
is an edge crack in the left edge of the plate. The aluminum 7075-T6 is used as the
4
material of plate in the previous case, the Young’s modulus E1  7.17 10 MPa , the

Poisson’s ratio   0.33 , and a plane strain state is assuming. The material of inclusion
18

is considered as carbon fiber reinforced composite and assume that the Young’s
modulus E2  2.1105 MPa , the Poisson’s is the same as the plate. The initial crack
length a0  10mm , and the uniformed load q  50 N / mm . The increment of
propagation a  1mm . Then the DUR method is used to solve this problem and only
the first iterative step needs to be calculated by full analysis method while others
iterative steps should be predicted efficiently by the DUR method.
q

120

60

30

30

60

R15

Inclusion

q

Fig. 13 The geometry of edge crack in a plate with a circular inclusion
The comparisons of displacement and stress are shown as Fig. 14 and Fig. 15, Fig. 15
presents the Von Mises stress results of iterative step 1, 20 and 40. It is obviously that
the result of DUR method is very close to the full analysis method. Furthermore, an
error list of some selected iterative steps are shown as Fig. 16. It can be found that the
DUR method is highly accurate because the maximum of error is 1.0 1010 . Moreover,
the computational costs of the DUR and full analysis are also listed in Tab. 3. It shows
that the efficiency of the DUR method is much higher than the full analysis method.
Moreover, the crack tip coordinates are also listed in Tab. 4 and the results of the DUR
and full analysis are matched.

19

Fig. 14 Displacement contour of edge crack in a plate with a circular inclusion

Fig. 15 Von Mises stress contour of edge crack in a plate with a circular inclusion
20

Fig. 16 The errors of each iterative step
Tab. 3 Performance comparison of edge crack in a plate with a circular inclusion
Computational Time/s

Average Errors

DUR

Full analysis

Displacement

Von Mises Stress

4.078

17.203

1.0534e-11

1.0951e-11

Tab. 4 Crack tip coordinates of edge crack in a plate with a circular inclusion
Iterative steps

0
1
5
10
15
20
25
30
35
40
45
49

Crack tip coordinates of DUR
Method

Crack tip coordinates of full analysis

X-coordinate

Y-coordinate

X-coordinate

Y-coordinate

0.000
10.000
10.987
14.941
19.828
24.647
29.401
34.135
38.956
43.911
48.908
53.908
55.908

60.000
60.000
59.842
59.240
58.188
56.857
55.308
53.698
52.383
51.754
51.594
51.552
51.545

0.000
10.000
10.987
14.941
19.828
24.647
29.401
34.135
38.956
43.911
48.908
53.908
55.908

60.000
60.000
59.842
59.240
58.188
56.857
55.308
53.698
52.383
51.754
51.594
51.552
51.545

21

4.3 Center crack in a plate with a circular inclusion and a hole
A center crack in a plate with a circular inclusion and a hole is considered as shown in
Fig. 17. The material of plate is aluminum 7075-T6, the Young’s modulus

E1  7.17 104 MPa , the Poisson’s ratio   0.33 , and a plane strain state is assuming.
The material of inclusion is considered as carbon fiber reinforced composite and
assume that the Young’s modulus E2  2.1105 MPa , the Poisson’s is the same as the
plate. The initial crack length a0  10mm , and the uniformed load q  50 N / mm . The
increment of propagation a  1mm . Then, the DUR method is used to solve this
problem and only the first iterative step needs to be calculated by the full analysis
method while other iterative steps should be predicted efficiently by the DUR method.
q

30

60

R10

120

Inclusion

30

60

R10

30

q

Fig. 17 The geometry of center crack in a plate with a circular inclusion and a hole
The displacement and Von Mises stress contours of iterative step 1, 10 and 24 are shown
as Fig. 18 and Fig. 19. It is obviously that the result of the DUR method is very close
to the full analysis method. Furthermore, an error list of some selected iterative steps
are shown as Fig. 20. It can be found that the DUR method is accurate because the
maximum of error is 4.0 1012 . Moreover, the computational costs of the DUR and
full analysis are also listed in Tab. 5. It shows that the efficiency of the DUR method is
22

much higher than full analysis method. Besides this, the crack tip coordinates are also
listed in Tab. 6 and the result of DUR and full analysis are matched.

Fig. 18 Displacement contour of center crack in a plate with a circular inclusion and a
hole

23

Fig. 19 Von Mises stress contour of center crack in a plate with a circular inclusion
and a hole

Fig. 20 The errors of each iterative step
Tab. 5 Performance comparison of center crack in a plate with a circular inclusion and
a hole
Computational Time/s

Average Errors

DUR

Full analysis

Displacement

Von Mises Stress

2.3438

8.4375

4.9236e-13

5.6032e-13

Tab. 6 Crack tip coordinates of center crack in a plate with a circular inclusion and a
hole
Iterative steps

Crack tip coordinates of DUR
Method
Crack tip 1

Crack tip 2
24

Crack tip coordinates of full analysis
Crack tip 1

Crack tip 2

0
1
5
10
15
20
24

(25.000, 60.000)
(24.012, 60.152)
(20.058, 60.511)
(15.107, 60.748)
(10.133, 61.037)
(5.150,61.077)
(1.163, 61.116)

(35.000, 60.000)
(35.988, 60.152)
(39.942, 60.510)
(44.894, 60.746)
(49.864, 61.032)
(54.858, 61.037)
(58.857, 61.019)

(25.000, 60.000)
(24.012, 60.152)
(20.058, 60.511)
(15.107, 60.748)
(10.133, 61.037)
(5.150,61.077)
(1.163, 61.116)

(35.000, 60.000)
(35.988, 60.152)
(39.942, 60.510)
(44.894, 60.746)
(49.864, 61.032)
(54.858, 61.037)
(58.857, 61.019)

4.4 Accuracy and efficiency comparison
Three numerical examples have been tested in this section and it can be found that the
DUR is an accurate and efficient and method. Moreover, a representative case has been
calculated by the DUR method under different computational scales from 1000 to
100,000 to fully investigate the efficiency of the DUR method. The log-log plots of
comparison results are shown in Fig. 21, Fig. 22 and the error analysis is also shown.

Fig. 21 The comparison of computational cost used in solving equilibrium equations

25

Fig. 22 The comparison of computational cost used in stiffness matrix updating
It can be found that the efficiency of the DUR method is much higher than the full
analysis method, and the advantage is more obvious for large scale problems. It is
obvious that the accuracy should be improved significantly with the increase of DOFs.
Accord to Fig. 21, it can be observed that the accuracy of the DUR method is very high,
and the DUR can be regarded as an exact method. Moreover, the computational cost of
stiffness matrix updating was also plotted, and Fig. 22 shows that the local updating
strategy largely reduces the computational cost of stiffness matrix updating than
traditional global updating strategy.

5 Conclusions
In this study, the DUR method is proposed for crack propagation. The DUR method
consists of three strategies: local stiffness matrix updating, decomposed updating
reanalysis, local Cholesky factorization updating strategies. Considering the
characteristic of local change of stiffness matrix during X-FEM iterative procedure, the
local stiffness matrix updating method can achieve the modified stiffness matrix quickly,
and the local Cholesky factorization updating strategy is used to guarantee the accuracy
and efficiency of the DUR method. More importantly, the decomposed updating
26

reanalysis strategy is suggested to improve the efficiency of solving the equilibrium
equations significantly. Therefore, the DUR method not only reduces the computational
cost of solving equilibrium equations but also saves computational cost of stiffness
matrix assembling.
Numerical examples show that the DUR method is accurate for the crack propagation.
For both the edge and center crack propagation, the accuracy of DUR method is very
high. The log-log plots show that the efficiency of DUR method is much more higher
than full analysis, and the advantage is more obvious for large scale problems.
Moreover, compared with other reanalysis method, the comparisons of stress between
DUR and full analysis are made, and the stress can be obtained accurately and
efficiently.
Acknowledgements This work has been supported by Project of the Key Program of
National Natural Science Foundation of China under the Grant Numbers 11572120,
National Key Research and Development Program of China 2017YFB0203701.

References
[1] Kanaun S. Discrete model of hydraulic fracture crack propagation in homogeneous isotropic
elastic media. International Journal of Engineering Science. 2017;110:1-14.
[2] Areias P, Msekh MA, Rabczuk T. Damage and fracture algorithm using the screened Poisson
equation and local remeshing. Engineering Fracture Mechanics. 2016;158:116-43.
[3] Branco R, Antunes FV, Costa JD. A review on 3D-FE adaptive remeshing techniques for crack
growth modelling. Engineering Fracture Mechanics. 2015;141:170-95.
[4] Santana E, Portela A. Dual boundary element analysis of fatigue crack growth, interaction and
linkup. Engineering Analysis with Boundary Elements. 2016;64:176-95.
[5] Tanaka S, Suzuki H, Sadamoto S, Imachi M, Bui TQ. Analysis of cracked shear deformable
plates by an effective meshfree plate formulation. Engineering Fracture Mechanics. 2015;144:14257.
[6] Khosravifard A, Hematiyan MR, Bui TQ, Do TV. Accurate and efficient analysis of stationary
and propagating crack problems by meshless methods. Theoretical and Applied Fracture Mechanics.
2017;87:21-34.
[7] Nguyen-Xuan H, Liu GR, Bordas S, Natarajan S, Rabczuk T. An adaptive singular ES-FEM
for mechanics problems with singular field of arbitrary order. Computer Methods in Applied
Mechanics and Engineering. 2013;253:252-73.
[8] Liu GR, Nourbakhshnia N, Zhang YW. A novel singular ES-FEM method for simulating
27

singular stress fields near the crack tips for linear fracture problems. Engineering Fracture
Mechanics. 2011;78(6):863-76.
[9] Liu Z, Zheng H, Sun C. A domain decomposition based method for two-dimensional linear
elastic fractures. Engineering Analysis with Boundary Elements. 2016;66:34-48.
[10] Zhang HH, Li LX, An XM, Ma GW. Numerical analysis of 2-D crack propagation problems
using the numerical manifold method. Engineering Analysis with Boundary Elements.
2010;34(1):41-50.
[11] Zeng Q, Liu Z, Xu D, Wang H, Zhuang Z. Modeling arbitrary crack propagation in coupled
shell/solid structures with X-FEM. International Journal for Numerical Methods in Engineering.
2016;106(12):1018-40.
[12] Belytschko T, Gracie R, Ventura G. A review of extended/generalized finite element methods
for material modeling. Modelling and Simulation in Materials Science and Engineering.
2009;17(4):043001.
[13] Chopp DL, Sukumar N. Fatigue crack propagation of multiple coplanar cracks with the
coupled extended finite element/fast marching method. International Journal of Engineering
Science. 2003;41(8):845-69.
[14] Cuomo M, Contrafatto L, Greco L. A variational model based on isogeometric interpolation
for the analysis of cracked bodies. International Journal of Engineering Science. 2014;80:173-88.
[15] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing.
International journal for numerical methods in engineering. 1999;45(5):601-20.
[16] Dolbow J, Belytschko T. A finite element method for crack growth without remeshing.
International journal for numerical methods in engineering. 1999;46(1):131-50.
[17] Belytschko T, Moës N, Usui S, Parimi C. Arbitrary discontinuities in finite elements.
International Journal for Numerical Methods in Engineering. 2001;50(4):993-1013.
[18] Sukumar N, Moës N, Moran B, Belytschko T. Extended finite element method for threedimensional crack modelling. International Journal for Numerical Methods in Engineering.
2000;48(11):1549-70.
[19] Stolarska M, Chopp D, Moës N, Belytschko T. Modelling crack growth by level sets in the
extended finite element method. International journal for numerical methods in Engineering.
2001;51(8):943-60.
[20] Xu D, Liu Z, Liu X, Zeng Q, Zhuang Z. Modeling of dynamic crack branching by enhanced
extended finite element method. Computational Mechanics. 2014;54(2):489-502.
[21] Song J-H, Areias PMA, Belytschko T. A method for dynamic crack and shear band
propagation with phantom nodes. International Journal for Numerical Methods in Engineering.
2006;67(6):868-93.
[22] Huynh D, Belytschko T. The extended finite element method for fracture in composite
materials. International Journal for Numerical Methods in Engineering. 2009;77(2):214-39.
[23] Ahmed A, van der Meer FP, Sluys LJ. A geometrically nonlinear discontinuous solid-like
shell element (DSLS) for thin shell structures. Computer Methods in Applied Mechanics and
Engineering. 2012;201-204:191-207.
[24] Zhuang Z, Cheng B. A novel enriched CB shell element method for simulating arbitrary crack
growth in pipes. Science China Physics, Mechanics and Astronomy. 2011;54(8):1520-31.
[25] Areias P, Belytschko T. Non‐linear analysis of shells with arbitrary evolving cracks using
XFEM. International Journal for Numerical Methods in Engineering. 2005;62(3):384-415.
28

[26] Zilian A, Legay A. The enriched space–time finite element method (EST) for simultaneous
solution of fluid–structure interaction. International Journal for Numerical Methods in Biomedical
Engineering. 2008;75(3).
[27] Sukumar N, Chopp DL, Moës N, Belytschko T. Modeling holes and inclusions by level sets
in the extended finite-element method. Computer methods in applied mechanics and engineering.
2001;190(46):6183-200.
[28] Zhuang Z, Cheng BB. Equilibrium state of mode-I sub-interfacial crack growth in bimaterials. International Journal of Fracture. 2011;170(1):27-36.
[29] Chessa J, Smolinski P, Belytschko T. The extended finite element method (XFEM) for
solidification problems. International Journal for Numerical Methods in Engineering.
2002;53(8):1959-77.
[30] Areias P, Belytschko T. Two‐scale shear band evolution by local partition of unity.
International Journal for Numerical Methods in Engineering. 2006;66(5):878-910.
[31] Belytschko T, Gracie R. On XFEM applications to dislocations and interfaces. International
Journal of Plasticity. 2007;23(10):1721-38.
[32] Yazid A, Abdelkader N, Abdelmadjid H. A state-of-the-art review of the X-FEM for
computational fracture mechanics. Applied Mathematical Modelling. 2009;33(12):4269-82.
[33] Abdelaziz Y, Hamouine A. A survey of the extended finite element. Computers & Structures.
2008;86(11-12):1141-51.
[34] Fries TP, Belytschko T. The extended/generalized finite element method: an overview of the
method and its applications. International Journal for Numerical Methods in Engineering.
2010;84(3):253-304.
[35] Woodbury MA. Inverting modified matrices. Memorandum report. 1950;42(106):336.
[36] Sherman J, Morrison WJ. Adjustment of an Inverse Matrix Corresponding to a Change in
One Element of a Given Matrix. The Annals of Mathematical Statistics. 1950;21(1):124-7.
[37] Song Q, Chen P, Sun S. An exact reanalysis algorithm for local non-topological high-rank
structural modifications in finite element analysis. Computers & Structures. 2014;143:60-72.
[38] Liu HF, Wu BS, Li ZG. Method of Updating the Cholesky Factorization for Structural
Reanalysis with Added Degrees of Freedom. Journal of Engineering Mechanics. 2014;140(2):38492.
[39] Huang G, Wang H, Li G. A reanalysis method for local modification and the application in
large-scale problems. Structural and Multidisciplinary Optimization. 2013;49(6):915-30.
[40] Haftka RT, Nachlas JA, Watson LT, Rizzo T, Desai R. Two-point constraint approximation in
structural optimization. Computer Methods in Applied Mechanics and Engineering.
1987;60(3):289-301.
[41] Li Z, Wu B. A preconditioned conjugate gradient approach to structural reanalysis for general
layout modifications. International journal for numerical methods in engineering. 2007;70(5):50522.
[42] Kirsch U. A unified reanalysis approach for structural analysis, design, and optimization.
Structural and Multidisciplinary Optimization. 2003;25(2):67-85.
[43] Chen SH, Yang ZJ. A universal method for structural static reanalysis of topological
modifications. International Journal for Numerical Methods in Engineering. 2004;61(5):673-86.
[44] Wu BS, Lim CW, Li ZG. A finite element algorithm for reanalysis of structures with added
degrees of freedom. Finite Elements in Analysis and Design. 2004;40(13-14):1791-801.
29

[45] Zheng SP, Wu BS, Li ZG. Vibration reanalysis based on block combined approximations
with shifting. Computers & Structures. 2015;149:72-80.
[46] Zuo W, Bai J, Yu J. Sensitivity reanalysis of static displacement using Taylor series expansion
and combined approximate method. Structural and Multidisciplinary Optimization. 2015:1-7.
[47] Zuo W, Xu T, Zhang H, Xu T. Fast structural optimization with frequency constraints by
genetic algorithm using adaptive eigenvalue reanalysis methods. Structural and Multidisciplinary
Optimization. 2011;43(6):799-810.
[48] Sun R, Liu D, Xu T, Zhang H, Zuo W. New Adaptive Technique of Kirsch Method for
Structural Reanalysis. AIAA Journal. 2014;52(3):486-95.
[49] He G, Wang H, Li E, Huang G, Li G. A multiple-GPU based parallel independent coefficient
reanalysis method and applications for vehicle design. Advances in Engineering Software.
2015;85:108-24.
[50] Materna D, Kalpakides VK. Nonlinear reanalysis for structural modifications based on
residual increment approximations. Computational Mechanics. 2016;57(1):1-18.
[51] Erdogan F, Sih G. On the crack extension in plates under plane loading and transverse shear.
Journal of basic engineering. 1963;85(4):519-27.
[52] Gravouil A, Moës N, Belytschko T. Non‐planar 3D crack growth by the extended finite
element and level sets—Part II: Level set update. International Journal for Numerical Methods in
Engineering. 2002;53(11):2569-86.
[53] Kirsch U. Design-Oriented Analysis of Structures. Dordrecht: Kluwer Academic Publishers;
2002.
[54] Giner E, Sukumar N, Tarancón JE, Fuenmayor FJ. An Abaqus implementation of the
extended finite element method. Engineering Fracture Mechanics. 2009;76(3):347-68.

30

