Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

57, 1, pp. 37-48, Warsaw 2019
DOI: 10.15632/jtam-pl.57.1.37

TIME INTEGRATION OF STOCHASTIC GENERALIZED EQUATIONS OF
MOTION USING SSFEM

Mariusz Poński

Czestochowa University of Technology, Department of Civil Engineering, Częstochowa, Poland

e-mail: mponski@bud.pcz.czest.pl

The paper develops an integration approach to stochastic nonlinear partial differential equ-
ations (SPDE’s) with parameters to be random fields. The methodology is based upon
assumption that random fields are from a special class of functions, and can be described
as a product of two functions with dependent and independent random variables. Such an
approach allows one to use Karhunen-Loève expansion directly, and themodified stochastic
spectral finite element method (SSFEM). It is assumed that a random field is stationary
and Gaussian while the autocovariance function is known. A numerical example of one-
dimensional heat waves analysis is shown.

Keywords: spectral stochastic finite element method, time integration, heat waves

1. Introduction

In the literature, one can find works describing SPDE’s solutions using SSFEM for stationa-
ry problems (Le Maitre and Knio, 2010; Matthies and Keese, 2005). There are also methods
handling transient problems such as theMonte Carlomethod (MC), perturbationmethod (Ka-
miński, 2013; Służalec, 2003), stochastic collocation method (SCM) (Acharjee and Zabaras,
2006; Babuška et al., 2007; Xiu andHesthaven, 2005). Thesemethods are widely used formany
problems, e.g. continuum mechanics, fluid dynamics, heat flow and have their advantages and
disadvantages (Stefanou, 2009; Xiu, 2010). As the leading in the literature, the SCMmethod is
listed due to the possibility of analyzing complex nonlinear problems. It reduces analysis time by
usingmultithreading and can be simply implemented (deterministic solver is treated as so-called
“black box”). Competitive to SCM is themethod of Stochastic Spectral Finite ElementMethod
proposed byGhanem and Spanos (2003). Thismethod is one of the so-called intrusivemethods,
which is very effective in solving linear problems, but requires building the source code from
scratch. An important disadvantage that is mentioned in many works is the coupling of equ-
ations that prevent the use of parallel solvers. This problemwas solved by applying the domain
decomposition method (Subber and Sarkar, 2014). Another important disadvantage mentioned
in theworks on numerical solutionwith SSFEMis the considerable difficulty of solving nonlinear
problems. Mathematical formulation of nonlinear stationary equations can be found in works
(Arregui-Mena et al., 2016; Ghosh et al., 2008; Hu et al., 2015;Matthies andKeese, 2005; Nouy,
2008; Nouy and LeMaitre, 2009; Stefanou et al., 2017; Xiu and Karniadakis, 2003; Zakian and
Khaji, 2016) rather than description of a general numerical approach. There is also no compre-
hensive solution to transient problems. This paper presents amethodology based upon a special
class of functions occurring in constitutive equations which can be described as a product of
two functions, respectively with dependent and independent variables. This approach allows one
to extend the applicability of SSFEM to solve wide range of nonstationary nonlinear stochastic
PDE’s.



38 M. Poński

2. Stochastic description

In this work, a modified SSFEM is used. This method is based upon notion of a random field.
The random field α(x,ω) (Acharjee and Zabaras, 2007; Xiu, 2010) (x ∈ D ⊂ R, ω ∈ Ω) is a
real valued measurable function which assigns a random variable α(ω) to each point x on a
fixed probability space (Ω,Z,P). HereΩ is the set of elementary events, Z is the σ-algebra and
P : Z → [0,1] is a probability measure. To obtain a computationally useful representation of
the process α(x,ω), it will be presented in the canonical form. Among various forms of such a
representation, a spectral representation -Karhunen-Loève expansionwill be adopted in further
considerations (Ghanem and Spanos, 2003). This expansion may be presented in the following
form

α(x,ω) = ξ0α(x)+
∞
∑

i=1

ξi(ω)
√

λifi(x) x∈ D, ω∈Ω (2.1)

Such a Karhunen-Loève expansion is truncated toM terms

α(x,ω)≈ ξ0α(x)+
M
∑

i=1

ξi(ω)
√

