In the present investigation, a temperature-dependent constitutive model for concrete damage plasticity is proposed using the concept of continuum thermodynamics. The essence of the constitutive model is in the definition of the relation between the mechanical, temperature, and damage fields using the Helmholtz free energy function. The Helmholtz free energy function is written as Equation 1: ψ ε e , α , d t , d c , T = ψ e ε e , d t , d c , T + ψ p α , T (1) where, ψ e represents elastic free energy, ψ p represents plastic hardening energy, ε e is the elastic strain tensor, α is the hardening variable, d t and d c are tensile and compressive damage variables, respectively, and T is temperature. The elastic free energy adopts a quadratic form. Based on the Clausius-Duhem inequality, the stress tensor expression can be derived as Equation 2: σ = ∂ ψ ∂ ε e = 1 − d C 0 T : ε e (2) where the damage variable d is determined by combining d t and d c according to the stress state, C 0 T is the temperature-dependent initial elastic stiffness tensor, ensuring the model satisfies the positive dissipation requirement of the second law of thermodynamics.

The proposed model differs from the conventional CDP model in ABAQUS in three key aspects: (1) thermodynamic consistency through Helmholtz free energy formulation ensuring positive dissipation at all temperatures; (2) temperature-dependent damage evolution with coupled thermal-mechanical degradation, whereas conventional CDP treats damage parameters as temperature-independent; (3) explicit thermal strain decomposition with temperature-varying expansion coefficients. Conventional CDP requires manual input of discrete stress-strain curves at each temperature without ensuring thermodynamic consistency between states.

The strain rate tensor takes an additive decomposition expression as Equation 3: ε · = ε · e + ε · p + ε · t h (3) Where, ε · e , ε · p and ε · t h represent elastic, plastic, and thermal strain rates, respectively. Thermal strain adopts an isotropic thermal expansion assumption as Equation 4: ε t h = α T T − T 0 I (4) where α T is the temperature-dependent linear expansion coefficient, T 0 is the reference temperature (20 °C), and I is the second-order identity tensor. The evolution of plastic strain follows a non-associated flow rule as Equation 5: ε · p = λ ˙ ∂ G ∂ σ (5) where, λ ˙ is the plastic multiplier and G σ , T is the plastic potential function, adopting a Drucker-Prager type dilation form as Equation 6: G = ϵ σ t 0 T tan ⁡ ψ 2 + q 2 − p ⁡ tan ⁡ ψ (6) where, q = 3 J 2 is the Mises equivalent stress, p = − tr σ / 3 is the hydrostatic pressure, ϵ is the eccentricity parameter, ψ is the dilation angle, and σ t 0 T is the temperature-dependent uniaxial tensile strength.

The yield criterion adopts a modified Drucker-Prager form as Equation 7: F = 1 1 − α q − 3 α p + β ε ∼ p l , T σ ^ max − γ − σ ^ max − σ c 0 T ≤ 0 (7) where, α and γ are material parameters, β is the hardening function related to temperature and equivalent plastic strain, σ ^ max is the maximum principal stress, · denotes the Macaulay bracket, and σ c 0 T is the temperature-dependent uniaxial compressive strength. The hardening variable is controlled through the equivalent plastic strain parameter ε ∼ p l , with its evolution rate as Equation 8: ε ∼ ˙ p l = λ ˙ 2 3 ∂ G ∂ σ : ∂ G ∂ σ (8)

Damage evolution is a key component of the model. The tensile damage variable d t and compressive damage variable d c describe stiffness degradation under tensile and compressive stress states, respectively. The damage evolution equations adopt an exponential softening form as Equation 9: d t = 1 − r t 0 T r t exp − 2 H t T G f T r t − r t 0 T d c = 1 − r c 0 T r c 1 − α c T r c − r c 0 T − r c 0 T exp 2 H c T r c − r c 0 T (9) where, r t and r c are damage driving forces in tensile and compressive states defined through plastic damage energy release rates; r t 0 T and r c 0 T are initial thresholds; H t T and H c T are softening moduli; G f T and G c T are temperature-dependent tensile and compressive fracture energies; and α c T is the compressive softening parameter. These parameters are all regularized through temperature reduction factors to ensure mesh objectivity of the damage process. Damage variables increase monotonically from zero (undamaged material) to unity (complete failure), and the damage process is irreversible.

The influence of temperature on concrete material properties is reflected through temperature dependence of material parameters. Compressive strength, tensile strength, and elastic modulus all adopt normalized temperature reduction coefficients as Equation 10: f c T = k c T · f c 20 ° C , f t T = k t T · f t 20 ° C , E c T = k E T · E c 20 ° C (10)

Temperature-dependent material properties are derived from EN 1992-1-2 recommendations (Standard, 2004) and validated against experimental data from Schneider (1988) and Phan and Carino (2002), obtained following RILEM TC 129-MHT testing protocols for heated concrete. Based on these data under ISO 834 standard fire conditions, this study establishes material performance degradation laws for normal-strength concrete in the range of 20 °C–800 °C. Compressive strength changes little below 200 °C, decreases rapidly between 200 °C and 600 °C, drops to 45% of ambient strength at 600 °C, and retains only about 20% at 800 °C. Tensile strength degrades more significantly, dropping to 40% of ambient value at 400 °C. Elastic modulus degradation trends similarly to strength. Detailed temperature-dependent material parameters are provided in Table 1.

