Combination of 3D solid finite elements with rotation-free beam elements for non-linear analysis of fiber reinforced polymer rebars

<!-- metadata commented in wiki content
==Combination of 3D solid finite elements with rotation-free beam elements for non-linear analysis of fiber reinforced polymer  rebars==

'''Francisco Zárate<math>^{1}</math>,  Thibault Villette<math>^{2}</math>, Xavier Martínez<math>^{1,3}</math>,  

Fernando Rastellini<math>^{1}</math>, David Cendón<math>^{4}</math>,  Carmen Andrade<math>^{1}</math> and Eugenio Oñate<math>^{1,3}</math>

<math>^{1}</math> Centre Internacional de Metodes Numerics a l'Enginyeria (CIMNE),

Gran Capitán s/n, 08034 Barcelona, Spain

zarate@cimne.upc.edu; frastellini@cimne.upc.edu; candrade@cimne.upc.edu; 

xmartinez@cimne.upc.edu; onate@cimne.upc.edu

<math>^{2}</math> Aramco UK,  thibault.villette@aramco.com

<math>^{3}</math> Universitat Politecnica de Catalunya (UPC), Jordi Girona 1-3

08034 Barcelona, Spain

<math>^{4}</math> Centro de Investigación en Materiales Estructurales (CIME)

Technical University of Madrid (UPM), Madrid, Spain

david.cendon.franco@upm.es'''
-->
==Abstract==

We present  a combination of three dimensional (3D) solid elements and rotation-free beam elements for non-linear analysis of fiber-reinforced polymer (FRP) of rebars. The matrix material is modelled by 3D solid elements while the fibers are modelled with rotation-free 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 elastic-damage model.  The efficiency and accuracy of the combined 3D-beam 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, rotation-free beam elements, polymer fiber rebars

==1. Introduction==

