Numerical modeling and design of precast prestressed UHPFRC I beams

: Ultra-high-performance fiber reinforced concrete (UHPFRC) is a new material developed to present superior properties, as high compressive strength (higher than 130 MPa), high durability, and satisfactory ductile behavior. This paper reports the procedure to design precast UHPFRC beams, subject to flexural loads. First, an I-girder AASHTO Type II was designed, and simulated. Next, I-beams with diverse depths and steel ratios were designed, and simulated considering a four-point bending load test. It was found that the classical design equations used to predict the strength bending moment (M rd ) showed good accuracy with the simulated models with a 4.5% error.

The notable advances of superplasticizers and mineral additions allowed the production of concretes with high volumetric amounts of fine particles, perfectly packed and with low water-cement ratios allowing the development of the precursor Ultra-High-Performance Concrete (UHPC) [6]- [8]. UHPC presents (i) high mechanical properties in the long term, (ii) low permeability, and (iii) long life cycle under aggressive environments [6], [9]. Therefore, UHPC must present a water-cement ratio of around 0.2 and a minimum compressive strength of 130 MPa at 28 days [10], [11]. However, the matrix densification and decrease of microstructure imperfections may induce brittle failure of the elements that present low ductility due to its reduced tensile strength. The mitigation of brittle failures, metallic microfibers are added to the mixture, improving the ductile behavior, and then UHPC is named UHPFRC [7].
Considering its elevated properties, durability, and appropriated workability, UHPFRC is an emergent material with notable applications in the precast industry [7], [12]- [14]. Given this context, this paper describes a procedure for designing UHPFRC elements subjected to bending loads. The analytical results obtained with the proposed equations are compared to numerical simulation results of beams constituted of UHPFRC, presenting good accuracy.

BENDING DESIGN
Fehling et al. [7], [15] proposed basic assumptions and analytical equations to design UHPC and UHPFRC elements subjected to flexural loads. Figure 1 presents the stress-strain distribution of the rectangular beam studied by Fehling et al. [7], [15] considering axial force (N sd) and bending moment (Msd): In Figure 1, x is the neutral axis position; Fcc is the concrete compressive resultant force located at a distance equal to x/3 from the top; Fft is the concrete tensile force located at the centroid of the parabolic area representing the stress distribution; and Fst is the resultant steel force.
The herein presented assumptions are considered to develop the design equations: • Bernoulli hypothesis: plane sections remain plane, and the deformed beam angles are small. • Triangular compressive stress distribution is considered with a linear stress-strain compressive response for the concrete until the failure (see e.g., Figure 1). • Tensile stress cannot be disregarded for UHPFRC, and its distribution can be considered parabolic (Figure 1a) or rectangular ( Figure 1b) The concrete tensile resultant force can be calculated by integrating the area of the stress distribution. Equation 1 indicates the resultant force considering a parabolic stress distribution. The resultant force is located at the centroid of the parabolic area from a distance equal to 0.56 (h -x) to the neutral axis [7], [15]. Equation 2 indicates the concrete tensile force considering a rectangular distribution. Comparing Equations 1 and 2, it can be observed an error smaller than 2.5%. Hence, for simplicity, this paper considers the rectangular distribution for concrete tensile stresses, according to Equation 2.
In Equation 2, h is the beam height, x is the depth of the neutral axis, b is the beam width, σcf,0d is the direct tensile strength.
Equation 3 gives the resultant compressive force (Fcc), and Equation 4 presents the resultant force applied to the prestressed cables (Fst) calculates using the stress-strain diagram of the steel cables.
The equilibrium equations in terms of the axial forces and bending moments are given by Equations 5 and 6, respectively, Finally, it is necessary to impose a strain compatibility relationship for the cross-section (Equation 7). Figure 1c shows the strain response along with the height of the section, considering the cables embedded in the concrete.
In Equation 7, (εst) is the tensile strain at the reinforcement and (εc) the compressive strain of the concrete. Finally, following the above-suggested steps, it is possible to design a beam with arbitrary dimensions (Figure 2), 1. Calculate the resultant axial forces: Fft, Fcc, Fsc e Fst; 2. Define the depth of the neutral axis (x) using an iterative algorithm to promote the equilibrium of bending moments in the cross-section; 3. Update resultant forces at concrete and steel cables Fcc and Fft; 4. Impose the compatibility relationship to calculate the strains along with the cross-section height; 5. Determine the stress in the reinforcements σs,t using the constitutive stress-strain relation; 6. Calculate the required reinforcement area (Ast); Where: k = 0.9 χ = 0.85 if the width of the cross-section decreases towards the tensile edge χ = 0.90 in general

