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

Andrea

Andrea

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

where:

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

In the above equation, we have the following terms:

  • describes diffusion;
  • describes convection.
model = createpde();

Let's consider a simple geometry:

model.importGeometry('PlateHolePlanar.stl');
figure();
pdegplot(model, 'EdgeLabels', 'on');

https://i1.wp.com/picciau.net/wp-content/uploads/2021/01/figure_0.png

MATLAB expects an equation of the form:

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

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);

https://i2.wp.com/picciau.net/wp-content/uploads/2021/01/figure_1.png

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);

https://i0.wp.com/picciau.net/wp-content/uploads/2021/01/figure_2.png

Visualise one particular snapshot

figure()
progressValue = 61;
numStep = round(numel(timeSteps)*progressValue/100);
iShapshot(model, result, timeSteps, numStep, initialConcentration);

https://i2.wp.com/picciau.net/wp-content/uploads/2021/01/figure_3.png

Export the solution as a video

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

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

Leave a comment

Leave a Reply

Your email address will not be published. Required fields are marked *

Comments