Determining the size of RVE for nonlinear

random composites in an incremental

computational homogenization framework

Abstract

In this paper, we address the issue of determining the size of a representative

volume element (RVE) in the case of nonlinear random composites with either
elastoplastic or elasto-viscoplastic phases. In such a case, the general form of the
effective constitutive behavior is not known in advance and the response must be
evaluated either by direct numerical computations on the RVE, either by an ap-
propriate approximation scheme. Previous methodologies for determining the size
of RVE usually rely on analyzing the convergence of the RVE response computed
numerically with respect to its size. In the present work, we analyze the convergence
of parameters related to an incremental homogenization scheme, with respect to (i)
the size of the RVE and (ii) to statistical convergence related to microstructure real-
izations. For that purpose, we combine an incremental homogenization method with
a statistical convergence analysis of parameters related to the matrix phase only.
The advantage is that the range of parameters to be identified is much narrower
than for a general empirical constitutive law. Once identified and the convergence
analysis performed with respect to both size of RVE and statistical realizations,
the macroscopic constitutive law can be readily used for structure calculations.
We illustrate the methodology by analyzing two-dimensional microstructures with
randomly distributed cylindrical elastic rigid fibers, embedded in a elastoplastic or
elasto-viscoplastic matrix. For these materials, the existence of an RVE is demon-
strated for sizes of RVE corresponding to 17 − 18 and 14 − 15 times the diameter
of the inclusions, respectively.

Key words: Computational homogenization, Size of RVE, nonlinear composites,

Incremental homogenization.

1 Introduction

A fundamental concept in homogenization of composites is that of Represen-

tative Volume Element (RVE), determining the domain associated with the
microstructure to evaluate the effective behavior of the heterogeneous materi-
als. If this notion is nowadays well defined for random linear composites, many

Preprint submitted to Elsevier Science December 16, 2014

issues remain when considering nonlinear heterogeneous materials. A review of
the different definitions of an RVE can be found in [9] for linear and nonlinear
materials. Several statistical processes have been proposed, including, among
others, the procedures of Kanit et al. [15], Gitman et al. [9], or Pelissou et al.
[20].

A first class of methods for studying the size of RVE includes approaches
based on analytical homogenization schemes, mainly restricted to linear cases.
These techniques (see e.g. [17,21,33]), have been used to consider spherical
or spheroidal-shaped inclusions and have been useful in some situations to
determine the size of the RVE with respect to the inclusions characteristic
size.

A second class of approaches, based on numerical methods such as the finite

element method (FEM) (see, e.g., among many others, [2,10,22,25,31]), uses
computations on a unit cell and allows determining the size of the RVE via
statistical analyses relying on numerical computations. These techniques have
been mainly applied in the linear case, and a few recent studies involve non-
linear heterogeneous materials. For linear composites, determining the size
of the RVE can be performed by analyzing the statistical convergence of ef-
fective material parameters with respect to the size of the unit cell. In [15],
Kanit et al. studied the linear thermal and elastic properties of random 3D
polycristalline microstructures. In [23], Ostoja-Starzewski et al. investigated
random polycristal microstructures made-up of cubic single crystals. Other
examples in elasticity can be found in [7,12,18,24,34]. Applications to cortical
bone, molecular dynamics models of polymer or porous media have been stud-
ied in [5,11,32]. In [28], new criteria to determine the size of RVE with random
elastic matrix have been proposed as well as estimates for RVE sizes. In [30], a
stochastic homogenization theory has been introduced for random anisotropic
elastic composites which cannot be described in terms of their constituents
and for which the standard methods cannot be applied, like cortical bones
or biological membranes. In [26], a method using the concept of periodiza-
tion of random media was used to estimate the effective properties of random
composites using small volumes.

For nonlinear compos, most of the proposed methodologies are based on an-
alyzing the convergence of the effective response (e.g. the effective stress),
computed numerically at one point of the loading curve, with respect to the
size [3,13,27]. More recent studies analyze the convergence of identified pa-
rameters related to an empirical macroscopic model with respect to the unit
cell size [16,20].

In the present paper, a new methodology to estimate the size of the RVE for
nonlinear composites is proposed. The technique combines statistical analysis
based on the convergence of identified macroscopic parameters and homog-

2
enization. An incremental homogenization technique [4] is employed, which
provides the effective constitutive behavior of the material. Then, the conver-
gence of parameters associated with the homogenization process is analyzed
with respect to the RVE size through an appropriate metric. The advantages
of the present framework are twofold: the size of the RVE for nonlinear ma-
terials like elastoplastic or elasto-viscoplastic ones can be estimated; (b) the
converged parameters can be used to efficiently compute the effective con-
stitutive law through the incremental homogenization scheme for structure
calculations.

