Control Dictionary

The zCFD control file is a python file that is executed at runtime. This provides a flexible way of specifying key parameters, enabling user defined functions and leveraging the rich array of python libraries available.

A python dictionary called ‘parameters’ provides the main interface to controlling zCFD

parameters = {....}

zCFD will validate the parameters dictionary specified at run time against a set of expected values for each key. In general, values that should be integers and floats will be coerced to the correct type. Values which should be strings, booleans, lists or dictionaries will cause an error if they are the incorrect type. The time marching scheme and solver/equation set chosen defines which dictionary keys are valid. A simple script is provided to allow the user to check the validity of their input parameters before submitting their job. See Input Validation.

Reference Settings

Define units

# Optional (default 'SI'):
# units for dimensional quantities
'units' : 'SI',

Scale mesh

# Optional (default [1.0,1.0,1.0]):
# Scale vector
'scale' : [1.0,1.0,1.0],
# Optional
# Scale function
'scale' : scale_func,

Example scale function

# Scale function
def scale_func(point):
    # Scale z coordinate by 100
    point[2] = point[2] * 100.0
    return {'coord' : point}

Reference state

# reference state
'reference' : 'IC_1',

See Initial Conditions

Partitioner

The ‘metis’ partitioner is used to split the computational mesh up so that the job can be run in parallel. Other partitioners will be added in future code releases

# Optional (default 'metis'):
# partitioner type
'partitioner' : 'metis',

Safe Mode

Safe mode turns on extra checks and provides diagnosis in the case of solver stability issues

# Optional (default False):
# Safe mode
'safe' : False,

Initialisation

Initial conditions are used to set the flow variable values in all cells at the start of a simulation. The initial conditions default to the reference state (above) if not specified

# Initial conditions
'initial' : 'IC_2',

See Initial Conditions. A function which defines the initial fields may be supplied

# Initial conditions
'initial' : {
              # Name of initial condition
              'name' : 'IC_1',
              # Optional: User defined function
              'func' : my_initialisation,
            },

Example User Defined Initialisation Function

def my_initialisation(**kwargs):

    # Dimensional primitive variables
    pressure = kwargs['pressure']
    temperature = kwargs['temperature']
    velocity = kwargs['velocity']
    wall_distance = kwargs['wall_distance']
    location = kwargs['location']

    if location[0] > 10.0:
      velocity[0] *= 2.0

    # Return a dictionary with user defined quantity.
    # Same schema as above
    return { 'velocity' : velocity }

To restart from a previous solution there are two options. A solve can be restarted on the same mesh using the options:

# Restart from previous solution
'restart' : True,
# Optional: Restart from a different case
'restart casename' : 'my_case',

Alternatively, a solution from another mesh can be mapped on to a new mesh using the options:

# Restart from previous solution on a different mesh
"interpolate restart": True,
# Required: Restart from a different case
"restart casename": "different_case",
# Required: Restart from a different mesh
"restart meshname": "different_mesh",
# Optional: Number of Laplace smoothing cycles to smooth initial mapped solution (default 0)
"solution smooth cycles": 0,

It is possible to restart the SA-neg solver from results generated using the Menter SST solver and vice versa, however by default the turbulence field(s) are reset to zero. An approximate initial solution for the SA-neg model can be generated from Menter-SST results by adding the key

# Generate approximate initial field for SA-neg
'approximate sa from sst results' : True,

Low Mach number preconditioner settings (Optional)

The preconditioner is a mathematical technique for improving the speed of the solver when the fluid flow is slow compared to the speed of sound. Use of the preconditioner does not alter the final converged solution produced by the solver. The preconditioner factor is used to improve the robustness of the solver in regions of very low speed. The preconditioner is controlled in the 'equations' dictionary (see Equations) by setting the 'precondition' key, for example when using the euler solver

'euler' : {
          # Optional (default False):
          # Use low speed mach preconditioner
          'precondition' : True,
        },

# Optional, Advanced: Preconditioner factor
'preconditioner' : {
                    'minimum mach number' : 0.5,
                   },

Time Marcher

zCFD is capable of performing both steady-state and time accurate simulations. In the first case the solution is marched in pseudo time using the dual time-stepping algorithm towards the steady-state solution. For time accurate simulations the solution can be either globally time-stepped or use dual time-stepping. The 'time marching' dictionary controls the selection of these schemes.

'time marching' : {....},

Time accurate (unsteady) simulation control

'time marching' : {
                   'unsteady' : {
                                 # Total time in seconds
                                 'total time' : 1.0,
                                 # Time step in seconds
                                 'time step' : 1.0,
                                 # Required only when dual time-stepping:
                                 # Dual time step time accuracy (options: 'first' or 'second' order)
                                 'order' : 'second',
                                 # Required only when dual time-stepping:
                                 # Number of pseudo time (steady) cycles to run before starting dual timestep time accurate simulation
                                 'start' : 3000,
                                 },
                   },

Solver scheme

There are three different time marching schemes available in zCFD: Explicit Euler, various orders of Runge Kutta (like Euler…) and the implicit Euler scheme. These can be run with dual time-stepping for both steady and unsteady simulations. Explicit schemes with no preconditioning can be used with global time-stepping for unsteady simulations. The choice of scheme and time stepping algorithm affects which dictionary keys are valid.

Runge Kutta

'time marching' : {
                  'scheme' : {
                              # Either 'euler' or 'runge kutta' or 'implicit euler'
                              'name' : 'runge kutta',
                              # Number of RK stages 'euler' or 1, 4, 5, 'rk third order tvd'
                              'stage': 5,
                              # Time-stepping scheme: 'local timestepping' for steady state or dual time step time accurate
                              # 'global timestepping' for time accurate
                              'kind' : 'local timestepping',
                              # With AMGX, the linearsolver can be run in single or double precision modes
                              'double precision linear solve': True,
                              },
                   },

When global time-stepping is specified the keys ['unsteady']['order'] and ['unsteady']['start'] are not valid. The CFL keys in in the CFL section below are also not valid as the time step is explicitly specified. When the ['unsteady']['total time'] and ['unsteady']['time step'] keys are set to be equal values and ['time marching']['scheme']['kind'] is set to 'local timestepping' the solver marches in pseudo time towards a steady state solution. If implicit euler is selected then the 'multigrid' keys must be removed. For large CFL numbers, it is recommended the Roe inviscid flux scheme is used.

CFL

The Courant-Friedrichs-Lewy (CFL) number controls the local pseudo time-step that the solver uses to reach a converged solution. The larger the CFL number, the faster the solver will run but the less stable it will be. The default values should be appropriate for most cases, but for lower-quality meshes or very complex geometries it may be necessary to use lower values (e.g. CFL of 1.0). In such cases it may also be helpful to turn off Multigrid (above) by setting the maximum number of meshes to 0. The CFL dictionary keys are only valid when using dual time-stepping.

'time marching' : {
                  # Default max CFL number for all equations
                  'cfl': 2.5
                  # Optional (default 'cfl' value):
                  # Override max CFL number for transported quantities
                  'cfl transport' : 1.5,
                  # Optional (default 'cfl' value):
                  # Override max CFL number for coarse meshes
                  'cfl coarse' : 2.0,
                  # Control of CFL for polynomial multigrid (highest and lowest order)
                  'multipolycfl' : [2.0,2.0],
                  # Optional: Ramp cfl
                  'cfl ramp factor': {
                                     # cfl = cfl_ini * growth ^ cycle
                                     'growth': 1.05,
                                     # min cfl
                                     'initial': 0.1,
                                     },
                  },

