Multiscale stochastic finite element method on random field modeling of geotechnical problems – a fast computing procedure

Abstract

The Green-function-based multiscale stochastic finite element method (MSFEM) has been formulated based on the stochastic variational principle. In this study a fast computing procedure based on the MSFEM is developed to solve random field geotechnical problems with a typical coefficient of variance less than 1. A unique fast computing advantage of the procedure enables computation performed only on those locations of interest, therefore saving a lot of computation. The numerical example on soil settlement shows that the procedure achieves significant computing efficiency compared with Monte Carlo method.

Keywords

multiscale / finite element / settlement / perturbation / random field / geotechnical

Introduction

To model uncertainty of spatial and/or temporal variations widely present in natural and engineering systems, a variety of stochastic finite element methods (SFEM) have been formulated based on the standard displacement-based finite elements. In Xu [ 1] by distinguishing a quasi-weak form from a weak form in both real and random space, a unifying framework of variational formulation is presented covering both the displacement-based SFEM and the recently proposed Green-function-based Multiscale Stochastic Finite Element Method (MSFEM) [ 2, 3]. The study [ 1] shows that Monte Carlo, perturbation, and weighted integral SFEMs correspond to the quasi-weak form, while the weak form results in spectral SFEM, pseudo-spectral SFEM, and GFB-SFEM. Within the unifying framework, dynamic problems are further addressed especially to demonstrate the unique feature of the GFB-SFEM on problems with inputs characterized as random fields or random processes [ 4].

In Xu [ 3] the MSFEM is formulated based on the stochastic variational principle. In this work a perturbation-type MSFEM is proposed on random field geotechnical problems typically characterized with coefficient of variance (c.o.v) less than unity. A unique fast computing advantage of the perturbation-type MSFEM enables computation performed only on those locations of interest, therefore saving a lot of computation resource spent on those not needed. For instance, in evaluation of a foundation settlement a major engineering interest is focused on the maximum settlement at the bottom of the foundation. A conventional probabilistic method such as Monte Carlo method obtains the probabilistic information not only at that particular location but all over the whole soil domain underneath. From this perspective, the new fast computing procedure of the MSFEM also emphasizes this type of computational saving as an important direction in developing efficient computational methods in engineering applications.

The perturbation-type MSFEM is first formulated in Section 2. In Section 3 the convergence of Neumann expansion is analyzed. The numerical example is provided in Section 4, with conclusion given at the end.

Formulation of perturbation-type MSFEM

A stochastic boundary value problem (BVP) of a geotechnical body in domain D is described via the following equations

(1) \{\begin{array}{ll}\nabla \cdot [L(x,\omega ){\nabla }^{s}u(x,\omega )]+f(x)=\mathrm{0}& inD\\ u_{\mathrm{0}}=\widehat{u}& on\partial D_u\\ (L_{\mathrm{0}}\nabla u_{\mathrm{0}})\cdot n=\widehat{t}& on\partial D_t\end{array},

where the displacement random field u(x,\omega ) results from the elastic moduli random field L(x,\omega ) , with ω denoting a sample in random space and the superscript s a symmetric operation, i.e., {\nabla }_j^s{u}_i\equiv e_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right) . Based on the multiscale method formulated in Xu [ 5], the stochastic BVP Eq. (1) can be decomposed into a slow-scale and a fast-scale BVP. The slow-scale BVP corresponds to the traditional deterministic model adopted in geotechnical design whereas soils are assumed spatially homogeneous, with the governing equations given as follows

