THE U N I V E R S I T Y OF M I C H IG A N College of Engineering Department of Nuclear Engineering Laboratory for Fluid Flow and Heat Transport Phenomena Report No. 07738-5-T NUMER ICAL AND PHOTOGRAPHIC STUD IES OF ASYMMETRIC BUBBLE COLLAPSE T. M. Mitchell C. L. Kling R. Cheesewright F. G. Hammitt Financial Support Provided by: National Science Foundation Grant No. GK-730 July 1967

ABSTRACT Recent investigations of the fundamental aspect of cavitation damage have indicated the importance of asymmetric bubble collapse. This investigation is an effort to theoretically and experimentally justify the role of asymmetric bubble collapse in cavitation damage. The theoretical investigation employs a combined Eulerian and Lagrangian computer code capable of following the collapse of a bubble much further than previous perturbation techniques. The experimental investigation is concerned with the development of a specialized spark camera with a maximum framing rate of 106 frames per second and a spatial resolution higher than any other camera presently available at this laboratory. ii

TABLE OF CONTENTS Page Abstract......................... ii Table of Contents...I............. *. iii List of Figures.....................iv Introduction............ 1 Theoretical Considerations...2 Experimental Considerations. 14 Bibliography.....................28 iii

LIST OF FIGURES Page l(a) Grid and Initial Marker Particle Location for MAC 5 Solution of the Spherically Symmetric (One-Dimensional) Problem; (b) Grid and Initial Marker Particle Locations for MAC Solution of the Axially Symmetric (Two-Dimensional) Problem. 2 Comparison of Radius as a Function of Time by Two 11 Methods: MAC and Simpson's Rule Numerical Integration. 3 Effect of Variable Time Interval on Accuracy of the 13 MAC Method. 4 Wedge and Spoiler in Plexiglas Section of Two- 15 Dimensional Venturi. 5 Dynafax Photographs Showing Separation Bubble. 16 6 Dynafax Photographs Showing Elimination of Separation 17 Bubble. 7 Schematic Diagram of Spark Camera. 22 8 Spark Camera Including Delay Network and High Voltage 23 Power Supply. 9 Spark Camera Photographs Showing Water Jet Impacting 25 on a Plexiglas Specimen. iv

THE ASYMMETRIC COLLAPSE OF CAVITATION BUBBLES Introduction: Over the last few years, interest in the fundamental aspect of cavitation which is concerned with the collapse of a single isolated bubble has moved away from the classic concept of symmetrical collapse exemplified by the work of Rayleigh (1), Hickling and Plesset (2), Ivany and Hammitt (3), etc. Following an initial suggestion by Kornfeld and Suvorov (4) and important work by Rattray (5), Naude and Ellis (6), and Shutler and Mesler (7), attention has been increasingly directed towards the consideration of asymmetric collapse phenomena. This interest in the asymmetric collapse mode is justified by general considerations of the collapse phenomena (e.g., Benjamin and Ellis (8)), as well as the experimental evidence of Ivany, Hammitt, and Mitchell (9), Shutler and Mesler (7), etc., and analytical considerations of the early stages of the collapse phenomena by Rattray (5), Naude' and Ellis (6), and Yeh (10), which have all shown a tendency for the cavitation bubble to become asymmetric during the collapse. These latter theoretical considerations have all made use of perturbation techniques based on the classical symmetric solution, and have therefore been limited to relatively small departures from spherical symmetry. One aspect of the present investigation is aimed at extending the theoretical considerations of the asymmetric collapse phenomena to much

