Saint-Venant torsion of orthotropic bars with inhomogeneous rectangular cross section

Composite Structures 92 (2010) 1449–1457

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Saint-Venant torsion of orthotropic bars with inhomogeneous rectangular cross section Xu Rongqiao a,*, He Jiansheng a, Chen Weiqiu b a b

Department of Civil Engineering, Zhejiang University, Zijingang Campus, Hangzhou 310058, China Department of Engineering Mechanics, Zhejiang University, Yuquan Campus, Hangzhou 310027, China

a r t i c l e

i n f o

Article history: Available online 5 November 2009 Keywords: Saint-Venant torsion Inhomogeneous rectangle Orthotropy Exact solution

a b s t r a c t This work presents an exact solution of the Saint-Venant torsion of a straight bar with an orthotropic and inhomogeneous rectangular cross section whose material properties obey the exponential law in one direction. An approximate solution is also obtained for the material properties being arbitrarily distributed in one direction by using a layerwise model, in which the inhomogeneous rectangle is simulated by a composite rectangle composed of multiple rigidly connected homogeneous rectangular regions. The warping function, stresses and torsional rigidity are analytically expressed in terms of Fourier series for both solutions, in which the hyperbolic functions are not directly employed to avoid numerical difficulties. Some numerical examples are finally presented to verify the proposed method and the parametrical study is also performed. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction The analysis of the torsion of a straight bar could be traced back to Coulomb, who investigated the torsion of a circular cylinder and invented a torsion balance to accurately measure extremely minute forces. Due to the torsion balance, he developed the well-known Coulomb’s law of the theory of attraction and repulsion between bodies of the electrical charge. Afterward, Navier investigated the torsion of a straight bar with a non-circular section using the torsion theory proposed by Coulomb but obtained self-contradictory results. Whereas, Saint-Venant successfully derived the governing equation of the torsion of the straight bar with an arbitrary section through a semi-inverse method. Prandtl then introduced the stress function of the Saint-Venant torsion and the method of membrane analogy [1]. Since these pioneering works established the theory of torsion and solved many problems in engineering application, the torsion of a straight bar became a classical problem in the theory of elasticity [2], which was also presented as a numerical example in a seminal paper about the finite element method by Courant [3]. Some analytical solutions of the homogeneous section with various shapes are available in the literatures [4,5]. The torsion of composite bar has attracted many researchers’ attention in the development and application of composite materials. Muskhelishvilli [6] presented the governing equation and boundary condition of the torsion of composite bars and its solution in Fourier series for composite section with two sub-rectangles. This solution was * Corresponding author. Tel.: +86 571 88208478; fax: +86 571 88208685. E-mail address: [email protected] (X. Rongqiao). 0263-8223/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2009.10.042

extended later for multiple rectangular composite section by Booker and Kitipornchai [7]. Kuo and Conway [8–11] analyzed the torsion of the composite sections of various shapes. Packham and Shail [12] extended their work on two-phase fluid to the torsion of composite bars. Ripton [13] investigated the torsional rigidity of composite section reinforced by fibers. Chen et al. [14] also analyzed exactly the torsion of composite bars. Apart from these analytical methods, numerical methods have also been employed to solve the torsion of the straight bars. Ely and Zienkiewicz [15] firstly solved the Poisson’s equation of the Prandtl’s stress function using finite difference method and investigated the rectangular section with and without holes. Herrmann [16] utilized the finite element method to calculate the warping function of the torsion of irregular sectional shapes. The boundary element method was applied to solve the boundary integral equation of the warping function of the torsion in Refs. [17–20]. Recently functionally graded materials (FGMs), whose properties continuously vary with spatial coordinates, have been developed to overcome the problems associated with interfaces in traditional composite materials due to the abrupt change of the materials properties [21]. Although Ely and Zienkiewicz [15] and Plunkett [22] presented the governing equation of the torsion of inhomogeneous material before the introduction of the conception of FMGs, it was paid little attention as there is no engineering significance at that time. Once the FGMs were fabricated and applied in engineering practice, Rooney and Ferrari [23,24] and Horgan and Chan [25] resumed the research on the torsion of FGM bars. More recently, Tarn and Chang [26] obtained the exact solution of the torsion of orthotropic inhomogeneous cylinders and also analyzed

1450

X. Rongqiao et al. / Composite Structures 92 (2010) 1449–1457

the end effect. All these reported works only focus on the circular cylinders, and the investigation of non-circular shapes is still scarce. This paper considers a straight bar with a rectangular section, whose material is orthotropic and its properties vary along one coordinate. An exact solution is obtained for the case in which the material properties vary in an exponential law. For the case of arbitrary variation of the material properties, a layerwise model is adopted and the FGM section is simulated by a composite section with multiple sub-rectangular regions, which are rigidly connected at the fictitious interfaces. The material properties are assumed constant in each sub-rectangular region and the governing equation of the warping function is derived as well as the corresponding boundary conditions and compatible conditions at the interfaces. An analytical solution in series form is obtained through the separation of variables method. For the numerical convenience, the hyperbolic functions are avoided in the solutions for the sake of numerical stability. Numerical examples and parametric studies are also carried out to validate the solutions.

3. Exact solution If the variation of the shear moduli obeys

 x Gx ¼ G0x exp k ; a

 x Gy ¼ a2 G0x exp k a

ð6Þ

in which the parameter k is called the gradient factor. Substituting it into Eq. (3) yields

  2 @2/ k @/ 2@ / þ a þ  y ¼0 @x2 @y2 a @x

ð7Þ

Another function w is introduced as [6]

w ¼ / þ xy

ð8Þ

and substituting it into Eq. (7) gives

  2 @2w k @w 2@ w þ a þ  2y ¼0 @x2 @y2 a @x

ð9Þ

