Nonlinear analysis of reinforced concrete structures using thin flat shell elements

: This paper presents the development of a nonlinear finite element analysis program for reinforced concrete structures, subject to monotonic loading, using thin flat shell finite elements. The element thickness is discretized in concrete and steel layers. It is adopted the Newton-Raphson method, considering a secant stiffness approach for the Material Nonlinear Analysis, based on the Modified Compression Field Model (MCFT), unlike the usual tangent stiffness approach. The original formulation was expanded to also consider the Geometric Nonlinear Analysis, through a Total Lagrangian Formulation. The program was validated through comparison with experimental results, for different structures. It was observed good agreement, besides adequate computational cost.


INTRODUCTION
In reinforced concrete structures design, the civil engineer analyzes the real situation based on simplifying hypotheses, so that the structural models used in the analysis are sufficiently accurate and safe, but still having adequate simplicity, for use on project office.
Besides that, the development of construction technology has allowed the achievement of complex structures, with large spans and highly slender elements. In these cases, some of usual simplifying hypotheses may no longer represent where the subscript indicates the iteration number where the parameter must be evaluated, { } is the global nodal displacements vector, { } is the external forces vector, and [ ] and { } are, respectively, the tangent stiffness matrix and the internal forces vector that can be obtained through the corresponding elements contributions, [ ] and { }. The iterative process continues until a stopping criterion is met, such as: the iteration number exceeds the maximum value or, for a given tolerance , the relative norm of the difference between the vectors { } +1 and { } , convergence criterion parameter , is small enough, namely Equation 2.
According to De Borst et al. [8], it is important to apply the external forces incrementally, otherwise, due to the material nonlinear behavior or numerical characteristics of the solution procedure, in very large load steps, it is possible to arise serious convergence problems or inappropriate results. Thus, the solution procedure adopted in this paper is called incremental-iterative, using a load control approach, where Equation 1 is applied iteratively in each incremental load step. Also, according to De Borst et al. [8], since the solution procedure tends to reach an equilibrium configuration, in most cases, which stiffness matrix was adopted in the iterative process is less relevant. Based on this and knowing the numerical stability, which is often observed in secant stiffness analysis, even though the convergence rate may be lower compared to tangent stiffness analysis [9], in this paper, it was used a secant stiffness approach, unlike Barrales [5] who adopted a tangent stiffness approach.

Quadrilateral thin flat layered shell element -QTFLS
In this paper, the Quadrilateral Thin Flat Layered Shell Element -QTFLS, proposed by Barrales [5] and Rojas et al. [10], was adopted. It is a combination of the Quadrilateral Layered Membrane Element with Drilling Degrees of Freedom (DOF) -QLMD [11], and the Discrete Kirchhoff Quadrilateral Element -DKQ [12], where the Kirchhoff's assumptions for thin plates are considered. Figure 1 represents this 4-nodes finite element, the discretization of the element thickness in layers and its 6 degrees of freedom per node: 2 in-plane translations, 1 in-plane rotation (drilling), 2 out-of-plane rotations and 1 translation perpendicular to the element plane.  [10] In the QTFLS element, both the displacement and deformation fields are established through the superposition of membrane and plate behaviors. According to Barrales [5], this approach has the advantage of allowing different shape functions for each behavior.
As and � L � are well-known in the technical community and its detailed development can be easily found in appropriate bibliographies [5], [10]- [13]. However, to contribute to the paper completeness, these parameters will be briefly discussed in the following subsections.

Quadrilateral Layered Membrane Element with Drilling DOF -QLMD
The QLMD membrane element uses a combination of linear shape functions, Equation 3, and cubic Hermite functions, Equation 4. This 4-node element has 3 degree of freedom per node (2 in-plane translations and 1 in-plane rotation, drilling).
1 ( ) = The membrane displacement field in natural coordinates ( , ) is given by the following interpolation: where [MN] is a matrix defined by the shape functions in Equations 3 and 4 and [Tr] is a transformation matrix to ensure the compatibility between the rotation DOF, where 1 to 4 and 1 to 4 are the element nodes local coordinates. Thereby, the kinematic matrix [ ] can be obtained according to Equation 6, where [J] −1 represents the inverse of the Jacobian matrix [J] (2×2) , which relates, through bilinear shape functions, the natural ( , ) and local ( , ) coordinate systems. The corresponding formulation can be found in finite elements introductory textbooks [14].

