Abstract

In the past three decades, one-dimensional (1D) thermal model was usually used to estimate the thermal responses of glass fiber-reinforced polymer (GFRP) materials and structures. However, the temperature gradient and mechanical degradation of whole cross sections cannot be accurately evaluated. To address this issue, a two-dimensional (2D) thermomechanical model was developed to predict the thermal and mechanical responses of rectangular GFRP tubes subjected to one-side ISO-834 fire exposure in this paper. The 2D governing heat transfer equations with thermal boundary conditions, discretized by alternating direction implicit (ADI) method, were solved by Gauss-Seidel iterative approach. Then the temperature-dependent mechanical responses were obtained by considering the elastic modulus degradation from glass transition and decomposition of resin. The temperatures and midspan deflections of available experimental results can be reasonably predicted. The overestimation of deflections could be attributed to the underestimation of bending stiffness. This model can also be extended to simulate the thermomechanical responses of beams and columns subjected to multiside fire loading, which may occur in real fire scenarios.

1. Introduction

Compared with traditional building materials, fiber-reinforced polymer (FRP) materials have many obvious advantages such as high strength-to-weight ratio, superior durability, good fatigue endurance, and rapid installation. Hence, the FRP composites have increasingly been used in civil engineering due to these excellent properties. However, the resin in the FRP composites undergoes a glass transition stage when subjected to the high temperatures around the glass transition temperatures (). The resin matrix changes from a rigid, glassy solid to a softer and viscoelastic material during the glass transition stage. The mechanical properties of FRP exhibit a significantly reduction at , which determines the limited use temperature for FRP composites [1]. The matrix transfers from compact viscous polymer to porous solid char and volatiles over the decomposition temperature () range. The mechanical properties of FRP composites dropped once again. Therefore, the usage of FRP composites in engineering can be hindered by the high degradation of mechanical performance at elevated and high temperatures.

In the past two decades, a large number of experimental studies have been conducted to evaluate the thermal and thermomechanical performances of FRP materials. In 1999, the fire resistance of GFRP sprinkler pipes with empty cavity, stagnant water, and flowing water was investigated by Davies and Dewhurst [2]. The failure time of the pipes with empty cavity and stagnant water was 1.5 min and 8.5 min, respectively, while the flowing water-cooled pipe remained structural integrity subjected to 120 min fire exposure. Inspired by the water-cooling concept of filament wound GFRP pipe used in [2], this concept, developed by Keller et al., was applied in the fire endurance of pultruded GFRP multicellular panels [3, 4]. The fire performance of GFRP panels was investigated in three different experimental parts: charring of GFRP laminates, fire endurance of liquid cooling moderately sized GFRP panels, and fire resistance of structural liquid cooling full-scale GFRP panels, respectively. In the first part, the pultruded GFRP laminates started burning at roughly 6 min, and the temperature of cold face reached at approximately 10 min [5]. In the second part, the temperature profiles through the thickness of lower face sheet were much lower than those of charring experiments [3]. In the third part, for the noncooled specimens, the cold face of lower face sheet reached and in 10 min and 57 min, respectively. The test results demonstrated that liquid cooling was an effective way to improve the fire resistance of pultruded GFRP components. Correia et al. [6] conducted structural fire endurance experiments on full-scale cellular GFRP columns subjected to one-side fire exposure. The noncooled column failed at 49 min due to global buckling, while structural function of the water-cooled columns could be maintained for two hours. In 2010, the fire resistances of noncooled and water-cooled pultruded GFRP beams with square hollow section were investigated by Bai et al. [7]. The bottom flange of the GFRP beams was subjected to ISO-834 fire. The fire resistance was over 120 min by using the water-cooled approach with a 72 mm/s flowing rate. But the noncooled specimen failed abruptly after about 38 min due to the kinking and buckling of top flange. In 2015, the fire performance of pultruded GFRP columns was investigated by Morgado et al. [8]. The columns with passive (calcium silicate boards) and active (water-cooling) fire protection systems were exposed to one-side (bottom flange) and three-side (bottom flange and two webs) fire, respectively. The test results showed that the efficiency of fire protection systems depended on the fire exposure sides. For one-side fire exposure, the water-cooling with a 72 mm/s flowing rate exhibited the best fire protection, which provided more than 120 min of fire resistance, while for three-side fire exposure, the column with calcium silicate boards showed the most effective protection, which provided roughly 40 min of fire resistance.

