FINITE ELEMENT ANALYSIS OF COUPLED HEAT AND MOISTURE TRANSFER IN CONCRETE SUBJECTED TO FIRE

Prévia do material em texto

FINITE ELEM ENT ANALYSIS OF COUPLED HEAT AND
MOISTURE TRANSFER IN CONCRETE SUBJECTED TO
FIRE
R. T. Tenchev, L. Y . Li, and J. A. Purkiss
School of Engineering and Applied Science, Aston University, Aston Triangle,
Birmingham, B4 7ET, United Kingdom
A system of coupled transient differential equations governing heat, mass transfer, and
pore pressure built up in porous media ( concrete) , subjected to intensive heating, is
derived. Water vapor and liquid water are considered separately in the mass transfer
formulation. The primary unknowns are temperature, water vapor content, and pore press-
ure of the gaseous mixture. A � nite element formulation and corresponding � owchart of
computat ions of all required data are presented. The numerical example solved represents
a cross section of a concrete column exposed to � re. The domain and time distributions
of temperature, pore pressure, water vapor, and liquid water content are presented. Com-
puted pore pressure is higher than those usually reported by other analytical studies.
The in� uence of some initial parameters ( permeability, initial water content, and porosity)
on maximum pore pressure is investigated.
INTRODUCTION
Concrete as a structural material comprises three different phasesösolid,
liquid (water), and gaseous mixture (water vapor and air). The solid phase consists
of cement paste and aggregate, both of which are porous. The cement paste contains
very small gel pores (diameter about 2 nm) and capillary pores (about 1 mm). The
products of hydration in the cement paste contain chemically bound water. The
gel pores hold water absorbed on their large surface. The capillary pores may be
fully ¢lled (saturated concrete) or partially ¢lled with water. The gel and capillary
water is usually called free or evaporable water. The gaseous phase exists in the
partially ¢lled capillary pores.
When a concrete structure is exposed to high external temperature, heat is
transferred within the concrete by both conduction and convection. The free water
begins to evaporate. At higher temperatures the chemically bound water is released
by dehydration and then evaporates. Concrete has low permeability and when
subjected to extensive heat, as is the case during ¢re, the rate of evaporation exceeds
the rate of vapor migration and this results in pore pressure buildup. The transport of
Numerical Heat Transfer, Part A, 39:685^710 , 2001
Copyright # 2001 Taylor & Francis
1040-7782 /01 $12.00 + .00
685
Received 3 March 2000; accepted 2 January 2001.
This work is a part of an EPSRC funded project (grant GR/M13565) on investigaton of concrete
spalling. The authors gratefully acknowledge this support.
Address correspondence to Dr. Rosen T. Tenchev, School of Engineering and Applied Science,
Aston University, Aston Triangle, Birmingham B4 7ET, UK. E-mail: tenchrt@aston.ac.uk
vapor and water is governed by the pressure gradient (the dominant factor),
moisture, and temperature gradients. The vapor is transferred in two directions:
(i) outward, where it escapes; (ii) inward, where it condenses, saturating the concrete
and blocking further movement in that direction.
The combined effect of pore pressure and thermal stresses may cause spalling
and reduction of its ¢re resistance. It is a very serious problem, especially for high
strength concrete because of its lower permeability.
The occurrence of spalling in concrete has been reported in the technical litera-
ture [1, 2]. It also has been seen in real structures: the Storebaelt Tunnel in Denmark,
NOM ENCLATURE
Ci speci¢c heat of phase i,
[J/(kg ¯C)], i ˆ A; V ; S; L
DAV diffusion coef¢cient of air or
vapor in the gaseous mixture
[m2/s]
_EL rate of evaporation of liquid
water per unit volume of
concrete [kg/(m3s)]
hq convective heat transfer
coef¢cient [J/(m2 s ¯K]
hqr hqr ˆ hq ‡ hr
Ji mass £ux of phase i, per unit
area of concrete [kg/(m2s)],
i ˆ A; V ; L
k thermal conductivity of concrete
[J/(m2s¯C/m)]
K speci¢c (intrinsic) permeability
of dry concrete [m2]
Ki relative permeability of phase i,
through dry concrete [1%],
i ˆ G; L
Pi partial pressure of phase i,
i ˆ A; V ; pressure of phase i,
i ˆ G; L, [Pa]
PSat saturation pressure of water
vapor [Pa]
por initial porosity
p¤or increase of porosity with
temperature
Ri gas constant for phase i, i=A, V
T temperature [¯C] or [¯K]
vi velocity of phase i, [m/s],
i ˆ G; L
ei volume fractons of phase i, [%],
i ˆ G; L; S; ei ˆ Vi=V
(V -volumes)
lD speci¢c heat of dehydraton of
bound water [J/kg]
lE speci¢c heat of evaporation
[J/kg]
b water vapor transfer coef¢cient
on boundary [m/s]
mi (dynamic) viscosity of phase i
[kg/(ms)], i ˆ A; V ; G; L
ri density of phase i ˆ L; S,
~ri mass of phase i per unit volume
of gaseous mixture [kg/m3],
i ˆ A; V
ri mass of phase i per unit volume
of concrete [kg/m3], i ˆ L; D
rCem mass of anhydrous cement per
unit volume of concrete [kg/m3]
rC effective heat capacity of
concrete [J/(m3 ¯C)]
rCv vector of energy transport by
£uid £ow [J/(m2 ¯C)]
j external humidity
& (overbar) variable de¢ned per unit volume
of concrete
~& (tilde) variable de¢ned per unit volume
of gaseous mixture
Subscripts
A dry air within the gaseous
mixture
D released bound water by
dehydration
G gaseous phase (mixture of air
and water vapor)
L liquid water phase
S dry skeleton
V water vapor
1 surrounding environment
Superscripts
e element number
n number of time step (n ˆ 0
initial ambient state)
686 R. T. TENCHEV ET AL.
where the spalling reduced the concrete structure to a mere 25% of its original
thickness during a ¢re in 1994 [3]; the Channel Tunnel, where on November 18,
1996, a ¢re of 10 h duration destroyed parts of the concrete tunnel rings (depth
of 0.45 m) by thermal spalling over a length of a few thousand meters and the dam-
age reached an average depth of 0.1 to 0.2 m, [4, 5].
There are two main hypotheses explaining thermal spalling:
(i) Spalling due to pressure buildup. The hypothesis is supported by investi-
gations showing that spalling is closely related to moisture content, con-
crete permeability, and porosity [6, 7, 8]. It is opposed by the facts
that predried concrete is still susceptible to spalling (as is the case with
the 10-year-old rings in the Channel Tunnel) and that prevailing data,
experimentally determined or theoretically predicted, for the magnitude
of the buildup pressure does not exceed 1 MPa for normal strength con-
crete and 2 to 3 MPa for high performance concrete [4, 5].
(ii) Spalling due to restrained thermal dilatation [4, 5, 9]. The restrained
thermal dilatation causes high compressive stresses parallel to the heated
surface. These stresses are released by brittle cracking (i.e., spalling) along
45¯ planes. Once a crack is formed (i.e., there is rapid expansion of the
void that can be occupied by the evaporated water), the pore pressure
drops down, and it can only play a secondary role in the growth of a larger
crack. However, pore pressure is the dominant factor in the so-called
explosive spalling [3].
The in£uence of spalling on the ¢re resistance of concrete structures was inves-
tigated in [10]. The correlation of analytical methods with observed experimental
behavior of reinforced concrete columns subjected to ¢re could be improved if
spalling is taken into account [11]. The correct prediction of spalling requires a
reliable solution to the complex problem of heat and moisture transfer in concrete.
Heat and moisture transfer in porous materials is a coupled phenomenon. The
migration of moisture (liquid water and water vapor) augments the conductive heat
transfer by convective heat transfer. The pore pressure gradient developed through
the process of heat transfer and water evaporation is an additional driving force
of mass transfer. The phenomenon is transient and strongly nonlinear because most
material properties of the solid, liquid, and gaseous phases are temperature
dependent.
Luikov [12, 13] was one of the ¢rst to study in detail, analytically as well as
experimentally,the heat and moisture transfer in porous materials. He stated that
moisture transfer results from gradients of moisture content, temperature, and
pressure. On the bases of mass conservation equations for each phase (mass £uxes
given by Fick’s law), the energy conservation equation (heat £ux given by Fourier’s
law), and Darcy’s law he derived a system of three coupled differential equations.
In the derivation of these equations two basic assumptions were made: (i) there
is thermal equilibrium between all phases; (ii) the mass of the liquid is equal to
the total mass involved in the mass transport, that is, the mass of vapor and air
is negligible. Another shortcoming in the formulation is that in the case of laminated
materials composed of inhomogeneous moist bodies in contact there occurs moisture
content discontinuity at the contact surfaces even at equilibrium. Based on the
TRANSFER IN CONCRETE SUBJECTED TO FIRE 687
principles of irreversible thermodynamics he proposed that a mass-transfer
potential, similar to the heat-transfer one (temperature), should be used as the
driving force for moisture transfer, rather than the moisture content. A drawback
of both formulations is that a number of transfer mechanisms are lumped together,
masking the individual differences between mechanisms and the dependence of each
on different controlling variables. To distinguish the differences between vapor and
liquid transfer, two-zone models were proposed [13^15]. The governing differential
equations for high velocity mass transfer, when the mean temperatures of the
dry body and the liquid within an in¢nitesimal volume are not equal, also are given
in [13].
Philip and deVries [16] and deVries [17] used an approach similar to Luikov’s
one, based on moisture content and temperature gradients in the presentation of
£uxes, neglecting pressure gradient but including gravitational contribution.
Whitaker [18] developed a formal theory of drying starting from point
equations in each phase and obtained volume averaged conservation equations.
Although his formulation does not differ substantially from that of Luikov, it offers,
through its averaging procedure, the most substantiative justi¢cation for the con-
tinuum mechanics approach to the modeling of heat and mass transfer in porous
materials.
Bazant and Thonguthai [19^21], on the basis of Luikov’s theory, proposed ¢eld
equations for the coupled heat and moisture transfer in concrete at high
temperatures. Making use of experimental data they simpli¢ed the expression of
the mass £ux (a function only of the gradient of the pressure) and presented some
semiempirical relationshipsösorption isoterms, dependence of permeability on
temperature, change of porosity because of dehydration, and so on.
More complex formulations for heat and mass transfer in concrete are given in
[7, 23^25].
Once the heat and moisture transfer in porous materials has been formulated
successfully in terms of a coupled system of partial differential equations, the ana-
lytical solution of these equations presents a major problem. Analytical solutions
of Luikov’s equations are available only for some simple geometric con¢gurations
with simple boundary conditions and constant material properties [12, 26, 27].
For the solution of more realistic problems approximate numerical techniques often
are required.
The ¢nite difference method was used for the one dimensional analysis of heat
and mass transfer in concrete exposed to ¢re [23]. The governing equations were
based on mass conservation of liquid water, water vapor, and gaseous mixture (water
vapor and air) and an energy conservation equation. A system of three coupled
differential equations was derived with the prime variables being temperature,
pressure, and mole fraction of water vapor. The numerical example solved showed
that the maximum pore pressure was at a considerable distance (about 0.1 m) from
the exposed surface and that spalling was unlikely, but no data for the material
properties used was given.
The ¢nite element method formulation for Luikov’s equations was used in
[28^32] for the solution of a low-rate heat and moisture transfer, when the effect
of pore pressure was neglected. The effect of pore pressure and the deformability
of the porous media was taken into account in [33]. Finite element formulations
688 R. T. TENCHEV ET AL.
of the model proposed by Baz­ant and Thonguthai for the analysis of heat and
moisture transfer in concrete are given in [19^22, 34^38].
The purpose of this article is to derive a comprehensive model for heat and
moisture transfer in concrete subjected to high temperatures and to give the corre-
sponding three-dimensoional ¢nite element formulation. Its most distinguishing
feature, when compared with other models, is that water vapor and liquid water
are treated separately.
ANALYTICAL STUDY
Governing Equations
The coupled system of differential equations for heat and mass transfer in a
porous material is derived from the following basic laws for mass and energy con-
servation:
Water conservation
@rL
@t|{z}
a
ˆ ¡r ¢ JL|‚‚‚‚{z‚‚‚‚}
b
¡ _EL|‚{z‚}
c
‡ @rD
@t|{z}
d
…1†
Term a represents the rate of change of liquid water rL in a unit volume of the porous
material, b is the mass of water transferred by convection, c is the mass lost by
evaporation, and d is the mass gained from dehydration of chemically bound water.
Water vapor conservation
@…eG ~rV †
@t
ˆ ¡r ¢ JV ‡ _EL …2†
Air conservation
@…eG ~rA†
@t
ˆ ¡r ¢ JA …3†
Energy conservation
…rC†
@T
@t|‚‚‚‚‚{z‚‚‚‚‚}
a
ˆ ¡r ¢ …¡krT †|‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚}
b
¡ …rCv† ¢ rT
|‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚}
c
¡ lE _EL|‚‚{z‚‚}
d
¡ lD
@rD
@t|‚‚‚{z‚‚‚}
e
…4†
Term a represents the rate of change of the accumulated energy in a unit volume, b is
the energy diffused by conduction (Fourier’s law), c is the energy transported by £uid
£ow (i.e., convection), d is the energy required for evaporation of liquid water, e is
the energy required for release of bound water by dehydration.
Fick’s law and total mass £ux expressions
JA ˆ eG ~rAvG ¡ eG ~rGDAVr
~r
A
~rG
…5†
TRANSFER IN CONCRETE SUBJECTED TO FIRE 689
JV ˆ eG ~rV vG ¡ eG ~rGDAVr
~rV
~rG
…6†
JL ˆ rLvL …7†
The £uxes are de¢ned per unit area of concrete. The mass concentrations (densities)
~ri are de¢ned per unit volume of gaseous mixture, and ri is de¢ned per unit volume
of concrete. The mass concentrations of air and water vapor per unit volume of
concrete are eG ~rA and eG ~rV , respectively. The diffusion of the adsorbed water
on the pore surface is ignored because it is assumed that this £ux is negligible.
Darcy’s law
vG ˆ ¡
KKG
mG
rPG …8a†
vL ˆ ¡
KKL
mL
rPL …8b†
It is assumed that PL ˆ PG and the in£uence of the capilary pressure is
accounted for by KL.
Ideal gas law
PA ˆ RA ~rAT …9a†
PV ˆ RV ~rV T …9b†
Partial pressures and partial densities
PG ˆ PA ‡ PV …10a†
~rG ˆ ~rA ‡ ~rV …10b†
Volume ratio of liquid water
eL ˆ
rL
rL
…11†
Porosity (the equation de¢nes the gas volume fractions eG)
por ˆ 0por ‡ p¤or ˆ eL ‡ eG …12†
The following relationship is required so that the number of unknowns is equal
to the number of available equations:
690 R. T. TENCHEV ET AL.
Sorption curves
rL ˆ rCemfr
PV
PSat
; T …13†
The sorption curves are derived by analytical reasoning and enhanced empirically to
give the best ¢t to experimental data [19^22]. Details are given in Appendix A,
Eq. (A3).
The following two experimental relationships also are used:
Release of bound water by dehydration
rD ˆ fD…T ; rCem† …14†
Increase of porosity
p¤or ˆ fp…T† …15†
The empirical function fp…T † represents the increase of porosity with temperature.
The latter is because of several factors: dehydration of bound water, chemical
decomposition of aggregates, thermal strains, and microcracking. If bound water
dehydration is only considered, then
p¤or ˆ ro=rL:
Energy transport by £uid £ow
rCv ˆ rLCLvL ‡ …~rVCV ‡ ~rACA†eGvG …16†
This term should be used when the thermal conductivity k is determined experimen-tally for absolutely dry concrete. Otherwise the effect of the convective transfer
already may be represented by the temperature dependence of k. Numerical tests
have shown that much better agreement between numerical and experimental results
is achieved if the term is neglected when k is given by [47]. So, further on, rCv ˆ 0.
The dependencies of other material properties on temperature or pressure are
given in Appendix A. Whether such dependencies exist will not in£uence the ¢nal
form of the differential equations.
In the formulation of the above equations and in the subsequent calculations
the following assumptions were made:
(i) There is thermal equilibrium between all phases within an in¢nitesimal
volume.
(ii) Water vapor, air, and their gaseous mixture behave as ideal gases.
(iii) The dependence of mass £uxes on temperature gradient (Soret effect) is
negligible according to the test data for moisture transfer [19^22].
(iv) There is no diffusion of bound water. It diffuses and evaporates only after
it is released as free water.
TRANSFER IN CONCRETE SUBJECTED TO FIRE 691
System of Differential Equations
From Eqs. (1)^(4) the following system is obtained:
@…eG ~rA†
dt
ˆ ¡r ¢ JA …17a†
@…eG ~rV †
@t
‡
@rL
@t
¡
@rD
@t
ˆ ¡r ¢ …JV ‡ JL† …17b†
…rC†
@T
@t
¡ lE
@rL
@t
‡ …lD ‡ lE†
@rD
@t
ˆ r…krT † ‡ lEr ¢ JL ¡ …rCv†rT …17c†
Equation (19b) results from summing up Eqs. (1)-(2). Equation (19c) is
obtained from Eq. (4) after eliminating _EL by means of Eq. (1).
The set T , PG, ~rV is chosen as the primary unknown variables. As a check for
the analytical derivation and software programming all analyses were repeated
for sets T , PA, PV and T , ~rA, ~rV . The results were almost identical.
The system of transient differential equations, after some algebraic
manipulations given in Appendix B, may be written in a form suitable for a ¢nite
element solution:
CTT
@T
@t
‡ CTP
@PG
@t
‡ CTV
@ ~rV
@t
ˆ r ¢ …KTTrT ‡ KTPrPG ‡ KTVr ~rV †
CAT
@T
@t
‡ CAP
@PG
@t
‡ CAV
@ ~rV
@t
ˆ r ¢ …KATrT ‡ KAPrPG ‡ KAVr~rV †
CMT
@T
@t
‡ CMP
@PG
@t
‡ CMV
@ ~rv
@t
ˆ r ¢ …KMTrT ‡ KMPrPG ‡ KMVr ~rV † …18†
where
CTT ˆ …rC† ‡ …lD ‡ lE†
@rD
@t
¡ lE
@rL
@T
CTP ˆ 0 CTV ˆ ¡lE
@rL
@ ~rV
…19a†
CAT ˆ ~rA
@ror
@T
¡
~rA
rL
@rL
@T
¡ eGPG
RAT 2
CAP ˆ
eG
RAT
CAV ˆ ¡eG
RV
RA
¡
~rA
rL
@rL
@ ~rV
…19b†
CMT ˆ 1 ¡
rV
rL
@rL
@T
‡ ~rA
@por
@T
¡ @rD
@T
CMP ˆ 0 CMV ˆ eG ‡ 1 ¡
~rV
rL
@rL
@ ~rV
…19c†
KTT ˆ keff ; KTP ˆ ¡l
rLKKL
mL
; KTV ˆ 0 …20a†
KAT ˆ ¡
DVAeG ~rV
~rGRAT
PG
T
; KAP ˆ
KKG
mG
eG ~rA ‡
DVAeG ~rV
~rGRAT
KAV ˆ ¡ ~rA ‡
RV
RA
rV
DVAeG
~rG
…20b†
KMT ˆ ¡KAT ; KMP ˆ
KKG
mG
eG ~rV ¡
DVAeG ~rV
~rGRAT
‡
~rLKKL
mL
; KMV ˆ ¡KAV …20c†
The sequence of computations of all coef¢cients is given in Appendix A.
692 R. T. TENCHEV ET AL.
Boundary Conditions
Gaseous mixture mass balance
JG ¢ n|‚‚{z‚‚}
a
¡ b… ~rG ¡ ~rG;1†|‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚}
b
ˆ 0 …21†
Term a is the £ux transferred to the boundary surface from within the body in the
direction of the normal and b is the substance dissipated into the surrounding
medium. This equation is used only to simplify the energy balance, Eq. (24).
Energy balance
k
@T
@n|‚{z‚}
a
¡ …HG ¡0 HG†JG ¢ n|‚‚‚‚‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚‚‚‚‚}
b
‡ …HG ¡0 HG†b…~rG ¡ ~rG;1†|‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚}
c
‡ lJL ¢ n|‚‚‚{z‚‚‚}
d
‡ hqr…T ¡ T1†|‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚}
e
ˆ 0
…22†
Term a is the heat energy transferred to the boundary surface from inside the body by
conduction. Term b is the heat energy accumulated in the gaseous mixture, which
enters the boundary surface from the inside of the concrete. Term c is the heat energy
accumulated in the gaseous mixture, that leaves from the boundary surface. Term d
is the heat energy used for evaporation of liquid water at the boundary. Term e
is the heat energy dissipated by convection and radiation to the surrounding
medium. The enthalpy of the gaseous mixture at ambient and current conditions
is 0HG and HG [J/kg], respectively.
The exposed surface of the column will be dry in a very short time after the start
of the ¢re (if not already dried during its service life), and the liquid £ux thus can be
omitted, that is,. JL ˆ 0. The ¢nal expression for the boundary condition reduces
to the one used in conventional thermal analysis:
k
@T
@n
‡ hqr…T ¡ T1† ˆ 0 …23†
Water vapor boundary condition
JV ¢ n ¡ b…~rV ¡ ~rV ;1† ˆ 0 …24†
The expression for the vapor £ux JV in terms of the prime unknowns is given in
Appendix B, Eqs. (B1e)^(B1f ). The gaseous pressures on the concrete surface and in
the surrounding environment are both equal to the atmospheric pressure, so
@PG=@n ˆ 0. The simultaneous solution of Eqs. (22)^(24) gives:
@T
@n
@ ~rV
@n
8
>><
>>:
9
>>=
>>;
ˆ KTT 0
KVT KVV
¡1 hqrT1
b ~rV;1|‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚}
a
¡ KTT 0
KVT KVV
¡1 hqr 0
0 b
T
~rV|‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚}
b
…25†
Term a contributes to the load vector; term b contributes to the stiffness matrix.
TRANSFER IN CONCRETE SUBJECTED TO FIRE 693
Pressure boundary condition
The pressure on the boundary surface is equal to the atmospheric pressure:
PG ˆ PG;1 …26†
FINITE ELEM ENT FORMULATION
The system of differential equations to be solved is
C
@ U
@t
¡ r ¢ …KrU † ‡ KVrU ˆ 0 …27†
with boundary conditions
@ U
@n
ˆ R1 ¡ Rk U ; PG ˆ PG;1 …28†
and initial conditions at t ˆ 0
U ˆ 0 U …29†
where the matrix notation is given in Appendix C.
The ¢nite element solution is based on a piecewise approximation of the
unknown vector function U ˆ
P
U e ˆ
P
Nue, where N is a matrix of polynomial
shape functions and the superscript refers to element e. The approximate solution
must minimize the integral of the weighted residual error. The minimization is
applied for the whole domain, but it is valid for each element Oe, too.
Xnel
eˆ1
Z
Oe
NT CeN
@ue
@t
¡ r…KerN†ue ‡ KeVrNu
e dO ˆ 0 …30†
with natural boundary conditions
@N
@n
ue ˆ Re1 ¡ R
e
kNu
e …31†
where nel is the total number of elements. Integration by parts results in a system of
¢rst-order differential equations:
Ĉ
@u
@t
‡ K̂u ˆ R̂ …32†
where
Ĉ ˆ
Xnel
eˆ1
Ĉe; K̂ ˆ
Xnel
eˆ1
…K̂e ‡ K̂eV ‡ K̂eR†; R̂ ˆ
Xnel
eˆ1
R̂e; u ˆ
Xnel
eˆ1
ue …33†
Ĉe ˆ
Z
Oe
NT CeN dO; K̂e ˆ
Z
Oe
rNT KerN dO; K̂eV ˆ
Z
Oe
NT KeVrN dO …34†
K̂eR ˆ
Z
Ge
NT KeRekN dG R̂
e ˆ
Z
Ge
NT KeRe1 dG …35†
694 R. T. TENCHEV ET AL.
and Ge is a part of the boundary surface where convection boundary conditions are
applied.
A ¢nite difference scheme is used for the time discretisation of Eq. (34), which
yields
K £ …nu† ˆn R …36†
where
K ˆ yK̂ ‡
1
Dt
Ĉ …37†
nR ˆ Ĉ
Dt
¡ …1 ¡ y†K̂
" #
£ …n¡1u† ‡ …1 ¡ y† £ …n¡1R̂† ‡ y £ …nR̂† …38†
and n is the number of the current time step, Dt is the magnitude of the time step, and
y is a prede¢ned coef¢cient (0, 1/2, 2/3, or 1 rendering the explicit, Crank^Nicolson,
Galerkin, and backward-difference methods, respectively).
NUM ERICAL EXAM PLE
The test problem is shown on Figure 1 and represents a rectangular cross sec-
tion ABCD of a concrete column exposed to ¢re. The problem is two dimensional
and so is the ¢nite element program developed by the authors. The task is simpli¢ed
and it is to compute the temperature, pore pressure, and water vapor content dis-
tribution along the horizontal axis of symmetry x, thus making the problem one
dimensional. The corresponding boundary conditions are that sides DA and CB
are exposed to convection heat, vapor transfer, and atmospheric pressure; sides
AB and CD are either isolated or it may be assumed that during the time period
under consideration the convection on these sides will not affect the results for
the x axis. Taking into account that y is an axis of symmetry, the domain to be
analyzed is the rectangle KLMN. The mesh is composed of one element in y direction
and a variable number of elements in the x direction. Side NK, which is part of side
DA, has the boundary conditions applied on side DA. The other sides are isolated,
that is,no boundary conditions in the FE analysis are speci¢ed. Side KL is 0.2 m
long. Side NK may have arbitrary length. In this study it is determined from
the general ¢nite element method recommendation that the generated elements
are square shaped. Numerical values for the rest of the input data and a £owchart
of the calculation of all variables are given in Appendix A. Results are presented
in Figures 2^6.
Four noded quadrilateral elements and the backward difference algorithm,
y ˆ 1 in Eqs. (37)^(38), are used. The in£uence of mesh density and time step
magnitude on the results is investigated. Meshes of 50, 100, 200, and 400 elements
and time steps varying from 30 sec to 0.2 sec are used. Temperature distribution
is hardly affected by mesh density or time step. Pore pressure and water vapor con-
tent show time oscillations for the mesh of 50 elements, which are, on average, under-
estimated. The oscillations may be reduced by increasing the mesh density or
decreasing the time step or both. Based on this comparative study, a mesh of
TRANSFER IN CONCRETE SUBJECTED TO FIRE 695
200 elements and a time step of 2 sec are chosen for the subsequent analyses. The
convection term, Eq. (16), when taken into account, does not cause spurious
oscillations and no special numerical treatment is required. The same observation
is made in [36].
The numerical results for the temperature distributions as a function of the
distance from the ¢re exposed surface, x, at times t ˆ 30 min and t ˆ 60 min and
the temperature distributions as function of time at a distance x ˆ 13 mm and
x ˆ 32 mm are shown in Figures 2a and 2b, respectively, by solid lines. The discrete
points correspond to test results, [51, 49]. The agreement is very good.
The numerical results for the pore pressure, PG, the mass content of water
vapor per unit volume of gas, ~rV , and per unit volume of concrete, eG ~rV , and
the mass content of liquid water, rL, are shown in Figure 3. Graphs t10 correspond
to time of t ˆ 10 min of exposure to ¢re; graphs x10 give the results at a distance
of x ˆ 10 min from the ¢re exposed surface. The same labeling pattern is applied
to the rest of the graphs.
The most important result is that high values for the pore pressure are com-
puted. It is most likely that, in combination with the thermal stresses, these will
cause spalling. Reported numerical values for pore pressure in concrete under
intense heating vary in wide margins: from less than 1 MPa, [14, 24]; about 3.5 MPa,
[48, 49], and up to 16 MPa [53].
Comments on the interaction between the plotted variables are given for
graphs t60 and x30. They are valid for the other graphs, too.
The zone of increased pore pressure (see Figure 3a), corresponds to the
increased content of water vapor, rV (see Figure 3b). It is because of the higher
rate of evaporation than the rate of mass transfer. The sudden drop in eG ~rV , (vapor
content per unit volume of concrete, see Figure 3c) is because of the drop in the
gas volume fraction, eG, and corresponds to the transition between wet and dry zone
(see the dotted line a in Figure 3d).
The distribution of liquid water (see Figure 3d), shows several well-de¢ned
zones: a dry zone between x ˆ 0 and line a; a zone of evaporation between lines
a and b; a zone of increased water content between lines b and c; and a zone of
Figure 1. Layout of test probelm ABCD is the cross section of a concrete column; KLMN is part meshed
and analysed.
696 R. T. TENCHEV ET AL.
initial state rL ˆ 0rL ˆ 60 kg=m3 between line c and z ˆ 0:1 m. The increased water
content is because of the condensation of the water vapor when it is transferred
into the cooler interior of the concrete. It corresponds to the so-called moisture clog
[5, 6]. A similar numerical result is presented in [34, 35, 49]. The ability to predict
the moisture clog is an advantage of the model when compared with other models
that fail to do so [24].
A rapid increase of PG as function of time at a given location (graph x30 in
Figure 3e) starts with the formation of the moisture clog (dotted line d going to
Figure 3h). The moisture clog blocks the vapour mass transfer toward the interior.
During the time interval between lines d and e the vapor content ~rV increases, which
causes an increase of PG, too. The pores are almost saturated with liquid water (see
Figure 3h). The gas volume fraction eG is very low, which results in low vapor content
per unit volume of concrete, eG ~rV , Figure 3g. The time when the maximum satu-
ration at x ˆ 30 mm occurs is 30 minutes and it corresponds to temperature about
200¯C (see Figure 2b). After that, reduction of rL follows because of evaporation.
The increase in eG ~rV is because of the increase of eG. After its maximum value,
PG decreases only slightly. After evaporation of all liquid water, the relatively high
value of PG is maintained by the thermal expansion of the gaseous mixture, as well
as by the vapour transferred from the interior.
The mass content of liquid water rL is computed from the sorption curves,
Eq. (A3), and not from the water conservation condition, Eq. (1), because the rate
of evaporation is not known. The computed domain and time distributions of
~rV are continuous (some oscillation may be observed for coarse mesh and large
time steps), while the sorption curves have a sharp discontinuity in the interval
0:96 µ PV =PSat µ 1:04. It represents the intense evaporation, and when this happens
(a ¢xed time) there is a discrepancy between the mass of water lost and the mass of
vapor gained. However, over a time interval, the equations of mass conservation,
Eqs. (1)^(3), will insure the required mass balance.
Figure 2. Distribution of temperature T [¯C]. Numerical results (solid line). (a) ‡ , £ test results, Ref.
[51]; (b) &, * test results, Ref. [49].
TRANSFER IN CONCRETE SUBJECTED TO FIRE 697
Figure 3. Distribution of pore pressure, PG [MPa]: liquid water content, rL [kg/m3]; vapor content per
unit volume of gas, eG ~rV [kg/m
3]; and per unit volume of concrete, ~rV [kg/m
3], at times
t10 ˆ 10 min, t30 ˆ 30 min, t60 ˆ 60 min; at locations x10 ˆ 0:01 m, x20 ˆ 0:02 m, x30 ˆ 0:03 m.
698 R. T. TENCHEV ET AL.
The maximum pore pressure, max PG, its location x, and the temperature
where it occurs are plotted as a function of time in Figure 4. There is a rapid increase
of max PG up to about 30 minutes after the beginning of the ¢re. At that time its
location is at x ˆ 0:02 m from the ¢re exposed surface and the temperature is about
220¯C. If spalling is caused by high pore pressure, this is the most likely time
and place where it will occur. This is in good agreement with experimental obser-
vations [41, 43]. After that time, max PG and the temperature where it occurs remain
almost constant. Its location moves steadily toward the interior where evaporable
water is still available.
The ratios between the current and initial contents of liquid water, water
vapour, and air are shown in Figure 5. As it may be expected, liquid water decreases
when evaporation takes place, and the total moisture (water plus vapor) decreases as
vapor escapes from the surface. Air decreases because it is displaced by vapour and
expands thermally. The total change of the moisture content may be measured exper-
imentally and thus may be used to check the accuracy of the presented analysis indi-
rectly.
The in£uence of permeability K , initial water content 0rL, and porosity por on
pore pressure PG is shown in Figure 6. The parameters are varied independently,
although there is an interconnection between them [39]. According to experimental
data [50, 51], permeability may vary by several orders of magnitude. Within this
variation permeability has the most pronounced in£uence on pore pressure, when
compared with other material properties.
Initial water content 0rL has a ¢xed range of variationöfrom zero (absolutely
dry concrete) to the value of porosity (fully saturated concrete). The higher the water
content, the higher the pore pressure, becausemore vapor can be produced. The
relationship is almost linear, and the total pore pressure increase is about 1.5 to
2 times. In absolutely dry concrete, pore pressure results from the evaporation
of bound water that has been released by dehydration as free water.
Lower porosity increases pore pressure because there is less room for the evap-
orated water to expand. The variation of porosity has a much more pronounced
Figure 4. Maximum pore pressure, max PG [MPa], its location, x [m], from the ¢re exposed surface and
the temperature T [¯C] at x as function of time.
TRANSFER IN CONCRETE SUBJECTED TO FIRE 699
in£uence when permeability is low; see graphs n and p in Figure 6. The worst com-
bination is low permeability, low porosity, and high water content.
CONCLUSIONS
A system of transient differential equations for coupled heat, mass transfer,
and pore pressure buildup has been derived based on well-established physical laws.
Water vapor content and liquid water content are treated separately. Several physi-
cal quantities are determined based on published experimental data and
Figure 5. Relative content of liquid water, rL=
0rL; water vapor, ~rV =
0 ~rV , and air, ~rA=
0 ~rA, after 60 minutes
of exposure to ¢re.
Figure 6. In£uence of permeability, K [m2], initial water content, 0rL [kg/m3], and porosity, por, on maxi-
mum pore pressure PG [MPa] after 20 minutes of ¢re exposure. a: 0rL ˆ 30, por ˆ 0:16; b: 0rL ˆ 60,
por ˆ 0:16; c: 0rL ˆ 30, por ˆ 0:08; d: 0rL ˆ 60, por ˆ 0:08; e: K ˆ 5 £ 10¡16 , por ˆ 0:16; f:
K ˆ 5 £ 10¡16 , por ˆ 0:08; g: K ˆ 5 £ 10¡17, por ˆ 0:16; h: K ˆ 5 £ 10¡17 , por ˆ 0:08; k: 0rL ˆ 30,
K ˆ 5 £ 10¡16 ; m: 0rL ˆ 60, K ˆ 5 £ 10¡16 ; n: 0rL ˆ 30, K ˆ 5 £ 10¡17 ; p: 0rL ˆ 60, K ˆ 5 £ 10¡17 .
700 R. T. TENCHEV ET AL.
semianalytical relationships. A detailed ¢nite element formulation and a correspond-
ing £owchart of computational procedures are presented. A two-dimensional ¢nite
element method program is developed and a numerical example (a cross section
of a concrete column under ¢re) has been presented. The computed results show
that signi¢cant pore pressure builds up in the interior, close to the ¢re exposed
surface, which agrees with published experimental data. Low permeability and high
initial water content increase the magnitude of the pore pressure. Future work will
investigate the joined effect of high pore pressure, thermal and mechanical stresses,
and the prediction of spalling in concrete exposed to ¢re.
REFERENCES
1. W. J. Copier, The Spalling of Normal Weight and Lightweight Concrete Exposed to Fire,
Fire Safety of Concrete Structures, ACI Publ. SP^80, American Concrete Institute (ACI),
Detroit, MI, pp 219^237, 1983.
2. H. L. Mathotra, Spalling of Concrete in Fires, Technical notes 118, Construction Industry
Research and Information Assn., London, 1994.
3. R. J. Connolly, The Spalling of Concrete, Fire Engineering Journal, pp. 38^40, Jan., 1997
4. F. J. Ulm, O. Coussy, and Z. P. Bazant, The ``Chunnel’’ Fire. I: Chemoplastic Softening in
Rapidly Heated Concrete, J. Engineering Mechanics, ASCE, vol. 125, no. 3, pp. 272^282,
1999.
5. F. J. Ulm, O. Coussy, and Z. P. Bazant, The ``Chunnel’’ Fire. II: Analyses of Concrete
Damage, J. Engineering Mechanics, ASCE, vol. 125, no. 3, pp. 283^289 , 1999.
6. T. Z. Harmathy, Effect of Moisture on the Fire Endurance of Building Materials, ASTM
No. 385, Philadelphia, pp. 74^95, 1965.
7. A. K. Abdel-Rahman and G. H. Ahmed, Computational Heat and Mass Transport in
Concrete Walls Exposed to Fire, Numer. Heat Transfer Part A, vol. 29, pp. 373^395, 1996.
8. Y. Anderberg, Spalling Phenomena of HPC and OC, in L. T .Phan, N. J. Carino, D.
Duthinh, and E. Garboczi (eds.), Proc., Int. Workshop on Fire Performance of
High-Strength Concrete, NIST Spec. Publ. 919, National Institute of Standards and
Technology, Gaithersburg, MD, pp. 69^73, 1997.
9. Z .P. Bazant, Analysis of Pore Pressure, Thermal Stresses and Fracture in Rapidly Heated
Concrete, in L. T. Phan, N. J. Carino, D. Duthinh, and E. Garboczi (eds.), Proc., Int.
Workshop on Fire Performance of High-Strength Concrete, NIST Spec. Publ. 919,
National Institute of Standards and Technology, Gaithersburg, MD pp. 155^164, 1997.
10. J. A. Purkiss and K. N. Mustapha, A Study of the Effect of Spalling on the Failure Modes
of Concrete Columns in a Fire, Concrete 95: Towards Better Concrete Structures,
Brisbane, Australia, pp. 263^272 , 1995.
11. J. A. Purkiss, W. A. Morris, and R. J. Connolly, Fire Resistance of Reinforced
ConcreteöCorrelation of Analytical Methods with Observed Experimental
Behaviour, in C. J. Franks and S. Grayson (eds.), Inter£am ’96, Interscience
Communications Ltd., London, pp. 531^541 , 1996.
12. A. V. Luikov, Heat and Mass Transfer in Capillary Porous Bodies, Pergamon, Oxford,
1966.
13. A. V. Luikov, Systems of Differential Equations of Heat and Mass Transfer in
Capillary-Porous Bodies (Review), Int. J. Heat Mass Transfer, vol. 18, pp. 1^14, 1975.
14. A. Dayan and E. L. Gluekler, Heat and Mass Transfer within an Intensely Heated
Concrete Slab, Int. J. Heat Mass Transfer, vol. 25, pp. 1461^1467, 1982.
15. P. Chen and D. C. T. Pei, A Mathematical Model of Drying Processes, Int. J. Heat Mass
Transfer, vol. 32, pp. 297^310 , 1989.
TRANSFER IN CONCRETE SUBJECTED TO FIRE 701
16. J. R. Philips and D. A. deVries, Moisture Movement in Porous Materials under
Temperature Gadients, Trans. Am. Geophys. Union, vol. 38, pp. 222^232, 1957.
17. D. A. deVries, Simultaneous Transfer of Heat and Moisture in Porous Materials, Trans.
Am. Geophys. Union, vol. 39, pp. 909^916, 1958.
18. S. Whitaker, Simultaneous Heat, Mass and Momentum Transfer in Porous Media: A
Theory of Drying, in J. P. Hartnett and T. F. Irvine (eds.), Advances in Heat
Transfer, Academic Press, New York, vol. 13, pp. 119^203, 1977.
19. Z. P. Bazant and W. Thonguthai, Pore Pressure and Drying of Concrete at High
Temperature, J. Eng. Mecnanics Division. ASCE, vol. 104, pp. 1059^1079 , 1978.
20. Z. P. Bazant and W. Thonguthai, Pore Pressure in Heated Concrete Walls: Theoretical
Prediction, Magazine of Concrete Research, vol. 31, no. 107, pp. 67^76, 1979.
21. Z. P. Bazant, J. C. Chern, and W. Thonguthai, Finite Element Program for Moisture and
Heat Transfer in Heated Concrete, Nuclear Engineering and Design, vol. 68, pp. 61^70,
1981
22. Z. P. Bazant and M. F. Kaplan, Concrete at High Temperatures, Longman, London, 1996
23. A. K. Abdel-Rahman and G. N. Ahmed, Computational Heat and Mass Transport in
Concrete Walls Exposed to Fire, Numer. Heat Transfer Part A, vol. 29, pp. 373^395, 1996.
24. G. N. Ahmed and J. P. Hurst, Coupled Heat and Mass Transport Phenomena in Siliceous
Aggregate Concrete Slabs Subjected to Fire, Fire and Materials, vol. 21, pp. 161^168,
1997.
25. J. Selih, A. C. M. Sousa and T. W. Bremmer, Moisture and Heat Flow in Concrete Walls
Exposed to Fire, J. of Engineering Mechanics, ASCE, vol. 120, pp. 2028^2043 , 1994.
26. M. D. Mikhailov and B. K. Shishedjiev, Temperature and Moisture Distributions during
Contact Drying of a Moist Porous Sheet, Int. J. Heat Mass Transfer, vol. 18, pp.
15^24, 1975.
27. M. D. Mikhailov, Exact Solution of Temperature and Moisture Distributions in a Porous
Half-Space with Moving Evaporation Front, Int. J. Heat Mass Transfer, vol. 18, pp.
797^804 , 1975.
28. R. W. Lewis, K. Morgan, and H. R. Thomas, The Nonlinear Modeling of Drying-Induced
Stresses in Porous Bodies, in A. S. Mujumdar (ed.), Advances in Drying, McGraw-Hill,
New York, 1983.
29. G.Comini and R. W. Lewis, A Numerical Solution of Two-Dimensional Problems
Involving Heat and Mass Transfer, Int. J. Heat Mass Transfer, vol. 19, pp.
1387^1392 , 1976.
30. H. R. Thomas, K. Morgan, and R. W. Lewis, A Fully Nonlinear Analysis of Heat and
Mass Transfer Problems in Porous Bodies’, Int. J. Numerical Methods in
Engineering , vol. 15, pp 1381^1393 , 1980.
31. R. W. Lewis, K. Morgan, R. H. Thomas, and K. N. Seetharamu, The Finite Element
Method in Heat Transfer Analysis, John Wiley & Sons,Chichester, 1996.
32. J. Irudayaraj and Y. Wu, Finite Element Analysis of Coupled Heat, Mass, and Pressure
Transfer in Porous Biomaterials, Numer. Heat Transfer Part A, vol. 26, pp. 337^350,
1994.
33. Y. Zhou, R.K.N.d Rajapakse, and J. Graham, Coupled Heat-Moisture-Air Ttransfer in
Deformable Unsaturated Media, J. Engineering Mechanics, vol. 124, pp. 1090^1099 ,
1998.
34. P. Majumdar and A. Marchertas, Heat Moisture Transport and Induced Stresses in
Porous Materials under Rapid Heating, Numer. Heat Transfer Part A, vol. 32, pp.
111^130 , 1997.
35. P. Majumdar, A. Gupta, and A. Marchertas, Moisture Propagation and Resulting Stress
in Heated Concrete Walls, Nuclear Engineering and Design, vol. 156, pp. 159^165, 1995.
702 R. T. TENCHEV ET AL.
36. Z. X. Gong and A. S. Mujumdar, A Two Dimensional Finite Element Model for
Kiln-Drying of Refractory Concrete, Drying Technology, vol. 13, pp. 585^605 , 1995.
37. Z. X. Gong, B. Song, and A. S. Mujumdar, Numerical Simulation of Drying of Refractory
Concrete, Drying Technology, vol. 9, pp. 479^500, 1991.
38. Z. X. Gong and A. S. Mujumdar, The In£uence of an Impermeable Surface on Pore Steam
Pressure during Drying of Refractory Concrete Slabs, Int. J. Heat Mass Transfer, vol. 38,
pp. 1297^1303 , 1995.
39. A. E. Scheidegger, The Physics of Flow through Porous Media, 3rd ed., University of
Toronto Press, 1974.
40. S. B. Nasrallah and P. Perre, Detailed Study of a Model of Heat and Mass Transfer during
Convective Drying of Porous Mmedia, Int. J. Heat Mass Transfer, vol. 31, pp. 957^967,
1988.
41. J. A. Purkiss, Fire Safety Engineering Design of Structures, Butterworth-Heinemann,
Oxford, 1996.
42. Y. A. Cengel, Heat Transfer: A Practical Approach, McGraw-Hill, New York, 1998.
43. Ir. W. J. Copier, The Spalling of Normal Weight and Lightweight Concrete on Exposure
to Fire, Heron, vol. 24, no.2, 1979.
44 X. Lin, W. Sun, and S. Y. N. Chan, Effect of Heating and Cooling Regimes on Residual
Strength and Microstructure of Normal and High-Performance Concrete, Cement and
Concrete Research, vol. 30, pp. 379^383, 2000.
45. F. Vodak, C. J. Drchalova, S. Hoskova, O. Kapickova, O. Michalko, P. Semerak, and J.
Toman, Thermophysical Properties of Concrete for Nuclear-Safety Related
Structures, Cement and Concrete Research, vol. 37, pp.415^426, 1997.
46. A. Atkison and A. K. Nickerson, The Diffusion of Ions through Water-Saturated
Cement, J Materials Science, vol. 19, pp. 3068^3078 , 1984.
47. PrEN 1992-1-2 (1st draft), Eurocode 2: Design of Concrete StructuresöPart 1.2: General
Rules Structural ¢re design, Comitë Europëan de Normalization (British Standard
Institution), London, pp. 66^67, Oct. 2000.
48. G. R. Consolazio, M. C. McVay, and J. W. Rish III, Measurment and Prediction of Pore
Pressure in Cement Mortar Subjected to Elevated Temperature, Int Workshop on Fire
Performance of High-Strength Concrete, NIST Gaithersburg, MD, Feb. 13^14, paper
B.8, pp 125^148, 1997.
49. G. N. Ahmed and J. P. Hurst, An Analytical Approach for Investigating the Cases of
Spalling of High-Strength Concrete at Elevated Temperatures, Int Workshop on Fire
Performance of High-Strength Concrete, NIST Gaithersburg, MD, Feb. 13^14, paper
B.6, pp 95^108, 1997.
50. G. L. England and N. Khoylou, Moisture Flow in Concrete under Steady State
Non-Uniform Temperature States: Experimental Observations and Theoretical
Modelling, Nuclear Engineering and Design, vol. 156, pp. 83^107, 1995.
51. A. S. El-Dieb and R. D. Hooton, Water-Permeability Measurement of High Performance
Concrete Using a High-Pressure Triaxial Cell, Cement and Concrete Research, vol. 25,
pp.1199^1208 , 1995.
52. Anonymous, Design and Detailing of Concrete Structures for Fire Resistance, Interim
Guidance by a Joint Committee of The Institution of Structural Engineersöthe
Concrete Society, April, 1978.
53. J. W. Ju and Y. Zhang, Axisymmetric Thermomechanical Constitutive and Damage
Modeling for Air¢eld Concrete Pavement under Transient High Temperature,
Mechanics of Materials, vol. 29, pp. 307^323, 1998.
TRANSFER IN CONCRETE SUBJECTED TO FIRE 703
APPENDIX A (FLOW CHART)
The computations are done according to the following input data and
£owchart.
STEP 0: Constants and initial values:
High strength concrete:
rCem ˆ 300 ‰kg=m3Š; rS ˆ 2400 ‰kg=m3Š; 0 ror ˆ 0:08
0K ˆ 5:0 £ 10¡17 ‰m2Š; equivalent to permeability 3 £ 10¡13 ‰m=sŠ; see Note:
hq ˆ 25 ‰W=…m2 ¯K Š; Ref: ‰41Š
Initial conditions:
0T ˆ 20 ‰¯CŠ; 0 PG ˆ PG;1 ˆ 0:1 MPa; 0 rL ˆ 60 ‰kg=m
3Š; j ˆ 0:8
0PV ˆ0 PSat…0T† 0 ~rV ˆ0 PV =…0TRV † 0rV;1 ˆ0 ~rV j
Thermal load:
T1 ˆ0 T ‡ 345 log10…8t ‡ 1†; t¡time [min]; (standard fire curve, ISO 834).
Constant material properties, Ref. [42]:
RA ˆ 287 ‰J=…kg¯K †Š; RV ˆ 461:5 ‰J=…kg¯K †Š; rL ˆ 1000 ‰kg=m
3Š
STEP 1: Computation of variables at time step n:
The values of n¡1T , n¡1PG, and n¡1 ~rV are known at time step n ¡ 1. They will be
used to compute the coef¢cients required for the computations at time step i. Further
on the left superscript n ¡ 1 will be omitted. Temperature is ¯C except where
explicitly stated otherwise.
PV ˆ RV ~rV T ; …T in ¯K † …A1†
Psat ˆ Psat…T† ¡ tabulated data ‰42Š …A2†
Sorption curves [19^22]:
rL
PV
PSat
ˆ
rCem
PV
PSat
0rL
rCem
1=m…T†
if
PV
PSat
µ 0:96
0rL 1 ‡ 0:12
PV
PSat
¡ 1:04 if PV
PSat
¶ 1:04
rL…0:96† ‡
PV
PSat
¡ 0:96 rL…1:04† ¡ rL…0:96†
0:08
if 0:96 <
PV
PSat
< 1:04
8
>>>>>>>>>><
>>>>>>>>>>:
…A3†
704 R. T. TENCHEV ET AL.
where
m…T † ˆ 1:04 ¡
…T ‡ 10†2
…T ‡ 10†2 ‡ 22:3…25 ‡ 10†2
In Refs. [19^22] water vapor and liquid water content are treated as one variable
(moisture). The case PV =PSat ¶ 1:04 refers to the condition of saturated concrete
when moisture occupies all the available pore space; hence the formulation in [19^22]
makes use of the porosity. In this article the pores may be partially ¢lled with water,
so the initial water content 0rL is used instead. A smoothed transition at the end
points of the interval 0:96 µ PV =PSat µ 1:04 is used, [35], which makes the slope
of the sorption curves continuous.
@rL
@T
ˆ rL…T ‡ DT † ¡ rL…T †
DT
…A4†
@rL
@ ~rV
ˆ rL…
~rV ‡ D ~rV † ¡ rL… ~rV †
D ~rV
…A5†
Eqs. (A1)^(A2) are included in the computation of @rL=@T , too.
DT ˆ 0:0001T ; D ~rV ˆ 0:0001 ~rV
rD ˆ rcem £
0 if T µ 200¯C
7:0 £ 10¡4…T ¡ 200† if 200¯C < T µ 300¯C
0:4 £ 10¡4…T ¡ 300† ‡ 0:07 if 300¯C < T µ 800¯C
0:09 if T > 800¯C
;
8
>><
>>:
…A6†
Eq. (A6) is approximation of a graph in Ref. [22].
por ˆ0 por £
1 if T < 100¯C
aT 3 ‡ bT 2 ‡ cT ‡ d if 100¯C µ T µ 800¯C
3 if T > 800¯C
8
<
: …A7†
A threefold increase of porosity is reported [44] when heating up to 800¯C. It is
assumed here that the increase is by a cubic function of temperature, such as
por, and its ¢rst derivative are continuous functions.
K ˆ
por
0por
2=3
£0 K …A8†
It is also assumed that porosity increases equally in all dimensions and permeability
is proportional to the increase of the cross section of the pores.
eL ˆ
rL
rL
…A9†
eG ˆ por ¡ eL …A10†
PA ˆ PG ¡ PV …A11†
TRANSFER IN CONCRETE SUBJECTED TO FIRE 705
~rA ˆ
PA
RAT
…T in ¯K † …A12†
~rG ˆ ~rA ‡ ~rV …A13†
DAV ˆ D
d
t2
Ref ‰46Š D ˆ 1:89 £
…273:15 ‡ T†2:072
PG
£ 10¡5 Ref : ‰42Š …A14†
D is the coef¢cient for atmospheric diffusion. The diffusion coef¢cient in the porous
concrete, DAV , is less than D for two reasons: (i) The diffusion paths are torturous
and not parallel to the concentraton gradient; (ii) the paths have varying cross sec-
tion and may become very constricted at certain points, where the effect of he walls
on the free molecular motion is not negligible. Approximate values for the toruosity
and constrictivity are, respectively, t ˆ 3, d ˆ 0:5. In Ref. [45] DAV & 5 £ 10¡17
[m2/s] is reported.
mV ˆ mV …T † mL ˆ mL…T † mA ˆ mA…T †
Cv ˆ CV …T† CL ˆ CL…T† CA ˆ CA…T†
9
=
; …A15†
according to tabulated data, Ref. [42]
mG ˆ
~rAmA ‡ ~rV mV
~rA ‡ ~rV
or mG ˆ 0 if ~rA ˆ ~rV ˆ 0 …A16†
Asumption: weighted average.
CS ˆ 900 ‡ 80
T
120
¡ 4 T
120
2
; Ref: ‰47Š …A17†
k ˆ 2:0 ¡ 0:24 T
120
‡ 0:012 T
1202
; Ref: ‰47Š …A18†
KG ˆ 1 ¡ S;
where saturation lof liquid water is S ˆ eL=por , Refs. [14, 39].
KL ˆ 0:01; …A20†
It is reported in [39] that the surface tension of the liquid water in the capillaries of a
porous media causes about 100 times a constant reduction of PL compared with PG.
This effect is accounted for by KL. An alternative treatment of capillary pressure is
given in Ref. [39^40].
rC ˆ rSCS ‡ rLCL ‡ eG ~rV CV ‡ eG ~rACA …A21†
lE ˆ lE…T†; tabulated data, Ref. [42] …A22†
lD ˆ 2400 £ 103 Ref : ‰22Š …A23†
hr ˆ es‰…T ‡ 273:15†2 ‡ …T1 ‡ 273:15†2Š £ ‰…T ‡ 273:15† ‡ …T1 ‡ 273:15†Š …A24†
706 R. T. TENCHEV ET AL.
e ˆ 0:6 emisivity of the surface; s ˆ 5:67 £ 10¡8 ‰W=…m2 ¯K 4†Š, Stefan^Boltzman
constant.
b ˆ
hq
…rC†Air
DAV …T1†
aAir
2=3
Chilton¡Colburn analogy, Ref. [42] …A25†
where …rC†Air and aAir are heat capacity and thermal diffusivity of air, respectively.
Computation of coefficients, Eqs. (21)¡…22† …A26†
STEP 2: Solve the system of differential equations : Obtain new values: iT , iPG, i ~rV .
STEP 3: Equilibrium iterations:
STEP 1 and STEP 2 are repeated for the same time step, but coef¢cients are
calculated on the basis of the last computed values for T , PG, ~rV . The equilibrium
iteratons are stopped if either convergence is achieved or the number of allowed
iterations, Niter, is reached (Niter ˆ 5 is used). The error for time step i, equilibrium
iteration j, and the convergence criteria are de¢ned as
DijF ˆ
Pnodes
kˆ1 …
i
j¡1Fk ¡ij Fk†
2
Pnodes
kˆ1 …
i
jFk†
2 < Err; …A27†
where F stands for each one of T , PG, ~rV and Err ˆ 10¡9.
STEP 4: New time step:
Time step increment Dt ˆ 2 sec. STEP 1 to STEP 3 are repeated until the
required total time for the ¢re simulation is reached.
Note. Moisture £ux (liquid water and water vapor considered as one
entityömoisture) often is de¢ned by J ˆ ¡…a=g†rP, where a is permeability in
m/s and g ˆ 9:81 m=s2 is the gravity acceleration used for dimensional con-
venience [19^22, 34^38] . When the liquid water mass £ux is given by
Eqs. (7)^(8b), that is, JL ˆ ¡rL…KKL†=…mL†rPL, the relationship between the
two measures of permeability is K ˆ …mL=…grLKL†a. Experimental values for
permeability vary in a very wide range, about four orders of magnitude [50, 51].
Permeability of a ˆ 3 £ 10¡13 m=s for high performance concrete with com-
pressive strength of 54 MPa is reported in [51]. It is equivalent to
K & 5 £ 10¡17 m2, because at ambient conditions mL ˆ 0:001 kg=…m:s†.
rL ˆ 0rL ˆ 60 kg=m3, and KL ˆ 0:01 (for the magnitude of the latter see the note
where KL is de¢ned). The value agrees with Ref. [48], where permeability of
8:3 £ 10¡17 m2 is used for cement mortar with higher porosity twice that of
the one used here.
TRANSFER IN CONCRETE SUBJECTED TO FIRE 707
APPENDIX B (DERIVATION OF EQUATIONS)
B1. Water vapor £ux:
¡JV ˆ
KKG
mG
eG ~rVrPG ‡ eG ~rGDVAr
~rV
~rG
…B1a†
r
~rV
~rG
ˆ …
~rA ‡ ~rV †r ~rV ¡ ~rVr… ~rA ‡ ~rV †
~r2G
ˆ
~rAr ~rV ¡ ~rVr~rA
~r2G
…B1b†
r ~rA ˆ r
PA
RAT
ˆ r PG ¡ PV
RAT
ˆ r PG
RAT
¡ RV
RA
~rV
ˆ
TrPG ¡ PGrT
RAT 2
¡
RV
RA
r ~rV
r ~rA ˆ ¡
PG
RAT2
rT ‡ 1
RAT
rPG ‡ ¡
RV
RA
r ~rV …B1c†
~rGr
~rV
~rG
ˆ
~rVPG
~rGRAT2
rT ‡ ¡
~rV
~rGRAT
rPG ‡
~rA
~rG
‡
~rV RV
~rGRA
r~rV …B1d†
¡JV ˆ
DVAeG ~rV PG
~rGRAT2
rT ‡ eG ~rV
KKG
mG
¡ DVA
~rGRAT
rPG
‡
DVAeG
~rG
~rA ‡
RV
RA
~rV r~rV …B1e†
¡JV ˆ KVTrT ‡ KVPrPG ‡ KVVr~rV …B1f †
B2. Air £ux:
¡JA ˆ
KKG
mG
~eG ~rArPG ‡ ~eG ~rGDAVr
~rA
~rG
…B2a†
r
~rA
~rG
ˆ r
~rG ¡ ~rV
~rG
ˆ ¡r
~rV
~rG
using Eq: …B1d† : …B2b†
¡JA ˆ ¡
DAV eG ~rV
~rGRAT
PG
T
rT ‡ KKG
mG
eG ~rA ‡
DAV eG ~rV
~rGRAT
rPG
‡ ¡ ~rA ¡
RV
RA
~rV
DAV eG
~rG
r ~rV …B2c†
¡JA ˆ KATrT ‡ KAPrPG ‡ KAVr ~rV …B2d†
708 R. T. TENCHEV ET AL.
B3. Water £ux:
¡JL ˆ rL
KKL
mL
rPL …B3a†
¡JL ˆ …0†rT ‡ KLPrPG ‡ …0†r ~rV …B3b†
B4. Moisture (liquid water and water vapour) £ux:
¡…JV ‡ JL† ˆ …KVT ‡ 0†rT ‡ …KVP ‡ KLP†rPG ‡ …KVV ‡ 0†r ~rV …B4a†
¡…JV ‡ JL† ˆ KMTrT ‡ KMPrPG ‡ KMVr ~rV …B4b†
B5. Time derivative:
@…eG ~rA†
@t
@…eG ~rA†
@t
ˆ eG
@ ~rA
@t
‡ ~rA
@eG
@t
…B5a†
~rA ˆ
PA
RAT
ˆ PG ¡ PV
RAT
ˆ PG ¡
~rV RV T
RAT
ˆ 1
RA
PG
T
¡ RV
RA
~rV …B5b†
@ ~rA
@t
ˆ
T @PG
@t ¡ PG
@T
@r
RAT 2
¡ RV
RA
@ ~rV
@t
ˆ ¡ PG
RAT2
@T
@t
‡ 1
RAT
@PG
@t
‡ ¡ RV
RA
@ ~rV
@t
…B5c†
@eG
@t
ˆ
@por
@T
¡
1
rL
@rL
@T
@T
@t
‡ …0†
@PG
@t
‡ ¡
1
rL
@rL
@ ~rV
@ ~rV
@t
…B5d†
@…eG ~rA†
@t
ˆ ~rA
@por
@T
¡
~rA
rL
@rL
@T
¡ eGPG
RAT 2
@T
@r
‡ eG
RAT
@PG
@t
‡ ¡eG
RV
RA
¡
~rA
rL
@rL
@ ~rV
@ ~rV
@t
@…eG ~rA†
@t
ˆ CAT
@T
@t
‡ CAP
@PG
@t
‡ CAV
@ ~rV
@t
…B5c†
B6. Time derivative:
@…eG ~rV †
@t
@…eG ~rV †
@t
ˆ eG
@ ~rV
@t
‡ ~rV
@eG
@t
…B6a†
~rV
@eG
@t
ˆ ~rV
@por
@T
¡
~rV
rL
@rL
@T
@T
@t
‡ …0†
@PG
@t
‡ ¡
~rV
rL
@rL
@ ~rV
@ ~rV
@t
…B6b†
@…eG ~rV †
@t
ˆ ~rV
@por
@T
¡
~rV
rL
@rL
@T
@T
@t
‡ …0†
@PG
@t
‡ eG ¡
~rV
rL
@rL
@ ~rV
@ ~rV
@t
…B6c†
TRANSFER IN CONCRETE SUBJECTED TO FIRE 709
B7. Left-hand side:
@…eG ~rV †
@r
‡ @rL
@t
¡ @rD
@t
@rL
@t
ˆ
@rL
@T
@T
@t
‡ …0†
@PG
@t
‡
@rL
@ ~rV
@ ~rV
@t
@rD
@t
ˆ
@rD
@T
@T
@t
…B7a†
@…eG ~rV †
@t
‡ @rL
@t
¡ @rD
@t
ˆ ~rV
@por
@T
‡ 1 ¡
~rV
rL
@rL
@T
¡ @rD
@T
@T
@t
‡ …0†
@PG
@t
‡ eG ‡ 1 ¡
rV
rL
@rL
@ ~rV
@ ~rV
@t
@…eG ~rV †
@t
‡ @rL
@t
¡ @rD
@t
ˆ CMT
@T
@t
‡ CMP
@PG
@t
‡ CMV
@ ~rV
@t
…B7b†
B8. Left-hand side: …rC†
@T
@t
‡ …lD ‡ lE†
@rD
@t
¡ lE
@rL
@t
…rC†
@T
@t
‡ …lD ‡ lE†
@rD
@t
¡ lE
@rL
@t
ˆ …rC† ‡ …lD ‡ lE†
@rD
@T
¡ ¡lErL
@rL
@T
@T
@t
‡ …0†
@PG
@t
¡ lErL
@rL
@T
@ ~rV
@t
…rC†
@T
@t
‡ …lD ‡ lE†
@rD
@t
¡ lE
@rL
@t
ˆ CTT
@T
@t
‡ CTP
@PG
@t
‡ CTV
@ ~rV
@t
APPENDIX C (FINITE ELEM ENT MATRIX NOTATION)
Ce ²
CTT CTP CTV
CAT CAP CAV
CMT CMP CMV
2
64
3
75 Ke ²
KTT KTP KTV
KAT KAP KAVr
KMT KMP KMV
2
64
3
75 …C1†
KeV ²
…rCv† 0 0
0 0 0
0 0 0
2
64
3
75 U ²
T
PG
~rV
8
><
>:
9
>=
>;
…C2†
Re1 ²
RT ;1
0
RV ;1
8
<
:
9
=
; where
RT ;1
RV ;1
ˆ KTT 0
KVT KVV
¡1 hqrT1
b ~rV;1
…C3†
ReK ˆ
KhT 0 Khr
0 0 0
KbT 0 Kbr
2
64
3
75
where
KhT Khr
KbT Kbr
ˆ
KTT 0
KVT KVV
¡1 hqr 0
0 b
…C4†
N ²
bNc 0 0
0 bNc 0
0 0 bNc
2
64
3
75
where bNc ˆ bN1…x; Z† N2…x; Z† ¢ ¢ ¢ Nn…x; Z†c …C5†
710 R. T. TENCHEV ET AL.
Copyright of Numerical Heat Transfer: Part A -- Applications is the property of Taylor &
Francis Ltd and its content may not be copied or emailed to multiple sites or posted to a
listserv without the copyright holder's express written permission. However, users may print,
download, or email articles for individual use.