INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING

Int. J. Numer. Meth. Engng 2007; 71:14661482

Published online 12 February 2007 in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/nme.2001
A multiscale projection method for macro/microcrack simulations
Stefan Loehnert
,
and Ted Belytschko
Department of Mechanical Engineering, Northwestern University, 2145 Sheridan Road,
Evanston, IL 60208, U.S.A.
SUMMARY
We present a new multiscale method for crack simulations. This approach is based on a two-scale
decomposition of the displacements and a projection to the coarse scale by using coarse scale test
functions. The extended nite element method (XFEM) is used to take into account macrocracks as well
as microcracks accurately. The transition of the eld variables between the different scales and the role
of the microeld in the coarse scale formulation are emphasized. The method is designed so that the ne
scale computation can be done independently of the coarse scale computation, which is very efcient and
ideal for parallelization. Several examples involving microcracks and macrocracks are given. It is shown
that the effect of crack shielding and amplication for crack growth analyses can be captured efciently.
Copyright q 2007 John Wiley & Sons, Ltd.
Received 8 December 2006; Revised 20 December 2006; Accepted 22 December 2006
KEY WORDS: multiscale; cracks; XFEM
1. INTRODUCTION
Multiscale methods offer great promise in modelling the effects of microstructural features such as
microcracks. Different concepts of multiscale methods are the variational multiscale methods [1],
the homogenized Dirichlet projection method [2, 3], the FE
2
method [4], domain decomposition
methods [57], multiscale methods based on homogenization techniques [8], and other concurrent
methods, e.g. [9]. When local effects in certain subdomains of a large structure need to be modelled,
multiscale methods are attractive since standard homogenization techniques, which in general are
very useful for the bulk of the structure, usually fail to predict localizing phenomena accurately.
Using a multiscale method in those domains overcomes this difculty since it is possible to

Correspondence to: Stefan Loehnert, Department of Mechanical Engineering, Northwestern University, 2145 Sheridan
Road, Evanston, IL 60208, U.S.A.

E-mail: s-loehnert@northwestern.edu
Contract/grant sponsor: Army Research Ofce; contract/grant number: W911NF-05-1-0049
Copyright q 2007 John Wiley & Sons, Ltd.
A MULTISCALE PROJECTION METHOD FOR MACRO/MICROCRACK SIMULATIONS 1467
predict the localizing phenomena with a detailed microstructural computation. A method using
superimposed meshes for the improvement of solutions in local areas applied to linear elastostatics
can be found in [10].
Two of the most common localization phenomena which can be observed in a variety of materials
are cracks and shear bands. For the investigation of crack propagation, stress or strain criteria in the
vicinity of the existing crack tip or the stress intensity factors of the crack tip are used to predict
the crack propagation and its direction. In many ceramic materials it can be observed that around
the crack tip of an existing macrocrack, microcracks nucleate and develop. These microcracks
have a signicant inuence on the stress eld in the vicinity of the macrocrack tip. They can
result in crack shielding or crack amplication and thus greatly inuence the propagation of the
macrocrack. A brute force, single-scale analysis of such problems would require an extreme local
adaptive mesh renement in the vicinity of each macrocrack tip. Rening to the needed scale
around such a local area would lead to poorly conditioned equations and high computational cost
which make such an approach awkward.
This paper presents a computationally efcient method for treating subscale (i.e. microscale)
behaviour for cracked bodies. The approach is based on exploiting the separation of scales typical
of microscale/macroscale problems. As shown here, such a scale separation enables the problem to
be treated in two parts: a macroscale model which embodies the overall effects of the microcracks,
and a microscale model which models detailed interaction of cracks. The novel contribution of
this paper is that we show that the acceptable accuracy can be obtained by using the coarse mesh
test functions in conjunction with the effects of the microcracks.
Our method is formulated within the context of the extended nite element method (XFEM)
[11, 12] which makes it possible to perform crack simulations with an arbitrary number of randomly
distributed cracks without remeshing whenever a crack propagates. XFEM has been applied to large
arrays of cracks at a single scale by Budyn et al. [13] and in microcrack/macrocrack strategies by
Guidault et al. [14]. The synthesis of XFEM with level sets we use here was proposed by Stolarska
et al. [15] and Belytschko et al. [16]. Improvements of XFEM for crack applications have been
proposed by Xiao and Karihaloo [17] and Laborde et al. [18]. Extensions to three dimensions are
reported by Gravouil et al. [19, 20] and Areias and Belytschko [21]. Applications to shear bands
are given by Areias and Belytschko [22] and Samaniego and Belytschko [23]. Applications to
dislocations are given in Ventura et al. [24] and Gracie et al. [25]. Thus, the method proposed
here is also applicable to localized shear bands and dislocations.
In Section 2, the multiscale method for crack simulations is derived. In particular, this involves the
transition of the eld variables between the different scales and consequently the decoupling of the
computations for different scales as well as the separation of the structural details for each scale. In
Section 3, the discretized equations and the proposed projection of the eld variables are presented.
Section 4 shows the solution procedure. In Section 5, the multiscale method is applied to several
tests for which analytical results can be found in literature and to some crack shielding/amplication
problems involving a larger number of randomly distributed microcracks around the crack tips of
a macrocrack. Section 6 gives some conclusions and an outlook to future work.
2. THE MULTISCALE STRATEGY
For simplicity, here we restrict ourselves to two scales even though the multiscale strategy described
in the following can be extended to an arbitrary number of scales. We assume that the behaviour
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
1468 S. LOEHNERT AND T. BELYTSCHKO
Figure 1. Fine scale (microscale) and coarse scale (macroscale).
of the problem, either because of its geometry or the evolution of the cracks, exhibits two scales
of behaviour. We denote by l
0
the scale of coarse scale features and by l
1
the scale of ne
scale features. Furthermore, we assume that the geometry and response are characterized by scale
separation, i.e. that the ne scale features are signicantly smaller than the coarse scale features,
i.e. l
0
l
1
.
Let the domain of the structure under consideration be
0
. We denote the boundary of
0
by *
0
and subdivide it into complementary boundaries *
0
t
and *
0
u
where the tractions and
displacements, respectively, are prescribed. The domain
0
consists of parts where a ne scale
analysis is not necessary and other parts
1

