Hygro-thermo-mechanical analysis of spalling in concrete walls at high
temperatures as a moving boundary problem
Michal Beneša , Radek Štefanb,∗
a

arXiv:1401.1152v3 [] 9 Jan 2015

b

Department of Mathematics,
Department of Concrete and Masonry Structures,
Faculty of Civil Engineering,
Czech Technical University in Prague,
Thákurova 7, 166 29 Prague 6, Czech Republic

Abstract
A mathematical model allowing coupled hygro-thermo-mechanical analysis of spalling in concrete
walls at high temperatures by means of the moving boundary problem is presented. A simplified
mechanical approach to account for effects of thermal stresses and pore pressure build-up on spalling
is incorporated into the model. The numerical algorithm based on finite element discretization in
space and the semi-implicit method for discretization in time is presented. The validity of the
developed model is carefully examined by a comparison between the experimentally determined
data stated in literature and the results obtained from the numerical simulation.
Keywords: concrete, high temperature, spalling, hygro-thermo-mechanical analysis, moving
boundary, thermal stress, pore pressure, finite element method

1. Introduction
Mathematical modelling seems to be an effective and powerful tool to simulate the heated
concrete behaviour (see e.g. [22]). Several models based on more or less general physical background
have been developed to simulate transport processes in heated concrete (see [27] and references
therein). All developed models build on a system of conservation of mass and energy, but differ
in the complexity of phase description of state of pore water as well as chemical reactions and
different physical mechanisms of coupled transport processes in a pore system. A descriptive
phenomenological approach was used to explore hygro-thermal processes in concrete exposed to
temperatures exceeding 100 ◦ C starting with the Bažant & Thonguthai model [5]. Here, the liquid
water and water vapour are treated as a single phase, moisture, and the evaporable water is assumed
to be formed by the capillary water only. The main advantages for usage of this approach is relative
simplicity and small number of parameters that can be obtained from experiments. However, in
such a way it is impossible to consider the effects of phase changes of water and the applicability
of the single-fluid-phase models for temperatures above the critical point of water is disputed.
These deficiencies led to the development of more detailed multi-phase description, see e.g. the
works of Gawin et al. [21] and Davie et al. [13] for specific examples. Coupled multi-phase models
reflect the multi-phase structure of concrete, interactions between phases, phase changes of fluids
and solids and non-linear couplings between thermal, hygral and mechanical processes. However,
such increase in complexity comes at the expense of a large number of model parameters, which
determination can be hardly obtained directly from experiments. Moreover, multi-phase models
are computationally expansive. Despite rapid progress in computer technologies, complex models
∗

Corresponding author. Tel.: +420 22435 4633; fax: +420 23333 5797
Email address: radek.stefan@fsv.cvut.cz (Radek Štefan)

Preprint submitted to International Journal of Heat and Mass Transfer

January 12, 2015

still exceed capabilities of recently developed numerical algorithms and computational hardware.
Therefore, in this paper we will adopt a pragmatic concept and consider the simplified model
obtained directly from the complex multi-phase description. Relative importance of thermodynamic
fluxes will be quantitatively evaluated at the level of material point and these results will allow
us to neglect less important transport phenomena without lost of capability to realistically predict
behavior of concrete at extremely high temperatures.
High-temperature exposure of concrete can lead to the risk of concrete spalling and, consequently, to the damage of the entire structure (see e.g. [25] and references therein, for the examples
see [34, 53, 54]). It is generally accepted that the spalling process in rapidly heated concrete is
caused by two main processes – increase of pore pressure end development of thermal stresses –
that may act separately or, which is more likely, in a combined way (see e.g. [38, Section 4]).
In literature, we can find many criteria to asses the risk and, in some cases, also the amount of
concrete spalling (for a brief summary of some of these criteria, see e.g. [25, 38]). In our previous
work [7], we have employed a heuristic engineering approach, originally proposed by Dwaikat and
Kodur [15, Section 3.3], in which the spalling is supposed to occur if the effective pore pressure
exceeds the temperature dependent tensile strength of concrete. In the present paper, we extend
this criterion in order to take into account not only the pore pressure (which seems to be not the
dominant mechanism of spalling, as observed by recent experimental investigations [31, 46] and
numerical simulations [48]) but also the thermal stress as a driving force of spalling.
The paper is organized as follows. In Section 2, we specify general thermodynamical and
mechanical assumptions on concrete as a porous multi-phase medium to obtain a reasonably simple
but still realistic model to predict hygro-thermal behavior of concrete at very high temperatures.
This Section is concluded by presentation of conservation of mass and thermal energy, carefully
derived in Appendix A by quantitative parameter analysis. In Sections 3 and 4, constitutive
relationships are discussed in details and suitable boundary conditions for description of transport
processes through the surfaces of the concrete wall are presented, respectively. Section 5 deals with
the problem of moving boundary for spalling simulation. Our approach is based on the combined
effect of pore pressure and thermal stresses in concrete under high temperature exposure. In
Section 6, the complex problem is formulated as a fully coupled system of highly nonlinear partial
differential equations (PDE’s) supplemented with appropriate boundary and initial conditions.
Based on the full FEM (in space) and semi-implicit (in time) discretization of the mentioned
system of PDE’s, an in-house research MATLAB code has been developed for the solution of
the system of nonlinear algebraic equations. The numerical algorithm is presented in Section 7.
Section 8 brings the complete list of material properties of moist concrete at high temperatures
used in the numerical model. Section 9 is the key part of the paper. Here we present validation
of the model by comparison of the numerical results with experiments reported in literature by
means of three examples. First two examples, performed by Kalifa et al. [36] and, more recently,
Mindeguia [44], examine the comparison of the measurements of pore pressure and temperature
distributions with the predicted numerical simulations based on the present model under specific
conditions excluding spalling phenomena. Finally, we apply the present model to investigate the
surface spalling of high strength concrete prismatic specimen under unidirectional heating by the
ISO 834 fire curve and compare the numerical results with experimental observations reported by
Mindeguia [44] and Mindeguia et al. [45, 46]. The summary of the outcomes achieved in this paper
as well as the general conclusions and recommendations for future research appear in Section 10.
2. Basic conservation equations
2.1. General assumptions
A number of assumptions have been adopted to develop the coupled hygro-thermo-mechanical
model for concrete spalling due to high temperatures exposure. Some of them of particular interest
2

are as follows:
• concrete is assumed as a multi-phase system consisting of different phases and components:
solid skeleton (composed of various chemical compounds and chemically bound water), liquid
phase (combined capillary and adsorbed water), gas phase (a mixture of dry air and water
vapour) (see e.g. [21]);
• the diffusive mass flux of water vapour is neglected and the adsorbed water diffusion is
assumed to be expressed by the liquid water relative permeability term, Krw , instead of a
separate term (see [10, 13]);
• the effects of variations of pressure of dry air is neglected. The mass of dry air in concrete is
considered to be much smaller than the mass of liquid water and vapour and, consequently,
the vapour pressure plays the crucial role (when compared to dry air) in spalling phenomena
(see [15]);
• above the critical point of water, only the vapour contribution to the mass transport is taken
into account (cf. [23]);
• the spalling of concrete in a heated wall is assumed to be caused by a combination of the
hygro-thermal stress due to the pore pressure build-up and the thermo-mechanical stress
resulting from the restrained thermal dilatation (see e.g. [38]);
• the wall is considered as fully mechanically restrained in the plane perpendicular to the wall
thickness (see [48, 58]);
• the stresses resulting from the external mechanical load as well as the self-weight of the wall
are not accounted for since their effect on the potential spalling may be considered to be
negligible compared with the effect of mechanical restraint (see [48, 58]);
• the hygro-thermal and the mechanical problems are coupled in the analysis.
2.2. Conservation laws
The mathematical model of moisture and heat transfer in concrete consists of the balance
equations governing the conservation of mass (moisture) and thermal energy (cf. equations (A.9)
and (A.23)).
The mass balance equation of moisture (liquid water and vapour):
∂
∂
∂md
(ηw ρw + ηv ρv ) +
(ηw ρw vw + ηv ρv vv ) =
;
∂t
∂x
∂t

(1)

energy conservation equation for moist concrete as the multi-phase system:
(ρcp )

∂θ
∂qc
∂θ
∂me
∂md
=−
− (ρcp v)
−
he −
hd ,
∂t
∂x
∂x
∂t
∂t

(2)

g
s
(ρcp ) = cw
p ρw ηw + cp ρg ηg + cp ρs ηs ,

(3)

where

(ρcp v) =

cw
p ρw ηw vw

+

cgp ρg ηg vg

(4)

and cgp ρg = cvp ρv + cap ρa . Governing equations (1) and (2) are carefully derived in Appendix A and
the meanings of all symbols are explained in a well arranged way in Appendix B. It should be
underlined that material coefficients of concrete and fluids that appear in governing equations (1)
and (2) depend in non-linear manner upon the primary unknowns – temperature and pore pressure,
which completely describe the state of concrete under thermal loading.
3

3. Constitutive relationships
Balance equations (1) and (2) are supplemented by an appropriate set of constitutive equations.
Moisture flux. As the constitutive equations for fluxes of fluid phases (liquid water and vapour)
the Darcy’s law is applied
KKrw ∂Pw
,
µw ∂x
KKrg ∂P
,
= −ρv
µg ∂x

ηw ρw vw = −ρw
η v ρ v vv

Pw = P − Pc ,

(5)
(6)

where Pw [Pa] is the pressure of liquid water, Pc [Pa] is the capillary pressure and P [Pa] represents
the pore pressure due to the water vapour. Further, K [m2 ] represents the intrinsic permeability,
Krw [−] and Krg [−] are the relative permeability of liquid water and relative permeability of gas,
µw [Pa s] and µg [Pa s] represent the liquid water and gas dynamic viscosity.
The equilibrium state of the capillary water with the water vapour is expressed in the form of
the Kelvin equation
 
θR
P
Pc = −ρw
ln
,
(7)
Mw
Ps
where Pc denotes the capillary pressure and the water vapour saturation pressure Ps [Pa] can be
calculated from the following formula as a function of temperature θ [K]


4042, 9
Ps (θ) = exp 23, 5771 −
.
(8)
θ − 37, 58
Water vapour is considered to behave as perfect gas, therefore, Clapeyron equation [21]
ρv =

P Mw
θR

(9)

is assumed as the state equation, where R [J mol−1 K−1 ] is the gas constant (8.3145 J mol−1 K−1 )
and Mw [kg mol−1 ] represents molar mass of water vapour.
Heat flux. For the heat flux induced by conduction the Fourier’s law is applied in the form
qc = −λc

∂θ
,
∂x

(10)

where λc [W m−1 K−1 ] represents the temperature and saturation dependent effective thermal conductivity of moist concrete.
Evaporation. In order to determine the amount of heat due to evaporation or, reversely, condensation processes, the water vapour conservation equation
∂me
∂(ηv ρv ) ∂ (ηv ρv vv )
+
=
∂t
∂x
∂t

(11)

needs to be incorporated into the energy conservation equation (2) which will be handled in Section
6.2.

4

Dehydration. Following Dal-Pont and Ehrlacher [11], Feraille-Fresnet et al. [19], the evolution of
mass source term md [kg m−3 ] related to the dehydration process is considered through the following
evolution law
∂md
1
= − (md − md,eq (θ)),
(12)
∂t
τ
where md,eq [kg m−3 ] is the mass of water released at the equilibrium according to temperature θ and
τ [s] is the characteristic time of mass loss governing the asymptotic evolution of the dehydration
process.
4. Boundary conditions
To describe coupled transport processes through the surfaces of the wall, one should prescribe
the appropriate boundary conditions across the boundary. Homogeneous Neumann conditions
∂θ
= 0,
∂x
ρw ηw vw + ρv ηv vv = 0
he ρv ηv vv − λc

are usually applied on the insulated surface of the wall. Boundary conditions of the form


∂θ
4
he ρv ηv vv − λc
nx = αc (θ − θ∞ ) + eσSB (θ4 − θ∞
),
∂x
(ρw ηw vw + ρv ηv vv )nx = βc (ρv − ρv∞ )

(13)
(14)

(15)
(16)

