TECHNICAL REFERENCE

SOLIDWORKS FLOW SIMULATION 20Ŋň

ENG
 Contents

1 Physical Capabilities of Flow Simulation . . . . . . . . . . . . . . . . . . . . . . 1

2 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
The Navier-Stokes Equations for Laminar and Turbulent Fluid Flows . . 3
Laminar/turbulent Boundary Layer Model. . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Mass Transfer in Fluid Mixtures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
Constitutive Laws and Thermophysical Properties . . . . . . . . . . . . . . . . . . . 9
Real Gases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
Compressible Liquids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
Non-Newtonian Liquids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Water vapor condensation and relative humidity . . . . . . . . . . . . . . . . . . . . . 18
Cavitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
Equilibrium cavitation model (for pre-defined water only) . . . . . . . . . . . . . . . . . . . . . . . . 19
Isothermal cavitation model (for user-defined liquids only) . . . . . . . . . . . . . . . . . . . . . . . 20
Conjugate Heat Transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
Joule Heating by Electric Current in Solids . . . . . . . . . . . . . . . . . . . . . . . . . 24
Radiation Heat Transfer Between Solids. . . . . . . . . . . . . . . . . . . . . . . . . . . 26
Environment and Solar Radiation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
Discrete Transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
General Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
The Discrete Transfer Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
Clustering and Ray Tracing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
Exchange Factor Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
Discrete Ordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

Flow Simulation 2020 Technical Reference i

 General Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .36
The Discrete Ordinates Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .36
Absorption and Transmission. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
Reflection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
Radiation Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
Radiative Surface and Radiation Source Types . . . . . . . . . . . . . . . . . . . . . . 45
Radiative Surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .45
Radiation Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .46
Simultaneous Use of Radiative Surface and Radiation Source Conditions. . . . . . . . . . . . .47
Viewing Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
Rotation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
Global Rotating Reference Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
Local rotating region(s) (Averaging) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
Local rotating region(s) (Sliding) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
Flows in Porous Media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
General Approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
Perforated Plates in Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 53
Two-phase (fluid + particles) Flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
Euler-Lagrange Approach for the Pre-Processor Simulation . . . . . . . . . . . . 54
Lagrangian Approach for the Post-Processor Simulation . . . . . . . . . . . . . . 58
Free Surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
HVAC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
Tracer Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
Local Mean Age . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
Comfort Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
3 Boundary Conditions and Engineering Devices. . . . . . . . . . . . . . . 71
Internal Flow Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
External Flow Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
Wall Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
Periodic Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
Heat Pipes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

ii Flow Simulation 2020 Technical Reference

 Thermal Joints. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
Two-resistor Components. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
Printed Circuit Boards . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
Thermoelectric Coolers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4 Numerical Solution Technique. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
Computational Mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
Two-Scales Wall Functions Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
Spatial Approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
Spatial Approximations at the Solid/fluid Interface. . . . . . . . . . . . . . . . . . . 88
Temporal Approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
Form of the Numerical Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
Methods to Resolve Linear Algebraic Systems. . . . . . . . . . . . . . . . . . . . . . 90
Iterative Methods for Nonsymmetrical Problems . . . . . . . . . . . . . . . . . . . . 90
Iterative Methods for Symmetric Problems . . . . . . . . . . . . . . . . . . . . . . . . . 90
Multigrid Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
Nested Iterations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
Criterion of mass equation convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
Criterion of energy equation convergence . . . . . . . . . . . . . . . . . . . . . . . . . . 92
Criterion of momentum equation convergence . . . . . . . . . . . . . . . . . . . . . . 93
Notes on the residuals behavior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
5 Mesh Settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
Types of Cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
Mesh Construction Stages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
Basic Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
Control Planes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
Contracting the Basic Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
Mesh Refinement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
Refining Cells by Type . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
Advanced Refinement at Interfaces Between Substances . . . . . . . . . . . . . 104

Flow Simulation 2020 Technical Reference iii

 Small Solid Features Refinement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .105
Curvature Refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .106
Small Solid Feature Refinement Level or Curvature Level. . . . . . . . . . . . . . . . . . . . . . . .108
Tolerance Refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .109
Channel Refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
Thin walls resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
Automatic Mesh Setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
Minimum Gap Size and Minimum Wall Thickness . . . . . . . . . . . . . . . . . . 118
Level of Initial Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
Advanced channel refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
Local Mesh Settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
Recommendations for Creating the Computational Mesh . . . . . . . . . . . 121
Solution-vs.-Mesh Convergence Investigation . . . . . . . . . . . . . . . . . . . . . 123
6 Calculation Control Options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
Finishing the Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
Refinement of the Computational Mesh During Calculation . . . . . . . . . 127
Flow Freezing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
What is Flow Freezing? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
How It Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
Flow Freezing in a Permanent Mode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .132
Flow Freezing in a Periodic Mode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .134
Radiation Freezing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
Radiation Freezing in a Periodic Mode . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
Radiation Freezing in an Automatic Mode. . . . . . . . . . . . . . . . . . . . . . . . . 136
7 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137

iv Flow Simulation 2020 Technical Reference

 1
Physical Capabilities of Flow Simulation

With Flow Simulation it is possible to study a wide range of fluid flow and heat transfer
phenomena that include the following:
• External and internal fluid flows;
• Steady-state and time-dependent fluid flows;
• Compressible gas and incompressible fluid flows;
• Subsonic, transonic, and supersonic gas flows
• Free, forced, and mixed convection;
• Fluid flows with boundary layers, including wall roughness effects;
• Laminar and turbulent fluid flows;
• Multi-species fluids and multi-component solids;
• Fluid flows in models with moving/rotating surfaces and/or parts;
• Heat conduction in fluid, solid and porous media with/without conjugate heat
transfer and/or contact heat resistance between solids and/or radiation heat transfer
between opaque solids (some solids can be considered transparent for radiation),
and/or volume (or surface) heat sources, e.g. due to Peltier effect, etc.

• Joule heating due to direct electric current in electrically conducting solids;1

• Various types of thermal conductivity in solid medium, i.e. isotropic, unidirectional,

biaxial/axisymmetrical, and orthotropic
• Fluid flows and heat transfer in porous media;
• Flows of non-Newtonian liquids;

1. Capabilities and features marked with are available for the Electronics Cooling
module users only.

Flow Simulation 2020 Technical Reference 1

 • Flows of compressible liquids;
• Real gases;
• Cavitation in incompressible water flows;
• Equilibrium volume condensation of water from steam and its influence on fluid
flow and heat transfer;
• Relative humidity in gases and mixtures of gases;
• Two-phase (fluid + particles) flows;
• Periodic boundary conditions.

2
 2
Governing Equations

This chapter provides the theoretical background for the basic governing equations for a
wide range of incompressible and compressible, laminar and turbulent fluid flows and for
transport phenomena such as heat transfer and chemical reactions.

The Navier-Stokes Equations for Laminar and

Turbulent Fluid Flows
Flow Simulation solves the Navier-Stokes equations, which are formulations of mass,
momentum and energy conservation laws for fluid flows. The equations are supplemented
by fluid state equations defining the nature of the fluid, and by empirical dependencies of
fluid density, viscosity and thermal conductivity on temperature. Inelastic non-Newtonian
fluids are considered by introducing a dependency of their dynamic viscosity on flow
shear rate and temperature, and compressible liquids are considered by introducing a
dependency of their density on pressure. A particular problem is finally specified by the
definition of its geometry, boundary and initial conditions.
Flow Simulation is capable of predicting both laminar and turbulent flows. Laminar
flows occur at low values of the Reynolds number, which is defined as the product of
representative scales of velocity and length divided by the kinematic viscosity. When the
Reynolds number exceeds a certain critical value, the flow becomes turbulent, i.e. flow
parameters start to fluctuate randomly.

Flow Simulation 2020 Technical Reference 3

The Navier-Stokes Equations for Laminar and Turbulent Fluid Flows

Most of the fluid flows encountered in engineering practice are turbulent, so Flow
Simulation was mainly developed to simulate and study turbulent flows. To predict
turbulent flows, the Favre-averaged Navier-Stokes equations are used, where
time-averaged effects of the flow turbulence on the flow parameters are considered,
whereas the other, i.e. large-scale, time-dependent phenomena are taken into account
directly. Through this procedure, extra terms known as the Reynolds stresses appear in the
equations for which additional information must be provided. To close this system of
equations, Flow Simulation employs transport equations for the turbulent kinetic energy
and its dissipation rate, the so-called k- model.
Flow Simulation employs one system of equations to describe both laminar and turbulent
flows. Moreover, transition from a laminar to turbulent state and/or vice versa is possible.
Flows in models with moving walls (without changing the model geometry) are computed
by specifying the corresponding boundary conditions. Flows in models with rotating parts
are computed in coordinate systems attached to the models rotating parts, i.e. rotating with
them, so the models' stationary parts must be axisymmetric with respect to the rotation
axis.
The conservation laws for mass, angular momentum and energy in the Cartesian
coordinate system rotating with angular velocity  about an axis passing through the
coordinate system's origin can be written in the conservation form as follows:
 
+ (u i ) = S Mp (2.1)
t  xi
  ui  p 
+ ( u i u j ) + = ( ij   ijR )  S i  S Iip i  1,2,3 (2.2)
t x j  xi  x j
 H ui H  p u
t

xi
=
xi
 
u j ( ij   ijR )  qi    ijR i    Si ui  S Hp  QH
t x j
(2.3)

u2 5  2r 2
H h + k   h m0 y m
2 3 2 m

4
 where u is the fluid velocity,  is the fluid density, S i is a mass-distributed external force
per unit mass due to a porous media resistance (Siporous), a buoyancy (Sigravity = - gi,
where gi is the gravitational acceleration component along the i-th coordinate direction),
and the coordinate system’s rotation (Sirotation), i.e., Si = Siporous + Sigravity + Sirotation,
p p p
h is the thermal enthalpy; S M ,
S I i , S H are additional interfacial exchange terms due to
Euler-Lagrange particle interaction; Q H is a heat source or sink per unit volume,  i j is
the viscous shear stress tensor, q i is the diffusive heat flux,  is an angular velocity of
the rotating coordinate system, r is the distance from a point to rotation axis in rotation
0
reference frame, k is kinetic energy of turbulence, h m is an individual thermal enthalpy of
the m-th mixture component, y m is a concentration of the m-th mixture component and it
is described by the equation (2.17). The subscripts are used to denote summation over the
three coordinate directions.
For calculations with the High Mach number flow option enabled, the following energy
equation is used:

 p
  u i  E  
 E   =  u (   R )  q   R  u i    Q ,
t
 
 xi xi
j ij  ij i ij
x j
 H
(2.4)

u2
E e ,
2
where e is the internal energy.
For Newtonian fluids the viscous shear stress tensor is defined as:

 ui u j2 u 
 ij       ij k  (2.5)
 x j xi 3 xk 

Following Boussinesq assumption, the Reynolds-stress tensor has the following form:

 ui u j 2 uk  2
 ijR  t     ij   k ij
 3
(2.6)

 jx xi 3 x k 

Flow Simulation 2020 Technical Reference 5

The Navier-Stokes Equations for Laminar and Turbulent Fluid Flows

Here  i j is the Kronecker delta function (it is equal to unity when i = j, and zero
otherwise),  is the dynamic viscosity coefficient,  t is the turbulent eddy viscosity
coefficient and k is the turbulent kinetic energy. Note that  t and k are zero for laminar
flows. In the frame of the k- turbulence model,  t is defined using two basic turbulence
properties, namely, the turbulent kinetic energy k and the turbulent dissipation ,

C  k 2
t  f  (2.7)


Here f  is a turbulent viscosity factor. It is defined by the expression

2  20,5  ,

f   1  exp 0.0165R y   1   
RT 
(2.8)


where R  k , R   k y
2


T

y

and y is the distance from the wall. This function allows us to take into account
laminar-turbulent transition.
Two additional transport equations are used to describe the turbulent kinetic energy and
dissipation,

 k   


u i k         t  k   Sk ,

(2.9)
t xi  xi    k   xi 

   


 u i         t     S ,

(2.10)
t  xi xi       xi 

where the source terms S k and S  are defined as

ui
Sk   ijR    t PB (2.11)
x j

  ui   2 . (2.12)
S  C1 f1 ijR  t CB PB   C 2 f 2
k  x j 
 k

6
 Here P B represents the turbulent generation due to buoyancy forces and can be written as

gi 1 
PB   (2.13)
 B  xi

where g i is the component of gravitational acceleration in direction x i , the constant

B = 0.9, and constant C B is defined as: CB = 1 when P B  0 , and 0 otherwise;
3
 0.05  ,
f1  1  
 f 
 f 2  1  exp  RT2  (2.14)
 

The constants C  , C  1 , C  2 ,  k ,   are defined empirically. In Flow Simulation the

following typical values are used:
C = 0.09, C1 = 1.44, C2 = 1.92,  = 1.3,
k = 1 (2.15)

Where Lewis number Le=1 the diffusive heat flux is defined as:

     h , i = 1, 2, 3.
q i    t  (2.16)

 Pr c  xi

Here the constant c = 0.9, Pr is the Prandtl number, and h is the thermal enthalpy.

These equations describe both laminar and turbulent flows. Moreover, transitions from
one case to another and back are possible. The parameters k and  t are zero for purely
laminar flows.

Laminar/turbulent Boundary Layer Model

A laminar/turbulent boundary layer model is used to describe flows in near-wall regions.
The model is based on the so-called Modified Wall Functions approach. This model is
employed to characterize laminar and turbulent flows near the walls, and to describe
transitions from laminar to turbulent flow and vice versa. The modified wall function uses
a Van Driest's profile instead of a logarithmic profile. If the size of the mesh cell near the
wall is more than the boundary layer thickness the integral boundary layer technology is
used. The model provides accurate velocity and temperature boundary conditions for the
above mentioned conservation equations.

Flow Simulation 2020 Technical Reference 7

The Navier-Stokes Equations for Laminar and Turbulent Fluid Flows

Mass Transfer in Fluid Mixtures

The mass transfer in fluid mixtures is governed by species conservation equations. The
equations that describe concentrations of mixture components can be written as

 y m  


u i y m     Dmn  Dmnt y m   S m , m  1,2,..., M
  (2.17)
t x i x i  x i 

t
Here Dmn, D m n are the molecular and turbulent matrices of diffusion, Sm is the rate of
production or consumption of the m-th component.
In case of Fick's diffusion law:

t t
Dmn  D   mn , Dmn   mn  (2.18)

The following obvious algebraic relation between species concentrations takes place:

 ym  1 . (2.19)
m

8
Constitutive Laws and Thermophysical Properties
The system of Navier-Stokes equations is supplemented by definitions of thermophysical
properties and state equations for the fluids. Flow Simulation provides simulations of gas
and liquid flows with density, viscosity, thermal conductivity, specific heats, and species
diffusivities as functions of pressure, temperature and species concentrations in fluid
mixtures, as well as equilibrium volume condensation of water from steam can be taken
into account when simulating steam flows.
Generally, the state equation of a fluid has the following form:
  f P,T ,, y , (2.20)
where y =(y1, ..., yn) is the concentration vector of the fluid mixture components.
Excluding special cases (see below subsections concerning Real Gases, Water vapor
condensation and relative humidity), gases are considered ideal, i.e. having the state
equation of the form
P ,
 (2.21)
RT
where R is the gas constant which is equal to the universal gas constant Runiv divided by
the fluid molecular mass M, or, for the mixtures of ideal gases,
yi
R  Runiv  , (2.22)
i Mi
where y i , i=1, 2, ..., n, are the mass fractions of mixture components, and M i is the
molecular mass of the i-th component.
Specific heat at constant pressure, as well as the thermophysical properties of the gases,
i.e. viscosity and thermal conductivity, are specified as functions of temperature. In
addition, proceeding from Eq. 2.21, each of such gases has constant specific heat ratio
Cp/Cv.
Excluding special cases (see below subsections Compressible Liquids, Non-Newtonian
Liquids), liquids are considered incompressible, i.e. the density of an individual liquid
depends only on temperature:
  f T  , (2.23)

and the state equation for a mixture of liquids is defined as

1
 y 
    i  (2.24)

 i i 
In Flow Simulation the specific volume, i.e. reciprocal to the density, is used in defining
the state of the liquid.

Flow Simulation 2020 Technical Reference 9

Constitutive Laws and Thermophysical Properties

The specific heat and the thermophysical properties of the liquid (i.e. viscosity and
thermal conductivity), are specified as functions of temperature.
The specific heat at constant pressure for a mixture of liquids is defined as
C p   yiC pi (2.25)
i

where C p i is the specific heat of the i-th component.

The dynamic viscosity for a mixture of liquids is defined as
 
  yiv i (2.26)
 i i
where y iv is the volume fraction, and i is the dynamic viscosity of the i-th
component.
The thermal conductivity for a mixture of liquids is defined as
1 2
 
    y i  i 2  (2.27)
 i 
where i is the thermal conductivity of the i-th component.

10
Real Gases
The state equation of ideal gas (2.21) become inaccurate at high pressures or in close
vicinity of the gas-liquid phase transition curve. Taking this into account, a real gas state
equation together with the related equations for thermodynamical and thermophysical
properties should be employed in such conditions. At present, this option may be used
only for a single gas, probably mixed with ideal gases.
In case of user-defined real gas, Flow Simulation uses a custom modification of the
Redlich-Kwong equation, that may be expressed in dimensionless form as follows:
 1 a·F 
Pr  Tr ·   (2.28)
  r  b  r ·( r  c) 
where Pr = P/Pc, Tr = T/Tc, r = Vr·Zc, Vr=V/Vc, F=Tr-1.5, Pc, Tc, and Vc are the
user-specified critical parameters of the gas, i.e. pressure, temperature, and specific
volume at the critical point, and Zc is the gas compressibility factor that also defines the a,
b, and c constants.
A particular case of equation (2.28) with Zc=1/3 (which in turn means that b=c) is the
original Redlich equation as given in Ref. 1.
Alternatively, one of the modifications (Ref. 1) taking into account the dependence of F
on temperature and the Pitzer acentricity factor  may be used:
• the Wilson modification

F  1  1.57  1.62  Tr1  1  (2.29)

• the Barnes-King modification


F  1  0 .9  1.21  Tr1.5  1  (2.30)

• the Soave modification

F
1
Tr
 
1  0.48  1.574  0.176 2 1  Tr0.5  2
(2.31)

The specific heat of real gas at constant pressure (Cp ) is determined as the combination of
the user-specified temperature-dependent "ideal gas" specific heat (Cpideal) and the
automatically calculated correction. The former is a polynomial with user-specified order
and coefficients. The specific heat at constant volume (Cv ) is calculated from Cp by
means of the state equation.
Likewise, the thermophysical properties are defined as a sum of user-specified "basic"
temperature dependency (which describes the corresponding property in extreme case of
low pressure) and the pressure-dependent correction which is calculated automatically.

Flow Simulation 2020 Technical Reference 11

 The basic dependency for dynamic viscosity  of the gas is specified in a power-law form:
 = a T n . The same property for liquid is specified either in a similar power-law form
a 1 T  1 n 
 = a T n or in an exponential form:  = 10 . As for the correction, it is given
according to the Jossi-Stiel-Thodos equation for non-polar gases or the Stiel-Thodos
equations for polar gases (see Ref. 1), judging by the user-specified attribute of polarity.
The basic dependencies for thermal conductivities  of the substance in gaseous and liquid
states are specified by the user either in linear  = a  n  T or in power-law  = a  T n
forms, and the correction is determined from the Stiel-Thodos equations (see Ref. 1).
All user-specified coefficients must be entered in SI unit system, except those for the
exponential form of dynamic viscosity of the liquid, which should be taken exclusively
from Ref. 1. These equations predict the user-defined real gases properties accurately in
the specified range of parameters, including both sub- and supercritical regions (at
pressure up to 1.1Pc).
In case of pre-defined real gas, the custom modification of the Redlich-Kwong equation of
the same form as Eq. 2.28 is used, with the distinction that the coefficients a, b, and c are
specified explicitly as dependencies on Tr in order to reproduce the gas-liquid phase
transition curve at P < Pc and the critical isochore at P > Pc with higher precision:

 1 a 
Pr  Tr ·   (2.32)

 r  b  r ·(  r  c ) 
where Pr = P/Pc, Tr = T/Tc, r = V·Pc/(R·Tc) are the reduced pressure, temperature, and
specific volume respectively, V is the molar volume, Pc and Tc are the critical pressure and
temperature respectively, and R is the gas constant. This equation predicts the pre-defined
real gases properties accurately in the specified range of parameters, including both sub-
and supercritical regions (at pressure up to 2Pc).
When the calculated (p, T) point drops out of the region bounded by the temperature and
pressure limits (zones 1 - 8 on Fig. 2.1) or gets beyond the gas-liquid phase transition
curve (zone 9 on Fig. 2.1), the corresponding warnings are issued and properties of the
real gas are extrapolated linearly.

12
 Fig. 2.1 Pressure-temperature phase diagram.

Supercritical

Liquid

Vapor

If a real gas mixes with ideal gases (at present, mixtures of multiple real gases are not
considered), properties of this mixture are determined as an average weighted with mass
or volume fractions:
N
v   yi vi (2.33)
i 1

where  is the mixture property (i.e., Cp, , or ), N is the total number of the mixture
gases (one of which is real and others are ideal), yi is the mass fraction (when calculating
Cp) or the volume fraction (when calculating  and ) of the i-th gas in the mixture.
The real gas model has the following limitations and assumptions:
• The precision of calculation of thermodynamic properties at nearly-critical
temperatures and supercritical pressures may be lowered to some extent in
comparison to other temperature ranges. Generally speaking, the calculations
involving user-defined real gases at supercritical pressures are not recommended.
• The user-defined dependencies describing the specific heat and transport properties
of the user-defined real gases should be applicable in the whole Tmin...Tmax range
(or, speaking about liquid, in the whole temperature range where the liquid exists).
• Tmin for user-defined real gas should be set at least 5...10 K higher than the triple
point of the substance in question.

Flow Simulation 2020 Technical Reference 13

Compressible Liquids
Compressible liquids whose density depends on pressure and temperature can be
considered within the following approximations:
• obeying the logarithmic law:

 BP 
   0 / 1  C  ln  , (2.34)
 B  P0 
where o is the liquid's density under the reference pressure Po, C, B are
coefficients, here o, C, and B can depend on temperature, P is the calculated
pressure;
• obeying the power law:
1/ n
 PB 
   0    , (2.35)
 P0  B 
where, in addition to the above-mentioned designations, n is a power index which
can be temperature dependent.
Also if the liquid dynamic viscosity depends on pressure, it can be considered within the
following approximation:
   0 e a P  P   , (2.36)

where, in addition to the above-mentioned designations, 0 is the liquid’s dynamic

viscosity under the reference pressure Po, a is a coefficient which can be temperature
dependent and P' = 105 Pa is a constant.

14
Non-Newtonian Liquids
Flow Simulation is capable of computing laminar flows of inelastic non-Newtonian
liquids. In this case the viscous shear stress tensor is defined, instead of Eq. 2.5, as
 u i u j 
 ij        , (2.37)
 x j xi 
where shear rate,
ui u j
  d ij 2  d ii  d jj , d ij   (2.38)
x j xi
and for specifying a viscosity function    the following five models of inelastic
non-Newtonian viscous liquids are available in Flow Simulation:
 Herschel-Bulkley model:
o ,
    K   n1  (2.39)

where K is the liquid's consistency coefficient, n is the liquid's power law index, and
o is the liquid's yield stress.
This model includes the following special cases:

• n = 1,  o = 0 describes Newtonian liquids, in this case K is the liquid's dynamic

viscosity;
• n = 1,  o > 0 describes the Bingham model of non-Newtonian liquids, featured
by a non-zero threshold yield stress (  o ) below of which the liquid behaves as a
solid, so to achieve a flow this threshold shear stress must be exceeded. (In Flow
Simulation this threshold is modeled by automatically equating K, named plastic
viscosity in this case, to a substantially high value at    o );

• 0 < n < 1,  o = 0 describes the power law model of shear-thinning

non-Newtonian liquids (see also below);

• n > 1,  o = 0 describes the power law model of shear-thickening non-Newtonian

liquids (see also below);
 Power-law model:

    K   n1 , (2.40)

in contrast to the Herschel-Bulkley model's special case, the  values are restricted:
min ≤ ≤ max;

Flow Simulation 2020 Technical Reference 15

  Carreau model:


    o     1  Kt   2  n1 / 2 , (2.41)

where:

•  is the liquid's dynamic viscosity at infinite shear rate, i.e., the minimum
dynamic viscosity;
• o is the liquid's dynamic viscosity at zero shear rate, i.e., the maximum
dynamic viscosity;
• Kt is the time constant;
• n is the power law index.
This model is a smooth version of the power law model with the  restrictions.
In all these three models described above, all parameters with the exception of the
dimensionless power law index can be temperature-dependent.
 Cross-WLF model (Cross-William-Landel-Ferry) is another modification of the
