Numerical modelling of acoustic streaming during the ultrasonic melt treatment of direct-chill (DC) casting

Acoustic streaming and its attendant effects in the sump of a direct-chill (DC) casting process are successfully predicted under ultrasonic treatment for the first time. The proposed numerical model couples acoustic cavitation, fluid flow, heat and species transfer, and solidification to predict the flow pattern, acoustic pressure, and temperature fields in the sump. The model is numerically stable with time steps of the order of 0.01 s and therefore computationally attractive for optimization studies necessitating simulation times of the order of a minute. The sump profile is altered by acoustic streaming, with the slurry region depressed along the centreline of the billet by a strong central jet. The temperature gradient in the transition zone is increased, potentially interfering with grain refinement. The cooling rate in the sump is also altered, thereby modifying the dendrite arm spacing of the as-cast billet. The relative position of the sonotrode affects the sump profile, with the sump depth decreased by around 5mm when the sonotrode is moved above the graphite ring level by 100mm. The acoustic streaming jet penetrates into the slurry zone and, as a result, the growth direction of dendritic grains in the off-centre position is altered.


Introduction
Ultrasonic melt treatment is applied to direct-chill (DC) casting for degassing, reducing the macrosegregation level, and refining the grain structure. Ultrasonic processing is performed by dipping one or several sonotrodes into the sump of the billet, as shown in Fig. 1. Eskin and Eskin attributed the grain refining effect of ultrasound to the activation of substrates by wetting, deagglomeration and dispersion of nucleating particles, and dendrite fragmentation, all associated with acoustic cavitation [1].
Numerical modelling of conventional DC casting is popular in the literature. DC casting models can be broadly classified as multiphase and continuum models. In multiphase models, the interfaces between the different phases are explicitly tracked. Ni and Beckermann [2] described a two-phase model based on volume averaging that enables the coupling between microscopic and macroscopic phenomena, but this approach is computationally expensive for optimization purposes.
Bennon and Incropera [3] avoided the requirement of tracking phase interfaces by adopting a continuum formulation that integrates the microscopic description of transport behaviour. Their work has been clarified by Prescott et al. [4] who re-derived the continuum momentum equation and established the need for accurate closure rules to the continuum model. Vreeman et al. [5] later incorporated the transport of free floating dendrites in the slurry to this continuum model and modelled macrosegregation in DC casting. In a separate paper, Vreeman et al. [6] obtained a good comparison between empirical temperature and sump profiles with their continuum model: this approach is the starting point of our Computational Fluid Dynamics (CFD) model. Due to the low casting speeds and slow flows due to natural convection, the flow models in the sump are all laminar.
Mean-field models that closely couple solidification growth at the microscopic level and the macroscopic multiphase flow have been recently developed, both in DC casting and other processes such as vacuum arc re-melting [7]. Heyvaert et al. [8] adopted a three-phase model of grain growth consisting of the solid, inter-dendritic liquid, and extra-dendritic liquid to model macrosegregation and the effect of grain refiners in DC casting. While their model provides a better understanding of coupling of microstructure and macrosegregation at the process scale, more work is required to adequately explain the macrosegregation mechanism. Tveito et al. [9] developed a simplified threephase model that is applicable to equiaxed solidification and established that grain morphology must be correctly described to obtain accurate macrosegregation predictions. These models are not currently considered in this work but will be implemented at a later stage to quantify the effect of ultrasound on macrosegregation.
While DC casting modelling has been successful in the literature, acoustic streaming modelling is more difficult, and even more so in the presence of acoustic cavitation and turbulence. Acoustic cavitation models are based upon the set of equations derived by van Wijngaarden [10] which developed a set of non-linear equations to model flow in bubbly liquids in the presence of moderate pressure oscillations. While van Wijngaarden developed these equations through physical reasoning, Caflisch et al. [11] mathematically re-derived these equations based on Foldy's approximation [12]. This model is valid at slow flow fields because it neglects convection and assume that the bubbles are disperse, i.e. is valid far from the sonotrode. Lebon et al. [13] used such a model to compute the acoustic pressures in water and aluminium and obtained a good agreement with measured values. However, this set of non-linear equations is computationally expensive to solve and require the solution of ordinary differential equations that describe bubble dynamics in each computational cell. This requirement makes the Caflisch model inadequate for optimization studies with the current computational power available. Recent advances to model the non-linear effect of bubbles on sound propagation has led to computationally tractable formulations. The starting point of these types of model is from the linear model of Commander and Prosperetti [14]. Their linearized model incorporated the mixture effect through a complex wave number in a Helmholtz equation. However, this model is applicable at low pressures where the bubbles oscillate linearly. Non-linear effects due to acoustic cavitation were considered by Louisnard [15] who used the Rayleigh-Plesset equation to derive the attenuation term in the complex wave. This term depends on pressure making the pressure propagation equation nonlinear. This method was extended by Jamshidi and Brenner who used the Keller-Miksis equation instead to consider compressibility in the bubble dynamics: this effect cannot be neglected because energy dissipation due to acoustic radiation is of the same order of magnitude as thermal dissipation [16]. The non-linear and linear models were compared by Dogan and Popov [17] who established that the non-linear models more adequately represent the effect of attenuation in bubbly liquids. While Louisnard [15] assumed that the real part of the wave number was given by the linear approximation, Trujillo [18] rigorously re-derived a non-linear model and validated his formulation at low pressure amplitudes. His derivation revealed that the real part of the wave number is related to the average acoustic energy while the imaginary part is related to average energy dissipation.
These recent development in non-linear pressure propagation theory has enabled the emergence of an acoustic streaming model that properly accounts for the effect of cavitation bubbles. Louisnard [19] coupled his non-linear model to the Reynolds-Averaged Navier Stokes equations and modelled turbulent flow in the presence of cavitation and obtained good agreement with a Particle Image Velocimetry (PIV) experiment in water. In this paper, we apply Louisnard's acoustic streaming model [19] coupled with Trujillo's non-linear Helmholtz equation [18] to the ultrasonic treatment in DC casting (USDC) of an AA6XXX series aluminium alloy to predict the acoustic streaming pattern in the sump. This model extends our previously validated model of acoustic streaming in water [20] by including the heat transfer and species conservation equations and considering the effect of acoustic radiation in non-linear pressure propagation. Results are presented for two configurations: with the sonotrode either submerged to the level of the graphite ring or positioned 100 mm above the graphite ring level, for which experimental measurements of grain size and dendrite arm spacing are available [21]. Both configurations result in larger temperature gradients across the reduced width of the transition region and larger cooling rates near the centre of the billet. Results show that the flow effect is weaker for the higher position, resulting in less melt penetration into the slurry.