The mathematical modeling of the thermal response of FRP composite under elevated and high temperatures can be traced back to the 1980s. In 1981, Griffis et al. [9] proposed a 1D finite difference model to simulate the temperature profiles of graphite epoxy coupons subjected to intense surface heating. The analytical thermal responses agreed well with the test results. A thermomechanical model was developed to predict the thermal response of graphite epoxy composites and wood [10]. The proposed model can reasonably predict the temperature profiles and strength degradations. A 1D transient thermal model based on chemical kinetics, developed by Henderson et al., was applied to predict the temperature responses of fiber-reinforced phenol-formaldehyde polymer composite [11]. Good agreement was found between the analytical and experimental temperature profiles. In 1995, a thermochemical model of Thick composite laminates subjected to hydrocarbon fire was presented by Gibson et al. [12]. The 1D finite difference method was adopted to solve the proposed model. In 1996, this model was improved by Looyeh et al. [13]. The finite element method was introduced to predict the thermal performance of polymer composite materials. In 2000, Dodds et al. [14] updated the model for simulating the composite laminates. A straightforward explicit finite difference format was applied to solve the model. In 2006, the apparent specific capacity, developed by Lattimer and Ouellette [15], was input into a heat transfer model to predict the temperature profile through E-glass/vinyl ester composite laminates. Gibson et al. [16, 17] investigated the postfire mechanical properties of polymer composites by combing the thermal model with Mouritz’s two-layer mechanical model [18]. In 2012, the thermal responses of FRP laminates subjected to three-point bending and one-side heat flux were experimental and numerical investigated by Gibson et al. [19]. More recently, Miano and Gibson [20, 21] introduced a simplified thermal model by using the apparent diffusivity (ATD) method. The temperature-dependent ATD can significantly simplify the computational procedures and improve the stability of the numerical solutions. Feih et al. [22, 23] presented the thermomechanical model to predict the tension and compression properties of FRP laminates in fire. A 1D model was further developed to simulate the thermomechanical responses of sandwich composites [24]. This model also validated and applied by Anjang et al. [25, 26] to investigate the in-fire and postfire mechanical properties of FRP sandwich composite structures with balsa wood core. In 2006, a thermochemical and thermomechanical models were introduced by Keller et al. [27]. The models were used to predict the structural response of water-cooled multicellular GFRP slabs. In 2007, Bai et al. [28] proposed the chemical kinetics-based thermophysical model by considering the decomposition of matrix resin. Based on the temperature-dependent thermophysical model, 1D thermal model with the effective specific capacity was developed [29]. The predicted temperature responses agreed well with the experimental results. In 2008, the mechanism-based models [30–32] were developed to investigate the mechanical behavior of FRP composites in fire. This model was further validated on the predicting of the time-dependent temperature responses, elastic modulus degradation, and time to failure of water-cooled full-scale GFRP cellular columns [6]. In 2011, Miano [21] developed a three-dimensional (3D) finite element (FE) thermal model with ATD to estimate the temperature profiles of carbon fiber-reinforced polymer (CFRP) wing box laminates. More recently, a 3D FE model was developed by Shi et al. [33] to investigate the coupled temperature-diffusion-deformation problem of silica/phenolic composite materials. The accuracy of this model was validated by comparing the measured temperatures and displacements with numerical results.

Among the various existing thermal models, the 1D finite difference method was most frequently used to predict the thermal response of FRP composites. However, the FRP structural components applied in building construction were subjected to multidimensional fire exposures. The temperature profile of heated zone of FRP composites cannot be predicted by the use of 1D model. Hence, to address this issue, a 2D thermal model was proposed in this study. The accuracy of the proposed model was validated by the existing experimental data of pultruded E-glass/polyester rectangular tube under four-point bending and fire from one side. The comparison indicated that the developed two-dimensional thermomechanical model can be reasonable to predict the thermomechanical responses of GFRP composites.

2. Mathematical Model

2.1. Modeling Assumptions

When subjected to fire, FRP composites undergo many complex thermal, physical, chemical, and structural failure processes [35]. The challenge to accurately modeling the temperature responses is to consider complex interaction of degradation processes. Previous studies [5, 28] showed that the thermal and mechanical processes were dominated by glass transition and thermal decomposition. The convective heat transfer of volatiles flow up through the decomposition front to the surface may have a small effect on the temperature profile in the FRP material [15, 21]. Hence, in this study, the 2D thermomechanical model was simplified by considering heat conduction, physical process (glass transition), chemical (decomposition), and mechanical degradations.

2.2. Thermophysical Properties Model

During fire processes, thermophysical properties include density, specific capacity, and thermal conductivity of virgin (nonchar) and char material.

