Inelastic nonuniform torsion of bars of doubly symmetric cross section by BEM

Computers and Structures 89 (2011) 2388–2401

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Inelastic nonuniform torsion of bars of doubly symmetric cross section by BEM E.J. Sapountzakis ⇑, V.J. Tsipiras School of Civil Engineering, National Technical University, Zografou Campus, GR-157 80 Athens, Greece

a r t i c l e

i n f o

Article history: Received 15 October 2010 Accepted 27 June 2011 Available online 20 July 2011 Keywords: Nonuniform torsion Inelastic torsion Warping shear stresses Distributed plasticity Boundary element method Bar

a b s t r a c t In this paper a boundary element method is developed for the inelastic nonuniform torsional problem of simply or multiply connected cylindrical bars of arbitrarily shaped doubly symmetric cross section taking into account the effect of warping shear stresses. The bar is subjected to arbitrarily distributed or concentrated torsional loading along its length, while its edges are subjected to the most general torsional boundary conditions. A displacement based formulation is developed and inelastic redistribution is modeled through a distributed plasticity model exploiting three dimensional material constitutive laws and numerical integration over the cross sections. An incremental – iterative solution strategy is adopted to restore global equilibrium along with an efficient iterative process to integrate the inelastic rate equations. Three boundary value problems with respect to the variable along the bar axis angle of twist, to the primary and to the secondary warping functions are formulated and solved employing the boundary element method. Numerical results are worked out to illustrate the method, demonstrate its efficiency and wherever possible its accuracy. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction In engineering practice we often come across the analysis of members of structures subjected to twisting moments. Curved bridges, ribbed plates subjected to eccentric loading or columns laid out irregularly in the interior of a plate due to functional requirements are most common examples. Moreover, design of bars and bar assemblages based on elastic analysis are most likely to be extremely conservative not only due to significant difference between initial yield and full plastification in a cross section, but also due to the unaccounted for yet significant reserves of strength that are not mobilized in redundant members until after inelastic redistribution takes place. Thus, material nonlinearity is important for investigating the ultimate strength of a bar that resists torsional loading, while distributed plasticity models are acknowledged in the literature [1–3] to capture more rigorously material nonlinearities than cross sectional stress resultant approaches [4] or lumped plasticity idealizations [5,6]. When an elastic bar is subjected to uniform torque arising from two concentrated torsional moments at its ends while the warping of the cross section is not restrained, the angle of twist per unit length remains constant along the bar. Under these conditions, the bar is leaded to uniform torsion and the well known primary (St. Venant) shear stress distribution arises [7]. When arbitrary torsional boundary conditions are applied either at the edges or ⇑ Corresponding author. E-mail addresses: [email protected] (E.J. Sapountzakis), tsipiras@gmail. com (V.J. Tsipiras). 0045-7949/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2011.06.009

at any other interior point of a bar due to construction requirements, this bar under the action of general twisting loading is leaded to nonuniform torsion and additional normal and secondary (warping) shear stresses arise [8], complicating the problem at hand. If inelastic effects are considered, especially through distributed plasticity formulations, then both the uniform [9,10] and the nonuniform [11,12] torsion problems require a much more rigorous analysis. Apart from research efforts in which bars are idealized with computationally demanding three dimensional [13] or shell [14] elements, several researchers proposed specialized beam elements to analyze bars under inelastic nonuniform torsion. Bathe and Wiener [14] employed Hermitian and isoparametric beam elements to model thin-walled I-beams using one element for each flange and one element for the web of the cross section. Wunderlich et al. [15] employed a power series numerical technique using an Updated Lagrangian formulation to study thin-walled beams under general loading conditions. Izzuddin and Lloyd Smith [16,17] followed an Eulerian approach employing cubic shape functions to study thin walled bars under general loading conditions. Chen and Trahair [18] and Pi and Trahair [12], using a mitre model [19] to describe the shear strain distribution over the cross section, developed finite element models for the inelastic analysis of thin-walled I-beams under torsion. Nie and Zhong [20] and Wang et al. [21] studied thin-walled beams under general loading conditions by employing an Updated Lagrangian description and the FE method. Alsafadie et al. [22] employed a corotational technique and a mixed FE formulation to analyze thin-walled cross section bars under general loading conditions.

E.J. Sapountzakis, V.J. Tsipiras / Computers and Structures 89 (2011) 2388–2401

The aforementioned contributions are applicable to cross sections of special geometry (e.g. thin-walled ones). A more general approach is presented by Baba and Kajita [11] who used a two-node, four-degree-of-freedom element for the longitudinal modeling and a four-node, 12-degree-of-freedom rectangular element for the cross sectional modeling of bars of bisymmetrical section under a total Lagrangian formulation, taking into account plasticity effects in the plane of the cross section as well. Moreover, in the more recent contributions of Battini and Pacoste [23], Nukala and White [1] and Gruttmann et al. [24] beams of arbitrarily shaped cross section under general loading conditions are analyzed. Battini and Pacoste [23], Nukala and White [1] use an independent warping parameter associated with the intensity of warping deformations to model nonuniform torsion. This parameter is eventually equated with the angle of twist per unit length in their analyses, thus only St. Venant shear stresses are taken into account. In all of the aforementioned research efforts warping shear stresses are not taken into account, with the exceptions of Wunderlich et al. [15], Gruttmann et al. [24], Nie and Zhong [20], Wang et al. [21] and Alsafadie et al. [22] who exploit a stress distribution arising from the introduction of an independent warping parameter in the displacement field of the bar. However, since this theory is analogous to the Timoshenko beam theory of shear – bending loading conditions, it does not satisfy local equilibrium equations under inelastic or even elastic conditions (for relevant discussions see for example Simo et al. [25] and Minghini et al. [26]). Moreover, the boundary element method (BEM) seems to be an alternative powerful tool for the solution of 3-D inelastic problems [27]. Beam problems are common cases of these in which the beam is modeled as 3-D inelastic solid and solved after discretizing its lateral surface with 2-D elements [28]. However, to the authors’ knowledge the BEM has not yet been employed for the numerical analysis of a 1d beam formulation of the aforementioned nonuniform torsional problem. In this paper a boundary element method is developed for the inelastic nonuniform torsional problem of simply or multiply connected cylindrical bars of arbitrarily shaped doubly symmetric cross section taking into account the effect of warping shear stresses. The bar is subjected to arbitrarily distributed or concentrated torsional loading along its length, while its edges are subjected to the most general torsional boundary conditions. A displacement based formulation is developed and inelastic redistribution is modeled through a distributed plasticity model exploiting three dimensional material constitutive laws and numerical integration over the cross sections. An incremental – iterative solution strategy is adopted to restore global equilibrium along with an efficient iterative process to integrate the inelastic rate equations [29]. Three boundary value problems with respect to the variable along the bar axis angle of twist, to the primary and to the secondary warping functions are formulated and solved employing the boundary element method [29]. The essential features and novel aspects of the present formulation compared with previous ones are summarized as follows. i. For the first time in the literature, warping shear stresses are taken into account in the problem at hand and their influence is quantified. Evaluation of both St. Venant and warping shear stresses is based on the solution of boundary value problems formulated by exploiting local equilibrium considerations under elastic conditions. ii. The formulation is a displacement based one taking into account inelastic redistribution along the bar axis by exploiting three dimensional material constitutive laws and numerical integration over the cross sections (distributed plasticity approach).

2389

