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