Efficient bounds for the Monte Carlo–Neumann solution of stochastic thermo-elasticity problems

Accepted Manuscript Efficient Bounds for the Monte Carlo – Neumann Solution of Stochastic Thermo-Elasticity Problems Cláudio R. Ávila da Silva Jr., André Teófilo Beck PII: DOI: Reference:

S0020-7683(14)00495-8 http://dx.doi.org/10.1016/j.ijsolstr.2014.12.025 SAS 8606

To appear in:

International Journal of Solids and Structures

Received Date: Revised Date:

7 April 2014 18 December 2014

Please cite this article as: Ávila da Silva, C.R. Jr., Beck, A.T., Efficient Bounds for the Monte Carlo – Neumann Solution of Stochastic Thermo-Elasticity Problems, International Journal of Solids and Structures (2015), doi: http://dx.doi.org/10.1016/j.ijsolstr.2014.12.025

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

1 2

Efficient Bounds for the Monte Carlo – Neumann

3

Solution of Stochastic Thermo-Elasticity Problems

4 5

Cláudio R. Ávila da Silva Jr., [email protected]

6

Federal University of Technology of Parana

7

Av. 7 de setembro, 3165, Curitiba, PR, Brazil.

8 9

André Teófilo Beck,[email protected]

10

Structural Engineering Department, EESC, University of São Paulo

11

Av. Trabalhador Sancarlense, 400, São Carlos, SP, Brazil.

12 13 14

Abstract: The numerical solution of stochastic thermo-elasticity problems can be computationally

15

demanding. In this article, a well-known property of the Neumann series is explored in order to

16

derive lower and upper bounds for expected value and second order moment of the stochastic

17

temperature and displacement responses. Uncertainties in axial stiffness and conductivity are

18

represented as parameterized stochastic processes. Monte Carlo simulation is employed to obtain a

19

few samples of the stochastic temperature and displacement fields, from which lower and upper

20

bounds of expected value and second order moment are computed. The proposed methodology is

21

applied to two linear one-dimensional thermo-elastic example problems. It is shown that accurate

22

and efficient bounds can be obtained, for a proper choice of operator norm, with as few as one or two

23

terms in the Neumann expansion. The Monte Carlo – Neumann bounding scheme proposed herein is

24

shown to be an efficient alternative for the solution of stochastic thermo-elasticity problems.

25 26

keywords: Neumann series; Monte Carlo simulation; thermo-elasticity; stochastic processes.

27 28

1.

INTRODUCTION

29

The last few decades have witnessed tremendous developments in the modeling of mechanical and

30

structural systems, due to advances in computational mechanics. Numerical methods such as finite

31

elements, finite difference, boundary elements, etc., have reached broad acceptability and wide

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

_____________________2

32

coverage of applications. New developments address the solution of complex, non-linear problems.

33

Multi-physics analyses allow investigation of new, unforeseen interaction effects between structures,

34

soils, fluids, thermo-dynamic and electric effects. Significant developments have also been recently

35

achieved in modeling uncertainty propagation through mechanical and other types of systems.

36

The Monte Carlo simulation method remains a popular, yet computationally expensive tool for

37

analyzing uncertainty propagation through mechanical systems. The computational cost of Monte

38

Carlo simulation can easily become prohibitive, for highly non-linear problems and complex

39

geometries. More efficient, intrusive methods have recently been developed, such as the stochastic

40

finite element method (Ghanem & Spanos, 1991) or stochastic Galerkin Method (Avila et al., 2010).

41

Intrusive methods have the inconveniency of requiring full re-programming of conventional finite

42

element software. Hence, non-intrusive Monte Carlo simulation methods remain popular in the

43

solution of stochastic mechanics problems.

44

In linear stochastic mechanics problems, the numerical solution of a differential equation is

45

replaced by the solution of a system of linear algebraic equations (stiffness matrix). In this context,

46

when Monte Carlo simulation is employed, for each system realization, the stiffness matrix needs to

47

be evaluated and inverted. Depending on the dimensions of the linear system, and the required

48

number of samples, this can become computationally intensive. For a linear operator in finite

49

dimensions, A :
n →
n , the inverse can be represented by the Neumann series, composed of

50

operators P :
n →
n related to A .In finite dimensions, linear operators are matrices. The

51

objective of using the Neumann series is to replace the matrix inversions by a truncated series

52

expansion. However, depending on the number of terms in the Neumann series, the number of

53

operations to be performed may become larger than required for the direct linear system solution,

54

when some interactive method is used to compute the inverse. Therefore, in this paper, properties of

55

the Neumann series are exploited in order to derive, formally, lower and upper bounds for each

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

_____________________3

56

realization of the numerical approximation of the stochastic system response. To improve efficiency

57

of the proposed scheme, equivalence between norms of finite dimension operators are considered.

58

First use of the Neumann series to solve stochastic problems in mechanics goes back to the late

59

eighties. Yamazaki (1989) employed the Neumann series to obtain samples of the displacement

60

response, for a plane elasticity problem with random Youngs modulus. Araújo & Awruch (1994)

61

derived the response on non-linear static and dynamic problems, with random mechanical properties.

62

Chakraborty & Dey (1996) obtained expected value and variance of displacement responses

63

considering uncertain geometries and mechanical properties. Chakraborty & Dey (1998) formulated

64

the dynamic bending of curved beams, with uncertain stiffness parameters. Lei & Qiu (2000)

65

employed the Neumann series to study uncertainty propagation in structures using movement

66

equations. Chakraborty & Sarkar (2000) estimated statistical moments of the transversal

67

displacement response of curved beams resting on Winkler foundations. Chakraborty &

68

Bhattacharyya (2002) obtained response statistics for linear tri-dimensional elasticity problems,

69

representing elastic properties as Gaussian processes of the continuous. Li et al. (2006) presented a

70

methodology to obtain solutions of linear systems, based on a discretization of the stochastic

71

problem. Schevenels et al. (2007) proposed a methodology based on Green functions, which was

72

compared to the Neumann series, in a wave propagation problem on random Winkler foundation. In

73

all references above, the Neumann series was employed in a conventional fashion, as an alternative

74

to solve the stochastic problem.

75

The present article advances the state of the art by exploring new convergence properties of the

76

Neumann series and by employing well-known properties (Golub & Van Loan, 2012) of the

77

Neumann series to derive formal lower and upper bounds for the realizations of a stochastic system

78

response.To improve efficiency of the proposed scheme, equivalence between norms of finite

79

dimension operators are considered. The developed scheme is illustrated by means of application to

80

two linear one-dimensional thermo-elastic problems. Uncertainties in membrane stiffness and

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

_____________________4

81

conductivity coefficients are represented as parameterized stochastic processes.

Monte Carlo

82

simulation is employed to obtain a few samples of the stochastic temperature and displacement

83

fields, from which lower and upper bounds of expected value and second order moment are

84

computed.

85 86

2.

LINEAR STOCHASTIC THERMO-ELASTICITY PROBLEM

87

The linear stochastic thermo-elasticity problem to be addressed herein is defined in a probability

88

space ( Ω, F , P ) , where " Ω " is the space of events, " F " is a σ -algebra of events and " P " is a

89

probability measure. The linear stochastic thermo-elasticityproblem in a bar consists in finding the

90

displacement and temperature ( u, T ) responses, such that:

(

91

)

 Find ( u, T ) ∈ L2 Ω, F , P; H 2 ( 0, l ) ∩ H 1 ( 0, l ) 2 , such that: ( ) 0    − d EA du + βT = f , dx  dx  d dT ∀ ( x, ω) ∈ ( 0, l ) × Ω, a. e.;  − dx κ dx = q,  u ( 0, ω) = u ( l , ω) = 0,  ∀ω∈ Ω; T ( 0, ω) = T ( l , ω) = 0,  

( (

)

)

(1)

92

where ( f , q ) are the mechanical and thermal loading terms, respectively; EA is the membrane or “unit

93

length”axial stiffness, hereafter simply referred to as stiffness; and κ and β are the thermal

94

conductivity and linear dilatation coefficients, respectively. The following hypotheses are required to

95

warrant existence and uniqueness of the solution:

96

 ∃ κ, κ ∈
+ , P ( x, ω ) ∈ ( 0, l ) × Ω : κ ( x, ω ) ∈ [ κ, κ ] = 1; H1 :  +  ∃λ , λ ∈
, P ( x, ω ) ∈ ( 0, l ) × Ω : EA ( x, ω ) ∈  λ , λ  = 1;

( (

∞

(

)

)

2

)