A lot of works were made to investigate these properties. The temperature-dependent density can be described by using an Arrhenius equation based on chemical kinetics [11]:where is the density; represents the preexponential factor; is the activation energy; is the universal gas constant; and is temperature.

The change of density can be described by (1). But in the existing literatures, specific capacity and thermal conductivity of virgin and char were usually expressed by polynomial fitting instead of analytical from. Hence, the parameters used in fitting have no clear physical meaning. Bai et al. [28] developed a model for predicting temperature-dependent thermophysical properties. In this model, a conversion degree of decomposition was introduced to characterize the pyrolysis process of polymer resin, as indicated in where denotes the decomposition degree and , , and are kinetic parameters, which can be derived from thermogravimetric analysis.

Based on decomposition degree , the temperature-dependent density can be obtainedwhere is the density of virgin composite and is the density of char material.

The apparent specific capacity can be given as [28]where is the apparent specific capacity, which increases due to endothermic phenomenon during decomposition process [35]; and are apparent specific capacity before and after decomposition respectively; is decomposition heat; and is mass fraction of undecomposed FRP materials, which can be calculated by [28]where () is the initial (final) mass of FRP material.

In this paper the main research object is the cross section of rectangular GFRP tubes; only the cross-sectional properties should be considered. Previous study [36] showed that thermal conductivity parallel to the axis of fiber (0°) was much higher than in the transverse (90°) and through thickness direction, while the thermal conductivity of transverse and through thickness direction were similar for glass fiber/polyester composites. Hence, the transverse thermal conductivity can be simplified by using the through thickness thermal conductivity. The thermal conductivity through the thickness of FRP composites can be expressed as [28]where denotes the through thickness thermal conductivity; the subscripts and represent the thermal conductivity before and after decomposition, respectively; and and can be obtained by using the inverse rule of mixtures:where , , and are the thermal conductivities of fiber, matrix, and volatile gas, respectively; , , and are the volume fractions of fiber, matrix, and volatile gas, respectively.

2.3. 2D Heat Transfer Model

2D transient heat transfer equation can be expressed aswhere and are the thermal conductivity in and directions, respectively. The term on the left side of (9) refers to the rate of change of internal energy, and the other two terms on the right side denote the net heat flux.

The boundary conditions applied in thermal analysis can be divided into three types. The first one is the Dirichlet boundary condition, which can be expressed aswhere and represent the 2D spatial Descartes coordinates; subscript bd denotes the boundaries of rectangular tube; and is the time-dependent specific temperature on the boundary. The second one is the prescribed heat flux, which can be expressed aswhere is the time-dependent heat flux at the boundaries. The last one is the convection between the boundary and external environment, which can be expressed aswhere denotes the convective coefficients and represents the time-dependent temperature of external environment.

In this study, convection and radiation between the boundaries and external/internal environment can be obtained from

2.4. Thermomechanical Model

Young’s modulus, , can be expressed as [31]where is Young’s modulus in glassy state and is Young’s modulus in leathery and rubbery states; and is the conversion degree of glass transition and can be expressed by [31]where , , and are kinetic parameters of glass transition for FRP composites. Based on beam theory, the midspan deflection of GFRP rectangular beam under four-point bending can be calculated by where denotes the time-dependent midspan deflection; is half of the total load; is the distance between the support and one load; is the span; and is the second moment of area of the section. The total bending stiffness of the GFRP beam, , can be calculated as the sum of the stiffness contribution of each elements [31]: where the subscripts w, tf, and denote the web, top flange, and bottom flange, respectively.

3. Solutions of Governing Equation

3.1. ADI Techniques

Equation (8) is a 2D nonlinear parabolic partial differential equation, which can be solved directly by explicit difference scheme. However, the numerical stability of this scheme is rather poor when solving the governing equation with time-dependent thermophysical properties. When heat transfer switches from 1D to 2D problem, the coefficient matrix of governing equation transfers from tridiagonal to nontridiagonal matrix by using the traditional implicit difference scheme. The formats of results are in a very complicated set of algebraic equations in two dimensions, which is difficult to solve [37]. Alternating direction implicit (ADI) format was an effective method for solving multidimensional heat transfer problems, which was first proposed by Peaceman and Rachford [38] in 1955 and developed by Douglas and Rachford [39]. By using ADI method, the 2D equation can be splitting into the same 1D implicit format with simpler tridiagonal matrix algorithm at each time step. In this study, ADI method was introduced and applied in this simulation of 2D thermal and mechanical responses of GFRP rectangular profiles.

3.2. Mesh

