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. Currently, capability exists to simulate a gas with number of species (nsp) ≤ 5.
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 (see input.config in source code)
- 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 - Turbulent eddy viscosity, [Pa-s]
- pressure - Specification of static pressure over a surface, [Pa]
- Used to compute mole fractions of each species of the gas with Dalton's Law of partial pressures in conjunction with reference mole fractions specified in solver.inp
- 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 (related to initial vibrational temperature, see input.config in source code), [K]
- initial scalar_1 - Initial value of turbulent eddy viscosity, [Pa-s]
- initial pressure - Static pressure of the gas, [Pa]
- For multi-species flows, this value is used in combination with solver.inp values to compute the mole fraction of each species
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.
Example of solver.inp file:
#SOLUTION CONTROL
#{
Equation of State: Compressible
Number of Timesteps: 1000
Time Step Size: 1e-8
Turbulence Model: No-Model # RANS-SA
# Limit instructions : switch min max -- change switch from zero to activate
Limit Temperature : 0 0 0 # also limits vibrational temperature
Limit u1 : 0 0 0
Limit u2 : 0 0 0
Limit u3 : 0 -1 1
Limit rho1 : 0 1e-20 0
Limit rho2 : 0 1e-20 0
Limit rho3 : 0 1e-20 0
Limit rho4 : 0 1e-20 0
Limit rho5 : 0 1e-20 0
Limit Scalar 1 : 0 0 0
#}
#OUTPUT CONTROL
#{
Number of Timesteps between Restarts: 100
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
#{
Viscous Control: Viscous #None
Shear Law: Wilke's Mixing Rule # ishear=1 => matflag(2,n)
Bulk Viscosity Law: Constant Bulk Viscosity # ibulk=0 => matflag(3,n)
Conductivity Law: Wilke's Mixing Rule # icond=1 => matflag(4,n)
#}
#REACTING FLOW
#{
Number of species: 1 # nsp
Species IDs: 1 2 3 4 5 # ispecIDs
# IDs numbered in order: N2=1,O2=2,NO=3,N=4,O=5
Inflow Concentrations: 0.78997 0.21 1e-5 1e-5 1e-5 # concinf
Allow reactions: False # ichem = 1 if True. right now can only be used for nsp=5
Temperature threshold: 2000 # Tth--below which, reactions ignored. 2000 K is from Chalot 1990
Equilibrium Tolerance: 1e-5 # chemtol (max species production rate for reactions to equilibrium in IC's/BC's)
Two Temperature coefficient: 0.5 # qta (Tvib**qta*T**(1-qta))
Vibrational Temperature BC: 1 # TvibBC. must be a positive number. at any BC with T set:
# if greater than 5, then Tvib = TvibBC
# else, then Tvib = TvibBC*T
Vibrational Temperature IC: -1 # TvibIC (if negative, Tvib = T initially. if positive, value here is used)
Restart from primitive file: 0 # 1 if restarting from file from primitive code
#}
#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
#{
Shock Sensor : False # i_ss_on = 1 if True
Shock Sensor Value Full Off : 0.05 # % change of temperature across element above which we ramp mu until
Shock Sensor Value Full On : 0.2 # percent change after which we hold at scale factor * mu
Shock Sensor Scale Factor : 100.0 # scale factor on mu or other sensor
Wall Distance to Shield Shock Sensor : -1 # The above won't be applied within this wall distance (-1 ignores this condition)
#}
#STEP SEQUENCE
#{
Step Construction : 0 1 0 1 # laminar
# Step Construction : 0 1 0 1 10 11 10 11 # turbulent
#}
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 evisc and dwal fields would require turbulence modeling to be turned on for the example below to work. This is controlled through the solver.inp options.
Example of flow.pht file
<?xml version="1.0" ?>
<PhastaMetaFile number_of_pieces="40">
<GeometryFileNamePattern pattern="ent_from_NAS/geombc.dat.%d"
has_piece_entry="1"
has_time_entry="0"/>
<FieldFileNamePattern pattern="ent_from_NAS/restart.%d.%d"
has_piece_entry="1"
has_time_entry="1"/>
<TimeSteps number_of_steps="9"
auto_generate_indices="1"
start_index="65800"
increment_index_by="1000"
start_value="0"
increment_value_by="1">
</TimeSteps>
<Fields number_of_fields="10">
<Field phasta_field_tag="solution"
paraview_field_tag="rho_1"
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_2"
start_index_in_phasta_array="1"
number_of_components="1"
data_dependency="0"
data_type="double"/>
<Field phasta_field_tag="solution"
paraview_field_tag="rho_3"
start_index_in_phasta_array="2"
number_of_components="1"
data_dependency="0"
data_type="double"/>
<Field phasta_field_tag="solution"
paraview_field_tag="rho_4"
start_index_in_phasta_array="3"
number_of_components="1"
data_dependency="0"
data_type="double"/>
<Field phasta_field_tag="solution"
paraview_field_tag="rho_5"
start_index_in_phasta_array="4"
number_of_components="1"
data_dependency="0"
data_type="double"/>
<Field phasta_field_tag="solution"
paraview_field_tag="velocity"
start_index_in_phasta_array="5"
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="8"
number_of_components="1"
data_dependency="0"
data_type="double"/>
<Field phasta_field_tag="solution"
paraview_field_tag="temp"
start_index_in_phasta_array="9"
number_of_components="1"
data_dependency="0"
data_type="double"/>
<Field phasta_field_tag="solution"
paraview_field_tag="evisc"
start_index_in_phasta_array="10"
number_of_components="1"
data_dependency="0"
data_type="double"/>
<Field phasta_field_tag="dwal"
paraview_field_tag="dwal"
start_index_in_phasta_array="0"
number_of_components="1"
data_dependency="0"
data_type="double"/>
</Fields>
</PhastaMetaFile>
