The influence of stress concentration on creep rupture at non-stationary loading

The influence of stress concentration on creep rupture at non-stationary loading

Int. Z mech. Sci., Vol. 19, pp. 521-531. Pergamon Press 1977. Printed in Great Britain THE I N F L U E N C E OF STRESS CONCENTRATION ON CREEP R U P ...

525KB Sizes 0 Downloads 48 Views

Int. Z mech. Sci., Vol. 19, pp. 521-531.

Pergamon Press 1977. Printed in Great Britain

THE I N F L U E N C E OF STRESS CONCENTRATION ON CREEP R U P T U R E AT N O N - S T A T I O N A R Y LOADING N. N. MALININ and A. A. NIGUIN Department of Strength of Materials and of Dynamics and Strength of Machines, Moscow Higher Technical School, Moscow, 2-ya Baumanskaya ul.5, U.S.S.R. (Received 17 Nonember 1975; and in revised form 15 February 1977) Abstract--A method is given for the calculation of creep rupture strength of parts containing stress concentrations. The creep theory takes account of the damage and anisotropy of materials. which arise during deformation. For the example of a plane specimen with a notch, results of calculations are compared with experimental data. The results of an investigation into the creep strength of the fir-tree root of a turbine blade at non-stationary loading is given. NOTATION

stress tensor strain tensor ~pt proportional limit E modulus of elasticity /z Poisson ratio epl = ~ra/E strain corresponding to a proportional limit reduced stress tensor e,j = * , d ~ , reduced strain tensor reduced stress deviator /3 .2 equivalent reduced stress

d,j

reduced reduced reduced reduced

active stress deviator at plasticity additional stress deviator at plasticity active stress deviator at creep additional stress deviator at creep

equivalent reduced active stress at creep "3

.1t2 equivalent reduced additional stress at creep #,

#

2

. ~n to

reduced average reduced reduced reduced thermal

normal stress in the plane with the normal v reduced normal stress elastic strain tensor plastic strain tensor creep strain tensor expansion

equivalent reduced creep strain tensor transgranular damage intergranular damage vector 1

&(e.)

1

a,(l + Ct#,')

P ( ~ o ) = ao + C o e o + Co#o2

a, C, Bo De, H, Fc, c, m, q. K,. K2, a~. C1, n, a0, Co, Co constants for the material in question T absolute temperature t time x, y, z axes of rectangular coordinates with an origine in a centre of gravity of an element {u} matrix of displacements of some point {8} matrix of displacements of nodes of an element A area of an element [D] matrix of mechanical properties {40} matrix of initial strains {p} matrix of distributed load {F} matrix of external contour forces [G] stiffness matrix 8o Kronecker delta AI area of section MS Vol. 19, No. 9-~B

~21

522

N . N . MALININ and A. A. NIGU1N 1. INTRODUCTION

COMPLICATED conditions of w o r k i n g of m a c h i n e parts with d e m a n d s f o r i n c r e a s e d

service life require more precise methods of strength calculation. The investigation of rupture in stress concentration zones is of considerable practical importance. When it is necessary to determine the service life of a component the basic question becomes not the calculation of stresses but the estimation of a rupture time. As most parts operate in non-stationary conditions, it is important, as is shown in Ref (1), to take into account the anisotropy arising during plastic and creep deformation. 2. BASIC EQUATIONS OF THE THEORY OF CREEP WITH ANISOTROPIC HARDENING AND DAMAGE Reduced stresses and strains are used, and it is assumed that when these variables are used, stress-strain curves at different temperatures coincide and there is a single plastic deformation equation of the form ~p = ~0(O). It is supposed that the processes of plasticity and creep are independent. For a description of instantaneous plastic deformation we use the theory of plasticitty with kinematic hardening. Suppose that the stress deviator is divided into two deviators: the active stress deviator and the additional stress deviator

~,~= ~j + p,j,

(1)

and during plastic flow the yield surface displaces as a rigid body. This displacement is characterized by the additional stress deviator. Such a model describes the Bauschinger effect. Results of the experiment show that the deformation anisotropy is not removed for closed strain paths. The dependence between the additional stress deviator and the plastic strain increment tensor is expressed

by~ 2

dp,j = ~ A, (de) d~,~,.

(2)

As the hardening is kinematic we have:

~,,d~,~ = 0. From this condition, using the associated flow rule and the dependence of (2), we obtain

Considering creep we shall apply a theory of creep with anisotropic hardening and damage. In Ref. (1) the basic principles of the theory of creep with anisotropic hardening were formulated. Bailey's idea that the first and second stages of creep may be explained on the basis of the interaction between strain hardening and annealing has been used? In accordance with this idea, strain hardening prevails in the first stage of creep and therefore the creep rate decreases and anisotropy arises. In the second stage the equilibrium between strain hardening and annealing is established and the creep rate is constant. The application of the theory of creep with anisotropic hardening makes it possible to describe the I~haviour of materials subject to variable loads more precisely than by means of the ordinary creep theories. To represent tertiary creep it is necessary to introduce an additional structural parameter which is a scalar quantity and measures the transgranular damage. In Ref. (4) to describe the transgranular damage an equation of the following form was offered: do, = P (6"o)d~..

(4)

Also the parameter o, was introduced into the creep equations in the form: 4

d~j~ = Bc exp (-TH) (l + co,m)sh d-~2d. 3-~dj 2 / H\ m X. Xi~ d~,, = ~Ao(~.)d~,,c - Oc exp ~,-y)(l + co, ) s h E 7 " dr.

