A 3D solid-like frame finite element applied to steel structures under high temperatures

Finite Elements in Analysis and Design 91 (2014) 68–83

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

A 3D solid-like frame finite element applied to steel structures under high temperatures Ronaldo Rigobello b, Humberto Breves Coda a,n, Jorge Munaiar Neto a a Structural Engineering Department, São Carlos School of Engineering, University of São Paulo, Av Trabalhador São Carlense, 400, ZIP-13560-590 São Carlos, SP, Brazil b Federal University of Technology – Paraná, Campus Campo Mourão, BR 369 – km 0.5, ZIP-87301-006 Campo Mourão, PR, Brazil

art ic l e i nf o

a b s t r a c t

Article history: Received 13 November 2013 Received in revised form 24 June 2014 Accepted 5 July 2014

This study presents the development of a finite element formulation based on positions to deal with steel structures under high temperatures. The nonlinear thermal analysis is done regarding cross sections using two-dimensional flat elements while the three dimensional structural members are modeled by physical and geometrical nonlinear solid-like frame elements of any order. The thermal analysis consists in achieving, along time, the temperature profile at frame elements cross sections that affects the constitutive relations and, consequently, residual stresses at the structural mechanical equilibrium. The proposed element presents a geometrically exact description for large displacements and rotations by means of translations and generalized vectors nodal parameters. The adopted three dimensional Thermo-elastoplastic constitutive relation is based on the von Mises yielding surface for which the temperature dependence is established and an alternative flow rule is developed. Several examples are presented in order to validate the formulation and to show its good results. & 2014 Elsevier B.V. All rights reserved.

Keywords: 3D Frame FEM based on positions Physical and geometrical non-linear analysis Structures at high temperatures Thermo-plasticity Fire

1. Introduction It is well known that the mechanical properties of steel degenerate significantly at elevated temperatures [1]. For example, at 700 1C the strength of steel is only 23% of its ambient temperature strength. At 800 1C the strength reduces to 11% and, at 900 1C, to 6% of its original value. Therefore, structural designers should be able to consider the high temperature effects in the analysis of structures, especially in case of steel structures in fire situation. The simplest way to design steel structures in fire situation is based on the verification of isolated structural members. This verification follows existent standards as described in [2–5], among others. However, it is well known that the structural element behavior is very dependent on its interaction with other elements, mainly in fire [5,6], and that isolated analysis may lead both to super estimated or sub estimated strengths [7]. This affirmation is corroborated by the conclusions of the Cardington full-scale fire tests [8]. Full scale tests are very expensive and complex, therefore the development of computational codes or packages to accurately describe the steel (or other material) structure behavior at high temperatures (fire situations) are of great importance.

n

Corresponding author. Tel.: þ 55 16 33739451; fax: þ55 16 33739482. E-mail address: [email protected] (H.B. Coda).

http://dx.doi.org/10.1016/j.finel.2014.07.005 0168-874X/& 2014 Elsevier B.V. All rights reserved.

There are some computational codes that can model structures under high temperatures. These codes (or packages) can be classified as generalists and specialists, most of them are based on the Finite Element Method. Among the generalists one may mention the ANSYS, ABAQUS and DIANA. Among the specialists one may mention SAFIR [9,10], VULCAN [11,12], ADAPTIC [13–15] and FEMFAN [16–18]. Most of existing codes when using frame elements adopts one dimensional stress strain relations together with simple frame elements based on Euler–Bernoulli kinematics. In order to reduce the computational effort, some researchers have used the plastic hinge concept instead of distributed plasticity [19–23]. Lien et al. [24] proposed the vector form intrinsic finite element method to analyze the behavior of steel frames in fire. In general they consider constant or linear temperature distribution over cross sections, i.e., few of them develop thermal analysis to determine the temperature profile over the cross section in order to make the thermo-mechanical coupling, as presented in this study. In this work we are interested to apply a more recent kinematic description, the generalized vector frame element Coda [25] and Coda and Paccola [26,27] (a solid-like formulation) to model three dimensional steel frames under fire. The interesting features of this formulation are the use of three-dimensional generalized constitutive model considering associative or non-associative plasticity and the generalized vector parameterization that avoids the use of large rotation schemes.

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

The proposed three dimensional elastic–plastic model for steel under high temperatures is achieved here transforming the one dimensional diagrams to a set of multi-linear stress–strain curves at the equivalent von-Mises space. When a change of temperature occurs a new stress–strain diagram necessarily takes place and the thermodynamic consistent internal force changes in order to keep the stress level on the “moving” diagram. The change of the internal force is achieved here directly by the nullity of the plastic objective yielding potential. Due to the internal force change, for the same load level, the iterative solution process is started determining a new equilibrium position with a consequent plastic strain evolution. The thermal analysis is done by means of a 2D finite element approach considering the nonlinear thermal behavior of materials associated with Neumann and Dirichlet conditions, including convection and radiation. The thermal formulation is nonlinear, i. e., considers the thermal properties as functions of temperature and is based on the works of Huang and Usmani [28] and Lewis et al. [29]. The 2D discretization of the thermal analysis is coincident with the discretization of the solid-like frame element cross section, making possible the direct transfer of the current temperature profile to the mechanical analysis. The engineering standard CEN EN 1993-1-2:2005 [30] is used to relate the material properties and temperature. The proposed thermo-mechanical model allows the use of any cross section, open or not, and introduces the consistent warping into solution. The nonlinear procedure is solved by the Newton– Raphson implicit technique. Several examples are shown in order to demonstrate the capability of the proposed formulation in solving thermomechanical elastoplastic problems, particularly, 3D steel frames under high temperatures.

X2

69

H2

P3

m

P2

P2

P4

H1 m

P1 P1 2

(-1,1)

(1,1)

(-1,1)

(1,1)

1

X1 Fig. 1. Two-dimensional solid element: (a) Classic mapping; (b) Vector Mapping.

Mapping

3 2

V21

1

V22

LM m

m

m

X21, X 22 , X 23

V11

m

m

V12

m

X11 , X 12 , X 13 Fig. 2. Tridimensional vector mapping scheme.

1 φ2 ðξ1 Þ ¼ ð1 þ ξ1 Þ 2 2. The solid-like finite element kinematics

The vector mapping is described by Eq. (3) where (X m ℓi ) are nodal (ℓ) coordinates (i) of the reference line (m) and the vectors (H ℓi ) are given by (4).

2.1. Solid-like mapping The purpose of this section is providing the basis of the adopted formulation for the solid-like finite element. This strategy has been proposed by Coda [25] for elastic applications and consists on describing solids by means of generalized unconstrained vectors. Differently from classical formulations, the continuous mapping is done through non-unitary vectors that do not keep the orthogonality after deformation, which implies the natural consideration of shear effects on frame elements and the possibility of distorting cross sections. Complete three-dimensional stress–strain relation is employed and, to avoid locking and improve the structure representation, warping and strain with linear variation along the transverse direction are considered. The resulting mapping represents a long solid and uses only vector variables in Euclidean space and, consequently, the values of stress, strain and position are objective measures, i.e., independent of the translations, rotations and equilibrium path. Fig. 1 shows a two-dimensional quadrilateral solid mapped from the dimensionless space to its shape, following both 2D classical solid procedure and the proposed vector procedure. To illustrate the generalized vector idea we make use of a linear flat membrane and then generalize it to comprise any order frame element. The adopted linear shape functions are given by Eq. (1) and (2), with dimensionless coordinates ξ1 and ξ2 A ½  1; 1. 1 φ1 ðξ1 Þ ¼ ð1  ξ1 Þ 2

ð2Þ

ð1Þ

m xi ðξ1 ; ξ2 ; X m ℓi ; H ℓi Þ ¼ φℓ ðξ1 ÞX ℓi þ

H ℓi φℓ ðξ1 Þ ξ2 2

for i ¼ 1; 2 and ℓ ¼ 1; 2 ð3Þ

Xm 1i ¼

X 1i þX 4i X 2i þX 3i X  X 1i X X 2i ; Xm ; H 1i ¼ 4i and H 2i ¼ 3i 2i ¼ 2 2 2 2 ð4Þ

Through the analysis of Eq. (4) it is evident that the vector mapping is done with the use of non-unitary vectors (H ℓi ) as parameters, which are non-orthogonal to the “reference line” that m connect points P m 1 and P 2 (see Fig. 1(b)). For frame elements, we assume constant height (h0) and change Eq. (4) to the mapping described by Eq. (5), where V ℓi is called a generalized vector. m xi ðξ1 ; ξ2 ; X m ℓi ; V ℓi Þ ¼ φℓ ðξ1 ÞX ℓi þ

h0 ξ V φ ðξ Þ where i ¼ 1; 2 and ℓ ¼ 1; 2 2 2 iℓ ℓ 1

ð5Þ

For simplicity, at the initial configuration the generalized vectors are adopted as unitary and orthogonal to the reference line. However, for the current configuration these vectors are unknown and are neither unitary nor orthogonal to the reference line. Finally, Eq. (6) presents the three-dimensional version of Eq. (5), whose mapping is illustrated by Fig. 2 for a two node