while the boundary conditions (4) and (5) become 2. Description of the problem Consider a straight bar with the rectangular cross section with height 2b and width a as shown in Fig. 1. A vertical side of the rectangle is taken as y-axis while x-axis is parallel to the horizontal direction. The origin of the coordinates is located at the mid point of the vertical side, and the material is assumed orthotropic and inhomogeneous. The shear moduli in xz and yz planes are denoted by Gx and Gy , respectively, both being functions of x only. According to Saint-Venant torsion theory, the displacements can be expressed as

u ¼ hzy;

v ¼ hzx;

w ¼ h/ðx; yÞ

ð1Þ

in which h is the twist angle of unit length, and / is the warping function. Substitution of Eq. (1) into stress displacement relationships, the two nonzero stress components are

syz

  @/ ¼ Gy h þx ; @y

sxz

  @/ ¼ Gx h y @x

ð2Þ

Since there are only two nonzero stress components and they are independent of z coordinate, the equilibrium equations in the x and y directions are satisfied immediately and the equilibrium equation in the z direction becomes

      @ @/ @ @/ Gx y þ Gy þx ¼0 @x @x @y @y

ð3Þ

@w ¼ 2y at x ¼ 0; a @x @w ¼ 0 at y ¼ b @y

@/  y ¼ 0 at x ¼ 0; a @x

ð4Þ

ð11Þ

To satisfy boundary condition (11), we assume [6] 1 X

w¼

fn ðxÞ sinðnn yÞ

ð12Þ

n¼1;3;5;...

in which fn ðxÞ is a series of undetermined functions and

nn ¼

np 2b

ðn ¼ 1; 3; 5; . . .Þ

ð13Þ

The term 2y in Eq. (9) can also be expressed by Fourier series, namely, 1 X

2y ¼

C n sinðnn yÞ

ð14Þ

n¼1;3;5;...

where the coefficient C n is given by n1 2

Cn ¼

16b ð1Þ p2 n2

ðn ¼ 1; 3; 5; . . .Þ

ð15Þ

Substitution of Eq. (12) into Eq. (9) yields 2

The traction-free condition on the lateral surface of the bar requires that

ð10Þ

d fn ðxÞ 2

dx

þ

k dfn ðxÞ k  a2 n2n fn ðxÞ ¼ C n a dx a

ðn ¼ 1; 3; 5; . . .Þ

Thus the function fn ðxÞ can be rendered as

fn ðxÞ ¼ An exp½b1n ðx  aÞ þ Bn expðb2n xÞ 

and

@/ þ x ¼ 0 at y ¼ b @y

ð16Þ

kC n aa2 n2n

ðn ¼ 1; 3; 5; . . .Þ ð17Þ

ð5Þ

where An and Bn are unknown parameters and can be determined from the boundary condition (11) and

y b1n

R

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k2 k ¼ þ a2 n2n  ; 4a2 2a

b2n

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k2 k ¼ þ a2 n2n þ 4a2 2a

ðn ¼ 1; 3; 5; . . .Þ ð18Þ

2b

x

a Fig. 1. The inhomogeneous rectangular section.

As shown in Eq. (17), the hyperbolic functions are not used intendedly to avoid the numerical difficulties. Alternatively, two exponential functions are adopted to guarantee that the independent variables of the exponential functions are always negative, assuring numerical convergence at large n.

1451

X. Rongqiao et al. / Composite Structures 92 (2010) 1449–1457

If the gradient factor k ¼ 0, the material becomes isotropic and Eqs. (16) and (17) reduce to the forms in the isotropic case. Substitution of Eqs. (17) and (14) into Eq. (11) yields

½1  expðab2n ÞC n ; b1n f1  exp½aðb1n þ b2n Þg ½1  expðab1n ÞC n Bn ¼ ðn ¼ 1; 3; 5; . . .Þ b2n fexp½aðb1n þ b2n Þ  1g

An ¼

ð19Þ

ui ¼ hzy;

Thus, the function w is determined and the warping function / is then obtained from Eq. (8). The torsional rigidity can be calculated by

     @/ @/ Gy þ x x  Gx  y y dxdy @y @x R  Z Z   x @w @w 0 2y2  y ¼ Gx þ a2 x exp k dxdy @x @y a R

Z Z 

syz;i ¼ Gyi h ð20Þ

Substituting of Eqs. (12) and (17) into Eq. (20) yields

  4 3 expðkÞ  1 D ¼ G0x ab  Dx þ 2a2 Dy 3 k

ð21Þ

ð24Þ

  @/i þx ; @y



sxz;i ¼ Gxi h

@/i y @x

 ð25Þ

since the shear moduli are assumed constant, the warping function /i satisfies 2 @ 2 /i 2 @ /i þ a ¼0 i @x2 @y2

1 X

n1 2

a2

ð1Þ

n¼1;3;5;...

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

  8b expðkÞ  expðab1n Þ expðk  ab2n Þ  1  2 2 An b1n þ Bn b2n b1n þ k=a b2n  k=a n p 2

ð22Þ and

Dy ¼

wi ¼ h/i ðx; yÞ

ð26Þ

Table 1 The torsional rigidity factor k for orthotropic and homogeneous rectangular section.

where

Dx ¼

v i ¼ hzx;

in which the subscription i means the corresponding physical quantity associated with the sub-rectangle Ri . The two nonzero stress components in this sub-rectangle are

(

  a2 An expðkÞ  expðab1n Þ expðkÞ  ab1n þ k ab1n þ k n¼1;3;5;...   2 a Bn expðk  ab2n Þ  1 exp ðk  ab2n Þ  þ k  ab2n k  ab2n   ) 2 4ab C n 1 1 ð23Þ  2 2 2 1 expðkÞ þ k k anp 1 X

n1 2

ð1Þ

2 3 4 5 6 7 8 9 10

4. Layerwise model and its solution * **