Discrete Kirchhoff Quadrilateral Element -DKQ
According to Barrales [5], in the DKQ plate element, proposed by Batoz and Tahar [12], initially, the deflection and rotation fields are established independently, and later, they are related by applying the Kirchhoff's assumptions in a discrete manner on the element edges. For this purpose, it is adopted 8-node serendipity isoparametric element shape functions for the rotation fields, Equation 7, and cubic function for the deflections along the edges.
Although this shape functions are related to an 8-node element, in the DKQ element development, using: coordinate transformations, applying Kirchhoff's assumptions in a discrete manner on the element nodes, especially in nodes 5 to 8, and other simplifications, it is possible to reduce the element node number to 4, even using the Equation 7 shape functions. Consequently, the DKQ is a 4-node element that has 3 degree of freedom per node (2 out-of-plane rotations and 1 translation perpendicular to the element plane).
where [Ψ] is a matrix developed based on the discrete application of Kirchhoff's assumptions, and the geometric coefficients , , , and are functions of the element nodes local coordinates, Table 1. Thereby, the kinematic matrix [ 0 ] can be obtained according to Equation 9:

Geometric Nonlinearity
In the present paper, the geometric nonlinearity is considered through a Total Lagrangian Formulation (TLF), where the problem is analyzed in terms of the structure undeformed configuration. The Von Karman's hypotheses for large deflections of plates are considered, as presented by Figueiras [7]: • The shell thickness is small compared to the other dimensions; • The shell transverse deflection is of the same order of magnitude as the thickness; • The slopes are small, | / | ≪ 1 and | / | ≪ 1; • The tangential displacements and are small enough to allow the nonlinear terms associated with these fields be disregarded; • All the strain vector components are small.
Based on these hypotheses, the nonlinear plate strain component can be written as: where the vector { } with the derivatives of shell transverse deflection can be evaluated based on the element plate displacements { } and considering a bilinear interpolation in this field, as used in the Jacobian matrix [ ], Batoz and Tahar [12] and Rojas et al. [10], Equation 11: Based on the vector { } components, it is possible to assemble the matrix [ ]. Thereby, it is possible to note that these two elements depend on the nodal displacements { }, and the multiplication between them results in a nonlinear relationship between { } and { }. Calculating the variation of { } with respect to { } we obtain the matrix [ ]:

Element stiffness matrix
Through the Total Lagrangian Formulation, the element stiffness matrix [ ] is defined by two components [ ] and [ ] which consider, respectively, the large displacements and the stresses acting on the structure: ) where [ ] is the material tangent stiffness matrix, being adopted in its place, in this paper, the material secant stiffness matrix, and [ ] is a matrix defined according to the acting stresses. The element volume integrals are evaluated numerically. In the element plane, it is adopted the Gauss quadrature [14]. As the shell thickness is discretized in concrete layers and steel layers, in the integral along the thickness, it is used a mixed approach between the presented by Zhang et al. [15], [16], Barrales [5] and Vasilescu [13], which considers the sum of each layer individual contribution.
Thereby, using Equations 13 and 14, the matrix [ ] is evaluated as: The parameters and |[ ] | represent, respectively, the integration weights of the Gauss points, and the determinant of the corresponding Jacobian matrix [ ] .
The concrete layers contributions to the [ ] matrix, evaluated for each Gauss point and each concrete layer, are given by: The parameter [ ] , is the material secant stiffness matrix, evaluated at the concrete layer in the Gauss points, discussed in section 3. The terms 1 , 2 and 3 arise from the numerical integration of the layer coordinate parameter in the matrices [ 0 ], component of [ ] matrix, and are calculated in a discrete way, for each concrete layer, as shown in Equation 17, where and represent, respectively, the analyzed layer top and bottom coordinates.
The steel layers contributions to the [ ] matrix, evaluated for each Gauss point and each steel layer, are given by: The parameter [ ] , is the material secant stiffness matrix, evaluated at the steel layer in the Gauss points, also discussed in section 3. The variable refers to the position of the steel layer axis. On the other hand, and are the corresponding layer thickness and reinforcement ratio. Similarly, the matrix [ ] can be defined numerically as: where , and represent the stresses acting on the concrete and steel layers. Note that the matrix [ ] contributes to [ ] only in the region associated with the plate degrees of freedom, a consequence of considering nonlinear deformations only in the nonlinear plate strain component { }.

