Non-Linear Three Dimensional Finite Elements for Composite Concrete Structures

The current investigation focused on the development of effective and suitable modelling of reinforced concrete component with and without strengthening. The modelling includes physical and constitutive models. New interface elements have been developed, while modified constitutive law have been applied and new computational algorithm is utilised. The new elements are the Trusslink element to model the interaction between concrete and reinforcement bars, the interface element between two plate bending elements and the interface element to represent the interfacial behaviour between FRP, steel plates and concrete. Nonlinear finite-element (FE) codes were developed with pre-processing. The programme was written using FORTRAN language. The accuracy and efficiency of the finite element programme were achieved by analyzing several examples from the literature. The application of the 3D FE code was further enhanced by carrying out the numerical analysis of the three dimensional finite element analysis of FRP strengthened RC beams, as well as the 3D non-linear finite element analysis of girder bridge. Acceptable distributions of slip, deflection, stresses in the concrete and FRP plate have also been found. These results show that the new elements are effective and appropriate to be used for structural component modelling.


INTRODUCTION
In the analysis of the composite structures, modelling of various structural elements interface (such as beam, truss, plate, membrane, solid element, etc.) should be carefully taken into account.Otherwise, accuracy will not be achieved.Moreover, these issues cannot be handled by using any of the above-mentioned elements.Interface translation and rotational degrees of freedom should be considered.Introducing suitable interface elements will solve this problem for specific behaviour.
Attempt was made by Ahmad (1987) to develop a new finite element called Ahmad's link to represent the interface between bar and concrete.However, the development and application of this particular element in RC component need multi-level generation of the link element but only a few reports are available.Hence, further improvement in the development and generation of this type of element is required.
Interface elements were used to model the bonding between FRP plates and concrete beams by Lu et al. (2006), Niu and Wu (2006) and Abdel Baky et al. (2007).In addition, the interface elements were also used to model the bonding between the FRP plates and concrete components such as slabs (Elsayed et al., 2007), curved concrete beams (Jalili & Chao, 2010) and masonry structures (Waleed et al., 2008;Milani, 2010).Ebead and Marzouk (2004) simulated the behaviour of the concrete-FRP interface using a very fine mesh to model the adhesive layer that is defined as a linear elastic material.However, Hu et al. (2004), Lundquist et al. (2005), as well as Santhakumar and Chandrasekaran (2004) who have studied the behaviour of retrofitted structures ignored the effects of the interfacial behaviour between the FRP plates and concrete for difficulty in the modelling of interface behaviour.Rodrigues et al.(2015) were modeled the interaction between steel bars and concrete through the use of interface finite elements with high aspect ratio and a damage model designed to describe the bond-slip behavior which proposed by Manzoli et al. (2012).
More recently Vinh (2014) presented an open source one/two dimensional programme to generate zero-thickness cohesive interface elements.Ali Zokaei et al. (2014) proposed a 2D plate element for the idealization of composite slabs.Interface elements are used to connect the unbonded layers within the composite slab and to model the interface between the bottom slab layer and the foundation layer.Khelifa et al. (2015) focused on the flexural behaviour of timber beams externally reinforced using Carbon Fibre-Reinforced Plastics (CFRP).A cohesive model in the Abaqus software was used to represent the interaction between two adherent surfaces (CFRP and timber).The mesoscale model for concrete was detailed by Eduardo et al. (2016) Attention is given to the strategy used for modelling the coarse aggregates and the presentation of the mesh fragmentation technique.The formulation of the interface solid finite element and the continuum tension damage model used to describe its behavior.
Meanwhile, non-linear analysis of the reinforced concrete slab-girder bridge has been carried out by several researchers (e.g.Song et al., 2002;Tianyu et al., 2005;Chansawat et al., 2006;Sulata, 2007, Tong Guo et al. 2016).Most of the investigations have focussed on some specific issues, such as cracking, crushing and yielding, with different elements and computational procedures.Moreover, the development and application of the interface elements, such as Truss-link or the interface between plate bending element and brick element to analyse (static or dynamic) the RC structures, have not been well reported in the literature.Hence, there is a need to combine all these important parameters in research programmes.
The current investigation focused on the development of effective and suitable modelling of reinforced concrete component with and without strengthening.The modelling includes physical and constitutive models.New interface elements have been developed, while modified constitutive law have been applied and new computational algorithm is utilized.The application of the 3D FE code was further enhanced by carrying out the numerical analysis of the, three dimensional finite element analysis of FRP strengthened RC beams, as well as the 3D non-linear finite element analysis of girder bridge.Acceptable distributions of slip, deflection, shear and axial stresses in the FRP plate have also been found.These results show that the new elements are effective and appropriate to be used for structural component modelling.

