Ultrasonic liquid metal processing: The essential role of cavitation bubbles in controlling acoustic streaming

The acoustic streaming behaviour below an ultrasonic sonotrode in water was predicted by numerical simulation and validated by experimental studies. The flow was calculated by solving the transient Reynolds-Averaged Navier-Stokes equations with a source term representing ultrasonic excitation implemented from the predictions of a nonlinear acoustic model. Comparisons with the measured flow field from Particle Image Velocimetry (PIV) water experiments revealed good agreement in both velocity magnitude and direction at two power settings, supporting the validity of the model for acoustic streaming in the presence of cavitating bubbles. Turbulent features measured by PIV were also recovered by the model. The model was then applied to the technologically important area of ultrasonic treatment of liquid aluminium, to achieve the prediction of acoustic streaming for the very first time that accounts for nonlinear pressure propagation in the presence of acoustic cavitation in the melt. Simulations show a strong dependence of the acoustic streaming flow direction on the cavitating bubble volume fraction, reflecting PIV observations. This has implications for the technological use of ultrasound in liquid metal processing.


Introduction
The complexity of non-linear acoustics and cavitation phenomena, the opaqueness, high temperatures and chemical reactivity of metallic melts hinder the study of ultrasonic melt processing. To improve the understanding of ultrasonic process effects on liquid metals, an extensive research program was undertaken over the past five years [1] to implement a numerical model that can realistically predict acoustic pressures and acoustic streaming in the melt in the presence of cavitation, and to validate the model experimentally.
Processing melts with ultrasound improves physical and mechanical properties of the treated metallic materials [2][3][4]. These beneficial effects of ultrasonic melt processing are attributed to acoustic cavitation and acoustic streaming before, and during, the solidification of the liquid metal [5]. While this technology has been successfully applied in the laboratory and at the pilot plant scale, further up-scaling requires multiple ultrasonic sources, which undermines the economic attraction, and imposes technological restrictions on widespread industrial adoption. To up-scale ultrasonic melt processing to the industrial scale effectively, there is a need to bridge the gap between understanding of ultrasonic fundamentals and the interaction between cavitation and larger-scale melt flow, with the aim of operating the ultrasound devices in a manner that minimises the number of ultrasonic sources while simultaneously maximising the volume of effectively treated melt [6].
Quantifying recirculation patterns and mass exchange between the cavitation zone and the rest of the melt bulk is crucial for understanding how to optimize ultrasonic melt processing because this effectively determines the treated volume. Recirculation of the liquid melt, over length scales that are matched to the dimensions of the melt volume, will also help to reduce temperature gradients in the melt, thereby promoting a preferred equiaxed grain structure [7]. In the direct-chill (DC) casting process that is used widely in the aluminium industry, acoustic streaming promotes forced convection which has been shown to decrease macrosegregation and to promote solid fragmentation and thus self-grain refinement [4,8]. Other beneficial effects include deagglomeration and wetting of inclusions [9], and their dispersion that increases the number of substrates available for heterogeneous nucleation of the solid, thereby promoting grain refinement.
Numerical models describing the dynamics of acoustic cavitation can synthesize existing empirical knowledge and provide a framework for understanding the complex mechanisms involved. Accurate prediction of acoustic streaming is however plagued by difficulties, especially the challenges in modelling nonlinear acoustic pressure propagation [10]. Therefore, it is not surprising that a model of acoustic streaming in the presence of acoustic cavitation has only recently appeared in the literature [11]. Earlier attempts of modelling acoustic cavitation and flows in liquid metals [12] were limited due to the use of a homogeneous cavitation model [13] that is only applicable to liquids containing vapour bubbles [11]. More accurate efforts in quantifying acoustic pressures in liquid aluminium [14,15] employ non-linear equations as suggested by van Wijngaarden [16] following previous successful implementation in water [17,18]. However, these models are computationally expensive, as the dynamics of bubbles in each computational cell have to be resolved [10,15]. Incorporating bubble motion in such models [19] results in additional computational complexity.
From the experimental perspective, there are relatively few quantitative studies on acoustic cavitation and streaming during ultrasonic processing. Acoustic streaming is known to form a conical or jet flow pattern, depending on the input acoustic power, with recirculating flows influencing the entrainment of bulk liquid and small solid inclusions back into the cavitation zone [20][21][22]. Using Particle Image Velocimetry (PIV) experiments in low temperature liquids (mostly water) and high-speed synchrotron X-ray radiography experiments, many groups have attempted to visualise and understand what is happening inside the cavitation zone of a liquid melt. Feng et al. [23] described the mechanisms of microstructural refinement based on real time acoustic flow observations of a sonicated Al-35Cu alloy melt. Mirihanage et al. [24] reported the flow velocity of the acoustic streaming by observing the bubble streamlines and local turbulent fluctuations of an Al-Cu metal matrix nano-composite melt. Bing et al. [25] investigated the effect of acoustic flow on the liquid-solid interface of a Bi-8%Zn alloy using X-ray radiography and image analysis: a toroidal flow pattern was observed, coaxial with the sonotrode, appearing in 2D as a clockwise vortex on the left and a counter-clockwise vortex on the right. The recorded velocities were in the range of 0.5-0.6 ms −1 . Overall, there remains a lack of a generic, validated approach that can describe behaviour across a range of conditions, materials and geometries that can guide industrial practise, with the ultimate goal of ultrasonic treatment implementation in casting processes.
Recent progress in the theory on non-linear sound propagation [26][27][28][29] resulted in an easier-to-solve nonlinear Helmholtz equation to quantify the acoustic pressure field. Louisnard extended his nonlinear model to account for acoustic streaming in the presence of cavitation with two-dimensional results comparing well with experiment [11]. In this paper, this acoustic streaming model is adapted to a finite volume computational fluid dynamics solver compiled in OpenFOAM [30] to predict the flow below the sonotrode in both water and liquid aluminium. The model is compared with PIV measurements in water sonicated at two different acoustic powers. Water is used for comparison, because it is deemed a suitable physical analogue to aluminium for studying ultrasonic melt processing [31]. The model is then applied to aluminium to predict the acoustic streaming pattern at different   transducer powers and the likely implications for liquid metal processing are discussed.

