2D Hypersonic Inlet Tutorial

This tutorial details how to set up, run, and analyze a simulation of a hypersonic jet inlet—an example. In this case, a high-speed compressible flow is used to simulate the complex shock-boundary layer interactions and recirculation phenomena present in hypersonic inlet designs.

In this problem a high-speed compressible flow enters from the right and encounters the hypersonic inlet geometry. Oblique shocks are created and a shock train is setup in the inlet area. 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.

#.. image:: media/shock_animation.gif # :align: center



Geometry and Boundary Conditions

The high speed flow enters the domain through the left-most boundary. It strikes the hypersonic inlet geometry and creates a series of oblique shocks which end up forming a shock train in the inlet. There are two outlet boundaries: one inside the inlet and the other one at the top-right of the domain. The top horizontal boundary is a fixed pressure boundary.

../../../_images/grid_geometry_diagram.png

Geometry hypersonic inlet. Credit: Christian Max Fischer

Grid

A mixed mesh 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 262,364 prismatic cells. An image of the region near the particle is shown below.

../../../_images/mesh_picture_1.png

A close-up of the grid near the leading edge of the geometry.

../../../_images/mesh_picture_2.png

A close-up of the grid near the inlet.

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 shock-over-sphere tutorial.

// This grid is in millimeters
// Geometry is from Fischer's thesis

// Parameters for inlet, these are angles(beta) and lengths(L)
beta_1 = 9;
L1 = 267.7;

beta_2 = 20.5;
L2 = 119.5;

beta_3 = 0.0;
L3 = 58.8;

beta_4 = -1;
L4 = 141.8;

beta_5 = 15;
L5 = 206.8;

intake_gap = 15.5; // Height of intake

domain_height = 350; // Height of domain
domain_upstream_length = 100;

// Convert angles to radians
beta_1 = beta_1 * Pi/180.0;
beta_2 = beta_2 * Pi/180.0;
beta_4 = beta_4 * Pi/180.0;
beta_5 = beta_5 * Pi/180.0;

y_1 = L1*Tan(beta_1);
y_2 = y_1 + L2*Tan(beta_2);
y_3 = y_2;
y_4 = y_3 + L4*Tan(beta_4);
y_5 = y_3 + intake_gap;
y_6 = y_5;
y_7 = y_6 + (L3+L4)*Tan(beta_5);

// Points rotate counterclockwise around the domain, roughly.
Point(1) = {0.0, 0.0, 0.0}; // Leading point of the inlet geometry
Point(2) = {L1, y_1, 0.0};
Point(3) = {L1+L2, y_2, 0.0};
Point(4) = {L1+L2+L3, y_3, 0.0};
Point(5) = {L1+L2+L3+L4, y_4, 0.0};
Point(6) = {L1+L2+L3+L4, y_6, 0.0};
Point(7) = {L1+L2+L3+L4-L5, y_5, 0.0};
Point(8) = {L1+L2+L3+L4, y_7, 0.0}; // outer outlet
Point(9) = {L1+L2+L3+L4, domain_height, 0.0};
Point(10) = {0.0-domain_upstream_length, 0.0, 0.0};
Point(11) = {0.0-domain_upstream_length, domain_height, 0.0};


// Extra points for refinement region
Point(12) = {260.4, 62.4, 0.0};
Point(13) = {0.0, 21, 0.0 };

// Connect the points
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6}; // Internal outlet
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 9}; // External outlet
Line(9) = {9, 11}; // Top boundary
Line(10) = {11, 10}; // Inlet
Line(11) = {10, 1};

// Extra lines for refinement region
Line(12) = {3,7};
Line(13) = {7, 12};
Line(14) = {12, 13};
Line(15) = {13, 1};
Line(16) = {2, 12};

// First refinement region
Curve Loop(1) = {15, 14, 16, 1};
Plane Surface(1) = {1};

// Second refinement region
Curve Loop(2) = {16, -13, -12, -2};
Plane Surface(2) = {2};

// Inlet region
Curve Loop(3) = {12, -6, -5, -4, -3};
Plane Surface(3) = {3};

// External surface
Curve Loop(4) = {-11, -10, -9, -8, -7, 13, 14, 15};
Plane Surface(4) = {4};

Transfinite Curve {12, 5, 16, 15} = 50; // Vertical lines of the refinement region

