0 Votes

Models - Calculation of airflow

Last modified by Walter Fornari on 2026/08/26 10:28

Back to Start | Next (Calculation of heat flow)

TABLE OF CONTENT

Calculation of airflow

The models

Calculation of airflow

The airflow in any section of the tunnel system is derived from balances of air mass flow and total pressure for the whole tunnel. The unsteady flow energy equation is the basis for the calculations and the energy losses (or gains) are modeled similarly to normal pipe flow calculations.

Loss coefficients are in some cases calculated by the program, but can also be given by the program user. The present lack of detailed knowledge for some of the airflow configurations has been overcome by continuity assumptions and comparisons to measurements of loss coefficients.

Rail tunnel section (RailSection)

In the tunnel section model the length is divided at nodes into segments. Each segment has constant cross section area and hydraulic diameter. The expressions of the different pressure contributions are in all essentials the same as those used in the SES program [SES User’s Manual 1997] and described in [Harris 1973] and [DOT 1997].

Note that in the computation of pressure changes due to friction and piston action we neglect the effect of vehicles side by side.

The different contributions from segments are

$∆p_{frict}=-\sum_{i}{\frac{\rho_i}{2}\frac{\lambda_{Ti}}{d_i}}\left[\left(l_i-\sum_{j} l_{V_{ij}}\right)v_i\left|v_i\right|-\sum_{j}{\frac{l_{V_{ij}}}{\left(1-\frac{A_{V_j}}{A_i}\right)^3}\left(\frac{A_{V_j}}{A_i}v_{V_{ij}}-v_i\right)\left|\frac{A_{V_j}}{A_i}v_{V_{ij}}-v_i\right|}\right]-\sum_{i}\sum_{k}{\frac{\rho_i\lambda_{s_k}}{2}v_i\left|v_i\right|}$

where

Δp_frict	pressure change due to friction against walls (Pa)	dPFrict
λ_Ti	Darcy friction factor for tunnel wall in segment i	lambda
ρ_i	mean density of air in tunnel segment i (kg/m³)	Rho
d_i	hydraulic diameter of tunnel segment i (m)	dTun
l_i	length of tunnel segment i (m)	lTun
l_{V_ij}	length of vehicles of type j in tunnel segment i (m)	lBody
v_i	mean air velocity in segment i, positive from left to right (m/s)	VAir
v_{V_ij}	velocity for vehicle j in segment i, positive from left to right (m/s)	speed
A_{V_j}	cross section area of vehicle type j (m²)	frontArea
A_i	tunnel area of segment i (m²)	aTun
λ_{s_k}	loss coefficient for sign type k	lambda2

$\Delta p_{stack}=\sum_{i}{\Delta p_{stack_i}=\sum_{i}\left(p_{stack_{i+1}}-p_{stack_i}\right)}$

\[ \Delta p_{\text{stack}} = \underbrace{\sum_{i} \Delta p_{\text{stack},i}}_{\text{Sum of individual pressure drops}} = \underbrace{\sum_{i} \left( p_{\text{stack},i+1} - p_{\text{stack},i} \right)}_{\text{Difference between consecutive pressures}} \]

Δp_stack	stack pressure rise (Pa)	dPStack
Δp_{stack_i}	stack pressure rise in segment i(Pa)	dPStack
p_{stack_i}₊₁	$p_{stack_i}\cdot\ e^{-g\cdot\mathrm{\Delta}h_i/(R\cdot(T_i+273.15))}$ = stack pressure in node i+1(Pa)	PStack
p_{stack_}₁	static pressure in left node of section (Pa)	PStack
Δh_i	change of altitude in segment i (m)	dH
g	9.80665 = acceleration of gravity (m/s²)
T_i	air temperature at node i (˚C)	TAir
R	287 = gas constant for air (J/kg/K)

If one is interested in the average velocity around the train in a tunnel segment, it would be:

$V_{real}=\frac{v_i-\sum_{j}{\frac{A_{V_j}}{A_i}v_{V_{ij}}}}{1-\sum_{j}\frac{A_{V_j}}{A_i}}$

This velocity takes the train body into account in the continuity equation and provides the average real velocity around the train. IDA Tunnel calculates this velocity for a maximum of two trains in a segment because the 1D approximation for this velocity around two or more trains is far from accurate. The importance of this physical aspect of velocity lies in its direction, which can be opposite to the direction of train movement.

$\Delta p_{\mathrm{fire\ }}=-FPD\cdot q_{\mathrm{fire\ }}/1000000\cdot\delta_i$

Δp_fire	pressure drop across fire (Pa)
FPD	pressure drop across fire per MW of fire power (Pa/MW)	FPD
q_fire	actual fire power (see 2.2.6) (W)
δ_i	if fire in segment i then if \|m_{dot_i}/A_i\| > 0.1 sign(m_{dot_i}) else sin(π/2·m_{dot_i} /A_i /0.1) else 0
m_{dot_i}	air mass flow in segment i (kg/s)	m_dot

$\Delta p_{\mathrm{fricv\ }}=\sum_{i}\hairsp\sum_{j}\hairsp\frac{\rho_i}{8}\frac{\lambda_{V_-j}\cdot\ l_{V_-ij}\cdot\ p_{V_-j}}{A_i\left(1-A_{V_-j}/A_i\right)^3}\left(v_{V_-ij}-v_i\right)\left|v_{V_-ij}-v_i\right|$

where $\lambda_{V_j}=\lambda_{VS_j}+\frac{c_{DTV_j}}{4l_{V_j}p_{V_j}}$

