Numerical simulation of composite steel-concrete alveolar beams: web-post buckling, Vierendeel and flexural mechanisms

Composite alveolar beams consist in the union of two structural systems largely employed in civil construction sector: the steel-concrete composite beams and the alveolar steel beams. Thus, its use allows their advantages to be enhanced, enabling to design even larger spans and to achieve more economical and sustainable solutions. Considering that Brazilian and international standards do not directly specify criteria for the analysis and design of these beams, in this paper it is presented the development and validation of an updated finite element model, using ANSYS software, capable of simulating different failure modes that may occur, such as web-post buckling, Vierendeel mechanism and flexural mechanism. The obtained results presented a good correlation with experimental results from previous works. After the model validation, the effect of the openings on the composite beam was investigated and discussed, and it was concluded that the web-post buckling may limit the structural gains on load capacity, so it is important to adopt opening patterns that enhance the resistance of the beam to this mode of failure.


INTRODUCTION
Steel-concrete composite beams represent an efficient solution to explore the best mechanical properties of these materials, especially in simply supported beams, where the concrete slab works basically under compression and the steel profile under tension. The joint work of these elements is ensured by shear connectors, which aim to restrict longitudinal slip and the vertical separation at the interface. This structural system provides an increase in the strength of the composite beam, when compared to the structural elements working independently, enabling the design of larger spans and the achievement of more economical solutions (Nie et al., 2004).
Alternatively, the employment of alveolar steel beams, which consist of steel profiles with sequential openings in the web, can be an advantageous solution to overcome larger spans, especially when subjected to distributed loads with a low or medium magnitude (Cimadevila et al., 2000). Its manufacturing process consists in cutting the original steel profile longitudinally in a certain pattern, resulting in two parts that can be repositioned and then welded together in a new configuration, in which the flanges are farther apart, as shown in Figure 1. Thus, with practically the same self-weight, the expanded profiles are produced with a greater moment of inertia and, consequently, greater flexural strength, resulting in a better performance under serviceability limit states. According to Badke-Neto et al. (2015), depending on the shape of the openings, the alveolar beams can be named as castellated beams (with hexagonal openings) or cellular beams (with circular openings). Furthermore, the use of steel alveolar beams presents other advantages, in addition to the structural ones, such as an appealing aesthetic aspect and the possibility of passing pipes through the openings, which allows the optimization of the ceiling height in buildings. However, in terms of structural behavior, these beams also present some disadvantages, such as the reduction in shear resistance and a worse performance when subjected to concentrated loads, in comparison to full web beams, often requiring reinforcements on the web, thus increasing their cost (Badke-Neto et al., 2015). Moreover, due to their peculiar geometry, specific modes of failure may occur and must be properly understood for their correct employment as reliable structural elements. Kerdal and Nethercot (1984) studied and listed these modes of collapse: (i) Formation of a Vierendeel mechanism; (ii) Buckling of web-post due to shear; (iii) Rupture of a welded joint in a web-post; (iv) Lateral-torsional buckling of an entire span; (v) Formation of a flexure mechanism (plastic hinge); and (vi) Buckling of web-post due to compression.
Aiming to enhance the advantages and to mitigate the disadvantages of these two structural systems, the composite alveolar beams are formed by composite steel-concrete beams with alveolar steel profiles. Therefore, its use allows the design of even larger spans, the achievement of even more economical solutions, the reduction in material consumption and, consequently, the reduction of the generated environmental impacts. In addition, it is expected that the composite action may improve the resistance of these beams to concentrated loads and to some failure modes, in comparison to alveolar steel beams working independently. In this context, Redwood (2000) outlined that once the composite action increases the resistances to flexural and Vierendeel mechanisms, there is an increased likelihood of web-post buckling in composite castellated sections. Thus, the resistance of these beams to failure modes involving local instabilities is a fundamental issue to be investigated.
The Brazilian standard NBR 8800 (2008), the Eurocodes EN 1993-1-1 (2005 and EN 1994EN -1-1 (2004, and the American standard ANSI/AISC 360-16 (2016) do not directly specify criteria for the analysis and design of composite alveolar beams considering their specific modes of failure. In the scientific literature there are some important formulations available, such as the works of Veríssimo et al. (2006) and Lawson and Hicks (2011), applicable to steel and composite beams with web openings, and, more recently, the AISC Steel Design Guide 31: Castellated and Cellular Beam Design (Fares et al., 2016). However, none of these publications have been recognized as an official standard. Therefore, advances in numerical simulation are essential for a better understanding of the structural behavior of composite alveolar beams and can be helpful for the validation and improvement of future standardizations. In the last decades, several researchers have developed finite element models for the analysis of composite beams (Gattesco, 1999;Queiroz et al., 2007;Marconcin et al., 2010;Tamayo et al., 2015;Dias et al., 2015;Liu et al., 2017;Katwal et al., 2018), alveolar steel beams (Ellobody, 2012;Durif et al., 2013;Erdal and Saka, 2013;Vieira, 2015;Wang et al., 2016;Oliveira et al., 2019;Shamass and Guarracino, 2020) and composite alveolar beams (Megharief, 1997;Müller et al., 2006;Gizejowski and Salah, 2008;Bake, 2010;Ferrari, 2013;Lawson et al., 2018). Since the approached numerical problem presents many nonlinearities, both physical (materials behavior) and geometric (profile buckling), different modeling strategies have been used, employing several commercial and authorial software products, with varying degrees of simplification.
In this context, the aims of this work are: • To develop and validate an updated finite element model for composite alveolar beams, using ANSYS software, version 19.2, with current-technology elements and compatible material models. This model should be able to simulate different modes of failure, especially the web-post buckling and the yielding mechanisms.
• To investigate and discuss the effect of the openings on the behavior of composite alveolar beams, by comparing the results of the validated model to the numerical results of similar composite beams without holes.

NUMERICAL MODEL
The following items present information about the element types, material models, boundary conditions, initial geometrical imperfections and solution stages of the developed finite element model.

Element types
The steel profile was modeled by 4-node quadrilateral shell elements, with six degrees of freedom per node (translations and rotations in x, y, z directions), named SHELL181 in ANSYS (2018). Both membrane and bending stiffnesses were considered. The element formulation is based on the work of several authors, including Bathe and Dvorkin (1986), and it uses the Reissner-Mindlin first-order shear-deformation theory, being suitable for analyzing thin to moderately thick shell structures. It is applicable to both linear and nonlinear problems, including those with large deformations. Thus, its formulation employs logarithmic strain and true stress measures rather than nominal engineering strain and stress, but for small deformations the difference between nominal and true values is negligible. The element can also contain several layers. In this work, a single and centralized layer was used, with five integration points along the thickness.
The concrete slab was modeled by 20-node hexahedral solid elements, with three degrees of freedom per node (translations in x, y, z directions), named SOLID186 in ANSYS (2018), whose formulation is based on Zienkiewicz et al. (2013). Since SOLID186 is classified as a current technology element, it is compatible with several current ANSYS features, such as the generation of embedded elements and the use of updated material models. In the present work, these elements were employed in their homogeneous form with full integration.
The shear connectors were modeled by nonlinear spring elements, named COMBIN39 in ANSYS (2018), which act in the longitudinal direction of the composite beam, relating the ux displacements of the nodes they are applied to. This relation respects the curve of shear force versus longitudinal slip, which can be experimentally obtained from push-out tests. These elements were applied to nodes at the interface, which are spaced by half the thickness of the top flange (tfs/2), because of the centralized positioning of the SHELL181 elements' cross-sections, as shown in Figure 2. Since this distance is very small, these nodes can be assumed as practically coincident. At the same time, regarding the transverse and vertical directions, the displacements uy and uz of these same nodes were perfectly coupled.
The reinforcement steel bars of the concrete slab were modeled by discrete embedded elements, named REINF264 in ANSYS (2018). These elements use the same nodes of the base elements SOLID186, even if their geometric position does not coincide with them, and present only axial stiffness. Therefore, the stiffnesses to bending, torsion and shear are neglected. A perfect interaction between the reinforcing elements (REINF264) and the concrete base element (SOLID186) is admitted, so there is no relative movement between them.
A new ANSYS functionality was used for the generation of the reinforcement embedded elements, named mesh-independent method. This method uses MESH200 elements, which do not directly contribute to the solution but determine the positions where REINF264 reinforcement elements are created, being classified as guide-only elements. Thus, it is possible to insert the positions of the reinforcement bars from the lines drawn in absolute coordinates, differently from the standard method, in which it is necessary to use relative coordinates in respect to the base elements, resulting in mesh dependence during the process of reinforcement generation. The steel-deck sheet was modeled by 8-node shell elements, with six degrees of freedom per node (translations and rotations in x, y, z directions), named SHELL281 in ANSYS (2018). A perfect interaction between the concrete slab and the steel-deck sheet was assumed, so, for this purpose, the lower slab faces were selected for the mesh generation of the shell elements. Therefore, once the SOLID186 elements have 8 nodes on each face, it was necessary to choose this 8-node shell element for the steel-deck sheet, instead of the 4-node shell element that was used for the steel profile.
In Figure 3, the mentioned elements are highlighted in the developed model to a generic composite cellular beam.

Figure 3
Finite elements used in the developed model.

Material models
The profile steel was modeled by von Mises yield criterion, with isotropic hardening, adopting the constitutive model proposed by Gattesco (1999), illustrated in Figure 4a. This model is divided into three stages of loading: (i) elastic-linear; (ii) yield plateau; (iii) hardening governed by parabolic curve, given by equation (1): where fy and fu are the steel yield and ultimate strengths, εy = fy /E is the yield strain, εh is the strain at the start of the hardening process, εu is the strain at ultimate stress, E is the steel modulus of elasticity and Eh is the tangent modulus of elasticity at the start of the hardening process. The steels of the reinforcement bars and the steel deck sheet were simplified as perfectly elastoplastic (Figure 4b). The connectors' behavior was represented by the real constants of the COMBIN39 nonlinear spring element, in which the points of the shear force versus slip curve, obtained experimentally via push-out tests, can be stored. In order to systematize the model, it was adopted the theoretical curve given by equation (2), proposed by Ollgaard et al. (1971), that is adjustable to the push-out test data of stud bolts.
where Q is the shear force acting in the connector; Qu is the ultimate shear force resisted by the connector; s is the slip; e = 2.718 is the Euler number; and m, n are curve fitting parameters. Two different models were used to simulate the concrete's behavior: 1. a customized model using the USERMAT subroutine; 2. an ANSYS native model, named DP-CONCRETE, recently made available by the software (since 17.0 version).
Once these two models were recently developed, in this work it was decided to use both and to compare the obtained results, aiming to increase the reliability of the numerical model.
The USERMAT is an ANSYS subroutine that can be programmed to customize a material model (Quevedo et al., 2018). The customized model adopted was developed by Lazzari et al. (2019), based on Ottosen four-parameter criterion (Ottosen and Ristinmaa, 2005). The constitutive relation used for compressive behavior is suggested by fib Model Code 2010 (FIB, 2013), as given in the equation (3), and is illustrated in Figure 5a. On the other hand, the constitutive relation for tensile behavior considers the tension stiffening effect and is illustrated in Figure 5b.
where fcm is the mean uniaxial compressive strength of concrete; σc is the compressive stress; εc is the compressive strain; η= εc/εc1; εc1 is the strain at the peak compressive stress; εc,lim is the ultimate compressive strain; κ= Eci/Ec1 is the plasticity number; Eci and Ec1 are the tangent and the secant modulus of elasticity of concrete, respectively.
Latin American Journal of Solids and Structures, 2020, 17(5), e289 6/28 In Figure 5, besides the symbols described for eq. (3), fctm is the mean uniaxial tensile strength of concrete; α is a reducing factor; and εctu is the limit strain value in tension. The parameters α = 0.6 and εctu = 0.001 were adopted, as suggested by Schmitz (2017). In addition to these values, the USERMAT model also uses the mean biaxial compressive strength of concrete (fc2m) to define the Ottosen surface. In the absence of experimental values, the fib Model Code 2010 (FIB, 2013) may be used to calculate some of these parameters.
The DP-CONCRETE model is a native ANSYS model, recently developed, and based on the combination of two yield surfaces. Once a single Drucker-Prager surface does not represent the large differences in tensile and compressive behavior of concrete, this model uses a first Drucker-Prager yield surface for compression, and a second one, that can be a Drucker-Prager or a Rankine surface, for tension and tension-compression. In this paper, it was used the combination Drucker-Prager with Rankine, that is represented in the two-dimensional principal stress system in Figure 6. The surfaces' formulations can be found in Chen (2007). In order to define these surfaces, it is necessary to inform the values of the uniaxial and the biaxial compressive strengths of concrete (named Rc and Rb in ANSYS), and of the concrete tensile strength (named T). In this work, Rc = fcm, Rb = fc2m and T = fctm were adopted. The hardening and softening laws for tension and compression depend on the chosen HSD model and start to act in a point after its stress state reaches the respective yield surface -until then, it presents linear elastic behavior. There are four different HSD models available in ANSYS, which associate the relative stress evolution (Ω) with the effective plastic strain (κ). In this research, the Linear HSD model was adopted, as illustrated in Figure 7, in which the highlighted parameters (Ωci, Ωcr, κcm, κcr, κtr) need to be inserted by the user. Further information about DP-CONCRETE and HSD models can be found in ANSYS (2018)   In order to approximate, as much as possible, the shape of the hardening and softening laws of the two models (Usermat and DP-Concrete), the values of Ωci = 0.4, Ωcr = 0.65, Ωtr = 0.02 and κtr = 0.001 were adopted, and the parameters κcm and κcr were calculated by the equations (4) and (5), so that the total strains (elastic added to the plastic ones) in compression are equal to εc1 when the concrete reaches the maximum stress, and to εc,lim when the concrete reaches the residual stress. In their turn, the values of εc1 and εc,lim can be obtained from the fib Model Code 2010 (FIB, 2013). The value Ωci = 0.4 means that the concrete presents linear elastic behavior until the effective compressive stress reaches 40% of fcm, thus a reduced modulus of elasticity (Ec) should be used for concrete. In this work, it was adopted Ec = 0.9Eci.

Boundary conditions
The beam was admitted as simply supported in the developed model, and the symmetry condition was used when applicable. In this case, the nodes at the first support, on the bottom flange of the steel profile, had the displacements in y and z restricted; and the nodes at the central cross-section of the beam had the displacements in x and the rotations around y and z restricted, as shown in Figure 8a. Otherwise, when the beam was not symmetrical, the nodes at the first support had the displacements in x, y and z restricted and the nodes at the second support had the displacements in y and z restricted, as shown in Figure 8b. A concentrated load was preferably applied by imposing y-displacements on the respective nodes, at the top face of the slab, once, as outlined by Queiroz et al. (2007), the displacement control may overcome convergence problems. However, as explained by these authors, in cases of distributed loads, or of a set of multiple concentrated loads, it is necessary to apply forces, since in these cases it is difficult to establish a relation between loads and the associated displacements, especially during the plastic range of behavior.

Stages and methods of solution
Aiming to add initial geometric imperfections to the steel profile through the combination of buckling modes, in order to simulate the failure modes involving local instabilities, the analysis was performed in four stages: Stage 1: Solution of a linear static analysis, in which a unit load was distributed among the same nodes of the load that is applied in the final analysis (stage 4).
Stage 2: Solution of an Eigen-Buckling analysis, which consists in a linear eigenvalue problem (Bathe, 2014) to determine the buckling modes and the load factors associated with the linear static analysis performed in stage 1.
Stage 3: Insertion of geometric imperfections to the steel profile by updating the model geometry, based on the buckling modes (eigenvectors) calculated in stage 2. The weighted combination of two buckling modes was applied (usually the first ones), with an amplitude of dg/600 each, where dg is the height of the expanded profile. This amplitude value was suggested by Bake (2010).
Stage 4: Solution of the nonlinear analysis. When the load was applied by imposing equivalent displacements, the full Newton-Raphson method with displacement control was used to solve the nonlinear problem. Otherwise, when forces were applied, the Arc-Length method (Riks, 1979) was adopted, aiming to capture the post-buckling behavior, Latin American Journal of Solids and Structures, 2020, 17(5), e289 8/28 since the Newton-Raphson method with force control does not support sudden changes in the structure stiffness that lead to negative force increments. Further explanations about these methods can be found in Bathe (2014). The main advantage of the Arc-Length method is that the load increment factor is treated as an extra unknown to the system, thus both positive and negative load increments are supported during the solution.

Comparison to other numerical models
In terms of composite beams modeling, the developed model resembles the strategy used by Queiroz et al. (2007), with the connectors being modeled by nonlinear spring elements, the concrete slab by solid elements and the steel profile by shell elements. These authors also employed ANSYS software, however the element types used by them (SOLID65 and SHELL43) are classified as legacy elements in the most current versions of the software (ANSYS, 2018). Therefore, it was decided to use the current-technology elements SOLID186 and SHELL181, as in Schmitz (2017). One important issue is that the mentioned 8-node solid element (SOLID65) is the only one compatible with the material model CONCRETE, based on Willam and Warnke formulation (Chen, 2007), which had been usually employed to the concrete modeling in ANSYS. However, besides the fact that this element is considered outdated, numerical instability problems associated with this material model are described by several researchers, who reported the necessity to disable its standard crushing option, substituting it by the von Mises criterion in order to improve convergence (Queiroz et al, 2007;David, 2007;Marconcin et al., 2010). This is the main reason why, in this paper, it was decided to use more recent material models (DP-CONCRETE and USERMAT), which are compatible with the element SOLID186.
Regarding the composite alveolar beams modeling itself, in Table 1 a comparison is made between the strategies employed by different works, in terms of the material models, elements types and geometric imperfections. It can be noted that the mentioned researches used others finite element software, thus the differences for the concrete's behavior modeling were already expected, considering the different libraries of materials available. Besides, the constitutive model used in the present research for the profile steel is also different, adopting nonlinear hardening.
In terms of the model construction, the steel-deck sheet was not modeled by Müller et al. (2006) and by Ferrari (2013). Beyond, the strategy used by Bake (2010) differs from the others, since the author modeled the slab by shell elements and admitted full interaction between the steel profile and the concrete slab. In terms of geometric imperfections, differences in the amplitudes and in the number of modes are observed. In this context, it is worth mentioning that, for the analysis of full web profiles, the Eurocode EN 1993-1-5 (2006 indicates that the amplitude of the imperfection must be equal to A/200, where A is the minor dimension of the rectangular panel. However, a priori, this value is not applicable to alveolar beams, reason why, in this paper, it was decided to use the value suggested by Bake (2010). Finally, it may be emphasized that the mentioned works in Table 1 were selected because they also simulated some of the beams used for Latin American Journal of Solids and Structures, 2020, 17(5), e289 9/28 the model validation, so they can be used for comparing results -reason why it is important to observe the differences between the adopted strategies.

VALIDATION EXAMPLES
The validation of the proposed numerical model was performed through the numerical simulation of experimental tests carried out by previous works, with the subsequent comparison between the results. For this purpose, two castellated beams tested by Hosain and Speirs (1973) and five composite cellular beams tested by Müller et al. (2006) and by Nadjai et al. (2007) were selected.

Experimental tests configurations
The castellated beams A-1 and G-1 geometries are presented in Figure 9. In both beams, the load was applied at the midspan. Instead of measuring the midspan deflection, the authors chose to measure the rotation at the second vertical support. Lateral supports have also been added at the specified positions. Further details about the tests can be found in Hosain and Speirs (1973).

Figure 9
Geometries of the beams A-1 and G-1, tested by Hosain and Speirs (1973).
The geometries of the composite cellular beams A1 and B1 are presented in Figure 10. Both have a free span of 450 cm, and were produced with the same concrete, steel reinforcement (mesh A142), steel-deck plate (HOLORIB 51/150) and connectors arrangement (one row, ф19 mm, 120 mm height, at each 15 cm). They differ in the geometries of the cellular profiles and in the points of load application. The expanded profile of beam B1 is asymmetric, once it had been generated from two distinct full web profiles. Stiffeners were added at the points of load application and at the supports. Further information about the tests can be found in Nadjai et al. (2007).
The geometries of the composite cellular beams 1 and 3 are presented in Figure 11. Beam 1 was tested twice: in the first test (named 1A) the beam has failed by web-post buckling due to shear, close to the second support. When the buckling started, the beam was unloaded, and then the web was stiffened in this location, through the sealing of the openings with a rigid bar. Next, a new test was carried out (named 1B). Both beams 1 and 3 have the same concrete slab geometry, steel reinforcement (0.4% in both directions), steel-deck sheet (HOLORIB 51/150) and connectors arrangement (one row, ф19 mm, 100 mm height, at each 15 cm). They have been loaded with four concentrated loads, in order to simulate a uniformly distributed load. The expanded profile of beam 1 is symmetric (generated from two IPE 400 profiles), and of the beam 3 is highly asymmetric, once it has been generated from two distinct steel profiles (HEB340 and IPE 300). Stiffeners were added at the points of load application and at the supports. Further information about the tests can be found in Müller et al. (2006).

Material properties
Some of the material properties were provided in the original papers, such as the profile steel strengths (fy, fu) and the concrete's mean uniaxial compressive strength (fcm). These properties were inserted in the numerical models, but many others were necessary to complete the input data. Since the results of push-out tests were not provided, the values of Qu, m and n were arbitrated based on the values used by Schmitz (2017) for connectors with the same diameter and height, which were calibrated with the push-out tests of Chapman and Balakrishnan (1964). For the profile steel, the arbitrated values for εh/εy and Eh were based on Gattesco (1999) and Sadowski et al. (2015). The concrete properties were calculated according to the fib Model Code 2010 (FIB, 2013), using the provided fcm of each beam. All the material properties inserted in the numerical model are presented in Table 2. The symbols that appear in this table were described in section 2 of this article. For the profile steel properties, subscripts w (for the web) and f (for the flanges) were added. For example, fyw represents the yield strength of the web. It should be noted that, since the steel profiles of beams B1 and 3 are asymmetrical, the steel strengths are presented for the lower and upper parts, respectively.

Numerical model implementation
The number of elements used in each numerical model is presented in Table 3. The symmetry condition was adopted only in models of the beams tested by Nadjai et al. (2007), as shown in Table 4. For the generation of SHELL181 elements, first the profile areas were subdivided into regular areas, aiming to create a mapped mesh (containing only quadrilateral elements). After, the profile mesh was generated using the maximum element size criterion, adopting the value of 3 cm for the beams tested by Hosain and Speirs (1973) and Nadjai et al. (2007), 3.5 cm for beam 3 and 4 cm for beam 1, both tested by Müller et al. (2006). The connectors were arranged in two rows, instead of in a single one as in the experimental tests. This adaptation was necessary due to the modeling strategy employed (with spring elements), since when using only one row it was observed that the slab had a strong tendency to rotate around the longitudinal axis, especially in the preliminary Eigen-Buckling analysis. Therefore, the values of Qu, shown in Table 2, were divided by two, in order to compensate the duplicated number of connectors, thus the global behavior of the composite beams was not changed. To exemplify the mesh refinement adopted, Figure 12 Table 4 presents a summary of the performed analyses and indicates the solution methods adopted for the final stage. Because of the load configuration, it was not possible to use displacement control for the beams tested by Müller et al. (2006), so the Arc-Length method was preferably adopted, aiming to capture their post-buckling behavior. However, since this method presented convergence difficulties when using the Usermat model for concrete, in these cases the Newton-Raphson method with force control was adopted, thus the convergence stopped after the maximum load was reached. It is worth mentioning that, in order to simulate the web stiffening of beam 1 that was performed in test 1B by Müller et al. (2006), a rigid element (named MPC184 in ANSYS) was used in the numerical model, as shown in Figure 13, with the direct elimination method, i.e. through the addition of compatibility equations that relate the translation and rotation degrees of freedom of the nodes the element is applied to.

VALIDATION RESULTS
The results obtained with the developed numerical model for the described examples are now presented, being compared to the experimental results and, when available, to numerical results from previous works.

Hosain and Speirs (1973) beams
The numerical and experimental results for the diagram of applied load versus rotation at the support, for beams A-1 and G-1, tested by Hosain and Speirs (1973), are presented in Figure 14. According to Hosain and Speirs (1973), the beginning of the failure process of beam A-1 occurred with the formation of a Vierendeel Mechanism, which started to develop, and, after the steel hardening, culminated in a local buckling around one of the central openings. The developed numerical model captured this behavior. The von Mises equivalent stresses for the applied load of 173 kN are presented in Figure 15, where it can be verified that the corners of the central openings are completely yielded and the deformed shape resembles a parallelogram, characteristics of the Vierendeel Mechanism.
On the other hand, in beam G-1, as the web-posts and the bases of the hexagonal openings are narrower, the secondary moments do not develop significantly in the tee-sections. Therefore, instead of the Vierendeel Mechanism formation, there was a pure plastic hinge formation (flexural mechanism), due to the principal moment, as reported by Hosain and Speirs (1973). The numerical model captured this failure mode, as shown in Figure 16, where the von Mises stresses for an applied load of 164 kN are presented: it can be noted that the entire region above and below the central openings are yielded. After the total development of the flexural mechanism, the beam presented local web-post buckling in the experimental test, which was also captured by the numerical model, as shown in Figure 17

Nadjai et al. (2007) beams
The numerical and experimental results for the diagram of applied load versus midspan deflection, for beams A1 and B1, tested by Nadjai et al. (2007), are presented in Figures 18 and 19, respectively. As it can be observed, the developed numerical model was able to well simulate the experimental behavior of these beams, whether using the Usermat model or the DP-Concrete model.
In the experimental tests, both beams failed by web-post buckling due to shear, as reported by Nadjai et al. (2007). The numerical models were able to capture this mode of collapse, as illustrated in Figure 20, that shows the displacements in z direction of beam A1, when uy = 1.4 cm. The insertion of initial geometric imperfections to the model through the combination of elastic buckling modes has proved to be an appropriate strategy to capture this failure mode.
Moreover, it can be noted that both models adopted for concrete -DP-Concrete and Usermat -presented quite similar results. Figure 21 shows the normal stresses in x direction on beam A1, when uy = 1.0 cm, obtained by both models. Except for small localized differences, the stresses distribution presented very similar patterns.

Figure 20
Web-post buckling due to shear in beam A1 (transversal displacements, in cm, when uy = 1.4 cm).

Figure 21
Beam A1: Comparison between the normal stresses in x direction, in kN/cm 2 , when uy=1.0 cm.
Regarding the steel profile, it was observed that, when the web-post buckling started (uy≈0.9 cm for beam A1 and uy≈1.0 cm for beam B1), the yielding had only started at concentrated points around the openings, but it did not develop farther, as shown in Figure 22, which represents the equivalent plastic strains for beam B1. Therefore, failure modes related to yielding (Vierendeel or flexural mechanisms) were not yet about to occur.

Müller et al. (2006) beams
The numerical and experimental results for the diagram of applied load versus midspan deflection, for beam 1 in tests 1A and 1B, performed by Müller et al. (2006), are presented in Figure 23. In test 1A diagrams, it can be noted that the results of the developed models presented an expressive correlation with the experimental result. Moreover, the results of Usermat and DP-Concrete models were almost coincident until buckling occurred. However, after reaching the peak load, the convergence of the Usermat model stopped, since the Newton-Raphson method with force control was used (because in this model the Arc-Length method had already presented convergence difficulties in the elastic range of behavior). The numerical models were able to capture the web-post buckling due to shear, as shown in Figure 24, where the displacements in z direction for the 725 kN load are shown (at the post-peak load range), which were obtained by the DP-Concrete model. Figure 25 presents a numerical-experimental comparison of the y-displacements evolution over the span for different load levels, illustrating the good correlation between numerical and experimental results.
In test 1B diagrams, it can also be observed a good correlation between the numerical and experimental results, although not as good as in test 1A. However, the evaluation of these results must be done with caution, since this second test was carried out after the first one was unloaded, therefore permanent plastic strains and concrete cracks may have interfered in the final experimental response, and were not considered in the numerical analysis. Anyway, the numerical model captured the web-post buckling that occurred at the other end of the beam, as shown in Figure 26, since the web-post that failed in the first test was stiffened in the test 1B.   The numerical and experimental results for the diagram of applied load versus midspan deflection, for beam 3, tested by Müller et al. (2006), are presented in Figure 27. In the experimental test, the failure started with the formation of a Vierendeel Mechanism, and then culminated in the buckling of the last web-post, as reported by Müller et al. (2006). The two proposed numerical models, with Usermat and DP-Concrete, presented results practically coincident with each other and with the numerical results of Bake (2010). In relation to the experimental test, the models presented a good correlation during the elastic range, but slightly overestimated the loads at the plastic range of behavior, after the formation of the Vierendeel Mechanism began. These small differences could possibly be justified by not considering the self-weight, which reduces the resistance of the beam to yielding, and/or by the arbitrated hardening parameters for profile steel, which were not provided by the authors.
In the Usermat model the convergence stopped after the maximum load was reached, as expected, once the Newton-Raphson method with force control was adopted. In the DP-Concrete model, although the Arc-Length method was used, the convergence stopped shortly after the maximum load was reached, at the beginning of the decreasing part of the curve, when the solution showed the tendency to converge to an elastic discharge, thus it was unable to capture post-buckling behavior. This convergence difficulty during the post-peak load range can be justified by the complexity of the numerical model, which contains a large number of elements, and by the behavior of the analyzed beam, which involves yielding and buckling in an irregular geometry, with an asymmetrical profile and different web thicknesses.
The numerical models were able to capture the formation of the Vierendeel Mechanism, which occurred in the corner of the last opening, at the upper part of the steel profile, where the thickness is thinner. This corner is localized at the non-concreted end of the beam, where the profile works isolated, therefore it is more susceptible to this failure mode. Figure 28 shows the von Mises stresses in that part of the beam, for the applied load of 516 kN. When the peak load was reached, the buckling of the last web-post was already quite developed and visible in the numerical models. Besides influencing the yielding, the different web thicknesses also changed the shape of the web-post buckling, with the transversal displacements occurring higher up, in the thinner part of the web. This difference was reported by Müller et al. (2006) and can also be observed numerically, as shown in Figure 29.

Summary and evaluation of validation results
In Tables 5 and 6, a comparison is made between the numerical results of this paper and the experimental results, in terms of failure modes and their respective loads, for the two castellated beams and five composite cellular beams analyzed, respectively. In these tables, the following abbreviation is adopted for the failure modes: formation of Vierendeel mechanism (VM), local buckling (LB), formation of flexural mechanism (FM) and web-post buckling (WPB). Regarding the castellated beams (Table 5), the numerical models captured the correct failure modes with excellent precision in their respective failure loads, presenting an absolute mean of the relative differences in relation to the experimental results equal to 1.91%.
With respect to the composite cellular beams (Table 6), it is initially emphasized that, despite Müller et al. (2006) have reported the formation of a Vierendeel mechanism before the web-post buckling in beam 3, the authors did not specify the load when it was formed in the experiment. In the numerical models, to determine the formation load of this mechanism, it was verified the load step in which the corner of the opening was completely yielded, as shown in Figure 28. Regarding the other results, referring to the buckling of the web-posts, it can also be observed a good precision of the numerical results, with absolute means of the relative differences in relation to the experimental results equal to 4.07%, for the models with DP-Concrete; and 3.02%, for the models with Usermat. Therefore, it can be concluded that the developed numerical models were capable to satisfactorily simulate the beams of the validation examples.

EFFECT OF THE OPENINGS
To evaluate the effect of web openings on the global behavior of the composite alveolar beams, similar composite beams without openings were also numerically simulated, aiming to compare their results. Two examples analyzed in the validation stage -beam A1 by Nadjai et al. (2007) and beam 1 by Müller et al. (2006) -were chosen for this comparison, since these beams have symmetrical expanded steel profiles. Finally, it was proposed an example of a beam with a larger span (11 m), subjected to a uniformly distributed load. In each of these cases, two additional composite beams were simulated: one with the original steel profile, before being expanded, and the other with a full web profile with the same dimensions as the expanded profile, but without holes.
Since the Usermat and DP-Concrete models presented similar results in the validation stage, it was decided to carry out this study only with the DP-Concrete model. The other characteristics of the models, such as mesh refinement, element types and material models, were kept similar to those previously described.
Latin American Journal of Solids and Structures, 2020, 17(5), e289 21/28 Figure 30 shows the results of the composite cellular beam A1 in comparison to the two corresponding composite beams with the full web profiles (UB 406x140x39 original profile and an expanded profile without holes). It can be verified that in this example the stiffness gain obtained with the expansion of the original profile is not significant, possibly due to its small free span (4.5 m, as shown in Figure 10), so that shear forces are considerable in comparison to the bending moments. In addition, it is observed that there is a decrease in the ultimate load capacity, since the composite cellular beam failed by web-post buckling due to shear, which occurred at the load of 351.92 kN (numerical) and 370.12 (experimental). Therefore, from a structural point of view, in this case there are no great advantages in expanding the steel profile, although in a real project this expansion could be considered for other reasons, such as the passage of ducts and pipes or the achievement of a lighter design.

Figure 30
Load-deflection curves: comparison between composite beams and composite cellular beam A1. Figure 31 shows the numerical results of the composite cellular beam 1 compared to the two corresponding composite beams with the full web profiles considered (IPE400 original profile and an expanded profile without holes). Unlike the previous example, in this case it is already possible to notice a greater gain in initial stiffness, which occurs because the span is larger (6.84 m, as shown in Figure 11), and the four concentrated loads are distributed along the span, simulating a uniformly distributed load. For example, Table 7 presents the loads corresponding to a midspan deflection of 27.4 mm (equal to L/250), where it can be observed that the profile expansion generates a gain of load equal to 25.42% in test 1A and 26.42% in test 1B. However, it can be once more verified that, in terms of ultimate behavior, the web-post buckling limited the possible gains in ultimate load capacity.

Figure 31
Load-deflection curves: comparison between composite beams and composite cellular beams 1A and 1B.

Beam with 11m span and subjected to a uniformly distributed load
The results of the last two examples show that with a greater free span the structural gains with the steel profile expansion become more significant. Thus, a natural sequence of studies would be the analysis of a beam with an even greater span. However, once experimental tests with larger spans were not found in literature, it was necessary to create a new example to be analyzed with the numerical model developed and validated in this paper.
Following this logic, it was proposed the case of a beam with 11m span, subjected to a uniformly distributed load. The other geometric characteristics and material properties were admitted, for simplification purposes, as being the same as those of beam 1 by Müller et al. (2006). IPE400 profile was used as the original steel profile, and expansions with different types of opening patterns were analyzed. Figure 32 shows the nomenclature adopted for the geometric parameters, and Table 8 presents the names adopted for the beams and their respective geometric parameters. Composite beams with full web profiles were named FW1 and FW2. The proposed composite castellated beams cover the most usual patterns, called Litzka (CA1), Anglo-Saxon (CA2) and Peiner (CA3). The proposed composite cellular beams have different ratios s/D0, assuming the values of 1.1 (CE1), 1.2 (CE2), 1.3 (CE3) and 1.4 (CE4), all of them with a diameter of 44 cm and expansion ratio k=1.5. The beam CE5 has a smaller diameter (38 cm), s/D0=1.5, and expansion ratio k=1.388, as in the pattern used by Müller et al. (2006) in the tests of beam 1. All beams were numerically analyzed considering the symmetry condition and using the Arc-Length method with force control.   Figure 33 presents the results obtained for the load-deflection curves of all cases considered in this example. As it can be seen, all beams with expanded profiles, except for beam CE1, had significant gains in initial stiffness. Table 9 Latin American Journal of Solids and Structures, 2020, 17(5), e289 23/28 presents the applied linear load values for a midspan deflection of 44 mm (equal to L/250), showing the mentioned gain. Due to the smaller expansion ratio k = 1.388 of beam CE5, it presented lower gains than CE4, CA1, CA2 and CA3 beams. At the same time, beam CE4, with an expansion ratio k=1.5, had a greater gain than beam CE5, but slightly less than the castellated beams. It has occurred because, despite having profiles with equal height (dg = 1.5x40 = 60 cm), the heights of the tee-sections are different due to the manufacturing process in double cut (cellular beams) or single cut (castellated beams), as shown in Figure 1. On the other hand, the composite castellated beams CA1, CA2 and CA3 presented similar results, once their expanded profiles have not only the same height dg, but also the same tee-section height ht.

Figure 33
Load-deflection curves for the beam with 11m span Regarding the ultimate behavior, beams CE1, CE2 and CE3 failed early by web-post buckling due to shear, thus in these three cases the gain in load capacity was limited. This failure mode occurred in theses beams due to the slenderness of their web-posts, which increases with the decrease of the s/D0 ratio, resulting in a greater susceptibility for the occurrence of buckling. Figure 34 shows the out-of-plane displacements for these three beams, at their respective ultimate loads (32.45 kN/m, 50.01 kN/m and 61.01 kN/m). Such results indicate that the two phenomena (web-post buckling and concrete cracking) may be related. Therefore, it is recommended to carry out more in-depth and specific studies to verify the influence of the reinforcement bars in this region on the resistance to this mode of failure.

Figure 35
Cracking of concrete slab at top face, above the first opening, in beam CE2 (equivalent plastic strains) With respect to the other composite alveolar beams results (CE4, CE5, CA1, CA2, CA3), it can be observed that there was a gain not only in the initial stiffness but also in the ultimate load, when compared to the composite beam FW1. It occurred because in these beams the profiles failed by the formation of a flexural mechanism, for which increasing the height of the profile also increases its load capacity. Table 10 shows the values of the linear load applied at the beginning of the yielding of these beams and at the total yield of the lower tee-section. In the cases of the beams CE4 and CE5, the profile web has also reached its ultimate strength fu at localized points below the central openings, indicating its possible rupture, so that the values of applied linear load when it occurred are also identified in Table 10. In these seven beams, the failure process of the numerical models occurred in three stages: (i) yielding of the steel profile, (ii) significant increase in the shear force and slip of the connectors and (iii) concrete crushing or, alternatively, possible rupture of the steel profile (in beams CE4 and CE5). However, it is noteworthy that, since in the developed model it was not programmed an explicit failure criterion for the connectors, they may present large slips in the numerical model without failing, which does not necessarily correspond to reality, once the real failure can occur in the connectors before the concrete's crushing. Therefore, Table 11 presents the values of applied linear loads and midspan deflections when the first (most loaded) and the twentieth connectors reach 98% of their ultimate shear forces, which is equivalent to a slip of 2.05 mm, considering the parameters m = 1.9 and n = 1.0 used in eq. (2). These values do not necessarily correspond to the failure of the connectors, since they are ductile, but indicate that their rupture may be close. Table 11 also indicates the linear loads and midspan deflections obtained for the first crushed point on the upper face of the concrete slab. From Table 11 data, it can be verified that in the composite alveolar beams the first connector takes longer to be loaded, indicating a lower shear flow at the interface, when compared to composite beams FW1 and FW2. On the other hand, it can be observed that the applied linear load difference between the total loading of the 1 st and 20 th connectors becomes smaller, that is, once the first connector is completely loaded, the others quickly also become.
In summary, from this example it can be verified that the expansion of the profile may result in a significant gain in initial stiffness, and may also increase the ultimate load, as long as the web-post buckling is not the predominant failure mode. Thus, when designing composite alveolar beams, it is important to adopt an opening pattern that presents greater resistance to this failure mode, for example with larger web-posts and greater s/D0 ratios.

CONCLUSIONS
A finite element model was developed in ANSYS software, version 19.2, using updated material models and element types. This model was capable to simulate different failure modes of castellated beams and composite cellular beams. The obtained numerical results presented a good correlation with experimental and numerical results from previous works. Regarding the implemented modeling strategy, the following conclusions are obtained: • The addition of initial geometric imperfections by the weighted combination of elastic buckling modes, obtained in a preliminary Eigen-Buckling analysis, has shown to be an efficient approach for web-post buckling simulation.

•
The two employed material models for concrete -Usermat and DP-Concrete -presented quite similar results, which increases their reliability, considering that both were recently developed. However, the DP-Concrete model presented greater numerical stability and fewer convergence problems, especially when the Arc-Length method was used.
• For the solution, the Newton-Raphson method with displacement control proved to be quite efficient, being capable to capture the post-buckling behavior of composite alveolar beams, however it is limited to certain loading patterns. On the other hand, the Arc-Length method with force control is suitable for a wider range of loadings, but it presented some convergence difficulties in a few examples. Also, Newton-Raphson method with force control can be used, but in this case the convergence will stop at the peak load without capturing the post-buckling behavior.
• With respect to the effect of the openings on the behavior of composite beams, it is concluded that: • The structural gains due to steel profile expansion become more evident in larger spans.
• Load capacity gains may be limited when the composite alveolar beam fails by web-post buckling. Thus, when designing this type of beam, it is important to ensure a great resistance to this mode of failure. On the other hand, when the failure mode is the formation of a flexural mechanism, the load capacity gain is influenced by the expansion ratio and the tee-section height.
• The concrete cracking above the web-post buckling indicates that these two phenomena may be related.

RECOMMENDATIONS
The authors recommend that the numerical model validated in this paper be used for comparisons to analytical formulations of composite alveolar beams, aiming to help the process of future standardizations.
Additionally, the following topics could be further investigated: • The possible relation between the slab reinforcement and the composite beam resistance to web-post buckling.
• The effect of the steel profile expansion on the connectors' behavior.