(nx = ±1) are of importance on the exposed side of the wall, where the terms on the right hand
side of (15) represent the heat energy dissipated by convection and radiation to the surrounding
medium and the term on the ride hand side of (16) represents a water vapour dissipated into the
surrounding medium.
5. Moving boundary for spalling simulation in concrete walls
5.1. Spalling criterion
Let us assume a concrete wall of a thickness ` exposed to fire on boundary x = `. The spalling
at position x ∈ (0, `) and time t occurs if (cf. [48, p. 613]; [51, Section 2.3])
F (fc (θ), ft (θ), σht (P, θ), σtm (θ)) > 1,

(17)

where F is a dimensionless failure parameter (failure function), fc and ft are the temperature
dependent uniaxial compressive and tensile strengths of concrete, respectively, σht and σtm are the
actual stresses in concrete caused by hygro-thermal and thermo-mechanical processes, respectively,
and P and θ are the pore pressure and temperature in concrete at position x and time t, respectively.
Note that the strengths of concrete (both compressive and tensile) are assumed to be positive
values while the strains and stresses are taken as positive in tension and negative in compression.
5.2. Hygro-thermal stress
Hygro-thermal stress in heated concrete is a tensile stress caused by the pore pressure buildup. There are several approaches to determine the hygro-thermal stress. These approaches may
differ both in the definition of the pore pressure and also in the manner in which the pore pressure
is converted into the hygro-thermal stress. The pore pressure may be assumed to be equal to
the vapour pressure (e.g. [15]), to the gas pressure, or can be calculated as the Bishop’s stress
[12, 21, 48, 51]. The hygro-thermal stress (i.e. the effective pore pressure) can be determined
5

from the pore pressure by the hollow spherical model [30], or multiplying respectively by a Biot’s
coefficient [48, 51] or by a concrete porosity [15].
Here, we adopt the approach proposed by Dwaikat and Kodur [15], in which
σht (P, θ) = P φ(θ),

(18)

where P is the pore pressure due to water vapour and φ is the temperature dependent concrete
porosity.
5.3. Thermo-mechanical stress
The thermo-mechanical stress in a heated concrete wall generally depends on the external
mechanical load applied on the wall, on its geometry (wall thickness, load eccentricity), material
properties (both hygro-thermal and mechanical), and boundary conditions (heating, mechanical
restraint).
In our approach, we follow a conservative assumption, also adopted by e.g. Msaad [48], Zeiml
et al. [58], in which the wall is supposed to be fully mechanically restrained in the plane perpendicular to the wall thickness. On the other hand, the stresses resulting from the external mechanical
load as well as the self-weight of the wall are not accounted for since their effect on the potential
spalling may be considered to be negligible when compared with the effect of mechanical restraint
(cf. [48, 58]).
In order to asses the thermal stress arising from the mechanical restraint, we have to focus on
the stress-strain conditions in the heated wall. It is widely accepted that the total strain in concrete
subjected to high temperatures may be decomposed into several parts which differ in their physical
meaning. Hence, we can write ([2, eq. (1)]; [3, eq. (4.2)]; [41, eq. (1)])
tot = θ (θ) + σ (σ, θ) + cr (σ, θ, t) + tr (σ, θ),

(19)

where tot is the total strain, θ is the free thermal strain, σ is the instantaneous stress-related
strain (which can be divided into the elastic part and the plastic part, e and p , respectively (see
e.g. [8]), cr is the creep strain, tr is the transient strain, σ is the stress, and t is the time.
As stated in e.g. [28], the creep strain, cr , is usually neglected. Moreover, the stress-dependent
strains can be assumed together (with or without σ or cr ) as the mechanical strain (see [28]) or
as the so called load induced thermal strain – LITS (see e.g. [2, 39, 56, 58]).
Here, we adopt a constitutive law proposed by Eurocode 2 [17], in which the stress is expressed
in terms of the total mechanical strain, m , that includes the transient strain implicitly, while the
creep strain is omitted [2, 9, 28], and hence, we can write
tot = θ (θ) + m (σ, θ).

(20)

Since the wall is supposed to be fully mechanically restrained, the total strain is equal to zero
and the mechanical strain can be expressed as
m = −θ (θ).

(21)

As mentioned above, the constitutive law given by Eurocode 2 [17] has the form
σ = L(m , θ).

(22)

In our case (see equation (21)), the stress can be expressed directly as a function of temperature
(the formulas provided by Eurocode 2 [17] for θ (θ) and L(m , θ) are stated in Section 8). It should
be noted that the stress in concrete determined by (22) belongs to the uniaxial conditions. For the
plane stress conditions, we get
1
σtm (θ) =
σ(θ),
(23)
1 − ν(θ)
where ν is the Poisson’s ratio (see Section 8), and σ is the uniaxial stress given by (22) (with m
determined by (21)).
6

5.4. Resulting stress conditions
Let us assume that the wall thickness, `, which is much smaller than the other dimensions of
of the wall, is parallel to the x-axis, and the wall surfaces (at positions of x = 0 and x = `) are
parallel to the plane y-z (cf. [48]). As mentioned above and as stated in [48], we can suppose
that the hygro-thermal stress acts as a tensile stress perpendicular to the wall surface (i.e. in
the x-direction) and the thermo-mechanical stress acts as a compressive plane stress in the planes
parallel to the wall surface (i.e. in the directions of y and z), see Figure 1. Hence, in terms of the
principal stresses, σ1 , σ2 , σ3 , and the normal stresses, σx , σy , σz , we may write [48]
σ1 = σx = σht (P, θ),

(24)

σ2 = σ3 = σy = σz = σtm (θ).

(25)



σtm
P
σht

θ

σtm

σht

z

σtm
y

x
Figure 1: Stress conditions in the wall

5.5. Failure function
As a failure function, we adopt the Menétrey–Willam triaxial failure criterion defined as [43,
eq. (6)]

 2



√
fc − ft2 eF
ρ 2
ρ
ξ
√ r(ϑ, eF ) + √
F =
1.5
,
(26)
+ 3
fc
fc ft eF + 1
6fc
3fc
where eF is the dimensionless eccentricity (0.5 < eF ≤ 1.0), which influences the shape of the
failure function (see e.g. [43, Fig. 3]) and can be expressed as a function of fc , ft and fb , with fb
being the biaxial compressive strength of concrete (see [33, eq. (21.28)]). Further, function r(ϑ, eF )
is given by Menétrey and Willam [43, eq. (1)] in the form
r(ϑ, eF ) =

4(1 − e2F ) cos2 ϑ + (2eF − 1)2
q
,
2(1 − e2F ) cos ϑ + (2eF − 1) 4(1 − e2F ) cos2 ϑ + 5e2F − 4eF
7

(27)

and ξ, ρ and ϑ are the Haigh–Westergaard coordinates (see e.g. [33, 43]).
The Menétrey–Willam failure function for the plane stress conditions (i.e. for σ3 = 0) is shown
in Figure 2.
In our case, in which σ1 = σht and σ2 = σ3 = σtm , we can write
1
ξ = √ (σht + 2σtm ),
3
r
2
ρ=
(σht − σtm )2 ,
3
σht − σtm
.
cos (3ϑ) = p
(σht − σtm )2

(28)
(29)
(30)

Since σtm ≤ 0 < σht , and 0.5 < eF ≤ 1.0, it is obvious that ϑ = 0, and r = 1/eF . Substituting
this in (26) and assuming the strengths of concrete as temperature dependent, we get the resulting
failure function used in our model (cf. [48, p. 613])

F =

σht (P, θ) − σtm (θ)
fc (θ)

2

fc (θ)2 − ft (θ)2
+
fc (θ)2 ft (θ)



2eF − 1
σht (P, θ) +
σtm (θ) ,
eF + 1

(31)

which is illustrated in Figure 3.
σ2

−fb −fc

σtm (θ)

ft
ft

F <1

F =1

ft (θ)

σ1

F >1
−fc
−fb

σht (P, θ)

−fb (θ)

Figure 2: Menétrey–Willam failure function [43]

Figure 3: Failure function used in our model

5.6. Evolution law for moving boundary due to spalling
In our regularized approach, the instantaneous spalling of concrete is approximated by rapid
continuous process in such a way that sheet of concrete is continuously removed from the wall and
the receding outer surface forms the moving boundary. In general, in one dimension, the external
loading to the concrete wall according to the boundary conditions (15) and (16) is prescribed on
the unknown ablating boundary x = `(t). For the spalling of the wall originally occupying the
space 0 ≤ x ≤ `0 , we propose the governing equation of the form
d`
`
= − [max(F ) − 1]+ ,
dt
γ

`(0) = `0 ,

(32)

where the position of the moving boundary ` has to be determined as a part of the solution.
Further, x+ = max(0, x), γ [s] represents the characteristic time of spalling process and max(F ) is
the maximal value of failure function (31) achieved within the wall (0, `) at actual time t.

8

6. Formulation of the problem
6.1. Modification of mass balance equation
Let us start with modification of the equation (1). Denote by m = ηw ρw + ηv ρv the total mass
of moisture (including liquid water and vapour). Incorporating the constitutive relations (5) and
(6) into the mass balance equation (1) yields



KKrg ∂P
KKrw
KKrw ∂Pc
∂md
∂m
∂
ρw
+ ρv
+
=
− ρw
.
(33)
∂t
∂x
µw
µg
∂x
µw ∂x
∂t
Here Pc = Pc (P, θ) via equation (A.14). After additional modification, the equation (33) can be
written in a general form


∂m ∂md
∂P
∂
∂θ
KmP
,
(34)
−
=
+ Kmθ
∂t
∂t
∂x
∂x
∂x
where
KmP
Kmθ



KKrg
∂Pc
KKrw
= ρw
1−
+ ρv
,
µw
∂P
µg
KKrw ∂Pc
.
= −ρw
µw ∂θ

(35)
(36)

6.2. Modification of energy conservation equation
Incorporating water vapour conservation equation into the term corresponding to the latent
heat of evaporation leads to the following modified energy balance equation


 ∂θ
∂θ
∂ (ηv ρv )
∂md
∂
∂θ
∂ (ηv ρv vv )
g
(ρcp ) + he
+ hd
=
λc
− he
− cw
. (37)
p ρw ηw vw + cp ρg ηg vg
∂t
∂t
∂t
∂x
∂x
∂x
∂x
Simple calculation yields
he

∂ (ηv ρv vv )
∂x

∂ (he ηv ρv vv )
∂he
− (ηv ρv vv )
∂x

 ∂x
KKrg ∂P
KKrg ∂he ∂P ∂θ
∂
= −
he ρv
+ ρv
.
∂x
µg ∂x
µg ∂θ ∂x ∂x
=

Incorporating the equation (38) into (37) reads

 

∂P
∂θ
∂md
∂
∂P
∂θ
∂P
∂θ ∂θ
MθP
+ Mθθ
+ hd
=
KθP
+ Kθθ
+ CθP
+ Cθθ
,
∂t
∂t
∂t
∂x
∂x
∂x
∂x
∂x ∂x

(38)

(39)

where
MθP

Mθθ =
KθP

=

Kθθ =
CθP

∂ (ηv ρv )
,
∂P
∂ (ηv ρv )
(ρcp ) + he
,
∂θ
KKrg
he ρv
,
µg
λc ,

 


KKrg
KKrw
∂Pc
∂he
w
v
a
cp ρw
1−
+
cp −
ρv + cp ρa
,
µw
∂P
∂θ
µg
KKrw ∂Pc
−cw
.
p ρw
µw ∂θ

= he

=

Cθθ =

9

(40)
(41)
(42)
(43)
(44)
(45)

6.3. Resulting model
The full mathematical model consists of the balance equations for moisture and energy, state
equation of pore water (moisture), governing equation for dehydration, evolution equation for
moving boundary due to spalling, the set of appropriate boundary and initial conditions specifying
the fields of pore pressure, temperature, mass of dehydrated water and initial thickness of the wall.
Moisture conservation equation:
∂
∂m ∂md
−
=
∂t
∂t
∂x



∂P
∂θ
KmP
+ Kmθ
∂x
∂x


;

(46)

energy conservation equation:
MθP

∂P
∂θ
∂md
∂
+ Mθθ
+ hd
=
∂t
∂t
∂t
∂x


KθP

∂P
∂θ
+ Kθθ
∂x
∂x





∂P
∂θ ∂θ
+ CθP
+ Cθθ
;
∂x
∂x ∂x

(47)

state equation of moisture:
m = ηw ρw + ηv ρv ;

(48)

∂md
1
= − (md − md,eq (θ));
∂t
τ

(49)

governing equation of dehydration:

evolution law for moving boundary due to spalling:
`
d`
= − [max(F ) − 1]+ ;
dt
γ
boundary conditions (at x = 0 and x = `):


∂P
∂θ
− KmP
+ Kmθ
nx = βc (ρv − ρv∞ ),
∂x
∂x


∂P
∂θ
4
− KθP
+ Kθθ
nx = αc (θ − θ∞ ) + eσSB (θ4 − θ∞
)
∂x
∂x

(50)

(51)
(52)

and initial conditions (at t = 0):
P = P0 ,

(53)

θ = θ0 ,

(54)

` = `0 ,

(55)

md = 0.

(56)

The unknowns in the model are the moisture content m, pore pressure P , temperature θ, mass
of dehydrated water md and the actual thickness of the wall `. Transport coefficients KmP , Kmθ ,
MθP , Mθθ , KθP , Kθθ , CθP , Cθθ , defined by (35) and (36) and (40)–(45), depend in non-linear
manner upon temperature θ and pressure P . Material non-linearities are described in detail in
Section 8.

10

7. FEM formulation and solution strategy
Let 0 = t0 < t1 < · · · < tN = T be an equidistant partitioning of time interval [0, T ] with
step ∆t. Set a fixed integer n such that 0 ≤ n < N . In what follows we abbreviate f (x, tn ) by
f n for any function f . The time discretization of the continuous model is accomplished through a
semi-implicit difference scheme


n+1
n+1
mn+1 − mn
∂
n ∂P
n ∂θ
KmP
=
+ Kmθ
∆t
∂x
∂x
∂x
mn+1
− mnd
d
,
 ∆t

n+1
n+1
∂
n ∂P
n ∂θ
KθP
=
+ Kθθ
∂x
∂x
∂x
 n

n
n
∂θ
∂θ
∂P
n
n
+ Cθθ
+ CθP
∂x
∂x
∂x
+

n
MθP

n+1 − θ n
P n+1 − P n
n θ
+ Mθθ
∆t
∆t

(57)

mn+1
− mnd
d
,
(58)
∆t
n+1
= ηw
ρw (θn+1 ) + ηvn+1 ρv (θn+1 , P n+1 ), (59)
−hnd

mn+1
mn+1
− mnd
d
∆t
n+1
`
− `n
∆t

n+1


n+1
n ∂P
n ∂θ
−nx KmP
+ Kmθ
∂x
∂x


n+1
n+1
n ∂θ
n ∂P
−nx KθP
+ Kθθ
∂x
∂x

x=0, x=`n+1

x=0, x=`n+1

1
= − (mnd − md,eq (θn )),
τ
`n+1
= −
[max(F n+1 ) − 1]+ ,
γ

(60)
(61)


n+1 n+1
= βc ρv (P n+1 , θn+1 ) − ρv (P∞
, θ∞ ) , (62)
n+1
= αc (θn+1 − θ∞
)


n+1 4
+eσSB (θn+1 )4 − (θ∞
) ,

(63)

n+1 = η (θ n+1 , P n+1 ) and η n+1 = η (θ n+1 , P n+1 ). Further, n = +1 for x = `n+1 and
where ηw
w
v
x
v
nx = −1 for x = 0.
Applying the Galerkin procedure to the mass and energy conservation equations leads to the
system of non-linear algebraic equations


1 n n+1
M x
− xn + Kn xn+1 + f n+1 (xn+1 ) = 0,
(64)
∆t
T
where xn+1 = mn+1 , θ n+1 , P n+1 , mn+1 (x) = N(x)mn+1 , θn+1 (x) = N(x)θ n+1 and P n+1 (x) =
N(x)P n+1 , stores the unknown nodal values of moisture, temperature and pore pressure at time
tn+1 , respectively. The constant matrices in (64) exhibit a block structure
 n



Mmm
0
0
0 Knmθ KnmP
Mnθθ MnθP  ,
Mn =  0
Kn =  0 Knθθ KnθP 
(65)
0
0
0
0
0
0
and the non-linear term reads as

n+1 )
f n+1
m (x
f n+1 (xn+1 ) =  f n+1
(xn+1 )  .
θ
n+1 )
f n+1
P (x


11

(66)

The non-linear system (64) is solved iteratively using the Newton’s method (see [52]).
The individual matrices M•• can be expressed as
Z

Mnmm

`n+1

N(x)T N(x) dx,

=

(67)

0
`n+1




∂ (ηv ρv )
(ρcp ) + he
N(x)
=
N(x) dx,
∂θ
0


Z `n+1
∂ (ηv ρv )
T
=
N(x) he
N(x) dx,
∂P
0
Z

Mnθθ
MnθP

T

(68)
(69)

whereas the blocks K•• attain the form
Knmθ
KnmP

`n+1

T 

KKrw ∂Pc ∂
∂
=
N(x)
−ρw
N(x) dx,
∂x
µw ∂θ ∂x
0



T 

Z `n+1 
KKrg ∂
∂
∂Pc
KKrw
1−
+ ρv
=
N(x)
ρw
N(x) dx,
∂x
µw
∂P
µg
∂x
0
Z



(70)

(71)
Knθθ
KnθP

`n+1

T 

∂
∂
N(x)
λc N(x) dx,
=
∂x
∂x
0
T 

Z `n+1 
KKrg ∂
∂
N(x)
he ρv
N(x) dx.
=
∂x
µg ∂x
0
Z



(72)
(73)

The non-linear terms f • are provided by
n+1
f n+1
) =
m (x

`n+1

Z

−

0

− mnd
mn+1
d
N(x)T dx
∆t
#x=`n+1

"
+
(xn+1 ) =
f n+1
θ

βc ρv (P
`n+1

Z

mn+1
hnd d

0

−

Z
0

n+1

`n+1

,θ

−
∆t

n+1

mnd

)−

n+1 n+1
ρv (P∞
, θ∞ )



T

,

N(x)

x=0

N(x)T dx


 n
n
n
∂θ
n ∂P
n ∂θ
CθP
+ Cθθ
N(x)T dx
∂x
∂x
∂x
#x=`n+1

"
+

(74)

n+1
n+1
αc (θn+1 − θ∞
) + eσSB (θn+1 )4 − (θ∞
)

n+1
n+1 n+1
f n+1
) = mn+1 − ηw
(θ
, P n+1 )ρw (θ n+1 )
P (x

−ηvn+1 (θ n+1 , P n+1 )ρv (θ n+1 , P n+1 ).

12


4

N(x)T

,

(75)

x=0

(76)

Numerical algorithm. The described numerical algorithm becomes:
Step 1. Set T , ∆t, N = T /∆t
Step 2. Set initial values P0 , θ0 , `0 , md (0) = 0
Step 3.
For n = 0, . . . , N − 1
tn = n∆t
update md,eq (θn ), F n+1
solve (61) to get `n+1
solve (60) to get mn+1
d
n+1 , θ n+1
update P∞
∞

update Mn , Kn , f n+1
solve (64) by Newton’s iteration procedure to get xn+1
end n
8. Material data for concrete at high temperatures
In our approach, concrete is assumed to be a homogeneous multi-phase system. Hence, most
of its material properties (such as free thermal strain, thermal conductivity, etc.) are treated as so
called effective or smeared characteristics. It means that for the solid skeleton, we do not distinguish
between the aggregates and the cement paste. Such macroscopic approach has been widely accepted
by the scientific community for investigation of hygro-thermo-mechanical and spalling behaviour
of heated concrete (see e.g. the works by Dwaikat and Kodur [15], Gawin et al. [25], Witek et al.
[57]). It is clear that by turning to a mesoscale level, one can obtain more realistic results, especially
when investigating the spalling phenomenon (see e.g. the recent works by Le [40], Zhao et al. [59]).
However, such approach is out of scope of the present paper.
Based on the literature review, the material properties of concrete and its components subjected
to high temperatures are assumed as follows.
The free thermal strain of concrete, θ [−], is defined in the temperature range of 293.15 K to
1473.15 K as [17, Section 3.3.1]
for siliceous aggregates concrete:

−1.8 × 10−4 + 9 × 10−6 θ + 2.3 × 10−11 θ3 for θ ≤ 973.15 K,
θ (θ) =
14 × 10−3
for θ > 973.15 K,
for calcareous aggregates concrete:

−1.2 × 10−4 + 6 × 10−6 θ + 1.4 × 10−11 θ3 for θ ≤ 1078.15 K,
θ (θ) =
12 × 10−3
for θ > 1078.15 K.

(77)

(78)

The constitutive law for concrete in compression, σ = L(m , θ), is adopted from Eurocode 2 in
the form [17, Section 3.2.2.1]

3 f (θ)


" m c

3 # for cu1 (θ) < m ≤ 0,
 −
m
σ(m , θ) =
(79)
c1 (θ) 2 +

c1 (θ)



0
for m ≤ cu1 (θ),
13

m [−]
-0.05 -0.04 -0.03 -0.02 -0.01

c1 (θ)

cu1 (θ)

1273.15 K
1073.15 K

m

00
-0.2
-0.4

873.15 K
-0.6

673.15 K
473.15 K
293.15 K

−fc (θ)

σ(θ)/fc,ref [−]

σ

-0.8
-1

Figure 4: Constitutive law for concrete in compression Figure 5: An example of the constitutive law for congiven by Eurocode 2 [17, Figure 3.1]
crete in compression at high temperatures [17]

where m is the total mechanical strain, fc is the compressive strength of concrete, c1 is the strain
corresponding to fc , and cu1 is the ultimate strain, see Figure 4 and Figure 5.
The temperature dependent parameters of the constitutive law (79), fc (θ), c1 (θ), and cu1 (θ),
for the normal strength concrete with the siliceous and calcareous aggregates, NSC-S and NSC-C,
respectively, as well as for three classes of high strength concrete (HSC-1, HSC-2, and HSC-3, see
Eurocode 2 [17, Section 6.1]) are stated in Eurocode 2 [17, Table 3.1, Table 6.1N] in the form of
tabulated data that are illustrated in Figures 6 and 7.
1

0
NSC-S
NSC-C
HSC-1
HSC-2
HSC-3

0.6

-0.01
c1 , cu1 [−]

fc /fc,ref [−]

0.8

0.4
0.2
0
273.15

c1 (θ)
cu1 (θ)

-0.02
-0.03
-0.04

673.15

1073.15
θ [K]

-0.05
273.15

1473.15

673.15

1073.15
θ [K]

1473.15

Figure 6: Reduction of the compressive strength of NSC Figure 7: Temperature evolution of the mechanical
and HSC given by Eurocode 2 [17]
strains of concrete according to Eurocode 2 [17]

The constitutive law of concrete in tension need not be defined in our approach. It is sufficient
to describe the tensile strength, ft [Pa], that can be assumed as [15, eq. (34)]

1
for θ ≤ 373.15 K,



(873.15 − θ)/500
for 373.15 K < θ ≤ 823.15 K,
(80)
ft (θ) = ft,ref ×
(1473.15 − θ)/6500 for 823.15 K < θ ≤ 1473.15 K,



0
for θ > 1473.15 K,
where ft,ref [Pa] is the reference tensile strength of concrete at the room temperature.
If the reference tensile strength, ft,ref [Pa], is not known, it can be estimated from the reference compressive strength of concrete at the room temperature, fc,ref [MPa]. In Eurocode 2 [16,
Table 3.1], the tensile strength of high strength concrete is defined as


fc,ref
ft,ref = 2.12 ln 1 +
[MPa].
(81)
10
14

The Young’s modulus of concrete, Ec [Pa], can be easily derived from (79) in the form [see e.g.
18, eq. (3.14)]
3fc (θ)
Ec (θ) =
.
(82)
2c1 (θ)
The Poisson’s ratio of concrete, ν [−], is usually taken as temperature independent (see e.g. [1,
p. 46]) since the data about its temperature dependency are ambiguous [49, Section 2.2.1.2]. Here,
we follow the assumption (which leads to conservative results) that the Poisson’s ratio increases
with increasing temperature. Based on the data experimentally determined by Mindeguia [44,
pp. 132–133] for high strength concrete (fc,ref ≈ 60 MPa), we assume that

for θ ≤ 293.15 K,
 0.2
0.2 + 0.5(θ − 293.15)/580 for 293.15 K < θ ≤ 873.15 K,
(83)
ν(θ) =

0.7
for θ > 873.15 K.
The eccentricity, eF [−], used in the failure function of concrete (Section 5.5) generally depends
on the type of concrete and probably also on temperature. Here, we assume constant (temperature
independent) value, eF = 0.505, derived from the data measured by He and Song [29, Table 2] for
high strength concrete (fc,ref ≈ 60 MPa) under biaxial compression tests at high temperatures, see
Figure 8, where fc (θ) is taken from He and Song [29, Table 2], ft (θ) is assumed according to (80),
with ft,ref = 4 MPa, and F (θ) is determined by (26), with σ3 = 0.
The porosity of concrete, φ [−], may be expressed as [21, eq. (41)]
φ(θ) = φref + Aφ (θ − θref ),

(84)

where φref [−] is the reference porosity (at the reference temperature, θref [K]), and Aφ [K−1 ] is a
concrete-type-depend constant (see [21, pp. 46–47]).
The thermal conductivity of concrete, λc [W m−1 K−1 ], is given by [21, eqs. (46–47)]


4 φ(θ) ρw (θ) Sw (P, θ)
λc (P, θ) = λd (θ) 1 +
,
(85)
(1 − φ(θ))ρs
with
λd (θ) = λd,ref [1 + Aλ (θ − θref )] ,

(86)

where λd,ref [W m−1 K−1 ] is the reference thermal conductivity of a dry concrete (at the reference
temperature, θref [K]), and Aλ [K−1 ] is an experimentally determined coefficient.
For the intrinsic permeability of concrete, K [m2 ], we adopt the Bary function recommended
by Davie et al. [14, eq. (13)]
K(P, θ) = Kref × 104D(P,θ) ,
(87)
where Kref [m2 ] is the reference permeability at the room temperature, and D [−] is the multiplicative damage parameter that can be defined as ([13, eq. (70)]; [24, eq. (54)])
D(P, θ) = Dm (P, θ) + Dθ (θ) − Dm (P, θ)Dθ (θ),

(88)

where Dm [−] is the mechanical damage parameter, which is assumed to be equal to the failure
parameter (31) in our approach, i.e.
Dm (P, θ) = F (P, θ),
and Dθ [−] is the thermal damage parameter that we define as (cf. [24, eq. (52)])


1
1 Ec (θ)
fc (θ)
ft (θ)
Dθ (θ) = 1 − (Dθ,Ec + Dθ,f c + Dθ,f t ) = 1 −
+
+
.
3
3 Ec,ref
fc,ref
ft,ref
15

(89)

(90)

It should be noted that usually, the thermal damage parameter is based only on the degradation
of the Young’s modulus (see [13, eqs. (29–30)]; [24, eq. (52)]). This is however not possible to
assume in our approach since we employ the constitutive model proposed by Eurocode 2 [17], see
(79), where the stress is expressed as a function of the total mechanical strain, m , instead of
the instantaneous stress-related strain, σ , and hence the resulting softening of the material with
increasing temperature is higher than in reality (see also [2, 28]). In order to eliminate this effect,
we propose to include not only the temperature dependent degradation of Ec but also the reduction
of the other mechanical properties (fc , ft ), which leads to good agreement of the resulting thermal
damage parameter with the data stated in literature, see Figure 9.
1

σ1 (θ) [MPa]
-90
-60
-30
Test
Model

0 1010
0

-60

Dθ [−]

-30

σ2 (θ) [MPa]

873.15 K

0.8

Dθ,Ec
Dθ,f c
Dθ,f t
Dθ
Dθ,[a]
Dθ,[b]

0.6
0.4
0.2
0
273.15

293.15 K

673.15

-90

1073.15
θ [K]

1473.15

Figure 8: Comparison of failure function (26) used in Figure 9: Thermal damage assumed in our model (Dθ )
our model with the experimental data measured by He and its comparison with Davie et al. [13, eq. (30)],
and Song [29, Table 2]
Dθ,[a] , and Gawin et al. [24, eq. (52)], Dθ,[b]

The gas relative permeability and the liquid water relative permeability, Krg [−] and Krw [−],
respectively, can be expressed by the formulas experimentally determined by Chung and Consolazio
[10, eqs. (7–9)]
Krg (P, θ) = 10Sw (P,θ) ψ(θ) − 10ψ(θ) Sw (P, θ),
(91)
Krw (P, θ) = 10(1−Sw (P,θ)) ψ(θ) − 10ψ(θ) (1 − Sw (P, θ)),

(92)

ψ(θ) = 0.05 − 22.5 φ(θ).

(93)

with
The gas dynamic viscosity, µg [Pa s], is approximated as [21, eq. (50)]

µg (P, θ) = µgv (θ) + [µga (θ) − µgv (θ)]

Pa
P + Pa

0.608
,

(94)

with
µgv (θ) = µgv,ref + αv (θ − θref )

(95)

µga (θ) = µga,ref + αa (θ − θref ) + βa (θ − θref )2 ,

(96)

µw (θ) = 0.6612(θ − 229)−1.532 .

(97)

and
where µgv,ref = 8.85 × 10−6 Pa s, αv = 3.53 × 10−8 Pa s K−1 , µga,ref = 17.17 × 10−6 Pa s, αa =
4.73 × 10−8 Pa s K−1 , βa = 2.22 × 10−11 Pa s K−2 , and θref = 273.15 K.
The liquid water dynamic viscosity, µw [Pa s], can be calculated from [21, eq. (51)]

16

The volume fraction of liquid (all free) water, ηw [−], is governed by the sorption isotherm
function proposed in [4, 5, 6]. After some modification related to the saturated region, P/Ps ≥ 1.0
(see [13, eq. (73)]), we can write
if θ ≤ θcr :



1/m(θ)
φ(θref )ρw (θref ) P
P
c



for
≤ 0.96,


ρ
(θ)
c
P
(θ)
P

w
s
s (θ)



 X
3
i
P
P
ηw (P, θ) =
ξ
(θ)
−
0.96
< 1.00,
for 0.96 <
i

Ps (θ)
Ps (θ)


i=0



P

 φ(θ)
for
≥ 1.00,
Ps (θ)
if θ > θcr :
ηw = 0,

with
m(θ) = 1.04 −

22.34(θref

(θ − 263.15)2
,
− 263.15)2 + (θ − 263.15)2

(98)

(99)

where c [kg m−3 ] is the mass of cement per unit volume of concrete, θref [K] is the reference
temperature (θref = 298.15 K), the term [φ(θref )ρw (θref )] expresses the saturation water content
at 298.15 K, and ξi (θ) are temperature dependent coefficients that ensure a smooth transition
between the unsaturated and saturated region, thus the sorption isotherm function and its first
derivative with respect to (P/Ps (θ)) are continuous (cf. [12, 13, 42, 55]). Therefore, we can write
ξ0 (θ) = ηw,0.96 (θ),
ξ1 (θ) =
ξ2 (θ) =
ξ3 (θ) =

(100)

0
ηw,0.96
(θ),

(101)
0
2ηw,0.96
(θ)

0
ηw,1.00
(θ)

+
3(ηw,1.00 (θ) − ηw,0.96 (θ))
−
,
2
0.04
0.04
0
0
(θ) + ηw,1.00
(θ)
2(ηw,0.96 (θ) − ηw,1.00 (θ)) ηw,0.96
+
,
3
2
0.04
0.04

where
c
ηw,0.96 (θ) =
ρw (θ)
ηw,1.00 (θ) = φ(θ),



1/m(θ)
φ(θref )ρw (θref )
0.96
,
c

(102)
(103)

(104)

and

dηw,0.96 (θ)
dηw,1.00 (θ)
0
and ηw,1.00
(θ) =
.
dθ
dθ
The degree of saturation with liquid water, Sw [−], is defined as
0
ηw,0.96
(θ) =

Sw (P, θ) =

ηw (P, θ)
.
φ(θ)

(105)

(106)

The mass of dehydrated water, md [kg m−3 ], is governed by equation (12), with ([19, eq. (8)];
[11, eq. (C-22)]; [1, eq. (1.9)])



7.5 378.15
θ − 378.15
md,eq (θ) =
meq
1 − exp −
H(θ − 378.15)
100
200



2 378.15
θ − 673.15
+
m
1 − exp −
H(θ − 673.15)
100 eq
10



1.5 378.15
θ − 813.15
+
m
1 − exp −
H(θ − 813.15),
(107)
100 eq
5
17

378.15 [kg m−3 ] is the equilibrium mass at 378.15 K and H is the Heaviside function.
where meq
In [11, p. 146], the following values of material parameters that belong to equations (12) and
(107) are stated: τ = 10800 s, m378.15
= 210 kg m−3 , which are also adopted in our simulations,
eq
since these parameters are usually not measured within the fire tests.
The density of solid skeleton, ρs [kg m−3 ], is, together with the concrete porosity end the mass
of water released into the pores by dehydration, governed by the solid mass conservation equation
(A.8) (cf. [21, pp. 47–48]). In our approach, the density of solid skeleton is assumed to be a
constant value. As obvious from numerical experiments, this simplification has negligible effect on
the obtained results.
The density of liquid water, ρw [kg m−3 ], can be expressed by the experimentally determined
formula, originally proposed by Furbish [20, eq. (4.79)] and simplified by Gawin et al. [23, eq. (35)]
(neglecting the water pressere dependence of ρw ), of the form

ρw (θ) =

5
X
i=0

5
 X
ai (θ − 273.15)i −1 × 107 +
bi (θ − 273.15)i
i=0

for θ ≤ θcr ,

(108)

where a0 = 4.8863 × 10−7 , a1 = −1.6528 × 10−9 , a2 = 1.8621 × 10−12 , a3 = 2.4266 × 10−13 , a4 =
−1.5996×10−15 , a5 = 3.3703×10−18 , b0 = 1.0213×103 , b1 = −7.7377×10−1 , b2 = 8.7696×10−3 , b3 =
−9.2118 × 10−5 , b4 = 3.3534 × 10−7 , b5 = −4.4034 × 10−10 .
The specific enthalpy of evaporation, he [J kg−1 ], can be expressed by the Watson formula ([21,
eq. (49)]; [26, eq. (28)]; [12, eq. (AI.27)])

2.672 × 105 (θcr − θ)0.38 for θ ≤ θcr ,
(109)
he (θ) =
0
for θ > θcr .
The specific enthalpy of dehydration, hd [J kg−1 ], is considered to be a constant value [12,
eq. (AI.26)]
hd = 2400 × 103 J kg−1 .
(110)

The specific heat capacity of dry air, cap [J kg−1 K−1 ], is determined by Davie et al. [12,
eq. (AI.28)]
cap (θ) = a θ3 + b θ2 + c θ + d,
(111)

where a = −9.84936701814735×10−8 , b = 3.56436257769861×10−4 , c = −1.21617923987757×10−1 ,
and d = 1.01250255216324 × 103 .
−1 K−1 ], is given by [12, eq. (AI.29)]
The specific heat capacity of liquid water, cw
p [J kg

b
aθ
w
cp (θ) = (2.4768 θ + 3368.2) +
for θ ≤ θcr ,
(112)
513.15
where a = 1.08542631988638 and b = 31.4447657616636.
The specific heat capacity of solid skeleton, csp [J kg−1 K−1 ], can be calculated from [12, eq. (AI.30)]


θ − 273.15
θ − 273.15 2
s
cp (θ) = 900 + 80
−4
.
(113)
120
120
For the specific heat capacity of water vapour, cvp [J kg−1 K−1 ], the following equation is proposed
by Davie et al. [12, eq. (AI.31)]


b
aθ

(7.1399 θ − 443) +
for θ ≤ θcr ,
cvp (θ) =
(114)
513.15

45821.01
for θ > θcr ,
where a = 1.13771502228162 and b = 29.4435287521143.
It should be noted that above the critical temperature (if θ > θcr ), there is no liquid water in
concrete (see equation (98)) and hence, the liquid water properties (ρw , cw
p ) need not be defined.
18

9. Numerical results and experiments
In order to validate the model presented in this paper, the results obtained by numerical simulations are compared with the experimentally determined data. Three experiments reported in
literature are closely investigated:
• the test of the hygro-thermal behaviour of high strength concrete (fc,ref ≈ 90 MPa) prismatic specimen 300 × 300 × 120 mm3 under unidirectional heating by the radiant heater of
a temperature of 600 ◦ C reported by Kalifa et al. [36] and related publications Kalifa et al.
[35, 37];
• the test of the hygro-thermal behaviour of high strength concrete (fc,ref ≈ 60 MPa) prismatic specimen 300 × 300 × 120 mm3 under unidirectional heating by the radiant heater of a
temperature of 600 ◦ C reported by Mindeguia [44] and related publications Mindeguia et al.
[45, 47];
• the test of the spalling behaviour of high strength concrete (fc,ref ≈ 60 MPa) prismatic
specimen 700 × 600 × 150 mm3 under unidirectional heating by the ISO 834 fire reported by
Mindeguia [44] and related publications Mindeguia et al. [45, 46].
Within the first two of the above tests, the pore pressure and the temperature propagation
through the specimen as well as its mass loss were recorded. This type of tests is therefore denoted
as the ”PTM test” by Mindeguia [44], which is also adopted in this paper. It should be noted that
no spalling was observed during the both PTM experiments and hence, these tests are employed
primarily for the validation of the hygro-thermal behaviour simulation while the spalling prediction
can be validated only on the qualitative level (whether it occurs or not). For that reason, the third
test is also analysed, since it enables to validate our model also with respect to the quantitative
prediction of spalling (i.e. the prediction of the amount of the potential spalling).
It should be noted that the usage of the model can provide only an approximate picture of the
experiments described above. The inaccuracies may arise mainly from the fact that:
• within the mechanical part of the model, a simplified constitutive law of concrete given by
Eurocode 2 [17] is utilized. In this model, the transient strain is included implicitly and the
creep strain is neglected (see Section 5.3);
• in the simulation, the specimen is assumed to be fully mechanically restrained in the plane
perpendicular to its heated surface. This assumption is however disputable in this cases. It
is not obvious, whether the ceramic blocks placed on the lateral sides of the small specimens
within the PTM tests (see Fig. 10) or the unheated parts of the large specimens within the
spalling test (see [44, Figure 170]) provide a sufficient level of restraint as assumed in the
model (probably not);
• the model is not able to capture the ”size effect” influencing the spalling behaviour of concrete
specimens - i.e. that for small specimens, the spalling is less likely to occur in comparison
with large specimens, as observed by recent experimental works [32, 45].
9.1. Simulation of the PTM test 1
The Kalifa’s experiments [35, 36, 37] have been accepted as a benchmark problem by many
authors focused on modelling of concrete subjected to high temperatures (cf. [1, Section III-3];
[13, Section 3.2]; [15, Section 5]; [22, Section 6.1]; [23, Section 6.2]; [39, Sections 4.3, 5.4]; [50,
Section 6.4]; [57, Section 3.1]).
In [37], we can find the high temperature thermal and hygral properties measured by Kalifa
et al. [37] for 4 types of concrete named as M30, M75C, M75SC and M100C. The M30 and M100
19

concretes were used for the subsequent PTM tests reported by Kalifa et al. [36]. Within this
experiment, the specimens were subjected to an unidirectional heating up to 450 ◦ C, 600 ◦ C, and
800 ◦ C using a radiant heater. Moreover, similar experiment – in this case for the specimens made
of C110 concrete (practically the same as the M100 concrete, only with a little bit higher amount
of cement and superplasticizer [35]) with polypropylene fibres of the content of (0–3) kg m−3 , were
made by Kalifa et al. [35] in order to assess the effect of fibres on the high-temperature behaviour
of high strength concrete.
For the validation of our model, the results measured by Kalifa et al. [36] for the M100 concrete
heated up to 600 ◦ C are employed. This case is chosen for the following reasons: our model is
primarily designed to study the behaviour of high strength concrete (such as M100); material
properties of the M100 concrete are described in detail in [37] (more comprehensive description
than for the C110 concrete in [35]); and the results obtained for the other heating conditions (up
to 450 ◦ C or 800 ◦ C) are not sufficiently reported in [36].
9.1.1. Experiment description
The tests performed by Kalifa et al. are described in detail in [36]. Therefore, only a brief
description of the experiment is stated below.
As mentioned above, two types of concrete (M30 and M100) and three different heating conditions (450 ◦ C, 600 ◦ C, and 800 ◦ C) were used by Kalifa et al. [36]. Hereafter, only the case of
M100 – 600 ◦ C is assumed. Material (thermal and hygral) properties of the M100 concrete have
been reported in [37] and are discussed in Section 9.1.2.
As shown in Figure 10, the 300 × 300 120 mm3 test specimen was subjected to an unidirectional
heating (in the 120-mm direction) by a radiant heater of the temperature of 600 ◦ C placed near its
surface for a heating period of 6 hours.
Radiant heater

Specimen
300 × 300 × 120 mm3

Ceramic
blocks
Pressure
Temperature

2 20 40
10 30 50
[mm]

Balance

Figure 10: Scheme of the test set-up (according to Kalifa et al. [36, Figure 2])

Within the test, the temperature, the pore pressure, and the mass loss of the specimen were
recorded. The temperature together with the pore pressure were measured using the combined
pressure-temperature gauges (see [35, Figure 3]) placed at the distance of (10, 20, 30, 40, and 50)
mm from the heated surface. The heated surface temperature was measured by a thermocouple
placed at the depth of 2 mm. Moreover, as obvious from the graphs listed in [36], the temperature
on the unexposed side was also recorded (probably by an additional thermocouple or by a surface
thermometer). The mass loss was monitored by a balance on which the specimen was placed during
the test (see Figure 10).

20

9.1.2. Modelling
In this section, the material properties of concrete, the initial and boundary conditions as
well as the discretization employed for the numerical modelling are discussed. The thermal and
hygral properties of concrete are assumed according to the relationships given in Section 8, with
the parameters determined from the data measured by Kalifa et al. [37]. All the parameters are
summarized in Table 1.
Table 1: Material properties and parameters used in our simulation of the PTM test 1

Parameter

Value

Unit

Type of concrete

HSC2-C

−

fc,ref
ft,ref
c
θref
φref
Aφ

91.8
4.9
414.8
293.15
0.0897
2.4457 × 10−5

MPa
MPa
kg m−3
K
−
K−1

ρs

2660

kg m−3

θref
λd,ref
Aλ
Kref

293.15
1.9759
−6.4215 × 10−4
1.3 × 10−20

K
W m−1 K−1
K−1
m2

Reference
Eurocode 2 [17];
Kalifa et al. [37]
Kalifa et al. [37, Table 1]
Eurocode 2 [16], see (81)
Kalifa et al. [37, Table 1]
determined from
Kalifa et al. [37, Table 2],
see Figure 11
determined from
Kalifa et al. [37, Table 4],
see Figure 13
determined from
Kalifa et al. [37, Table 7],
see Figure 14
−

The concrete of the test specimens (named as M100 concrete in [36, 37]) is assumed to be a
high strength concrete, class 2 (see Eurocode 2 [17, Section 6.1] and Section 8), with calcareous
aggregates, which can be denoted as ”HSC2-C”. This classification determines the free thermal
strain of concrete and the temperature dependent reduction of its compressive strength according
to Eurocode 2 [17], as described in Section 8.
The reference compressive strength of concrete at the room temperature is taken from Kalifa
et al. [37, Table 1]. The reference tensile strength of concrete at the room temperature is not
specified in [36, 37] and hence, it is calculated by (81).
The mass of cement per unit volume of concrete is taken from Kalifa et al. [37, Table 1].
For the porosity of concrete, the parameters of equation (84) are determined by a linear regression of the data stated in [37, Table 2], see Figure 11.
Assuming the above values of c and φ(θ), the resulting degree of saturation with liquid water
(106), based on the sorption isotherms (98), is illustrated in Figure 12.
As mentioned in Section 8, the density of solid skeleton, ρs , is assumed to be a constant value
in our approach. This value can be estimated from the data measured by Kalifa et al. [37, Table 4]
for the apparent density. The corresponding apparent density assumed in our model, that can be
expressed as (cf. [24, eq. (26)])
ρ = ρw ηw + ρg ηg + ρs ηs = ρw (φ Sw ) + (ρv + ρa )[φ(1 − Sw )] + ρs (1 − φ),

(115)

appears in Figure 13.
The temperature dependence of the thermal conductivity of concrete was measured by Kalifa
et al. [37, Table 7] at a dry state. The parameters of equation (86) are determined by a linear
regression of the measured data, see Figure 14 (RH = 0). The resulting thermal conductivity of
concrete assumed in our model is shown in Figure 14.
The intrinsic permeability of concrete was tested by Kalifa et al. [37, Table 6] on concrete
samples dried at 105 ◦ C and then heated up to several temperature levels. However, the reference
21

0.14

1.2

0.8
Sw [−]

φ [−]

θ = 293.15
θ = 373.15
θ = 473.15
θ = 573.15

1.0

Test
Model

0.12
0.10

K
K
K
K

0.6
0.4

0.08

0.2
0.06
273.15

473.15

673.15
θ [K]

0

873.15

0

0.2

0.4

0.6 0.8
RH [−]

1

1.2

Figure 11: Porosity of concrete measured by Kalifa Figure 12: Degree of saturation with liquid water aset al. [37] (points) and assumed in our model (line)
sumed in our model

2.4

2550

ρ [kg m−3 ]

2500
2450

λc [W m−1 K−1 ]

Test
RH = 0
RH = 0.3
RH = 0.6
RH = 0.9

2400

Test
RH = 0
RH = 0.3
RH = 0.6
RH = 0.9

2.2
2.0
1.8
1.6
1.4

2350
273.15

473.15

673.15
θ [K]

1.2
273.15

873.15

473.15

673.15
θ [K]

873.15

Figure 13: Apparent density of concrete measured by Figure 14: Thermal conductivity of concrete measured
Kalifa et al. [37] at a dry state (points) and assumed in by Kalifa et al. [37] at a dry state (points) and assumed
our model (lines)
in our model (lines)

permeability of undamaged concrete at the room temperature is not reported in [37] and hence it
is determined by a trial-and-error method in our simulation (cf. [57, p. 277]).
Scheme of the analysed problem is displayed in Figure 15. The spatial discretization is performed
with the use of linear 1-D elements. In total, 120 elements are employed – 30 elements in the interval
x ∈ (0, `/2)), 30 elements for x ∈ (`/2, 3`/4) and 60 elements in the interval x ∈ (3`/4, `). For the
time discretization, the time step is set to ∆t = 1 s. The characteristic time of spalling is assumed
as γ = 10 s.
x
θ∞

Δx(t)

Radiant
heater

(t)

θ∞ (t)

P∞

P∞
Figure 15: Scheme of the analysed problem

The initial conditions are assumed as P0 = 1.9039 × 103 Pa, θ0 = 293.15 K, `0 = 0.12 m. For
these values, the initial saturation with liquid water Sw,0 = 0.77, which is the value reported by
Kalifa et al. [36, Table 2].
As stated above, a radiant heater of a temperature of 600 ◦ C was employed for the experiment.
22

However, as mentioned by e.g. Gawin et al. [23, pp. 554–556], the temperature of the ambient air
on the heated side, which needs to be defined for the boundary conditions assumed in the model
(see Section 4), is lower than the temperature of the heater. For our simulation, the ambient
temperature on the heated side, θ∞ [K], has been determined by a trial-and-error method as (see
Figure 16)


 293.15 + t 410
for t ≤ 300 s,
300
θ∞ (t) =
(116)
35

 703.15 + (t − 300)
for t > 300 s,
21300
where t [s] is the time of heating.
773.15
673.15
θ∞ [K]

573.15
473.15
373.15
273.15

0

1

2

3
4
t [hours]

5

6

Figure 16: Air temperature on the heated side of the specimen used for the simulation of the PTM test 1

The boundary conditions are summarized in Table 2.
Table 2: Boundary conditions parameters for simulation of PTM test 1

Variable
P∞
θ∞
αc
eσSB
βc

Value
Unexposed side
Exposed side
P0
P0
θ0
θ∞ (t), see (116)
4
20
0.7 × 5.67 × 10−8 0.7 × 5.67 × 10−8
0.009
0.019

Unit
Pa
K
W m−2 K−1
W m−2 K−4
m s−1

9.1.3. Discussion
The resulting pore pressure, temperature and mass loss evolutions determined by our simulation
are compared with the experimentally measured data in Figures 17–19 (note that the distances in
millimeters stated in the figures are, in accordance with Kalifa et al. [36], measured from the heated
surface).
The pore pressure evolution in the depth of 30 mm from the heated surface is not displayed in
Figure 17 since it was not correctly measured by Kalifa et al. [36], probably due to the damage of
the gauge (see [36, p. 1923]). On the other hand, the pore pressure evolution in the depths of (10,
20, and 50) mm can be verified more precisely because the investigated PTM test was duplicated
and the data from the both measurements are stated in [36, Figure 6], see Figures 20–22.
From Figures 17–22, it is obvious that in this case, despite of the simplifications mentioned in
the introductory part of this section, the present model predicts the hygro-thermal behaviour of
the heated concrete specimen on the sufficient level of accuracy. It should be however mentioned
that some of the potential inaccuracies arising from the simplifications embedded in the model are
23

4
3.5

40 mm
20 mm

P [MPa]

3
2.5

50 mm
Test
Model

10 mm

2
1.5
1
0.5
0
0

1

2

3
t [hours]

4

5

6

Figure 17: Pore pressure evolution measured by Kalifa et al. [36, Figure 4] (dashed lines) and determined by the
present model (solid lines)

probably eliminated by the appropriate setting of the reference permeability of undamaged concrete
at the room temperature (fitted by a trial-and-error method).
The spalling behaviour is also simulated correctly since, as mentioned above, during the test, no
spalling was observed by Kalifa et al. [36], which is in accordance with our simulation, see Figure 23
– the maximal value of failure function (31) achieved within the specimen during the test period
does not exceeds the value of 1. The evolution of the maximal values of the damage parameter and
its components within the specimen during the test period determined by our model is shown in
Figure 24.
For illustration, the spatial and time distributions of the primary unknowns of the present model
(except the thickness of the specimen, which is constant in this case) for the Kalifa’s PTM test 1
[36] are shown in Figures 25 and 26.

24

773.15

573.15

2 mm
10 mm
20 mm
30 mm
40 mm
50 mm

473.15

120 mm

θ [K]

673.15

373.15
Test
Model
273.15

0

1

2

3
t [hours]

4

5

6

Figure 18: Temperature evolution measured by Kalifa et al. [36, Figure 4] (dashed lines) and determined by the
present model (solid lines)

4

-1

3
Test
Model

-2

P [MPa]

Mass loss [%]

0

-3
-4

Test 1
Test 2
Model

2
1

0

2

t [hours]

4

0
0

6

2

t [hours]

4

6

Figure 19: Mass loss measured by Kalifa et al. [36, Figure 20: Pore pressure development measured by
Figure 7]) (dashed line) and determined by the present Kalifa et al. [36, Figure 6] and determined by the
model (solid line)
present model at the depth of 10 mm from the heated
surface
4

3
P [MPa]

3
P [MPa]

4

Test 1
Test 2
Model

2
1
0
0

Test 1
Test 2
Model

2
1

2

t [hours]

4

0
0

6

Figure 21: Pore pressure development measured by
Kalifa et al. [36, Figure 6] and determined by the
present model at the depth of 20 mm from the heated
surface

2

t [hours]

4

6

Figure 22: Pore pressure development measured by
Kalifa et al. [36, Figure 6] and determined by the
present model at the depth of 50 mm from the heated
surface

25

1

1
0.8
max(D) [−]

max(F ) [−]

0.8
0.6
0.4
0.2

0.6
0.4

max(Dθ )
max(Dm )
max(D)

0.2

0
0

2

t [hours]

4

0
0

6

Figure 23: Evolution of the failure parameter

4

1
0
0

60
90
x [mm]

1
0
6

120
90

0.2 h
0.5 h
1h
2h
3h
4h
5h
6h

θ [K]
473.15

60

2
t [hours]

673.15

273.15

2

120

873.15

30
0 0

30

60
90
x [mm]

x [mm]

θ [K]
873.15
673.15
473.15
273.15

120

6

90

4

0

6

3

4

30

4

4
P [MPa]

P [MPa]

2

t [hours]

Figure 24: Evolution of the damage parameter

0.2 h
0.5 h
1h
2h
3h
4h
5h
6h

3

2

60

2

120

t [hours]

30
0 0

Figure 25: Distribution of P and θ for the analysed PTM test 1

26

x [mm]

0.2 h
0.5 h
1h
2h
3h
4h
5h
6h

50
25
0
0

md [kg m−3 ]

75
50
25
0
6

120
90

4

30

60
90
x [mm]

t [hours]

0.2 h
0.5 h
1h
2h
3h
4h
5h
6h

9
6
3

60

2

120

12

0
0

100
m [kg m−3 ]

75

30
0 0

9
6
3
0
6

120
90

4

30

60
90
x [mm]

x [mm]

12
md [kg m−3 ]

m [kg m−3 ]

100

60

2

120

t [hours]

30
0 0

Figure 26: Distribution of m and md for the analysed PTM test 1

27

x [mm]

9.2. Simulation of the PTM test 2
9.2.1. Experiment description
The PTM test procedure proposed by Kalifa et al. [36] was adopted by Mindeguia as a part of
an extended experimental program reported in [44, 45, 46, 47]. Within this program, the material
properties investigation, hygro-thermal behaviour measurements (the PTM tests) and the spalling
experiments were performed (a detailed summary of the experiments included in the program
appears in [44, Section 1]). Hereafter, one of the PTM tests performed by Mindeguia [44] is closely
investigated – the PTM test of a high strength concrete (denoted as ”B60” in [44]) prismatic
specimens (conditioned in ”Air 1” according to Mindeguia [44]) heated by a radiant heater of the
temperature of 600 ◦ C (the type of heating denoted as ”modéré” in [44]) placed near its surface for
a heating period of 5 hours. The test set-up is the same as described for the Kalifa’s experiment
(Section 9.1.1). Note that (i) within the investigated PTM test 2, no spalling was observed by
Mindeguia [44], and (ii) the test was duplicated and the data from the both measurements of P
are stated in [44], see Section 9.2.3.
9.2.2. Modelling
The thermal and hygral properties of concrete are assumed according to the relationships given
in Section 8, with the parameters determined from the data measured by Mindeguia [44]. All the
parameters are summarized in Table 3.
Table 3: Material properties and parameters used in our simulation of the PTM test 2

Parameter

Value

Unit

Type of concrete

HSC1-C

−

fc,ref
ft,ref
c
θref
φref
Aφ

61.0
3.76
550
293.15
0.1027
1.0624 × 10−4

MPa
MPa
kg m−3
K
−
K−1

ρs

2660

kg m−3

θref
λd,ref
Aλ
Kref

293.15
2.0153
−9.8533 × 10−4
4.0 × 10−20

K
W m−1 K−1
K−1
m2

Reference
Eurocode 2 [17];
Mindeguia [44]
Mindeguia [44, Table 24]
Mindeguia [44, Table
determined from
Mindeguia [44, Table
see Figure 27
determined from
Mindeguia [44, Table
see Figure 28
determined from
Mindeguia [44, Table
see Figure 29
−

14]
25],

26],

28],

The concrete of the test specimens (denoted as C60 concrete in [44]) is assumed to be a high
strength concrete, class 1 (see Eurocode 2 [17, Section 6.1] and Section 8), with calcareous aggregates, which can be denoted as ”HSC1-C”. This classification determines the free thermal strain
of concrete and the temperature dependent reduction of its compressive strength according to
Eurocode 2 [17], as described in Section 8.
The reference compressive and tensile strengths of concrete at the room temperature are taken
from Mindeguia [44, Table 24].
The mass of cement per unit volume of concrete is taken from Mindeguia [44, Table 14].
For the porosity of concrete, the parameters of equation (84) are determined by a linear regression of the data stated in [44, Table 24], see Figure 27.
The density of solid skeleton, ρs , is estimated from the data measured by Mindeguia [44, Table 26]. The corresponding apparent density assumed in our model appears in Figure 28.
The parameters of equation (86) for the thermal conductivity of dry concrete are determined
by a linear regression of the data stated in [44, Table 28] (excluding the value for the temperature
28

0.18

2500

ρ [kg m−3 ]

φ [−]

0.16
0.14
Test
Model

0.12
0.10
273.15

Test
RH = 0
RH = 0.3
RH = 0.6
RH = 0.9

2450

473.15

673.15
θ [K]

2400
2350
2300
2250
273.15

873.15

473.15

673.15
θ [K]

873.15

Figure 27: Porosity of concrete measured by Mindeguia Figure 28: Apparent density of concrete measured by
[44] (points) and assumed in our model (line)
Mindeguia [44] (points) and assumed in our model
(lines)

2.6
2.4
2.2
2.0
1.8
1.6
1.4
1.2
1.0
0.8
273.15

773.15
Test
RH = 0
RH = 0.3
RH = 0.6
RH = 0.9

673.15
573.15

θ∞ [K]

λc [W m−1 K−1 ]

of 20 ◦ C), see Figure 29 (RH = 0). The resulting thermal conductivity of concrete assumed in our
model is shown in Figure 29.

473.15
373.15

473.15

673.15
θ [K]

273.15

873.15

0

1

2
3
t [hours]

4

5

Figure 29: Thermal conductivity of concrete measured Figure 30: Air temperature on the heated side of the
by Mindeguia [44] (points) and assumed in our model specimen used for the simulation of the PTM test 2
(lines)

The intrinsic permeability of concrete was tested by Mindeguia [44] on concrete samples dried
at 80 ◦ C and then heated up to several temperature levels. However, the reference permeability of
undamaged concrete at the room temperature is not reported in [44] and hence it is determined by
a trial-and-error method in our simulation (cf. [57, p. 277]).
The scheme of the analysed problem, the spatial discretization, the time discretization, and the
characteristic time of spalling are identical as in the previous example (see Section 9.1.2).
The initial conditions are assumed as P0 = 1.9194 × 103 Pa, θ0 = 293.15 K, `0 = 0.12 m. For
these values, the initial saturation with liquid water Sw,0 = 0.78, which is the value reported by
Mindeguia [44, Table 18].
The ambient temperature on the heated side, θ∞ [K], has been determined by a trial-and-error
method as (see Figure 30)


 293.15 + t 380
for t ≤ 300 s,
300
(117)
θ∞ (t) =
50

 673.15 + (t − 300)
for t > 300 s,
17700

29

where t [s] is the time of heating, and it is practically the same as for the Kalifa’s PTM test 1, cf.
Figures 16 and 30.
The boundary conditions are summarized in Table 4.
Table 4: Boundary conditions parameters for simulation of PTM test 2

Variable
P∞
θ∞
αc
eσ
βc

Value
Unexposed side
Exposed side
P0
P0
θ0
θ∞ (t), see (117)
4
20
0.7 × 5.67 × 10−8 0.7 × 5.67 × 10−8
0.009
0.019

Unit
Pa
K
W m−2 K−1
W m−2 K−4
m s−1

9.2.3. Discussion
The resulting pore pressure and temperature evolutions determined by our simulation are compared with the experimentally measured data in Figures 31 and 32 (note that the distances in
millimeters stated in the figures are, in accordance with Mindeguia [44], measured from the heated
surface). The more detailed illustration of the pore pressure prediction and its comparison with
the test data is shown in Figure 33.
3.5
3

Test 1
Test 2
Model
30 mm

P [MPa]

2.5
2

40 mm 50 mm

20 mm
10 mm

1.5
1
0.5
0
0

1

2

t [hours]

3

4

5

Figure 31: Pore pressure evolution measured by Mindeguia [44, Figure 154c] (dashed lines, dotted lines) and determined by the present model (solid lines)

From Figures 31–33, it is obvious that also in this case, the present model provides an accurate
prediction of the hygro-thermal behaviour of the heated concrete specimen. Similarly as in the
previous case, the potential inaccuracies arising from the simplifications of the model (see the
introductory part of this section) are probably excluded by the appropriate setting of the reference
permeability of undamaged concrete at the room temperature (fitted by a trial-and-error method).
The spalling behaviour is also simulated correctly since, as mentioned above, during the test, no
spalling was observed by Mindeguia [44], which is in accordance with our simulation, see Figure 34
– the maximal value of failure function (31) achieved within the specimen during the test period
does not exceeds the value of 1. The evolution of the maximal values of the damage parameter and
its components within the specimen during the test period determined by our model is shown in
Figure 35.

30

673.15

10
20
30
40
50

θ [K]

573.15

mm
mm
mm
mm
mm

473.15
120 mm
373.15
Test
Model
273.15

0

1

2

t [hours]

3

4

5

Figure 32: Temperature evolution measured by Mindeguia [44, Figure 157c] (dashed lines) and determined by the
present model (solid lines)

9.3. Simulation of the spalling test
9.3.1. Experiment description
The spalling experiments were performed by Mindeguia [44] on the 700×600×150 mm3 prismatic
specimens made of various types of concrete and exposed to various fire scenarios. Hereafter, the
spalling test of a high strength concrete (denoted as ”B60” in [44]) sample (conditioned in ”Air 2”
according to Mindeguia [44]) heated by the ISO 834 fire (in the 150-mm direction) for 60 minutes.
This test was duplicated and the data from the both measurements of the spalling depths are stated
in [44]. The specimen was placed on the top of the furnace and it was exposed to heating on the
area of 420 × 600 mm2 (see [44, Figure 170]). As mentioned by Mindeguia [44], Mindeguia et al.
[46], during the test, the specimen was not subjected to any mechanical load.
Within the test, the temperature and the pore pressure propagations through the specimen were
recorded. After the termination of heating, the spalling depths on the exposed surface were closely
measured which enables to determine the maximal and average depths of spalling and depict the
well-arranged ”spalling maps” (see [44, 46]).
The detailed description of the spalling experiments performed by Mindeguia can be found in
[44, 46].
9.3.2. Modelling
The thermal and hygral properties of the B60 concrete as well as its hygro-thermal behaviour
are described in Section 9.2.2.
Scheme of the analysed problem is displayed in Figure 36. The spatial discretization is performed
with the use of linear 1-D elements. In total, 160 elements are employed – 40 elements in the interval
x ∈ (0, `/2)), 40 elements for x ∈ (`/2, 3`/4) and 80 elements in the interval x ∈ (3`/4, `). For
the time discretization as well for the characteristic time of spalling, three values are assumed:
∆t = 1 s, ∆t = 0.5 s, ∆t = 0.1 s, and γ = 1 s, γ = 10 s, and γ = 100 s, respectively, which leads
to 9 numerical experiments discussed in Section 9.3.3.
The initial conditions are assumed as P0 = 1.9194 × 103 Pa, θ0 = 293.15 K, `0 = 0.15 m. For
these values, the initial saturation with liquid water Sw,0 = 0.78, which is the average value of the
two initial saturations reported by Mindeguia [44, Table 20].
The ambient temperature on the heated side is governed by the ISO 834 fire curve
θ∞ (t) = 293.15 + 345 log(8t + 1),
31

(118)

(10 mm)

3

2
1.5
1
0.5

2
1.5
1
0.5

0
0

2
t [hours]

0
0

4

(40 mm)

3

2.5
P [MPa]

2
1.5
1

Test 1
Test 2
Model

0.5
0
0

2
t [hours]

2
t [hours]

2

4

(50 mm)

3

2.5
P [MPa]

Test 1
Test 2
Model

2.5

Test 1
Test 2
Model

P [MPa]

P [MPa]

2.5

(20 mm)

3

Test 1
Test 2
Model

1.5
1
0.5
0
0

4

2
t [hours]

4

Figure 33: Pore pressure evolution measured by Mindeguia [44, Figure 157c] and determined by the present model
at the depth of (10, 20, 40, 50) mm from the heated surface

where t [minutes] is the time of heating.
The boundary conditions are summarized in Table 5.
Table 5: Boundary conditions parameters for simulation of spalling test

Variable
P∞
θ∞
αc
eσ
βc

Value
Unexposed side
Exposed side
P0
P0
θ0
θ∞ (t), see (118)
4
25
0.7 × 5.67 × 10−8 0.7 × 5.67 × 10−8
0.009
0.019

Unit
Pa
K
W m−2 K−1
W m−2 K−4
m s−1

9.3.3. Discussion
In this case, the spalling of concrete occurred during the test, which is also correctly simulated
by the present model. The reduction of the specimen thickness obtained by the calculations is
shown in Figure 37.
The spatial and time distributions of the primary unknowns of the present model are shown in
Figure 38.
It is obvious that the model predicted a large amount of concrete spalling – the spalling depth
at the end of heating reaches the values from 85 mm to 99 mm, depending on ∆t and γ assumed
for the calculation. These values of spalling depths are about 2.5 to 4 times as high as measured

32

1

1
0.8
max(D) [−]

max (F ) [−]

0.8
0.6
0.4
0.2
0
0

0.6
0.4

max(Dθ )
max(Dm )
max(D)

0.2
2
t [hours]

0
0

4

Figure 34: Evolution of the failure parameter

4

Figure 35: Evolution of the damage parameter

x
θ∞

2
t [hours]

Δx(t)

ISO 834
ﬁre

(t)

θ∞ (t)

P∞

P∞
Figure 36: Scheme of the analysed problem

by Mindeguia [44, Table 36] the maximal spalling depths of 34 mm and 23 mm for the first and for
the second test are reported, respectively.
Moreover, the values of pore pressure determined by the simulation are also significantly higher
than measured within the experiment. The maximal value of pore pressure obtained by the simulation is about 2.5 MPa (see Figure 38); the maximal value of pore pressure measured by Mindeguia
is 0.63 MPa (see [44, Table 38]).
The overestimation of pore pressures indicates that in such a rapid heating conditions, the
specimen was damaged more severely (which leads to increase of the permeability of concrete) than
has been predicted by the model. The fact that the increasing heating rate leads to decreasing
pore pressures within the spacemen was observed by Mindeguia et al. [45, Section 3.2]. In our
simulations, this trend is not so evident. When comparing the pore pressure distributions within
the specimen heated by a radiant heater of a temperature of 600 ◦ C (see Figures 31, 33) with
the pore pressure in the specimen made of the same concrete exposed to the ISO 834 fire (see
Figure 38), we get practically identical values (about 2.5 MPa). In our opinion, there are three
possible explanations of this fact:
• the Bary function (87) adopted in our model from Davie et al. [14, eq. (13)] is not valid for
such a rapid heating conditions as the ISO 834 fire exposure;
• the mechanical damage parameter influencing the increase of concrete permeability cannot
be simply assumed to be equal to the failure parametr assumed in our spalling criterion, as
we supposed in equation (89);
• the value of reference permeability of undamaged concrete at the room temperature that has
been set in our simulation of the spalling test is not correct. It has been adopted from the
previous PTM test of small specimen made of the same concrete. Within the PTM test,
concrete permeability has been fitted by a trial-and error method and hence, it is possible
that it covers indirectly some of the factors related to the PTM experiment of small samples.
However, these factors don’t need to be relevant for the spalling test of large specimens.
33

150

 [mm]

120
90
60
30
0
0

Δt = 1 s, γ = 1 s
Δt = 1 s, γ = 10 s
Δt = 1 s, γ = 100 s
Δt = 0.5 s, γ = 1 s
Δt = 0.5 s, γ = 10 s
Δt = 0.5 s, γ = 100 s
Δt = 0.1 s, γ = 1 s
Δt = 0.1 s, γ = 10 s
Δt = 0.1 s, γ = 100 s
10

20

30
t [minutes]

150

 [mm]

 [mm]

149.5
149
148.5
148
16

16.5

17
17.5
t [minutes]

18

40

70
68
66
64
62
60
58
56
54
52
50
55

50

56

60

57
58
t [minutes]

59

60

Figure 37: Reduction of the specimen thickness due to concrete spalling for the analysed experiment [44] obtained
by the present model

The detailed investigation of this phenomenon needs further research and it is out of scope of
the present paper.
The overestimation of the specimen thickness reduction is consequently caused by the overestimation of pore pressures as discussed in the previous paragraph, and probably also by the
assumption adopted in the model that the specimen is fully mechanically restrained in the plane
perpendicular to its thickness. In the investigated experiment, the only restraint arose from the
fixing the heated concrete by the unheated lateral parts and upper layers of the specimen (see [32]
and references therein). It should be noted that the overestimation of the spalling process has
been supposed since the assumption of the fully restrained conditions is generally considered as
conservative. Moreover, in the real conditions, the structure usually is mechanically restrained and
hence, in such a case, the model can provide more accurate prediction.
Nevertheless, numerical experiments demonstrate a computational stability of the present model.
From Figure 37, it is obvious that the model is not extremely sensitive on the value of time step.
It is also evident that by increasing the characteristic time of spalling as well as by decreing the
time step, the spalling proces becomes smoother.
10. Conclusions
In the paper, we proposed a one-dimensional coupled model for simulation of hygro-thermomechanical behaviour of concrete walls exposed to high temperatures including the prediction of
concrete spalling. The transport-processes-related part of the model was derived from the multi34

P [MPa]
2.5

150

90

1.5

60

1

30

0.5

0
0

10

20 30 40
t [minutes]

50

120

40
20

20 30 40
t [minutes]

50

573.15

60

473.15

30

373.15
10

20 30 40
t [minutes]

50

60

md [kg m−3 ]
0.8
0.6

60

10

90

120

90

0
0

673.15

150

60

30

120

0
0

60

m [kg m−3 ]
80

150

x [mm]

x [mm]

2

x [mm]

x [mm]

120

θ [K]

150

90
0.4

60

0.2

30
0
0

60

10

20 30 40
t [minutes]

50

60

0

Figure 38: Spatial and time distribution of the primary unknowns for the analysed experiment [44] obtained by the
present model (for ∆t = 0.1 s and γ = 1 s)

phase formulations of the laws of conservation of mass and energy. Based on the detailed quantitative parameter analysis, the less important transport phenomena were neglected in order to obtain
a simplified but still robust and realistic model of transport processes in heated concrete, which can
be expressed in terms of relatively small number of material parameters that can be easily obtained
from material tests. The main simplifications employed in our approach consisted in (i) neglecting
the diffusive mass flux of water vapour, (ii) removing the separate term describing the effect of
diffusion of adsorbed water and including this effect into the liquid water relative permeability,
Krw , and (iii) ignoring the effects of variations of pressure of dry air. The permeability of concrete
was supposed to be influenced both by the thermal and mechanical damage of the material, which
can be expressed by a multiplicative total damage parameter.
The spalling of concrete was assumed to be caused by a combination of the hygro-thermal stress
due to the pore pressure build-up and the thermo-mechanical stress resulting from the restrained
thermal dilatation. We employed a simplified mechanical approach based on the assumption that
the wall is fully mechanically restrained and the effects of stresses resulting from the external
mechanical load as well as the self-weight of the wall on the potential spalling are neglected.
An evolution law for moving boundary due to spalling was proposed. It enables to simulate the
instantaneous spalling of concrete as a continuous process, which contributes to the computational
stability of the numerical solution.
For the model represented by a system of partial differential equations and appropriate boundary
and initial conditions, the finite element discretization in space and the semi-implicit discretization
in time were employed. The resulting numerical algorithm was incorporated into an in-house
MATLAB code that was applied for numerical experiments in order to validate the present model.

35

The validity of the model was verified by comparing the results obtained by the numerical
simulations with the data measured within real high-temperature experiments. Two types of tests
were investigated – the tests of hygro-thermal behaviour of high-strength concrete specimens 300 ×
300 × 120 mm3 heated by a radiant heater of temperature of 600 ◦ C performed by Kalifa et al. [36]
and Mindeguia [44] and the spalling test of high-strength concrete specimens 700 × 600 × 150 mm3
exposed to the ISO 834 fire performed by Mindeguia [44]. All the material properties assumed
for the simulations, except the reference intrinsic permeability of undamaged concrete at the room
temperature, Kref , were adopted from the data measured by Kalifa et al. [36] and Mindeguia [44],
respectively. The values of Kref were not reported by the authors of the tests (Kalifa et al. [36]
and Mindeguia [44] reported only the values of the intrinsic permeability dried at 105 ◦ C or 80 ◦ C,
respectively, and then heated up to several temperature levels).
By comparing the calculated and measured data, the following findings can be drawn:
• the model is able to simulate the hygro-thermal behaviour of heated concrete on an appropriate level of accuracy, the temperature and pore pressure distributions obtained by the present
model are in close agreement with the data measured by Kalifa et al. [36] and Mindeguia [44];
• for the investigated examples, the model provided a correct qualitative prediction of concrete
spalling (whether it occurred or not);
• the quantitative prediction of spalling for the Mindeguia’s spalling test [44] was overestimated
(the values of spalling depths determined by the model at the end of the test period are
about 2.5 to 4 times as high as measured by Mindeguia [44]). This result arises from the
overestimation of pore pressures (as discussed in detail in Section 9.3.3) and probably also
from the assumption of the full mechanical restraint of the specimen. This assumption has
been adopted in the model but it is probably not fulfilled in the investigated test setup by
Mindeguia [44]. It can be supposed that in real structures (with some level of restraint), the
model can provide more accurate prediction.
For future research, it could be recommended to focuss on:
• determination of the reference intrinsic permeability of undamaged concrete at the room
temperature;
• validation of the present model on a real structure exposed to fire or on the results of full-scale
fire experiments;
• extension of the present model for two-dimensional problems and improvement of the mechanical part of the model in order to describe various types of loading and restrain conditions.

Appendix A. Physical mechanisms of transport processes in concrete at high temperatures
In this Appendix we address, on the quantitative level, contributions of different transport
mechanisms in concrete exposed to ambient high-temperatures. To this goal, we will first summarize
modelling assumptions and state variables employed in the derivation of complete multi-phase
formulations, together with the dependence of material constants on the state variables. Then, by
expressing the thermodynamic fluxes in terms of the same state variables, we will obtain closedform expressions for the corresponding transport coefficients. Their relative importance will be
evaluated at the level of material point and these results will allow us to explore the domains of
applicability and the connection of the presented simplified model against detailed multi-phase
descriptions.
36

Appendix A.1. Conservation of mass
The local form of the macroscopic mass balance equation of the phase α is [24]
Dα ρα
+ ρα ∇ · v α = Mα .
Dt

(A.1)

ρα = ηα ρα

(A.2)

Here,
represents the phase averaged density, ηα [-] is the volume fraction of the α phase and ρα stands
for the intrinsic phase averaged density. Note that
X
ηα = 1.
(A.3)
α

Further, v α [m s−1 ] is the velocity of α phase and Mα is the volumetric mass source. Using (A.2)
the equation (A.1) can be rewritten as
∂(ηα ρα )
+ ∇ · (ηα ρα v α ) = Mα .
∂t

(A.4)

Based on multi-phase modelling, the moist concrete is considered as a multi-phase porous
material consisted of solid skeleton with pores filled by liquid water and gas, which is a mixture
of dry air and water vapour. Let us write the balance equation (A.4) for each phase of concrete
mixture:
liquid water conservation equation:
∂(ηw ρw )
∂me ∂md
+ ∇ · (ηw ρw v w ) = −
+
;
∂t
∂t
∂t

(A.5)

water vapour conservation equation:
∂(ηv ρv )
∂me
+ ∇ · (ηv ρv v v ) =
;
∂t
∂t

(A.6)

∂(ηa ρa )
+ ∇ · (ηa ρa v a ) = 0;
∂t

(A.7)

dry air conservation equation:

solid mass conservation equation:
∂(ηs ρs )
∂md
+ ∇ · (ηs ρs v s ) = −
,
∂t
∂t

(A.8)

where md [kg m−3 ] is the mass source term related to the dehydration process and me [kg m−3 ] is
the vapour mass source caused by the liquid water evaporation.
The mass balances of the liquid water and of the vapour, summed together to eliminate the
source term related to phase changes (evaporation or condensation), form the mass balance equation
of moisture (liquid water and vapour):
∂
∂md
(ηw ρw + ηv ρv ) + ∇ · (ηw ρw v w + ηv ρv v v ) =
.
∂t
∂t

37

(A.9)

Appendix A.2. Moisture flux
The total moisture flux in moist concrete can be expressed as
J M = ηw ρw v w + ηg ρv v v ,
| {z } | {z }
Jw
Jv

(A.10)

where J w represents the liquid (capillary and adsorbed) water mass flux and J v is the water vapour
mass flux.
The chosen primary state variables are temperature θ, water vapour pressure Pv and dry air
pressure Pa . We express the thermodynamic fluxes in terms of state variables θ, Pv and Pa to
obtain the corresponding transport coefficients related to different physical mechanisms of moisture
transfer.
Liquid water mass flux. Different mechanisms governing the liquid (free) water mass flux need to be
distinguished. In particular, above the solid saturation point the transport of liquid water consists of
capillary water flows driven by the capillary pressure gradient. Otherwise, the physically adsorbed
water flows diffuses due to the saturation gradient. Consequently, mathematical expression of liquid
water (adsorbed or capillary) mass flux reads as (cf. [12])


SB
KKrw
SB
Jw = − 1 −
ρw
∇(Pg − Pc ) −
ρw DB ∇SB ,
(A.11)
Sw
µw
Sw
where the degree of saturation with adsorbed water SB is defined as

for Sw ≤ Sssp ,
 Sw
SB =

Sssp for Sw > Sssp ,

(A.12)

Sw represents the degree of saturation with liquid water and Sssp is the solid saturation point. Db
is the bound water diffusion tensor (see [21, eq. (39)]).
Note that the saturation degree with liquid water Sw is defined as
Sw =

ηw
φ

[−].

(A.13)

Here the volume fraction of liquid water ηw in concrete can be calculated via the sorption isotherms
(introduced by Bažant et al. [4, 5, 6]) as a function of temperature θ and vapour pressure Pv .
The equilibrium state of the capillary water with the water vapour is expressed in the form
corresponding to the Kelvin equation
 
Rθ
Pv
Pc (Pv , θ) = −ρw
ln
,
(A.14)
Mw
Ps
where Pc denotes the capillary pressure and the water vapour saturation pressure Ps can be calculated from the following formula as a function of temperature θ [K]


4042.9
.
(A.15)
Ps (θ) = exp 23.5771 −
θ − 37.58
Finally, the liquid water mass fluxes can be expressed in the closed-forms with respect to
gradients of primary state variables:
Adsorbed water diffusion (Sw ≤ Sssp , SB = Sw ):
ηw ρw v w = −ρw DB ∇Sw
∂Sw
∂Sw
= −ρw DB
∇θ;
∇Pv − ρw DB
∂Pv
∂θ
38

(A.16)

Capillary water flow (Sw > Sssp , SB = Sssp ):




KKrw
KKrw
SB
SB
ρw
∇Pg + 1 −
ρw
∇Pc
ηw ρw v w = − 1 −
Sw
µw
Sw
µw






SB
KKrw
SB
KKrw ∂Pc
∂Pc
=− 1−
ρw
∇Pg + 1 −
ρw
∇Pv +
∇θ
Sw
µw
Sw
µw
∂Pv
∂θ



SB
KKrw
KKrw ∂Pc
= 1−
−ρw
+ ρw
∇Pv
Sw
µw
µw ∂Pv




KKrw
SB
KKrw ∂Pc
SB
ρw
∇Pa + 1 −
ρw
− 1−
∇θ.
Sw
µw
Sw
µw ∂θ

(A.17)

Water vapour mass flux. Decomposition of the water vapour velocity into the diffusional (v v − v g )
and advectional (v g ) components yields
J v = ηg ρv v g + ηg ρv (v v − v g ).

(A.18)

The advective flux may be desribed by Darcy’s law in the form (vapour flow)
KKrg
∇Pg
µg
KKrg
KKrg
∇Pv − ρv
∇Pa .
= −ρv
µg
µg

ηg ρv v g = −ρv

(A.19)

For the diffusive mass flux of water vapour, Fick’s law (vapour diffusion) is applied in the form
 
Ma Mw
Pv
v
J d = ηg ρv (v v − v g ) = −ρg
Def f ∇
2
Mg
Pg
Ma M w P a
Ma M w P v
=−
Def f ∇Pv +
Def f ∇Pa . (A.20)
θR(Pv Mw + Pa Ma )
θR(Pv Mw + Pa Ma )
Total moisture flux. The total moisture flux can be expressed in terms of gradients of state variables
Pv , Pa and θ to obtain

ηw ρw v w + ηg ρv v v






KKrg
SB
KKrw
∂Pc

− 1−
= −ρv
ρw
1−
µg
Sw
µw
∂Pv

|
{z
}
| {z }
liquid water f low

vapour f low


SB
∂SB
− ρ w DB
S
∂Pv
| w {z
}

adsorbed water dif f usion


Ma Mw P a

−
Def f  ∇Pv
θR(Pv Mw + Pa Ma )

|
{z
}
vapour dif f usion









KKrg
SB
KKrw
Ma M w P v


+ −ρv
ρw
Def f  ∇Pa
− 1−
+
µg
Sw
µw
θR(Pv Mw + Pa Ma )


{z
}
{z
}|
| {z } |
liquid water f low

vapour f low

vapour dif f usion







SB
KKrw ∂Pc

+ 1 −
ρw
Sw
µw ∂θ

|
{z
}
liquid water f low

SB
∂SB
− ρw DB
S
∂θ
| w {z
}



 ∇θ.


adsorbed water dif f usion

(A.21)
39

10
10

-15

II

-25

1

10
300

I

450

[s]

[s]

III
-5

0

10

-15

IV

-25

1

10
300

0.5
600
θ [K]

10

-5

RH [−]

450

0.5
600
θ [K]

0

RH [−]

Figure A.39: Moisture fluxes related to ∇Pv representing the different moisture transfer mechanisms – (I) vapour
flow, (II) vapour diffusion, (III) liquid water flow, (IV) adsorbed water diffusion

Figures A.39 and A.40 present the dependence of transport coefficients upon the temperature
and relative humidity at constant atmospheric pressure (Pa = 101325 Pa) related to the gradient of
vapour pressure ∇Pv . Material properties of concrete are taken from Gawin et al. [21]. The moisture
flux according to gradient of vapour pressure, as derived in (A.21), consists of Darcian vapour flow,
vapour diffusion and Darcian flow of capillary water (provided Sw > Sssp ) or diffusion of adsorbed
water (provided Sw ≤ Sssp ). Figures A.39 and A.40 show that with increasing temperature Darcian
vapour flow plays a dominant role (approximately above 200◦ C). This observation is in accordance
with Dwaikat and Kodur [15] where it is assumed that only the vapour transport takes place during
the temperature exposure of HSC. At lower temperatures, the adsorbed water diffusion or Darcian
water flow (depending on degree of saturation) are of importance. Moreover, at the temperatures
between 100 ∼ 200 ◦ C in the relatively small RH-range, approximately up to 10 %, the vapour
diffusion may be the dominant factor of moisture transport. However, following our observation
based on numerical experiments, neglecting vapour diffusion has negligible effects on the spatial
position of the amplitude of the peak and the peak value of the pore pressure.
In order to simplify the model as much as possible and minimize the number of input parameters,
against complete multi-phase formulation, where the adsorbed and capillary water are assumed
separately, we follow the approach introduced by Chung and Consolazio [10], where the diffusion of
adsorbed water is considered to be accounted for within the liquid water relative permeability Krw
(see [13]). The main assumption in our model (46)–(56) is that the moisture transport in HSC at
high temperature related to the dry air pressure gradient is negligible when compared to the other
causes, in particular, the vapour pressure and temperature gradients, respectively (see [15]).
Appendix A.3. Heat energy conservation equation
The local form of the macroscopic energy balance equation of the phase α is [24]
ρα cαp

Dα θα
+ ∇ · q αc = Qα + E α − H α Mα .
Dt

(A.22)

Here, cαp is the specific isobaric heat, q αc the heat flux due to conduction, Qα the volumetric heat
source, E α represents energy exchange with other phases and H α is the specific enthalpy of the
phase α.
Summing up the energy balances for all phases (α ≡ a, P
w, v, s) of the multi-phase system,
neglecting internal sources Qα and taking into consideration α E α = 0 one arrives at the final
form of the

40

10

-10

10

-15

10

-20

10

-25

10

-5

10

-10

10

-15

10

-20

10

-25

θ = 293.15 K

[s]

-5

0

0.2

I

III

II

IV

0.4
0.6
RH [−]

0.8

1

θ = 393.15 K

[s]

[s]
[s]

10

0

0.2

I

III

II

IV

0.4
0.6
RH [−]

0.8

1

θ = 493.15 K

10

-5

10

-10

10

-15

10

-20

10

-25

10

-5

10

-10

10

-15

10

-20

10

-25

0

0.2

I

III

II

IV

0.4
0.6
RH [−]

0.8

1

θ = 593.15 K

0

0.2

I

III

II

IV

0.4
0.6
RH [−]

0.8

1

Figure A.40: Moisture fluxes related to ∇Pv representing the different moisture transfer mechanisms – (I) vapour
flow, (II) vapour diffusion, (III) liquid water flow, (IV) adsorbed water diffusion

energy conservation equation for moist concrete as the multi-phase system:
(ρcp )

∂θ
∂me
∂md
= −∇ · q c − (ρcp v) · ∇θ −
he −
hd ,
∂t
∂t
∂t

(A.23)

where
v
a
s
(ρcp ) = cw
p ρw ηw + cp ρv ηv + cp ρa ηa + cp ρs ηs ,

(ρcp v) =
qc =

v
cw
p ρw ηw v w + cp ρv ηv v v
v
a
q sc + q w
c + qc + qc ,
v
w

+

cap ρa ηa v a

+

csp ρs ηs v s ,

he = H − H ,
w

s

hd = H − H .

(A.24)
(A.25)
(A.26)
(A.27)
(A.28)

v a
s
−1 K−1 ] are the specific heats at constant pressure of liquid water, water
Here, cw
p , cp , cp and cp [J kg
v
a
s
vapour, dry air and solid matrix, q w
c , q c , q c and q c represent heat fluxes due to conduction corresponding to liquid water, water vapour, dry air and solid matrix. Further, he [J kg−1 ] represents
the enthalpy of evaporation per unit mass, while hd [J kg−1 ] is the enthalpy of dehydration per unit
mass. Finally, H w , H v and H s [J kg−1 ] are specific enthalpies of liquid water, water vapour and
chemically bound water.
Neglecting the impact of diffusional flows of dry air and vapour in the gas mixture on heat
transfer by convection yields

cvp ρv ηv v v + cap ρa ηa v a = (cvp ρv ηv + cap ρa ηa )v g = cgp ρg ηg v g ,
where cgp is the specific heat at constant pressure of gas mixture.
41

(A.29)

Appendix B. Nomenclature
Symbol
αc
βc
tot
θ
σ
cr
tr
ηs
ηw
ηg
φ
ρs
ρg
ρa
ρv
ρv∞
ρw
µw
µg
σ
σht
σtm
σSB
θ
θ∞
λc
τ

Unit
[W m−2 K−1 ]
[m s−1 ]
[−]
[−]
[−]
[−]
[−]
[−]
[−]
[−]
[−]
[kg m−3 ]
[kg m−3 ]
[kg m−3 ]
[kg m−3 ]
[kg m−3 ]
[kg m−3 ]
[Pa s]
[Pa s]
[Pa]
[Pa]
[Pa]
[W m−2 K−1 ]
[K]
[K]
[W m−1 K−1 ]
[s]

c
cw
p
cvp
cap
csp
cgp
e
fc
ft
F
he
hd
Hw
Hv
Hs
Jw
Jv
JM
K

kg m−3
[J kg−1 K−1 ]
[J kg−1 K−1 ]
[J kg−1 K−1 ]
[J kg−1 K−1 ]
[J kg−1 K−1 ]
[−]
[Pa]
[Pa]
[−]
[J kg−1 ]
[J kg−1 ]
[J kg−1 ]
[J kg−1 ]
[J kg−1 ]
[kg m−2 s−1 ]
[kg m−2 s−1 ]
[kg m−2 s−1 ]
[m2 ]

Description
convective heat transfer coefficient
convective mass transfer coefficient
total strain
free thermal strain
instantaneous stress-related strain
creep strain
transient strain
volume fraction of the solid microstructure
volume fraction of liquid phase
volume fraction of gas phase
porosity
density of the solid microstructure
gas phase density
dry air phase density
vapour phase density
ambient vapour phase density
liquid phase density
liquid water dynamic viscosity
gas dynamic viscosity
stress
stress caused by hygro-thermal processes
stress caused by thermo-mechanical processes
Stefan-Boltzmann constant
absolute temperature
ambient absolute temperature
effective thermal conductivity of moist concrete
characteristic time of mass loss governing the asymptotic evolution of the dehydration process
mass of cement per m3 of concrete
specific heat at constant pressure of liquid water
specific heat at constant pressure of water vapour
specific heat at constant pressure of dry air
specific heat at constant pressure of solid phase
specific heat at constant pressure of gas mixture
emissivity of the interface
uniaxial compressive strength of concrete
uniaxial tensile strength of concrete
failure function
enthalpy of evaporation per unit mass
enthalpy of dehydration per unit mass
specific enthalpy of liquid water
specific enthalpy of water vapour
specific enthalpy of chemically bound water
liquid water mass flux
water vapour mass flux
moisture flux
intrinsic permeability
42

Krw
Krg
`
md
md,eq
me
Mw
Ma
P
Pc
Pv
Pa
Pg
Pw
Ps
qc
R
RH
Sw
SB
Sssp
vg
vw
va
vv

