We present a combination of three dimensional (3D) solid elements and rotationfree beam elements for nonlinear analysis of fiberreinforced polymer (FRP) of rebars. The matrix material is modelled by 3D solid elements while the fibers are modelled with rotationfree beam elements. The absence of rotation variables in the beam elements allows the straightforward coupling of 3D solid and beam elements using a formulation with displacement nodal degrees of freedom only. Both solid and beam elements are formulated in an updated Lagrangian description. The behavior of the matrix and the fiber material are modelled with an elasticdamage model. The efficiency and accuracy of the combined 3Dbeam element formulation are verified in examples of application to the analysis of FRP rebars up to fracture in axial, bending and shear tests for which experimental results are available.
Keywords: 3D solid finite elements, rotationfree beam elements, polymer fiber rebars
There is an increasing interest in the use of nonmetallic materials, such as polymers, in the building and construction (B&C) sector as an alternative to standard materials, such as steel and cement [1,2]. These new materials have many advantages, compared to traditional ones, in terms of strengthweight and stiffnessweigh ratios, corrosion resistance to chemical agents, and fatigue performance, to name a few. All these advantages can extend the lifeframe of structures and minimize the maintenance required, thus contributing to improve the sustainability of the sector. Current research efforts aim to develop resistant, sustainable and costeffective materials incorporating polymers for a number of applications in building and construction.
Polymers and composite materials can be used in B&C in many forms: as external reinforcements of concrete beams and slabs [3,4], as short fibers embedded in concrete [5], or as pultruded elements which can be used as rebars for concrete reinforcement [6] or as beams elements to construct different structures such as bridges [7], or lighthouses [8] (Figures 1a and 1b), to name some examples.
 
Figure 1. Applications of pultruded beams in construction [7,8] 
The use of pultruded FRP as reinforcement in concrete has several advantages compared to steel reinforcement such as corrosion resistance, high longitudinal (unidirectional) strength (2 or 3 times higher than steel for equivalent bar diameter), high fatigue endurance (which depends on the type of reinforcing fiber and bar), magnetic transparency, lightweight (25% of steel weight) and low thermal and electric conductivity (for glass and aramid fibers) [9]. Glass fibers are the common choice due to their good performance and reduced cost, compared to other synthetic fibers that have improved properties but are more expensive, such as carbon. Research efforts are currently carried out for developing new families of rebars with synthetic fibers, as reinforcement for concrete structures [10,11]. This type of synthetic fiber rebars could be an alternative to current FRP rebars using glass fibers, and even to steel rebars for specific applications.
Pultrusion beams are obtained by pulling the reinforcement fibers, which have previously impregnated with a matrix resin, through a heated die where the composite is cured [12]. The shape of the die defines the shape of the final beam produced. Due to the fabrication process, the final composite contains mainly fibers aligned with the beam axis, and the efficiency of the method allows having fiber participations that can reach 80% of the total composite weight content [13].
This paper presents a new computational approach for assessing the nonlinear behavior of pultruded composites beams using a combination of 3D solid finite elements with rotationfree beam elements. These last ones are included in the model to characterize the longitudinal performance of the composite.
The formulation presented is of special interest for pultruded beams, as in these elements fibers are mainly oriented in the longitudinal direction, which simplifies the distribution of the rotationfree beams in the numerical model. As it will be shown, the adaptability of the beam elements in the model geometry allows characterizing complex failure mechanisms such as shear failure, which depends highly on the tensile forces taken by the fibers once the matrix has failed. To do so, the geometrical nonlinear effects in the solid and beam elements are taken into account by using an updated Lagrangian description, allowing for large displacements effects that can occur in specific loading cases, such as in shear tests and in postfailure situations. The fact that both the 3D solid elements and the rotationfree beam element are formulated in terms of displacement degrees of freedom (DOFs) only facilitates the integration of both types of elements in a unified formulation. An elasticdamage model is used for characterizing the behavior of both the matrix and the fibers.
The content of the paper is as follows. In the next section we present the overall finite element approach for nonlinear analysis of FRP rebars combining solid elements and rotationfree beam elements. Then we describe the 3D solid finite element formulation for modelling the behavior of the matrix material. This is followed by the finite element formulation for the fibers using rotationfree beam elements based on a cellvertex approach. The coupling of 3Dsolid elements and rotationfree beam elements is detailed. The finite element formulation is validated in the nonlinear analysis up to failure of FRP bar specimens formed by an epoxy matrix and carbon fibers tested under axial, bending and shear loads. Numerical results are compared with experimental results for the same tests. Excellent agreement is obtained in the comparison.
Figure 2 depicts the modelling strategy for analysis of FRP rebars with the FEM [14,15]. The procedure proposed dissociates the longitudinal contribution of fibers to the composite strength and stiffness, from their contribution in the transversal and shear directions. This is done by modeling a transverse reinforced resin, that accounts for the resin and the transverse and shear performance of the fibers, with standard 3D solid elements; and by modeling the longitudinal contribution of fibers to the composite with a collection of beams embedded in the 3D solid model.
Figure 2. Finite element modelling of FRP rebar 
Both the 3D solid elements and the beam elements are formulated using an updated Lagrangian description in order to account for geometrically nonlinear effects in the deformation of the FRP bars. Following standard finite element procedures, the discretized equilibrium equation for the bar structure, expressed in the current configuration at time , can be written in the following quasistatic form [15]