Acoustic field model
The Caflisch equations describe acoustic propagation in a bubbly liquid represented as a continuous medium [17]. Following conservation of mass and momentum, we obtain where is the (pure) liquid density, c is the speed of sound, p a is the acoustic pressure, u is the velocity, and is the bubble volume fraction given by is the bubble volume, and N is the bubble density (number of bubbles per unit volume). To complete the Caflisch model, the Rayleigh-Plesset equation [32,33] describes the bubble dynamics where the bubble pressure is given by with p g,0 being the gas pressure at the equilibrium radius R 0 . Assuming an adiabatic gas, the polytropic exponent 1.4 = , p v is the vapour pressure, is the interfacial tension between the gas and the liquid, µ is the dynamic viscosity of the liquid, p 0 is the pressure at infinity (set to atmospheric pressure) and p A is the pressure amplitude due to the excitation source whose angular frequency is During rapid pulsation, a cavitating bubble undergoes rapid changes in pressure p b and temperature. Heat transfer between the gas and liquid during these pulsations results in energy dissipation. A nonlinear model can be obtained by taking into account gas continuity and energy conservation without explicitly solving for them [26,28,34]. We follow Louisnard's approach [26] to obtain the following dissipation functions: where T f 1 = is one period, th is the period-averaged heat loss from the bubble and v is the period-averaged dissipation by viscous friction.
Denoting the harmonic part of p a as Pe , the complex amplitude P approximately fulfils the nonlinear Helmholtz equation [26,27] P K P 0, where the real and imaginary parts of K 2 are given by The resonant frequency 0 is given by [34,35] and D is the gas diffusivity. Table 3 Boundary conditions.
No slip Normal gradient = 0 Pa m −1 kqRWallFunction (zero gradient boundary condition), omega wall function [40], nutkWallFunction (turbulent kinematic viscosity condition when using wall functions, calculated from turbulent kinetic energy) [30] Fixed The bubble density is assumed to follow the step function [26] where the Blake threshold is P p 1 B S S 0 Note that the bubble fraction is not explicitly computed in this model. Instead, the acoustic pressure is estimated by solving equation (7), from which the acoustic streaming force can be inferred. Also, the model is concerned with the ultrasonic processing of light alloy melts or analogues, and not in simulating the different process of bubble entrainment (e.g. due to degassing). Therefore, the current model only considers the effect of internally generated cavitation bubbles.

