Jtam.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

3, 39, 2001

STOCHASTIC RELIABILITY IN CONTACT PROBLEMS FOR

COMPOSITES WITH SPHERICAL PARTICLES

Marcin Kamiński

Division of Mechanics of Materials, Technical University of Łódź

e-mail: marcin@kmm-lx.p.lodz.pl

The application of the stochastic perturbation technique in the classical
Hertz contact problemof a compositewith spherical particle is presented
in the paper. The perturbation approach is used in conjunctionwith the
Weibull-Second Order Third Moment (W-SOTM) reliability model in
numerical analysis of the composite. All computational experiments are
carried out using symbolic computation packageMAPLE, thatmakes it
possible to visualize the deterministic sensitivity gradients as well as to
show probabilisticmoments of contact stresses. The entire methodology
can be applied to analytical solutions of contact problems with a more
complicated contact surface geometry and, on the other hand, to com-
putational implementation in the Finite or Boundary Element Method
programs for simulations of some contact phenomena.

Key words: reliability analysis, contact problems, composite materials,
stochastic perturbationmethod

1. Introduction

Usually, contact problems are formulated for homogeneous bodies, but in
general, they are very important for composites as well. Contact stresses can
be analyzed in the case of composites reinforced with spherical particles and
particles with a more general shape such as ellipsoidal or paraboloidal one,
for instance. Moreover, contact problems in composites can be analyzed to
simulate delamination effects or interface defects between various components
in the laminates andfiber-reinforcedcomposites.As it canbeexpected, contact
stresses dependon interrelations betweenmaterial properties of the contacting
composite components. The analysis is definitelymore complicated in the case
of a randomly definedmaterial and/or geometrical parameters – not only the



540 M.Kamiński

expected values, but also higher order probabilistic moments can be decisive
for the final results of the contact problem.

At the same time, significant progress is being observed in the field of sto-
chastic mechanics (see Grigoriu, 2000). Probabilistic numerical experiments
with Monte-Carlo simulation (MCS) (see Kamiński, 2001; Kleiber and Hien,
1992) and with Stochastic Finite Element Method (SFEM) (see Kleiber and
Hien, 1992; Liu et al., 1986) prove that usually the second order perturba-
tion analysis is exact enough, while the second probabilisticmoment approach
is not satisfactory, considering the fact that in most cases the third order
(skewness) parameters differ from 0. Additionally, the stochastic perturba-
tion method gives essential computation time savings in comparison with the
MCStechnique and the stochastic spectralmodeling (seeGhanemandSpanos,
1992, 1997). Furthermore, even if a Gaussian random variable is very precise-
ly digitally simulated, the third order probabilistic moments and coefficients
slightly differ from 0. That is why, the SecondOrder ThirdMoment (SOTM)
stochastic analysis is proposed below, in addition to the reliability study of
contact problems based on the Weibull probability density function (PDF).
It should be underlined that the Weibull probability density function is fre-
quently used in probabilistic analysis of fracture phenomena in homogeneous
structures and composites (see Beyerlein and Phoenix, 1997; Lund, 1998).

On the other hand, it is known that the reliability analysis has beenmode-
led using the First Order Reliability Method (FORM) and the Second Order
Reliability Method (SORM) (see Brandt et al., 1984; Ghanem and Spanos,
1992). These approaches do not enable one to include the influence of the
third order probabilistic characteristics, necessary for the Weibull statistics,
and that is why the Third Order Reliability Method (TORM) is proposed on
the basis of the SOTM stochastic perturbation method.

To illustrate the entire approach, the stochastic perturbation reliability
analysis of a linear elastic contact is carried out for a composite reinforced
with spherical particles. Since the solution for the deterministic problem is
knownandhas beenworked out analytically, the probabilistic analysis ismade
using symbolic computations package MAPLE (see Char et al., 1992). Com-
putations of the reliability limit function and probabilistic moments of the
contact stresses as well as some numerical sensitivity studies together with
visualization of all computed functions are carried out by the use of this pac-
kage. This methodology can be successfully applied to randomization of all
contact problems in reliability studies, where contact stresses are described
by the closed form equations. Otherwise, computational implementations of
theStochasticFinite (seeBuczkowski andKleiber, 1999) orBoundaryElement



Stochastic reliability in contact problems... 541