Δp_fricv	pressure change due to friction against vehicles (Pa)	dPFricv
λ_{V_j}	skin friction coefficient for vehicles of type j	lambda
λ_{VS_j}	skin friction coefficient related to viscous drag for vehicles of type j	lambdas
c_{DTV_j}	drag coefficient weighted total truck area of vehicle type j (m²)	cdTruck
l_{V_j}	length of vehicles of type j (m)	length
p_{V_j}	perimeter of vehicles of type j (m)	Perim

$\Delta p_{\mathrm{piston\ }}=&\sum_{i}\hairsp\sum_{j}\hairsp\frac{\rho_i}{2}A_{V_-j}\left\{\left[\frac{c_{D{FV}_-j}}{A_i}+\frac{\left(2A_i-A_{V_-j}\right)}{\left(A_i-A_{V_-j}\right)^2}\right]N_{\mathrm{front_\ }ij}+\right. \left.\left[c_{DBV_-j}\frac{A_i}{\left(A_i-A_{V_-j}\right)^2}-\frac{2}{A_i-A_{V_-j}}\right]N_{back_ij}\right\}\left(v_{V_-ij}-v_i\right)\left|v_{V_-ij}-v_i\right|$

where $c_{DBV_j}=\frac{0.029}{\sqrt{0.5\lambda_{VS_j}l_{V_j}p_{V_j}/(4A_{V_j})}}$

Δp_piston	piston pressure rise (Pa)	dPPiston
c_{DBV_j}	drag coefficient at back end of vehicle type j	cdBack
c_{DFV_j}	drag coefficient at front end of vehicle type j	cdFront
N_{front_ij}	# of front ends (see below) of vehicle type j in segment i	nFront
N_{back_ij}	# of back ends (see below) of vehicle type j in segment i	nBack

front end is here defined as the end, either physical front or back, which has headwind

Thus the aerodynamic drag force on each vehicle becomes (omitting indices for vehicle and segment)

$F_{drag}=F_{fricv}+F_{piston}$

where $F_{\mathrm{fricv\ }}=\frac{\rho}{8}\frac{\lambda_V\cdot\ l_V\cdot\ p_V}{\left(1-A_V/A\right)^3}\left(v_V-v\right)\left|v_V-v\right|$

$F_{\mathrm{piston\ }}&=\frac{\rho}{2}A\cdot A_V\left\{\frac{c_{DFV}}{A}+\frac{\left(2A-A_V\right)}{\left(A-A_V\right)^2}+c_{DBV}\frac{A}{\left(A-A_V\right)^2}-\frac{2}{A-A_V}\right\}\left(v_V-v\right)\left|v_V-v\right|$

$=\frac{\rho}{2}A_V\left\{c_{DFV}+\frac{c_{DBV}}{\left(1-A_V/A\right)^2}-\frac{A_V}{A\left(1-A_V/A\right)^2}\right\}\left(v_V-v\right)\left|v_V-v\right|$

The contributions from nodes, which cause discontinuities in pressure profile along the section, are

$\Delta p_{\mathrm{area\ }}=-\sum_{i}\hairsp\frac{\rho_i}{2}\cdot\zeta_i\cdot\ v_i^2$ dPArea

where ρ_i= density of air at node i (kg/m³) Rho

v_i = air velocity at node i (with area A_i), positive from left to right (m/s) VAir

ζ_i = resistance factor at node i, depending on areas, shape of transition, and direction of flow. This factor is expressed in terms of expansion/contraction loss coefficients and areas of surrounding segments. Loss coefficients are related to smallest segment area. If areas are equal contraction loss coefficient is valid for positive direction and expansion loss coefficient for negative direction.

$\zeta_i=1-(1-k_{exp_i})\cdot(A_i/A_{i-1})^2$ if A_i / A_i-₁ > 1 and v_i > 0 kPos

$\zeta_i=1-(1+k_{cnt_i})\cdot(A_i/A_{i-1})^2$ if A_i / A_i-₁ > 1 and v_i < 0 kNeg

$\zeta_i=1+k_{cnt_i}-(A_i/A_{i-1})^2$ if A_i / A_i-₁ ≤ 1 and v_i > 0 kPos

$\zeta_i=1-k_{exp_i}-(A_i/A_{i-1})^2$ if A_i / A_i-₁ ≤ 1 and v_i < 0 kNeg

k_{exp_i} = expansion loss coefficient at node i kExp

k_{cnt_i} = contraction loss coefficient at node i kCntr

A_i = cross section area at node i = cross section area of segment to the right (m²)

${\Delta p}_{fan}=\sum_{i}\hairsp\sum_{j}\hairsp\delta_{ij}\cdot\rho_i\cdot\ k_{f_-j}\cdot\left(sign\funcapply\left(c_{act_-j}\right)\cdot\ v_{fan_-j}-v_i\right)\cdot\left|c_{act_-j}\cdot\ v_{f_{fan_-}j}\right|$ thrust control

$=\sum_{i}\hairsp\sum_{j}\hairsp\delta_{ij}\cdot\rho_i\cdot\ k_{f_-j}\cdot\left(c_{act_-j}\cdot\ v_{fan-_-}-v_i\right)\cdot\left|c_{act_-j}\cdot\ v_{f_{fan_-j}}\right|$ speed control