70

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

Fig. 4. (a) Initial configuration; and (b) current configuration.

The pair of vectors V 1i and V 2i in Fig. 4 represents an orthonormal basis for the initial cross section configuration, while the pair of vectors g 1i and g 2i is the generic base for the current cross section configuration. For this reason, Eqs. (7) and (8) are rewritten as:

Fig. 3. Mapping of initial and current configurations.

element. m xi ðξ1 ; ξ2 ; ξ3 ; X m ℓi ; V ℓi Þ ¼ φℓ ðξ1 ÞX ℓi þ

1 h0

2

ξ2 V 1iℓ φℓ ðξ1 Þ þ

2 h0

2

ξ3 V 2iℓ φℓ ðξ1 Þ

ð6Þ

for i ¼ 1; 2; 3 and ℓ ¼ 1; 2.

Eq. (6) represents the initial mapping of the analyzed frame element. To describe the current position one changes coordinates x by coordinates y and vectors V 1ℓi and V 2ℓi by the so called generalized vectors g 1ℓi and g 2ℓi . Fig. 3 illustrates a typical mapping from initial to current configuration for a curved solid-like bar element with cubic approximation (in the longitudinal direction) and flat cross section. Although this figure shows a schematic cubic element, the developed and applied element allows the choice of any approximation order. In Fig. 3, f and A represent the Lagrangian description of the change of configuration (deformation) and its gradient, respec0 1 tively. Similarly, ðf ; f Þ and ðA0 ; A1 Þ are the mappings and their gradients of initial (0) and current (1) configurations from the auxiliary dimensionless space, written as: 1

0

1 fi

¼

2

h0 h ξ V 1 φ ðξ Þ þ 0 ξ3 V 2iℓ φℓ ðξ1 Þ 2 2 iℓ ℓ 1 2

1 2 h0 h ξ2 g 1iℓ φℓ ðξ1 Þ þ 0 ξ3 g 2iℓ φℓ ðξ1 Þ yi ¼ φℓ ðξ1 ÞY m ℓi þ

2

ð9Þ

1

ð10Þ

j j 1 j j 2 f i ¼ yi ¼ φℓ ðξ1 ÞY m ℓi þ α ψ ðγ 1 ; γ 2 Þg iℓ φℓ ðξ1 Þ þ β ψ ðγ 1 ; γ 2 Þg iℓ φℓ ðξ1 Þ

in which γ1 and γ2 are ordinary dimensionless coordinates of triangular finite element, αj and βj (known values) are coordinates of the points of the auxiliary mesh of the cross section and ψ j ðγ 1 ; γ 2 Þ are the triangular element shape functions.

2.2. Deformation mapping

f i ¼ xi ¼ φℓ ðξ1 ÞX m ℓi þ

0

1 2 j j j j f i ¼ xi ¼ φℓ ðξ1 ÞX m ℓi þ α ψ ðγ 1 ; γ 2 ÞV iℓ φℓ ðξ1 Þ þ β ψ ðγ 1 ; γ 2 ÞV iℓ φℓ ðξ1 Þ

2

ð7Þ

2.3. Modes of warping and linear strain variation for 3D stress–strain consideration By adopting a complete three-dimensional constitutive law and Reissner kinematics a volume locking appears. This is due to the influence of Poisson's ratio and unbalanced approaches to normal and shear strains [25]. To correct this problem, the current configuration kinematics is enhanced by the introduction of linear strain behavior along transverse directions and the consideration of warping effects. From this reasoning Eq. (10) is rewritten as: 1

j j 1 j j 2 f i ¼ yi ¼ φℓ ðξ1 ÞY m ℓi þ α ψ ðγ 1 ; γ 2 Þg iℓ φℓ ðξ1 Þ þ β ψ ðγ 1 ; γ 2 Þg iℓ φℓ ðξ1 Þ

þ f½αj ψ j ðγ 1 ; γ 2 Þ2 ðΛ1ℓ φℓ ðξ1 ÞÞgg 1iℓ φℓ ðξ1 Þ þ f½βj ψ j ðγ 1 ; γ 2 Þ2 ðΛ2ℓ φℓ ðξ1 ÞÞgg 2iℓ φℓ ðξ1 Þ 1 2 þ ½Del  ψ  ðγ 1 ; γ 2 Þζ ijk ½g jℓ φℓ ðξ1 Þ½g kℓ φℓ ðξ1 Þφℓ ðξ1 ÞW ℓ

Λ1ℓ

ð8Þ

In Eqs. (7) and (8) the lower index ℓ is no longer limited to two nodes and i varies from 1 to 3. As mentioned before, vectors g 1i , g 2i and the reference line generate the bar element. However, the geometry of the cross section can be improved to consider different shapes and materi1 2 als. In this sense it is necessary to replace the constants h0 and h0 by a discretization which, multiplied by the corresponding base vectors V 1i and V 2i , in the initial configuration, results in the desired cross section geometry. This procedure comprises a cross section deformation that depends only on the transformations V 1i -g 1i and V 2i -g 2i . To do this discretization, we adopt a two-dimensional mapping with cubic isoparametric finite elements, as sketched in Fig. 4. It is important noting that this discretization is also used to solve the nonlinear thermal problem to find the temperature profile of elements for fire analysis, see Section 6.

ð11Þ

Λ2ℓ

and are the intensity of the transverse strain rate in which and ½αj ψ j ðγ 1 ; γ 2 Þ2 and ½βj ψ j ðγ 1 ; γ 2 Þ2 are the associated enhancement modes. The last term in Eq. (11) represents warping, where Del  are the modal warping parameters related to the  points of each triangular finite element (superscript el) of the cross section discretization. This modal warping is solved for a unitary rotation by a preprocessor described by [26]. The warping direction in each solid-like finite element is given by the vector product between g 1i and g 2i given by ni ðξ1 Þ ¼ ζ ijk g 1j ðξ1 Þg 2k ðξ1 Þ ¼ ζ ijk ½g 1jℓ φℓ ðξ1 Þ½g 2kℓ φℓ ðξ1 Þ

ð12Þ

in which ζ ijk is the well-known Levi–Cevita permutation tensor. The warping intensity W ℓ is the new degree of freedom (unknown) approximated along the finite element as the other variables. Eq. (11) represents an objective vector mapping from the dimensionless space to the current configuration of the body. This mapping results in 12 degrees of freedom per node that are: three translations, six vector components, two intensity values of the strain rate in the transverse directions and the warping intensity.

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

71

2.4. The deformation gradient

of temperature and the nominal thermal dilatation coefficient as:

! It is not necessary to explicitly know f to achieve the deformation gradient of the total change of configuration function, !0 !1 but to know f given by Eq. (9) and f given as a trial by Eq. (11). From these two mappings one writes the following gradients

αG ðθÞ ¼ αðθÞ þ αðθÞ2 ðθ  θref Þ=2

0

A0ij ¼

1

∂f i ∂f and A1ij ¼ i ∂ξj ∂ξj

ð13Þ

in which γ1 and γ2 are called ξ2 and ξ3 respectively. From these values the deformation gradient given by [31,32] holds as, A ¼ A1 UðA0 Þ  1 or Aij ¼ A1ik B0kj

ð14Þ

where B0 ¼ ðA0 Þ  1 . From the deformation gradient one writes the Cauchy–Green stretch and the Green strain as: C ¼ At UA ¼ ðA0 Þ  t UAt1 U A1 U ðA0 Þ  1 or C ij ¼ B0γi A1βγ A1βα B0αj

ð15Þ

1 1 E ¼ ½C  I or Eij ¼ ½C ij δij  2 2

ð16Þ

in which θref is a reference temperature value. It is important to note that when plastic applications take place the Green strain given in Eq. (18) is considered as the elastic part of the total train, i.e., E e . As the Green strain has been written as a function of nodal positions, the strain energy stored in the structure is written as: Z ! Ueð Y Þ ¼ ue dV 0 ð21Þ V0

where V0 is the initial volume of the analyzed structure [25] and ! Y is a general nodal position vector including all 12 degrees of freedom per node. 3.2. Total potential energy In order to write the total potential energy of the mechanical ! system Πð Y Þ one sums the strain energy, the potential energy of concentrated and distributed forces, resulting: Z Z ! ! ! ! ! !m ! ue ð Y ÞdV 0  F U Y  q U y ð Y ÞdS0 ð22Þ Πð Y Þ ¼ V0

3. Elastic procedure In order to give a brief, but complete, description of the developed thermo-mechanical solution, this section shows the basic steps of the elastic geometrical nonlinear solution procedure. First of all it is necessary to detach that the elastic constitutive relation is dependent of temperature; however this dependence is quite simple as the temperature profile is introduced as an independent variable (related to time) to be achieved by the thermal part of the method. Moreover, we consider negligible the heat generation resulting from the mechanical dissipation for quasi-static fire analysis. For applications that this heat generation is important (as high speed impact) the reader is invited to consult the work of Carrazedo and Coda [33] and others cited therein. In this section, the elastic Saint-Venant–Kirchhoff constitutive relation is adopted to relate the second Piola-Kirchhoff stress and the Green strain for frame elements, including thermal dilatation strains. Moreover the nonlinear equilibrium equation is achieved by the principle of stationary energy. Plasticity is introduced in Section 4 and the thermo-mechanical procedure in Sections 6 and 7. 3.1. Frame element strain energy:

ð20Þ

S0

! ! in which F is the external nodal force vector and q is the general distributed force vector written as function of nodal values, given by: qi ¼ ϕℓ ðξÞQ ℓi

ð23Þ

!m ! In Eq. (22), y ð Y Þ is the current position at the reference line ! of frame elements, written as function of nodal positions Y (see Eqs. (7) and (8)) and dS0 is the infinitesimal length of the curved frame element. In order to find the equilibrium configuration one applies the principium of stationary potential energy on Eq. (24), using as parameters the nodal positions, as: Z δΠ ¼

! ! ! ∂ue ∂E : !dV 0 Uδ Y  F Uδ Y  V 0 ∂E ∂Y

Z

!m ! ! ∂y q U ! dS0 U δ Y ¼ 0 S0 ∂Y

ð24Þ

The first term of the first integral of Eq. (24) is the second PiolaKirchhoff stress. Therefore, in order to simplify the understanding of the physical nonlinear procedure shown in next section, one rewrites Eq. (24) as: Z

! ! ! ∂E S : !dV 0 U δ Y  F U δ Y  V0 ∂Y

Z

!m ! ! ∂y q U ! dS0 Uδ Y ¼ 0 S0 ∂Y

As mentioned before the specific strain energy to be adopted here is the Saint-Venant–Kirchhoff one, that relates the Green strain and the second Piola-Kirchhoff stress in a linear way,

δΠ ¼

1 ue ¼ E:ℭ:E 2

The understanding of Eq. (25) can be further improved defining the first integral as the internal nodal force and the last integral as the equivalent nodal force of distributed forces [34], resulting:

in which E is the constitutive tensor. mechanical analysis mechanical (m) and strains, results:

ð17Þ Green strain tensor and ℭ is the elastic It is worth to mention that for thermothe Green strain tensor is constituted by a a thermal (th) part, that, for relative small

E ¼ E m  E th

ð18Þ

ð19Þ

in which I is the identity tensor of second order and αG is the thermal dilatation coefficient for Green strain, given as a function

ð26Þ

in which L is the matrix that transforms the distributed forces into ! nodal equivalent ones. Due to the arbitrariness of δ Y Eq. (26) results into the geometrical nonlinear equilibrium equation, as: ! !int !int ! F þ M  F LU Q ¼ 0

For a given temperature (θ) one writes: E th ¼ αG ðθÞðθ θref ÞI

! ! !int ! ! ! δΠ ¼ F Uδ Y  F Uδ Y  L U Q U δ Y ¼ 0

ð25Þ

ð27Þ

The Newton–Raphson procedure is used here to solve the nonlinear equilibrium, which is described in Section 5, including the physical nonlinear behavior of materials.

72

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

in which Ee and Ep are, respectively, the elastic and plastic parts of E. Moreover Ee comprise the thermal strain as described by Eq. (19). Therefore, the free energy potential can be written as:

4. Inelastic procedure 4.1. Equilibrium equation The difference between the elastic and inelastic procedures is the way the potential energy of external forces is transferred to the deformed body and dissipated during the loading process. In this case, instead of using the elastic strain energy, one uses the Helmholtz free energy potential, that usually does not have an explicit expression; however, its variation is known [35]. Following this reasoning the total potential energy is written as: ! Πð Y ; θÞ ¼

Z V0

! ! ! Ψ ðEðθ; Y Þ; ℭðθÞ; ℌðθÞ; aÞdV 0  F U Y 

Z S0

! !m ! q U y ð Y ÞdS0

ð28Þ where Ψ is the free energy potential, written as function of the ! Green strain tensor Eðθ; Y Þ, the elastic constitutive tensor ℭðθÞ, the hardening parameter tensor ℌðθÞ and the thermodynamic parameter a. It is important to note that this free energy potential is similar to the isothermal one as we decide that the heat generated by the mechanical dissipation is negligible in the analyzed problem. Therefore, the dependence of temperature is omitted in the passages where it does not have importance and opportunely resumed. The variation of the total potential energy is null at the equilibrium position, i.e.: Z δΠ ¼

! ! ! ∂Ψ ∂E : !dV 0 Uδ Y  F Uδ Y  ∂E V0 ∂Y

Z

!m ! ! ∂y q U ! dS0 U δ Y ¼ 0 S0 ∂Y

ð29Þ

in which parameter a is not present due to its intrinsic relation with E and S to be shown in item 4.2. Even for inelastic problems, the energy conjugate of the free energy potential regarding Green strain is the second PiolaKirchhoff stress tensor. In this way, Eq. (29) can be written exactly as Eq. (25). However, with a physical nonlinear meaning, i.e., while the passage from Eqs. (24) and (25) is done by a simple differentiation of the quadratic potential shown by Eq. (17), now it is necessary the definition of inelastic (plastic for instance) constitutive relations to describe the material behavior and its evolution rule. Moreover, this elastoplastic relation should be adapted to comprise imposed temperature changes. 4.2. Elastoplastic constitutive relation In this section, the constitutive elastoplastic relation is described. One may classify the developed constitutive relation as non-associative with an unknown plastic potential for which its normal direction (plastic flow) is defined by the trial stress direction. The theoretical basis of this constitutive relation is given by the study of volume change in plastic phase given by Botta et al. [36]. Moreover, we assume a multi-linear stress strain diagram (in the von-Mises equivalent space) without a limit for the number of elastoplastic patches, which enables an accurate reproduction of any stress strain relation. 4.2.1. Basic definitions Although, the developed displacements are high, the strain level present in our applications is relatively small. Therefore, the Green strain approximates the linear strain and the second PiolaKirchhoff stress can be used in place of the Cauchy stress. Following this reasoning we adopt the additive strain decomposition, as: E ¼ Ee þ Ep

ð30Þ

1 1 Ψ ðE  E p ; ℭ; ℌ; aÞ ¼ E e : ℭ : E e þ a:H:a 2 2

ð31Þ

where ℭ is the elastic constitutive tensor and ℌ (in a general representation ) is a second order tensor that includes the hardening parameters. The temperature dependence is omitted as it is considered an independent variable. The second Piola–Kirchhoff stress and the thermodynamic force χ are written as: S¼

∂Ψ ¼ ℭ : E e ¼ ℭ : ðE  E p Þ ∂E

χ¼

∂Ψ ¼  ℌ:a ∂a

ð32Þ ð33Þ

At this point it is important to mention that χ as a function of ℌ is dependent of temperature and this dependence is solved in Section 7. As mentioned after Eq. (29), in order to eliminate a from equilibrium equation it is necessary to relate Ep and a. In classical formulations this is done introducing a plastic potential FðS; χÞ [37] from which one writes: ∂F p E_ ¼ λ_ ∂S

ð34Þ

∂F a_ ¼ λ_ ∂χ

ð35Þ

p where λ is the plastic multiplier. From Eqs. (34) and (35) E_ and a_ are related when one eliminates the plastic multiplier that satisfies the Khun–Tucker conditions:

λ_ Z0; f r 0; λ_ f ¼ 0

ð36Þ

where f ðS; χÞ is the yielding (failure) surface. Particularly, in this work we adopt the Von-Mises surface, i.e., pffiffiffiffi f ¼ J 2  χ Sp ð37Þ in which the equivalent von-Mises stress t

σ ¼ J 2 ¼ 1=2½S : S 

ð38Þ

is the second invariant of the deviatory stress S, and pffiffiffi pffiffiffi Sp ¼ σ p = 3 ¼ f p = 3

ð39Þ

is the initial size of the Von-Mises surface, in which σ p is the yielding uniaxial stress. In expression (37) χ is assumed to be scalar resulting in scalar values for α and ℌ in Eqs. (33) and (35). Moreover, isotropic hardening is adopted. Although the equivalent strain is not used at this stage of the solution, see Section 7, it is worth mentioning that it is calculated following a similar way as the equivalent stress, i.e., t

ε ¼ 1=2½E:E 

ð40Þ

Although not mentioned σ p and Sp depends on temperature. Therefore the yielding surface size and evolution depend upon temperature, this dependence will be solved in Section 7. Instead of introducing a plastic potential which gradient imposes the plastic flow direction, see Eq. (34), we prefer to define the plastic flow from the current stress state (trial in practical situations) as indicated by [36], that is: ∂F p ¼ η ¼ D :S ∂S

)

p _ E_ ¼ λη

ð41Þ

and ∂F ¼ 1 ∂χ

)

λ_ ¼  α_

ð42Þ

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

where G p p D ¼ pffiffiffiffiðℭ Þ  1 J2

ð43Þ

p

73

Introducing (50) and (52) into Eq. (47) and operating ℭ : ΔE P one writes:   k ð53Þ S ¼ t S ðn þ 1Þ  Δλ0 t S ðn þ 1Þ þ P trðt S ðn þ 1Þ ÞI 3k