Method (seeBurczyński, 1995) are to bemade in order to get general approxi-
mative probabilistic solutions to contact problemsof composites. Furthermore,
theW-SOTMnumerical approach to the stochastic reliability, stochastic con-
tact modeling and related aspect of analytical computations can be applied
and explored in various areas of modern engineering, especially in the field of
composite materials (see Fish, 1999).

2. Contact problem of a particle reinforced composite

Let us consider a contact phenomenon between two linear elastic isotropic
regions characterized by the Young’s moduli E1, E2 and Poisson’s ratios
ν1, ν2. Let us assume that the regions have spherical shapes with radii R1
and R2, respectively, and that the contact is considered at the point C, as
it is shown in Fig.1. A 3D view of the particle-reinforced composite plane
cross-section is shown in Fig.2.

Fig. 1. Geometry of contact surface

Let us observe that the contact problem is axisymmetric with respect to
the vertical axis introduced at the center of a spherical particle and at the
bottom of this sphere (Fig.1). It is assumed that there is no pressure between
the composite constituents and, therefore, the contact appears only at the
point C. The distance between the points on the contacting surfaces and the
plane transverse to the vertical axis of both surfaces is assumed to be small
and can be described, according to Timoshenko and Goodier (1951), as

z2−z1 =
r2(R1−R2)

2R1R2
(2.1)



542 M.Kamiński

Fig. 2. 3D cross-section view of particle-reinforced composite

where r denotes the distance between the points M, N and the symmetry
axis introduced at C. If the composite components are loaded by the vertical
force P acting along the vertical axis at the point C, then some local strains
are induced in the neighborhood of this point. They are a result of the contact
phenomenon on a small circular surface (contact area). Assuming that the
radii R1 and R2 of the composite constituents are sufficiently greater than
the radius of the contact area, then the results of the Bussinesq problem of
a linear elastic halfspace loaded by a concentrated force can be adopted here
(see Timoshenko and Goodier, 1951). For this purpose let us denote by w1
the vertical displacement induced by the local vertical strain of the point M
belonging to the matrix; and by w2 – the corresponding displacement of the
point N in the vertical direction. Finally, assuming that the tangential plane
at the point C remains unmovable during local compression, the close-up of
the two points M and N can be expressed by some real η as

η=α− (w1+w2)=
r2(R1−R2)

2R1R2
(2.2)

If M and N belong to the contact area, their displacements wi for i=1,2
can be written as (see Timoshenko and Goodier, 1951)

wi =
1−ν2i
πEi

∫∫

q dsdϕ (2.3)

which follows the symmetry of the pressure intensity q and the corresponding
local strains with respect to the vertical axis at the point C. The integration
is carried out over the entire contact surface, therefore

(k1+k2)

∫∫

q dsdϕ=α−
r2(R1−R2)

2R1R2
(2.4)



Stochastic reliability in contact problems... 543

where k1 and k2 can be obtained from Eq. (2.3). Now, we are looking for
such an expression for q, so that the above equation could be fulfilled. It can
be obtained for the pressure distribution on the contact surface represented
by the coordinates of a hemisphere with the radius a constructed on that
surface. If q0 is taken as the pressure at the point C, then one can show that

∫

q ds=
qo
a
A (2.5)

where

A=
π

a
(a2−r2 sin2ϕ) (2.6)

what gives

π(k1+k2)q0
a

π

2
∫

0

(a2−r2 sin2ϕ) dϕ=α−
r2(R1−R2)

2R1R2
(2.7)

Finally, the parameters a and α can be found as

a=
3

√

3π

4

P(k1+k2)R1R2
R1−R2

(2.8)

α=
3

√

9π2

16

P2(k1+k2)
2(R1−R2)

R1R2

which gives the maximal pressure on the contact surface

q=
3P

2πa2
(2.9)

Then, the normal stress can be defined as

σz =−

a
∫

0

3qr dr z3(r2+z2)
−5

2 = q
[

−1+z3(a2+z2)
3

2

]

(2.10)

Let us note that the shear stresses are equal to 0, which results from the
spherical symmetryof the reinforcingparticle, however in thecase of ellipsoidal
reinforcement the shear stresses differ from 0 (see Timoshenko and Goodier,
1951).

Themain purpose of further analysis is to determine probabilistic charac-
teristics of themaximal contact stresses as well as contact surface geometrical



544 M.Kamiński

parameters. Since the spherical particle surrounding thematrix is considered,
let us assume that the difference R1−R2 is smaller than R2.