iii. The cross section is an arbitrarily shaped thin- or thick walled doubly symmetric one. The formulation does not stand on the assumptions of a thin walled structure and therefore the cross section’s torsional and warping rigidities are evaluated ‘‘exactly’’ in a numerical sense. iv. An incremental – iterative solution strategy is adopted to restore global equilibrium of the bar. Integration of the inelastic rate equations is performed for each monitoring station with an efficient iterative process and stress resultants are obtained employing incremental strains. v. The developed procedure retains most of the advantages of a BEM solution over a pure domain discretization method, although it requires domain discretization to the longitudinal problem, having the following features.  Contrary to classical displacement based FE formulations, shape functions are not needed to approximate the angle of twist along the bar.  Cross sectional discretization is employed exclusively for numerical integration of domain integrals.  Finite differences of plastic terms are required exclusively if the effect of warping shear stresses is taken into account. 2. Statement of the problem 2.1. Displacements, strains, stresses Consider a prismatic bar of length l (Fig. 1) with an arbitrarily shaped doubly symmetric constant cross section, occupying the two dimensional multiply connected region X of the y, z plane bounded by the Cj(j = 1, 2, . . . , K) boundary curves, which are piecewise smooth, i.e. they may have a finite number of corners. In Fig. 1a Syz is an arbitrary coordinate system through the cross section’s shear center. The normal stress–strain relationship for the material is assumed to be elastic–plastic-strain hardening with initial modulus of elasticity, shear modulus and yield stress E, G and rY0, respectively. The bar is subjected to the combined action

(a)

(b)

Fig. 1. Prismatic bar of an arbitrarily shaped doubly symmetric cross section occupying region X (a) subjected to torsional loading (b).

2390

E.J. Sapountzakis, V.J. Tsipiras / Computers and Structures 89 (2011) 2388–2401

of arbitrarily distributed or concentrated twisting mt = mt(x) and warping mw = mw(x) moments acting in the x-direction (Fig. 1b). Under the aforementioned loading, the displacement field of the bar taking into account warping shear stresses is assumed to be given as S uðx; y; zÞ ¼ h0x ðxÞ/PS ðy; zÞ þ h000 x ðxÞ/S ðy; zÞ;

v ðx; zÞ ¼ zhx ðxÞ

ð1aÞ

wðx; yÞ ¼ yhx ðxÞ;

ð1b; cÞ

where u, v, w are the axial and transverse bar displacement components with respect to the Syz system of axes; h0x ðxÞ denotes the rate of change of the angle of twist hx(x) regarded as the torsional curvature, while h000 x ðxÞ is the third order derivative of the angle of twist with respect to x; /PS , /SS are the primary and secondary warping functions, respectively with respect to the shear center S [8]. Substituting Eq. (1) in the well known three-dimensional linear strain–displacement relations, the non vanishing (total) strain resultants are obtained as

exx ¼ h00x /PS ; cxy ¼

h0x

ð2aÞ

! @/PS @/SS  z þ h000 x @y @y

cxz ¼

h0x

! @/PS @/SS þ y þ h000 ; x @z @z ð2b; cÞ

where in Eq. (2a) the higher order contribution of the secondary warping function to the normal strain component is neglected [8]. The first terms of Eqs. (2b), (2c) are related to the primary shear stress distribution accounting for uniform torsion, while the last ones are related to the warping shear stress distribution accounting for nonuniform torsion. Considering strains to be small, employing the Cauchy stress tensor and assuming an isotropic and homogeneous material without exhibiting any damage during its plastification, the stress rates are defined in terms of the strain ones as [31,32]

8 9 2 38 el 9 E 0 0 > > < drxx > = < dexx > = 6 7 dcel ; dsxy ¼ 4 0 G 0 5 xy > > > : ; : el > ; 0 0 G dsxz dcxz

ð3Þ

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi     r2xx þ 3 s2xy þ s2xz  rY epleq ¼ 0;

ð4Þ

where rY is the yield stress of the material and epl eq is the equivalent plastic strain [33]. Moreover, the plastic modulus h is defined as h ¼ drY =depl eq . As long as the material remains elastic or elastic unloading occurs, the total strain components are equal to the elastic strain ones, thus the stress rates are expressed as (Eq. (3))

8 9 2 9 38 E 0 0 > > < drxx > = < dexx > = 6 7 dsxy ¼ 4 0 G 0 5 dcxy : > > > > : ; : ; dcxz 0 0 G dsxz |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl}

ð5Þ

½Del 

If plastic flow occurs then the total strain components are equal to sum of the elastic and plastic strain ones. Following standard procedures [33,34], the plastic strain components may be expressed with respect to the total ones, thus the stress rates–total strain rates relations are resolved as

ð6Þ

½Delpl 

with

  h  i c11 ¼ E hr2e þ 9G s2xy þ s2xz ; c ¼ hr2e þ Er2xx þ 9G s2xy þ s2xz ð7a; bÞ c21 ¼ 3EGrxx sxy

c31 ¼ 3EGrxx sxz ;

  c22 ¼ G hr2e þ Er2xx þ 9Gs2xz c32 ¼ 9G2 sxy sxz ; h i c33 ¼ G hr2e þ Er2xx þ 9Gs2xy

re ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   r2xx þ 3 s2xy þ s2xz :

ð7c; dÞ ð7e; fÞ ð7g; hÞ

By setting h = 0 (perfectly plastic case) in the above relations, the constitutive matrix presented by Baba and Kajita [11] is obtained, while if one of the shear stress components (along with the corresponding strain one) is dropped out, the constitutive relations presented by Chen and Trahair [18] are also precisely recovered. 2.2. Primary /PS and secondary /SS warping functions The primary and secondary warping functions are resolved by exploiting local equilibrium considerations along the longitudinal x-axis (together with the associated boundary condition at the lateral surface of the bar) under elastic conditions. Thus, /PS , /SS are evaluated by solving the boundary value problems [35,8]

r2 /PS ¼ 0 in X E G

@/PS ¼ zny  ynz @n

r2 /SS ¼  /PS ðy; zÞ in X

where d() denotes infinitesimal incremental quantities over time (rates) and the superscript el denotes the elastic part of the strain components. A Von Mises yield criterion, an associated flow rule and an isotropic hardening rule for the material are considered. The yield condition is described with the expression

f ¼

8 9 9 2 38 sym: > > < drxx > = 1 c11 < dexx > = 6 7 dsxy ¼ 4 c21 c22 5 dcxy > > > > c : ; : ; c31 c32 c33 dcxz dsxz |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

on Cj ;

@/SS ¼ 0 on Cj ; @n

ð8a; bÞ

ð9a; bÞ

where »2 = o2/ oy2 + o2/oz2 is the Laplace operator and o/on denotes the directional derivative normal to the boundary C. Since the problems at hand have Neumann type boundary conditions, each of the evaluated warping functions contains an integration constant (parallel displacement of the cross section along the bar axis). The evaluation of these constants is presented in [8]. Apparently, local equilibrium is violated when Eqs. (8) and (9) are introduced in Eq. (6) to perform stress calculations and inelastic redistribution is modeled only through the angle of twist hx(x). Eqs. (8) and (9) exhibit less accuracy as plastic deformations throughout the cross section increase. Introducing appropriate plastic terms in Eq. (8) to satisfy local equilibrium has already been achieved by Wagner and Gruttmann [9] and Sapountzakis and Tsipiras [34] for the uniform torsion problem under small and large displacements, respectively and by Baba and Kajita [11] for the geometrically nonlinear nonuniform torsion problem of bars of bisymmetrical cross section. However, none of these research efforts accounts for warping shear stresses, while it is pointed out in the literature [23] that such efforts lead to a much higher computational cost. An interesting and simple solution to the same problem is presented by Billinghurst et al. [19] and by Chen and Trahair [18] who employ a ‘‘mitre model’’ to distribute the primary shear stresses along the cross section in such a manner so as to resemble the shear stress distribution of bars under uniform torsion just before plastic collapse. However, this approach is limited in thin-walled bisymmetrical I-shaped cross sections, while it is expected to be less accurate in the elastic part of the bar and during the early phases of plastic loading.