The paper is organized as follows. In section 2, the incremental homogenization

method for nonlinear materials is reviewed. In section 3, the proposed statis-
tical process is presented. In section 4, the technique is applied to the analysis
of random nonlinear fiber-composites whose matrix has either a elastoplastic
or an elasto-viscoplastic behavior.

2 Review of the incremental homogenization method

Incremental homogenization schemes are extensions of the formulation pro-

posed by Hill [14] in which macroscopic stress and strains are related by a
constitutive law in the form:

σ̇(t) = C(t) : ε̇(t) (1)

where σ̇ and ε̇(t) are rates of macroscopic stress and strain, respectively and
C(t) is the tangent operator depending on both the actual strain state and on
the history of loading. For the linearized problem, it is possible to apply the
superposition principle and to compute the effective tangent operator C(t) at
each iteration, by the knowledge of the nonlinear constitutive law in each phase
and of the strain at the former iteration. Let Ω be a domain in Rd associated to
the RVE, d being the space dimension. The boundary of Ω is denoted by ∂Ω.
We assume that Ω contains two phases, the matrix, associated to a domain Ω0
and the inclusion, defined in a domain Ω1 . The interfaces between the phases
are assumed to be perfectly bounded.

For a given analytical homogenization scheme (e.g. self-consistent, Mori-Tanaka,

see e.g. [6]), given the strain increment ∆ε prescribed on the RVE at time tn ,
it is possible to evaluate the tangent modulus associated to the nonlinear con-
stitutive laws in each phase (see e.g. [4,1,29]), which are used to compute the
effective tangent modulus at time tn+1 . For example, in the scheme proposed
by Doghri et al. [4], the problem is to determine at time tn+1 the average strain
in the inclusions, which is evaluated iteratively. Let ⟨∆ε⟩Ω = ∆ε be the∫ strain
increment, denoted at time tn by ⟨∆εn ⟩Ω = ∆εn . The symbol ⟨.⟩ω = |ω| 1
ω (.)dΩ

3
is the spatial averaging over a domain ω. The matrix and inclusions are associ-
ated to the indices 0 and 1, respectively. The average strain in the matrix and
in inclusions are denoted by ⟨∆ε⟩Ω0 and ⟨∆ε⟩Ω1 , respectively. For a new in-
crement of ⟨∆εn ⟩Ω1 , the average strain in the matrix ⟨∆ε⟩Ω0 can be computed
by

∆ε = ⟨∆ε⟩Ω0 (1 − f1 ) + f1 ⟨∆ε⟩Ω1 , (2)

where f1 denotes the volume fraction of inclusion. Using the expression of

the tangent modulus computed in each phase, the Eshelby tensor E can be
evaluated (see e.g. [6]). We can then compute the concentration tensor Bϵ to
relate the average strain in each phase to the macroscopic strain by

⟨∆ε⟩Ω0 = [f1 Bϵ + (1 − f1 )I]−1 : ε, (3)

⟨∆ε⟩Ω1 = Bϵ : [f1 Bϵ + (1 − f1 )I]−1 : ε, (4)

where Bϵ is expressed by
{ [ ]}−1
Bϵ = I + E : C−1
0 : C1 − I , (5)

where C0 and C1 are the tangent moduli associated to the nonlinear constitu-
tive laws in each phase. Eventually, the effective tangent modulus is expressed
by

C = [f1 C1 : Bϵ + (1 − f1 )C0 ] : [f1 Bϵ + (1 − f1 )I]−1 . (6)

For a Mori-Tanaka scheme, the algorithm proposed by Doghri et al. in [4], is

described as follows. Let [tn , tn+1 ] a time interval. Given εn , ∆εn and internal
variables in phases at time tn , the problem consists in determining the overall
stress σ n+1 and the tangent modulus Cn+α , where n + α denotes the time
tn+α = tn + α∆t. The different steps of the algorithm are summarized as
follows:

• Initialization ⟨ε⟩Ω1 = ∆ε.

• WHILE ∥R∥ > T OL , where T OL is a numerical tolerance:
(1) Iteration k (upper index k is omitted for the sake of simplicity)
(2) Given ⟨εn ⟩Ω1 and ⟨∆ε⟩Ω1 , compute the tangent modulus C1 ≡ C1(n+1)
from the nonlinear constitutive law of the inclusion.
(3) Compute the average strain in the matrix by (2):
∆ε − f1 ⟨∆ε⟩Ω1
⟨∆ε⟩Ω0 = . (7)
1 − f1
(4) Given ⟨εn ⟩Ω0 et ⟨∆ε⟩Ω0 , compute the tangent modulus C0 ≡ C0(n+1) from
the nonlinear constitutive law of the matrix.
0 from C0 (see a justification and details in
(5) Extract the isotropic part Ciso
[4]).

