TCNEQ Version
From PHASTA Wiki
Contents
Background
The following information relates to the use of the thermochemical nonequilibrium (TCNEQ) version of PHASTA written in terms of entropy variables. The reader is referred to the following for additional information.
- F. Chalot, T.J.R. Hughes, and F. Shakib, "Symmetrization of Conservation Laws with Entropy for High-Temperature Hypersonic Computations," Computing Systems in Engineering, 1(2-4):495–521, 1990.
- J. Pointer, "Influence of Interpolation Variables and Discontinuity Capturing Operators on Inviscid Hypersonic Flow Simulations Using a Stabilized Continuous Galerkin Solver," Ph.D. dissertation, University of Colorado, Boulder, CO, 2022.
Pre-Processing
In this section, details of the meshing and model attributes are provided. For simulation cases where the gas is modeled as a single species, the scalar quantities for boundary and initial conditions are not required. Currently, capability exists to simulate a gas with number of species (nsp) ≤ 5. When 2 ≤ nsp ≤ 5, the scalar quantities are used to govern the composition of the gas.
Meshing
Within the Simmodeler utility, the mesh can either be created or loaded from an existing .cas file. Below are steps for loading a mesh from a .cas file:
- Launch Simmodeler (for this example, SimModeler7.0-190604 is used)
- File > Import Discrete Data > (select .cas file to import) > (keep defaults and click OK) > (select YES to keep volume mesh)
- Save .sms and .smd files
- Attributes can now be assigned to the model as normal
Boundary Conditions
Below are the recognized boundary conditions that can be applied for the current version:
- comp1/comp2/comp3 - Specification of one/two/three components of velocity, [m/s]
- temperature - Specification of translational-rotational temperature, [K]. By default, vibrational temperature is held in equilibrium with this value and nonequilibrium is controlled through simulation inputs.
- surfID - When value is set to 702, the boundary is treated as a slip wall. If using this option, include a boundary layer mesh along the surface to ensure the wall normal direction is accurately computed.
- scalar_1 - Mole fraction of species 2 of the gas
- scalar_2 - Mole fraction of species 3 of the gas
- scalar_3 - Mole fraction of species 4 of the gas
- scalar_4 - Mole fraction of species 5 of the gas
- pressure - Specification of static pressure over a surface, [Pa]
- Used to compute mole fraction of species 1 of the gas with Dalton's Law of partial pressures and subtracting the summation of the other mole fractions from a value of 1
- heat flux - set to zero for adiabatic wall boundary condition
Initial Conditions
Below are the required initial conditions for the current version:
- initial velocity - Components and magnitude of flow velocity, [m/s]
- If a supersonic outlet condition is used, set such that flow is initialized Mach > 1
- initial temperature - Value used to set translational-rotational temperature, [K]
- initial scalar_1 - Initial value of species 2 mole fraction
- initial scalar_2 - Initial value of species 3 mole fraction
- initial scalar_3 - Initial value of species 4 mole fraction
- initial scalar_4 - Initial value of species 5 mole fraction
- initial pressure - Static pressure of the gas, [Pa]
- For multi-species flows, this value is used in combination with the initial scalar values to compute the mole fraction of species 1
Simulation Inputs
Below is an example of the input script for the current version of the code. Capability is included for handling multi-species flows up to number of species nsp equal to 5. This feature has been tested and shown to work for five species mass conservation equations at one time. However, when chemical species production governed by the finite-rate chemistry module is active, the solution becomes unstable. Additional work on this feature is needed.
Example of solver.inp file:
# PHASTA_HYP Version 1 Input File
#SOLUTION CONTROL
#{
Equation of State: Compressible
Number of Timesteps: 100
Time Step Size: 3.2339656949004e-08
Limit Density: 0 0.01 0.1 # solution limiting on variables [switch, min, max]
Limit u1: 0 0. 2.8e3
Limit u2: 0 0 0
Limit u3: 0 0 0
Limit Temperature: 0 230 3500 # also limits vibrational temperature
#}
#OUTPUT CONTROL
#{
Number of Timesteps between Restarts: 250
Print Error Indicators: True # shock error stored in column 6, DC factor \nu stored in column 10
Error Indicator Threshold: 0.01 # err > thresh*err_max is flagged as 1 (i.e. identified for refinement)
# --> smaller values = narrower flagged region along shock
Number of Error Smoothing Iterations: 0 # ierrsmooth
Load and set 3D IC: False # load the flowfield from a file as the initial condition
Position Tolerance on IC Load: 1e-7 # sets the tolerance for matching node locations while loading the initial condition
#}
#MATERIAL CONTROL
#{
Viscosity: 1.95508431704028e-5 # dynamic viscosity (only used if nsp.eq.99 for air mixture)
Thermal Conductivity: 26.6843390135759e-3 # only used if nsp.eq.99 for air mixture
Viscous Control: None #Viscous #None
#}
#REACTING FLOW
#{
Number of species: 1 # nsp
#.......currently only allowing 1<=nsp<=5.
#.........and make sure scalars passed from simmodeler are in correct slots
#
Species IDs: 1 # specIDs, length of array must equal
# nsp (IDs numbered in order:
# N2,O2,NO,N,O
# ID:99 for Air molecule
# e.g. Species IDs: 1 3 <<< would give N2 and NO
Ref Entropy Conditions: 1e3 230 230 0 #[P0,T0,T0vib,S0]
Ref Entropy Mole Frac: 1 0 0 0 0 # composition of gas used as reference entropy condition
# >> must sum to zero
Allow reactions: False # chemical reactions, ichem = 1 if True
Chemical heat release: False # chemical heat release, iqtot = 1 if True
Limit on reaction step: 0.00001 # rlim (limits change in species cs per step)
Tolerance to global time: 0.01 # ttol, chem solver is advanced in time until time diff < ttol*dt_global
Temperature threshold: 500 # Tth (below which, reactions ignored)
Reaction solver MIN steps: 5 # nstepmin, minimum number of time steps
Reaction solver MAX steps: 100 # nstepmax, maximum number of time steps
Two Temperature coefficient: 0.5 # qta (Tvib**qta*T**(1-qta))
Exclude vib energy: True # ivib0 = 1 if True
Exclude vib source: True # ivibS0 = 1 if True
Tvib BC Ratio: 1.0 # at any BC with T set, Tvib = tvibBC * T
Vibrational Temperature IC: -1 # set negative value to force Tvib = T, otherwise positive value set as IC value
#}
#LINEAR SOLVER
#
Solver Type: GMRES sparse
Number of GMRES Sweeps per Solve: 1 # replaces nGMRES
Minimum Number of Iterations per Nonlinear Iteration: 10 # minIters
Number of Krylov Vectors per GMRES Sweep: 100 # replaces Kspace
Tolerance on Momentum Equations: 0.01 # epstol(1), affects etol for Hessenberg problem
#}
#DISCRETIZATION CONTROL
#{
Weak Form: SUPG # alternate is Galerkin only for compressible
Time Integration Rule: First Order # 1st Order sets rinf(1) -1
Tau Matrix: Matrix-Ent-Adv
Include Viscous Correction in Stabilization: False # if p=1 idiff=1
# if p=2 idiff=2
Tau Time Constant: 1.0
Tau C Scale Factor: 1.0 # taucfct best value depends
Number of Elements Per Block: 64 #ibksiz
#}
#DISCONTINUITY CAPTURING
#{
Discontinuity Capturing: DC-quadratic # Current Options: DC-mallet, DC-minimum, DC-quadratic, DC-yzbeta
Multiplier for DC factor: 1 # scales DC variable in e3DC
Discontinuity Capturing Scheme: 1 # 0: discontinuous, 1: continuous (L2 projection)
Include Source Term in DC: 0 # 1: sets idcSRC to 1
Write DCqpt: 0 # if > 0, writes out data at a quadrature point iDCqpt
#----Parameters for YZBeta DC operator ----
Beta Value: 1 # 1: smoother , 2: sharper, 12: compromise between 1 and 2
YZB Farfield Conditions: 1e5 2119 10 10 300 # [Pressure, X-Vel, Y-Vel, Z-Vel, Temperature]
YZB Farfield Mole Frac: 1 0 0 0 0 # mole fractions at reference condition
# [xN2,xO2,xNO,xN,xO]
# must sum to 1, must be length 5
Include Umod Term: 1 # 0: no, 1: yes
Mach Adjustment Bm Value: 1 # 0: off,1: smoother shock, 2: sharper shock
Mach Adjustment Bj Value: 6 # 0: off,1: smoother shock, 2: sharper shock
Include Time Term in Z: 1 # 0: no, 1: yes
#------------------------------------------
#}
#STEP SEQUENCE
#{
Step Construction : 0 1 0 1
#}
Post-Processing
An example of the flow.pht file is provided to demonstrate the ordering of the variables that can be viewed in Paraview. Note that both the errors and DCqpt fields would need to be output to the restart file for the example below to work. Saving these two fields is controlled through the solver.inp options.
Example of flow.pht file
<?xml version="1.0" ?>
<PhastaMetaFile number_of_pieces="24">
<GeometryFileNamePattern pattern="24-procs_case/geombc.dat.%d"
has_piece_entry="1"
has_time_entry="0"/>
<FieldFileNamePattern pattern="24-procs_case/restart.%d.%d"
has_piece_entry="1"
has_time_entry="1"/>
<TimeSteps number_of_steps="1"
auto_generate_indices="1"
start_index="2010"
increment_index_by="50"
start_value="0"
increment_value_by="1">
</TimeSteps>
<Fields number_of_fields="9">
<Field phasta_field_tag="solution"
paraview_field_tag="rho_N2"
start_index_in_phasta_array="0"
number_of_components="1"
data_dependency="0"
data_type="double"/>
<Field phasta_field_tag="solution"
paraview_field_tag="rho_O2"
start_index_in_phasta_array="1"
number_of_components="1"
data_dependency="0"
data_type="double"/>
<Field phasta_field_tag="solution"
paraview_field_tag="velocity"
start_index_in_phasta_array="2"
number_of_components="3"
data_dependency="0"
data_type="double"/>
<Field phasta_field_tag="solution"
paraview_field_tag="temp_vib"
start_index_in_phasta_array="5"
number_of_components="1"
data_dependency="0"
data_type="double"/>
<Field phasta_field_tag="solution"
paraview_field_tag="temperature"
start_index_in_phasta_array="6"
number_of_components="1"
data_dependency="0"
data_type="double"/>
<Field phasta_field_tag="errors"
paraview_field_tag="nu"
start_index_in_phasta_array="9"
number_of_components="1"
data_dependency="0"
data_type="double"/>
<Field phasta_field_tag="DCqpt"
paraview_field_tag="DCqpt1"
start_index_in_phasta_array="0"
number_of_components="1"
data_dependency="1"
data_type="double"/>
</Fields>
</PhastaMetaFile>