(1) 
with

(2) 
In Eqs.(1) and (2), is the vector of residual forces, is the vector of nodal displacement, is the vector of internal nodal forces, is the strain matrix relating the strains with the nodal displacements, is the Cauchy stress vector, is the vector of equivalent nodal forces, is the volume of the FRP bar and the upper left index denotes the time .
The vector of internal forces is computed by sum of the contributions from the solid elements and the beam elements, i.e.

(3) 
where indices and denote the contribution of the 3D solid elements and the beam elements to vector , respectively.
The equilibrium solution at time is found by solving Eq.(1) using a standard NewtonRaphson procedure briefly outlined below.
For each iteration of the NewtonRaphson scheme [15] we perform the following computations:
Computation of displacement increments

(4) 
with

(5) 
Displacement update

(6) 
Stress update

(7) 

(8) 

(9) 
Residual force vector computation

(10) 
Convergence check
The iteration stop when

(11) 
where is a prescribed error value. In our computations we have taken .
In Eq.(4) is the tangent stiffness matrix. We have computed this matrix as

(12) 
where indices and denote the contribution of the solid elements and the beam elements to the tangent stiffness matrix, respectively. These contributions are computed as

(13) 
where and with denote the deformation matrix and the tangent constitutive matrix for the solid elements and the beam elements, respectively at time .
The computation of these matrices is detailed in the next sections.
Let us consider a FRP rebar discretized using 3Dsolid elements of nodes. The displacement field within the element can be expressed in terms of the nodal displacements in the standard finite element form as [14,15]

(14) 
where is the displacement vector, and are the displacements along the cartesian axes , respectively, and are the shape functions matrix and the nodal displacement vector for the element given by

(15) 
where

(16) 
are the nodal components of matrix and vector , and is the shape function matrix of node .
The strains within the element can be readily expressed in terms of the nodal displacements as

(17) 
where

(18) 
where and are the strain matrices for the element and the th node, respectively.
The stressstrain relationship is written in incremental form (neglecting initial strains and initial stressed) as

(19) 
where

(20) 
is the Cauchy stress vector and is the tangent constitutive matrix. In our computations we have modelled the matrix material with an elastoplastic model based on a DruckerPrager yield functions and an elasticplasticdamage model. The damage model is described in the next section.
The expressions of the strain matrix and the Cauchy stress vector are used for computing the nodal forces (Eq.(2)) and the tangent constitutive matrix.
The validation examples presented in the paper have been solved using threelinear 8noded hexahedral elements [14].
The material model defined for the solid elements has to reproduce the performance of what has been called transverse reinforced resin. This is, a resin with improved properties in its transversal and shear directions, due to the presence of fibers. The transversal directions are defined as those directions perpendicular to the fiber longitudinal orientation.
As stated by Rastellini et al. when defining the serial parallel mixing theory [17] and used by Martinez et al. in the simplified model presented in [16], the contribution of the fibers to the composite in the transversal and shear directions can be captured by applying the inverse mixing theory. This states that the constitutive matrix of the composite () can be obtained by the inverse combination of the stiffness matrix of its constituents ( and ), proportionally to their volumetric participation ( and ), being and matrix and fiber, respectively. This is