λifi(x) x∈ D, ω∈Ω (2.2)

In equation (2.1), {ξi(ω)}∞i=1 is a set of othonormal independent Gaussian random variables
with mean ξ0 = 1 and standard deviations equal to one, (·) is the expected value operator.
Constants {λi}∞i=1 and deterministic functions {fi(x)}∞i=1 are the eigenvalues and eigenfunctions
of the covariance kernel
∫

D

Ckernel(x1,x2)fi(x2) dx2 =λifi(x1) i∈ N = {1,2, . . .} (2.3)

The polynomial chaos (Ghanem and Spanos, 2003; LeMaitre andKnio, 2010) representation of
the random variableU(ω), truncated to P terms, can be written as

U(ω)≈U(ξ)=
P
∑

i=0

uiΦi(ξ) (2.4)

where {Φi(ξ)}Pi=0 denotes polynomial chaoses (Ghanem and Spanos, 2003), and {ui}Pi=0 are
coefficients of the expansion.The coefficientP describedby the expression (GhanemandSpanos,
2003)

P =1+
p
∑

s=1

2

s!

s−1
∏

r=0

(M+r) (2.5)

is the total number of polynomial chaoses used in the expansion, excluding the zero-th order
term,withpdenoting theorderof polynomial chaoses (detailed descriptionof polynomial chaoses
can be found in Ghanem and Spanos (2003), LeMaitre and Knio (2010)).

3. Governing equations formulation of SSFEM – an approach for a special class
of equations

3.1. Method of solution of the stochastic problem (SSFEM)

First step in solving SPDE’s, after finite element discretisation, is the stochastic randomfield
discretisation. Let K denotes the number of nodes of the discretized domain. The equation of



Time integration of stochastic generalized equations... 39

motion with given initial and boundary conditions, which is analyzed, can be written in a well
knownmatrix form (Bathe, 1996)

M(x,ω,q(t,ω))q̈(t,ω)+C(x,ω,q(t,ω))q̇(t,ω)+K(x,ω,q(t,ω))q(t,ω) =F(x,ω,q(t,ω))

(3.1)

where t ∈ T denotes time, x ∈ D × D × D def= D3 ⊂ R3 denotes space variables, q(t,ω) is
a generalized K× 1 displacement vector and q̇(t,ω), q̈(t,ω) is the first and second derivative
in time. Suppose that for the K ×K matrices M,C,K and K × 1 vector F separation of
dependent and independent variables can be made, which allows one to use Karhunen-Loève
expansion directly, e.g.

K(x,ω,q(t,ω)) =K(x,ω)fK(q(t,ω)) (3.2)

where fK(q(t,ω)) is a real valued Riemann integrable function on a suitable space andK(x,ω)
is a K×K matrix. Such an approach limits applicability of this method to problems in which
constitutive equations or nonlinear boundary conditions can bewritten as a product of functions
of independent and dependent variables, e.g. k(x,q(x, t)) = ka(x)kb(q(x, t)). Many physical
relations can be written in this way or can be reduced to either k(x,q(x, t)) = kb(q(x, t)) or
k(x,q(x, t)) = ka(x).
Assume that the discretized function of generalized displacement at each node has represen-

tations in polynomial chaos

qk(t,ω)≈
P
∑

i=0

(qspect(t))k,iΦi(ξ(ω)) k=1,2, . . . ,K (3.3)

(index k denotes node number) and let it be derived its first and second order derivative with
respect to the time

q̇k(t,ω)≈
P
∑

i=0

(q̇spect(t))k,iΦi(ξ(ω)) k=1,2, . . . ,K

q̈k(t,ω)≈
P
∑

i=0

(q̈spect(t))k,iΦi(ξ(ω)) k=1,2, . . . ,K

(3.4)

The vector of nodal generalized displacement can be written as

q(t,ω)≈q(t,ξ)=qmatrixspect (t)Φnode(ξ) (3.5)

where qmatrixspect (t) is aK× (P +1) matrix built from spectral coefficients and

Φnode(ξ)=
[

Φ0(ξ) Φ1(ξ) · · · ΦP(ξ)
]T

(3.6)