Acoustic streaming model
Denoting the harmonic part of acoustic pressure p as Pe ( ) iωt R , the complex amplitude P is approximately described by the nonlinear Helmholtz equation [19] ∇ + = P K P 0, 2 2 (1) where the real and imaginary parts of K 2 are given by where ω is the angular frequency, c is the speed of sound in the pure liquid, and the terms A and B derived by Trujillo as [18] as either ρ is the pure liquid density. τ denotes time within one period between [0, 2π]. The bubble volume fraction β is given by where R is the bubble radius and = V πR 4 3 3 is the bubble volume. The bubble density N is assumed to follow the step function where the Blake threshold = ⎡ ⎣ being the bubble equilibrium radius.
The terms A and B are estimated from the solution of a bubble dynamics equation. This work follows the approach of Trujillo [18] by adopting the Keller-Miksis equation to account for dissipation due to acoustic radiation: this effect is considered to yield a periodic solution to the evolution of bubble dynamics: σ is the surface tension between the liquid and gas phases, μ is the dynamic viscosity of the liquid, p 0 is the pressure at infinity (set to atmospheric pressure), A is the pressure amplitude (normalized by p 0 ) of the excitation source of angular frequency ω, and p v is the vapour pressure.
The gas pressure p g is evaluated by solving the differential equation which takes into account the effect of heat transfer during bubble pulsation [16,22]. k is the heat conductivity of the hydrogen gas. The gas pressure at the equilibrium radius R 0 , denoted by p g,0 , is used as the initial value for Eq. (11). Assuming adiabatic pulsation, the polytropic exponent is = γ 1.4, the ratio of specific heats. Since the vapour pressure of aluminium at its melting point is 0.000012 Pa [23], aluminium vapour bubble formation is highly unlikely and the vapour pressure for modelling purposes can be approximated as zero. Therefore, no vapour transfer equation is coupled with the Keller-Miksis equation. Also, rectified diffusion of hydrogen bubbles is a slow process [1] as evidenced by empirical observation of stably cavitating bubbles by X-ray radiography [24]. Therefore, the transfer of hydrogen into cavitating bubbles is also neglected.
Employing the method of Toegel et al. [25], the temperature gradient at the bubble surface, required in the evaluation of Eq. (11), is approximated linearly as where D is the diffusivity of the gas [16]. The temperature of the liquid bulk ∞ T is approximated as the inlet temperature. The temperature of the gas inside the bubble, T, is evaluated using the first law of thermodynamics s the thermal diffusion length and C v is the specific heat capacity of the gas.