It is difficult to solve Eq. (3) exactly if the shear moduli Gx ðxÞ and Gy ðxÞ vary arbitrarily with respect to the coordinate x. For this purpose, the original rectangular section is simulated by a composite section of multiple sub-rectangles as shown in Fig. 2. The number of the sub-rectangles is denoted by m and the material of the each sub-rectangle is assumed homogeneous. Namely, In sub-rectangle i, the material properties are assumed constant and equal to those at the center of the sub-rectangle, denoted by Gxi and Gyi . Let x0 ¼ 0 and xm ¼ a are the left and right boundaries, respectively, of the original rectangular section while x ¼ xi ði ¼ 1; 2; . . . ; m  1Þ denote the interfaces between the sub-rectangles. It can be

y

R1

Rm −1

R2

Rm

s ¼ 0:001*

s ¼ 0:1

s ¼ 0:2

s ¼ 0:4

s ¼ 0:6

s ¼ 0:8

s¼1

0.0333 0.0667 0.1000 0.1333 0.1666 0.1999 0.2332 0.2665 0.2998 0.3331

0.0327 0.0648 0.0965 0.1280 0.1592 0.1902 0.2210 0.2516 0.2821 0.3123 (0.312)** 0.6072 0.8908 1.1653 1.4318 1.6912 1.9442 2.1913 2.4328 2.6690

0.0320 0.0629 0.0931 0.1227 0.1518 0.1805 0.2087 0.2366 0.2641 0.2913 (0.291) 0.5478 0.7817 0.9975 1.1977 1.3845 1.5592 1.7232 1.8775 2.0229

0.0307 0.0591 0.0862 0.1121 0.1370 0.1609 0.1841 0.2066 0.2283 0.2494 (0.249) 0.4308 0.5726 0.6869 0.7813 0.8607 0.9284 0.9870 1.0382 1.0834

0.0293 0.0554 0.0793 0.1015 0.1222 0.1416 0.1598 0.1770 0.1932 0.2086

0.0280 0.0516 0.0724 0.0910 0.1077 0.1229 0.1367 0.1493 0.1610 0.1717

0.3272 0.4055 0.4614 0.5036 0.5366 0.5633 0.5854 0.6040 0.6199

0.2468 0.2899 0.3183 0.3387 0.3540 0.3661 0.3760 0.3842 0.3912

0.0267 0.0479 0.0657 0.0809 0.0941 0.1057 0.1159 0.1250 0.1332 0.1406 (0.1406) 0.1882 0.2131 0.2287 0.2395 0.2476 0.2540 0.2591 0.2633 0.2669

0.6661 0.9989 1.3317 1.6643 1.9969 2.3294 2.6619 2.9943 3.3266

s denotes the width-to-height ratio of the cross section, i.e., s ¼ a=ð2bÞ. The data in parenthesis are obtained by Timoshenko and Goodier [2].

0.68

The torsion rigidity factor k =D/(2ba3G0x)

D¼

expected that the behavior of the composite section will converge to the original section with the increase of the number of the subrectangles m. Muskhelishvilli [6] has given the solution of the torsion of a composite section with two isotropic rectangular regions and the present paper extends his work to a composite section with multiple orthotropic rectangular regions. We also assume the displacements in a sub-rectangle Ri to be

0.67

layer model exact solution

0.66

0.65

2

Gx = Gx0 exp(2 x a ) G y = 4Gx

0.64

0.63

2

x 0.62

x0

x1

x2

xm −1

xm

Fig. 2. The layerwise model of inhomogeneous section.

0

10

20

30 40 50 60 70 The number of the sub-rectangle (m)

80

90

100

Fig. 3. The convergence of the torsional rigidity factor of the layerwise model.

1452

X. Rongqiao et al. / Composite Structures 92 (2010) 1449–1457 Table 2 Non-dimensional torsion constant and maximum resultant shear stress of composite cross-section.

y

G1 =G2

G2

G1

1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0

a x

a/2

a/2

smax

It Present

Ref. [20]

Ref. [6]

Present

Ref. [20]

Ref. [6]

0.1406 0.1970 0.2394 0.2764 0.3105 0.3431 0.3746 0.4055 0.4360 0.4661

0.1405 0.1969 0.2394 0.2763 0.3105 0.3431 0.3746 0.4055 0.4360 0.4661

0.1406 0.1970 0.2395 0.2764 0.3105 0.3431 0.3746 0.4055 0.4360 0.4661

0.6755 1.2108 1.7110 2.1972 2.6764 3.1516 3.6243 4.0954 4.5652 5.0343

0.6751 1.2101 1.7102 2.1962 2.6752 3.1502 3.6227 4.0936 4.5633 5.0321

0.6753 1.210 1.711 2.197 2.676 3.151 3.623 4.094 4.564 5.033

Fig. 4. Composite square cross-section consisting of two isotropic materials.

@ 2 wi @ 2 wi þ a2i ¼0 2 @x @y2

in which a2i ¼ Gyi =Gxi . If the material is isotropic, i.e., Gxi ¼ Gyi and a2i ¼ 1, Eq. (26) becomes Laplacian. Since the lateral surface of the bar is traction-free, we have

@/1  y ¼ 0 at x ¼ 0 @x @/m  y ¼ 0 at x ¼ a @x @/i þ x ¼ 0 ði ¼ 1; 2; . . . mÞ at y ¼ b @y

@w1 ¼ 2y at x ¼ 0 @x @wm ¼ 2y at x ¼ a @x @wi ¼ 0 ði ¼ 1; 2; . . . ; mÞ @y

ð27Þ ð28Þ ð29Þ

at y ¼ b

@wi @w  Gx;iþ1 iþ1 ¼ 2ðGxi  Gx;iþ1 Þy; @x @x ði ¼ 1; 2; . . . ; m  1Þ

ð35Þ

wiþ1 ¼ wi

Now the problem is reduced to find a warping function /i , which must satisfy simultaneously the governing Eq. (26) and the boundary conditions (27)–(29) with the continuity conditions (30) at the interfaces. In a similar manner, a new function wi is introduced

wi ¼

ð36Þ

1 X

fin ðxÞ sinðnn yÞ ði ¼ 1; 2; . . . ; mÞ