Δp_fan		fan pressure rise (Pa)	dPFan
δ_ij		1 if fans j is in node i, otherwise 0
k_{f_j}		N_{fan_j} · k_{fan_j} · A_{fan_j} / A_i = efficiency factor of fans j	kF
N_{fan_j}		number of jet fans j at node i	mult
A_{fan_j}		cross section area of jet fan j (m²)	area
k_{fan_j}		pressure rise coefficient (installation factor) of jet fan j
v_{f_j}		c_{act_j}· v_{fan_j} = outlet velocity of fan j, positive from left to right (m/s)
c_{act_j}		-τ·dc_{act_j}/dt + c_{tr_j} = delayed fan capacity control	cntrAct
c_{tr_j}		fan capacity control in interval [-1, 1], -1 = right to left action	cntr
τ		time constant for startup (s)	tau
v_{fan_j}		jet speed of fan j (m/s)	v
Q_{f_j}	=	N_{fan_j} · A_{fan_j} ·ρ_i· v_{fan_j} ² / Tp_{fan_j} _i· \| c_{act_j} \| if thrust control N_{fan_j} · A_{fan_j} ·ρ_i· v_{fan_j} ² / Tp_{fan_j} _i· c_{act_j} ²; if speed control power consumption of fans j (W)	PowerCons [kW]
Tp_{fan_j}		total performance (thrust/fan power ratio) of fan j (N/W).	TotPerform [N/kW]

The air density ρ_i at node i is assumed to satisfy the perfect gas law

$\rho_i=\frac{p_i}{R\cdot(T_i+273.15)}$

where

R	287 = gas constant for air (J/kg/K)
p_i	static pressure at node i (Pa)
T_i	air temperature at node i (˚C).	TAir

The air density ρ_i in a segment i satisfies the gas law with a pressure which is the mean of the pressures at the end nodes of the segment.

The relation between mass flow m_{dot_i}, positive from left to right, and the velocity v_i at node i or in a segment i, is given by

$v_i=\frac{m_{dot_i}}{\rho_iA_i}$ .

The tunnel area is initially computed at the nodes between the segments. As tunnel area of a segment the area at the left node of the segment is used.

Road tunnel section (RoadSection)

The road section model is very much alike the rail section model. The main difference occurs because there is another type of traffic: car traffic instead of train traffic. We will just list the differences.

${\Delta p}_{frict}=-\sum_{i}\hairsp\frac{\rho_i}{2}\frac{\lambda_T}{d_i}l_iv_i\left|v_i\right|$

where

Δp_frict	pressure change due to friction against walls (Pa)	dPFrict
λ_T	Darcy friction factor for tunnel wall	lambda
ρ_i	mean density of air in tunnel segment i (kg/m³)	Rho
d_i	hydraulic diameter of tunnel segment i (m)	dTun
l_i	length of tunnel segment i (m)	lTun
v_i	mean air velocity in segment i, positive from left to right (m/s)	VAir

${\Delta p}_{fricv}=0$

where

Δp_fricv		pressure change due to friction against vehicles (Pa)

${\Delta p}_{piston}=\sum_{i}\hairsp\sum_{j}\hairsp\frac{\rho_i}{2}\frac{A_{V_-j}}{A_i}\frac{c_{D{FV}_-j}}{\left(1-A_{V_-j}/A_i\right)^2}N_{\mathrm{front_\ }ij}\left(v_{V_-ij}-v_i\right)\left|v_{V_-ij}-v_i\right|$

where

Δp_piston	piston pressure rise (Pa)	dPPiston
c_{DFV_j}	drag coefficient at front end of vehicle type j	cdFront
N_{front_ij}	# of front ends of vehicle type j in segment i	nFront
A_{V_j}	cross section area of vehicle type j (m²)	frontArea
A_i	tunnel area of segment i (m²).	aTun

Platform (Platform)

In the platform the mass conservation is given by

$V\cdot\frac{d\rho_{dry}}{dt}=\sum_{j} m_{dot_j}+m_{dot_{leak}}$

where

ρ_dry	density of dry air = ρ ·(1 + w)^-1 (kg/m³)	RhoDry
ρ	density of air (kg/ m³)	Rho
w	humidity ratio of air	xHum
m_{dot_j}	mass flow from tunnel j (kg/s)	tunnel.m_dot
m_{dot_leak}	mass flow from leak (kg/s)	leak.m_dot
V	volume of platform (m³).	VEncl

The static pressure is given by

p= ρ·R·(T+273.15)

where

T		air temperature in platform (˚C).	TAir

Plenum (Plenum)

In the Plenum the mass conservation is given by

$V\cdot\frac{d\rho_{dry}}{dt}=\sum_{j} m_{dot_j}$

where

ρ_dry	density of dry air =  ·(1 + w)^-1 (kg/m³)	RhoDry
ρ	density of air (kg/ m³)	Rho
w	humidity ratio of air ( = 0 in IDA RTV)	xHum
m_{dot_j}	mass flow from tunnel j (kg/s)	tunnel.m_dot
V	volume of Plenum (m³).	VEncl

The static pressure is given by

p= ρ·R·(T+273.15)

where

T		air temperature in Plenum (˚C).	TAir

Tunnel entry (Entry)

The tunnel entry is modeled as a flanged entrance for airflow moving into the tunnel and an abrupt exit for airflow moving out of the tunnel. The reduction of wind pressure due to a vehicle in the tunnel entry is neglected.

Assuming given total pressure p_out outside the entry the Bernoulli equations become

$p_{out}=p_{tunnel}+\frac{\rho_{mean}}{2}\cdot(1+k_c)\cdot v_{tunnel}^2$

for airflow entering the tunnel and

$p_{out}=p_{tunnel}+\frac{\rho_{mean}}{2}\cdot(1-k_d)\cdot v_{tunnel}^2$

for airflow leaving the tunnel,

where

