Multi-Scale Failure For Heterogeneous Materials: Link With Morphological Modeling - Complas Xi
Multi-Scale Failure For Heterogeneous Materials: Link With Morphological Modeling - Complas Xi
Multi-Scale Failure For Heterogeneous Materials: Link With Morphological Modeling - Complas Xi
2
e
xy/Lc
where L
c
is the correlation length.
The orthogonal decomposition of Gaussian correlated random elds theory stipulates
[8] that mean zero Gaussian eld with continuous covariance function (such as C) can be
written as follows
(x, w) =
n=1
n
(x)
n
(w), (1)
where
n
(w) are zero mean, unit variance Gaussian random variables, and
n
(x) are
functions on M determined by the covariance function C. It is worth noting that eq.(1)
2
E. Roubin, M. Bogdan and J.B. Colliat
allows for stochastic - w - and spatial -x - variables separation. Therefore, implementing
this framework comes to put the eort in the determination of the spatial functions
n
(x).
The Karhunen-Lo`eve decomposition is based on the previous orthogonal decomposi-
tion. It allows us to determine these spatial functions
n
(x) for simple compact M in
R
N
. Demonstration can be found in [9] that they can be determined by rst solving the
following eigenvalues problem (known as Fredholm problem):
_
M
C(x, y)(y)dy = (x) (2)
where and are respectively the eigenvalues and eigenvector and then by setting
n
(x) =
n
(x). Theoretically, an innite sum is needed to dene the exact random
eld in eq.(1). For the numerical implementation made here, a Finite Element Method is
used to solve a discretized Fredholm problem. Therefore, using a nite set of eigenvalues
and eigenvectors, the following troncated Karhunen-Lo`eve decomposition eq.(3) denes
an approximative realization of the underlying random eld.
(x, w) =
m
n=1
_
n
(w)
n
(x). (3)
The fact that stochastic and spatial variables are still separated is an essential result
for any numerical implementation. Indeed, once the m couples {
n
;
n
} of a certain
correlated random eld are determined, the generation of a realization comes to generate
a set of independent Gaussian variables (which only requires a random number generator).
Moreover, the same couples can be used to produce any other realizations of the same
eld.
The precision of this method, involving full squared matrix eigenvalues problem, is
quickly limitated by the memory storage when one deals with multi-dimensional random
elds of large size. The turning bands projectional method has been developped by Math-
eron [7] in order to reduce the amount of numerical ressources. The idea is to generate
several one-dimensional realizations of random elds to produce a multi-dimensional one.
The algorithm below explains this projectional method with details.
Let M be the discreted multi-dimensional bounded region where the nal realization
will be created. Several lines have to be generated (we shall call L their number) with
one arbitrary intersection point 0 and an uniform distribution of directions over the unit
ball (see Fig.1).
Let z(, w
i
), i = 1..L be the L realizations of a one-dimensional correlated random eld
generated over the L lines. For each point N on M, the value of the multi-dimensional
realization is the average of the one-dimensional realization values at the projection of N
on each line i:
(N, w) =
1
L
L
i=1
z(
N
i
, w
i
) (4)
3
E. Roubin, M. Bogdan and J.B. Colliat
0
z(
Ni
)
line i
u
i
N
M
N
i
i
x
Ni
Figure 1: Schematic representation of the turning band method (from [10])
In this paper, the application of the method is made for three-dimensional random
elds. The key of this method is the link between the three-dimensional covariance func-
tion C and the equivalent one-dimensional covariance function C
1
we need to generate
the L realizations. Let C(r) be as above (with r = x y). Following [7] we have
C
1
(r) =
d
dr
(rC(r)) =
2
_
1
2r
2
L
2
c
_
e
r
2
/L
2
c
(5)
3 Excursion Set
We call an excursion set the morphology of a subset of a bounded region dened by
thresholding a realization of a random eld. It allows us to create a set of random shapes.
Let be a realization of (x, w) : M R
N
R dene as above and u R a chosen
threshold. The underlying excursion set A
u
is dened by the points of M where the values
of are above u (eq.(6)).
A
u
A
u
(, M) {x M : (x) u} (6)
This principle, applied for M R is shown on Fig.2.
In our case, random elds will be yield in a three dimensional space (M R
3
) and
therefore dene three-dimensional excursion sets. The two excursions represented in Fig.3
are made from the same realization with two dierent threshold values. It is clear that, by
changing this value, a large range of varied morphologies can be generated. This exemple
shows that low values of u produce excursions mainly made of handles with high volume
fraction, giving a sponge like topology (Fig.3(a)), whereas high values of u produce
excursion made of several connected components with a lower volume fraction (Fig.3(b)).
In order to provide a global description of the resulting morphology, the Lipschitz-
Killing curvatures, hereafter LKCs, are choosen. In a N-dimensional space N + 1 LKCs
can be dened where each can be thought of measures of the j-dimensional sizes of
4
E. Roubin, M. Bogdan and J.B. Colliat
x
u
A
u M
Figure 2: Schematic representation of a one-dimensional excursion set A
u
(a) Low threshold - sponge
topology
(b) High threshold - meatball
topology
Figure 3: Eect of threshold value on tri-dimensional excursion topology
A
u
. In our three-dimensional case, the four LKCs, denoted by L
j
, j = 0..3, provide both
geometrical - L
1
, L
2
, L
3
- and topological - L
0
- descriptions of the morphology A
u
. They
are dened by:
- L
3
(A
u
) is the three dimensional volume of A
u
.
- L
2
(A
u
) is half the surface area of A
u
.
- L
1
(A
u
) is twice the caliper diameter of A
u
.
- L
0
(A
u
) is the Euler characteristic of A
u
, which contrary to the other LKCs is a
topological measure. In three-dimension, it can be calculated by:
L
0
(A
u
) = #{connected components in A
u
}#{handles in A
u
}+#{holes in A
u
}
For exemple, a ball or a cube are topologicaly identical (Euler characteristic L
0
= 1)
but dier from a hollow ball (L
0
= 2) or a ring torus (L
0
= 0).
5
E. Roubin, M. Bogdan and J.B. Colliat
Following [2], a probabilistic link has been made between excursion set properties and
random eld thresholding parameters giving an explicit formulae for the expectation of
the LKCs - E{L
i
(A
u
(, M))}. It is not the purpose of this paper to give details on these
formulae, however, full proof and details can be found in [9]. The only idea one need
to remember to go through this paper is that this theory gives a new tool helping us to
predict all the geometrical and topological properties of an excursion set from the random
eld characteristics and the threshold - , L
c
, u -. These relations have been made explicit
for (x, w) as above on a cube M =
3
i=1
[0; T]:
_
_
E{L
0
}(A
u
) =
_
2
2
2
T
3
L
3
c
_
u
2
2
1
_
+
3
2
2
3/2
T
2
L
2
c
u
+
3
2
2
T
Lc
_
e
u
2
/2L
2
c
+
_
u
_
E{L
1
}(A
u
) =
_
2
3/2
T
3
L
2
c
u
+
3
2
4
T
2
Lc
_
e
u
2
/2L
2
c
+ 3T
_
u
_
E{L
2
}(A
u
) =
T
3
Lc
e
u
2
/2L
2
c
+ 3T
2
_
u
_
E{L
3
}(A
u
) = T
3
_
u
_
(7)
Fig.4(a) and Fig.4(b) represent respectively the Euler characteristic and the volume
fraction - directly linked with the fourth LKC by E{L
3
}(A
u
)/T
3
- of excursion sets of
(x, w) for u from 20 to 20.
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
-20 -15 -10 -5 0 5 10 15 20
V
o
l
u
m
e
f
r
a
c
t
i
o
n
[
-
]
Threshold [-]
(a)
-80
-60
-40
-20
0
20
40
60
-20 -15 -10 -5 0 5 10 15 20
E
u
l
e
r
C
a
r
a
c
t
e
r
i
s
t
i
c
[
-
]
Threshold [-]
(b)
Figure 4: LKCs of excursion sets of Gaussian random eld in term of threshold values.
Expected values of LKCs provided by (7), + Numerical values calculated from one
realization of (x, w).
The constant decreasing shape of the volume fraction curve in term of u clearly reects
the eect of the threshold level on the size of A
u
. Even if more peculiar, the Euler
characteristic curve shape reects also easily the eect of the threshold on excursion sets
topology. For values of u lower than the lowest value of , the Euler characteristic is
the one of the full cube (L
0
= 1). By increasing u, several holes appear, counting in
positive for the Euler characteristic (L
0
> 1). Then, the expansion of the holes starts to
form handles which lead to a sponge like topology (L
0
< 0). By increasing u even more,
6
E. Roubin, M. Bogdan and J.B. Colliat
handles disappear forming a meatball like topology of connected components (L
0
> 0).
Finally, the Euler characteristic decreases to L
0
= 0 when no more connected components
remain.
From the comparison between theorical values and measures on one realization, we can
point out that the variability of the numerical generation is very low. Therefore, although
eq.(7) gives only expectations of LKCs, for this range of excursion sets we can assume
that V{L
i
(A
u
)} 1.
So far, we have seen the eect of the threshold value on excursion sets. But one needs
to remember that, according to eq.(7), both variance and covariance length of (x, w)
aect the morphology as well. Understanding the full behaviour of these equations is a
key point for anyone who wants to make excursion set modeling.
4 Application of the modeling framework on concrete like material
The material is represented as an heterogeneous material with two phases. One phase
(aggregates) is represented by an excursion set of a correlated random eld while its second
phase (concrete) is represented by its complementary. Therefore in this part, the eort
will be put in a realistic representation of the aggregates phase. We keep only three
relevant characteristics from the four LKCs: the volume fraction V
v
, the volumic surface
area S and the number of agregates N which are respetively linked with L
3
, L
2
and L
0
.
Thought V
v
and S can be directly estimated, attention must be taken when it comes to
N. Indeed the Euler characteristic does not indicate the number of aggregates for every
topology. In our case, the meatball topology has to be targeted and it is only once we
assume that the excursion set is free from holes and handles that N can be estimated by
L
0
. In this specic kind of topology: N #{connected components} = L
0
.
Once the three characteristics (N, S, V
v
) of the phase are set, the generation of the
underlying excursion set rely on nding a solution for (u, , L
c
) that satisfy the following
system:
_
_
_
E{L
3
}(u, ) = V
v
T
3
E{L
2
}(u, , L
c
) =
1
2
ST
3
E{L
0
}(u, , L
c
) = N
(8)
Due to the intrinsic non linearity of eq.(7), depending on the dierent values of (N, S, V
v
)
(especially for meatball topology - N 1) the problem eq.(8) do not always have a
solution. For exemple, we can clearly see on Fig.4 that we can not expect N to be up-
per than 40 while keeping a high volume fraction (V
v
> 40%). Which in our case of
concrete like material modeling leads to a major issue. So far, the more realistic solution
for meatball topology we get with this framework allows us to represent an aggregate
phase with a maximum of 15% volume fraction.
Until now, the framewok has been presented considering Gaussian random elds. But
estimation of LKCs for excursion set can also be worked out considering Gaussian related
elds. The application of this paper is made using a chi-square distribution with k degrees
7
E. Roubin, M. Bogdan and J.B. Colliat
of freedom -
2
k
-. Realizations of such elds can be seen as sum of k independent squared
realizations of a correlated Gaussian random eld. Let be a realization of such eld and
i
, i = 1..k be k realizations of the Gaussian eld (x, w) described above. We have :
=
k
i=1
2
i
(9)
Although similar to eq.(7), the use of a
2
k
distribution add the parameter k to the
system eq.(8). With such eld, the nearest solution is found for k = 1 and enable us to
double the previous volume fraction V
vmax
30%. Fig.5 shows a two-dimensional slice of
excursions from a Gaussian realization and a
2
1
made from the same realization. Fig.5(b),
being the excursion from the squared realization of the excursion Fig.5(a), shows clearly
that, for the same threshold, it is natural to expect the volume fraction to double between
excursions of Gaussian and
2
1
random elds.
(a) Gaussian realization - (b)
2
1
realization - =
2
Figure 5: Comparison between Gaussian and
2
1
excursion sets for the same threshold value.
The
2
1
distribution remains the more suitable solution for meatball topology and high
volume fraction morphology we found.
5 FE model for heterogeneous material - Application to hydration process
modeling
The approach made here relies on a spatial truss, to model pattern of heterogeneities.
The choice of a not adapted meshing process is made here thus, the spatial positions of
nodes are not constrained by the morphology. Therefore, both gemetrical and mechanical
properties have to be handle inside some interface elements. These cut elements are split
into two parts, each having dierent elastic properties by enhancing them with strain
(weak) discontinuities [11]. An elementary enhancements method (E-FEM) method for
kinematic enhancement of Finite Element using the Hu-Washizu variational formulation is
used here. For example, if we consider a two-phase material (inclusions within a matrix),
8
E. Roubin, M. Bogdan and J.B. Colliat
three sets of elements are needed: those entirely in the matrix, those entirely in the
inclusions, and those which are split between both (cut elements). To calculate these
elements repartition, a projection of the previous excursion set is made onto the truss.
In order to illustrate this linear framework, a simple hydration process of concrete like
material modeling has been implemented. Considering a simplistic version of the Powers
and Brownyard hydration model [12], with only three phases: unreacted cement, hydration
products (including gel water) and free water, the volume fraction of each one of them
can be calculated according to the following equations:
_
_
p =
w/c
w/c+w/c
V
anh
= (1 p)(1 )
V
h
= 2.12(1 p)
V
w
= 1 V
h
V
anh
(10)
where p is the initial porosity, the hydration degree and V
anh
, V
h
, V
w
respectively the
volume fractions of anhydrous cement, hydration products and water.
(a) = 0.1 (b) = 0.5 (c) = 1
Figure 6: Projection of excusion set shapes on FE truss for dierent hydration degrees.
water, hydration products, anydrous cement
As explained previously, thresholding a random eld with a scalar allows to create a
two phase material. One can easily imagine, that a second threshold, with a dierent
value, will allow to create an additional phase, concentrical to the rst one. Therefore,
setting two thresholds will allow us to create a three phase material. Thus, for dierent
hydration degrees, each phases volume fraction is known and can be linked to the random
elds thresholds u
i
(equation eq.(7)). Eventually, the initial morphology is set up by one
threshold (two phases: water and unhydrated cement), and then, for a growing hydration
degree, two thresholds are calculated and applied to the random eld, creating a three
phase material (water, unhydrated cement and hydration products).
Within this framework, macroscopic material characteristics like Young modulus can
be estimated over a given hydration degree with simple tension tests. The following
9
E. Roubin, M. Bogdan and J.B. Colliat
characteristics have been chosen E
anh
= 135 000 MPa, E
h
= 25 000 MPa and E
w
=
1 MPa.
0
5000
10000
15000
20000
25000
30000
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Y
o
u
n
g
m
o
d
u
l
u
s
[
M
P
a
]
Hydration degree [-]
Figure 7: Young modulus of a concrete like material for dierent hydration degrees.
Fig.7 shows that the continuous growing of the macroscopic Young modulus over hy-
dration degree is well handled by this FE representation. A slight raising of the slope can
be seen after = 0.4.
6 FE models with embedded discontinuities
In addition to the geometrical representation of heterogeneities, displacement (strong)
discontinuities are also introduced in the elements, in order to model a non-linear softening
response based on failure quasi-brittle. These discontinuities represent micro-cracks that
can occurs in both phases as well as at the interfaces (debonding). Details of this FE
numerical implementation can be found in [3].
A other simple tension test is presented here. Material properties are dened according
to Tab.1.
Table 1: Material properties
Matrix Inclusions Interface
E = 10GPa 70GPa
u
= 3MPa 3MPa
Gf = 11J/m
2
11J/m
2
Two remarks are worthy of attention. The rst is that the interface is of rigid-brittle
type. The second is that we choosed for inclusions to remains in the linear elastic regime.
10
E. Roubin, M. Bogdan and J.B. Colliat
(a) Displacement eld and crack pat-
tern at last time step
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04
L
o
a
d
[
k
N
]
Displacement [mm]
(b) Load vs imposed displacement
Figure 8: Results for simple tension
The cracking pattern is shown on Fig.8(a) where two zones are splited by a macroscopic
crack (represented by means of the broken elements). Fig.8(b) shows the macroscopic
load vs imposed displacement curve where three steps can be seen. First, a linear part
where no failure occurs. Then, with the apparition of several microscopic cracks, we can
observe a yield behaviour. Finaly, the softening part begin when the localisation of these
microscopic cracks creates a macrosopic one.
7 Concluding remarks
This communication presents a rst attempt to create a sequential multi-scale frame-
work where morphology of heterogeneous material is dened by excursion sets of correlated
random elds. Though, eorts still have to be made in order to generate more realistic
morphologies, advantages have been shown through two examples. We can also add that
this framework is well adapted to other problematics related with concrete like materials
such as the eect of morphological variability on macroscopic behaviour. Indeed, the
use of both Karhunen-Lo`eve decomposition and non-adapted meshes allows fast compu-
tations, limiting the growing amount of numerical ressources needed when dealing with
large sets of morphologies. Futhermore, being able to represent broken elements by means
of a strong discontinuity in the FE method allows calculations of permeability or diusion
in such damaged materials [13].
REFERENCES
[1] Feyel, F. and Chaboche, J.-L. Multi-scale non-linear FE analysis of composite struc-
ture: damage and ber size eects. Revue europeenne des
Elements Finis : NUM-
DAM00 isues (2001) 10:449472
11
E. Roubin, M. Bogdan and J.B. Colliat
[2] Adler, R.J. Some new random eld tools for spatial analysis. Stochastic Environmen-
tal Research and Risk Assessment (2008) 22:809822
[3] Benkemoun, N., Hautefeuille, M., Colliat, J.-B., Ibrahimbegovic, A. Failure of het-
erogeneous materials: 3D meso-scale FR models with embedded discontinuities. Int
J. Num. Meth. in Engng. (2010) 82:16711688
[4] Sukumar, N., Chapp, D.L., Moes, N. and Belytshcko, T. Modeling holes and in-
clusions by level sets in the extended nite element method. Computer Methods in
Applied Mechanics and Engineering (2001) 190:61836200
[5] Oliver J. Modelling strong discontinuities in solid mechanics via strain softening
constitutive equations. International Journal for Numerical Methods in Engineering
(1996) 39:35753623
[6] Loeve, M. Probability Theory. Graduated Texts in Mathematics (1978) Vol. II, 4th
ed..
[7] Matheron, G. The intrinsic random functions and their applications. Advances in
Applied Probability (1973) 5:439468
[8] Karhunen, K.
Uber lineare Methoden in der Wahrscheinlichkeitsrechnung. Ann.
Acad. Sci. Fennicae. Ser. A. I. Math.-Phys (1947) 37:179
[9] Adler, R.J. and Taylor, J.E. Random Fields and Geometry (2007) Springer, Boston
[10] Mantoglou, A. and Wilson, J.L The Turning Bands Method for Simulation of Ran-
dom Fields using Line Generation with a Spectral Method Water Resources Research
(1982) Vol.II, 2:129149
[11] Ortiz, M., Leroy, Y. and Needleman, A. A Finite Element method for localized failure
analysis Computer Methods in Applied Mechanics and Enginnering (1987) 61:189
214
[12] Powers, T.C. and Brownyard, T.L. Studies of the physical properties of hardened
Portland cement paste. J. Am. Concr. Inst., (1947). 43101132, 249336, 469505,
549602, 669712, 845880, 933992.
[13] Jourdain, X., Colliat, J.-B., De Sa, C., Benboudjema, F. and Gatuingt, F. Upscal-
ing permeability for fractured concrete : meso-macro numerical approach within a
sequential framework submitted
12