Modified Rhie-Chow/PISO Algorithm for Collocated Variable Finite Volume Porous Media Flow Solvers

Fluid flow through porous media is fundamental to many natural and industrial processes, such as groundwater flows, filtration and chemical and biomass processing. In order to efficiently simulate these processes and predict their performances, robust mathematical and numerical models are of high importance. While the equations governing the flow in porous media are readily specified using the method of volume averaging [1] in conjunction with suitable closure models, it often remains challenging to obtain physically acceptable solutions in the vicinity of fluid-porous interfaces, where the porosity is discontinuous. Without special care in the algorithm development, such discontinuity may yield spurious oscillations in the solution variables. This is especially true for high Reynolds number flows and/or low porosity porous media, for which the jump in flow resistance and/or porosity is high. The occurrence of spurious oscillations is particularly pronounced when segregated algorithms are applied, in which the velocity and pressure equations are decoupled and an iterative solution process is required. Betchen et al. [2] and DeGroot et al. [3] proposed numerical schemes that avoid spurious oscillations at sharp interfaces for both structured and unstructured grids using a collocated variable, finite volume block-coupled solver, which solves pressure and velocity simultaneously. Their rather complex schemes, which treat fluid-porous interfaces consistently, work well for high Reynolds number flows, but require special treatment of the fluid-porous interfaces. In contrast to the works of Betchen et al. [2] and DeGroot et al. [3], this work focuses on more common segregated algorithms for collocated variable finite volume methods for incompressible porous media flow. A modified Rhie-Chow / PISO (Pressure-Implicit with Splitting of Operators) algorithm is proposed, which by construction avoids the development of spurious oscillations near sudden changes in effective properties. The algorithm is based on a distribution over nearby grid cells of the discontinuous flow resistance, in a similar way as was proposed by Mencinger et al. [4] for discontinuous body force fields. In order not to create time step restrictions for high Reynolds number flows or low porosities, for which the flow resistance is high, the flow resistance term is treated implicitly. The proposed algorithm is successfully verified against published data for the velocity, cf. Figure 1, and pressure fields for the flow through a porous plug at different Darcy numbers. The proposed algorithm method is also compared to a traditional algorithm for high flow resistance, to demonstrate its robustness.