Tutorial 26 – Adaptive Mesh Refinement

Introduction

This tutorial is an extension of the previous tutorial so it is highly recommended to first go thorugh the Tutorial 25 – Compressible flow over a bump. In this tutorial we will focus on mesh adaptation. We want to adapt mesh in regions where high pressure jump occurs. By doing so we will be able to get better and sharper visualization of shock wave.

Adaptive mesh refinement

The whole process of adapting the mesh is automatic. It is performed with the function AMR that uses a specific objective function, here weightedPressureJump that serves as a criterion where to adapt the grid. As the name suggests, it finds the regions of large pressure gradients and denotes them to be refined. There is one crucial step in doing AMR – after it one needs to call clear imposeCfdBoundaryCondition2D. That results from the fact that there is a static structure generate when the imposeCfdBoundaryCondition2D is called for the first time. This structure is strictly connected to the mesh topology. Hence, changing the mesh (by AMR for instance) would result in some con icts. Only a total clearing out of the BC function and fresh call of imposeCfdBoundaryCondition2D results in the correct underlying data.

Code

clear;
clear functions;
clc;


% Read and display mesh
[p,e,t]=importMeshGmsh('bump.msh');
% [p,e,t]=convertMeshToSecondOrder(p,e,t);
figure(4)
displayMesh2D(p,t);

nnodes = size(p,2);

% constant specific heat ratio of gas
refValues.kappa = 1.4;

% specific heat of gas (you may provide an arbitrary value or [] if you
% are disinterested in the temperature field (no call to
% compressibleSolution(...,'T') later in the code)
refValues.R = 287;
% reference pressure value
refValues.p = 101325;
% reference density value
refValues.rho = 1.225;

% initial field values
uini = 561.485; %x velocity given so that the Mach number at inlet is 1.65
vini = 0;
pini = 101325;
rhoini = 1.225;

% Boundary IDs explicitly stated for clarity
wallIds = [10 12];
inletIds = 9;
outletIds = 11;


% Init solution and set convergence settings
[u, convergence] = initCompressibleSolution(p,e,wallIds,refValues,rhoini,[uini,vini],pini);
nstep = 200;
maxres = 1e-5;


% The solver is strongly convergent, we suggest using large Courant numbers to eliminate oscillations.
% To obtain more accurate solutions use adaptive mesh refinement (see 'AMR')
dt = courantNumber(p,t,u,5);
du=[];


% Iterate nonlinear terms
for step = 1:nstep

    %adapt mesh at time steps 35, 45 55
    if ismember(step,[15, 25, 35,45,55])
        disp('Adapting mesh based on weighted pressure gradients... ')
        [p,e,t,u] = AMR(p,e,t,u,weightedPressureJump(p,t,u));
        disp('DONE')
        figure(4)
        displayMesh2D(p,t)
        clear imposeCfdBoundaryCondition2D
    end
    % Assemble matrix and right hand side
    [K,F] = assembleEulerProblem2D(p,t,u,refValues.kappa,dt);

    % Impose boundary conditions (wall BCs must be imposed at the end, for all walls in one function call)
    [K,F] = imposeCfdBoundaryCondition2D(p,e,t,K,F,inletIds,'farfield', u, refValues, rhoini, [uini, vini], pini);
    %The outlet BC will not actually be imposed since the flow is supersonic in the whole domain.
    [K,F] = imposeCfdBoundaryCondition2D(p,e,t,K,F,outletIds,'compressibleOutlet', u, refValues, pini);
    [K,F] = imposeCfdBoundaryCondition2D(p,e,t,K,F,wallIds,'slipwall');

   % Solve for increment in flow variables and update solution
    du=K\F;

    u=u+du;

    % Compute and plot residuals
    [stop, convergence] = computeCompressibleResiduals(K,F,du,size(p),convergence,maxres);
    plotCompressibleResiduals(convergence,2);

    % Break if converged
    if (stop)
        disp('solution converged');
        break;
    end

    disp(step)
    if(mod(step,5)==0)
        %plot the results
        figure(3)
        subplot(2,1,1)
        displaySolution2D(p,t,compressibleSolution(p,u,refValues,'p'),'pressure field');
        subplot(2,1,2)
        displaySolution2D(p,t,compressibleSolution(p,u,refValues,'Ma'),'Mach number field');
        drawnow;
    end

end
% Final results should look like below. You can see the mesh after a series of mesh adaptation steps below.

% Lastly, below is attached pressure and Mach number field obtained without refining mesh. Difference in shock
% wave resolution is easily observed.