This parameter is treated as the input design parameter in further sen-
sitivity analysis. The characteristics mentioned above are necessary in final
stochastic reliability computations and, considering complexity of equations
describing the reliability parameters, the second order stochastic perturbation
method is proposed.

3. Stochastic perturbation approach

Let H be some Hilbert’s space defined on the domain D with the real
values set R. Let (Ω,σ,P) denote the probability space, x∈D and θ ∈ Ω,
while Θ : Ω → R. Let us consider the differential operator L(x;θ) defined
on the product H×Θ, whose coefficients exhibit random fluctuations with
respect to some independent variables. The random coefficients of the diffe-
rential operator are restricted to be second order random processes. Then,
we illustrate the stochastic second order third central probabilistic moment
approach using a solution to the following equation

L(x;ω)u(x;ω)= f(x;ω) (3.1)

being amathematical description of a physical problem.The randomresponse
of the system is denoted here by u(x;ω), while its admissible excitation by
f(x;ω).

As it is known, analytical solutions to such a class of partial differential
equations are available for some specific cases only, and that is why a general
methodof extending such solutionswouldbeverypowerful.Variousmathema-
tical methods can be proposed here, some of them are listed by Ghanem and
Spanos (1997), Kleiber and Hien (1992), however, in this paper, the second
order perturbation third central probabilisticmoment approach is applied (see
Peng et al., 1998). The second order perturbation follows a traditional appro-
ach in this area (and the lack of convergence studieswith respect to theTaylor
expansion order). The third probabilistic momentmethod reflects the need of
modeling of unsymmetric random variables.

Let us denote for this purpose the vector of input random fields by
{br(x;ω)} and its probability density by p(br) and p(br,bs) respectively,
where r,s= 1,2, ...,R are indices of the random variables. Then, the expec-



Stochastic reliability in contact problems... 545

ted value of br is defined as

E[br] =

+∞
∫

−∞

brp(br) dbr (3.2)

and the covariance as

Cov(br,bs)=

+∞
∫

−∞

+∞
∫

−∞

(br−E[br])(bs−E[bs])p(br,bs) dbrdbs (3.3)

Next, the skewness parameter S(ui) is calculated

S(u)=
1

σ3(u)

+∞
∫

−∞

(u−E[u])3pR(b) db (3.4)

In further applications, the Weibull distribution is used with the probability
density function defined as

pR =











β

λ

(x−x

λ

)β−1
exp
[

−
(x−x

λ

)β]

for x>x

0 for x<x

(3.5)

where
β – Weibull shape parameter
λ – scale parameter
x – location parameter, which indicates the smallest value of the

random variable x for which the probability density function
is positive.

Considering this definition, theWeibull PDF is used for general mechani-
cal applications, wheremany randomvariables must be nonnegative (Young’s
modulus, some geometrical parameters, for instance) and especially in com-
posites failure and fatigue modeling (see Beyerlein and Phoenix, 1997). Let
us note that if a discrete representation of the random variable b(x;θ) is
used, then statistical estimators can be applied to approximate probabilistic
moments of this variable of any order.
Afterwards, all material and physical parameters of Ω and all state func-

tions being random fields are represented by the use of the stochastic second
order expansion via the Taylor series as follows

F
(

b(x;θ);x
)

=F0
(

b0(x;θ);x
)

+εF ,rs
(

b(x;θ);x
)

∆br+
ε2

2
F ,rs
(

b(x;θ);x
)

∆br∆bs

(3.6)



546 M.Kamiński

where ε is a given small perturbation, ε∆br denotes the first order varia-
tion of ∆br about its expected value E[br] and F(n)

(

b(x;θ);x
)

represents
the nth order partial derivatives with respect to the random variables, while
F(n)
(

b0(x;θ);x
)

denotes the nth order partial derivative calculated at the
expected value of the input random variables vector. Let us note that the se-
cond order equation is obtained here bymultiplying the R-variate probability
density function pR(b1,b2, ...,bR) by the ε

2-terms and integrating over the
random field domain b(x;θ).
Applying the stochastic second order Taylor expansion to Eq. (3.1) and

equating terms of the same order for ε=1 yields:
— 0th order equations

L0
(

b0(x;θ);x
)

u0
(

b0(x;θ);x
)

= f0
(

b0(x;θ);x
)

(3.7)

— 1st order equations (for r=1, ...,R)

L0
(

b0(x;θ);x
)

