A micromechanics-enhanced finite element formulation for
modelling heterogeneous materials
Jan Nováka,b,c , Lukasz Kaczmarczyka , Peter Grassla , Jan Zemanb , Chris J. Pearcea,∗

arXiv:1103.5633v2 [cond-mat.mtrl-sci] 15 Jun 2011

a

School of Engineering, University of Glasgow, Glasgow G12 8QQ, UK
Faculty of Civil Engineering, Czech Technical University in Prague, Thákurova 7,
166 29 Praha 6, Czech Republic
c
School of Civil and Environmental Engineering, University of New South Wales, NSW 2052,
Sydney, Australia
b

Abstract
In the analysis of composite materials with heterogeneous microstructures, full resolution of the heterogeneities using classical numerical approaches can be computationally prohibitive. This paper presents a micromechanics-enhanced finite element
formulation that accurately captures the mechanical behaviour of heterogeneous materials in a computationally efficient manner. The strategy exploits analytical solutions derived by Eshelby for ellipsoidal inclusions in order to determine the mechanical perturbation fields as a result of the underlying heterogeneities. Approximation
functions for these perturbation fields are then incorporated into a finite element
formulation to augment those of the macroscopic fields. A significant feature of
this approach is that the finite element mesh does not explicitly resolve the heterogeneities and that no additional degrees of freedom are introduced. In this paper,
hybrid-Trefftz stress finite elements are utilised and performance of the proposed formulation is demonstrated with numerical examples. The method is restricted here
to elastic particulate composites with ellipsoidal inclusions but it has been designed
to be extensible to a wider class of materials comprising arbitrary shaped inclusions.
Keywords: Micromechanics; Equivalent inclusion method; Eshelby’s solution;
Heterogeneous materials; Hybrid-stress finite elements; Displacement perturbations

∗

Corresponding author. Tel.: +44-141-330-5207
Email addresses: novakj@cml.fsv.cvut.cz (Jan Novák),
Lukasz.Kaczmarczyk@glasgow.ac.uk (Lukasz Kaczmarczyk), Peter.Grassl@glasgow.ac.uk
(Peter Grassl), zemanj@cml.fsv.cvut.cz (Jan Zeman), Chris.Pearce@glasgow.ac.uk (Chris J.
Pearce)

Preprint submitted to Computer Methods in Applied Mechanics and Engineering

June 16, 2011

1. Introduction
In the analysis of materials with complex microstructures, full resolution of
the heterogeneities using classical numerical approaches such as the Finite Element
method can be computationally prohibitive. To overcome this, one option is to model
the macroscale problem using equivalent properties; however, this can lead to a critical loss of information about the finer scale behaviour and poor understanding of the
heterogeneities’ influence on the macroscale response. Numerical approaches such
as computational homogenization (often called FE2 ) provide an alternative strategy [1, 2]. These techniques comprise nested Finite Element analyses, where each
macroscopic material point response is determined via the numerical solution of an
RVE subject to the macroscopic strains. Although such approaches have significant
potential for certain classes of problems, they are still computationally demanding
and are restricted to situations involving clear separation of scales.
The objective of this work is to develop a Finite Element formulation for modelling
the macroscopic mechanical problem that is enhanced to capture the influence of the
underlying heterogeneities. In our approach, the Finite Element mesh is not required
to explicitly resolve the heterogeneities. Closed-form expressions for the perturbation
of the mechanical fields due to the presence of the heterogeneities are determined
and these are then utilised to enhance the Finite Element formulation.
The ability to capture the effect of microstructural features independently of the
underlying finite element mesh has been an ongoing challenge in computational mechanics research. Partition of Unity methods [3, 4, 5, 6] provide a potential solution
to this problem, without mesh refinement, by extending a given solution space with
additional functions and has been successfully applied to problems such as cracks
and material interfaces. The application of this approach in the context of the current work will be briefly discussed in this paper, whereby the closed-form solutions
derived for the mechanical perturbation fields are used to extend the classical finite
element method. However, it will be shown that there are some disadvantages to this
approach for the particular problem at hand and an alternative approach, centred
on the Hybrid-Trefftz stress element formulation [7], represents the main focus of
this paper. This method does not result in additional degrees of freedom, although
it does involve an additional, albeit relatively minor, computational overhead.
The heterogeneities, although currently restricted to simple shapes (ellipsoids),
can be randomly sized and randomly distributed without reference to the finite element mesh. Therefore, the proposed approach has the potential to be applicable to a wide range of composite materials, such as fibre reinforced composites [8],
porous media [9, 10], functionally graded materials [11], etc. Moreover, it can be

2

extended to general inclusion shapes by evaluating the perturbation functions numerically [12, 13].
The paper is structured as follows. The methodology of the proposed strategy is
described in Section 2. Construction of the perturbation approximation functions
for Finite Element Analysis is derived in Section 3. The implementation into the
Hybrid-Trefftz stress element formulation containing an arbitrary number of inclusions is presented in Section 4. Section 5 comprises examples demonstrating the
model’s performance. Finally we present the conclusions as well as a discussion on
future research directions. An appendix is included that highlights some important,
but rather technical, aspects of the proposed technique in order to keep the paper
self-contained.
2. Micromechanics approach
This section outlines the strategy to calculate the perturbation of mechanical
fields due to a heterogeneous microstructure, exploiting the Equivalent Inclusion
Method [14] in conjunction with analytical micromechanics. Our objective is to
convert the heterogeneous problem into an equivalent homogeneous problem and to
derive analytical expressions for the perturbations of the stress, strain and displacement fields that we can then utilise within a finite element formulation.
Consider a body consisting of clearly distinguishable heterogeneities in a matrix
(Fig. 1a) subjected to a displacement u and traction t field. The stiffness of such a
u (x)

, t (x)

Γ

Ω
C1 (x)

6
C (x)

y

Ω4
C (x)

≈

Ω3

Ω

Ω :C

Γ

Ω
C5 (x)
5

C

C0
ε2τ

0

ε4τ

Ω2

Ω4

C0
ε τ3
Ω3

C0

Ω5

6

C0
ε τ5

4

0

x

ε τ1

C3
(x)

Ω

, t (x)

Ω1

Ω2
C2 (x)

1

6

u (x)

C 0 ετ

5

y

0

0

x

(a)

Ω :C

0

(b)

