Multi-Physics Modeling of Fire-Induced Damage in High-Performance Concrete

The response of high-performance concrete (HPC) to rapid heating due to exposure to fire with air temperatures exceeding 700°C is analyzed. The analysis focuses on coupled thermal-mechanical-transport processes at the mesoscale in the first 10 minutes of exposure during which heating rates on the order of 10 °C/minute and temperature gradients on the order of 15 °C/cm are involved. The driving forces for damage considered are thermomechanical stresses and internal pore pressure resulting from the expansion of water within the material. The HPC is modeled as a two-phase composite consisting of a cementitious matrix and a population of quartz aggregate particles. Mesostructures with aggregate sizes of 400-1600 µm and aggregate volume fractions of 10-30% are considered. To capture the development of stresses and pore pressure, the cementitious matrix is modeled using a coupled thermal-mechanical-transport formulation and the aggregate is modeled using a thermal-mechanical formulation. Simulations show that the composition of the mesostructures significantly influences the time and spatial distribution of damage. Materials with smaller aggregate sizes and the lower effective permeability are found to exhibit more rapid property degradation. The time to failure and depth of thermal spall are quantified as functions of structural variables. This framework and the mesostructure-response relations obtained serves as a tool for the design of HPC that are more resistant to fire-induced damage.


INTRODUCTION
Concrete is a porous, reactive, composite material. Chemical transformations occurring during heating include the evaporation of free water; the dehydration and subsequent evaporation of bound water in calcium silicate hydrates (C-S-H), ettringite, and gypsum; the dehydroxylation of calcium hydroxide; and the decarbonatation of calcium carbonate [1]. Solid-state phase transitions can also occur, including the transformation of quartz aggregate from its triagonal α phase to its hexagonal β phase and the transformation of steel from its body-centered cubic ferrite phase to its face-centered cubic austenite phase. The constituents of concrete have differing coefficients of thermal expansion which lead to short-range thermal stresses. Additionally, large thermal gradients arising from the low thermal conductivity of the material lead to long-range stresses. The combination of chemically and thermally induced stresses when concrete is exposed to fire can lead to property degradation, damage, and failure. Obviously, the response of concrete to rapid heating is complex and multi-physics in nature and the thermal, mechanical and transport processes governing the response are coupled.
Hertz first reported the explosion of silica fume concrete during rapid heating [2]. Experiments on high-performance concretes (HPCs) have demonstrated that the build-up of internal pore pressure along with other factors lead to explosive spallation [3,4]. It is postulated that the thermal expansion of liquid and gaseous water and their transport towards the inner part of specimens play a significant role in the build-up of gas pressure. This occurs through a quasi-saturated layer that precedes the drying front and acts as a moisture clog. Internal measurements show that pore pressures reached prior to spall in HPCs range from 1 MPa to 5 MPa [3,4] which are on the order of the tensile strength of the materials. This behavior is observed at temperatures of approximately 200-270°C, which are reached in the first 20 minutes of fire exposure [3,4].
The coupled thermo-poro-mechanical constitutive behavior of homogenous concrete has been computationally modeled [5][6][7][8]. The framework has been used to predict the time and location of thermal spall during heating [9,10]. Models have idealized concrete as a threeor four-phase homogenized media consisting of a solid cementitious skeleton, liquid water, water vapor, and dry air. Conservations of mass, energy, and linear momentum provide the governing equations for these models. High-temperature phenomena such as the triple point of water have also been considered [11]. Homogenized capillary porosity and cracks are considered as the two primary networks for fluid transport that relieve high localized pressures.
Thermally-induced damage in cement-based composites was first thermo-mechanically modeled at the mesoscale by Fu et al. to capture cracking around a single heterogeneity [12] or multiple heterogeneities [13]. The authors studied the influence of the heterogeneous constituents' coefficients of thermal expansion, size scale, and spatial distribution on stresses and cracking. Willam et al. studied thermo-mechanical coupling and degradation, accounting for volumetric and deviatoric damage evolution [14]. Grondin et. al. implemented a mesoscale model to calculate homogenized thermo-mechanical properties and degradation [15]. Grassl and Pierce implemented a thermo-mechanical lattice model to capture not only bulk phases, but also the interfacial transition zone (ITZ) between the cementitious matrix and aggregate [16]. They used the model to predict cracking and structural stability under compressive loading and uniaxial constraint. Recent efforts by Xotta et. al. have utilized three-dimensional models to study concrete repair after thermal exposure [17] and the effect of aggregate on the visco-damage response of concrete [18]. Their hygral-thermal study concerned the transport response under uniform heating of concrete systems with several large aggregates. The size of the models allowed for explicit modeling of the interfacial transition zone. So far, no mesoscale study on the thermo-poro-mechanical behavior and damage development has been reported.
Here, we develop a three-dimensional framework for analyzing the thermo-poromechanical response of HPC mesostructures. Mesostructures are comprised of a population of quartz aggregate and a cementitious matrix. The matrix is a homogenized multi-phase medium consisting of a cementitious skeleton, pores, and liquid water. The analysis focuses on the effect of mesostructure on material response. The primary structural parameters considered are aggregate volume fraction, aggregate size, and permeability of the cementitious matrix. The aggregate volume fraction is varied from 10% to 30% and the aggregate size is varied from 400 µm to 1600 µm. The permeability is varied between 10 -19 m 2 and 10 -21 m 2 , reflecting a range of values obtained from experiments on ordinary and high-performance concretes [19,20]. Simulations carried out concern the first 10 minutes of exposure to a standard ASTM fire in which air temperatures exceed 700 °C. The results are used to formulate empirical relations which can be used to design concrete mixes that are more fire-resistant.

