# A Pseudospectral Approach for Kirchhoff Plate Bending Problems with Uncertainties

Authors: Ling Guo; Jianguo Huang
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

This paper proposes a pseudospectral approach for the Kirchhoff plate bending problem with uncertainties. The Karhunen-Loève expansion is used to transform the original problem to a stochastic fourth-order PDE depending only on a finite number of random variables. For the latter problem, its exact solution is approximated by a gPC expansion, with the coefficients obtained by the sparse grid method. The main novelty of the method is that it can be carried out in parallel directly while keeping the high accuracy and fast convergence of the gPC expansion. Several numerical results are performed to show the accuracy and performance of the method.


---

## Body

## 1. Introduction

Elastic plates are basic elements in structural mechanics, which are frequently used in aviation, aerospace, mechanical structures, civil engineering, and so on [1]. To date, there have developed many computational methods, including finite difference, finite element, and boundary element, to numerically solve their deformation in the deterministic setting. However, many physical parameters and quantities involved, such as the Poisson ratio, Young’s modulus, and applied loads, are uncertain in real-world applications. Hence, it is very important to investigate the deformation of elastic plates in the stochastic setting. Depending on the amount of information available, uncertainties can be described in several ways, for example, fuzzy set theory, worst-case scenario analysis, and probabilistic setting. We refer the reader to [2] and the references therein for details along this line. In this paper, we focus on probabilistic description of uncertainties, where the uncertain input data is modeled as random variables or random field with a given spatial correlation structure. Thus, the mathematical models of plates with uncertainties are governed by stochastic differential equations with random input data.As far as we know, great efforts have been made to devise computational methods for plate structures with uncertainties over the past several decades. Monte Carlo sampling (MCS) is the most commonly used method. Since the approach only requires repetitive executions of deterministic solvers, it is very easy to implement in real applications. The main limitation of MCS is its slow convergence [3]. The most popular nonsampling methods is the perturbation methods, where random fields are expanded via Taylor series around their mean and truncated at certain order. Such technique has been extensively used to obtain the random transversal displacements of plates with random properties [4–8]. However, these methods are efficient only when the magnitude of uncertainties is not very large.A recently developed method involving generalized polynomial chaos (gPC) has now become one of the most widely used methods for solving stochastic partial differential equations. gPC started with the seminal work on PC (Polynomial Chaos) by Ghanem and coworkers. Motivated by the theory of Winer-Hermite homogeneous chaos, Spanos, and Ghanem employed Hermite polynomials as orthogonal basis to represent random process, leading to a great success in numerical solutions of many engineering problems [9]. To alleviate the difficulty in the use of Hermite polynomials for non-Gaussian process, the generalized polynomial chaos (gPC) was proposed by Xiu and Karniadakis in [10]. In the gPC setting, relevant orthogonal polynomials are taken as basis functions in view of the probability distribution of random inputs. To apply the method for solving stochastic differential equations with random inputs, the key step is to determine the expansion coefficients of the gPC expansion of the solution. In [11], da Silva and Beck used gPC and the Galerkin method, named stochastic Galerkin method, to obtain the stochastic displacement response of Kirchhoff plates on Winkler foundations. However, the use of rapidly growing number of basis functions will deteriorate the efficiency of the stochastic Galerkin methods significantly.Very recently, there is a surge of interests in high-order stochastic collocation method (SC) following the works of [12]. By introducing “sparse grid” from multivariate interpolation analysis, the SC method enjoys the power of Monte Carlo methods as well as Stochastic Galerkin methods. It was shown in [13] (see also [14]) that SC converges exponentially fast with respect to the number of points employed for each random input variable. Moreover, this method leads to the solution of uncoupled deterministic problems, just like MCS, even when the input data depend nonlinearly on the driving random variables. Therefore, the solution process is highly parallelizable. In the earlier high-order SC formulation, the basis functions are Lagrange polynomials defined by nodes. A more practical “pseudospectral” approach was suggested in [15], which is easier to manipulate in practice than the Lagrange interpolation approach.In this paper, we are going to apply the pseudospectral approach mentioned above to obtain approximate solutions of stochastic Kirchhoff plates with uncertainties. To this end, we seek an approximate solution in the form of gPC expansion and evaluate the gPC coefficients by solving deterministic plate bending problems defined over a set of collocation nodes-sparse grids built upon either Clenshaw-Curtis abscissas and the Gaussian abscissas. On each collocation node, the deterministic plate bending problem is computed by the Morley element method. In order to assess the efficiency of our method proposed, several numerical results are performed. The convergence rate is discussed inL2(𝒟) norm on different levels of sparse grids, and the first- and second-order moments of the approximated displacement obtained by our method are also compared with the ones obtained by Monte Carlo simulation.The rest of the paper is organized as follows. In Section2, the mathematical model for the stochastic Kirchhoff plate with uncertainties is given. Then we briefly review the gPC expansion and give the pseudospectral approach in Section 3. Some numerical results are performed to illustrate the efficiency of the method in the last section.


## 2. Mathematical Formulation of the Problem

### 2.1. Formulation of the Stochastic Kirchhoff Plate Problem