ℭ is a fourth order constitutive tensor similar to a typical p elastic constitutive tensor and D is related with its inverse. In order to write η in a concisely form one appeals to some basic elastic operations, that is:

with

S ¼ ℭ : E ¼ 2GE þ k trðEÞI

Splitting t S ðn þ 1Þ in its deviatory and hydrostatic parts, Eq. (53) is rewritten as:   S ðn þ 1Þ Þ k S ¼ t S ðn þ 1Þ ð1  Δλ0 Þ þ trðt 1  P Δλ0 I ð55Þ 3 k

E ¼ D:S ¼

ð44Þ

S trðSÞ þ I 2G 9k

ð45Þ

where k is the bulk modulus and G is the shear modulus. Therefore, ! G S trðSÞ p þ pI η ¼ D :S ¼ pffiffiffiffi ð46Þ J 2 2G 9k p

p

In expression (46) k ¼ E=½3ð1  2νp Þ means that ℭ is similar to C less the bulk modulus that depends on the plastic Poisson. Following this strategy, when νp ffi 0:5 isochoric plasticity takes place and when ν ¼ νp the plastic flow occurs in the same direction p of elastic flow. In this work we adopt ℭ ¼ ℭ, defined by the second derivative of expression (17) regarding strain. This choice releases any possible locking due to the adopted Reissner kinematic. 4.2.2. Numerical evolution of variables Instead of following classical procedure, for which the previous equations are solved in an infinitesimal way, from basic definitions and considering the hardening parameter H constant by parts (multi-linear plasticity) we write an implicit (numerical) closed solution for the constitutive model variables. We assume a typical interval ½t n ; t n þ 1  and adopt an elastic trial, written as: t

E ðn þ 1Þ -Total stress elastic trial

t

E pðn þ 1Þ

t

λðn þ 1Þ ¼ λn ¼  an -Internal variable trial

t

χ ðn þ 1Þ ¼ χ n -Thermodynamic internal force trial

¼ E pn -Accumulated

ð47Þ

S ðn þ 1Þ S ðn þ 1Þ Þ þ trðt p IÞ 2G 9k

ð49Þ

ΔE ¼ Δλη

ð50Þ

E pðn þ 1Þ ¼ t E pðn þ 1Þ þ ΔE p

ð51Þ

χ ðn þ 1Þ ¼ χ ðn þ 1Þ þ ℌΔλ

or Δλ ¼

t

ð52Þ

f ðn þ 1Þ ðG þ ℌÞ

ð58Þ

Finally, the searched plastic strain is found using (48), (50) and (58), that is, 0 1 ΔE P ¼

t

f ðn þ 1Þ Bt S ðn þ 1Þ trðt S Þ C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffiA @ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P t ðG þ ℌÞ 2 t J þ 9k J 2ðn þ 1Þ

ð59Þ

2ðn þ 1Þ

To find the Hessian matrix for the solution technique one uses the differential process. The equivalent expressions to (58) and (59) are: S : ℭ : dE pffiffiffiffi S : ℭ : η þ 2 J2 ℌ

dE ¼

t

t

ð57Þ

p

If f ðn þ 1Þ r 0 the step is elastic and the trial variables are the t correct ones. However, if f ðn þ 1Þ 4 0 then the step is elastoplastic p and ΔE must be fount to guaranty the nullity of f ðn þ 1Þ . This situation is solved here in a closed form, because ℌ is considered constant by pats. Therefore, from Eq. (46) one writes:

p

From Eqs. (56) and (54) one finds: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ðG þ ℌÞΔλ þ t J 2ðn þ 1Þ  t χ ðn þ 1Þ  Sp ¼ 0

and

where the trial stress t S is known. Using this value we calculate the trial yielding potential: qffiffiffiffiffiffi t f ðn þ 1Þ ¼ t J 2  t χ  Sp ð48Þ

t

ð56Þ

plastic strainðassumed as trialÞ

S ðn þ 1Þ ¼ ℭ:ðt E ðn þ 1Þ  t E pn  ΔE p Þ ¼ t S ðn þ 1Þ  ℭ:ΔE p

ð54Þ

Calculating f ðn þ 1Þ and imposing its nullity, results: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f ðn þ 1Þ ¼ ð1 Δλ0 Þ t J 2ðn þ 1Þ  t χ ðn þ 1Þ  ℌΔλ Sp ¼ 0

dλ ¼

Using these variables we calculate the stress level considering ΔE p as the main unknown, as follows:

G η ¼ Dp : t S ðn þ 1Þ ¼ pffiffiffiffi J2

G Δλ0 ¼ Δλqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t J 2ðn þ 1Þ

ð60Þ

! S : ℭ : dE pffiffiffiffi η S : ℭ : η þ 2 J2 ℌ

ð61Þ

Finally, the tangent constitutive tensor is: ! ∂S ∂2 Ψ ℭ:ηS:ℭ ep pffiffiffiffi ¼ ¼ ℭ ¼ ℭ  ℭp ℭ ¼ ∂E ∂E∂E S : ℭ : ηþ 2 J 2 ℌ

ð62Þ

This tensor is not symmetric, as the plastic potential (not explicitly written) is not the same as the yielding (failure) surface. In order to avoid asymmetry, in the numerical procedure, we use only the diagonal of ℭp to assemble the Hessian. This quasitangent matrix accelerates the convergence of the solution process, see Section 5.

5. Newton–Raphson procedure Knowing the elastoplastic constitutive model, to start the solution process one rewrites the equilibrium Eq. (29) removing ! the arbitrary variation of positions (δ Y ), as: !! g ðY Þ¼

Z

! ∂Ψ ∂E : !dV 0  F  ∂E V0 ∂Y

Z

!m ! ! ∂y q U ! dS0 ¼ 0 S0 ∂Y

in which ∂Ψ =∂E is the second Piola–Kirchhoff stress, Eq. (32).

ð63Þ

74

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

Remembering that the process is nonlinear, Eq. (63) is ! expanded in Taylor series from a trial position Y 0 , i.e.:  ! ! ! ∂ g  !! !! ΔY ¼ 0 ð64Þ g ð Y Þ ffi g ð Y 0 Þ þ ! !  ∂Y Y 0 Solving the linear system of Eq. (64) one finds the position ! ! ! ! ! correction, applied as Y ¼ Y þΔ Y until jΔ Y j=j X jo tol, in which ! tol is the tolerance in positions and X is the initial nodal positions of the body. As the applied forces are conservative the derivative of ! g regarding positions results: ! ! ∂g ∂E ∂2 Ψ ∂E ∂2 E ¼ : : þ S : ð65Þ ! ! ! ! ! ∂Y ∂ Y ∂E∂E ∂ Y ∂Y ∂Y in which ∂2 Ψ =∂E∂E is given by Eq. (62). The derivatives of the Green strain regarding positions can be seen elsewhere [25–27].

For a given time interval Δt and observing the scheme of Fig. 5 one writes

6. Thermal analysis As mentioned in introduction the main objective of this study is the mechanical behavior of steel structures under high temperatures (fire). Therefore, this section briefly describes the transient thermal analysis used to generate the temperature profiles at the cross sections of structural members. Most of equations, developments and implementations are based on Huang and Usmani [28] and Lewis et al. [29]. In the next section, knowing the temperature profile, the elastoplastic constitutive relation is adapted to comprise temperature changes. In this study, thermal analyses are based on the classical heat diffusion (or conduction) equation in solids, ρc

dθ λ∇2 θ  q ¼ 0 dt

∂θn þ ϕ θn þ 1  θn  ∂t Δt

In Eq. (67) [C] is the heat capacity matrix, [K] is the thermal stiffness matrix and {F} is the vector of nodal heat flux. The matrices [C], [K] and vector {F} are known while {θ} is the vector of unknown parameters (temperature). The nonlinear characteristic of Eq. (67) is the dependence on temperature of matrices [C] and [K] which becomes evident in the solution process. To obtain the numerical solution of Eq. (67) it is necessary to discretize the differential operator involving the transient term ∂θ=∂t. Here it is done using a numerical approximation based on the Finite Difference Method, with a strategy similar to that presented in [28,29]. Fig. 5 illustrates a typical temperature variation during a time interval [tn;tn þ 1]. Using Taylor series, it is possible to describe the variation of temperature at time n þ ϕ as ∂θn þ ϕ Δt ∂ θn þ ϕ þφ þ::: ∂t 2 ∂t 2 2

2

ð68Þ

ð69Þ

The temperature at the instant t n þ ϕ can be written as θn þ ϕ ¼ ϕθn þ 1 þð1 ϕÞθn Writing Eq. (67) for instant t n þ ϕ ,   ∂θ þ ½Kn þ ϕ fθgn þ ϕ ¼ fFgn þ ϕ ½Cn þ ϕ ∂t n þ ϕ

ð70Þ

ð71Þ

