Numerical simulation of steel-concrete composite beams: updated strategies of finite element modeling

abstract: The use of steel-concrete composite beams allows the best properties of these materials to be explored, resulting in more economical solutions. Many researchers have studied the behavior of composite beams from different strategies of numerical modeling, and some of these are presented in this article. In this context, the present work proposes the construction of a tridimensional numerical model using ANSYS software, version 19.2, with current-technology elements and compatible material models. For the simulation of concrete behavior, two models have been used: the first, denominated DP-CONCRETE, is a native ANSYS model, available in the more recent versions of this software; and the second, denominated USERMAT, is a customizable model that was developed based on Ottosen criterion. The results obtained with these models for the analyzed beams presented a good correlation with the experimental results and with numerical results from previous works.


INTRODUCTION
The use of steel-concrete composite beams has been widespread in the civil construction sector, as it allows the best mechanical properties of the involved materials to be explored. In these beams, a steel profile is connected to a concrete slab through shear connectors, which aim to restrict the longitudinal slip and the vertical separation at the interface [1], ensuring the joint work of the elements. This configuration allows to use these materials more rationally, especially in simply supported beams, where the concrete works basically under compression and the steel profile under tension. The result of this rationalization is an increase in the load carrying capacity of the composite beam, when compared to the materials working separately, allowing the design of larger spans and/or the achievement of more economical structural solutions [2], [3].
There are different possible shapes for the shear connectors, and the most common are known as stud bolts. The level of interaction between the steel profile and the concrete slab is a fundamental issue in the design and analysis of composite beams, and can be classified according to two different criteria: strength or stiffness. The Brazilian code NBR 8800 [4] uses the first of them, defining that "complete interaction" occurs between steel and concrete when the connectors present a design strength equal to or greater than the design strength of the steel profile under tension or the concrete slab under compression, whichever is less; and, otherwise, "partial interaction" occurs. Also adopting the strength criterion, but in other words, the Portuguese version of Eurocode NP EN 1994-1-1 [5] and the original Eurocode 4 [6] use the terms "partial shear connection" and "full shear connection", respectively, to distinguish situations in which an increase in the resistance of the longitudinal shear connection increases, or not, the design bending resistance of the composite beam. However, it is worth mentioning that, in the specialized literature, the term "partial interaction" is often used to characterize a deformable connection between two structural elements [7], i.e., the stiffness criterion is adopted. In fact, according to Queiroz et al. [8], a composite beam presents "full interaction" when the connectors can be admitted as infinitely stiff, thus the relative slip between the steel profile and the concrete slab next to the connections is negligible; and that presents "partial interaction" when the deformation of the connectors is significant, resulting in slips that cannot be neglected. It is noteworthy that, in practice, in cases of full connection (or complete interaction, strength criterion), it is common to adopt, for design purposes, the hypothesis of full interaction (stiffness criterion) -reason why the cited codes only use the strength criterion in their classifications. Anyway, the real behavior of a connector can be represented by a slip versus shear force curve (δxQ), which can be obtained, along with the ultimate shear force (Qu), from a push-out test, standardized by Eurocode 4 [6]. Therefore, numerical models that incorporate the results of this test to the connectors assume the hypothesis of partial interaction (stiffness criterion) and, therefore, can simulate more realistically the behavior of composite beams, for both partial and full shear connection (strength criterion).
In the first applications, the steel profile of the composite beam was designed to support all the acting loads, assuming that the profile and the slab worked independently. Due to the scarcity of steel in the post II World War period, engineers began to design these beams considering the contribution of the concrete slab, starting the movement for systematic research on the subject [3]. Between the 1960s and the 1970s, several experimental studies were carried out [9]- [11]. More recently, studies using the finite element method for the numerical analysis of these beams have been multiplied, with the construction of models that contemplate the partial interaction [1], [8], [12]- [16].
There are different ways to create a finite element model for the analysis of composite beams. The one-dimensional models are those in which the steel profile and the concrete slab are modeled as bar elements, being connected by spring elements to simulate the shear connectors. An example of this approach is the work of Gattesco [12]. The threedimensional models are those that use shell elements to model the steel profile, and may present different approaches for modeling the concrete slab and shear connectors.
Tamayo et al. [13], Tamayo [14] and Dias [15] developed a computational code at CEMACOM/UFRGS for the analysis of composite beams, considering a three-dimensional finite element model in which the concrete slab was modeled by 8-node degenerated shell elements, the shear connectors by bar elements with penalties and the steel profile by shell elements. ANSYS software was employed by Kotinda [1] and by Marconcin et al. [16], with the steel profile being modeled by shell elements (named SHELL43), the concrete slab by hexahedral elements (SOLID65), the shear connectors by bar elements (BEAM189), and using pairs of contact elements (TARGE170 and CONTA173) at the interface between steel and concrete. Queiroz et al. [8] also used ANSYS and modeled the steel profile by shell elements (SHELL43), the concrete slab by hexahedral elements (SOLID65) and the connectors by non-linear spring elements (COMBIN39). Figure 1 illustrates the different models mentioned.
The 8-node element SOLID65, employed by the mentioned works, received notoriety for being the only element that can be associated with the constitutive model CONCRETE, made available by ANSYS, and based on the formulation of Willam and Warnke [17]. However, although the software still reads the codes that use this element, it is now classified as a "legacy element" by the ANSYS documentation, version 19.2 [18], being incompatible with several current functionalities, such as the use of embedded elements for reinforcement. In addition, numerical instability problems associated with this model were described by several researchers, who reported the need to disable the crushing option to improve convergence [1], [8].  [1], [8], [12]- [15].
In this context, the present work proposes the construction of a three-dimensional finite element model for the numerical simulation of composite beams, with ANSYS software, version 19.2, using elements classified as current technology by the software documentation [18]. In this model, the applicability of two different approaches for the simulation of concrete behavior will be evaluated: (i) a customized model by the subroutine USERMAT, based on the Ottosen criterion [19]; and (ii) the DP-CONCRETE model, available in the most recent versions of ANSYS. The validation of the developed model will be carried out from the numerical simulation of two composite beams tested experimentally by Chapman and Balakrishnan [9], which were also numerically simulated by previous works [1], [8], [12], [14], [15].