MESOSCALE FINITE ELEMENT MODEL
The choice of constituent size and volume fraction is motivated by the characteristics of real HPC mesostructures [21]. A micrograph of a representative mesostructure obtained using backscattered electron diffraction under a scanning electron microscope is shown in Figure  1 [21]. The quartz aggregate particles are gray, steel fibers are white, and pores are black. This HPC sample was prepared with Lafarge Ductal® premix, a water-reducing superplasticizer, and metallic fibers. The diameter of the aggregate phase ranges up to 800 µm, this phase accounts for approximately 40% of the total volume of the material. The steel fibers have a diameter of 200 µm and a volume fraction of 2%. The pores are typically 200-1,000 µm in size and account for 2% of the structure by volume. This sample was cured at high-temperature with high-humidity to maximize hydration. Reinforcing steel rebar and fibers are not considered here and porosity is homogenized into the matrix.
To provide a systematic sampling of structural variation, two-phase mesostructures are computationally generated with randomly placed aggregate particles. Each sample measures 20 mm × 20 mm × 20 mm. The volume fraction of quartz aggregate is 10%, 20%, or 30% and the aggregate diameter is 400 µm, 800 µm, 1600 µm, or a combination of these three diameters in equal volumetric proportions. Particles are spherical in shape and morphology changes are not considered. The resulting twelve structural variants are shown in Figure 2. The structures are mapped onto a regular, hexahedral finite element mesh and imported into Abaqus/Standard for finite element method calculations. Each mesh consists of 125,000, 8noded, hexahedral finite elements (type C3D8T or C3D8PT) and 2,500, 18-noded, infinite elements (type CIN3D18R). The resulting mesh resolution is 400 µm per element. The mesh with boundary and loading conditions applied to each face is shown in Figure 3.   Periodic boundary conditions (PBCs) are applied on faces 3 and 4 and faces 5 and 6 of the model. The application of the PBCs is natural and straightforward. For all nodes with the same x 1 and x 3 coordinates on faces 3 and 4 and the same x 2 and x 3 coordinates on faces 5 and 6, the degree of freedom vectors, u 3 and u 4 and u 5 and u 6 , respectively, are constrained such that, where the degree of freedom vector on face i is , with u 1 , u 2 , and u 3 being the nodal displacements in the x 1 , x 2 , and x 3 directions, respectively, u T being the nodal temperature, and u P being the nodal pore pressure. The resulting deformation, temperature, and pressure patterns on opposite faces of the model are thus, identical.
To model each sample as part of the surface portion of a larger structure, infinite elements with linear-elastic properties are placed on face 2. These elements account for the effects of the rest of the larger structure and obviate the need of using more strict constraints such as zero-displacement boundary conditions. On this face, adjacent finite and infinite elements are assigned identical linear-elastic material properties. The infinite medium is modeled using interpolation functions which decay by ξ 3 -3 , where ξ is the local coordinate system of the element, such that as ξ 3 → ∞ the displacement and strain within the element tends to zero. The prescribed infinite elements can only model linear material behavior. This treatment is considered a reasonable approximation because significantly far away from the heated surface damage does not occur. In addition to the infinite elements, surface flux and flow conditions are specified at face 2 to account for diffusive and convective transport into the remaining part of the larger medium. Specifically, the heat flux, q, at any node on face 2, in the -x 3 direction, is (1) where h is the heat transfer coefficient of the surface, u T 2 is the nodal temperature and θ°is a prescribed far-field temperature which is taken as 20 °C in the calculations carried. The heat transfer coefficient on the surface is chosen to be 1.5 W m -2 K -1 to simulate the effects of heat transfer with a concrete medium. The assumption that the far-field temperature remains constant is less accurate at long times, but is quite reasonable for the first 10 minutes of exposure modeled here for distances beyond 0.5 m from the exposed surface. Fluid flow through Face 2 in the -x 3 direction follows the relation of (2) where v is the fluid velocity, k is the seepage coefficient of the surface, u P 2 is the nodal pore pressure, and P 0 is the far-field pore pressure. An intrinsic permeability of κ = 1×10 -19 m 2 is used for the matrix. A value of 2.6×10 -15 m 3 N -1 s -1 is then calculated for the seepage coefficient, k, as recommended in the Abaqus Documentation [22]. The far-field pressure, P 0 , is chosen to be 0 MPa. This assumption is valid for the first several minutes of analysis, similar to the far-field temperature assumption stated earlier.
Exposure to convective heating from fire occurs on face 1. The temperature-time relation used in the ASTM E119 standard [23] is used to prescribe the temperature of the heated air in contact with the structure. Heat transfer into the mesostructure occurs through a film on face 1, where the heat flux is determined by Equation (2). On this surface, θ 0 is the ASTM prescribed air temperature and h is 8 W m -2 K -1 , the heat transfer coefficient between air and concrete. Mechanically, face 1 is traction-free and has a pore pressure of 0 MPa.
The initial conditions are as follows. There are no applied pre-stresses, the initial temperature is 20 °C, the cementitious matrix is 100% saturated, and there is an initial void ratio of 0.20. The void ratio, e, is defined as, where n is the porosity. A summary of all boundary and initial conditions is given in Table  1. Periodic displacement boundary -3 Periodic temperature boundary -3 Periodic pressure boundary -3 5 Periodic displacement boundary -6 Periodic temperature boundary -6 Periodic pressure boundary - 6 6 Periodic displacement boundary -5 Periodic temperature boundary -5 Periodic pressure boundary -5