4
 (6) Compute the Eshelby tensor E for Ciso 0
(7) Compute C0 and C1 at time tn+α :
Ci(n+α) = (1 − α)Ci(n) + αCi(n+1) i = 0, 1, α ∈ ]0, 1] . (8)
(8) Compute the concentration tensor B by ϵ
{ [ ]}−1
Bϵ = I + E : C−1 0(n+α) : C1(n+α) − I . (9)
(9) Check compatibility of the average strain in the inclusion by computing
the residual:
R = Bϵ : [f1 Bϵ + (1 + f1 )I]−1 : ∆ε − ⟨∆ε⟩Ω1 . (10)
(10) IF ∥R∥ ≤ T OL THEN END of iterations
(11) ELSE GO TO (1) with the new average strain in the inclusion
⟨∆ε⟩k+1
Ω1 = ⟨∆ε⟩Ω1 + ξR, ξ ∈ ]0, 1] .
k
(11)
• At convergence, compute the effective tangent modulus Cn+α and the macro-
scopic stress by:
[ ]
Cn+α = f1 C1(n+α) : Bϵ + (1 − f1 )C0(n+α) : [f1 Bϵ + (1 − f1 )I]−1 , (12)
∆σ = Cn+α : ∆ε, (13)
σ n+α = (1 − α)σ n + α∆σ. (14)

3 Proposed methodology to determine the size of the RVE

3.1 Description of the procedure

The procedure to determine the RVE’s size is described as follows.

{ }K
We define a set of volumes χV = V k , where K is determined by the
k=1
tolerance criterion defined in (18)-(22), and V 1 the minimal size of the unit
cell. The volume is given by V = L2 × H, L denoting the length of the square
domain associated with Ω(V ) and H has a unitary [Link] boundary of
Ω(V ) is denoted by ∂Ω(V )

For each volume V k , we assume that the geometry of the microstructure is

modeled by a Rm − valued random variable with probability distribution Pξ
for which its support is a subset of Rm . The random possible geometrical
parameters are represented by the random vector ξ. The independent real-
izations χξ = {ξ r }R
r=1 of random vector ξ are generated here by a uniform
distribution of inclusions centers within Ω, with rejecting trial positions in-
ducing intersection with previously introduced inclusions. An illustration of
one realization of the microstructure is provided in Fig. 1. The constitutive
law in each phase is nonlinear. Then, R realizations of the microstructure mor-
phology are generated, R being determined by the criterion defined in (23).
For each realization r, a direct finite element computation is conducted, by
prescribing a macroscopic strain load history ε(t) on the RVE through periodic

5
 Figure 1. Example of realization of the microstructure.
boundary conditions:

u(t) = ε(t)x + ũ(x, t) on ∂Ω(V k ), (15)

with ũ(x, t) a periodic fluctuation. Defining a time-stepping T = [t1 , t2 , ..., T ]

and solving the nonlinear problem for t ∈ [0, T ], the effective stress response
is computed by

r1 ∫
k
σ(V , ξ , t) = k σ(V k , ξ r , t, x)dΩ (16)
V Ω(V )
k

where σ(V k , ξ r , t, x) is the local stress in the RVE at time t, for a volume V k ,
and for the realization r. This solution constitutes the reference solution for the
set {k, r}. Then, choosing an appropriate analytical constitutive law defined
by P coefficients α = {α1 , α2 , ..., αP }, the incremental procedure described
in the previous section is conducted, where the behavior of the inclusions is
determined by fixed values of coefficients, but where the behavior of the matrix
is associated with the coefficients α. The macroscopic stress response can be
efficiently evaluated by (14) and is denoted by σ̂(V k , ξ r , t).

Then, the coefficients α(ξ r , V k ) are adjusted so as to minimize the error func-
tion:
∫ T 2
e= σ(V k , ξ r , t) − σ̂(V k , ξ r , t) dt, (17)
0

with ∥e∥2 = e : e.

For one given volume V k , the operation is repeated for R realizations, until a
convergence criterion, defined in the next section, is reached. The associated
converged parameters related to the analytical constitutive law of the matrix
in the incremental homogenization scheme are denoted by α(V k ).

Then, the volume is increased to V k+1 and the former steps are repeated.
Finally, the convergence is checked with respect to the volume k, as defined in

6
 250

200

σ̄11 (MPa)
150

100
Incremental homogenization
method with fixed coefficients
50 Identification of an empirical law
Reference solution
Proposed method
0
0 0.005 0.01 0.015
ε̄11

Figure 2. Identification of the numerical parameters α by three approaches for an

