This tutorial is intended for the full version of the toolbox.
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  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.
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:
where stands for permeability and 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 ):
with – surface tension, – contact angle and – 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 – 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:
with at the bottom boundary (contact with water surface).
Futhermore, by exploiting definition of , we can easily define residence time of liquid at any particular location in material matrix with the following time integral:
Easy determination of the residence time is important for simulation of swelling porous media, where permeability can be updated in every time step of the simulation.
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 – corrected permeatbility coefficient as follows: in area occupied by liquid and in the dry region 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:
The above formula provides that we get permeability equal to in liquid region and times greater in the dry region. Ideally, should approach infinity, which yields very big permeability and minimal pressure loss. However, big values of can lead to numerical problems, so there is a trade-off between required accuracy and numerical stability. In practice, it seems sufficient, to choose such that permeability of the dry region is an order of magnitude greater than that of wet region.
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]');
 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