power-law model for shear-thinning liquids. It takes into account the effects of the
temperature T:
 0 T , p 
 T , , p   1 n 
  T , p    
(2.42)
1  0 
 * 
where:

 A1 T  T*  
 
•  0 T   D 1  e  A 2  T  T*  
is the zero-shear viscosity or the
'Newtonian limit' in which the viscosity approaches a constant at very low shear
rates;
• T* = D2 is the glass-transition temperature;
• n is the power-law index in the high shear rate regime;
• * is the critical stress level at the transition to shear thinning;
• A1, A2, D1 and D2 are the data-fitted constants.
 Second Order model:

ln   , T   C1  C 2 ln   C 3 ln 2   C 4T  C 5T ln   C 6T 2 (2.43)

where Ci, i = 1…6 are the user-specified coefficients.

16
 The six coefficients Ci can be found by taking the two input viscosity sets i , i 
at T1 and T2 respectively, and operating a simultaneous least-squares
multi-regression minimization. Thus, this equation represents the master curve for
the viscosity representation inside a narrow range of processing temperatures
between the lower limit T1 and the upper limit T2.
Outside of this range, or if the range becomes too wide, its validity falls
progressively down.
The minimum shear rate, below which the viscosity is considered constant, can be
determined automatically as the point where    reaches its maximum value.
The minimum shear rate can also be manually specified by user. The maximum
shear rate, after which the viscosity is considered constant, can be specified by user
also.

 The Viscosity table model defines the liquid's viscosity  , T  by linear
interpolation or polynomial approximation of the user-specified tabular dependencies
of the viscosity  on the shear rate  at the various temperatures T.
Coefficients of the 2nd and 3rd order polynomials are automatically determined with
the least squares method.
The 3rd order polynomials are preferable for shear-thickening liquids which viscosity
rises intensively when the shear rate approaches zero.
In addition, the user can specify the minimum and maximum shear rates, outside of
which the viscosity  is constant. The minimum shear rate for 2nd order polynomials
can be determined automatically as the point where    reaches its maximum value.

Flow Simulation 2020 Technical Reference 17

Water vapor condensation and relative humidity
If the gas whose flow is computed includes steam, Flow Simulation can predict an
equilibrium volume condensation of water from this steam (without any surface
condensation) taking into account the corresponding changes of the steam temperature,
density, enthalpy, specific heat, and sonic velocity.
In accordance with the equilibrium approach, local mass fraction of the condensed water
in the local total mass of the steam and the condensed water is determined from the local
temperature of the fluid, pressure, and, if a multi-component fluid is considered, the local
mass fraction of the steam. Since this model implies an equilibrium conditions, the
condensation has no history, i.e. it is a local fluid property only.
In addition, it is assumed that:
• the volume of the condensed water is neglected, i.e. considered zero, so this
prediction works properly only if the volume fraction of the condensed water does
not exceed 5%,
• The temperature and pressure in the phase transition area should be within the
following ranges: 200.15 K < T < 583.15 K and 103 Pa < P < 107 Pa.
Flow Simulation also allows you to consider the relative Humidity of the gas or mixture of
gases. This allows you to analyze engineering problems where the condensation of water
vapor contained in the air (or other gas), or, more generally speaking, where any
differences in physical properties of wet and dry air play an important role.

18
Cavitation
A liquid subjected to a low pressure below some threshold ruptures and forms vaporous
cavities. More specifically, when the local pressure at a point in the liquid falls below the
saturation pressure at the local temperature, the liquid undergoes phase transition and form
cavities filled with the vapor with an addition of gas that has been dissolved in the liquid.
This phenomenon is called cavitation.
If the characteristic time-scale of the vapor formation process is much less than the
characteristic time-scale of the liquid flow, the cavitation process occurs in conditions
close to thermodynamic equilibrium. Therefore, the equilibrium approach can be used to
simulate the liquid-vapor phase transition.
The following models of cavitation are available in Flow Simulation:
• Equilibrium cavitation model (for pre-defined water only):
This model is based on a homogeneous equilibrium model of cavitation, which
provides a simple technique for analyzing two-phase flows, and is available for
pre-defined water only. It has the capability to account the thermal effects (for
example, localized boiling of water due to intense heating).
• Isothermal cavitation model (for user-defined liquids only):
This model is based on the approach considering isothermal, two-phase flows.

Equilibrium cavitation model (for pre-defined water only)

The homogeneous equilibrium approach is employed. It is applicable for a variety of
important industrial processes.
The fluid is assumed to be a homogeneous gas-liquid mixture with the gaseous phase
consisting of the vapor and non-condensable (dissolved) gas. The vapor mass fraction is
defined at the local equilibrium thermodynamic conditions. The dissolved gas mass
fraction is a constant, which can be modified by user.
The velocities and temperatures of the gaseous (including vapor and non-condensable gas)
and liquid phases are assumed to be the same.
The following assumptions and limitations are made in the cavitation models:
• In the Equilibrium Cavitation Model, cavitation is currently available only for the
pre-defined water (when defining the project fluids you should select Water from
the list of Pre-Defined liquids).
• Cavitation process is assumed to be equilibrium.
• The mixture is assumed to be a homogeneous two-phase fluid., the velocities and
temperatures of the gaseous (including vapor and non-condensable gas) and liquid
phases are assumed to be the same.
• The model does not describe the detailed structure of the cavitation area, i.e.
parameters of individual vapor bubbles and migration of bubbles.

Flow Simulation 2020 Technical Reference 19

 • The temperature and pressure in the phase transition area should be within the
following ranges: 280.15 K < T < 583.15 K, 103 Pa < P < 107 Pa.
• The volume fraction of vapor is limited by 0.95. The parameters of the flow at the
initial and boundary conditions must satisfy this requirement.
• The Cavitation option is not applicable if you calculate a fluid flow in the model
without flow openings (inlet and outlet).
• The fluid region where cavitation occurs should be well resolved by the
computational mesh.
• If the calculation has finished or has been stopped and the Cavitation option has
been enabled or disabled, the calculation cannot be resumed or continued and must
be restarted from the beginning.
The density of the gas-liquid mixture is calculated as:
1 RunivT R Tz (T , P)
 , v  yg  (1  yg  yv )vl (T , P)  yv univ v (2.44)
v P g Pv
where v is the specific volume of the gas-liquid mixture, vl is the specific volume of
liquid, zv (T,P) is the vapor compressibility ratio, Runiv is the universal gas constant, P
is the local static pressure, T is the local temperature, yv is the mass fraction of vapor, v
is the molar mass of vapor, yg is the mass fraction of the non-condensable gas; g is the
molar mass of the non-condensable gas. The properties of the dissolved non-condensable
gas are set to be equal to those of air. By default, the mass fraction of non-condensable gas
is set to 10-4. This is a typical model value appropriated in most cases but it can be
modified by the user in the range of 10-3…10-5.
The mass fraction of vapor yv is computed numerically from the following non-linear
equation for the full enthalpy gas-liquid mixture:
2
I v 5
H  yg hg (T, P)  (1 yg  yv )hl (T, P)  yv hv (T, P)  C  k (2.45)
2 2
where temperature of the mixture T is a function of pressure P and yv. Here hg, hl, hv are
the enthalpies of non-condensable gas, liquid and vapor, respectively, k is the turbulent
energy,
I C  ( u x ) 2  ( u y ) 2  ( u z ) 2 is the squared impulse.

Isothermal cavitation model (for user-defined liquids only)

This model provides a capability to analyze two-phase flows of industrial liquids which
thermophysical properties are not described in details.

20
 The basis of the model is formed from the assumption that the cavitation process occurs in
conditions close to thermodynamic equilibrium. Moreover, it can be assumed that, in most
cases, this process is isothermal, since the medium temperature changes due to phase
transition are insignificant.
Therefore, analyses assume the existence of a barotropic relation for the two-phase
homogeneous mixture. The barotropic state law assumes that the fluid pressure is a
function of fluid density only. The multi-phase fluid density is calculated from the
barotropic law as presented at the diagram below.
Fig. 2.2 Density-pressure phase diagram.

Density
Compressible
liquid

Two-phase
area

Ideal gas

PEV PEL Pressure

The density of the gas-liquid mixture is calculated as:

P  PE T0   g
  , PEV  P  PEL (2.46)
Runiv T0 yg
where Runiv is the universal gas constant, P is the local static pressure, T0 is the local
L
temperature, PE is the saturation pressure of liquid at T0, PE is the local static pressure
V
at which the vapor appears, PE is the local static pressure at which the liquid turns into
the vapor completely, yg is the mass fraction of the non-condensable gas; g is the molar
mass of the non-condensable gas.
V
When the pressure is below PE , a liquid disappears from the mixture and the fluid is
treated as an ideal gas with the molar mass 0. In this case the density of the gas mixture
is calculated as:
1
P  yg 1 yg 
    , P  PEV
(2.47)
Runiv T0   g 0 


Flow Simulation 2020 Technical Reference 21

 L
When the pressure is above the saturation pressure PE , the fluid density equals the
compressible liquid density:

P  PEL
  L
E , P  PEL (2.48)
a2
where  E is the liquid density at the saturation pressure PE , a is the sonic velocity.
L L

In the Isothermal Cavitation Model, the mass fraction of the dissolved (noncondensable)
gas (yg) is a variable value. The four gases can be used as a dissolved gas: Air, Carbon
dioxide, Helium and Methane, in the Equilibrium Cavitation Model Air is used as a
dissolved gas. By default, the mass fraction of the dissolved gas is set to 10-4. This is a
typical value under normal conditions and appropriate in most cases but it can be modified
by the user in the range of 10-2...10-6.
The following additional assumptions and limitations are made in this model:
• The process temperature is constant and the thermal effects are not considered.
• In the Isothermal Cavitation Model, cavitation is currently available only for the
user-defined liquids.
This model requires a minimal number of setting fluid parameters, such as density, molar
mass, saturation pressure, and dynamic viscosity at the local temperature T0. The required
data are available in the literature for the most industrial liquids, such as gasolines, diesel
fuel, mineral and synthetic oils.

22
Conjugate Heat Transfer
Flow Simulation allows to predict simultaneous heat transfer in solid and fluid media with
energy exchange between them. Heat transfer in fluids is described by the energy
conservation equation (2.3) where the heat flux is defined by (2.16). The phenomenon of
anisotropic heat conductivity in solid media is described by the following equation:

e   T 
  i   QH , (2.49)
t xi  xi 
where e is the specific internal energy, e = c·T, c is specific heat, QH is specific heat
release (or absorption) per unit volume, and i are the eigenvalues of the thermal
conductivity tensor. It is supposed that the heat conductivity tensor is diagonal in the
considered coordinate system. For isotropic medium 1 = 2 = 3 = 
If a solid consists of several solids attached to each other, then the thermal contact
resistances between them (on their contact surfaces), specified in the Engineering database
in the form of contact conductance (as m2·K/W), can be taken into account when
calculating the heat conduction in solids. As a result, a solid temperature step appears on
the contact surfaces. In the same manner, i.e. as a thermal contact resistance, a very thin
layer of another material between solids or on a solid in contact with fluid can be taken
into account when calculating the heat conduction in solids, but it is specified by the
material of this layer (its thermal conductivity taken from the Engineering database) and
thickness. The surface heat source (sink) due to Peltier effect may also be considered (see
"Thermoelectric Coolers" on page 79).
The energy exchange between the fluid and solid media is calculated via the heat flux in
the direction normal to the solid/fluid interface taking into account the solid surface
temperature and the fluid boundary layer characteristics, if necessary.
If a solid medium is porous with a fluid flowing through it, then a conjugate heat transfer
problem in this porous-solid/fluid medium can be solved also in the manner described
below (see "Flows in Porous Media" on page 51).

Flow Simulation 2020 Technical Reference 23

Conjugate Heat Transfer

Joule Heating by Electric Current in Solids

This feature is available for the Electronics Cooling module users only.

Flow Simulation is able to calculate steady-state (and quasi time-dependent) direct electric
current in electroconductive solids. In presence of the electric current, the corresponding
specific Joule heat QJ [W/m3] is released and included in QH of heat transfer equation
(2.49) (see "Conjugate Heat Transfer" on page 23). In the case of isotropic material QJ
is
QJ = r·i2, (2.50)

where r is the solids' electrical resistivity [·m] (it can be temperature-dependent) and i is
the electric current density [A/m2].
The electric current density vector

 1  1   1  
i   , , , (2.51)
 r11 x1 r22 x 2 r33 x 3 
is determined via the electric potential  [V]. To obtain the electric potential , Flow
Simulation utilizes the steady-state Laplace equation

  1  
   0 (2.52)
x i  rii  x i 
Here rii is the temperature-dependent electrical resistivity in the i-th coordinate direction.
Transient electric problems with boundary conditions depending on time are considered as
quasi-steady-state. In this case the steady-state problem for potential is solved at each time
step, not taking into account transient electrical processes in solids.
The Laplace equation is solved numerically in the computational subdomain (it may be a
part of the overall computational domain) of electroconductive materials. This
computational subdomain automatically excludes dielectric solids and fluid areas inside.
The total electric current in normal direction over a surface In [A] or electric potential
 [V] may be specified by user as boundary conditions for the problem. These conditions
may be imposed on surfaces between fluid/electroconductive solid, electroconductive
solid/electroconductive solid, dielectric solid/electroconductive solid, and outer solid
surfaces. If no electrical boundary conditions are specified by user, the In = 0 boundary
condition is automatically specified by default on bounding surfaces.
A surface between electroconductive solids in the computational subdomain is either
considered zero-resistance (default) or the electric contact resistance is specified on it. The
resistance value is either given explicitly or calculated from the given material and its
thickness.

24
 A contact resistance specified on a surface implies that the current passing through it
produces the corresponding Joule heating, which yields the following surface heat source
QJS [W/m2]:
QJS = in· (2.53)

where in is the electric current density normal to the surface and  is the electric
potential drop at this surface.
The anisotropic electrical resistivity of electroconductive solids can be anisotropic, i.e.
specified by its components in the coordinate system's directions ri , i = 1, 2, 3. The
isotropic/anisotropic type of material is specified for electrical resistivity and thermal
conductivity simultaneously, i.e. so that their main axes coincide.

Flow Simulation 2020 Technical Reference 25

Radiation Heat Transfer Between Solids

Radiation Heat Transfer Between Solids

In addition to heat conduction in solids, Flow Simulation is capable of calculating
radiation heat transfer between solids. If necessary, a heat radiation from the
computational domain far-field boundaries or the model openings can also be defined and
considered either as thermal, i.e. by specifying the boundaries emissivity and temperature,
or/and as a solar radiation defined by the specified location (on the surface of the Earth)
and time (including date) or by a constant or time-dependent direction and intensity.
Moreover, it is assumed that the project fluids neither emit nor absorb heat radiation (i.e.,
they are transparent to the heat radiation).
In Flow Simulation the following radiation models are employed, each of them has its
own limitations and benefits:
 Discrete Transfer. Its main idea can be described as followed: the radiation leaving
the surface element in a certain range of solid angles can be approximated by a single
ray. The radiation heat is transferred along a series of rays emanating from the radiative
surfaces only. Rays are then traced as they traverse through the fluid and transparent
solid bodies until it hits another radiative surface. This approach, usually called "ray
tracing", allows "exchange factors" to be calculated as the fractions of the total
radiation energy emitted from one of the radiative surfaces that is intercepted by other
radiative surfaces (this quantity is a discrete analog of view factors). "Exchange
factors" between radiative surface mesh elements are calculated at initial stage of
solver, it allows to form matrix of coefficients for a system of linear equations which is
solved on each iteration.

 Discrete Ordinates. This model solves the radiative transfer equation for a finite
number of discrete solid angles, each associated with a vector direction. This method
allows the solution of radiation in absorptive (semi-transparent) media and to model
spectral dependencies. The accuracy of the solution depends on the number of discrete
directions used.1
These models have the following capabilities:
Table 1: Capabilities of radiation models

Discrete Discrete
Transfer Ordinates

Reflection o1 o2

1. Capabilities and features marked with are available for the HVAC module users
only.

26
 Discrete Discrete
Transfer Ordinates

Refraction v o2
Spectral dependencies × v
(wavelength dependency)

Absorption in × v
semi-transparent solids

Geometric optics (focusing on v3 o2

lenses, sharp shadows)

Smooth distribution of the o4 v

radiative heat fluxes

Quadratic Quadratic
dependence on dependence on
View factor Discretization
Parameters on which the resolution level level;
calculation time depends Linear
dependence on
Number of
bands

v – supported
o – supported with limitations
× – not supported
1 – Discrete Transfer is capable for simulating only diffusive reflection or 100% specular
reflection if defined by the Radiative Surface of Symmetry type.
2
– Discretization of space with the ordinates is applicable for majority of tasks though
the method is known as one providing moderate solution for cases where geometric
optics is strong and essential, for example the solar radiation focusing on lenses.
3
– Geometric optics (focusing on lenses, sharp shadows) is accurately simulated with a
directional radiation source (i.e. Solar radiation source) in case the Forward Solar
ray tracing is used. By default the Backward Solar ray tracing is used. For the differ-
ence between Forward and Backward ray tracing see "Exchange Factor Calculation"
on page 34.
4
– Increasing number of rays makes distribution of fluxes less noisy often at the expense
of longer computational time.
When deciding which radiation model to use, consider the following:

Flow Simulation 2020 Technical Reference 27

Radiation Heat Transfer Between Solids

• The Discrete Ordinates model is recommended for most cases with low-temperature
variations and/or with semi-transparent media. At the same time in most cases with
low-temperatures, it is not necessary to consider spectral dependencies (i.e. Number
of bands can be set to 0). Also the Discretization level can be set to 2 or 3 to get an
acceptable accuracy. If the geometric-optical effects of the shadows are significant,
it is recommended to increase the Discretization level to improve the geometrical
characteristics of the model.
If that is not enough, the using of Discrete Transfer model is recommended.
• For cases with high-temperature gradients (i.e. with high power heat sources such as
incandescent lamps, furnaces or other heat sources of temperatures grater then
1000K) and geometric-optical effects of the shadows, focusing and etc. it is
recommended to use the Discrete Transfer model.
But this model does not allow to simulate the absorption and/or spectral
dependencies.

28
Environment and Solar Radiation
Environmental and solar radiation can be applied to external and internal problems. In
fact, the environment radiation is the non-directional energy flux generated by the walls of
an imaginary huge "room" that surrounds the body. This flux has predefined radiation
parameters.
In contrast to the environment radiation, the solar radiation is modeled by the directional
energy flux. Therefore, the solar radiation is defined via its power flow (intensity) and its
directional vector. In addition to the solar radiation from the computational domain
boundaries, a solar radiation source emitting directional radiation can be specified.
The Sun directional vector can be calculated using the following formula:
 
 
 
S  E  sin  s  N  cos  s  cos  s  Z  sin  s (2.54)

where
    
Z is the zenith direction, N is the north direction, E  N  Z is the east direction;
s is the solar elevation angle. That is, the angle between the direction of the geometric
center of the sun's apparent disk and the (idealized) horizon;
s is the solar azimuth angle. That is, the angle from due north in a clockwise direction.
The solar elevation angle s can be calculated, to a good approximation, using the
following formula:
sin  s  cos h cos  cos   sin  sin  (2.55)

where  is the local latitude,  is the current Sun declination,

t
h 2   is the hour angle of the present time t[s] (for example, at solar noon
86400
12:00AM, t = 43200 s and h = 0).

The solar azimuth angle s can be calculated, to a good approximation, using the
following formulas:

sin h cos 
sin  s  
cos  s (2.56)

sin  cos   cos h cos  sin 

cos  s 
cos  s

Flow Simulation 2020 Technical Reference 29

 The current Sun declination  is approximately given by the following expression:

 360 
  23.44  cos   N  10  (2.57)
 365 
where N is the day of the year beginning with N=0 at midnight Coordinated Universal
Time as January 1 begins (i.e. the days part of the ordinal date -1). The number 10, in
(N+10), is the approximate number of days after the December solstice to January 1. The
number 23.44 is the Earth's axial tilt.
The solar radiation intensity can be calculated using the following formula:
Qs0
q s  s    1  n  (2.58)
1  C m m1

1.029
where Cm=0.092 and m1  are empirical values;
sins  0.029
n is the cloudiness, i.e. the presence of clouds in the sky and ranges from 0 (no clouds in
the sky) to 1 (all the sky is covered by clouds);

Qs0 = 1360 W/m2 is the solar constant (i.e. the energy reaching the Earth from the Sun at
the top of the atmosphere).

30
Discrete Transfer
The Discrete Transfer model is a method that has its origin in the flux model but also
exhibits features of the Hottel zone and Monte Carlo techniques.

General Assumptions
 The heat radiation from the solid surfaces, both the emitted and reflected, is assumed
diffuse (except for the symmetry radiative surface type), i.e. obeying the Lambert law,
according to which the radiation intensity per unit projected source area per unit solid
angle is the same in all directions.
 The radiative solid surfaces which are not specified as a blackbody or whitebody are
assumed an ideal graybody, i.e. having a continuous emissive power spectrum similar
to that of blackbody, so their monochromatic emissivity is independent of the emission
wavelength. For certain materials with certain surface conditions, the graybody
emissivity can depend on the surface temperature.
 The solar radiation is absorbed and reflected by surfaces independently from thermal
radiation from all other heat radiation sources.
 The propagating heat radiation passes through a solid specified as radiation transparent
without any absorption. A solid can be specified as transparent to the solar radiation
only, or transparent to the thermal radiation from all sources except the solar radiation,
or transparent to both types of radiation, thermal and solar.

The Discrete Transfer Equations

In a general case, the surfaces participating in the heat radiation (hereinafter radiative
surfaces) can emit, absorb, and reflect a heat radiation (both the solar and thermal). The
heat radiation leaving a radiative surface or radiation source is defined as:

QTout      T 4  A  1     QTin for thermal radiation,

where:  is the surface emissivity,  is the Stefan-Boltzmann constant, T is the
temperature of the surface (     T 4 is the heat flux radiated by this surface in

accordance with the Stefan-Boltzmann law), A is the radiative surface area, QTin is the
incident thermal radiation arriving at this surface;
and


QSout  1     QSin  QSsource  for solar radiation,
in source
where: Q S is the incident solar radiation arriving at this surface, QS is the radiation
arriving at this surface from solar sources.

Flow Simulation 2020 Technical Reference 31

 The net radiation Qnet being the difference between the heat radiation Qout leaving this
surface and the incident heat radiation Qin arriving at it:

  
Qnet  Qout  Qin  QTout  QSout  QTin  QSin ; 
are calculated for each of the surfaces participating in the radiation heat transfer.

Clustering and Ray Tracing

In order to reduce memory requirements, the problem of determining the leaving and net
heat radiation fluxes is solved using a discrete ray Monte-Carlo approach consisting of the
following main elements:
 To reduce the number of radiation rays and, therefore, the required calculation time
and resources, the computational mesh cells containing faces approximating the
radiative surfaces are joined in clusters by a special procedure that takes into account
the face area and angle between normal and face in each partial cell. The cells
intersected by boundaries between radiative surfaces of different emissivity are
considered as belonging to one of these surfaces and cannot be combined in one
cluster. This procedure is executed after constructing the computational mesh before
the calculation and after each solution-adaptive mesh refinement, if any.
 For each of the clusters, the hemisphere governed by the ray’s origin and the normal to
the face at this origin is evenly divided into several nearly equal solid angles generated
by several zenith (latitudinal) angles (at least 3 within the 0...90° range, including
the zero zenith angle of the normal to the face) and several azimuth (longitudinal)
angles (at least 12 within the 0...360° range).
A radiation ray is emitted in each of the solid angles in a direction that is defined
randomly within this solid angle. Each ray is traced through the fluid and transparent
solid bodies until it intercepts the computational domain’s boundary or a cluster
belonging to another radiative surface, thus defining a ‘target’ cluster. Since the
radiation heat is transferred along these rays only, their number and arrangement
govern the accuracy of calculating the radiation heat coming from one radiative surface

32
 to another (naturally, the net heat radiated by a radiative surface does not depend on
number of these rays).
Polar ray

Fig.2.3 Definition of rays emitted from cluster.

The total number of rays emitted from a cluster is calculates as:

N  m  1  n  1,
where m is the number of different latitude values for the rays (including the polar
ray), n is the number of different longitude values (n = 2 for 2D case).
The value of m is defined directly by the View factor resolution level which can be
changed by the user via the Calculation Control Options dialog box. The value of n
depends on m as follows: n  m  4 . The higher the View factor resolution level, the
better the accuracy of the radiation heat transfer calculation, but the calculation time
and required computer resources increase significantly when high values of View
factor resolution level are specified.
To increase the accuracy of heat radiation calculation, the number of radiation rays
emitted from each cluster can be increased automatically during the calculation,
depending on the surface temperature and emissivity, to equalize the radiation heat
emitting through the solid angles.
 When a radiation ray intercepts a cluster of other radiative surfaces, the radiation heat