Let𝒟 be a bounded polygonal domain in ℝ2, in which every point is assigned by Cartesian coordinates x={x1,x2}T. Let (Ω,ℱ,P) be a complete probability space, where Ω is the sample space, ℱ⊂2Ω the minimal σ-algebra of elementary events, and P:ℱ→[0,1] a probability measure over ℱ. When there is no confusion caused, we will simply write Ω for (Ω,ℱ,P). Then the clamped Kirchhoff plate bending problem in the stochastic setting is to find a deflection field u:Ω×𝒟→ℝ such that the following equations hold in the sense of P, almost everywhere in Ω:
(2.1)-MIJ,IJ(u(ω,x))=f(ω,x)inD,u(ω,x)=∂nu(ω,x)=0on⁡∂D,
where n denotes the unit outward normal to ∂𝒟,(2.2)MIJ(u(ω,x)):=D(ω,x)((1-ν)KIJ(u(ω,x))+νKLL(u(ω,x))δIJ),KIJ(u(ω,x)):=-∂IJu(ω,x)=-∂2u(ω,x)∂xI∂xJ,I,J=1,2,
with D(ω,x) being the rigid flexibility of the plate and δIJ being the usual Kronecker delta, and f(ω,x) denotes the load force. We mention in passing that we use Einstein’s summation convention, throughout the paper, whereby the summation is implied when a subscript I, J, or L appears exactly twice.To ensure the existence and uniqueness of the solution, we need the following assumptions [16], that is:(2.3)∃D̲,D¯∈(0,+∞)suchthatP(ω∈Ω:D(ω,x)∈[D̲,D¯]∀x∈D¯)=1,
and the right-hand side f in (2.1) satisfies(2.4)∫Ω∫Df2(ω,x)dxdP(ω)<+∞which implies∫Df2(ω,x)dx<+∞a.s.


### 2.2. Representation of the Uncertainty

To solve (2.1) numerically, one of the most important steps is to parameterize the input uncertainty by a set of finite independent random variables. Such a procedure is often achieved via a certain type of decomposition which can approximate the target random process with desired accuracy. One of the choices is the Karhunen-Loève type expansion [17], which is based on the spectral decomposition of the covariance function of the input random process. Following such decomposition and assuming that the random inputs can be characterized by N random variables, we can write the random inputs in abstract form, for example,(2.5)D(ω,x)=DN(Y1(ω),…,YN(ω),x),f(ω,x)=fN(Y1(ω),…,YN(ω),x),
where N is a certain positive integer, {Yi}i=1N are real-valued random variables with zero mean value and unit variance. Hence, following the Doob-Dynkin lemma [18], the solution u of the plate bending problem (2.1) can be described by the same set of random variables {Yi(ω)}i=1N, that is, u(ω,x)=uN(Y1(ω),…,YN(ω),x). By characterizing the probability space with N random variables, we have effectively reduced the infinite-dimensional probability space to an N-dimensional space.In this paper, we take the customary approach (see, e.g, [16, 19]) of assuming that the random inputs are already characterized by a set of mutually independent random variables with satisfactory accuracy. Let us now assume that {Yi}i=1N are independent random variables with probability density functions ρi:Γi→ℝ+ and their images Γi:=Yi(Ω) are bounded intervals in ℝ for i=1,…,N. Then(2.6)ρ(y)=∏i=1Nρi(Yi)∀y∈Γ,
is the joint probability density function of y=(Y1(ω),…,YN(ω)) with the support(2.7)Γ:=∏i=1NΓi⊂RN.This allows us to rewrite the stochastic plate bending problem (2.1) as an (N+2)-dimensional differential equation in a strong form
(2.8)-MIJ,IJ(uN(y,x))=fN(y,x)∀(y,x)∈Γ×D,uN(y,x)=∂nuN(y,x)=0∀(y,x)∈Γ×∂D,
where N is the dimension of the random space Γ.Similar to the deterministic problems, we can define a (finite-dimensional) subspaceVΓ⊂Lρ2(Γ), the space of all square integrable function in Γ with respect to the measure ρ(y)dy, and find uNV(y,x)∈VΓ(y) in the following weak form:
(2.9)∫Γ-ρ(y)MIJ,IJ(uNV(y,x))v(y)dy=∫Γρ(y)fNv(y)dy∀v(y)∈VΓ,x∈D,∫Γρ(y)uNV(y,x)v(y)dy=∫Γρ(y)∂nuNV(y,x)v(y)dy=0∀v(y)∈VΓ,x∈D.A typical approach to getting numerical solution via (2.9) is to employ a stochastic Galerkin approach, which was investigated in [11]. The stochastic Galerkin methods offer fast convergence when the solution is sufficiently smooth in the random space. However, if N is very large, the rapidly growing number of basis functions would deteriorate the efficiency of the method considerably.


## 2.1. Formulation of the Stochastic Kirchhoff Plate Problem

