Chaos–Galerkin solution of stochastic Timoshenko bending problems

Computers and Structures 89 (2011) 599–611

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Chaos–Galerkin solution of stochastic Timoshenko bending problems Cláudio R. Ávila da Silva Jr. a, André Teófilo Beck b,⇑ a b

Department of Mechanical Engineering, Federal University of Technology of Paraná, Avenida Sete de Setembro, 3165, 80230-901 Curitiba, PR, Brazil Department of Structural Engineering, EESC, University of São Paulo, Av. Trabalhador Sancarlense, 400, 13566-590 São Carlos, SP, Brazil

a r t i c l e

i n f o

Article history: Received 23 April 2010 Accepted 10 January 2011

Keywords: Timoshenko beam Galerkin method Wiener–Askey scheme Parameterized random fields Monte Carlo simulation

a b s t r a c t This paper presents an accurate and efficient solution for the random transverse and angular displacement fields of uncertain Timoshenko beams. Approximate, numerical solutions are obtained using the Galerkin method and chaos polynomials. The Chaos–Galerkin scheme is constructed by respecting the theoretical conditions for existence and uniqueness of the solution. Numerical results show fast convergence to the exact solution, at excellent accuracies. The developed Chaos–Galerkin scheme accurately approximates the complete cumulative distribution function of the displacement responses. The Chaos–Galerkin scheme developed herein is a theoretically sound and efficient method for the solution of stochastic problems in engineering.  2011 Elsevier Ltd. All rights reserved.

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 different problems of mechanics. Applications to fields as diverse as soil mechanics [1], solid mechanics [2], diffusion processes [3], tribology [4] and fracture mechanics [5,6] are just a few examples. The analysis of stochastic engineering systems has received new impulse with use of finite element methods to obtain response statistics. In early approaches, finite element solutions were obtained via Monte Carlo simulation: samples of random system response were computed based on samples of random system parameters. Perturbation and Galerkin methods were used in this context [7]. These methods allowed representation of uncertainty in system parameters by means of random fields. At the end of the 80’s, Spanos and Ghanem [8] presented a novel approach for the solution of stochastic mechanics problems, using the finite element method. The space of approximate solutions was built using the finite element method and chaos polynomials. These polynomials form a complete system in L2 ðX; F ; PÞ, where ðX; F ; PÞ is a probability space. New theoretical insights into the problem were given by Babuska et al. [9], who presented a stochastic version of the Lax– Milgram lemma [34]. In that paper, the authors presented limitations to the modeling of uncertainty via Gaussian processes. They ⇑ Corresponding author. Tel.: +55 16 3373 9460; fax: +55 16 3373 9482. E-mail addresses: [email protected] (C.R. Ávila da Silva), [email protected] (A.T. Beck). 0045-7949/$ - see front matter  2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2011.01.002

showed that for certain problems of mechanics, use of Gaussian processes can lead to loss of coercivity of the bi-linear form associated to the stochastic problem. Moreover, Babuska and Chatzipantelidis [10] present theoretical results which relate the need of limited variance of random parameters to existence and uniqueness of the solution. Difficulties were indeed encountered in the study of Silva Jr. [11], and resulted in non-convergence of the solution for the bending of plates with random parameters. The lack of convergence was due to the choice of a Gaussian process to represent the uncertainty in strictly positive parameters of the system. The lack of convergence also affects solutions based on perturbation or simulation methods. Despite being an elementary fact (that one cannot use Gaussian processes to represent the uncertainty in strictly positive system parameters), this has been overlooked by several research papers published in the nineties and even recently. Just to provide some examples, Liu and Liu [12] studied the spectral response of concrete structures with uncertainties in material properties and ambient temperature modeled as Gaussian processes; Ghanem and Brzkala [1] used the Galerkin method in a problem of uncertain interface between two continuous media modeled as Gaussian; Anders and Hori [13] solved non-linear problems in mechanics, modeling uncertainty in material properties using Gaussian processes; Charmpis [14] solved stochastic shell problems modeling uncertain Young’s modulus and shell thickness as Gaussian processes; finally Elman and Furnival [3] solved diffusion problems, using Gaussian processes to represent the uncertainty in diffusion coefficients. In many cases, the authors acknowledge the flaw of using Gaussian processes to model strictly positive properties, justifying it by different heuristics and by use of small variances for these processes. Indeed, in these flawed applications,

600

C.R. Ávila da Silva Jr., A.T. Beck / Computers and Structures 89 (2011) 599–611

coefficients of variation of the Gaussian processes, in numerical examples, are always kept very small, in order to avoid convergence problems. These justifications, however, do not overcome the theoretical results about existence and uniqueness of the solution, as argued in [9,10]. In contrast, in the present paper the Lax– Milgram lemma is used to state the conditions for existence and uniqueness of the solution to bending of Timoshenko beams with uncertain material and geometrical properties. The theoretical and numerical solution spaces are constructed by respecting these conditions. In a paper by Xiu and Karniadakis [15] the Wiener–Askey scheme was presented. This scheme represents a family of polynomials which generate dense probability spaces with limited and unlimited support. These polynomial families were called Generalized Chaos Polynomials. The scheme increases the possibilities for modeling of uncertain system parameters. In recent years, large efforts are being made at representing uncertainty in stochastic systems via non-Gaussian processes. The stochastic beam bending problem has been studied by several authors. Vanmarcke and Grigoriu [16] studied the bending of Timoshenko beams with random shear modulus. Elishakoff et al. [17] employed the theory of mean square calculus to construct a solution to the boundary value problem of bending with stochastic bending modulus. Ghanem and Spanos [18] used the Galerkin method and the Karhunen–Loeve series to represent uncertainty in the bending modulus by means of a Gaussian process. Chakraborty and Sarkar [19] 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. Chakraborty and Bhattacharyya [20] solved beam and spatial frame problems, modeling the uncertainty in Young’s modulus as three-dimensional random fields. Impollonia and Elishakoff [21] employed mean square calculus to obtain classical solutions for the covariance of transverse displacements of Timoshenko beams with uncertain shear modulus. Silva Jr. et al. [22] used the Galerkin method and generalized chaos polynomials to obtain numerical solutions for the displacement response of Euler–Bernoulli beams with uncertain material and geometrical properties. Singh and Zheng [23] obtained the response of prismatic beams with random material properties under different deterministic loads and boundary conditions using linear finite element analysis. Kamin´ski [24] employed a perturbation-based stochastic finite element solution using polynomial response functions in static analysis of Timoshenko beams. Although presenting numerical solutions to stochastic beam bending problems, none of the papers just referenced explicitly addresses matters of existence and uniqueness of the solutions. In particular, there is no record of the use of Galerkin method and generalized chaos polynomials in the construction of approximated solutions to Timoshenko beam bending problems. In the present paper, the Galerkin method is used to obtain approximate solutions for the bending of Timoshenko beams with random parameters. Using the Lax–Milgram lemma [9], a brief but novel study about the existence and uniqueness of theoretical solutions is presented, for a system of ordinary, second-order differential equations with random coefficients. A Galerkin solution is constructed by respecting the conditions for existence and uniqueness of the theoretical solution. Uncertainties in material and geometrical properties are represented by parameterized random fields, indexed in uniform random variables. Generalized Chaos polynomials of the Wiener–Askey scheme are used to construct the approximate solution space. The Timoshenko bending problem is formulated in Section 2. Representation of uncertain parameters via parameterized random fields is presented in Section 3. Section 4 shows construction of the approximate Galerkin solution. In Section 5 it is shown how the

