A mixed finite element for three-dimensional nonlinear analysis of steel frames

Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545 www.elsevier.com/locate/cma

A mixed finite element for three-dimensional nonlinear analysis of steel frames Phani Kumar V.V. Nukala a, Donald W. White a

b,*

Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831-6359, USA School of Civil and Environmental Engineering, Georgia Institute of Technology, Atlanta, GA 30332-0355, USA

b

Received 16 March 2003; received in revised form 24 November 2003; accepted 21 January 2004

Abstract This paper develops and demonstrates a mixed beam finite element that can represent the inelastic three-dimensional stability behavior of frames composed of open-walled cross-section members. The kinematics of deformation of the element include finite rotation and warping of the cross-section due to torsion. The material inelasticity is based on a two-space model that includes the effect of shear stresses due to uniform torsion in addition to normal stresses due to axial force, biaxial bending, and bimoment. Initial geometric imperfections and residual stresses are accommodated. The formulation is based on a two-field (displacement and generalized stress) Hellinger–Reissner (HR) variational principle. The interpolation of the generalized stresses (i.e., cross-section stress resultants) along the element length is based on the geometrically-exact nonlinear governing differential equations of equilibrium. This leads to significant improvements in the accuracy over conventional displacement-based elements in problems involving highly nonlinear variations in the curvature along the member lengths due to geometric nonlinearity and/or distributed plasticity effects. Several numerical examples involving two- and three-dimensional elastic and inelastic stability behavior are presented to illustrate the capabilities of the formulation.
2004 Elsevier B.V. All rights reserved. Keywords: Mixed finite element; Three-dimensional nonlinear analysis; Steel frames

1. Introduction In the last several decades, considerable progress has been made on the three-dimensional nonlinear analysis of structural members and framing systems subjected to large displacements and rotations. Notable contributions in this area include, but are not limited to [1–15]. Many of these contributions have focused only on geometrically nonlinear analysis. For purposes of computational efficiency, many formulations that address material nonlinearity have been based on concentrated plasticity idealizations, that

*

Corresponding author.

0045-7825/$ - see front matter
2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.01.029

2508

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

is, the plasticity is modeled within zero-length plastic hinges [16,17]. Although these types of idealizations are often sufficient for accurate analysis of the overall inelastic stability behavior of framing systems, generally they are inadequate for prediction of limit states such as inelastic lateral buckling of beams and torsional-flexural buckling of beam-columns. In addition to the effects of initial geometric imperfections and geometric nonlinearity, distributed plasticity effects––including initial residual stresses––must be modeled in some way for an accurate representation of the general inelastic stability behavior of steel framing members. Some progress has been made in recent research on direct modeling of the material nonlinearity directly in terms of cross-section stress resultants [18]. Such models require numerical integration only along the element length, and thus eliminate the high computational cost of numerical integration over the element cross-section. However, stress resultant plasticity models that can predict reliably the detailed coupling of bending, axial, and torsional deformations due to partial cross-section plasticity, the influence of shear stresses due to torsion on the axial and bending behavior, and the effects of nonproportional cycling do not exist at the present time. Therefore, stress-space constitutive models and numerical integration both over the cross-section and along the element lengths are still required for rigorous assessment of three-dimensional inelastic member stability behavior. This type of analysis is often referred to as a plastic-zone analysis or distributed plasticity approach [10,19–22]. Since the section curvatures and axial strains can vary in a highly nonlinear fashion along the length of an inelastic structural member, a large number of displacement-based beam finite elements are typically required to represent the distributed plasticity behavior accurately. For example, Izzuddin and Smith [23,24] utilize from two to ten elements along the member lengths for analysis using their element, depending on the problem being addressed. Teh and Clarke [25] have found that three to five of their elements are required per member to obtain accurate solutions. This is because the approximations of the axial strains and curvatures of a displacement-based element are constrained by the element’s assumed displacement field. Mixed elements have tremendous potential to alleviate this problem, particularly if the independent assumed stress resultant fields are specified to satisfy the governing differential equations of equilibrium of the deformed member convergence requirements and For example, in a first-order distributed plasticity analysis of framing members subjected only to end-forces, the exact variation of the moments along the element length is linear. If a mixed element formulation is used that is based on a linear interpolation of the moments, then a model can be established in which the differential equations of equilibrium are satisfied section-by-section along the length of the element. An improved representation of the curvatures along inelastic element lengths is achieved by the use of this force interpolation, and the resulting element exhibits improved coarse mesh accuracy relative to comparable displacement-based formulations [20]. This paper presents a three-dimensional, geometric-nonlinear, distributed plasticity beam element based on the two-field (displacement and generalized stress) Hellinger–Reissner (HR) variational principle. This element is targeted particularly at analysis of the three-dimensional inelastic stability behavior of members of open-walled cross-section. The kinematics of deformation include finite rotation, and warping of the cross-section due to torsion. The inelastic material behavior is modeled using a von-Mises representation that captures the coupling between the shear stresses due to uniform torsion and the normal stresses due to axial force, biaxial bending, and bimoment. The proper consistent linearization of the governing variational principle to obtain the element tangent stiffnesses, and state determination algorithms is outlined. The organization of the paper is as follows. In Section 2, the three-dimensional beam formulation is summarized. Section 3 describes the finite element discretization and the corresponding state determination procedure. Numerical examples that illustrate the capabilities of the element are presented in Section 4. A detailed investigation of four alternative consistent state determination algorithms for this and other similar mixed beam element formulations is presented in a companion paper [26].

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2509

2. Formulation 2.1. Kinematic description In this work, we adopt a specialized form of the kinematic hypothesis presented by Simo and Vu-Quoc [8] (see Fig. 1). The reference configuration ðB
R3 Þ is chosen as a straight beam of length equal to the length along the initial (straight or curved) geometry of the beam. Curved geometry and/or initial geometric imperfections are modeled as an initial displacement effect, the details of which are addressed within the subsequent developments. Let S denote the coordinate along the axis of the beam within the reference configuration. We define an orthonormal material reference frame, at a typical beam cross-section located at S, by fO; E1 ðSÞE2 ðSÞE3 ðSÞg such that E3 is directed along the axis of the beam within the reference configuration. Furthermore, let the orthonormal spatial fixed frame be fe1 ; e2 ; e3 g such that ei ¼ diI EI . Let fsI ðSÞ; I ¼ 1; 2; 3g represent the orthonormal basis vectors of a moving frame attached to the crosssection located at S in the reference configuration. This orthonormal moving reference frame coincides with the fixed material reference frame in the reference configuration, i.e., sI ðSÞ ¼ EI ðSÞ. The finite rotation of the cross-section is defined by specifying the orientation of the moving reference basis sI with respect to the fixed (material) reference basis EI as sI ¼ KðSÞEI ¼ KiI ðSÞei ;

ð1Þ

where K is the orthogonal rotation tensor. The deformed configuration of a beam in three-dimensions is specified by the displacement of the reference axis, the finite rotation of the cross-sections, and finally a superposed out-of-plane warping relative to the cross-section plane. Thus, the position vector x, which describes the location of a material point in the deformed configuration that was located at X ¼ ðX1 ; X2 ; SÞ in the reference configuration, is given by x ¼ /ðX1 ; X2 ; SÞ ¼ /0 ðSÞ þ Xa sa ðSÞ þ -ðX1 ; X2 ÞpðSÞs3 ðSÞ where a ¼ 1; 2;

s1

g

s3

ΦB Φ E1

e1

O

e3

B

E3

S Fig. 1. Kinematic description of the beam.

3

ð2Þ

2510

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

where /0 ðSÞ ¼ uE1 þ vE2 þ ðS þ wÞE3 :

ð3Þ

The vector /0 ðSÞ describes the position vector of a point on the beam axis in the deformed configuration, which was initially located at S in the reference configuration. The beam axis can be taken at any location within the cross-section, although the description of the elastic response is simplified by locating the beam axis at the centroid of the cross-section. Also, -ðX1 ; X2 Þ denotes the warping function of the cross-section and pðSÞ represents the warping amplitude, as in [8]. For nonprismatic sections, a general warping function -ðX1 ; X2 ; SÞ, that changes along the length of the element may be used. In the above expression for /0 ðSÞ, u; v, and w are the displacements of the reference point within a cross-section located at S, in the directions of the material bases E1 , E2 and E3 , respectively. Thus the deformed configuration of a three-dimensional beam UðSÞ ¼ Uð/0 ðSÞ; KðSÞ; pðSÞÞ is uniquely determined by specifying the position vector /0 ðSÞ, the finite rotation transformation KðSÞ and the warping amplitude pðSÞ. Although the general kinematic model of Simo and Vu-Quoc [8] includes beam-shear effects, beam-shear deformations are not significant in many types of frame members. It is traditionally accepted that beamshear effects tend to be negligible in members that have lengths greater than about 10 times their largest cross-section dimension [27–29]. Beam-shear deformations are assumed to be negligible in the element developed here. This leads to significant simplification of the element kinematic equations, and avoids potential problems associated with shear locking in framing members with lengths significantly larger than their cross-section dimensions. When beam-shear deformations are included within the formulation, the cross-section displacements and rotations are interpolated independently with C 0 continuous interpolations. On the other hand, when shear deformations are not included, Hermitian interpolations (C 1 continuous) are used for the displacements. A higher order interpolation for the displacements typically increases the coarse-mesh accuracy of the solutions. Furthermore, warping shear strains due to nonuniform torsion (i.e., bi-shear strains [8]) are neglected in this work. These shear strains tend to be an order of magnitude less than those due to uniform torsion in framing members of ordinary proportions. Also, neglect of warping shear is consistent with the assumption that the beam-shear effects are small. The effects of nonuniform torsion are manifested predominantly as axial strains in most open-section frame members. In this work, it is assumed that the interaction of axial and warping shear stresses is not significant, just as it is assumed that the interaction of axial and beam-shear stresses is negligible. However, rigorous assessment of inelastic uniform torsion and its coupling with bending, axial, and warping deformations requires the direct consideration of the associated St. Venant shear strains. In this study, the shear strains due to uniform torsion are included based on the mitre model [30,31]. In [10,30,31], this model is calibrated specifically for I-sections and is shown to provide good accuracy. 2.2. Finite rotation parametrization In the literature, various parameters have been chosen to represent the orthogonal transformation operator KðSÞ. These include Euler angles, quaternion or Euler–Rodrigues parameters, the conformal rotation vector, the rotational vector and the rotational pseudo-vector [6,32]. By using the rotational vector parametrization, the orthogonal transformation operator K can be represented as
 32 2
h sin
h 1 4 sin 2 5 2
 K¼Iþ
Hþ H
h 2 h

ð4Þ

2

¼I þHþ

1 2 H þ    ¼ expðHÞ; 2!

ð5Þ

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2511

where the rotational axial vector
h is defined as
h ¼
h^n, and
h is the magnitude of the rotation about an axis ^ specified by the unit vector n. Also, H represents the skew-symmetric tensor corresponding to the axial vector
h, such that Hh ¼
h  h 8h 2 R3 , i.e., in matrix form 8 9 2 3
h2 0 
h3
> = < h1 > 6
7 ð6Þ ½H ¼ 4 h3 hg ¼
h2 : 0 
h1 5 f
> ; :
>
1
2 h h3 h 0 The variation dK of the orthogonal finite rotation tensor can be expressed as dK ¼