u,r
(

b0(x;θ);x
)

= f ,r
(

b0(x;θ);x
)

−L,r
(

b0(x;θ);x
)

u0
(

b0(x;θ);x
)

(3.8)
— 2nd order equations (for r,s=1, ...,R)

L0
(

b0(x;θ);x
)

u(2)
(

b0(x;θ);x
)

=
{

f ,rs
(

b0(x;θ);x
)

−
(3.9)

−L,rs
(

b0(x;θ);x
)

u0
(

b0(x;θ);x
)

+2L,r
(

b0(x;θ);x
)

u,s
(

b0(x;θ);x
)

}

Cov(br,bs)

Equation (3.7) is solved for u0 to obtain the probabilistic solution for the
considered equilibrium problem, they Eq. (3.8) – for first order terms of u,r

and, finally, Eq. (3.9) – for u(2). The two probabilistic moment characteri-
zation of all the state functions for the boundary value problem starts from
the expected value of the stress tensor components. Using definition (3.2) and
substituting the second order expansion it can be found (seeKleiber andHien,
1992), that

E
[

u[b(x;θ);x]
]

=u0[b(x;θ);x]+
1

2
u,rs[b(x;θ);x]S

,rs
b

(3.10)

Next, the first order cross-covariances for the solution components are de-
rived as follows (see Kamiński, 2001; Kleiber and Hien, 1992)

Cov
(

u[b(x(1);θ);x(1)];u[b(x(2);θ);x(2)]
)

=Su(x
(1);x(2))=

=

+∞
∫

−∞

{

u[b(x(1);θ);x(1)]−E
[

u[b(x(1);θ);x(1)]
]

}

· (3.11)

·
{

u[b(x(2);θ);x(2)]−E
[

u[b(x(2);θ);x(2)]
]

}

pR
(

b(x;θ)
)

db



Stochastic reliability in contact problems... 547

which gives

Su(x
(1);x(2);θ)=u,r(x(1);θ)u,s(x(2);θ)Srsb (x

(1);x(2);θ) (3.12)

Analogously, the third order probabilistic moments can be derived in the sto-
chastic second order perturbation approach.

4. Weibull-Second Order Third Moment reliability approach

Starting from the moments and the stochastic perturbation methodology
presented in the previous section, we compute the first three probabilistic cha-
racteristics of the vertical stresses E[σz(x;ω)], var

(

σz(x;ω)
)

and S
(

σz(x;ω)
)

as follows (see Timoshenko and Goodier, 1951)

E[σz] = σ
0
z +
1

2

n
∑

i=1

∂2σz
∂b2i
σ2(bi)

(4.1)

σ2(σz) = σ
2
z +

n
∑

i=1

[(∂σz
∂bi

)2
+σz
∂2σz
∂b2i

]

σ2(bi)+

+
n
∑

i=1

∂σz
∂bi

∂2σz
∂b2i
S(bi)σ

3(bi)−E
2[σz]

and

S(σz) =
{

σ3z +
3

2

n
∑

i=1

[

2σz
(∂σz
∂bi

)2
+σ2z
∂2σz
∂b2i

]

σ2(bi)+

+
n
∑

i=1

[(∂σz
∂bi

)3
+3σz

∂σz
∂bi

∂2σz
∂b2i

]

S(bi)σ
3(bi)− (4.2)

− E3[σz]−3E[σz]σ
2(σz)

} 1

σ3(σz)

Having computed the first three probabilistic moments of contact stresses
(expected values, standard deviations and skewness coefficients), the random
field of the limit function g(z;ω) is to be proposed. Usually, it can be in-
troduced as the difference between the allowable and actual stresses σz(z;ω)
induced in the composite

g(z;ω) =σall(ω)−σz(z;ω) (4.3)



548 M.Kamiński

Let us underline that the allowable stresses are most frequently analyzed as
random variables in the interior of statistically homogeneous materials, whe-
reas the actual stresses are random fields. That is why the computational
analysis is carried out in Section 5.3 for a specific value of the vertical coordi-
nate z. The random variable of the allowable stresses σall(ω) is specified by
the use of the first three probabilistic moments E[σall(ω)], var(σall(ω)) and
S(σall(ω)). Then, the corresponding probabilistic characteristics of the limit
function are calculated as

E[g] = g0+
1

2

n
∑

i=1

∂2g

∂b2i
σ2(bi)

