)..r:.-::‘ '

nun .m.
Au L «uxiu

guinzrg,

,..

. JAN. .3.-
,.. . 4hr-Vu
gm -

‘.. n.\~L

Lu.“ a
mum".

‘~- ‘-:£G=‘~.:i r’

l

_-.«.A.4.w.
It“ w.

x4 . A
551‘ Nunl‘
‘ . A. y“

’ ‘4', .
L‘A.\n.‘l."».‘ ; .
m.- w. . -~u.K.~"-nr.. m.
u "“‘IJ. ‘ . k .|.I.". ‘ ’ .
nun. ‘_ .\\A 2
. M L

x
‘1: -
..«1. "1. ' 1
~‘ ‘-§< . A , . .. >
h A . I l ‘ I. lA. V i 4
“1‘ - ~ m. h . ‘\ '
.. “ “i ‘ . 1‘03. ,
' .a .. ‘ _~ ‘ m M u..v....m‘|,«;-.-.
. Ad .‘n. \‘»‘<1.I.\.-I AV '.
--.. , t\'J n-.-.~. :1.
K K. mu, .4. W‘.
~
-.

‘:'=::.-n.. ‘
a; ' .Ictn‘u A

.u figk
~5~1.I- 4..

.>‘.. .-.-. “mm. \.
.'~~.~.\- 4\.<\"<« .mJ.‘
xv...» r. - . .u
,_. “-

—.~.-.u'..1r, ,
.V w.“ -.._.
ma. ,
:51.- M u
.\.S..‘x¢_\A~.‘,._~<-5 {
‘ f ‘ :.' 1‘“
~ .

. w =90.
Km}.

.‘ _ x.
Lst'. .n‘-.
. “um-4.;

tux-xx u v 1.1

mm.

-. -..‘ ‘ a
" uva

in. . v
.\-‘m
"‘1"

'-..—‘

mu”
' F ., ”:1‘,

, 5..

I,‘ 5'":

"'3“ I
.‘

.

...,. N
c .v

.1‘1‘fll

. u. .H—z
4 gr.

Myra:

:4'

" ‘7. .x ‘

.n

 

 

 

 

This is to certify that the

dissertation entitled

COMPUTATIONAL METHODS IN STOCHASTIC

MICROMECHANICS OF HETEROGENOUS SOLIDS

presented by 1

KHALID IBRAHIM ALZEBDEH

has been accepted towards fulfillment
of the requirements for

Ph.D. Mechanics

 

degree in

 

 

Date 8—4-94

 

MSU is an Affirmative Action/Equal Opportunity Instizution 0-12771

 

 

 

 

 

llllll

Till

@fi 1 0

\l\\

129

27

Till“

ll

 

 

_ LIBRARY
MWhigan State
Unwersnty

 

 

“km“ from Your record.

h
c
.B
m
e
m
m
m
X
0
Bu
.NHr
us
am
RF
mm
E0
V
WA
0
PT

DATE DUE DATE DUE DATE DUE

 

COMPUTATIONAL METHODS IN STOCHASTIC
MICROMECHANICS OF HETEROGENOUS SOLIDS

By

Khalid Ibrahim Alzebdeh

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Materials Science and Mechanics

1 994

ABSTRACT

COMPUTATIONAL METHODS IN STOCHASTIC
MICROMECHANICS OF HETEROGENOUS SOLIDS

By

Khalid Ibrahim Alzebdeh

Due to the spatial variability in heterogenous (random) materials at microscopic level,
its effective properties at various length scales have to be determined in a stochastic (prob-
abilistic) fashion. In this dissertation, such a determination is based on a discrete random
modeling of microstructure. Three problems of stochastic nature, are presented and meth-
odologies of treatment are proposed. Solutions are performed via computer simulations.
In the first problem, the effective elastic moduli of Delaunay networks, modeling two—
phase granular media are calculated. Combined with a Delaunay network, two spring
models are used to represent interactions between particles. In the first model, central
interaction is taken into account, while in the second one both central and angular interac-
tions are considered. Results of numerical simulations are used to identify that self-

consistent model which most closely approximates effective elastic properties of two-

phase Delaunay networks.

In the second problem, a micromechanics-based stochastic finite element method is
developed to account for the variability in material properties at micro level. The method
is illustrated through an out-of-plane elasticity problem of a membrane with a microstruc-

ture of a spatially random inclusion-matrix under a deterministic load. The key concept

ii

introduced here is a random meso scale continuum model. It is found that two bounds on
the material properties and in turn on the global response have to be considered in the

analysis.

In the last problem, the effective thermal conductivity of functionally graded heteroge-
neous interphases between fiber and matrix in such composites is determined. The topol-
ogy of microstructure is taken as a mosaic or a random chessboard where both phases
have locally isotropic properties. The resulting meso-continuum model of the interphase
is used to calculate the effective macroscopic properties (transverse conductivity or,
equivalently, axial shear modulus) of such composite materials. This problem requires the
treatment of several length scales; the fine interphase microstructure, its meso-continuum
representation, the fiber size and the macroscale level at which the effective properties are

defined.

iii

DEDICATION

In memory of my father

iv

ACKNOWLEDGMENTS

I wish to acknowledge all who assisted me in this undertaking and completion of this
dissertation. I am indebted to Dr. Martin Ostoja-Starzewski, my adviser, for his valuable
time, assistance and encouragement. His kind consideration and understanding have been

an incentive for the completion of this dissertation.

Sincere gratitude is extended to the other members of my guidance committee Dr.
Ivona Jasiuk, Dr. Michael Thorpe, and Dr. Dashin Liu for their contribution, helpful dis-
cussions, and criticism. Also, I would like to thank Wei Wang who provided me with

some of the figures of chapter five.

Finally, I would like to thank my mother and pray for my deceased father who raised
me and devoted their lives to provide me with happiness. Also, my heartfelt appreciation
is extended to all the members of my family for their constant encouragement and their
moral support. My special gratitude goes to my wife Suha for her patience, care, and
encouragement in each step I made towards the completion of this research. I do not for—
get my son Osama whose smiles and childish looks gave me the courage to accomplish

my goals.

TABLE OF CONTENTS

SUBJECT ...................................................................................................................... page
LIST OF FIGURES ......................................................................................................... viii
CHAPTER 1 INTRODUCTION ..................................................................................... 1
CHAPTER 2 BACKGROUND ....................................................................................... 6
2.1 Random Field Concept ........................................................................................ 6
2.2 Definition of Effective Properties ...................................................................... 11

CHAPTER 3 MICROMECHANICALLY—BASED STOCHASTIC FINITE

ELEMENTS ............................................................................................. 16
3.1 Introduction ....................................................................................................... 16
3.2 The Finite Element Formulation ........................................................................ 23
3.3 The Micromechanics Basis ............................................................................... 27
3.4 Methodologies of Solutions .............................................................................. 37
3.4.1 Exact Calculation Method ....................................................................... 37
3.4.2 Covariance-Based Method ..................................................................... 40
3.5 Numerical Results and Discussion .................................................................... 42
CHAPTER 4 EFFECTIVE ELASTIC MODULI OF DELAUNAY NETWORKS ....... 53
4.1 Introduction ....................................................................................................... 53
4.2 Problem Statement ............................................................................................ 57
4.3 Models of Random Microstructures .................................................................. 61
4.3.1 Regular Triangular Networks (Lattices) ................................................. 61

vi

TABLE OF CONTENTS (CONT)

SUBJECT ....................................................................................................................... page
4.3.2 Dual Voronoi-Delaunay Networks ......................................................... 63
4.3.3 Central Spring Model ............................................................................. 69

4.3.4 General Spring Model (Combined Central and Angular Spring

Models) .................................................................................................. 72

4.4 Methodology of Solution ................................................................................... 75
4.4.1 Spider-Inclusion Analogy ...................................................................... 75
4.4.2 Periodic Boundary Conditions ............................................................... 78
4.4.3 Numerical Solution Implementation ...................................................... 81

4.5 Numerical Results .............................................................................................. 83
4.5.1 Central Interaction Model ..................................................................... 83
4.5.2 Combined Central and Angular Interactions Models ............................ 92

CHAPTER 5 HEREROGENOUS INTERFACES IN FIBER-MATRIX

COMPOSITES ......................................................................................... 99

5.1 Introduction ........................................................................................................ 99

5.2 Problem Statement ........................................................................................... 101

5.3 Solution Procedure ........................................................................................... 104

5.3.1 Passage from the Micro to Meso Level ................................................ 104

5.3.2 Passage from the Meso to Macro Level ............................................... 107

5.4 Discussion of Results ....................................................................................... 111

CHAPTER 6 CONCLUSIONS .................................................................................... 118

BIBLIOGRAPHY ............................................................................................................ 1 23
vii

LIST OF FIGURES

FIGURE ......................................................................................................................... page
2.1- Realizations of 2-D fiber-matrix composite taken from its sample space ................... 8
3.1- Triangular finite element with its three out-of-plane degrees of freedom .................. 24

3.2- A window of size 5 = L/d in (a) Voronoi microstructure (b) matrix-inclusion
composite .................................................................................................................. 28

-r
3.3- A schematic 8-dependence of upper and lower bounds C ; and (5;) ,
respectively. ............................................................................................................... 34

3.4- Microstructure of a finite element discretized by a fine finite difference scheme (a)
one grain cell 1 with spring model (b). ...................................................................... 39

3.5 - A membrane of a matrix-inclusion composite with a finite element mesh of
resolution 8 ................................................................................................................ 43

3.6- Ensemble average “upper” and “lower” responses, normalized over the deterministic

case Ci = Cm = 1. 5“) and 50!) correspond to C; and C; , respectively ................. 44

3.7- Graph of standard deviation of the “upper” and “lower” responses ........................... 45

3.8- The domain, boundary conditions, and finite element discretization of the
domain ........................................................................................................................ 48

3.9- Bounds on response (reference temperature of comer point) as a function

of 5 ............................................................................................................................. 49
3.10- Standard deviation of the upper bound (for reference temperature) as a function
of 5. ......................................................................................................................... 51
311- CPU time of FE and MM calculations plotted as a function of 5. ........................... 52
4.1- Random assemblage of particles from granular material with spring
interaction model. ...................................................................................................... 58
4.2- Triangular (regular) network of size 7x7 with a hexagonal unit cell. ........................ 62
4.3 - Unit cell with six springs modeling the stiffness of the original continuum. ............ 63
viii

LIST OF FIGURES (CONT)

FIGURE ......................................................................................................................... page

4.4- Formation process of a Voronoi graph. Poisson points (a) expand into disks at the
same rate (b). No disk may impinge on the area of another (c), hence, they evolve

into area-filling polygons (d). .................................................................................... 65
4.5- A Voronoi tessellation in 2-D plane (200 random Poisson points in
a unit square) ............................................................................................................. 66
4.6 - Delaunay network dual to the Voronoi tessellation graph of Fig. 4.5
(200 random Poisson points in a unit square) ........................................................... 67
4.7- Central spring model for an in-plane elasticity problem (a), Voronoi unit cell
with m springs incident into it (b) ............................................................................. 70
4.8 - Combined central and angular spring model. ........................................................... 73

4.9 - A Voronoi tessellation graph and its dual Delaunay network (a),
“spider-inclusion” analogy (b). ............................................................................... 77

4.10 - A periodic Delaunay network generated by duplication of Poisson points
of one square into other 9 squares ........................................................................... 80

4.11- Dependence of the Poisson’s ratio v on p for a one-phase Delaunay network

with all spring constants assigned according to k = 1". ........................................... 84
4.12- Effective bulk modulus at contrast 10 for Delaunay and regular networks;

along with various self-consistent models and Voigt and Reuss bounds. ............... 87
4.13- Effective shear modulus at contrast 10 for Delaunay and regular networks;

along with various self-consistent models and Voigt and Reuss bounds. ............... 87
4.14- Effective bulk modulus of Delaunay and regular networks at contrast 100 ............. 88
4.15- Effective bulk modulus of Delaunay and regular networks at contrast _

1,000 ......................................................................................................................... 88
4.16- Effective bulk modulus of Delaunay and regular networks at contrast

10,000 ....................................................................................................................... 89
4.17- The reciprocal of effective bulk modulus at contrast 100. ....................................... 90
4.18- The reciprocal of effective bulk modulus at contrast 1000. ..................................... 91
4.19- The reciprocal of effective bulk modulus at contrast 10,000. .................................. 91

ix

LIST OF FIGURES (CONT)

FIGURE ......................................................................................................................... page

4.20- Effective bulk modulus of Delaunay network at contrast 10 for different
values of r. ................................................................................................................ 95

4.21- Effective shear modulus of Delaunay network at contrast 10 for different
values of r. ................................................................................................................ 95

4.22- Effective Young’s modulus of Delaunay network at contrast 10 for different
values of r. ............................................................................................................. 96

4.23- Effective Poisson’s ratio of Delaunay network at contrast 10 for different
values of r. ................................................................................................................ 96

4.24- Effective bulk modulus from numerical simulations and self-consistent models

results for r = 10 ....................................................................................................... 97
4.25- Effective bulk modulus from numerical simulations and self-consistent models
results for r = 50. .................................................................................................... 97
4.26- Effective bulk modulus from numerical simulations and self-consistent models
results for r = 100. .................................................................................................... 98
5.1- Sketch of a transverse cross-section of a unidirectional fiber-reinforced composite
with heterogeneous interphase. ................................................................................ 102
5.2- Sketch of a fiber-interphase-matrix system. and two models of
microstructure .......................................................................................................... 1 05
5.3- Sketch of Composite Cylinders Assemblage model ................................................. 107
54- Variation of conductivity in anisotropic interphase under essential and
natural boundary conditions for 5 = 10. ............................................................... 114
5.5- Variation of conductivity in anisotropic interphase under essential and
natural boundary conditions for 5 = 20. ................................................................ 115
5.6- Variation of conductivity in anisotropic interphase under essential and
natural boundary conditions for 5 = 40. ................................................................ 116
5.7- Conductivity of composite with anisotropic interphases vs. volume
fraction f of inclusions for 5 = 20 ........................................................................... 117
X

 

CHAPTER 1

INTRODUCTION

Solution of any engineering problem involves two types of information: input and out-
put. In mechanical problems the input includes specifying a system (material properties),
a loading (source of disturbance) and boundary conditions on the system. Output, which
depends on the input, is nothing but the response of the system under the specified distur-
bance and boundary conditions. In spite of high quality of control processes in manufac-
turing of engineering parts, a fully deterministic description of their properties is hard or
even impossible. Also, because it is difficult to measure (or describe) disturbances and
boundary conditions, the exact values of these quantities are usually unknown. Therefore,
in an analysis one must admit uncertain nature of the problem which necessitates the use
of stochastic treatments rather than deterministic ones. Stochastic (probabilistic) analysis
in its broadest sense refers to the explicit treatment of uncertainty in any quantity entering
the corresponding deterministic analysis. Difference in the Young’s modulus along a
beam axis, eccentricity of mass density on a disc of a turbomachinary rotor, a change in a
boundary shape or condition, wind-induced forces acting on a structure, etc. are all exam-
ples of randomness in typical structural engineering problems. In this dissertation, sto-
chastic treatment is devoted only to variability in systems (i.e. material properties), while

randomness entering in other parameters is beyond its scope.

It is needless to mention that the necessity for such a stochastic treatment is highly
needed in particular for heterogenous materials. Composites, a special form of heteroge-
nous materials, consist of two or more constituents (phases) which may be perfectly

bonded together. Fiber-reinforced plastics are an example of two-phase materials, while

 

polycrystals are another example of a heterogeneous solid with an infinite number of
phases where each crystal orientation is associated with one phase. Composite materials
are classified according to different factors such as geometry and distribution of phases
and type of bonding between these phases upon which overall properties are determined.
Due to uncertainty in manufacturing and/or variability (randomness) in its constituents,
the composite properties vary at micro-scale from one point to another and possess a local
directional dependence. Such properties are termed as locally inhomogeneous and aniso-
tropic, respectively. Thus, for random materials, calculation of effective properties at the

macro level should be extracted based on their microstructures.

In recent decades, the subject of effective pr0perties calculations has become of inter-
est to many researchers. Many investigations are in progress, but the subject is not well
established yet where more studies are needed. The sub-branch of mechanics which deals
with aspects at the level of microstructures is called “micromechanics”. Several analyti-
cal micromechanics-based methods to compute the effective materials properties like elas-
tic, thermal and others are available in literature. Basically, these methods can be

classified according to Hashin [1] into three categories:

1) Direct (Exact) approach: where exact calculations of effective properties for some geo-
metrical models of composite materials are performed in terms of averages. It implies that
the averaged microfield quantities satisfy the phase governing differential equations, the
phase interface conditions and the external boundary conditions on the composite. Actual
direct computation is extremely difficult due to the necessity of satisfying the interface

conditions between phases. Therefore, it must be restricted to simple models.

2) Variational approach: In a certain sense, it is more powerful than the direct approach

since it leads to bounds on the effective properties when exact calculation is not possible.

 

 

In particular, results for composites with irregular phase geometry of partial information

can be obtained only by this approach.

3) Approximations methods: they have various forms and are independent of a model or a
theory. Perhaps, the most primitive form is to assume that the property has an expression
with undetermined parameters which might be fitted by an experimental data. Within the
dilute limit, another form of approximation is to analyze a model for the a composite
material on the basis of assumptions that are in principle incorrect, with the hope that the
error introduced will be minimized. Effective medium theories may be categorized under

this approach.

In this dissertation, stochastic treatments of medium’s variability (i.e. in material prop-
erties) will be presented through several engineering problems. Considering of variability
of microstructure, we seek meso-continuum models which provide a basis for calculations
of effective macroscopic responses. Micromechanics-based methods along with various
random discrete models of microstructures are used for the passage to the meso-contin-
uum approximations. These methods which involve computer simulations and numerical
analysis will be established and used in various problems. For sake of clarity in presenta-
tion, we focus on simple problems in elasticity as well as thermal problems like two-
dimensional (2-D) conductivity. Namely, the problems that will be addressed and treated

here are:

'm 'l- ifinilmnnli'
The bulk of finite element work and its applications to date have concentrated on deter-
ministic analysis, where a single problem is modeled and solved numerically. For many
applications this is inadequate due to the non-unifonn properties of the media or in other

parameters at specific length scales. Therefore, there is a great need for stochastic tech-

 

 

niques to incorporate such non-uniformity. Thus, Stochastic Finite Element Methods
(SFEM) have recently become an active area of research. Here we attempt to combine
two crucial methodologies developed to deal with complex problems of modern engineer-
ing: finite element analysis and probabilistic analysis. Chapter three of this dissertation
is concerned with stochastic finite element analysis through solutions of elasticity and
conductivity problems. A method based on micromechanics is developed to accommo-
date any variation in material prOperty at the micro-level over the domain of finite ele-
ments. Discretization of element’s random field using a digital-image based method in a

form of finite difference scheme is used to link micro to meso properties.

i ' Ii D l n n rk :

In the context of granular media, most previous works in this area were concerned with
specific parameters or phenomena (e.g. contact forces between grains and shear flow of
granular media) [2]. Macroscopic constitutive modeling was based on primitive models
where particles were represented as round disks with one contact point between them
through which force can be transmitted, e.g. see [2]. Unlike composites, there was a lack
of effective medium theories for effective properties prediction. These observations pro-
vide a motivation for the present study. A discrete random model developed with the help
of a well-established graph model, so-called Delaunay graphs, is used to represent micro-
structures of granular media. For a two-phase granular medium modeled by a Delaunay
network, computer simulations are carried out to calculate its effective elastic properties.
Results are used to establish a special effective medium which takes into account several
parameters. Details are presented in chapter four. It is worthwhile to mention here that
such development will furnish a foundation for a further research work in this area.

_iA' H! H 01“u"10_n -.'. 10 '11-. °. 01100"

In fiber-matrix composites, random mixture of fibers and matrix grains leads to a cre-

ation of a new zone between them called “interphase”[3,4]. This zone plays an important
role in affecting the overall properties of the composite. In literature, a few analytical
studies dealing with the interface problem can be found, see for example [4—12]. Actually,
most of these studies focused on interfaces’ elastic properties. Even for thermal proper-
ties, pure approximated treatments with some specific assumptions were adapted [4].
Therefore, it is highly needed to tackle the problem with a real consideration of micro-
structure of the interface. To that end, random field modeling of the two-phase interface
microstructures is developed in order to calculate its effective thermal conductivity (axial
shear modulus). In chapter five, we develop two meso-continuum approximation of a
heterogeneous interface based also, as in the finite element problem, on a digital-image
based method and computer simulation procedure. Then, we proceed by introducing
effective conductivity results in a specific effective medium theory, Composite Cylinder

Assemblage model, to obtain the overall thermal conductivity of the whole composite.

In this dissertation, the sequence of presentation will be as follows: in Chapter two, we
discuss some of the random field concepts and how they can be applied to heterogenous
materials to calculate their effective elastic and thermal properties. Chapter three is con-
cerned with the stochastic finite element problem in which we propose methodologies for
treatment for materials’ variability over the domain of the finite element. In chapter four
we describe how numerical implementation can be used to calculate effective elastic mod-
uli of Delaunay network type of microstructures and identify proper effective medium the-
ories for two-phase granular media. Chapter five deals with the linkage between
different length scale calculations of thermal conductivity of heterogenous interfaces in
fiber-matrix composites. Finally, in Chapter six we present some conclusions and recom-

