Numerical formulation for nonlinear analysis of concrete and steel shells with deformable connections

A two-dimensional interface finite element capable of associating flat shell elements positioned one above the other was developed. The implemented interface element can physically simulate the contact between the flat shell elements and connect the reference planes of the shell elements above and below it. The formulation presented allows consideration of nonlinear behavior for the deformable connection as well as for the concrete and steel materials that make up the shell structure. One of the practical applications analyzed in this research is the numerical simulation of composite floors formed by a reinforced concrete slab connected to steel beams through a de-formable connection. In this case, the concrete slab and the steel beams are discretized by flat shell elements and the deformable connection is discretized by two-dimensional interface elements. Experimental and numerical results from literature were used to validate the implemented elements. In the two examples analyzed, the results obtained for the displacements were close, with the difference, in the first case, being associated with uncertainties during the experimental test and in the second, the difference in theories used in the formulation of the implemented elements.


Introduction
A formulation of a two-dimensional interface finite element that can numerically simulate a deformable connection between structural elements and be analyzed by flat shell elements is presented.Practical civil construction components, such as composite floors, beams, and pillars, have a deformable connection between different materials.Among these, composite floors formed by a reinforced concrete slab connected to steel beams by shear connectors are the most common.In general, the structural analysis and design for this element are based on simplified models that can generate significant errors (e.g., by using bar elements to represent the elements).Therefore, the effect of the variation of the normal stress along the width of the concrete slab or shear lag is not considered, and if torsion effects on the beam are significant, it is necessary to accurately represent the warping.One way to better represent the mechanical behavior of this structural element and an optimized design is to model the concrete slab and steel beams by flat shell elements and to model the deformable connection by interface elements.
The interface element was initially developed to work in conjunction with two-and three-dimensional elements that represent a thin layer of material or the contact between two distinct materials, such as the case of soil-structure interaction.The first study on interface elements was by Goodman et al. (1968).In that study, the interface element was used to simulate the slip and separation between two bodies in contact.
To represent two distinct materials, Kalikin and Li, (1995) used twodimensional interface elements with zero thickness to model the problem of contact between the soil and a shallow foundation.Subsequently, Carol et al. (2001) used an interface element of zero thickness to analyze the cracking process in concrete elements.
For the analysis of the ultimate capacity of a composite beam with a deformable connection, Sousa and Silva (2007) developed an interface element capable of simulating longitudinal slip and vertical uplift at the contact between materials of a composite section.The implemented interface element works in conjunction with one-dimensional beam elements implemented by considering Euler-Bernoulli beam theory, the physical nonlinearity of materials, and the possibility of generic cross sections.
Sousa and Silva (2009) presented a family of interface elements for the numerical analysis of composite beams with longitudinal deformable connections.The proposal of these elements includes formulations to be employed with Euler-Bernoulli beam theory, Tymoshenko's beam theory, different numbers of nodes in the elements, and different integration processes for the displacements along the element.
Sousa and Silva (2010) presented an analytical solution for composite beams with multiple layers.This solution was used to verify the ability of the interface element to simulate composite beam problems with more than one shear plane.
The interface element implemented in the present study is associated with flat shell elements above and below it.According to Batoz et al. (2010), the first analysis of a shell structure using the finite element method was performed using a set of flat shell elements to approximate the true shell shape.Owing to the simplicity of the formulation.computational efficiency, and application flexibility, shell elements are extensively used in practice.
According to Hughes (1987), Reissner-Mindlin's plate theory, which includes transverse shear strain, has promoted the development of several interpolation schemes of nodal displacements, because, in this case, translations and rotations are interpolated independently.Because of this, shell elements have recently been obtained based on Reissner-Mindlin's theory, and these are superior to the elements obtained according to Kirchhoff's classic plate theory in the numerical analyses of thick shells.In the case of thin shells, shear strain tends to be very small and may lead to problems in the evaluation of these elements in the analysis using Reissner-Mindlin's theory (the shear locking effect).To avoid this, a slightly higher discretization is recommended when thin shells are obtained, as well as reduced numerical integration of the shear strain.
One of the practical applications analyzed in this study is the numerical simulation of composite floors formed by a reinforced concrete slab connected to steel beams through a deformable connection.Despite the prodigious use of composite floors in construction, few numerical and experimental studies on composite floors can be found in literature.A vast majority of studies focus on simplification of the problem, modeling the floor as a composite beam.One of the most relevant articles for the study of composite floors is that of Nie et al. (2008), which introduces a new definition of the effective width for the final resistance calculations of composite beams subjected to a bending moment using a commonly accepted rectangular stress distribution.Izzuddin et al. (2004) developed a shell element to simulate concrete slabs reinforced by cold-rolled steel sheets.The formulation, as in the present article, is based on Reissner-Mindlin's plate theory, but with a modification that makes the proposed element capable of simulating orthotropic geometry and the discontinuity of the material between the adjacent ribs.
In this study, a two-dimensional interface finite element capable of associating flat shell elements positioned one above the other was developed.The implemented interface element can physically simulate the contact between the flat shell elements and connect the reference planes of the shell elements above and below it or side by side.The formulation presented allows consideration of nonlinear behavior for the deformable connection as well as for the concrete and steel materials that make up the shell structure.
One of the practical applications analyzed herein is the numerical simulation of composite floors formed by a reinforced concrete slab connected to steel beams through a deformable connection.In this case, the concrete slab and the steel beams are discretized by flat shell elements and the deformable connection is discretized by two-dimensional interface elements.
Experimental and numerical results from the literature were used to validate the implemented elements.Two examples were analyzed, and their results demonstrated the efficiency of the implemented element.

