An interface element for numerical analysis of flat plate/shell elements with deformable connection

Abstract

The use of numerical methods such as the finite elements method for solving structural analysis problems is becoming ever more efficient. An interface element capable of associating flat plate/shell element in three different combinations has been developed. In a structural analysis, when only the flat plate/shell’s finite elements are used in the discretization, some problems concerning dimensional variation of the transversal section, as well as overlapping of areas can occur. . The aforementioned problems can be solved by the developed interface elements that can also simulate a possible deformable connection that is existent in the association of materials with different characteristics, such as the composite steel-concrete elements. One of the applications of the finite elements developed is for the numerical simulation of composite steel-concrete elements, such as the composite beam formed by a reinforced concrete slab attached to a steel beam using a deformable connection. In this case, the concrete slab and the steel beam are discretized by the flat plate/shell element, and the deformable connection done by the interface elements. In the validation of the implemented elements, we used numerical and experimental results found in the literature, and analytical solutions considering the classical plate theory.

Keywords
Interface element; flat plate/shell element; composite beam; deformable connection

1 INTRODUCTION

In Brazil, the use of steel in civil construction is rapidly increasing, especially in the manufacturing of columns and beams for multi-story buildings. The presence of composite elements formed by the mechanical association between a concrete surface element and steel beams is becoming popular in civil construction. Usually, this mechanical association is made by connectors welded to the steel element and surrounded by the concrete element. The composite structural elements that are mostly used in the civil construction are the composite floors, composite beams and composite columns. Among them, the composite beam stands out, which is a simplification of the composite floor considering the concept of effective width of the concrete slab. Even in situations of isolated composite beams, the problem of normal tension variation along the width of the concrete slab can generate significant errors when evaluating this element using the beam theory. This project presents the formulation of a interface element that has the function of numerically simulating a deformable connection between structural elements that can be analyzed by the flat plate/shell elements, for example, the modeling of the composite steel-concrete beam problem can be done by using flat plate/shell elements for the concrete slab and the steel profile, and interface elements for the deformable connection between them.

Approximately 45 years ago, the finite interface elements have been studied by several researchers. These elements were initially developed to represent a thin layer of material or to simulate contact between distinct materials such as the soil-structure interaction studied in the work of Kaliakin and Li (1995) . The first project that was found in literature on interface elements was proposed by Goodman et al. (1968) . The authors implemented the GTB interface element, capable of simulating the slip and separation between two bodies in contact.

It has been known the one-dimensional beam elements have been used extensively in the analysis of composite beams with partial interaction ( Salari and Spacone, 2001 ; Dall’ Asta and Zona, 2004a; Dall’ Asta and Zona, 2004b). Considering the possibility of generic cross sections and the physical non-linearity of the cross-section’s materials, Silva (2006) , and Sousa Jr and Silva (2007), developed an interface element capable of simulating vertical separation and slip at the interface of contact between the materials of the composite section. These interface element works in conjunction with the one-dimensional beam elements implemented considering the Euler-Bernoulli beam theory.

Sousa Jr and Silva (2009), present a family of interface elements for numerical analysis of composite beams with a longitudinal deformable connection. In the model idealized by the authors, the developed formulations were based on the Euler-Bernoulli beam theory, Timoshenko beam theory, different amounts of nodes on the elements, and different integration processes of the displacements along the element. This way, the authors were able to identify, for example, which elements avoid slip locking problems in the numerical analysis of composite beams considering their partial interaction.

Sousa Jr and Silva (2010), used the interface element associated with the one-dimensional beam element for numerical analysis of multilayered beams formed by the association of beams of different materials and cross-sections, linked by deformable connections which allow for relative horizontal slip between the layers.

Silva (2010) , implemented an interface element capable of simulating deformable connections between the distinctive materials and associating a plate element and a one-dimensional beam element with degrees of freedom compatible with the plate element. For the physical and geometric non-linear analysis of the concrete slab, the author uses a formulation similar to that of the Huang et al. (1991) . Different from the cases of composite beams, the interface element presents relative longitudinal, transversal and vertical displacements for this particular problem.

The interface element implemented in this project works together with a flat plate/shell elements, so the formulation of the flat plate/shell element is also presented in the paper. There are large numbers of publications describing the characteristics of structural elements forming shells. Pereira (2005) defines a shell element as a body where one of the dimensions is much smaller than the other two. This fact allows it to reduce a three-dimensional problem into a two-dimensional problem, and thus, the displacement of any point inside the shell can be expressed in terms of the components of the mean surface displacement. This surface will be neutral, that is, it does not present a membrane behavior, in the case of homogeneous shells with only perpendicular load to the plane tangent to the middle surface.

For the development of the equations using the shell elements, two theories can be used. Kirchoff's classic theory does not take into account the shear deformation, and Reissner/Mindlin's plate theory, which according to Hughes Thomas (1987) , in addition to including transverse shear deformation, has paved the way for various interpolation schemes of nodal displacements, since, in this case, translations and rotations are interpolated independently. The results using the Reissner/Mindlin plate theory have been shown to be superior to the elements obtained according to Kirchoff's classic plate theory in the numerical analysis of thick shells, some papers present formulations of this element ( Bathe, 1982 ; Dvorkin and Bathe, 1984 ; Bathe and Dvorkin, 2005 ). As for the analysis of thin shells, the shear deformation tends to be very small and may lead to problems in the evaluation in the elements based on the Reissner/Mindlin theory, a so called shear locking effect. To avoid this, a slightly higher discretization is recommended when having thin shells, or the use of shell elements based on the Kirchoff's theory, such as the elements implemented by Batoz and Tahar (1982) and Batoz et al. (2000) .

2 NUMERICAL FORMULATION