Let𝒟 be a bounded polygonal domain in ℝ2, in which every point is assigned by Cartesian coordinates x={x1,x2}T. Let (Ω,ℱ,P) be a complete probability space, where Ω is the sample space, ℱ⊂2Ω the minimal σ-algebra of elementary events, and P:ℱ→[0,1] a probability measure over ℱ. When there is no confusion caused, we will simply write Ω for (Ω,ℱ,P). Then the clamped Kirchhoff plate bending problem in the stochastic setting is to find a deflection field u:Ω×𝒟→ℝ such that the following equations hold in the sense of P, almost everywhere in Ω:
(2.1)-MIJ,IJ(u(ω,x))=f(ω,x)inD,u(ω,x)=∂nu(ω,x)=0on⁡∂D,
where n denotes the unit outward normal to ∂𝒟,(2.2)MIJ(u(ω,x)):=D(ω,x)((1-ν)KIJ(u(ω,x))+νKLL(u(ω,x))δIJ),KIJ(u(ω,x)):=-∂IJu(ω,x)=-∂2u(ω,x)∂xI∂xJ,I,J=1,2,
with D(ω,x) being the rigid flexibility of the plate and δIJ being the usual Kronecker delta, and f(ω,x) denotes the load force. We mention in passing that we use Einstein’s summation convention, throughout the paper, whereby the summation is implied when a subscript I, J, or L appears exactly twice.To ensure the existence and uniqueness of the solution, we need the following assumptions [16], that is:(2.3)∃D̲,D¯∈(0,+∞)suchthatP(ω∈Ω:D(ω,x)∈[D̲,D¯]∀x∈D¯)=1,
and the right-hand side f in (2.1) satisfies(2.4)∫Ω∫Df2(ω,x)dxdP(ω)<+∞which implies∫Df2(ω,x)dx<+∞a.s.


## 2.2. Representation of the Uncertainty

To solve (2.1) numerically, one of the most important steps is to parameterize the input uncertainty by a set of finite independent random variables. Such a procedure is often achieved via a certain type of decomposition which can approximate the target random process with desired accuracy. One of the choices is the Karhunen-Loève type expansion [17], which is based on the spectral decomposition of the covariance function of the input random process. Following such decomposition and assuming that the random inputs can be characterized by N random variables, we can write the random inputs in abstract form, for example,(2.5)D(ω,x)=DN(Y1(ω),…,YN(ω),x),f(ω,x)=fN(Y1(ω),…,YN(ω),x),
where N is a certain positive integer, {Yi}i=1N are real-valued random variables with zero mean value and unit variance. Hence, following the Doob-Dynkin lemma [18], the solution u of the plate bending problem (2.1) can be described by the same set of random variables {Yi(ω)}i=1N, that is, u(ω,x)=uN(Y1(ω),…,YN(ω),x). By characterizing the probability space with N random variables, we have effectively reduced the infinite-dimensional probability space to an N-dimensional space.In this paper, we take the customary approach (see, e.g, [16, 19]) of assuming that the random inputs are already characterized by a set of mutually independent random variables with satisfactory accuracy. Let us now assume that {Yi}i=1N are independent random variables with probability density functions ρi:Γi→ℝ+ and their images Γi:=Yi(Ω) are bounded intervals in ℝ for i=1,…,N. Then(2.6)ρ(y)=∏i=1Nρi(Yi)∀y∈Γ,
is the joint probability density function of y=(Y1(ω),…,YN(ω)) with the support(2.7)Γ:=∏i=1NΓi⊂RN.This allows us to rewrite the stochastic plate bending problem (2.1) as an (N+2)-dimensional differential equation in a strong form
(2.8)-MIJ,IJ(uN(y,x))=fN(y,x)∀(y,x)∈Γ×D,uN(y,x)=∂nuN(y,x)=0∀(y,x)∈Γ×∂D,
where N is the dimension of the random space Γ.Similar to the deterministic problems, we can define a (finite-dimensional) subspaceVΓ⊂Lρ2(Γ), the space of all square integrable function in Γ with respect to the measure ρ(y)dy, and find uNV(y,x)∈VΓ(y) in the following weak form:
(2.9)∫Γ-ρ(y)MIJ,IJ(uNV(y,x))v(y)dy=∫Γρ(y)fNv(y)dy∀v(y)∈VΓ,x∈D,∫Γρ(y)uNV(y,x)v(y)dy=∫Γρ(y)∂nuNV(y,x)v(y)dy=0∀v(y)∈VΓ,x∈D.A typical approach to getting numerical solution via (2.9) is to employ a stochastic Galerkin approach, which was investigated in [11]. The stochastic Galerkin methods offer fast convergence when the solution is sufficiently smooth in the random space. However, if N is very large, the rapidly growing number of basis functions would deteriorate the efficiency of the method considerably.


## 3. Pseudospectral Approach

In collocation methods, we require the solution to satisfy (2.8) only at a discrete set of points in the corresponding random space. Unlike Monte Carlo sampling, our focus here is on the approaches that utilize polynomial approximation theory to strategically locate the nodes to gain accuracy.


### 3.1. Generalized Polynomial Chaos Expansion