carried by this ray is evenly distributed over the area of this cluster. The same
procedure is performed if several radiation rays hit the same cluster. To smooth a
possible non-uniformity of the incident radiation heat distribution on a radiative
surface, a fraction of the radiation heat arriving with rays at a cluster can be transferred
to the neighboring clusters also. In addition, small fluctuations are smoothed by the
heat conduction in solid regions.

Flow Simulation 2020 Technical Reference 33

 Exchange Factor Calculation
The "exchange factor" Fij between two clusters is the fraction of the total radiation energy
emitted from one of the clusters that is intercepted by other clusters. Hence, given a point j
at the centre of a cluster on the radiative surface, a radiation ray is fired from j for each one
of the specified directions. Let i be the impingement point, starting from the origin j, the
ray is traced to the i.

in
Then, the incident radiation at point i, Qi , is calculated by adding up the contributions
due to all the rays that reach point i:

Q Tiin  F Q
j
ij
out
Tj is the incident thermal radiation;

Q Siin  F
j
ij Q Sjout is the incident solar radiation.

out
The radiation leaving from the i-th cluster, Qi , is then computed as a sum of the
in
reflected portion of Q i and the emissive power of the surface:

QTiout  1   i QTiin   i Ti 4 Ai for thermal radiation;

 
QSiout  1   Si   QSiin  QSisource for solar radiation.
source
The radiation arriving at i-th cluster from solar sources, Q Si , is calculated as:

QSisource  ISk  FSki  ki for the backward solar ray tracing method;
k

QSisource  QSray  N i for the forward solar ray tracing method.

"Backward" means that rays are emitted from the radiative surfaces. Each solar radiation
source produces one ray that follows the directional vector backward. After it reaches the
outer boundary or the surface having appropriate radiation boundary condition, the
 
exchange factor can be estimated as FSki  nk , ni   Ai . The ISk is the solar radiation
intensity of the k-th solar source. The value  ki equals 1 if the i-th cluster is visible for the
k-th solar source, and 0 if not.

34
 "Forward" means that the rays are emitted from the radiation source. The radiation emitted
by all solar radiation sources can be distributed evenly over all rays traced from radiation
ray
sources. The radiation transferred by each ray, QS , can be estimated by summing up the
source
radiation emitted by all solar sources, Q Sk , and then divided by the number of emitted

Q source
Sk
rays, Nrays: QSray  k .
N rays
The value Ni is the number of rays intercepted the i-th cluster. The forward method is
recommended in case of refraction and reflection radiative surfaces.
The solar ray tracing method and the number of rays traced from radiation sources can be
changed by the user via the Calculation Control Options dialog box.
Please note that for the sake of simplicity the set of equations described here defines the
radiation heat transfer between clusters only and does not take into account the outer
boundaries radiation and diffusive radiation sources. In the full set of equations these
sources are also considered.

Flow Simulation 2020 Technical Reference 35

Discrete Ordinates
The Discrete Ordinates model is a way to approximately solve the radiative transfer
equation (RTE) by discretizing both the xyz-domain and the angular variables that specify
the direction of radiation.

This feature is available for the HVAC module users only.

General Assumptions
 The whole 4 directional domain at any location within the computational domain is
discretized into the specified number of equal solid angles.
 Radiation absorptive (semi-transparent) solids absorb and emit heat radiation in
accordance with the specified solid material absorption coefficient. Scattering is not
considered.
 Surfaces of opaque solids absorb the incident heat radiation in accordance with their
specified emissivity coefficients, the rest incident radiation is reflected specularly or
diffusively, or both specularly and diffusively, in accordance with the specified
specularity coefficient.
 Radiation absorptive solids reflect radiation specularly, the radiation is refracted in
accordance with the specified refraction indices of the solid and adjacent medium
(another radiation absorptive solid, or a transparent solid or fluid, which refraction
index is always considered as equal to 1). The refraction index value cannot exceed 4.

The Discrete Ordinates Method

In the discrete ordinates method the radiation transfer equation is solved for a set of

discrete directions s representing the directional domain of 4 at any position within the

computational domain defined by the position vector r . The directional domain is broken
down into the specified number of equal solid angles or directions. The total number of
directions is defined as:
RL  RL  1
N ord  8  (2.59)
2
where RL is the Discretization level specified by the user. Within each direction the
radiation intensity is considered constant.
In the absence of scattering the radiation transfer equation can be written as:
 
dI s , r 
ds
   

   n 2  I b r   I s , r  (2.60)

36
 T 4
where I is the radiation intensity per solid angle, Ib  is the blackbody radiation

intensity,  is the medium absorption coefficient, n is the refraction index.
At the surfaces of opaque solids the incident heat radiation is absorbed depending on the
specified emissivity coefficient, the rest of the incident radiation is reflected. The opaque
surfaces can also emit heat radiation diffusively in accordance with the surface
temperature and the specified emissivity coefficient.
Since all fluids are considered as transparent to heat radiation, the heat radiation
propagates through them, as well as through transparent solids, without any interaction
with them. However, as the heat radiation is traced through the computational domain by
the discrete ordinates method, the "false scattering" effect caused by the discretization
inaccuracies can appear - it is similar to the "numerical diffusion" effect in fluid flow
calculations. Creating a finer computational mesh allows to reduce this effect.
When the radiation propagates from a small source to a significant distance, the "ray
effect" can be encountered as a result of breaking down the directional domain into several
directions. As radiation in these directions is traced from the source, it begin to
demonstrate a ray-like behavior because more and more cells fall within the same
direction with the distance and the radiation intensity distribution within the direction
becomes non-uniform. This effect must be considered when decreasing the cell size to
avoid the "false scattering" effect, so the discretization level must be increased if the finer
mesh is used.

Flow Simulation 2020 Technical Reference 37

Absorption and Transmission
This feature is available for the HVAC module users only.

The absorptive (semi-transparent) solids absorb heat radiation in accordance with the
specified absorption coefficient.
This option is available for the Discrete Ordinates model only.
In the absence of emitting the radiation transfer equation (2.60) can be written along axis x
as:
dI  x
   I  x (2.61)
dx
The radiation intensity continuously decreases with distance x that the light traverses:
I  x   I 0 e  x (2.62)

where I0 is the initial intensity of radiation and , the absorption coefficient, is the
characteristic of particular material.
At the interface between two transparent or semi-transparent media with different
refractive indices the incident radiation changes its direction in accordance with the
Snell’s law:

sin2 n1
 (2.63)
sin1 n2
where n1 and n2 are the refractive indices of the first and second medium (n is always
equal to 1 for all fluids and for a fully transparent solid material, if its radiation properties
are not specified in the Engineering Database) and 1 and 2 are the incident and
refraction angles correspondingly.
Fig. 2.4 Radiation reflection and transmission at the interface of absorptive media.

I0
1 IR

medium 1
medium 2

2
IT

38
 The radiation reflection at the interface between two transparent or semi-transparent
media with different refractive indices follows the Fresnel’s relation for unpolarized light:

1  tan 2     sin 2      (2.64)

  2 1 2  2 1 2 
2  tan 1   2  sin 1   2  
where  is the reflectivity which define the fraction of reflected radiation.
The fraction of radiation transmitted through the interface is defined as:
 1  (2.65)

If the light is normal to the interface, then:

2
n n  4 n1 n 2
   1 2    (2.66)
 n1  n2  n1  n 2 2
The phenomena of absorption, reflection, and transmission may be applied to the passage
of light through a semi-transparent solid plate. In case of the incident light is normal to the
solid surface, the transmitted intensity can be estimated as shown in Fig. 2.5.
Fig. 2.5 Radiation transmission through a transparent medium.
I’R=I0
I0

I’T=I0(1-)
l
IA=I0(1-)e-l I"R=I0(1-)e-l

IT=I0(1-)2e-l

For an incident beam of intensity I0 that impinges on the front surface of a

semi-transparent solid plate of thickness l and absorption coefficient  the total reflected
intensity IR from the back face can be estimated as:

 1   2 e 2l 
I R  I 0  1  2 2l 
(2.67)

 1  e 

Flow Simulation 2020 Technical Reference 39

 The total transmitted intensity IT at the back face can be estimated as:

IT  I 0
1   2 el  I 0 1    e l
2 (2.68)
2 2l
1  e
The total absorbed intensity IA in the solid plate can be estimated as:

 1   e l 
I A  I 0  I R  I T  I 0 1   1   (2.69)
 1  e l 
So the absorption coefficient  can be expressed as:

1  I  I 1    
   ln  A 0  (2.70)
l  I A   I 0 1    
Thus, the fraction of incident light that is transmitted through a semi-transparent material
depends on the losses that are incurred by absorption and reflection.
So in accordance with the equation (2.68) the absorption coefficient  can be
approximately expressed as:

 IT 
 
1  I0 
   ln (2.71)
l  1   2 
 
 

40
Reflection
At the surfaces of opaque solids the incident heat radiation is absorbed depending on the
specified emissivity coefficient, the rest of the incident radiation is reflected by one of the
following ways or their combination:
• diffusive (Lambertian) reflection,
• diffusive and specular (ideal) reflection which available for the Discrete Ordinates
model only,
The reflection can be expressed as a transfer function that records for a given incoming
direction, the amount of light that is reflected in a certain outgoing direction. The function
is called BRDF (bi-directional reflectance distribution function). This is function of four
angles, two incident (i,i) and two reflected (r,r).
Fig.2.6 Radiation reflection at the surface.
t
ligh
ent
id
Inc n

i
Re

r
fle
cte

i=sinididi
d

d
lig
ht

d r

i

r

The BRDF function bd (i, i, r, r) is defined in terms of incident and reflected
radiance by the following integral:
2  2
Lr  r , r     L  ,   , ,
i i i bd i i r , r cos  i sin  i d i d i (2.72)
0 0

where:
 is the azimuthal angle measured about the surface normal,
Lr(r,r) is the reflected radiance [W/sr/m2],
Li(i,i) is the incident radiance,

Flow Simulation 2020 Technical Reference 41

 bd (i, i, r, r) is the BRDF [sr-1].
Depending on the surface reflective properties the BRDF can be estimated as follows:
• diffusive (Lambertian) surface: light is equally likely to be reflected in any output
direction, independent of the incoming direction.
Fig.2.7 Diffuse reflection.

Incident light n
Reflected light

The BRDF of Lambertian surface is constant:

d
bd  i , i ,  r ,  r   (2.73)

where d is the diffuse reflectance.
• ideal specular surface: the surface that reflects all the light coming from the
direction (i,i) into the direction (i,i+).
Fig.2.8 Ideal specular (mirror) reflection.

Incident light n Reflected light

i r

42
 So the BRDF is zero everywhere except where: r=i and r=i+. Thus the
BRDF is a delta function at direction of ideal mirror reflection:
   i   r   i     r 
 bd  i ,  i , r ,  r   s (2.74)
cos  i sin  i

where s is the specular reflectance.

In case of diffusive and specular (ideal) reflection the BRDF is defined as follows:
    i   r   i     r 
 bd  i ,  i ,  r ,  r   d  s (2.75)
 cos  i sin  i

•
The surface specularity (s) defines the fraction of ideal specular reflected radiation. And
the diffusively reflected fraction is determined as (d = 1 - s).

Flow Simulation 2020 Technical Reference 43

Radiation Spectrum
This feature is available for the HVAC module users only.

Spectrum can be specified for radiation from the computational domain boundaries and
for radiation sources set on surfaces of opaque solids (representing openings).
This option is available for the Discrete Ordinates model only.
In case of using the Discrete Ordinates model, a multiband (Discrete Bands) approach is
used to model spectral dependencies in radiative heat transfer.
The radiation spectrum is considered as consisting of several bands, which edges are
specified by the user. Properties of radiation sources, surfaces and materials are
considered constant within each band. The wavelength-dependent properties of solid
materials are averaged over the specified spectrum bands, so it is recommended to specify
the band edges at the wavelengths, at which the material properties change substantially.
If the radiation spectrum is considered, the equation (2.60) takes the following form:

dIλi s, r 
ds
  

 kλi  n2  Ib λi r   I λi s, r  (2.76)

where Iλi is the heat radiation intensity in the i-th spectrum band, Ib λi is the intensity of
the blackbody radiation in the i-th spectrum band, kλi is the specified the medium
absorption coefficient in the i-th spectrum band.

44
Radiative Surface and Radiation Source Types
Radiative surfaces allow to consider the selected surfaces in contact with the fluid as a
radiative surface for a conjugate heat transfer problem. And radiation sources allow to
consider heating of the model components and solid bodies due to the radiation coming
into the internal space of the model through openings.

Radiative Surfaces
Table 2: Radiative Surface and Radiation Source Parameters.

Radiative surface or
Specified values Dependent values
radiation source type

Wall T, , solar, s  = 1-,  = 

Wall to ambient solar = 1-solar
Wall to environment wall d = 1-s
Symmetry (ideal reflection) No parameters

Absorbent wall (no ray starting) No parameters

Non-radiating wall (no rays) No parameters

A radiative surface can be specified on the surface of either an opaque or a

semi-transparent and a transparent solid.
For the Wall, Wall to ambient and Wall to environment wall boundary condition, the
program gets T from the current results set.
The rays are emitted only from surfaces and boundaries on which the Wall or radiation
source conditions are applied.
Wall to ambient and Wall to environment wall reproduces the most elementary
phenomenon among the radiation effects. The walls with this condition are treated as
black body and do not interact with any other surfaces. They can only exhaust energy into
the space that surrounds the computational domain. Heat flux from such surface could be
calculated as:


Q    T 4  Tout
4
 (2.77)

where T out is the temperature of the environmental radiation. In case of the Wall to
environment wall boundary condition, a particular temperature of the environmental
radiation can be specified at the wall boundary.

Flow Simulation 2020 Technical Reference 45

 When a ray reaches the surface of a radiation source or the Wall to ambient (or Wall to
environment wall) surface, it disappears. All energy carried by this ray also dies away.

The Wall to ambient and Wall to environment wall surface type of radiative
surfaces is not consistent with the Discrete ordinates method, and is substituted in
the calculation with the Wall.1
The Symmetry boundary condition forces the walls to which it is applied to reflect rays as
an ideal mirror.
Surfaces with the specified Absorbent wall boundary condition are taken into account
during the calculation but they can act as absorptive surfaces only. This wall type takes all
heat from the radiation that reaches it and does not emit any heat.
Non-radiating boundary condition removes specific surfaces from the radiation heat
transfer analysis, so they do not affect the results.

The Absorbent wall and Non-radiating surface types of radiative surfaces are
not consistent with the Discrete ordinates method, and are substituted in the
calculation with Whitebody wall, which emissivity is considered equal to 0 (i.e. to
that of whitebody), so that the surface fully reflects all the incident radiation (in
accordance with the Lambert’s law) and does not emit any heat by itself.

Radiation Sources

Diffusive radiation source Power (Q) or Intensity (I) or T and/or spectruma

Solar radiation source n, Power (Q) or Intensity (I) or Tb

a. The spectrum definition is available for the HVAC module users only.
b. The spectrum of a solar radiation source is always the pre-defined Daylight
Spectrum (available for the HVAC module users only).

A radiation source can only be specified on the surface of an opaque solid. Note that when
calculating net radiation rate on such surface, radiation source is not accounted. All the
incident radiation disappears without absorption or reflection on the radiation source
surface, unless there is a radiative surface condition defined on the same surface.

1. Capabilities and features marked with are available for the HVAC module users
only.

46
 The solar radiation source makes the wall to emit radiation like the outer solar radiation. It
is specified by the direction vector and power or intensity or temperature. The solar
radiation at the computational domain boundaries can be specified not only by the
direction vector and intensity, but also by the location (on the surface of the Earth) and
time.

Simultaneous Use of Radiative Surface and Radiation Source Conditions

The user can specify a radiative surface and radiation source on the same surface of an
opaque solid. In this case the total heat radiation leaving the surface is determined as:

q  1   qi      T 4  q source (2.78)

where:  is the emissivity of the radiative surface, qi is the incident heat radiation arriving
on the surface, T is the surface temperature and qsource is the heat radiation emitted by the
radiation source. Please note that the radiative surface emissivity coefficient and
temperature do not influence the heat radiation emitted by the radiation source defined on
the same surface - the radiation source properties are specified independently.

Viewing Results
The main result of the radiation heat transfer calculation is the solids’ surface or internal
temperatures. But these temperatures are also affected by heat conduction in solids and
solid/fluid heat transfer. To see the results of radiation heat transfer calculation only, the
user can view the Leaving radiant flux and distribution of Net radiant flux over the
selected radiative surfaces at Surface Plots. User can also see the maximum, minimum,
and average values of these parameters as well as the Leaving radiation rate and Net
radiation rate as an integral over the selected surfaces in Surface Parameters. All these
parameters can be viewed separately for the solar radiation and radiation from solar
radiation sources (solar radiation) and the radiation from all other heat radiation sources
(thermal radiation).
When the absorptive (semi-transparent) solids are considered, the additional parameters
such as Absorption volume radiant flux, Net volume radiant flow and Net volume
radiant flux become available both for the solar and thermal radiation, as well as for the
total radiation and heat flux.
If the spectral dependencies are specified in terms of discrete bands, the distribution of
radiant heat fluxes (Net and Leaving) over the selected radiative surfaces can be shown in
each band.

Flow Simulation 2020 Technical Reference 47

Rotation

Rotation

Global Rotating Reference Frame

The rotation of the coordinate system is taken into account via the following
mass-distributed force:

S irotation  2 eijk  j u k    2 ri (2.79)

where eijk is the Levy-Civita symbols (function),  is the angular velocity of the
rotation, r is the vector coming to the point under consideration from the nearest point
lying on the rotation axis.

Local rotating region(s) (Averaging)

This option is employed for calculating time-dependent (transient) or steady-state flows in
regions surrounding rotating solids, which are not bodies of revolution (e.g. impellers,
mixers, propellers, etc), when a single global rotating reference cannot be employed. For
example, local rotating regions can be used in analysis of the fluid flow in the model
including several components rotating over different axes and/or at different speeds or if
the computational domain has a non-axisymmetrical (with respect to the rotating
component) outer solid/fluid interface.
In accordance with the employed approach, each rotating solid component is surrounded
by an axisymmetrical (with respect to the component's rotation axis) Rotating region,
which has its own coordinate system rotating together with the component. If the model
includes several rotating solid components having different rotation axes, the rotating
regions surrounding these components must not intersect with each other. The fluid flow
equations in the stationary (non-rotating) regions of the computational domain are solved
in the inertial (non-rotating) Cartesian Global Coordinate System. The influence of the
rotation on the flow is taken into account in each of the rotating coordinate systems.

48
 To connect solutions obtained within the rotating regions and in the non-rotating part of
the computational domain, special internal boundary conditions are set automatically at
the fluid boundaries of the rotating regions. To correctly set these boundary conditions and
the rotating coordinate system, a rotating region must always be a body of revolution. The
rotating region’s boundaries are sliced into rings of equal width as shown on the Fig.2.9.
Then the values of flow parameters, transferred as boundary conditions from the adjacent
fluid regions, are averaged circumferentially over each of these rings.
Computational domain or fluid subdomain
Flow parameters are calculated in the inertial Global Coordinate System

Local rotating region

Flow parameters are calculated
in the local rotating coordinate
system

Rotation axis
Flow parameters are
averaged over these rings

Fig.2.9 The flow parameters averaging on the stator-rotor interface.

To solve the problem, an iterative procedure of adjusting the flow solutions in the rotating
regions and in the adjacent non-rotating regions, therefore in the entire computational
domain, is performed with relaxations.
Please note that even in the case of time-dependent (transient) analysis the flow
parameters within the rotating regions are calculated using a steady-state approach and
averaged on the rotating regions' boundaries as described above.

Local rotating region(s) (Sliding)

This option is employed for calculating time-dependent (transient) flows in cases where
the rotor-stator interaction is strong, when a time-accurate solution for rotor-stator
interaction (rather than a time-averaged solution) is desired.
The Sliding approach assumes that the flow field is unsteady. This assumption allows to
obtain more accurate simulation than by using the Averaging approach, however, that
because this approach requires an unsteady numerical solution, it is computationally more
demanding than the Averaging approach.

Flow Simulation 2020 Technical Reference 49

Rotation

In the Sliding technique rotor and stator cells zones are connected with each other through
"sliding interface". During the calculation, zones linked through "sliding interface" remain
in contact with each other (i.e., the rotor cells “slide” relative to stator cells along the
interface boundary in discrete steps as shown in Figure 2.10).

Initial position of the rotor and stator cells Rotor cells slides with respect to the stator cells
 
( = t)

Fig.2.10 Rotor-stator grid interface

The Sliding approach has the following limitations and/or assumptions:

• The solid body motion is not taken into account in heat transfer analysis.

50
Flows in Porous Media

General Approach
Porous media are treated in Flow Simulation as distributed resistances to fluid flow, so
they can not occupy the whole fluid region or fill the dead-end holes. In addition, if the
Heat conduction in solids option is switched on, the heat transfer between the porous
solid matrix and the fluid flowing through it is also considered. Therefore, the porous
matrix act on the fluid flowing through it via the Si, Siui, and (if heat conduction in solids
is considered) QH terms in Eqs. (2.2) and (2.1), whose components related to porosity are
defined as:
S iporous   k  ij u j , (2.80)

QHporosity    TP  T,  (2.81)

where k is the resistance vector of the porous medium (see below),  is the user-defined
volumetric porous matrix/fluid heat transfer coefficient which can dependent on the flow
velocity, Tp is the temperature of the porous matrix, T is temperature of the fluid flowing
through the matrix, and the other designations are given in Section 1. In addition, the fluid
density in Eqs. (2.1)-(2.3) is multiplied by the porosity n of the porous medium, which is
the volume fraction of the interconnected pores with respect to the total medium volume.
In the employed porous medium model turbulence disappears within a porous medium
and the flow becomes laminar.
If the heat conduction in porous matrix is considered, then, in addition to solving
Eqs. (2.1)-(2.3) describing fluid flow in porous medium, an Eq. (2.49) describing the heat
conduction in solids is also considered within the porous medium. In this equation the
source QH due to heat transfer between the porous matrix and the fluid is defined in the
same manner as in Eq. (2.81), but with the opposite sign. The values of  and c for the
porous matrix may differ from those of the corresponding bulk solid material and hence
must be specified independently. Density of the solid material is multiplied by the solid
volume fraction in the porous matrix, i.e. by (1-n).
Thermal conductivity of the porous matrix can be specified as anisotropic in the same
manner as for the solid material.
The conjugate heat transfer problem in a porous medium is solved under the following
restrictions:
• heat conduction in a porous medium not filled with a fluid is not considered,
• porous media are considered transparent for radiation heat transfer,
• heat sources in the porous matrix can be specified in the forms of heat generation
rate or volumetric heat generation rate only; heat sources in a form of constant or
time-dependent temperature can not be specified.

Flow Simulation 2020 Technical Reference 51

Flows in Porous Media

To perform a calculation in Flow Simulation, you have to specify the following porous
medium properties: the effective porosity of the porous medium, defined as the volume
fraction of the interconnected pores with respect to the total medium volume. Later on, the
permeability type of the porous medium must be chosen among the following:
• isotropic (i.e., the medium permeability is independent of direction),
• unidirectional (i.e., the medium is permeable in one direction only),
• axisymmetrical (i.e., the medium permeability is fully governed by its axial and
transversal components with respect to a specified direction),
• orthotropic (i.e., the general case, when the medium permeability varies with
direction and is fully governed by its three components determined along three
principal directions).
Then you have to specify some constants needed to determine the porous medium
resistance to fluid flow, i.e., vector k defined as k = - gradP  V  , where P, , and V are
fluid pressure, density, and velocity, respectively. It is calculated according to one of the
following formulae:

• k = P S m
  L , where P is the pressure difference between the opposite sides of
a sample parallelepiped porous body, m  is the mass flow rate through the body, S
and L are the body cross-sectional area and length in the selected direction,
respectively. You can specify P as a function of m  , whereas S and L are constants.
Instead of mass flow rate you can specify volume flow rate, v. In this case Flow
Simulation calculates m = v   . All these values do not specify the porous body for
the calculation, but its resistance k only.

• k = P V  L    , where P is the pressure difference between the opposite sides of

a sample parallelepiped porous body, V and  are the fluid velocity and density,
respectively, L is the body length in the selected direction. You can specify P as a
function of V, whereas L is constant, since V and  are calculated.

• k =  AV  B  , where V is the fluid velocity, A and B are constants,  is the fluid
density. Here, only A and B are specified, since V and  are calculated.

•  
k = 32     D2 , where  and  are the fluid dynamic viscosity and density, D is
the reference pore size determined experimentally,  is the porous medium's
porosity. Here, only D and  are specified, since  and  are calculated.

•  
k =  2  D 2  f Re , differing from the previous formula by the f(Re) factor,
yielding a more general formula. Here, in addition to D and , f(Re) as a formula
dependency is specified: f Re = Re   Re  1  , where (Re) is the flow
resistance coefficient of a narrow channel.

52
 The specified reference pore size is also used to calculate the turbulence dissipation after
the porous medium. In case of the two first formulae, the specified pore size is used to
calculate the turbulence dissipation after the porous medium only. By default, it is set to
0.00001 m.To define a certain porous body, specify both the body position in the model
and, if the porous medium has a unidirectional or axisymmetrical permeability, the
reference directions in the porous body.

Perforated Plates in Boundary Conditions

A perforated plate is a particular case of porous media treated in Flow Simulation as
distributed hydraulic resistances. Flow Simulation allows to specify a perforated plate as a
special feature imposed on a pressure-opening- or fan-type flow boundary condition (see
"Internal Flow Boundary Conditions" on page 71) via the Perforated Plates option.
With the use of this feature, the user specifies a perforated plate via its porosity  (defined
as the fraction of the plate area covered by holes) and the holes' shape and size (holes may
be either round with specified diameter, or rectangular with specified width and height, or
regular polygons with specified side length and number of vertices). Then it can be
assigned to any of the above mentioned flow boundary conditions as a distributed
hydraulic resistance yielding the additional pressure drop at this boundary
1
 p    ·  ·u 2 , (2.82)