QUARTZ AGGREGATE
Diffusive heat conduction in the quartz aggregate is modeled using Fourier's Law such that the heat flux, q, is where m is the conductivity tensor of the material and ∇θ is the spatial gradient of the temperature field. Under the assumption of isotropic heat conduction, m = mI and the flux simplifies to q = -m∇θ, where is the isotropic conductance of the material and I is the identity tensor. Thermal equilibrium is then stated as (5) where ρ and c are the density and heat capacity of the material, respectively, ∆ is the Laplacian operator, and q s is the surface flux per unit area through S. The quartz aggregate is considered to be linear-elastic; dislocation glide in quartz crystals is not a prevalent deformation mechanism until temperatures exceed 800 °C [24]. A simplified material model can thus be used to provide an adequate account of the material response. Thermal and mechanical responses are coupled such that the stresses in the quartz aggregate are determined as (6) where L is the fourth-order elastic stiffness tensor, ε is the infinitesimal strain tensor, α and . Air time-temperature relation specified in the ASTM E119 standard and used in the analyses [24] K are the coefficient of thermal expansion and bulk modulus of the material, respectively, and θ and θ 0 are the current and reference temperatures, respectively. In the absence of body forces and surface tractions, equilibrium can then be stated as (7) For an isotropic, elastic material, L is completely defined by two parameters. In the case of quartz, the elastic modulus, E, and the Poisson Ratio, ν, are 95 GPa and 0.10, respectively. Equations (6) and (7) are stated in the weak form and solved for nodal displacements and temperatures. The coefficient of thermal expansion is chosen to be constant, at 15×10 -6 °C -1 . The density of quartz is 2.65 g cm -3 . The thermal conductivity is 15 W m -1 °C -1 at room temperature. This value decreases linearly to 4.0 W m -1 °C -1 at the α-β transition temperature of 573 °C and remains constant at that value thereafter.