elastoplastic behavior, and volume fraction f = 0.5.

the next section. The volume associated to convergence with respect to both
sizes K and realizations R gives an estimate for the RVE size, denoted by
V = V K,R .

3.2 Convergence criteria

3.2.1 Convergence of empirical parameters for a given volume

For a given volume V k , the convergence of the P parameters αj (V k ), j =

1, ..., P is determined by a simple criterion:

2Dαj
ϵrel = √ , (18)
Mαj R

where R denotes the number of independent realizations for the volume V k

fixed, for r. Let χjα be the vector collecting the R independent realizations of
the parameter αj :
{ }
χjα = αj1 (V k ), αj2 (V k ), ..., αjR (V k ) . (19)

In (19), Mαj and Dαj denote mean and standard deviations of αj :

[ ] 1 ∑R
Mαj χjα = αjr (V k ) ≡ µj , (20)
R r=1

7
 √ [(
[ ] )2 ]
Dαj χjα = Mαj χjα − µj . (21)

The convergence criterion is chosen such that

ϵrel ≤ ϵ∗rel , (22)

where ϵ∗rel is a given tolerance error. In the present work, we choose ϵ∗rel = 0.05.

3.2.2 Convergence of empirical parameters with respect to the volume

Different approaches have been proposed previously to determine convergence

criteria related to the convergence with respect to the volume of the RVE, and
associated estimates for the RVE volume (see Kanit et al. [15], Gitman et al.
[9], or Pelissou et al. [20]). In [15], the estimation of the variance with respect
to the volume size was modeled by a power law, requiring the identification of
a metric associated with the covariance. In the present work, as the material
is fully nonlinear, the convergence of α(V ) is simply estimated by a numerical
convergence analysis with respect to V as follows:

α(V k+1 ) − α(V k )

<δ (23)
∥α(V k+1 )∥

where δ is a tolerance parameter.

Remark a more straightforward approach could have been alternatively pro-

posed, consisting in directly identifying the parameters of a chosen empirical
law for the homogenized composite, without employing the incremental ho-
mogenization procedure to obtain σ̂(V k , ξ r , t). However, we have found two
drawbacks for this option: (a) first, without any knowledge about the ef-
fective constitutive law, the range of parameters used for the identification
can be very large, leading to heavy computations associated to minimiz-
ing (17); (b) a classical constitutive law can lead to inaccurate results with
respect to the reference solution, as illustrated below in Fig. 2. Another
option would be to analyze the macroscopic stress response of the compos-
ite by applying directly the incremental homogenization procedure, and by
using a constitutive law with fixed values for the matrix. As shown below,
this can also lead to inaccurate results with respect to the reference so-
lution. For illustration, we compare in Fig. 2 the effective stress response
obtained by a direct computation by FEM (reference solution) for the non-
linear composite described in section 4 with the following estimates: i) the
stress response obtained by directly applying the incremental procedure;
ii) the stress response obtained by identifying an empirical constitutive law
without the incremental homogenization scheme and iii), the stress response
obtained by the proposed corrected incremental homogenization scheme. We

8
 can note that the last solution is the closest to the reference solution. The
case corresponds to the case of elastic inclusions in an elastoplastic matrix,
for a volume fraction f = 0.5 and material parameters provided in Table 1,
except that the Young moduli of inclusions and matrix are XX and YY, re-
spectively. Furthermore, this last strategy is also the most efficient, because
in that case the range for the parameters to be identified can be set much
narrower.

4 Numerical examples

In the following two examples, a 2D RVE is considered, containing randomly

distributed cylindrical fibers, as depicted in figure 1. The fibers cannot in-
terpenetrate and the inclusions are supposed to be elastic, while the matrix
has a nonlinear behavior. We study two cases: an elastoplastic matrix and
a elasto-viscoplastic matrix. For each case, two volume fractions are studied:
f = 0.3 and f = 0.5. For a given number of fibers and a given volume frac-
tion, we define the size of the RVE√as the length of one edge of the square
domain containing the fibers, L = (πN )/f D/2, D being the diameters of
the inclusions.

4.1 Composite with elastoplastic matrix

In this first example, the matrix is assumed to have an elastoplastic behav-

ior, described by J2 -flow theory with nonlinear isotropic hardening. In that
context, the constitutive law is given by:

σ = C : (ε − εp ) (24)

where εp is the plastic strain( and where)C the elastic tensor, assumed to be
isotropic, i.e. C = κ1 ⊗ 1 + µ I − 13 1 ⊗ 1 , with κ the bulk modulus expressed
by κ = E/(3(1 − 2ν)) and µ the shear coefficient expressed by µ = E/2(1 + ν).

In the context of the J2 theory, the response is assumed to be linear if