For steel reinforcement, this study adopts a temperature-dependent elastoplastic constitutive model. The stress-strain relationship of HRB400 grade hot-rolled ribbed steel at elevated temperatures is characterized by a multilinear hardening model, including elastic stage, yield plateau, and strain hardening stage. The yield strength, ultimate strength, and elastic modulus of steel also adopt temperature reduction factors as Equation 11: f y T = k s y T · f y 20 ° C , f u T = k s u T · f u 20 ° C , E s T = k s E T · E s 20 ° C (11)

Steel strength degrades faster than concrete at elevated temperatures, with yield strength dropping to approximately 70% of ambient value at 400 °C, 40% at 600 °C, and only 10% at 800 °C. Meanwhile, steel ductility increases at elevated temperatures, with ultimate strain increasing from approximately 10% at ambient temperature to 15% at 600 °C. The thermal expansion coefficient of steel remains essentially constant in the range of 20 °C–800 °C, with a value of 1.4 × 10 − 5   / C ° . While creep and relaxation effects of steel under long-duration fire cannot be neglected, this study focuses on standard fire scenarios within 2 h, where time-dependent effects are relatively minor and therefore not explicitly considered in the model.

Through this series of temperature-dependent constitutive relationships, the model can accurately describe the stress-strain response, strength degradation, stiffness loss, and thermal expansion behavior of concrete and steel from ambient temperature to 800 °C at the material point level, providing a reliable material basis for structural-level fire performance analysis.

The temperature-dependent CDP model is integrated into the finite element software ABAQUS through the user material subroutine VUMAT, enabling thermal-mechanical coupled analysis of structures under fire conditions (Aliş et al., 2022). The strategy employed in the problem is a kind of sequential coupling, where the transient heat transfer calculation is done first to determine the temperature field distribution in the structure, followed by the structural calculation done using temperature history as an external field loading in the structure, updating properties accordingly at each time step. The temperature diffusion time is relatively long, but it is much shorter compared to the structural time, making it possible to use an approximation where the temperature field can be modeled using a unidirectional external loading on the mechanical system. The time step in heat transfer calculation is in the range 1–5 s for an accurate calculation since it is done for transient responses in temperature variations, while the time step in the structural calculation starts with 10 s since structural variation is a relatively slower process than temperature variation. Time step in structural calculation dynamically varies according to its convergence level with a maximum time step no more than 60 s and no less than 0.1 s.

Within the VUMAT routine, constitutive integration is performed using an incremental-iterative procedure whose flow is illustrated in Figure 1. The procedure begins by reading the total strain increment, the preceding stress state, internal variables, the present temperature, and the preceding temperature from the main program. After that, strain decomposition is done by determining the value of the thermal strain increment and then subtracting it from the total strain increment to get the mechanical strain increment. An elastic trial step is executed next, assuming the current increment is a purely elastic response and calculating the trial stress state using the temperature-updated undamaged stiffness. The algorithm then checks whether the trial stress satisfies the yield condition, with yield tolerance set to one millionth of the compressive strength at the current temperature. If the trial stress lies inside the yield surface, the step represents a purely elastic response and proceeds directly to damage update; if it exceeds the yield surface, plastic correction is required.

During the plastic correction phase, a back-mapping algorithm is employed to solve the consistency condition through Newton-Raphson iteration, bringing the stress back to the yield surface. The iterative process simultaneously updates the stress, plastic multiplier, plastic strain, and hardening variables, and a consistent tangent stiffness matrix is used to accelerate convergence. The iterative convergence criteria require the residual norm to be less than 1.0 × 10 − 8 and the relative error of the stress increment to be less than 1.0 × 10 − 6 . The maximum number of iterations is set to 50. If convergence fails, the time step is reduced and the calculation is repeated. After convergence, the plastic strain and hardening variables are updated to provide a basis for subsequent damage evolution calculations. The entire plastic correction process ensures that the updated stress state strictly lies on the yield surface at the current temperature, satisfying the plastic consistency condition.

Damage evolution calculates the damage driving force based on the updated elastic strain, tracking the historical maximum values of tensile and compressive damage, respectively. The tensile and compressive damage variables are updated according to the damage evolution equation, and the damage increment uses implicit backward Euler integration to ensure numerical stability. In order to maintain numerical stability during the softerening stage, a method for regularizing the fracture energy is considered. Fracture energy is considered to be the energy for total fracture work done over a unit area, which is a function that reduces with an increase in temperature. By associating fracture energy with the characteristic length of the element, the softening modulus is decoupled from the mesh size, avoiding mesh sensitivity issues. Furthermore, a viscosity regularization technique is introduced to stabilize the numerical response of the material in the post-peak softening region. The viscosity parameter is taken as one ten-thousandth of the time increment. This viscosity regularization introduces a small strain rate effect in the quasi-static analysis, effectively suppressing numerical oscillations during the softening process.

Finally, the actual stress is calculated, the effective elastic stress is reduced by the damage variable, and the consistent tangent stiffness matrix is updated and output to the main program for assembly of the global stiffness matrix. The model was verified at the material point level by simulating uniaxial tension and compression tests at different fixed temperatures. The simulation results are in full agreement with the stress-strain curves used during calibration, confirming the correctness of the constitutive integral. The model meets the thermodynamic consistency requirements, and the global energy dissipation is verified to be always positive by monitoring the energy balance, without generating spurious energy creation. Although an explicit version of VUMAT was also developed for potential dynamic fire scenarios, this paper mainly uses implicit static analysis, which is more suitable for quasi-static fire problems. The temperature-dependent CDP model is seamlessly integrated with finite element analysis. By updating the material stiffness and strength as a function of local temperature, the finite element simulation can capture the gradual weakening process of the concrete structure during the development of the fire.

