TCNEQ Version

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:

1. Launch Simmodeler (for this example, SimModeler7.0-190604 is used)
2. File > Import Discrete Data > (select .cas file to import) > (keep defaults and click OK) > (select YES to keep volume mesh)
3. Save .sms and .smd files
4. 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>