as density, t as time and xi , xj as space coordinates. The R(T) = ¼ T (1 – T) stands for Air Movement Within Enclosed Road-Objects with Contra-Traffica CFD-Investigation
75
reaction rate[24] where the reciprocal value of reaction time-scale is represented by 1/T and
is here again a constant average density. Temperature T will be used as expression for reaction-progress-variable as well, whose purpose is to distinguished burned, unburned and
partially burned state, providing facile interpretation of flame propagation. The term g
i
T
denotes buoyancy, treated according to the Boussinesq approximation, where T is showing the difference between local and reference temperature. The symbol g denotes the gravity and
is coefficient of the thermal linear deformation. The model for the unknown Reynolds
stress[25], ij, s is related to the local strain rate:
( ) ( )
ij
ij N
ij T (4)
where we distinguish between Newtonean stress ( ) 2
ij N
i
S j featuring molecular
viscosity; and turbulent Reynolds stress ( ) 2
ij T
T i
S j featuring turbulent viscosity
recognising again i
S j as local strain stress rate reads:
1
v
v
Sij
j
i
(5)
2 x
x
j
i
2
k
Here we recognise the turbulent viscosity:
T
C
( 6 )
The applied k- model[26] is a two-equation eddy viscosity model [27, 28] uses transport
equations for these two variables[29]. One of these equations governs the distribution through the field of k, the local kinetic energy of the fluctuating motion. The other one is explaining a turbulence characteristic of different dimensions, the energy dissipation rate [30].
2
2
k
k
k
2
C
k C
P
N (7)
t
d
v
2
k
C
C kP
d
2
C N C
1
3 v
2 (8)
t
k
The d
P , the production term for the squared shear frequency is:
T
d
P
V
V 2
1
(9)
2
where T
V
V is the 2-norm matrix. The constants are given: C1 = 0.126, C2 = 1.92, Cµ =
0.09, C= 0.07. The variable N denotes the Brunt-Väisälä frequency, with the N2, related
following squared buoyancy:
v g
2
N v
(10)
t
0
z
Here is constant 1
t
and following expressions of eddy coefficients vv , v should include
the stability parameters to account for the turbulence damping in the stratified fluid flows: 76
Applied Computational Fluid Dynamics
2
k
2
k
v S
S
v
u
v ,
(11)
v
b
b
0.108 0.0229
S
N
,
0.177
S
(12)
u
1 0.47 0.0275
b
N
N
1
0.043 N
where b =10-6 and the following stability coefficients N and C3 can be expressed as: 2
k
2
N
N
2
2
C
N (13)
3
0.4;
0
2
C N
3
1;
0
During these investigations, the wall shear stress was obtained from the logarithmic law of
wall (“wall function”) for the distance from the wall, y, the von Karman constant =0,41 and v
1
the constant C=9,793:
v y
ln
C
(14)
v
v
Energy transfer within the computational domain of such non-premixed turbulent flow in
tunnel cavity is governed by next equation:
H vH
t H Q (15)
t
c
h
p
where contribution to the energy budget, characterised through turbulent thermal
conductivity t is performed through the heat conduction and diffusion combined in term
t H , accompanied by the heat-increase component
Q due to the chemical reaction
c
h
p
rate R of a species n and it’s molar mass M:
0
h
Q n
h
n
R (16)
n Mn
In the governing energy equation is total enthalpy H based onto particular enthalpy of the
nth species Hn and it’s mass fraction Yn
H n
Y Hn (17)
n
with
0
H h
n
n
hn (18)
having so species n characterised through it’s thermal capacity cp,n and nth enthalpy
T
h
T
n
cp, ndT with referent temperature
293,15
ref
K that makes an impact on the
ref
T
formation enthalpy 0 n
h r
T
.
ef , n
Air Movement Within Enclosed Road-Objects with Contra-Traffica CFD-Investigation
77
The radiative transfer equation (RTE) for an absorbing, emitting and scattering medium at
the position r in the direction s reads:
4
4
dI( r , s)
2 T
( a ) I( r , s) ax
s
I( r , s '
) ( s, s ')
d '
s
(19)
ds
4 0
The applied Discrete Ordinate approach for treating the radiation considers the RTE in the
direction s as a field equation, allowing the modelling of non-gray radiation, using the gray-
band model. So, for the spectral intensity I( r , s) the RTE can be written as: 4
2
(
s
I( r , s) s) (
a )
s I ( r , s )
a x Ib
I( r, s')( s, s')
d ' (20)
4 0
Here is r a position vector and s a direction vector. I denotes the radiation intesity, i.e.
radiation intensity I for a certain wawe-lenght . The radiation intensity for the black body reads Ib . Scattering coefficient is represented by s and Stefan-Bolzman-constant reads =5.672x10-8W/m²K4. stands for phase function, T is local temperature and x is
refraction index. Path length is s and s ' denotes direction vector. ' stands for solid angle and a is absorption coefficient, i.e. absorption coefficient for a particular wave-length
a .
So will the RTE be integrated over each wave-length-interval resulting in transport
equation for entity Icontaining the radiant energy of the wave-length . The total
intensity I( r , s) in each direction s at some position r is computed as:
I( r , s) I r s
, ( , )
n
n (21)
n
The weighted-sum-of grey-gases model (WSGGM), that was applied in this study for
handling the irradiance in the complicated phenomenon of the confined combustion[31]
presents a reasonable compromise between the oversimplified grey gas model and a
complete model which takes into account particular absorption bands. Within the Discrete
Ordinate approach, where the soot generation[32] was described by the Magnussen
model[28], the basic assumption of the WSGGM means that the emissivity E of the any nth
gray gas characterised by pSUM, the sum of itś partial pressures p and the absorption
within the computational domain over a distance d is presented by:
i
E a
T
n S
p UMd
E, n 1
e
(22)
n0
Here
n pd
E
a , n denotes the emissivity weighting factor for the nth gray gas and 1
e
is
emissivity of the nth gray gas. The temperature dependence of E
a , n is commonly
J
approximated by the relation:
n1
a
E, n
Eb, nT
(23)
n1
where polynomial coefficient E
b , n , depending on the gas temperature emissivity is obtained
experimentally. The coefficients E
b , n and n are slowly changing for broad range of the
78
Applied Computational Fluid Dynamics
pressure ( 0.001 bar p 10 bar ) and temperature ( 600 K T 2400 K ) which simplifies J
relation (22) to:
E a
E, n n S
p UM
(24)
n1
J
ln(1 E)
for
4
a a p ; s 10 m
4
E, n n
and a
; s
10 m (25)
n1
s
The combustion – the chemistry development is explained by fast chemistry assumption
including the prePDF [33] and in the ideal stoichiometric conditions the reaction runs as
foolows:
ignition
energy
C H
O
N
CO H O
N
7
16
11 2 3,76 2
7
2
8 2
41,36 2
H (26)
However, in a case of unfulfilled combustion, soot is influencing the radiation absorption.
Employed generalised model estimates the effect of the soot onto the radiative heat transfer by determining an effective absorption coefficient for soot. In this case, absorption
coefficient is a sum of the one for the pure gas ag and that one for the pure soot s a :
a a
s g
s
ag (27)
The ag is obtained from presented WSGGM approach and
a b b T
s
1
1
m
T
2000 (28)
with
m²
b
4
1
1232,4
and
4,8 10
1
kg
T
b
, where
K
m denotes soot density.
3. Numerical approach
Discretizing the governing equations in both space and time[28, 34] while choosing the
numerical method of the standard finite volumes [22, 34, 35], the CFD-tool is “adapted” for
the time-dependent modelling of the explored event. The spatial discretization of time-
dependent equations within here-employed segregated solution method are linearized in an
implicit way. Since variables described the flow by their value at the centre of a single mesh-cell in explored computational domain, the system of linear equations occurs, containing
mentioned unknown variable at the cell centre as well as their unknown values in
neighbour cells. A scalar transport-equation within such method is also used to discretize the momentum equations; in the same mode for the pressure-field (if cell-face mass-fluxes
were known). This applies for velocity-field as well. In case that the pressure-field and face mass fluxes are not known, applied commercial software of FLUENT [28] uses a co-located scheme, whereby the local values for pressure and velocity are stored at centres of a single mesh-cell. A need for interfacial values includes an application of an interpolation scheme to compute pressure and velocity out of mesh-cell values. The integration over the observed
control-volume, where is a mesh-cell in a computational domain, can be performed yielding
the discretised equation for the mass-flux through a control-surface – a mesh-face on the
mentioned mesh-volume. Procedure of the segregated solver that follows as a mechanism in
this research, the continuity equation, is used as an equation for pressure as well[28].
Air Movement Within Enclosed Road-Objects with Contra-Traffica CFD-Investigation
79
However, the pressure does not appear for an incompressible flow in explicit way, since the
density is not directly related to pressure. However, through the application of the SIMPLE-
algorithm-family[34] during this investigation (within the commercial software package of
FLUENT) the won values for pressure were introduced into the continuity equation and in this way supported the pressure-velocity coupling. This was done through the NITA-
algorithm for the unsteady.. Further, the SIMPLE uses a relationship between velocity and pressure corrections to enforce mass conservation and to obtain the pressure field. So,
executing these numerical steps, the equations can express the state for each other cell in the computational grid, the mesh. This will result in a set of algebraic equations with a
coefficient-matrix. In this way the segregated solver is “updating” variable-field by
considering all the cells of the domain in the same time, solving the governing equations
sequentially (segregated one from another). What follows are the iterations of the solution-
loop before one wins the converged solution. Subsequently, the next field of another
variable will be solved by again considering entire cells at the same time, etc. Finally,
determining of the boundary conditions relied on some previous done studies [19, 36] that
were basis for the estimation.
4. Objects of investigation - Two different tunnel-shapes
Todayś major two shapes of the underground traffic (road) facilities that are having their
cross-section as a “horse-shoe” and rectangular shape – we explored here. And according to
the final use, these traffic covered object, were taken as very frequently used ones.
Additionally, those tunnels (one, having bifurcation – is going under Slovenian Capital of
Ljubljana, and another is situated in a curvy valley of the river Bosna, between Bosnian cities of Travnik and Zenica) are “offering” interesting geometry for a CFD-based investigation;
such is air-movement and the movements of air-pollutants and other gaseous products (in
accidental situations).
4.1 The “Sentvid” tunnel
As mentioned, the northern tube of tunnel “Sentvid” (undergoing the city of Ljubljana) was
investigated with the CFD-approach, where the “chimney-effect” was expected (coming