mendations for further future research work for of the three problems.

 

CHAPTER2

BACKGROUND

2.1 Random Field Concept

In elementary probability theory, random variable X((D), (0 E Q is defined as a real
function mapping a sample space 9 into the real line R. Usually, in experiments, out-
comes are described in terms of random variables as uncertain quantities to be observed.
Similarly, it is meaningful to associate with the term “random field” a giant laboratory
where an investigator can do numerous experiments. Observing the outcomes of all the
experiments in the laboratory and observing the realization of the random field are both

equivalent.

In the terminology of random field theory, the location of each experiment setup in the
laboratory is identified by a set of “coordinates” or “parameters”. In an n-dimensional
space, a vector t = (t1, t2, t3, ....... , tn) whose elements, are the coordinates or parameters
corresponds to the locations. A random field is also called a random (or stochastic) pro-
cess, where field indicates that the parameters space is multi-dimensional, while random

variable depends on a single coordinate, usually time.

A random field may be defined as a family of random variables X(t) depending on more
than one deterministic parameter. The collective outcome of all the experiments compris-
ing the random field is denoted by X(t), the realization of the random field. The outcome
of the random variable associated with any location I is referred as the state of random

field at that location. Thus, the random field is denoted as “discrete state” or

 

 

 

“continuous state” one depending on whether the range of random variables is discrete or
continuous. Depending upon the locations of the observation points, random field may be
called more specifically [13]:

1) a random series: if observations are made at discrete points on a time axis,

2) a lattice process: if observation are made at the nodes or sites of a lattice in space,

3) a continuous random function: if observations are made at all points along one or more
spatial coordinates (x1, x2, x3) and/or the time axis (t),

4) a random partition of space: if a discrete random variable is observed at every point in
space,

5) a random point process: if points (or, in reality small objects) are located in a random

pattern in space.

Of course, these different descriptions are often closely related. For example, a contin-
uous parameter random field X(t) may be sampled (observed) at a lattice of locations in
the parameter space. At these locations, one may choose to assign one of two values, 1 or
-1. This converts a continuous state field into a binary (two-phase) field. A regular point
pattern can be made to fill the entire space by assigning to each point its own polyhedron
For n = 2, the vector connecting adjacent points is bisected by a plane, and the space
enclosed by all such planes is associated with the point it contains. If the pattern of points
is random, this same construction yields so-called “Voronoi cells”. By assigning a value

to each cell, one obtains a discrete-state, continuous-parameter random field.

The foregoing discussion is about the random fields in general. Here, more specifically,
and as an application in the mechanics field, random heterogeneous materials (compos-
ites) are considered. Due to uncertainty in manufacturing, any actual composite must be
treated as a random medium. In statistical sense, the random medium B is defined as:

B = mm»; (0 6 O} which is a family of deterministic specimens 8(0)), (1) is one realiza-

 

 

tion from the sample space 0 (see Fig. 2.1).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0. 0 0 .0” 0 .0 0 .
0 0 0 0
0 0 0 0
0 0 0 . 0. . 0 .0
___I__
0. 0.0. 0 . J 0 .0
.0 0 . .0 . .0.0.0 ........ etc.
0 0 0 0 0. . . 0
__0J _I__0___

 

 

 

 

 

 

 

Fig. 2.1 Realizations of 2-D fiber-matrix composite taken from its sample space.

For simplicity and clarity, we consider a random medium of 8(0)) which is a two-
dimensional random composite of fiber-matrix type as shown by Fig. 2.1, where the ran-
domness is in the location of the fibers. Thus, any property of the random composite B,
say Young’s modulus E as an example, is regarded as a random process or random field, E
= {E((D); (t) E 9} over the domain of the specimens where E((D) describes Young’s mod-

ulus of one specific realization. The mean or ensemble average of B may be calculated as
(E) = JE(w)p(w)dco (2.1)
Q

where p((0) is the probability density function of E over the sample space 9. For a dis-
crete arrangement of phases, the probability function P, ()5) of finding phase r at x is
given by
P, (at) = (E, (15)) = IE, (r, (0);? ((0) dc» (2.2)
o

 

Here, E (1;) is the indicator function which is defined as

{I if it lies in r phase}

E (x) = . . (2.3)
r {0 If r; [res elsewhere}

Similarly, the probability of finding simultaneously phase r at )5 and phase 5 at .3’ is

P,, (2!. 15') = (15.0.6) E, (5')) (2.4)
Probabilities involving more points or phases follow the same pattern.
The auto-correlation function of Young’s modulus E ()5) is defined as
n n

(E (2!) E (20> = 2‘. 2 E,E,P,, (at. r) (2.5)

r=ls=l

where E, and Es are the Young’s moduli of the two phases r and s, respectively, and the

E ()5) is the Young’s modulus of the composite at point )5.

The material occupied by a random region of volume V is defined to be Statistically
Homogeneous (SH) by the requirement that P, (.3) be invariant with respect to any trans-

lation inside the body. In mathematical notation, it can be stated as:

V (r. r') e V (2.6)

Polycrystal is an example of a material which is statistically homogeneous. This property
is called strict stationary in the theory of stochastic (random) processes. The statically
independence (or no long-range order) property of statistical homogeneity is attained if we

have. on large volume of the region, the condition

10

Paar) =P,(;c>P,(;c') as (as-41*” (27)

This condition is of importance in the analysis of random fields where we note that it is not
satisfied for a material with periodic microstructure. On the other hand, the random com-
posite is described as being statistically isotropic, if the property over the body domain

remains the same in all directions irrespective of any rotation in space.

Suppose now that the dimensions of the region V occupied by the composite are infi-
nitely large compared with typical microstructure dimensions, or microscale, such as grain
size or mean inclusion spacing. If the specimen size is sufficiently large compared to the
size of heterogeneity, then the effective property can be defined based on the volume aver-

age over this specimen. Thus, the volume average of E over one realization is
E = I E (1;) dV (2.8)
v

For statistically homogeneous media, it is reasonable and of much help to make the
ergodic assumption which asserts that the volume average of one specific E((D) over V
equals the ensemble average of the random function E (x, 0)) at any location it. Mathe-

matically, the ergodicity condition is stated as
Eaéjeupm = <E>sjE(r.m>p(w)dw (2.9)
v Q

where E is the volume average over a fixed (D, while (E) is the ensemble average over all
specimens at fixed 3;. By making such an assumption, the statistics of Young’s modulus

(E) can be determined from a single realization E(0)).

11

2.2 Definition of Effective Properties of Heterogenous Media

Based on these stochastic (statistical) concepts of random fields theory and on a proper
model of the microstructure, any property of the random composite (or heterogenous
solid) can be defined. In the present dissertation two mechanical pr0perties are of interest;
the elastic moduli and thermal conductivity. The definition of effective elastic properties

is given here, while the analogy to the thermal conductivity will be presented at the end of

this chapter.

Consider one realization 8(0)) taken from the random body shown by Fig. 2.1 with
finite volume V and boundary 5V. Hooke’s law at any point inside the domain of the body
is
011(5) = CW (1:) 8a (2:) (2.10)
where Cijkl is a fourth-order tensor characterized by a random distribution over the body
domain. Depending on the type of boundary condition (displacement, traction or a combi-

nation of the two) applied over 8V and body forces specified over V, then the problem is

to solve the differential equation in the form of
00,].(5) +fi(g) = 0 (2.11)

to determine 9 (it) , § ()5) and y ()5) . If the moduli of (2.10) vary on the microscale, its
solution will provide some locally averaged information. But, we seek an overall consti-

tutive equation of the body at the macroscale level in the form of

9 = 9% (2.12)

12

where Qeff is the tensor of effective moduli at macroscopic level of the body. Definitions
and derivations of C” form much of the goals of Micromechanics. The determination of
such a constitutive equation may be achieved based on either the volume averaging over
the body domain, or on the ensemble averaging over the body’s sample space of speci-

mens (realizations) O.

For the sake of defining effective moduli (3”) based on volume averages, we assume

a particular specimen B(0)) to be subjected to a displacement boundary condition given by

u. = 8.x. xje 8V (2.13)

where Bil-is the uniform strain field over V if 8(0)) is a homogeneous body which means
that C is a constant tensor field. But if C is varying over the domain, and hence, g and g

tensors, then with help of Gaussian theorem the average strain may expressed as
1

from which and after solving the differential equation (2.11) we can define the effeCIive

moduli by

“’2 (2.15)

10
ll
10

where g, g is the average value of stress and strain fields, respectively over the domain

the body.

Also average value of energy density field U can be calculated as

13

(2.16)

l
- Rica-80W

By using the Gaussian theorem and since Oi“- = 0, the average energy density (2.16) can

be expressed after applying the boundary condition (2.13) as

g‘ffs (2.17)

U=-—jous gems g

___ .1.

2
from which C“ can be extracted. Thus, the moduli may be defined either from the aver-
age stress and strain or from the average energy density if the boundary condition (2.13) is

imposed.

On the other hand, the specimen B(0)) may be subjected to a traction boundary condi-

tion in the form of

t. = Unn- n. on 3V (2.18)

which generates a uniform stress field over V if the body is homogeneous. Here, 6i} is the
average stress and nj is the unit normal of the boundary. As before, if the domain is heter-

ogenous, the stresses can be averaged over the body volume as

1

Vioija’v (2.19)
Similarly, the average constitutive equation in the form of (2.15) is obtained after solving
for the average strain field from the governing differential equation. Alternatively, with
the help of equation (2.19) and after substituting the boundary condition (2.18) in (2.16), it
yields

14

U = Z—IVJ/O‘feifdv = $0,112,}. = égs‘ffg (2.20)
from which flea = ($61.le can be extracted. Thus, again the overall moduli may be defined
from the average stress and strain or from the average energy density for this type of
boundary condition. Indeed, either boundary value problem defined by (2.13) or (2.18)
yields different values for cm at some scale smaller than the dimension of the body. Also,
solution of such boundary value problems provides rigorous definitions for the overall
moduli for any particular inhomogeneous body 3(0)) and are useful in predicting other

behaviors of composite bodies.

On the other hand, elastic moduli QC” may be defined based on the ensemble average.
For a given boundary value problem with displacement or stress controlled boundary con-
ditions, we seek a solution 6 (g), e (15) , u (x) at a fixed 1; e V for each 8(0)) realization

so that, the effective constitutive law can be defined as
(o) = c"”<e> (2.21)

here ( ) is the ensemble average over the sample space 52 (recall (2.1)). Further applica-
tion of ergodicity condition leads to the ensemble averaged version of (2.21), in which the

mean energy density (U) at )5 is expressible in the form of
1 l 1
(U) = 5058) = -2-<o><8> = 5<8>Qeff<8> (2.22)

from which QC“ can be extracted. Hill conditions are assumed in order to obtain 2.22
which can be stated as:

i) the body B is ergodic,

ii) the number of grains --> infinity,
iii) distribution of volume stress sources, if present. are not correlated to the local elastic

moduli.

The forgoing discussion which focused on elastic properties is still applied for any other
property, say thermal conductivity, with the following analogous conceptual relations in

terminology [ l]

 

 

 

 

 

 

 

 

Elasticity Conductivity
u,- - displacement T- temperature
8 ij - strain T J - temperature gradient
O’ij - stress qi - flux
ti - traction qn- normal flux component
Cijkl' elastic moduli Ki} - conductivity
Sim - compliance Rij - resistivity

 

 

 

 

The mathematics of the conductivity problem is considerably easier than that of elastic-
ity since vectors take place of second rank tensors. In addition, scalar Laplace equation

takes place of vectorial elasticity Naiver displacement equation.

 

CHAPTER 3

MICROMECHANICALLY-BASED STOCHASTIC FINITE
ELEMENTS

3.1 Introduction

Stochastic finite element methods (SFEM) in its broadest sense refer to the explicit
treatment of uncertainty in any quantity entering the corresponding deterministic analysis.
In such an analysis randomness may arise from three sources: distortion of the boundary,
fluctuation in the external force fields, and heterogeneity of the system (medium). SFEM
have begun to be used in solving mechanical problems of system stochasticity only in the
last decade. Researchers in this field have attempted to combine two types of analysis
developed to deal with complicated problems: finite element analysis and stochastic anal-
ysis. However, solving such problems, which involve continuous random variation in
medium’s properties, is not an easy task even for simple cases. Analytical solutions in
SFEM are available only in few simple cases, generally limited to one-dimensional prob-
lems [17]. Useful analytical tools for performing any structural analysis involving uncer-

tain properties are provided by the theory of random fields which is a branch of the

probability theory.

Several approaches for establishing SFEM are known. Theoretical background, advan-
tages and disadvantages of various stochastic finite element methods can be found in a
comprehensive survey conducted by Brenner [14]. In this survey, the author discusses
various representations of stochastic fields such as the point discretization, the local aver-

aging, the interpolation, and series expansion methods. In literature, most of stochastic

l6

l7

finite element solutions were performed based on one of these different techniques of ran-
dom field representations. Perhaps the perturbation approach, a series expansion-based
method, is the most common one which was developed by many researchers. The method
is based on expressing random functions as sums of deterministic and random compo-
nents, which are then substituted in the descretized equations transforming the random
system equations to a theoretically infinite number of deterministic equations. Depending
on the number of terms considered in the expansion, on (or two), the method is termed as
first (or second)-order approximation. First-order expansion form was utilized, first time,
by Cambou [15] in seeking finite element solutions of linear static problems with system
stochasticity. For solution of similar problems, Hisada & Nakagiri [16] utilized a second-
order perturbation approximation. In another simple case, Vanmarcke [17] developed a
stochastic finite element procedure for solving problems in which the physical properties
exhibit one-dimensional spatial variation. The procedure was based on the second-order
statistics of the random parameter where randomness was represented by a local spatial
variation over the element and characterized by three parameters; the mean, the standard
deviation and the scale of fluctuations. For, multi-degree-of-freedom structures, Shinozo-
uka and Yamazaki [18] performed a static finite element analysis where Young’s modulus
or Poisson ratio were assumed to vary randomly according to Gaussian (or non-Gaussian)
distributions. Results with different stochastic finite element techniques were used; first
and second order perturbation techniques, Neumann expansion, an expansion Monte
Carlo simulation and direct Monte Carlo simulation methods were obtained. These results
were found to be very close for small values of a coefficient of variation for the random

parameter.

Linear partial derivative is also a particular form of the perturbation method where the
first and second order correction terms are expressed as partial derivative of the response

in terms of a defined random parameter. Baecher and Ingra [19] used this form in solving

 

 

18

geomechanical problems to calculate structural settlement in soils with stationary random
properties. The covariance matrix of the total settlement was expressed in terms of partial
derivatives of the displacement vector with respect to element moduli. For a comprehen-
sive review of perturbation techniques-based stochastic finite element methods the reader
may refer to a survey conducted by Benaroya and Rehak [20]. The major contribution of
their review was in illustrating the similarities and differences between various methods
based on the first and second-order moment calculations. They concluded that all pertur-
bation methods rely on the same theoretical foundations with only differences in the

numerical implementation.

On the other hand, many researchers concluded that due to the discrete nature of a finite
element procedure, stochastic functions of spatial variation must be descretized. Then,
choosing a finite element size depends on the statistical correlation of the random field -
mesh should be fine enough for accuracy purposes but large with respect to the correlation
length. Der Kiureghian [21] proposed to use two meshes; the first one for finite element
computational efficiency based on structural topology, the other for random field discreti-
zation on the correlation length scale. In this regard, discretization of the random field can
be found in Vanmarcke [17] where he defined a scale of fluctuations and suggested it to be
used in the finite element analysis. Liu, et al [22] discussed the discretization of a random
field issue in stochastic finite element analysis. Based on a second moment method, they
transformed the correlated random variables to uncorrrelated ones which leads to a reduc-
tion in the computational time of analysis. Such treatment was presented for linear and

nonlinear continua with inhomegenous random fields.

As was mentioned earlier, the random global stiffness matrix of a structure may be split

into two parts, a deterministic part (l_(0) and a random one US.) as follows

 

19

N‘
K = 21g“) = sour (3.1)
e = 1
To solve the static problem
(Ko'l'K’k)! = E (32)

an inverse of the stiffness matrix is needed, so several methods were proposed to deal with
this inversion problem. In the weighted integral method, proposed by Deodatis [23] and
others, the elements of fluctuating part of stiffness matrix are assumed to be random vari-

ables and expressed in an integral form. Then, these can be decomposed into weighted
(e)

integrals Xi“) and the associated coefficient matrices AK i .

Neumann series expansion method has been treated by many researchers (e.g. Shino-
zoka, Yamazaki , Spanos, Ghanem and Deodatis). Application of the perturbation tech-
nique in stochastic finite element analysis yields a stiffness matrix expanded in the form of
Taylor series. Responses are also required to be expanded in a Taylor series expansion
form. The homogeneous chaos expansion method was applied to stochastic finite ele-
ments by Ghanem and Spanos [24] where they used a Karhunen-Loéve expansion to
replace random variables involved in the stiffness matrix. Other methods to handle sto-
chasticity in finite element analysis can be found. Among these is a method in which the
probability density of the response is expressed in the form of a Fokker-Planck equation.
Thus, there is only a need to solve a deterministic equation using a deterministic finite ele-
ment method. Also, the Monte Carlo simulation method consists of solving deterministic
governing equations for some chosen realizations of the random parameter. Point discret-
ization-based methods which use successive iteration simulations are useful in estimating

the first and second moment of responses. Despite the fact that the method is less accurate

20

than the Monte Carlo method, it is applicable to a more complicated class of problems.
Although reliability methods are complicated, they have received an attention of numer-
ous researchers recently. More information can be found in a paper by Vanmarcke et a1
[25] in which the authors conducted a survey on the development of stochastic finite ele-
ment analysis. Discussion about several methods such as digital simulation, Monte Carlo
simulation, perturbation techniques and reliability analysis was presented. Their work
contained a discretization of the random parameter (vector) for a material property whose
covariance matrix depends on the finite element mesh. Analytical tools to convey first and
second order information about homogenous random fields were described. Since statisti-
cal information such as spatial correlation about materials properties was difficult to

obtain, it was assumed to be the form

R (Ax) = [W (3.3)

where Ax is the distance between two points, and 7k is an adjustable constant. Similarly,
the uncertainty in a parameter, 2 in the stiffness matrix, was represented by a small random

variable defined as

z=2(1+or) (3.4)

in which Otis a variation parameter and 2 is the corresponding deterministic value.

The computational time of a stochastic finite element procedure (CPU-time), which
determines its efficiency and practicality, has become an important aspect. As an exam-
ple, the computational time of the direct Monte Carlo method, which has proven to be the
most powerful approach, severely limits its applicability. In a deterministic finite element

solution, the CPU time depends on the procedure of taking the inverse of the stiffness

 

 

21

matrix. Therefore, in a more complicated problem of system stochasticity, any solution
procedure should concern the CPU time issue. Many researchers in developing their pro-

cedures, considered such an aspect and reported data about it (e. g. see [25]).

Thus, the different methods used in incorporating the random effects in material proper-

ties can be summarized and listed as:

1) Perturbation methed: where the random system is replaced by a theoretically infinite
number of identical deterministic systems each of which depends on solution of lower

order equations. Up to second-order, the response in a static problem is expressed as
2
{U} = {U0}+e{U]}+e {U2} (3.5)

where 8 « 1

2) WM: which, as mentioned earlier is based on a Neumann series

for the inverse of a random operator [K((r))], takes the following form

[k(w)] = [1—P(m)+P2(m)—-P3(m) +...)[(K)]“
Pa») = 1<K)1"1k(w)1

(3.6)

This method reduced CPU time when stochastic problem is solved by a Monte Carlo sim-

ulation.

3) Weighted integxal method: This method focuses on the determination of the random
stiffness matrix [k(0))]. The idea, which was introduced in [26] and presented through an
elastic plate problem, was to start with a locally isotropic random field, say, Young’s mod-

ulus, and assign it to all the finite elements according to

22

iEULw) = ‘<E>[1+‘f(r.w)] <‘f(sco)> = o (3.7)
then, all the element of [K(0))] are calculated as
[iK (w)] = ['20] +iX0(w) [A‘KO] (3.8)

' i . . . . . l . . .
where [‘Ko] and [A K0] are deterrnrnrstrc matrices, while X015 a random vanable grven

as

ix0 (w) = (91., co) d’A (3.9)
I

4) Spectral methed: The coefficients of expansion become correlated in the process of rep-
resenting the random function by a Fourier series. In order to retain the uncorrelateness
while obtaining the desired orthogonality of random coefficients, Karhunen-Loéve expan-

sion is used for that purpose.

It is worthwhile to mention that all above works may be considered as reasonable in
problems involving only small variations in material properties. However, all of them
lack a correct formulation of medium’s variability, which is due to purely hypothetical
assumptions concerning the random field description of a material. That is in the area of
linear elastic structures, all studies relied on a stochastic interpretation of the locally iso-