Flat shell element
Figure 1 shows a composite floor formed by a concrete slab connected to two steel beams.In the discretization of the composite floor, a two-dimensional interface element is used to make the connection between the flat shell elements and to simulate a possible deformable connection in the contact between the concrete slab and the steel beams.The flat shell finite element implemented for the nonlinear analysis of a steel and reinforced concrete shell has nine nodes and five degrees of freedom per node, as shown in Figure 2. Because of the independence between translation and rotational degrees of freedom, the formulation described in this section accounts for shear deformation; therefore, it is ap-plicable to thick shells.In the case of thin shells, attention should be given to possible numerical errors resulting from shear locking.To avoid this, we can refine the finite element mesh and Ouro Preto, 74(3), 297-307, jul Using the kinematic assumptions of Reissner-Mindlin's plate theory, we obtain the displacement equations given by: where u 0 , v 0 , and w 0 represent the translations of the reference plane of the flat shell element in the x, y, and z directions, respectively; θ x and θ y are the rotations of the sections in relation to the x and axes, respectively; and z is the position of the fiber along the thickness of the flat shell element.To simplify the notation, the zero superscript is omitted in the following equations.
In the definition of the strain-displacement equations, the Green-La-grange relations were used by considering Von Karman's hypothesis, which implies that the derivatives of u and v in relation to y, and z are small and can be neglected, while the variation of w can be also be neglected: Geometric nonlinearity is observed in Eqs. ( 4) and ( 5): Another nonlinearity of the problem is defined by the stress-strain relationship of steel and concrete for a uniaxial stress state.For concrete under compression, the curve defined by the CEB-FIP model code (2010) was used.For the concrete under tension, the stress-strain curve of Figure 3 was adopted, as suggested by Rots et al. (1984) and also used by Huang et al. (1999Huang et al. ( , 2003a)).use reduced integration in the evaluation of shear strain.
In the analysis of physical nonlinearity, a multilayer element is used and specific characteristics are considered for each layer.These include, for example, different mechanical properties and independent stress-strain relation-ships.For other characteristics and more details of the flat shell element, consult Huang et al. (1999Huang et al. ( , 2003aHuang et al. ( , 2003b) ) and Silva and Dias (2018a).
For the particular case of isotropic material, observed when the principal strain is outside the failure region, we have and the D xy matrix of Eq. ( 10) reduces to the traditional form of Reissner-Mindlin's plate theory.
In the definition of the internal force vector and tangent stiffness matrix, it is assumed that the displacements and rotations have quadratic variations along the element and can be written in terms of the nodal displacements.For details on how to obtain the numerical formulation of this element, consult Silva and Dias (2018a).
where E 1 and E 2 are given by tangents to the material stress-strain curve at points e = e α and G α = 0.5E α / (1+v) for α =1, 2. The stiffness matrix in the direction of the x-y axes (D xy ) can be obtained from D 12 and the angle φ of rotation of the principal axes in relation to the x-y axes: where Herein, e tu = 10e tr was adopted, and, for the concrete tension strength, the relationship f t = 0.3321 √ f c (ASCE, 1982) was used, where f c is the concrete compressive strength in MPa.For the steel, the stressstrain curve defined by linear segments was adopted, where values of the limits of stresses and strain are presented in the application examples.
For the physical and geometric nonlinear analysis, an incremental method with displacement control was used.At each step of the incremental method, a linear material with a modulus of elasticity given by the tangent to the strain-strain curve was considered.In this manner, the stress-strain relationship can be obtained by using Hooke's law for the problem analyzed.Following the constitutive matrix for materials under stress levels inside and outside, the failure region was defined.
The materials exhibited orthotropic behavior inside the failure region (after cracking or crushing of concrete and after yielding of the steel); that is, they exhibited different characteristics for each principal direction.By considering the layers in the plane stress state, the principal directions were calculated, as indicated herein by subscripts 1 and 2, where 1 indicates the direction of the greater principal strain.If the principal strains (e 1 and e 2 ) are within the failure region, the materials are considered orthotropic with the stress-strain relationship decoupled to the principal directions; in this way, the constitutive matrix of the material is given by: In Eqs. ( 20)-( 22), u 0 a , v 0 a , and w 0 a represent the translational displacements of the midplane (or any reference plane) of the flat shell elements; θ xa and θ ya represent the rotations with respect to the x and y axes, which are independent of the position along the thickness, the hypothesis of the maintenance of the plane section; and h 1 and h 2 are the thicknesses of the flat shell elements below and above, respectively, the interface element.In the following equations, the index 0 will be omitted to simplify the notation.
Because E Sb , E Vb , and E Nb are the deformable connection stiffnesses of the interface element in the direction of the relative displacements s l , s t , and s v , the forces per unit area that emerge at this interface are S b = E Sb s l , V b = E Vb s t , and N b = E Nb s v .Applying a compatible virtual deformation field to an interface element of Figure 4, we have, through the principle of virtual work, In the finite element approximation based on displacements, the displacement equations are approximated by the shape functions associated with the nodal displacements (q).Because the lower and upper degrees of freedom are independently interpolated and the order of these followed the same order as that for the flat shell element shown in Section 2.1, the shape functions for the two-dimensional interface element are the same as those adopted for the flat shell element.
Analogous to the flat shell element, we arrive at the internal force vector and stiffness matrix of the two-dimensional interface element given by Eq. ( 25).For the element in question, shape functions given by quadratic polynomials represented by the column vector Φ of nine terms were adopted: Substituting the variation of the relative displacements in Eq. ( 23) of the principle of virtual work, we have s l (x, y) = u 0 2 (x, y) -u 0 1 (x, y) -h 2 /2 θ y2 (x, y) -h 1 /2 θ y1 (x, y) s t (x, y) = v 0 1 (x, y) -v 0 2 (x, y) -h 2 /2 θ x2 (x, y) + h 1 /2 θ x1 (x, y) (24) The two-dimensional interface element has, in addition to the function of physically simulating the contact between the flat shell elements, the function of connecting the reference surfaces of the flat shell elements above and below it.For this, reference surfaces are defined, which in this case are taken as the mid-plane of the element to be discretized into flat shell elements.In this way, the thickness or distance between the upper and lower nodes of the two-dimensional element is taken as the sum of half of the upper and lower thickness of the elements connected by the two-dimensional element.Although this thickness is not zero, its actual physical thickness is always zero.
The interface element physically The 18-node two-dimensional interface element is responsible for simulating the deformable connection and making the connection between 9-node flat shell elements, as shown in Figure 4.The orientation of the degrees of freedom follows analogously to that of the flat shell element of Section 2.1.

Two-dimensional interface element
S b , V b , and N b can be obtained directly from the force per unit area versus the relative displacement curve.The derivatives of these force terms in relation to nodal displacements are given by where O is a column vector with nine null terms.Nie et al. (2008) numerically and experimentally evaluated a floor formed by a reinforced concrete slab connected by shear connectors to five steel beams.In their model, the load was applied through three hydraulic jacks in increments of 2 kN until rupture.During the test, displacements, strain, and slip were monitored.Figure 5 shows the composite floor with the geometric parameters of the slab, steel beams, reinforcing bars, and shear connector.Load P is applied to the three longitudinal beams and divided into four points on each beam.In this work, the slab and the I-shaped steel profile were simulated by a nine-node flat shell element, and the deformable connection was represented by a two-dimensional interface element of 18 nodes, as shown in Figure 6.

Steel-concrete composite floor
Among the results provided by Nie et al. (2008), the numerical and experimental load-displacement curves of the composite floor are shown.In the analysis using the elements implemented in this study, f c = 30.3MPa, e c2 = 0.2%, e cu = 0.35%, and v = 0.17 for the concrete.For the reinforcement bars and steel beam, a stress-strain curve defined by elastic-perfectly plastic-hardening tri-linear model was used.The four points that define the tensioned part of the curve are: (0;0), (0.1432%; 295 MPa), (0.23%; 295 MPa) e (2%; 448 MPa) (steel beams), and (0;0), (0.1845%; 380 MPa), (0.2%; 380 MPa) e (4.5%; 478 MPa) (reinforcement bars).For both steels E s = 206000 MPa, v = 0.30, and compression behavior equal in traction was admitted.For the longitudinal and transverse stiffness of the deformable connection between the steel beams and the concrete slab, the equation of Ollgaard et al. (1971) or Aribert (1992) was used: S b = S bu (1 -e C1Sl ) C2 , with C 1 = 0.7 mm -1 and C 2 = 0.4.The ultimate strength of the connector can be obtained through the equations described in the technical norms and textbooks on the subject.For stud connectors Figure 9 shows the discretization of the composite beam in flat shell elements and a two-dimensional interface.From the figure, it is observed that the concrete slab was discretized by 20 flat shell elements.The steel beam had four divisions in the longitudinal direction, and, for each division, there were two flat shell elements for the flange and one for the web.For the two-dimensional interface elements, four elements were used.
It is observed that the result obtained by using the implemented elements is relatively close to that obtained from the experimental test, demonstrating the effectiveness of the implemented elements.
In this example, the elastic line of a fixed composite beam with a 10 m span was subjected to a concentrated force of 50 kN at the midspan is analyzed.The cross section of the composite beam was composed of a concrete slab, with E c = 13 GPa and v = 0.2, connected to an I-shaped steel profile (IPE 300), with E s = 200 GPa and v = 0.30. Figure 8 shows the cross section with the dimensions.   in the contact area).For the vertical stiffness at the contact interface, a high stiffness (E Nb = 10 9 kPa ) was considered, that is, the total interaction for the vertical uplift.