2
where  is the fluid density, u is the fluid velocity inside the plate's holes, and  is the
perforated plate's hydraulic resistance, calculated according to Ref. 2 from the plate
porosity  and the Reynolds number of the flow in the holes:
Re = ·u·Dh/, (2.83)

where Dh = 4F  is the hole hydraulic diameter defined via the area of a single hole F
and its perimeter , and  is the fluid's dynamic viscosity. Please note that the  (Re, )
dependence taken from Ref. 2 and employed in Flow Simulation is valid for non-swirled,
normal-to-plate flows only.

Flow Simulation 2020 Technical Reference 53

Two-phase (fluid + particles) Flows

Two-phase (fluid + particles) Flows

In Flow Simulation two approaches are employed to simulate fluid-particles flows:
• the Euler-Lagrange approach allows to two-phase gas-particles flows by taking into
account the influence of the discrete second phase (particles) on the first continuous
phase (fluid flow);
• the Lagrangian approach allows to calculate two-phase flows as a motion of
particles in a steady-state flow field where the influence of the particles on the fluid
flow (including its temperature) is negligible.

Euler-Lagrange Approach for the Pre-Processor Simulation

Flow Simulation has capability to calculate two-phase gas-particles flows by using the
hybrid Euler-Lagrange approach combining the continuum description of the gas phase
with the statistically significant sample of individual particles. This involves dividing the
particles into n groups at their injection and then computing subsequent trajectories of the
individual droplets representing these groups in the fluid flow.
The discrete phase may represent non-evaporating solid or liquid particles or evaporating
droplets. Evaporating droplets can be represented by chemically inert substances (water
vapor or Real Gases).
A fundamental assumption made in this model is that discrete phase is sufficiently dilute
that particle-particle interactions and the effects of the particle volume fraction on the fluid
flow are negligible. In practice, it means that the spray occupies a fairly low volume
fraction, usually less than 10-12%, even though high mass loading is acceptable. The
present model does not account for the effects due to particle breakup and coalescence
processes, which might be significant in a dense spray situation and the effects due to
particles radiation. The spray-wall interaction is simulated that when the particle reaches a
wall surface, it disappears.
The two-phase flow field is assumed to be steady-state.
The continuous phase is determined by the system of transport equations (2.1), (2.2), (2.3)
p p p
and (2.17). The droplet-gas interaction is described by source terms S M , S I i , S H and

S mp  v which are taken into account in determining the mass, momentum, energy and
vapor component exchange between the two phases. These source terms may be estimated
from the Lagrangian single droplet behavior calculation and the initial mass flow rate of
the group represented by the single droplet.

54
 The particles parameters are determined by numerical integration of the system of
equations:

 dx pi
  u pi
dt
 du
 pi  u pi
 dt
 dm (2.84)
 p
 m p
 dt
 dT p
  Tp
 dt

where mp is the particle mass, u pi is the particle velocity vector components, Tp is the
particle temperature.
The model describing the droplet evaporation in a multicomponent gas is based on the
continuum description of resistance and heat and mass exchange of single isolated
spherical droplet (Ref. 18). It is assumed that the droplet consists of a single spices liquid
and its temperature is spatially uniform in the droplet. This model is not intended to
calculate the evaporation and growth of droplet in its pure vapor.
The right-hand side of the droplet velocity components equation is defined as follows:

3 C D  Re
u pi  u i  u pi  (2.85)
4  p d p2

where ui is the freestream velocity component and  is the gas phase viscosity, dp and p
are the particle diameter and density, CD is the particle resistance coefficient and Re is the
Reynolds number.
The particle resistance coefficient is calculated as the standard drag coefficient (Ref. 18):

 27Re 0.84 , Re  80
CD   (2.86)
0.271Re
0.217
, 80  Re  104
The right-hand side of the droplet mass equation is defined as follows:

 0.278 Re1 2 Sc1 3 

m p  m 0p 1  12 
  
(2.87)
 1  1.232 Re Sc
43


Flow Simulation 2020 Technical Reference 55

Two-phase (fluid + particles) Flows

where Sc is the Schmidt number and m 0p is the mass flux of the evaporating droplet at
Re = 0.

 1  yv 
m 0p  2d p  D ln  
s 
(2.88)
 1  yv 

where  and D are the gas phase density and vapor diffusion coefficient, yv is the

vapor mass fraction in free stream and y vs is the vapor mass fraction at the droplet
surface, and it is determined, knowing the liquid surface temperature and the pressure,
from the vapor pressure characteristics of the liquid.
The right-hand side of the droplet temperature equation is defined as follows:

m
Tp 
6h
T  Tp   p L (2.89)
 l C pl d p m p C pl

where l and Cpl are the droplet density and specific heat, h is the heat transfer
coefficient, L is the heat of evaporation.

d 3p  l
mp  (2.90)
6
The heat transfer coefficient is calculated as follows:

 0.278 Re 1 2 Pr 1 3 
h  h 0 1  
  
12
(2.91)
 1  1.232 RePr 4 3 

where Pr is the Prandtl number and h0 is the heat transfer coefficient at Re = 0.

m 0p C pv
 d p2
h 
0
(2.92)
 m 0p C pv 
exp   1
 2 d p  
 

where  and C pv are the thermal conductivity and the vapor specific heat.

56
 According to the model describing the droplets behavior in a multicomponent gas, effects
due to changes in thermophysical gas properties across the droplet boundary layer are
taken into account approximately. For this purpose averaged values of temperature and
mass fractions are used. The properties of the gas phase and vapor are evaluated at the
average conditions defined as follows:

T  T p  1   T y m   y ms  1    y m (2.93)

where p and  are the temperatures at the droplet surface and in free stream,

respectively; ym and y ms are the mass fractions of gas phase component in free stream
and at the droplet surface, respectively; and  is the empirical factor that can vary within
0 <  < 1 and which is selected to obtain agreement between model predictions and
existing measurements, wherever possible. By default, the empirical factor is set to 2/3.

Flow Simulation 2020 Technical Reference 57

Two-phase (fluid + particles) Flows

Lagrangian Approach for the Post-Processor Simulation

In addition to simulating two-phase gas-particle flows by taking into account the influence
of the discrete second phase (particles) on the first continuous phase (fluid flow), Flow
Simulation allows to calculate two-phase flows as a motion of particles in a steady-state
flow field where the influence of the particles on the fluid flow (including its temperature)
is negligible. Generally, in this case the particles mass flow rate should be lower than
about 30% of the fluid mass flow rate. To simulate dilute two-phase flows, where flows of
gases or liquids contaminated with particles, the Lagrangian approach is used.
The particles of a specified (liquid or solid) material and constant mass are assumed to be
spherical. The particles trajectories are determined by numerical integration of the
equation:

dU p
U  U  16  d
1   
mp   dReC D p
3
P  m p g (2.94)
dt 8

where mp is the particle mass, U p is the particle velocity vector, d is the particle

diameter, CD is the particle resistance coefficient, g is the gravitational acceleration
vector, and Re is the Reynolds number based on the particle diameter, the relative
velocity, and freestream density and viscosity. The particle resistance coefficient is
calculated with Henderson’s formula (Ref. 3), derived for continuum and rarefied,
subsonic and supersonic, laminar, transient, and turbulent flows over the particles, and
taking into account the temperature difference between the fluid and the particle. The
particles’ rotation, their interaction with each other, Brown motion and additional forces
are not taken into account.
For subsonic flow (M≤1):

1
  Tp 
  3 .65  1 .53 

T exp   0 .247 Re 
C Dsub  24  Re  S  4 .33    
  Tp  S  
  1  0 .353 
 T 

 exp   

 0 .5 M   4 .5  0 .38 0 .03 Re  0 .48 Re  
 0 .1M 2  0 .2 M 8  
(2.95)

 Re   1  0 .03 Re  0 .48 Re 
  M 
 1  exp     0 .6 S
  Re  

58
 For supersonic flow (M≥1.75):

 M   1 
0.5 0.5
0.34 2 1.058  T p 
0.9  2  1.86  2  2     
M  Re   S  S   T  S 4 
C Dsuper   (2.96)
0.5
 M 
1  1.86 
 Re 
For the flow regimes with Mach between 1 and 1.75, a linear interpolation is to be used:

C D  C Dsub 1, Re  
4
3

M   1 C Dsuper 1 .75, Re    C Dsub 1, Re   (2.97)

where CDsub 1, Re is the drag coefficient calculated using the correlation for subsonic flow

with M=1, CDsuper 1.75, Re  is the drag coefficient calculated using the correlation for
supersonic flow with M=1.75, M is the Mach number based on the relative velocity

between the particle and the fluid flow, S  M is the molecular speed ratio,  is the
2
specific heat ratio, T is the fluid temperature in freestream, Tp is the particle temperature,
and subscript ∞ refers to freestream conditions.
Thermal energy equation for a particle is given by:

dTp
mpC p  d  k  NuT  Tp  (2.98)
dt
where Cp is the specific heat of particle, Tp is the particle temperature, k is the thermal
conductivity of fluid, and Nu is the Nusselt number. The particle/fluid heat transfer
coefficient is calculated with the formula proposed in Ref. 4:

2  0.459Re0.55 Pr 0.33
Nu 
 
(2.99)
M
1  3.42 2  0.459Re0.55 Pr 0.33
RePr
If necessary, the gravity is taken into account. Since the particle mass is assumed constant,
the particles cooled or heated by the surrounding fluid change their size.

Flow Simulation 2020 Technical Reference 59

Two-phase (fluid + particles) Flows

The interaction of particles with the model surfaces is taken into account by specifying
either full absorption of the particles (that is typical for liquid droplets impinging on
surfaces at low or moderate velocities) or ideal or non-ideal reflection (that is typical for
solid particles). The ideal reflection denotes that in the impinging plane defined by the
particle velocity vector and the surface normal at the impingement point, the particle
velocity component tangent to surface is conserved, whereas the particle velocity
component normal to surface changes its sign. A non-ideal reflection is specified by the
two particle velocity restitution (reflection) coefficients, en and eτ, determining values of
these particle velocity components after reflection, V2,n and V2,τ, as their ratio to the ones
before the impingement, V1,n and V1,τ:
V2 ,n V2 ,
en  e  (2.100)
V1,n V1,

As a result of particles impingement on a solid surface, the total erosion mass rate,
RΣerosion, and the total accretion mass rate, RΣaccretion, are determined as follows:
N
R erosion    K i  V pb i  f 1 i (  p i )  f 2 i ( d p i )dm p i (2.101)
i 1 M p i

N
R accretion  M pi (2.102)
i 1
where:
N is the number of fractions of particles specified by user as injections in Flow Simulation
(the user may specify several fractions of particles, also called injections, so that the
particle properties at inlet, i.e. temperature, velocity, diameter, mass flow rate, and
material, are constant within one fraction),
i is the fraction number,
Mp i is the mass impinging on the model walls in unit time for the i-th particle fraction,
Ki is the impingement erosion coefficient specified by user for the i-th particle fraction,
Vp i is the impingement velocity for the i-th particle fraction,
b is the user-specified velocity exponent (b = 2 is recommended),
f1 i (αp i) is the user-specified dimensionless function of particle impingement angle αp i,
f2 i (dp i) is the user-specified dimensionless function of particle diameter dp i.

60
Free Surface
Flow Simulation allows to model the flow of two immiscible fluids with a free surface.
Liquids are considered as immiscible if they are completely insoluble in each other. A free
surface is an interface between immiscible fluids, for example, a liquid and a gas (any pair
of fluids belonging to gases, or liquids or non-Newtonian liquids are allowed excluding
gas-gas contact).
However, any phase transitions (including humidity, condensation, cavitation), rotation
problems, surface tension and boundary layer on an interface between immiscible fluids
are not allowed in this version.
Free surfaces are modeled with the Volume of Fluid (VOF) technique by solving a single
set of momentum equations and tracking the volume fraction of each of the fluids
throughout the domain.
The VOF technique is based on the concept of a fluid volume fraction q (q = 0..Nq-1),
which must have a value between 0 and 1. In a two-phase system, for example, in mesh
cells (control volumes) of liquid 0 = 0 and 1 = 1, while in cells of gas 0 = 1 and
1 = 0. The location of a free surface is where q changes from 0 to 1. In each control
volume, the volume fractions of all phases sum to unity:
N q 1


q 0
q 1 (2.103)

where Nq is the number of immiscible fluids (component phases).

Flow Simulation allows to model a two-phase system with Nq = 2 only. Pairs of liquid-gas
or liquid-liquid are available.
The density in each cell is either purely representative of one of the phases or a mixture of
the phases, depending upon the volume fraction values:
N q 1
 
q 0
q q (2.104)

The absence of immiscible fluids mixing means that there is no convective transfer of
mass, momentum, energy and other parameters through the free surface. Also there is no
diffusion transfer of the mixture components. However there are diffusion fluxes for
momentum, energy, and turbulent parameters. It is assumed that velocities, shear stresses,
static pressures, temperatures and heat fluxes are equal at the interface between
immiscible fluids.

Flow Simulation 2020 Technical Reference 61

Free Surface

The terms grad  q  are assumed to be negligible within each fluid. For a gas phase
changes in density q are mainly related to temperature variations, not pressure:

max  q 
 3, q = 0 .. N q  1
min  q 
(2.105)

It is assumed that a gas phase flow is low-compressible, that is Mach number M < 0.3.
At the free surface the ratio between immiscible fluids densities (or viscosities) may be
large, e.g., for water to air the ratio between densities is about 103 and for Non-Newtonian
liquid to air the ratio between viscosities may be up to 1010.
The volume fraction equation is written as follows:
 q  q q  qui
   0, q  0..Nq  1 (2.106)
t q t i xi
Then as a consequence of the volume fraction equations, the conservation law for mass
has the following form:
N q 1
  q  q 
     u i  0

q 0   q  t  i  x i
(2.107)

The momentum equation can be written as:

u i   q  q 
N q 1
  u i         u i u j    p 
t 
q 0  q t  j x j x i
(2.108)

  
 ij   ijR  S i , i  0 ,1, 2
j x j

The energy equation can be written as:

N q 1
  q  q 

H
 H       u i H 
t 
q  0   q t 
 i xi
(2.109)
 p u
 
u j ( ij   ijR )  qi    ijR i     S i u i  QH
ij xi t ij x j i

Each immiscible fluid is one substance, not a mixture.

62
 The state equation of immiscible fluids has the following form:
N q 1
hq q  q
h 
q 0 
(2.110)

The other fluid properties (viscosity, heat conductivity and heat capacity) are defined in
similar way:
N q 1

 
q 0
q q

N q 1

 
q 0
q q (2.111)

N q 1
 q q
Cp  q 0 
C p ,q

The free surface position can be restored at any time for visualization. The restored free
surface position is not used in the numerical algorithm.

Flow Simulation 2020 Technical Reference 63

HVAC

HVAC

Tracer Study
This feature is available for the HVAC module users only.

If a gaseous (or liquid) substance (e.g. a contaminant) diffuses in a gaseous (or liquid)
fluid (if this fluid flows and carries this substance, this fluid is usually named the carrier
fluid) and this substance’s mass fraction y in the carrier fluid is too small, i.e. y<<1, so it
can not influence the carrier fluid flow’s properties (velocity, pressure, temperature), then
distribution of this substance over the computational domain due to carrying it by the fluid
flow and its diffusion in this fluid can be calculated with the Tracer Study option.
According to this option, this substance’s diffusion is calculated in the previously
calculated steady-state or unsteady carrier fluid flow by solving the following Tracer
Study equation taking the substance’s concentration non-uniformity and the carrier fluid
flow’s pressure gradient (for gaseous fluids only) into account:

y   RT  μ μ  y 
  yui    t
 
t xi  pm  Pr  Le Prt  Let  xi  (2.112)

 m1m2  yv1  y  μ μt  p 
2 
   
xi m  p  Pr  Le Prt  Let  xi 
where
 is the carrier fluid and the substance’s mixture density (since y<<1,  can be
considered as the carrier fluid’s density),
t is time,
xi is the i-th component of the used coordinate system,
ui is the i-th component of the carrier fluid’s velocity (the substance has the same
velocity),
p is the carrier fluid’s static pressure,
R is the universal gas constant,
m is the carrier fluid and the substance’s mixture molar mass,
m1 is the substance’s molar mass,
m2 is the carrier fluid’s molar mass,
v1 is the substance’s specific volume,
 is the carrier fluid’s laminar viscosity,
t is the carrier fluid’s turbulent viscosity,
Pr, Prt are the carrier fluid’s laminar and turbulent Prandtl numbers,

64
 Le, Let are the carrier fluid’s laminar and turbulent Lewis numbers.
The Tracer Study equation is solved after finishing the carrier fluid flow calculation, i.e. in
the postprocessor, either as steady-state or as unsteady (time-dependent, with its own
time).
The substance’s physical properties required to solve the Tracer Study equation are
specified in the Engineering Database as the substance’s molar mass (m1), binary (this
 t
substance in the carrier fluid) diffusion coefficient D  ( Dt   D is
Pr Le Prt  Let
assumed) or Lewis number Le (Let = Le, Prt =Pr, t =  are assumed), and saturation
pressure curve vs. temperature (optional).
The Tracer Study equation is solved in the computational domain (or its subdomain) with
the user-specified boundary conditions, initial conditions, and volume sources of the
substance carried by the fluid flow. The substance’s boundary conditions are
superimposed on the carrier fluid flow boundary conditions and consist of the substance’s
surface sources and wall condition specified on the user-selected surfaces. The substance’s
surface sources are specified by the substance’s mass fraction, or mass flow rate, or
specific (i.e. per unit surface area) mass flow rate, or evaporation on the user-selected
walls and/or by the substance’s mass fraction on the user-selected inlet flow openings. At
that the substance’s evaporation can be specified (by enabling the Liquid Surface option)
only for those substances whose saturation pressure curve has been specified in the
Engineering Database. The substance’s wall condition is the substance’s condensation
calculated on the user-selected model walls (only for those substances whose saturation
pressure curve has been specified in the Engineering Database).
The substance’s volume sources are specified by the substance’s mass fraction, or mass
generation rate (total over the source volume), or volumetric (i.e. per unit volume) mass
generation rate in the user-selected fluid volumes.
To solve the Tracer Study equation the user must also specify the substance’s initial
conditions in the form of the substance’s mass fraction (it may be zero) distribution over
the computational domain (or its subdomain).
As a result of the Tracer Study calculation you can see the substance’s mass fraction, mass
flow rates (through the user-selected surfaces or in the user-selected volumes), CRE and
LAQI (see page 69 and page 70) distributions over the computational domain (or its
subdomain).
The Tracer Study can be performed simultaneously for several substances in the same
fluid flow carrying them.

Flow Simulation 2020 Technical Reference 65

HVAC

Local Mean Age

Local mean age (LMA) is the average time for fluid to travel from the selected inlet
opening to the point with taking both the velocity and the diffusion into account. It is
determined by solving the following equation:
3
    t   
 x   u i  
     x   
(2.113)
i 1 i   t  i 
where xi is the i-th coordinate,  is the density, ui is the i-th velocity component,  is the
dynamic viscosity coefficient, t is the turbulent eddy viscosity coefficient,  and t are
the laminar and turbulent Schmidt numbers. The equation is solved under the = 0
boundary condition on the inlet opening.
Dimensionless LMA is the LMA divided by the V/Q ratio, where V is the computational
domain fluid volume, Q is the volume flow rate of the air entering this fluid volume.
LACI (Local Air Change Index) is equal to 1/(Dimensionless LMA), i.e. reciprocal to the
Dimensionless LMA.
By default, the calculation of these parameters is disabled. You can enable the calculation
of LMA, Dimensionless LMA and LACI in the Calculation Control Options.

Comfort Parameters
These parameters are available for the HVAC module users only.

Flow Simulation has the capability to predict the general thermal sensation, degree of
discomfort (thermal dissatisfaction) of people exposed to moderate thermal environments
and estimate air quality by calculating comfort criteria (Ref. 5, Ref. 6, Ref. 7). These
criteria are used when designing occupied spaces and their HVAC systems and are
intended to determine whether environmental conditions are acceptable in terms of
general thermal comfort and air quality or represent discomfort. The calculation of the
comfort criteria assumes that the analyzed fluid is Air.
Mean Radiant Temperature (MRT) is the uniform surface temperature of an imaginary
black enclosure in which an occupant would exchange the same amount of radiant heat as
in the actual non-uniform space.
The mean radiant temperature Tr is defined as follows:
1 1
Tr4 
4  I diffuse ()d 
4
I sun
(2.114)

where Idiffuse is the intensity of the diffuse (thermal) radiation (W/m2/rad), Isun is the
intensity of the solar radiation (W/m2),  is the Stefan-Boltzmann constant.

66
 To calculate the Mean Radiant Temperature, it is assumed that the emissivity of all the
surfaces within the computational domain equals to unity.
Operative Temperature is the uniform temperature of an imaginary black enclosure, in
which an occupant would exchange the same amount of heat by radiation plus convection
as in the actual non-uniform environment.
The operative temperature Tc is defined as follows (Ref. 6):

Tr  T 10V
Tc  (2.115)
1  10V
where Tr is the mean radiant temperature (°C), T is the fluid temperature (°C), V is the
fluid velocity (m/s).
Predicted Mean Vote (PMV) is an index that predicts the mean value of the votes of a
large group of persons on the 7-point thermal sensation scale, based on the heat balance of
the human body. Thermal balance is obtained when the internal heat production in the
body is equal to the loss of heat to the environment. In a moderate environment, the
human thermoregulatory system will automatically attempt to modify skin temperature
and sweat secretion to maintain heat balance (Ref. 7).

cold cool slightly neutral slightly warm hot

cool warm

-3 -2 -1 0 +1 +2 +3

The PMV is defined as follows:

PMV  0.303e 0.036 M  0.028
M  W   3.05 10 3 5733  6.99M  W   pa  
 
 0.42M  W   58.15 1.7  10 M 5867  pa   0.0014M 34  Ta 
5
(2.116)


 
 3.96  10 f cl Tcl  273  Tr  273  f cl hc Tcl  Ta 
8 4 4 

where
Tcl  35.7  0.028( M  W )
 
(2.117)

 I cl 3.96 10 8 f cl (Tcl  273) 4  (Tr  273) 4  f cl hc (Tcl  Ta )

Flow Simulation 2020 Technical Reference 67

HVAC


hc  max 2.38(Tcl  Ta ) 0.25 ,12.1 V 
 m2
 1.00  1.29 I cl , for I cl  0.078 (2.118)
kW
f cl   2
1.05  0.645I , for I  0.078 m
 cl cl