Figure 1: Principle of Equivalent Inclusion Method: a) composite body with inclusions, b) homogeneous reference body with additional equivalent eigenstrains

material is decomposed as follows [14, 13]
C = C0 + V C∗
3

(1)

P
i
0
where C0 is the stiffness tensor of a homogeneous matrix Ω0 and C∗ = N
i [C − C ]
is due to the presence of N inclusions. C∗ is nonzero only within the domain Ω =
Ω1 ∪ · · · ∪ ΩN , so that

0 in Ω0
(2)
V =
1 in Ω
As a result of the heterogeneities, the mechanical fields (displacement, strains, stresses)
experience a perturbation for which we will derive closed-form expressions based on
analytical micromechanics. Symbolically, we can express the decomposition of the
mechanical fields as follows:
u = u0 + u∗ , ε = ε0 + ε∗ , σ = σ 0 + σ ∗

(3)

where, the superscript ·0 indicates the macroscopic component of the fields in the absence of heterogeneities and superscript ·∗ indicates the perturbation (or microscopic)
component due to the presence of the heterogeneities. It is worth noting that, traditionally, in analytical micromechanics, the macroscopic fields are assumed to be
uniform across the domain, e.g. [15, 16]. Here it is assumed that they can be position
dependent functions of the Neumann and Dirichlet boundary conditions.
The perturbation fields are determined by employing the equivalent inclusion
method for a single heterogeneity embedded in a matrix and then extended here
for multiple heterogeneities. In the equivalent inclusion method, the heterogeneous
solid is replaced by an equivalent homogeneous solid with uniform material stiffness
C0 everywhere (Fig. 1a, b) and suitable stress-free eigenstrains ετi applied in the
inclusions Ωi so that the homogeneous equivalent solid has the same mechanical
fields as the original heterogeneous solid.
2.1. Equivalent inclusion method for single heterogeneity problem
Consider first a single heterogeneity embedded in a matrix. Following Eshelby’s
fundamental work [17], this problem can be decomposed into two problems of known
solution and then assembled back via superposition [14, 17], see Fig. 2. In brief,
the solution of the inhomogeneity problem requires the determination of the transformation eigenstrain ετ that induces the identical local mechanical response as the
original heterogeneous body. Due to the absence of other inclusions, no strain or
stress concentrations are induced and thus ετ remains constant [17].
In the original heterogeneous body, the stress state can be expressed as
σ = σ 0 + σ ∗ = C : [ε0 + ε∗ ]

4

(4)

u 0(x) ,q0(x)

(a)

u (x) , t (x)

Γ
x

8

ετ
Ω: C 0
Ω0: C 0
−

8

−

u 0(x) ,q0(x)

everywhere: C 0
−

+

y

8

8

−

8

8

Ω0: C 0
−

≡

x

u (x) , t (x)

Ω: C1 (x)

y

8

8

u (x) , t (x)

−

8

x

u (x) , t (x)

u (x) , t (x)

Γ

y

8

8

8

u (x) , t (x)

u 0(x) ,q0(x)

u (x) , t (x)

u (x) , t (x)

(b)

u 0(x) ,q0(x)

(c)

Figure 2: Equivalent Inclusion Method: a) inhomogeneity problem, b) infinite homogeneous body,
c) homogeneous inclusion problem

In the homogeneous solid, we add a stress-free eigenstrain ετ inside the domain of
the inclusion, which has the same material stiffness C0 as the matrix such that
σ = C0 : [ε0 + ε∗ − ετ ]

(5)

It should be noted that ετ = 0 in the matrix.
Given that the macroscopic stress is σ 0 = C0 : ε0 , it can be see from equating Eqs. (4) and (5), that the stress perturbation in the homogeneous solid can be
expressed as
σ ∗ = C0 : [ε∗ − ετ ]
(6)
Furthermore, equating Eqs. (4) and (5) also results in the following expression:
C : [ε0 + ε∗ ] = C0 : [ε0 + ε∗ − ετ ]

(7)

where the transformation eigenstrain is as yet unknown. Eshelby’s solution of the
homogeneous inclusion problem [17], relates the eigenstrain to the perturbation strain
as follows
ε∗ = S : ετ
(8)
where S denotes the Eshelby tensor and is a function of the heterogeneity’s geometry
and the material stiffness of both the matrix and heterogeneity. Substituting this
expression into Eq. (7) and rearranging yields




C − C 0 : ε0 = C 0 : S − C : S − C 0 : ετ
(9)
This can be recast in a compact form to give an expression for the eigenstrain ετ
which depends on the homogeneous strain ε0 , the material stiffness of both the matrix
and heterogeneity and the Eshelby tensor as
ετ = B : ε0
5

(10)

where the tensor B is provided by:

−1
B = − C∗ : S + C0
: C∗

(11)

Once the transformation eigenstrain has been determined, the stress perturbation
can be computed from Eq. (6) in the form
σ ∗ = C0 : [S − I] : B : ε0

(12)

It can be seen that the stress perturbation depends on stiffness of the different
material phases, the macroscopic strain field and the geometry of the heterogeneity.
This closed-form expression for the stress perturbation is at the heart of the proposed
finite element enrichment to be discussed later in this paper. It is also useful to derive
an expression for the displacement perturbation field as follows
u∗ = L : ετ = L : B : ε0

(13)

where the operator L is a third order tensor, mapping ετ → u∗ . For the sake
of conciseness and to keep the paper self-contained, the detailed derivation of this
operator can be found in Appendix Appendix A.
2.2. Self-compatibility algorithm for multiple inclusions
In the case of multiple inclusions, the mechanical perturbation fields within individual inclusion domains are no longer uniformly distributed as a result of their
mutual interaction. Here, we account for this approximately by assuming that the
eigenfields to be piecewice constant within each inclusion. Thus, the perturbations
are determined from the Eshelby solution for each individual inclusion, as described
above, plus an iterative self-compatibility procedure (Tab. 1) to ensure that the
solution correctly reflects the influence of multiple heterogeneities.
This procedure iteratively enforces compatibility (see eg. [18] for further reference) between the imposed macroscopic strain and the average eigenstrain inside
any given inclusion i, so as to account for the influence of the remaining inclusions
N \i. An iterative algorithm has been chosen because a closed form solution for
multiple inclusions does not exist and a numerical solution would be prohibitively
expensive [18].
First, the eigenstrain ετi for each inclusion i is calculated (Eq. 10) without reference to the other inclusions (Line 2). Next, the associated perturbation strain ε∗i
for each inclusion i is evaluated (Eq. 8) at the centre of all other inclusions (Line
3). The mutual interaction of inclusions is then taken into account via a correction
of the eigenstrain (∆ετi ). For each inclusion i, this correction is calculated (Line 8)
6

