Nonlinear analysis of reinforced concrete structures using thin flat shell elements

Análise não linear de estruturas de concreto armado usando elementos finitos de cascas finas e planas

abstract:

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.

Keywords:
reinforced concrete; finite elements; nonlinear analysis; shells; plates

resumo:

O presente artigo apresenta o desenvolvimento de uma ferramenta para a análise não-linear de estruturas de concreto armado, sujeitas a carregamentos monotônicos, utilizando o elemento finito de cascas finas e planas, o qual, é discretizado, ao longo da sua espessura, em lamelas de concreto e camadas de aço. É utilizado o método de Newton-Raphson, adotando, na consideração da não-linearidade física, a rigidez secante do material, baseada no modelo do campo de compressão modificado, no lugar da abordagem via rigidez tangente. A formulação original do elemento foi expandida para considerar também a não-linearidade geométrica, através de uma formulação Lagrangiana total. A validação da ferramenta é via comparação com resultados experimentais da literatura, para diversas estruturas, onde pode ser observada boa aderência além de adequado custo computacional.

Palavras-chave:
concreto armado; elementos finitos; análise não-linear; cascas; placas

1 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 the actual structural behavior, due to the increase of relevant nonlinear effects associated with the material response, such as concrete cracking (material nonlinearity) or large displacements (geometric nonlinearity).

Fortunately, the modern computer has allowed the use of sophisticated structural models, once considered unfeasible for practical applications, to gain space in the market, including nonlinear finite element analysis (NL-FEA). It leads the technical community to a constant review of the structural models used, always seeking to associate a safe design with the productivity resulting from the most efficient technologies available.

Several structural practical applications can be analyzed through shell models, such as: shear walls, shell roofs, water tanks or other storage structures. The analysis complexity and the computer development has been stimulating the search for numerical tools to solve this problem (shell finite elements). The challenge becomes greater when the structure material is reinforced concrete (RC), where the material behavior plays a crucial role in the construction response. Consequently, the material constitutive models are a determining point for a satisfactory analysis. Among the most common approaches for shell element formulation, one can mention degenerate shell elements, which are based on three-dimensional equilibrium equations, and shell elements developed by the superposition of membrane and plate elements.

Following the first option, Luu et al. [¹] proposed the CSMM-based shell element for reinforced concrete structures, for material nonlinear analysis (MNA). It uses the smeared crack theory Cyclic Softened Membrane Model (CSMM), created at the University of Houston. This 8-node degenerate shell element has 40 degrees of freedom (DOF): 3 translations and 2 rotations per node, where each nodal rotation follows a specific nodal coordinate system. Its nonlinear analysis procedure uses a Newton-Raphson approach (tangent stiffness).

Another notorious tool that also makes use of degenerate shell elements is VecTor4. This software considers both MNA and Geometric Nonlinear Analysis (GNA), through a Total Lagrangian Formulation (TLF). It was developed at the University of Toronto and its material model is based on the Modified Compression Field Theory (MCFT) [²], [³]. The VecTor4 quadrilateral 9-node shell element has 42 DOF: 3 translation (all nodes) and two rotations (only at the edge nodes). Again, the nodal rotations follow specific nodal coordinate systems, defined in each node. It is an element that combines in its formulation Langrangian and Serendipity shape functions, being therefore also called heterosis [⁴]. Its nonlinear analysis procedure differs from Luu et al. [¹] and uses a direct secant stiffness approach.

Following the other mentioned option in shell element development (superposing membrane and plate elements), Barrales [⁵] proposed a simple and efficient quadrilateral thin fat layered shell element (QTFLS), for MNA. This 4-node element has 6 DOF per node (3 translations and 3 rotations). Its nonlinear analysis procedure also uses a Newton-Raphson approach (tangent stiffness). An interesting point about this formulation is that, unlike degenerate shell elements, this element explicitly has the nodal in-plane rotation DOF (drilling). Therefore, it is not necessary to use nodal coordinate systems to assemble the elements rotation DOF. According to Silva and Horowitz [⁶], when modeling U-Shaped RC shear walls using degenerate elements, such as VecTor4, special attention is required in the rotations DOF compatibility, in the L-connection elements (through nodal coordinate systems). The Barrales [⁵] element requires only local and global coordinate systems, an attractive feature. As usual in RC Shell structures NL-FEA, all discussed elements have its thickness discretized in concrete and steel layers, to properly consider the internal stresses variation along the thickness.

This paper expands the Barrales [⁵] QTFLS shell element formulation to consider both MNA and GNA, using a TLF [⁷]. The formulation is implemented in a NL-FEA program for RC structures, subject to monotonic loads. The element thickness is also discretized in concrete and steel layers. It is adopted the Newton-Raphson method, but considering a secant stiffness approach, using the basis of the Modified Compression Field Model (MCFT), unlike the original tangent stiffness approach. The program was validated through comparison with experimental and numerical results [¹], [⁴], for different structures. It was observed good agreement, besides adequate computational cost.

2 NONLINEAR FINITE ELEMENT PROCEDURE

This section presents the material and geometric nonlinear finite element procedure used for analysis of reinforced concrete structures.

2.1 Incremental-iterative procedure

In nonlinear analysis, the set of equilibrium equations can be obtained through the principle of virtual works (PVW). Based on these equations, and discretizing the structure in finite elements, usually, the problem can be solved iteratively, using the Newton-Raphson method, Equation 1:

${\left\{d\right\}}_{n+1}={\left\{d\right\}}_n+{\left[K_G\right]}_n^{-1}\left(\left\{F_{ext}\right\}-\left\{F_{int}\right\}\right)$ (1)

where the subscript $n$ indicates the iteration number where the parameter must be evaluated, $\left\{d\right\}$ is the global nodal displacements vector, $\left\{F_{ext}\right\}$ is the external forces vector, and $\left[K_G\right]$ and $\left\{F_{int}\right\}$ are, respectively, the tangent stiffness matrix and the internal forces vector that can be obtained through the corresponding elements contributions, $\left[k_e\right]$ and $\left\{f_e\right\}$. The iterative process continues until a stopping criterion is met, such as: the iteration number exceeds the maximum value or, for a given tolerance $tol$, the relative norm of the difference between the vectors ${\left\{d\right\}}_{n+1}$ and ${\left\{d\right\}}_n$, convergence criterion parameter $error$, is small enough, namely Equation 2.

$error=\mathrm{}\frac{\left|{\left\{d\right\}}_{n+1}-{\left\{d\right\}}_n\right|}{\mathrm{}\left|{\left\{d\right\}}_{n+1}\right|}<tol$ (2)

According to De Borst et al. [⁸], 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. [⁸], 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 [⁹], in this paper, it was used a secant stiffness approach, unlike Barrales [⁵] who adopted a tangent stiffness approach.

2.2 Quadrilateral thin flat layered shell element - QTFLS

In this paper, the Quadrilateral Thin Flat Layered Shell Element - QTFLS, proposed by Barrales [⁵] and Rojas et al. [¹⁰], was adopted. It is a combination of the Quadrilateral Layered Membrane Element with Drilling Degrees of Freedom (DOF) - QLMD [¹¹], and the Discrete Kirchhoff Quadrilateral Element - DKQ [¹²], 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.

Figure 1
Quadrilateral Thin Flat Layered Shell Element - QTFLS [¹⁰]

In the QTFLS element, both the displacement and deformation fields are established through the superposition of membrane and plate behaviors. According to Barrales [⁵], this approach has the advantage of allowing different shape functions for each behavior.

