0 Votes

Other models

Last modified by omid on 2025/03/27 10:20

Previous (Air quality) | Back to Start | Next (References)

TABLE OF CONTENT

Other models
Tunnel water flow model: condensation, evaporation, seepage and drainage
Critical velocity
Mold growth index
Ice build-up
Train motion
Fire throttling effect
Back Layering length for tunnel fire under longitudinal ventilation
Inventory of models

Other models

Tunnel water flow model: condensation, evaporation, seepage and drainage

IDATunnelTheory_work -2024-1-O2_html_5ac1e76ee645a636.gif

In IDA Tunnel, water may enter the tunnel via condensation on tunnel walls, seepage through the wall or through explicitly defined water or vapor sources (an option for tunnel Equipment). Wall seepage is considered to be homogeneously radial, flowing into the tunnel from all directions and the advected heat will influence the temperature profile into the tunnel wall.

Once inside the tunnel, entering water will form a free water table, which may evaporate into tunnel air. The size of the water table may be controlled by the user through the parameters DesignmDrain and fWet. Non-evaporated water is regarded to be immediately drained from each tunnel segment. When the flow, mDrain, of this water equals the given parameter DesignmDrain, the area of the water table is regarded to be fWet times the tunnel inner wall area. At smaller flows, the wetted area tends linearly to zero as mDrain approaches zero. At flows larger than DesignmDrain, the wetted area stays constant. To define a wetted area which is independent of actual flow, select DesignmDrain to be considerably smaller than actual mDrain (but non-zero).

The drain flow gets three contributions

$m_{drain}=m_{seep}+m_{cond}+m_{condAir}$

where

m_drain	drain flow per tunnel length (kg/sm)	mDrain
m_seep	water seepage per tunnel length out of wall (kg/sm)	mSeep
m_cond	condensation/evaporation flow per tunnel length at wall (kg/sm)	Humf
m_condAir	condensation in air per tunnel length (not in platform) (kg/sm).	mAirCond

$m_{cond}=U_{conv}/c_{pAir}\cdot p_{act}\cdot f_{mol}\cdot(C-C_{surf})$

where

C	humidity fraction (m³/m³)
C∙ f_mol	humidity ratio (kg/kg)	xHum
C_surf	humidity fraction on wall (m³/m³)
C_surf f_mol	surface humidity ratio (temperature and pressure dep.) (kg/kg)	xHumSurf
U_conv	convective heat transfer coefficient (see sec. 2.2.1) (W/ m²˚C)	UConv
c_pAir	specific heat of air (J/kg˚C)	CP_AIR
p_act	actual wet perimeter (m)
f_mol	0.62198 = ratio of molecular masses of water and air.

$m_{condAir}=c_{cond}\cdot f_{mol}\cdot A\cdot m a x{(}0,C-C_{sat})$

where

C_sat	saturated humidity fraction in air (m³/m³)
C_sat∙f_mol	saturated humidity ratio (kg/kg)	xHumSat
c_cond	coefficient for air condensation (kg/sm³)	CondensCoeff
A	tunnel cross section area of (m²).	aTun

The actual perimeter used p_actin the computations of m_cond is

p_act		if C > C_surf then p else p_wet (m)

where p_wet is given by

$\tau_{wet}\cdot\frac{dp_{wet}}{dt}+p_{wet}=f_{wet}\cdot p\cdot m i n{(}\frac{m_{drain}}{m_{design}},1)$ if m_drain > 0

= 0 otherwise

p_wet	wet perimeter (m)	wetPerim
p	physical perimeter (m)
τ_wet	time constant (s)	TauWet
f_wet	fraction of perimeter which is wet at design drain flow	fWet
m_design	design drain flow/meter (kg/sm).	designmDrain

Critical velocity

The minimum steady-state velocity of the ventilation air moving towards the fire, the critical velocity, needed for removing the smoke and preventing backlayering is computed according to the Kennedy formula. In the tunnel section models this is done by solving the system

$V_C=K_1\cdot K_g\cdot\left(\frac{g\cdot h\cdot Q_{fire}}{\rho_{appr}\cdot c_{pAir}\cdot A\cdot(T_f+273.15)}\right)^{1/3}$

$T_f=\frac{Q_{fire}}{\rho_{appr}\cdot c_{pAir}\cdot A\cdot V_C}+T_{appr}$

