3D SSAM Gen-5 Fighter Tutorial

This tutorial covers the basics of setting up and running the case of a fifth generation high-performance aircraft. This geometry is part of the Sydney Standard Aircraft Models (SSAM) that are openly available to the research community. 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.

This tutorial is based on a configuration published with the Gen-5 SSAM geometry in a paper that can be found here. The experiments in the paper examine the lift and drag coefficients on the model of the aircraft at varying angles of attack.

Geometry and Boundary Conditions

The geometry and boundary conditions for this case are shown below. All boundaries except the exit, symmetry, and plane boundaries are specified with a fixed velocity in order to set the angle of attack of the aircraft. The exit boundary is assigned a fixed pressure and the plane boundary is assigned as an no-slip surface.

Geometry and boundary conditions for the 3D SSAM fighter tutorial domain.

The inlet velocity magnitude was 20 m/s. The angle of attack of the jet was 10 degrees for this case. The pressure of the outlet was set to be 101325 Pascals.

Grid

A surface grid for the plane with 229,342 quadrilateral faces was first generated using Capstone. A selection of images for the plane surface mesh are shown below. This surface mesh and the surface meshes for the external boundaries were passed to the AFLR volumetric mesh generator. The initial boundary layer spacing on the plane surface was set to 1.0e-5 meters. The final volume mesh comprised 14,157,595 cells: 3,304,171 hexahedral elements, 198,460 pentagonal elements (5 sides), 844,114 pentagonal elements (6 sides), and 9,810,850 tetrahedral elements.

The surface mesh of the front of the plane.

The surface mesh of the wing of the plane.

The surface mesh of the rear of the plane.

The surface mesh of the underside of the plane.

Run Control File Setup for Laminar Simulation

A simulation was first performed using laminar flow in order to obtain an initial condition for the turbulent simulation. The run control file for the laminar simulation is shown below.

Run control file for the laminar case.