As usual in finite element analysis, the strain vector $\left\{\mathrm{\epsilon }\right\}$ can be related to the element displacement vector $\left\{{\mathrm{d}}_{\mathrm{e}}\right\}$ through the kinematic matrix $\left[\mathrm{B}\mathrm{\text{'}}\right]$. Equation 2 expands this relationship by superposing the linear membrane component $\left\{{\mathrm{\epsilon }}_{\mathrm{m}}\right\}$, the linear plate component $\left\{{\mathrm{\epsilon }}_{\mathrm{b}}^0\right\}$ and the nonlinear plate component $\left\{{\mathrm{\epsilon }}_{\mathrm{b}}^{\mathrm{L}}\right\}$:

$\left\{\mathrm{\epsilon }\right\}=\left[B\mathrm{\text{'}}\right]\left\{d_e\right\}=\left\{{\epsilon }_m\right\}+\left\{{\epsilon }_b^0\right\}+\left\{{\epsilon }_b^L\right\}=\left[B_m\right]\left\{d_e^m\right\}+z_L\left[B_b^0\right]\left\{d_e^b\right\}+\frac{1}{2}\left[B_b^L\right]\left\{d_e^b\right\}$ (2)

where $\left[B_m\right]$ and $\left[B_b^0\right]$ are the kinematic matrices that represent the linear relationship between the membrane $\left\{d_e^m\right\}$ and plate $\left\{d_e^b\right\}$ element displacements and the corresponding strain component $\left\{{\epsilon }_m\right\}$ and $\left\{{\epsilon }_b^0\right\}$. On the other hand, the matrix $\left[B_b^L\right]$ is related to the consideration of the geometric nonlinearity of the problem, through the nonlinear plate strain component $\left\{{\epsilon }_b^L\right\}$, discussed in section 2.5. The $z_L$ parameter refers to the layer $z$ local coordinate. The formulation of the matrices $\left[B_m\right]$, $\left[B_b^0\right]$ and $\left[B_b^{\mathrm{L}}\right]$ are well-known in the technical community and its detailed development can be easily found in appropriate bibliographies [⁵], [¹⁰]–[¹³]. However, to contribute to the paper completeness, these parameters will be briefly discussed in the following subsections.

2.3 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).

$\begin{array}{cc}{M}_1\left(\eta \right)=\frac{1}{2}(1-\eta )\mathrm{}& \begin{array}{cc}\mathrm{}& {\mathrm{}M}_2\left(\eta \right)=\frac{1}{2}(1+\eta )\mathrm{}\end{array}\end{array}$ (3)

$\begin{array}{cc}{N}_1\left(\eta \right)=\frac{1}{2}-\frac{3}{4}\eta +\frac{1}{4}{\eta }^3& \begin{array}{cc}\mathrm{}& N_2\left(\eta \right)=+\frac{1}{4}-\frac{1}{4}\eta -\frac{1}{4}{\eta }^2+\frac{1}{4}{\eta }^3\end{array}\end{array}$ (4a)

$\begin{array}{cc}\mathrm{}& \mathrm{}\\ N_3\left(\eta \right)=\frac{1}{2}+\frac{3}{4}\eta -\frac{1}{4}{\eta }^3& \begin{array}{cc}\mathrm{}& N_4\left(\eta \right)=-\frac{1}{4}-\frac{1}{4}\eta +\frac{1}{4}{\eta }^2+\frac{1}{4}{\eta }^3\end{array}\end{array}$ (4b)

The membrane displacement field in natural coordinates $(\xi ,\eta )$ is given by the following interpolation:

$\left\{\begin{array}{c}{u}_m\\ {\nu }_m\end{array}\right\}={\left[MN\right]}_{(2\times 16)}{\left[Tr\right]}_{(16\times 12)}{\left\{d_e^m\right\}}_{(12\times 1)}$ (5a)

$\left[MN\right]={\left[\begin{array}{cc}\begin{array}{c}{M}_1\left(\xi \right)N_1\left(\eta \right)\\ \begin{array}{c}0\\ \begin{array}{c}{-M}_1\left(\xi \right)N_2\left(\eta \right)\\ \begin{array}{c}0\\ \begin{array}{c}{M}_2\left(\xi \right)N_1\left(\eta \right)\\ \begin{array}{c}0\\ \begin{array}{c}{-M}_2\left(\xi \right)N_2\left(\eta \right)\\ \begin{array}{c}0\\ M_2\left(\xi \right)N_3\left(\eta \right)\\ \begin{array}{c}0\\ {-M}_2\left(\xi \right)N_4\left(\eta \right)\\ \begin{array}{c}0\\ M_1\left(\xi \right)N_3\left(\eta \right)\\ \begin{array}{c}0\\ \begin{array}{c}{-M}_1\left(\xi \right)N_4\left(\eta \right)\\ 0\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}0\\ \begin{array}{c}{M}_1\left(\eta \right)N_1\left(\xi \right)\\ \begin{array}{c}0\\ \begin{array}{c}{M}_1\left(\eta \right)N_2\left(\xi \right)\\ \begin{array}{c}0\\ \begin{array}{c}{M}_1\left(\eta \right)N_3\left(\xi \right)\\ \begin{array}{c}0\\ \begin{array}{c}{M}_1\left(\eta \right)N_4\left(\xi \right)\\ 0\\ \begin{array}{c}{M}_2\left(\eta \right)N_3\left(\xi \right)\\ 0\\ \begin{array}{c}{M}_2\left(\eta \right)N_4\left(\xi \right)\\ 0\\ \begin{array}{c}{M}_2\left(\eta \right)N_1\left(\xi \right)\\ \begin{array}{c}0\\ M_2\left(\eta \right)N_2\left(\xi \right)\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]}^T$ (5b)

$\left[Tr\right]=\frac{1}{2}\left[\begin{array}{cccccccccccc}2& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 2& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& y_4-y_1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& x_2-x_1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 2& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 2& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& y_3-y_2& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& x_2-x_1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 2& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 2& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& y_3-y_2& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& x_3-x_4& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 2& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 2& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& y_4-y_1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& x_3-x_4\end{array}\right]$ (5c)

$\left\{d_e^m\right\}={\left\{\begin{array}{cccccccccccc}{u}_1& v_1& {\theta }_{z1}& u_2& v_2& {\theta }_{z2}& u_3& v_3& {\theta }_{z3}& u_4& v_4& {\theta }_{z4}\end{array}\right\}}^T$ (5d)

where $\left[\mathrm{M}\mathrm{N}\right]$ is a matrix defined by the shape functions in Equations 3 and 4 and $\left[\mathrm{T}\mathrm{r}\right]$ is a transformation matrix to ensure the compatibility between the rotation DOF, where $x_1$ to $x_4$ and $y_1$ to $y_4$ are the element nodes local coordinates. Thereby, the kinematic matrix $\left[B_m\right]$ can be obtained according to Equation 6, where ${\left[\mathrm{J}\right]}^{-1}$ represents the inverse of the Jacobian matrix ${\left[\mathrm{J}\right]}_{(2\times 2)}$, which relates, through bilinear shape functions, the natural $(\xi ,\eta )$ and local $(x,y)$ coordinate systems. The corresponding formulation can be found in finite elements introductory textbooks [¹⁴].

$\left\{{\epsilon }_m\right\}=\left\{\begin{array}{c}\frac{\partial u_m}{\partial x}\\ \frac{\partial {\nu }_m}{\partial y}\\ \frac{\partial u_m}{\partial y}+\frac{\partial {\nu }_m}{\partial x}\end{array}\right\}=\left[B_m\right]\left\{d_{em}\right\}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 1& 1& 0\end{array}\right]\left[\begin{array}{cc}{\left[J\right]}^{-1}& \left[0\right]\\ \left[0\right]& {\left[J\right]}^{-1}\end{array}\right]\left[\begin{array}{c}\frac{{\partial MN}_{1,i}}{\partial \xi }\\ \frac{{\partial MN}_{1,i}}{\partial \eta }\\ \frac{{\partial MN}_{2,i}}{\partial \xi }\\ \frac{{\partial MN}_{2,i}}{\partial \eta }\end{array}\right]\left[Tr\right]\left\{d_e^m\right\}\mathrm{}$ (6)

2.4 Discrete Kirchhoff Quadrilateral Element - DKQ

According to Barrales [⁵], in the DKQ plate element, proposed by Batoz and Tahar [¹²], 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.

${\psi }_i\left(\xi ,\eta \right)=-\frac{1}{4}\left(1+{\xi }_i\xi \right)\left(1+{\eta }_i\eta \right)\left(1-{\xi }_i\xi -{\eta }_i\eta \right)\mathrm{}for\mathrm{}i=1,\mathrm{}2,\mathrm{}3\mathrm{}and\mathrm{}4$ (7a)

${\psi }_k\left(\xi ,\eta \right)=\frac{1}{2}\left(1-{\xi }^2\right)\left(1+{\eta }_k\eta \right)\mathrm{}for\mathrm{}i=5\mathrm{}and\mathrm{}6$ (7b)

${\psi }_k\left(\xi ,\eta \right)=\frac{1}{2}\left(1+{\xi }_k\xi \right)\left(1-{\eta }^2\right)\mathrm{}for\mathrm{}i=7\mathrm{}and\mathrm{}8$ (7c)

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).

According to Rojas et al. [¹⁰], the middle surface normal rotation field of the plate element, in natural coordinates $(\xi ,\eta )$, is given by the following interpolation:

$\left\{\begin{array}{c}{\beta }_x\\ {\beta }_y\end{array}\right\}={\left[\mathrm{\Psi }\right]}_{(2\times 12)}{\left\{d_e^b\right\}}_{(12\times 1)}$ (8a)

$\left[\mathrm{\Psi }\right]={\left[\begin{array}{rr}\mathrm{}3/2\left[{\psi }_5\left(\xi ,\eta \right)a_5-{\psi }_8\left(\xi ,\eta \right)a_8\right]& 3/2\left[{\psi }_5\left(\xi ,\eta \right)d_5-{\psi }_8\left(\xi ,\eta \right)d_8\right]\\ \mathrm{}{\psi }_5\left(\xi ,\eta \right)b_5+{\psi }_8\left(\xi ,\eta \right)b_{8\mathrm{}}& -{\psi }_1\left(\xi ,\eta \right)+{\psi }_5\left(\xi ,\eta \right)e_5+{\psi }_8\left(\xi ,\eta \right)e_8\mathrm{}\\ {\psi }_1\left(\xi ,\eta \right)-{\psi }_5\left(\xi ,\eta \right)c_5-{\psi }_8\left(\xi ,\eta \right)c_8& \mathrm{}{-\psi }_5\left(\xi ,\eta \right)b_5-{\psi }_8\left(\xi ,\eta \right)b_{8\mathrm{}}\\ 3/2\left[{\psi }_6\left(\xi ,\eta \right)a_6-{\psi }_5\left(\xi ,\eta \right)a_5\right]& 3/2\left[{\psi }_6\left(\xi ,\eta \right)d_6-{\psi }_5\left(\xi ,\eta \right)d_5\right]\\ \mathrm{}{\psi }_6\left(\xi ,\eta \right)b_6+{\psi }_5\left(\xi ,\eta \right)b_5\mathrm{}& -{\psi }_2\left(\xi ,\eta \right)+{\psi }_6\left(\xi ,\eta \right)e_6+{\psi }_5\left(\xi ,\eta \right)e_{5\mathrm{}}\\ {\psi }_2\left(\xi ,\eta \right)-{\psi }_6\left(\xi ,\eta \right)c_6-{\psi }_5\left(\xi ,\eta \right)c_5& \mathrm{}{-\psi }_6\left(\xi ,\eta \right)b_6-{\psi }_5\left(\xi ,\eta \right)b_5\mathrm{}\\ \mathrm{}3/2\left[{\psi }_7\left(\xi ,\eta \right)a_7-{\psi }_6\left(\xi ,\eta \right)a_6\right]& 3/2\left[{\psi }_7\left(\xi ,\eta \right)d_7-{\psi }_6\left(\xi ,\eta \right)d_6\right]\\ \mathrm{}{\psi }_7\left(\xi ,\eta \right)b_7+{\psi }_6\left(\xi ,\eta \right)b_{6\mathrm{}}& -{\psi }_3\left(\xi ,\eta \right)+{\psi }_7\left(\xi ,\eta \right)e_7+{\psi }_6\left(\xi ,\eta \right)e_{6\mathrm{}}\\ {\psi }_3\left(\xi ,\eta \right)-{\psi }_7\left(\xi ,\eta \right)c_7-{\psi }_6\left(\xi ,\eta \right)c_6& \mathrm{}{-\psi }_7\left(\xi ,\eta \right)b_7-{\psi }_6\left(\xi ,\eta \right)b_{6\mathrm{}}\\ 3/2\left[{\psi }_8\left(\xi ,\eta \right)a_8-{\psi }_7\left(\xi ,\eta \right)a_7\right]& 3/2\left[{\psi }_8\left(\xi ,\eta \right)d_8-{\psi }_7\left(\xi ,\eta \right)d_7\right]\\ {\mathrm{}\psi }_8\left(\xi ,\eta \right)b_8+{\psi }_7\left(\xi ,\eta \right)b_{7\mathrm{}}& {-{\psi }_4\left(\xi ,\eta \right)+\psi }_8\left(\xi ,\eta \right)e_8+{\psi }_7\left(\xi ,\eta \right)e_{7\mathrm{}}\\ {\psi }_4\left(\xi ,\eta \right)-{\psi }_8\left(\xi ,\eta \right)c_8-{\psi }_7\left(\xi ,\eta \right)c_7& {\mathrm{}-\psi }_8\left(\xi ,\eta \right)b_8-{\psi }_7\left(\xi ,\eta \right)b_7\mathrm{}\end{array}\right]}^T$ (8b)

$\left\{d_e^b\right\}={\left\{\begin{array}{cccccccccccc}{w}_1& {\theta }_{x1}& {\theta }_{y1}& w_2& {\theta }_{x2}& {\theta }_{y2}& w_3& {\theta }_{x3}& {\theta }_{y3}& w_4& {\theta }_{x4}& {\theta }_{y4}\end{array}\right\}}^T$ (8c)

where $\left[\mathrm{\Psi }\right]$ is a matrix developed based on the discrete application of Kirchhoff's assumptions, and the geometric coefficients $a_k$, $b_k$, $c_k$, $d_k$ and $e_k$ are functions of the element nodes local coordinates, Table 1. Thereby, the kinematic matrix $\left[B_b^0\right]$ can be obtained according to Equation 9:

Table 1
Geometric coefficients.

Coefficient	Equation*
$a_k$	$-\frac{(x_i-x_j)}{{(x_i-x_j)}^2+{(y_i-y_j)}^2}$
$b_k$	$\frac{3(x_i-x_j)(y_i-y_j)}{{4(x_i-x_j)}^2+{4(y_i-y_j)}^2}$
$c_k$	$\frac{{(x_i-x_j)}^2/4-{(y_i-y_j)}^2/2}{{(x_i-x_j)}^2+{(y_i-y_j)}^2}$
$d_k$	$-\frac{(y_i-y_j)}{{(x_i-x_j)}^2+{(y_i-y_j)}^2}$
$e_k$	$\frac{{-(x_i-x_j)}^2/2+{(y_i-y_j)}^2/4}{{(x_i-x_j)}^2+{(y_i-y_j)}^2}$

• ^*
  The indexes (k,i,j) can be related as 5,1,2, 6,2,3, 7,3,4 and 8,4,1.

$\left\{{\epsilon }_b^0\right\}=z_L\left\{\begin{array}{c}\frac{\partial {\beta }_x}{\partial x}\\ \frac{\partial {\beta }_y}{\partial y}\\ \frac{\partial {\beta }_x}{\partial y}+\frac{\partial {\beta }_y}{\partial x}\end{array}\right\}=z_L\left[B_b^0\right]\left\{d_{eb}\right\}=z_L\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 1& 1& 0\end{array}\right]\left[\begin{array}{cc}{\left[J\right]}^{-1}& \left[0\right]\\ \left[0\right]& {\left[J\right]}^{-1}\end{array}\right]\left[\begin{array}{c}\frac{{\partial \mathrm{\Psi }}_{1,i}}{\partial \xi }\\ \frac{\partial {\mathrm{\Psi }}_{1,i}}{\partial \eta }\\ \frac{\partial {\mathrm{\Psi }}_{2,i}}{\partial \xi }\\ \frac{\partial {\mathrm{\Psi }}_{2,i}}{\partial \eta }\end{array}\right]\left\{d_e^b\right\}\mathrm{}$ (9)

2.5 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 [⁷]:

• The shell thickness $t$ is small compared to the other dimensions;

• The shell transverse deflection $w$ is of the same order of magnitude as the thickness;

• The slopes are small, $\left|\partial w/\partial x\right|\ll 1$ and $\left|\partial w/\partial y\right|\ll 1$;

• The tangential displacements $u$ and $v$ 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:

$\left\{{\epsilon }_b^L\right\}=\left\{\begin{array}{c}{\frac{1}{2}\left(\frac{\partial w}{\partial x}\right)}^2\\ {\frac{1}{2}\left(\frac{\partial w}{\partial y}\right)}^2\\ \frac{\partial w}{\partial x}\frac{\partial w}{\partial y}\end{array}\right\}=\frac{1}{2}\left[\begin{array}{cc}\frac{\partial w}{\partial x}& 0\\ 0& \frac{\partial w}{\partial y}\\ \frac{\partial w}{\partial y}& \frac{\partial w}{\partial x}\end{array}\right]\left\{\begin{array}{c}\frac{\partial w}{\partial x}\\ \frac{\partial w}{\partial y}\end{array}\right\}=\frac{1}{2}\left[A\right]\left\{\theta \right\}\mathrm{}$ (10)

where the vector $\left\{\theta \right\}$ with the derivatives of shell transverse deflection $w$ can be evaluated based on the element plate displacements $\left\{d_e^m\right\}$ and considering a bilinear interpolation in this field, as used in the Jacobian matrix $\left[J\right]$, Batoz and Tahar [¹²] and Rojaset al. [¹⁰], Equation 11:

$\left\{\theta \right\}=\left\{\begin{array}{c}\frac{\partial w}{\partial x}\\ \frac{\partial w}{\partial y}\end{array}\right\}={\left[J\right]}^{-1}\left\{\begin{array}{c}\frac{\partial w}{\partial \xi }\\ \frac{\partial w}{\partial \eta }\end{array}\right\}={\left[G\right]}_{(2\times 12)}{\left\{d_e^b\right\}}_{(12\times 1)}\mathrm{}$ (11)

Based on the vector $\left\{\theta \right\}$ components, it is possible to assemble the matrix $\left[A\right]$. Thereby, it is possible to note that these two elements depend on the nodal displacements $\left\{d_e^b\right\}$, and the multiplication between them results in a nonlinear relationship between $\left\{{\epsilon }_b^L\right\}$ and $\left\{d_e^b\right\}$. Calculating the variation of $\left\{{\epsilon }_b^L\right\}$ with respect to $\left\{d_e^m\right\}$ we obtain the matrix $\left[B_b^L\right]$:

$\partial \left\{{\epsilon }_b^L\right\}=\frac{1}{2}\partial \left(\left[A\right]\right)\left\{\theta \right\}+\frac{1}{2}\left[A\right]\partial \left(\left\{\theta \right\}\right)=\left[A\right]\partial \left(\left\{\theta \right\}\right)=\left[A\right]\left[G\right]\partial \left\{d_e^b\right\}=\left[B_b^L\right]\partial \left\{d_e^b\right\}$ (12)

Unlike the matrices $\left[B_m\right]$, $\left[B_b^0\right]$ and $\left[G\right]$ which remain constants through the analysis, the matrix $\left[B_b^L\right]$ depends on the nodal displacements $\left\{d_e^b\right\}$ and needs to be updated in each iteration.

The incremental kinematic matrix $\left[B\right]$, which differs from the kinematic matrix $\left[B\text{'}\right]$ due to the problem linearization process, can be written as:

${\left[B\right]}_{(3\times 24)}=\left[\begin{array}{cc}{\left[B_m\right]}_{(3\times 12)}& z_L{\left[B_b^0\right]}_{(3\times 12)}+{\left[B_b^L\right]}_{(3\times 12)}\end{array}\right]$ (13)

2.6 Element stiffness matrix

Through the Total Lagrangian Formulation, the element stiffness matrix $\left[k_e\right]$ is defined by two components $\left[k_e^L\right]$ and $\left[k_e^{\sigma }\right]$ which consider, respectively, the large displacements and the stresses acting on the structure:

$\left[k_e\right]\mathrm{}=\left[k_e^L\right]+\left[k_e^{\sigma }\right]={\int }_{V_e}^{\mathrm{}}{\left[B\right]}^T\left[D\right]\left[B\right]d{V}_e+{\int }_{V_e}^{\mathrm{}}{\left[G\right]}^T\left[M\right]\left[G\right]d{V}_e$ (14)

where $\left[D\right]$ is the material tangent stiffness matrix, being adopted in its place, in this paper, the material secant stiffness matrix, and $\left[M\right]$ is a matrix defined according to the acting stresses.

The element volume $V_e$ integrals are evaluated numerically. In the element plane, it is adopted the Gauss quadrature [¹⁴]. As the shell thickness is discretized in $n_c$ concrete layers and $n_s$ steel layers, in the integral along the thickness, it is used a mixed approach between the presented by Zhang et al. [¹⁵], [¹⁶], Barrales [⁵] and Vasilescu [¹³], which considers the sum of each layer individual contribution.

Thereby, using Equations 13 and 14, the matrix $\left[k_e^L\right]$ is evaluated as:

$\left[k_e^L\right]=\sum _{k=1}^{n_g}{p}_k\left|{\left[J\right]}_k\right|\left(\sum _{i=1}^{n_c}\left[\begin{array}{cc}{\left[k_{c.mm}^L\right]}_{k,i}& {\left[k_{c.mb}^L\right]}_{k,i}\\ {\left[k_{cbm}^L\right]}_{k,i}& {\left[k_{c.bb}^L\right]}_{k,i}\end{array}\right]+\sum _{j=1}^{n_s}\left[\begin{array}{cc}{\left[k_{s.mm}^L\right]}_{k,j}& {\left[k_{s.mb}^L\right]}_{k,j}\\ {\left[k_{s.bm}^L\right]}_{k,j}& {\left[k_{s.bb}^L\right]}_{k,j}\end{array}\right]\right)$ (15)

The parameters $p_k$ and $\left|{\left[J\right]}_k\right|$ represent, respectively, the integration weights of the $n_g$ Gauss points, and the determinant of the corresponding Jacobian matrix ${\left[J\right]}_k$.

The concrete layers contributions to the $\left[k_e^L\right]$ matrix, evaluated for each $k$ Gauss point and each $i$ concrete layer, are given by:

${\left[k_{c.mm}^L\right]}_{k,i}={z_{c1}}_i{\left[B_m^0\right]}_i^T{\left[D_c\right]}_{k,i}{\left[B_m^0\right]}_i$ (16a)

${\left[k_{c.mb}^L\right]}_{k,i}={z_{c2}}_i{\left[B_m^0\right]}_i^T{\left[D_c\right]}_{k,i}{\left[B_b^0\right]}_i+{z_{c1}}_i{\left[B_m^0\right]}_i^T{\left[D_c\right]}_{k,i}{\left[B_b^L\right]}_i$ (16b)

${\left[k_{c.bm}^L\right]}_{k,i}={z_{c2}}_i{\left[B_b^0\right]}_i^T{\left[D_c\right]}_{k,i}{\left[B_m^0\right]}_i+{z_{c1}}_i{\left[B_m^L\right]}_i^T{\left[D_c\right]}_{k,i}{\left[B_m^0\right]}_i$ (16c)

$\begin{array}{c}{\left[k_{c.bb}^L\right]}_{k,i}={z_{c3}}_i{\left[B_b^0\right]}_i^T{\left[D_c\right]}_{k,i}{\left[B_b^0\right]}_i+{z_{c2}}_i{\left[B_b^0\right]}_i^T{\left[D_c\right]}_{k,i}{\left[B_b^L\right]}_i+{z_{c2}}_i{\left[B_b^L\right]}_i^T{\left[D_c\right]}_{k,i}{\left[B_b^0\right]}_i+{z_{c1}}_i{\left[B_b^L\right]}_i^T{\left[D_c\right]}_{k,i}{\left[B_b^L\right]}_i\\ \end{array}$ (16d)

The parameter ${\left[D_c\right]}_{k,i}$ is the material secant stiffness matrix, evaluated at the $i$ concrete layer in the $k$ Gauss points, discussed in section 3. The terms ${z_{c1}}_i$, ${z_{c2}}_i$ and ${z_{c3}}_i$ arise from the numerical integration of the layer coordinate parameter $z_L$ in the matrices $\left[B_b^0\right]$, component of $\left[B\right]$ matrix, and are calculated in a discrete way, for each concrete layer, as shown in Equation 17, where ${z_c}_i^{top}$ and ${z_c}_i^{bot}$ represent, respectively, the analyzed layer top and bottom coordinates.

${z_{c1}}_i=\left[{z_c}_i^{top}-{z_c}_i^{bot}\right]$ (17a)

${z_{c2}}_i=\frac{1}{2}\left[{\left({z_c}_i^{top}\right)}^2-{\left({z_c}_i^{bot}\right)}^2\right]$ (17b)

${z_{c3}}_i=\frac{1}{3}\left[{\left({z_c}_i^{top}\right)}^3-{\left({z_c}_i^{bot}\right)}^3\right]$ (17c)

The steel layers contributions to the $\left[k_e^L\right]$ matrix, evaluated for each $k$ Gauss point and each $j$ steel layer, are given by:

${\left[k_{s.mm}^L\right]}_{k,j}={{\rho }_s}_j{t_s}_j{\left[B_m^0\right]}_j^T{\left[D_s\right]}_{k,j}{\left[B_m^0\right]}_j$ (18a)

${\left[k_{s.mb}^L\right]}_{k,j}={{\rho }_s}_j{t_s}_j{z_s}_j{\left[B_m^0\right]}_j^T{\left[D_s\right]}_{k,j}{\left[B_b^0\right]}_j+{{\rho }_s}_j{t_s}_j{\left[B_m^0\right]}_j^T{\left[D_s\right]}_{k,j}{\left[B_b^L\right]}_j$ (18b)

${\left[k_{s.bm}^L\right]}_{k,j}={{\rho }_s}_j{t_s}_j{z_s}_j{\left[B_b^0\right]}_j^T{\left[D_s\right]}_{k,j}{\left[B_m^0\right]}_j+{{\rho }_s}_j{t_s}_j{\left[B_m^L\right]}_j^T{\left[D_s\right]}_{k,j}{\left[B_m^0\right]}_j$ (18c)

$\begin{array}{c}{\left[k_{s.bb}^L\right]}_{k,j}={{\rho }_s}_j{t_s}_j{z_s}_j^2{\left[B_b^0\right]}_j^T{\left[D_s\right]}_{k,j}{\left[B_b^0\right]}_j+{{\rho }_s}_j{t_s}_j{z_s}_j{\left[B_b^0\right]}_j^T{\left[D_s\right]}_{k,j}{\left[B_b^L\right]}_j+\\ \mathrm{}{{\rho }_s}_j{t_s}_j{z_s}_j{\left[B_b^L\right]}_j^T{\left[D_s\right]}_{k,j}{\left[B_b^0\right]}_j+{{\rho }_s}_j{t_s}_j{\left[B_b^L\right]}_j^T{\left[D_s\right]}_{k,j}{\left[B_b^L\right]}_j\end{array}$ (18d)

The parameter ${\left[D_s\right]}_{k,i}$ is the material secant stiffness matrix, evaluated at the $j$ steel layer in the $k$ Gauss points, also discussed in section 3. The variable $z_s$ refers to the position of the $j$ steel layer axis. On the other hand, ${t_s}_j$ and ${{\rho }_s}_j$ are the corresponding layer thickness and reinforcement ratio.

Similarly, the matrix $\left[k_e^{\sigma }\right]$ can be defined numerically as:

$\left[k_e^{\sigma }\right]={\int }_V^{\mathrm{}}{\left[G\right]}^T\left[M\right]\left[G\right]dV=\left[\begin{array}{cc}{\left[0\right]}_{(12\times 12)}& {\left[0\right]}_{(12\times 12)}\\ {\left[0\right]}_{(12\times 12)}& \sum _{k=1}^{n_g}{p}_k\left|{\left[J\right]}_k\right|\left(\sum _{i=1}^{n_c}{\left[k_c^{\sigma }\right]}_{k,i}+\sum _{j=1}^{n_s}{\left[k_s^{\sigma }\right]}_{k,j}\right)\end{array}\right]$ (19a)

${\left[k_c^{\sigma }\right]}_{k,i}={\left[G\right]}_i^T{\left[M_c\right]}_{k,i}{\left[G\right]}_i={{z_{c1}}_i\left[G\right]}_i^T\left[\begin{array}{cc}{{\sigma }_{xx}^c}_i& {{\tau }_{xy}^c}_i\\ {{\tau }_{xy}^c}_i& {{\sigma }_{yy}^c}_i\end{array}\right]{\left[G\right]}_i$ (19b)

${\left[k_s^{\sigma }\right]}_{k,j}={\left[G\right]}_j^T{\left[M_s\right]}_{k,j}{\left[G\right]}_j={{{\rho }_s}_j{t_s}_j\left[G\right]}_j^T\left[\begin{array}{cc}{{\sigma }_{xx}^s}_j& {{\tau }_{xy}^s}_j\\ {{\tau }_{xy}^s}_j& {{\sigma }_{yy}^s}_j\end{array}\right]{\left[G\right]}_j$ (19c)

where ${\sigma }_{xx}$, ${\sigma }_{yy}$ and ${\tau }_{xy}$ represent the stresses acting on the concrete and steel layers. Note that the matrix $\left[k_e^{\sigma }\right]$ contributes to $\left[k_e\right]$ only in the region associated with the plate degrees of freedom, a consequence of considering nonlinear deformations only in the nonlinear plate strain component $\left\{{\epsilon }_b^L\right\}$.

2.7 Internal forces vector

The internal forces vector $\left\{F_{int}\right\}$ unlike the external forces vector $\left\{F_{ext}\right\}$ is not constant. It depends on the structure nodal displacement vector $\left\{d\right\}$ and must be updated at each iteration. The corresponding element contribution $\left\{f_e\right\}$ can be defined by Equation 20a and calculated numerically by Equation 20b.

$\left\{f_e\right\}={\int }_{V_e}^{\mathrm{}}{\left[B\right]}^T\left\{\sigma \right\}d{V}_e={\int }_{V_e}^{\mathrm{}}{\left[B\right]}^T\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\tau }_{xy}\end{array}\right\}d{V}_e$ (20a)

$\left\{f_e\right\}=\sum _{k=1}^{n_g}{p}_k\left|{\left[J\right]}_k\right|\left(\sum _{i=1}^{n_c}{\left[B_c\right]}_{k,i}^T\left\{\begin{array}{c}{{\sigma }_{xx}^c}_i\\ {{\sigma }_{yy}^c}_i\\ {{\tau }_{xy}^c}_i\end{array}\right\}+\sum _{j=1}^{n_s}{\left[B_s\right]}_{k,j}^T\left\{\begin{array}{c}{{\sigma }_{xx}^s}_j\\ {{\sigma }_{yy}^s}_j\\ {{\tau }_{xy}^s}_j\end{array}\right\}\right)$ (20b)

where the matrices ${\left[B_c\right]}_{k,i}$ and ${\left[B_s\right]}_{k,j}$ can be calculated as:

${\left[B_c\right]}_{k,i}=\left[\begin{array}{cc}{{z_{c1}}_i\left[B_m\right]}_k& {z_{c2}}_i{\left[B_b^0\right]}_k+{{z_{c1}}_i\left[B_b^L\right]}_k\end{array}\right]$ (21a)

${\left[B_s\right]}_{k,j}=\left[\begin{array}{cc}{{{\rho }_s}_j{t_s}_j\left[B_m\right]}_k& {{\rho }_s}_j{t_s}_j{z_s}_j{\left[B_b^0\right]}_k+{{{\rho }_s}_j{t_s}_j\left[B_b^L\right]}_k\end{array}\right]$ (21b)

2.8 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.

Figure 2 details the implemented program. As this figure illustrates, in the first load step first iteration, the concrete constitutive matrix can be initialized based on a linear-elastic model, Equation. (4.22), where ${\mathrm{\upsilon }}_0$ is the initial Poisson ratio and ${\mathrm{E}}_{\mathrm{c}}$ is the Modulus of elasticity of concrete.

Figure 2
Nonlinear analysis procedure

$\left[D_{linear}\right]=\frac{E_c}{1-{\upsilon }_0^2}\left[\begin{array}{ccc}1& {\upsilon }_0& 0\\ {\upsilon }_0& 1& 0\\ 0& 0& \frac{1-{\upsilon }_0}{2}\end{array}\right]$ (22)

3 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 [²] [³], [¹⁷]. In addition, this section also details the formulation of the concrete $\left[D_c\right]$ and steel $\left[D_s\right]$ layers secant stiffness matrices, according to Vecchio [³]. 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) [²].