DC casting model
A continuum formulation is used to present the DC casting problem. The mass conservation equation is where u is the velocity of the liquid phase. The energy balance equation is where = h C T p is the enthalpy, κ is thermal conductivity, T is temperature, C p is specific heat capacity, L f is latent heat of fusion, and f l is the volume fraction of liquid. The source term in Eq. (15) is due to phase change [26].
The species conservation equation is given by where u s is the velocity of the solid shell which is set as the casting speed, C s is the concentration of species s, and D l s is the diffusivity of species s in the liquid. The liquid concentration is calculated using the lever rule where k p is a binary partition coefficient. It is conventional to divide the transition region (between liquidus and solidus) into a slurry (above the coherency isotherm) and a mush (below the coherency isotherm) [5]. The momentum conservation equation in the liquid and slurry region ( ≤ ≤ g f 1 c l ) is given by where = + μ μ μ t lm , is the effective viscosity, p is pressure, and g is the acceleration due to gravity. f represents the force driving acoustic streaming, and g c is the liquid fraction reflecting the coherency, where = ∇ v ρω P is the acoustic velocity that is estimated by solving the equation for sound propagation only [27]. Once this acoustic velocity is estimated, the flow velocity u can be calculated by solving Eq. (18). The overbar indicates that the values are obtained from averaging over a period of the acoustic bubble. The buoyancy term is evaluated in assuming the Boussinesq approximation, i.e.
where β T is the thermal expansion coefficient and β s is the solution expansion coefficient for species s.
In the slurry region, the viscosity is modified to simulate flow with resistance due to the presence of the grains where F μ is a switching function and A c is a crystal constant [28].
In the mushy zone and solid regions ( ≤ ≤ f g 0 l c ), the momentum conservation equation is given by where K is the permeability coefficient. The last term of Eq. (22) is a Carman-Kozeny source term that accounts for the resistance to flow in the mushy/solid region. The − k ω shear stress transport (SST) model is used for closure: The turbulent viscosity is given by k is the kinetic energy of turbulence. ω t is the dissipation rate. The turbulent model coefficients are identical to those in the original reference [29].

Geometry of ultrasonic melt treatment in DC model
Fig . 1 illustrates the ultrasonic treatment process. A 24 mm Ø sonotrode is inserted along the axis of a 155 mm Ø DC casting mould and introduces power ultrasound in the sump. The inlet is set to 50 mm above the sonotrode tip: Reese [30] established that the temperature is stratified in the sump, and as such the inlet temperature can be fixed to the liquidus temperature of the melt and the feeding does not have to be explicitly modelled. The graphite ring depth is 40 mm, followed by a 17 mm deep water cooled aluminium mould. The water jet flow rate is 60 L/min with an average temperature of 20°C. The sonotrode is connected to a magnetostrictive transducer supplying 2.0 kW of power at = f 17.3 kHz: this corresponds to an estimated amplitude of = y 20 µm peak-to-peak. The sonotrode is located at two positions in the sump: (i) submerged to the level of the graphite ring and (ii) at 100 mm above the graphite ring level. This setup corresponds to the experiment described in [21].

Material properties and model parameters
The material properties for an A6060 alloy (composition in Table 1 and solutal properties in Table 2) were calculated using the National Physical Laboratory's (NPL) Virtual Measurement System (VMS) [31] and the Pro-CAST material property calculator [32]. For the specific heat capacity, the values in the transition region were calculated using the method of mixtures. Other parameters are listed in Table 3.

Numerical implementation
The finite volume solver buoyantPimpleFoam from the open source package OpenFOAM 5.x [33] was modified as described in this section. The discretization schemes and solver control parameters are listed in Table 4. The boundary conditions for the model are presented in Fig. 2 and Table 5. The sonotrode is sufficiently far from the solidification front so that the expected intensive cavitation zone does not interact with growing dendrites.

Secondary cooling heat transfer boundary condition
The heat transfer coefficient at the mould and water-cooled mould is estimated from Rohsenow's formula [34] and entered as a table at the wall. The tabular boundary condition is implemented by modifying the compressible: externalWallHeatFluxTemperatureFvPatchScalarField class to read an interpolation table. Material properties for water are given in Table 6.
The heat transfer coefficient is evaluated as where the heat transfer at incipient boiling is given by The forced-convection heat transfer coefficient is evaluated as [36] ⎜ where = Pr C μ k p is the Prandtl number.
The effect of nucleate boiling is accounted for by

