The influence of cutting edge radius on the machined surface of brittle materials in simulated orthogonal machining

The influence of cutting edge radius on the machined surface of brittle materials in simulated orthogonal machining

Inl. J. Msch. Tools M~rmfact. Vol. 36. No. 9. pp. 971-983. 1996 Copyright © 1996 Elsevier Science Lid Printed in Groat Britain. All rigbt~ reserved 08...

596KB Sizes 0 Downloads 36 Views

Inl. J. Msch. Tools M~rmfact. Vol. 36. No. 9. pp. 971-983. 1996 Copyright © 1996 Elsevier Science Lid Printed in Groat Britain. All rigbt~ reserved 0890--6955/96,515.00 + .00

Pergamon

PII: S0890-6955(96)00019.3

THE INFLUENCE OF CUTTING EDGE RADIUS ON THE M A C H I N E D S U R F A C E O F B R I T T L E M A T E R I A L S IN S I M U L A T E D ORTHOGONAL

MACHINING

JOO-HYUN KIM* (Original received 3 October 1995; in final form I 1 January 1996)

Abstract--When cutting brittle materials, many surface cracks appear on the fiat surface due to tensile stresses induced by the cutting tool. These cracks propagate into the substrate and will affect the fracture properties of a body. Crack spacings and the directions of crack propagation into glass were calculated numerically by applying the finite element method and using linear-elastic fracture mechanics. The calculated crack spacings were in the range of the experimental results. Stress intensity factors and crack extension angles were dependent on the edge radius and the load on the tool, and from these two factors the possible directions of crack propagation were calculated. The calculated propagation directions were in good agreement with the actual crack propagation directions. Copyright © 1996 Elsevier Science Ltd

I.

INTRODUCTION

Although difficult to machine, it is becoming more feasible to cut most of the brittle materials in current use. Brittle materials such as glass and ceramics have been successfully cut and high quality surfaces have been obtained at small depths of cut. However, there has been very little research on the surface formation mechanisms of these materials in cutting. This is mainly because the mechanism of the material removal process has not yet been sufficiently clarified. Many researchers have been interested in the theoretical mechanics of machining which cannot be applied to the study of machined surface. Complete modelling of the mechanics of cutting for the purpose of analysing machined surface is probably not possible at this time. However, a universal event and a characteristic of cutting is that the flank face or heel of the cutting tool slides along the finished surface, imposing very high stresses upon that surface. The modelling problem can be simplified by studying only this sliding or rubbing action, ignoring the sequence and pathway of strain and fracture ahead of the cutting tool and the machine characteristics. The work that follows assumes the geometry of a cylinder sliding in a direction perpendicular to its axis on a flat plate. There is ample evidence that new tools have rounded rather than perfectly sharp edges, and that worn tools slide against the finished surface on a narrow flat surface referred to as a wear land on the flank face of the tool. A cylinder sliding on a flat surface has a more amenable geometry for calculation purposes when compared to a flat wear land. In the following work the cylinder will not cut, but rather it will only slide along the flat surface. The cutting action and the nature of the approaching surface in real cutting is too complicated for calculation purposes at this time. In two-dimensional geometry sliding, the imposed normal and tangential stresses can be taken to be distributed elliptically as shown in Fig. 1. A tensile stress is induced in the flat surface on either side of a slider due to the normal load, but an additional tensile stress arises behind the slider from the tangential stress (Fig. 2) [ 1]. This stress is frequently sufficient enough to cause fracture on the flat surface behind the slider. The important point of the existence of such cracks is that they extend into the flat surface at a sufficient depth to influence the fracture properties of the body. Periodic surface cracks often appear in the flat surface of a brittle material behind a heavily loaded slider. Neither the spacing between the cracks nor the depth of these cracks

*National Industrial Technology Institute, 2 Joongang, Kwacheon, Kyunggi, 427-010, Korea. 971

Joo-Hyun Kim

972

sliding direction

slider

N ~

fiat

---~

-q~

qo

>

--'-'~

surface

Fig. I. The load distribution of a semi-infiniteelastic plate in two-dimensional sliding contact (Po = maximum normal stress, q(~= maximum tangential stress, 2a = width of contact).

i

T

T

1

. . . .

~

~

.i~lo

"

"

.]

--al~"q

gJ t',,'<....7...6~

~

/?t

I

//

1

";

1

[

1

1

\"

t ]