f (σ, p) = J2 (σ) − σY − R(p) ≤ 0 (25)

( )
where J2 (σ) = 23 s : s , s = σ − 13 T r(σ), T r(.) being the trace operator, σY is
the yield stress, R(p) is the hardening stress and p the cumulated equivalent
plastic strain express by:
∫ t
p(t) = ṗ(τ )dτ (26)
0

9
 Matrix Inclusion

Young modulus G (MPa) 15000 300000

Poisson coefficient ν 0.2 0.45
Yield Stress σY (MPa) 80
Hardening modulus σ∞ (MPa) 60
Hardening exponent m 75
Linear hardening modulus k (MPa) 20
Table 1
Material parameters for the matrix and the inclusions for the composite with elast-
plastic model.

( )1/2
2 p
with ṗ = 3
ε̇ : ε̇p .

The evolution of the plastic strain εp is given by the normality rule:

∂f
ε̇p = ṗ (27)
∂σ

In this work, the hardening strain is assumed to be in the form:

( )
R(p) = kp + σ∞ 1 − e−mα (28)

where k is the linear hardening modulus and σ∞ is the hardening modulus,

and m is hardening exponent. The numerical parameters of the model are
indicated in Table 1.

As described in the previous sections, Finite Element computations are per-

formed by prescribing a tensile loading ε = ε11 (t)e1 ⊗ e1 on the RVE, in order
to provide a reference solution. Illustrations of some von Mises stress fields
for cases involving 4 and 144 inclusions, with volume fractions f = 0.3 and
f = 0.5 are provided in Figs. 3 and 4. We have observed highly localized shear
bands for small values of L/D. However, for much larger values of L/D, these
bands exist but are stopped by surrounding inclusions and remain limited.
This induces a lack of ergodicity of the process, which requires high sizes of
RVE (ratio L/D) to obtain converged values of parameters, as shown in the
following. This phenomenon has also been shown in [8].

The procedure described in section 3 is carried out for the two cases, f = 0.3
and f = 0.5. The statistical convergence of the parameters σY , σ∞ , m and k is
depicted in Figs. 5-6. The respective errors ϵ∗rel have been chosen respectively
as ϵ∗rel = 0.02, ϵ∗rel = 0.02, ϵ∗rel = 0.03 and ϵ∗rel = 0.04.

10
 (a) (b)
Figure 3. Elastoplastic composite: von Mises stress field in the case of volume frac-
tion f = 0.3: (a) 4 inclusions, corresponding to L/D =3.23; (b) 144 inclusions,
corresponding to L/D = 19.41

(a) (b)
Figure 4. Elastoplastic composite: von Mises stress field in the case of volume frac-
tion f = 0.5: (a) 4 inclusions, , corresponding to L/D =2.5066; (b) 144 inclusions,
corresponding to L/D = 15.03

80 70

70 60
60
50
σ∞ (MPa)
σY (MPa)

50
40 Mean
Mean
40 Standard deviation
Standard deviation
30
30
20
20

10 10

0 0
0 20 40 60 80 100 0 20 40 60 80 100
R R

(a) (b)
Figure 5. Elastoplastic composite, f = 0.3: statistical convergence of parameters:
(a) σY and (b) σ∞ for a fixed size of RVE L/D = 3.23.

In figures 7-10, the convergence of the parameters with respect to the RVE
size L/D is depicted. For each value of L/D, we have reported the converged
value with respect to the number of realization R.

11
 100 25

80 20

60 15
Mean Mean
m

k
Standard deviation Standard deviation
40 10

20 5

0 0
0 20 40 60 80 100 0 20 40 60 80
R R

(a) (b)
Figure 6. Elastoplastic composite, f = 0.3: statistical convergence of parameters:
(a) m and (b) k for a fixed size of RVE L/D = 3.23.

24 22

22 20
k (MPa)

k (MPa)

20 18

18 16

16 14
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16
L/D L/D

(a) (b)
Figure 7. Elastoplastic composite: convergence of the parameter k with respect to
the volume of the RVE, or equivalently with respect to L/D; (a) volume fraction
f = 0.3; (b) volume fraction f = 0.5.

76 67.5

74
σY (MPa)

σY (MPa)

65.5

72

70 63.5
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16
L/D L/D

(a) (b)
Figure 8. Elastoplastic composite: convergence of the parameter σY with respect to
the volume of the RVE, or equivalently with respect to L/D; (a) volume fraction
f = 0.3; (b) volume fraction f = 0.5.

12
 96
84

82

80
m

m
94
78

76

74

72 92
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16
L/D L/D

(a) (b)
Figure 9. Elastoplastic composite: convergence of the parameter m with respect to
the volume of the RVE, or equivalently with respect to L/D; (a) volume fraction
f = 0.3; (b) volume fraction f = 0.5.