[−]
[−]
[m]
[kg m−3 ]
[kg m−3 ]
[kg m−3 ]
[kg kmol−1 ]
[kg kmol−1 ]
[Pa]
[Pa]
[Pa]
[Pa]
[Pa]
[Pa]
[Pa]
[W m−2 ]
[J kmol−1 K−1 ]
[−]
[−]
[−]
[−]
[m s−1 ]
[m s−1 ]
[m s−1 ]
[m s−1 ]

relative permeability of liquid water
relative permeability of gas
thickness of a concrete wall
mass source term related to the dehydration process
mass of water released at the equilibrium
vapour mass source caused by the liquid water evaporation
molar mass of liquid water
molar mass of dry air
pore pressure due to water vapour (P = Pv )
capillary pressure
water vapour pressure
dry air pressure
gas pressure
pressure of liquid water
water vapour saturation pressure
heat flux vector
gas constant
relative humidity (RH = Pv /Ps )
degree of saturation with liquid water
degree of saturation with adsorbed water
solid saturation point
velocity of gaseous phase
velocity of liquid phase
velocity of dry air
velocity of vapour

Acknowledgement. This research has been supported by the Czech Science Foundation, project
GAČR 13-18652S (the first and the second author), and by the Grant Agency of the Czech Technical
University in Prague, project SGS13/038/ OHK1/1T/11 (the second author). The support is
gratefully acknowledged.
References
References
[1] Alnajim, A., 2004. Modelisation et simulation du comportement du beton sous hautes temperatures par une approche thermo-hygro-mecanique couplee. application a des situations accidentelles. Ph.D. thesis, Universite de Merne la Vallee.
[2] Anderberg, Y., 2008. The impact of various material models on structural fire behaviour
prediction. In: Tan, K., Kodur, V., Tan, T. (Eds.), Proceedings of the Fifth International
Conference on Structures in Fire (SiF’08). Nanyang Technological Uniersity, Singapore, pp.
253–265.
[3] Anderberg, Y., Thelandersson, S., 1976. Stress and deformation characteristics of concrete at
high temperatures. 2. Experimental investigation and material behaviour model. Lund Institute of Technology.
[4] Bažant, Z., Chern, J.-C., Thonguthai, W., 1981. Finite element program for moisture and heat
transfer in heated concrete. Nucl. Eng. Des. 68 (1), 61–70.