Elastic line
Figure 7 illustrates the results obtained from the numerical analyses using the elements developed in the present study and results of the experimental model of Nie et al. (2008).
The influence of longitudinal stiffness variation on vertical displacement along the length of the beam was evaluated.The values of longitudinal stiffness used were E Sb = 10 3 kPa and E Sb = 10 7 kPa Because transverse and vertical displacements at the contact interface were not allowed, values for the stiffnesses of E Vb = 10 7 kPa and E Nb = 10 7 kPa respectively, were adopted.
In the following figures, Shell9+int6 is the result for the composite beam analyzed by using the nine-node flat shell element and the six-node interface element implemented by Silva and Dias (2018b).Shell9+Int18 is the result obtained by simulating the composite beam through the 9-node shell element and the 18-node interface element implemented in this study.
BeamTQ+IntTQ is the result obtained by simulating the steel beam with the three-node beam element (BeamTQ) and the six-node interface element (IntTQ); both are elements presented in Sousa and Silva (2007), wherein the beam theory of Timoshenko was used in the formulation and quadratic interpolation for translation and rotation.REF is the result obtained by Brighenti and Bottoli (2014), which simulates the beam with beam elements based on Euler-Bernoulli beam theory.
The responses for the three analyses (Shell9+Int18, Shell9+Int6, and BeamTQ+IntTQ) are shown in Figure 10, and they are practically the same.The difference observed when compared with the response indicated by REF in Figure 10 is due to the difference between the theories used in the simulation of the problem, because Shell9+Int6, Shell9+Int18, and BeamTQ+IntTQ account for the shear strain, which is not considered in the REF analysis.The effects of shear lag (normal stress variation along the slab width) and poisson, which are verified when simulating the problem using flat shell elements and do not appear in the beam analysis, do not cause this difference because the responses are practically identical for the Timoshenko beam analysis (BeamTQ+IntTQ) and for the Reisner-Mindlin's plate analysis (Shell9+Int6 or Shell9+Int18).
Figure 11 shows the variation along the width of the concrete slab of normal stress obtained in the most compressed fiber.It can be seen from the figure that the results obtained by Shell9+Int6 and Shell9+Int18 analyses are very close to each other.The response obtained from the BeamTQ+IntTQ analysis does not exhibit normal stress variation along the width because it simulates the problem by considering beam theory with flexion in only one plane.Figure 11 shows that, although the plate analysis reveals a shear lag effect, the areas bound by the curves and the horizontal axis are close together, generating the same contribution for the concrete slab in the beam and plate analyses.
Similar to Figure 10, it can be seen from Figure 12 that both the 6-and 18-node interface elements give almost identical results, but they present a different relationship from that of the BeamTQ+IntTQ analyses.In this case,    Figure 13 shows that the area bound by the curves and the horizontal axis for the plate analysis is larger than the area found using Timoshenko beam analysis (for the same load).Therefore, greater stiffness is found when the problem is simulated through the beam theory of Timoshenko, which generated a smaller displacement, as shown in Figure 12.This effect, which also manifests when using Euler-Bernoulli beam theory, generated a smaller displacement in the REF analysis to compensate for the difference in the response between the flat shell elements and the Euler-Bernoulli beam element shown in Figure 10.That is, a reduction of the concrete slab width should be made so that the displacement and the area defined by the region bound by the beam theory curve and the horizontal axis are the same as in the plate analysis.Therefore, it is concluded that beam analysis may generate an overestimation of the bending stiffness of the composite beam.