2391

E.J. Sapountzakis, V.J. Tsipiras / Computers and Structures 89 (2011) 2388–2401

2.3. Equations of global equilibrium To establish global equilibrium equations, the principle of virtual work neglecting body forces is employed

Z V

ðrxx dexx þ sxy dcxy þ sxz dcxz Þ dV ¼

Z

ðt x du þ ty dv þ tz dwÞ dF;

F

where d() denotes virtual quantities, tx, ty, tz are the components of the traction vector, V is the volume and F is the surface area of the bar including the end cross sections. In the elastic case, shear stresses in the left hand side of Eqs. 000 0 (10) yield four sets of terms depending on h0x dh0x , h0x dh000 x , hx dhx and 000 000 hx dhx , respectively. The first term is related to primary shear stresses, the second and third terms are proved to vanish by exploiting Eqs. (8) and (9) and the homogeneity of the constitutive matrix, while the last term, related to warping shear stresses, yields a sixth order derivative of hx in the (Euler–Lagrange) governing equation of the bar. By dropping this term, the well known fourth order governing differential equation of the bar under nonuniform torsion is obtained [35]. Since throughout the present work, governing differential equations of up to the fourth order are considered, in the elastic range warping shear stresses do not affect global equilibrium. However, they are still evaluated through Eqs. (2)–(4), following the procedure presented in Sapountzakis and Mokos [8]. In the inelastic case, the incremental version of Eqs. (10) neglecting body forces

Z h

ðrxx þ Drxx ÞdDexx þ ðsxy þ Dsxy ÞdDcxy þ ðsxz þ Dsxz ÞdDcxz V Z dV ¼ ðt x dDu þ ty dDv þ tz dDwÞ dF

i ð11Þ

F

is employed, where D() denotes incremental quantities (over time). The incremental strains are obtained employing Eqs. (2). As in the elastic case, a sixth order governing differential equation is formulated due to the third order derivative h000 x appearing in Eq. (2b, c). By dropping terms of fifth order or higher, a fourth order governing differential equation is once again formulated. However, contrary to the elastic case, the inhomogeneity of the constitutive Eq. (6) (dependence on the longitudinal coordinate x) and the property of the fully populated constitutive matrix lead to nonvanishing terms of fourth or lower order related to warping shear stresses. Since all these terms depend on /SS , they are neglected in the formulation in order to reduce the computational cost of obtaining the tangent stiffness matrix of the bar. It is pointed out that this process affects only the convergence rate of the algorithm and not the equilibrium path. Moreover, the virtual terms dDcxy, dDcxz, dDu depend on the third order derivative dDh000 x which will eventually lead to the requirement of an additional (third) boundary condition apart from the well-known torsional and warping ones. In order to avoid the formulation of a third-order shear deformation theory, terms of dDh000 x are neglected in Eq. (11), resulting in an approximate solution of the problem at hand. Nevertheless, it is pointed out that warping shear stresses still affect the global equilibrium of the bar (as it will be subsequently presented in detail), without being completely neglected due to the aforementioned simplifications. The stress resultants of the bar are defined as

SM t ¼

X

!

sxy

ð13a; bÞ

X

ð10Þ

Z "

! !# Z " @/PS @/PS DSMt ¼ Dsxy  z þ Dsxz þ y dX @y @z X Z DSMw ¼ Drxx /PS dX;

@/PS @/PS  z þ sxz þy @y @z

!# dX SM w ¼

Z X

By exploiting Eqs. (1), (2), (12) and (13) and the aforementioned simplifications, the governing equation of the bar is obtained after some algebra through Eq. (11) as stiffness matrix  incremental displacements

zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ 2 dDSM t d DSMw  þ 2 dx dx

out-of-balance forces

zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ 2 dmw ðxÞ dSMt d SM w ¼ mt ðxÞ þ þ ;  2 dx ffl} dx dx ffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl externally applied loading

where SMt and SMw correspond to the internal twisting and warping moments of the bar, respectively and the relevant incremental stress resultants are defined as

internal stress resultants

along with its corresponding boundary conditions



a1 DSMt 

 dDSM w dSM w þ mw þ a2 Dhx ¼ a3  a1 SM t  ; dx dx ð15aÞ

b1 DSM w þ b2 Dh0x ¼ b3  b1 SMw ;

ð15bÞ

where ai, bi (i = 1, 2, 3) are functions specified at the bar ends. The boundary conditions (15) are the most general ones for the problem at hand, including also the elastic support. It is apparent that all types of the conventional boundary conditions (clamped, simply supported, free or guided edge) may be derived from Eq. (15) by specifying appropriately the functions ai and bi (e.g. for a clamped edge it is a2 = b2 = 1, a1 = a3 = b1 = b3 = 0). Finally, the expressions of the externally applied loading (mt, mw) with respect to the components of the traction vector (or the Cauchy stress tensor) may be easily resolved by virtue of Eq. (11) as

mt ðxÞ ¼

Z



 t y ðzÞ þ yt z ds mw ðxÞ ¼ 

C

Z

C

t x /PS ds:

ð16a; bÞ

Eqs. (14) and (15) must be brought into a more convenient form in order to employ the boundary element technique described in the next section to solve the problem numerically. The inhomogeneous elastoplastic constitutive matrix [Del–pl] of Eq. (6) can be always additively split into a homogeneous part (being the elastic constitutive matrix [Del] of Eq. (5)) and a remaining inhomogeneous one [Dpl]. Employing Eqs. (5), (6), and (22) the incremental stress resultants (Eq. (13)) are written as

DSMt ¼ GIt Dh0x þ DSMpl t

DSM w ¼ EC S Dh00x þ DSMpl w;

ð17a; bÞ

where It, CS are the torsion (St. Venant) and warping constants, respectively, depending exclusively on the geometry of the cross pl section and defined in Ref. [35], while DSMpl t , DSM w are expressed as 00 0 DSMpl t ¼ I1 Dhx þ I 2 Dhx

00 0 DSM pl w ¼ I 3 Dhx þ I1 Dhx ;

ð18a; bÞ

In the above equations, I1, I2, I3 are cross sectional properties depending on the geometry of the cross section and the progress of its plastic region given as

I1 ¼

rxx /PS dX; ð12a; bÞ

ð14Þ

I2 ¼

! !# Z " @/PS @/PS pl cpl  z þ c þ y /PS dX; 21 31 @y @z X Z X

2

!2 ! ! @/PS @/PS @/PS pl  z þ 2c32 z þy @y @y @z !2 3 @/PS þ y 5 dX; @z

ð19aÞ

4cpl 22

þ cpl 33

ð19bÞ

2392

I3 ¼

E.J. Sapountzakis, V.J. Tsipiras / Computers and Structures 89 (2011) 2388–2401

Z X

P 2 cpl dX; 11 /S

