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

Simulação numérica de vigas mistas de aço e concreto: estratégias atuais de modelagem por elementos finitos

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.

Keywords:
steel-concrete composite beams; numerical modeling; finite element method

resumo:

O emprego de vigas mistas de aço e concreto permite que sejam exploradas as melhores características mecânicas desses materiais, resultando em soluções estruturais mais econômicas. Vários pesquisadores estudaram o comportamento de vigas mistas a partir de diferentes estratégias de modelagem numérica, algumas das quais são apresentadas neste artigo. Dentro desse contexto, propõe-se a construção de um modelo numérico tridimensional com o software ANSYS, versão 19.2, utilizando elementos finitos atuais e modelos de materiais compatíveis com eles. Para a simulação do comportamento do concreto, foram utilizados dois modelos, o primeiro denominado DP-CONCRETE, disponibilizado nas versões mais recentes do ANSYS; e o segundo customizado com a interface USERMAT, e que tem por base o critério de Ottosen. Os resultados obtidos para as vigas analisadas apresentaram boa correlação com os resultados experimentais e com resultados numéricos de trabalhos anteriores.

Palavras-chave:
vigas mistas de aço e concreto; modelagem numérica; elementos finitos

1 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 [¹], 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 [²], [³].

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 [⁴] 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 [⁵] and the original Eurocode 4 [⁶] 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 [⁷], i.e., the stiffness criterion is adopted. In fact, according to Queirozet al. [⁸], 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 (Q_u), from a push-out test, standardized by Eurocode 4 [⁶]. 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 [³]. Between the 1960s and the 1970s, several experimental studies were carried out [⁹]–[¹¹]. 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 [¹], [⁸], [¹²]–[¹⁶].

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 [¹²]. The three-dimensional 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.

Tamayoet al. [¹³], Tamayo [¹⁴] and Dias [¹⁵] 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 [¹] and by Marconcinet al. [¹⁶], 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. Queirozet al. [⁸] 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.

Figure 1
Numerical models developed by previous works [¹], [⁸], [¹²]–[¹⁵].

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 [¹⁷]. 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 [¹⁸], 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 [¹], [⁸].

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 [¹⁸]. 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 [¹⁹]; 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 [⁹], which were also numerically simulated by previous works [¹], [⁸], [¹²], [¹⁴], [¹⁵].

2 NUMERICAL MODEL

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

2.1 Element types