MATRIX
The cementitious matrix is modeled as an isotropic, temperature-dependent thermo-poroelastic, pressure-dependent thermo-poro-plastic, homogenized continuum comprised of a solid skeleton, pores, and fluid. Accounting for the presence of fluid flowing through the matrix, the thermal equilibrium over volume V with boundary surface S in the absence of heat generation is (8) where ρ w and c w , and ρ s and c s are the density and heat capacity of the fluid and solid, respectively, v w is the fluid velocity, and q s is the conductive heat flux through S. Fluid flow in the matrix is modeled using Forchheimer's Law; which states that high flow velocities reduce the effective permeability of the medium. For low velocities, it closely approximates the more well-known Darcy's Law. The Forchheimer flow equation gives the volumetric flow vector, f, under conditions of no gravitational effects as (9) where γ w and p w are the specific weight and pore pressure of the fluid. The effective permeability is defined as k -= k(1+λ v w ) -1 , where k is the permeability tensor of the matrix, λ is a velocity coefficient dependent on the void ratio, and v w is the fluid velocity vector. Conservation of pore fluid mass requires that the change in mass within V be equal to the mass exiting S, such that (10) where n is the unit normal to S, J is the Jacobian of the material deformation defined as J≡det(dy/dx)=dV/dV 0 , y and dV are the material coordinates and differential volume in the current configuration, respectively, and x and dV 0 are the corresponding coordinates and volume in the reference configuration, respectively. The effective stress in the porous medium is defined as (11) where d t is a dimensionless, scalar measure of isotropic damage representing the reduction of the effective load-carrying area within the volume due to cracking, and ε′ is the inelastic strain tensor. In the absence of applied body forces and surface tractions, equilibrium can be stated as, The magnitude of p w is determined solely by strain incompatibility between the thermally expanding and slightly compressible pore fluid and the thermally expanding and deformable solid matrix. Thermal equilibrium, fluid mass continuity, and mechanical equilibrium provide the coupled governing differential equations, Equations (9), (11), and (14). The equations are then stated in a weak form and solved for nodal displacements, temperatures, and pore pressures.
The inelastic deformation of the cementitious matrix is allowed to evolve in an anisotropic manner to capture both tensile damage and dilatation during pressure-dependent plastic flow. The mathematical development of the constitutive model follows that of Lubliner et. al. [25]. The plastic-damage model assumes a non-associative plastic flow potential in the form of (13) where e is eccentricity; σ t0 = σ t0 (θ) is the uniaxial tensile stress at failure as a function of temperature; ψ is the dilation angle; q = √  3 2 S:S is the Mises equivalent stress, and p=- 1 3 tr(σ) is the hydrostatic pressure.
Stress in the solid skeleton is defined as The tensile and compressive equivalent plastic strain rates, ε · t pl and ε · c pl , are defined as (15) where u=x-X and F= ∂x ∂X are the maximum and minimum eigenvalues of the plastic strain rate tensor ε ·pl , respectively. The parameter r(σ ) is a stress weight factor valued between zero and one and defined as where σ i is the i th principal stress and 〈•〉 is the McCauley bracket such that 〈σ 〉= 1 2 (σ +|σ |). Tensile cracking data is input into the model to calculate the amount of recoverable elastic strain, permanent plastic strain, and unrecoverable cracking strain or damage at each time step.
A pressure-dependent yield criterion is chosen such that , where ω and is dimensionless material constant; ; σ t is the tensile cohesion stress; and σ c is the compressive cohesion stress. The parameter ω is defined as (18) where σ b0 is the initial yield stress for equi-biaxial compression and σ c0 is the initial yield stress for uniaxial compression. The uniaxial stress-strain behavior of the cementitious matrix phase in compression and tension is shown in Figure 5 [26].

NUMERICAL RESULTS
The first 10 minutes of fire exposure are simulated and the representative thermal, mechanical, and hygral response is presented, followed by an analysis of the effect of mesostructure on the damage response.