To validate the model’s predictive capabilities at the structural level, a three-dimensional finite element model of a typical reinforced concrete simply supported beam was constructed. The beam dimensions for Cases 1–3 were 200 mm width, 400 mm height, and 3,600 mm span; for Cases 4–5, the geometry was adjusted to 200 × 350 × 3,200 mm following the respective experimental configurations (Table 2). Conventional reinforcement was employed: four 16 mm diameter longitudinal bars were arranged in the tension zone, two 12 mm diameter longitudinal bars were arranged in the compression zone, and stirrups were 8 mm diameter bars spaced 200 mm apart. The concrete strength grade was C30, with a standard axial compressive strength of 30 MPa and an elastic modulus of 30 GPa. HRB400 grade hot-rolled ribbed steel bars were used as the reinforcement, with a yield strength of 400 MPa and an elastic modulus of 200 GPa. The concrete cover thickness was 35 mm, meeting conventional fire protection design requirements. The detailed geometry of the finite element model is shown in Figure 2.

In the finite element model, the concrete was modeled using an eight-node hexahedral temperature-displacement coupling element, C3D8T, with each node representing temperature and three translational degrees of freedom. The element mesh size was kept around 25 mm, with refinement to 15 mm near the rebar and on the exposed surface to improve the accuracy of capturing temperature gradients. To verify the rationality of the mesh size, a mesh convergence analysis was performed using four mesh sizes: 35, 25, 15, and 10 mm. Results showed that when the mesh size was reduced from 25 to 15 mm, the change in mid-span deflection and fire resistance time was less than 3%, specifically, fire resistance times for mesh sizes of 35, 25, 15, and 10 mm are 122, 126, 127, and 127.5 min respectively, confirming mesh convergence. While the calculation time increased by approximately 150%. When further improved to 10 mm, it entails an improvement of less than 1% in performance, though the computational time is doubled. Taking into consideration factors such as accuracy and computational time, a mesh size of 25 mm for general regions with further meshing to 15 mm for critical regions is adopted. The reinforcement is modeled using embedded constraints assuming perfect bond, which is justified for ribbed bars (HRB400) when rebar temperatures remain below 400 °C–500 °C. For the baseline 35 mm cover, this temperature threshold is reached at approximately 90 min. Literature indicates bond degradation effects are less pronounced for ribbed bars than smooth bars, and the 2.9% prediction error against experimental data confirms acceptable accuracy for fire durations within 2 h. The rebar temperature is transferred from the surrounding concrete by conduction. Due to the rebar’s low heat capacity and high thermal conductivity, its temperature approaches that of the concrete at the same location.

The boundary conditions are set according to the mechanical characteristics of the simply supported beam: one end constrains three translational degrees of freedom to simulate a hinged support, and the other end constrains vertical and lateral displacements to allow longitudinal free sliding to simulate a roller support. The load is applied in the form of four-point bending, and concentrated forces are applied at each of the three equal points of the span to produce a pure bending area in the middle. The load ratio is defined as the ratio of the load acting during a fire to the ultimate load at room temperature. In this study, 0.5 is taken as the typical working condition, and other working conditions include 0.3 and 0.7 to examine the influence of load levels. The fire action adopts the ISO 834 standard temperature-time curve, which is expressed as Equation 12: T fire = 20 + 345 ⁡ log 10 8 t + 1 (12) where T fire is the furnace temperature (°C) and t is the time (min). The bottom and two side surfaces of the beam were exposed to the fire, while the top surface was not exposed to the fire, simulating the typical fire conditions of floor beams in a fire.

A heat transfer analysis was performed independently before the structural analysis. The calculation of the thermal parameters for concrete and steel, with values for thermal conductivity, specific heat capacity, and density, followed the recommendations in Eurocode 2, using temperature dependence where applicable. The thermal conductivity is a function of temperature, decreasing with an increase in temperature from 1.6 W/(m·K) for room temperature to 0.8 W/(m·K) for 800 °C. The specific heat capacity reaches its maximum between 100 °C and 200 °C because of the water evaporation effect; it then increases with an increase in temperature. Both convection and radiation heat transfer coefficients for the surface in contact with the fire environment were 25 W/(m²·K) and 0.7 respectively. The heat transfer coefficient between the unexposed top surface and the environment was 9 W/(m²·K) to simulate natural convection. The initial temperature was set to 20 °C, and the heat transfer analysis was performed in transient mode with an adaptive time step to ensure convergence and accuracy. After the heat transfer analysis, the nodal temperature histories were exported and imported as predefined fields into the structural analysis. Structural analysis also uses implicit dynamics steps to solve, taking into account geometric nonlinearity and material nonlinearity, and the time increment step is automatically adjusted according to load and temperature changes.

The model was validated using three sets of fire test data from reinforced concrete beams reported in the literature. A total of five validation cases were selected from these sources. Table 2 summarizes the experimental configurations and material properties, with selection criteria including: (1) ISO 834 standard fire exposure; (2) comparable geometric scale; (3) documented deflection-time history and failure time; (4) normal-strength concrete with siliceous aggregate.