d d K
j ¼ ½expð
dWÞKj
¼0 ¼ dWK; d

¼0 d


ð7Þ

where dW is a skew-symmetric tensor (with a corresponding axial vector dw) and is an element of the tangent space at K, TK ðSOð3ÞÞ, i.e., dW is an infinitesimal rotation tensor superposed on the existing rotation K. Hence, dW is a spatial object defined as dW ¼ dKKT :

ð8Þ

Similarly, the variation of the finite rotation tensor can be constructed in terms of the material object dH as dK ¼

d d K
j ¼ ½K expð
dHÞj
¼0 ¼ KdH; d

¼0 d


ð9Þ

where dH is a skew-symmetric tensor (with a corresponding axial vector dh) belonging to the tangent space at identity, TI ðSOð3ÞÞ, and is expressed as dH ¼ KT dK:

ð10Þ

The spatial and material objects dW and dH are related as follows:  dW ¼ KdHKT ; dw ¼ Kdh : dH ¼ KT dWK; dh ¼ KT dw

ð11Þ

As explained above, the finite rotation tensor K can be parameterized in terms of various quantities. The following description is particularly useful as an intermediate step in the derivations for elements that are based on the assumption of zero beam-shear deformations. This description involves the definition of two ^ and one angle ^c that establishes the orientation of an axis N b about which the angles of rotation #^ and /, ^ rotation # occurs (see Fig. 2). By using these parameters, a finite rotation can be defined as:

Fig. 2. Finite rotation transformation.

2512

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

b which is perpendicular to E3 ðSÞ and • an initial rotation through an angle #^ is performed about the axis N, b s3 ðSÞ, where the angle between N and E2 is denoted by the symbol ^c, followed by ^ about the current ‘‘longitudinal’’ unit vector s3 ðSÞ. • a rotation / The same finite rotation transformation matrix is obtained if the rotations are performed in the ^ about the fixed material reference vector E3 ðSÞ, followed by a reverse order as an initial rotation / ^ b , which is perpendicular to the plane containing E3 ðSÞ and rotation through an angle # about the axis N s3 ðSÞ. The finite rotation tensor can be expressed in matrix form in terms of the above three parameters as 2 3 lx ly lz ½KiI  ¼ 4 mx my mz 5; ð12Þ nx ny nz where

9 ^ sin ^c þ cosð^c  /Þ ^ cos #^ cos ^c > lx ¼ sinð^c  /Þ > > > ^ sin ^c þ sinð^c  /Þ ^ cos #^ cos ^c > > ly ¼ cosð^c  /Þ > > > > ^ lz ¼ cos ^c sin # > > > ^ ^ ^ > mx ¼ sinð^c  /Þ cos ^c þ cosð^c  /Þ cos # sin ^c > = ^ ^ ^ my ¼ cosð^c  /Þ cos ^c þ sinð^c  /Þ cos # sin ^c > > > > mz ¼ sin ^c sin #^ > > > ^ > nx ¼ sin #^ cosð^c  /Þ > > > > ^ ^ > ny ¼ sin # sinð^c  /Þ > > ; ^ nz ¼ cos #

ð13Þ

are the direction cosines of sI relative to EI . If the variation of the finite rotation tensor dK, expressed in terms of the pseudo-vector like parameters ^ T , is related to the variation of the material object dh such that both parametrizations lead ^ d^c; d/g dv ¼ fd#; to the same perturbed finite rotation tensor K
, i.e., K
ðv þ
dvÞ ¼ expðHÞ expð
dHÞ

ð14Þ

the following result is obtained: dh ¼ Tdv ! dh ðvÞdv;

ð15Þ

where 2

3 ^ cosð^c  /Þ ^ sin #^ 0 sinð^c  /Þ ^ ^ sin #^ 0 5: ð16Þ Tdv ! dh ðvÞ ¼ 4 cosð^c  /Þ sinð^c  /Þ ^ 0 ð1 cos #Þ 1 ^ T to represent the finite rotation transformation, the issues of ^ ^c; /g By using the parametrization v ¼ f#; consistent linearization and updating of the rotation variables are greatly simplified, since v and Dv belong to the same vector space. Similar finite rotation parametrizations have been adopted in [6,11]. In a similar manner to the above, the relationship between the variation dw and dv is obtained as dw ¼ TH dv ! dw ðvÞdv;

ð17Þ

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

where

3 sin ^c cos ^c sin #^ cos ^c sin #^ 4 cos ^c sin ^c sin #^ sin ^c sin #^ 5: TH dv ! dw ðvÞ ¼ ^ 0 ð1 cos #Þ cos #^

2513

2

ð18Þ

It should be noted that the transformations Tdv ! dh and TH dv ! dw are singular for values of # ¼ kp where k ¼ 0; 1; 2; 3; . . ., since detðTdv ! dh Þ ¼ detðTH dv ! dw Þ ¼ sin #:

ð19Þ

H1 With the exception of these singularities, the inverse mappings T1 dv ! dh and Tdv ! dw are well defined and can be obtained by simply inverting the respective matrices. The commutative diagram shown in Fig. 3 summarizes the relations given by Eqs. (11), (15) and (17) for representing the variation of the finite rotation tensor.

Remark 1. From the commutative diagram shown in Fig. 3, we obtain T H 1 dw ¼ TH dv ! dw dv ¼ Tdv ! dw Tdv ! dh K dw ¼ Idw

ð20Þ

1 which leads to the relation TH dv ! dw Tdv ! dh ¼ K. This relationship can be verified by simply multiplying the matrices.

The element proposed in this paper is similar to previously proposed C 1 -type elements for eliminating the so-called membrane/shear/warping locking in the straight and curved beam formulations [10,33,34]. More recently, similar approach has been used in [35–38] to make the beam formulation fully consistent with the existing shell formulations. In these C 1 -type formulations, the warping as well as the bending rotational degrees of freedom are associated with the spatial derivatives of the three reference-line translations ^ In the following, we derive the relationship between the fu; v; wg, and the cross-section twist rotation /. rotations and the displacement derivatives. 2.2.1. Relationship between rotations and displacement derivatives Fig. 4 illustrates the relationship between a gage length AB along the beam axis in the reference configuration and the corresponding length within a deformed configuration. Based on the assumption that zero beam-shear deformation, the axial strain of the beam axis,
, is related to the displacement derivatives by the equation 2

2

ð1 þ
Þ ¼ u02 þ v02 þ ð1 þ w0 Þ ;

Fig. 3. Finite rotation parametrizations––commutative diagram.

ð21Þ

2514

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

Fig. 4. Geometric representation of deformation.

where the ‘‘0’’ symbol denotes differentiation with respect to the coordinate S. Furthermore, based on the ^ T of the previous section can be expressed in terms ^ ^c; /g geometry shown in Fig. 4, the parameter set v ¼ f#; of the displacement derivatives as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðu02 þ v02 Þ 1 þ w0 ; ð22Þ and sin#^ ¼ cos #^ ¼ 1þ
1þ
u0 cos ^c ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðu02 þ v02 Þ

and

v0 sin^c ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : ðu02 þ v02 Þ

ð23Þ

By use of these relations, the variation of the material rotation axial vector, dh, along with the variation of the constraint condition (Eq. (21)), can be expressed in terms of the variation of the displacement deriv^ as follows: atives and the angle / fdhg ¼ ½Ca ½Ta fdag ¼ ½Tda ! dh fdag;

ð24Þ

where fdhg ¼ f dh1

dh2

fdag ¼ f du0

dv0

d
g

ð25Þ

^ g; dw0 d/ 3 0 07 7 15 0

ð26Þ

dh3

and 2

0 1 61 0 ½Ca  ¼ 6 40 0 0 0 and

0 0 0 1

2

lx 1 6 l y 6 ½Ta  ¼ ð1 þ
Þ 4 ð1 þ
Þlz kmz

mx my ð1 þ
Þmz klz

ð27Þ

nx ny ð1 þ
Þnz 0

3 0 0 7 7; 0 5 ð1 þ
Þ

ð28Þ

1 . In the above set of equations, dh1 , dh2 and dh3 are the components of the material vector where k ¼ 1þcos #^

dh, and d
is the variation of axial strain at the reference point in the cross-section. It should be noted that the transformation from fdag to fdhg is objective (one-to-one and onto) in the sense that the matrix ½Tda ! dh  is nonsingular since detðTda ! dh Þ ¼

1 2

ð1 þ
Þ

:

ð29Þ

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2515

2.3. Generalized axial strain Following the standard procedures illustrated in [4–6,8,11], it is possible to derive the following strain measures for the beam element, c ¼ /00  s3 ;

ð30Þ

C ¼ KT ð/00  s3 Þ ¼ KT /00  E3 ;

ð31Þ

where c and C denote spatial and material strain measures, respectively, and /00 is the derivative of the position vector of the beam axis with respect to the coordinate S. Based on the constraint of zero shear deformations due to bending, we can write /00 ¼ u0 E1 þ v0 E2 þ ð1 þ w0 ÞE3 ¼ ð1 þ
Þðlz E1 þ mz E2 þ nz E3 Þ ¼ ð1 þ
Þs3 :

ð32Þ ð33Þ

As a result, the spatial and material strains can be written as c ¼ /00  s3 ¼
s3 ;

ð34Þ

C ¼ KT ð/00  s3 Þ ¼
E3 :

ð35Þ

The variation of the generalized strain measures c (Eq. (30)) and C (Eq. (31)) can be expressed as 9 od/0 > >  dw  s3 > > oS > > > = T od/0 T o/0 þ ðK Þ  dh : dC ¼ K oS oS > >   > > > od/0 o/ > > ;  dw  0 ¼ KT oS oS dc ¼

ð36Þ

However, based on the assumption of zero shear deformations due to bending (see Eqs. (34) and (35)), the variation of the axial strain at the location of the reference axis can be written as dC ¼ d
E3 ;

ð37Þ

where d
¼ lz du0 þ mz dv0 þ nz dw0 :

ð38Þ

2.4. Curvatures By differentiating the identity KT K ¼ KKT ¼ I with respect to S [4–6,8,11], the spatial and material curvature tensors (XðSÞ and KðSÞ, respectively) are given by XðSÞ ¼

oKðSÞ T K ðSÞ; oS

ð39Þ

oKðSÞ ; oS

ð40Þ

KðSÞ ¼ KT ðSÞ

2516

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

where oK ¼ XK ¼ KK: oS

ð41Þ

Based on the procedure presented earlier in Section 2.2, the spatial and material curvatures can also be expressed in terms of v0 as 0 x ¼ TH dv ! dw ðvÞv ;

ð42Þ

j ¼ Tdv ! dh ðvÞv0 ;

ð43Þ

where ‘‘0’’ denotes the derivative with respect to the coordinate S, and j and x denote the material and spatial curvature axial vectors, respectively. In terms of the displacement derivatives, the components of the material curvature axial vector and the derivative of axial strain with respect to S can be represented together as fjg ¼ ½Ca ½Ta fa0 g ¼ ½Tda ! dh fa0 g;