from the inverse of the inclusion’s Eshelby tensor S−1
and the perturbation strains
i
of all other inclusions, evaluated at the centroid of inclusion i. The perturbation
strain resulting from inclusion j at the centre of inclusion i is denoted as ε∗i,j . This is
demonstrated in Fig. 3 for a two inclusion problem in 1D. The eigenstrain correction
∆ετi is then used to calculate the correction to the perturbation strains (Line 10).
The algorithm continues until a small Euclidean norm between the last two iterations of the total eigenstrains is achieved. At convergence, the corresponding stress
and displacement perturbations are recalculated from the corrected transformation
eigenstrains.
It is worthwhile noting that the algorithm does not depend on a particular sequence of inclusions, as follows from the elastic reciprocity theorem [18, and references therein]. Moreover, since the perturbation fields are calculated for the entire
macroscopic domain (no RVE is considered), stress admissibility and strain compatibility, in the sense of macro-micro field relations, are fulfilled a priori.
The computational complexity of the Self-compatibility algorithm is O(N 2 ). However, this can be improved by taking into account only those inclusions which have
a non-negligible influence to the inclusion of interest i. Preliminary studies have
shown this to give a significant computational speed-up and will be reported in a
future paper.

1
2
3
4
5
6
7
8
9
10
11
12

Self Compatibility Algorithm (ε0i , Bi , Si , S−1
i , N)
For (i ≤ N )
ετi = Bi : ε0i (Eq. 10)
ε∗i = Si : ετi (Eq. 8)
Set ∆ε∗i = ε∗i
EndFor
Do
For (i ≤ N )
P
−1
∗
(Eq. 8)
∆ετi = N
j\i Si : ∆εi,j
τ
τ
τ
εi = εi + ∆εi
∆ε∗i = Si : ∆ετi
EndFor

PN
∗
While
i k∆εi k > η

Table 1: Self-compatibility algorithm. Note, that η stands for an acceptable tolerance.

7

Figure 3: Principle of self-compatibility algorithm for double inclusion problem in 1D.

3. Construction of perturbation approximation functions for FEA
The above methodology can be utilised to formulate an enhanced Finite Element
formulation. The primary task is to determine appropriate approximation functions
for the mechanical perturbation fields u∗ , ε∗ and σ ∗ based on the analytical micromechanics developed above and which can then augment the standard macroscopic field
approximations. It should be noted that the Voigt-Mandel notation is exclusively
used in the forthcoming text.
The perturbation field approximation functions are determined a priori as a linear
combination of the perturbation fields evaluated analytically for six load cases, with
self-equilibrium enforced by means of the self-compatibility algorithm outlined above
(Tab. 1). Each load case corresponds to a unit component of the macroscopic strain
vector 0i , i = 1, . . . , 6 and the resulting analytically determined stress, strain and
displacement perturbation fields are arranged, column-by-column into s∗6×6 , e∗6×6 and
u∗3×6 matrices, respectively:
s∗ =



1 ∗

σ

...
u∗ =

6 ∗

σ





, e∗ =

1 ∗

u



1 ∗

 ...

. . . 6 u∗

6 ∗





(14)
(15)

where the left superscript refers to a specific load case, 1 to 6.
3.1. Partition of unity method
Partition of Unity (PU) Methods (for example the eXtended Finite Element
Method) extend the underlying basis functions used for interpolating the displacement field by adding an appropriate set of additional functions. Following [19, 20] it
has been shown that the displacement field u(x) within an element can be interpolated by
n
X

u(x) =
Ni (x)ai + Ni Nγ (x)bi
(16)
i=1

8

where n is the number of nodes per element, Ni (x) = N i (x)I is the standard matrix
of element shape functions for node i, I is the identity matrix and ai the standard
displacement degrees of freedom at node i. Nγ is a matrix containing the additional
basis terms and bi are the associated additional degrees of freedom at node i. It is
important to recognise that the element shape functions form a partition of unity,
i.e.
n
X
N i (x) = 1
(17)
i=1

The six analytically derived displacement perturbation functions contained in u∗ can
be used as the additional functions Nγ to augment the standard basis functions. Thus

 1 ∗
u1 (x) . . . 6 u∗1 (x)
0
...
0
0
...
0
1 ∗
 (18)
0
...
0
0
...
0
u2 (x) . . . 6 u∗2 (x)
Nγ (x) = 
1 ∗
0
...
0
0
...
0
u3 (x) . . . 6 u∗3 (x)
With this at hand, the PU-based finite element formulation can be derived, see for
example [21]. PU methods are particularly useful in problems where the extension of
the basis functions is introduced on a node by node basis, so that additional degrees
of freedom are only introduced at nodes where the basis is extended. One obvious
example of such a local feature that can be modelled in this way is discrete cracks [21].
However, this favourable property is not exploited here because we wish to model a
large number of heterogeneities throughout the domain. In 3D problems, there are
3 standard displacement degrees of freedom per node; this would be extended by an
additional 18 degrees of freedom per node and per heterogeneity with the proposed
approach.
It is also worth noting that for standard finite elements, the volume integration
of the discrete system of equations is relatively straightforward. However, extension
of the basis functions to include the perturbation functions in Eq. (18) makes this
process significantly more arduous. For these reasons, an alternative Finite Element
approach using Hybrid Trefftz Stress elements [22, 7] is considered where the standard
basis function is not extended, as with PU methods, but enhanced such that no
additional degrees of freedom are introduced. This is described in the next section.
4. Hybrid-Trefftz stress element formulation
In this section a finite element formulation based on an enhancement of a hybridTrefftz stress (HTS) element formulation [7] is presented.

9

t (x)

Γe

Γσe
u3

Ωe
z

x

u1

Γue

u(x)

u2

y
x

Figure 4: Elastic body representing HTS element

The problem requires a solution to the displacement u and stress σ fields as a
result of given boundary displacements u and tractions q on Γeu and Γeσ , respectively.
The displacement and stress fields must fulfil the following governing equations:
LT σ
Lu
σ
Nσ
u