NUMERICAL MODEL
The following items present information about the element types, material models and boundary conditions used in the developed finite element model.

Element types
The steel profile was modeled by quadrilateral shell elements, named SHELL181 [18], with four nodes and six degrees of freedom per node (x, y, z translations and rotations), considering both membrane and bending stiffnesses. Its formulation is based on the work of several authors, including Bathe and Dvorkin [20] and MacNeal and Harder [21], and uses Reissner-Mindlin first-order shear-deformation theory. It is applicable to linear and nonlinear problems, including large deformations and rotations. For this reason, the formulation uses logarithmic strain and true stress measures rather than nominal engineering strain and stress. For small deformations, the difference between nominal and true values is negligible. The element can also contain several layers, but in this paper a single and centralized layer was used, with five integration points along the thickness.
The concrete slab was modeled by hexahedral elements with twenty nodes and three degrees of freedom per node (translations in x, y and z), named SOLID186 [18]. In the present work these elements were used in their homogeneous form with full integration. Its formulation is based on Zienkiewicz et al. [22]. Because it is a current-technology element, SOLID186 is compatible with several current ANSYS features, such as the generation of embedded elements and the use of new material models, e.g. DP-CONCRETE.
The connectors were modeled by nonlinear spring elements, named COMBIN39 [18], acting in the longitudinal direction of the beam. Thus, the relative slip between the nodes of the steel profile and the concrete slab in this direction (ux) are governed by points of the force versus slip curve of the element, whose data are obtained from push-out tests. When acting on a single degree of freedom per node, these elements should preferably be applied to coincident nodes [18]. In this work, however, they were applied to nodes spaced by half the thickness of the superior flange (tfs/2), because of the model geometry and the centralized positioning of the profile shell elements cross-section, as shown in Figure 2. However, once this distance is very small, it can be assumed that these nodes are practically coincident. At the same time, in the transverse and vertical directions, the couplings of the displacements uy and uz of these same nodes are performed. Thus, the COMBIN39 element acts only in the longitudinal direction of the beam, and compatibility equations are responsible for simulating the behavior of connectors in the other directions, as in the model proposed by Queiroz et al. [8]. The reinforcement bars were modeled by discrete embedded elements, named REINF264 [18], which are suitable for simulating steel bars. These elements use the same nodes of the base elements SOLID186, even if their geometric position does not coincide with them. The REINF264 element presents only axial stiffness, thus the stiffnesses to bending, torsion and shear are neglected. A perfect interaction between the reinforcing element and the concrete base element is admitted, so there is no relative movement between them [18].
For the generation of reinforcement embedded elements, a new ANSYS functionality was used, denominated meshindependent method. This method uses MESH200 elements, which are only guide elements, thus do not directly contribute to the solution, but determine the positions where REINF264 reinforcement elements are created [18]. Therefore, it is possible to insert the positions of the reinforcement bars from the lines drawn in absolute coordinates, unlike the standard method, in which it is necessary to use relative coordinates in respect to the base elements, generating mesh dependence.
The elements SHELL181, SOLID186 and REINF264 are illustrated in Figure 3.

