An improved hybrid-stress element approach to torsion of shafts

Computers and Structures 71 (1999) 535±563

An improved hybrid-stress element approach to torsion of shafts Q.Z. Xiao a, 1, B.L. Karihaloo b,*, Z.R. Li a, F.W. Williams b a

Department of Modern Mechanics, University of Science and Technology of China, Hefei, 230026, People's Republic of China b Cardi€ School of Engineering, Cardi€ University, Cardi€, CF2 3TB, UK Received 30 April 1998; accepted 12 November 1998

Abstract A four-noded torsional element with only four stress parameters is developed based on the Hellinger±Reissner principle by optimizing element trial stresses through the introduction of incompatible displacements. Two boundary hybrid-stress elements with six stress parameters are also developed to meet the traction-free condition on one or two adjacent sides of a four-noded element, respectively. Numerical results show that these elements give accurate results for the warping displacement, the angle of twist per unit length, as well as the shear stresses for shafts with a polygonal cross-section. Comparison with alternative displacement and stress elements shows that the present elements give results of the same accuracy with a coarser mesh. The elements are shown to be insensitive to mesh distortion and to exhibit rapid convergence with h-re®nement. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: St Venant's torsion of shaft; Optimization of hybrid-stress element; Hybrid-stress torsional element; Boundary hybridstress element

1. Introduction The solution of the St Venant torsion problem of shafts may be regarded as a settled question in theory. There are indeed solutions available for shafts of various cross sections [1,2]. Unfortunately, it is dicult to obtain analytical solutions for shafts with complicated cross sections, for example, those with a slot, cutout or ®llet at a re-entrant corner. In order to solve St Venant's problem for such shafts, a lot of simpli®cations have to be made, thus rendering the solutions of dubious practical value. There is thus a need for robust

* Corresponding author. Tel.: +44-1222-874-934; fax: +441222-874-597. E-mail address: karihaloob@cardi€.ac.uk (B.L. Karihaloo) 1 Currently Academic Visitor at Cardi€ School of Engineering, Cardi€ University, Cardi€ CF2 3TB, UK.

and accurate numerical methods for the solution of St Venant's torsion of shafts in order ®rst to examine the validity of already published results for various complicated cross sections and then to solve more complicated problems of practical interest. Courant [3] was the ®rst to use a triangular element based on the complementary energy principle to calculate the St Venant torsion stress function. Since then, the ®nite element methodology has advanced greatly [4±7]. It has become the most widely employed numerical technique for the solution of problems in solid mechanics. Studies on St Venant's torsion of shafts have been carried out by many researchers. Herrmann [8], Kawai and Yoshimura [9], and Kruhula and Lauterbach [10] used the displacement formulation in terms of the warping function; Synge and Cahill [11], Zienkiewicz and Cheung [12], Feldmann [13], Moan [14], Golley and Cashman [15], Heise [16], and Bathe

0045-7949/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 8 ) 0 0 2 4 5 - 4

536

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

[6] used the stress approach based on the Prandtl stress function; Yamada et al. [17,18], and Noor and Andersen [19] used the hybrid stress and the mixed approaches, respectively. Desai discusses these ®nite element approaches in his book [20]. Based on the Prandtl stress function, Jaswon and Ponter [21], and Mukherjee and Morjaria [22] developed the boundary element approach, and Zielinski and Zienkiewicz [23] the Tre€tz method. Also based on the Prandtl stress function, some researchers have recently applied the di€erential quadrature method to the torsion of bars [24] and shafts [25±27]. The boundary element, the Tre€tz and the di€erential quadrature methods mentioned above are all stress-based approaches. The authors using displacement or stress-based approaches usually solve the problem for a ®xed applied angle of twist. In practical cases, the angle of twist under the prescribed twisting moment is unknown, so that the relation between the twisting moment and the warping/ stress function is an additional equation to be solved simultaneously by the displacement/stress approach. This requires complicated integrations over the crosssection. In addition, the traction-free conditions on the boundary cannot be correctly enforced by the displacement approach, and the warping function or the warping displacement cannot be easily determined by the stress approach. It is also dicult to obtain accurate stresses by di€erentiating the warping or the stress function. The unknowns of most interest in the torsion of shafts are the torsional rigidity, the angle of twist per unit length, the warping displacement, and the stresses in the cross-section. For determining these unknowns simultaneously and with high accuracy, hybrid ®nite elements are necessary. Based on the classical formulation of the hybrid stress element [28] using the modi®ed complementary functional, Yamada et al. [17,18] introduced a triangular hybrid stress element in which the warping function, the angle of twist per unit length and the stress function are chosen as the independent unknowns. This overcomes most of the diculties mentioned above, except that the stresses are obtained by di€erentiating the stress function. It has been found that a more convenient and e€ective way to formulate the hybrid stress elements is an optimization approach [7,29] based on the Hellinger±Reissner functional [30]. Wu and Pian [7] have developed a systematic and practical approach to constructing and optimizing hybrid stress elements. In theory, by using the hybrid stress model based on the Hellinger± Reissner principle, all the above unknowns in the torsion of shafts including the stresses can be simultaneously obtained. In addition, the traction-free conditions on the boundary can be easily enforced, so that accurate distribution of stresses in the cross-section can be expected.

In this paper, we will further develop the hybrid stress approach to the solution of torsion of shafts. This development will take two forms. First, by optimizing the hybrid stress element through the introduction of incompatible displacement ®elds, a new fournoded hybrid-stress element for the torsion of shafts will be introduced. Secondly, two special boundary hybrid-stress elements will be introduced to meet exactly the traction-free boundary constraints. The accuracy of these new elements will be judged by solving problems for which the analytical solutions are known. The elements will then be used to solve shafts of complicated cross sections for which no exact solutions are possible. Details of mesh distortion sensitivity, convergence rate and comparison with alternative ®nite element methods are given in the Appendix.

2. Hellinger±Reissner principle and hybrid-stress approach Consider a uniform shaft of arbitrary cross-section twisted by couples applied at the ends. Without loss of generality, the origin of co-ordinates is taken at the left end cross-section, with the x- and y-axes as the principal axes of inertia, and the z-axis along the axis of the shaft and pointing to its other end. According to St Venant's theory of torsion [1,2,31], the displacement components ui, i = 1, 2, 3, are u1 ˆ ÿyzy u2 ˆ yzx u3 ˆ w ˆ w…x,y†

…1†

where y represents the angle of twist per unit length (clockwise about the z-axis). With the assumed displacements, Eq. (1), the nonzero components of strains are 3 2 @   ÿy 7 6 @x gxz 7 6 gˆ ˆ Du u, Du ˆ 6 7, gyz @ 5 4 x …2† @y   w uˆ y Corresponding to Eq. (2), the non-zero components of stresses are     t Cs 0 t ˆ xz ˆ Cg, C ˆ …3† tyz 0 Cs where Cs is the modulus of rigidity.

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