\{\begin{array}{lll}\nabla \cdot \left[L_0{\nabla }^s{u}_0\left(x\right)\right]+f\left(x\right)=0\hfill & inD,\hfill & \left(2\mathrm{a}\right)\hfill \\ u_0=\widehat{u}\hfill & on\partial D_u,\hfill & \left(2\mathrm{b}\right)\hfill \\ \left(L_0\nabla u_0\right)\cdot n=\widehat{t}\hfill & on\partial D_t,\hfill & \left(2\mathrm{c}\right)\hfill \end{array}

which can be solved using the standard finite element method with the matrix form given as

(3) \mathbf{K}\mathbf{U}=\mathbf{F}.

Note the displacement boundary condition Eq. (2b) is accounted for in the stiffness matrix K of Eq. (3) by using the penalty scheme.

The fast-scale BVP accounts for spatial variability of soil properties, such as Young’s modulus, and the equations are given as

(4) \{\begin{array}{ll}\nabla \cdot [L_{\mathrm{0}}{\nabla }^s{u}^{\ast }(x)]+\nabla \cdot p=\mathrm{0}\hfill & inD\hfill \\ u^{\ast }=\mathrm{0}\hfill & on\partial D_u\hfill \\ (L_{\mathrm{0}}\nabla u^{\ast })\cdot n=\mathrm{0}\hfill & on\partial D_t\hfill \\ p\cdot n=\mathrm{0}\hfill & on\partial D_t\hfill \end{array},

where the so-called stress polarization p is defined as p=(L(x,\omega )-L_{\mathrm{0}})({\nabla }^s{u}_{\mathrm{0}}+{\nabla }^s{u}^{\ast }) .

In Xu [ 3, 5], the stress polarization is treated as a primary unknown to be solved by using an even-order stochastic variational principle and an odd-order principle, respectively. In this study to achieve fast computing, we rather go to the Green function solution of Eq. (4) directly by using divergence theorem that is expressed as [ 5]

u^{\ast }(x)={\displaystyle \underset{D}{\int }G(x,x\prime )(\nabla \cdot p(x\prime ))}dx\prime =-{\displaystyle \underset{D}{\int }(x,x\prime )p(x\prime )}dx\prime =-{\displaystyle \underset{D}{\int }(x,x\prime )\mathrm{\Delta }L(x\prime )(e_{\mathrm{0}}(x\prime )+e^{\ast }(x\prime ))}dx\prime ,

which further yields

(5) e^{\ast }(x)=-{\displaystyle \underset{D}{\int }\mathrm{\Gamma }(x,x\prime )\mathrm{\Delta }L(x\prime )[e_{\mathrm{0}}(x\prime )+e^{\ast }(x\prime )]}dx\prime ,

where \mathrm{\Delta }L=L(x,\omega )-L_{\mathrm{0}} , and the modified Green functions

(6) {}_{ikl}(x,x\prime )={\displaystyle \frac{\partial G_{ik}(x,x\prime )}{\partial x_l^{\prime }}},

(7) {\mathrm{\Gamma }}_{ijkl}(x,x\prime )={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{{\partial }^{\mathrm{2}}{G}_{ik}(x,x\prime )}{\partial x_j\partial x_l^{\prime }}}+{\displaystyle \frac{\partial G_{jk}(x,x\prime )}{\partial x_i\partial x_l^{\prime }}}\right].

Note in the above decomposition, u₀, e₀ and u*, e* represent the slow scale deformation and fast scale (i.e., perturbation) deformation, respectively, with the total displacement u = u₀ + u*, and the total strain e = e₀ + e*.

The solution of the fast-scale BVP Eq. (4) for the perturbation displacement therefore can be concisely expressed in Neumann expansion as

(8) u^{\ast }=-\mathrm{\Delta }L{e}_{\mathrm{0}}+\mathrm{\Delta }L(\mathrm{\Gamma }\mathrm{\Delta }L{e}_{\mathrm{0}})-\cdots \cdots

The Green functions are computed using finite element method, with the formula given as Ref. [5,6 Shen and Xu, 2010]

G_{ik}(x,x\prime )={\displaystyle \sum _p^N{\displaystyle \sum _q^N{\psi }_p(x)G_{pq}^{ik}{\psi }_q(x\prime )},}

{}_{ikl}(x,x\prime )={\displaystyle \sum _p^N{\displaystyle \sum _q^N{\psi }_p(x)G_{pq}^{ik}{\nabla }_l{\psi }_q(x\prime )},}

{\mathrm{\Gamma }}_{ijkl}(x,x\prime )={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \sum _p^N{\displaystyle \sum _q^N({\nabla }_j{\psi }_p(x)G_{pq}^{ik}{\nabla }_l{\psi }_q(x\prime )+{\nabla }_i{\psi }_p(x)G_{pq}^{jk}{\nabla }_l{\psi }_q(x\prime ))}},

where \psi denotes the finite element shape function, and N the total number of nodes. The numerical Green function \mathbf{G}={\mathbf{K}}^{-\mathrm{1}} is the inverse of the standard finite element stiffness matrix of the slow-scale BVP Eq. (3), with the subscript pq indicating two nodes and the superscript ik the axial directions of the nodal values on p and q. Note the inverse matrix is not necessarily computed explicitly, and for a large system by using the iteration method only the matrix K is required.

The Neumann expansion similar to Eq. (8) was given early in Kröner [ 7] in homogenization of heterogeneous materials. By employing the rigorous proof presented in Xu [ 2] for decomposition of the stochastic BVP Eq. (1) with general boundary conditions, hereby Eq. (8) is established for any BVP characterized with general boundary conditions. In next sections the concise notation used in Eq. (8) are followed.

Convergence of the perturbation method

In deterministic computation using a particular periodic sample L(x) , the convergence of Eq. (8) has been shown in Ref. [ 8- 10]. For stochastic computation, a convergence study was numerically given in Xu and Brady [ 11] using Fourier-Galerkin method. For stochastic finite element computation, the probabilistic information of Eq. (8) can be evaluated via probabilistic moments, e.g., the mean and the mean square values are given as

(9a) \overline{u^{\ast }}=-\overline{\mathrm{\Delta }L}{e}_{\mathrm{0}}+\overline{\mathrm{\Delta }L\mathrm{\Gamma }\mathrm{\Delta }L}{e}_{\mathrm{0}}-\cdots

(9b) \overline{{(u^{\ast })}^{\mathrm{2}}}=\overline{(\mathrm{\Delta }L{e}_{\mathrm{0}})(\mathrm{\Delta }L{e}_{\mathrm{0}})}-\mathrm{2}\overline{(\mathrm{\Delta }L{e}_{\mathrm{0}})(\mathrm{\Delta }L\mathrm{\Gamma }\mathrm{\Delta }L{e}_{\mathrm{0}})}+\cdots

where the overbar denotes the ensemble average.

For a geotechnical field, we can experimentally characterize the Lame’s constant l and shear modulus m as two correlated random fields, i.e.,

(10) \lambda (x,\omega )=\overline{\lambda }(x)+{\sigma }_{\lambda }Z(x,\omega ),

(11) \mu (x,\omega )=\overline{\mu }(x)+{\sigma }_{\mu }Z(x,\omega ),

where Z is the underlying zero mean unit variance random field. The elastic moduli is therefore expressed in terms of Z as

(12) L_{ijkl}(x,\omega )=\left[\overline{\lambda }{\delta }_{ij}{\delta }_{kl}+\overline{\mu }({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})\right]+[{\sigma }_{\lambda }{\delta }_{ij}{\delta }_{kl}+{\sigma }_{\mu }({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})]Z(x,\omega ),

or in a matrix form

(13) L(x,\omega )=\overline{L}+Z(x,\omega )L_{\sigma },

with

(14) \overline{L}=\left[\begin{array}{llllll}\mathrm{2}\overline{\mu }+\overline{\lambda }& \overline{\lambda }& \overline{\lambda }& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \overline{\lambda }& \mathrm{2}\overline{\mu }+\overline{\lambda }& \overline{\lambda }& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \overline{\lambda }& \overline{\lambda }& \mathrm{2}\overline{\mu }+\overline{\lambda }& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ 0& 0& 0& \overline{\mu }& 0& 0\\ 0& 0& 0& 0& \overline{\mu }& \mathrm{0}\\ 0& 0& 0& \mathrm{0}& \mathrm{0}& \overline{\mu }\end{array}\right],

(15) L_{\sigma }=\left[\begin{array}{llllll}\mathrm{2}{\sigma }_{\mu }+{\sigma }_{\lambda }& {\sigma }_{\lambda }& {\sigma }_{\lambda }& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\sigma }_{\lambda }& \mathrm{2}{\sigma }_{\mu }+{\sigma }_{\lambda }& {\sigma }_{\lambda }& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\sigma }_{\lambda }& {\sigma }_{\lambda }& \mathrm{2}{\sigma }_{\mu }+{\sigma }_{\lambda }& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ 0& 0& 0& {\sigma }_{\mu }& \mathrm{0}& \mathrm{0}\\ 0& 0& 0& 0& {\sigma }_{\mu }& \mathrm{0}\\ 0& 0& 0& 0& 0& {\sigma }_{\mu }\end{array}\right].

