Solving the advection-diffusion equation in a couple of lines of code

en , MATLAB , Numerical solvers , PDE · 30 Jan 2021

Let's implement the advection equation in the form (1):

$\frac{\partial u}{\partial t} - \nabla \cdot (D \nabla u) + \nabla \cdot (\mathbf{v} u) = R$

where:

$u$ is the species concentration, and the variable in our equation;
$D$ is the diffusivity;
$v$ is the velocity field;
$R$ describes sources or sinks, and may be a function of and other parameters.

In the above equation, we have the following terms:

$\nabla \times (D \nabla u)$ describes diffusion;
$\nabla \times (v u)$ describes convection.

model = createpde();

Let's consider a simple geometry:

model.importGeometry('PlateHolePlanar.stl');

figure();

pdegplot(model, 'EdgeLabels', 'on');

MATLAB expects an equation of the form:

$m \frac{\partial^2 u}{\partial t^2} + d \frac{\partial u}{\partial t} - \nabla \cdot (c \nabla u) + au = f$

so we’ll have to reorder (1) to:

$\frac{\partial u}{\partial t} - \nabla \cdot (D \nabla u) = R - \nabla \cdot (\mathbf{v} u)$

diffusivity = 1e-2;
velocity = [0.1, 0.1];
source = 0;
specifyCoefficients(model, ...
    "m", 0, ...
    "d", 1, ...
    "c", diffusivity, ...
    "a", 0, ...
    "f", @(x, s) iRightHandSide(x, s, velocity, source));

Apply an initial condition on the central circle

initialConcentration = 1;
model.setInitialConditions(0);
model.setInitialConditions(initialConcentration, 'edge', 5);

Apply boundary conditions:

model.applyBoundaryCondition('neumann', ...
    'edge', 1:4, 'q', 0, "g", 0);

Generate the mesh and show what it'll look like

mesh = generateMesh(model, "Hmax", 0.5);

figure();

pdeplot(model);

Resolve the PDE:

timeSteps = iTimeSteps(60, 10);
result = solvepde(model, timeSteps);

Show trend at some points of the solution

figure();

subplot(2, 2, 1);

iShapshot(model, result, timeSteps, 1, initialConcentration);

subplot(2, 2, 2);

iShapshot(model, result, timeSteps, round(numel(timeSteps)*1/3), initialConcentration);

subplot(2, 2, 3);

iShapshot(model, result, timeSteps, round(numel(timeSteps)*2/3), initialConcentration);

subplot(2, 2, 4);

iShapshot(model, result, timeSteps, round(numel(timeSteps)*0.99), initialConcentration);

Visualise one particular snapshot

figure()

progressValue = 61;

numStep = round(numel(timeSteps)*progressValue/100);

iShapshot(model, result, timeSteps, numStep, initialConcentration);

Export the solution as a video

iExportResultsAsVideo("result.mp4", model, result, timeSteps, initialConcentration);

Warning: The video's width and height has been padded to be a multiple of two as required by the H.264 codec.

Helper functions

function f = iRightHandSide(~, state, velocity, source)
ux = state.ux;
uy = state.uy;
v = velocity;
convection = ux*v(1)+uy*v(2);
f = source - convection ;
end
function ts = iTimeSteps(sec, varargin)
if nargin==1
    stepsPerSecond = 60;
else
    stepsPerSecond = varargin{1};
end
nSteps = sec*stepsPerSecond;
ts = linspace(0, sec, nSteps);
end
function iShapshot(model, solution, timeSteps, stepNumber, maxConcentration)
currSolution = solution.NodalSolution(:, stepNumber);
pdeplot(model, ...
    'XYData', currSolution, ...
    'ColorMap', 'jet', ...
    'ColorBar', "off");
title("Time " + compose("%2.3f", timeSteps(stepNumber)));
colorbar();
caxis(gca(), [0, maxConcentration]);
end
function iExportResultsAsVideo(fileName, model, result, timeSteps, maxConcentration)
h = figure();
numSteps = numel(timeSteps);
% create the video writer
writerObj = VideoWriter(fileName, 'MPEG-4');
writerObj.FrameRate = 60;
open(writerObj);
for k = 1:numSteps
    iShapshot(model, result, timeSteps, k, maxConcentration);
    % Capture the plot as an image
    frame = getframe(h);
    writeVideo(writerObj, frame);
end
close(writerObj);
end