tropic constitutive law

oi]. = A(£,w)5;j€kk+211(&w)€;j (3.10)

where the Lame constants 71., [.1 are taken to be random fields. Particularly, Young’s mod-

ulus was assumed to be a random field of Gaussian type, with the Poisson’s ratio taken as

23

constant. On the other hand, in recent treatments of the SFE problem, Ostoja-Starzewski
[27, 28] suggested a more correct approach; variation in any material property should be
considered at micro level and he proposed a micromechanically-based stochastic finite
difference method. Accordingly, in this dissertation, a similar treatment of medium’s vari-
ability in stochastic finite element [29] will be introduced. The treatment takes into
account the real nature of medium’s properties without any hypothetical assumption.
Thus, we develop, for the first time, a micromechanically-based stochastic finite element
method to handle the variability in any material property over each finite element. A ran-
dom field continuum model based on micromechanics will be presented and used as input
of a specific boundary value problem, say, stochastic finite element method. Due to the
discreteness of materials on the micro level, their non-deterministic constitutive response

has to be considered at a meso level (i.e. scale of a finite element).
3.2 The Finite Element Formulation

To briefly and clearly present our micromechanically-based method, the basic steps in
setting up the finite element equations for out-of-plane displacements of a membrane are
outlined. However, this formulation is the same for all other problems governed by the

same differential equation in the form of

Bu __
(CU-(2r. (0)3) — f(2.c) (3.11)

J

3’.
8x,-

in which,

Q ()5, 0)) : elastic moduli of the membrane 2-D second rank tensor. As will explained in
next section, C is taken as a random field of sample space S).

u: out-of-plane displacement,

f ()5) : a deterministic uniformly distributed load.

24

0): one realization taken from the sample space 0.

Examples of these problems are: in-plane thermal (electrical, dielectric) conductivity, tor-
sion of rods, diffusion, permeability and any other potential field problem. On the other
hand, the same methodology is applicable to a more general class of problems in which
materials have elastic as well as inelastic microstructures, e.g. plane stress, plane strain,

and three-dimensional problems.

A right-sided triangular finite element, shown in Fig. 3.1, is considered in the present
formulation. However, the key point is that the micromechanical meso-cell which will be
introduced in the next section and called “window” corresponds to one (or more) finite

element cell(s).

u3

x2

X1

x3

 

 

U1 . U2

Fig. 3.1 Triangular finite element with its three out-of—plane degrees of freedom.
Now, the nodal out-of-plane displacement vector for the i-th element is
{‘UUDH = [u1 u2 1.13] (3.12)

While, the displacement at any point inside the element is given in terms of the nodal dis-

placement vector as

25

{u(2c,w)} = [111, ~,N,] {‘u(w)} (3.13)
where [N] is the shape functions matrix given by

N. = %(ai+bixl +cix2) where i = 1,2,3 (3.14)

l

in which A is the element area and ai, bi, c, are coefficients defined in terms of the nodal
coordinates. Thus, the displacement inside the element is taken to be linearly varying in

terms of the element’s nodal displacement values. Its gradient vector may expressed as

gr{u(r.w)} = 181(‘u(co)} (3.15)
where 5;“ (1." (1)) b b b
gr{u(2c.w)} = x .181 = [1 2 3] (316)
a .
511050)) Cr 62 C3

in which [B] is the gradient matrix. The stiffness matrix of each finite element can be cal-

culated from

[‘K(w)] =]181T[‘C(w)]181dA
4 (3.17)

= 181T["6(w)]131A
where PC ((0)) is stiffness tensor of the i-th element of one realization of the random
medium. Since we are dealing with heterogenous media, [C] possesses a variation in its
material properties at micro level. Therefore, it is very difficult (or even impossible) in
some cases to be obtained by the conventional methods, say from experimental measure-

ments. Thus, more sophisticated techniques are highly needed in this regard. Random

26

field description of the material at micro level and Micromechanics-based method are the
basis of [C] determination (see next section for more details). This approach revealed that
a unique deterministic value of [C] should not be assumed. Instead it can be bounded in a
statistical sense by two values C; or 92 which are upper and lower bounds on the
actual [C], where 5 is a parameter related to the size of the element and e, n stand for
essential and natural boundary conditions, respectively. In a case of a homogeneous iso-
tropic medium, which corresponds to the deterministic case, a unique value of [C] is given
by

8

£5 = Q; = k[; (1)] ,k=constant (3.18)

At this stage, we synthesize the global stiffness matrix from those of the sub-elements

over the entire region as

[K(wn" = 2E‘K(w)1€ or 81mm)" = 2 [‘K(w)]' (3.19)

i=1 i=l

depending on which (I; or Q; is used. Here, (ne) is the total number of finite elements

into which the domain is descretized.

Finally, after applying the specified boundary conditions in the problem, the unknown
quantities (responses) are obtained by solving a system of linear algebraic equations (say,

by Gaussian elimination method) in the form of
81K(w)1’{U(w)1’ = {F} g1K(w)1"1U(w)1" =1F} (3.20)

Here, [U((D)] is the global response vector which represents the random out-of-plane dis-

placements in case of the membrane problem. While Ue provides a stiffer statistical solu-

27

tion and Un a softer one, the actual solution U lies between these two bounds. The

effective global force vector {F}, is assumed to be deterministic in this problem.

3.3 The Micromechanical Basis

The random medium (random microstructure) concept, as already discussed in chapter
2, plays a fundamental role in this micromechanical formulation. Commonly, in stochastic
mechanics random medium is taken as B = (8(0)); 0) E Q}, i.e. a set of deterministic
media, where 0) indicates one specimen and Q is the underlying sample (probability)
space. Formally, Q is equipped with a G-algebra F and a probability distribrrtion P. In an
experimental setting, P may be specified by a set of stereological measurements, while in a
theoretical setting P is usually specified by a chosen model of a microstructure. For certain
types of materials, e.g. a polycrystalline material, a Voronoi tessellation provides a simple
model (Fig.3.2-a), while in case of a matrix-inclusion composite, 8(0)) may be taken as a
planar Poisson point process in which overlapping disks are excluded (Fig. 3.2-b). In case
of a Voronoi tessellation, each cell, centered at a Poison’s point 3;, is assumed to be occu-
pied by a homogeneous continuum governed by a stiffness tensor Q (g, 0)) following a
space-homogeneous probability distribution P(Q). The disks are assumed to be occupied
by a homogeneous isotropic continuum of one kind, while the matrix by a continuum of
another. Any microstructure used here is assumed to be linear elastic, and the statistics of
micro scale properties are space-homogeneous and ergodic. For simplicity and clarity of
presentation, the discussion is conducted in two-dimensions (2-D); a generalization to 3-D

is quite straightforward.

 

 

 

 

 

 

 

 

(a)
Fig. 3.2 A window of size 8 = L/d in (a) Voronoi microstructure (b) matrix-inclusion com-

posite.

In both models, all the phases are assumed to satisfy the so-called ellipticity conditions:

But, [3 > 0 such that for any g the following inequalities hold
are S §§§ S Beg (3.21)
Thus, we have two realistic ergodic media models described by random fields
C = {magmas B;coe 0} (3.22)

with piecewise-constant realizations. This latter property is an obstacle to employing the

governing equations of continuum elasticity, which require that the stiffness fields be dif-

29

ferentiable. Thus, there is a need for another continuum model - one that possibly looses
some information due to a “smearing-out” procedure, but is sufficiently differentiable and

grasps the meso-level behavior.

The window concept is introduced as a square-shaped domain taken from the random
medium 8(0)), as shown in Fig. 3.2. It is important here to note that the window is equiv-

alent to one (or more) finite element cell(s). The window size is introduced as

(3.23)

0’1
11
Qll‘

which defines a non-dimensional parameter 5 typically greater than 1, specifying the scale
L of observation (or measurement) relative to a typical microscale d (i.e. grain size) of the
material microstructure. The window may be placed arbitrarily in the domain of 8(0)) with
smallest size 5 = 1; i.e. scale of crystal or inclusion. In view of the fact that either one of
the microstructures in Fig. 3.2 is a result of a certain random point process in plane, the win-
dow bounds a random microstructure 85 = (85(0)); 0) 6 Q l, where 85(0)) is a single win-
dow realization from a given specimen 8(0)). In order to define the effective moduli of 85,
we postulate the existence of an effective homogeneous continuum 8°°m(0)) of the same
volume V (i.e. area in 2-D) equivalent to 85. This equivalence is established by equating
the potential energy U (or complementary energy U") of the effective medium to that of a
microstructure 85(0)) under same boundary conditions. The boundary 885 is subjected to
two types of boundary conditions (analogous to those given by (2.13) and (2.18)):

1) Deformation-controlled (essential) boundary conditions

u = 8.01. on 885 (3.24)

I l

U . . .
where e rs a grven constant tensor and 885 15 the boundary of 85.

30

2) Stress-controlled (natural) boundary conditions
I = Gin. on 885 (3.25)

where 90 is a given constant tensor and ni is the unit normal of boundary. Thus, the essen-
tial and natural boundary conditions define two different inhomogeneous tensor fields at the
scale 8 with continuous realizations, which lead to two random meso-continuum approxi-
mations: B‘s = {8°5(0)); 0) e 52 },or B”5 = [8"5(0)); 0) e Q }, respectively. Since for any 8
> 1, 85 is a random rather than a deterministic medium, the window of size 8 plays the role
of a Representative Volume Element (RVE) of these two continuous random medium mod-
els. Then two boundary value problems should be solved in a continuum setting under the
boundary conditions (3.24) and (3.25), respectively. From solution of the boundary value
problem with condition (3.24) for potential energy U (B 5) ,the effective stiffness tensor of

any specific body 35(0)) is obtained as
Q; = Q; (co) (3.26)

superscript e indicates essential boundary conditions. Based upon such determination, the

continuum type constitutive law can be found as
s (w) = 921mg" (327)

which points to a random nature of the resulting stress field; overbar indicates a volume
averaging (i.e area averaging in 2-D). It has to be pointed out that the surface traction is
random inhomogeneous on 385(0)) , with the fluctuations disappearing in the limit as
8 —) co. On the other hand, the solution of the boundary value problem for complementary

energy U* (85) under condition (3.25) results in a following effective compliance tensor

 

31

$2 = .5200) (3.28)

in which n indicates natural boundary condition. Upon averaging over the strain field

another constitutive law is obtained in the form of
e (co) = s; (w) 9 (3.29)

which also points to a random nature of the resulting strain field, and the presence of ran-
dom fluctuations in the displacements u, on the window boundary. Due to the heterogeneity

of the microstructure 85(0)), the effective stiffness tensor given by
-1
a; (w) = [3; (1.1)] (3.30)

is, in general, different from Q; (0)) obtained under essential condition for any given finite
8. Also, both these moduli are, generally, anisotropic with the nature of anisotropy being
dependent on chosen 85(0)). The scatter in £22 and C: is strongest on the scale 8 = l, and
it tends to 0 as 8 —> oo . Since we have two random tensor fields CE and g; at our dis-
posal, this continuum approximation is nonunique. These conclusions indicate that a
locally isotropic unique random tensor field should not be assumed. Principles of minimum
potential and complementary energies can be used to show that the two values satisfy the

inequality
cguo) segue) (3.31)

The present definition of two inhomogeneous tensor fields given by (3.26) and (3.28)

are analogous to a moving locally averaged random field. Although it is conceptually sim-

32

ilar to the procedure of local averaging in the theory of random fields [13] applied to a sin-
gle realization Q (0)) ;(u e 9, but it is not the same. It becomes the same in case of a 1-D
model only when applied to compliance. In two and three dimensions computational
mechanics, methods have to be implemented in a Monte Carlo sense to find the energies
and, hence, the effective moduli of finite windows and their probability distributions
Pkgé) and Pig”. Accordingly, the autocorrelation (autocovariance) functions of Cij
were found to be 8-dependent [17] in contradistinction to the procedures employed by oth-
ers - the “weighted integral method” [26], or the “spectral representation method” [24] -

which have no connection to the material microstructure.

. . . . . . e
In VICW of the spatial homogeneity of microstructure’s statrstrcs, Q5 and Q; converge
as 8 tends to infinity. This defines a detemrinistic continuum 8del for a single specimen

8(0))
c“‘"(co) = 91(0)) = (21(0)) (3.32)

whereby the window of an infinite extent plays the role of an RVE of deterministic elastic-
ity theory. In other words, it is at 8 --) oo that the uniqueness of the constitutive law is
obtained. The assumption of ergodicity of the microstructure implies that the effective

moduli of the medium can be determined from one realization as
gd"((o) = g‘” Vcoe Q (3.33)

where C“ denotes the effective stiffness tensor of a deterministic continuum correspond-
ing to the scale 8 -—> oo . The characteristics of isotropy and homogeneity of Qefi depend on
the given medium and may be ensured by a choice of a proper model. Here, (3.26 ) and

(3.28) are isotropic due to the spatial homogeneity and isotropy of the underlying Poisson

 

33

point process and the spatial homogeneity of Pg).

Also principles of minimum potential and complementary energies can be used to obtain

a hierarchy of scale dependent bounds on the effective stiffness tensor chr

-1
(ck) - <5?) s (3;) s (.55.) s cc”: (cg) s (cg) s (Q?) a c w > s
(3.34)
This is equivalent, by inversion, to a hierarchy of bounds on the effective compliance §eff

= (geffyl

s” s (5?) 2 (53;) 2 (3;) 2 SW2 ((2 2 (93" a (6)" vs > a (3.33)
in which ( ) denotes ensemble averaging, and C R , C V are the Reuss and Voigt bounds,
respectively. These two bounds correspond to calculations of effective moduli of 8(0)) at
size 8 = 1 involving volume fraction information only. The derivations of such bounds are
based on solving two boundary value problems for minimum energy under specific types
of boundary conditions. Mathematically, Reuss and Voigt bounds correspond to arithmetic

and harmonic averages on C“, respectively.

To present graphically the discussed behavior, the eigenvalues of Q for both cases (essen-

tial and natural conditions) are plotted as a function of 8 as shown in Fig. 3.3.

34

CV

Eigenvalue
of Q 8
(c5)
1

 

cefi= c... =c"...

QR

 

 

 

 

0017

II
E

-1
Fig. 3.3 A schematic 8-dependence of upper and lower bounds (Q2) and (5;) , respec-

tively.

In the foregoing discussion, three measuring levels were introduced. Explicitly, they are:
microscale 8 = I, mesoscale 8mm, and macroscale 8M where 8M = macroscopic dimension

of the body 8(0)). These measuring scales satisfy the inequality

micro < 8mesa < 8M (3'36)
with the difference in behavior of Q (elastic moduli) at each level may be summarized as:
i) Microscale (8mm): represents the scale at the heterogeneity level (8 = 1). The difference
between C; and C; values is maximum which corresponds to Voigt and Reuss bounds

with highest scatter.

 

35

ii) Mesoscale (8mm): an arbitrary chosen scale for window (finite element cell) at which
we must deal with two values of Q (generally locally anisotropic) with some finite fluctua-
dons

iii) Macroscale ( 8M): corresponds to the dimension of the body, at which both bounds con-
verge to a unique deterministic value with scatter being converge to zero. If 8meso is on
order of 8M, one is forced to deal with spatial fluctuations at the macro scale. Thus, a sta-

tistical rather than a deterministic continuum approximation is applicable.

Since two different random anisotropic continua result, a given boundary value problem
must be solved under both types of boundary conditions (essential and natural) which
yields two global responses. In stochastic finite element analysis, the fact that the element
has a finite size implies that the fluctuations, have to be admitted at the meso level (size of

finite element).

The displacement finite element formulation is not fully consistent with the microme-
chanical lower (stress condition) bound calculation. This is due to the fact that in the upper
bound case the boundary displacements are random. Therefore, a need arises for modifi-
cation in the micromechanics computation to overcome this inconsistency. One possibility
is to perform a stress finite element analysis, which has as it is known, its own drawbacks.
Therefore, an alternative approach, which is adapted here, is to carry out an ensemble aver-

aging of Hooke’s law
§ (0)) = 52 (w) 9 (3.37)

to obtain

(*5) = <§§>o (3.38)

36

which results in a uniform strain field and linear displacements on the boundary of each i-
th finite element. Since the surface tractions are now linear on the element boundary and
there is no randomness present, this can now be used in a deterministic displacement for-

-r
mulation provided we replace 90 by (g) = (52;) (g) and e by (5).

Of course the upper bound on uac‘ is obtained by either averaging {ue;0) e O} or, more

simply, by employing (3.27) above with an ensemble averaged form of Hooke’s law as
(9) = (£980 (3.39)

Finally, it is interesting to make a note about the finite-size scaling, i.e. about the behav-
-r
ior of C; and (5;) , as function of 8, f(8). It follows from our computer calculations that

a fit can be achieved with a function

f(8) = 2265” (3.40)

where, a, b and p are constants. A fit by f (8) = a8” was tried, but could not cover the
entire tested range from 8 = 1 through 200. This subject, including finite-size scaling of

higher moments of Q; and .5; , is being presently investigated by M. Ostoja-Starzewski.

Let us note here that, effective average moduli calculated under periodic boundary con-

. 0 o -1 c o o
drtrons Ire between (9;) and (.52) . However, we did not investigate the 8-dependence
in this case, as we are not interested in random microstructures with finite scale periodicity

in our work on SFEM - such microstructures are not relevant here.

 

37

3.4 Methodologies of Solution

The foregoing observations in the previous section suggest a methodology of solution for
stochastic boundary value problems. When a statistical rather than a deterministic contin-
uum approximation is applicable, a solution of a specific problem is then conducted with a
certain choice of meso scale 8mso of RVE of the statistical continuum. However, since two
alternative definitions of boundary conditions are possible - deformation-controlled and
stress-controlled - then two different random anisotropic continua result. Thus, a given
boundary value problem may then be solved analytically or numerically, say, by stochastic
finite elements to find the upper and lower bounds on response according as Q; and g;
are used. However, for global response calculation, two distinct methods, which have sim-
ilar features, termed as “Exact calculation method” and “Covariance-based method” are
proposed here. In the first method, first and second moments of the random response
(mean and variance) are calculated in a Monte Carlo sense. Therefore, each random field
over the finite element is considered where its elastic tensors are calculated based on a suf-
ficient number of realizations. In the second method, the averages and correlation between
elements’ properties are employed to generate their elastic tensors. In such a treatment.
variation is admitted in the material properties only i.e. not in the response since the latter
is obtained by performing a deterministic finite element solution. The second approach
may be considered as approximate but it is more powerful since it dramatically decreases

the computational time.
3.4.1 Exact calculation method
It is termed “exact” since random microstructure of each finite element is discretized by

a very fine finite difference mesh for the purpose of calculating its material property in a

monte Carlo sense. Monte Carlo, which is a common statistical method in probabilistic

 

38

mechanics, is employed as a means of solving problems numerically through sampling
experiments. The problem may be posed in either a probabilistic or a deterministic form.
In a probabilistic setting, solution procedure via the exact method for response based on a

Monte Carlo method may be summarized in the following steps:

1- Qenemn'gg Qt Qne regiizgg’gn at the random micrgstructurg

We generate one realization 8(0)) of the random microstructure of the material at hand.
Actually we mean here a proper discrete model of the microstructure which will be
employed in the computation. For example, in case of a matrix-inclusion composite, round
shaped disks are thrown on a two-dimensional plane according to a restricted Poisson point

process and a specific volume fraction.

For a specific 8, we divide the domain into 2 x (n x n) finite elements for the purpose of

obtaining the global responses ue (15, 0)) and u" ()5, 0)) , of all nodes; 11 is the number of

square windows in X] or x2 direction.

3- M iergsrruemre egleulgtigns;

For each window (meso scale element) the "upper" moduli Q; (0)) and "lower" moduli
Q; (0)) are calculated based on the discrete random model introduced in step 1 and com-
puter simulations. Here we adopt a specific model; a random chessboard (with square lat-
tice) network where the window’s microstructure is descretized by a fine finite-difference
(as shown in Fig. 3.4) type of mesh representing the continuous phases of matrix and inclu-
sions, at'the microscale. Each node represents a grain cell from one phase with its stiffness
being modelled by a linear spring of either km or ki stiffness constant depending on which
phase it belongs to - matrix or inclusion, respectively. The interaction between neighbor-

ing nodes (phases) is modeled by these linear springs which are connected in series.

39

Assignment of phases to nodes is performed randomly according to a specified probability

distribution function - uniform probability density function is adopted here.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(b)

 

 

 

 

 

 

 

 

 

 

(a)

Fig. 3.4 Microstructure of a finite element discretized by a fine finite difference scheme (a),

one grain cell 1 with spring model (b).

Then, an effective Spring constant which models the interaction any two neighboring

nodes (cells) 1 and 2 is calculated as

—1
16"” = (i + i) (3.41)

Once the model is established, we apply essential (or natural) boundary conditions on the
window’s boundary. Consequently, solving a numerical boundary value problem for min-
imum energy stored in the network, under say essential condition, the elastic moduli
Q; (0)) can be extracted for 8(0)). Conjugate gradient method is used here as a minimiza-