Note in the above expression, the concept of engineering shear strain is adopted in the strain tensor e₀. The following term is therefore obtained by using the auto-correlation function of the elastic moduli, as

(16a) \overline{\mathrm{\Delta }L(x)\mathrm{\Gamma }(x,x\prime )\mathrm{\Delta }L(x\prime )}=L_{\sigma }\mathrm{\Gamma }(x,x\prime )L_{\sigma }{\rho }_{\mathrm{2}}(x-x\prime ),

(16b) \overline{\mathrm{\Delta }L(x)(x,x\prime )\mathrm{\Delta }L(x\prime )}=L_{\sigma }(x,x\prime )L_{\sigma }{\rho }_{\mathrm{2}}(x-x\prime ),

where r₂ denotes coefficient of covariance of the random field Z.

Choose the reference moduli of the slow scale bvp Eq. (2) as L_{\mathrm{0}}=\overline{L} . The first term on the right hand side of Eq. (9a) vanishes, and the second term becomes

(17) \overline{\mathrm{\Delta }L\mathrm{\Gamma }\mathrm{\Delta }L}{e}_{\mathrm{0}}=L_{\sigma }\mathrm{\Gamma }{L}_{\sigma }{\rho }_{\mathrm{2}}{e}_{\mathrm{0}}.