{
// Grid file information
grid_file_info: <Lref=1m>

boundary_conditions: <
 inlet=incompressibleInlet(v=[19.69615506 m/s, 0.0 m/s, 3.472963553 m/s]),
 outlet=fixedPressureOutlet(p=1 atm),
 plane=noslip,
 rear=incompressibleInlet(v=[19.69615506 m/s, 0.0 m/s, 3.472963553 m/s]),
 bot=incompressibleInlet(v=[19.69615506 m/s, 0.0 m/s, 3.472963553 m/s]),
 top=incompressibleInlet(v=[19.69615506 m/s, 0.0 m/s, 3.472963553 m/s]),
 symmetry=symmetry
>

// Initial conditions
initialCondition: <rho=1.225 kg/m/m/m, p=1 atm, v=0.0 m/s>

// Flow properties
flowRegime: laminar
flowCompressibility: incompressible

// Transport properties
transport_model: const_viscosity
mu: 1.789e-05

// Time-stepping
timeIntegrator: BDF
timeStep: 1.0e+30  // Value used for steady-state simulations
numTimeSteps: 1001
convergenceTolerance: 1.0e-30
maxIterationsPerTimeStep: 1

// Inviscid flux
inviscidFlux: FOU

// Gradient limiting
limiter: venkatakrishnan

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

// Diagnostics variables
// turnoverTimeScale: used to normalize residual turn-over time (TT) values
//   Plane length is 0.75 m. Velocity is 20 m/s. Time scale is (0.75/20).
diagnostics: <turnoverTimeScale=3.75e-02>

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

// Pressure equation
pressureCorrectionEquationOptions: <linearSolver=HYPRE, relaxationFactor=0.1, maxIterations=20>
  pressureBasedMethod: SIMPLEC
  PclipMin: 80000.0
  PclipMax: 140000.0

// Printing, plotting, and restart parameters
print_freq: 10
plot_freq: 1000
plot_modulo: 0
plot_output: pResidualTT
restart_freq: 1000
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. This tutorial is concerned with the 10 degree angle of attack configuration. To simulate with this without having to rotate the plane, the inflow boundary velocity is adjusted to direct the inflow at an angle relative to the plane. The specification of the velocity in the cartesian coordinate system [19.69615506 m/s, 0.0 m/s, 3.472963553 m/s] is computed from the relation V = V_inf * [cos(α), 0, sin(α)], where V_inf = 20 m/s and α = 10 degrees.

2. This laminar simulation is run to setup an initial pressure and velocity field for the subsequent turbulent simulation. The initial condition is set to a zero velocity flow and constant pressure of 101325 Pascals. The simulation is run for a sufficient number of time steps to allow the flow to develop and the initial transient to damp out.

3. The hypreSolverOutput option is set to on to provide information about the convergence of the HYPRE linear solver for the pressure-correction equation. This information can be used to adjust the linear solver parameters if necessary. Information can be found by doing a tail -f debug/debug terminal command when running on a single processor, or a tail-f debug/debug.0 command when running in parallel. This option should be turned off for long simulations with many thousands of time steps to prevent the debug files from becoming excessively large.

4. Since we are only interested in getting an initial approximate condition for the pressure and velocity fields, we can set the inviscidFlux to FOU. This will give us a baseline spatially low-order solution that we can then restart from with higher order spatial flux schemes.

5. Relaxation factors for the governing equations have been set to relatively low values to overcome the poor initial condition.

Running the Laminar Simulation

The grid file that is required to run the laminar and turbulent cases is located in the tutorial main directory.

• case.vog: the grid

Also in this directory is a sub-directory called laminar/ , which contains the run control file for the laminar case. To run the laminar case, symbolically link to the grid file in the tutorial main directory by executing the following command from within the laminar/ directory. This is done to save disk space when working with large grid files.

ln -s ../case.vog .

For this case, the grid has around 14 million cells. Stream has good parallel efficiency down to 5,000 cells per process, so you will want to run using less than 2800 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. This case was run on 80 processes and took about 1.5 hours. To run the case, enter the laminar/ directory and execute Stream using the following command.

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

The laminar 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. In this case we are only writing one restart solution at the final timestep. Once the laminar case is completed, the laminar solution can be used as an initial condition for the turbulent case. Details about restarting a turbulent case from a laminar solution are discussed in the turbulent case section.

Laminar Results

The laminar simulation is only being used as an initial condition for the turbulent simulation, so we do not show any results here. If you wish to view the laminar solution, you can do so by using the extract utility.

The extract utility is used to view the results of a simulation. For this case you can view the simulation at the 1,000th time step. In the run directory, use the following command: extract -vtk case 1000 v P . 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 1000, and the variables to extract are the velocity(v) and pressure(P). Other variables are not needed because this is an incompressible laminar simulation.

Run Control File Setup for Turbulent Simulation

The turbulent simulation restarts from the laminar solution that was obtained in the previous section. The run control file for the turbulent simulation is shown below.

Run control file for the turbulent case.

{
// Grid file information
grid_file_info: <Lref=1m>

boundary_conditions: <
 inlet=incompressibleInlet(v=[19.69615506 m/s, 0.0 m/s, 3.472963553 m/s], k=1.307e-5, omega=88.3),
 outlet=fixedPressureOutlet(p=1 atm),
 plane=noslip,
 rear=incompressibleInlet(v=[19.69615506 m/s, 0.0 m/s, 3.472963553 m/s], k=1.307e-5, omega=88.3),
 bot=incompressibleInlet(v=[19.69615506 m/s, 0.0 m/s, 3.472963553 m/s], k=1.307e-5, omega=88.3),
 top=incompressibleInlet(v=[19.69615506 m/s, 0.0 m/s, 3.472963553 m/s], k=1.307e-5, omega=88.3),
 symmetry=symmetry
>

// Initial conditions
initialCondition: <rho=1.225 kg/m/m/m, p=1 atm, v=0.0 m/s, k=1.307e-5, omega=88.3>

// Flow properties
flowRegime: turbulent
flowCompressibility: incompressible

// Transport properties
transport_model: const_viscosity
mu: 1.789e-05

// Time-stepping
timeIntegrator: BDF
timeStep: 1.0e+30
numTimeSteps: 3001
convergenceTolerance: 1.0e-30
maxIterationsPerTimeStep: 1

// Inviscid flux
inviscidFlux: SOU
turbulenceInviscidFlux: SOU

// Gradient limiting
limiter: venkatakrishnan

// Linear solver options
linearSolverTolerance: 1.0e-03
hypreSolverName: AMG
hypreStrongThreshold: 0.9
hypreSolverOutput: on

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

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

// Pressure equation
pressureCorrectionEquationOptions: <linearSolver=HYPRE, relaxationFactor=0.3, maxIterations=20>
  pressureBasedMethod: SIMPLEC
  PclipMin: 80000.0
  PclipMax: 140000.0

// Turbulence equation
turbulenceEquationOptions: <model=menterSST2003, linearSolver=SGS, relaxationFactor=0.5, maxIterations=5, crossDiffusion=standard>
  kFreestream: 1.307e-5
  omegaFreestream: 88.3

// Printing, plotting, and restart parameters
print_freq: 10
plot_freq: 4000
plot_modulo: 0
plot_output: pResidualTT, viscosityRatio, kResidualTT, omegaResidualTT, k, omega
restart_freq: 4000
restart_modulo: 0



}

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

1. While this case is being restarted from a laminar solution, note that we must still specify an initial condition in the run control file to circumvent Stream’s run control file consistency checking coding from returning an error. Upon restart, however, the values specified in the variable initialCondition will be ignored, except for the new turbulent variables (k and omega).

2. For the inviscid fluxes, we have chosen SOU for both the inviscidFlux and turbulentInviscidFlux options to achieve a high spatial accuracy.

3. For the pressure-correction equation, one should always use the HYPRE linear solver. The parameters in the linear solver section apply only to the HYPRE solver and are not used by the SGS linear solver which is being employed here for the momentum and turbulence equations. Here we select the Algebraic Multi-Grid (AMG) solver, which is superior in general to the GMRES solver in HYPRE. The hypreStrongThreshold parameter is set to a relatively high value (range 0-1,default is 0.5) to further minimize the work per linear solver iteration.

Running the Turbulent Simulation

As mentioned before, the grid file that is required to run the turbulent case is located in the tutorial main directory.

• case.vog: the grid

Also in this directory is a sub-directory called turbulent/ , which contains the run control file for this case.

To run the turbulent case, once again, make a symbolic link to the case.vog file from within turbulent/ directory. The command to make the symbolic link is shown below.

ln -s ../case.vog .

In addition, you will create a new directory in the turbulent/ directory called restart/. This directory is where you will place a set of laminar restart files for the turbulent simulation to restart from. For this case, the laminar simulation was run to 1,000 timesteps and a set of restart files were written at the 1,000th timestep. Go to the laminar/ directory and perform a copy like the one shown below.

cp restart/*1000* ../turbulent/restart/

For this case, the grid has around 14 million cells. Stream has parallel efficiency down to 5,000 cells per processor, so you will want to run using less than 2800 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. This case was run on 80 processors and took about 2.5 hours. Enter the turbulent/ directory and execute Stream using the following command.

• mpirun -np 80 <path_to_stream_exec> --scheduleoutput -q solution case 1000 >& run.log_1000 &

The turbulent case will run an additional number of timesteps from the restart as specified by the numTimeSteps variable in the case.vars file. The case will thus terminate at timestep 4,000. 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. The addition of the integer argument ``1000 in the command above tells Stream to use the solution in the restart/ directory as the initial condition. In the command above, the file that the log information is redirected to is named run.log_1000. This suffix is used to differentiate the log file from other log files that are generated when a user is restarting cases. The suffix notation prevents log file information from previous simulations from being overwritten. The naming convention is not required, but is a good practice to follow when using Stream.

Note: For this case, we are interested in examining the lift and drag coefficients on the plane. These values are computed from the forces on the plane, which are written out to specific files in the output/ directory. They are written out at the frequency specified by the print_freq variable in run control file.

Turbulent Results

The extract utility is used to view the results of a simulation. For this case we view the simulation at the 4,000th time step. In the run directory, use the following command: extract -vtk case 4000 v P k omega pResidualTT . This will generate a solution directory with files in the VTK format which can be opened by Paraview. The case name is case, the time step is 4000, and the variables to extract are the velocity(v), pressure(P), turbulent kinetic energy(k), turbulent specific dissipation rate(omega), and pressure residual turn-over time(pResidualTT).

Shown below are select plots of the solution.

Contour plot of the pressure field on the surface of the plane.

Stream tracer plot over the wing of the plane.

Vorticity magnitude contours at 0.1 meter spacing down the body of the plane.

The lift and the drag coefficients were computed from the x-component and z-component boundary force data in the flux_plane.dat file which is located in the output/ directory.

Note: The lift, by convention, is perpendicular to the flow direction and the drag is parallel to the flow direction. The transformation required to compute the lift and drag from the x-component and z-component force data is shown below. The angle of attack is represented by \(\alpha\).

\[ \begin{align}\begin{aligned}F_L = F_z cos(\alpha) - F_x sin(\alpha)\\F_D = F_z sin(\alpha) + F_x cos(\alpha)\end{aligned}\end{align} \]

The lift (\(F_L\)) and drag (\(F_D\)) forces can be non-dimensionalized to obtain the lift and drag coefficients (\(C_L\) and \(C_D\)) using the following formulas:

\[ \begin{align}\begin{aligned}C_L = \frac{F_L}{0.5 \rho_{\infty} V_{\infty}^2 S}\\C_D = \frac{F_D}{0.5 \rho_{\infty} V_{\infty}^2 S}\end{aligned}\end{align} \]

where \(S\) is the reference area of the plane (0.1091 m²), \(V_{\infty}\) is the free stream velocity (20 m/s), and \(\rho_{\infty}\) is the density (1.225 \(kg/m^3\)).

To estimate the timestep at which the simulation could be stopped and the coefficients extracted, the Stream boundary flux plotting utility was used. This script should be run in the simulation directory as it searches for an output/ directory. The name of the boundary (plane) should be passed to the utility, as shown in the command below.

Stream/bin/plot_bc_fluxes plane

The plot_bc_fluxes script generates an output directory with image plots of the fluxes versus iteration number. For this case we ran the simulation further than the 4000 timesteps in the run control file in order to see when the forces would stabilize. The x-component of the force on the plane stabilized at around 4000 timesteps as shown in the boundary flux plot below.

X-component force on the plane surface.

The lift and drag coefficients were computed from the forces on the plane at the final time step of the simulation. The results below show a suite of simulations (including this 10 degree angle-of-attack case) that cover a range of angles of attack between -5 and 20 degrees. The experimental data from the SSAM paper is included in the plots for comparison.

Comparison of lift and drag coefficients from Stream simulations and experimental data from the SSAM paper.