In the finite-dimensional random spaceΓ, the gPC expansion seeks to approximate a random function via orthogonal polynomials of random variables. Let us define one-dimensional orthogonal polynomial spaces with respect to the measure ρi(Yi)dYi in Γi,(3.1)Wi,mi:={v:Γi→R:v∈span⁡{ϕm(Yi)}m=0mi},i=1,…,N,
where {ϕm(Yi)} are a set of orthogonal polynomials satisfying the orthogonality conditions(3.2)∫Γiρi(Yi)ϕm(Yi)ϕn(Yi)dYi=hm2δmn,
with normalization factors hm2=∫Γiρiϕm2dYi. With proper scaling, one can always normalize the bases such that hm2≡1 and this will be adopted throughout this paper. The probability density function ρi(Yi) in the above orthogonality relation (3.2) serves as a role of integration weight, which in turn defines the type of orthogonal polynomials {ϕm}. We refer this to [10, 14] for more details.The correspondingN-variate orthogonal polynomial space in Γ is defined by(3.3)WNP=⨂m∈ΛWi,mi,
where the tensor product is taken over all possible combinations of the multi-index m=(m1,…,mN) satisfying m∈Λ:={α:|α|=∑i=1Nmi≤P}. Thus, WNP is the space of N-variate orthogonal polynomials of total degree at most P, and the number of basis polynomials in (3.3) is dim(WNP)=(P+N)!/P!N!. Let {Φm(y)} be the N-variate orthonormal polynomials from WNP. They are constructed as products of a sequence of univariate polynomials in each direction of Yi, i=1,…,N, that is,(3.4)Φm(y)=∏m∈Λϕmi(Yi),
where mi is the order of the polynomial ϕ(Yi) in Yi for 1≤i≤d. For the N-variate polynomial Φm, each index 1≤m≤M=dim(WNP) corresponding to a unique sequence (m1,…,mN)∈Λ; they are also satisfying the orthogonality conditions(3.5)∫ΓΦm(y)Φn(y)ρ(y)dy=δmn,1≤m,n≤M.The finite-order, Pth-order, gPC approximation of the stochastic transverse displacementuN(y,x) in (2.8) is(3.6)PNPu:=uNP(y,x)=∑m=1Mûm(x)Φm(y),M=(N+PN),
where ℙNP denotes the orthogonal projection operator from Lρ2(Γ) onto WNP and {ûm} are the Fourier coefficients defined by(3.7)ûm(x)=∫ΓuN(y,x)Φm(y)ρ(y)dy=E[uN(y,x)Φm(y)].When a sufficient accurate gPC approximation (3.6) is available, one has in fact an analytical representation of uN(y,x) in terms of the random inputs y. Therefore, practically all statistical information can be retrieved in a straightforward manner. For example, the mean solution of the transverse displacement is(3.8)E(uN(y,x))≈E(uNP(y,x))=∫(∑m=1MûmΦm(y))ρ(y)dy=û1.
The variance of the solution can be obviously approximated as(3.9)Var⁡(uN(y,x))≈E[(uNP(y,x)-E(uNP(y,x)))2]=∑m=2M[ûm2].


### 3.2. Pseudospectral Approach

In the pseudospectral approach, we again seek an approximate solution in the form of gPC expansion, similar to (3.6), that is, for any x∈𝒟,(3.10)INPu:=wNP(y,x)=∑m=1Mŵm(x)Φm(y),
where 𝕀NP denotes the other projection operator from Lρ2(Γ) onto WNP and the expansion coefficients are determined as(3.11)ŵm=∑j=1QuN(yj,x)Φm(yj)αj,m=1,…,M.
Here {yj,αj}j=1Q are a set of nodes and weights, where yj=(yj1,…,yjN) and αj stand for the jth node and its associated weights, respectively. uN(yj,x) is again the deterministic solution of (2.8) with fixed yj. The choice of the nodes and weights should be made such that(3.12)UQ(f):=∑j=1Qf(yj)αj
is an approximation to the integral ∫Γf(y)ρ(y)dy for sufficiently smooth functions f(y).


### 3.3. Points Selection-Sparse Grids

The selection of nodes is the key ingredient in the above pseudospectral method. It is essential that the nodal set is a good cubature rule such that multiple integrals can be well approximated by a weighted discrete sum in the form of (3.12). It is straightforward in one-dimensional spaces N=1 and the optimal choice is usually the Gauss quadratures. The challenge is in multidimensional spaces, especially for large dimensions N≫1. We employ the sparse grids in our approach and refer to [14, 15] for more details along this line.Sparse grids were first proposed in [20]. It is a linear combination of product formulas, which is chosen in such a way that an interpolation property for N=1 is preserved for N>1 as much as possible. The sparse grids consist of much less number of nodes than that of the full tensor grids and were found to be particular useful in solving random differential equations by stochastic collocation approach in high-dimensional random space in [12, 21]. In [21], Nobile et al. established strong error estimates for the fully discrete solution using Lq norms and the computational efficiency demonstrates algebraic convergence with respect to the total number of collocation points.For every directioni=1,…,N, suppose we can construct a good one-dimensional interpolation operator 𝒰i:(3.13)Ui(f)=∑j=1mif(yji)⋅αji,
based on nodal sets Θi={y1i,…,ymii}⊂Γi. Following this formula, the Smolyak algorithm is given by(3.14)UQ(f):=A(J,N)=∑J+1≤|i|≤J+N(-1)J+N-|i|(N-1J+N-|i|)⋅(Ui1⊗⋅⋅⋅⊗UiN),
where |i|=(i1,…,iN)∈ℕN. To compute 𝒜(J,N), we only need to evaluate function on the sparse grid(3.15)H(J,N)=⋃J+1≤|i|≤J+N(Θi1×⋅⋅⋅×ΘiN).In this paper, both of theClenshaw-Curtis abscissas and Gaussian abscissas are used in the construction of the smolyak formula (3.13)–(3.15).Clenshaw-Curtis Abscissas
These abscissas are the extrema of Chebyshev polynomials, and for any choice ofmi>1,
(3.16)yji=-cos⁡(π(j-1)mi-1),j=1,…,mi.
In addition, one sets y1i=0 if mi=1 and lets the number of abscissas mi in each level to grow according to the following formula:
(3.17)m1=1,mi=2i-1+1fori>1.
With this particular choice, one obtains nested sets of abscissas, that is, Θi⊂Θi+1 and subsequently ℋ(J,N)⊂ℋ(J+1,N).Gaussian Abscissas
We also propose to use Gaussian abscissas, that is, the zeros of the orthogonal polynomials with respect to a weightρ. However, these Gaussian abscissas are in general not nested. Regardless, as in the Clenshaw-Curtis case, we choose the number mi of abscissas that are used by 𝒰i as in (3.17). The natural choice of the weight should be the probability density function ρ of the random variables Yi(ω) for all i.


