Jtam.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

47, 3, pp. 537-551, Warsaw 2009

APPLICATION OF THE MULTISCALE FEM TO THE

MODELING OF NONLINEAR MULTIPHASE MATERIALS

Sandra Ilic

Klaus Hackl

Ruhr-University of Bochum, Institute of Mechanics, Germany

e-mail: sandra.ilic@rub.de; klaus.hackl@rub.de

This contribution is concernedwith themodeling of compositematerials
and, particularly, with the application of the multiscale finite element
method for that purpose. The method is a result of combining homoge-
nization theory with the finite elementmethod and is based on the idea
of splitting the simulation of a heterogeneous body into two tasks: the
first one is the modeling of the actual body and the second one the mo-
deling of the representative volume element, the material sample whose
analysis replaces themissing effective constitutive law.The connectionof
these two simulation levels is achieved by introducing the Hill macroho-
mogeneity condition which requires the equality of the macropower and
the volumeaverageof themicropower.Themethodhas the advantage to
be applicable for simulation of materials with very different microstruc-
ture types. This is illustrated by the examples concerned with effective
behavior of two- and three-phase composite materials.

Key words: multiscale FEM, composite materials, homogenization

1. Introduction

Themodeling of compositematerials has already a long history. It starts with
the well known solutions of Voigt (1889) andReuss (1929) which are even to-
day often used because of their simplicity. These solutions are pretty intuitive:
Voigt’s approach states that the effective material tensor is just a volume ave-
rage of the elasticity tensor of the original composite while the Reuss solution
claims the same but for the compliance tensor. A further fundamental contri-
bution to this field has been made by Hill (1952, 1963, 1972). Among others,
he introduced the idea of energy bounds and shown that the solutions of Voigt
and Reuss can be obtained from the upper and lower energy bounds for the



538 S. Ilic, K. Hackl

assumedhomogeneous state of strain and stress.Hill’s work openedmanynew
possibilities for the modeling of composite materials, with a tendency to find
the narrowest possible energy bounds.This ideawas also followed in theworks
of Hashin and Shtrikman (1962a,b, 1963), proposing energy bounds based on
the consideration of a comparison homogeneous body. Such an approach was
also suggested in contributions considering the modeling of nonlinear compo-
site materials (Castañeda, 1991, 1992; Talbot andWillis, 1985).

Finally, the complexity of analytical solutions as well as the rush develop-
ment of computer technology gave rise to the strong development of numerical
methods appropriate for themodeling of composites. Somemethodswhich are
currently used most frequently are those of domain decomposition, adaptive
hierarchical modeling and the multiscale FEM.

The essential idea of the domain decomposition method (Zohdi andWrig-
gers, 1999; Zohdi et al., 2001) is to choose a homogeneous substitutematerial
used to calculate the approximative deformation of the whole body, the so-
called regularized problem. Furthermore, the body is decomposed into seve-
ral subdomains where the real material structure is considered. The bounda-
ry conditions at this level are the deformations calculated in the regularized
step. The final result is obtained as a superposition of the regularized solution
(step I) and local perturbations of the regularized solution for each particular
subdomain (step II). From the aspect of numerical costs, the method is very
convenient as a big system of equations is split into several small ones.

The method of hierarchical adaptivity (Oden and Zohdi, 1997; Zohdi et
al., 1993) also starts from the idea of regularization. Here again a substitute
material, not necessarily homogeneous, is chosen. The domain is decomposed
in several partitions, in each of them the error due to the regularization is
calculated and compared with a prescribed tolerance. For the subdomains
where the error bound is exceeded, a new finer subpartitioning is carried out
and more precise data on the material microstructure are introduced. The
process is iteratively continueduntil the error remains in thepermitteddomain
for all subpartitions.

A basic disadvantage of both numerical methodsmentioned is that in spi-
te of the even high-level partitioning, the influences of the very fine micro-
structure cannot be caught. This disadvantage is eliminated in the approaches
belonging to the group of homogenization methods which has been developed
for the so-called statistically uniform materials. This kind of materials has a
crucial property that a corresponding representative volume element (RVE)
can bedefined for them.AnRVE represents the smallestmaterial sample con-
taining enough information about the material microstructure such that the