H2 : ( q, f ) ∈ L Ω, F , P; ( L2 ( 0, l ) ) .

(2)

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

_____________________5

97

The values ( κ, κ ) and ( λ , λ ) are limits of the coefficients, which can be based on physical and/or

98

mathematical restrictions. Hypothesis H1 ensures that the thermal conductivity and axial stiffness

99

coefficients are positive, and uniformly limited in probability (Babuska et al., 2005). Hypothesis H2

100

ensures that mechanical and thermal loading terms have finite variance. From hypothesis H1 and H2,

101

the Lax-Milgram lemma ensures existence and uniqueness of the solutions, for random samples of

102

thermal conductivity and axial stiffness.

103

Approximated numerical solutions are obtained from the abstract variational problem (AVP)

104

associated to Eq. (1). The AVP is defined in V = H 01 ( 0, l ) × H 01 ( 0, l ) . For the kth thermo-mechanical

105

sample EA ( x, ξ ( ωk ) ) , κ ( x, ζ ( ω′k ) ) , the AVP is given by:

106

For fixed {ξ ( ωk ) , ζ ( ω′k )} , find ( u, T ) ∈ V such that :  a ( u , v ) = c ( T , v ) + l ( v ) ,  ∀ ( v, θ ) ∈ V ; b (T , θ ) = h ( θ ) ,

107

{

}

where a ( ⋅, ⋅) , b ( ⋅, ⋅) are bi-linear forms and l ( ⋅) , h (⋅) and c (T , ⋅) are linear functionals, defined as: L  dv ′ , a u v = ( )  ∫0 EA du dx dx ( x, ζ ( ωk ) ) dx,  L  d (β T ) x, ζ ( ω′ ) .v x dx, c (T , v ) = − ∫ dx ( k ) ( )  0  L l v = f .v x dx, ( ) ∫0 ( )( )   L  h θ = ( )  ∫0 ( h.θ)( x ) dx. 

(

108

(3)

)

(4)

109

Eq. (3) represents a system of variational equations. Eq. (3a) is related to the mechanical

110

problem(Eq. 1a), whose solution, for the kth thermo-mechanical sample, is the kth realization of the

111

displacement process. Eq. (3b) is the variational equation associated to the heat transfer problem (Eq.

112

1b), whose solution is the temperature field.Following the Doob-Dynkin lemma (Rao, 2010), thekth

113

realization

114

(u = u ( x, ω , ω′ ) and T = T ( x, ω ) ) . From the AVP and for the kth realization, numerical solutions

of k

ik

the

displacement k

and

temperature

processes

depend

on

{ωk , ω′k }

:

_____________________6

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

115

are obtained via Galerkin method. From the ensemble of realizations, expected value and second

116

order moment of the approximated numerical solution can be evaluated. Importantly, the numerical

117

solutions for the displacement process depend on the solution for the temperature field.

118 119

3.

120

The Galerkin method is employed herein to obtain approximated numerical solutions to the

121

realizations of displacement and temperature fields. The approximation space is obtained as

122

Vmn = span {{ϕ1 ,… , ϕ m } , {ψ1 , … , ψ n }} , using a subset of V, V = span {φi }i∈

123

generality, the thermo-mechanical loading terms are considered deterministic in this paper, and the

124

uncertainty in conductivity and stiffness coefficients is represented using parameterized stochastic

125

process:

126

GALERKIN METHOD

(

V

) . Without loss of

Mκ   κ ( x, ξ ( ω) ) = µ κ ( x ) + ∑ ϑi ( x ) ξi ( ω);  i =1  M EA  EA ( x, ζ ( ω) ) = µ ( x ) + γ ( x ) ζ ( ω ); ∑ EA i i  i =1

(5)

127

where µ κ ( ⋅) and µ EA ( ⋅ ) are the expected values of conductivity and stiffness coefficients, respectively;

128

ϑi , γ i ∈ C0 ( 0, l ) ∩ C

129

ζ ( ω) = {ζi ( ω)}i =1 are random vectors, whose entries are statistically independent random

130

variables, with the following properties:

2

( 0, l )

are

deterministic

functions.

Mκ

ξ ( ω) = {ξi ( ω)}i=1

and

M EA

131

 E [ ξi ] = 0 ∧ P ( ω∈ Ω : ξi ( ω ) ∈ Γi ) = 1, ∀i ∈ {1,…, M κ } ,   E [ ζ i ] = 0 ∧ P ( ω∈ Ω : ζ i ( ω ) ∈ Π i ) = 1, ∀i ∈ {1,…, M EA} ,

(6)

132

where E [⋅] is the expected value operator. In Eq. (5), Γi and Πi are the images of the random

133

variables, i.e. Γ i = ξ i ( Ω ) and Π i = ζ i ( Ω ) with Γ i = [ ai , bi ] ⊂
, Γ i = bi − ai < ∞, ∀i ∈ {1,…, M κ } , and

134

Π i = [ ci , d i ] ⊂
, Π i = di − ci < ∞ , ∀i ∈ {1,… , M ΕΑ } . The image of random vector is ξ : Ω → Γ ,

_____________________7

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

135

Mκ

M

with Γ ⊂
N , which in terms of {Γi }i=1κ , is given by Γ = ∏ Γi . With a proper choice of sets Γ i and i =1

136

Π i and using the uncertainty representation in Eq. (5), one ensures that the stochastic processes"

137

κ ( ⋅, ⋅) "and" EA (⋅, ⋅ ) "will be positive and uniformly limited in probability, satisfying hypothesis

138

H1.The restriction on these coefficients ensures that there will be no violations of physical principles

139

such as total potential energy; hence warranting existence and uniqueness of the solution for the

140

stochastic system. The interested reader may consult the following references for the details on how

141

existence and uniqueness of the solutions are warranted: (Silva Jr. et al. 2009; Silva Jr. et al. 2010;

142

Beck & Silva Jr. 2011; Silva Jr. & Beck, 2011).

143 144

Following the Doob-Dynkin lemma, the displacement and temperaturefields depend on random variables ( ξ ( ω ) , ζ ( ω ) ) :

(

)

145

T ( x, ω) = T ( x, ξ ( ω ) ) = T x, ξ1 ( ω) ,…, ξ M ( ω) ; κ   u ( x, ω) = u ( x, ζ ( ω ) , ξ ( ω) ) = u x, ζ1 ( ω) ,…, ζ MΕΑ ( ω ) , ξ1 ( ω ) ,…, ξ Mκ ( ω ) .

146

For the kth realization of random variables ζ k ( ω ) , ξ k ( ω ) , the Galerkin method is employed

147

(

(

)

to derive approximated numerical solutions, ( um , Tn ) , to the AVP problem in Eq. (3): m  u x u i ζ k ( ω ) , ξ k ( ω ) ϕi ( x ) , , ζ ω , ξ ω = ( ) ( ) k k  m  i =1  n T x , ξ ( ω ) = Ti ξ k ( ω ) ψ i ( x ) , k  n i =1

) ∑ (

(

148

)

) ∑ (

(

)

(7)

)

n

m

i =1

i =1

149

where {Ti ( ξ k ( ω ) )} and {ui ( ζ k ( ω ) , ξ k ( ω ) )} are the coefficients of the linear combination, to be

150

determined. Eq. (7) represents one realization of the approximated displacement and temperature

151

fields for one realization of random vector ( ξ ( ω ) , ζ ( ω ) ) . For the kth realization, the approximated

152

variational problem is obtained as:

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

_____________________8

 For fixed {ζ k , ξ k } , find ( um ,Tn ) ∈Vmn , such that :   a ( um , v ) = c (Tn , v ) + l ( v ) ,  ∀ ( v, θ) ∈Vmn . b (Tn , θ) = h ( θ ) ,

153

(8)

154

By replacing the discrete solution in Eq. (7) into the AVP (Eq. 3), a linear system of algebraic

155

equations is obtained: n

m

 For fixed {ζ k , ξ k } , find ( U ( ζ k , ξ k ) , T ( ξ k ) ) ∈
×
such that :   R ( ζ k ) U ( ζ k , ξ k ) + DT ( ξ k ) = F;  ( CT )( ξ k ) = Q;

156

(9)

157

where R ( ζ k ) ∈ M m (
) , C ( ξ k ) ∈ M n (
) and D ∈ M m×n (
) are stiffness, conductivity and thermic

158