where

V_C	critical velocity (m/s)	VC
T_f	average temperature of fire site gases (˚C)	Tf
T_appr	temperature of approaching air (˚C)
Q_fire	heat from fire directly to air (W)
ρ_appr	density of approaching air (kg/m³)
c_pAir	specific heat of air (J/kg˚C)	CP_AIR
A	cross section area of fire segment (m²) (NFPA 502-2014 and 502-2017) annulus area, not including the area blockage, such as vehicles (m²) (NFPA 502-2023)
h	tunnel height in fire segment, default is hydraulic diameter but if parameter VCHeight set > 0 this value is used (m) (NFPA 502-2014 and 502-2017) height from the base of the fire to the tunnel ceiling at the fire, default is hydraulic diameter but if parameter VCHeight set > 0 this value is used (m) (NFPA 502-2023)
g	9.80665 = acceleration of gravity (m/s²)
K₁	=Froude number factor = Fr_c^-1/3 is a function of Q_fire according to NFPA 502-2017. =Froude number factor = 4.5^-1/3= 0.606 according to NFPA 502-2014 and NFPA 502-2023.
Fr_c	Froude number.
K_g	1 + 0.0374 ∙ \|min(0, f_slope)\|^0.8 = grade factor	Kg
f_slope	slope in direction of moving air (%).

The subject of critical velocity remains an area of ongoing research, and the suggested formulas have been revised in every published edition of NFPA 502 over the last decade. These formulas will likely change again in future editions of NFPA 502. Therefore, it is the user’s responsibility to stay informed about the applicability and limitations of the various critical velocity formulas. Some guidance on their applicability and limitations can be found in Annex D of the NFPA 502-2023 version.

The critical velocity according to (NFPA 502-2014) is presented by VC0 in the list of variables.

The critical velocity according to (NFPA 502-2017) is presented by VC in the list of variables.

The critical velocity according to (NFPA 502-2023) is presented by VC1 in the list of variables.

As output signal we compute the quotient between velocity of approaching air v_appr (with sign) and critical velocity V_Cand VC1

V_CSignal	=	v_appr / V_C .	VCSignal

V_CSignal1	=	v_appr / V_C1 .	VCSignal1

Mold growth index

A mathematical model for mold growth on wooden materials [Hukka, Viitanen 1999] has been implemented. The growth is described by the mold index, which is used in experiments for visual inspections and is quantified as follows:

no growth
some growth detected only with microscope
moderate growth detected with microscope (cover more than 10 %)
some growth detected visually
visually detected coverage more than 10 %
visually detected coverage more than 50 %
visually detected coverage 100%.

The model consists of a differential equation (with constraints) where effects of exposure time, temperature, relative humidity and dry periods are taken into account. Our interpretation of the model is described below where the mold growth index is set equal to zero initially.

if T < 0.1 or RH < RH_crit then

