. «1::
.1’( Y... .
I...

a h 9' f
eiiftfif.
5.9-2“?
flfg‘fi:
Jhtsrzgtf

mgu“ :. ;
iAf 4.: 1.7"
flcfihfil

.3

5.8.4.? .0.

r»... an" . in;
gamma. 1.
a}...

§ ,
:-

.!. . or

as»). 5..

.3

y...

‘ : i . . .2

«Hub...» lam! . . V ‘. . page. :5

w. .04.. . 1v..- 233$)

unu.1,,.« . ‘ an». nun

w i in $10!. 3'
z W! 5...:

.3949
4.1;: ‘
.7 X5».li..11b
,1? {sins}. y
«‘

{:25 o. v \
, . .i
vi-.. . . . . $23.31;

3h!!!” . . . . . V , ‘ :1 tannin.

v»...~ul: - u ‘55:.
;y.nunnv.,‘.uaq xx
1...“; sl .xr...‘

h" 1
145:3

4»... {Mn
x. .2, . Runway...
., 2:3.

om"

¢

5
ill? #31: . i 0 v
«V5.3. h. u .t‘ 10.1%.4‘
1‘50» 1 u! $3. , , 3!:
.3... O... ’39.!” I. 5

“1:! mm“;

5.1 H3.~..+A.

..13.:Loh

0.1.5.0 \
a). )1

3.0 ~.v\ A

I'll. XIKAA‘! £~¢ :l

’- .vo..;ct!c 5|
Al}‘\

0...

JfimmNnfih...’

Wag? ‘ ‘ ‘ I V . . ‘
. , ‘ awfimimwrfirn ‘ ¢ .

.. . ‘ , __ , _ ‘ mm,

 

7)

lll’lli’lli! lllllllllllllilllllll’liill

3 1293 014101830

This is to certify that the

dissertation entitled
ANALYTICAL AND COMPUTATIONAL

APPROACHES

IN MECHANICS OF COMPOSITE MATERIALS
presentedby

Peiying Sheng

has been accepted towards fulfillment
of the requirements for

:94, D degree in ”8:24 a k {If (AL—«jiméé‘l cruj ) I

 

CZQDaQMméjébhuLQ
”a Ma Ma. Lam. W.

Major {tofessor

 

 

Date 7"”(6 251/7775-
J 7

MS U i: an Affirmative Action/Equal Opportunity Institution 0-12771

 

 

 

A. _V.___,_..._. —_-...u_—F—M_ I I

 

LIBRARY
Michigan State
University

 

 

 

PLACE N RETURN BOX to remove this checkout from you: record.
TO AVOID FINES Mum on or More data duo.

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU lsAn Aflimativo Wl—m-d—W Opponunity Institution
Wanna

 

 

ANALYTICAL AND COMPUTATIONAL
APPROACHES
IN MECHANICS OF COMPOSITE MATERIALS

By
Peiying Sheng

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

1995

ABSTRACT

ANALYTICAL AND COMPUTATIONAL APPROACHES
IN MECHANICS OF COMPOSITE MATERIALS
By
Peiying Sheng

This thesis consists of two parts. In the first part, we analyze elastic fields of a half
space having a spherical inclusion. The half-space is subjected to either a remote shear
stress parallel to its plane boundary or to a uniform shear eigenstrain in the inclusion. The
interface between the inclusion and the surrounding matrix is either perfectly bonded or is
allowed to slide (slip) without friction. We obtain an analytical solution using displace-
ment potentials in forms of integrals and infinite series.

In the second part, we study the effective moduli and damage formation in out-of-
plane elasticity (or, equivalently, two-dimensional conductivity) of matrix-inclusion com-
posite materials with either randomly or periodically distributed inclusions (fibers) by a
computational approach based on a finite difference spring—network scheme. Damage evo-
lution is simulated by sequentially removing/breaking springs (bonds) in this lattice in
accordance with the evolving state of stress/strain concentrations. The composite systems
are characterized by two parameters: stiffness ratio and strength ratio of both phases. In
particular we investigate the following aspects: basic classification of effective constitu-
tive responses, geometric patterns of damage, varying degrees of randomness of the incl-u-

sions’ arrangements, and mesh resolution of continuum phases.

CONTENTS

LIST OF FIGURES ................................................................................................. VI

CHAPTER 1 A SPHERICAL INCLUSION IN AN ELASTIC

HALF-SPACE UNDER SHEAR
1 INTRODUCTION ..................................................................... l
l .1 INTRODUCTION ....................................................................... l
1.2 PROBLEM STATEMENT .......................................................... 3
2 METHOD OF SOLUTION ......................................................
2.] ELASTICITY THEORY IN THREE-DIMENSIONS ...............
2.2 UNDISTURBED FIELD ............................................................. 7
2.3 HARMONIC POTENTIALS ...................................................... 9
2.4 SOLUTION ................................................................................. ll
3 RESULTS AND DISCUSSION ................................................ 28
REFERENCES ........................................................................... 53
APPENDIX ................................................................................ 57

CHAPTER 2 EFFECTIVE PROPERTIES AND DAMAGE
FORMATION IN RANDOM COMPOSITE
MATERIALS: A COMPUTATIONAL APPROACH

1 INTRODUCTION ...................................................................... 71
2 COMPOSITE MICROSTRUCTURE AND THE
MESO-CONTINUUM MODELS ....... 75
2.1 RANDOM COMPOSITE MODEL ............................................. 75
2.2 TWO SCALE-DEPENDENT RANDOM MESO-
CONTINUUM FIELDS ............................................................... 76

IV

3.1
3.2
3.3

3.4

4.1
4.2

5.1
5.2
5.3

FINITE-DIFFEREN CE MODEL ............................................. 8O

LAYOUT ...................................................................................... 80
MODEL FOR PLANE ELASTICITY ......................................... 81
APPLICATION OF THE MODEL TO

PLANE ELASTICITY ................................................................. 88
APPLICATION TO CONDUCT IVIT Y ....................................... 94
EFFECTIVE CONDUCTIVITY ............................................... 97
METHOD OF SOLUTION .......................................................... 97
RESULTS ..................................................................................... 99
DAMAGE AND FRACTURE SIMULATION

OF COMPOSITES ..................................................................... 105
METHOD OF SOLUTION .......................................................... 105
RESULTS ..................................................................................... 107
CLOSURE .................................................................................... 113
REFERENCES ............................................................................ 115

LIST OF FIGURES

CHAPTER 1

Fig. l Spherical inhomogeneity in a half—space.
Fig. 2 Stress owww) versus angle (p for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when P = 100 and 6 = 0° for perfect

bonding and remote shear loading case.

Fig. 3 Stress Owww) versus angle (p for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when I" = 2 and 9 = 0° for perfect bond-

ing and remote shear loading case.

Fig. 4 Stress owww) versus angle (p for different radii a = 02 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when P = 0.5 and 6%: 0° for perfect

bonding and remote shear loading case.

Fig. 5 Stress owww) versus angle (p for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when I‘ = 0.01 and 6 = 0° for perfect

bonding and remote shear loading case.

Fig. 6 Stress 6W) (OW) versus angle (p for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when P = 100 and 6 = 0° for sliding

and remote shear loading case.

VI

Fig. 7 Stress owww) versus angle (p for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when F = 2 and 0 = 0° for sliding and

remote shear loading case.

Fig. 8 Stress mew) versus angle (p for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when F = 0.5 and 0 = 0° for sliding and

remote Shear loading case.

Fig. 9 Stress owww) versus angle (p for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when F = 0.01 and 0 = 0° for sliding

and remote shear loading case.

Fig.10 Stress 6w versus angle (9 for perfect bonding and Sliding cases when 1‘ -.- 0,5 ,
0 = 0° and different radii a = 0.2 (solid line), 0.5 (dashed line) and 0.8 (dash-

dot line) for remote shear loading.

Fig.1] Jump in displacement [u‘p] versus angle (p for different radii a = 0.2 (solid
line), 0.5 (dashed line) and 0.8 (dashdot line) when F = 100, 0 = 0° for pure

sliding case and remote shear loading.

Fig.12 Jump in displacement ("<91 versus angle (p for different radii a = 0.2 (solid
line), 0.5 (dashed line) and 0.8 (dashdot line) when F = 2 , 0 = 0° for pure

sliding case and remote shear loading.

Fig.13 Jump in displacement [u‘p] versus angle (p for different radii a = 0.2 (solid
line), 0.5 (dashed line) and 0.8 (dashdot line) when F = 0.5 , 0 = 0° for pure

sliding case and remote shear loading.

VII

Fig.14 Jump in displacement [uq’] versus angle (p for different radii a = 0.2 (solid
line), 0.5 (dashed line) and 0.8 (dashdot line) when F = 0.01 , 0 = 00 for pure

sliding case and remote shear loading.

Fig.15 Stress Oxx(6xx) along z-axis for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when F = 100 for perfect bonding case and

remote shear loading.

Fig.16 Stress oxx(6xx) along z-axis for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when F = 2 for perfect bonding case and

remote shear loading.

Fig.17 Stress Oxx(6xx) along z-axis for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when F = 0.5 for perfect bonding case and

remote shear loading.

Fig.18 Stress Gxx(6xx) along z-axis for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when F = 0.01 for perfect bonding case and

remote shear loading.

Fig. 19 Stress (5M (6“) at several points along z-axis versus radius of inclusion a for
both perfect bonding (dashed lines) and pure sliding (solid lines) when F = 100

(remote shear loading).

Fig.20 Stress Gum”) along z-axis for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when F = 100 for perfect bonding and the

eigenstrain case.

VIII

ix

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

14

15

l6

17

18

19

20.

21

22

Fig. 23

Fig. 24

Effective conductivity, at window size L = 64.
Effective conductivity versus volume fraction, at window size L = 64.

Effective conductivity and Hashin second order bounds versus volume fraction,

Ci/Cm = 0.02 and 50.

Effective conductivity and Hashin second order bounds versus volume fraction,

Ci/Cm = 0.1 and 10.

Effective conductivity and Hashin second order bounds versus volume fraction,

Ci/Cm = 0.2 and 5.

Effective conductivity and Hashin second order bounds versus volume fraction,

c:‘/cm = 0.5 and 2.

Effective conductivity, at window size L = 128.
Effective conductivity versus volume fraction, at window size L = 128.
Two kinds of constitutive responses.

Damage pattern for random arrangement of inclusion,

ei/em= 0.5, Ci/Cm = 0.1.

Effective constitutive laws for the case of

29/6“: 0.5, Ci/Cm = 0.1.

Fig. 25

Fig. 26

Fig. 27

Fig. 28

Fig. 29

Fig. 30

Fig.3]

Fig. 32

Fig. 33

Damage pattern for random arrangement of inclusion,

e‘/e'“= 0.5, Ci/Cm = 0.1.

Damage pattern for random arrangement of inclusion,

9.92“”: 0.5, Ci/Cm = 0.1.

Damage pattern for periodic arrangement of inclusion,

29/5“: 0.5, Ci/Cm = 0.1.

Damage pattern for random arrangement Of inclusion,

eVe'“: 2, Ci/Cm = 0.2.

Damage pattern for random arrangement of inclusion,

e‘/e"‘= 2, Ci/Cm = 0.2.

Damage pattern for random arrangement of inclusion,

292'“: 2, Ci/C =02.

Damage pattern for periodic arrangement of inclusion,

e‘/e"‘= 2, Ci/Cm = 0.2.

Effective constitutive laws for the case of

21/8“: 2, Ci/Cm = 0.2.

Damage pattern for random arrangement of inclusion,

Fig. 34

Fig. 35

Fig. 36

Fig 37

Fig. 38

Fig. 39

Fig. 40

Fig. 41

Fig. 42

892‘“: 10, Ci/Cm = 0.1.

Effective constitutive laws for the case of

ei/em= 10, Ci/Cm = 0.1.

Damage pattern for random arrangement of inclusion,

e‘/e"‘= 10, Ci/Cm = 0.1.

Damage pattern for random arrangement of inclusion,

e‘/e"‘= 10, ci/c’"=0.1 ,L=127.

Damage map of effective constitutive laws, also showing responses of inclusions

and matrix.
Damage map of effective constitutive laws.
Damage map of damage pattern.

Effect of fluctuation of inclusion arrangement on damage pattern,

raj/em: 2, Ci/Cm = 0.2.

Effect of fluctuation of inclusion arrangement on stress-strain laws,

ei/emz 2, Ci/Cm = 0.2.

Effect of fluctuation of inclusion arrangement on damage pattern,

XII

Fig. 43

Fig. 44

Fig. 45

Fig. 46

Fig. 47

Fig. 48

Fig. 49

Fig. 50

e‘/e"‘= 2, chm = 5.

Effect of fluctuation of inclusion arrangement on damage pattern,

292'“: 0.5, Ci/c’" = 0.1.

Effect of fluctuation of matrix material on damage pattern,

ei/em= 2, Ci/Cm = 0.2.

Effect of fluctuation of matrix material on constitutive law,

e‘/e'“= 2, Ci/Cm = 0.2.

Effect of fluctuation of matrix material on damage pattern,

292’“: 0.5, Ci/Cm = 0.1.

Effect of fluctuation of matrix material on constitutive law,

e‘/e"‘= 0.5, c‘/C"‘ = 0.1.

Effect of fluctuation of matrix material on damage pattern.

ei/Em= 2, (9/6“ = 5.

Effect of fluctuation of matrix material on constitutive law,

e‘/e"‘= 2, Ci/Cm = 5.

Effective constitutive curve, (ii/em: 1, Ci/Cm = 0.1 and 20 samples.

XIII

Fig. 51

Fig. 52

Fig. 53

Fig. 54

Fig. 55

Fig. 56

Fig. 57

Fig. 58

Fig. 59

Fig. 60

Probability distribution function F (omax) for the case of Ei/emz 1.0,

Ci/Cm = 0.1.
Probability distribution function F (emax) for the case of
e‘/e"‘= 1, C‘/c’" = 0.1.
Effective constitutive curve, ei/emz 0.1, Ci/Cm = 10 and 20 samples.

Probability distribution function F (Gmax) for the case of

eVe’": 0.1, Ci/Cm = 10.

Probability distribution function F (Emax) for the case of

e‘/e"‘= 0.1, Ci/Cm = 10.

Effective constitutive curve. Ei/Emz 10, Ci/Cm = 10 and 20 samples.

Probability distribution function F (Gmax) for the case of

e‘/e"‘= 10. c‘/c"‘ = 10.

Probability distribution function F (E ) for the case of

max

93/2”: 10, chm = 10.
Effective constitutive curve, iii/8m: 1, Ci/Cm = 10 and 20 samples.

Probability distribution function F (0' ) for the case of

max

XIV

Fig. 61

Fig. 62

Fig. 63

Fig. 64

Fig. 66

Fig. 67

Fig. 68

Fig. 69

e‘/e"‘= 1.0, ci/c'“ = 10.

Probability distribution function F (emax) for the case of

e‘/e"‘= 1, C'/C = 10.
Effective constitutive curve, ei/em= 1o, Ci/c'“=0.1.

Probability distribution function F (0’ ) for the case of

max

si/em= 1o, (:i/cm = 0.1.

Probability distribution function F (emax) for the case of

ei/e'“: 10, Ci/Cm = 0.1.
Effective constitutive curve, (Si/em: 0.1, Ci/Cm = l and 20 samples.

Probability distribution function F (0’ ) for the case of

max

e‘/e'“= 0.1, Ci/Cm = 1.

Probability distribution function F (8 ) for the case of

max

e‘/e"‘= 0.1, ci/C’“ = 1.
Effective constitutive curve, Ei/emz 0.1, Ci/Cm = 0.1 and 20 samples.

Probability distribution function F (omax) for the case of

ei/e‘“: 0.1, Ci/Cm = 0.1.

XV

Fig. 70 Probability distribution function F (e ) for the case of

max
e‘/e"‘= 0.1, ci/C'“ = 0.1.
Fig. 71 Damage map, constitutive curves.

Fig. 72 Damage map, distribution function F (omax) .

Fig. 73 Damage map, distribution function F (Emax) .

XVI

CHAPTER 1

A SPHERICAL INCLUSION
IN AN ELASTIC HALF-SPACE
UNDER SHEAR

1. INTRODUCTION

1.1 INTRODUCTION

When an inclusion is present in a matrix and a loading is applied, the elastic fields
are disturbed in the vicinity of the inclusion. These stresses depend on a number of fac-
tors: the geometry and location of the inclusion, the elastic constants of the inclusion and
the matrix, the boundary conditions at the matrix—inclusion interface, the loading, and
other.

Inclusion problems have been addressed by many researchers. For a review of liter-
ature in this area see Mura (1987). Most of these studies, however, considered the cases
when the inclusion is embedded in an infinitely extended material and the matrix—inclu-
sion interface is perfectly bonded.

The elasticity problem of a half-space with a spherical (or spheroidal) inclusion or
cavity has been studied by several researchers. Among them, Tsuchida and Nakahara

(1970, 1972) solved the problem of a semi-infinite elastic body with a spherical cavity

subjected to a remote all-around tension or a uniform pressure on either the surface of cav-
ity or on the plane boundary. Tsutsui and Saito (1973) solved the problem of a semi-infi-
nite body containing a perfectly bonded spherical inhomogeneity under all around tension,
while Tsuchida and Mura (1983) considered a similar problem involving a spheroidal
inhomogeneity. Other related papers are due to Atsumi and Itou (1974), Tsuchida et al.
(1973), Tsuchida et al. (1982), and Yu and Sanday (1990). The problem of an ellipsoidal
inclusion, having same elastic constants as the matrix and subjected to dilatational strains
was solved by Seo and Mura (1979), and the same problem but involving a spherical
inclusion was studied by Mindlin and Cheng (1950) and Wachtman and Dundurs (1971).
All the above works were restricted to axisymmetric cases. Tsuchida and Nakahara
(1974), using a combination of Boussinesq, Neuber and Dougall displacement potentials,
solved an asymmetric problem of a spherical cavity in a half space subjected to either a
uniaxial tension or a pure shear. Recently, Jasiuk et al. (1991) solved the problem of a
half space containing a sliding spheroidal inclusion under either an axisymmetric remote
tension or an eigenstrain loading.

The corresponding two dimensional problems involving a half-plane include: a hole
under a uniaxial tension solved by Jeffery (1920) and Mindlin (1948), an inclusion sub—
jected to an eigenstrain loading considered by Richardson (1969), and a circular perfectly
bonded inclusion under a remote uniaxial tension studied by Saleme (1958) and Shioya
(1967). Also, Lee et al. (1992) addressed the case of a sliding circular inclusion in a half-
plane under either a remote uniaxial tension or an eigenstrain loading.

In this thesis, we investigate the joint effect of a traction free-plane boundary and a
sliding spherical inclusion embedded in a half-space. The inclusion is subjected to either
a remote pure shear stress parallel to the plane boundary or a pure shear eigenstrain in the
inclusion. The interface between the inclusion and the matrix is either perfectly bonded or
allows slip without friction (shear tractions are zero) while maintaining continuity of nor-

mal displacement and tractions. This study is related to earlier works involving a spheroi-

dal sliding inclusion under shear embedded in an infinite matrix (Jasiuk et al., 1987;

Sheng,1992)

1.2 PROBLEM STATEMENT

We consider a semi-infinite elastic material containing a spherical inclusion of radius
0, having different elastic constants than those of the matrix, as shown in Fig. 1. The load-
ing is either a uniform pure shear stress, parallel to the plane boundary, applied at infinity,
or an eigenstrain loading of shear type in the inclusion.

These loading conditions can be expressed as

o'er = —ny = P0 (x, y —> w)
(1)
SH=—£),y=8 (in Q)

t

it at
where on, 6),), are stresses, Sn, 8),), are e1genstrams, and p0 and e are constants. In
our notation we use a bar to denote quantities in the inclusion.
In the solution we use all three coordinate systems, Cartesian (x, y, z), cylindri-

cal (r, 0, z) and spherical (R, 0, (p) . The relations among these coordinate systems are

x = rcosB = Rsintpcose

y = rsin0 = RsintpsinG (2)

N
ll

Rcoscp

We let the origin of coordinates be at the center of the spherical inclusion and the
positive direction of the z-axis be downward and, without the loss of generality, we take

the perpendicular distance from the origin of the inclusion to the traction-free surface as

unity so that the plane boundary is located at z = -1.
The boundary conditions are as follows:

1. tractions at infinity (x, y —> oo) are

Ox = ‘0),).| = P0 (3)

for a remote shear loading, and vanishing tractions at infinity for an eigenstrain case,

2. the traction-free condition on the plane surface (2 = -1)

 

0.4:“, = 0
oz, ”-1 = O (4)
ozelzhl = O

3. and either perfect bonding boundary conditions at the inclusion-matrix interface

(r=a)

”Rh“. = ”Rik“.

“elk“. " uelkea

“(leza = amine

"MIR=a = ”MIR” ‘5)
“Rel“. ”Reign

0' =6
R¢R=a R¢R=a

or frictionless sliding at the interface with no separation in the normal direction (pure slid-

ing)

ORRIR=0 = ORR|R=a
GRGIRM = O (6)
6R9|R=a = O

6,“le = a _

CRwlR = a —

In the above expressions, at. and GU represent displacements and stresses, respectively.

2. METHOD OF SOLUTION

2.1 ELASTICITY THEORY IN THREE-DIMENSIONS

The problem considered in this thesis is solved by using the displacement poten-

tials approach. If no body forces are accounted for, the displacement equations of equilib-

rium for elasticity in three-dimensions are

V24+112vXVe = O

 

(7)

where 15 = displacement vector
v = Poisson’s ratio

e = dilatation, V01!

2 2 2
V2: 8 + a + 82,Laplacian operator

8x2 ayz dz

3 e a . . .. . .
V: —i +—' +—k , where i , ' , k are unit vectors in x, v. z direction res ec-
ax ay/ az J ' p
tively.

The general solution of the boundary value problem of the partial differential equations of
equilibrium can be given by harmonic displacement potentials.

In order to construct a solution to our problem we use a combination of six harmonic
displacement potential functions, namely (1)0, 11)], $2, (113, $4 and 13. Among them,
the set of 430, 4)], $2, $3 is due to Papkovich and Neuber, the set of $0, $3 and k3 to
Boussinesq, and the set of (110, (b4 and 2.3 to Dougall. For reference, we include the
eXpressions for the stress and displacement fields resulting from these potentials in the
Appendix A. It should be pointed out that the combination of the potentials chosen here

is not unique. Other combinations, consisting of at least three potentials, may also give

the same solution. We choose our combination for a mathematical convenience.

2.2 UNDISTURBED FIELD

According to a superposition principle in elasticity, for the applied remote loading
case the stress and displacement field in the matrix can be considered as a sum of two
parts, the undisturbed field in the absence of the inclusion caused only by the applied load-
ing at infinity and the field due to the disturbance by the inclusion. Similarly, for the
eigenstrain case, the stress and displacement field in the inclusion can be considered as a
sum of two parts, the nonelastic field when the inclusion is allowed to deform freely due to
the eigenstrain without any constrain from the matrix and a field resulting from elastic
strains caused by the presence of matrix.

For the matrix, the potential function which gives the elastic field due to the remote

shear loading (3) in the absence of inclusion is

l 2 7 1 2 2
4’0 = i170 (x ‘11") = pO-OR P2(i~1)00529 (8)

where u = coscp and P: (11) is the associated Legendre’s function of the first kind of

order n and degree m. The displacements and stresses corresponding to (8) are
ZGuR = ngRPg (u) c0520 = pORSIIlszCOSZB

2Gue = pORP§(u) c0820 = —pORsintpcos20

-3sintp

1

2Gu¢ - -6p0RP§’(u) simpcosZO = pORsintpcos¢c0520

oRR = £1061”; (p) c0520 = Posin2tpeos29

669 = —p0cos 20

Q
I

- pocosztpcos 20 (9)

(Ike = —3—S-i-;]Tpp0P§ (u) sin20 = —posintpsin20

Ge —p0cosq)sin 20

‘P

1p()P§’ (u) sincpcosZB = posintpcoscpc0520

where the prime refers to the derivative with respect to u and G is the shear modulus of
the matrix.
For the eigenstrain loading, the undisturbed stresses in the inclusion are zero and the

displacements. derived from displacement-strain relations, are as follows

a; = %€*RP:(u) c0529
-* l .1. 2 .
“e = —3£—18 RP2(u) sm29 _ (10)

_ 1 . 2,
u = male RP2 (u) c0520

wherefl = sinq)

2.3 HARMONIC POTENTIALS

The main concern in this investigation is to choose the proper harmonic potential
functions so that the specified boundary conditions can be satisfied by the corresponding
stresses and displacements. Generally, the harmonic potential functions which satisfy the
displacement equation of equilibrium in the forms of spherical and cylindrical coordinates
are related to Legendre functions and Bessel functions.

For the matrix, the potentials accounting for the disturbance due to the presence of

inclusion are

2 CnR_("+ ”P: (u) c0520

n=2

4’6

0° — n+1 1
$1 = 2 DnR ( )Pn(u)c050

n=l

¢2 = _ Z DnR-("+1)P;(u)sin0 (11)

n=l

0° -(n l) 2
zEnR + Pn(u)c0520

n=2

9-
L»
11

$4 = - Z n—_"1R’("+”Pn(u) c0520
n=2

While the following potentials allow to satisfy the traction-free condition (4)

¢0 = Jw3(l)12(}tr)e’kC0529dl
0

(b, = jw4(x)11(xr)e-keosed7t
O

10

115 = -jw4(x)11(xr)e-izsinedx (12)
0

4’3 = pulls (3.)12 (lr) e‘kcos20dk
0

2.3 = [(16 (2.) J2 (xr) e—ksinzedx
0
In the above equations C", D" and 15,, are the unknown constants, W3, W4, W5 and “’6
are the unknown functions which will be determined from the boundary conditions, and
Jn( 1r) is the Bessel function of the first kind of n order.

For the inclusion region we choose the following displacement potentials

Z C,R"Pf; (p) c0520

(D0 =
n=2
¢1 = Z DanPifll) cosB
n=l
62 = — 2 DanPi ()J.) sine (13)
n=l
63 = X EnR"Pfi(tt) c0520
n=2

The displacements and stresses can be derived from potentials (11)-(13) by using
relations given in the Appendix A. We combine these fields with the fields given in (9)
and/or (10) and find the total fields.

Note that for the remote shear loading the potential function (8) yields the stresses
0” = 43),), = p0 at infinity, while the other stress components are zero. The stresses

derived from the potentials (11) and (12) vanish at infinity. Therefore, the total stresses

11

satisfy the boundary conditions at infinity (3). For the eigenstrain case for the matrix
region we only use the potentials (11) and (12). These satisfy the traction-free boundary

conditions at infinity automatically.

2.4 SOLUTION

Potentials (11)-(12) are expressed in spherical and cylindrical coordinates, respec-
tively. In order to satisfy the traction-free boundary condition (4) on the surface 2 = -1 it is
convenient to use cylindrical coordinates. With the aid of the relation,

P__’,,"__ (11)

n+1

=(( -—1———n_)n - —-;)—!Jl”]m (kr)e7‘zd}t (z<0) (14)
R