The n-th term is similarly given as

(18) \overline{\mathrm{\Delta }L{(\mathrm{\Gamma }\mathrm{\Delta }L)}^{n-\mathrm{1}}}{e}_{\mathrm{0}}=L_{\sigma }{(\mathrm{\Gamma }{L}_{\sigma })}^{n-\mathrm{1}}{\rho }_n{e}_{\mathrm{0}},

where {\rho }_n denotes the n-th order correlation function of the random field Z, i.e., {\rho }_n(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_n)=\overline{Z(x_{\mathrm{1}})Z(x_{\mathrm{2}})\cdots Z(x_n)} . The right hand side of Eq. (18) can be bounded from above by replacing {\rho }_n with a constant C(n) that is a function of n, as follows:

(19) L_{\sigma }{(\mathrm{\Gamma }{L}_{\sigma })}^{n-\mathrm{1}}{\rho }_n{e}_{\mathrm{0}}\le C(n)\cdot L_{\sigma }{(\mathrm{\Gamma }{L}_{\sigma })}^{n-\mathrm{1}}{e}_{\mathrm{0}}.

To evaluate the right hand side of inequality Eq. (19), we need the following theorem that is applicable for any BVP subject to arbitrary boundary conditions.

THEOREM The operation of the modified Green function \mathrm{\Gamma } on the self-equilibrated stress L_{\mathrm{0}}^{}{e}_{\mathrm{0}} of the slow-scale BVP Eq. (2) results in the slow-scale strain solution e₀ itself, i.e., \mathrm{\Gamma }{L}_{\mathrm{0}}{e}_{\mathrm{0}}=e_{\mathrm{0}} .

Proof:

Construct a bvp as shown below

(20) \begin{array}{ll}\nabla \cdot [L_{\mathrm{0}}{\nabla }^s(u\prime (x)+u_{\mathrm{0}}(x))]+f(x)=\mathrm{0}\hfill & inD\hfill \\ u\prime =\mathrm{0};u_{\mathrm{0}}=\widehat{u}\hfill & on\partial D_u\hfill \\ (L_{\mathrm{0}}{\nabla }^{s}u\prime )\cdot n=\mathrm{0};(L_{\mathrm{0}}{\nabla }^s{u}_{\mathrm{0}})\cdot n=\widehat{t}\hfill & on\partial D_t\hfill \end{array}.

Compared with the slow scale BVP Eq. (2), the BVP Eq. (20) clearly has the solution u\prime (x)+u_{\mathrm{0}}(x) identical to u_{\mathrm{0}}(x) , and therefore we simply have u\prime (x)=\mathrm{0} . By using the Green function solution, u\prime (x) can be alternatively expressed from Eq. (20) as

(21) u\prime =G(\nabla \cdot (L_{\mathrm{0}}{e}_{\mathrm{0}})+f),

