《高等传热学》课程教学资源（文献资料）A stabilized stochastic finite element second-order projection method for modeling natural convection in random porous media.pdf_大学文库

A stabilized stochastic finite element second-order projection method for modeling natural convection in random porous media Xiang Ma, Nicholas Zabaras * Materials Process Design and Control Laboratory, Sibley School of Mechanical and Aerospace Engineering, 101 Frank H.T. Rhodes Hall, Cornell University, Ithaca, NY 14853-3801, USA article info Article history: Received 15 April 2007 Received in revised form 6 March 2008 Accepted 5 June 2008 Available online 21 June 2008 Keywords: Natural convection Random porous media Stochastic finite element method Stochastic projection method Sparse grid collocation Stabilized finite element method abstract We consider natural convection in flow saturated porous media with random porosity. The porosity is treated as a random field and a stochastic finite element method is developed. The stochastic projection method is considered for the solution of the high-dimensional stochastic Navier–Stokes equations since it leads to the uncoupling of the velocity and pressure degrees of freedom. Because of the porosity dependence of the pressure gradient term in the governing flow equations, one cannot use the first-order projection method. A stabilized stochastic finite element second-order projection method is presented based on a pressure gradient projection. A two-dimensional stochastic problem with moderate and large variation in the random porosity field is examined and the results are compared with Monte-Carlo and sparse grid (Smolyak) collocation approaches. Excellent agreement between these results indicates the effectiveness and accuracy of the proposed methodology. 2008 Elsevier Inc. All rights reserved. 1. Introduction Fluid flow through porous media is an ubiquitous process occurring in various applications such as fluidized beds, solidification of alloys, geothermal energy systems and oil recovery. The analysis of flow through a medium with deterministic porosity has been well studied [1]. However, in practice, only limited statistical information is available regarding the structure and material properties of the medium. These statistics are easily extracted and reconstructed from experimental data. The porosity can thus be conveniently described by random fields. This enables us to develop a methodology that treats the porosity as input uncertainty and analyzes the propagation of this uncertainty through the governing equations of thermal and flow transport. In this context, a number of methods have been proposed [2–7] that however are limited to small fluctuations and do not provide higher-order statistics of the solution. A more effective technique is the spectral stochastic finite element method (SSFEM) [8]. In this method, the random field is discretized directly, i.e., uncertainty is treated as an additional dimension along with space and time and a field variable is expanded along the uncertain dimension using suitable expansions. The most widely used expansions are the Karhunen–Loève expansion (KLE) and the polynomial chaos expansion (PCE). In PCE, Gaussian random variables are used with Hermite polynomials to represent a second-order stochastic process. This was extended to the generalized polynomial chaos expansion (GPCE), which uses the Wiener–Askey orthogonal polynomials [9]. The SSFEM has been widely used in fluid flow simulations [10–12]. 0021-9991/$ - see front matter 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2008.06.008 * Corresponding author. Tel.: +1 607 255 9104; fax: +1 607 255 1222. E-mail address: zabaras@cornell.edu (N. Zabaras). URL: http://mpdc.mae.cornell.edu/ (N. Zabaras). Journal of Computational Physics 227 (2008) 8448–8471 Contents lists available at ScienceDirect Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

During the past few years, a number of works have been reported using the GPCE method driven by a K–L expansion to discretize the random porosity [13–18]. Others used the so called KL-based moment-equation approach [19–21]. Most of these works focused on either thermal diffusion or flow motion. There is no detailed analysis as of today that accounts for both thermal and momentum transport (e.g., natural convection) in such random porous media. In order to reduce the problem complexity and decouple the calculation of velocity and pressure, a stochastic projection method was developed based on the first-order projection method and was applied successfully to natural convection in a closed cavity [22–25]. However, in porous media flow, due to the porosity dependence of the pressure gradient term in the momentum equation, we cannot impose the divergence-free constraint as is the case in the first-order projection method [26]. Thus, in order to model uncertainty propagation in natural convection in random heterogeneous porous media we need to extend the stochastic projection formulation to the second-order projection approach. This is one of the primary contributions of this paper. The projection method for the incompressible Navier–Stokes equations, also known as the fractional step method or operator splitting method, has attracted widespread popularity [27]. The reason for this lies on the uncoupling of the velocity and pressure computation. It was first introduced as the first-order projection method (also called non-incremental pressure-correction method) [28]. Later, it was extended to the second-order scheme (also called incremental pressure-correction scheme) in which part of the pressure gradient is kept in the momentum equation [29,30]. These techniques either employ spectral or finite differences techniques and it is not straightforward to extend them to a finite element interpolation. As it is known, the approximation spaces for velocity and pressure must a priori satisfy the inf–sup condition, otherwise, there will be a severe node-to-node spatial oscillation in the pressure field [31]. The first-order projection method has some pressure stability control which depends on the time step size. However, there is a severe node-to-node spatial pressure oscillation for a second-order scheme if we do not satisfy the inf–sup condition (e.g., by using a mixed finite element formulation for velocity and pressure). In order to utilize the advantage of the incremental projection method which retains the optimal space approximation property of the finite element and allows equal-order finite element interpolation, a pressure stabilized finite element second-order projection formulation for the incompressible Navier–Stokes equations has been developed [32–35]. This method mimics the stabilizing effect of the first-order projection method. It consists of introducing the projection of pressure gradient and adding the difference between the Laplacian of the pressure and the divergence of this new field to the incompressibility equation. Following this idea, we developed a framework of pressure stabilized stochastic second-order projection method. This paper is organized as follows: In Section 2, a brief review of the deterministic problem definition and second-order projection method is given. A framework for representing stochastic processes is presented in Section 3. Various issues related to modeling the uncertainties in heterogeneous porous media are detailed in Section 4. Example problems are presented in Section 5. Results are compared with those obtained through Monte-Carlo and the sparse grid collocation approach discussed in [36]. Finally, concluding remarks and future suggestions are given in Section 6. 2. Deterministic problem formulation 2.1. Problem statement A generalized non-Darcian porous medium model for natural convective flow has been developed in [1] that includes linear and non-linear matrix drag components as well as the inertial and viscous forces within the fluid. In [37], a similar model was utilized using a volume-averaged method. Here, we just present the governing equations. For detailed derivation, the interested reader may refer to the above papers. Consider a d-dimensional bounded domain D Rd with a boundary oDd S oDn. Dirichlet boundary conditions are applied on oDd, while Neumann boundary conditions are applied on oDn. The problem consists of finding the velocity v, pressure p and temperature h such that the following non-dimensional governing equations are satisfied: ov ot þ v rv ¼ Pr Da ð1 Þ 2 2 v 1:75kvkð1 Þ ð150DaÞ 1=2 2 v þ Prr2 v rp PrRaheg; ð1Þ r v ¼ 0; ð2Þ oh ot þ v rh ¼ r2 h; ð3Þ where is the porosity of the medium and eg is the unit vector in the direction of gravity. The other important non-dimensional parameters are the Prandtl number Pr, Darcy number Da and the thermal Rayleigh number Ra. Also, we assume the Boussinesq approximation is satisfied and that appropriate boundary conditions are imposed. 2.2. Pressure stabilized formulation We will discuss the pressure stabilized second-order projection method formulation based on the pressure projection. For detailed discussion and derivation, the interesting reader may refer to [33,35]. The method consists in adding to the X. Ma, N. Zabaras / Journal of Computational Physics 227 (2008) 8448–8471 8449

incompressibility equation the divergence of the difference between the pressure gradient and its projection onto the velocity space, both multiplied by algorithmic parameters defined element-wise. We take these parameters as in [34] se :¼ c1 m h2 e !2 þ c2 kvhk he 2 2 4 3 5 1=2 ; ð4Þ where m is the viscosity, he is the local size of the element e, kvhk is the local velocity and c1 and c2 are algorithmic constants, which we take as c1 ¼ 4 and c2 ¼ 2 for linear elements and c1 ¼ 16 and c2 ¼ 4 for quadratic elements. Having introduced these parameters, the continuity equation is modified as follows: sr2 p þ sr p þ r v ¼ 0; ð5Þ rp þ p ¼ 0; ð6Þ where s is the stabilized parameter as discussed before and the new auxiliary variable p is the projection of the pressure gradient rp onto the velocity space. Eq. (5) is the modified continuity equation. 2.3. Stabilized second-order projection method Let us denote with superscript n the value of each variable at the end of the nth time step and with Dt the time step. The second-order projection method corresponding to Eqs. (1), (3), (5), (6) consists of the following three major steps [35]: Step 1. Solve for the intermediate velocity vnþ1=2 in the momentum equation: 1 Dt ðvnþ1=2 vnÞ þ 2 vn rvn vn1 rvn1 ¼ Pr Da ð1 Þ 2 2 ð2vn vn1Þ 1:75kvnkð1 Þ ð150DaÞ 1=2 2 ð2vn vn1Þ þ Prr2 vnþ1=2 rpn PrRahn eg: ð7Þ In the first-order projection method, the pressure gradient term rpn is neglected in this equation. It is important to understand that the intermediate velocity vnþ1=2 does not satisfy the continuity (divergence-free) constraint. Thus, we employ the following projection step. Step 2. Projection step ðvnþ1 vnþ1=2Þ Dt þ rðpnþ1 pnÞ ¼ 0; ð8Þ sr2 pnþ1 þ sr pn þ r vnþ1 ¼ 0: ð9Þ A common approach to avoid using a mixed finite element interpolation is to use the pressure Poisson equation. In particular, we take the divergence of Eq. (8) and make use of Eq. (9) to yield ðDt þ sÞr2 pnþ1 ¼ Dtr2 pn þ r vnþ1=2 þ sr pn: ð10Þ In this equation, pn; vnþ1=2 and pn are known quantities from the previous step. Thus, we can solve for pnþ1. Note that the boundary condition for this equation are homogeneous Neumann boundary conditions [28]. In addition, the pressure at a given point (here at (0, 0)) is fixed at zero value. Step 3. Finally, we update pnþ1, velocity vnþ1 and temperature hnþ1 according to pnþ1 ¼ rpnþ1; ð11Þ vnþ1 ¼ vnþ1=2 Dtðrpnþ1 rpnÞ; ð12Þ hnþ1 hn Dt þ 2vn rhn vn1 rhn1 ¼ r2 hnþ1 : ð13Þ Note that Eqs. (7) and (10)–(13) provide a complete solution to the problem and are fully-decoupled. Thus, we can first solve for the intermediate velocity vnþ1=2 from Eq. (7) and then solve the pressure from Eq. (10). Finally, we update pnþ1, velocity vnþ1 and temperature hnþ1 according to Eqs. (11)–(13). Each of them is a rather simple equation that can be easily solved using the finite element method. This enables us to introduce next the GPCE formulation of this problem. 3. Representation of stochastic processes Consider a complete probability space ðX; F;PÞ with sample space X which corresponds to the outcomes of some experiments, the r-algebra F of subsets of X (these subsets are called events) and the probability measure P on F. The uncertainty 8450 X. Ma, N. Zabaras / Journal of Computational Physics 227 (2008) 8448–8471

is modeled as a second-order stochastic process, which is a process that has finite variance. In subsequent sections, any quantity with an x-dependence represents a stochastic quantity and x 2 X. Theoretically, the stochastic process can be represented as a random variable at each spatial and temporal location. Therefore, we require an infinity number of random variables to completely characterize a stochastic process. This poses a numerical challenge in modeling uncertainty in physical quantities that have spatio-temporal variations, hence necessitating the need for a reduced-order representation. In this section, we will consider two most popular ways of approximating a second-order stochastic process using a truncated spectral expansion comprising of a few random variables: (a) approximation by the Karhunen–Loève expansion (K–L) [8] and (b) approximation by the generalized polynomial chaos expansion (GPCE) [9,10]. 3.1. Karhunen–Loève expansion The most common approach for representing input uncertainty is the K–L expansion, which is optimal in the sense that the mean-square error of the finite representation of the process ðx;xÞ (in our problem random porosity) is minimized. Let us denote the process by ðx;xÞ and its correlation function by Rhhðx; yÞ, where x and y are spatial coordinates. By definition, the correlation function is real, symmetric and positive definite. All its eigenfunctions are mutually orthonormal and form a complete set spanning the function space to which ðx;xÞ belongs. The K–L expansion then takes the following form: ðx;xÞ ¼ ðxÞ þX1 i¼1 ffiffiffiffi ki p fiðxÞniðxÞ; ð14Þ where ðxÞ denotes the mean of the random process, and n ¼ fniðxÞg1 i¼1 forms a set of uncorrelated random variables. If the process is a Gaussian process, then they are standard identically independent Nð0; 1Þ Gaussian random variables. Also, fiðxÞ and ki are the eigenfunctions and eigenvalues of the correlation function, respectively. They are the solutions from the following eigenvalue problem: Z Rhhðx; yÞfiðyÞdy ¼ kifiðxÞ: ð15Þ We can solve this equation numerically [8,38]. The K–L expansion provides an effective way to represent the input uncertainties which have a spatial variation such as material properties when their correlation structure is known. Also note that we always truncate the expansion into finite number of terms. The number of expansion terms represents the stochastic dimensions used in the problem. The final form of the truncated K–L expansion is as follows: ðx;xÞ ¼ ðxÞ þXM i¼1 ffiffiffiffi ki p fiðxÞniðxÞ; ð16Þ where M is the truncated number of terms, which denotes the stochastic dimension in the problem. We use the scalable library SLEPc for eigenvalue problem computations [39]. SLEPc is based on the PETSc [40] data structure and it employs the MPI standard for message-passing communication for parallel high performance computation. 3.2. Generalized polynomial chaos expansion The polynomial chaos (PC) expansion technique for representation of L2-random processes was originally described in [41]. This constituted representing the random process as an expansion in terms of Hermite polynomials in the random space (a trial basis for L2ðXÞ). An extension to the original polynomial chaos expansion technique was introduced in [9], where hypergeometric orthogonal polynomials from the Askey series in the random space were used as a trial basis for L2ðXÞ. We represent the general second-order random process XðxÞ, which is a process with finite variance, as XðxÞ ¼ a0C0 þX1 i1¼1 ai1C1ðni1 ðxÞÞ þ þX1 i1¼1 Xin1 in¼1 ai1i2...inCnðni1 ðxÞ; ... ; nin ðxÞÞ þ ; ð17Þ where Cnðni1 ðxÞ; ... ; nin ðxÞÞ denote the Wiener–Askey polynomial chaos of order n in terms of the uncorrelated random vector n :¼ ðni1 ðxÞ; ... ; nin ðxÞÞ. In the original polynomial chaos, fCng are multi-dimensional Hermite polynomials and n are orthonormal standard Gaussian random variables. In the GPCE, however, fCng and n are inter-related through the joint PDF of n. For example, gamma distribution corresponds to Laguerre polynomials and uniform distribution corresponds to Legendre polynomials [9]. Since in our problem the porosity is modeled as Gaussian random filed, Hermite polynomials are used in the PC expansion. For notational convenience, Eq. (17) can be rewritten as XðxÞ ¼ X1 j¼0 a^jWjðnÞ; ð18Þ X. Ma, N. Zabaras / Journal of Computational Physics 227 (2008) 8448–8471 8451

where, the equality is interpreted in the L2ðXÞ sense and there is a one-to-one correspondence between Cnðni1 ðxÞ; ... ; nin ðxÞÞ and WjðnÞ. Since each type of polynomial in the Askey-series forms a complete basis for L2ðXÞ, we can expect the GPCE to converge to any L2 random process in the mean-square sense. The orthogonality relation of the Wiener–Askey polynomial chaos takes the form hWiWji¼hW2 i idij; ð19Þ where dij is the Kronecker delta and h; i denotes the ensemble average, which is the inner product in the Hilbert space of the variables n, hfðnÞgðnÞi ¼ Z fðnÞgðnÞWðnÞdn: ð20Þ Here, WðnÞ is the weighting function corresponding to the Wiener–Askey polynomial chaos basis Wj [9]. Note that, some types of orthogonal polynomials from the Askey scheme have weighting functions the same as the probability function of certain types of random distributions. For example, the weighting function of the p-dimensional Hermite polynomial is just the probability density function of multivariate standard normal distribution, i.e., WðnÞ ¼ 1 ffiffiffiffiffiffiffiffiffiffiffiffi ð2pÞ p p e1 2nTn : ð21Þ In practice, we then choose the type of independent variables n in the polynomials fWjðnÞg according to the type of random distribution and truncate the expansion at finite term P, i.e., XðxÞ ¼ XP j¼0 a^jWjðnÞ: ð22Þ The total number of expansion terms is (P + 1) and is determined by the dimension (M) of random vector n and the highest order (n) of the polynomials fWjg: P þ 1 ¼ ðn þ MÞ! n!M! : ð23Þ We choose a suitable order of the GPCE to capture strong non-linear dependence of the solution process on the input uncertainty (uncertainty quantification or uncertainty propagation process). Remark 1. The truncated GPCE expansion is characterized by the stochastic dimension and the order of the expansion. The stochastic dimension is determined by the number of terms M in the truncated K–L expansion of the input random processes. Since the accuracy of the truncated GPCE depends on the order of the expansion, we require techniques to determine the optimal truncation order. We use the weak-Cauchy convergence criterion for this purpose [13]. Let the guess for optimal order be q. Then we construct an order m GPCE, where m ¼ q þ 1; q þ 2. In the criterion, we require that the L2 norm of the difference in the two approximations be negligible. Note the convergence should hold point-wise and these checks are made a priori in order to determine the optimal order of the GPCE. 4. Stochastic finite element method formulation In this section, we will present the complete stochastic finite element formulation for this problem. 4.1. Non-polynomial function evaluations of stochastic spectral expansion First, let us determine the spectral expansion of the product of the form c ¼ ab ¼ PP i¼0ciWiðnÞ, where a and b are given by aðnÞ ¼ XP i¼0 aiWiðnÞ; bðnÞ ¼ XP i¼0 biWiðnÞ: ð24Þ We want to find the coefficients ck of the expression cðnÞ ¼ XP k¼0 ckWkðnÞ ¼ XP i¼0 XP j¼0 aibjWiðnÞWjðnÞ: ð25Þ Following the method introduced in [42], we perform a Galerkin projection onto the polynomial orthogonal basis and use the orthogonality of the basis discussed in the previous section. Then the expression of the coefficients can be found as ck ¼ XP i¼0 XP j¼0 hWiWjWki hW2 k i aibj; ð26Þ 8452 X. Ma, N. Zabaras / Journal of Computational Physics 227 (2008) 8448–8471

where the expectation value h; i can be evaluated by any numerical integration rule. Next, we consider a general non-linear function gðx;Þ, where is the random porosity. We need to express this function as gðx;Þ ¼ XP i¼0 giWi; ð27Þ where gi is the expansion coefficient onto the polynomial basis and the porosity ðx;xÞ is written here based on the K–L expansion as follows: ðx;xÞ ¼ ðxÞ þXM i¼1 ffiffiffiffi ki p fiðxÞniðxÞ ¼ XP i¼0 iðxÞWiðxÞ: ð28Þ Here, 0ðxÞ ¼ ðxÞ, iðxÞ ¼ ffiffiffiffi ki p fiðxÞ, for i ¼ 1; ... ; M and iðxÞ ¼ 0, for i > M. This is because the first M þ 1 term Hermite polynomials are just W0 ¼ 1; W1 ¼ n1ðxÞ; ... ; WM ¼ nMðxÞ. By writing the K–L expansion as Eq. (28) instead of Eq. (16), it is easy to formulate and perform the polynomial chaos calculations. Using the same method as before (i.e., performing a Galerkin projection onto each basis element), we can obtain from Eq. (27) the following: g x; XP i¼0 iWi !; Wj * + ¼ gjhW2 j i: ð29Þ Thus, we obtain gj ¼ hgðx; PP i¼0iWiÞ; Wji hW2 j i : ð30Þ This expression may be evaluated using quadrature rule as discussed in [43] or a non-quadrature-based approach such as the integration, Taylor series and sampling approach discussed in [42]. Here, we employ a Monte-Carlo based sampling approach based on the Latin-Hypercube sampling (LHS) strategy [44]. In this case, we first generate samples of uncorrelated standard normal variables n using LHS. The idea of LHS is to subdivide the stochastic support space of the joint PDF of n into N subintervals along each stochastic dimension and to ensure that one sample of n lies in each subinterval. For each sample, we then calculate the value of the integrand in the numerator of Eq. (30). Summing all the values, the expectation value h; i is just the arithmetic mean of these realizations. Also, the value of hW2 j i can be pre-computed using quadrature rule. 4.2. GPCE-based formulation By now, we have developed all the tools we need to formulate natural convection in random porous media. In the stochastic natural convection problem, the input uncertainties are due to the Gaussian random field of porosity. Note that since there are non-linear functions ofðx; xÞ in the governing equation (1), we need to first express them in the polynomial basis using the method discussed in the previous section: 1 ðx;xÞ ¼ XP i¼0 ^iWi; ð31Þ ð1 ðx;xÞÞ2 ðx;xÞ 2 ¼ XP i¼0 iWi; ð32Þ 1 ðx;xÞ ðx;xÞ 2 ¼ XP i¼0 ~iWi; ð33Þ where P is the number of expansion terms determined by the stochastic dimension and expansion order as in Eq. (23). Thus, the stochastic problem is to find the stochastic functions that describe the velocity field vðx;t;xÞ : D ½0; T X ! Rd , the pressure field pðx;t;xÞ : D ½0; T X ! R and the temperature field hðx;t;xÞ : D ½0; T X ! R, such that the following equations are satisfied: ovðx;t;xÞ ot þ vðx;t;xÞ ðx;t;xÞ rvðx;t;xÞ¼ Pr Da ð1 ðx;t;xÞÞ2 ðx;t;xÞ 2 vðx;t;xÞ 1:75kvðx;t;xÞkð1 ðx;t;xÞÞ ð150DaÞ 1=2 ðx;t;xÞ 2 vðx;t;xÞ þ Prr2 vðx;t;xÞ ðx;t;xÞrpðx;t;xÞ ðx;t;xÞPrRahðx;t;xÞeg; ð34Þ r vðx;t;xÞ ¼ 0; ð35Þ ohðx;t;xÞ ot þ vðx;t;xÞ rhðx;t;xÞ ¼ r2 hðx;t;xÞ: ð36Þ X. Ma, N. Zabaras / Journal of Computational Physics 227 (2008) 8448–8471 8453

Since the input uncertainty is taken as a Gaussian random field, we use Hermite polynomials to represent the solution: vðx;t;xÞ ¼ XP i¼0 viðx;tÞWiðnÞ; pðx;t;xÞ ¼ XP i¼0 piðx;tÞWiðnÞ; hðx;t;xÞ ¼ XP i¼0 hiðx;tÞWiðnÞ; pðx;t;xÞ ¼ XP i¼0 piðx;tÞWiðnÞ: ð37Þ Substitution of Eqs. (28), (31)–(33) and (37) into the stabilized second-order projection method formulation Eqs. (7) and (10)–(13), results in the following: 1 Dt XP i¼0 ðvnþ1=2 i vn i ÞWi þXP i¼0 XP j¼0 XP l¼0 ^lð2vn i rvn j vn1 i rvn1 j ÞWiWjWl ¼ Pr Da XP i¼0 XP j¼0 ið2vn j vn1 j ÞWiWj þ PrXP i¼0 r2 ðvnþ1=2 i ÞWi 1:75kvnk ð150DaÞ 1=2 XP i¼0 XP j¼0 ~ið2vn j vn1 j ÞWiWj XP i¼0 XP j¼0 irpn j WiWj PrRaXP i¼0 XP j¼0 ihn j WiWjeg; ð38Þ ðDt þ sÞr2XP i¼0 pnþ1 i Wi ¼ Dt XP i¼0 r2 pn i Wi þXP i¼0 r vnþ1=2 i Wi þ s XP i¼0 r pn i Wi; ð39Þ XP i¼0 pnþ1 i Wi ¼ XP i¼0 rpnþ1 i Wi; ð40Þ 1 Dt XP i¼0 ðvnþ1 i vnþ1=2 i ÞWi þXP i¼0 rpnþ1 i Wi XP i¼0 rpn i Wi ¼ 0; ð41Þ 1 Dt XP i¼0 ðhnþ1 i hn i ÞWi þXP i¼0 XP j¼0 ð2vn i rhn j vn1 i rhn1 j ÞWiWj ¼ XP i¼0 r2 hnþ1 i Wi: ð42Þ Then performing a Galerkin projection of each equation by h; Wki [8,10], and using the orthogonality of the polynomial basis Eq. (19), we obtain 1 Dt ðvnþ1=2 k vn k Þ þ hWiWjWlWki hW2 k i XP i¼0 XP j¼0 XP l¼0 ^lð2vn i rvn j vn1 i rvn1 j Þ ¼ Pr Da XP i¼0 XP j¼0 ið2vn j vn1 j Þeijk 1:75kvnk ð150DaÞ 1=2 XP i¼0 XP j¼0 ~ið2vn j vn1 j Þeijk þ Prr2 vnþ1=2 k XP i¼0 XP j¼0 irpn j eijk PrRaXP i¼0 XP j¼0 ihn j eijkeg; ð43Þ ðDt þ skÞr2 pnþ1 k ¼ Dtr2 pn k þ r vnþ1=2 k þ skr pn k ; ð44Þ pnþ1 k ¼ rpnþ1 k ; ð45Þ 1 Dt ðvnþ1 k vnþ1=2 k Þ þ rpnþ1 k rpn k ¼ 0; ð46Þ 1 Dt ðhnþ1 k hn k Þ þXP i¼0 XP j¼0 ð2vn i rhn j vn1 i rhn1 j Þeijk ¼ r2 hnþ1 k ; ð47Þ where eijk ¼ hWiWjWki hW2 k i , k ¼ 0; 1; ... ; P. This results in ðP þ 1Þð3d þ 2Þ decoupled deterministic equations. Remark 2. In Eq. (43), it is time consuming to evaluate the fourth-order product term 2 hW2 k i PP i¼0 PP j¼0 PP l¼0^lvn i rvn j hWiWjWlWki directly using the method discussed in Section 4.1. To simplify this calculation, we introduce an auxiliary random variable as follows: 8454 X. Ma, N. Zabaras / Journal of Computational Physics 227 (2008) 8448–8471