is a vector built from the polynomial chaoses.
Substituting appropriate derivatives of equation (3.5) to equation (3.1) and using (3.2),

(3.5)-(3.6), the following equation

M(x,ω)fM(q
matrix
spect (t)Φnode(ξ))q̈

matrix
spect (t)Φnode(ξ)

+C(x,ω)fC(q
matrix
spect (t)Φnode(ξ))q̇

matrix
spect (t)Φnode(ξ)

+K(x,ω)fK(q
matrix
spect (t)Φnode(ξ))q

matrix
spect (t)Φnode(ξ)=F(x,ω)fF(q

matrix
spect (t))

(3.7)

can be obtained.



40 M. Poński

Let us represent thematricesM(x,ω),C(x,ω),K(x,ω),F(x,ω) using theKarhunen-Loève
expansion, e.g.

K(x,ω)≈K(x,ξ)= (K(x))0ξ0+
M
∑

i=1

ξi(ω)(K
0(x))i (3.8)

where (·)0 denotes a matrix computed for the mean value of the process and (·)0 denotes a
matrix built of the shape functions and Karhunen-Loève expansion terms (e.g. (K0(x))i =
∫

D
∇N(
√
λifi(x))∇NT dx, where N is a vector built of test functions from the Sobolev space

H1(0, lelem) where lelem is a finite element length).
In order to formulate a suitable system of equations, let us represent thematrix of stochastic

eigenmodes of the solution as a vector qmatrixspect →qvectorspect where

qvectorspect (t)=
[[

q00(t) q01(t) · · · q0P(t)
] [

· · · qK0(t) · · · qKP(t)
]]T

(3.9)

where P, as before, is the total number of polynomial chaoses used in the expansion and K is
the total number of nodes in the FEM solution. The global vector of polynomial chaoses takes
the form

Φ(ξ)=
[

(Φnode(ξ))0 (Φnode(ξ))1 · · · (Φnode(ξ))K
]T

(3.10)

where the vectors of polynomial chaoses are the same for each node

(Φnode(ξ))0 =(Φnode(ξ))1 = . . .=(Φnode(ξ))K (3.11)

After substitution of appropriate expansions (3.8) of the matrices M(x,ω), C(x,ω), K(x,ω),
F(x,ω), equation (3.9) and (3.10) into (3.7), then multiplying by Φ(ξ) and averaging with
respect to the random space

〈

Φ(ξ)
(

(M(x))0ξ0+
MM
∑

iM=1

ξiM(M
0(x))iM

)

fMT

(

(qvectorspect (t))
TΦ(ξ)

)(

q̈vectorspect (t)
)T
Φ(ξ)

〉

+
〈

Φ(ξ)
(

(C(x))0ξ0+
MC
∑

iC=1

ξiC(C
0(x))iC

)

fCT

(

(qvectorspect (t))
TΦ(ξ)

)(

q̇vectorspect (t)
)T
Φ(ξ)

〉

+
〈

Φ(ξ)
(

(K(x))0ξ0+
MK
∑

iK=1

ξiK(K
0(x))iK

)

fKT

(

(qvectorspect (t))
TΦ(ξ)

)(

qvectorspect (t)
)T
Φ(ξ)

〉

=
〈

Φ(ξ)
(

(F(x))0ξ0+
MF
∑

iF=1

ξiF(F
0(x))iF

)

fF
(

(qvectorspect (t))
TΦ(ξ)

)〉

(3.12)

can be obtained, where

〈·〉=
∫

Ω

(·) dP(ω) (3.13)

andMM,MC,MK,MF denote the numbers of terms in the Karhunen-Loève expansion.
Finally, the set of (P +1)K nonlinear deterministic equations of the SSFEM method is

obtained

Mexpand(x,qvectorspect (t))q̈
vector
spect (t)+C

expand(x,qvectorspect (t))q̇
vector
spect (t)

+Kexpand(x,qvectorspect (t))q
vector
spect (t)=F

expand(x,qvectorspect (t))
(3.14)



Time integration of stochastic generalized equations... 41

where, for example, the generalized stiffness matrix (matrices Mexpand and Cexpand can be
written in the same way)

Kexpand(x,qvectorspect (t))

=
〈