CASES OF STUDY
In this section, the numerical and experimental results showed by Graybeal [16] are used to validate the proposed method for designing UHPFRC elements (section 2). The example consists of an I-beam PCI AASHTO II constituted by UHPFRC with 26 prestressed steel reinforcements and subjected to a bending test, according to Figure 3a and 3b. This section is commonly applied to bridge structures using C70 concrete. The Graybeal's experimental results were used to calibrate a Finite Elements Model (FEM) using Concrete Damage Plasticity (CDP) to describe the non-linear behavior of the material, presented in Section 5.1.
Moreover, aiming to test the capability of the analytical equations proposed by Fehling et al. [7], [15], the strength of prestressed sections are predicted using a numerical FEM model and analytical model. These beams present heights h = 500 mm, 400 mm and 300 mm, two reinforcement strands (ϕ = 12,5 mm) placed in the inferior flange, constant width b = 300 mm and constant thickness e = 50 mm, according to Figure 3c. The spans of the beams are 3 meters, and the point load application and boundary conditions are described in Figure 3d. The UHPFRC characteristics used in this analytical-numerical model were adopted, according to Krahl et al. [17], [18].  Table 1 presents the different material characteristics for UHPFRC produced by Krahl and Graybeal. It is relevant the fact of the Graybeal's concrete presents divergency in uniaxial tensile behavior, i.e.,: the value of 9MPa was presented by the author as the direct tensile strength, performed in dog bone samples [16]; and the value of 15.9MPa was presented as the calibrated uniaxial response of numerical model [16]. In this paper, we are showing the results obtained, considering these two different tensile strengths to analyze. Hence, the first girder "Complete Graybeal model" was simulated using the tensile strength of 15.9 MPa and considering the tensile and compressive uniaxial behaviors; and "Simplified Graybeal model" was modeled using the tensile strength of 9 MPa and a simplified uniaxial tensile law based in the direct tensile tests. *d' is the distance of the top and bottom face until steel strands gravity center. ** The value of 9 MPa is obtained by direct tensile test [16]. ***The value of 15.9 Mpa is obtained by numerical-experimental calibration [16] 4 NUMERICAL SIMULATION

Graybeal's girder constitutive model
The constitutive tensile law based on the calibrated tensile and compressive uniaxial behaviors given by Graybeal [16] is according to Figure 4a and b, with 15.9 MPa tensile strength. The second girder, "Simplified Graybeal model", was modeled using the tensile strength of 9 MPa and a simplified uniaxial tensile law based on direct tensile tests; see Figure 4b simplified model. The constitutive model of Concrete Damage Plasticity (implemented in ABAQUS CAE Simulia) was applied to simulate both girders. The compressive and tensile damage evolution (i.e., dt and dc) are obtained by the methodology proposed by Birtel and Mark [19], through Equations 8 and 9: In Equations 8 and 9 is the tensile stress of the concrete; c σ is the compressive stress of the concrete; Ec is the elastic modulus of the intact material; bc = 0.7 and bt = 0.57 are empirical parameters related to damage evolution [19]; and are plastic strains, defined as , and pl in in c b ε ε ε = is given by the difference between total and elastic deformation (σ/Ec). Table 2 presents the elasticity and plasticity parameters adopted.
ASTM 270-ksi steel was used in prestressed strands, with the constitutive law given by Figure 5 [20]. A prestress load of 885 MPa was adopted. The prestressed strands were considered totally embedded into the concrete. The steel elastic modulus is Es = 197 GPa.