There is an increasing interest in the use of non-metallic materials, such as polymers, in the building and construction (B&C) sector as an alternative to standard materials, such as steel and cement <span id='citeF-1'></span><span id='citeF-2'></span>[[#cite-1|[1,2]]]. These new materials have many advantages, compared to traditional ones, in terms of strength-weight and stiffness-weigh ratios, corrosion resistance to chemical agents, and fatigue performance, to name a few. All these advantages can extend the life-frame 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 cost-effective 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 <span id='citeF-3'></span><span id='citeF-4'></span>[[#cite-3|[3,4]]], as short fibers embedded in concrete <span id='citeF-5'></span>[[#cite-5|[5]]], or as pultruded elements which can be used as rebars for concrete reinforcement <span id='citeF-6'></span>[[#cite-6|[6]]] or as beams elements to construct different structures such as bridges <span id='citeF-7'></span>[[#cite-7|[7]]], or lighthouses <span id='citeF-8'></span>[[#cite-8|[8]]] ([[#img-1|Figures 1]]a and [[#img-1|1]]b), to name some examples.

<div id='img-1a'></div>
<div id='img-1b'></div>
<div id='img-1'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|align="center" | 
{|style="margin: 1em auto 0.1em auto;width:auto;" 
|+
|-
|style="text-align: center;padding:10px;"| [[Image:Draft_Onate_870101189-Figure1-intro-a.png|393px|]]
|style="text-align: center;padding:10px;"| [[Image:Draft_Onate_870101189-Figure1-intro-b.png|121px|]]
|-
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(a) 
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(b) 
|}
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figure 1'''. Applications of pultruded beams in construction <span id='citeF-7'></span><span id='citeF-8'></span>[[#cite-7|[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) <span id='citeF-9'></span>[[#cite-9|[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 <span id='citeF-10'></span><span id='citeF-11'></span>[[#cite-10|[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 <span id='citeF-12'></span>[[#cite-12|[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 <span id='citeF-13'></span>[[#cite-13|[13]]].

This paper presents a new computational approach for assessing the non-linear behavior of pultruded composites beams using a combination of 3D solid finite elements with rotation-free 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 rotation-free 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 non-linear 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 post-failure situations. The fact that both the 3D solid elements and the rotation-free 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 elastic-damage 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 non-linear analysis of FRP rebars combining solid elements and rotation-free 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 rotation-free beam elements based on a cell-vertex approach. The coupling of 3D-solid elements and rotation-free beam elements is detailed. The finite element formulation is validated in the non-linear 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.

==2. Combination of 3D solid elements and rotation-free elements for analysis of polymer fiber rebars==

[[#img-2|Figure 2]] depicts the modelling  strategy for analysis of FRP rebars with the FEM <span id='citeF-14'></span><span id='citeF-15'></span>[[#cite-14|[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.

<div id='img-2'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Onate_870101189-Figure1.png|600px|Finite element modelling  of FRP rebar]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''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 non-linear 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 <math display="inline">t</math>, can be written in the following quasi-static form <span id='citeF-15'></span>[[#cite-15|[15]]]

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t {\boldsymbol r} ({}^t {\boldsymbol a}):= {}^t {\boldsymbol p} ({}^t {\boldsymbol a})-{}^t {\boldsymbol f}=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

with

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t {\boldsymbol p} ({}^t {\boldsymbol a})= \int _{{}^t V} {}^t {\boldsymbol B}{}^t \boldsymbol{\sigma }dV </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

In Eqs.[[#eq-1|(1)]] and [[#eq-2|(2)]], <math display="inline">{\boldsymbol r}</math> is the vector of residual forces, <math display="inline">{\boldsymbol a}</math> is the vector of nodal displacement, <math display="inline">{\boldsymbol p}</math> is the vector of internal nodal forces, <math display="inline">{\boldsymbol B}</math> is the strain matrix relating the strains with the nodal displacements, <math display="inline">\boldsymbol{\sigma }</math> is the Cauchy stress vector, <math display="inline">{\boldsymbol f}</math> is the vector of equivalent nodal forces, <math display="inline">V</math> is the volume of the FRP bar and the upper left index <math display="inline">t</math> denotes the time <math display="inline">t</math>.

The vector of internal forces <math display="inline">{}^t {\boldsymbol p}</math> is computed by sum of the contributions from the solid elements and the beam elements, i.e.

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t {\boldsymbol p}({}^t {\boldsymbol a})= {}^t {\boldsymbol p}_s ({}^t {\boldsymbol a}) + {}^t {\boldsymbol p}_b({}^t {\boldsymbol a}) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where indices <math display="inline">s</math> and <math display="inline">b</math> denote the contribution of the 3D solid elements and the beam elements to vector <math display="inline">{}^t {\boldsymbol p}</math>, respectively.

The equilibrium solution at time <math display="inline">t +\Delta t</math> is found by solving Eq.[[#eq-1|(1)]] using a standard Newton-Raphson procedure briefly outlined below.

For each iteration of the Newton-Raphson scheme <span id='citeF-15'></span>[[#cite-15|[15]]] we perform the following computations:

Computation of displacement increments

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t {\boldsymbol K}_T \Delta {\boldsymbol a}^i = - {}^t {\boldsymbol r}^i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

with

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t {\boldsymbol r}^i = {}^t {\boldsymbol p} ({}^{t+\Delta t} {\boldsymbol a}^i) - {}^{t+\Delta t} {\boldsymbol f} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

Displacement update

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^{t+\Delta t} {\boldsymbol a}^{i+1} =  {}^{t+\Delta t} {\boldsymbol a}^i + \Delta  {\boldsymbol a}^i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

Stress update

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta \boldsymbol{\varepsilon }^i =  {}^t {\boldsymbol B}   \Delta  {\boldsymbol a}^i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta \boldsymbol{\sigma }^i =  {\boldsymbol D}_T    \Delta \boldsymbol{\varepsilon }^i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^{t+\Delta t} \boldsymbol{\sigma }^{i+1} = {}^{t+\Delta t} \boldsymbol{\sigma }^i + \Delta \boldsymbol{\sigma }^i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

Residual force vector computation

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^{t+\Delta t} {\boldsymbol r}^{i+1} =  \int _{{}^{t+\Delta t}V} {}^{t+\Delta t}{\boldsymbol B} {}^{t+\Delta t} \boldsymbol{\sigma }^{i+1}dV - {}^{t+\Delta t} {\boldsymbol f} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

Convergence check

The iteration stop when

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Vert {}^{t+\Delta t} r^{i+1}\Vert < \boldsymbol{\varepsilon } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

where <math display="inline">\boldsymbol{\varepsilon }</math> is a prescribed error value. In our computations we have taken <math display="inline"> \varepsilon = 10^{-3}</math>.

In Eq.[[#eq-4|(4)]] <math display="inline">{}^t {\boldsymbol K}_T </math> is the tangent stiffness matrix. We have computed this matrix as

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t {\boldsymbol K}_T = {}^t {\boldsymbol K}_{T_s} + {}^t {\boldsymbol K}_{T_b} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

where indices <math display="inline">s</math> and <math display="inline">b</math> denote the contribution of the solid elements and the beam elements to the tangent stiffness matrix, respectively. These contributions are computed as

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t {\boldsymbol K}_{T_r} = \int _{V_r} {}^t{\boldsymbol B}_r^T {}^t {\boldsymbol D}_{T_r} {}^t{\boldsymbol B}_r dV\, , \qquad \mbox{with }\,\,\, r=s,b </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where <math display="inline">{}^t{\boldsymbol B}_r</math> and <math display="inline">{}^t {\boldsymbol D}_{T_r} </math> with <math display="inline">r=s,b</math> denote the deformation matrix and the tangent constitutive matrix for the solid elements and the beam elements, respectively at time <math display="inline">t</math>.

The computation of these matrices is detailed in the next sections.

==3. Formulation of 3D solid elements for the transverse reinforced resin material==

===3.1 Definition of displacements, strains and stresses===

Let us consider a FRP rebar discretized using 3D-solid elements of <math display="inline">n</math> nodes. The displacement field within the element can be expressed in terms of the nodal displacements in the standard finite element form as <span id='citeF-14'></span><span id='citeF-15'></span>[[#cite-14|[14,15]]]

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol u} =[u,v,w]^T = {\boldsymbol N}^e {\boldsymbol a}^e </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

where <math display="inline">{\boldsymbol u}</math> is the displacement vector, <math display="inline">u, v </math> and <math display="inline">w</math> are the displacements along the cartesian axes <math display="inline">x,y,z</math>, respectively, <math display="inline">{\boldsymbol N}^e</math> and <math display="inline">{\boldsymbol a}^e</math> are the shape functions matrix and the nodal displacement vector for the element <math display="inline">e</math> given by

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol N}^e = [{\boldsymbol N}^e_1,{\boldsymbol N}^e_2, \cdots , {\boldsymbol N}^e_n] \quad , \quad              {\boldsymbol a}^e = \left\{\begin{matrix}{\boldsymbol a}^e_1\\{\boldsymbol a}^e_2\\\vdots \\{\boldsymbol a}^e_n\end{matrix}\right\} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol N}^e_i = N_i^e \left[\begin{matrix}1 &0&0\\ 0&1&0\\ 0&0&1\end{matrix}\right]\quad , \quad              {\boldsymbol a}^e = \left\{\begin{matrix}u_i\\v_i\\w_i\end{matrix}\right\} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

are the nodal components of matrix <math display="inline">{\boldsymbol N}^e</math> and vector <math display="inline">{\boldsymbol a}^e</math>, and <math display="inline">N_i^e</math> is the shape function matrix of node <math display="inline">i</math>.

The strains within the element can be readily expressed in terms of the nodal displacements as

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\varepsilon } ={\boldsymbol B}^e {\boldsymbol a}^e </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

where

<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\varepsilon } = \left\{\begin{matrix}\varepsilon _x\\ \varepsilon _y\\ \varepsilon _z \\ \gamma _{xy}\\\gamma _{xz}\\\gamma _{yz}   \end{matrix} \right\}= \left\{\begin{matrix}\displaystyle \frac{\partial u}{\partial x}\\ \displaystyle \frac{\partial v}{\partial y}\\[.3cm]\displaystyle \frac{\partial w}{\partial z}\\[.3cm] \displaystyle \frac{\partial u}{\partial y} + \displaystyle \frac{\partial v}{\partial x}\\[.3cm] \displaystyle \frac{\partial u}{\partial z} + \displaystyle \frac{\partial w}{\partial x}\\[.3cm] \displaystyle \frac{\partial v}{\partial z} + \displaystyle \frac{\partial w}{\partial y}\end{matrix} \right\} \quad ,\quad {\boldsymbol B}^e = [{\boldsymbol B}_1^e, {\boldsymbol B}_2^e, \cdots , {\boldsymbol B}_n^e] \quad ,\quad  {\boldsymbol B}_i = \left[\begin{matrix}\displaystyle{{\partial N_i} \over {\partial x} } & 0 & 0\\[.3cm] 0 & \displaystyle { {\partial N_i} \over {\partial y} } & 0\\[.3cm] 0 & 0 & \displaystyle{ {\partial N_i} \over {\partial z} }\\[.3cm] \displaystyle{ {\partial N_i} \over {\partial y} } & \displaystyle{ {\partial N_i} \over  {\partial x} } & 0\\[.4cm] \displaystyle{{\partial N_i} \over {\partial z} } & 0 & \displaystyle{ {\partial N_i} \over {\partial x} }\\[.4cm] 0 & \displaystyle { {\partial N_i} \over {\partial z} } & \displaystyle{ {\partial N_i} \over  { \partial y} }\\\end{matrix}\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

where <math display="inline">{\boldsymbol B}^e</math> and <math display="inline">{\boldsymbol B}^e_i</math> are the strain matrices for the element and the <math display="inline">i</math>th node, respectively.

The stress-strain relationship is written in incremental form (neglecting initial strains and initial stressed) as

<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>d\boldsymbol{\sigma } = {\boldsymbol D}_T d\boldsymbol{\varepsilon } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

where

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\sigma } = [\sigma _x, \sigma _y, \sigma _z , \tau _{xy}, \tau _{xz},\tau _{yz}]^T </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

is the Cauchy stress vector and <math display="inline">{\boldsymbol D}_T</math> is the tangent constitutive matrix. In our computations we have modelled the matrix material with an elastoplastic model based on a Drucker-Prager yield functions and an elastic-plastic-damage 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.[[#eq-2|(2)]]) and the tangent constitutive matrix.

The validation examples presented in the paper have been solved using three-linear 8-noded hexahedral elements <span id='citeF-14'></span>[[#cite-14|[14]]].

===3.2 Material model for solid elements===

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 <span id='citeF-17'></span>[[#cite-17|[17]]] and used by Martinez et al. in the simplified model presented in <span id='citeF-16'></span>[[#cite-16|[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 (<math display="inline">\hat {\boldsymbol D}</math>) can be obtained by the inverse combination of the stiffness matrix of its constituents (<math display="inline">{\boldsymbol D}_m</math> and <math display="inline">{\boldsymbol D}_f</math>), proportionally to their volumetric participation (<math display="inline">V_m</math> and <math display="inline">V_f</math>), being <math display="inline">m</math> and <math display="inline">f</math> matrix and fiber, respectively. This is

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hat{\boldsymbol D}^{-1}=V_m {\boldsymbol D}^{-1}_m + V_f {\boldsymbol D}^{-1}_f </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

The constitutive matrices <math display="inline">{\boldsymbol D}_m</math> and <math display="inline">{\boldsymbol D}_f</math> are shown in Eq.[[#eq-22|(22)]]

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol D}_i = \left[\begin{matrix}H_i & F_i & F_i & 0 & 0 &0 \\   & H_i &F_i  & 0 & 0 & 0 \\ &  &  H_i & 0 & 0 & 0 \\ &  &  &  G_i  & 0 & 0 \\  & & &  & G_i & 0\\ Sym.& & & & & G_i \end{matrix}\right]\, , \quad i=m,f </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

with

<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>H_i= \displaystyle \frac{E_i (1-\nu _i)}{(1+\nu _i)(1-2\nu _i)}\, ,  \quad  F_i = \displaystyle \frac{E_i \nu _i}{(1+\nu _i)(1-2\nu _i)} \, ,  \quad  G_i = \displaystyle \frac{(1-2\nu _i)}{2(1-\nu _i)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

We note that the  constitutive matrix <math display="inline">\hat{\boldsymbol D}</math> accounts for  the effect of the fiber matrix <math display="inline">{\boldsymbol D}_f</math> where <math display="inline">H_f</math> 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 <math display="inline">{\boldsymbol D}_f</math> and <math display="inline">{\boldsymbol D}_m</math> is needed in order to accurately reproduce the shear stiffness  in the resulting matrix <math display="inline">\hat{\boldsymbol D}</math>.  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 <span id='citeF-17'></span>[[#cite-17|[17]]]. An advantage of this approach versus the standard serial-parallel 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.

Non-linear effects in the response of the material are modelled by a simple isotropic damage model <span id='citeF-18'></span><span id='citeF-19'></span>[[#cite-18|[18,19]]]. According to this, Eq.[[#eq-21|(21)]] is modified as

<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hat{\boldsymbol D}^{-1} = V_m [{\boldsymbol D}_m(1-d_m)]^{-1} + V_f {\boldsymbol D}_f(1-d_f)]^{-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

where <math display="inline">\hat{\boldsymbol D}</math> is the damaged constitutive matrix for the solid element and the damage parameters <math display="inline">d_i (i=m,f)</math> are an scalar numbers between 0 and 1  defined as

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>d_i=\frac{2\left({\sigma }_{0_i}-{\sigma }_I\right)\left(E_i G_i-h{\sigma }_{0_i} \right)}{(2E_iG_i-h{\sigma }_{0_i}){\sigma }_I}\quad \quad i=m,f </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

In Eq.[[#eq-25|(25)]] <math display="inline">\sigma _{0_i}</math> is the elastic limit, <math display="inline">\sigma _I</math> is the maximum principal stress, <math display="inline">G_i</math> is the fracture energy, <math display="inline">E_i</math> is the Young's modulus and index <math display="inline">i</math> (<math display="inline">i=m,f</math>) refers to the matrix or fiber material, respectively. In addition, <math display="inline">h</math> 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 shear-dominated 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. <math display="inline">\hat{\boldsymbol D}\to 0</math>) 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 <math display="inline">\hat{\boldsymbol D} </math> is important in these zones. This distinct phenomenon has been modelled by limiting the maximum value of the damage parameters in Eq.[[#eq-25|(25)]] in these regions as

<span id="eq-26"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>0\le d_i\le 0.5 \quad \quad i=m,f </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

==4. Formulation of  rotation-free beam elements for the fibers==

The elimination of the rotations as nodal degrees of freedom in bending-type elements (such as beams, plates and shells) leads to the so-called ''rotation-free'' elements <span id='citeF-20'></span><span id='citeF-21'></span><span id='citeF-22'></span>[[#cite-20|[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 <span id='citeF-20'></span>[[#cite-20|[20]]]. The idea evolved in the last decades of the 20th century and lead to the formulation of different families of rotation-free elements for analysis of plates and shells. Oñate and Zárate <span id='citeF-20'></span>[[#cite-20|[20]]] developed a unified formulation for this type of rotation-free plate and shell elements combining ideas from the finite volume method with the FEM. Some of the rotation-free shell elements have enjoyed much popularity for non-linear analysis of shells structures under dynamic and impact loads and sheet stamping processes using explicit and implicit dynamic formulations <span id='citeF-23'></span>[[#cite-23|[23]]].

The use of rotation-free beam elements has been less popular in practice. The derivation of this type of elements based on Euler-Bernoulli beam theory can be found in <span id='citeF-14'></span><span id='citeF-22'></span><span id='citeF-24'></span>[[#cite-14|[14,22,24]]]. Extensions of these elements to account for shear deformation effect was reported by Zarate and Oñate <span id='citeF-21'></span>[[#cite-21|[21]]].

A distinct feature of all rotation-free 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 rotation-free Euler-Bernoulli beam elements formulated using an updated Lagrangian description and a  ''cell-vertex'' approach <span id='citeF-14'></span><span id='citeF-20'></span>[[#cite-14|[14,20]]]. We present below details of the element formulation.

===4.1 Straight cell-vertex rotation-free beam elements===

A particular class of rotation-free beam element  can be derived by computing  the curvature ''at each node'' using a  finite difference scheme. The resulting element is termed CVB (for Cell-Vertex Beam element) <span id='citeF-14'></span><span id='citeF-20'></span>[[#cite-14|[14,20]]].

For the sake of simplicity we will show the derivation of the stiffness equation of a  rotation-free straight beam element that bends in the <math display="inline">x-z</math> plane ([[#img-3|Figure 3]]). A similar process can be followed for deriving the stiffness equations for bending in the <math display="inline">x-y</math> plane. Due to the particular one-dimensional (1D) features of the rebars considered in this work, we will neglect  torsional effects in the fibers.

<div id='img-3'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Onate_870101189-Figure2.png|600px|Modelling  strategy for analysis of FRP with the FEM]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''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 <math display="inline">t</math>. 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 <math display="inline">i</math>  is  computed as

<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\kappa _i ={1\over l_i} \left[\left({\partial w\over \partial  x}\right)_{\!\scriptscriptstyle R_i }-   \left({\partial w\over \partial x}\right)_{\!\scriptscriptstyle L_i}\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

where <math display="inline">l_i={1\over 2} (l^e+ l^{e-1})</math> and subscripts <math display="inline">R_i</math> and <math display="inline">L_i</math> denote the midpoints of the elements located at the right and left of node <math display="inline">i</math> (Fig. [[#img-4|4]]). The curvature <math display="inline">\kappa _i</math> is assumed to be constant in the control domain <math display="inline">l_i</math> assigned to node <math display="inline">i</math> ([[#img-4|Figure 4]]).

<div id='img-4'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Onate_870101189-Fig1_15.png|440px|CVB element. Control domains around two nodes i and i+1]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figure 4'''. CVB element. Control domains around two nodes <math>i</math> and <math>i+1</math>
|}


The rotations <math display="inline">\left(\displaystyle {\partial w\over \partial      x}\right)_{\!\scriptscriptstyle R_i}</math> and <math display="inline">\left(\displaystyle {\partial w\over \partial x}\right)_{\!\scriptscriptstyle L_i}</math> are  expressed in terms of the nodal deflections as

<span id="eq-28"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left({\partial w\over \partial x}\right)_{\scriptscriptstyle R}={w_{\scriptscriptstyle i+1}-w_{\scriptscriptstyle i}\over l^e}\quad , \quad \left({\partial w\over \partial x}\right)_{\scriptscriptstyle L}={w_{\scriptscriptstyle i}-w_{\scriptscriptstyle i-1}\over l^{e-1}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

Substituting Eqs.[[#eq-28|(28)]] into [[#eq-27|(27)]] gives

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\kappa _i ={2\over  l^e l^{e-1}(l^e+ l^{e-1})}[l^e ,-(l^e+ l^{e-1}),l^{e-1},0] \left\{\begin{matrix}w_{\scriptscriptstyle i-1}\\ w_{\scriptscriptstyle i}\\ w_{\scriptscriptstyle  i+1}\\w_{\scriptscriptstyle i+2}\end{matrix}\right\}={\boldsymbol  B}_i \bar{\boldsymbol w}^{(e)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

where

<span id="eq-30"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol B}_i = {2\over  l^e l^{e-1}(l^e+ l^{e-1})}[l^e ,-(l^e+ l^{e-1}),l^{e-1},0]~~ ,~~ \bar{\boldsymbol w}^{(e)} = \left\{\begin{matrix}w_{\scriptscriptstyle  i-1}\\ w_{\scriptscriptstyle i}\\ w_{\scriptscriptstyle  i+1}\\w_{\scriptscriptstyle i+2}\end{matrix}\right\} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

Similarly, the curvature  at node <math display="inline">i+1</math> is found as

<span id="eq-31"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\kappa _{\scriptscriptstyle i+1}={2\over (l^e+ l^{e+1})} \left[{w_{\scriptscriptstyle i+2}-w_{\scriptscriptstyle  i+1}\over l^{e+1}}-   {w_{\scriptscriptstyle i+1}-w_{\scriptscriptstyle i}\over l^e}\right]=   {\boldsymbol B}_{\scriptscriptstyle i+1}\bar{\boldsymbol w}^{(e)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

with

<span id="eq-32"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol B}_{i+1} = {1\over l^e l^{e+1}(l^e+ l^{e+1})}[0,l^{e+1},-(l^e+ l^{e+1}),l^e] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

The internal virtual work over the element <math display="inline">e</math> is obtained by adding up the contributions from the two control domains <math display="inline">l_i</math> and <math display="inline">l_{i+1}</math> as

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta U^{(e)} = \int _{{l_i\over 2}} \delta  \kappa _i (E_f I)  \kappa _i dx + \int _{{l_{i+1}\over 2}} \delta  \kappa _{\scriptscriptstyle i+1} (E_f I) \kappa _{\scriptscriptstyle i+1}dx </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

<span id="eq-34"></span>

Substituting Eqs.([[#eq-27|27]]) and ([[#eq-31|31]]) into ([[#eq-33|33]]) and following the standard assembly process of the FEM  yields the element stiffness matrix as

<span id="eq-34.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol K}^{(e)} = {\boldsymbol K}^{(e)}_i+{\boldsymbol  K}^{(e)}_{\scriptscriptstyle i+1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34.a)
|}

where

<span id="eq-34.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}\displaystyle {\boldsymbol K}^{e}_i = \int _{{l_i\over 2}} {\boldsymbol B}_i^T (E_f I) {\boldsymbol B}_i dx = {l_i\over 2} {\boldsymbol B}_i^{T} (E_f I)^{(e)} {\boldsymbol B}_i \\[.5cm] \displaystyle {\boldsymbol K}^{e}_{\scriptscriptstyle i+1} = \int _{{l_{i+1}\over 2}} {\boldsymbol B}_{i+1}^{T} (E_f I) {\boldsymbol    B}_{i+1} dx = {l_{\scriptscriptstyle i+1}\over 2}{\boldsymbol B}_{i+1}^{T} (E_f I)^{(e)} {\boldsymbol    B}_{\scriptscriptstyle {i+1}} \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34.b)
|}

In Eqs.([[#eq-33|33]]) and ([[#eq-34|34]]) <math display="inline">E_f</math> and <math display="inline">I</math> 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 <math display="inline">E_f</math> and <math display="inline">I</math> are constant over the control domains.

Note that, as it is usual in rotation-free beam, plate and shell elements,  the stiffness matrix of an element involves the nodal deflections of the adjacent elements.

The global equivalent nodal force <math display="inline">{f}_i</math> is computed

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{f}_i={P}_{z_i} + {f}_{z_i} {(l^e + l^{e+1})\over 2} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35.a)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{f}_{z_i} ={1\over 2} \left[f_z^{(e)}+ {f}_z^{(e+1)} \right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35.b)
|}

where <math display="inline">{P}_{z_i}</math> is the external nodal force acting at node <math display="inline">i</math> and and <math display="inline"> {f}_z^{(e)}</math> and <math display="inline">{f}_z^{(e+1)}</math> are uniformly distributed loads acting in the vertical direction over elements <math display="inline">e</math> and <math display="inline">e+1</math>, respectively.

===4.2 Boundary conditions===

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 <span id='citeF-14'></span>[[#cite-14|[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 ([[#img-5|Figure 5]]a)
<span id="eq-36"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol K}^{(e)} = {\boldsymbol K}^{(e)}_{\scriptscriptstyle i+1} \quad  \hbox{if }\quad \kappa _i =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36.a)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol K}^{(e)} = {\boldsymbol K}^{(e)}_{\scriptscriptstyle i} \quad  \hbox{if }\quad \kappa _{i+1} =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36.b)
|}

<div id='img-5'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:65%;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Onate_870101189-Fig1_16.png|480px]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figure 5'''. CVB element. (a) Boundary condition of zero curvature at a free or simply supported node <math>B</math>. (b) Control domain for a clamped or symmetry node
|}


The curvature at a ''clamped or symmetry left-end node'' <math display="inline">i</math> is  computed as ([[#img-5|Figure 5]]b)

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\kappa _i ={2\over  l^e}{(w_{\scriptscriptstyle i+1}-w_{\scriptscriptstyle i})\over  l^e} ={2\over (l^e)^2} [1,-1,1,0]\left\{\begin{matrix}w_{\scriptscriptstyle i-1}\\ w_{\scriptscriptstyle i}\\ w_{\scriptscriptstyle i+1}\\ w_{\scriptscriptstyle i+2}\end{matrix}\right\}={\boldsymbol B}_i \bar {\boldsymbol w}^{(e)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

where <math display="inline">w_{\scriptscriptstyle i-1}</math> is an auxiliary fictitious deflection.

Similarly, for a right-end clamped or symmetry node ([[#img-5|Figure 5]]b)

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\kappa _{\scriptscriptstyle i+1} = {2\over  l^e}{(w_{\scriptscriptstyle i+1}-w_{\scriptscriptstyle i})\over  l^e}={2\over    (l^e)^2}[0,-1,1,1]\left\{\begin{matrix}w_{\scriptscriptstyle i-1}\\ w_{\scriptscriptstyle i}\\ w_{\scriptscriptstyle i+1}\\ w_{\scriptscriptstyle  i+2}\end{matrix}\right\}={\boldsymbol B}_{\scriptscriptstyle  i+1}\bar {\boldsymbol w}_i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

where <math display="inline">w_{\scriptscriptstyle i+2}</math> is an auxiliary fictitious deflection.

The element stiffness matrix  is computed in both cases by Eq.[[#eq-36|(36)]] with <math display="inline">{\boldsymbol B}_i</math> or <math display="inline">{\boldsymbol B}_{\scriptscriptstyle i+1}</math> equal to zero, as appropriate. The fictitious nodal deflections <math display="inline">w_{i-1}</math> or <math display="inline">w_{i+2}</math> (and indeed <math display="inline">w_i</math>) are prescribed to a zero value when solving the global system of equations.

The above rotation-free formulation can be extended to account for the bending behavior of the beam in the <math display="inline">y-z</math> plane in a straightforward manner. The nodal degrees of freedom are the vertical displacement along the <math display="inline">y</math> and <math display="inline">z</math> axes <span id='citeF-14'></span>[[#cite-14|[14]]].

The extension of the bending formulation of rotation-free beam element to account for transverse shear deformation effect is possible. The details can be found in <span id='citeF-21'></span>[[#cite-21|[21]]].

===4.3 Material model for the rotation-free beam element===

The material behavior of the rotation-free CVB element is simply characterized by the Young modulus. Accounting for non-linear effects in the fiber material under axial stresses has been modelled with a standard one-dimensional elastic-damage constitutive law,  i.e.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hat E_f = E_f (1-d_f^b)</math>
|}
|}

where <math display="inline">\hat E_f</math> and <math display="inline">E_f</math> are the Young modulus of the fiber in the undamaged and damaged states, respectively.

The evolution of the damage parameter <math display="inline">d_f^b</math> for the beam element is governed the maximum deformation of the fiber at failure, i.e.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>d_f^b =1 + \frac{\sigma _{0_f} (\varepsilon - 2\varepsilon _{0_f} + \varepsilon _{u_f})}{\varepsilon E_f (\varepsilon _{0_f} + \varepsilon _{u_f})}</math>
|}
|}

where <math display="inline">\varepsilon </math> is the axial deformation of the fiber, <math display="inline">\varepsilon _{0_f}</math> is the deformation at the elastic limit, <math display="inline">\sigma _{0_f}</math>, and <math display="inline">\varepsilon _{u_f}</math> is the deformation at the ultimate (fracture) strength value.

===4.4 Extension to curved rotation-free CVB elements===

A curved rebar will be modelled by a collection of straight rotation-free CVB elements. [[#img-6|Figure 6]] shows a curved rebars under bending loads in a plane <math display="inline">x-z</math>, discretized with three CVB elements.

<div id='img-6'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:70%;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Onate_870101189-Fig9_26.png|480px]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figure 6'''. Cell-vertex rotation-free  beam element. Control domains for computing the average  curvature. A, B and C are element mid-points
|}


The stiffness matrix is obtained as the sum of the axial stiffness matrix (<math display="inline">{\boldsymbol K}_A</math>) and the bending stiffness matrix (<math display="inline">{\boldsymbol K}_B</math>).

====4.4.1 Axial stiffness matrix====

The axial strain matrix <math display="inline">{\boldsymbol B}_a</math> is obtained as

<span id="eq-39"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\varepsilon _a = \frac{\partial u'}{\partial x'}= \lambda _s= {\partial u' \over \partial x'} = \frac{1}{l^e}[-1,1] \left\{\begin{matrix}{u'}^{e}_{i}\\ {u'}^{e}_{i+1}\end{matrix} \right\}=\underbrace{\frac{1}{l^e} [{\boldsymbol 0},-{\boldsymbol t}^e,{\boldsymbol t}^e,{\boldsymbol 0}]}_{{\boldsymbol B}_{a}}{\boldsymbol a}^e =  {\boldsymbol B}_{a} {\boldsymbol a}^e </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

with <math display="inline">{\boldsymbol t}^e = [\cos \phi ^e ,\sin \phi ^e]</math> and <math display="inline">()'</math> denotes variables in the local coordinates system.

The axial stiffness matrix is computed as

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol K}_A^{e} =(E_f I)^{(e)} {\boldsymbol B}_a^T {\boldsymbol B}_a l^e   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