### 3.4. Deterministic Plate Bending Solvers

Since at each collocation point, the deterministic plate bending problem (2.8) is a fourth-order partial differential equation in space variables, the corresponding conforming finite-element space should be C1-smooth at least which will result in an expensive computational cost [22, 23]. Here, we use the nonconforming Morley element to solve it. Compared to the conforming element methods, its computational cost alleviates greatly. Let 𝒯h be the shape of regular triangulation of 𝒟 such that each element K is an open triangle of size h. Then, we construct the usual Morley element space as follows (cf. [23]). For each K∈𝒯h with {pi}i=13 and {mi}i=13 as its three vertices and midpoints, respectively, the shape function space is chosen as P2(K) equipped with the nodal variables(3.18)ΣK={v(pi),∂nv(mi),1≤i≤3},
where Pk(K) denotes the space of all polynomials with the total order no more than k on K. We also use p or m to represent a vertex or a midpoint of a triangle in 𝒯h. Thus the Morley nonconforming element space related to H02(𝒟) is given by(3.19)VhM(D):={v∈L2(D);v|K∈P2(K)∀K∈Th,v|K(p)=v|K′(p),∂nv|K(m)=∂nv|K′(m)∀p,m∈∂K∩∂K′,∀K,K′∈Th,v(p)=∂nv(m)=0∀p,m∈∂DL2}.Then the nonconforming Morley element method at each collocation pointsyj for (2.8) is to find uNh(yj,x)∈VhM, such that for j=1,2,…,M,(3.20)∫DMIJ(uNh(yj,x))KIJ(vh)dx=∫DfN(yj,x)vhdx∀vh∈VhM.


### 3.5. Summary of Our Algorithm

To sum up, our pseudospectral algorithm for the stochastic plate bending problem consists of the following steps:(1)
choose a collocation sets{yj,αj} in space Γ according to the Clenshaw-Curtis abscissas or Gaussian abscissas described in Section 3.3;(2)
for eachj=1,…,M, solve problem (2.8) with a fixed parameter yj=(yj1,…,yjN) via the Morley element method (3.20);(3)
evaluate the approximate gPC expansion coefficients (3.14) and finally construct the N-variate, Pth-order gPC approximation (3.13).


## 3.1. Generalized Polynomial Chaos Expansion

In the finite-dimensional random spaceΓ, the gPC expansion seeks to approximate a random function via orthogonal polynomials of random variables. Let us define one-dimensional orthogonal polynomial spaces with respect to the measure ρi(Yi)dYi in Γi,(3.1)Wi,mi:={v:Γi→R:v∈span⁡{ϕm(Yi)}m=0mi},i=1,…,N,
where {ϕm(Yi)} are a set of orthogonal polynomials satisfying the orthogonality conditions(3.2)∫Γiρi(Yi)ϕm(Yi)ϕn(Yi)dYi=hm2δmn,
with normalization factors hm2=∫Γiρiϕm2dYi. With proper scaling, one can always normalize the bases such that hm2≡1 and this will be adopted throughout this paper. The probability density function ρi(Yi) in the above orthogonality relation (3.2) serves as a role of integration weight, which in turn defines the type of orthogonal polynomials {ϕm}. We refer this to [10, 14] for more details.The correspondingN-variate orthogonal polynomial space in Γ is defined by(3.3)WNP=⨂m∈ΛWi,mi,
where the tensor product is taken over all possible combinations of the multi-index m=(m1,…,mN) satisfying m∈Λ:={α:|α|=∑i=1Nmi≤P}. Thus, WNP is the space of N-variate orthogonal polynomials of total degree at most P, and the number of basis polynomials in (3.3) is dim(WNP)=(P+N)!/P!N!. Let {Φm(y)} be the N-variate orthonormal polynomials from WNP. They are constructed as products of a sequence of univariate polynomials in each direction of Yi, i=1,…,N, that is,(3.4)Φm(y)=∏m∈Λϕmi(Yi),
where mi is the order of the polynomial ϕ(Yi) in Yi for 1≤i≤d. For the N-variate polynomial Φm, each index 1≤m≤M=dim(WNP) corresponding to a unique sequence (m1,…,mN)∈Λ; they are also satisfying the orthogonality conditions(3.5)∫ΓΦm(y)Φn(y)ρ(y)dy=δmn,1≤m,n≤M.The finite-order, Pth-order, gPC approximation of the stochastic transverse displacementuN(y,x) in (2.8) is(3.6)PNPu:=uNP(y,x)=∑m=1Mûm(x)Φm(y),M=(N+PN),
where ℙNP denotes the orthogonal projection operator from Lρ2(Γ) onto WNP and {ûm} are the Fourier coefficients defined by(3.7)ûm(x)=∫ΓuN(y,x)Φm(y)ρ(y)dy=E[uN(y,x)Φm(y)].When a sufficient accurate gPC approximation (3.6) is available, one has in fact an analytical representation of uN(y,x) in terms of the random inputs y. Therefore, practically all statistical information can be retrieved in a straightforward manner. For example, the mean solution of the transverse displacement is(3.8)E(uN(y,x))≈E(uNP(y,x))=∫(∑m=1MûmΦm(y))ρ(y)dy=û1.
The variance of the solution can be obviously approximated as(3.9)Var⁡(uN(y,x))≈E[(uNP(y,x)-E(uNP(y,x)))2]=∑m=2M[ûm2].