where fin ðxÞ denotes a series of functions with respect to x and the parameter nn is defined in Eq. (13). Apparently, the function wi in Eq. (37) satisfies the boundary condition (35). Substitution of Eq. (37) into (32), we have

The original governing Eq. (26) and the boundary/continuity conditions (27)–(30) become

0.06 0.04 0.02 0 -0.02 -0.04 -0.06 0.5

y

0 -0.5

ð37Þ

n¼1;3;5;...

ð31Þ

φM

at x ¼ xi

Thus if we can find the function wi , which satisfies the governing Eq. (32) and the boundary conditions (33)–(35) with the continuity conditions (36) at the interfaces, the problem will be solved. Assume the function wi in the form [6]

ð30Þ

wi ¼ /i þ xy

ð34Þ

Gxi

    @/i @/iþ1 Gxi  y ¼ Gx;iþ1 y at x ¼ xi @x @x

ði ¼ 1; 2; . . . ; m  1Þ

ð33Þ

and

And at the interfaces between sub-rectangles, the compatibility of the displacement and equilibrium of shear stress in the normal direction of the interfaces yield

/i ¼ /iþ1 ;

ð32Þ

0

0.2

0.6

0.4

0.8

x

 M of the composite cross section. Fig. 5. The distribution of the dimensional warping function /

1

1453

X. Rongqiao et al. / Composite Structures 92 (2010) 1449–1457

τ ( aG2θ ) 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0.5

y

0 -0.5

0.2

0

0.6

0.4

1

0.8

x

Fig. 6. The distribution of the dimensionless constant shear stress s=ðaG2 hÞ.

Table 3 The torsional rigidity k of isotropic and inhomogeneous rectangular section   Gax ¼ 3G0x . k

s ¼ 0:001*

s ¼ 0:1

s ¼ 0:2

s ¼ 0:4

s ¼ 0:6

s ¼ 0:8

s¼1

0 2 4 6 8 10 15 20 25 30 35 40 50 100 1000

0.9994 0.5330 0.4917 0.4719 0.4566 0.4441 0.4210 0.4055 0.3945 0.3863 0.3801 0.3751 0.3676 0.3510 0.3333

0.9370 0.4996 0.4586 0.4393 0.4249 0.4133 0.3922 0.3781 0.3682 0.3608 0.3551 0.3505 0.3438 0.3286 0.3124

0.8739 0.4659 0.4252 0.4064 0.3929 0.3823 0.3631 0.3505 0.3416 0.3350 0.3298 0.3258 0.3197 0.3061 0.2914

0.7481 0.3985 0.3586 0.3409 0.3292 0.3204 0.3052 0.2954 0.2886 0.2835 0.2795 0.2764 0.2717 0.2610 0.2494

0.6258 0.3332 0.2948 0.2784 0.2685 0.2614 0.2498 0.2426 0.2376 0.2339 0.2310 0.2287 0.2253 0.2174 0.2087

0.5152 0.2742 0.2386 0.2240 0.2156 0.2100 0.2013 0.1961 0.1925 0.1899 0.1879 0.1863 0.1838 0.1782 0.1718

0.4217 0.2245 0.1926 0.1798 0.1728 0.1684 0.1618 0.1580 0.1555 0.1536 0.1522 0.1510 0.1493 0.1453 0.1406

Table 5 The torsional   rigidity factor k of isotropic and inhomogeneous rectangular section Gax ¼ 7G0x .

* *

s denotes the width-to-height ratio of the cross section, i.e., s ¼ a=ð2bÞ.

Table 4 The torsional   rigidity factor k of isotropic and inhomogeneous rectangular section Gax ¼ 5G0x . k

s ¼ 0:001*

s ¼ 0:1

s ¼ 0:2

s ¼ 0:4

s ¼ 0:6

s ¼ 0:8

s¼1

0 2 4 6 8 10 15 20 25 30 35 40 50 100 1000

1.6657 0.6758 0.5939 0.5652 0.5443 0.5266 0.4916 0.4664 0.4478 0.4335 0.4223 0.4133 0.3997 0.3682 0.3334

1.5616 0.6350 0.5538 0.5249 0.5047 0.4880 0.4559 0.4331 0.4163 0.4034 0.3933 0.3852 0.3728 0.3443 0.3126

1.4566 0.5939 0.5133 0.4842 0.4647 0.4491 0.4199 0.3995 0.3845 0.3730 0.3640 0.3568 0.3458 0.3202 0.2915

1.2468 0.5116 0.4325 0.4031 0.3851 0.3717 0.3482 0.3325 0.3211 0.3125 0.3057 0.3001 0.2918 0.2721 0.2495

1.0431 0.4315 0.3552 0.3265 0.3101 0.2988 0.2806 0.2691 0.2609 0.2547 0.2499 0.2459 0.2399 0.2256 0.2087

0.8587 0.3584 0.2875 0.2607 0.2462 0.2368 0.2228 0.2145 0.2087 0.2044 0.2011 0.1983 0.1942 0.1841 0.1718

0.7029 0.2959 0.2322 0.2081 0.1956 0.1879 0.1769 0.1708 0.1667 0.1637 0.1614 0.1595 0.1566 0.1495 0.1406

s denotes the width-to-height ratio of the cross section, i.e., s ¼ a=ð2bÞ.

s ¼ 0:001*

s ¼ 0:1

s ¼ 0:2

s ¼ 0:4

s ¼ 0:6

s ¼ 0:8

s¼1

0 2 4 6 8 10 15 20 25 30 35 40 50 100 1000

2.3320 0.7996 0.6705 0.6341 0.6106 0.5907 0.5496 0.5184 0.4944 0.4756 0.4606 0.4483 0.4295 0.3848 0.3335

2.1862 0.7530 0.6258 0.5881 0.5650 0.5460 0.5080 0.4798 0.4582 0.4413 0.4277 0.4167 0.3997 0.3594 0.3127

2.0392 0.7060 0.5806 0.5423 0.5190 0.5009 0.4661 0.4408 0.4216 0.4066 0.3946 0.3848 0.3697 0.3337 0.2916