====4.4.2 Bending stiffness matrix====

The procedure for computing the bending stiffness matrix is an extension of that followed for  straight CVB elements in Section 4.2. <span id="eq-41"></span>

The  curvature at node <math display="inline">i</math> is computed in terms of the gradients of the normal deflection at points <math display="inline">A</math> and <math display="inline">B</math> ([[#img-6|Figure 6]])

<span id="eq-41.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\kappa _{i}= \frac{2}{l^e + l^{e-1}}\left[\left({\partial w' \over \partial x'}\right)_B-\left({\partial w' \over \partial x'}\right)_A \right]= \frac{2}{l^e + l^{e-1}}\left[\frac{{w'}^{e}_{i+1}{-w_i'}^e}{l^e} - \frac{{w'}^{e-1}_{i}-w_{i-1'}^{e-1}}{l^{e-1}}\right]=</math>
|-
| style="text-align: center;" | <math> =\underbrace{\frac{1}{l^e l^{e-1}(l^e +l^{e-1})}[l^{e},-l^{e},-l^{e-1},l^{e-1}]}_ {{\bar{\boldsymbol B}'}_{b_i}}\left\{\begin{matrix}{w'}^{e-1}_{i-1}\\{w'}^{e-1}_{i} \\ {w'}^{e}_{i+1}\\{ w'}^{e}_{i+1} \end{matrix}\right\}={\boldsymbol B'}_{b_i} \mathbf{w}^{'}={\boldsymbol B}_{b_i}{\boldsymbol a}^e </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41.a)
|}

In Eqs.([[#eq-41.a|41.a]]) <math display="inline">{w'}_j</math> are the nodal deflections in the local axes of the element defined in [[#img-7|Figure 7]] and

<span id="eq-41.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol B}_{b_i}={\boldsymbol B'}_{b_i}{\boldsymbol T} \quad \hbox{with  }~ {\boldsymbol T}= \left[\begin{matrix}{\boldsymbol n}^{e-1} & {\boldsymbol 0} &  {\boldsymbol 0}& {\boldsymbol 0}\\ {\boldsymbol 0} &{\boldsymbol n}^{e-1} & {\boldsymbol 0} & {\boldsymbol 0}\\ {\boldsymbol 0} & {\boldsymbol 0}& {\boldsymbol n}^{e} &{\boldsymbol 0}\\ {\boldsymbol 0} & {\boldsymbol 0}& {\boldsymbol 0}& {\boldsymbol n}^{e} \end{matrix} \right]  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41.b)
|}

and

<span id="eq-41.c"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol a}^e = \left\{\begin{matrix}{\boldsymbol a}_{i-1}\\{\boldsymbol a}_i\\{\boldsymbol a}_{i+1}\\{\boldsymbol a}_{i+2} \end{matrix} \right\}\quad \hbox{with}\quad {\boldsymbol a}_i=\left\{\begin{matrix}u_i\\ w_i \end{matrix} \right\} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41.c)
|}

<div id='img-7'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Onate_870101189-Figure3.png|600px]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figure 7'''. Local coordinate system of curved rotation-free CVB element bending in the <math>xz</math> plane
|}


Similarly, for node <math display="inline">i+1</math>

<span id="eq-42.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\kappa _{i+1}= \frac{2}{l^e + l^{e+1}}\left[\left({\partial w' \over \partial x'}\right)_C- \left({\partial w' \over \partial x'}\right)_B \right]= \frac{2}{l^e + l^{e+1}}\left[\frac{{w'}^{e+1}_{i+2}{-w'}^{e}_{i+1}}{l^{e+1}} - \frac{{w'}^{e}_{i+1}-w_{i'}^{e}}{l^{e}}\right]=</math>
|-
| style="text-align: center;" | <math> =\underbrace{\frac{2}{l^e l^{e+1}(l^e +l^{e+1})}[l^{e+1},-l^{e+1},-l^{e},l^{e}]}_ {\bar{\boldsymbol B}_{b_{i+1}}}\left\{\begin{matrix}{w'}^{e}_i\\{w'}^{e}_{i+1} \\ {w'}^{e+1}_{i+1}\\{ w'}^{e+1}_{i+2}\end{matrix}\right\}={\boldsymbol B'}_{b_{i+1}} \mathbf{w}^{'}= {\boldsymbol B}_{b_{i+1}}{\boldsymbol a}^e  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42.a)
|}

with

<span id="eq-42.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol B}_{b_{i+1}}={\boldsymbol B'}_{b_{i+1}}{\boldsymbol T}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42.b)
|}

and

<span id="eq-42.c"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol T}= \left[\begin{matrix}{\boldsymbol n}^{e-1} & {\boldsymbol 0} & {\boldsymbol 0}& {\boldsymbol 0}\\ {\boldsymbol 0} &{\boldsymbol n}^{e-1} &{\boldsymbol 0}& {\boldsymbol 0}\\ {\boldsymbol 0} & {\boldsymbol 0}&{\boldsymbol n}^{e+1} &{\boldsymbol 0}\\ {\boldsymbol 0} & {\boldsymbol 0}& {\boldsymbol 0}&{\boldsymbol n}^{e+1} \end{matrix} \right]\quad ,\quad {\boldsymbol n}^e =[-\sin \phi ^e, \cos \phi ^e]\quad ,\quad {\boldsymbol 0}=[0,0] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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 ([[#img-6|Figure 6]]), i.e.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol K}^{e}_B =  l^e\Bigl[{\boldsymbol B}_i^T (E_f I) {\boldsymbol B}_i  + {\boldsymbol B}_{i+1}^T (E_f I){\boldsymbol B}_{i+1} \Bigr] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

The total stiffness matrix for the element is finally obtained by

<span id="eq-44"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol K}^e = {\boldsymbol K}_A^e + {\boldsymbol K}_B^e </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

where <math display="inline">{\boldsymbol K}_A</math> is given by Eq.[[#eq-40|(40)]].

We note that the size of <math display="inline">{\boldsymbol K}^{(e)}</math> for the CVB element bending in the <math display="inline">x-z</math> plane is <math display="inline">8\times 8</math>, as it involves the two DOFs of nodes <math display="inline">i-1,i,i+1</math> and <math display="inline">i+2</math>. A similar procedure will be followed for the case of bending in the plane <math display="inline">y-z</math> <span id='citeF-14'></span>[[#cite-14|[14]]].

The assembly process follows the general rule of the FEM <span id='citeF-14'></span>[[#cite-14|[14]]].

==5. Coupling of 3D solid elements and rotation-free beam elements==

Fig. [[#img-8|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 (8-noded hexahedral are shown). The fibers contained within each solid element are assumed to be concentrated in one-dimensional (1D) packs of fibers distributed within the 3D solid elements. Fig. [[#img-8|8]] shows the concentration of fibers in four 1D packs placed at the <math display="inline">2\times 2</math> Gauss points locations in each cross section <span id='citeF-14'></span>[[#cite-14|[14]]]. Each fiber pack is modelled as a rotation-free beam element with a cross-section area equivalent to the amount of fibers considered in the pack.

The displacements of the rotation-free beam element nodes are related to the nodal displacement of the 3D solid element by a simple interpolation.

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Onate_870101189-Figure4.png|600px|Patch of 3D solid and rotation-free beam elements]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Patch of 3D solid and rotation-free beam elements
|}

The discretized equilibrium equations for the coupled solid-beam system can be written as explained in Section [[#2 Combination of 3D solid elements and rotation-free elements for analysis of polymer fiber rebars|2]]. The numerical response of the FRP bars for different incremental load up to failure is computed using the standard Newton-Raphson iterative procedure described in Eqs.[[#eq-4|(4)]]-[[#eq-13|(13)]].

==6. Examples of application==

===6.1 Composite material properties===

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 three-point 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 ECR-Glass and a matrix of DERAKANE<math display="inline">\_</math>8084<math display="inline">\_</math>Epoxy<math display="inline">\_</math>Vinyl<math display="inline">\_</math>ester.

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%

===6.2 Axial tension test===

Fig. [[#img-9|9]] shows the bar modelled with 960 8-noded 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.

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Onate_870101189-Figure_8.png|540px|Tensile test geometry: 3D solid and RFB elements mesh]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Tensile test geometry: 3D solid and RFB elements mesh
|}

Fig. [[#img-10|10]] shows a typical bar tested at UPM after failure.

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Onate_870101189-Fig6_2.png|270px|Tensile test of the PFR performed at UPM]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Tensile test of the PFR performed at UPM
|}

The numerical results obtained are shown in  Figs. [[#img-11|11]] and [[#img-12|12]]. The displacement field obtained presents a linear distribution, until the rebar reaches brittle failure.

<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Onate_870101189-Fig6_3.png|600px|Displacement field [m] in deformed geometry]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Displacement field [m] in deformed geometry
|}

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Onate_870101189-Fig6_4.png|574px|Axial tension tests stress-strain curves. Numerical and experimental values]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' Axial tension tests stress-strain curves. Numerical and experimental values
|}

The numerical load-displacement curve obtained is converted into a stress-strain curve to be compared with the experimental results. A good correlation can be seen between the numerical and experimental values (Fig. [[#img-12|12]]).

Table [[#table-1|1]] shows the good correlation between the numerical and experimental stiffness and strength magnitudes.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 Tensile test. Comparison of numerical and experimental results
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  '''Property''' 
| style="border-left: 2px solid;border-right: 2px solid;" | '''Lab tests''' 
| style="border-left: 2px solid;border-right: 2px solid;" | '''Numerical''' 
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  Tensile strength (MPa) 
| style="border-left: 2px solid;border-right: 2px solid;" | 1150 <math>\boldsymbol{\pm }</math> 34  (MPa) 
| style="border-left: 2px solid;border-right: 2px solid;" | 1173.33 (MPa) 
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  Strain to max strength 
| style="border-left: 2px solid;border-right: 2px solid;" | 2.3 <math>\pm </math> 0.1 (%) 
| style="border-left: 2px solid;border-right: 2px solid;" | 2.22 (%) 
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  Young Modulus (GPa) 
| style="border-left: 2px solid;border-right: 2px solid;" | 55 <math>\pm </math> 4 (GPa) 
| style="border-left: 2px solid;border-right: 2px solid;" | 53.03 (GPa) 

|}

===6.3 Bending test===

The experiment performed for measuring the bending properties, corresponds to the standard three-point bending test showed in Fig. [[#img-14|14]]. The tests were carried out on bars of full circular section, instead of using cut specimens as indicated in the ASTMD-4476 norm. Fig. [[#img-13|13]] shows the bar modelled with 5400 8-noded hexahedral elements (6327 nodes) and 6156 rotation-free 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 Fig. [[#img-13|13]]. To avoid problems with the stress concentration at the loading points, a small part of the beam (the blue zone in Fig. [[#img-13|13]]) has been defined with an unbreakable material. The elastic properties of this material are the same as for the rest of the bar.

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Onate_870101189-Fig6_5.png|600px|Three point bending test for PFR rebar. Dimensions, FEM mesh and boundary conditions]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Three point bending test for PFR rebar. Dimensions, FEM mesh and boundary conditions
|}

Fig. [[#img-14|14]] shows a snapshot of the laboratory test for the FRP rebar in a three-point bending configuration. Note the rollers used as a pin-supports of the beam.

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Onate_870101189-Fig6_6.png|534px|Three-point bending test of FRP rebars performed at UPM]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14:''' Three-point bending test of FRP rebars performed at UPM
|}

Fig. [[#img-15|15]] shows the strain distribution over the deformed bar   (Fig. [[#img-15|15]]a) and the mesh of hexahedra (Fig. [[#img-15|15]]b) in a load state close to the onset of damage. The damage distribution for each element is displayed in Fig. [[#img-16|16]].

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Onate_870101189-Fig6_7.png|600px|Three-point bending tests. Numerical results. Deformed bar. (a) Strain distribution for  RFB elements. (b) Strain distribution in  3D solid elements]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Three-point bending tests. Numerical results. Deformed bar. (a) Strain distribution for  RFB elements. (b) Strain distribution in  3D solid elements
|}

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Onate_870101189-Figure_15.png|600px|Contour field of damage variable in longitudinal and transverse sections on the deformed geometry (1x)]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Contour field of damage variable in longitudinal and transverse sections on the deformed geometry (1x)
|}