Material Models
The profile steel was modeled by von Mises yield criterion with isotropic hardening. As hardening law, the constitutive model proposed by Gattesco [12], shown in Figure 4a, was adopted. 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, is the yield strain, εh is the strain at the initial hardening, εu is the strain at the ultimate stress, E is the modulus of elasticity and Eh is the tangent modulus of elasticity at the start of the hardening stage.
On the other hand, the steel of the reinforcement bars was simplified as perfectly elastoplastic, as shown in Figure 4b. Differently from the other materials, the steel properties of the connectors are not inserted into a material model itself, but rather through real constants of the COMBIN39 spring element. These real constants provide points for the force versus slip curve, obtained experimentally via push-out test. To systematize the model, the theoretical curve given by Equation 2 was adopted, as proposed by Ollgaard et al. [23], adjustable to the push-out test data of stud bolts.

( )
where Q is the shear force acting in the connector, in kN; Qu is the ultimate shear force resisted by the connector, in kN; s is the slip, in mm; e = 2.718 is the Euler number; and m, n are curve fitting parameters, in mm -1 and dimensionless, respectively.
For each analyzed example, two numerical models were developed, differing in the way they simulate the concrete's behavior: (i) via a customized model with the USERMAT subroutine; (ii) via an ANSYS native model, called DP-CONCRETE, recently made available by the software (since 17.0 version).
USERMAT is an ANSYS subroutine that can be programmed to customize a material model [18], [24]. The customized model used in the present work was developed by Lazzari et al. [25]. The rupture surface adopted was proposed by Ottosen [19], and is given by Equations 3, 4 and 5.
( ) ( ) , , cos where I1 is the first stress invariant; J2 and J3 are the second and third deviatoric stress invariants; fcm is the mean compressive strength of concrete; θ is the similarity angle; and α, β, c1, c2 are material parameters, automatically calculated by USERMAT from equations provided by the fib2010 model code [26], which uses as input data the values of the mean uniaxial compressive strength of concrete (fcm), biaxial compressive strength (fc2m), tensile strength (fctm), and triaxial compressive strength, the latter being defined by a point located on the compressive meridian and described by the stresses σcom and τcom.
For the compressive behavior, the hardening law adopted is given by Equation 6, illustrated in Figure 5a, as suggested by fib2010 model code [26].
where σc is the compressive stress; εc is the compressive strain of concrete; εc1 is the strain at the maximum compressive stress; εc,lim is the ultimate compressive strain; Eci is the tangent modulus of elasticity of concrete; Ec1 is the secant modulus of elasticity of concrete from the origin to the peak compressive stress; For the tensile behavior, a model of distributed cracks was adopted, considering the tension stiffening effect through the softening law shown in Figure 5b. Initially, the concrete is admitted as linear elastic until the tensile strength (fctm) is reached. After the cracking occurs, the softening is governed by a decreasing line that intersects the vertical axis at the value of α.fctm and the horizontal axis at the limit strain value (εctu). In the present work, the parameters α = 0.6 and εctu = 0.001 were adopted, as suggested by Martinelli [27]. The customized model also contemplates the reduction of shear stresses transfer in the crack, considering the aggregate interlock mechanism and the reinforcement pin effect, from the reduction of the cracked concrete's transversal elasticity modulus. The second model used to simulate the concrete's behavior, denominated DP-CONCRETE, is a composition of two distinct surfaces: a Drucker-Prager yield surface for compression, and a second surface, which may be Drucker-Prager or Rankine, for tension and tension-compression. It is necessary to use one of these combinations since a single Drucker-Prager surface does not represent the large differences in tensile and compressive behavior of concrete [18]. Figure 6 illustrates the two possible compositions, in the two-dimensional principal stresses system. Drucker-Prager (adapted from [18]).
In this paper, the combination of Drucker-Prager (compression) with Rankine (tension) was adopted. These surfaces are defined, respectively, by Equations 7 and 8. The parameters βc and σYc of Equation 7 are calculated with Equations 9 and 10.
( ) where Rc and Rb are the uniaxial and the biaxial compressive strengths of concrete, respectively; T is the concrete tensile strength; Ωc and Ωt are hardening and softening functions in compression and in tension. In this work, In function of the HSD (hardening, softening and dilatation) model adopted, the functions Ωc and Ωt assume a certain shape. With these models it is possible to simulate approximately the phenomena of cracking, in tensile behavior, and crushing, in compressive behavior, from increments of plastic strains related to the hardening and softening rules. ANSYS provides four types of HSD models. In this work, the "linear model" was adopted, which is governed by Equations 11 and 12, for compression and tension, respectively. These functions are illustrated in Figure 7.
where κ is the effective plastic strain; κcm is the plastic strain when the maximum compressive stress is reached; κcr is the ultimate plastic strain in compression; Ωci is the relative tension at the start of the plastification; Ωcr is the residual relative stress in compression; κtr is the plastic strain when the residual tensile stress is reached; Ωtr is the residual relative stress in tension. 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 with Equations 13 and 14, 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. , Given the fact that in the DP-CONCRETE model the cracking is simulated by increments of plastic strains, it is recommended to use elastic materials at stress concentration points due to external loads, such as at the supports, to avoid an unreal level of strains. In this work, the vertical supports are applied directly to the steel profile, so they do not constitute a problem in this regard. However, the connectors, being represented by discrete spring elements, impose concentrated shear loads on the steel-concrete interface. Thus, to avoid the excessive increase of plastic deformations in this region, it was decided to add a fictitious plate on the lower face of the slab, along the interface, formed by SHELL281 8-node shell elements, with 1 mm thickness, as shown in Figure 8. The plate material is admitted as linear elastic and has the same modulus of elasticity as the concrete.