70 66

68 64
σ∞ (MPa)

σ∞ (MPa)

66 62

64 60

62 58
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16
L/D L/D

(a) (b)
Figure 10. Elastoplastic composite: convergence of the parameter σ∞ with respect
to the volume of the RVE, or equivalently with respect to L/D; (a) volume fraction
f = 0.3; (b) volume fraction f = 0.5.

In the case of the volume fraction f = 0.3, we can observe a convergence of

the different parameters for RVE size of roughly L/D = 17 − 18. In the case
f = 0.5, the RVE size is roughly L/D = 14 − 15.

4.2 composite with Elasto-viscoplastic matrix

In this second example, the matrix is assumed to have an elasto-viscoplastic

behavior. The total strain ε is assumed to be the sum of an elastic part and
a viscoplastic part εvp as:

ε = εe + εvp . (29)

13
 Matrix Inclusion

Young modulus G (MPa) 15000 300000

Poisson coefficient ϑ 0.2 0.45
Yield Stress σY (MPa) 80
Hardening modulus σ∞ (MPa) 60
Hardening exponent m 75
Linear hardening modulus k (MPa) 20
Viscoplastic coefficient η 2
Viscoplastic exponent s 100
Table 2
Material parameters for the matrix and the inclusions for the composite with elasto-
viscoplastic model

The viscous effects are only considered in the plasticity domain. The evolution
of the viscoplastic strain is given by the flow rule:

∂f
ε̇vp = ṗ (30)
∂σ
and the stress satisfies (25). In this work, a current yield stress Norton law
[19] has been chosen to describe the evolution of ṗ = dp
dt
:

( )s
σY f
ṗ = (31)
η σY + σ(p)

where η is the viscoplastic Norton coefficient. This parameter indicates the

viscoplastic sensitivity of the material to the strain rate, s is the viscoplastic
exponent and f is the viscoplastic part of the stress, defined as :

f = J2 (σ) − σY − R(p). (32)

The numerical values for these parameters are indicated in Table 2.

Here again, the procedure described in section 3 is carried out. Illustrations

of von Mises stress fields obtained during the computation of the reference
solution are depicted in figures 11 and 12, corresponding to f = 0.3, f = 0.5,
respectively, for cases involving N = 4 and N = 144 inclusions.

14
 (a) (b)
Figure 11. Elasto-viscoplastic composite: von Mises stress field in the case of volume
fraction f = 0.3: (a) 4 inclusions, corresponding to L/D =3.23; (b) 144 inclusions,
corresponding to L/D = 19.41

(a) (b)
Figure 12. Elasto-viscoplastic composite: von Mises stress field in the case of vol-
ume fraction f = 0.5: (a) 4 inclusions, , corresponding to L/D =2.5066; (b) 144
inclusions, corresponding to L/D = 15.03

In this case of an elasto-viscoplastic composite, we also observe shear bands,

which propagate outside of the cell for small sizes L/D and are blocked for
larger sizes. However, these bands are much less localized in the viscoplastic
case and induce a faster convergence with respect to the size L/D, and thus
leading to smaller RVE sizes.

In figures 13-18, the convergence of the parameters with respect to L/D is

depicted, as in the previous example.

We can note that the convergence of the hardening mudulus k and of the
viscoplastic coefficient η is quite fast with respect to the size L/D, while the
other coefficients converge much slower. Eventually, the RVE size L/D for
which a convergence is observed with respect to all parameters is in this case
about L/D = 14 − 15. The faster convergence in the case of composites with
elasto-viscoplastic matrix can be explained by less localization in the matrix
as compared to the elastoplastic case.

15
 22 22

20
k (MPa)

k (MPa)
20

18

18 16
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16
L/D L/D

(a) (b)
Figure 13. Elasto-viscoplastic composite: convergence of the parameter k with re-
spect to the volume of the RVE, or equivalently with respect to L/D; (a) volume
fraction f = 0.3; (b) volume fraction f = 0.5.

68 68

66
σY (MPa)

σY (MPa)

66

64

64 62
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16
L/D L/D

(a) (b)
Figure 14. Elasto-viscoplastic composite: convergence of the parameter σY with
respect to the volume of the RVE, or equivalently with respect to L/D; (a) volume
fraction f = 0.3; (b) volume fraction f = 0.5.

76 64

74
62
σ∞ (MPa)

σ∞ (MPa)

72

60
70

68 58
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16
L/D L/D

(a) (b)
Figure 15. Elasto-viscoplastic composite: convergence of the parameter σ∞ with
respect to the volume of the RVE, or equivalently with respect to L/D; (a) volume
fraction f = 0.3; (b) volume fraction f = 0.5.