(21) 
The constitutive matrices and are shown in Eq.(22)

(22) 
with

(23) 
We note that the constitutive matrix accounts for the effect of the fiber matrix where has a null value in the fiber direction (in practice a small value is chosen to avoid numeric instability). The stiffness in the fiber direction is provided by the RFB elements.
The distinct form of matrices and is needed in order to accurately reproduce the shear stiffness in the resulting matrix . In this manner, the combination of solid elements and RFB elements can reproduce a parallel mixed material behavior along the longitudinal direction of the RFB elements (deformations in the fiber direction are compatible between all materials) and a serial mixed material behavior in the orthogonal direction to the fibers (stresses in the normal fiber direction are compatible between all materials), according to the inverse mixing theory [17]. An advantage of this approach versus the standard serialparallel material mixing model is the compatibility between the stresses in the concrete and the fibers in the orthogonal direction to the fibers, which facilitates the convergence of the nonlinear iterative process.
Nonlinear effects in the response of the material are modelled by a simple isotropic damage model [18,19]. According to this, Eq.(21) is modified as

(24) 
where is the damaged constitutive matrix for the solid element and the damage parameters are an scalar numbers between 0 and 1 defined as

(25) 
In Eq.(25) is the elastic limit, is the maximum principal stress, is the fracture energy, is the Young's modulus and index () refers to the matrix or fiber material, respectively. In addition, is the characteristic length of the element.
The above constitutive model can accurately reproduce axial and pure bending stress states in which the RFB elements are mainly subjected to tensile forces. The model has proved to work well in sheardominated regions in which two different stress situations appear in the shear zone and in the region where the loads are applied. In the shear zone the matrix material is soon fully damaged (i.e. ) and the stiffness is mainly contributed by the fibers as they rotate and change their orientation. In zones where the loads are applied in a direction orthogonal to the fibers, the fibers are agglomerated in packages that allow the transmission of compression stresses to the resin material, inducing damage. The contribution of matrix is important in these zones. This distinct phenomenon has been modelled by limiting the maximum value of the damage parameters in Eq.(25) in these regions as

(26) 
The elimination of the rotations as nodal degrees of freedom in bendingtype elements (such as beams, plates and shells) leads to the socalled rotationfree elements [20,21,22]. Rotations can be simply eliminated by expressing the curvature field in a node, or a point within an element, in terms of the displacements of nodes contained in an appropriate patch of elements surrounding the node or the point. This procedure emanates from the efforts of finite difference practitioners for analysis of thin plates based on Kirchhoff plate theory using displacement variables only [20]. The idea evolved in the last decades of the 20th century and lead to the formulation of different families of rotationfree elements for analysis of plates and shells. Oñate and Zárate [20] developed a unified formulation for this type of rotationfree plate and shell elements combining ideas from the finite volume method with the FEM. Some of the rotationfree shell elements have enjoyed much popularity for nonlinear analysis of shells structures under dynamic and impact loads and sheet stamping processes using explicit and implicit dynamic formulations [23].
The use of rotationfree beam elements has been less popular in practice. The derivation of this type of elements based on EulerBernoulli beam theory can be found in [14,22,24]. Extensions of these elements to account for shear deformation effect was reported by Zarate and Oñate [21].
A distinct feature of all rotationfree elements is that the discretized equilibrium equations can be expressed in terms of displacement degrees of freedom only. This makes these elements advantageous for their coupling with solid elements, that are also formulated in terms of nodal displacement DOFs, in a straightforward manner. In this work we will model the mechanical behavior of the fibers with rotationfree EulerBernoulli beam elements formulated using an updated Lagrangian description and a cellvertex approach [14,20]. We present below details of the element formulation.
A particular class of rotationfree beam element can be derived by computing the curvature at each node using a finite difference scheme. The resulting element is termed CVB (for CellVertex Beam element) [14,20].
For the sake of simplicity we will show the derivation of the stiffness equation of a rotationfree straight beam element that bends in the plane (Figure 3). A similar process can be followed for deriving the stiffness equations for bending in the plane. Due to the particular onedimensional (1D) features of the rebars considered in this work, we will neglect torsional effects in the fibers.
Figure 3. Modelling strategy for analysis of FRP with the FEM 
Let us formulate the deformation of a straight beam element for a current configuration at time . The bending behavior of the beam is modelled by considering the control domain formed by half of the lengths of the elements adjacent to a node. The curvature at node is computed as