kW
where
M is the metabolic rate (W/m2 of the body area). It is the rate of transformation of
chemical energy into heat and mechanical work by metabolic activities within an
organism (by default it is set to 70 W/m2);
W is the external work (W/m2 of the body area). It accounts for the effective
mechanical power (by default it is set to 0 W/m2);
Icl is the clothing thermal resistance (m2K/W). It is the resistance to sensible heat
transfer provided by a clothing ensemble. The definition of clothing insulation
relates to heat transfer from the whole body and, thus, also includes the uncovered
parts of the body, such as head and hands. The typical values of thermal resistance
for a certain clothing ensemble can be found in Ref. 7 (by default it is set to
0.11 m2K/W);
fcl is the ratio of clothed surface area to nude surface area;
Ta is the air temperature (°C);
Tr is the mean radiant temperature (°C);
V is the relative air velocity (m/s);
pa is the water vapor partial pressure (Pa) calculated in accordance with the
saturation curve, the air temperature Ta and the Relative humidity (see "Water
vapor condensation and relative humidity" on page 18);
hc is the convective heat transfer coefficient (W/m2/K);
Tcl is the clothing surface temperature (°C).
Predicted Percent Dissatisfied (PPD) is an index that provides information on thermal
discomfort or thermal dissatisfaction by predicting the percentage of people likely to feel
too warm or too cool in a given environment (Ref. 7). It can be obtained from the PMV
using the following equation:

PPD  100 - 95 exp - 0.03353 PMV 4 - 0.2179 PMV 2  (2.119)

68
 Draught Rate is the percentage of people predicted to be bothered by draught. The
draught rate DR is defined as follows:

DR  34  T a ,l v a ,l  0 .05  0 .37  v  Tu  3 .14 

0 . 62 (2.120)
a ,l

For v a ,l  0 .05 m/s : v a ,l  0 .05 m / s

For T a ,l  34 C : T a ,l  34 C
For DR  100 % : DR  100 %
where
Ta,l is the local air temperature (°C), 20 °C to 26 °C;
va,l is the local mean air velocity (m/s), < 0.5 m/s;
Tu is the local turbulence intensity (%), 10 % to 60 %. If the Laminar Only flow
type is selected, Tu equals 40 %.
This model applies to people at light, mainly sedentary activity with a thermal sensation
for the whole body close to neutral and for prediction of draught at the neck. At the level
of arms and feet, the model could overestimate the predicted draught rate. The sensation
of draught is lower at activities higher than sedentary (> 1.2 met) and for people feeling
warmer than neutral.
Draft Temperature is the difference in temperature between any point in the occupied
zone and the control condition. "Draft" is defined as any localized feeling of coolness or
warmth of any portion of the body due to both air movement and air temperature, with
humidity and radiation considered constant (Ref. 5).
The draft temperature Td is defined as follows:
Td  T -Tm - 7.6553 (V -0.1524 ) (2.121)

where:
T is the local fluid temperature (°C);
Tm is the average fluid temperature within the control space (°C);
V is the local fluid velocity (m/s).
Air Diffusion Performance Index (ADPI) is the percentage of the space in which the air
speed is less than 0.35 m/s and the Draft Temperature falls between -1.7 ºC and 1.1 ºC
(Ref. 5).
If the Draft Temperature or ADPI are calculated in Volume Parameters, then the “control
space” will correspond to the specified volume region. In all other cases the “control
space” corresponds to the whole computational domain.
Contaminant Removal Effectiveness (CRE) is an index that provides information on the
effectiveness of a ventilation system in removing contaminated air from the whole space.
For a perfect mixing system CRE=1. Values above 1 are good, values below 1 are poor.

Flow Simulation 2020 Technical Reference 69

HVAC

This parameter is only available if more than one fluid is present in the control space.
Its value is defined as follows:
Ce
CRE  (2.122)
C
where:
Ce is the bulk average mass fraction of the contaminant calculated over all faces
through which the fluids outflow from the computational domain;
<C> is the bulk average mass fraction of the contaminant calculated over the whole
computational domain.
Local Air Quality Index (LAQI) is an index that provides information on the effectiveness
of a ventilation system in removing contaminated air from a point. For a perfect mixing
system LAQI =1. For other systems, the higher the value at a point, the better the
capability of the ventilation system in removing contaminated air from that point. This
parameter is only available if more than one fluid is present in the control space. Its value
is defined as follows:
Ce
LAQI  (2.123)
C
where:
Ce is the bulk average mass fraction of the contaminant calculated over all faces
through which the fluids outflow from the computational domain;
C is the mass fraction of the contaminant at a point.
The Flow Angle shows the flow deviation from the design flow direction. Consider one of
the axis of selected coordinate system as the design flow direction, the results can then be
viewed as the deviation from the design. Typically, flow angles of less than 15° might be
considered as good.
The flow angle components are defined as follows:
Angle (X)  arccos(V x /V)
Angle (Y)  arccos(V y /V) (2.124)

w Angle (Z)  arccos(V z /V)

where Vx, Vy, Vz are the X, Y, and Z components of the fluid velocity and V is the
absolute value of the fluid velocity vector.

70
 3
Boundary Conditions and Engineering Devices

This chapter describes flow boundaries conditions implemented in Flow Simulation and a
wide range of the engineering device models specifying conditions at the model surfaces
and in solid components.

Internal Flow Boundary Conditions

For internal flows, i.e., flows inside models, Flow Simulation offers the following two
options of specifying the flow boundary conditions: manually at the model inlets and
outlets (i.e. model openings), or to specify them by transferring the results obtained in
another Flow Simulation calculation in the same coordinate system (if necessary, the
calculation can be performed with another model, the only requirement is the flow regions
at the boundaries must coincide).
With the first option, all the model openings are classified into "pressure" openings,
"flow" openings, and "fans", depending on the flow boundary conditions which you
intend to specify on them.
A "pressure" opening boundary condition, which can be static pressure, or total pressure,
or environment pressure is imposed in the general case when the flow direction and/or
magnitude at the model opening are not known a priori, so they are to be calculated as part
of the solution. Which of these parameters is specified depends on which one of them is
known. In most cases the static pressure is not known, whereas if the opening connects the
computational domain to an external space with known pressure, the total pressure at the
opening is known. The Environment pressure condition is interpreted by Flow Simulation
as a total pressure for incoming flows and as a static pressure for outgoing flows. If,
during calculation, a vortex crosses an opening with the Environment pressure condition
specified at it, this pressure considered as the total pressure at the part of opening through
which the flow enters the model and as the static pressure at the part of opening through
which the flow leaves the model.

Flow Simulation 2020 Technical Reference 71

External Flow Boundary Conditions

Note that when inlet flow occurs at the "pressure" opening, the temperature, fluid mixture
composition and turbulence parameters have to be specified also.
A "flow" opening boundary condition is imposed when dynamic flow properties (i.e., the
flow direction and mass/volume flow rate or velocity/ Mach number) are known at the
opening. If the flow enters the model, then the inlet temperature, fluid mixture
composition and turbulence parameters must be specified also. The pressure at the
opening will be determined as part of the solution. For supersonic flows the inlet pressure
must be specified also.
A "fan" condition simulates a fan installed at a model opening. In this case the
dependency of volume flow rate on pressure drop over the fan is prescribed at the opening.
These dependencies are commonly provided in the technical documentation for the fans
being simulated.
With the second option, you specify the boundary conditions by transferring the results
obtained in another Flow Simulation calculation in the same coordinate system. If
necessary, the calculation can be performed with another model, the only requirement is
the flow regions at the boundaries must coincide. At that, you select the created boundary
conditions’ type: either as for external flows (so-called "ambient" conditions, see the next
Section), or as for "pressure" or "flow" openings, see above. If a conjugate heat transfer
problem is solved, the temperature at the part of the boundary lying in a solid body is
transferred from the other calculation.
Naturally, the flow boundary conditions specified for an internal flow problem with the
first and/or second options must be physically consistent with each other, so it is expedient
to specify at least one "pressure"-type boundary condition and at least one "flow"-type
boundary condition, if not only "ambient" boundary conditions are specified.
If one or several non-intersecting axisymmetric rotating regions (local rotating reference
frames) are specified, the flow parameters are transferred from the adjacent fluid regions
and circumferentially averaged over rotating regions’ boundaries as boundary conditions.

External Flow Boundary Conditions

For external problems such as flow over an aircraft or building, the parameters of the
external incoming flow (so-called "ambient" conditions) must be defined. Namely the
velocity, pressure, temperature, fluid mixture composition and turbulence parameters must
be specified. Evidently, during the calculation they can be partly violated at the flow
boundary lying downstream of the model.

72
Wall Boundary Conditions
In Flow Simulation the default velocity boundary condition at solid walls corresponds to
the well-known no-slip condition. The solid walls are also considered to be impermeable.
In addition, the wall surface's translation and/or rotation (without changing the model's
geometry) can be specified. If a calculation is performed in a rotating coordinate system,
then some of the wall surfaces can be specified as stationary, i.e. a backward rotation in
this coordinate system (without changing the model geometry). Flow Simulation also
provides the "Ideal Wall" condition that corresponds to the well-known slip condition. For
example, Ideal Walls can be used to model planes of flow symmetry.
If the flow of non-Newtonian liquids is considered, then the following slip condition at all
solid walls can be specified for each non-Newtonian liquid in the project separately: if the
shear stress  exceeds the yield stress value 0,slip, then
a slip velocity vslip is determined from vslip  C1 (τ  τ0,slip )C2 , where C1 and C2, as well
as 0,slip, are specified by user.
If conjugate heat transfer in fluid and solid media is not considered, one of the following
boundary conditions can be imposed at solid walls: either the wall temperature
T  Tw , (3.1)
or the heat flux,
q  qw (3.2)
being positive for heat flows from fluid to solid, equal to zero for adiabatic
(heat-insulated) walls, and negative for heat flows from solid to fluid.
When considering conjugate heat transfer in fluid and solid media, the heat exchange
between fluid and solid is calculated by Flow Simulation, so heat wall boundary
conditions are not specified at the walls.

Periodic Boundary Conditions

The "periodicity"condition may be applied if the model consist of identical geometrical
features arranged in periodic linear order. Periodic boundary conditions are specified at
the pair of computational domain boundaries for the selected direction in which a
geometrical feature or a group of features repeats regularly over distance. Periodic
boundary conditions allows to reduce the analysis time by calculating the fluid flow only
for a small group of identical geometrical features or even just for one feature, but taking
into account influence of other identical features in the pattern. Please note that the
number of basic mesh cells along the direction in which the "periodicity"condition is
applied must be no less than five.

Flow Simulation 2020 Technical Reference 73

Heat Pipes

Heat Pipes

This feature is available for the Electronics Cooling module users only.

"Heat Pipe" is a device (of arbitrary form: a tube, a plate, etc.) transferring heat from its
hotter surface to its colder surface due to evaporating a liquid (water, etc.) and condensing
its vapor in this device's inner hollow. This liquid evaporates near this hollow’s hotter
surface and its vapor condensates near this hollow’s colder surface. The condensed liquid
then returns to the hotter end due to gravity or by a wick. If such a device operates under
its design conditions, the heat transfer rate provided by this device exceeds substantially
any value achievable by a solid body of similar dimensions.
In Flow Simulation a "Heat Pipe" is modeled simplistically as an extremely
heat-conducting body. To model real efficiency of a heat pipe, a non-zero thermal
resistance can be assigned to it.

Thermal Joints
A thermal joint model implemented in Flow Simulation can be used to simplify the heat
transfer simulation from one (hot) surface S1 to the other (cold) surface S2 with the actual
thermal joint geometry excluded from the analysis. Here, to simulate heat transfer
between these surfaces, the value of contact resistance  (or its reciprocal, heat transfer
coefficient ) is used. The resulting heat transfer rate from the hot surface to the cold
surface is calculated as Q   T1  T 2  , where T1 and T2 are the temperatures of the
surfaces S1 and S2.
Notice that once you select the faces to specify thermal joint between them, these faces
become thermally insulated in respect to the surrounding medium and only participate in
the heat exchange between each other.

74
Two-resistor Components

This feature is available for the Electronics Cooling module users only.

A two-resistor model implemented in Flow Simulation can be used for simplified solving
of heat transfer problems in an electronic device with small electronic packages (chips,
etc.). A small package is considered as a flat solid plate, which is mounted on the printed
circuit board (see Fig. 3.1). The compact model consists of three nodes: Junction, Case
and Board. These are connected together by two thermal resistors which are the
user-specified values of the junction-to-board JB (from the junction to the board on
which it is mounted) and junction-to-case JC (from the junction to the top surface of the
package) thermal resistances (in K/W in SI).
 JB  TJ  TB  PH (3.3)

 JC  TJ  TC  PH (3.4)

where TJ is the junction temperature, TB is the board temperature and TC is the case
temperature, measured at the package top surface, PH is the heat generation rate which
produced the change in junction temperature and applied at the Junction.
The Board is considered to be in direct thermal contact with the local environment
immediately below the footprint of the package; normally the printed circuit board
(PCB).
The Case is considered to be in direct thermal contact with the local environment
immediately above the top of the package (normally air, or a thermal interface material
used in conjunction with a heat sink).
The Junction plate is modeled as a high conductivity body with heat-insulating side walls,
so the only surfaces through which heat can be exchanged with the environment are the
top and bottom surfaces
Fig. 3.1 Flow Simulation representation of the two-resistor model.

Top Plate (Case, qJC)

High Conductivity
Block (Junction)

PCB
Bottom Plate (Board,JB)
Adiabatic walls

Flow Simulation 2020 Technical Reference 75

Printed Circuit Boards

Printed Circuit Boards

This feature is available for the Electronics Cooling module users only.

Printed Circuit Board (PCB) is concerned as a special case of solid material having
anisotropic thermal conductivity. The integral characteristics of a PCB, i.e. its effective
density, specific heat, and components of thermal conductivity, are calculated on the basis
of the PCB structure.
To define a PCB characteristics, you have to specify density, specific heat, and thermal
conductivity for both dielectric (denoted by 'D' index) and conductor (denoted by 'C'
index) materials of the PCB, and describe its internal structure in one of the provided
types:
• Conductor Volume Fraction type requires the specification of the volume fraction
of conductor material in the PCB A. The in-plane (planar) conductivity Kin and the
through-plane (normal) conductivity Kthro are calculated as follows:

A A
1
A  A  1
Kin   KC  1    KD ,  100  100 (3.5)
100  100 K thro KC KD
The effective density  and the effective specific heat C of the board are calculated
as follows:

C  VC  D  VD CC  C  VC  CD   D  VD
 , C (3.6)
V  V

A
where V is the total volume of the PCB, VC  V  is the volume of the
100

conductor material in the PCB, V D  V  1 
A 
 is the volume of the dielectric
 100 
material in the PCB;
• Board Mass type allows to calculate the fraction of conductor material A in PCB by
using the PCB total mass M and PCB total volume V:

  C
A  100 (3.7)
C  D

76
 The effective density  is calculated as follows:
M
  (3.8)
V
The in-plane (planar) conductivity Kin and the through-plane (normal) conductivity
Kthro are calculated in the same manner;
• Layer Definition type implies that you must specify the PCB total thickness t and
the number of conducting layers nC. Also for each conducting layer you must
specify the volume fraction of conductor material in the layer Ai and the layer
thickness tCi . The in-plane (planar) conductivity Kin and the through-plane
(normal) conductivity Kthro are calculated as follows:

nC
 Ai  Ai    nC 
 t Ci 
100
 KC  

1 
100
 D    t  tCi   KD

 K
Kin 
i 1    i1  (3.9)
t
t
K thro  (3.10)
Ai    
nC
 Ai 
nC
  1    t  
100    i 1 
 t Ci

 100 

i 1
tCi 
K

KD 

KD
 C

 
The effective density  and the effective specific heat C of the board are calculated
as follows:

nC
 Ai  Ai    nC 
 tCi 
100
 C  1 
 100 
   D    t  tCi    D

i 1    i1  (3.11)
t
nC
 Ai  Ai    nC 
 tCi   C CC 1  D CD  t tCi  D CD
C
i1 100  100   i1  (3.12)
 t
The fraction of conductor material A in PCB is calculated as follows:

Flow Simulation 2020 Technical Reference 77

Printed Circuit Boards

nC
Ai
t Ci 
100
A i 1
 100 (3.13)
t

78
Thermoelectric Coolers
Thermoelectric cooler (TEC) is a flat sandwich consisting of two plates covering a circuit
of p-n semiconductor junctions inside. When a direct electric current (DC) i runs through
this circuit, in accordance with the Peltier effect the a  i  Tc heat, where a is the Seebeck
coefficient, Tc is the TEC's "cold" surface temperature, is pumped from the TEC's "cold"
surface to its "hot" surface (the "cold" and "hot" sides are determined from the DC
direction). This heat pumping is naturally accompanied by the Joule (ohmic) heat release
at both the TEC surfaces and the heat transfer from the hotter side to the colder (reverse to
i2
the Peltier effect). The ohmic heat release is defined as R , where R is the TEC's
2
electric resistance, while the heat transfer is defined as k·T, where k is the TEC's thermal
conductivity, T  Th  Tc , Th is the TEC's "hot" surface temperature. The net heat
transferred from the TEC's "cold" surface to its "hot" surface, Qc, is equal to

Qc  a  i  Tc  R  i 2 / 2  k  T (3.14)

Correspondingly, the net heat released at the TEC's "hot" surface, Qh, is equal to

Qh  a  i  Th  R  i 2 / 2  k  T (3.15)

Flow Simulation 2020 Technical Reference 79

Thermoelectric Coolers

In Flow Simulation a TEC is specified by selecting a flat plate (box) in the model,
assigning its "hot" face, and applying one of the TECs already defined by user in the
Engineering Database.

Heat Sink

Hot Side

TEC
Cold Side

Object being cooled

Fig.3.2 Flow Simulation representation of the Thermoelectric Cooler model.

The following characteristics of TEC are specified in the Engineering Database:

• Maximum pumped heat is the maximum heat Qc,max transferred at this imax from
the "cold" face to the "hot" face at temperature difference between these faces
T = 0;
• Maximum temperature drop is the maximum temperature difference Tmax
between the "hot" and "cold" faces attained at the heat transferred between these
faces Qc = 0;

• Maximum current strength is the maximum DC current, imax;

• Maximum voltage is the voltage Vmax corresponding to imax.

80
 All of these characteristics are specified for two Th values, in accordance with the
information usually provided by the TEC suppliers. Proceeding from these characteristics,
the a(T), R(T), and k(T) linear functions are determined. The functional boundary
conditions are specified automatically on the TEC's "cold" and "hot" surfaces, which must
be free from other boundary conditions.
The temperature solution inside the TEC and on its surfaces is obtained using a special
procedure differing from the standard Flow Simulation calculation procedure for heat
conduction in solids.
The TEC's "hot" face must be in contact with other solids, i.e it must not be in contact with
any fluid. In addition, it is required that the obtained TEC solution, i.e. Th and T, lie
within the TEC's operating range specified by its manufacturer.

Flow Simulation 2020 Technical Reference 81

Thermoelectric Coolers

82
 4
Numerical Solution Technique

The numerical solution technique employed in Flow Simulation is robust and reliable, so
it does not require any user knowledge about the computational mesh and the numerical
methods employed. But sometimes, if the model and/or the problem being solved are too
complicated, so that the Flow Simulation standard numerical solution technique requires
extremely high computer resources (memory and/or CPU time) which are not available, it
is expedient to employ Flow Simulation options which allow the adjustment of the
automatically specified values of parameters governing the numerical solution technique.
To employ these options properly and successfully, take into account the information
presented below about Flow Simulation’ numerical solution technique.
Briefly, Flow Simulation solves the governing equations with a discrete numerical
technique based on the finite volume (FV) method. Cartesian rectangular coordinate
system is used. To obtain space discretization, the axis-oriented rectangular grid is used
far from a geometry boundary. Thus, the control volumes (i.e. mesh cells) are rectangular
parallelepipeds. Near the geometry boundary Cartesian cut cells approach is used.
According to this approach, the near-boundary mesh is obtained from the original
background Cartesian mesh by cutting original parallelepiped cells that intersect the
geometry (see Fig.4.1). Consequently, the near-boundary cells are polyhedrons with both
axis-oriented and arbitrary oriented plane faces in this case. Thus, Flow Simulation
combines advantages of approaches based on regular grids and ones with highly accurate
representation of geometry boundaries.
Also the local refinement of mesh is used in Flow Simulation to take into account
geometry and solution peculiarities. It is usually exploited at the solid/fluid interface, in
regions of high gradients, and so on.
All physical parameters are referred to control volume mass centers. Following the FV
approach, the direct discretization of the integral form of the conservation laws is used. It
guarantees that the basic quantities mass, momentum and energy remain conserved in the
discrete representation.

Flow Simulation 2020 Technical Reference 83

Computational Mesh

The spatial derivatives are approximated with implicit difference operators of

second-order accuracy. The time derivatives are approximated with an implicit first-order
Euler scheme. The numerical viscosity generated by the discretization error of the scheme
is small enough and allows obtaining adequately accurate results in the most practical
cases.

Computational Mesh
Flow Simulation computational approach is based on locally refined rectangular mesh
near geometry boundaries. The mesh cells are rectangular parallelepipeds with faces
orthogonal to the specified axes of the Cartesian coordinate system. However, near the
boundary mesh cells are more complex. The near-boundary cells are portions of the
original parallelepiped cells that cut by the geometry boundary. The curved geometry
surface is approximated by set of polygons which vertexes are surface’s intersection
points with the cells’ edges. These flat polygons cut the original parallelepiped cells. Thus,
the resulting near-boundary cells are polyhedrons with both axis-oriented and arbitrary
oriented plane faces in this case (see Fig.4.1). The original parallelepiped cells containing
boundary are split into several control volumes that are referred to only one fluid or solid
medium. In the simplest case there are only two control volumes in the parallelepiped, one
is solid and another is fluid.

84
 Fig.4.1 Computational mesh near the solid/fluid interface.

centers of fluid control volumes

centers of solid control volumes

original curved geometry boundary

The rectangular computational domain is automatically constructed (may be changed

manually), so it encloses the solid body and has the boundary planes orthogonal to the
specified axes of the Cartesian coordinate system. Then, the computational mesh is
constructed in the following several stages.
First of all, a basic mesh is constructed. For that, the computational domain is divided into
slices by the basic mesh planes, which are evidently orthogonal to the axes of the
Cartesian coordinate system. The user can specify the number and spacing of these planes
along each of the axes. The so-called control planes whose position is specified by user
can be among these planes also. The basic mesh is determined solely by the computational
domain and does not depend on the solid/fluid interfaces.

Flow Simulation 2020 Technical Reference 85

Computational Mesh

Then, the basic mesh cells intersecting with the solid/fluid interface are split uniformly
into smaller cells in order to capture the solid/fluid interface with mesh cells of the
specified size (with respect to the basic mesh cells). The following procedure is employed:
each of the basic mesh cells intersecting with the solid/fluid interface is split uniformly
into 8 child cells; each of the child cells intersecting with the interface is in turn split into 8
cells of next level, and so on, until the specified cell size is attained.
At the next stage of meshing, the mesh obtained at the solid/fluid interface with the
previous procedure is refined (i.e. the cells are split further or probably merged) in
accordance with the solid/fluid interface curvature. The criterion to be satisfied is
established as follows: the maximum angle between the normals to the surface inside one
cell should not exceeds certain threshold, otherwise the cell is split into 8 cells.
Finally, the mesh obtained with these procedures is refined in the computational domain to
satisfy the so-called narrow channel criterion: for each cell lying at the solid/fluid
interface, the number of the mesh cells lying in the fluid region along the line normal to
the solid/fluid interface and starting from the center of this cell must not be less than the
criterion value. Otherwise each of the mesh cells on this line is split into 8 child cells.
As a result of all these meshing procedures, a locally refined rectangular computational
mesh is obtained. The final set of mesh cells that includes parallelepiped as well as more
complex polyhedrons is used for approximation of the governing equations.
Since all the above-mentioned meshing procedures are performed before the calculation,
the obtained mesh is unable to resolve all the solution features well. To overcome this
disadvantage, the computational mesh can be refined further at the specified moments
during the calculation in accordance with the solution spatial gradients (both in fluid and
in solid, see Online Help for details). As a result, in the low-gradient regions the cells are
merged, whereas in the high-gradient regions the cells are split. The moments of the
computational mesh refinement during the calculation are prescribed either automatically
or manually.
Note that in some cases the fine mesh may be unnecessary to give results of the required
accuracy.

Two-Scales Wall Functions Model

The Two-Scales Wall Functions (2SWF) model in Flow Simulation consist of two
approaches to coupling the boundary layer calculation with the main flow properties:
• To describe boundary layers on a fine mesh (the number of cells across a boundary
layer is 6 or greater) the “thick-boundary-layer” approach is used. In this approach
the calculation of parameters of laminar boundary layer is doing via Navier-Stokes
equations and for the turbulent boundary layer is performed by modification of well
known wall function approach. However instead of the classical approach where the
logarithmic velocity profile is used in EFD technology the full profile proposed by
Van Driest (Ref. 11) is used. All other assumptions are similar ones to the classical
wall functions approach.

86
 • The “thin-boundary-layer” approach is used to describe flow on a coarse mesh (the
number of cells across a boundary layer is 4 or less). In this approach the Prandtl
boundary layer equations already integrated along the normal to the wall from 0 (at
the wall) to the dynamic boundary layer thickness  are solved along a fluid
streamline covering the walls. If the boundary layer is laminar these equations are
solved with a method of successive approximations based on the Shvetz trial
functions technology (Ref. 10). If the boundary layer is turbulent or transitional
between laminar and turbulent, a generalization of this method to such boundary
layers employing the Van Driest hypothesis about the mixing length in turbulent
boundary layers is used (Ref. 11).
In intermediate cases, a compilation of the two above approaches is used, ensuring a
smooth transition between the two models as the mesh is refined, or as the boundary layer
thickens along a surface.
In addition to 2SWF model mentioned above, the “thin-channel” approach is used to
describe flow through narrow slots (including flow in pipes, plane and circular Couette
flows) on a coarse mesh (the number of cells across a slot is 7 or less). According to this
approach, the shear stress and heat flux near the wall are calculated by using the
approximations based on the experimental data.
By default, an appropriate boundary-layer approach is selected automatically according to
the computational mesh.
In the most cases all of these approaches provide the good accuracy, even on a coarse
mesh. However, in some cases when the appropriate boundary-layer approach is selected
automatically and the computational mesh is rather fine, the solution accuracy may fall
off. The reasons for the accuracy decrease are that the mesh is excessively fine to apply
the thin-boundary-layer, but it is insufficiently fine to resolve boundary layers and apply
the thick-boundary-layer. The further refinement of the computational mesh improves the
accuracy gradually.