dilatation matrices, respectively. The entries of these matrices are defined as: L  d ϕi d ϕ j  R ( ζ k ) =  rij ( ζ k )  m×m , rij ( ζ k ) = ∫ EA ( x, ζ ( ωk ) ) dx dx dx; 0  L  d ψi d ψ j C ( ξ k ) = cij ( ξ k )  n×n , cij ( ξ k ) = ∫ κ ( x, ξ ( ωk ) ) dx dx dx; 0  L  dϕ  D =  d ij  , d ij = ∫ βψ i j ( x ) dx.   dx m×n  0

( (

159

(

160

163

)

(10)

)

Solution of the linear system of equations in Eq. (9) is given by: −1

U ( ζ k , ξk ) = ( R ( ζ k ) ) F ( ξk ) = H ( ζ k ) F ( ξk ) ,

161 162

)

(11)

where −1

F ( ξk ) = F − D ( C ( ξk ) ) Q ,

(12) t

164

Vector F ( ξ k ) =  f1 ( ξ k ) ,..., f m ( ξ k )  is the thermo-mechanical loading vector. In this vector the

165

thermo-mechanical coupling is evident. In Eq. (11), matrix H ( ζ k ) =  hij ( ζ k ) 

166

stiffness matrix, R ( ζ k ) ,and U ( ζ k , ξ k ) = u1 ( ζ k , ξ k ) ,…, um ( ζ k , ξ k )  is a vector with entries given as

167

the linear combinations:

t

m×m

is the inverse of the

_____________________9

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

m

ui ( ζ k , ξ k ) = ∑ hij ( ζ k , ξ k ) f j .

168

(13)

j =1

169

From Eqs. (7a) and (13), the approximated numerical solution, for the kth realization of the thermo-

170

mechanical loading, is given by: m

m

um ( x, ζ k , ξ k ) = ∑ ∑ hij ( ζ k ) f j ( x, ξ k ) ϕi ( x ) = F ( ξ k )⋅ ( H ( ξ k ) Φ ( x ) ) .

171

i =1 j =1

(

)

(14)

172

An analogous process can be followed to obtain the approximated temperature solution, for the kth

173

realization of the thermo-mechanical loading.

174

In the solution by direct Monte Carlo simulation, estimates of expected value and second

175

order moment of the responses are obtained from the ensemble of realizations. In Eqs. (13) and (14),

176

one observes that each system response evaluation requires one inversion of the stiffness matrix,

177

R (⋅) . Computation of this inverse should be avoided, due to the large computational cost involved.

178

For solution of the linear system (Eq. 9), iterative methods such as Jacobi, Gauss-Seidell or

179

Conjugated Gradients are generally employed (Kincaid & Cheney, 2002). However, the cost of

180

direct Monte Carlo simulation can still become prohibitive, for instance for non-linear problems

181

involving large number of samples. Under such conditions, one alternative to reduce the

182

computational cost is use of the Neumann series.

183 184

4.

THE NEUMANN SERIES

185

In this paper, the Neumann series is employed for a linear operator defined in a finite dimensional

186

space. In this case, the mathematical object of the operator is a matrix A :
n →
n . For the

187

parameterized uncertainty model employed in this paper (Eq. 4), the stiffness and conductivity

188

matrices, for the kth thermo-mechanical realization EA ( x, ζ k ( ω) ) , κ ( x, ξ k ( ω) ) , can be decomposed

189

as follows:

{

}

_____________________10

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

 R ( ζ k ) = R 0 + ∆R ( ζ k )    C ( ξ k ) = C0 + ∆C ( ξ k ) 

190

⇒

R ( ζ k ) = R 0 ( I + PR ( ζ k ) ) ,  C ( ξ k ) = C 0 ( I + PC ( ξ k ) ) ,

−1

(15)

−1

−1

191

where PR ( ζ k ) = ( R 0 ) ∆R ( ζ k ) and PC ( ζ k ) = ( C0 ) ∆C ( ζ k ) . The entries of matrices ( R 0 ) and

192

( C0 )

193

Eq. (15) in Eq. (11), formally, a solution U = U ( ζ k , ξ k ) is given by:

194

T ( ξ ) = ( I + P ( ξ ) ) −1 T ; k C k 0   −1 −1 −1  U ( ζ k , ξ k ) = ( I + PR ( ζ k ) ) U 0 − ( R 0 ) D ( I + PC ( ξ k ) ) T0 ; 

−1

are obtained as the expected values of the stiffness and conductivity coefficients. Substituting

)

(

−1

(16)

t

−1

t

195

with T0 = ( C0 ) Q , U0 = ( R 0 ) F , and with T0 = T10 ,…, Tn0  and U 0 = u10 , …, um0  . From the

196

Neumann series theorem, if 0 < PC ( ξk ) , PR ( ζ k ) < 1 then ( I + PR ( ζ k ) ) and ( I + PC ( ξk ) ) exist,

197

and can be represented by the following series (Golub & Van Loan, 2012):

198

−1

0 < PC ( ξ k )  0 < PR ( ζ k )

∞ i −1 −1  ∃ I + P ξ ∧ I + P ξ = ( ) ( ) ( ) ( ) ( PC ( ξk ) ) ; ∑ C k C k  < 1   i =1 ⇒  ∞ i −1 −1 < 1  ∃( I + PR ( ξ k ) ) ∧ ( I + PR ( ζ k ) ) = ∑ ( PR ( ζ k ) ) .  i =1

−1

(17)

199

Moreover, the theorem provides an upper bound, for the summation in Neumann series, Eq. (17), of

200

matrices ( I + PC ( ξk ) ) and ( I + PR ( ζ k ) ) . These bounds are given by:

201

202

203

204

−1

−1

−1 1   ( I + PC ( ξ k ) ) ≤ 1 − P ( ξ ) ; C k  .  1  ( I + P ( ζ ) ) −1 ≤ . R k  1 − P ζ ( ) R k 

(18)

From Eqs. (15) and (17), one obtains: ∞ −1 −1 i  −1 C ξ = I + P ξ C = ( ) ( ) ( ) ( ) ( ) ( PC ( ξk ) ) ( C0 )−1 ; ∑ k C k 0   i =1  ∞ −1 −1 − 1 ( R ( ζ ) ) = ( I + P ( ζ ) ) ( R ) = ( P ( ζ ) )i ( R )−1 . ∑ k R k 0 R k 0  i =1

Replacing Eq. (19) in Eq. (13), one obtains:

(19)

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

n n ∞  pC( qij) ( ξ k ) T j0 ψ i ( x ) ∑∑ Tn ( x, ξ k ) = i =1 j =1 q =1  ∞  q T = ( T0 ) ∑ ( PC ( ξ k ) ) Ψ ( x ) ;  q =1   m m ∞ u x, ξ , ζ = pR( qij) ( ζ k ) f j x, ξ k ϕi ( x ) k ∑∑ k  m i =1 j =1 q =1  ∞ T q  = F ξ ( ) ( ) ( PR ( ζ k ) ) Φ ( x ) ; ∑ k  q =1 

_____________________11

∑

205

) ∑

