Engineering Fracture Mechmics Vol. 34, No. 516,pp. 1097-l 107, 1989 Printed in Great Britain.
CRACK BRANCHING: ON MACROSCOPIC
0013-7944/89 $3.00 + 0.00 0 1989 Pergamon Press pk.
A “TWIN-CRACK” MODEL BASED ENERGY FRACTURE CRITERIA
P. S. THEGCARIS, N. P. ANDRIANOPOULOS
and S. K. KOURKOULIS
Department of Engineering Science, Section of Mechanics, Athens National Technical University, P.O. Box 77230, Athens GrI75-10, Greece Abatraet-An improved approach for the prediction of the branching angle of running cracks is presented. It is based on the assumption that from the early steps of propagation, ahead of the running main crack-tip a tuft of microcracks is formed, from which two symmetric microcracks prevail over the others at the final steps before macroscopic branching. Further expansion of these two microcracks, which can be predicted by existing macroscopic fracture criteria, rest&s in macroscopic branching. The unknown initial orientation of the two dominant microcracks represents the stochastic part of the whole branching phenomenon, while the main crack geometry, velocity and external load refer to the deterministic part. The results obtained by applying the T-criterion for the expected branching angle show a satisfactory agreement with existing experimental data.
1. INTRODUCTION IT IS WELL-ESTABLISHED that existing macroscopic fracture criteria can adequately predict the conditions for crack initiation in the case of static loading. Both experimentally measured quantities, i.e. the critical fracture load and the expected angle of crack propagation, agree in general with theoretical predictions. This is not true in the case of dynamic loading. Experimental data, especially for the angle of crack-branching, deviate considerably from the corresponding theoretical prediction@-51. Usually, this disagreement is addressed to the unstable nature of the crack propagation, absent in static loadings, which, then, can be affected by small-scale irregularities of the material, variations of its mechanical properties, etc. To the authors’ opinion an additional and important factor should be put with all previously mentioned factors. Indeed, in dynamic problems the fracture criteria are usually applied to a single running crack and it is supposed that the crack bifurcates when a suitable mechanical quantity exhibits two separate extreme values, the direction of which coincides with the expected direction of each crack-branch. This approach was firstly introduced by Yoffe[6]. However, it is proved experimentally that, at the tip of a running ductile crack, a tuft of microcracks is always presentf7,8]. This cloud of microdefects consists of many microcracks of various lengths and directions. Indeed, extensive experimental evidence with slowly running cracks under mode I loading observed in the waning electron microscope has revealed that for ductile cracking in polycarbonate and metals the procedure of blunting and kink initiation followed a sequence of blunting of the crack tip to form a flat crack front rounded-off at the corners followed by the creation of either an asymmetric kink at the one extremity of the crack front or a double rather unsymmetric kink of parallel microcracks. Figure l(c) of Ref.[7] indicates clearly the one kink situation, whereas Fig. l(a) of this paper shows the two-kink case. As the dynamic loading precedes, the crack front is stuffed with several microcracks and kinks as shown in Fig. l(b). It is then reasonable to assume that, at the final steps before macroscopic crack propagation two microcracks predominate over all the others and run almost parallel with roughly the same velocity at least for some small distance, thus not violating the macroscopic impression of a single running crack with extremely irregular lips. Hence, the problem of crack-branching is to describe and predict the conditions necessary to cause mutual deviation of the two predominant microcracks, lying on both sides of the main crack axis, to obtain macroscopic crack-branching. Then, for ductile fracture, any criterion is preferable to apply in a model with two running singularities, instead of a single one. This concept was recognized as a reasonable way to escape the inaccuracies arised from the Y~~-~j~e application of various fracture criteria, and there exist 1097
1098
P. S. THEOCARIS
et al.
either experimental[9] or theoretical[9-I I] recent approaches, where crack branching is studied in pre-branched crack geometries. However, the above-mentioned studies are limited to the point of locating geometries, where the initial crack-branches are expected to expand straightforward[9] or the mode II stress intensity factor is zero[9-111. In both cases it is implied that the branching criterion is K,, = 0, as it was, initially, proposed by Bilby et aL[9]. If it is correct, one could expect that the condition K,, = 0 must be valid in case of predicting the angle of propagation of a single inclined crack loaded statically. However, the variety of existing fracture criteria convinces that such a simple and smart criterion is nonexistent and those existing are only rough approximations. Consequently, the necessity of applying a more capable fracture criterion in case of crack branching is obvious. Then, it will be shown in the sequel that the alternative approach of the twig-croci model is a convenient method to solve the problem of branching and it is superior when compared to the Yo@-mode. Thus, it constitutes an improvement of a similar technique[l2,13], applied to cases of already bifurcated cracks. For the application of the twin-crack model the corresponding stress-field around the crack tip must be known. The values of the dynamic stress intensity factors, Kr and Z&, for either branch of the crack are imputed separately by applying a suitable n~e~cal method[ 11,141 and using the appropriate correction factors for the dynamic case. Then, these values were introduced into the dynamic stress field expressions given in[15, 161. The numerical method incorporated in this analysis has been proved stable and accurate, especially for extreme geometries, which are of great importance for the present case. The aim of this paper is to predict the branching angle, by using the above-referenced numerical method and the twig-crack model in the macroscopic fracture criteria, and to compare the predictions obtained with existing experimental data. 2. THEORETICAL
ANALYSIS
2.1. The stress $eld The results published in[l 1, 141 are used for the crack configuration shown in Fig. l(c). This configuration consists of a main crack, OA, and two branches, OB and OC, at different angles with respect to the main-crack axis. The thin infinite plate containing the cracks is loaded at infinity by a static load, normal to the main-crack axis, under generalized plane stress conditions. The values of stress intensity factors at the three tips, A, B and C, are computed by solving numerically a system of three complex singular integral equations with Cauchy-type kernels with 17 integration points. In Fig. 2 the values of & = K,,, and K1,,B= - K,,, normalized to Ki,* of the main crack are plotted vs branch inclination c;o= (ps = - qc for two crack-length ratios b/a = c/a = 0.25 and b/a = c/a = 0.025. It is seen from these figures that both quantities for both geometries possess maximum values at small angles cp. Consequently, at these small angles the stress field is stronger than at any other neighbou~ng angle. Also, Kr,,s and Ku, change sign at ip w 11” for b/a = 0.25 and at cp N 17.5” for b/u = 0.025. This zeroing of K,, implies that there exists a geometry of cracks at which the local stress field is symmetric with respect to the axis of each crack branch, In addition, at angles lower than those of this local symmetry, K,,s and K,,, are bigger than K,,*, the stress intensity factor of the tip A (Fig. 2). This superiority of K,,* and Kr,c over Kr,, implies that in crack configurations like those of Fig. 1 crack propagation is easier in the form of a do~bie s~gu~urity (tips B and C), than in the form of a single one (tip A), the same being not true for large angles (p and lengths 6, c[12]. Finally, it must be noted that the value of the angle of local symmetry, where &r,B = J&c = 0, depends, only, on the crack geometry and not on its velocity. 2.2. Description of the twin-crack mo&l approach At the final steps before macroscopic branching, it is assumed that two microcracks from the cloud surrounding the running main crack tip predominate over all the other@]. These two cracks may be approximated in the model to have the same length (b = c), and the same inclination (cp,,= - cpc) with respect to main crack axis (Fig. lc). For arbitrary (but small) values of the lengths and inclinations of the microcracks, the respective values of stress intensity factors, K,,, = K,,c and
Crack branching
A
Fig. l.(a), (b) The initial steps of mi~o-bmnchi~g.
(c) The respective crack ~n~gumtion.
1099
Crack branching
1101
-0
i-0 I
Fig. 2. Values of stress intensity factors for the tips B, C, reduced to mode I stress intensity factor of tip A for mode I and mode II.
- Knc, are computed, as previously described. With the values of stress intensity factors = known, the stress field for the running microcracks is computed by means of the well known Freund solution in the form given in Refs[ 17, 181. Finally, a fracture criterion can be applied to predict the expected angle of propagation, O,, for each branch. Due to symmetry, the expected angles of propagation, O,, are equal and of opposite signs for the two branches, Referring to Fig. 3(a), the angle, &, of the crack branching is given by:
K 11,B
&, = cp + sin-’
(l-b2/12 sin2 Q1’* - f cos 0,
11
as it can be obtained from elementary calculations. In real cases, where 1>>b = c, eq. (1) reduces to:
(al
(bl
Fig. 3.(a) Definition of symbols I, 8,, @,,.(b) The two regions of branch survival or coalescence.
(1)
P. S. THECXARIS
1102
et al.
In the sequel, the simplified eq. (2) will be used to define the angle of crack branching. At this point the following remarks are worth making: (i) The value of branching angle, 8,, depends on 1, i.e. it depends on the distance from the crack-tip, 0, where all measurements are performed[5]. (ii) Depending on the relative magnitudes of 40and f$, angle tfbobtains values which may result to an early absorption of the expanding branch by the main crack, which is then disappearing. Namely, from eq. (2>, B,, vanishes when cp + 0, = 0. In this case the two branches run parallel to each other. When $ becomes negative (u, + 0, < 0), it leads to an early coalescence of the two branches (Fig. 3b). These unsuccessful branching attempts may be correlated with the commonly observed hackle zones, appearing just before a successful branching[$ 81. (iii) Macroscopic branching will occur when the characteristic for each fracture criterion quantity reaches a critical level, C,, depending on the material, provided that the existing couple of microcracks has the appropriate orientation (angle 40) to avoid coalescence. For a given stress at infinity, G,, this critical branching condition can be achieved by many combinations of microcrack lengths b, orientations tp and velocities u. Hence, a function f (rr,, b/a, 9, Y) = C, governs the branching phenomenon. The exact form of this function depends on the fracture criterion assumed as valid. However, the twin-crackmodel implies that there could not be a unique critical quantity, like the stress intensity factor or the crack velocity, responsible for branching. Instead, a critical surface in the 4-dimensional space of ((Jo, b/a, q, V) is characterizing the phenomenon. Each one of the above four quantities has a low limit, necessary, but not sufficient to cause alone branching. This remark agrees well with experimental evidence[fi]. It may be useless to repeat that according to Yoff-mode it is assumed, contrary to the present concept, that a critical crack velocity (or, equivalently, a critical stress intensity factor) is sufficient to cause branching. 2.3. Application of the twig-crack model The above described twin-crackmodel was applied to three welLknown macroscopic fracture criteria, namely the o,-criterion[l9], the S-criterion [20] and the T_criterion[21]. In this application, a crack-length ratio b/a = c/a = 0.025 was selected. It was checked that for ratios lower than 0.025 the angle, OP,of branch propagation remains practically constant, the same being also true for the other characteristic angles of the problem (Fig. 4). Thus, the predictions for the branching angles, 6,, are settled, implying that for branch lengths b = c less than 0.025 of the main crack length, a, the lengths of the two microcracks do not affect strongly the branching angle. In Fig. 5 the polar distribution of the values of o,, S and Tv is plotted for an initial branch inclination rp = - 7” and a crack velocity Y = 0.5c,, where c2 is the shear-wave velocity. It is seen bla
(=ctal --m
Fig. 4. Influence of the ratio {b/a) on the characteristic angles of ~irn~ k; , local s~metry, and branch propagation.
Fig. 5. Polar distribution of u#e-, S-, T’quantities cp = -7” and v -0.5~~.
for
Craclc branching
1103
in this figure that ce and TV exhibit a maximum, whereas S presents a minimum, at the direction where the branch is expected to propagate. These extreme values are plotted in Fig. 6 vs angle, cp, of initial branch orientation, for the same crack velocity u = 0.5~~.It is seen in this figure that the S-criterion shows a maximum minimorum at (o e 17.5” and, hence, two microcracks having an initial orientation f 17.5” are prone to continue running and to cause macroscopic branching. The T-criterion has a local rna~rn~ at angle Q,55:7” and a global one at CJJ x 40” and the oOcriterion, also, presents a local maximum at v, x 11” and a global one at cpx 48”. Hence, the two latter criteria predict two optimum microcrack orientations for branching, instead of the unique one predicted by the S-criterion. However, the global maximum corresponds to states of an early branch absorption and hence it may be discarded. The two maxima of TV and a, are separated by a relative local minimum value, which in the case of the T-criterion lines at cpx 18” and in the case of the a,-criterion at Q x 17.5’ (Fig. 6). Figure 7 presents the influence of the crack velocity on the behaviour of the extreme values of the dilatational strain energy density, TV, for three different crack velocities v = 0.35c,, OSc, and o = 0.65~~.The following remarks are valid: (i) The global maximum remains constant and independent of the crack velocity. Its value equals: 2 (1 - 2v) GA4*=~~Ql!
(3)
where E and v are the elastic modulus and the Poisson ratio and a, is the yield stress of the material. (ii) The global maximum appears at bigger q-angles as the crack velocity increases. (iii) The local maxima and minima remain eseentially at the same location (40x constant). (iv) The magnitudes of these two local extrema decrease rapidly as the crack velocity increases, but with different rates. Namely, for ti = 0.5 c2 the local maximum is only 1% higher than the t.00
0.98
0.92
Fig. 6. The extreme values of Q-, S-, Tvquantities vs branch inclination 0 reduced to the respective maximummaximorum value for v = 0.5 c,.
Fig. 7. Them extreme values of TVvs bran& inclination q reduced to the constant maximum-maximorum value, for crack velocitiesv = 0.35 c,, v = 0.50 c,, v = 0.65 c,.
1104
P. S. THEOCARIS
et al.
1.20
0.98
t 2 0.94 E
A Global maximum 0 x
Local maximum Local minimum
0.b
Fig. 8. Variation of the three extrema of Tv vs reduced crack velocity.
corresponding minimum and for u = 0.65 c2 it becomes 4.5% higher. Then, as the crack velocity increases the local maximum becomes sharper. A cumulative presentation of the influence of the crack velocity on the magnitude of the above three extrema of TV is given in Fig. 8. It is evident in this figure that the minimum velocity to engender branching is Omi,k: 0.2 c2, where a slight difference between the three extreme values begins to emerge. For any initial branch orientation, cp, there exists according to all fracture criteria an extreme value of their respective characteristic quantity (ae, S or TV), whose direction coincides with the expected branch expansion. The direction, 8,, of branch propagation is plotted in Fig. 9, vs branch
M-
01;
t
o^ d -2o-
-601
I
I 60
9. Angle of branch propagation vs initial branch inclination cp according to a,-, S-, T-criteria.
Fig. 10. Angle of branch propagation vs initial branch inclination cp,according to Tsriterion, for crack velocities v=0.35c2, v=0.SOc,,v=0.65c2.
Crack branching
0.15
1105
0.65
0.40 v/c2 -
Fig. IL Branching angie f&, vs reduced crack velocity, according to the Yq@-rno& and twin-crack model versions of T-criterion, and experimental data.
orientation, tp, for the three above-mentioned criteria and for a representative crack velocity v =osc,. It is seen in this figure that the angle of branch propagation, e,,, may take either positive or negative values, depending on cp. Eventual branching occurs when a coalescence between the main crack and its branches is avoided, i.e. when cp + 0, > 0 (eq. 2). In the specific case of Fig. 9, all three curves 0, =f(p) lie above the limiting line of cp + 6, = 0, for 9 > 3” and then branching is possible for any branch orientation satisfying this condition. From all the orientations, q, two are the most probable, according to the og and the T-criteria, and one according to the S-criterion. These most probable orientations are coinciding with the respective maxima shown in Fig. 6. Thus, for v = 0.5 c,, according to the a,-criterion, the optimum branch orientations are 9 = 11’ and cp = 48” (Fig. 6), and the respective angles of propagation are from Fig. 9, 0, = 14” and 0, = - 46”, resulting in a branc~ng angle, according to eq. (Z), equal t0 eb = 25” or 6),,= 2”. Similarly, in the case of T-criterion we have cp = 7” or cp = 40”, and 0, = 16” or 8, = - 38”, resulting in e,, = 23” or &, = 2”. Finally, the S-criterion gives cp = 17.5”, 0, = 0” and hence eb = 17.5”. Figure 10 presents the influence of crack velocity, v, on the value of the angle 0, of branch propagation according to the T-criterion. This angle is absolutely increasing with crack velocity in an accelerating manner. It is worth mentioning that, for v = 0.65 c,, parts of the respective curve 0, =f(rp) lie below the line cp + 19,= 0, implying that branching isunsuccessful. In addition, all the three curves of Fig. 10, pass through the same point of local symmetry (rp, 0,) = (17.5”, 0’) where K,, = 0. This point, as it has been already mentioned, depends only on the crack configuration and it is independent of the crack velocity. Figure 11 shows the prediction of T-criterion for the expected branching angle, e,,, vs crack velocity, v. In this figure experimental points have been taken from Refs[l-51. For comparison purposes, the predictions for the same quantity according to the Yofl-mode of the T-criterion are, also, given. In the same figure the twin-crack predictions of the T-criterion are represented by three curves, corresponding to initial branch inclinations cp = 5”, cp = 7” and cp = 11”. From these three curves the middle one (for cp = 7”) corresponds to the local maximum value of TV (see Fig. 7), whereas the other two cases correspond to a value of Tv slightly lesser than its maximum value. The three curves of T-criterion predictions constitute a gradually narrowing band, as crack velocity increases. This narrowing band is clear in Fig. 12, where it was assumed that the branches may propagate even when TV has a value only 0.05% lesser than its maximum value TV,,,,=. In this case two angles, t9,, of branch propagation are defined from Fig. 7, corresponding to the same TV-value (TV = 0.9~5Tv.~~) on both sides of TV,,,. The difference, A@,, between the two extreme values EFM 34.5/6-c
P. S. THEOCARIS
1106
0.15
et al.
0.65
0.40
090
v/c2 +
Fig. 12. Width of the divergence zone of propagation angle 0,, for TV= 0.9995Tv,,, velocity.
vs reduced crack
of 0, varies from A(_$, x 12” for v z 0.2 c2, to A8, < 0.5” for v > 0.7 c2. Consequently, a strong scattering in 6, (and similarly in 0,) values may be caused by slight variations in the Tv-values, especially in low velocities. Returning to Fig. 11, one may conclude that the twin-crack model agrees better with experiments, as compared to the Yofi-mode. For example experimental points of Ref.[l] are unexplained by the Yofl-mode, introduced in any of the three criteria[5]. 3. CONCLUSIONS An improved model for the branching phenomenon of cracks under dynamic propagation is described in the present paper. This model reflects the experimental fact that at the vicinity of the tip of a running crack a cloud or tuft of microcracks always exists in dynamic fractures. To overcome theoretical complexities in the computation of the respective stress fields it was assumed that in the model only two equal microcracks, symmetrically placed to the main crack axis, are always dominating all other microcracks in the tuft. Application of a fracture criterion in either of the new two branch-crack tips permits the evaluation of the expected angle of the microcrack propagation. The two microcracks may progressively deviate from each other or they coalesce, depending on the geometry of the crack zone and the crack velocity. In the former case the macroscopic phenomenon of crack branching is achieved. Consequently, there do not exist definite branching points. The crack from the early steps of its propagation is in a state of microscopic branching. Furthermore, the lack of sharp branching conditions implies that the common assumption of a single critical quantity, necessary for branching to occur is not realistic. It seems much more reasonable to accept the assumption of a critical branching surface, in the form of a function f(a, , b/u. q, v) = C,, depending on the type of the stress field, the crack geometry and the crack velocity. The parameters of this function can be classified in two groups. In the first group belong the applied stress, cm, and the crack velocity, v, which are definite macroscopic quantities, representing the deterministic character of the phenomenon. In the second group belong the microcrack parameters b/a and cp, which represent the stochastic part of the phenomenon of branching. It was proved that from the latter two parameters the microcrack length, b/a, is of minor importance, concerning the value of the macroscopic branching angle. On the contrary, the initial branch orientation, cp, affects strongly the branching angle, not only by means of eqs (1) or (2) but, also, by introducing an uncertainty component in the value of 0,, as it is shown in Fig. 12. The formation of microcracks at the vicinity of the tip of a running crack depends on the real state of stress of the non-homogeneous and non-isotropic, imperfect material and, consequently, the value of angle, cp, can be a good representative of this state, since microcracks are sensitive to local irregularities. Anyhow, angle cp is restricted by the state of stress. It cannot be very small, because in this case the two microcracks coalesce, nor it can be very large, because in such geometries the local stress field weakens and a single crack is more prone to branch.
Crack branching
1107
Hence, a band of acceptable values of the q-angle exists, a fact which means that the branching phenomenon is not completely stochastic and, on the other hand, that the expected angles of crack branching cannot be uniquely defined to have a definite value, but, rather, they belong to a band of predictions. This band of branching-angle predictions has a decreasing width as the macroscopic parameters (i.e. velocity Y, etc.) are increasing. Really, the stochastic character of branching dominates when its deterministic counterpart is weak and vice versa. This is clearly shown in Fig. 12. Finally, experimental data gathered from existing literature show a good agreement with the twin-crack model presented here. This agreement supports the hypothesis that an initial asymmetric orientation of microcracks may be responsible for the macroscopic crack kinking and arrest, observed in many types of fractures in reality. Acknowledgemenrs-The results of the research programme presented in this paper constitute a part of the Dr Eng. thesis, submitted to the department of Engng Science of the Athens Nat. Technical University by S. K.,Kourkoulis. This research program was supported by the University Research Fund.
REFERENCES VI A. S. Kobayashi, B. G. Wade, W. B. Bradley and S. T. Chiu, Crack branching in Homalite-100 sheets. Engng Fracture Me& 6, 81-92 (1974).
PI F. P. Bullen, F. Henderson and H. L. Wain, Crack branching in heavily drawn chromium. Phil. Mug. 21, 689-699 (1970). [31 F. P. Bowden, J. H. Brunton, J. E. Field and A. P. Heyes, Controlled fracture of brittle solids and interruption of electric current. Nature, Lond. 216, 3842 (1967). 141 A. J. Carlsson, On the mechanism of brittle fracture propagation. Trans. R. Insr. Technol. Stockholm, Sweden, 805, l-39 (1963). PI P. S. Theocaris, N. P. Andrianopoulos and S. K. Kourkoulis, Brittle curving and branching under high dynamic loading. Res Mechunicu, submitted. WI E. H. Yoffe, The moving Griffith crack. Phil. Mug. 42, 739-750 (1951). (71 P. S. Theocaris, Blunting phenomena in cracked ductile plates under mixed mode conditions. Engng Frucrure Me&. 31, 255-270 (1988).
181 K. Ravi Chandar and W. G. Knauss, An experimental investigation into dynamic fracture: II Microstructural aspects. Inr. J. Fracture 26, 65-80 (1984).
[91 J. F. Kalthoff, On the propagation direction of bifurcated cracks, in Dynamic Crack Propagation (Edited by G. C. Sih), pp. 449-458. Noordhoff, Leyden (1973). See also: B. A. Bilby, G. E. Cardew and I. C. Howard, Stress intensity factors at the tips of kinked and forked cracks. Fracture, 3, 197-200 (1977). J. P. Dempsey and P. Burgers, Dynamic crack branching in brittle solids. Inr. J. Fracture 27, 203-213 (1985). it;; P. S. Theocaris, Asymmetric branching of cracks. J. uppf. Mech. 44, 61 l-618 (1977). VI N. P. Andrianopoulos, Propagation predictions in bifurcated cracks by energy methods. 1st U.S.A.-Greece Symposium of Mixed Mode Crack Propagation, pp. 99-108 (1980). iI31 P. S. Theocaris and N. P. Andrianopoulos, Propagation predictions for branched cracks by the strain energy density criterion. J. Mech. phys. Solids 30, 23-36 (1982). P. S. Theocaris, Branched cracks at small angles (An Addendum). Engng Fracture Mech. 17, 361-366 (1982). it:; L. B. Freund, Crack propagation in an elastic solid subjected to general loading-I. Constant rate of extension. J. Mech. Phys. Solids 20, 129-140 (1972).
WI L. B. Freund, Crack propagation in an elastic solid subjected to general loading-II. Nonuniform rate of extension. J. Mech. Phys. Solids 20, 141-152 (1972). u71 M. Ramulu and A. S. Kobayashi, Strain energy density fracture criterion in elastodynamic mixed mode crack propagation. Ennnn Fracture Mech. 18, 1087-1098 (1988). WI M. Ramulu and A. S. Kobayashi, Strain energy density fracture criterion in elastodynamic mixed mode crack propagation (Errata). Engng Fracture Mech. 19, 583 (1984). iI91 F. Erdogan and G. C. Sih, On the crack extension in plates under plane loading and transverse shear. J. Bus. Engng 8SD, 519-527 (1963). [20] G. C. Sih, Strain energy density factor applied to mixed-mode crack problems. Inf. J. Fracture Mech. 10, 305-321 (1974). [21] P. S. Theocaris and N. P. Andrianopoulos, The Griffith-Growan fracture theory revisited: The T-Criterion. Inr. J. Mech. Sci. 27, 783-801 (1985). (Received 26 January 1989)