(4.4)

σ2(g) = (g0)2+
n
∑

i=1

[(∂g

∂bi

)2
+g0
∂2g

∂b2i

]

σ2(bi)+

+
n
∑

i=1

∂g

∂bi

∂2g

∂b2i
S(bi)σ

3(bi)−E
2[g]

and

S(g) =
{

(g0)3+
3

2

n
∑

i=1

[

2g0
(∂g

∂bi

)2
+(g0)2

∂2g

∂b2i

]

σ2(bi)+

+
n
∑

i=1

[(∂g

∂bi

)3
+3g0

∂g

∂bi

∂2g

∂b2i

]

S(bi)σ
3(bi)− (4.5)

− E3[g]−3E[g]σ2(g)
} 1

σ3(g)

By inserting the limit state function g from Eq. (4.3) into Eqs (4.4) and
(4.5) and assuming that the randomvariable of the allowable stresses and the
random field of the actual stresses are uncorrelated, one obtaines

E[g] = σ0all−σ
0
z −
1

2

n
∑

i=1

∂2σz
∂b2i
σ2(bi)

(4.6)

σ2(g) = (σ0all−σ
0
z)
2+

n
∑

i=1

[(∂σz
∂bi

)2
− (σ0all−σ

0
z)
∂2σz
∂b2i

]

σ2(bi)+

+
n
∑

i=1

∂σz
∂bi

∂2σz
∂b2i
S(bi)σ

3(bi)−E
2[σall−σz]

and



Stochastic reliability in contact problems... 549

S(g) =
{

(σ0all−σ
0
z)
3+
3

2

n
∑

i=1

[

2(σ0allσ
0
z)
(∂σz
∂bi

)2
+(σ0all−σ

0
z)
2∂
2σz
∂b2i

]

σ2(bi)+

+
n
∑

i=1

[(∂g

∂bi

)3
+3g0

∂g

∂bi

∂2g

∂b2i

]

S(bi)σ
3(bi)− (4.7)

− E3[σ0all−σ
0
z]−3E[σ

0
all −σ

0
z]σ
2(σ0all−σ

0
z)
} 1

σ3(σ0
all
−σ0z)

Comparing the second order second moment (SOSM) approach with the
second order third moment (SOTM) approach, it is seen that the expected
values are described by exactly the same equation, while the standard devia-
tions (or variances) have some extra components resulting from the skewness
of the analyzed PDF; the third order parameter of the output PDF is taken
into account in the SOTM-based analysis, only.
Finally, the parameters x(g), λ(g) and β(g) of the equivalent Weibull

PDF are derived from the following system of equations

E[g] = λ(g)Γ
(

1+
1

β(g)

)

+x(g)

σ2(g) = λ2(g)
[

Γ
(

1+
2

β(g)

)

−Γ2
(

1+
1

β(g)

)]

(4.8)

S(g) = λ3(g)
[

Γ
(

1+
3

β(g)

)

−3Γ
(

1+
2

β(g)

)

Γ
(

1+
1

β(g)

)] 1

σ2(g)

+ 2Γ3
(

1+
1

β(g)

) 1

σ2(g)

where the Gamma function is defined as

Γ(x)=























∞
∫

0

e−ttx−1 dt for x> 0

lim
n→∞

n!nx−1

x(x+1)(x+2)...(x+n−1)
for x∈R

(4.9)

It should bementioned that approach based on the symbolic computations is
themost efficientmethodof solving these equations andof computingprobabi-
listic parameters of the equivalentWeibull distribution. Finally, the structural
reliability index R of the limit function g is calculated from the following
formula

R(g)= exp
[

−
(

−
x(g)

β(g)

)λ(g)]

(4.10)

Values of this index should behave like the classical probability index, i.e. be
greater or equal to 0 and less or equal to 1.



550 M.Kamiński

5. Computational experiments

Computational experiments are conducted by the use of system MAPLE
for symbolic computations (seeChar et al., 1992), where the stochastic second
order perturbation method in W-SOTM reliability analysis of contact pro-
blems has been implemented. The entire analysis is divided into three groups
of essentially various numerical examples:

1. Deterministic analysis and sensitivity study (see Haug et al., 1986) of a
contact problemwith respect to the vertical spatial coordinate

2. Numericalmodeling by randomizingmost of the input parameters based
on the stochastic second order perturbation

3. Stochastic reliability modeling according to the Weibull Second Order
ThirdMoment approach (see Peng et al., 1998).