THERMAL RESPONSE Calculated surface temperatures agree with measurements performed by Jenson and van
Horn on HPC specimens in experiments [27], providing a validation for the methodology of calculating the heat flux at the exposed surface. The calculated interior temperatures range from 55°C to 65°C, consistent with thermocouple measurements obtained at 20 mm to 25 mm from the exposed surface in experiments [1,27]. The evolution of temperature within mesostructures with 10%, and 20%, and 30% aggregate are shown in Figure 6(a), Figure 6(b), and Figure 6(c), respectively. Mesostructures with equal volume fractions of each aggregate size are denoted "Sieve" in the plots and from hereon. The average thermal gradients within the structures are found to range from 8 °C cm -1 to 19 °C cm -1 . The magnitude of the gradients is observed to increase with time and vary with mesostructure characteristics. The average thermal gradient within the mesostructures as a function of aggregate size and volume fraction is shown in Figure 7. Structures with higher aggregate content are found to have lower thermal gradients. This result is logical considering that the conductivity of the aggregate is higher than that of the matrix. It is also observed that the average diameter of the aggregate has a secondary effect on the average temperature gradient. Specifically, lower average diameter aggregates decrease the average thermal gradient. This effect will be analyzed in more detail later along with internal pore pressures.  Mesostructures with low aggregate contents have higher thermal gradients due to the higher diffusive conductivity of the aggregate phase compared to that of the cementitious matrix phase.

THERMO-MECHANICAL STRESSES
The average σ 11 in the matrix as a function of depth into the structure during heating is shown in Figure 8(a), Figure 8(b), and Figure 8(c) for mesostructures with 10%, 20%, and 30% aggregate, respectively. Compressive stresses occur near the surface of the structure; a transition to tensile stresses is observed in the interior. The addition of aggregate alleviates stresses in the matrix phase. Specifically, an increase from 10% to 30% in volume fraction results in an approximate 15% reduction in the average compressive stress in the matrix. The stresses first develop around the aggregate, primarily because of the differences in thermal expansion between the matrix and the aggregate. As time progresses, regions of higher compressive stresses grow away from the aggregate. This development toward longer-range, more homogeneous distributions of stresses is related to the confinement the mesostructure experiences. Due to mesostructural level isotropy in the x 1x 2 plane, the results for σ 11 and σ 22 are essentially the same, as expected.
Similarly, the average σ 11 as a function of depth is calculated in the aggregate phase. The evolution of these stresses can be found in Figure 10. The stress component σ 11 in the aggregate phase is observed to be compressive throughout the analysis and does not show a compression-tension transition point near the center of the structure as the σ 11 in the matrix phase does. The magnitude of the average compressive stress in the aggregate is also higher than that in the matrix by 20-25%. The stress fields in the matrix agree with the results of prior numerical studies. However, the results here suggest that tensile σ 11 and σ 22 in the aggregate and across the matrix-aggregate interface contribute to damage evolution and spall. The distribution of σ 33 , the out-of-plane stress component, is found to differ significantly from those of the in-plane stresses. The evolution of σ 33 within the matrix can be found in Figure 11. In general, this stress component is lower in magnitude than σ 11 and σ 22 . σ 33 is typically compressive near the surface and tensile in the interior. The transition from compression to tension is near the midpoint of the specimen depth. Structures with larger aggregates are observed to generate higher (more positive) σ 33 than those with smaller aggregates. Conversely, σ 33 in the aggregate is tensile throughout the heating process and structures with smaller aggregates have higher (more tensile) σ 33 .

PORE PRESSURE
The evolution of internal pore pressure is shown in Figure 13. For all cases, the pore pressure is less than 1 MPa, in agreement with experimental measurements [4]. A strong correlation between pore pressure and the volume fraction of the aggregate phase is observed. Specifically, higher aggregate volume fractions lead to higher internal pore pressures. Naturally, mesostructures with increasingly fewer impermeable phases have higher overall permeability and faster fluid flow. Aggregate size has only a minor effect on fluid transport, with larger aggregates of 1600 µm diameter producing approximately 2% to 5% lower pore pressure values than smaller aggregate of 400 µm diameter. Thus, it can be assumed that transport is largely governed by the effective permeability of the matrix.
The observation that larger aggregate lead to lower pore pressures and higher temperature gradients is explained by examining the mean free path between aggregate particles in the x 1 direction. The lineal-Path correlation function profiles for all structures analyzed are shown in Figure 14. The lineal-path correlation function, L(x) I , is the probability of a particle traveling a distance x in the x 1 direction without encountering an aggregate particle. This probability is a function of x. Structures with smaller mean aggregate sizes have more finely distributed aggregate and are therefore have lower probabilities of allowing a particle to travel a certain distance x without encountering the aggregate phase. The lower average thermal gradients can be explained with similar consideration. Thermal energy propagating in a mesostructure containing more and finer aggregate has a higher probability of encountering a more highly conductive aggregate phase than in a mesostructure with the same volume fraction of larger aggregate. For this reason, mesostructures with smaller aggregates have lower effective permeability and higher effective conductivity at a given volume fraction.