and using approximations (69) and (70) results   ½Cn þ ϕ ½Cn þ ϕ þ ϕ½Kn þ ϕ fθgn þ 1 ¼ þ ðϕ  1Þ½Kn þ ϕ fθgn þfFgn þ ϕ Δt Δt ð72Þ

ð66Þ

In which c is the specific heat, ρ is the density, λ is the thermal conductivity, q is the amount heat generated in the material per unit volume, θ is temperature (Kelvin) and t is time. Considering the discretization of a domain (in this case a cross section at the reference configuration) by finite elements (Fig. 4(a)) and applying the method of weighted residuals with consideration of the appropriate boundary conditions, one achieves Eq. (67) which describes, in matrix form, the balance of heat flow inside the considered domain at a given instant.   ∂θ ½C þ ½Kfθg ¼ fF g ð67Þ ∂t

θn þ ϕ ¼ θn þϕΔt

Fig. 5. Temperature variation in a time step.

The system of Eq. (72) is compacted according to Eq. (73) where the terms ½K~ n þ ϕ and fF~ gn þ ϕ are given by Eqs. (74) and (75), respectively. ½K~ n þ ϕ fθgn þ 1 ¼ fF~ gn þ ϕ ½K~ n þ ϕ ¼ fF~ gn þ ϕ ¼

½Cn þ ϕ þ φ½Kn þ ϕ Δt 

½Cn þ ϕ þðϕ  1Þ½Kn þ ϕ fθgn þ fFgn þ ϕ Δt

ð73Þ ð74Þ ð75Þ

In order to solve Eq. (73) one needs to choose the value of ϕ, which defines the integration scheme. According Cook et al. [38], for values of ϕ Z 0:5 the solution of linear or nonlinear problems are unconditionally stable and when it is assumed to be 2/3 the Galerkin scheme is followed. At elevated temperatures one must consider the thermal properties of materials as a function of temperature, such as specific heat, density and thermal conductivity. Furthermore, the thermal flow at boundary of structural members is highly dependent of the typical actions of fire. This aspect gives to the thermal problem a nonlinear characteristic, detached in the temperature dependence of ½Kn þ ϕ , ½Cn þ ϕ and fFgn þ ϕ , in Eq. (72). It is obvious that depending on temperature, these values will also depend on time. To solve this nonlinear system of equation the Newton– Raphson procedure is employed as follows. Eq. (73) is rewritten considering a trial temperature, as fψ gin þ 1 ¼ fF~ gin þ ϕ  ½K~ in þ ϕ fθgiðtrialÞ nþ1 fψgin þ 1

ð76Þ

is the residuum of the solution for the trial in which temperature. Eq. (76) is expanded in Taylor series up to the first

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

75

order imposing the nullity of the final residuum, i.e., fψginþþ11 ¼ fψgin þ 1 þ ∇fψgin þ 1 fΔθgin þ 1 ffi0

ð77Þ

Solving the linear system of Eq. (77) a correction for temperature is achieved, fΔθgin þ 1 ¼ ð∇fψgin þ 1 Þ  1 fψgin þ 1

ð78Þ

improving the trial temperature by i fθginþþ1ðtrialÞ ¼ fθgiðtrialÞ 1 n þ 1 þ fΔθgn þ 1

ð79Þ

For a time step this procedure is repeated until fΔθgin þ 1 can be considered negligible.

fψgin þ 1

or Fig. 7. Stress–strain relation for carbon steel at elevated temperatures (CEN EN 1993-1-2:2005 [30]).

7. Thermo-mechanical coupling Section 4 describes the elastoplastic constitutive relation, Section 5 indicates how the nonlinear mechanical problem is solved and Section 6 briefly describes how the temperature profile is achieved in a bar element cross section. In this item, considering known the temperature profile, the constitutive relation is extended to include the temperature influence. 7.1. General aspects of stress–strain relations under high temperatures As mentioned before, when temperature changes the mechanical properties of material change, as elastic modulus, hardening parameter and strength levels. Taking as reference the Eurocode CEN EN1993-1-2:2005 [30], the steel reduction factors for yielding strength fy, elastic limit fp and elastic modulus E as a function of temperature can be described by Fig. 6. Following [30] a general diagram of uniaxial stress–strain relation for carbon steel can be seen in Fig. 7. Particularly, as an example, for steel with yielding limit of 250 MPa and strain up to 2% some stress–strain relations achieved using [30] can be seen in Fig. 8. As mentioned before, we assume a multi-linear stress strain diagram without a limit for the number of elastoplastic patches, which enables an accurate reproduction of any stress strain relation. The most important point of the diagram that should be precisely indicated is ðε; f p Þ that defines the beginning of plastic evolution, as it should be in accordance with the elastic modulus definition. 7.2. Thermo-mechanical coupling of the 3D constitutive model First of all, the uniaxial diagrams of Figs. 7 and 8 are transformed to equivalent stress–strain relations (based on Eqs. (39)

pffiffiffi and (40)) simply dividing it by 3. Therefore, at an imposed load level and temperature the structure is at equilibrium and one knows, from Eqs. (39), (40), (51) and (52), the equivalent stress (σðθÞ), the equivalent strain (εðθÞ), the developed equivalent plastic strain (λðθÞ) and the thermodynamic internal force (χðθÞ). The last two values store the history of plastic evolution along the loading process, including previous temperature changes. The equivalent strain already contains the thermal dilatation related to temperature changes. If, for this load level, a temperature change occurs our task is to find the new values σðθ þ ΔθÞ, εðθ þ ΔθÞ, λðθ þ ΔθÞ and χðθ þΔθÞ. As it is a nonlinear process, one should remember that before and after the change of temperature these values should respect the restriction: f ðθÞ ¼ σðθÞ  χðθÞ  Sp ðθÞ r0

ð80Þ

and

1.0

f ðθ þ ΔθÞ ¼ σðθ þ ΔθÞ  χðθ þ ΔθÞ  Sp ðθ þ ΔθÞ r0

0.8

Reduction factor

Fig. 8. Graphic representation of stress–strain relation for steel with f y ¼ 250 MPa and strain up to 2%. using CEN EN 1993-1-2:2005 [30].

ð81Þ

with 0.6

_ ¼0 λ_ Z 0 and λf

ð82Þ

0.4 0.2 0.0 0

200

400

600

800

1000

1200

Temperature (ºC) Fig. 6. Reduction factors for the stress–strain relation of carbon steel at elevated temperatures.

The procedure to be followed is the same as described in Section 4.2.2. (variable evolution) with the difference that, before starting the calculation of Δλ the stress–strain curve is changed and some starting (trial) variables should p beffiffiffi calculated. The new value Sp ðθ þΔθÞ ¼ f p ðθ þΔθÞ= 3 is known from the changing stress–strain diagram, see previous item. The internal force produced by the history of hardening is no longer valid and to achieve a new plausible trial value, for this

76

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

Displacement at the midspan (mm)

0 -5 -10 -15 -20 -25

Yu et al. [12] - ABAQUS - Solid element Yu et al. [12] - VULCAN - Brick element Yu et al. [12] - VULCAN - Beam element ANSYS - Beam 189 element SYSAF - Solid-like element

-30 -35 -40 0

100

200

300

400

500

600

700

Temperature (ºC) Fig. 11. Displacement at the midspan. Fig. 9. Coupling procedure – schematic description.

to consult [39]. For applications that heat generation is important (as high speed impact) the reader is invited to consult [33]. Our code will be called SYSAF all along this section. Only the last two examples employ the nonlinear thermal solver, as there is a lack of references to perform comparisons. 8.1. Simply supported beam under actions of heating and cooling Fig. 10. Simple supported beam geometry [12].

variable, the following auxiliary equation is proposed: σðθ þ ΔθÞ  t χ ðθ þ ΔθÞ  Sp ðθ þ ΔθÞ ¼ 0

ð83Þ

in which σðθ þ ΔθÞ is the equivalent stress level achieved using the equivalent stress–strain diagram (Fig. 9) for ðθ þΔθÞ considering εðθÞ, resulting: t

χ ðθ þ ΔθÞ ¼ σðθ þΔθÞ  Sp ðθ þ ΔθÞ

ð84Þ

Now, using Eqs. (18) and (40), the total equivalent strain is changed to comprise the variation of thermal strain. Using Eq. (47) and the updated strain (elastic part), the trial stress value is changed only in its elastic part, as t

S ðn þ 1Þ ¼ ℭ:ðt E ðn þ 1Þ  t E pn Þ

ð85Þ

resulting in the trial equivalent stress σ ðθ þΔθÞ, using (38). Finally the new evolution equation, the same as Eq. (55), is written as t

f ðn þ 1Þ ¼ ð1  Δλ0 Þt σ ðθ þ ΔθÞðn þ 1Þ  t χ ðθ þ ΔθÞðn þ 1Þ  ℌðθ þΔθÞΔλ Sp ðθ þ ΔθÞ ¼ 0

ð86Þ

