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

Joumal of Computational Physics 227(0)44-847 Contents lists available at ScienceDirect Journal of Computational Physics ELSEVIER journal homepage:www.elsevier.com/locate/jcp A stabilized stochastic finite element second-order projection method for modeling natural convection in random porous media Xiang Ma,Nicholas Zabaras' Uiversiry tc USA ARTICLE INFO ABSTRACT m 6 March 2008 of the velocity and ng flow not use the first-orde essure gradient pro n.A two-di roble ate and nte-Carl suits are compare of the proposed Elsevier nc.All rights reserved. 1.Introduction d directly.i.e..uncertainty is treated as a 9.The SSFEM has been widely used in fluid flow simulations 10-12. Wiener-A +160725512z mae comelledu/(N.Zabaras) o2l078niS08ro2aeBoacrseiend

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, solid￾ification 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 struc￾ture 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

X Ma N.Zabares/Joumal of Computational Physics 27(0)844-471 849 During the past few years.a number of works have been reported using the GPCE method driven by a K-L expansion to thermal and momentum transport (natural convectio)nsuch random porous media.n order to reduce the problem exity and de 22 momentum e quation.weCame formulation to the secon ental pressure-c spectralor finite differences techniques and it is not straightforward to extend them toa finite element interpolation.Asit ity and p ressure) the which retains ibility tion following this idea we developed a fram orkof pressure stabilized stochastic second-order presented 0 2.Deterministic problem formulation 2.1.Problem statement Ageneralized non-Darcian porous medium model for natural convective flow has been developed in [1]that includes lin ear and non-l inar ma interested reade may refer to the ove paper ain DcRwith a bou irichlet mpacameoeRsae (150Da' 7.v=0. (2) 2+vT0=% (3) Darcy numb d tha 2.2.Pressure stabilized formulation

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 uncer￾tainty 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 oper￾ator 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-cor￾rection 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 veloc￾ity 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 sec￾ond-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 gradi￾ent and adding the difference between the Laplacian of the pressure and the divergence of this new field to the incompress￾ibility 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 sec￾ond-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 ap￾proach 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 lin￾ear 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-dimen￾sional 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

8450 X Ma,N.Zabaras foumal of Computational Physics 227(2008)8448-8471 incompressibility ation the div ce of the differ 1-12 (4 where is the viscosity,h is the local size of the elemente.is the local velocity and c and cz are algorithmic constants. -tV2p+tV,r+又.v=0. -Vp+π=0. (6) d the 2.3.Stabilized second-order projection method Step 1.Solve for the intermediate velocityin the momentum equation: +PrV2v+R-eVp"-cPrRaleg. Step 2.Projection step (8) -+tV.+V.=0. A2oeaChcmc2尚gacomDgsoeepmeeRaaaquieanunk (+t)2p+1=2+,v+n+t, (10) Step 3.Finally.we update,velocity y+and temperature according to 1=了p+1 (11) v1=viR -At(Vp1-Vp"). (12) +1- A +2.了0--1.了0m-1=2+1 (13) pled.Thus.we can firs 3.Representation of stochastic processes Consider a complete probability space(1.F.P)with sample space which corresponds to the outcomes of some exp iments.the a-algebraFof subsets of(these subsets are called events)and the probability measure PonF.The uncertainty

incompressibility equation the divergence of the difference between the pressure gradient and its projection onto the veloc￾ity 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 under￾stand 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 partic￾ular, 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 exper￾iments, 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

X Ma N.Zabares/Joumal of Computational Physics 27(0)844-471 is modeled as a second-order stochastic process,which is a process that has finite variance.In subsequent sections,any e rep ()the Karhunen-Loevep(K)ad( approximation by the generalized polynomial chaos expansion(GPCE)9.10 3.1.Karhunen-Loeve expansion us denote the process by()byR).whereandyare spatia coordinates By x.)+)(). (14 a set of of the rv They are therom the Ra.(x.y)fi(y)dy if(x) (15 We can solve this equation numerically (838)The K-L expansion provides an effective way to represent the input uncer tainties which have a spatial variatic The final form of the truncated K-Lexpansion is as follows x-W+之Vx@ (16) where M is the truncated number of terms,which denotes the stochastic dimension in the problem age-p 32.Generalized polynomial chaos expansion 黑o四nd where hypergeometric orthogonal polynomials from the Askey series in the random space were used as a trial basis fo ,e pe the era-orerdmpehicwth he a「(a(@a(oj+ (17) where I(().....()denote the Wiener-Askey polynomial chaos of order n in terms of the uncorrelated random vec me For notational convenience.Eq.(17)can be rewritten as (18)

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 rep￾resented 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 phys￾ical 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 spec￾tral 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 fol￾lowing 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 uncer￾tainties 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 struc￾ture 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 vec￾tor 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

8452 X Ma,N.Zabaras foumal of Computational Physics 227(2008)8448-8471 where.the equality is interpreted in the l (o)sense and there is a one-to-one con spondence between r(( (平平》=（Ψ (19) 20 W( (21) X@=∑aΨ) (22) ()ndmntee P+1- (23) ygmciCEeeppammatctkaimpmmcp威nm Remark 1.The truncated GPCE expansion is characterized by the stochastic dimension and the order of the expansion.The 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 formcab-).where a andbare given by a0=工aΨ，(.b(8=∑bΨ，( (24 We want to find the coefficients cy of the expression (25) -含20 (26)

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 uncer￾tainty (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

XMa N.Zabares/Joumal f Computational Physics()44-471 8453 where the expectation value(can be evaluated by any numerical integration rule. Next.we consider a general non-linear function g(x.).where e is the random porosity.We need to express this function (27) 28) Here.))))ornd)oMThis is beca heterm Hermite poly- nomials are justo()By writing the k-Le jection onto each basis element)we can obtain from Eq.(27)the following: (2）-明 29 Thus,we obtain 哥-gK62 (30) ( This expression may be evaluated usng discusd-quraturebased approach suchas the ased on the Latin-Hypercubestatethiscase wer be The ide. subdivide the space of then calculate the value of the integrand in the numerator of Eq.(30).Summing all the values,the expectation value ()is just the arithmetic mean of these realizations.Also.the value of(can be pre-computed using quadrature rule. 4.2.GPCE-based formulation 院 By no nedia.n the sto cussed in the pre 议可. 1 31 --含 32) e(Xω) (33) (X.0 where Pis the number of expansion terms determined by the stochastic dimension and expansion order as in Eq.(23) hus.the stochastic problem is to fir the stochas ch tha e loe aand (150Da)e(x,t,) +PrV2v(x.t.@)-c(x.t.w)Vp(x.t.@)-c(x.t.@)PrRa0(x.t.w)eg. (34 7vx,t,=0. 65 .vx.t.)-V..)-vo.t.) (36)

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 poly￾nomials 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 pro￾jection 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 sub￾intervals 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 sto￾chastic 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

8454 XMa,N.Zabaras/Joumal of Computational Physics 7(0)8448-8471 Since the input uncertainty is taken as a Gaussian random field.we use Hermite polynomials to represent: z..-...) xto-xΨ⑤.t,o-4xΨ（⑤ (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: 2n-+会会会42wW-g网m 品2会-w+容92名w-g (38) w+加玄以"男=A三+会+会 (39 "-g (40) 立名w1-9%+玄-名防=0 2w-明+22四-呀m-云g。 (42 (0dtheroonity of thepoly -0,得言若名e听- -高含名-名2a网-+n-名名 -prRaeneu (43) 0-0 (At+t)p"AtV+ 4) 暖=Vp吸l (45) 心1-+暖1-哦=0， (46) -用+客呀时g (47) wheree0.1...P.This results in (P+1)(3d+)decoupled deterministic equations. Remark 2.In Eq.(43)it is time consuming to evaluate the fourth-order product term ()directly using the method discussed in Section 4.1.To simplify this calculation,we introduce an auxiliary random variable as follows:

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

X Ma N.Zabares/Joumal of Computational Physics 27(0)844-471 8455 d-宫4-名名新m (48) such that 4-22am (49 =0.u==0 0=1 u=v=0 Free fluid Porous Medium Fg.Schematic of atural saturated variable porous medium .Comparison of th

d ¼ XP m¼0 dmWm ¼ XP i¼0 XP l¼0 ^lvn i WiWl; ð48Þ such that dm ¼ XP i¼0 XP l¼0 ^lvn i elim: ð49Þ Free fluid Porous Medium = 1 u = v =0 = 0 u = v =0 θ θ = 0, u = v =0 n ∂ θ ∂ = 0, u = v =0 n ∂ θ ∂ l ε =1 w ε Fig. 1. Schematic of natural convection in a fluid saturated variable porous medium. 0. 0 1 .2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 -1.73 -1.73 -1.73 -0.09 -0.09 -0.09 -0.09 -0.01 -0.01 -0.01 -0.01 -0.01 -5.48 -9 -8 -7 -6 -5 -4 -4 -3 -3 -2 -2 -1 -1 1 1 2 2 3 4 5 Fig. 2. Comparison of the streamlines (left) and isotherms (right). Top row: first-order projection method. Bottom row: second-order projection method. X. Ma, N. Zabaras / Journal of Computational Physics 227 (2008) 8448–8471 8455

8456 X.Ma,N.Zabaras/Joumal of Computational Physics 227 (2008)8448-8471 Soourterm of interest can be reduced to).This form is now easier to calculate,since it ony evolves third-order product terms. Remark3.In the on-linerdrate -5∑o∑o(2-1e we assume that the magnitude of the velocity is determined by the mea 60010010006 Fig.3.Contour results obtained with the second-order projection method.Top row:u and v velocity component contours.Bottom row:pressure contours 00=0.u-v-0 01 u=v-0 E(w) 77777777 777777777 Fg4.Schematic of stchastic natural vectin heterogeneous porous medium

So our term of interest can be reduced to 2 hW2 k i PP j¼0 PP m¼0dm rvn j hWjWmWki. This form is now easier to calculate, since it only evolves third-order product terms. Remark 3. In the non-linear drag term  1:75kvnk ð150DaÞ 1=2 PP i¼0 PP j¼0 ~ið2vn j  vn1 j Þeijk, we assume that the magnitude of the velocity kvnk is determined by the mean velocity vn 0. Fig. 3. Contour results obtained with the second-order projection method. Top row: u and v velocity component contours. Bottom row: pressure contours. = 1 u = v =0 = 0 u = v =0 θ θ = 0, u = v =0 n ∂ θ ∂ = 0, u = v =0 n ∂ θ ∂ ε ( ) w Fig. 4. Schematic of stochastic natural convection in a heterogeneous porous medium. 8456 X. Ma, N. Zabaras / Journal of Computational Physics 227 (2008) 8448–8471

X Ma N.Zabares/Joumal of Computational Physics 27(0)844-471 8457 Remark4.From E.(44).since the PCE terms are decoupled,the stabilized parameteris determined by the kth coefficie g.5.Eigenspe 100 right:first 20 terms) g 6.Eig 14h

Remark 4. From Eq. (44), since the PCE terms are decoupled, the stabilized parameter s is determined by the kth coefficient of the spectral expansion of the velocity components. So we denote it as sk to emphasize that it is a function of the kth coef- ficient vk according to Eq. (4). Index Eigenvalue 20 40 60 80 100 0.1 0.2 0.3 0.4 0.5 5 10 15 20 0.1 0.2 0.3 0.4 0.5 Index Eigenvalue Fig. 5. Eigenspectrum of the correlation kernel (left: first 100 terms; right: first 20 terms). x y 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0.006 0.004 0.002 0 -0.002 -0.004 -0.006 -0.008 -0.01 x y 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0.004 0.003 0.002 0.001 0 -0.001 -0.002 -0.003 x y 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0.004 0.003 0.002 0.001 0 -0.001 -0.002 -0.003 -0.004 -0.005 x y 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0.027 0.026 0.025 0.024 0.023 0.022 0.021 0.02 0.019 0.018 0.017 Fig. 6. Eigenmodes of the porosity correlation field: (a) the 1st eigenmode, (b) the 5th eigenmode, (c) the 11th eigenmode and (d) the 14th eigenmode. X. Ma, N. Zabaras / Journal of Computational Physics 227 (2008) 8448–8471 8457