1.7455 0.6122 0.4906 0.4505 0.4276 0.4112 0.3827 0.3632 0.3488 0.3376 0.3286 0.3213 0.3100 0.2826 0.2496

1.4603 0.5204 0.4044 0.3639 0.3421 0.3277 0.3049 0.2905 0.2802 0.2723 0.2660 0.2609 0.2529 0.2333 0.2088

1.2021 0.4357 0.3286 0.2901 0.2702 0.2577 0.2396 0.2291 0.2218 0.2164 0.2121 0.2086 0.2031 0.1896 0.1719

0.9840 0.3625 0.2666 0.2315 0.2138 0.2032 0.1887 0.1808 0.1757 0.1719 0.1689 0.1665 0.1628 0.1534 0.1407

s denotes the width-to-height ratio of the cross section, i.e., s ¼ a=ð2bÞ.

Table 6 The torsional   rigidity factor k of isotropic and inhomogeneous rectangular section Gax ¼ 9G0x .

* *

k

k

s ¼ 0:001*

s ¼ 0:1

s ¼ 0:2

s ¼ 0:4

s ¼ 0:6

s ¼ 0:8

s¼1

0 2 4 6 8 10 15 20 25 30 35 40 50 100 1000

2.9983 0.9147 0.7333 0.6883 0.6630 0.6422 0.5982 0.5633 0.5356 0.5135 0.4954 0.4805 0.4573 0.4007 0.3337

2.8109 0.8630 0.6853 0.6387 0.6129 0.5926 0.5516 0.5199 0.4950 0.4751 0.4589 0.4456 0.4247 0.3738 0.3128

2.6218 0.8109 0.6368 0.5886 0.5623 0.5425 0.5046 0.4761 0.4540 0.4364 0.4222 0.4103 0.3919 0.3467 0.2917

2.2443 0.7068 0.5402 0.4890 0.4619 0.4432 0.4111 0.3891 0.3725 0.3595 0.3489 0.3402 0.3265 0.2925 0.2497

1.8775 0.6046 0.4474 0.3952 0.3684 0.3511 0.3247 0.3083 0.2965 0.2874 0.2800 0.2740 0.2645 0.2405 0.2089

1.5456 0.5096 0.3655 0.3153 0.2902 0.2749 0.2533 0.2411 0.2327 0.2265 0.2215 0.2174 0.2110 0.1946 0.1719

1.2652 0.4266 0.2980 0.2521 0.2294 0.2161 0.1983 0.1891 0.1830 0.1787 0.1752 0.1724 0.1681 0.1569 0.1407

s denotes the width-to-height ratio of the cross section, i.e., s ¼ a=ð2bÞ.

Its solution can be obtained readily as 2

d fin ðxÞ 2

dx

 a2i n2n fin ðxÞ ¼ 0 ði ¼ 1; 2; . . . ; m; n ¼ 1; 3; 5; . . .Þ

ð38Þ

fin ðxÞ ¼ Ain exp½ai nn ðx  aÞ þ Bin expðai nn xÞ ði ¼ 1; 2; . . . ; m; n ¼ 1; 3; 5; . . .Þ

ð39Þ

1454

X. Rongqiao et al. / Composite Structures 92 (2010) 1449–1457

 of an isotropic square section. Fig. 7. The distribution of the shear stress factor s

 of an isotropic rectangular section ðs ¼ 0:5Þ. Fig. 8. The distribution of the shear stress factor s

where Ain and Bin are undetermined constants dependent on the boundary and continuity conditions. Here the hyperbolic functions are not adopted again in Eq. (39) to avoid numerical difficulty. The unknown constants Ain and Bin can be determined from Eqs. (33), (34) and (36). Noticing Eq. (14), we have

ai Gxi exp½ai nn ðxi  aÞAin  ai Gxi expðai nn xi ÞBin  aiþ1 Gx;iþ1  exp½aiþ1 nn ðxi  aÞAiþ1;n þ aiþ1 Gx;iþ1 expðaiþ1 nn xi ÞBiþ1;n ¼ ðGxi  Gx;iþ1 ÞC n =nn

ði ¼ 1; 2; . . . ; m  1; n ¼ 1; 3; 5; . . .Þ

ð42Þ

and

Ain exp½ai nn ðxi  aÞ þ Bin expðai nn xi Þ  Aiþ1;n exp½aiþ1 nn ðxi  aÞ

A1n expða1 nn aÞ  B1n ¼ C n =ða1 nn Þ ðn ¼ 1; 3; 5; . . .Þ

ð40Þ

 Biþ1;n expðaiþ1 nn xi Þ ¼ 0 ði ¼ 1; 2; . . . ; m  1; n ¼ 1; 3; 5; . . .Þ

Amn  Bmn expðam nn aÞ ¼ C n =ðam nn Þ ðn ¼ 1; 3; 5; . . .Þ

ð41Þ

ð43Þ

1455

X. Rongqiao et al. / Composite Structures 92 (2010) 1449–1457 Table 7 The torsional rigidity factor k of orthotropic and inhomogeneous square section   Gax ¼ 3G0x .

Table 10 The torsional rigidity factor k of orthotropic and inhomogeneous square section   Gax ¼ 9G0x .

k

a2 ¼ 0:1

a2 ¼ 0:3

a2 ¼ 0:5

a2 ¼ 0:7

a2 ¼ 3

a2 ¼ 5

a2 ¼ 7

k

a2 ¼ 0:1

a2 ¼ 0:3

a2 ¼ 0:5

a2 ¼ 0:7

a2 ¼ 3

a2 ¼ 5

a2 ¼ 7

0 2 4 6 8 10 15 20 25 30 35 40 50 100 1000

0.0801 0.0427 0.0386 0.0368 0.0356 0.0346 0.0329 0.0318 0.0311 0.0305 0.0301 0.0297 0.0292 0.0280 0.0267