Φ(ξ)
(

(K(x))0ξ0+
MK
∑

iK=1

ξi(K
0(x))iK

)

fK
(

(qvectorspect (t))
TΦ(ξ)

)

(Φ(ξ))T
〉 (3.15)

and the generalized force vector

Fexpand(x,qvectorspect (t))=
〈

Φ(ξ)
(

(F(x))0ξ0+
MF
∑

iF=1

ξi(F
0(x))iF

)

fF
(

(qvectorspect (t))
TΦ(ξ)

)〉

(3.16)

For fr(q
TΦ(ξ)) = 1, r = M,C,K,F equation (3.14) is linear. Moreover, the terms

〈ξiΦ(ξ)(Φ(ξ))T〉 and 〈ξiΦ(ξ)〉 of the above stated matrices obtained from integration over the
random space have a lot of zero entries (Ghanem and Spanos, 2003). In addition, this terms
may be determined in advance and only once.
In the proposed solution of nonlinear equation (3.14), the matrices in this equation ha-

ve to be numerically integrated both in the iteration step and time step, over a random and
geometric space. This is due to nonlinear functions fr (where r = M,C,K,F) of the depen-
dent variable q appearing in the parts 〈ξiΦ(ξ)fg(qTΦ(ξ))(Φ(ξ))T〉 (where g = M,C,K) and
〈ξiΦ(ξ)fF(qTΦ(ξ))〉 of these matrices. Methods of evaluating the above mentioned inner pro-
duct for different types of nonlinearities of the functions fr can be found in work of Le Maitre
and Knio (2010).
For a complete presentation of the stochastic process q(x, t,ω), the covariance matrix is

determined (Ghanem and Spanos, 2003)

Covqq(t)=
P
∑

j=0

〈Φj(ξ)Φj(ξ)〉(qvectorspect (t))j[(qvectorspect (t))j]T (3.17)

where Covqq = {(Covqq)i,j}Pi,j=1, the j-th eigenmode of the vector qvectorspect (t) has been written
as (qvectorspect (t))j. Therefore, the expected value may be calculated

E(qvectorspect (t))= (q
vector
spect (t))0 (3.18)

and the variance

Varq(t)i =Covqq(t)i,i i=1,2, . . . ,P (3.19)

Matrix equation (3.14) is a deterministic system of nonlinear equations, therefore, in order to
solve it, one ofmethods of direct integration, for exampleNewmarkmethod, can beused (Bathe,
1996).

3.2. Application of SSFEM to non-classical stochastic heat conduction constitutive model
– heat waves analysis

Themostwidely usedmodel formany engineering problems is the classic equation of Fourier
(Fourier, 1822), which can be represented as a function of a random field. The expression for
the heat flux can be written

qF(x, t,ω) =−kF(x,ω,T)∇T(x, t,ω) x∈ D3 ⊂ R3, ω∈Ω, t∈ T (3.20)

where kF(x,ω,T) is Fourier thermal conductivity.



42 M. Poński

Because of theanomalies associatedwith theFouriermodel (Vernotte, 1958) andthepresence
of finite speedpropagation of heat, aCattaneomodel that takes into account heat flux relaxation
has been introduced (Cattaneo, 1948; Służalec, 2003)

τ∂tqC(x, t,ω)+qC(x, t,ω) =−kC(x,ω,T)∇T(x, t,ω) x∈ D3 ⊂ R3, ω∈Ω, t∈ T
(3.21)

where τ represents the relaxation time and kC(x,ω,T) is the Cattaneo thermal conductivity.
TheJeffreys typemodel is another heat conduction constitutivemodel ofwhich theCattaneo

model and a Fourier-like diffusive model are subcases that can be obtained from this model
(Joseph and Preziosi, 1989; Straughan, 2011; Tamma and Zhou, 1998; Ván and Fülöp, 2012)

τ∂tq(x, t,ω)+q(x, t,ω) =−k(x,ω,T)[∇T(x, t,ω)−K(x,ω,T)∂t(∇T(x, t,ω))] (3.22)

where x∈ D3 ⊂ R3, ω∈Ω, t∈ T, and

k(x,ω,T) = kF(x,ω,T)+kC(x,ω,T) (3.23)

and