p_out	p_{height +}p_wind= total pressure outside tunnel relative reference level (Pa)	P1
p_tunnel	static pressure inside tunnel relative reference level (Pa)	tunnel.P
p_height	$p_{amb}\cdot(e^{-g\cdot h/(R\cdot(T_{amb}+273.15))}-1)$ ambient static pressure relative reference level (Pa)	PHeight
p_wind	wind pressure outside tunnel (Pa)	PWindCalc
p_amb	ambient pressure at reference level (Pa)	amb.P
ρ_mean	mean air density at entry (kg/m³)	Rho
k_c	loss coefficient for ingoing flow (into tunnel), default = 0.5	kPos
k_d	loss coefficient for outgoing flow (out of tunnel), default = 1	kNeg
v_tunnel	air velocity in tunnel inside entrance (m/s)	VAir
T_amb	ambient temperature (˚C)	amb.T
h	height relative reference level (m)		hTun
g	9.80665 = acceleration of gravity (m/s²)
R	287 = gas constant for air (J/kg/K).

Tunnel portal (Portal)

Openings of tunnels to stations/plenums are modeled in the same way as tunnel entries regarding airflows. For airflow entering the tunnel Bernoulli equation become

$p_{station}=p_{tunnel}+\frac{\rho_{mean}}{2}\cdot(1+k_c)\cdot v_{tunnel}^2-\mathrm{\Delta}p_{height}$

and for airflow entering the station/plenum

$p_{station}=p_{tunnel}+\frac{\rho_{mean}}{2}\cdot(1-k_d)\cdot v_{tunnel}^2-\mathrm{\Delta}p_{height}$

where

p_station	total pressure outside tunnel at floor level relative reference level (Pa)	station.P
p_amb	ambient at reference level (Pa)	amb.P
Δp_height	${(p}_{station}+\ p_{amb})\cdot(e^{-g\cdot h/(R\cdot(T+273.15))}-1)$ static pressure relative station/plenum floor level (Pa)	dPHeight
ρ_mean	mean air density at portal (kg/m³)	Rho
k_c	loss coefficient for flow entering tunnel, default = 0.5	kPos
k_d	loss coefficient for flow leaving tunnel, default = 1	kNeg
v_tunnel	air velocity in tunnel (m/s)	VAir
T	temperature in station/plenum (˚C)	station.T
h	height to tunnel from station floor (m)		hTun
g	9.80665 = acceleration of gravity (m/s²)
R	287 = gas constant for air (J/kg/K).

Recirculation (Recirc)

A model for the situation where two tunnel entries are so close that airflowing out from one tunnel can re circulate into the other has been formulated. The model has two parameters, maximum percentage flowing from each entry to the other. The mass flows from one tunnel to the other is given by

$m_{dot_21}=max{(}0,min{(}-m_{dot_1},m_{dot_2}\cdot f_{21}))$

$m_{dot_12}=max{(}0,min{(}-m_{dot_2},m_{dot_1}\cdot f_{12}))$

and from ambient by

$m_{dot_a1}=max{(}0,-m_{dot_1}-m_{dot_21})$

$m_{dot_a2}=max{(}0,-m_{dot_2}-m_{dot_12})$

where

m_{dot_}₁	mass flow in tunnel 1 (positive out of tunnel) (kg/s)	tunnel1.m_dot
m_{dot_}₂	mass flow in tunnel 2 (positive out of tunnel) (kg/s)	tunnel2.m_dot
m_{dot_}₂₁	mass flow recirculating from tunnel 2 to tunnel 1 (kg/s)	m_dot2_1
m_{dot_}₁₂	mass flow recirculating from tunnel 1 to tunnel 2 (kg/s)	m_dot1_2
m_{dot_a}₁	mass flow from ambient to tunnel 1 (kg/s)	m_dota_1
m_{dot_a}₂	mass flow from ambient to tunnel 2 (kg/s)	m_dota_2
f₂₁	maximum fraction of mass flow recirculating from tunnel 2 to tunnel 1	part2_1 [%]
f₁₂	maximum fraction of mass flow recirculating from tunnel 1 to tunnel 2.	part1_2 [%]

Tunnel junction (Branch)

The models for the converging flows and diverging flows in the branching component are the ones derived by [Gardel 1957], used in [Wallin 1994], and valid for a wye consisting of a straight uniform section with a connected branch, with a possibly different cross section area.

Shape1

One can distinguish between six cases of the flow pattern, four of which are covered by the original model by Gardel: converging flows (2, 3) to 1 and (3, 1) to 2 and diverging flows 1 to (2, 3) and 2 to (3, 1). The remaining two cases, converging to branch 3 and diverging from branch 3, are not described there but have been estimated by continuity to the other four cases and fitting to measurements (see below).

Denoting properties of flows at different locations by subscripts i, j and l, Bernoulli’s equations for the flow paths are written as follows.

Converging flow, (i, j) to l:

$p_i+\frac{\rho_i}{2}\cdot v_i^2=p_l+\frac{\rho_l}{2}\cdot(1+k_{cil})\cdot v_l^2$

$p_j+\frac{\rho_j}{2}\cdot v_j^2=p_l+\frac{\rho_l}{2}\cdot(1+k_{cjl})\cdot v_l^2$

Diverging flow, i to (j, l):

$p_i+\frac{\rho_i}{2}\cdot(1-k_{dij})\cdot v_i^2=p_j+\frac{\rho_j}{2}\cdot v_j^2$

$p_i+\frac{\rho_i}{2}\cdot(1-k_{dil})\cdot v_i^2=p_l+\frac{\rho_l}{2}\cdot v_l^2$

where