Application of the multiscale FEM... 539

analysis of it yields the same effective material parameters independent from
its shape. The notion of the RVE is rather important as the application of the
homogenization concept assumes that the ratio of the characteristic size of
the RVE and of the simulated body has to tend to zero. Themultiscale finite
element method (Hackl and Ilic, 2005; Ilic, 2008; Ilic and Hackl, 2004, 2005,
2006; Miehe et al., 2002a,b,c; Schröder, 2000), which is the main topic of this
contribution, also belongs to the group of homogenization methods, and its
concept and application within the paper will be presented as follows.

Thebasic definitions and theorems of thehomogenization theoryare expla-
ined in Section 2. The combination of this method with the finite element
method (FEM) is demonstrated on a particular example concerned with the
nonlinear elastic materials with the potential proposed by Simo, Taylor and
Pister (Chap. 3). In Section 4, we consider two numerical examples. The first
example investigates the effective material parameters of a material with cir-
cular microinclusions periodically ordered in the matrix material. The second
example intends to simulate materials with the transit region around the inc-
lusion representing amixture of thematrix and inclusionmaterial. The paper
ends with several conclusions and a brief outlook.

2. Homogenization theory

As indicated in the introduction, the concept of homogenization requires two
problems to be simultaneously studied: themodeling of themacroscopic body
and the analysis of an RVE corresponding to its material structure. Of cour-
se, the quantities related to different problems are dependent on each other.
Moreover, according to Hill (1972), the macroscopic quantities are dependent
only on the microquantities occurring at the boundary of the RVE. This can
be easily explained if the RVE is understood as a small piece cut from the
macroscopic body. Then the influences of the body on the cut part are just
those occurring at the boundary of the RVE and not those inside it. Another
motivation is provided by the fact that laboratory investigations are carried
out on the boundary of the sample, being an RVE in the present case.

As we are particularly interested in the theory of finite deformations, the
definitions of the macroscopic first Piola-Kirchhoff stress tensor P and defor-
mation gradient F will be elucidated especially in the following. Note that
hereafter the overbar symbol will be used to denote the quantities related to
the macrolevel while those without overbar symbol relate to the microlevel.
For the definition of themacroscopic quantities, an arbitraryRVE Bwith the



540 S. Ilic, K. Hackl

boundary ∂B is considered (Fig.1). Such an RVE may also contain one or
more voids and their boundaries are denoted by G in Fig.1. It is also possible
to study the singular surfaces inside theRVE, but in that casemostly the zero
jump of the deformations is assumed which leads to the same definitions as
they are presented here (Schröder, 2000).

Fig. 1. An arbitrary RVE with a cavity

In order to define the first Piola-Kirchhoffmacrostress tensor, the analysis
is started by observing that the volume average of itsmicroscopic counterpart
over the volume V of the RVE B is

〈P〉=
1

V

∫

B

P dV (2.1)

Here the volume average is denoted by the symbol 〈·〉. Bearing in mind that
P = 0 on the boundary parts without load (such as the boundaries of the
intrinsic voids G) and that GradX = I, Eq. (2.1) can be transformed into
the following boundary integral

1

V

∫

B

P dV =
1

V

∫

∂B

(PN)⊗X dA=
1

V

∫

∂B

T ⊗X dA (2.2)

Here T represents the tractions, i.e. it holds T =PN, X is thepositionvector
and N the normal vector on the surface of the boundary, all of them in the
reference configuration. The notation ”Grad” represents differentiation with
respect to the coordinates in the reference coordinate system. As the volume
average of the microscopic stress tensor can be rewritten as an expression
depending only on quantities occurring at the boundary of the RVE (2.2), the
volume average can be directly assumed as the definition of the macroscopic
Piola-Kirchhoff stress tensor

P= 〈P〉=
1

V

∫

B

P dV (2.3)



Application of the multiscale FEM... 541

A similar procedure repeated for themicrodeformation gradient yields its vo-
lume average in the form

1

V

∫

B

F dV =
1

V

(∫

∂B

x⊗N dA+

∫

G

x⊗N dA
)

(2.4)

