Today, we are showing a simple script which can be run in our QuickerSim CFD Toolbox to solve an underground fluid flow in MATLAB or simply to simulate fluid flow in porous media. We exploit the Darcy law here. The complete pdf tutorial for CFD problems including porous media (but not only) can be downloaded from the homepage of this fluid dynamics toolbox.
The resulting pressure and velocity fields can be plotted with the displaySolution2D function.
Example CFD code to solve porous media flow in MATLAB:
% Flow in porous media in MATLAB (QuickerSim CFD Toolbox for MATLAB Version 1.0)<br>
% Import mesh and convert to second order mesh
[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');
% Extract ids of boundary
inletNodes = extractNodeIdsOnEdges(e,8);
outletNodes = extractNodeIdsOnEdges(e,9);
% 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(inletNodes,:) = 0;
M(inletNodes,inletNodes) = eye(length(inletNodes));
F(inletNodes) = 1000*9.81*5.5;
% Set zero pressure at outlet
M(outletNodes,:) = 0;
M(outletNodes,outletNodes) = eye(length(outletNodes));
F(outletNodes) = 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]');