p_i	static pressure at location i (Pa)	tun<i>.P, i=1,2,3
ρ_i	p_i/(R·(T+273.15)) = density at location i (Pa)	Rho<i>
v_i	velocity at location i (positive inwards) (m/s)	VAir<i>
k_cil	loss coefficient for converging flow from i to l	K<il>
k_dil	loss coefficient for diverging flow from i to l	K<il>
T	temperature (˚C).	T

Depending on the flows the loss coefficients are computed as follows:

Converging flow, (2, 3) to 1, case 1:

$k_{c21}=0.03(1-\frac{Q_3}{Q_1})^2+(2-\frac{A_3}{A_1})\frac{Q_3}{Q_1}(1-\frac{Q_3}{Q_1})-(\frac{Q_3}{Q_1})^2[1+(1.62-\sqrt r)(\frac{cos{(}\pi-\phi)}{A_3/A_1}-1)-0.38(1-\frac{A_3}{A_1})]$ K21

$k_{c31}=-0.92(1-\frac{Q_3}{Q_1})^2+(2-\frac{A_3}{A_1})\frac{Q_3}{Q_1}(1-\frac{Q_3}{Q_1})-\ (\frac{Q_3}{Q_1})^2[(1.2-\sqrt r)(\frac{cos{(}\pi-\phi)}{A_3/A_1}-1)-(1-\frac{A_3}{A_1})\frac{cos{(}\pi-\phi)}{A_3/A_1}+0.8(1-\frac{1}{(A_3/A_1)^2})]$ K31

Converging flow, (3, 1) to 2, case 2:

k_c₁₂= formula for k_c₂₁ with subscripts 1 and 2 exchanged K12

and π - φ replaced by φ

k_c₃₂= formula for k_c₃₁ with subscript 1 replaced by 2 and π - φ by φ K32

Diverging flow, 1 to (2, 3), case 4:

$k_{d12}=0.03(1-\frac{Q_3}{Q_1})^2-0.2\frac{Q_3}{Q_1}(1-\frac{Q_3}{Q_1})+0.35(\frac{Q_3}{Q_1})^2$ K12

$k_{d13}=0.95(1-\frac{Q_3}{Q_1})^2+0.4\frac{Q_3}{Q_1}(1-\frac{Q_3}{Q_1})(1+\frac{A_1}{A_3})cot{(}\frac{\phi}{2})+(\frac{Q_3}{Q_1})^2(1-0.9\sqrt{\frac{r}{A_3/A_1}})[1.3cot{(}\frac{\phi}{2})-0.3+\frac{0.4-0.1A_3/A_1}{(A_3/A_1)^2}]$ K13

Diverging flow, 2 to (3, 1), case 5:

k_d₂₁= formula for k_d₁₂ with subscripts 1 and 2 exchanged K21

k_d23= formula for k_d₁₃ with subscript 1 replaced by 2 and φ by π - φ K23

where

Q_i	absolute value of volume flow at location i (m³/s)	\|VFAir<i>\|
A_i	cross section area at location i (A₁= A₂) (m²)	aTun<il>
φ	angle between straight section 1 and branch 3, φ [π/4, 3π/4]	Fi [degree]
r	radius of rounding at branch to main / diameter of main.	r

To estimate the remaining loss coefficients, i.e. converging to branch 3 (case 3) and diverging from branch 3 (case 6), we utilize that the solutions should be continuous when flows change direction. This is the case for transitions between the above described cases. E.g. the transition between converging flow case 2 ((3, 1) to 2) and diverging flow case 4 (1 to (2, 3)), i.e. when Q₃ = 0 we get by combining the above Bernoulli’s equations that

$\frac{\rho_l}{2}\cdot\left.k_{d12}\right|_{Q_3=0}=\frac{\rho_2}{2}\cdot\left.k_{c12}\right|_{Q_3=0}$

$\frac{\rho_l}{2}\cdot\left.k_{d13}\right|_{Q_3=0}=\frac{\rho_2}{2}\cdot(\left.k_{c12}\right|_{Q_3=0}-\left.k_{c32}\right|_{Q_3=0})$

should be satisfied for continuity. If the densities are equal, ρ₁= ρ₂, we see by inspecting the expressions for these loss coefficients that this is the case.

Now we require the same continuity properties for transitions where case 3 and case 6 (converging to and diverging from branch 3) are involved. From transition case 3 – case 4 (Q₂ = 0) we get assuming that ρ₁= ρ₃

$\left.k_{c13}\right|_{Q_2=0}=(\frac{A_3}{A_1})^2\cdot\left.k_{d13}\right|_{Q_2=0}\ =(\frac{A_3}{A_1})^2(1-0.9\sqrt{\frac{r}{A_3/A_1}})[1.3cot{(}\frac{\phi}{2})-0.3+\frac{0.4-0.1A_3/A_1}{(A_3/A_1)^2}]$

$\left.k_{c23}\right|_{Q_2=0}=\left.k_{c13}\right|_{Q_2=0}-(\frac{A_3}{A_1})^2\left.k_{d12}\right|_{Q_2=0}=\left.k_{c13}\right|_{Q_2=0}-(\frac{A_3}{A_1})^20.35$

and from transition case 3 – case 5 (Q₁ = 0) we get assuming that ρ₁= ρ₃

k_c₂₃|_Q1=0 = formula for k_c₁₃|_Q2=0 with subscripts 1 and 2 exchanged and φ replaced by π - φ

k_c₁₃|_Q1=0 = formula for k_c₂₃|_Q2=0 with subscript 2 replaced by 1.

In an analogous manner we derive coefficient for case 6 (3 to (1, 2)). From transition case 2 – case 6 (Q₁ = 0) we get assuming that ρ₁= ρ₃