16
 112 104
110
108 102
106
104 100
102
98
100
s

s
98
96
96
94 94
92
90 92
88
86 90
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16
L/D L/D

(a) (b)
Figure 16. Elasto-viscoplastic composite: convergence of the parameter s with re-
spect to the volume of the RVE, or equivalently with respect to L/D; (a) volume
fraction f = 0.3; (b) volume fraction f = 0.5.

86 92

90
84
88

86
m

82
84

82
80
80

78 78
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16
L/D L/D

(a) (b)
Figure 17. Elasto-viscoplastic composite: convergence of the parameter m with re-
spect to the volume of the RVE, or equivalently with respect to L/D; (a) volume
fraction f = 0.3; (b) volume fraction f = 0.5.

2.8 3.6
3.4
2.6 3.2
3
2.4 2.8
2.6
2.2 2.4
2.2
η

2 2
1.8
1.8 1.6
1.4
1.6 1.2
1
1.4 0.8
0.6
1.2 0.4
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16
L/D L/D

(a) (b)
Figure 18. Elasto-viscoplastic composite: convergence of the parameter η with re-
spect to the volume of the RVE, or equivalently with respect to L/D; (a) volume
fraction f = 0.3; (b) volume fraction f = 0.5.

17
5 Conclusion

A new procedure to determine the size of RVE for nonlinear random materials
has been proposed. The technique is based on analyzing the convergence of
parameters used in an incremental homogenization framework [4] and identi-
fied by means of FEM nonlinear finite element computations, with respect to
both statistical realization of the random microstructure and to the size of the
unit cell. Two cases were studied, involving fiber-composites with hard elastic
cylindrical fibers and elastoplastic and elasto-viscoplastic matrix, for different
volume fractions. We have found a convergence of the parameters leading to
the existence of an RVE for these cases, with a size estimated as 17-18 and
13-15 times the diameter of the inclusions for the elastoplastic and viscoplas-
tic cases, respectively. The faster convergence in the case of composites with
elasto-viscoplastic matrix can be explained by less localization in the matrix as
compared to the elastoplastic case. Once the RVE size is identified, the related
parameter used in the incremental homogenization framework provide a very
accurate homogenization model, useable in efficient structure calculations.

6 Aknowledgements

The financial support this work enjoys from Ile-de France Region, BpiFrance,
DGE and ILMAB Team is gratefully aknowledged.

References

[1] W.L. Azoti, A. Tchalla, Y. Koutsawa, A. Makradi, G. Rauchs, S. Bellouetar,

and H. Zahrouni. Mean-field constitutive modeling of elasto-plastic composites
using two (2) incremental formulations. Compos. Struct., 105:256–262, 2013.
[2] V. Bulsara, R. Talreja, and J. Qu. Damage initiation under transverse loading
of unidirectional composites with arbitrarily distributed fibers. Compos. Sci.
Technol., 59:673–682, 1999.
[3] J. Cugnoni and M. Galli. Representative volume element size of
elastoplastic and elastoviscoplastic particle-reinforced composites with random
microstructure. Comput. Model. Eng. Sci., 66(2):165–185, 2010.
[4] I. Doghri and A. Ouaar. Homogenization of two-phase elasto-plastic composite
materials and structures, study of tangent operators, cyclic plasticity and
numerical algorithms. Int. J. Sol. Struct., 40:1681–1712, 2003.
[5] X. Du and M. Ostoja-Starzewski. On the size of representative volume element
for darcy law in random media. Proc. R. Soc. A, 462:2949–2963, 2006.

18
[6] G. Dvorak. Micromechanics of Composites Materials. Springer, New-York,
2013.

[7] A. Elvin and S.S. Sunder. Microcracking due to grain boundary sliding in
polycrystalline ice under uniaxial compression. Acta Mater., 44(1):43–56, 1996.

[8] F. Fritzen, S. Forest, T. Böhlke, D. Kondo, and T. Kanit. Computational

homogenization of elasto-plastic porous metals. Int. J. Plast., 29(0):102 – 119,
2012.

[9] I.M. Gitman, H. Askes, and L.J. Sluys. Representative volume: Existence and
size determination. Eng. Frac. Mech., 74(16):2518 – 2534, 2007.

[10] C. Gonzalez and J. LLorca. Mechanical behavior of unidirectional fiber-

reinforced polymers under transverse compression: Microscopic mechanisms and
modeling. Compos. Sci. Technol., 67(13):2795 – 2806, 2007.

[11] Q. Grimal, K. Raum, A. Gerisch, and P. Laugier. A determination of the