and solved, as did previously, finding Δλ and all new variables to be updated in the Newton–Raphson iterative process of the mechanical equilibrium. It is worth noting that the hardening parameter in Eq. (86) is also known directly from the new equivalent stress–strain diagram. Moreover, after solving (86) a plastic strain develops resulting in a mechanical accommodation due to the reduction of the thermodynamic internal force promoted by the change in stress–strain diagram. It is important to note that for cooling the procedure is the same, however in most situations, the internal force grows and the step is elastic. Fig. 9 shows a representation of the steps described in this item.

8. Examples In this section some high temperature examples are used to show the possibilities of the formulation, as well as to demonstrate its good results when compared with other computational codes and references. For validation examples at room temperature comprising inelastic analyses of steel frames the reader is invited

A simply supported beam, Fig. 10, is analyzed under actions of heating and cooling using the CEN EN 1993-1-2:2005 [30] temperature dependence of the constitutive model and relative thermal elongation for steel at high temperatures. The achieved results are compared with the ones presented by Yu et al. [12] using solid elements, via ABAQUS and VULCAN softwares. This beam is heated from 20 1C to 670 1C and then cooled back to 20 1C. The elastic modulus is E ¼ 210 GPa, the reference yielding stress f y ¼ 355 MPa and the adopted Poisson is ν ¼ 0:0. Yu et al. [12] solved this problem using both solid and frame elements. They adopted 1920 hexahedral solid finite elements (16 for a cross section) and three quadratic beam finite elements. In our work we adopt a temperature step of 2.5 1C, four cubic finite elements and eighth triangular cubic “elements” for thermal considerations and elastoplastic treatment of cross sections. Moreover, we solved the same example using the ANSYS package using ten quadratic BEAM189 elements. Results of midspan displacements versus temperature are compared in Fig. 11. As one can see, results achieved by [12] are in reasonable accordance to ABAQUS, however our results compare very well with ANSYS. Next examples compare SYSAF with experimental results and other computational codes confirming the good behavior of our formulation to model steel frame structures under fire. 8.2. Beam subject to uniform temperature test Fig. 12 shows a typical beam tested by Rubert and Schaumann [40] with the objective of verifying the behavior of structures subjected simultaneously to mechanical and thermal actions. Izzuddin et al. [14] present numerical results using their computational code (ADAPTIC) to model a Rubert and Schaumann [40] beam. We did the same using our code (SYSAF). The cross section of the beam is of type IPE 80 and is discretized as depicted in Fig. 12(b). The beam is discretized by four cubic finite elements. Four load levels, proportional to the collapse load at 20 1C, are imposed to the beam. The proportionality constants are 0.2, 0.5, 0.7 and 0.85. At 20 1C the elastic modulus is E ¼ 210 GPa, the reference yielding stress is f y ¼ 399 MPa and the adopted Poisson is ν ¼ 0:0. The thermal dilatation coefficient is constant and equal to 14  10  6 1C  1. The elastic modulus (kE;θ ), the proportionality limit (kp;θ ) and yielding limit (ky;θ ) reduction factors are given by

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

77

F w 114 cm

Fig. 12. (a) Structural scheme; (b) Cross section discretization.

Table 1 Reduction factors from [40]. θa ðo CÞ

kE;θ

ky;θ

kp;θ

20 100 200 300 400 500 600 650 700 800 900 1000

1.0000 1.0000 0.9000 0.8000 0.7000 0.6000 0.5000 0.4500 0.4000 0.3000 0.2000 0.1000

1.0000 1.0000 1.0000 1.0000 1.0000 0.7333 0.4667 0.3333 0.2000 0.1333 0.0667 0.0000

1.0000 1.0000 0.9000 0.7000 0.6000 0.5000 0.2333 0.1000 0.0857 0.0571 0.0286 0.0000

60 0.2 - test [40] 0.5 - test [40] 0.7 - test [40] 0.85 - test [40] 0.2 - ADAPTIC [14] 0.5 - ADAPTIC [14] 0.7 - ADAPTIC [14] 0.85 - ADAPTIC [14] 0.2 - SYSAF 0.5 - SYSAF 0.7 - SYSAF 0.85 - SYSAF

Displacement w (mm)

50 40 30 20

0 100

200

Table 2 Adopted physical parameters. Configuration

L (cm)

h (cm)

fy,0

F1 (kN)

F2 (kN)

EHR3 EGR1 ZSR1

124 122 120

117 117 118

38.2 38.2 35.5

112 65 74

28 2.5 2.85

Table 3 Collapse temperatures (1C).

10

0

Fig. 14. Configuration of tested frames [40].

300

400

500

600

700

800

Temperature (ºC)

Configuration

Experimental

ADAPTIC

dif (%)

SYSAF

dif (%)

EHR3 EGR1 ZSR1

475 515 547

452 489 514

 4,8  5,0  6,0

455 491 515

 4,2  4,7  5,9

Fig. 13. Midspan displacements versus temperature.

Table 1 following [40,14]. The stress–strain diagrams are built following CEN EN 1993-1-2:2005 [30]. The test displacement at the midspan of the beam by [40], regarding temperature and each load level, are depicted in Fig. 13 together with the results achieved by SYSAF and ADAPTIC [14]. As one may note the results achieved by the developed code compare well with the ones achieved by ADAPTIC and reproduce the experimental behavior of the beam. 8.3. Steel frame tests Three series of scaled steel frames, called EHR, EGR and ZSR (Fig. 14), tested by and Rubert and Schaumann [40] are considered. All cross sections are of the type IPE 80. Frame series EHR and EGR are uniformly heated while only the left cell of frame ZSR is heated. These cases were numerically analyzed by [14,16,23,24,40] among others. In order to make comparisons our code (SYSAF) is compared with experimental values given by Rubert and Schaumann [40] and the numerical results presented by Izzuddin et al. [14] using their package ADAPTIC. In Table 2 one finds all required physical parameters. It should be noted that fy,0 means the yielding limit at 20 1C.

At 20 1C the elastic modulus is E ¼ 210 GPa, the adopted Poisson is ν ¼ 0:0 and the thermal dilatation coefficient is constant and equal to 14  10  6 1C  1. The elastic modulus (kE;θ ), the proportionality limit (kp;θ ) and yielding limit (ky;θ ) reduction factors are given by Table 1 and the stress–strain diagrams are built following CEN EN 1993-1-2:2005 [30]. We employed SYSAF with four cubic frame elements and the cross section discretization presented in Fig. 12(b). The experimental and numerical collapse temperatures are presented in Table 3. Both SYSAF and ADAPTIC achieved collapse temperatures lower than the experimental ones. Moreover relative differences are less than 6%. The experimental displacements [40], regarding temperature and load level, are depicted in Figs. 15–17 together with the results achieved by SYSAF and ADAPTIC [14]. As it can be seen our formulation gives good results for these cases when compared to experimental and numerical results achieved by [14] using DAPTIC. 8.4. Two storey three-dimensional frame The three-dimensional frame, Fig. 18, was firstly analyzed by Souza Júnior and Creus [23] using the package SAFIR. All sections are of the type H 150  150  7  10 and the material is steel with

78

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

Displacement (mm)

40

u2 - test [40] w4 - test [40] u2 - ADAPTIC [14] w4 - ADAPTIC [14] u2 - SYSAF w4 - SYSAF

30

20

10

0 0

100

200

300

400

500

Temperature (ºC) Fig. 15. Experimental and numerical results for EHR3.

u1 - test [40] u3 - test [40] u5 - test [40] u1 - ADAPTIC [14] u3 - ADAPTIC [14] u5 - ADAPTIC [14] u1 - SYSAF u3 - SYSAF u5 - SYSAF

30

20

Fig. 18. Two storey three-dimensional frame and boundary conditions. 80

10

70

0 0

100

200

300

400

500

600

Temperature (ºC ) Fig. 16. Experimental and numerical results for EGR1.

Displacement (mm)

Displacement (mm)

40

H1 - SAFIR [23]

60

H1 - SYSAF

50

H1 - ANSYS H2 - SAFIR [23]

40

H2 - SYSAF

30

H2 - ANSYS

20 10

60

0 u1 - test [40] u2 - test [40] u1 - ADAPTIC [14] u2 - ADAPTIC [14] u1 - SYSAF u2 - SYSAF

Displacement (mm)

50 40 30

-10 0

100

200

300

400

500

600

Temperature (ºC) Fig. 19. Displacement at points H1 and H2 as a function of temperature.

20

8.5. Steel beam-columns under fire

10 0 0

100

200

300

400

500

600

Temperature (ºC) Fig. 17. Experimental and numerical results for ZSR1.

f y ¼ 325 MPa. The load F ¼ 250 kN is applied before heating. The constitutive relation and relative thermal elongation follow CEN EN 1993-1-2:2005 [30] recommendations. The fire is simulated by a uniform heating of all first floor frame elements. Each member is discretized by four cubic frame elements. The cross section discretization follows the same pattern described in Fig. 12(b). Displacements at points H1 and H2 (Fig. 18) achieved by our formulation (SYSAF), ANSYS (BEAM189) and SAFIR [23] are compared in Fig. 19. As it can be seen all codes give similar results (less than 2.5% of difference for collapse temperature), confirming that our formulation is adequate to solve steel structures under high temperatures and uniform heating conditions. Moreover, our formulation applies three-dimensional constitutive relation, which indicates a more general approach for frame elements. Next examples show cases for which the cross section thermal analyses are required.