K(x,ω,T) =
τkF(x,ω,T)

k(x,ω,T)

is the so-called retardation time.
When the retardation time, K = 0, the Jeffreys model is reduced to the Cattaneo model.

When selecting, K = τ, the Jeffreys model only degenerates to a Fourier-like diffusive model
with relaxation (Tamma and Zhou, 1998).
Zhou and co-workers (Tamma and Zhou, 1998) introduced C-process and F-process models

which are a linear combination of the Fourier andCattaneomodels. The basic assumption is the
simultaneous occurrence of a fast process based on equation (3.21) and a slow process related
to equation (3.22). This model is a generalization of the above stated relations. The equations
describing connection between heat flux and temperature have the following form (Tamma and
Zhou, 1998)

qCF(x, t,ω) =qF(x, t,ω)+qC(x, t,ω) x∈ D3 ⊂ R3, ω∈Ω, t∈ T (3.24)

where

qC(x, t,ω)+ τ∂tqC(x, t,ω)=−(1−FT(x,ω,T))k(x,ω,T)∇T(x, t,ω)
qF(x, t,ω) =−FT(x,ω,T)k(x,ω,T)∇T(x,ω)

FT(x,ω,T) =
kF(x,ω,T)

kF(x,ω,T)+kC(x,ω,T)

k(x,ω,T) = kF(x,ω,T)+kC(x,ω,T)

(3.25)

After substitution of equation (3.25)1,2 into (3.24), it can be obtained

qCF(x, t,ω)+ τ∂tqCF(x, t,ω)

=−[k(x,ω,T)∇T(x, t,ω)+ τ∂t(k(x,ω,T)FT(x,ω,T)∇T(x, t,ω))]
(3.26)

IndexesF andC in the conductivity coefficient and in the heat flux vector respectively refer to
the Fouriermodel with infinite propagation speed of wave and to theCattaneomodel, occurring
simultaneously. Also the model number FT ∈ [0,1] has been introduced. It can be seen that
for FT ∈ [0,1] the Jeffrey model is obtained, for FT = 1 the Fourier one, and for FT = 0 the
equation is reduced to the Cattaneo model.



Time integration of stochastic generalized equations... 43

4. Stochastic nonlinear 1D transport equation

Toobtain equations describing the flowof heat in a rigid conductor, the following energy balance
equation (for clarity the shorthand notation is adopted T ≡T(x,t,ω))

C(x,ω,T)∂tT +∇·qCF =0 x∈ D, ω∈Ω, t∈ T (4.1)

whereC(x,ω,T) = ρc(x,ω,T) is combined with one of the model constitutive equation, initial
conditionT(x,t=0,ω)=T0 and a suitable boundary condition on ∂D. The equation associated
with C- and F-process models (which can be reduced to the above mentionedmodels) is stated
below

τ∂t(C(x,ω,T)∂tT)+C(x,ω,T)∂tT −∂x(kF(x,ω,T)∂xT)− τ∂t(∂x(kF(x,ω,T)∂xT))
−∂x((1−FT(x,ω,T))k(x,ω,T)∂xT)= 0

(4.2)

5. Governing equations formulation of the Galerkin Finite Element Method
(GFEM) – C- and F-processes model

Thefirst step to solve the stochastic problem is discretisation of a deterministic space by the use
of the finite element method in the Galerkin approach (Bathe, 1996). To this end, the response
in terms of the temperature field is approximated by the expression

T(x, t,ω) = (T(t,ω))TN(x) (5.1)

whereN(x) is a vector built of the test functions (as defined in Section 3.1), T(t,ω) is a vector
representing the discretized temperature field and the upper indexTdenotes transposition. The
final equations after appropriate transformations read

M(T)T̈+C(T)Ṫ+K(T)T=F(T) (5.2)

where Ṫ and T̈ denote the first and second order temperature derivative with respect to time.
The individualmatrices in equation (5.2) canbeobtained throughastandardFEmethod(Bathe,
1996).

6. Numerical example

Distributions of temperature statistical moments as functions of time of the considered model
for a thin steel sheet (similar to the work (Al-Nimr, 1997)) (Fig. 1) heated by a sudden heat
impulse (Ván and Fülöp, 2012) will be analyzed. For this purpose and for further analysis, the
following data will be adopted (Joseph and Preziosi, 1989; Ván and Fülöp, 2012).