Cycles

For steady-state simulations, the number of pseudo time cycles is the same as the number of steps that the solver should use to reach a converged solution. Note that the solver uses local pseudo time-stepping (the time-step varies according to local conditions) so any intermediate solution is not necessarily time-accurate.

For unsteady (time-accurate) simulations using ‘dual time-stepping’ to advance the solution in time the number of pseudo time cycles determines the number of inner iterations that the solver uses to converge each real time step. When global time-stepping is specified the cycles key is not valid.

'time marching' : {
                  # Required only for dual time-stepping:
                  # Number of pseudo time cycles
                  'cycles' : 5000,
},

Multigrid

The mesh is automatically coarsened by merging cells on successive layers in order to accelerate solver convergence. This does not alter the accuracy of the solution on the finest (original) mesh. Advanced users can control the number of geometric multigrid levels and the ‘prolongation’ of quantities from coarse to fine meshes.

'time marching' : {
                # Maximum number of meshes (including fine mesh)
                'multigrid' : 10,
                # Optional:
                # Number of multigrid cycles before solving on fine mesh only
                'multigrid cycles' : 5000,
                # Optional (default 0.75), Advanced:
                # Prolongation factor
                'prolong factor' : 0.75,
                # Optional (default 0.3), Advanced:
                # Prolongation factor for transported quantities
                'prolong transport factor' : 0.3,
                },

Polynomial Multigrid

The polynomial basis on which the solution is computed can be successively coarsened, thus allowing the solution to be evolved quicker due to weakened stability restrictions at lower polynomial orders. This does not alter the accuracy of the solution on the highest polynomial order.

'time marching' : {
                   # Optional (default False):
                   # Switch on polynomial multigrid
                   'multipoly' : True,
                   },

Equations

The equations key selects the equation set to be solved as well as the finite volume solver or the high order strong form Discontinuous Galerkin/Flux Reconstruction solver. For each equation type the valid keys are listed below.

# Governing equations to be used
# one of: RANS, euler, viscous, LES, DGRANS, DGviscous, DGLES, DGCAA
'equations' : 'RANS',

Compressible Euler flow is inviscid (no viscosity and hence no turbulence). The compressible Euler equations are appropriate when modelling flows where momentum significantly dominates viscosity - for example at very high speed. The computational mesh used for Euler flow does not need to resolve the flow detail in the boundary layer and hence will generally have far fewer cells than the corresponding viscous mesh would have.

'euler' : {
            # Spatial accuracy (options: first, second)
            'order' : 'second',
            # Optional (default 0):
            # Run first order spatial accuracy for the first x cycles
            'first order cycles': 100,
            # Optional (default 'vanalbada'):
            # MUSCL limiter (options: vanalbada)
            'limiter' : 'vanalbada',
            # Optional (default off)
            # Cycle on which to freeze the limiter values
            'freeze limiter cycle': 500,
            # Optional (default False):
            # Use low speed mach preconditioner
            'precondition' : True,
            # Optional (default False):
            # Use linear gradients
            'linear gradients' : False,
            # Optional (default 'HLLC'):
            # Scheme for inviscid flux: HLLC, Rusanov or Roe
            'Inviscid Flux Scheme': 'HLLC',
          },

The viscous (laminar) equations model flow that is viscous but not turbulent. The Reynolds number of a flow regime determines whether or not the flow will be turbulent. The computational mesh for a viscous flow does have to resolve the boundary layer, but the solver will run faster as fewer equations are being included.

'viscous' : {
              # Spatial accuracy (options: first, second)
              'order' : 'second',
              # Optional (default 0):
              # Run first order spatial accuracy for the first x cycles
              'first order cycles': 100,
              # Optional (default 'vanalbada'):
              # MUSCL limiter (options: vanalbada)
              'limiter' : 'vanalbada',
              # Optional (default False):
              # Use low speed mach preconditioner
              'precondition' : True,
              # Optional (default False):
              # Use linear gradients
              'linear gradients' : False,
              # Optional (default 'HLLC'):
              # Scheme for inviscid flux: HLLC or Rusanov
              'Inviscid Flux Scheme': 'HLLC',
            },

The fully turbulent (Reynolds Averaged Navier-Stokes Equations)

'RANS' : {
            # Spatial accuracy (options: first, second, euler_second)
            'order' : 'second',
            # Optional (default 0):
            # Run first order spatial accuracy for the first x cycles
            'first order cycles': 100,
            # Optional (default 'vanalbada'):
            # MUSCL limiter (options: vanalbada)
            'limiter' : 'vanalbada',
            # Optional (default off)
            # Cycle on which to freeze the limiter values
            'freeze limiter cycle': 500,
            # Optional (default False):
            # Use low speed mach preconditioner
            'precondition' : True,
            # Optional (default False):
            # Use linear gradients
            'linear gradients' : False,
            # Optional (default False):
            # Use least squares gradients
            'leastsq gradients' : False,
            # Optional (default 'HLLC'):
            # Scheme for inviscid flux: HLLC, Rusanov, Roe or Roe low diss
            'Inviscid Flux Scheme': 'HLLC',
            # Optional (default 'NONE')
            # Roe low dissipation sensor type: 'FD' (wall distance based), 'NTS' (vorticity based), 'DES' (DES blending function based) or 'NONE'
            'roe low dissipation sensor': 'DES',
            # Optional (default 0.05)
            # Minimum value of the Roe low dissipation sensor
            'roe dissipation sensor minimum': 0.05,
            # Turbulence
            'turbulence' : {
                              # turbulence model (options: 'sst', 'sa-neg', 'sst-transition')
                              'model' : 'sst',
                              # Optional (default 'none'):
                              # LES model (options 'none', 'DES', 'DDES', 'SAS')
                              'les' : 'none',
                              # Optional (default 0.09):
                              # betastar turbulence closure constant
                              'betastar' : 0.09,
                              # Optional (default True):
                              # turn off mu_t limiter
                              'limit mut' : False,
                              # Optional (default CDES=0.65, CDES_kw=0.78, CDES_keps=0.61):
                              # DES constants
                              'CDES': 0.65
                              'CDES_kw': 0.78,
                              'CDES_keps': 0.61,
                              # Optional (default Cd1=20.0, Cd2=3, Cw= 0.15):
                              # (I)DDES constants
                              'Cd1': 20.0,
                              'Cd2': 3,
                              'Cw': 0.15,
                              # Optional when using 'sst' (default 0):
                              # Menter SST production term: 0=SST-V, 1=incompressible, 2=SST
                              'production': 0,
                              # Optional when using 'sa-neg' (default False):
                              # Use sa-neg rotation correction
                              'rotation correction': True,
                              # Optional when using sa-neg (default 0.0):
                              # Limit the sa-neg gradient based on the maximum value between neighbouring cells.
                              # k = 0.0 corresponds to no limiting, k = 1.0 maximum limiting.
                              'limit gradient k': 0.5
                              # Option when using 'sst-transition'
                              # Transition model constants
                              'ca1': 2.0,
                              'ca2': 0.06,
                              'ce1': 1.0,
                              'ce2': 50.0,
                              'cthetat': 0.03,
                              'sigmagamma': 1.0,
                              'sigmathetat': 2.0,
                              # Option when using 'sst-transition'
                              # Turn transition separation correction model on/off
                              'separation correction': True,
                              # Option when using 'sst' or 'sa-neg'
                              # Use the non-linear Quadratic Constitutive Relation
                              # Spalart, P. R., "Strategies for Turbulence Modelling and Simulation,"
                              # International Journal of Heat and Fluid Flow, Vol. 21, 2000, pp. 252-263,
                              # https://doi.org/10.1016/S0142-727X(00)00007-2
                              'qcr': True
                            },
           },

