Non-linear elastic non-uniform torsion of bars of arbitrary cross-section by BEM

ARTICLE IN PRESS International Journal of Non-Linear Mechanics 45 (2010) 63–74

Contents lists available at ScienceDirect

International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm

Non-linear elastic non-uniform torsion of bars of arbitrary cross-section by BEM E.J. Sapountzakis , V.J. Tsipiras Institute of Structural Analysis, School of Civil Engineering, National Technical University of Athens, Zografou Campus, GR-157 80 Athens, Greece

a r t i c l e in fo

abstract

Article history: Received 7 April 2008 Received in revised form 16 August 2009 Accepted 14 September 2009

In this paper an elastic non-uniform torsion analysis of simply or multiply connected cylindrical bars of arbitrary cross-section taking into account the effect of geometric non-linearity is presented employing the boundary-element(BE) method. The torque–rotation relationship is computed based on the finitedisplacement (finite-rotation) theory, that is the transverse displacement components are expressed so as to be valid for large rotations and the longitudinal normal strain includes the second-order geometric non-linear term often described as the ‘‘Wagner strain’’. The proposed formulation does not stand on the assumption of a thin-walled structure and therefore the cross-section’s torsional rigidity is evaluated exactly without using the so-called Saint-Venant’s torsional constant. The torsional rigidity of the cross-section is evaluated directly employing the primary warping function of the cross-section depending on its shape. Three boundary-value problems with respect to the variable along the beam axis angle of twist, to the primary and to the secondary warping functions are formulated. The first one, employing the Analog Equation Method (a BEM-based method), yields a system of non-linear equations from which the angle of twist is computed by an iterative process. The remaining two problems are solved employing a pure BE method. Numerical results are presented to illustrate the method and demonstrate its efficiency and accuracy. The developed procedure retains most of the advantages of a BEM solution over a pure domain discretization method, although it requires domain discretization. & 2009 Elsevier Ltd. All rights reserved.

Keywords: Non-uniform torsion Non-linear analysis Shear stresses Warping Bar Beam Twist Boundary-element method Shear center Wagner strain

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. When arbitrary torsional boundary conditions are applied either at the edges or at any other interior point of the bar due to construction requirements, this bar under the action of general twisting loading is led to non-uniform torsion. Besides, since thinwalled open sections have low torsional stiffness, the torsional deformations can be of such magnitudes that it is not adequate to treat the angles of cross-section rotation as small. When finite twist rotation angles are considered, the elastic non-uniform torsion problem becomes non-linear. Moreover, this problem becomes much more complicated in the case the cross-section’s centroid does not coincide with its shear center (asymmetric beams), leading to the formulation of a flexural–torsional coupled

 Corresponding author. Tel.: + 302107721718; fax: + 302107721720.

E-mail addresses: [email protected] (E.J. Sapountzakis), [email protected] (V.J. Tsipiras). 0020-7462/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijnonlinmec.2009.09.003

problem. The extensive use of the aforementioned structural elements necessitates a reliable and accurate analysis of bars of arbitrary cross-sections subjected to torsional loading taking into account geometric non-linearity. Though several researchers have dealt either the linear nonuniform torsional behaviour of beams [1–4] or with the nonlinear uniform torsional behaviour of doubly symmetric beams [5,6], to the author’s knowledge very little work has been done on the corresponding non-linear non-uniform torsional problem of arbitrary cross-section beams. Tso and Ghobarah [7] and Attard in [8,9] have presented a set of displacement relationships for a straight prismatic thin-walled open beam applicable to situations where displacements are finite, cross-section do not distort, strains are small and flexural displacements are small to moderate while cross-sectional twist can be large. The presented numerical examples in these studies concern only uniform torsion of either mono- or doubly symmetric cross-sections. Finally, Trahair [10], employing the finite-element method and presenting examples of only doubly symmetric cross-sections, and Mohri et al.[11], employing similar equations to those established by Attard [8] and presenting examples of either doubly symmetric crosssections subjected to non-uniform torsion or buckling or postbuckling behaviour of arbitrary cross-section beams, also analyze the non-linear non-uniform torsional problem. Nevertheless, all of

ARTICLE IN PRESS 64

E.J. Sapountzakis, V.J. Tsipiras / International Journal of Non-Linear Mechanics 45 (2010) 63–74

the aforementioned studies, which are the only ones considering finite angles of twist in asymmetrical bars (and taking into account all of the arising non-linear terms), are not general since they are restricted to thin-walled beams. Moreover, to the authors’ knowledge the boundary-element(BE) method has not yet been used for numerical analysis of the aforementioned problem. In this paper an elastic non-uniform torsion analysis of simply or multiply connected cylindrical bars of arbitrary cross-section taking into account the effect of geometric non-linearity is presented employing the boundary-element method. The torque–rotation relationship is computed based on the finitedisplacement (finite-rotation) theory, that is the transverse displacement components are expressed so as to be valid for large rotations and the longitudinal normal strain includes the second-order geometric non-linear term often described as the ‘‘Wagner strain’’. The proposed formulation does not stand on the assumption of a thin-walled structure and therefore the crosssection’s torsional rigidity is evaluated exactly without using the so-called Saint-Venant’s torsional constant. The torsional rigidity of the cross-section is evaluated directly employing the primary warping function of the cross-section [3] depending on its shape. Three boundary-value problems with respect to the variable along the beam axis angle of twist, to the primary and to the secondary warping functions are formulated. The first one, employing the Analog Equation Method [12], a BEM-based method, yields a system of non-linear equations from which the angle of twist is computed by an iterative process. The rest two problems are solved employing a pure BEM-based method. The proposed formulation procedure is based on the assumption of no local or lateral torsional buckling or distortion and includes the following essential features and novel aspects compared with previous ones: i. large deflections and rotations are taken into account, that is the strain-displacement relationships contain higher-order displacement terms; ii. for the first time in the literature, the present formulation is applicable to arbitrarily shaped thin or thick-walled crosssections occupying simple or multiple connected domains; iii. for the first time in the literature, both the linear and the nonlinear parts of the secondary warping shear stresses arising from the non-uniform torsion are evaluated as well as all the linear and non-linear stress resultants, including the nonlinear warping moment; iv. the presented formulation does not stand on the assumption of a thin-walled structure and therefore the cross-section’s torsional rigidity is evaluated exactly without using the socalled Saint-Venant’s torsional constant and v. the developed procedure retains most of the advantages of a BEM solution over a pure domain discretization method, although it requires domain discretization. Numerical results are presented to illustrate the method and demonstrate its efficiency and accuracy.