(

206 207

(

(20)

)

where

( P ( ξ )) C

k

q

q =  pC( ij) ( ξ k ) 

n× n

∧

( P ( ζ )) R

k

q

q =  pR( ij) ( ζ k ) 

m× m

.

(21)

208

In order to obtain the approximated solutions ( um , Tn ) , it is necessary to truncate the infinite series

209

appearing in Eq. (17). Hence, approximated solutions ( ums , Tnr ) are obtained as: n n r  q T x , ξ = pC( ij) ( ξ k ) T j0 ψ i ( x ) ∑∑  nr ( k ) i =1 j =1 q =1  r  q T = ( T0 ) ∑ ( PC ( ξ k ) ) Ψ ( x ) ;  q =1   m m s u x , ξ , ζ = pR( qij) ( ζ k ) f j x, ξ k ϕi ( x ) k ∑∑ k  ms i =1 j =1 q =1  s T q  = F ξ ( ) ( ) ( PR ( ζ k ) ) Φ ( x ) . ∑ k  q =1 

∑

210

(

) ∑

(

(22)

)

211

The requirement that PC ( ξk ) , PR ( ζ k ) < 1 is a well-known limitation of the Neumann series. The

212

authors are unaware of theoretical results warranting the above for general problems. It is known,

213

however, that larger variance in the random coefficients increases the above norms, eventually

214

leading to non-convergence of the series. Hence, there are practical limitations on the parameter

215

variances that can be handled using the Neumann series and the solutions presented in this article.

216

The following references can be consulted for limitations on the use of the Neumann series:

217

(Yamazaki 1989; Araújo & Awruch 1994; Chakraborty & Dey 1996; Chakraborty & Dey 1998; Lei

_____________________12

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

218

& Qiu 2000; Chakraborty & Sarkar 2000; Chakraborty & Bhattacharyya 2002; Li et al. 2006;

219

Schevenels et al. 2007).

220 221 222

5.

LOWER AND UPPER BOUNDS FOR REALIZATIONS OF THE STOCHASTIC DISPLACEMENT AND TEMPERATURE FIELDS

223

In this section the main contribution of this article is presented: based on properties of the Neumann

224

expansion and classical results of algebra, a methodology is proposed to derive lower and upper

225

bounds for samples of the displacement and temperature fields, ( um ,Tn ) . The main ideas are: to

226

avoid direct computation of the linear system of equations (Eq. 9) and to obtain an approximated

227

solution, ( ums , Tnr ) , using a truly small (n=0 or n=1) number of terms in the Neumann expansion.

228

If the set of thermo-mechanical load samples,

229

0 < PC ( ξi ) , PR ( ζ i ) < 1, ∀i ∈{1, ..., N } ; then the numerical approximations to the displacement and

230

temperature fields can be represented in terms of the Neumann series using matrices ( I + PR ( ζ k ) )

231

and ( I + PC ( ξk ) ) , Eq. (22).For the kth realization of the displacement and temperature fields, and for

232

a fixed x ∈ [ 0, l ] , the distance between response realizations obtained by direct Monte Carlo ( um , Tn )

233

, Eq. (7), and by the Neumann series ( ums , Tnr ) , Eq. (22), can be estimated as:

234

235

{κ ( x, ζ ( ω ) ) , EA ( x, ξ ( ω) )} i

i

i =1

, are such that

−1

−1

∞ q  T − T x , ξ ≤ PC ( ξ k ) U 0 Ψ ( x ) , ( ) ( ) ∑ k  m mr q = r +1   ∞ q  ( u − u ) ( x, ζ , ξ ) ≤ PR ( ζ k ) F ( ξ k ) Φ ( x ) . ∑ m ms k k  q = s +1

Based on properties of the Neumann series, Eq. (18), the terms

(23)

∞

∑ q = r +1

236

N

can be evaluated as:

PC ( ξ k )

q

and

∞

∑ q = s +1

PR ( ζ k )

q

_____________________13

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

r +1  ∞ PC ( ξ k ) q  ∑ PC ( ξ k ) = ; 1 − PC ( ξ k )  q =r +1  s +1 PR ( ζ k )  ∞ q .  ∑ PR ( ζ k ) = 1 − PR ( ζ k )  q =r +1

237

(24)

238

Substituting Eq. (24) in Eq. (23), the distance between the direct Monte Carlo andthe Monte Carlo –

239

Neumann solutions become:  (Tm − Tmr )( x, ξ k ) ≤ CT ( r , ξ k ) T0 Ψ ( x ) , ∀ ( x , ξ k ) ∈ [ 0, l ] × Γ;   ( u m − ums )( x, ζ k , ξ k ) ≤ CU ( s , ζ k )  U 0  +CT ( r , ξ k ) T0  Φ ( x ) , ∀ ( x, ζ k , ξ k ) ∈ [ 0, l ] × Γ × Π ; 

240

241

242

(25)

where r +1  PC ( ξ k ) C T ( r , ξ k ) = ; 1 − PC ( ξ k )   s +1  PR ( ζ k ) . C U ( s , ζ k ) = 1 − PR ( ζ k ) 

(26)

243

Coefficients C T (⋅, ⋅ ) and CU ( ⋅, ⋅ ) , for the kth thermo-mechanical sample, depend on the number

244

of terms ("r" and "s") used in the Neumann series for matrices ( I + PC ( ξk ) ) and ( I + PR ( ζ k ) ) ,and

245

on the numerical values of PC ( ξ k ) and PR ( ζ k ) . FromEqs.(25) and(26) one obtains the lower and

246

upper bounds, for the kthsample, of the displacement and temperature fields:

247

ϖ  ( x, ξ k ) ≤ Tn ( x, ξ k ) ≤ θ ( x , ξ k ) , ∀ ( x , ξ k ) ∈ [ 0, l ] × Γ;   α ( x, ζ k , ξ k ) ≤ um ( x , ζ k , ξ k ) ≤ β ( x , ζ k , ξ k ) , ∀ ( x, ζ k , ξ k ) ∈ [ 0, l ] × Γ × Π;

−1

−1

(27)

248

where ϖ ( ⋅, ⋅) and θ ( ⋅, ⋅) are lower and upper bounds, respectively, for the temperature field; and α (⋅, ⋅)

249

and β ( ⋅, ⋅) are lower and upper bounds, respectively, for the displacement field. The quotas

250

( ϖ, θ, α, β ) are defined as:

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

ϖ ( x, ξ k ) = Tmr ( x, ξ k ) + χ ( r , ξ k ) T0 Φ ( x ) ,  θ ( x, ξ k ) = Tmr ( x, ξ k ) + λ ( r , ξ k ) T0 Φ ( x ) , ∀ ( x, ξ k ) ∈ [ 0, l ] × Γ  α ( x, ζ k , ξ k ) = ums ( x, ζ k , ξ k ) + τ ( s, ζ k )  U 0 +λ ( r , ξ k ) T0  Ψ ( x ) ,  β ( x, ζ k , ξ k ) = ums ( x, ζ k , ξ k ) + υ ( s , ζ k )  U 0  +λ ( r , ξ k ) T0  Ψ ( x ) , ∀ ( r , s , ζ k , ξ k ) ∈ [ 0, l ] × Γ × Π ; 

251

252

_____________________14

(28)

where  χ ( r , ξ k ) = −λ ( r , ξ k ) ;   τ ( s , ζ k ) = −υ ( s , ζ k ) ;   λ ( r , ξ k ) = CT ( r , ξ k ) ; υ s, ζ = C s, ζ . U ( k )  ( k)

253

(29a) (29b) (29c) (29d)

254 255

Unfortunately, there is a “catch” in this formulation. In order to compute coefficients " λ" and

256

" υ" in Eq. (29), one needs to solve for the norm P , which may involve high computational cost

257

and hence turn the procedure useless. In order to circumvent this difficulty, we resort to the

258

equivalence between matrix norms (Golub & Van Loan, 2012): cη P

259

η

≤ P ≤ Cη P η ,

(30)

260

where 0 < cη ≤ Cη < +∞ ,are equivalence constants. Based on this equivalence (Eq. 30),and on Eqs.

261

(26), (29c) and (29d), five alternatives are proposed to evaluate P and the coefficients " λ" and

262

" υ" :

_____________________15

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

 n +1  m. PC ( ξ k ) 1 λ1 ( n, ξ k ) = ; 1 − m . PC ( ξ k ) 1   n +1  PC ( ξ k ) 2 λ 2 ( n, ξ k ) = ; 1 − PC ( ξ k ) 2   n +1  m . P ξ ( ) C k  ∞ ;  λ 3 ( n, ξ k ) = 1 − m. PC ( ξ k ) ∞   n +1 PC ( ξ k ) F  ;  λ 4 ( n, ξ k ) = 1 − PC ( ξ k ) F   n +1  m. max pCij ( ξ k ) i,j  λ ( n, ξ ) = ; k  5 1 − m. max p ξ ( ) Cij k  i,j

 n +1  m . PR ( ζ k ) 1  υ1 ( n, ζ k ) = ; 1 − m. PR ( ζ k ) 1   n +1  PR ( ζ k ) 2  υ 2 ( n, ζ k ) = ; 1 − PR ( ζ k ) 2   n +1  m . P ζ ( ) R k  ∞ ;  υ3 ( n, ζ k ) = 1 − m . PR ( ζ k ) ∞   n +1 PR ( ζ k ) F  ;  υ 4 ( n, ζ k ) = 1 − PR ( ζ k ) F   n +1  m. max pRij ( ξ k ) i, j  υ ( n, ζ ) = ; k  5 1 − m. max p ξ ( ) Rij k  i,j

 m   PC ( ξ k ) 1 = max  ∑ pCij ( ξ k )  ; 1≤ j ≤ m  i =1    m   pCij ( ξ k )  ;  PC ( ξ k ) ∞ = max  ∑ 1≤ i ≤ m  j =1    m m 2  P (ξ ) = pCij ( ξ k ) ; ∑ ∑ C k F  i =1 j =1

 m   PR ( ζ k ) 1 = max  ∑ pRij ( ζ k )  ; 1≤ j ≤ m  i =1    m   pRij ( ζ k )  ;  PR ( ζ k ) ∞ = max  ∑ 1≤ i ≤ m  j =1    m m 2  P (ζ ) = pRij ( ζ k ) . ∑ ∑ R k F  i =1 j =1

(

)

(

263

)

(

{

(

264

265

{

}) })

(

)

(

)

(

{

(

{

}) })

( 31a ) ( 31b ) ( 31c ) ( 31d ) ( 31e )

where

(32)

266 267

Within this paper, for simplicity, the norms in Eq. (31) are referred to as follows: ⋅ 1 is called

268

norm one, ⋅

269

maximum norm.

2

is norm two, ⋅ ∞ is the infinity norm, ⋅

F

is the Frobenius norm and max { ⋅ } is the

270

The definition of coefficients " λ" and " υ" are based on estimates about the sum in Neumann

271

series; hence the equivalence norm, ⋅ η , in Eq. (31), must comply with the hypothesis in Eq. (17),

272

that is: Cη P

273 274

η

< 1.

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

275 276

6.

_____________________16

LOWER AND UPPER BOUNDS FOR THE EXPECTED VALUE AND SECOND ORDER MOMENT OF DISPLACEMENT AND TEMPERATURE FIELDS

277

From the ensemble of computed lower and upper bounds, one can compute lower and upper bounds

278

for the expected value and for the second order moment of the displacement and temperature fields.

279

Let

({u ( x ζ (ω ) ξ (ω ))} m

,

j

,

j

N j =1

{ ( ( ))}

, Tn x, ξ ω j

N j =1

, ∀x ∈ [0, l ]

) be the ensemble of realizations of

280

numerical approximations of the displacement and temperature responses. Then, estimates for

281

expected value and second order moment are defined as: N  1 µ x = um x , ζ ( ω j ) , ξ ( ω j ) , ( ) ( ) u ∑ N  m i =1  N  1 µ x = ( ) ( )  Tn N ∑ Tn x, ξ ( ω j ) , ∀x ∈ [ 0, l ] ; i =1   N  ( 2) 1 µ um ( x, y ) = ( N −1 ) ∑ um x , ζ ( ω j ) , ξ ( ω j ) um y , ζ ( ω j ) , ξ ( ω j ) ,  j =1  N 2  ( 2) µ 1 x , y = ) ( N −1 ) ∑ Tn x, ξ ( ω j ) Tn y, ξ ( ω j ) , ∀ ( x, y ) ∈ [0, l ] .  Tn ( j =1 

(

)

(

282

)

(

) (

(

) (

(33)

)

)

283

Inserting Eqs. (27)-(29) in Eq. (33), one obtains lower and upper bounds for the expected value and

284

second order moment of the displacement and temperature fields:

285

   µ ϖi ( x ) ≤ µTn ( x ) ≤ µθi ( x ) ,    µ αi ( x ) ≤ µum ( x ) ≤ µβi ( x ) , ∀x ∈ [ 0, l ] , i = 1, ...,5 ;   ( 2)  ( 2)  ( 2) µ ϖi ( x, y ) ≤ µTn ( x, y ) ≤ µ θi ( x, y ) ,   ( 2) 2  ( 2)  ( 2) µ αi ( x, y ) ≤ µum ( x, y ) ≤ µβi ( x, y ) , ∀ ( x, y ) ∈ [ 0, l ] , i = 1, ...,5.

(34)

286

In the following, direct Monte Carlo simulation is employed to evaluate the performance of the

287

proposed lower and upper bounds, in the solution of two linear thermo-elastic problems. In order to

288

evaluate the accuracy of the proposed lower and upper bounds, deviation (error) functions need to be

289

established. The relative deviation, with respect to direct Monte Carlo simulation, of the “k” order

290

moment lower and upper bounds, εµ ( ) k

( ϖi ,αi )

, εµ ( k )

(θi ,βi )

:  0, l  →


, are given by:

_____________________17

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

k   µ(( α)i ,ϖi )  100 % . 1 − ( )   µ ((uk )m ,Tn ) ( x )  , ∀ x ∈ 0, l ; εµ ( k ) ( x ) =  ( αi ,ϖi ) 0, ∀ x ∈ 0, l , i = 1, ...,5 ;     µ((βk )i ,θi )  (100% ) . 1 − µ ( k ) ( x ) , ∀ x ∈ 0, l ; (um ,Tn ) εµ ( k ) ( x ) =    (βi ,θi )  0, ∀ x ∈ 0, l , i = 1, ...,5 . 

( ) { }

291

(35)

( )

{ }

292 293

7.

NUMERICAL RESULTS

294

In this section, numerical results are presented for two examples involving a bar subject to a

295

deterministic thermo-mechanical loading, with uncertainty in stiffness and thermal conductivity. The

296

problem is illustrated in Figure 1. The deterministic thermo-mechanical loading, q ( x ) , f ( x ) , is

297

given by:

(

)

 f ( x ) = 8400.( 2 x − 1) x kPa.m, ∀x ∈ ( 0,1) ;   q ( x ) = 50 kW.m, ∀x ∈ ( 0,1) .

298

299

For

both

300

(β = 23,5.10

examples, −6

the linear

dilatation

(36) coefficient

is

deterministic,

and

given

by

K −1 ) . Thebar is clamped at both ends and the length is one meter.

301

Modelling of the uncertainty in stiffness and conductivity coefficients is made by

302

parameterized stochastic processes. The expected value and coefficient of variation for these

303

processes are given by: 1 δ EA ( x ) = 10 , ∀x ∈ [ 0,1];  1 δ κ ( x ) = 10 , ∀x ∈ [ 0,1] .

304

µ EA ( x ) = 7.106 N.m 2 , ∀x ∈ [ 0,1] ;  µ κ ( x ) = 250 kW mK , ∀x ∈ [ 0,1] ;

305

In

306

{ξ ( ω )}

307

employed. Properties of these vectors are given in Eq. (6). All random variables ( ζ k , ξk ) have

j

N j =1

order

to

characterize

= {ξ1 ( ω j ) , ξ2 ( ω j ) , ξ3 ( ω j ) , ξ 4 ( ω j )}

N j =1

and

the

uncertainty,

{ζ ( ω )} j

N j =1

(37) random

= {ζ1 ( ω j ) , ζ 2 ( ω j ) , ζ 3 ( ω j ) , ζ 4 ( ω j )}

vectors N j =1

are

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

_____________________18

308

uniform distribution. All computer simulations are performed in MATLAB. In order to reduce

309

spurious correlations, Latin Hypercube Sampling is employed (Olsson &Sandberg, 2002).

310

The reference solution is direct Monte Carlo simulation. In this solution, the kth realization of

311

the displacement and temperature fields is obtained from the kth thermo-mechanical sample, which

312

consists in inserting in Eq. (5) the sampled values of vectors ( ζ k = ζ ( ωk ) , ξ k = ξ ( ωk )) ; evaluate from

313

Eq. (10) the stiffness and conductivity matrices; solve the linear system of equations (Eq. 9) for

314

(U (ζ

k

, ξ k ) , T ( ξ k ) ) and compute the coefficients of the linear combinations in Eq. (7).

315

A total of thirty thousand samples are considered (N=30,000). Figure 2 shows convergence

316

  plots of expected value and second order moment, µum , σu2m , with respect to the number of samples

317

(N),for the random variable displacement u m

(

)

( 2l , ζ, ξ ) , obtained by fixing x = 2l .

318

One observes in Figure 2 that for N ≥ 15, 000 realizations, the expected value and second

319

 order moment of mid-spam displacements stabilize at values µ um ( 12 ) ≅ −0.00863 m and

320

 σ u2m ( 12 ) ≅ 2.03524 × 10−7 m 2 . Figure 2 shows that using N=30,000 samples is sufficient to obtain

321

reliable results via direct Monte Carlo simulation. Note that a large number of samples is used only

322

to exclude effects of statistical sampling errors. For larger, realistic problems, a smaller number of

323

samples would be employed, following the usual compromise between accuracy and computational

324

effort. For any number of samples, the proposed bounds will be more efficient to compute than a

325

brute force MC simulation.

326 327

7.1 Example 1: Random stiffness (EA)

328

In this example, the uncertainty is assumed on the“unitary length” axial stiffness coefficient,

329

EA : [ 0,1] × Ω →  a, a  , which is modeled as a parameterized stochastic process:

330

EA ( x, ζ ( ω) ) = µ EA + (

3 2

2

) .σ ∑ ζ ( ω) cos ( k π4x ) + ζ ( ω) sin ( k π4x ) , EA

k =1

2.k −1

2.k

(38)

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

_____________________19

331

where σ EA is the standard deviation of the stiffness coefficient.In this example, the conductivity

332

coefficient is deterministic and equal to the mean (Eq. 37). Hence, the displacement field is a

333

function of the random outcomes of the stiffness coefficient only.

334

Figure 3 shows relative deviations of lower and upper bounds for the expected value and

335

second order moment of the displacement field, Eq. (22). The deviations are computed by comparing

336

the proposed Monte Carlo – Neumann bounds with the direct Monte Carlo solution, following Eq.

337

(35). One observes in Fig. 3 that the most accurate bounds are obtained when coefficient " λ " is

338

given by Eq. (31b).The largest deviations are obtained for coefficient " λ " evaluated by the

339

maximum norm, ( máx ⋅ ) , given by Eq. (31e). One can also observe that the relative deviations are

340

larger for second order moment then for expected values.

341

Table 1 shows numerical values of the relative deviations in expected value and second order

342

moment, computed via Monte Carlo-Neumann using different norms, for n=0 and n=1. The table

343

shows maximum deviation values within the domain, max { ⋅ } , and deviations computed for

344

x = 107 m . One observes that the deviations in second order moment are larger than deviations in

345

expected values. Also, one observes that using n=1 in the Monte Carlo-Neumann solution

346

significantly reduces the relative deviations. The largest deviations are obtained when coefficient

347

" λ " is computed using the maximum norm, Eq. (31e). When the bounds are evaluated by norm 2,

348

Eq. (31b), the smallest deviations are obtained.

x∈[ 0,1]

349

Using norm 2 leads to the most accurate results, but there is a computational penalty. Table 2

350

compares the computation time taken to obtain the direct Monte Carlo solutions ( TMC ), and the

351

solutions via Monte Carlo – Neumann with different norms (T1 ,…, T5 ) , Eq. (31), and n values. One

352

observes that there is a gain in computational time for all Monte Carlo – Neumann solutions, when

353

compared to direct Monte Carlo ( TMC ). However, this gain is minimum when norm 2 is used to

_____________________20

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

( ⋅ ).

354

compute coefficient " λ " . The largest gain in computational time was obtained using norm 1,

355

This norm, however, leads to larger deviations, as observed in Table 1. Therefore, there is a

356

compromise between the accuracy of the results and the computational time to achieve it, using

357

different norms. A good compromise between efficiency and accuracy is obtained, for instance,

358

using the Frobenius norm (norm 4, Eq. 31d). When the bounds are evaluated via Neumann series

359

with n=0, the computation time is zero, (TSN = 0 ) , since ( I + PR ( ⋅) ) ≈ I . Also, it is noted that

360

computing the bounds for n=2 or larger leads to computation times larger than brute Monte Carlo.

1

−1

361 362

7.2 Example 2: Random conductivity coefficient ( κ )

363

In this example, uncertainty over the conductivity coefficient is considered, and modelled as a

364

parameterized stochastic process κ : [ 0,1] × Ω →  b, b  , given by:

365

κ ( x, ξ ( ω) ) = µ κ +

3

( ) .σκ ∑ ξ ( ω) cos ( kπ4x ) + ξ ( ω) sin ( k π4x ) , 3 3

k =1

2.k −1

2.k

(39)

366

where σ κ is the standard deviation. In this example, the axial stiffness is deterministic and equal to

367

the mean (Eq. 37). However, the displacement field is a function of the random outcomes of the

368

conductivity.

369

Figure 4 shows relative deviations of lower and upper bounds for the expected value and

370

second order moment of the displacement field, Eq. (22). As for Example 1, one observes that the

371

relative deviations are larger for second order moment then for expected values. However, for

372

Example 2 the deviations are significantly smaller. Again, the most accurate bounds are obtained

373

when coefficient " λ " is evaluated from Eq. (31b). The largest deviations are obtained for coefficient

374

" λ " evaluated by the maximum norm, norm 5, given by Eq. (31e).

375

Table 3 shows numerical values of the relative deviations in expected value and second order

376

moment, computed via Monte Carlo-Neumann using different norms for n=0 and n=1. The table

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

_____________________21

377

shows maximum deviation values within the domain, max { ⋅ } , and computed for x = 107 m . One

378

observes that the deviations in second order moment are larger than deviations in expected values,

379

but also that deviations are smaller for example 2 than for example 1. This indicates that when

380

uncertainty is on the conductivity coefficient “ κ ”, propagation of the uncertainty through the

381

thermo-elastic system is smaller.

x∈[ 0,1]

382

Computation times for the direct Monte Carlo solution and for the Monte Carlo-Neumann

383

solutions using different norms are shown in Table 2. Computation times are smaller for all Monte

384

Carlo-Neumann solutions, using the bounds proposed herein, than for direct Monte Carlo. The most

385

accurate bounds are obtained when coefficient " λ " is computed using norm 2 (Table 3), Eq. (31b),

386

but this is also the most costly norm to compute (Table 2). The largest deviations are obtained for

387

norm 5, Eq. (31e). The fastest norm to compute is norm 1, but deviations are large. Also for this

388

example, a good compromise between efficiency and accuracy is obtained by using the Frobenius

389

norm (norm 4, Eq. 31d).

390 391

7.3 Example 3: Generic Linear Thermo-elastic Problem

392

In this example, solution to a generic, linear, one-dimensional thermo-elastic problem is addressed,

393

using the Neumann bounds proposed in this paper. The discretized solution to a general one-

394

dimensional linear problem in thermo-elasticity, with uncertainty in the mechanical or thermo-

395

physical parameters, can be cast as (reproduced from Eq. 9): n

396

m

 For fixed {ζ k , ξ k } ,find ( U ( ζ k , ξ k ) , T ( ξ k ) ) ∈
×
such that :   R ( ζ k ) U ( ζ k , ξ k ) + DT ( ξ k ) = F;  ( CT )( ξ k ) = Q.

(40)

397

Clearly, similar formulations can be obtained for two and three-dimensional thermo-elastic

398

problems.

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

_____________________22

399

In this example, lower and upper bounds are obtained for the kth realization of the displacement

400

and temperature fields, when the thermal conductivity parameter is assumed random (similar to

401

example 2). The distance between displacement vectors U n and U , computed using the Neumann

402

series with n terms ( n = 0 or n = 1) and Monte Carlo simulation, respectively, is computed:

403

U n − U ≤ λ R −1D T0

U = U + λ ( R )−1 D T 1; n 0  ⇒  −1 U = U n + χ ( R ) D T0 1; 

(41)

T

404

where 1 = [1,…,1] ∈
m , λ and χ are the coefficients defined in Eqs. (31d) and (29a); and

405

T U= [ u1 ,…, u m ] and U= u1,…, u m  are the lower and upper bound vectors of U= [u1 ,…, u m ] . The

406

expected values and second order moments are computed for components of vectors U , U and U .

407

Relative deviations in expected values and second order moments are computed, similar to Eq. (35):

T

T



408

µ(u ,ui )  i εµ (ui ,ui ) = (100% ) . 1 − µui , i = 1, ..., m;    µ(2u) ,u ( i i) ε  , i = 1, ..., m. = 100 % . 1 − ( )  µu(2i)  µ((2u)i ,ui )

409

Results are computed for the same numerical values of the previous examples. The matrices

410

and vectors in Eq. (40) are those computed for example 2 above. In this example, only the

411

Frobenious norm (Eq. 31d) is used. Results are computed for 30000 samples. Results are presented

412

in Tables 4 and 5, when vector U n is computed using Eq. (41) and for n=0 and n=1 terms in the

413

Neumann Series, respectively, for a system with dimension m=3.

(42)

414

In Tables 4 and 5 one observes that smaller deviations are obtained for the components of

415

vectors U and U when the Neumann series solution U n is computed with n=1, in comparison to

416

n=0. For both approximations (n=0 and n=1 terms), the relative deviation for expected value is

417

smaller than for second order moment.

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

Finally, it is observed that accurate lower and upper bounds are obtained using the proposed

418 419

_____________________23

methodology.

420 421

8.

CONCLUDING REMARKS

422

In this paper, a well-known property of the Neumann series was explored to derive lower and upper

423

bounds for expected value and second order moment of the stochastic temperature and displacement

424

responses, for thermo-elasticity problems. To improve efficiency of the proposed scheme,

425

equivalence between norms of finite dimension operators were considered. The uncertainties in axial

426

stiffness and conductivity coefficients were represented as parameterized stochastic processes. The

427

developed scheme was illustrated by means of application to two specific linear one-dimensional

428

thermo-elastic problems and to a general linear thermo-elastic problem. Solutions employing the

429

proposed Monte Carlo – Neumann scheme were computed for five equivalent norms using one or

430

two terms in the Neumann expansion (n=0 or n=1). Deviations of the proposed Monte Carlo –

431

Neumann bounds were computed by comparing to a reference solution, computed by direct Monte

432

Carlo simulation. It was shown that accurate bounds can be obtained using the proposed Monte Carlo

433

– Neumann scheme, using as few as one or two terms in the Neumann expansion (n=0 or n=1).

434

Using more terms in the Neumann solution led to computational times larger than brute force Monte

435

Carlo. The most accurate bounds were obtained with norm 2

436

computationally expensive bounds to compute, leading to little gain in comparison to direct Monte

437

Carlo simulation. A good compromise between accuracy and efficiency was obtained by using the

438

Frobenius norm

439

holds for other problems as well. The proposed methodology can be readily extended to non-linear

440

problems, in conjunction with different linearization techniques, such as Newton-Rhapson. Recall

441

that the proposed methodology completely circumvents solution of the linear system. The proposed

( ⋅ ) , but these were also the most 2

( ⋅ ) , for both problems studied herein. It remains to be investigated weather this F

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

_____________________24

442

Monte Carlo – Neumann bounding scheme proposed herein is a viable alternative for the solution of

443

uncertainty propagation problems in mechanics.

444 445

ACKNOWLEDGEMENTS

446

The authors gratefully acknowledge the Brazilian National Research Council (CNPq) for sponsoring

447

this research.

448 449

REFERENCES

450 451

AdhikariS. (2011): A reduced spectral function approach for the stochastic finite element analysisComputer Methods in Applied Mechanics and Engineering 200, 1804-1821.

452 453

Araújo J.M., AwruchA.M., (1994): On stochastic finite elements for structural analysisComputers & Structures 52, 461-469.

454 455

Ávila S Jr. CR, Beck AT, Mantovani G, Azikri H. (2010) “Galerkin Solution of Stochastic Beam Bending on Winkler Foundations.” Computer Modeling in Engineering & Science; v. 1729, p. 1-31.

456 457 458

Babuska I; Tempone R, Zouraris GE. (2005) “Solving elliptic boundary value problems with uncertain coefficients by the finite element method: the stochastic formulation.” Computer Methods in Applied Mechanics and Engineering, v. 194, n. 12-16, p. 1251-1294.

459 460

Beck AT, Silva Jr. CRA, 2011: Timoshenko versus Euler beam theory: Pitfalls of a deterministic approach. Structural Safety 33, 19-25.

461 462

Bhattacharyya B., ChakrabortyS., NE-MCS technique for stochastic structural response sensitivity, Computer Methods in Applied Mechanics and Engineering, v. 191, n. 49–50, 2002, p. 5631-5645.

463

Brenner SC, Scott LR. (1994), The Mathematical Theory of Finite Element Methods. Springer-Verlag.

464 465

Chakraborty S., DeyS.S. (1995): Stochastic finite element method for spatial distribution of material properties and external loading, Computers & Structures 55, 41-45.

466 467

Chakraborty S., DeyS.S. (1996): Stochastic finite element simulation of random structure on uncertain foundation under random loading,International Journal of Mechanical Sciences 38, 1209-1218.

468 469

Chakraborty S., DeyS.S. (1998) A stochastic finite element dynamic analysis of structures with uncertain parametersInternational Journal of Mechanical Sciences, v. 40, n. 11,p. 1071-1087.

470 471

Chakraborty S., SarkarS. K., (2000): Analysis of a curved beam on uncertain elastic foundationFinite Elements in Analysis and Design 36, 73-82.

472 473

Chakraborty S., BhattacharyyaB. (2002) An efficient 3D stochastic finite element method, International Journal of Solids and Structures, v. 39, n. 9, p. 2465-2475.

474

Ghanem R. &Spanos P.D., 1991: “Stochastic Finite Elements: A Spectral Approach”, Dover, NY.

475

Golub GH, Van LoanCF, (2012) Matrix Computations, Johns Hopkins Studies in the Mathematical Sciences.

476 477

Hyuk-Chun Noh, Taehyo Park(2006) Monte Carlo simulation-compatible stochastic field for application to expansion-based stochastic finite element methodComputers & Structures, v. 84, n. 31–32,p. 2363-2372.

478

Kincaid D Cheney W (2002). Numerical Analysis: mathematics of scientific computing, 3a ed., AMS.

_____________________25

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

479 480

Lei Z., QiuC. (2000)Neumann dynamic stochastic finite element method of vibration for structures with stochastic parameters to random excitationComputers & Structures, v. 77, n. 6, p. 651-657.

481 482

Li C.F., Feng Y.T., OwenD.R.J. (2006) Explicit solution to the stochastic system of linear algebraic equations, Computer Methods in Applied Mechanics and Engineering, v. 195, n. 44–47, p. 6560-6576.

483 484

Rahman S., XuH. A univariate dimension-reduction method for multi-dimensional integration in stochastic mechanicsProbabilistic Engineering Mechanics, v. 19, n. 4, October 2004, p. 393-408

485

Rao MM, Swift JR. (2010), Probability Theory with Applications. Springer; 2nd ed.

486 487

Schevenels M., Lombaert G., Degrande G., ClouteauD.The wave propagation in a beam on a random elastic foundation, Probabilistic Engineering Mechanics, v. 22, n. 2, 2007, p. 150-158.

488 489

Silva Jr. CRA, Beck AT, Rosa E, 2009: Solution of the Stochastic Beam Bending Problem by Galerkin Method and the Askey-Wiener Scheme. Latin American Journal of Solids and Structures 6, 51 - 72.

490 491

Silva Jr. CRA, Beck AT, 2010: Bending of stochastic Kirchhoff plates on Winkler foundations via the Galerkin method and the Askey-Wiener scheme. Probabilistic Engineering Mechanics 25, 172 - 182.

492 493

Silva Jr. CRA, Beck AT, 2011: Chaos-Galerkin Solution of Stochastic Timoshenko Bending Problems. Computers & Structures 89, 599 - 611.

494 495

Olsson AMJ,SandbergGE.(2002) “Latin Hypercube Sampling for Analysis.”Journal of Engineering Mechanics; v. 128, n. 1, p. 121-125.

496

Yoshida K. (1978) Functional analysisSpringer Berlin.

Stochastic

Finite

Element

497 498 499

FIGURE CAPTIONS

500 501 502

Figure 1: Bar subject to thermo-mechanical loading.

503

Figure 2: Convergence of Monte Carlo simulation results for the mean and second order moment of

504

mid-spam displacement ( x =

1 2

).

505 506 507

Figure 3: Relative error functions for expected value εµ , εµ , i = 2, 4,5 (above) and for second αi

βi

order moment ε  ( 2) , ε  ( 2) , i = 2,4,5 (below). µ µ α i

β i

508 509 510

Figure 4: Relative error functions for expected value εµ , εµ , i = 2, 4,5 (above) and for second αi

order moment ε  ( 2) , ε  ( 2) , i = 2,4,5 (below). µ µ α i

511

β i

βi

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

_____________________26

512

Table 1: Relative deviations in expected value and second order moment, computed via

513

Monte Carlo-Neumann using different norms, for n=0 and n=1 (example 1). Expected Value Norm

n

{

max εµ x∈[0,1]

αi

( x)

}%

εµ α ( 107 ) %

  max  εµ ( x )  % x∈[ 0,1] β  i 

εµ ( 107 ) % β

i

i

n=0

21.376641981258

21.373505564049

20.8641793659111

20.432999912677

n=1

4.8426006707301

3.9089626093706

3.90896260937064

3.9089626093706

n=0

9.1602776186399

9.1566538674697

8.60358798967606

8.2161482160989

n=1

1.3272218427113

0.3935837813525

0.393583781352569

0.3935837813525

n=0

20.147182045325

20.143996582624

19.6302684051304

19.203490931252

n=1

4.4101538187991

3.4765157574408

3.47651575744081

3.4765157574408

n=0

13.011962790467

13.008492690329

12.4692174519237

12.067987038957

n=1

2.2876171624128

1.3539791010550

1.35397910105506

1.3539791010550

n=0

25.397922645337

25.394946644645

24.9000183084938

24.454440993273

n=1

6.8925102542441

5.9588721928854

5.95887219288547

5.9588721928854

1

2

3

4

5

Second order moment Norm

n

  máx  ε (2) ( x )  % x∈[0,1] µ  αi 

ε  (2) ( 107 ) % µα

i

    máx  ε  ( 2 ) ( x )  % x∈[0,1] µ  βi 

ε  (2) ( 107 ) % µβ

i

n=0

37.081563811368

37.069940729798

46.4299932763907

45.756212384637

n=1

9.2513852286976

9.2449949210784

8.80652854150028

8.0888113292324

n=0

17.664604039865

17.649394010035

17.3090127747284

16.769229335595

n=1

2.6400606545767

2.6324881720284

1.47959804281452

0.7760429863682

n=0

35.262203886142

35.250244709251

43.3638269045107

42.704154691749

n=1

8.4507182779713

8.444194588629

7.89082593054284

7.1746868377654

n=0

24.056554641698

24.042525414504

26.2192420780234

25.638459062547

n=1

4.4776555447964

4.4704254596371

3.44971098949154

2.7427176030510

n=0

41.254493957511

41.243641752284

58.2450439965844

57.516897213540

n=1

12.529144976908

12.523374151606

13.6581286827352

12.930982466756

1

2

3 4

5

514

_____________________27

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

515

Table 2: Computation (clock) time to obtain the direct Monte Carlo ( TMC ) and Monte Carlo –

516

Neumann solutions, using different norms ( Ti ) and n values. Norm Example

1

2

517 518

TMC

n

TSN

T1

T2

T3

T4

T5

n=0

0

138.14

692.44

157.03

249.97

291.63

n=1

172.88

180.65

624.29

194.49

228.81

360.17

n=0

0

144.31

728.09

213.80

227.87

282.97

n=1

138.29

182.86

649.34

205.72

188.08

229.37

874.90

962.12

_____________________28

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

519

Table 3: Relative deviations in expected value and second order moment, computed via

520

Monte Carlo-Neumann using different norms, for n=0 and n=1 (example 2). Expected Value Norm

n

{

máx εµ x∈[ 0,1]

αi

( x)

}

/ 10−6%

εµ α ( 107 ) / 10

−6

%

i

{

máx εµ x∈[0,1]

βi

( x)

}/

−6

10 %

εµ ( 107 ) / 10 β

−7

%

i

n=0

5.965287275

3.06264025255

5.996997143

3.04079494917

n=1

1.2408067129

0.65210992163

1.272172078

0.63007532524

n=0

2.470117976

1.2793176518

2.5018263238

1.25747812163

n=1

0.231092134

0.1369211410

0.2624580774

0.11489609264

n=0

5.614841958

2.88382173696

5.6465513820

2.86199264287

n=1

1.117715320

0.58930338386

1.1490821849

0.56727622599

n=0

3.576681006

1.84390924751

3.6083905197

1.82207937626

n=1

0.506273856

0.27733393359

0.5376409539

0.25530133562

n=0

7.122282219

3.6529597124

7.1539917545

3.63112961920

n=1

1.824845496

0.9500960018

1.8562125170

0.92806184959

máx εµ ( ) ( x ) .10 −6 %

{

εµ ( 2) ( 107 ) .10−7%

1

2

3

4

5

Second order moment

Norm

}

n

  máx  εµ ( 2) ( x )  .10 −6% x∈[ 0,1]  αi 

n=0

11.930573495799

6.1252893868157

11.9939945975034

6.08158168269

n=1

2.4816110055780

1.3042086299819

2.54434617819044

1.26015109458

n=0

4.9402325341674

2.5586547325318

5.00365273659043

2.51494691738

n=1

0.4621820814953

0.2738492765885

0.52491733182336

0.22979196323

n=0

11.229682539415

5.7676776687998

11.2931035856079

5.72397085285

n=1

2.2354297524174

1.1786015496895

2.29816508046099

1.13454357020

n=0

7.1533608370089

3.6878445852650

7.21678135029436

3.64413610398

n=1

1.0125464022792

0.5546582082516

1.07528157489156

0.51059978467

n=0

14.244562474008

7.3059411853648

14.3079841752325

7.26223281510

n=1

3.6496900501781

1.9001895612547

3.71242554475515

1.85613280301

εµ ( ) ( 107 ) .10 2

−7

%

αi

x∈[ 0,1]

2

βi

βi

1

2

3

4

5

521 522

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

523 524

_____________________29

Table 4: Expected value, second order moment and relative deviations for the three components of vectors U , U and U , for Neumann solution with n=0.  µ ×10 −2 [ m ]

2 µ( ) ×10−5  m2 

εµ ×10−5 [ %]

ε  ( 2) × 10−5 [ %]

u1 u1

0.001466262072327

0.214999613128218

0.020610053348483

0.041220113264531

0.001466262069305

0.214999612241987

-

-

u1

0.001466262066284

0.214999611356002

0.020604473589063

0.041208701937292

u2 u2

-0.000967546030022

0.093617652609471

0.031485462149996

0.062971061402464

-0.000967546033068

0.093617653198992

-

-

u2

-0.000967546036065

0.093617653778908

0.030972714389449

0.061945219481465

u3 u3

0.000260860942038

0.006805069943771

0.115827050762482

0.231653894656217

0.000260860939017

0.006805069786129

-

-

u3

0.000260860935995

0.006805069628478

0.115833291365495

0.231666709321112

Component

µ

525 526 527

Table 5: Expected value, second order moment and relative deviations for the three components of vectors U , U and U , for Neumann solution with n=1.

 µ ×10−2 [ m]

2 µ( ) ×10−5  m2 

εµ ×10−5 [ %]

ε ( 2) ×10−5 [%]

u1 u1

0.001466262069513

0.214999612302978

0.141474777914224

0.28294792664581

0.001466262069305

0.214999612242144

-

-

u1

0.001466262069098

0.214999612181242

0.141628136263513

0.28326389020342

u2 u2

-0.000967546032836

0.093617653154025

0.240002365948933

0.47999175184744

-0.000967546033068

0.093617653198960

-

-

u2

-0.000967546033251

0.093617653234355

0.189035970506294

0.37807087591520

u3 u3

0.000260860939224

0.006805069796954

0.795742137209016

1.59147557974853

0.000260860939017

0.006805069786124

-

-

u3

0.000260860938809

0.006805069775296

0.795559470207738

1.59112103225871

Component

528 529 530

µ

figure 1

figure 2 a

−3

1.4

x 10

m

µu (1/2)

1.35 1.3 1.25 1.2 1.15 1.1 0

0.5

1

1.5

N

2

2.5

3 4

x 10

figure 2 b

−8

2

x 10

m

σ2u (1/2)

1.5 1 0.5 0 0

0.5

1

1.5

N

2

2.5

3 4

x 10

figure 3 a

εµ

10

εµ

β

β

β

4

εµ

εµ

α

5

α

2

εµ

α

4

5

i

5

0

α

i

β

εµ (x), εµ (x) [%]

2

εµ

−5

−10 0

0.1

0.2

0.3

0.4

0.5

x

0.6

0.7

0.8

0.9

1

figure 3 b

15

εµ(2)

εµ(2)

εµ(2)

2

4

5

α

α

α

εµ(2)

εµ(2)

β

β

2

4

εµ(2) β

5

i

5

α

i

β

εµ(2)(x), εµ(2)(x) [%]

10

0 −5 −10 −15 0

0.1

0.2

0.3

0.4

0.5

x

0.6

0.7

0.8

0.9

1

figure 4 a

−6

2

x 10

εα

εα

4

5

εβ

2

εβ

4

εβ

5

1 0

i

i

εα (x), εβ (x), [%]

2

εα

−1 −2 0

0.2

0.4

0.6 x

0.8

1

figure 4 b

−6

4

x 10

ε(2) α

ε(2) α

4

5

ε(2) β 2

ε(2) β 4

ε(2) β 5

2 0

i

i

(2) ε(2) ε (x), [%] α β

2

ε(2) α

−2 −4 0

0.2

0.4

0.6 x

0.8

1

Avila & Beck: Efficient Bounds for the MC-Neumann Sol. of Stochastic Thermo-Elasticity Problems

_____________________30

531

•

Solution of random thermo-elastic problems addressed;

532

•

Well-known result about Neumann series upper bound explored;

533

•

Bounds for realizations of stochastic system response proposed;

534

•

Accuracy and efficiency shown to depend on choice of matrix norm;

535

•

Good compromise between accuracy and efficiency obtained with Frobenius norm.

536 537 538