## 3.2. Pseudospectral Approach

In the pseudospectral approach, we again seek an approximate solution in the form of gPC expansion, similar to (3.6), that is, for any x∈𝒟,(3.10)INPu:=wNP(y,x)=∑m=1Mŵm(x)Φm(y),
where 𝕀NP denotes the other projection operator from Lρ2(Γ) onto WNP and the expansion coefficients are determined as(3.11)ŵm=∑j=1QuN(yj,x)Φm(yj)αj,m=1,…,M.
Here {yj,αj}j=1Q are a set of nodes and weights, where yj=(yj1,…,yjN) and αj stand for the jth node and its associated weights, respectively. uN(yj,x) is again the deterministic solution of (2.8) with fixed yj. The choice of the nodes and weights should be made such that(3.12)UQ(f):=∑j=1Qf(yj)αj
is an approximation to the integral ∫Γf(y)ρ(y)dy for sufficiently smooth functions f(y).


## 3.3. Points Selection-Sparse Grids

The selection of nodes is the key ingredient in the above pseudospectral method. It is essential that the nodal set is a good cubature rule such that multiple integrals can be well approximated by a weighted discrete sum in the form of (3.12). It is straightforward in one-dimensional spaces N=1 and the optimal choice is usually the Gauss quadratures. The challenge is in multidimensional spaces, especially for large dimensions N≫1. We employ the sparse grids in our approach and refer to [14, 15] for more details along this line.Sparse grids were first proposed in [20]. It is a linear combination of product formulas, which is chosen in such a way that an interpolation property for N=1 is preserved for N>1 as much as possible. The sparse grids consist of much less number of nodes than that of the full tensor grids and were found to be particular useful in solving random differential equations by stochastic collocation approach in high-dimensional random space in [12, 21]. In [21], Nobile et al. established strong error estimates for the fully discrete solution using Lq norms and the computational efficiency demonstrates algebraic convergence with respect to the total number of collocation points.For every directioni=1,…,N, suppose we can construct a good one-dimensional interpolation operator 𝒰i:(3.13)Ui(f)=∑j=1mif(yji)⋅αji,
based on nodal sets Θi={y1i,…,ymii}⊂Γi. Following this formula, the Smolyak algorithm is given by(3.14)UQ(f):=A(J,N)=∑J+1≤|i|≤J+N(-1)J+N-|i|(N-1J+N-|i|)⋅(Ui1⊗⋅⋅⋅⊗UiN),
where |i|=(i1,…,iN)∈ℕN. To compute 𝒜(J,N), we only need to evaluate function on the sparse grid(3.15)H(J,N)=⋃J+1≤|i|≤J+N(Θi1×⋅⋅⋅×ΘiN).In this paper, both of theClenshaw-Curtis abscissas and Gaussian abscissas are used in the construction of the smolyak formula (3.13)–(3.15).Clenshaw-Curtis Abscissas
These abscissas are the extrema of Chebyshev polynomials, and for any choice ofmi>1,
(3.16)yji=-cos⁡(π(j-1)mi-1),j=1,…,mi.
In addition, one sets y1i=0 if mi=1 and lets the number of abscissas mi in each level to grow according to the following formula:
(3.17)m1=1,mi=2i-1+1fori>1.
With this particular choice, one obtains nested sets of abscissas, that is, Θi⊂Θi+1 and subsequently ℋ(J,N)⊂ℋ(J+1,N).Gaussian Abscissas
We also propose to use Gaussian abscissas, that is, the zeros of the orthogonal polynomials with respect to a weightρ. However, these Gaussian abscissas are in general not nested. Regardless, as in the Clenshaw-Curtis case, we choose the number mi of abscissas that are used by 𝒰i as in (3.17). The natural choice of the weight should be the probability density function ρ of the random variables Yi(ω) for all i.


## 3.4. Deterministic Plate Bending Solvers

Since at each collocation point, the deterministic plate bending problem (2.8) is a fourth-order partial differential equation in space variables, the corresponding conforming finite-element space should be C1-smooth at least which will result in an expensive computational cost [22, 23]. Here, we use the nonconforming Morley element to solve it. Compared to the conforming element methods, its computational cost alleviates greatly. Let 𝒯h be the shape of regular triangulation of 𝒟 such that each element K is an open triangle of size h. Then, we construct the usual Morley element space as follows (cf. [23]). For each K∈𝒯h with {pi}i=13 and {mi}i=13 as its three vertices and midpoints, respectively, the shape function space is chosen as P2(K) equipped with the nodal variables(3.18)ΣK={v(pi),∂nv(mi),1≤i≤3},
where Pk(K) denotes the space of all polynomials with the total order no more than k on K. We also use p or m to represent a vertex or a midpoint of a triangle in 𝒯h. Thus the Morley nonconforming element space related to H02(𝒟) is given by(3.19)VhM(D):={v∈L2(D);v|K∈P2(K)∀K∈Th,v|K(p)=v|K′(p),∂nv|K(m)=∂nv|K′(m)∀p,m∈∂K∩∂K′,∀K,K′∈Th,v(p)=∂nv(m)=0∀p,m∈∂DL2}.Then the nonconforming Morley element method at each collocation pointsyj for (2.8) is to find uNh(yj,x)∈VhM, such that for j=1,2,…,M,(3.20)∫DMIJ(uNh(yj,x))KIJ(vh)dx=∫DfN(yj,x)vhdx∀vh∈VhM.


## 3.5. Summary of Our Algorithm

To sum up, our pseudospectral algorithm for the stochastic plate bending problem consists of the following steps:(1)
choose a collocation sets{yj,αj} in space Γ according to the Clenshaw-Curtis abscissas or Gaussian abscissas described in Section 3.3;(2)
for eachj=1,…,M, solve problem (2.8) with a fixed parameter yj=(yj1,…,yjN) via the Morley element method (3.20);(3)
evaluate the approximate gPC expansion coefficients (3.14) and finally construct the N-variate, Pth-order gPC approximation (3.13).


## 4. Numerical Results

This section illustrates the computational performance of our pseudospectral method for the stochastic Kirchhoff plate bending problem (2.1). Uncertainty is first imposed on Young’s modulus of the plate and then on the Poisson ratio. In both cases, uncertainty is modeled by parameterized stochastic process. The plate domain is 𝒟=[0,1]2 and it is subjected to a deterministic load f(ω,x)=cos(x1)sin(x2). The plate thickness is h=0.02. The deterministic plate bending problem is solved by the Morley element method introduced in Section 3.4 and the space domain 𝒟 is partitioned into 648 triangles with 1369 unknowns.


### 4.1. Elastic Plate with Random Young’s Modulus