2. Statement of the problem 2.1. Displacements, strains and stresses Let us consider a prismatic bar of length l (Fig. 1), of constant arbitrary cross-section of area A. The homogeneous isotropic and linearly elastic material of the bar’s cross-section, with modulus of elasticity E, shear modulus G and Poisson’s ratio n occupies the two-dimensional multiply connected region O of the y–z plane

Fig. 1. Prismatic bar subjected to arbitrarily distributed conservative twisting and warping moments (a) with a cross-section of arbitrary shape occupying the twodimensional region O (b).

and is bounded by the Gj(j= 1, 2 ,y, K) boundary curves, which are piecewise smooth, i.e. they may have a finite number of corners. In Fig. 1b CYZ is the principal coordinate system through the crosssection’s centroid C, while yC, zC are its coordinates with respect to Syz system of axes through the cross-section’s shear center S. The bar is subjected to the combined action of the arbitrarily distributed or concentrated conservative twisting mt =mt(x) and warping mw =mw(x) moments acting in the x direction (Fig. 1a). Under the aforementioned loading the displacement field of the bar with respect to the Syz system of axes for large twisting rotations and small bending ones is given as u ¼ um ðxÞ þ yY ðxÞðz  zC Þ P

S

yZ ðxÞðy  yC Þ þ y0 x ðxÞfS ðy; zÞ þ fS ðx; y; zÞ

ð1aÞ

v ¼ vS ðxÞ  zsin yx ðxÞ  yð1  cos yx ðxÞÞ

ð1bÞ

w ¼ wS ðxÞ þysin yx ðxÞ  zð1  cos yx ðxÞÞ

ð1cÞ

where the transverse displacement components v, w are valid for large rotations [13]; yY, yZ are the angles of rotation due to bending of the cross-section with respect to its centroid; y0 x(x) denotes the rate of change of the angle of twist yx regarded as the torsional curvature; vS(x), wS(x) are the transverse displacement components of the shear center S; fSP and fSS are the primary and secondary warping functions with respect to the shear center S, respectively [3]; and um(x) is an ‘‘average’’ axial displacement of the cross-section of the bar, which will be later explained. Substituting Eqs. (1a)–(1c) in the non-linear (Green) strain– displacement relations of the non-vanishing strains "   2  2 # @u 1 @u 2 @v @w ð2aÞ exx ¼ þ þ þ @x 2 @x @x @x

gxy ¼

@u @v @u @u @v @v @w @w þ þ þ þ @y @x @y @x @y @x @y @x

ð2bÞ

gxz ¼

@w @u @u @u @v @v @w @w þ þ þ þ @x @z @z @x @z @x @z @x

ð2cÞ

ARTICLE IN PRESS E.J. Sapountzakis, V.J. Tsipiras / International Journal of Non-Linear Mechanics 45 (2010) 63–74

assuming moderate large deflections ((qu/qx)2 5 qu/qx, (qu/qx)(qu/ qy)5 (qv/qx)+(qu/qy), (qu/qx) 5(qw/qx)+(qu/qz)) and having in mind that ignoring shear deformation effect the angles of rotation due to bending, yY, yZ, approximate the slope between the longitudinal axis in the deformed state and in the initial vertical and horizontal planes, respectively, the following relations are obtained:

yY ðxÞ ¼ v0 S sin yx  w0 S cos yx 0

ð3aÞ

0

yZ ðxÞ ¼ v S cos yx þ w S sin yx

1 y x ðyC yY þ zC yZ Þ þ ½ðv0 S Þ2 þ ðw0 S Þ2 þ ðy2 þ z2 Þðy0 x Þ2  2

gxy ¼ y x

! P S @fS @fS z þ @y @y

gxz ¼ y0 x

! P S @fS @fS þy þ @z @z

ð4aÞ

ð4bÞ

ð4cÞ

kY ðxÞ ¼ v00 S ðxÞsin yx  w00 S ðxÞcos yx

ð5aÞ

kZ ðxÞ ¼ v00 S ðxÞcos yx þ w00 S ðxÞsin yx

ð5bÞ

while the second-order geometrically non-linear term in the right hand side of Eq. (4a) (y 2 + z 2) (yx0 )2/2 is often described as the ‘‘Wagner strain’’ [10]. It is worth noting here that in obtaining Eq. (4a) the rate of change of the secondary warping function fSS, that is the normal stress arising due to the secondary shear one due to warping [14], has been ignored. Considering strains to be small, employing the second Piola–Kirchhoff stress tensor and assuming an isotropic and homogeneous material for zero Poisson ratio, the stress components are defined in terms of the strain ones as 8 9 2 9 38 e = E 0 0 > > < Sxx > = < xx > 7 Sxy ¼ 6 ð6Þ 4 0 G 0 5 gxy > > :S > ; :g > ; 0 0 G xz xz

Sxz ¼ Gy0 x

! P S @fS @f þy þG S @z @z

E G

@fS ¼0 @n

ð10Þ

E G

ð11Þ

on G

ð12Þ

where r2 =q2/qy2 + q2/qz2 is the Laplace operator and q/qn denotes the directional derivative normal to the boundary G. Thus, the primary fSP and the secondary fSS warping functions will be evaluated from the solution of the Neumann problems described by the governing equations (9) and (11) inside the two-dimensional multiply connected region O, subjected to the boundary conditions (10) and (12) on its boundary G, respectively. It is worth noting here that the evaluated warping functions due to the solution of the corresponding Neumann problems contain an integration constant (parallel displacement of the crosssection along the bar axis), which can be obtained from the conditions that Z fPS dO ¼ 0 ð13Þ O

Z O

fSS dO ¼ 0

ð14Þ

Moreover, in the case where the origin O of the coordinates is a point of the plane other than the shear center S, the primary warping function with respect to this point fOP is first established from the Neumann problem (9) and (10) substituting fSP by fOP. Using the evaluated warping function fOP, fSP is then established using the transformation given in [16], while the coordinates of the shear center S with respect to the Oy z system of coordinates are established from Eq. (13) and the conditions Z fPS y dO ¼ 0 ð15aÞ

O

ð7aÞ

ð7bÞ

ð7cÞ

2.2. Equations of local equilibrium Applying the stress components (7a)–(7c) in the first equation of equilibrium and neglecting the body forces [15]: @Sxx @Sxy @Sxz þ þ ¼0 @x @y @z

on G

and for the secondary warping function fSS as

Z

P

Sxx ¼ E½u0 m þ kY ðz  zC Þ  kZ ðy  yC Þ þ y00 x fS ðy; zÞ E y0 x ðyC yY þ zC yZ Þ þ ½ðv0 S Þ2 þ ðw0 S Þ2 þ ðy2 þ z2 Þðy0 x Þ2  2

