2D Prandtl-Meyer Tutorial

This tutorial will cover the basics of setting up and running a Prandtl-Meyer expansion fan problem. The Prandtl-Meyer expansion fan is a situation where a supersonic changes direction in such a way that the flow bends away from itself. In this problem a compressible flow with a Mach number of 2.5 flows left to right and bends around the corner, which produces the Prandtl-Meyer expansion fan. 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.

../../../_images/prandtl_velocity.png

Geometry and Boundary Conditions

The domain is a simple rectangular geometry with an edge that bends away at 15 degrees. A supersonic flow with a Mach number of 2.5 enters from the left boundary and exits through an extrapolated pressure outlet boundary on the right.

../../../_images/domain_diagram.png

Geometry and boundary conditions for the Prandtl-Meyer expansion fan tutorial.

Grid

A mixed quadrilateral cell grid was first generated using GMSH.. This was then extruded into a single-cell thick grid in the z-direction. The grid is comprised of 50,544 hexahedral cells. An image of of the grid is shown below.

../../../_images/prandtl_grid.png

Grid for the Prandtl-Meyer expansion fan 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 forward-step tutorial.
// This file is drawn in units of feet

// Geometry parameters
length = 1.1;
theta = 15;  // Bend angle of the domain downward in degrees


theta = theta * (Pi/180.0);  // Convert the degrees to radians

// Define points
Point(1) = {0.0, 0.0, 0.0};
Point(2) = {0.0, 2*Tan(theta)*length, 0.0};

Point(3) = {0.0, 1.0, 0.0};
Point(4) = {1.1, 1.0, 0.0};
Point(5) = {2.2, 1.0, 0.0};


Point(6) = {2.2, 2*Tan(theta)*length, 0.0};  // point for refinement region, bends twice as much upwards as the geometry
Point(7) = {2.2, -Tan(theta)*length, 0.0};

Point(8) = {1.1, 0.0, 0.0};
Point(9) = {1.1, 2*Tan(theta)*length, 0.0};


Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};

Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 1};
Line(9) = {2, 9};
Line(10) = {4, 9};
Line(11) = {9, 6};
Line(12) = {8, 9};

// Define surfaces now
Curve Loop(1) = {1, 9, -12, 8};
Plane Surface(1) = {1};

Curve Loop(2) = {2, 3, 10, -9};
Plane Surface(2) = {2};

Curve Loop(3) = {-10, 4, 5, -11};
Plane Surface(3) = {3};

Curve Loop(4) = {12, 11, 6, 7};
Plane Surface(4) = {4};


//---- Structured mesh specification -------------------
Transfinite Line {1, 12, -6} =  201 Using Progression 1 ; // Lower left/right boundaries
Transfinite Line {2, 10, 5} = 17 ; // Upper left/right boundaries

Transfinite Line {8, -9, -3} = 35 Using Progression 1.03;  // Upstream top/bottom boundaries
Transfinite Line {-7, 11, 4} =  201 Using Progression 1 ;  // Downstream top/bottom boundaries

Transfinite Surface {1, 2, 3,4};
Recombine Surface {1, 2, 3, 4};
//----------------------------------------------------


width = -0.05;
Extrude {0, 0, width} {
Surface{1,2,3,4};
Layers{1};
Recombine;
}

Physical Surface("Inlet", 1) = {21, 43};
Physical Surface("Top", 2) = {47, 69};
Physical Surface("Outlet", 3) = {73, 95};
Physical Surface("Wall", 4) = {33, 99};
Physical Surface("Symmetry", 5) = {1, 2, 3, 4, 56, 34, 100, 78};

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

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 -ft 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, Lref=1m, pieSlice>

boundary_conditions: <
 Inlet=supersonicInlet(v=2874.0322 ft/s, p=12 psia, T=550 R),
 Outlet=extrapolatedPressureOutlet,
 Wall=slip(),
 Top=slip(),
 Symmetry=symmetry()
>

// Initial conditions
initialCondition: <T=550 R, p=12 psia, v=2874.0322 ft/s>

// Flow properties
flowRegime: laminar
flowCompressibility: compressible

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

// Time-stepping
timeIntegrator: BDF
timeStep: 1.0e-5
numTimeSteps: 501
convergenceTolerance: 1.0e-30
maxIterationsPerTimeStep: 10

// Inviscid flux
inviscidFlux: SOU

// Gradient limiting
limiter: venkatakrishnan

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

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

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

// Pressure equation
pressureCorrectionEquationOptions: <linearSolver=HYPRE,relaxationFactor=0.3,maxIterations=50>
  pressureBasedMethod: SIMPLEC

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


// Printing, plotting and restart parameters.
print_freq: 100
plot_freq: 50
plot_modulo: 0
plot_output: pResidualTT
restart_freq: 250
restart_modulo: 0

}

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 50,000 cells. Stream has parallel efficiency down to 5,000 cells per processor, so you will want to run using less than 10 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 10 <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 10,000th time step. In the run directory extract the solution using: extract -vtk case 10000 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 time step is 10000, and the variables to extract are the pressure(P), velocity(v), density(r), temperature(t), and mach number(m).

../../../_images/prandtl_mach_result.png

Contour plot of the magnitude of the mach number in the flow field at timestep 10,000.