Pulsed heating can bemodeled as an internal heat source (Bargmann andFavata, 2014) with
various time characteristics or as external boundary conditions. Let the heating function for the
mixed boundary condition (convection-radiation) takes the form

Tpulse(t)=T0+asin(bt)exp(−ct) (6.1)

with parameters T0 =293.15K, a=15 ·104K, b=π/10, c=50.
One dimensional region of the sample (modeled as a bar) is divided into 20 elements. The

adoptedheating time is equal to tmax =0.6·10−10 swith60 timesteps (time increment∆t=0.01·
10−10). Cattaneo thermal conductivity has been considered as a random function independent



44 M. Poński

Fig. 1. Schematic of a thin steel plate heated by a pulse

of temperature (stationary and Gaussian process), kC(x,ω,T) = kC(x,ω) with the covariance
kernel

Ckernel(x1,x2)=σ
2
kC
exp
(−|x1−x2|

b

)

(6.2)

where the coefficient of variation σ2kC and correlation length b are stated in Table 1.

Table 1.Parameters used in analysis

Parameter Value

Heat capacity per unit volume
c=434.0J/(kgK)

(deterministic)

Density (deterministic) ρ=7850.0kg/m3

Fourier thermal conductivity
kF =54W/(mK)

(deterministic)

Cattaneo thermal conductivity 〈kC(x,ω)〉=1210.0 ·1010W/(mKs)
(random –mean)

Cattaneo thermal conductivity σ2kC =2500.0 ·10
10W/(mKs)

(random – coefficient of variance)

Relaxation time τ =20.0 ·10−12 s
Heat convection coefficient αc =9.0W/(m

2K)
(deterministic)

Emissivity (deterministic) εr =0.625

Stefan-Boltzmann constant σB =5.67 ·10−8W/(m2K4)
Thickness L=0.005m

Correlation length b=0.001m

Using the Karhunen-Loève expansion, the Cattaneo thermal conductivity can be written in
the form

kC(x,ω)≈〈kC(x,ω)〉+
M
∑

i=1

ξi(ω)
√

λifi(x) (6.3)



Time integration of stochastic generalized equations... 45

Therefore, the randommatrix present in equation (3.14) can be expressed as

F
expand
x=0 (x=0,T

vector
spect (t))

=
〈

Φ(ξ)
(

(Fx=0(x=0))0ξ0+
MF
∑

iF=1

ξi(F
0(x=0))iF

)

fF
(

(Tvectorspect (t))
TΦ(ξ)

)〉 (6.4)

This matrix takes non-zero values for the boundary node (x = 0). The matrices (F(x))0 = 0,
(F0(x))iF = 0 are generally equal to zero, only for x = 0 (Fx=0(x = 0))0 = I. For the mixed
boundary condition (convection-radiation) it can be assumed that (parameters in Table 1)

fF
(

(Tvectorspect (t))
TΦ(ξ)

)

=αc
(

Tpulse(t)−
(

(Tvectorspect (t))
TΦ(ξ)

))

+εrσB
[

Tpulse(t)
4−
(

(Tvectorspect (t))
TΦ(ξ)

)4] (6.5)

Because the Cattaneo thermal conductivity is independent of the temperature function fK
included in the matrix

Kexpand(x,Tvectorspect (t))

=
〈

Φ(ξ)
(

(K(x))0ξ0+
MK
∑

iK=1

ξiK(K
0(x))iK

)

fK
(

(Tvectorspect (t))
TΦ(ξ)

)

(Φ(ξ))T
〉 (6.6)

(6.6) can bewritten as fK(·)= 1, and thematrices (K(x))0, (K0(x))iK can be determined from

(K(x))0 =

∫

D

∇NkF∇NT dx+
∫

D

∇N〈kC(x,ω)〉∇NT dx

(K0(x))iK =

∫

D

∇N(
√

λiKfiK(x))∇N
T dx

(6.7)

the proposedmodified SSFEMhas been compared to theMonte Carlomethod using theC- and
F-processmodels. As the relevant set of points, the heated surface has been chosen. As shown in
Fig. 2, there is a good correlation between themethods. SSFEM is giving smaller values for the