=0
=
= C
=t
=u

in Ωe
in Ωe
in Ωe
on Γeσ
on Γeu

. . . Cauchy equilibrium equation
. . . strain-displacement relationship
. . . constitutive equation
. . . static boundary conditions
. . . kinematic boundary conditions

(19)

where σ and  are the column matrix representation of the second order stress and
strain tensor, respectively, u represents the displacement vector, C is the matrix
representation of fourth order stiffness tensor and finally u and t represent the applied
displacements and tractions, respectively. The gradient operator L and the matrix
of directional cosines N of the outward normal to element boundary Γe have the
following forms [23]


∂/∂x
0
0
 0


∂/∂y
0 


n
0
0
0
n
n
x
z
y
 0
0
∂/∂z 
 , N =  0 ny 0 nz 0 nx 
L=
(20)
 ∂/∂y ∂/∂x
0 


0 0 nz ny nx 0
 0
∂/∂z ∂/∂y 
∂/∂z
0
∂/∂x
4.1. Stress, strain and displacement approximations
The macroscopic stress field within the HTS element is approximated as
σ 0 = S0v v in Ωe
10

(21)

where v is the vector of generalised stress degrees of freedom, S0v denotes the matrix of
stress approximation functions chosen so as to automatically satisfy the equilibrium
conditions Eq. (19)1,4 . Thus,
LT S0v v = 0 in Ωe
(22)
and
t = NS0v v on Γe

(23)

where t represents the traction vector induced by the macroscopic stress approximation field.
The macroscopic strain and displacement fields are expressed analogously to
Eq. (21) as
(24)
0 = E0v v and u0 = U0v v
where E0v and U0v are directly associated with the stress approximation by means of
the compatibility equation (19)2 and constitutive equation (19)3 as
S0v = C0 E0v = C0 LU0v

(25)

Since the stress approximation functions S0v are typically polynomial functions, the
integration of E0v to get U0v is relatively straightforward.
Rather than extend the solution space to capture the influence of the heterogeneities, as was briefly described in Section 3.1, here we enhance the macroscopic
approximations to include the influence of the heterogeneities, thereby not increasing the number of unknowns. The total stress (macroscopic plus perturbation) field
within the HTS element is approximated, following Eq. (3), as

σ = σ 0 + σ ∗ = S0v + S∗v v in Ωe

(26)

where S∗v is the perturbation counterpart to S0v . S∗v can be constructed from Eq. (14):
σ ∗ = s∗ 0

(27)

where s∗ is the set of analytically defined stress perturbations for six load cases,
each one representing a unit component of the macroscopic strain Eq. (14). For
the purposes of constructing S∗v , we approximate the macroscopic strain field as
constant within each finite element and computed as the volume average of the
actual macroscopic strain field. From Eq. (24),
Z
Z
1
1
0,ave
0
e

= e
 dΩ = e
Ev dΩe v = Eave
(28)
v v
|Ω | Ωe
|Ω | Ωe
11

Substituting 0,ave for 0 into Eq. (27) leads to
σ ∗ = s∗ Eave
v v

(29)

Thus, the matrix of stress perturbation approximation functions is:
S∗v = s∗ Eave
v

(30)

Analogously, the approximation of total strain and displacement fields within the
element domain is given by


(31)
 = E0v + E∗v v and u = U0v + U∗v v
where the perturbation approximation matrices U∗v and E∗v are, as with their macroscopic counterparts, directly associated with the stress approximation as
S∗v = C0 E∗v = C0 LU∗v

(32)

It is worthwhile noting that the stress perturbation fields and the corresponding
traction perturbation fields, approximated as σ ∗ = S∗v v and t∗ = NS∗v v, remain in
self-equilibrium.
4.2. Static boundary conditions
Contrary to general condition in Eq. (19)4 , the equilibrium on the element traction
boundary is imposed only in the weighted residual sense as:
Z



W1T N σ 0 + σ ∗ − t dΓe = 0
(33)
Γeσ

along with W1 representing the matrix of weighting functions. Replacing the total
stress field in Eq. (33) by its approximation from Eq. (26), the traction boundary
condition becomes
Z
Z

T
0
∗
e
W1 N Sv + Sv v dΓ =
W1T t dΓe
(34)
Γeσ

Γeσ

4.3. Kinematic boundary conditions
Compatibility inside the element domain Ωe is also enforced in a weighted residual
sense, such that:
Z

W2T 0 + ∗ − Lu dΩe = 0
(35)
Ωe

12

Next, utilising integration by parts and applying Green’s theorem to W2T Lu, Eq. (35)
results in
Z
Z

e
∗
0
T
(LT W2 )T u dΩe
W2 Ev + Ev v dΩ +
e
e
Ω
Z
ZΩ
T
e
(NW2 ) u dΓ =
(NW2 )T u dΓe
−
(36)
Γeσ

Γeu

With the current formulation, it is not necessarily possible to find a solution to both
Eqs. (34) & (36) that satisfies both traction and kinematic boundary conditions acting on the element boundary. As a consequence, Eq. (36) is relaxed by introducing an
additional and independent approximation of displacements on the element traction
boundary:
uΓ = UΓ q in Γeσ
(37)
Here, q stands for the set of displacement unknowns and UΓ is the matrix of boundary
displacement approximation functions. Such a formulation of the stress element leads
to a hybrid approach [7, 22].
Given the above consideration, introducing Eq. (37) into Eq. (36) results in
Z
Z

T
0
∗
e
W2 Ev + Ev v dΩ +
(LT W2 )T u dΩe
e
Ωe
Z
ZΩ
T
e
(NW2 )T u dΓe
(NW2 ) uΓ dΓ =
−
(38)
Γeu

Γeσ

4.4. Weighting functions
In order to achieve an energy-consistent formulation, it is required that all weighted
terms within the integrals defined above have the dimension of work. The weighting
functions then directly follow from the integrands in Eq. (33) and Eq. (35) representing the increment of internal work of strains within Ωe and external work of tractions
on Γeσ , respectively. These functions thus admit the following forms:
W2 = S0v

W1 = UΓ ,

(39)

First, introducing Eq. (39)2 into Eq. (38) and taking into account condition (22)
yields
Z
Z
Z



T
0 T
0
∗
e
0 T
e
(40)
Sv
Ev + Ev v dΩ −
NSv UΓ q dΓ =
NS0v u dΓe
Ωe

