2D Supersonic Bump Tutorial

This tutorial will cover the basics of setting up and running the problem of a compressible supersonic flow moving over a bump.

The files that are needed to run this case can be extracted from a gzipped archive located here. The files will extract into a directory structure that mirrors this case’s location in the tutorials repository.

Geometry and Boundary Conditions

This case is set up as a 2D problem. A high-speed flow from left to right is developed. The supersonic flow encounters the bump and oblique shockwaves are generated.

Geometry and boundary conditions for the tutorial case.

Grid

A triangular cell grid was generated using GMSH. This grid is a single-cell thick 3D mesh with an extrusion in the z-direction. The grid is comprised of 23,611 prismatic (extruded triangle) cells. An image of the region near the particle is shown below.

The mesh that is used for this tutorial.

To generate the mesh yourself, if you have GMSH installed, you can use the following .geo control file to generate the mesh.

GMSH .geo control file for the case.

// This GMSH file is drawn in units of millimeters
// The bump has a chord length of 1000mm and the height of the bump is 4.2% of the bump chord length.

// This grid is drawn in millimeters
height_factor = 0.042; // Height of bump as a fraction of chord length
chord = 1000; // Chord length of bump
bump_height = height_factor*chord; // Height of bump

extrude_thickness = 100; // Thickness of the domain

// Center of the bump
CX = 0.0; // X coordinate of leading edge of bump
CY = 0.0; // Y coordinate of leading edge of bump

// Outer domain parameters
upstream_spacing = 1000; // Distance from the bump to the inlet
downstream_spacing = 2000; // Distance from the bump to the outlet
upper_spacing = 1000; // Distance from the bump to the top wall


// Points for the bump
Point(1) = {CX, CY, 0}; // Leading edge of bump
Point(2) = {CX + 0.5*chord, CY+bump_height, 0}; // Upper point of bump
Point(3) = {CX + chord, CY, 0}; // Trailing edge of bump

// This is the center of a circle that creates a 30 degree angle with respect to the vertical
Radius = (chord*chord + 4*bump_height*bump_height) / (8*bump_height); // Radius of the circle
Printf("Radius = %f", Radius);
Point(4) = {CX + 0.5*chord, CY-(Radius-bump_height), 0}; // Center point used to create the bump arc

// Points for rectangular outer domain
Point(5) = {CX - upstream_spacing, CY, 0}; // Upstream lower point
Point(6) = {CX - upstream_spacing, CY + upper_spacing, 0};   // Upstream upper point
Point(7) = {CX + downstream_spacing, CY + upper_spacing, 0}; // Downstream upper point
Point(8) = {CX + downstream_spacing, CY, 0}; // Downstream lower point

Circle(1) = {1, 4, 2}; // Upstream arc of bump
Circle(2) = {3, 4, 2}; // Downstream arc of bump

// Lines connecting outer domain
Line(3) = {1, 5}; // Bottom upstream
Line(4) = {5, 6}; // Upstream vertical line (inlet)
Line(5) = {6, 7}; // Top wall
Line(6) = {7, 8}; // Downstream vertical line (outlet)
Line(7) = {8, 3}; // Bottom downstream


// Surface
Curve Loop(1) = {3, 4, 5, 6, 7, 2, -1};
Plane Surface(1) = {1};


//---- Sizing field for uniform unstructured mesh near bump----
// Create a Distance field
DistanceField = 1;
Field[DistanceField] = Distance;
Field[DistanceField].CurvesList = {1, 2};
Field[DistanceField].Sampling = 100;

// Apply a Threshold field to control mesh size growth
ThresholdField = 2;
Field[ThresholdField] = Threshold;
Field[ThresholdField].IField = DistanceField; // Use the Distance field
Field[ThresholdField].LcMin = 5; // Minimum cell size
Field[ThresholdField].LcMax = 35; // Maximum cell size
Field[ThresholdField].DistMin = 0; // Distance for min cell size
Field[ThresholdField].DistMax = 1000; // Distance for max cell size

// Set the Background Field
Background Field = ThresholdField;
//-------------------------------------------------------------------

// Extrude the geometry
Extrude {0, 0, -extrude_thickness} {
   Surface{1};
   Layers{1};
   Recombine;
}

// These numbers are extracted manually with the GUI after the extrusion
Physical Surface("inlet", 1) = {23};
Physical Surface("bump", 2) = {43, 39};
Physical Surface("top_wall", 3) = {27};
Physical Surface("symmetry", 4) = {19, 35};
Physical Surface("outlet", 5) = {31};
Physical Surface("front", 6) = {1};
Physical Surface("back", 7) = {44};