The three sources include Gernay et al. (2013), Dwaikat and Kodur (2009), and Kodur and Banerji (2024), covering variations in load ratio, cover thickness, and reinforcement ratio. These test results offered critical information on temperature distribution, time-dependent curve of the mid-span deflection, and time to failure for the beams in the test, making it a valid basis for checking the correctness of the finite element simulation model. The concretes’ strengths, reinforcement, and temperature rise in fires during the tests were more or less similar to the simulation conditions in this simulation, thereby making it valid for verification purposes. The parameters for the finite element simulation model are listed in Table 3 below.

Six design parameters were identified for the sensitivity analyses in evaluating the effect of design parameters on the fire resistance of reinforced concrete beams. The identified parameters were load ratio, thickness of concrete cover, ratio of longitudinal reinforcement, section heights, number of fire-exposed surfaces, and the angle of parameter expansion for models. The ratio λ is changed between 0.3, 0.5, and 0.7, expressing low, medium, and high levels of loading, respectively. The thickness c is changed between 25, 35, and 45 mm, reflecting the permissible values according to the code standards for different structures. The ratio ρ, also referred to as the ratio of the area of tensile reinforcement to the effective area, is changed between 1.0%, 1.5%, and 2.0% for under-reinforced, normally reinforced, and over-reinforced structures, respectively. The section heights were changed between 350 and 450 mm for an investigation into differences in cross-sectional dimensions. Fire exposure from either three sides for floors, independent beams, columns, or four sides for uninterrupted beams is considered in these analyses. The dilation angle ψ is a CDP model parameter defining the ratio of volumetric to shear plastic strain in the p-q stress plane (Equation 6). It governs concrete dilatancy under multiaxial stress states, with typical values ranging from 20 ° to 40 ° for normal-strength concrete (Vermeer and de Borst, 1984). The selected range (20 °, 30 °, 40 °) encompasses commonly adopted values in structural fire analysis to assess model sensitivity to this parameter.

To further quantify the relative importance of each parameter on the response, the Pearson correlation coefficient was used for statistical analysis. The correlation coefficient was calculated based on the parameter-response data pairs of all simulated conditions. The absolute value of the Pearson coefficient r is defined as a strong correlation above 0.6, a moderate correlation between 0.3 and 0.6, and a weak correlation below 0.3. A positive coefficient indicates that an increase in the parameter leads to an increase in the response, while a negative coefficient indicates the opposite relationship. The main evaluation indicator is the fire resistance time, which is defined as the duration from the onset of the fire to the failure of the beam. In addition, the characteristics of the failure mode are also recorded, such as whether it is ductile failure dominated by yielding of the steel bars or brittle failure dominated by crushing of the concrete. At high temperatures, the failure mode may change from plastic hinge failure at room temperature to sudden failure after softening of the steel bars at high temperature, or rapid collapse of the concrete in the compression zone under the combined action of high temperature and high stress.

The operating conditions for the sensitivity analysis were summarized in a table format. The baseline condition served as the reference group, and each parameter variation condition was compared with the baseline condition. Besides univariate analyses, an orthogonal experimental design was also adopted to explore the interaction effects among parameters. Unlike the other five parameters which are physical design variables, the dilation angle is a numerical constitutive parameter that primarily affects local stress redistribution rather than global fire resistance. Therefore, the dilation angle was excluded from the orthogonal design to focus on practically meaningful design variables. The orthogonal experiment used an L16 (4^5) orthogonal matrix, choosing factors such as load ratio, cover, reinforcement ratio, section height, and number of visible surfaces, with each factor having four levels. In orthogonal experiments, it is possible to effectively traverse the entire factor space using only a limited number of experiments, making it possible to study the variance significance levels and main effects for different factors. The finite element simulation cases comprised: (1) single-factor analysis with 13 cases (1 baseline case plus 2 additional levels for each of 6 parameters); (2) L16 (4^5) orthogonal array with 16 cases; and (3) 10 supplementary Latin hypercube sampling points for surrogate model training. After removing four duplicate configurations, a total of 35 unique finite element simulations were conducted (13 + 16 + 10 − 4 = 35), which is acceptable from a computational viewpoint for understanding quantitative relationships between parameters and responses.

Based on the validated finite element model and parameter sensitivity analysis results, this study established a framework for the optimized design of reinforced concrete beam sections. The goal of the optimization problem is to minimize the beam’s material usage while meeting fire resistance requirements, thereby reducing construction costs and environmental impact. The design variable vector x = b , h , A s , c T consists of four independent variables: beam width b , beam height h , tensile reinforcement area A s , and concrete cover thickness c . These variables range from b ∈ 200 , 400 mm, h ∈ 300 , 500 mm, A s ∈ 600 , 2500 mm², and c ∈ 25 , 50 mm, respectively. The objective function of the single-objective optimization problem is defined as the beam’s cross-sectional area f x = b × h , aiming to minimize material usage. Constraints ensure the feasibility and safety of the design, including: fire resistance constraint g 1 x : t f i r e x ≥ 120 min, requiring the fire resistance time to be no less than 2 h; deformation limit constraint g 2 x : u max x ≤ L / 30 , limiting the maximum deflection during a fire; reinforcement ratio constraint g 3 x : 0.8 % ≤ ρ x ≤ 2.5 % , ensuring that the reinforcement ratio is within a reasonable range; protective layer thickness constraint g 4 x : c ≥ 25 mm, meeting the basic requirements of the specification; and geometric dimension constraint g 5 x : h / b ≥ 1.5 , avoiding overly flat sections. The fire resistance time constraint is an implicit constraint and needs to be evaluated through finite element simulation or proxy model. Table 4 summarizes the complete multi-objective optimization problem formulation, including design variables, objectives, and constraints.