.alo

i

i

~

i

i

~

i

i

i

Fig. 2. Contours of constant Om~lpo due to surface traction (om,~= maximum principal stress, Po = maximum normal stress).

appear to have been the subject of analysis, particularly to see if these two quantities might be related. Work of this type is described in this paper. The stress in a flat surface behind a slider was calculated using forces and material properties that would result in a crack. The direction and rate of crack propagation as the slider advanced was then determined by fracture mechanics. At some point the crack ceased to grow because of a weakening stress field whereas the tensile stress behind the slider increased until a new crack began. Crack propagation can be analysed by calculating stress intensity factors at the crack tip. Several different methods have been used to calculate the stress intensity factors of an existing vertical crack in sliding. Keer et al. [2] used the numerical method, Jahanmir et al. [3] used the photoelastic method, and Kimura and Shima [4] used the finite element method. Because all of the above works considered only a single existing crack, a full description of the crack path and the spacings between cracks were not analysed. Dubourg and Villechaise [5] and Dubourg et al. [6] developed a semi-analytical model of multiple fatigue cracks in studies of the mutual interactions between cracks but did not analyse the crack path down into the subsurface. In order to completely analyse crack propagation, stress intensity factors at the incrementing crack tip must be used. Mathematical analysis of an entire fracture field in which there are incremental changes of direction is difficult or impossible because of the complex geometry. Therefore a finite element analysis was used to calculate the stress intensity factors for this complex geometry and stress system.

Influence of Cutting Edge Radius on the Machined Surface of Brittle Materials

973

2. EXPERIMENTALWORK Sliding tests were performed using cylindrical tool tips of several different radii in order to study the influence of the radius on crack pattern productions. The cylindrical surface of the tool tip was pressed and slid in a direction perpendicular to the axis of the cylinder. All tests were done with a constant depth of indentation of 0.05 mm. Abrasive paper was used to produce the desired tip radius and the radius of each edge was measured with an optical microscope. The tool tips were frequently reground to maintain the desired edge radius during sliding. Glass was selected as a good approximation to a homogeneous, isotropic and linearelastic ideal brittle material in order that it was possible to compare experimental results with the numerical analysis. Small glass plates were attached to aluminum alloy plates with epoxy, and then each aluminum alloy plate was clamped to the vice on the shaper. The width of the cylindrical tool tip was 2 mm as shown in Fig. 3. The vertical and horizontal components of the sliding force were measured by strain-gauged transducers. The sliding speed was very low, 8 m rain-~, so that the heating from the sliding action would not contribute to any localized heating of the glass. Periodic surface cracks on the glass were observed as shown in Fig. 4 and these were tool (slider) attached with epoxy

Fig. 3. Sliding test on glass.

sliding direction Fig. 4. Periodic surface cracks on glass (edge radius -- 0.44 ram, sliding speed = 8 m min-~).

974

Joo-HyunKim

measured by an optical microscope. Table 1 shows the sliding force components and crack spacings for sliders of various edge radius, where F, is the component normal to the glass surface and F, is the component tangential to it. These sliding force components will be used in numerically analysing the crack spacing on the surface and the crack path down into the subsurface. After sliding of the tool, the glass plates were detached from the aluminum alloy plates by heating sufficiently to soften the epoxy, but not high enough to affect the glass. Then the glass plates were broken along the sliding direction by scoring the underside of the glass plate with a glass cutter. The cross-section of the glass was observed by an optical microscope. The direction of crack propagation in the glass depended on the edge radius of the tool when maintaining constant depth of indentation, because the different edge radii produced different normal and tangential forces, which affect the stress pattern under the surface. Figure 5 shows three types of crack propagation of the surface cracks where the cracks propagate: in the direction of sliding (Type A); in the direction opposite to the sliding direction (Type B); or straight and vertically (Type C). In these experiments, some edges produced both Type A and B while others always produced Type B. For example, the cross-sections of glass produced by an edge radius of 0.23 mm showed both Type A and Table I. Loads and crack spacings for each radius of tool tip Edge radius (ram)

F, (N)

F, (N)

Crack spacing (I.tm)

130 175 191 204 258 345

67 66 80 92 131 155

3-9 4-15 6-22 9-24 15-28 37-55

0.05 0.11 0.16 0.23 0.44 0.76

\\

\ (Type A) f

// (TypeB) v