As a practical application for the elements implemented in this work, a composite floor formed by a concrete slab connected to two steel beams is presented in Figure 1 . As can be seen in the figure, the composite floor is discretized by shell and interface element. The shell elements simulate the behavior of the concrete slab and steel beams while the interface element connects the shell elements and simulates a possible deformable connection in the contact between the concrete slabs and the steel beams.

Figure 1
Composite beam discretized by plane of the shells and interface elements

2.1 Flat plate/shell element

In this section, the formulation for the flat plate/shell finite element with nine nodes and five degrees of freedom per node at local level is presented, as shown in Figure 2 . For the development of the formulation, the consideration was applied to the Reissner/ Mindlin theory which defines the Eqs. (1) to (3) for the displacements. This theory considers deformation energy due to the shear stresses, so the rotations of the sections are different of the derived of w with respect to x and y. In the case of thin plates when these derivatives are almost equal to the rotations, numerical errors can occur due to compatibility problems between the functions that interpolate the different displacements, this effect is known in the literature as shear locking.

Figure 2
Nine node shell elements

In Eqs. (1) to (3) , $u^0$ , $v^0$ and $w^0$ represent the translocation displacements of the mean surface (or any reference surface) of the flat plate/shell element in the x, y and z directions. ${\theta }_x$ and ${\theta }_y$ represent the rotations of the sections with respect to axes x and y, while z is the position of the fiber in relation to the average surface along the thickness of the element where it is desired to evaluate the displacements. To facilitate the notation, subscript zero will be omitted in the following equations.

$$u\left(x,y,z\right)=u^0\left(x,y\right)+z{\theta }_y\left(x,y\right)$$ (1)

$$v\left(x,y,z\right)=v^0\left(x,y\right)-z{\theta }_x\left(x,y\right)$$ (2)

$$w\left(x,y,z\right)=w^0\left(x,y\right)$$ (3)

For determining the deformation equations of the flat plate/shell element, the displacement relation of Green-Lagrange was used considering the Von Karman hypothesis. Thus the derivatives of u and v in relation to x, y and z are rather small, which can be neglected, as well as the variation of w with z. Applying a compatible virtual deformation field to the deformable flat plate/shell element, we have the Eq. (4) for the internal virtual work.

$$\begin{array}{l}\delta W_{\mathrm{int}}={\displaystyle \underset{V}{\iiint }[\left(\delta u_{,x}+z\delta {\theta }_{y,x}\right){\sigma }_x+\left(\delta v_{,y}-z\delta {\theta }_{x,y}\right){\sigma }_y-\left(\delta u_{,y}+z\delta {\theta }_{y,y}+\delta v_{,x}-z\delta {\theta }_{x,x}\right){\tau }_{x,y}}\\ +\left(\delta {\theta }_y+\delta w_{,x}\right){\tau }_{xz}+\left(\delta w_{,y}-\delta {\theta }_x\right){\tau }_{yz}]dV\end{array}$$ (4)

Defining the efforts N_x, N_y , M_y, M_x as shown in Eqs. (5) to (8) and $N_{xy}={\displaystyle \underset{h}{\int }{\tau }_{xy}}dz$ , $N_{xz}=k{\displaystyle \underset{h}{\int }{\tau }_{xz}}dz$ , $Q_{yz}={\displaystyle \underset{h}{\int }{\tau }_{yz}}dz$ , and $M_{xy}={\displaystyle \underset{h}{\int }z{\tau }_{xy}}dz$ , and knowing that the displacements are constant along the thickness of the flat plate/shell element, we get, from Eq. (4) , the Eq. (9) for the internal virtual work of evaluated element.

$$N_x={\displaystyle \underset{h}{\int }{\sigma }_x^c}dz+{\displaystyle \sum _{i=1}^{nx}\frac{A_{xi}}{S_{xi}}\left({\sigma }_x^s\left(z_i\right)-{\sigma }_x^c\left(z_i\right)\right)}$$ (5)

$$N_y={\displaystyle \underset{h}{\int }{\sigma }_y^c}dz+{\displaystyle \sum _{i=1}^{ny}\frac{A_{yi}}{S_{yi}}\left({\sigma }_y^s\left(z_i\right)-{\sigma }_y^c\left(z_i\right)\right)}$$ (6)

$$M_y={\displaystyle \underset{h}{\int }z{\sigma }_x^c}dz+{\displaystyle \sum _{i=1}^{nx}\frac{A_{xi}}{S_{xi}}{z}_i\left({\sigma }_x^s\left(z_i\right)-{\sigma }_x^c\left(z_i\right)\right)}$$ (7)

$$M_x={\displaystyle \underset{h}{\int }z{\sigma }_y^c}dz+{\displaystyle \sum _{i=1}^{ny}\frac{A_{yi}}{S_{yi}}{z}_i\left({\sigma }_y^s\left(z_i\right)-{\sigma }_y^c\left(z_i\right)\right)}$$ (8)

For the flat plate/shell element to cover the reinforced concrete slabs, it is necessary that in the definition of the internal efforts, Eqs. (5) to (8) , consider the contribution to the reinforcement bars represented in these equations by their summation, where nx is the number of layers of bars in the section in the direction of x; S_xi the spacing of the bars along the x direction; A_xi, the area of bar in direction x; ${\sigma }_x^s\left(z_i\right)$ and ${\sigma }_x^c\left(z_i\right)$ the stress within the steel and concrete in the center of the layer of bars arranged in the direction of x. The same is true for the y direction.

$$\begin{array}{l}\delta W_{\mathrm{int}}={\displaystyle \underset{A}{\iint }(\delta u_{,x}{N}_x+\delta {\theta }_{y,x}{M}_y+\delta v_{,y}{N}_y-\delta {\theta }_{x,y}{M}_x+\delta u_{,y}{N}_{xy}+\delta {\theta }_{y,y}{M}_{xy}+}\\ +\delta v_{,x}{N}_{xy}-\delta {\theta }_{x,x}{M}_{xy}+\delta {\theta }_y{Q}_{xz}+\delta w_{,x}{Q}_{xz}-\delta {\theta }_x{Q}_{yz}+\delta w_{,x}{Q}_{yz})dA\end{array}$$ (9)