Fig. 2. Temperature mean value and relative error between SSFEM (3rd order polynomial chaos,
2nd order Karhunen-Loève expansion) andMonte Carlo (5000 samples) solution in function of time at

the heated surface



46 M. Poński

standard deviation than MC (Fig. 3), which is typical for this method (Ghanem and Spanos,
2003). It can be seen that the biggest relative error for the standard deviation occurs in time
nodes with smallest values. The relative error of temperature (Fig. 2) between the mean value
obtained from the SSFEMandMonte Carlo solution is small with the extremum not exceeding
2.5 ·10−8. Computations have been performed for P =3 order of polynomial chaos andM =2
order of Karhunen-Loève expansion.

Fig. 3. Temperature standard deviation and relative error between SSFEM (3rd order polynomial chaos,
2nd order Karhunen-Loève expansion) andMonte Carlo (5000 samples) solution in function of time at

the heated surface

Fig. 4. Temperature eigenmodes obtained from SSFEM solution in function of time at the heated
surface – 1st and 2nd order of polynomial chaos, 2nd order Karhunen-Loève expansion

Figure 4 illustrates solutions for successive orders of polynomial chaos. As can be seen, the
biggest influence is the first order of the expansion on the basis of which it can be concluded that
the improvement of the convergence of statistical moments by increasing the order of expansion
in polynomial chaos will be small.



Time integration of stochastic generalized equations... 47

7. Concluding remarks

The paper develops an approach to analysis stochastic nonlinear partial differential equations
(SPDE’s). As an example of stochastic analysis, the heat waves equation has been shown. The
C-F-processes constitutive model has been chosen for the analysis. It can be reduced to the
Fourier, Cattaneo and Jeffery types of models. Amodified SSFEM, which consists in the sepa-
ration of dependent and independent variables in themainmatrices, has been proposed to solve
nonlinear governing equations. The modified SSFEM has been compared to the Monte Carlo
method. The comparison has shown that the proposed method works with nonlinear problems
well and for equation (4.1) the solutions generated by the SSFEM method are convergent to
solutions generated by theMonte Carlomethod due to the first and second statistical moment.
The analysis has revealed that the largest difference in the results obtained from the SSFEM
andMCmethod is generated in time nodeswith the smallest standard deviation (localminima).
Comparison of the results from themethods has aimed at demonstrating compliance rather than
efficiency or time consumption of the SSFEM. In order to check time consumption of the me-
thod in relation toMCor SCM, one should use the domain decompositionmethod (Subber and
Sarkar, 2014) andmethods of reducing the integration time of the main matrices (e.g. Smolyak
sparse grid method (Smolyak, 1963)) which would allow one to use parallel processing.

References

1. Acharjee S., Zabaras N., 2006, Uncertainty propagation in finite deformations – a spectral
stochastic Lagrangian approach,Computer Methods in Applied Mechanics and Engineering, 195,
19, 2289-2312

2. Acharjee S., Zabaras N., 2007, A non-intrusive stochastic Galerkin approach for modeling
uncertainty propagation in deformation processes, Computational Stochastic Mechanics, 85, 5,
244-254

3. Al-Nimr, M., 1997, Heat transfer mechanisms during short-duration laser heating of thin metal
films, International Journal of Thermophysics, 18, 5, 1257-1268

4. Arregui-Mena J.D., Margetts L., Mummery P.M., 2016, Practical application of the sto-
chastic finite element method,Archives of Computational Methods in Engineering, 23, 1, 171-190

5. Babuška I., Nobile F., Tempone R., 2007, A stochastic collocation method for elliptic par-
tial differential equations with random input data, SIAM Journal on Numerical Analysis, 45, 3,
1005-1034

6. Bargmann S., Favata A., 2014, Continuummechanical modeling of laser-pulsed heating in po-
lycrystals: a multi-physics problem of coupling diffusion, mechanics, and thermal waves,ZAMM –
Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Me-
chanik, 94, 6, 487-498

7. Bathe K.J., 1996,Finite Element Procedures, Prentice Hall, New Jersey