Constitutive model (I-beams)
For the three simulated I-beams (Figure 3c and 3d), the UHPFRC constitutive law written in terms of the inelastic strain is applied according to Figures 6a and 6b, and Krahl et al. [17], [18]. Plasticity and elasticity parameters are given in Table 2; the design parameters and prestress load are presented in Table 1.

Boundary conditions, mesh, and load
All beams were considered simply supported, and a 3D 8-nodes solid elements mesh was applied (C3D8R) to represent UHPFRC, with 3 degrees of freedom each node. Graybeal's girders were modeled with a 50 mm mesh size for transversal direction and 150 mm for longitudinal axis. I-beams were modeled using a 40 mm mesh size for longitudinal direction, and 20 mm mesh size for transversal discretization. 3D truss finite elements were used to represent the prestressed strands with a 100 mm mesh size. The gravity load was considered as -9.81 m/s 2 , and concrete density equals to 2500 kg/m 3 .  [17], [18] A displacement control test was used for the three studied I-beams and for the Simplified Graybeal model. For the Complete Graybeal model, the same experimental methodology given in Graybeal [16] was performed, i.e., initial force control until 18 kN, and posterior displacements control until failure. Therefore, the following cases were studied in this paper: (I) Complete Graybeal model: force and displacements control, with constitutive law calibrated according to experimental results of Graybeal [16]; (II) Simplified Graybeal model: displacements control and simplified constitutive law based in the direct tensile tests; (III) I-beams: displacement control and constitutive law was given by Krahl et al. [17], [18]. Figure 7a presents the comparison between the numerical results and experimental behavior given by Graybeal [16]. The force-displacement curve obtained by the complete Graybeal model presented a 5% maximum error (gray area in the graph). The estimated maximum force was 3.17% higher in the numerical model. The displacements prediction was 5.6% higher in the numerical model. For the Simplified Graybeal model, a maximum error of 10% was achieved in the force x displacement curve comparing the numerical and experimental results.
Relevant parameters are the tensile and compressive damage indexes (Figure 7b and 7c) over the load evolution. Yang et al. [21] showed that damage indexes are used to predict the cracking zones, which can be useful for the identification of the failure modes. For the Graybeal's girders bending failures were observed, presenting predominant tensile damage indexes in the middle of the spans. When the peak force is reached, both models presented tensile damage indexes of 90% (Figure 7c), with excessive strains in the reinforcements, characterizing the flexural failure.