ð44Þ

where fjg ¼ f j1

j2

j3

v00

w00

T


0 g41

ð45Þ

and fa0 g ¼

n

u00

^0 /

oT 41

ð46Þ

:

The linearization of the spatial and material curvature axial vectors can be obtained as [4–6,8,11] 9 odw odh > >  x  dw ¼ K dx ¼ oS oS = : ð47Þ > odh > T odw ; dj ¼ K ¼ þ j  dh oS oS In terms of displacement derivatives (higher), the variation of the terms in fjg (see Eqs. (44) and (45)) are obtained after some manipulation as d
d
fjg ¼ ½Ca ½½Ta fda0 g þ ½Sa fdhg  fjg ð1 þ
Þ ð1 þ
Þ d
d
¼ ½Ca ½½Ta fda0 g þ ½Sa ½Ca ½Ta fdag  fjg ¼ ½Tca fda0 g þ ½Rca fdag  fjg; ð1 þ
Þ ð1 þ
Þ

fdjg ¼ ½Ca ½½Ta fda0 g þ ½dTa fa0 g 

ð48Þ

where 2 6 6 ½Sa  ¼ 6 4

0
0 j1 00 ^ þ v00 sin /Þ ^ kðu cos /


0

j1

0 j2 j2 0 00 00 ^ ^ kðv cos /  u sin /Þ 0

0

3

0 7 7 7 ð1 þ
Þ 5

ð49Þ

1

and ½Rca  ¼ ½Ca ½Sa ½Ca ½Ta :

ð50Þ

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2517

2.5. Deformation gradient The deformation gradient tensor, which relates an arbitrary infinitesimal material vector dX at X in the reference configuration with a vector dx at x in the deformed configuration, is defined as F ¼ xr ¼

o/ðX1 ; X2 ; SÞ  EI ; oXI

ð51Þ

where Malvern’s convention [39] is adopted for the gradient operator r. By making use of the relations, o/ðX1 ; X2 ; SÞ ¼ sa þ -;a pðSÞs3 ; oXa

ð52Þ

o/ðX1 ; X2 ; SÞ ¼ /00 þ Xa x  sa þ ð-pÞ0 s3 þ ð-pÞx  s3 ; oS

ð53Þ

where the symbols ðÞ0 and ðÞ;a denote F ¼ K þ ps3  r- þ

½ð/00

o oS

and

o oXa

for a ¼ 1; 2, the deformation gradient can be written as

 s3 Þ þ x  ð/  /0 Þ þ ð-pÞ0 s3   E3 0

¼ K½I þ pE3  r- þ fC þ j  KT ð/  /0 Þ þ ð-pÞ E3 g  E3 :

ð54Þ

2.6. Green–Lagrange strain The Green–Lagrange strain tensor, which is conjugate to the second Piola–Kirchhoff stress, can be expressed in terms of deformation gradient as E ¼ 12ðFT F  IÞ:

ð55Þ

Furthermore, the components of the Green–Lagrange strain tensor may be evaluated as EIJ ¼ EI  E  EJ :

ð56Þ

It should be noted that EI and EJ refer to the orthonormal basis vectors of the fixed material reference frame in the reference configuration, and E refers to the Green–Lagrange strain tensor. In the following we assume small strains but arbitrarily large displacements, rotations and twist. If the terms
, j1 , j2 and ð-pÞ0 are assumed to be of higher order than the term j3 , such that any products of the first four terms and any products of these terms with j3 can be neglected, then the strain components E13 , E23 , and E33 can be expressed as   1 oE13 ¼ p ; ð57Þ X2 j3 þ 2 oX1   1 oE23 ¼ X1 j3 þ p ; ð58Þ 2 oX2   1 2 0 2 2 E33 ¼
 X1 j2 þ X2 j1 þ ð-pÞ þ ðX1 þ X2 Þj3 : ð59Þ 2 In the above equation for E33 , the terms X1 j2 þ X2 j1 represent the normal strain due to bending, and the 0 term ð-pÞ represents the normal strain due to warping, and the term 12 ðX12 þ X22 Þj23 gives a normal strain due to twisting. This last term produces a coupling between the axial and twist behavior and is commonly referred to in the literature as the Wagner effect [10,40].

2518

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

Even though the above kinematic equations could be used to incorporate the shear strains due to bending, the assumption of /00 ¼ s3 made in deriving E13 and E23 restricts this extension. Hence, the above expressions for E13 and E23 represent only the shear strains due to uniform torsion. As noted in Section 2.1, the shear strains due to uniform torsion must be considered for accurate analysis of the inelastic torsionalflexural response of steel frame members within a distributed plasticity or plastic-zone approach. However, warping shear strains due to nonuniform torsion are neglected in this work as these shear strains tend to be an order of magnitude less than those due to uniform torsion in framing members of ordinary proportions. In this work, the warping amplitude p in Eqs. (57)–(59) is assumed equal to j3 . This assumption of p ¼ j3 is equivalent to the assumption of zero warping shear strains (cw ¼ ðj3  pÞ) due to nonuniform torsion [41]. In this work, the above expressions for E13 and E23 are further simplified to those given by the mitre model [10,30,31] for doubly-symmetric I-cross-sections. The mitre model has been shown [10,30,31] to provide an accurate representation of the shear strains due to uniform torsion, especially in the fully plasticregime. The elastic torque-twist relationship is predicted with high accuracy, and the inelastic relationship with reasonable accuracy. The model is also accurate in predicting the inelastic behavior of steel members when uniform torsion acts in combination with axial force and bending moment. Based on this model, the shear strain cp (¼ 2E23 or 2E13 ) due to uniform torsion is approximated as (see Fig. 5) cp ¼  j3 ;

ð60Þ

where  denotes the generalized coordinate associated with the shear strain due to uniform torsion. For a typical I-section, the function  can be expressed as [30,31]    ¼ 2t jX2 j  12 ðd  tf Þ sgnðX2 Þ; ð61Þ  ¼ 2tðjX1 j  12 ðbf  tf ÞÞsgnðX1 Þ

ð62Þ

bf

tw d

X1 X2 tf

bf

tf

Fig. 5. Geometric representation of mitre model.

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2519

along the long and short edges of the flanges, respectively, and ð63Þ

 ¼ 2tX1

for the web. In the above equations, d represents the depth of the cross-section, tf represents the thickness of the flanges and bf represents the width of the flanges. Also, sgnðÞ represents the sign of the argument ðÞ, and t is the elastic stiffness correction factor given by t¼

J 2bf tf3 3

þ

ðd 

2tf Þtw3 3

t3

ðd  2tf Þðtf  tw Þtw2 ðtf3  2w Þtf þ  4 6

;

ð64Þ

where tw is the thickness of the web and J is the St. Venant (uniform) torsional constant of the section [30,31]. Remark 2. It should be emphasized that the Eqs. (61)–(63) can be obtained directly from the Eqs. (57) and (58) simply by substituting the contour warping function associated with the doubly-symmetric I-sections [29,40]. In such a derivation, X1 and X2 refer to the cross-sectional coordinates of a material point. On the otherhand, if we invoke a vlasov type kinematics [27] for thin-walled beam cross-sections, then the Eqs. (61)–(63) can be derived from the Eqs. (57) and (58) by simply considering that the warping function - is
[29,41]. In such a
and the thickness warping function, x the sum of the contour warping function, x, derivation, it should be noted that X1 and X2 refer to the contour cross-sectional coordinates. In general, this second methodology is more accurate and hence extensively used for thin-walled cross-sections. For doubly-symmetric I-cross-sections composed of rectangular plates, both these methods are equivalent (see Ref. [29, pp. 93–95]). The above expressions for strains E33 and cp at a point can also be written in terms of the generalized coordinates and generalized strains at a cross-section as E ¼ Yk;

ð65Þ

where ½Y denotes the matrix of generalized coordinates and is given by 

1 ½Y ¼ 0

X1 0

X2 0

0

ðX12 þ X22 Þ 0

0 



k denotes the generalized strains at a cross-section, expressed as

k ¼
j2 j1 j03 12 j23 j3 :

ð66Þ

ð67Þ

By taking the variation of Eq. (65), we can write dE ¼ Ydk:

ð68Þ

Thus, the variation of strains at a material point is expressed in terms of the variation of generalized strains at a cross-section. Remark 3. As noted in Section 2.1, curved geometry and/or initial geometric imperfections are addressed within the formulation as an initial displacement effect. In simple terms, the strain at a point is taken as the difference in the Green–Lagrange strains associated with the current deformed geometry and the initial curved or imperfect geometry. Residual stresses within the initial deformed geometry are also specified.

2520

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2.7. Mechanical stress power The stress power may be expressed as Z }¼ P : F_ dV ;

ð69Þ

B

where P is the first Piola–Kirchhoff stress tensor, defined as P ¼ t0a  Ea þ t03  E3

ð70Þ

and P : F_ ¼ tr½PF_ T . Differentiating the Eq. (54), we obtain _ 3  E3  F_ ¼ WF þ K½C_ þ j_  KT ð/  /0 Þ þ -pE 0 _ þ pK½E 3  r- þ -ðj  E3 Þ  E3 þ - E3  E3 ;

ð71Þ

where W is the skew-symmetric spatial angular velocity tensor for which w the corresponding axial vector. Substituting Eqs. (70) and (71) into Eq. (69), and simplifying the result, we have Z o o _ dS }¼ ðn  c þ m  x þ Bp_ 0 þ N- pÞ ð72Þ L Z 0 _ dS; ðN  C_ þ M  j_ þ Bp_ 0 þ N- pÞ ¼ ð73Þ L0

where the spatial stress resultants that are conjugate to the spatial strains and curvatures are given by Z t03 dA0 ; ð74Þ n¼ A0

m¼

Z

ð/  /0 Þ  t03 dA0 ;

ð75Þ

A0

 Z B ¼ s3 

-t03 dA0

;

ð76Þ

A0

N - ¼ s3 

Z

ð-;a t0a þ -t03  x þ -0 t03 Þ dA0 :

ð77Þ

A0

The term n is the spatial cross-section force vector, m is the cross-section spatial moment vector, B is the o bimoment, and N- ¼ Tw , where Tw is the warping (nonuniform) torque. The symbol ð Þ denotes the corotating objective rate measured by an observer fixed in the moving frame sI , that is o

ð Þ ¼

o ðÞ  w  ðÞ: ot

ð78Þ

~, is equal to the second Based on the assumption of small strains, the rotated Cauchy stress tensor, r Piola–Kirchhoff stress tensor, S, up to the first-order (see the discussion in Section 5 of Ref. [8]). That is, ~ ¼ KT rK  S  KT P; r

ð79Þ

where r is the Cauchy stress tensor. Therefore, the material stress resultants that are work conjugate to the material strains and curvatures can be expressed as Z T ~3 dA0 ; N¼K n¼ ð80Þ r A0

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

M ¼ KT m ¼

Z

~3 dA0 ; KT ð/  /0 Þ  r

2521

ð81Þ