8. Cattaneo M., 1948, Sulla Conduzione de Calor,Mathematics and Physics Seminar, 3, 3, 83-101

9. Fourier J., 1822,Theorie Analytique de la Chaleur, Chez Firmin Didot, Paris

10. Ghanem R.G., Spanos P.D., 2003, Stochastic Finite Elements: A Spectral Approach, Dover
Publications, Mineola, USA

11. Ghosh D., Avery P., Farhat C., 2008, A method to solve spectral stochastic finite element
problems for large-scale systems, International Journal for Numerical Methods in Engineering, 00,
1-6

12. Hu J., Jin S., Xiu„ D., 2015, A stochastic Galerkinmethod for Hamilton-Jacobi equations with
uncertainty, SIAM Journal on Scientific Computing, 37, 5, A2246-A2269



48 M. Poński

13. Joseph D.D., Preziosi L., 1989, Heat waves,Reviews of Modern Physics, 61, 1, 41-73

14. KamińskiM., 2013,The Stochastic PerturbationMethod for Computational Mechanics, JohnWi-
ley & Sons, Chichester

15. Le Maitre O.P., Knio O.M., 2010, Spectral Methods for Uncertainty Quantification: with Ap-
plications to Computational Fluid Dynamics, Springer Science &BusinessMedia, Doredrecht

16. Matthies H.G., Keese A., 2005, Galerkin methods for linear and nonlinear elliptic stochastic
partial differential equations, Computer Methods in Applied Mechanics and Engineering, 194, 2,
1295-1331

17. Nouy A., 2008, Generalized spectral decomposition method for solving stochastic finite element
equations: invariant subspace problem and dedicated algorithms, Computer Methods in Applied
Mechanics and Engineering, 197, 51, 4718-4736

18. Nouy A., Le Maitre O.P., 2009, Generalized spectral decomposition for stochastic nonlinear
problems, Journal of Computational Physics, 228, 1, 202-235

19. Służalec A., 2003, Thermal waves propagation in porous material undergoing thermal loading,
International Journal of Heat and Mass Transfer, 46, 9, 1607-161

20. Smolyak S.A., 1963,Quadrature and interpolation formulas for tensor products of certain classes
of functions,Doklady Akademii Nauk SSSR, 4, 240-243

21. Stefanou G., 2009, The stochastic finite element method: past, present and future, Computer
Methods in Applied Mechanics and Engineering, 198, 9, 1031-1051

22. Stefanou G., Savvas D., Papadrakakis M., 2017, Stochastic finite element analysis of com-
posite structures based on mesoscale random fields of material properties, Computer Methods in
Applied Mechanics and Engineering, 326, 319-337

23. Straughan B., 2011,Heat Waves, Springer, NewYork

24. SubberW., SarkarA., 2014,A domain decompositionmethod of stochasticPDEs:An iterative
solution techniques using a two-level scalable preconditioner, Journal of Computational Physics,
257, 298-317

25. Tamma K.K., Zhou X., 1998, Macroscale and microscale thermal transport and thermo-
-mechanical interactions: some noteworthy perspectives, Journal of Thermal Stresses, 21, 3-4,
405-449

26. Ván P., Fülöp T., 2012, Universality in heat conduction theory: weakly nonlocal thermodyna-
mics,Annalen der Physik, 524, 8, 470-478

27. Vernotte P., 1958, Les paradoxes de la theorie continue de l’equation de la chaleur, Comptes
Rendus, 246, 3154-3155

28. XiuD., 2010,NumericalMethods for Stochastic Computations: A SpectralMethodApproach, Prin-
ceton University Press, Princeton

29. Xiu D., Hesthaven J.S., 2005, High-order collocation methods for differential equations with
random inputs, SIAM Journal on Scientific Computing, 27, 3, 1118-1139

30. Xiu D., Karniadakis G.E., 2003, A new stochastic approach to transient heat conduction mo-
deling with uncertainty, International Journal of Heat and Mass Transfer, 46, 24, 4681-4693

31. Zakian P., Khaji N., 2016, A novel stochastic-spectral finite element method for analysis of
elastodynamic problems in the time domain,Meccanica, 51, 4, 893-920

Manuscript received January 31, 2018; accepted for print June 17, 2018