In this example, the log Young’s ModulusY(ω,x)=log(EN(ω,x)-100) is modeled as a parameterized stochastic process as in [13, page 2336]:(4.1)Y(ω,x)=1+Y1(ω)(πL2)1/2+∑i=2Nζifi(x)Yi(ω),
where(4.2)ζi:=(πL)1/2exp⁡(-(⌊i/2⌋πL)28),ifi>1,fi(x):={sin⁡⁡(⌊i/2⌋πx1Lp),ifieven,cos⁡⁡(⌊i/2⌋πx1Lp),ifiodd.
Here {Yi}i=1∞ are supposed to be independent and uniformly distributed in the interval [-1,1]. Thus a family of Legendre polynomials is used to construct space Φm. For x1∈[0,1], let Lc be a desired physical correlation length for Young’s Modulus E. Then the parameter Lp and L are Lp=max{1,2Lc} and L=Lc/Lp, respectively.We first study the convergence of the Smolyak sparse grid approximation and investigate the behavior when the levelJ in the Smolyak formula is increased linearly with a fixed P=4. The sparse grids are built upon either Clenshaw-Curtis abscissasorGaussian abscissas. To estimate the computational error in the Jth level we approximate ∥E[wNP,J(y,x)-wNP,J+1(y,x)]∥ as in [21], where wNP,J(y,x) is the solution obtained by the Jth level sparse grids. On the other hand, to investigate the performance of the Smolyak approximation with fixed gPC polynomial order by varying the correlation length we also include the cases where Lc=1/64 and Lc=1/2 for both N=4 and N=10. The computational results are shown in Figures 1 and 2. The results reveal that for Lc=1/64, the error decreases (sub)exponentially, as the level J increases. However, for Lc=1/2, the convergence rate slows down. We also observe that the convergence rate is dimension dependent and slightly deteriorates as N increases.Figure 1
The rate of convergence of the pseudospectral approach with different correlation lengthsLc=1/64 and Lc=1/2 with fixed P=4 for N=4.Figure 2
The rate of convergence of the pseudospectral approach with different correlation lengthsLc=1/64 and Lc=1/2 with fixed P=4 for N=10.We next compare the computational results, obtained by different orders of chaos polynomialsP∈{1,2,3}, with the standard Monte Carlo Finite Element method (MCM) to illustrate the accuracy of our method. Figure 3 shows the expected value and variance of the stochastic displacement at x2=1/2. The results show that the accuracy of the pseudospectral method increases with increasing polynomial order.Figure 3
Comparisons of the expectation and variance of the deflectionu at x2=1/2, derived from MCM and pseudospectral approach, for Lc=1/64, with fixed J=3.


### 4.2. Elastic Plate with Random Poisson Ratio

In this example, Young’s modulus is assumed to be deterministic and the poisson ratio is represented by the following parameterized stochastic process:(4.3)ν(ω,x)=0.3+0.05*(Y1(ω)cos⁡(πx1)+Y2(ω)sin⁡(πx1)+Y3(ω)cos⁡(πx2)+Y4(ω)cos⁡(πx2)),
where {Yi}i=14 are supposed to be independent and uniformly distributed in the interval [-1,1]. To evaluate the error of approximate solutions, relative error functions in expected value and variance are defined as following:(4.4)Errormean=|EwNP(x1,x2)-Eumc(x1,x2)||Eumc(x1,x2)|,Errorvar=|Var⁡wNP(x1,x2)-Var⁡umc(x1,x2)||Var⁡umc(x1,x2)|,
where EwNP and VarwNP are the mean and variance obtained by pseudospectral solution, Eumc and Varumc are obtained via the Monte Carlo simulation, based on 5000 samples.Figure4 shows the expected value and variance of the stochastic displacement obtained via different polynomial orders P∈{1,2,3} at the centerline x2=1/2 with fixed sparse gird level J=3. Figure 5 shows the relative error defined above at x2=1/2. The computational results show that the accuracy of our method with increasing polynomial order.Figure 4
Comparisons of the expectation and variance of the deflectionu at x2=1/2, derived from MCM and pseudospectral approach with fixed J=3.Figure 5
Relative error in expected value and variance of the deflectionu at x2=1/2, derived from MC and pseudospectral approach with fixed J=3.


## 4.1. Elastic Plate with Random Young’s Modulus

In this example, the log Young’s ModulusY(ω,x)=log(EN(ω,x)-100) is modeled as a parameterized stochastic process as in [13, page 2336]:(4.1)Y(ω,x)=1+Y1(ω)(πL2)1/2+∑i=2Nζifi(x)Yi(ω),
where(4.2)ζi:=(πL)1/2exp⁡(-(⌊i/2⌋πL)28),ifi>1,fi(x):={sin⁡⁡(⌊i/2⌋πx1Lp),ifieven,cos⁡⁡(⌊i/2⌋πx1Lp),ifiodd.
Here {Yi}i=1∞ are supposed to be independent and uniformly distributed in the interval [-1,1]. Thus a family of Legendre polynomials is used to construct space Φm. For x1∈[0,1], let Lc be a desired physical correlation length for Young’s Modulus E. Then the parameter Lp and L are Lp=max{1,2Lc} and L=Lc/Lp, respectively.We first study the convergence of the Smolyak sparse grid approximation and investigate the behavior when the levelJ in the Smolyak formula is increased linearly with a fixed P=4. The sparse grids are built upon either Clenshaw-Curtis abscissasorGaussian abscissas. To estimate the computational error in the Jth level we approximate ∥E[wNP,J(y,x)-wNP,J+1(y,x)]∥ as in [21], where wNP,J(y,x) is the solution obtained by the Jth level sparse grids. On the other hand, to investigate the performance of the Smolyak approximation with fixed gPC polynomial order by varying the correlation length we also include the cases where Lc=1/64 and Lc=1/2 for both N=4 and N=10. The computational results are shown in Figures 1 and 2. The results reveal that for Lc=1/64, the error decreases (sub)exponentially, as the level J increases. However, for Lc=1/2, the convergence rate slows down. We also observe that the convergence rate is dimension dependent and slightly deteriorates as N increases.Figure 1
The rate of convergence of the pseudospectral approach with different correlation lengthsLc=1/64 and Lc=1/2 with fixed P=4 for N=4.Figure 2
The rate of convergence of the pseudospectral approach with different correlation lengthsLc=1/64 and Lc=1/2 with fixed P=4 for N=10.We next compare the computational results, obtained by different orders of chaos polynomialsP∈{1,2,3}, with the standard Monte Carlo Finite Element method (MCM) to illustrate the accuracy of our method. Figure 3 shows the expected value and variance of the stochastic displacement at x2=1/2. The results show that the accuracy of the pseudospectral method increases with increasing polynomial order.Figure 3
Comparisons of the expectation and variance of the deflectionu at x2=1/2, derived from MCM and pseudospectral approach, for Lc=1/64, with fixed J=3.


## 4.2. Elastic Plate with Random Poisson Ratio

In this example, Young’s modulus is assumed to be deterministic and the poisson ratio is represented by the following parameterized stochastic process:(4.3)ν(ω,x)=0.3+0.05*(Y1(ω)cos⁡(πx1)+Y2(ω)sin⁡(πx1)+Y3(ω)cos⁡(πx2)+Y4(ω)cos⁡(πx2)),
where {Yi}i=14 are supposed to be independent and uniformly distributed in the interval [-1,1]. To evaluate the error of approximate solutions, relative error functions in expected value and variance are defined as following:(4.4)Errormean=|EwNP(x1,x2)-Eumc(x1,x2)||Eumc(x1,x2)|,Errorvar=|Var⁡wNP(x1,x2)-Var⁡umc(x1,x2)||Var⁡umc(x1,x2)|,
where EwNP and VarwNP are the mean and variance obtained by pseudospectral solution, Eumc and Varumc are obtained via the Monte Carlo simulation, based on 5000 samples.Figure4 shows the expected value and variance of the stochastic displacement obtained via different polynomial orders P∈{1,2,3} at the centerline x2=1/2 with fixed sparse gird level J=3. Figure 5 shows the relative error defined above at x2=1/2. The computational results show that the accuracy of our method with increasing polynomial order.Figure 4
Comparisons of the expectation and variance of the deflectionu at x2=1/2, derived from MCM and pseudospectral approach with fixed J=3.Figure 5
Relative error in expected value and variance of the deflectionu at x2=1/2, derived from MC and pseudospectral approach with fixed J=3.



---

*Source: 750605-2012-05-22.xml*