// Upper horizontal wall of the inlet region
Transfinite Curve {6} = 151; // This number must be the sum of the numbers below, minus 1
Transfinite Curve {3} = 51; // First part of the bottom of the inlet region
Transfinite Curve {4} = 101; // Second part of the bottom of the inlet region

Transfinite Curve {1, 14} = 400; // Horizontal lines of the first refinement region
Transfinite Curve {2, 13} = 200; // Horizontal lines of the second refinement region

Transfinite Surface {1} = {1, 13, 12, 2}; // Four Points of the first refinement region
Transfinite Surface {2} = {2, 12, 7, 3}; // Four Points of the second refinement region
Transfinite Surface {3} = {3, 7, 6, 5}; // Four Points of the inlet region

Recombine Surface {1, 2, 3};


width = 100;
extrude() = Extrude {0, 0, -width} {
  Surface{1, 2, 3, 4};
  Layers{1};
  Recombine;
};


// Mesh refinement-----------------
Field[0] = Distance;
//Field[0].SurfacesList = {1, 2};
Field[0].CurvesList = {15, 14, 13};
Field[0].Sampling = 2000;

Field[1] = Threshold;
Field[1].InField = 0;
Field[1].DistMax = 40;
Field[1].DistMin = 0;
Field[1].SizeMax = 10;
Field[1].SizeMin = 0.8;
Field[1].Sigmoid = 1;

Background Field = 1;
//-----------------------------------


Physical Surface("inlet", 1) = {104};
Physical Surface("walls", 2) = {29, 59, 86, 82, 74};
Physical Surface("external_outlet", 3) = {112};
Physical Surface("internal_outlet", 4) = {78};
Physical Surface("symmetry", 5) = {1, 2, 3, 4, 129, 38, 60, 87, 100, 116};
Physical Surface("farfield", 6) = {108};

Physical Volume("internal", 1) = {1, 2, 3, 4};

Mesh 3;
Mesh.Smoothing = 100;

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.
{
// Grid file information
grid_file_info: <file_type=VOG, twoDimensions>

// Mach = 7.7 
boundary_conditions: <
  inlet=supersonicInlet(v=1725.64 m/s, p= 460 Pa, T= 125 K),
  external_outlet=extrapolatedPressureOutlet,
  internal_outlet=extrapolatedPressureOutlet,
  walls=symmetry,
  farfield = fixedPressureOutlet(p=460 Pa,entrainment),
  symmetry=symmetry // Front and Back z-direction symmetry planes
>

// Initial conditions
initialCondition: <v=1725.64 m/s, p= 460 Pa, T= 125 K>

// Flow properties
flowRegime: laminar
flowCompressibility: compressible

// Equation of state and transport properties
chemistry_model: air_1s0r
transport_model: const_viscosity
mu: 1.0e-12
kcond: 1.0e-3

// Time-stepping
timeIntegrator: BDF
timeStep: 1.0e-6
numTimeSteps: 1201
convergenceTolerance: 1.0e-30
maxIterationsPerTimeStep: 11

// Inviscid flux
inviscidFlux: SLAU
slau_reconstruct: 3

// Gradient limiting
limiter: venkatakrishnan
multidimensional_limiter: mlp

// Linear solver options
linearSolverTolerance: 1.0e-03
hypreSolverName: AMG
hypreStrongThreshold: 0.5

// Diagnostics variables
//  turnoverTimeScale: used to normalize residual turn-over time (TT) values
diagnostics: <turnoverTimeScale=4.0e-4>

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

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

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

// Printing, plotting, and restart parameters
print_freq: 100
plot_freq: 100
plot_modulo: 0
plot_output: pResidualTT
restart_freq: 500
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 strong 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 strong 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_1s0r: 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_1s0r.mdl file.
// Model for Air as an ideal gas
species = {
  _Air   = < m=28.89, n=2.5, href=0, sref=0, Tref=298.0, Pref=10000.0, mf=1.0 > ;
} ;

reactions = {
} ;

For this case, the grid has around 75,000 cells. Stream has parallel efficiency down to 5,000 cells per processor, so you will want to run using less than 15 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 15 <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 1,000th time step. In the run directory extract the solution using: extract -vtk case 1000 P v r t . 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 1000, and the variables to extract are the pressure(P), velocity(v), density(r), and temperature(t). Below is a contour plot of the mach number field showing the oblique shocks and the shock train within the inlet.

../../../_images/mach_plot_timestep_1000.png

Contour plot of the Mach number field at the 1,000th time step.