we can express potentials (11) in terms of cylindrical coordinates as

(to = 2 CnIX"J2(kr)e7‘Zc0520d)t
n = 2 0
(1)] = z DJN’JI (Kr) ekzcosedl
n=1 0
(1)2 = “'2 D41”, (Xr) elzsinedk (15)
n - 2 0

$3 = Z EnJ'l"J2(Xr)e7‘Zcos29d2.
n = 2 0

 

(1)4 = .. 2 "hillkflz (1r) eAzCOSZBdA
n=2

where

 

~ _ (—1)"
C"‘Cn(n—2)!

(—1)"
'D"(n—1)!
(16)

 

aw)
I

(4)"
" E"(n-2)!

D1)

 

Then we substitute the potentials (00, (1)1, (1)2, (1)3, 64 and X3into the boundary

conditions (4) at the plane surface (2 = -1) as follows.

 

  

2
a a 8 a a 8
zzr-———O+[ 2% ¢ :lcose+2v1—gl s1n0+ [r— fizz-Way] sine

(“29.--1 = ‘57.: ”a7 36 322
2
a a a a a
+2v132cose+z ¢3—2(l—v)—¢3—r-—¢i—2(2—v);¢4=0
dz Brdz dz

r80 8:3

2 2

a a 8 a a

1 ¢O+ (b cos0+(l --2v) ¢ sin0+— (b2 sin0- (1 -2V)a—(:2 9059
dz 3082

 

(639):=-1 = r8082 +2308
512636 82¢ 36 32;.
z 3 1 3 l 4 1 4 3 _
+;aeaz ( ‘V2 )rae 2‘ar'a'6' “2“ “”736 Traz ' O “7)
2
8% 2‘12) 31¢ $2 3422 ,
(017)z=—1 - m+[r———-rarZa —-—(1—2v)a—z ]cos€)+[r:——:az —(1—2v)a—Z s1n6
2 2
22¢. 21. 12¢. M. ai. 1223
'“Mla +7363?-

+2r832(1_2V)§+2_r8—-02 az2

Note that the potential function (8) gives no tractions on the plane surface 2 — -1

Next we substitute the potentials (12) and (15) into the above equations (17) and express

boundary conditions on the traction-free plane (4) as

13

JR 2 {(A,‘,,2."+2(1—2v)la),,7t"—l—(}t+2-2v)lA§,,7t."_l }e-22k
0 n=2

+ {1113(2) +21!de + (—}t+2—2v)w5(}t) —2vf)1e-2x} )JZOtr)

+ {—b,xe-22+w,(x)}r1,(xr) ]e2x2dx = o

I“ Z {c,,x"+2(1—2v)D,,>.""—(x+1—2v)1‘5,,x"" }e'2"
0 n=2

—{\y3(}t) +2va40t) + (—}t+1—2v)w5(}t) +2lee_2}‘} ) 2JZ’Otr) (18)

2). °° .
+ 7[ (1 — 2v) 2 D,,?.""‘ .‘J-wé (7») ]J2 (Ar)

n=2
+(-g+r12){blle”n+w4(l)}Jz(7i.r) :lei‘d}. = 0

N

j[[— z {C,,A"+2(1—2v)D,,x""——(x+1—2v)i:,,x"" }
n=2

0

x -2 2
+ “113(k) 43%de + (-}t+1—2v)w5(}t) +2vD1e "1} )-%Jz(}tr)

+[—(1—2v) z ann—le-21+2vi)1e'2l+ 27VW40‘) +W6(A)]
n=2

szlz'Otr) ]d7\. = 0

Then, we let the coefficients of terms involving Jl (Ar), 12 (Ar) and J2' (Xr) be

zero and use the following relations:

14

(a) the Fourier-Bessel integral

f(x) = jjyrf<t)1,(yt)1,(yx)dydt (s>-1.x>0) (19)
00

(b) the transformation from the Legendre function to Bessel function

[mm (Ar) .22sz (z < 0) (20)
0

P301) __ (4).-..
Run " (n—m)!

 

(c) the relation between Bessel and modified Bessel function

co

5+1 b! S-IK ’b
j x -,+1J5(xb)dr= y! NO) (21)
Dixhyz) 2r(t+1)

 

 

where K 5 (z) is a modified Bessel function of the second kind and F(n) is a Gamma

function,

((1) the definition of the modified Bessel function
K1 /2 (A) = 191/2 (7.) = Jig-‘6" (22)
(e) the recursion formula for the Gamma function
F(x+ 1) = xF(x) (x>0) (23)

Then, the unknown functions W3: W4, W5 and 1116, which satisfy the boundary con-

dition at the plane surface (4), become

1113(2) =2(-x+1-2v)zb,e‘”‘+ Z [(—2x+3—4v)é,,x
n=2

15

°+2(1-2v)(—2x+3—4v)bn—2{2(1_v)(1-26)-)2mnW-1 —22

2. = J) x “2*
“I“ ) ‘ e (24)

—2(-2t+ 1— 2v)D1e‘2"

v5 (2»)

00

+ 2 [—2cnx—4(1—2v)b,, + (2x+3—4v)i~:,, 1x

n—l ~21
8

1116(2) = (1—22) EDNKHK”

n = 2
Next, we use the relation (16) and express the unknown functions
W3, W4, W52 and “’6 in terms of unknown constants C", D" and E", which are now the
only unknowns in the harmonic potential functions (1 l) and (12) and will be determined
from the boundary conditions at the inclusion—matrix interface (5) or (6).
It is convenient to use the spherical coordinate system to satisfy the boundary condi-
tions at the inclusion’s interface. We use the following relation to transform potentials

from Bessel function to Legendre’s function.

00

1,, (me 2...}: (—1>"”"—("—f-)%P,,(m (25)

Then, we rewrite potentials (12) in the spherical coordinates as follows

2 éanP: (u) c0520
n = 2

9-
o
11

Z nanPl (11) c050

n=1

6-
ll

16

62 = — Z nanP; (11) sin0 (26)
n =1
(1)3 = 2 CanPiut) cos20

n=2

2., = 2 KanP:(}.l) sin20

n=2

where

 

ee _}‘."
§n=IW3(7‘-) (+) d?)
0

(n 2)!

n

 

 

_ .. (4»)
"n ‘ "JV/40‘) (n+1)!d}"
,3 I (27)
— (_A)n+
Q" ‘ ’i‘VSO‘) (n+2)!‘0‘
0

n

_ .. <4»)
Kn ‘ “iwbm (n+2)!d}‘
0

 

After substituting 1413, W4, 1415, and W6’ given in equations (24), into (27) and using the

formula
J-e-bedx -_- F—(ij’l—ll (c>0, b>—l) (28)
c
0
we have

2
_ 2(1-2v) (1—2v) ]
E-‘n-—2I:Y0,n+2+ n+2 Y0,n+1+(n+]) (n+2)70,n D1

 

+ z [2(n+3)Ym—2,n+3+ (3_4V)Ym—2,n+2] Cm
m=2

17

 

3— 4v
+ 22(1_2v)[n2+27m—1,n+l (n+2) (n+1)Ym—1,mn]D

m=2

 

2(1—V) (l—2V)
+ m§22[ n+2 Ym—2,n+l—(n+3)Ym—2,n+3:lEm

 

 

T"n = Dly0,n+l
(29)
1-2v
5n = '2[70,n+2+_n:'2” 70.n+1]D1
°° 4 1—2v
+ Z 2(n+3)Ym-— 2,mn+3C - 2 ( )Ym—l.n+le
n+2 1
m=2 m=2
—2 [2(n+3)Ym—2,n+3_(3_4V)Ym-2,n+2]Em
m=2
l—2v
’ —;2(n+2) (n+1)Ym-1.an
where
_ (—1)”*"(p+q)!
1M — (30)

p!q! 2p+q+l

Then, we express displacements and stresses in terms of the potentials (11), (13) and
(26) by using the relations given in the Appendix A. Finally we substitute the displace-
ments and stresses, given by the potentials (11)-(13), and the undisturbed ones (9) in the
matrix and/or non-elastic ones (10) in the inclusion, into equations (5) and use the recur-
sion formula for Legendre functions given in the Appendix B. Then, the boundary condi-
tions at the interface, for the case of perfect bonding, become:

u =17
RR=a RR=a

l8

5—4v 2 0° C
——3 Dlpz(p)+2[{xmcn+KnD D+I+K§n_lE_l+KE E

2 ,n+1 n n n,n+1 n+1
0 n=2

K
+K§,n§n+K2,n—lnn-l +KI'21,n+lnn+l +Kr§,n—1§n—l +Kr§,n+l§n+l+Kn,nKn,n}

1 C '
-I:{Kn,ncn+Krli),n—1Dn-l+Kr€n+an+l+Kt€n-1En—1
K5 E P2 = 1 —2G * P2 31
+ n,n+1 run-+1} "(i-l) 3 po 8 a 2“” ( )
“9! = 179'

n n,n+l n+1

Z[{Hinc"+H'l')‘"‘1D"'1+HRDJI+1DN+1+HIEII—1E —l+HE E +
n=2

+H§ §n+Hrun—lnn—l +Hn

m" mn+lnn+

1+H§,n—1§n—1+H§ 1Cn+l+HrinKn

71,714"

] C
K _ D D E
+Hn,n+2Kn+2}_I“{Hn,ntn+Hn,n—1Dn-1+Hn,n+1Dn+1+Hn.n—1E

n-l

1 . 1
+H£n+1En+1}]P3(ll)=6(P0-ZG€)0P§(M) (32)

u¢|R=a = i-lq)|R=a

E
n-1+Rn,n+lEn+1

co C
2 [{Rn,nCn+RnD,n-1Dn—l +RnD,n+an+l+Rr11§,n-1E
n=2

+RE,n§n+RI?.n-lnn—l +1?"

mn+1nn+l

+R§,n—1Cn-l "I'RC,n-+l§n+l+RK K

n n,n n

1 C D D F
+Rrfin+2Kn+2 }—f~{Rn,nC +Rn,n-1Dn—1+Rn,n+1Dn+l +Ruin-1E

n n-l

19

+Rr€n+1En+l } ]P5’(p)

11.1: n n,n+1

+ 2 [{RanDn+RE E +R2nnn+R§an+RK 1c"+1 }
n=2 .

1 . 3 I at I
_1_.{R,I,2,,D,,+R,{{,,E,,} ]P,;(11) = goo—268 >an (11) (33>
GRR|R=0 = ORRlea

E
n-1+Ln.n+lEn+l

2 S-v 2 0° C
—(—3-T-)-D1P2(tt) + Z [{L C +LD Dn+l+Ln€n_lE

n,n n n,n+l
0 n=2

+L§ €11+Lii|,n-lnn-l+Lr'11,n+lnn+l +Lr§,n—l§n-l+Lr§,n+lcn+l+LK K }

II, n n, in II, n

_{L:nCn+LnD,n—1Dn-l+Lrli),n+an+l+LnE,n—1En—1+L£n+lEn+l } :] P3(u)

= ——P§(u) (34)

+5E E

n,n n n n,n+l n+1 n,n-1 n—l

5—4v 2 C
—3—D,P2(tt) + z [{s C +sgn_,D _1+sD D

C
K K D D E
+5 K +5 n+2Kn+2}—{SN.HCH+SII.n—1Dn-l+Sn,n+1Dn+1+Sn.n—1E

n,n n n, n—l

P .
+S£n+lEn+l } ]P,§(u) = goPgm) (35)

GR‘pIR = a = ORWIR = a

5—4v 2 1 " c
3 D1P2(u) +3 2 [{Tn.nCn+Trlt),n—1Dn—1+TnD,n+an+l+Tr€n-1E
n=2

n—l
60

+71; E +1I".TingnflkTr'11,r1—1r'rz-—l-+-Tr‘i‘,n+lnn+l+TIC1",n-1Cn-1+TC

n,n+l n n.n+l§n+l

C
K K D D E
+Tn,nKn+Tn,n-i.2Kn-i.2}"’{Tn,nz:n+Tn,r1—an-1'I'Tn,n+1Dn-i-l+Tn,n—lEn—l

+735n+15n+t}]P3’<u>+%2[17f.0.+T£.E.+7:1,.n.+T§,.c.+
n=2
p , .
Trin+lxn+l }—{7f.nDn+TIlzi:nEn } ]P5(ll) = 3013301) (36)

where the primed quantities are the first derivatives with respect to u and 1" = G/ G . The

coefficients are defined as

 

 

C _ n+1
KIM — —an+2
_ 2(1—2v)
K’lgn‘tl - nan+2
KE __ (n+3—4v)(n—2)
“-1 (2n—l)a”
KE ("+5-4V) (n+3)

 

n.n+1- (2n+3)an+2

21

_n—4+4v
2n—l

n-l

_n—2+4v
2n+3

an+l

 

= (n—4+4v) (n-2)a,,_l
2n—l

= (n-2+4v) (n+3)a

 

 

n+1
2n+3
4a""
1 _ C
—an+2_Rn»"
1—2v _ D

 

= (2n-1)a" - "’"‘]

 

_ 1—2v _ D
-—(2n+3)a”+2_ "’"H
"—2 _ E

 

R
(2n—1)a" "'"71

n+3 E

- R
(2n+3)a"+2 "’"H

 

n—l _
—a —Rs~n

(37)

(38)

 

 

 

mn—l

 

22

mn+2

1-2v _
”Rn—1:57:10" '=R"
1-2\r ,
”2n“ =-§7,‘:§a”1 =13"
_ (n-2) - _
Hg’n_1——2n_1a" 1-
n+3
H§,n+1 =_2n+3an+1=
H" = Mun—l — R"
M 2n—l ‘ .
(n+3) (n+4)

K = n+1 = K
”M” 2n+3 a R
RD _ (n+1)(l—2v)

mn — — n+1

(n—l)a
_4(1-v)
R£fl+l — an+l
R2,n—1=’2(1‘V)an
REM] =4(1-V)a"
REM] =n(n+3)a"