Γeσ

Γeu

13

Second, introducing Eq. (39)1 into the traction boundary condition (34) yields
Z
Z
e
∗
0
T
UTΓ t dΓe
UΓ N(Sv + Sv )v dΓ =
(41)
Γeσ

Γeσ

Combining Eqs. (40) & (41) results in a coupled system of linear equations that can
be expressed in compact form as

  

F
−AT
v
pu
=
(42)
−(A + A∗ )
0
q
−pσ
where the submatrices on the left-hand side follow from Eqs (40, 41) and are given
by the following integrals
Z
Z


T 0

0 T
0
∗
e
Sv
Ev + Ev dΩ =
N S0v
Uv + U∗v dΓe
(43)
F =
e
e
ZΩ
ZΓ
UΓ NS0v dΓe and A∗ =
UΓ NS∗v dΓe
(44)
A =
Γeσ

Γeσ

and for the terms on the right-hand side it holds
Z
Z

0 T
e
pu =
NSv u dΓ , and pσ =
Γeu

Γeσ

UTΓ t dΓe

(45)

Note that Eq. (43) illustrates that the F matrix can be evaluated via a boundary
rather than volume integral. Thus all terms in Eq. (42) can be evaluated using
boundary integrals only.
The size of the system of equations to be solved simultaneously can be reduced
via static condensation, representing a significant reduction in computational effort.
First, from the first equation in Eq. (42), the generalised stress degrees of freedom v
are expressed in terms of the displacement degrees of freedom q as
v = F−1 (pu + AT q)
This is then substituted into the second equation of Eq. (42) to yield:


(A + A∗ )F−1 AT q = pσ − (A + A∗ )F−1 pu

(46)
(47)

This sparse system of equations is then solved for the displacement degrees of freedom
q. Subsequently, the stress degrees of freedom v can then be calculated on an elementby-element basis.
Our implementation of these HTS elements for composite materials (C-HTS elements) utilises displacement degrees of freedom that are associated with element
faces rather than vertices. This has the advantage that the bandwidth of the stiffness
matrix is minimised, as is interprocessor communication.
14

5. Numerical Examples
The performance of the key components of the proposed strategy (micromechanical solution, self-compatibility algorithm, finite element analysis convergence, etc.),
in terms of efficiency and accuracy, have been explored through two numerical examples.
5.1. Three ellipsoidal inclusions in matrix
This example comprises three ellipsoidal inclusions embedded in a cube of matrix.
The geometry of the problem is illustrated in Fig. 5 and details of the ellipsoids are
given in Table 2, including the semi-axes’ dimensions, centroid coordinates and Euler
angles φ, ν and ψ, which are successive rotations of the semi-axes a1 , a2 and a3 about
global coordinate axes z, x and z, respectively. The cube has side lengths of 600.
The displacement boundary conditions were prescribed on faces x = 300, y = 300
and z = 300 as ux = 0, uy = 0 and uz = 0, respectively. The remaining faces at
x = −300, y = −300 and z = −300 were subject to uniform normal unit tractions.
The Young’s modulus for the homogeneous matrix was chosen as E = 1 and for the
heterogeneities as E = 2. Poisson’s ratio was chosen as ν = 0.1 for both materials.
All units are consistent.
It is worthwhile noting that the small contrast in stiffness between the two materials was chosen deliberately to maximise the extent of the perturbation fields
emanating from the heterogeneities. Large contrasts in stiffness between the matrix and heterogeneities lead to perturbation fields that decay rapidly with distance
from the heterogeneities. The close proximity of the three ellipsoidal heterogeneities
to each other was also chosen deliberately in order to demonstrate the ability of
the formulation to capture the interaction of multiple heterogeneities. Furthermore,
the close proximity of one of the ellipsoids to a traction boundary was chosen to
demonstrate the ability of the formulation to capture boundary effects.
Incl.
1
2
3

Centroid coordinates
x
y
z
-48.07 78.27
14.81
16.45 178.64 -154.51
127.93 -65.94 -27.32

Semiaxes dim.
a1
a2
a3
50
75 100
50 100 75
100 75
50

Euler
φ
74.21
37.27
46.74

angles of
ν
48.44
22.27
11.17

max ai
ψ
-48.07
-25.51
-26.30

Table 2: Topology and geometry of ellipsoidal inclusions of triple inclusion problem

The problem was analysed using two three-dimensional finite element meshes
with different densities, comprising C-HTS elements. The coarse mesh comprised
15

z
y
x

Figure 5: Geometry and topology of triple inclusion problem

24 elements and 540 DOFs (Fig. 6a) and the second, refined, mesh comprised 192
elements and 3888 DOFs (Fig. 6b). Results from the two enhanced finite element

(a)

(b)

Figure 6: Triple inclusion task discretization by C-HTS elements: a) Coarse mesh with 24 enhanced
elements (540 DOFs) b) Finer mesh with 1,536 enhanced elements (29,376 DOFs)

analyses are plotted in the yz-plane (at x = 0). Fig. 7a and Fig. 8a show the two
meshes in this plane. The σyy stress component for both analyses are shown in
Fig. 7b and Fig. 8b.
In addition, a reference finite element analysis of the same problem was undertaken for comparison sake. The reference analysis utilised HTS elements but without
the proposed enhancement. Unlike the other two analyses, the reference analysis
utilised a mesh that explicitly resolved the three ellipsoidal heterogeneities and comprised 309, 406 tetrahedrons with 5, 596, 776 DOFs (Fig. 9a). The corresponding
mesh and stress results of the reference analysis are shown in Fig. 9. Comparison of
the stress results from the enhanced formulation and the reference analysis leads to
the relative error plots shown in Fig. 7c and Fig. 8c. It can be seen that even the
16

very coarse mesh with the enhanced formulation results in good agreement.
Further comparison of the stress results is shown in Fig. 10, where it can been
seen that along the traction boundary, the finer mesh of enhanced elements is able
to capture the imposed constant stress field more accurately than the very much
finer reference mesh.
The perturbation fields are based on the assumption of

(a)

(b)

(c)

Figure 7: Coarse mesh solution: a) Enhanced finite element mesh in yz plane at x = 0, b) σyy in
ref
2
yz plane, c) error calculated as (σyy
− σyy )2 /σyy

(a)

