Numerical simulation of the experimental behavior of RC beams at elevated temperatures

Introduction Reinforced concrete (RC) is globally the predominant construction material of choice as structural components. Structures during their service life in addition to dead, live, wind and earthquake loads, may encounter elevated temperatures resulted from fire incidents. Therefore, resistance of RC members against fire is becoming an important parameter in recent design code such as ASCE [1], Eurocode2 [2], BS 8110–2 [3], ACI 216.1 [4] and AS 3600 [5]. The overall behavior of RC exposed to fire is dependent on mechanical and thermal properties of concrete and steel that fluctuate as a function of temperature. Concrete as a material can resist elevated temperatures resulted from fire events [6] and does provide some sort of thermal insulation to protect steel bar due to low conductivity, slower strength degradation at elevated temperatures and high heat capacity [4–6]. Fire safety of structural members are determined in term of duration with respect to structural stability and integrity [7, 8] which is the behavior of RC members exposed to temperature–time curves per standard fire tests as defined by ASTM E119 [9], EN, 1992–1-2 [2], ISO 834–1 [10], BS 476–20 [11]. Most current practices are more of a visual descriptive approach that is more concern with the adequate steel bar minimum cover and member size while ignoring a more meticulous approach of understanding the thermal and mechanical behavior of RC members exposed to elevated temperatures [12]. Abstract


Experimental program
In order to gain an understanding of the thermal effects on reinforced concrete structural elements, a study was conducted to experimentally evaluate the durability of concrete cylinders exposed to thermal conditions. Thus, the following objectives were perused: -To examine the response of concrete cylinders before and after thermal exposure due to dry heat. -To determine the strength reduction factors related to various thermal exposure levels. -To identify the different possible concrete failure modes as related to thermal exposure levels.
For the purpose of this experimental study, two batches of fifteen and twelve standard cylinders (150 by 300 mm) were of casted natural aggregates from coarse limestone aggregates and quartz sand. The specimens of the first batch (B1) and the second batch (B2) were divided respectively into five and four series of dry-heat and lab ambient temperature series. All irregularities found on the concrete surface of the cylinders were removed and the cylinder surfaces were well cleaned. Three specimens were tested at each temperature level with the highest and lowest temperature limits set as 260 °C and 25 °C, respectively. The highest temperature was limited to 260 °C due to physical and practical constraints. The specimens were placed and maintained in an oven at set temperatures for a period of 24 h prior to testing, thus guaranteeing that the cylinders are preserving the overall constant temperature while being tested. The strength of the test specimens was obtained according to ASTM C39 [7]. Instrumentation and testing procedure for all cylinders were kept the same for all types of thermal exposures. The tested specimens were removed from exposure and moved quickly to the UTM machine where the compression test was conducted in order to avoid as little as possible of thermal differential. A compressometer-extensometer was mounted on each concrete cylinder prior to placing it in the UTM, where it was positioned between 2 removable neoprene caps, protected with cardboard in order to shield the neoprene the heated specimens. Compression test was conducted on the specimens in accordance with ASTM C39 [25]. The data acquisition system was setup to take readings from the displacement sensors every 5 s in synchrony with load data acquisition and collection process.

Experimental results
The compressive strength experimental results are displayed in Table 1. The results clearly indicate the strength of concrete had decreased with the increase of temperature. Stress strain figures were generated using the data collected from the automatic data acquisition system incorporating the corresponding collected loads values. Figure 1 represents a single graph of the results of all Batch 1 tested specimens combined. Every curve represents the average result of three samples at the same temperature. Figure 2a shows the average of ultimate strength clearly displaying that the ultimate concrete strength decreases with the increase of temperature exposure. Reduction of compressive strength of normal strength concrete at elevated temperature was confirmed by other references such as Eurocode [2] and ASCE [1]. Figures 2b illustrates the reduction factor as obtained from test data and comparison of them to relations that proposed by Eurocode [2] and ASCE [1]. According to ASCE [1], there is no variation in the reduction factor in the range of 20 °C to 450 °C, but at higher temperatures, there is a decrease. In the Eurocode [2], the reduction factor is different in term of concrete material. Although concrete stress-strain curve is not linear, but it is assumed as being nearly linear for a compressive strength lower than forty percent of ultimate compressive strength [2,26]. Thus, the slope of the curve is the modulus of elasticity. In this study, a linear behavior for concrete up to a strength of about 100 MPa was assumed. Figure 3 presents the modulus of elasticity of the specimens at various temperatures.