This example was presented by Landesmann et al. [22] and is the simulation of a pin-ended IPE-360 beam-column exposed to a standard time–temperature curve ISO 834 [41]. The two analyzed models are depicted in Fig. 20, including the geometry and both mechanical and thermal loads. The load combination (N and M) is fraction of the strength capacity of the column at 20 1C, i.e., Ny20 (compression) and Mp20 (bending). In the first situation the beam-column is heated at all faces (Fig. 20(a)) while in the second situation only three faces of the column are heated (Fig. 20(b)). Following [22] the last situation represents a facade column where one of the flanges is partially protected from fire. The adopted cross section discretization and the achieved temperature profile at 200 s, 500 s and 800 s of analysis using SYSAF are described, respectively, in Figs. 21–23. At 20 1C the material presents E ¼ 205 GPa and f y ¼ 250 MPa. The stress–strain diagram and relative thermal elongation elongation follow CEN EN 1993-1-2:2005 [30] recomendations. The achieved results by Landesmann et al. [22] using SAFIR and SAAFE are compared with our results at Fig. 24. It is important to note that SAAFE adopts concentrated plasticity, while SAFIR and SYSAF adopt distributed plasticity. One can see in Fig. 24 that our results are in good conformity with the ones presented by [22] using SAFIR, i.e., distributed plasticity.

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

8.6. Vogel's portal frame under fire The Vogel's portal frame [42] is solved here for a fire situation proposed by Landesmann et al. [22]. The original loading combination was applied at different load ratios ψ. The structure is exposed to heating according the ISO 834 curve [41]. All the members were assumed to be fire protected at the external faces (Fig. 25). We adopted four cubic finite elements for each structural member. The thermal analysis is carried out using the discretizations depicted in Fig. 26(a) and (b) for cross sections HEB 300 and HEA 340. Fig. 26 also shows the temperature fields at 750 s.

79

The relative thermal elongation and stress–strain relation follow CEN EN 1993-1-2:2005 [30]. Fig. 27 compares our results (SYSAF) with the ones presented by Landesmann et al. [22] (SAFIR) for load ratios ψ equal to 0.2, 0.4, 0.6 and 0.8. As one can see the results are in good agreement between formulations which implies that both formulations may be used to model structures under high temperatures. Our results present a more stiff behavior when compared to the one presented by Landesmann et al. [22] for low load factors. However, it is not possible to conclude that this difference is associated with difference in formulations, as the problem is very dependent on the large amount of variables and boundary conditions prescribed. However, we believe that using a 3D formulation, as ours, more realistic results can be achieved in more general analysis.

8.7. Steel beam submitted to non-uniform heating

Fig. 20. Pin-ended beam–column models: (a) Four faces heated section; (b) Three faces heated section. From: Landesmann et al. [22].

In order to explore the possibilities of the developed formulation we solve a typical case that demand thermal analysis for determining the temperature field of the steel profile during heating, i.e., temperature field cannot be determined by means of simplified equations. Moreover, a convergence analysis is performed regarding mechanical and thermal behavior of the solution indicating the fast convergence of the integration procedure. A steel beam with cross section H 200  100  5.5  8 embedded in a concrete slab (Fig. 28) submitted to fire is analyzed. Only the bottom surfaces of the beam and slab are exposed to effects of ambient heating according the ISO 834 standard fire curve [41]. Heat exchange to an ambient at 20 1C is also considered at the top of slab. The resultant emissivity for steel and concrete is assumed 0,7. The heat transfer convection coefficient is assumed 25 W=m2 o C at the surfaces exposed to fire and 4 W=m2 o C at the top of slab. The concrete is considered non-structural, acting as a heat sink material. The thermal properties of steel and concrete are in accordance with CEN EN 1994-1-2:2005 [43]. At 20 1C a simplified stress–strain diagram is adopted with elastic modulus E ¼ 210 GPa, reference yielding stress f y ¼ 355 MPa and

Fig. 21. Temperature field at 200 s (SYSAF). (a) Four sides heated and (b) three sides heated.

80

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

Fig. 22. Temperature field at 500 s (SYSAF). (a) Four sides heated and (b) three sides heated.

Fig. 23. Temperature field at 800 s (SYSAF). (a) Four sides heated and (b) three sides heated.

Poisson ratio ν ¼ 0:0. Perfect elastoplasticity is adopted until a strain of 0.15, when a softening is imposed reducing stress to zero at a strain of 0.2.

8.7.1. Mechanical convergence analysis Only two cubic elements are able to capture the exact elastic bending moment distribution of the proposed static scheme. However, for the elastoplastic analysis the integration points, responsible to capture residual stresses, should be near to the

center of the beam. Therefore, we used four finite elements, two with 195 cm length near the ends and two with 5 cm near the center to model the structure. No sensible changes in results are achieved when varying the central elements lengths. The very good convergence of the plastic behavior can be seen in Fig. 29 for which two cross section discretizations are tested, one very poor with 18 triangular divisions, see Fig. 30a, and other with 62 triangular divisions, see Fig. 30b. As one can see, the very poor discretization presents practically the same results as the regular one, characterizing convergence. Moreover, the maximum

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

difference in force among each other is less than 0.27%. Thirty seven displacement steps of 1 mm each are imposed at the center of the beam.

Midspan horizontal displacement (cm)

SAFIR - 4 sides [22] SAAFE - 4 sides [22]

4

SYSAF - 4 sides SAFIR - 3 sides [22]

3

The achieved collapse load is 69:8 kN, 1.7% less than a reference value (71 kN) calculated using the uniaxial yielding stress, revealing the influence of the shear stress in the collapse. 8.7.2. Thermal convergence analysis The same pattern of steel cross section discretizations, used for the mechanical convergence analysis, are used to solve the thermal problem, now including the presence of concrete acting as a heat sink material, as shown in Fig. 31. Fig. 32 shows the ISO 834 standard fire curve, which determines the ambient temperature evolution, and temperature of points A and B (see fig. 28) along time for both discretizations. The maximum time for thermal analysis (240 min) was adopted considering the double of the maximum required fire resistance time found in design codes. As one can see no sensible difference

6 5

81

SAAFE - 3 sides [22] SYSAF - 3 sides

2 1 0 0

200

400

600

800

1000

10

Fire elapsed time (s)

Fig. 24. Midspan horizontal displacement versus time – results comparison. Displacement Δ (cm)

= 0.2, 0.4, 0.6, 0.8 E = 205 GPa y = 235 MPa F = · 2800 kN H = ·35 kN

8

F

ψ = 0.8

ψ = 0.6

ψ = 0.4

ψ = 0.2

6

4 SAFIR [22]

2

HEA 340

SYSAF

0 0

250

500

750

1000

1250

1500

1750

Fig. 27. Vogel's portal frame under fire – Lateral displacement at the top (Δ).

220 mm

5m

H E B 30 0

H EB 300

Fire elapsed time (s)

B

A

Force

4.0 m

4m Fig. 25. Vogel's portal frame under fire as in Landesmann et al. [22].

Fig. 28. Beam submitted to non-uniform heating.

Fig. 26. Discretizations and temperature fields at 750 s. (a) HEB 300 and (b) HEA 340.

82

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

among the presented results for both discretizations is achieved, revealing the very good convergence of the thermal problem when using the adopted cubic finite element.

75

Based on the achieved results, it is evident that the difference between the top and bottom temperatures has a large dependence regarding time. Which means that, for this case, the use of a simplified temperature profile during the heating is not an option, especially when considering the initial time, until 120 min, the maximum found in design codes.

70 65 60 55

Force (kN)

50 45 40 35

Reference value 62 elements 18 elements

30 25 20 15 10 5 0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

Displacement (cm) Fig. 29. Convergence of the collapse load.

8.7.3. Thermo-elastoplastic analysis As simplified temperatures profile does not apply for this case, we present a complete thermo-mechanical analysis to serve as a reference for further works and also to perform a final convergence analysis. In order to do so, the CEN EN 1993-1-2:2005 [30] temperature dependence of the constitutive model and relative thermal elongation for steel at high temperatures are considered. From these considerations, the temperature profile for the steel cross section achieved by the previous thermal analysis is used to model the simple supported beam indicated in Fig. 28 subjected to a transverse load of F ¼ 60 kN. Fig. 33 shows the vertical displacement versus time and indicates that the collapse occurs in 25 min. Moreover, the convergence is also characterized as the use of the very poor discretization presents practically the same result as the one achieved by the final discretization of the thermal and mechanical convergence analyses.

1200 1100 1000

Temperature (oC)

900 800 700

ISO 834 Point A - 84 elements Point B - 84 elements Point A - 350 elements Point B - 350 elements

600 500 400 300 200 100 0 0

20

40

60

80

100 120 140 160 180 200 220 240

Time (min) Fig. 30. Adopted cross section discretizations. (a) 18 elements. (b) 62 elements.