Analytical design of Graybeal's girder
The analytical equations showed in section 2 are used to estimate the maximum bending moment of Graybeal's girders. Table 1 presents the material characteristics applied to prestress force and section parameters. Figure 8 presents a geometrical law in terms of the depth of the neutral axis and the compressed area, used in Equations 13 to 16. For the analytical design, better results were obtained with the tensile strength of 9 MPa.
The maximum strain in the compressive behavior is given by Equation 10, following the French Association of Civil Engineering (AFGC) [22]: In Equation 10, fctm is the tensile strength, which could be higher than cracking stress when UHPFRC presents strain-hardening behavior; fcm is the average compressive strength; εc0d is the strain written in terms of the elastic limit (fcm/Eci), and Ec is the young modulus initial tangent.
The calculus of εcud allows the determination of the limit between strain domains 2 and 3 of the Brazilian code ABNT NBR 6118:2014 (i.e., concrete with a strain of 5.8 ‰ and steel strands with 10 ‰). The determination of the neutral axis between dominium 2 and 3 (i.e., xlim) is given by Equation 11: The transversal section was divided into strips of 0.01 m heigh to account for the cumulative areas. Through the geometrical law, the geometrical law, it was possible to determine the values of Ac = 0.076 m 2 (compressed area), and At = 0.154 m 2 (tensile area) for the strain domain limit 2 and 3 (xlim); the force components for this configuration are expressed by Equation 12 to 15:  In Equations 12-15, Acc is the concrete compressed area, nbar is the total number of strands, Atot is the total area of the cross-section, fp is the stress considering the total stress at strands, fc is the compressive strength; fct is the direct tensile strength; εprev is the previous strain in the prestressed strands, equals to 4.5 ‰ (for fpi = 885 MPa), fp is the stress in the strands. Table 3 presents the initial parameters obtained with Equations 13-16 using x = xlim, establishing the equilibrium of forces in the cross-section, it is possible to determine the imbalanced vector ΔR that governates the position of the neutral axis to ensure the balance until ΔR approaches zero. In this way, for each iteration of ΔR, Acc, Ast, εc, and εst are updated. The values of Rcc (force in concrete), Rsc (force in compressed reinforcements), Rst (force in tensile reinforcement), and Rct (force in tensile concrete) for the initial and final ΔR are presented in Table 3. Evaluating the neutral axis final position, it was possible to verify that the prestressed strands are working with strains equals to εst = εpnd + 10 ‰ = 15.02 ‰, following the hypothesis of domain 2. The resistant bending moment of cross-section is obtained by Equation 16: In Equation 16, Yi are the distances of gravity center of prestressed steel strands to the neutral axis of the beam; Finally, the strength bending moment of 4527 kN ⋅ m was obtained through the analytical design equations. This value is very close to the experimental value obtained by Graybeal [16] and the numerical simulation developed in this paper (4318 kN ⋅ m), with a 4.7% error. Figure 9 presents the force-displacement of I-beams (see the simulation condition in Figure 3c and 3d). It is possible to verify the high ductile behavior of the beams, showing high displaceable capacity until the total loss of strength. Figures 9b and 9c show the tensile and compressive damage distribution for the I-beam section with h = 400 mm and b = 300 mm. It can be noted a high level of tensile damage in the inferior zone at the center of the span, characterizing the bending failure mode. Nevertheless, it is also possible to detect the presence of tensile damage at the diagonals around the supports, showing the influence of the bending-shear composed failure mode, probably due to the minor span/height relation (L/h) in comparison with Graybeal's girder.

I-beams
Equations presented in Section 2 were used to predict the maximum bending moment of the designed sections. Figure 9d presents the comparison of analytical bending moments and the values obtained through the numerical simulation. There is a correspondence between the numerical and analytical results that can be achieved, showing the accuracy of the proposed model given by Fehling et al. [7], [15] and adapted by this paper.

CONCLUDING AND REMARKS
This paper approaches the numerical modeling and design of prestressed UHPFRC I-beams subjected to flexural tests. The main aspects can be highlighted: • The potential of concrete damage plasticity constitutive model was demonstrated to application in UHPFRC prestressed beams, presenting correspondence with the experimental results of Graybeal [16] for the PCI AASHTO SII beams; • The strength bending moments obtained by the analytical design equations showed to be accurate to predict the experimental and numerical results, with 5% error; • A qualitative analysis shows the high strength and satisfactory ductile behavior of beams, with high energy dissipation before the failure; • Due to the high span/height relation, the Graybeal's girders presented a typical flexural failure. For the I-beams studied, with reduced L/h coefficient, it was possible to observe a composed bending-shear failure mode; • The numerical simulation developed in this paper presented high accuracy, showing an error of 5% in the prediction of strength for Graybeal's girders and I-beams. • The simplified Graybeal model presented an error of 10%, estimating the maximum force strength. This difference evidences the latent influence of experimental variability between uniaxial sample tests and real structure behavior. The Brazilian code ABNT NBR 6118:2014 does not present the prescriptions for to design of prestressed and nonprestressed UHPFRC sections. In this way, this paper proposes a simple procedure to be applied at the design of prestressed beams subjected to bending loads.