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. 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:

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 (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>