DAMAGE EVOLUTION
The evolution of averaged sub-surface damage as a function of depth within all mesostructures considered is shown in Figure 14. The damage is highly inhomogeneous spatially and is, therefore, averaged over planes parallel to the exposure surface. The highest values are observed in mesostructures containing fewer, larger aggregate. The mesostructure Figure 13. Lineal-Path correlation function profiles for all mesostructures analyzed. Mesostructures with (a) 10% volume fraction of aggregate, (b) 20% volume fraction of aggregate and (c) 30% volume fraction of aggregate. denotes the probability of a fluid particle's traveling a distance x in the x 1 direction without encountering the aggregate phase. decreases as aggregate volume fraction increases and as the mean aggregate size decreases.
with 10% 1,600 µm diameter aggregates has the highest level of damage. Note that this mesostructure does not show the highest stresses in Figure 9, Figure 10, Figure 11, or Figure  12, nor the highest pore pressures in Figure 13. This structure does show the highest thermal gradients. This result shows that damage during rapid heating prior to the evaporation of water is primarily driven by bulk thermal gradients, not necessarily the absolute levels of stress. Figure 15 shows the distributions of tensile damage in mesostructures containing 3.33% 400 µm, 3.33% 800 µm, and 3.33% 1600 µm aggregate. Damage primarily nucleates at the matrix-aggregate interfaces and propagates in the x 1 -x 2 plane. It can be seen that tensile σ 33 is a primary driving force for damage development.

MATRIX PERMEABILITY AND MESOSTRUCTURAL EFFECT ON THE DEVELOPMENT OF CRACKS
Three levels of intrinsic permeability of the cementitious matrix (1×10 -21 m 2 , 1×10 -20 m 2 , and 1×10 -19 m 2 ) are considered using the mesostructures shown in Figure 2. These values span the range of intrinsic permeability measured in ordinary and high-performance concretes [19,20]. While the average damage may be a more reliable measure of dissipation through fracture or loss in bulk stiffness, it does not accurately quantify or predict the initiation, coalescence, or propagation of cracks. Cracks can either serve as pathways for pressurized water if they exist prior to the evaporation of free water and the dehydration and evaporation of bound water, or serve as sites for spallation. Clearly, the development of cracks prior to the evaporation of water is important in mitigating the effects of pressurized water, preventing spallation, and maintaining structural stability. Three critical levels (0.1, 0.2, and 0.3) of element damage, d t , are chosen. It is observed that these values are first reached along planes parallel to the heated surface. The initiation time of a crack plane is defined as the time when the area of a network of adjacent elements, all meeting the critical damage level criterion, exceeds 1 mm 2 . These planes subsequently coalesce near aggregates, especially large aggregates in the case of the "Sieve" mesostructures. The average depth and time of initiation of these planes are shown for all 36 mesostructures considered in Figure 16 and Figure 17, respectively. To account for the statistical variations between structures and the resulting crack planes at each structural setting, five instantiations of each structural variant were generated and used. When the permeability of the matrix is higher than 1×10 -20 m 2 , no cracking is seen in mesostructures containing only small aggregates. In contrast, structures consisting primarily of large or a mixture of large and small aggregates show spall cracks at depths of 10-13 mm within the first 8 minutes of fire exposure. From Figure 16 and Figure 17, it can be seen that reduction in matrix permeability, increase in aggregate volume fraction, or reduction in aggregate size causes cracking to occur closer to the exposure surface. An increase in aggregate volume fraction from 10% to 30% causes the depth of spall crack initiation to decrease by 70%, from 8.99 mm to 2.70 mm. Reducing the aggregate size from 1600 µm to 400 µm results in a 27% average reduction in crack depth, from 15.9 mm to 11.6 mm. The corresponding time to cracking also decreases. Shallow cracks, generated in mesostructures with 30% aggregate volume fraction, occur within the first 3 minutes of fire exposure. Conversely, mesostructures with 10% aggregate volume fraction require 7-9 minutes of exposure to initiate a crack. Additionally, it was observed that for some combinations of aggregate volume fraction, size, and matrix permeability, no cracking occurred. The depth and time of spall crack initiation closely correlate with variation in the mesostructures. To quantify the relation, a scaling law between the depth of spall, r s , the permeability of the matrix, κ, the aggregate diameter, d a , and the aggregate volume fraction, v f , is developed. This relation can be stated as (19) where α, β, γ, and δ 0 are empirical material parameters depending on mesostructural variables. A similar relation is obtained for the time to spall crack initiation, i.e., (20) where α, β, γ, and t* are empirical material parameters depending on mesostructural variables and is a time scaling factor. Table 2 shows the values of the parameters for the above relations. One insight gained from analyzing the scaling parameters is on the influence of aggregate size on the depth of spall. Specifically, mesostructures with larger aggregate show spalls closer to the exposed surface. level of at least 0.001 (this level is much lower than the damage levels for the occurrence of spall considered earlier). The tracer analysis yields the full three-dimensional geometry of the crack network and salient characteristics such as the location and length of cracks. The results of the tracing analysis are shown in Figure 18. Large aggregates, particularly at high volume fractions, induce the longest crack networks. For example, in the most extreme case (30% volume fraction of 1.6 mm diameter aggregates) the crack length is 9.6 mm.
Mesostructures with small aggregates show only localized damage which is limited to areas on the order of 1-2 aggregate diameters along vertical planes. The coalescence of damaged regions prior to the evaporation of water provides pathways for the relief of pressurized gas.