where x, in contrast to X, denotes the coordinates in the current configura-
tion. Obviously, the volume average cannot be expressed dependent only on a
single boundary integral over the external boundary of the RVE as the defor-
mations at the internal boundaries are, in general, different from zero. For this
reason we take the first term at the right hand side of (2.4) as the definition
of the macrodeformation tensor

F :=
1

V

∫

∂B

x⊗N dA=
1

V

(∫

B

F dV −

∫

G

x⊗N dA
)

(2.5)

Obviously, in the case of the continuous composite materials, the macrode-
formation tensor can be equalized with the volume average of its microscopic
counterpart.
Given the definitions of the macroscopic quantities, the last component

that needs to be specified are the boundary conditions at the microlevel. To
this end, an additional requirement concerning the connection of macro and
micro power is introduced. This is known as theHill condition and states that
macropower must be equal to the volume average of the micropower (Hill,
1963, 1972; Huet, 1990)

P : Ḟ=
1

V

∫

B

P : Ḟ dV (2.6)

An alternative form of the Hill condition (2.6) can be derived bymeans of the
Gauss theorem

1

V

∫

∂B

(T −PN)(ẋ− ḞX) dA=0 (2.7)

This expression is especially convenient as it allows one to directly define two
types of boundary conditions satisfying (2.6) and (2.7), respectively. These
conditions are known as the static boundary conditions

T =PN on ∂B (2.8)

and the kinematic boundary conditions

x=FX on ∂B (2.9)



542 S. Ilic, K. Hackl

andbecause of their simplicity they are regularly applied in laboratory investi-
gations. A further possible solution to (2.7) is derived by using the assumption
that the real deformation at the microscale is equal to the macrodeformation
superimposed with microfluctuations w̃

x=FX+ w̃ (2.10)

Substitution of this assumption into (2.7) results in the conclusion that this
condition is also satisfied in the case that themicrofluctuations w̃ are periodic
and the tractions T antiperiodic at the periodic boundary of the RVE

w̃
+ = w̃− T+ =−T− (2.11)

Here the superscripts ”+” and ”−” distinguish the quantities related to the
boundary parts whose normal vectors are oppositely oriented: N+ =−N−.
The investigations show that this type of boundary conditions leads to very
good results even if the microstructure is not perfectly periodic (Zohdi and
Wriggers, 2005). An important advantage of using assumption (2.10) is the
simple relation of the microdeformation tensor: it is just a superposition of
the macrodeformation and the microfluctuation gradient

F= Gradx=F+ Gradw̃=F+ F̃ (2.12)

3. FEM implementation

Let us summarize the data of the boundary value problems (BVPs) that are
to be solved within the scope of homogenization theory. The first BVP is
concerned with the simulation of themacroscopic body. Here the geometry as
well as the loads are defined. Thematerial properties and thematerial tensor
or directly the stress tensors are obtained from the microscale calculations.
The secondBVP considers the behavior of theRVE for themacrodeformation
obtained as the input from the macroscale, and the periodic (alternatively
static or kinematic) boundary conditions derived from theHill condition have
tobe satisfied.Thegeometry including thematerial characteristics at this level
is known. Obviously, after introducing the Hill macrohomogeneity condition,
both BVPs are complete and any standard method can be used for their
solution. Our particular choice considers a combination of homogenization
theory with FEM, known as themultiscale or FE2 method. The advantage of
this approach is that it is easily adaptable for simulating materials with very



Application of the multiscale FEM... 543

different microstructures. This applies to the geometry of the RVE as well as
to the constitutive laws of the component materials.
Within this contribution we will consider a case where the material beha-

vior at both scales is determined by the following potential in the three-field
description (Simo and Hughes, 1997)

Π(u,Θ,p)=

∫

V

[Ψvol(Θ)+Ψdev(C
∗(u))+p(J(u)−Θ)] dV +Πext (3.1)

Its advantage in comparisonwith the standard formulation is the applicability
to cases of nearly incompressible materials. This is achieved by introducing
the secondary variables Θ and p, permitting the potential to be split into a
deviatoric and a volumetric part. The new variables Θ and p represent the
volumetric change and hydrostatic pressure respectively, C∗ is the isochoric
rightCauchyGreen deformation tensor defined as C∗ = J−