tion algorithm for the energy.

4O

4- F' .
For each finite element, the moduli Q; (0)) or Q; (0)) are taken as the output of the above
micromechanical calculations. Feeding them into a finite element code and performing two

finite element analyses, two different responses are obtained.

Performing the above calculations (steps 1-4) for a sufficient number of realizations, Monte
Carlo sampling, the global response for each realization is calculated. Upon ensemble
averaging over the number of realizations, the upper and the lower bounds on the global

. . . . . 6
response are obtained depending on Wthh elastrc tensor rs used, 95 or 92.

In the above methodology of solution, an exact micromechanical solution over m
finite element was conducted. Then, a deterministic finite element solution for global
response was obtained for each realization. However, following the modification in the
micromechanical formulation when the natural condition is used, as presented previously,
the same steps of exact method still have to be implemented but in a different order. In
other words, effective moduli over one finite element are calculated via Monte Carlo
method. Since we deal with a statistically homogeneous medium, this result is then appli-
cable for all other finite elements. Henceforth, it may be fed into a deterministic finite ele-

ment analysis which yields a deterministic response.

3.4.2 Covariance-based method

In this method averages and auto-correlations between elements’ elastic moduli are used
to generate the random field of all finite elements’ elastic moduli. The method to achieve
such a generation is based on a simulation of random vectors technique (see e.g. Elishakoff

[30]), which may be summarized as follows:

41

Treating a finite element of size 5 as an RVE of the medium, the mean value of elastic
stiffness tensor of one typical element and its variance-covariance matrix are calculated via
the micromechanics-based method described in section 3.4.1. Furthermore, the method
relies on an observation that correlations between C tensors of any two finite elements are
known. In fact, as demonstrated by Ostoja-Starzewski [31], the correlations between any
two adjacent elements are almost zero, while those for any other pair are zero due to the
lack of spatial “memory” in our disk-matrix composite on scales 5 greater than the mini-
mum separation distance. Because _(_3 is a symmetric tensor, it is written here in a vector

form as
T

where {C} is the ensemble average (mean) of elastic stiffness tensor and V is the variance-
covariance matrix. Next, we decompose the variance-covariance matrix to upper and lower

triangular matrices [L], [U]T, respectively. This decomposition is expressed as

[L] 1U]T = (V1 (3.43)

where [1.]3x3 is given by
1/2 1/2 2 1/2
L“ = v“ L21 = v12/V11 L22 = (sz‘vrz/Vrl) (3-44)

A A A

. T
Then, to generate the i-th element’s elastic moduli { C} = [C1, C22 C12], we use the

relationship

{613... = 1L1,,, (2);... + (6)3... (3.45)

42

where, {Z} is a vector of random numbers taken from an independent standard normal dis-
tribution N(0,1). Thus, we generate a realization of a random vector for each finite element
in a sequential fashion. The corresponding elastic moduli tensor Q (0)) is obtained by
inserting the generated random numbers in the formula given by (3.45). Once all elements’
elastic moduli are assigned, a deterministic finite element solution is performed to obtain

the global response.
3.5 Numerical Results and Discussion

The micromechanics-based stochastic finite element is used here to solve two problems;
an elastostatic membrane and a steady state thermal conduction. In both problems, the
microstructure is taken as a matrix—inclusion composite with its microstructure modeled
under the following specifications:

1) both phases are locally isotropic homogenous continua;

ii) inclusions are round non-overlapping disks of the diameter d = 5A1, where Al is a unit
length,

iii) the distribution of the inclusions in the matrix corresponds to planar Poisson process
with a specific minimal distance between any two adjacent points. The reason of choosing

such a restricted process is to avoid the stress concentration problem.
I) Membrane problem:

The out-of-plane deformation u (x1, x2) will be obtained by applying the stochastic finite
element method corresponding first to a random continuum approximation 82 , and then to
82, at a given scale 8. Thus, Ci] in (3.11) are components, and realizations, of the random
tensor fields 9; or 92- The membrane problem is being expressed in the form of a Poisson

equation (3.11), under Dirichlet boundary conditions u = 0 on 88 of a square-shaped

43

body domain of the size 400Alx400Al, subjected to a uniform force distribution f = 10°2/
(Al)2 (see Fig. 3.5). Inclusions are generated with a minimal distance between any two
neighboring inclusions of140% of inclusion’s diameter. Inclusions occupy 10% of the

whole volume.

 

 

 

 

 

 

 

 

 

 

 

 

Fig. 3.5 A membrane of a matrix-inclusion composite with a finite element mesh of reso-

lution 8 = 20.

To take into account the randomness in the inclusions’ distribution, ten different 8(0))’s
are generated and the solution for the global response is obtained. As a measure of the glo-
bal response of 8(0)), we choose the volume V (0)) under the membrane, so that we obtain

a set of 10 ‘upper’ estimates V (e) (0)) and a set of 10 ‘lower’ estimates V (n) (0)) . Figure

44

3.6 shows the ensemble averages (V M) and (V (n) ) for three different values of 8 plotted
as functions of increasing stiffness Ci of the inclusion with Cm kept constant. Here Cm
(stiffness of matrix phase) is taken as unity for simplicity, so that Ci becomes a so-called
contrast, and all the plots are nomralized by V“) = V (n) at Ci = Cm = 1, which is the

purely deterministic case.

 

 

 

 

1.0000
0.9500- ' _
- a)
E \ . -
2 0.9000- \ ._
° \
> \
\\\
a) 0.8500- \Qé‘:~~-___ 6‘”) = 4 ..
.5! ~:~:,=: _:'§""-- 5(0) = 8
B _ “*—?.~ 50") = 20
E 0.8000- 5") = 20 _
0 5(9) = a
2
0.7500- ' 5“” = 4 _
0.7000 , , , , .
1.0 6.0 1 1.0 16.0 21.0 26.0 31.0
c3 '

Fig. 3.6 Ensemble average “upper” and “lower” responses, normalized over the deterrnin-

istic case Ci = Cm = 1. 8“) and 8(“) correspond to C; and Q2 , respectively.

As may be expected, the membrane deformation under essential boundary condition is
smaller than that under the natural boundary condition, since the moduli Q; are always

. fl . . . . .
stiffer than the 96' On the other hand, an increase 1n the stiffness of the rnclusrons leads to

45

stiffening of the membrane, and hence, to a reduction of the deformation and consequently,
of V“) and V“). Therefore, for a fixed 8, both curves are decreasing monotonically and
diverging away from Ci = Clm = 1, with increasing C. In addition, for a fixed contrast, we
see that the two responses (bounds) get closer with increasing 8, and have a tendency to do
converge to a unique value as 8 -—> oo , which corresponds to the deterministic case. How-
ever, we should note that this limit can only be attained in an approximate sense, providing
the fluctuations in stiffness disappear for such large finite elements. Thus, it is expected
that, for a fixed contrast, the scatter in the responses (as measured here by the standard devi-
ation) decreases when 8 goes up, (see Fig. 3.7). On the other hand, for any fixed 8, it

increases with the contrast (i.e. microstruCtural randomness) increasing.

 

 

 

 

7.00
- ~ . 5‘91 = 4
'3 6.00- —
o
,_ _ (e) =
x 5.00- - 5 8 -
8
:3;- 4.00- -
02
> , 5(0) = 4
8 3.00- ’,./” ..
'0 .3”
'8 2,00— ’/:.____________..—-- 5(9) = 20 ..
g // .___ ________ 8W — 20
'2; 1.00- / ---"'”" -
0’00 l I I I I I
1.0 6.0 11.0 16.0 21.0 26.0 31.0
Ci

Fig. 3.7 Graph of standard deviation of the “upper” and “lower” responses.

46

II) Heat conduction problem:

In this section a problem of a heterogenous medium, undergoing a steady state heat con-

duction, is considered. The local constitutive relationship at a point x in the domain is

given by

q; = —Kij ($3 (1)) TJ- (3-46)

where,

qi: heat flux,

Kij: second rank thermal conductivity tensor,

TJ: temperature gradient,

0): one realization taken from the sample space Q.

According to our micromechanical formulation, we have two constitutive relationships in

the form of

g = —I_{§YT g = —I_("5YT (3.47)

over each finite element. Thus, all our concepts and observations presented in section 3.3
still hold here. Therefore, two stochastic finite element problems have to be solved to find

the upper and lower bounds on the actual response.

Our solution of this problem takes into account several features:

i) Both exact and covariance-based methods of solution are adapted in stochastic finite

element analysis.

ii) For lower bound problem, the compatibility modification in the micromechanical com-
putation is incorporated. Accordingly, averaging are conducted for the compliance tensor

and then a deterministic finite element solution is sought.

47

iii) Unlike the membrane problem, window, i.e. meso element cell, will be equivalent to

one triangular finite element.

With these features taken into account, we focus on the response as a function of the
element size 5. The CPU time for both finite element and micromechanics computations
will be measured and compared. To achieve these objectives, we consider a unit square
domain of the medium in which the inclusions are distributed randomly in the matrix with
a contrast of 10 according to a minimal distance between disks centers being 140% of inclu-
sion’s diameter and volume fraction of 0.20. The domain is subjected from above to a uni-
form source of heat with unit magnitude. The boundary conditions are specified as two

insulated neighboring edges, with the other two are flux free (see Fig. 3.8).

48

 

dT/dx = O

 

 

 

 

Fig. 3.8 The domain, boundary conditions, and finite element discretization of the

domain with 6 = 10; a plot at 402x402 pixel resolution.

The temperature at the corner point (the origin of the coordinate system) is used as a
measured reference in our present problem. Firsr, considering the case of a homogenous
microstructure (without inclusions), a deterministic finite element solution is found to be T
= 0.2913, which is very accurate compared to the series solution (T = 0.2965) and Reddy’s
solution (T = 0.3013) [32]. In the random medium problem, we study the effeCI of finite
element size (5: the ratio of meso element size to micro element size) on the results. Two

sets of bounds on temperature based on exact and covariance-based calculations are

49

obtained and plotted as function of 5 (see Fig. 3.9). Six values of 5 are considered 2, 4, 8,

10, 16, 20 and fifty realizations are used in the ensemble averaging.

0.4

 

 

0.3 -

 

Temp
‘1
'9
t
t
t
L

 

 

------------------- I—--------------
02— .>17 """"""" 4: —---'k. + J

 

0.1 — + Exact, essential
—0— Exact, natural
-—-I--- Covariance. essential

“0"" Covariance, natural

 

 

 

 

 

 

0.0 I I I I I I l I I l I I I I I I I I I
0 5 10 15 20

Fig. 3.9 Bounds on response (reference temperature of corner point) as a function of 5.

Examining the results of Figure 3.9, we can observe that the difference between borh
bounds decreases as 5 increases. It is with increasing 5 that both bounds converge to a con-
stant value which corresponds to the deterministic solution. The window’s (finite element)
size plays an important role in both micromechanical and finite element calculations. For
example, in the micromechanical computations, it is of interest to have a large size of win-
dow in order to decrease the fluctuations in the material property. On the other hand, finer
finite element mesh leads to a better accuracy in response results. Comparing the two sets

Of bounds, under the exact method and the under random assignment method, the former

50

set apperas to be softer than the latter one. The difference is not that significant which
makes the latter one more reasonable to use. As expected, both methods provide conver-
gent bounds on response with 5 increasing. However, they are in some disagreement; this
is due to the fact that the covariance-based method is an approximate one since:

i) the elastic moduli are generated according to a Gaussian distribution based on the first
and second moment information; actual statistics of C; and S; are not really Gaussian
[31].

ii) the elastic moduli are generated assuming zero correlations between contiguous finite
elements; this is not exactly true as some spatial correlation exits in a non-overlapping ran-

dom disk field.

Since randomness in the response is involved in the upper bound only, standard devia-
tion has a meaning here. Thus, it is plotted as a function of 5 as shown in Fig. 3.10 which
shows that scatter decreases as 5 increases. It should be noted that neither the full conver-
gence at 8 —> cc of upper and lower bounds in Fi g. 3.9, nor the complete decrease of scatter
at 8 -—> oo in Fig. 3.10 can be attained since the finite element has to remain finite relative

to a heterogeneity (i.e. microscale) in a finite element problem.

51

 

St. Dev. x 10 e 3

 

 

 

 

Fig. 3.10 Standard deviation of the upper bound (for reference temperature) as a function

of 5.

CPU time is a very important criterion in determining the practicality of any computa-
tional technique. Here, we measure and compare the CPU time for both micromechanics
(MM) and finite element (FE) calculations based on one realization as shown in Figure
3.11. From that graph, we can see that the FE time decreases dramatically with the increase
of 5, while the MM time is increasing very rapidly. This is due to a decrease in the number
of finite elements, but at the same time, an increase in the element size (the number of res-
olutions increases). This indicates that an optimal 5 may be deducted which is found here
to be approximately 3. On the other hand, the advantage of using the covariance-based
method rather than the exact one is that it reduces the nricromechanics calculations time by

two order of magnitude. Nevertheless, the CPU time of either one of our methods (exact

CPU time, minute

52

or covariance-based) is much less (1) than that of the problem of exact solution using very
fine resolutions of a finite difference scheme over the entire domain. In such a case, CPU
time was found to be 443 minutes, whereby the Conjugate Gradient method was used as

the minimization procedure. This comparison shows that our methodologies are efficient

and practical in numerical solution of Stochastic boundary value problems

 

 

—
"a
c
-’

 

 

 

 

 

 

 

Fig. 3.11 CPU time of FE and MM calculations plotted as a function of 5.

CHAPTER 4

EFFECTIVE ELASTIC MODULI OF DELAUNAY NETWORKS

4.1 Introduction

Variation in material properties at micro level as well as random geometry of micro-
structure have a great impact on the effective (macroscopic) properties of heterogenous
media. Particularly, we are concerned here with two-phase granular media. Deterrnina-
tion of constitutive laws for such media is one of the primary goals of micromechanics.
Besides, Micromechanics gives a better understanding of the relationship between the glo-
bal behavior of such materials and micro behavior. In this context, one can see that differ-
ent means of study were used: analytical, numerical computations, and experiments with
real or model material. We refer the reader to a special issue of Mechanics of Materials

Journal for a collection of recent advances in the granular materials area [2].

Naturally, behavior of granular materials is influenced by a multiplicity of different
parameters. Anisotropy, contact points between particles, type of packing (i.e regular or
random) and displacement and rotation of individual particles at micro level are some
examples of these parameters. Incorporating one or more of these parameters, many
researchers have attempted to derive constitutive models at macroscopic level. Perhaps
the statistical approach is more frequently used than other methods. As an example, a sta-
tistical method was employed in defining the stress tensor for particulate media [33] and
hence in finding the effective stiffnesses for assemblies of granules with linear elastic
interactions [34]. In this regard, many studies modeled the grains at micro level as circu-
lar (spherical) disks with one contact point between any two neighboring grains by which

53

54

force can be transmitted from one to anorher. Among those, Satake’s [35] work who pre-
sented a formulation of new mechanical quantities where void and contact points were
introduced. While Oda [36] was concerned with the importance of fabric on mechanical
response of soils, Satake [37] concentrated on the analogy between governing relations of
discrete and continuous media, as well as on relating the anisotropy of microstructure with

stress and strain.

The contact between grains at micro level received a great attention of many workers in
this field. Based on statistical approaches, a number of studies have been attempted along
this line of approach, see for example work [38 - 40]. These works were limited to condi-
tions of small deformation without sliding at inter-particle contacts. Chang [41] dealt with
translation and rotation of discrete particles which yield the micro scale discontinuity due
to sliding and separation at contact points. At micro-element level, stress and strain were
introduced so that a stress-strain relationship based on contact force transmission was
derived. On the other hand, Oda et. al [33] introduced several measures of granular fabric
for a random assembly of spherical granules. Therefore, parameters which characterize the
granular fabric were utilized in a description of the overall macroscopic stresses and
deformations. In a general frame of continuum mechanics, Sidoroff et. al [42] considered
a description of the contact force distribution in granular media, taking into account the
induced anisotropy. Based on contact points between particles, Mehrabadi et al [34]
examined the macroscopic stresses of a granular mass and, hence, defined these stresses

in terms of its microsnuctural characteristics like the fabric.

An alternative approach to analytical methods of solution in granular media is that of
numerical simulations. If used in a proper way and checked by real physical experiments,
computer simulations can provide a wide range of valuable information not easily avail-

able from real experiments. This type of approach has great advantages in the micro level

55

analysis of granular materials: motion of each individual particle and mechanical state of
each contact point can be followed; boundary conditions, material properties, etc. can be
prescribed exactly. Recently, Rothenburg and Bathurst [43] used numerical simulations to
study a medium of initially isotropic plane assemblies of elliptical-shaped particles. They
focused in their work on the strength characteristics of real sand and micromechanical fea-
tures of behavior. In another study [39], they presented a micromechanical analysis of
plane granular assemblies interacting according to linear contact force interparticle rela-
tionship. Explicit expression for elastic constants in terms of rrricrostructure’s parameters
were derived and verified by computer simulations for dense isotropic assemblies. How-
ever, one drawback of the computer simulation method is that its power is limited by the
capacity of computers which deal with a limited number of particles. This fact motivated
some experimental line of approach which can be found in literature. Gourves [44], as an
example, used various experimental techniques to measure displacement of each grain,

stresses and other parameters.

Recent decades have seen an increasing use of the concept of graph theory in the anal-
ysis of solids with heterogenous microstructures. Different kinds of networks whose
geometry ranges from regular to random were adapted in modeling these microstructures.
Motivated by the lack of models which connect random field models to underlying micro-
structures, Ostoja-Starzewski [45] showed how a graph representation may be used to ana-
lyze the constitutive laws of granular-type materials with the inclusion of micro scale
randomness. For an arbitrary polygonal convex-shaped cells in a graph, the author
derived constitutive coefficients for linear elastic interactions between grains in which
randomness was incorporated per edge. Davis and Deresiewiecz [46] used a graphical
approach to study the compressibility and force transmission in a granular material mod-
eled as a two-dimensional random packing of like spheres in elastic contacts. They repre-

sented the packing by a stochastic planar graph, with the nodes of graph taken as the

56

centers of spheres and the branches as contacts between adjacent spheres. Then, the sto-
chastic graph was replaced by a lattice whose each branch has a random stiffness modulus

assigned to it.

It is well known that the Voronoi tessellation, is an important model in different fields
of science. It was suggested and improved by several workers like Meijering, Cahn and
Prager based on the idea of subdividing a space into random regions where each region
corresponds to one Poisson point. Many workers in the random media field adapted the
Voronoi model in their studies. Meijering [47] introduced the Voronoi tessellation as a
model of crystals aggregate with randomly distributed nucleation sites. Uniform growth
of these sites, represented by randomly distributed points in space, resulted in space filling
Voronoi polyhedra. Hatfield [48] generated random bond networks from a two-dimen-
sional Voronoi tessellation as an attempt to model the structure and conductivity of micro-
emulsions. For more discussion about Voronoi (and other types) tessellation the reader is
encouraged to refer to Winterfield’s dissertation [49] entitled “Percolation and Conduction
in Disordered Media” in which he presented it as a graphical model for chaotic compos-
ites. In this dissertation, a further discussion about generation of a Voronoi graph will be

presented in light of what is given in [49].

Dual to the Voronoi tessellation graph is a so called Delaunay triangulation, i.e. trian-
gles or there analogous higher dimensions. A Delaunay triangulation network is con-
structed by joining a pair of contiguous Poisson points by a line which is bisected by the
in-common edge of the Voronoi cells containing the two points. Linking all Poisson
points according to this mentioned criteria makes up random triangles. Both networks
(Voronoi and Delaunay) are generated according to a stationary Poisson process in R2
plane. More details about construction of Voronoi and Delaunay networks will be pre-

sented in next section. Delaunay graphs (networks) have been used in modeling the inter-

57

actions at micro level of solids. Consequently, effective properties, like elastic, thermal, at
macroscopic level may investigated. For linear elastic Delaunay networks, Ostoja-Starze-
wski and Wang [50,51] calculated the stiffness matrix in terms of the first and second
moment based on computer simulation. For a homogeneous medium, bounds on macro-

scopic moduli of Delaunay networks were found as a function of some of its parameters.

In the present study, we conduct a study for a special class of two-phase granular media
with a topology of Delaunay networks. A special feature of this study is that it establishes
a linkage between different concepts: discrete graph models, random field models (sto-
chastic concepts) and numerical (computer) simulations which are used as a methodology
of solution. Our goal here is to identify a proper effective medium theory (self-consistent

model) based on these computer simulation results.

4.2 Problem Statement