Sxy ¼ Gy x

P

@fS ¼ zny  ynz @n

O

or, employing Eqs. (4a)–(4c) as

! P S @fS @f  z þG S @y @y

ð9Þ

S

where the curvature components kY, kz are given from the following relations:

0

= 2 fPS ¼ 0 in O

P = 2 fSS ¼  ym ðy2 þ z2 Þy0x y00x x fS ðy; zÞ 

exx ¼ u0 m þ kY ðz  zC Þ  kZ ðy  yC Þ þ y00 x fPS ðy; zÞ

0

parts of both the Eq. (8) and the traction vector on the free surface of the bar to vanish, we obtain the following Neumann problems for the primary warping function fSP as

ð3bÞ

and the non-vanishing strain resultants are given as

0

65

fPS z dO ¼ 0

ð15bÞ

that define point S as the point for which the axial and bending stress resultants arising from an angle of twist at this point vanish in a geometrically linear analysis. It is worth noting here that the shear center axis of an asymmetric cross-section bar subjected to pure torsion is not a rectilinear one but is assigned from the displacement field uS(x), vS(x), wS(x). Based on the conditions (13) and (14) the meaning of the ‘‘average’’ axial displacement of the cross-section of the bar can now be explained by the relation Z  Z u dO ¼ um ðxÞA þ yY ðxÞ z dO  zC A O O Z  Z Z y dO  yC A þ y0 x fPS dO þ fSS dO yZ ðxÞ O

O

O

ð16Þ ð8Þ

ignoring the independent of warping terms such as E u00 m or warping terms due to shear such as E(w000 S(z  zC)+ v000 S(y yC)) and requiring both the primary and the secondary due to warping

leads to R um ðxÞ ¼

O u dO

A

ð17Þ

ARTICLE IN PRESS 66

E.J. Sapountzakis, V.J. Tsipiras / International Journal of Non-Linear Mechanics 45 (2010) 63–74

that is um generally does not coincide either with the axial displacement of the cross-section’s centroid uC or with that of the shear center uS. 2.3. Equations of global equilibrium

Z

Neglecting body forces the principle of virtual work yields Z ðSxx dexx þ Sxy dgxy þ Sxz dgxz ÞdV ¼ ðtx du þ ty dv þ tz dwÞdA

V

Z

CS ¼

b2 ¼

1 2IYY

ð18Þ

O

Z

MY ¼

Sxx ðz  zC ÞdO

ð19bÞ

O

Z

MZ ¼ 

Sxx ðy  yC ÞdO

ð19cÞ

O

"

Z

MtP ¼

SPxy

O

Mw ¼ 

Z

! !# P P @fS @fS  z þSPxz þy dO @y @z P

O

Sxx fS dO

ð19dÞ

ð19eÞ

where MtP is the primary twisting moment [3] resulting from the primary shear stress distribution SxyP, SxzP and Mw is the warping moment due to the torsional curvature. Substituting Eqs. (7a)– (7c) into (19a)–(19e) the stress resultants are obtained as     1 IP ðv0 S Þ2 þðw0 S Þ2 þ ðy0 x Þ2  y0 x ðyC yY þ zC yZ Þ N ¼ EA u0 m þ A 2 ð20aÞ MY ¼ EIYY ½kY þ b2 ðy0 x Þ2 

ð20bÞ

MZ ¼ EIZZ ½kZ  b1 ðy0 x Þ2 

ð20cÞ

MtP

ð20dÞ

GIt yx0

¼



Mw ¼  ECS yx00 þ

Uw 0 2 ðyx Þ 2CS

r 2 dO

ð21bÞ

O

Z

IYY ¼

2

ðz  zC Þ dO

ð21cÞ

ðy  yC Þ2 dO

ð21dÞ

O

IZZ ¼

Z O

It ¼

Z

P

y2 þ z2 þ y O

Z O

P

!

@fS @f  z S dO @z @y

ð21eÞ

ð21fÞ

Z

r 2 ðz  zC ÞdO

ð22bÞ

O

fPS r2 dO

ð22cÞ

It is worth noting here that the aforementioned stress resultants refer to the directions of the cross-section at its deformed configuration (they take into account the cross-sections’ rotations), since they have been defined with respect to the second Piola–Kirchhoff stress tensor. Moreover, though the direction of the axial force N is tangential to the deformed centroid axis, it is not necessarily applied to the cross-section’s centroid, since as it is already mentioned the ‘‘average’’ axial displacement um is not necessarily the centroid’s displacement uC [8]. In deriving the relations (20a)–(20e) the properties of the principal system of axes and the orthogonality conditions of the primary warping function fSP with respect to the shear center S have been employed. Substituting the stress components given in Eqs. (7a)–(7c) and the strain resultants given in Eqs. (4a)–(4c) to the principle of virtual work (Eq. (18)) the equation of torque equilibrium of the bar is obtained as NyC yZ y0 x þ NzC yY y0 x  MZ kU þ MY kZ   d 1  MtP þ EIn ðy0 x Þ3 þ Cy0 x  NyC yY  NzC yZ dx 2 

d2 M w d ¼ mt ðxÞ þ ½mw ðxÞ dx dx2

ð23Þ

and the corresponding boundary conditions are written as   1 ðMw Þ0 þ MtP þ EIn ðy0 x Þ3 þ Cy0 x  NyC yY  NzC yZ  Mt dyx ¼ 0 2 ð24aÞ 

  Mw þMw dy0 x ¼ 0

ð24bÞ

where the ‘‘higher-order torsion constant’’ In and the variable C are given as

ð20eÞ

O

IP ¼

Uw ¼

In ¼ IPP 



where the area A, the polar moment of inertia IP with respect to the shear center S, the moments of inertia IYY, IZZ with respect to the cross-section’s centroid, the torsion constant It and the warping constant CS with respect to the shear center S are given as Z A¼ dO ð21aÞ Z

ðfS Þ2 dO

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi with r ¼ y2 þ z2 , while the geometric constants b1, b2 and Uw are given as Z 1 b1 ¼ r 2 ðy  yC ÞdO ð22aÞ 2IZZ O

A

where d(  ) denotes virtual quantities, tx, ty, tz are the components of the traction vector with respect to the undeformed surface of the bar including the cross-section ends, A is the surface area and V the volume of the bar at its initial undeformed state. Moreover, the stress resultants of the bar are given as Z Sxx dO ð19aÞ N¼

P

O

C¼N

IP2 U2 2 2  4b1 IZZ  4b2 IYY  w A CS

