# Application of Galerkin Method to Kirchhoff Plates Stochastic Bending Problem

Authors: Cláudio Roberto Ávila da Silva Júnior; Milton Kist; Marcelo Borges dos Santos
Journal: ISRN Applied Mathematics
 (2014)

---

## Abstract

In this paper, the Galerkin method is used to obtain approximate solutions for Kirchhoff plates stochastic bending problem with uncertainty over plates flexural rigidity coefficient. The uncertainty in the rigidity coefficient is represented by means of parameterized stochastic processes. A theorem of Lax-Milgram type, about existence and uniqueness of the theoretical solutions, is presented and used in selection of the approximate solution space. The Wiener-Askey scheme of generalized polynomials chaos (gPC) is used to model the stochastic behavior of the displacement solutions. The performance of the approximate Galerkin solution scheme developed herein is evaluated by comparing first and second order moments of the approximate solution with the same moments evaluated from Monte Carlo simulation. Rapid convergence of approximate Galerkin's solution to the first and second order moments is observed, for the problems studied herein. Results also show that using the developed Galerkin's scheme one gets adequate estimates for accrued probability function to a random variable generated by the stochastic process of displacement.


---

## Body

## 1. Introduction

The field of stochastic mechanics has been subject of extensive research and significant developments in recent years. Stochastic mechanics incorporates the modeling of randomness or uncertainty in the mathematical formulation of mechanics problems. This is in contrast to the more established field of structural reliability, where uncertainty and randomness are also addressed, but where problem solutions are obtained mainly based on deterministic mechanics models. The stochastic analysis applied to engineering structural systems, in addition to providing more information, can lead in some applications to more knowledge on the behavior of mathematic models used in presence of uncertainty.The analysis of stochastic engineering systems has received new impulse with use of finite element methods to obtain response statistics. Initially, finite element solutions were combined with the Monte Carlo method, and samples of random system response were obtained. This methodology is based on samples generation and, from there, getting the system’s set of realizations and, then, to get the estimates of responses and histograms statistics with the set of achievements. This procedure demands, due to samples dimension, much computing effort and time. In order to minimize this inconvenience, but following the same approach, the work of Yamazaki et al. [1] one may be quoted, which used Neumann’s series along with Monte Carlo’s methods and finite elements to reduce samples dimension, needed to determine statistical estimates of response; thus one may quote also the work of Araújo and Awruch [2], who applied the technique to structural problems subject to static and dynamic load and the geometrical nonlinearity with uncertainty on mechanical properties of materials.An alternative to get the numerical solutions to stochastic problems without samples generation and realization, emerged from an association of perturbation and finite elements methods, [2–5]. In this method, uncertainty on parameters of the system is represented through an expansion to the second order and through Galerkin’s type method equations are gotten for the stochastic solution process. One may quote Sobczyk’s [6] works with this methodology, in order to study the free vibrations problem from a rectangular plate with Young’s module modeled as a random field; Nakagiri et al. [4] formulated the problem with self-values for thin plates composite materials, with uncertainty on fibers’ direction; Ramu [7] used Galerkin’s technique to study the free vibration problem with Young’s module and mass distribution modeled by random fields; Kaminski [5] uses the principle of minimum complementary energy and a second order Taylor’s series for uncertainty representation in mechanical properties of the problem and to get the stochastic process of tensions. This methodology presents good performance for problems in which uncertainty presents small dispersion (
μ
/
σ
≤
0.1
) in system’s parameters.At the end of the 80s, Spanos and Ghanem [8] used the Galerkin finite element method to solve a stochastic beam bending problem, where Young’s modulus was modeled as a Gaussian stochastic process. The space of approximate solutions was built using the finite element method and chaos polynomials. These polynomials form a complete system in L
2
(
Ω
,
F
,
P
)
=
Ψ
¯
L
2
(
Ω
,
F
,
P
), where Ψ
=
span
⁡
〈
{
ψ
i
}
i
=
0
∞
〉 is the space generated by the chaos polynomials and (
Ω
,
F
,
P
) is a probability space. The following works stand out with this methodology: Ghanem and Spanos [9] used the Galerkin method and the Karhunem-Loeve series to represent uncertainty in the bending modulus by means of a Gaussian process. The ideas presented in this study were innovative and represented a new method to solve stochastic problems. Particularly, the construction of approximation space through Galerkin method with polynomials chaos provides that the stochastic approach to have similarities with the deterministic approach, thus becoming a motivating factor for using this methodology.Babuška et al. [10] presented a stochastic version of the Lax-Milgram lemma. The authors showed that, for certain problems of mechanics, use of Gaussian processes can lead to loss of coercivity of the bilinear form associated with the stochastic problem. This difficulty was indeed encountered in the study of da Silva Jr. [11] and resulted in nonconvergence of the solution for the bending of plates with random parameters. This nonconvergence was due to the choice of a Gaussian process to represent the uncertainty in some parameters of the system. This failure to converge also affects solutions based on perturbation or simulation methods. Although the message of Babuška et al. was assimilated by some researchers [12–15], one finds very recent papers where Gaussian processes are still used to represent the uncertainty in strictly positive mechanical properties [16–18].In recent years, much effort is being addressed at representing uncertainty in stochastic engineering systems via non-Gaussian processes. In the paper by Xiu and Karniadakis [19] the Askey-Wiener scheme was presented. This scheme represents a family of polynomials, denominated as generalized polynomials chaos (gPC's), which generate dense probability spaces with limited and unlimited support. This increases the possibilities for uncertain system parameter modeling. The use of gPC's has increased, from engineering applications standpoint, since they were used firstly by Ghanem; while one can quote several works, in several areas of engineering, such as those regarding plates bending problem especially, there are many works using the perturbation method [20–24]; however only few used gPC's. Chen and Soares [16] used Kharunen-Loeve’s expansion and the gPC’s to get numerical solutions for the plates bending problem in composite material; Vanmarcke and Grigoriu [25] studied the bending of Timoshenko beams with random shear modulus; Elishakoff et al. [26] employed the theory of mean square calculus to construct a solution to the boundary value problem of beam bending with stochastic bending modulus; Chakraborty and Sarkar [27] used the Neumann series and Monte Carlo simulation to obtain statistical moments of the displacements of curved beams, with uncertainty in the elasticity modulus of the foundation. Although they present numerical solutions for stochastic beam problems, none of the papers referenced above address the matter of existence and uniqueness of the solutions.In the present paper, the Galerkin method is used to obtain approximated, numerical solutions for Kirchhoff plates bending problem with uncertainty in flexural rigidity coefficient. Uncertainties in the flexural rigidity coefficient are represented by parameterized random processes [28]. A major contribution from this work is presented by a theoretical result in the way of a theorem on existence and uniqueness of the solution for plates stochastic bending problem. This result uses the Lax-Milgram lemma [10]. This study subsidizes and gives consistency to the Galerkin solution scheme developed herein, avoiding the flaws of some formulations mentioned earlier [16–18]. From this point, one uses Galerkin method to get the approximate solutions (numerical) for the displacement stochastic process. The tensioning product between spaces of finite dimension generates the space of approximate solutions. From the isomorphism and density properties assures that the approximate solutions spaces is dense in the theoretical solutions space. A subspace of the Askey-Wiener’s scheme [19, 29] is used to represent the stochastic behavior of the solution. Performance of the Galerkin solution developed herein is evaluated by comparing first and second order moments of random beam displacement responses with the same statistics evaluated via Monte Carlo simulation. One compare still the estimates for the accrued probability function, and from Monte Carlo stimulation.


## 2. Kirchhoff Plates Stochastic Bending Problem

In this section, are presented the strong and weak formulations of the problem of stochastic bending of Kirchhoff plates. At the end of this section, the Lax-Milgram lemma is used to present a proof of existence and uniqueness of the solution. Let(
Ω
,
F
,
P
) be a probability space, where Ω is a sample space, F is an σ-algebra, and P is a probability measure. The strong form of the Kirchhoff’s plates bending problem is
(1)
Δ
(
α
(
x
,
ω
)
·
Δ
u
)
=
q
,
∀
(
x
,
ω
)
∈
D
×
Ω
,
a
.
e
.
;
u
(
x
,
ω
)
=
0
,
Δ
u
(
x
,
ω
)
=
0
,
∀
(
x
,
ω
)
∈
Γ
×
Ω
,
where u is stochastic processes of the transversal plate displacement and q is a load term. The term α
(
·
,
·
) is the plate flexural rigidity coefficient, which is given by
(2)
α
(
x
,
ω
)
=
(
E
·
t
3
12
·
(
1
-
ν
2
)
)
(
x
,
ω
)
,
where E is Young’s modulus, t is plate thickness, and ν is Poisson’s coefficient.For the analysis of existence and uniqueness of the response, the following hypotheses are considered:(3)
(
H
1
)
∃
α
_
,
α
¯
∈
R
+
∖
{
0
}
:
P
(
ω
∈
Ω
:
α
(
x
,
ω
)
∈
[
α
_
,
α
¯
]
,
:
∃
α
_
,
H
1
α
¯
∈
R
+
∖
{
0
}
:
P
w
∀
x
∈
D
[
α
_
,
α
¯
]
)
=
1
;
(
H
2
)
α
∈
L
2
(
Ω
,
F
,
P
;
H
2
(
D
)
)
;
(
H
3
)
q
∈
L
2
(
Ω
,
F
,
P
;
L
2
(
D
)
)
.
Hypothesis (H1) ensures that the plate flexural rigidity coefficient is positively defined and uniformly limited in probability [10], whereas the hypothesis (H2) is that the strong form is well-defined. Hypothesis (H3) ensures that the stochastic load process has finite variance. These hypotheses are necessary for the application of the Lax-Milgram lemma, which is used in the sequence to demonstrate the existence and uniqueness of the solution.


### 2.1. Existence and Uniqueness of the Solution

In this section, a brief theoretical study of existence and uniqueness of the solution for Kirchhoff’s plate bending stochastic problem with uncertainty on the bending rigidity is discussed. Results presented in this section are based on classical results from functional analysis and from theory of Sobolev spaces [10, 30, 31]. The study requires definition of stochastic Sobolev spaces, tensor product, and density between distribution spaces and L
p spaces. In order to study existence and uniqueness, the abstract variational problem associated with the strong form (see (1)) needs to be defined. The stochastic Sobolev space, where the solution to the stochastic beam bending problem is constructed, is V
=
L
2
(
Ω
,
F
,
P
;
Q
). For fixed ω
∈
Ω, one has(4a)
Q
=
{
u
(
·
,
ω
)
∈
H
2
(
D
)
∣
u
(
x
,
ω
)
=
0
∧
Δ
u
(
x
,
ω
)
=
0
,
∀
(
x
,
ω
)
∈
Γ
{
u
(
·
,
ω
)
∈
H
2
(
D
)
∣
u
(
x
,
ω
)
=
0
∧
Δ
u
(
x
,
ω
)
=
0
,
}
;
hence
(4b)
V
=
{
∫
Ω
∥
u
(
ω
)
∥
H
2
(
D
)
2
d
P
(
ω
)
u
:
D
×
Ω
⟶
R
∣
u
is
mensurable
and
∫
Ω
∥
u
(
ω
)
∥
H
2
(
D
)
2
d
P
(
ω
)
<
+
∞
}
,where ∥
·
∥
H
2
(
D
) is the H
2
(
D
) norm. Equation (4b) means that, for fixed ω
∈
Ω, u
(
·
,
ω
)
∈
Q. On the other hand, for fixed x
∈
Γ, u
(
x
,
·
)
∈
L
2
(
Ω
,
F
,
P
) is a random variable. Defining the tensor product between g
∈
L
2
(
Ω
,
F
,
P
) and w
∈
Q as u
∶
=
w
·
g [32, 33], one notes that, for fixed ω
∈
Ω,
(5)
u
(
·
,
ω
)
=
w
(
·
)
·
g
(
ω
)
∈
Q
,
whereas for a fixed x
∈
D,
(6)
u
(
x
,
·
)
=
w
(
x
)
·
g
(
·
)
∈
L
2
(
Ω
,
F
,
P
)
.
In order to obtain the isomorphism between V and L
2
(
Ω
,
F
,
P
)
⊗
Q, it is necessary to redefine the differential operator Δ
ω
:
L
2
(
Ω
,
F
,
P
)
⊗
Q
→
L
2
(
Ω
,
F
,
P
)
⊗
L
2
(
D
). This operator acts over an element u
∈
V in the following way [34]:
(7)
Δ
ω
u
=
(
Δ
w
)
(
x
)
·
g
(
ω
)
,
where α
∈
N and α
≤
2. From the definitions presented in (5) to (6), the isomorphism between V and L
2
(
Ω
,
F
,
P
)
⊗
Q is obtained:
(8)
V
=
L
2
(
Ω
,
F
,
P
;
Q
)
L
2
(
Ω
,
F
,
P
;
Q
)
≃
L
2
(
Ω
,
F
,
P
)
⊗
Q
⇓
V
≃
L
2
(
Ω
,
F
,
P
)
⊗
Q
.
V is a Hilbert space, with inner product (
·
,
·
)
V
:
V
×
V
→
R, defined as
(9)
(
u
,
v
)
V
=
∫
Ω
∫
D
(
Δ
ω
u
·
Δ
ω
v
)
(
x
,
ω
)
d
x
d
P
(
ω
)
.
The inner product defined in (9) induces the V norm ∥
u
∥
V
=
(
u
,
u
)
V
1
/
2, following Kinderlehrer and Stampacchia [35]. The bilinear form a
:
V
×
V
→
R is defined as
(10)
a
(
u
,
v
)
=
∫
Ω
∫
D
(
α
·
Δ
ω
u
·
Δ
ω
v
)
(
x
,
ω
)
d
x
d
P
(
ω
)
.
The bilinear form a
(
·
,
·
) is obtained by multiplying (1) by v
∈
V and integrating in D
×
Ω, using the product rule for differentiation and the definition of space V given in (4b). In order to preserve the equality, the right-hand side of (1) becomes
(11)
f
(
v
)
=
∫
Ω
∫
D
(
q
·
v
)
(
x
,
ω
)
d
x
d
P
(
ω
)
.
It is easy to see that (11) defines a linear functional f
:
V
→
R. Using the Cauchy-Schwartz inequality, it can be shown that this functional is well-defined; hence f
∈
V
′, where V
′ is the dual of V. Hence, the abstract variational problem associated with (1) can be written as
(12)
Find
u
∈
V
such
that
a
(
u
,
v
)
=
f
(
v
)
,
∀
v
∈
V
.Theorem 1 (existence and uniqueness).
Letα and q be such that (H1)–(H3) ((2) and (3)) are satisfied. Then, a solution to the problem defined in (12) exists, and it is unique in V.Proof.
The proof of existence and uniqueness uses the Lax-Milgram lemma [10, 36]. It is necessary to show that the bilinear form (10) of the problem defined in (12) is continuous and coercive. In order to do that, one makes use of hypothesis (H1) of limited probability (3) of the Cauchy-Schwartz inequality [30]:(a)
continuity(13)
|
a
(
u
,
v
)
|
≤
α
¯
[
∫
Ω
∫
D
(
Δ
ω
u
)
2
(
x
,
ω
)
d
x
d
P
(
ω
)
]
1
/
2
×
[
∫
Ω
∫
D
(
Δ
ω
v
)
2
(
x
,
ω
)
d
x
d
P
(
ω
)
]
1
/
2
≤
C
∥
u
∥
V
∥
v
∥
V
,
where C
=
α
¯;
(b)
coercivity(14)
a
(
u
,
u
)
≥
α
_
∫
Ω
∫
D
(
Δ
ω
u
)
2
(
x
,
ω
)
d
x
d
P
(
ω
)
≥
c
∥
u
∥
V
2
,
where c
=
α
_.Following (11), the functional f
:
V
→
R is linear and well-defined; hence, in view of the Lax-Milgram lemma, it is guaranteed that the problem defined in (12) has unique solution and continuous dependency on the data [10, 36].Theoretical solutions to the abstract variational problem in (1), associated with the problem in (12), are found in V
≃
L
2
(
Ω
,
F
,
P
)
⊗
Q. Numerical solutions are obtained in less abstract spaces: continuous functions of class C
2, sequentially dense in Q, and a family of generalized chaos polynomials, belonging to the Askey-Wiener scheme, are used.


## 2.1. Existence and Uniqueness of the Solution

In this section, a brief theoretical study of existence and uniqueness of the solution for Kirchhoff’s plate bending stochastic problem with uncertainty on the bending rigidity is discussed. Results presented in this section are based on classical results from functional analysis and from theory of Sobolev spaces [10, 30, 31]. The study requires definition of stochastic Sobolev spaces, tensor product, and density between distribution spaces and L
p spaces. In order to study existence and uniqueness, the abstract variational problem associated with the strong form (see (1)) needs to be defined. The stochastic Sobolev space, where the solution to the stochastic beam bending problem is constructed, is V
=
L
2
(
Ω
,
F
,
P
;
Q
). For fixed ω
∈
Ω, one has(4a)
Q
=
{
u
(
·
,
ω
)
∈
H
2
(
D
)
∣
u
(
x
,
ω
)
=
0
∧
Δ
u
(
x
,
ω
)
=
0
,
∀
(
x
,
ω
)
∈
Γ
{
u
(
·
,
ω
)
∈
H
2
(
D
)
∣
u
(
x
,
ω
)
=
0
∧
Δ
u
(
x
,
ω
)
=
0
,
}
;
hence
(4b)
V
=
{
∫
Ω
∥
u
(
ω
)
∥
H
2
(
D
)
2
d
P
(
ω
)
u
:
D
×
Ω
⟶
R
∣
u
is
mensurable
and
∫
Ω
∥
u
(
ω
)
∥
H
2
(
D
)
2
d
P
(
ω
)
<
+
∞
}
,where ∥
·
∥
H
2
(
D
) is the H
2
(
D
) norm. Equation (4b) means that, for fixed ω
∈
Ω, u
(
·
,
ω
)
∈
Q. On the other hand, for fixed x
∈
Γ, u
(
x
,
·
)
∈
L
2
(
Ω
,
F
,
P
) is a random variable. Defining the tensor product between g
∈
L
2
(
Ω
,
F
,
P
) and w
∈
Q as u
∶
=
w
·
g [32, 33], one notes that, for fixed ω
∈
Ω,
(5)
u
(
·
,
ω
)
=
w
(
·
)
·
g
(
ω
)
∈
Q
,
whereas for a fixed x
∈
D,
(6)
u
(
x
,
·
)
=
w
(
x
)
·
g
(
·
)
∈
L
2
(
Ω
,
F
,
P
)
.
In order to obtain the isomorphism between V and L
2
(
Ω
,
F
,
P
)
⊗
Q, it is necessary to redefine the differential operator Δ
ω
:
L
2
(
Ω
,
F
,
P
)
⊗
Q
→
L
2
(
Ω
,
F
,
P
)
⊗
L
2
(
D
). This operator acts over an element u
∈
V in the following way [34]:
(7)
Δ
ω
u
=
(
Δ
w
)
(
x
)
·
g
(
ω
)
,
where α
∈
N and α
≤
2. From the definitions presented in (5) to (6), the isomorphism between V and L
2
(
Ω
,
F
,
P
)
⊗
Q is obtained:
(8)
V
=
L
2
(
Ω
,
F
,
P
;
Q
)
L
2
(
Ω
,
F
,
P
;
Q
)
≃
L
2
(
Ω
,
F
,
P
)
⊗
Q
⇓
V
≃
L
2
(
Ω
,
F
,
P
)
⊗
Q
.
V is a Hilbert space, with inner product (
·
,
·
)
V
:
V
×
V
→
R, defined as
(9)
(
u
,
v
)
V
=
∫
Ω
∫
D
(
Δ
ω
u
·
Δ
ω
v
)
(
x
,
ω
)
d
x
d
P
(
ω
)
.
The inner product defined in (9) induces the V norm ∥
u
∥
V
=
(
u
,
u
)
V
1
/
2, following Kinderlehrer and Stampacchia [35]. The bilinear form a
:
V
×
V
→
R is defined as
(10)
a
(
u
,
v
)
=
∫
Ω
∫
D
(
α
·
Δ
ω
u
·
Δ
ω
v
)
(
x
,
ω
)
d
x
d
P
(
ω
)
.
The bilinear form a
(
·
,
·
) is obtained by multiplying (1) by v
∈
V and integrating in D
×
Ω, using the product rule for differentiation and the definition of space V given in (4b). In order to preserve the equality, the right-hand side of (1) becomes
(11)
f
(
v
)
=
∫
Ω
∫
D
(
q
·
v
)
(
x
,
ω
)
d
x
d
P
(
ω
)
.
It is easy to see that (11) defines a linear functional f
:
V
→
R. Using the Cauchy-Schwartz inequality, it can be shown that this functional is well-defined; hence f
∈
V
′, where V
′ is the dual of V. Hence, the abstract variational problem associated with (1) can be written as
(12)
Find
u
∈
V
such
that
a
(
u
,
v
)
=
f
(
v
)
,
∀
v
∈
V
.Theorem 1 (existence and uniqueness).
Letα and q be such that (H1)–(H3) ((2) and (3)) are satisfied. Then, a solution to the problem defined in (12) exists, and it is unique in V.Proof.
The proof of existence and uniqueness uses the Lax-Milgram lemma [10, 36]. It is necessary to show that the bilinear form (10) of the problem defined in (12) is continuous and coercive. In order to do that, one makes use of hypothesis (H1) of limited probability (3) of the Cauchy-Schwartz inequality [30]:(a)
continuity(13)
|
a
(
u
,
v
)
|
≤
α
¯
[
∫
Ω
∫
D
(
Δ
ω
u
)
2
(
x
,
ω
)
d
x
d
P
(
ω
)
]
1
/
2
×
[
∫
Ω
∫
D
(
Δ
ω
v
)
2
(
x
,
ω
)
d
x
d
P
(
ω
)
]
1
/
2
≤
C
∥
u
∥
V
∥
v
∥
V
,
where C
=
α
¯;
(b)
coercivity(14)
a
(
u
,
u
)
≥
α
_
∫
Ω
∫
D
(
Δ
ω
u
)
2
(
x
,
ω
)
d
x
d
P
(
ω
)
≥
c
∥
u
∥
V
2
,
where c
=
α
_.Following (11), the functional f
:
V
→
R is linear and well-defined; hence, in view of the Lax-Milgram lemma, it is guaranteed that the problem defined in (12) has unique solution and continuous dependency on the data [10, 36].Theoretical solutions to the abstract variational problem in (1), associated with the problem in (12), are found in V
≃
L
2
(
Ω
,
F
,
P
)
⊗
Q. Numerical solutions are obtained in less abstract spaces: continuous functions of class C
2, sequentially dense in Q, and a family of generalized chaos polynomials, belonging to the Askey-Wiener scheme, are used.


## 3. Uncertainty Representation

In most engineering problems, complete statistical information about uncertainties is not available. Sometimes, the first and second moments are the only information available. The probability distribution function is defined based on experience or heuristically. In order to apply Galerkin’s method, an explicit representation of the uncertainty is necessary. Given the incomplete information about the probability distribution of a given parameter, a hypothesis of finite dimensional noise is assumed, following, for example, [37–39]. This implies that the uncertainty over a given input parameter ϑ
:
D
×
Ω
→
R
+ will be represented in terms of a finite set of random variables:
(15)
ϑ
(
x
,
ω
)
=
ϑ
(
x
,
ξ
(
ω
)
)
=
ϑ
(
x
,
ξ
1
(
ω
)
,
…
,
ξ
N
(
ω
)
)
.
From this hypothesis, the uncertainty in beam and foundation stiffness coefficients are modeled via parameterized stochastic processes. These are defined from a linear combination of deterministic functions and random variables [28]:
(16)
ϑ
(
x
,
ω
)
=
μ
ϑ
(
x
)
+
∑
i
=
1
N
φ
i
(
x
)
ξ
i
(
ω
)
=
μ
ϑ
(
x
)
+
Φ
t
(
x
)
ξ
(
ω
)
,
where μ
ϑ
(
·
) is the expected value of random process ϑ
(
·
,
·
) and Φ
:
D
→
R
N is a vector-valued function with terms φ
i
∈
C
0
(
D
)
∩
C
2
(
D
), ∀
i
∈
{
1
,
…
,
N
}. ξ
(
ω
)
=
{
ξ
i
(
ω
)
}
i
=
1
N is a vector of independent random variables, such that
(17)
E
[
ξ
i
]
=
0
,
∀
i
∈
{
1
,
…
,
N
}
;
P
(
ω
∈
Ω
:
ξ
i
(
ω
)
∈
Γ
i
)
=
1
,
∀
i
∈
{
1
,
…
,
N
}
,
where E
[
·
] is the expected value operator. In (17), Γ
i is the image of random variable ξ
i; that is, Γ
i
=
ξ
i
(
Ω
), with Γ
i
=
[
a
i
,
b
i
]
⊂
R, |
Γ
i
|
=
|
b
i
-
a
i
|
<
∞, ∀
i
∈
{
1
,
…
,
N
}, limited. In this form, the image of random vector ξ
:
Ω
→
Γ, with Γ
⊂
R
N and in terms of {
Γ
i
}
i
=
1
N, is given by Γ
=
∏
i
=
1
N
Γ
i. Since the random variables are independent, the joint probability density is given by
(18)
ρ
(
ξ
(
ω
)
)
=
∏
i
=
1
N
ρ
i
(
ξ
i
)
,
where ρ
i
(
·
) is the marginal probability density function of random variable ξ
i. Hence, the probability measure d
P
(
·
) is defined as
(19)
d
P
(
ξ
(
ω
)
)
=
∏
i
=
1
N
ρ
i
(
ξ
i
)
d
ξ
i
.
From the measure and integration theory [40], one knows that the probability measure defined in (19) is the product measure obtained from the product between probability measure spaces associated with the random variables ξ
(
ω
)
=
{
ξ
i
(
ω
)
}
i
=
1
N, with ξ
i
:
Ω
→
Γ
i.From the Doob-Dynkin lemma [41], the transversal displacement random process will be a function of random variables ξ
(
ω
)
=
{
ξ
i
(
ω
)
}
i
=
1
N; hence
(20)
u
(
x
,
ω
)
=
u
(
x
,
ξ
(
ω
)
)
=
u
(
x
,
ξ
1
(
ω
)
,
…
,
ξ
N
(
ω
)
)
.
In this paper, polynomials of the Askey-Wiener scheme are used to construct the problems solution space.


### 3.1. The Askey-Wiener Scheme

The Askey-Wiener scheme is a generalization of chaos polynomials, also known as Wiener-chaos. Chaos polynomials were proposed by Wiener [42] to study statistical mechanics of gases. Xiu and Karniadakis [19] have shown the close relationship between results presented by Wiener [42] and Askey and Wilson [43] for the representation of stochastic processes by orthogonal polynomials. Xiu and Karniadakis [19] extended the studies of Ghanem and Spanos [9] and Ogura [44] for polynomials belonging to the Askey-Wiener scheme.The Cameron-Martin theorem [29] shows that Wiener-Askey polynomials form a base for a dense subspace of second order random variables L
2
(
Ω
,
F
,
P
). Let H
⊆
L
2
(
Ω
,
F
,
P
) be a separable Gaussian Hilbert space and let H
=
span
⁡
〈
{
ξ
i
}
i
=
1
∞
〉 be an orthonormal base of Gaussian random variables. The Wiener-Askey scheme allows a representation of any second order random variable u
∈
L
2
(
Ω
,
Σ
(
H
)
,
P
), where Σ
(
H
) is the σ-algebra generated by H [45]. Let P
n
(
H
) be the vector space spanned by all polynomials of order less than n:
(21)
P
n
(
H
)

=

{
Γ
(
{
ξ
i
}
i
=
1
N
)
:
Γ
is
the
polynomial
of
degree
≤
n
;
ξ
i
∈
H
,
∀
i
=
1
,
…
,
N
;
N
<
∞
{
Γ
(
{
ξ
i
}
i
=
1
N
)
:
Γ
is
the
polynomial
of
degree
≤
n
;
}
,
with
(22)
H
:
0
:
=
P
¯
0
(
H
)
,
H
:
n
:
=
P
¯
n
(
H
)
∩
P
¯
n
-
1
(
H
)
⊥
,
where P
¯
n is the closure of P
n in L
2
(
Ω
,
F
,
P
). As shown by Janson [46], for L
2
(
Ω
,
Σ
(
H
)
,
P
), the following orthogonal decomposition is admitted:
(23)
L
2
(
Ω
,
Σ
(
H
)
,
P
)
=
⨁
n
=
0
∞
H
:
n
:
.
Equation (23) is an orthogonal decomposition of L
2
(
Ω
,
F
,
P
), known as Wiener-chaos decomposition or simply chaos decomposition. One application of this decomposition is the representation of an element X
∈
L
2
(
Ω
,
F
,
P
) in terms of elements X
n
∈
H
:
n
::
(24)
X
=
∑
n
=
0
∞
X
n
.
Equation (24) represents an important result for the approximation theory applied to stochastic systems. Solution of a stochastic system is expressed as a nonlinear function in terms Gaussian random variables. This function is expanded in terms of chaos polynomials as
(25)
u
i
(
ω
)
=
∑
p
≥
0
∑
n
1
+
⋯
+
n
r
=
p
∑
i
1
,
…
,
i
r
u
i
1
⋯
i
r
j
1
⋯
j
r
ψ
p
(
ξ
i
1
(
ω
)
,
…
,
ξ
i
r
(
ω
)
)
,
where ψ
p is the chaos polynomial of order “p” and u
i
1
⋯
i
r
j
1
⋯
j
r are the polynomial coefficients. The superindex refers to the number of occurrences of ξ
i
k
(
ω
). Chaos polynomials of order “p” are formed by a Hermite polynomial, in standard Gaussian variables ξ
(
ω
)
=
{
ξ
i
k
(
ω
)
}
k
=
1
r of order less than “p.” Introducing a mapping in the sets of indexes {
i
k
}
k
=
1
r and {
j
k
}
k
=
1
r, (25) can be rewritten as
(26)
u
(
ω
)
=
∑
λ
∈
Λ
u
λ
ψ
λ
(
ξ
(
ω
)
)
,
where λ is a multi-index, Λ
⊂
N is a set of natural numbers with compact support, {
ψ
λ
}
λ
∈
Λ are chaos polynomials, {
u
λ
}
λ
∈
Λ are coefficients of the linear combination, and ξ
:
Ω
→
R
N is a vector of random variables. In this paper, chaos polynomials are multidimensional Hermite polynomials:
(27)
ψ
λ
(
ξ
(
ω
)
)
=
∏
m
=
1
∞
h
λ
m
(
ξ
m
(
ω
)
)
,
where h
λ
m
(
·
) is a Hermite polynomial defined in random variable ξ
m. The inner product between polynomials ψ
i and ψ
j in L
2
(
Ω
,
F
,
P
) is defined as
(28)
(
ψ
i
,
ψ
j
)
L
2
(
Ω
,
F
,
P
)
=
∫
Ω
(
ψ
i
·
ψ
j
)
(
ξ
(
ω
)
)
d
P
(
ω
)
,
where d
P is a probability measure. These polynomials form a complete orthonormal system with respect to the probability measure, with the following properties:
(29)
ψ
0
=
1
,
(
ψ
i
,
ψ
j
)
L
2
(
Ω
,
F
,
P
)
=
δ
i
j
,
∀
i
,
j
∈
N
.
It is important to observe that in (28) the polynomials are orthogonal with respect to the standard Gaussian density function of vector ξ. The convergence rate is exponential for the case where the random variable is Gaussian. For other random variables the convergence rate is smaller.The proposal of the Wiener-Askey scheme is to extend the result presented in (25) to other types of polynomials. In analogy to (21), taking P
n
(
H
)
=
span
⁡
〈
{
ψ
i
}
i
=
1
N
〉, with H a separable Hilbert space of finite variance random variables, one has that P
=
⋃
n
∈
N
P
n
(
H
) is a family of polynomials of the Wiener-Askey scheme, generating a complete orthogonal system in L
2
(
Ω
,
F
,
P
). The Askey-Wiener scheme represents a family of subspaces generated by orthogonal polynomials obtained from ordinary differential equations [19]. Among them, the Hermite, Laguerre, Jacobi, and Legendre polynomials can be cited. Every subspace generated by these polynomials is a complete system in L
2
(
Ω
,
F
,
P
). The orthogonality between the polynomials is defined with respect to a weight function, which is identical to the probability density function of a certain random variable. For example, the Gaussian density function is used as weight function to obtain the orthogonality between Hermite polynomials. Table 1 shows the correspondence between subsets of polynomials of the Askey-Wiener scheme and the corresponding probability density functions.Table 1
Correspondence between some random variables and polynomials of the Askey-Wiener scheme.
Random variable
Polynomial
Weight function
Support
Gaussian
Hermite
e
-
|
ξ
|
2
/
2
R
Gama
Laguerre
1
Γ
(
α
+
1
)
ξ
α
e
-
ξ
[
0
,
+
∞
)
Beta
Jacobi
2
-
(
α
+
β
+
1
)
Γ
(
α
+
β
+
2
)
Γ
(
α
+
1
)
Γ
(
β
+
1
)
(
1
-
ξ
)
α
(
1
+
ξ
)
β
e
-
ξ
[
a
,
b
]
Uniform
Legendre
1
b
-
a
[
a
,
b
]


## 3.1. The Askey-Wiener Scheme

The Askey-Wiener scheme is a generalization of chaos polynomials, also known as Wiener-chaos. Chaos polynomials were proposed by Wiener [42] to study statistical mechanics of gases. Xiu and Karniadakis [19] have shown the close relationship between results presented by Wiener [42] and Askey and Wilson [43] for the representation of stochastic processes by orthogonal polynomials. Xiu and Karniadakis [19] extended the studies of Ghanem and Spanos [9] and Ogura [44] for polynomials belonging to the Askey-Wiener scheme.The Cameron-Martin theorem [29] shows that Wiener-Askey polynomials form a base for a dense subspace of second order random variables L
2
(
Ω
,
F
,
P
). Let H
⊆
L
2
(
Ω
,
F
,
P
) be a separable Gaussian Hilbert space and let H
=
span
⁡
〈
{
ξ
i
}
i
=
1
∞
〉 be an orthonormal base of Gaussian random variables. The Wiener-Askey scheme allows a representation of any second order random variable u
∈
L
2
(
Ω
,
Σ
(
H
)
,
P
), where Σ
(
H
) is the σ-algebra generated by H [45]. Let P
n
(
H
) be the vector space spanned by all polynomials of order less than n:
(21)
P
n
(
H
)

=

{
Γ
(
{
ξ
i
}
i
=
1
N
)
:
Γ
is
the
polynomial
of
degree
≤
n
;
ξ
i
∈
H
,
∀
i
=
1
,
…
,
N
;
N
<
∞
{
Γ
(
{
ξ
i
}
i
=
1
N
)
:
Γ
is
the
polynomial
of
degree
≤
n
;
}
,
with
(22)
H
:
0
:
=
P
¯
0
(
H
)
,
H
:
n
:
=
P
¯
n
(
H
)
∩
P
¯
n
-
1
(
H
)
⊥
,
where P
¯
n is the closure of P
n in L
2
(
Ω
,
F
,
P
). As shown by Janson [46], for L
2
(
Ω
,
Σ
(
H
)
,
P
), the following orthogonal decomposition is admitted:
(23)
L
2
(
Ω
,
Σ
(
H
)
,
P
)
=
⨁
n
=
0
∞
H
:
n
:
.
Equation (23) is an orthogonal decomposition of L
2
(
Ω
,
F
,
P
), known as Wiener-chaos decomposition or simply chaos decomposition. One application of this decomposition is the representation of an element X
∈
L
2
(
Ω
,
F
,
P
) in terms of elements X
n
∈
H
:
n
::
(24)
X
=
∑
n
=
0
∞
X
n
.
Equation (24) represents an important result for the approximation theory applied to stochastic systems. Solution of a stochastic system is expressed as a nonlinear function in terms Gaussian random variables. This function is expanded in terms of chaos polynomials as
(25)
u
i
(
ω
)
=
∑
p
≥
0
∑
n
1
+
⋯
+
n
r
=
p
∑
i
1
,
…
,
i
r
u
i
1
⋯
i
r
j
1
⋯
j
r
ψ
p
(
ξ
i
1
(
ω
)
,
…
,
ξ
i
r
(
ω
)
)
,
where ψ
p is the chaos polynomial of order “p” and u
i
1
⋯
i
r
j
1
⋯
j
r are the polynomial coefficients. The superindex refers to the number of occurrences of ξ
i
k
(
ω
). Chaos polynomials of order “p” are formed by a Hermite polynomial, in standard Gaussian variables ξ
(
ω
)
=
{
ξ
i
k
(
ω
)
}
k
=
1
r of order less than “p.” Introducing a mapping in the sets of indexes {
i
k
}
k
=
1
r and {
j
k
}
k
=
1
r, (25) can be rewritten as
(26)
u
(
ω
)
=
∑
λ
∈
Λ
u
λ
ψ
λ
(
ξ
(
ω
)
)
,
where λ is a multi-index, Λ
⊂
N is a set of natural numbers with compact support, {
ψ
λ
}
λ
∈
Λ are chaos polynomials, {
u
λ
}
λ
∈
Λ are coefficients of the linear combination, and ξ
:
Ω
→
R
N is a vector of random variables. In this paper, chaos polynomials are multidimensional Hermite polynomials:
(27)
ψ
λ
(
ξ
(
ω
)
)
=
∏
m
=
1
∞
h
λ
m
(
ξ
m
(
ω
)
)
,
where h
λ
m
(
·
) is a Hermite polynomial defined in random variable ξ
m. The inner product between polynomials ψ
i and ψ
j in L
2
(
Ω
,
F
,
P
) is defined as
(28)
(
ψ
i
,
ψ
j
)
L
2
(
Ω
,
F
,
P
)
=
∫
Ω
(
ψ
i
·
ψ
j
)
(
ξ
(
ω
)
)
d
P
(
ω
)
,
where d
P is a probability measure. These polynomials form a complete orthonormal system with respect to the probability measure, with the following properties:
(29)
ψ
0
=
1
,
(
ψ
i
,
ψ
j
)
L
2
(
Ω
,
F
,
P
)
=
δ
i
j
,
∀
i
,
j
∈
N
.
It is important to observe that in (28) the polynomials are orthogonal with respect to the standard Gaussian density function of vector ξ. The convergence rate is exponential for the case where the random variable is Gaussian. For other random variables the convergence rate is smaller.The proposal of the Wiener-Askey scheme is to extend the result presented in (25) to other types of polynomials. In analogy to (21), taking P
n
(
H
)
=
span
⁡
〈
{
ψ
i
}
i
=
1
N
〉, with H a separable Hilbert space of finite variance random variables, one has that P
=
⋃
n
∈
N
P
n
(
H
) is a family of polynomials of the Wiener-Askey scheme, generating a complete orthogonal system in L
2
(
Ω
,
F
,
P
). The Askey-Wiener scheme represents a family of subspaces generated by orthogonal polynomials obtained from ordinary differential equations [19]. Among them, the Hermite, Laguerre, Jacobi, and Legendre polynomials can be cited. Every subspace generated by these polynomials is a complete system in L
2
(
Ω
,
F
,
P
). The orthogonality between the polynomials is defined with respect to a weight function, which is identical to the probability density function of a certain random variable. For example, the Gaussian density function is used as weight function to obtain the orthogonality between Hermite polynomials. Table 1 shows the correspondence between subsets of polynomials of the Askey-Wiener scheme and the corresponding probability density functions.Table 1
Correspondence between some random variables and polynomials of the Askey-Wiener scheme.
Random variable
Polynomial
Weight function
Support
Gaussian
Hermite
e
-
|
ξ
|
2
/
2
R
Gama
Laguerre
1
Γ
(
α
+
1
)
ξ
α
e
-
ξ
[
0
,
+
∞
)
Beta
Jacobi
2
-
(
α
+
β
+
1
)
Γ
(
α
+
β
+
2
)
Γ
(
α
+
1
)
Γ
(
β
+
1
)
(
1
-
ξ
)
α
(
1
+
ξ
)
β
e
-
ξ
[
a
,
b
]
Uniform
Legendre
1
b
-
a
[
a
,
b
]


## 4. Galerkin Method

The Galerkin method is used in this paper to solve the stochastic beam bending problem with uncertainty in the beam and foundation stiffness coefficients. In order to develop numerical solutions which are compatible with the conditions for existence and uniqueness of the theoretical solution, results presented by Besold [32] and Ryan [47] are used. An element from a space isomorph to the space obtained via tensor product, V
≃
L
2
(
Ω
,
F
,
P
)
⊗
Q, can be represented from elements of separable spaces, dense in spaces L
2
(
Ω
,
F
,
P
), and Q. Hence, the strategy to construct approximated, numerical solutions and to use Galerkin’s method is to use bases of finite dimensions, but dense in L
2
(
Ω
,
F
,
P
) and Q. Due to the simplicity of the spatial domain for the present problem, no spatial discretizations are employed. Hence, functions employed in construction of the responses are defined in the whole problem domain. For more complex spatial domains, special techniques like finite element, boundary element, or finite differences would have to be employed.It is proposed that approximated solutions to the stochastic displacement response of the beam have the following form:(30)
u
(
x
,
ξ
(
ω
)
)
=
∑
i
=
1
∞
u
i
δ
i
(
x
,
ξ
(
ω
)
)
,
where u
i
∈
R and ∀
i
∈
N are coefficients to be determined and δ
i
∈
V are the test functions. Numerical solutions to the variational problem defined in (12) will be obtained. Hence, it becomes necessary to define spaces less abstract than those defined earlier, but without compromising the existence and uniqueness of the solution. Consider two complete orthogonal systems Φ
=
span
⁡
〈
{
ϕ
i
∈
C
2
(
D
,
R
)
∣
ϕ
i
(
x
)
=
0
∧
Δ
ϕ
i
(
x
) = 0
,
in
x
∈
Γ
and
∀
i
∈
N
}
〉 and Ψ
=
span
⁡
〈
{
ψ
i
}
i
=
1
∞
〉, such that Φ
¯
Q
=
Q, Ψ
¯
L
2
(
Ω
,
F
,
P
)
=
L
2
(
Ω
,
F
,
P
), [29], and define the tensor product between Φ and S as
(31)
(
ϕ
⊗
ψ
)
i
(
x
,
ξ
(
ω
)
)
=
ϕ
j
(
x
)
·
ψ
k
(
ξ
(
ω
)
)
,
000000000000000
with
(
j
,
k
)
∈
N
×
N
.
To simplify the notation, we will use δ
i
=
(
ϕ
⊗
ψ
)
i. Since approximated numerical solutions are derived in this paper, the solution space has finite dimensions. This implies truncation of the complete orthogonal systems Φ and Ψ. Hence one has
(32)
Φ
m
=
span
⁡
〈
{
ϕ
i
∈
C
2
(
D
,
R
)
∣
ϕ
i
(
x
)
=
0
∧
Δ
ϕ
i
(
x
)
=
0
,
in
x
∈
Γ
,
∀
i
∈
{
1
,
…
,
m
}
ϕ
i
}

{
ϕ
i
∈
C
2
(
D
,
R
)
〉
,
and Ψ
n
=
span
⁡
〈
{
ψ
i
}
i
=
1
n
〉, which results in V
M
=
Φ
m
⊗
Ψ
n. The dimension of Ψ
n
(
n
=
dim
⁡
〈
Ψ
n
〉
) depends on the dimension of the random variable vector ξ
(
ω
) and on the order of chaos polynomials. Let “s” be the dimension of random vector ξ
(
ω
) and “p” the order of chaos polynomials; then the dimension of Ψ
n is given by
(33)
n
=
(
s
+
p
)
!
s
!
·
p
!
.
Since dim
⁡
〈
Φ
n
〉
<
∞ and dim
⁡
〈
Ψ
n
〉
<
∞ (following [47]), one has that the dimension of the approximation space, V
M, is given by
(34)
M
=
dim
⁡
(
V
M
)
=
dim
⁡
(
Φ
m
⊗
Ψ
n
)
=
dim
⁡
(
Φ
m
)
·
dim
⁡
(
Ψ
n
)
=
m
·
n
.With the above definitions and results (31), it is proposed that numerical solutions are obtained from truncation of the series expressed in (30) at the Mth term:
(35)
u
M
(
x
,
ξ
(
ω
)
)
=
∑
i
=
1
M
u
i
δ
i
(
x
,
ξ
(
ω
)
)
=
∑
i
=
1
M
u
i
(
ϕ
⊗
ψ
)
i
(
x
,
ξ
(
ω
)
)
.
Substituting (35) in (12), one arrives at the approximated variational problem:
(36)
Find
{
u
i
}
i
=
1
M
∈
R
M
such
that
,
∑
i
=
1
M
a
(
δ
i
,
δ
j
)
u
i
=
f
(
δ
j
)
,
∀
δ
j
∈
V
M
,
where
(37)
a
(
δ
i
,
δ
j
)
=
∫
Ω
∫
D
(
α
·
Δ
ω
δ
i
·
Δ
ω
δ
i
)
(
x
,
ω
)
d
x
d
P
(
ω
)
;
f
(
δ
j
)
=
∫
Ω
∫
D
(
q
·
δ
j
)
(
x
,
ω
)
d
x
d
P
(
ω
)
.
The approximated variational problem (36) consists in finding the coefficients of the linear combination expressed in (35). Using a vector-matrix representation, the system of linear algebraic equations defined in (36) can be written as
(38)
K
U
=
F
,
where K
∈
M
M
(
R
) is the stiffness matrix, U
=
{
u
i
}
i
=
1
M is the displacement vector, and F
=
{
f
i
}
i
=
1
M is the loading vector. Elements of the stiffness matrix are defined as
(39)
K
=
[
k
i
j
]
M
×
M
,
k
i
j
=
∫
Ω
∫
D
(
α
·
Δ
ω
δ
i
·
Δ
ω
δ
j
)
(
x
,
ω
)
d
x
d
P
(
ω
)
.
The load vector is given by
(40)
F
=
{
f
j
}
j
=
1
M
,
f
j
=
∫
Ω
∫
D
(
q
·
δ
j
)
(
x
,
ω
)
d
x
d
P
(
ω
)
.
The sparseness of the stiffness matrixes for Example 1 (to be presented) is shown in Figure 1.Sparsity pattern of stiffness matrix of Example 1, (a) form = 2, n = 5, and p = 1; (b) for m = 2, n = 35, and p = 3.
(a)
(b)Remember that “p” is the degree of the chaos polynomial. The stiffness matrix corresponding to Figure 1(a) has dimension 10 and n
z = 42 (number of nonzero elements), whereas the matrix corresponding to Figure 1(b) has dimensions 70 and n
z
=
1344. The conditioning number (n
C) corresponding to these two matrixes is n
C = 19.70 and n
C = 25.79, respectively. These numbers show that the conditioning number increases with the dimension of the approximation space (
V
M
=
Φ
m
⊗
Ψ
n
). This can lead to ill-conditioning of the of the stiffness matrix and hence to the loss of accuracy of the approximated solution. For Example 2—to be presented—the sparsity pattern is similar. However, for Example 1, case (b) presents a different pattern, with n
z
=
100 and n
C
=
32.82 (for p = 2) and n
z
=
4592 and n
C
=
71.74 (for p
=
3). In this paper, uncertainties in beam and foundation stiffness coefficients are modeled using Legendre polynomials to construct space Ψ
n, defined in the variables ξ
(
ω
)
=
{
ξ
1
(
ω
)
,
ξ
2
(
ω
)
,
ξ
3
(
ω
)
,
ξ
4
(
ω
)
}
∈
Γ
1
×
Γ
2
×
Γ
3
×
Γ
4
=
[
-
1,1
]
4.


## 5. Statistical Moments and Estimates

In this section, the evaluation of first and second order moments, from the approximated Galerkin solution, is presented. In an evaluation of the performance of the developed Galerkin scheme, these moments are compared with the same moments computed through Monte Carlo simulation. In order to reduce spurious correlations between the samples, Latin Hypercube Sampling (LHS) is used in the simulations [48].The statistical moment ofkth order of a random variable u
(
x
,
·
)
∈
L
2
(
Ω
,
F
,
P
), generated by the displacement stochastic process (35) by fixing x
∈
D, taking the kth power of this random variable and integrating with respect to its probability measure,
(41)
μ
u
M
k
(
x
)
=
∫
Ω
u
M
k
(
x
,
ξ
(
ω
)
)
d
P
(
ξ
(
ω
)
)
=
∑
i
1
⋯
∑
i
k
︷
k
times
u
i
1
×
⋯
×
u
i
k
(
ϕ
i
1
×
⋯
×
ϕ
i
k
)
(
x
)
×
∫
Ω
(
ψ
i
1
×
⋯
×
ψ
i
k
)
(
ξ
(
ω
)
)
d
P
(
ξ
(
ω
)
)
.
The integration term d
P
(
·
) is a probability measure, defined in (19). From (41), μ
u
M
k
(
·
) is given by
(42)
μ
u
M
k
(
x
)
=
∑
i
1
,
j
1
⋯
∑
i
k
,
j
k
︷
k
times
(
u
i
1
ϕ
i
1
×
⋯
×
u
i
k
ϕ
i
k
)
×
(
x
)
E
[
ψ
i
1
,
…
,
ψ
i
k
]
,
with
(43)
E
[
ψ
i
1
,
…
,
ψ
i
k
]
=
∫
a
1
b
1
⋯
∫
a
N
b
N
(
ψ
j
1
×
⋯
×
ψ
j
k
)
(
ξ
(
ω
)
)
ρ
1
(
ξ
1
)
=
=
=
=
×
⋯
×
ρ
N
(
ξ
N
)
×
d
ξ
1
(
ω
)
×
⋯
×
d
ξ
N
(
ω
)
.
The integrals in (43) are called iterated integrals. The first order statistical moment, or expected value, of the stochastic displacement process evaluated at a point x
∈
D is
(44)
μ
u
M
(
x
)
=
∑
i
=
1
m
u
(
i
-
1
)
·
n
+
1
ϕ
i
(
x
)
.
The variance of this displacement is
(45)
σ
u
M
2
(
x
)
=
∑
i
=
1
n
∑
j
=
1
n
∑
k
=
2
m
u
i
k
u
j
k
(
ϕ
i
·
ϕ
j
)
(
x
)
.In this paper, Monte Carlo simulation is used to estimate the statistical moments of the displacement stochastic process and of random variables generated. Monte Carlo simulation consists in generating “N” samples for flexural rigidity coefficient {
α
(
x
,
ξ
1
(
ω
j
)
,
…
,
ξ
n
(
ω
j
)
)
}
j
=
1
N, and from this point to undertake the solution of the problem, which defined below
(46)
(
K
U
)
(
ξ
(
ω
j
)
)
=
F
.
The jth-realization of the system defined in (46) consists in getting the vector U
(
ξ
(
ω
j
)
)
=
{
u
k
j
(
ξ
(
ω
j
)
)
}
k
=
1
m. To that it is necessary to evaluate stiffness matrix K
:
Ω
→
M
m
(
R
) for the jth sample of α
=
α
(
x
,
ξ
(
ω
j
)
), evaluated in ξ
(
ω
j
)
=
{
ξ
i
(
ω
j
)
}
i
=
1
N. Rigidity matrix is set by Galerkin method. For the jth-sample of α
=
α
(
x
,
ξ
(
ω
j
)
), has its inputs defined by
(47)
K
(
ξ
(
ω
j
)
)
=
[
k
p
q
(
ξ
(
ω
j
)
)
]
m
×
m
,
k
p
q
(
ξ
(
ω
j
)
)
=
∫
D
α
(
x
,
ξ
(
ω
j
)
)
·
(
Δ
ϕ
p
·
Δ
ϕ
q
)
(
x
,
ω
j
)
d
x
.
It is important to mention that rigidity matrix is set by the Galerkin method in the sense of a deterministic problem for each sample of the flexural rigidity coefficient. By solving (46) the jth-realization of u
(
x
,
ξ
(
ω
j
)
) is given by
(48)
u
(
x
,
ξ
(
ω
j
)
)
=
∑
k
=
1
m
u
k
j
(
ξ
(
ω
j
)
)
ϕ
k
(
x
)
=
U
t
(
ξ
(
ω
j
)
)
Φ
(
x
)
,
where Φ
:
D
→
R
m is a vector-valued function whose components {
ϕ
k
}
k
=
1
m are functions. The Monte Carlo estimates for expected value and variance are
(49)
μ
⌢
u
(
x
)
=
1
N
∑
i
=
1
N
u
(
x
,
ξ
(
ω
j
)
)
;
σ
⌢
u
2
(
x
)
=
1
N
-
1
∑
j
=
1
N
[
u
(
x
,
ξ
(
ω
j
)
)
-
μ
⌢
u
(
x
)
]
2
,
where u
(
x
,
ξ
(
ω
j
)
) is the jth-realization of the displacement stochastic process u
=
u
(
x
,
ξ
(
ω
)
).In the paper, in presented examples the estimates of accrued probability function of a random variableu
(
z
,
ξ
(
ω
)
) generated when fixed in the displacement stochastic process of z
∈
D position are presented. Thus, the estimate for random variable u
(
z
,
ξ
(
ω
)
) generated through numerical solutions for the displacement stochastic process is obtained by fixing z
∈
D and generating numerical values for the vector of random variables ξ
(
ω
j
)
=
{
ξ
1
(
ω
j
)
,
…
,
ξ
n
(
ω
j
)
} and evaluating for the jth-sample of this vector, values for u
(
z
,
ξ
(
ω
j
)
) in (35),
(50)
u
M
(
z
,
ξ
(
ω
j
)
)
=
∑
i
=
1
M
u
i
δ
i
(
z
,
ξ
(
ω
j
)
)
=
∑
i
=
1
M
u
i
(
ϕ
⊗
ψ
)
i
(
z
,
ξ
(
ω
j
)
)
.
Samples of random vector {
ξ
(
ω
j
)
}
j
=
1
N
=
{
ξ
1
(
ω
j
)
,
…
,
ξ
n
(
ω
j
)
}
j
=
1
N are obtained from direct Monte Carlo simulation, using Latin Hypercube Sampling [48]. Twenty thousand samples (N = 20,000) are used in the solution of each problem. The estimates for the accrued probability function of random variables u
(
z
,
ξ
(
ω
)
) and u
M
(
z
,
ξ
(
ω
)
) generated by Monte Carlo simulation and by the numerical solutions, (48) and (50), respectively, are given by
(51)
F
u
M
⌢
(
z
)
=
(
1
N
)
∑
i
=
1
N
1
u
M
(
ω
i
)
;
F
u
⌢
(
z
)
=
(
1
N
)
∑
i
=
1
N
1
u
(
ω
i
)
,
whereas 1
u
M
,
1
u
:
Ω
→
{
0,1
} are the characteristic functions of the random variable generated by (50) and (48), respectively, and defined as follows:
(52)
1
u
(
ω
)
=
{
1
,
ω
∈
B
;
0
,
ω
∉
B
;
∧
1
u
M
(
ω
)
=
{
1
,
ω
∈
B
;
0
,
ω
∉
B
,
characteristic function of the set B.


## 6. Numerical Examples

In this section, two numerical examples for Kirchhoff plates bending problem with uncertainty in bending rigidity are presented. Uncertainty is, in the first example, on the thickness of the plate, and in the second the uncertainty is on Poisson’s module. For all mentioned examples, numerical solutions for a plate singly supported with the geometric domainD
=
{
(
x
,
y
)
∈
R
2
∣
0
<
x
<
1
∧
0
<
y
<
1
} are obtained, subjected to a load distributed q
(
x
,
y
)
=
1 kPa·m, ∀
(
x
,
y
)
∈
D. In the examples, the following mean-value functions are considered for the thickness of the plate and Poisson’s coefficient:
(53)
μ
t
(
x
,
y
)
=
1
100
m
,
∀
(
x
,
y
)
∈
D
;
μ
ν
(
x
,
y
)
=
1
3
,
∀
(
x
,
y
)
∈
D
.The performance of numerical solutions obtained via Galerkin method will be evaluated in terms of the approximation of the first and second order statistical moments (expected value and variance) in relation to its respective estimates obtained via Monte Carlo simulation. The first and second order statistical moments of the numerical solution obtained via Galerkin method are compared with the same moments evaluated via Monte Carlo simulation. For the better exposition of the numerical results, in the statistical moments of the displacement stochastic process is evaluated for a restriction in the dominium of the problem,(
x
,
1
/
2
,
ω
)
⊂
D
×
Ω. Comparisons and evaluations made in face of the Monte Carlo simulation in terms of expected value and variance will be carried out in the stochastic process u
(
x
,
1
/
2
,
ω
). In order to evaluate the accuracy of the Galerkin solutions, the functions relative error in mean, ε
μ
u
:
[
0,1
]
→
R
+, and relative error in variance, ε
σ
u
2
:
[
0,1
]
→
R
+, are defined as
(54)
ε
μ
u
(
x
)
=
{
(
100
%
)
×
|
(
1
-
μ
u
μ
⌢
u
)
(
x
,
1
2
)
|
,
∀
x
∈
(
0,1
)
;
0
,
∀
x
∈
{
0,1
}
;
ε
σ
u
2
(
x
)
=
{
(
100
%
)
×
|
(
1
-
σ
u
2
σ
⌢
u
2
)
(
x
,
1
2
)
|
,
∀
x
∈
(
0,1
)
;
0
,
∀
x
∈
{
0,1
}
,
where μ
u and σ
u
2 are the Galerkin-based expected value and variance, respectively, and μ
⌢
u and σ
⌢
u
2 are the Monte Carlo estimates of the same moments. The accuracy of developed Galerkin solutions, in approximating the expected value and variance of the displacement response, is evaluated for numerical solutions using chaos polynomials of order p
∈
{
1,2
,
3
}. Monte Carlo estimates of expected value and variance (
μ
⌢
u
,
σ
⌢
u
2
) are evaluated from (43) and (44), using 20,000 samples. One still compares the estimates for the accrued probability function obtained via Galerkin method and Monte Carlo simulation. Numerical results presented in this paper were obtained in a HP Pavilion personal computer, running a MATLAB computational code.


### 6.1. Example 1: Plate Thickness

In this example, the uncertainty is present in the thickness of plate being modeled through a parameterized stochastic processt
:
D
×
Ω
→
R
+, defined by,
(55)
t
(
x
,
y
,
ξ
(
ω
)
)
=
μ
t
+
3
·
σ
t
[
ξ
1
(
ω
)
cos
⁡
(
x
l
)
+
ξ
2
(
ω
)
sin
(
x
2
l
)
=
μ
t
+
3
·
σ
t
w
w
+
ξ
3
(
ω
)
cos
⁡
(
y
3
l
)
+
ξ
4
(
ω
)
sin
(
y
4
l
)
]
,
where μ
t is the expected value, σ
t is the standard deviation of plate thickness, and ξ
(
ω
)
=
{
ξ
i
(
ω
)
}
i
=
1
4 is a vector of uniform and independent random variables.Figures2(a) and 2(b) show the evolution of the Monte Carlo estimates for expected value and variance, (
μ
⌢
u
(
z
)
,
σ
⌢
u
2
(
z
)
) as the number of samples, for a random variable u
(
z
,
ω
) generated, by fixing the position z
=
(
1
/
2
,
1
/
2
). In Figures 2(a) and 2(b) it can be observed that for the estimates expected value and random variable u
(
z
,
ω
), present small variation for N > 10,000.Convergence of Monte Carlo simulation results for average and variance for a random variableu
(
z
,
ω
).
(a)
(b)Figure3 presents results for the expected value of the displacement stochastic process u
(
x
,
1
/
2
,
ω
), obtained through Monte Carlo simulation, and by means of the Galerkin solutions of order p
∈
{
1,2
,
3
}. The difference between the curves is imperceptible, showing that even a Galerkin solution of order p = 1 is already acceptable.Figure 3
Expected value function.Figure4 shows the relative error function in expected value (49a) for the numerical solutions obtained through Galerkin method using polynomials of order p
∈
{
1,2
,
3
}. Here it can be observed that, although the error is already small for p = 1, Figure 3, the Galerkin solutions approximate themselves to estimate for the expected value of u
(
x
,
1
/
2
,
ω
), while increasing polynomial orders.Figure 4
Relative error function in expected value.Figure5 shows the graphics of variance functions of the u
(
x
,
1
/
2
,
ω
) stochastic process, for the estimate obtained by Monte Carlo simulation and for those obtained from approximated solutions gotten via Galerkin method. The variance function obtained from the numerical solutions uses (45) for the approximation of the displacement stochastic process. One notices that approximations improve their performance from p = 2. In Figure 6 the graphics of the related error functions in variance for the u
(
x
,
1
/
2
,
ω
) stochastic process based in numerical solutions obtained with polynomials chaos with p
∈
{
1,2
,
3
} level.Figure 5
Variance function.Figure 6
Relative error function in variance.Figure7 presents the graphics for the accrued probability function estimate for a random variable u
(
z
,
ω
). This random variable is generated by fixing the position z
=
(
1
/
2
,
1
/
2
)
∈
D. into the displacement stochastic process. The accrued probability function estimates are obtained from numerical solutions via Galerkin method and determined through generation and realization of 20,000 samples of the random variable vector {
ξ
(
ω
j
)
}
j
=
1
N
=
{
ξ
1
(
ω
j
)
,
ξ
2
(
ω
j
)
,
ξ
3
(
ω
j
)
,
ξ
4
(
ω
j
)
}
j
=
1
N in the displacement stochastic process in z
=
(
1
/
2
,
1
/
2
).Figure 7
Graphics of accrued probability functions estimates of the random variableu
(
z
,
ω
).One may notice in Figure7 that, from p = 2, a suitable approximation between the estimates of probability distribution functions via Galerkin method and Monte Carlo simulation.


#### 6.1.1. Summary of Results for Example 1

Table2 summarizes results of expected value, variance, and corresponding relative errors for the random variable u
(
z
,
ω
), obtained by fixing z
=
(
1
/
2
,
1
/
2
) in the random process displacement, for Example 1. Results are presented for approximated solutions with p
∈
{
1,2
,
3
}. Monte Carlo estimates of expected value and variance are
(56)
μ
⌢
u
(
z
)
=
-
0.000217520750865259
m
;
σ
⌢
u
2
(
z
)
=
1.5897011842709
×
1
0
-
9
m
2
.
Comparing the expected value and variance obtained via Galerkin solutions with simulation results, one notes that the approximated solutions are always smaller than the corresponding Monte Carlo estimates.Table 2
Summary of numerical results of Example  1: expected value, variance, relative errors in expected value, and variance evaluated atz
=
(
1
2
,
1
2
).
p
μ
u
(
z
)
[
m
]
σ
u
2
(
z
)
×
1
0
-
9
[
m
2
]
ε
μ
u
(
z
)
[
%
]
ε
σ
u
2
(
z
)
[
%
]
1
−0.000216110689910037
1.42915224414275
0.648242041098541
10.0993156271553
2
−0.000216313424602775
1.56854000498604
0.555039580215695
1.33114194631272
3
−0.000216317483992034
1.57368391783567
0.553173372397411
1.0075646035687One notes in Table2 that the expected value of the random variable u
(
z
,
ω
) presents good approximation from p
=
1. In this case the relative error function in expected value assumes numerical values under 1%. The same behavior is not noticed in approximation of the random variable variance u
(
z
,
ω
). For P
=
1, the relative error function in variance points to a deviation above 10% between the variance obtained through numerical solutions and the estimate via Monte Carlo simulation. In Table 2, one may note that, from p = 1 to p = 3, the relative error function in expected value presented a reduction of 17.2% while the relative error function in variance presents a reduction of 900.2%. This proved that, in the approximation of variance, the chaos polynomials presented a higher convergence rate in relation to approximation for the expected value.


### 6.2. Example 2: Poisson’s Coefficient

In this example, the uncertainty on Poisson’s coefficient for plate bending rigidity is expressed in (2). Again the parameterized stochastic process is used to model the uncertainty for Poisson’s coefficient ν
:
D
×
Ω
→
R
+ being defined as
(57)
ν
(
x
,
y
,
ξ
(
ω
)
)
=
μ
ν
+
3
·
σ
ν
[
ξ
1
(
ω
)
cos
⁡
(
x
l
)
+
ξ
2
(
ω
)
sin
(
x
2
l
)
μ
ν
3
σ
ν
+
3
e
w
+
ξ
3
(
ω
)
cos
⁡
(
y
3
l
)
+
ξ
4
(
ω
)
sin
(
y
4
l
)
]
,
where μ
ν is the expected value and σ
ν is the standard deviation.Figure8 presents the expected value function in the displacement stochastic process u
(
x
,
1
/
2
,
ω
) obtained from numerical solutions via Galerkin method and Monte Carlo simulation. The expected value functions obtained from numerical solutions approximate suitably of the estimate for the expected value function. One notes that with the numerical solutions obtained with polynomials chaos with p = 1 present already good approximation with its respective estimate. Comparing Figures 3 and 8, one notes that the numerical solution present good performance in the approximation of the expected value function.Figure 8
Expected value function.The graphics of relative error functions in expected value for the stochastic processu
(
x
,
1
/
2
,
ω
). These functions are evaluated from the estimate for the expected value via Monte Carlo simulation and with numerical solution obtained through Galerkin method with p
∈
{
1,2
,
3
}. Values assumed by the relative error function in expected value reflects the good approximation of the expected value function obtained via Galerkin method and the estimate of the expected value function via Monte Carlo simulation. Comparing Figures 4 and 9, one verifies that approximations for the expected value through numerical solutions presented better performance for this example.Figure 9
Relative error function in expected value.The graphics of the variance functions of the stochastic process obtained from numerical solutions via Galerkin method is the estimated variance obtained via Monte Carlo simulation are presented in Figure10. One notes that the graphics for the variance functions determined through numerical solutions present values above the estimate via Monte Carlo simulation. This behavior is due to the variance estimate nature.Figure 10
Variance function.Figure11 presents the graphics of the relative error functions in variance of the stochastic process u
(
x
,
1
/
2
,
ω
) for the solution with p
∈
{
1,2
,
3
}. One notices that the graphics of the relative error functions in variance are accrued for different value of “p.” This indicates that suitable approximations for the variance function are obtained from p = 1. Comparing the graphics of relative error function in variance presented in Figures 6 and 11, one verifies that for p = 1 and p = 2 the approximation for variance is better for Example 2.Figure 11
Relative error function in variance.Figure12 presents the graphics of estimates for the accrued probability functions for the random variable u
(
z
,
ω
). An excellent agreement between the estimates of the accrued probability functions can be observed, showing that the developed Galerkin solution is accurate not only for the first and second moments, but also for the whole probability distribution of the response.Figure 12
Graphics for estimates of the accrued probability functions of the random variableu
(
z
,
ω
).


#### 6.2.1. Summary of Results for Example 2

Table3 summarizes results of expected value, variance, and corresponding relative errors for the random variable obtained by fixing z
=
(
1
/
2
,
1
/
2
) in the random process displacement, for Example 2. Results are presented for approximated solutions with p
∈
{
1,2
,
3
}. Monte Carlo estimates of expected value and variance are
(58)
μ
⌢
u
(
z
)
=
-
0.000211203852750318
m
;
σ
⌢
u
2
(
z
)
=
1.3813208499463
×
1
0
-
11
m
2
.
Comparing the estimates for variance, equations (56) and (58), one verifies that greater dispersion in the random variable (
z
,
ω
) was observed for Example 1. Table 3 presents the numerical results for expected value, variance, error function in expected value, and variance for the random variable u
(
z
,
ω
).Table 3
Summary of numerical results of Example  2: expected value, variance, relative errors in expected value, and variance evaluated atz
=
(
1
2
,
1
2
).
p
μ
u
(
z
)
[
m
]
σ
u
2
(
z
)
×
1
0
-
11
[
m
2
]
ε
μ
u
(
z
)
[
%
]
ε
σ
u
2
(
z
)
[
%
]
1
−0.000211202194750868
1.39814560400783
0.000785023297583971
1.21801926483494
2
−0.000211202310577466
1.39813703070284
0.000730182161605587
1.21739860490758
3
−0.000211202310577466
1.39813702629232
0.000730182158833529
1.21739828560978One notices, in Table3, that approximations for the expected value and variance of the random variable presented good results from p = 1. Thus, relative error functions in expected value and in variance present small variation for different values for “p.” One verifies that the values for variance obtained through numerical solutions present higher values in relation to its respective estimate. By comparing the values of relative error functions in expected value presented in Tables 2 and 3, one notices that Example 2 presented the best approximations for expected value, which the best approximations for variance were seen for Example 1. Differently from what was seen in Example 1, reduction of numerical values of the relative error functions in expected value and variance noticed between p = 1 and p = 3 was lower. In this example, for p = 1 and p = 3, one noticed a reduction of 7.5% in assumed values for the relative error function in expected value, while for the relative error function in variance there is a reduction of 0.05%. By comparing these results with those presented in Example 1,one verifies that convergence rates for expected value and variance were lower for Example 2. Particularly in Example 2, convergence rate for the expected value was higher than for variance, a behavior in opposition to what was noticed for Example 1.


## 6.1. Example 1: Plate Thickness

In this example, the uncertainty is present in the thickness of plate being modeled through a parameterized stochastic processt
:
D
×
Ω
→
R
+, defined by,
(55)
t
(
x
,
y
,
ξ
(
ω
)
)
=
μ
t
+
3
·
σ
t
[
ξ
1
(
ω
)
cos
⁡
(
x
l
)
+
ξ
2
(
ω
)
sin
(
x
2
l
)
=
μ
t
+
3
·
σ
t
w
w
+
ξ
3
(
ω
)
cos
⁡
(
y
3
l
)
+
ξ
4
(
ω
)
sin
(
y
4
l
)
]
,
where μ
t is the expected value, σ
t is the standard deviation of plate thickness, and ξ
(
ω
)
=
{
ξ
i
(
ω
)
}
i
=
1
4 is a vector of uniform and independent random variables.Figures2(a) and 2(b) show the evolution of the Monte Carlo estimates for expected value and variance, (
μ
⌢
u
(
z
)
,
σ
⌢
u
2
(
z
)
) as the number of samples, for a random variable u
(
z
,
ω
) generated, by fixing the position z
=
(
1
/
2
,
1
/
2
). In Figures 2(a) and 2(b) it can be observed that for the estimates expected value and random variable u
(
z
,
ω
), present small variation for N > 10,000.Convergence of Monte Carlo simulation results for average and variance for a random variableu
(
z
,
ω
).
(a)
(b)Figure3 presents results for the expected value of the displacement stochastic process u
(
x
,
1
/
2
,
ω
), obtained through Monte Carlo simulation, and by means of the Galerkin solutions of order p
∈
{
1,2
,
3
}. The difference between the curves is imperceptible, showing that even a Galerkin solution of order p = 1 is already acceptable.Figure 3
Expected value function.Figure4 shows the relative error function in expected value (49a) for the numerical solutions obtained through Galerkin method using polynomials of order p
∈
{
1,2
,
3
}. Here it can be observed that, although the error is already small for p = 1, Figure 3, the Galerkin solutions approximate themselves to estimate for the expected value of u
(
x
,
1
/
2
,
ω
), while increasing polynomial orders.Figure 4
Relative error function in expected value.Figure5 shows the graphics of variance functions of the u
(
x
,
1
/
2
,
ω
) stochastic process, for the estimate obtained by Monte Carlo simulation and for those obtained from approximated solutions gotten via Galerkin method. The variance function obtained from the numerical solutions uses (45) for the approximation of the displacement stochastic process. One notices that approximations improve their performance from p = 2. In Figure 6 the graphics of the related error functions in variance for the u
(
x
,
1
/
2
,
ω
) stochastic process based in numerical solutions obtained with polynomials chaos with p
∈
{
1,2
,
3
} level.Figure 5
Variance function.Figure 6
Relative error function in variance.Figure7 presents the graphics for the accrued probability function estimate for a random variable u
(
z
,
ω
). This random variable is generated by fixing the position z
=
(
1
/
2
,
1
/
2
)
∈
D. into the displacement stochastic process. The accrued probability function estimates are obtained from numerical solutions via Galerkin method and determined through generation and realization of 20,000 samples of the random variable vector {
ξ
(
ω
j
)
}
j
=
1
N
=
{
ξ
1
(
ω
j
)
,
ξ
2
(
ω
j
)
,
ξ
3
(
ω
j
)
,
ξ
4
(
ω
j
)
}
j
=
1
N in the displacement stochastic process in z
=
(
1
/
2
,
1
/
2
).Figure 7
Graphics of accrued probability functions estimates of the random variableu
(
z
,
ω
).One may notice in Figure7 that, from p = 2, a suitable approximation between the estimates of probability distribution functions via Galerkin method and Monte Carlo simulation.


### 6.1.1. Summary of Results for Example 1

Table2 summarizes results of expected value, variance, and corresponding relative errors for the random variable u
(
z
,
ω
), obtained by fixing z
=
(
1
/
2
,
1
/
2
) in the random process displacement, for Example 1. Results are presented for approximated solutions with p
∈
{
1,2
,
3
}. Monte Carlo estimates of expected value and variance are
(56)
μ
⌢
u
(
z
)
=
-
0.000217520750865259
m
;
σ
⌢
u
2
(
z
)
=
1.5897011842709
×
1
0
-
9
m
2
.
Comparing the expected value and variance obtained via Galerkin solutions with simulation results, one notes that the approximated solutions are always smaller than the corresponding Monte Carlo estimates.Table 2
Summary of numerical results of Example  1: expected value, variance, relative errors in expected value, and variance evaluated atz
=
(
1
2
,
1
2
).
p
μ
u
(
z
)
[
m
]
σ
u
2
(
z
)
×
1
0
-
9
[
m
2
]
ε
μ
u
(
z
)
[
%
]
ε
σ
u
2
(
z
)
[
%
]
1
−0.000216110689910037
1.42915224414275
0.648242041098541
10.0993156271553
2
−0.000216313424602775
1.56854000498604
0.555039580215695
1.33114194631272
3
−0.000216317483992034
1.57368391783567
0.553173372397411
1.0075646035687One notes in Table2 that the expected value of the random variable u
(
z
,
ω
) presents good approximation from p
=
1. In this case the relative error function in expected value assumes numerical values under 1%. The same behavior is not noticed in approximation of the random variable variance u
(
z
,
ω
). For P
=
1, the relative error function in variance points to a deviation above 10% between the variance obtained through numerical solutions and the estimate via Monte Carlo simulation. In Table 2, one may note that, from p = 1 to p = 3, the relative error function in expected value presented a reduction of 17.2% while the relative error function in variance presents a reduction of 900.2%. This proved that, in the approximation of variance, the chaos polynomials presented a higher convergence rate in relation to approximation for the expected value.


## 6.1.1. Summary of Results for Example 1

Table2 summarizes results of expected value, variance, and corresponding relative errors for the random variable u
(
z
,
ω
), obtained by fixing z
=
(
1
/
2
,
1
/
2
) in the random process displacement, for Example 1. Results are presented for approximated solutions with p
∈
{
1,2
,
3
}. Monte Carlo estimates of expected value and variance are
(56)
μ
⌢
u
(
z
)
=
-
0.000217520750865259
m
;
σ
⌢
u
2
(
z
)
=
1.5897011842709
×
1
0
-
9
m
2
.
Comparing the expected value and variance obtained via Galerkin solutions with simulation results, one notes that the approximated solutions are always smaller than the corresponding Monte Carlo estimates.Table 2
Summary of numerical results of Example  1: expected value, variance, relative errors in expected value, and variance evaluated atz
=
(
1
2
,
1
2
).
p
μ
u
(
z
)
[
m
]
σ
u
2
(
z
)
×
1
0
-
9
[
m
2
]
ε
μ
u
(
z
)
[
%
]
ε
σ
u
2
(
z
)
[
%
]
1
−0.000216110689910037
1.42915224414275
0.648242041098541
10.0993156271553
2
−0.000216313424602775
1.56854000498604
0.555039580215695
1.33114194631272
3
−0.000216317483992034
1.57368391783567
0.553173372397411
1.0075646035687One notes in Table2 that the expected value of the random variable u
(
z
,
ω
) presents good approximation from p
=
1. In this case the relative error function in expected value assumes numerical values under 1%. The same behavior is not noticed in approximation of the random variable variance u
(
z
,
ω
). For P
=
1, the relative error function in variance points to a deviation above 10% between the variance obtained through numerical solutions and the estimate via Monte Carlo simulation. In Table 2, one may note that, from p = 1 to p = 3, the relative error function in expected value presented a reduction of 17.2% while the relative error function in variance presents a reduction of 900.2%. This proved that, in the approximation of variance, the chaos polynomials presented a higher convergence rate in relation to approximation for the expected value.


## 6.2. Example 2: Poisson’s Coefficient

In this example, the uncertainty on Poisson’s coefficient for plate bending rigidity is expressed in (2). Again the parameterized stochastic process is used to model the uncertainty for Poisson’s coefficient ν
:
D
×
Ω
→
R
+ being defined as
(57)
ν
(
x
,
y
,
ξ
(
ω
)
)
=
μ
ν
+
3
·
σ
ν
[
ξ
1
(
ω
)
cos
⁡
(
x
l
)
+
ξ
2
(
ω
)
sin
(
x
2
l
)
μ
ν
3
σ
ν
+
3
e
w
+
ξ
3
(
ω
)
cos
⁡
(
y
3
l
)
+
ξ
4
(
ω
)
sin
(
y
4
l
)
]
,
where μ
ν is the expected value and σ
ν is the standard deviation.Figure8 presents the expected value function in the displacement stochastic process u
(
x
,
1
/
2
,
ω
) obtained from numerical solutions via Galerkin method and Monte Carlo simulation. The expected value functions obtained from numerical solutions approximate suitably of the estimate for the expected value function. One notes that with the numerical solutions obtained with polynomials chaos with p = 1 present already good approximation with its respective estimate. Comparing Figures 3 and 8, one notes that the numerical solution present good performance in the approximation of the expected value function.Figure 8
Expected value function.The graphics of relative error functions in expected value for the stochastic processu
(
x
,
1
/
2
,
ω
). These functions are evaluated from the estimate for the expected value via Monte Carlo simulation and with numerical solution obtained through Galerkin method with p
∈
{
1,2
,
3
}. Values assumed by the relative error function in expected value reflects the good approximation of the expected value function obtained via Galerkin method and the estimate of the expected value function via Monte Carlo simulation. Comparing Figures 4 and 9, one verifies that approximations for the expected value through numerical solutions presented better performance for this example.Figure 9
Relative error function in expected value.The graphics of the variance functions of the stochastic process obtained from numerical solutions via Galerkin method is the estimated variance obtained via Monte Carlo simulation are presented in Figure10. One notes that the graphics for the variance functions determined through numerical solutions present values above the estimate via Monte Carlo simulation. This behavior is due to the variance estimate nature.Figure 10
Variance function.Figure11 presents the graphics of the relative error functions in variance of the stochastic process u
(
x
,
1
/
2
,
ω
) for the solution with p
∈
{
1,2
,
3
}. One notices that the graphics of the relative error functions in variance are accrued for different value of “p.” This indicates that suitable approximations for the variance function are obtained from p = 1. Comparing the graphics of relative error function in variance presented in Figures 6 and 11, one verifies that for p = 1 and p = 2 the approximation for variance is better for Example 2.Figure 11
Relative error function in variance.Figure12 presents the graphics of estimates for the accrued probability functions for the random variable u
(
z
,
ω
). An excellent agreement between the estimates of the accrued probability functions can be observed, showing that the developed Galerkin solution is accurate not only for the first and second moments, but also for the whole probability distribution of the response.Figure 12
Graphics for estimates of the accrued probability functions of the random variableu
(
z
,
ω
).


### 6.2.1. Summary of Results for Example 2

Table3 summarizes results of expected value, variance, and corresponding relative errors for the random variable obtained by fixing z
=
(
1
/
2
,
1
/
2
) in the random process displacement, for Example 2. Results are presented for approximated solutions with p
∈
{
1,2
,
3
}. Monte Carlo estimates of expected value and variance are
(58)
μ
⌢
u
(
z
)
=
-
0.000211203852750318
m
;
σ
⌢
u
2
(
z
)
=
1.3813208499463
×
1
0
-
11
m
2
.
Comparing the estimates for variance, equations (56) and (58), one verifies that greater dispersion in the random variable (
z
,
ω
) was observed for Example 1. Table 3 presents the numerical results for expected value, variance, error function in expected value, and variance for the random variable u
(
z
,
ω
).Table 3
Summary of numerical results of Example  2: expected value, variance, relative errors in expected value, and variance evaluated atz
=
(
1
2
,
1
2
).
p
μ
u
(
z
)
[
m
]
σ
u
2
(
z
)
×
1
0
-
11
[
m
2
]
ε
μ
u
(
z
)
[
%
]
ε
σ
u
2
(
z
)
[
%
]
1
−0.000211202194750868
1.39814560400783
0.000785023297583971
1.21801926483494
2
−0.000211202310577466
1.39813703070284
0.000730182161605587
1.21739860490758
3
−0.000211202310577466
1.39813702629232
0.000730182158833529
1.21739828560978One notices, in Table3, that approximations for the expected value and variance of the random variable presented good results from p = 1. Thus, relative error functions in expected value and in variance present small variation for different values for “p.” One verifies that the values for variance obtained through numerical solutions present higher values in relation to its respective estimate. By comparing the values of relative error functions in expected value presented in Tables 2 and 3, one notices that Example 2 presented the best approximations for expected value, which the best approximations for variance were seen for Example 1. Differently from what was seen in Example 1, reduction of numerical values of the relative error functions in expected value and variance noticed between p = 1 and p = 3 was lower. In this example, for p = 1 and p = 3, one noticed a reduction of 7.5% in assumed values for the relative error function in expected value, while for the relative error function in variance there is a reduction of 0.05%. By comparing these results with those presented in Example 1,one verifies that convergence rates for expected value and variance were lower for Example 2. Particularly in Example 2, convergence rate for the expected value was higher than for variance, a behavior in opposition to what was noticed for Example 1.


## 6.2.1. Summary of Results for Example 2

Table3 summarizes results of expected value, variance, and corresponding relative errors for the random variable obtained by fixing z
=
(
1
/
2
,
1
/
2
) in the random process displacement, for Example 2. Results are presented for approximated solutions with p
∈
{
1,2
,
3
}. Monte Carlo estimates of expected value and variance are
(58)
μ
⌢
u
(
z
)
=
-
0.000211203852750318
m
;
σ
⌢
u
2
(
z
)
=
1.3813208499463
×
1
0
-
11
m
2
.
Comparing the estimates for variance, equations (56) and (58), one verifies that greater dispersion in the random variable (
z
,
ω
) was observed for Example 1. Table 3 presents the numerical results for expected value, variance, error function in expected value, and variance for the random variable u
(
z
,
ω
).Table 3
Summary of numerical results of Example  2: expected value, variance, relative errors in expected value, and variance evaluated atz
=
(
1
2
,
1
2
).
p
μ
u
(
z
)
[
m
]
σ
u
2
(
z
)
×
1
0
-
11
[
m
2
]
ε
μ
u
(
z
)
[
%
]
ε
σ
u
2
(
z
)
[
%
]
1
−0.000211202194750868
1.39814560400783
0.000785023297583971
1.21801926483494
2
−0.000211202310577466
1.39813703070284
0.000730182161605587
1.21739860490758
3
−0.000211202310577466
1.39813702629232
0.000730182158833529
1.21739828560978One notices, in Table3, that approximations for the expected value and variance of the random variable presented good results from p = 1. Thus, relative error functions in expected value and in variance present small variation for different values for “p.” One verifies that the values for variance obtained through numerical solutions present higher values in relation to its respective estimate. By comparing the values of relative error functions in expected value presented in Tables 2 and 3, one notices that Example 2 presented the best approximations for expected value, which the best approximations for variance were seen for Example 1. Differently from what was seen in Example 1, reduction of numerical values of the relative error functions in expected value and variance noticed between p = 1 and p = 3 was lower. In this example, for p = 1 and p = 3, one noticed a reduction of 7.5% in assumed values for the relative error function in expected value, while for the relative error function in variance there is a reduction of 0.05%. By comparing these results with those presented in Example 1,one verifies that convergence rates for expected value and variance were lower for Example 2. Particularly in Example 2, convergence rate for the expected value was higher than for variance, a behavior in opposition to what was noticed for Example 1.


## 7. Conclusions

In this paper theoretical and numerical results are presented for the Kirchhoff plate bending problem with uncertainty in plate bending rigidity. The uncertainty on bending rigidity and modeled through a parameterized stochastic process. A theoretical result is presented under the format of a theorem on the existence and uniqueness of theoretical solutions for this problem. This result is based on Lax-Milgram theorem. From this point, one uses Galerkin method to obtain numerical solutions for displacement stochastic process. The space for approximate solutions is constructed through the tensor product between the space generated by chaos generalized polynomials, derived from Askey-Wiener scheme, and a conventional space of continuous functions immersed inV
=
L
2
(
Ω
,
F
,
P
;
Q
). Both spaces have finite dimensions, but, due isomorphism and disunity due to construction via tensor product, one has that for approximate solutions the space is sequentially dense in the space of theoretical solutions. The construction of the approximation space meets the necessary hypothesis for the existence and uniqueness theorem. The methodology proposed for obtaining numerical solutions for the displacement stochastic process is evaluated for the two examples. In both, uncertainty is on the plate bending rigidity. In the first example, one considers that uncertainty is on thickness of the plate, while in the second example, uncertainty is on Poisson’s coefficient. In all examples, the mathematical model used to describe uncertainty is a parameterized stochastic process. To measure the performance of numerical solutions obtained through Galerkin method in view of the Monte Carlo simulation are defined relative error function in expected value and variance. These functions have the estimates for the expected value and variance of a stochastic process in their definition. One compares also the estimates obtained from the numerical solutions and by Monte Carlo simulation for the accrued probability function for a random variable generated by the displacement stochastic process. The estimates for the expected value, variance, and accrued probability function were obtained through realization of 20,000 samples for the parameter that presented uncertainty. Generally, one may verify that the expected value is the better approximated by the numerical solutions that variance. With p = 1 for gPC’s, good approximations for the expected value have been obtained already. Additionally, as one increases the level of chaos generalized polynomials, approximations toward expected value and variance are improved. Another feature to be highlighted is uncertainty propagation through Kirchhoff plates model; observed variances in numerical solution were lower than those used in uncertainties on the thickness of plates and Poisson’s coefficient. It is important to mention that this behavior is intrinsic to Kirchhoff plates model, as it was seen both in variances obtained from numerical solutions and in estimates determined by Monte Carlo simulation. Example 1 presented a higher variance; the gPC’s level increase had greater influence on the approximations for the expected value and variance. For this example, it was verified that relative error functions in expected value and variance assumed numerical values lower than the measure that would increase the gPC’s level. The same behavior was noticed in Example 2, but with lower sensitiveness, that is, relative error function in expected value and variance decreased in their numerical values in lower proportion. By comparing the relative error functions in variance it is verified that Example 2 presented better approximation for variance than Example 1. Reductions of relative error functions in variance related to the level of polynomials chaos were higher than for Example 2. In both examples, the estimate for the accrued probability function obtained from numerical solutions presented good approximation of the estimate obtained by Monte Carlo simulation. Approximation between estimates improved with the increase of gPC’s level. This paper presented a basic theoretical result, and yet to be exploited in many applications that, for simpler they are, in the sense of deterministic analysis are still open when uncertainties are considered in the model. The (developed) Galerkin solution developed herein, using polynomials chaos of the Askey-Wiener scheme, was shown to be an accurate and efficient solution for presented examples for Kirchhoff plates stochastic bending problem. It was verified that increase in gPC’s level improved approximations for the expected value, variance, and the estimate of the density function of accrued probability. This shows that increase of approximation space dimension leads to better approximations for the displacement stochastic process. Therefore, if the conditions for existence and uniqueness of the solution are respected a proper family of chaos polynomials is selected, one may obtain good approximation for the displacement stochastic process with a suitable dimension. Thus, Galerkin method becomes an effective alternative to methods based on simulation.



---

*Source: 604368-2014-05-15.xml*