0
where the ne scale behaviour has great inuence
on the coarse scale behaviour (see Figure 1). It is assumed that along the boundaries *
1
, the
uctuations of the eld variables due to the microstructure are negligible. This strongly affects
the choice of the area within which the detailed microstructure has to be resolved accurately. On
the coarsest scale, only cracks which are longer than a typical nite element size are explicitly
considered. In the ne scale model, the microcracks and all the discontinuities on coarser scales
are taken into account.
We consider a body subjected to body forces f so that the equilibrium equation is
div(r) + f =0 (1)
In addition, the traction boundary conditions on *
0
t
and the crack surfaces must be satised.
r n =t on *
0
t
(2)
r n =t
c
on

i

i
D
(3)
where n is the unit normal to the surface pointing out from the body. The coarse scale model
includes macrocracks which correspond to discontinuities

0
D
in the displacement eld, while
the microcracks in
1
corresponds to discontinuities

1
D
. For convenience, we let
i
D
=

i
D
,
i =0, 1.
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
A MULTISCALE PROJECTION METHOD FOR MACRO/MICROCRACK SIMULATIONS 1469
The macro displacement (coarse scale displacement) is given by
u
0
(X) =u
0
C
(X) + u
0
D
(X) (4)
where u
0
C
(X) is the continuous part of the macro displacement and u
0
D
(X) is the discontinuous
part of the macro displacement. We require that u
0
C
U
0
C
and u
0
D
U
0
D
where
U
0
C
={u|u H
1
(
0
), u= u on *
0
u
} (5)
U
0
D
={u|u H
1
(
0
\
0
D
), u=0 on *
0
u
} (6)
In
1
, the displacement eld is given by
u
1
=u
0
+ u
1
(7)
where the micro displacement eld u
1
is given by
u
1
(X) = u
1
C
(X) + u
1
D
(X) (8)
where u
1
C
and u
1
D
are, respectively, continuous and discontinuous parts of the micro displacement
eld, u
1
C
U
1
C
and u
1
D
U
1
D
and
U
1
C
={u|u H
1
(
0
), u=0 on *
1
, u=0 in
0
\
1
} (9)
U
1
D
={u|u H
1
(
0
\{
0
D

1
D
}), u=0 on *
1
, u=0 in
0
\
1
} (10)
We include both cohesive cracks and cracks with singular crack tips in the formulation, so that
crack bridging can be treated, cf. Bao and Suo [26].
The weak form on
0
will be based on the macro displacement eld, i.e. it will resolve
equilibrium on the coarse level. For this purpose, we employ as test functions the functions from
the macro displacement space U
0
C
U
0
D
. The weak form is obtained by multiplying the equilibrium
equation by the test function g
0
(X)
_