To address multi-objective optimization requirements, this study simultaneously considers the conflicting objectives of material usage and fire resistance. The multi-objective optimization problem is formulated as: minimize  F x = f 1 x , f 2 x T $ , w h e r e $ f 1 x = b × h represents the material usage and f 2 x = − t f i r e x represents the negative fire resistance time. By converting the negative fire resistance time into a minimization problem, both objective functions can be solved within a unified framework. Each point on the Pareto optimal frontier represents a design solution that achieves the optimal trade-off between the two objectives. Designers can select appropriate design points from the frontier based on their specific project preferences for cost and safety.

The optimization algorithm uses a genetic algorithm for single-objective problems and the NSGA-II algorithm for multi-objective problems. The genetic algorithm parameters were set as follows: population size 50, 100 iterations, crossover probability 0.8, and mutation probability 0.1. Tournament selection was used for selection, arithmetic crossover was used for crossover, and Gaussian perturbations were applied for mutation with adaptive decay over the number of iterations. The NSGA-II algorithm maintained the Pareto-optimal solution set using fast non-dominated sorting and crowding distance. The population size was 100 and the iterations were 200. The algorithm divided the population into different non-dominated hierarchies, prioritizing individuals at the highest level. Within each hierarchy, individuals were selected based on crowding distance to maintain a uniform distribution at the frontier. An elite retention strategy ensured that the optimal solution was not lost in each generation.

To reduce computational costs, the optimization process incorporated Kriging surrogate model technology. While direct finite element simulations would take 4–6 h to evaluate each candidate design, the surrogate model uses Gaussian process regression to establish a probabilistic mapping between design variables and responses, allowing a single evaluation to take only 0.01 s. The Kriging model was trained on 35 finite element results using Latin hypercube sampling, approximately following the 10 × d guideline for four design variables. Leave-one-out cross-validation (R ² = 0.92, RMSE = 4.7%, maximum absolute error = 11.2%) confirmed adequate surrogate accuracy; learning curve analysis showed R ² stabilizes above 0.90 when n > 25. The surrogate reduces per-evaluation time from 4.2 h to 0.01 s, and final Pareto solutions were validated using high-fidelity finite element analysis. This surrogate modeling strategy aligns with recent AI-based prediction frameworks for structural fire analysis (Ali et al., 2024; 2025a; 2025b).

Model validation was conducted at both the material and structural levels. Material-level validation first examined the model’s ability to predict the uniaxial compressive stress-strain curves of concrete at different temperatures. Figure 3 compares the experimental and simulated curves at four temperatures: 20 °C, 200 °C, 400 °C, and 600 °C. At room temperature (20 °C), the model accurately captures the initial elastic phase, nonlinear rise phase, peak strength, and softening phase of concrete. The peak stress is 30.2 MPa, highly consistent with the nominal strength of 30 MPa, and the corresponding peak strain is approximately 0.002. The downward slope of the softening phase reflects the brittle nature of concrete, and the model predictions are consistent with the experimental data. At 200 °C, the peak strength decreases slightly to approximately 27.0 MPa, approximately 90% of the room temperature strength, indicating that concrete performance degrades slowly at low temperatures. At 400 °C, the peak strength decreases to 22.5 MPa, approximately 75% of the room temperature strength, while the peak strain increases to 0.0028, indicating improved ductility at elevated temperatures. At 600 °C, the strength further decreased to 13.5 MPa, only 45% of the room temperature strength, with a peak strain of 0.0038. The slope of the softening section increased, indicating that high-temperature concrete rapidly loses its bearing capacity after the peak. The model’s response predictions at all temperatures agreed well with the experimental data, with an average prediction error of 7.9%, validating the rationality of the temperature-dependent damage evolution law. At the material level, the proposed model achieves 7.9% average error versus 12%–15% for conventional CDP with temperature-independent damage parameters, particularly improving prediction accuracy in the 400 °C–600 °C range where damage-temperature coupling is most significant.

Structural verification focused on the deflection response and temperature distribution of reinforced concrete beams subjected to fire. Figure 4a compares the mid-span deflection development under three different conditions with experimental data, clearly demonstrating the complete response of the beam under fire. The first condition is the baseline condition, with a load ratio of 0.5, a 35 mm cover, and fire on three sides. The deflection curve is generally characterized by an S-type curve feature. At the early stage (0–20 min), the deflection is small and even positive, due to the upward deflection resulted from the heat expansion effect countering the downward deflection effect caused by the weight load. At the middle stage (20–100 min), the deflection increases faster, due to the reduced mechanical properties under high temperature, thus weakening the strengthened beam’s weight-bearing capability. At the late stage (100–125 min), the deflection increases quickly, implying that the strengthened beam is close to its critical state in terms of weight-bearing capability. Finally, at 124 min, the deflection exceeded the L/30 criteria, thus leading to a failure. The simulation took 126 min for the beam to fail with only a 1.6% error in comparison with the experimental result, which occurred at 124 min, making it highly accurate. In the second test condition, the thickness of the protective layer is extended to 45 mm, with other test parameters remaining unchanged. There is a right-shift in the deflection curve, with the time to failure lengthening to 152 min, with an error of 2.7% in comparison with the experimental result, which occurred at 148 min. The protective layer’s thickness increase resulted in a marked slowing down in the heating rate for the steel bar, thus increasing the fire resistance time accordingly. In the third test condition, the weight ratio is raised to 0.7, with the thickness of the protective layer reduced back to 35 mm. The deflection curve becomes steeper, with the time to failure shortening to 82 min, with an error of 5.1% in comparison with the experimental result, which occurred at 78 min. The strengthened beam’s stress reserve is considerably small for high weight ratios, thus any slight weakening in the mechanical properties would easily cause a system failure.

