Tutorial 11 – Wicking of Porous Medium

Posted on July 17, 2017 - Porus Media, Tutorials

This tutorial is intended for the full version of the toolbox.

Introduction

This tutorial is intended to show one of the possible approaches to simulate wicking process in QuickerSim CFD Toolbox. The tutorial makes use of several concepts and equations presented in [1] and is loosely coupled to the paper, since we introduce a couple of customized concepts, which make it easier to implement the underlying mathematical model in numerical environment.

Problem Definition

We are interested in simulating wicking of a liquid absorbing material. We will study this process in a 2-D geometry shown below

At the beginning of simulation we assume that the bottom edge of the material is immersed in water (water surface exactly touches lower edge of geometry). Futhermore, water flow in the porous material can be modelled by means of Darcy’s law, leading to the following equation:

    $$\nabla \cdot (K \nabla P) = 0$$

where $K$ stands for permeability and $P$ for difference between actual pressure value and hydrostatic pressure head. The flow is driven by capillary pressure acting at the liquid front. Capillary pressure can be described by (see [1]):

    $$p_c = \frac{2 \gamma cos(\theta)}{R_c}$$

with $\gamma$ – surface tension, $\theta$ – contact angle and $R_c$ – capillary radius.

Proper boundary conditions for the pressure problem are the following:

  • Bottom boundary: Atmospheric pressure
  • Liquid front: Atmospheric pressure minus capillary pressure

In addtion to the concepts mentioned above, we define some auxiliary fields. First, we define $\phi$ – a field variable specifying water fraction (it is equal to 1 in the part of domain already occupied by liquid and 0 in the dry material). This concept allows us to track evolution of wicking process by observing the liquid fraction field. Evolution of liquid fraction is governed by the following advection equation:

    $$\frac{\partial \phi}{\partial t} + \vec{v} \nabla \phi = 0$$

with $\phi = 1$ at the bottom boundary (contact with water surface).

Futhermore, by exploiting definition of $\phi$, we can easily define residence time of liquid at any particular location in material matrix with the following time integral:

    $$t_{res} = \int_{t=0}^{t=T}{\phi d \tau}$$

Easy determination of the residence time is important for simulation of swelling porous media, where permeability $K$ can be updated in every time step of the simulation.

Numerical Approach

A major numerical difficulty stems from the fact that region occupied by liquid changes from one time step to another. One, very difficult possibility would be to generate a seperate domain mesh for every time step of the simulation. To avoid this tiduous and difficult in practical realization task, we keep one constant domain mesh forever, apply the capillary pressure at the top of the domain and only modify the governing equation resulting from Darcy’s law to correct our mistake that we do by this approximation. Let us summarize the mistake: Capillary pressure should be applied at the liquid front, we impose it at the top of the domain. Pressure loss in a porous media flow should only be present in the region occupied by liquid. Obviously, region of dry material with no liquid flow does not generate drag. To correct our error, we need to make apply a correction to the pereability coefficient. Let us define lowercase $k$ – corrected permeatbility coefficient as follows: $k = K$ in area occupied by liquid and in the dry region $k$ is much bigger to generate only minimal pressure losses in the region, where no pressure loss should take place. In practice, we make the following definition by linear interpolation between wet and dry region:

    $$k = \phi K + (1-\phi) \alpha K$$

The above formula provides that we get permeability equal to $K$ in liquid region and $\alpha$ times greater in the dry region. Ideally, $\alpha$ should approach infinity, which yields very big permeability and minimal pressure loss. However, big values of $\alpha$ can lead to numerical problems, so there is a trade-off between required accuracy and numerical stability. In practice, it seems sufficient, to choose $\alpha$ such that permeability of the dry region is an order of magnitude greater than that of wet region.

Code

In this section we present the whole MATLAB script, which solves the problem in QuickerSim CFD Toolbox. We provide most of the description in comments in the code.

clc;
clear;


% Import, translate mesh by 5 units in y-direction and display mesh
[p,e,t] = importMeshGmsh('porous.msh');
p(2,:) = p(2,:) + 5; % p(2,:) is a vector with y-coordinates of all nodes
displayMesh2D(p,t);


% Define fields
phead = initScalarSolution(p,0);
phi = initScalarSolution(p,0);
resTime = initScalarSolution(p,0);


% Define timstepping parameters (number of steps and time step size)
dt = 20;
nstep = 80;


% Initialize phi with 1 at the bottom edge (Physical Line id = 9)
bottomNodes = extractNodeIdsOnEdges(e,9);
phi(bottomNodes) = 1;


% Define physical properties
rho = 1000; % water density
g = 9.81; % graviational acceleration
gamma = 0.0723; % surface tension
theta = 0; % contact angle
Rc = 5e-6; % capillary radius
pc = 2*gamma*cos(theta)/Rc; % capillary pressure
K = 5e-8; % porous zone permeability
alpha = 10; % permeability of the dry zone one order of magnitude greater


% Init movie with evolution of liquid fraction
movieObj = initMovie('liquidFraction.avi',16);


% Precompute number of nodes and mass matrix for transient advection
% problem
nPnodes = size(p,2);
M = assembleMassMatrix2D(p,t,1);


% Begin timestepping
for step = 1:nstep
    time = step*dt;

    % Solve pressure field problem with corrected permeability
    k = K*phi + (1-phi)*alpha*K;
    [D,F] = assembleDiffusionMatrix2D(p,t,k);

    % Impose boundary conditions
    [D,F] = imposeScalarBoundaryCondition2D(p,e,D,F,9,'value',0);
    [D,F] = imposeScalarBoundaryCondition2D(p,e,D,F,8,'value',-pc);

    % Solve pressure problem
    phead = D\F;

    % Compute velocity field from Darcy's law
    pvel = phead;
    vel = -repmat(k,1,2).*solutionGradient2D(p,t,pvel);

    % Make a correction to the velocity field (pressure field in this
    % formulation does not account for hydrostatic head - hence, y-velocity
    % field must be adjusted by -K*rho*g)
    vel(:,2) = vel(:,2)-K*rho*g;

    % Solve advection problem for phi
    C = assembleScalarConvectionMatrix2D(p,t,1e-3,vel(:,1),2*vel(:,2),'supgDoublyAsymptotic');
    C = C + M/dt;
    Fphi = M*phi/dt;

    % Define phi = 1 at the bottom edge (contact with water container)
    [C,Fphi] = imposeScalarBoundaryCondition2D(p,e,C,Fphi,9,'value',1);

    % Solve equations
    phi = C\Fphi;

    % Increase numerical stability by clipping values of phi in physically
    % permissible range
    phi(phi>1) = 1;
    phi(phi<0) = 0;

    % Update residence time
    resTime = resTime + phi*dt;

    % Shoot movie frame
    figure(2)
    displaySolution2D(p,t,phi,'\phi [-]',[0 1]);
    drawnow;
    writeVideo(movieObj, getframe(gcf));

    % Report simulation progress (display time step number) - uncomment
    % disp(step)
end
Display solution

close(movieObj);

figure(2)
displaySolution2D(p,t,phi,'Liquid Fraction \phi [-]',[0 1]);


figure(3)
displaySolution2D(p,t,phead,'Pressure Field');


% Display residence time
figure(4)
displaySolution2D(p,t,resTime,'Residence Time [s]');

References

[1] R. Masoodi, H. Tan, K. M. Pillai, “Numerical Simulation of Liquid Absorption in Paper-Like Swelling Porous Media”, American Institute of Chemical Engineering Journal, Vol 58, No 8, DOI: 10.1002/aic.12759

Tags