Numerical simulation
The performance of RC structures at elevated temperatures is very intricate. As temperature increases, mechanical, physical, and thermal properties of concrete and steel reinforcement considerably change. Empirical and computational procedures [12,[27][28][29] have been comprehensively utilized in this study of RC structures exposed to elevated temperatures or accidental fire. While investigating the effects of elevated temperatures, buildings codes [30] usually do provide a variation of equation for the calculation of fire ratings of structures. Nevertheless, these equations do not incorporate  the thermo-mechanical performance of a structure during fire [28]. Meanwhile, highly sophisticated non-linear mechanical numerical 2D or 3D finite element method [27][28][29] were being developed incorporating all possible effects of elevated temperatures on RC. However, these methods were computationally very intensive and did not accurately reflect the constitutive laws incorporating the effects of elevated temperature. Empirical and numerical procedures are widely utilized for the estimation of the capacity of a RC structure to resist fire, which are limited to simple cases. Finite element methods based concrete damaged plasticity models have been employed to perform deformation analysis of bending RC members. The models guarantee higher accuracy of deflection predictions in comparison to the design code methods. The model incorporates efficiently a combination of accuracy and simplicity through the utilization of the layer section model. In this study, an effort has been made to extend application of the layer section model to strain analysis of RC bending members exposed to elevated temperatures incorporating non-linear physical and thermo-mechanical materials properties. An advanced finite element model was developed to evaluate the stress-strain state, load bearing capacity and failure time. The material properties of both concrete and reinforcing steel bars critically affects the mechanical and thermal behavior of RC elements. Numerous studies are available to comprehend the mechanical and thermal properties of the reinforcing steel [1,2,31,32]. Meanwhile, there are three concrete material models which are commonly being used in the numerical simulation for concrete damage: (1) Smeared crack concrete model, (2) Brittle crack concrete model, and (3) Concrete damaged plasticity model. In this study, the third model was selected because it can represent the inelastic behavior of concrete both in tension and in compression including damage behavior. In this model, it is assumed that "tensile cracking" and "compressive crushing" are the two main mechanism of concrete failure. In addition, this model has stiffness recovery ability in cyclic loading i.e., changing tensile stress to compressive stress at an element or reverse. The failure criterion is not defined in this model, so cracking or eliminate of an element are not possible. Nevertheless, concrete damage plasticity model can show the location of cracks and their direction. In addition, to avoid extensive failure, it is recommended to use adaptive meshing [33].
In fact, concrete damage plasticity failure surface is a modification of the Drucker-Prager hypothesis [34] that assumes that failure is determined by non-dilatational strain energy and the failure surface is a cone shape in the stress space. The ratio of the distances between the hydrostatic axis and respectively the compression meridian and the tension meridian in the deviatoric cross section is referred to as K c with a value higher than 0.5 and can increase to 1. According to experimental results from triaxial-stress tests, the K c value for a mean normal zero stress amounts to 0.6 and slowly increases with a decreasing stress and it is recommended to assume K c = 0.667 in concrete damage plasticity model [34].
In this model, the ratio of tensile strength to compressive strength is introduced by plastic potential eccentricity (ε). Experimental results indicate that the plane's meridians in the stress space are almost hyperbola. The parameter ε is a representation of the rate of approach of the plastic potential hyperbola to its asymptote. In the Drucker-Prager hypothesis, the meridian plane assumes a straight line so the ε becomes zero but Page 6 of 17 Issa and Izadifard Adv. Model. and Simul. in Eng. Sci.
(2021) 8:12 in concrete damage plasticity model it is recommended to use a value of 0.1 for ε. The concrete internal friction angle is an important parameter that effects the material dilation when the concrete is under compound stress. Theoretically, this parameter (ψ) indicates the angle of inclination of the failure surface towards the hydrostatic axis. Many of researchers have recommended assuming ψ = 36° to ψ = 40° for concrete numerical simulation [32,34]. The ratio of the strength in the biaxial state compared to the strength in the uniaxial state (f b0 /f c0 ) should be defined in this model and is normally assigned a value of 1.16. The compressive stress-strain curve of a concrete could be defined based on uniaxial compression test. The initial part of this curve σ c is less than 0.5σ cu and for that portion, the curve is almost linear. The initial tangential modulus of elasticity E is obtained exactly by experimental data or approximately by one of the equations in Table 2.
The strain analogous to the ultimate compressive strength (σ cu ) can be exactly obtained by experimental tests or approximately by [35]: The descending portion of the curve ends at the ultimate compressive strain of concrete (ε cu ) that is given by [36]: Hsu and Hsu [35] used the following equations to calculate the compressive stress valve (σ c ) between the 0.5 σ u in ascending portion and the 0.3 σ u in the descending portion of stress-strain curve: The parameter β is given by: In the absence of the compressive load, a part of the strain remains in the elements as a residual strain. The tangent of return line in stress-strain curve (Fig. 4) depends on the compressive damage parameter (0 < d c < 1) that is defined as the ratio between the inelastic strain and total strain. In undamaged concrete and fully damaged concrete d c assumes 1 and 0 respectively. The following relation is used to evaluate at d c : According to Fig. 4, the total compressive strain (ε c ), elastic strain corresponding to the undamaged concrete (ε el 0c ), inelastic strain (ε in c ) and plastic strain (ε pl c ) are given by: The plastic strain values are neither negative nor decreasing with increasing compressive stress [32]. Considered for this study σ cu = 24.6 MPa, E 0 = 25,700 MPa (Table 3), β = 2.289 (Eq. 4), ε 0 = 0.0017 (experimental result). The utilized calculated values of other compressive parameters using Eqs. 2 to 6 and are presented in Table 3.
The maximum tensile strength of concrete could be determined by a direct tension test or cylinder splitting. Due to the absence of test results, the approximate relationship provided ACI-318 [14] was employed: Fig. 4 Modified tension stiffening model [36]  ABAQUS Manual [32] recommends the use of a tension-stiffening and strain-softening model. In this model, the total tension strain (ε t ), elastic strain analogous to the undamaged concrete (ε el 0t ), cracking strain (ε ck t ) and plastic strain (ε pl t ) are given by: where d t is the tensile damage parameter and it is zero for undamaged concrete and one for fully damaged concrete. Nayal and Rasheed [36] as presented in Fig. 4 have modified this model. In this new model, the stress-strain curve is divided in four parts: (1) the linear elastic relation from (0, 0) to (ε cr , σ t0 ), (2) the sudden drop at ε cr from σ t0 to 0.8 σ t0 , (3) primary cracking stage from (ε cr , 0.8 σ t0 ) to (4 ε cr , 0.45σ t0 ) and (4) secondary cracking stage from (4 ε cr , 0.45σ t0 ) to (10 ε cr , 0).
To avoid numerical errors, it is necessary to change the drop part to a sharp slope and where the curve ends to a stress more than zero. In this study, values considered were: E 0 = 25,700 MPa (Table 4), σ t0 = 2.73 MPa (Eq. 7) and ε cr = σ t0 / E 0 = 0.000106. Values of other parameters for concrete tension stiffening are calculated using Eq. 8 and presented in Table 4.
There are numerous studies concerning the values of Poisson's ratio for various types of concrete exposed to diverse situations. Marechal [37] conducted a few experimental tests to determine the variation of the modulus of elasticity and Poisson's ratio with temperature. Also Elghazouli and Izzuddin [38] proposed a model to determine the values. Based on those references [37,38], the Poisson's ratio of concrete was taken as a constant value of 0.20 from 20 °C to 150 °C. Beyond the 150 °C temperature, the Poisson's ratio was assumed to decrease bilinearly to 0.1 and with a zero value from 400 °C to 1200 °C.

Concrete cylinders
The experimental concrete cylinders of 150 mm in diameter and 300 mm in height were modeled using one eightf of the specimen due to symmentry in loading and geoemetry. Because of symmetry, a quarter of a cylinder with 150 mm length was modeled in ABAQUS/Explicit using C3D8R elements [32]. Abaqus/Explicit is a finite element analysis product that is particularly well-suited to simulate brief transient dynamic events. The C3D8R element is a general-purpose linear brick element, with reduced integration (1 integration point). The element has eight nodes and three translational degrees of freedom at each node. The model boundary conditions were defined complying with symmetry conditions. In the model, either force or displacement control methods could be employed to impose pressure on the samples. In the force control method, a uniform increasing pressure is imposed on the top face of specimens and in the displacement control method the top face is gradually load until failure occurs. Although both methods were utilized, but in order to more accurately simulate the experimental conditions, the displacement control method was employed in this study. The results of the numerical analysis was recorded at various sub steps and were plotted in the form of a stress-strain graph. In the numerical analysis, the element size is very important in determining the accuracy and computational time required. The smaller mesh size results in an increase in the number of elements and thus leads to an increase accuracy of results and computaional time. A mesh size effect study is displayed in Fig. 5-a and b. Thus it was concluded that if every radius of the cylinder was divided into eight partsthat would produce accurate results within reasonable computational time. A finite element code was implemented in the ABAQUS software to compute axial and diametrical strain versus axial stress. The material model of concrete at elevated (2021) 8:12 temperatures, the geometry of model, meshing and boundary conditions, the top face of model pull down gradually was employed. The cyclinderical diameterical elongation and the reaction force were computed at several height intervals. Figure 5-c represents the results of all models at different temperatures combined in one graph. It shows the ultimate compressive strength and the ultimate axial and diametrical strain of concrete decrease with temperature increasing resulting in the decrease of concrete ductility.

Reinforcing steel bars
The other elementary component was the longitudinal reinforcing bars having the main task of transferring the normal forces and thus, were modeled as a three-dimensional truss element (T3D2 [32]). Hence, eliminating the necessity of modeling the strain hardening and the elasto-plastic phenomena. The steel mechanical properties included a modulus of elasticity of 208 GPa and yield stress of 498 MPa in tension and 459 MPa in compression. In order to model the steel bars accounting for elevated temperatures, fire curves are required. Conflagration in a structure causes temperature increasing of compartment. Thus, time-temperature curve of Eurocodes [2], ISO [10], ASTM [9], or BS [11] could be used to determine the temperature profile. In Eurocodes, generally two types of fire curves are obtainable, namely a nominal curve and a parametric curve. The nominal curves just include the ascending heating phase, but the parametric curves include both heating and cooling phases. Fire proceeding continuity depends on the building utilization, type and number of combustible materials, internal fire extinguishment systems and distance of fire stations. If the fire extinguishing process is to be short in time, decreasing in the structural element's strength is not expected, but if it is for a long time, one can expect the beams and columns strength to diminish. Since the main purpose of this study is obtain the behavior of reinforced concrete at elevated temperatures, the heating phase is the focus of this study. Hence, the standard fire curve was utilized where the temperature increases as a logarithm function (Eq. 9) and the cooling phase is ignored.
where t and θ g are time (min) and temperature of fire compartment ( o C), respectively. An increase in temperature results in the decrease in stiffness and strength of the steel bars. The curves displaying the decrease in young's modulus, proportional limit and yielding strength of steel are presented in Fig. 6. Steel properties at elevated temperatures were obtaining by multiplying the steel properties at room temperatures by the reduction factors illustrated in Fig. 7.
In the present finite element model, the density of steel was set at 7850 kg/m 3 . Thermal conductivity and specific heat capacity of steel are temperature dependent, and they are utilized as specified in EN 1993-1-2 [39]. Also based on EN 1992-1-2 [2], the simulation of the constitutive model for steel at elevated temperatures included a total strain of two parts: the free thermal strain ε th and the stress-induced strain ε σ. Reinforced concrete RC beams were modeled utilizing the ABAQUS/Explicit routine employing the C3D8R hexahedral elements for concrete and T3D2 truss elements for the reinforcing steel bars. Bond interaction between the reinforcing steel and the surrounding concrete makes the tensile behavior of cracked concrete more complicated. Effect of elevated temperature on bond between reinforcing steel and concrete is not significant [40,41] and it was modeled accordingly. Seabold [42] conducted an experimental study on fifteen reinforced concrete beams. Eight beams were subjected to blast loading in a shock chamber and seven beams were exposed to a uniform distributed static pressure applied on the top face. The dimensions of the beams were 7.75 inches in width, 15 inches in depth and a length of 154 inches. Two longitudinal steel bars at top and bottom faces reinforced the beam. There were 2 inches center-to-center vertical stirrups at the right beam region, but the spacing of the stirrups increased at the opposite region to realize shear cracks and damage to this region. One of the beams, designated as WE7, was selected due to the identical concrete properties as the tested concrete cylinders. Figure 7-a and b display the RC cross sectional details and materials properties of the WE7 beam exposed to static loading [41]. In the numerical model a distributed uniform increasing pressure was imposed gradually on the top face of the beam until it was damaged and destroyed. The pattern of shear and bending cracking in the beam under static load is clearly visible in Fig. 7-c. Figure 7-d displays a comparison of the midspan deflection confirming that the utilized RC model will lead to accurate and reliable results.
There are three modes of heat transfer: convection, radiation and conduction. In the fire compartment, heat fluxes flow to the outermost surfaces of the beam through convection and radiation, whereas within the concrete body, heat transfer occurs by conduction. Two fire scenarios were considered: 1) the beam engulfed in fire (all faces of beam exposed to fire), 2) three sides of beam (under the cell) exposed to fire and the top face of the beam (on the upper floor) exposed to room temperature. The effects on the deflection of the beam due to elevated temperatures were investigated at three time intervals of 626, 1248 and 1864s posteriori to the start of the fire incident. Tables 5 and 6 display the values of temperatures for the two fire scenarios at various locations of the beam at the indicated time.
In the model, the sections of the meshed beam were assigned temperature related properties at the mesh element position. Furthermore, the properties of steel bars were modified according to the temperature profile. Figure 9 presents the deflection of mid span at various time intervals after the initiation of the fire event and clearly indicates that the beam exposed to fire becomes more flexible with the progress of time and thus displaying an increase in deflection with the progression of heat towards the center of its beam cross section. For a given strength capacity, as the duration of fire exposure increases, the beam is expected to fail at a given constant loading condition. In fact, when the rotation of two supports at the ends of beam reaches a value of 0.04 rad, the bearing capacity of beam would have reached the ultimate limit and thus resulting in failure. In the other words, when the mid span deflection reaches two percent of the beam length, the beam would fail. This issue is coherent with the results illustrated in Fig. 9. The selected beam could bear about 85 Psi pressure at room temperature. When the beam was engulfed with fire,  In the case where three sides of the beam were exposed to fire, the ultimate pressure decreased to 84.6, 84 and 82.4. Figure 10 displays the reduction of the beam load carrying capacity in terms of fire duration.