Boundary Conditions
The models were developed considering symmetry, thus only half of the beam was modeled. The nodes located near the first support, at the inferior flange of the profile, had the displacements in y and z constrained; and all the nodes located at the midspan had the displacements in x and rotation in y and z constrained. The total applied load on the experiment was divided by two, because of the symmetry, and then it was distributed evenly on the central nodes on the superior face of the concrete slab. These boundary conditions are highlighted in Figure 9. For the solution of the nonlinear system, it was adopted the full Newton-Raphson method.

Comparison to other numerical models
In Table 1 the main characteristics of the developed models are presented, and they are compared to models from previous works [1], [8], [12]- [15]. It can be observed that the employed strategy, in terms of geometry and element types, approaches the strategy used by [8]: the main differences are in the adopted material models for the concrete and in the choice of more current elements. However, it maintains the same constructive logic, using hexahedral solid elements for the slab, spring elements for the connectors and shell elements for the steel profile. In the present work, discrete embedded elements were used for modeling the reinforcement steel bars of the concrete slab, differently from other works that employed ANSYS [1], [8], in which layers of distributed reinforcement were adopted, once the element SOLID65, used by them, is not compatible with discrete embedded elements. In the CEMACOM model [13]- [15], in which the slab is modeled by degenerated shell elements, the reinforcement is also distributed in layers. For the constitutive relation of the reinforcement steel [1], [8] adopted perfect elastoplastic models, as in the present work, while [13]- [15] considered linear hardening.

ANALYZED EXAMPLES
Two beams that were experimentally tested by Chapman and Balakrishnan [9], denominated A2 and E1, were numerically analyzed in this article. The selection of these beams is justified for comparison of results, once they were also numerically analyzed by previous works [1], [8], [12], [14], [15]. The geometry of the two beams are the same, and the differences between them are in the connectors disposition and in the materials data, as presented in Figure 10 and in Table 2. More details about the tests can be obtained in Chapman and Balakrishnan [9].
The numerical model of beam A2 is illustrated in Figure 9. The concrete's data that were not provided by Chapman and Balakrishnan [9] were calculated according to the model code fib2010 [26], using the mean compressive strength. For the Usermat model it was used the initial tangent modulus of elasticity (Eci), while for the DP-Concrete model it was used the reduced modulus of elasticity ( , once the latter admits an elastic linear stage in compression until reaching the stress of 0.4.fcm. The parameters m and n for the connectors curve, given in equation [2], were adjusted according to the results obtained from push-out tests, which were provided by Chapman and Balakrishnan [9], as shown in Figure 11. All adopted values for connectors, steel and concrete are presented in Table 2. Figure 10. Geometry of the beams tested by Chapman and Balakrishnan [9].

