Nonuniform shear warping effect in the analysis of composite beams by BEM

Engineering Structures 76 (2014) 215–234

Contents lists available at ScienceDirect

Engineering Structures journal homepage: www.elsevier.com/locate/engstruct

Nonuniform shear warping effect in the analysis of composite beams by BEM I.C. Dikaros, E.J. Sapountzakis ⇑ 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 3 December 2013 Revised 19 May 2014 Accepted 8 July 2014

Keywords: Nonuniform shear warping Shear lag Shear deformation Composite beams Boundary element method

a b s t r a c t In this paper a general formulation for the nonuniform shear warping analysis of composite beams of arbitrary simply or multiply connected cross section with at least one axis of symmetry, under general boundary conditions is presented. The beam is subjected to the combined action of arbitrarily distributed or concentrated transverse loading passing through the shear center of the cross section, as well as to bending and warping moments. The nonuniform shear warping distribution is taken into account by employing one independent warping parameter multiplying a shear warping function, which is obtained by solving a corresponding boundary value problem, formulated exploiting the longitudinal local equilibrium equation. By taking into account the aforementioned nonuniform shear warping distribution, the developed formulation is capable of capturing ‘‘shear lag’’ effect on beams under flexure. Three boundary value problems are formulated with respect to the displacement and rotation components as well as to the independent warping parameter and solved using the Analog Equation Method, a Boundary Element Method based technique. The warping function and the geometric constants including the additional ones due to warping are evaluated employing a pure BEM approach. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction In engineering practice the analysis of beam members is frequently encountered. In the vast majority of these cases the assumptions of the well established Euler–Bernoulli theory play a dominant role on the structural analysis. Nevertheless, in cases where shear deformation is not negligible, this theory fails to give acceptable results and Timoshenko beam theory has to be employed. However, even within the context of the latter theory the assumption that plane cross sections remain plane after deformation is maintained and the warping due to shear is taken into account in an indirect manner through appropriate shear correction factors which can be employed with various procedures [1–3]. Nevertheless, in many structural members (e.g. beams of box-shaped cross sections, thin-walled beams or beams made of materials weak in shear), the shear effect becomes more complicated, since warping is more intense and its nonuniform distribution is responsible for the well-known shear lag phenomenon, which is associated with the significant modification of normal stress distribution, especially near the joint of the various cross sectional components [4]. As a consequence, elementary beam the⇑ Corresponding author. E-mail addresses: [email protected] (I.C. Dikaros), cvsapoun@central. ntua.gr (E.J. Sapountzakis). http://dx.doi.org/10.1016/j.engstruct.2014.07.009 0141-0296/Ó 2014 Elsevier Ltd. All rights reserved.

ories cannot predict accurately the actual structural behavior of the beam and may underestimate significantly the arising displacement and stress component values. Therefore, it is necessary to include the effects of nonuniform distribution of shear warping in the analysis. This kind of analysis becomes more complicated when beams composed by multiple materials (composite beams) are involved, which are used in many structural applications in an attempt to achieve a better distribution of rigidity and weight by exploiting the favorable characteristics of various structural materials. Shear lag phenomenon as a consequence of nonuniform shear warping distribution in box-shaped and folded structural members has been pointed out several decades ago and many efforts have been made for its analysis e.g. by Reissner [5] employing energy method approach or by Malcolm and Redwood [6] and Moffatt and Dowling [7] using finite elements [8]. In up to date regulations, in order to simplify the analysis as far as nonuniform shear warping of beams is concerned, the ‘‘effective’’ breadth concept [9–11] is employed. However, this simplifying approach may fail to capture satisfactorily the actual structural behavior of the member, since the influence of shear lag phenomenon is not constant along the beam length, while apart from the geometrical configuration of the cross section it depends also on the type of loading [12]. Many approaches towards analyzing shear lag phenomenon and estimating the effective breadth of box girders and stiffened

216

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

plates has been based either on folded plate theory assumptions [13–15] or on beam-plate models [16,17]. Moreover, more recently refined models based on shell or solid finite elements have been developed for the investigation of this phenomenon [18–20]. Nevertheless, it is true that these formulations require sophisticated and computationally demanding modeling. Hence, including the nonuniform shear warping in beam element formulations based on so-called ‘‘Higher-Order Beam Theories’’ [21] is of increased interest due to their important advantages over more refined approaches. Beam elements are practical, computationally efficient and offer better insight of the structural phenomena since they permit their isolation and their independent investigation. Furthermore, due to easy parameterization of all necessary data, beam elements are more convenient for parametric analyses than more refined models which, in most cases, require the construction of multiple models. During the past years various beam formulations for the nonuniform shear warping analysis of prismatic beams have been proposed, most of which employ the assumptions of thin-walled beam theory [22]. In the majority of these studies the shear warping is taken into account by ad-hoc assumptions using pre-defined polynomial curves to express the warping of the flanges of homogeneous or laminate composite box-shaped cross sections [4,23– 31]. In some of these studies the shear deformation of the webs is not taken into account, while the aforementioned assumptions generate restrictions as far as the cross sectional shape and the loading direction are concerned. Thin-walled homogeneous cross sections without the above simplifications have been examined in [32–36] by means of analytical solutions or the Finite Element Method (FEM). The nonuniform shear warping analysis of beams without adhering to the simplifying assumptions of thin-walled structures has received limited amount of literature. Among these studies, Massonnet [37] analyzed shear lag effect in general cross sections employing the simplifying assumption of a constant normal stress correction due to nonuniform shear warping. Ie and Kosmatka [12] performed nonuniform shear warping analysis employing firstorder warping functions and presented results for homogeneous beams of rectangular cross sections. Hoffmann [38] analyzed the nonuniform shear problem for rectangular and T-shaped homogeneous cross sections. Park et al. [39] presented a FEM formulation accounting for the nonuniform shear warping of general cross sections; however they focused exclusively on thin-walled beams. El Fatmi [40,41], El Fatmi and Ghazouani [21] and Ghazouani and El Fatmi [42,43] examined the case of the St. Venant beam problem including nonuniform shear and torsional warping of arbitrary homogeneous and composite cross sections. Finally, very recently Le Corvec and Filippou [44] presented a FEM formulation for elastic and elastoplastic analysis of beams of general cross section taking into account shear and torsional warping by employing a number of warping degrees of freedom on the cross sectional domain, while Ferradi et al. [45] developed a FEM formulation for the analysis of beams of arbitrary cross section including a multitude of shear and torsional warping functions. However, both of these latter formulations are developed for homogeneous cross sections. Moreover, to the authors’ knowledge a BEM procedure for the analysis of the nonuniform shear warping problem of composite beams has never been reported in the literature. The objective of this paper is to present a general boundary element formulation for the nonuniform shear warping analysis of composite beams of simply or multiply connected cross section, under general boundary conditions. It is noted that in the case of arbitrary cross sections several couplings between flexure and torsion may occur making the analysis more complicated [21]. Consequently, in the present study cross sections having at least one axis of symmetry will be studied. However, the proposed method can

be extended in order to account for any arbitrary cross sectional configuration with main purpose the application to framed structures. The composite beam (thin- or thick-walled) consists of materials in contact, each of which can surround a finite number of inclusions. The materials have different elasticity and shear moduli and are firmly bonded together. The beam is subjected to the combined action of arbitrarily distributed or concentrated transverse loading passing through the shear center of the cross section (torsionless bending), as well as to bending and warping moments. The nonuniform shear warping distribution is taken into account employing an independent warping parameter multiplying a shear warping function, which is obtained by solving a corresponding boundary value problem, formulated exploiting the longitudinal local equilibrium equation. By taking into account the aforementioned nonuniform shear warping distribution, the developed formulation is capable of capturing ‘‘shear lag’’ effect on beams under flexure. Three boundary value problems are formulated with respect to the displacement and rotation components as well as to the independent warping parameter and solved using the Analog Equation Method (AEM) [46], a BEM based method. The warping function, the additional geometric constants due to warping and the elementary ones are evaluated employing a pure BEM approach, i.e. only boundary discretization of the cross section is used. After the establishment of the kinematical components, the aforementioned boundary discretization can be employed to evaluate through BEM the normal and shear stress components at any arbitrary point of each cross section. The essential features and novel aspects of the present formulation are summarized as follows: i. The cross section is a monosymmetric arbitrarily shaped thin- or thick-walled composite one. The formulation does not stand on the assumption of a thin-walled structure and therefore the cross section’s warping rigidities are evaluated ‘‘exactly’’ in a numerical sense. ii. The beam is subjected to arbitrary external transverse torsionless loading and is supported by the most general boundary conditions including elastic support or restraint. iii. The proposed formulation is suitable for the investigation of shear lag effects in beams of closed or open cross sections. iv. The proposed method employs a BEM approach (requiring boundary discretization for the cross sectional analysis) resulting in line or parabolic elements instead of area elements of the FEM solutions (requiring the whole cross section to be discretized into triangular or quadrilateral area elements), while a small number of line elements are required to achieve high accuracy. v. The use of BEM permits the effective computation of derivatives of the field functions (e.g. stresses) which is very important during the nonuniform shear warping analysis of beams. Numerical examples of great practical interest are worked out to demonstrate the efficiency and the accuracy of the developed method compared with FEM results. In the aforementioned examples, the effects arising in the nonuniform shear warping analysis of composite beams are illustrated.

2. Statement of the problem 2.1. Global equations of equilibrium Consider a prismatic beam of length l, of arbitrarily shaped composite cross section having at least one axis of symmetry (Z axis), consisting of materials in contact, each of which can sur-

217

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

round a finite number of inclusions. The materials, occupying the simply or multiply connected regions Xm (m = 1, 2, . . . , M) of the YZ plane (Fig. 1), are firmly bonded together and are assumed homogeneous, isotropic and linearly elastic with modulus of elasticity Em, shear modulus Gm and Poisson ratio mm. Let also the boundaries of the non-intersecting regions Xm be denoted by CXm (m = 1, 2, . . . , M). These boundary curves are piecewise smooth, i.e. they may have a finite number of corners. In Fig. 1 CXYZ is the principal bending coordinate system through the cross section’s centroid C, while ZS is the coordinate of the cross section’s shear center S, with respect to the CXYZ system. The beam is subjected to the combined action of the arbitrarily distributed or concentrated transverse loading pz = pz(X) along Z direction, passing through the shear center S (torsionless bending), bending moment mY = mY(X) along Y direction (Fig. 1), as well as warping moment muP ¼ muP ðXÞ, which will be later defined. CY