43

[5] Bažant, Z., Thonguthai, W., 1978. Pore pressure and drying of concrete at high temperature.
Proc. ASCE, J. Eng. Mech. Div. 104 (5), 1059–1079.
[6] Bažant, Z., Thonguthai, W., 1979. Pore pressure in heated concrete walls: theoretical prediction. Mag. Concr. Res. 31 (107), 67–76.
[7] Beneš, M., Štefan, R., Zeman, J., 2013. Analysis of coupled transport phenomena in concrete
at elevated temperatures. Appl. Math. Comput. 219, 7262–7274.
[8] Bratina, S., Čas, B., Saje, M., Planinc, I., 2005. Numerical modelling of behaviour of reinforced
concrete columns in fire and comparison with eurocode 2. International J. Solids Struct. 42 (2122), 5715 – 5733.
[9] Carstensen, J., 2011. Material modelling of reinforced concrete at elevated temperatures. Master’s thesis, Technical University of Denmark.
[10] Chung, J., Consolazio, G., 2005. Numerical modeling of transport phenomena in reinforced
concrete exposed to elevated temperatures. Cement Concrete Res. 35 (3), 597–608.
[11] Dal-Pont, S., Ehrlacher, A., 2004. Numerical and experimental analysis of chemical dehydration, heat and mass transfers in a concrete hollow cylinder submitted to high temperatures.
Int. J. Heat Mass Tran. 47 (1), 135–147.
[12] Davie, C., Pearce, C., Bićanić, N., 2006. Coupled heat and moisture transport in concrete
at elevated temperatures-effects of capillary pressure and adsorbed water. Numer. Heat Tr.
A-Appl. 49 (8), 733–763.
[13] Davie, C., Pearce, C., Bićanić, N., 2010. A fully generalised, coupled, multi-phase, hygrothermo-mechanical model for concrete. Mater. Struct. 43 (1), 13–33.
[14] Davie, C., Pearce, C., Bićanić, N., 2012. Aspects of permeability in modelling of concrete
exposed to high temperatures. Transport Porous Med. 95 (3), 627–646.
[15] Dwaikat, M., Kodur, V., 2009. Hydrothermal model for predicting fire-induced spalling in
concrete structural systems. Fire Safety J. 44 (3), 425–434.
[16] Eurocode 2, 2004. Design of concrete structures – Part 1-1: General rules and rules for buildings. CEN.
[17] Eurocode 2, 2004. Design of concrete structures – Part 1-2: General rules – Structural fire
design. CEN.
[18] Fellinger, J., 2004. Shear and anchorage behaviour of fire exposed hollow core slabs. Ph.D.
thesis, TU Delft.
[19] Feraille-Fresnet, A., Tamagny, P., Ehrlacher, A., Sercombe, J., 2003. Thermo-hydro-chemical
modelling of a porous medium submitted to high temperature: An application to an axisymmetrical structure. Math. Comput. Model. 37 (5), 641–650.
[20] Furbish, D., 1997. Fluid Physics in Geology. An Introduction of Fluid Motions on Earth’s
Surface and Within Its Crust. Oxford University Press.
[21] Gawin, D., Majorana, C., Schrefler, B., 1999. Numerical analysis of hygro-thermal behaviour
and damage of concrete at high temperature. Mech. Cohes.-Frict. Mat. 4 (1), 37–74.
[22] Gawin, D., Pesavento, F., 2012. An overview of modeling cement based materials at elevated
temperatures with mechanics of multi-phase porous media. Fire Technol. 48 (3), 753–793.
44