which by using divergence theorem becomes

(22) u\prime =-L_{\mathrm{0}}{e}_{\mathrm{0}}+G\left(\widehat{t}+f\right),

and

(23) e\prime =-\mathrm{\Gamma }{L}_{\mathrm{0}}{e}_{\mathrm{0}}+\prime \left(\widehat{t}+f\right),

with

(24) {\prime }_{ijk}(x,x\prime )={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{\partial G_{ik}(x,x\prime )}{\partial x_j}}+{\displaystyle \frac{\partial G_{jk}(x,x\prime )}{\partial x_i}}\right].

Since u\prime =e\prime =\mathrm{0} , and

(25) e_{\mathrm{0}}=\prime \left(\widehat{t}+f\right),

Equation (23) yields

(26) \mathrm{\Gamma }{L}_{\mathrm{0}}{e}_{\mathrm{0}}=e_{\mathrm{0}}.

By using the result Eq. (26), the right hand side of Eq. (19) becomes

(27) C(n)\cdot L_{\sigma }{(\mathrm{\Gamma }{L}_{\sigma })}^{n-\mathrm{1}}{e}_{\mathrm{0}}\le C(n){\gamma }^{n-\mathrm{1}}\cdot L_{\sigma }{e}_{\mathrm{0}},

where the coefficient of variance

(28) \gamma =\max \left\{{\displaystyle \frac{{\sigma }_{\mu }}{\overline{\mu }}},{\displaystyle \frac{{\sigma }_{\lambda }}{\overline{\lambda }}}\right\}.

Since Eq. (22) shows that

(29) L_{\mathrm{0}}{e}_{\mathrm{0}}=u_{\mathrm{0}},

we have

(30) C(n){\gamma }^{n-\mathrm{1}}\cdot L_{\sigma }{e}_{\mathrm{0}}\le C(n){\gamma }^n\cdot u_{\mathrm{0}},

and the inequality Eq. (19) becomes

(31) L_{\sigma }{(\mathrm{\Gamma }{L}_{\sigma })}^{n-\mathrm{1}}{\rho }_n{e}_{\mathrm{0}}\le C(n){\gamma }^n\cdot u_{\mathrm{0}}.

In general the increase of C with n is slower than {\gamma }^n when \gamma is smaller than 1. In geotechnical fields, the coefficient of variance (c.o.v) is normally smaller than 1, and we expect (31) and therefore the Neumann expansion Eq. (9) to converge. In Section 4 we show that for a soil with the c.o.v. of the Young’s modulus being 40%, even the first nonzero term of Eq. (9) is sufficient to capture the mean settlement accurately.

With respect to the variance

(32) \overline{u^{\mathrm{2}}}-{\overline{u}}^{\mathrm{2}}=\overline{{(u^{\ast })}^{\mathrm{2}}}-{\left(\overline{u^{\ast }}\right)}^{\mathrm{2}}=\overline{(\mathrm{\Delta }L{e}_{\mathrm{0}})(\mathrm{\Delta }L{e}_{\mathrm{0}})}-\mathrm{2}\overline{(\mathrm{\Delta }L{e}_{\mathrm{0}})(\mathrm{\Delta }L\mathrm{\Gamma }\mathrm{\Delta }L{e}_{\mathrm{0}})}+\cdots

the terms on the right hand side can be similarly given in terms of correlation functions, e.g., the first term is expressed as

(33) \overline{(\mathrm{\Delta }L{e}_{\mathrm{0}})(\mathrm{\Delta }L{e}_{\mathrm{0}})}=(L_{\sigma }{e}_{\mathrm{0}})({\rho }_{\mathrm{2}}{L}_{\sigma }{e}_{\mathrm{0}}).