(27) 
where and subscripts and denote the midpoints of the elements located at the right and left of node (Fig. 4). The curvature is assumed to be constant in the control domain assigned to node (Figure 4).
Figure 4. CVB element. Control domains around two nodes and 
The rotations and are expressed in terms of the nodal deflections as

(28) 
Substituting Eqs.(28) into (27) gives

(29) 
where

(30) 
Similarly, the curvature at node is found as

(31) 
with

(32) 
The internal virtual work over the element is obtained by adding up the contributions from the two control domains and as

(33) 
Substituting Eqs.(27) and (31) into (33) and following the standard assembly process of the FEM yields the element stiffness matrix as

(34.a) 
where

(34.b) 
In Eqs.(33) and (34) and are the Young modulus of the fiber material and the moment of inertia of the fiber section with respect the bending axis, respectively. In the computation of the above integrals, we have assumed that and are constant over the control domains.
Note that, as it is usual in rotationfree beam, plate and shell elements, the stiffness matrix of an element involves the nodal deflections of the adjacent elements.
The global equivalent nodal force is computed

(35.a) 
with

(35.b) 
where is the external nodal force acting at node and and and are uniformly distributed loads acting in the vertical direction over elements and , respectively.
The prescribed values for the deflection are imposed when solving the global system of equations, as usual. We briefly explain next how to compute the curvature at a free end and at simply supported and clamped nodes. For details see [14].
The curvature at a free end node and at a simply supported node is zero. This condition is implemented by neglecting the contribution of the boundary node to the element stiffness matrix. For a boundary element with a free or simply supported node the stiffness matrix is (Figure 5a)

(36.a) 

(36.b) 
Figure 5. CVB element. (a) Boundary condition of zero curvature at a free or simply supported node . (b) Control domain for a clamped or symmetry node 
The curvature at a clamped or symmetry leftend node is computed as (Figure 5b)

(37) 
where is an auxiliary fictitious deflection.
Similarly, for a rightend clamped or symmetry node (Figure 5b)

(38) 
where is an auxiliary fictitious deflection.
The element stiffness matrix is computed in both cases by Eq.(36) with or equal to zero, as appropriate. The fictitious nodal deflections or (and indeed ) are prescribed to a zero value when solving the global system of equations.
The above rotationfree formulation can be extended to account for the bending behavior of the beam in the plane in a straightforward manner. The nodal degrees of freedom are the vertical displacement along the and axes [14].
The extension of the bending formulation of rotationfree beam element to account for transverse shear deformation effect is possible. The details can be found in [21].
The material behavior of the rotationfree CVB element is simply characterized by the Young modulus. Accounting for nonlinear effects in the fiber material under axial stresses has been modelled with a standard onedimensional elasticdamage constitutive law, i.e.

where and are the Young modulus of the fiber in the undamaged and damaged states, respectively.
The evolution of the damage parameter for the beam element is governed the maximum deformation of the fiber at failure, i.e.