ð19cÞ

uðnÞ ¼

Z

cpl ij ðx; y; zÞ

where ¼ (i = 1, 2, 3, j = 1, 2, 3) are the coefficients of the [Dpl] matrix. Substituting Eqs. (17) in Eqs. (14) and (15), the following boundary value problem is obtained 00 EC S Dh0000 x  GI t Dhx 

ð24Þ where u⁄ is the fundamental solution given in [35]. Since ECS is independent of x, Eq. (24) can be written as

EC S uðnÞ ¼

2

dDSM pl d DSMpl t w þ 2 dx dx

dmw dSM t d SMw þ  2 dx dx dx

a1 GIt Dh0x  EC S Dh000x þ DSMplt 

inside the bar; dDSMpl w þ mw dx

þ a2 Dhx ð21aÞ

ð21bÞ

This boundary value problem by virtue of Eq. (18) can be rewritten as

I01 Dh00x 2

dmw dSMt d SM w þ  2 dx dx dx

4

dx

K4 ðrÞdx 3

2

d u

d u

inside the bar;



a1 GIt Dh0x  EC S Dh000x þ I2 Dh0x  I03 Dh00x  I3 Dh000x  I01 Dh0x þ mw  dSMw ; þ a2 Dhx ¼ a3  a1 SMt  dx

ð22Þ 

2

2

dDSM pl d DSM pl t w þ mt  2 dx 0 dx dx ! 2 dmw dSM t d SMw dx þ þ  2 dx dx dx " #l 3 2 d u d u du  EC S K4 ðrÞ 3  K3 ðrÞ 2 þ K2 ðrÞ  K1 ðrÞu : ð26Þ dx dx dx 0

EC S

du ¼ dn

Z 0

þ ð23aÞ þ

 b1 EC S Dh00x þ I3 Dh00x þ I1 Dh0x þ b2 Dh0x ¼ b3  b1 SMw at the bar ends x ¼ 0; l;

Z

l

K4 ðrÞ GIt

d u 2

þ

Having in mind that the lowest order derivative of u contained in Eq. (22) is the first one, Eq. (26) is differentiated once with respect to n and after carrying out several integrations by parts yields

0 00 0 00 0 00 00 000 EC S Dh0000 x  GI t Dhx  I 2 Dhx  I2 Dhx þ I 3 Dhx þ 2I3 Dhx

þ

"

d u

where r = x  n, x, n points of the bar, while the kernels Kj(r) (j = 1, 2, 3, 4) are also given in [35]. Solving Eq. (20) with respect to EC S Dh0000 x and substituting the result in Eq. (25), the following integral representation is obtained:

EC S uðnÞ ¼

at the bar ends x ¼ 0; l:

þ

EC S

ð25Þ !

  0 b1 EC S Dh00x þ DSM pl w þ b2 Dhx ¼ b3  b1 SM w

I001 Dh0x

4

l

#l du  EC S K4 ðrÞ 3  K3 ðrÞ 2 þ K2 ðrÞ  K1 ðrÞu ; dx dx dx 0

ð20Þ

 dSMw ¼ a3  a1 SMt  ; dx

I3 Dh0000 x

Z 0

2

¼ mt þ

¼ mt þ

d u

0

cpl ij

þ

" #l 3 2 @u d u @ 2 u du @ 3 u  d u  u dx  u  þ u ; 4 3 @x dx2 @x2 dx @x3 dx dx 0 4

l

ð23bÞ

where I01 , I02 , I03 , I001 , I003 are the derivatives with respect to x of the cross sectional properties of Eqs. (19). 3. Integral representations – numerical solution

þ

l

K2 ðrÞGIt Z Z Z

l 0

du dx þ dx

Z

l

0

K1 ðrÞDSMpl w dx 

l

K2 ðrÞmw dx þ 0

Z

K2 ðrÞDSM pl t dx

Z

l

K3 ðrÞmt dx

0 l

K2 ðrÞSM t dx

0 l

K1 ðrÞSMw dx 0

"

þ K3 ðrÞGIt

du dDSMpl w  K3 ðrÞDSMpl t þ K3 ðrÞ dx dx

 K2 ðrÞDSMpl w  K3 ðrÞSM t þ K3 ðrÞ

!#l du K3 ðrÞ 3  K2 ðrÞ 2 þ K1 ðrÞ dx dx dx 0 3

3.1. Integral representations for the angle of twist hx According to the precedent analysis, the inelastic nonuniform torsional problem of bars reduces to establishing the displacement component Dhx(x) having continuous derivatives up to the fourth order with respect to x and satisfying the boundary value problem described by the governing differential Eq. (22) along the bar and the boundary conditions (23) at the bar ends x = 0, l. The main difficulty of the problem at hand is the inhomogeneity of its governing equation. This boundary value problem (Eqs. (22) and (23)) is solved employing the BEM [30], as this is developed in [35,36] for the solution of a fourth order differential equation with constant coefficients, after modifying it as follows. The motivation to use this particular technique is justified from the intention to retain the advantages of a BEM solution over a domain approach, while keeping the finite differences required to solve the problem to a minimum and avoiding shape functions to approximate the displacement component Dhx(x). According to this method, let u(x) = Dhx(x) be the sought solution of the problem. The solution of the fourth order differential 4 4 equation d u=dx ¼ Dh0000 x is given in integral form as [35]

 K2 ðrÞSM w þEC S

dSM w dx

d u

2

d u

ð27Þ Eq. (27) may be used as a basis for the numerical solution of the problem at hand. However, it is observed that derivatives of plastic terms and stress resultants at the bar ends appear in the aforementioned equation. Finite differences required to resolve such terms may be avoided by expressing the boundary terms of Eq. (27) more conveniently as

du dDSM pl w  K3 ðrÞDSMpl  K2 ðrÞDSMpl t þ K3 ðrÞ w dx dx dSM w  K3 ðrÞSM t þ K3 ðrÞ  K2 ðrÞSM w dx ! 3 2 d u d u du þ EC S K3 ðrÞ 3  K2 ðrÞ 2 þ K1 ðrÞ dx dx dx

K3 ðrÞGIt

du tot   K3 ðrÞ DSM tot t þ SM t dx  K2 ðrÞðDSM w þ SMw Þ;

¼ EC S K1 ðrÞ

ð28Þ

2393

E.J. Sapountzakis, V.J. Tsipiras / Computers and Structures 89 (2011) 2388–2401 tot where DSMtot are the incremental and internal total twisting t , SM t moments at the bar ends, respectively, which by virtue of Eqs. (17) and (28) are defined as

DSM tot t ¼ DSM t 

dDSMw dx

SM tot t ¼ SM t 

dSM w : dx

ð29a; bÞ

Observing Eqs. (27) and (28), it is deduced that the integral representation (27) has been brought into a convenient form to establish a numerical computation of du/dx without any approximation of u(x) along the bar. Thus, the interval (0, l) is divided into L elements, on each of which du/dx is assumed to vary according to a certain law (constant, linear, parabolic, etc). The linear element assumption is employed here (Fig. 2) as the numerical implementation is simple and the obtained results are very good. This technique leads to differentiation of shape functions only for the approximation of d2u/dx2 appearing in Eq. (27) (through Eq. (18)). Employing the aforementioned procedure and a collocation technique, a set of L + 1 algebraic equations is obtained with respect to L + 7 unknowns, namely the values of (du/dx)i (i = 2, 3, . . . , L) at the L  1 internal nodal points and the boundary values of uj, (du/dx)j, ðDSMtot t Þj , (DSMw)j (j = 1, L + 1) at the bar ends n1 = 0, nL+1 = l (Fig. 2)). Two additional algebraic equations are obtained by applying the integral representation (26) at the bar ends n = 0, l after employing integration by parts, similar with those performed for Eq. (27). These L + 3 equations along with the four boundary conditions (Eq. (15)) yield a linear system of L + 7 simultaneous algebraic equations