3.1 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) [¹⁷], which through the softening coefficient ${\beta }_d$, models the material strength loss due to transversal tension, Equation 24.

${\sigma }_{cc}=\left\{\begin{array}{c}-f_p\left[2\left(\frac{{\epsilon }_{cc}}{{\epsilon }_p}\right)-{\left(\frac{{\epsilon }_{cc}}{{\epsilon }_p}\right)}^2\right],\mathrm{}if\mathrm{}\left|{\epsilon }_{cc}\right|<2\left|{\epsilon }_p\right|\\ \mathrm{}0,\mathrm{}if\mathrm{}\left|{\epsilon }_{cc}\right|<2\left|{\epsilon }_p\right|\end{array}\right.$ (23)

$f_p={\beta }_d{f}_c$ (24a)

${\epsilon }_p={{\beta }_d\epsilon }_0$ (24b)

${\beta }_d=\frac{1}{1+C_s{C}_d}\le 1$ (24c)

$C_d=\left\{\begin{array}{c}0,ifmin(-{\epsilon }_{c1}/{\epsilon }_{c2},400)\ge 0.28\\ 0.35{\left[\mathrm{m}\mathrm{i}\mathrm{n}(-{\epsilon }_{c1}/{\epsilon }_{c2},\mathrm{}400)-0.28\right]}^{0.80},ifmin(-{\epsilon }_{c1}/{\epsilon }_{c2},400)\ge 0.28\end{array}\right.$ (24d)

where ${\sigma }_{cc}$ and ${\epsilon }_{cc}$ are, respectively, the average compressive stress and strain. The variable $f_c$ is the cylinder compressive strength (at 28 days), and ${\epsilon }_0$ is the corresponding peak strain. The factor $C_d$ represents the influence of the relationship between the principal tensile ${\epsilon }_{c1}$ and compression ${\epsilon }_{c2}$ strains in concrete. The factor $C_s$ is equal to 0.55 when the slip deformations in the cracks are considered in the model (subsection 3.6), and 1.0 otherwise.

3.2 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. [¹⁷], 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 [³], the magnitude of the average principal tensile stress ${\sigma }_{ct}$ must be limited by the remaining steel resistant capacity ${\sigma }_{ct}^{max}$.

${\sigma }_{ct}=\left\{\begin{array}{c}{E}_c{\epsilon }_{ct},\mathrm{}if\mathrm{}{\epsilon }_{ct}<{\epsilon }_{cr}\\ \frac{f_{cr}}{1+c_t{\epsilon }_{ct}},\mathrm{}if\mathrm{}{\epsilon }_{ct}\ge {\epsilon }_{cr}\end{array}\le {\sigma }_{ct}^{max}\right.$ (25a)

$c_t=\frac{2.2}{4\frac{{\rho }_x}{{\varphi }_x}\left|\mathrm{c}\mathrm{o}\mathrm{s}(\theta -{\alpha }_x)\right|+4\frac{{\rho }_y}{{\varphi }_y}\left|\mathrm{c}\mathrm{o}\mathrm{s}(\theta -{\alpha }_y)\right|}$ (25b)

${\sigma }_{ct}^{max}={\rho }_x(f_{scrx}-f_{sx}){\mathrm{c}\mathrm{o}\mathrm{s}(\theta -{\alpha }_x)}^2+{\rho }_y(f_{scry}-f_{sy}){\mathrm{c}\mathrm{o}\mathrm{s}(\theta -{\alpha }_y)}^2$ (25c)

where ${\epsilon }_{ct}$ is the concrete average principal tensile strain. The parameters $f_{cr}$ and ${\epsilon }_{cr}$ are, respectively, the crack stress and strain. The coefficient $c_t$ refers to the influence of the reinforcement on the stiffening, and it is obtained based on: the steel ratios ${\rho }_x$ and ${\rho }_y$; the reinforcement rebar nominal diameter ${\varphi }_x$ and ${\varphi }_y$; and the angles $\theta$, ${\alpha }_x$ and ${\alpha }_y$ that illustrate the principal system direction and the reinforcements orientation. The parameters $f_{sx}$ and $f_{sy}$ are the reinforcements average stresses, while $f_{scrx}$ and $f_{scry}$ represent these same parameters evaluated at the crack surface. The Modulus of elasticity of concrete $E_c$ illustrates the ratio between $f_{cr}$ and ${\epsilon }_{cr}$, and can be estimated as ${2f}_c/\left|{\epsilon }_0\right|$.

When the reinforcement properties are directly assigned to the finite element, like in reinforced concrete membranes analysis programs [¹⁸], 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 [⁴] presented a solution to this problem, based on CEB-FIP [¹⁹], 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 ${\mathrm{\sigma }}_{\mathrm{c}1}$ and ${\mathrm{\sigma }}_{\mathrm{c}2}$ 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.

3.3 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.

$f_s=\left\{\begin{array}{c}\frac{f_y}{{\epsilon }_{sy}}{\epsilon }_s,if\left|{\epsilon }_s\right|<{\epsilon }_{sy}\\ f_{y}sign\left({\epsilon }_s\right),if\left|{\epsilon }_s\right|\ge {\epsilon }_{sy}\end{array}\right.$ (26)

$f_s=\left\{\begin{array}{c}\frac{f_{sy}}{{\epsilon }_{sy}}{\epsilon }_s,if\left|{\epsilon }_s\right|<{\epsilon }_{sy}\\ \left[f_y+\frac{f_u-f_y}{{\epsilon }_{su}-{\epsilon }_{sy}}({\epsilon }_s-{\epsilon }_{sy})\right]sign\left({\epsilon }_s\right),if\left|{\epsilon }_s\right|\ge {\epsilon }_{sy}\end{array}\right.$ (27)

where $f_y$ and ${\epsilon }_{sy}$ are the yield stress and strain and $f_s$ and ${\epsilon }_s$ are the reinforcement average stress and strain. The parameters $f_u$ and ${\epsilon }_{su}$ 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.

3.4 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 [¹⁸], based on the work of Vecchio [²⁰], Kupfer et al. [²¹], Richart et al. [²²] and Wong et al. [¹⁷], using the enhancement factors $K_{c1}$ and $K_{c2}$:

$K_{c1}\left({\sigma }_{c2}\right)=1+0.92\left(-\frac{{\sigma }_{c2}}{f_c}\right)-0.76{\left(-\frac{{\sigma }_{c2}}{f_c}\right)}^2+4.1\frac{{\rho }_{sz}{f}_{sz}}{f_c}$ (28a)

$K_{c2}\left({\sigma }_{c1}\right)=1+0.92\left(-\frac{{\sigma }_{c1}}{f_c}\right)-0.76{\left(-\frac{{\sigma }_{c1}}{f_c}\right)}^2+4.1\frac{{\rho }_{sz}{f}_{sz}}{f_c}$ (28b)

where ${\rho }_{sz}$ and $f_{sz}$ are the reinforcement ratio and its stress in the out-of-plane direction. The stress $f_{sz}$ can be evaluated by applying the concrete strain in the corresponding direction ${\epsilon }_{cz}$ 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:

${\epsilon }_{cz}=\frac{E_{cn}}{E_{cn}+{\rho }_{sz}{E}_{sz}}(-v_{12}{\epsilon }_{c2}-v_{21}{\epsilon }_{c1})$ (29)

where $E_{cn}$ and $E_{sz}$ are the concrete and steel modulus of elasticity in the out-of-plane direction, where $E_{cn}$ is considered equal to $\left|2f_c/{\epsilon }_0\right|$. The parameter $v_{12}$ represents the Poisson ratio that relates the strain in the 1-direction due to the stress in 2-direction [²⁰]. The parameter $v_{21}$ is defined similarly. The model adopted for calculating the Poisson coefficients $v_{12}$ and $v_{21}$ is detailed in the following subsection. The factors $K_{c1}$ and $K_{c2}$ are used to determine the peak stress and strain in the principal system (1-2):

$f_{p1}=K_{c1}{f}_c$ (30a)

$f_{p2}=K_{c2}{f}_c$ (30b)

${\epsilon }_{p1}={(3K}_{c1}-2){\epsilon }_0$ (31a)

${\epsilon }_{p2}={(3K}_{c2}-2){\epsilon }_0$ (31b)

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:

$f\left({\sigma }_{c1},{\sigma }_{c2}\right)=\left\{\begin{array}{c}{\sigma }_{c1}-{\sigma }_{cc}(f_{p1},{\epsilon }_{p1},{\epsilon }_{c1})\\ {\sigma }_{c2}-{\sigma }_{cc}(f_{p2},{\epsilon }_{p2},{\epsilon }_{c2})\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\end{array}\right\}$ (32)

This solution strategy was implemented in the developed program, where optimization functions from SciPy [²³] 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).

3.5 Poisson ratio and lateral strains

According to Vecchio [²⁰], 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 $v_0$ increase, also disregards this parameter when the concrete presents, in the transverse direction, a principal tensile strain greater than the cracking strain ${\epsilon }_{cr}$.

$v_{12}=\left\{\begin{array}{c}0,if{\epsilon }_{c2}\ge 0and{\epsilon }_{c2}\ge {\epsilon }_{cr}\\ {\mathrm{}v}_0,if{\epsilon }_{c2}\ge 0and{\epsilon }_{c2}<{\epsilon }_{cr}\\ v_0,if{\epsilon }_{c2}<0and\left|{\epsilon }_{c2}\right|<\left|{\epsilon }_0\right|/2\\ v_0\left[1+1.5{\left(\frac{2{\epsilon }_{c2}}{{\epsilon }_0}-1\right)}^2\right]\ngtr 0.5,if{\epsilon }_{c2}<0and\left|{\epsilon }_{c2}\right|\ge \left|{\epsilon }_0\right|/2\end{array}\right.$ (33)

The concrete principal lateral strain in 1-direction ${\epsilon }_{c01}$ is given by Equation 34. It is important to note that the corresponding strain in 2-direction ${\epsilon }_{c02}$ and its Poisson ratio $v_{21}$ can be obtained through Equations 33 and (34) switching the indexes.

${\epsilon }_{c01}=-v_{21}{\epsilon }_{c2}$ (34)

The concrete principal lateral strain vector $\left\{{\epsilon }_{c0}^{1-2}\right\}$ can be transformed to the corresponding cartesian system vector $\left\{{\epsilon }_{c0}^{}\right\}$, Equation 35, using the rotation matrix $\left[T\right]$, Equation 36.

$\left\{{\epsilon }_{c0}^{\mathrm{}}\right\}=\left[T(-\theta )\right]\left\{{\epsilon }_{c0}^{1-2}\right\}=\left[T(-\theta )\right]{\left\{\begin{array}{ccc}{\epsilon }_{c01}& {\epsilon }_{c02}& 0\end{array}\right\}}^T$ (35)

$\left[T\left(\theta \right)\right]=\left[\begin{array}{ccc}{\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)}^2& {\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)}^2& \mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\\ {\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)}^2& {\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)}^2& -\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\\ -2\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)& 2\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)& {\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)}^2-{\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)}^2\end{array}\right]$ (36)