IP Uw  2MZ b1 þ 2MY b2  Mw A CS

ð25aÞ

ð25bÞ

while Mt and Mw in the boundary conditions (24a) and (24b) are the externally applied conservative twisting and warping moR ments at the bar ends, respectively, and IPP = Or4 dO (in Eq. (25a)). It is worth noting here that in deriving the equation of torque equilibrium of the bar the secondary shear stress distribution has been ignored [14]. Considering that the axial force and the bending moments vanish along the bar N¼0

ð26aÞ

MY ¼ 0

ð26bÞ

MZ ¼ 0

ð26cÞ

the governing equation of the non-uniform non-linear torsional problem is written as ECS yx0000 

3 d EIn2 ðyx0 Þ2 yx00  GIt yx00 ¼ mt ðxÞ þ ½mw ðxÞ 2 dx

ð27Þ

ARTICLE IN PRESS E.J. Sapountzakis, V.J. Tsipiras / International Journal of Non-Linear Mechanics 45 (2010) 63–74

subjected to the following boundary conditions at the bar ends:

a1 Mt þ a2 yx ¼ a3

ð28aÞ

b1 Mw þ b2 y0 x ¼ b3

ð28bÞ

where In2 is defined as In2 ¼ In þ

2 Uw CS

1 EIn2 ðyx0 Þ3 2   Uw 0 2 ðy x Þ Mw ¼  ECS y00 x þ 2CS

d4 yx ¼ qðxÞ dx4

ð33Þ

ð30aÞ

ð30bÞ

yx ðxÞ ¼

Mt and Mw are the twisting and warping moments at the boundary of the bar, respectively, given as m

According to this formulation, if yx is the sought solution of the boundary-value problem described by Eqs. (27), (28a) and (28b), differentiating this function four times yields

Eq. (33) indicates that the solution of the original problem can be obtained as the angle of twist of a bar with unit warping rigidity and vanishing torsional and higher-order rigidity subjected to a torsional fictitious load q(x) under the same boundary conditions. The fictitious load is unknown. However, it can be established using BEM as follows. The solution of Eq. (33) is given in integral form as

ð29Þ

Mt ¼  ECS yx þGIt yx0 þ

67

while ai, bi (i= 1, 2, 3) are functions specified at the boundary of the bar. The boundary conditions (28a) and (28b) are the most general linear torsional boundary conditions for the beam problem including also the elastic support. It is apparent that all types of conventional torsional boundary conditions (clamped, simply supported, free or guided edge) can be derived form these equations by specifying appropriately the functions ai and bi (e.g. for a clamped edge it is a2 = b2 =1, a1 = a3 = b1 = b3 = 0). The solution of the boundary-value problem described by Eqs. (27), (28a) and (28b) for the evaluation of the angle of twist yx assumes that the warping CS, the torsion It and the geometric Uw constants defined from Eqs. (21e), (21f) and (22c) are already established. Eqs. (21e), (21f) and (22c) indicate that the evaluation of the aforementioned constants presumes that the primary warping function fSP at any interior point of the domain O of the cross-section of the bar is evaluated after solving the boundaryvalue problem described by Eqs. (9) and (10). Once yx is established, the secondary warping function fSS at any interior point of the domain O of the cross-section of the bar is evaluated after solving the boundary-value problem described by Eqs. (11) and (12). Subsequently the second Piola–Kirchhoff stress components are evaluated employing Eqs. (7a)–(7c). Moreover, using the evaluated angle of twist yx, the transverse displacement components of the shear center vS(x), wS(x) can be established integrating the following relations

Z

 l    3 dyx d2 yx d2 yx dyx d3 yx   d yx  qyx dx  yx  þ y x dx dx2 dx3 dx2 dx dx3 0 0 l

ð34Þ

where y*x is the fundamental solution, which is given as

yx ¼

1 3 l ð2 þ jrj3  3jrj2 Þ 12

ð35Þ

with r = r/l and r = x  x, x, x points of the bar, which is a particular singular solution of the equation 

d4 yx ¼ dðx; xÞ dx4

ð36Þ

Employing Eq. (35) the integral representation (34) can be written as

yx ðxÞ ¼

Z

 l d3 yx d2 yx dyx qL4 ðrÞdx  L4 ðrÞ 3 þ L3 ðrÞ 2 þ L2 ðrÞ þ L1 ðrÞyx dx dx dx 0 0 l

ð37Þ

v00 S ¼  ðy0 x Þ2 ðb2 sin yx  b1 cos yx Þ

ð31aÞ

where the kernels Li(r), (i =1, 2, 3, 4) are given in [17]. Notice that in Eq. (37) for the line integral it is r =x  x, x, x points inside the bar, whereas for the remaining terms r =x  z, x inside the bar and z at the bar ends 0 and l. Differentiating Eq. (37) results in the integral representations of the derivatives of the angle of twist yx as given in [17]. The integral representation (37) and the one of y0 x written for bar ends 0, l together with the boundary conditions (28a) and (28b) can be employed to express the unknown boundary quantities yx, y0 x, y00 x and y0000 x in terms of q. This is accomplished numerically. If L is the number of nodal points along the bar axis, this procedure yields the following set of non-linear equations:

w00 S ¼ ðy0 x Þ2 ðb2 cos yx þ b1 sin yx Þ

ð31bÞ

½e11 fyx g þ ½e12 fyx0 g þ ½e13 fy 0 x g þ ½e14 fy x g ¼ fa3 g

obtained from the solution of the linear system of equations arising from the substitution of Eqs. (20b)–(26c) after employing Eqs. (5a) and (5b). Finally, using the evaluated derivatives of the transverse displacement components v0 S, w0 S the ‘‘average’’ axial displacement of the cross-section of the bar um can be established after substitution of Eqs. (20a)–(26a) as   1 IP u0 m ¼  ð32Þ ðv0 S Þ2 þ ðw0 S Þ2 þ ðy0 x Þ2 þ y0 x ðyC yY þzC yZ Þ A 2 and subsequent integration.

3. Integral representations—numerical solution 3.1. For the angle of twist yx The numerical solution of the non-linear boundary-value problem described by Eqs. (27), (28a) and (28b) will be accomplished by the Analog Equation Method [12]. This method is applied for the aforementioned problem using the formulation presented in [17].

3

m

2

½e21 fy0 x g þ ½e22 fy0 x g þ ½e23 fy00 x g ¼ fb3 g

ð38aÞ ð38bÞ

m

½e31 fyx g þ ½e32 fyx0 g þ ½e33 fyx00 g þ ½e34 fy x g ¼ ½f1 fqg m