first and second order moments of the response are computed from the approximate Galerkin solution. The Galerkin scheme developed herein is applied to two example problems in Section 6. Some conclusions are presented in Section 7. 2. Bending of Timoshenko beams In this section, the strong and weak formulations of the problem of stochastic bending of Timoshenko beams are presented. The Lax–Milgram lemma is used to present a proof of existence and uniqueness of the solution. The strong form of the stochastic Timoshenko beam bending problem is given as,

8   d > ða d/ Þ þ b dw  / ¼ 0; > dx dx dx > >     wð0; x Þ ¼ wðl; xÞ ¼ 0; > > > : /ð0; xÞ ¼ /ðl; xÞ ¼ 0; 8x 2 X;

ð1Þ

where a = E.I and b = G.A are the bending and shear stiffness, respectively, X is the sample space, w is the transverse displacement field, / is the angular displacement field and f is a load term. In order to guarantee existence and uniqueness of the solution, the following hypotheses are necessary:

( H1 :

 2 Rþ n f0g : Pðfx 2 X : aðx; xÞ 2 ½a; a  ; 8x 2 ½0; lgÞ ¼ 1; 9a; a   8x 2 ½0; lgÞ ¼ 1; 9b; b 2 Rþ n f0g : Pðfx 2 X : bðx; xÞ 2 ½b; b;

H2 : f 2 L2 ðX; F ; P; L2 ð0; lÞÞ: ð2Þ Hypothesis H1 ensures that the beam stiffness modulus is positivedefinite and uniformly limited in probability [9]. Moreover, if there  ; x 2 ½0; lg and/or B ¼ fx 2 X : exists A ¼ fx 2 X : aðx; xÞ R ½a; a  x 2 ½0; lg, then PðAÞ ¼ 0 and/or PðBÞ ¼ 0. Hypothesis bðx; xÞ R ½b; b; H2 ensures that the stochastic load process has finite variance. These hypotheses are necessary for the application of the Lax– Milgram lemma [9], 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 of stochastic Timoshenko beam bending problems is presented. A stochastic version [9] of the Lax–Milgram lemma [34] is used for this purpose. For operators with derivatives of order greater than two, no such study has been found in the literature. Classical results from functional analysis are used in this study [9,25,26]. The study requires definition of stochastic Sobolev spaces, tensorial product and density between distribution spaces and Lp spaces. The study also requires definition of the abstract variational problem associated to the strong form (Eq. (1)). The stochastic Sobolev space where the solution to the stochastic 2 beam bending problem is constructed is V ¼ L2 ðX; F ; P; H10 ð0; lÞ Þ, such that,  ( )  w and / are measurable and  R V ¼ ðw; /Þ : ð0; lÞ  X ! R2  R : 2 2  X kwðxÞkH1 ð0;lÞ dPðxÞ; X k/ðxÞkH1 ð0;lÞ dPðxÞ < þ1 ð3Þ

Expression (3) means that, for x 2 X fixed, ðwð; xÞ; /ð; xÞÞ 2 H10 ð0; lÞ  H10 ð0; lÞ, whereas for x 2 (0, l) fixed, wðx; Þ; /ðx; Þ 2 L2 ðX; F ; PÞ. Defining the tensorial product between v 2 L2 ðX; F ; PÞ and ðh; #Þ 2 H10 ð0; lÞ  H10 ð0; lÞ as (w, /) = (v.h, v.#), one has, for fixed x 2 X,

(

wð; xÞ ¼ v ðÞ:hðxÞ 2 H10 ð0; lÞ; /ð; xÞ ¼ v ðÞ:#ðxÞ 2 H10 ð0; lÞ;

601

C.R. Ávila da Silva Jr., A.T. Beck / Computers and Structures 89 (2011) 599–611

whereas for fixed x 2 (0, l),

(

From hypothesis H1, which state the continuity and coervicity of the bilinear form, and from the Lax–Milgram lemma [9], it is guaranteed that the problem defined in Eq. (7) has unique solution and continuous dependency on the data (a, b, f) [9,28]. This yields theoretical support for the developments to follow. For the following numerical computations, it is more appropriate to express the abstract variational problem as,

2

wðx; Þ ¼ v ðxÞ:hðÞ 2 L ðX; F ; PÞ; /ðx; Þ ¼ v ðxÞ:#ðÞ 2 L2 ðX; F ; PÞ:

Hence, one has

 2  2 L2 ðX; F ; P; H10 ð0; lÞ Þ ’ L2 ðX; F ; PÞ  H10 ð0; lÞ  2 ) V ’ L2 ðX; F ; PÞ  H10 ð0; lÞ : It is also necessary to redefine the differential operator for the space obtained via tensorial product. The operator Dmx : V ! L2 ðX; F ; PÞ ðL2 ð0; lÞÞ2 , acts over an element w 2 V in the following way [27],

 m  d v m ðxÞ:hðxÞ; dx

Dmx w :

ð4Þ

where m 2 N and m 6 2. Let u, v 2 V, with u = (w, /) and v = (h, t). The inner product ð; ÞV : V  V ! R is defined as:

ðu; v ÞV ¼

Z Z

l

½ðDx /:Dx tÞ þ ðDx w  /Þ:ðDx h  tÞðx; xÞdxdPðxÞ:

0

X

ð5Þ 1=2 ÞV .

This inner product induces the V norm kukV ¼ ðu; v The bilinear form a : V  V ! R, associated to the problem defined in Eq. (1), is defined as:

aðu; v Þ ¼

Z Z

l 0

X

½ða:Dx /:Dx tÞ þ b:ðDx w  /Þ:ðDhx  tÞðx; xÞdxdPðxÞ: ð6Þ

Find u 2 V such that aðu; v Þ ¼ ‘ðv Þ;

ð7Þ

8v 2 V;

where ‘ : V ! R is a linear functional given by,

‘ðv Þ ¼

Z Z X

l

ðf :v Þðx; xÞdxdPðxÞ:

ð8Þ

0

From hypothesis H1 (Eq. (2)), one can show that the bilinear form has the following properties: a. continuity jaðu; v Þj 6

Z Z X

l

  jðDx /:Dx tÞðx; xÞj þ b:jððD ½a x w  /Þ

0

þ

X

l 0

This alternative formulation of the abstract variational problem does not eliminate the coupling between the fields w and /. 3. Uncertainty representation via Wiener–Askey scheme In order to apply Galerkin’s method, an explicit representation of the uncertainty is necessary. In this paper, uncertainties on Young’s modulus and on beam cross-section height are modeled via parameterized random fields. These are defined as a linear combination of deterministic functions and random variables [29],

jðx; xÞ ¼ lj ðxÞ þ

N X

g i ðxÞni ðxÞ;

ð12Þ

P n ðHÞ ¼

9 #1=2 "Z #1=2 = l jðDx w  /Þðx; xÞj2 dx : jðDx h  tÞðx; xÞj2 dx dPðxÞ ; 0

h 6 C: kDx /kL2 ðX;F ;PÞL2 ð0;lÞ :kDx tkL2 ðX;F ;PÞL2 ð0;lÞ

i

þ Dx w  /kL2 ðX;F ;PÞL2 ð0;lÞ :kDx h  tkL2 ðX;F ;PÞL2 ð0;lÞ 6 C:kukV :kv kV ;

Z Z

ð9Þ

H:0: ¼ P 0 ðHÞ; H:n: ¼ P n ðHÞ \ P n1 ðHÞ? ;

ð15Þ

Hence, random variable u can be represented by the series expansion:

l

½ða:Dx /:Dx /Þ þ b:ðDx w  /Þ

uðxÞ ¼

X

uk wk ðnðxÞÞ;

ð16Þ

k2I

0

c:kuk2V ;

ð14Þ 2

where P n is the closure of P n in L ðX; F ; PÞ. For L ðX; RðHÞ; PÞ, the following orthogonal decomposition is admitted [31]: 1

:ðDx w  /Þðx; xÞdxdPðxÞ P where c = min{a, b}.

with

n¼0

:ðDx w  /Þðx; xÞdxdPðxÞ Z Z l ½ðDx /:Dx /Þ þ ðDx w  /Þ Pc X

)

ð13Þ

L2 ðX; RðHÞ; PÞ ¼  H:n: :

0

X

Cðfni gNi¼1 Þ : C is the polynomial of degree 6 n; ni 2 H; 8i ¼ 1; . . . ; N; N < 1

2

  ; bg. where C ¼ maxfa b. coercivity

aðu; uÞ P

where lj() is the expected value of field j (, ), gi 2 C0(0, l) \ C1(0, l), "i 2 {1, . . . , N} and fni gNi¼1 are statistically independent random variables. The Wiener–Askey scheme is a generalization of chaos polynomials, also known as Wiener–chaos. Chaos polynomials were proposed by Wiener [27] to study statistical mechanics of gases. Xiu and Karniadakis [15] extended the studies of Spanos and Ghanem [8] and Ogura [28] for polynomials belonging to the Wiener–Askey scheme [30], for the representation of random fields by orthogonal polynomials. The Cameron–Martin theorem [30] shows that Wiener–Askey polynomials form a base for a dense subspace of second order random variables L2 ðX; F ; PÞ. Let H # L2 ðX; F ; PÞ be a separable Gaussian Hilbert space and H ¼ span½fni g1 i¼1  be an orthonormal base of Gaussian random variables. The Wiener–Askey scheme allows a representation of any second order random variable u 2 L2 ðX; RðHÞ; PÞ, where RðHÞ is the r-algebra generated by H. Let P n ðHÞ be the vector space spanned by all polynomials of order less than n:

(

:ðDx h  tÞÞðx; xÞjdxdPðxÞ 8 " #1=2 "Z #1=2
ð11Þ

i¼1

The abstract variational problem associated to the strong form (Eq. (1)) is defined as follows:



8 Find ðw;/Þ 2 V such that > > > : R R l ða:D /:D tÞðx; xÞdxdPðxÞ ¼ R R l ðb:ðD w/Þ:tÞðx; xÞdxdPðxÞ; 8ðh; tÞ 2 V: x x x X 0 X 0

ð10Þ

where k is a multi-index, I is a set of natural numbers with compact support, fwk gk2I are chaos polynomials and fuk gk2I are coefficients of the linear combination. In this case chaos polynomials are multi-dimensional Hermite polynomials:

602

C.R. Ávila da Silva Jr., A.T. Beck / Computers and Structures 89 (2011) 599–611

wk ðnðxÞÞ ¼

1 Y

hkm ðnm ðxÞÞ;

ð17Þ

m¼1

where hkm ðÞ is a Hermite polynomial defined in random variable nm. The inner product between polynomials wi and wj in L2 ðX; F ; PÞ is defined as

ðwi ; wj ÞL2 ðX;F ;PÞ ¼

Z

ðwi :wj ÞðnðxÞÞdPðxÞ;

ð18Þ

ð23Þ

X

where dP is a probability measure. These polynomials form a complete ortho-normal system with respect to the probability measure, with the following properties:

w0 ¼ 1;

8 M > > > Find fwi ; /i gi¼1 such that > > h i 8 9 > R R > l > > > > ðb:Dx ti :Dx tj Þðx; xÞdxdPðxÞ wi = > R Rl PM < X 0 > > > i¼1 > hR R i > ¼  X 0 ðf :tj Þðx; xÞdxdPðxÞ; < > l : þ X 0 ðb:tj :Dx ti Þðx; xÞdxdPðxÞ /i ; > > > i > PM hR R l > > > > i¼1 X 0 ða:Dx ti :Dx tj þ b:ti :tj Þðx; xÞdxdPðxÞ /i > > > i > > PM hR R l > : ¼ i¼1 X 0 ðb:Dx ti :tj Þðx; xÞdxdPðxÞ wi ; 8tj 2 V M :

ðwi ; wj ÞL2 ðX;F ;PÞ ¼ dij ;

8i; j 2 N:

ð19Þ

The approximated variational problem (Eq. (23)) consists in finding the coefficients of the linear combination expressed in Eq. (22). Using a vector–matrix representation, the system of linear algebraic equations defined in Eq. (23) can be written as

8 Find w; / 2 RM such that > > < Aw þ B/ ¼ F; > > : Cw ¼ D/;

It is important to observe that in Eq. (19) the polynomials are orthogonal with respect to the standard Gaussian density function of vector n. The proposal of the Wiener–Askey scheme is to extend the result presented in Eq. (16) to other types of polynomials. In analogy to Eq. (13), taking P n ðHÞ ¼ span½fwi gNi¼1 , with H a separable Hilbert space of finite variance random variables, one has that S P ¼ n2N P n ðHÞ is a family of polynomials of the Wiener–Askey scheme, generating a complete orthogonal system in L2 ðX; F ; PÞ.

where A; B; C; D 2 MM ðRÞ. The solution to the linear system given in Eq. (24) is:

4. Galerkin method

Elements of these matrices are:

The Galerkin method is used in this paper to build an approximate solution to stochastic beam bending problems. It is proposed that approximated solutions have the following form:

(

wðx; xÞ ¼

P1

/ðx; xÞ ¼

P1

t

i¼1 wi i ðx;

xÞ;

ð20Þ

i¼1 /i ti ðx; xÞ;

2

W

1 Let W = C0(0, l) \ C1(0, l), with W H0 ð0;lÞ  2 1 ¼ ðH0 ð0; lÞÞ . Consider two complete orthogonal systems   span f i g1 and ¼ span fwi g1 with i 2 W; i 2 N i¼1 , i¼1 , 2 L2 ðX;F ;PÞ

ðX;F ;PÞ

¼ L2 ðX; F ; PÞ.

H10 ð0;lÞ

U¼

u

u

8

W

¼ L ðX; F ; PÞ. Now define the tensorial product such that, W between U and W as:

ðu  wÞi ðx; xÞ ¼ uj ðxÞ:wk ðxÞ; with ðj; kÞ 2 N2 :

ð21Þ

To simplify the notation, we will use ti = (u  w)i. Since approximated numerical solutions are constructed, the solution space has finite dimensions. This implies truncation of the complete orthogonal systems U and W. Hence, one has Um ¼ span½fui gm i¼1  and Wn ¼ span½fwi gni¼1 , which results in VM = Um  Wn, where M = dim(VM) = m.n. With the above definitions and results, it is proposed that numerical solutions are obtained from truncation of the series expressed in Eq. (21) at the Mth term:

(

wM ðx; xÞ ¼ /M ðx; xÞ ¼

PM

t xÞ; i¼1 /i ti ðx; xÞ: i¼1 wi i ðx;

PM

/ ¼ ðAC1 D þ BÞ1 F; w ¼ C1 DðAC1 D þ BÞ1 F:

ð25Þ

8 R Rl > A ¼ ½aij MM ; aij ¼ X 0 ða:Dx ti :Dx tj Þðx; xÞdxdPðxÞ; > > > > R Rl > > < B ¼ ½bij MM ; bij ¼ X 0 ðb:tj :Dx ti Þðx; xÞdxdPðxÞ; Rl R > > C ¼ ½cij MM ; cij ¼ X 0 ðb:Dx ti :tj Þðx; xÞdxdPðxÞ; > > > > > R Rl : D ¼ ½dij MM ; dij ¼ X 0 ða:Dx ti :Dx tj þ b:ti :tj Þðx; xÞdxdPðxÞ: ð26Þ

where wi ; /i 2 R; 8i 2 N are coefficients to be determined and ti 2 V are the test functions. Numerical solutions to the variational problem defined in Eq. (11) 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. From the theorem of Cameron–Martin [32] one has PL

(

ð24Þ

ð22Þ

Substituting Eq. (23) in Eq. (11), one arrives at the approximated variational problem:

The load vector is given by

F ¼ ffi gM i¼1 ;

fi ¼ 

Z Z X

l

ðf :ti Þðx; xÞdxdPðxÞ:

ð27Þ

0

Fig. 1 shows the sparsity pattern of matrices A, B, C and D, for the numerical examples 1 and 2 (to be presented), with p = 4. It is observed that the sparsity pattern structure of matrix A varies significantly from example 1 (random elasticity modulus) to example 2 (random beam cross-section height). The same is not true for matrices B, C and D.The conditioning number (nC) for matrix (AC1D + B)1, in example 1, is nC = 9172.70 for p = 1 and nC = 15962.83 with p = 4. In example 2, conditioning numbers for the same matrix are nC = 7029.73 for p = 1 and nC = 8248.15 for p = 4. Hence, it is observed that matrix (AC1D + B)1 is better conditioned for example 2. For example 1, the conditioning numbers are more sensitive to variations in p. For both examples and for the different values of p, the system (AC1D + B) is always positive-definite. 5. Statistical moments Numerical solutions to be obtained are defined in VM. Interest lies in the statistical moments of the stochastic displacement response. In this section, it is shown how the first and second order moments are evaluated from the numerical solution. The statistical moment of kth order of a random variable w(x, ) is obtained, for a fixed point x 2 [0, l], by taking the kth power of the displacement and integrating with respect to it’s probability measure:

603

C.R. Ávila da Silva Jr., A.T. Beck / Computers and Structures 89 (2011) 599–611

Fig. 1. Sparsity pattern of matrices A, B, C and D for examples 1 (left) and 2 (right), for p = 4.

lkwM ðxÞ ¼

Z

6. Numerical examples

wkM ðx; nðxÞÞdPðnðxÞÞ

X

k times

zfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflffl{ X X ¼  wi1      wik ð/i1      /ik ÞðxÞ i1 ik Z  ðwi1      wik ÞðnðxÞÞdPðnðxÞÞ:

ð28Þ

X

The integration term dP() is a probability measure defined as,

dPðnðxÞÞ ¼

N Y

qi ðni Þdni ðxÞ;

ð29Þ

i¼1

where qi : ni ! R is the probability density function of random variable ni. From the measure and integration theory [31], one knows that the probability measure defined in Eq. (28) is the product measure obtained from the product between probability measure spaces associated to the random variables nðxÞ ¼ fni ðxÞgNi¼1 , with ni:X ? [ai, bi]. The probability measure defined in Eq. (29) one has, k times

l

k wM ðxÞ

zfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflffl{ X X ¼  ðwi1 ui1      wik uik ÞðxÞ  wi1 ;    ; wik i i 1

k

ð30Þ with

 wi1 ;    ; wik ¼

Z

b1

a1



Z

bN

aN

ðwi1      wik ÞðnðxÞÞ

 q1 ðn1 Þ      qN ðnN Þ  dn1 ðxÞ      dnN ðxÞ:

ð31Þ

The integrals in Eq. (31) are called iterated integrals. The first order statistical moment, or expected value, of the stochastic displacement field evaluated at a point x 2 [0, l] is:

lwM ðxÞ ¼

m X

wði1Þ:nþ1 ui ðxÞ:

ð32Þ

i¼1

The second-order statistical moment (the variance) of this displacement is

r2wM ðxÞ ¼

n X n X m X i¼1

wik wjk ðui  uj ÞðxÞ:

ð33Þ

j¼1 k¼2

In the following numerical examples, the statistical moments defined in Eqs. (32) and (33) are evaluated and compared with the same moments obtained via Monte Carlo simulation.

In this section, two numerical examples of the stochastic Timoshenko beam bending problem are presented. The purpose of these examples is to show that use of the Galerkin method and of generalized chaos polynomials, originated from the Wiener–Askey scheme, result in an efficient strategy to approximate the numerical solution of stochastic Timoshenko bending problems. In the first example, uncertainty in the Young’s modulus is considered. In the second example, uncertain beam height is considered. In both cases, uncertainty is modeled by parameterized random fields. In both examples, the beam is clamped at both ends,  the span equals half meter l ¼ 12 m , the cross-section is rectangu1 1 lar with b ¼ 25 m and h ¼ 20 m, and the beam is subject to a deterministic, uniformly distributed loading, given by f ðxÞ ¼ 100   KPa:m;8x 2 0; 12 . Fig. 2 shows a sketch of the beam studied in both examples. In the examples presented herein, the choice of representation for parameter uncertainties is guided by mathematical soundness and simplicity, rather than on fitting to actual, measured data. The performance of Galerkin solutions is evaluated with respect to the ability of this method to represent first and second order moments of angular and transverse displacement responses. Accuracy is measured with respect to results of Monte Carlo simulation. To evaluate the error of approximated solutions, relative error functions in expected value and variance are defined as: 8 _ < el ðxÞ ¼ jðl  l w ÞðxÞj w w  ^ _ : el ðxÞ ¼ jðl  l ÞðxÞj; 8x 2 0; 12 ; / / /

8   _ > < er2 ðxÞ ¼ ðr2w  r2w ÞðxÞ w

 _ > : er2 ðxÞ ¼ jðr2/  r2/ ÞðxÞj; 8x 2 0; 12 ; /

ð34Þ

where (lw, l/) and ðr2w ; r2/ Þ, are the value _Galerkin-based  _ expected  _ _ and variance, respectively, and lw ; l/ and r2w ; r2/ are the Monte Carlo estimates of the same moments. In this paper, a family of Legendre polynomials is used to construct space Wn, defined by two independent, uniform random variables (nrv = 2). Numerical solutions are obtained for m = 3(m = dim(Um)) and for different orders of chaos polynomials, with p 2 {1, 2, 3}. The size of the chaos polynomial basis (Wn) ber v Þ! comes n 2 {3, 6, 10}, since n ¼ ðpþn . This results in numerical solup!nrv ! tions with M 2 {9, 18, 30} coefficients to determine.

604

C.R. Ávila da Silva Jr., A.T. Beck / Computers and Structures 89 (2011) 599–611

f ( x ) = 100 KPa.m

l h

b Fig. 2. (a) Clamped beam studied in examples 1 and 2 (b) Details of the cross-section.

Numerical results presented in this paper were obtained in a personal computer, running a MATLAB computational code. 6.1. Example 1: uncertain Young’s modulus In this example, the materials Young’s modulus is modeled as a parameterized random field,

pffiffiffi 3:rE ðxÞ½n1 ðxÞ cosð2xÞ þ n2 ðxÞ sinð2xÞ;   1 8ðx; xÞ 2 0;  X; 2

Eðx; xÞ ¼ lE ðxÞ þ

ð35Þ

 where lE ðxÞ ¼ 210 GPa; 8x 2 0; 12 is the mean value, rE is the standard deviation and {n1, n2} are orthogonal uniform random  1  variables. Numerical solutions are obtained for rE ðxÞ ¼ 10 :lE ðxÞ;  8x 2 0; 12 . The covariance function of Young’s modulus is shown in Fig. 3, where it can be observed that Young’s modulus is a strictly stationary random field. The covariance function shown in Fig. 3 is obtained in exact form from Eq. (35) and from the orthogonality property of random variables {n1, n2}. Results obtained via Monte Carlo simulation are used as reference and are shown first. Fig. 4 shows all the ten thousand samples or realizations (NS = 10.000) of the transverse and angular displacements of the beam. Fig. 5 shows the expected value (left) and absolute error in expected value (right) of random transverse displacements obtained for Galerkin solutions with p 2 {1, 2, 3}. It is observed that approximated solutions converge very fast, and monotonically, to the exact solution. A good approximation of the expected value is already obtained for p = 1, with a relative error of the order of 0.3%. Similar results are shown in Fig. 6, in terms of the expected value (left) and absolute error in expected value (right) of angular displacements, obtained via Galerkin method with p 2 {1, 2, 3}. It is observed that convergence of Galerkin solutions to the exact

result is monotonous and very fast. Even for p = 1, the absolute error is already very small. Fig. 7 shows the covariance functions of transverse displacements (left) and angular displacement (right) random fields. There are no observable differences between the Galerkin and Monte Carlo estimates of these functions. Convergence of Galerkin solutions in variance of transverse displacements is presented in Fig. 8. The variance obtained via simulation and via Galerkin method, for p 2 {1, 2, 3}, is shown left. The absolute error in variance is shown right. It is observed that convergence is fast and monotonous, and that a good approximation of the variance is already obtained for p = 2. Similar results are presented in Fig. 9, for the variance of angular displacements. Fig. 9 shows the variance (left) and absolute error in variance (right) obtained via Galerkin method for p 2 {1, 2, 3}. It is observed that convergence of Galerkin solutions to the exact result is monotonous and very fast. For p = 2, the absolute error in variance is already very small. Comparing Figs. 8 and 9 with Figs. 5 and 6, it can be observed that the convergence rate is faster for the expected value than for the variance. For the expected value, a very good approximation is already obtained for p = 1. An approximation for variance, with the same level of accuracy, is only obtained for p P 2. The fast and monotonous convergence rate of Galerkin solutions, as observed, is not limited to the first and second moments of the random displacement fields. Fig. 10 shows that this excellent convergence rate is also obtained for the probability distribution functions of the random displacement responses. Fig. 10 shows histograms (above) and cumulative distribution functions (below)   of the random variables transverse displacement wð1Þ ¼ w 14 ; x 1  4 (left) and angular displacement /ð 1 Þ ¼ / 10 ; x (right), for a Galer10 kin solution with p = 3. It is observed that the Galerkin solutions converge to the exact probability distribution function of the random response. In Fig. 10, Monte Carlo histograms are obtained directly from ten thousand samples of the beams displacement responses. The Galerkin histograms (and cumulative distribution functions) are obtained by sampling the random variables {n1, n2} which form the chaos polynomial basis. Very good agreement between the histograms and cumulative distribution functions can be observed in Fig. 10. Table 1 summarizes results of expected value, variance and corresponding absolute errors for the random variables wð1Þ and /ð 1 Þ . 4 10 Monte Carlo estimates of expected value and variance were obtained as: _

lw

¼ 0:00356222482526572 m;

_

¼ 1:79131108179577  106 m2 ;

ð14Þ

r2w

ð14Þ

_

l/

¼ 0:0214559516927593 rad;

_

¼ 6:69005154692174  105 rad :

ð101 Þ

r2/ Fig. 3. Covariance function of Young’s modulus, example 1.

ð101 Þ

2

C.R. Ávila da Silva Jr., A.T. Beck / Computers and Structures 89 (2011) 599–611

Fig. 4. Ten thousand samples of the transverse (left) and angular (right) displacements of the beam.

Fig. 5. Expected value of transverse displacement (left) and absolute error in expected value (right).

Fig. 6. Expected value of angular displacements (left) and absolute error in expected value (right).

605

606

C.R. Ávila da Silva Jr., A.T. Beck / Computers and Structures 89 (2011) 599–611

Fig. 7. Covariance functions of transverse (left) and angular (right) displacements.

Fig. 8. Variance of transverse displacements (left) and absolute error in variance (right).

Fig. 9. Variance of angular displacements (left) and absolute error in variance (right).

Comparing these results, it is possible to observe that the coefficient of variation (c.o.v.) of random variable /ð 1 Þ is larger than that of 10 variable wð1Þ . 4

Comparing the expected value and variance obtained via Galerkin solutions with simulation results, one notes that both approximated solutions are always smaller than the corresponding Monte

607

C.R. Ávila da Silva Jr., A.T. Beck / Computers and Structures 89 (2011) 599–611

Fig. 10. Histograms (above) and cumulative distribution functions (below) of the random variables transverse displacement wð1Þ (left) and angular displacement /ð 1 Þ (right). 4

10

Table 1 Summary of numerical results for example 1: expected value, variance, and absolute errors in expected value and variance of random variables wðLÞ and /ð LÞ . 5

2

p

wð L Þ 2

L

lwM ð2LÞ

r2wM ðLÞ

elw

1 2 3

0.00350582949397029 0.00355247654498115 0.00355975330839172

1.11836164058010  106 1.61106357505796  106 1.74574892254426  106

5.63953312954367  105 9.74828028457416  106 2.47151687400755  106

p

/ð L Þ

r2/ ðLÞ M 5

el/

2

L 2

6.72949441215671  107 1.80247506737815  107 4.55621592515130  108

5

l/M ð5LÞ 1 2 3

er2w

2

5

4.16924437705643  10 6.01051253833884  105 6.51715708574915  105

0.0211120220213612 0.0213962873837220 0.0214408039100452

Carlo estimates. It can also be observed in the tabulated values that the convergence rate is faster and the absolute errors are smaller for random variable wð1Þ than for /ð 1 Þ . 4

10

6.2. Example 2: random beam cross-section height In this example, random height of the beams cross-section is considered, and modeled as a parameterized random field:

pffiffiffi 3:rh ðxÞ:½n1 ðxÞ cosð2xÞ þ n2 ðxÞ sinð2xÞ;   1 8ðx; xÞ 2 0;  X; 2

hðx; xÞ ¼ lh ðxÞ þ

ð36Þ

 1 where lh ðxÞ ¼ 201 m; is the mean value, rh ðxÞ ¼ 1 8x 2 0; 2 1 l ðxÞ; 8 x 2 0; : is the standard deviation and {n1, n2} are indeh 10 2 pendent and uniformly distributed random variables. As in Example 1, the height is a strictly stationary random field. Results of Monte Carlo simulation are shown first. Fig. 11 shows all 10.000 sampled realizations of the transverse and angular displacements of the beam, for Example 2. Fig. 12 shows the convergence of Monte Carlo simulation results, in terms of the number of samples N, for the mean and variance of random variables wð1Þ and /ð 1 Þ . The figure shows that 4 10 results are stable, and independent of sample size, for NS> 2000. Fig. 13 shows the expected value (left) and absolute error in expected value (right) of random transverse displacements obtained for Galerkin solutions with p 2 {1, 2, 3}. It is observed that approximated solutions converge very fast, and monotonically, to the exact solution. Very good approximation of the expected value is already obtained for p = 1.Comparing Figs. 5 and 13, it is observed that convergence is faster and the error is smaller for the case of

L

er2

5

/

4

3.43929671398107  10 5.96643090372932  105 1.51477827141332  105

L 5

2.52080716986531  105 6.79539008582901  106 1.72894461172592  106

random beam height, in comparison to random Young’s modulus. This behavior can be explained by the way uncertainties propagate through the system of equations. Referring to Eqs. (1), (25) and (26), it can be observed that w 1/E andw (h2 + 1)/(h5 + h3 + h). Hence, since (x2 + 1)/(x5 + x3 + x) 6 1/x for all x > 0, the error is smaller and convergence is faster for beam height uncertainty. This is a characteristic of the Timoshenko beam. For the Euler–Bernoulli beam, the opposite behavior was observed by the authors [22]: the error is smaller and convergence is faster for random Young’s modulus. These differences in uncertainty propagation are further discussed in Ref. [33]. Fig. 14 shows the covariance functions of transverse displacement (left) and angular displacement (right) random fields. There are no observable differences between the Galerkin and Monte Carlo estimates of these functions. Comparing Figs. 7 and 14, it is noted that the variance is smaller for Example 2 than for Example 1, for both displacement fields. Fig. 15 shows convergence in variance of transverse displacements of Galerkin solutions with p 2 {1, 2, 3}. Fig. 15 shows the variance (left) and absolute error in variance (right). It is observed that convergence is fast and monotonous, and that good approximation of the variance is already obtained for p = 2. Fig. 16 shows similar results for the variance (left) and absolute error in variance (right) of angular displacements, obtained via Galerkin method with p 2 {1, 2, 3}. It is observed that convergence of Galerkin solutions to the exact result is monotonous and very fast. For p = 2, the absolute error in variance is already very small. Fig. 17 shows histograms (above) and cumulative distribution functions (below) of the random variables transverse displacement  1  wð1Þ ¼ w 14 ; x (left) and angular displacement /ð 1 Þ ¼ / 10 ;x 4 10 (right), for a Galerkin solution with p = 3. The figure shows that

608

C.R. Ávila da Silva Jr., A.T. Beck / Computers and Structures 89 (2011) 599–611

Fig. 11. Ten thousand samples of the transverse (left) and angular (right) displacements of the beam.

Fig. 12. Convergence of simulation results in mean (left) and variance (right) of random variables wð1Þ (above) and /ð 1 Þ (below). 4

Fig. 13. Expected value of transverse displacement (left) and absolute error in expected value (right).

10

C.R. Ávila da Silva Jr., A.T. Beck / Computers and Structures 89 (2011) 599–611

609

Fig. 14. Covariance functions of transverse (left) and angular (right) displacements.

Fig. 15. Variance of transverse displacements (left) and absolute error in variance (right).

the excellent convergence rate of Galerkin solutions, observed in terms of mean and variance of the displacements fields, is also obtained for the probability distribution functions of these fields. In Fig. 17, the Monte Carlo histograms are obtained directly from

the ten thousand samples of the beams displacement responses. The Galerkin histograms (and cumulative distribution functions) are obtained by sampling random variables {n1, n2} of the chaos polynomial basis.

Fig. 16. Variance of angular displacements (left) and absolute error in variance (right).

610

C.R. Ávila da Silva Jr., A.T. Beck / Computers and Structures 89 (2011) 599–611

Fig. 17. Histograms (above) and cumulative distribution functions (below) of random variables transverse displacement wð1Þ (left) and angular displacement /ð 1 Þ (right). 4 10

Table 2 Summary of numerical results for Example 2: expected value, variance, and absolute errors in expected value and variance of random variables wð LÞ and /ðLÞ . 2 5 p

wð L Þ 2

L

lwM ð2LÞ

r2wM ðLÞ

elw

1 2 3

0.00307737877021987 0.00308451861810926 0.00308468958449653

3.47079383834219  107 2.29898684812457  108 2.41182013409204  108

5.38905094815993  106 1.75079694122874  106 1.92176332849481  106

p

/ð L Þ

r2/ ðLÞ M 5

el/

2

7

2.09700852631451  10 9.59288438400093  107 1.00395733221219  106

0.0183588246035917 0.0184025728073331 0.0184036352449307

Table 2 summarizes results of expected value, variance and corresponding absolute errors for random variables wð1Þ and /ð 1 Þ . 4 10 Monte Carlo estimates of expected value and variance were obtained as: _

lw

ð14Þ

_

r2w

ð14Þ

_

l/

ð101 Þ

_

r2/

¼ 0:00308276782116803 m; ¼ 2:41821247039281  108 m2 ; ¼ 0:0183919751968186 rad; 2

ð101 Þ

L 2

2.07113308655859  108 1.19225622268237  109 6.39233630077205  1011

5

l/M ð5LÞ 1 2 3

er2w

2

¼ 1:01233159433793  106 rad :

Comparing the estimates for expected value and variance of random variables wð1Þ and /ð 1 Þ , obtained for Examples 1 and 2, it can be ob4 10 served that expected value and variance are larger for Example 1. Since the c.o.v. of both input random fields (Young’s modulus and beam height) were the same, it is observed that the Timoshenko beam model is more sensitive to uncertainties in elastic properties than in dimensional cross-section properties. For Euler–Bernoulli beams, an opposite behavior was observed by the authors [22]: Euler–Bernoulli beams are more sensitive to uncertainties in cross-section height than to uncertainties in elastic properties. These observations have led to interesting conclusions about the comparative uncertainty propagation properties of the two beam models [33]. Comparing the tabulated values for expected value and variance obtained via Galerkin solutions, it is noted that convergence is

L

er2

5

/

5

3.31505932268522  10 1.05976105144838  105 1.16600481121047  105

L 5

8.02630741706474  107 5.30431559378329  108 8.37426212573846  109

faster and the absolute error is smaller for random variable wð1Þ 4 than for /ð 1 Þ , as was observed for Example 1. 10

7. Conclusions In this paper, theoretical and practical results for bending of stochastic Timoshenko beams were presented. The Lax–Milgram lemma [9] was used to state the conditions for existence and uniqueness of the solution. By respecting these conditions, a robust Galerkin scheme was developed to obtain approximated solutions to the transverse and angular displacements of Timoshenko beams. The space of approximated solutions was built by tensor product between deterministic functions and generalized chaos polynomials. A family of Legendre polynomials, derived from the Wiener– Askey scheme, was used to construct the solution space. The uncertainty in input parameters was represented via parameterized random fields, respecting the necessary conditions for existence and uniqueness of the solution. The Galerkin scheme developed herein was applied in the solution of two example problems. It was observed that uncertainty propagation through the system is larger for the case of random Young’s modulus than for random beam height. For the same c.o.v. of the input parameters, the dispersion of displacement results was much larger for random Young’s modulus. Interestingly, an opposite behavior was observed by the authors when studying bending of stochastic Euler–Bernoulli beams. These observations have led to interesting conclusions about the comparative uncertainty propagation properties of Timoshenko and Euler–Bernoulli beam models [33].

C.R. Ávila da Silva Jr., A.T. Beck / Computers and Structures 89 (2011) 599–611

For both examples, the Galerkin scheme presented very fast convergence rates to the first and second order moments of transverse and angular displacement fields. The mean value of both fields was accurately predicted with a chaos basis of polynomials of first order (p = 1). For the variance of both fields, a higher polynomial order was required. However, accurate results where already obtained for p = 2. The convergence rate of the developed scheme is so good, that the complete probability distribution functions of transverse and angular displacement fields was accurately represented by a Galerkin solution using third-order polynomials (p = 3). The scheme developed herein is shown to be a theoretically sound and efficient method for the solution of stochastic problems in engineering. Acknowledgments Sponsorship of this research project by the São Paulo State Foundation for Research – FAPESP (Grant No. 2008/10366-4) and by the National Council for Research and Development – CNPq (Grant No. 301679/2009-6) is greatly acknowledged. The first author also acknowledges the Federal University of Technology of Paraná, Curitiba, for the leave permit to the University of São Paulo. Comments by the anonymous reviewers have greatly improved the article and are therefore specially acknowledged. References [1] Ghanem R, Brzakala W. Stochastic finite-element analysis of soil layers with random interface. J Eng Mech 1996;122(4):361–9. [2] Silva Jr CRA, Beck AT. Bending of stochastic Kirchhoff plates on Winkler foundations via the Galerkin method and the Askey–Wiener scheme. Prob Eng Mech 2010;25:172–82. doi:10.1016/j.probengmech.2009.10.002. [3] Elman H, Furnival D. Solving the stochastic steady state diffusion problem using multigrid. IMA J Numer Anal 2007(27):675–88. [4] Silva Jr CRA, Pintaude G, Al-Qureshi HA, Alves MK. An Application of Mean Square Calculus to Sliding Wear. J Appl Mech 2010;77(2):14–21. [5] Rahman S. Probabilistic fracture mechanics: J-estimation and finite element methods. Eng Fract Mech 2001;68(1):107–25. [6] Zhao JP, Huang WL, Dai SH. The applications of 2D elastic stochastic finite element method in the field of fracture mechanics. Int J Press Vessels Pip 1997;71(2):169–73. [7] Araújo JM, Awruch AM. On stochastic finite elements for structural analysis. Comput Struct 1994;52(3):461–9. [8] Spanos PD, Ghanem R. Stochastic finite element expansion for media random. J Eng Mech 1989;125(1):26–40. [9] Babuska I, Tempone R, Zouraris GE. Solving elliptic boundary value problems with uncertain coefficients by the finite element method: the stochastic formulation. Comput Methods Appl Mech Eng 2005;194(12–16):1251–94. [10] Babuska I, Chatzipantelidis P. On solving elliptic stochastic partial differential equations. Comput Methods Appl Mech Eng 2002;191(37–38):4093–122.

611

[11] Silva Jr., CRA, Application of the Galerkin method to stochastic bending of Kirchhoff Plates. Doctoral Thesis, Department of Mechanical Engineering, Federal University of Santa Catarina, Florianópolis, SC, Brazil (in Portuguese), 2004. [12] Liu N, Liu TG. Spectral stochastic finite element analysis of periodic random thermal creep stress in concrete. Eng Struct 1996;18(9):669–74. [13] Anders M, Hori M. Stochastic finite element method for elasto-plastic body. Int J Numer Methods Eng 1999;46(11):1897–916. [14] Charmpis DC. Incomplete factorization preconditioners for the iterative solution of Stochastic Finite Element equations. Comput Struct 2010;88: 178–88. [15] Xiu D, Karniadakis GE. The Wiener–Askey Polynomial Chaos for Stochastic Differential Equations. SIAM J Sci Comput 2002;24(2):619–44. [16] Vanmarcke EH, Grigoriu M. Stochastic finite element analysis of simple beams. J Eng Mech 1983;109(5):1203–14. [17] Elishakoff I, Ren YJ, Shinozuka M. Some exact solutions for the bending of beams with spatially stochastic stiffness. Int J Solids Struct 1995;32(16): 2315–27. [18] Ghanem R, Spanos PD. Stochastic finite elements: a spectral approach. NY: Dover; 1991. [19] Chakraborty S, Sarkar SK. Analysis of a curved beam on uncertain elastic foundation. Finite Elem Anal Des 2000;36(1):73–82. [20] Chakraborty S, Bhattacharyya B. An efficient 3D stochastic finite element method. Int J Solids Struct 2002;39(9):2465–75. [21] Impollonia N, Elishakoff I. Exact and approximate solutions, and variational principles for stochastic shear beams under deterministic loading. Int J Solids Struct 1998;35(24):3151–64. [22] Silva Jr CRA, Beck AT, Rosa E. Solution of the Stochastic Beam Bending Problem by Galerkin Method and the Askey–Wiener Scheme. Lat Am J Solids Struct 2009;6:51–72. [23] Singh AV, Zheng Y. On finite element analysis of beams with random material properties. Prob Eng Mech 2003;18(4):273–8. [24] Kamin´ski M. Perturbation-based stochastic finite element method using polynomial response function for the elastic beams. Mech Res Commun 2009;36(3):381–90. [25] Matthies HG, Keese A. Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations. Comput Methods Appl Mech Eng 2005;194(12–16):1295–331. [26] Brenner SC, Scott LR. The mathematical theory of finite element methods. Springer-Verlag; 1994. [27] Wiener N. The homogeneous chaos. Am J Math 1938;60:897–936. [28] Ogura H. Orthogonal functionals of the Poisson process. IEEE Trans Inform Theory 1972;18(4):473–81. [29] Grigoriu M. Applied non-gaussian processes: examples, theory, simulation, linear random vibration, and matlab solutions. Prentice Hall; 1995. [30] Askey R, Wilson J. Some Basic Hypergeometric Polynomials that Generalize Jacobi Polynomials. Memoir American mathematical society, vol. 319. Providence, RI: AMS; 1985. [31] Jason S. Gaussian hilbert spaces. Cambridge University Press; 1997. [32] Cameron RH, Martin WT. The orthogonal development of nonlinear functionals in series of Fourier–Hermite functionals. Ann Math 1947;48: 385–92. [33] Beck AT, Silva Jr CRA. Timoshenko versus Euler beam theory: Pitfalls of a deterministic approach. Struct Saf 2011;33:19–25. doi:10.1016/j.strusafe. 2010.04.006. [34] Lax PD, Milgram AN. Parabolic equations: Contributions to the theory of partial differential equations. Annals of mathematical studies, vol. 33. Princeton University Press; 1954. pp. 167–190.