Figure 4b illustrates the verification results for the temperature distribution along the height of the cross-section at 60 min. The finite element calculation modeled a beam with the following cross-sectional area measurements: 200 by 400 mm. The bottom, two sides, but no top, were subjected to fire exposure. The temperature gradients obtained match well with the experimental values, thus validating the heat transfer simulation results. The surface temperature exposed to fire reached approximately 840 °C, close to the furnace temperature corresponding to the ISO 834 curve. The temperature gradient along the cross-section height was approximately 600 °C at 50 mm from the bottom, approximately 450 °C at 100 mm, approximately 300 °C at the center, and only approximately 200 °C at the top. This non-uniform temperature gradient leads to complex thermal strain and stress distributions within the cross-section, which is the root cause of the complexity of the structural behavior under fire. The temperature at the rebar location depends on the cover thickness. At 35 mm from the bottom, the center of the bottom rebar reached approximately 650 °C, entering the region of significant rebar strength degradation.

The evolution of the temperature field in the beam section is shown in Figure 5. Figures 5a,b show the temperature contours of the beam section after 60 and 120 min of fire exposure, respectively. At 60 min, the temperature of the exposed surface reached approximately 900 °C. The temperature gradient along the cross-section height was observed, with temperatures reaching approximately 600 °C at 50 mm from the bottom, approximately 450 °C at 100 mm, approximately 300 °C at the center, and only approximately 200 °C at the top. The bottom reinforcement temperature reached approximately 650 °C, entering a region where reinforcement strength had significantly degraded. At 120 min, the temperature increased further and penetrated deeper, reaching approximately 1,000 °C on the exposed surface, 680 °C at the bottom reinforcement temperature, approximately 550 °C at 100 mm from the bottom, and approximately 400 °C at the center and 300 °C at the top. The simulated temperature distribution agrees well with the experimental measurement points, validating the accuracy of the heat transfer model. This non-uniform temperature gradient distribution leads to complex thermal strain and stress fields within the cross section, which is the fundamental reason for the complexity of structural response under fire.

The quantitative results of the fire resistance time validation are summarized in Table 5. In addition to the three groups mentioned above, two additional groups of beams with different reinforcement ratios are included. As shown in Table 5, the fire resistance time of the validation beam with a reinforcement ratio of 1.0% (load ratio 0.6, cover 30 mm) was 108 min, compared to the simulated value of 110 min, with an error of 1.9%. The fire resistance time of the validation beam with a reinforcement ratio of 2.0% (load ratio 0.55, cover 30 mm) was 98 min, compared to the simulated value of 95 min, with an error of −3.1%. The average prediction error for the five conditions was 2.9%, with a standard deviation of 3.1%, demonstrating the model’s good predictive stability under different parameter combinations. It should be noted that these conditions were derived from experimental studies in different literature, with varying design parameters (such as load ratio and cover thickness). The purpose of this study was to validate the model’s applicability under a variety of conditions. The single-factor effect of reinforcement ratio on fire resistance will be systematically discussed in the parameter sensitivity analysis in Section 3.2. Furthermore, the simulation reproduced the failure mode observed in the tests: plastic elongation after yielding of the tensile reinforcement, concrete cracking at the bottom of the mid-span, and ultimately, concrete crushing or reinforcement failure in the compression zone, leading to loss of the beam’s load-bearing capacity. Combined validation results at both the material and structural levels confirmed the accuracy and reliability of the temperature-dependent CDP model and its finite element implementation, laying the foundation for subsequent parameter analysis and optimized design.

To demonstrate the advantages of the proposed model, additional simulations were conducted using the conventional ABAQUS built-in CDP model. As shown in Table 6, the proposed model achieves an average prediction error of 2.9%, compared to 11.4% for conventional CDP. This improvement is attributed to the thermodynamically consistent damage-temperature coupling, which captures accelerated stiffness degradation in the 400 °C–600 °C range more accurately.

Parameter sensitivity analysis helped in understanding the effect of different design parameters on the fire-resistance ratio of reinforced concrete beams. Six significant parameters, namely, load ratio, cover thickness, reinforcement ratio, section height, number of fire-exposed sides, and expansion angle, were investigated using 30 finite element analyses with an orthogonal experimental design. The single-factor variation method was used in the investigation, where each factor is altered while others were held at their base values (load ratio 0.5, cover thickness 35 mm, reinforcement ratio 1.5%, section height 400 mm, number of fire-exposed sides 3, and expansion angle 30 °) to effectively establish the independent contribution made by each parameter to the fire-resistance ratio value. Fire-resistance time is adopted here as the principal responses in the investigation, whereby it is measured from the time fire is initiated to the time the structure collapses. Other factors, namely, deflection in the middle span in excess of 1/30 times the span length, with a loss in load-carrying capacity to less than 50% of original values, were also considered for structural collapses in the study. The base case fire-resistance time of 126 min matches validation Case 1 (Section 3.1), confirming model reliability.