A0

 B¼

s3 

Z

-t03 dA0



 ¼

E3 

A0

N- ¼ E3 

Z

Z

-~ r3 dA0 ;

ð82Þ

A0

~a þ -~ ~3 Þ dA0 ; ð-;a r r3  j þ - 0 r

ð83Þ

A0

~3 ¼ r ~E3 , and A0 is the area of the cross-section in the reference geometry. It should be noted that where r the bimoment B is work-conjugate to p0 , which is taken as j03 , in this work. As mentioned in Section 2.1, beam-shear deformations and warping shear deformations are not considered in this work as they are expected to be negligible in framing members of ordinary proportions. Based on this assumption, and the assumption that the terms
, j1 , j2 and ð-pÞ0 are of higher order than the term j3 such that any products of the first four terms and any products of these terms with j3 can be neglected, then the mechanical stress power (Eq. (73)) may be expressed as Z }¼ ðk_  RÞ dS; ð84Þ L0 T

where R ¼ fN3 ; M2 ; M1 ; B; W ; Tsv g represents the generalized stress resultants at a cross-section, given by Z ~ dA0 R¼ YT r ð85Þ A0

and ½Y is given by Eq. (66). In the above description, N3 is the E3 component of the vector N, M1 and M2 are the components of M in the E1 and E2 directions, respectively. The symbols W and Tsv ¼ ðM3  Tw Þ refer to the Wagner term and the St. Venant torque, which are defined later. Substituting Eqs. (36) and (47) into the Eq. (72) and using the Eq. (78) to obtain the corotating objective rates, the governing equations of equilibrium in the spatial system can be stated as [8,29,42] 9 on > > þ
n¼0 > > oS > > = om o/0 ð86Þ
¼0 ; þ nþm > oS oS > > > > oB > ;  N- þ B ¼ 0 oS
, m
and B represent the applied spatial distributed force and moment vectors and the bimoment per where n unit length, respectively. These equations, and Eqs. (80)–(82) are utilized in Section 3.1 to develop equilibrium based functions for interpolation of the element generalized stresses. 2.8. Governing equations: weak form Consider a beam ðB
R3 Þ that occupies a volume V0 ¼ A0  L0 with the boundary S0 in its reference configuration. Let Sr0 and Su0 denote the segments of S0 over which the tractions and displacements are prescribed such that S0 ¼ Sr0 þ Su0 . The deformation of the body from its reference configuration is described by the displacement field u : B ! R3 , which is assumed to be prescribed as u ¼
u on Su0 . Further, assume that the nominal traction vector is prescribed as t0 ¼
t0 on Sr0 . Let U and G denote the spaces of the generalized displacement (e.g., displacements, rotations, and warping degree of freedom) and generalized stress (i.e., stress resultant) fields, respectively. Similarly, let Ud

2522

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

and Gd denote the spaces of the admissible generalized displacement and generalized stress variations. The displacements are assumed to be smooth within the element and compatible between adjacent elements, satisfying the essential (kinematic) boundary conditions on Su0 . Since the beam kinematics are those of Bernoulli–Euler beam theory with superposed torsional-warping of the cross-sections, the generalized displacements and the admissible displacement variations belong to the following spaces: ) n u and S 2 L0 g U ¼ fu 2 ½H 2 ðSÞ dim jujSu ¼
0 ; ð87Þ n Ud ¼ fdu 2 ½H 2 ðSÞ dim jdujSu ¼ 0 and S 2 L0 g 0

where ndim ¼ 3 for the displacement vector and the rotation vector, and ndim ¼ 1 for the torsional warping amplitude, and H 2 ðSÞ denotes the Sobolev space of functions [43]. The generalized stress variation is continuous within the element, but is not interelement compatible for the formulation presented in this paper. Thus, the physical and admissible variations of generalized stress belong to the square integrable function spaces denoted by G ¼ Gd ¼ fsjs 2 ½L2 ðSÞnstr g;

ð88Þ

where nstr ¼ 6 for the current beam formulation (axial force, two components of moment, bimoment, uniform torque and Wagner force). Given the above definitions, the two-field Hellinger–Reissner HR(u; R) variational principle can be expressed in terms the of generalized quantities as Z Z T ^ ^ dS ¼ 0 HRðu; R; du; dRÞ ¼ dk R dS  Gext  dRT ðk  kÞ ð89Þ L0

L0

^ in Eq. (89) represents the compatible for all variations ðdu; dRÞ 2 Ud  Gd , and ðu; RÞ 2 U  G. The term k ^ generalized strain field (derived from the displacements), and dk is the variation in this strain. Furthermore, k represents the generalized strain field that is derived from the generalized stress field (R). The term Gext in Eq. (89) is defined as Z Z 0
Gext ¼ q0 b0  du dV0 : ð90Þ t  du dS þ Sr0

V0

For a beam element, this expression can be further reduced to Z
 dw þ BdpÞ dS: Gext ¼ ð
n  d/0 þ m

ð91Þ

L0

^ leads to Linearization of the HR(u; R; du; dR) functional (Eq. (89)) about a configuration (^u; R) b du; dRÞ þ DHRð^u; R; b du; dRÞ; L½HRðu; R; du; dRÞ ¼ HRð^ u; R; where b du; dRÞ ¼ DHRð^ u; R;

Z

ð92Þ

^  DR þ Ddk ^  RÞ dS ðdk

L0



Z

^ þ DdR  ðk  kÞÞ ^ dS  DGext : ðdR  ðDk  DkÞ

ð93Þ

L0

2.9. Two-space von-Mises constitutive representation Based on the kinematic assumptions outlined in the previous subsections, the state of stress at a general location within the element cross-section includes the axial stress due to axial force, biaxial bending and

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2523

bimoment, and the shear stress due to uniform torsion. For the examples presented in this paper, a simple J2 material with multi-linear kinematic hardening is assumed. Reduction of the corresponding six-space constitutive relations onto the required two-dimensional stress-constrained subspace imposes nontrivial constraints on the return mapping algorithms, since the stresses are constrained to remain within the stressconstrained subspace. In this work, we adopt the strategy proposed for plane stress problems by Simo and Taylor [44] to derive the return mapping algorithm and consistent tangent modular matrix for the twodimensional stress-constrained subspace defined by the beam kinematics. Accordingly, if the projection matrix P of Ref. [44] is replaced with   1 2 0 ð94Þ P¼ 3 0 6 the appropriate two-space constitutive equations are obtained. The reader is referred to [44] for a detailed description of the return mapping algorithm and the consistent tangent formulation for stress-constrained subspaces.

3. Finite element discretization Snelm Based on the finite element discretization V0 ¼ e¼1 V0e and the finite-dimensional approximation subspaces for the generalized displacements and the generalized stresses, such that 8ðu; RÞ 2 U  G, we have  Uh ¼ fuh 2 H 2 ðSÞjuh ¼ N rg ; ð95Þ Gh ¼ fRh 2 L2 ðSÞjRh ¼ Pbg where Uh
U and Gh
G. The symbol h denotes the association of Uh and Uh with the discretization of the domain, where h is the characteristic element length scale. In the present study, cubic Hermitian ^ The cubic fields for u and v are interpolation is used for all the displacement fields (u; v; w, and /). necessary to capture the bending response of end-loaded beams with good accuracy. The cubic field for w is necessary for these displacements to be represented with comparable accuracy to that for u and v within the total-Lagrangian description, and to avoid locking behavior in the thin-membrane limit. In a mixed-formulation, independent interpolations are chosen for both the primary and dual field variables, while in displacement-based formulations, the dual variable is obtained a posteriori by differentiation. The choice of the interpolation fields plays an important role in accurately representing the behavior of constrained problems (e.g., incompressible flow), and in achieving good convergence characteristics. It should be noted that both mixed and displacement-based formulations lead to identical solutions if the interpolations for the primary and dual variables are chosen such that the limitation principles are satisfied [45–47]. In addition, the order of the interpolation functions for the dual variables should be selected such that the element accuracy is improved in cases where the displacement interpolation is not exact or may not correspond to the element deformation (e.g., locking problems). Also, the interpolation functions, P, must be of sufficient order (to satisfy patch test) and the interpolation space must satisfy certain stability constraints such that the element will not have spurious zero-energy modes (i.e., the element stiffness matrix is of sufficient rank) [48]. In the present study, the selected generalized stress interpolation is dependent on the displacements, such that the geometric nonlinear effects are captured more accurately within each element. Thus, the generalized stress interpolations are nonlinear functions of both the dual (force) and primary (displacement) variables as explained below.

2524

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

In a two-dimensional end-loaded beam, the variation of bending moment along the length of the member is composed of a primary moment, which varies linearly, and a secondary moment (typically referred to in structural engineering as the P  d moment), associated with the geometric nonlinearity. Therefore, if this P  d moment is accounted for within the generalized stress interpolation, the interpolation matrix P does not belong to a linear vector space. It depends on the displacements in the current deformed configuration. An important result of such an interpolation is that the variations and the physical fields for the generalized stresses do not belong to the same approximation subspaces (see Eq. (96) below). However, in the limit that the geometric nonlinearity vanishes, i.e., in the geometrically linear case, Pb reduces to a linear function of the dual variables b. In such instances, the actual fields (for end loaded elements) and the variations belong to the same approximation subspaces. The approximation subspaces for the variations in the displacements and the generalized stresses may be expressed as Uhd ¼ fduh 2 H 2 ðSÞjduh ¼ N drg Ghd ¼ fdRh 2 L2 ðSÞjdRh ¼ Pdb þ dPbg

) :

ð96Þ

3.1. Generalized stress interpolation functions, P In this section, the generalized stress interpolations are developed based on the governing equations of equilibrium. Although, it is not necessary that the generalized stress interpolation functions P satisfy the equilibrium equations section-by-section along the length of the element, better results are obtained if equilibrium is satisfied. In the beam element developed in this research, compatibility of the force fields is not enforced at the element boundaries (i.e., at the nodes). Correspondingly, the dual variables are condensed out at the element level. The governing equations of equilibrium in terms of generalized stresses are given by Eq. (86) in the fixed (spatial) reference frame. By integrating the first of these equations, we can write n ¼ ðnA þ
nS0 Þ;

ð97Þ

where
nS0

¼

Z

S


n dS

ð98Þ

0

and nA is one of the generalized stress vectors at S ¼ 0. The components of nA are expressed as T nA ¼ fb1 ; b2 ; b3 g . Since cubic Hermitian interpolation is used for the axial displacement in this work, the axial force variation along the length predicted by the displacement interpolation is quadratic for elastic element behavior. Therefore, the corresponding force term in R (i.e., n3 of n) must vary quadratically along the length to avoid spurious modes in the element stiffness matrix. If we assume the following member loading
n ¼ ðb4 þ 2b5 SÞe3

ð99Þ

then
nS0

¼

Z 0

S


n dS ¼ ðb4 S þ b5 S 2 Þe3 :

ð100Þ

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2525

Thus, we obtain

9 n1 ¼ b1 = n2 ¼ b2 ; n3 ¼ b3  b4 S  b5 S 2

ð101Þ