2

3C, and J =detF
is the Jacobian. For a complete specification of the material behavior, the
Helmholtz free energy corresponding to the Neo-Hook material is chosen

Ψ =Ψdev +Ψvol =
1

2
µ(trC−3)+K(J lnJ−J+1) (3.2)

where K and µ are the bulk and shear moduli, respectively.
Since the potential is fixed,we can now focus on theFEM implementation.

To this end the second variation of (3.1) is needed, and using the standard
tools of variational calculus (Boonet andWood, 2000; SimoandHughes, 1997)
it can be shown that the necessary expression has the form
∫

V

Gradδu : [Grad∆u · (Sdev +Svol)] dV +

+

∫

V

(Grad
⊤
δu ·F) : (Adev +Avol) : (F

⊤
· Grad∆u) dV +

+

∫

V

(Grad
⊤
δu ·F) :JC

−1
dV
( 1
V

∂2Ψvol

∂Θ
2

)∫

V

JC
−1
: (F
⊤
· Grad∆u)dV +

+∆δuΠ
ext
=−δuΠ

res
(3.3)

u=u0, δu=0, ∆u=0 on ∂Bu

t= t0 on ∂Bt

Recall that the second variation in this concept represents the linearized form
of the condition that the actual state of the deformationminimizes the poten-
tial. In (3.3), the followingnotation isused: ∆denotes the increment, Sdev and



544 S. Ilic, K. Hackl

Svol are thedeviatoric andvolumetricpartsof the secondPiola-Kirchhoff stress
tensor, defined in a standard way as

Sdev =2
∂Ψdev(C)

∂C
Svol =2

∂Ψvol(C)

∂C
(3.4)

and Adev and Avol are the corresponding elasticity tensors

Adev =2
∂Sdev

∂C
=4
∂2Ψdev(C)

∂C
2

(3.5)

Avol =2
∂Svol

∂C
=4
∂2Ψvol(C)

∂C
2

The derivative of the pressure has to be calculated according to

∂2Ψvol(C)

∂Θ
2
=
∂p

∂J
=
2

3
J
∂p

∂C
:C
∗

(3.6)

As expression (3.3) exactly corresponds to the BVP at the macroscale, the
overbar symbol is directly introduced. The particularity of the macroscale
problem is that the quantities dependent on Ψ(C), such as the stresses and
their derivatives, cannotbecalculateddirectlybutusingthemicroscale results,

expression (2.3) and relation S=F
−1
P. For the calculation of Adev,Avol and

∂p/∂C, the numerical derivatives must be used in addition.
The same material description (3.1) is used at the microlevel so that for

the macrodeformation imposed from the macrolevel (in the calculation it is
treated as the value from the previous step of NewtonRaphson iteration) and
neglecting the external load, the microlevel BVP becomes

∫

V

Gradδw̃ : [Grad∆w̃ · (Sdev +Svol)] dV +

+

∫

V

(Grad⊤δw̃ ·F) : (Adev +Avol) : (F
⊤ · Grad∆w̃) dV +

+

∫

V

(Grad⊤δw̃ ·F) : JC−1 dV
( 1
V

∂2Ψvol
∂Θ2

)∫

V

JC−1 : (F⊤·Grad∆w̃) dV+

+δ
w̃
Πresmac =−δw̃Π

res (3.7)

w̃
+ = w̃−, δw̃=0, ∆w̃=0 on ∂B

δ
w̃
Πresmac = f(F)



Application of the multiscale FEM... 545

In contrast to (3.3), this expression depends on the microfluctuations w̃.
Periodic boundary conditions as a consequence of the Hill macrohomogene-
ity conditions are assumed. Of course static or kinematic boundary condi-
tions are possible alternatives but due to the periodic microstructure of the
materials, which will be simulated here, periodic boundary conditions are
chosen.
The further procedure includes the standard steps of FE implementation

(Hughes, 1987; Zienkiewicz and Taylor, 2000). To this end, the displacement
function and its derivatives are approximated bymeans of the shape functions
and nodal values

u=Nûe

(3.8)
Gradu=Bûe B= GradN

e

and analogous expressions hold also for the increment and variation of displa-
cements