Internal forces vector
The internal forces vector { } unlike the external forces vector { } is not constant. It depends on the structure nodal displacement vector { } and must be updated at each iteration. The corresponding element contribution { } can be defined by Equation 20a and calculated numerically by Equation 20b.
where the matrices [ ] , and [ ] , can be calculated as:

Program implementation
The program was developed using a simple imperative concept with repetition statements, to perform the necessary calculations on all elements and Gaussian points, as well as for each load step. It is possible to apply the presented formulation using other strategies, such as the Object-Oriented Programming (OOP) to ensure greater code reusability.

MATERIAL CONSTITUTIVE MODELS
This section presents, in a concise way, the material constitutive models implemented in the developed computational program, since they are known formulations and widely discussed in the technical literature [2] [3], [17]. In addition, this section also details the formulation of the concrete [ ] and steel [ ] layers secant stiffness matrices, according to Vecchio [3]. In this study, the cracked concrete, despite its evident discrete nature, is modeled as a homogeneous orthotropic material, through the concept of mean stress and strain evaluated in regions containing several cracks, following the basis of the Modified Compression Field Theory (MCFT) [2].

Concrete in compression
Concrete in compression is modeled using a combination between the well-known Hognestad parabola, Equation 23, for both pre-peak and post-peak behavior, and the Vecchio 1992-A model (e1/e2-Form) [17], which through the softening coefficient , models the material strength loss due to transversal tension, Equation 24.
where and are, respectively, the average compressive stress and strain. The variable is the cylinder compressive strength (at 28 days), and 0 is the corresponding peak strain. The factor represents the influence of the relationship between the principal tensile 1 and compression 2 strains in concrete. The factor is equal to 0.55 when the slip deformations in the cracks are considered in the model (subsection 3.6), and 1.0 otherwise.

Concrete in tension
The concrete tensile model is defined by two distinct behaviors: pre-cracking and post-cracking, Equation 25. According to Wong et al. [17], after cracking, the reinforced concrete leaves the linear-elastic behavior, and the concrete tensile stresses tend to zero, on the crack surface, while it can present considerable values, between cracks, due to steel interaction. The Modified Bentz 2003 tensile stiffening model can represent this behavior. Furthermore, according to Vecchio [3], the magnitude of the average principal tensile stress must be limited by the remaining steel resistant capacity .
where is the concrete average principal tensile strain. The parameters and are, respectively, the crack stress and strain. The coefficient refers to the influence of the reinforcement on the stiffening, and it is obtained based on: the steel ratios and ; the reinforcement rebar nominal diameter and ; and the angles , and that illustrate the principal system direction and the reinforcements orientation. The parameters and are the reinforcements average stresses, while and represent these same parameters evaluated at the crack surface. The Modulus of elasticity of concrete illustrates the ratio between and , and can be estimated as 2 /| 0 |. When the reinforcement properties are directly assigned to the finite element, like in reinforced concrete membranes analysis programs [18], for example, the application of equation Equation 25 is direct. However, in the present study, the thickness of the shell element is discretized in layers and the constitutive model must be applied separately in each one. In this case, it is necessary to define which reinforcement parameters should be considered in each concrete layer. Hrynyk [4] presented a solution to this problem, based on CEB-FIP [19], which consists in defining a rebar tensile stiffening influence area equal to 7.5 times the rebar diameter. Thus, the concrete layers located within the rebar influence area are considered stiffened by this reinforcement. However, it is important to note that this evaluation must be done for all concrete layer-steel layer combinations, where a concrete layer could be stiffened by more than one steel layer (or none). Therefore, it is possible to see that the adoption of this criterion in the tensile stiffening model, along the shell thickness, tends to obtain a more realistic structural response. In the concrete layer plane state analysis, the principal stresses σ c1 and σ c2 can be evaluated with either the tensile model, Equation 25, or the compression model, Equation 23, depending on the layer plane state: biaxial tension, biaxial compression or tension-compression state.