as the components of the spatial stress resultant (force) vector n in the fixed (spatial) reference directions. Similarly, by integrating the moment equilibrium equation (Eq. (86)), one obtains
Sn0  m
S0 ; m ¼ mA þ ð/0  nA Þ þ m

ð102Þ

where mA is the cross-section moment vector at S ¼ 0,   Z S S S
n0 ¼ /0 
m ð/0 
nÞ dS n0 

ð103Þ

0

and
S0 m

¼

Z

S

ð104Þ


dS: m 0 T

The components of mA are expressed as mA ¼ fb6 ; b7 ; b8 g . Also, the term /0 in Eq. (103) is defined by Eq. (3). Since a cubic field is chosen for the twist interpolation in this work, the Saint-Venant torque predicted by the assumed displacement field is quadratic within an elastic element. As a result, the moment term m3 in R must vary quadratically in S for the element stiffness to have sufficient rank. If we assume the member loading,
¼ ðb9 þ 2b10 SÞe3 m then
S0 ¼ m

Z

ð105Þ

S


dS ¼ ðb9 S þ b10 S 2 Þe3 : m

ð106Þ

0

Finally, after defining the following quantities 9 RS Iu ¼ 0 u dS > > > RS > Isu ¼ 0 uS dS = RS Iv ¼ 0 v dS > > > > RS ; Isv ¼ 0 vS dS

ð107Þ

we can write the components of the spatial stress couple (i.e., moment) vector m in the fixed (spatial) reference directions as 9 m1 ¼ b6 þ ðS þ wÞn2  vn3  b4 Iv  2b5 Isv = ð108Þ m2 ¼ b7  ðS þ wÞn1 þ un3 þ b4 Iu þ 2b5 Isu : ; m3 ¼ b8 þ vn1  un2  b9 S  b10 S 2 It should be emphasized that a quadratic function is used for both the axial force and Saint-Venant torque interpolations to eliminate spurious zero-energy modes even though a constant term would be sufficient to satisfy equilibrium for a straight axially loaded member in the absence of external element loading. This is required to give the element stiffness matrix sufficient rank. The constitutive law must be based on material components of the generalized stresses. Thus, it is necessary to transform the above spatial generalized stress interpolation into a material generalized stress interpolation. This transformation is given by Eq. (86) for the total force and by

2526

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

DN ¼ KT Dn þ DKT n ¼ KT Dn  DHKT n ¼ KT Dn þ ½NDh

ð109Þ

for the increment in the force, where ½N is the skew-symmetric operation for which N ¼ fN1 ; N2 ; N3 gT is the corresponding axial vector. Similarly, the increment in the material stress couple is given by DM ¼ KT Dm þ ½MDh;

ð110Þ T

where the ½M is the skew-symmetric operation for which M ¼ fM1 ; M2 ; M3 g is the corresponding axial vector. The interpolation of the bimoment B along the length of the element is assumed to be the same form as the exact elastic, geometrically linear variation of the bimoment within a prismatic elastic member, with X1 and X2 taken as the principal axes of the cross-section. This results in improved accuracy compared to the linear variation obtained from the displacement interpolation, and has been suggested by Aalberg [49] as providing an accurate representation of the inelastic variation of the bimoment along the length of inelastic members. Thus B is expressed as B ¼ Nw1 b11 þ Nw2 b12 ;

ð111Þ

where b11 and b12 are the element end values of B at S ¼ 0 and L0 , respectively,     kz cosh k kz sinh Nw1 ¼  cosh  ; ð112Þ l sinh k l    sinh kzl Nw2 ¼  ð113Þ sinh k qffiffiffiffiffiffi GJ , where G and E are the elastic shear and Young’s moduli of the material, J the St. Venant and k ¼ L0 EC w torison constant of the cross-section, and Cw the cross-section warping constant. As a result, the interpolation for the warping torque Tw is obtained as Tw ¼

dB ¼ Nw0 1 b11 þ Nw0 2 b12 : dS

ð114Þ

Given Tw , the St. Venant torque is calculated as Tsv ¼ M3  Tw . Finally, a constant interpolation is used for the Wagner stress resultant (W ). W ¼ b13 :

ð115Þ

In summary, we obtain the interpolation functions for the material generalized stresses by combining Eqs. (101), (108), (111), (114) and (115) with Eqs. (80)–(82). The resulting equation is represented symbolically by Eq. (95), where b is the 13 · 1 vector of the generalized stress parameters and R is the 6 · 1 vector of the generalized stresses given by Eq. (85). 3.2. Incremental equations for generalized stress and strain Following standard procedures, the incremental generalized strains can be expressed as ^ ¼ BDr; Dk

ð116Þ

where B is the discretized strain-displacement operator obtained from Eqs. (38) and (48). Also, the incremental material generalized stresses are obtained as DR ¼ PDb þ WDr

ð117Þ

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2527

and the corresponding generalized strains are given by Dk ¼ D1 cp DR;

ð118Þ

where WDr ¼ DPb

ð119Þ

and Z

Dcp ¼

ST Ccp S dA:

ð120Þ

A

In the above equation, Ccp represents the material consistent tangent matrix at each Gauss point within the cross-section, while Dcp is the cross-section consistent tangent modular matrix. Also, in Eq. (119), DP is a function of incremental nodal displacements Dr. In Eq. (119), DPb is rearranged into the matrix W containing only total generalized stress terms multiplied by the vector of the incremental nodal displacements Dr. 3.3. Discretized form of HR(u; R) The discrete form of the HR(u; R) principle leads to the following system of nonlinear equations: g¼