where is the axial deformation of the fiber, is the deformation at the elastic limit, , and is the deformation at the ultimate (fracture) strength value.
A curved rebar will be modelled by a collection of straight rotationfree CVB elements. Figure 6 shows a curved rebars under bending loads in a plane , discretized with three CVB elements.
Figure 6. Cellvertex rotationfree beam element. Control domains for computing the average curvature. A, B and C are element midpoints 
The stiffness matrix is obtained as the sum of the axial stiffness matrix () and the bending stiffness matrix ().
The axial strain matrix is obtained as

(39) 
with and denotes variables in the local coordinates system.
The axial stiffness matrix is computed as

(40) 
The procedure for computing the bending stiffness matrix is an extension of that followed for straight CVB elements in Section 4.2.
The curvature at node is computed in terms of the gradients of the normal deflection at points and (Figure 6)

(41.a) 
In Eqs.(41.a) are the nodal deflections in the local axes of the element defined in Figure 7 and

(41.b) 
and

(41.c) 
Figure 7. Local coordinate system of curved rotationfree CVB element bending in the plane 
Similarly, for node

(42.a) 
with

(42.b) 
and

(42.c) 
The bending stiffness matrix for the CVB element is expressed as the sum of the stiffness contributions from the two control subdomains 1 and 2 that form the element (Figure 6), i.e.

(43) 
The total stiffness matrix for the element is finally obtained by

(44) 
where is given by Eq.(40).
We note that the size of for the CVB element bending in the plane is , as it involves the two DOFs of nodes and . A similar procedure will be followed for the case of bending in the plane [14].
The assembly process follows the general rule of the FEM [14].
Figure 8 shows a conceptual representation of a FRP bar of rectangular section containing of a number of fibers embedded in the (bulk) matrix material. The bar geometry is discretized into 3D solid elements (8noded hexahedral are shown). The fibers contained within each solid element are assumed to be concentrated in onedimensional (1D) packs of fibers distributed within the 3D solid elements. Figure 8 shows the concentration of fibers in four 1D packs placed at the Gauss points locations in each cross section [14]. Each fiber pack is modelled as a rotationfree beam element with a crosssection area equivalent to the amount of fibers considered in the pack.
The displacements of the rotationfree beam element nodes are related to the nodal displacement of the 3D solid element by a simple interpolation.
Figure 8. Patch of 3D solid and rotationfree beam elements 
The discretized equilibrium equations for the coupled solidbeam system can be written as explained in Section 2. The numerical response of the FRP bars for different incremental load up to failure is computed using the standard NewtonRaphson iterative procedure described in Eqs.(4)(13).
The three examples shown hereafter demonstrate the good behavior of the formulation described for analysis of FRP rebars. The examples include the analysis of a FRP rebar in three standard laboratory tests: an axial test, a threepoint bending test and a shear test. The numerical results are compared with results from laboratory tests carried out at the experimental lab for structural materials of the Centro de Investigación en Materiales Estructurales (CIME) of the Technical University of Madrid (UPM), using rebars made of ECRGlass and a matrix of DERAKANE8084EpoxyVinylester.
The rebar nominal diameter is 10mm, including helicoidal ribs. As these ribs are manufactured by machining the originally pultruded rebars, the effective diameter of the bars is reduced to 9.58 mm, on average. This diameter is the one adopted for creating the model of the bar geometries.
The mechanical parameters of the fibers are: Young modulus: 74.61 GPa, Poisson ratio: 0.22, Elastic limit (tension): 1,653 MPa, Elastic limit (compression): 1,322 MPa, Fracture energy: 126,778 J/m2, Deformation at failure: 0.0277, Volume fraction: 70%
The mechanical parameters of the matrix material are: Young modulus: 3.17 GPa, Poisson coefficient: 0.38, Elastic limit: 72 MPa, Fracture energy: 45,657 J/m2, Deformation at failure: 0.11, Volume fraction 30%
Figure 9 shows the bar modelled with 960 8noded hexahedral elements and 1197 nodes. The fibers have been modelled using 1140 RFB elements. The length of the sample is 100 mm long with 9.58 mm diameter. Symmetry conditions have been considered by prescribing a zero axial displacement in the normal direction to the surface at a bar end, while the other end has a prescribed axial displacement condition to induce tensile stresses in the bar. The bar axis has the displacements in the normal plane yz constrained to zero.
Figure 9. Tensile test geometry: 3D solid and RFB elements mesh 
Figure 10 shows a typical bar tested at UPM after failure.
Figure 10. Tensile test of the PFR performed at UPM 
The numerical results obtained are shown in Figures 11 and 12. The displacement field obtained presents a linear distribution, until the rebar reaches brittle failure.
Figure 11. Displacement field [m] in deformed geometry 
Figure 12. Axial tension tests stressstrain curves. Numerical and experimental values 
The numerical loaddisplacement curve obtained is converted into a stressstrain curve to be compared with the experimental results. A good correlation can be seen between the numerical and experimental values (Figure 12).
Table 1 shows the good correlation between the numerical and experimental stiffness and strength magnitudes.
Property  Lab tests  Numerical 