Steel in tension or compression
In this paper, two steel constitutive models were implemented, both for tension and compression: a simple perfect elastic-plastic curve, Equation 26, and a bilinear curve that considers the material hardening after it reaches the yield condition, Equation 27.
where and are the yield stress and strain and and are the reinforcement average stress and strain. The parameters and are the corresponding steel ultimate stress and strain. The second model was implemented to allow the program to capture the post-yielding behavior, in the Polak shells problem, subsection 4.2.

Concrete confinement
Unlike the concrete softening effect, section 3.1, when this material is submitted to a biaxial (or triaxial) compression state, there is a confinement effect that tends to increase its resistant capacity. In the present study, this was modeled in a similar way to what was presented by Silva [18], based on the work of Vecchio [20], Kupfer et al. [21], Richart et al. [22] and Wong et al. [17], using the enhancement factors 1 and 2 : where and are the reinforcement ratio and its stress in the out-of-plane direction. The stress can be evaluated by applying the concrete strain in the corresponding direction in the steel constitutive model. In this paper, this strain is calculated in a simplified way regardless of whether the reinforcement yields or not as: where and are the concrete and steel modulus of elasticity in the out-of-plane direction, where is considered equal to |2 / 0 |. The parameter 12 represents the Poisson ratio that relates the strain in the 1-direction due to the stress in 2-direction [20]. The parameter 21 is defined similarly. The model adopted for calculating the Poisson coefficients 12 and 21 is detailed in the following subsection. The factors 1 and 2 are used to determine the peak stress and strain in the principal system (1-2): The peak stresses and strains can be used to calculate the concrete average principal compressive stresses. Analyzing the set of equations described in this subsection, it is possible to see the evaluation of the concrete behavior in biaxial compression as a nonlinear system of equations with two equations and two variables: This solution strategy was implemented in the developed program, where optimization functions from SciPy [23] scientific computational library were applied. It was observed good results, in addition to an adequate convergence, even applying a simple initial estimate for the solution (coordinate system origin).

Poisson ratio and lateral strains
According to Vecchio [20], lateral strains related to the Poisson ratio can be relevant for the reinforced concrete structures behavior, especially near failure. The Equation 33 represents the Poisson ratio model adopted in this paper. This model, in addition to considering the initial Poisson ratio 0 increase, also disregards this parameter when the concrete presents, in the transverse direction, a principal tensile strain greater than the cracking strain .
The concrete principal lateral strain in 1-direction 01 is given by Equation 34. It is important to note that the corresponding strain in 2-direction 02 and its Poisson ratio 21 can be obtained through Equations 33 and (34) switching the indexes.