CY

Under the action of the aforementioned general loading and of possible restraints, the beam is leaded to flexure with varying bending moment along its length (nonuniform bending). This bending moment represents the distribution of normal stresses due to bending (primary normal stresses rPxx ) at a cross section. Due to the aforementioned bending moment variation along the length of the beam, shear stresses arise on horizontal sections of an infinitesimal beam element (Fig. 2) equilibrating the variation of normal stresses due to bending. Cauchy principle dictates that corresponding shear stresses arise on the plane of the cross section as well. If the assumption that plane sections remain plane after deformation (Euler–Bernoulli or Timoshenko beam theories) is maintained, the arising shear stresses obtain a uniform distribution over the section. However, this distribution violates local equilibrium since the requirement of vanishing tractions sxn on the lateral surface of the beam is not satisfied. Thus, the aforementioned shear stresses exhibit a nonuniform distribution over the cross section’s domain so that both local equilibrium and vanishing tractions sxn on the lateral surface of the beam are satisfied. These nonuniform shear stresses will be referred to as primary (or St.Venant) shear

(a)

pz mY

y

mϕ P

Y

CY

z Z

l S C

C: Centroid S: Shear center

x X

pz ( X )

pz ( X ) P P σ xx + dσ xx

P σ xx

P τ xz P σ xx

P P σ xx + dσ xx

Y

X ≡x

P τ zx

Z

Z≡z

dx

x

Fig. 2. Development of primary shear stresses





sPxz over a rectangular cross section.

stresses sPxy ; sPxz and lead the cross section to warp (Fig. 3b). Furthermore, due to the nonuniform character of this warping along the beam length a secondary normal stress distribution rSxx is developed (Fig. 3c). This normal stress distribution is responsible for the well-known shear lag phenomenon and it is taken into account by employing an independent warping parameter multiplying the warping function, which depends on the cross section geometry. The nonuniform distribution of the secondary normal stresses along the length of the beam results in the development of the secondary shear stresses sSxy ; sSxz (Fig. 3d), which equilibrate the variation of the rSxx at an infinitesimal beam element. The sequence of stress generation previously described is schematically illustrated in Table 1. In order to take into account the aforementioned nonuniform warping behavior, in the present study an additional degree of freedom (DOF) is added to the well-known ones of the elementary beam theory. This additional DOF is a warping parameter gY(X) describing the ‘‘intensity’’ of the cross sectional warping along the length of the beam. This latter warping is defined by the warping function uPCY ðY; ZÞ, depending only on the cross sectional configuration. Thus, the ‘‘actual’’ deformed configuration of the cross section due to warping (Fig. 3b) is given as gY ðXÞuPCY ðY; ZÞ at any position along the X axis. The corresponding stress resultant of the aforementioned additional DOF is a warping moment MuP CY (bimoment) along the beam length, arising from the secondary S normal stress distribution rxx the variation of which is equilibrated by a corresponding secondary stress resultant Q Sz arising from the secondary sSxz ; sSxy shear stress distributions. The aforementioned bimoment M uP due to shear warping constitutes an additional CY ‘‘higher order’’ stress resultant, which is developed in the nonuniform shear theory, in analogy to the warping moment (bimoment) of the nonuniform torsion theory [47]. Within the context of the above considerations, the displacement components of an arbitrary point of the beam are given as

 ðX; Y; ZÞ ¼ u  P ðX; Y; ZÞ þ u  S ðX; Y; ZÞ u ¼ hY ðXÞZ þ gY ðXÞuPCY ðY; ZÞ

(b)

ð1aÞ

t n a

Γ Ωm = Γ= Ω=

En, Gn r = q − P

k Γi i =1 Ωm M Γ m =1 Ω m M Ω m=1 m

y

E1, G1

S

ZS Y

Γ0Ω1

Γ1Ω1

C 2 ΓΩ 1 z

Ωm

v ðX; Y; ZÞ ¼ 0

ω

Ω1

Γ3Ω1 Em, Gm

Z Fig. 1. Prismatic beam under flexural loading (a) of a monosymmetric composite cross-section occupying the two dimensional region X (b).

 wðX; Y; ZÞ ¼ wðXÞ

ð1b; cÞ

; w  are the total axial and transverse beam displacement where u components with respect to the CXYZ centroidal system of axes; P ; u  S are the corresponding primary and secondary axial displaceu ment components; w(X) describes the deflection of the centroid C, while hY(X) is the angle of rotation due to bending about the centroidal Y axis. In Eq. (1a) uPCY ðY; ZÞ is the shear warping function with respect to the centroid C, while gY(X) is the dimensionless parameter introduced to take into account the nonuniform distribution of warping due to shear (‘‘intensity’’ of the cross sectional warping). It is worth here noting that by employing hY(X) instead of w,x(X) as measure of the angle of bending rotation, the inconsistency of

218

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

P σ xx

(a)

P τ xz

=

+

(b)

Primary shear stress distribution

Primary normal stresses due to bending P ( σ xx = EθY , x Z )

P P = Gγ ZP ΦCY ( τ xz ,z )

S σ xx S τ xz

(c)

(d)

Secondary shear stress distribution S P = Gγ ZSϕCY ( τ xz ,z )

Secondary normal stresses due to warping S P ) ( σ xx = EηY , xϕCY

Fig. 3. Sequence of stress generation over the height of a rectangular cross section of a beam under flexure. Table 1 Sequence of stress generation in nonuniform shear problem.

rPxx

equilibrium of variation









! P @ sP @ rP xy @ sxz xx @x þ @y þ @z ¼0



sPxy sPxz



Generation





!

rSxx

of warping

equilibrium of variation









! S @ sS @ rS xy @ sxz xx @x þ @y þ @z ¼0



sSxy sSxz



where cPZ ¼ w;x þ hY ; cSZ ¼ gY  w;x  hY are considered as primary and secondary ‘average’ shear strain quantities and

UPCY ¼ Z þ uPCY the vanishing of shear stresses in boundary conditions [39] is avoided. Employing the strain–displacement relations of the threedimensional elasticity for small displacements, employing the displacement expressions of Eq. (1) and applying the Hooke’s stress-strain law, the non-vanishing components of the Cauchy stress tensor in the regions Xm (m = 1, 2, . . . , M) are obtained as

h    i  P;x þ u  S;x ðrxx Þm ¼ Em ðexx Þm ¼ Em u m m h i   ¼ Em hY;x Z þ gY;x uPCY m ; m ¼ 1; 2; . . . ; M

 ;y Þm  ðsxy Þm ¼ Gm ðcxy Þm ¼ Gm ½ðv ;x Þm þ ðu   P ¼ Gm gY uCY;y ; m ¼ 1; 2; . . . ; M m    ;x Þm þ ðu  ;z Þm ðsxz Þm ¼ Gm ðcxz Þm ¼ Gm ðw h i ¼ Gm ðw;x þ hY Þ þ gY uPCY;z ; m ¼ 1; 2; . . . ; M

ð2aÞ

ð2bÞ

secondary

In order to establish the differential equations of equilibrium, the principle of virtual work is employed as

Z V

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

ð2cÞ

Z

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

ð5Þ

Lat

In the above equation, d() denotes virtual quantities, V is the volume of the beam in the initial undeformed state and F stands for the lateral surface of the beam including end cross sections. Moreover, tx, ty, tz are the components of the traction vector applied on the lateral surface of the beam including the end cross sections. The stress resultants of the beam can be defined as MY ¼

M X

ðMY Þm ¼

m¼1

where (),i denotes differentiation with respect to i. Taking into account the definition of primary and secondary normal and shear stresses given in the above description of stress generation sequence (Table 1) the non-vanishing components of the stress tensor (2) can be written as   ðrxx Þm ¼ Em hY;x Z þEm gY;x uPCY m ; m ¼ 1;2;...;M ð3aÞ |fflfflfflffl{zfflfflfflffl} |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} primary secondary     ð3bÞ ðsxy Þm ¼ Gm cPZ UPCY;y þGm cSZ uPCY;y ; m ¼ 1;2;...;M m m |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} primary secondary     ð3cÞ ðsxz Þm ¼ Gm cPZ UPCY;z þGm cSZ uPCY;z ; m ¼ 1;2;...;M m m |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} primary

ð4Þ

M uP ¼ CY

Q Pz ¼

M  X

M uP

CY

m¼1 M  X

M  X

Q Pz

m¼1

Q Sz



 m

 m

¼

m

¼



rxx ZdX

M Z X

m¼1

h

Xm

M Z X m¼1



Xm

m¼1

M Z X

¼

ð6aÞ

Xm

m¼1

m¼1

Q Sz ¼

M Z X

Xm

rxx uPCY dX

ð6bÞ

    i ðsxy Þm UPCY;y þ ðsxz Þm UPCY;z dX m

h

m

ð6cÞ

    i ðsxy Þm uPCY;y þ ðsxz Þm uPCY;z dX ð6dÞ m

m

where MY is the bending moment and M uP is the warping moment CY (bimoment). Q jz ðj ¼ P; SÞ are the primary and secondary parts of total shear force Qz. It is noted that the secondary shear force is also referred to as bishear stress resultant [32] since it equilibrates the corresponding warping moment (bimoment). Moreover, the geometric constants of the beam can be obtained by the following definitions

219

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

Z M M X X Em ðAÞm ¼ km dX E Xm m¼1 ref m¼1 Z M M X X Em SY ¼ ðSY Þm ¼ km ZdX E Xm m¼1 ref m¼1 Z M M X X Em SZ ¼ ðSZ Þm ¼ km YdX E Xm m¼1 ref m¼1 Z M M X X Em IYY ¼ ðIYY Þm ¼ km Z 2 dX E Xm m¼1 ref m¼1 Z M M X X Em IZZ ¼ ðIZZ Þm ¼ km Y 2 dX E ref X m m¼1 m¼1 Z M M X X Em IYZ ¼ ðIYZ Þm ¼ km YZdX E Xm m¼1 ref m¼1 Z M M  X X Em  SuP ¼ SuP ¼ km uPCY dX CY CY m E Xm m¼1 ref m¼1 Z M M X X Em Iij ¼ ðIij Þm ¼ km ðiÞm ðjÞm dX; i ¼ Z; uPCY ; j ¼ uPCY E Xm m¼1 ref m¼1 A¼

Dij ¼