(5)

Thus the transgranular damage arising during the first and second stages begins to influence creep and in the third stage the creep rate increases sharply. Besides this, during all the three stages of creep, intergranular damage arises. As cracking on grain boundaries occurs more readily in the presence of tensile stresses and for complex loading principal planes can rotate, it is natural to assume that intergranular damage is a vector quantity. A value of this vector varies from 0 at the initial moment of time to 1 at rupture. Supposing that the accumulation of intergranular

The influence of stress concentration on creep rupture at non-stationary loading

523

damage takes place during deformation of the material and neglecting the influence of grain boundary cracking on creep, a formula for the vector of the intergranular damage was obtained: 4 Id¢,[ =

qDcexp

---

(1 + c~ m) exp (Kl6, + K22,) dt

It is supposed that when the parameter of the transgranular damage or the vector of the intergranular damage at some point in an arbitrary direction reaches a critical value an instantaneous rupture takes place. The movement of the front of a crack is not investigated. In describing the results of creep rupture tests of materials under multiaxial states of stress the criteria of a maximum shear stress and a maximum principal stress are used as usual. Both criteria can be taken from the equation (6) if we neglect the third stage of creep. In more general cases the damage at a given temperature is a function of 6~, 6 , 60 and i.~. Thus the equations (2) to (6) describe plastic and creep deformation as well as creep rupture at non-stationary loading. In order to use these equations it is necessary to have only an instantaneous stress-strain curve, creep curves and creep rupture curves. 3. SOLUTION OF P L A N E STRAIN PROBLEM It is convenient to use the finite element method for investigations of the state of stress in components of complex shapeJ We propose to use only triangular elements. Such elements have frequently been used for solving plastic and creep problems. In this case we can assume that strains and stresses within each element are constant and therefore the onset of plastic flow there is no necessity to determine a boundary between elastic and plastic regions within the element. Let us consider the plane strain problem at constant axial strain. For a trian4gular plane element with the nodes i, k, k displacement increments at an arbitrary point of the element can be expressed in terms of increments of nodal displacements in the following way /d~,i

n . t /, {du}={d~}ffi[N]{dS}=[[Ni][NI][N~']]J/ ddSj'

/dS / [dSk,,J

where

[Ni]=a~+c~x+c~y[l ~] 2A

L0

'

and 2

c~ -- Yi - yk,

a~ = ~ a ,

cj = xk - xi.

The connection between strain increments and increments of displacements of nodes of an element are obtained from the previous equation in the form: [d,x ]'

{de} = ~d~, ~ - [Bl{dS} =

[[B,I[BII[B~II{dS},

(7)

[dTx,J where

s l[c, c~o]

[ l]=~[O

t.c,

c~ "

Then stress increments in an element will be given by {do'} =

Idor"l do', = [D] ({d~} - {d~.t). [d~z,J

(8)

In the generalcase, disiributed load may act on an element and external forces may be applied to the nodes, thus f

y

and { ~ =

I F

524

N.N.

MALININ and A. A. NIGUIN

From equilibrium conditions we obtain [G] {dS} = {dF} + {dF}p + {dF},e,

(9)

where f {dF}p = JA [N]T {dp } dA,

f {dF}, e= JA [B]r [D] {dEo} d~,

are the matrices of external force increments, equivalent distributed load and initial strain increments and [G] --- fA [B]r [D] [B] dA. is the stiffness matrix. The equation (9) is the system of linear algebraic equations for increments of displacements of nodes. Knowing them it is possible from expressions (7) and (8) to determine strain and stress increments. Then the basic problem is that of determining the matrices [D] and (d~o}, which are defined by the mechanical properties of the material. The tensor of total strain increments can be presented in the form: d0 d ~ = d ~ + d~ip + d~i ~ + - - , ept where the elastic strain increment tensor, taking into account that Ept and It depend on temperature, is given by d~q. = ( l + i t ) [ d O o _ 8 ~ j

L

3it d 6 o + /

l+it

3it \dEpl] /#,j - a , j - - a - o / - - / + \ l+it / ep~J

('~o - 3a,#o) dit.

For plane strain when de33 is a constant 1 d~33 = ~ [(1t -- P1133) d # . + (it - P2233)d@22-- 2p1233d~12 - de330 + d~33]. I + P3333 where { 9 ~

Piiu= 4

- -

SijSkt when when

S,=I

and

Sqd~>0

S, < l

or

Sii d~i <~O,

and dept dO dg33o = de33c + [#33 - It(flu + #z2)] - - - (~Ii + ~22) dit + - - . ept %t Then the physical equation, connecting strain and stress increments are.

,,

d~,yj

[d~,yj

Ld~,yo)

,

L~,y)

where [A] is the symmetric matrix (3 × 3) with elements: (it - p.33) 2 - - , 1 + P3333

A.=l+p..

Aa2 = A21 =

-

It + p..~

(it - P m~) (it - P~33) 1 + P3333

l

A , 3 = A31 = 2 [p,,,2 + (it - PlI33).______PP1233 1 + P3333 J ' A22 = 1+P2222

(lit -- P2233)2 - - . 1 + P3333

A23= A32 = 2 [p,222 + ( i t - P2233)P,233] 1 + P3333 J' p2233 ) . A33 = 2 (1 + It + 2Pro2 - 2 l + P3333/

(1o)

The influence of stress concentration on creep rupture at non-stationary loading

525

Elements of matrices {d~0} and {#} are given by, depj dO ~ -- P1133d~33o, d~x0 = d~H~ + [6.H - ]J, (6.22 + 6.33)] - - - (6.22+ 6.33) d/~ + - - + Epl P3333

Epl l +

dO /t - P2233j d~yo = d~22~+ [6.22- / ~ (6.11+ 6"3;)]depl _ (6.11+ 6.33)d/~ + - - + .-------- u~330, • pj opt 1 + P3333

d3'x'°=2Fd~nc+(l+lt)6.12dePl+6.12dlt-L Pl233 j_ ] ~pl 1 + P3333 0e330]" -- P1133 ~'x

~J" -- P2233

epl(1-t-P3333)'

eY

2pt233

¢pt(l+p3333) '

ex'= ~pl( 1 +P3333)

Matrices contained in the expression (8) can be determined from the equation (10) in the following way: [D] =

E[A] -~,

{de0}]= epa{d~0} and

{e0}2 = ¢pl{~}.

To determine d¢33 it is necessary to solve twice the system of equation (9), the first time calculating the displacement increments {d8h from external forces {dF} and {dp} and from initial strain increments {de0h, and the second time computing displacements {8}2 from initial strains {e0}2. Having determined from the expressions (7) and (8) {d~r}mand {~}2 and supposing there to be an absence of external force in the direction of axis z, we obtain

fA,(do'~dAl )l d¢33

f~. (¢rz)2d~j Then the total stress increments in the solid are given by, {do~} = {do'}! + d¢33{0"}2.

In solving the problem it is necessary to divide the total working life of a component into a series of time intervals during which stresses and strains change insignificantly. Then replacing infinitely small increments of functions by finite values we obtain the change of the state of stress and strain for a given interval of time. The right-hand part of equation (9) for a given regime of loading can usually be written in the form: {dF} + {dF}p + {dF}~0 =

{~}

dr.

In this case the dependence of the stress increments on the time increment is linear. Therefore it is possible not to take beforehand a time interval and to do the calculation for a single interval. Further, its value will be determined from the condition that a maximum stress increment will not exceed a given value (0.02-0.05~rp~). This makes it possible to decrease the number of calculation steps while keeping the necessary accuracy. 4. A CREEP STRENGTH CALCULATION FOR A SPECIMENS WITH STRESS CONCENTRATION WITH A COMPARISON OF THEORETICAL AND E X P E R I M E N T A L DATA The method described above was used for calculations for plane specimens with a notch made from a Nimonic at constant load and loading in steps at a temperature of 700~C. The shape of specimens is shown in Fig. 1. The thickness of plane specimens were given such values that the places of stress concentration

.~_~/R0.36 FIG. 1. A notched specimen.

526

N . N . I~L~JJNIN and A. A. NIGUIN

were in a plane strain condition. For this it was necessary that the thickness of a specimen was 2.5-3.0 times greater than the least width.6 Therefore the thickness of the specimens was made 9 mm. Stress-strain, creep and creep rupture curves, obtained by experiments, are depicted in Figs. 2-4. In Fig. 4 experimental and theoretical creep rupture curves of specimens with a notch are also presented. For notched specimens the nominal stress is the ratio of a tensile force to the area of the minimum cross-section. These stresses received theoretically and experimentally for the same time of rupture differ by no more than 6%. Such a difference corresponds to the scatter typical of experimental data. For stepped loading two cases were considered: the increasing of nominal stresses from 226 to 294 MNIm e and decreasing them from 294 to 245 MNIm ~. The loading of a specimen was altered when the time of the test was equal to one half of the rupture time at initial load, i.e. 30 hr at 226 MN/m 2 and 4 hr at 294 MN/m 2. 1200

8O0

o

4OO

0

4

8

12

E, %

FIG. 2. Stress-strain diagram. (~6

0.4

0.2

I I0

! 20

[ 50

I 40

50

t, hr

FIG. 3. Creep curves at 700°C. 600

200 ~

o I

2

I

4

I I0

I 20

J 40

I I00

J 200

400

t, hr

experiment (O smooth specimens, O notched FIG. 4. Creep-rupture curves at 700~C: specimens), - . . . . calculation (notched specimens).

The influence of stress concentration on creep rupture at non-stationary loading

527

Diagrams of normal stresses in the least section of a specimen, drawn as a result of the calculation for stepped loading in the first reghne at different moments in time is depicted in Fig. 5. Immediately after initial loading the largest stress appears at the root of a notch. Later during the process of creep the peaks of the stresses smooth down somewhat and are displaced into the inner part of a specimen, where the first cracks are initiated. Experimental and theoretical rupture times for stepped loading are given in Table 1. As follows from Table 1, the agreement of the results in the first case are better than in the second case. It can be explained that the ratio of the time of the second step of loading to the whole time of the test in the first case is considerably smaller than in the second case. As the time of the first step of loading in the experimental and theoretical investigations are the same it is evedent that the difference in the time of rupture in the first case must be less than in the second case. It follows from Fig. 4 that the slope of the creep rupture curves with respect to the time axis is small and therefore the difference in the rupture times even by a factor of four corresponds to a difference in rupture stresses of only 10-15%. As far as the rupture time of the specimen of complicated shape under the non-stationary loading is considered, the difference between theoretical and experimental date is satisfactory. It follows from this comparison that the methods described above can be used for creep calculations of machine components with stress concentrations subjected to non-stationary loading. 5. CALCULATIONS FOR THE FIR-TREE ROOT OF A TURBINE BLADE The fir-tree root is one of the most highly stressed elements operating at high temperature and stress concentration. Usually in calculations on fir-tree roots two problems are considered: an investigation of the distribution of tooth forces without taking into account stress concentration and then an evaluation of stress concentration for a given distribution of forces.~'s The method is correct for stresses below the elastic limit but strictly speaking is not true over the elastic limit and with creep taking place. The calculation for fir-tree roots within the elastic limit by the finite element method has been previously reported. 9 Simultaneously the contact problem was solved and stress concentration was investigated. In the present article, taking into account stress concentration, the plastic and creep contact problem for the fir-tree root is solved. In previous work 9 it was shown that the distribution of tooth forces and stresses does not significantly depend upon friction forces arising at the contact surfaces on the teeth. Therefore friction in the contact zones will not be taken into consideration. This permits the formulation of a symmetric stiffness matrix, thus decreasing to a considerable extent the volume of calculation. Let us consider two neighbouring elements of a disc and a blade (Fig. 6). As the elements can displace

TABLE I Stresses

(MN/m 2)

orI = 226;

N specimen t experimental" hours

~r, = 294

1

2

3

45.5

38.6

42.0

t theoretical: hours

err = 294; ~r, = 245 4

1

2

3 6 - 5 1 2 - 5 18.3

33.0

3

4

8.9

20.0

44.0

, 2k

]o

0.8

L32.9

0.4

I~ 0

i

0.4.

I

x,0.8mrn

I

1.2

t

~32.9 t ~'z t,30 'o3O29 ~'x I 1.6

FIG. 5. Diagrams of normals stresses in the least section of a notched specimen at loading in step (increasing nominal stresses after 30 hr from 226 to 294 MN/m2).

528

N.N.

MALININ and A. A. NtGUtN

k

/

/ L,~

X

FIG. 6. Contact elements. freely along the contact line it is possible to write down the following expressions for increments o! displacements of nodes in places where contact occurs: (dSi, - dB.y) cos ¢ = (dS~ - dS,~) sin ¢, and (dSjy - dS,~,) cos ¢ = (d/~/~ - dS~,) sin ¢.

(11)

Using equation (9), we obtain the increment of the contact force in the nodes i and n, acting perpendicular to the contact line, dP/n ffi

1 cos
[dFnyp + dFn~o - (G.y,,x d6,,x + G,,y.y dSn~,+ O,,y,,,xdBmx + Gnymy dS,ny + On,t~ dSt~ + anyty dSly)].

,

The force P~. is external for the elements of a disc and a blade, considered separately. Substituting dP~ instead of d F in the equation (9) written for the projection on the disc and the root of a blade, using the dependences (11) and putting together the equations of equilibrium of the elements of a disc and of a blade, we obtain the new stiffness matrix [G] and the matrix of the increments of external forces {dF} for the whole connection. The line and the column, for which in matrix [G] the element G.y.y was a diagonal one is absent in the matrix [G]. The c o ~ e s p o n d i n g element in the matrix {dF} is also absent. The remainder of the terms of the matrices [G] and {dF}, which change by comparison with the elements of the matrices [G] and {dF} in the result taking into account the contact of the nodes i and n will be given by, G~

= G,.~ + tg2 ~. G.,.,

G~j, = G ~ y - tg ,p - G.y,.y,

G',i~ = G,,~ - tg ~ "

G.,.,.

dt~r,~ = -- tg q~(G.y,,~ + tg ~ . G.ym,), d,,,.~ = G.ymx + tg q~" G.y,.y,

d~.= = G~.x + tg ~ ( G~.y + G.,.x + tg cp " G.,.y ), d,~l~ = G,~t, + tg cp- G.,t~, O,~,y = G,~t, + tg q~" G.y,,, d/~= = d F ~ + dF,~ o- tg ¢ (dF.y, + dF.~,), d~y = dFiyo + d F ~ o+ dF.,p + dFn~,

dP~ = dF~p + d F ~ , + tg q~(dF.yp + d F ~ ) . Having obtained dependences in the investigation of the contact of the nodes ] and m it can be shown that the matrix [(~] remains symmetric.

The influence of stress concentration on creep rupture at non-stationary loading

529

Thus, having obtained the matrices of stiffness and matrices of increments of external load separately for the projection on the disc and the root of a blade which by definition take into account plasticity and creep, it is possible to construct the full matrices of stiffness and the matrices of increments of external load for the complete body of a fir-tree root. In the process of the calculation it is necessary after each time step to determine the forces in all the nodes which are in contact and to calculate the displacements of the nodes, which are able to enter into contact with each other. If the force in some contacting node is negative then it is necessary to consider these nodes independently of one another at the next step of calculation. This method of calulation allows the creep contact problem to be solved when the contact area is unknown beforehand. In such a way it is possible to determine what teeth or what parts of teeth are in contact at any moment of time as well as yielding a distribution of forces on the teeth or on the sides of a tooth. For example, the fir-tree root of a turbine blade shown in Fig. 7 was considered. The alternation of rotational velocity of the disc and the temperature at the points A, B, C, D in the loading cycle are depicted in Fig. 8. In the calculations it is assumed that the temperature both in the root of a blade and in the projection on the disc varies linearly with distance from the axis of rotation. As shown by calculations the force distribution in a root, and also stresses and creep strength, depend to a considerable extent on the initial clearances between the teeth. Figs. 9 and 10 show the alternation of forces on fir-tree root teeth and maximum circumferential stresses (in relative coordinates) in the fillets of the teeth of a disc and a blade for two initial cycles of loading, for the condition of a close fit between all teeth before commencement of working. As shown in Fig. 10 maximum stresses arise in the fillet Dj of the projection on the disc. This place defines the creep strength of the whole connection. For unloading the teeth I of a root it is necessary to provide initial clearances between the teeth ID and IB and between the teeth lid and liB. With initial clearances between the teeth ID and Is and between the teeth lid and lib creep strength can increase seven times and rupture takes place in the fillet Dm of the projection on the disc. With initial clearances between the teeth lid and lib and between the teeth IIID and IIIB creep strength of a fir-tree root can decrease by a factor of twelve as compared with the root without initial clearances. It follows from Fig. 10 that a blade root is stressed to a considerable extent less than the projection on the disc. It may be explained that the force at point F is more than at point E (see Fig. 7). Therefore having provided differently inclined contact planes EF of a blade root and the projection on the disc at their

B

F

F

FIG. 7. Fir-tree root.

530

N.N.

IV[AUNIN and A. A. NIGUIN 800

/

,? T8

.oIV

/Te

n

F

600

f

x

T

I0

II

I

o

io

40OO

E

200O/

"~

I

,

400

k.r

200

i

20 30 t, min FIG. 8. Cycle of loading of a fir-tree root.

0.-

--

40

/

~ ~"-~ Ir

E

tO

20

30

40

50

t , min

FIG. 9. D/stribution of fir-tree root tooth forces at the a b s e n c e of initial clearances.

Dz Dn 0.8

E Dm D~ 0.4

Bz

I 10

I 20

30

40

50

t, min

FIG. 10. M a x i m u m s t r e s s e s in fir-tree root tooth fillets at the a b s e n c e of initial clearances.

contact in point E, it is possible at the e x p e n s e of the loading of a blade root to unload the projection on the disc, i.e. to p r o m o t e a more equal strength construction. In a given case in the provision of an initial angle of opening at the contact zone of the tooth pairs, the creep strength of a fir-tree root can be increased s e v e n times. T h u s , having calculated a fir-tree root at different initial clearances a m o n g teeth and angles o f opening of a joint, it is possible to c h o o s e s u c h tolerances for the d i m e n s i o n s of a root as will best suit a given service life.

The influence of stress concentration on creep rupture at non-stationary loading

531

6. CONCLUSION

A worked out method of creep strength calculation allows one to evaluate the rupture time of parts with stress concentrations at non-stationary loading. In the present article this method is applied to calculations for a plane specimen with a notch and to the fir-tree roots of turbine blades. Theoretical and experimental data of the creep strength of the plane specimens with a notch agree satisfactorily. On the basis of the calculation for the fir-tree root recommendations about tolerances for dimensions of the joint are given. The basic relations obtained in the present work may be used for calculations for any components working under creep conditions. REFERENCES N. N. MALImN and G. M. KHArm~SKY, Int. J. mech. Sci. 14, 235 (1972). R. A. ARUTYUNYAN and A. A. VAKULENKO, IZV. Akad. Nauk SSSR. Mech. 4, 53 (1965). R. W. BAILEY, Proc. Inst. mech. Engrs. 131, 131 (1935). G. M. KHADJINSKY, IZV. Akad. Nauk SSSR. M T T 6, 29 (1971). O. C. ZmNKIEWICZ, The Finite Element Method in Engineering Science. McGraw-Hill, London (1971). V. S. JouKowsIo, In Problems of Strength in Mechanical Engineering (in Russian), Akad. Nauk S S S R 2, 54 (1959). 7. A. N. GRUBIN, Strength Calculations of Fir-tree Root of Turbine Blades (in Russian), Machinostroeni~ 0970). 8. A. A. NIC;UIN, lzv. Vuz. Mach. 11, 38 (1972). 9. S. K. CHAN and I. S. TUBA, Int. J. mech. Sci. 13, 7 (1971).

1. 2. 3. 4. 5. 6.