RESULTS AND DISCUSSION
Figures 12 and 13 present the diagram of the applied load versus vertical displacement at midspan for the beams A2 and E1, respectively. Besides the obtained results from the developed model, the results obtained by previous works [1], [8], [12], [14], [15] and the experimental results [9] are also presented.  It is observed that both of developed models in this work properly simulated the global behavior of beams A2 and E1, with the DP-Concrete model achieving the best results for beam A2, and the Usermat model for beam E1. In comparison to the models from previous works, it is verified that, in both analyzed beams, the Usermat model presents results relatively near to the obtained by Queiroz et al. [8] and by Kotinda [1]. These authors also employed ANSYS, despite the facts that they used different finite elements for profile and slab, a different material model for concrete, and, in the case of Kotinda [1], a different strategy for the connectors modeling. The CEMACOM model [14], [15] presented a greater stiffness in the plastic regime, with smaller displacements for the same load levels, but simulated with good precision the experimental failure loads.
The models with Usermat and with DP-Concrete presented almost identical results in the elastic regime, and differed in the final loading stage, when the concrete reaches high stresses. This difference of behavior in high compression stresses was already expected, once, after the maximum compressive stress is reached, the DP-Concrete model presents linear softening, while the Usermat model presents a curvy softening, as illustrated in Figure 14, that makes a comparison between the constitutive relation in compression used by the two models for the concrete of beam A2. Therefore, the Usermat model is stiffer in the plastic regime. It can be observed, in Figure 12, that the small differences between the obtained results for beam A2 start at the applied load of 300 kN and become more significant from the applied load of 350 kN. When analyzing the equivalent plastic deformation in DP-Concrete model, it was verified that the compressed elements at the midspan started the softening process from the applied load of 343 kN. In Figure 15, the top views of beam A2 for different load levels are presented, and it is indicated the points in which the equivalent plastic strains exceed the value of the plastic strain associated to the maximum stress on the compression constitutive model (κcm).
It can be observed that the bigger the number of elements in softening, the more significant becomes the difference between the results of the models. However, it is important to note that other factors may influence the response, for example the shape of the yield surfaces in compression (Drucker-Prager or Ottosen) and the cracking in the regions of the slab under tension. In fact, as shown in Figure 12, differences in the results already exist before the compressive softening, for the applied load between 300 kN and 343 kN. Nevertheless, these initial differences are small, and practically insignificant in comparison to the differences after the softening. Therefore, it is possible to conclude that the smaller stiffness of the DP-Concrete model in the plastic regime can be mainly justified due to the differences between the softening curves in compression.