0.1971 0.1050 0.0933 0.0883 0.0851 0.0829 0.0791 0.0768 0.0751 0.0739 0.0730 0.0722 0.0711 0.0686 0.0657

0.2823 0.1503 0.1317 0.1240 0.1195 0.1163 0.1114 0.1084 0.1063 0.1047 0.1035 0.1026 0.1011 0.0978 0.0941

0.3477 0.1851 0.1606 0.1506 0.1449 0.1411 0.1354 0.1319 0.1296 0.1279 0.1265 0.1255 0.1239 0.1202 0.1159

0.6392 0.3417 0.2846 0.2624 0.2511 0.2444 0.2357 0.2312 0.2285 0.2266 0.2251 0.2240 0.2223 0.2182 0.2131

0.7186 0.3856 0.3187 0.2925 0.2794 0.2717 0.2618 0.2571 0.2543 0.2524 0.2510 0.2499 0.2483 0.2445 0.2396

0.7619 0.4100 0.3379 0.3095 0.2952 0.2868 0.2762 0.2712 0.2683 0.2664 0.2650 0.2640 0.2624 0.2588 0.2540

0 2 4 6 8 10 15 20 25 30 35 40 50 100 1000

0.2402 0.0750 0.0581 0.0531 0.0504 0.0484 0.0450 0.0425 0.0406 0.0392 0.0379 0.0369 0.0354 0.0315 0.0267

0.5913 0.1893 0.1413 0.1256 0.1175 0.1122 0.1038 0.0985 0.0946 0.0916 0.0892 0.0872 0.0840 0.0761 0.0658

0.8470 0.2762 0.2009 0.1752 0.1621 0.1540 0.1422 0.1352 0.1303 0.1266 0.1236 0.1211 0.1173 0.1075 0.0942

1.0431 0.3454 0.2464 0.2118 0.1944 0.1840 0.1694 0.1613 0.1558 0.1517 0.1485 0.1458 0.1416 0.1309 0.1160

1.9177 0.6871 0.4567 0.3704 0.3285 0.3047 0.2758 0.2627 0.2553 0.2503 0.2468 0.2440 0.2399 0.2297 0.2132

2.1559 0.7927 0.5224 0.4188 0.3678 0.3388 0.3038 0.2885 0.2801 0.2748 0.2711 0.2683 0.2644 0.2551 0.2397

2.2857 0.8531 0.5619 0.4483 0.3918 0.3594 0.3202 0.3031 0.2939 0.2882 0.2843 0.2814 0.2775 0.2687 0.2541

Table 8 The torsional rigidity factor k of orthotropic and inhomogeneous square section   Gax ¼ 5G0x . k

a2 ¼ 0:1

a2 ¼ 0:3

a2 ¼ 0:5

a2 ¼ 0:7

a2 ¼ 3

a2 ¼ 5

a2 ¼ 7

0 2 4 6 8 10 15 20 25 30 35 40 50 100 1000

0.1335 0.0546 0.0466 0.0437 0.0418 0.0404 0.0378 0.0360 0.0348 0.0338 0.0330 0.0324 0.0314 0.0292 0.0267

0.3285 0.1356 0.1124 0.1037 0.0986 0.0951 0.0893 0.0855 0.0828 0.0808 0.0792 0.0779 0.0759 0.0712 0.0657

0.4705 0.1956 0.1587 0.1448 0.1371 0.1320 0.1241 0.1193 0.1159 0.1133 0.1113 0.1097 0.1072 0.1013 0.0942

0.5795 0.2422 0.1935 0.1750 0.1651 0.1588 0.1495 0.1440 0.1402 0.1374 0.1352 0.1334 0.1307 0.1241 0.1160

1.0654 0.4602 0.3457 0.3020 0.2803 0.2677 0.2518 0.2441 0.2395 0.2364 0.2340 0.2322 0.2294 0.2226 0.2132

1.1977 0.5239 0.3895 0.3375 0.3117 0.2967 0.2782 0.2697 0.2648 0.2616 0.2593 0.2575 0.2549 0.2486 0.2396

1.2698 0.5598 0.4149 0.3582 0.3298 0.3134 0.2930 0.2838 0.2786 0.2753 0.2729 0.2712 0.2686 0.2627 0.2540

Table 9 The torsional rigidity factor k of orthotropic and inhomogeneous square section   Gax ¼ 7G0x . k

a2 ¼ 0:1

a2 ¼ 0:3

a2 ¼ 0:5

a2 ¼ 0:7

a2 ¼ 3

a2 ¼ 5

a2 ¼ 7

0 2 4 6 8 10 15 20 25 30 35 40 50 100 1000

0.1868 0.0651 0.0528 0.0489 0.0466 0.0449 0.0417 0.0396 0.0379 0.0366 0.0356 0.0348 0.0335 0.0304 0.0267

0.4599 0.1632 0.1279 0.1157 0.1090 0.1045 0.0973 0.0926 0.0892 0.0866 0.0845 0.0828 0.0802 0.0737 0.0658

0.6588 0.2369 0.1811 0.1613 0.1508 0.1441 0.1340 0.1280 0.1237 0.1205 0.1179 0.1158 0.1126 0.1045 0.0942

0.8113 0.2949 0.2214 0.1948 0.1811 0.1726 0.1604 0.1535 0.1487 0.1452 0.1423 0.1401 0.1365 0.1276 0.1160

1.4916 0.5748 0.4025 0.3374 0.3055 0.2873 0.2647 0.2543 0.2482 0.2440 0.2410 0.2387 0.2351 0.2264 0.2132

1.6768 0.6591 0.4569 0.3791 0.3407 0.3187 0.2918 0.2799 0.2732 0.2688 0.2658 0.2635 0.2601 0.2521 0.2397

1.7778 0.7071 0.4892 0.4041 0.3616 0.3372 0.3073 0.2941 0.2869 0.2823 0.2791 0.2768 0.2735 0.2659 0.2541