later stages of collapse than is possible with perturbation techniques. With this object in view, a study has been made of some of the powerful techniques which have been developed over the last few years for solving complex fluid mechanics problems using a high speed digital computer. The first part of this report describes the results of this study which led to the choice of the MAC (Marker and Cell) method as being the most suitable for the problem in hand. This method was first applied to the classic problem of symmetric collapse in order to gain a feeling for the problems involved in its application to problems of bubble collapse. The results of this preliminary investigation and their implications with regard to the study of the asymmetric collapse, which is presently under way, are reported. Another aspect of this present investigation which is concerned with obtaining more accurate and more detailed experimental data on the asymmetric collapse phenomena is contained in the second part of this report. This part is mainly concerned with the development of a special camera having sufficient spatial and temporal resolution to record the details of the bubble collapse on film. Theoretical Considerations: In selecting a suitable technique for the calculation of the asymmetric mode of bubble collapse, consideration was given to the long term objectives. These long term aims include the desirability of being able to consider the effect of such factors as bubble gas content, fluid viscosity,

arbitrary pressure gradients, and the presence of a solid wall, whereas the short term objective is merely to study the asymmetric collapse of a vapour bubble in an infinite volume of perfect fluid, subject to a step change in fluid pressure. The selection of the most suitable technique is at first a matter of elimination. The non-linearity of the partial differential equations together with the free surface boundary condition precludes a direct analytic solution. Numerical solution of the difference equations with an Eulerian (stationary) grid cannot handle the free surface condition. A purely Lagrangian numerical solution, on the other hand, will not adequately describe a fluid flow in the presence of large distortions such as may be anticipated in the asymmetric collapse problem. A combination of the Eulerian and Lagrangian approaches could offer some of the advantages of both methods while overcoming some of the disadvantages. One example of the combined approach is the Particle-inCell (PIC) method (11), developed by Harlow at Los Alamos. This method, which was initially developed for use in compressible flow problems, was discarded because it requires that no part of the flow be far subsonic, a condition which is untrue for most of the bubble collapse problem. The eventual choice for the appropriate method came to be the Markerand-Cell (MAC) method developed by Harlow and Welch at Los Alamos (12). This method appears to have worked very well for solution of a variety of problems involving incompressible liquids with free surfaces. It requires a high-speed computer with extensive storage capacity.