Spatial Approximations
The cell-centered finite volume (FV) method is used to obtain conservative
approximations of the governing equations on locally refined grid that consists of
parallelepipeds and more complex polyhedrons near boundary. Following the FV method,
the governing equations are integrated over a control volume which is a mesh cell, and
then are approximated. All basic variables are referred to mass centers of control volumes.
These cell-centered values are used for approximations.

Flow Simulation 2020 Technical Reference 87

Spatial Approximations

The integral conservation laws may be represented in the form of the cell volume and
surface integral equation:


t 
Udv   F ds   Qdv (4.1)

are replaced by the discrete form


Uv    F  S  Qv (4.2)
t cell faces

The fluxes F are approximated in accordance with a type of the related faces. The faces
are classified and divided into two groups in accordance with their properties: if they are
axis-oriented and common for two adjacent control volumes or they are arbitrary oriented
boundary faces. Approximations are constructed in different way for these two cases.
Second order approximations are used for the faces that are common for two adjacent
control volumes. For convective fluxes the upwind approach is used. Nonlinear
approximations with limiters (Ref. 8) are used to provide monotonic discrete solutions.
For diffusive terms the central approximation is used. The approximations are treated
implicitly.
In Flow Simulation, especially consistent approximations for the convective terms, div
and grad operators are employed in order to derive a discrete problem that maintains the
fundamental properties of the parent differential problem in addition to the usual
properties of mass, momentum and energy conservation.

Spatial Approximations at the Solid/fluid Interface

On arbitrary oriented boundary faces, the fluxes are approximated in accordance with
boundary conditions and taking into account the curved geometry boundary. To calculate
accurately boundary fluxes on solid/fluid interface, sophisticated Flow Simulation
boundary layer model is used.
Conjugate heat transfer in fluid and solid is calculated at once as one problem, without
splitting into two explicitly related to each other problems.

88
Temporal Approximations
Time-implicit approximations of the continuity and convection/diffusion equations (for
momentum, temperature, etc.) are used together with an operator-splitting technique (Ref.
12, Ref. 13, and Ref. 14). This technique is used to efficiently resolve the problem of
pressure-velocity decoupling. Following the SIMPLE-like approach (Ref. 15), an elliptic
type discrete pressure equation is derived by algebraic transformations of the originally
derived discrete equations for mass and momentum, and taking into account the boundary
conditions for velocity.

Form of the Numerical Algorithm

Let index 'n' denotes the time-level, and '*' denotes intermediate values of the flow
parameters. The following numerical algorithm is employed to calculate flow parameters
on time-level (n+1) using known values on time-level (n):

U* - U n
t
 
 Ah U n , p n U*  S n ,
(4.3)

Lhp 
 
divh u*

1 *   n
,
(4.4)

t t t
* = (pn+p,T*,y*),

u n 1  u*  t  gradh p ,
(4.5)

pn 1  pn  p , (4.6)

T n1  T * ,  n1   * ,  n1   * , yn1  y* , (4.7)

 
 n 1   p n 1 , T n1 , y n 1 . (4.8)

Here U = (u,T, , , y)Tis the full set of basic variables excluding pressure p,
u=(u1,u2,u3)T is the velocity vector, y = (y1, y2, ..., yM)T is the vector of component
concentrations in fluid mixtures, and p = pn+1 - pn is an auxiliary variable that is called
a pressure correction. These parameters are discrete functions stored at cell centers. They
are to be calculated using the discrete equations (4.3)-(4.8) that approximate the governing
differential equations. In equations (4.3)-(4.8) Ah, divh, gradh and
Lh = divhgradh are discrete operators that approximate the corresponding differential
operators with second order accuracy.

Flow Simulation 2020 Technical Reference 89

Methods to Resolve Linear Algebraic Systems

Equation (4.3) corresponds to the first step of the algorithm when fully implicit discrete
convection/diffusion equations are solved to obtain the intermediate values of momentum
and the final values of turbulent parameters, temperature, and species concentrations.
The elliptic type equation (4.4) is used to calculate the pressure correction p. This
equation is defined in such a way that the final momentum field un+1 calculated from
(4.3) satisfies the discrete fully implicit continuity equation. Final values of the flow
parameters are defined by equations (4.5)-(4.8).

Methods to Resolve Linear Algebraic Systems

Iterative Methods for Nonsymmetrical Problems

To solve the asymmetric systems of linear equations that arise from approximations of
momentum, temperature and species equations (4.3), a preconditioned generalized
conjugate gradient method (Ref. 16) is used. Incomplete LU factorization is used for
preconditioning.

Iterative Methods for Symmetric Problems

To solve symmetric algebraic problem for pressure-correction (4.4), an original
double-preconditioned iterative procedure is used. It is based on a specially developed
multigrid method (Ref. 17).

Multigrid Method
The multigrid method is a convenient acceleration technique which can greatly decrease
the solution time. Basic features of the multigrid algorithm are as follows. Based on the
given mesh, a sequence of grids (grid levels) are constructed, with a decreasing number of
nodes. On every such grid, the residual of the associated system of algebraic equations is
restricted onto the coarser grid level, forming the right hand side of the system on that
grid. When the solution on the coarse grid is computed, it is interpolated to the finer grid
and used there as a correction to the result of the previous iteration. After that, several
smoothing iterations are performed. This procedure is applied repeatedly on every grid
level until the corresponding iteration meets the stopping criteria.
The coefficients of the linear algebraic systems associated with the grid are computed
once and stored.

90
Nested Iterations
Employment of an operator-splitting technique is one of basic numerical approaches in
Flow Simulation. This means that not the whole set of non-linear governing equations
with non-linear boundary conditions are solved simultaneously. Instead of this, equation
by equation is solved in some linearized form with linearized boundary conditions.
In order to provide more robust and accurate numerical method, nested iterations are
introduced on each time step. These iterations allow treating mutual influence of flow
parameters and non-linear terms in equations and boundary conditions.
These iterations are incorporated within each time step and repeated until a solver
convergence is obtained or until a Maximum number of nested iterations is reached or
each conservation law convergence condition is satisfied.
Let the governing equations are treated by using the nested iterations on the time span
[t,t+t]. Let k be a number of the nested iteration, the k-th iteration is finished and the
(k+1)-th iteration is starting.

Criterion of mass equation convergence

Let 0 be the density at time t and k is calculated approximation for the density at time
k
t+t. While performing the (k+1)-th nested iteration, the residual error rmass for the mass
equation is calculated and defined as follows:
k
rmass k
= rmass,cv
1
= r

k
mass,cv , (4.9)

where 
means summing over all control volumes in the domain (sub-domain) Ω and
the residual error at the control volume cv is:

r k

 k
 0 
V cv   j kfcv ,
(4.10)
mass , cv
t f cv

where jmass, fcv is the mass flux at the control volume face fcv .

k k
To make decision if rmass is small enough, a reference residual value rmass , ref is defined
as:

, ref   rmass,cv ,ref   jmass, f cv  2 jmass,cv

k k k
rmass (4.11)
  f cv

Flow Simulation 2020 Technical Reference 91

Nested Iterations

So a sum of absolute values of inlet and outlet fluxes at control volume faces can be
considered as doubled mass flux through the control volume cv.
At the end of the (k+1)-th nested iteration, a stopping criteria are checked. For the mass
equation it is:
k
rmass   mass  rmass
k
,ref ,
(4.12)

and the stopping criteria (4.12) can be treated as:

k
rmass ,cv
k
  mass (4.13)

2 jmass ,cv

Criterion of energy equation convergence

k
The residual error renergy for the energy equation is defined as follows:
k
renergy k
= renergy,cv =  renergy,cv
k
, (4.14)
1


where the residual error at the control volume cv is:

r k
 0 H k
H0   
Vcv  LkH ,Conv  Diff H k - f Hk
(4.15)
energy ,cv
t
For the energy equation, the calculation of reference residual value is based on the use of
approximation coefficients. The discrete operator

  c
LkH ,Conv  Diff H k  j  H kj (4.16)
j cv

approximates convection and diffusion terms of the equation at the control volume cv.
Here ωcv is a stencil of the operator L at the control volume cv, c0 is a coefficient for H at

the control volume under consideration and c j  j0  are coefficients for other control
volumes from the stencil. Note, that the coefficients of the approximation are:

 
 h 

c j ~ u n  S (4.17)
 hcv  fcv
 
 2 
where un is the velocity component normal to the face and hcv is the characteristic size of
the control volume.

92
 k
The explicit member f H contains approximations of other terms of the equation and the

heat source at the control volume face equals q  Vcv .

k
The reference residual value renergy , ref is introduces as:

  0  k  0 H 0  (4.18)
,ref   renergy,cv,ref   max c0 
k
renergy k   Vcv , max c j   H  f Hk  Vcv
t j , j 0  t 
     

hcv hcv2
In common case of t   , the first summand is the reference value of the
un H 
energy flux through the control volume cv. The second summand is the reference value of
the energy flux in/out of the control volume due to the heat source q  Vcv .

,cv ,ref 
  
k
In case of non-large Δt: renergy t 0 .

k
Thus, renergy,cv,ref is treated as the reference value of the energy flux jenergy ,cv through
the control volume. Finally, the convergence condition for the energy equation is:
k
renergy   energy  renergy
k
,ref ,
(4.19)

and the stopping criteria (4.19) can be treated as:

k
renergy,cv
k
  energy (4.20)
jenergy,cv

Criterion of momentum equation convergence

k
In a way that is quite similar to the energy equation, the residual error rmomentum for the
momentum equation is defined as follows:
k
rmomentum k
= rmomentum,c v 1 =  r

k
momentum 0 ,cv
 rmomentum
k
1 ,cv
 rmomentum
k
2 ,cv
 (4.21)
where the residual error for i-th component of momentum at the control volume cv is:

m k
 mi0    (4.22)
, cv  Vb  Lmi ,Conv  Diff mik - f mki , i  0,1,2
k i
rmomentum
t

Flow Simulation 2020 Technical Reference 93

Nested Iterations

Here:

Lkmi ,Conv  Diff mik    c ij  mijk (4.23)

j cv

where:

 
   t 
cij ~  un  S (4.24)
 hcv  fcv
 
 2 
k
where the member f contains approximations of other terms including a pressure
mi

gradient and forces.

k
The reference residual value rmomentum , ref is defined as:
k
rmomentum , ref =  r
k
momentum 0 ,cv , ref  rmomentum
k
1 ,cv , ref
 rmomentum
k
2 ,cv , ref
 (4.25)

where
  V  m i0 
rk
momentum i ,cv , ref    max  c i 0  b , max c ij   m  f mi 
k k
Vb 

(4.26)
    t j , j  0  t 

k
Thus, rmomentum ,cv ,ref is treated as the reference value of the momentum flux

jmomentum, cv through the control volume. Finally, the convergence condition for the
momentum equation is:
k
rmomentum   momentum  rmomentum
k
, ref , (4.27)

and the stopping criteria (4.27) can be treated as:

k
rmomentum ,cv
k
  momentum (4.28)
jmomentum ,cv

94
Notes on the residuals behavior
In a particular case of complex flow it is not evident what the time scale is for unsteady
processes. The residuals do not go down if the time step in calculations is comparable or
greater than the time scale of unsteady processes.
In any case, any behavior of residuals does not guarantee accuracy of a discrete solution
due to the fact that there are still approximation errors.
Note, that residuals do not tell whether the solution is accurate. Low residuals do not
automatically mean an accurate solution, and high residuals do not automatically mean a
wrong solution. Only investigation of physical parameters convergence on a sequence of
decreasing time steps tells about the actual accuracy of the discrete solution.

Flow Simulation 2020 Technical Reference 95

Nested Iterations

96
 5
Mesh Settings

Introduction
This chapter provides the fundamentals of working with the Flow Simulation
computational mesh, describes the mesh generation procedure and explains the use of
parameters governing both automatically and manually controlled meshes.
First, let us introduce a set of definitions.
Flow Simulation considers the real model created in SOLIDWORKS and automatically
generates a rectangular computational mesh in the Computational Domain distinguishing
the fluid and solid domains.
The corresponding Computational Domain is generated in the form of a rectangular
parallelepiped enclosing the model for both the 3D analysis and 2D analysis. Its
boundaries are parallel to the Global Coordinate System planes. For External flows, the
computational domain’s boundary planes are automatically distanced from the model. For
Internal flows, the computational domain’s boundary planes automatically envelop either
the entire model, if Heat Conduction in Solids is considered, or if Heat Conduction in
Solids is not considered, the model’s flow passage only.

In the mesh generation process, the computational domain is divided into uniform
rectangular parallelepiped-shaped cells, which form a so-called basic mesh. Then, using
information about the model geometry, the specified boundary conditions and goals Flow
Simulation further constructs the mesh by means of various refinements, i.e. splitting of
the basic mesh cells into smaller cells, thus better representing the model and fluid
regions. The mesh, which the calculation starts from, so-called initial mesh, is fully
defined by the generated basic mesh and the refinement settings.

Flow Simulation 2020 Technical Reference 97

Types of Cells

Each refinement has its criterion and level. The refinement criterion denotes which cells
have to be split, and the refinement level denotes the smallest size, which the cells can be
split to. Regardless of the refinement considered, the smallest cell size is always defined
with respect to the basic mesh cell size so the constructed basic mesh is of great
importance for the resulting computational mesh.
The main types of refinements are:
 Cell Type Refinement

 Channel Refinement

 Advanced Refinements:
• Small Solid Features Refinement
• Curvature Refinement
• Tolerance Refinement
In addition, the following type of refinements can be invoked locally:
 Equidistant Refinement