(b)

(c)

Figure 8: Finer mesh solution: a) Enhanced finite element mesh in yz plane at x = 0, b) σyy in yz
ref
2
plane, c) error calculated as (σyy
− σyy )2 /σyy

a heterogeneity in an infinite medium but the enhanced formulation still exhibits
convergence in the regions strongly influenced by the traction boundary.
5.2. L-shaped specimen
The proposed modelling strategy is also demonstrated on an example with a large
number of inclusion. A 3D L-shaped specimen with fully fixed boundary conditions
on the right surface of the right-hand arm and normal traction applied on the top
17

z

y

x

(a)

(b)

z

y

(c)

(d)

Figure 9: Reference analysis: a) discretization of entire body containing 309,406 tetrahedra, b) mesh
refinement on surface of heterogeneities, c) Reference finite element mesh in yz plane at x = 0, d)
ref
in yz plane.
σyy

(a)

(b)

Figure 10: Detailed 3D-plots of σyy stress concentrations due to the boundary effects in yz plane
at x = 0: Comparison of solution from reference analysis with a) coarse mesh C-HTS solution with
24 enhanced elements b) refined mesh C-HTS solution with 1,536 enhanced elements

18

surface of the left-hand arm is analysed, see Fig. 11. The length of the plate is 300 in
both x and y direction, the depth is 150 in z direction. The Young moduli were chosen
as E = 1 and E = 2 for matrix and inclusion respectively. Poisson’s ratio was ν = 0.1
for both material phases. The microstructure comprised 2,523 spherical inclusions
varying in size between 4 and 8 with a uniform spatial distribution (Fig. 11b). All
units are consistent.
The solution for three different mesh densities (Fig. 11c, d & e) are compared
in Fig. 12. These results are plotted on the x-y mid-plane. Fig. 12(a–d) shows
a plot of the σxy stress component for the three different meshes. These results
show that the complex stress distribution resulting from the heterogeneities can be
captured and that the solution is converging with mesh refinement. Further local
mesh refinement near corners and stress concentrations is possible, although this
was not undertaken in this case. As an estimate of the computational overhead of
the enhanced formulation for this particular problem on 16 processors, we note that
the total solution for 30, 007 C-HTS elements was 597s, whereas the problem on the
identical mesh of HTS elements for the equivalent homogeneous problem consumed
1.6s of computer time. This represents a large increase of computational time in
comparison to the homogenous problem. However, it should be noted that not all
aspects of the solution procedure have been parallelised. It should also be noted
that it was not possible to obtain the reference solution for this problem by means
of conventional FEA due to the excessive number of degrees of freedom associated
with a mesh required to resolve all of the heterogeneities. The mesh generation
itself was not possible using our currently available software and hardware facilities.
Furthermore, comparison with a PUM based solution, was also not undertaken due
to the complexity of the problem.
6. Conclusions
A new micromechanics-enhanced finite element formulation has been presented
for modelling the influence of a large number of heterogeneities in composite materials
in a computationally efficient manner. The strategy exploits closed form solutions
derived by Eshelby for ellipsoidal inclusions in order to determine the mechanical
perturbation fields as a result of the underlying heterogeneities. Approximation
functions for these perturbation fields are then incorporated into a finite element
formulation to augment those of the macroscopic fields. A significant feature of
this approach is that the finite element mesh does not explicitly resolve the heterogeneities, although the resulting solution still explicitly accounts for their presence.
In contrast with traditional homogenization approaches, this method does not rely
19

(a)

(c)

(b)

(d)

(e)

Figure 11: L-shaped specimen a) geometry, b) microstructure comprising 2,523 spherical inclusions,
c) coarse mesh comprising 1074 tetrahedral elements, d) medium mesh comprising 8772 elements
and e) fine mesh comprising 30007 elements.

on separation of scale and does not suffer from loss of information due to averaging
or localization.
The proposed technique has been implemented into a hybrid-Trefftz stress (HTS)
element formulation and it has been shown that the resulting enhanced elements (CHTS) require significantly fewer degrees of freedom to capture the detailed mechanical response compared to standard finite elements. The paper also outlines how the
proposed micromechanics approach could be used within a Partition-of-Unity (PoU)
formulation, although we conclude that this does not fully exploit the advantages of
PoU methods and that the proposed hybrid-Trefftz formulation is most appropriate.
A self-compatibility algorithm is used to determine the mutual interactions between inclusions, assuming that the eigenstrain fields are uniform within the domain
of each inclusion. It was found that even for topologies exhibiting extremely small
distances between the inclusions, this assumption is sufficient. Furthermore, it has

20

(a)

(b)

(c)

(d)

Figure 12: Solution of L-shaped specimen on x-y mid-plane. a) microstructure. Plots of σxy
resulting from b) coarse mesh, c) medium mesh and d) fine mesh.

been shown that boundary effects, that are not accounted for by the classical micromechanical solution due to the assumption of an infinite medium, can be captured
through local mesh refinement.
21

We have implemented this formulation into our FE code that is optimized for
parallel computing. Additional parallelization of the micromechanical aspects of
the formulation needs to be investigated for increased efficiency. Further research
is required in order to incorporate other improvements such as nonuniform eigenstrains [14], debonding effects [24] and inclusions of arbitrary shape by evaluating
the perturbation functions numerically [12, 13].
Acknowledgements. Funding by the Glasgow Research Partnership in Engineering
(GRPE) under project “Multi-scale modelling of fibre reinforced composites” is gratefully acknowledged. This research was partially supported by the Czech Science
Foundation through project GAČR 103/09/P490 and by the Ministry of Education,
Youth and Sports of the Czech Republic trough project MSM 684077003. All the
analyses were performed by means of YAFFEMS FE code. For more details we
refer to the code’s homepage at http://code.google.com/p/yaffems.
Appendix A. Detailed solution of perturbation displacements
The displacement perturbation field in an infinite homogeneous material due to a
uniform eigenstrain ετij applied to an ellipsoidal region Ω is provided by the following
integral equation [14, Eq. (11.30)]
u∗i =

1
[Ψjk,jki − 2νΦkk,i − 4(1 − ν)Φik,k ]
8π(1 − ν)

(A.1)

where ν denotes the Poisson’s ratio and the elliptic potentials Ψij and Φij are defined
as [14, Eq. (11.32)]
Z
Z
1
τ
0
0
τ
τ
Ψij = εij
|x − x | dx = εij ψ, and Φij = εij
dx0 = ετij φ
(A.2)
0
Ω
Ω |x − x |
The integrals φ and ψ in Eq. (A.2) are the harmonic and bi-harmonic potentials
respectively. Note that in Eq. (A.1) and thereafter, standard index notation is
employed, together with the generalised summation convection due to Mura [14].
Thus, a P
repeated index is summed according to the Einstein summation rule (e.g.
ai bij = 3i=1 ai bij ), whereas a non-repeated
upper-case index equals to the lowerP3
case equivalent (e.g. ai bi cIj = i=1 ai bi cij ). The symbol ajk,i denotes the partial
derivative of ajk with respect to the coordinate xi .
By combining Eq. (A.1) and Eq. (A.2), we obtain
u∗i =


 τ
1
εjk ψ,jki − 2νδjk ετjk φ,i − 4(1 − ν)δij ετjk φ,k
8π(1 − ν)
22

(A.3)

Similarly to Eshelby’s approach [17], the displacement perturbations are expressed
in compact form:
u∗i = Lεijk ετjk ,

Lεijk =

1
[ψ,jki − 2νδjk φ,i − 4(1 − ν)δij φ,k ]
8π(1 − ν)

(A.4)

where the third-order operator Lijk maps a transformation eigenstrain ετjk to the displacement perturbation field u∗i . It is therefore analogous to the well-known Eshelby
tensor [17], which relates a transformation eigenstrain to the strain perturbation
field.
The operator Lijk can be conveniently expressed in terms of the Ferrers-Dyson
elliptic integrals, e.g. [14, Eq. (11.36)]
Z ∞
ds
I(λ) = 2πa1 a2 a3
,
∆(s)
λ
Z ∞
ds
,
Ii (λ) = 2πa1 a2 a3
2
(ai + s)∆(s)
λ
Z ∞
ds
Iij (λ) = 2πa1 a2 a3
(A.5)
2
(ai + s)(a2j + s)∆(s)
λ
where ai stands for the i-th semi-axis of ellipsoid Ω and ∆(s) is obtained from
3
Y
∆ (s) =
(ai + s)2
2

(A.6)

i=1

The variable λ is the largest positive root of equation [14, Eq. (11.37)]
xi xi
=1
(aI + λ)2

(A.7)

Notice that λ is generally position dependent and non-zero for the points xi placed
outside the inclusion domain Ω, hence called the exterior points. Contrary, λ = 0
for interior points.
All integrals in (A.5) admit a closed-form expression in terms of the LegendreJacobi integrals of the first and second kind, defined as a fuction of an auxiliary angle
θ, [14, Eq. (12.17)]. It is worth noting that its definition via [14, Eq (11.18)]
s
a2
θ = sin−1 1 − 32
(A.8)
a1
23

is valid for interior points only and not everywhere as stated in [14]. Thus, it needs
to be replaced with a general formula:
s
a21 − a23
(A.9)
θ = sin−1
a21 + λ
available, e.g. in [25, 26]. Moreover, the following identity [14, Eq. (11.40.4)]
[xn xn Ii...jN (λ)],p = 2xp Ii...jP + Ii...j,p (λ)

(A.10)

will be repeatedly proved useful in the sequel.
It follows from Eq. (A.4) that to express operator Lεijk , we need to evaluate the
first-order derivatives of the potential φ and the third-order derivatives of ψ. To this
end, we start with expressions [14, Eq. (11.38)]

and

φ(λ) = 12 [I(λ) − xn xn IN (λ)]

(A.11)

n
o
ψ,i (λ) = 21 xi 2φ(λ) − a2I [II (λ) − xn xn IIN (λ)] = 12 xi Q(λ)

(A.12)

Employing Eq. (A.10), the first derivative of φ becomes
n
o
φ,i (λ) = 12 I,i (λ) − [xn xn IN (λ)],i = 12 [I,i (λ) − 2xi II (λ) − I,i (λ)] = −xi II (λ)
(A.13)
The third derivative of potential ψ is expressed from Eq. (A.12) as
ψ,ijk (λ) = 12 [δij Q,k (λ) + δik Q,j (λ) + xi Q,jk (λ)]

(A.14)

With the help of Eqs. (A.13) and (A.10), the term Q,j can be evaluated from


Q,j (λ) = 2φ,j (λ) − a2I [II (λ) − xn xn IIN (λ)],j = 2xj a2I IIJ (λ) − IJ (λ)
(A.15)
This provides the second derivatives of Q in the form
n 


o
Q,jk (λ) = 2 δjk a2I IIJ (λ) − IJ (λ) + xj a2I IIJ,k (λ) − IJ,k (λ)

(A.16)

After utilising the derivatives of Q(λ) and re-ordering the indices, Eq. (A.14) becomes




ψjki (λ) = xi δjk a2J IJI (λ) − II (λ) + xk δji a2J IJK (λ) − IK (λ)




+ xj δki a2J IJK (λ) − IK (λ) + xj xk a2J IJK,i (λ) − IK,i (λ)
(A.17)
24

with IIJ,k and II,j provided by [14, Eqs. (11.40.1, 11.40)]
Ii...jk,p (λ) =
λ,p =

(a2i

−2πa1 a2 a3
λ,p ,
+ λ) . . . (a2j + λ)(a2k + λ)∆(λ)

2xp (a2I + λ)2
a2P + λ xi xi

(A.18)

Now we are in a position to evaluate Lεijk in terms of the Ferrers-Dyson integrals.
Introducing Eqs. (A.13) and (A.17) into Eq. (A.4) gives




[8π(1 − ν)] Lεijk = xi δjk a2J IJI (λ) − II (λ) + (xk δji + xj δki ) a2J IJK (λ) − IK (λ)


+ xj xk a2J IJK,i (λ) − IK,i (λ) + 2νδjk xi II (λ) + 4(1 − ν)δij xk IK (λ)
(A.19)
Since λ = 0 for the points inside the inclusion, recall Eq. (A.7), all derivatives of Ii...j
vanish as well. Therefore, Eq. (A.19) yields




[8π(1 − ν)] Lε,int = xi δjk a2J IJI (λ) − II (λ) + (xk δji + xj δki ) a2J IJK (λ) − IK (λ)
+ νδjk xi II (λ) + 4(1 − ν)δij xk IK (λ)
(A.20)
and Eq. (A.4) receives its final form
Lεijk = Lε,int
ijk +