DISCUSSION
During heating of concrete, the macro-scale stiffness and tensile strength decrease even at temperatures below 100 °C. The degradation of properties has been linked to micro-cracking [28]. The results obtained through numerical simulations in this paper are consistent with the results obtained from experiments. The constitutive relations used are able to track this degradation during rapid heating of the thermo-poro-mechanical matrix and thermomechanical aggregate. While the temperature regimes considered here are lower than the boiling point of water and the typical point of thermal spall, the damage created in the first minutes of rapid heating serves as transportation networks for hot fluids and gases at higher temperatures. Structures with larger aggregate are especially more susceptible to damage, with the spall planes closer to the exposed surface of the structure. The formation of nearsurface cracks may serve as an important mechanism for avoiding explosive spallation and maintaining structural stability during a fire by providing channels for water vapor transport. In a previous research, the permeability of the cementitious skeleton was taken to be dependent on both chemical porosity and mechanical damage, since both act as networks for fluid flow [30]. Since the effects of damage and chemical porosity on the permeability of Figure 18. The longest crack network length in the direction perpendicular to the exposed surface as a function of mesostructure parameters. Large aggregates lead to longer crack lengths.
the skeleton are not well-understood, only estimates can be made concerning the degree of their influences. The evaluation of these internal state variables and the quantification of their influences on pore fluid transport are important. This task will be undertaken in future analyses.
It is important to note that the interfaces between the matrix and aggregate are assumed to be perfect here. The interfacial bonding between the matrix and aggregate is known to be stronger and the thickness of the ITZ is known to be thinner in HPC relative to those in ordinary concrete. Thus, perfect thermal and mechanical property assumptions may be reasonable approximations here. If the ITZ is large, or the bonding between the matrix and aggregate is weak, special treatment of the interface may be necessary. The analyses of σ 11 in Figures 9-10 and σ 33 in Figures 11-12 indicate that large stress gradients exist between the two phases. Interfacial failure between the two phases may very well be the dominant damage initiation mechanism during rapid heating. After the interface failure, cracks will propagate into the matrix, a progress governed by the thermo-poro-mechanical stress state. Additionally, the interfacial failure sites can serve as vital pathways for the flow of high pressure fluid and gas. These phenomena should be considered in future studies in which the evaporation and dehydration of water is considered.

CONCLUSIONS
The thermo-poro-mechanical response of HPCs at the mesoscale to high-rate heating from fire is computationally studied. The effect of aggregate size and volume fraction on thermal gradients, thermo-mechanical stresses, pore pressure, and damage is assessed. Structures with higher aggregate volume fractions are found to have lower thermal gradients, lower tensile stresses in the aggregate phase, and higher tensile stresses in the matrix phase. Higher matrix phase stresses lead to the development of matrix damage, which is found to nucleate in the vicinity of larger aggregates. The damage networks generated in this low temperature regime are found to be important for relieving high pressures in the fluids and gases that arise when free water vaporizes and bound water dissociates and subsequently vaporizes. The model and the analysis presented here represent a first step in understanding the mesostructural evolution in high-performance concretes up until the time of spall.