Choice of the Element
In reinforced concrete material the external load is applied to the concrete and the reinforcing bars only receive its part of the load from the surrounding concrete by bond.In the composite structures, the bond between different components of reinforced concrete member has a fundamental role and its neglect is conducted to poor structural response.In the following sections, suitable elements for modelling the concrete, bars, FRP or steel plate and new interface elements are chosen and developed into a realistic model of RC structures with and without strengthening.
The detail for modelling of each element proposed for finite element analysis and physical modelling of Reinforced concrete beam is discussed briefly.

Concrete Element
To represent modelling the concrete, twenty seven node Lagrangian brick element (Figure 1), is a powerful element for 3D finite element modelling of concrete.The element has three degrees of freedom in each node.The element is from lagrangian family.The shape functions for the 27 node Lagrangian brick element are generated using zero method and products of three Lagrange polynomials which reported by Kohnehpooshi et al (2010).
where: u , v , w are displacements in x, y and z directions and θ , θ are rotations in the x and y directions.This element, can allow for the transverse shear deformations.Rotations, θ and θ can be expressed as: where: ϕ , ϕ are average shear rotations of the mid surface normal.More detail and formulation of this element has already been presented by Viladkar and co. Workers (1994).Figure 2 shows the mentioned element.Truss-link element which presented by Kohnehpooshi et al. (2010) is an element that assembled three dimensional truss and linkage elements as a one element as shown in Figure 4.The two linkage elements, each one with two nodes will be assembled with one 3D truss element.Subsequently, the element has four nodes with 3 DOF in each node.Nodes i and j are bar element nodes, and each pair of nodes i, i' and j, j' are linkage element nodes as shown in Figure 4.Only one spring element in the direction of the slip is considered.In the other two directions the same displacements for concrete and steel are used with selecting high values of Young's Modulus (E) in these direc-

Vertical spring
Lateral spring Z' (r, s,t) Latin American Journal of Solids and Structures 14 (2017) 398-421 tions, then concrete and steel have same deflections in these directions.More detail and formulation of this element has already been presented by Kohnehpooshi et al. (2010).

Plate Bending Interface Element
This newly plate bending interface element was developed to connect two plate bending elements.It has 16 nodes with 5 DOF in each node.This element connects 8 nodes of the top plate and the bottom plates with identical coordinates, as shown in Figure 5. Physically, the element does not exist, but its mechanical action for each node is represented by 5 springs which include 3 orthogonal springs that are connected in the horizontal, vertical and lateral directions, as well as 2 rotational springs to the plate elements.The horizontal and lateral spring represents the bond stiffness and also acts as a bond between two adherends material.The vertical and two rotational bonds can also be applied if the bond-slip relation of the vertical deformation and two rotations can be found for the adhesive material.The explicit form of K, for each link node is given as follow.
, , , , , , , : are direction cosines of the local axis with respect to the global axes , , , , : are the bond-slip modules in the five directions.These are obtained using an idealized form of the bond-slip curves between the two plates: where b   is the incremental bond stress and S  is the slip from a specified bond-slip curve.
Eq. 3 is the stiffness matrix of the link element to connect the top and bottom continues in just one node.However, there are 8 nodes in the top and the bottom of interface elements.Therefore, the stiffness matrix for the plate interface element includes 8 link elements stiffness matrices that will become 80×80 stiffness matrix, as follows: More detail of this element has already been presented by Kohnehpooshi and co. Workers (2010b).

The Interface Element Between the Brick and Plate Bending Elements
This interface element was developed to connect the plate and brick elements together such as simulation of the bond between FRP or steel plate and concrete elements, and it has 16 nodes, as shown in Figure 6.It connects 8 nodes of the top brick element and 8 nodes of the bottom plate with identical coordinates.Meanwhile, its mechanical action for the nodes adjacent to the plate is represented by 5 springs, including 3 orthogonal springs that are connected in the horizontal, vertical and lateral directions and 2 rotational springs linked to the plate's elements.The mechanical action for the nodes adjacent to the brick is represented by 3 orthogonal springs that are connected in the horizontal, vertical and lateral directions.Derivation of the stiffness matrix follows the step used in deriving the plate bending interface element.The difference lies only in the degrees of freedom of the nodes that are adjacent to the brick element.More detail and formulation of this element has already been presented by Kohnehpooshi et al. (2011).