$\left.k_{d32}\right|_{Q_1=0}=(\frac{A_3}{A_2})^2\cdot\left.k_{c32}\right|_{Q_1=0}\ =-(\frac{A_3}{A_2})^2[(1.2-\sqrt r)(\frac{cos{(}\phi)}{A_3/A_2}-1)-(1-\frac{A_3}{A_2})\frac{cos{(}\phi)}{A_3/A_2}+0.8(1-\frac{1}{(A_3/A_2)^2})]$

$\left.k_{d31}\right|_{Q_1=0}=\left.k_{d32}\right|_{Q_1=0}-(\frac{A_3}{A_2})^2\left.k_{c12}\right|_{Q_1=0}=\left.k_{d32}\right|_{Q_1=0}+(\frac{A_3}{A_2})^2[1+(1.62-\sqrt r)(\frac{cos{(}\phi)}{A_3/A_2}-1)-0.38(1-\frac{A_3}{A_2})]$

and from transition case 1 – case 6 (Q₂ = 0) we get assuming that ρ₂= ρ₃

k_d₃₁|_Q2=0 = formula for k_d₃₂|_Q1=0 with subscripts 1 and 2 exchanged and φ replaced by π - φ

k_d₃₂|_Q2=0 = formula for k_d₃₁|_Q1=0 with subscript 2 replaced by 1 and φ by π - φ.

The loss coefficients of the remaining two flow cases can now be approximated.

Converging flow, (1, 2) to 3, case 3:

$k_{c13}=\left.k_{c13}\right|_{Q_2=0}\frac{Q_1}{Q_3}+\left.k_{c13}\right|_{Q_1=0}(1-\frac{Q_1}{Q_3})$ K13

$k_{c23}=\left.k_{c23}\right|_{Q_2=0}\frac{Q_1}{Q_3}+\left.k_{c23}\right|_{Q_1=0}(1-\frac{Q_1}{Q_3})$ K23

Diverging flow, 3 to (1, 2), case 6:

$k_{d31}=\left.k_{d31}\right|_{Q_1=0}\frac{Q_2}{Q_3}+\left.k_{d31}\right|_{Q_2=0}(1-\frac{Q_2}{Q_3})$ K31

$k_{d32}=\left.k_{d32}\right|_{Q_1=0}\frac{Q_2}{Q_3}+\left.k_{d32}\right|_{Q_2=0}(1-\frac{Q_2}{Q_3})$ K32

When the cross-sectional area changes between segments 1 and 2, the abovementioned equations are modified by replacing A₁ or A₂ with 0.5(A₁+ A₂). This is a modification to Gardel’s equations.

Supply or exhaust fan to ambient (FanDuct)

The supply or exhaust duct is an extension of the branch with the model described below connected to branch 3 and to ambient. The model consists of a fan and optional damper. At no external Δp and f_cntr = 1, the given flow (V_{f_Nom}) is obtained. When f_cntr = 0 and DamperAtFanOff = Open, air will flow according to external Δp (including stack effect) with a system characteristic through the point p_head and V_{f_Nom}. When f_cntr = 0 and DamperAtFanOff = Closed, there is no airflow.

$\mathrm{\ 10}\frac{dV_f}{dt}+V_f=V_{f_Nom}sign(\mathrm{\Delta p})\cdot\sqrt{|\mathrm{\Delta p}|/|p_{head}|}$ if f_cntr > 0 or DamperAtFanOff = Open

=0 otherwise

where

V_f	volume flow, negative for supply and positive for exhaust (m³/s)	VFAir
V_{f_Nom}	nominal volume airflow at no external pressure difference (m³/s)	VFNom
Δp	$p_{tun}+p_{rise}-p_{height}-f_{cntr}^2\cdot p_n\cdot\rho/\rho_{nom}$ = total static pressure difference	dPtot
p_tun	static pressure at tunnel relative ambient pressure at reference level (Pa)	tunnel.P
p_rise	$-(p_{amb}+p_{height})\cdot(e^{-g\cdot(h_{tun}-h_{amb})/(R\cdot(T+273.15))}-1)$ = static pressure rise in duct (Pa)	PRise
p_height	$p_{amb}\cdot(e^{-g\cdot h_{amb}/(R\cdot(T_{amb}+273.15))}-1)$ = ambient static pressure at opening relative ambient pressure at reference level (Pa)	PHeight
p_n	if supply p_head else -p_head (Pa)	p_n
p_head	static fan pressure at nominal flow (Pa)	pHead
f_cntr	fan capacity control, 0 for no fan dp, 1 for p_head (Pa)	fanCntr
p_amb	ambient pressure at reference level (Pa)	amb.P
T_amb	ambient temperature (˚C)	amb.T
h_amb	level of ambient opening (grille) relative to reference level (m)	hAmb
h_tun	tunnel level relative to reference level (m)	hTun
ρ	if V_f < 0 then ρ_amb else ρ_tunnel = air density (kg/m³)	Rho
T	if V_f < 0 then ρ_amb else ρ_tunnel = air temperature (˚C)	TAir
T_tunnel	temperature at branch to tunnel (˚C)	tunnel.T
ρ_nom	air density at which the fan performance is defined (kg/m³)	RhoNom
g	9.80665 = acceleration of gravity (m/s²)
R	287 = gas constant for air (J/kg/K)
Q_sup	$V_f\cdot\|f_{cntr}^2\cdot p_n\cdot\rho/\rho_{nom}\|/e_{ff}$ = power supplied(W)	PowerSup kW]
e_ff	fan efficiency.	effic

Leak (Leak)

In a leak the relation between mass flow and pressure difference is given by

