Tutorial 25 – Supersonic Flow with a Shockwave

Posted on July 17, 2017 - Tutorials

This tutorial is intended for the full version of the toolbox.

Introduction

In this tutorial we present how to simulate inviscid (i.e. $\mu = 0$ ) compressible flows in 2D. In order to do so, we have to solve Euler problem and continuity equation.

    $$\frac{\partial u}{\partial t}+(u\cdot\nabla)u=-\frac{1}{\rho}\nabla p$$

    $$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho u)=0$$

Those equations are appended by the equation for energy ( $e = c_V T$ ) and the equation of state (perfect gas).

    $$\frac{\partial e}{\partial t}+\nabla\cdot(ue)=-\frac{p}{\rho}\nabla\cdot u$$

    $$p=\rho RT$$

Geometry

We will investigate flow over a bump. Geometry is shown below.

We will prescribe ‘slipwall’ boundary condition on lines with physical ids 10 and 12. Remember that the flow has no viscosity, thus there is slip on those walls. We will prescribe ‘farfield’ boundary condition at the inlet (physical id 9). Finally, ‘compressibleOutlet’ boundary condition is the utlet BC (id 11). In this simulaiton the flow is actually supersonic and the ‘compressibleOutlet’ is unnecessary (no infor- mation is traveling upstream), but it is recommended to always use this boundary condition in compressible flow simulations.

Code

clc;
clear;


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


nnodes = size(p,2);

% Notice that we do not update the mesh to the 2nd order! Simply speaking, 1st order schemes are more stable
% in calculations. You can try to run the 2nd order mesh, however.

% Reference values
%
%The solver requires reference values in order to nondimensionalize the
%variables. We recommend the reference values be the same as far field/free
%stream values. The structure fields must be named as below.


%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 tutorial 26)
dt = courantNumber(p,t,u,5);
du=[];


% Iterate nonlinear terms
for step = 1:nstep

    clear CfdBoundaryCondition2D

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

    % The resulting system of equations is symmetric positively-definite,
    % for large problems you may consider using the numerically faster
    % preconditioned conjugate gradient method. If you choose to do so,
    % please use the template below or see the pcg help page for more
    % information.
    %
    %     L = speye(4*nnodes); %in this example we use no preconditioner
    %     du=pcg(K,F,1e-8,1e4,L,L',du);

    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

   %plot the results
    figure(1);
    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

% After solution is converged your pressure field and Mach number field
% should look like below.

 

Tags