Brick Interface Element
Brick interface element represents the bond modelling between two brick elements.Once again, the interface element has 16 nodes but with 3 DOF in each node.Moreover, there is no rotation degree of freedom.This interface element can also be used to model the interface between the membrane shell elements with 3 translations DOF.More detail and formulation of this element has already been presented by Kohnehpooshi (2012).

Constitutive Modelling for Bond-Link, Concrete, Steel and CFRP Elements
A few constitutive models of bond stress-slip between the bar and concrete had been developed in the past, as for example: Eligehausen et al. (1983) and Harjli et al. (1995) Non-linear behaviour of the RC components is modelled by taking into account the material nonlinearity of the reinforced concrete structures, which is mainly associated with its compressive yielding.The concrete is assumed to behave according to an elasto-plastic isotropic constitutive relationship under compression.The reinforcing steel is assumed to be elastic-perfectly plastic in tension and compression.The tensile behaviour of the CFRP is assumed as linear elastic until failure.

Constitutive Modelling for the Bond-Slip Between Concrete and FRP Plates
A few models can be found for the bond-slip model between the FRP plates and concrete in the works by several researchers such as Lu et al. (2005) and Sayed-Ahmed (2009).One of the most accurate bond stress-slip is proposed by Lu et al. (2005).In their study, the behaviour of the FRP/concrete interface was simulated as a relationship between the local shear stress, τ and the relative displacement, s.Three different bond-slip relations have been suggested and are referred to as: (i) precise, (ii) simplified and (iii) bilinear models.In this study, the bilinear model shown in Figure 8 was adopted for its simplicity and little difference in the result as compared to the other two models, as shown in Lu et al. (2005).If max is the maximum bond stress and s0 is the corresponding slip, then:

Ottosen Four-Parameter Model Formulation
A failure criterion proposed by Ottosen (1977), known as the four-parameter model, was proposed in this study for formulation and implementation of the developed code as a yield criterion.This failure surface contains all the three stress invariants with the following characteristics (Bangash, 2001).The surface is smooth and convex with curved meridians, as determined by the constants a and b.It is open in the negative direction of the hydrostatic axis.The trade in the deviatoric plane changes from an almost triangular to a circular shape with the increasing hydrostatic pressure, as defined by λ on the deviatoric plane.The surface is in good agreement with the experimental results over a wide range of stress states, including those where tensile stresses occur.Ottosen defines the analytical failure surface containing all the above characteristics in the following form: , , 1 0 (9) Uniaxial compressive cylinder strength for concrete Uniaxial tensile strength for concrete As for a>0, b>0, the meridians become curvy, smooth and convex, while the surface opens in a negative direction of the hydrostatic axis.When a=0 and λ=constant, this criterion becomes closer to the Drucker-Prager criterion.When a=0, b=0 and λ=constant, it represents the von Mises criterion

COMPUTATIONAL PROCEDURES
The standard finite element procedure has been used to compute strains and stresses at the element Gauss points.For any load step, an iterative procedure has been used as well.For the J th iteration, the resulting incremental nodal displacement vector dUj, which is obtained through global tangent stiffness , is: accumulate the displacements: where, ∆ is the residual force vector and is accumulated displacement.and are the total displacement at the end of and 1 iteration.Meanwhile, is the incremental nodal displacement during the iteration.The total strain vector after the iteration is: where, is the strain at the end of the 1 iteration The incremental stress vector during the j th iteration is evaluated as follows: where is the material matrix based on the state of strain/stress at the chosen gauss point at the end of 1 iteration.Compute total stress vector, , where is the total stresses at the end of 1 and incremental stresses during the iteration.While adopting the incremental iterative technique, specification on convergence criteria is necessary for the termination of the iterations at particular incremental load step.In the present study, the convergence criterion to terminate the iteration process in the load step is based on the incremental displacements and it is expressed as follows:

Development of Computer Programmes
The computational procedures described in the previous sections have been written in the FORTRAN computer language working, under power station environment.The original finite element programme by Noorzaei ( 2007) was significantly modified.Eighteen new subroutines integrated with the finite element programme, see Figure 9.The new version of the finite element programme, called the three-dimensional finite element analysis of the girder bridges (3DFEMGB).
The new 3DFEMGB programme has the following features: (a) Is able to perform a linear elastic analysis of RC members (b) Elasto-Plastic analysis of RC components, (c) Frontal solver with multi-element and multi degree of freedom features.(d) Pre-processing facilities (e) Ottosen for the parameter yield criterion has been used in the programme (f) The new mentioned interface elements include: Linkage element; Truss-Linkage element; Interface element between plate bending and brick elements; Interface element between two plate bending elements; Interface element between two 20 noded Lagrangian brick elements Meanwhile, the inclusion of the above elements involves the implementation of the computation algorithm of the programme.An iterative tangent-stiffness procedure under incremental load was developed to perform a non-linear analysis for reinforced concrete structures.A flowchart of the procedure is shown in Figure 10.In the finite element analysis of structures, if the number of the element increases, it will be difficult to generate the element geometry input data without using any geometry generator tool, and for this reason, a mesh generation programme is written in the FORTRAN language.The programme has the ability to generate 3D solid elements and 2D elements.The outputs of the programme include:

Non-Linear Analysis of Reinforced Concrete Beam
Based on the non-linear modelling using Ottosen four-parameter model formulation presented in Sub-section 2.4, the computational procedures were carried out.Meanwhile, the output of the developed computational programme has been verified with the experimental work by Karihaloo (1990) and modelled by Ashour (1993).The beams were in the size of 100 mm in width by 150 mm in depth and 1600 mm in span.Beam 1 had one 12 mm reinforcing bar and beam 2 had two 12 mm reinforcement bars.The modulus of elasticity for steel reinforcement and concrete were Es= 210000 N/mm 2 and Ec = 30800 N/mm 2 , respectively.In the present study, the beam was modelled by nine elements in the x direction and three elements in the y direction (height of beam), as shown in Figure 11

Results
The load-deflection behaviour of the developed model, together with the experimental results by Karihaloo (1990) and model by Ashour (1993), is shown in Figure 12.Two beams [namely, Beam 1 that is reinforced with one 12 mm diameter bar as in Figure 12 A comparison of the results of the experimental and modelled reinforced concrete beams with non-linear models based on Ottosen, Mohr-Coulomb and Drucker-Prager yield criteria showed a reasonable agreement.The Ottosen model is the closest to that of Ashour's (1993), but the experimental result is almost similar to those of the Mohr-Coulomb yield criteria.

Non-Linear Analysis of the FRP Strengthened Reinforced Concrete Beams
The developed computational programme is further verified in this section based on the experimental results of the reinforced concrete beams strengthened with FRP plates by Rahimi and Hutchinson (2001), and the numerical studies by Carlos et al. (2006).All three yield criteria, Ottosen, Mohr and Drucker-prager, are also verified.The interface of plate bending and brick element as discussed in Section 2.1.5 is used in this model.Rahimi and Hutchinson (2001) tested three types of beam notated as A, B and C, in which the variable parameters are internal and external reinforcements.Only beams from series B are referred to validate the interface element between plate bending and brick element.There are eight type B beams with different external reinforcement and plate thickness apart from the two control beams.Details of the internal reinforcements are shown in Figure 13(a).The properties of the materials are described in Table 1, while the Type of external reinforcement and plate thickness are presented in Table 2.Additional details can be obtained from Rahimi and Hutchinson (2001).These beams are modelled using the three-dimensional finite element with internal reinforcement bars positioned at 140 mm apart in the width direction and 90 mm apart in the depth direc- The experimental results of type B beams by Rahimi and Hutchinson (2001), was used to verify the developed model.The developed model computed load-deflection curve of type B beams using three different yield criteria notated as Ottosen, Mohr and Drunker-Prager in Figures 14(a) to 14(d).The developed models of the reinforced concrete beams showed a reasonable agreement with those of the experimental results by Rahimi and Hutchinson (2001) and also the numerical studies by Carlos et al. (2006).The results for all the numerical studies are close enough for all the B1-B8 types, and these prove the ability of the model to capture the behaviour in linear and non-linear stages.In this study again, model that used Ottosen criteria, shows more deflection as compared to the other models, and this proves that the Ottosen model is a bit flexible compare to, Mohr and Drunker-Prager for concrete modelling.As shown in these figures, there are distributions of the slip between the concrete and FRP plates to justify the shear stress occurred between the two.In the region of maximum shear, the beams have maximum slip and in the region of zero shear stress, zero slip is recorded.Almost the same distribution of slip is computed in elastic stage, regardless of the thickness and properties of the FRP material.This is because slip is related to the amount of shear stress in the region between the plate and concrete, but it is not related to the type of plates.

Application to the Prefabricated Concrete Bridge Girder
Based on validation of the linear and non-linear finite element code for the reinforced concrete beam and the finite element interface developed code, the application of the developed model on the bridge girder is presented below.
Application of the developed code was made on a prefabricated highway bridge girder tested by Tianyu Xiang et al. (2005).The girder is a continuous reinforced concrete beam with four spans at 160 cm centre to centre.However, in the analysis of a prefabricated component, the girder was analysed as simply supported.The arrangement of the reinforcement of each prefabricated girder is given in Table 3, while the basic material properties are presented in Table 4 and the geometrical configuration of the bridge and dimension of the girder is shown in Figure 16.This bridge was designed according to the Specifications for Design of Reinforced Concrete and Pre-stressed Concrete Highway Bridges and Culverts (JTJ023-85) of PR China.The total dead load is 17KN/m for each girder, and the design vehicle load is shown in Figure 17.By using the G-M method (Johnson, 1980), the vehicle distribution loads of all beams were calculated.For simplicity and safety, the design vehicle load for all the beams was adopted for the maximum transverse distribution load, and the corresponding vehicle load distribution coefficient was set at 0.313.The service limit state and ultimate limit state were also analyzed.For the service limit state, the load combination is given as D+V, where D and V are dead load and vehicle load, respectively.The load combination for the ultimate limit state is expressed as 1.2D+1.1V.The prefabricated girder specimen was modelled with the three-dimensional finite element as shown in Figure 18.The total numbers of different elements used in the physical model is shown in the above-mentioned figure.The beam size is 1600 mm flange width by 1300 mm total depth.Meanwhile, the thickness of the web is 180 mm and the thickness of flange changes from 80 mm in the free end to 140 mm in the junction of the web and flange, with a clear distance between the supports of 19500 mm.

Results
In the ultimate limit state analysis, the variations of the vertical deflection, the relation between the maximum compressive stress of concrete and the maximum tensile stress of reinforcement bars, as well as the slip between bar and concrete in the fifth reinforcement layer and the moment at midspan are plotted in Figures 19 to 22    As shown in Figures 19 and 20, the variations of the vertical deflection, as well as the maximum tensile stress of reinforcement bars, are almost similar to the results for all the models, and this is sufficient to verify the ability of the proposed Ottosen four-parameter model.Figure 21 shows the results of the relationship between the maximum compressive stress of the concrete using three yield criteria have only a slight difference with that of Tianyu's work (2005), as they are in a close range.Meanwhile, the slip results from Figure 22 show the ability of the Truss-link element to model the interfacial behaviour between the bars and the concrete.As expected, the model, using the Ottosen yield criterion, showed higher value of deflections, and therefore, using this model
plate bending element with 5 DOF (3 translation and 2 rotation) is used to model the steel or FRP plates.The nodal displacements at any node i, are:

Figure 6 :
Figure 6: Interface Element between the Brick and Plate Bending Elements.

If
b and b , being the widths of the concrete prism and the FRP plate respectively, G f is the interfacial fracture energy and f is the concrete tensile strength.

K
first invariant of the stress tensor (10) J 2 = the second invariant of the stress deviator tensor and K , a and b are material parameters 0 K 1

Figure 9 :
Figure 9: Flowchart of the Modified Finite Element Programme.

Figure 10 :
Figure 10: Flowchart of the Non-linear Programme.
Property: Element Type, Nodal coordinates, Boundary conditions: GDATA Sets to zero arrays required for accumulation of data: ZERO Check Errors in the programmes: DOCTOR Read: Load Data Generate Statically Equivalent Nodal Forces: LDATA Generate elasticity matrix [D], B and Stiffness matrix [K]: STIFM Loop on no. of load increment Call shape function routine: SFR Stress level.IN-Determine the flow vector, YlELD and below.A new truss-link element was used to model the bond-slip between the concrete and the steel.In the developed code, apart from the Ottosen analytical yield criterion, two other yield criteria by Mohr-Coulomb and Drucker-Prager have been used for both beams.
(a), and Beam 2 with two 12 mm reinforcement bars in Figure 12(b) are presented.

Figure 14 :
Figure 14: Load vs Mid-Span deflection of Rahimi, Carlos and the present study using different Yield Criteria.

Figure 15 :
Figure 15: Slip (mm) between the FRP Plate and Concrete in the Adhesive Region of the Strengthened Beams.

Figure 16 :
Figure 16: Geometrical Configuration of the Bridge (all dimensions are in cm).Tianyu Xiang et al. (2005).

Figure 17 :
Figure 17: The longitudinal and Transverse Assignment of Vehicle Load (all dimensions are in cm).

Figure 21 :
Figure 21: Response of the Maximum Concrete Compressive Stress to Moment at Mid-Span.

Figure 22 :
Figure 22: Slip between the Maximum Stressed Reinforcement Bar (Fifth layer) and Concrete.

Table 1 :
Properties of the materials used.

Table 2 :
The Geometric properties of external reinforcement in beams type B.