[23] Gawin, D., Pesavento, F., Schrefler, B., 2002. Modelling of hygro-thermal behaviour and damage of concrete at temperature above the critical point of water. Int. J. Numer. Anal. Met.
26 (6), 537–562.
[24] Gawin, D., Pesavento, F., Schrefler, B., 2003. Modelling of hygro-thermal behaviour of concrete
at high temperature with thermo-chemical and mechanical material degradation. Comput.
Methods Appl. Mech. Engrg. 192 (13-14), 1731–1771.
[25] Gawin, D., Pesavento, F., Schrefler, B., 2006. Towards prediction of the thermal spalling risk
through a multi-phase porous media model of concrete. Comput. Methods Appl. Mech. Engrg.
195 (41-43), 57075729.
[26] Gawin, D., Pesavento, F., Schrefler, B., 2011. What physical phenomena can be neglected when
modelling concrete at high temperature? A comparative study. Part 1: Physical phenomena
and mathematical model. Int. J. Solids Struct. 48 (13), 1927–1944.
[27] Gawin, D., Pesavento, F., Schrefler, B., 2011. What physical phenomena can be neglected
when modelling concrete at high temperature? A comparative study. Part 2: Comparison
between models. Int. J. Solids Struct. 48 (13), 1945–1961.
[28] Gernay, T., Franssen, J.-M., 2010. Consideration of transient creep in the eurocode constitutive
model for concrete in the fire situation. In: Kodur, V., Franssen, J.-M. (Eds.), Proceedings of
the Sixth International Conference Structures in Fire. DEStech Publications, Lancaster, PA,
pp. 784–791.
[29] He, Z., Song, Y., 2008. Failure mode and constitutive model of plain high-strength highperformance concrete under biaxial compression after exposure to high temperatures. Acta
Mech. Solida Sin. 21 (2), 149–159.
[30] Ichikawa, Y., England, G., 2004. Prediction of moisture migration and pore pressure build-up
in concrete at high temperatures. Nucl. Eng. and Des. 228 (1-3), 245–259.
[31] Jansson, R., Boström, L., 2010. The influence of pressure in the pore system on fire spalling
of concrete. Fire Technol. 46 (1), 217–230.
[32] Jansson, R., Boström, L., 2013. Factors influencing fire spalling of self compacting concrete.
Mater. Struct. 46 (10), 1683–1694.
[33] Jirásek, M., Bažant, Z., 2002. Inelastic Analysis of Structures. Wiley.
[34] Ju, J., Zhang, Y., 1998. A thermomechanical model for airfield concrete pavement under
transient high temperature loadings. Int. J. Damage Mech. 7 (1), 24–46.
[35] Kalifa, P., Chene, G., Galle, C., 2001. High-temperature behaviour of hpc with polypropylene
fibres: From spalling to microstructure. Cement Concrete Res. 31 (10), 1487–1499.
[36] Kalifa, P., Menneteau, F., Quenard, D., 2000. Spalling and pore pressure in hpc at high
temperatures. Cement Concrete Res. 30 (12), 1915–1927.
[37] Kalifa, P., Tsimbrovska, M., Baroghel-Bouny, V., 1998. High-performance concrete at elevated
temperatures – an extensive experimental investigation of thermal and hygral properties and
microstructure. In: Aı̈tcin, P.-C., Delagrave, Y. (Eds.), International Symposium on HighPerformance and Reactive Powder Concretes. Université de Sherbrooke, Sherbrooke, Canada,
pp. 259–279.