Consider a two-phase random granular medium (see Fig. 4.1) in 2-D plane, whose
microstructure is modelled by a Voronoi planar tessellation graph. The central interac-
tions between grains at micro level are modeled by Delaunay graph (network) in which
each vertex represents a unit cell (grain) of the medium, while the edges represent the lin-
ear elastic central interactions which are modelled as linear springs. The shear interactions
between neighboring cells are modeled by linear angular springs connecting two consecu-
tive bonds in the Delaunay graph. A discussion about various models of microstructure
can be found in the next section. Based on the two spring models, we seek effective elas-
tic moduli of periodic Delaunay networks at macroscopic level. Although properties at
micro level are locally anisotropic, statistics of Delaunay networks shows that macroscop-
ically they are homogeneous and isotropic. Therefore, two effective constants (for exam-

ple, K, [1.) are a sufficient information to provide.

58

The problem of determination of effective moduli of Delaunay networks falls into a
general category of linear transport problems. The randomness in those networks may be
per edge or per vertex. The first case is easier to handle and has been solved in a number
of conductivity and elasticity problems [50, 51, 52]. The second case is especially hard in
the case of a topologically disordered microstructure lacking periodic geometry. Such is

clearly the situation with Delaunay networks.

 

Fig. 4.1 Random assemblage of particles from granular material with spring interaction

model.

In this study, we are interested in using the second approach (randomness per vertex) to
determine the effective moduli of a two-phase granular medium with Delaunay network
type of microstructure, where the two phases are mixed at random. Vertices in the net-
works are assigned in a random way to be either of type 1 (soft) or type 2 (stiff) with the
volume (i.e. area) fraction C] and C2, respectively [53, 54]. More specifically, type I (or
type 2) vertex implies that the halves of all edges incident onto it are soft (or stiff). In par-

ticular, the stiffness of any half length 1/2 ofa given edge is given by 2ki/l, where i = l. 2.

59

One edge of Delauanay graph represents two springs (which belong to two neighboring

grains) connected in series. Thus, its effective stiffness constants are calculated according

to

eff 1 1 ’1
= — — 4.1
k (k1 + k2) ( )

The stiffness ratio r = k2/k1 is called a contrast. When r = 1, the model reduces to the
homogeneous one-phase case. Therefore, the randomness is present not only in the geom-
etry of Voronoi cells, but also in the distribution of the two phases. For the sake of com-
parison, the effective moduli of the medium will be determined using regular (triangular)
networks of approximately the same size number of vertices. Macroscopically, elastic
moduli of s triangular lattice modeling a two-phase medium are isotropic. For a homoge-
neous case, an exact analytical formula can be derived considering only one unit cell.
Based on these two models of microstructure (Delaunay and regular), we seek the specifi-
cation of Hooke’s law as

e e
or} = 71 8050+ 20 51,- (4.2)

where, he and [.1e are the effective Lame constants.

Numerical results of effective moduli are obtained via computer simulations for
Delaunay networks using two models; (1) central and (2) both central and angular elastic

springs. For regular triangular networks only the central spring model will be utilized.

On the other hand, we attempt to identify a proper effective medium theory for analyti-
cal calculations of effective elastic moduli of such two-phase Delaunay networks. Our

computer simulation results are used to determine which of several different self-consis-

60

tent models can provide the best possible approximation to effective Hooke’s law.

Following reference [55], we employ two self-consistent approximations (S CA):

i) a Symmetric Self-Consistent Approximation (SCA-S), which treats the medium as a

mixture of N phases without distinguishing any one as a matrix.

ii) an Asymmetric Self-Consistent Approximation (SCA-A), which treats the medium as a

system of N-l inclusions embedded in a matrix. The SCA-S is specified by the following

equations

N N
XCiPz(K-Ki) = 0 2400—11,.) = 0 (4.3)
‘=1 i=1

while the SCA-A is given by

(K1 — 19-)

1 1 N 1 l N (Ill-14;) '
s = E, 1* 220107" 11 = ii}- 1’ 2240—7,— (4'4)
1: I:

where the matrix phase (i.e. i = 1) serves as a reference. In the above coefficients Pi and

 

 

 

Qiaregivenby
P 1 Z ZiYiil-sY—l 2 PHZ" 2 x *1 4
i-[ + i—Xi+Y‘. 1+3 ] Q‘ — ‘X,+Yl.+( +Ai_ i) (’5)
where (s) is the ellipse aspect ratio and
A.s(1+v) A.(1—v) C.(1+v)
Xi = + ._‘_____2._.. Y‘, : ..._‘___ ZI : —‘———— (4.6)
(1+5) 2 2
p. K- _.
4,=—‘ C,=—‘—1 1):" “ (4.7)
p K ic+u

61

To apply the self-consistent models to our media, we perform a mapping between the
Delaunay network and an inclusion-matrix system where the former is viewed as if it is a
field of inclusions, rather than vertices connected by elastic edges, see details in next sec-
tion. The main advantage of establishing such effective medium theories is that it gives

approximate results without conducting computer-intensive calculations.

4.3 Models of Random Microstructures

Macrostructures of a planar two-phase granular materials are represented by two differ-
ent discrete models; Voronoi graph and regular hexagonal lattice. While Voronoi model
captures the randomness in microstructure’s geometry, hexagonal lattice is a proper model
for special types of microstructures, like that of ceramics. In both models, each cell repre-
sents one particle at micro level and interactions between grains are modeled by Delaunay
and regular triangular networks, corresponding to Voronoi and hexagonal graphs, respec-
tively. Since the regular lattice is a well known model, it will be used just for the sake of

comparison with our Delaunay results.

In the Delaunay network, two spring models will be used; central, and combined cen-
tral and angular (noncentral) spring model. Thus, an approximate continuum model for
heterogeneous medium is defined from these discrete models. In what follows, we give a

brief description of each model.

4.3.1 Regular Triangular Networks (Lattices)

Regular triangular network is a well established model which was used by many work-

ers, to represent macrostructures of composite and crystalline materials, e.g. see [56-59].

In this model, material in a continuous sense has cells of hexagonal shape which have

62

same size and orientation and cover whole area of the medium. When all cells of a lattice
have the same property, the medium is termed as “homogenous”, while if it changes from
one cell to another it is called “heterogenous”. Randomness in phases’ distribution may
be incorporated by random assignments of springs’ stffinesses. Representing each cell
(hexagon) by one vertex (node) at its center, and interaction between cells as bonds con-
necting all vertices, a regular triangular network, named as lattice model, is obtained (see
Fig. 4.2). Clearly, it is equivalent to a finite difference scheme in the sense that it is a dis-

cretized representation of the original continuum (i.e. nodes and bonds).

 

 

 

 

 

 

 

 

 

Fig. 4.2 Triangular (regular) network of size 7x7 with a hexagonal unit cell.

Consider a single cell, as shown in Fig. 4.3, of the discrete model; we describe the orig-
inal continuum in terms of six normal linear springs (Note: we assume here central inter-
actions only). These springs control the stretching of the span OL, i.e. the deformation of
the cell. For a two-phase medium, each hexagonal-shaped cell springs may be either of 1
type or 2. To answer the question of what material this kind of model can describe, we
need the relation between the stiffness of springs and elastic moduli of the continuum it
models. MacrOSCOpically, triangular lattice model yields isotropic properties of a random

two phase medium. Further discussion can be found in sub-section 4.3.3.

63

 

Fig. 4.3 Unit cell with six springs modeling the stiffness of the original continuum.
4.3.2 Dual Voronoi -Delaunay Networks

Tessellation of a plane means subdividing it into smaller regions with a certain shape
like square, honeycomb, triangle etc. that cover the whole plane without any overlapping.
Two types of tessellation may be distinguished in literature; regular and random. One type
of random tessellation in 2-D plane is the Voronoi tessellation in which a plane is subdi-
vided into cells of convex polygons. The two-dimensional tessellation has appeared in
the literature under the names of Dirichlet [60] and Thiessen [61] as well as Voronoi. For-
mally, Voronoi planar tessellation graph G (E,V) consists of two sets V - the set of vertices
of cells and E - the set of its boundary edges (see Fig. 4.6). Generalized to N dimensions,
the definition of the Voronoi tessellation as presented by Winterfeld [38] may be summa-

rized as:
Let P1 (x1), P2 (x2), ...Pk (xk) be an infinite set of randomly distributed Poisson
points, according to a specific distribution, in N-dimensional space. The region inside the

polytope of Pi is the set of points x closer to Pi than PJ-

|.§—)§j| < lx—xil [:tj (4.8)

64

Points on an edge of a Voronoi polygon are equidistant from two Poisson points Pi and Pk.
lat—26,] = [at-3,] <|2c—r,| l¢i.k (4.9)

An edge lies on the perpendicular bisector of the segment connecting Pi to Pk. Vertices of

a Voronoi polygon are equidistant from three Poisson points.

[at-2,] = lat—2,1 = |i-J£,|<|Js-r,,,| m¢i,k.l (4.10)

A vertex is the intersection of the perpendicular bisectors of the segments connecting P, to

Pk and Pi to P1.

The process to generate a two dimensional tessellation network, shown in Fig. 4.4 may

described as follows:

Starting with a random distribution of Poisson points in an area, all points are simulta-
neously allowed to expand into circular disks at the same rate. The constraint - no disk is
allowed to impinge on the area of another - is imposed during the formation process.
Thus, it must deform upon further expansion with all disks evolving into polygons of
which each edge is equidistant from two Poisson points. Another characteristics of the
Voronoi graph is that any circle centered at any point inside the Voronoi polygon with a
radius equals to the distance of this point to the polygon’s Poisson point contains no other
points inside [47]. In addition, the circle centered on an edge of the polygon passes
through two Poisson points while the circle centered on a vertex passes through three

Poisson points.

 

 

*9.” qn-um—ns-r" ,. ,_.

65

 

 

 

 

 

 

 

 

(a) (b)

 

 

 

 

 

 

 

 

 

(C) (d)

Fig. 4.4 Formation process of a Voronoi graph. Poisson points (a) expand into disks at
the same rate (b). No disk may impinge on the area of another (0), hence, they evolve into

area-filling polygons (d).

 

 

 

 

 

 

 

 

 

 

Poisson points in a un iiiiiiiii

 

 

 

68

We employ Voronoi networks as two-dimensional models of multi-phase granular
materials. The Voronoi tessellation determines a structure of simply connected grains that
fill an area. Voronoi polygons are the building blocks of a random material which impart
to it a random geometry and topology. The random topology results from the variation in:
coordination number of one site to the other in the Voronoi graph, the random location of
Poisson points, and the random sequence of labeling polygons as either one of the two
phase. The random geometry results from this labeling and the geometrical statistics of
the tessellation process. From the Voronoi tessellation, a useful model of two-phase gran-
ular media is obtained by selecting at random a given number of polygons and labeling
them as a certain phase and the rest as the other phase. The network representation of the
Voronoi tessellation is dual to a so called Delaunay network G(V, E) which consists of two
sets - V set of sites (vertices) representing polygons of the Voronoi graph and E - set of
bonds (edges) connecting the contiguous polygons. This network representation assists in
determining the effective properties of the Voronoi granular media. Once the medium is
created, its structure is fixed and its elastic moduli can be determined via numerical com-
putation after a proper spring model is introduced into the network which models interac-

tions between particles.

To perform a computer simulation, the Delaunay network is translated digitally in a
form of three matrices stored in different files. First one, called a connectivity matrix,
contains (N) rows where N is the number of Poisson points of the network. In each row
the first entry is the Poisson point (vertex) label and in the second the coordination number
(# of nearest neighbors vertices) while the rest of entries represent the nearest neighbors to
this vertex. Statistics of Delaunay networks shows that the average coordination number
per vertex is six, recall Euler’s theorem [62]. The numbering scheme used for the matrix
is the order in which the polygons are constructed, i.e., row #1 represents the first polygon

constructed, #2 the second, etc. In the other two files the X-coordinates and Y-coordinates

69

of all Poisson points are stored, respectively. Such generation of Voronoi-Delaunay net-
works is implemented by a computer program written in FORTRAN and run on a SPARC

10 - SUN station.

4.3.3 - Central Spring Model

In a generally disordered network, any edge, which represents the central interaction
between two adjacent cells is modeled as two linear springs connected together in series,
each of constant Oti, where i = 1 or 2 in a two-phase medium. Analogously, this model is
applicable to a special case of the regular triangular lattice. Considering a cell, in a shape
of polygon, as shown in Fig. 4.7 - b, its mechanics is simulated by m central (normal)
springs which are connected at its vertex. We should notice here that all springs incident
into one cell posses the same stiffness since we assign stiffness to cells (vertices) not to
edges (springs). Nodes of cells can deform freely in a 2-D plane i.e. have two-degrees of
freedom (ul, u2). The m-th edge (OL) makes an angle 9 with respect to the coordinate
axis X1 with its rotation being negligible. Also, since we are dealing with infinitesimal

linear plane elasticity problem, higher order terms in energy expressions are neglected.

The relation between elastic moduli of the continuum and spring stuffiness of the dis-
crete model is established via equivalence setup in which energies of both models are
equal. For the continuum model, the elastic energy of a single cell is expressed as

1 . .
E = ECijkleijekI/l 1,), k = 1, 2 (4.11)
where,
Cijkl : forth order tensor of elastic moduli (stiffness tensor),

EU: a second order linear strain tensor,

70

A: area of the unit cell.

 

 

X
2‘
L/2
u2
I, L
x “1
' 0

O >

X1

(a)
(b)

Fig. 4.7 Central spring model for an in-plane elasticity problem (a), Voronoi unit cell with

m springs incident into it (b).

To obtain an energy expression for the discrete model, Fig. 4.7 - a is useful where OL
represents the m-th spring of stiffness 0t. Without loss of generality, let the length of the
spring be included in its stiffness constant. i.e spring constant is inversely proportional to
its length (k ~ 1/1) . Assuming that O is a fixed point and the in-plane infinitesimal dis-
placements of point L are 111 and uz in X1 and X2 directions, respectively, then, the total

energy, stored in all of springs, is given by

n n

_ m_8[j8k[ mmmm
E— 2E -37 6111.11.51,“ (4.12)

m=l m=1

where, 11 = c050, I2 = sin0. Comparing expression (4.11) and (4.12), the elastic moduli

of the discrete spring model are found as

 

 

Cir/d = in: 2 am’inlitlin’i" (4.13)

1 .5

For the triangular regular lattice model A = 2J3 and I1 = 3,12 = —2—. It is obvious
that expression (4.13) for the elastic moduli satisfies the essential requirements of symme-
try in stress and strain tensors. Result presented above was found by many workers, e.g.
see Holnicki and Rogula [63] for regular lattice which was extended by Ostoja-Starzewski

to a general case of disordered geometry [45].

Suppose a triangular lattice of central spring model is used to model a homogeneous
isotropic medium e.g. [64, 65]. Thus, in any cell, all the six normal springs have same

stiffnesses; in such a case, (4.13) can be rewritten as

C3“: 53-”; 1”“1, (4.14)

Substituting 11 = 0.5 and 12 = #2 , into above, elastic moduli can be evaluated as

3
C1111 = C11 = §./301

J3
C1122 = C12 = Ya (4.15)
./3
C1212 = C66 = ‘g’a
which satisfy all the following conditions of isotropy
C1111 = C2222
C1112 = C2212 = O (4.16)

l
C1212 = 5(C1111" C1122)

72

In a such a 2-D case, we express bulk (and shear) moduli in terms of stiffness matrix com-

ponents as

1c+u = C11 [1 = C66 (4.17)

so that a pair of independent elastic constants are obtained as

J3 ”/3 (4.18)

K=—a li=—8—a

Substituting the elastic moduli of expression (4.18) into the following relationship,

1:
(3

K' — ‘2 (4.19)

V=K+1J._C;;

 

we obtain a Poisson ratio (V) equal to 1/3. This result which also may be obtained form an
analytical solution for one unit cell indicates that this model can be used only for a special
case of two-dimensional medium of V being 1/3. To present a medium of a general isot-
ropy, two independent elastic constants i.e. E and V, another parameter like angular

springs should be added to the model as will be discussed next.

4.3.4 - General Spring Model - Combined Central and Angular Spring Models

As shown by Fig. 4.8, an angular linear spring linking two adjacent normal springs in
order to constrain their free rotation, is introduced into the network. For each unit cell, we
have (n) angular springs, where n is total number of normal springs. The stiffness of the
m-th spring, denoted by 5‘“, takes into account the magnitude of angle associated to it. We
notice here that all angular springs incident into a cell have the same stiffness and there is

no interaction between springs of different cells.

73

112

 

 

 

Fig, 4.8 Combined central and angular spring model.

Writing an expression for the elastic energy due to angular change (AG) (i.e. angular

spring deformation) between two bonds as shown in Fig. 4.8 in the form of

E’" = épmlAelz (4.20)

”1+1 m
0

where, AO = A —A0

After carrying out some manipulations and making use of the stress-and strain symmetry

property in stress and strain tensors, the total energy can be expressed in terms of the dis-

74

crete model parameters. Again, from energy equivalence between those of discrete and

continuum models we obtain the elastic moduli as

(1,8111%Bm“)s..z;'tz;"t-(t'"+1““)1".it"

_Bm5ikl;rllm+ll;n+ll;n+ Bmén+1171r+117
K -Bm0ikl:1?1:+1l7+l+Bmén+14nlzzl§n+l )

C2.“ = (4.21)

 

 

where superscript (m-l) equals 11 when m = 1. Here n is the total number of springs in the

cell. Deriviation of above formula can be found in an unpublished work of Ref. [66].

Since, we are dealing with an elastic behavior, superposition of elastic moduli given by
(4.11) and (4.21) yields the elastic moduli of a cell of the continuum expressed in terms of

the stiffness constants 01 and B of normal and angular springs, respectively.

We consider a homogeneous isotropic medium modeled by a triangular lattice. All nor-
mal springs are assigned same stiffness (01), and angular springs same stiffness (B). In a
continuum sense, elastic bulk and shear moduli of a cell are expressed in terms of the dis-

crete model’s parameters (01 and B) as follows

. = 5383)

1 (30+28) (4.22)

“-573

Eqn. (4.22) implies that the presence of the angular springs in the discrete model has no

13

effect on bulk modulus (i.e. K = 0). Substitution of these values into the formula of

two-dimensional Poisson’s ratio given by (4.19) yields

0—30

 

75

Thus, the above equation suggests the following range of Poisson’s ratio that can be cov-

ered by the present model

(4.24)

4.4 Methodology of Solution
4.4.1 Spider-Inclusion Analogy

The Voronoi tessellation can be viewed as a planar system of continua, e. g. a polycrys-
talline solid. If we observe that the self-consistent method is used to predict effective
(macroscopic) moduli of polycrystalline media, based on approximations of a typical crys-
tal by an elliptical inclusion, such moduli of Delaunay networks can be determined pro-
viding we set up a suitable equivalence between both systems. On the other hand, the
condition that self-consistent methods are only applicable to macroscopically isotropic
medium is also consistent with our Delaunay network system since it is statistically an iso-
tropic one. Due to the fact that a typical Voronoi cell is somewhat elongated, it seems rea-
sonable to approximate it by an ellipse of aspect ratio (s). Actually, the computational
results which will be obtained by computer simulations in a range of volume fraction
between 0.0 and 1.0 will be used to determine the closest value of s, the aspect ratio of
ellipse, such that the self-consistent model fits our exact results. To that end, we define a
“spider” to be a Delaunay vertex and those halves of edges incident onto it which lie in
this vertex’ Voronoi cell. Now, with the help of Fig. 4.9- b below we call this equivalence
a “spider-inclusion analogy”. The precise nature of this analogy depends on the actual

inclusion model, that is on the choice of shape, and the type of inclusion-matrix interface

76

(i.e. perfectly bonded or imperfect). In addition, the results also depend on whether a

symmetric or asymmertic form of the self-consistent method is used.

In summary, the two self-consistent models, summarized before, will be tested against
the exact calculations of effective moduli carried out by computer simulations of very

large Delaunay networks.

 

\r1. 4‘ M
(afifivflW/w M
1.3V?
fik‘W @

Ir . Q
«bal/'4' .-
57/29 ..
56%»? E M
h" evol-

' 9,
sron analogy.

78

4.4.2 Periodic boundary conditions

Periodic boundary conditions have been used in many fields of science and proven to
be an excellent way to represent various problems of infinite nature. For Delaunay net-
works, it means that the geometry is repeated with periodicity L in x1 and x2 directions.
One window of Delaun-ay network with periodic boundary conditions may be generated as

follows:

First we generate n Poisson points from a planar homogeneous Poisson process of in a
square L x L. Then, we duplicate these Poisson points of the square into other two squares
one on the right and one on the left. The resulting three squares are duplicated twice in the
x2 direction: above and below. Thus, we end up with a total of 9 sub-squares forming a
new big square of size 3Lx3L and 9n Poisson points. Then, the process of generating a
Voronoi graph along with its dual Delaunay graph, as it was described before, is applied
over the new big square. Thus, generated Delaunay graph (as shown in Fig. 4.10) yields
in the middle square, an approximately squared-shaped part of the Delaunay network,
whose boundary BB is shown in bold lines; this is called a periodic Delaunay network. It
is next used for the purpose of calculating the effective moduli of the medium. It is of
importance to note that periodicity is meant here not only in the geometry of the network
but also in the physical nature of the boundary condition. That is, physical periodic

boundary conditions may be specified as

ui(.§+[~.) = 11, (1g) +€ijxj Vite BB

1(8) = —t(8+é) V-EE 58

(4.25)

where 8,]. is the macroscopic strain, ti is the traction on 88, and L. = Le, with e‘. being a

~45:

unit base vector.

79

Periodic boundary conditions represent a certain combination of the essential and the
natural boundary conditions given by (2.13) and (2.18), since the periodicity of the net-
work imposes also a periodicity in the forces acting between the neighboring vertices.
Thus, the effective elastic moduli of a network with periodic boundary conditions must lie
between those under essential and natural boundary conditions; actually they are closer to
those obtained under the essential boundary conditions. Hence, simulations of periodic
Delaunay networks give the closest estimates of the effective properties in a computer

simulation.

8O

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Fig. 4.10 A periodic Delaunay network generated by duplication of Poisson points of one

square into other 9 squares.

81

4.4.3 Numerical Solution Implementation

The method of solution is described here in a general sense so that it can be applicable
to any network regardless of its geometry or a model of interactions between grains used
along with it. The procedure to calculate the effective moduli (K, [1) of the generated peri-

odic Delaunay network is conducted in the following steps:

1) We generate one realization (specimen) of the Voronoi and its dual Delaunay networks
(note: a large size of network gives a better presentation of an infinite medium). Once a

realization is created we conduct the following operations on it:

i) Random assignment of the two phases type 1 (and 2) to all network’s vertices according

to a given volume fraction C1 (and C2) without any spatial correlations.

ii) Subjecting the boundary of the network to a uniform in-plane extension, corresponding
to a macroscopic strain 811 = 822 , and then calculating the total energy U(l) set in the net-
work as a sum all the energies of all its edges (bonds) which are modeled as linear springs.
It is important to note that the boundary edges are to be calculated only once, although
they show up twice in Fig. 4.10. Accordingly, bulk modulus of this realization can be

extracted from the energy calculations.

iii) Subjecting the boundary of the network to a uniform extension in x1 direction and a
uniform compression in x2 direction, corresponding to a simple shear 812, in a coordinate
system rotated by 45°, and then calculating the total energy U0) as before. This assign-
ment of boundary condition yields the shear modulus of the realization directly.