½K T fDdg ¼ fbext g  fbint g;

ð30Þ

where [KT] is a known tangent generalized stiffness matrix, {Dd} is a generalized incremental unknown vector given as

     du du du du du fDdgT ¼  u1 uLþ1 dx 2 dx 3 dx L dx 1 dx Lþ1 o



 DSM tot DSMtot ðDSM w Þ1 ðDSMw ÞLþ1 ; ð31Þ t t 1 Lþ1 while {bext}, {bint} are known vectors representing all the terms related to the externally applied loading and the internal stress resultants, respectively. After solving the system of Eqs. (30), a post-processing step is required to obtain the derivatives d2u/dx2, d3u/dx3 at any nodal point ni (i = 1, 2, . . . , L + 1) which are essential to the developed solution algorithm. Differentiating Eq. (27) with respect to n, the following integral representations are obtained 2

EC S

d u 2

dn

Z

Z l du K1 ðrÞDSM pl dx  t dx dx 0 0 Z l Z l ðnÞ þ K ðrÞm dx  K1 ðrÞmw dx  DSM pl 2 t w

¼



K1 ðrÞGIt

0

Z

0

l

K1 ðrÞSMt dx  SM w ðnÞ

ð32aÞ 3

3

dn

Z l Z l du dx  K3 ðrÞDSMpl dx  K2 ðrÞDSMpl t w dx dx 0 0 0 Z l Z l Z l þ K4 ðrÞmt dx  K3 ðrÞmw dx  K3 ðrÞSM t dx 0 0 0

 Z l du K2 ðrÞSMw dx þ EC S K2 ðrÞ þ K1 ðrÞu  dx 0 l

tot  þ SM K ðrÞ ð D SM þ SM Þ ð33Þ þ K4 ðrÞ DSMtot þ 3 w w t t

EC S uðnÞ ¼ 

Z

l

K3 ðrÞGIt

0

The presented integral representations may also be employed to compute u(x) and its derivatives at any interior point of the bar other than ni (i = 2, 3, . . ., L). 3.2. Integral representations for the primary /PS and secondary /SS warping functions The evaluation of the primary warping function /PS (Eq. (8)) and of its derivatives with respect to y and z at any interior point for the calculation of stresses (Eqs. (2), (5) and (6)) or the axial displacement component (Eqs. (1a)) is accomplished using BEM as this is presented in [35]. The numerical solution of the boundary value problem described by Eq. (9) for the evaluation of the secondary warping function /SS (and of its derivatives with respect to y and z) is also accomplished using BEM as this is presented in [8,30]. 3.3. Incremental-iterative solution algorithm

 l tot þ K2 ðrÞðDSM tot t þ SM t Þ þ K1 ðrÞðDSM w þ SM w Þ 0 ;

d u

However it is pointed out that if warping shear stresses are ignored, d3u/dx3 does not need to be evaluated. Finally, the angle of twist u(x) = Dhx(x) may be computed by reformulating Eq. (26) as

l

0

EC S

Fig. 2. Discretization of the beam interval into linear elements, distribution of the nodal points and approximation of several quantities.

Z l du dDSMpl w þ DSM pl K1 ðrÞmt dx þ mw ðnÞ  t ðnÞ  dn dn 0

dSM w  tot l þ SM t ðnÞ  : ð32bÞ þ K1 ðrÞ DSM tot t þ SM t 0 dn

¼ GIt

Observing Eqs. (32b), it is concluded that the calculation of d3u/dx3 presumes the computation of derivatives related to plastic terms such as dSMw/dx, I01 , I03 (required to resolve dDSMpl w =dx). The evaluation of these derivatives results after approximating them with appropriate central, forward or backward finite differences [30]).

The modified or the fully nonlinear stiffness methods are usually employed for the incremental-iterative solution strategy of inelastic analysis of bars, that is the tangent stiffness matrix is kept constant throughout the whole increment or it is recomputed after each iteration, respectively. Apart from these methods, the initial stiffness method has also been implemented in the present study, since (i) it requires exclusively BEM computations to obtain the stiffness matrix (which is computed only once and stored in the beginning of the algorithm) and (ii) no convergence difficulties have arisen. Moreover, load control over the incremental steps is used and load stations are chosen according to the load history and convergence requirements. Usually, internal stress resultants and the externally applied loading are not in equilibrium, hence the out-of-balance forces do not vanish, thus an iterative process is initiated during an incremental step to restore equilibrium. Thus, using the subscript m to denote the load step, the superscript l to denote the iterative cycle and the symbol d() to denote iterative quantities, the lth iteration of the mth load step of the

2394

E.J. Sapountzakis, V.J. Tsipiras / Computers and Structures 89 (2011) 2388–2401

Fig. 3. Flowchart of the incremental-iterative solution algorithm. (See above mentioned references for further information.)

incremental-iterative solution algorithm of the modified stiffness method is outlined in the flowchart of Fig. 3 (full description of a similar algorithm is presented in [34]). 3.4. Computational aspects of the incremental - iterative solution algorithm (Fig. 3) The monitoring cross sections are located at the nodal points ni (i = 1, 2, . . ., L + 1) of the bar (Fig. 2). Internal stress resultants ½SMt ðni Þlm (i = 1, 2, . . ., L + 1), ½SM w ðni Þlm (i = 2, . . ., L) are evaluated

by employing a two-dimensional numerical integration scheme to approximate the domain integrals of Eqs. (12a, b) (step (vi)). In this study, the bar’s monitoring cross sections are divided into a number of triangular or quadrilateral cells and standard twodimensional Gauss quadrature rules are employed in each cell. Thus, the monitoring stations of each cross section coincide with the Gauss points of its cells. Exact patch between adjacent cells is not required due to the combined use of BEM to compute /PS , /SS and the fact that only domain integrals are approximated without performing any structural discretization. Moreover, the same

2395

E.J. Sapountzakis, V.J. Tsipiras / Computers and Structures 89 (2011) 2388–2401

tral and backward finite differences exploiting the values of SMw at the monitoring cross sections. Rl Rl  line integrals of the form 0 Ki ðrÞSMt dx, 0 Kj ðrÞSMw dx (i = 1, 2, 3, 4, j = 1, 2, 3) (eqns (27), (32a), (33)). A numerical integration scheme must be employed to resolve these integrals since stress resultants are not known in the whole bar interval (0, l). The most effective scheme that has been implemented is a semianalytical one, according to which SMt, SMw vary on an element k (k = 1, 2, . . ., L) (Fig. 2) of the bar interval following the same law that is used to approximate Dh0x (see Section 3.1). This leads to the integration of kernels being products of functions Ki(r) and two-node linear shape functions, thus it is performed analytically without any difficulty.