½e42 fyx0 g þ ½e43 fyx00 g þ ½e44 fy x g ¼ ½f2 fqg

ð38cÞ ð38dÞ

in which [e1i], (i=1, 2, 3, 4), [e2j], (j=1, 2, 3) are 2  2 matrices including the nodal values of the functions a1, a2, b1, b2 of Eqs. (28a) and (28b), respectively, and [eij], (i=3, 4, j=1, 2, 3, 4) are square 2  2 known coefficient matrices resulting from the values of the kernels Li at the bar ends; {a3}, {b3} are 2  1 known column matrices including the boundary values of the functions a3, b3 in Eqs. (28a) and (28b) and [fi] (i=1, 2) are 2  L rectangular known matrices originating from the integration of the kernels on the axis of the bar. Finally, {yx}, {y0 x}, {y0 x2}, {y0 x3}, {y00 x} and {y000 x} are vectors including the two unknown nodal values of the respective boundary quantities and {q} is a vector including the L unknown nodal values of the fictitious load. Elimination of the linear boundary quantities yx, y00 x and y000 x from Eq. (38a)–(38d) yields 2

3

½e1 fy0 x g þ ½e2 fy0 x g þ ½e3 fy0 x g þfe4 g ¼ ½f fqg

ð39Þ

ARTICLE IN PRESS 68

E.J. Sapountzakis, V.J. Tsipiras / International Journal of Non-Linear Mechanics 45 (2010) 63–74

in which [ei], (i=1, 2, 3) are 2  2 known matrices; {e4} is a 2  1 known column matrix and [f] is a 2  L rectangular known matrix. The discretized counterparts of Eq. (37) and of the integral representations of the derivatives of yx when applied to all nodal points in the interior of the bar and after the elimination of the linear boundary quantities yx, y00 x and y000 x yield 2

3

fYx g ¼ ½Bfqg þ ½C1 fy0 x g þ ½C2 fy0 x g þ ½C3 fy0 x g 2

ð40aÞ 3

fY0 x g ¼ ½B0 fqg þ½C 0 1 fy0 x g þ ½C 0 2 fy0 x g þ ½C 0 3 fy0 x g 2

ð40bÞ 3

fY00 x g ¼ ½B00 fqg þ ½C 00 1 fy0 x g þ ½C 00 2 fy0 x g þ ½C 00 3 fy0 x g 2

3

fYx g ¼ ½Bfqg þ ½C1m fy0 x g þ½C2m fy0 x g þ ½C3m fy0 x g

ð40cÞ ð40dÞ

where [B], [B0 ], [B00 ], [B000 ] are known L  L coefficient matrices and [Ci], [C0 i], [C00 i], [C000 i], (i= 1, 2, 3) are L  2 also known matrices. The aforementioned elimination is performed using Eqs. (38a)–(38d) for homogeneous boundary conditions (28a,b) (a3 = b3 = 0). For non-homogeneous boundary conditions, an additive constant vector will appear in the right hand side of Eqs. (40a)–(40d). Finally, applying Eq. (27) to the L nodal points in the interior of the bar the following non-linear system of equations with respect to q and to the boundary quantities y0 x is obtained:  n o n o



2 3 ð½ECS   GIt ½B00 Þ q  GIt ½C 00 1  y0 x þ½C 00 2  y0 x þ½C 00 3  y0 x 3  EIn2 fB0 1 gfqgfB0 1 gfqgfB00 1 gfqg 2

...

fB0 2 gfqgfB0 2 gfqgfB00 2 gfqg

fB0 L gfqgfB0 L gfqgfB00 L gfqg

T

n

n 2o n 3o o 3  EIn2 Fð½B0 ; ½B00 ; ½C 0 1 ; ½C 0 2 ; ½C 0 3 ; ½C 00 1 ; ½C 00 2 ; ½C 00 3 ; y0 x ; y0 x ; y0 x ; q Þ 2 ¼ fmt g þ fm0 w g

ð41Þ

where [ECS] is a diagonal L  L matrix including the values of the warping rigidity at the L nodal points in the interior of the bar; {B0 i}, {B00 i} are L  1 column matrices including the values of the ith row of the corresponding matrices; {mt} and {m0 w} are vectors with L elements including the values of the fictitious loading, the torsional one and the derivative of the warping loading. Eqs. (41) and (39) constitute a non-linear system of equations with respect to q and y0 x. This system is solved using iterative numerical methods, such as the two-term acceleration method [18–19]. Following the solution of the non-linear system of equations and the evaluation of the fictitious load q and the boundary quantities y0 x, the rest of the boundary quantities are evaluated employing Eqs. (38a)–(38d), while the angle of twist and its derivatives in the interior of the bar are obtained using Eqs. (40a)–(40d). 3.2. For the primary warping function fSP The integral representations and the numerical solution for the evaluation of the angle of twist yx presented in the previous section assume that the warping CS, the torsion It and the geometric Uw constants from Eqs. (21e), (21f) and (22c) are already established. Eqs. (21e), (21f) and (22c) indicate that the evaluation of the aforementioned constants presumes that the primary warping function fSP and its derivatives with respect to y and z at any interior point of the domain O of the cross-section of the bar are known. Once fSP and its derivatives are established, CS and It constants are evaluated by converting the domain integrals into line integrals along the boundary using the corresponding relations presented in Sapountzakis and Mokos [4], as well as Uw, b1 and b2 using the following relations:  Z  1 3 1 4 P ð42aÞ ðy ny þz3 nz ÞfS  ðy þ z4 Þðzny  ynz Þ ds Uw ¼ 12 G 3

b1 ¼

1 2IZZ

Z  4  Z y yz3 ny þ nz ds  yC 3 G 4



ð42bÞ

G

  y3 z3 ny þ nz ds 3 3

b2 ¼

1 2IYY

Z  3  Z y z z4 ny þ nz ds  zC 4 3 G

 3   y z3 ny þ nz ds 3 G 3

ð42cÞ

along with constant boundary elements for the approximation of these line integrals. Moreover, the evaluation of the primary warping function fSP with respect to the shear center S and its derivatives with respect to y and z at any interior point for the calculation of stress resultants (Eqs. (7b) and (7c)) is accomplished using BEM as presented in Sapountzakis [16]. Finally, since the non-uniform non-linear torsional problem of composite bars is solved by the BEM, the domain integrals in Eqs. (21a)–(21d) have to be converted to boundary line integrals. This can be achieved using integration by parts, the Gauss theorem and the Green identity [14].