Slip strain model
The Disturbed Stress Field Model (DSFM) is an extension of the well-known Modified Compression Field Theory (MCFT) [2], which admits disagreements between the stress and strain principal systems, through the consideration of crack slip strains. The DSFM was proposed by Vecchio [3] to solve some MCFT drawbacks in calculate the strength and stiffness of high or low reinforcement ratio elements.
The crack surface shear stress can be evaluated applying equilibrium in the reinforced concrete element: According to Vecchio [3] due to this shear stress there is a rigid body local slip along the crack (slip displacement ) which causes slip strains { }. These strains must be considered in the model additionally to the principal strains related to the material constitutive response. According to Vecchio [3], the displacement can be calculated as: where represents the average crack width, which can be estimated from the average crack spacing , Equation 39, based on the nominal crack spacings in x and y ( = = 50 ).
When the average crack width is greater than or equal to 5 mm, the concrete principal compressive stress 2 is considered equal to zero. The parameter is cubic compressive strength, adopted as: = /0.85. The crack slip shear strain is evaluated as the ratio between and . Based on Mohr's circle coordinate transformations, the slip strain vector in the Cartesian system { } can be written as: The formulation described above represents an overview of the slip strain vector { } calculation procedure. However, it is important to present some additional details about the shear stresses along the crack surfaces . Although this parameter can be obtained by Equation 37, according to Silva [18], it should be limited to the maximum shear stress that can be resisted on the crack by aggregate interlock, Equation 41.
where is the aggregate size (adopted as 25mm). Furthermore, according to Equation 37, to evaluate the shear stress it is previously necessary to calculate the reinforcement stresses in the crack surface and . Although these parameters can be easily obtained using the steel constitutive models, section 3.3, the evaluation of the corresponding steel strains, in the crack surface, and , necessary to obtain these stresses, is not immediate. Through the reinforced concrete element equilibrium in the crack surface, there is a steel stress (and strain) increment, due to the concrete tensile stress absence. Thus, Vecchio [3] proposes that the steel strains, in the cracks and , should be determined through the sum of the corresponding average strain and and the local incremental strain contribution in 1-direction, Δ 1 : Thus, considering the concrete in tension and steel constitutive models presented and the formulation described above, it is possible to formulate a nonlinear equation whose solution is the incremental strain Δ 1 , Equation 43. Again, it was used SciPy [23] scientific computing library optimization functions.
Finally, it is important to note that the slip theory is only applied to cracked concrete. Thus, it must be verified whether the concrete principal tensile strain 1 is greater than the crack strain . Otherwise, the slip strain vector { } can be computed as a null vector.

Material secant stiffness matrices
In the previous subsections, it was presented the lateral strain vector { 0 } and the slip strain vector { } calculation procedure. In order to ensure the concrete secant stiffness matrix symmetry and the associated benefits, according to Vecchio [3] and Silva [18], the concrete total strain vector { } is defined by three distinct components: { }, { 0 } and { }, where, { } represents the elastic strains due to stress, which can be evaluated, as: Once the vector { } are obtained, it is possible to estimate the inclination of the principal strains , Equation 45, and the concrete principal strains: Through Equations 44 and (46) and the models presented in the previous subsections, it is possible to see the nonlinear relationship between the vectors { }, { 0 } and { }. In the present paper, this problem was solved iteratively, as illustrated in Figure 3. In the cracked reinforced concrete element, the concrete stress vector { } can be related to the corresponding strain vector { } as: where the concrete layer secant stiffness matrix [ ] in the Cartesian system is obtained by a coordinate transformation of the corresponding matrix in the principal system [ 1−2 ], which is defined based on the concrete secant moduli � 1 and � 2 , Equations 48 and 49.

PROGRAM VALIDATION AND DISCUSSIONS
This section presents the developed computer program validation through comparison with experimental and some numerical results [1], [4], available in the technical literature, for different structures. The load-displacement curves presented were created by the program, while the nodal displacements diagrams and the average internal forces diagrams, in the Gauss points, were obtained using the software Paraview [24], a Python module called PyEVTK [25] and additional codes written by the author, in the same language. Other useful structure diagrams for practical applications, like principal stresses, reinforcement stresses and crack pattern, are features that have not been implemented yet in the presented code. However, it can be done using the same approach mentioned above. In fact, any node or element property, in each load step, can be represented this way, in the program post-processor.

Cervenka deep beam
Initially, to evaluate the program performance in material nonlinear membrane problems, the deep beam W2 tested by Cervenka [26], Figure 5, was analyzed. All the plate degrees of freedom have been fixed. It was adopted 3 steel layers (2 horizontal and 1 vertical) to model the reinforcement, where its positions were defined according to the experiment. It was considered the steel perfect elasto-plastic model. The external load was applied using an initial load of 40 kN, in addition to 85 increments equal to 0.88 kN. The stopping criteria tolerance , associated with the increment displacement criterion described in subsection 2.1, was equal to 1%, . Figure 5 illustrates the results obtained, using the problem symmetry. Figure 4 also shows a comparison between the obtained load-displacement curve (ydisplacement at the deep beam bottom midpoint), the experiment and a literature numerical response [27]. It was observed an adequate structural behavior and a good accuracy to the experimental data. The processing time was about 6 minutes, using a processor: Intel Core i7-5500U CPU @ 2.40GHz.