Fluid flow model
Fluid flow is modelled using the incompressible Reynolds-Averaged Navier-Stokes (RANS) equations closed with the kshear stress transport model [36]: where p is the flow pressure, µ/ = is the kinematic viscosity, and f represents the force driving acoustic streaming.
where v P = is the acoustic streaming velocity [11]. The overbar indicates that the values are obtained from averaging over a period of the acoustic bubble (specifically in the calculation of th and v . .
The turbulent viscosity is given by

Geometry
The physical arrangement corresponds to the experimental configuration of particle image velocimetry (PIV) experiment in a water tank, illustrated in Fig. 1. A 20 kHz transducer (Hielscher UP1000) sonicated water in a rectangular glass tank of base area 300 mm × 200 mm with water height 150 mm. The cylindrical sonotrode with a diameter 40 mm was immersed 20 mm below the water surface. Two conditions were considered: when the transducer was operating at 50% power corresponding to a null-peak amplitude of 4.25 µm, and 100% power corresponding to an amplitude of 8.5 µm.
A TSI (USA) PIV system measured velocities in the plane of a narrow two-dimensional 10 mm × 10 mm laser-illuminated light sheet, or "window", using multiple, serial double short Nd-YAG laser flashes with "frozen" images of the flow filled with glass hollow sphere seeded particles (∼10 μm in diameter) recorded perpendicular to the light sheet at a rate of 5 pair of images per second. Image pairs were then analysed using TSI software to retrieve 2D velocity vector maps which were time averaged over typically 100 pairs. Further details of the PIV arrangement can be found elsewhere [37].

Numerical implementation
The numerical simulations are performed in a modified version of the solver pisoFoam of the open source library OpenFOAM version 4.x. Table 1 lists the material properties that are used to solve the Rayleigh-Plesset equation (4) and the coupled acoustic-flow equations (7) and (12)(13)(15)(16). Table 2 lists the discretization schemes and solver control parameters for the simulations. A summary of the solver algorithm is detailed below: 1. The Rayleigh-Plesset equation is solved for a bubble of equilibrium radius R 0 for a range of pressures which include the pressure below the sonotrode at the operating power. This pressure is estimated from measurements in similar experiments using a calibrated cavitometer [38]. 2. The values of th and v at the operating pressure are obtained from the integrals (5-6) over a period of oscillation.
The following steps are looped over the run time. Even if steady state is expected at very low bubble volume fractions, solving transient equations adds more stability due to the transient term acting as an inertial relaxation term in the flow equations.
3. For each time step, the harmonic component of the acoustic pressure equation (7) is solved as follows: a. Equation (7) is split into two equations for the real and imaginary b. Equations (18) and (19) are solved sequentially using the finite volume method. Note that convergence is achieved only with a suitable use of preconditioners for both equations, i.e. Simplified Diagonal-based Incomplete Cholesky preconditioner (DIC) with DIC smoothing followed by Gauss-Seidel (DICGaussSeidel).
c. The computed values of acoustic pressure are used in the source term (14) 4. If activated, the momentum equations (13) are solved in the momentum predictor step. 5. Pressure and velocity are solved for using the Pressure-Implicit Splitting of Operators (PISO) algorithm [39]. 6. The turbulence model is solved during the last stage of the time loop. Fig. 2 and Table 3 show the boundary conditions of the simulation. The surfaces of the tank are rigid walls at which the no-slip condition is imposed. The free surface is modelled as a surface of constant pressure. The pressure at the sonotrode is prescribed by specifying the expected pressure normal gradient as specified in Table 3 [28]. The pressure below the sonotrode was indirectly prescribed by fixing the pressure gradient at the surface as v y n 2 = . The walls of the sonotrode were presumed to be rigid walls.

Application of the acoustic streaming model to water
The numerical model was solved in the two-dimensional plane of the laser sheet shown in Fig. 1. Comparison with experiment was principally qualitative as a full three-dimensional simulation is required for quantitative comparison, which is the focus of further work, together with the use of large-eddy simulations (LES) to model turbulence more accurately [41,42]. Nevertheless, a qualitative comparison is helpful to understand, for example, a surprising upward flow in PIV experiments at low power, as shown in Fig. 3(a). At 50% of the maximum sonotrode power, the net average flow along the axis of the sonotrode was upwards, with recirculation patterns around the axis. At full power, the acoustic streaming forces the liquid at high speed towards the bottom of the tank (which is the streaming pattern usually reported). Fig. 4 shows the recording of the flow field with the PIV at 2 s intervals. These 6 snapshots are among the 100 vector file recordings used to calculate the average flow pattern presented in Fig. 3(a). The flow is always upwards. Note that, while the sonotrode is oscillating up and down, the fluid below the sonotrode surface will also be pushed upwards and downwards with the forcing frequency of the sonotrode. However, the resolution of the PIV equipment cannot capture this fast instantaneous motion and what is recorded is the net average flow. This net average flow is of particular significance in ultrasonic melt processing where ultrasound is continuously fed into the melt for time frames of the order of minutes. Fig. 5 illustrates the predicted flow pattern as a function of bubble volume fraction, with velocities averaged over 10 s and 20 s of simulation time. This relatively large time simulation time is required at low powers to ensure that the solution for the acoustic streaming pattern is compatible with the time-averaged approach also applied to the PIV data. The whole 2D domain is presented in these results. The hatched rectangle at the top of the contours represents the sonotrode. The dotted rectangular window represents the PIV windows in which direct comparison with the experiments is possible. The left contours  represent velocities while the right contours represent the flow pressure p. The flow direction is represented using streamlines overlaid on the pressure contours.

Acoustic streaming pattern at 50% power
The change in bubble volume fraction in the code was implemented by setting different values of bubble number densitiesN 0 in Eq. (11). The results in Fig. 5 show a strong dependence on the predicted acoustic pattern on the bubble fraction in the computational domain. Bubble volume fractions 0.01% < gave good qualitative agreement with the experiment, with a maximum velocity of 0.04 m s −1 below the sonotrode surface, a net upflow along the sonotrode axis, and turbulent recirculation at the sides. The location of the maximum velocity is predicted to be slightly lower than the experiments, but the maximum is of the correct order of magnitude. At larger bubble volume fractions 0.01% > , the numerical simulations became unstable and diverged to unrealistic solutions: these velocities are not presented here. These results demonstrate one shortcoming of the model: the large sensitivity of the model to bubble number density.
The flow direction can be understood by the flow pressure pattern: zones of low pressure are predicted in the middle of the PIV window, i.e. 40 mm below the surface of the sonotrode, and the flow is forced towards these zones from the bottom of the tank. There is also a predicted net low (flow) pressure zone just under the sonotrode (located 40 mm below the sonotrode surface) as indicated by the white contours in Fig. 5 (right) explaining why the net flow below the sonotrode is upwards. However, a high-pressure zone also pushes some flow down between 20 and 35 mm below the sonotrode: this was not observed in the experiments and is a short-coming of this 2D model. This infers a competition between the in the vicinity high and low-pressure zones that regulate the flow direction.
At 0.01% < , the acoustic pressure pattern did not deviate significantly from predictions of the linear Helmholtz equation, as shown in Figs. 6(b) and 7. The pressure pattern was almost identical to that without bubbles i.e. 0% = .
There is a slight difference between Fig. 7(a)-(c): this is better visualized in Fig. 6(b). The three lines do not perfectly overlap in the pressure line plot, but the difference is very low compared with alpha = 0.01%. However, increasing the bubble density to 0.01% = slightly supressed the pressure field, as shown in Fig. 6(b). While the acoustic pressure solver is still robust and forces a 'converged' solution to a pressure profile, the increase in acoustic dissipation results in lower predicted acoustic pressures, but these lead to larger source terms in the momentum equation, thereby making the solver unstable. show good agreement with the PIV measurements in Fig. 3(b). The flow pattern consists of a strong downwards jet with a maximum velocity of 0.24 m s −1 , with the same order of magnitude as the recorded net maximum velocities. The time evolution of one such predicted structure for the case of 0.001% = is shown in Fig. 11. Recirculations were present at the side of the axial jet, and there the conical flow pattern was established just below the sonotrode. Increasing the bubble density further resulted in chaotic, fluctuating flow patterns that indicate numerical divergence.

Acoustic streaming pattern at 100% sonotrode power
Note that the net downflow can be explained by the large high pressure zone sandwiched between the low pressure zones at around 40 mm below the sonotrode. This large relative flow pressure value is not present in the 50% sonotrode case and prevents the net upflow from the bottom of the tank and forces the liquid down, consistent with what is observed in the PIV experiments. As in the water experiment at 50% sonotrode power, large kinks in velocity appear in corresponding large pressure gradient regions.
The acoustic pressure patterns showed a large deviation from linear (Helmholtz) predictions as the bubble fraction increased, as shown in Fig. 10(b) and Fig. 12. For 0.01% = , the acoustic pressures were no longer accurately described by a linear Helmholtz equation and the non-linear terms became dominant. This demonstrates the need for accurate non-linear pressure models to study acoustic cavitation adequately at high forcing amplitudes in the presence of a large fraction of bubbles.

Application of the model to the ultrasonic treatment of molten aluminium
The changes in acoustic streaming pattern due to the presence of bubbles are now investigated by application of the model for the case of the ultrasonic processing of liquid aluminium. The simulations are repeated using the same conditions as the water tank simulations at 50% sonotrode power, but with the liquid replaced by molten aluminium, mimicking the treatment of the alloy just prior to casting and solidification. At low bubble fractions, the results were similar to those of water: a net upward flow along the axis of the sonotrode. On increasing the bubble fraction to larger values, the flow begins to reverse as shown in Fig. 13(h). This effect is also demonstrated in the line plot of the melt y-velocity in Fig. 14(a).
Figs. 14(b) and 15 illustrate the acoustic pressure profile for molten aluminium. At lower bubble fractions, the non-linear contribution to the pressure equation was small and did not affect the acoustic pressure predictions significantly, but pressure was significantly increased as the bubble fraction increased.
The results in aluminium also demonstrate competing interplays in the flow profile inside the melt. On one hand, large pressures are generally considered to be required to de-agglomerate particle clusters inside the bulk [1,43]. However, according to the results here, these large pressures are accompanied by a strong induced downflow that pushes melt away from the sonotrode and the active cavitation zone, and higher velocities that will reduce the residence time of liquid in the energetic zone close to the sonotrode tip. Upward flow, at lower power, below the sonotrode could be attractive in this regard to extend the melt residence time in the cavitation zone. Also, operating at lower power and velocity may help preserve the integrity of the free, top surface of the liquid, which is known to be critical if entrainment of surface oxide [44], which always forms on liquid aluminium, into the bulk of the melt is to be minimised. If ultrasonic processing were applied to the melt sump in direct chill casting of aluminium alloys (that typically contain zinc, copper and/or magnesium) [45], the results of either upward or downward flow on the final microstructure and macrosegregation might be investigated, noting that downflow would "push" hot, relatively dilute liquid down the billet centreline quickly, while upward flow would "suck" relatively cool, more concentrated liquid toward the melt bulk more slowly. The effects on the final grain structure and distribution of alloying elements under these conditions are difficult to predict, but tuning between upward and downward flow may add an extra degree of freedom to industrial practice and the range of final structures that might be contrived.

Conclusions
Acoustic streaming in water has been predicted using Louisnard's non-linear propagation model [26] coupled with transient Reynold-Averaged Navier-Stokes equations, closed with the k-SST model. The coupled set of equations have been solved stably using a finite volume method and the results compare favourably with velocity vector fields measured by particular image velocimetry in water at two ultrasonic powers (at 50% and 100% of available sonotrode power).
A net upward flow on the model centreline was predicted at low sonication power and observed by experiment, which has previously been difficult to predict accurately because the effect of cavitation on the acoustic streaming pattern has been neglected. The model shows a high sensitivity of the acoustic flow pattern (both in magnitude and direction) to the cavitation bubble volume fraction , which is implicated as a critical parameter for useful model predictions.
The model was subsequently applied to liquid aluminium to achieve the prediction of acoustic streaming that considers the effect of cavitating bubbles for the first time. Low bubble volume fractions at low sonication power result in an upward centreline flow, which may be attractive in processing liquid metal in continuous processes. The availability of either upward or downward centreline flow by model and experiment according to sonotrode power may introduce additional useful flexibility on how ultrasonic processing may be applied to liquid metal processing during casting, and in other sono-technologies.