Solution of the nonlinear Helmholtz equation
Special consideration is given to the implementation of the solution of the Helmholtz Eq. (1), which is not straight-forward: 1. The Keller-Miksis equation (10), gas pressure equation (11), and first law of thermodynamics equation (13) are solved for a cavitating hydrogen bubble with equilibrium radius = R 5 0 µm for a range of pressures including the pressure below the sonotrode at the operating power. The material properties for the hydrogen gas are given in Table 7. This pressure is estimated from measurements in experiments featuring a similar transducer using a calibrated hightemperature cavitometer [37].

Results and discussion
An axisymmetric model of the USDC setup described in Fig. 1 is run in the custom OpenFOAM solver whose implementation was described in the previous section. This model has been implemented in version 5.x. The model is first run without the sonotrode, but with the same operating conditions, to obtain the initial conditions for the USDC simulations. The results from this conventional DC casting simulation are also compared with the USDC results.
The period-averaging assumption to calculate A and B makes the proposed model computationally cheap, since this stage negates the prohibitive use of very fine time steps that are required in van Wijngaarden-type models. The current results have been run with time steps of the order of 0.01 s, resulting in run times of around 9 h until convergence on a 16 core (3.0 GHz) CPU. These run times make this model attractive for optimization and uncertainty quantification studies.

Evaluation of attenuation terms
The coupled Eqs. (10), (11), and (13) are solved for hydrogen bubbles of equilibrium radius 5 µm cavitating due to sinusoidal forcing signals of frequency 17.3 kHz and amplitudes in the range ≤ ≤ A 0. 3 10. A bubble density of = N 0 1.9 × 10 9 is assumed, corresponding to the range of bubble fraction β for which the numerical model is stable [20]. This system of non-linear equations is solved using the ODE solver supplied by the SciPy Python library [38]. Because of the stiffness of the problem, the Adams/BDF method with automatic stiffness detection and switching [39] is used.
A and B are evaluated using Eqs. (4)- (7). The integrals are evaluated in the last cycle: to ensure that the evaluation is precisely performed in the last cycle, the ODE solution is output at 400 regularly spaced intervals per cycle and only the solution in the last cycle comprising the 400 last values of R and Ṙis used. The variation of the attenuation terms, normalized by bubble volume fraction, with forcing amplitude A is given in Fig. 3. The required values of A and B for the ultrasonic processing simulation are interpolated using splines of order 3 assuming A = 2.4, the expected average pressure under the horn based on previous measurements in aluminium processing [13].

Mesh convergence analysis
The converged mesh density for the numerical simulations is determined by running the DC casting model for three different mesh densities and evaluating an estimate of the continuum solution of the   temperature field using the Richardson Extrapolation method. Fig. 4 shows the grid refinement analysis for the case using the predicted centre line temperature: the solution evaluation on a mesh of 20,736 cells, corresponding to an average cell length of 1 mm, is grid independent. This mesh density is used for all the results presented in this work.

Treatment with the sonotrode aligned with the graphite ring
Fig . 5 shows the effect of ultrasonic treatment on the sump when the sonotrode is submerged down to the level of the graphite ring (as in Fig. 2) by comparing the sump profile in conventional DC casting (left) with the modified profile with acoustic streaming (right). A strong central acoustic streaming jet (see Fig. 6 left) depresses the liquidus and shortens the transition zone in the centre of the billets, thereby drastically increasing the temperature gradient in the transition region. The flow is opposite to the natural convection direction and the melt flows upwards, parallel to the solidification front towards the mould. The transition zone is also depressed at the sides, leaving a larger  with amplitude y assumed to be 10 µm.

Table 6
Water properties at saturation temperature [36].  Table 7 Hydrogen gas bubble properties [1]. Unavailable properties were approximated by those of air.
The acoustic pressure decays exponentially below the sonotrode, as shown in Fig. 6 (right) and Fig. 7. This prediction is consistent with pressure measurements in a crucible [40]. Due to acoustic shielding and energy dissipation of the acoustic wave, most of the ultrasound energy is consumed inside the cavitation zone [41]. With the cavitation zone being active only below the sonotrode, this implies that the acoustic streaming pattern could be the main mechanism behind the grain morphology modification in the sump.
In Fig. 9, the cooling rate T, defined here as ∇ − u u T·( ) s , inside the slightly elevated sump is larger in the centre of the billet due to the large speed jet impinging directly from under the sonotrode and the larger temperature gradient in the shortened transition zone. This increase implies an increase in the solidification velocity. These large cooling rates result in smaller dendrite arm spacing λ [28] compared with the DC casting case without ultrasound.