Figure 1 presents the “C type” profiles by considering the symmetry of shape and boundary conditions of a rectangular GFRP tube. The lower-left corner was set as the coordinate origin. In general, the solution domain was meshed into elements, where is the number of elements in direction. represents number of elements through the thickness. denotes the total width (height) of the GFRP tube; then the space interval can be written as or . Similarly, the time step was used for simulating program, where t represents the total fire exposure time, and is the number of time step.

Figure 1 

Thermal boundary conditions of GFRP rectangular tube.

3.3. Finite Difference Scheme

The ADI difference method splits the time step from to into two substeps: Step I: and Step II: . In Step I, implicit scheme was used in direction while explicit scheme was applied in direction. In Step II, implicit scheme was then applied in direction while explicit scheme was applied in direction. For example, the difference schemes of all notes inside the domain were derived as follows:where the subscripts , denote the coordinates; , represent the space interval in and direction, respectively. As indicated in (17), based on thermophysical properties and temperatures of the last time step , the temperature in time step can be calculated by using of Gauss-Seidel iterative approach.

The above only presents the notes inside the calculation domain without any external heat source. Therefore, the differential equations related to various boundary conditions must be derived. There are a total of 16 types (8 points and 8 sides) of boundary conditions in this model. The detailed expressions of ADI difference format were given in Appendices A~O, and the corresponding schematic diagram was given in Figure 2.

Figure 2 

Schematic diagram of thermal boundary conditions.

The difference format of decomposition degree can be expressed aswhere is the current time step and denotes the previous time step.

Based on the decomposition degree at the last time step , the difference schemes of thermophysical properties at time step can be derived as

The difference schemes of temperature-dependent Young’s modulus at time step can be expressed as

The schemes of the total bending stiffness at time step can be expressed as where and were the number of layers of the webs and flanges, respectively. Therefore the midspan at each time step can be calculated:

4. Validation and Discussion

4.1. Basic Model

As described in Section 1, Bai et al. [7] conducted the fire structural experiment of pultruded E-glass/polyester composite square beam subjected to ISO-834 fire, and a 37 min fire resistance was reached. In this paper, this experiment is used to validate the 2D numerical model developed in Section 3.

For boundary conditions, the hot face of bottom flange was set as prescribed ISO-834 curve because the temperature progressions of hot face were very close to ISO-834 fire curve. The outer faces of webs were set as adiabatic condition. As shown in Figure 1, convection and radiation play an important role in the rest boundaries, and the heat convective coefficients at the cold face of bottom flange were estimated as 2 W/m²K by fitting the temperature profiles at the cold face of the bottom flange [34]. The convective coefficients at inner faces of web and top flange were set as 0 by minimalizing the difference between modeling temperatures and test results for webs and top flanges. This indicated that convection can be omitted. In fact, the temperature of air in the inner cavity was higher than ambient temperature (20°C), and the convection was weak due to the close temperature for both webs and hot gas in the cavity. The effects of heat radiation, which was described by Stephan-Boltzmann law, were considered for all inner faces and out face of top flanges. The Stephan-Boltzmann constant was given in Table 1 and the solid emissivity can be seen from (23), which was adopted from [27] for modeling thermal response of a very similar E-glass fiber-reinforced polyester composite. The convection coefficient of at the outer face of upper flange can be obtained based on hydrodynamics proposed by Tracy (see (24)). The detailed modeling parameters were summarized in Table 1.where is the thermal conductivity of air (0.03 W/m K); is the acceleration of gravity (9.81 m/s²); is the Prandtl number (0.722); denotes the volumetric coefficient of thermal expansion of air (3.43 × 10^-3K⁻¹); and represents the kinetic viscosity of air (1.57 × 10^-5m²/s) [27].

4.2. Validation and Comparison

Figure 3 shows the comparisons of temperature progressions at different depths from the hot face of bottom flange. During the first 5 min, the heating rate of oven was much higher than the ISO-834 curve, which was used as the temperature of hot face, and this led to an underestimation of temperatures from 5 min to 20 min. This may attribute to the temperature lag resulting from the heat transfer from the oven to GFRP beam. In addition, the forming charred layers at the bottom flange may act as an insulator which producing a high experimental temperature due to the accumulation of hot gas in the pyrolysis zone. This insulation effect can be implemented in the model by adding a convection term of pyrolysis gas in (8) and a momentum equation (using Darcy’s law) to simulate the gas transport in the porous charred layers. A complex computational fluid dynamics (CFD) model should be developed and solved. Nevertheless, for the duration from 0 min to 5 min and 20 min to 40 min, the modeling results agreed well with the test values. The temperature gradients at 5 min, 10 min, 20 min, and 40 min were presented in Figure 4. A large temperature gradient through the thickness of bottom flange was found at 5 min and 10 min; however, the gradient decreased from 10 min to 40 min. This indicated that a more even temperature distribution was found through the thickness of bottom flange, and the temperatures of hot face and inner face were approximately 750°C and 900°C at 40 min (Figure 3). The elastic modulus profiles at 5 min, 10 min, 20 min, and 40 min were given in Figure 5. At 5 min, the modulus was 5.8 GPa and the glass transition was finished. Besides, the modulus of web also was also degraded. At 40 min, the modulus of bottom flange was almost lost and the bottom flange was fully decomposed with the lower part of web was in rubbery state.