3.3. For the secondary warping function fSS Numerical solution of the boundary-value problem described by Eqs. (11) and (12) for the evaluation of the secondary warping function fSS will be accomplished using BEM [20]. This method is applied for the aforementioned problem using the formulation presented in [3]. According to this formulation, the Green identity  Z Z  @j @C j ds ð43Þ ðC r2 j  j r2 CÞdO ¼ C @n @n O G when applied to the secondary warping function fSS and to the fundamental solution

C¼

1 ln rðP; Q Þ; 2p

P; Q A O

ð44Þ

which is a particular singular solution of the equation

r2 C ¼ dðP; Q Þ yields after exploiting Eqs. (11) and (12) Z E 000 P efSS ðPÞ ¼  yx ln r fS dO G O Z Z  E cos a  yx0 yx00 ln rðy2 þz2 ÞdO þ fSS ðqÞ ds G r O G

ð45Þ

ð46Þ

a being the anticlockwise angle between r and n vectors, denoted as a ¼ rb;n (as shown in Fig. 1b); r = |P  q|, P, QAO, qAG and e = 2p, p or 0 depending on whether the point P is inside the domain O, P p on the boundary G or P outside O, respectively. It is worth noting here that the boundary has been assumed to be smooth at the point pAG. Applying once more the Green identity given by Eq. (43) for the following pair of functions   1 2 1 2 C¼ r ðln r  1Þ; r C ¼ ln r ð47aÞ 8p 2p j ¼ fPS C¼

1 4 ðy þ z4 Þ; ðr2 C ¼ y2 þz2 Þ 12

j ¼ ln r

ð47bÞ ð48aÞ ð48bÞ

the domain integrals of Eq. (46) can be converted into line integrals along the boundary of the cross-section and the integral

ARTICLE IN PRESS E.J. Sapountzakis, V.J. Tsipiras / International Journal of Non-Linear Mechanics 45 (2010) 63–74

4. Numerical examples

representation (46) is written as fSS ðPÞ ¼

e

Z E P  yx ðfS ðqÞð2lnr  1Þrcos a  ðzny  ynz Þðlnr  1Þr2 Þds 4G G  4   Z  y þz4P E 1 3 1 4 cos a ds  yx0 yx00 e P þ ðy ny þz3 nz Þlnr  ðy þz4 Þ 12 G 3 12 r G Z cos a þ fSS ðqÞ ð49Þ ds r G

The values of the function fSS inside the domain O can be established from the integral representation (49) if fSS is known on the boundary G. Eq. (49) written for the boundary points of the domain O (e = p) constitutes a system of simultaneous linear equations that can be solved with respect to the unknown boundary quantities fSS. Thus, using constant boundary elements to approximate the line integrals along the boundary and a collocation technique the following linear system of simultaneous algebraic equations is established: ½HS fFS g ¼ fF S g

ð50Þ

where fFS gT ¼

n

S

S

ðfS Þ1

ðfS Þ2

...

S

ðfS ÞN

o

ð51Þ

are the values of the boundary quantities fSS at the nodal points of the N boundary elements. Moreover, in Eq. (50) [HS] and {FS} are square N  N and column N  1 known coefficient matrices, respectively. Finally, the derivatives of fSS with respect to y and z at any interior point for the calculation of the stress resultants (Eqs. (7b) and (7c)) are computed by differentiating the integral representation (49) of the function fSS as S

@fS ðPÞ E m ¼ y @y 8pG x

Z G

P ðfS ðqÞð2cos

o cos a þ ð2lnr  1Þny Þ  ðzny  ynz Þð2lnr  1Þr cos oÞds

   Z  E 0 00 y3P 1 1 3 cos o 1 4 cosðo  aÞ ðy ny þ z3 nz Þ ðy þ z4 Þ  þ yx yx ds 2 G 12 3 2p G 3 r r Z 1 cosðo  aÞ þ fS ðqÞ ds ð52aÞ 2p G S r2



S

@fS ðPÞ E m ¼ y @z 8pG x

Z G

69

P

ðfS ðqÞð2sin o cos a þ ð2lnr  1Þnz Þ  ðzny  ynz Þð2lnr  1Þr sin oÞds

  Z  E 0 00 z3P 1 1 3 sin o 1 4 sinðo  aÞ ðy ny þ z3 nz Þ ðy þ z4 Þ  þ yx yx ds G 12 3 2p G 3 r r2 Z 1 sinðo  aÞ S þ f ðqÞ ds ð52bÞ 2p G S r2





with a ¼ rb;n; r =|P q|, PAO, qAG and o ¼ yb;r.

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 the efficiency, wherever possible accuracy and range of applications of the developed method. 4.1. Example 1 In the first example, for comparison reasons, a simply supported torsion member (E=2.0  108 kN/m2, G=8.0  107 kN/m2) of length l=1.0 m and of narrow rectangular cross-section 10  200 mm2 (It =6.6488  104 mm4, CS =5.499  107 mm6, In2 =1.7778  1010 mm6) subjected to a uniformly distributed torque per unit length mt has been studied (300 boundary elements and 50 nodal points as depicted in Fig. 2). In Fig. 3 the variation of the angle of twist yx at the midspan with respect to the applied external torque per unit length mt for zero value of the warping constant CS =0 for both the linear and the non-linear cases is presented as compared with the values obtained from a FEM solution (It =6.6667  104 mm4, In2 =1.7778  1010 mm6) [10]. Moreover, in Table 1 these values of yx are presented as compared with those obtained using the ‘‘correct’’ value of the warping constant. The accuracy of the proposed numerical procedure is demonstrated, the minor discrepancy of the results due to the vanishing warping constant for this narrow cross-section is observed, while the stiffening of the structure due to the geometrical non-linearity is remarked. 4.2. Example 2 In the second example, again for comparison reasons, the tested by Gregory [21-22] 901 monosymmetric angle section of flange midline length b= 14.61 mm and thickness t = 0.9601 mm cantilever bar of length l= 177.8 mm (E=8.966  107 kN/m2, G=3.113  107 kN/m2) twisted by a pure concentrated torque Mt applied at its free end has been studied (400 boundary elements and 50 nodal points as depicted in Fig. 4). In Table 2 the constants of the examined bar are shown as compared with those presented in [9]. Moreover, in Fig. 5 the shear center movement of the end cross-section for the cases of restrained warping at both bar ends, at one end or free warping at both ends is presented (vS(0)= wS(0) =0, v0 S(0)= w0 S(0) =0) as compared with the measured path by Gregory’s experiment [21,22] and the predicted one by a finite-element formulation [9], in which warping is free at both bar ends (uniform torsion). From this figure, it is observed that the accuracy of the obtained results employing the proposed method compared with those of the Gregory’s experiment and with those

Fig. 2. Cross-sectional (a) and longitudinal (b) mesh of the beam of example 1.