quadrature rules are employed to compute the cross sectional properties ½I1 ðni Þlm , ½I2 ðni Þlm , ½I3 ðni Þlm (i = 1, 2, . . ., L + 1) as well (step 3). The obtained stress resultants ½SM t ðni Þlm , ½SMw ðni Þlm (i = 1, 2, . . ., l L + 1), ½SMtot t ðnj Þm (j = 1, L + 1) from steps (vi), (vii) are employed to evaluate (a) the vector fbint glm representing all the terms of Eqs. (15), (27) and (33) related to stress resultants and (b) terms of Eqs. (32) related to stress resultants required to perform step (ii) for the next iteration l + 1 (step (viii)). Apart from elementary computations, this step requires the resolution of:  dSMw/dx at the bar ends and at the internal nodal points (Eq. (32b)). This is performed employing appropriate forward, cen-

Table 1 Twisting moment versus angle of twist at x = l/2 along with ultimate twisting moment M ut undertaken by the bar of Example 1 for various longitudinal discretization schemes. Ignoring warping shear stresses Number of elements Mt(l/2) (kipin) 8.851 77.887 88.507 101.784

14 hx ðl=2Þ 0.647 5.225 7.125 10.356 M ut ðkipinÞ 122.140

Considering warping shear stresses

30

42

14

30

42

0.645 5.844 7.310 10.771

0.645 5.887 7.345 10.837

0.647 5.356 7.915 14.424

0.645 5.861 7.439 12.290

0.645 5.890 7.598 12.142

122.140

123.025

105.324

111.519

113.290

125

Mt (kip*in)

100

Torque - Angle of twist at x=l/2 75

perfectly plastic - ignoring warping shear stresses perfectly plastic - with warping shear stresses strain hardening - ignoring warping shear stresses strain hardening - with warping shear stresses

50

25

0 0

5

10

15

20

25

30

Angle of twist (degrees)

Mt (Nm)

Fig. 4. Torque – angle of twist curves at x = l/2 for two cases of analysis and two cases of materials of the bar of Example 1.

34 32 30 28 26 24 22 20 18 16 14 12 10 8 6 4 2 0

Angle of twist at x=l present study (BEM) - without warping shear stresses present study (BEM) - with warping shear stresses

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

Angle of twist (rad) Fig. 5. Torque – angle of twist curves at x = l of the bar of Example 2.

0.4

0.45

2396

E.J. Sapountzakis, V.J. Tsipiras / Computers and Structures 89 (2011) 2388–2401

Moreover, the requirements of various computations of step 4 are similar to the ones described above. The integrals appearing in [KT]m+1 are also treated with a semi-analytical scheme, according to which I1, I2, I3 vary on an element k (k = 1, 2, . . . , L) (Fig. 2) of the bar interval following the same law that is used to approximate Dh0x (see Section 3.1). The resolution of terms of eqns (15), (27), (33) and (32a,b) related to externally applied loading (step 5) requires the computaRl Rl tion of line integrals of the form 0 Ki ðrÞmt dx, 0 Kj ðrÞmw dx (i = 1, 2, 3, 4, j = 1, 2, 3). Since the distributions of mt, mw are usually prescribed in codes and regulations with simple analytical relations, these integrals are evaluated analytically, demonstrating the efficiency of the developed numerical procedure (e.g. concentrated loads may be treated using the Dirac function, without adhering to any simplifications). It is finally worth noting that step (ii) of the solution algorithm should be enriched properly by applying Eq. (33), if values of hx(x) at any interior point of the bar are monitored.

4. Numerical examples On the basis of the analytical and numerical procedures presented in the previous sections, a computer program has been written and representative examples have been studied to demonstrate wherever possible the accuracy, the efficiency and the range of applications of the developed method. Example 1. In the first example, for comparison purposes an open thin-walled I-shaped (total height h = 0.1524 m, total width b = 0.1509 m, flange width tf = 0.0122 m, web width tw = 0.0080 m) cross section bar (E = 213,400 MN/m2, G = 80,000 MN/ m2, rY0 = 285 MN/m2) of length l = 1.93 m subjected to a monotonically increasing concentrated twisting moment at its midpoint and restrained against twisting at both ends (free warping boundary conditions) has been studied, employing 500 boundary elements, 64 quadrilateral cells and a 3  3 Gauss integration scheme for each cell (cross sectional discretization). Two material

600

Mt (kN*m)

500

400

Angle of twist at x=2.4m Without warping shear stresses With warping shear stresses

300

200

100

0

0.000

0.004

0.008

0.012

0.016

0.020

Angle of twist (rad) Fig. 6. Torque – angle of twist at x = 2.4 m curves for two cases of analysis of the first load case of Example 3.

Table 2 Spread of plasticity at selected cross sections of the bar of Example 2 subjected to M ext ¼ 31:0 N m. t x=0 With warping shear stresses

Without warping shear stresses

x = l/3

x = 2l/3

x=l

E.J. Sapountzakis, V.J. Tsipiras / Computers and Structures 89 (2011) 2388–2401

2397

700 600

Mt (kN*m)

500

Angle of twist at x=2.8m

400

Present method - Without warping shear stresses Present method - With warping shear stresses 3-D FEM solution [39]

300 200 100 0 0.000

0.005

0.010

0.015

0.020

0.025

0.030

Angle of twist (rad) Fig. 7. Torque – angle of twist at x = 2.8 m curves for two cases of analysis and a 3-D FEM solution [39] of the second load case of Example 3.

Number of plastified Gauss points

(a) 600 550 500 450

Mt = 220 kNm Without warping shear stresses With warping shear stresses

400 350 300 250 200 150 100 50 0 2.30 2.35 2.40 2.45 2.50 2.55 2.60 2.65 2.70 2.75 2.80 2.85 2.90 2.95 3.00

(b)

650

Number of plastified Gauss points

x(m)

600 550 500

Mt = 300 kNm Without warping shear stresses With warping shear stresses

450 400 350 300 250 200 150 100 50 0 2.30 2.35 2.40 2.45 2.50 2.55 2.60 2.65 2.70 2.75 2.80 2.85 2.90 2.95 3.00

x(m) Fig. 8. Number of plastified Gauss points at selected cross sections along the bar of Example 3 subjected to the first load case with torque levels (a) Mext ¼ 220:0 kN m and (b) t M ext ¼ 300:0 kN m. t

2398

E.J. Sapountzakis, V.J. Tsipiras / Computers and Structures 89 (2011) 2388–2401

cases have been analyzed, namely an elastic-perfectly plastic (Et = 0) one and an elastoplastic-strain hardening (Et = 6000 MN/ m2) one. To demonstrate convergence characteristics of the developed numerical procedure (elastic-perfectly plastic case), in Table 1 pairs of applied twisting moment and angle of twist values at the midpoint of the bar are presented taking into account or ignoring warping shear stresses, for three longitudinal discretization schemes. The first torque level (Mt(l/2) = 8.851 kipin) of the table corresponds to an elastic behavior, while the remaining ones refer to inelastic response. Moreover, in this table the ultimate twisting moment M ut that can be undertaken by the bar (plastic collapse load) is also presented for the aforementioned longitudinal discretization schemes. From this table, it is concluded that the scheme with only 14 elements yields acceptable results when warping shear stresses are ignored (especially in predicting M ut value), while finer discretizations should be employed if these stresses are considered in the analysis. In Fig. 4, the torque – angle of twist curves (42 linear longitudinal elements) at the midpoint of the bar are presented for two cases of analysis (considering or ignoring warping shear stresses) and two cases of materials (elastic–perfectly plastic, elastoplastic–strain hardening) as compared with those obtained from FEM solutions for thin-walled sections [17,18], ignoring warping shear stresses and experimental results [38]. From this figure a very good agreement between the results of the proposed method and those of the literature is observed. From Table 1 and Fig. 4 it is also concluded that warping shear stresses increase the torsional flexibility