iii) Noting that the energy of a two-dimensional linear elastic continuum of volume V (=

L2) is

82

v 1
U = 2[wait-1317+ 211(8080- — iaiieifl] (4.26)

By equating the energy of the discrete model to that of continuum model, we calculate the

(area) bulk modulus as

20
K = “I (4.27)

2
W811 + 822)

 

Similarly, the (area) shear modulus can be calculated as

20
11 = __(2>_2 (4.28)
V9512)

It may be recalled from 2-D elasticity that bulk and shear moduli can be written as [8]

E E

=2(1-v) ”2 2(1+v) (429)

K

where E is the area Young’s modulus and V is the area Poisson’s ratio which can be

expressed in a 2—D plane as [67]

 

v = K_ II (4.30)

K + p
2) Carrying out steps i) through step iii) in Monte Carlo sense for a sufficient number of
realizations, say it, the effective bulk and shear moduli of the medium are obtained in the

sense of ensemble averaging which may expressed as

ff 2 K1 2 “t
[(e : (K) : :7:— “eff = (u) = :71... (4.31)

83

The outlined procedure of calculations is carried out via computer simulations where
three computer programs were developed and written in FORTRAN for this purpose.
While, one program is used to generate the geometry of the periodic Delaunay network, a
second one does the mechanics part of calculations where the energy of the discrete model
is minimized using a conjugate gradient algorithm [68]. As a result, the effective moduli
can be extracted. In the case of a regular network, the third program is developed to do

both tasks of the first two programs simultaneously.
4.5 Numerical Results

As stated earlier, we are going to calculate the effective elastic moduli of two systems:
(1) periodic Delaunay network and (2) periodic regular triangular network. The main
focus is given to Delaunay medium under two cases: (1) homogeneous case and (2) heter-
ogenous case. Two sets of numerical results are obtained corresponding to: (1) the central
interaction between grains and (2) the general model with both central and noncentral

interactions are considered.

4.5.1 Central Interaction Model
i1 Homogeneous medium:

It is a case of networks representing a granular media which are made of entirely type 1
(or type 2) vertices, i.e a homogeneous medium. Poisson’s ratio in such a homogeneous
case is found to be 0.38. The spring constant k of an edge of length 1 has been adapted
according to a series model (k (l) = 1'1). This observation motivates an investigation of V
for a Delaunay network in which springs are assigned according to a more general rule k
(1) = 1", where the exponent p is the same for all springs. Fig. 4.11 illustrates the depen-
dence of V on p from which we observe that the maximum Poisson’s ratio is reached when

p = -1. This value makes much sense for a spring model, since it assumes that the stiffness

84

of a spring is inversely proportional to its length. Therefore, it will be adapted throughout
all calculations. Since in the central spring model, the only input parameter is k (spring
constant), a special isotropy (here, V = 0.38) is yielded. For a homogeneous triangular lat-
tice Poisson’s ratio obtained from numerical simulations is 1/3, as it can be verified from a

unit cell spring model.

' 0.4 _

 

0.3-

3 0.2—

0.1—

 

 

0.0

 

Fig. 4.11 Dependence of the Poisson’s ratio v on p for a one-phase Delaunay network

with all spring constants assigned according to k :19.

ii H r n m i m

It is a non-trivial case where two phases of a granular medium are mixed at random.
Calculations of effective moduli are performed for c2 ranging from 0.0 to 1.0 varying by
0.1 increments for several contrasts; 10,100,1000,10,000. For contrast equal to 10, results
of computer simulations varying continuously are shown in Fig. 4.12 for bulk modulus

and in Fig. 4.13 for shear modulus (indicated by black triangles). The results of various

85

self-consistent models for (N = 2, i.e two-phase medium) with aspect ratio s = 2 are shown
also on the same graphs. We tried several ratios and, SCA-S with s = 2 gave the optimal
approximation to Delaunay network results for the entire range of contrast (stiffness ratio).
On the other hand, SCA-A method seems to be too ‘stiff’ for Delaunay networks but good
for regular triangular networks. Results of regular triangular network of approximately
the same size, indicated by black circles, show that the effective stiffness at c2 between 0.0
and 1.0 are stiffer than those of the Delaunay networks. Also, we observe that the random
geometry leads to a softening relative to a regular structure in the two-phase medium case.
The same observation was made in the case of one phase medium [50]. In Figs. 4.12 and
4.13 we show results of circular inclusions with springy interphase where the tangential
stiffness (spring constant) of the springy interphase equals to 7. Due to the drawback -
lack of representation of numerical results — it will not be used beyond this point. Voigt
and Reuss bounds at contrast 10 are plotted in Figs. 4.12 and 4.13 for bulk and shear mod-
uli, respectively. Since they are so widely spaced, they shall be dropped also in any fur-

ther calculation.

Then we proceed to perform calculations for bulk and shear moduli (for Delaunay and
regular triangular networks) at contrasts 100, 1,000, 10,000 (see Figs. 4.14-4.16). Aspect
ratio (5 = 2) in SCA-S model continues to be the best approximation to our numerical sim-
ulation results for Delaunay networks. For regular geometry, simulation results fit very

well the self-consistent model of circular inclusions with perfectly bonded interfaces.

On the other hand, the results of Delaunay networks are used to test the quality of the
self-consistent approximations in two limiting cases: networks with a finite area fraction
of one phase being either holes or rigid inclusions. These two cases are characterized by a
critical point at which the system makes a transition: from elastic to floppy - in case of per-

colation of holes, and from elastic to rigid - in case of percolation of rigid inclusions.

86

With our numerical simulation results of Delaunay networks, we approach both cases by
an assignment of a contrast increasing from 10, to 100, 1000 and 10,000, rather by direct
simulations of partially void or partially rigid networks. Thus, Fig. 4.14 shows the results
for K normalized by K1 at contrast 100, Fig. 4.15 at contrast 1000, and finally, Fig. 4.16 at
contrast 10,000. Graphs for [.1 normalized by [11 are exactly the same and, hence, they are
not shown here. Clearly, since we are plotting the resulting effective K and u on a scale
where K/K2 is finite, it is at contrast 10,000, shown in Fig. 4.16 that we get the case of par-

tially void networks. This is a situation of percolation of holes.

The results of the three self-consistent models; SCA-S, SCA-A and circular inclusion
are also plotted in Figs. 4.14, 4.15 and 4.16. It is found that SCA-S ellipses model at
aspect ratio 8 = 2 continues to be the best approximation here except in the range 02 = 0.5
to 0.6 which is a situation of high fluctuations. On the other hand, the circular inclusion
model with perfect bonding remains an excellent approximation to the regular network
results except between same points: C2 = 0.5 to 0.6. In fact, the self-consistent model pre-
dicts the transition in response at 2/3 volume fraction of phase 2, while the simulations do
so around 50%. It is interesting to note that 50% is the critical probability of Bernoulli
percolation of sites on such networks [69]. Also, it is important to observe that while the
self-consistent model with circular inclusion is off by 16.66%, other effective medium the-

ories of composites with perfect bonding give percolation at 100% [70].

87

 

 

 

 

  
  
 
  

 

 

 

 

 

10.0
—— p. b. (circles)
9'07 — — p. b. (symmetric ellipses)
8 0— - ' p. b. (asymmetric ellipses)
' p. b. (springy interface)
7.0- + Reuss bound
+ Voigt bound
50.. A exact (random geometry)
2" 0 exact (regular geometry)
\ 5.0-
52
4.0-
3.0-
2.0-
/-’=
l .0
0.0 I I T I I I I l I U I I I T I I F? I I
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Fig. 4.12 Effective bulk modulus at contrast 10 for Delaunay and regular networks; along

with various self-consistent models and Voigt and Reuss bounds.

 

 

 

 

 

 

 

 

10.0
9 O — p. b. (circles) g
' — — p. b. (symmetric ellipses)
8 0- .— ' - p. b. (asymmetric ellipses)
' ' ' ' ' p. b. (springy interface)
7.0_ + Reuss bound /"
+ Voigt bound /.'
/ a
5.0— exact (random geometry) " 3
1" exact (regular geometry) /"
/.
E 5.0-i /."
4.0- "
3.0- '-
2.0- _ =
1.0 “‘ ‘ '
0.0 T I I l I l I l f I I I I I I I I
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Fig. 4.13 Effective shear modulus at contrast 10 for Delaunay and regular networks,

along with various self-consistent models and Voigt and Reuss bounds.

88

 

 

 

 

 

 

 

 

 

100.0

90.0— — p. b. (circles) /

'— — p. b. (symmetric ellipses) /
80-0“ — ' " p. b. (asymmetric ellipses) ' /
70 O A exact (random geometry) //
' 0 exact (regular geometry) f
60.0— . //
" /
53 /
50.0— .
E / /
. /
40.0 '7‘
/ /
30.0— . ,
/
. /
20.0- / /
/
10.0— ‘ '/ /_
+— = 3 .
0.0J.=$———l——r r I r I r I l I r I F I l I r

Fig. 4.14 Effective bulk modulus of Delaunay and regular networks at contrast 100.

 

 

 

 

 

 

 

 

1000.0
900.0- — p. b. (circles) . /
— - p. b. (symmetric ellipses) /
800.0- — - - p. b. (asymmetric ellipses) /
exact (random geometry) /
700.0- /
0 exact (regular geometry) 4
soo.o-_~ /
- / /
’3 500 o ' /
\ . — /
>2 - /
4oo.o— ./ /
/ ’/
300.01 / /
._ / /
200.0 /’ /
/
A - ‘ / z
0.0 I fri? r T if T r I I I y l , i —r i
00 01 0.2 0.3 04 05 06 0.7 08 09 10

Fig. 4.15 Effective bulk modulus of Delaunay and regular networks at contrast 1,000.

89

 

 

 

 

 

 

 

 

10000.0
9000‘}: _ p. b. (circles)
— — p. b. (symmetric ellipses) //
8000.0- — - - p. b. (asymmetric ellipses) /
exact random eomet
7000.0- ( g ry) /
0 exact (regular geometry) ’/
6000.04 . /
s H
2 5000.04 /' /

' ' /

4000.0- /
' /
3000.0— . 7
, /
2000.0- /
1000.0— ./
. /
4 /
0.0'Jh—fl—F—H—hl—H l l 1 w l u
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Fig. 4.16 Effective bulk modulus of Delaunay and regular networks at contrast 10,000.

 

90

The case of rigid inclusions is approached by plotting KIIK for a very high contrast of
10,000 on a scale where K1 is finite. Thus, Fig. 4.17 gives the case. of contrast 100, Fig.
4.18 contrast 1000, and, finally, Fig. 4.19 contrast 10,000. Simulation results of both
Delaunay and triangular networks showed that SCA-A provides now a better fit, with a
change in aspect ratio: 3 = 2 in Fig. 4.17, s = 4 in Fig. 4.18. However, the final case of Fig.
4.19 cannot be approximated by any single aspect ratio in the entire range of (>2; this figure
shows the case of s = 8, which was found to provide the best overall fit. The transition in

the self-consistent model of circular inclusion occurs again at volume fraction 2/3.

In all simulations carried out here, at any specific c2, we employed 20 samples with
1,000 vertices each. This was the best we could achieve given the computer resources
available to us. The strongest scatter was encountered at 02 = 0.2, 0.3, 0.4 and 0.5 in Fig.
4.19, which motivated a setup of 10 samples with 2,000 vertices each. As a result, only

the regular triangular network at c2 = 0.3 appears tolbe off.

 

 

 

 

 

 

 

 

1.0
0'9‘ —- p. b. (circles)
0.8- — — p. b. (symmetric ellipses)
\ — ° - p. b. (asymmetric ellipses)
0.7- A exact (random geometry)
\ 0 exact (re ular eomet
0.5— \ 9 9 W)
2 A
\20-5- \
\\
0.4— s“
0.3— ‘\
0.2- 8‘
A \.
. — \\
0 l . \\. \
\ \
0.0 r I r I r I r I r I . l , l I l I I
0.0 0.1 0.2 0.3 0.4 0.5 0.5 0.7 0.8 0.9 1.0

Fig. 4.17 The reciprocal-of effective bulk modulus at contrast 100.

 

 

 

 

 

 

 

 

1.0
0.9— \ —.- p. b. (circles)
\ . .

O 8 - - p. b. (symmetric ellipses)

' \ — ' - p. b. (asymmetric ellipses)

0.7— \.\ . A exact (random geometry)

\§ 0 exact (regular geometry)
\\

0.6— .\
2 \\
\_o.5- ‘\\
5a e\

A\ \
0.4— . \
\ \
‘ \
0.3— \0 \
s
\
\ \
0.2- \
. . \
0.1- \ \\
A ‘ \ . 2 \ \ \
000 I I I l I I I I I l I I T- I I ~? I ? I
0.0 0.1 0.2 0.3 0.4 0.5 0 5 O 7 0.8 O 9 1 O

 

 

   
  
   
  

/ p. b. (circles)
------ p. b. (symmetric ellipses)
- ------ p. b. (asymmetric ellipses)

A exact (random geometry)

 

. exact (regular geometry)

 

 

 

 

Fig. 4.19 The reciprocal of effective bulk modulus at contrast 10,000.

 

92

4.5.2 Combined Central and Angular Interactions Models

So far, all calculations were carried out based on a central interaction between grains.
In this sub-section, results are obtained based on a consideration of two types of interac-
tions; central and angular. Hence, the general spring model introduced into Delaunay net-
works in such a way that in addition to each edge being a normal spring, an angular spring
connects each two neighboring edges to model angular interactions. Therefore, now, a
new parameter, r = knlka ratio is needed to be specified, kn is normal spring stiffness and
k3 is angular spring stiffness. As mentioned earlier, this model can represent a medium of

general isotropic elastic properties (i.e. two independent elastic constants, here K and It).

On the other hand, inclusion of the angular (shear) interactions in a Delaunay network
needs more efforts since computer time needed to run a single case becomes very large. In
addition, a higher number of realizations is essential to kill fluctuatibn in output results
especially for higher contrast. Due to limitation in space and time on any computer facil-
ity available to us, all calculations are performed at contrast 10 (note here that we used
both same contrasts in normal and shear interactions). To examine the effect of a change
in r on the effective elastic moduli of Delaunay networks, six values of r are chosen; 0.1,
1.0, 10, 50,100. When r = 00, the model reduces to the central spring one and all pervious
results are valid, a case which was verified here, while r = 0 means that only pure shear
interactions (no central interactions) exist between grains. In latter case, Poisson’s ratio is
found to be equal -1. Hence, the present model can predict effective properties of a
medium whose Poisson’s ratio ranges from -1 to 0.38. Proceeding to effective moduli cal-
culations, Fig. 4.20 shows the results of K normalized over K1 for a range of volume frac-
tion (C2) ranges from 0.0 to 1.0 with an increasing increment of 0.1. Examining these
curves one can observe that an increase in r parameter leads to a softening in effective bulk

and shear moduli. This observation is consistent with that one might expect - introduc-

93

ing angular springs into a network increases its stiffness. Fig. 4.21 shows the results of
effective [1. as a function of C2 while same behavior continues to be observed here with
only one exception - a high fluctuation in results at small values of r (= 0.1, 1.0) is
observed. These cases corresponds to a network which is very strong in shear but very

weak in normal interaction.

To present the results in form of two other independent elastic constants, E, V are calcu-
lated from K and It and plotted versus volume fraction (C2) as shown in Figs. 4.22 and
4.23, respectively. Effective Young’s modulus is decreasing as r increases (the same
behavior as before), but the problem of high irregularity at small values of r has disap-
peared. For effective Poisson’s ratio, the variation with respect to volume fraction (C2)
appears to be very weak near the extreme values of r (i. e. at r = 0.1 V has approximately
the same value of 0.98). At r = 00, (shown by the curve on top) we recover results
obtained by using the central spring model, which possesses also a weak variation. At
moderate values of r, a stronger variation which takes a form of a concave-up curve can be
seen for the cases at r =1, 10, 50, 100. It is of interest to note that it is approximately at r
= 20 when V equals 0. A note is made here of an investigation by Snyder et a1 [71] which

reports the dependence of V on the volume fraction for continuum-type composites.

Now we examine the best approximation among the three self-consistent models, SCA-
S, SCA-A and circular inclusion to our simulation results at three cases; r = 10, 50, 100.
Since these ratios represents practical cases of materials, so they are chosen here. Pure
shear or close to it is of a pure academic interest and does not posses any practical basis,
hence, it is ignored here. For such a case, we could not find any self-consistent model
among the three which may fit our results, which strengthen our decision - that it should
be ignored. Fig. 4.24 shows the simulation results of K for Delaunay network (normalized

over K1) along with the three models of self-consistent method. It is the circular inclusion

94

with perfect interface model that fits our results excellently (see Fig. 4.24). In the con-
trary, for both cases r = 50, 100, as shown in Figs. 4.25 and 4.26, the SCA-A inclusion
model is found to be the best approximation to its corresponding simulation results,
respectively. One should remember here that at r = 00, which corresponds to the case of
using the central spring model only, it was found that the SCA-S model was the best

approximation to the simulation results.

95

 

 

 
 
  
  
  
   
   

rs 0.1
r-1.0
r- 10.0
r= 50.0
r- 100.0
r- lnflnlty

liiiii

 

 

 

 

 

 

Fig. 4.20, Effective bulk modulus of Delaunay network at contrast 10 for different

values of r.

 

16

 

 

 

 

 

 

 

 

Fig. 4.21 Effective shear modulus of Delaunay network at contrast 10 for difi‘erent

values of r.

ll-

96

 

 

 

 

 

 

4o 2
-e— r-0.1 a
d -8— ram
+ r-10.0
304
+ r-50.0
—9i(- ”100.0
+ r-lnflnlty
$30-

10-

 

 

 

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 . 1.0

Fig. 4.22 Effective Young’s modulus of Delaunay network at contrast 10 for different val—

ues of r. ..

 

1.0 l 1

 

 

 

 

 

 

 

 

 

 

.04 -«
-0.6 -‘
.0.8 £11 8 3 48 S 8 C} G #8 8 DJ
' I I r T vfi T fi'T t T r T t I r T r I“ r
0.0 0.1 0.2 0.3 0.4 0.5 0.3 0.7 0.8 0.9 1.0
C

Fig. 4.23 Effective Poisson’s ratio of Delaunay network at contrast 10 for different

values of r.

97

 

 

10
9 .I I Slmulatlon

J -— P. e. Clrcle
a 7 -— SCA-S .//

 

 

 

 

 

 

 

Fig. 4.24 Effective bulk modulus from numerical simulations and self-consistent models

results for r = 10.

 

 

 

9 _] I Slmulatlon ’77
~ P. a. Circle
A I

a —— sacs / /

 

 

 

 

 

 

 

Fig. 4.25 Effective bulk modulus from numerical simulations and self-consistent models

results for r = 50.

98

 

 

9 .. I Simulatlon
P. B. Clrcle

SCA-S x
,/

 

 

 

 

 

 

 

 

Fig. 4.26 Effective bulk modulus from-numerical simulations and self-consistent models

results for r = 100.

CHAPTER 5

EFFECTIVE CONDUCTIVITY OF HETEROGENOUS INTERFACES
IN FIBER-MATRIX COMPOSITES

5.1 Introduction

The transition zone between matrix and inclusion (or fiber) in inclusion-matrix com-
posites usually forms due to bonding processes during manufacturing [3]. Diffusion,
nucleation and chemical reactions are some examples of these processes. Thus, it is a
region where the inclusion and matrix grains are chemically and mechanically combined.
Two terms most commonly used to describe this region are: interface and interphase.
When the transition zone is a two dimensional boundary separating matrix and inclusion,
it is called “interface”. Whereas, when it is assumed to be an additional phase between
matrix and inclusion or a third layer, it is called “interphase”. Most frequently, theoretical
studies in this area represent this zone as a two dimensional layer or as a three dimensional
region of certain microstructure, i.e. an interphase [3], which has distinct properties. Inter-
faces in composites influence the local fields and their effective properties. Therefore,
composites with inhomogeneous interfaces have attracted attention of several researchers
recently, i.e. see [72]. Works in this field include studies of stress fields due to thermal [5,
6, 7] and mechanical loadings and the evaluation of effective elastic moduli of composites
with inhomogeneous interfaces [8 9, 10]. However, in these studies the interphase region
is assumed to be isotmpic with Young’s modulus, as an example varying linearly, as a

power law, or as a polynomial, and in general Poisson’s ratio is taken as constant for sim-

99

100

plicity purpose. Anisotropy of the interphase is considered in [11, 12] where the anisotro-

pic constants vary as a power law.

All of these works were based on effective medium theories and no consideration was
given to the microstructure of interphase. In the present study, this complicated problem
is approached in a new way which is hinged on a discrete random model for the micro-
structure of interphase. The microstructure is represented as a zone of two randomly
interpenetrating phases with statistics distribution to be linearly along its thickness. We
consider two specific models of the random microsu'ucture: a fine-grained model with a
topology of a random chessboard, and a coarser-grained model with a topology of a
Voronoi planar tessellation, whose cells are occupied by either one of the two phases.
Although we mention these two models for random microstructures, our calculations will

performed based on the first model only (i.e. random chessboard with square pixels).

Next, the window concept is used as a “Representative Volume Element (RVE)” taken
from the interphase region. In order to calculate the effective properties of a window, the
same idea adapted in previous chapters is implemented here. Namely, at a given resolution
5 of the window, we apply two types of boundary conditions: deformation-controlled and
stress-controlled on the boundary of the window. Equivalency setup between energy of
two systems - discrete model for window’s microstructure and a meso—continuum approx-
imation - yields two bounds on the effective properties of the interphase. Results showed
that properties of interphase, are locally anisotropic and are given as a continuous differ-
entiable tensor function which may be approximated by polynomials. On the other hand,
it is important to note here that there is an alternative method, based on a “star-triangle

transformation” (Y --> A) for calculating an effective conductivity of a lattice due to Frank

101

and Lob [73]. However, our own program appears to be entirely sufficient and up to the
task of computing all the components of Keff of the interphase under specific - essential or
natural - boundary conditions. These results are fed into an effective medium theory to
predict the effective transverse conductivity of a composite with random distribution of
inclusions within the matrix [69]. In particularly, for the sake of simplicity, the Composite
Cylinders Assemblage model (CCA) [12] is employed in the present study to achieve this
purpose. A special feature of our study is its linkage of behavior at various length scales
of the composite body. While random field models of microstructure are established,
computer simulations are used to calculate effective properties of interphase at meso scale,
then, the CCA model is employed to calculate the effective (macroscopic) properties of

composites.

Finally, the governing equations for conductivity problem in 2-D are the same as for
out-of-plane elasticity, so that whatever outlined above is applicable also to an axial shear

modulus of a unidirectional composite.

5.2 Problem Statement

The main objective of the present study is to compute the effective transverse thermal
conductivity of the inhomogeneous interphase zone in a fiber-matrix composite. The inter-
face is explicitly taken as a finite thickness zone of two component phases. A particular
unidirectional fiber reinforced composite with continuous and aligned fibers, embedded in
a matrix, which are in shape of circular cylinders is considered. Each fiber-matrix inter-
face is a region of two randomly mixed phases, matrix and fiber (see Fig. 5.1). Both fiber
and matrix are assumed to be homogeneous and isotropic, that is, they are described by

two constant isotropic conductivity tensors. A random fluctuation is present in the com-

102

posite from two sources: (i) the random mixture of both fiber and matrix cells in the inter-

phase region and (ii) the random distribution of fibers in the matrix at macro level.

 

Fig. 5.1 Sketch of a transverse cross-section of a unidirectional fiber-reinforced compos-

ite with heterogeneous interphase.

To account for such randomness, four different length scales are introduced here:

i) 5micro = 1 - scale of heterogeneity,
ii) Simerphase = (b — a)/d - scale of the interphase,
iii) Smeso = L/d - scale of window taken from interphase,

iv) 5mm - macroscopic dimension of the composite body.

For a homogeneous isotropic interphase undergoing a steady state of heat conduction,

the constitutive law given by Laplace’s equation is

103

qt = —K..T (5.1)

where,

qi- heat flux,

T’j- temperature gradient,

Kij- thermal conductivity tensor in 2-D, given by

Ki. = K[1 0] where K = constant (5.2)
’ 0 1

but since the interphase region is heterogeneous and a random mixture of two phases, the

problem then turns to find an effective constitutive law analogues to equation 5.1 as
a; = 43” T J- (5.3)

where I_(_"”’ff is the effective conductivity of interphase, over-bars indicates volume averag-
ing. Calculations are conducted over a window of size 5. Postulating that the window is
equivalent to a RVE taken from the interphase, helps in finding Eeff. But since it is
obtained when window size goes to infinity (5—> 00), then, it is impossible in our problem
to be predicted by a unique value. Instead, compatible with our observations made in
chapter 3, only bounds on Keg may be obtained. The micro-meso linkage (i.e. between
discrete variation of conductivity at micro level and meso-continuum approximations)

will be performed via computer simulation. Details are provided in the next section.

Then, the two bounds on interphase’s effective conductivity at meso level will be used
to estimate the effective conductivity of the composite containing such interphases.
Here, we use a continuum mechanics approach to establish a linkage between meso and

macro properties. CCA model is utilized for this purpose.

104

5.3 Solution Procedure
5.3.1 Passage from the Micro to Meso Level

We employ a cylindrical coordinate system r, 0, 2 such that the z-axis is parallel to the
axes of the fibers and the r—9 plane is a transverse plane which becomes the plane of our
2-D problem (Fig. 5.2- a). The random mixture of two phases (m—matrix, and f-fiber) in
the interphase is represented by a random chess board (Fig. 5.2-b) or by a Voronoi mosaic
(Fig. 5.2-0). In both models cells are occupied at random by either one of the phases
according to a so-called indicator function:

1 if re V
X(E,C0) = r = (ne) we 9 (5.4)

0 if r e V
where: Vs- domain occupied by a phase, 5 = m, f.
Q- sample space,

(1)- one realization of the random interphase.

An established model of description of random media [74, 75] is used here, whereby B
= {B((0); to e Q}is a random microstructure, which was presented in chapter 2. Two
typical realizations of B((t)) for both models of the interphase are shown in Fig. 5.3- b and
c. In what follows we characterize X in terms of a probability distribution P {x (r) }

which is assumed to satisfy the end conditions
P(x(a0) =1) = 0 and P(X(b0) =1) = 1 (55)

where a0 is the radius of fiber, while b0 is the radius of interphase. The fiber and matrix

phases are both taken as locally isotropic and hence equations given by (5.5) implies that

105

(”1)

Kg) = 8‘.ij atr<a0, and [fig-m) = 8in atr>b0.

   

l.-.

C)

Fig. 5.2 Sketch of a fiber-interphase-matrix system, and two models of microstructure.

The variation in properties of interphase at micro level is piecewise-constant. Thus, the
highly heterogeneous interphase su'ucture — discrete (non-differentiable) and random -
poses a major obstacle in the analysis. The meso-continuum approximation, K ((1.0 (r, 0) ,
where i denotes the interphase, provides a solution for such a problem. The meso contin-
uum is obtained by treating the window as if it is an RVE characterizing the material
response K 5.0 (r, 0) as a transformation from the temperature gradient VT = [3 into the
heat flux 2], or vice versa. Hereinafter, overbears indicate volume averages over the vol-

ume (i.e. area) of RVE.

106

In order to define the effective conductivity of the window (RVE), the same methodol-
ogy as the one used earlier in the stochastic finite element problem is implemented here

based on a computer simulation. For sake of clarity, it is briefly summarized here as:

i) Application of two types of boundary conditions on window’s boundary:
a) essential (displacement-controlled): T = Bjxj ,
b) natural (stress-controlled): qknk = aknk.