LC _(n+U(n+D

mn — an+3
LD 2(1-2V)(n+2)

mn+l

nan+3

(39)

23

__ (n—2) (n(n+3) —2v)

 

 

E

Lmn—l — (2n__1)an+l

LE _ (n+2)(n+3) (n+5—4V)
run+l - (2n_+3)an+3

La = n(n—l)a"‘2

(n— 1) (n—4+4v)an_2

 

 

 

 

"’"‘l 211—]
L“ __(n+1)(n—2)—2v ,,
an+l - 2n-F3

_(n—l)(n—2) (n—4+4v) _
LEM'“ 2n—l an 2
C = (n+3)((n—2)(n+1)—-2V) n
L"t"+‘ 2n+3 a

LK = 4(n—l)a"’2

n. n

SC = 2(n+2) ___ TC

mn mn
an+3

 

 

D _ (1—2v) (n+1) _ D

Sn.n—l — (Zn-l)an+] - Tn,n-1

D _ (1-2V) (n—3)(n+2) D

Sn n+1 ‘ = T +1
’ n(2n+3)a"+3 "’"

SE _ 2(72-2) ("+2-2V) = TE

 

n,n—l - (2n_])an+] ILR-I

(40)

24

_ 2(n+3) (n+4—-2V)

 

 

 

 

 

 

 

 

5E _ _
n,n+l (2n+3)an+3 - Tini-I
t
sin = —2(n—1)a"‘2 = TS."
2(1—vn)
$71 =____ —2-
n,n-1 211-] an ’77:;1—1 (41)
2(l—V(n+2))
S“ =— —
n,n+1 2714-3 an—Tg,n+1
2 -2 -
Sg‘n-l _ (n )2(n 3+2V) (In-2 = TC
n_1 n,n-—l
2(n+3)(n—l+2
5C =_ V) _
n,n+l 2714-3 an — Tr§,n+1
Sf,"=—(n(n+2)+4-(n_2)2(n+2))a"‘2=T"
211-] "v"
S" =n(n+3)(n+4) _
n,n+2 2n+3 an — Thin-+2
TD ____(1—2v)(n+1)(n+2)
"’" 2(n—l)a"+2
TE _ 2(1—v)(n+2)
"in an+2
721 =— - — -'
,,,,, (l V) (n l)a" (42)

T§,n=2(1-v)(n—1)a"'1

-l
Trin+1=n(n )2(n+3)a"’l

 

25

All terms with a bar in equations (31)-(36) can be obtained from the corresponding terms
without a bar, given in eqns. (37)-(42), by replacing i, 11, C and v with C, D, E and V.
For example, KC K, n_ 1, KHDH 1, Kfn_1 and K5 M1 can be obtained from K5,”

K"

n, n’
n,n_1,K‘1n,n+1 KCH’MI and Kg,“ + 1, respectively, by replacing v with V.

Using eqns. (29), we reduce the unknown constants in equations (31)-(36), repre-
senting the boundary conditions on the interface, to only C", D", E", C", D" and En.
Finally, we equate the coefficients of P3 (1.1) and PS’ ( u) on both sides of eqns. (31)-(36)
for each n from n=2 to n -—> oo , and obtain an infinite set of algebraic equations. Each of
these sets contains six equations and six unknowns. For the numerical solution we trun-
cate this infinite set of equations at n = N. Therefore, there are 6N equations and 6N
unknowns to be solved. We truncate the series such that the boundary conditions (5) or
(6) are satisfied to at least three significant figures. After these constants are evaluated,
the stresses and displacements are known everywhere in the matrix and in the inclusion.

We follow the same procedure for the sliding interface case, but use the sliding
boundary conditions (6) instead of (5). The sliding boundary conditions at the interface

include eqns. (31) and (34)—(36) and

ORGIRza = O