ARTICLE IN PRESS 70

E.J. Sapountzakis, V.J. Tsipiras / International Journal of Non-Linear Mechanics 45 (2010) 63–74

80

Table 2 Constants of the bar of example 2.

A (mm2) IYY (mm4) IZZ (mm4) IP (mm4) IPP (mm6) b1 (mm) yC (mm) It (mm4) CS (mm6) In2 (mm6)

40

3

mt(10 Ntmm/mm)

60

20

Present study

Attard [9]

28.06 999.2 252.2 1991 2.544  105  10.22 5.135 8.766 152.154 7.784  103

28.05 998.0 249.5 1996 2.556  105  10.33 5.165 8.620 – 7.102  103

Present study - Warping restrained at both ends Present study - Warping restrained, warping free

Present study with CS=0 Trahair (2005)

Present study - Warping free at both ends

Linear case

0

0 0

0.2

0.4 0.6 Angleof twist at x=l/2(rad)

0.8

0

1

4

8

12

y (mm)

Fig. 3. Torque–twist curves for the simply supported beam of example 1.

Gregory’s experiment

mt

CS = 5.499  107 mm6

CS = 0

10 20 30 40 50 60 70

0.2114 0.3618 0.4741 0.5643 0.6403 0.7063 0.7650

0.2142 0.3661 0.4798 0.5712 0.6483 0.7153 0.7748

-4

Attard’s FEM solution

z (mm)

Table 1 Angle of twist yx (rad) at the midspan of the simply supported beam of example 1 for various values of the uniformly distributed torque per unit length mt.

-8 Measured position at 75° twist -12 Fig. 5. Path of the shear center of the end cross-section for the cantilever beam of example 2.

Mt is presented for the cases of restrained warping at both bar ends, at one end or free warping at both ends. For all boundary conditions at the bar ends the stiffening of the structure arising from the geometrical non-linearity is observed, while this stiffening is slightly increased with the restraint of warping at one bar end and a bit more at both bar ends. 4.3. Example 3

Fig. 4. Cross-sectional (a) and longitudinal (b) mesh of the beam of example 2.

of the finite-element formulation [9] is very good, while the influence of the applied torsional boundary conditions at the bar ends is not significant. Finally, in Fig. 6 the variation of the angle of twist yx at the bar end with respect to the applied external torque

In the third example, in order to demonstrate the efficiency and the range of applications of the developed method, the nonsymmetric steel L-beam (E =2.1 108 kN/m2, G=8.0769  107 kN/m2) shown in Fig. 7a, of length l= 1.0 m clamped at its one end, torsionally elastically supported (k =200 kNm/rad) with free warping at the other one and twisted by a pure concentrated torque Mt applied at the midspan of the bar has been studied (400 boundary elements and 50 nodal points as depicted in Figs. 7b and c, respectively). In Table 3 the constants of the examined bar (y is the angle of the principal coordinate system through the crosssection’s centroid C, as shown in Fig. 7a) and in Fig. 8 the variation of the angle of twist yx at a point just before the midspan of the bar with respect to the applied external torque Mt are presented verifying again the stiffening of the structure arising from the

ARTICLE IN PRESS E.J. Sapountzakis, V.J. Tsipiras / International Journal of Non-Linear Mechanics 45 (2010) 63–74

Table 3 Constants of the bar of example 3.

Nonlinear - Warping restrained at both ends Nonlinear - Warping restrained, warping free

A (cm2) IYY (cm4) IZZ (cm4) IP (cm4) IPP (cm6) b1 (cm) b2 (cm) y (rad) yC (cm) zC (cm) It (cm4) CS (cm6) Uw (cm6) In (cm6) In2 (cm6)

Nonlinear - Warping free at both ends Linear - Warping free at both ends 8000

6000

Mt(Ntmm)

71

25.00 723.593 132.198 1445.518 168733.198 8.192 3.887 0.430 3.662 3.190 8.390 119.890  38.340 5937.214 5949.475

4000 200 Geometrically nonlinear case Geometrically linear case 2000

0 0

1

2 3 Angleof twist at x = l(rad)

4

Fig. 6. Torque–twist curves for the cantilever beam of example 2.

Mt(kNm)

160

120

80

40 t=1cm 0

h=15.5cm

0 Z

1 1.5 Angleof twist at x=l/2(rad)

2

2.5

Fig. 8. Torque–twist curves for the elastically supported beam of example 3.

θ

Y

C

S

0.5

t=1cm b=10.5cm

Fig. 7. L-shaped cross-section of unequal legs (a), cross-sectional (b) and longitudinal (c) mesh of the beam of example 3.

Table 4 Maximum and minimum values of the angle of twist yx and its derivatives along the bar of example 3.

yx (rad) y0 x (rad/cm) y00 x (rad/cm2) y000 x (rad/cm3)

Maximum values

Minimum values

2.1893 0.0482 0.0162 0.0026

0.0000  0.0401  0.0130  0.0041

geometrical nonlinearity. For the final loading step of this curve Mt = 170 kNm, in Table 4 the extreme values of the angle of twist yx and its derivatives along the bar are presented. Moreover, in Table 5 the distributions of the normal stresses Sxx and the flow of the resultant shear stress Sxn together with their maximum values at the clamped end and at the midspan of the bar are presented for Mt = 20 kNm. In Fig. 9 the resultant of the primary shear stress distribution SxyP, SxzP, that is the primary twisting moment MtP (second term of Eq. (30a) or Eq. (20d)), the resultant of the linear part of the secondary shear stress distribution SxyS, SxzS, that is the linear part of the secondary twisting moment MtS (first term of Eq. (30a)), the resultant of the non-linear part of the secondary

ARTICLE IN PRESS 72

E.J. Sapountzakis, V.J. Tsipiras / International Journal of Non-Linear Mechanics 45 (2010) 63–74

Table 5 Distributions of the normal stresses Sxx (a, b) and flow of the resultant shear stress Sxn (c, d) together with their maximum values at the clamped end and at the midspan of the bar (Mt = 20 kNm) of example 3. Clamped end

shear stress distribution SxyS, SxzS, that is the non-linear part of the secondary twisting moment MtS (third term of Eq. (30a)), the total twisting moment Mt (Eq. (30a)) and its total linear part (the first two terms of Eq. (30a)) are presented along the bar for

Midspan of the bar

Mt =80 kNm, while for the same loading step, in Fig. 10 the linear (first term of Eq. (30b)), the non-linear (second term of Eq. (30b)) and the total warping moment Mw, which is the resultant of the warping normal stress distribution Sxx, are also