During the calculation, the initial mesh can be refined further using the
 Solution-Adaptive Refinement (see "Refinement of the Computational Mesh During
Calculation" on page 127).

Though it depends on a refinement which criterion or level is available for user control,
we will consider all of them (except for the Solution-Adaptive Refinement) to give you a
comprehensive understanding of how the Flow Simulation meshing works.

Types of Cells
Any Flow Simulation calculation is performed in a rectangular parallelepiped-shaped
computational domain which boundaries are orthogonal to the axes of the Cartesian
Global Coordinate System. A computational mesh splits the computational domain with a
set of planes orthogonal to the Cartesian Global Coordinate System's axes to form
rectangular parallelepipeds called cells. The original parallelepiped cells containing
boundary are split into several parts that are referred to only one fluid or solid medium.
The resulting computational mesh consists of cells of the following types:
• Fluid cells are the cells located entirely in the fluid.
• Solid cells are the cells located entirely in the solid.
• Solid-fluid boundary cells are the cells which are partly in the solid and partly in the
fluid.
As an illustration let us look at the generated computational mesh (Fig. 5.1).

98
 Fig. 5.1 The computational mesh cells of different types.
Fluid cells Solid-Fluid boundary cells

Solid cells

Mesh Construction Stages

The mesh is constructed in the following steps:
 Construction of the basic mesh taking into account the user-specified Control Planes
and the respective values of cells number (or cells size) and cell size ratios.
 Refinement of all fluid and solid mesh cells up to the user-specified level.

 Resolving of the interface between substances to resolve the relatively small solid
features and solid/solid interface, tolerance and curvature refinement of the mesh at a
solid/fluid boundaries to resolve the interface curvature (e.g. small-radius surfaces of
revolution, etc).
 Channels refinement, that is the refinement of the mesh in narrow channels taking into
account the respective user-specified settings.
Note: The mesh at this stage is called the initial mesh. The initial mesh implies the
complete basic mesh with the resolution of the solid/fluid (as well as solid/solid)
interface by the small solid features refinements, the curvature and channels
refinement also taking into account the local mesh settings.
After each of these stages is passed, the number of cells is increased to some extent.
Tip: If you switch on or off heat conduction in solids, or add/modify solid materials, you
should rebuild the mesh.
In Flow Simulation you can control the following parameters and options which govern
the computational mesh:
1 Nx, the number of basic mesh cells (zero level cells) along the X axis of the Global
Coordinate System. 1 ≤ Nx ≤ 100 000

Flow Simulation 2020 Technical Reference 99

Basic Mesh

2 Ny, the number of basic mesh cells (zero level cells) along the Y axis of the Global
Coordinate System. 1 ≤ Ny ≤ 100 000.
3 Nz, the number of basic mesh cells (zero level cells) along the Z axis of the Global
Coordinate System. 1 ≤ Nz ≤ 100 000.
4 Control planes. By adding and relocating them you can contract and/or stretch the
basic mesh in the specified directions and regions. In case of the 3D analysis, six
control planes coincident with the computational domain's boundaries are always
present in any project.
5 Channels refinement: Characteristic number of cells across a channel (Nch), Maximum
Channels refinement level (0 ≤ Lch ≤ 9), The minimum and maximum height of
channels (Hmin and Hmax) to be refined.
6 Small solid features refinement level (0 ≤ Lssf ≤ 9).
7 Curvature refinement level (0 ≤ Lcur ≤ 9).

8 Curvature refinement criterion (0 ≤ Ccur ≤ ).

9 Tolerance refinement level (0 ≤ Ltol ≤ 9).
10 Tolerance refinement criterion (0 ≤ C tol).
11 Display refinement level. The Display Refinement Level option is intended to inspect
the probable refined initial mesh but not to change the specified refinement settings
that will be applied. The specified refinement level is measured from the basic mesh.
These options are described in more details below.

Basic Mesh
The basic mesh is a mesh of zero level cells. It is constructed for the whole computational
domain at the beginning of the meshing process and formed by dividing the computational
domain into slices by parallel planes which are orthogonal to the Global Coordinate
System’s axes. The basic mesh can be defined either by the number of basic mesh cells or
by the average size of the basic mesh cells along each coordinate direction.
In case of 2D calculation (i.e. if you select the 2D plane flow option in the Computational
Domain dialog box) only one basic mesh cell is generated automatically along the
eliminated direction. By default Flow Simulation constructs each cell as close to cubic
shape as possible.
Note: The number of basic mesh cells could be one or two less than the user-defined
number (Nx, Ny, Nz). There is no limitation on a cell oblongness or aspect ratio, but you
should carefully check the calculation results in all cases for the absence of too oblong or
stretched cells.

100
Control Planes
The Control Planes option is a powerful tool for creating an optimal computational mesh,
and the user should certainly become acquainted with this tool if he is interested in
optimal meshes resulting in higher accuracy and decreasing the CPU time and required
computer memory. Control planes allow you to contract the basic mesh locally to resolve a
particular region of interest by a denser mesh and stretch the basic mesh to avoid
excessively dense meshes in other regions.
All of the control planes are orthogonal to the Global Coordinate System’s axes. The
Computational domain boundary planes (X min, X max, Y min, Y max, Z min, Z max) are
among the Control Planes by default.

Contracting the Basic Mesh

Using control planes you may contract the basic mesh in the regions of interest. To do this,
you need to set control planes surrounding the region and assign the proper Ratio values to
the respective intervals. The cell sizes on the interval are changed gradually so that the
proportion between the first and the last cells of the interval is close (but not necessarily
equal) to the entered Ratio value. Negative values of the ratio correspond to the reverse
order of cell size increase. Alternatively, you may explicitly set the Number of cells for
each interval, in which case the Ratio value becomes mandatory. For example, assume
that there are two control planes Plane1 and Plane2 (see Fig. 5.2) and the ratio on the
interval between them is set to 2. Then the basic mesh cells adjacent to the Plane1 will be
approximately two times longer than the basic mesh cells adjacent to the Plane2.

Flow Simulation 2020 Technical Reference 101

Basic Mesh

Fig. 5.2 Specifying custom control planes.

Plane 4 Default control plane

Interval 3:
number of cells=3 (automatic)
Plane 3 ratio=1.5
Custom control plane

Interval 2:
number of cells=12 (manual)
ratio=1

Plane 2 Custom control plane

Interval 1:
number of cells=12 (automatic)
ratio=2

Plane 1 Default control plane

Use of control planes is especially recommended for external analyses, where the
computational domain may be substantially larger than the model.
Fig. 5.3 Uniform mesh, default control planes Fig. 5.4 Two custom control planes in each
only. direction.

In the Fig. 5.4 two custom control planes are set around the body in each coordinate
direction with the ratio set to 5 and -5, respectively.

102
Mesh Refinement
Refinement is a process of splitting a rectangular computational mesh cell into eight cells
by three orthogonal planes that divide the cell's edges in halves. The non-split initial cells
that compose the basic mesh are called basic cells or zero level cells. Cells obtained by the
first splitting of the basic cells are called first level cells, the next splitting produces
second level cells, and so on. The maximum level of splitting is nine. A ninth level cell is
89 times smaller in volume than the basic cell.
Tip: During the solution-adaptive meshing the cells can be refined and merged. See
"Refinement of the Computational Mesh During Calculation" on page 127.

The following Cell Mating rule is applied to the processes of refinement: the levels of two
neighboring cells (i.e. cells having a common face) can only be the same or differ by one,
so that, say, a fifth level cell can have only neighboring cells of fourth, fifth, or sixth level.
This rule has the highest priority.

The fourth-level red cells appearing Fig. 5.5 Fluid cell refinement due to the Cell Mating
after resolving the cog cause the rule.
neighboring cells to be split up to
third level (yellow cells), that, in
turn, causes the subsequent
refinement producing second level
cells (green cells) and first level cells
(blue cells). The white zero level cell
(basic mesh cell) remains unsplit
since it borders on first level cells
only, thus satisfying the rule.

Note: The Cell Mating rule is strict and has higher priority than the other cell operations.
The rule is also enforced for the cells that are entirely in a solid.

Refining Cells by Type

The refinement level of cells of a specific type (the combination of fluid cells, and solid
cells) denotes the minimum level to which the corresponding cells must be split if it
doesn’t contradict the Cell Mating rule.
Note: The minimum level means the lower bound to which it is obligatory to split cells,
though the cells can be split further if it is required to satisfy the other criteria such as
Small solid features refinement, Curvature refinement, Channels refinement or
Tolerance refinement.

Flow Simulation 2020 Technical Reference 103

Mesh Refinement

Advanced Refinement at Interfaces Between Substances

Different interface types (solid/fluid or solid/solid) are checked on different refinement
criteria, namely:
• small solid features criterion for solid/fluid, solid/solid, solid/porous and
porous/fluid interfaces;
• curvature refinement criterion for solid/fluid, solid/porous and porous/fluid
interfaces;
• tolerance refinement criterion for solid/fluid and solid/porous interfaces;
• narrow channel refinement criterion for solid/fluid and solid/porous interfaces.
Whereas the specified refinement levels are equally applied to any interface type.
Tip: The solid/solid interface is the interfaces between different solid materials. If two
solids of the same material contact to each other, they are combined into one solid.
To resolve the interface between substances, the refined mesh is constructed in the
following steps:
1 Analyze triangulation in each basic mesh cell at the interfaces between different
substances (such as solid/fluid and solid/solid interfaces) in order to find the maximum
angle between normals to the triangles which compose the interface within the cell.
Depending on the maximum angle found, the decision whether to split the cell or not is
made in accordance with the specified Curvature Criterion (Ccur), Tolerance
Criterion (Ctol) and the Small Solid Feature Refinement criterion (Cssf), also called
the 120-degree criterion.
2 Analyze the distances between the opposite walls of each flow passage in the direction
normal to wall. Depending on the distances found, the decision whether to split the cell
or not is made in accordance with the specified Characteristic Number of Cells
Across Channel (Nch), and Minimum (Maximum) Height of Channel (Hmin and
Hmax).
3 The cells that marked to be split as described in items 1-2, are split if the specified
Small Solid Feature Refinement Level (Lssf), Curvature Level (Lcur), Tolerance
Level (Ltol) and Maximum Channel Refinement Level (Lch) are not reached.

Tip: If a cell belongs to a local initial mesh area, then the corresponding local
refinement levels will be applied (see "Local Mesh Settings" on page 120).
4 If a basic mesh cell is split, the resulting child cells are analyzed as described in items
1-3, and split further, if necessary. The cell splitting will proceed until the interface
resolution satisfies the Small Solid Feature Refinement criterion (Cssf), the specified
Ccur, Ctol, Nch, Hmin and Hmax or the corresponding level of splitting reaches its
specified value.
Tip: The specified levels of splitting denote the maximum admissible splitting, i.e. they
show to which level a basic mesh cell can be split if it is required for resolving the
solid/fluid interface within the cell.

104
 5 The operations 1 to 4 are applied for the next basic mesh cell and so on, taking into
account the following Cell Mating rule: two neighboring cells can be only cells which
levels are similar or differ by one. This rule has the highest priority as it is necessary
for simplifying numerical algorithm in solver.

Small Solid Features Refinement

The procedure of resolving small solid features refines only the cells where the solid/fluid
and solid/solid interface curvature is too high: the maximum angle between the normals to
a solid surface inside the cell is strictly greater than 120, i.e. the solid surface has a
protrusion within the cell. Such cells are split until the Small Solid Feature Refinement
Level (Lssf) of splitting mesh cells is achieved.

To make this Small Solid Feature Refinement criterion (Cssf), also called the 120-degree
criterion, easier to understand, let us consider simple small solid features of planar faces
only. The normal to triangles that form the planar face is normal to the planar face too.
Therefore, instead of considering the normals to the triangles we can consider normals to
faces, or better the angle between faces.
Fig. 5.6 Lssf = 1, Lcur = 0, Lch = 0.

In Fig. 5.6 the cells with the cogs of 150 and 60 degrees were not split by the small solid
features refinement because the maximum angles between the faces (i.e. between normals
to the triangles enclosed by the cell) are 30 and 120, respectively. If the angle between
the normals becomes greater than 120 (121 for the 59-cog) then the cell is split. The
cell with the square spike surely has to be split because the lateral faces of the spike have
their normals at the angle of 180, thus satisfying the 120-degree criterion.
Note: Rectangular corners (like in the rightmost cell) do not satisfy the criterion and
therefore will not be resolved by the small solid features refinement.
Remember that if the Channel refinement is enabled, the maximum level to which the
small solid features refinement can split the cells is set as the maximum level from the
specified Small Solid Feature Refinement Level (Lssf) and Maximum Channel
Refinement Level (Lch). In other words, if the Channel refinement is enabled, the Lssf
has no effect if it is smaller than the Lch.

Flow Simulation 2020 Technical Reference 105

Mesh Refinement

Fig. 5.7 Lssf = 0, Lcur = 0, Lch = 1.

From Fig. 5.7 it is clear that the cells are split by the 120-degree criterion up to the first
level, as defined by the narrow channel refinement level.
For the information about how the Lch influences the narrow channel refinement see
"Channel Refinement" on page 110.

Curvature Refinement
The Curvature refinement level is the maximum level to which the cells will be split
during refinement of the computational mesh until the curvature of the solid/fluid
interface within the cell becomes lower than the specified Curvature criterion (Ccur).
The curvature refinement procedure has the following stages:
1 Each solid surface is triangulated: Flow Simulation gets triangles that make up the
surfaces.
Note: The performance settings do not govern the triangulation performance.
2 A local (for each cell) interface curvature is determined as the maximum angle
between the normals to the triangles within the cell.
3 If this angle exceeds the specified Ccur, and the curvature refinement level is not
reached then the cell is split.
The curvature refinement works in the same manner as the small solid features refinement
with the difference that the critical angle between the normals can be specified by the user
(in radians) as Curvature refinement criterion (Ccur). Here, the smaller the criterion, the
better resolution of the solid curvature. To give more precise and descriptive explanation,
the following table presents several Ccur values together with the corresponding angles
between normals and the angles between planar faces.
Table 1: Influence of the curvature criterion on the solid curvature resolution.

Curvature
0.3491 0.5062 0.5411 0.6982 1.0472 1.5708 2.0944 3.1416
criterion, rad

' between 20 29 31 40 60 90 120 180

normals, [degrees]

106
 Curvature
0.3491 0.5062 0.5411 0.6982 1.0472 1.5708 2.0944 3.1416
criterion, rad

 between faces, 160 151 149 140 120 90 60 0

[degrees]

The table states that if the Ccur is equal to 0.5062 rad, then all the cells where the angle
between normals to the surface-forming triangles is more than 29 degrees will be split.
Fig. 5.8 Lssf = 0, Lch = 0, Lcur = 1, Fig. 5.9 Lssf = 0, Lch = 0, Lcur = 1,
Ccur = 0.5411 rad (' = 31°,  = 149°) Ccur = 0.5062 rad (' = 29°,  = 151°)

You can see that the curvature criterion set to 0.5062 rad splits the cells with the
151-degrees cog.
Note: The curvature refinement works exactly as the small solid features refinement when
the curvature criterion is equal to 2.0944 rad (120°).
However, the default curvature criterion values are small enough to resolve obtuse angles
and curvature well. Increasing the curvature criterion is reasonable if you want to avoid
superfluous refinement but it is recommended that you try different criteria to find the
most appropriate one.
The curvature refinement is a powerful tool, so that the competent usage of it allows you
to obtain proper and optimal computational mesh. Look at the following illustrations to
the curvature refinement by the example of a sphere.

Flow Simulation 2020 Technical Reference 107

Mesh Refinement

You can see that the curvature criterion set to 0.317 rad and 0.1 rad splits the cells up to
the first level only. In Fig. 5.11 and Fig. 5.12 the cells with the cogs of 162 were split. In
Fig. 5.13 the cells with the cogs of 174 were split.

Fig. 5.10Lcur = 0; Fig. 5.11 Lcur = 1; Ccur = 0.317 rad (' = 18°);
Total number of cells is 64. Total number of cells is 120.

Fig. 5.12Lcur = 2; Ccur = 0.317 rad (' = 18°); Fig. 5.13Lcur = 2; Ccur = 0.1 rad (' = 6°);
Total number of cells is 120. Total number of cells is 148.

Small Solid Feature Refinement Level or Curvature Level

Why is it necessary to have two criteria? As you can see, the curvature refinement has
higher priority than the small solid features refinement if the curvature criterion is smaller
than 120°.
Note: The Flow Simulation-specified values of the curvature criterion are always smaller
than 120°. In other words, if you did not set the Ccur greater than 120° and if the Lssf and
Lch are smaller than the Lcur, then the small solid feature refinement would be idle.

108
 Nevertheless, the advantage of the small solid features refinement is that being sensitive to
relatively small geometry features it does not “notice” the large-scale curvatures, thus
avoiding refinements in the entire computational domain but resolving only the areas of
small features. At the same time, the curvature refinement can be used to resolve the
large-scale curvatures. So both the refinements have their own coverage providing a
flexible tool for creating an optimal mesh.

Tolerance Refinement
Tolerance refinement allows you to control how well (with what tolerance) mesh
polygons approximate the real interface. The tolerance refinement may affect the same
cells that were affected by the small solid features refinement and the curvature
refinement. It resolves the interface's curvature more effectively than the small solid
features refinement, and, in contrast to the curvature refinement, discerns small and large
features of equal curvature, thus avoiding refinements in regions of less importance (see
images below).
Fig. 5.14Lcur = 3; Ccur = 0.1 (' = 6°,  = 174°).

Fig. 5.15Ltol = 3; Ctol = 0.1 mm.

Flow Simulation 2020 Technical Reference 109

Mesh Refinement

Any surface is approximated by a set of polygons which vertices are surface's intersection
points with the cells' edges. This approach accurately represents flat faces though
curvature surfaces are approximated with some deviations (e.g. a circle can be
approximated by a polygon). The tolerance refinement criterion controls the precision of
this approximation. A cell will be split if the distance (h) between the outermost
interface's point within the cell and the polygon approximating this interface is larger than
the specified criterion value.
Fig. 5.16Tolerance refinement vs. Small Solid Features and Curvature refinements
(refinement level is set to 1).

Tolerance Refinement Small Solid Feature Refinement

(refinement occurs regardless of the
Ctol = 0.1 mm Ctol = 0.08 mm feature’s size)

Curvature Refinement
Ctol = 0.1 mm Ctol = 0.03 mm (refinement occurs regardless of the
Refines cells only if the solid part cut by curvature only)
the polygon is large enough (h > 0.1 mm)

Channel Refinement
After the primary mesh has been created, the channel refinement is put in action. The
channel refinement is applied to each flow passage within the computational domain (or a
region, in case that local mesh settings are specified) unless you specify for Flow
Simulation to ignore passages of a specified height.
The Narrow Channels term is conventional and used for the definition of the flow
passages of the model in the direction normal to the solid/fluid interface. Regardless of the
real solid curvature, the mesh approximation is that the solid boundary is always
represented by a set of flat elements, which nodes are the points where the model
intersects with the cell edges. Thus, whatever the model geometry, there is always a flat
element within a solid-fluid boundary cell and the normal to this element denotes the
direction normal to the solid/fluid interface for this solid-fluid boundary cell.

110
 The basic concept of narrow channel refinement is to resolve the narrow channels with a
sufficient number of cells to provide a reasonable level of solution accuracy. It is
especially important to have narrow channels resolved in analyses of low Reynolds
numbers or analyses with long channels, i.e. in such analyses where the boundary layer
thickness becomes comparable to the size of the solid-fluid boundary cells where the layer
is developed.
There are two ways to specify the channels refinement:
• Number of Cells mode governs the procedure of mesh refining in the model’s
narrow channels by specifying the number of mesh cells across model’s flow
passages and restricting the refinement level;
• Refinement Level mode allows to define an uniform mesh across each model’s
flow passage by specifying the refinement level as a tabular dependency on the
channel height.
The narrow channel settings available in the Number of Cells mode are the following:
• Maximum Channel Refinement Level – the maximum level of cells refinement in
narrow channels with respect to the basic mesh cell.
• Characteristic Number of Cells Across Channel – the number of cells (including
cells lying at the solid/fluid interface) that Flow Simulation will attempt to set
across the model flow passages in the direction normal to the solid/fluid interface. If
possible, the number of cells across narrow channels will be equal to the specified
characteristic number, otherwise it will be as close to it as possible. The
Characteristic Number of Cells Across Channel (Nch) and the Maximum
Channel Refinement Level (Lch) both influence the mesh in narrow channels in the
following manner: the basic mesh in narrow channels will be split to have Nch
number per channel, if the resulting cells satisfy the specified Lch. In other words,
whatever the specified Nch, the smallest possible cell in a narrow channel is 8L times
smaller in volume (or 2L times smaller in each linear dimension) than the basic
mesh cell. This is necessary to avoid undesirable mesh splitting in very fine
channels that may cause the number of cells to increase to an unreasonable value.
• Minimum Height of Channel, Maximum Height of Channel – the minimum and
maximum bounds for the height outside of which a flow passage will not be
considered as a narrow channel and thus will not be refined by the narrow channel
resolution procedure.
For example, if you specify the minimum and maximum height of narrow channels, the
cells will be split only in those fluid regions where the distance between the opposite walls
of the flow passage in the direction normal to wall lies between the specified minimum
and maximum heights.
The Number of Cells channel refinement mode operates as follows:
1 For each solid-fluid boundary cell Flow Simulation calculates the “local” narrow
channel width as the distance between this solid-fluid boundary cell and the next

Flow Simulation 2020 Technical Reference 111

Mesh Refinement

solid-fluid boundary cell found on the line normal to the solid/fluid interface of this
cell (i.e. normal to the flat surface element located in the cell).
Tip: If the line normal to the solid/fluid interface crosses a local initial mesh area, then
the corresponding local narrow channel refinement settings is applied to the cells in
this direction.
2 If the distance value falls within the range defined by the Minimum height of channel
(Hmin) and Maximum height of channel (Hmax) options, the number of cells per this
interval is calculated including both cells lying at the solid/fluid interface and taking
into account which portion of each cell is in fluid.
3 More precisely, the number of cells across the channel (i.e. on the interval between the
two solid-fluid boundary cells) is calculated as N = Nf + np1 + np2, where Nf is the
number of fluid cells on the interval, and np1 and np2 are the fluid portions of the both
solid-fluid boundary cells. This value is compared with the specified Characteristic
number of cells across a channel (Nch). If N is less than the specified Nch then the cells
on this interval are split. For example, on Fig. 5.17 Nf = 2, np1 = np2 = 0.4, and
N = 2+0.4+0.4 = 2.8 which is less than the criterion Nch = 3. On Fig. 5.18 the
solid-fluid boundary cells are split, so that the fluid portions of the newly-formed
solid-fluid boundary cells are np1 = np2 = 0.9, and the criterion is satisfied (N > Nch).
Fig. 5.17Lch = 2; Nch = 3; Fig. 5.18Lch = 3; Nch = 3;
N = 2.8 < Nch N = 3.8 > Nch

Note: Like in the other refinements, the Maximum channel refinement level (Lch)
denotes the maximum level to which the cells can be split to satisfy the Nch criterion.
The Lch has higher priority than the Nch, so the refinement will proceed until the Nch
criterion is satisfied or all the cells reach the Lch.
The channel refinement is symmetrical with respect to the midpoint of the interval and
proceeds from the both ending solid-fluid boundary cells towards the midpoint. Since
the actual number of cells across narrow channels can be greater by 1 than the specified
characteristic number.
Fig. 5.19Nch = 5; Lch = 1 Fig. 5.20Nch = 5; Lch = 3

112
 In Fig. 5.19, the specified Characteristic number of cells across a channel (Nch) is 5
but only two cells were generated since the Maximum channel refinement level (Lch)
of one allows only basic mesh cells and first-level cells to be generated.
In Fig. 5.20, the specified Maximum channel refinement level (Lch) is high enough to
allow 5 cells to be placed across the channel, but there are 6 cells across the channel
due to the symmetry requirement of the channel refinement.
If the channels height is variable along the channel length, pay attention to the Nch
criterion and make sure that the proper value of the Nch criterion is satisfied along the
entire channel.
Fig. 5.21Hmax = 60 mm; Lch = 4; Nch = 10 Fig. 5.22Hmax = 60 mm; Lch = 4; Nch = 25

In Fig. 5.21, the specified Characteristic number of cells across a channel (Nch) is not
enough to reach the specified Maximum channel refinement level (Lch) across the
entire channel. Only the cells located in the narrowest end of the channel reach the
specified Lch.
In Fig. 5.22, the specified Characteristic number of cells across a channel (Nch) is high
enough to reach the specified Lch across the entire channel.
Alternatively, to generate an uniform mesh across a channel, you can specify the
refinement level as a tabular dependency on the channel height by using the
Refinement Level mode.
Fig. 5.23Number of Cells (Hmax = 80 mm): Fig. 5.24Refinement Level:
Lch = 3; Nch = 10 H ≤ 81 mm: Lch = 3

Flow Simulation 2020 Technical Reference 113

Mesh Refinement

4 Next, for all the fluid cells within the entire computational domain the following Fluid
Cell Leveling procedure is applied: if a fluid cell is located between two cells of higher
level, it is split to be equalized with the level of neighboring smaller cells.
Although the settings that produce an optimal mesh depends on a particular task, here are
some ’rule-of-thumb’ recommendations for narrow channel settings:
1 Set the number of cells across a channel to a minimum of 5.
2 Use the minimum and maximum heights of narrow channels to concentrate on the
regions of interest.
3 If possible, avoid setting high values for the narrow channels refinement level, since it
may cause a significant increase in the number of cells where it is not necessary.
Fig. 5.25Lssf = 3; Channels refinement is disabled.

114
 Fig. 5.26 Lssf = 3; Channels refinement is on: Nch = 5, Lch = 2.

Fig. 5.27Lssf = 3; Channels refinement is on: Nch = 5, Lch = 5.

Thin walls resolution

In contrast to the narrow channels, thin walls can be resolved without the mesh refinement
inside the wall, since the both sides of the thin wall may reside in the same cell. Therefore,
the amount of cells needed to resolve a thin wall is generally lower than the number of
cells needed to properly resolve a channel of the same width. See Fig. 5.28 - 5.30
illustrating the thin walls resolution technology and its limitations.

Flow Simulation 2020 Technical Reference 115

Thin walls resolution

Fig. 5.28One mesh cell can contain more than one fluid and/or solid volume; during calculation each volume
has an individual set of parameters depending on its type (fluid or solid).

Solid 2 Fluid 1

Solid 1
Fluid 2

Fig. 5.29If the wall thickness is greater than the basic mesh cell's size across the wall or if the wall creates only
one fluid volume in the cell, then the opposite sides of the wall will not lay within the same cell. Such walls are
resolved with two or more cells across.

Fig. 5.30The edges of thin walls ending within a mesh cell may be trimmed in certain cases. These mesh cells
are called Trimmed cells.

Model geometry Meshed geometry

Trimmed edge

Trimmed cell

116
Automatic Mesh Setting
As noted above, Flow Simulation mesh generator uses a system of refinement criteria,
each with their own values and refinement levels. However, in order to further simplify
the process of defining the mesh a method to automatically define the refinement criteria,
values and levels for a case Automatic Parameters Definition (APD) technology has been
implemented.
The main idea of the APD technology is to define both the basic mesh and refinement
settings from the physical definition of the model’s boundary conditions, etc. together
with the most general information about the geometry.
APD uses the following input parameters:
1 Problem type (Internal\External, Compressible\Incompressible, 3D\2D)
2 Minimum gap size (hgap) and minimum wall thickness (hwall)
3 Level of initial mesh (1 ≤ Lini ≤ 7)
4 Computational model bounding box (Bmodel)
5 Fluid region bounding box (Bfluid)
6 Symmetry settings applied to the computational domain boundaries
7 Middle size of the computational model (Hmean)
8 Wind direction for external flows (ewind).

The Level of initial mesh (Lini) is specified by user. The default values of Minimum gap
size (hgap) and Minimum wall thickness (hwall) are based on the model geometry and
defined automatically. However, they can be also set manually. The values of other
parameters are set automatically by Flow Simulation.
Also you can enable the Advanced channel refinement option to improve the narrow
channel refinement strategy.
APD output:
1 Basic mesh settings:
a) Computational domain size (HCD)
b) Control planes sets for x,y,z direction and cells ratios at these planes
c) Number of cells for x,y,z direction (Nx, Ny, Nz).
2 Refinement settings:
a) Small solid feature refinement level (LSSF)
b) Tolerance refinement level (Ltol) and tolerance criterion (Ctol)
c) Maximum channel refinement level (Lch) and the number of cells per channel
height (Nch=Ngap).

Flow Simulation 2020 Technical Reference 117

Automatic Mesh Setting

Minimum Gap Size and Minimum Wall Thickness

Before creating the initial mesh, Flow Simulation automatically determines the Minimum
gap size and the Minimum wall thickness for the walls contacting the fluid with both
sides. This is required for resolving the essential geometrical features of the model with
the computational mesh. Flow Simulation creates the initial mesh so that the number of
the mesh cells along the normal to the model surface must not be less than a certain
number, if the distance along this normal from this surface to the opposite wall is not less
than the minimum gap size (hgap). This insures that a flow passage or gap with the width
larger than the specified minimum gap size will be resolved with a certain number of cells
across it.
The minimum wall thickness (hwall) governs the resolution of the sharp edges such as tips
of thin fins and it does not influence the meshing if it is equal to or greater than the
minimum gap size (hgap). The default value of hwall is calculated using information about
the model size, the Computational Domain, volume sources, initial conditions, surface
sources and surface goals.
The default value of hgap is the smallest size of the surfaces with the boundary conditions
specified, such as the model inlet and outlet openings in an internal analysis, as well as the
surfaces and volumes with the heat sources, local initial conditions, surface and/or volume
goals and some of the other conditions and features. If set manually, hgap should be set to
the smaller of the size of smallest flow passage to be resolved by the mesh, and the size of
smallest geometric object in the flow that needs to be resolved. At that, it is necessary to
remember that the number of the computational mesh cells generated to resolve the model
geometrical features depends on the specified result resolution level.

Level of Initial Mesh

The Level of initial mesh governs the mean number of cells (Nmean) per the model’s
middle size (Hmean) and the number of cells (Ngap) per the smallest flow passage/channel
height, which specified by the Minimum gap size value (hgap).

First of all, the geometry resolution coefficient (Kres) is defined as the ratio of the model’s
middle size (Hmean) to the minimum gap size (hgap). If it is smaller than 2.5, the number of
cells per the model size (Nmean) and the number of cells per the minimum gap size (Ngap)
are specified by using the table below. Note, that the real number of cells per a gap is
restricted by the maximum allowed channel refinement level (L* ). ch

118
 Table 4: APD Mesh Settings.

Nmean

External
Level of
initial Incompres Compressi
*
mesh Internal sible ble Ngap L ch

1 5 2 3 2 0
2 7 3 5 3 1

3 10 5 7 5 2

4 14 7 10 7 2
5 20 10 14 10 2

6 28 14 20 14 4

7 40 20 28 20 4

See Fig. 5.31 - 5.32 illustrating how the APD technology works
Fig. 5.31Lini = 4; hgap = 18 mm Fig. 5.32Lini = 4; hgap = 4 mm

If Kres is greater than 2.5, the Nmean value is multiplied by the normalized coefficient
Knorm depended on the Kres. Also in these cases the maximum allowed channel refinement
level (L* ) cannot provide a sufficient number of cells per a gap. So if the K >>1, it is
ch res
recommended to enable the Advanced channel refinement option.

Flow Simulation 2020 Technical Reference 119

Local Mesh Settings

Advanced channel refinement

When checked, an improved narrow channel refinement strategy is effective which
ensures that narrow channels/passages are resolved by a sufficient number of cells to
predict the flow and heat transfer phenomena including boundary layers with high
accuracy. The default Maximum channel refinement level (Lch) is set to the value which
provides the number of cells across a gap, no less than Ngap. Otherwise, it depends on the
Level of initial mesh (Lini) and cannot be greater than four, resulting in the restriction of
the minimum size of cells across the model’s flow passage in the
normal-to-solid/fluid-interface directions. The consequence may be a significant increase
of the number of cells up to one or more orders of magnitude.

Local Mesh Settings

The local mesh settings option is one more tool to help create an optimal mesh. Use of
local mesh settings is especially beneficial if you are interested in resolving a particular
region within a complex model.
The local mesh settings influence only the initial mesh and do not affect the basic mesh in
the local area, but are basic mesh sensitive in that all refinement levels are set with respect
to the basic mesh cell.
The local mesh settings are applied to the cells intersected with the local mesh region
which can be represented by a component, face, edge, vertex, or by simple geometric
shape such as a cube, cylinder or sphere. By using Equidistant Refinement option, you
can specify settings to refine the cells within the region bounded by surfaces that are
equidistant (keep the same distance) from the selected objects. You can also apply local
mesh settings to fluid regions and solid bodies. To apply the local mesh settings to a fluid
region you need to specify this region as a solid part or subassembly and then disable this
component in the Component Control dialog.
If a cell intersects with different local mesh setting regions, the refinement settings in this
cell will be used to achieve the maximum refinement.
The local mesh settings have higher priority over the initial mesh settings. Therefore, the
local mesh cells will be split to the specified local refinement levels regardless of the
general Lssf, Lcur and Lch (specified in the Global Mesh dialog). This, however, may
cause refinement of cells located outside of the local region due to imposing the Cell
Mating rule.

120
 Fig. 5.33The local mesh settings used: Two narrow channels are refined to have 10 cells across them.

Recommendations for Creating the Computational

Mesh
1 First, define an appropriate problem type: Internal\External,
Compressible\Incompressible, 3D\2D.
2 In case of Internal flow analysis, first of all, using information about the model
geometry, estimate the ratios of overall model size to sizes of its smallest parts (or to
heights of the narrowest flow passages).
a) Note, that if you use some of engineering methods and techniques implemented in
Flow Simulation (such as “thin-channel” and “thin-boundary-layer” approaches,
see "Two-Scales Wall Functions Model" on page 86), in the most cases these
approaches provide the good accuracy, even on a coarse mesh.
b) If the estimated ratios are not very large, it is recommended to create the mesh
using the default (Automatic) mesh settings. Start with the Level of initial mesh of
3 if both the model geometry and the flow field are relatively smooth. For more
complex problems we recommend to perform the calculation at the result
resolution level of 4 or 5. On this stage it is important to recognize the appropriate
value of the Minimum gap size which will provide you with the suitable mesh.
The default value of the Minimum gap size is calculated using information about
the overall model dimensions, the Computational Domain size, and area of
surfaces where conditions (boundary conditions, sources, etc.) and goals are
specified.
c) If the estimated ratio of overall model size to heights of the flow passages is rather
large, closely analyze the obtained automatic mesh, paying attention to the total
numbers of cells, resolution of the regions of interest and narrow channels. If the
automatic mesh does not satisfy you and changing of the Minimum gap size value

Flow Simulation 2020 Technical Reference 121

Recommendations for Creating the Computational Mesh