Z n[ elm L0

e¼1

¼

n[ elm

BT R dS 

fFeint  Feext g ¼

e¼1

V¼

Z

Z

^ dS  Fe WT ½k  k ext

 ¼0

L0 n[ elm e¼1

ge ¼ Fint  Fext ¼ 0

ð121Þ

^  kÞ dS ¼ 0 PT ðk

L0

U ¼ ðRR  RÞ ¼ 0: The first of the above equations (g ¼ 0) represents the global discretized equilibrium equations, the second of these equations (V ¼ 0) represents the (local) intra-element compatibility between the strains derived from the displacement fields and the strains derived from the interpolation of the element generalized stresses, and the third of these equations (U ¼ 0) represents the (local) cross-section constitutive law. Although we use boldface capital letters V and U for local element compatibility and cross-sectional constitutive relations, respectively, these variables indeed denote vector quantities. This is an exception toSthe nelm convention of using boldface capital letters for second-order tensor or matrix quantities. In Eq. (121), e¼1 e e denotes the finite element assembly operator, Fint is the element internal force vector, and Fext is the element contribution of the external force vector, including any directly applied nodal loads. The variables Fint and Fext denoteSthe global internal and external force vectors, respectively. It should be noted that the assembly elm operator ne¼1 implicitly handles the mapping between the element nodal displacement and the global displacement degrees of freedom. The above nonlinear system of equations (Eq. (121)) are solved using Newton iteration. Since interelement compatibility is not enforced in the generalized stress variables, the nonlinear discretized straindisplacement compatibility equation V ¼ 0 can be solved iteratively at the element level for every global

2528

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

iteration. Similarly, the nonlinear constitutive equation U ¼ 0 can be solved iteratively at the crosssection level for every element level iteration. Alternatively, linearized equations may be employed at the section and/or the element levels, and the resulting section and/or element residuals may be transferred to the equations at the higher level or levels. The reader is referred to [26,50] for more detailed discussion of the consistent linearization of the element equations, and the development of the corresponding state determination procedures. 3.4. Tangent stiffness matrix Consider the incremental form of HR(u; R) (see Eq. (93)). After substituting the discretized strain-displacement operator and the interpolation functions, this equation can be expressed as 9 T > fdrg fGT fDbg þ GTw fDrg þ Kg fDrgg > > > T T T fdrg f½Hpw  fDbg þ Hw fDrg  Gw fDrg þ Mk fDbgg = ¼ 0; ð122Þ T > fdbg fHfDbg þ Hpw fDrg  GfDrg þ Mk fDrgg > > > ; T T T fdrg fDfg  fdbg fc1  c2 g  fdrg fc3  c4 g where G¼

R

Gw ¼

L0

PT B dS

R

L0

WT B dS

9 > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > =

R Hpw ¼ L0 PT D1 cp W dS R T 1 Hw ¼ L0 W Dcp W dS R H ¼ L0 PT D1 cp P dS R ^T R dS fdrgT Kg fDrg ¼ L0 Ddk : R ^ dS > > fdbgT Mk fDrg ¼ L0 DdRT fk  kg > > > R > > > c1 ¼ L0 PT k dS > > > > R > > T^ > c2 ¼ L0 P k dS > > > > R > T > > c3 ¼ L0 W k dS > > > > R > T^ > > c4 ¼ L0 W k dS > > > > R ; e T fDfg ¼ Fext  L0 B R dS

ð123Þ

For Eq. (122) to be valid for any virtual set of fdrg and fdbg, the coefficients of fdrg and fdbg in this equation must be equal to zero. In matrix form, the resulting linearized equations can be expressed as   T   Dr Df  ðc3  c4 Þ ðGw þ Gw  Hw þ Kg Þ ðG  Hpw  Mk ÞT ð124Þ ¼ ðc1  c2 Þ Db ðG  Hpw  Mk Þ H from which fDbg can be obtained as fDbg ¼ H1 fðG  Hpw  Mk ÞfDrg  ðc1  c2 Þg:

ð125Þ

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2529

Substituting Eq. (125) into Eq. (124), we have ½Ke fDrg ¼ fDpg;

ð126Þ

where ½Ke  ¼ ½G  Hpw  Mk T H1 ½G  Hpw  Mk  þ ½GTw þ Gw  Hw þ Kg 

ð127Þ

fDpg ¼ fDfg þ fc3  c4 g þ ½G  Hpw  Mk T H1 fc1  c2 g:

ð128Þ

and

3.5. Joint compatibility within structural framing systems At a rigid joint that undergoes a rotation in three-dimensions, the orientations of the moving reference axes are one and the same for all the members (elements) meeting at that joint. Thus, the three independent quantities that are used in parametrizing the finite rotation transformation of the joint coordinate frame are ^ T . It should be ^ ^c; /g continuous across the joint. In the present formulation, these quantities are v ¼ f#; noted that the quantity d
used in transforming fdag to fdhg (Eq. (24)) is not continuous at a joint. Hence the displacement derivatives that are used in the earlier sections as state variables are not in general compatible at a joint. Due to the above factor, the state variables fdrg need to be transformed to an alternate state variable set fdrh g before assembling the global stiffness. This transformation may be written as fdrg ¼ ½T1 drh ! dr fdrh g;

ð129Þ

where ( fdrg ¼

du1

( fdrh g ¼ and

du1

2

dv1 du2 dv1 du2

I33 6 043 6 6 013 1 ½Tdrh ! dr  ¼ 6 6 033 6 4 043 013

dw1 dv2

du01 dw2

dw1 dv2

dv01 du02

dw01 dv02

^ d/ 1 dw02

dhx1 dw2

dhy1 dhx2

dhz1 dhy2

T1 da ! dhð1Þ ð44Þ 014 034 044 014

031 041 1 031 041 0

033 043 013 I33 043 013

034

d
1 dhz2

^0 d/ 1 ^ d/ 2

)T ^0 d/ 2

^0 d/ 1 d
2

034 044 014 034

T1 da ! dhð2Þ 44 014

;

ð130Þ

161

)T ^0 d/ 2

ð131Þ 161

3 031 041 7 7 0 7 7: 031 7 7 041 5

ð132Þ

1

In the above equation, Tda ! dhð1Þ and Tda ! dhð2Þ are evaluated at nodes 1 and 2, respectively, using the Eqs. (24)–(28). Based on Eq. (129), the transformed stiffness is obtained as 1 H ½Krh  ¼ ½TT drh ! dr ½Ke ½Tdrh ! dr  þ ½H ;

ð133Þ

where ½HI fDrh g ¼ ½DTT drh ! dr fQg

ð134Þ

2530

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

and fQg is the set of work conjugate forces to the dof set fDrg. Also, it should be noted that the degrees of freedom d
1 ; d
2 that appear in Eq. (131) are statically condensed out at the element level prior to assembly of the global equations. In handling the warping degrees of freedom at beam-to-column joints in structural frameworks, the following simplifying assumptions are employed in this work. If a physical member or members are con^ 0 is imposed between the elements representing the member on tinuous through a joint, then continuity of d/ each side of the joint. However, these derivatives (warping deformations) are assumed to be uncoupled

24

120

E = 7.2 × 10 A=6 I=2

6

120 Fig. 6. Snap through of a hinged right-angle frame.

2.5

4 x 10

2.5

1.5

1.5

1

1

Load P

2

Load P

2

0.5

0.5 0

0 –0.5

–0.5

___ Displacement method _._. Mixed method

–1

___ Displacement method _._. Mixed method

–1 –1.5

–1.5 0

(a)

4 x 10

10

20

30

40

50

60

70

Vertical displacement

80

90

100

0

(b)

10

20

30

40

50

60

70

80

90

100

Horizontal displacement

Fig. 7. Snap through of a hinged right-angle frame: (a) load versus vertical displacement under applied load; (b) load versus horizontal displacement under applied load.

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2531

between any members that meet at an angle other than 0 or 180 at a structural joint. This assumption is similar that adopted by Yang and McGuire [14,15]. 3.6. State determination algorithm In general, state determination algorithms for mixed finite element formulations involving material and geometric nonlinearity are more complex than those in displacement-based formulations. This is due to the presence of more than one independent field variable in the variational principle. Also, the structure of the constitutive theory plays an important role as many constitutive theories consider strain as an independent variable and stress as a dependent variable. In a standard displacement-based formulation, the state determination is a strain driven problem, i.e., the stresses are derived from the strains. This is in sharp contrast with stress-based mixed finite element formulations where the displacements and stresses are independent variables. In the examples presented in this paper, the nonlinear system of equations at the global, element and section levels (Eq. (121)) are each solved separately using Newton’s method. The local element

P 10

E = 6.895 x 10 N/m 2 10

G = 2.652 x 10

N/m2

-4

2 A = 6.452 x 10 m

R

α

-8

4 I = 4.162 x 10 m

R = 2.54 m o

α = 215

10

___ Displacement method _._ Mixed method

EI

Normalized Load

PR 2

8

v/R 6

w/R 4

2

0

–2 0

0.2

0.4

0.6

0.8

1

1.2

Normalized displacements

Fig. 8. Clamped-hinged deep circular arch subjected to point load.

1.4

2532

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

Table 1 Normalized buckling load versus analysis method Model

Normalized buckling load

Present (displacement element) Present (mixed element) DaDeppo and Schmidt [56] Wood and Zienkiewicz [57] Simo and Vu-Quoc [5] Borri and Bottasso [58] Franchi and Montelaghi [60]

9.00 8.966 8.97 9.24 9.053 9.07 8.984

PR2 EI

Z 1 00

R

=1

1 X

45

E = 10

o

P

7

R = 100

Y

Fig. 9. Initial geometry of 45 cantilever bend.

Table 2 Tip geometry (x; y; z) for 45 bend, initially ðx; y; zÞ ¼ ð29:29; 0:0; 70:71Þ Model

Load level 300

450

Present (displacement element) Present (mixed element) Bathe and Bolourchi [3] Simo and Vu-Quoc [5] Cardona and Geradin [6] Crisfield [7] Ibrahimbegovic et al. [12] Avello [61] Borri and Bottasso [58] Franchi and Montelaghi [60]

22.28, 40.09, 58.81 22.28, 40.11, 58.79 22.5, 39.5, 59.2 22.33, 40.08, 58.84 22.14, 40.35, 58.64 22.16, 40.53, 58.53 – 22.14, 40.65, 58.66 22.38, 39.33, 59.24 22.33, 40.11, 58.86

18.58, 18.56, – 18.62, 18.38, 18.43, – 18.23, 18.68, 18.62,

600 48.37, 52.30 48.41, 52.25 48.39, 52.32 48.59, 52.11 48.79, 51.93 49.31, 51.84 47.51, 52.88 48.43, 52.32

15.78, 53.33, 47.23 15.74, 53.38, 47.16 15.9, 53.4, 47.2 15.79, 53.37, 47.23 15.55, 53.50, 47.04 15.61, 53.71, 46.84 15.62, 53.50, 47.01 15.26, 54.54, 46.48 15.88, 52.37, 47.93 15.80, 53.40, 47.24

strain-displacement compatibility equation V ¼ 0 is solved iteratively at the element level for every global iteration. Similarly, the local cross-section constitutive equation U ¼ 0 is solved iteratively at the crosssection level for every element level iteration. In the companion paper [26], this algorithm is referred to as the N–N algorithm. The reader is referred to [26] for a detailed discussion of four alternative variationally consistent state determination algorithms and their relative efficiency for elastoplastic analysis using mixed finite element formulations.

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2533

Table 3 Convergence rate for the 45 bend Iteration

Residual norm

2nd load step 0 1 2 3 4 5 6 7 8 9

1.000 · 102 5.229 · 105 1.489 · 104 2.549 · 104 4.510 · 102 5.606 · 104 1.704 · 102 9.360 · 101 5.729 · 103 2.485 · 106

6th load step 0 1 2 3 4 5 6

1.000 · 102 4.666 · 104 3.899 · 102 9.945 · 101 3.256 · 101 3.745 · 103 9.091 · 106

30

0.6 X

240

Z Y

240

E = 71240

P

Fig. 10. Initial geometry of right-angle frame.

4. Numerical examples In the following, numerical examples that verify the correctness of the proposed formulation (through a number of two- and three-dimensional elastic benchmark problems), and which verify its accuracy in

2534

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

Z

X

Z 1.43 1.43 X

Y

(b)

(a)

Applied Load (lb)

1. 5

1

___ Displacement method _._ Mixed method 0. 5

o

Crisfield’s solution

0 0

10

20

30 40 50 Lateral tip displacement (in )

60

70

(c)

Fig. 11. Lateral buckling of a cantilever right-angle frame: (a) perspective view of final deformed shape; (b) projection of final deformed shape onto the X –Z plane; (c) applied load versus lateral tip displacement of the free end.

representing three-dimensional inelastic stability behavior are presented. The solutions obtained from the mixed-formulation are also compared to those obtained from a displacement-based formulation having the same beam kinematics. The reader is referred to [50] for details on the formulation of displacementbased beam element. Unless noted otherwise, in the following examples, full Newton iteration is used at the global level in conjunction with a minimum residual displacement iterative strategy [51] and Crisfield’s automatic arclength incrementation procedure [52] to trace the equilibrium path. At the structural level, convergence of the finite element solution is based on the Euclidean norm of the out-of-balance force vector as well as the incremental displacement vector. Similarly, the Euclidean norm of the incremental force parameters and of the unbalanced generalized stresses are used for the convergence criteria at element and section levels, respectively. Convergence tolerances equal to 103 , 108 , and 1012 are used at structural, element and section levels, respectively.

4.1. Buckling of Lee’s hinged right-angle frame This example is concerned with the loss of stability and the general large deflection response of the rightangle frame shown in Fig. 6 [5,53,54]. Ten elements, five in each member, are used to model the frame. The

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2535

6.11"

P 0.33" X

6.45"

L ey ex

E = 30000 ksi

Y W 6x25 Top: ex = 1.67" ey = 2.95"

Est = 100 ksi

σy = 36 ksi

0.45" Bottom: e x = 1.66" e y = 2.95"

90 80 70

(- u)

Load P (kips)

60 50

___ Displacement method (2 elements) 40

_._ Mixed method (2 elements) o Birnstiel test data (specimen #4)

30 20 10 0 0

0.5

1

1.5

2

2.5

Out-of-plane displacement (inches)

Fig. 12. Birnstiel’s test no. 4: (a) test setup; (b) cross-section dimensions; (c) applied axial load versus mid-span out-of-plane displacement.

90

90

80

80 70

(- v)

60 50

___ Displacement method (2 elements)

40

_._ Mixed method (2 elements)

30

60 50 40

_._ Mixed method (2 elements)

20

10

10 0.1

0.2

0.3 0.4 0.5 0.6 Inplane displacement (inches)

___ Displacement method (2 elements)

30

o Birnstiel test data (specimen #4)

20

0 0

(a)

Load P (kips)

Load P (kips)

70

0.7

0 0

0.8

(b)

0.005

0.01 0.015 0.02 Central twist rotation (rad)

0.025

Fig. 13. Birnstiel’s test no. 4: (a) applied axial load versus mid-span in-plane displacement; (b) applied axial load versus central twist rotation.

2536

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

T0

N0

u = 0 v = 0

140.45 mm

H 7.3 mm

Z, w

2030 mm

2090 mm

140.9 mm

φz = 0 u = 0

11.4 mm

Y, v

HEB 140

u = 0 X, u

T0

v = 0 w = 0

E = 210

3

× 10 N/mm2

E st = 500 N/mm2

σy = 279 N/mm2 Fig. 14. Aalberg’s tests, load and support conditions.

initial load increment and the target number of iterations are chosen as 2000 and 5, respectively. The load– deflection curves for the entire post-buckling range are shown in Fig. 7. Using the displacement-based finite element, the entire post-buckling path was traced in 54 load increments. On the other hand, the mixed based finite element traced the entire post-buckling path in 40 load increments. The buckling load for the frame is obtained as 18,616.

4.2. Clamped-hinged deep circular arch This example (see Fig. 8) has been considered by a number of authors [5,55–60] and the exact solution based on the Kirchoff–Love theory is given by DaDeppo and Schmidt [56]. The finite element mesh consists of 40 elements. A target number of five iterations is chosen for the problem. The solution is obtained in 53 load increments using the mixed element and in 73 load increments using the displacement-based element. Convergence to within 103 tolerance is achieved at an average of six iterations per load increment. The post-buckling load–displacement curve is shown in Fig. 8. 2 Table 1 compares the present solution for the normalized buckling load PR with that given by other EI authors.

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2537

HEB 140 (a) 6

Torsional moment (kNm)

5

4

3

___ Experimental (Aalberg) _ _ Displacement method _._ Mixed method

2

1

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

Twist rotation (rad)

(b) Fig. 15. Uniform torsion test on HEB 140 beam: (a) test setup; (b) torsional moment versus twist rotation.

4.3. Forty-five degree cantilever bend This example presents the geometrically nonlinear analysis of a cantilever bend placed in the horizontal plane with a vertical load applied at the free end as shown in Fig. 9. Since this example involves large threedimensional rotations and an initially curved geometry, it has been considered by a number of authors [3,5– 7,11–12,58–61]. Eight elements are used. Table 2 compares the present solution with those given by number of authors for the position of the tip at loads of 300, 450 and 600. Solutions are obtained by using eight equal load increments of 75. This problem is also solved by using six equal load increments of 100 by a number of authors. An average of seven iterations per load increment are observed when six equal load increments are applied. The residual norm indicating the convergence rate with number of iterations is shown in Table 3. 4.4. Cantilever right-angle frame subjected to end load This example also involves a full three-dimensional response and considers the buckling and postbuckling response of a right-angle frame when subject to transverse load P at the free end (see Fig. 10). The frame is subjected to an in-plane end load in the z direction. A small perturbation load of 2 · 104 is applied at the free end normal to the plane of the frame. The frame is modeled using ten elements, with five elements in each leg. The present solution was obtained using the arc-length method and the calculation was stopped after 25 load steps. The buckling load is found to be 1.086 and load versus out-of-plane deflection curves obtained with the displacement-based and mixed based element are compared with Crisfield’s [7] solution in Fig. 11.

2538

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

HEB 140 (a) 7

Torsional moment (kNm)

6

5

4

___ Experimental (Aalberg) 3

_ _ Displacement method (2 elements) _._ Mixed method (2 elements)

2

1

0 0

0.2

0.4

0.6

0.8

1

1.2

Twist rotation (rad)

(b) Fig. 16. Nonuniform torsion test on HEB 140 beam: (a) test setup; (b) torsional moment versus twist rotation.

4.5. Birnstiel’s tests In this example, the finite element solutions are compared to the Birnstiel’s [40] test data for Specimen No. 4, a beam-column subjected to biaxial bending from an eccentric axial load. The test setup and the cross-sectional dimensions of the beam-column are shown in Fig. 12. The beam-column specimen is assumed to have initial sinusoidal imperfections with amplitudes of 0.01 in the transverse direction, and a twist imperfection of magnitude 0.004 rad. The Galambos and Ketter residual stress pattern [62] is assumed. Horizontal displacements and twist rotations are restrained at the top and bottom of beam-column. Warping is also assumed to be restrained at the top and bottom of the member. Two elements are used in the finite element analysis. The load versus mid-span out-of-plane displacement curve is shown in Fig. 12. Also, the load versus mid-span in-plane displacement and load versus central twist rotations are shown in Fig. 13. Good agreement is observed between the test data and the finite element solutions. 4.6. Aalberg’s tests In this subsection, finite element solutions are compared to test results obtained by Aalberg [49] with a variety of loading conditions. The test setup, loads and support conditions are shown in Fig. 14. The beam-columns are assumed to be initially straight and to have zero residual stresses. For the uniform and nonuniform torsion cases, one element (either displacement-based or mixed) is used to model ^ and the out-of-plane displacement u are restrained half of the length of the member. The twist rotation / z

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2539

at the middle of the beam-column (see Fig. 14). The following tests are considered on this member geometry: • Uniform and nonuniform torsion tests: The load–deflection curves for uniform and nonuniform torsion cases are shown in Figs. 15 and 16, respectively. From the figures, it can be seen that the mitre model has the ability to account for the flange tip shear stresses due to uniform torsion. • Combined bending and torsion tests: The finite element solutions are compared to the results for test H-2-MT [49] in this case. The loads and support conditions are shown in Fig. 17. Four elements are used to model the full length of the beamcolumn. The applied loads and displacements are normalized with respect to their yield values as shown below: N¼

N Ny


¼ w

M¼

M My


v ¼

T ¼

T Ty

9 w> > > wy > > > > > = v >

ð135Þ

vy > ; > > > > ^ >
^ / > > > /¼ ; ^ /y

T

T H HEB 140 (a)

1.6

1.4

T

φ

Normalized load

1.2

v

M 1

φ

M

0.8

___ Displacement method (4 elements) 0.6

_._ Mixed method (4 elements) +, o, x Experimental data (Aalberg)

0.4

0.2

0

0

1

2

3

4

5

6

Normalized displacement s

(b) Fig. 17. Bending and torsion (MT) combination: (a) test setup; (b) normalized load versus displacement curve.

2540

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

where the yield values for HEB 140 specimen are given by 9 Ny ¼ 1198 kN wy ¼ 2:85 mm > = My ¼ 58:9 kN m vy ¼ 7:3 mm : ; ^ ¼ 4:95 > Ty ¼ 2:65 kN m / y

ð136Þ

The normalized load versus normalized displacement curves are plotted in Fig. 17. It is seen that the finite element solutions agree closely with the test results. • Combined axial load and torsion test: The finite element solutions are compared with the results for test H-9-NT [49] in this study. Two elements are used to model the full length of the beam-column. The normalized torque versus normalized rotation curves are shown in Fig. 18. Again, the finite element solutions agree closely with the test results. 4.7. Portal sway frame in strong-axis bending This frame is proportioned to represent a realistic nonredundant framing system. As shown in Fig. 19, column AB is the only member of the lateral resisting system, and is rotationally restrained both at its top and bottom, thus simulating typical reversed-bending curvature conditions in a multi-story frame. Both two- and three-dimensional plastic-zone analyses are presented for the frame. An initial out-of-plumbness L L of 500 as well as an initial out-of-straightness of 1000 in the plane of the frame are assumed. The sinusoidal out-of-straightness is placed to the right in the form of a half-sine wave. For three-dimensional analysis, a

T

N

N

T

10

12

HEB 140 (a) 2 1.8 1.6

Normalized Torque

1.4 1.2 1

N = 0.50

0.8

___ Displacement method (2 elements)

0.6

_._ Mixed method (2 elements) o Experimental data (Aalberg)

0.4 0.2 0 0

2

4

6

8

Normalized rotation

(b)

Fig. 18. Axial force and torsion (NT) combination: (a) test setup; (b) normalized load versus displacement curve.

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2541

Fig. 19. Nonredundant sway frame in strong-axis bending.

L sinusoidal initial sweep with a maximum out-of-straightness of 1000 at the mid-length is assumed in the outof-plane direction. The Galambos and Ketter [62] residual stress pattern is assumed. For the three-dimensional analysis, column AB is braced at top and bottom in the out-of-plane direction. End twist rotations are assumed to be fixed, but the out-of-plane rotations and the warping displacements at the column ends are modeled as unrestrained. Ten displacement-based finite elements or two mixed based finite elements are used to model column AB. Two- and three-dimensional strength interaction diagrams between the normalized axial load PPy and HL the lateral load 2M , generated by multiple analyses at different ratios of HP , is shown in Fig. 19. The strength p predicted by a combined set of nonlinear in-plane and out-of-plane beam-column design interaction curves from AS4100 [63] is also shown in the figure. From the design interaction curves, from the analysis-based strength curves, and from the load–displacement behavior in the different analysis solutions, it can be observed that the failure mode of member AB is predominantly in-plane for axial loads up to approximately the level associated with the out-of-plane strength of the member as a concentrically-loaded column. However, as P approaches the out-of-plane column strength, the strength predicted by the three-dimensional analysis deviates significantly from the analysis-based in-plane strength. It can be observed that the

2542

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

0.25

H Hp

0.2

0.15

___ 2D analysis, imperfect 0.1

_ _ 3D analysis, imperfect 0.05

0 0

0.001

0.002

0.003

0.004

0.005

0.006

0.007

0.008

0.009

0.01

∆ Lc

Fig. 20. Load versus deflection curves for

P Py

¼ 0:67.

three-dimensional analysis-based strength interaction curve matches well with the combined AS4100 design curves. This is significant since the AS4100 interaction formula tends to be more liberal than comparable curves in other steel design standards for cases similar to the one studied here (i.e., wide-flange beamcolumns with bdf  1 and subjected to reversed-curvature bending). An example load–deflection plot, corresponding to the application of PPy ¼ 0:67 followed by monotonic lateral loading through the limit load of the structure with PPy held constant, is shown in Fig. 20. The HL normalized horizontal load that can be supported by column AB (HHp ¼ 2M ) is sensitive to changes in PPy at p this axial load level. Nevertheless, the ductility of the beam-column in the three-dimensional failure mode is still good. That is, the post-peak load shedding by the column is rather mild with increasing lateral displacements, similar to that obtained from the comparable two-dimensional analysis.

4.8. Portal space frame The space frame shown in Fig. 21 is considered to demonstrate the applicability of the proposed elements to three-dimensional second-order inelastic analysis [64]. All the members in the frame are assumed to be initially straight. The Galambos and Ketter residual stress pattern [62] is assumed in all the members. As noted previously, the warping degree of freedom is assumed to be continuous if the finite elements meeting at a node are co-linear. On the other hand, the warping degree of freedom is assumed to be discontinuous at joints where the members intersect at an angle. Two mixed elements are used to model each member of the frame. The resulting load–deflection curve is shown in Fig. 21.

5. Conclusions An HR(u; R) mixed finite element formulation of a beam element that can represent the three-dimensional inelastic stability behavior of members of open-walled cross-section is presented. The beam kine-

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545 2.375 H

2543

2.875 H

L = 12 ft. 1.625 H

L = 12 ft.

2.125 H

18WF60 10WF60

H

1 in. = 25.4 mm 2

1 ksi = 6.9116 N/mm 1 kip = 4.459 kN E = 29000 ksi

Est = 290 ksi

σ y = 36 ksi

90 80

Lateral Load (kips)

70 60 50 ___ Mixed Method (2 elements per member)

40

_o_ Morino second-order analysis _x_ Ovunc second-order analysis

30 20 10 0 0

1

2

3

4

5

6

7

8

9

Horizontal displacement (in)

Fig. 21. Square space frame: load versus deflection curve.

matics include finite rotation and warping of the cross-section due to torsion. Material inelasticity is modeled using a two-space von-Mises model including the effect of shear stresses due to uniform torsion and the normal stresses due to axial force, biaxial bending, and bimoment. Numerical examples representing a wide range of two- and three-dimensional elastic and inelastic stability behavior are presented.

References [1] J.H. Argyris, H. Balmer, J.St. Doltsinis, P.C. Dunne, M. Haase, M. Kleiber, G.A. Malejannakis, H.P. Mlejnek, D.W. Scharpf, Finite element method––the natural approach, Comput. Methods Appl. Mech. Engrg. 17/18 (1979) 1–106. [2] J.H. Argyris, O. Hilpert, G.A. Malejannakis, D.W. Scharpf, On the geometric stiffness of a beam in space––a consistent v.w. approach, Comput. Methods Appl. Mech. Engrg. 20 (1979) 105–131. [3] K.J. Bathe, S. Bolourchi, Large displacement analysis of three-dimensional beam structures, Int. J. Numer. Methods Engrg. 14 (1979) 961–986. [4] J.C. Simo, A finite strain beam formulation. The three-dimensional dynamic problem. Part I, Comput. Methods Appl. Mech. Engrg. 49 (1985) 55–70.

2544

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

[5] J.C. Simo, L. Vu-Quoc, A three-dimensional finite-strain rod model. Part II: computational aspects, Comput. Methods Appl. Mech. Engrg. 58 (1986) 79–116. [6] A. Cardona, M. Geradin, A beam finite element non-linear theory with finite rotations, Int. J. Numer. Methods Engrg. 26 (1988) 2403–2438. [7] M.A. Crisfield, A consistent co-rotational formulation for nonlinear, three-dimensional, beam elements, Comput. Methods Appl. Mech. Engrg. 81 (1990) 131. [8] J.C. Simo, L. Vu-Quoc, A geometrically-exact rod model incorporating shear, torsion-warping deformation, Int. J. Solids Struct. 27 (3) (1991) 371–393. [9] A.F. Salleb, A.S. Gendy, Shear-flexible models for spatial buckling of thin-walled curved beams, Int. J. Numer. Methods Engrg. 31 (1991) 729–757. [10] Y.L. Pi, N.S. Trahair, Nonlinear inelastic analysis of steel beam-columns––theory, Research Report, School of Civil and Mining Engineering, The University of Sydney, NSW, Australia, November 1992. [11] A. Ibrahimbegovic, F. Frey, I. Kozar, Computational aspects of vector-like parametrization of three-dimensional finite rotations, Int. J. Numer. Methods Engrg. 38 (1995) 3653–3673. [12] A. Ibrahimbegovic, On finite element implementation of geometrically nonlinear Reissner’s beam theory: three-dimensional curved beam elements, Comput. Methods Appl. Mech. Engrg. 122 (1995) 11–26. [13] G. Jelenic, M. Saje, A kinematically exact space finite strain beam model - finite element formulation by generalized virtual work principle, Comput. Methods Appl. Mech. Engrg. 120 (1995) 131–161. [14] Y.B. Yang, W. McGuire, Stiffness matrix for geometric nonlinear analysis, J. Struct. Engrg.-ASCE 112 (4) (1986) 853–877. [15] Y.B. Yang, W. McGuire, Joint rotation and geometric nonlinear analysis, J. Struct. Engrg.-ASCE 112 (4) (1986) 878–905. [16] J.G. Orbison, W. McGuire, J.F. Abel, Yield surface applications in nonlinear steel frame analysis, Comput. Methods Appl. Mech. Engrg. 33 (1–3) (1982) 557–573. [17] Y. Zhao, Modelling of inelastic cyclic behavior of members, connections, joint panels of steel frames, Ph.D. thesis, Cornell University, 1993. [18] M. Attalla, G. Deierlein, W. McGuire, Spread of plasticity: quasi-plastic-hinge approach, J. Struct. Engrg. 120 (8) (1994) 2451– 2473. [19] M.S. Park, B.C. Lee, Geometrically nonlinear and elastoplastic three-dimensional shear flexible beam element of von-Mises-type hardening material, Int. J. Numer. Methods Engrg. 39 (1996) 383–408. [20] F. Taucer, E. Spacone, F.C. Filippou, A fiber beam-column element for seismic response analysis of reinforced concrete structures, Ucb/eerc-91/17, University of California, Berkeley, December 1991. [21] J.-M. Battini, C. Pacoste, Co-rotational beam elements with warping effects in instability problems, Comput. Methods Appl. Mech. Engrg. 191 (2002) 1755–1789. [22] J.-M. Battini, C. Pacoste, Plastic instability of beam structures using co-rotational elements, Comput. Methods Appl. Mech. Engrg. 191 (2002) 5811–5831. [23] B.A. Izzuddin, D.L. Smith, Large-displacement analysis of elastoplastic thin-walled frames. I. Formulation, implementation, J. Struct. Engrg.-ASCE 122 (8) (1996) 905–914. [24] B.A. Izzuddin, D.L. Smith, Large-displacement analysis of elastoplastic thin-walled frames. II. Verification and application, J. Struct. Engrg.-ASCE 122 (8) (1996) 915–925. [25] L.H. Teh, M.J. Clarke, Plastic-zone analysis of 3d steel frames using beam elements, J. Struct. Engrg.-ASCE 125 (11) (1999) 1328– 1337. [26] P.K.V.V. Nukala, D.W. White, Variationally consistent state determination algorithms for nonlinear analysis using mixed finite element formulations, Comput. Methods Appl. Mech. Engrg., in press. [27] V.Z. Vlasov, Thin walled elastic bars, Israel Program for Scientific Translation, Jerusalem, 1961. [28] S.P. Timoshenko, J.M. Gere, Theory of Elastic Stability, McGraw-Hill, New York, 1961. [29] A. Gjelsvik, The Theory of Thin Walled Bars, John Wiley and Sons, New York, 1981. [30] A. Billinghurst, J.R.L. Williams, G. Chen, N.S. Trahair, Inelastic uniform torsion of steel members, J. Struct. Mech. 8 (January) (1991) 61–92. [31] G. Chen, N.S. Trahair, Inelastic nonuniform torsion of steel i-beams, Research Report, r647, School of Civil and Mining Engineering, The University of Sydney, NSW, Australia, December 1991. [32] K.W. Spring, Euler parameters and the use of quaternion algebra in the manipulation of finite rotations: a review, Mech. Mach. Theory 21 (1986) 365–373. [33] Y.B. Yang, S.R. Kuo, Y.D. Cheng, Curved beam element for nonlinear analysis, J. Struct. Engrg.-ASCE 115 (1989) 840–855. [34] E.N. Dvorkin, D. Celentano, A. Cuitino, G. Gioai, A vlasov beam element, Comput. Struct. 33 (1) (1989) 187–196. [35] F. Gruttmann, R. Sauer, W. Wagner, A geometrical nonlinear eccentric 3d-beam element with arbitrary cross-sections, Comput. Methods Appl. Mech. Engrg. 160 (1998) 383–400. [36] M. Li, The finite deformation theory for beam, plate and shell. Part I: the two-dimensional beam theory, Comput. Methods Appl. Mech. Engrg. 146 (1997) 53–63.

P.K.V.V. Nukala, D.W. White / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2507–2545

2545

[37] M. Li, The finite deformation theory for beam, plate and shell structures. Part II: the kinematic model and the Green–Lagrange strains, Comput. Methods Appl. Mech. Engrg. 156 (1998) 247–257. [38] M. Li, The finite deformation theory for beam, plate and shell structures. Part III: the three-dimensional beam theory and Fe formulation, Comput. Methods Appl. Mech. Engrg. 162 (1998) 287–300. [39] L.E. Malvern, Introduction to the Mechanics of a Continuous Medium, Prentice-Hall, 1969. [40] W.F. Chen, T. Atsuta, in: Theory of Beam-Columns, vol. 2, McGraw Hill, 1977. [41] A.S. Gendy, A.F. Saleeb, T.Y.P. Chang, Generalized thin-walled beam models for flexural-torsional analysis, Comput. Struct. 42 (4) (1992) 531–550. [42] S.S. Antman, The theory of rods, in: Handbuch der Physik, vol. VIa/2, Springer, Berlin, 1972. [43] T.J.R. Hughes, The Finite Element Method, Prentice-Hall, 1987. [44] J.C. Simo, R.L. Taylor, A return mapping algorithm for plane stress elastoplasticity, Int. J. Numer. Methods Engrg. 22 (1986) 649–670. [45] H. Stolarski, T. Belytschko, Limitation principles for mixed finite elements based on the Hu-Washizu variational formulation, Comput. Methods Appl. Mech. Engrg. 60 (1987) 195–216. [46] D.N. Arnold, Mixed finite element methods for elliptic problems, Comput. Methods Appl. Mech. Engrg. 82 (1990) 281–300. [47] G. Alfano, F. Marotti de Sciarra, Mixed finite element formulations and related limitation principles: a general treatment, Comput. Methods Appl. Mech. Engrg. 138 (1996) 105–130. [48] J.C. Simo, M.S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Int. J. Numer. Methods Engrg. 29 (1990) 1595–1638. [49] A. Aalberg, An experimental study of beam-columns subjected to combined torsion, bending and axial actions, Ph.D. thesis, The Norwegian Institute of Technology, July 1995. [50] P.K.V.V. Nukala, Three-dimensional second-order inelastic analysis of steel frames, Ph.D. thesis, Purdue University, West Lafayette, IN 47907, May 1997. [51] M.J. Clarke, G.J. Hancock, A study of incremental-iterative strategies for non-linear analyses, Int. J. Numer. Methods Engrg. 29 (1990) 1365–1391. [52] M.A. Crisfield, Non-linear Finite Element Analysis for Solids and Structures, John Wiley and Sons, 1991. [53] J.H. Argyris, Sp. Symeonidis, Nonlinear finite element analysis of elastic systems under nonconservative loading-natural formulation. Part I. Quasistatic problems, Comput. Methods Appl. Mech. Engrg. 26 (1981) 75–123. [54] J.H. Argyris, Sp. Symeonidis, A sequel to: nonlinear finite element analysis of elastic systems under nonconservative loadingnatural formulation. Part I. Quasistatic problems, Comput. Methods Appl. Mech. Engrg. 26 (1981) 377–383. [55] A.K. Noor, J.M. Peters, Mixed models, selective/reduced integration displacement models for nonlinear analysis of curved beams, Int. J. Numer. Methods Engrg. 17 (1981) 615–631. [56] D.A. DaDeppo, R. Schmidt, Instability of clamped-hinged circular arches subjected to a point load, Trans. ASME 27 (1975) 894– 896. [57] R.D. Wood, O.C. Zienkiewicz, Geometrically nonlinear finite element analysis of beams, frames, arches and axisymmetric shells, Comput. Struct. 7 (1977) 725–735. [58] M. Borri, C. Bottasso, An intrinsic beam model based on a helicoidal approximation––Part I: formulation, Int. J. Numer. Methods Engrg. 37 (1994) 2267–2290. [59] M. Borri, C. Bottasso, An intrinsic beam model based on a helicoidal approximation––Part II: linearization and finite element implementation, Int. J. Numer. Methods Engrg. 37 (1994) 2291–2310. [60] C.G. Franchi, F. Montelaghi, A weak–weak formulation for large displacements beam statics: a finite volumes approximation, Int. J. Numer. Methods Engrg. 39 (1996) 584–604. [61] A. Avello, J. Garcia de Jalon, E. Bayo, Dynamics of flexible multibody systems using Cartesian co-ordinates and large displacement theory, Int. J. Numer. Methods Engrg. 32 (1991) 1543–1563. [62] T. Galambos, R. Ketter, Columns under combined bending and thrust, J. Engrg. Mech.-ASCE 85 (EM2) (1959) 1–30. [63] Australian Institute of Steel Construction, Standards Australia, AS4100-1990 Steel Structures, 1990. [64] S. Morino, Analysis of space frames, Ph.D. thesis, Lehigh University, 1970.