Tensile strength (MPa)  1150 34 (MPa)  1173.33 (MPa) 
Strain to max strength  2.3 0.1 (%)  2.22 (%) 
Young Modulus (GPa)  55 4 (GPa)  53.03 (GPa) 
The experiment performed for measuring the bending properties, corresponds to the standard threepoint bending test showed in Figure 14. The tests were carried out on bars of full circular section, instead of using cut specimens as indicated in the ASTMD4476 norm. Figure 13 shows the bar modelled with 5400 8noded hexahedral elements (6327 nodes) and 6156 rotationfree beam elements. The length of the bar sample is 100 mm long with 9.58 mm diameter. Symmetry and boundary conditions have been considered, as shown in Figure 13. To avoid problems with the stress concentration at the loading points, a small part of the beam (the blue zone in Figure 13) has been defined with an unbreakable material. The elastic properties of this material are the same as for the rest of the bar.
Figure 13. Three point bending test for PFR rebar. Dimensions, FEM mesh and boundary conditions 
Figure 14 shows a snapshot of the laboratory test for the FRP rebar in a threepoint bending configuration. Note the rollers used as a pinsupports of the beam.
Figure 14. Threepoint bending test of FRP rebars performed at UPM 
Figure 15 shows the strain distribution over the deformed bar (Figure 15a) and the mesh of hexahedra (Figure 15b) in a load state close to the onset of damage. The damage distribution for each element is displayed in Figure 16.
Figure 15. Threepoint bending tests. Numerical results. Deformed bar. (a) Strain distribution for RFB elements. (b) Strain distribution in 3D solid elements 
Figure 16. Contour field of damage variable in longitudinal and transverse sections on the deformed geometry (1x) 
In Figure 17, the numerical loaddeflection curve is shown together with that obtained in the laboratory for comparison purposes. The maximum loading force and the flexural stiffness captured by the numerical simulation are in good agreement with the experimental results.
Figure 17. Threepoint bending tests. Forcedeflection curves for shearing testing. Numerical results compared with laboratory tests 
At the laboratory, samples of 240mm of total length were cut for the 10mm diameter full cylinder tests. The span of the tests was set to 200mm. Tests were carried out under displacements control, with a crosshead displacement velocity of 3mm/min. Figure 18 shows the test setup and one specimen after failure.
Figure 18. Threepoint bending tests. Failed specimen of PFR rebar 
The experiment carried out for evaluating the shear properties, corresponds to the ASTM7617 standard test showed in Figure 19 [25]. For the numerical analysis, the symmetry of the geometry and the imposed displacement have been taken into account as shown in Figure 20a in order to speed up the calculations. Figure 20b displays the bar modelled with 5225 8noded hexahedral elements (4608 nodes) and 5016 rotationfree beam elements. The length of the sample is 100 mm long with 9.58 mm diameter.
 
Figure 19. Shear test. Double shear fixture according to ASTM D7617 standard [25] 
 