The value of the load ratio is the key factor influencing the fire-resistance time for beams. As the load ratio increases from 0.3 to 0.7, the fire-resistance time for the beam reduces sharply from 178 to 75 min, resulting in a 58% decrease (Figure 6a). This strong negative correlation reflects diminished stress reserves at higher load ratios, where minor strength degradation triggers structural failure. At a load ratio of 0.3, the beam is in a state of reduced stress during the early stages of the fire, whereby even at 600 °C–700 °C, the nominal structural strengths of the steel as well as the concrete would remain intact to resist the loads. At a load ratio of 0.7, the beam is in a state where it reaches almost 70% of its loading capacity at ordinary temperature levels. When the temperature reaches around 400 °C, the nominal strengths of the steel reinforcement would reduce by almost 20%, leading to an absence of margins of safety. Additionally, it is also discerned from the analysis that the critical value for the load ratio is fixed at 0.55. At values higher than 0.55, the nominal increase in the fire-resistance time with respect to the load ratio increases with evidence from the drastically higher slope in Figure 6a.

Cover thickness has the significant impact on fire-resistance time in an inverted manner. As the thickness of the cover increases from 25 mm to 45 mm, the fire-resistance time rises from 85 to 165 min, thereby rising by 94% (Figure 6b). The underlying physical phenomenon for these positive results is the heat insulation offered by the cover thickness. An increase in the thickness of the concrete cover results in an increase in the distance of heat travel from the surface to the rebar core, with the additionally low heat conduction offered by the concrete further retarding the heat travel phenomenon. After 60 min into the fire, the temperature in the rebar core under the 25 mm thick cover has risen to 720 °C, in comparison with only 480 °C in the rebar core under 45 mm thick covers. This difference in temperature, with readings reaching 240 °C, resulted in the rebar in the 25 mm thick cover having only 40% of its strength in room temperature values, in comparison with values over 70% in rebar with 45 mm thick covers. There is an optimum thickness for covers, though. When the cover thickness exceeds 50 mm, the steam pressure generated by the high-temperature concrete can cause the surface layer to crack and fall off.

The reinforcement ratio shows a weak negative correlation with the fire resistance limit time. When the reinforcement ratio increases from 1.0% to 2.0%, the fire resistance limit time decreases by approximately 16%, from 138 to 116 min (Figure 6c). This nonlinear relationship reflects a competition between two mechanisms: while a higher reinforcement ratio provides a greater load reserve at room temperature, more steel bars mean a larger cross-sectional area. The high thermal conductivity of steel bars accelerates temperature diffusion within the cross-section, causing the strength of steel bars to decay faster than that of concrete. Therefore, a high reinforcement ratio makes the structure more dependent on steel bars, and once the steel bars soften, the beam will fail rapidly. The interplay of positive and negative effects results in a relatively small net effect of the reinforcement ratio on the fire resistance time. Analysis indicates that the optimal reinforcement ratio range is approximately 1.0%–1.4%. Increasing the cross-section height has a significant positive effect on the fire resistance time. When the cross-section height increases from 350 to 450 mm, the fire resistance time increases by 49%, from 95 to 142 min (Figure 6d). This increase is attributed to the fact that the taller cross-section provides a larger moment arm and prolongs the time it takes for the concrete in the compression zone to reach high temperatures. The effect of the number of exposed surfaces was as expected. Fire resistance lasted approximately 25%–30% longer in a three-sided fire than in a four-sided fire (Figure 6e). In a three-sided fire, the top, unexposed surface acts as a heat dissipation channel and provides strength reserve for the compression zone. The expansion angle, a model parameter, had no significant effect on fire resistance within a range of 20 °–40 ° (Figure 6f), with the difference within 5%, confirming the robustness of the model within a reasonable parameter range.

Quantitative correlation analysis assessed the strength of association between each parameter and fire resistance time. Pearson correlation coefficients showed that the load ratio (r = −0.82) and cover thickness (r = +0.78) were strongly correlated with fire resistance time (p < 0.01), while section height (r = +0.51) and reinforcement ratio (r = −0.35) showed moderate correlations (p < 0.05), and the influence of expansion angle was negligible (r = +0.12, p > 0.1). An orthogonal analysis of variance further confirmed the order of parameter importance. Based on the F statistic, cover thickness (F = 42.3) and load ratio (F = 38.7) were the most significant factors, followed by section height (F = 18.5) and reinforcement ratio (F = 9.2). The number of exposed surfaces (F = 7.8) and expansion angle (F = 1.5) had a minor influence. The cover and load ratio together explained 71% of the total variance and were the key variables controlling the fire resistance performance, while other parameters accounted for 29%.

This section uses a genetic algorithm and the NSGA-II algorithm for single-objective and multi-objective optimization, respectively, to minimize material usage while satisfying the fire resistance time constraint. Conventional design methods typically employ conservative cross-sectional dimensions to ensure compliance with the 2-h fire resistance rating requirement. As shown in Figure 7, the conventional design (Conv.) has a cross-sectional area of 150,000 mm² (300 mm wide × 500 mm high), is equipped with four 20 mm diameter steel bars, and has a 30 mm cover. Finite element analysis verifies that its fire resistance time is 125 min. This optimization algorithm significantly reduces material usage while maintaining the same fire resistance constraint and provides designers with a variety of performance-cost trade-offs.