ARTICLE IN PRESS E.J. Sapountzakis, V.J. Tsipiras / International Journal of Non-Linear Mechanics 45 (2010) 63–74

Mt=80kNm Angle of twist Lateral deflection at y direction - vs Lateral deflection at z direction - ws

1.6

1.2

40

0.8

40

20

y,z(cm)

Mt(kNm)

50

0

30

0.4

20 10

0

0 0

-20

-40 0

20

40

60

80

100

x(cm)

2

1

2

0

-1 Mt=80kNm Total Warping moment Linear warping moment Nonlinear warping moment

-2

-3 0

20

40

60 x (cm)

80

20

40 60 x(cm)

80

100

Fig. 11. Transverse displacement components of the shear center together with the variation of the angle of twist, along the bar of example 3.

torque–rotation relationship is computed based on the finitedisplacement (finite-rotation) theory, that is the transverse displacement components are expressed so as to be valid for large rotations and the longitudinal normal strain includes the second-order geometric non-linear term often described as the ‘‘Wagner strain’’. The main conclusions that can be drawn from this investigation are

Fig. 9. Twisting moments along the bar of example 3.

Mw(kNm )

Angle of twist (rad)

Mt=80kNm Total torque Total linear torque Uniform torque Warping torque Nonlinear torque

60

73

100

Fig. 10. Linear, non-linear and total warping moment Mw along the bar of example 3.

presented. From both of these figures the role of the non-linear part of the secondary twisting moment as well as that of the nonlinear warping moment is demonstrated. Finally, in Fig. 11 the transverse displacement components of the shear center vS, wS with respect to the undeformed shear center axis together with the variation of the angle of twist yx along the bar (vS(0) =wS(0) =0, v0 S(0)= w0 S(0) = 0) are presented for Mt = 80 kNm, demonstrating the lateral displacement of the shear center axis due to the geometrical non-linearity.

5. Concluding remarks In this paper an elastic non-uniform torsion analysis of simply or multiply connected cylindrical bars of arbitrary cross-section taking into account the effect of geometric non-linearity is presented employing the boundary-element method. The

a. accurate results are obtained using a relatively small number of boundary and domain elements along the cross-section and bar axis, respectively; b. geometrical non-linearity leads to stiffening of the structure and eventually its better response against torsional loading; c. as observed, restrained warping boundary conditions lead to slightly greater torsional resistance than those of free warping for both the geometrically linear and non-linear cases; d. bars of asymmetrical cross-section exhibit lateral displacement of their shear center axis under pure torsional loading, due to the geometrical non-linearity; e. the non-linear Wagner torque can reach significant values locally along the bar axis while the non-linear warping moment can be neglected in most cases; f. the secondary shear stress distribution and magnitude is influenced by the geometrical non-linearity and g. the developed procedure retains most of the advantages of a BEM solution over a FEM approach, although it requires domain discretization.

References [1] M.Z. Ahmed, F.E. Weisgerber, Torsion constant for matrix analysis of structures including warping effect, Int. J. Solids Struct. 33 (1996) 361–374. ¨ [2] F. Gruttmann, W. Wagner, R. Sauer, Zur Berechnung von Wolbfunktion und Torsionskennwerten beliebiger Stabquerschnitte mit der Methode der Finiten Elemente, Bauingenieur 73 (3) (1998) 138–143. [3] E.J. Sapountzakis, V.G. Mokos, Warping shear stresses in non-uniform torsion by BEM, Comput. Mech. 30 (2) (2003) 131–142. [4] E.J. Sapountzakis, V.G. Mokos, Warping shear stresses in non-uniform torsion of composite bars by BEM, Comput. Methods Appl. Mech. Eng. 192 (2003) 4337–4353. [5] D.G. Ashwell, The axis of distortion of a twist prism, Phil. Mag. 42 (1951) 820–832. [6] M. Gregory, The bending and shortening effect of pure torque, Austn J. Appl. Sci. 11 (1960) 209–216. [7] W.K. Tso, A.A. Ghobarah, A non-linear thin-walled beam theory, International Journal. J. Mech. Sci. 13 (1971) 1025–1038.

ARTICLE IN PRESS 74

E.J. Sapountzakis, V.J. Tsipiras / International Journal of Non-Linear Mechanics 45 (2010) 63–74

[8] M.M. Attard, Non-linear theory of non-uniform torsion of thin-walled open beams, Thin-Walled Struct. 4 (1986) 101–134. [9] M.M. Attard, I.J. Somervaille, Non-linear analysis of thin-walled open beams, Comput. Struct. 25 (3) (1987) 437–443. [10] N.S. Trahair, Non-linear elastic non-uniform torsion, J. Struct. Eng. 131 (7) (2005) 1135–1142. [11] F. Mohri, N. Damil, M.M. Potier-Ferry, Large torsion finite element model for thin-walled beams, Comput. Struct. 86 (7–8) (2008) 671–683. [12] J.T. Katsikadelis, The analog equation method. A boundary-only integral equation method for non-linear static and dynamic problems in general bodies, Theor. Appl. Mech. 27 (2002) 13–38. [13] G. Chen, N.S. Trahair, Inelastic non-uniform torsion of steel I-—beams, J. Constr. Steel. Res. 23 (1992) 189–207. [14] V.G. Mokos, Contribution to the generalized theory of composite beams structures by the boundary element method, Ph.D. Thesis for the PhD degree at the, NTUA, National Technical University of Athens 2007 (in Greek), 2007.

[15] K. Washizu, Variational Methods in Elasticity and Plasticity, Pergamon Press, Oxford, 1975. [16] E.J. Sapountzakis, Solution of non-uniform torsion of bars by an integral equation method, Comput. Struct. 77 (2000) 659–667. [17] E.J. Sapountzakis, V.G. Mokos, Non-uniform torsion of bars of variable cross section, Comput. Struct. 82 (9–10) (2004) 703–715. [18] E. Isaacson, H.B. Keller, Analysis of Numerical Methods, John Wiley and Sons, New York, 1966. [19] E.J. Sapountzakis, J.T. Katsikadelis, Unilaterally supported plates on elastic foundations by the boundary element method, J. Appl. Mech. Trans. ASME 59 (1992) 580–586. [20] J.T. Katsikadelis, Boundary Elements: Theory and Applications, Elsevier, Amsterdam-–London, 2002. [21] M. Gregory, A non-linear bending effect when certain unsymmetrical sections are subjected to a pure torque, Aust. J. Appl. Sci. 11 (1960) 33–48. [22] M. Gregory, The bending and shortening effect of pure torsion, Aust. J. Appl. Sci. 11 (1960) 209–216.