3.6 Slip strain model

The Disturbed Stress Field Model (DSFM) is an extension of the well-known Modified Compression Field Theory (MCFT) [²], which admits disagreements between the stress and strain principal systems, through the consideration of crack slip strains. The DSFM was proposed by Vecchio [³] to solve some MCFT drawbacks in calculate the strength and stiffness of high or low reinforcement ratio elements.

The crack surface shear stress $v_c$ can be evaluated applying equilibrium in the reinforced concrete element:

$\begin{array}{c}{v}_c={\rho }_x\left(f_{scrx}-f_{sx}\right)\cos \left(\theta -{\alpha }_x\right)\sin \left(\theta -{\alpha }_x\right)+{\rho }_y\left(f_{scry}-f_{sy}\right)\cos \left(\theta -{\alpha }_y\right)\sin \left(\theta -{\alpha }_y\right)\\ \end{array}$ (37)

According to Vecchio [³] due to this shear stress there is a rigid body local slip along the crack (slip displacement ${\delta }_s$) which causes slip strains $\left\{{\epsilon }^s\right\}$. These strains must be considered in the model additionally to the principal strains related to the material constitutive response. According to Vecchio [³], the displacement ${\delta }_s$ can be calculated as:

${\delta }_s=\frac{v_c}{1.8w_{cr}^{-0.8}+(0.234w_{cr}^{-0.707}-0.20)f_{cc}}$ (38)