minimum sizes of representative volume elements for the prediction of cortical
bone elastic properties. Biomech. Model. Mechanobiol., 10:925–937, 2011.

[12] A.A. Gusev. Representative volume element size for elastic composites: a
numerical study. J. Mech. Phys. Solids, 45:1449 –1459, 1997.

[13] C. Heinrich, M. Aldridge, A.S. Wineman, J. Kieffer, A.M. Waas, and

K. Shahwan. The influence of the representative volume element (rve) size on
the homogenized response of cured fiber composites. Modelling Simul. Mater.
Sci. Eng., 20:075007, 2012.

[14] R. Hill. A self-consistent mechanics of composites materials. J. Mech. Phys.

Solids, 13:213–222, 1965.

[15] T. Kanit, S. Forest, I. Galliet, V. Mounoury, and D. Jeulin. Determination of

the size of the representative volume element for random composites: statistical
and numerical approach. Int. J. Solids Struct., 40(13-14):3647 – 3679, 2003.

[16] Z. F. Khisaeva and M. Ostoja-Starzewski. On the size of rve in finite elasticity

of random composites. J. elasticity, 37, 2006.

[17] I. Monetto and W.J. Drugan. A micromechanics-based non local constitutive

equation and minimum rve size estimates for random elastic composites
containing aligned spheroidal heterogeneities. J. Mech. Phys. Solids, 57:1578–
1595, 2009.

[18] M. Ostoja-Starzewski. Material spatial randomness: from statistical to

representative volume element. Probab. Eng. Mech., 21:112–132, 2006.

[19] N.S. Ottosen and M. Ristinmaa. Mechanics of Constitutive Modeling. Elsevier

Science Ltd, 2005.

[20] C. Pelissou, J. Baccou, Y. Monerie, and F. Perales. Determination of the size

of the representative volume element for random quasi-brittle composites. Int.
J. Solids Struct., 46(14-15):2842 – 2855, 2009.

19
[21] V. Pensée and Q.-C. He. Generalized self-consistent estimation of the apparent
isotropic elastic moduli and minimum representative volume element size of
heterogeneous media. Int. J. Sol. Struct., 44(7):2225–2243, 2007.
[22] J.F. Rakow and A.M. Waas. Size effects in metal foam cores for sandwich
structures. AIAA J., 42:7, 2004.
[23] I.S. Ranganathan and M. Ostoja-Starzewski. Scaling function, anisotropy and
the size of rve in elastic random polycrystals. J. Mech. Phys. Solids, 56:2773–
2791, 2008.
[24] Z.Y. Ren and Q.S. Zheng. A quantitative study of minimum sizes of
representative volume elements of cubic polycrystals - numerical experiments.
J. Mech. Phys. Solids, 50:881–893, 2002.
[25] S. Swaminathan S and S. Ghosh. Statistically equivalent representative volume
elements for unidirectional composite microstructures: Ii. with interfacial
debonding. J. Compos. Mater., 49:605–621, 2006.
[26] K. Sab and B. Nedjar. Periodization of random media and representative
volume element size for linear composites. C. R. Mech., 333:187–195, 2005.
[27] A. Salahouelhadj and H. Haddadi. Estimation of the size of the rve for isotropic
copper polycrystals by using elastic-plastic finite element homogenisation.
Computat. Mater. Sci., 48(3):447 – 455, 2010.
[28] M. Salmi, F. Auslender, M. Bornert, and M. Fogli. Various estimates of
representative volume element sizes based on a statistical analysis of the
apparent behavior of random linear composites. C. R. Meca., 340:230–246,
2012.
[29] J.C. Simo and T.J.R. Hugues. Computational Inelasticity. Springer-Verlag,
1998.
[30] C. Soize. Tensor-valued random fields for meso-scale stochastic model of
anisotropic elastic microstructure and probabilistic analysis of representative
volume element size. Probab. Eng. Mech., 23:307–323, 2008.
[31] S. Swaminathan, S. Ghosh, and N.J. Pagano. Statistically equivalent
representative volume elements for unidirectional composite microstructures:
I. without damage. J. Compos. Mater., 40:583–604, 2006.
[32] P.K. Valavala, G.M. Odegard, and E.C. Aifantis. Influence of representative
volume element size on predicted elastic properties of polymer materials.
Modelling Simul. Mater. Sci. Eng., 17:045004, 2009.
[33] W.J. Drugan J.R. Willis. A micromechanics-based nonlocal constitutive
equation and estimates of representative volume element size for elastic
composites. J. Mech. Phys. Solids, 44(4):497–524, 1996.
[34] T.I. Zhodi and P. Wriggers. On the sensitivity of homogenized material
responses at infinitesimal and finite strains. Commun. Numer. Methods Eng.,
16:657–670, 2000.

20