$\mathrm{\Delta p}\cdot(C_d\cdot A\cdot f_{open})^2=1/(2\cdot\rho)\cdot|m_{dot}|\cdot m_{dot}$

where

m_dot	mass flow to ambient (kg/s)	station.m_dot
C_d	discharge coefficient	cd
A	opening area (m²)	ar
f_open	opening control, 1 for fully open and 0 for closed	partOpen
ρ	mean air density in leak (kg/m³)	Rho
Δp	$p_{platf}+p_{rise}-p_{height}-p_{wind}$ = total static pressure difference (Pa)	dPtot
p_platf	static pressure at platform relative ambient pressure at reference level (Pa)	station.P
p_rise	-g ·ρ_station· h_leak = stack pressure rise to leak (Pa)
p_height	$p_{amb}\cdot(e^{-g\cdot(h_{floor}+h_{leak})/(R\cdot(T_{amb}+273.15))}-1)$ = ambient static pressure at leak relative ambient pressure at reference level (Pa)	PHeight
p_wind	wind pressure outside leak (Pa)	PWindCalc
h_leak	leak height from floor of platform (m)	hLeak
h_floor	platform level relative to reference level (m).	hFloor

Given volume airflow (VFGiven)

The volume airflow is specified by a parameter V_{f_Nom} and achieved through a differential equation with a time constant (10 sec.)

$\tau\cdot\frac{dm_{dot}}{dt}+m_{dot}=f_{cntr}\cdot\rho\cdot V_{f_Nom}$

where

m_dot	mass flow (kg/s)	inlet.m_dot
V_{f_Nom       }	specified volume airflow ( m³/s)	VFNom
f_cntr	control signal in interval [-1 : 1], -1 for reverse action	cntr
ρ	air density (kg/ m³)	Rho_tun1
τ	time constant (s).	tau

Filter (Filter)

The volume airflow is given by an input variable V_{f_Filter} and achieved through a differential equation with a time constant (10 sec.)

$\ 10\cdot\frac{dm_{dot}}{dt}+m_{dot}=\rho\cdot V_{f_{Filter}}$

where

m_dot	mass flow (kg/s)	inlet.m_dot
V_{f_Filter}	specified volume airflow ( m³/s)	VFFilter
ρ	air density (kg/ m³).	Rho_tun1

Damper (Damper)

In the damper the pressure difference is described as a function of the mass flow. At fully open position the flow is governed by a loss coefficient and at closed by an Equivalent Leakage Area. The damper can be defined as asymmetric. In that case a different Equivalent Leakage Area is given for air flow in the opposite direction.

$\mathrm{\Delta p}=sign(m_{dot})\cdot|m_{dot}/c_{Act}|^{1/n_{Act}}$

where

m_dot	mass flow (kg/s)	left.m_dot
n_Act	if powerLawLeakage then $0.5+f_{closed}^k(n-0.5)$ else 0.5	nAct
c_Act	if m_dot > 0 then $c_{open}+f_{closed}^k(c_{closed}-c_{open})$ else $c_{open}+f_{closed}^k(c_{closR2L}-c_{open})$	cAct
f_closed	- τ_open ·df_closed/dt + (1 - f_openSig) = delayed opening/closing	partClosed
f_openSig	opening control, 1 for fully open and 0 for closed	partOpen
τ_open	time constant for opening/closing (s)	tau
n	power law exponent at closed	n
c	power law coefficient at closed	c
k	integer exponent in linear combinations	k
c_open	$\sqrt{2\cdot\rho_{air}/k_{loss}}\cdot A$	c_open
c_closed	if powerLawLeakage then c else $\sqrt{2\cdot\rho_{air}}\cdot A_{el}\cdot C_d$	c_closed
c_closR2L	if powerLawLeakage or symmetric then c_closed else $\sqrt{2\cdot\rho_{air}}\cdot A_{elR2L}\cdot C_d$	c_closR2L
C_d	discharge coefficient	cd
Δp	static pressure difference between left and right (Pa)	dPtot
ρ_air	density of air (kg/m³)	Rho_air
k_loss	loss coefficient at fully open	kLoss
A	cross section area (m²)	tArea
A_el	equivalent leakage area (m²)	ela
A_elR2L	equivalent leakage area for air flow from right to left if not symmetric (m²).	elaR2L

Passage (Passage)

A model for parallel and/or serial doors for passenger passage has been implemented. Both stochastic and intermittent passenger flows are handled. We assume that a door section consists of a user defined number of doors and in the case of two serial door sections that these are far enough apart so that in the statistical treatment they can be considered as independent.

The relation between pressure difference and air mass flow through a passage is described by the equation

$\mathrm{\Delta p}=sign(m_{dot})\cdot|m_{dot}/c_{Act}|^2$

where

Δp	static pressure difference between left and right (Pa)	dPtot
m_dot	mass flow, positive from left to right (kg/s)	left.m_dot
c_Act	coefficient (see below).	cAct

For stochastic passenger flows we assume that the number of open doors in a door section is binomially distributed, i.e.

$f(i,p,n)=\left(\begin{matrix}n\\i\\\end{matrix}\right)\cdot p^i\cdot(1-p)^{n-i}$

where

f(i,p,n)		probability that i doors are open
n		number of doors in door section
p		time fraction when a door is open.

For parallel doors (a left door section) we get

$c_{ActS}=\mathbit{f}_L^T\cdot\mathbit{c}_L=(1-p_L)\cdot c_{closedL}+p_L\cdot c_{openL}$

where

