Tutorial 4 – SUPG Stabilization

Posted on July 13, 2017 - Tutorials

Introduction

This tutorial shows a very easy example how to solve a 2-D advection-diffusion-reaction problem on a rectangular geometry at high Peclet numbers. It is intended to give basic introduction to using SUPG stabilization of advection diffusion problems in QuickerSim CFD Toolbox.

Let us solve the following problem

    $$\overrightarrow{a}\nabla u=\nabla\cdot(\nu\nabla u)+f$$

on a rectangle $[0,1]\times[0,0.5]$ with source term $f=1$ everywhere and the following natural and essential boundary conditions:

    $$u_{x=0}=0$$

    $$u_{x=1}=0$$

    $$\left(\frac{\partial u}{\partial n}\right)_{y=0}=0$$

    $$\left(\frac{\partial u}{\partial n}\right)_{y=0.7}=0$$

Note that the above problem definition makes it identical to the corresponding 1-D advection-diffusion problem and the solution of our 2-D problem is constant with respect to y-coordinate and identical to the 1-D problem solution with respect to x-axis. We make this choice on purpose to be able to compute error norms of the numerical solution with respect to the analytical 1-D solution given in the following form:

    $$u(x)=1+\frac{1}{a}\left(x-\frac{1-\exp(\gamma x)}{1-\exp(\gamma)}\right)$$

with

    $$\gamma=\frac{a}{\nu}$$

Let us also define mesh Peclet number as $Pe=\frac{ah}{\nu}$ . Mesh Peclet number expresses the ratio of total strength of convection on the mesh with given elemnt size h to the diffusion strength, adenotes advection velocity (in our 2-D case we will simply choose $\overrightarrow{a}=[a,0]$ everywhere).

Geometry

The Figure below shows ids of the left, right, bottom and top boundary, which we will need to use, when specifying boundary conditions.

Code

% Clear MATLAB Command Window and Workspace
clc;
clear;

% Define problem parameters
a = 1;
nu = 0.1;
f = 1;

% Read and dispaly mesh
[p,e,t] = importMeshGmsh('rectangleDomain.msh');
[p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t);
displayMesh2D(p,t);


% Create advection velocity vectors with ax = a and ay = 0
ax = a*ones(nVnodes,1);
ay = zeros(nVnodes,1);

% Assemble advection and diffusion matrix and compute right hand side
D = assembleDiffusionMatrix2D(p,t,nu);
C = assembleScalarConvectionMatrix2D(p,t,nu,ax,ay,'nosupg');
F = assembleScalarSourceTerm2D(p,t,f);

% We test now solving advection diffusion problem without SUPG ('nosupg'
% flag above). Final problem matrix is a sum of D and C
M = D + C;

% Apply boundary conditions
%
% Let us impose essential (Dirichlet) boundary conditions on left (id = 10)
% and right (id = 8) boundary. Since, bottom and top boundaries have
% homogenous Neumann conditions (natural in finite element formulation)
% we do not need to modify anything for the top and bottom walls.

[M,F] = imposeScalarBoundaryCondition2D(p,e,M,F,10,'value',0);
[M,F] = imposeScalarBoundaryCondition2D(p,e,M,F,8,'value',0);


% Solve equations
u = M\F;


% Display solution
figure(2)
displaySolution2D(p,t,u,'u-values');


% Display 1-D solution
figure(3)
scatter(p(1,:), u);
title('u-values');
xlabel('x')
ylabel('u')


% Compute analytical solution
x = p(1,:);
ua = 1+(1/a)*(x-(1-exp(a/nu*x)/(1-exp(a/nu))));

% Plot both numerical and analytical 1-D solution
figure(4)
scatter(p(1,:),u);
hold on;
scatter(p(1,:),ua);
legend('Numerical solution','Analytical solution','location','NorthWest');
title('u-values');
xlabel('x');
ylabel('u');
hold off;

 

Tags