Grad∆u=B∆ûe Gradδu=Bδûe (3.9)

Here N is amatrix containing the shape functionsand B is amatrix containing
the derivatives of the shape functions with respect to the physical coordinates
X. A hat symbol denotes the nodal values and e indicates that all the DOFs
of an element are considered. The implementation of (3.8) and (3.9) into (3.7)
yields the element stiffness matrix K

e

K
e
=

∫

V
e

G
⊤
(Sdev +Svol)G dV +

∫

V
e

(B
⊤
F) : (Adev +Avol) : (F

⊤
B) dV +

(3.10)

+

∫

V
e

(B
⊤
F) : JC

−1
dV
( 1
V

∂2Ψvol

∂Θ
2

)∫

V
e

JC
−1
: (F
⊤
B) dV

with matrix G being defined in the following way
∫

V
e

Gradδu : [Grad∆u · (Sdev +Svol)] dV =

(3.11)

+δû

∫

V
e

[G
⊤
(Sdev +Svol)G] dV∆û

Of course, an approximation of the fluctuations w̃ at the microlevel is also
required, but this procedure is analogous to steps (3.8)-(3.11) and yields the



546 S. Ilic, K. Hackl

same results, so that it will not be repeated. The approximation formacrode-
formations (3.8)1 and its microscopic counterpart can be chosen in a different
manner dependent on the type of problems to be solved. In the current si-
mulations, the P0Q1 element with bilinear shape functions (Hughes, 1987;
Zienkiewicz and Taylor, 2000) is taken at both levels.

4. Numerical examples

Our first example examines materials whose microstructure is determined by
a square RVE containing a circular inclusion (Fig.2). The square has unit
side length 2d = 1mm and the inclusion of radius r = 0.125mm is placed
centrally. Thematerial parameters of thematrixmaterial are fixed: theYoung
modulus amounts to 1000N/mm2 and the Poisson ratio is 0.3. The Young
modulus of the inclusions is varied in different tests, and it takes values in
the range of 0-1000N/mm2. The Poisson ratio of the inclusion is constant:
ν =0.3. The plane strain state is assumed.

Fig. 2. Microlevel: RVE of two-phasematerial with circular inclusion

Themicrostructure being defined, a tension test of a square plate is consi-
dered at themacroscale. The plate has the dimensions 40×40cm, but due to
symmetry consideration of one fourth of the plate is satisfactory (Fig.3).Here,
in addition, the vertical displacements at the lower boundary and horizontal
displacements at the left boundarymust be constrained. Themesh consists of
10×10 elements with four Gauss points, and to each of these points an RVE
is associated. The load chosen is p=1kN/cm.

The final results of the tension tests are the displacements of the uncon-
strained boundaries, further used for calculation of the effective material pa-
rameters presented in Fig.4. Here both diagrams show a smooth increase of



Application of the multiscale FEM... 547

Fig. 3. Macrolevel: Tension test of a plate

the effective material parameters with respect to a change of Young’s modu-
lus of the inclusion. The effective Youngmodulus takes values in the range of
886.7-1000N/mm2 and Poisson’s ratio 0.29-0.3. Notice that the values obta-
ined for the zero Young’s modulus of inclusion (EI =0) correspond tomicro-
porous material while the other limit case (EI = 1000N/mm

2) corresponds
to the homogeneous material as the matrix and the inclusion have the same
properties.

Fig. 4. Change of Young’s modulus and Poisson’s ratio with respect to Young’s
modulus of the inclusion. Two-phasematerial is considered

The second example studies a more complex geometry. The RVE is still
a unit square with a circular inclusion of radius 0.125mm, but a third phase
is introduced in addition (Fig.5). This phase forms a ring with the internal
radius r=0.125mm and external radius r=0.25mm. The intention here is
to simulate a belt consisting of amixture of thematrix and inclusionmaterial.
Accordingly, the Young modulus of the third phase takes values in the range
of 300-1000N/mm2, where the first value is typical for the inclusion and the



548 S. Ilic, K. Hackl