Considering the finite element approximation based on displacements, the equations of the displacements are approximated by shape functions associated to the nodal displacements of each element, as shown in Eq. (10) . Being that rotations and translations are interpolated independently, shape functions are given by quadric polynomials in relation to parametric coordinates 𝜉 and η, and these functions are represented by the column vector Φ of nine terms.

$$\left\{\begin{array}{l}u\\ v\\ w\\ {\theta }_x\\ {\theta }_y\end{array}\right\}=\left[\begin{array}{ccccc}{\Phi }^T& O^T& O^T& O^T& O^T\\ & {\Phi }^T& O^T& O^T& O^T\\ & & {\Phi }^T& O^T& O^T\\ & sim& & {\Phi }^T& O^T\\ & & & & {\Phi }^T\end{array}\right]q$$ (10)

In Eq. (10), q is a column vector with its terms given by forty-five nodal displacements and O is a null column vector with nine elements. Being the displacements u, v, w, θ_x and θ_y functions of nodal displacements, their variational can be written from the expression $\delta a=\delta q^T\left(\frac{\partial a}{\partial q}\right)$ , where $a$ it must be replaced by u, v, w, θ_x and θ_y. Substituting these variational into Eq. (9) and equating internal virtual work with external virtual working given by δW _ext=δq^T f_ext, we get the expression f_int -f_ext=0, where f_int is the vector of internal forces given by Eq. (11) .

$$f_{int}=\left[\begin{array}{c}{N}_x{\Phi }_{,x}+N_{xy}{\Phi }_{,y}\\ N_{xy}{\Phi }_{,x}+N_y{\Phi }_{,y}\\ Q_{xz}{\Phi }_{,x}+Q_{yz}{\Phi }_{,y}\\ -M_{xy}{\Phi }_{,x}-M_x{\Phi }_{,y}-Q_{yx}\Phi \\ M_y{\Phi }_{,x}+M_{xy}{\Phi }_{,y}+Q_{xz}\Phi \end{array}\right]$$ (11)

Using the Newton-Rapshon method in the solution of the problem f _int-f_ext=0 it is necessary to determine the derivative of this expression in relation to the nodal displacements, thus obtaining the tangent stiffness matrix. Being f_ext constant in relation to nodal displacements the tangent stiffness is given by $K=\frac{\partial f_{int}}{\partial q}$ .

$$K={\displaystyle \underset{A}{\iint }\left[\begin{array}{l}{\Phi }_{,x}{\left(\frac{\partial N_x}{\partial q}\right)}^T+{\Phi }_{,y}{\left(\frac{\partial N_{xy}}{\partial q}\right)}^T\\ {\Phi }_{,x}{\left(\frac{\partial N_{xy}}{\partial q}\right)}^T+{\Phi }_{,y}{\left(\frac{\partial N_y}{\partial q}\right)}^T\\ {\Phi }_{,x}{\left(\frac{\partial Q_{xz}}{\partial q}\right)}^T+{\Phi }_{,y}{\left(\frac{\partial Q_{yz}}{\partial q}\right)}^T\\ -{\Phi }_{,x}{\left(\frac{\partial M_{xy}}{\partial q}\right)}^T-{\Phi }_{,y}{\left(\frac{\partial M_x}{\partial q}\right)}^T-\Phi {\left(\frac{\partial Q_{yz}}{\partial q}\right)}^T\\ {\Phi }_{,x}{\left(\frac{\partial M_y}{\partial q}\right)}^T+M_{xy}{\Phi }_{,y}{\left(\frac{\partial M_{xy}}{\partial q}\right)}^T+\Phi {\left(\frac{\partial Q_{xz}}{\partial q}\right)}^T\end{array}\right]}dA$$ (12)

To determine the internal efforts and their derivatives that arise in the vector of internal forces and in the stiffness matrix, the reference surface given by the average surface of the flat plate/shell element was considered. The axial internal efforts (N_x ) and its derivative are shown in Eqs. (13) and (14) , respectively. All other internal efforts and their derivatives can be obtained in an analogous way.

$$N_x=q^T\left(\begin{array}{c}{\Phi }_{,x}\left(\lambda +2G\right)h\\ {\Phi }_{,x}\lambda h\\ O\\ O\\ O\end{array}\right)+{\displaystyle \sum _{i=1}^{nx}\frac{A_{xi}}{S_{xi}}\left(E_s-E_c\right)\left(\begin{array}{c}{\Phi }_{,x}\\ O\\ O\\ O\\ z_i{\Phi }_{,x}\end{array}\right)}$$ (13)

$$\frac{\partial N_x}{\partial q}=\left(\begin{array}{c}{\Phi }_{,x}\left(\lambda +2G\right)h\\ {\Phi }_{,x}\lambda h\\ O\\ O\\ O\end{array}\right)+{\displaystyle \sum _{i=1}^{nx}\frac{A_{xi}}{S_{xi}}\left(E_s-E_c\right)\left(\begin{array}{c}{\Phi }_{,x}\\ O\\ O\\ O\\ z_i{\Phi }_{,x}\end{array}\right)}$$ (14)

2.2 Interface Element

In this section, the formulation of the finite interface element capable of simulating the deformable connection between the flat plate/shell elements is presented. The flat plate/shell elements can be associated in three different ways, as shown in Figure 3 .

Figure 3
The types of flat shell elements in the association.

Like the flat plate/shell element, the implemented interface element locally displays 5 degrees of freedom per node. However, the possibility of simulation of structures formed by non- coplanar planes, which requires that of the sixth degree of freedom to be inserted to assemble the global stiffness matrix of the problem. In the Figure 4 , the sequence of the degrees of freedom for the formation of the nodal displacement vector is shown.