where $w_{cr}^{}$ represents the average crack width, which can be estimated from the average crack spacing $s_{cr}^{}$, Equation 39, based on the nominal crack spacings in x and y ($s_x^{}=s_y^{}=50mm$).

$w_{cr}={\epsilon }_{c1}{s}_{cr}={\epsilon }_{c1}\frac{1}{\sin \left(\theta \right)/s_x+\cos \left(\theta \right)/s_y}$ (39)

When the average crack width is greater than or equal to 5 mm, the concrete principal compressive stress ${\sigma }_{c2}$ is considered equal to zero. The parameter $f_{cc}$ is cubic compressive strength, adopted as: $f_{cc}=f_c/0.85$. The crack slip shear strain ${\gamma }_s$ is evaluated as the ratio between ${\delta }_s$ and $s_{cr}$.

Based on Mohr's circle coordinate transformations, the slip strain vector in the Cartesian system $\left\{{\epsilon }^s\right\}$ can be written as:

$\left\{{\epsilon }^s\right\}\mathrm{}={\gamma }_s\left\{\begin{array}{c}-0.5\sin \left(2\theta \right)\\ 0.5\sin \left(2\theta \right)\\ \cos \left(2\theta \right)\end{array}\right\}$ (40)

The formulation described above represents an overview of the slip strain vector $\left\{{\epsilon }^s\right\}$ calculation procedure. However, it is important to present some additional details about the shear stresses along the crack surfaces $v_c$. Although this parameter can be obtained by Equation 37, according to Silva [¹⁸], it should be limited to the maximum shear stress that can be resisted on the crack by aggregate interlock, Equation 41.