In Fig. [[#img-17|17]], the numerical load-deflection 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.

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. Fig. [[#img-18|18]] shows the test set-up and one  specimen after failure.

<div id='img-17'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Onate_870101189-Fig6_9.png|586px|Three-point bending tests. Force-deflection curves for shearing testing. Numerical results compared with laboratory tests]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' Three-point bending tests. Force-deflection curves for shearing testing. Numerical results compared with laboratory tests
|}

<div id='img-18'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Onate_870101189-Fig6_10.png|600px|Three-point bending tests. Failed specimen of PFR rebar]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 18:''' Three-point bending tests. Failed specimen of PFR rebar
|}

===6.4 Shear test===

The experiment carried out for evaluating the shear properties, corresponds to the ASTM-7617 standard test showed in Fig. [[#img-19|19]] <span id='citeF-25'></span>[[#cite-25|[25]]]. For the numerical analysis, the symmetry of the geometry and the imposed displacement have been taken into account as shown in Fig. [[#img-20|20]]a in order to speed up the calculations. Fig. [[#img-20|20]]b displays the bar modelled with  5225 8-noded hexahedral elements (4608 nodes) and 5016 rotation-free beam elements. The length of the sample is 100 mm long with 9.58 mm diameter.

<div id='img-19'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"
|-
| style="padding:10px;" | [[Image:Draft_Onate_870101189-Fig6_11a.png|200px|]]
| style="padding:10px;" | [[Image:Draft_Onate_870101189-Fig6_11b.png|176px|]]
|-
| colspan="2" style="padding:10px;" |[[Image:Draft_Onate_870101189-Fig6_11c.png|350px|Shear test.  Double shear fixture  according to ASTM D-7617 standard <span id='citeF-25'></span>[[#cite-25|[25]]]]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 19:''' Shear test.  Double shear fixture  according to ASTM D-7617 standard <span id='citeF-25'></span>[[#cite-25|[25]]]
|}

<div id='img-20'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"
|-
|style="padding:10px;" |[[Image:Draft_Onate_870101189-Fig6_12a.png|600px|]]
|style="padding:10px;" |[[Image:Draft_Onate_870101189-Fig6_12b.png|600px|Shear test geometry. (a) Boundary conditions. (b) mesh and dimensions]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 20:''' Shear test geometry. (a) Boundary conditions. (b) mesh and dimensions
|}

In Fi. [[#img-20|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.

Fig. [[#img-21|21]] shows a representative configuration of the numerical results. Fig. [[#img-21|21]]a displays the deformed bar and the localized damage in the resin material at the shear zone. Fig. [[#img-21|21]]b is a close-up of the shear zone were the RFB are presented. The curvature of these elements due to the shear deformation is depicted.

Fig. [[#img-22|22]] shows the numerical load-deflection curve superimposed to a set of curves obtained in the laboratory tests for comparison purposes. In the curves the non-linear 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 non-linear, proves the validity of the proposed procedure to simulate pultruded fiber reinforced polymer rebars.

<div id='img-21'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Onate_870101189-Fig6_13.png|600px|Shear test. Contour field of damage variable in the matrix. Deformed geometry (1x)]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 21:''' Shear test. Contour field of damage variable in the matrix. Deformed geometry (1x)
|}

<div id='img-22'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Onate_870101189-Figure_21.png|600px|Shear test. Force-deflection curves. Numerical results compared with laboratory tests]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 22:''' Shear test. Force-deflection curves. Numerical results compared with laboratory tests
|}

==7. Concluding remarks==

The combination of 3D solid finite elements with rotation-free 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 3D-1D 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 rotation-free beam elements can be implemented in more sophisticated formulations for modelling fiber reinforced materials, such as that proposed in <span id='citeF-26'></span>[[#cite-26|[26]]].

==Data Availability Statement==

Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

==Acknowledgements==

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” (CEX2018-000797-S).

==References==

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''  Pacheco-Torgal F., Jalali S., Labrincha J., John V.M. Eco-efficient concrete. Woodhead Publishing Series in Civil and Structural Engineering. 2013 ISBN: 9780857098993.

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]'''  Bank L.C. Composites for construction: Structural design with FRP materials. John Wiley & Sons, Inc.2006. ISBN 9780471681267

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]'''   Marí A.R., Oller E.,  Bairán J.M. Predicting the response of FRP-strengthened reinforced-concrete flexural members with nonlinear evolutive analysis models. J. Compos. Constr., 15(5), 799–809, 2011.

<div id="cite-4"></div>
'''[[#citeF-4|[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(15-16), 1604-1618, 2008. doi:10.1016/j.compstruc.2008.01.007

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''  Nogales A.,  Tosić N.,  de la Fuente A. Rotation and moment redistribution capacity of fiber-reinforced concrete beams: Parametric analysis and code compliance. Struct. Concr., 23(1), 220-239, 2021.

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]'''  Pultron Composites. [Online]. Available: https://www.pultron.com/products/features/. [Accessed: 17-Apr-2020].

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]'''  Lleida GFRP Pedestrian Bridge. [Online]. Available: http://www.pedelta.com/lleida-gfrp-pedestrian-bridge-p-52-en. [Accessed: 17-Apr-2020].

<div id="cite-8"></div>
'''[[#citeF-8|[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.

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]'''   GangaRao H.V.S.,  Taly N.,  Vijay P.V. Reinforced concrete design with FRP composites. 1st Edition, CRC Press, 2019.

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]'''  Sharifianjazi F., Zeydi P.,  Bazli M.,  Esmaeilkhanian A., Rahmani R., Bazli  L.,  Khaksar S. Fibre-reinforced polymer reinforced concrete members under elevated temperatures: a review on structural performance.  Polymers, 14(3), 472, 2022. https://doi.org/10.3390/polym14030472

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]'''  Sagar B.,  Sivakumar M.V.N.  Performance evaluation of basalt fibre-reinforced polymer rebars in structural concrete members–a review. Innovative Infrastructure Solutions, 6(2), 1-18, 2021. https://doi.org/10.1007/s41062-020-00452-2

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]'''   Strau S.,  Senz A.,  Ellinger J., Comparison of the processing of epoxy resins in pultrusion with open bath impregnation and closed-injection pultrusion. J. Compos. Sci., 3(3), 87, 2019.

<div id="cite-13"></div>
'''[[#citeF-13|[13]]]'''  Starr T.F. Pultrusion for engineers. Woodhead Publishing Series in Composites Science and Engineering, pp. 336, 2000. ISBN: 978-1-85573-425-8

<div id="cite-14"></div>
'''[[#citeF-14|[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. Springer-CIMNE.

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]'''  O.C Zienkiewicz and R.L. Taylor. ''The Finite Element Method for Solids and Structural Mechanics''. 6th edition, Elsevier, 2005

<div id="cite-16"></div>
'''[[#citeF-16|[16]]]'''  Martínez X., Rastellini F., Oller S., Flores F. and Oñate E.  Computationally optimized formulation for the simulation of composite materials and delamination failures. Composites Part B: Engineering,  42(2), 134-144, 2011. 10.1016/j.compositesb.2010.09.013

<div id="cite-17"></div>
'''[[#citeF-17|[17]]]'''  Rastellini F., Oller O., Salomón O., Oñate E. Composite materials non-linear modelling for long fibre-reinforced laminates: Continuum basis, computational aspects and validations. ''Computers & Structures'', Vol. 86(9), 879&#8211;896, 2008. https://doi.org/10.1016/j.compstruc.2007.04.009.

<div id="cite-18"></div>
'''[[#citeF-18|[18]]]'''  Oller S., Oñate E., Oliver J. and Lubliner J.  Finite element non-linear analysis of concrete structures using a plastic-damage model. Int. Journal of Engineering Fracture Mechanics,  35(1/2/3), 219-231, 1990. 10.1016/0013-7944(90)90200-Z

<div id="cite-19"></div>
'''[[#citeF-19|[19]]]'''  Lubliner J., Oller S., Oliver J. and Oñate E.  A plastic damage model for concrete. Int. Journal of Solids and Structures,  25(3), 299-326, 1989. 10.1016/0020-7683(89)90050-4

<div id="cite-20"></div>
'''[[#citeF-20|[20]]]'''  Oñate E., Zárate F. Rotation-free triangular plate and shell elements. ''Int. J. Numer. Meth. Engng.'', 47(1-3), 557&#8211;603, 2000.

<div id="cite-21"></div>
'''[[#citeF-21|[21]]]'''  Oñate E., Zárate F. Extended rotation-free plate and beam elements with shear deformation effects. ''Int. J. Numer. Meth. Engng.'', Vol. 82(2), 196&#8211;227, 2010.

<div id="cite-22"></div>
'''[[#citeF-22|[22]]]'''  Jovicevic J. and Oñate E.  Analysis of beams and shells using a rotation - free finite element - finite volume formulation. Monograph M43,  CIMNE, ISBN: 84-89925-36-4, Barcelona, 1999.    Analysis of beams and shell using a rotation - free finite element - finite volume formulation - Jovicevic Onate 2018a - Scipedia.

<div id="cite-23"></div>
'''[[#citeF-23|[23]]]'''  Oñate E. and Flores F.G.  Advances in the formulation of the rotation-free basic shell triangle. Computer Methods in Applied Mechanics and Engineering,  194 (21-24), 2406-2443, 2005. 10.1016/j.cma.2004.07.039

<div id="cite-24"></div>
'''[[#citeF-24|[24]]]'''  Flores F. and Oñate E.  Rotation-free finite element for the non-linear analysis of beams and axisymmetric shells. Computer Methods in Applied Mechanics and Engineering,  195(41-43), 5297-5315, 2006. 10.1016/j.cma.2005.08.021

<div id="cite-25"></div>
'''[[#citeF-25|[25]]]'''  ASTM-International, ASTM D7617/D7617M - 11(2017) Standard Test Method for Transverse Shear Strength of Fiber-reinforced Polymer Matrix Composite Bars, West Conshohocken, PA, 2017.

<div id="cite-26"></div>
'''[[#citeF-26|[26]]]'''  Khristenko U., Schu S., Krüger M., Schmidt F., Wohlmuth B. and Hesch C. Multidimensional coupling: A variationally consistent approach to fiber-reinforced materials. Computer Methods in Applied Mechanics and Engineering, 382:113869, 2021. 10.1016/j.cma.2021.113869