As a closing remark of this section, below we discuss the physical ground for the assumption of an isotropic elastic random soil medium that certainly is convenient for illustration of the GFB multiscale method. In Ostoja-Starzewski [ 12, 13] it has been rigorously shown that a random field of the isotropic smooth elastic continuum is questionable given the statistically homogeneous microstructure. In geotechnical engineering a soil medium however presents much complexity and there are a number of factors (e.g., particle types, water/air content and pressure, etc.) besides particle volume fraction that affect the homogenized equivalent elastic moduli of a local window or so-called statistical volume element (SVE). The variations of particle volume fraction and particle types are normally anisotropic and have the correlation lengths several orders of magnitude larger than soil particle size. With a SVE size appropriately chosen that is much larger than soil particle size, and sufficiently smaller than the above correlation lengths, it is theoretically possible that the homogenized SVE moduli can be locally isotropic while varying smoothly and randomly from one SVE to the next. The anisotropy of a soil medium is reflected by the anisotropic correlation among SVEs, e.g., in this study the correlation function models a layered soil with vertical variation only.

In the case of a highly anisotropic soil, the homogenized moduli for any SVE size chosen are always anisotropic, which is considered uncommon in soils though. The elastic tensor in Eq. (9) will become anisotropic, and the correlation among tensors may be complicate dependent on the engineering approximations made about the random field. A possible approach to circumvent such an issue is to directly use the correlation function of porosity or particle types rather than of the moduli, supposed a relation between the moduli and porosity or particle types is given. In this way Eq. (9) can be conveniently evaluated, and as long as the corresponding c.o.v of the elastic moduli for an anisotropic soil is smaller than 1, von Neumann series converge and Eq. (19) can be obtained.

Numerical example

In this example the plane strain soil is 30 m deep and 120 m wide. The foundation located at the top center of the soil is 10 m wide imposing a pressure 0.2 MPa onto the soil underneath. The Young’s modulus is assumed to be a lognormal random field along the depth only, which has a mean value 50 MPa with a 40% c.o.v. The Poisson ratio is 0.3. The bottom of the soil is fixed in the depth direction, and is free in the horizontal direction. Due to the symmetry, only the right half of the soil is modeled as shown in Fig. 1, whereas the fluctuation of the Young’s modulus is also illustrated for a particular sample of the soil.

The Young’s modulus lognormal random field E is given in terms of the standard zero mean unit variance Gaussian field Y as

(34) E(x,\theta )=\exp (a+{\sigma }_{Y}Y(x,\theta )).

By letting E(x,\theta )=\overline{E}+{\sigma }_{E}Z(x,\theta ), we have the mean and variance of the Young’s modulus

(35) \overline{E}=\exp \left(a+{\displaystyle \frac{{\sigma }_Y^{\mathrm{2}}}{\mathrm{2}}}\right),

(36) {\sigma }_E^{\mathrm{2}}=\exp (\mathrm{2}a+{\sigma }_Y^{\mathrm{2}})(\exp ({\sigma }_Y^{\mathrm{2}})-\mathrm{1}),

and therefore the c.o.v = \sqrt{\exp ({\sigma }_Y^{\mathrm{2}})-\mathrm{1}} . The auto-correlation function of the centered lognormal random field Z can be expressed in terms of {\rho }_{\mathrm{2}}^Y as

(37) {\rho }_{\mathrm{2}}^Z(x)={\displaystyle \frac{\exp ({\rho }_{\mathrm{2}}^Y(x))-\mathrm{1}}{\exp ({\sigma }_Y^{\mathrm{2}})-\mathrm{1}}}.

The auto-correlation function of the underlying Gaussian field Y is specified as

(38) {\rho }_{\mathrm{2}}^Y(x)=e^{-\frac{\left|x\right|}{l_c}},

with l_c the correlation length.

The slow-scale BVP using the mean value of the Young’s modulus is solved by using the standard finite element method as shown in Fig. 2 with a 36 × 48 non-uniform mesh. This “apparent” settlement of the foundation for the slow-scale bvp is found to be 5.67 cm, which corresponds to geotechnical design considering only mean values of soil properties.

In the fast-scale BVP computation, the second term of Eq. (9a) provides an additional settlement. To verify the accuracy of the fast-scale BVP computation, we compare numerical results of four meshes of linear quadrilateral finite elements with the semi-analytical result Eq. (39) below when the correlation length l_c becomes infinitely large. Equation (17) therefore becomes

(39) \overline{\mathrm{\Delta }L\mathrm{\Gamma }\mathrm{\Delta }L}{e}_{\mathrm{0}}=L_{\sigma }\mathrm{\Gamma }{L}_{\sigma }{\rho }_{\mathrm{2}}{e}_{\mathrm{0}}=L_{\sigma }\mathrm{\Gamma }{L}_{\sigma }{e}_{\mathrm{0}}={\gamma }^{\mathrm{2}}{u}_{\mathrm{0}},