In view of the non-zero components of stresses in Eq. (3), the equations of equilibrium reduce to @ txz @ tyz ‡ ˆ0 @x @y

…4†

Consider now a cross-section A of the shaft, with boundary @A and a unit outward normal n~(`, m, n ). On the surface of the shaft, which is free from external forces and has normals perpendicular to the z-axis, the boundary conditions reduce to txz ` ‡ tyz m ˆ 0

…5†

The applied torque at the ends of the twisted shaft is … Mt ˆ …tyz x ÿ txz y†dA …6† A

The boundary-value problem, Eqs. (2)±(6), is equivalent to the stationary condition of the Hellinger± Reissner functional [31]  …  1 T ÿ1 T PHR …t,u† ˆ ÿ t C t ‡ t Du u dA ÿ yMt 2 A dPHR …t,u† ˆ 0

…7†

Based on the Hellinger±Reissner principle, Eq. (7), we will ®rst consider a hybrid-stress formulation in which the warping displacement w is related to nodal values qw via the shape functions Nw     qw Nw 0 , qˆ u ˆ Nq, N ˆ …8† 0 1 y

537

the sti€ness matrix of element ``e'' K e ˆ G T H ÿ1 G

…12†

and the discretized equations of equilibrium   X 0 K eq ˆ Mt e

…13a†

Rewriting the overall sti€ness matrix of the system into blocks   X Kww Kwy Kˆ Ke ˆ Kyw Kyy e and denoting the vector of total nodal warping displacements as Qw, Eq. (13a) can be rewritten as      Kww Kwy Qw 0 …13b† ˆ Kyw Kyy y Mt The solution of the system of Eq. (13b) gives the discretized warping displacements and angles of twist per unit length. Qw ˆ ÿK ÿ1 ww Kwy y,

…Kyy ÿ Kyw K ÿ1 ww Kwy †y ˆ Mt

where (Kyy ÿ KywK ÿ1wwKwy ) is the torsional rigidity. The distribution of warping displacements can then be obtained from Eq. (8), and the stress parameters and the distribution of stresses from Eqs. (11) and (9). Note however that the warping displacement of at least one node has been prescribed to suppress rigid body motion.

and the stress is related to stress parameters b via stress interpolation function j t ˆ jb

…9†

Note that as y is constant in any cross-section, it is unnecessary to interpolate it. Substituting Eqs. (8) and (9) into the energy function in Eq. (7), and denoting the characteristic matrices of element ``e'' by … … Hˆ jT C ÿ1 jdA, G ˆ jT Du NdA Ae

w ˆ wq ‡ wl ˆ Nw qw ‡ Nl l

…14†

condition …11†

The optimization condition (OPC) on stress presented in [7,29] reduces to the following for the torsion

where A e represents the area of element ``e'', we have   X 1 PHR …b,q† ˆ ÿ bT Hb ‡ bT Gq ÿ yMt …10† 2 e

b ˆ H ÿ1 Gq

Let us append an incompatible term to the warping displacement in Eq. (8) as follows:

in which the compatible displacement wq is related to the nodal values qw via the shape functions Nw, while the incompatible term wl is related to the element inner parameters l via the shape functions Nl. The assumed stress is divided into lower and higher order parts:   b t ˆ jb ˆ ‰jI jII Š I …15† bII

Ae

Making use of the stationary dPHR(b,q ) = 0, Eq. (10) gives

3. A four-noded hybrid stress element for torsion of shafts

538

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

of shafts: … @ Ae

where

wTl ‰` m Štds ˆ 0

…16a†

where @A e is the boundary of element ``e''. Physically, Eq. (16a) means that along the element boundary stresses do not work on the incompatible part of the displacement. Making use of Eqs. (14) and (15), the OPC on stress corresponding to Eq. (16a) is … @ Ae

N Tl ‰`

 b jII Š I ds ˆ 0 bII

Ni ˆ

jI ˆ

m Š‰jI

…16b†

N2

N3



…17†

where (xi, Zi ) represent the isoparametric co-ordinates of point i with the global co-ordinates (xi, yi ) = 1, 2, 3, 4. The additional shape functions for the incompatible displacement in Eq. (14) are  Nl ˆ 1 ÿ x2

1 ÿ Z2



2

b2 6 ÿb 3 6 bII ˆ 6 4 b2 a1

Fig. 1. A four-noded element.

0 , x

jII ˆ



x 0 , 0 Z



8 9 b > > > = < 1> .. , bI ˆ . > > > ; : > b4

a2 b3 a2 ÿ a1

0 b1 a1

3 a3 b3 7 7 7bI 5 0

…20†

Eliminating bII from Eqs. (19) using (20), we have the following optimized trial stresses b2 61 ÿ b x 3 6 t ˆ 6 4 b2 Z a1

a2 x b3 a2 1ÿ Z a1

Z b1 Z a1

38 9 a3 b > > x > = < 1> b3 7 7 b2 7 b > 5> > ; : 3> x b4

…21a†

or, equivalently, 

a1 b3 ÿ a1 b2 x a1 a2 x a1 Z a3 x a1 b3 ÿ a2 b3 Z b1 Z b3 x b2 b3 Z 8 9 b > > > < 1> = b2 ˆ j b b > > > : 3> ; b4

t ˆ

…18†

…19†





Substituting Eqs. (18) and (19) into Eq. (16b), gives

The primary assumed stresses take the form 8 9 b > >  > < 1> = 1 0 Z 0 x 0 .. tˆ . > 0 1 0 x 0 Z > > : > ; b6   b ˆ ‰jI jII Š I bII

b5 b6

2

N4 Š,

1 …1 ‡ xxi †…1 ‡ ZZi † 4

1 0 Z 0 1 0

bII ˆ



We will use the OPC, Eq. (16b), for constructing a four-noded hybrid-stress element for the torsion of shafts. Referring to the four-noded element shown in Fig. 1, the shape functions for compatible warping displacement in Eqs. (8) or (14) are Nw ˆ ‰N1





…21b†

The coecients ai and bi (i = 1, 2, 3) are dependent on the element nodal co-ordinates as follows: 2

a1 4 a2 a3

2 2 3 x 3 b1 ÿ1 1 1 ÿ1 6 1 1 x2 b2 5 ˆ 4 1 ÿ1 1 ÿ1 56 4 x3 4 ÿ1 ÿ1 1 1 b3 x4

3 y1 y2 7 7 y3 5 y4

In the element formulation discussed in Section 2, we now use the assumed displacement (8), in which Nw is de®ned by Eq. (17), and the optimized trial stresses Eq. (21b). The resulting element which has four stress parameters is designated H4. 2  2 Gauss quadrature is employed for the element formulation. In each element, the number of independent nodal displacements is nd = dim(q ) = 5, the degree of rigid body motion = 1, while the number of independent stress parameters is nb = dim(b ) = 4. It follows therefore that nb = nd ÿ 1, so that we have ensured the best parameter matching for hybrid-stress elements [7,29].

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

539

4. Hybrid-stress elements for satisfying traction-free boundary constraints We now consider two special hybrid-stress elements for satisfying the traction-free conditions on the boundary. One of the elements ensures a traction-free side, while the second has two adjacent traction-free sides. We will discuss each of them in turn in the following. 4.1. A four-noded element with a traction-free side The primary trial stresses that satisfy the equilibrium equation have the form of a second order complete polynomial 2 6 tˆ4

1 0

x

0

0 1 ÿy x

y

x2

0

0 ÿ2xy x 2

xy 1 ÿ y2 2

38 9 b > > y2 > = < 1> 7 . 5 .. > 0 > > ; : > b9

…22†

Let the equation of the traction-free side be `x ‡ my ˆ p

…23†

where (`, m ) are direction cosines of the unit outward normal to the side, and p is a constant. The tractionfree condition to be satis®ed is Eq. (5). We now proceed as follows: write x (or y ) as a function of y (or x ) in Eq. (23), and substitute this function into Eq. (22), and then the latter into Eq. (5). In order to meet the constraint (5) at every point whose co-ordinates satisfy (23), we will have three equations in b. Thus, three b parameters will be linearly dependent on the remaining ones. Eliminating these three b parameters from Eq. (22), we have 1. For `$0, uyx = (m/`) and px = ( p/`), 2

ÿuyx ÿpx ‡ x ‡ 2uyx y ÿuyx px ‡ u2yx y

6 6 6 6 6  t ˆ6 6 ÿp2x ‡ 4uyx px y ‡ x 2 ÿ 3u2yx y2 6 6 6 ÿuyx p2x ‡ 2u2yx px y ÿ u3yx y2 4 ÿpx y ‡ xy ‡ 1:5uyx y2

1 ÿy x

3T

7 7 7 7 7 ÿ2xy 7 7 7 7 2 x 7 5 2 ÿ0:5y

8 9 b > > > < 1> = .. . > > > > : ; b6

…24a†

2. For m$0, uxy = (`/m ) and py = ( p/m ), 2