5.1. Example 1: Deterministic sensitivity of contact problem in a com-

posite

Deterministic analysis and sensitivity of contact stresses in a two compo-
nent composite with spherical particles is verified with respect to the verti-
cal spatial coordinate. The following data are adopted for the computational
analysis: E2 = 2.0 · 10

9, ν1 = 0.3, ν2 = 0.2, R2 = 0.18, P = 10.0 · 10
5,

α=E1/E2 =2.0, β=R1/R2 =1.001−1.01.

Fig. 3. Contact stresses in the spherical particle-reinforced composite



Stochastic reliability in contact problems... 551

Fig. 4. Sensitivity of contact stresses to vertical spatial coordinate z

The computational analysis of the vertical contact stresses and their sen-
sitivity gradients dσz/dz with respect to the spatial coordinate is presented
in Fig.3 and Fig.4, this coordinate and the radii ratio β are marked on the
vertical axes of these figures.We observe that the compressive contact stresses
aremost sensitive to the parameter β for its value tending to ∞ (progressive
separation of the reinforcing particle from the surroundingmatrix).Moreover,
it is visible that the results is quite nonsensitive with respect to the vertical
coordinate zwhat cannot bepredicted fromthedefinition of σz, cf Eq.(2.10).
One of the main benefits flowing from making use of the package MAPLE is
visualization of variations of the stresses and their sensitivity gradients, which
can be studied in these figures.

5.2. Example 2: Stochastic modeling of the contact phenomenon

All the input parameters of the analyzed contact problem are now treated
as random variables: Young’s moduli and Poisson’s ratios of the composite
components as well as their radii. The deterministically calculated vertical
stresses in the contact area are compared below with the expected values,
standard deviations and probabilistic envelope values of the vertical contact
stresses for z=0.018 and z=−0.5. The standard deviations of the variables
are taken in the range of 10% of the corresponding expectations; all variables
are assumed to be uncorrelated.

The computed deterministic contact stresses are shown in Fig.5, Fig.9,



552 M.Kamiński

Fig. 5. Deterministic contact stresses; z=0.018

Fig. 6. Expected values of contact stresses; z=0.018



Stochastic reliability in contact problems... 553

Fig. 7. Standard deviations of contact stresses; z=0.018

Fig. 8. Probabilistic envelope of contact stresses; z=0.018



554 M.Kamiński

Fig. 9. Deterministic contact stresses; z=−0.5

Fig. 10. Expected values of contact stresses; z=−0.5



Stochastic reliability in contact problems... 555

Fig. 11. Standard deviations of contact stresses; z=−0.5

Fig. 12. Probabilistic envelope of contact stresses; z=−0.5



556 M.Kamiński

Fig. 13. Deterministic contact stresses; z=−2.0

Fig. 14. Expected values of contact stresses; z=−2.0



Stochastic reliability in contact problems... 557

Fig. 15. Standard deviations of contact stresses; z=−2.0

Fig. 16. Probabilistic envelope of contact stresses; z=−2.0



558 M.Kamiński

Fig.13 for the particle center z = 0.018, for some point belonging to the
matrix z = −0.5 and z = −2.0, respectively. The expected values of the
contact stresses are shown in Fig.6, Fig.10, Fig.14 for the same points, the
standarddeviations are shown inFig.7, Fig.11, Fig.15,while the probabilistic
envelopes for the surfaces of these stresses are presented in Fig.8, Fig.12 and
Fig.16. The vertical contact stresses are marked on the vertical axes; the
horizontal axes define the reinforcement ratio of the composite α and the
ratio between theparticle and the surroundingmatrix radii β. All the surfaces
shown in these figures have the same character and variability with respect to
the input parameters α and β, apart from the standard deviations plots.
Analyzing Figures 5, 6, 9, 10, 13 and 14, it can be seen that the expected

values of the contact stress surfaces are quite close to those obtained from
the corresponding deterministic analyses. Essential differences are observed
betweenFigures 5-8, 9-12 and13-16, where the probabilistic envelopes of these
stresses are shown. These envelopes are determined for a particular z on the
basis of the results presented in Figures 6, 7, 10, 11, 14 and 15 as

Env
(

f(z);z
)

=E[f(z);z]−3σ
(

f(z);z
)

(5.1)