Figure 3 

Comparisons of temperature progressions of bottom flange.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 4 

Temperature distribution of GFRP tube: (a) 5 min, (b) 10 min, (c) 20 min, and (d) 40 min.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 5 

Elastic modulus distribution of GFRP tube: (a) 5 min, (b) 10 min, (c) 20 min, and (d) 40 min.

The time-dependent elastic modulus profiles of web can be seen from Figure 6(a), and the decreasing modulus of elements resulted in an upper migration of neutral axis (Figure 6(b)). In the fire duration, the distance from the neutral axis at hot face was increasing from 50 mm to 75.5 mm while the one at outer face of top flange was decreasing from 50 mm to 36.5 mm. As a result, the time-dependent second moment of area of each element and bending stiffness was shown in Figures 6(c) and 6(d). A 78% reduction of bending stiffness was found during the whole fire durations. Figure 7 shows comparison of the midspan deflection progressions between the 2D modeling results, test values and 1D modeling results from [34]. As can be seen from Figure 7, the 1D thermomechanical model developed by the Bai et al. can accurately predict the midspan deflections during the first 12 min. However, due to the limitation of 1D model, the thermal responses of webs cannot be predicted; therefore mechanical responses after the onset of web degradation were not presented. The 2D modeling results agreed well with the test values in the first 8 min; however, overestimation was found from 8 min to 37 min. This may be attributed to two main reasons. First, it was assumed that the total bending stiffness was calculated as the sum of independent individual element. However, in fact each element in the cross section was constrained by the adjacent elements. Therefore this led to underestimation of total bending stiffness. In addition, here it was assumed that the whole span was subjected to the same loading and exhibited the same stiffness degradation, but only 0.95 m (64% of the span of 1.51 m) was directly exposed to the oven [34]. Therefore the overall bending stiffness was underestimated and the deflections were overestimated.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 6 

Mechanical responses of GFRP tube: (a) elastic modulus, (b) neutral axis, (c) second moment of area, and (d) bending stiffness.

Figure 7 

Comparison of mechanical responses of GFRP tube.

5. Conclusions

In this study, a 2D model to predict the thermomechanical responses of rectangular GFRP tubes were developed and solved by ADI finite difference scheme. The modeling results were compared with the temperature and mechanical responses from experiment on GFRP rectangular beam under the four-point bending and one-side fire exposure. The main conclusion remarks were drawn:(1)A 2D thermomechanical model based on finite difference method was developed and solved by Gauss-Seidel iterative approach. The detailed ADI schemes considering complex boundary conditions were also presented, and the numerical results were stable due to the stability of this difference format.(2)The proposed 2D thermal model can reasonably predict the temperatures of GFRP rectangular tube. Based on the calculated temperature, the temperature-dependent conversion degree of glass transition and decomposition as well as thermophysical properties of the whole cross section can also be obtained.(3)The evolution of elastic modulus, neutral axis, second moment of area, and bending stiffness can also be captured by the thermomechanical model, and midspan deflections of GFRP beam under four-point bending and fire exposure from one side can be reasonably predicted.(4)The overestimation of the midspan deflections may be resulted from the underestimation of bending stiffness and this may be attributed to two reasons: (a) the individual element in the cross section, in fact, was enhanced by the constrain of the adjacent elements; (b) it was assumed that the whole span was subjected to the same loading and exhibited the same stiffness degradation, but only part of the span was exposed to the oven directly [34], and this also make a contribution of underestimation of stiffness.(5)The 2D thermomechanical model can also be extended to simulate the thermomechanical responses of beam subjected to four-point bending and three-side fire loading. The 2D thermal also could be used to model the thermal responses of GFRP tubular columns subjected to multiside fire exposure, and this makes it possible to provide a base for estimating the mechanical performance of the GFRP tubular column under mechanical and multiside fire loading.

Appendix

A.

B.

C.

D.