The multi-objective optimization problem for using materials and fire resistance time is solved using the NSGA-II algorithm. After an iterative search procedure, 23 solutions in the Pareto front were obtained (Figure 7). The front has a nonlinear property with a convex region towards the origin, indicating the trade-off between the two objective values. There are three regions in the front: the high-efficiency region (cross-sectional area 100,000–115,000 mm²), which is a steep area, indicating a small increase in material usage offset by high gains in fire resistance time; the optimal trade-off region (115,000–140,000 mm²), indicating that an acceptable level is attained between usage and efficiency; and the low-efficiency region (140,000–175,000 mm²), which is a flat region indicating large amounts of usage yielding small increases in efficiency. The bottom left corner corresponds to the minimal usage solution (Design A with 102,900 mm² cross-sectional area and 110 min fire resistance time), which only satisfies minimal requirements for usage. The top right corner corresponds to maximal fire resistance time (cross-sectional area 175,000 mm² with 175 min fire resistance time). The front is continuous with optimal solutions providing designers with an apt choice according to project importance factors for safety in regard to assigned economical bounds.

Three samples from the Pareto front were picked for in-depth investigation (Table 7). The three representative designs were selected based on their positions on the Pareto front: Design A at the minimum-material extreme representing cost-priority applications; Design B at the knee point (maximum curvature) representing optimal cost-performance trade-off; Design C at the high-safety region for critical infrastructure requiring enhanced fire resistance margins. Design A (Minimum Material) emphasizes maximum saving on material. Its cross-sectional dimensions are 245 by 420 mm, reinforcement is 3φ20 mm, while its cover is 38 mm. The area of its cross-section is 102,900 mm², with a material saving percentage of 31.4% more on material saving efficiency compared to traditional designs. Design B (Recommended), which is at the middle position on the curve for maximum efficiency, is the optimized mix between saving on materials versus optimized performance. Its dimensions for its cross-sectional area consist of 260 by 440 mm, with 3φ22 mm rebar, while its cover is 40 mm thick. The area for its cross-sectional dimension is 114,400 mm², with 23.7% more saving on its material compared to conventional designs. Design C (Maximum Safety), at the other end on the curve for maximum safety, gives maximum emphasize on safety considerations over other factors. Its dimensions for its cross-sectional area consist of 285 by 480 mm, with 4φ20 mm rebar, while its cover is 45 mm thick. The area for its cross-sectional dimension is 136,800 mm², with an 8.8% saving on its material compared to conventional structural designs. All three designs gradually increase the cover thickness from 38 to 45 mm, a key parameter for delaying the temperature rise of the steel bars and improving fire resistance. In engineering practice, designers can use the Pareto front as a decision-making tool to select the design from the front that best meets the project requirements based on the owner’s preference for cost and safety. Recommended Design B is located at the critical inflection point of the front. Moving it to the left saves material but sacrifices an excessive safety margin, while moving it to the right increases cost but does not significantly improve performance. Therefore, it is the most cost-effective option.

A complete finite element verification was conducted on three representative design schemes, and the results are shown in Figure 8. The temperature distribution cloud diagram in Figure 8a shows that at 120 min, the temperature of the bottom fire-exposed surface exceeds 900 °C, while the top remains between 200 °C–300 °C, exhibiting a significant temperature gradient. The temperature of the steel bars is controlled by the thickness of the protective layer. The deflection-time curves in Figure 8b all exhibit the typical three-stage characteristics of a slight upward arch in the early stage, accelerated downward deflection in the middle stage, and a sharp increase in the late stage. As shown in Figure 8b, the Kriging surrogate model predicted that Designs A, B, and C would reach the L/30 deflection limit at 108.5, 127.8, and 152.3 min, respectively. High-fidelity FE verification confirmed fire resistance times of 110, 128, and 150 min for the three designs, yielding prediction errors of −1.4%, −0.2%, and +1.5%, respectively, with an average absolute error of only 1.0%. This close agreement validates the accuracy of the surrogate model and the reliability of the optimization procedure. From Figure 8c, it can be seen that the temperature history for the steel bar in Design C with 45 mm thick cover is 600 °C at 120 min, while for Designs B with 40 mm thick cover, it is 620 °C; for Design A with 38 mm thick cover, it is 650 °C. Additionally, it is also clear that for every 5 mm increase in thickness, the temperature falls by 30 °C.

The rolling over limit for each of the three designs corresponds to reaching the deflection limit, in accordance with the fire resistance law for reinforced concrete beams. Design A achieves 31.4% material savings with 108.5 min fire resistance, suitable for areas with reduced fire protection requirements. Design B (23.7% savings, 127.8 min) is recommended for general residential structures. Design C (8.8% savings, 152.3 min) suits critical structures and densely populated areas. The design optimization method presented in the research yields simultaneous benefits over traditional approaches in the following areas: Material saving improvement up to 31.4%, accuracy in prediction improved to 3%, and instead of a single conservative design, the set of design solutions is enlarged to encompass a Pareto-optimal solution set. With the aid of the Pareto front, designers can simultaneously pursue safe and economically optimal designs for structures. The proposed theoretical structure can be extended for other types of components, henceforth providing technical aid for performance-based fire protection design approaches in practice.