2 [Sincn+SEn-1Dn-l +SnD,n+an+l+SIEn—1En—I+Sr€n+lEn+1 11’3“” =0

(43)

6R¢|R=a = 0

Dn+l+TIEn—1En—l +TSI€n+lEn+11P3’(u)

n,n+l

é22[Tn,nD,nnCn+T-1Dn-1+TD

(44)

n,nn

+£X [if D +75 fling-(it)

26

In the calculations, for the case of a perfectly bonded inclusion, we set the constant
D1 to zero since D1 and C2 play the same role in the potentials and only one needs to be
kept. Furthermore, we take one of the constants D2 and E 2 as zero for the case of the per-
fectly bonded inclusion, and set both constants to zero for the sliding inclusion case.

In particular, when the matrix is infinite, the analytical solution involves only a finite
number of terms in potential functions. For the case of an infinite body containing a spher-
ical cavity and subjected to uniform shear stresses (1'm = ”on = p0 at infinity, the con-
stants become

-_‘__ 5
C2‘2(7—5v)“

5 3

2 - “W0 (45)

while the other terms vanish.

For the case of an infinite body containing a spherical inhomogeneity subjected to

 

uniform shear stresses on = ‘ny = p0 at infinity with the pure sliding interface, the
solution is
E _ 701‘ (1 —v) _1_
3' (l7—l9v)(7+59)+4(7—5v)(7-4V)az
C2 = —(l + 2%)(1253
(46)
C4 _ ’4—vE3

27

 

3
__3a é 2(7+SV) 2 )
D" 7—5v(6+ 7 “E3
a2 3
C2=§Z( -2(1+V)DI)

with the other terms vanishing. This solution in a form of finite series is expected from

the work of Ghahremani (1980). The solution for perfectly bonded spherical inclusion in

an infinite space is also expressed in terms of finite series as shown by Goodier (1933).
Similarly, the solution for the spherical inclusion in an infinite space with an eigen-

strain loading case involves finite series for both perfect bonding and sliding cases.

3. RESULTS AND DISCUSSION

We carry out the numerical computations for various radii of inclusion a ranging
from 0.2 to 0.99 and for the different ratios of shear modulus of the inclusion to the matrix
I‘ = G/ G. We take the Poisson’s ratio as 0.3 for both the inclusion and the matrix, for
simplicity. We give the numerical results for the case of the uniform shear loading in
Figs. 2-19 and for the uniform shear eigenstrain in Fig. 20.

Fig. 2 shows stresses 6W and CW along the surface of the perfectly bonded inclu-
sion when the radius of inclusion a = 0.2, 0.5 and 0.8 (plotted in solid, dashed and dashdot
lines respectively), 8 = 00, and F = 100. Recall that the inclusion is located a unit dis-
tance from the free surface so the larger the radius a the closer is the inclusion to the free
surface. As expected, the radius of the inclusion has a small influence on the stresses CW
and CW when the angle (p is less than 900, and the larger effect when (p is greater than
900 which corresponds to the region near the free surface. Note that the stress 6W in the
inclusion is more influenced by the free surface than the corresponding stress component
0W in the matrix.

Figs. 3 - 5 show stresses 6W and aw in the same cases as in Fig. 2, but the ratio of
shear moduli is 2, 0.5, 0.0] respectively. Once again, the radius of the inclusion has a
small influence on the stresses CW and 6‘94) when the angle tp is less than 900, and the
larger effect when (p is greater than 900. Further more, the following can be observed as
the ratio of shear moduli decreases. The influence of the free surface on stress aw in the

matrix becomes larger and larger and on stress 6w smaller and smaller. The stress Ow

almost vanishes as I“ = 001 which is close to the case of a cavity so that the obtained

28

29

result is reasonable. The stress 6W in the matrix increases around the points (p = 00
and 1800 while decreases around point (p = 900, and it will be compressive when the
ratio F is less then 1, or say, the inclusion is weaker than the matrix.

Figs. 6 — 9 show stresses CW and CW along the surface of the sliding inclusion
when the radius of inclusion a = 0.2, 0.5 and 0.8 (plotted in solid, dashed and dashdot lines
respectively), 6 = 00. The corresponding ratios of shear moduli F are 100, 2, 0.5, 0.01.
As expected, similar to the case of perfectly bonded inclusion, the radius of the inclusion
has a small influence on the stresses 6W and CW when the angle (p is less than 900, and
the larger effect when (p is greater than 900 which corresponds to the region near the free
surface. Note that the stress 6W in the inclusion is more influenced by the free surface

than the corresponding stress component 0' in the matrix when F takes big values,

‘9‘?
while the stress 6W in the inclusion is less influenced by the free surface than the corre-

sponding stress component 6 in the matrix when F takes small values. As I“ = 0,01

‘9‘?
which is close to the case of a cavity, the stress 6W in the inclusion goes towards zero.
The stress aw in the matrix increases around the points (p = 00 and 1800 as the ratios
of shear modulus F decrease, however, it decreases around point (p = 900, and is always
compressive there regardless of value of ratio, F .

Fig. 10 compares the stress 6W versus the angle (p for the cases of perfectly
bonded and sliding interfaces when F = 2. We observe that the stress increases for the
sliding case and decreases for the perfect case as (p increases from 00 to 900, and vice
versa from 900 to 1800, so that 6W reaches its maximum value at (p = 00 and the mini—
mum value at (p = 900 for the case of perfectly bonded interface, while the minimum at
(p = 00 and the maximum at (p = 900 for the case of sliding interface. The maximum
value of CW in the inclusion for the sliding case is about 2.7 pO for F = 2 and increases
as F goes up. The similar situation occurs for OW for the case of perfectly bonded inter-
face when F increases. For example, when F = 100 the maximum of stress a in the

(PIP

matrix occurs at q) = 900 but when F = 2 at (p = 00. In this case the stress OW is much

30

lower than that of the sliding case and has a value of 0.8 to 0.9pO depending on the radius
of the inclusion (i.e. the distance from the free surface). Note that the traction-free sur-
face has a higher influence on OW for the sliding case. As the ratio F becomes small, say
0.01, which is close to the case of a cavity, the stresses ow and CW coincide in cases of
both perfectly bonded and sliding inclusion.

Figs. 11 - l4 illustrate the jump of the tangential displacement [14¢] along the inter-
face for various radii of the inclusion a = 0.2, 0.5 and 0.8 (plotted in solid, dashed and
dashdot lines respectively) and 6 = 00 for the sliding case. The same as before, the corre-
sponding ratios of shear modulus F are 100, 2, 0.5, 0.01. Note that the jump of the dis-
placement is higher for (p 2 9O0 (i.e. near the traction-free surface) and increases when the
inclusion is closer to the free surface. This is true for all ratios of shear modulus F,
according to the results from our sample computations. This is expected since the inclu-
sion can deform more freely near the free surface.

Figs. 15 - 18 illustrates the stresses on and on along the z-axis from 2 = -l to z = 1
when the inclusion is perfectly bonded. The radius of the inclusion takes on the values
a = 0.2, 0.5 and 0.8 (plotted in solid, dashed and dashdot lines respectively). The corre-
sponding ratios of shear modulus F are 100, 2, 0.5, 0.01. First, we consider the case in
which the ratio F is larger than 1, i.e., the case of stiffer inclusion. The stress on in the
inclusion decreases from 2 = a to z = -a when the inclusion is close to the free surface
( a 2 0.5 ). The slope of 6“ goes to zero when the radius is small, which implies that for
the situation of 0 << 1, equivalent to that of an infinite body containing the inclusion, the
stress on is uniform which is expected from Eshelby’s (1957) solution. The curves have
similar forms for other F > 1, but have bottom up image for F < 1. Note from Fig. 18 that
when F is very small (close to the case of a cavity), there is large stress concentration at
the interface near the free surface. When F = 1 the curve is a straight line, as expected,
since this is the case of a homogeneous material.

Fig. 19 shows the stresses on and 6er at selected points P, M1, 11,12, and M2 along

31

the z-axis (see Fig.1) when the radius a of the inclusion varies continuously from 0.2 to
0.99 for both the perfectly bonded (dashed curves) and sliding (solid curves) cases and
F = 100. The point P is on the free plane surface, M1 and M2 on the side of the matrix at
the interface, 11 and 12 on the side of the inclusion at the interface. It can be seen that the
free plane surface will contribute significantly to the stresses on and on at the points P,
M1 and 11, but will have a very small effect at points 12 and M2, which are away from the
plane boundary.

For the case of the perfectly bonded interface, Oxx at M1 and M2 in the matrix is
very small (almost zero for F = 100), but on. at 11 and 12 in the inclusion is high (larger
than 2pO for F = 100). This is expected because when the inclusion is much stiffer than
the matrix it carries most of the loading. The stress oM at the point P decreases when the
radius a increases. The inverse situation occurs when the inclusion is softer than the
matrix.

For the case of the sliding inclusion, when the inclusion is not close to the free plane
surface, the stresses (3m and on at these five points are not greater than the applied stress.
But if the radius of the inclusion is close to unity, there is a high stress concentration in
on at the points P and M1 and it increases as the radius a increases. Also, the stresses are
very high in the extreme situation of a cavity addressed by Tsuchida and Nakahara (1974).

Fig. 20 illustrates the variation of stress on and on along 2 axis from 2 = -1 to z =
1 when the inclusion is perfectly bonded at the interface for the shear eigenstrain loading
case. The radius of the inclusion has the value of 0.2, 0.5 and 0.8, and F = 2. It is inter-
esting to note that comparing the results to a similar case for the shear stress loading prob-
lem with the same F, Fig. 7 has an inverse image. However it is similar to one for the
case of a uniform shear loading when F = 0.5.

In conclusion, in this paper we solve the problem of an inclusion in a half-space
under an asymmetric loading, either remote or eigenstrain loading of shear type. We study

the joint effect of traction-free surface and the boundary conditions at the matrix-inclusion

interface, which is either perfectly bonded or allows slip. We show that both the free sur-
face and sliding may alter significantly the stress fields and may cause high stress concen-

trations.

33

 

 

 

 

 

 

 

 

P0

Fig. l Spherical inhomogeneity in a half-space.

34

 

6W (6w) /p0
2,5 I I r T r r l 1
2p“'.\: ’d-‘l
v“ 'I
\ I
A I,
I .....
0"
.1
13r- ./ d
.1
r
I.
(I
C
1. 99 .
s 5.
~.
\. ‘
\. ‘.c\
5 'Qt
0'5,- .x "t‘. "
\ \
\ I .\I\\
\x I ‘. ‘\\
a \
\ O / ' ‘
“W ‘P‘p / ‘
~\ 1 '''''
ov- '.‘s-_.r.,
J l

 

 

 

4,5 l J l l 1 ' l
0 20 40 60 80 100 120 140 160 180

Fig.2 Stress 0' 6 versus angle (9 for different radii a = 0.2 (solid line), 0.5

¢¢( ww)
(dashed line) and 0.8 (dashdot line) when r = 100 and e = 0° for perfect

bonding and remote shear loading case.

 

 

35

6‘1"? (6WD) /p0

1.4 I I T I I I I I

 

O
........

 

 

)—
l-
r.

l

20 40 60 80 100 120 140 160 180

 

-0.2
0

Fig. 3 Stress 6W“? versus angle (p for different radii a = 0.2 (solid line), 0.5

‘9‘?)
(dashed line) and 0.8 (dashdot line) when F = 2 and 9 = O0 for perfect

bonding and remote shear loading case.

36

CW (6w) / p o
1.5 . . . . . . . .

 

I
‘I
I

1.4

I
g
l

1.2

¢¢

'
l

0.8

I
n

0.4

I
l

0.2

§>
60
I

 

 

 

 

-0.4 I l l l l l
0 2O 4O 60 80 100 120 140 160 180

Fig. 4 Stress 0 versus angle (p for different radii a = 0.2, (solid line), 0.5

99(699)
(dashed line) and 0.8 (dashdot line) when F = 0.5 and 9 = O0 for perfect

bonding and remote shear loading case.

37

 

 

 

 

 

 

Fig.5 Stress 0'

(C

W W) versus angle (p for different radii a = 0.2 (solid line), 0.5

(dashed line) and 0.8 (dashdot line) when F = 0.01 and 9 = O0 for perfect

bonding and remote shear loading case.

38

 

 

 

 

 

'o 20 40 so to no 120 140 150 no

(P

Fig. 6 Stress OWN versus angle (p for different radii a = 0.2 (solid line), 0.5

(PCP)
(dashed line) and 0.8 (dashdot line) when F = 100 and 6 = 00 for sliding and

remote shear loading case.

 

39

 

 

 

 

 

'o 2o 4o so so 100 no 14o 160 no

Fig. 7 Stress o (6 versus angle (p for different radii a = 0.2 (solid line), 0.5

$9 $9)
(dashed line) and 0.8 (dashdot line) when F = 2 and 9 = O0 for sliding and

remote shear loading case.

40

 

 

 

 

o o /
<p<p( (mp) p o
2 I l T I T T I I
0"!
.I
.’ J’.
I I
1.5- o ,v .
(PRO ’
1- '~ _ , -
\ , .’.
' I
0 \
o.s~ ; .
I ~_
I "v. .
r. 7
,-
o- ’ .
u(‘
-05- . 69¢ -4
1 l l J l l l l l

'0 2o 4o so so no no no no no

Fig. 8 Stress OW“? versus angle (p for different radii a = 0.2 (solid line), 0.5

WP)
(dashed line) and 0.8 (dashdot line) when F = 0.5 and 6 = O0 for sliding and

remote shear loading case.

41

 

 

 

 

 

'o 2o so so so no no no no no

 

Fig. 9 Stress ow“) versus angle (p for different radii a = 0.2 (solid line), 0.5

(Pl?)
(dashed line) and 0.8 (dashdot line) when F = 0.01 and 6 = O0 for sliding

and remote shear loading case.

42

GW/ p o

 

2.5 *- .

pure sliding

1.3 _ . q

0.5 ~
perfect bonding

 

 

 

Fig.10 Stress CW versus angle (p for perfect bonding and sliding cases when
I“ = 0,5, 9 = O0 and different radii a = 0.2 (solid line), 0.5 (dashed line) and

0.8 (dashdot line) for remote shear loading

43

 

G[u(p]/p0
62 - . . I . . . r

on x. .

 

-0J5

I

o
. 4"
I

 

 

 

20 40 6O 80 100 120 140 160 180

Fig.1] Jump in displacement [uq’] versus angle (p for different radii a = 0.2 (solid
line), 0.5 (dashed line) and 0.8 (dashdot line) when F = 100 , 9 = O0 for pure

sliding case and remote shear loading.

44

 

0.25 r T . . . . . .

o.2- , ’ """"" ‘ . -
I \.

on- , -

 

 

 

 

-02 - ‘ \ .r -
43.25 1 l l l l l l l
O 20 40 60 80 100 120 140 160 1110

Fig.12 Jump in displacement [u‘p] versus angle (p for different radii a = 0.2 (solid
line), 0.5 (dashed line) and 0.8 (dashdot line) when F = 2 , 9 = O0 for pure

sliding case and remote shear loading.

45

 

ooooo

 

I
J
I

l

on '_ ‘ .’

------

 

 

 

-0. 15 I I I I I J I I
0 2O 40 60 80 100 120 140 160 180

Fig.13 Jump in displacement [144)] versus angle (p for different radii a = 0.2 (solid
line), 0.5 (dashed line) and 0.8 (dashdot line) when F = 0.5 , 9 = 00 for pure

sliding case and remote shear loading.

 

46

 

0.15 I r 1 r I I r I
d -
.l s
I \
I \
I \
’ O
01» I \ -
o . I
I s
I ’ ’’’’ -\\ \
’ ’ \ ‘
I
\
I ” \ s

 

 

 

 

'01” ’- .s \\ ‘ ll .1. "
\ ‘ I’ ’
‘ \\\ l’ I
u\ \~~~ -' I, ’
n\ - ’
-0.1 *- ‘. ' 4
\ I
K. .1.
$.15 l l I I l l I I
0 20 40 60 80 100 120 140 I60 180

(P

Fig.14 Jump in displacement [um] versus angle (p for different radii a = 0.2 (solid
line), 0.5 (dashed line) and 0.8 (dashdot line) when r = 0.01 , 9 = 0° for

pure sliding case and remote shear loading.

47

OXX (61.x) /p0

 

 

 

 

    

 

 

 

 

 

 

2.5 I I I I I I I I I
1- ----------------------- _ f:fJ--.JO--P"c -
2 1 --------- l-""' I 1
| -‘-.P'- I .
.N” I
P.-. I .
- ..... 1 : .
I.‘ | '
1 1 1 f
, 1 | I
I I I i
1'5.- I I I o '1
, l ' 1
l I .
g . . .
I l '
' l
| I .
l | l I
5 l i .
1 l 1 l 1T
I 1 - ,
a“ . I fl— L3,.
“~L I I ,1
I‘ I I ’ o
' \s I I If I
! ‘\ 1 1 (I '
' ~ 1 l ' i ,
0.5- I ‘t | ' I! i ’Ifl
i \ I I f , .
'\ . ll It I f
.\ I ‘I I) i ’
', ' ll 11 .-
. E I1 I! '3
‘0' ll I! "
\. I‘ u I.
0 'l 1 ' 1 1 I ‘ 1 i
-1 0 8 -0 6 4) 4 41.2 0 0.2 O 4 0.6 0.8 1

Fig.15 Stress Gum“) along z-axis for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when F = 100 for perfect bonding case and

remote shear loading.

48

0.. (6..) /p0

 

1.4 I I I I I I I I I

 

---------------------------
--

P ”.-
I-.

u--
.....
d.

 

 

       

0.9 I-

0.8 '-

 

’

0.7 - ’I

 

 

-l
‘t

 

 

 

0.6 l l l l l l l l l
-l -0.II -0.6 «0.4 -0.2 0 0.2 0.4 0.6 0.8 1

Fig.16 Stress cum”) along z-axis for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when P = 2 for perfect bonding case and

remote shear loading.

49

.011 (cxx) /p0

 

1.5 I I I I r I I I I

1.4 -

 

‘-

 

13 " I

12-

 

 

0.8 r

I

0.? .........

 

 

I
I
I
l
I
I
I
I
I
I
I
I
I
I
I
I
l
I
I
I
I
I
I
I
I
I
I
I
l .
I

 

 

 

 

0.6 l l l l l l l l l

Fig.17 Stress onwxx) along z-axis for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when P = 0.5 for perfect bonding case and

remote shear loading.

50

611(61X) /p0

 

       

 

 

 

 

 

 

 

 

2.5 "I I I I I I I I I
I
[I
I'!
i I
I' !
2" l I ..
r. I 'I I. '.
, I I: I 3'.
i. I I II I I.
p ' I \
E I . I -\
I I I I .‘
1.5 I' i ’ : : I .-
I " l l a I ‘
i x I I ‘. g
i” I I \.. I
I-’.’i I I ‘i‘q~~
4' I i I
1 i I I g
i I I !
I I I I
I I I !
i I : I
0.5 - i I . ! -
i ' I I
. I I .
I ' ' I
I I I !
I I I I
o I ------- l-°-I'-—!--—- ; — .L J -., ..... -1'
-1 4).! -0.6 0.4 -0.2 0 0.2 0.4 0.6 0.8 1
Z

Fig.18 Stress Gum”) along z-axis for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when 1" = 0.01 for perfect bonding case

and remote shear loading.

51

 

 

 

\__'_;——— aP
M2 "'

 

 

 

 

l2 FIJI/Ix
l1 /
0.5“ I/P 'l
I
I
I
0...;ar-M="-----------------------------:
M2
.0.5 l I l l l I l
1 1.5 2 2.5 3 3.5 4 4.5 5

1/a

Fig. 19 Stress on (On) at several points along z-axis versus radius of inclusion a for
both perfect bonding (dashed lines) and pure sliding (solid lines) when l" = 100

- (remote shear loading).

52

on; (ch) / (202*)

 

008 I I I

-

I
' x

0.4

I
\

0.2

\
—*_——_—_—_—

 

 

 

I

-0.2

05.
5".

-0.6"

_+—-_-—__-“

 

 

 

F.““I-o-

 

 

 

 

‘0.8 l l l
-1 . .

Fig.20 Stress Gum”) along z-axis for different radii a = 0.2 (solid line), 0.5
(dashed line) and 0.8 (dashdot line) when P = 100 for perfect bonding and the

eigenstrain case.

REFERENCES

Atsumi, A. and Itou, S., 1974, “Stresses in a transversely isotropic half space having a

spherical cavity,” ASME Journal of Applied Mechanics, lb]. 41, pp. 708-712.
Eshelby, J .D., 1957, “The determination of the elastic field of an ellipoidal inclusion, and
related problems,” Proceedings of Royal Society, London, Series A, Vol. 241, pp.

376-396.

Ghahremani, F., 1980, “Effect of grain boundary sliding on anelasticity of polycrystals,”
International Journal of Solids and Structures, Vol. 16, pp. 825-845.

Goodier, J.N., 1933, “Concentration of stress around spherical and cylindrical inclusions

and flaws,” ASME Transactions, Vol. 55, pp. A-39-44.

Jasiuk, 1., Tsuchida, E. and Mura, T., 1991, “Spheroidal sliding inclusion in an elastic half-

space,” Applied Mechanics Reviews, Vol. 44, No. 11, pp. 143-149.

Jasiuk, I., Tsuchida, E. and Mura, T., 1987, “The sliding inclusion under shear,” Intema-

tional Journal of Solids and Structures, Vol. 23, No. 10, pp. 1373-1385.

Jeffery, 6.8., 1920, “Plane stress and plane strain in bipolar coordinates,’ Transactions of

Royal Society, Series A, Vol. 221, pp. 265-293.

53

54

Lee, M., Jasiuk, I. and Tsuchida,E., 1992, “The sliding circular inclusion in an elastic half-

plane,” ASME Journal of Applied Mechanics,” Vol. 59, pp. 857-864.

Mindlin, RD, 1948, “Stress distribution around a hole near the edge of a plate under ten-

sion,” Proceedings of the Society for Experimental Stress Analysis, Vol. 5, pp. 56-

68.

Mindlin, RD. (and Cheng, DH, 1950, “Thermoelastic stress in the semi-infinite solid,

Journal of Applied Physics, Vol. 21, pp. 931-933.

Mura, T., 1987, Micromechanics of Defects in Solids, 2nd. ed. Martinus Nijhoff Publish-

ers, Dordrecht.

Richardson, M.K., 1969, “Interference stresses in a half plane containing an elastic disk of

the same material,” ASME Journal of Applied Mechanics, Vol. 36, pp. 128-130.

Saleme, E.M., 1958, “Stress distribution around a circular inclusion in a semi-infinite elas-

tic plate,” ASME Journal of Applied Mechanics, Vol. 25, pp. 129-135.
Seo, K. and Mura, T., 1979, “The elastic field in a half-space due to ellipsoidal inclusions
with uniform dilatational eigenstrains,” ASME Journal of Applied Mechanics, Vol.

46. pp. 568-572.

Sheng, P.Y., 1992, The Stress Field of a Spheroidal Sliding Inclusion under Shear, M.S.

Dissertation, Michigan State University.

Shioya, S., 1967, “On a semi-infinite thin plate with a circular inclusion under uniform

I,

55

tension,” Bulletin of the JSME, Vol. 10, pp. 1-9.

Tsuchida, E. and Mura, T., 1983, “The stress field in an elastic half space having a spheroi-
dal inhomogeneity under all-around tension parallel to the plane boundary,” ASME
Journal of Applied Mechanics, Vol. 50, pp. 807-816.

Tsuchida, E. and I. Nakahara, I. 1970, “Three-dimensional stress concentration around a
spherical cavity in a semi-infinite elastic body,” Bulletin of the JSME, Vol.13,
No.58, pp. 499-508.

Tsuchida, E. and Nakahara, 1., 1972, “Stresses in a semi-infinite body subjected to uniform
pressure on the surface of a cavity and the plane boundary,” Bulletin of the JSME,

Vol.15, No.79, pp. 1-10.

Tsuchida, E. and Nakahara, I., 1974, “Stress concentration around a spherical cavity in a
semi-infinite body under uniaxial tension,” Bulletin of the JSME, Vol. 17, No. 1 l 2,

pp.1207-1217.

Tsuchida, E, I. Nakahara, I., and M.Kodama, M., 1973, “Stresses in a semi-infinite body
having a spherical cavity compressed by rigid body, “ Proceedings of the 23rd
Japan National Congress for Applied Mechanics, Vol. 23, pp. 531-546.

Tsuchida, E., Saito, Y., Nakahara, I., and Kodama, M., 1982, “Stress concentration around
a prolate spheroidal cavity in a semi-infinite elastic body under all-around ten-

sion,” Bulletin of the JSME, Vol. 25, No. 202, pp. 493-500.

Tsutsui, S. and Saito, K., 1973, “On the effect of a free surface and a spherical inhomoge-

56

neity on stress fields of a semi-infinite medium under axisymmetric tension,” Pro-
ceedings of the 23rd Japan National Congress for Applied Mechanics, Vol. 23, pp.
547-560.

Wachtman, J .B. and Dundurs, J., 1971, “Large localized surface stresses caused by ther-

mal expansion anisotropy,” Journal of the American Ceramic Society, Vol. 54, pp.

525-526.

Yu, H.Y. and Sanday, SC, 1990, “Axisymmetric inclusion in a half space”, ASME Jour-
nal of Applied Mechanics, Vol. 57, pp. 74-77.

APPENDIX A

displacement and stress fields in spherical coordinates

for (be

 

 

a
ZGuR = 3%)
_ l 8%
20“" _ Rsintpfi
13%
ZCJ’l/flp = Ea—q—l
2
B (to
ORR = 8R3

2.
o 1 a ¢o+ laid} cottpa_¢0
69 stinchaez RBR R2 3‘9

 

2
_ 18% 13%
w - [Sigh—HEB?

 

 

 

 

2 -
O' _ la $0 ] a¢o 1
R9 RBRBO R239 sintp
2 -

_ a ¢0 t 3% l
09‘? _ 868$ C0 (pfi stintp
32¢ 64»

o __l o 1 0

 

(PR - Ramp RZBIp

57

58

for (I1

3
ZGuR = [Rsintpgc-gl — (3 -4v)sin(p¢1]cose

8¢
cosG + (3— 4v) (bl sin0

261.49 = a—e-

a
26% = [sintp—£1 — (3 —4v) cos<p¢1:| c050

 

 

 

aw
_ _2__,, act I 32¢ a¢ Lowery]
ORR — R_s——intp 8—0 sin0+ Rsintpfi— ——2(l v)sin(p§§l- -vZ R 8 0059
_ 2(1—v)a¢1 . [_ 1 324’1 a¢1costpa¢1:l
Gee — Rsintp 3—0 smO-I- Rs sin MPB—Gz +(1 —2v)sin(pa—-R + (1—2v)—R-% cosB
_ 2v all] [ 32¢] 132+¢1 costp_ a‘l’l
6W — RT—intpgé sinG— Rsintpg-k—HI- Rsintpa_92 -I-—(3 2v v) RB—tp

69¢
+(1+2v)sintpa—R]cosfl

2
8(1) 8¢13¢1
(Ike (1— —2v)aR sin0+Ia———Rae— 2(l—v )Ra—O]c050

13¢ ¢1 34>,
(l —Il?2v) #sin0+— RIB—7238-2 2(1—v)cot(p§§ cosG

09¢

2
3 <1) 3(1)] s_in(p 4’1
OIPR = [Sifflpm—(I- 2V) COS(paRl— “2(1 V) R —§;6 C059

59

for I),

a
2GuR= [Rsintpa—%2— (3— 4v)sintp¢2]sin0

a
= izsinO—(3 — 4v) ¢2cose

ZGu6 39

. 3‘92 .
20u,p = [smog-6 — (3—4v) costpoz] smO

2
2v a¢2 a 4’2 a‘1’2c_ostpa‘l’2 .
ORR — -R——sintp-a—0 2cosO+ [Rsintpa—zz— 2(1— v) sintpa—Rz— —v2 R _56 sm0

 

 

a a a B
2(1—v) ¢Zcose+[ l ”)2 +(1 —2v) smog—(f?2 + (1 —2v )L;S¢T:2]sin6

099 = Rsintp 30 Rs s1n (p862
aI 32¢ a2¢ BI
_ _2__v 2 _1_ c__osq)_ 2
CW - Rsincp 8—0 C089 [RSimPg—z+ Rsintpae2 +(3 2v v)— Ratp
3%
+ (l +2v) sinIp-a— 2:Isine
3$2 4’2 la¢2 .
0R9 = —(1—2v)aR 28ine+[:R—28_6—2(1_V)R-BE s1n0

59

2
la¢2 . 1 a (1’2 3% .
(p = —(1—2v)§§6 51n9+-[———2(1—v) coups-6 sm0

R 8684)

. a 4’2 8% sintpa¢2 ..
5ka = [gnome-W41—2v)cos<p§l—a —2(1—v)—R—% smO

60

for ¢3

3%
RCOS‘P8R3- (3—4v)coscp(113

ZGuR =
3%
26149 = cotcpfi

a
2Gu = cos1p3—¢3+ (3-4v)s1n(p(1)3

a Q) 8413 sin<pa¢3
ORR = RCOS¢-a-ITZ§-2(1—V)C05¢fi ‘12va

8413
- (l—2v) cosip-a—R-I-

G
<9
<1:

I

 

8% 32¢ sin 8%
(SW: (1— 2v)c031p§— +C(:(Pa _23 2(]_ ) Rwfi
2
3 _c__ot(p_ 8%
CR9 = COtwm—ZU _v)__ Ra—e
cotcp82¢3 1 [(1_2V) _szcplidja
“W ‘ me a9

8(1) 34) cos a¢3
M (l -2v)sin¢fi3+coscp9—Ifi6-2(1_ v) _RLPBTP

Q
ll

6]

for (114

. 9‘94
ZGuR = —s1n(p3—(5 —4(l —v) cosqxi)4

ZGue -

I
O

, 3% .
2Gu = RsmcpS-E +4(1—v)sm(p¢4

8% a 2¢4 sin (113%

O'RR —--2(2— V)Coscpa—R _Sinwa—M +2va

8¢4_ sin _(_pa¢4 4]

0'99 = (1-2v)[coscp (pa—R- R 81p

2

8(3) 8 R ' at)
= (1-2\1)coscp3§4 + sin¢51%+ (3 — 2v) 212.9354

2
_ 1 13434 v.cot__(_pa_¢4
O-"’9—2[R86841p+4(1 R 39
2
—1—a¢“ +4(1—v)la<p4
69¢ ‘ 2 M89 R55

2 2
l E) 4’4 8‘94 simp ¢4
6ch — §[:Rsin(paR2 +4(l —v) sincpa—R— ——(p:— (5—4v)

 

51mm”-

ZGuR

26149

7
«Gu‘p

0' :

(pR

- comp

62

for >13

2ax3
R36

(”13 COS ([3813]

= —2[sm(p-3—§ + 7%

2 8X3

7 1823 1823
‘ Rama—1356

 

 

_2 1813 +coma-7x; 1 37:3
RBRBG R2 8681p stinzgpae

7 compa—x} _COI2(.p8_)"3
“ R2 888$ R2 as

 

 

 

2 2 2

8 A3 18KB cotcpa A3 2C0t(pal3 ] a 7L3

Slnw —— + —— — + _ + .
8R3 RBR R BRBq) R3 aq) R251n2¢392

2 2 2
a A3_*_2cotq)_a_)_‘_3_simpa A3 + (2-sin21p) 313+ 2cotcp 8 A3
3R2 R BI R aRaI stincp Bcp Risinchae2

 

 

 

 

2 2
COt‘Pa 7‘3 + _1_8 A3 2C0t£pal3
R BRBB R38981p R3 86

 

 

APPENDIX B

displacement and stress fields in cylindrical coordinates

for 410
20“, = .33"
,.
8C)
26MB : £360
8
2Gu: = £0
0_ : 20
" ar-
0’ : lfi0+ifl
9" rar r3882
324
c__ = ——,)-O
~~ 9;...
2
G 13% 13%

I = $5755

 

_ la ‘1’

“6: ‘ raeaz
3240
05’ 2 8rd;

63

64
for I]

7 a¢l
.014 = r— — (3—4v)¢ c059
r Br 1

aI .
20149: 56 c058+ (3-4v)¢lsm8

80]
2014:: r— C088

8:
_ 2va¢1 . ab] Be] \
_ 7-8—6 s1n6+ |:ra—-—:2 —2v-a—r] cosG

_2(1—v)3¢1. 1M1 3‘1) .
Gee-+55 sm8+[-— 8T): +(l— —2v)ar c056

_ 2va¢1, 8— 8:6), 1311+ 3_ 3%] 9
G" — 756 sin r872 +— ra—e—z +( 2v)§; cos

- 1—2 84) 6 23¢ 1 18¢" 9
Gr6—( War sin + a———ae— 2( _V);§é cos

2
II, II.
59; = (l—2v)§-_- sin8+3——l—ea c036

— 811" 1 2 84, ‘6
Ozr— ra—ra'F—( - V)-a—E COS

 

in

65

for (1)2

[ ad)w ]
2011f: r:—;"—(3—4V)(i)2 sinG

 

am,

20149 = E‘sinG—(3-4v)¢2cose
a w

2014: = ra:"sin6

 

aII 32% 69¢
01: = —g—Y—‘c056+ [r "—2v-.— +1119

 

~ r 89 8:3 dr
_ ac, a2 a
099 = _Méfi-cosm [15% + (I —2v)aq: 2Jsine
224 924 32¢ 34»
_ _7-V 2C1 2 1 2 _ _2
or, — r —B—§C Os —[razz + r862 + (3 2V)ar :Icose
8‘92 82% 13¢2 .
Ore = —(l—2V)8—r cos9+ [ma—2U —v);fi smO
a a2
06; = —(l-2v)a—(:)2c039479—922z sine

6in II
on = I’rarT—E (l —2v)a—z ]sin0

66

for (113
343
20Hr = Z;-
za¢3
ZGUB = Pé—e—

BI
7—3—(3—4v) I3

 

 

2611: = :d”
_ 2‘43 7 aI3
0" ' “arz "' 8—:
2
_ 7 ‘86} za¢3+£a¢3
099 - —-VCOS(p-a—E +;§; r2563
32% 8%
6:: = 3.8:: —2(1-V)a—:
_ RI; z III
“re - Era—e779
32¢ aI
_ z 3 l__}
06: ‘ ;—aea‘: (1‘2V)rae
32¢ aI
— 7___£ _ _ __3
0v ’ ”are; U Mar

67

for (114
26 - 8%
ur —- r—a—Z-
2Gu6 = O

864
20%) = —r-a—r —4(1 —V) (1)4

2
a (1)4 6(1),;
on = r8r8: + (1—2v)3—Z-

 

'7 8%
099 = (1_~v)-é_i
3244 a 4
(5~ = —r—_-8r82 2(2 )EL
_ 13194
6’9 _ 2888:
_ 1 $4 18 4
“9 ’[arae 2“ ) as]
a 4 32¢ a
_ _l___4 __4_ _ _4
OJ — 2’862 razz 2(1 v)8r

68

for X3
2313
2014f — ;86
8X3
20119 — _2'a—r'
2Gu:=()

1813 1823
“I - 2 791579736

1313 1523
“II - ‘2 fire—972%

<5,_=0

2113 3323 2313
“r9 ‘ a— *?$ 72597

813
“I: — "37);?
2
18 A3
0’ :

3’ 7’888;

69

APPENDIX C

definitions of Legendre’s function and recursion formulas

 

 

PM -,,‘——‘-’,-,(x'-I)"
2 nldx
P:'(x) = (i—fl) denm
1'
Pm(r) = ———l——[(n—m+l)Pm (x) +(n+m)Pm (x)]
n ' (212+l).\‘ n+1 n—l
Pm (x) = —l——[(2n+ l)xP"’(x) — (n+m) Pm (x)]
n+1 "_m+1 n n—l
d "I 1 "I m
—P (x) = , [(n+1)xPn(x)—(n—m+l)Pn+l(x)]
dx n I l—x")
d2 1 d 2
-—3P:(x) : 7 [(2,1'3—jPZ1U)—{n(n+l)— m 2}an(x)]
dx l—x“ x l—x
2 ”2 m—l 1 m m
(l-X) P" (Y)='2_71—fi P,,+](x)—Pn_l(x)]
2 1/2 m—l ] m ,

2 ”2 m— l m
(l—x) P" l(.\') = [Pn+](X)-XP;,"(X)]

n+m

 

7O

de._ ‘21/2m+i
(l—x )al’nh) — (i-l) P"

(x) — meZ2 (x)

CHAPTER 2

EFFECTIVE PROPERTIES AND DAMAGE
FORMATION IN RANDOM COMPOSITE
MATERIALS: A COMPUTATIONAL APPROACH

1. INTRODUCTION

Industrial applications of composite materials are constantly growing, necessitat-
ing an improved understanding of their thermomechanical response. The heterogeneous
nature of composites indicates that their mechanical properties have to be understood basi-
cally at the microscopic level - this is the task of micromechanics.

Most research conducted in micromechanics concerns the local fields and effective

properties of composites and is based on the simplest case: a single inclusion problem.
Some analytical solutions were obtained. The early research involving heterogeneous
materials can be found from the classical textbooks of elasticity where some simple inclu-
sion problems, especially the hole problem are mentioned. For example, one of them writ-
ten by Muskhelishvili (1953) shows the solution of the plane problem containing a
CirCular inclusion under several loadings. A famous paper dealing with a single inclusion

is due to Eshelby (1957), who found that the stresses in an ellipsoidal and perfectly

b0nded inclusion are uniform under a constant loading.

7]

72

It is obvious that even for a single inclusion case, the micromechanics problem is
very complex. There is a large number of factors which influence the local and global
fields. These include the material constants of the constituents, the shape and relative size
of reinforcement (inclusion), the geometric arrangement, the boundary condition at the
inclusion-matrix interfaces, and others. As a result, only a few analytical solutions have
been obtained for some simpler cases.

Composites containing many inclusions, usually distributed in an irregular fash-
ion, pose great challenges making those analytical solutions practically impossible. For
this case, many theories or methods, such as generalized self-consistent scheme, self-con-
sistent scheme, differential scheme and Mori-Tanaka method are presented to predict
approximate effective elastic properties of composites, based on the fundamental solution
of a single inclusion. Generally. these methods give effective properties as a function of
the properties of each phase and the volume fractions of phases, and are applicable to rela-
tively simple condition, such as two phases and perfectly bonded interfaces, although
some solutions dealing with more than two phases have been reported. Comprehensive
reviews of literature in micromechanics are given by Christensen (1979), Hashin (1983),
Mura (1982).

Although the above mentioned methods show quite good approximation of effec-
tive elastic moduli of composites, solving for the local fields, and hence, the strength and
damage characteristics still poses a challenge. Here, we take note of two approximation
concepts - the periodic assumption of structure in composite materials and the numerical
solution for a discrete model of composites. First, geometric distribution of fibers (inclu-
Sions) in the composites is often approximated by assuming a periodic arrangement. This
Simplifies significantly the complicated problem and enables one to solve for the local and
effective fields by considering a unit cell. Regarding the numerical calculation, one
Choice, of course, would be the finite element method which is powerful and can handle

arbitrary geometry and properties of constituents. However there is another algorithm,

73

involving a so-called spring network, which is simpler and more efficient with regard to
time and expense in calculation.

The lattice algorithm is based on the idea introduced by Kirkwood (1939) and
Keating (1966) for the long wave method whereby the force on an atom is calculated in
terms of two nearest neighbor force constants 0t and [3 (central and noncentral, or
stretching and bending force, respectively) and atomic displacements. The relation
between the force constants and the elastic coefficients is obtained by comparing the
potential energy at macroscopic level and microscopic level. When the method is trans-
lated from physics to micromechanics, the force constants in the former become stiff-
nesses of springs in a spring-network, which is a discrete representation of the original
continuum model. The spring network model of composites can have a topology of trian-
gle, square, honeycomb, and other lattices. The algorithm is actually equivalent to a finite
difference method.

More and more, the lattice algorithm has been used for composite mechanics
because of the advancing computer technology. In most cases the triangular lattice was
used; Snyder, Garboczi and Day (1992), and Day et al. (1992) worked on effective elastic
response; Chen and Thorpe (1992) compared the elastic moduli of glass/epoxy fiber rein-
forced composites with the second and third order bonds; Sahimi and Goddard (1986),

and Beale and Srolovitz (1988) investigated elastic fracture in randomly inhomogeneous
materials; Chan et al. (1992), using honeycomb lattice, gave the elastic moduli and tensile
fracture stress as a function of the number of removed bonds in the network (randomly
distributed defects in materials).

It is interesting to note that Holnicki-Szulc and Rogula (1979) used this idea
inversely, i.e., to simulate a discrete engineering structure with a continuum model. In this
Way, they simplified the analysis of mechanics of a truss structure.

There are many other references in the literature taking advantage of the lattice

network. However, to our knowledge, no study has been conducted on the local fields and

74

the damage initiation and propagation in composites with randomly distributed inclusions.

In this second part of the thesis, we study the effect of random geometry in com-
posites on their effective properties and damage by using the above mentioned lattice
algorithm. The relation between elastic moduli of continuous material model and stiffness
of springs in the corresponding discrete model is derived in details in terms of tensors.
This is done for both two-dimensional elasticity and conductivity. Then, problems involv-
ing effective two-dimensional conductivity (equivalent to the out-of—plane elasticity) of
composite materials with either randomly or periodically distributed inclusion are solved
and the results compared with Hashin’s bounds on these properties. Finally, we study the
initiation and propagation of microdamage, which are simulated by breaking springs
(bonds) according to local stress concentrations in the spring network modelling the com-
posite material. The main focus is on the effects of random geometry on damage forma-
tion in composite materials; this will be studied for various parameters of constituents,
such as elastic moduli and fracture criteria, different system size and different random

configurations of composites.

2. COMPOSITE MICROSTRUCTURE AND THE
MESO-CONTINUUM MODELS

2.1 THE RANDOM COMPOSITE MODEL

By a random composite we understand a set B = {B((o); (o e (2} of deterministic
media 8(a)), where 0) indicates one specimen (realization), and Q is an underlying sample
(probability) space (Ostoja-Starzewski 1993a). All specimens B(0)) occupy the same
domain in x1,x2-plane. whereby we employ a two-dimensional (2-D) setting for the clarity
of presentation. Formally speaking. Q is equipped with a G-algebra F and a probability
distribution P. We specify P as a Poisson process with inhibition, which ensures that there
is no overlap of inclusions shaped as round disks. In particular, we use a sequential Pois-
son process and assume the disks to be occupied by a homogeneous isotropic continuum
of one kind, while the matrix by a continuum of another kind.

In setting up of the model we assume both phases to satisfy the so-called ellipticity

conditions: 3", n > 0 such that for any g the following inequalities hold for all the phases
C§§S§C§Sn§§ (1)

AS a result, we have a realistic ergodic medium model without holes and rigid inclusions
described by a random field Q = { Q (5, (o) ;)~c e B;(o e Q } with piecewise-constant real-
izations. It is clear that this piecewise-constant nature of stiffness fields is an obstacle with
employing the governing equations of continuum elasticity, which require that the stiff-
ness fields be once differentiable. Thus, there is a need for another, so-called, mesa-con-

‘immm model - one that possibly loses some information due to a “smearing - out”

76

procedure, but is sufficiently differentiable.

2.2 TWO SCALE-DEPENDENT RANDOM MESO-CONTINUUM
FIELDS

First, with the help of Fig. 1, we introduce a square-shaped window of scale

_E
5—d (2)

Equation (2) defines a nondimensional parameter 5 _>_ 1 specifying the scale L of observa-
tion (and/or measurement) relative to a typical microscale d (i.e. grain size) of the material
structure. 5 = 1 is the smallest scale we consider: scale of a fiber. In view of the fact that
the composite is a random medium, the window bounds a random microstructure 85 =

{85((0);(o e Q }, where 85(0)) is a single realization from a given specimen 8(a)).

 

 

 

 

 

Fig.1 Square - shaped window

A continuum-type constitutive law is obtained by postulating the existence of an effective

. I . I .
homogeneous continuum B?" ((0) of the same volume V5 (i.e. area in 2-D). whose

77

potential energy E, or complementary energy E*, under given uniform boundary condi-
tions equals that of a microstructure B5(o)) under the same boundary conditions. These

are of two basic types:

i) displacement-controlled (essential)

u (x) = $19.79. on 385 (3)

0 . . .
where g 18 a given constant tensor and 885 IS the boundary of B5,

ii) stress—controlled (natural)

I = Ogn. on 885 (4)

0 . .
where 9' IS a given constant tensor.

Boundary condition (3) results in an effective random stiffness tensor Egan).

with the constitutive law being stated as

9(a)) = cgtmg" <5)

(lo-dependence in (5) points to a random nature of the resulting stress field and of the effec-
tive stiffness tensor, with the fluctuations disappearing in the limit 5 —-> oo .
Condition (4) results in a following random operator form, involving a random

Compliance tensor. of constitutive law

Em) = 52(0))9 (6)

Hereinafter, superscripts '3 and " stand for essential and natural conditions, respectively.
Also, we use here the same type of notation for out-of-plane elasticity and conductivity,
Whereby g and g are vectors, Cij is an elastic stiffness (conductivity), and Sij is an elastic

compliance (resistivity).

78

Following (Ostoja—Starzewski, 1992, 1993a) we list here these principal observa-
tions:

1. Due to the heterogeneity of the microstructure B5((o), the inverse

dado) = [EMT1 ‘ (7)

is for any finite 8, in general, different from Qg obtained under essential conditions.

2. Q; (0)) and Q2 (0)) satisfy an inequality

gym) 3 gym) (8)

Hereinafter, for two second-rank tensors A and B, an order relation B 5 A means

11.81.1931,.ij V13: 0 (9)

3. In view of the spatial homogeneity of microstructure’s statistics, C; (to) and Q; ((1))
converge as 5 tends to infinity: this defines a deterministic continuum Bdet for a single

specimen 8(a))

C“"’(w> = ciao) = aim) (10)

whereby the window of infinite extent plays the role of an RVE of deterministic elasticity
(Conductivity) theory; in other words, it is at 5 —> oo that the invertibility of the constitu-
tive law is obtained. '

4. Ergodicity of the microstructure implies that

Cd“((o) = Qeff Vwefl (11)

where Q” is the effective response tensor (independent of 0)) of a homogeneous medium.

5- At any finite 5 both response tensors are, in general, anisotropic, with the nature of

79

anisotropy dependent on any specific 85(0)).

6. Since the window may be placed arbitrarily in the domain of B((:)), the essential and
natural boundary conditions define two different inhomogeneous tensor fields at the scale
5 with continuous realizations, which lead to two basic random continuum approxima-
tions: B; = {Bi-) ((0) ;(o e Q } and Hg = {B}: (to) , (o e Q }, respectively. Accordingly, a
window of size 5 may be considered as an RVE of these two random continuum models.
7. Our definition of two inhomogeneous tensor fields is conceptually similar - but not the
same — to the procedure of local averaging (Vanmarcke, 1983) in the theory of random
fields applied to a single realization C((o); me Q ; it becomes the same in case of a l—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 ener-
gies and, hence, the effective constants of finite windows and their probability distribu-
tions PL £3sz and PL Q2} Similarly, the autocorrelation (autocovariance) functions may be
determined (see also Ostoja-Starzewski and Wang, 1989 and 1990).

8. Principles of minimum potential and complementary energies can be used to obtain a

hierarchy of scale-dependent bounds on the effective stiffness tensor _Cfff

CR2 (5")"1 E (s’l’fl s <59" 3 <59" s 6% (cg) 5(Cg.)s(C‘1’)a CV (12)
5’<5

9. Since two different random anisotropic continua result, a given boundary value prob-

lem must then be solved to find the upper and lower bounds on response of composites.

3. THE FINITE-DIFFERENCE MODEL

3.1 LAYOUT

The nature of a composite microstructure - i.e. presence of many inclusions distrib-
uted at random in the specimen - precludes any simple analytical solution of the local
fields for this problem and of the crack propagation study in such a system. Therefore we
adopt a method analogous to those used in modelling of stiffness and/or failure of atomic
lattices (see Day et al., 1992; Beale and Srolovitz, 1988) and conduct the simulations in two
dimensions. Such investigations are based, in the simplest setup, on simulations of break—
down of a lattice with central-force interactions, in which some bonds are initially taken out
to represent the disorder, and from which cracks then develop in the course of simulation
(see Section 4 below).

To use the‘discrete model, a two-dimensional triangle lattice (Fig.2) is arranged to
cover the area occupied by the material. Suppose the material consists of cells of hexagon
shaded in Fig.1. All the cells have the same size, shape and orientation, and are centered

at vertices of the lattice. However, the material for each cell maybe different.

 

 

 

 

 

 

 

 

Fig.2 Two - dimensional discrete triangle lattice configuration

80

We consider a single cell (Fig.3), as a discrete model and describe the original con-
tinuum with six normal springs 0t" and six angular springs (or watch spring) [3". To
answer what material problems this kind of model can describe, we need to find the rela—

tion between the stiffness of springs and the elastic moduli of the represented continuum.

 

Fig.3 A discrete spring model of a single hexagon cell

3.2 MODEL FOR PLANE ELASTICITY
3.2.1 Alpha Model

As mentioned above, the model consists of straight and angular springs, but we
prefer to discuss them separately for the sake of clarity and simplicity. Consider a cell of
isotropic material as shown in Fig.3. A discretized model consisting of six straight springs
Which connect at the center of hexagon will be used to simulate the original material of
Continuous model. All springs are linear elastic and can rotate around the center freely.
The nth spring has an angle 6" with the coordinate axis XI and stiffness 0t", so the model
is called an at - model.

The relation of elastic properties between the continuum and the discrete model

Can be derived from the equivalence of energy of both models. For the continuum model,

the elastic energy E of the hexagon is expressed as

82

E = écjueijem i,j,k=1,2 (13)

where, C is the fourth order tensor of elastic moduli, the stiffness tensor. A is the area of

the hexagon cell. 81.]. is the linear strain tensor of the second order.

.,A c j

 

 

Fig.4 Normal spring

To get the energy expression for the discrete model, we note that 0L represents nth
spring of stiffness 0t" and length I (Fig. 4). X I and X2 are the rectangular Cartesian coor-
dinates. N is directed along the axis of the spring. We assume that point 0 is fixed and the
displacements of point L are u" , and u”?

The force P in the spring is

n n n
P (I "N v (14)

and

I!

”N = 1:14; (15)

where I] = cosB, I2 = sine are the directions of the spring.

Using (14 ) and (15). the strain energy in the spring can be written as

83

En__l_Pnun
2

1 nnun

=§an1i1ju Jul.

2 (16)

I
20L 11.11.211.128]

[antinnn

=2a 1111, U

The total energy of all six springs with the same length I is

12

6 6
2 15": [5’22 a"1’.'1’.'1"1"e ek, (17)

Alij
n=l

Comparing expressions (17) and (13) . the elastic moduli of a material represented

by the discrete spring model are found as

13 6
Cal!“ : “1;"le anlnlnlnln (18)
The area of the hexagon is
A = 27:312 (19)

Substituting the expression (19) into (18), we get

{7321a a"1’.'1’.’1”1" (20)

3n=l

CijkI:

It is obvious that this expression for elastic moduli satisfies the requirements of the
Syrnmetry of stress and strain tensors which is necessary. The result obtained here is

exactly the same as given by Holnicki - Szulc and Rogula (1979).

84

3.2.2 Beta Model
Again consider the single cell given in Fig.3, but only with six linear elastic angu-
lar springs used to constrain the free-rotation of the (1 springs. The stiffness of the nth

angular spring is denoted by [3", so the model is called a [3 - model.

K” /”

[J
E
N

Ni
a

 

 

Fig.5 Angle change, AB

Using Fig. 5, we want to find what is the angle change A8 when the point L
moves to the point L. through an infinitesimal displacement. Let f be the unit vector

along 0L and [be the length of 0L. It is seen that

ixa==ixi’=I-Aé (2n
Equation (21) leads to

-ixa=1gypj=emelr an

" l
A6=7 1 jptp

where, ski}. is the permutation tensor and a”, is a strain tensor, 1', j, p = l, 2, I1 = cosG

85

and I2 = sine.

The angle change between two segments is measured by (Fig. 6)

 

 

~ ~n l ~n
A4) = A9 + -A9 (23)
“2
n + 1 9n + 1 n
‘1’"
B 6n xl
0

Fig.6 Angular spring

Substituting equation (22) into (23). we get,

~ n + l n + l n n
An = emejpili 1p — I, 1],) (24)
For a single angle spring B". the strain energy due to Ad) is
I l r 2
E' = ,1) WI (25)

Substituting equation (24) into (25), we get

2
l n n+1 n+1 nn
Efl=§B {€,,,€,-,,(l,- 1,, ‘55)} (26’

86

Expanding expression (26) and using the identity

5. 8. —5. 5. (27)

..E :
8li km 1p M 1q JP

where 51.]. is Kronecker delta, and taking advantage of permutation function of 8i)" we

finally get,
_ B" ll n ll ’1 n
— 3%“ [PI 11p] 11.1]. Ikll
Olklz+lfjn+lfn+11n+l- ln+lljn+lln+lln+l (28)
-25Mggn+lf*‘f'+ 21f1"+|fl*‘/"}eueu

The total energy in a cell with six angular springs is

21%;u511nff—1f1wh"

(k p jpl
n=l
+8.A]’l+11"+11n+11n+l—l’-1+]l"+lln+lln+l (29)
i'p
25dq¢+lP*U"+21M+IF*U"}g;H

Where superscript (n+1) takes the value of 1 when n=6.

Equation (29) is really an expression of the strain energy of the unit cell of discrete
model. But to get the expression for elastic moduli, the symmetry with respect to i and j, k
and l, and ij and k1 should be satisfied. That means some subscripts in equation (29)
Should be permuted and the expression of the energy needs to be reconstructed so that it
Can be symmetric. There is no difficulty to work it out, but the desired expression will be

tedious for the general case.

For the cases we consider, i.e., the cases of two-dimensional isotropy and orthot-

87

mm. Only a few permutations are needed. In this case, equation (29) can be rewritten as

follows

E = g 2 {113" +Bn-1Jaikzynjzglf—[BHBn-‘szzyfizf

n=l

_Brzsik[;];z+ll;+l[7+Bnlylf+llz+117 (30)

n+ln+l n+1

-B”8iklnl'.*l 1, +B"I?HIJ’.*IZI, ”if“

P1P

The elastic energy of a continuous cell should be

Ezlcl’

2 queue/(IA (31)

Comparing equation (31) and (30), the elastic constants of the equivalent contin-

uum model can be simulated by the discrete lattice model according to

6

B_l (n "_j n nn(n —-)"n""

CW .. Z 2 { [3 +13 1 siklpznj1p1,— B +3" ' liljlkll
n=l

_an 1"1;+11;“17+ang‘ljwllf'zf (32')

i k p

n+ln+l n+1 nn+l

—B”5 ["121 1, +8”, 1121,}, }

ikpj P

where superscript (n— 1) takes the value of 1 when n=l.

3-2.3 The General Model
By superposition of equation (20) and (32), the elastic moduli of a cell of contin-

uUm are expressed in terms of stiffness of general discrete 0t — [3 model as follows

88

6
l
Cijkl— Enlanz
6
l n "_ n n _ n nnn
w; ([3 +13 Inlkllljsiklgljlgl" —(B +5" Illiljlkl,

n n+1 n n+ln
-B"5.-klp1j”'1p 17+B"lil}'.’+‘lk 1,

(33)

_Bnalklpljnl;+lll [n+1 +[3nln+llnln[n+l}

where A = 2../3I2

3.3 APPLICATION OF THE MODEL FOR PLANE ELASTICITY

3.3.1 Application to orthotropic material
For a two dimensional orthotropic material, there should be two axes about Wthl’l

the material is symmetric. To satisfy this condition, we take the coordinate axes as the
Symmetric ones, and choose (Fig.7)

(34)

(35)

 

Fig.7 Discrete model of a two dimensional orthotropic cell

89

Bl = B3 = B4 = B6 (36)
133 = 85 (37)

Using directly expression (33) with the conditions (34)-(37), the elastic moduli

corresponding to the or‘hotropic discrete model are evaluated as

 

 

C1112 = C2212 = O (38)

C1111 = fiUZal +iazj+ [13(331 +359] (39)
C1122 = 23/§[3a2_113(gfll+% 2)] (40)
C222: = $[3a3+l]—2(381 4352):] (41)
C1212 = 23/3[3a3+;§§81] ' (42)

Where subscripts in the tensor are simplified notation as general convention. The obtained

tensor of elastic moduli is symmetric.

90

3.3.2 Application to isotropy
3.3.2.1 Using both at and B model

For isotropic material, we assign all straight and angular springs the same values
of stiffness, respectively. In this case, the elastic moduli of a cell in (33) degenerate as fol-

lows

CU“- fix, I nlnlnln

3n=l

B n n n n n n
+2 2 {26,Ipl ,1p121,1,'11,1, (43)

n=l

_Blklzln + 11p

n+ln n+ln

1,271+II'HII, I,

n+ln+l n+1 nn+l

—.'511’.’I 1, +1. ["1], }

”\‘PJ 1P

Assume that the straight and angular springs are distributed uniformly, the length
of I is taken as a unit, the corresponding relation between the elastic moduli of the contin-

uum model and the stiffnesses of the discrete model is determined as

C1112 = C1211 = C2212 = C1222 = O (44)
C1111 = C2222 = 517.3(§a+112§ ) (45)
C1122 = C2211 = 2—1fi(ia‘;—2§B) ' (46)

C1212 = 2—j/_—3(ga+212-29-, ) (47)

In view of the listed results, it is obvious that the obtained elastic moduli in terms

91

of stiffness of the discrete model satisfy the requirements of symmetry due to stress and

strain tensors and strain energy. It is also observed that the condition

1
C1212 = §(C1111'C1122) (48)

is satisfied, so that there are only two independent constants and the proposed model can
really represent an isoti xpic continuum of isotropy.

Using the relation,

K+u = C” (49)

u = C66 (50)

3

l 3 l9
(awn—:13) ‘52)

n=—
2f3

Where [.1 is the shear modulus and K is the two-dimensional bulk modulus. It may be

r10ted that the bulk modulus due to the angular springs is identically zero, i. e.,

KB = 0 (53)

Th at implies the presence of angular springs never change the area of material.

The formula for a two-dimensional Poisson’s ratio is as follows,

V:——:— (54)

92

Substituting equations (51) and (52) into (54), we get

al2—3B

30:]2 + 3B (55)

From equation (55), we find the range of Poisson’s ratio which can be covered with this

model, i.e.,

i
8

(56)

q;

l
8

1:19 Ql‘m

Wlfl

Note that this is a two dimensional Poisson‘s ratio.

3.3.2.2. Using only the on - model
In this modeling without angular spring, the straight springs are chosen such that
all the six have the same length and stiffness, and all the angles between any adjacent ones

are equal to each other. In this case, equation (8) can be rewritten as

6
2 1f1f1217 (57)
n =1

01
Cum " 517-;
From equation (57), the elastic moduli are evaluated. The results show all the con-

ditions of isotropic elastic properties are satisfied, such as

C1111 = C2222 (58)
C1112 = C2212 = 0' (59)

l
C1212 = §(C1111‘C1122) (60)

93

The elastic constants are listed as

C1111: C11: “8—“ (6])
3

C1122 = C12 = 1g“ (62)
3

C1212 _ C66 = %a (63)

K = £01 (64)
4
u = iga (65)

where t1 is the shear modulus and K is two-dimensional bulk modulus.

The same expression as in equation (57) was mentioned by Holnicki-Szulc (1979).

The results equivalent to equations (64) and (65) were given by Keating (1966) and Kirk—
wood (1939), but the stiffness of springs which they took is half of what we take here so
that the coefficients in their expression would be twice the ones in equations (64) and (65).

As we did above, if we substitute equations (64) and (65) into (54), we get
V = 1/ 3, which means the above-mentioned model can be used only for the case of two-
dimensional Poisson’s ratio being 1/3.

To get Poisson‘s ratio different from 1/3, the parameters, such as stiffness of
SPrings or angles between springs may be changed instead of adding angular springs.
Chen and Thorpe (1992) employed the equilateral triangular lattice without angle springs,
but assigned three different stiffnesses to the straight springs. the range of Poisson ratio is
f r 0m 1/3 to 1. Thus, by combination of Thorpe’s model and ours, Poisson’s ratio in the

f“ 11 range for 2-dimensions, from -1 to 1, can be obtained.

94

3.4 APPLICATION TO CONDUCTIVITY

We now consider the two-dimensional conductivity problem governed by
Laplace’s equation. The model can also be used for several other problems, such as out -

of - plane elasticity, membrane problem, etc., which have the same governing equation.

3.4.] Triangle lattice model
Again we refer to Fig.2 , to the model with the normal springs. Following the same
procedure we used for the two-dimensional elastic problem but at lower level of tensor,

we get

a"(e,l;’)b (66)

Where it is out of plane displacement (temperature), 8 is gradient of displacement (gradi-

ent of temperature). The total energy of all six springs is

6

6
N l "I?"
E = 2 E = 52 a I,l,e,ej (67)
n=1 n=l

For the continuous model. the elastic energy E of the hexagon is expressed as

l or
E = ECU-8,8114 (68)

whore C ,j is the tensor of elastic moduli (conductivity). Comparing expressions (67) and
(68) and substituting the area of the hexagon (19), the elastic moduli represented by the

d1 Screte spring model are found as

6
C1". 1 2019717 (69)

Uzfi 1}

Assume the straight springs are chosen such that all six have the same length and

n=l

stiffness, and all the angles between any adjacent ones are equal to each other. In this case,

equation (69) can be evaluated as

_ J3
C11 C22 " ’2’“ (70)
CD = C2, = o (71)

It is obvious that obtained conductivity tensor satisfies the requirements of symme-

try and indicates the isotropy of the model.

3.4.2 Square Lattice Model

Fig. 8 shows the two—dimensional square lattice which is to cover the area which
materials occupy. Suppose the material consists of cells of square shaded in the figure. All
the cells have the same size, shape and orientation, and are centered at vertices of the lat-
tice and fill the lattice. However, the material for each cell maybe different.

We consider a single cell (Fig.9), as a discrete model, and describe the original
continuum with four straight springs 01" of the same length. Following the same procedure
as for the triangular lattice, we get the conductivity of the original continuum in terms of

the stiffness of the discrete model as

4
1 n ,
CU = 4 2 a If],2 (72)

n=l

96

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Fig.8 Square lattice Fig.9 Discrete square cell

Let all four stiffnesses 01 be equal to each other. The expression (72) is evaluated
as

C = C22 = ~01 (73)
C12 = C2, = 0 (74)

An orthotropic model can also be established by using square lattice if we set

1 2 3 4 . .
a , 0t equal to (1 , on respectwely. As a result, the model W111 be represented by

(75)

N1—
Q

to

(3
to
M

II

(76)

Nl—t

C,, = C2, = 0 (77)

4. EFFECTIVE CONDUCTIVITY

4.1 METHOD OF SOLUTION

To study the effective conductivity of a composite material, the above-mentioned
square and triangular lattices are used to discretize its domain. We consider two - phase
materials, i.e., composites containing inclusions, which are arranged either periodically or
randomly in the matrix. Regarding the latter, the structure of composite material is created
by generating, in a Monte Carlo sense, pseudo-random numbers for coordinates of the
centers of inclusions. A rounding off of these numbers is done to set the center of inclu-
sion exactly on the grid points. Then the stiffness of all springs (or, say, bonds) are
assigned according to whether they fall in the matrix or in the inclusion. Any spring strad-
dling the circular boundaries of inclusions has its stiffness ( 01 ) assigned according to a
series spring system weighted by the partial lengths ( 1,, 1,, ) of the springs that belong to

the respective domains, i.e.. by the harmonic average:

I] [11 "I
01 = [—- + —— (78)
01, 101,,
where I = total spring length
01,, 01,, = stiffness of partial spring.
The two-dimensional conductivity and the out - of - plane elasticity problems that
we are concerned with in this investigation have the same governing equations mathemat-
ically. For the two problems, the constitutive equation is

0,- (2) = C,,-(2) 8,- (2) (79)

97

98

where 0', is either heat flux for the conductivity case or stress for the out - of - plane elas-
ticity case. 8, is temperature gradient for the conductivity and strain for the out-of-plane
elasticity. C ,1. is the conductivity or, out-of—plane shear modulus for the two cases, respec-
tively. All of these are random variables depending on their local position.

In order to solve problems related to inhomogeneity, ‘effective modulus’ should be
defined. The formula used here can be found from a general textbook so that it is shown
here without more details. Based on the stochastic concept of random field theory, there
are two kinds of ‘average’, volume average and ensemble average. Using the concept of
volume average, we take a special specimen subject to a displacement boundary (essen-
tial) condition as follows

u = 5.x. (80)

where u is the displacement and temperature for out - of - plane elasticity and conductiv-
ity cases, respectively. x, is the position vector, 5, is the average strain or intensity which

is defined by

1
e, = AISidA (81)
A

Furthermore, the effective modulus can be defined as

_ _ eff—
o, — C 1'1 J. (82)
In this case, the potential energy is
E = 52c”; (83)

On the other hand, using the discrete model in section 2, the same potential energy

can be obtained by simply summing the energy in each spring as

01,,)(A11 (84)

~Nl-‘

where 01, is the stiffness of the 1111 spring in the model, and A11, is its deformation. n is
the number of springs in the model.
Equating the energy in (82) and (83), the effective modulus can be found. The fol-

lowing three tests are carried out to evaluate C H, C 22 and C12 respectively.

1) taking 8, = 1 , 82:0 leads to (85)
C11 = 2El (86)

ii)taking 8, = 0,122 = 1 leads to (87)
C22 = 2E2 (88)
iii)taking E, = 1,352 = 1 leads (89)

(90)

where E,, E2 and E, are the energies in the three tests. respectively.

4.2 RESULTS

We examine the effects of several factors on the conductivity. These include the
arrangement of inclusions, the type and size of lattice, the ratio of conductivity of inclu-
sion to matrix, the volume fraction of inclusions and randomness, and compare the
obtained results with their Hashin second - order bounds.

Specimens have inclusions arranged either periodically or randomly. In both two

100

cases, all the boundaries are periodic. The samples of composite configuration are shown
in Fig. 10. Only essential (displacement) boundary problems are solved. Following issues
are discussed,

(1) Since we represent a continuum system by a discrete lattice, the first important
question is how large should the diameter of the inclusion relative to the mesh spacing be
in order to adequately approximate the continuous material. We explore this issue by con-
sidering a classical case of effective conductivity of a material with a dilute concentration

of circular hole for which an exact solution is known to have the form

C"ff = (1 -—<.Rf)c"‘ (91)

where C‘ff is the effective conductivity of a composite, Cm is the conductivity of a
matrix, f is the volume fraction of the inclusion to the composite and :91 is a constant
which equals to 2 for a hole in dilute case.

A specimen of size 114x114 with square lattices is used to test the effect of the size
of an inclusion on accuracy of solution. The diameter of the single cavity in the matrix
changes from 2 to 20 spacings, the corresponding conductivities are obtained and the con-
stant 91 is calculated from equation (78). All the related data are given in Table 1 below.
It is seen that the inclusion having a very small diameter d = 2 or 4 does not represent the
continuum system well. As we increase the diameter, the accuracy improves and 9i con-
verges to the analytical solution, although there is a fluctuation. However, for large diam-
eters we have two competing effects - the larger the diameter the better the approximation,
but, at the same time, the dilute assumption ceases to hold. Note that we have computer
limitations on the lattice size. Finally we compromise by choosing d=l4 where the volume
fraction, close to 1 percent, corresponds to a dilute system. The effect of the hole size on
the effective conductivity is also shown in Fig. l 1. In Fig. 11 and the following Figs. 12 -

15 and 20 - 21, the calculation points are connected simply by the straight lines for clarity.

101

No regression has been made.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Table l:
dimer 1.211121 9‘ 92:55?

2 .00024 0.9997

4 .00097 1.7373

6 .00218 2.0356 +1.5
8 .00387 1.9783 -1.5
10 .00604 1.9078 -5.0
12 .00870 2.0541 +2.5
14 .01185 1.9848 -1 .0
16 .01547 2.0032 +0.1
18 .01958 2.0193 +1.0
20 .02417 1.9914 -0.5

 

 

(ii) In addition to (i), the mesh dependence is examined as follows. We vary the
size of specimens and inclusions simultaneously, but keep the volume fraction of the
inclusions unchanged. A volume fraction of 0.349 is used here. The sizes of the window
are L = 64, 127, 253 and 379, and the corresponding diameters of the inclusions are d =
14, 28, 56 and 84, respectively. 8 is fixed. In other words, we consider the dependence
of effective conductivities on the mesh. The investigation is conducted for various ratios,
0.0, 0.01, 0.1, 0.2, 0.5, l, 2, 5 and 10, of the conductivity of inclusions to those of matri-
ces. Fig.12 shows the relation between the ratio and the effective conductivity which is
also normalized by the conductivity of matrix. The results indicate that the effective con-

ductivities from the window size 64 and 127 coincide well, and the effective conductivi-

102

ties go up when the window size increases after that. The bigger is the contrast Ci/ Cm
in conductivities, the stronger is the scale dependence.

(iii) As mentioned in section 2, the conductivity is scale-dependent. To understand
this effect, we vary the size of window with L = 64, 127, 253 and 379. The diameter of
inclusion is still taken as 14 in order to ensure the accuracy of solution. The correspond-
ing values of scale 5 are 4.5, 9, 18 and 27, respectively. No bigger specimen can be taken
because of the computer limitation.

Fig. 13 shows the relation between the ratio of the conductivity of inclusions to
those of matrices and the effective conductivity which is also normalized by the conduc-
tivity of matrix. The results show that the effective conductivities coincide well for the
case of Ci/ Cm <10. The effective conductivities go down when the ratio L/ d increases.

Figs. 14 - 15 and Figs. 20 - 21 show the conductivity for the cases of L=64 and
128, respectively. Fig. 14 and Fig. 20 give the dimensionless effective conductivity (nor-
malized by the conductivity of matrix), C eff / Cm , versus the ratio of the conductivity of
the inclusion to the matrix, C i/ C' n and the volume fraction of the inclusion to the com-
posite, f, in two - dimensions. The volume fractions are taken as 0.03879, 0.1939, 0.3491,
0.4266, 0.5042 and 0.5818, respectively. The upper figure is the shaded effective conduc-
tivity surface and the lower one is the intersection of the above surface with the planes
where the volume fraction is fixed. To see more details, Fig. 15 and Fig. 21 give the
curves of the effective conductivity versus the volume fraction with ratio of the conductiv-
ity of the inclusion to the matrix being 100, 50, 10, 5, 2, 1, 0.5, 0.2, 0.1, 0.02 and 0.01,
respectively. In the figures, the curves for 11 cases are plotted in various forms of solid,
dashed or dashdot lines for the sake of clarity. The difference of conductivities obtained in
. these cases is very small, for most cases, in a range of 1 percent relatively, so that they
appear to be weakly scale-dependent for the range of the ratio of the conductivity of the
inclusion to one of the matrix, 0.1 - 10, which is more practically applicable to compos-

ites. In other words, the effective conductivity obtained from different 5 scale but contain-

103

ing inclusions of diameter 14 are comparable to each other.

(iv) Regarding the stochastic aspects, 10 samples are randomly chosen for each
case. The mean value and variance of the results obtained are calculated. The results from
two loading conditions, 8, = 1, £2 = 0 and £1 = 0, £2 = l are also averaged to elimi-
nate the macroscopic anisotropy because of the limited window size. The distribution of
the conductivity with different random arrangements of inclusions can be found in Figs.
16 - 19 which show the effective conductivities corresponding to various volume (area,
here) fractions for a certain conductivity ratio, C'/ Cm . In the figures, the dot marker
stands for the case of L=64 and the circular marker for 1.2127. The scatter of the results is
obvious but not significant, and the relative variance, for most cases, does not exceed 0.06.
It is noted that randomness plays a more important role when the ratio of conductivity of
the inclusion to the matrix is far from a unity, such as the case Cl/Cm = 0.02 in Fig. 16.

(v) As a special case, comparing to random ones, the effective conductivity of
composites with periodically arranged inclusions are evaluated for the same volume frac-
tion and stiffness ratio. Fig. 24 and 32 display for two different stiffness ratios the stress-
strain curves of a periodic composite and of several random composites. It is seen that
the effective conductivity for the periodic composite is within the range of scatter of the
effective conductivities for the random composites, i.e., although the stress-strain curves
appear different. the initial slopes, i.e. conductivities, are almost the same. It should be
noted, however, that this result has been obtained for a rather restrictive model: a periodic
window with the dimension only about 3.5 times larger than the diameter of the inclusion.

Ceff calculated via CE; and C; (see Section 2.2).

This is not the case with

(vi) The range of volume fraction of inclusions is from zero to about 0.6, the max-
imum value that can be obtained in the case of non-overlapping inclusions. For the above
test, the volume fractions are taken as 0.03879, 0.1939, 0.3491, 0.4266, 0.5042 and
0.5818, respectively. Figs. 16-19, where the solid lines are the Hashin second upper and

lower bounds for the effective conductivity of composite in two-dimensions, show that the

104

effective conductivities for various volume fractions are consistent with the Hashin sec-
ond - order bounds.

(vii) A sample is also made of for the ratio of conductivities of the inclusion to the
matrix from 0.01 to 100 which covers a range wide enough to be applicable to most prac-
tical composite materials. The results show that all the values of the effective conductivity
fall within their feasible range enclosed by Hashin second—order bounds in Figs. 16-19.
Furthermore, the values approach the lower bounds as the ratio is larger than 1, and the
upper bounds for the ratio less than 1. This implies the model lattice has a little lower
stiffness, since according to the theory in (12), all conductivity should be near the upper
bounds

(viii) In order to understand the sensitivity of conductivities to the type of lattice,
both triangle and square lattice are used. The results from the two different types of dis-
crete models are quite consistent for the same composite material. No more comparison of
lattice types will be mentioned, it is not the main concern in our study.

In this section, the effect on the conductivity of many facts, such as the size of
mesh, the ratio of the window dimension to the diameter of inclusions, the composite con-
figurations (randomly and periodically arranged inclusions), the ratio of the conductivity
of inclusion to the conductivity of matrix, the volume fraction of inclusions and the lattice
shapes are examined. The inclusion of diameter 14 and the window of 64x64 are used to
get reasonable conductivities with feasible expenses on HP or SUN work station. The cal-

culated conductivities fall within Hashin second order bounds.

5. DAMAGE AND FRACTURE SIMULATION OF
COMPOSITES

5.1 METHOD OF SOLUTION

The investigation of damage and fracture in composite material is based on a sim-
ple and straight - forward idea, i.e., simulation of breakdown of a spring lattice model, in
which heterogeneity is introduced by, for example, taking out some springs to represent
the disorder and from which cracks then develop in the course of simulation.

In this part of thesis, we assume that the microstructural response of composites is
elastic-brittle. In other words, the springs in the lattice will have linear elastic property
until they break according to the established criterion. For the two phase material, there
are two kinds of springs representing matrix and inclusion, respectively. Therefore, the
composite material can be characterized by four parameters. Of them, Cm and Cl are
stiffnesses of springs for inclusions and matrices, and 8:: and 8:, are critical strains at
which springs will break. In this notation, the superscripts m and 1' stand for the matrix and
inclusion, respectively, and cr for ‘critical’. The global response of a composite is, how-
ever, a result of cooperative spatio-temporal phenomena of breaking along with the stress
and strain redistribution in the microstructure - a process which is simulated as described
below.

The fine mesh model of the matrix-inclusion composite forms the basis of compu- '
tational simulation of damage formation. These are conducted in the following steps.

(1) We subject the boundary of a square—shaped lattice to kinematic boundary con-

105

106

ditions

u = ijj j = 1,2 (92)

where g = (8,, 82) is the macroscopic (applied) strain (for out of plane elasticity case

where E E E3, and £2 a £32 ) or intensity (for conductivity case). In the numerical

l

examples in this paper we take 8, = E and £2 = 0 .

(ii) We simulate the increasing loading conditions through raising E by small
increment A8 in every next run of the simulation.

(iii) In every run we use a relaxation algorithm for the square lattice (conjugate
gradient algorithm for the triangle lattice) to solve for the equilibrium of the lattice, and
next, we search for all springs, in either matrix or inclusions, whose strains exceed the
local fracture criterion

m m i i

8 >8 and 8 >8 (93)

cr CI’

For any spring straddling the circular boundaries of inclusions having a weighted stiffness
assigned, a harmonic averaged fracture criterion is set up also and a criterion similar to
(93) is used. If the above criterion is met, given springs will be removed from the lattice -
thus representing a crack - and the microscopic strain E is increased according to step ii).
(iv) The increase of E by A8 is conducted by first unloading the entire lattice, and
then reloading it by strain E + As .
(v) We repeat the above steps until a continuous crack is formed through the whole

specimen (state of crack percolation).

108

the terminology of Section 2, a single specimen B5((0) at 5 E 4.5 and o) = (1)], whose'
overall stress-strain response is given in Fig. 24. Thus, a question arises as to how strong _
is the scatter among specimens of the same type. To that end we present damage evolution
patterns (Figs. 25 and 26) and corresponding stress-strain curves (Fig. 24) for two other
configurations indexed by (02 and (03, respectively. These three stress-strain curves are
marked with a, b and c in the figure, respectively. Fig. 24 shows also the response curves
of the inclusion material and matrix material per se.

Following observations are now in order:

- Scatter in effective moduli in the pre-damage region, as shown in the last chapter, is
much smaller than that which sets in after damage starts taking place.

- As expected, qualitatively same damage phenomena occur. A big space of 20 samples
will be presented later for statistical study.

— The stress-strain curves are all of the same type - i.e. pre-peak - albeit the scatter in
their shapes and the strain - to - failure is significant.

- Presence of inclusions - both random and periodic - has a strongly weakening effect
on the overall strength of a composite as compared with the strength of the inclusion mate-
rial or matrix material; this is due to the strain/stress concentrations.

- Spatial randomness in the distribution of inclusions has a strongly weakening effect on
the overall strength of a composite as compared with the strength of an idealized periodic
composite, Fig. 24; note that this observation is highly dependent on the particular mate-
rial parameters - for example, the case Ci/C’" = 5 and sir/E: = 5 exhibiteda much
smaller effect.

Up until now the randomness was due to nonuniform geometric arrangements,
while the properties of matrix and inclusions were taken as homogeneous - in reality each
material possesses some spatial random inhomogeneity in elastic moduli and strength, and
thus the damage patterns and the effective response of an actual periodic system would

differ from the one in Fig. 24 because of damage localization; those effects are examined

109

later.

(iii) The issue of material microstructure randomness effects on the overall com-
posite properties can further be clarified by a comparison with a composite of exactly
same volume fraction (i.e. 0.349) but a perfectly periodic distribution of inclusions, the
results from which configuration are shown in Fig. 24 and Fig. 27. The former contains
the entire stress-strain response of such a composite, and the latter illustrates the evolution
of damage.

Once again, in view of Fig. 24, it is seen that the effective conductivities in both
random and periodic cases are almost identical. However, the stress-strain curve for the
periodic case extends much longer than for the random cases, and there are some differ-
ences between the random cases themselves. The smaller strength in the random case is
caused by the stress concentration due to the nonuniform distribution of inclusions. This
can also be observed from Fig. 27 where the deformation appears uniform, and the dam-
age happens periodically in the configuration at same time so that the composite can
undertake higher loading than in the random case. Also, the fact that the breaks occur uni—
formly in Fig. 27 for the case of a composite with the periodic arrangement of inclusions
shows our algorithm and program work well.

(iv) In order to better understand the effect of randomness of special arrangements
of inclusions on the overall mechanical properties, the other case with the different stiff-
ness ratio and strength ratio, C‘/C’" = 2 and air/8:: = 0.2, is illustrated here. Figs. 28 - 30
and Fig. 31 are the damage evolution graphs for three random samples and the periodic
one. The corresponding stress-strain curves are shown in Fig. 32. From these graphs the
following conclusions are drawn.

The stress-strain curves are all of the same type, although the scatter between them
in their strain-failure is noticeable. But different from the case in the last section, the
effective response of the composite, both random and periodic, is of the elastic-brittle type

with feature of one peak only in the curve.

110

As expected, qualitatively same damage phenomena occur, but the idealized peri-
odic composite leads to an overestimate of the effective strength: geometric randomness
has a weakening effect. However, different from the last case, here the damage forms ini-
tially in the denser region due to stress concentrations in their vicinity, stiff inclusions tend
to form links carrying relatively more load. Microcracking then spreads across the speci-
men whereby the random character of the damage pattern reflects the heterogeneity of the
microstructure.

Presence of inclusion, both random and periodic, has a weakening effect on the
overall strength of a composite as compared with the strength of the inclusion material or
matrix material, due to the stress and strain concentration

(v) Comparing to the case in section (ii), the effect of changing sir/8:: while
keeping the stiffness ratio 076"" fixed is shown in Figs. 33 and 34; the first one gives the
damage evolution pattern, while the second presents the stress-strain curve; note that the
same geometry as in Fig. 23 is being used. We see that:

- Although the inclusions are softer than the matrix, just as in the previously studied case,
their now higher strength causes damage to initiate in the matrix.

- Cracking occurs within the inclusions at a later stage, and is seen to be needed for a total
breakdown of a composite. Comparing with the case in the last section, the different
strength ratio leads to the different damage processes.

- The stress-strain response is of the pre-peak type.

(vi) A question may be raised at this point as to the effect of a small change in com-
posite’s geometry on both the damage pattern and the overall stress-strain response. Such
a change is present in the composite of Fig. 35 as compared to the original one of Fig. 33.
It is seen from the ensuing damage evolution that such small microstructural variations can
lead to some differences in crack patterns and thus to nonnegligible changes in stress-strain
curves as shown in Fig. 34 in which the curves for the two different configuration are
marked by d and e, respectively.

(vii) At this stage we return to the aspect of mesh dependence in our simulation

method. Some light on this very broad and important issue is shed by simulations of

lll

exactly same composite as the one of Fig. 35 at double the mesh spacing per inclusion as
presented in Fig. 36. The corresponding stress-strain curve is marked by f. A comparison
of these two figures indicates that cracks form and evolve in almost same locations at two
different mesh resolutions. However, this is a delicate issue which requires a further study;
see also (de Borst et a1, 1993). The stress-strain response curves given in Fig.34 show that
the qualitative features are preserved under a change of mesh.

(viii) Furthermore, we are concerned with the effects of intrinsic randomness of
geometric arrangements of inclusions on the overall stress-strain response of composite in
the entire parameter space which is specified by two parameters, as used previously: the
stiffness ratio Ci/C'" and strain-to-failure (strength) ratio sir/8:: . In the following we call
the plane defined by these two parameters a ‘damage plane’, and a representation of a given
property, such as effective stress-strain curves, a ‘damage map’.

Fig. 37 shows the damage map of stress-strain curves of composites and the related
matrix and inclusion materials for the all nine combination of the two parameters which
are valued as 0.1, 1 or 10. Fig.38 shows the same damage map without the curves for
matrix and inclusion materials so that a bigger scale can be used to get more details. It is
observed that from the left of the map to the right, from the bottom to the top, i.e., as the
strength ratio and stiffness ratio increase, the type stress-strain curve changes from the pre-
peak to the elastic-brittle one. Correspondingly, the break initiates from inside the inclu-
sion to inside the matrix. In the former case, cracks initiate in the weak inclusion, then stop
growing until a high level of load is reached due to the difference between the properties of
matrix and inclusion materials. The cracking procedure is discontinuous so that the stress-
strain curves are of the pre—peak type. In the latter case, cracks develop through the matrix
material continuously and the stress- strain curve is of the post-peak (or elastic—brittle) type.

Fig. 39 is the corresponding damage pattern map. It is clear that from the left side,
fluctuation of location of the figure to the right side, from the bottom to the top, the crack
initiates from inside the inclusion to inside the matrix.

(ix) As stated earlier, models of uniform materials are usually used for analysis,

however, real materials are always heterogeneous due to fluctuation of geometry or pr0p-

112

erties of material’s constituents. The influence of geometry is considered here.

Figs. 40 - 41 illustrate the effect of fluctuation of inclusion location on the damage
pattern map and the stress-strain curve for the parameters sir/8:; = 2 and Ci/C'" = 0.2,
respectively. In Fi g. 40, the left column represents the evolution of the idealized periodic
arrangement of inclusions, while the right column the case of that the inclusion at the third
column and the third row is moved up one unit. Although the cracks initiate in the matrix
in both cases, only one crack, instead of many in the idealized case, forms in the vicinity of
the inclusion which is moved. It indicates the damage pattern is quite different due to the

stress concentration.

Fig.4] shows the stress-strain curve also fluctuates due to the concentration and the

composite is weakened. ,

Figs. 42 - 43 illustrate the cases with two different sets of parameters: air/8:1, = 2,
Ci/C" = 5 and sir/8:: = 0.5, 076'" = 0.1. Similarly as before, the cracks are not uni-
formly distributed any more. It should be noted that for the first two case is cracks always
occur in the matrix no matter if the location of the inclusion moved, but for the last case
the crack which initiates in the inclusion for the idealized composite, will start in the
matrix responding to the small change of the inclusion location.

(x) In addition to geometry, we examine the effect of fluctuation of matrix materi-
als on the damage pattern and stress-strain curves which are shown in Figs. 44 - 49. At this
point, we assume that the strength of matrix deviates a little from the uniform status as fol-

lowing

a : 8C7+AECT

_ (l r) (94)
_ Ecr +1?

where, 8" = strength of uniform material
A8,, = strength fluctuation

r = uniform distribution random number, [-0.5, 0.5]

N = weight factor, 100 in this case.

The parameter pairs chosen for the material fluctuation are: sir/8:: = 2, 076’" =
0.2; sir/8:: = 2, C‘/C’"= 5 and sir/8:; = 0.5, C‘/C’" = 0.1 for three different cases. Com-
paring with the uniforrn materials, the figures show that the damage of composites is not
uniform any more due to the stress concentration caused by the existing fluctuation. The
appearance and location of cracks are very much similar to those without material fluctua-
tion, but the number of them reduce to one. Figs. 45, 47 and 49 indicate the stress-strain
responses with or_without material fluctuation are quite similar to each other, the stiffness
is unchanged and only the strength of composites is reduced.

(xi) All previous studies are conducted by using one or a few samples for each case
of randomly arranged inclusions. The problem to be deal with is: what is the response of
stress-strain relation to the random configuration of inclusions? For this purpose, the Q
space containing 20 samples is employed for statistical study. The numerical characteris-
tics are also included. Figs. 50, 53, 56, 59, 62, 65 and 68 are assembly of effective consti-
tutive responses. Figs. 51, 54, 57, 60, 63, 66, 69 and 52, 55, 58, 61, 64, 67, 70 give the
distribution function of maximum stresses and strain which are normalized by the corre-
sponding matrix properties, respectively.

Fig. 71 gives the damage map of constitutive curves. The pair of parameters,
Sir/8:: and 076"” are used. The results from the 20 samples for each different pair of
parameter show the type of constitutive curves may change by chance. Fig. 72 and 73 give
the damage map of maximum stress distribution function and maximum strain distribution

function, respectively.

5.3 CLOSURE

As expected, the statistical analysis reported here shows that the conductivities of
composites are less dependent on the randomness of configuration. In contrast to this, big
deviations of maximum stress or strain are found, which suggests that the stochastic study

is necessary and important to determine properties of heterogeneous materials.

114

We would like to conclude with the following statement. The discrete lattice and
continuum models have quite different - if not opposing, and hence mutually complement-
ing - features. Thus, for example, the discrete lattice approach as presented here enables
one to capture the global aspects of randomness, while the continuum based approach is
suited to conduct a detailed study of interaction of a crack with a single inclusion, see e. g.
Li and Chudnovsky (1993). The damage initiation and propagation in composite materials
cannot be adequately addressed without the incorporation of geometric and material ran-
domness, and thus the present study illustrates a need as well as a possible basis for a sto-
chastic continuum damage mechanics of such systems. The damage information obtained
in this study needs to be compared with experiments and results of other numerical meth-

ods.

REFERENCES

Beale, PD. and Srolovits, D.J., 1988, “Elastic fracture in random materials,” Phys. Rev.,

Vol. B37, pp. 5500-5507.

Brockenbrough, J.R., Suresh, S. and Wienecke, H.A., 1991 , “Deformation of metal-matrix
composites with continuous fibers: geometrical effects of fiber distribution and
shape,” Acta Metal]. Mater, Vol. 39 (5), pp. 735-752.

Carmeliet, .l. and De Borst, R., 1994, “Stochastic approaches for damage evolution in
quasi-brittle materials,” Probabilities and Materials - Tests, Models and Applica-

tions, NATO ASI Series E - Vol.269, pp. 491-503, Kluwer, Dordrecht.

Chan. 8., Kim, 8.6. and Duxbury, P.M., 1992, “Elasticity and fracture of tubular honey-

combs containing random vacancies,” Unpublished work.

Chen, J. and Thorpe, M.F., 1992, “Fiber aligned glass/epoxy composites with second and
third order bounds,” Unpublished work.

Christensen, R.M., 1979, Mechanics of Composite Materials, J. Wiley & Sons, New York.

115

116

de Borst, R., Pamin, J., Schellekens, J.C.J. and Sluys, L.J., 1993, “Continuum models for
localized failure,” Proc. IUTAM Symp. on Fracture of Brittle Disordered Materi-

als.

Day, A.R., Snyder, K.A., Garboczi, EJ. and Thorpe, M.F., 1992, “The elastic moduli of a
sheet containing circular holes,” J. Mech. Phys. Solids, Vol. 40(5), pp. 1031-1051.

Duxbury, P.M., Kim, 8.6., and Leath, PL, 1994, “Size effect and statistics of fracture in
random materials.” Mat. Sci. and Eng. A, A176, pp. 25-31.

Duxbury, P.M., 1990, “Breakdown of diluted and hierarchical systems,” in statistical
Models for the Fracture of Disordered Media, Herrmann, HI. and Roux, S., eds.,

Elsevier/North-Holland, pp. 189-228

Eshelby, .1 .D., 1957, “The determination of the field of an ellipsoidal inclusion and related

problems,” Proc. Royal Soc. Lond. A, Vol. 241, pp. 376-396.

Hashin, Z., 1983, “Analysis of composite materials - A survey,” J. Appl. Mech., Vol. 50,
pp. 481-505.

Hen'mann, HI. and S. Roux, S., eds., 1990, Statistical Models for Fracture of Disordered

Media, Elsevier.

Holnicki-Szulc, J. and Rogula, D., 1979, “Nonlocal, continuum models of large engineer-

ing structures,” Arch. Mech., Vol. 31 (6), pp. 793-802.

Huet, C., 1990, “Application of variational concepts to size effects in elastic heteroge-

117

neous bodies," J. Mech. Phys Solids, 38 (6), pp. 813-841.

Jagota, A. and Bennison, 8.1., 1993, “Spring-network and finite-element models for elas-
ticity and fracture,” in Breakdown and Non-Linearity in Soft Condensed Matter,

Bardhan, K.K., Chakrabarti, 8K. and Hansen, A.,eds., Springer-Verlag, in print.

Keating, RN, 1966, “Effect of invariance requirements on the elastic strain energy of
crystals with application to the diamond structure,” Phys. Rev., Vol. 145, pp. 637-

645.
Khang, B., Batrouni, G.G., Redner, S., de Arcangelis, L. and Herrmann, H.J., 1988, “Elec-
trical breakdown in a fuse network with random, continuously distributed breaking

strengths,” Phys. Rev., Vol. B37, pp. 7625-7637.

Kirkwood, 1.6., 1939, “The skeletal modes of vibration of long chain molecules,” J.

Chem. Phys, v01. 7, pp. 506-509.

Krajcinovic, D. and Lemaitre, 1., 1987, Continuum Damage Mechanics - Theory and

Applications, CISM Lectures, Springer-Verlag.

Lemaitre, J., 1992, A Course on Damage Mechanics, Springer-Verlag, Berlin.

Li, R. and Chudnovsky, A., 1993, “Variation of the energy release rate as a crack
approaches and passes through an elastic inclusion,” Int]. J. Fracture, Vol. 59, pp.

R69-R74.

Maugin, G.A., 1992, The Thermomechanics of Plasticity and Fracture, Cambridge Uni-

118

versity Press.
Mura, T., 1987, Micromechanics of Defects in Solids, Martinus Nijhoff Publishers, Hague.

Murzewski, 1., 1957, “A statistical theorey of brittle quasi-homogeneous solids,” Proc. 9th

IUTAM Congress, Vol. 5, pp. 313-320.

Muskhelishrili, N .I., 1953, Some Basic Problems of the Mathematical Theory of Elasticity,
Noordhoff.

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

pp.6l-80.

Ostoja-Starzewski, M. and Wang, C., 1990, “Linear elasticity of planar Delaunay net-
works - II: Voigt and Reuss bounds, and modification for centroids,” Acta Mech,

Vol. 84. pp.47-61.

Ostoja-Starzewski, M., 1992, “Random fields and processes in mechanics of granular
materials,” in Advances in Micromechanics of Granular Materials (Proc. 2nd US-
Japan Seminar on Micromechanics of Granular Materials; H. Shen et al, Eds.);
Stud. App]. Mech., Vol. 31, pp. 71-80, Elsevier; extended version in Mech. Mater.,
Vol. 16 (1-2), pp. 55-64, 1993.

Ostoja-Starzewski, M., 1993a, “Micromechanics as a basis of random elastic continuum

approximations,” Probabilistic Engng Mech., Vol. 31 (2), pp. 107-114.

119

Ostoja-Starzewski, M., 1993b, “Micromechanics as a basis of stochastic finite elements
and differences - an overview,” App]. Mech. Rev. (Special Issue: Mechanics Pan-
America 1993, M.R.M. Crespo da Silva and C.E.N. Mazzilli, eds.), Vol. 46 (11,
Part 2), pp. 8136-8147.

Ostoja-Starzewski, M., Alzebdeh, K. and Jasiuk, 1., 1994, “Linear elasticity of planar

Delaunay networks - III: Self-consistent approximations,” Acta Mech., in press.

Ostoja-Starzewski, M. and Jasiuk, I. eds., 1994, Appl. Mech. Rev. (Special Issue: Micro-
mechanics of Random Media), Vol. 47 (1, Part 2).

Ostoja-Starzewski, P.Y. Sheng and I. Jasiuk, 1994, “Influence of random geometry on
effective properties and damage formation in composite materials,” ASME Journal

of Engineering Materials and Technology, Vol. 116, pp. 384-391.

Ostoja—Starzewski, P.Y. Sheng and l. Jasiuk, 1994, “Micromechanics as a basis of stochas-
tic continuum damage mechanics,” in Material Instabilities: Theory and Applica-
tion (R.C. Batra and HM. Zbib, Eds), AMD-Vol.183, MD-Vol.50, ASME-Press,
pp. 1 31-141.

Ostoja-Starzewski, P.Y. Sheng and I. Jasiuk, 1994, “Influence of random fiber arrange-
ment on crack propagation in brittle-matrix composites,” Proc. 4th Conference on

Brittle-Matrix Composites (A.M. Brandt, ED. ), Warsaw, Poland, pp. 200-207.

Sahimi, M. and Goddard, J.D., 1986, “Elastic percolation models for cohesive mechanical
failure in heterogeneous systems,” Phys. Rev, Vol. B33, pp. 7848-7851.

120

Schlangen, E., 1993, “Experimental and numerical analysis of fracture processes in con-

crete,” Ph.D. Dissertation, Technological University of Delft.

Snyder, K.A., Garboczi, 13.1., and Day, A.R., 1992, “The elastic moduli of simple two-
dimensional isotropic composites: Computer simulation and effective medium the-

ory,” J. Appl. Phys, Vol. 72 (12), pp. 5948-5955.

Vanmarcke, EH, 1983, Random Fields: Analysis and Synthesis, MIT Press, Cambridge,
MA.

Volkov, SD, 1960. Statistical Theory of Strenth of Materials, Mashgiz, Moscow.

Weibull, W., 1938, A statistical theory of the strenth of materials, Ingen. Vetenskap. Akad.
Handl, Vol. 15, pp. 1-45.

Ziegler, H. and Wehrli, Ch., 1987, “The derivation of constitutive relations from the free
energy and the disspation functions,” in Adv. App]. Mech. 25, pp. 182-238, Aca-

demic Press.

121

 

 

 

 

 

Fig. 10 Composite configuration: periodic and random arragement of inclusions

122

 

1.1 1 1 1 1

511/2

 

 

 

 

0‘ l l 1 l
0 0.005 0.01 0.015 0.02 0.025

Fig. 11 Influence of hole size on effective conductivity

 

 

window size

o----64

+----127
x----253
*—--—379

 

 

 

-1 II

10 10

9519...

Fig.12 Mesh dependence of effective conductivity

124

I I IIIIIII I I IIIIIIT I I IIIIIII I I IIIII

 

 

 

l lJllIJll l l llLllIl l l lllllll l l llllll

 

1

10 to" 10' 10 ~ 1112

Ci/Cm

Fig. 13 Delta dependence of effective conductivity

125

 

 

Fig. 14 Effective conductivity, at window size L=64

126

 

Ceff

 

 

 

 

 

 

Fig. 15

 

Effective conductivity versus volume fraction, at window size L=64

127

 

so I I I I 1' T I I T

45"

0L=64
50 OL=127

mm
q
a
9319..

 

 

 

 

 

9.19..

 

 

 

 

 

Fig. 16 Effective conductivity and Hashin second order bounds versus

volume fraction, C i/ Cm = 0.02 and 50

 

Cé’ff

 

Cm

128

 

 

 

 

 

 

oL=64
oL=127
10
I
DE 0: of» QJB Of
V1
.L=64
oL=127
Ci
C’"
03 a; a} o} 0}

V1

1
0.. 0.9

Fig. 17 Effective conductivity and Hashin second order bounds versus

volume fraction, Ci/Cm = 0.1 and 10

 

 

129

 

 

 

 

 

 

 

 

Fig. 18 Effective conductivity and Hashin second order bounds versus

volume fraction, c‘/c’" = 0.2 and 5

 

130

 

 

 

 

 

 

 

 

2 r t t T r
I" h _ o L = 64
a,” C' _2 oL=127
.1- . 6’"
(ff 6 7
C
a, -
Cm
¢ .-
..,, _
.3 _
. D D 1;! 0:8 054 of» 0.5 D? 1' 0:8 0.19
V1
0'” - o L = 64
as r O L = 127 i
—C— = 0.5
0.36 )- Cm
Ceff as -
m 0.!!! 1'
C
D. f '-
DIvb "
0.1:- ’
D.bb '-
u. . . . . . . . . .
D Q 2 0.8 0.4 [lb R5 'D.( 0.3 n,
Vr
Fig. 19 Effective conductivity and Hashin second order bounds versus

volume fraction, C 1./ Cm = 0.5 and 2

 

 

131

 

 

100

 

Fig. 20 Effective conductivity, at window size L = 128

132

 

I

3.5

Ceff

 

0.5-

 

 

 

 

Fig. 21 Effective conductivity versus volume fraction, at window size L = 128

133

 

 

Post-peak type

 

 

Pre-peak type

Fig. 22 Two kinds of constitutive responses

134

 

Fig. 23 Damage pattern for random arrangement of inclusion,

s‘/e"’= 0.5, c‘/c’" = 0.1

135

 

 

 

 

 

 

 

0.35 I I I I I I
0.3 r .
, periodic
matr1x
0.25 - I .
G /' / 1
—m r'J/A' l
Ger l.’ i 1 I
0.2 - ,5 " I | .
,1” I l
/ z ' l
. l
/’//P i 1
11.1 5 - 1-31” - 1 -
//l I 1
motor)
. ' I
01 - ’1', a :c -
" 1
l
l
0.05 ' inclusion I ‘
|
I
. l
o ' l l l l - l
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
8
—m'
8.,

Fig. 24 Effective constitutive laws for the case of 2'72“: 0.5, we" = 0.1

136

 

Fig. 25 Damage pattern for random arrangement of inclusion,

292'": 0.5, c‘rc'" = 0.1

137

 

Fig. 26 Damage pattern for random arrangement of inclusion,

s‘/e’"= 0.5, C‘rc’" = 0.1

138

 

Fig. 27 Damage pattern for random arrangement of inclusion,

e’/e"'= 0.5, Ci/C'" = 0.1

139

 

Fig. 28 Damage pattern for random arrangement of inclusion,

5'72”“: 2, C‘/C’" = 0.2

 

Fig. 29 Damage pattern for random arrangement of inclusion,

e‘/e"'= 2, d/C’" = 0.2

141

 

Fig. 30 Damage pattern for random arrangement of inclusion,

eye": 2, C‘/C’" = 0.2

142

 

Fig. 31 Damage pattern for random arrangement of inclusion,

sin-1m: 2, CJ/Cm = 0.2

143

Q

 

”I

 

Q

or

0.0 r 1 1 1 1 .J 1
I
I

0.5

I

0 .4
matrix.

I

I

0.2

0.1

I
5 o
a. c
1 B.
c. o
:1

 

 

 

0
_ol'
I
I.
I.
,-
0..
0
."

 

0 0.1 02 0.3 0.4

 

Fig. 32 Effective constitutive laws for the case of

e‘/e’"= 2, Ci/C'" = 0.2

 

Fig. 33 Damage pattem for random arrangement of inclusion,

eye: 10, C'VC'" = 0.1

145

 

 

0‘ 1 H 011011 1 1 1 1 1
,1
0.35 '" In" ‘ '21" | .1
.Il Id
.’ |
o 3 1- 1/ 1’ i ..
' I ' /,_ .
11‘, 1.1 | /
i v1 /| / /|/ : .4: /|
O i' ' 1 l / - . .’ I-’| V I l
7 25 I. (A 3.4., 1 O / 54:" 1’ I [’1 1 '1
CC, .1 : / / I, .- v‘ , l
r ’ | :
°2 ! 1 «0111)
I
0.15 " '

 

 

 

 

 

l
01" 1 7
incl . 1 W
i I
015- 1 : -
l
o l J l 1 l l I | 1
0 02 0.4 0.6 00 1 12 1.4 1.6
8
'71
8.,

1.9

Fig. 34 Effective constitutive laws for the case of ei/em= 10, Cl/Cm = 0.1

146

 

Fig. 35 Damage pattern for random arrangement of inclusion,

2'72“: 10, Ci/C'" = 0.1

147

 

Fig. 36 Damage pattern for random arrangement of inclusion

e‘/e’"= 10, C‘/C’"=o.1 ,L=127

148

C1/C'"

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

m m m
o/ocr G/GC, o/ocr
1 101 r“ r 1 ‘00.._.-__ ”'W’W“
. 1 1 1
. : ; 1 1
. ‘ L 1
10 .5. 1; C * m 51 1 50 1
’ e w 1 1
'1 1 1 _ (1'9”? m """ m C
o 5 1 0 .5 1 0 5 10
m m m
E/ECI‘ E/ECV E/ECV
m m m
o/ocr o/ocr o/ocr
1,; 1 l I, 1
1 m .‘/ - 4
l l i i , ' ~ 1
1 1 1 1 'm . 1
1 /__, £— — ~——' 1 . C
0' 5 1 o 5 1 o 5 IO
m m m
8/8
C, e/ecr e/ecr
m "I m
o/ocr o/ocr o/ocr
1. 1’[ 1 1.1 TIn
m 1 m * 1 f§ i _-
‘. 1 g: '
5 5 .5 '1. ,
0.1 J 1‘ :
r 1‘ ,1 C
C l C 1 1,"11’11 '.
l ”rm/’1 1 1’14“" '1 1
.,~._/”_~ i L." ~ a . , " a- _ - - _. _
o 5 1 0 .5 1 0' 5 10
e/em "1 "'
cr E/ec, E/Ecr
0.1 1 10 1' / m
8cr 8cr

Fig 37 Damage map of effective constitutive laws, also showing

response of inclusions and matrix

o/ocr

l49

o/ocr

Cl/Cm

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

r
m C
11 u E 4 11J
A .- /
x . 8
/. u
/. u //
/. o //
/ ..
., /
// l ,/
/ . //
/ r
/ m C
/ . O
. m H mm m ".1”an W -1 111111..
.1
m C
. 1 1 u a 1 11
a.
/ . ac /.
/ .m / a
/ , .
./., fl D
/.r A."
/, a
a" ..
.. .6 /
.. ,,
// .
Jh r
, . mCC
F 1 L 1 L 11 I. 1 1111
“uunuuuuuw. /.""..uunnnr.
G
r
. m C
a. r E
xv /
E 11-1 -i-
,x
u r ,
/ m C
o -- ,
_u.uu..ua /.r..pn....m-

 

 

 

 

e/ecr
e/ecr
i / m
r ECI'

8C

10

OCT

 

 

O/

 

 

e/ccr
e/ccr

"1
cr

 

 

o/o

m
cr

e/e,

e/ 8
Fig. 38 Damage map of effective constitutive laws

0.]

OCI'

 

 

I.-
III' 1'

G/

 

0.1

ISO

0.]

 

 

0.1

1
i m
Scr/Scr

Fig. 39 Damage map of damage patterns

151

 

Fig. 40 Effect of fluctuation of inclusion arrangement on damage pattern,

s'Ve’”: 2.0, c1/c’" = 0.2

152

o/kn,

 

0.6 I I I I I I I r

0.5- 4

03*

0.1r

 

 

 

 

m

e/eu

Fig. 4] Effect of fluctuation of inclusion arrangement on stress-strain laws,

21/2’": 2.0 C1/C'" = 0.2

153

 

Fig. 42 Effect of fluctuation of inclusion arrangement on damage pattern,

6713’": 2, Ci/C'" = 5

154

 

Fig. 43 Effect of fluctuation of inclusion arrangement on damage pattern,

31/2“: 0.5, C1/c’" = 0.1

155

 

Fig. 44 Effect of fluctuation of matrix material on damage pattern,

2'72“: 2, C1/c’" = 0.2

 

 

 

 

 

"I
o/ 0.,
0.6 I I I I I I I T
uniform matrix
0.5 — ..
. . l
_ - _ - fluctuated matrix material .
I
l
0.4 ~ 1 -
I
I
l
l
0.3 - , .
l
l
|
|
0 2 ” | "
|
l
I
l
0.1 r- l 1
l
I
I
o 1 l 1 l l l l 1
o 01 02 03 on as as 07 no an
m
e/ 8.,

Fig. 45 Effect of fluctuation of matrix material on constitutive law,

2'72”: 2, C1/c’" = 0.2

157

 

Fig. 46 Effect of fluctuation of matrix material on damage pattern,

151/81": 0.5, C1/C'" = 0.1

158

 

 

0/63
035 1
0.3

0.25

0.2

0.15

0.1

0.05

Fig. 47 Effect of fluctuation of matrix material on constitutive law,

 

_ _ _ _ fluctuated matrix material

uniform matrix

 

 

 

0.1 0.2 0 8 0.4 05 0.6 0.7
m
e/a,

i m
8/2

0.5, C1/C’" = 0.1

159

 

Fig. 48 Effect of fluctuation of matrix material on damage pattern,

(ii/em: 2, Ci/C'"=5

o/o”

00

0]

06

05

0!

03

02

01

0
0

"I

160

 

 

uniform matrix

 

_ _ _ - fluctuated matrix material

1 L L l l

 

 

 

005 01 015 02 025

03

Fig. 49 Effect of fluctuation of matrix material on constitutive law,

31/2“: 2, Ci/c’" = 5

o/oc,

161

 

0.35

I

0.25 P

I

0.2

0.15"

0.1

0.05 -

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.2 0.4

0.6

Fig. 50 Effective constitutive curve,

and 20 samples

8

i/Em

éw

8/8

1, Ci/C'" = 0.1

1.2

m
Cf

162

 

 

 

 

 

F(omm/o:1r)
1 . . . . . 1 1
0.0“ _
0.3- a
0.7- a
0.6- -
0.5- ..
0.4“ -1
0.3~ -
0.2% -
0.1- .

 

 

 

0
0.22 0.24 0.26 0.26 0.3 0.32 0.34 0% 0.30

1 m
Ornax/Gcr

Fig. 51 Probability distribution function F (6mm) for the case of

9.1/3“: 1, C1/Cm = 0.1

/8

max C r

163

 

0.0"

0.01

0.7-

0.6-

0.3~

0.1-

 

 

 

 

0.6

0.7 0.3

0.0

Fig. 52 Probability distribution function F (8mm) for the case of

i m
e/e

1, Ci/c’" = 0.1

164

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

o/o:

0.5 . . . . . T . .

1145- 4
0.4— -
1135- 1 -
0.3r 11‘1 1 ‘
o.25~ ~1 .
0.2- / -

. V/

0.15- " -
11.1 - -
005- D/I'l -

Fig.53 Effective constitutive curve, 81/€m= 0.1, C1/ Cm = 10

and 20 samples

0.0

0.3

0.7

0.6

0.5

0.4

0.3

0.2

0.1

165

 

I

 

 

 

 

 

 

 

F. .

l l l

 

 

 

 

 

l

l

l

 

l

 

0
0.33 0.30 0.

4

0.41 0.42 0.43

0.44

0.45

0.46

0.47 0.43

m
Gmax/Gcr

Fig. 54 Probability distribution function F (6mm) for the case of

i m
8/8

0.1, Ci/C'" =10

/8

max Cf

l

166

 

0.0"

0.7-

0.6"

0.4”

0.3»-

 

 

l

 

 

0
0.2

0.25 0.3 0.35

0.4

0.45

0.5

m
emax/Scr

Fig. 55 Probability distribution function F( 8mm) for the case of

2’72”: 0.1, Ci/C’" = 10

167

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/02
0.7 l l I l l
0.6“ .
0.5 ' .4
0.4 ~ ”/1/ -
u
0.3“ J '
r"
0.2 ' ‘
W./ //
l .55”
01 _ w u ”JAM ‘
.i-PFFFFFFV]
0 l L 1 J 7 1
0 0.1 - 02 0.3 ' 0.4 0.5 0.6
e/eZ

Fig. 56 Effective constitutive curve,

and 20 Samples

9372’": 10, C’/C’"= 10

168

F(o mm)

max Cf

 

0.0“

0.8-

0.6"

05"

0.3-

0.2’

 

 

 

 

0.4 0.45 0.5 0.55 0.6 0.65

m
GmaX/GCI'

Fig. 57 Probability Distribution function F<Gmax) for the case of

2’72“: 10, C’/C’"= 10

169

m
max/EU)

 

00*

0.7~

0.5"

0.4“

0.2”

0.1"

 

 

 

 

o J l l l l l l

0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6

m
Emax/Ecr

Fig. 58 Probability distribution function F( Ema): ) for the case of

2’72“: 10, (3’76" = 10

170

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

m
o/oc,
0.0 l I I l I l T
0.7- .
0.6- _
A
05- L -
04- .
.4
,4
0.3” 4
0.2- ‘ -
L4
3L
o . . . . Fl .1 m ‘. I
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.0
m
e/ec,

Fig. 59 Effective constitutive curve, €i/€m= 1, Ci/Cm =10

and 20 samples

171

 

0.3-

0.6"-

0.5-

0.3-

0.2-

 

 

1 1 J l l 1

0
0.5 0.55 0.6 0.65 0.7 0.75 0.3 0.05

 

 

m
Gmax/ OCT

Fig. 60 Probability distribution function Rom“) for the case of

2'73””: 1, C’/C’"= 10

172

F(e a“)

 

 

 

 

o I l l l l l l l

0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.0 0.05

 

m
emax/Ecr

Fig. 61 Probability distribution function F( 8mm ) for the case of

2'72“: 1, C‘/C’" = 10

 

173

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0/02
as t r I 1 I I . 1' .
0.45 - a
0.4 ~ -
0.35 ~ .
0.3 - .
0.25 - .
0.2 ~ fl .
,.
0.15 - «WE tx”J/M .
0.1 r ~
0.05 - 4
o i . i i . i i . .

 

Fig.6.? Effective constitutive curve, ei/em= 10, Ci/Cm = 0.1

and 20 samples

174

 

 

0.0-

0.7-

0.6-

T

0.5

04'-

0.3-

0.11-

 

 

 

° 1 J 1 1 1 1 1 l

0.41 0.42 0.43 0.44 0.45 0.46 0.47 0.40 0.40 05

 

(Sm/02".

Fig. 63 Probability distribution function F(omax) for the case of

e’/e’"= 10, Ci/c’" = 0.1

175

 

 

 

 

 

 

 

.4

1.0 1 9
m
8max/8N

Fig. 64 Probability distribution function F( em“) for the case of

e‘/e’"= 10, 676’" = 0.1

176

<1/o"
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1

0.05

 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.0 I

m

e/ea

Fig. 65 Effective constitutive curve, (ii/em: 0.1, Cl/C’" = l

and 20 samples

177

 

 

F( (Sm/0:)

1 T i . I I . I

I
0.0- i
0.0- .
0.7- 4
0.6- a
0.5» .
0.4- .
o.a~ .
o.2~ .
0.1- .

i i i . i . .

 

 

 

0
0.33 0.34 0.35 0.06 0.37 0.30 0.30 0.4 0.41

Gmax/Gcr

Fig. 66 Probability distribution function F(omax) for the case of

673'": 0.1, C‘/c’" = 1

178

 

 

 

 

 

 

F( Ema/8:1,)
1 I I I I I I
0.0 ~ ~
0.0 - 1
0.7 - .
0.6 - .
0.5 - a
0.4 - .
0.0 - .
0.2 r .1
0.1 - a
3.05 077 033 070 0.05 070 0.05 i 1.05
amt/8:2,

Fig. 67 Probability distribution function F( 8mm) for the case of

Eli/Em: 0.1, Ci/C'"=1

179

 

 

 

m
o/oc,
0.25 I I I I I T I T l
0.2 - -
0.15 - -
1
0.1 F .
0.05 - .

 

 

 

 

 

 

 

 

0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.5 0.4 0.45 0.5

m

e/cc,

Fig. 68 Effective constitutive curve, e’/e’"= 0.1, C’/c’" = 0.1

and 20 samples

F(o Kym)

max cr

00

00

07

06

05

04

00

02

01

0
015 01

180

 

 

l

1 1

l

l

 

4 l

 

Fig. 69 Probability distribution function F(on

6 017

i m
3/8

010 040

0.1, (“76" = 0.1

02

021

tax

022 023 024

Gmax/Gcr

) for the case of

181

 

 

 

 

 

F1 tam/8:1)

' ' ' ' ' ' ' 1
0.0 - '
0.0 ~ '
0.7 - .
0.6 - _
0.5 - -
0.4 - .
0.3 ~ .
02 ~ -
0.1 ~ -

3.32 0.04 0.36 0.30 0.14 0.42 0.44 0.46

m
Emax/Scr

Fig. 70 Probability distribution function F( em“) for the case of

e‘/e’"= 0.1, Ci/c’" = 0.1

182

‘35.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

C /
m m m
o/ocr o/ocr o/ocr
.50~*-- ~— MW“-J_ .80' - . - .70
. ,a’ .
I /'
10 j ,2
1t /
1 .
0 fl ’1 ‘ 0 0 1 JV 1
0 .45 O 0 .60
m m
e/ecr e/ecr
m m m
o/ocr o/ocr o/ocr
'01 1.0 /
/
1 . /
.//
I 1 , ; /
[ / 1 ~ / 1
t / i /
0 'L,.’ 0
0 1.0 0 1.0
m m
e/ecr e/ecr
m m m
o/ocr o/ocr o/ocr
0.1

 

 

 

 

 

 

 

 

 

0.1

H

Fig. 71 Damage map, constitutive curves

 

 

 

 

 

 

 

 

 

183

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

C’ /C’"
m m
F [Om/0.. F[loma_or/ F [Om/0 2",]
‘F __ _JH .. .
H . JJIJ ff
[_
10 ' f ' .
’ .1
: JH
o: of ___s o J
. 0.5 .85 0.40 "1.65
m m
G/Ocr G/Ocr G/Ocr
m
F [Gina/<5 I F [Gmax/ mCi’] F [Gmax/Gc
If ’F r:..'.__. ‘ l” 1’ fl "*1”?
1 1r“ '1 f ’ ‘
.r‘ . .
1 i if I ;
1 j i :
0 'rr 0: 01 l
.33 .41 O 1.2 0 ”('12
m m
o/ocr o/ocr G/GC,
m m m
F [omax/Gc F [Gmax/Gc F [Gmax/Oc
1 '_+ ______I lrj__1 l r___.'
1" f" ,._,l
J Ir ,
0.1 » I“ r] a4
_ _J y I I
1 J _r'
.1 F Jr J .
Of - 0’ 1" 0 I
.16 .24 0.22 .36 0.41 .50
m m - "1
o/ocr o/oc, G/Ocr
i m
0.1 1 10
8cr/fic‘r

Fig.72 Damage map, distribution function F(omax )

184

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

i m
C / C
m m
F[£max/8:‘nr] F[€max/£Cr:| F[£max/€Cr]
l _I—F—J—1 1 1‘ ___H—J I» :J ‘
; 1 ; 4H; 3 ’ HJ
10 1 1 . 1 r5 i . 1
0’ rr‘ 1 o T— 4 0 f!— T
20 m so .40 m 85 20 moo
e/ecr e/ecr e/ecr
"1 m m
F[Emax/8cr F[Emax/Ecr F[8maX/€Cr
l v 1 . l P f
1 J i 1 1
1 -—J" - 1 1
r H .1 ’
1 t
{.1 . 1 1 I 1 4
o . 0' L o '
.65 ,,} .05 o 1.2 0 ,,.‘ ~2
E/Ecr e/EZ: E/Ecr
m m
F[€max/€Cr:l F [ENIGX/ECI’] F[8max/8:‘nr]
1 'I .———-————}1 l ; _ ___ 1 > J
1 “7 ’ 1 f . . I” .
1 1 I 3L 1— 3 I —-H i
0.1 1 f 1 1 {I i . H]
l j- | 4
r‘ 11 irr 1 ’ Ill 1
o J} 01 _ -.,,.-_-__ -__.____ 1
.32 ”1.46 .5 m 1.1 1 2 ”1'9
e/ccr e/ecr s/ecr
8cr/ecr

Fig. 73 Damage map, distribution function F( 8

max )

 

'—
.I