1 x y

6 6 6 6 6  t ˆ6 2 6x 6 xy 6 4 y2

ÿuxy py ÿ 2uxy x ÿ y ÿuxy py ‡ u2xy x

3T

7 7 7 7 7 7 2 2py x ÿ 2xy ÿ 3uxy x 7 2 2 2 27 0:5py ÿ 2uxy py x ‡ 1:5uxy x ÿ 0:5y 7 5 ÿuxy p2y ‡ 2u2xy py x ÿ u3xy x 2

8 9 b > > > < 1> = .. . > > > : > ; b6

…24b†

Note that as truncation errors during computations are inevitable, care must be taken to avoid division by a very small number. Thus, it is advisable to use Eq. (24a), when v`vrvmv, and Eq. (24b) otherwise.

540

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

4.2. A four-noded element with two adjacent traction-free sides The primary trial stresses are now chosen as  tˆ

1 0 x 0 1 0

0 y 0 x 0 y

2

x 0

0 x2

xy 0 0 xy

2

y 0

9 8 b > > > = < 1 > 0 .. 2 . > y > > > ; : b12

…25†

Let the equations of the two traction-free sides have the same form as Eq. (23), with (`, m ) and (`, m ) the direction cosines of the unit outward normals to these two sides, respectively, and p and p the respective constants. The traction-free constraint for each side is again Eq. (5). Following the procedure of Section 4.1, we get three additional equations in b. In order to meet the constraint (5) at every point on each of the two sides, we will thus have six equations relating the 12 bs. Thus, six b parameters will be linearly dependent on the remaining six. Eliminating these six b parameters from Eq. (25), we get the required trial stress. In view of the fact that the traction-free sides are adjacent, ` and ` cannot vanish simultaneously. Therefore, without loss of generality, we assume `$0, and denote uyx = (m/`), px = ( p/`). Now we consider all remaining practical situations. 1. For `$0, uyx = (m/`) and px = (p/`), 2

c7 u yx ÿ p x ‡ x

ÿc7 ‡ y

3T

7 6 c3 7 6 c3 c8 ÿ u yx p x ÿ u yx uyx y ÿ ‡ x ‡ c1 y 7 6 Duyx 7 6 7 6 7 6 c 3 2 7 6 c2 c7 u yx ÿ p 2x ‡ 2…u yx p x ÿ c3 c8 †y ‡ x 2 ‡ u yx uyx y2 ÿc c ‡ 2 y ÿ c y 2 7 1 7  6 Du yx 7 6  t ˆ6 7  b1 c c 4 5 7 6 2 2 2 2 2 2 ‡2 y ‡ x ÿ c6 y 7 6 c4 c8 ÿ u yx p x ‡ u yx …ÿu yx ‡ c6 †y ‡ 2…u yx p x ÿ c5 c8 †y ÿ 7 6 Duyx Duyx 7 6 7 6 …ÿp x ‡ c7 u yx †y ‡ xy ÿc7 y ‡ y2 7 6 7 6 5 4 c 3 …ÿu yx p x ‡ c3 c8 †y ÿ u yx uyx y2 ÿ y ‡ xy ‡ c1 y2 Duyx






b6

T

…26a† where Duyx ˆ uyx ÿ u yx , Dpx ˆ px ÿ p x , c1 ˆ uyx ‡ u yx , c2 ˆ px ‡ p x , c3 ˆ uyx px ÿ u yx p x , c4 ˆ uyx p2x ÿ u yx p 2x , c5 ˆ u2yx px ÿ u 2yx p x , c6 ˆ u2yx ‡ uyx u yx ‡ u 2yx , c7 ˆ

u yx Dpx , c8 ˆ Duyx Duyx

2. For m = 0, m$0, uxy = (`/m ) and py = (p/m ), 2

ÿpx ‡ x 6 0 6 6 6 ÿp2 ‡ x 2 6 x 6 t ˆ 6 6 …x ÿ px †y 6 6 6 0 6 4 0

u xy … px ÿ x† ÿp y ‡ u xy x ‡ y

3T

7 7 7 7 2 2 u xy … px ÿ x † 7 h i7 7 2 u xy p y px ÿ …p y ‡ u xy px †x ‡ u xy x 7 7 7 7 2 ÿp y x ‡ u xy x ‡ xy 7 5 2 2 2 2 ÿp y ‡ 2u xy p y x ÿ u xy x ‡ y

8 9 b > > > < 1> = .. > . > > > : ; b6

…26b†

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

541

