Tutorial 3 – Flow in Porous Media

Posted on July 13, 2017 - Tutorials

Introduction

Problem definition In this tutorial we show a very easy example how to solve fluid flow in porous media with QuickerSim Toolbox for MATLAB®. Let us first define the problem. We deal with an incompressible fluid flow with pressure field p and velocity field $\overrightarrow{u}$ subject to incompressibility constraint $\nabla\cdot\overrightarrow{u}=0$ . We exploit Darcy’s law as the mathematical model of flow in porous zone. Darcy’s law states that the local velocity vector is defined by:

    $$\overrightarrow{u}=-\frac{\kappa}{\mu}\nabla p$$

where $\kappa$ stands for permeability and $\mu$ for fluid dynamic viscosity. Substitution of the second equation into divergence constraint yields the following Laplace’s (diffusion) equation for pressure field:

    $$\nabla\cdot(\frac{\kappa}{\mu}\nabla p)=0$$

Boundary conditions

We need to supplement the above equation with appropriate pressure boundary conditions. When simulating flows in porous media it is reasonably to define a zero level reference pressure at outlet boundary, thus having:

    $$p_{outlet}=0$$

a certain pressure head at inlet (if we typically assume that the flow takes place in gravitational field from top to bottom through a porous zone and the total fluid pressure head can be assumed equal to the hydrostatic pressure). Hence, we have:

    $$p_{inlet}=\rho_{fluid}gh$$

Eventually, we assume free flow (symmetry boundary condition or in other words no viscous tension) at the side walls of the geometry. This boundary condition is coincident with the homogeneous Neumann condition for pressure, i.e.:

    $$\left(\frac{\partial p}{\partial n}\right)_{sidewalls}=0$$

Code

Clear MATLAB Command Window and Workspace

clc;
clear;

% Read and convert mesh to second order
[p,e,t] = importMeshGmsh('porous.msh');
[p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t);

% Define water dynamic viscosity [Pa*s]
mu = 0.001;


% Generate random permeability field of sand
K = (1e-9)*rand(nVnodes,1);
figure(1)
displaySolution2D(p,t,K,'Porous media permeability');


% Assemble global problem matrix accounting for Darcy law
[M, F] = assembleDiffusionMatrix2D(p,t,K);


% Set high water pressure head at inlet (rho*g*h) with h about 5.5 m
[M,F] = imposeScalarBoundaryCondition2D(p,e,M,F,8,'value',1000*9.81*5.5);


% Set zero pressure at outlet
[M,F] = imposeScalarBoundaryCondition2D(p,e,M,F,9,'value',0);


% Solve for pressure field
pressure = M\F;


% Display pressure field
figure(2)
displaySolution2D(p,t,pressure,'Pressure field');


% Compute horizontal vx and vertical vy velocity
gradp = solutionGradient2D(p,t,pressure);
vx = -K/mu.*gradp(:,1);
vy = -K/mu.*gradp(:,2);

% Compute velocity magnitude
vmag = sqrt(vx.^2+vy.^2);

% Display velocity field
figure(3)
displaySolution2D(p,t,vmag,'Velocity magnitude [m/s]');

Tags