45

[38] Khoury, A., Anderberg, Y., 2000. Fire safety design – Concrete spalling review. Swedish National Road Administration.
[39] Kukla, K., 2010. Concrete at high temperatures – hygro-thermo-mechanical degradation of
concrete. Ph.D. thesis, University of Glasgow.
[40] Le, T. T. H., 2011. Étude multi-échelle du comportement thermo-hydro-mécanique des
matériaux cimentaires. approche morphologique pour la prise en compte de la mésostructure.
Ph.D. thesis, Université Paris-Est.
[41] Li, L., Purkiss, J., 2005. Stress–strain constitutive equations of concrete material at elevated
temperatures. Fire Saf. J.l 40 (7), 669 – 686.
[42] Majumdar, P., Gupta, A., Marchertas, A., 1995. Moisture propagation and resulting stress in
heated concrete walls. Nucl. Eng. Des. 156 (1-2), 147–158.
[43] Menétrey, P., Willam, K., 1995. Triaxial failure criterion for concrete and its generalization.
ACI Struct. J. 92 (3), 311–318.
[44] Mindeguia, J.-C., 2009. Contribution expérimentale a la compréhension des risques d’instabilité
thermique des bétons. Ph.D. thesis, UPPA.
[45] Mindeguia, J.-C., Pimienta, P., Carré, H., Borderie, C. L., 2009. Experimental study on the
contribution of pore vapour pressure to the thermal instability risk of concrete. In: Dehn, F.,
Koenders, E. A. B. (Eds.), International Workshop on Concrete Spalling due to Fire Exposure.
MFPA Institute, Leipzig, Germany.
[46] Mindeguia, J.-C., Pimienta, P., Carré, H., Borderie, C. L., 2013. Experimental analysis of
concrete spalling due to fire exposure. Eur. J. Environ. Civ. En. 17 (6), 453–466.
[47] Mindeguia, J.-C., Pimienta, P., Noumowé, A., Kanema, M., 2010. Temperature, pore pressure
and mass variation of concrete subjected to high temperature – experimental and numerical
discussion on spalling risk. Cem. Concr. Res. 40 (3), 477–487.
[48] Msaad, Y., 2007. Comparison between hydraulic and thermal spalling in heated concrete based
on numerical modeling. J. Eng. Mech. 133 (6), 608615.
[49] Naus, D., 2010. A compilation of elevated temperature concrete material property data and
information for use in assessments of nuclear power plant reinforced concrete structures. Tech.
rep., Office of Nuclear Regulatory Research.
[50] Pesavento, F., 2000. Non-linear modelling of concrete as multiphase porous material in high
temperature conditions. Ph.D. thesis, Università degli Studi di Padova.
[51] Phan, M., 2012. Modélisation explicite de l’écaillage sous incendie du béton: Approche thermohydro-mécanique avec des conditions aux limites évolutives. Ph.D. thesis, Université Paris-Est.
[52] Reddy, J., 2004. An Introduction to Nonlinear Finite Element Analysis. Oxford University
Press, USA.
[53] Schrefler, B., Brunello, P., Gawin, D., Majorana, C., Pesavento, F., 2002. Concrete at high
temperature with application to tunnel fire. Comput. Mech. 29 (1), 43–51.
[54] Shekarchi, M., Debicki, G., Billard, Y., Coudert, L., 2003. Heat and mass transfer of high
performance concrete for reactor containment under severe accident conditions. Fire Technol.
39 (1), 63–71.
46

[55] Tenchev, R., Li, L., Purkiss, J., 2001. Finite element analysis of coupled heat and moisture
transfer in concrete subjected to fire. Numer. Heat Tr. A-Appl. 39 (7), 685–710.
[56] Terro, M., 1998. Numerical modeling of the behavior of concrete structures in fire. ACI Struct.
J. 95 (2), 183–192.
[57] Witek, A., Gawin, D., Pesavento, F., Schrefler, B., 2007. Finite element analysis of various
methods for protection of concrete structures against spalling during fire. Comput. Mech.
39 (3), 271–292.
[58] Zeiml, M., Lackner, R., Mang, H., 2008. Experimental insight into spalling behavior of concrete
tunnel linings under fire loading. Acta Geotechnica 3 (4), 295–308.
[59] Zhao, J., Zheng, J.-J., Peng, G.-F., van Breugel, K., 2014. A meso-level investigation into the
explosive spalling mechanism of high-performance concrete under fire exposure. Cem. Concr.
Res. 65 (0), 64–75.

47