3. For `=0, m$0 and py = (p/m ), 2

ÿpx ‡ x ‡ uyx y uyx …p y ÿ y†

6 6 6 6 ÿp2x ‡ 2uyx px y ‡ x 2 ÿ u2yx y2 6 6  t ˆ6 ÿpx y ‡ xy ‡ uyx y2 6 6 h i 6 6 uyx px p y ÿ … px ‡ uyx p y †y ‡ uyx y2 6 4 uyx …p 2y ÿ y2 †

3T

0

ÿp y ‡ y 7 7 7 7 0 7 7 7 0 7 7 7 x… y ÿ p y † 7 7 5 y2 ÿ p 2y

8 9 b > > > < 1> = .. . > > > : > ; b6

…26c†

We can select v`vrv`v, and use Eq. (26a) if v`vrvmv. If v`v R vmv, then we use Eq. (26b) or Eq. (26c). Sometimes, it may be necessary to choose suitable global co-ordinates. The boundary hybrid-stress elements with one or two adjacent traction-free sides have each six stress parameters, and are therefore designated BH6. 3  3 Gauss quadrature is employed for the element formulation. Theoretical analysis and results of eigenvalue check show that both H4 and BH6 are rank-sucient. If instead of Eq. (25), one chose a third-order complete polynomial for the primary trial stress, i.e. 2 6 tˆ4

1 0

x

0

0 1 ÿy x

y

x2

0

0 ÿ2xy x

2

xy 1 ÿ y2 2

y2 0

x3

x 2y

0 2

ÿ3x y

x

3

2

ÿxy

xy2 1 ÿ y3 3

y3 0

3 7  5  b1






b14

T

an alternative boundary element meeting the traction-free constraints on two adjacent sides of a four-noded element would follow. This element also has six stress parameters, but an eigenvalue analysis shows it has a zero energy mode (ZEM), i.e. the element is rank-de®cient.

5. Numerical results and discussion In this section, we solve the torsion of shafts with an equilateral triangular or a square section, ®rst by using H4 only and then by using H4, together with BH6 (denoted H4 + BH6), and compare the approximate

solutions with known analytical solutions to judge the accuracy of the hybrid-stress elements. We will follow this by using H4 + BH6 for the solution of a shaft whose square cross-section has a slot, and of a shaft with Z-section. For the latter geometry, we will study the e€ect of introducing a ®llet in the re-entrant

Table 1 Results for torsion of an equilateral triangular shafta H4 Point 1 2 3 4 5 6 7 Torsional rigidity a

txz tyz w txz tyz w txz tyz w txz tyz

H4 + BH6

Mesh 2(b)

Mesh 2(c)

Mesh 2(b)

Mesh 2(c)

Exact

ÿ0.521/0.521 ÿ12.550 ÿ0.960 ÿ5.065 ÿ2.915 0.951 ÿ10.58/ ÿ 11.09 6.686/5.827 ÿ0.955 0.0 5.835 0.0424

ÿ0.349/0.349 ÿ13.210 ÿ0.944 ÿ3.267 ÿ1.866 0.934 ÿ11.29/ ÿ 11.69 6.926/6.283 ÿ0.937 0.0 3.749 0.0395

0.0 ÿ12.727 ÿ0.854 0.0 0.0 0.846 ÿ11.02/ ÿ 11.00 6.347/6.379 ÿ0.849 0.0 0.0 0.0406

0.0 ÿ12.876 ÿ0.865 0.0 0.0 0.856 ÿ11.15/ ÿ 11.16 6.422/6.427 ÿ0.858 0.0 0.0 0.0390

0.0 ÿ12.990 ÿ0.938 0.0 0.0 0.938 ÿ11.25 6.495 ÿ0.938 0.0 0.0 0.0385

Two values in a cell separated by a slash correspond to results in two adjacent elements.

542

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

corner. Without loss of generality, the modulus of rigidity is assumed to be Cs = 1.0, and the torque applied at each end of the shaft is Mt = 1.0. In all the ®gures to follow that illustrate the stress distribution (e.g. Fig. 3b), a line segment represents the distribution within an element. 5.1. Example 1. Torsion of an equilateral triangular shaft The exact solution for the torsion of a shaft of equilateral triangular section (length of side 2/Z3), shown in Fig. 2(a), is [1,2,32]: Torsional rigidity ˆ 0:0385, Warping displacement w ˆ y… y3 ÿ 3x 2 y†=2 Shear stresses:   txz ˆ ÿCs y… y ‡ 3xy†, tyz ˆ Cs y x ÿ 1:5…x 2 ÿ y2 Two discretized meshes will be used here. The one shown in Fig. 2(b) has 12 elements and 19 nodes, while the other, shown in Fig. 2(c), has 48 elements and 61 nodes. The latter mesh is obtained by dividing each element in the former into four sub-elements. Numerical results for the torsional rigidity, warping displacements and shear stresses at the sample points shown in Fig. 2(b), are listed in Table 1 and compared with the above exact solution. The comparison shows that the hybrid-stress element H4 gives very accurate results for torsional rigidity and warping displacement even with a very coarse mesh. The use of H4 + BH6 gives even better results for torsional rigidity and stresses, but the accuracy of the warping displacement is somewhat reduced. In the following, we will give detailed results obtained using the mesh of Fig. 2(c). Fig. 3 shows the results along the side DC: Fig. 3(a) shows the distribution of warping displacement, and Fig. 3(b) and (c) the shear stress tyz obtained by using the elements H4 and H4 + BH6, respectively. Along this side, the shear stress txz must vanish. The use of the element combination H4 + BH6 gives the exact txz, while the use of the H4 element alone produces the maximum error of ÿ3.267 at the point C. Fig. 4 shows the results along the side CB: Fig. 4(a) shows the warping displacement, Fig. 4(b) and (c) the shear stress txz, and Fig. 4(d) and (e) the shear stress tyz. The results shown in Fig. 4(b) and (d) are obtained by H4, while the results in Fig. 4(c) and (e) are obtained by H4 + BH6.

Fig. 2. Geometry and discretized meshes of equilateral triangular pro®le. (a) Equilateral triangular pro®le. (b) Mesh with 12 elements. (c) Mesh with 48 elements.

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

543

Fig. 3. Results along the side DC. (a) Warping displacement along the side DC. (b) Approximate tyz by H4 along the side DC; (± ± ±) polynomial smoothing of H4 results. (c) Approximate tyz by H4 + BH6 along the side DC; (± ± ±) polynomial smoothing of H4 + BH6 results.

544

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

Along the line of symmetry DB ( y = 0, Fig. 2(c)), the warping displacement and stress txz must vanish. In the approximate numerical solution the warping along this line nearly vanishes, whereas the maximum error in txz is less than 0.4, irrespective of whether the H4 element is used alone or in combination with BH6. The approximate tyz obtained by H4 and by H4 + BH6 are shown in Fig. 5(a) and (b), respectively. Figs. 3(a) and 4(a) showing the warping displacements on sides DC and CB con®rm that accurate results can be obtained by using H4 or H4 + BH6. In fact, the results obtained by using H4 alone are better. This is not surprising, given the fact that this element has only four stress parameters and satis®es best the parameter matching criterion for hybrid stress elements. A comparison of Figs. 3(b), 4(b), 4(d) with 3(c), 4(c), and 4(e), respectively, shows on the other hand that by using the BH6 in combination with H4 the traction-free conditions along the boundary can be better met. The approximate stress distribution obtained by H4 + BH6 and smoothed by polynomial ®tting is almost identical to the exact solution. Fig. 5(a) and (b) shows that the boundary element BH6 improves the stress distribution not only on the boundary, but also in the interior of the cross-section.

5.2. Example 2. Torsion of a square shaft Let us consider now a shaft of square cross-section (side = 2), as shown in Fig. 6(a). The analytical solution to this problem is [1,2]: Torsional rigidity ˆ

   16 192 p 1 3p Cs 1 ÿ tanh ‡ 5 tanh 3 p5 2 3 2

ˆ 2:2506 Warping displacement 2

0

p B sinh 2 y 6 32 p B w ˆ y6 p sin 2 4xy ÿ p3 @ cosh 2 13 3p 7 1 sinh 2 y 3p C ÿ xC7 sin 3p 27 2 A5 cosh 2

Shear stresses

quadquad 0

txz

p sinh y 16 B 2 cos p B ˆ ÿCs y 2 @ p p 2 cosh 2 1 3p sinh y 1 2 cos 3p x C C ÿ 3p 9 2 A cosh 2 2

0

p cosh y 6 16 B 2 sin p 6 B tyz ˆ Cs y42x ÿ 2 @ p p 2 cosh 2 13 3p y cosh 7 1 2 sin 3p x C C7 ÿ 3p 9 2 A5 cosh 2 Two meshes, one consisting of 16 elements and 25 nodes, and the other of 143 elements and 168 nodes, are shown in Fig. 6(b) and (c). The computational results for the torsional rigidity, warping displacements, and shear stresses at the sample points shown in Fig. 6(b), are listed in Table 2, and compared with the exact solution. The warping displacement and txz at points 3±7 obtained both by the use of H4 and of H4 + BH6 are identically zero, as in the exact solution, and are therefore not listed in the table. It is clear that the hybrid-stress element H4 either alone or in combination H4 + BH6 gives accurate values for the torsional rigidity, warping function and shear stresses even with a very coarse mesh. The use of the boundary element BH6 improves the accuracy of torsional rigidity and stresses. In the following, we will give detailed results obtained using the mesh of Fig. 6(c). Due to symmetry, the distributions of all ®elds on the four sides are similar, so that it is sucient to plot the computational results along one side only, say x = ÿ 1, as in Fig. 7. Fig. 7(a) shows that the approximate solutions for the warping function obtained by H4 alone or in combination H4 + BH6 are very accurate. Fig. 7(b) shows that tyz in each element obtained by H4 is constant, whereas Fig. 7(c) clearly shows that the more reasonable distribution of tyz is picked up when H4 is used in combination with BH6. The use of H4 + BH6 produces the exact solution for txz on x = ÿ 1, namely txz = 0, whereas if H4 only is used, the maximum error in txz is 0.207 at the corner. On the line of symmetry, y = 0, the warping displacement and shear stress txz must both vanish. By using H4 only, the computed warping displacement is less

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

545

Table 2 Results for torsion of a square shafta H4 Point 1 2 3 4 5 6 7 Torsional rigidity a

w txz tyz w txz tyz tyz tyz tyz tyz tyz

H4 + BH6

Mesh 6(b)

Mesh 6(c)

Mesh 6(b)

Mesh 6(c)

Exact

0.0019 0.2895 ÿ0.3037 0.0629 0.0777 ÿ0.3037/ ÿ 0.5522 ÿ0.5522 ÿ0.2196 ÿ0.0327 0.1867 0.5494 2.3747

ÿ0.0021 0.2074 ÿ0.1612 0.0600 0.0491 ÿ0.4667/ ÿ 0.5310 ÿ0.5957 ÿ0.2461 ÿ0.0352 0.2146 0.5968 2.2666

ÿ0.0026 0.0 0.0 0.0550 0.0 ÿ0.3144/ ÿ 0.4785 ÿ0.5974 ÿ0.1991/ ÿ 0.2260 ÿ0.0331 0.1936/0.1640 0.5954 2.2989

ÿ0.0039 0.0 0.0 0.0574 0.0 ÿ0.5307/ ÿ 0.4785 ÿ0.5913 ÿ0.2472 ÿ0.0356 0.2151 0.5943 2.2589

0.0068 0.0 0.0 0.0606 0.0 ÿ0.5091 ÿ0.6001 ÿ0.2487 ÿ0.0363 0.2157 0.6001 2.2506

Two values in a cell separated by a slash correspond to results in two adjacent elements.

than 0.49  10ÿ4, whereas H4 + BH6 gives the warping displacement that is less than 0.32  10ÿ4. The computed txz is less than 0.1  10ÿ5 both by using H4 and H4 + BH6. The non-vanishing component of the shear stresses tyz is shown in Fig. 8. The value shown is at the midside obtained by averaging the end values. The averaging does not introduce an error because the stress interpolation functions of H4 element are linear. The error in the BH6 is very small and can be ignored. Fig. 8 also con®rms the fact that a hybrid stress element gives the most accurate stresses in the middle of the element. The above two examples illustrated that the use of the hybrid stress element H4 produces very accurate

results for the torsional rigidity, warping function, and the stresses. The inclusion of special boundary hybridstress element BH6 further improves the computed stress distributions. Having gained con®dence in the accuracy of the developed hybrid-stress elements, let us now solve the problem of torsion of a shaft of square section with a slot and of a shaft with Z-section, for which no exact solutions are possible. 5.3. Example 3. A square cross-section with a slot The geometry and mesh of a shaft of square crosssection (side = 2) with a square slot (side = 0.4) are shown in Fig. 9. The mesh consists of 134 elements

546

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

Fig. 4. Results along the side CB. (a) Warping displacement along the side CB. (b) Approximate txz by H4 along the side CB; (± ± ±) polynomial smoothing of H4 results. (c) Approximate txz by H4 + BH6 along the side CB; (± ± ±) polynomial smoothing of H4 + BH6 results. (d) Approximate tyz by H4 along the side CB; (± ± ±) polynomial smoothing of H4 results. (e) Approximate tyz by H4 + BH6 along the side CB; (± ± ±) polynomial smoothing of H4 + BH6 results.

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

Fig. 5. Results along the line DB ( y = 0). (a) Approximate tyz by H4 along the side DB; (





) polynomial smoothing of H4 results. (b) Approximate tyz by H4 + BH6 along the side DB; (





) polynomial smoothing of H4 + BH6 results.

547

Fig. 6. Geometry and discretized meshes for a shaft of square pro®le. (a) Square pro®le. (b) Mesh with 16 elements. (c) Mesh with 143 elements.

548

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

Fig. 7. Numerical results on the side x = ÿ 1. (a) Warping on the side x = ÿ 1. (b) Approximate tyz by H4 on the side x = ÿ 1. (c) Approximate tyz by H4 + BH6 on the side x = ÿ 1.

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

549

Fig. 8. Approximate tyz by H4 and H4 + BH6 on the side y = 0.

Fig. 10. Numerical results on the line x = ÿ 1. (a) Approximate warping function along x = ÿ 1 by H4 + BH6 with or without a slot. (b) Approximate tyz along x = ÿ 1 by H4 + BH6 with or without a slot; (± ± ±) polynomial smoothing.

Fig. 9. Geometry and discretized mesh of a square section with a slot. (a) A square pro®le with a slot. (b) Discretized mesh.

and 162 nodes. The computed torsional rigidity is 1.7223 which when compared with a section without the slot (example 2 above) shows that the torsional rigidity has decreased by 23.47%. The warping displacement and the shear stress tyz on the line x = ÿ 1 are shown in Fig. 10(a) and (b), respectively. For comparison, the results for a square section without a slot are also shown. The warping function is almost una€ected by the slot, whereas the maximum value of tyz shows an increase of about 24.15% when the slot is present. The computed txz equals the exact value (txz = 0) both with and without a slot.

550

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

0.59 at D. When a slot is present, it of course vanishes on CD. 5.4. Example 4. Torsion of shaft of Z-section

Fig. 11. Approximate tyz along y = 0 by H4 + BH6 with or without a slot.

On the line y = 0, the warping function for a square pro®le without a slot vanishes, as in example 2 above. For a square pro®le with a slot, numerical results show that the warping displacement also vanishes for x < 0. However, for x>0, it reaches an almost constant value 0.7  10ÿ2. Numerical results also show that the absolute value of txz is less than 0.1  10ÿ5. The computed tyz is shown in Fig. 11. It is evident that there is a stress concentration at the edge of the slot; the stress concentration factor exceeds 3. It is known that the presence of a small circular slot in a shaft of circular cross-section doubles the maximum tangential stress. The results of Fig. 11 show that, as was to be expected, a rectangular slot creates a more severe stress concentration. For practical purposes, we will compute results along the edges BC and CD of the slot. These are plotted in Fig. 12. Again for comparison, the corresponding results along these lines in a square pro®le without a slot are also plotted. In the latter case, the warping displacement and tyz are almost constant on BC as seen from Fig. 12(a) and (b), while txz changes linearly from 0.063 at B to ÿ0.063 at C. In the presence of a slot, txz of course vanishes on BC, but the warping displacement changes rapidly, although its average value would seem to remain constant (Fig 12. (a)). The shear stress tyz along BC increases considerably (Fig. 12(b)), causing a serious stress concentration at the corner (stress concentration factorc4.5). Computed results on the edge CD are shown in Fig. 12(c) and (d). Fig. 12(c) shows that the warping displacement increases considerably due to the slot and that there is a serious stress concentration in txz at C (Fig. 12(d)). In the absence of the slot, tyz varies almost linearly on CD from 0.28 at C to

The pro®les of a shaft of Z-section with or without ®llets at re-entrant corners and the corresponding discretized meshes are shown in Fig. 13. An approximate analytical formula based on the theory of elasticity gives the torsional rigidity of the shaft with the crosssection shown in Fig. 13(a) as 56.333. The present computational approach gives the torsional rigidity 52.711 and 55.821 for the pro®les shown in Fig. 13(a) and (c), respectively. The torsional rigidity increases when the re-entrant corners are ®lleted. Computed results for the warping function and the non-zero components of shear stress on the boundaries AB, BC and DA1 are shown in Figs. 14±16, respectively. In order to get a better understanding of the distributions of the respective ®elds in the strip, the results on the line x = ÿ 3 are given in Fig. 17. Figs. 14(a), 15(a), 16(a) and 17(a) clearly show that the warping function varies almost linearly through the thickness and along the length, and that the ®llets decrease warping. A comparison of Fig. 14(b) with Fig. 17(c) shows that tyz varies almost parabolically through the thickness and decreases rapidly away from the edge. Fig. 17(b) shows that the main component of the shear stress (txz ) varies linearly through the thickness in the regions away from the edge and the reentrant corner. However, Figs. 15(b) and 16(b) show that txz varies signi®cantly along the length for the particular ratio of the thickness to length of 1/3, chosen for this example. However, the ®llets reduce the stress concentration near the re-entrant corner, as can be seen from Fig. 16(b). In the solution of the torsion problem of simply connected thin-walled sections made up of rectangular strips by the theory of elasticity, the e€ect of ®llets at re-entrant corners is neglected and it is assumed that only tangential stresses act in the strip, parallel to its long edges, which vary linearly with its thickness but are independent of the longitudinal coordinate. The results of Figs. 15(b) and 16(b) show that these assumptions are not strictly valid, so that care should be exercised in using the approximate formula available in the literature. The computed results along the side CD are plotted in Fig. 18. For comparison, the results on the lines y = 2 and y = 0 are also given in Figs. 19 and 20. These results clearly show that, in the narrow rectangular strip in a simply connected thin-walled section (the ratio of thickness to length is 1/4 in this example), the warping function varies linearly both in the thickness (x ) and length ( y ) directions, and that tyz varies linearly through the thickness but is constant along the length, except in the regions near the corners. Thus,

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

551

Fig. 12. Warping, txz and tyz along the edges of the slot. (a) Warping along BC. (b) tyz along BC. (c) Warping along CD. (d) txz along CD.

552

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

Fig. 13. Geometry and discretized meshes of Z-section with or without ®llets. (a) Geometry of the Z-section. (b) Discretized mesh for (a) (160 elements, 205 nodes). (c) Geometry of the Z-section with ®llets. (d) Discretized mesh for (c) (152 elements, 195 nodes).

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

Fig. 14. Numerical results on the side AB. (r) Results corresponding to Fig. 13(a) and (b). (R) Results corresponding to Fig. 13(c) and (d). (a) Warping displacement. (b) Distribution of tyz; (





) polynomial smoothing of R.

553

Fig. 15. Numerical results on the side BC. (r) Results corresponding to Fig. 13(a) and (b). (R) Results corresponding to Fig. 13(c) and (d). (a) Warping displacement. (b) Distribution of txz.

554

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

Fig. 16. Numerical results on the side DA1. (r) Results corresponding to Fig. 13(a) and (b). (R) Results corresponding to Fig. 13(c) and (d). (a) Warping displacement. (b) Distribution of txz.

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

555

Fig. 17. Numerical results on the side 1±1 (x = ÿ 3). (r) Results corresponding to Fig. 13(a) and (b). (R) Results corresponding to Fig. 13(c) and (d). (a) Warping displacement. (b) Distribution of txz. (c) Distribution of tyz (102); (





) polynomial smoothing of R.

556

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

Fig. 18. Numerical results on the side CD. (r) Results corresponding to Fig. 13(a) and (b). (R) Results corresponding to Fig. 13(c) and (d). (a) Warping displacement. (b) Distribution of tyz.

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

557

Fig. 19. Numerical results on the line 2±2 ( y = 2). (r) Results corresponding to Fig. 13(a) and (b). (R) Results corresponding to Fig. 13(c) and (d). (a) Warping displacement. (b) Distribution of txz (103); (± ± ±) polynomial smoothing of R. (c) Distribution of tyz.

558

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

Fig. 20. Numerical results for distribution of tyz on the line y = 0. (r) Results corresponding to Fig. 13(a) and (b). (R) Results corresponding to Fig. 13(c) and (d).

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

for a thin strip (ratio of thickness to length R1/4), the numerical results con®rm the validity of assumptions made in the approximate solutions of the theory of elasticity. Fig. 19(b) shows that txz is very small compared with tyz, but it varies in a very complicated manner through the thickness and along the length. In the middle part of the strip, the presence of ®llets has no e€ect on the warping function and stress tyz. Near the corners, however, the ®llet a€ects both w and tyz, which vary in a complicated manner. There is high stress concentration at the corner in the absence of the ®llet. The computed results along the line D1C are given in Fig. 21. The warping function varies nearly linearly

559

in the regions remote from the outer corner, the stresses txz and tyz also vary nearly linearly, except near the outer corner, where the variation is rather complicated. The ®llet has little e€ect on the warping function and stresses. 6. Conclusions The hybrid-stress approach developed in the paper gives accurate numerical results for the torsional rigidity, warping function and shear stresses of shafts with a polygonal cross-section. Special boundary hybridstress elements that enforce traction-free constraints on

Fig. 21. Numerical results on the side D1C. (r) Results corresponding to Fig. 13(a) and (b). (R) Results corresponding to Fig. 13(c) and (d). (a) Warping displacement. (b) Distribution of txz; (





) polynomial smoothing of R. (c) Distribution of tyz; (± ± ±) polynomial smoothing of R.

560

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

the boundary further improve the accuracy of the computed stresses. In principle, the approach is also applicable to smooth cross sections, but in practice, this would require a very ®ne mesh to approximate the smooth boundary curves. As expected, the computed results show that a rectangular slot reduces signi®cantly the torsional rigidity of the section and causes serious stress concentration. The numerical solutions of the torsion problems of simply connected thin-walled sections made up of rectangular strips show that the assumptions made in the approximate analytical solutions based on the theory of elasticity are reasonable, if the ratio of strip thickness to its length does not exceed 1/4; the assumptions are no longer valid if the ratio exceeds 1/3. A comparison with alternative ®nite elements available in the literature shows that the hybrid stress element BH6 developed in the paper gives the most accurate results for the torsional rigidity, as well as the shear stresses. Moveover, the use of H4 in combination with BH6 has the advantage of giving results of the same accuracy as the alternative displacement and stress elements but with a coarser mesh. It has been shown that both H4 and BH6 are insensitive to mesh distortion, and exhibit rapid converge with h-re®nement. These elements have also been found to be very accurate for the solution of crack problems in ®nite bodies in which the traction-free conditions have to be enforced on the cracks, as well as the unloaded boundaries [33].

Acknowledgements Financial assistance from the National Natural Science Foundation of China under grant number 19772051 and the Cardi€ Advanced Chinese Engineering Center of the University of Wales Cardi€

is gratefully acknowledged. Q.Z. Xiao and Z.R. Li would also like to thank Professor C.C. Wu. The revised paper greatly bene®ted from the reviewer's comments on the original paper. Appendix A In the following, we will ®rst present some comparative studies on torsion analysis using the present method and alternative solution techniques. We will then test the sensitivity of the present method to mesh distortion. Finally, we will study the pattern of its convergence with h-re®nement. A.1. Comparison with alternative ®nite element methods We will compare the torsional rigidity of the square shaft shown in Fig. 6(a) computed by the present hybrid stress elements, H4 and BH6, with the alternative ®nite element solutions reported in [20]. The alternative methods include three constant strain or stress triangular elements, the displacement triangular element, the stress triangular element and the hybrid stress triangular element introduced by Yamada et al. [17,18], and the stress quadrilateral element. As the mesh used in the present computation is very coarse, the errors in the computed shear stresses can be expected to be large.We will therefore compare only the computed torsional rigidity. Numerical results listed in Table A1 show that H4 gives similar results to the three constant strain or stress triangular elements in this example. The numerical solutions with a ®ner mesh are closer to the exact valuesÐthe error in the torsional rigidity by the stress triangular element is reduced from 21.0% to 14.0% when the number of elements is doubled. The higher order stress quadrilateral element which satis®es the traction-free condition on the boundary gives improved results. The most

Table A1 The error in the computed torsional rigidity by di€erent ®nite element methods Method

Elements

Nodes

Error in the torsional rigidity (%)a

Displacement triangular element Stress triangular elementb Stress triangular elementc Stress quadrilateral elementd Hybrid stress triangular element H4 BH6

4 16 32 16 4 4 4

5 13 25 25 5 9 9

ÿ18.5 21.0 14.0 9.0 21.0 ÿ18.5 ÿ1.2

a

Error = 100 (exact value ÿ computed value )/exact value. A quarter of the shaft, i.e. four elements and ®ve nodes, is used in the computation. c An eighth of the shaft, i.e. four elements and six nodes, is used in the computation. d A quarter of the shaft, i.e. four elements and nine nodes, is used in the computation. b

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

accurate result is provided by the hybrid stress element BH6 developed in the present work. Yamada et al. [17,18] also computed the torsion of the equilateral triangular shaft shown in Fig. 2(a) using a ®ne mesh. Due to symmetry, only one-sixth of the triangle was discretized in their computations using 120 elements and 79 nodes, as shown in Fig. 11±13(a) in Ref. [20]. In our computations by comparison, the entire triangle is divided into only 48 elements and 61 nodes, as shown in Fig. 2(c). Comparison of our results plotted in Figs. 3 and 4 with their results listed in Tables 11±1 and 11±2, and plotted in Figs. 11±13(b) and 11±16 in Ref. [20], show that the values of warping obtained by the displacement triangular element, the hybrid stress triangular element, and hybrid elements H4 and H4 + BH6 are close to the exact solution. The accuracy of the shear stresses computed by the hybrid stress triangular element and H4 + BH6 is

561

similar and better than that given by the displacement method. The advantage of using the higher order elements H4 and BH6 together, with the latter satisfying exactly the traction-free condition on the boundary, over other elements is that accurate results can be obtained with a coarse mesh. A.2. Mesh distortion sensitivity The square shaft shown in Fig. A1(a) is employed to test the mesh distortion sensitivity of the present method based on the H4 and BH6 elements. The coordinates of point A give the distortion e of the mesh relative to a regular mesh, shown by the broken lines. For convenience, the modulus of rigidity is assumed to be Cs = 1.0, and the torque applied at each end of the shaft is Mt = 1.0. The variations of the computed torsional rigidity and of the non-vanishing component of the shear stresses tyz at point B with the distorted

Fig. A1. Model for test on mesh distortion sensitivity. (





) Exact solution. (r) Results by H4. (R) Results by BH6. (a) Geometry and ®nite element mesh. (b) Variation of the computed torsional rigidity with the distorted distance e. (c) Variation of tyz at point B with the distorted distance e.

562

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563

Fig. A2. Convergence of the computed torsional rigidity and shear stress with mesh re®nement. (





) Exact solution. (r) Results by H4. (R) Results by H4 + BH6. (a) The computed torsional rigidity of the equilateral triangular shaft. (b) tyz at point D in Fig. 2(c). (c) The computed torsional rigidity of the square shaft. (d) tyz at point 7 in Fig. 6(b).

distance e are plotted in Figs. A1(b) and (c). Numerical results show clearly that both H4 and BH6 are insensitive to mesh distortion. In fact, in examples 1 and 4 in the body of the text, distorted elements have been used in the computation. A.3. Test of convergence with h-re®nement In order to test the pattern of convergence of the energy error norm with h-re®nement, shafts of an equilateral triangular and of a square section discussed in examples 1 and 2, respectively, are revisited. The computed torsional rigidity and the non-vanishing component of the shear stresses tyz at point D in Fig. 2(c) and at point 7 in Fig. 6(b) vary with the re®nement of mesh, or equivalently, with an increase in the number of nodes. The variations are plotted in Fig. A2. Numerical results show clearly that not only the torsional rigidity (i.e. the energy of the system), but also the stresses converge rapidly with h-re®nement. Generally speaking, the use of H4 together with BH6 gives better convergence than the use of H4 alone.

References [1] Timoshenko SP, Goodier JN. Theory of elasticity. 3rd ed. New York: McGraw-Hill, 1970. p. 291±353. [2] Novozhilov VV. Theory of elasticity. London: Pergamon, 1961. p. 282±335. [3] Courant R. Variational methods for the solution of problems of equilibrium and vibration. Bull Am Math Soc 1943;49:1±23. [4] Zienkiewicz OC, Taylor RL. The ®nite element method, 4th ed. New York: McGraw-Hill, 1989. [5] Cook RD. Concepts and applications of ®nite element analysis, 2nd ed. New York: Wiley, 1981. [6] Bathe KJ. Finite element procedures in engineering analysis. Englewood Cli€s, NJ: Prentice-Hall, 1982. p. 420±423. [7] Wu CC, Pian THH. Incompatible numerical analysis and hybrid element method. Beijing: Science Press, 1997. [8] Herrmann LR. Elastic torsional analysis of irregular shapes. J Engng Mech Div, ASCE 1965;91:11±19. [9] Kawai T, Yoshimura N. Finite element analysis on the torsion of a bar with uniform cross-section. Seisan Kenkyu 1968;20:246±8. . [10] Kruhula JL, Lauterbach GF. A ®nite element solution for Saint Venant torsion. AIAA J 1969;7:2200±3.

Q.Z. Xiao et al. / Computers and Structures 71 (1999) 535±563 [11] Synge JL, Cahill WF. The torsion of a hollow square. Q Appl Math 1957;55:217±24. [12] Zienkiewicz OC, Cheung YK. Finite elements in the solution of ®eld problems. Engineer 1965;220:507±10. [13] Feldmann JM. Direct numerical determination of stresses in elastic solids as illustrated by the torsion problem. Int J Solids Struct 1968;4:675±88. [14] Moan T. Finite element stress ®eld solution of the problem of Saint Venant torsion. Int J Numer Meth Engng 1973;5:455±8. [15] Golley BW, Cashman JD. The torsional sti€ness of prismatic members with sections divisible into rectangles. Int J Numer Meth Engng 1974;8:671±5. [16] Heise U. A ®nite element analysis in polar co-ordinates of the Saint Venant torsion problem. Int J Numer Meth Engng 1974;8:713±29. [17] Yamada Y, Nakagiri S, Takatsuka K. Analysis of SaintVenant torsion problem by a hybrid stress model. In: Proceedings of US±Japan Seminar on Matrix Methods of Structural Analysis and Design, Tokyo, 1969 August. [18] Yamada Y, Nakagiri S, Takatsuka K. Elastic±plastic analysis of Saint Venant torsion problem by a hybrid stress model. Int J Numer Meth Engng 1972;5:193±207. [19] Noor AK, Andersen CM. Mixed isoparameteric elements for Saint Venant torsion. Comput Meth Appl Mech Engng 1975;6:195±218. [20] Desai CS. Elementary ®nite element method. Englewood Cli€s, NJ: Prentice-Hall, 1979. [21] Jaswon MA, Ponter AR. An integral equation solution of the torsion problem. Proc R Soc Lond 1963;273A:237±46. [22] Mukherjee S, Morjaria M. Comparison of boundary element and ®nite element methods in the inelastic torsion

[23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

563

of prismatic shafts. Int J Numer Meth Engng 1981;17:1576±88. Zielinski AP, Zienkiewicz OC. Generalized ®nite element analysis with T-complete boundary solution functions. Int J Numer Meth Engng 1985;21:509±28. Chen CN. The warping torsion bar model of the di€erential quadrature element method. Comput Struct 1998;66:249±57. Lam SSE. Application of the di€erential quadrature method to two-dimensional problems with arbitrary geometry. Comput Struct 1993;47:459±64. Zhong HZ. Elastic torsional analysis of prismatic shafts by di€erential quadrature method. Commun Num Meth Engng 1998;14:195±208. Zhong HZ, He YH. Solution of Poisson and Laplace equations by quadrilateral quadrature element. Int J Solids Struct 1998;35:2805±19. Pian THH. Derivation of element sti€ness matrices by assumed stress distributions. AIAA J 1964;2:1333±6. Wu CC, Bu¯er H. Multivariable ®nite elements: consistency and optimization. Sci China (A) 1991;34:284±99. Pian THH, Chen DP. Alternative ways for formulation of hybrid stress elements. Int J Numer Meth Engng 1982;18:1679±84. Washizu K. Variational methods in elasticity and plasticity, 2nd ed.. New York: Pergamon, 1975. Karihaloo BL, Hemp WS. Optimum shapes for given torsional and bending rigidities. Proc R Soc Lond 1987;A409:67±77. Xiao QZ, Karihaloo BL, Williams FW, Application of improved penalty-equilibrium hybrid stress element method to crack problems, Engng Fract. Mech. (in press).