Boundary Element Method for Transport Phenomena in Porous Medium

L. Škerget*leo@uni-mb.si; R. Jecl†renata.jecl@uni-mb.si    * University of Maribor, Faculty of Mechanical Engineering, Smetanova 17, 2000 Maribor, Slovenia
† University of Maribor, Faculty of Civil Engineering, Smetanova 17, 2000 Maribor, Slovenia

2.1 INTRODUCTION

Fluid transport phenomena in porous medium refers to the processes related to and accompanied with the transport of fluid momentum, mass and heat, through the given porous medium. These processes which are encountered in many different branches of science and technology, e.g., hydrology, geomechanics, and civil, petroleum, chemical and mechanical engineering, etc. are commonly subject to theoretical treatments which are based upon the methods traditionally developed in classical fluid dynamics. Over recent decades, fluid flows in porous medium have been studied both experimentally and theoretically. Different numerical methods were used for obtaining the solutions of some transport phenomena in porous media, e.g., the finite-difference method (FDM), finite element method (FEM), finite volume method (FVM), as well as the boundary element method (BEM). The main comparative advantage of the BEM, the application of which requires the given partial differential equation to be mathematically transformed into the equivalent integral equation representation, which is latter to be discretized, over the discrete approximative methods is demonstrated in cases where this procedure results in boundary integral equations only. This turns out to be possible only for potential problems, e.g., inviscid fluid flow, heat conduction, etc. In general, the procedure results in boundary domain integral equations and therefore several techniques were developed to extend the classical BEM. The dual reciprocity boundary element method (DRBEM) represents one of the possibilities for transforming the domain integrals into a finite series of boundary integrals, see for example Blobner et al. (2000) and Pérez-Gavilán and Aliabadi (2000). The key point of the DRBEM is the approximation of the field in the domain by a set of global approximation functions and the subsequent representation of the domain integrals of these global functions by boundary integrals. The discretization of the domain is represented only by grid points (poles of global approximation functions) in contrast to FDM meshes. However, the discretization of the geometry and fields on the boundary is piecewise polygonal, which gives the method greater flexibility over the FDM methods in coping with the boundary quantities. In the DRBEM all calculations reduce to the evaluation of boundary integrals only. Another more recent extension of the BEM is the so-called boundary domain integral method (BDIM), see Škerget et al. (1989, 1999) and Jecl et al. (2001). Here, the integral equations are given in terms of the variables on the integration boundaries as well as within the domain of the integration. This situation arises when we are dealing with strongly nonlinear problems with prevailing domain-based effects, for example diffusion–convection problems. The Navier–Stokes equations are commonly used as a framework for the solution of such a problems since they provide a mathematical model of physical conservation laws of mass, momentum and energy. The velocity–vorticity formulation of these equations results in the computational decoupling of the kinematics and kinetics of the fluid motion from the pressure computation, see Wu (1982). Since the pressure does not appear explicitly in the field functions conservation equations, the difficulty connected with the computation of the boundary pressure values is avoided. The main advantage of the BDIM, as compared to the classical domain type numerical techniques, is that it offers an effective way of dealing with boundary conditions on the solid walls when solving the vorticity equation. Namely, the boundary vorticity in the BDIM is computed directly from the kinematic part of the computation and not through the use of some approximate formulae. One of the few drawbacks of the BDIM are considerable CPU time and memory requirements, but they can be drastically reduced by the use of a subdomain technique, see Hriberšek and Škerget (1996). Convection-dominated fluid flows suffer from numerical instabilities. In domain-type numerical techniques upwinding schemes of different orders are used to suppress such instabilities while in BDIM the problem can be avoided by the use of Green’s functions of the appropriate linear differential operators which results in a very stable and accurate numerical description of coupled diffusion–convective problems. There are no oscillations in the numerical solutions, which would have to be eliminated by using some artificial techniques, e.g., upwinding, artificial viscosity, as is the case with other approximation methods.

2.2 GOVERNING EQUATIONS

Due to the general complexity of the fluid transport process in porous medium, our work is based on a simplified mathematical model in which it is assumed that:

• the solid phase is homogeneous, non-deformable, and does not interact chemically with respect to the fluid,

• the fluid is single phase and Newtonian; its density does not depend on pressure variations, but only on variations of the temperature,

• the two average temperatures, Ts for the solid phase and Tf for the fluid phase are assumed to be identical and the porous medium is in thermodynamic equilibrium, that means it is described by a single equation for the average temperature T = Ts = Tf,

• no heat sources or sinks exist in the fluid; thermal radiation and Rayleigh dissipation are negligible,

• the natural convection effect is considered by using the Boussinesq approximation, where the temperature influence on the density is considered only in the term describing the body force, while in all the other terms the density is assumed to be constant.

Under these assumptions, having in mind that we are dealing with time dependent incompressible viscous fluid flow through porous medium, we may write the macroscopic set of equations, commonly called the modified Navier–Stokes equations or the Brinkman equations, as follows:

• continuity equation

vj∂xj=0,

• momentum equation

ϕ∂vi∂t+1ϕ2∂vjvi∂xj=−1ρ∂P∂xj+Fgi−γKvi⏟Darcylaw+γϕ∂2vi∂xj∂xj⏟Brinkmanextension,

• energy equation

∂T∂t+∂vjT∂xj=∂∂xja∂T∂xj,

where the vector field functions vi, gi and xi are the filtration velocity, gravity and position, respectively. The scalar quantity P = p − ρgiri is the modified pressure, while ρ and ϕ stand for the density and porosity. The quantities γ, K, a and σ are the kinematic viscosity coefficient, permeability, thermal diffusivity and the heat capacity ratio, which is defined by the expression

=ϕ+ρscp,sρcp1−ϕ,

where the subscript s denotes the solid part of porous medium. The following general expression for the density variation can be used: ρ = ρ0 (1 + F), where the function F is frequently defined as...