which is semi-analytical since the slow-scale solution u₀ is obtained numerically using standard finite element method. In Table 1 the four meshes are made non-uniformly, with the coarsest one 6x8 shown in Fig. 3. By taking the mesh size ratios approximately as 1, 2/3, 1/2, and 1/3, as shown in Fig. 4, for two-dimensional linear finite elements the error for the mean value computation Eq. (17) is approximately quadratic to the mesh size h, i.e., Error~O(h²).

With the above meshing verification, the additional settlement of the fast-scale BVP for the case l_c = 10 m is solved by using the four meshes, as shown in Table 2, which contributes to approximately further 14% increase. The total settlement is very close to the Monte Carlo result obtained with 10,000 samples (Fig. 5). It is noted that, in this case the mean value become stabilized after the number of Monte Carlo samples reaches about 6000, as shown in Fig. 5. In the fast-scale BVP computation, a mesh coarser than that of Monte Carlo simulation (36 × 48) can be used since the correlation function is much smoother than heterogeneities in a soil sample. Finally we remark that the computing time of the MSFEM is significantly less than the Monte Carlo method, in the 9 × 12 case equivalent to approximately 3% of computing time for the Monte Carlo method where about 2000 samples is required to achieve the same level of accuracy.

Conclusion

A fast computing procedure is developed for the MSFEM on geotechnical problems. The example shows that using mean value of the soil elastic properties underestimates the settlement, about 14% in this case. It is remarked that MSFEM can achieve further computing saving by using numerical truncation scheme on Green function and the autocorrelation function. Further results will be reported on the computing saving and computation of variance. Regarding higher moments the current method seems not competitive against Monte Carlo method yet. However, with mean and variance efficiently obtained it opens a way to quickly estimate the reliability bound of a geotechnical system.

References

[1] Xu X F. Quasi-weak and weak formulation of stochastic finite elements on static and dynamic problems: A unifying framework. Probabilistic Engineering Mechanics, 2012, 28: 103–109

[2] Xu X F. Generalized variational principles for uncertainty quantification of boundary value problems of random heterogeneous materials. Journal of Engineering Mechanics, 2009, 135(10): 1180–1188

[3] Xu X F. Stochastic computation based on orthogonal expansion of random fields. Computer Methods in Applied Mechanics and Engineering, 2011, 200(41–44): 2871–2881

[4] Xu X F, Stefanou G. Convolved orthogonal expansions for uncertainty propagation: Application to random vibration problems. International Journal for Uncertainty Quantification, 2012, 2(4): 383–395

[5] Xu X F, Chen X, Shen L. A green-function-based multiscale method for uncertainty quantification of finite body random heterogeneous materials. Computers & Structures, 2009, 87(21–22): 1416–1426

[6] Shen L, Xu X F. Multiscale stochastic finite element modeling of random elastic heterogeneous materials. Computational Mechanics, 2010, 45(6): 607–621

[7] Kröner E. Statistical continuum mechanics. Wien: Springer-Verlag, 1972

[8] Fokin G. Iteration method in the theory of nonhomogeneous dielectrics. Physica Status Solidi. B, Basic Research, 1982, 111(1): 281–288

[9] Eyre D J, Milton G W. A fast numerical scheme for computing the response of composites using grid refinement. Eur Phys J AP, 1999, 6(1): 41–47

[10] Michel J C, Moulinec H, Suquet P. Effective properties of composite materials with periodic microstructure: A computational approach. Computer Methods in Applied Mechanics and Engineering, 1999, 172(1–4): 109–143

[11] Xu X F, Graham-Brady L. Computational stochastic homogenization of random media elliptic problems using Fourier Galerkin method. Finite Elements in Analysis and Design, 2006, 42(7): 613–622

[12] Ostoja-Starzewski M. Stochastic finite elements: Where is the physics? Theoretical and Applied Mechanics, 2011, 38(4): 379–396

[13] Ostoja-Starzewski M. On the admissibility of an isotropic, smooth elastic continuum. Arch Mech, 2005, 57(4): 345–355

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg