High-order Spectral Stochastic Finite Element Analysis of Stochastic Elliptical Partial Differential Equations

This study presents an experiment of improving the performance of spectral stochastic finite element method using high-order elements. This experiment is implemented through a two-dimensional spectral stochastic finite element formulation of an elliptic partial differential equation having stochastic coefficients. Deriving this spectral stochastic finite element formulation couples a two-dimensional deterministic finite element formulation of an elliptic partial differential equation with generalized polynomial chaos expansions of stochastic coefficients. Further inspection of the performance of resulting spectral stochastic finite element formulation with adopting linear and quadratic (9-node or 8-node) quadrilateral elements finds that more accurate standard deviations of unknowns are surprisingly predicted using quad-ratic quadrilateral elements, especially under high autocorrelation function values of stochastic coefficients. In addition, creating spectral stochastic finite element results using quadratic quadrilateral elements is not unacceptably time-consuming. Therefore, this study concludes that adopting high-order elements can be a lower-cost method to improve the performance of spectral stochastic finite element method.


Introduction
A stochastic partial differential equation is a partial differential equation having stochastic coefficients or forcing terms.Problems expressed as stochastic partial differential equations include such as population dynamics and elastostatics with the uncertainty in the spatial variability of mechanical properties.The spectral stochastic finite element method [1] may be the most popular numerical tools for solving stochastic partial differential equation.Briefly, deriving a spectral stochastic finite element formulation couples a deterministic finite element formulation and representations of those stochastic coefficients and forcing terms.These representations of stochastic forcing terms and coefficients are derived by such as generalized polynomial chaos [2] and Karhunen-Loève expansions [3].
Numerous spectral stochastic finite element formulations are available for some branches of science and engineering.References [4,5] are two recent examples.Nevertheless, the performance of spectral stochastic fi-nite element method is not always satisfactory.For example, the spectral stochastic finite element method predicts less accurate mean values or standard deviations of random fields than the spectral stochastic meshless local Petrov-Galerkin method does [6].Similar experiences bring about the motive for improving the performance of spectral stochastic finite element method.
Since a spectral stochastic finite element formulation contains a deterministic finite element solution and random field representations, improving the performance of spectral stochastic finite element method may be through adopting more accurate deterministic finite element formulation or random field representations.Available studies (e.g.[1]), which were devoted to evaluate the performance of spectral stochastic finite element method, seem to focus on the latter method.Experiences of applying a spectral stochastic finite element formulation using high-order elements seem to be unavailable.After a web search, only Ngah and Young (2007) [7] had adopted quadratic quadrilateral elements to predict the performance of composite structures.Thus, this study focuses on applying the spectral stochastic finite element method using high-order elements and evaluating the accuracy of corresponding spectral stochastic finite element results.A two-dimensional elliptical partial differential equation having stochastic coefficients is chosen as the model problem.A spectral stochastic finite element formulation of such a stochastic partial differential equation is derived by coupling a two-dimensional deterministic finite element formulation of an elliptical partial differential equation with generalized polynomial chaos expansions of stochastic coefficients.Two benchmark problems having deterministic analytical solutions are then introduced to test the performance of resulting spectral stochastic finite element formulation with adopting linear and high-order finite elements.
The remainder of this study is organized in three sections.Section 2 presents the two-dimensional spectral stochastic finite element formulation of an elliptical partial differential equation having stochastic coefficients.Section 3 evaluates the performance of resulting spectral stochastic finite element formulation.Based on this performance evaluation, Section 4 draws some conclusion.

Spectral Stochastic Finite Element Formulation
Suppose x and y are the spatial coordinates,  is an event in the probability space,  is a problem domain, and  is its boundary.
where u is the unknown (such as temperature and stream function), on in which 0 and n are two known functions, U U Q  is the essential boundary,  T is the natural boundary, q n is the secondary variable, n x and n y are the components of a unit vector n normal to  T , and .
To derive a spectral stochastic finite element formulation of (1), a two-dimensional deterministic finite element formulation of an elliptical partial differential equation is needed.Suppose w is the test function.The weak form of (1) is where the superscript e denotes an element.Simplifying (3) by the divergence theorem yields where  e is the boundary of  e , and ds is the arc length of an infinitesimal line along the boundary  e .Next, approximate u over each  e by Equation ( 6) can be summarized in the form as Moreover, (7) in which N T is the total number of nodes in the problem domain .Since and f are dependent upon x, y, and , (10) is not ready for use.The generalized polynomial chaos expansion is chosen to estimate the distribution of and f.Similar manipulating the published study [6], generalized polynomial chaos expansions of and f are defined by n P n P   , P is the highest order of Ψ, and n is the total number of uncorrelated random variables.
For facilitating the computation of (11), Ψ 0 , , and 0 f are; respectively, set to 1, mean values of a ij , a 00 and f; respectively.Furthermore, computing coefficients in which   is computed as follows: If f and g are two functions,   is computed by 1) Continuous case: 2) Discrete case: , where  are the weighting functions.In (1), the succeeding study focuses on the continuous random fields; therefore, Table 1 [6] lists examples of orthogonal polynomials, statistical distributions and weighting functions to generate Table 1.Examples of polynomials and corresponding weighting functions and statistical distributions for generating the generalized polynomial chaos [6].
Copyright © 2013 SciRes.AM Similarly, substituting (11) into (8b) yields Meanwhile, the generalized polynomial chaos expansion of u is Modifying ( 10) with ( 14) and ( 16) results in Requiring the residual resulting from a finite representation of u (i.e.truncating 1 ) to be orthogonal to the approximation space spanned by , Solving ( 18) can obtain . Accumulating the resulting J û values can construct the generalized polynomial chaos expansion of u.

Results and Discussions
Two benchmark problems are introduced to evaluate the performance of (18) with adopting linear and high-order elements.The first benchmark problem involves a heat conduction problem over a square region.The second problem involves the transverse deflection of a square membrane.

Heat Conduction over a Square Region
As outlined by Figure 1, suppose a heat conduction problem over a unit square region in where  is the gradient vector.The origin of x and y coordinates locates at the lower left corner.The boundaries and 0 x  0 y  are insulated.The other boundaries 1 x  and 1 y  are maintained at zero temperature.In addition, the square region is subjected to a uniform heat generation.
To predict the temperature T, the governing equation is where k is the thermal conductivity and 0 is the rate of uniform heat generation.In addition, the boundary conditions are Similarly manipulating (7), the deterministic finite element formulation of (19) is where . However, the succeeding study considers that the thermal conductivity k is stochastic.Assume the distribution of k is described by where k  is the mean value of k and this k  value is independent upon x and y.Meanwhile,


is a zero-mean, scalar, homogeneous random field with its autocorrelation function   equal to where d 1 and d 2 denote the correlation lengths, S k is the standard deviation of thermal conductivity and   , x y represent two points of the square region.
To compare the predicted temperature T with adopting linear or high-order elements, essential data is provided below 1) Define the problem domain  as and 0 x  1 0 1 y   .
Copyright © 2013 SciRes.AM 2) Generate two cases of the finite element discretizations.As outlined by Figure 2, the left sub-figure denotes the first case which is composed of 25 nodes and 16 linear quadrilateral elements.The right sub-figure denotes the second case which is composed of 81 nodes and 16 quadratic (9-node) quadrilateral elements.
3) Consulting with Table 1 and the previous study [6], experiment to represent k by the Lauguerre polynomial chaos.If the accuracy of corresponding spectral stochastic finite element results is unsatisfactory, apply another type of the generalized polynomial chaos.10 , 16, 1, N q is the total number of quadrature points in a finite element.Furthermore, in an attempt of quantifying the errors between Monte Carlo simulation and spectral stochastic finite element results, two error estimators   and  S are defined below in which the subscripts M and S denote the Monte Carlo simulation and spectral stochastic finite element method; respectively.In (28), choosing is better checked by observing Monte Carlo simulation results with respect to different sample value.Figures 3(a 2, are set according to (28).Similar practices are followed in the following.
Since we can easily expect that adopting quadratic (9-node) elements can obtain more accurate deterministic finite element results, Figure 4

Transverse Deflection of a Square Membrane
Suppose a membrane occupies a region and and its edges are fixed.Initially, the membrane is stretched so that the tension a in the membrane is uniform and that tension a is so large that it is not appreciably altered when the membrane is deflected by a distributed normal force To predict the transverse deflection u of the membrane, the governing equation is Due to symmetry, only one quadrant of the membrane is analyzed.The boundary conditions are Meanwhile, the deterministic finite element formulation of (30) is similar to (7) except that components e ij K and e i F are; respectively Nevertheless, the succeeding study assumes the tension a varies according to the following uniform distribution: where  a and S a are; respectively, the mean value and standard deviation of tension a, To compare the predicted deflection u with accounting for a stochastic tension a, essential data is provided below 1) Still define  as 0 x 1   and 0 1 y   .
2) Generate two cases of the finite element discretizations.As outlined by Figure 8, the left sub-figure denotes the first case which is composed of 81 nodes and 64 linear quadrilateral elements.The right sub-figure denotes the second case which is composed of.65 nodes and 16 quadratic (8-node) quadrilateral elements.
3) Consulting with Table 1, represent a by the Legendre polynomial chaos.
4) Generate Monte Carlo simulation results by first sampling the tension a according to (31).Each sample of the tension a is then substituted into the aforementioned analytical solution of (29) to predict a sample of the deflection u.Similarly manipulating (26), the mean value  u and standard deviation S u of deflection u are equal to

Closure
As introduced in Section 1, the most popular numerical tool for solving stochastic partial differential equations may be the spectral stochastic finite element method.Numerous resources are available for generating spectral stochastic finite element results.As compared to the Monte Carlo simulation, applying the spectral stochastic finite element method doesn't require sampling the existing random fields sufficiently; thus, creating spectral stochastic finite element results is usually time-saving.


Probably since deterministic analytical solutions are usually unavailable for producing Monte Carlo simulation results, the accuracy of spectral stochastic is not often discussed and linear elements were usually adopted in applying this stochastic numerical method.However, the succeeding study demonstrates that adopting highorder (e.g.quadratic) elements can improve the performance of spectral stochastic finite element method.The previous section has shown that adopting linear elements has the danger of obtaining uncontrollable errors between Monte Carlo simulation and spectral stochastic finite element results.Whereas, adopting quadratic (9-node or 8-node) elements to apply the spectral stochastic finite element method stably produces more accurate predicted On the other hand, comparing Figures 10(a  mean values and standard deviations under high autocorrelation function values of existing stochastic coefficients ranges.In addition, the time spent to apply the spectral stochastic finite element method with using quadratic quadrilateral elements is not unacceptably time-consuming.
In conclusion, replacing linear elements with highorder elements to apply the spectral stochastic finite element method can be as a low-cost method to improve the performance of this stochastic numerical method.


denote multi-dimensional uncorrelated random variables having zero mean and unit variance (for facilitating the computation of mean values and standard deviations of 00 and f), N PC is equal to ,

C
in which   is the ensemble average.For example,

Figure 1 .
Figure 1.A heat conduction problem over a square region (not to scale).

4 )
Generate Monte Carlo simulation results to serve as the accuracy standard in evaluating the accuracy of spectral stochastic finite element results.A Monte Carlo simulation is implemented by first sampling the thermal conductivity k according to (24).Each sample of the thermal conductivity k is then substituted into (23) to predict a sample of the temperature T. If a sufficient amount of samples of k are created, the corresponding mean value  T and standard deviation S T of temperature T will approach their exact values.These Monte Carlo simulation-based  T and S T are computed by (e.g.[1the total number of samples for implementing the Monte Carlo simulation, and the subscript j denotes T predicted using the j-th sample of thermal conductivity k.Meanwhile, similarly manipulating (16), the generalized polynomial chaos expansion of tempera-N ture T is and the corresponding spectral stoelement-based predicted  T and S T are computed by (e.g.[1])  

Figure 2 .
Figure 2. Finite element discretizations for analyzing the first benchmark problem (not to scale).

Figure 3 .Figure 2 .
Figure 3. Variation of Monte Carlo simulation results with respect to different values: (a) mean value; (b) standard deviation (First benchmark problem).N sample

Figure 4 .
Figure 4. Comparison of accuracy of spectral stochastic finite element results with respect to linear and quadratic quadrilateral elements: (a) mean value; (b) standard deviation (First benchmark problem,  = mean value, S = standard deviation).

Figure 2 .values with respect to Figure 2 .
ting quadratic quadrilateral elements consequently predict more accurate mean values  T , since the correis surprising.Adopting quadratic (9-node) elements can also obtain more accurate predicted standard deviation S T values, since the corresponding smaller.Meanwhile, Table2indicates that the time spent to generate spectral stochastic finite element results is not unacceptably time-consuming.However, we may argue that Figures4(a) and (b) only outline the effects of spacings of any two connecting nodes on the accuracy of spectral stochastic finite elements.As the nodal distribution becomes denser, obtaining more accurate spectral stochastic finite element results may be expected.To figure out this argument, the problem domain  is re-discretized by using 64   8 8  linear quadrilateral elements and 81 equally-spaced nodes.Based on the resulting finite element discretization and Figure 2(b), Table3inspects variation of to these two cases of finite element discretization.Inspecting Table3finds that adopting quadratic (9node) elements still produces smaller is, the improvement of accuracy of spectral stochastic finite element results in Figure 4(b) is not due to the denser nodal distribution in the right sub-figure of Figure 2. In fact, comparing Table 3 and Figures 4(a) and (b) finds that adopting a denser nodal distribution but still using linear elements only slightly improve the accuracy of spectral stochastic finite element results.Furthermore, increase the k k S  value from 0.12 to 0.32.Figures 5(a) and (b) present the corresponding Next, change the d 2 value from 1.0 to 2.0 (but revert the k k S  value to 0.12).Figures 6(a) and (b) present the corresponding The incentive of plotting Figures 5(a) and (b) comes from the published study [6] that high autocorrelation function values of stochastic coefficients (i.e.  values) has apparent effects on the accuracy of spectral stochastic finite element results.Further inspection of Figures 5(a) and (b) finds that adopting linear quadrilateral elements to apply the spectral stochastic finite element method

Figure 5 .Figure 6 .
Figure 5.Comparison of accuracy of spectral stochastic finite element results with respect to . k k S 0 32 

Figure 7
Figure 7 further illustrates the layout of problem domain  and boundary conditions.If the tension a is deterministic, the analytical solution of u is

Figure 7 .
Figure 7. Transverse deflection of a square membrane (not to scale).

5 )Figure 8 .
Figure 8. Finite element discretizations for analyzing the second benchmark problem (not to scale).

Figures 9 (
Figures 9(a) and (b) report that choosing is still reasonable.As the sample value is larger than 10 6 , the Monte Carlo simulation-based predicted 6 sample 10 N  N  u x y   and remain approximately fixed.Meanwhile, we still find from Figures 10(a) and (b) that adopting quadratic (8-node) elements is conducive to improving the accuracy of spectral stochastic finite elmentbased predicted , since the corresponding S value is smaller.Note that the total number of nodes in the left sub-figure of Figure 8 is more than the total number of nodes in the right sub-figure of the same figure.Since similar results are found in Table3and Figures4(a) and (b), adopting more nodes but still using the same element type consequently improve the accuracy of corresponding spectral stochastic finite element results slightly.

Figure 9 .
Figure 9. Variation of Monte Carlo simulation results with respect to different N sample values: (a) mean value; (b) standard deviation (Second benchmark problem).

Figure 10 .Figure 11 .
Figure 10.Comparison of accuracy of spectral stochastic finite element results with respect to linear and quadratic quadrilateral elements: (a) mean value; (b) standard deviation (Second benchmark problem,  = mean value, S = standard deviation).

Table 2 . Comparison of the time spent to generate spectral stochastic finite element results adopting linear and quad- ratic (9 nodes) quadrilaterals.
*On a MacBook Pro with an Intel Core i5 Processor.