Treatment with the sonotrode positioned 100 mm above the graphite ring
With the sonotrode elevated 100 mm above the graphite ring level, the flow pattern is similarly reversed at the centre, as shown in Fig. 10. However, the sump becomes somewhat shallower as compared with the sonotrode at the graphite ring level (Fig. 11). Therefore, with a weaker penetration jet, solidification occurs faster than with the sonotrode at the lower position. The temperature gradients are of the same order of magnitude (Fig. 12).
When the sonotrode is located 100 mm above the graphite ring level, the predicted cooling rates are lower than those obtained when the sonotrode is aligned with the graphite ring level, as shown in Fig. 13. This implies that the dendrite arm spacing at the centre of the billet will be larger when the horn in the elevated position, but still smaller than in the reference case.

Verification of simulation results with the microstructure observations
A cast billet of an AA6XXX-series alloy was cut in the horizontal direction in the middle of the cast length. The observed samples are taken from the off-centre and half-radius locations (schematically illustrated in Fig. 1). Details of the experiment can be found elsewhere [21].
Left parts of Figs. 15 and 16 show the grain morphology at half radius of the billet when the ultrasonic sonotrode is located at the graphite ring level (Fig. 15) and 100 mm above the graphite ring level, respectively (Fig. 16). When the sonotrode is located at the lowest position in the hot top, closer from the solidification front, a non-dendritic, very fine grain structure is obtained at half radius. The strong liquid flow penetrating into the semi-solid region predicted by our simulations in Fig. 5 could explain that dendrite fragments and free grains can be transported from the centre of the sump toward the half radius of the billet resulting in the small grain size. Both cases show a clear grain refinement compared with the reference cast without ultrasound (Fig. 14). Grain refinement is more pronounced in Fig. 15, consistent with the stronger velocities, deeper penetration into the slurry region (which facilitates fragmentation and transport of fragments) and higher cooling rates (that facilitate structure refinement) demonstrated in the numerical model for the lower position of the sonotrode (Fig. 9).
Note the interesting dendritic growth towards the sonotrode position in Fig. 15 (right). The numerical simulations predicted the acoustic streaming pattern as shown in Fig. 5. When the strong jet penetrates the transition zone, hot melt is brought into the solidification front from the sonotrode, lowering the liquidus position in the sump and dramatically increasing the temperature gradient at the solidification front.
Studies of dendritic grains growing under the conditions of forced flow has been explored experimentally, as well as numerically by phase field simulations [42]. They conclude that dendrites will grow upstream towards the fluid flow, developing elongated dendritic grains. There is a clear shift from unconstrained to constrained growth [43], leading to more elongated grains, when the sonotrode is at its lowest position thus increasing the temperature gradient at the solid-liquid interface. In situ X-ray investigation of aluminium alloy solidification showed that when the cooling rate is not high enough, grains tend to elongate in the direction of the temperature gradient [44].
The dendrites shown in Fig. 15 grew preferentially towards the incoming melt flow created by the acoustic streaming. However, this was  not pronounced when the sonotrode was located 100 mm higher, with the incoming jet and the flow parallel to the solidification front being considerably weaker than when the sonotrode was placed at the graphite ring level (Fig. 11). The higher temperature gradient in the centre of the billet would also reduce the solidification time at this location in the billet, generating significantly smaller secondary dendrite arm spacing λ 2 than the reference cast as shown in Fig. 17 (measurements from the centre of the billet) and elsewhere [21]. These dendrite arm spacing measurements are in agreement with the prediction of larger cooling rates in the sump in Figs. 9 and 13.

Conclusions
A novel acoustic streaming model applicable to the ultrasonic melt processing integrated in the direct-chill casting process has been presented for the first time. The model is numerically stable, representative of the physics of acoustic cavitation in the melt, and is in qualitative agreement with the grain morphology modification of ultrasonic treatment of billets:   1. The model is computationally attractive: the period averaging assumption negates the use of fine time steps as is required in van Wijngaarden type models. This makes the proposed model suitable for optimization studies. 2. Ultrasound modifies the sump profile by depressing the slurry region along the axis of the billet and pushing the slurry sideways along the solidification front. This contributes to the transport of floating grains towards the melt and increases the temperature gradient in the phase transition region.
3. The sonotrode position affects the acoustic streaming pattern and grain morphology accordingly. A higher sonotrode position within the hot top position leads to a slightly raised sump profile, with no elongated grains at the centre of the billet compared with the sonotrode aligned with the graphite ring level.

Data availability statement
The processed data required to reproduce these findings are available to download from https://dx.doi.org/10.17633/rd.brunel. 7610924.

Appendix A. Supplementary data
Supplementary data to this article can be found online at https://