second one for the matrix material. The former limit case corresponds to the
two-phase material with an inclusion of radius 0.25mm and the later limit
case to the two-phase material with the inclusion of radius 0.125mm (as in
the previous example). Themacroscale problem remains unchanged: a tension
test of the plate results in displacements used for calculation of the effective
material parameters (Fig.3). The values calculated for the new type of micro-
structure are presented in Fig.6. Here, again, both diagrams show a smooth,
gradual increase of the effective parameters.

Fig. 5. Microlevel: RVE of a three-phase material with a circular inclusion and
mixture belt

Fig. 6. Change of Young’s modulus and Poisson’s ratio with respect to Young’s
modulus of the ring. Three-phasematerial is considered

5. Conclusions

Thepresent contribution is focused on an approach to themodeling of compo-
sites known as the multiscale FEM. In order to demonstrate the concept and
implementation of thismethod,we chose a particularmaterial with the poten-
tial in three-field description and the free Helmholtz energy corresponding to



Application of the multiscale FEM... 549

the Neo-Hookmaterial. For such a combination, two illustrative examples are
considered. The first shows the change of the effective material parameters in
the case of a material with periodic microstructure with a circular inclusion
whereby thematerial parameters of the inclusion are varied.The second exam-
ple also considers the periodicmicrostructure with a circular inclusion, but in
addition a narrow belt around this inclusion is introduced. In this test group,
Young’s modulus of the belt varies taking values from the range limited by
the values characteristic for the inclusion and matrix material. The resulting
diagrams for both examples show a gradual increase of the effective material
parameters.

The simulations presented endorse that the multiscale FEM can be relia-
bly applied to the modeling of composite materials with complex microstruc-
ture. This is a strong motivation for further concentrating on the modeling
of practice-oriented macro- andmicroscale problems such as estimation of the
effective properties of damaged materials or simulation of composites consi-
sting of solid and fluid fractions, which is typical for biomechanics. The 3D
simulations as well as the implementation of the static and kinematic bo-
undary conditions at the microlevel are some other challenging topics to be
investigated.

References

1. Boonet J., Wood R.D., 2000, Nonlinear Continuum Mechanics for Finite
Element Analysis, Cambridge University Press

2. Castañeda P.P., 1991, The effective mechanical properties of nonlinear iso-
tropic composites, J. Mech. Phys. Solids., 39, 1, 45-71

3. Castañeda P.P., 1992, New variational principles in plasticity and their ap-
plication to composite materials, J. Mech. Phys. Solids., 40, 8, 1757-1788

4. Hackl K. Ilic S., 2005, Solution-precipitation creep – continuummechanical
formulation andmicromechanical modelling,Arch. Appl. Mech., 74, 773-779

5. Hashin Z., Shtrikman S., 1962a,A variational approach to the theory of the
elastic behaviour of polycrystals, J. Mech. Phys. Solids, 10, 343-352

6. Hashin Z., Shtrikman S., 1962b,On some variational principles in anisotro-
pic and nonhomogeneous elasticity, J. Mech. Phys. Solids, 10, 335-342

7. Hashin Z., Shtrikman S., 1963, A variational approach to the theory of the
elastic behaviour of multiphase materials, J. Mech. Phys. Solids, 11, 127-140



550 S. Ilic, K. Hackl

8. Hill R., 1952, The elastic behaviour of a crystalline aggregate, Proc. Phys.
Soc., London, A, 65, 349-354

9. Hill R., 1963, Elastic properties of reinforced solids: some theoretical princi-
ples, J. Mech. Phys. Solids, 11, 357-372

10. HillR., 1972,Onconstitutivemacro-variables forheterogeneoussolidsatfinite
strain,Proc. R. Soc. Lond., A, 326, 131-147

11. Huet C., 1990, Application of variational concepts to size effects in elastic
heterogeneous bodies, J. Mech. Phys. Solids, 38, 6, 813-841

12. Hughes T.J.R., 1987,The Finite Element Method, Prentice Hall, Englewood
Cliffs, NJ

13. Ilic S., 2008,Application of theMultiscale FEM to theModeling of Composite
Materials, Ph.D. Thesis, Ruhr-University Bochum

14. Ilic S., Hackl K., 2004, Homogenisation of random composites via the mul-
tiscale finite-element method,PAMM, 4, 326-327