Example 2. In the second example, also for comparison purposes an open thin-walled I-shaped cross section (total height h = 0.02794 m, total width b = 0.0254 m, flange and web width tf = tw = 0.00254 m) cantilever of length l = 0.254 m subjected to a monotonically increasing concentrated twisting moment at its right end and restrained against twisting and warping at its left end has been studied, employing 30 linear longitudinal elements, 500 boundary elements, 60 quadrilateral cells and a 3  3 Gauss integration scheme for each cell (cross sectional discretization). The material is considered to be elastic–perfectly plastic with E = 206,850 MN/m2, G = 79,558 MN/m2, rY0 = 206.850 MN/m2. In Fig. 5, the torque – angle of twist curves at the bar’s right end are presented considering or ignoring warping shear stresses, as compared with those obtained from a FEM solution for thin-walled sections [18] ignoring warping shear stresses. A very good agreement is once again verified between the corresponding curve of the present study and the one of Chen and Trahair [18], especially in the inelastic regime, as well as with the one of Bathe and Wiener [14]. Moreover, it is concluded that warping shear stresses slightly decrease the torsional rigidity and the ultimate torsional loading that the bar can undertake. In order to demonstrate the influence of warping shear stresses on kinematical components, the variation of the angle of twist along the bar axis is depicted in Fig. 6 for several inelastic torque levels, considering or ignoring warping shear

650

Number of plastified Gauss points

(a)

and decrease the ultimate torsional loading that the bar can undertake.

Mt = 220 kNm Without warping shear stresses With warping shear stresses

600 550 500 450 400 350 300 250 200 150 100 50

0 2.700 2.725 2.750 2.775 2.800 2.825 2.850 2.875 2.900 2.925 2.950 2.975 3.000

x(m)

Mt = 300 kNm

650

Number of plastified Gauss points

(b)

Without warping shear stresses With warping shear stresses

600 550 500 450 400 350 300 250 200 150 100 50

0 2.700 2.725 2.750 2.775 2.800 2.825 2.850 2.875 2.900 2.925 2.950 2.975 3.000

x(m) Fig. 9. Number of plastified Gauss points at selected cross sections along the bar of Example 3 subjected to the second load case with torque levels (a) M ext ¼ 220:0 kN m and t (b) M ext ¼ 300:0 kN m. t

2399

E.J. Sapountzakis, V.J. Tsipiras / Computers and Structures 89 (2011) 2388–2401

stresses. The conclusions already drawn from Fig. 5 and Example 1 are once again verified. Finally, in Table 2 the spread of plasticity at selected cross sections is presented for an inelastic torque level (Mext ¼ 31:0 N m), considering or ignoring warping shear stresses. t From this table both the influence of the aforementioned stresses especially at the clamped edge and the good agreement with the corresponding plasticity patterns obtained from a FEM solution [14] ignoring warping shear stresses are once more verified.

Example 3. In order to demonstrate the range of applications of the developed method, in the third example a thick-walled rectangular cross section (h = 0.60 m, b = 0.30 m) bar of length l = 3.0 m subjected to two load cases of a monotonically increasing concentrated twisting moment near its right end and restrained against twisting at both ends (free warping) has been studied, employing 40 linear longitudinal elements, 400 boundary elements, 72 quadrilateral cells and a 3  3 Gauss integration scheme

Table 3 Cross sectional distribution of Von Mises stresses rVM (103 MPa) at x = 2.55 m considering (a, b) or ignoring (c, d) warping shear stresses of the bar of Example 3 subjected to two inelastic torque levels of the first load case. M ext ¼ 350:0 kN m t

Mext ¼ 400:0 kN m t

0.02 0.019 0.018 0.017 0.016 0.015 0.014 0.013 0.012 0.011 0.01 0.009 0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001

0.02 0.019 0.018 0.017 0.016 0.015 0.014 0.013 0.012 0.011 0.01 0.009 0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001

(a)

0.02 0.019 0.018 0.017 0.016 0.015 0.014 0.013 0.012 0.011 0.01 0.009 0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001

(b)

0.02 0.019 0.018 0.017 0.016 0.015 0.014 0.013 0.012 0.011 0.01 0.009 0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001

(c)

(d)

2400

E.J. Sapountzakis, V.J. Tsipiras / Computers and Structures 89 (2011) 2388–2401

for each cell (cross sectional discretization). The material is considered to be elastic–perfectly plastic with E = 32318.4 MN/ m2, G = 13,466 MN/m2, rY0 = 20 MN/m2. In the first load case the twisting moment is eccentrically placed at x = 2.4 m, while in Fig. 7 the torque – angle of twist curves at the same point are presented considering or ignoring warping shear stresses verifying the conclusions already drawn from the previous examples concerning the effect of warping shear stresses. To demonstrate the amplified role in some cases of these stresses, in the second load case the concentrated twisting moment is placed more closely at the bar’s right end (x = 2.8 m). In Fig. 8 the torque – angle of twist curves at the same point are presented considering or ignoring warping shear stresses, indicating that the effects of these stresses are more pronounced in this load case rather than the first one. Moreover, in this figure, the torque – angle of twist curve (along with the ultimate torsional loading) obtained from a 3-D FEM solution [39] employing 4320 solid (brick) elements is also presented for comparison purposes. A good agreement is observed between this solution and the analysis considering warping shear stresses, especially in predicting the torque – angle of twist curve before plastic collapse. In Figs. 9 and 10 the number of plastified Gauss points at selected (near the bar’s right end) cross sections situated at the nodal points along the bar axis are presented for two inelastic torque levels (Mext ¼ 220:0 kN m, M ext ¼ 300:0 kN m) of the first and second t t load cases, respectively. From these figures, it is easily concluded that warping shear stresses affect inelastic redistribution patterns of bars under nonuniform torsion. In this case, it is observed that if warping shear stresses are taken into account, inelastic redistribution is initiated at positions where external twisting moments are placed, in which normal and warping shear stresses get their maximum values (Figs. 9 and 10a). At higher torque levels, plasticity spreads towards the bar’s right end where primary shear stresses mainly occur (Figs. 9 and 10b). If warping shear stresses are ignored, inelastic redistribution is initiated at the bar’s right end, where primary shear stresses get their maximum values. At higher torque levels, plasticity spreads towards positions where external twisting moments are placed. To demonstrate the spread of plasticity throughout the cross section, in Table 3 the distribution of Von Mises stresses rVM ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2xx þ 3ðs2xy þ s2xz Þ at the cross section with x = 2.55 m are presented for two inelastic torque levels of the first load case, namely M ext ¼ 350:0 kN m and M ext ¼ 400:0 kN m, considering or ignoring t t warping shear stresses, respectively.

5. Concluding remarks In this paper, a boundary element method is developed for the inelastic nonuniform torsional problem of simply or multiply connected cylindrical bars of arbitrarily shaped doubly symmetric cross section taking into account the effect of warping shear stresses. The main conclusions that can be drawn from this investigation are a. The numerical technique presented in this investigation is well suited for computer aided analysis of cylindrical bars of arbitrarily shaped doubly symmetric cross section, supported by the most general torsional boundary conditions and subjected to the action of arbitrarily distributed or concentrated torsional loading. The present formulation may be extended to more general cross section geometry and loading conditions by exploiting techniques such as the ones employed in [40–42].