Figure 1 -
Figure 1 -Composite floor discretized by flat shell and interface elements.

Figure 3 -
Figure 3 -Stress-strain relationship for concrete in tension.
Figure -4 Two-dimensional interface element associated with flat shell elements.

Figure 6 -
Figure 6 -Steel-concrete composite slab discretized in flat shell and interface elements.

Figure 8 -
Figure 8 -Cross section of the mixed section (in units of mm).

Figure 9 -
Figure 9 -Discretization of the fixed composite beam.

Figure 10 -
Figure 10 -Elastic line of composite fixed-fixed beam with E Sb = 10 3 kPa.

Figure 11 -
Figure 11 -Stress variation in the most compressed fiber along the concrete slab width for E Sb = 10 3 kPa.

Figure 12 -
Figure 12 -Elastic line of fixed composite beam with E Sb = 10 7 kPa.

Figure 13 -
Figure 13 -Stress variation in the most compressed fiber along the concrete slab width for E Sb = 10 7 kPa.

Formulations
of the flat shell and twodimensional interface finite elements are presented for the numerical analysis of shell structures with possible deformable connections.The most practical case would be the composite floors formed by reinforced concrete slabs associated with steel beams by means of a deformable connection.Two examples were analyzed to verify the potential usefulness of the elements implemented in this study.Some other examples have been evaluated by other researchers using numerical methods, commercial software, and experimental tests.The results demonstrate that the implemented elements are a useful tool for structural analysis in some practical situations.The authors would like to thank the Federal University of Ouro Preto, PROPEC, CNPq, and FAPEMIG for their financial support.