Polak shells
To evaluate the program performance in material nonlinear plate and shell problems, three reinforced concrete shells (SM1, SM2 and SM3) tested by Polak and Vecchio [28] were analyzed. The specimens' characteristics, in addition to the results obtained, in comparison with literature numerical solutions [1], [4], are illustrated in Figure 6. It was adopted 4 steel layers (2 horizontal and 2 vertical), in addition to an out-of-plane reinforcement, to model the structures rebars, where its positions were defined according to the experiment. In all three cases, the load increments number was about 90. The stopping criteria tolerance was equal to 1%, subsection 2.1. The finite element mesh used contains 8x8 elements, and its thickness were discretized into 10 concrete layers. In the problem analysis, the program was not able to find a post-yielding equilibrated response, in less than 100 iterations (default maximum iterations number). However, during the study, it was observed that, considering the bilinear steel constitutive model, Equation 27, and disregarding the confinement (subsection 3.4) and slip strain models (subsection 3.6), consequently a simpler set of constitutive models, the tool found an equilibrated configuration, between 20 and 50 iterations. After these considerations, again, a good agreement was observed between the developed program and the experimental results. The total processing time was 15 minutes, where most of this time refer to the post-yielding behavior.

Geometric nonlinear plates analysis
Three rectangular linear-elastic plates, subjected to a uniform load , presented by Figueiras [7], were analyzed to verify the geometric nonlinear model implemented. The plates span were equal to 6m and its thickness ℎ was assumed as 0.15m. The material elastic modulus and Poisson ratio were adopted, respectively, as 30 / 2 and 0.316. The difference between the three structures was the boundary conditions applied to the edges: clamped, simply supported (horizontally fixed) and simply supported (horizontally free). The finite element mesh used contains 6x6 elements, and its thickness were discretized into 10 concrete layers. In all three cases, the number of load increments was equal to 100. The stopping criteria tolerance was equal to 10 −5 . Figure 7 shows the load parameter ( . 4 / . ℎ 4 ) versus central displacement parameter ( /ℎ) curves obtained in this study, where it is possible to observe good accuracy when compared with literature analytical [29] and numerically [7] solutions. The processing time was about 12 minutes.

Slender shear wall
The last problem analyzed was a slender shear wall. The structure geometry was adapted from the literature [6], [30], [31] to reach a wall slenderness ratio equal to 90. It was considered both the material and the geometric nonlinearities The problem characteristics are illustrated in Figure 8. The thickness was discretized into 10 concrete layers. The shear wall is simply supported. The external loads were applied in 10 increments and stopping criteria tolerance was equal to 1%. The results obtained in the developed program were compared with VecTor 4 structural analysis software in Figure 8c. The two tools obtained similar results. However, it is important to emphasize that a proper validation must be performed using experimental results. VecTor 4 was adopted as a reference given the lack slender shear walls test results, like Figure 7a. The total processing time was close to 10 minutes.

CONCLUSIONS
This paper presents the development of a nonlinear finite element analysis program for reinforced concrete structures, subject to monotonic loading, using thin flat shell finite elements QTFLS [5]. The material nonlinear analysis considered a secant stiffness approach, based on the Modified Compression Field Theory (MCFT) [2]. The element original formulation was expanded to also consider the problem geometric, through a Total Lagrangian Formulation [7]. Based on this study, the following was observed: • Reinforced concrete structures nonlinear analysis, based on the Newton-Raphson method, using the materials secant stiffness matrices and a Total Lagrangian Formulation can be considered an attractive approach, according to the results accuracy and the computational cost; • However, the fact that part of the formulation needed to be adjusted to the program be able to find a post-yielding equilibrated response, in some problems (subsections 4.2), shows the difficulty present in the development of a tool with wide range of potential applications, as also exposed by Figueira [7]. • The slender shear wall validation, subsection 4.4, indicates the necessity to conduct experimental programs for this type of structure, in order not only to obtain a better understanding of the construction behavior, but also to produce test data to enable the development of more accuracy computational tools.