where, T- temperature, B- temperature gradient, Ei- heat flux, and n- unit normal on bound-

ary.

ii) Energy calculations based on the discrete model under both (i-a) and (i-b) conditions for
one realization taken from the random microstructure of interphase. As mentioned earlier,
energy equivalence principle yields two values of conductivity; K; (0)) and

152(0)) = (R; (0)) )1, respectively.

iii) Statistical analysis (in Monte Carlo sense) over a sufficient number of realizations.
Thus, two (upper and lower) bounds on the effective conductivity of interphase given by
either (If (1')) or (K (r)) are obtained, where ( ) denotes ensemble averages. It turns
out that both tensors are orthotropic as can be seen from Fig. 5.2- b and c by noting a gra-
dient in the radial direction. Indeed, the 09-component was found larger that the rr-com-
ponent. Three features of the present model in contradistinction to the commonly assumed
isotropic meso—phase models [72] are:

i) the tensors posses an anisotropic character,

ii) a sigmoidal transition curve without a kink at r = a0,

iii) an effectively larger thickness of the interphase: a = a0 -L/2 and b = b0 +L/2 instead of

3.0 and 00.

107

5.3.1 Passage from the Mesa to Macro Level

So far we obtained two bounds on the effective conductivity of the heterogenous inter-
phase. Now, admitting the interphase as a distinct zone, the two bounds are fed into the
CCA model of Hashin and Rosen [12] to find fie“ of the whole composite. In this model
the actual unidirectional fiber-matrix composite is represented as a set of composite cylin-
ders which completely fill the space, (see Fig. 5.3). Each composite cylinder consists of a
cylindrical fiber enclosed in two concentric cylinders, the inner one representing the inter-
phase and the outer one the matrix. In addition, each composite cylinder has a similar
geometry such that a/b and a/c are constant, where a, b and c are - in the meso-continuum
approximation - the factious radius of the fiber, the factious outside radius of the inter-
phase and the factious outside radius of the matrix (of the composite cylinder), respec-

tively.

 

Fig. 5.3 Sketch of a Composite Cylinders Assemblage model.

108

In CCA model, first we calculate the effective conductivity of fig cylinder by replacing
it by an equivalent homogeneous cylinder with unknown effective conductivity property.
To achieve that, we can apply either essential or natural boundary conditions on the
boundary of the cylinder; but as was found By Hashin and Rosen [12], both yield identical
results. Therefore, we choose to apply a uniform temperature gradient field B such that

the temperature on boundary of the composite cylinder is

T|r___c = B~c~ 0059 (5.6)

in which c is the radius of the composite cylinder. The governing equation for our problem
is heat equation which for a case of no coupling between mechanical and thermal fields in

steady state condition is
4m. = KZITJI: = 0 s = f.m (5.7)

for two homogeneous regions, of fiber and matrix. For inhomogenous locally anisotropic

interphase eqn. (5.7) is
qk.k = (195;)ka = 0 (5.8)

Solution of equation (5.7), which is a Laplace’s equation in a polar coordinate (k = r, l =

0), it yields the temperature fields in the fiber as

Tm =Amr'cose 0<r<a (5.9)

and in the matrix as

109

1‘“) = [A(m)r+B(m)/r] cose b<r<c (5.10)

For the interphase region, eqn. (5.8) takes the following form in the polar coordinate

2 (i) (i) (i) 2 (i) 2
g1+ [819’ +1K(‘)]§I laK'e 3T Kre a T KrO a T __

3,2 Br r ” (Br-'75 W4, r 8100+ r2 392 — 0 (5'11)

 

(i)
KT?

in which a S r .<_ b with the K?) (r) tensor being either (K! (r)) or (K; (r) ). Recogniz-
ing that they are both orthotropic i.e. K g)“: = K 5;)" = 0, and they are varying in the
radial direction only, we develop polynomial fits for K (:2 e, K (i) e, K (:2 n and K ,9) n in the

rr

form ofk0+k1r+k2r2+k3r3+

Plugging the polynomial fits of (say the upper bound) into (5.11) and using the separa-

tion of variables technique, its series solution is obtained in the form of

N N
T = [Am 2 dn(r—r0)n+B(‘) 2en(r—r0)n]cose (5.12)
n=0 n=0

which involves two new constants A“) and Ba), here i stands for interphase and N for the

number of terms in the series.

Finally, we evaluate the five unknown constants A“), A“), B0), A0“), BC“) by using the

boundary conditions

(J) (i)
T = T
U) (i) } at r = a (5.13)
q = q
To) = T(m)
(i) ( ) at r = b (5.14)
m
‘1 = q

110

T0") = B-c- c050 at r = c (5.15)

where we note that

(3)2143)

(S)
q 8r

= _K s =f,i,m (5.16)
This permits the determination of the effective conductivity Eeff from equating the flux
on the boundary of our composite cylinder (10") (c) with the that on the boundary of the

equivalent homogeneous cylinder.

61W (6) = -Keffl3c089 (5.17)
which yields
ff ff ff K (m) (m) 3"")
K2, = 6HK‘ K" = — A —— (5.18)
l3 62

As proven by Hashin and Rosen [12], this results applies for a RVE (Representative Vol-
ume element) taken from the composite medium. Based on such a result, eqn. (5.18) is

truly our effective conductivity of the whole composite.

As pointed out in section 5.3, we note again that we actually conduct this procedure

twice, for K (1) equal to either (15;) or (Kg).

 

111

5.4 Discussion of Results

Based on the foregoing theoretical formulation, results are obtained for a composite
having a random chessboard interphase. Interphase under consideration characterized by
two parameters: contrast ratio Km/ K0“) = 10, and a0 = 1.0, b0 = 1.1 which implies that the
interphase thickness is 10% of fiber radius. We assume a linear distribution for fiber

grains along the interphase thickness as
P(x (5) = 0) = Ar+B (5.19)

where A and B are found from the end conditions (5.5) such that A = -l0.0 and B = 11.0.
Choosing the above linear distribution is due to the lack of experimental information.
Thus, we have a random mixture of two phases with concentration of fiber phase decreas-
ing gradually (here linearly) in the radial direction along the thickness of interphase. In
this study, the interphase is simulated by 100 squares (pixels). Three windows of 10x10,
20x20, and 40x40 pixels corresponding to 5 = 10, 20, 40, respectively from interphase are
considered. For each window of size 5 taken from a specific B((t)), numerical simulations,
i.e. solving a boundary value problem numerically, yield I_(e or (R2)-1 (R is the resisriv-
ity tensor) according to which boundary conditions; essential or natural is used, respec-
tively. For n radial points, ensemble averaging over N realizations in Monte Carlo sense
yields two bounds (Kg (r)) and (K? (r)) on the effective thermal conductivity of inter-
phase. We mention here, that although both Kre and Rre are, in general, different from
zero for any given realization 8(0)), the ensemble averages of these two components are
zero due to the assumption of axisymmetry of P (x (5) ) . Besides, (Kr) and (K99) were

found to be different which indicates that the interphase possesses an orthotropic property.

 

112

Averages for KH and R99 are shown in Figs. 5.4 - 5.6 for the three different values of 5

according to the following specifications:

 

 

 

 

' 5 =10 6: 20 8:40
11 8 13 13
a .995 .99 .98
b 1.105 1.11 1.12
N 1000 100 20

 

 

 

 

 

 

Polynomial approximation for each case is shown on the same graph along with its corre-
sponding output simulation results. These polynomial fits were developed by Wang [77],
and are reproduced here for the completeness of presentation. Examining these curves,
we can make the following observations:

i) difference between bounds becomes smaller as 5 increases. Also, as expected, results
obtained under essential boundary condition are stiffer than those under natural boundary
condition, where the former one is an upper bound on interphase conductivity while the
latter one is a lower bound.

ii) the polynomial functions shows a non zero slope at r = a - as it should not be the case -
this is due to a limitation in the fitting process itself. A gradual decrease through the
whole thickness of interphase can be observed where no jump or sharp gradient exists. It
is due to the nature of method used in this study that such a behavior occurs, i.e. the win-

dow concept is employed here. As a result, the thickness of the interphase becomes (b - a)

nor (b0 - a0).

iii) larger 5 leads to a higher anisotropy, while for the smallest 5 (= 10) the anisotropy is

almost non—existent.

113

iv) since scatter decreases as 5 increases, smaller number of realization (N) is needed for

high values of 5.

With the meso-continuum approximation at hand, we proceed to the calculation of
macrosc0pic effective moduli via the CCA model. Both bounds are used as an input for
K ,2) in this model, and thus will obtain two bounds on the macroscopic conductivity of
the whole composite [72]. In Fig. 5.7, we show only effective conductivity for the case 5
= 20, in which two graphs are plotted corresponding to (K?) and (EQ) as a function of
volume fraction of fibers. It is clear how both bounds on interphase conductivity produce

very close bounds on the macrosc0pic conductivity of the composite.

Finally, in this study, we showed how properties of heterogeneous interphase affect
composites’ overall properties. We used two different procedures of calculation at both
interphase and macro levels; random field modeling and computer simulations, and CCA

model, respectively.

U4

 
 
 
 
  
   

 

 

5 .
e
K00
1 Polynomial simulation of K e
4" rr
n
['1' 2 Polynomial simulation of K n
3 Polynomial simulation of K 8 9
3 .1-
: Pol nomial simulation of K n
. .1) y 90
U
2 u
1 a.
00
0 i i i i i l 4
0.995 1.00 1.02 1.04 1.06 1.08 1.10 1.105
r

Fig. 5.4 Variation of conductivity in anisotr0pic interphase under essential and natural

boundary conditions for 5 = 10.

115

    
  
   
  
   

 

 

l Polynomial simulation of K gr
11
KS]. 2 Polynomial simulation of Krr
4 «-
e
K09
rr
3 Polynomial simulation of K89
3 .
K90 4 P 1 ' 1 ' l ' f Kn
(I) o ynomra srmu ation o 99
i1
2 -.
1 ..
0 1 i i 'v 1 i i
0.99 1.00 1.02 1.04 1.06 1.08 1.10 1.11

Fig. 5.5 Variation of conduCtivity in anisotropic interphase under essential and natural

boundary conditions for 5 = 20.

H6

   
  
  
  
   

1 Polynomial simulation of K51.

n
2 Polynomial simulation of K rr
4 ._ 3 Polynomial simulation of K 8 9

, n
4 Polynomial simulatron of K 0 0

[KB
3“ rr

 

O 1 1 1 1 i 1 1

I I I

0.98 1.00 1.02 1.04 1.06 1.08 1.10 1.12

 

Fig. 5.6 Variation of conductivity in anisotropic interphase under essential and natural

boundary conditions for 5 = 40.

117

1.6 ..

   
 

1 Upper bound

1‘4 " 2 Lower bound

 

1.2 ..

 

1.0 ..

0.8 .-

 

l l l l l l
l I

0.0 0.1 0.2 0.3 0.4 0.5 0.6 1

Volume fraction f ,

Fig. 5.7 Conductivity of a composite with anisotropic interphases vs. volume fraction f of inclu-

sions for 5 = 20.

 

CHAPTER 6

CONCLUSIONS

Three engineering problems in micromechanics which have in common a stochastic
nature have been presented and methodologies of treatment were outlined. Due to vari-
ability in microstructures of special materials (heterogenous materials), effective (macro-
SCOpic) properties calculation should be admitted at micro level. Methodologies of
solution had the same basis: establishing a linkage between responses at different length
scales of the medium - micro, meso and macro - where randomness was incorporated
with the help of a random field theory. Introducing random discrete models at micro level
of the medium yields the meso-continuum models from which macroscopic responses can
be identified. In order to present clearly the basic ideas associated with the proposed
methods, analysis was conducted only in the context of linear elasticity (out-of-plane), and
thermal problems. However, for more complicated analyses such as nonlinear elastic, vis-
coelastic or elasto—plastic responses the underlying theme still holds. In this dissertation, a
finite element method based on micromechanics was proposed for solution of stochastic
boundary value problems. In another problem, a random discrete model, Delaunay net-
work, was employed to model random microstructures of granular media. Based on such
a modeling, effective elastic properties of the Delaunay network were calculated via com-
puter simulations. Finally, heterogenous interfaces in fiber-matrix composites were mod-
eled so that meso-continuum approximations were obtained and allowed the

determination of all effective properties of the composite.

1n the stochastic finite element problem, we developed a random continuum model

using a micromechanics approach, and then used it as input for solution of a specific

118

119

boundary value problem. Thus, we developed for the first time, a micromechanically
based stochastic finite element method. The link between the microstructure and finite
element method was provided by a so-called window which corresponds to a single finite
element cell. The window concept played a key role in the finite element formulations
since it yielded a specification of constitutive laws of a random meso—continuum. The
membrane example problem through which the method was presented here was a simple
but generic one - it illustrated the essential features of a micromechanics-based approach
in grasping the spatial heterogeneity of materials in stochastic finite element analyses. We
conclude that two random fields have to be considered in describing the microstructure
and bounding the response of the heterogeneous medium, with the choice of both fields
being dictated by the choice of a finite element mesh. Also, we note that the finer the finite
element mesh, the higher is the computational cost of a finite element solution, but the
lower is the computational cost of micromechanical specification of the material. Investi-

gation of both times as a function of the length scale 5, gave an optimal value of about 3.

It follows that by a relatively straightforward modification, a finite element program
can be made to incorporate a meso-mechanical input. Such an input may, in principle, be
provided in three different ways:

i) a direct computational simulation of the microstructure in each window (cell), repeated
11 times in a Monte Carlo sense for various realizations [29],