The filtered Large Eddy Simulation (LES) equations

'LES' : {
              # Spatial accuracy (options: first, second)
              'order' : 'second',
              # Optional (default 'vanalbada'):
              # MUSCL limiter (options: vanalbada)
              'limiter' : 'vanalbada',
              # Optional (default False):
              # Use low speed mach preconditioner
              'precondition' : True,
              # Optional (default False):
              # Use linear gradients
              'linear gradients' : False,
              # Optional (default False):
              # Use least squares gradients
              'leastsq gradients' : False,
              # Optional (default 'HLLC'):
              # Scheme for inviscid flux: HLLC, Rusanov, Roe or Roe low diss
              'Inviscid Flux Scheme': 'HLLC',
              # Optional (default 'NONE')
              # Roe low dissipation sensor type: 'FD' (wall distance based), 'NTS' (vorticity based), 'DES' (DES blending function based) or 'NONE'
              'roe low dissipation sensor': 'DES',
              # Optional (default 0.05)
              # Minimum value of the Roe low dissipation sensor
              'roe dissipation sensor minimum': 0.05,
              'turbulence' : {
                              # LES model (options 'none', 'WALE')
                              'les' : 'none',
                            },

            },

When using the high order strong form Discontinuous Galerkin/Flux Reconstruction solver the keys in the equations dictionary are the same for each equation set except that the 'linear gradients' and 'limiter' keys are not valid and additional keys are available in the equations dictionary. These additional keys are given below

'DG...' : {
               # Spatial polynomial order 0,1,2,3
               'order' : 2,
               # Optional (default 0.0)
               'c11 stability parameter': 0.0,
               # Optional (default 0.5)
               'LDG upwind parameter': 0.5,
               # Optional (default False)
               'Shock Sensing': False,
               'Shock Sensing k': 1.0,
               'Shock Sensing Viscosity Scale': 1.0,
               # Variable used for shock sensing, one of: 'density', 'temperature', 'mach', 'turbulence'
               'Shock Sensing Variable': 'density',
            },

Hence all the available keys for 'DGviscous' are a combination of the 'viscous' keys and the 'DG...' keys. The same logic applies to the 'DGRANS', and 'DGLES' equation sets.

'DGviscous' : {
               # Spatial polynomial order 0,1,2,3
               'order' : 2,
               # Optional (default False)
               # Use low speed mach preconditioner
               'precondition' : True,
               # Optional (default 'HLLC')
               # scheme for inviscid flux: HLLC or Rusanov
               'Inviscid Flux Scheme': 'HLLC',
               # Optional (default 0.0)
               'c11 stability parameter': 0.0,
               # Optional (default 0.5)
               'LDG upwind parameter': 0.5,
               # Optional (default False)
               'Shock Sensing': False,
               'Shock Sensing k': 1.0,
               'Shock Sensing Viscosity Scale': 1.0,
               # Variable used for shock sensing, one of: 'density', 'temperature', 'mach', 'turbulence'
               'Shock Sensing Variable': 'density',
               # Optional (default False)
               'inviscid flow': False,
               # Optional (default False)
               'freeze diffusion during rk stages': False
            },

DG Nodal locations

The location of the DG solution points must be specified in the parameters dictionary. The definitions first need to be imported by including this import statement at the start of the control dictionary

from zcfd.solvers.utils.DGNodalLocations import *

To use the default values

'Nodal Locations' : nodal_locations_default['Nodal Locations']

Or select from the available types

'Nodal Locations' : {
                      # Options line_evenly_spaced, line_gauss_lobatto or line_gauss_legendre_lobatto
                      'Line':  line_gauss_legendre_lobatto,
                      # Options tet_evenly_spaced, tet_shunn_ham
                      'Tetrahedron': tet_evenly_spaced,
                      # Options tri_evenly_spaced, tri_shunn_ham
                      'Tri' : tri_evenly_spaced,
                    },

Material Specification

The user can specify any fluid material properties by following the (default) scheme for ‘air’:

'material' : 'air',

Options

'air' : {
          'gamma' : 1.4,
          'gas constant' : 287.0,
          'Sutherlands const': 110.4,
          'Prandtl No' : 0.72,
          'Turbulent Prandtl No' : 0.9,
        },

Initial Conditions

The initial condition properties are defined using consecutively numbered blocks

'IC_1' : {....},
'IC_2' : {....},
'IC_3' : {....},

Each block can contain the following options

# Required: Static temperature in Kelvin
'temperature': 293.0,
# Required: Static pressure in Pascals
'pressure':101325.0,
# Fluid velocity
'V': {
        # Velocity vector
        'vector' : [1.0,0.0,0.0],
        # Optional: specifies velocity magnitude
        'Mach' : 0.20,
      },

Dynamic (shear, absolute or molecular) viscosity should be defined at the static temperature previously specified. This can be specified either as a dimensional quantity or by a Reynolds number and reference length

# Dynamic viscosity in dimensional units
'viscosity' : 1.83e-5,

or

# Reynolds number
'Reynolds No' : 5.0e6,
# Reference length
'Reference Length' : 1.0,