15. Ilic S., Hackl K., 2005, Solution-precipitation creep –micromechanicalmo-
delling and numerical results,PAMM, 5, 277-278

16. IlicS.,HacklK., 2006,MultiscaleFEMinmodellingof solution-precipitation
creep,PAMM, 6, 483-484

17. Miehe C., Schotte J., Lambrecht M., 2002a, Homogenisation of inelastic
solid materials at finite strains based on incremental minimization principles,
J. Mech. Phys. Solids, 50, 2123-2167

18. Miehe C., Schröder J., Bayreuther C., 2002b, On the homogenisation
analysis of composite materials based on discretized fluctuations on themicro-
structure,Acta Mechanica, 155, 1-16

19. Miehe C., Schröder J., Becker M., 2002c, Computational homogeniza-
tion analysis in finite elasticity: material and structural instabilities on the
micro- and macro-scales of periodic composites and their interaction, Comp.
Met. Appl. Mech. Eng., 191, 4971-5005

20. Oden J.T., Zohdi T.I., 1997, Analysis and adaptive modeling of highly he-
terogeneous elastic structures,Comp. Met. Appl. Mech. Eng., 148, 367-391

21. Reuss A., 1929, Calculation of flow limits ofmixed crystals on the basis of the
plasticity of mono-crystals,Z. Angew. Math. Mech., 9, 49-58

22. Schröder J., 2000, Homogenisierungsmethoden der nichtlinearen Kontinu-
umsmechanik unter Beachtung von Stabilitäts Problemen, Habilitationsshrift,
Universität Stuttgart

23. Simo J.C., HughesT.J.R., 1997,Computational Inelasticity, SpringerVerlag

24. TalbotD.R.S.,Willis J.R., 1985,Variational principles for inhomogeneous
non-linear media, IMA-Journal of Applied Mathematics, 35, 39-54



Application of the multiscale FEM... 551

25. Voigt W., 1889, Über die Beziehung zwischen den beiden Elasti-
citätskonstanten isotroper Körper,Ann. Phys., Leipzig, 38, 3, 573-587

26. Zienkiewicz O.C., Taylor R.L., 2000, The Finite Element Method,
Butterworth-Heinemann

27. Zohdi T.I., OdenJ.T., Rodin G.J., 1993, Hierarchical modeling of hetero-
geneous bodies,Comp. Met. Appl. Mech. Eng., 138, 273-298

28. Zohdi T.I.,Wriggers P., 1999, A domain decompositionmethod for bodies
with heterogeneousmicrostructure based on thematerial regularization, Int. J.
Sol. Struct., 36, 2507-2525

29. Zohdi T.I.,Wriggers P., 2005, Introduction to Computational Micromecha-
nics, SpringerSeries in:LectureNotes inApplied andComputationalMechanics,
20

30. ZohdiT.I.,WriggersP.,HuetC., 2001,Amethod of substructuring large-
scale computationalmicromechanical problems,Comp.Met. Appl. Mech. Eng.,
190, 13, 5639-5656

Zastosowanie wieloskalowej MES do modelowania nieliniowych

wieloskładnikowych materiałów

Streszczenie

Pracadotyczymodelowaniamateriałówkompozytowych, aw szczególności zasto-
sowania wieloskalowej metody elementów skończonych do tego celu. Metoda ta jest
kombinacjąMES oraz teorii homogenizacji i opiera się na podziale zadania symulacji
niejednorodnego ciała na dwa poziomy: modelowania właściwego tego ciała i mode-
lowania reprezentatywnego elementumateriału próbki, którego analizama uzupełnić
brakujące efektywne równanie konstytutywne materiału. Połączenie obydwu pozio-
mów symulacji jest osiągalne po wprowadzeniu warunku homogenizacji Hilla, który
narzuca równośćmocy w skali makro z mocąmikroskalową uśrednioną objętościowo.
Zaletąmetody jest jej aplikacyjność do analizymateriałów o bardzo różnorodnej mi-
krostrukturze.Zilustrowanotowpracynaprzykładziebadańefektywnegozachowania
dwu- i trójskładnikowych kompozytów.

Manuscript received February 23, 2009; accepted for print March 13, 2009