$v_c\le v_c^{max}=\frac{0.18\sqrt{f_c}}{0.31+24w_{cr}/(a_g+26)}$ (41)

where $a_g$ is the aggregate size (adopted as 25mm). Furthermore, according to Equation 37, to evaluate the shear stress $v_c$ it is previously necessary to calculate the reinforcement stresses in the crack surface $f_{scrx}$ and $f_{scry}$. 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, ${\epsilon }_{scrx}$ and ${\epsilon }_{scry}$, 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 [³] proposes that the steel strains, in the cracks ${\epsilon }_{scrx}$ and ${\epsilon }_{scry}$, should be determined through the sum of the corresponding average strain ${\epsilon }_{sx}$ and ${\epsilon }_{sy}$ and the local incremental strain contribution in 1-direction, ${\mathrm{\Delta }\epsilon }_{1cr}$:

${\epsilon }_{scrx}={\epsilon }_{sx}+{\mathrm{\Delta }\epsilon }_{1cr}{\mathrm{c}\mathrm{o}\mathrm{s}(\theta -{\alpha }_x)}^2$ (42a)

${\epsilon }_{scry}={\epsilon }_{sy}+{\mathrm{\Delta }\epsilon }_{1cr}{\mathrm{c}\mathrm{o}\mathrm{s}(\theta -{\alpha }_y)}^2$ (42b)

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 ${\mathrm{\Delta }\epsilon }_{1cr}$, Equation 43. Again, it was used SciPy [²³] scientific computing library optimization functions.