Turbulence intensity is defined as the ratio of velocity fluctuations \(u^{'}\) to the mean flow velocity. A turbulence intensity of 1% is considered low and greater than 10% is considered high.

# Turbulence intensity %
'turbulence intensity': 0.01,
# Turbulence intensity for sustainment %
'ambient turbulence intensity': 0.01,

The eddy viscosity ratio \((\mu_t/\mu)\) varies depending type of flow. For external flows this ratio varies from 0.1 to 1 (wind tunnel 1 to 10)

For internal flows there is greater dependence on Reynolds number as the largest eddies in the flow are limited by the characteristic lengths of the geometry (e.g. The height of the channel or diameter of the pipe). Typical values are:

Re 3000 5000 10,000 15,000 20,000 > 100,000
eddy 11.6 16.5 26.7 34.0 50.1 100
# Eddy viscosity ratio
'eddy viscosity ratio': 0.1,
# Eddy viscosity ratio for sustainment
'ambient eddy viscosity ratio': 0.1,

The user can also provide functions to specify a ‘wall-function’ - or the turbulence viscosity profile near a boundary. For example an atmospheric boundary layer (ABL) could be specified like this:

'profile' : {
             'ABL' : {
                       'roughness length' : 0.0003,
                       'friction velocity' : 0.4,
                       'surface layer height' : -1.0,
                       'Monin-Obukhov length' : -1.0,
                       # Non dimensional TKE/friction velocity**2
                       'TKE' : 0.928,
                       # ground level (optional if not set wall distance is used)
                       'z0'  : -0.75,
                      },
            },
'profile' : {
             'field' : 'inflow_field.vtp',
             # Localise field using wall distance rather that z coordinate
             'use wall distance' : True,
            },

Note

The conditions in the VTK file are specified by node based arrays with names ‘Pressure’, ‘Temperature’, ‘Velocity’, ‘TI’ and ‘EddyViscosity’. Note the field will override the conditions specified previously therefore the user can specify only the conditions that are different from default.

Certain conditions are specified relative to a reference set of conditions

'reference' : 'IC_1',
# total pressure/reference static pressure
# This can be a single value or VTK field file name with array name 'TotalPressureRatio' or function
'total pressure ratio' : 1.0,
# total temperature/reference static temperature
# This can be a single value or VTK field file name with array name 'TotalTemperatureRatio' or function
'total temperature ratio' : 1.0,
# Mach number
'mach' : 0.5,
# Direction vector or function
'vector' : [1.0,0.0,0.0],
'reference' : 'IC_1',
# static pressure/reference static pressure
# This can be a single value or VTK field file name with array name 'StaticPressureRatio' or function
'static pressure ratio' : 1.0,
'reference' : 'IC_1',
# Mass flow ratio
'mass flow ratio' : 1.0,
# Or mass flow rate
'mass flow rate' : 10.0
# Example pressure ratio function
def total_pressure_ratio(x, y, z):

    ratio = 1.0
    if x < 10.0:
        ratio = 0.9

    return ratio

Note

\(W^* =\frac{\dot{m}}{\rho V}\)

Driven Initial Conditions

Another way of prescribing a farfield initial condition is by defining a ‘driving function’, which on evaluation returns an initial condition dictionary. The driving function is re-evaluated at zCFD’s reporting frequency, meaning that the condition (e.g. inlet velocity) can be programmed to change as the simulation progresses. For a driving function lift_target_ic_func, the IC should be set up as follows in the parameters dictionary.

'IC_2' : lift_target_ic_func

At startup, the function is only provided with the key word arguments RealTimeStep=0, Cycle=0. However, the driving function is subsequently provided with key word arguments containing all the data in the ‘_report.csv’ file, evaluated at the previous timestep. For any driven initial condition, every value in the initial condition dictionary is also appended to the _report.csv file.

This means that the initial condition can be programmed to respond to the flow. An example is shown below, where the driven IC will continually adjust inlet angle of attack until a specified lift coefficient is acheived.

If the driven initial condition has been programmed to respond to changes in the flow, it is important to ensure that the flow has had enough time to propogate from the farfield through to the region of interest and settle before readjusting the driven initial condition. For simulations where the farfield is distant from the geometry, local and implicit timestepping will likely give the fastest propogation of adjusted farfield conditions through the domain.

def lift_target_ic_func(**kwargs):
   alpha_init = 0 # Initial angle of attack
   lift_target= 0.5 # Targeted lift coefficient
   update_period = 1000 # How often alpha should be updated
   flow_settling_period = 1500

   assumed_d_cL_d_alpha = 0.1 # Used to estimate alpha which would give required lift
   relaxation_factor = 1.15 # <1: under-relaxation, >1: over-relaxation


   if 'lift_target_alpha' in kwargs.keys(): # If timestep > 1
       alpha_current = kwargs['lift_target_alpha']
       F_xyz = [kwargs['wall_Fx'], kwargs['wall_Fy'], kwargs['wall_Fz']]
       F_LDS = zutil.rotate_vector(F_xyz, alpha_current, 0.0)
       lift_current = F_LDS[2] # Lift coefficient, having rotated to account for alpha
   else:
       lift_current = 0 # placeholder value

   if kwargs['Cycle'] < flow_settling_period:
       alpha_new = alpha_init
   else:
       if (kwargs['Cycle'] % update_period) == 0:
           d_cL_required = lift_target - lift_current
           d_alpha = relaxation_factor / assumed_d_cL_d_alpha *  d_cL_required
           alpha_new = alpha_current + d_alpha
           print("Alpha updated from {} to {}".format(alpha_current, alpha_new), flush = True)
       else:
           alpha_new = alpha_current

   dict_out = {
       "temperature": 300,
       "pressure": 101325.0,
       'v': {
           'mach' : 0.15,
           'vector' : zutil.vector_from_angle(alpha_new, 0.0)
           },
       "reynolds no": 6.0e6,
       "eddy viscosity ratio": 1.0,
       "monitor": { # Monitor entries can be floats or lists of floats
                    "prefix": "lift_target",
                    "alpha": alpha_new,   # Extra keys can be added to a driven IC dictionary,
                    "lift": lift_current,  # allowing monitoring of key parameters flow in the _report.csv
                  },
   }

   return dict_out

Note

Driven initial conditions cannot be used as reference conditions, to initialise the flow or within a restarted simulation.

Boundary Conditions

Boundary condition properties are defined using consecutively numbered blocks

'BC_1' : {....},
'BC_2' : {....},

zCFD will automatically detect zone types and numbers in a number of mesh formats, and assign appropriate boundary conditions. The type tags follow the Fluent convention (wall = 3, etc), and if present no further information is required. Alternatively, the mesh format may contain explicitly numbered zones (which can be determined by inspecting the mesh). In this case, the user can specify the list of zone numbers for each boundary condition ‘type’ and ‘kind’ (see below).

Wall

Wall boundaries in zCFD can be either slip walls, no-slip walls, or use automatic wall functions. No-slip walls are appropriate when the mesh is able to fully resolve the boundary layer (i.e \((y^{+} \leq 1)\)) otherwise automatic wall functions should be used. This behaviour is chosen by setting 'kind' in the wall specification below. Optionally it is possible to set wall velocity, roughness and temperature.

'BC_1' : {
          # Zone type tag
          'ref' : 3,
          # Optional: Specific zone boundary condition override
          'zone' : [0,1,2,3],
          # Required: Boundary condition type
          'type' : 'wall',
          # Type of wall, one of 'slip', 'noslip', 'wallfunction'
          'kind' : 'slip',
          # Optional: Roughness specification
          'roughness' : {
                        # Type of roughness specification (option: height or length)
                        'type' : 'height',
                        # Constant roughness
                        'scalar' : 0.001,
                        # Roughness field specified as a VTK file
                        'field' : 'roughness.vtp',
                        },
          # Optional: Wall velocity specification either 'linear' or 'rotating' dictionary supplied
          'V' : {
                'linear' : {
                           # Velocity vector
                           'vector' : [1.0,0.0,0.0],
                           # Optional: specifies velocity magnitude
                           'Mach' : 0.20,
                           },
                'rotating' : {
                             # Rotational velocity in rad/s
                             'omega' : 2.0,
                             # Rotation axis
                             'axis' : [1.0,0.0,0.0],
                             # Rotation origin
                             'origin' : [0.0,0.0,0.0],
                             },
                },
          # Optional: Wall temperature specification either a constant 'scalar' value or a 'profile' from file
          'temperature' : {
                          # Temperature in Kelvin
                          'scalar' : 280.0,
                          # Temperature field specified as a VTK file
                          'field' : 'temperate.vtp',
                          },
          },

Note

The roughness at each boundary face is set by finding the nearest point to the face centre on the supplied VTK file with the roughness value looked up in a node based scalar array called ‘Roughness’

Note

The temperature at each boundary face is set by finding the nearest point to the face centre on the supplied VTK file with the temperature value looked up in a node based scalar array called ‘Temperature’

Farfield

The farfield boundary condition can automatically determine whether the flow is locally an inflow or an outflow by solving a Riemann equation. The conditions on the farfield boundary are set by referring to an 'IC_' block. In addition to this an atmospheric boundary layer profile or a turbulence profile may be set.

'BC_2' : {
         # Zone type tag
         'ref' : 9,
         # Optional: Specific zone boundary condition override
         'zone' : [0,1,2,3],
         # Required: Boundary condition type
         'type' : 'farfield',
         # Required: Kind of farfield options (riemann, pressure or supersonic)
         'kind' : 'riemann',
         # Required: Farfield conditions
         'condition' : 'IC_1',
         # Optional: Atmospheric Boundary Layer
         'profile' : {
                     'ABL' : {
                             # Roughness length
                             'roughness length' : 0.05,
                             # Fiction velocity
                             'friction velocity' : 2.0,
                             # Surface Layer Heigh
                             'surface layer height' : 1000,
                             'Monin-Obukhov length' : 2.0,
                             # Turbulent kinetic energy
                             'TKE' : 1.0,
                             # Ground Level
                             'z0' : 0.0,
                              },
                      },
         # Optional: Specify a turbulence profile
         'turbulence' : {
                        'length scale' : 'filename.vtp',
                        'reynolds tensor' : 'filename.vtp',
                        },
         },

Inflow

zCFD can specify two kinds of inflow boundary condition: pressure or massflow, selected using the 'kind' key. When using the pressure inflow the boundary condition is specified by a total pressure ratio and a total temperature ratio that needs to be defined by the condition this refers to. When using the massflow condition a total temperature ratio and mass flow ratio or a mass flow rate is required. Note the mass flow ratio is defined relative to the condition set in the 'reference' key. By default the flow direction is normal to the boundary but either a single vector or a python function can be used to specify the inflow normal vector of each face on the surface.

'BC_3' : {
         # Zone type tag
         'ref' : 4,
         # Optional: Specific zone boundary condition override
         'zone' : [0,1,2,3],
         # Required: Boundary condition type
         'type' : 'inflow',
         # Required: Kind of inflow, either 'default' for pressure or 'massflow'
         'kind' : 'default',
         # Required: Reference inflow conditions
         'condition' : 'IC_2',
         },

Example Pressure Inflow

An example setup for a pressure inflow is given below:

'IC_1' : {
         'temperature' : 293.0,
         'pressure' : 101325.0,
         'V': {
              'vector' : [1.0,0.0,0.0],
              'Mach' : 0.20,
              },
         'Reynolds No' : 1.0e6,
         'Reference Length' : 1.0,
         'turbulence intensity' : 0.1,
         'eddy viscosity ratio' : 1.0,
         },
'IC_2' : {
         'reference' : 'IC_1',
         'total temperature ratio' : 1.0,
         'total pressure ratio' : 1.0,
         # Optional - specify mach number to set
         # static pressure at inflow from total pressure
         'mach' : 0.1,
         },
'BC_3' : {
         'ref' : 4,
         'zone' : [0,1,2,3],
         'type' : 'inflow',
         'kind' : 'default',
         'condition' : 'IC_2',
         },

Example Massflow Inflow

An example setup for a massflow inflow is given below:

def inflow_direction(x, y, z):
    # x, y, and z are the coordinates of the face centre
    cen = [1.0, 0.0, 0.0]
    y1 = y - cen[1]
    z1 = z - cen[2]
    rad = math.sqrt(z1*z1 + y1*y1)
    phi = math.atan2(y1,z1)
    angle = -0.24
    x2 = math.cos(angle)
    y2 = math.sin(phi) * math.sin(angle)
    z2 = math.cos(phi) * math.sin(angle)
    return {'v1': x2, 'v2': y2, 'v3': z2}

'IC_1' : {
         'temperature' : 293.0,
         'pressure' : 101325.0,
         'V': {
              'vector' : [1.0,0.0,0.0],
              'Mach' : 0.20,
              },
         'Reynolds No' : 1.0e6,
         'Reference Length' : 1.0,
         'turbulence intensity' : 0.1,
         'eddy viscosity ratio' : 1.0,
         },
'IC_2' : {
         'reference' : 'IC_1',
         'total temperature ratio' : 1.0,
         'mass flow rate' : 10.0,
         'vector' : inflow_direction
         },
'BC_3' : {
         'ref' : 4,
         'zone' : [0,1,2,3],
         'type' : 'inflow',
         'kind' : 'massflow',
         'condition' : 'IC_2',
         },

Outflow

Pressure or massflow outflow boundaries are specified in a similar manner to the Inflow boundary. They can be of :code`’kind’` 'pressure', 'radial pressure gradient' or 'massflow'. When either pressure boundary is set the 'IC_' conditions it refers to must set a static temperature ratio, when 'massflow' is specified a mass flow ratio or mass flow rate must be set. Note the mass flow ratio is defined relative to the condition set in the 'reference' key.

'BC_4' : {
         # Zone type tag
         'ref' : 5,
         # Optional: Specific zone boundary condition override
         'zone' : [0,1,2,3],
         # Boundary condition type
         'type' : 'outflow',
         # Kind of outflow 'pressure, 'massflow' or 'radial pressure gradient'
         'kind' : 'pressure',
         # Required only when using 'radial pressure gradient'
         # Radius at which the static pressure ratio is defined
         'reference radius' : 1.0,
         # Outflow conditions
         'condition' : 'IC_2',
         },

Symmetry

'BC_5' : {
         # Zone type tag
         'ref' : 7,
         # Optional: Specific zone boundary condition override
         'zone' : [0,1,2,3],
         # Required: Boundary condition type
         'type' : 'symmetry',
         },

Periodic

This boundary condition needs to be specified in pairs with opposing transformations. The transforms specified should map each zone onto each other.

'BC_6' : {
         # Required: Specific zone index
         'zone' : [1],
         # Required: Type of boundary condition
         'type' : 'periodic',
         # Required: Either 'rotated' or 'linear' periodicity
         'kind' : {
                  # Rotated periodic settings
                  'rotated' : {
                              'theta' : math.radians(120.0),
                              'axis' : [1.0,0.0,0.0],
                              'origin' : [-1.0,0.0,0.0],
                              },
                  # Linear periodic settings
                  'linear' : {
                             'vector' : [1.0,0.0,0.0],
                             },
                },

Fluid Structure Interaction

zCFD has the cability of simulating Fluid Structure Interaction (FSI) using both a modal model, as well as any generic full FEA program through the use of our generic coupling scheme.

General FSI

The general FSI capabilities of zCFD allow the solver to expose solution variables at run time, in order to allow for a user to couple the solver to any structural model. The behaviour of the coupling scheme can be modified within the open source python layer, by editing /bin/zcfd/utils/coupling/genericfsi.py. The default file contains 4 hooks definining points within the code when FSI actions can be performed- post intialisation, start of the real time cycle, post advance and post solve. In addition the file contains examples of all the functions available to the user for data communication, and how they should be executed in regards the rank of the CPU. In general any user written code for communication should be performed on the rank 0 CPU to avoid issues with parallelisation. The generic FSI scheme can be setup using the fsi key in the control dictionary.

'fsi': {
    # Required: The fsi driver to launch into, for generic FSI use the key 'generic'
    'type' : 'generic',
    # Required: The wall zones in the mesh which correspond to the wetted surface of the structural model
    'zone' : [1, 2, 3],
    # Optional: Zones in the mesh which remain stationary
    'fixed zones' : [4, 5, 6],
    # Required: Type of mesh transform to use: 'rbf multiscale', 'rbf' or 'idw'
    'transform type' : 'rbf multiscale',
    # Required (RBF only): Support radius of the RBFs
    'alpha' : 0.02,
    # Required (RBF only): Basis function for the rbf: 'c0', 'c2', 'c4', 'tps', 'gaussian' (default = 'c2')
    'basis' : 'c2',
    # Required (Multiscale RBF only): Fraction of surface points in the RBF base set
    'base point fraction': 0.1,
    # Required (IDW only): number of nearest neighbours to
    'n nearest' : 250,
    # Required (IDW only): power of the interpolation weight in the far field, at the boundary and the blending function between them
    'power' : [3.0, 5.0, 8.0],
    # Required (IDW only): distance over which to decay the deformation
    'deformation distance' : 100.0,
    # Required (IDW only): Stiffness of the blending function used to decay the deformation (0.0-1.0)
    'blending stiffness' : 0.5
    # Optional: extra user variables to pass to solver through validation script
    'user variables' : {'data path' : '/path/to/data'}
},

Fluid Zones

The fluidic zone properties are defined using consecutively numbered blocks like

'FZ_1' : {....},
'FZ_2' : {....},
'FZ_3' : {....},

For actuator disk zones

'FZ_1':{
        'type':'disc',
        'def':'T38-248.75.vtp',
        'thrust coefficient':0.84,
        'tip speed ratio':6.0,
        'centre':[-159.34009325,-2161.73165187,70.0],
        'up':[0.0,0.0,1.0],
        'normal':[-1.0,0.0,0.0],
        'inner radius':2.0,
        'outer radius':40.0,
        # Location of reference conditions used to calculate thrust from C
        'reference point': [1.0,1.0,1.0],
},

For tree canopies in terrain wind models

'FZ_1':{
        'type':'canopy',
        # Use def when the vtp file specifies a closed volumetric region defining the canopy
        'def':'canopy_dacosta.vtp',
        # Use field when the vtp specifies a set of points on the surface (see below for details)
        'field':'canopy_field.vtp',
        'func':lad_function,
        'cd':0.25,
        'beta_p':0.17,
        'beta_d':3.37,
        'Ceps_4':0.9,
        'Ceps_5':0.9,
},

where the parameters cd, beta_p,beta_d, Ceps_4 and Ceps_5 are defined in the literature (for example Desmond, 2014) and the leaf area density (LAD) is typically specified using a function

def lad_function(cell_centre_list):
    lai = 5.0  # for example following Da Costa 2007
    h_can = 20.0
    lad_list = []
    for cell in cell_centre_list:
        wall_distance = cell[3]
        lad = (np.interp(wall_distance, a_z[0] * h_can, a_z[1]) * lai / h_can,)
        lad_list.append(lad)

    return lad_list

and the variation in leaf area density is a function of height

a_z = (
        np.asarray([0.0,0.43,0.1,0.45,0.2,0.56,0.3,0.74,0.4,1.10,
                  0.5,1.35,0.6,1.48,0.7,1.47,0.8,1.35,0.9,1.01,1.0,0.00]).reshape(11, 2).T
      )

Note

When using the field specification the height of the canopy at each boundary face is set by finding the nearest point to the face centre on the supplied VTK file with the height value looked up in a node based scalar array called ‘Height’

For porous media

'FZ_1':{
        'type':'porous',
        # Use def when the vtp file specifies a closed volumetric region defining the canopy
        'def':'porous_region.vtp',
        # Use zone when the porous region is defined in the mesh as a cell zone
        'zone': [124, 234],
        'alpha': 0.1,
        'c2':0.0,
},

Mesh Transformations

Radial Basis Functions (RBFs) or Inverse Distance Weighted (IDW) interpolation can be used to deform the volume mesh according to the motion of some or all of the boundary nodes. zCFD can perform an arbitrary number of transforms which are specified in consecutively numbered blocks.

'TR_1' : {...},
'TR_2' : {...},

Two kinds of RBF transformation are available: the “Allen et al. greedy” method and the “multiscale” method. For the “greedy” method

'TR_1' : {
     # Required: 'rbf multiscale', 'rbf' or 'idw'
     'type' : 'rbf',
     # Required: list of surface zones to be deformed
     'zone' : [4],
     # Required: function describing surface deformation
     'func' : transformFunc,
     # Required: support radius of RBFs
     'alpha' : 200.0,
     # Optional: type of RBF - one of c0, c2, gaussian, tps (default: c2)
     'basis' : 'c2',
     # Required: maximum error tolerance for greedy
     'tol' : 0.005,
     },

For the “multiscale” method

'TR_1' : {
         # Required: 'rbf multiscale', 'rbf' or 'idw'
         'type' : 'rbf multiscale',
         # Required: faction of surface points to be solved in the base set
         'base point fraction' : 0.1,
         # Required: list of surface zones to be deformed
         'zone' : [4],
         # Required: function describing surface deformation
         'func' : transformFunc,
         # Required: support radius of RBFs in the base set
         'alpha' : 200.0,
         # Optional: type of RBF - one of c0, c2, gaussian, tps (default: c2)
         'basis' : 'c2', # type of RBF - one of c0, c2, gaussian, tps
         },

For the IDW transform

'TR_1' : {
         # Required: 'rbf multiscale', 'rbf' or 'idw'
         'type' : 'idw',
         # Required: number of nearest neighbours to
         'n nearest' : 250,
         # Required: list of surface zones to be deformed
         'zone' : [4],
         # Required: function describing surface deformation
         'func' : transformFunc,
         # Required: power of the interpolation weight in the far field, at the boundary and the blending function between them
         'power' : [3.0, 5.0, 8.0],
         # Required: distance over which to decay the deformation
         'deformation distance': 100.0,
         # Required: Stiffness of the blending function used to decay the deformation (0.0-1.0)
         'blending stiffness': 0.5
         },

An example of a function to deform a surface is given below

def transWing(x,y,z):
 ymax = 1.0
 ymin = 0.1
 yNorm =(y-ymin)/(ymax-ymin)
 x2 = x
 y2 = y
 z2 = z
 if y > 0.1 and y < 1.1 and x < 1.1 and x > -0.1 and z > -0.1 and z < 0.1:
     z2 = z - 0.5*(20.0*yNorm**3 + yNorm**2 + 10.0*yNorm)/31.0

 return {'v1' : x2, 'v2' : y2, 'v3' : z2}

Reporting

In addition to standard flow field outputs (see below), zCFD can provide information at monitor points in the flow domain, and integrated forces across all parallel partitions, any number of 'MR_', 'FR_' or 'MF_' blocks may be spcified.

'report' : {
            # Report frequency
            'frequency' : 1,
            # Extract specified variable at fixed locations
            'monitor' : {
                        # Consecutively numbered blocks
                        'MR_1' : {
                                 # Required: Output name
                                 'name' : 'monitor_1',
                                 # Required: Location
                                 'point' : [1.0, 1.0, 1.0],
                                 # Required: Variables to be extracted
                                 'variables' : ['V','ti'],
                                 },
                        },

            # Report force coefficient in grid axis as well as using user defined transform
            'forces' : {
                       # Consecutively numbered blocks
                       'FR_1' : {
                                # Required: Output name
                                'name' : 'force_1',
                                # Required: Zones to be included
                                'zone' : [1, 2, 3],
                                # Transformation function
                                'transform' : my_transform,
                                # Optional (default 1.0):
                                'reference area' : 0.112032,
                                # Optional (default 1.0):
                                'reference length' : 1.0,
                                # Optional (default [0.0, 0.0, 0.0]):
                                'reference point' : [0.0, 0.0, 0.0],
                                # Optional (default 0.0): Calculate forces using a specified reference pressure rather than absolute
                                'reference pressure': 0.0
                                },
                       },

            # Report massflow ratio
            'massflow' : {
                         # Consecutively numbered blocks
                         'MF_1' : {
                                  # Required: Output name
                                  'name' : 'massflow_1',
                                  # Required: Zones to be included
                                  'zone' : [1, 2 ,3],
                                  },
                         },
        },

Example Transformation Function

A transformation function my be supplied to the force report block to transform the force into a different coordinate system, this a particularly useful for reporting lift and drag which are often rotated from the solver coordinate system.

# Angle of attack
alpha = 10.0
# Transform into wind axis
def my_transform(x,y,z):
    v = [x,y,z]
    v =  zutil.rotate_vector(v,alpha,0.0)
    return {'v1' : v[0], 'v2' : v[1], 'v3' : v[2]}

Overset meshes

Within zCFD in finite volume mode it is possible to run with an arbitrary number of meshes overlapping one another. It is also possible to run a single high order DG mesh (simulation) overlapping a single finite volume background mesh. To activate the overset mesh capability, the user simply strings mesh and case names after the -p and -c options respectively. There should be a python dictionary for each mesh with the appropriate boundary conditions for that mesh defined in it. It should be noted that the ordering of the mesh and case names indicates superiority where the second mesh overlaps the first, and the third overlaps the second and first, and so on. Any mesh which overlaps another must have the boundaries defined where solution information is to be from a lower mesh. This is specified as follows.

'BC_1': {
          'zone': [5],
          'type': 'overset'}
        }

With finite volume mode, the overset solver can also perform simulations with rotating domains overlapping fixed background domains. It is assumed the geometry of the rotating domains is such that there is no new covering or uncovering of cells in the background mesh. Multiple rotating domains can also be defined. If an overset mesh is to rotate, the entire mesh must rotate with no fixed volumes and dual time-stepping time marching must be used. The rotating mesh / domain must be specified in fluid zone 1 as follows:

'FZ_1': {
    "type": "rotating",
    "zone": [0],
    "axis": [0, 0, -1],
    "origin": [0, 0, 0],
    "omega": 100,
},

To control the order of accuracy when mapping between overset meshes, the user can specify the order of accuracy of the interpolation. First (0) or second (1) order mapping can be chosen. First order mapping will use the solution in the cell containing the location requiring data, whereas second order mapping will use an inverse distance weighted interpolation using neighbouring information as well.

'mapping' : {
                # Order of accuracy of mapping.  Defaults to 0.
                'idw from cells order' : 0,
                },

Sliding meshes

Sliding meshes (where interfaces on each mesh slide past one another, for a moving rotor simulation for example) can be simulated in a similar fashion to overset meshes. In fact, sliding meshes are a simple extension of the overset mesh philosophy however no mesh is defined as a background mesh (with cells blanked off). For two meshes sliding across one another, we simply set the sliding boundary as 'overset' in both python dictionaries as described above for BC_1 in the Overset meshes section. The only additional requirement is that to indicate to the solver that the first mesh and solver dictionary combination are not background overset, we must add the following key value pair to the first solver dictionary.

'create overset halos' : True,

The 'create overset halos' keyword instructs the solver to generate an overset halo cell centre by projecting along the face normal a distance equal to the distance from the overset boundary face to the adjacent cell centre. This fictitious cell centre is used for mapping data.

Aero-acoustics

Broadband aero-acoustic simulations can be performed by selecting the DGCAA solver. Aero-acoustic sources are currently limited to stochastic noise sources generated by the FRPM method. These sources are typically driven by a steady state RANS simulation. Several FRPM domains can be defined within the aero-acoustic simulation by adding additional 'FRPM_2', 'FRPM_3', etc., dictionaries. With a single FRPM domain, the dictionary must be labelled 'FRPM_1'.

'DGCAA' : {
             # Order of polynomial used to spatially represent
             # the solution within an element
             'order' : 2,

             # If refraction effects are of interest, these can be included
             # by mapping a CFD solution on to the acoustic domain.  Defaults to True.
             # (Optional)
             'Map CFD to CAA mesh' : True,

             # Option to run the FRPM solver only, without running the acoustic propagation
             # solver.  Defaults to False. (Optional)
             'FRPM solver only' : False,

             # Specify mesh file to be used when mapping RANS simulation to CAA mesh
             # (and / or FRPM mesh)
             'RANS mesh name' : "rans_mesh.h5",

             # Specify solution file to be used when mapping RANS simulation to CAA mesh
             # (and / or FRPM mesh)
             'RANS case name' : "rans_case.py",

             # Specify if FRPM sources are to be used in the acoustic simulation.
             # Defaults to True.  (Optional)
             'frpm sources' : True,

             # Specify after how many cycles microphone data should be extracted.  Note this
             # is cycles rather then Hz.
             'microphone output frequency' : 10,

             # FRPM dictionary.  Specifies FRPM parameters.  Full details below.
             'FRPM' : {
             # See dictionary contents below
             }
}
'FRPM_1' : {
             # The FRPM solver can run on a subset of MPI processes
             # Choose number of MPI processes to run FRPM (Optional)
             'number of frpm mpi processes' : 1,

             # The FRPM can use freestream meanflow or flow from a RANS solution
             # Set to True (default) if RANS solution to be used to drive FRPM (Optional)
             'Map CFD to FRPM mesh' : True,

             # Mesh spacing to be used for implicitly defined FRPM domain
             # Ideally should be less than 0.5 the minimum turbulence integral
             # length scale of interest
             'FRPM Spacing' : 1.0,

             # Constant of proportionality for turbulence integral length scale
             # computation from RANS.  Defaults to 1.0. (Optional)
             'FRPM turbulence integral length scale parameter': 1.0,

             # The user can inspect the RANS solution and choose a turbulence integral
             # length scale best suited to match the RANS simulation.  A uniform
             # integral length scale allows much faster source computations.  Defaults
             # to True (recommended).  Otherwise the turbulence integral length scale
             # will be taken from the mapped RANS solution. (Optional)
             'use constant integral length scale': True,

             # Used in conjunction with 'use constant integral length scale' set to True
             # User specified turbulence integral length scale.
             'FRPM Turbulence Integral Length Scale' : 0.5,

             # Number of cells in each direction for implicity defined FRPM mesh
             # If using parallel FRPM, x-direction should be the largest value for
             # optimum parallel efficiency
             'FRPM Cart Num Mesh Cells' : [10, 10, 10],

             # Frequency with which the FRPM sources should be updated.  Acoustic sources
             # at time levels between updates will be linearly interpolated.  Set to -1 if
             # solver should determine (recommended).  Set to 1 if FRPM should be updated
             # each acoustic solver time step
             'FRPM March Frequency' : 1,

             # Vector the implicitly defined FRPM mesh should be translated to correctly
             # position sources in acoustic domain
             'FRPM Domain Translate' : [0.0, 0.0, 0.0],

             # Vector the implicitly defined FRPM mesh should be rotated to correctly
             # position sources in acoustic domain
             'FRPM Domain Rotate (Deg)' : [0.0, 0.0, 0.0],

             # Physical distance sources should be blended over from the
             # [min-X, max-X, min-Y, max-Y, min-Z, max-Z] cartesian boundaries
             # Values of 0.0 indicate no blending of sources
             'FRPM Blend Sources From Side' : [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],

             # Physical distance from walls acoustic sources should be blended over to zero.
             # Values of 0.0 or less indicates no blending using wall distance
             # (default = -1.0).  The first half of the blending distance has no decay, the
             # second half blends sources linearly to zero.
             'wall distance blend over': -1.0,

             # The FRPM method can reconstruct turbulence fields from many turbulence
             # spectrum models.  Typically Gauss (default) is used. It is also possible to
             # specify Liepmann.  Liepmann is associated with greater computational cost.
             # (Optional)
             'FRPM turbulence spectrum': "gauss",

             # Used in conjunction with the Liepmann spectrum.  Defaults to 1 (Gauss
             # spectrum).  The Liepmann spectrum is obtained by a series of Gauss spectra,
             # the number to be used is specified here (1 to 10). (Optional)
             'FRPM number of discrete filters': 1,


             # Used in conjunction with the Liepmann spectrum.  Defines FRPM mesh refinement
             # as a function of turbulence integral length scale from smaller scales.  0.5 is
             # recommended if using Liepmann. (Optional)
             'FRPM length scale resolution factor': 0.5,

             # Used in conjunction with the Liepmann spectrum.  Maximum turbulence integral
             # length scale to be modelled (multiple of length scale from RANS).
             # Defaults to 5.0. (Optional)
             'FRPM maximum length scale factor': 1.0,

             # Used in conjunction with the Liepmann spectrum.  Minimum turbulence integral
             # length scale to be modelled (multiple of length scale from RANS).
             # Defaults to 0.1. (Optional)
             'FRPM minimum length scale factor': 0.2,

             # Run FRPM solver a number of cycles before starting the acoustic solver.
             # Defaults to 0.
             'FRPM Flush Particle Cycles' : 0,

             # Smooth the turbulence kinetic energy field.  If the RANS solution mapped
             # on to the FRPM domain is discontinuous, that can cause problems with the
             # computation of the acoustic sources.  Ideally a finer RANS mesh should be used,
             # however if this is not possible, Laplace smoothing can be used.  There is
             # potential for the accuracy of the results to deteriorate when using this.
             # Defaults to 0.  (Optional)
                 'FRPM Mapped TKE Smoothing iterations' : 5,

             # Relaxation for turbulence kinetic energy smoothing.  Defaults to 0.2.
             # (Optional)
             'FRPM Mapped TKE Smoothing relaxation' : 0.2,

             # There are two options to map acoustic sources from the FRPM mesh to the
             # acoustic mesh.  Firstly linearly interpolate based on the location of
             # the acoustic mesh solution points within the FRPM domain.  Secondly use
             # in inverse distance weighted mapping based on sources in the neighbourhood
             # of the acoustic solver solution points.  Defaults to True (Recommended).
             # (Optional)
             'FRPM inverse distance source mapping' : True,

             # Used in conjunction with inverse distance weighted source mapping.
             # Distance to specify neighbourhood to take sources from. Recommend
             # at least twice the turbulence integral length scale.  Note, this increases
             # memory usage.
             'frpm inverse distance source mapping distance' : 6.0,

             # List used to specify which boundaries should be used for reseeding particled
             # once they have left the domain [XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX]
             # Defaults to all True.
             'FRPM inflow boundaries': [True, True, True, True, True, True]
}

Output

For solver efficiency, zCFD outputs the raw flow field data (plus any user-defined variables) to each parallel partition without re-combining the output files. The 'write output' dictionary controls solver output.

'write output' : {
                  # Required: Output format: 'vtk', 'ensight', 'native'
                  'format' : 'vtk',
                  # Required: Variables to output on each boundary type
                  'surface variables': ['V','p'],
                  # Required: Field variables to be output
                  'volume variables' : ['V','p'],
                  # Optional: Volume interpolation - interpolates cell values to points in supplied meshes
                  #           File types can be vtp, vtu or binary stl
                  'volume interpolate' : ['mesh.vtp', 'mesh.vtu', 'mesh.stl'],
                  # Optional: Paraview catalyst scripts
                  'scripts' : ['paraview_catalyst1.py','paraview_catalyst2.py'],
                  # Optional: If downstream processes need variables named using a specific convention a naming alias dictionary can be supplied
                  'variable_name_alias' : { "V" : "VELOCITY", },
                  # Required: Output frequency
                  "frequency": {
                      # Output visualisation files every 3 cycles (default = 1000000)
                      "volume data": 3,
                      "surface data": 3,
                      "volume interpolate": 3,
                      "checkpoint": 3,
                      # Start output every 3 cycles after 1 cycle (default = 1)
                      "volume data start": 1
                      "surface data start": 1,
                      "volume interpolate start": 1,
                      "checkpoint start": 1,
                  },
                  # Optional for DG: Output high order surface data on zones
                  'high order surface list' : [1],
                  # Optional (default False):
                  # Don't output vtk volume data
                  'no volume vtk' : False
                  # Optional
                  # Calculate mean and RMS data when running unsteady
                  # Adds 'V_avg", "V_rms", "T_avg", "T_rms", "p_avg" and "p_rms" to the available output variables.
                  'calculate averge and rms': True,
                  # Start time averaging at a specific time step
                  'average start time cycle': 1000
                },

Output Variables

Variable Name Alias Definition
temperature T, t Temperature in (\(K\))
temperaturegrad   Temperature gradient vector
potentialtemperature   Potential temperature in (\(K\))
pressure p Pressure (\(Pa\))
gauge_pressure   Pressure relative to the reference value (\(Pa\))
density rho Density (\(kg/m^3\))
velocity V, v Velocity vector (\(m/s\))
velocitygradient   Velocity gradient tensor
cp   Pressure coefficient
totalcp   Total presure coefficient
mach m Mach number
viscosity mu Dynamic viscosity (\(Pa/s\))
kinematicviscosity nu Kinematic viscosity (\(m^2/s\))
vorticity   Vorticity vector (\(1/s\))
Qcriterion   Q criterion
reynoldstress   Reynolds stress tensor
helicity   Helicity
enstrophy    
ek    
lesregion    
turbulenceintensity ti Turbulence intensity
turbulencekineticenergy   Turbulence kinetic energy (\(J/kg\))
turbulenceeddyfrequency   Turbulence eddy frequency (\(1/s\))
eddy   Eddy viscosity (\(Pa/s\))
cell_velocity    
cellzone    
cellorder    
centre   Cell centre location
viscouslengthscale   Minimum distance between cell and its neighbours
walldistance   Wall distance (\(m\))
sponge   Sponge layer damping coefficient
menterf1   Menter SST blending function

Surface Only Quantities

Variable Name Alias Definition
pressureforce   Pressure component of force on a face
pressuremoment   Pressure component of moment on a face
pressuremomentx   Component of pressure moment around x axis
pressuremomenty   Component of pressure moment around y axis
pressuremomentz   Component of pressure moment around z axis
frictionforce   Skin friction component of force on a face
frictionmoment   Skin friction component of moment on a face
frictionmomentx   Component of skin friction moment around x axis
frictionmomenty   Component of skin friction moment around y axis
frictionmomentz   Component of skin friction moment around z axis
roughness   Wall roughness height
yplus   Non-dinemsional wall distance
zone   Boundary zone number
cf   Skin friction coefficient
facecentre   Face centre location
frictionvelocity ut Friction velocity at the wall