c_ActS	the value of c_Act in case of stochastic passenger flow
f_L	vector with components f(i, p_LnL_{) , i = 0, n}L
n_L	number of doors in left door section	nl
c_L	vector [ c_closedL ,{ c_openL· i/ n_L , i = 1, n_L }]^T
p_L	time fraction when a door in left section is open (see below)	fOpenL
c_closedL	$\sqrt{2\cdot\rho_{air}}\cdot A_{elL}\cdot C_d$ = coefficient when all left section doors closed
c_openL	$\sqrt{2\cdot\rho_{air}/k_{lossL}}\cdot A$ = coefficient when all left section doors open
C_d	discharge coefficient	cd
ρ_air	density of air (kg/m³)	Rho_air
k_lossL	loss coefficient at fully open left door section	kLossL
A	cross section area (m²)	tArea
A_elL	equivalent leakage area in left door section (m²).	ela

In the case with serial door sections (left and right door section) the expression for c_ActS becomes, with an analogous notation for the right door section (R for right)

$c_{ActS}=\mathbit{f}_L^T\cdot\mathbit{C}_{LR}\cdot\mathbit{f}_R$

where

C_LR		matrix with $C_{LR}(i,j)=1/\sqrt{c_L^{-2}(i)+c_R^{-2}(j)}$ , i = 0, n_L, j = 0, n_R
f_R		vector with components f(i, p_RnR_{) , i = 0, n}R
c_R		vector [ c_closedR ,{ c_openR· i/ n_R , i = 1, n_R }]^T .

The probability that a door is open (time fraction when a door is open) is given by

$p_L=1-e^{-\tau_L\cdot q_{stoch}\cdot s_\tau/3600/n_L}$

$p_R=1-e^{-\tau_R\cdot q_{stoch}\cdot s_\tau/3600/n_R}$

where

τ_L	door opening time for single passenger in left door section (s)	tauL
τ_R	door opening time for single passenger in right door section (s)	tauR
q_stoch	q_sL + q_sR stochastic passenger flow at rush hour (passengers/h)	passStoch
q_sL	stochastic passenger flow from left at rush hour (passengers/h)	passLeftStoch
q_sR	stochastic passenger flow from right at rush hour (passengers/h)	passRightStoch
s_τ	- τ ·ds_τ /dt + s = delayed passenger flow fraction	passFrtau
s(t)	passenger flow fraction relative rush hour passenger flow	passFract
τ	time constant (s)	tau
t	time (s).

For intermittent passenger flow and parallel doors (left door section) we have

$c_{ActP}=c_{openL}$

$p_{pulse}=min(1,q_{pulse}\cdot s_\tau/(q_{max}-q_{stoch}\cdot s_\tau))$

where

c_ActP	the value of c_Act during the intermittent passenger flow
p_pulse	time fraction when all door sections are open for intermittent passenger flow	fOpenPulse
q_pulse	q_pL + q_pR = intermittent passenger flow at rush hour (passengers/h)	passPulse
q_pL	intermittent passenger flow from left at rush hour (passengers/h)	passLeftPulse
q_pR	intermittent passenger flow from right at rush hour (passengers/h)	passRightPulse
q_max	user-defined passage capacity (passengers/h).	passMax

The equation for p_pulse takes into account that the intermittent flow is augmented by the simultaneous stochastic flow.

For serial door sections (left and right door sections) the expression for c_ActP becomes

$c_{ActP}=1/\sqrt{c_{openL}^{-2}+c_{openR}^{-2}}$ .

The coefficient c_Act for combination of the stochastic and intermittent passenger flows is evaluated as the mean value of c_ActS and c_ActP:

$c_{Act}=(1-p_{pulse})\cdot c_{ActS}+p_{pulse}\cdot c_{ActP}$ .

Axial fan (AxialFan)

The axial fan is described by a fan curve given as a spline, where pressure rise is given as a function of air volume flow

$\ \mathrm{\Delta p}=f_{spline}(V_f/f_{vAct})\cdot|f_{vAct}|\cdot f_{vAct}\cdot\rho/\rho_{nom}$

where

Δp	static pressure difference between inlet and outlet ( Pa)	dPtot
V_f	air volume flow ( m³/s)	VFAir
f_spline(V_f)	spline giving pressure rise as a function of air volume flow (Pa)	spline
f_vAct	- τ ·df_vAct/dt + f_vRel = delayed f_vRel	speedAct
f_vRel	quotient fan speed to nominal fan speed in interval [0 : 1]	speedRel
τ	time constant (s)	tau
ρ	air density (kg/m³)
ρ_nom	air density at which the fan performance is defined (kg/m³)	RhoNom
Q_sup	max(0, V_{f ;}·Δp )/ e_ff = power supplied(W)	PowerSup [kW]
e_ff	e_spline(V_f/ f_vAct) = fan efficiency	effic
e_spline(V_f)	spline giving fan efficiency as a function of air volume flow.	spline

Saccardo nozzle (SaccarNoz)

In the Saccardo nozzle the momentum of the supplied air is taken into account. Furthermore the pressure at location (3) of the supplied air is assumed to be the same as at the incoming (1).

Shape1

$(p_1-p_2)\cdot A=m_{dot_2}\cdot v_2-m_{dot_1}\cdot v_1-m_{dot_3}\cdot v_3\cdot c o s{(}\phi)\cdot k_{mx}$

p_3=p_1

where

m_{dot_i}	mass flow at location i (positive inward) (kg/s)	tuni.m_dot
v_i	velocity at location i (m/s)	VAiri
p_i	static pressure at location i (Pa)	tuni.P
φ	angle between nozzle and main (rad)	Fi [degree]
k_mx	momentum exchange coefficient	Kmx
A	cross section area of main (m²).	aTun12