Fig. 32. Temperature of points A and B along time – thermal convergence analysis.

Fig. 31. Cross section and concrete discretizations – thermal convergence analysis. (a) 84 elements. and (b) 350 elements.

R. Rigobello et al. / Finite Elements in Analysis and Design 91 (2014) 68–83

Displacement at the midspan (cm)

0 -2 -4 -6 -8 -10 -12

Cross section - 18 elements Cross section - 62 elements

-14 -16 -18 0

5

10

15

20

25

30

Time (min) Fig. 33. Displacement at the midspan versus time in the heating process.

9. Conclusions In this study we successfully apply a new geometrical nonlinear formulation, based on positions, to model three-dimensional frame structures under high temperatures (fire). The formulation is geometrically exact and does not use finite rotation, but generalized vectors, to represent large displacements and rotations for frame analysis. A two-dimensional nonlinear thermal module of analysis is developed and implemented to solve the cross section temperature profile. The three-dimensional thermo-mechanical coupling is developed adjusting, at the equivalent stress space, the stress–strain diagrams following standardized reduction factors. Moreover the adjustment of internal force, respecting the yielding surface at any temperature, is proposed and successfully tested. Numerical implementations have been carried out resulting in the SYSAF package. Various examples are solved, comparing results with experimental and other numerical methods, confirming that the proposed formulation is able to solve steel structures in fire. Acknowledgments Authors would like to thanks FAPESP (São Paulo Research Foundation – Brazil) and CNPq (National Council for Scientific and Technological Development – Brazil) for the financial support of this research. References [1] R.R. Sun, Z. Huang, I.W. Burgess, Progressive collapse analysis of steel structures under fire conditions, Eng. Struct. 34 (2012) 400–413. [2] J.-M. Franssen, V. kodur, R. Zaharia, Designing steel structures for fire safety, CRC Press, Leiden, 2009. [3] J.A. Purkis, Fire Safety Engineering – Design of Structures, second ed., Butterworth-Heinemann, Oxford, 2007. [4] P. Vila Real, Incêndio em Estruturas Metálicas – Cálculo Estrutural, Edições Orion, Mafra, 2003. [5] Y.C. Wang, Steel and Composite Structures, - Behaviour and Design for Fire Safety, Spon Press, London, 2002. [6] Y.C. Wang, Advances in research on fire engineering of steel structures, P. I. Civil Eng –Civ. En 162 (3) (2009) 129–135. [7] Steel Construction Industry Forum, Investigation of Broadgate Phase 8 Fire, The Steel Construction Institute, Ascot, 1991. [8] T. Lennon, D. Moore, Client report: results and observations from full-scale fire test at BRE CardingtonBuilding Research Establishment, Client Report Number 215-741 2004. [9] J.-M. Franssen, SAFIR – a thermal/structural program modelling structures under fire, Eng. J. AISC 42 (3) (2005) 143–158. [10] J.-M. Franssen, User's Manual for SFIR 2007 – A Computer Program for Analysis of Structures Subjected to Fire, University of Liège, Liège, 2007.

83

[11] Z. Huang, I.W. Burgess, R.J. Plank, 3D modeling of beam-columns with general cross-sections in fire, in: Proceedings of the 3rd International Workshop – Structures in Fire, SiF'04 (2004), pp. 323–334. [12] C. Yu, Z. Huang, I.W. Burgess, R.J. Plank, Development and validation of 3D composite structural elements at elevated temperatures, J. Struct. Eng. – ASCE 136 (3) (2010) 275–284. [13] B.A. Izzuddin, A.S. Elnashai, ADAPTIC: A Program for the Adaptive Dynamic Analysis of Space Frames, Tech Report No. ESEE-89/7, Imperial College, London, 1989. [14] B.A. Izzuddin, L. Song, A.S. Elnashai, P.J. Dowling, An integrated adaptive environment for fire and explosion analysis of steel frames – Part II: verification and application, J. Constr. Steel Res. 53 (1) (2000) 87–111. [15] L. Song, B.A. Izzuddin, A.S. Elnashai, A.S., P.J. Dowling, An integrated adaptive environment for fire and explosion analysis of steel frames – Part I: analytical models, J. Constr. Steel Res. 53 (1) (2000) 63–85. [16] K.H. Tan, S.K. Ting, Z.F. Huang, Visco–elasto-plastic analysis of steel frames in fire, J. Struct. Eng. – ASCE 128 (1) (2002) 105–114. [17] Z.F. Huang, K.H. Tan, Structural responses of axially restrained steel beams with semirigid moment connection in fire, J. Struct. Eng. – ASCE 131 (4) (2005) 541–551. [18] Z.F. Huang, K.H. Tan, S.K. Ting, Heating rate and boundary restraint effects on fire resistance of steel columns with creep, Eng. Struct. 28 (2006) 805–817. [19] W.S. Toh, T.C. Fung, K.H. Tan, Fire resistance of steel frames using classical and numerical methods, J. Struct. Eng. – ASCE 127 (7) (2001) 829–838. [20] C.K. Iu, S.L. Chan, A simulation-based large deflection and inelastic analysis of steel frames under fire, J. Constr. Steel Res. 60 (2004) 1495–1524. [21] C.K. Iu, S.L. Chan, X.X. Zha, Nonlinear pre-fire and post-fire analysis of steel frames, Eng. Struct 27 (2005) 1689–1702. [22] A. Landesmann, E.M. Batista, J.L.D. Alves, Implementation of advanced analysis method for steel framed structures under fire conditions, Fire Saf. J. 40 (4) (2005) 339–366. [23] V. Souza Junior, G.J. Creus, Simplified elastoplastic analysis of general frames on Fire, Eng. Struct. 29 (2007) 511–518. [24] K.H. Lien, Y.J. Chiou, R.Z. Wang, P.A. Hsiao, Vector form intrinsic finite element analysis of nonlinear behavior of steel structures exposed to fire, Eng. Struct. 32 (2010) 80–92. [25] H.B. Coda, A solid-like FEM for geometrically non-linear 3D frames, Comput. Method Appl. M. 198 (47-48) (2009) 3712–3722. [26] H.B. Coda, R.R. Paccola, Improved finite element for 3D laminate frame analysis including warping for any cross section, Appl. Math. Model. 34 (2010) 1107–1137. [27] H.B. Coda, R.R. Paccola, A FEM procedure based on positions and unconstrained vectors applied to non-linear dynamic of 3D frames, Finite Elem. Anal. Des. 42 (4) (2011) 319–333. [28] H.C. Huang, A.S. Usmani, Finite Element Analysis for Heat Transfer, SpringerVerlag, London, 1994. [29] R.W. Lewis, P. Nithiarasu, K.N. Seetharamu, Fundamentals of the Finite Element Method for Heat and Fluid Flow, John Wiley & Sons, Chichester, 2004. [30] European Committee for Standardization – CEN, EN 1993-1-2:2005 Eurocode 3 – Design of Steel Structures – Part 1–2: General rules – Structural Fire Design, Brussels, 2005. [31] J. Bonet, R.D. Wood, J. Mahaney, P. Heywood, Finite element analysis of air supported membrane structures, Comput. Method Appl. M. 190 (5–7) (2000) 579–595. [32] H.B. Coda, R.R. Paccola, An alternative positional FEM formulation for geometrically nonlinear analysis of shells: Curved triangular isoparametric elements, Comput. Mech. 40 (1) (2007) 185–200. [33] R. Carrazedo, H.B. Coda, Alternative positional FEM applied to thermomechanical impact of truss structures, Finite Elem. Anal. Des. 46 (11) (2010) 1008–1016. [34] H.B. Coda, Two dimensional analysis of inflatable structures by the positional FEM, Latin Am. J. Solids Struct. 6 (3) (2009) 187–212. [35] C. Lanczos, Variational Principles of Mechanics, University of Toronto Press, Toronto, 1970. [36] A.S. Botta, R.R. Paccola, W.S. Venturini, H.B. Coda, A discussion on volume change in the plastic phase, Commun. Numer. Meth. En. 24 (2008) 1149–1162. [37] J.C. Simo, T.J.R. Hughes, Computational Inelasticity, Springer-Verlag, New York, 1998. [38] R. Cook, D. Malkus, M. Plesha, Concepts and applications of finite element analysis, third ed., John Willey & Sons, New York, 1989. [39] R. Rigobello, H.B. Coda, J. Munaiar Neto, Inelastic analysis of steel frames with a solid-like finite element, J. Construct. Steel Res. 86 (2013) 140–152. [40] A. Rubert, P. Schaumann, Structural steel and plane frame assemblies under fire action, Fire Saf. J. 10 (1986) 173–184. [41] International Organization for Standardization – ISO 834-1 – Fire-resistance tests – elements of building construction – Part 1: general requirements, Geneva, 1999. [42] U. Vogel, Calibrating Frames, Stahlbau 10 (1985) 295–301. [43] European Committee for Standardization – CEN, EN 1994-1-2: 2005 Eurocode 4 – Design of composite steel and concrete structures – Part 1–2: General Rules – Structural Fire Design, Brussels, 2005.