$f\left({\mathrm{\Delta }\epsilon }_{1cr}\right)={\sigma }_{c1}^{\mathrm{}}-{\rho }_x\left(f_{scrx}-f_{sx}\right){\cos \left(\theta -{\alpha }_x\right)}^2-{\rho }_y\left(f_{scry}-f_{sy}\right){\cos \left(\theta -{\alpha }_y\right)}^2=0$ (43)

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 ${\epsilon }_{c1}$ is greater than the crack strain ${\epsilon }_{cr}$. Otherwise, the slip strain vector $\left\{{\epsilon }^s\right\}$ can be computed as a null vector.

3.7 Material secant stiffness matrices

In the previous subsections, it was presented the lateral strain vector $\left\{{\epsilon }_{c0}^{}\right\}$ and the slip strain vector $\left\{{\epsilon }^s\right\}$ calculation procedure. In order to ensure the concrete secant stiffness matrix symmetry and the associated benefits, according to Vecchio [³] and Silva [¹⁸], the concrete total strain vector $\left\{\epsilon \right\}$ is defined by three distinct components: $\left\{{\epsilon }_c\right\}$, $\left\{{\epsilon }_{c0}^{}\right\}$ and $\left\{{\epsilon }^s\right\}$, where, $\left\{{\epsilon }_c\right\}$ represents the elastic strains due to stress, which can be evaluated, as:

$\left\{{\epsilon }_c\right\}=\left\{\epsilon \right\}-\left\{{\epsilon }_{c0}^{\mathrm{}}\right\}-\left\{{\epsilon }^s\right\}{=\left\{\begin{array}{ccc}{\epsilon }_{cx}& {\epsilon }_{cy}& {\gamma }_{cxy}\end{array}\right\}}^T$ (44)

Once the vector $\left\{{\epsilon }_c\right\}$ are obtained, it is possible to estimate the inclination of the principal strains $\theta$, Equation 45, and the concrete principal strains:

$\theta =0.5{tan}^{-1}\left(\frac{{\gamma }_{cxy}}{{\epsilon }_{cx}-{\epsilon }_{cy}}\right)$ (45)

${\epsilon }_{c1},{\epsilon }_{c2}=\frac{{\epsilon }_{cx}+{\epsilon }_{cy}}{2}\pm \frac{1}{2}\sqrt{{({\epsilon }_{cx}-{\epsilon }_{cy})}^2+{{\gamma }_{cxy}}^2}$ (46)

Through Equations 44 and (46) and the models presented in the previous subsections, it is possible to see the nonlinear relationship between the vectors $\left\{{\epsilon }_c\right\}$, $\left\{{\epsilon }_{c0}^{}\right\}$ and $\left\{{\epsilon }^s\right\}$. In the present paper, this problem was solved iteratively, as illustrated in Figure 3. In the cracked reinforced concrete element, the concrete stress vector $\left\{{\sigma }_c\right\}$ can be related to the corresponding strain vector $\left\{{\epsilon }_c\right\}$ as:

$\left\{{\sigma }_c\right\}\mathrm{}=\left[D_c\right]\left\{{\epsilon }_c\right\}$ (47)

where the concrete layer secant stiffness matrix $\left[D_c\right]$ in the Cartesian system is obtained by a coordinate transformation of the corresponding matrix in the principal system $\left[D_c^{1-2}\right]$, which is defined based on the concrete secant moduli ${\stackrel{-}{E}}_{c1}$ and ${\stackrel{-}{E}}_{c2}$, Equations 48 and 49.

Figure 3
Concrete layer secant constitutive model

$\left[D_c\right]\mathrm{}={\left[T\left(\theta \right)\right]}^T\left[D_c^{1-2}\right]\left[T\left(\theta \right)\right]={\left[T\left(\theta \right)\right]}^T\left[\begin{array}{ccc}{\stackrel{-}{E}}_{c1}& 0& 0\\ 0& {\stackrel{-}{E}}_{c2}& 0\\ 0& 0& \frac{{\stackrel{-}{E}}_{c1}{\stackrel{-}{E}}_{c2}}{{\stackrel{-}{E}}_{c1}+{\stackrel{-}{E}}_{c2}}\end{array}\right]\left[T\left(\theta \right)\right]$ (48)

${\stackrel{-}{E}}_{c1}={\sigma }_{c1}/{\epsilon }_{c1}\mathrm{}{\stackrel{-}{E}}_{c2}={\sigma }_{c2}/{\epsilon }_{c2}$ (49)

Like Equation 48, the steel secant stiffness matrix $\left[D_s\right]$ (in $x$ or $y$-direction) can be calculated as shown in Equation 50.