- 4 - With the eventual solution of the asymmetric problem in mind, we have first applied the MAC method to the problem of the spherical collapse of a cavitation bubble (Rayleigh's problem), because an analytic solution already exists for this problem and we can thus determine the characteristics of this technique as applied to bubble collapse problems and the effect of such parameters as grid sizes and time interval. The MAC method starts with an Eulerian grid as shown in Figure 1 (a). Field variables such as pressure and velocity are directly associated with the cells of the grid. In addition to the grid, a series of marker particles are assigned to the liquid; these particles do not participate in the dynamics of the problem and are essentially "massless", but they are necessary to take care of the free surface movement. Although the Los Alamos work suggests that marker particles be located throughout the liquid initially, experience with the symmetric collapse problem has shown that a considerable saving in computer time and storage can be obtained by starting with particles located only along the free surface. This arrangement is possible because of the condition (see Lamb (13) ) that fluid particles initially on the free boundary always remain on the free boundary. In the spherically symmetric case, the above considerations lead to the use of one single particle (marked A in Figure 1 (a) ) to define the free boundary. In the asymmetric case it is intended to start with particles located along the free boundary as in Figure 1 (b).

r LIQUID-VOID INTERFACE 2488 Fig. l(a) —Grid and Initial Fig. l(b) —Grid and Initial Marker Particle Location for MAC Marker Particle Locations for MAC Solution of the Spherically Sym- Solution of the Axially Symmetric metric (One-Dimensional) Problem. (Two-DimensiGnal) Problem.

To explain the technique in more detail, it is convenient to take as the starting point of the problem, the hydrodynamic equations for continuity and conservation of momentum, in one-dimensional spherical coordinates: DD(JZg) - 0 (1) -l + (. V /t = - CP (2) where u = the radial velocity in the liquid, P = p is the liquid pressure, and p is the liquid density (constant for the incompressible problem). The boundary conditions for the problem are: P_ pO, As; — oo - (3) P=O fob Ro R t (4) where R is the radius of the bubble wall; R (o) - Ro(5) R o) O (65 'elf() Q (6) Equation (2) can be rewritten as -a 2- - D<-~ D ~P (7)

-7 -Harlow and Welch point out that this transformation is essential for rigorous momentum conversation after the equations are written in difference form. The cells of the grid system are then labeled with an index i; further, the P-quantity (=p/p) is defined at the center of each cell as Pj and the velocities are defined on the cell-boundaries as Ui-1/2 and Ui+l/2. The definition of these quantities at any other locations increases the complexity of the calculation and brings in values from "far-distant" points, thus contributing to increased inaccuracy. It is then possible to write equation (1) in difference form as: i — -' — (8) ' An h= All 1 ' -- ",,. Equation (7) similarly becomes (r,* _____ (9) The superscript n+l indicates that the quantity is evaluated at t = (n+l) t while the absence of a superscript indicates that the evaluation is at t = n t. The quantity Ui is defined as 1/2(Ui+/2+il-1/2). If, in addition, equation (8) is rewritten for Din+l and (9) is rewritten for uin 12 and the resulting four equations are combined, one obtains

8-,2 cu142 (Di D1).L. ( ( AL- g Equations (9) and (10) provide the basis of the computational method. Knowing the pressure distribution at any time t, equation (9) is solved to give the velocity distribution at t+S t; equation (10) can then be solved, assumed Din+l= O,to give the pressure distribution at time t+% t. Two questions arise in connection with this scheme of computation. The first is concerned with the initial conditions to be imposed at t = 0. Strictly this should be P = POO for r > Ro at t = 0. Because of the computational scheme, however, the pressure discontinuity initially at the bubble wall would require a finite time to propagate out through the liquid and thus would violate the assumption of an incompressible fluid. In the symmetric case this difficulty has been overcome by the use initially of the pressure distribution as given by the Rayleigh solution. In the asymmetric case it is felt that it will be necessary to begin with either the Rayleigh solution or perhaps the results of a perturbation solution (Rattray, Yeh, etc.) for small degrees of asymmetry. The second question concerns the boundary condition at the outer extremity of the grid system. The true boundary condition in the Rayleigh

problem is that P is specified at infinity. Since we are forced to use a grid system of finite extent, it is not possible to employ this condition. To overcome this difficulty, the pressure predicted by the Rayleigh solution has been used as the boundary condition at the outer edge of the grid system. It is proposed to use the same procedure in the asymmetric collapse problem. After solving equations (9) and (10) for the velocity and pressure distributions for time t + 3 t, the marker particle is moved. The movement is accomplished by looking at the present location of the particle, interpolating between the velocities at the boundaries of the cell in which the particle is located to find the particle velocity, and then moving the particle a distance equal to this velocity times S t. The procedure then is to scan the grid to see if the particle has moved into another cell and, if so, to incorporate the new cell into the succeeding Eulerian calculations. After printing any necessary outputdata, the cycle of calculations is repeated. The primary goal of the study of the spherical bubble problem by the MAC method has been the investigation of the effects of different grid sizes and time intervals. 25 and 40 cell grids were used for two different sets of calculations. Use of the 40 cell grid reduced the error in the initial stages of the calculation by about 1/3 from that with the 25 cell grid, at a cost of a 20% increase in computing time.

- 10 - At later stages in the collapse, however, the error improvement disappeared, although this was masked by the constant time step error explained below. Two problems were apparent when constant time steps of 5 and 10 us were utilized. In the late stages of collapse, the bubble wall velocity diverged rapidly from the expected Rayleigh value. This difficulty arose apparently because of stability considerations violated when the particle moved across a large part of a cell in one time interval. This difficulty was overcome by varying the time interval so as to impose a maximum on the particle movement in one time cycle of 10/o of a single cell width. The other difficulty concerned the error in the early stages of the collapse. This error is dependent on both grid size and time increment, but a method of reducing this error without involving a large increase in computing time, has not yet been achieved. The error does, however, appear to die out with time and so possibly may be ignored in the asymmetric case. A check on the accuracy of the computed results can be made by comparing the results with corresponding values calculated from Rayleigh's analytical solution. Figure 2 shows a plot of the bubble wall radius vs. time. The solid line indicates the Rayleigh solution; the data points are those obtained from the Marker-and-Cell program. The MAC method

100 200 TI00 sS 000 400 ~iO0 600 700 Boo 900 1000 INITIAL CONDITIONS: POi atntm, R0 a cm IO~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1.0 RO A MAC COMPUTATION NUM. INTEGRATION 24,9 Fig. 2.*Comparison of of Time by Two Methods: MsAC ard Simpson'c Rule?fumerical Integration,

- 12 - data virtually coincides with the Rayleigh results except for a slight divergence at the very end of collapse. Figure 3 shows another comparison of the two methods expressed in terms of the percentage error as a function of time (the percentage error is defined here as the difference between the Rayleigh velocity for the present bubble radius and the MAC velocity divided by the Rayleigh velocity). The curves of Figure 3 show the improvement arising from using a variable time step rather than one of constant value and also the problem of the initial error. It is hoped that the following conclusions drawn from the study of the spherically symmetric case can be applied to the axially-symmetric two-dimensional case: 1. Marker particles need only be located initially along the liquidcavity interface. 2. A variable time step is necessary for accuracy. 3. For greatest accuracy, some initally correct pressure distribution should be imposed upon the liquid. 4. Approximately 25 cells radially are sufficient to give an accuracy of the order of 5%. It is hoped that this can be improved further by correcting the initial error problem. 5. The pressure in the one cell necessarily external to the grid system must be controlled by an appropriate expression, possibly the correct Rayleigh pressure as was used in the spherically symmetric case.

TIME/(/ps) 0 100 200 300 400 500 600 700 800 900 1000 P E G 8, R -12 R I VARIABLE TIME INTERVAL PERCENTAGE ERRORRA G COMPUTER -20 (STARTING AT 10 S) VRAYLEIGH A Zt- 10 Us INITIAL CONDITIONS: Po I atm, Ro- I cm 2490 Fig. 3. —Effect of Variable Time Interval on Accuracy of the MAC Method.

- 14 - Experimental Considerations: The objective of this part of the investigation is to determine the detailed mechanism of the asymmetric mode of cavitation bubble collapse. The asymmetry is introduced into the experimental system by the presence of a solid wall. The previous experimental study undertaken with a Fastax camera has been discussed in reference (14). Since this initial investigation, several modifications have been made. The major modification to the two-dimensional venturi used in this investigation has been the attachment of a spoiler to the adjustable bolt shown in Figure 4. By turning the bolt, the area of the spoiler perpendicular to the flow can be changed; thereby, partially controlling the division of flow at the tip of the wedge. Experimentation, while taking pictures with both the Fastax and Dynafax cameras, showed that an unequal division of flow at the tip of the wedge caused the separation bubble shown in Figure 5. By placing the spoiler on the same side of the wedge as the separation bubble, the flow division at the tip could be equalized and the bubble nearly eliminated as shown in Figure 6. Additional modifications of the venturi have been undertaken which included gluing the four sides together and installing a plexiglas window to facilitate adjustment and removal of the wedge. These modifications serve the two-fold purpose of eliminating the leaks associated with the multitude of seals previously required and of reducing the time required for wedge adjustment by a factor of ten.

- 15 -Wedge Pivot' F1nw/ PivotX TOP VIEW Adjustable bolt Top plate Throat Diffuser ISpoiler Flo w... — Field of view Bottom plate 249 SIDE VIEW Fig. 4. —Wedge and Spoiler in Plexiglas Section of Two-Dimensional Venturi.

Mo 0'I qS2uaq oeoS 'omvei ad spuoOasoJDTW 08 'axqqnfg uoTjvjdoS SUT&~RS sqdvaSolOqd xvJVuaaG —' *S a........::::....................^..........a............ a" '""BB X~~~~~~~~~~~~~~~~~~~~i"'0::i''00 00:'ll; i i i ii i: 0; * |~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~:::;::::~1:iii~i.:!i~~ 'B::::R:::.::::':.................0;,:,,v~~g;:i.-i: ti::... -.: -..... -.R'..',R -fi..ii,~~!?,:: ~.........: i.:.:.. 9..........I.........i,..,,,~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~.-:.............

17 oil~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~it ii ~. ~........... -iii "~ ~~..........?i ~~ ~ EA 2 T7%::- - --.:.............- - ~. Z'......... ~ Fig. 6. —Dynafax Photographs Showing seconds erra, S cale- L e n g t h -.0 cm.

18 Fig. 6. — (Cont.)

19 Fig. 6. —(Cont.) 1i

- 20 - By the time tests were completed showing that the flow problem had been resolved, the conclusion was reached that the Dynafax camera,although offering a three fold increase in framing rate (25,000 frame per second maximum) over the Fastax camera, still had neither the framing rate nor the resolution required. Results obtained both with the Fastax and the Dynafax cameras show that a bubble lifetime in the two-dimensional venturi is on the order of 3 ms, while the final collapse phase occurs in both cases in the time between frames, or less then 40 ps. Therefore, a camera capable of photographing a bubble during this final collapse is required. It may also be seen from the small size of the bubble during the important final stages of collapse, that a camera with improved spatial resolution, as compared with the Fastax and Dynafax cameras, is required. Only a few types of cameras are capable of producing several photographs within a 40 ps period, i.e., framing rates of 105 frames per second or greater. Ultra-high speed rotating mirror cameras are capable of framing rates in the millions of frames per second with exposure times on the order of the reciprocal of the framing rate. 35 mm film is generally used, limiting negative image size and therefore resolution. Unless a continuous access type of rotating mirror camera is employed, syncronization of the camera with the random occurrence of the cavitation bubbles of this investigation would be possible. Spark cameras are capable of framing rates in the millions of frames per second if a reliable delay

- 21 - network can be constructed. Exposure times are equal to the duration of the spark light source and are generally a fraction of a microsecond. The only other common camera capable of such high framing rates is the image converter camera. Both exposure time and framing intervals may be measured in nanoseconds giving extremely high framing rates on the order of 108 frames per second. However only a very small number of frames, about 5, can be obtained and the spatial resolution is an order of magnitude less that than obtainable with rotating mirror or spark cameras. A very high speed rotating mirror camera, having continuous access, has been ordered for use in this investigation, but it will not be available for at least six months. In order that the work may proceed during this period, a specialized spark camera is being constructed. A schematic of the camera is shown in Figure 7 and a picture of the physical setup is shown in Figure 8. The camera uses a series of spark light sources mounted in a circle. These spark light sources are controlled by a twelve stage delay network. The network has an initial delay that is continuously variable from 2 to 1000 ps with each successive delay variable from 1 to 100 ys in seven steps. The effective framing rate is from 104 to 106 frames per second. The pulse width (or exposure time) of the discharging spark light source is approximately 0.2 us.

LIGHT SOURCES CAMERA LENSES 000 00 0 0o0 0~~~~~~~~~~~0 ~ ~ ~ ~ SECTION A-A SECTION B-B SPARK LIGHT CENTRAL WORKING CAMERA FILM SOURCE A LENS rSECTION B PLANE "-OPTICAL BENCH A B 2494 Fig. 7. —Schematic Diagram of Spark Camera.

Fig. 8.-Spark Camera Including Delay Net. work and High Voltage Power Supply.

- 24 - The light from each light source is focused by the large central lens onto one of the camera lenses. There is a camera lens for each light source and the optical system is arranged so that the camera lenses and light sources are mounted at the same radius. In turn, the camera lenses focus the working section on the stationary film plane. The resulting negative has a circle of images on a radius slightly greater than the camera lenses. The camera is capable of holding 8 inch x 10 inch sheet film in order to obtain negative images as large as possible. A representative sequence of pictures from the initial use of this camera shows a liquid jet emerging from an orifice and impacting on a plexiglas specimen (Figure 9). The high resolution and short exposure time of this camera are plainly evident. The major limitation of this camera, for our work, is the low number of frames per run, 12. This limitation is imposed by the physical size of the spark sources used and the need to make the negative images as large as possible. It may be possible to partially overcome this disadvantage by providing some means of detecting a bubble in its latter stages of collapse and then starting the camera with this detector. Other limitations of this camera not pertinent to our work are that its use is restricted to subjects which are non-luminous and can be back lighted.

- 25 Fig. 9. —Spark Camera Photographs Showing Water Jet Impacting on a Plexiglas Specimen. Time in Microseconds Measured from Frame Number 1, Scale Length 0.5 cm.

- 26 - Fig. 9. — (Cont.)

- 27 - The high speed rotating mirror camera which will eventually be available for this work will have a maximum framing rate of 2 x 106 frames per second and a total of 82 frames per sequence. The image size will be somewhat smaller than that obtained with the spark camera, causing a corresponding reduction in resolution, but this is not expected to be a serious limitation.

- 28 - BIBLIOGRAPHY 1. Lord Rayleigh, "On the Pressure Developed in a Liquid During the Collapse of a Spherical Cavity," Phil. Maq., 34, (1917) 94-98. 2. Hickling, R. and Plesset, M., "Collapse and Rebound of a Spherical Bubble in Water," The Physics of Fluids, 7, No. 1, (1964) 7-14. 3. Ivany, R. D. and Hammitt, F. G., "Cavitation Bubble Collapse in Viscous, Compressible Liquids - Numerical Analysis," J. Basic Enq., Trans. ASME, D, 87, (1965) 977-985. 4. Kornfeld, M. and Suvorov, L., "On the Destructive Action of Cavitation," J. Applied Physics, 15, No. 6, (1944) 495-506. 5. Rattray, M., "Perturbation Effects in Cavitation Bubble Dynamics," Cal. Inst. of Tech., Ph.D. Thesis, 1952. 6. Naude, C. F. and A. T. Ellis, "On the Mechanism of Cavitation Damage by Nonhemispherical Cavities Collapsing in Contact with a Solid Boundary," Trans. ASME, J. Basic Enq., (1961) 648-656. 7. Shutler, N. D. and Mesler, R. B., A Photoqgraphic Study of the Dynamics and Damage Capabilities of Bubbles Collapsinq Near Solid Boundaries, Trans. ASME, J. Basic Eng., D, 87, (1965) 511-517. 8. Benjamin, T. B. and Ellis, A. T., "The Collapse of Cavitation Bubbles and the Pressures Thereby Produced Against Solid Boundaries," Trans. Royal Society, A, 260, (1966) 221-240. 9. Ivany, R. D., Hammitt, F. G., and Mitchell, T. M., "Cavitation Bubble Collapse Observations in a Venturi," J. Basic Eng., Trans ASME, D, 88, (1966) 649-657. 10. Yeh, H., "The Dynamics of Gas Bubbles Moving in Liquids with Pressure Gradient," Ph.D. Thesis, University of Michigan, 1967. 11. Harlow, F. H., "The Particle-in-Cell Computing Method for Fluid Dynamics," Methods in Computational Physics, 3, (1964) 319-343. 12. Welch, J. E., Harlow, T. H., et al, "The MAC Method: A Computing Technique for Solving Viscous Incompressible Transient Fluid - Flow Problems Involving Free Surfaces," Los Alamos Scientific Laboratory Report LA-3425 (March, 1966).

- 29 - Lamb, Sir Horace, Hydrodynamics, Dover Publ., 6th Ed., New York, (1945) 7. Mitchell, T. M., Kling, C. L., and Hammitt, F. G., "Bubble Collapse Behavior Around a Sharp Wedge in a Diffuser," Department of Nuclear Engineering Laboratory for Fluid Flow and Heat Transport Phenomena, University of Michigan, Report No. 07738-T-1, July, 1966.