Analysis of crushing
Both of analyzed beams failed by concrete crushing in the experimental tests. The failure loads registered by Chapman and Balakrishnan [9] were 441 kN and 510 kN for beams A2 and E1, respectively.
In Usermat, it was programed the generation of a text file that lists the elements and their respective crushed integration points at each load substep. In this model, the crushing criterion is presented in relation to deformations. When the equivalent strain exceeds the ultimate strain in compression, adopted as 3.5/1000 in this paper, it admits that the integration point has crushed, and its stresses are zeroed. Analyzing the generated text file, it is verified that the first integration points were crushed at the loads of 382 kN and 450 kN for beams A2 and E1, respectively. These loads are highlighted in Figures 12 and 13. In subsequent substeps, a bigger number of elements and integration points have crushed, and the iterative process continues until ANSYS cannot reach the balance between internal and external forces, so the solution is interrupted. The loads reached at the end of convergence were 458 kN and 518 kN, for beams A2 and E1, respectively, as it is presented in Table 3. It can be observed, for beams A2 and E1, that the crushing failure loads obtained by the experimental tests are between the numerical values of the first crushed integration point and the end of convergence.
In DP-Concrete model, on the other hand, it is not possible to program the generation of a text file to control de number of crushed integration points at each substep, once this model is native from ANSYS, thus it cannot be customized. An alternative to control the crushing is to do a graphical analysis of the equivalent plastic strain evolution on the compressed elements. In this paper, it was adopted a simplified criterion that when the equivalent plastic strain of each point exceeds the ultimate plastic strain in compression (κcr), the crushing starts. Top views of beam A2 slab are presented in Figure 16, where the evolution of crushed points is illustrated for different load levels: 363 kN (first points are crushed), 379 kN, 395 kN and 447 kN (end of convergence). As can be observed, despite reaching the load of 447 kN at the last calculated substep, at this stage the crushing is well advanced, with many crushed elements. Therefore, it can be concluded that the numerical model indicates that the failure of beam A2 may occur before, for a load between 363 kN and 447 kN. The experimental failure load, equal to 441 kN, is contained in this interval. The same analysis was performed for beam E1, and it was verified that the first points were crushed at the load of 403 kN and that the end of convergence occurred at 549 kN. Again, the experimental failure load, in this case equal to 510 kN, is contained in this interval. The mentioned results are presented in Table 3. It is worth mentioning that, for both analyzed beams, the load at the end of convergence was bigger than the experimental failure load, which attests the importance to perform a critical analysis of the strain's evolution in the slab. Also, it is worth mentioning that the end of convergence can be considerably influenced by numerical parameters adopted, for example: error tolerance, number of substeps, number of iterations at each substep, among others -one more reason why the failure load obtained by a numerical model must be carefully analyzed.  The convergence can also stop due to the failure of other component of the structure. Therefore, beyond verifying the strains in the concrete slab, it is recommendable to analyze the evolution of shear forces in the connectors and the stresses in the steel profile, in order to certify that the model really captured the crushing of concrete. Regarding the connectors, it was verified that in all analyses the maximum shear forces did not exceed 55% and 89% of the ultimate shear forces (Qu) obtained by push-out tests for beams A2 and E1, respectively. Regarding the steel profile, despite the plastification had developed significantly in the region near the midspan at the end of convergence in both beams, it was verified that the von Mises stresses did not exceed, in any analyzed case, the value of 80% of the steel ultimate strength. For illustration, Figure 17 presents the von Mises stresses in the steel profile of beam E1, obtained by the numerical model with Usermat, for the failure load 518 kN: the maximum stress is 34.14 kN/cm 2 , while the ultimate strength of the steel is 45.3 kN/cm 2 . Therefore, it is possible to conclude that the models in fact are capturing the concrete crushing, and not the failure of other components.  Figure 18 presents the curves of applied load versus vertical displacement at midspan obtained for beam E1, with and without the fictitious plate, composed by elastic material, added to the inferior face of the slab, at the interface between steel and concrete. It can be observed that the plate practically does not change the global behavior of the Usermat model, but it has significant importance in the DP-Concrete model. This fact occurs because, without the plate, the plastic strains in the cracked concrete near the connectors increase in an excessive way, in function of the stress concentrations caused by the concentrated shear forces imposed by the discrete spring elements. In this case, the stresses are not correctly distributed between profile and slab, and the composite beam stiffness decreases considerably when its behavior becomes nonlinear. The same does not occur with the Usermat model, in which the crack is simulated by changing the stiffness matrices of the cracked finite elements.

CONCLUSIONS
Two numerical models were developed in ANSYS to analyze composite beams, one with a material model recently made available by the software, denominated DP-Concrete, and the other with a customized model programed in the interface Usermat, developed by Lazzari et al. [25], considering the Ottosen criterion [19] for the concrete modeling. Both were developed with finite elements classified as current technology by the ANSYS documentation, version 19.2 [18]. The obtained results for the beams A2 and E1 of Chapman and Balakrishnan [9] were quite satisfactory, presenting a good correlation with the experimental results, as with the results by numerical analysis of previous works [1], [8], [12], [14], [15].
Two alternatives were proposed to carefully analyze the concrete crushing, instead of simply considering the load at the end of convergence as the failure load of the composite beam. In the customized model with Usermat it was programmed the generation of a text file that lists the crushed integration points at each calculated substep. Thus, it is possible to verify when the first point is crushed. Based on a similar idea, the start of crushing in DP-Concrete model was analyzed with a graphical investigation of the equivalent plastic strains at each substep. The delimitated intervals between the start of crushing and the end of convergence, obtained with both numerical models, contain the respective experimental failure loads, for both analyzed beams. It was verified that the other components of the composite beams did not reach their respective failure stresses. Therefore, it is possible to conclude that the numerical models were able to capture the concrete crushing phenomenon.
The DP-Concrete model presented good results when a fictitious elastic plate was added at the interface between steel and concrete, aiming to avoid the excessive increase of plastic strains due to the cracking in this region. Without this plate, the quality of results decreased considerably. It was also verified that the same plate, when added to the Usermat model, did not change the results significantly.
At last, it is concluded that the proposed strategies of finite element modeling are proper to the numerical simulation of composite beams. The possibility of using two different material models to simulate the concrete's behavior allows the comparison of results between them, increasing the reliability of the analyses.