0
div(r(u
0
+ u
1
)) g
0
d=
_

0
f g
0
d +
_
*
0
t
t g
0
d* +
_

0
D
t
c
'g
0
D
( d (11)
where
g
0
(X) =g
0
C
(X) + g
0
D
(X) (12)
and g
0
C


U
0
C
, g
0
D


U
0
D
with

U
0
C
={g|g H
1
(
0
), g =0 on *
0
u
} (13)

U
0
D
={g|g H
1
(
0
\
0
D
), g =0 on *
0
u
} (14)
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
1470 S. LOEHNERT AND T. BELYTSCHKO
Integration by parts yields
_

0
r(u
0
+ u
1
) : grad
sym
(g
0
) d
=
_

0
f g
0
d +
_
*
0
t g
0
d* +
_

0
D
t
c
'g
0
D
( d g
0
C


U
0
C
, g
0
D


U
0
D
(15)
where t
c
is the cohesive traction across the cracks.
Note that the microcrack cohesive tractions are eliminated from the equilibrium equation because
the macroscale test eld is continuous across the microcracks, i.e. the discontinuities
1
D
. Thus, the
use of the coarse scale test functions eliminates the role of the cohesive tractions on the microcracks
from the coarse scale equilibrium equations. However, their effect is still included in the sense
that u
1
depends on the microcracks. When the cracks are traction free, they do not appear in the
coarse scale work expression, but it is easy to verify that the EulerLagrange equations (15) are
_

0
r(u
0
+ u
1
) : grad
sym
(g
0
) d
=
_

0
f g
0
d +
_
*
0
t g
0
d* g
0
C


U
0
C
, g
0
D


U
0
D
(16)
with r n=0 on
0
D
.
For the purpose of obtaining the ne scale equilibrium equations, we dene the micro test
functions
g
1
(X) =g
1
C
(X) + g
1
D
(X) (17)
where g
1
C
U
1
C
and g
1
D
U
1
D
. We then multiply the equilibrium equation by these test functions,
which yields, after integrating by parts, that
_

1
r(u
0
+ u
1
) : grad
sym
(g
1
) d
=
_

1
f g
1
d +
_

0
D
t
c
'g
1
D
( d +
_

1
D
t
c
'g
1
D
( d (18)
Note that cohesive forces of both the coarse scale and the ne scale cohesive cracks enter the ne
scale equations. As for the macroscale, when the cohesive forces vanish, the last terms do not
appear. In this case, the EulerLagrange equations are
_

1
r(u
0
+ u
1
) : grad
sym
(g
1
) d=
_

1
f g
1
d (19)
with natural boundary conditions r n=0 on
0
D

1
D
, i.e. the tractions vanish on the macrocrack
and microcrack surfaces. In the above, we did not account for crack closure, which would entail
either changing (15) and (18) to variational inequalities or checking for the absence of crack
closure after the computation. We did the latter.
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
A MULTISCALE PROJECTION METHOD FOR MACRO/MICROCRACK SIMULATIONS 1471
Figure 2. Coarse scale mesh and ne scale mesh.
3. DISCRETIZATION
The domain
0
is subdivided into elements
0
e
. The subdomain
1
is subdivided into elements
that are congruent to the elements
0
e
, so that for any element
0
e
on the coarse scale there exists a
set of elements S
e
so that
0
e
=

S
e

(see Figure 2). The ne scale and coarse scale domains

will be treated separately, with the ne scale solution only providing the ne scale stress eld
to the coarse scale equations, see Equation (15). Instead of discretizing the eld u
1
, we directly
discretize u
1
. Following Mo es et al. [12] the displacement approximations can then be written in
the form
u
0
h
=
n
0
n

I =1
N
0
I
_
u
0
I
+
n
enr

j =1
f
j
a
0
j I
_
(20)
u
1
h
=
n
1
n

I =1
N
1
I
_
u
1
I
+
n
enr

j =1
f
j
a
1
j I
_
(21)
where f
j
are the enrichment functions chosen to account for crack behaviour, N
0
I
the shape
functions for the coarse scale, N
1
I
the shape functions for the ne scale, u
0
I
and u
1
I
the standard
nodal displacement degrees of freedom for the coarse and ne scale, respectively, and a
0
j I
and a
1
j I
are the coarse scale and ne scale enrichment parameters.
The nodes belonging to elements within which a crack ends are enriched with the crack tip
enrichment functions developed by Fleming et al. [27]
f
1
(r, ) =

r sin
_

2
_
(22)
f
2
(r, ) =

r cos
_

2
_
(23)
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
1472 S. LOEHNERT AND T. BELYTSCHKO
f
3
(r, ) =

r sin
_

2
_
sin() (24)
f
4
(r, ) =

r cos
_

2
_
sin() (25)
where
r =
_
(x x
tip
)
2
+ (y y
tip
)
2
(26)
=arctan
_
y y
tip
x x
tip
_
(27)
and is the angle of the crack with the x-axis at the crack tip which is located at (x
tip
, y
tip
). If the
element is completely intersected by a crack, the element nodes are enriched with the Heaviside
step function
f
1
(n) = H((n)) =
_
+1, (n)0
1, (n)<0
(28)
where (n) is a level set function so that (n) =0 gives the geometry of the crack and H() is the
Heaviside step function. The methods of Bellec and Dolbow [28] are used for very short cracks.
To simplify the notation the displacement eld is written as
u
i
h
=
n
i
n

I =1

N
i
I
u
i
I
(29)
where

N
i
I
=

N
i
I
N
i
I
f
1
(X)
.
.
.
N
i
I
f
n
enr
(X)

, u
i
I
=

u
i
I
a
i
1I
.
.
.
a
i
n
enr
I

(30)
where the superscript i indicates the respective scale.
The test functions g
i
are discretized in the same way as the displacements
g
i
h
=
n
i
n

I =1
N
i
I
_
g
i
I
+
n
enr

j =1
f
j
b
i
j I
_
(31)
Applying the same notation as for the displacements gives
g
i
h
=
n
i
n

I =1

N
i
I
g
i
I
(32)
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
A MULTISCALE PROJECTION METHOD FOR MACRO/MICROCRACK SIMULATIONS 1473
with
g
i
I
=

g
i
I
b
i
1I
.
.
.
b
i
n
enr
I

(33)
The symmetric part of the gradient of the test function is given by
grad
sym
(g
i
h
) =
n
i
n

I =1
B
i
I
g
i
I
(34)
where B
i
I
=grad
sym

N
i
I
. Accordingly, the discrete weak form of balance of momentum on the
coarse scale becomes
n
0
n

I =1
( g
0
I
)
T

__

0
B
0
I
T
: r(u
0
h
+ u
1
h
) d
. ,, .
f
int

N
0
I
f
h
d
_
*
0

N
0
I
t
h
d*
_

0
D

N
0
I
t
c
h
d
_
=0 (35)
The integration of the nodal internal forces f
int
is carried out over the entire domain including
the domain
1
. While in the domain
0
\
1
the stresses only depend on the displacement eld
u
0
h
since in
0
\
1
we have u
1
h
0, in the domain
1
the stresses depend on the displacement
eld u
1
h
=u
0
h
+ u
1
h
which are computed from the ne scale model by solving the discretized form
of (18)
n
1
n

I =1
( g
1
I
)
T

_
_

1
B
1
I
T
: r(u
1
h
) d
_

N
1
I
f
h
d
_

0
D

1
D

N
1
I
t
c
h
d
_
=0 (36)
for the displacement eld u
1
h
.
In solving the standard macroscale problem, standard Gauss quadrature methods are used in

0
\
1
except around the macrocrack, where element subdivision as described in Mo es et al. [12] is
employed. In the macroscale problem,
0

1
is integrated by taking advantage of the congruence
of the two meshes. Thus, we employ
_

0
e
() d=

S
e
_

() d (37)
where we use Gauss quadrature for each element S
e
. The subdivision of the microscale
elements around discontinuities is also carried over in the quadrature, though not indicated above.
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
1474 S. LOEHNERT AND T. BELYTSCHKO
Figure 3. Enrichments along the boundary of the ne scale domain.
For the ne scale problem, the boundary conditions (i.e. continuity conditions) are
u
1
h
=u
0
h
on *
1
(38)
Although the meshes are congruent, the enriched functions of the macro and micro elds differ.
Therefore, the continuity condition is enforced by a least-squares projection
_

1
_

N
1
I
u
1
I

N
0
J
u
0
J
_

N
1
M
g
1
M
_
d=0 (39)
Note that in (39) the integral is carried out over a subdomain

1
of the microstructure
contiguous to the interface. This is necessary to avoid non-uniquenesses of u
1
h
which can occur
when the integral is performed only over the boundary of the microstructure. Figure 3 shows
the reason for such a situation. Here, due to the macrocrack, within the ne scale domain nodes
A and B on the boundary are enriched with the jump enrichment, although no crack intersects
the boundary between the two enriched nodes. On the boundary of the ne scale domain, the
standard degrees of freedom and the enriched degrees of freedom of those nodes are linearly
dependent. Therefore, it is necessary to integrate over the domain of the ne scale element C,
which takes into account the geometry of the crack within the microstructure, to eliminate linear
dependence. Accordingly, the domain

1
is chosen to be a strip of a width of one element along
the boundary.
This projection method requires solving another linear equation system. However, this system
is much smaller than the equation system for the microscale nite element problem. Additionally,
it is possible to apply a lumping strategy for the coefcient matrix such that this operation is still
quite cheap.
Close to the crack tip of a macrocrack, the displacement eld has high gradients. For arbitrary
crack problems with microcracks around a macrocrack tip, the approximation of the displacement
eld by the macroscopic ansatz may be relatively coarse. The macroscopic solution will then be
stiffer than the solution of the microstructural problem, not only because it cannot accurately capture
the discontinuities within the microstructure, which weaken the effective stiffness of the material,
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
A MULTISCALE PROJECTION METHOD FOR MACRO/MICROCRACK SIMULATIONS 1475
but also because it simply contains a much lower number of degrees of freedom. Therefore, to
avoid an excessively stiff boundary of the microstructure domain this domain should be big enough
so that the uctuations due to the microstructure and the non-linearities due to the macrocrack are
relatively small on the boundary. This will be investigated in Section 5.
4. SOLUTION PROCEDURE
The solution procedure consists of a series of iterations in which the ne scale domain is solved
independently, followed by solution of the coarse scale equations for the entire domain. The
essence of the procedure for a linear elastic material without cohesive cracks is shown in Figure 4.
The method can easily be extended to non-linear material models and nite deformation theory
by replacing the solution step 4 by a Newton iteration step and using a geometrically consistent
linearization of the weak form and a consistent material tangent operator D in solution step 2.
However, for simplicity for the test cases investigated in Section 5 we restrict ourselves to isotropic
linear elasticity.
1. Initialize
k = 0, u
0,k
h
= 0, u
1,k
h
= 0,
0,k
h
= 0,
1,k
h
= 0
2. Solve macroscale iteration step
K
0
u
0,k
h
= f
0
ext
f
0,k
int
for u
0,k
h
with
K
0
IJ
=
0
B
0
I
T
D B
0
J
d and
f
0,k
int,I
0
B
0
I
T
: u
0,k
h
+ u
1,k
h
d +
0

N
0
I
f
h
d
set u
0,k +1
h
= u
0,k
h
+ u
0,k
h
3. Project boundary conditions (continuity conditions) from the macroscale solution onto the
boundary of the microscale domain by solving
P u
1,k +1
h
= r
k
h
for u
1,k +1
h
on the boundary
1
with
P
IJ
=

N
1
I
T

N
1
J
d and r
k
I
=

1

N
1
I
M

N
0
M
u
0,k +1
M
d
4. Solve microscale iteration step
K
1
u
1,k +1
h
= f
1,k +1
ext
f
1
int
for u
1,k +1
h
with
K
1
IJ
=
1
B
1
I
T
D B
1
J
d , f
1
int,I
=
1

N
1
I
f
h
d and
f
1,k +1
ext,I
= K
1
IJ
u
1,k +1
J
for all nodes J on the boundary of the microscale domain
5. If || u
0,k
h
|| > tol
0
or ||u
1,k +1
h
u
1,k
h
|| > tol
1
then set k k + 1 and goto step 2.
Figure 4. Solution procedure of the multiscale method.
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
1476 S. LOEHNERT AND T. BELYTSCHKO
5. NUMERICAL STUDIES
In this section, we examine a number of numerical test cases to study the different parameters
inuencing the results of the described multiscale method. In particular, those parameters are the
mesh resolution of the coarse and ne scale meshes and the domain size for the microstructure. For
all computations only square-shaped standard displacement elements are used. The stress intensity
factors, which are the quantities we investigate in all the numerical tests, are computed using a
domain form of the interaction integral as proposed in [11, 12, 29].
5.1. Inuence of the coarse scale mesh resolution
The inuence of the resolution of the coarse scale meshes on the result is investigated by considering
crack shielding of the problem shown in Figure 5. The body is loaded by prescribing mode 1
displacement boundary conditions along all edges. This example has been studied in Reference [30]
using analytical methods. The geometry is shown in Figure 5. In the coarse scale problem only
the large edge crack is considered explicitly. The inuence of the two small cracks (b =
2
135
a) is
captured only implicitly through the ne scale solution. To keep the inuence of the size of the ne
scale domain (the darkly shaded area in Figure 5) negligible, the radius r around the macrocrack tip
within which all elements are dened as multiscale elements is chosen to be constant (r =0.165a)
for all mesh resolutions. Additionally, the ne scale mesh resolution is kept constant for all coarse
scale mesh resolutions. This means that mainly the number of coarse scale elements within the
ne scale domain will change in this convergence study. The number of ne scale elements
remains almost constant. We choose the coarse scale mesh to consist of 9 9, 15 15, 27 27,
45 45 and 81 81 elements, respectively. The ne scale mesh is chosen to have a mesh size
of h =
1
1215
a.
Figure 5. Sketch of a crack shielding test: left, coarse scale mesh; and right, ne scale mesh.
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
A MULTISCALE PROJECTION METHOD FOR MACRO/MICROCRACK SIMULATIONS 1477
Figure 6. Stress intensity factor K
I
of the macrocrack tip depending on the coarse scale mesh resolution.
The results in Figure 6 show that there is hardly any inuence of the coarse scale mesh
resolution on the resulting stress intensity factor of the macrocrack. Even for a very coarse mesh
resolution the values for the stress intensity factor deviate only in the third digit. For comparison
purposes, it should be noted that the stress intensity factor of the macrocrack in the absence of
the microcracks is 1.0, so there is substantial crack shielding. The results shown in Figure 6
indicate that the far eld solution is almost independent of the detailed microstructure in the
vicinity of the macrocrack tip and that our two-scale approach can effectively capture the effects
of the microstructure. It also shows that the size of the ne scale domain is large enough for this
test and that the projection of the displacement eld on the boundary of the ne scale domain
works well.
5.2. Inuence of the ne scale mesh resolution
For the investigation of the inuence of the ne scale mesh resolution on the accuracy of the result
we use the same crack shielding example (see Figure 5) but with b =
2
27
a. The mesh resolution of
the coarse scale mesh is kept constant at 15 15 elements, and the radius r for the ne scale domain
is set to r =0.165a again. The ne scale element size varies between h =
1
45
a and h =
1
1215
a which
corresponds to a range between 225 and 164 025 elements. The results for mode 1 stress intensity
factor of the macrocrack tip are shown in Figure 7. Even for the coarsest ne scale mesh the
results are quite good. The deviation of the coarse mesh solutions from the solution for the nest
mesh resolution is less than 5%.
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
1478 S. LOEHNERT AND T. BELYTSCHKO
Figure 7. Stress intensity factor K
I
of the macrocrack tip depending on the ne scale mesh resolution.
5.3. The domain size of the microstructural domain
The size of the domain in which a microstructural computation needs to be performed has to be
chosen so that the uctuations of the eld variables on the boundary of this domain are negligible.
This means that the boundary of the ne scale domain has to be far away from the zone of inuence
of the microstructure on the macrocrack. However, in order to reduce the computational cost, it
is desirable to reduce the size of the ne scale domain as much as possible without losing too
much accuracy. Here, the domain is chosen to be of circular shape with the centre at the tip of the
macrocrack. To investigate the inuence of the ne scale domain size, the radius of the circle is
changed. The geometry for this test again is given in Figure 5 with b =
2
243
a. The coarse and ne
scale mesh element sizes are xed. The coarse scale mesh consists of 81 81 elements and
the ne scale element size is set to h =
1
2187
a. This rather ne discretization is chosen to
minimize the inuence of the mesh resolution. Figure 8 shows mode 1 stress intensity factor
as a function of the normalized radius r/a. It can be seen that even for the smallest ne scale
domain, for which the radius r is not even twice as big as the length of the microcracks, the error
of the stress intensity factor is less than 1%. For a slightly bigger radius, the error even decreases
to less than 12. This shows that the uctuations in the eld variables due to the presence of the
microcracks rapidly decay with the distance to the microcracks.
5.4. Mixed mode multiple crack problem
We consider a crack problem with many microcracks and a macrocrack under mixed mode loading.
The purpose of this example is to show the applicability of the method to random microstructures
and arbitrary loading conditions. A sketch of the geometry and loading conditions is shown
in Figure 9. The 114 microcracks, which are randomly distributed around the crack tip of the
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
A MULTISCALE PROJECTION METHOD FOR MACRO/MICROCRACK SIMULATIONS 1479
Figure 8. Stress intensity factor K
I
of the macrocrack tip depending on the radius of the ne scale domain.
Figure 9. Sketch of a mixed mode multiple crack problem.
macrocrack, are of length
1
600
a. The darkly shaded area in Figure 9 is chosen to be the ne scale
domain and is chosen to be a square of (
3
49
a) (
3
49
a). The macrocrack is at an angle of =15.64

and ends at the centre of the plate. On the coarse scale level a mesh of 49 49 elements is used,
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
1480 S. LOEHNERT AND T. BELYTSCHKO
Figure 10. Deformed conguration of the ne scale domain.
Figure 11.
yy
stress component of the ne scale domain.
while on the ne scale 261 261 elements are used to capture the quite non-smooth displacement
and stress eld in the vicinity of the macrocrack tip accurately. In Figure 10 is displayed the
deformed ne scale domain. One can see that the macro- and microcrack behaviour is captured
correctly in the whole domain and on the boundary of the ne scale domain. None of the mode 1
stress intensity factors of the microcracks is negative which indicates that crack closure does not
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
A MULTISCALE PROJECTION METHOD FOR MACRO/MICROCRACK SIMULATIONS 1481
occur. Figure 11 shows the
yy
component of the stresses in the ne scale domain. The mode 1
stress intensity factor for the macrocrack computed by the macro/micro model is K
I
=19.836. For
a coarse scale mesh resolution of 147 147 elements, the stress intensity factor for the macrocrack
is 19.924. In the absence of the microcracks, the stress intensity factor is K
I
=28.163. We have
not been able to do a single-scale computation for this problem because the microcracks are very
short compared to the scale of the structure.
6. CONCLUSIONS
A method has been presented for the multiscale analysis of bodies with extensive microstructure
on a scale signicantly smaller than features of the macrostructure. This scale separation makes it
possible to solve the microproblem separately from the macroproblem, thus providing substantial
savings in computational cost. The scale separation is affected by using as test functions the coarse
scale functions on the macroscale, but retaining the microscale features in the stress eld. In this
paper, the focus has been on microcrack/macrocrack interactions. As shown in the development,
only the cohesive tractions from the macrocracks appear in the coarse scale equations, whereas
the macrocrack and microcrack cohesive tractions appear in the ne scale equations. The proposed
multiscale method for crack simulations considerably improves the efciency of simulations of
microcrack/macrocrack interaction problems.
We have employed the extended nite element method for crack modelling. This facilitates
the treatment of complex crack patterns. However, the method is applicable to other multiscale
approaches, such as combinations of molecular and continuum mechanics.
The numerical examples demonstrate the excellent convergence behaviour of the method in the
interplay of coarse scale and ne scale behaviour, namely, the stress intensity factor as inuenced
by the microcracks. The effectiveness of the method has been demonstrated entirely in terms of
examples. At this time, we have no mathematical basis for how well and under what circumstances
the method works. However, it is apparent from the examples that when there is signicant
separation of scales, a two-scale method such as proposed here is quite effective. Additionally, this
method is applicable to a large variety of other problems, such as dislocations and shear bands.
An advantage of this method is the ease of its implementation, for it requires little modication
of standard nite element software. The method can also easily be extended to three-dimensional
problems and other applications such as heterogeneities.
ACKNOWLEDGEMENTS
The support of the Army Research Ofce under grant W911NF-05-1-0049 is gratefully acknowledged.
REFERENCES
1. Hughes TJR. Multiscale phenomena: Greens functions, the Dirichlet-to-Neumann formulation, subgrid scale
models, bubbles and the origins of stabilized methods. Computer Methods in Applied Mechanics and Engineering
1995; 127:387401.
2. Zohdi TI, Oden JT, Rodin GJ. Hierarchical modeling of heterogeneous bodies. Computer Methods in Applied
Mechanics and Engineering 1996; 138:273298.
3. Oden JT, Zohdi TI. Analysis and adaptive modeling of highly heterogeneous elastic structures. Computer Methods
in Applied Mechanics and Engineering 1997; 148:367391.
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme
1482 S. LOEHNERT AND T. BELYTSCHKO
4. Feyel F, Chaboche JL. Fe
2
multiscale approach for modelling the elastoviscoplastic behaviour of long bre sic/ti
composite materials. Computer Methods in Applied Mechanics and Engineering 2000; 183:309330.
5. Farhat C. A Lagrange multiplier based divide and conquer nite element algorithm. Journal of Computers in
System Engineering 1991; 2:149156.
6. Xiao SP, Belytschko T. A bridging domain method for coupling continua with molecular dynamics. Computer
Methods in Applied Mechanics and Engineering 2004; 193:16451669.
7. Ben Dhia H, Rateau G. The Arlequin method as a exible engineering design tool. International Journal for
Numerical Methods in Engineering 2005; 62:14421462.
8. Terada K, Kikuchi N. A class of general algorithms for multi-scale analyses of heterogeneous media. Computer
Methods in Applied Mechanics and Engineering 2001; 190:54275464.
9. Karpov EG, Wagner GJ, Liu WK. A Greens function approach to deriving non-reecting boundary condition
in molecular dynamics simulations. International Journal for Numerical Methods in Engineering 2005; 62:
12501262.
10. Fish J, Markolefas S, Guttal R, Nayak P. On adaptive multilevel superposition of nite-element meshes for linear
elastostatics. Applied Numerical Mathematics 1994; 14:135164.
11. Belytschko T, Black T. Elastic crack growth in nite elements with minimal remeshing. International Journal
for Numerical Methods in Engineering 1999; 45(5):601620.
12. Mo es N, Dolbow J, Belytschko T. A nite element method for crack growth without remeshing. International
Journal for Numerical Methods in Engineering 1999; 46(1):131150.
13. Budyn E, Zi G, Mo es N, Belytschko T. A method for multiple crack growth in brittle materials without remeshing.
International Journal for Numerical Methods in Engineering 2004; 61:17411770.
14. Guidault PA, Allix O, Champeney L, Cornault C. A multiscale extended nite element method for crack
propagation. Computer Methods in Applied Mechanics and Engineering 2006, in press.
15. Stolarska M, Chopp DL, Mo es N, Belytschko T. Modelling crack growth by level sets in the extended nite
element method. International Journal for Numerical Methods in Engineering 2001; 51:943960.
16. Belytschko T, Mo es N, Usui S, Parimi C. Arbitrary discontinuities in nite elements. International Journal for
Numerical Methods in Engineering 2001; 50:9931013.
17. Xiao QZ, Karihaloo BL. Improving the accuracy of XFEM crack tip elds using higher order quadrature and
statically admissible stress recovery. International Journal for Numerical Methods in Engineering 2006; 66:
13781410.
18. Laborde P, Pommier J, Renard Y, Sala un M. Higher-order extended nite element method for cracked domains.
International Journal for Numerical Methods in Engineering 2005; 64:354381.
19. Mo es N, Gravouil A, Belytschko T. Non-planar 3d crack growth by the extended nite element and level sets
Part I: Mechanical model. International Journal for Numerical Methods in Engineering 2002; 53:25492568.
20. Gravouil A, Mo es N, Belytschko T. Non-planar 3d crack growth by the extended nite element and level sets
Part II: Level set update. International Journal for Numerical Methods in Engineering 2002; 53:25692586.
21. Areias PMA, Belytschko T. Analysis of three-dimensional crack initiation and propagation using the extended
nite element method. International Journal for Numerical Methods in Engineering 2005; 63:760788.
22. Areias PMA, Belytschko T. Two-scale shear band evolution by local partition of unity. International Journal for
Numerical Methods in Engineering 2006; 66:878910.
23. Samaniego E, Belytschko T. A continuumdiscontinuum modelling of shear bands. International Journal for
Numerical Methods in Engineering 2005; 62:18571872.
24. Ventura G, Moran B, Belytschko T. Dislocations by partition of unity. International Journal for Numerical
Methods in Engineering 2005; 62:14631487.
25. Gracie R, Ventura G, Belytschko T. A new fast nite element method for dislocations based on interior
discontinuities. International Journal for Numerical Methods in Engineering 2006. DOI: 10.1002/nme.1896
26. Bao G, Suo Z. Remarks on crack bridging concepts. Applied Mechanics Reviews 1992; 45:355366.
27. Fleming M, Chu YA, Moran B, Belytschko T. Enriched element-free Galerkin methods for crack tip elds.
International Journal for Numerical Methods in Engineering 1997; 40:14831504.
28. Bellec J, Dolbow J. A note on enrichment functions for modelling crack nucleation. Communications in Numerical
Methods in Engineering 2003; 19:921932.
29. Moran B, Shih CF. Crack tip and associated domain integrals from momentum and energy balance. Engineering
Fracture Mechanics 1987; 27(6):615642.
30. Rose LRF. Microcrack interaction with a main crack. International Journal of Fracture 1986; 31:233242.
Copyright q 2007 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:14661482
DOI: 10.1002/nme