==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]]] (Figs. [[#img-1|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="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Onate_870101189-Figure1-intro-a.png|393px|]]

|[[Image:Draft_Onate_870101189-Figure1-intro-b.png|121px|]]

|- style="text-align: center; font-size: 75%;"

| (a) 

| (b) 

|- style="text-align: center; font-size: 75%;"

| colspan="2" | '''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==

​

Fig. [[#img-2|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="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Onate_870101189-Figure1.png|600px|Finite element modelling  of FRP rebar]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''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 (Fig. [[#img-3|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="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Onate_870101189-Figure2.png|600px|Modelling  strategy for analysis of FRP with the FEM]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''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> (Fig. [[#img-4|4]]).

​

<div id='img-4'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Onate_870101189-Fig1_15.png|440px|CVB element. Control domains around two nodes i and i+1]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''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]]].

​

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

​

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 (Fig. [[#img-5|5]]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+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="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Onate_870101189-Fig1_16.png|600px|CVB element. (a) Boundary condition of zero curvature at a free or   simply supported node B. (b) Control domain for a clamped or symmetry node]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''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 (Fig. [[#img-5|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 (Fig. [[#img-5|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. Fig. [[#img-6|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="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Onate_870101189-Fig9_26.png|570px|Cell-vertex rotation-free  beam element. Control domains for computing the average  curvature. A, B and C are element mid-points]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''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> (Fig. [[#img-6|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 Fig. [[#img-7|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="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Onate_870101189-Figure3.png|600px|Local coordinate system of curved rotation-free CVB element bending in the xz plane]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''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>