The steel profile was modeled by quadrilateral shell elements, named SHELL181 [¹⁸], 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 [²⁰] and MacNeal and Harder [²¹], 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 [¹⁸]. In the present work these elements were used in their homogeneous form with full integration. Its formulation is based on Zienkiewiczet al. [²²]. 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 [¹⁸], 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 (u_x) 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 [¹⁸]. In this work, however, they were applied to nodes spaced by half the thickness of the superior flange (t_fs/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 u_y and u_z 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 Queirozet al. [⁸].

Figure 2
Detail of COMBIN39 nonlinear spring element for connectors modeling.

The reinforcement bars were modeled by discrete embedded elements, named REINF264 [¹⁸], 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 [¹⁸].

For the generation of reinforcement embedded elements, a new ANSYS functionality was used, denominated mesh-independent 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 [¹⁸]. 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.

Figure 3
Finite elements adopted for modeling the steel profile (SHELL181), the concrete slab (SOLID186) and the embedded reinforcement (REINF264) [¹⁸].

2.2 Material Models

The profile steel was modeled by von Mises yield criterion with isotropic hardening. As hardening law, the constitutive model proposed by Gattesco [¹²], 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.

$\sigma =f_y+E_h\cdot \left(\epsilon -{\epsilon }_h\right)\cdot \left(1-\frac{E_h\cdot \left(\epsilon -{\epsilon }_h\right)}{4\cdot \left(f_u-f_y\right)}\right)$ (1)

where f_y and f_u are the steel yield and ultimate strengths, ${\epsilon }_y=f_y/E$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 E_h is the tangent modulus of elasticity at the start of the hardening stage.

Figure 4
Constitutive relations: (a) profile steel; (b) reinforcement steel.

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. [²³], adjustable to the push-out test data of stud bolts.

$Q=Q_u\cdot {\left(1-e^{-ms}\right)}^n$ (2)

where Q is the shear force acting in the connector, in kN; Q_u 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:

1. i

   via a customized model with the USERMAT subroutine;

2. 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 [¹⁸], [²⁴]. The customized model used in the present work was developed by Lazzariet al. [²⁵]. The rupture surface adopted was proposed by Ottosen [¹⁹], and is given by Equations 3, 4 and 5.

$f\left(I_1,J_2,\cos \left(3\theta \right)\right)=\alpha \cdot \frac{J_2}{f_{cm}^2}+\lambda \cdot \frac{\sqrt{J_2}}{f_{cm}}+\beta \cdot \frac{I_1}{f_{cm}}-1=0$ (3)

$\lambda =\left\{\begin{array}{c}{c}_1\cdot \cos \left(\frac{1}{3}\arccos \left(c_2\cdot \cos \left(3\theta \right)\right)\right),if\cos \left(3\theta \right)\ge 0\\ c_1\cdot \cos \left(\frac{\pi }{3}-\frac{1}{3}\arccos \left(-c_2\cdot \cos \left(3\theta \right)\right)\right),if\cos \left(3\theta \right)\le 0\end{array}\right.$ (4)

$\cos \left(3\theta \right)=\frac{3\sqrt{3}}{2}\cdot \frac{J_3}{\sqrt{J_2^3}}$ (5)

where I₁ is the first stress invariant; J₂ and J₃ are the second and third deviatoric stress invariants; f_cm is the mean compressive strength of concrete; θ is the similarity angle; and α, β, c₁, c₂ are material parameters, automatically calculated by USERMAT from equations provided by the fib2010 model code [²⁶], which uses as input data the values of the mean uniaxial compressive strength of concrete (f_cm), biaxial compressive strength (f_c2m), tensile strength (f_ctm), 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 [²⁶].

$\frac{{\sigma }_c}{f_{cm}}=\frac{{\eta }^2-k\cdot \eta }{1+\left(k-2\right)\cdot \eta }$ (6)

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; E_ci is the tangent modulus of elasticity of concrete; E_c1 is the secant modulus of elasticity of concrete from the origin to the peak compressive stress; $k=E_{ci}/E_{c1}$is the plastic number; and $\eta ={\epsilon }_c/{\epsilon }_{c1}$.

Figure 5
Constitutive relations in the Usermat model: (a) compression; (b) tension.

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 (f_ctm) is reached. After the cracking occurs, the softening is governed by a decreasing line that intersects the vertical axis at the value of α.f_ctm 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 [²⁷]. 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 [¹⁸]. Figure 6 illustrates the two possible compositions, in the two-dimensional principal stresses system.

Figure 6
Combinations Drucker-Prager with Rankine or Drucker-Prager with Drucker-Prager (adapted from [¹⁸]).

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.

$f_{DPc}=\frac{{\beta }_c}{3}\cdot I_1+\sqrt{J_2}-{\sigma }_{Yc}\cdot {\Omega }_c=0$ (7)

$f_R=I_1+2\cdot \sqrt{3\cdot J_2}\cdot cos\left(\theta \right)-3\cdot T\cdot {\Omega }_t=0$ (8)

${\beta }_c=\frac{\sqrt{3}\cdot (R_b-R_c)}{2\cdot R_b-R_c}$ (9)

${\sigma }_{Yc}=\frac{R_b\cdot R_c}{\sqrt{3}\cdot \left(2\cdot R_b-R_c\right)}$ (10)

where R_c and R_b 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, $T=f_{ctm}$, $R_c=f_{cm}$and $R_b=f_{c2m}$ were adopted.

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.

Figure 7
Functions of the linear HSD model: (a) compression; (b) tension [¹⁸].

${\Omega }_c=\left\{\begin{array}{c}{\Omega }_{ci}+\left(1-{\Omega }_{ci}\right)\cdot \sqrt{2\cdot \frac{\kappa }{{\kappa }_{cm}}-\frac{{\kappa }^2}{{\kappa }_{cm}^2}},if0<\kappa \le {\kappa }_{cm}\\ \begin{array}{c}\\ 1-\frac{1-{\Omega }_{cr}}{{\kappa }_{cr}-{\kappa }_{cm}}\cdot \left(\kappa -{\kappa }_{cm}\right),if{\kappa }_{cm}<\kappa \le {\kappa }_{cr}\\ \end{array}\\ {\Omega }_{cr},if\kappa \ge {\kappa }_{cr}\end{array}\right.$ (11)

${\Omega }_t=\left\{\begin{array}{c}1-\left(1-{\Omega }_{tr}\right)\cdot \frac{\kappa }{{\kappa }_{tr}},if0<\kappa <{\kappa }_{tr}\\ \\ {\Omega }_{tr},if\kappa \ge {\kappa }_{tr}\end{array}\right.$ (12)

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 Equation 13 and Equation 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.

${\kappa }_{cm}={\epsilon }_{c1}-\frac{f_{cm}}{E_c}$ (13)

${\kappa }_{cr}={\epsilon }_{c,lim}-\frac{{\Omega }_{cr}\cdot f_{cm}}{E_c}$ (14)

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.

Figure 8
Fictitious plate of elastic material added at the interface between steel and concrete in the DP-Concrete model.

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

Figure 9
Numerical model of beam A2 with the adopted boundary conditions.

2.4 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 [¹], [⁸], [¹²]–[¹⁵]. It can be observed that the employed strategy, in terms of geometry and element types, approaches the strategy used by [⁸]: 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.

Table 1
Comparison between numerical models.

Model	Software	Steel profile	Connectors	Concrete slab
Gattesco [12]	Authorial	Bar elements. Nonlinear hardening.	Spring elements. Curve from push-out test.	Bar elements. Nonlinear constitutive model.
CEMACOM [13]–[15]	Authorial	Shell elements. von Mises criterion with linear hardening.	Bar elements. Curve from push-out test.	Degenerated shell elements. Modified Drucker-Prager model.
Kotinda [1]	ANSYS	Shell elements SHELL43. von Mises criterion with nonlinear hardening [12].	Bar elements BEAM189. von Mises criterion with linear hardening.	Hexahedral elements SOLID65. Model concrete [17].
Queirozet al. [8]	ANSYS	Shell elements SHELL43. von Mises criterion with nonlinear hardening [12].	Spring elements COMBIN39. Curve from push-out test.	Hexahedral elements SOLID65. Model concrete [17].
Present analysis	ANSYS	Shell elements SHELL181. von Mises criterion with nonlinear hardening [12].	Spring elements COMBIN39. Curve from push-out test.	Hexahedral elements SOLID186. Models: (i) Usermat (ii) DP-Concrete

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 [¹], [⁸], 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 [¹³]–[¹⁵], 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 [¹], [⁸] adopted perfect elastoplastic models, as in the present work, while [¹³]–[¹⁵] considered linear hardening.

3 ANALYZED EXAMPLES

Two beams that were experimentally tested by Chapman and Balakrishnan [⁹], 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 [¹], [⁸], [¹²], [¹⁴], [¹⁵]. 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 [⁹].

Figure 10
Geometry of the beams tested by Chapman and Balakrishnan [⁹].

Table 2
Materials data.

	Beam		A2	E1	Unit	Source
Connectors	Total number of connectors		76	100	-	Chapman and Balakrishnan [9]
	Number of rows		2	2	-	Chapman and Balakrishnan [9]
	Height		101.6	50.8	mm	Chapman and Balakrishnan [9]
	Diameter		19	12.7	mm	Chapman and Balakrishnan [9]
	Ultimate shear force (Qu)		120.6	63.76	kN	Chapman and Balakrishnan [9]
	Parameter m		1.2	1.15	1/mm	Ollgaardet al. [23]
	Parameter n		0.55	0.5	-	Ollgaardet al. [23]
Steel profile	Modulus of elasticity (E)	Flanges	20216	20520	kN/cm2	Chapman and Balakrishnan [9]
		Web	20368	19912	kN/cm2	Chapman and Balakrishnan [9]
	Yield strength (fy)	Flanges	23	24,5	kN/cm2	Chapman and Balakrishnan [9]
		Web	23.9	29.2	kN/cm2	Chapman and Balakrishnan [9]
	Ultimate strength (fu)	Flanges	41	45.8	kN/cm2	Chapman and Balakrishnan [9]
		Web	41.8	45.3	kN/cm2	Chapman and Balakrishnan [9]
	Ratio εh/εy	Flanges	4.9	2.2	-	Chapman and Balakrishnan [9]
		Web	7.6	1	-	Chapman and Balakrishnan [9]
	Tangent modulus (Eh)	Flanges	350	350	kN/cm2	Kotinda [1]
		Web	350	350	kN/cm2	Kotinda [1]
Concrete slab	Uniaxial compressive strength (fcm)		2.69	3.27	kN/cm2	Chapman and Balakrishnan [9]
	Uniaxial tensile strength (fctm)		0.21	0.25	kN/cm2	FIB [26]
	Biaxial compressive strength (fc2m)		3.16	3.82	kN/cm2	FIB [26]
	Stress σcom (on the compressive meridian)		-24	-24	kN/cm2	FIB [26]
	Stress τcom (on the compressive meridian)		15.3	15.1	kN/cm2	FIB [26]
	Initial tangent modulus of elasticity (Eci)		2990	3191	kN/cm2	FIB [26]
	Reduced modulus of elasticity (Ec)		2575	2790	kN/cm2	FIB [26]
	Secant modulus of elasticity (Ec1)		1300	1480	kN/cm2	FIB [26]
	Strain at the max. compressive stress (εc1)		2.1	2.2	‰	FIB [26]
	Ultimate compressive strain (εc,lim)		3.5	3.5	‰	FIB [26]
Slab Reinforcement	Modulus of elasticity (E)		20500	20500	kN/cm2	Kotinda [1]
	Yield strength (fy)		32	32	kN/cm2	Kotinda [1]

The numerical model of beam A2 is illustrated in Figure 9. The concrete's data that were not provided by Chapman and Balakrishnan [⁹] were calculated according to the model code fib2010 [²⁶], using the mean compressive strength. For the Usermat model it was used the initial tangent modulus of elasticity (E_ci), while for the DP-Concrete model it was used the reduced modulus of elasticity ($E_c={\alpha }_i.E_{ci}$), once the latter admits an elastic linear stage in compression until reaching the stress of 0.4.f_cm. The parameters m and n for the connectors curve, given in equation [²], were adjusted according to the results obtained from push-out tests, which were provided by Chapman and Balakrishnan [⁹], as shown in Figure 11. All adopted values for connectors, steel and concrete are presented in Table 2.

Figure 11
Curve fitting of Equation 2 to the experimental results of the push-out tests.

4 RESULTS AND DISCUSSION

Figure 12 and Figure 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 [¹], [⁸], [¹²], [¹⁴], [¹⁵] and the experimental results [⁹] are also presented.

Figure 12
Load-deflection curve – beam A2.

Figure 13
Load-deflection curve – beam E1.

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 Queirozet al. [⁸] and by Kotinda [¹]. 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 [¹], a different strategy for the connectors modeling. The CEMACOM model [¹⁴], [¹⁵] 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).

Figure 14
Comparison between the constitutive relations for compression used by the models for beam A2.

Figure 15
Top views of beam A2 concrete slab midspan: evolution of softening obtained with the DP-Concrete model for different load levels.

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.

4.1 Analysis of crushing

Both of analyzed beams failed by concrete crushing in the experimental tests. The failure loads registered by Chapman and Balakrishnan [⁹] 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 Figure 12 and Figure 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.

Table 3
Crushing loads.

Model		Beam A2		Beam E1
		Load P [kN]	(P-Pexp)/Pexp	Load P [kN]	(P-Pexp)/Pexp
Experimental	Pexp	441	-	510	-
Numerical (Usermat)	P1	382	-13.38%	450	-11.76%
	P2	458	3.85%	518	1.57%
Numerical (DP-Concrete)	P1	363	-17.69%	403	-20.98%
	P2	447	1.36%	549	7.65%

• Where: P₁ = load at the first crushed point; P₂ = load at 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.

Figure 16
Top views of beam A2 concrete slab midspan: evolution of crushing obtained with the DP-Concrete model for different load levels.

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 (Q_u) 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², while the ultimate strength of the steel is 45.3 kN/cm². Therefore, it is possible to conclude that the models in fact are capturing the concrete crushing, and not the failure of other components.

Figure 17
Von Mises stresses in beam E1 at the applied load of 518 kN (Usermat model).

4.2 Importance of adding a fictitious elastic plate in the DP-Concrete model

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.

Figure 18
Influence of the fictitious elastic plate at the interface on the beam E1 results.

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.

5 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 Lazzariet al. [²⁵], considering the Ottosen criterion [¹⁹] for the concrete modeling. Both were developed with finite elements classified as current technology by the ANSYS documentation, version 19.2 [¹⁸]. The obtained results for the beams A2 and E1 of Chapman and Balakrishnan [⁹] were quite satisfactory, presenting a good correlation with the experimental results, as with the results by numerical analysis of previous works [¹], [⁸], [¹²], [¹⁴], [¹⁵].

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.

ACKNOWLEDGEMENTS

The authors are grateful to CAPES and CNPq for the financial support for this research, and to CEMACOM/UFRGS for providing the infrastructure for the development of this work.

• Financial support: This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq). Finance Code 001
• How to cite: M. E. Benincá, R. J. Schmitz, and I. B. Morsch, “Numerical simulation of steel-concrete composite beams: updated strategies of finite element modeling”, Rev. IBRACON Estrut. Mater., vol. 13, no. 5, e13510, 2020, https://doi.org/10.1590/S1983-41952020000500010

REFERENCES

• ¹
  T. I. Kotinda, “Modelagem numérica de vigas mistas aço-concreto simplesmente apoiadas: ênfase ao estudo interface laje-viga,” M.S. thesis, Esc. Eng. São Carlos, Univ. São Paulo, São Carlos, 2006. http://dx.doi.org/10.11606/D.18.2006.tde-22062006-110301
  » http://dx.doi.org/10.11606/D.18.2006.tde-22062006-110301
• ²
  J. Nie, V. Xiao, and L. Chen, "Experimental studies on shear strength of steel-concrete composite beams," J. Struct. Eng., vol. 130, no. 8, pp. 1206–1213, 2004, http://dx.doi.org/10.1061/(ASCE)0733-9445(2004)130:8(1206)
  » http://dx.doi.org/10.1061/(ASCE)0733-9445(2004)130:8(1206)
• ³
  W. Pfeil and M. Pfeil, Estruturas de Aço – Dimensionamento Prático de Acordo com a NBR 8800:2008, 8th ed. Rio de Janeiro: LTC, 2013.
• ⁴
  Associação Brasileira de Normas Técnicas, Projeto de Estruturas de Aço e de Estruturas Mistas de aço e Concreto de Edifícios, ABNT NBR 8800, 2008.
• ⁵
  Instituto Português da Qualidade, Projecto de Estruturas Mistas Aço-Betão – Parte 1-1 – Regras Gerais e Regras para Edifícios, NP EN 1994-1-1, 2011.
• ⁶
  European Committee for Standardization, Eurocode 4: Design of Composite Steel and Concrete Structures – Part 1-1: General Rules and Rules for Buildings, CEN EN 1994-1-1, 2004.
• ⁷
  A. R. Silva, “Análise numérica de elementos estruturais com interação parcial,” Ph.D. dissertation, Univ. Fed. Ouro Preto,, Ouro Preto, 2010. [Online]. Available: http://www.repositorio.ufop.br/handle/123456789/2254
  » http://www.repositorio.ufop.br/handle/123456789/2254
• ⁸
  F. D. Queiroz, P. C. G. S. Vellasco, and D. A. Nethercot, "Finite element modelling of composite beams with full and partial shear connection," J. Construct. Steel Res., vol. 63, no. 4, pp. 505–521, 2007, http://dx.doi.org/10.1016/j.jcsr.2006.06.003
  » http://dx.doi.org/10.1016/j.jcsr.2006.06.003
• ⁹
  J. C. Chapman and S. Balakrishnan, "Experiments on composite beams," Struct. Eng., vol. 42, no. 11, pp. 369–383, 1964.
• ¹⁰
  R. G. Slutter and G. C. J. Driscoll, "Flexural strength of steel-concrete composite beams," J. Struct. Div., vol. 91, no. 2, pp. 71–99, 1965.
• ¹¹
  J. A. Grant, R. Slutter, and J. W. Fisher, "Composite beams with formed steel deck," Eng. J., vol. 14, pp. 24–43, 1977.
• ¹²
  N. Gattesco, "Analytical modeling of nonlinear behavior of composite beams with deformable connection," J. Construct. Steel Res., vol. 52, no. 2, pp. 195–218, 1999, http://dx.doi.org/10.1016/S0143-974X(99)00026-7
  » http://dx.doi.org/10.1016/S0143-974X(99)00026-7
• ¹³
  J. L. P. Tamayo, I. B. Morsch, and A. M. Awruch, "Short-time numerical analysis of steel–concrete composite beams," J. Braz. Soc. Mech. Sci. Eng., vol. 37, no. 4, pp. 1097–1109, 2015, http://dx.doi.org/10.1007/s40430-014-0237-9
  » http://dx.doi.org/10.1007/s40430-014-0237-9
• ¹⁴
  J. L. P. Tamayo, “Análise numérica de vigas mistas pelo método dos elementos finitos,” M.S. thesis, Prog. Pós-grad. Eng. Civ., Univ. Fed. Rio Grande do Sul, Porto Alegre, 2011. [Online]. Available: https://hdl.handle.net/10183/30117
  » https://hdl.handle.net/10183/30117
• ¹⁵
  M. M. Dias, “Análise numérica de vigas mistas aço-concreto pelo método dos elementos finitos: efeitos de longa duração,” M.S. thesis, Prog. Pós-grad. Eng. Civ., Univ. Fed. Rio Grande do Sul, Porto Alegre, 2013. [Online]. Available: https://hdl.handle.net/10183/77726
  » https://hdl.handle.net/10183/77726
• ¹⁶
  L. R. Marconcin, R. D. Machado, and M. A. Marino, "Numerical modeling of steel-concrete composite beams," IBRACON Struct. Mater. J., vol. 3, no. 4, pp. 449–476, 2010, http://dx.doi.org/10.1590/S1983-41952010000400006
  » http://dx.doi.org/10.1590/S1983-41952010000400006
• ¹⁷
  K. J. Willam and E. P. Warnke, “Constitutive model for the triaxial behavior of concrete,” in IABSE Seminar on Concrete Structures Subjected to Triaxial Stresses, Bergamo, Italy, 1974.
• ¹⁸
  ANSYS Inc., ANSYS Help System – Version 19.2, 2018.
• ¹⁹
  N. S. Ottosen, "A failure criterion for concrete," J. Eng. Mech. Dvision, vol. 103, no. 4, pp. 527–535, 1977.
• ²⁰
  K. J. Bathe and E. N. Dvorkin, "A formulation of general shell elements: the use of mixed interpolation of tensorial components," Int. J. Numer. Methods Eng., vol. 22, no. 3, pp. 697–722, 1986, http://dx.doi.org/10.1002/nme.1620220312
  » http://dx.doi.org/10.1002/nme.1620220312
• ²¹
  R. H. MacNeal and R. L. Harder, "A refined four-noded membrane element with rotational degrees of freedom," Comput. Struc., vol. 28, no. 1, pp. 75–84, 1988, http://dx.doi.org/10.1016/0045-7949(88)90094-6
  » http://dx.doi.org/10.1016/0045-7949(88)90094-6
• ²²
  O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, 7th ed. Oxford, UK: Elsevier, 2013.
• ²³
  J. G. Ollgaard, R. G. Slutter, and J. W. Fisher, "Shear strength of stud connectors in lightweight and normalweight concrete," Eng. J., vol. 8, pp. 55–64, 1971. [Online]. Available: https://preserve.lehigh.edu/engr-civil-environmental-fritz-lab-reports/2010
  » https://preserve.lehigh.edu/engr-civil-environmental-fritz-lab-reports/2010
• ²⁴
  F. P. M. Quevedo, R. J. Schmitz, I. B. Morsch, A. Campos Filho, and D. Bernaud, "Customization of a software of finite elements to analysis of concrete structures: long-term effects," IBRACON Struct. Mater. J., vol. 11, no. 4, pp. 696–718, 2018, http://dx.doi.org/10.1590/s1983-41952018000400005
  » http://dx.doi.org/10.1590/s1983-41952018000400005
• ²⁵
  P. M. Lazzari, A. Campos Fo., B. M. Lazzari, A. R. Pacheco, and R. R. S. Gomes, "Numerical simulation of the constructive steps of a cable-stayed bridge using ANSYS," Struct. Eng. Mech., vol. 69, no. 3, pp. 269–281, 2019., http://dx.doi.org/10.12989/sem.2019.69.3.269
  » http://dx.doi.org/10.12989/sem.2019.69.3.269
• ²⁶
  Fédération Internationale du Béton, Fib Model Code for Concrete Structures 2010, 2013.
• ²⁷
  M. Martinelli, “Modelagem de situações de punção em lajes de concreto armado através do método dos elementos finitos,” M.S. thesis, Prog. Pós-grad. Eng. Civ., Univ. Fed. Rio Grande do Sul, Porto Alegre, 2003. [Online]. Available: https://hdl.handle.net/10183/2299
  » https://hdl.handle.net/10183/2299

Edited by

Editors: José Marcio Calixto, José Luiz Antunes de Oliveira e Sousa, Guilherme Aris Parsekian

Publication Dates

• Publication in this collection
  16 Oct 2020
• Date of issue
  2020

History

• Received
  15 Apr 2019
• Accepted
  08 Feb 2020