Physical Volume(1) = {1};

Mesh 3;
Save "case.msh";

Once the mesh is generated, you can use the following Loci utility to create a vog file from the msh file.

msh2vog -mm case

This will create a case.vog file that can be used as the grid file for the simulation.

Run Control File Setup

The run control file for the simulation is shown below.

Run control file for the case.

// Inviscid compressible bump flow.
{

// Grid file information
grid_file_info: <file_type=VOG>

// Boundary conditions
boundary_conditions:
<
 inlet=supersonicInlet(v=[1943.6 ft/s, 0.0, 0.0], p=39.329 psi, T=558.77 R),
 outlet=extrapolatedPressureOutlet,
 bump=symmetry, 
 symmetry=symmetry,
 top_wall=symmetry,
 front=symmetry,
 back=symmetry
>

// Initial condition
initialCondition: <v=1943.6 ft/s, p=39.329 psi, T=558.77 R>

// Flow properties
flowRegime: laminar
flowCompressibility: compressible

// Equation of state and transport properties
chemistry_model: air
thermodynamic_model: curve_fit
transport_model: const_viscosity
mu: 1.0e-10
kcond: 1.0e-10

// Time-integration parameters
timeIntegrator: BDF
timeStep: 5.0e-04
numTimeSteps: 401
convergenceTolerance: 1.0e-19
maxIterationsPerTimeStep: 11

// Inviscid flux 
inviscidFlux: SLAU
slau_reconstruct: 3

// Gradient limiting
limiter: venkatakrishnan
multidimensional_limiter: mlp

// Linear solver options
linearSolverTolerance: 1.0e-04
hypreSolverName: AMG

// Momentum equation
momentumEquationOptions: <linearSolver=SGS, relaxationFactor=0.5, maxIterations=5>

// Pressure correction equation 
pressureCorrectionEquationOptions: <linearSolver=HYPRE, relaxationFactor=0.5, maxIterations=20>
  pressureBasedMethod: SIMPLEC

// Energy equation
energyEquationOptions: <linearSolver=SGS, relaxationFactor=0.5, maxIterations=5, form=totalEnthalpy>

// Printing, plotting, and restart parameters
print_freq: 100
plot_freq: 100
plot_modulo: 0
plot_output: h, vResidual, pResidual, hResidual
restart_freq: 100
restart_modulo: 0

}

In the following, we briefly discuss some of the key choices made in setting the run control file parameters for this case.

1. Due to the flow being supersonic with the prescence of shockwaves, the inviscidFlux variable is set to the SLAU scheme. The SLAU scheme is a second-order accurate flux-vector splitting scheme that is based on the Roe scheme, and it is designed to handle shockwaves.

2. The flow is essentially inviscid, so the viscosity and thermal conductivity are set to very low values. The flow is also compressible, so the flowCompressibility variable is set to compressible.

Running the Simulation

In the tutorial main directory are three files that are required to run the case.

• case.vog: the grid

• case.vars: the run control file

• air.mdl: the .mdl file that contains the species definition for air that is used in this case.

The .mdl file is shown below for reference. More details about this file can be found in the Stream user guide.

The air.mdl file.

// Model for Air as an ideal gas. Note that the unit for cp are in
// J/(mol*K).

species = {
_Air = < m=28.0189, n=2.5, href=0, sref=0, Tref=1.0, Pref=101325.0, mf=1.0,
cp=[1.0,poly(29.100709729,0.0,0.0,0.0,0.0),10000.0] > ;
} ;

reactions = {
} ;

For this case, the grid has around 23,611 cells. Stream has parallel efficiency down to 5,000 cells per processor, so you will want to run using less than 5 processes for the value passed to the -np argument in the mpirun command shown below. Using fewer processes is fine and will only result in an increased time for the solution to be generated. Execute the code using the following command.

mpirun -np 2 <path_to_stream_exec> --scheduleoutput -q solution case >& run.log_0 &

The case will run until the number of timesteps specified in the case.vars file is reached. Restart checkpoint solutions will be written to the restart/ directory at the number of timesteps specified in the case.vars file for the restart_freq variable.

Results

The extract utility is used to view the results of a simulation. For this case we view the simulation at the 400th time step. In the run directory extract the solution using: extract -vtk case 400 P v r t m . This will generate a solution directory with files in the Paraview VTK format which can be opened by Paraview. The case name is case and the timestep is 400, and the variables to extract are the pressure(P), velocity(v), density(r), temperature(t), and Mach number(m). Below is a contour plot of the pressure and Mach number fields showing the steady shockwave pattern that the bump produces in the channel.

Contour plot of the pressure and Mach number fields at the 400th time step.