Let us note that Eq. (5.1) is frequently used in Stochastic Finite Element
computations and stochastic fatigue analysis (see Liu et al., 1986). The values
of the probabilistic envelope surfaces are significantly smaller than the corre-
sponding values obtained from deterministic analysis, which means that the
computational analysis based on the stochastic perturbation in more restric-
tive than that of the classical model. The same concerns the corresponding
expected values. All the surfaces combined in the probability envelope show
that the vertical stresses tend to 0 for the reinforcement ratio approaching 1
and the matrix assuming the radius of the spherical particle. Comparing the
deterministic and stochastic results, one can clearly see that the contact stres-
ses are most sensitive to the vertical spatial coordinate.
Analyzing Figures 5-8, 10-12 and 14-16 in terms of variations of the con-

tact stresses with respect to the composite reinforcement ratio, it is observed
that the greatest sensitivity appears for α → 1, which means that the gre-
atest variations of the examined probabilistic stresses are obtained for the
homogeneous contact problem. The sensitivity analysis requires further nu-
merical examinations aimed at interrelations between Poisson’s ratios of both
composite components.

5.3. Example 3: Reliability analysis for the particle reinforced composite

The computational study on structural reliability, proposed in the the-



Stochastic reliability in contact problems... 559

oretical considerations, is the main subject of the next example. The set of
input data together with their probabilistic characteristics is given in Table 1
for the same composite contact problem as before. The Weibull probability
density function (PDF) of the limit function is determined together with its
probabilistic moments of up to 3rd order (cf Table 1) obtained from symbolic
computational solution to nonlinear equations system (4.8).

Table 1.Probabilistic input data for reliability index computations

Parameter Value

E2 2.0 ·10
6

ν1 0.3

ν2 0.2

R2 1.8

P 5.0 ·102

z −0.018

σ(E2) 0.2 ·10
6

S(E2) 0.0

σ(R2) 0.018

S(R2) 0.0

σall −4.0 ·10
5

α(σall) 10.0

β(σall) 1.01

E[g] −211378.33

σ(g) 38213.61838 (α=0.18)

S(g) 5.158577

First, it can be seen that even for a relatively small input coefficient of
variation of the input parameters (not greater than 0.1), the randomness level
of the output function is about 18% of the relevant expected value. That is
why the proposed third moment approach is more accurate for the analyzed
contact problem. Furthermore, we observe that even for the input skewnesses
equal to 0, the corresponding third order probabilistic characteristics differ
from 0, which reflects the differences in algebraic combinations of characteri-
stics of lower order, cf Eq. (4.7). In further analysis it is necessary to verify



560 M.Kamiński

the sensitivity (both in the deterministic and stochastic context, see Haug et
al. (1986), Kleiber and Hien (1992)) of the outputWeibull PDF probabilistic
momentswith respect to all inputmechanical and geometrical parameters and
their random characteristics. At the same time, the cross-correlation function
of the contact stresses Cov(σ(z1),σ(z2)) can be symbolically computed using
the programMAPLE.

6. Conclusions

• The above presented analysis reveals various sources of randomness and
stochasticity in contact problemsof the spherical particle reinforced com-
posites. In comparison to the second order second probabilistic moment
approach, third order probabilistic moments of both input and output
parameters are analyzed. It is demonstrated that even for the skewnesses
of the inputs equal to 0, the output third order probabilisticmoments in
reliability studies slightly differ from 0. It results from the main idea of
the SOTM approach and from the interrelations between probabilistic
characteristics of the lower order.Furthermore, it is observed that thede-
terministic values of the state functions are quite close to the computed
expected values. They are considerably greater and well approximated
by their probabilistic envelopes, which confirms the usefulness of these
envelopes in various stochastic numerical experiments.

• Themost interesting extension of this study would be introducing:

– randomness of non-spherical (ellipsoidal) contact surface

– more realistic incremental Stochastic Finite or Boundary Element
Method (SFEM, see Kamiński (2001), Kleiber and Hien, (1992) or
SBEM, see Burczyński (1995), respectively) analysis of nonlinear
geometry of the contacting surface, see Buczkowski and Kleiber
(1999).

• Application of the computational W-SOTM reliability study to various
numerical analyses of composites would be interesting, too. Neglecting
the relatively simple character of the deterministic contact problem, the
geometrical sensitivity of the contact stresses values is decisive for this
analysis, both in the deterministic and stochastic cases. Considering the
above, one can conclude that the stochastic second order perturbation
analysis in conjunction with mathematical symbolic computations is a