ii) an assignment of medium’s properties according to a random meso—continuum [27],
iii) an analytical derivation of the stochastic constitutive laws at the meso-scale. This

yields a finite element version based on homogenization theory.

The first part of calculations in this problem were obtained in the sense of method 1).
However, in light of method ii), a modification in the analysis was achieved through link-

ing the random fields of all finite elements together. Variance-covariance statistical

120

parameters were used to establish such a treatment by which a major advantage was
achieved. CPU time of analysis was reduced by two orders of magnitude, which gave the
method a much needed practicality and efficiency. Concerning results, covariance-based
finite element analysis for problem with mixed boundary conditions gave softer results
compared with exact analysis performed earlier. On the other hand, CPU time measure-
ments of our recent methods showed that it is comparatively more elegant than a direct
exact analysis of the whole domain. The ratio between former time of analysis to latter one

was found to be very small.

Finally, it is of need to extent some of stochastic finite element procedures to cope with

nonlinear and dynamic analysis as well as apply it to mechanics pr0blems like stochastic

continuum damage mechanics.

In the second problem, our objectives were: (1) to introduce a discrete modeling for
random microstructures of a special class of granular media; namely, Delaunay network
model of microstructures combined with two different spring models for central interac—
tions between particles at micro level and both central and noncentral interactions were
employed; (2) to calculate the effective elastic properties of two-phase granular media
with periodic boundary condition imposed on window’s boundaries which represents a
RVE taken from the random medium; (3) to identify a self-consistent model for a matrix-
inclusion composite which could adequately give bulk and shear moduli of two phase
Delaunay networks. Various such models were investigated against our numerical simula-

tions with the following results:

1) Using a central spring model:

i) For moderate and high contrasts, the ellipse SCA-S model provides the best fit to

numerical simulation results of Delaunay networks. The circular model is excellent for

121

regular triangular networks.

ii) For very high contrasts corresponding to the percolation of holes, these two models con-
tinue to be adequate except near the critical point.

iii) The self-consistent model, in a symmetric formulation, of inclusions with springy inter-
phases is only adequate in fitting the bulk modulus. The asymmeuic formulation of this
model was apparently investigated here for the first time, and with a negative result: the

model breaks down for volume fractions of inclusions exceeding 50%.

2i nrlrinm fr hcenlnnlinrin

i) Due to the high variability involved in responses especially at high contrasts, contrast 10
was only considered. To decrease the fluctuations, a larger size of Delaunay network as well
as a higher number of realizations are highly needed.

ii) As the ratio of stiffnesses of two springs (r = knlka) increases, softer effective bulk and
shear moduli are obtained. When r —-> 0, a case of pure shear interactions, the moduli were
found to increase. Effective Poisson’s ratio, corresponding to a homogeneous medium case,
in these two limiting cases were found 0.38 and -1, respectively. It was observed that Pois-
son’s ratio equals 0 approximately at r = 20. In 2-D plane, the range of Poisson’s ratio, as is
well-known, is between -1 and 1. One drawback of the present Kirkwood type model is its
inability to cover this whole range.

iii) For r = 10, the best fit to our numerical results among the three self-consistent models is

the circular inclusion model, while at r = 50, 100, it was the SCA-A with aspect ratio 4.

On the application side, results of this study will form basis for derivation of effective
constitutive laws of granular materials in a number of more complex situations where more
work is needed (see c. g. [76])such as:

i) Hertzian contacts need to be taken into account,

ii) infinitesimal elastic waves propagate under conditions of constant macrosc0pic strain,

122
iii) elasto-plastic deformations take place.

While, the three cases above still lack effective medium theories, numerical simulations
such as those presented here may be employed. Thus, for example, introducing angular
springs between contiguous edges, one could adapt the Delaunay network to simulate the
elastic-plastic transition and the incipient plastic flow in granular media - these results
should be used to identity the best self-consistent model and calibrate it in the fashion

demonstrated here.

The need for meso-continuum approximations arises in another problem of mechanics;
inhomogeneous interphases between fibers and matrix in fiber-matrix composites. Intro-
ducing several length scales, we predicted the effective transverse conductivity (and axial
shear modulus) of a unidirectional fiber reinforced composites with inhomegenous inter-
phases. We presented a new treatment of heterogenous microstructures of interphases in
which materials were functionally graded. The highly heterogenous nature of interphase
necessitated introduction of a meso-continuum concept. The choice of a length scale of an
RVE of the latter is quite arbitrary over a wide range of values. Unlike the passage to an
infinite scale, as commonly targeted in other problems of micromechanics, the finite size
of interphase precluded this possibility. As a result, two random fields bounding the effec-
tive response at any given 5 were introduced. This provided a micromechanics—based der-
ivation of the profile of material characteristics across the interphase. Both random fields
were averaged in an ensemble sense to provide two meso-continuum approximations for
input to the effective medium calculations, that were based here on the CCA model. The
convergent nature of resulting bounds on effective modulus Cerf, with increasing 5, sug-
gested a 5-dependent hierarchy of bounds. Such an observation provided a clarification

about the choice of 5 to be used.

123

BIBLIOGRAPHY

[1] Hashin, 2., “Analysis of composite materials,” A survey, ASME J. Appl. Mech. Vol.
50. PP. 481- 505 (1983).

[2] Nemat-Nasser, N. (Editor-in Chief) and Weertman, J ., Mechanics of Materials Jour-
nal, a special issue on Mechanics of granular materials, 16(1-2), pp 1-248.

[3] Drzal, L.T., “Composite interphase characterization”, SAMPE Journal September/
October, Vol. 7, pp. 7—13 (1983).

[4] Jasiuk, I. and Kouider, M.W., “The effect of an inhomogeneous interphase on the clas-
tic constants of transversely isotropic composites”, Mechanics of Materials, Vol. 15,
pp.53-63, 1993.

[5] Sottos, N.R., McCullough and Guceri, S. I.”Thermal stresses due to property gradi-
ents at the fiber-matrix interface,” Mechanics of Composite Materials and Structures (ed.
J. N. Reddy and J. L. Teply), pp. 11-20, ASME, New York (1989).

[6] Jayaraman, K. and Reifsnider, K. L., “Micromechanical stress analysis of continuous-
fiber composites with local material property gradients,” Achievements in Composites in
Japan and the United States (ed. A. Kobayashi), pp. 421-428, Korkon Shoin Co., Tokyo,
Japan (1990).

[7] J ayaraman, K. and Reifsnider, K. L., “Residual stresses in a composite with continu-
ously varying Young’s modulus in the fiber/matrix interphase,” J. Campos. Mater, Vol.
26. pp. 770-791 (1992).

[8] Theocaris, PS, The Concepts of Mesophase in Composites, Springer Verlag, Berlin
(1987).

[9] Sideridis, E., “ The in-plane shear modulus of fiber reinforced composites as defined
by the concept of interphase,” Compos. Sci. Tech, Vol. 31, pp. 35-53 (1988).

[10] Jayaraman, K. and Reifsnider, K. L., “Local stress fields in a unidirectional fiber
reinforced composite with non-homogeneous interphase region: formulation,” Adv. Com-

pos. Lett., Vol. 1, pp. 54-57 (1992).

[11] Jayaraman, K. and Reifsnider, K. L., “Local stress fields in a unidirectional fiber
reinforced composite with non-homogeneous interphase region: solution,” Adv. Compos.
Lett., Vol. 1, pp. 84-87 (1992).

124

[12] Hashion, Z., and Rosen, B. W., “The elastic moduli of fiber-reinforced materials,” J.
Appl. Mech., Vol. 31, pp. 223-232 (1964).

[13] Vanmarcke, E., Random Fields: Analysis and Synthesis, MIT Press, Cambridge, MA
(1983).

[14] Brenner, C.E., Stochastic Finite Element Methods (Literature Review),”, Internal
Working Report, No. pp. 35-91 (1991).

[15] Cambou, 3., “Application of first-order uncertainty analysis in the finite element
method in linear elasticity, Proceeding of the 2nd International Conference on Application
of Statistics and Probability in Soil and Structural Engineering, Aachen, West Germany,
Deutche Gesellschaft fur Grd-und ev, Essen, FRG, pp. 67-87 (1975).

[16] Hisada, T. and Nakagiri, 8., “Stochastic finite element method developed for struc-
tural safety and reliability,” Proceeding of the 3rd International Conference on Structural
Safety and Reliability, Elsevier, Amsterdam, pp. 395-408, (1981).

[17] Vanmarcke, E. and Grigoriu, M., “Stochastic finite element analysis of simple
beams,” Journal of the Engineering Mechanics Division, ASCE, Vol. 109(5), pp. 1203- 14
(1983).

[18] Shinozuka, M. and Yamazaki, F., “Stochastic finite element analysis: an introduc-
tion,” in: ST. Ariaratnam, G.l. Schueller and I. Elishakoff, eds., Stochastic Structural
Dynamics: Progress in Theory and Applications, pp. 241-291, Elsevier Applied Sci., Lon-
don, (1988).

[19] Baecher, GB. and Ingra, T. S., “Stochastic FEM in settlement prediction,” J. Geo-
tech. Engng. Div. ASCE, Vol. 107 (GT104), pp. 449-463 (1981).

[20] Benaroya, H. and Rehak, M., “Finite element methods in probabilistic sauctural
analysis: A selective review,” Appl. Mech. Rev. Vol. 41(5), pp. 201-213 (1988).

[21] Der Kiureghian, A. and Ke, J .B., “The Stochastic finite element method in structural
reliability,” Probabilistic Engineering Mechanics, Vol. 3, No. 2, pp. 83-91 (1988).

[22] Lui, W.K., Belytschko, T. and Mani, A., “Random field finite elements,” Int. J. Num.
Meth. Engng., Vol. 23. PP. 1831-1845 (1986).

[23] Deodatis, G., “Bounds on response variability of stochastic finite element systems,”
Journal of Engineering Mechanics, Vol. 116, No.3 pp. 565-585 (1990).

[24] Ghanem, R., and Spanos, P.D., Stochastic Finite Elements: a Spectral Approach,
Springer-Verlag, Berlin, 1991.

125

[25] Vanmarcke, E., Shinozuka, M., Nakagiri, S., Schueller, 6.1. and Grigoriu, M., “Ran-
dom fields and stochastic finite elements,” Structural Safety, Vol. 3, pp. 143-166 (1986).

[26] Deodatis, 0., Wall, W. and Shinozuka, M., “Analysis of two-dimensional stochastic
systems by the weighted integral method,” Computational Stochastic Mechanics (P.D.
Spanos and CA. Brebbia, Eds.), CMP-Elsevier, pp. 395-406 (1991).

[27] Ostoja-Starzewski M., Micromechanics as a basis of random elastic continuum
approximations,” Probab. Eng. Mech., Vol. 8(2), pp. 107-114 (1993).

[28] Ostoja-Starzewski, M., “ Micromechanics as a basis of stochastic finite elements and
differences - an overview,” Appl. Mech. Rev. (Special Issue: Mechanics Pan-America
1993), Vol. 46(11, Part2), Pp. 136-147 (1993).

[29] Alzebdeh, K. and Ostoja-Starzewski, M., "Micromechanically-based stochastic finite
elements,” Finite Elements in Analysis and Design, Vol. 15, pp. 35-41 (1993).

[30] Elishakoff, I., Probabilistic Methods in the Theory of Structures, J. Wiley & Sons,
New York (1983).

[31] Ostoja-Starzewski, M., “ Micromechanics as a basis of continuum random fields,”
Appl. Mech. Rev. (Special Issue: Micromechanics of Random Media), Vol. 47(1, Part2),
pp. 221-230 (1994).

[32] Reddy, J. N., Applied Functional Analysis and Variational Methods in Engineering,
McGraw-Hill, New York (1986).

[33] Oda, M., Nemat-Nasser, S., Mehrabadi, M. M., “A statistical study of fabric in a ran-
dom assembly of spherical granules,” Int. J. Numer. Anal. Methods Geomech, Vol. 6, pp.
77-94 (1982).

[34] Mehrabadi, M. M., Nemat-Nasser, S., Oda, M., “On statistical description of stress
and fabric in granular materials,” Int. J. Numer. Anal. Methods Geomech, Vol. 6, pp. 95-
108 (1982).

[35] Satake, M., “New formulation of graph-theoretical approach in the mechanics of
granular materials,” J. Mech. Mater. (Special Issue: Mechanics of Granular materials),
Vol. 16, pp. 65-72 (1993).

[36] Oda, M., “Initial fabrics and their relation to mechanical properties of granular mate-
rial,”, Soils and Foundations, Vol. 12(2), pp. 17-36 (1972).

[37] Statakc, M., “Constitution of mechanics of granular materials through graph repre-
sentation,” Theoretical Applied Mechanics, Vol. 26, University of Tokyo Press, pp. 257-
266 (1978).

126

[38] Jenkins, J.T., “Volume change in small strain axisymmetric deformations of a granu-
lar material, in: M. satake and J .T. Jenkins, eds., Micromechanics of Granular Materials,
Elsevier, Amsterdam pp. 143-152 (1988).

[39] Bathurst, RJ. and Rothenberg, L., “Micromechanical aspects of isotropic granular
assemblies with linear contact interactions, J. Appl. Mech. ASME, Vol. 55, 17-23

[40] Chang, CS, “Micromechanical modeling of constitutive equation for granular mate-
rial,”, in: M. satake and J.T. Jenkins, eds., Micromechanics of Granular Materials,
Elsevier, Amsterdam pp. 271-278 (1988).

[41] Chang, CS, “Micromechanical modeling of deformation and failure for granulates
with frictional contacts,” J. Mech. Mater. (Special Issue: Mechanics of Granular materi-
als), Vol. 16. pp. 13-24 (1993).

[42] Sidoroff, F., Cambou, B. and Mahboubi, A., “Contact force distribution in granular
media,” J. Mech. Mater. (Special Issue: Mechanics of Granular materials), Vol. 16, pp.
83-89 (1993).

[43] Rothenburg, L. and Bathurst, R.J., “Influence of particle eccentricity on microme-
chanical behavior of granular materials,” J. Mech. Mater. (Special Issue: Mechanics of
Granular materials), Vol. 16. PP. 141-152 (1993).

[44] Gourves, R., “Application of the Schneebli model in the study of micromechanics of
granular media,” J. Mech. Mater. (Special Issue: Mechanics of Granular materials), Vol.
16, PP. 125-131 (1993).

[45] Ostoja-Starzewski, M., Graph approach to the constitutive modeling of heterogenous
solids,” Mech. Res. Comm, Vol. 14(4), pp. 225-262 (1987).

[46] Davis, R. A. and Deresiewicz, l-l., “A discrete probabilistic model for mechanical
response of a granular medium”, Acta Mechanica, Vol. 27, pp. 69-89 (1977).

[47] Meijering, J. L., Phillips Res. Rep., Vol. 8, 270- (1953).
[48] Hatfield, J. C., Ph.D. Thesis, University of Minnesota (1978).
[49] Winterfield, P. H., Ph.D. Thesis, University of Minnesota (1981).

[50] Ostoja-Starzewski, M. and Wang, C., “Linear elasticity of planar Delaunay net-
works: random field characterization of effective moduli”, Acta Mech., Vol. 80, 61-80
(1989).

[51] Ostoja-Starzewski M. and Wang, C., “Linear elasticity of planar Delaunay networks.
Part II: Voigt and Reuss bounds, and modification for centroids”, Acta Mech., Vol. 84, 47-
61 (1990).

127

[52] Shechao Feng, Thorpe, M.F. and Garboczi, E., Effective medium theory of percola-
tion on central-force elastic networks, Phys. Rev. B 31, 276-280 (1985).

[53] Ostoja-Starzewski, M., Alzebdeh, K. and J asiuk, 1., "Elastic moduli of Delaunay net-
works: exact results versus effective medium theories," MEET'N'93—First Joint ASME-
ASCE-SES meeting, Charlottesville, VA, June 1993; in Homogenization and Constitutive
Modelling for Heterogeneous Materials (CS. Chang and J .W.Ju, Eds.) ASME-AMD, Vol.
166, pp 79-85.

[54] Ostoja-Starzewski, M., Alzebdeh, K. and Jasiuk, 1., “Linear elasticity of planar
Delaunay networks III: Self-consistent approximations”, Acta Mechanica, accepted for
publication.

[55] Thorpe, M.F. and Sen, P.N., “Elastic moduli of two-dimensional composite continua
with elliptical inclusions”, J. Acoust. Soc. Am., Vol. 77 (5), 1674-1680 (1985).

[56] Shechao F., Thorpe, M. F. and Garboczi, E. J ., “Effective -medium theory of percola-
tion on central-force elastic networks”, Phys. Rev. B, Vol. 31, pp 7276-7281 (1985).

[57] Garboczi, E. J. and Thorpe, M. F., “Effective -medium theory of percolation on cen-
tral-force elastic networks. 11. Further results”, Phys. Rev. B, Vol. 31, pp 276-280 (1985).

[58] Garboczi, E. J. and Thorpe, M. F., “Effective -medium theory of percolation on cen-
tral-force elastic networks. Ill. The superelastic problem”, Phys. Rev. B, Vol. 33, pp 3289-
3294 (1986).

[59] Day, A., Snyder, K., Garboczi E. and Thorpe, M., “The elastic moduli of a sheet con-
taining circular holes,” J. Mech. Phys. Solids, Vol. 40(5), pp 1031-1051 (1992).

[60] Roger, C. A., Packing and Covering, Cambridge Univ. Press, London (1964).
[61] Rhynsburger, D., Geographical Analysis, Vol. 5, pp. 133- (1973).

[62] Stoyan D. Kendall W.S. and Mecke J ., Stochastic Geometry and its Applications,
John Willey & Sons, New York (1987).

[63] Holnicki-Szulc, J. and Rogula, D., “Nonlocal continuum models of large engineering
structures,” Arc. Mech., Vol. 31(6), pp 793-802 (1979).

[64] Keating, P.N., “Effect of invariance requirements on the elastic strain energy of crys-
tals with application on the diamond structure,” Phys. Rev., Vol. 145, pp 637- 645 (1966).

[65] Kirkwood, J .G., “The skeletal modes of vibration of long chain molecules,” J. Chem.
Phys, Vol. 7, pp 506-509 (1939).

[66] Peiying, S., A comprehensive Exam Report, Michigan State University (1994).

128

[67] Thorpe, M.F. and Jasiuk, 1., “New results in the theory of elasticity for two-dimen-
sional composites”, Proc. Roy. Soc. London, Vol. A438, 531-544 (1992).

[68] Press, W. H. , Teukolsky, S.A., Vetterling WT. and Flannery B.P., Numerical Recipes
in Fortran, Cambridge University Press, New York (1986).

[69] Ziman, J. M., Models of Disorder, Cambridge University Press, 1979

[70] Jun, S. and Jasiuk, 1., Elastic moduli of 2D cbmposites with sliding inclusions - a
comparison of effective medium theories, Int. J. Solids Struct., to be published (1993).

[71] Snyder, K.A., Garboczi E.J. and Day A.R., “The elastic moduli of simple two-
dimensional isotropic composites: Computer simulation and effective medium theory,” J.
Appl. Phys, Vol. 72 (12), pp 5948-5955 (1992).

[72] Ostoja-Starzewski, M., Jasiuk, 1., Wang, W. and Alzebdeh, K., “Composites with
functionally graded interphases: meso-continuum concept and effective properties,” sub-
mitted to Journal of the Mechanics and Physics of Solids.

[73] Frank, D. J. and Lobb, C. J ., “ Highly efficient algorithm for percolative transport
studies in two dimensions,” Phys. Rev B, Vol. 37(1), pp. 302-307 (1988).

[74] Willis, J.R., “Variational and related methods for the overall properties of compos-
ites, Adv. Appl. Mech., Vol. 21, pp. 2-78 (1981).

[75] Sobczyk, K. Stochastic Wave Propagation. Elsevier, Amsterdam.
[76] Shen, H.H., Satake, M., Mehrabadi, M., Chang, CS. and Campbell, C.S., Eds.,

Advances in Micromechanics of Granular Materials, Stud. Appl. Mech., Vol. 31, Elsevier
(1992).

[77] Wang, W., unpublished notes, 1994.

 

rv- LIBRARIES

”lilillljlllll“

5133