Figure 20. Shear test geometry. (a) Boundary conditions. (b) Mesh and dimensions 
In Figure 20 a gap of 0.001 m can be observed between the loaded and supported zones of the specimen. This zone is necessary to capture accurately the shear forces and to avoid stress concentration in the numerical model.
Figure 21 shows a representative configuration of the numerical results. Figure 21a displays the deformed bar and the localized damage in the resin material at the shear zone. Figure 21b is a closeup of the shear zone were the RFB are presented. The curvature of these elements due to the shear deformation is depicted.
Figure 22 shows the numerical loaddeflection curve superimposed to a set of curves obtained in the laboratory tests for comparison purposes. In the curves the nonlinear load phenomena of the fiber reorientation can be clearly observed. For an applied load slightly larger than 5kN there is an initial plateau in the experimental tests, also reproduced by the numerical model, that corresponds to the cracking of the matrix. At this point there is a fiber reorientation that provides the stiffness increase and the final composite failure, due to fiber failure, for a load close to 35kN. The good agreement obtained between the numerical and the experimental results, especially in this case, highly nonlinear, proves the validity of the proposed procedure to simulate pultruded fiber reinforced polymer rebars.
Figure 21. Shear test. Contour field of damage variable in the matrix. Deformed geometry (1x) 
Figure 22. Shear test. Forcedeflection curves. Numerical results compared with laboratory tests 
The combination of 3D solid finite elements with rotationfree beam elements has shown to be adequate for modelling axial, bending and shear deformation states in fiber reinforced polymer rebars. The fact that both types of elements are formulated in terms of displacement DOFs only simplifies the combination of the best features of both elements for modelling the 3D behavior of the matrix material and the 1D behavior of the fibers, as well as the coupling interactions during complex deformation modes.
The good agreement of the numerical results for several benchmark laboratory tests with experimental data for the same tests have, validated the accuracy of the combined 3D1D model presented. The combined model can be therefore be used confidently in practical applications of the fiber polymer rebars as reinforcement of concrete structures.
The approach presented in this work for combining 3D solid elements with rotationfree beam elements can be implemented in more sophisticated formulations for modelling fiber reinforced materials, such as that proposed in [26].
Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.
This research was carried out as a part of the project “Development of polymer fiber reinforced rebars” funded by Saudi Aramco. The authors also acknowledge the financial support from the Spanish Ministry of Economy and Competitiveness, through the “Severo Ochoa Programme for Centres of Excellence in R&D” (CEX2018000797S).
[1] PachecoTorgal F., Jalali S., Labrincha J., John V.M. Ecoefficient concrete. Woodhead Publishing Series in Civil and Structural Engineering, pp. 624, 2013.
[2] Bank L.C. Composites for construction: Structural design with FRP materials. John Wiley & Sons, Inc., pp. 560, 2006.
[3] Marí A.R., Oller E., Bairán J.M. Predicting the response of FRPstrengthened reinforcedconcrete flexural members with nonlinear evolutive analysis models. J. Compos. Constr., 15(5):799–809, 2011.
[4] Martínez X., Oller S., Rastellini F., Barbat A.H. A numerical procedure for simulating RC structures reinforced with FRP using the serial/parallel mixing theory. Computers & Structures, 86(1516):16041618, 2008. doi:10.1016/j.compstruc.2008.01.007
[5] Nogales A., Tosić N., de la Fuente A. Rotation and moment redistribution capacity of fiberreinforced concrete beams: Parametric analysis and code compliance. Struct. Concr., 23(1):220239, 2021.
[6] Pultron Composites [Online]. Available: https://www.pultron.com/products/features/. [Accessed: 17Apr2020].
[7] Lleida GFRP Pedestrian Bridge [Online]. Available: http://www.pedelta.com/lleidagfrppedestrianbridgep52en. [Accessed: 17Apr2020].
[8] Hernández S., et al. Design and development of an innovative structural system based on composites applied on the constrction of maritme lighthouses with lower maintenance requeriments (FAROCOMP). Mater. Compuestos, 3(2):54–58, 2019.
[9] GangaRao H.V.S., Taly N., Vijay P.V. Reinforced concrete design with FRP composites. 1st Edition, CRC Press, pp. 400, 2019.
[10] Sharifianjazi F., Zeydi P., Bazli M., Esmaeilkhanian A., Rahmani R., Bazli L., Khaksar S. Fibrereinforced polymer reinforced concrete members under elevated temperatures: A review on structural performance. Polymers, 14(3):472, 2022. 10.3390/polym14030472
[11] Sagar B., Sivakumar M.V.N. Performance evaluation of basalt fibrereinforced polymer rebars in structural concrete members–a review. Innovative Infrastructure Solutions, 6(2):118, 2021. 10.1007/s41062020004522
[12] Strau S., Senz A., Ellinger J. Comparison of the processing of epoxy resins in pultrusion with open bath impregnation and closedinjection pultrusion. J. Compos. Sci., 3(3):87, 2019.
[13] Starr T.F. Pultrusion for engineers. Woodhead Publishing Series in Composites Science and Engineering, pp. 336, 2000.
[14] Oñate E. Structural analysis with the finite element method. Linear statics. Vol 1: Basis and solids, 2009. Vol. 2: Beams, plates and shells, 2013, SpringerCIMNE.
[15] Zienkiewicz O.C., Taylor R.L. The finite element method for solids and structural mechanics. 6th edition, Elsevier, 2005.
[16] Martínez X., Rastellini F., Oller S., Flores F., Oñate E. Computationally optimized formulation for the simulation of composite materials and delamination failures. Composites Part B: Engineering, 42(2):134144, 2011. 10.1016/j.compositesb.2010.09.013
[17] Rastellini F., Oller O., Salomón O., Oñate E. Composite materials nonlinear modelling for long fibrereinforced laminates: Continuum basis, computational aspects and validations. Computers & Structures, 86(9):879896, 2008. 10.1016/j.compstruc.2007.04.009
[18] Oller S., Oñate E., Oliver J., Lubliner J. Finite element nonlinear analysis of concrete structures using a plasticdamage model. Int. Journal of Engineering Fracture Mechanics, 35(1/2/3):219231, 1990. 10.1016/00137944(90)90200Z
[19] Lubliner J., Oller S., Oliver J., Oñate E. A plastic damage model for concrete. Int. Journal of Solids and Structures, 25(3):299326, 1989. 10.1016/00207683(89)900504
[20] Oñate E., Zárate F. Rotationfree triangular plate and shell elements. Int. J. Numer. Meth. Engng., 47(13):557603, 2000. <557::AIDNME784>3.0.CO;29 10.1002/(SICI)10970207(20000110/30)47:1/3<557::AIDNME784>3.0.CO;29
[21] Oñate E., Zárate F. Extended rotationfree plate and beam elements with shear deformation effects. Int. J. Numer. Meth. Engng., 83(2):196227, 2010. 10.1002/nme.2836
[22] Jovicevic J., Oñate E. Analysis of beams and shells using a rotation  free finite element  finite volume formulation. Monograph M43, CIMNE, ISBN: 8489925364, Barcelona, 1999.
[23] Oñate E., Flores F.G. Advances in the formulation of the rotationfree basic shell triangle. Computer Methods in Applied Mechanics and Engineering, 194(2124):24062443, 2005. 10.1016/j.cma.2004.07.039
[24] Flores F., Oñate E. Rotationfree finite element for the nonlinear analysis of beams and axisymmetric shells. Computer Methods in Applied Mechanics and Engineering, 195(4143):52975315, 2006. 10.1016/j.cma.2005.08.021
[25] ASTMInternational, ASTM D7617/D7617M  11(2017) Standard Test Method for Transverse Shear Strength of Fiberreinforced Polymer Matrix Composite Bars, West Conshohocken, PA, 2017.
[26] Khristenko U., Schuß S., Krüger M., Schmidt F., Wohlmuth B., Hesch C. Multidimensional coupling: A variationally consistent approach to fiberreinforced materials. Computer Methods in Applied Mechanics and Engineering, 382:113869, 2021. 10.1016/j.cma.2021.113869
Published on 04/05/23
Accepted on 04/05/23
Submitted on 16/03/23
Volume 39, Issue 2, 2023
DOI: 10.23967/j.rimni.2023.05.001
Licence: CC BYNCSA license
Are you one of the authors of this document?