(Type C) Fig. 5. Three types of crack propagation.

Influence of Cutting Edge Radius on the Machined Surface of Brittle Materials

975

B crack propagation, as observed in Figs 6 and 7, but cross-sections produced by an edge radius of 0.76 mm showed only Type B crack propagation, as observed in Fig. 8. 3.

FINITE ELEMENT MODELLING

When the loads move over the fiat surface, there is a large tensile stress behind the slider, which is a maximum at the surface and decays with depth [1]. However, once a

sliding direction Fig. 6. Cross-section of glass showing crack propagation in the same direction as the sliding (edge radius = 0.23 ram).

sliding direction Fig. 7. Cross-section of glass showing crack propagation in the direction opposite to sliding (edge radius : 0.23 mm).

976

Joo-Hyun Kim

sliding direction Fig. 8. Cross-section of glass shov, ing crack propagation in the direction opposite to sliding (edge radius = 0.76 mm).

crack forms the stress diminishes between the back edge of the slider and the crack. This stress increases as the slider moves away from the crack, but in a manner or rate that depends on the stress field in the substrate, which in turn depends on the depth and direction of the propagated crack. Once the crack propagates a particular distance, it will stop as the load moves away. Then a new crack initiates at the next site of high tensile stress on the surface behind the slider. In this model, a vertical crack from the surface was initially assumed in order to apply fracture mechanics. The normal and tangential forces developed by the slider motion lead to a mixture of mode I and II crack extensions of the assumed vertical crack, which may cause the initial crack to propagate at some angle other than vertical, as shown in Fig. 9. Stress intensity factors were calculated for the zone at the crack tip, from which the direction of the next increment of cracking may be estimated. Further increments of crack propagation may be analysed by attempting to follow a path where the stress intensity factors remain greater than the critical stress intensity factor, Kc, of the material. Measurements were made of the Kc of glass, and it was found to cover a range of values from 4.5 to 5.5 MPa m ~2. These values, when inserted into the analysis, produced a range of possible crack extension directions. r

A

d!

Fig. 9. Crack propagation with crack extension angle(d = initial depth of crack, 0~ = crack extension angle, b = distance between slider and crack, a = half width of contact).

Influenceof Cutting Edge Radiuson the MachinedSurfaceof Brittle Materials

977

The analysis began with an assumed crack depth as well as direction. Figure 10 shows the flow chart for determining the depth of the initial crack. If the possible crack extension angle at the initial crack tip produced by the analysis was 5 ° or less, a longer crack was assumed in order to avoid tedious new mesh generation. The choice of an acceptable angle (e.g. 5 ° ) surely introduced an uncertainty in the results, but this uncertainly was considerably less than that which resulted from the range of values of Kc for glass taken from the tests. If the next calculated angle exceeded 30 ° , the initial crack depth was abandoned and a shorter one was assumed, on the principle that directions were not likely to change drastically during propagation. The real cracks from experiments did not show a sudden change of the direction for crack extension, but rather a gradual change, thus supporting this idea. This process was also used to determine the length of extended cracks (Fig. 10).