do not give the desired effect you can proceed with the Manual mesh with the
enabled Channels refinement.
d) If the Channels refinement does not allow to resolve the area of interest, try to use
the Local mesh settings. In the most cases it makes sense to specify the Refining
Cells by type. Note, that it is strictly recommended to construct a mesh as uniform
as possible, especially in the high-gradient flow regions. So define the local mesh
region so that its boundaries are not intersect the high-gradient flow regions.
e) To capture the relatively small solid features or to resolve the boundary between
substances (fluid/solid, fluid/porous, porous/solid interfaces or boundary between
different solids) you can use either the Local mesh settings or the Advanced
Refinement. The Local mesh settings are preferable if there are few small features
required to be resolved.
3 In case of External flow analysis, it makes sense to create the mesh using the default
(Automatic) mesh settings.
a) Start with the Level of initial mesh of 3 if both the model geometry and the flow
field are relatively smooth. For more complex problems we recommend, first of
all, to perform the calculation at the result resolution level of 4 or 5.
b) As described above, note, that if you use some of Flow Simulation engineering
techniques, in the most cases these approaches provide the good accuracy, even on
a coarse mesh.
c) Closely analyze the obtained automatic mesh, paying attention to the total numbers
of cells and resolution of the regions of interest. If the automatic mesh does not
satisfy you, you can proceed with the Manual mesh with custom Control Planes
and/or specify the Local mesh settings.
d) Use Control Planes to optimize the Basic mesh by setting custom control planes
surrounding the region of interest and assigning the proper Ratio values to the
respective intervals. Note, that it is strictly recommended to construct a mesh as
uniform as possible, especially in the high-gradient flow regions.
e) If needed, try to use the Local mesh settings to resolve the area of interest. It
makes sense to specify the Refining Cells by type or the Equidistant Refinement
depending on the selected region.
4 In the high-gradient flow regions, it makes sense to use the Solution-Adaptive
Refinement. At the beginning perform calculations with the disabled
Solution-Adaptive Refinement to obtain the flow pattern on a coarse mesh, then
perform calculations with enabled Solution-Adaptive Refinement to achieve the
prescribed solution accuracy.
The problem of resolving a model with the computational mesh is always model-specific.
In general, a denser mesh will provide better accuracy but you should tend to create an
optimal mesh and to avoid redundant refinement.
When performing a calculation, try different mesh settings and analyze the obtained
results carefully in order to understand whether it is necessary to refine the mesh or a
coarser resolution is acceptable for the desired accuracy.

122
Solution-vs.-Mesh Convergence Investigation
As it is mentioned in "Numerical Solution Technique" on page 83, Flow Simulation
solves the governing equations numerically, in a discrete form on a computational mesh.
The accuracy of solution of the mathematical problem depends on how well the
computational mesh resolves the regions of a non-linear behavior in the problem. To
provide a good accuracy, the mesh has to be rather fine in these regions. Moreover, the
usual way of estimating the accuracy of the solution consists in obtaining solutions on
several different meshes, from coarse to fine. So, if beginning from some mesh in this set,
the difference in the physical parameters of interest between the solutions obtained on the
finer and coarser meshes becomes negligible from the standpoint of the solved problem
(the solution flattens), then the accuracy of the solution of the mathematical problem is
considered to be attained, since the so-called solution mesh convergence is attained.
The formal method of establishing the mesh convergence requires a curve of a critical
result parameter (typically local or integral value) in a specific location, to be plotted
against some measure of mesh density. Naturally, the Flow Simulation mesh convergence
curve may looks non-monotonic. The reason for such behavior is that a number of
engineering methods and techniques are implemented in Flow Simulation and different
engineering techniques or their combinations are used automatically as the mesh gets finer
or coarser (for example, see "Two-Scales Wall Functions Model" on page 86).
Fig. 5.34The critical result parameter vs.the total number of computational mesh cells.

Flow Simulation 2020 Technical Reference 123

Solution-vs.-Mesh Convergence Investigation

Therefore, the strategy of establishing mesh convergence with Flow Simulation consists,
first of all, in performing several calculations on the same basic project (with the same
model, inside the same computational domain, and with similar boundary and initial
conditions) varying only the computational mesh and controlling if the same engineering
techniques are enabled on the given mesh. Note moreover, that using engineering methods
and techniques in Flow Simulation allows to obtain acceptable accuracy on coarser
meshes as compared with traditional CFD codes.

124
 6
Calculation Control Options

The Calculation Control Options dialog box introduced into Flow Simulation allows you
to control:
• conditions of finishing the calculation,
• saving of the results during the calculation,
• refinement of the computational mesh during the calculation,
• freezing the flow calculation,
• time step for a time-dependent analysis,
• number of rays traced from the surface if radiating heat transfer is enabled.
This dialog box is accessible both before the calculation and during the calculation. In the
last case the new-made settings are applied to the current calculation starting from the
next iteration.
The main information on employing the options of Finishing the calculation and
Refining the computational mesh during calculation is presented in this document.

Finishing the Calculation

Flow Simulation solves the time-dependent set of equations for all problems, including
steady-state cases. For such cases it is necessary to recognize the moment when a
steady-state solution is attained and therefore the calculation should be finished. A set of
independent finishing conditions offered by Flow Simulation allow the user to select the
most appropriate conditions and criteria on when to stop the calculation. The following
finishing conditions are offered by Flow Simulation:
• maximum number of refinements;
• maximum number of iterations;

Flow Simulation 2020 Technical Reference 125

Finishing the Calculation

• maximum physical time (for time-dependent problems only);

• maximum CPU time;
• maximum number of travels;
• convergence of the Goals.
Note: Travel is the number of iterations required for the propagation of a disturbance
over the whole computational domain. Current number of iterations per one travel is
presented in the Info box of the Calculation monitor.
In Flow Simulation you can select the finishing conditions that are most appropriate from
your viewpoint to solve the problem under consideration, and specify their values. For the
latter two conditions (i.e., for the maximum number of travels and the Goals convergence
settings) Flow Simulation presents their default values (details are described below). You
can also select the superposition mode for multiple finishing conditions in the Finish
Conditions value cell: either to finish the calculation when all the selected finishing
conditions are satisfied or when at least one of them is satisfied.
In any case, information on the finishing conditions due to which the calculation has
finished is shown in the Monitor’s Log box.
The Goals convergence finishing condition is complex since it consists of satisfying all the
specified Goals criteria. A specified Goal criterion includes a specified dispersion, which
is the difference between the maximum and minimum values of the Goal, and a specified
analysis interval over which this difference (i.e., the dispersion) is determined. The
interval is taken from the last iteration rearwards and is the same for all specified Goals.
The analysis interval is applied after an automatically specified initial calculation period
(in travels), and, if refinement of the computational mesh during calculation is enabled,
after an automatically or manually specified relaxation period (in travels or in iterations)
since the last mesh refinement is reached. As soon as the Goal dispersion obtained in the
calculation becomes lower than the specified dispersion, the Goal is considered
converged. As soon as all Goals included in the Goals convergence finishing condition (by
selecting them in the On/Off column) have converged, this condition is considered
satisfied. The Goals not included into the Goals convergence finishing condition are used
for information only, i.e., with no influence on the calculation finishing conditions.
Let us consider the Flow Simulation default values for the maximum number of travels
and the Goals convergence settings in detail. These default (recommended by Flow
Simulation) values depend on the Result resolution level either specified in the Wizard or
changed by pressing the Reset button in the Calculation Control Options dialog box. For
higher Result resolution levels the finishing conditions are tighter.
The default maximum number of travels depends on
• the type of the specified Goal (i.e., dynamic or diffusive, see below);
• the specified Result resolution level;

126
 • the problem's type (i.e., incompressible liquid or compressible gas, low or high
Mach number gas flow, time-dependent or steady-state).
Note: The Dynamic goals are: Static Pressure, Dynamic Pressure, Total Pressure,
Mass Flow Rate, Forces, Volume Flow Rate, and Velocity.
The Diffusive goals are: Temperature, Density, Mass in Volume, Heat flux, Heat
transfer rate, Concentrations, Mass Flow Rate of species, and Volume Flow Rate of
species.
The default Goals convergence settings are the default analysis interval, which is shown
in the Finish tab of the Calculation Control Options dialog box, and the default Goals
criterion dispersion values, which are not shown in the Calculation Control Options
dialog box, but, instead, are shown in the Monitor’s Goal Table or Goal Plot table (in the
Criteria column), since they depend on the values of the Goal physical parameter
calculated in the computational domain, and therefore are not known before the
calculation and, moreover, can change during it. In contrast, the Goals criterion dispersion
values specified manually do not change during the calculation.
As for the automatically specified initial calculation period (measured in travels), it
depends on the problem type, the Goal type, and the specified Result resolution level.
Note:
• the manually specified analysis interval for the Goals convergence finishing criteria
must be substantially longer than the typical period of the flow field oscillation (if it
occurs);
• the Goals determined on solid/fluid interfaces or model openings, as well as the
Post-processor Surface Parameters, yield the most accurate and correct numerical
information on flow or solid parameters, especially integral ones;
• Global Goals yield the most reliable information on flow or solid parameters, although
they may be too general;
• the CPU time depends slightly on the number of the specified Goals, but, in some cases,
vary substantially in the case of presence of a Surface Goal;
• Surface and Volume Goals provide exactly the same information that may be obtained
via the Surface and Volume Parameters Post-processor features, respectively.

Refinement of the Computational Mesh During

Calculation
Solution-adaptive refinement of the computational mesh during the calculation is a
process of splitting the computational mesh cells in areas where the calculation error
(specifically, the local truncation error (LTE)) may be sufficiently large and merging of
the computational mesh cells in areas where the calculation error is definitely small.

Flow Simulation 2020 Technical Reference 127

Refinement of the Computational Mesh During Calculation

In the solution-adaptive refinement the computational mesh cells are split until the
specified Refinement level is satisfied. The Refinement level specifies how much times
the initial mesh cells can be split to achieve the solution-adaptive refinement criteria, and
thus governs the minimum computational mesh cell size. If specified Refinement level
allows, the re-meshing adaptation cycles are performed with a sequentially increasing
refinement level. If the Refinement level is already reached, the re-meshing occurs
anyway, but at a constant refinement level.
If any local initial mesh is specified in the project, you can control the solution-adaptive
refinement in these regions independently.
Fig. 6.1 Refinement in the regions where the local initial meshes are specified.
Global Domain level = 1
Local Mesh 1 level = 2 Local Mesh 2 level = Use Global Local Mesh 3 level = Disabled

To locate regions of the computational domain that need mesh refinement, it is necessary
to analyze the solution obtained with the mesh that existed at the last refinement cycle.
Indicator functions are used to locate regions where the solution requires a mesh
refinement to reduce the local truncation errors. The output of an indicator function is used
to determine if the mesh cell must be split or merged or is adequately resolved.
The indicator function for the momentum conservation law is defined as follows:

 n   n 
C m    , foam   h  S ,   1    1   foam (6.1)
 n min   n min 

128
 where S  0.5 Sij  Sij is a convolution of the strain rate tensor, n is the total number of
adjacent cells, nmin is the number of the coordinate directions (plus or minus) in which at
least one cell is adjacent cell, foam is the cell level gap displacement factor, which shifts
the region where the mesh cells must be split to the area where the local truncation error is
definitely small, h is the characteristic cell size defined as follows:

u i

i j j x j
hj
h  (6.2)
u i

i j j x j

where u is the fluid velocity.

The indicator functions for the mass, energy and species conservation laws are defined as
follows:
1 1 1
  1  2 
2
  1 T  2 
2
  y  2  2

C       hi   , Ce     
 hi   , C y     hi   (6.3)
 i   xi    i  T xi    i  xi  

The mesh cell are split if the following condition is satisfied:


Cm   split
m
or C    split or C e   split
e
or C y   split
y

and the mesh refinement level is less than the current maximum refinement level.
The child mesh cells are merged if all following condition is satisfied for each child cell:

Cm   merge
m
and C    merge

and C e   merge
e
and C y   merge
y

and if the refinement level difference between the resulting merged mesh cells and their
neighbors will not exceed 1.
All limiting values (split and merge) are determined automatically at each refinement
cycle by using mesh histograms.
The solution-adaptive refinement can dramatically increase the number of cells so that the
available computer resources (physical RAM) will not be enough for running the
calculation. To limit the total number of cells the Approximate Maximum Cells value can
be specified. If the maximum cells number is exceeded, the number of cells will be limited
for the last refinement cycle in such a way that the LTE, increasing at the refinement level
gap surface, is minimized.
For a transient analysis the following three strategies are available:
• Periodic refinement;

Flow Simulation 2020 Technical Reference 129

Refinement of the Computational Mesh During Calculation

• Tabular Refinement;
• Manual Only refinement.
In the first two strategies the refinement moment is known beforehand. The solution
gradients are analyzed over iterations belonging to the Relaxation interval, which is
calculated from the current moment rearwards. As the result, only steady-state gradients
are considered. The default length of the Relaxation interval can be adjusted manually.
On the other hand, the analysis must not continue with the same relaxation interval
defined from the start of the calculation, in order to avoid considering the initial highly
unsteady period. Therefore, a period of at least two relaxation intervals is recommended
before the first refinement. If the first assigned refinement is scheduled in a shorter term
from the beginning, the period over which the gradients are analyzed is shortened
accordingly, so that in an extreme case it can be as short as one current iteration. If you
initiate a refinement manually within this period, the gradients are analyzed in one current
iteration only. Naturally, such a short period gives not very reliable gradients and hence
may result in an inadequate solution or excessive CPU time and memory requirements.
The Fig. 6.2 illustrates this concept. Here, the letter r denotes the relaxation interval. This
figure involves both the Periodic and Tabular refinements. Case 1 is the recommended
normal approach. In the Case 2 the first refinement is too close to the starting point of the
calculation, so the gradients are analyzed over the shorter interval (which could even be
reduced to just one current iteration in an extreme case). Case 3 is a particular case when
the refinement is initiated manually just before a previously assigned refinement. As the
result, the manual refinement is well-defined, since the gradients have been analyzed over
almost the entire relaxation interval, but on the other hand, the previously assigned
refinement is performed on the substantially shorter interval, and therefore its action can
be incorrect. Thus, Case 3 demonstrates the possible error of performing manual and
previously assigned refinements concurrently.

130
 Fig. 6.2 Refinement strategy.

Collecting of the statistics is prohibited

Monitoring interval where statistics are collected
Refinement

Case 1

r r

0 Ref. point 1 Ref. point 2

Case 2 Case 3

r r1 r2
r

0 Ref. point 1 Manual ref. Auto ref.

The mesh refinement performed during the calculation is idling and the warning messages
appear in the Monitor window in the following cases:
• Refinement canceled: the limit of maximum number of cells was achieved. This
message warns that the Approximate Maximum Cells is reached.
• Refinement canceled: refinement-unrefinement criteria is not satisfied. This
message indicates that there are no mesh cells where the conditions of splitting of
merging are satisfied.
• Refinement canceled: the limit of maximum refinement level was achieved in
field gradients. If there are mesh cells where the conditions of splitting or merging
are satisfied, but the Refinement level has been already reached, since the field
gradients have not been changed, this message appears. If the field gradients are
changed and the mesh refinement is performed again, the mesh refinement will
result in splitting or merging cells.

Flow Simulation 2020 Technical Reference 131

Flow Freezing

Flow Freezing

What is Flow Freezing?

Sometimes it is necessary to solve a problem that deals with different processes
developing at substantially different rates. If the difference in rates is substantial (10 times
or higher) then the CPU time required to solve the problem is governed almost exclusively
by the slower process. To reduce the CPU time, a reasonable approach is to stop the
calculation of the fastest process (which is fully developed by that time and does not
change further) and use its results to continue the calculation of the slower processes. Such
an approach is called “freezing”.
In the case of problems solved with Flow Simulation the processes of convective mass,
momentum, and energy transport are the fastest processes to develop and to converge,
whereas the processes of mass, momentum, and energy transfer by diffusion are the
slowest ones. Accordingly, Flow Simulation offers the “Flow Freezing” option that allow
you to freeze, or fix, the pressure and velocity field while continuing the calculation of
temperature and composition. This option is especially useful in solving steady-state
problems involving diffusion processes that are important from the user’s viewpoint, e.g.
species or heat propagation in dead zones of the flow. Time-dependent analyses with
nearly steady-state velocity fields and diffusion processes developing with time are also
examples of this class of problems. As a result, the CPU time for solving such problems
can be substantially reduced by applying the Flow Freezing option.
Flow Simulation treats Flow Freezing for the High Mach number flows differently. All
flow parameters are frozen, but the temperature of the solid is calculated using these fixed
parameters at the outer of the boundary layer and user defined time step. Temperature
change on the solid's surface and relevant variation of the heat flows are accounted in the
boundary layer only. It is impossible (and makes no sense) to use Flow Freezing for
calculation of concentration propagation in the High Mach number flow. If custom time
step is not specified, the steady-state temperature of solid will be reached in one time step
assumed to be infinite.

How It Works
To access the Flow Freezing option, open the Calculation Control Options dialog box,
then the Solving tab. This option has three modes: Disabled (by default), Periodic, and
Permanent.

Flow Freezing in a Permanent Mode

As an example of applying the Flow Freezing option, let us consider a plane flow (2D)
problem of heating the vortex core in a vessel (Fig. 6.3).

132
 Fig. 6.3 Heating the vortex core in a vessel.

At the beginning the entire fluid region is filled with a cold (T=300 K) liquid. A hot
(T=400K) liquid enters the vessel through the lower channel (the upper channel is the
exit). As a result, a vortex with a cold core is developed in the vessel. The vortex core
temperature is changed mainly due to heat diffusion. To measure it, a small body is placed
at the vortex center and disabled in the Component Control dialog box, so that it is
treated by Flow Simulation as a fluid region. Its minimum temperature (i.e., the minimum
fluid temperature in this region) is the Volume Goal of the calculation.
First of all, let us consider Flow Freezing operating in the Permanent mode. The only
user-specified parameter in Permanent mode is the starting moment of enabling the Flow
Freezing option. Until this moment the calculation runs in a usual manner. After this
moment the fluid velocity field becomes frozen, i.e., it is no longer calculated, but is taken
from the last iteration performed just before the Flow Freezing Start moment. For the
remainder of the run only the equations’ terms concerning heat conduction and diffusion
are calculated. As a result, the CPU time required per iteration is reduced.
The starting moment of the Flow Freezing option should be set not too early in order to let
the flow field to fully develop. As a rule, an initial period of not less than 0.25 travels is
required to satisfy this condition. In most problems the 0.5 travel initial period is
sufficient, but there are problems that require a longer initial period.
Tip: The Flow Freezing Start moment, as well as other parameters of the Calculation
Control Options dialog box can be changed during a calculation.
Tip: As soon as the Flow Freezing option is invoked, only the slowest processes are
calculated. As a result, the convergence and finishing criteria can become non-optimal.
Therefore, to avoid obtaining incorrect results when enabling the Flow Freezing option, it
is recommended to increase the maximum number of travels specified at the Finish tab of
the Calculation Control Options dialog box by 1.5…5 times compared to the number that
was set automatically or required for the calculation performed without the Flow
Freezing option.

Flow Simulation 2020 Technical Reference 133

Flow Freezing

When first solving the problem under consideration we set the maximum number of
travels to 10. The calculation performed without applying the Flow Freezing option then
required about 10 travels to reach the convergence of the project Goal (the steady-state
minimum fluid temperature in the vortex core). However, the steady-state fluid velocity
field was reached in about 0.5 travels, i.e., substantially earlier. So, by applying the Flow
Freezing option in the Permanent mode (just after 0.5 travels) the same calculation
requires substantially less time on the same computer to reach the convergence of the
project Goal.
If it is necessary to perform several calculations with the same fluid velocity field, but
different temperatures and/or species concentrations, it is expedient to first calculate this
fluid velocity field without applying the Flow Freezing option. Then, clone the Flow
Simulation project into several projects (including copying the calculation results), make
the required changes to these projects, and perform the remaining calculations for these
projects using the calculated results as initial conditions and applying the Flow Freezing
option in the Permanent mode with a zero Start period.
Tip: If you forget to use the calculated results as initial conditions, then the saved fluid
velocity field will be lost in the cloned project, so the project must be created again. To use
the calculated results as initial conditions for the current project, select the Transferred
type of Parameter definition for the initial conditions in the General Settings dialog box.

Flow Freezing in a Periodic Mode

In some problems the flow field depends on temperature (or species concentrations), so
both the velocity and the temperature (concentrations) change simultaneously throughout
the calculation. Nevertheless, since they change in a different manner, i.e., the velocity
field changes faster than the temperature (concentrations) field, therefore approaching its
steady state solution earlier, the Flow Freezing option can be used in a Periodic mode to
reduce the CPU time required for solving such problems. The Periodic mode of the Flow
Freezing option consists of calculating the velocity field not in each of the iterations (time
steps), but periodically for a number of iterations specified in No freezing (iterations)
after a period of freezing specified in the Freezing (iterations) (see Fig. 6.4) The
temperatures and concentrations are calculated in each iteration. Examples include
channel flows with specified mass flow rates and pressures, so the fluid density and,
therefore, velocity depend on the fluid temperature, or flows involving free convection,
where due to the buoyancy the hot fluid rises, so the velocity field depends on the fluid
temperature.
Fig. 6.4 The Periodic mode of the Flow Freezing.

No Freezing
Start

Iterations

Freezing

134
 As an example, let us consider a 3D external problem of an air jet outflow from a body
face into still air (see Fig. 6.5, in which the jet outflow face is marked by a red line). Here,
the wire frame is the computational domain. The other body seen in this figure is
introduced and disabled in the Component Control dialog box (so it is a fluid region) in
order to see the air temperature averaged over its face (the project Goal), depending on the
air temperature specified at the jet outflow face.
Fig. 6.5 Air jet outflow from a body face into a still air.

This problem is solved in several stages. At the first stage, the calculation is performed for
the cold (T = 300 K, which is equal to the environment temperature) air jet. Then we clone
the project including copying the results. Next, we set the outlet air temperature to T = 400
K, specify the Periodic mode of the Flow Freezing option by its Start moment of 0.25
travels (in order for the heat to have time to propagate along the jet to the measuring face)
and under Duration specify 10 as both the Freezing (iterations) and No freezing
(iterations) values. Then perform the calculation on the same computational mesh with
the Take previous results option in the Run box. As a result, the calculation with flow
freezing takes less CPU time than the similar calculation without the Flow Freezing option
enabled.

Radiation Freezing
Radiation calculation might be CPU consuming and radiation process might develop
faster than fluid flow. To reduce the CPU time, the radiation calculation can be stopped
and you can use the previous iteration values of all radiation parameters to continue the
calculation.
To access the Radiation Freezing option, open the Calculation Control Options dialog
box, then the Solving tab. This option has three modes: Disabled (by default), Periodic,
and Auto.

Flow Simulation 2020 Technical Reference 135

Radiation Freezing

Radiation Freezing in a Periodic Mode

The Periodic freezing strategy allows to calculate the flow field and the radiation
simultaneously throughout the calculation. Also this strategy allows to specify (in the
iterations) the Periodicity, which means that the radiation calculates every specified
iterations only, which allows to reduce the CPU time required for solving such problems.

Radiation Freezing in an Automatic Mode

The Automatic freezing strategy means that the radiation calculates every freezing step
only and this strategy is similar to the Periodic mode, however the Automatic mode allows
to vary the freezing period automatically and it depends on the selected radiation model:
• If the Discrete Transfer model is selected, the freezing step can be automatically
varied from 3 to 10 iterations.
• If the Discrete Ordinates model is selected, the freezing strategy is disabled.
Also in case of the Time-dependent (transient) or Heat conduction in solids only
analysis, the freezing strategy is disabled too.

136
 References

1 Reid R.C., Prausnitz J.M., Poling B.E. (1987). The properties of gases and liquids, 4th
edition, McGraw-Hill Inc., NY, USA.
2 Idelchik, I.E. (1986). Handbook of Hydraulic Resistance, 2nd edition, Hemisphere,
New York, USA.
3 Henderson, C.B. Drag Coefficients of Spheres in Continuum and Rarefied Flows.
AIAA Journal, v.14, No.6, 1976.
4 Carlson, D.J. and Hoglund, R.F. Particle Drag and Heat Transfer in Rocket Nozzles.
AIAA Journal, v.2, No.11, 1964.
5 ASHRAE Handbook–2001 Fundamentals.
6 ISO 7726:1998, Ergonomics of the Thermal Environment – Instruments for Measuring
Physical Quantities.
7 ISO 7730:2005, Ergonomics of the thermal environment – Analytical determination
and interpretation of thermal comfort using calculation of the PMV and PPD indices
and local thermal comfort criteria.
8 Roache, P.J., (1998) Technical Reference of Computational Fluid Dynamics, Hermosa
Publishers, Albuquerque, New Mexico, USA.
9 Hirsch, C., (1988). Numerical computation of internal and external flows. John Wiley
and Sons, Chichester.
10 Ginzburg, I. P., (1970). Theory of Drag and Heat Transfer. Leningrad, LGU (in
Russian).
11 Van Driest, E.R. On Turbulent Flow Near a Wall. Journal of the Aeronautical Science,
v.23, No.10, p.1007. 1956.
12 Glowinski, R., Le Tallec, P. (1989). Augmented Lagrangian Methods and
Operator-Splitting Methods in Nonlinear Mechanics. SIAM, Philadelphia.
13 Marchuk, G.I., (1982). Methods of Numerical Mathematics, Springer-Verlag, Berlin.
14 Samarskii, A.A., (1989). Theory of Difference Schemes, Nauka, Moscow (in Russian).

Flow Simulation 2020 Technical Reference 137

 15 Patankar, S.V., (1980). Numerical Heat Transfer and Fluid Flow, Hemisphere,
Washington, D.C.
16 Saad, Y. (1996). Iterative methods for sparse linear systems, PWS Publishing
Company, Boston.
17 Hackbusch, W. (1985). Multi-grid Methods and Applications, Springer-Verlag, NY,
USA.
18 Faeth. G. M. (1983) Evaporation and Combustion of Sprays. Progress in Energy and
Combustion Science, Pergamon Press, 9:1–76.

138