ð7aÞ ð7bÞ ð7cÞ ð7dÞ

Q z;x ¼ pz

ð7fÞ

where the externally applied loads are related to the components of the traction vector applied on the lateral surface of the beam tx, ty, tz as

ð7gÞ

pz ðXÞ ¼

i ¼ Z; U

u

P CY ;

P CY ;

j¼U

u

t z ds mY ðXÞ ¼

C

Z

C

 M uP

CY

;x

t x Zds muP ðXÞ ¼

CY

Z

CY

C

ð12a; b; cÞ

t x uPCY ds ð13a; b; cÞ

ð7hÞ

Employing expressions (9), (10a) and (11), the governing differential equilibrium equations for the beam in terms of its kinematical components can be written as

 Gref Aðw;xx þ hY;x Þ þ Gref ASZ gY;x ¼ pz

P CY

ð7iÞ

where r  (),yiY + (),ziZ is the gradient operator and iY, iZ the unit vectors along Y, Z axes, respectively. A is the transformed area of the cross section, SY, SZ are the first moments of inertia, and IYY, IZZ, IYZ are the second moments and the product of inertia, respectively, of the composite cross section. Moreover, SuP is the first CY warping moment of inertia, IuP uP is the warping constant and   CY CY P P Dij i; j ¼ Z; UCY ; uCY are geometric constants associated with the resistance of the cross section in shear. Finally, Eref, Gref are the modulus of elasticity and shear modulus of a reference material, respectively. It is worth here noting that any material of the composite cross section can be used as reference material. Moreover, the weighted elastic and shear moduli of the m-th material have the same value km, by neglecting the influence of Poisson ratio. Employing definitions (7), having in mind that the coordinate system is the principal bending centroidal one and the orthogonality conditions of the warping function, which will be discussed in section 2.2, the following relations are obtained

SY ¼ SZ ¼ SuP ¼ IYZ ¼ IZuP ¼ 0 CY

Z

 MY;x þ Q z ¼ mY

 Q Sz ¼ muP

ð7eÞ

Z M M X X Gm ðDij Þm ¼ km ½rðiÞm  rðjÞm dX; G Xm m¼1 ref m¼1 P CY ;

corresponds to shear rigidity of Timoshenko beam theory. Thus, the simplified notation Gref APZ could be adopted for this quantity. Similarly, Gref DuP uP refers to secondary shear rigidity due to nonCY CY uniform shear warping and can be denoted as Gref ASZ . In what follows, in order to maintain the compatibility with classical notations, the above simplified symbols will be employed. Using the expressions of strain components, the definitions of the stress resultants (Eq. (6)) and applying the principle of virtual work (Eq. (5)), the differential equations of equilibrium of the beam can be derived as

ð8Þ

CY

P CY

Furthermore, employing the definition of warping function U (Eq. (4)) and its boundary condition which will be presented in section 2.2 and applying the Gauss-Green theorem the following relation is obtained

gY ¼ mY gY;xx þ Gref ASZ ðgY  w;x  hY Þ ¼ muPCY

 Eref IYY hY;xx þ Gref Aðw;x þ hY Þ   Eref IuP

uP CY CY

Gref ASZ

ð14aÞ ð14bÞ ð14cÞ

The above differential equations (Eq. (14)) are subjected to the corresponding boundary conditions of the problem at hand, which are given as

a1 w þ a2 V bz ¼ a3 a 1 hY þ a 2 MbY ¼ a 3 a~ 1 gY þ a~ 2 MbuPCY ¼ a~ 3

ð15aÞ ð15bÞ ð15cÞ