4. NUMERICALPROCEDURE In combined modes I and II loading, the stresses around the crack tip can be derived as 1 0 [ 203 1 ~o = ~,/2Itr cos ~ Kicos ~ - ~ Kn sin 0 x,~ -

1

~

2~/2~r

cos

0

(1)

[Klsin0 + KH(3cos0 - 1)]

where r, 0 is a right-handed polar coordinate system located at the crack tip and 0 = 0 is the continuation of the line of the crack into the material. The stress intensity factors, /(i and Kn, for modes I and II were obtained by the crack opening displacement method coupled with quarter-point elements [7, 8].

Assumethe initial depth of crack

~-~

calculationofthe angle,0c

Assumelongercrack if[0e[ ~5 T

~

[

No

Assumesiorter crack ff[0e[ 230 I

No I Yes

( Assume the lengthofthe extendedcrack [

vl-Icalculationofthe angle,Ocrk I

Assumelongerlength J if leeI <5

I

[

Assumeshorterlength if loci >30

No

T

Fig. 10. Flowchart for the prediction of crack propagation.

978

Joo-Hyun Kim

The maximum principal stress criterion [9] states that crack growth will occur in the direction perpendicular to the maximum principal stress. This criterion implies that the mode II stress intensity factor becomes zero as the crack extends. From Equation (1) the stress Oo is a principal stress if "t~o--0. Therefore, the crack extension angle 0¢ can be found from Kx sin 0~ + K,(3 cos 0c - 1) = 0

(2)

which gives

/(KII2

0¢ 1 K I 1 tan 2 - 4 K . + 4 V \ ~ ]

+8

(3)

Once the new crack is initiated at some angle from the tip of the old crack, it will extend when the stress intensity factors at the tip of a small branch exceed the critical value. These stress intensity factors K,* and K." can be calculated in terms of the stress intensity factors KI and K. of the main crack and the crack extension angles 0c, by the following equations [I 0]: =

c, ,K, + c,:K,i

and K~. = c2tKi + c22K.

(4)

where c,,=~

3cos~-+cos

c,2-

4

sin~-+sin-~¢

1(0c

c2, = ~ sin ~- + sin and c:2=~

cos~-+3cos

If K~" exceeds the critical stress intensity factor, the crack will propagate with the calculated crack extension angle (Fig. 9). It is obvious that K." becomes zero because 0¢ is derived from the maximum principal stress criterion implying that the crack follows a path with zero value of mode II stress intensity factor. With the extended crack, the next possible crack propagation at the new crack tip can be estimated. Sliding was simulated by translating the distribution of Hertzian normal and tangential stresses (Fig. 1) over a finite element mesh. The finite element mesh had eight-node quadrilateral elements for a semi-infinite elastic two-dimensional solid with an initial vertical crack. The finite element solutions were obtained using the multi-purpose finite element package ABAQUS [l 1]. Material for this application was glass with elastic modulus = 70 GPa and Poisson's ratio = 0.17.

Influence of Cutting Edge Radius on the Machined Surface of Brittle Materials 5.

5.1.

979

RESULTS A N D DISCUSSION

Crack spacing

In order to study the influence of the depth of the initial vertical crack on the initiation of a new surface crack, finite element meshes with vertical cracks of different depths were generated. The variations of maximum principal stress as the load moves away from the vertical cracks for different initial depths are shown in Fig. 11 when r = 0.23 mm and in Fig. 12 when r = 0.76 mm, in which r = edge radius of tool tip. In these figures, the vertical crack is at b/a = 1.0 where b and a are as defined in Fig. 9. Calculated crack spacings were expected to produce a range of values because of the uncertainty in the calculation of the final crack depth and the angle of extension. However, the crack spacing was found to be strongly dependent on crack depth only for the smaller ranges of crack depth as shown in Figs 11 and 12. When r = 0.23 mm, an initial crack deeper than 40 I.tm appears not to influence the initiation of a new crack, and for the case of r = 0.76 mm, an initial crack deeper than 80 lam does not appear to influence the initiation of a new crack. This critical depth remained uncertain until after some experimental observations. Then it was found that all of the cracks were much longer than the critical crack length. This explained why the crack spacings were much more uniform than the depth and angle of cracks. If the breaking stress of these glasses is known, an estimate can be made of the crack spacing between periodic surface cracks from the figures. Using three-point bending tests,

3°I 2.5

~

~

~

-~

2.0

i

1.5

o

d~lO

o 1.0

=

0.6

0.0 1.0

1.4

1.8

2.2

2.6

3.0

3.4

3.8

4.2

4.6

b/a

Fig. 1 I. M a x i m u m principal stress on surface behind the slider for r = 0.23 mm(d = initial depth of crack). 2.5

2.0

¢h

1.5

1.0

0.5

0.0 1.0

1.4

1.8

2.2

2.6

3.0

3.4

3.8

4.2

4.6

b/a Fig. 12. M a x i m u m principal stress on surface behind the slider for r = 0.76 m m (d = initial depth of crack).

980

loo-Hyun Kim

the breaking stress of these glasses was found to be 110 + 20 MPa. These values are very small compared with the values in Figs I l and 12, because flaws are always present in glass and the fracture of glass in a bending test originates at the largest flaw within a wide region. However, in sliding the high stresses behind a slider are more localized and a high stress is required to initiate a new crack. The cracks will initiate when the maximum tensile stress behind the slider reaches the tensile strength of glass when there are only very small flaws. The breaking stress of glass fibres increases rapidly as the size of the fibre is decreased due to the decreasing probability of the occurrence of suitable flaws. In this sense the strength of glass with very small flaws may have similar values to the large breaking stress of fine fibres. The breaking stress of fine fibres is known to be in the range 0.07-2.8 GPa [12]. The results in Figs I l and 12 can be rationalized by using the experimental results from Table l: a breaking stress of 1.5 GPa results in reasonable crack spacings. With this value of breaking stress and considering that all real cracks are deeper than 80 I.tm, it can be seen from Figs i 1 and 12 that b/a = 1.5 for r = 0.23 ram, and b/a = 2.0 for r = 0.76 ram. Therefore the crack spacing (b - a) = l0 I.tm for r = 0.23 ram, and 48 lira for r = 0.76 mm. These values are in the range of crack spacing for each radius as shown in Table I.

5.2. Directions of crack propagation K~ and K . were calculated by the crack opening displacement method coupled with quarter-point elements for six different radii of cylindrical tool tips as the load on the surface of glass moved away from the initial vertical crack. From K~ and K,, 0c, and Kx" were calculated by Equations (3) and (4), and are shown in Figs 13 and 14, respectively. In these figures the initial depth of crack, d, was assumed to be 5 ~m. With an initial depth of crack of 5 ~tm, the calculated values of K~" for all edges was relatively large when the load began to move. If K~" exceeds the critical stress intensity factor of glass at the crack tip, the crack will propagate. The load position at which the crack will extend, b, may be found from Fig. 14 for a given value of the stress intensity factor: the value of b/a at which Kj" first exceeds the estimate of critical stress intensity factor describes the position of crack propagation. The crack extension angle may then be determined for this position from Fig. 13. Yet because of the uncertainty in the value of Kc for glass due to the fact that the actual value of Kc of a material may vary with position, it is difficult to estimate the exact load position at which the crack will extend. However, it was expected that all crack propagations would be completed by the time the slider had moved to a position where K~', within the range of critical stress intensity factors, was a maximum. That is, no further propagation would take place as K," decreased. Figure 15 shows the depth of the initial crack and the possible directions of crack extension, estimated from the method in Fig. 10, for r = 0.23 and 0.76 ram. In order to 10

0

-10

-20 x -30 a

r-0.23 r-0,44 r,,0,76

-40

-50 1.0

Fig.

I

I

1

1

1

1

I

1

1.4

1.8

2.2

2.6

3.0

3.4

3.8

,1.2

4.6

b/a 13. Crack extension angle for different edge radii (initial depth of crack d = 5 gin).

Influence of Cutting Edge Radius on the Machined Surface of Brittle Materials

981

0.40 0.35 0.30 ea

o

0.25 0.20 0.15 0.10

X

0.05 1.0

r.0.76

J

I

I

1

I

1

1

J

1.4

1.8

2.2

2.6

3.0

3.4

3.8

4.2

4.6

bla Fig. 14. Kn" for different edge radii (initial depth of crack d = 5 Bm). r = 0.23 mm

Id

= 5 Bm

Oe = - 10 to 10 deg

r = 0.76 m m

d = 10 Bm

0¢ = - 5 dog Fig. 15. Reasonable initial depth of crack and possible directions of crack extension.

study further propagations of the cracks, new finite element meshes for the extended cracks were generated, where the planes of the extended cracks were simply rotated around the tip of the old crack to make the desired crack extension angles. Because the load position at which the initial vertical crack extends was not known, the load was moved from the original position, that is, at b/a = 1. The process used for the determination of the initial depth of crack was applied to determine the length of the extended crack. For the two cases shown in Fig. 15, the following further crack propagations were expected. (a) r = 0.23 mm and d = 5 Bm. The possible angles, 0c, for this case were found to be in the range of - 1 0 to I0 ° with d = 5 lam. New meshes having the extended crack with

982

Joo-Hyun Kim

the different angles of - 10, - 5 , 0, 5, and 10° were generated. As an example, considering the angle of - 10° where the length of the extended crack was estimated to be 5 p.m, the range of the next possible crack extension at the new crack tip was from bla = 1.0 to 1.6. In this range 0~ varied from - 3 0 to 0 °. In this way the next possible range of crack extension angles at the tip of the extended crack was calculated. (b) r = 0.76 mm and d = 10 p.m. The only possible 0c for this case was found to be - 5 ° with d = 10 p.m. By using new meshes having the extended crack with the angle of - 5 ° where the length of the extended crack was estimated to be 10 p.m, 0c and Kj ° values at the new crack tip were calculated. It was found that the next possible angle of 0c for the extended crack varied from - 2 2 to 5 ° . The crack propagation behaviour based on the calculated results for the two cases discussed above are shown in Fig. 16. In the case of r = 0.23 mm, the crack may propagate in several different directions; that is, some cracks propagate in the direction of sliding, but others in a direction of opposite to the sliding. However, in the case of r = 0.76 mm, the directions of crack propagation were always opposite to the sliding direction. These results were in very good agreement with the experimental results shown in Figs 5-7.

6. CONCLUSIONS (1) In order to compare the numerical analysis results with experimental results, glass was used in the sliding test. Periodic surface cracks on the surface behind the slider were observed, and it was noted that the spacing between the surface cracks increased with the edge radius when the depth of indentation was constant. The normal and tangential forces increased as the edge radius increased. (2) Photographs of the cross-section of glass showed that the propagation path of the surface crack depended on the radius of the slider and the load. Some cracks propa-

r = 0.23 ram

r = 0.23 mm

~

I

O=-lOdeg

r = 0.23 nun

~

r = 0.23ram

0 = 0 deg

r ffi0.23mm ~

+ lOdeg

O=-5deg

t

= + 5 deg

[r ffi0.76 m m ] I

Of- 5deg

Fig. 16. The possible directionsof crack extension(r = edge radius, 0 = first calculatedpossibleangle).

Influence of Cutting Edge Radius on the Machined Surface of Brittle Materials

983

gated in the sliding direction (Type A) and others in a direction opposite to sliding (Type B). (3) The periodic occurrence of surface cracks could be explained by the gradual increase in the maximum principal stress on the surface behind the slider as the slider moved away from a previous crack. Approximate crack spacings were calculated and these compared well with the experimental results. (4) A method was used to calculate possible directions of crack propagation. Ideally an analysis could be done to yield a fairly precise answer of crack shape by assuming many very short cracks, then considering only those solutions within a range of angles, e.g. 0.5 ° or less. This would require a very lengthy procedure. The present method was considered valid because of the great uncertainty arising from the range of values of Kc of glass, and indeed the range of uncertainty of the calculated results explained the range of angle of crack extension seen in the experiments or in practice with a blunt slider sliding along a brittle material. REFERENCES [1] J. O. Smith and C. K. Liu, Stresses due to tangential and normal loads on an elastic solid with application to some contact stress problems, J. Appl. Mech. 20, 157 (1953). [2] L. M. Keer, M. D. Bryant and G. K. Haritos, Subsurface and surface cracking due to Hertzian contact, Z Lubri. Technol. 104, 347 (1982). [3] S. Jahanmir, J. W. Dally and Y. M. Chen, Fracture mechanics analysis of near-surface cracks, Proc. of the Japan Int. Tribology Conf., Nagoya, p. 581 (1990). [4] Y. Kimura and M. Shima, Stress intensity at surface cracks in sliding wear--parametric study, Proc. of the Japan Int. Tribology Conf., Nagoya, p. 569 (1990). [5] M. C. Dubourg and B. Villechaise, Analysis of multiple fatigue cracks--part I: theory, J. Tribol. 114, 455 (1992). [6] M.C. Dubourg, M. Godet and B. Villechalse, Analysis of multiple fatigue cracks--part If: results,J. Tribol. 114, 462 (1992). [7] C. F. Shih, H. G. Delorenzi and M. D. German, Crack extension modeling with singular quadratic isoparametric elements, Int. J. Fracture 12, 647 0976). [8] V. Mufti, S. Valliappan and I. K. Lee, Stress intensity factor using quarter point element, J. Engng Mech. III, 203 (1985). [9] F. Erdogan and G. C. Sih, On the crack extension in plates under plane loading and transverse shear, J. Basic Engng 85, 519 (1963). [10] B. Cotterell and J. R. Rice, Slightly curved or kinked cracks, Int. J. Fracture 16, 155 0980). [11] H. D. Hibbit, B. L. Karlsson and E. P. Sorensen, A B A Q U S : a General Purpose Finite Element Code, Version 4.6, Hibbit, Karlsson and Sorensen, Inc., Providence, RI (1987). [12] G. W. McLellan and E. B. Shand, Glass Engineering Handbook, 3rd edn, chap. 6. McGraw-Hill, New York (1984).