$\frac{dM}{dt}=\left\{\begin{matrix}-0.032\\0\\-0.016\\\end{matrix}\begin{matrix}\\\\\\\end{matrix}\begin{matrix}\\\\\\\end{matrix}\begin{matrix}when\\"\\"\\\end{matrix}\begin{matrix}\\\\\\\end{matrix}\begin{matrix}t-t_1\le0.25\\0.25<t-t_1\le1\\1<t-t_1\\\end{matrix}\right.$

M = max(M, 0)

else

if M < 1then

$\frac{dM}{dt}=t_m^{-1}$

else

$\frac{dM}{dt}=t_m^{-1}\cdot k_1\cdot k_2$

endif

M = min(M, M_max)

endif

where

M	mold growth index	mold
t	time (days)
t₁	beginning of dry period (T < 0.1 or RH < RH_crit) (days)
T	air temperature (˚C)	TAir
RH	relative humidity (%)	relHum
RH_crit	–0.00267·T³ + 0.16·T² – 3.13·T + 100 if T < 20 else 80 (%) critical relative humidity favorable to initiation of mold growth
M_max	$1+7\cdot\frac{RH_{crit}-RH}{RH_{crit}-100}-2\cdot\left(\frac{RH_{crit}-RH}{RH_{crit}-100}\right)^2$ = largest possible mold growth
t_m	7·exp(–0.68·ln(T) – 13.9·ln(RH) +0.14·W – 0.33·SQ + 66.02) = response time for initiation of mold growth (days)
t_v	7·exp(–0.74·ln(T) – 12.72·ln(RH) + 0.06·W + 61.50) = response time for visual appearance of mold growth (days)
k₁	2/( t_v / t_m – 1) if M > 1 else 1
k₂	max(0, 1 – exp(2.3·(M – M_max )))
W	wood species (0 = pine (used), 1 = spruce)
SQ	surface quality (0 = resawn (used), 1 = original kiln-dried).

Ice build-up

Ice build-up is computed using an enthalpy formulation. We assume that the source water comes from the drain flow described in section 2.4. Further we assume that the ice forms half spheres at every Δx meter.

$\frac{dM_{ice}}{dt}=m_{watIn}-m_{watOut}$

$\frac{d(M_{ice}\cdot H_{ice})}{dt}=m_{watIn}\cdot H(T_{watIn})-m_{watOut}\cdot H(T_{watOut})+\alpha_{air}\cdot(T_{air}-T_{ice})+\alpha_{wall}\cdot(T_{wall}-T_{ice})$

where

M_ice	mass of ice per tunnel length ≥ 0 (kg/m)	mIce
m_watIn	m_drain · p_wet/( f_wet·p) = drain flow available for ice build-up per tunnel length ≥ 0 (see 2.4) (kg/sm)
m_drain	drain flow per tunnel length ≥ 0 (see 2.4) (kg/sm)	mDrain
p_wet	wet perimeter (m)	wetPerim
p	physical perimeter (m)
f_wet	fraction of perimeter which is wet at design drain flow	fWet
m_watOut	water leaving per tunnel length ≥ 0 (kg/sm)
H_ice	enthalpy of ice (J/kg)
T_ice	temperature of ice (˚C)	TIce
H(T_watIn)	enthalpy of drain flow (J/kg)
T_watIn	temperature of drain flow (˚C)
H(T_watOut)	enthalpy of leaving water (J/kg)
T_watOut	temperature of leaving water (˚C)
T_air	temperature of tunnel air (˚C)	TAir
T_wall	temperature of tunnel wall (˚C)	TWall
α_air	2·π·r²/(p·Δx ) · α = heat transfer coefficient between ice and air (W/˚C m)
α_wall	π·r²/((3·r/8)·Δx) · k_ice = heat transfer coefficient between ice and wall (W/˚C m)
α	heat transfer coefficient between wall and air (see 2.2.4) (W/˚C m)	Alfa
k_ice	2.2 = thermal conductivity of ice (W/˚C m)
Δx	distance between ice masses (m)	dxIce
r	(3· M_ice ·Δx/(2·π·ρ_ice))^1/3 = radius of spheres (m)
ρ_ice	917 = density of ice (kg/m³).

The enthalpy function H(T) is approximated with linear functions within the two phases and a jump at the melting temperature such that the enthalpy of ice at melting is zero:

dH/dT = c_pWat = 4187 (J/˚C kg) in the liquid phase (T > 0)

dH/dT = c_pIce = 2060 in the solid phase (T < 0)

ΔH = 3.34·10⁵(J/kg) latent heat of solidification/melting (T = 0).

We have two equations but five unknowns (M_ice, m_watOut, H_ice, T_watOut, T_watIn), so three assumptions are needed to get a well defined problem. Three cases can be defined. Either there is no ice (case 0), water and ice (case 1) or just ice (case 2). In case 0 we set M_ice = 0 and H_ice = 0. In the other two cases three assumptions are made.

T_watOut = 0
T_watIn = max(0, T_wall)
H_ice = 0 in case 1 and m_watOut = 0 in case 2.

In case 1 we get

$\frac{dM_{ice}}{dt}=m_{watIn}-m_{watOut}$

$0=m_{watIn}\cdot H(T_{watIn})-m_{watOut}\cdot\mathrm{\Delta H}+\alpha_{air}\cdot T_{air}+\alpha_{wall}\cdot T_{wall}$

to solve for M_ice and m_watOut and in case 2

$\frac{dM_{ice}}{dt}=m_{watIn}$

$\frac{d(M_{ice}\cdot H_{ice})}{dt}=m_{watIn}\cdot H(T_{watIn})+\alpha_{air}\cdot(T_{air}-T_{ice})+\alpha_{wall}\cdot(T_{wall}-T_{ice})$

for M_ice and H_ice(T_ice = H_ice/c_pIce).

Set $s=m_{watIn}\cdot H(T_{watIn})+\alpha_{air}\cdot T_{air}+\alpha_{wall}\cdot T_{wall}$

then the cases are executed as follows to generate continuous functions M_ice and H_ice:

if s < 0 and (M_ice > 0 or m_watIn > 0) then

case 2 (M_ice > 0, H_ice ≤ 0)

else if s < 0 or (M_ice = 0 and s > m_watIn·ΔH) then

case 0 (M_ice = 0, H_ice;= 0)

else

case 1 (M_ice ≥ 0, H_ice= 0)

endif.

Train motion

Our model for the train motion is based on the requirement that it can be computed before starting the simulation. This implies that some factors affecting the train motion, e.g. the drag force, are omitted. We believe these factors have a minor influence on the train motion. They are still computed during the simulation and affect the train power consumption. The main factors determining the motion in our model are:

Maximum allowed acceleration and retardation of train max_acc, max_ret
Maximum available power of train max_pow
Mass of train totalMass, Km
Rolling resistance of train A, B
Speed at nodes between route segments rSeg.v0, rSeg.v1
Maximum allowed speed in route segments rSeg.vMax
Coasting rSeg.vCoast
Route profile (slopes along route).

A route is divided into route segments. A default segmentation is defined by the presence of platforms (or sections marked to be platforms) along the route. However, the number of and size of segments can be manually edited (e.g., to limit the speed on some section). Route segments are defined by a number, start position, length, maximum speed, end speed and dwell time (at stations).

A starting train will initially accelerate according to max_acc. If the necessary power, increasing with speed, reaches max_pow then the acceleration will be limited accordingly. Slopes along the route also affect the necessary power and can thus prevent the train from reaching maximum allowed speed or decrease train speed.

If the train speed should decrease at the end of a segment (i.e. the end speed is lower than the actual maximum speed) the deceleration is modelled as follows. If coasting is allowed the train coasts until the speed decreases to the minimum coasting speed, then it starts breaking. If coasting is not allowed or the train speed is lower than the minimum coasting speed, the train is only breaking. The breaking is done with retardation max_ret. Slope and rolling resistance affect the coasting.

On each route segment the basic equations for determining the train motion are:

$\frac{dx_{veh}}{dt}=v_{veh}$

$\frac{dv_{veh}}{dt}=a_{veh}$

$F_{loc}=F_{drag}+F_{roll}+M_{veh}\cdot g\cdot s i n{(}f_{slope})+km_{veh}\cdot M_{veh}\cdot a_{veh}$

with constraints

$v_{veh}\le v_{max}$

$& a_{veh}\left\{\begin{matrix}\le a_{max}&acceleration\\=-r_{max}&retardation\\\end{matrix}\right.$

$F_{loc}\cdot v_{veh}\le P_{max}$

$F_{loc}=\begin{matrix}0&coasting\\\end{matrix}$

and boundary conditions

$v_{veh}=v_0$ at start of route segment

$v_{veh}=v_{end}$ at end of route segment

where

x_veh	position of vehicle (m/s)	train.x
v_veh	velocity of vehicle (m/s)	train.x_dot
a_veh	acceleration of vehicle (m/s²)	train.x_dot_dot
v_max	Maximum velocity in route segment (m/s)	rSeg.vMax
a_max	Maximum acceleration of vehicle (m/s²)	max_acc
r_max	Maximum retardation of vehicle (m/s²)	max_ret
P_max	Maximum available power of train (W)	max_power
v₀	velocity of vehicle at start of route segment (m/s)	rSeg.v0, rSeg.v1
v_end	velocity of vehicle at end of route segment (m/s)	rSeg.v1
F_loc	locomotive force (N)	LocF [kN]
F_drag	aerodynamic drag force (see 2.1.1) (N), ignored when calculating train motion	DragF [kN]
F_roll	(A_veh+ B_veh· v_veh) · M_veh = rolling resistance force (N) only used when coasting, otherwise ignored	RollF [kN]
g	9.80665 = acceleration of gravity (m/s²)
km_veh	kinetic energy coefficient	Km
f_slope	mean slope of tunnel	slope
M_veh	train total mass (incl. passengers) (kg)	totalMass
A_veh	coefficient (m/ s²)	A
B_veh	coefficient (s^-1).	B

Note that the drag force F_drag is neglected here and that the rolling resistance force F_roll is included only for coasting computations.

The conditions on the velocity at nodes between route segments need not be completely satisfied. If limiting acceleration and power constraints prevents the train velocity at ends of segments from reaching the specified value it is accepted if it is smaller than specified velocity. Larger end velocities are not allowed, e.g. a train has to stop at a station, and will cause an error message (except the first segment of the route where the error will be fixed by decreasing the start velocity).

The solution to this problem for the train motion in a route is saved in three parameter vectors containing time, position and speed.

Fire throttling effect

The throttling effect refers to a tunnel fire's tendency to impede airflow. Certain sources of aerodynamic resistance, influenced by temperature, contribute to an overall increase in tunnel resistance as the fire's heat release rate rises. This effect primarily stems from the expansion of heated air and combustion products downstream of the fire, resulting in decreased static pressure at the fire location due to temperature elevation, air expansion, and acceleration within the tunnel. Subsequently, the static pressure increases owing to airflow deceleration caused by cooling at the tunnel wall downstream of the fire. Additionally, there is an augmented airflow resistance due to temperature stratification downstream of the fire. This stratification effect varies depending on the case and can typically be accurately determined using 3D CFD simulations. Among these effects, the first, the pressure drop occurring in a fire, can be analytically calculated using momentum and mass balance for the fire control volume.

The pressure drop is computed based on the following equation:

${\mathrm{\Delta p}}_{fire}=\rho_c{u_c}^2\left(1-\frac{T_h}{T_c}\right)$

where c represents the upstream of the fire and h represents its downstream (Figure).

By selecting the auto calculation option, the fire pressure drop is automatically calculated using the above formula, which depends on the air velocity and fire power.

Back Layering length for tunnel fire under longitudinal ventilation

The concept of backlayering length was initially proposed in the context of a horizontal tunnel fire under longitudinal ventilation. In a longitudinally ventilated tunnel, if the ventilation velocity is smaller than the so-called critical velocity, buoyant smoke can oppose the forced longitudinal flow, eventually leading to the formation of a steady-state backlayering flow in the upstream region of the fire source.

The back layering length is computed using the following equation.

${L_b}^\ast=\frac{L_b}{H}=19.8\ ln(\frac{V_{cr}}{V})$

Here, V represents the air velocity approaching the fire, and Vcr denotes the critical velocity of the fire. H stands for the tunnel height, and if it is unspecified, the tunnel hydraulic diameter is utilized. Mathematically, the equation yields a negative value when the air speed exceeds the critical velocity. However, in practical terms, this does not indicate smoke back layering. Therefore, in IDA Tunnel, we set the backlayering length (BL) equal to zero under this condition. It's worth noting that this equation is derived from the study by Hu et al. (2008).

Inventory of models

Table 1 shows the central models that have been developed for the program and some of their parameters.

Table 1: Modelica Models in IDA Tunnel Library

Model name	Function	Parameters
Entry	Entry portal	Tunnel area, loss coefficient
Branch	Junction	Tunnel areas in main and branch, angle between main and branch
RailSection, RoadSection	Tunnel /shaft section	Friction factor, discretisation, tunnel area, hydraulic diameter, loss coefficients at area changes, jet fan data, equipment data, fire data
Platform	Platform	Size dimensions, number of occupants
Plenum	Plenum	Size dimensions
Passage	Passenger passage	Number of gates (1 or 2), doors/gate, stochastic and intermittent passenger flows
Leak	Leak to ambient from platform	Opening area, height above floor
Portal	Portal between platform/plenum and tunnel	Loss coefficients
Recirc	Recirculation between tunnel entries	Percentage recirculating
FanDuct	Air supply/exhaust	Tunnel areas, loss coefficients
AxialFan	Axial fan	Fan curve (spline)
Damper	Damper	Loss coefficients
VFGiven	Given flow	Air volume flow
Filter	Filter particles	Air volume flow and filter efficiency
Heater	Heater	Temperature setpoint
SaccarNoz	Saccardo nozzle	Angle, momentum exchange coefficient