- Email: [email protected]
The line element-less method analysis of orthotropic beam for the De Saint Venant torsion problem
ARTICLE IN PRESS International Journal of Mechanical Sciences 52 (2010) 43–55
Contents lists available at ScienceDirect
International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci
The line element-less method analysis of orthotropic beam for the De Saint Venant torsion problem Roberta Santoro Dipartimento di Ingegneria Strutturale, Aerospaziale e Geotecnica DISAG, Universita degli Studi di Palermo, Viale delle Scienze, 90128 Palermo, Italy
a r t i c l e in fo
abstract
Article history: Received 3 February 2009 Received in revised form 22 September 2009 Accepted 7 October 2009 Available online 13 October 2009
This paper deals with the extension of a novel numerical technique, labelled line element-less method (LEM), in order to provide approximate solutions of the De Saint Venant torsion problem for orthotropic beams having simply and multiply connected cross-section. A suitable transformation of coordinates allows to take full advantage of the theory of analytic complex functions as in the isotropic case. A complex potential function analytic in all the transformed domain whose real and imaginary parts are related to the shear stress components and to the orthotropic ratio is introduced and expanded in the double-ended Laurent series involving harmonic polynomials. An element-free weak-form procedure has been proposed imposing that the square of the net flux of the shear stress across the border is minimized with respect to the series coefficients. Numerical implementation of the LEM results in system of linear algebraic equations involving symmetric and positive-definite matrices. All the integrals are transferred into the boundary without requiring any discretization neither in the domain nor in the contour. The technique provides the complete shear stress field as shown by some numerical applications in order to assess the efficiency and the accuracy of the method to handle shear stress problems in presence of orthotropic material. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Orthotropic material Analytic functions Torsion problem
1. Introduction Several materials used in construction show a sharp difference in elastic properties for different directions. Examples of such materials exhibiting elastic anisotropic behaviour are provided by natural wood, composites and synthetic materials used in aircraft construction. To design anisotropic members undergoing elastic strains, it is necessary to know how to determine the stresses and strains in an anisotropic body theoretically, i.e. to solve torsion problems of the theory of elasticity for anisotropic body as presented in this paper. Benchmark solutions for torsion and flexure problems in the analysis of orthotropic beams are represented by the series expansion solutions due to Lekhnitskii [1] and Tolf [2], respectively. The commonly used techniques to pursue numerical solution of the orthotropic De Saint Venant problem are based on the finite element method (FEM) and boundary element method (BEM). A finite element approach has been proposed in Herrmann [3] and Mason and Herrmann [4] for the two-dimensional torsion and flexure analysis of beams with arbitrary cross-sections restricted to the homogeneous isotropic case. This method was also
Tel.: + 39 091 6568406; fax: +39 091 6568407.
E-mail addresses: [email protected], [email protected]. 0020-7403/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmecsci.2009.10.002
extended to provide solutions for beams made by anisotropic material by Kosmatka and Dong [5]. A Ritz-based power series method was used for the solution of the flexure-torsion problem of a general anisotropic material with arbitrary cross-section as reported in Ref. [6]. Both the above methods require the whole cross-section to be discretized into area elements. Although the solving system obtained via FEM involves symmetric and positivedefinite matrices providing simple and immediate solutions, a large number of elements is required for a good accuracy; furthermore, the limitation with respect to the shape of the elements leads to cumbersome meshing process. In the Ritz-based methods the elements of the mesh are less limited with respect to their shape because they are created only to aid the numerical integration on the cross-sectional area but the integration requires a large number of integration points because of the high order of the power series. On the other hand for its features BEM allows a more accurate evaluation of the stress field over the cross-section. Advantage of the latter method is related with an easier discretization of the boundary resulting in a relatively small number of line elements with the possibility of addressing the computational resources to restricted regions of interest. The use of boundary integral method for the torsion problem for homogeneous and isotropic material has been proposed by Jawson and Ponter [7] by using a simple constant interpolation and by Chou and Mohr [8] and Katsikadelis and Sapountzakis [9] to handle composite materials. In Sapount-
ARTICLE IN PRESS 44
R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
zakis et al. [10–13] the BEM is employed to solve the non-uniform torsion problem of composite bars of arbitrary variable cross-sections. In Ref. [14] a more refined analysis is presented with a three-node isoparametric boundary element to allow an accurate geometric representation of cross-sections having curved boundaries. In the latter paper the stress field over the crosssection is described at points quite far from the boundary probably due to the difficulties implied by numerical integration. Solution of De Saint Venant torsion problem obtained via FEM and BEM has also been used to evaluate torsion factors [15]. Within the orthotropic analysis Gracia and Doblare [16] solved an optimization problem for sections under torsion. A boundary element model in the analysis of De Saint Venant orthotropic problem is presented in [17], where the differential equations governing the shear stress field are converted in three Neumann problems which form the basis of the boundary integral formulation and a quadratic B-spline approximation represents the boundary variables. The application of complex variables in the boundary element methods allows to obtain the solution for two conjugate functions at once. These related methods intend to solve boundary problems of Laplace equation and distinguish in complex polynomial method (CPM) [18,19] and complex variable boundary element method (CVBEM) [20,21]. Both the above methods [22] provide solution of two-dimensional potential problems as the torsion of elastic rods making use of complex variables by generalization of the Cauchy integral formula into a boundary integral equation method. The CVBEM was also proposed for non-uniform torsion with isotropic homogeneous material [23] and for pure torsion of orthotropic prismatic bar [24]. Recently Di Paola et al. [25] propose to solve the problem of pure torsion for the De Saint Venant cylinder by using a novel approach labelled line element-less method (LEM). As in the CVBEM and CPM the method firstly works in terms of a potential function. In order to avoid to calculate both warping function and its harmonic conjugate a new potential function involving directly the shear stresses has been introduced. This complex function, holomorphic in all the domain, allows to take full advantage of the double-ended Laurent series involving harmonic polynomials. The selected expression of this analytic complex function is such that governing equations of De Saint Venant elastic problem, namely equilibrium and compatibility equations, remain fulfilled in the cross-section domain. While in the BEM theory the shape functions are the fundamental ones, in the proposed method the shape functions are harmonic polynomials that satisfy the field equations in the whole domain. The coefficients of the Laurent series are then evaluated by imposing a weak-form condition that the square value of the total net flux across the border remains minimum with respect to parameters expansion under the condition that the static equivalence is satisfied. Numerical implementation of LEM results in systems of linear algebraic equations in positive-defined and symmetric matrices solving only contours integral. The method provides exact solutions when existing and very accurate results by using very few series terms in the other cases. In this paper the above explained method is extended to the case of orthotropic material by means of a suitable transformation of coordinates [26] with the introduction of a complex potential function analytic in all the domain whose real and imaginary parts are related to the shear stress components and to the orthotropic ratio. Some numerical applications to simply and multiply connected domains have been reported to demonstrate the accuracy and the efficiency of the proposed LEM to handle shear stress problem in presence of orthotropic material.
2. Introductory remarks By definition, a material having two or three mutually orthogonal planes of elastic symmetry at each point (or, what is the same thing, two or three mutually perpendicular principal directions) is said to be orthogonally anisotropic or, more shortly, orthotropic. Thus, its mechanical properties are, in general, different along the directions of each of the axes. One common example of an orthotropic material with two axes of symmetry would be a polymer reinforced by parallel glass or graphite fibres. The strength and stiffness of such a composite material will usually be greater in a direction parallel to the fibres than in the transverse direction. Such materials are sometimes called transverse isotropic. A familiar example of an orthotropic material with three mutually perpendicular axes is provided by natural wood, in which the properties (such as strength and stiffness) along its grain and in each of the two perpendicular directions are different. Orthotropic materials require 9 independent variables (i.e. elastic constants) in their constitutive matrices. In contrast, a material without any planes of symmetry is fully anisotropic and requires 21 elastic constants, whereas a material with an infinite number of symmetry planes (i.e. every plane is a plane of symmetry) is isotropic, and requires only 2 elastic constants. The 9 independent elastic constants in orthotropic constitutive equations are comprised of three Young’s moduli Ex, Ey, Ez, the three Poisson’s ratios nzx, nzy, nxy, and the three shear moduli Gzx, Gzy, Gxy. By choosing the axes normal to the three planes of symmetry, the compliance matrix takes the form ex 1=Ex nyx =Ey nzx =Ez 0 0 0 sx e 1=Ey nzy =Ez 0 0 0 sy y nxy =Ex ez nxz =Ex nyz =Ey 1=Ez 0 0 0 sz g ¼ 0 0 0 1=Gyz 0 0 tyz yz gxz 0 t 0 0 0 0 1=G xz xz 0 0 0 0 0 1=Gxy txy gxy ð1Þ It is worth to be
nij Ei
¼
nji
ð2Þ
Ej
where generally nij a nji. Assuming sx = sy = txy =0 as in the De Saint Venant problem the constitutive equations become:
ex ¼
nzx Ez
gxy ¼ 0;
sz ; gyz ¼
ey ¼ tyz Gyz
;
nzy Ez
sz ;
gxz ¼
ez ¼
1 sz Ez
txz
ð3aÞ
ð3bÞ
Gxz
The equilibrium equations are the same of the isotropic case @tzx ¼ 0; @z
@tzy ¼ 0; @z
@tzx @tzy @sz þ þ ¼0 @x @y @z
ð4a2cÞ
The compatibility equations become 2
@ ex @2 ey @2 ex @2 ez þ 2 ¼ 0; þ 2 ¼ 0; 2 @y @x @z2 @x @2 ex 1 @ @gzx @gyz @2 ey ¼ ; ¼ @y@z @x @x@z 2 @x @y @2 ez ¼0 @x@y
@2 ey @2 ez þ 2 ¼0 @z2 @y 1 @ @gyz @gxz ; @y 2 @y @x ð5Þ
ARTICLE IN PRESS R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
and bearing in mind Eqs. (3) and (4) can be rewritten in terms of stresses as follows:
mzx
@2 sz @2 sz þ mzy ¼ 0; 2 @y @x2
@2 sz ¼ 0; @x2
mzx @2 sz Ez
mzy Ez
@2 sz ¼ 0; @y2 !
@2 sz ¼0 @x@y
1 1 @2 tzx 1 @2 tyz ; ¼ @y@z 2 Gxz @x@y Gyz @x2 ! @2 sz 1 1 @2 tzy 1 @2 txz ¼ @x@z 2 Gyz @x@y Gxz @y2
ð6Þ
ðtzy x tzx yÞ dO ¼ Mz
Let us consider a beam of arbitrary cross-section having three mutually orthogonal planes of elastic symmetry twisted by a moment Mz. It is assumed that the axis of the beam is perpendicular to one of these planes. As in the isotropic case, we let the axis of the beam coincide with the z-axis and choose the x- and y-axis in one end of the beam oriented as the material directions of orthotropy (see Fig. 1). It is worth to note that in the most general case the two longitudinal planes of elastic symmetry, namely xz and yz, do not necessarily intersect the plane of the section at a given point along lines parallel to the principal axes of the cross-section [26,27]. The shear moduli associated with the z- and x-axis and z- and y-axis are denoted by Gzx and Gzy, respectively. We assume, as in the isotropic case, that the components of displacement in the directions of x, y and z axes denoted as u, v and w, respectively, are given by the formulas: w ¼ yoðx; yÞ
ð7a2cÞ
where y is the twist angle for unit length of the bar and o(x, y) is the warping function associated with this problem. It is easy to verify [26] that in the orthotropic case, the torsion function o no longer satisfies Laplace equation in the domain O of the crosssection. As in the isotropic case, the non-vanishing stress components should satisfy the free-traction boundary conditions on the contour of the cross-section, namely
tzx nx þ tzy ny ¼ sT n ¼ 0 on C
ð9Þ
has to be satisfied. A simplified compact formulation of the De Saint Venant orthotropic torsion problem can be obtained following Sokolnikoff [26] by the introduction of new independent variables x and Z defined as follows:
x ¼ Rx;
3. Formulation of the pure torsion problem for an orthotropic beam
v ¼ yzx;
Moreover the static equivalence Z O
Once equilibrium, compatibility and constitutive law are reported for the orthotropic De Saint Venant beam in the next section the problem of pure torsion will be presented.
u ¼ yzy;
45
ð8Þ
Z¼y
ð10a; bÞ
where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R ¼ Gzy =Gzz
ð11Þ
is defined as the orthotropic ratio. As it will be evident in the following the variable transformation in Eq. (10) allows to formulate the torsion orthotropic problem as a Neumann or Dirichlet problem in the transformed domain where the governing equations are of the same form as that appearing in the study of torsion of isotropic cylinder. The shear stress components in the transformed domain can be easily written in terms of warping function: @o @o ð12a; bÞ tzx ¼ yGzx R Z ; tzy ¼ yRGzx R þx @x @Z and the third equilibrium equation, namely div s ¼ 0, becomes: @2 oðx; ZÞ @x
2
þ
@2 oðx; ZÞ ¼ r2 oðx; ZÞ ¼ 0 @Z2
in O
ð13Þ
It follows that o(x, Z) is an harmonic function in the transformed domain O whose equation is changed in f ðx; ZÞ ¼ 0. Making the corresponding change of variables in the boundary condition (8) it returns:
tzx n x þ tzy n y ¼ sT n ¼ 0 on C
ð14Þ
or in terms of warping function ! @o @f @o @f @f @f R ¼Z x þ @x @x @Z @Z @x @Z
on C
ð15Þ
where plainly C is the transformed boundary. It is worth to note that the components of the normal n to the boundary C are proportional to R@f =@xand @f =@Z, respectively. Hence the solution of the torsion problem for a cylinder of no-isotropic material whose cross-section has contour C is reduced to the solution of the torsion problem for an isotropic beam whose cross-section has a different boundary C, defined by equation f ðx; ZÞ ¼ 0 (see Fig. 2).
where s ¼ j tzx tzy j and C represents the boundary of the crosssection having equation f(x, y)= 0; then the components nx and ny of the outward normal vector n to the boundary C are related to the equation of the boundary by nx ¼ @f =@x and ny ¼ @f =@y. T
C
x Mz x
u
w
v y
z
nx n
C
C
x
zx zy
n
n
ny y
Fig. 1. De Saint Venant cylinder under torsion: (a) displacement field and (b) stress field.
y Fig. 2. Domain transformation.
ARTICLE IN PRESS 46
R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
As is customary the shear stress components can be written likewise in terms of Prandtl function @c @c ð16a; bÞ tzx ¼ ; tzy ¼ R @Z @x Relations between partial derivatives of warping and Prandtl functions in terms of new coordinates x and Z are given by @c @o @c @o ð17a; bÞ ¼ yGzx R Z ; ¼ yGzx R þx @Z @x @x @Z It follows that the Prandtl function c satisfies the Poisson equation in the form: @2 cðx; ZÞ 2
@x
þ
@2 cðx; ZÞ ¼ r2 cðx; ZÞ ¼ 2Gzx y in O @Z2
ð18Þ
while the corresponding boundary condition becomes cðx; ZÞ ¼ const on the boundary C. It may be also calculated [26] the torsional rigidity of a noisotropic cylinder in terms of the torsional rigidity of the corresponding isotropic cylinder by a substitution of Eqs. (12) in the expression of the resultant moment (see Eq. (9)) taking into account the coordinate transformation in Eqs. (10) and that the integration now extends over the region O bounded by the curve C.
analysis, it may be expanded in the double-ended Laurent series as F ðzÞ ¼
þ1 X
ak ðz z0 Þk ;
z ¼ x þ iZ
ð19Þ
where
wx ðx; ZÞ ¼ Re½F ðzÞ ¼ Rtzx ðx; ZÞ þRGzx yZ
ð20aÞ
wy ðx; ZÞ ¼ Im½F ðzÞ ¼ tzy ðx; ZÞ þ RGzx yx
ð20bÞ
F ðzÞ ¼ wx ðx; ZÞ þ iwy ðx; ZÞ ¼
bk Qk Þ þi
1 X
ðak Qk þbk Pk Þ
R
in O
ð22a; bÞ
that is the complex potential F ðzÞ in terms of wx(x, Z) and wy(x, Z) as defined in Eq. (19) satisfies both equilibrium (div s ¼ 0 in O Þ and compatibility in O . Following the steps developed for the isotropic case, since F ðzÞ is analytic in all the domain O taking full advantage of complex
ð25Þ
Thus, due to the relations (20) between sðx; ZÞ and F ðzÞ an approximation for the shear stresses tzx and tzy can be expressed in terms of truncated Laurent series of harmonic polynomials as follows: sup 1 X ða P bk Qk Þ Gzx yZ R k ¼ inf k k
By letting Pinf ðx; ZÞ ^ pðx; ZÞ ¼ P0 ðx; ZÞ ; ^ Psup ðx; ZÞ binf ^ b ¼ b0 ^ bsup
The Cauchy–Riemann conditions ð@wx =@x ¼ @wy =@Z; @wx =@Z ¼ @wy =@xÞ) for the functions wx and wy above defined, lead to the fulfillment of equilibrium and compatibility as
ðak Pk
k ¼ 1 k a 1
r2 wx ðx; ZÞ ¼ r2 tzx ðx; ZÞ ¼ 0 r2 wy ðx; ZÞ ¼ r2 tzy ðx; ZÞ ¼ 0
ð21Þ
1 X k ¼ 1 k a 1
tzx ðx; ZÞ ¼
@tzx ðx; ZÞ @tzy ðx; ZÞ þ ¼0 @Z @x @tzx ðx; ZÞ @tzy ðx; ZÞ R ¼ 2RGzx y @x @Z
ð24Þ
terms of harmonic polynomials in the form:
It is worth to note that as in the isotropic case (see [25]) the real and imaginary parts of the complex potential function F ðzÞ are related directly to the shear stresses; moreover the orthotropic ratio R appears in the expressions (20).Functions wx ðx; ZÞ and wy ðx; ZÞ defined in Eqs. (20) are harmonic in the transformed domain O in O
Qk ðx; ZÞ ¼ Imðx þ iZÞk
By letting in Eq. (23) ak ¼ ak þibk ðak ; bk A RÞ and assuming x0 = Z0 = 0 the complex potential function can be now expanded in
Starting from the condition that the shear stress components tzx and tzy are harmonic functions in the transformed domain, a complex potential function denoted as F ðzÞ can be introduced in the form F ðzÞ ¼ wx ðx; ZÞ þiwy ðx; ZÞ;
ð23Þ
where (z z0)k with k4 0 is analytic in the whole domain, while for ko0 it is also analytic with exception of the point z0, that is Pþ1 k called a pole of order k. In Eq. (23) the series k ¼ 0 ak ðz z0 Þ is called regular part and it is capable to express any analytic P k function everywhere. While the summation 1 k ¼ 1 ak ðz z0 Þ is called principal part and accounts for singularities in z0. As it will be evident in the following numerical studies, the regular part of the Laurent series, coinciding with the Taylor series, is sufficient to provide a correct solution for the case of simply connected crosssection. Inspection of multiply connected cross-sections shows as the principal part of the Laurent series with the exception of the term k= 1 cannot be neglected in order to find a correct shear stress distribution. The presence of the term k = 1 in the expansion returns multi-valued functions for both warping and Prandtl functions that is meaningless from a physical point of view. Powers (z)k will be denoted as Pk +iQk where Pk and Qk are the 2 2 so-called harmonic polynomials ðr Pk ¼ 0; r Qk ¼ 08kÞand for the case under exam they are defined in the transformed domain as follows: Pk ðx; ZÞ ¼ Reðx þ iZÞk ;
4. Line element-less method in the solution of torsion problem for orthotropic material
ak ; z0 ¼ x0 þ iZ0 A C
k ¼ 1
k a 1
tzy ðx; ZÞ ¼
sup X
ðak Qk þ bk Pk Þ þ RGzx yx
ð26a; bÞ
k ¼ inf k a 1
Qinf ðx; ZÞ ^ qðx; ZÞ ¼ Q0 ðx; ZÞ ; ^ Qsup ðx; ZÞ
ainf ^ a ¼ a0 ^ asup
;
ð27a2dÞ
Eqs. (26a,b) can be rewritten in a compact form as
sðx; ZÞ ¼ LDðx; ZÞw þGzx yg
ð28Þ
ARTICLE IN PRESS R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
where pT Dðx; ZÞ ¼ qT
qT ; pT
a w ¼ ; b
1=R L ¼ 0
0 ; 1
Z g¼ Rx ð29a2dÞ
As for the isotropic case, Eq. (28) completely fulfills equilibrium and compatibility equations in the transformed domain; however, since the series in Eqs. (26a,b) are truncated retaining the terms included between inf and sup then the boundary condition sT n ¼ 0 may not be fulfilled in each point of the new contour line C . In order to satisfy the free-stress boundary condition sT n ¼ 0 on C in a weak-form, according to the LEM procedure, the unknown series coefficients w can be evaluated by minimizing a functional in the form I I ð30Þ Iðw; yÞ ¼ ðsT nÞ2 dC ¼ ðtzx n x þ tzy n y Þ2 dC C
C
that is the squared value net flux of the shear stress vector s through the boundary C under the condition of static equilibrium equation: 8 H Iðw; yÞ ¼ ðsT nÞ2 dC ¼ min > > > w;y < C ð31a; bÞ subjected to > > R T > : O s g dO ¼ M z By making use of Green’s lemma and the properties of harmonic polynomials, the domain integral in Eq. (31b) extended to the cross-section O can be converted into a line integral on C as widely explained in Appendix A. Then the constrained minimum problem described in Eqs. (31a,b) may be solved by using the Lagrange multiplier method by performing the minimization of the free enlarged ~ y; lÞ expressed in the following form: functional Iðw; I T ~ Iðw; y; lÞ ¼ ðsT nÞ2 dC þ lðh w þðGzx =RÞyI p Mz Þ ¼ min ð32Þ w;y;l
C
where l is the Lagrange multiplier. Substitution of Eq. (28) into (32) leads to I ~ Iðw; y; lÞ ¼ ððwT DT LT þ Gzx ygT Þnn T ðLDw þGzx ygÞÞ dC C T
þ lðh w þ ðGzx =RÞyI p Mz Þ ¼ min w;y;l
ð33Þ
By performing variations with respect to the unknown coefficients w, Gzxy and l, respectively, we get @I~ ¼ Q w þ Gzx yq þ lh ¼ 0 @w
ð34aÞ
@I~ ¼ q T w þ Gzx yH þ lI P =R ¼ 0 @ðGzx yÞ
ð34bÞ
T @I~ ¼ h w þðGzx =RÞyI P Mz ¼ 0 @l
ð34cÞ
where I Q ¼ 2 ðDT LT nn T LDÞ dC ; C
I H ¼ 2 ðgT nn T gÞ dC C
I q ¼ 2 ðDT LT nn T gÞ dC ; C
ð35a2cÞ
47
It is worth to note that the problem as expressed in Eq. (33) parallels the torsion problem in the isotropic case solved by LEM in Di Paola et al. [25]. Differences rise up in the presence of the orthotropic ratio and the operators Q ; q and H are over-signed to bear in mind that all the integrals are performed in the transformed domain. Eq. (34a) in terms of the vector w is rewritten w ¼ Gzx yQ
1
q lQ
1
h
By inserting Eq. (36) in Eqs. (34b) and (34c) we get 8 < Gzx yðH q T Q 1 qÞ þ lðI P =R q T Q 1 hÞ ¼ 0 : G yðI =R h T Q 1 qÞ þ lðhT Q 1 hÞ ¼ M zx z P
ð36Þ
ð37Þ
Eqs. (37) represent a linear algebraic system in the unknowns Gzxy and l. Once the solution of Eq. (37) is obtained, Eq. (36) provides the Laurent series coefficients and the shear stress distribution (Eq. (28)) in the transformed domain O . In Appendix B expressions for warping function as well as for Prandtl function are reported. Following the procedure above explained the torsion problem has been solved in the plane xZ: in order to find the stress field as well as the warping and Prandtl functions of the orthotropic crosssection under study we have to pass through an inverse variable transformation taking into account Eqs. (10a,b). Once an explanation on the line element-less method has been presented, some insight into the main differences between LEM and the complex polynomial method is discussed. First of all the coordinates transformation expressed in Eqs. (10) allows to extend the CPM as well as the LEM to the case of orthotropic material. As reported in the mentioned Refs. [18,19], the complex-valued analytical function used in the CPM is composed as usual by the warping function and its harmonic conjugate. The values of this potential function are approximated by the values of the complex polynomial pk(z) on the contour of the cross-section. The order of pk(z) is equal to (k 1) where 2k nodal points are specified along the contour of the cross-section. As a result the solution of the boundary problem of the Laplace equation comes down to the determination of the coefficients of the complex polynomial pk(z) from the system of the linear equation by substituting nodal values of real and imaginary parts of the potential function in Re[pk(z)]and Im[pk(z)], respectively. As stated before, in the LEM as well as in the CPM a polynomial approximation is performed and the series coefficients derive from the solution of a linear algebraic equation but some differences between the two method rise up. Firstly, the potential function presented in the LEM is in terms of shear stresses allowing a straightforward extension to the combined torsion–flexure problem. The weak-form procedure proposed in the LEM allows to exceed the limit of the CPM which is sensitive to the nodal point placement where instability can occur for irregular nodal collocation on the boundary for even high polynomial orders. Moreover, the truncation order of the series is related on the number of nodal points fixed on the contour while in the LEM also very few series terms assents to reach accurate results. By using Lagrange multiplier method the static equivalence condition is incorporated in the solving functional and value of the unitary twist angle is derived together with the unknown series coefficients; the CPM needs to perform a domain integral once the system in the unknown coefficients is solved. At last since the CPM model assumes that the potential function is analytic on the boundary and in the domain, the polynomial expansion is performed in a Taylor series coinciding with the regular part of the Laurent series that is no polar
ARTICLE IN PRESS 48
R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
singularities can be considered. The deduced consequence is that the CPM cannot be applied to handle hollow sections or section with circular grooves whose solutions are provided by LEM application as shown in the following.
a
x
In this section some numerical studies have been reported in order to validate the accuracy and the correctness of the new formulation of the line element-less method here reported. Numerical applications have been conducted on simply and multiply connected domains. As previously pointed out it is worth to note that for the case of simply connected cross-sections the regular part of the Laurent series is sufficient to provide a correct solution and for all the cases whose the shear stress distribution is already known the LEM returns the exact solution. The LEM shows high versatility also in the case of material having a marked orthotropy. 5.1. Elliptical cross-section The variable transformation in Eqs. (10a,b) changes the boundary of an elliptical cross-section with semi-axes a and b in another elliptical domain whose semi-axes are aR and b, respectively. Thus, as for the isotropic case, the LEM applied in the transformed domain returns an exact solution. For any order of truncation sup, all the coefficients of the Laurent series are equal zero with exception of b1. Table 1 summarizes results in terms of unknown series coefficients, Lagrange multiplier l, shear stress field s, Prandtl c(x, y) and warping o(x, y) functions, respectively for the cross-section under exam. It is worth to note as the stress components are independent of the elastic constants and clash with the isotropic case (see [1]). Values of torsional stiffness and stresses at points A and B, as pointed out in the cross-section depicted in Table 1, are reported in Table 2. These values have been obtained by LEM and CVBEM with 52 nodes for full boundary discretization [24] and they are Table 1 Results for shear stresses for an elliptical cross-section. Cross-section geometry
B y Fig. 3. Rectangular cross-section.
Table 3 Comparison of results for an orthotropic rectangular cross-section (d = Rab 1 = 4).
LEM CVBEM [24] Lekhnitskii [1]
b1 = (a2R2 b2)Mz/a3R3pR 0 (a2R2 + b2)Mz/a3R3GzxpR2 2Mzy/pab3 2Mzx/pab3 (a2y2 + b2x2)/a3b3p (b2 a2R2)xy/(b2 + a2R2) a3b3pGzxR2/(b2 + a2R2)
Table 2 Comparison of results for an orthotropic elliptical cross-section ðR ¼
LEM CVBEM [24] Lekhnitskii [1]
pffiffiffi 2Þ.
Mz/y
jtðBÞ zx j=Gzx ya
jtðAÞ zy j=Gzx ya
0.761a4Gzx 0.7618a4Gzx 0.7616a4Gzx
1.0909 1.0825 1.0909
0.7273 0.7170 0.7273
16Mz/ya4Gzx
2jtðBÞ zx j=Gzx ya
2jtðAÞ zy j=Gzx ya
1.906 1.897 1.897
1.4946 1.4851 1.4963
3.5901 3.2028 3.3433
compared with the available analytical results [1]. LEM results validate the robustness of the method in returning exact solutions when they are known in analytic form as in the case of an elliptical cross-section.
5.2. Rectangular cross-section A rectangular cross-section of sides a and b as depicted in Fig. 3 changes in the transformed domain in another rectangle of different length sides namely Ra and b. For this section results obtained by LEM have been compared with existing results provided by a series expansion similar to that known for the isotropic case [27,28]. In particular, Lekhnitskii [1] suggests three factors namely b, k1, k2 in form of numerical series depending only on the ratio of the elastic moduli R and the ratio of the sides namely the factor d = Rab 1. By using these factors very simple formulas for the torsional rigidity and the maximum shearing stress at points A and B (see Fig. 3) can be expressed as follows: Mz =y ¼ Gzx ab3 b;
Series coefficients a 0 Lagrange multiplier l Unitary twist angle y Shear stresses tzx(x, y) Shear stresses tzy(x, y) Prandtl function c(x, y) Warping function o(x, y) Torsional rigidity Mz/y
b
A
5. Numerical applications
2 tðAÞ zx ¼ Mz k1 =ab ;
2 tðBÞ zy ¼ RMz k2 =ab
ð38Þ
and considered for a comparison. LEM applied to a rectangular cross-section returns all the coefficients of the Laurent series equal zero with exception of bk 8k ¼ 1; 3; 5; . . .. Numerical applications have been conducted for two different values of the factor d (d = 4,1) and for both cases the order of truncation has been selected sup = 8. First let us consider the case where the lengths of the sides and the orthotropic ratio are selected a= 4, b =3 and R= 3, respectively (d = 4). For the above orthotropic data LEM results are compared with the CVBEM results available in [24] and those obtained via Lekhnitskii series expansion in terms of torsional rigidity and shear stresses at the middles of the sides namely point A and B (Table 3). For the same orthotropic data the coefficients b, k1 and k2 have been computed for increasing order of truncation sup in the LEM and then compared with values provided by a boundary method [17] and the tabular data in [1] (Table 4). In Figs. 4 and 5 stresses due to torsion and contour lines of warping functions are plotted, respectively, for the isotropic (R= 1) and orthotropic (R= 3) cases under a unit torsion moment.
ARTICLE IN PRESS R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
49
It is worth to note that while the shearing stress in an isotropic bar is maximum at the middles of the longer sides, in an orthotropic bar, the points of maximum shearing stress may be at the middles either of the longer or of the shorter sides, depending on d = Rab 1. Let, for example, the lengths of the sides be a =4 and b= 2 (a = 2b) and let the orthotropic ratio be R= 2, that is d= 4. For
these selected values the maximum stress occurs in the point A at the middle of the shorter side b. As expected in the isotropic (R= 1) rectangular cross-section with the same geometry the maximum shear stress occurs at the middle of the longer side a. This difference in the stress distribution is evident in the stress maps depicted in Fig. 6.
Table 4 Comparison of BEM and LEM results for an orthotropic rectangular cross-section (d =Rab 1 = 4).
Table 5 Comparison of LEM and Lekhnitskii results for an orthotropic rectangular crosssection (d = Rab 1 = 1).
BEM [17]
LEM
Nel
sup
10 8 20 12 40 16 Lekhnitskii [1]
k1
b
k2
LEM
k1
b
bk a 0
BEM
LEM
BEM
LEM
BEM
LEM
Sup
bk a 0
4 6 8
0.270 0.280 0.281 0.281
0.282 0.281 0.281
3.702 3.561 3.552 3.550
3.529 3.518 3.557
2.556 2.694 2.686 2.644
2.826 2.647 2.610
8 2 12 3 16 4 Lekhnitskii [1]
k2
LEM
Err (%)
LEM
Err (%)
LEM
Err (%)
0.140 0.141 0.141 0.141
0.4 0.3 0.3
4.754 4.827 4.793 4.803
1 0.5 0.2
4.754 4.827 4.793 4.803
1 0.5 0.2
Fig. 4. Stresses due to torsion: (a) isotropic case R= 1 and (b) orthotropic case R = 3.
Fig. 5. Contour lines of warping function: (a) isotropic case R =1 and (b) orthotropic case R= 3.
Fig. 6. Stress map (a) isotropic case R= 1; (b) orthotropic case R =2.
ARTICLE IN PRESS 50
R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
As well for d = 1 the rectangular section changes in a square cross-section of length side aR = b whose all the Laurent coefficients are zero with exception of b3 ; b7 ; b11 ; . . . and a good accuracy is reached with a lower number of series coefficients involved in the solution. In Table 5 a comparison between LEM results and the tabular Lekhnitskii benchmark solution in terms of b, k1 and k2 is reported.
5.3. Annular cross-section The annular cross-section depicted in Fig. 7 in the transformed domain changes in a hollow elliptical cross-section whose semiaxes are ReR and Re for the external contour and RiR and Ri for the internal contour, respectively. For the annular cross-section under torque Mz the LEM (for x0 ¼ Z0 ¼ 0) returns only one coefficient of the Laurent series different from zero, namely b1, whatever the orders of truncation
Re
x
Ri
inf and sup are. It has to be emphasized that this is an exact solution since sT n ¼ 0 proves to be locally fulfilled in the external and internal contours whatever the value of the orthotropic ratio R is. Table 6 summarizes results in terms of unknown series coefficients, Lagrange multiplier l, shear stress field s, Prandtl c(x, y) and warping o(x, y) functions, respectively, where Re and Ri are the internal and external radius, respectively, and IP ¼ pðR2e R2i Þ=2. As for the elliptic cross-section, the stress distribution is independent of the value of the orthotropic ratio R and clash with the stress field provided for the isotropic case. Otherwise a closer inspection in the expression of the warping function o(x, y) (see Table 6) shows its dependence on the value of the orthotropic ratio R: by setting R= 1 (isotropic case) the cross-section remains plane and the well-known result for isotropic circular and annular cross-sections is restored. In Fig. 8 warping function and its contour lines are depicted for the case Gzx = 4Gzy. It is worth to note that although the annular cross-section is a multiply connected cross-section no polar singularities are involved in the solution and the regular part of the Laurent series is sufficient in order to provide a correct shear stress distribution in presence of a torsion moment as for a simply connected circular cross-section. Moreover by assuming Ri = 0 the coefficient b1 in Table 6 exactly coincides with the same coefficient in Table 1 by letting a = b=Re in the elliptical cross-section.
y
R Fig. 7. Annular cross-section.
Ri i x
Table 6 Results for shear stresses for an annular cross-section. Series coefficients a 0 Lagrange multiplier l Unitary twist angle y Shear stresses tzx(x, y) Shear stresses tzy(x, y) Prandtl function c(x, y) Warping function o(x, y) Torsional rigidity Mz/y
b1 = (R2 1)Mz/2IP(Re2 Ri2)R 0 (R2 + 1)Mz/2GzxIPR2 Mzy/IP(Re2 +Ri2) Mzx/IP(Re2 + Ri2) Mz(x2 + y2)/2IP(Re2 + Ri2) + const (R2 1)xy/(R2 + 1) 2pGzxIP(Re2 + Ri2)R2/(1 + R2)
Re 2
y
Fig. 9. Annular cross-section.
Fig. 8. Annular cross-section: (a) warping function and (b) its contour lines ðMz ¼ 1; Re ¼ 2; Ri ¼ Re =2; R ¼ 1=2Þ.
ARTICLE IN PRESS R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
5.4. Hollow circular cross-section If the centre of the hollow does not coincide with the centre of the external contour as depicted in Fig. 9 the pole of the principal part of the Laurent series moves from the point of coordinates x0 ¼ Z0 ¼ 0 to the point with coordinates x0 ¼ Re R=2 and Z0 ¼ 0
51
that is the centre of the elliptical hollow in the transformed domain. The numerical application studies the orthotropic case of Gzy =2Gzx with Re =2 and Ri =Re/4. Due to the symmetry LEM returns all the coefficients ak = 08k and bk a08k; in this case the solution is not exact because the condition sT n ¼ 0 is not locally
Fig. 10. Hollow circular cross-section: (a) Prandtl function and (b) its contour line ðMz ¼ 1; Re ¼ 2; Ri ¼ Re =4; R ¼
pffiffiffi 2Þ.
Fig. 11. Hollow circular cross-section ðMz ¼ 1; Re ¼ 2; Ri ¼ Re =4Þ: stress vector contour lines (a) isotropic R =1 and (b) orthotropic case R ¼
Fig. 12. Hollow circular cross-section ðMz ¼ 1; Re ¼ 2; Ri ¼ Re =4Þ: warping contour lines (a) isotropic R= 1 and (b) orthotropic case R ¼
pffiffiffi 2.
pffiffiffi 2.
ARTICLE IN PRESS 52
R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
fulfilled in all points of the contour. It is worth to note that the solution for a hollow elliptical cross-section, that is the shape assumed by the cross-section in the transformed domain, requires a computational effort greater than that required by the solution for a hollow circular cross-section. It follows that a greater number of polynomials is required to reach a good accuracy in the solution of the orthotropic case that is the lower order of truncation has been brought from inf =3 for the isotropic case to inf = 9 in the orthotropic case. In both cases the higher order of the series has been selected sup= 3. In Fig. 10 the Prandtl function for the case under exam is depicted with its contour lines. In Figs. 11 and 12 a comparison with the isotropic solution (R= 1) is reported in terms of stress vector and warping contour lines, respectively.
y take the following expressions a2 ¼ 0:07318Mz R2 =b and y ¼ 0:747597Mz R=b4 Gzx . Two different combinations of mechanical and geometrical parameters have been considered both satisfying the condition R ¼ a1 b ¼ a 1 b. In Figs. 14 and 15 the modulus of the stress vector functions and their contour lines are depicted for R= 0.5 and a = 2 and for R=2 and a =0.5, respectively, once b has been selected equal 1. In Figs. 16 and 17 Prandtl functions and their contour lines are depicted for the two data combinations under exam. For both cases LEM provides exact solutions as follows from the observation that the Prandtl functions take a constant value on the contour and their external contour lines move tangential along the boundary of the cross-section.
5.5. Elliptical cross-section with elliptical groove Let us consider an elliptical cross-section with an elliptical groove as illustrated in Fig. 13a. The coordinate transformation leads the original cross-section in a cross-section whose external contour is an ellipse of semi-axes aR and b, respectively; likewise the elliptical groove changes in an ellipse of semi-axes aR and b. It is worth to note that for a suitable choice of mechanical and geometrical data that is R ¼ a1 b ¼ a 1 b the cross-section in the transformed domain becomes a circular cross-section with a circular groove (see Fig. 13b) whose exact analytical solution exists [28]. In the numerical application the latter case is pursued. Also in this case the principal part of the Laurent series is essential in order to get a correct torsion solution and the pole is selected in the centre of the groove. Whatever the orders of series truncation are the only series coefficient different from zero is a 2 as expected. By selecting the qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi angle and g ¼ p=10 with b ¼ b sin g2 þð1 cos gÞ2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a ¼ a sin g2 þ ð1 cos gÞ2 the coefficient a 2 and the twist angle
aR = b
a
b x
b
6. Conclusions The extension of the novel numerical technique, labelled line element-less method (LEM), is presented in order to provide approximate solutions of the De Saint Venant torsion problem for orthotropic beams with simply and multiply connected crosssection. A suitable transformation of coordinates allows to restore the torsion problem as that studied in the isotropic case. A complex potential function analytic in all the domain whose real and imaginary parts are related to the shear stress components and to the orthotropic ratio has been introduced and expanded in the double-ended Laurent series involving harmonic polynomials. The potential function returns the domain governing equations, namely equilibrium and compatibility, while in order to satisfy the missing condition on the boundary an element-free weakform procedure has been proposed. By imposing that the square of the net flux of the shear stress across the border is minimum with respect to the series coefficients under the static equivalence condition the unknown series coefficients are thus provided. Numerical implementation of the LEM results in system of linear algebraic equations involving symmetric and positive-definite matrices. All the integrals are transferred into the boundary without requiring any discretization neither in the domain nor in the contour. The technique provides the complete shear stress field as shown by various numerical applications in order to assess the efficiency and the accuracy of the method to handle shear stress problems in presence of orthotropic material.
γγ
b
b Acknowledgement
y
a
aR = b
Fig. 13. (a) Elliptical cross-section with elliptical groove and (b) circular crosssection with circular groove.
The author wishes to thank Professor Mario Di Paola for stimulating discussions and suggestions on the subject of the paper.
Fig. 14. Elliptical cross-section with elliptical groove (R =0.5, b = 1, a= 2): (a) modulus of stress vector and (b) its contour lines.
ARTICLE IN PRESS R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
53
Fig. 15. Elliptical cross-section with elliptical groove (R = 2, b= 1, a= 0.5): (a) modulus of stress vector and (b) its contour lines.
Fig. 16. Elliptical cross-section with elliptical groove (R= 0.5, b =1, a= 2): (a) Prandtl function and (b) its contour lines.
where I P represents the polar inertia moment in the transformed domain O respect to the axes x and Z
Appendix A The static equivalence of the shear stresses as expressed in Eq. (31b) in terms of harmonic polynomials is fulfilled requiring that:
O
ðA:1Þ
O
or Z sup 1 X ak 2 R k ¼ inf k a 1
2
ðx þ Z Þ dO ¼
I
x3 3
O
nx þ
Z3 3
! ny
dC
ðA:3Þ
and
Z sup 1 X 1 þ 2 bk ðxPk þ ZQk Þ dO þ Gzx yI P ¼ Mz R R k ¼ inf k a 1
2
C
Z sup 1 X ak ðxQk ZPk Þ dO 2 R k ¼ inf k a 1
IP ¼
Z
Z sup 1 X div u k dO þ 2 bk R k ¼ inf O k a 1
1 div v k dO þ Gzx yI P ¼ Mz R O ðA:2Þ
ZPk þ 1 ðx x0 ; Z Z Þ=ðkþ 1Þ 0 u k ¼ xPk þ 1 ðx x0 ; Z Z Þ=ðk þ 1Þ ; 0 ZQk þ 1 ðx x0 ; Z Z Þ=ðk þ 1Þ 0 vk ¼ xQk þ 1 ðx x0 ; Z Z0 Þ=ðkþ 1Þ
ðA:4a; bÞ
Then, according to the Green lemma, the integrals in Eq. (A.2) may be performed as line boundary integrals as
ARTICLE IN PRESS 54
R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
Fig. 17. Elliptical cross-section with elliptical groove (R= 2, b = 1, a= 0.5): (a) Prandtl function and (b) its contour lines.
follows: sup X
1 a R2 k ¼ inf k
I
k a 1
¼
u Tk n dC þ
sup X
1 b R2 k ¼ inf k k a 1
C
sup X
sup X
ak c k þ
k ¼ inf k a 1
bk d k þ
k ¼ inf k a 1
I
v Tk n dC þ
form in the transformed domain: sup 1 X ak b Qk þ 1 ðx; ZÞ þ k Pk þ 1 ðx; ZÞ cðx; ZÞ ¼ kþ 1 R k ¼ inf kþ 1
Gzx yI p R
k a 1
C
Gzx yI p ¼ Mz R
ðA:5Þ
where ck ¼
1 R2
I
u Tk n dC ;
dk ¼
C
1 R2
I
v Tk n dC
ðA:6a; bÞ References
C
By collecting c k and d k in the vector h as T
h ¼ jc inf c 0 c sup d inf d 0 d sup j
ðA:7Þ
Eq. (A.5) can be expressed as follows: T
h w þ ðGzx =RÞyI p ¼ Mz
1 2 ðB:2Þ yGzx ðx þ Z2 Þ 2 Once the Prandtl function is known in the transformed domain, the simple inverse variable transformation (10a,b) allows to provide the Prandtl function for the cross-section under exam in the orthotropic case.
ðA:8Þ
Appendix B Bearing in mind the relationship between shear stresses and warping function in Eqs. (12a,b) and their series expansion in Eqs. (26a,b) and by taking advantage of the harmonic polynomial rules the warping function can be expressed in the transformed domain in terms of harmonic polynomials as follows: sup X 1 ak b Pk þ 1 ðx; ZÞ k Qk þ 1 ðx; ZÞ oðx; ZÞ ¼ ðB:1Þ 2 kþ1 yGzx R k ¼ inf k þ1 k a 1
Once the warping function is provided in the transformed domain, the simple inverse variable transformation (10a,b) allows to provide the warping function for the cross-section under exam in the orthotropic case. Similarly bearing in mind the relationship between shear stresses and Prandtl function in Eqs. (16a,b) and their series expansion in Eqs. (26a,b) the Prandtl function takes the following
[1] Lekhnitskii SG. Theory of elasticity of an anisotropic body. Moscow: Mir Publishers; 1981. [2] Tolf G. St. Venant bending of an orthotropic beam. Composite Structures 1985;4:1–14. [3] Herrmann RL. Elastic torsional analysis of irregular shapes. ASCE Journal of Engineering Mechanics Division 1965;91(EM6):11–9. [4] Mason WE, Herrmann RL. Elastic shear analysis of general prismatic shaped beams. ASCE Journal of Engineering Mechanics Division 1968;94(EM4): 965–983. [5] Kosmatka JB, Dong SB. Saint-Venant solutions for prismatic anisotropic beams. International Journal of Solids and Structures 1991;28(7):917–38. [6] Kosmatka JB. Flexure-torsion behavior of prismatic beams, part I: section properties via power series. AIAA Journal 1993;31(1):170–9. [7] Jawson MA, Ponter AR. An integral equation solution of the torsion problem. Proceedings of the Royal Society of London 1963;275(A):237–46. [8] Chou SI, Mohr JA. A boundary integral method for torsion of composite shafts. Res Mechanica 1990;29:41–56. [9] Katsikadelis JT, Sapountzakis EJ. Torsion of composite bars by boundary element method. ASCE Journal of the Engineering Mechanics Division 1985;111(9):1197–210. [10] Sapountzakis EJ, Mokos VG. Nonuniform torsion of composite bars by boundary element method. Journal of Engineering Mechanics, ASCE 2001;127(9):945–53. [11] Sapountzakis EJ. Nonuniform torsion of multi-material composite bars by the boundary element method. Computers and Structures 2001;79(32):2805–16. [12] Sapountzakis EJ, Mokos VG. Warping shear stresses in nonuniform torsion of composite bars by BEM. Computer Methods in Applied Mechanics and Engineering 2003;192(39–40):4337–53. [13] Sapountzakis EJ, Mokos VG. Nonuniform torsion of composite bars of variable thickness by BEM. International Journal of Solids and Structures 2004;41(7):1753–71. [14] Friedman Z, Kosmatka JB. Torsion and flexure of a prismatic isotropic beam using the boundary element method. Computers & Structures 2000;74: 479–494.
ARTICLE IN PRESS R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
[15] Petrolo S, Casciaro R. 3D beam element based on Saint Vena nt’s rod theory. Computers & Structures 2004;82:2471–81. [16] Gracia L, Doblare M. Shape optimization of elastic orthotropic shafts under torsion by using boundary elements. Composite Structures 1988;30(6): 1281–1291. [17] Gaspari D, Aristodemo M. Torsion and flexure analysis of orthotropic beams by a boundary element model. Engineering Analysis with Boundary Elements 2005;29:850–8. [18] Hromadka II TV, Guymon GL. Complex polynomial approximation of the Laplace equation. Journal of Hydraulic Engineering ASCE 1984;110(3):329–39. [19] Poler AC, Bohannon AW, Schowalter SJ, Hromadka II TV. Using the complex polynomial method with mathematica to model problems involving the Laplace and Poisson equations. Published online in Wiley InterScience /www.interscience.wiley.com) doi:10.1002/num.20365S. [20] Hromadka II TV, Lai Ch. The complex variable boundary element method in engineering analysis. New York: Springer; 1987. [21] Hromadka II TV, Whitley RJ. Advances in the complex variable boundary element method. London: Springer-Verlag; 1998.
55
[22] Aleynikov SM, Stromov AV. Comparison of complex methods for numerical solutions of boundary problems of the Laplace equation. Engineering Analysis with Boundary Elements 2004;28:615–22. [23] Hromadka II TV, Pardoeon G. Application of the CVBEM to nonuniform St. Venant torsion. Computer Methods in Applied Mechanics and Engineering 1985;53:149–61. [24] Dumir PC, Kumar R. Complex variable boundary element method for torsion of anisotropic bars. Applied Mathematical Modelling 1993;17: 80–88. [25] Di Paola M, Pirrotta A, Santoro R. Line element-less method (LEM) for beam torsion solution (TRULY NO-MESH METHOD). Acta Mechanica 2008;195(1–4):349–64. [26] Sokolnikoff IS. Mathematical theory of elasticity. New York, NY: McGrawHill; 1956. [27] Muskhelishvili NI. Some basic problems of the mathematical theory of elasticity. Groningen, The Netherlands: P. Noordhoff Ltd.; 1953. [28] Timoshenko SP, Goodier JN. Theory of elasticity. New York, NY: McGraw-Hill; 1970.
Contents lists available at ScienceDirect
International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci
The line element-less method analysis of orthotropic beam for the De Saint Venant torsion problem Roberta Santoro Dipartimento di Ingegneria Strutturale, Aerospaziale e Geotecnica DISAG, Universita degli Studi di Palermo, Viale delle Scienze, 90128 Palermo, Italy
a r t i c l e in fo
abstract
Article history: Received 3 February 2009 Received in revised form 22 September 2009 Accepted 7 October 2009 Available online 13 October 2009
This paper deals with the extension of a novel numerical technique, labelled line element-less method (LEM), in order to provide approximate solutions of the De Saint Venant torsion problem for orthotropic beams having simply and multiply connected cross-section. A suitable transformation of coordinates allows to take full advantage of the theory of analytic complex functions as in the isotropic case. A complex potential function analytic in all the transformed domain whose real and imaginary parts are related to the shear stress components and to the orthotropic ratio is introduced and expanded in the double-ended Laurent series involving harmonic polynomials. An element-free weak-form procedure has been proposed imposing that the square of the net flux of the shear stress across the border is minimized with respect to the series coefficients. Numerical implementation of the LEM results in system of linear algebraic equations involving symmetric and positive-definite matrices. All the integrals are transferred into the boundary without requiring any discretization neither in the domain nor in the contour. The technique provides the complete shear stress field as shown by some numerical applications in order to assess the efficiency and the accuracy of the method to handle shear stress problems in presence of orthotropic material. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Orthotropic material Analytic functions Torsion problem
1. Introduction Several materials used in construction show a sharp difference in elastic properties for different directions. Examples of such materials exhibiting elastic anisotropic behaviour are provided by natural wood, composites and synthetic materials used in aircraft construction. To design anisotropic members undergoing elastic strains, it is necessary to know how to determine the stresses and strains in an anisotropic body theoretically, i.e. to solve torsion problems of the theory of elasticity for anisotropic body as presented in this paper. Benchmark solutions for torsion and flexure problems in the analysis of orthotropic beams are represented by the series expansion solutions due to Lekhnitskii [1] and Tolf [2], respectively. The commonly used techniques to pursue numerical solution of the orthotropic De Saint Venant problem are based on the finite element method (FEM) and boundary element method (BEM). A finite element approach has been proposed in Herrmann [3] and Mason and Herrmann [4] for the two-dimensional torsion and flexure analysis of beams with arbitrary cross-sections restricted to the homogeneous isotropic case. This method was also
Tel.: + 39 091 6568406; fax: +39 091 6568407.
E-mail addresses: [email protected], [email protected]. 0020-7403/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmecsci.2009.10.002
extended to provide solutions for beams made by anisotropic material by Kosmatka and Dong [5]. A Ritz-based power series method was used for the solution of the flexure-torsion problem of a general anisotropic material with arbitrary cross-section as reported in Ref. [6]. Both the above methods require the whole cross-section to be discretized into area elements. Although the solving system obtained via FEM involves symmetric and positivedefinite matrices providing simple and immediate solutions, a large number of elements is required for a good accuracy; furthermore, the limitation with respect to the shape of the elements leads to cumbersome meshing process. In the Ritz-based methods the elements of the mesh are less limited with respect to their shape because they are created only to aid the numerical integration on the cross-sectional area but the integration requires a large number of integration points because of the high order of the power series. On the other hand for its features BEM allows a more accurate evaluation of the stress field over the cross-section. Advantage of the latter method is related with an easier discretization of the boundary resulting in a relatively small number of line elements with the possibility of addressing the computational resources to restricted regions of interest. The use of boundary integral method for the torsion problem for homogeneous and isotropic material has been proposed by Jawson and Ponter [7] by using a simple constant interpolation and by Chou and Mohr [8] and Katsikadelis and Sapountzakis [9] to handle composite materials. In Sapount-
ARTICLE IN PRESS 44
R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
zakis et al. [10–13] the BEM is employed to solve the non-uniform torsion problem of composite bars of arbitrary variable cross-sections. In Ref. [14] a more refined analysis is presented with a three-node isoparametric boundary element to allow an accurate geometric representation of cross-sections having curved boundaries. In the latter paper the stress field over the crosssection is described at points quite far from the boundary probably due to the difficulties implied by numerical integration. Solution of De Saint Venant torsion problem obtained via FEM and BEM has also been used to evaluate torsion factors [15]. Within the orthotropic analysis Gracia and Doblare [16] solved an optimization problem for sections under torsion. A boundary element model in the analysis of De Saint Venant orthotropic problem is presented in [17], where the differential equations governing the shear stress field are converted in three Neumann problems which form the basis of the boundary integral formulation and a quadratic B-spline approximation represents the boundary variables. The application of complex variables in the boundary element methods allows to obtain the solution for two conjugate functions at once. These related methods intend to solve boundary problems of Laplace equation and distinguish in complex polynomial method (CPM) [18,19] and complex variable boundary element method (CVBEM) [20,21]. Both the above methods [22] provide solution of two-dimensional potential problems as the torsion of elastic rods making use of complex variables by generalization of the Cauchy integral formula into a boundary integral equation method. The CVBEM was also proposed for non-uniform torsion with isotropic homogeneous material [23] and for pure torsion of orthotropic prismatic bar [24]. Recently Di Paola et al. [25] propose to solve the problem of pure torsion for the De Saint Venant cylinder by using a novel approach labelled line element-less method (LEM). As in the CVBEM and CPM the method firstly works in terms of a potential function. In order to avoid to calculate both warping function and its harmonic conjugate a new potential function involving directly the shear stresses has been introduced. This complex function, holomorphic in all the domain, allows to take full advantage of the double-ended Laurent series involving harmonic polynomials. The selected expression of this analytic complex function is such that governing equations of De Saint Venant elastic problem, namely equilibrium and compatibility equations, remain fulfilled in the cross-section domain. While in the BEM theory the shape functions are the fundamental ones, in the proposed method the shape functions are harmonic polynomials that satisfy the field equations in the whole domain. The coefficients of the Laurent series are then evaluated by imposing a weak-form condition that the square value of the total net flux across the border remains minimum with respect to parameters expansion under the condition that the static equivalence is satisfied. Numerical implementation of LEM results in systems of linear algebraic equations in positive-defined and symmetric matrices solving only contours integral. The method provides exact solutions when existing and very accurate results by using very few series terms in the other cases. In this paper the above explained method is extended to the case of orthotropic material by means of a suitable transformation of coordinates [26] with the introduction of a complex potential function analytic in all the domain whose real and imaginary parts are related to the shear stress components and to the orthotropic ratio. Some numerical applications to simply and multiply connected domains have been reported to demonstrate the accuracy and the efficiency of the proposed LEM to handle shear stress problem in presence of orthotropic material.
2. Introductory remarks By definition, a material having two or three mutually orthogonal planes of elastic symmetry at each point (or, what is the same thing, two or three mutually perpendicular principal directions) is said to be orthogonally anisotropic or, more shortly, orthotropic. Thus, its mechanical properties are, in general, different along the directions of each of the axes. One common example of an orthotropic material with two axes of symmetry would be a polymer reinforced by parallel glass or graphite fibres. The strength and stiffness of such a composite material will usually be greater in a direction parallel to the fibres than in the transverse direction. Such materials are sometimes called transverse isotropic. A familiar example of an orthotropic material with three mutually perpendicular axes is provided by natural wood, in which the properties (such as strength and stiffness) along its grain and in each of the two perpendicular directions are different. Orthotropic materials require 9 independent variables (i.e. elastic constants) in their constitutive matrices. In contrast, a material without any planes of symmetry is fully anisotropic and requires 21 elastic constants, whereas a material with an infinite number of symmetry planes (i.e. every plane is a plane of symmetry) is isotropic, and requires only 2 elastic constants. The 9 independent elastic constants in orthotropic constitutive equations are comprised of three Young’s moduli Ex, Ey, Ez, the three Poisson’s ratios nzx, nzy, nxy, and the three shear moduli Gzx, Gzy, Gxy. By choosing the axes normal to the three planes of symmetry, the compliance matrix takes the form ex 1=Ex nyx =Ey nzx =Ez 0 0 0 sx e 1=Ey nzy =Ez 0 0 0 sy y nxy =Ex ez nxz =Ex nyz =Ey 1=Ez 0 0 0 sz g ¼ 0 0 0 1=Gyz 0 0 tyz yz gxz 0 t 0 0 0 0 1=G xz xz 0 0 0 0 0 1=Gxy txy gxy ð1Þ It is worth to be
nij Ei
¼
nji
ð2Þ
Ej
where generally nij a nji. Assuming sx = sy = txy =0 as in the De Saint Venant problem the constitutive equations become:
ex ¼
nzx Ez
gxy ¼ 0;
sz ; gyz ¼
ey ¼ tyz Gyz
;
nzy Ez
sz ;
gxz ¼
ez ¼
1 sz Ez
txz
ð3aÞ
ð3bÞ
Gxz
The equilibrium equations are the same of the isotropic case @tzx ¼ 0; @z
@tzy ¼ 0; @z
@tzx @tzy @sz þ þ ¼0 @x @y @z
ð4a2cÞ
The compatibility equations become 2
@ ex @2 ey @2 ex @2 ez þ 2 ¼ 0; þ 2 ¼ 0; 2 @y @x @z2 @x @2 ex 1 @ @gzx @gyz @2 ey ¼ ; ¼ @y@z @x @x@z 2 @x @y @2 ez ¼0 @x@y
@2 ey @2 ez þ 2 ¼0 @z2 @y 1 @ @gyz @gxz ; @y 2 @y @x ð5Þ
ARTICLE IN PRESS R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
and bearing in mind Eqs. (3) and (4) can be rewritten in terms of stresses as follows:
mzx
@2 sz @2 sz þ mzy ¼ 0; 2 @y @x2
@2 sz ¼ 0; @x2
mzx @2 sz Ez
mzy Ez
@2 sz ¼ 0; @y2 !
@2 sz ¼0 @x@y
1 1 @2 tzx 1 @2 tyz ; ¼ @y@z 2 Gxz @x@y Gyz @x2 ! @2 sz 1 1 @2 tzy 1 @2 txz ¼ @x@z 2 Gyz @x@y Gxz @y2
ð6Þ
ðtzy x tzx yÞ dO ¼ Mz
Let us consider a beam of arbitrary cross-section having three mutually orthogonal planes of elastic symmetry twisted by a moment Mz. It is assumed that the axis of the beam is perpendicular to one of these planes. As in the isotropic case, we let the axis of the beam coincide with the z-axis and choose the x- and y-axis in one end of the beam oriented as the material directions of orthotropy (see Fig. 1). It is worth to note that in the most general case the two longitudinal planes of elastic symmetry, namely xz and yz, do not necessarily intersect the plane of the section at a given point along lines parallel to the principal axes of the cross-section [26,27]. The shear moduli associated with the z- and x-axis and z- and y-axis are denoted by Gzx and Gzy, respectively. We assume, as in the isotropic case, that the components of displacement in the directions of x, y and z axes denoted as u, v and w, respectively, are given by the formulas: w ¼ yoðx; yÞ
ð7a2cÞ
where y is the twist angle for unit length of the bar and o(x, y) is the warping function associated with this problem. It is easy to verify [26] that in the orthotropic case, the torsion function o no longer satisfies Laplace equation in the domain O of the crosssection. As in the isotropic case, the non-vanishing stress components should satisfy the free-traction boundary conditions on the contour of the cross-section, namely
tzx nx þ tzy ny ¼ sT n ¼ 0 on C
ð9Þ
has to be satisfied. A simplified compact formulation of the De Saint Venant orthotropic torsion problem can be obtained following Sokolnikoff [26] by the introduction of new independent variables x and Z defined as follows:
x ¼ Rx;
3. Formulation of the pure torsion problem for an orthotropic beam
v ¼ yzx;
Moreover the static equivalence Z O
Once equilibrium, compatibility and constitutive law are reported for the orthotropic De Saint Venant beam in the next section the problem of pure torsion will be presented.
u ¼ yzy;
45
ð8Þ
Z¼y
ð10a; bÞ
where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R ¼ Gzy =Gzz
ð11Þ
is defined as the orthotropic ratio. As it will be evident in the following the variable transformation in Eq. (10) allows to formulate the torsion orthotropic problem as a Neumann or Dirichlet problem in the transformed domain where the governing equations are of the same form as that appearing in the study of torsion of isotropic cylinder. The shear stress components in the transformed domain can be easily written in terms of warping function: @o @o ð12a; bÞ tzx ¼ yGzx R Z ; tzy ¼ yRGzx R þx @x @Z and the third equilibrium equation, namely div s ¼ 0, becomes: @2 oðx; ZÞ @x
2
þ
@2 oðx; ZÞ ¼ r2 oðx; ZÞ ¼ 0 @Z2
in O
ð13Þ
It follows that o(x, Z) is an harmonic function in the transformed domain O whose equation is changed in f ðx; ZÞ ¼ 0. Making the corresponding change of variables in the boundary condition (8) it returns:
tzx n x þ tzy n y ¼ sT n ¼ 0 on C
ð14Þ
or in terms of warping function ! @o @f @o @f @f @f R ¼Z x þ @x @x @Z @Z @x @Z
on C
ð15Þ
where plainly C is the transformed boundary. It is worth to note that the components of the normal n to the boundary C are proportional to R@f =@xand @f =@Z, respectively. Hence the solution of the torsion problem for a cylinder of no-isotropic material whose cross-section has contour C is reduced to the solution of the torsion problem for an isotropic beam whose cross-section has a different boundary C, defined by equation f ðx; ZÞ ¼ 0 (see Fig. 2).
where s ¼ j tzx tzy j and C represents the boundary of the crosssection having equation f(x, y)= 0; then the components nx and ny of the outward normal vector n to the boundary C are related to the equation of the boundary by nx ¼ @f =@x and ny ¼ @f =@y. T
C
x Mz x
u
w
v y
z
nx n
C
C
x
zx zy
n
n
ny y
Fig. 1. De Saint Venant cylinder under torsion: (a) displacement field and (b) stress field.
y Fig. 2. Domain transformation.
ARTICLE IN PRESS 46
R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
As is customary the shear stress components can be written likewise in terms of Prandtl function @c @c ð16a; bÞ tzx ¼ ; tzy ¼ R @Z @x Relations between partial derivatives of warping and Prandtl functions in terms of new coordinates x and Z are given by @c @o @c @o ð17a; bÞ ¼ yGzx R Z ; ¼ yGzx R þx @Z @x @x @Z It follows that the Prandtl function c satisfies the Poisson equation in the form: @2 cðx; ZÞ 2
@x
þ
@2 cðx; ZÞ ¼ r2 cðx; ZÞ ¼ 2Gzx y in O @Z2
ð18Þ
while the corresponding boundary condition becomes cðx; ZÞ ¼ const on the boundary C. It may be also calculated [26] the torsional rigidity of a noisotropic cylinder in terms of the torsional rigidity of the corresponding isotropic cylinder by a substitution of Eqs. (12) in the expression of the resultant moment (see Eq. (9)) taking into account the coordinate transformation in Eqs. (10) and that the integration now extends over the region O bounded by the curve C.
analysis, it may be expanded in the double-ended Laurent series as F ðzÞ ¼
þ1 X
ak ðz z0 Þk ;
z ¼ x þ iZ
ð19Þ
where
wx ðx; ZÞ ¼ Re½F ðzÞ ¼ Rtzx ðx; ZÞ þRGzx yZ
ð20aÞ
wy ðx; ZÞ ¼ Im½F ðzÞ ¼ tzy ðx; ZÞ þ RGzx yx
ð20bÞ
F ðzÞ ¼ wx ðx; ZÞ þ iwy ðx; ZÞ ¼
bk Qk Þ þi
1 X
ðak Qk þbk Pk Þ
R
in O
ð22a; bÞ
that is the complex potential F ðzÞ in terms of wx(x, Z) and wy(x, Z) as defined in Eq. (19) satisfies both equilibrium (div s ¼ 0 in O Þ and compatibility in O . Following the steps developed for the isotropic case, since F ðzÞ is analytic in all the domain O taking full advantage of complex
ð25Þ
Thus, due to the relations (20) between sðx; ZÞ and F ðzÞ an approximation for the shear stresses tzx and tzy can be expressed in terms of truncated Laurent series of harmonic polynomials as follows: sup 1 X ða P bk Qk Þ Gzx yZ R k ¼ inf k k
By letting Pinf ðx; ZÞ ^ pðx; ZÞ ¼ P0 ðx; ZÞ ; ^ Psup ðx; ZÞ binf ^ b ¼ b0 ^ bsup
The Cauchy–Riemann conditions ð@wx =@x ¼ @wy =@Z; @wx =@Z ¼ @wy =@xÞ) for the functions wx and wy above defined, lead to the fulfillment of equilibrium and compatibility as
ðak Pk
k ¼ 1 k a 1
r2 wx ðx; ZÞ ¼ r2 tzx ðx; ZÞ ¼ 0 r2 wy ðx; ZÞ ¼ r2 tzy ðx; ZÞ ¼ 0
ð21Þ
1 X k ¼ 1 k a 1
tzx ðx; ZÞ ¼
@tzx ðx; ZÞ @tzy ðx; ZÞ þ ¼0 @Z @x @tzx ðx; ZÞ @tzy ðx; ZÞ R ¼ 2RGzx y @x @Z
ð24Þ
terms of harmonic polynomials in the form:
It is worth to note that as in the isotropic case (see [25]) the real and imaginary parts of the complex potential function F ðzÞ are related directly to the shear stresses; moreover the orthotropic ratio R appears in the expressions (20).Functions wx ðx; ZÞ and wy ðx; ZÞ defined in Eqs. (20) are harmonic in the transformed domain O in O
Qk ðx; ZÞ ¼ Imðx þ iZÞk
By letting in Eq. (23) ak ¼ ak þibk ðak ; bk A RÞ and assuming x0 = Z0 = 0 the complex potential function can be now expanded in
Starting from the condition that the shear stress components tzx and tzy are harmonic functions in the transformed domain, a complex potential function denoted as F ðzÞ can be introduced in the form F ðzÞ ¼ wx ðx; ZÞ þiwy ðx; ZÞ;
ð23Þ
where (z z0)k with k4 0 is analytic in the whole domain, while for ko0 it is also analytic with exception of the point z0, that is Pþ1 k called a pole of order k. In Eq. (23) the series k ¼ 0 ak ðz z0 Þ is called regular part and it is capable to express any analytic P k function everywhere. While the summation 1 k ¼ 1 ak ðz z0 Þ is called principal part and accounts for singularities in z0. As it will be evident in the following numerical studies, the regular part of the Laurent series, coinciding with the Taylor series, is sufficient to provide a correct solution for the case of simply connected crosssection. Inspection of multiply connected cross-sections shows as the principal part of the Laurent series with the exception of the term k= 1 cannot be neglected in order to find a correct shear stress distribution. The presence of the term k = 1 in the expansion returns multi-valued functions for both warping and Prandtl functions that is meaningless from a physical point of view. Powers (z)k will be denoted as Pk +iQk where Pk and Qk are the 2 2 so-called harmonic polynomials ðr Pk ¼ 0; r Qk ¼ 08kÞand for the case under exam they are defined in the transformed domain as follows: Pk ðx; ZÞ ¼ Reðx þ iZÞk ;
4. Line element-less method in the solution of torsion problem for orthotropic material
ak ; z0 ¼ x0 þ iZ0 A C
k ¼ 1
k a 1
tzy ðx; ZÞ ¼
sup X
ðak Qk þ bk Pk Þ þ RGzx yx
ð26a; bÞ
k ¼ inf k a 1
Qinf ðx; ZÞ ^ qðx; ZÞ ¼ Q0 ðx; ZÞ ; ^ Qsup ðx; ZÞ
ainf ^ a ¼ a0 ^ asup
;
ð27a2dÞ
Eqs. (26a,b) can be rewritten in a compact form as
sðx; ZÞ ¼ LDðx; ZÞw þGzx yg
ð28Þ
ARTICLE IN PRESS R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
where pT Dðx; ZÞ ¼ qT
qT ; pT
a w ¼ ; b
1=R L ¼ 0
0 ; 1
Z g¼ Rx ð29a2dÞ
As for the isotropic case, Eq. (28) completely fulfills equilibrium and compatibility equations in the transformed domain; however, since the series in Eqs. (26a,b) are truncated retaining the terms included between inf and sup then the boundary condition sT n ¼ 0 may not be fulfilled in each point of the new contour line C . In order to satisfy the free-stress boundary condition sT n ¼ 0 on C in a weak-form, according to the LEM procedure, the unknown series coefficients w can be evaluated by minimizing a functional in the form I I ð30Þ Iðw; yÞ ¼ ðsT nÞ2 dC ¼ ðtzx n x þ tzy n y Þ2 dC C
C
that is the squared value net flux of the shear stress vector s through the boundary C under the condition of static equilibrium equation: 8 H Iðw; yÞ ¼ ðsT nÞ2 dC ¼ min > > > w;y < C ð31a; bÞ subjected to > > R T > : O s g dO ¼ M z By making use of Green’s lemma and the properties of harmonic polynomials, the domain integral in Eq. (31b) extended to the cross-section O can be converted into a line integral on C as widely explained in Appendix A. Then the constrained minimum problem described in Eqs. (31a,b) may be solved by using the Lagrange multiplier method by performing the minimization of the free enlarged ~ y; lÞ expressed in the following form: functional Iðw; I T ~ Iðw; y; lÞ ¼ ðsT nÞ2 dC þ lðh w þðGzx =RÞyI p Mz Þ ¼ min ð32Þ w;y;l
C
where l is the Lagrange multiplier. Substitution of Eq. (28) into (32) leads to I ~ Iðw; y; lÞ ¼ ððwT DT LT þ Gzx ygT Þnn T ðLDw þGzx ygÞÞ dC C T
þ lðh w þ ðGzx =RÞyI p Mz Þ ¼ min w;y;l
ð33Þ
By performing variations with respect to the unknown coefficients w, Gzxy and l, respectively, we get @I~ ¼ Q w þ Gzx yq þ lh ¼ 0 @w
ð34aÞ
@I~ ¼ q T w þ Gzx yH þ lI P =R ¼ 0 @ðGzx yÞ
ð34bÞ
T @I~ ¼ h w þðGzx =RÞyI P Mz ¼ 0 @l
ð34cÞ
where I Q ¼ 2 ðDT LT nn T LDÞ dC ; C
I H ¼ 2 ðgT nn T gÞ dC C
I q ¼ 2 ðDT LT nn T gÞ dC ; C
ð35a2cÞ
47
It is worth to note that the problem as expressed in Eq. (33) parallels the torsion problem in the isotropic case solved by LEM in Di Paola et al. [25]. Differences rise up in the presence of the orthotropic ratio and the operators Q ; q and H are over-signed to bear in mind that all the integrals are performed in the transformed domain. Eq. (34a) in terms of the vector w is rewritten w ¼ Gzx yQ
1
q lQ
1
h
By inserting Eq. (36) in Eqs. (34b) and (34c) we get 8 < Gzx yðH q T Q 1 qÞ þ lðI P =R q T Q 1 hÞ ¼ 0 : G yðI =R h T Q 1 qÞ þ lðhT Q 1 hÞ ¼ M zx z P
ð36Þ
ð37Þ
Eqs. (37) represent a linear algebraic system in the unknowns Gzxy and l. Once the solution of Eq. (37) is obtained, Eq. (36) provides the Laurent series coefficients and the shear stress distribution (Eq. (28)) in the transformed domain O . In Appendix B expressions for warping function as well as for Prandtl function are reported. Following the procedure above explained the torsion problem has been solved in the plane xZ: in order to find the stress field as well as the warping and Prandtl functions of the orthotropic crosssection under study we have to pass through an inverse variable transformation taking into account Eqs. (10a,b). Once an explanation on the line element-less method has been presented, some insight into the main differences between LEM and the complex polynomial method is discussed. First of all the coordinates transformation expressed in Eqs. (10) allows to extend the CPM as well as the LEM to the case of orthotropic material. As reported in the mentioned Refs. [18,19], the complex-valued analytical function used in the CPM is composed as usual by the warping function and its harmonic conjugate. The values of this potential function are approximated by the values of the complex polynomial pk(z) on the contour of the cross-section. The order of pk(z) is equal to (k 1) where 2k nodal points are specified along the contour of the cross-section. As a result the solution of the boundary problem of the Laplace equation comes down to the determination of the coefficients of the complex polynomial pk(z) from the system of the linear equation by substituting nodal values of real and imaginary parts of the potential function in Re[pk(z)]and Im[pk(z)], respectively. As stated before, in the LEM as well as in the CPM a polynomial approximation is performed and the series coefficients derive from the solution of a linear algebraic equation but some differences between the two method rise up. Firstly, the potential function presented in the LEM is in terms of shear stresses allowing a straightforward extension to the combined torsion–flexure problem. The weak-form procedure proposed in the LEM allows to exceed the limit of the CPM which is sensitive to the nodal point placement where instability can occur for irregular nodal collocation on the boundary for even high polynomial orders. Moreover, the truncation order of the series is related on the number of nodal points fixed on the contour while in the LEM also very few series terms assents to reach accurate results. By using Lagrange multiplier method the static equivalence condition is incorporated in the solving functional and value of the unitary twist angle is derived together with the unknown series coefficients; the CPM needs to perform a domain integral once the system in the unknown coefficients is solved. At last since the CPM model assumes that the potential function is analytic on the boundary and in the domain, the polynomial expansion is performed in a Taylor series coinciding with the regular part of the Laurent series that is no polar
ARTICLE IN PRESS 48
R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
singularities can be considered. The deduced consequence is that the CPM cannot be applied to handle hollow sections or section with circular grooves whose solutions are provided by LEM application as shown in the following.
a
x
In this section some numerical studies have been reported in order to validate the accuracy and the correctness of the new formulation of the line element-less method here reported. Numerical applications have been conducted on simply and multiply connected domains. As previously pointed out it is worth to note that for the case of simply connected cross-sections the regular part of the Laurent series is sufficient to provide a correct solution and for all the cases whose the shear stress distribution is already known the LEM returns the exact solution. The LEM shows high versatility also in the case of material having a marked orthotropy. 5.1. Elliptical cross-section The variable transformation in Eqs. (10a,b) changes the boundary of an elliptical cross-section with semi-axes a and b in another elliptical domain whose semi-axes are aR and b, respectively. Thus, as for the isotropic case, the LEM applied in the transformed domain returns an exact solution. For any order of truncation sup, all the coefficients of the Laurent series are equal zero with exception of b1. Table 1 summarizes results in terms of unknown series coefficients, Lagrange multiplier l, shear stress field s, Prandtl c(x, y) and warping o(x, y) functions, respectively for the cross-section under exam. It is worth to note as the stress components are independent of the elastic constants and clash with the isotropic case (see [1]). Values of torsional stiffness and stresses at points A and B, as pointed out in the cross-section depicted in Table 1, are reported in Table 2. These values have been obtained by LEM and CVBEM with 52 nodes for full boundary discretization [24] and they are Table 1 Results for shear stresses for an elliptical cross-section. Cross-section geometry
B y Fig. 3. Rectangular cross-section.
Table 3 Comparison of results for an orthotropic rectangular cross-section (d = Rab 1 = 4).
LEM CVBEM [24] Lekhnitskii [1]
b1 = (a2R2 b2)Mz/a3R3pR 0 (a2R2 + b2)Mz/a3R3GzxpR2 2Mzy/pab3 2Mzx/pab3 (a2y2 + b2x2)/a3b3p (b2 a2R2)xy/(b2 + a2R2) a3b3pGzxR2/(b2 + a2R2)
Table 2 Comparison of results for an orthotropic elliptical cross-section ðR ¼
LEM CVBEM [24] Lekhnitskii [1]
pffiffiffi 2Þ.
Mz/y
jtðBÞ zx j=Gzx ya
jtðAÞ zy j=Gzx ya
0.761a4Gzx 0.7618a4Gzx 0.7616a4Gzx
1.0909 1.0825 1.0909
0.7273 0.7170 0.7273
16Mz/ya4Gzx
2jtðBÞ zx j=Gzx ya
2jtðAÞ zy j=Gzx ya
1.906 1.897 1.897
1.4946 1.4851 1.4963
3.5901 3.2028 3.3433
compared with the available analytical results [1]. LEM results validate the robustness of the method in returning exact solutions when they are known in analytic form as in the case of an elliptical cross-section.
5.2. Rectangular cross-section A rectangular cross-section of sides a and b as depicted in Fig. 3 changes in the transformed domain in another rectangle of different length sides namely Ra and b. For this section results obtained by LEM have been compared with existing results provided by a series expansion similar to that known for the isotropic case [27,28]. In particular, Lekhnitskii [1] suggests three factors namely b, k1, k2 in form of numerical series depending only on the ratio of the elastic moduli R and the ratio of the sides namely the factor d = Rab 1. By using these factors very simple formulas for the torsional rigidity and the maximum shearing stress at points A and B (see Fig. 3) can be expressed as follows: Mz =y ¼ Gzx ab3 b;
Series coefficients a 0 Lagrange multiplier l Unitary twist angle y Shear stresses tzx(x, y) Shear stresses tzy(x, y) Prandtl function c(x, y) Warping function o(x, y) Torsional rigidity Mz/y
b
A
5. Numerical applications
2 tðAÞ zx ¼ Mz k1 =ab ;
2 tðBÞ zy ¼ RMz k2 =ab
ð38Þ
and considered for a comparison. LEM applied to a rectangular cross-section returns all the coefficients of the Laurent series equal zero with exception of bk 8k ¼ 1; 3; 5; . . .. Numerical applications have been conducted for two different values of the factor d (d = 4,1) and for both cases the order of truncation has been selected sup = 8. First let us consider the case where the lengths of the sides and the orthotropic ratio are selected a= 4, b =3 and R= 3, respectively (d = 4). For the above orthotropic data LEM results are compared with the CVBEM results available in [24] and those obtained via Lekhnitskii series expansion in terms of torsional rigidity and shear stresses at the middles of the sides namely point A and B (Table 3). For the same orthotropic data the coefficients b, k1 and k2 have been computed for increasing order of truncation sup in the LEM and then compared with values provided by a boundary method [17] and the tabular data in [1] (Table 4). In Figs. 4 and 5 stresses due to torsion and contour lines of warping functions are plotted, respectively, for the isotropic (R= 1) and orthotropic (R= 3) cases under a unit torsion moment.
ARTICLE IN PRESS R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
49
It is worth to note that while the shearing stress in an isotropic bar is maximum at the middles of the longer sides, in an orthotropic bar, the points of maximum shearing stress may be at the middles either of the longer or of the shorter sides, depending on d = Rab 1. Let, for example, the lengths of the sides be a =4 and b= 2 (a = 2b) and let the orthotropic ratio be R= 2, that is d= 4. For
these selected values the maximum stress occurs in the point A at the middle of the shorter side b. As expected in the isotropic (R= 1) rectangular cross-section with the same geometry the maximum shear stress occurs at the middle of the longer side a. This difference in the stress distribution is evident in the stress maps depicted in Fig. 6.
Table 4 Comparison of BEM and LEM results for an orthotropic rectangular cross-section (d =Rab 1 = 4).
Table 5 Comparison of LEM and Lekhnitskii results for an orthotropic rectangular crosssection (d = Rab 1 = 1).
BEM [17]
LEM
Nel
sup
10 8 20 12 40 16 Lekhnitskii [1]
k1
b
k2
LEM
k1
b
bk a 0
BEM
LEM
BEM
LEM
BEM
LEM
Sup
bk a 0
4 6 8
0.270 0.280 0.281 0.281
0.282 0.281 0.281
3.702 3.561 3.552 3.550
3.529 3.518 3.557
2.556 2.694 2.686 2.644
2.826 2.647 2.610
8 2 12 3 16 4 Lekhnitskii [1]
k2
LEM
Err (%)
LEM
Err (%)
LEM
Err (%)
0.140 0.141 0.141 0.141
0.4 0.3 0.3
4.754 4.827 4.793 4.803
1 0.5 0.2
4.754 4.827 4.793 4.803
1 0.5 0.2
Fig. 4. Stresses due to torsion: (a) isotropic case R= 1 and (b) orthotropic case R = 3.
Fig. 5. Contour lines of warping function: (a) isotropic case R =1 and (b) orthotropic case R= 3.
Fig. 6. Stress map (a) isotropic case R= 1; (b) orthotropic case R =2.
ARTICLE IN PRESS 50
R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
As well for d = 1 the rectangular section changes in a square cross-section of length side aR = b whose all the Laurent coefficients are zero with exception of b3 ; b7 ; b11 ; . . . and a good accuracy is reached with a lower number of series coefficients involved in the solution. In Table 5 a comparison between LEM results and the tabular Lekhnitskii benchmark solution in terms of b, k1 and k2 is reported.
5.3. Annular cross-section The annular cross-section depicted in Fig. 7 in the transformed domain changes in a hollow elliptical cross-section whose semiaxes are ReR and Re for the external contour and RiR and Ri for the internal contour, respectively. For the annular cross-section under torque Mz the LEM (for x0 ¼ Z0 ¼ 0) returns only one coefficient of the Laurent series different from zero, namely b1, whatever the orders of truncation
Re
x
Ri
inf and sup are. It has to be emphasized that this is an exact solution since sT n ¼ 0 proves to be locally fulfilled in the external and internal contours whatever the value of the orthotropic ratio R is. Table 6 summarizes results in terms of unknown series coefficients, Lagrange multiplier l, shear stress field s, Prandtl c(x, y) and warping o(x, y) functions, respectively, where Re and Ri are the internal and external radius, respectively, and IP ¼ pðR2e R2i Þ=2. As for the elliptic cross-section, the stress distribution is independent of the value of the orthotropic ratio R and clash with the stress field provided for the isotropic case. Otherwise a closer inspection in the expression of the warping function o(x, y) (see Table 6) shows its dependence on the value of the orthotropic ratio R: by setting R= 1 (isotropic case) the cross-section remains plane and the well-known result for isotropic circular and annular cross-sections is restored. In Fig. 8 warping function and its contour lines are depicted for the case Gzx = 4Gzy. It is worth to note that although the annular cross-section is a multiply connected cross-section no polar singularities are involved in the solution and the regular part of the Laurent series is sufficient in order to provide a correct shear stress distribution in presence of a torsion moment as for a simply connected circular cross-section. Moreover by assuming Ri = 0 the coefficient b1 in Table 6 exactly coincides with the same coefficient in Table 1 by letting a = b=Re in the elliptical cross-section.
y
R Fig. 7. Annular cross-section.
Ri i x
Table 6 Results for shear stresses for an annular cross-section. Series coefficients a 0 Lagrange multiplier l Unitary twist angle y Shear stresses tzx(x, y) Shear stresses tzy(x, y) Prandtl function c(x, y) Warping function o(x, y) Torsional rigidity Mz/y
b1 = (R2 1)Mz/2IP(Re2 Ri2)R 0 (R2 + 1)Mz/2GzxIPR2 Mzy/IP(Re2 +Ri2) Mzx/IP(Re2 + Ri2) Mz(x2 + y2)/2IP(Re2 + Ri2) + const (R2 1)xy/(R2 + 1) 2pGzxIP(Re2 + Ri2)R2/(1 + R2)
Re 2
y
Fig. 9. Annular cross-section.
Fig. 8. Annular cross-section: (a) warping function and (b) its contour lines ðMz ¼ 1; Re ¼ 2; Ri ¼ Re =2; R ¼ 1=2Þ.
ARTICLE IN PRESS R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
5.4. Hollow circular cross-section If the centre of the hollow does not coincide with the centre of the external contour as depicted in Fig. 9 the pole of the principal part of the Laurent series moves from the point of coordinates x0 ¼ Z0 ¼ 0 to the point with coordinates x0 ¼ Re R=2 and Z0 ¼ 0
51
that is the centre of the elliptical hollow in the transformed domain. The numerical application studies the orthotropic case of Gzy =2Gzx with Re =2 and Ri =Re/4. Due to the symmetry LEM returns all the coefficients ak = 08k and bk a08k; in this case the solution is not exact because the condition sT n ¼ 0 is not locally
Fig. 10. Hollow circular cross-section: (a) Prandtl function and (b) its contour line ðMz ¼ 1; Re ¼ 2; Ri ¼ Re =4; R ¼
pffiffiffi 2Þ.
Fig. 11. Hollow circular cross-section ðMz ¼ 1; Re ¼ 2; Ri ¼ Re =4Þ: stress vector contour lines (a) isotropic R =1 and (b) orthotropic case R ¼
Fig. 12. Hollow circular cross-section ðMz ¼ 1; Re ¼ 2; Ri ¼ Re =4Þ: warping contour lines (a) isotropic R= 1 and (b) orthotropic case R ¼
pffiffiffi 2.
pffiffiffi 2.
ARTICLE IN PRESS 52
R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
fulfilled in all points of the contour. It is worth to note that the solution for a hollow elliptical cross-section, that is the shape assumed by the cross-section in the transformed domain, requires a computational effort greater than that required by the solution for a hollow circular cross-section. It follows that a greater number of polynomials is required to reach a good accuracy in the solution of the orthotropic case that is the lower order of truncation has been brought from inf =3 for the isotropic case to inf = 9 in the orthotropic case. In both cases the higher order of the series has been selected sup= 3. In Fig. 10 the Prandtl function for the case under exam is depicted with its contour lines. In Figs. 11 and 12 a comparison with the isotropic solution (R= 1) is reported in terms of stress vector and warping contour lines, respectively.
y take the following expressions a2 ¼ 0:07318Mz R2 =b and y ¼ 0:747597Mz R=b4 Gzx . Two different combinations of mechanical and geometrical parameters have been considered both satisfying the condition R ¼ a1 b ¼ a 1 b. In Figs. 14 and 15 the modulus of the stress vector functions and their contour lines are depicted for R= 0.5 and a = 2 and for R=2 and a =0.5, respectively, once b has been selected equal 1. In Figs. 16 and 17 Prandtl functions and their contour lines are depicted for the two data combinations under exam. For both cases LEM provides exact solutions as follows from the observation that the Prandtl functions take a constant value on the contour and their external contour lines move tangential along the boundary of the cross-section.
5.5. Elliptical cross-section with elliptical groove Let us consider an elliptical cross-section with an elliptical groove as illustrated in Fig. 13a. The coordinate transformation leads the original cross-section in a cross-section whose external contour is an ellipse of semi-axes aR and b, respectively; likewise the elliptical groove changes in an ellipse of semi-axes aR and b. It is worth to note that for a suitable choice of mechanical and geometrical data that is R ¼ a1 b ¼ a 1 b the cross-section in the transformed domain becomes a circular cross-section with a circular groove (see Fig. 13b) whose exact analytical solution exists [28]. In the numerical application the latter case is pursued. Also in this case the principal part of the Laurent series is essential in order to get a correct torsion solution and the pole is selected in the centre of the groove. Whatever the orders of series truncation are the only series coefficient different from zero is a 2 as expected. By selecting the qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi angle and g ¼ p=10 with b ¼ b sin g2 þð1 cos gÞ2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a ¼ a sin g2 þ ð1 cos gÞ2 the coefficient a 2 and the twist angle
aR = b
a
b x
b
6. Conclusions The extension of the novel numerical technique, labelled line element-less method (LEM), is presented in order to provide approximate solutions of the De Saint Venant torsion problem for orthotropic beams with simply and multiply connected crosssection. A suitable transformation of coordinates allows to restore the torsion problem as that studied in the isotropic case. A complex potential function analytic in all the domain whose real and imaginary parts are related to the shear stress components and to the orthotropic ratio has been introduced and expanded in the double-ended Laurent series involving harmonic polynomials. The potential function returns the domain governing equations, namely equilibrium and compatibility, while in order to satisfy the missing condition on the boundary an element-free weakform procedure has been proposed. By imposing that the square of the net flux of the shear stress across the border is minimum with respect to the series coefficients under the static equivalence condition the unknown series coefficients are thus provided. Numerical implementation of the LEM results in system of linear algebraic equations involving symmetric and positive-definite matrices. All the integrals are transferred into the boundary without requiring any discretization neither in the domain nor in the contour. The technique provides the complete shear stress field as shown by various numerical applications in order to assess the efficiency and the accuracy of the method to handle shear stress problems in presence of orthotropic material.
γγ
b
b Acknowledgement
y
a
aR = b
Fig. 13. (a) Elliptical cross-section with elliptical groove and (b) circular crosssection with circular groove.
The author wishes to thank Professor Mario Di Paola for stimulating discussions and suggestions on the subject of the paper.
Fig. 14. Elliptical cross-section with elliptical groove (R =0.5, b = 1, a= 2): (a) modulus of stress vector and (b) its contour lines.
ARTICLE IN PRESS R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
53
Fig. 15. Elliptical cross-section with elliptical groove (R = 2, b= 1, a= 0.5): (a) modulus of stress vector and (b) its contour lines.
Fig. 16. Elliptical cross-section with elliptical groove (R= 0.5, b =1, a= 2): (a) Prandtl function and (b) its contour lines.
where I P represents the polar inertia moment in the transformed domain O respect to the axes x and Z
Appendix A The static equivalence of the shear stresses as expressed in Eq. (31b) in terms of harmonic polynomials is fulfilled requiring that:
O
ðA:1Þ
O
or Z sup 1 X ak 2 R k ¼ inf k a 1
2
ðx þ Z Þ dO ¼
I
x3 3
O
nx þ
Z3 3
! ny
dC
ðA:3Þ
and
Z sup 1 X 1 þ 2 bk ðxPk þ ZQk Þ dO þ Gzx yI P ¼ Mz R R k ¼ inf k a 1
2
C
Z sup 1 X ak ðxQk ZPk Þ dO 2 R k ¼ inf k a 1
IP ¼
Z
Z sup 1 X div u k dO þ 2 bk R k ¼ inf O k a 1
1 div v k dO þ Gzx yI P ¼ Mz R O ðA:2Þ
ZPk þ 1 ðx x0 ; Z Z Þ=ðkþ 1Þ 0 u k ¼ xPk þ 1 ðx x0 ; Z Z Þ=ðk þ 1Þ ; 0 ZQk þ 1 ðx x0 ; Z Z Þ=ðk þ 1Þ 0 vk ¼ xQk þ 1 ðx x0 ; Z Z0 Þ=ðkþ 1Þ
ðA:4a; bÞ
Then, according to the Green lemma, the integrals in Eq. (A.2) may be performed as line boundary integrals as
ARTICLE IN PRESS 54
R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
Fig. 17. Elliptical cross-section with elliptical groove (R= 2, b = 1, a= 0.5): (a) Prandtl function and (b) its contour lines.
follows: sup X
1 a R2 k ¼ inf k
I
k a 1
¼
u Tk n dC þ
sup X
1 b R2 k ¼ inf k k a 1
C
sup X
sup X
ak c k þ
k ¼ inf k a 1
bk d k þ
k ¼ inf k a 1
I
v Tk n dC þ
form in the transformed domain: sup 1 X ak b Qk þ 1 ðx; ZÞ þ k Pk þ 1 ðx; ZÞ cðx; ZÞ ¼ kþ 1 R k ¼ inf kþ 1
Gzx yI p R
k a 1
C
Gzx yI p ¼ Mz R
ðA:5Þ
where ck ¼
1 R2
I
u Tk n dC ;
dk ¼
C
1 R2
I
v Tk n dC
ðA:6a; bÞ References
C
By collecting c k and d k in the vector h as T
h ¼ jc inf c 0 c sup d inf d 0 d sup j
ðA:7Þ
Eq. (A.5) can be expressed as follows: T
h w þ ðGzx =RÞyI p ¼ Mz
1 2 ðB:2Þ yGzx ðx þ Z2 Þ 2 Once the Prandtl function is known in the transformed domain, the simple inverse variable transformation (10a,b) allows to provide the Prandtl function for the cross-section under exam in the orthotropic case.
ðA:8Þ
Appendix B Bearing in mind the relationship between shear stresses and warping function in Eqs. (12a,b) and their series expansion in Eqs. (26a,b) and by taking advantage of the harmonic polynomial rules the warping function can be expressed in the transformed domain in terms of harmonic polynomials as follows: sup X 1 ak b Pk þ 1 ðx; ZÞ k Qk þ 1 ðx; ZÞ oðx; ZÞ ¼ ðB:1Þ 2 kþ1 yGzx R k ¼ inf k þ1 k a 1
Once the warping function is provided in the transformed domain, the simple inverse variable transformation (10a,b) allows to provide the warping function for the cross-section under exam in the orthotropic case. Similarly bearing in mind the relationship between shear stresses and Prandtl function in Eqs. (16a,b) and their series expansion in Eqs. (26a,b) the Prandtl function takes the following
[1] Lekhnitskii SG. Theory of elasticity of an anisotropic body. Moscow: Mir Publishers; 1981. [2] Tolf G. St. Venant bending of an orthotropic beam. Composite Structures 1985;4:1–14. [3] Herrmann RL. Elastic torsional analysis of irregular shapes. ASCE Journal of Engineering Mechanics Division 1965;91(EM6):11–9. [4] Mason WE, Herrmann RL. Elastic shear analysis of general prismatic shaped beams. ASCE Journal of Engineering Mechanics Division 1968;94(EM4): 965–983. [5] Kosmatka JB, Dong SB. Saint-Venant solutions for prismatic anisotropic beams. International Journal of Solids and Structures 1991;28(7):917–38. [6] Kosmatka JB. Flexure-torsion behavior of prismatic beams, part I: section properties via power series. AIAA Journal 1993;31(1):170–9. [7] Jawson MA, Ponter AR. An integral equation solution of the torsion problem. Proceedings of the Royal Society of London 1963;275(A):237–46. [8] Chou SI, Mohr JA. A boundary integral method for torsion of composite shafts. Res Mechanica 1990;29:41–56. [9] Katsikadelis JT, Sapountzakis EJ. Torsion of composite bars by boundary element method. ASCE Journal of the Engineering Mechanics Division 1985;111(9):1197–210. [10] Sapountzakis EJ, Mokos VG. Nonuniform torsion of composite bars by boundary element method. Journal of Engineering Mechanics, ASCE 2001;127(9):945–53. [11] Sapountzakis EJ. Nonuniform torsion of multi-material composite bars by the boundary element method. Computers and Structures 2001;79(32):2805–16. [12] Sapountzakis EJ, Mokos VG. Warping shear stresses in nonuniform torsion of composite bars by BEM. Computer Methods in Applied Mechanics and Engineering 2003;192(39–40):4337–53. [13] Sapountzakis EJ, Mokos VG. Nonuniform torsion of composite bars of variable thickness by BEM. International Journal of Solids and Structures 2004;41(7):1753–71. [14] Friedman Z, Kosmatka JB. Torsion and flexure of a prismatic isotropic beam using the boundary element method. Computers & Structures 2000;74: 479–494.
ARTICLE IN PRESS R. Santoro / International Journal of Mechanical Sciences 52 (2010) 43–55
[15] Petrolo S, Casciaro R. 3D beam element based on Saint Vena nt’s rod theory. Computers & Structures 2004;82:2471–81. [16] Gracia L, Doblare M. Shape optimization of elastic orthotropic shafts under torsion by using boundary elements. Composite Structures 1988;30(6): 1281–1291. [17] Gaspari D, Aristodemo M. Torsion and flexure analysis of orthotropic beams by a boundary element model. Engineering Analysis with Boundary Elements 2005;29:850–8. [18] Hromadka II TV, Guymon GL. Complex polynomial approximation of the Laplace equation. Journal of Hydraulic Engineering ASCE 1984;110(3):329–39. [19] Poler AC, Bohannon AW, Schowalter SJ, Hromadka II TV. Using the complex polynomial method with mathematica to model problems involving the Laplace and Poisson equations. Published online in Wiley InterScience /www.interscience.wiley.com) doi:10.1002/num.20365S. [20] Hromadka II TV, Lai Ch. The complex variable boundary element method in engineering analysis. New York: Springer; 1987. [21] Hromadka II TV, Whitley RJ. Advances in the complex variable boundary element method. London: Springer-Verlag; 1998.
55
[22] Aleynikov SM, Stromov AV. Comparison of complex methods for numerical solutions of boundary problems of the Laplace equation. Engineering Analysis with Boundary Elements 2004;28:615–22. [23] Hromadka II TV, Pardoeon G. Application of the CVBEM to nonuniform St. Venant torsion. Computer Methods in Applied Mechanics and Engineering 1985;53:149–61. [24] Dumir PC, Kumar R. Complex variable boundary element method for torsion of anisotropic bars. Applied Mathematical Modelling 1993;17: 80–88. [25] Di Paola M, Pirrotta A, Santoro R. Line element-less method (LEM) for beam torsion solution (TRULY NO-MESH METHOD). Acta Mechanica 2008;195(1–4):349–64. [26] Sokolnikoff IS. Mathematical theory of elasticity. New York, NY: McGrawHill; 1956. [27] Muskhelishvili NI. Some basic problems of the mathematical theory of elasticity. Groningen, The Netherlands: P. Noordhoff Ltd.; 1953. [28] Timoshenko SP, Goodier JN. Theory of elasticity. New York, NY: McGraw-Hill; 1970.