$\left[D_s\right]\mathrm{}={\left[T\left({\alpha }_s\right)\right]}^T\left[D_s^{\mathrm{\text{'}}}\right]\left[T\left({\alpha }_s\right)\right]={\left[T\left({\alpha }_s\right)\right]}^T\left[\begin{array}{ccc}{f}_s/{\epsilon }_s& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\left[T\left({\alpha }_s\right)\right]$ (50)

where the angle ${\alpha }_{sj}$ represents the reinforcement direction. Finally, the reinforcement layer stress vector $\left\{{\sigma }_s\right\}$ is calculated by:

$\left\{{\sigma }_s\right\}\mathrm{}=\left[T\left({\alpha }_s\right)\right]\left\{{\sigma }_s^{\mathrm{\text{'}}}\right\}=\left[T\left({\alpha }_s\right)\right]{\left\{\begin{array}{ccc}{f}_s& 0& 0\end{array}\right\}}^T$ (51)

Figures 3 and 4 detail the constitutive models’ implementation in the developed computer program, to obtain the secant stiffness matrices $\left[D_c\right]$ and $\left[D_s\right]$, and the stress vectors $\left\{{\sigma }_c\right\}$ and $\left\{{\sigma }_s\right\}$, on each layer.

Figure 4
Reinforcement layer secant constitutive model

4 PROGRAM VALIDATION AND DISCUSSIONS

This section presents the developed computer program validation through comparison with experimental and some numerical results [¹], [⁴], 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 [²⁴], a Python module called PyEVTK [²⁵] 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.

4.1 Cervenka deep beam

Initially, to evaluate the program performance in material nonlinear membrane problems, the deep beam W2 tested by Cervenka [²⁶], 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 $tol$, 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 (y-displacement at the deep beam bottom midpoint), the experiment and a literature numerical response [²⁷]. 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.

Figure 5
Cervenka deep beam

4.2 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 [²⁸] were analyzed. The specimens’ characteristics, in addition to the results obtained, in comparison with literature numerical solutions [¹], [⁴], 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 $tol$ 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.

Figure 6
Polak shells

4.3 Geometric nonlinear plates analysis

Three rectangular linear-elastic plates, subjected to a uniform load $q$, presented by Figueiras [⁷], were analyzed to verify the geometric nonlinear model implemented. The plates span $l$ were equal to 6m and its thickness $h$ was assumed as 0.15m. The material elastic modulus $E$ and Poisson ratio were adopted, respectively, as $30kN{/m}^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 $tol$ was equal to ${10}^{-5}$. Figure 7 shows the load parameter ($q.l^4/E.h^4$) versus central displacement parameter ($w/h$) curves obtained in this study, where it is possible to observe good accuracy when compared with literature analytical [²⁹] and numerically [⁷] solutions. The processing time was about 12 minutes.

Figure 7
Geometrically nonlinear plate analysis

4.4 Slender shear wall

The last problem analyzed was a slender shear wall. The structure geometry was adapted from the literature [⁶], [³⁰], [³¹] 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 $tol$ 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.

Figure 8
Slender shear wall

5 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 [⁵]. The material nonlinear analysis considered a secant stiffness approach, based on the Modified Compression Field Theory (MCFT) [²]. The element original formulation was expanded to also consider the problem geometric, through a Total Lagrangian Formulation [⁷]. 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 [⁷].

• 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.

• Financial support: None.
• Data Availability: The data that support the findings of this study are available from the corresponding author, J. R. B. Silva, upon reasonable request.
• How to cite: J. R. B. Silva and B. Horowitz, “Nonlinear analysis of reinforced concrete structures using thin flat shell elements,” Rev. IBRACON Estrut. Mater., vol. 15, no. 4, e15407, 2022, https://doi.org/10.1590/S1983-41952022000400007

REFERENCES

• ¹
  C. Luu, Y. Mo, and T. T. Hsu, "Development of csmm-based shell element for reinforced concrete structures," Eng. Struct., vol. 132, pp. 778–790, 2017.
• ²
  F. J. Vecchio and M. P. Collins, "The modified compression field theory for reinforced concrete elements subjected to shear," ACI Struct. J., vol. 83, no. 6, pp. 925–933, 1986.
• ³
  F. J. Vecchio, "Disturbed stress field model for reinforced concrete: formulation," J. Struct. Eng., vol. 126, no. 9, pp. 1070–1077, 2000.
• ⁴
  T. D. Hrynyk, “Behaviour and modelling of reinforced concrete slabs and shells under static and dynamic loads,” Ph.D. dissertation, Dept. Civ. Eng., Univ. Toronto, Toronto, Canada, 2013.
• ⁵
  F. R. Barrales, “Development of a nonlinear quadrilateral layered membrane element with drilling degrees of freedom and a nonlinear quadrilateral thin flat layered shell element for the modeling of reinforced concrete walls,” Ph.D. dissertation, Fac. USC Grad. Sch., Univ. Southern California, EUA, 2012.
• ⁶
  J. R. B. Silva and B. Horowitz, "Nonlinear finite element analysis of reinforced concrete shear walls," Rev. IBRACON Estrut. Mater., vol. 13, no. 6, pp. e13603, 2020.
• ⁷
  J. A. Figueiras, “Ultimate load analysis of anisotropic and reinforced concrete plates and shells,” Ph.D. dissertation, Dept. Civ. Eng., Univ. Coll. Swansea, Swansea, 1983.
• ⁸
  R. De Borst et al., Non-Linear Finite Element Analysis of Solids and Structures, 2nd ed. Chichester: Wiley, 2012.
• ⁹
  H. R. V. Goudarzi, “Nonlinear dynamic analysis of reinforced concrete frames under extreme loadings”, Ph.D dissertation, Sch. Civ. Environ. Eng., The Univ. New South Wales, Sydney, Australia, 2009.
• ¹⁰
  F. Rojas, J. Anderson, and L. Massone, "A nonlinear quadrilateral thin fat layered shell element for the modeling of reinforced concrete wall structures," Bull. Earthquake Eng., vol. 17, no. 12, pp. 6491–6513, 2019.
• ¹¹
  F. Rojas, J. Anderson, and L. Massone, "A nonlinear quadrilateral layered membrane element with drilling degrees of freedom for the modeling of reinforced concrete walls," Eng. Struct., vol. 124, pp. 521–538, 2016.
• ¹²
  J. L. Batoz and M. B. Tahar, "Evaluation of a new quadrilateral thin plate bending element," Int. J. Numer. Methods Eng., vol. 18, pp. 1655–1677, 1982.
• ¹³
  A. Vasilescu, “Analysis of geometrically nonlinear and softening response of thin structures by a new facet shell element,” M.S. thesis, Dept. Civ. Environ. Eng., Carleton Univ., Ottawa, Ontario, 2000.
• ¹⁴
  L. E. Vaz, Método dos Elementos Finitos em Análise de Estruturas Rio de Janeiro: Elsevier, 2011.
• ¹⁵
  Y. X. Zhang, M. A. Bradford, and R. I. Gilbert, "A layered cylindrical quadrilateral shell element for nonlinear analysis of RC plate structures," Adv. Eng. Softw., vol. 38, pp. 488–500, 2007.
• ¹⁶
  Y. X. Zhang, M. A. Bradford, and R. I. Gilbert, "A layered shear-flexural plate/shell element using Timoshenko beam functions for nonlinear analysis of reinforced concrete plates," Finite Elem. Anal. Des., vol. 43, pp. 888–900, 2007.
• ¹⁷
  P. S. Wong, H. Trommels, and F. J. Vecchio, VecTor2 and FormWorks User’s Manual. Technical Report, 2nd ed. Toronto: Dept. Civ. Eng., Univ. Toronto, 2012.
• ¹⁸
  L. M. T. Silva, “Análise não-linear de estruturas de concreto armado submetido ao estado plano de tensões,” M.S. thesis, Dept. Civ. Environ. Eng., Fed. Univ. Pernambuco, Recife, 2019.
• ¹⁹
  Comité EURO-International du Béton, Model Code for Concrete Structures – Design Code, 1990.
• ²⁰
  F. J. Vecchio, "Finite element modeling of concrete expansion and confinement," J. Struct. Eng., vol. 126, no. 9, pp. 1070–1077, 1992.
• ²¹
  H. Kupfer, H. K. Hilsdorf, and H. Rusch, "Behavior of concrete under biaxial stresses," J. Proc., vol. 66, no. 8, pp. 656–666, 1969.
• ²²
  F. Richart, A. Brandtzeag, and R. Brown, A Study of the Failure of Concrete Under Combined Compressive Stresses (Bulletin 185). Univ. Illinois Eng. Exp. Stn., p. 104, 1928.
• ²³
  SciPy. "Python-based ecosystem of open-source software for mathematics, science, and engineering." scipy.org (accessed Sept. 23, 2021).
• ²⁴
  J. Ahrens, B. Geveci, and C. Law, “ParaView: an end-user tool for large data visualization,” in Visualization Handbook, C. D. Hansen and C. R. Johnson, Eds., Amsterdam: Elsevier, 2005.
• ²⁵
  P. Herrera. “Evtk - Export Vtk.” https:// pypi.org/ project/ pyevtk/ (accessed Feb. 26, 2021).
• ²⁶
  V. Cervenka, “Inelastic finite element analysis of reinforced concrete panels,” Ph.D. dissertation, Univ. Colorado, Colorado, EUA, 1970.
• ²⁷
  F. J. Vecchio, "Reinforced concrete membrane element formulations," J. Struct. Eng., vol. 116, no. 3, pp. 730–750, 1990.
• ²⁸
  M. A. Polak and F. J. Vecchio, "Reinforced concrete shell elements subjected to bending and membrane loads," ACI Struct. J., vol. 91, no. 3, pp. 261–268, 1994.
• ²⁹
  S. Levy, Square Plate with Clamped Edges Under Normal Pressure Producing Large Deflections (DACA, Tech. Note 847). 1942.
• ³⁰
  R. L. S. França and A. E. Kimura, “Resultados de recentes pesquisas para o dimensionamento das armaduras longitudinal e transversal em pilares-parede,” in 9° Enc. Nac. Eng. Con. Estrut., 2006.
• ³¹
  A. E. Kimura, Cálculo de Pilares de Concreto Armado – Introdução, Visão Geral & Exemplos São Paulo: Assoc. Bras. Eng. Con. Estrut., 2016.

Edited by

Editors: Osvaldo Manzoli, Guilherme Aris Parsekian.

Publication Dates

• Publication in this collection
  28 Feb 2022
• Date of issue
  2022

History

• Received
  04 Oct 2021
• Accepted
  02 Jan 2022