Conclusions
In this paper, the experimental results of twenty-seven concrete specimens at various constant temperatures are explored and numerically simulated utilizing an accurate three-dimensional finite element model of concrete properties at room and elevated temperatures. The results of the numerical simulation are compared with experimental data and incorporated in the mechanical and thermal analysis of RC beams exposed to fire. Heat transfer in RC beams is analyzed exposed to two scenarios, first the beam fully engulfed with fire and then only three sides of the beam exposed to fire. Determined are the temperature profiles at different sections of the beam and the overall mechanical response of the beam. The effect of the three-dimensional non-uniformity on the concrete material and as well as on the structure was not known and it could cause additional thermal stresses and could lead to the occurrence of spalling. Another shortcoming in the material modeling was that the concrete specimens were heated unloaded. Studies however showed that both compressive strength and elastic modulus reduce far less with increasing temperature when the concrete was loaded before and during heating [43].
In the constitutive models, it was determined that the structural analysis of heated concrete must include transient strain, or its deviations transient creep strain or load induced thermal strain. Nevertheless, the concrete material model presented in the Eurocode ignores these two effects. However, in some cases, the calculated failure time according to these models were verified to predict experimental results fairly well. This is attributed to the fact that structural elements are typically exposed to fire from underneath. The heat will thus only affect the tension zones of the members whereas transient strain solely acts in compression. Furthermore, very large strains on the descending branch were incorporated to compensate for the transient strain [44]. The strains ε c1 and ε cu are much larger than their respective published literature values. Nonetheless, this study presents a reliable numerical model that can accurately predict the stiffness and strength degradation of RC beam at elevated temperatures. Thus, the goal to develop a numerical model that reflects the true behavior of RC structures when exposed to fire, both from a material and structural point of view was realized. This model could be utilized to investigate and predict the behavior of RC structures and thus provide a valuable tool to induce changes in new RC fire design codes and methods.