at the beam ends x = 0, l, where the reaction forces and moments Vbz, MbY, MbuP are given by relations (11), when applied at the beam CY ends. k ; a ~ k ðk ¼ 1; 2; 3Þ are functions specified at the Finally, ak ; a boundaries of the beam (x = 0, l). The boundary conditions (15) are the most general boundary conditions 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) can be derived from these equations by specifying appropriately these functions (e.g. for a clamped edge it is a1 ¼ a 1 ¼ a~ 1 ¼ 1; a2 ¼ a3 ¼ a 2 ¼ a 3 ¼ a~ 2 ¼ a~ 3 ¼ 0Þ. 2.2. Warping function uPCY

ð10a; b; cÞ

The analysis described in the previous section presumes that the shear warping function uPCY ðy; zÞ is already established. This can be accomplished by establishing UPCY ðy; zÞ (Fig. 4a). This latter warping function can be evaluated by exploiting the local equilibrium equation in the longitudinal direction and the corresponding boundary condition as

Employing Eqs. (3), (6)–(10) the expressions of the stress resultants in terms of the kinematical components can be obtained as

ðrxx;x Þm þ ðsxy;y Þm þ ðsxz;z Þm ¼ 0 ð16aÞ on the free surface of the beam ð16bÞ ðsxy Þm nY þ ðsxz Þm nZ ¼ 0

DuP

CY

uPCY

¼ DZuP

ð9Þ

CY

Consequently, employing Eqs. (4) and   i ¼ UPCY ; j ¼ Z; UPCY ; uPCY can be written as

DUP

CY

UPCY

¼ A  DuP

CY

M Y ¼ Eref IYY hY;x Q Pz ¼ Gref DUP

UP CY CY

DZUP ¼ DUP

uPCY

CY

CY

M uP ¼ Eref IuP CY

CY

UPCY

DUP

CY

(9),

uPCY

¼0

uPCY gY;x

constants

Dij

ð11a; bÞ

cPZ Q Sz ¼ Gref DuPCY uPCY cSZ

ð11c; dÞ

Having in mind Eqs. (9) and (10b), it is worth noting that the given definitions of the primary (6c) and secondary (6d) parts of the shear PM R  P  P force coincide  with  the classical ones Q z ¼ m¼1 Xm sxz m dXP; PM R S S Q z ¼ m¼1 Xm sxz m dX. Moreover, Gref DUP UP multiplying cZ CY

CY

ðsxy Þm nY þ ðsxz Þm nZ ¼ ðsxy Þn nY  ðsxz Þn nZ on the interfaces ð16cÞ where nY, nZ are the direction cosines with respect to Y,Z axes. Examining UPCY , stress field (3) is employed, ignoring at this stage   the secondary shear stresses sSxy ; sSxz and the normal stress due   to warping rSxx . According to the stress generation sequence described above (Table 1), primary shear stress components   sPxy ; sPxz equilibrate normal stresses due to bending rPxx (St. Venant

220

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

(a)

(b)

Fig. 4. Warping functions U

P CY

problem of uniform shear). Hence, employing (16), (11a,c), having in mind that equilibrium Eq. (12b) is valid (with the absence of secondary shear stress counterparts) and considering that mY = 0 in order to avoid evaluating load-dependent warping functions, a boundary value problem is formulated as



r2 UPCY



m

n

ð17aÞ 3. Integral representations – numerical solution

ð17bÞ ð17cÞ

  where UPCY ¼  IYY =APZ UPCY and (),n indicates the derivative with respect to the outward normal vector to the boundary n (directional derivative). After establishing UPCY from Eq. (17), the primary shear area APZ can be evaluated employing definition (7i) as

ðIYY Þ2 DUP UP CY

ð18Þ

CY

At this point the shear deformation coefficient aZ of Timoshenko beam theory (and the corresponding shear correction factor jZ) can be obtained (for comparison purposes at the numerical examples section) following the definition given in [3] and taking into account relation (18) as

aZ ¼

1

jZ

¼

constant indicating a parallel displacement of the cross section along the beam axis. The evaluation of this constant is performed as presented in [47]. After its evaluation, this constant is subtracted from the warping function leading to the satisfaction of the orthogR  P  P onality condition h1; uPCY i  SuP ¼ M m¼1 km Xm uCY m dX ¼ 0. CY

¼Z

 m ¼0 on the free surface of the beam Gm UPCY;n  m   Gm UPCY;n ¼ Gn UPCY;n on the interfaces

APZ ¼

(a) and warping functions uPCY (b).

Z  M 2  2  A X A P P dX ¼ P k U þ U m CY;y CY;z I2YY m¼1 AZ Xm

According to the precedent analysis, the nonuniform shear problem of composite beams reduces in establishing the components w(x), hY(x), gY(x) having continuous derivatives up to the second order with respect to x at the interval (0, l) and up to the first order at x = 0, l, satisfying the boundary value problem described by the coupled governing differential equations of equilibrium (Eq. (14)) along the beam and the boundary conditions (Eq. (15)) at the beam ends x = 0, l. Eqs. (14) and (15) are solved using the Analog Equation Method [46]. According to this method, let w(x), hY(x) and gY (x) be the sought solutions of the aforementioned problem. Setting as u1(x) = w(x), u2(x) = hY(x) and u3(x) = gY(x) and differentiating with respect to x these functions two times, respectively, yields

ui;xx ¼ qi ðxÞ; ð19Þ

Moreover, warping function uPCY (Fig. 4b) defining normal stress pattern due to warping can be established through Eq. (4). It is remarked that employing (4) and (18) and applying the GaussGreen theorem, it can be proved that uPCY fulfills the orthogonality  P  R P condition hZ; uPCY i  IZuP ¼ M m¼1 km Xm Z uCY m dX ¼ 0, which indiCY

cates that normal stresses due to uPCY do not generate additional bending moments. It is remarked that boundary value problem   R P P Eq. (17) has solution, since the condition M ds m¼1 CX Gm UCY;n m

3.1. Evaluation of kinematical components

m

¼ 0 ðm; n ¼ 1; . . . ; MÞ [48] is satisfied. It is also worth here noting that, since the problem (17) has a Neumann type boundary condition, the warping function contains an arbitrary integration

ði ¼ 1; . . . ; 3Þ

ð20Þ

Eq. (20) indicate that the solution of Eqs. (14) and (15) can be established by solving Eq. (20) under the same boundary conditions (Eq. (15)), provided that the fictitious load distributions qi (i = 1, . . . , 3) are first established. These distributions can be determined using BEM. The solutions of Eq. (20) are given in integral form as [49]

ui ðxÞ ¼

Z

l

u ðg; xÞqi ðgÞdg

0

h if¼l  u ðf; xÞui;f ðfÞ  u;f ðf; xÞui ðfÞ ; f¼0

i ¼ 1; . . . ; 3

ð21Þ

with x being a source point (x 2 (0, l)) and g an arbitrary moving point in the interval (0, l) of the beam. When point g reaches the beam ends 0, l, it is denoted as f. u⁄(g, x) is the fundamental solution given as

221

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

Xi: Collocation point

where q = r/l. Notice that in Eq. (24), for the line integral it is r = g  x, x, g points inside the beam, whereas for the rest of the terms it is r = f  x, x inside the beam, f at the beam ends 0, l. Differentiating Eq. (24) with respect to x (x 2 (0, l)), results in the integral representations of the derivatives of ui as

2

qj

1

qj 2

1

qi

qi

q(X )

i-th element

j-th element X

1

X=0

Xi

Xj

Semi-discontinuous element

Fig. 5. Discretization of the beam in linear continuous elements.

Z

A

pz=-0.1kN/m h=1.0 m

X

C

Y

l=3.0 m

b=0.1 m

(a)

(b)

with r = g  x, which is a particular solution of the equation

ð23Þ

where d(g  x) is the Dirac function in one dimension. Employing Eq. (22), the integral representation (21) can be written as

ui ðxÞ ¼

Z 0

l

2

E12 E22 0

ð22Þ

u;xx ¼ dðg  xÞ



K2 ðrÞqi ðgÞdg  K2 ðrÞui;f ðfÞ  K1 ðrÞui ðfÞ

i ¼ 1; . . . ; 3

f¼l f¼0

where the kernels Kj(r) (j = 1, 2) are given as

ð25a; bÞ

i ¼ 1; . . . ; 3

i ¼ 1; . . . ; 3

ð26aÞ ð26bÞ

8 9 0 > > > > > > > > 38 9 > a3 > > > > > E13 > d > > 1 < = <0 > = 7 0 5 d2 ¼ > >  > a3 > : ; > > > > E33 d3 > > > >0 > > > > > : > ; a~ 3

ð27Þ

where       F E1 E2 F E1 E2 F E1 E2 E22 ¼ E33 ¼ E11 ¼ 0 D21 D22 0 D41 D42 0 D63 D64

ð28a;b;cÞ

; ð24Þ

1 1 1 1 K1 ðrÞ ¼ sgnr ¼ sgnq K2 ðrÞ ¼ j r j¼ l j q j 2 2 2 2

K1 ðrÞqi ðgÞdg þ ½K1 ðrÞui;f ðfÞf¼l f¼0 ;

The integral representations (26a) written for the beam ends 0, l together with the boundary conditions (15) can be employed to express the unknown boundary quantities ui, ui,x in terms of qi (i = 1, . . . , 3). This is accomplished numerically as follows. The interval (0, l) is divided into L elements, on which qi (x) is assumed to vary according to a certain law. In this study, linear continuous elements are employed (Fig. 5). It is noted here that since fictitious loads qi (i = 1, . . . , 3) are not defined on the boundary points, special care is needed for the extreme elements in order to avoid the coincidence of the first and the last node with the beam ends (Fig. 5). In this case appropriate shape functions for discontinuous elements [50] have to be employed for the integration of the kernels. Thus, the following set of 12 algebraic equations is obtained

E11 6 4 0 0

Fig. 6. Side view (a) and cross section (b) of the beam of example 1.

u ðg; xÞ ¼ 1=2 j r j

l

ui;xx ðxÞ ¼ qi ðxÞ;

l

Z

Z 0

X=l

r

1

ui;x ðxÞ ¼ 

E12 ¼



0

0

0 D23

0 D24



E13 ¼



0

0

0 D27

0



D28

ð28d; eÞ

The matrices D21 to D64, are 2  2 rectangular known matrices j ; a ~ j (j = 1,2) of Eq. (15); including the values of the functions aj ; a a3, a 3 ; a~ 3 are 2  1 known column matrices including the boundary 3 ; a ~ 3 of Eq. (15); E1, E2 are rectangular values of the functions a3, a

0.45

0.09

0.4

0.08

0.35 0.3

0.07

0.25 0.2

0.05

0.15

0.03

0.1 0.05

0.02

0

0

-0.05

-0.01

-0.1

-0.02

-0.15 -0.2

-0.03

-0.25

-0.05

-0.3 -0.35

-0.06

-0.4

-0.08

-0.45

-0.09

0.06 0.04

0.01

-0.04

-0.07

(ΦCYP )max = 4.1667 E − 01 m

(ϕCYP )max = 8.3574E − 02 m

(a)

(b)

Fig. 7. Distribution of UPCY (a) and uPCY (b) over the cross sectional domain of the beam of example 1.

222

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

2  2 known coefficient matrices resulting from the values of kernels at the beam ends; F is 2  (L + 1) rectangular known matrix originating from the integration of kernels on the axis of the beam. Finally



di ¼

qi

^i u

^ i;x u

T

;

ð29Þ

ði ¼ 1; . . . ; 3Þ

Table 2 Geometric constants of the beam of example 1. A = 1.0000E01 m2

jZ = 8.3331E01

APZ ¼ 8:3331E  02 m2

IYY = 8.3333E03 m4

ASZ ¼ 1:6667E  02 m2

IuP uP ¼ 9:9201E  05 m4 CY CY

are generalized unknown vectors, where

  ^ i ¼ f ui ð0Þ ui ðlÞ gT ; ði ¼ 1;...;3Þ u ^ i;x ¼ ui;x ð0Þ ui;x ðlÞ T ; ði ¼ 1;...;3Þ u ð30a;bÞ are vectors including the two unknown boundary values of the  T respective boundary quantities and qi ¼ q1i q2i    qLþ1 ; i ði ¼ 1; . . . ; 3Þ are vectors including the L + 1 unknown nodal values of the fictitious loads. Discretization of the integral representations of the unknown quantities ui (i = 1, . . . , 3) inside the beam (x 2 (0, l)) and application to the L + 1 collocation nodal points yields

ui ¼ Hdi ; ði ¼ 1; . . . ; 3Þ ui;x ¼ H;x di ; ði ¼ 1; . . . ; 3Þ

ð31a; bÞ

where H,H,x are (L + 1)  (L + 5) known coefficient matrices. Applying Eq. (14) to the L + 1 collocation points and employing Eq. (32), 3  (L + 1) algebraic equations of equilibrium are formulated as

Kd ¼ f

ð32Þ

Fig. 8. Deflection w(X) of the clamped beam of example 1.

where K, f are generalized stiffness matrix and force vector, respecT tively, while d ¼ ½ d1 d2 d3  . Eq. (32) with Eq. (27) form a system of 3  (L + 1) + 12 equations with respect to the generalized unknown vector d. Solving the linear system for the fictitious load ^i; u ^ i;x distributions qi and the unknown boundary quantities u (i = 1, . . . ,3), the displacements and their derivatives in the interior of the beam are computed using Eq. (31). 3.2. Evaluation of the warping function UPCY The numerical solution for the evaluation of the kinematical components presented in the previous section assumes that the constants defined in Eq. (7h,i) are already computed. According to the analysis presented in Section 2.2, the evaluation of the aforementioned constants reduces in establishing warping function UPCY at any interior point of each domain Xm (m = 1, . . . ,M) of the cross section having continuous partial derivatives up to the second order with respect to Y, Z, satisfying the boundary value problem described by Eq. (17). The evaluation of the warping function UPCY is accomplished using BEM within the context of the method of subdomains [50,51]. According to this method, the second Green’s identity

Z

h

Xm

¼

   i GðP;Q Þ r2 UPCY ðQ Þ  r2 GðP;Q Þ UPCY ðQ Þ m dX

Z



m

h

CXm



GðP;qÞ UPCY;n ðqÞ

 m

  i  G;n UPCY ðqÞ m ds; m ¼ 1;2;...;M

when applied on every subdomain m = 1, 2, . . . , M to the warping function UPCY and to the fundamental solution

1 ln rðP; Q Þ 2p

ð34Þ

ð35Þ

where d(P, Q) is the Dirac function in two dimensions with P(y, z), Q(n, g) (when Q is located on the boundary CXm it is denoted as q) being arbitrary points of the domain and after converting the

ð36Þ

where the parameter e = e(P) = 1, 1/2 or 0 depending on whether point P is in the domain Xm, on the boundary CXm (P  p) or outside 1=2

Xm, respectively, r ¼ rðP; Q Þ ¼j Q  P j¼ ½ðn  xÞ2 þ ðg  yÞ2  , nZ is the direction cosine with respect to Z axis. Furthermore, kernels Ki(r) (i = 1, . . . , 4) are given as

r ;n K2 ðrÞ ¼ ln r r

1 1 1 rr;n K4 ðrÞ ¼ ðln r  1Þr 2 ln r  K3 ðrÞ ¼ 2 2 4

K1 ðrÞ ¼ 

which is a singular particular solution of the equation

r2 GðP; QÞ ¼ dðP; QÞ

remaining domain integral into a line one along the boundary CXm [48,3], yields Z Z h    i    1 1 ds; e UPCY m ¼ ðK3 Z  K4 nZ Þds K1 UPCY m þ K2 UPCY;n m 2p CXm 2 p C Xm

m ¼ 1;2;...;M

ð33Þ

GðP; Q Þ ¼

Fig. 9. Axial displacement components at the point A of the clamped beam of example 1.

ð37a; bÞ

ð37c; dÞ

It is noted that UPCY and its derivatives with respect to Y, Z at every point P of the domain Xm can be obtained by employing relation

223

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234 Table 3  PA ðXÞ ðmÞ; u  SA ðXÞ ðmÞ and u  A ðXÞ ðmÞ of the beam of example 1. Extreme values of the kinematical components wðXÞ ðmÞ; hY ðXÞ ðradÞ; u

(w)min (h )  YP max  u  AS max  A max u  A Þmax ðu

Present study

FEM Euler–Bernoulli beam model [55]

FEM Timoshenko beam model [55]

FEM solid model [55]

2.8010E05 1.2818E05 6.4090E06

1.2656E05 1.2932E05 6.4660E06

2.8543E05 1.2932E05 6.4660E06

2.8021E05 – –

1.2535E06

–

–

–

7.4454E06

6.4660E06

6.4660E06

7.3524E06

0.5

1.8 1.4 1 0.6

0

0.2

(a)

-0.2 -0.6 -1

-0.5

-1.4

0

0.5

1

2

1.5

(τ xz )max = 1.760 kN

2.5

3

-1.8

m2

(b)

(τ xz )max = 1.669 kN

m2

0.5

2 1.6 1.2 0.8 0.4

0

0

(c)

0

-0.4 -0.8

-0.5

-1.2

0

0.5

1

1.5

2

(τ xz )max = 1.800 kN

2.5

3

-1.6 -2

m2

Fig. 10. Contour lines of sxz (kN/m2) over the mid-plane of the clamped beam of example 1, employing the present study (a), a FEM solution using 4800 solid elements (b) and Timoshenko beam theory (c).

(36) and its derivatives with respect to Y, Z [3]. The unknown quan    tities UPCY ðpÞ m and UPCY;n ðpÞ on the boundary can be evaluated m

from the solution of a boundary integral equation on the boundary CXm , which is derived working as follows. Considering a point p and a point q lying on the boundary CXm (m = 1, 2, . . . , M), Eq. (36) may be written as



p UPCY

 m

¼

Z CXm

m ¼ 1; 2; . . . ; M

ðK3 Z  K4 nZ Þds 

Z CXm

h

  i   K1 UPCY m þ K2 UPCY;n ds; m

ð38Þ

In order to approximate the integrals appearing in (38), the bound  ary CXm is divided into Nm boundary elements, on which UPCY m ,   UPCY;n are assumed to vary according to a certain law (constant, m

linear, parabolic, etc.). Thus, using constant boundary elements and the collocation technique, the following system of linear algebraic equations is obtained for each subdomain as

b ¼ f m  Gm ð U b ;n Þ ; Hm ð UÞ m m

ðm ¼ 1; . . . ; MÞ

ð39Þ

where Hm, Gm are Nm  Nm known coefficient matrices originating from the integration of the kernels of Eq. (38) over the boundary

224

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

0.5

0

-0.5

0

0.5

1

1.5

2

(σ VM )max = 7.396 kN

2.5

3

7.5 7 6.5 6 5.5 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0

(a)

m2

(b)

(σ VM )max = 6.822 kN

m2

0.5

0

-0.5

0

0.5

1

1.5

2

(σ VM )max = 5.474 kN

m

2.5

3

5.4 5 4.6 4.2 3.8 3.4 3 2.6 2.2 1.8 1.4 1 0.6 0.2

(c)

2

Fig. 11. Contour lines of rVM (kN/m2) over the mid-plane of the clamped beam of example 1, employing the present study (a), a FEM solution using 4800 solid elements (b) and Timoshenko beam theory (c).

curve. It is noted that for the computation of these integrals, standard Gauss–Legendre numerical integration scheme has been employed since it facilitates implementation without compromising the accuracy. fm is a Nm  1 known vector containing the values of the first integral in the right hand side of Eq. (38). Moreover

b ¼ ð UÞ m

h

b ;n Þ ¼ ðU m

UPCY



1



m

UPCY;n

UPCY

1

2 m



m



UPCY;n

2 m



UPCY



 N m iT m



UPCY;n

ð40aÞ Nm T

ð40bÞ

m

where G1m, G2m are Nm  (Nm  Km) and Nm  Kmcoefficient   matrib int , U b fr ces, respectively, originating from Gm, while U are the ;n ;n m m vectors containing the values of the directional derivative of warping function on the Km interface nodes and on the Nm  Km free nodes, respectively, of the m-th boundary. Assembling appropriately Eq. (41) for all of the subdomains of the composite cross section, yields

( ½ H G2 

b U b U int

) b fr ¼ f  G1 U ;n

ð42Þ

;n

are vectors containing the values of the warping function and its directional derivative, respectively, on every nodal point of the boundary element mesh of region m. It is noted here that in an interface between two adjacent subdomains, both the value of the warping function and of its directional derivative are unknown. Thus, considering the number of boundary elements of the possible interfaces of the m-th boundary as Km 6 Nm, Eq. (39) are written as 8   9 > > b < U =   b fr ; ðm ¼ 1; . . . ; MÞ m ½ Hm G2m   ¼ f m  G1m U ð41Þ ;n > m b int > : U ; ;n m

which is a system of N linear algebraic equations with respect to P N + K unknown boundary quantities, where N ¼ M m¼1 N m , PM K ¼ m¼1 K m . Hence, K additional equations are necessary, which arise from the collocation of the displacement continuity condition and the condition of equilibrium of tractions (Eq. (17c)) [50] on the interface nodes as

b int Þ ¼ ð U b int Þ ðU n  m    int b b int Gm U ;n ¼ Gn U ;n m

n

ð43aÞ ð43bÞ

225

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

(a)

(b)

Fig. 12. Distribution of rxx (a) and rVM (b) along the height of the cross section of the clamped beam of example 1 at X = 0.034 m.

Z

Z

pz=-0.1kN/cm

Pz=-4.7kN X

X l=47.0cm

l=47.0cm

(a)

(b)

b=7.2 cm

h=6.0 cm

(c) Fig. 13. Load case I (a), Load case II (b) and cross section (c) of the beam of example 2.

Thus employing Eq. (43), the total system of equations can be written as

( Htot

b int U b int U ;n

) ¼ f tot

ð44Þ

whereh Htot is a (Ni+T K)  (N + K) known coefficient matrix and b fr 0 , with 0 being a K  1 zero vector. The soluf tot ¼ f  G1 U ;n tion of system (44) yields all the unknown boundary quantities. It is worth here noting that in order to avoid singular behavior during the inversion of the coefficient matrix Htot, due to the fact that a Neumann type boundary value problem is solved, the regularization

technique presented in [52] is employed. Subsequently the warping function UPCY and its derivatives can be evaluated on any arbitrary point of the cross sectional domain employing the integral representation (36) and its derivatives and the already developed boundary element mesh. In order to avoid near singular behavior when an internal point approaches the boundary, the element subdivision technique is employed [50]. In case of monosymmetric cross sections, the coordinate ZS of the shear center S of the cross section, lying on the symmetry axis, can be obtained employing the procedure given in [3]. In order to maintain the pure boundary character of the formulation, all the geometric constants defined in Eq. (7) are evaluated

226

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

(ΦCYP )max = 2.7342E + 00 cm

(ϕCYP )max = 1.5983470E + 00 cm

(a)

(b)

Fig. 14. Distribution of UPCY (a) and uPCY (b) over the cross sectional domain of the beam of example 2.

Table 4 Geometric constants of the beam of example 2.

SuCz ¼

A = 2.4480E+01 cm2

jy = 1.6446E01

APZ ¼ 4:0260E þ 00 cm2

IYY = 1.4668E+02 cm4

ASZ ¼ 2:0454E þ 01 cm2

IuP uP ¼ 1:0220E þ 01 cm4 CY CY

Iyy ¼ Izz ¼ Iyz ¼

I uP

"



1 zny þ ynz z2 yds 2 CXm m¼1 # Z Z 1 1 þ ðyny þ znz ÞðuCz Þm ds  ðy2 þ z2 ÞðuCz;n Þm ds ð45dÞ 2 C Xm 4 CXm M X

km 

1 4Izz

Z M X km m¼1

CXm

m¼1

CXm

Z M X km M 1X km 4 m¼1

uP CY CY

¼

M X

Z

ðz2 yny Þds

ð45eÞ



 y2 znz ds

ð45fÞ

ðy2 zny þ z2 ynz Þds

ð45gÞ

Z

CXm

2

1 APZ km 4 I 840 YY m¼1

!2 Z C Xm

Z 7 nZ ds þ

1 APZ 120 IYY

Z

Z5

C Xm



uPCY;n

 m

ds

Z Z   1 APZ 1 APZ Z 4 nZ uPCY m ds  ðFÞm Z 2 nZ ds 24 IYY CXm 2 IYY CXm ! # Z Z   1 APZ 3  P   ðFÞm uPCY;n ds þ Z þ uCY m ðF ;n Þm ds ð45hÞ m C Xm CXm 6 IYY



Fig. 15. Deflection w(X) of the cantilever beam of example 2.

DUP

UP CY CY

by converting the domain integrals into line ones along the boundary using integration by parts, the Gauss theorem and the Green identity as

A¼

M 1X km 2 m¼1

Sy ¼

Z

M X

1 km 2 m¼1 M X

1 Sz ¼ km 2 m¼1

ðyny þ znz Þds

ð45aÞ

CXm

Z Z

z2 nz ds

ð45bÞ

CXm

CXm

( Z

1 1 2 1 2 3 Y þ Z Z nZ ds 8 CXm 3 5 m¼1 Z   1 2  ðY þ Z 2 ÞZ UPCY;n ds m 8 CXm )   Z  P  1 1 2 2 ZYnY þ ðY þ 3Z ÞnZ UCY m ds þ 8 CXm 4 ¼

M X

km

and using constant boundary elements for the approximation of these line integrals. In Eq. (45), ny, nz are the direction cosines, while F(y, z) is an auxiliary function defined as the solution of the following Neumann problem

r2 F ¼ uPCY y2 ny ds

ð45cÞ

ð45iÞ

Gm ðF ;n Þm ¼ 0

ð46aÞ on the free surface of the beam

Gm ðF ;n Þm ¼ Gn ðF ;n Þn

on the interfaces

ð46bÞ ð46cÞ

227

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234 Table 5  PA ðXÞ ðcmÞ; u  SA ðXÞ ðcmÞ and u  A ðXÞ ðcmÞ of the beam of example 2. Extreme values of the kinematical components wðXÞ ðcmÞ; hY ðXÞ ðradÞ; u

Load case I (w)max (h )  YP max  A max u 

  SA max u  A Þmax ðu

Load case II (w)max (h )  YP max  u  AS max  A max u  A Þmax ðu

Present study

FEM Euler–Bernoulli beam model [55]

FEM Timshenko beam model [55]

FEM solid model [55]

1.6934E+00 4.1172E02 9.3608E02

1.46184E+00 4.1472E02 9.4290E02

1.7238E+00 4.1472E02 9.4290E02

1.6467E+00 – –

9.3361E03

–

–

–

9.5330E02

9.4290E02

9.4290E02

9.5811E02

4.4083E+00 1.2437E01 2.8277E01

3.8979E+00 1.2440E01 2.8283E01

4.4219E+00 1.2440E01 2.8283E01

4.2990E+00 – –

1.2286E02

–

–

–

2.9505E01

2.8283E01

2.8283E01

2.9900E01

Upper surface 3 2.7

10

Lower surface

15 10

2.4 5

2.1 1.8

0

1.5 1.2

-5

5 0 -5

0.9 -10

0.6 0.3 0

5

10

15

20

25

(σVM )max = 2.947 kN

30

35

40

45

0

cm 2 (pres. study)

(σVM )max = 2.801 kN

cm 2 (FEM)

-10 -15

0

5

10

15

(σVM )max

20

25

30

35

40

45

4.25 4 3.75 3.5 3.25 3 2.75 2.5 2.25 2 1.75 1.5 1.25 1 0.75 0.5 0.25 0

2

= 4.394 kN cm (pres. study)

(σVM )max = 4.792 kN

cm 2 (FEM)

Fig. 16. Contour lines of rVM (kN/cm2) over the upper and lower surface of the cantilever beam of example 2, employing the present study (62 line elements) and a FEM solution (7040 solid elements) for load case I.

The establishment of the auxiliary function F is accomplished using BEM [47,53], employing the numerical procedure presented above. By establishing DUP UP and employing Eq. (18), APZ ; ASZ ¼ A  APZ can CY CY be directly computed. 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 the efficiency, wherever possible the accuracy and the range

of applications of the developed method. The numerical results have been obtained employing up to 63 nodal points (longitudinal discretization) and up to 900 boundary elements (cross sectional discretization). However it is noted that in all of the treated examples satisfactory accuracy could be also achieved with coarser discretization. 4.1. Example 1 For comparison reasons, in the first example a homogenous beam of a narrow rectangular cross section having its ends

228

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

(a)

(b)

Fig. 17. Absolute values of rxx along the boundary of the cross section of the cantilever beam of example 2 at X = 0 (a) and at X = 30 cm (b) (values in parentheses have been obtained from a FEM solution using 7040 solid elements).

and in Table 2 the geometric constants of the cross section are presented. In Fig. 8 the displacement w(X) of the beam is presented as compared with the ones obtained from FEM solution [54,55], employing Euler–Bernoulli beam elements (22 elements), Timoshenko beam elements (22 elements) and 8-noded hexahedral solid finite elements (4800 elements). From this figure the accuracy of the results obtained with the presented model especially comparing with the solid one (which is considered as the most accurate) is demonstrated, while the inability of the Euler–Bernoulli model to capture accurately the response of the beam can be observed. Furthermore, in order to demonstrate the influence of shear warping, in Fig. 9 the axial displacement components  PA ; u  SA ; u  A at the point A of the cross section (Fig. 6b) given by u the relations

 PA ðXÞ ¼ Z A hY ðxÞ u

ð47aÞ

 SA ðXÞ ¼ gY ðXÞuPCY ðY A ; Z A Þ u

ð47bÞ

 ðX; Y A ; Z A Þ ¼ u PA ðXÞ þ u SA ðXÞ  A ðXÞ  u u

ð47cÞ

for YA = 0 and ZA = 0.5 m are presented. From this figure the discrepancies arising employing elementary beam theories and the accuracy of the presented method can be verified, while in Table 3 the extreme values of the above kinematical components are presented. Moreover, in order to examine the arising stress distribution, in Figs. 10 and 11 the shear stress component sxz and the Von Mises equivalent stress given as Fig. 18. Distribution of shear lag factor k along the length of the cantilever beam of example 2 for load case I (a) and load case II (b).

clamped, as this is shown in Fig. 6, (E  Eref = 2.0  105 kN/m2, G  Gref = 8.333  104 kN/m2, l = 3.0 m) subjected to a uniformly distributed load pz = 0.1 kN/m, is examined. Due to the fact that the examined beam is a ‘‘short’’ one, significant influence of shear deformation is anticipated. In Fig. 7 the warping functions UPCY , uPCY

rVM ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2xx þ 3s2xy þ 3s2xz

ð48Þ

respectively, over the mid-plane of the beam as compared to the aforementioned solid FEM and Timoshenko beam model results are presented. Finally, in Fig. 12 the rxx and rVM profile along the midline of the cross section at X = 0.034 m are shown. From these figures, it can be easily verified that Timoshenko beam theory leads to significant discrepancies in the estimation of stress distributions, while the accuracy of the presented method as compared to the

229

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

Fig. 19. Distribution of shear lag factor k along the upper plate at X = 30 cm of the cantilever beam of example 2 for load case I.

Z

pz=-50.0 kN/m

(a) X l=30.0 m

A

%z b

0.30

E2 , G2

A

0.50

0.25

0.70

0.12

2.34

2.00

E1, G1

%zC zS

C C S S

O

(b)

h

zS %zC

B 1.11

0.36 0.08

%y

B

Fig. 20. Side view (a) and cross section (b) of the beam of example 3.

(ϕCYP )max = 5.2884E − 01 m

(ϕCYP )max = 6.0763E − 01 m

(a)

(b)

Fig. 21. Distribution of uPCY of the trapezoidal (a) and the rectangular (b) configuration of the composite cross section of example 3.

solid FEM model results is satisfactory. Hence, it can be concluded that the developed formulation is advantageous over more refined approaches, such as solid models, since it reduces significantly

computational cost (e.g. 22 line elements require approximately 0.07 s of analysis, while 4800 solid elements require approximately 6 s), it eliminates modeling effort, permits isolation of structural

230

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

phenomena and results interpretation and facilitates parametric analyses without compromising accuracy. Furthermore, the superiority of the presented beam formulation over Euler–Bernoulli or Timoshenko elements, which exhibit the same computational effort, is also verified by yielding results closer to solid FEM model. 4.2. Example 2

4.3. Example 3

For comparison reasons, a homogeneous cantilever beam of a monosymmetric box-shaped cross section, as this is shown in Fig. 13, (E  Eref = 284.49 kN/cm2, G  Gref = 101.6036 kN/cm2, l = 47.0 cm) under two load cases, is examined. More specifically, the beam is subjected to a uniformly distributed load pz = 0.1 kN/cm (case I, Fig. 13a) or to a concentrated force Pz = pz  l = 4.7 kN applied at its tip cross section (case II, Fig. 13b). In Fig. 14 the warping functions UPCY ; uPCY and in Table 4 the geometric constants of the cross section are presented. In Fig. 15 the displacement w(X) of the beam for load cases I, II is presented as compared with the ones obtained from FEM solutions [54,55], employing Euler–Bernoulli beam elements (62 elements), Timoshenko beam elements (62 elements) and 8-noded hexahedral solid finite elements (7040 elements). From this figure the significant influence of shear deformation in the deflection magnitude can be observed, while in Table 5 the extreme values of the kine Pj ðXÞ; u  Sj ðXÞ and u  j ðXÞ (Eq. (47)) matical components wðXÞ; hY ðXÞ; u at point j = A(yA = 7.5 cm, zA = 2.3412 cm) (Fig. 13c) for load cases I, II are given. Moreover, in Fig. 16 the distribution of Von Mises stress rVM (Eq. (48)) on the upper and lower surface of the beam for load case I, is presented as compared with the one obtained from the aforementioned solid FEM solution [55], and in Fig. 17 the j rxxj distribution along the boundary curve of the cross sections at X = 0 and X = 30 cm for load case I are presented. From these figures the significant alteration in normal stress distribution due to shear lag phenomenon is apparent, while the occurrence of the so-called ‘‘negative shear lag’’ phenomenon [18,24,36] can be observed. From these figures and table the accuracy of the present study is once again verified. For the treated example, results available in the literature could be also retrieved. In order to compare the present study with the beam formulations presented in [24,36], the ‘‘shear lag factor’’ [36] is defined as

 xxj k ¼ rxxj =r

ð49Þ

where rxxj is the normal stress value at a characteristic point j of the  xxj is the corresponding stress value at the same cross section and r point according to the engineering beam theory, given as

r xxj ¼ ðMy =Iyy Þzj

X = 30 cm is presented. From this figure, the influence of the flange thickness on the distribution of shear lag factor can be observed. From all of the above figures, the convergence of the present study results to the numerical and experimental ones can be once again verified.

In order to demonstrate the range of applications of the developed method, in this final example a composite beam (l = 30 m, E1 = 2.1  108 kN/m2, G1 = 8.0769  107 kN/m2, E2  Eref = 3.2  107 kN/m2, G2  Gref = 1.3333  107 kN/m2) of a monosymmetric box-shaped cross section having its left end clamped and its right one simply supported, as this is shown in Fig. 20 is examined. The beam is subjected to a uniformly distributed load pz = 50 kN/m (Fig. 20a). In this study, three values of width b of the concrete plate are examined (Fig. 20b). Furthermore, in order to examine the influence of nonuniform shear warping in different cross sectional shapes, two configurations of the steel box, a trapezoidal (configuration I) and a rectangular one (configuration II) (Fig. 20b) are examined, while the shape of the concrete plate

Table 6 Geometric constants of the beams of example 3. b = 10 m

b = 12 m

b = 14 m

Configuration I (trapezoidal box) A (m2) 8.4452E+00 1.923225396 AP ðm2 Þ

1.0095E+01 1.7034E+00

1.1745E+01 1.4808E+00

6.521974604

8.3917E+00

1.0264E+01

7.4918E+00 1.5879E01

1.0426E+01 3.9752E01

1.3171E+01 9.2851E01

1.6226E+00 5.5625E01 2.2773E01

1.5084E+00 7.1815E01 1.6874E01

1.4262E+00 7.7324E01 1.2608E01

Configuration II (rectangular box) h (m) 1.7069E+00 A (m2) 8.9812E+00 1.5994E+00 AP ðm2 Þ

1.7990E+00 1.0728E+01 1.4605E+00

1.8497E+00 1.2431E+01 1.3017E+00

7.3818 E+00

9.2675E+00

1.1129E+01

7.4918E+00 3.5915E01

1.0426E+01 7.5591E01

1.3171E+01 1.3641E+00

1.5703E+00 4.3976E01 1.7808E01

1.4581E+00 5.4207E01 1.3614E01

1.3803E+00 5.7886E01 1.0471E01

Z ASZ

ðm2 Þ IYY (m4) IuP uP ðm4 Þ CY

CY

~zC ðmÞ ZS (m)

jZ

Z

ASZ ðm2 Þ IYY (m4) IuP uP ðm4 Þ CY

CY

~zC ðmÞ ZS (m)

jZ

ð50Þ

where zj is the z coordinate of point j. In the present study, for comparison reasons, points A, B (Fig. 13c) located at the joint between the upper plate and the web are selected. It is worth noting that point A lies on the upper boundary curve of the cross section while B refers to the mid-line of the upper plate. In Fig. 18 the distribution of shear lag factor k along the beam length for load cases I, II is presented as compared with FEM solutions [24,36], variational approach [24] and experimental results [24]. It is worth here noting that the numerical results presented in [24,36] have been obtained employing thin-walled assumptions; hence, the variation of warping along the thickness of the members (webs, flanges) of the cross section cannot be taken into account. The substantially different behavior of shear lag factor between both load cases can be observed. More specifically, in load case I the uniformly distributed load leads to the generation of negative shear lag (region where k < 1) while, concentrated force does not have the same effect. Furthermore, in Fig. 19 the distribution of shear lag factor along the half of the upper plate of the cross section for load case I at

Fig. 22. Deflection w(X) of the beam of example 3 for b = 12 m.

231

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

remains exactly the same. It is noted here that the height h of the web of the rectangular box is chosen so that both of the cross sections have equal second moments of inertia IYY. In Fig. 21 the warping functions uPCY for both configurations I, II for b = 12 m and in Table 6 the geometric constants of all of the cross sections are given. In Fig. 22 the displacement w(X) of the beam for b = 12 m is presented as compared with the ones obtained from the analytical integration [56] of the Euler–Bernoulli and Timoshenko differential equilibrium equations given in [57]. For the solution of the aforementioned differential equations the values of A, IYY and jZ presented in Table 6 were employed. In this figure, the significant influence of shear deformation in the deflection magnitude can be observed. Furthermore, it is noted that the Euler–Bernoulli solution is independent of the cross sectional configuration since IYY is the same for both cases; however when shear deformation is involved, each configuration yields different deflection curves. Moreover, the discrepancies of Timoshenko beam theory from the present one are noticed. Furthermore, in Fig. 23 the distribution of jrxxj along the boundary curve of the cross sections and the shear stress s vectors on the cross sectional domain, at X = 0.484 m for configurations I, II and b = 12 m are presented. In this figure the significant alteration in normal stress distribution due to shear lag phenomenon can be once again observed. Finally, in Fig. 24 the shear force t(X) generated on the interface between each steel flange and the concrete plate of configuration I for b = 12 m is presented as compared to the one calculated by the approximate formula often employed for the design of shear connectors in composite structures

Fig. 24. Shear interface force t(X) between steel flange and concrete plate along the beam of example 3 for configuration I and b = 12 m.

tðXÞ ¼

Scut Y Q z ðXÞ IYY

ð51Þ

where Scut Y is the first moment of inertia of the concrete plate. In this figure, it can be observed that in a region near the fixed support, the

(a)

σ xx max = 6658.7259 kN m 2

σ xx max = 6458.7980 kN m 2 1

1

0

0

(b) -1

-1

-6

-5

50

-4

-3

-2

-1

150 250 350 450 550 650 750 850 950 1050 1150 1250

(τ )max = 1263.4000 kN

m2

0

1 100

200

2 300

400

3 500

600

4 700

800

(τ )max = 1118.2000 kN

5

6

900 1000 1100

m2

Fig. 23. Absolute values of rxx (kN/m2) along the boundary of the cross section (a) and shear stress vectors (b) of the composite beam of example 3 at X = 0.484 m for b = 12 m.

232

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

Table 7  PA ðXÞ ðmÞ; u  SA ðXÞ ðmÞ; u  A ðXÞ ðmÞ; u  PB ðXÞ ðmÞ; u  SB ðXÞ ðmÞ and u  B ðXÞ ðmÞ of the beam of example 3. Extreme values of the kinematical components wðXÞ ðmÞ; hY ðXÞ ðradÞ; u Configuration I

Configuration II

b = 10 m

b = 12 m

b = 14 m

b = 10 m

b = 12 m

b = 14 m

(w)max (hY)max (hY)min

1.1322E03 7.5787E05 1.2536E04

8.9924E04 5.31572E05 9.3269E05

7.8739E04 4.10109E05 7.6300E05

1.1675E03 7.51751E05 1.2646E04

9.2771E04 5.27389E05 9.4012E05

8.1091E04 4.07134E05 7.6817E05

Axial displacement  P  u  AP max  A min u  S  u  AS max  A min u  A Þmax ðu  A Þmin ðu

at point A 7.5588E05

5.9090E05

4.8957E05

7.8915E05

6.1278E05

5.0471E05

1.2503E04

1.0368E04

9.1083E05

1.3275E04

1.0923E04

9.5227E05

Axial displacement  P  u  BP max  B min u  S  u  BS max  B min u  B Þmax ðu  B Þmin ðu

at point B 2.0342E04

1.4069E04

1.0882E04

1.6151E04

1.1818E04

1.2298E04

8.0182E05

5.8492E05

9.6010E05

6.6297E05

1.2269E05

1.4477E05

1.7420E05

1.3772E05

1.5366E05

1.7971E05

8.1602E06

1.0065E05

1.2796E05

9.5109E06

1.1111E05

1.3647E05

8.4203E05 1.3319E04

7.0163E05 1.1374E04

6.3687E05 1.0388E04

8.9285E05 1.4226E04

7.3888E05 1.2035E04

6.6415E05 1.0887E04 9.44887E05 5.0079E05

4.1472E06

7.9757E06

1.2453E05

9.9473E06

1.4643E05

1.9087E05

6.2356E06

1.1472E05

1.6952E05

1.4404E05

2.0251E05

2.5134E05

2.0757E04 1.2734E04

1.48662E04 8.8853E05

1.2128E04 7.2730E05

1.7146E04 1.0680E04

1.3282E04 8.2955E05

1.1358E04 7.2601E05

above formula is inaccurate leading to significant underestimation of the interface force. It is worth noting that force t(X) corresponds to the case of rigid shear connection since in the present study materials are considered firmly bonded together. Taking into account the deformability of this connection would lead to reduced values of interface force t(X) and larger values of deflection w(x) as also observed in the study of Sapountzakis and Mokos [58] which deals with the deformability of shear connection in plates stiffened by parallel beams. It is also noted that it is possible to extend the developed beam formulation so as to account for this effect as well. Furthermore, in Table 7 the extreme values of the kinematical  Pj ðXÞ; u  Sj ðXÞ and u  j ðXÞ (Eq. (47)) and in components wðXÞ; hY ðXÞ; u Table 8 the values of stress components rxxj, rVMj and rPxxj ; rSxxj corresponding to the primary and secondary normal stress distribution (Eq. (3a)), respectively, at j = A, B points (Fig. 20b)

and X = 0, X = 18.871 m, are given. In Fig. 25 the distribution of secondary normal stress component rSxxA at the point A of the cross section (Fig. 20b) along the beam length for configurations I, II and various values of b is presented. In Fig. 26 the distribution of  SA ðXÞ along the beam length the axial displacement component u for configurations I, II and various values of b is also shown. From these figures and tables, the increment of warping displacements at the point A, located above the joint of the steel box and the concrete plate, as b increases can be verified; it can also be observed that configuration II leads to larger values of warping displacement at the joint. As far as the stress magnitudes are concerned, it can be verified that in the region of positive bending moment, normal stress rSxx due to warping increases with the increment of b. It can also be observed that configuration II leads to increased rSxx as compared to configuration I. In the region of negative bending moment the opposite behavior occurs.

Table 8 Values of stress components rxx (kN/m2) and rVM (kN/m2) at the points A, B of the cross section of the beam of example 3. Configuration I b = 10 m

Configuration II b = 12 m

b = 14 m

b = 10 m

b = 12 m

b = 14 m

X = 0, stress at point A 7.1001E+02 rPxx 8.4215E+02 rSxx

5.6240E+02 6.5757E+02

4.7297E+02 5.3714E+02

7.4544E+02 6.3201E+02

5.8660E+02 5.1236E+02

4.8976E+02 4.5763E+02

rxx rVM

1.5522E+03 1.5612E+03

1.2200E+03 1.2286E+03

1.0101E+03 1.0182E+03

1.3774E+03 1.3871E+03

1.0990E+03 1.1078E+03

9.4739E+02 9.5527E+02

X = 0, stress at point B 7.5807E+03 2.8088E+03

5.0082E+03 3.4196E+03

3.7084E+03 3.4302E+03

5.9516E+03 4.3379E+03

4.1649E+03 4.4273E+03

3.1891E+03 4.2003E+03

rxx rVM

1.0389E+04 1.0439E+04

8.4278E+03 8.4735E+03

7.1386E+03 7.1807E+03

1.0290E+04 1.0337E+04

8.5922E+03 8.6339E+03

7.3894E+03 7.4266E+03

X = 18.871 m, stress at point A 4.2426E+02 2.4110E+01

3.4281E+02 3.0464E+01

2.9376E+02 4.0354E+01

4.4715E+02 2.8851E+01

3.5877E+02 3.6261E+01

3.0537E+02 4.44661E+01

rxx rVM

4.4837E+02 4.4837E+02

3.7327E+02 3.7327E+02

3.3412E+02 3.3412E+02

4.7600E+02 4.7600E+02

3.9503E+02 3.9503E+02

3.4983E+02 3.4983E+02

X = 18.871 m, stress at point B 4.5298E+03 8.0412E+01

3.0527E+03 1.5842E+02

2.3033E+03 2.5770E+02

3.5701E+03 1.9802E+02

2.5472E+03 3.1132E+02

1.9884E+03 4.0813E+02

rxx rVM

3.2111E+03 3.2113E+03

2.5610E+03 2.5614E+03

3.7681E+03 3.7681E+03

2.8586E+03 2.8586E+03

2.3966E+03 2.3968E+03

rPxx rSxx

rPxx rSxx

rPxx rSxx

4.6102E+03 4.6103E+03

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

Fig. 25. Secondary normal stress component rSxxA ðXÞ at point A of the cross section of the beam of example 3 for configurations I, II and various values of width b.

233

 It permits isolation of structural phenomena and results interpretation (quantities such as rotation, warping parameter, stress resultants, etc. are also evaluated in contrast to solid model which yields only displacements and stress components).  It allows straightforward model handling (boundary conditions and external loading are easily applied).  It facilitates parametric analyses (solid modeling often requires construction of multiple models). b. In the examined examples, the influence of nonuniform shear warping on the distribution of normal stress was apparent, while different type of loading led to a significantly different effect on box-shaped cross sections known as ‘‘negative shear lag’’. c. In the examined box-shaped composite cross sections, the increment of flange width resulted in the increment of axial warping displacements at the joint between plate and web, while the shape of the box also altered their magnitude. More specifically, in the case of rectangular configuration larger values of axial warping displacements at the joint were generated, as compared with the ones computed in the case of trapezoidal one. d. As far as stress distribution in the examined box-shaped composite beams is concerned, in the region of positive moments, the increment of plate width resulted in larger normal stress values at the joint, while in the case of rectangular configuration, further increase was observed. The opposite behavior occurred in the region near the fixed support.

Acknowledgements This research has been co-financed by the European Union (European Social Fund – ESF) and Greek National Funds through the Operational Program ‘‘Education and Lifelong Learning’’ of the National Strategic Reference Framework (NSRF) – Research Funding Program: THALES: Reinforcement of the interdisciplinary and/or inter-institutional research and innovation. References SA ðXÞ at point A of the cross section of the Fig. 26. Axial displacement component u beam of example 3 for configurations I, II and various values of width b.

5. Concluding remarks In this paper, a boundary element solution is developed for the nonuniform shear warping analysis of composite beams of arbitrary monosymmetric cross section. The beam is subjected to arbitrarily distributed or concentrated transverse loads, passing through the shear center of the cross section as well as bending and warping moments, while its ends are supported by the most general boundary conditions. The main conclusions that can be drawn from this investigation are: a. The developed numerical procedure retains most of the advantages of a BEM solution over a FEM approach. The proposed beam formulation is capable of yielding results of high accuracy, as verified by comparing with solid FEM models and experimental results, with minimum computational cost. Its advantageous character over more refined approaches is also enhanced by the following:  The developed beam formulation reduces significantly modeling effort (solid models require cumbersome preprocessing even in relatively simple cases).

[1] Schramm U, Kitis L, Kang W, Pilkey WD. On the shear deformation coefficient in beam theory. Finite Elem Anal Des 1994;16:141–62. [2] Nouri T, Gay D. Shear stress in orthotropic composite beams. Int J Eng Sci 1994;32(10):1647–67. [3] Sapountzakis EJ, Mokos VG. A displacement solution to the transverse shear loading of composite beams by BEM. Comput, Mater, Continua 2009;10(1):1–39. [4] Razaqpur AG, Li HG. A finite element with exact shape functions for shear lag analysis in multi-cell box girders. Comput Struct 1991;39(1):155–63. [5] Reissner E. Analysis of shear lag in box beams by the principle of minimum potential energy. Q Appl Math 1946;41:268–78. [6] Malcolm DJ, Redwood RG. Shear lag in stiffened box-girders. J Struct Div ASCE 1970;96(ST7):1403–15. [7] Moffatt KR, Dowling PJ. Shear lag in steel box-girder bridges. Struct Eng 1975;53:439–48. [8] Murray NW. Introduction to the theory of thin-walled structures. Oxford University Press; 1986. [9] Eurocode 3: design of steel structures – part 1.5: plated structural elements. European Committee for Standardization, prEN 1993-1-5; 2004. [10] Eurocode 4: design of composite steel and concrete structures – part 1.1: general rules and rules for buildings. European Committee for Standardization, prEN 1994-1-1; 2004. [11] Eurocode 4: design of composite steel and concrete structures – part 2: general rules and rules for bridges. European Committee for Standardization, prEN 1994-2; 2004. [12] Ie CA, Kosmatka JB. On the analysis of prismatic beams using first-order warping functions. Int J Solids Struct 1992;29(7):879–91. [13] Kristek V, Evans HR, Ahmad MKM. A shear lag analysis for composite box girders. J Construct Steel Res 1990;16:1–21. [14] Lopez-Anido R, GangaRao HVS. Warping solution for shear lag in thin-walled orthotropic composite beams. J Eng Mech 1996;122(5):449–57.

234

I.C. Dikaros, E.J. Sapountzakis / Engineering Structures 76 (2014) 215–234

[15] Tahan N, Pavlovic´ MN, Kotsovos MD. Shear-lag revisited: the use of single fourier series for determining the effective breadth in plated structures. Comput Struct 1997;63(4):759–167. [16] Adekola AO. On shear lag effects in orthotropic composite beams. Int J Solids Struct 1974;10:735–54. [17] Katsikadelis JT, Sapountzakis EJ. A realistic estimation of the effective breadth of ribbed plates. Int J Solids Struct 2002;39:897–910. [18] Lee SC, Yoo CH, Yoon DY. Analysis of shear lag anomaly in box girders. J Struct Eng 2002;128(11):1379–86. [19] Sa-nguanmanasak J, Chaisomphob T, Yamaguchi E. Stress concentration due to shear lag in continuous box girders. Eng Struct 2007;29:1414–21. [20] Gupta PK, Singh KK, Mishra A. Parametric study on behaviour of box-girder bridges using finite element method, technical note. Asian J Civil Eng (Build Housing) 2010;11(1):135–48. [21] El Fatmi R, Ghazouani N. Higher order composite beam theory built on SaintVenant’s solution. Part-I: Theoretical developments. Compos Struct 2011;93:557–66. [22] Vlasov V. Thin-walled elastic beams. Israel Program for Scientific Translations, Jerusalem; 1963. [23] Dezi L, Mentrasti L. Nonuniform bending-stress distribution (shear lag). J Struct Eng 1985;111(12):2675–90. [24] Chang ST, Zheng FZ. Negative shear lag in cantilever box girder with constant depth. J Struct Eng 1987;113(1):20–35. [25] Luo QZ, Li QS. Shear lag of thin-walled curved box girder bridges. J Eng Mech 2000;126(10):1111–4. [26] Wu Y, Zhu Y, Lai Y, Pan W. Analysis of shear lag and shear deformation effects in laminated composite box beams under bending loads. Compos Struct 2002;55:147–56. [27] Luo QZ, Tang J, Li QS. Shear lag analysis of beam-columns. Eng Struct 2003;25:1131–8. [28] Wu Y, Liu S, Zhu Y, Lai Y. Matrix analysis of shear lag and shear deformation effects in thin-walled box beams. J Eng Mech 2003;129(8):944–50. [29] Wu Y, Lai Y, Zhang X, Zhu Y. A finite beam element for analyzing shear lag and shear deformation effects in composite-laminated box girders. Comput Struct 2004;82:763–71. [30] Zhou SJ. Finite beam element considering shear-lag effect in box girder. J Eng Mech 2010:1115–22. [31] Gara F, Ranzi G, Leoni G. Simplified method of analysis accounting for shearlag effects in composite bridge decks. J Construct Steel Res 2011;67:1684–97. [32] Hjelmstad KD. Warping effects in transverse bending of thin-walled beams. J Eng Mech 1987;113(6):907–24. [33] Koo KK, Cheung YK. Mixed variational formulation for thin-walled beams with shear lag. J Eng Mech 1989;15(10):2271–86. [34] Laudiero F, Savoia M. Shear strain effects in flexure and torsion of thin-wailed beams with open or closed cross-section. Thin-Wall Struct 1990;10:87–119. [35] Koo KK, Wu XS. Shear lag analysis for thin-walled members by displacement method. Thin-Wall Struct 1992;13:337–54. [36] Prokic´ A. A new finite element for analysis of shear lag. Comput Struct 2002;80:1011–24. [37] Massonnet CE. A new approach (including shear lag) to elementary mechanics of materials. Int J Solids Struct 1983;19(1):33–54.

[38] Hofmann TJ. Beitrag zur verfeinerten Balkentheorie. Dissertation, Bericht 15, Institut für Baustatik, Universität Stuttgart; 1992. [39] Park SW, Fuji D, Fujitani Y. A finite element analysis of discontinuous thinwalled beams considering nonuniform shear warping deformation. Comput Struct 1997;65(1):17–27. [40] El Fatmi R. Non-uniform warping including the effects of torsion and shear forces. Part-I: A general beam theory. Int J Solids Struct 2007;44:5912–29. [41] El Fatmi R. Non-uniform warping including the effects of torsion and shear forces. Part-II: Analytical and numerical applications. Int J Solids Struct 2007;44:5930–52. [42] Ghazouani N, El Fatmi R. Extension of the non-uniform warping theory to an orthotropic composite beam. C R Mec 2010;338:704–11. [43] Ghazouani N, El Fatmi R. Higher order composite beam theory built on SaintVenant’s solution. Part-II: Built-in effects influence on the behavior of endloaded cantilever beams. Compos Struct 2011;93:567–81. [44] Le Corvec V, Filippou FC. Enhanced 3D fiber beam-column element with warping displacements. In: Proc. of the 3rd international conference on computational methods in structural dynamics and earthquake engineering COMPDYN; 2011. [45] Ferradi MK, Cespedes X, Arquie M. A higher order beam finite element with warping eigenmodes. Eng Struct 2013;46:748–62. [46] Katsikadelis JT. The analog equation method. A boundary-only integral equation method for nonlinear static and dynamic problems in general bodies. Theor Appl Mech 2002;27:13–38. [47] Sapountzakis EJ, Mokos VG. Warping shear stresses in nonuniform torsion of composite bars by BEM. Comput Methods Appl Mech Eng 2003;192:4337–53. [48] Muskhelishvili NI. Some basic problems of the mathematical theory of elasticity. P. Noordhoff Ltd; 1963. [49] Sapountzakis EJ, Katsikadelis JT. Analysis of plates reinforced with beams. Comput Mech 2000;26:66–74. [50] Katsikadelis JT. Boundary elements: theory and applications. AmsterdamLondon: Elsevier; 2002. [51] Beer G, Smith I, Duenser Ch. The boundary element method with programming – for engineers and scentists. Wien, New York: Springer; 2008. [52] Lutz E, Ye W, Mukherjee S. Elimination of rigid body modes from discretized boundary integral equations. lnt J Solids Struct 1998;35(33):4427–36. [53] Sapountzakis EJ, Mokos VG. Nonuniform torsion of composite bars by boundary element method. J Eng Mech, Am Soc Civil Eng (ASCE) 2001;127:945–53. [54] Siemens PLM Software Inc., NX Nastran User’s Guide; 2008. [55] FEMAP for Windows. Finite element modeling and post-processing software. Help System Index, Version 10; 2008. [56] Maplesoft, a division of Waterloo Maple Inc., Maple User Manual, Maple 12.0; 2008. [57] Armenakas AE. Advanced mechanics of materials and applied elasticity. New York: Taylor & Francis Group; 2006. [58] Sapountzakis EJ, Mokos VG. An improved model for the analysis of plates stiffened by parallel beams with deformable connection. Comput Struct 2008;86:2166–81.