For each n, Eqs. (40)–(43) can be solved simultaneously to get Ain and Bin . As long as the function wi is obtained, the warping function /i ¼ wi  xy can be determined. The shear stresses can also be obtained from Eq. (25). The torsional rigidity of the section is

     @/ @/ þ x x  Gxz  y y dxdy Gyz @y @x R  Z Z  m X @/ @/  Gxi a2i x2 þ y2 þ a2i x i  y i dxdy @y @x Ri i¼1  Z Z  m X @w @w ¼ Gxi 2y2 þ a2i x i  y i dxdy @y @x R i i¼1

D¼

Z Z 

ð44Þ

Using Eq. (39), we have

D

m X i¼1

( " #) 2 1 X n1 4 8b 3 xi 2 2 Gxi Dxi b þ ð1Þ 2ai F in  2 2 fin ðxÞjxi1 3 n p n¼1;3;5;... ð45Þ

where Dxi ¼ xi  xi1 denotes the width of the ith sub-rectangle and

F in ¼

Ain ðai nn Þ2

fðai nn xi  1Þ exp½ai nn ðxi  aÞ  ðai nn xi1  1Þ

 exp½ai nn ðxi1  aÞg 

Bin ðai nn Þ2

½ðai nn xi þ 1Þ expðai nn xi Þ

 ðai nn xi1 þ 1Þ expðai nn xi1 Þ

ð46Þ

If the torque applied on the ends of the bar is denoted by M, the twist angle of unit length is

h¼

M D

ð47Þ

5. Numerical example Firstly the torsional rigidity of a homogeneous and orthotropic rectangular section is investigated for different size ratio s ¼ a=ð2bÞ and shear modulus ratio a, and the results are tabulated in Table 1. The number of the item of the infinite series is dependent on the relative round error, which is set to be 106 in all examples. We find the convergent number increases with the decrease of the size ratio s in our calculation. In Table 1, the row corresponding to a2 ¼ 1 is for isotropic materials and the results in parenthesis are obtained by Timoshenko and Goodier [2]. To demonstrate the convergence of the layerwise model, the torsional rigidity of an orthotropic and inhomogeneous rectangular section is computed using two different methods. The resulting tor3 sional rigidity factor k ¼ D=ð2ba G0x Þ versus m is plotted in Fig. 3, in which the solid line is obtained from the layerwise model and the dashed line from the exact solution. It is shown that the solid line approaches to the dashed line with the increase of m. Some numerical experiments show that the layer number m = 100 can give the results accurate enough (relative round error less than 106 Þ and therefore the layer number m is taken 100 in this example. Secondly, we consider a square composite cross-section (as shown in Fig. 4) consisting of two symmetrically placed rectangular isotropic materials in contact, which has been studied by Sapountzakis [20]. Table 2 lists the torsional constant It ¼ D=ða4 G2 Þ and max ¼ the non-dimensional maximum resultant shear stress s smax =ðaG2 hÞ, in which smax is the maximum of the resultant shear stress, namely,

1456

X. Rongqiao et al. / Composite Structures 92 (2010) 1449–1457

 of an orthotropic square section ða ¼ 2Þ. Fig. 9. The distribution of the shear stress factor s

smax ¼ maxðsÞ; s ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2xz þ s2yz

ð48Þ

It is shown that the numerical results agree well with those obtained by Sapountzakis [20] and Muskhelishvilli [6]. Fig. 5 shows the distribution of the dimensionless warping  M ¼ / =a2 for G1 =G2 ¼ 2, in which / denotes the new function / M M warping function with respect to the shear center of the considered composite cross-section and can be transformed from the original warping function /, namely [20]

/M ¼ / þ xM y

ð49Þ

where xM is the x-coordinate of the shear center. The distribution of the dimensionless resultant shear stress s=ðaG2 hÞ for G1 =G2 ¼ 2 is also plotted in Fig. 6. It is shown that the maximum of the resultant shear stress occurred at the midpoint of the left boundary of the composite cross section and the shear stress vanishes at the four corners of the section. Thirdly, the material is considered to be isotropic and inhomogeneous, and the distribution of shear modulus obeys the following classical power law,

Gx ðxÞ ¼ G0x ½1  ðx=aÞk  þ Gax ðx=aÞk

ð50Þ

in which G0x and Gax are the shear moduli at the left ðx ¼ 0Þ and right ðx ¼ aÞ sides, respectively. When the gradient factor k is zero or infinite, the material becomes homogeneous, with the shear modulus

Gax or G0x . Apparently, G0x ¼ Gax also leads to a homogeneous material for any k. In this case only the layerwise model ðm ¼ 100Þ can be used to investigate the torsional rigidity and shear stress. Tables 3–6 list the torsional rigidity factor k for Gax ¼ 3G0x ; Gax ¼ 5G0x ; Gax ¼ 7G0x and Gax ¼ 9G0x for different gradient factors k and size ratios s. The torsional rigidity factor k decreases with the increase of the parameters k and s. In Table 3, the values of k at the second row ðk ¼ 0Þ are three times large as those at the last row ðk ¼ 1000Þ since these two rows correspond to homogeneous materials with shear modulus Gax and G0x , respectively. The similar results are found in Tables 4–6 except that the multiple becomes 5, 7 and 9. Fig. 7 shows the distribution of the non-dimensional shear stress factor s, which is defined as

s ¼ s=ðaG0x hÞ

ð51Þ Gax

3G0x .

for k ¼ 0; 1; 2; 3 and ¼ Fig. 7a shows that the maximum shear stress occurs at the midpoint of the four sides of the square section for homogeneous materials ðk ¼ 0Þ. The zero stress occurs at the four corners and the center of the square section. With the increase of gradient factor k, the location of the maximum shear stress in the y direction is retained at the midpoint of the vertical side since the material properties is homogeneous in the y direction. However, the location of the maximum shear stress in the x direction shifts to the side with a larger magnitude of the shear modulus and so does the location of the zero stress inside the section. Fig. 8 shows the distribution of the total shear stress factor of a rectangular section of

X. Rongqiao et al. / Composite Structures 92 (2010) 1449–1457

1457

 of an orthotropic rectangular section ða ¼ 2; s ¼ 0:5Þ. Fig. 10. The distribution of the shear stress factor s

size ratio s ¼ 0:5 and a similar distribution of the shear stress can be obtained. Finally, the material of the section is considered to be orthotropic and the shear modulus varies in the same law as described in Eq. (50). Tables 7–10 list the torsional rigidity factor k of the square section for Gax ¼ 3G0x ; Gax ¼ 5G0x ; Gax ¼ 7G0x and Gax ¼ 9G0x . It is shown the factor k also decreases with the increase of the parameter k, as indicated for isotropic materials in Tables 3–6. Meanwhile, the factor k increases with the parameter a. In Tables 7–10, we also find that the values of k at the second row ðk ¼ 0Þ are 3, 5, 7 and 9 times large as those at the last row ðk ¼ 1000Þ. Figs. 9 and 10 give the dis for a ¼ 2. Since the material is tribution of the shear stress factor s orthotropic the shear stresses in the x and y directions are not equal even for a square section of homogeneous material (Fig. 9a) and the maximum shear stress in the y direction is greater than the one in the x direction as Gyz > Gxz ða > 1Þ. For the section of inhomogeneous material, the location of the maximum shear stress in the x direction offsets to the side with a larger shear modulus as well as the location of zero shear stress inside the section. It is similar to those for an isotropic section. 6. Conclusions The Saint-Venant torsion of a straight bar with a rectangular section is described in this paper. The material can be orthotropic and inhomogeneous or functionally graded in one direction. The exact solution is obtained for a special case where the variation of the shear modulus obeys the exponential law. The approximate solution is also presented for an arbitrary variation in one direction of the shear modulus using the layerwise model, in which the original section is represented by a composite section with multiple subrectangular regions. The warping function, shear stress and torsional rigidity are obtained in an infinite series through the variable separation method. In this solution, the hyperbolic functions are not used purposely to avoid the numerical difficulties. Finally through the numerical examples, it is clearly shown that this proposed method is very effective for the Saint-Venant torsion of rectangular sections with orthotropic and inhomogeneous or functionally graded material. Acknowledgement This work is supported by the National Natural Science Foundation of China (Nos. 10872180, 10725210).

References [1] Timoshenko SP. History of strength of materials. New York: McGraw-Hill; 1953. [2] Timoshenko SP, Goodier JN. Theory of elasticity. 3rd ed. New York: McGrawHill; 1970. [3] Courant R. Variational methods for the solution of problems of equilibrium and vibration. Bull Am Math Soc 1943;49(1):1–23. [4] Timoshenko SP. Strength of materials. 3rd ed. Berkshire (England): Van Nostrand; 1956. [5] Quinlan PM. The torsion of an irregular polygon. Proc Royal Soc Lond Ser A Math Phys Sci 1964;282(1389):208–27. [6] Muskhelishvilli NI. Some basic problems of the mathematical theory of elasticity. Groningen (Holland): P. Noordhoff; 1953. [7] Booker JR, Kitipornchai S. Torsion of multilayered rectangular section. J Eng Mech Div ASCE 1971;97(EM5):1451–68. [8] Kuo YM, Conway HD. The torsion of composite tubes and cylinders. Int J Solids Struct 1973;9(12):1553–65. [9] Kuo YM, Conway HD. Torsion of cylinders with multiple reinforcement. J Eng Mech Div ASCE 1974;100(EM2):221–34. [10] Kuo YM, Conway HD. Torsion of composite rhombus cylinder. J Appl Mech 1974;31(2):302–3. [11] Kuo YM, Conway HD. Torsion of reinforced square cylinder. J Eng Mech Div ASCE 1980;106(EM6):1341–7. [12] Packham BA, Shail R. St. Venant torsion of composite cylinders. J Elasticity 1978;8(4):393–407. [13] Ripton R. Optimal fiber configurations for maximum torsional rigidity. Arch Ration Mech Anal 1998;144(1):79–106. [14] Chen T, Benveniste Y, Chuang PC. Exact solutions in torsion of composite bars: thickly coated neutral inhomogeneities and composite cylinder assemblages. Proc Math Phys Eng Sci 2002;458(2023):1719–59. [15] Ely JF, Zienkiewicz OC. Torsion of compound bars – a relaxation solution. Int J Mech Sci 1960;1(4):356–65. [16] Herrmann LR. Elastic torsional analysis of irregular shapes. J Eng Mech Div ASCE 1965(EM6):11–9. [17] Jaswon MA, Ponter AR. An integral equation solution of the torsion problem. Proc Roy Soc Lond Ser A Math Phys Sci 1963;273(1353):237–46. [18] Kasikadelis JT, Sapountzakis EJ. Torsion of composite bars by boundary element method. J Eng Mech 1986;111(9):1197–210. [19] Sapountzakis EJ. Solution of non-uniform torsion of bars by an integral equation method. Comput Struct 2000;77(6):659–67. [20] Sapountzakis EJ. Nonuniform torsion of multi-material composite bars by the boundary element method. Comput Struct 2001;79(32):2805–16. [21] Koizumi M. The concept of FGM. Ceram Trans Function Grad Mater 1993;34(1):3–10. [22] Plunkett R. Torsion of inhomogeneous elastic prismatic bars. J Eng Ind 1965;87:391–2. [23] Rooney FJ, Ferrari M. Torsion and flexure of inhomogeneous elements. Compos Eng 1995;5(7):901–11. [24] Rooney FJ, Ferrari M. On the St. Venant problems for inhomogeneous circular bars. J Appl Mech 1999;66(2):32–44. [25] Horgan CO, Chan AM. Torsion of functionally graded isotropic linearly elastic bars. J Elasticity 1999;52(2):181–99. [26] Tarn JG, Chang HH. Torsion of cylindrically orthotropic elastic circular bars with radial inhomogeneity: some exact solutions and end effects. Int J Solids Struct 2008;45(1):303–19.