Figure 4
Degrees of Freedom to interface element

At the contact interface between the plane elements, longitudinal ( $s_l$ ), transversal ( $s_t$ ), and vertical relative displacement ( $s_v$ ), along the directions x, z and y can occur, as shown in Figure 4 . The interface element physically represents only the contact between these elements. The equations of the relative displacements in this contact, given by Eqs. (15) - (17) , are obtained from the displacement equations for the flat plate/shell elements below (α=1) and above (α=2) the element interface.

$$s_l\left(x\right)=u_2^0\left(x\right)-u_1^0\left(x\right)+\raisebox{1ex}{$h_2$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{\theta }_{z2}\left(x\right)+\raisebox{1ex}{$h_1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{\theta }_{z1}\left(x\right)$$ (15)

$$s_t\left(x\right)=w_1^0\left(x\right)-w_2^0\left(x\right)+\raisebox{1ex}{$h_2$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{\theta }_{x2}\left(x\right)+\raisebox{1ex}{$h_1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{\theta }_{x1}\left(x\right)$$ (16)

$$s_v\left(x\right)=v_2^0\left(x\right)-v_1^0\left(x\right)+z\left[{\theta }_{x1}\left(x\right)-{\theta }_{x1}\left(x\right)\right]$$ (17)

In the Eqs. (15) - (17), $u_{\alpha }^0$ , $v_{\alpha }^0$ , and $w_{\alpha }^0$ represent the translation displacements of the mean surface (or any reference surface) on the flat plate/shell elements, whereas ${\theta }_{x\alpha }$ and ${\theta }_{z\alpha }$ represent the rotations with respect to the x and z axes, which are independent of the position along the thickness, which is the hypothesis for maintaining the plane section. h₁ and h ₂ refer to the position of the contact interface and depend on the type of interface element.

Being $E_{S_b}$ , $E_{V_b}$ and $E_{N_b}$ the deformation stiffness of interface element in the direction of the relative displacements $s_l$ , $s_t$ , e $s_v$ , then the forces per unit length that arise at this interface are $S_b=E_{S_b}{s}_l$ , $V_b=E_{V_b}{s}_t$ and $N_b=E_{N_b}{s}_v$ . Applying a compatible virtual deformation field to an interface element of Figure 4 , we arrive at Eq. (18) for the virtual works of the internal forces.

$$\delta W_{\mathrm{int}}={\displaystyle \underset{A}{\int }\left(S_b\delta s_l+V_b\delta s_t+N_b\delta s_v\right)dA}$$ (18)

In Eq. (18), A is the contact area between the flat plate/shell elements, $s_l$ and $s_t$ vary only along the axis of the interface element, while $s_v$ also varies along the z axis. In the following formulation, the z axis will be considered to vary from b₁ to b₂. Figure 5 illustrates these values for the type 3 interface element. These Values are important when vertical separation is allowed in the contact interface. In these cases, they prevent the interpenetration between the elements above and below the contact interface along its width, in other words, in the z direction shown in the Figure 5 .

Figure 5
Width of the interface element verifying vertical separation

Considering $N_b^1={\displaystyle {\int }_b{N}_{b}dz}$ and $N_b^2={\displaystyle {\int }_b{N}_{b}zdz}$ to be the vertical forces along with the width (z-axis) at the contact interface and substituting the variational of the relative displacements in Eq. (18) of the virtual work we have

$$\begin{array}{l}\delta W_{\mathrm{int}}=[S_b(\delta u_2-\delta u_1)+V_b(\delta w_1-\delta w_2)+N_b^1(\delta v_2-\delta v_1)+\raisebox{1ex}{$h_1$}\!\left/ \!\raisebox{-1ex}{$s$}\right.V_b\delta {\theta }_{x1}+\\ +N_b^2\delta {\theta }_{x1}+\raisebox{1ex}{$h_1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.S_b\delta {\theta }_{z1}+\left(\raisebox{1ex}{$h_2$}\!\left/ \!\raisebox{-1ex}{$2$}\right.V_b-N_b^2\right)\delta {\theta }_{x2}+\raisebox{1ex}{$h_2$}\!\left/ \!\raisebox{-1ex}{$2$}\right.S_b\delta {\theta }_{z2}]dx\end{array}$$ (19)

In the finite element approach based on displacements, the displacement equations are approximated by shape functions associated with the nodal displacements. For the element in question, we have adopted shape functions of quadratic polynomials represented by the column vector Φ of three terms. As can be seen in equations (15) to (17) the relative displacements are obtained from a combination of the axial displacements and rotations of the flat plate/shell elements connected thereto. Equal interpolating functions are used for axial displacements and rotation, which allows a good representation of problems whose relative displacements are very small. Other combinations of interpolating functions, such as linear and quadratic for axial displacement and rotation, can generate numerical errors when the stiffness of the connection is too great, this effect is known as slip locking in the literature.

From the interpolation functions, it is possible to determine the approximate equation of the displacements associated to the nodal displacements q, and from these we can determine their variational, thus:

$$\begin{array}{l}\delta W_{int}=\delta q^T{\displaystyle \underset{L}{\int }[S_b\left(\frac{\partial u_2}{\partial q}-\frac{\partial u_1}{\partial q}\right)+N_b^1}\left(\frac{\partial v_2}{\partial q}-\frac{\partial v_1}{\partial q}\right)+V_b\left(\frac{\partial w_1}{\partial q}-\frac{\partial w_2}{\partial q}\right)+\\ +\left(\frac{h_1}{2}{V}_b+N_b^2\right)\frac{\partial {\theta }_{x1}}{\partial q}+\frac{h_1}{2}{S}_b\frac{\partial {\theta }_{z1}}{\partial q}+\left(\frac{h_2}{2}{V}_b-N_b^2\right)\frac{\partial {\theta }_{x2}}{\partial q}+\frac{h_2}{2}{S}_b\frac{\partial {\theta }_{z2}}{\partial q}]dx\end{array}$$ (20)

Analogously to the flat plate/shell element, we arrive at the internal force vector and the tangent stiffness matrix to the one-dimensional interface element.

$$f_{int=}\left\{\begin{array}{c}-S_b\Phi \\ -N_b^1\Phi \\ V_b\Phi \\ \left(\frac{h_1}{2}{V}_b+N_b^2\right)\Phi \\ \frac{h_1}{2}{S}_b\Phi \\ S_b\Phi \\ N_b^1\Phi \\ -V_b\Phi \\ \left(\frac{h_2}{2}{V}_b+N_b^2\right)\Phi \\ \frac{h_2}{2}{S}_b\Phi \end{array}\right\},K=\left[\begin{array}{c}-\Phi {\left(\frac{\partial S_b}{\partial q}\right)}^T\\ -\Phi {\left(\frac{\partial N_b^1}{\partial q}\right)}^T\\ \Phi {\left(\frac{\partial V_b}{\partial q}\right)}^T\\ \Phi \left[\frac{h_1}{2}{\left(\frac{\partial V_b}{\partial q}\right)}^T+{\left(\frac{\partial N_b^2}{\partial q}\right)}^T\right]\\ \frac{h_1}{2}\Phi {\left(\frac{\partial S_b}{\partial q}\right)}^T\\ \Phi {\left(\frac{\partial S_b}{\partial q}\right)}^T\\ \Phi {\left(\frac{\partial N_b^1}{\partial q}\right)}^T\\ -\Phi {\left(\frac{\partial V_b}{\partial q}\right)}^T\\ \Phi \left[\frac{h_2}{2}{\left(\frac{\partial V_b}{\partial q}\right)}^T-{\left(\frac{\partial N_b^2}{\partial q}\right)}^T\right]\\ \frac{h_2}{2}\Phi {\left(\frac{\partial S_b}{\partial q}\right)}^T\end{array}\right]$$ (21)

In Eq. (22) shows the derivatives of the forces in relation to the nodal displacements. In this expression, O is a column vector with three null terms. For the calculation of the integrals that appear in these expressions, a bilinear curve for the relation vertical force per unit length ( N_b) versus the vertical relative displacement (s_v ) will be considered, as shown in Figure 6 . To avoid the inter-penetrability between the flat plate/shell elements connected by the interface element, in other words, s_v negative, a very high stiffness is considered for the negative part of the bi-linear curve.

Figure 6
The relationship of the vertical force per unit length versus vertical relative displacement

$$\frac{\partial S_b}{\partial q}=\left\{\begin{array}{c}-\Phi E{S}_b\\ O\\ O\\ O\\ \frac{h_1}{2}\Phi E{S}_b\\ \Phi E{S}_b\\ O\\ O\\ O\\ \frac{h_1}{2}\Phi E{S}_b\end{array}\right\},\frac{\partial V_b}{\partial q}=\left\{\begin{array}{c}O\\ O\\ \Phi E{V}_b\\ \frac{h_1}{2}\Phi E{V}_b\\ O\\ O\\ O\\ -\Phi E{V}_b\\ \frac{h_2}{2}\Phi E{V}_b\\ O\end{array}\right\},\frac{\partial N_b^1}{\partial q}=\left\{\begin{array}{c}O\\ -\Phi {\displaystyle \underset{-b1}{\overset{b2}{\int }}{E}_{N_b}dz}\\ O\\ \Phi {\displaystyle \underset{-b1}{\overset{b2}{\int }}{E}_{N_b}zdz}\\ O\\ O\\ \Phi {\displaystyle \underset{-b1}{\overset{b2}{\int }}{E}_{N_b}dz}\\ O\\ -\Phi {\displaystyle \underset{-b1}{\overset{b2}{\int }}{E}_{N_b}zdz}\\ O\end{array}\right\}$$ _e $$\frac{\partial N_b^2}{\partial q}=\left\{\begin{array}{c}O\\ -\Phi {\displaystyle \underset{-b1}{\overset{b2}{\int }}{E}_{N_b}zdz}\\ O\\ \Phi {\displaystyle \underset{-b1}{\overset{b2}{\int }}{E}_{N_b}{z}^{2}dz}\\ O\\ O\\ \Phi {\displaystyle \underset{-b1}{\overset{b2}{\int }}{E}_{N_b}zdz}\\ O\\ -\Phi {\displaystyle \underset{-b1}{\overset{b2}{\int }}{E}_{N_b}{z}^{2}dz}\\ O\end{array}\right\}$$ (22)

3 EXAMPLES

3.1 Plate with full interaction

In this example, a rectangular plate of 3.5 by 4.8m with 10cm thick is considered. Two cases will be assessed in terms of support for this slab and three different types of loading, as presented in Table 1 below.

Table 1
Types of support conditions and loads

	Support Conditions		Type of Load
Case	Size of 3.5m	Size of 4.8m
1	supported	supported	1
2	supported	supported	2
3	supported	supported	3
4	fixed	supported	2
5	fixed	supported	3

In Table 1 , the first load is given by a concentrated load of 150kN geometrically centered on the plate. The second load is given by a linear load 35kN/m distributed parallel to the side of 3.5m in the Middle plate. The third load is given by a load of 10kN/m² uniformly distributed across the plate.

In numerical analysis, the plate is simulated in two ways. At first, only rectangular flat plate/shell element of nine nodes is used with a discretization of 4 by 4 totaling 16 elements. In the second one, the plate has its thickness divided in half and the two halves are simulated by the rectangular flat plate/shell elements of nine nodes, being these two halves connected by interface elements of type 3. That is, the 16 flat plate/shell elements with a thickness of 5cm of the top plate are connected to 16 flat plate/shell elements with a thickness of 5 cm of the bottom plate through 32 interface elements (4 by 4 rows linear elements), as shown in Figure 7 .

Figure 7
Discretization in shell planes and interface elements: (a) Omitting the upper shell elements, (b) not omitting

In the analysis using the flat plate/shell and interface elements is considered a very high connection stiffness represented by the interface element. Thus, in this case, it is expected that the response to this analysis is very close to the answer when using only the 10cm thick plate element. Numerical responses are compared with the analytical solutions obtained considering the classical theory of plates ( Timoshenko and Woinowsky-Krieger, 1959 ). In both analyses, analytical and numerical, plate material with 25GPa elasticity module and Poisson coefficient of 0.20 was considered. The following are the analytical solutions for different cases of support conditions and loading of the plate shown in Table 1 . The Eqs. (23) and (24) present the analytical solution to Case 1. Being D the plate bending stiffness, P the concentrated load applied, while a and b are the dimensions of the plate and the variables x and y are the point coordinates on the plate.

$$w\left(x,y\right)=\frac{4P}{{\pi }^{4}Dab}{\displaystyle \sum _{m=\mathrm{1,3,5}}^{\infty }{\displaystyle \sum _{n=\mathrm{1,3,5}}^{\infty }\frac{{\left(-1\right)}^{\frac{m+n}{2}-1}}{mn{\left(\frac{m^2}{a^2}+\frac{n^2}{b^2}\right)}^2}\sin }}\frac{m\pi x}{a}\sin \frac{n\pi \left(y+0.5b\right)}{b}$$ (23)

$$w\left(\frac{a}{2},y\right)=\frac{4P}{{\pi }^{4}Dab}{\displaystyle \sum _{m=\mathrm{1,3,5}}^{\infty }{\displaystyle \sum _{n=\mathrm{1,3,5}}^{\infty }\frac{1}{mn{\left(\frac{m^2}{a^2}+\frac{n^2}{b^2}\right)}^2}}}$$ (24)

For the other cases and loads the middle plane deformed of the plate is given by the Eq. (25) where the parameters K, A_m, B _m and Cm are presented in Table 2 .

Table 2
Parameters for calculation of the deformed plate

Case	K	$A_m$	$B_m$	$C_m$
2	$q_l$	$-\frac{\sin \frac{m\pi }{2}(2+{\alpha }_m\tanh {\alpha }_m)}{2m^4\cosh {\alpha }_m}$	$\frac{\sin \frac{m\pi }{2}}{2m^4\cosh {\alpha }_m}$	$\frac{\sin \frac{m\pi }{2}}{m^4}$
3	$\frac{2q_0}{\pi }$	$-\frac{(2+{\alpha }_m\tanh {\alpha }_m)}{2m^5\cosh {\alpha }_m}$	$\frac{1}{2m^5\cosh {\alpha }_m}$	$\frac{1}{m^5}$
4	$q_l$	$-\frac{2\sin \frac{m\pi }{2}\cosh {\alpha }_m(\tanh {\alpha }_m+{\alpha }_m)}{m^4(\sinh 2{\alpha }_m+2{\alpha }_m)}$	$\frac{2\sin \frac{m\pi }{2}\sinh {\alpha }_m}{m^4(\sinh 2{\alpha }_m+2{\alpha }_m)}$	$\frac{\sin \frac{m\pi }{2}}{m^4}$
5	$\frac{2q_0}{\pi }$	$-\frac{2\cosh {\alpha }_m(\tanh {\alpha }_m+{\alpha }_m)}{m^5(\sinh 2{\alpha }_m+2{\alpha }_m)}$	$\frac{2\sinh {\alpha }_m}{m^5(\sinh 2{\alpha }_m+2{\alpha }_m)}$	$\frac{1}{m^5}$

$$w\left(x,y\right)=\frac{2a{K}^3}{{\pi }^{4}D}{\displaystyle \sum _{n=\mathrm{1,3,5}}^{\infty }\left[A_m\cosh \left(\frac{m\pi y}{a}\right)+B_m\sinh \left(\frac{m\pi y}{a}\right)+C_m\right]}\sin \frac{m\pi x}{a}$$ (25)

For a rectangular plate with a = 3.5m and b = 4.5m and considering four significant digits in the calculation of the displacement, the values given in Table 3 are for the maximum displacement on the plate considering the different support conditions and loadings.

Table 3
Maximum displacement on plate

Case	1	2	3	4	5
Maximum displacement	$0.1994\frac{P}{D}$	$0.4803\frac{q_l}{D}$	$1.0343\frac{q_0}{D}$	$0.3128\frac{q_l}{D}$	$0.6614\frac{q_0}{D}$

In Table 3 , P, q_l and q₀ are obtained from the different types of loading considered, that is P = 150kN, q_l= 35kN/m, q₀ = 10kN/m². While D is already the plate bending stiffness and is calculated using the material mechanical properties and the thickness of plate; doing this, D = 2170.139kNm.

Table 4 shows the values of the maximum displacement on the plate obtained from numerical and analytical analyses. In this Table, analysis 1 is the numerical analysis using only the flat plate/shell elements. Analysis 2 is the numerical analysis using flat plate/shell and interface elements, considering the total interaction of the connection. Analytical integration implies that the integrals of the area necessary for the calculation of the flat plate/shell element stiffness matrix are calculated analytically, while the reduced integration implies that these integrals are calculated using numerical integration using a reduced number of integration points to the terms relating to shear.

Table 4
Maximum displacement on plate (cm)

Case	Analysis 1		Analysis 2		Analytical
	int. analytical	int. reduced	int. analytical	int. reduced
1	1.13	1.33	0.93	1.10	1.38
2	0.78	0.85	0.62	0.71	0.78
3	0.47	0.49	0.36	0.40	0.48
4	0.54	0.63	0.52	0.61	0.51
5	0.27	0.31	0.26	0.30	0.30

As per Table 4 , we observed that the results obtained for the different analyses are near to one another, as was expected, since for analysis of the interface elements was considered a total interaction on the connection, namely, rigidity connection between the two halves of the plate. The difference found is due to an inevitably different discretization for the two analyses and also to the fact that for analysis 1, an average surface of the plate does not suffer membrane effect just the bending effect, while in analysis 2 the average surfaces of the halves plate suffering membrane and bending effects.

Reduced integration tends to prevent shear locking when the plate thickness is too small. This locking occurs because the interpolation functions are the same for rotations and translations, such as the shear deformation is given by the difference of the rotation with the derivative of the translation, whereby the non-compatibility of the interpolating functions of these two shall generate spurious values in the evaluation of the deformation, when this tends to have small values. In Table 4 there is a greater influence of this effect in case 1. This occurred because in this case, the shear has greater influence in the analysis, soon errors in the determination are more evident in the responses.

In Table 4 , we came to realize that the numerical solutions are quite close to the analytical solutions. There is a slight difference that can be explained due to the difference of theories used in the implementation for the flat plate/shell element and the classic theory of plates. While one ignores the deformation energy from the shearing force, the other considers this energy. It is exactly this difference in the theories that permits verification that Case 1 has a greater shearing influence, since, as shown in Table 4 , Case 1 has a greater difference in the numerical analysis solution.

3.2 Plate with free interaction

In this example cases 1, 2 and 3 of the previous example will be considered, but now the plate will be formed by the overlap of two plates of thickness of 5cm, with different plate bending stiffness and free interaction. The upper and lower plate bending stiffness is D _s and D_i, respectively. The contact between the two plates is considered complete in the vertical direction and free in the transverse direction; in other words, the slip in the contact is free while there is no uplifted vertical separation between the plates, therefore, they will have even vertical displacement and curvature.

The analytical solution for these three cases that considers two overlapping plates with free slip between them is determined in this work by utilizing the technique described below.

The loading applied to the upper plate is absorbed by deformation of that plate, and part is transmitted to the lower plate and absorbed by its deformation. To define the loading portion of the upper plate that is transmitted to the lower plate, the condition of equal displacement in the geometric center of both plates will be used.

The transfer of the load from the top plate to the bottom plate occurs more intensely in the center of the plate and remains very small in the contour thereof. Thus, it is assumed that the loading is concentrated in the center of the plate. Therefore, for Case 1, the concentrated force P_i transmitted to the lower plate is obtained from Eq. (26) .

$$\frac{4\left(P-P_i\right)}{{\pi }^4{D}_{s}ab}=\frac{4P_i}{{\pi }^4{D}_{s}ab}\to P_i=\frac{D_i}{D_s+D_i}P$$ (26)

Following the same reasoning for Eq. (26), we arrive at Table 5 for the maximum displacement in the plate for the other two cases.

Table 5
Maximum displacement for free interaction plate

Case	1	2	3
Analysis 2	$0.1994\frac{P}{D_s+D_i}$	$0.4803\frac{q_l}{D_s+D_i}$	$1.0343\frac{q_0}{D_s+D_i}$

In Table 5 , P, q_l and q₀ are obtained from the different types of loading considered, where, P =150kN, q_l= 35kN/m, and q₀ = 10kN/m². Assuming for the upper plate, E= 25000MPa and ν= 0.2, and for the lower plate E= 70000MPa and ν= 0.28, we have: D_s= 271.27kNm and D_i = 791.20kNm.

Table 6 below shows the maximum displacement values ​​in the plate with free interaction obtained from the numerical and analytical analyzes. This Table is analogous to Table 4 , so for some additional information, consult the descriptive text in Table 4 .

Table 6
Maximum displacement on the plate with free interaction (cm)

Case	Analysis 2		Analytical
	int. analítica	int. reduzida
1	2.28	2.63	2.81
2	1.57	1.74	1.58
3	0.96	1.01	0.97

From Table 6 it can be seen that the results of the maximum displacement observed using the elements implemented in this project and the analytical solutions are very close to one another. As in the previous example, the small difference can be explained due to the difference of the theories used in the implementation of the flat plate/shell element and the classical plate theory. Also, as in the previous example, the difference in the theories regarding whether to consider or not the energy of deformation due to the shear, shows that Case 1 suffers greater influence of the shear, since it is verified in Table 6 that the greater difference of the analytical solution with the numerical is given by Case 1.

3.3 Plate with Vertical Separation

In this example, the same plate from the previous example formed by the overlapping of two boards of 5cm thickness is analyzed. As in the previous example, the contact between the two plates is free in the longitudinal and transversal directions. However, in this example, the vertical separation between the two plates is allowed, while inter-penetrability must be prevented. This situation is obtained by assuming a bi-linear curve for the relationship between the relative displacements and the force per unit of length arising in that direction. A low stiffness value is allowed (E_NB=10^-3KPa ) in the case of positive relative displacement; in other words, separation between the plates, and a high stiffness (E_NB=10⁹KPa ) for negative relative displacements, which represents inter-penetration between the plates.

A distributed load of 10kN/m² is applied across the lower plate. As the loading is applied to the underside of the rectangular plate and null stiffness is admitted in the direction of vertical separation between the plates, then only the lower plate is expected to deform. Table 7 presents the results of the displacements considering numerical analysis using the flat palte/shell elements and the interface element, and the analytical solution. In this case the analytical solution is obtained admitting only the bottom plate working.

Table 7
Maximum displacement on the plate with vertical separation (cm)

ENb (KPa) sv < 0	ENb (KPa) sv > 0	EVb (KPa)	ESb (KPa)	Displacement Lower	Displacement Up
109	10-3	10-3	10-3	1.28	0.009
	Analytical			δ = 1.31

It can be seen from Table 7 that the analytical solution considering only the lower plate working is practically equal to the displacement obtained from the numerical analysis for the lower plate. The displacement on the top plate is very small, showing that only the bottom plate is working. Figure 8 presents the deformation of the plate considering the free vertical interaction. For the sake of better visibility of the upper and lower plates deformation, and knowing that the analyzed plate has two axes of symmetry, only a quarter of the plate is shown in the figure. As expected, it is observed from Figure 8 that the top plate does not deform.

Figure 8
Deformed plate for free vertical interaction

3.4. Composite Beam

Sousa Jr and Silva (2007) numerically simulate a steel-concrete composite beam with partial interaction using beam elements for the concrete section, as well as the steel profile and an interface element to simulate the deformable connection and associate the beam elements. The composite beam analyzed has a cross section formed by a rectangular section of concrete measuring 550 by 120mm and a steel profile of type I, with height of 300mm, width of the lower and upper flange of 150mm, thickness of the web of 10mm, and thickness of the flanges of 16mm. The mechanical properties of steel and concrete are given by: E_c = 30500MPa, ν_c= 0.2, E_s = 20000MPa, and ν_s= 0.3. The composite beam is bi-supported with a span of 6m subjected to a uniformly distributed load of 80kN/m.

The same composite beam of Sousa Jr and Silva (2007) is analyzed in this example by using the elements implemented in this work. As shown in Figure 9 , the concrete slab is discretized in 4 plane shell elements of nine nodes, called here SHELL9. The steel profile flange and web are also discretized by 4 SHELL9 elements. The connection between the plane shell elements of the concrete slab and the upper flange of the steel profile is made by 4 interface elements, called here INTSHELL6. It is assumed that there is total interaction in the vertical and transverse direction, whereas in the longitudinal direction, a stiffness of 9482kN/m² is considered. This stiffness is associated to the parameter αL=1 commonly used in literature to identify the degree of connection of composite beams with partial interaction. For further details of this parameter, see Sousa Jr and Silva (2007).

Figure 9
Composite discrete beam in plane of shell elements and interface

The graph of Figure 10 below shows the variation of the displacement along the composite beam. In this graph, 'reference' refers to the result obtained by Sousa Jr and Silva (2007) and shell9+intshell6 to the result obtained in this work. As it can be seen from Figure 10 , the results are practically coincident.

Figure 10
Variation of the displacement (cm) along the span (m) of the composite beam

4 CONCLUSION

The efficiency of the flat plate/shell and interface finite element formulations developed and implemented for the numerical simulation of surface structures with deformable connection were duly proven with results obtained in numerical and analytical examples using the classical plate theory.

• http://dx.doi.org/10.1590/1679-78254640

References

• Bathe, K. J. (1982). Finite element procedures in engineering analysis. New Jersey, EUA: Prentice-Hall.
• Bathe, K. J., Dvorkin, E. N., (2005). A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation. International Journal for Numerical Methods in Engineering 21: 367-383.
• Batoz, J. L, Hammad, F., Zheng, C., Zhong, W., (2000). On the linear analysis of plates and shells using a new-16 degrees of freedom flat shell element. Computers and Structures 78: 11–20.
• Batoz, J. L., Tahar, M. B., (1982). Evaluation of a new quadrilateral thin plate bending element. International Journal for Numerical Methods in Engineering 18: 1655-1677.
• Dall’Asta, A., Zona, A., (2004a). Comparison and validation of displacement and mixed elements for the non-linear analysis of composite beams. Engineering Computations 82: 2117–2130.
• Dall’Asta, A., Zona, A., (2004b). Three-field mixed formulation for the non-linear analysis of composite beams with deformable shear connection. Finite Elements in Analysis and Design 40: 425–448.
• Dvorkin, E. N., Bathe, K. J., (1984). A continuum mechanics based four‐node shell element for general nonlinear analysis. Engineering Computations 1: 77-88.
• Goodman RE., Taylor RL, Brekke TL, (1968). A model for the mechanics of jointed rock. Journal of the Soil Mechanics and Foundations Division 99: 637-659.
• Huang Z, Burgess IW, Plank RJ, (1991). Nonlinear Analysis of Reinforced Concrete Slabs Subjected to Fire. ACI Structural Journal 96: 127-135.
• Hughes Thomas JR, (1987). The Finite Element Method. Englewood Cliffs, New Jersey, EUA: Prentice-Hall.
• Kaliakin VN, Li J., (1995). Insight into deficiencies associated with commonly used zero thickness interface elements. Computers and Geotechnics 17: 225-252.
• Pereira, J. C. (2005). Curso de Estabilidade Estrutural. Florianópolis, Santa Catarina, Brasil. (in Portuguese).
• Salari, M. R., Spacone, E., (2001). Finite element formulations of one-dimensional elements with bond-slip. Journal of Structural Engineering 23: 815–826.
• Silva, A. R., (2006). Análise Numérica de Vigas Mistas com Interação Parcial, M. Sc. Dissertation (in Portuguese) Federal University of Ouro Preto, MG, Brazil.
• Silva, A. R., (2010). Análise Numérica de Elementos Estruturais com Interação Parcial. D. Sc. Thesis (in Portuguese), Federal University of Ouro Preto, Brazil.
• Sousa Jr. J. B. M. and Silva, A. R., (2007) Nonlinear analysis of partially connected composite beams using interface elements. Finite Elements in Analysis and Design 43: 954-964.
• Sousa Jr., J. B. M., Silva, A. R., (2009). A family of interface elements for the analysis of composite beams with interlayer slip. Finite Elements in Analysis and Design 45: 305-314.
• Sousa Jr, J. B. M., Silva, A. R., (2010). Analytical and numerical analysis of multilayered beams with interlayer slip. Engineering Structures 32: 1671-1680.
• Timoshenko, S., Woinowsky-Krieger S., (1959). Theory of plates and Shells. New York, EUA: McGraw-Hill.

Publication Dates

• Publication in this collection
  2018

History

• Received
  02 Nov 2017
• Reviewed
  06 Feb 2018
• Accepted
  09 Feb 2018