1
xj xk a2J IJK,i (λ) − IK,i (λ)
8π(1 − ν)

(A.21)

For implementation purposes, it is worth noting that Eq. (A.4) admits the Voigt
representation:
 τ 


 ε11

τ 
 ∗  



ε



L111 L122 L133 L112 L123 L113 
 u1 
 22
τ 
ε
∗
33
u
=  L211 L222 L233 L212 L223 L213 
(A.22)
2ετ12 
 2∗ 



u3
L311 L322 L333 L312 L323 L313 


2ετ23 



 τ 

2ε13
References
[1] F. Feyel, J.-L. Chaboche, FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials, Computer Methods in Applied Mechanics and Engineering 183 (3-4) (2000) 309 – 330.
25

[2] M. Geers, V. Kouznetsova, W. Brekelmans, Multi-scale computational homogenization: Trends and challenges, Journal of Computational and Applied Mathematics 234 (2010) 2175 – 2182.
[3] N. Sukumar, D. Chopp, N. Moës, T. Belytschko, Modeling holes and inclusions by level sets in the extended finite-element method, Computer Methods
in Applied Mechanics and Engineering 190 (46–47) (2001) 6183–6200.
[4] N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth
without remeshing, International Journal for Numerical Methods in Engineering
46 (1) (1999) 131–150.
[5] N. Moës, M. Cloirec, P. Cartraud, J.-F. Remacle, A computational approach
to handle complex microstructure geometries, Computer Methods in Applied
Mechanics and Engineering 192 (28–30) (2003) 3163–3177.
[6] G. Wells, L. Sluys, R. D. Borst, Simulating the propagation of displacement
discontinuities in a regularised strain-softening medium, International Journal
for Numerical Methods in Engineering 53 (5) (2002) 1235–1256.
[7] L. Kaczmarczyk, C. J. Pearce, A corotational hybrid-Trefftz stress formulation
for modelling cohesive cracks, Computer Methods in Applied Mechanics and
Engineering 198 (15-16) (2009) 1298 – 1310. doi:DOI:10.1016/j.cma.2008.
11.018.
[8] P. Kabele, Multiscale framework for modeling of fracture in high performance
fiber reinforced cementitious composites, Engineering Fracture Mechanics 74 (12) (2007) 194 – 209, fracture of Concrete Materials and Structures. doi:DOI:
10.1016/j.engfracmech.2006.01.020.
[9] P. Nicolaou, S. Semiatin, A hybrid micromechanical-macroscopic model for the
analysis of the tensile behavior of cavitating materials, Metallurgical and Materials Transactions A 35 (13) (2004) 1141–1149.
[10] A. Fritsch, C. Hellmich, L. Dormieux, Ductile sliding between mineral crystals
followed by rupture of collagen crosslinks: Experimentally supported micromechanical explanation of bone strength, Journal of Theoretical Biology 260 (2)
(2009) 230 – 252. doi:DOI:10.1016/j.jtbi.2009.05.021.
[11] Z. Sharif-Khodaei, J. Zeman, Microstructure-based modeling of elastic functionally graded materials: One dimensional case, Journal of Mechanics of Materials
and Structures 3 (2008) 1773–1796. arXiv:0802.0511.
26

[12] V. Maz’ya, G. Schmidt, Approximate approximations, American Mathematical
Society, Providence, RI, 2007.
[13] J. Novák, Calculation of elastic stresses and strains inside a medium with multiple isolated inclusions, in: M. Papadrakakis, B. Topping (Eds.), Proceedings of
the Sixth International Conference on Engineering Computational Technology,
Stirlingshire, UK, 2008, p. 16 pp, paper 127. doi:10.4203/ccp.89.127.
[14] T. Mura, Micromechanics of Defects in Solids., Martinus Nijhoff Publishers, P.
O. Box 163, 3300 AD Dordrecht, The Netherlands, 1987. 587.
[15] J. Vorel, M. Šejnoha, Evaluation of homogenized thermal conductivities of imperfect carbon-carbon textile composites using the mori-tanaka method, Structural Engineering and Mechanics 33 (4) (2009) 429–446.
[16] B. Pichler, S. Scheiner, C. Hellmich, From micron-sized needle-shaped hydrates
to meter-sized shotcrete tunnel shells: micromechanical upscaling of stiffness
and strength of hydrating shotcrete, Acta Geotechnica 3 (4) (2008) 273–294.
[17] J. D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion,
and related problems, Proceedings of the Royal Society of London. Series A,
Mathematical and Physical Sciences 241 (1226) (1957) 376–396.
[18] B. Pichler, C. Hellmich, et al., Estimation of influence tensors for eigenstressed
multiphase elastic media with nonaligned inclusion phases of arbitrary ellipsoidal
shape, Journal of Engineering Mechanics 136 (2010) 1043–1053.
[19] J. Melenk, I. Babuška, The partition of unity finite element method: Basic theory and applications, Computer Methods in Applied Mechanics and Engineering
139 (1–4) (1996) 289–314.
[20] I. Babuška, J. Melenk, The partition of unity method, International Journal for
Numerical Methods in Engineering 40 (4) (1997) 727–758.
[21] G. Wells, Discontinuous modelling of strain localisation and failure, Ph.D. thesis,
Delft University of Technology (2001).
[22] J. Teixeira de Freitas, Formulation of elastostatic hybrid-trefftz stress elements,
Computer Methods in Applied Mechanics and Engineering 153 (1998) 127 –
151.

27

[23] Z. Bittnar, J. Šejnoha, Numerical methods in structural mechanics, American
Society of Civil Engineers, 1996.
[24] J. Ju, Y. Ko, Micromechanical Elastoplastic Damage Modeling of Progressive
Interfacial Arc Debonding for Fiber Reinforced Composites, International Journal of Damage Mechanics 17 (4) (2008) 307.
[25] J. D. Eshelby, The elastic field outside an ellipsoidal inclusion, Proceedings of
the Royal Society of London. Series A, Mathematical and Physical Sciences
(1959) 561–569.
[26] M. Rahman, On the Newtonian potentials of heterogeneous ellipsoids and elliptical discs, Proceedings: Mathematical, Physical and Engineering Sciences
457 (2013) (2001) 2227–2250.

28