Stochastic reliability in contact problems... 561

powerful stochastic computational tool. However, the limitations on the
input randomness level, typical for such an analysis, must be fulfilled.

References

1. Beyerlein I.J., Phoenix S.L., 1997, Statistics of Fracture for an Elastic
Notched Composite Lamina Containing Weibull Fibers, Engng. Fract. Mech.,
57, 2/3, 267-299

2. Brandt A., Jendo S., Marks W., 1984, Probabilistic Approach to
Reliability-BasedOptimum Structural Design,Engng. Trans., 32, 1, 57-74

3. Buczkowski R., Kleiber M., 1999, A Stochastic Model of Rough Surfaces
for Finite ElementContactAnalysis,Comput. Meth. Appl. Mech. Engng., 169,
43-59

4. Burczyński T., 1995, Boundary Element Method in Mechanics (in Polish),
WNT, Warsaw

5. CharB.W., GeddesK.O., GonnetG.H., LeongB.L.,MonaganM.B.,
Watt S.M., 1992,First Leaves: ATutorial Introduction toMaple V, Springer-
Verlag

6. Fish J., edit., 1999, Computational Advances in Modeling Composites and
HeterogeneousMaterials,Comput. Meth. Appl. Mech. Engng., 172

7. GhanemR.G., Spanos P.D., 1997, Spectral Techniques for Stochastic Finite
Elements,Arch. Comput. Meth. Engng., 4, 1, 63-100

8. Ghanem R.G., Spanos P.D., 1992, Spectral Stochastic Finite-Element For-
mulation for Reliability Analysis, J. Engng. Mech. ASCE, 117, 10, 2351-2372

9. Grigoriu M., 2000, StochasticMechanics, Int. J. Sol. Struct., 37, 197-214

10. Haug E.J., Choi K.K., Komkov V., 1986, Design Sensitivity Analysis of
Structural Systems, Ser. Math. Sci. Engng., Academic Press

11. Kamiński M., 2001, Stochastic Second Order Perturbation Approach to the
Stress-Based Finite ElementMethod, Int. J. Sol. Struct., 38, 21, 3831-3852

12. Kleiber M., Hien T.D., 1992, Stochastic Finite Element Method, Wiley

13. LiuW.K.,BelytschkoT.,ManiA., 1986,ProbabilisticFinite Elements for
Nonlinear Structural Dynamics,Comput. Meth. Appl. Mech. Engng., 56, 61-81

14. LundE., 1998, ShapeOptimizationUsingWeibull Statistics of Brittle Failure,
Struct. Optimiz., 15, 208-214



562 M.Kamiński

15. Peng X.Q., Geng L., Liyan W., Liu G.R., Yam K.Y., 1998, A Stochastic
FiniteElementMethod forFatigueReliabilityAnalysis ofGearTeethSubjected
to Bending,Comput. Mech., 21, 253-261

16. TimoshenkoS.,GoodierS.N., 1951,Theory of Elasticity,McGraw-Hill, Inc.

Stochastyczna analiza niezawodności w zagadnieniu kontaktowym

kompozytu ze sferycznymi inkluzjami

Streszczenie

Tematem niniejszej pracy jest zastosowanie stochastycznej metody perturbacji
w analizie niezawodności materiałów kompozytowych ze sferycznymi inkluzjami.
W analizie tej wykorzystano stochastyczną metodę perturbacji drugiego rzędu
uwzględniając momenty probabilistyczne rzędu trzeciego. Zastosowane podejście
umożliwiło potraktowanie wybranych parametrów materiałowych i geometrycznych
jako zmienne losowe i wyznaczenie wartości oczekiwanych oraz odchyleń standardo-
wychnaprężeńkontaktowych.Obliczenianumerycznewykonanoprzypomocyprogra-
mu matematycznego do obliczeń symbolicznych MAPLE. Dzięki jego zastosowaniu
wykonano wizualizację wybranych elementów losowych, a także gradientów wrażli-
wości naprężeń kontaktowych. Zaproponowane w pracy podejście można rozszerzyć
na rozwiązania analityczne zagadnień kontaktowych o bardziej skomlikowanej geo-
metrii, a także w odpowiednich implementacjach numerycznych Metody Elementów
Skończonych lubMetody ElementówBrzegowych.

Manuscript received February 8, 2001; accepted for print April 27, 2001