b. Both St. Venant and warping shear stresses are evaluated from the solution of two boundary value problems formulated under elastic conditions. c. Accurate results are obtained using a relatively small number of nodal points across the longitudinal axis. d. Warping shear stresses slightly decrease the torsional rigidity and the ultimate torsional loading that bars can undertake under nonuniform torsion. e. Warping shear stresses have an influence on inelastic redistribution patterns of bars under nonuniform torsion. f. The developed procedure retains most of the advantages of a BEM solution over a domain approach.

Acknowledgements The authors would like to thank the Senator Committee of Basic Research of the National Technical University of Athens, Programme ‘‘PEVE-2008’’, R.C. No: 65 for the financial support of this work. References [1] Nukala P, White D. A mixed finite element for three-dimensional nonlinear analysis of steel frames. Comput Methods Appl Mech Eng 2004;193:2507–45. [2] Teh L, Clarke M. Plastic-zone analysis of 3D steel frames using beam elements. J Struct Eng 1999;125:1328–37. [3] Saritas A, Filippou FC. Frame element for metallic shear-yielding members under cyclic loading. J Struct Eng 2009;135:1115–23. [4] Attalla MR, Deierlein GG, McGuire W. Spread of plasticity: quasi-plastic-hinge approach. J Struct Eng 1994;120:2451–73. [5] Orbison JG, McGuire W, Abel JF. Yield surface applications in nonlinear steel frame analysis. Comput Methods Appl Mech Eng 1982;33:557–73. [6] Ngo-Huu C, Kim S, Oh J. Nonlinear analysis of space steel frames using fiber plastic hinge concept. Eng Struct 2007;29:649–57. [7] Timoshenko SP, Goodier JN. Theory of elasticity. New York, NY: McGraw-Hill; 1970. [8] Sapountzakis EJ, Mokos VG. Warping shear stresses in nonuniform torsion by BEM. Comput Mech 2003;30:131–42. [9] Wagner W, Gruttmann F. Finite element analysis of Saint–Venant torsion problem with exact integration of the elastic–plastic constitutive equations. Comput Methods Appl Mech Eng 2001;190:3831–48. [10] Nadai A. Theory of flow and fracture of solids. New York, NY: McGraw-Hill; 1954. [11] Baba S, Kajita T. Plastic analysis of torsion of a prismatic beam. Int J Numer Methods Eng 1982;18:927–44. [12] Pi YL, Trahair NS. Inelastic torsion of steel I-beams. J Struct Eng, ASCE 1995;121:609–20. [13] May IM, Al-Shaarbaf AS. Elasto-plastic analysis of torsion using threedimensional finite element model. Comput Struct 1989;33:67–678. [14] Bathe K-J, Wiener PM. On elastic–plastic analysis of I-beams in bending and torsion. Comput Struct 1983;17:711–8. [15] Wunderlich W, Obrecht H, Schroedter V. Nonlinear analysis and elastic–plastic load-carrying behaviour of thin-walled spatial beam structures with warping constraints. Int J Numer Methods Eng 1986;22:671–95. [16] Izzuddin BA, Lloyd Smith D. Large-displacement analysis of elastoplastic thinwalled frames. I: Formulation and implementation. J Struct Eng 1996;122:905–13. [17] Izzuddin BA, LloydSmith D. Large-displacement analysis of elastoplastic thinwalled frames. II: Verification and application. J Struct Eng 1996;122:915–25. [18] Chen G, Trahair NS. Inelastic nonuniform torsion of steel I-beams. J Constructional Steel Res 1992;23:189–207. [19] Bilinghurst A, Williams JRL, Chen G, Trahair NS. Inelastic uniform torsion of steel members. Comput Struct 1992;42:887–94. [20] Nie G, Zhong Z. The elasto-plastic and geometrically nonlinear finite element model of space beam considering restraint torsion. Key Eng Mater 2007;340– 341 I:335–40. [21] Wang X, Yang Q, Zhang Q. A new beam element for analyzing geometrical and physical nonlinearity. Acta Mech Sinica/Lixue Xuebao 2010;26:605–15. [22] Alsafadie R, Hjiaj M, Battini J-M. Three-dimensional formulation of a mixed corotational thin-walled beam element incorporating shear and warping deformation. Thin-Walled Struct 2011;49:523–33. [23] Battini J-M, Pacoste C. Plastic instability of beam structures using co-rotational elements. Comput Methods Appl Mech Eng 2002;191:5811–31. [24] Gruttmann F, Sauer R, Wagner W. Theory and numerics of three-dimensional beams with elastoplastic material behaviour. Int J Numer Methods Eng 2000;48:1675–702. [25] Simo JC, Hjelmstad KD, Taylor RL. Numerical formulations of elastoviscoplastic response of beams accounting for the effect of shear. Comput Methods Appl Mech Eng 1984;42:301–30.

E.J. Sapountzakis, V.J. Tsipiras / Computers and Structures 89 (2011) 2388–2401 [26] Minghini F, Tullini N, Laudiero F. Locking-free finite elements for shear deformable orthotropic thin-walled beams. Int J Numer Methods Eng 2007;72:808–34. [27] Banerjee PK. The boundary element methods in engineering. London, United Kingdom: McGraw-Hill; 1994. [28] Beskos D, Maier G, editors. Boundary element advances in solid mechanics. Udine, Italy: Springer-Verlag; 2003. [29] Ortiz M, Simo J. Analysis of a new class of integration algorithms for elastoplastic constitutive relations. Int J Numer Methods Eng 1986;23:353–66. [30] Katsikadelis JT. Boundary elements: theory and applications. Elsevier: Amsterdam-London, United Kingdom; 2002. [31] Vlasov V. Thin-walled elastic beams. Israel Program for Scientific Translations, Jerusalem; 1963. [32] Armenakas AE. Advanced mechanics of materials and applied elasticity. New York: Taylor & Francis Group; 2006. [33] Crisfield MA. Non-linear finite element analysis of solids and structures: essentials. vol. 1, New York, USA: John Wiley and Sons. [34] Sapountzakis EJ, Tsipiras VJ. Nonlinear inelastic uniform torsion of bars by BEM. Comput Mech 2008;42:77–94.

2401

[35] Sapountzakis EJ. Solution of non-uniform torsion of bars by an integral equation method. Comput Struct 2000;77:659–67. [36] Sapountzakis EJ. Nonuniform torsion of multi-material composite bars by the boundary element method. Comput Struct 2001;79:2805–16. [37] Isaacson E, Keller H. Analysis of numerical methods. New York, USA: Wiley; 1966. [38] Farwell Jr CR, Galambos TV. Nonuniform torsion of steel beams in inelastic range. J Struct Eng ASCE 1969;95:2813–29. [39] NX Nastran User’s Guide, Siemens PLM Software Inc.; 2007. [40] Sapountzakis EJ, Dourakopoulos JA. Flexural–Torsional nonlinear analysis of Timoshenko beam-column of arbitrary cross section by BEM. Comput Mater Continua 2010;18:121–53. [41] Dourakopoulos JA, Sapountzakis EJ. Postbuckling analysis of beams of arbitrary cross section using BEM. Eng Struct 2010;32:3713–24. [42] Sapountzakis EJ, Dikaros IC. Non-linear flexural–torsional dynamic analysis of beams of arbitrary cross section by BEM. Int J Non-Linear Mech 2011;46:782–94.