Scholarly article on topic 'Numerical modelling of ultrasonic waves in a bubbly Newtonian liquid using a high-order acoustic cavitation model'

Numerical modelling of ultrasonic waves in a bubbly Newtonian liquid using a high-order acoustic cavitation model Academic research paper on "Chemical engineering"

Share paper
Academic journal
Ultrasonics Sonochemistry
OECD Field of science
{"Acoustic cavitation" / "Numerical acoustics" / "Ultrasonic wave propagation" / "Ultrasonic melt processing" / "Light metal alloys"}

Abstract of research paper on Chemical engineering, author of scientific article — G.S. Bruno Lebon, I. Tzanakis, G. Djambazov, K. Pericleous, D.G. Eskin

Abstract To address difficulties in treating large volumes of liquid metal with ultrasound, a fundamental study of acoustic cavitation in liquid aluminium, expressed in an experimentally validated numerical model, is presented in this paper. To improve the understanding of the cavitation process, a non-linear acoustic model is validated against reference water pressure measurements from acoustic waves produced by an immersed horn. A high-order method is used to discretize the wave equation in both space and time. These discretized equations are coupled to the Rayleigh-Plesset equation using two different time scales to couple the bubble and flow scales, resulting in a stable, fast, and reasonably accurate method for the prediction of acoustic pressures in cavitating liquids. This method is then applied to the context of treatment of liquid aluminium, where it predicts that the most intense cavitation activity is localised below the vibrating horn and estimates the acoustic decay below the sonotrode with reasonable qualitative agreement with experimental data.

Academic research paper on topic "Numerical modelling of ultrasonic waves in a bubbly Newtonian liquid using a high-order acoustic cavitation model"

Accepted Manuscript

Numerical modelling of ultrasonic waves in a bubbly Newtonian liquid using a high-order acoustic cavitation model

G.S. Bruno Lebon, I. Tzanakis, G. Djambazov, K. Pericleous, D.G. Eskin



S1350-4177(17)30084-6 ULTSON 3568

To appear in:

Ultrasonics Sonochemistry

Received Date: Revised Date: Accepted Date:

27 May 2016 14 December 2016 22 February 2017

Please cite this article as: G.S.B. Lebon, I. Tzanakis, G. Djambazov, K. Pericleous, D.G. Eskin, Numerical modelling of ultrasonic waves in a bubbly Newtonian liquid using a high-order acoustic cavitation model, Ultrasonics Sonochemistry (2017), doi:

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Numerical modelling of ultrasonic waves in a bubbly Newtonian liquid using a high-order acoustic cavitation model

G. S. Bruno Lebon1, I. Tzanakis23, G. Djambazov1, K. Pericleous1, D. G. Eskin34

1 Computational Science and Engineering Group (CSEG), University of Greenwich, 30 Park Row, London, SE10 9ET, United Kingdom

2 Brunel Centre for Advanced Solidification Technology (BCAST), Brunel University London, Uxbridge, Middlesex, UB8 3PH, United Kingdom

3 Faculty of Technology, Design and Environment, Oxford Brookes University, Wheatley Campus, Wheatley, OX33 1HX, United Kingdom

4 Smart Materials and Technologies Institute (SMTI), Tomsk State University, Tomsk, 634050, Russia



To address

in treating large volumes of liquid metal with ultrasound, a fundamental

pro me

study of acoustic cavitation in liquid aluminium, expressed in an experimentally validated numerical model, is presented in this paper. To improve the understanding of the cavitation

cess, a non-linear acoustic model is validated against reference water pressure measurements from acoustic waves produced by an immersed horn. A high-order method is used to discretize the wave equation in both space and time. These discretized equations are coupled to the Rayleigh-Plesset equation using two different time scales to couple the bubble and flow scales, resulting in a stable, fast, and reasonably accurate method for the prediction of acoustic pressures in cavitating liquids. This method is then applied to the context of

treatment of liquid aluminium, where it predicts that the most intense cavitation activity is localised below the vibrating horn and estimates the acoustic decay below the sonotrode with reasonable qualitative agreement with experimental data.

Keywords: Acoustic cavitation, Numerical acoustics, Ultrasonic wave propagation, Ultrasonic

melt processing, Light metal alloys

1. Introduction

Significant improvements in the quality and properties of metallic materials are observed when treating them near their liquidus temperature [1 -3]: the beneficial effects of the treatment include the degassing of dissolved gases, improved wetting, activating inclusions by cleaning the solid-liquid interface, enhancing nucleation, and refining the structure of the solidified metal [1, 4]. These improvements are primarily attributed to acoustic cavitation [5]; the term "cavitation" here follows the definition of Neppiras [6] and is restricted to cases involving the formation, expansion, pulsation, and collapse of existing cavities and bubble nuclei. However, treating large volumes of liquid metal, as is required by industrial processes such as continuous casting, is still problematic: the process is time-consuming and volume-limited so it can currently be applied only to a fixed volume of melt in a crucible. To circumvent these difficulties and facilitate the transfer of this promising technology to industry, a fundamental study of melt cavitation treatment is required [7].

Nastac [8] used the 'full cavitation model' [9] developed for hydrodynamic cavitation to model ation structure evolution in an alloy in the presence of ultrasonic stirring, while uting the acoustic field analytically from the Helmholtz reduced wave equation. An mproved version of this model, based on the Keller-Miksis equation [10] and including a turbulent source term arising from the collapse of cavitating bubbles, has been proposed by the authors [11, 12] to model a moving liquid metal volume in a launder. However, the use of a homogeneous cavitation model, e.g. the 'full cavitation model', for acoustic cavitation is questionable. Also, the acoustic solver used in [11, 12] was second order in space and prone

to numerical diffusion; hence a higher-order model [13] is desired to improve the accuracy of the acoustic field prediction. The presence of bubbles significantly alters the ultrasonic wave propagation in the melt, and this influence must be accurately quantified to understand the effect of the acoustic field on a volume of treated metal. In this endeavour, we are proposing a macroscopic cavitation model coupled with a high-order acoustic solver, using reference experiments in water [14] for validation.

A plethora of empirical observations of acoustic cavitation in water is available in the literature. For example, observations of streamers and acoustic Lichtenberg figures have been recorded [15]. The conical bubble structure below the radiating surface of the sonotrode has been observed and studied by many authors [14, 16-19]. The tendency of bubbles to form clusters after collapsing has been observed with a high-speed camera [20]. Pressure measurements with calibrated hydrophones in water under ultrasonic treatment are also available [14, 21].

Alongside this empirical evidence, there exists a series of models attempting to explain and reproduce the bubble cloud behaviour numerically. One class of such models is based on a set of non-linear equations proposed by van Wijngaarden [22] to model wave propagation in a bubbly liquid. Caflisch etal. [23] re-derived this set of non-linear differential equations from the microscopic motion equations of a large number of bubbles. Commander and Prosperetti wrote an extensive review of pressure wave propagation models in bubbly liquids [24] following the insight given by Caflisch et al. [23]. A simplified version of the Caflisch equations has recently been used by Louisnard [25, 26] to model bubble structures below a sonotrode in ter. Other recent advances include the work of Dahnke et al. [27] who modelled the )ustic pressure field in sonochemical reactors with an inhomogeneous density distribution. Vanhille et al. applied and extended a model consisting of a coupled linear non-dissipative wave and volume variation equations [28] based on the Rayleigh-Plesset equation [29] to model the nonlinear propagation of ultrasonic waves in water-air bubble mixtures [30-33]. Servant et al. [34-36] considered the Bjernkes forces [37] in their model of mono- and dual-frequency sono-reactors.

Tudela et al. completed a more recent review [38] of the state of the research and outlined the need to model the nonlinear nature of the problem. They highlighted that the Caflisch-type equations have certain drawbacks for the simulation of the effect of bubbles on strong acoustic fields due to: the non-linear nature of the problem, the limits of using assumptions on bubble sizes and distribution, the assumption of low bubble volume fractions in the derivation of the model, and the applicability only in cases where bubble resonance plays a negligible role [38]. Moreover, this class of models requires extremely small time steps for acoustic pressures higher than the Blake threshold, making it unattractive for the design of experiment simulations that seek optimum parameters to enhance cavitation activity. Despite these drawbacks, resolving the complex coupling between the void fraction and acoustic pressure field is necessary and this class of numerical models is therefore unavoidable in acoustic cavitation modelling.

upled with a i

In this paper, a high-order acoustic model coupled with a cavitation model is presented, followed by validation against acoustic pressures measured in water [14] and then applied to the treatment of aluminium in a cru

crucible. 2. Theory

From the conservation of mass and momentum, and using dp/ dp = c2, the governing equations for sound propagation in a moving fluid are the Navier-Stokes equations in perturbation form:

'bation ь

дР , дР , 2 dvi г !л \

dr+vidrL + ldr = irl + Fi (2)

dt J dxi p dxi dxi i

p denotes pressure and vj are the velocities. p is the liquid density. S contains mass sources, such as due to vibrating solid surfaces and growing (or collapsing) bubbles. In the

model presented in this paper, an additional source term, pc2 d0/ dt, is added to the source term S of equation (1) to account for the acoustic pressure waves induced by the collapse of bubbles, and conversely the sink of acoustic pressure during the creation of bubbles [30].

The forcing terms Fi are usually set to zero for most practical acoustics problems [13] and contain acoustic velocity sources due to a vibrating surface. Ignoring dissipation due to

d"D dv'

viscosity (drij/ dxj term), the convection terms vj ^ and vj ^, and considering a constant

speed of sound, equations (1) and (2) reduce to the standard, linear Helmholtz equation. However, these assumptions are not accurate for modelling acoustic cavitation.

In this implementation, the viscosity and convection terms are retained: this makes equations (1) and (2) fully coupled and non-linear, unlike the linearized cases in [13, 39]. The effect of the flow on pressure predictions is modelled by including the convection terms. The speed of sound in the liquid is given by c = jK/p. However, variations of density and bulk modulus in the bubbly liquid lead to numerical instability due to the discontinuity in derivatives along the saturation curve that separates single phase and two-phase domains [40]. These numerical instabilities can be avoided by treating the speed of sound as a constant, thereby restricting the accuracy of the method to void fractions smaller than 1 % [23]. This assumption is applicable to liquid metals where the bubbles originate from dissolved hydrogen.

The gas volume fraction 0 = ^ nn0R3, where n0 is the number of stationary bubbles of radius

R per unit volume, is calculated from the solution to the Rayleigh-Plesset equation [29] which overns the dynamics of a bubble in the presence of a strong acoustic field:


rr + r2 = 2s, (3)

where ps is given by

Ps = P„(.t) + Pv-Y-4-f-p(.t), (4)

and pb is the pressure inside the bubble. pv is the vapour pressure in the bubble. a is the surface tension between the gas and liquid interface. p is the dynamic viscosity of the liquid.

The use of the Rayleigh-Plesset equation to derive the gas volume fraction assumes the following:

1. The internal pressure of the bubbles is homogeneous, since the i negligible.

2. The bubbles remain roughly spherical. Due to the large value of surface tension of an interface of hydrogen with liquid aluminium, bubbles observed during melt cavitation are small, in the region of 10-100 |im in radius [41].

3. Modelling

3.1. Wave equations discretization

The discretization method of Djambazov etal. [13] is used to solve equations (1) and (2). The computational meshes are fully staggered, as described in [13] and illustrated in Figure 1. Fully staggering pressures and velocities allows the formulation of a fully explicit, stable, second order accurate scheme [42]. The accuracy can then be extended to higher orders by allowing the scheme to become implicit provided it retains a strong diagonal dominance to ensure fast convergence [39]. The computational domain is divided into regular cells, with the scalar quantities pressure and bubble volume fraction stored in cell centres, and velocity components stored at cell faces in the middle of each time step. In this formulation, curves are represented by a castellated mesh, i.e. with a bitmap grid-like structure. The DRP (dispersion-relation-preserving) scheme [43] is used with its differentiation and temporal integration steps for the convection integrals only.

ertia of the gas is

3.1.1. Spatial derivatives

olume i

Figure 1. Discretisation of pressure and velocities on two-dimensional finite volume cells [39].

The terms pc2 and are evaluated with 6-point stencils for the spatial derivatives. At

OXj p ox^

solid boundaries, values are mirrored to provide the missing points, and the free surface is modelled as a layer of cells with fixed acoustic pressure p = 0 Pa. For radiative boundary conditions, the derivatives are extrapolated from previous time step values [13].

In all other cells, the first derivative of a function f with respect to x is expressed as

The coefficients aj have been optimized by Djambazov [13] to make the scheme exact to the fourth order in space. These coefficients are provided in the appendix.

Since the fluids considered here are Newtonian, the second derivatives of velocity are required to determine drij/ dxj. Since viscosity is added as a source term in F, the viscous forces are computed at each iteration: the second derivatives of velocity vt along xj are puted as

d2Vj if n

3.1.2. Temporal integration

The terms pc2^Vj and are computed in the middle of each time step and are then used

OXj p ox^

to update pressure and velocities at the end of their time steps, using the following approximate integration

j \t f{t)dt = AtZ3m=0 bmf(-mAt) (7)

with the coefficients bm chosen to make the scheme third order accurate in time [13]. These coefficients are provided in the appendix.

3.2. Adaptive time stepping for bubble dynamics equation

The Rayleigh-Plesset equation (3) is solved using the Runge-Kutta method with an adaptive time step h< At evaluated as follows [44]:

1. The Jacobian matrix \J] is calculated as

pfl dR

\ dR dR

\ dR dR


2. The eigenvalues X1 and X2 of the amplification matrix [A] = 1 + h\j] are evaluated. The maximum deviation from 1.0 is calculated as S = max (\A1 -1\,\X2 - 1|).

3. The time step h is resized so that S becomes close to a target value denoted by Smax bility, Smax is set to a small value (0.001) so that the maximum Lyapunov

xponent is 1.0 in practice.

Depending on the stage of the bubble cycle, h can vary between 10 s to 10 s. 3.3. Iterative procedure

Starting from an initial guess for the solved variables, the coupled equations (1), (2), and (3) are solved in each time step as follows:

1. The acoustic pressure is solved from equation (1).

2. The Rayleigh-Plesset equation (3) is solved in a separate time scale according to the procedure outlined in Section 3.2. The initial bubble radius and interface velocity in each cell are taken from the previous time step and their values updated with smaller time increments until the flow time step value is reached. The gas volume fractit calculated from the new radii.

3. The velocity components of the pressure perturbation are solved from equation (2) at the end of the time step.

4. The solution is advanced to the next time step and procedures 1 to 3 are repeated until the last time step is reached.

3.4. Material properties

Table 1: Material properties for water and aluminium [1, 45, 46]. Surface tension with air

interface for water. Hydrogen interface for aluminium.

adaptive on

Material Water Aluminium

Temperature (°C) 20 700

Density pi (kg m-3) 1000 2375

Dynamic viscosity p (10-3 Pa s) 1.004 1.0

Kinematic viscosity v (10-6 m2 s-1) 1.0 0.42

Speed of sound c (m s-1) 1482 4600

Surface tension a (N m-1) 0.079 0.860

Vapour pressure pv (kPa) 2.2 negligible

Bulk modulus K (GPa) 2.15 41.2

The material properties used in the numerical simulations are listed in Table 1. The gas phase is assumed to be adiabatic in each case and therefore k = y = 1.4.

3.5. Initial and boundary conditions

The liquid is initially unperturbed (constant hydrostatic pressure and all velocity components set to 0) and contains n0 bubbles of radius R0 per unit volume. In liquid metals, an initial number of nuclei is always assumed since cavitation is attributed to both the

hydrogen-containing inclusions and the dissolved hydrogen that is released from aluminium when the local pressure decreases [1]. The vapour pressure of aluminium at its melting point is 0.000012 Pa [47] and therefore vapour bubbles are unlikely to form in the liquid bulk [48]. Based on the numerical values of acoustic simulations from the literature [30, 31], the number of bubbles per unit volume (bubble density) is n0 = 1 x 1011 m-3 and the initial radius used in water was in the range R0 = [1, 10] ^m. This corresponds to an average distance of 22 radii between bubbles in the extreme case R0 = 10 ^m. This separation is long enough to prevent the motion of bubbles due to the effect of secondary Bjerknes forces [49]. For aluminium, n0 = 1 x 1011 m-3 and R0 = 1 ^m.

A sinusoidal pressure signal is indirectly prescribed below the sonotrode (see Figure 2) by specifying the acoustic velocity at the sonotrode surface. The acoustic velocity amplitude is calculated from the displacement amplitude A of the sonotrode as

v = 2 nfA (9)

per boundary is a free surface from which a 180° phase shift occurs upon reflection of ie acoustic wave: this is approximated by setting p = 0 Pa in the top row of computational cells (representing the atmosphere above the interface). All other boundaries, including the sonotrode walls, are fully reflective to sound and are modelled using the mirroring technique from [13]. Radiative boundary conditions are used to approximate absorbent boundaries. The derivatives in transparent cells at the edges of the domain are updated using a second order interpolation on a 3-point stencil [39].

3.6. Geometry and mesh

The geometry of a water vessel is shown in Figure 2 and corresponds to the setup from Campos-Pozuelo etal. [14]. The liquid depth is 18 cm. The radiating surface of the sonotrode, vibrating at 20 kHz, is 1 cm below the free surface. The sonotrode radius is 3.5 cm. The left and right boundaries are fully reflective to sound, as are the sonotrode walls. The bottom boundary consisting of an absorbent material is modelled as a transparent boundary. The top

boundary is a free surface. The velocity of the horn is provided in [14] and pre: measurements 4 cm below the sonotrode axis are used for validation.


3 2: Geo

undary. Th assure

Free surface boundary

Reflective boundary

<- 18 cm ->

Transparent boundary

Figure 2: Geometry of water vessel in the experimental setup from [14]. The origin (black dot) is taken as the point of intersection between the liquid free surface and the axis of the sonotrode. The hydrophone position (clear dot) is 4 cm below the sonotrode surface.

Three mesh densities were used in the simulation and the bubble dynamics at the monitoring point, corresponding to the hydrophone location in [14], were found to be independent of the grid size Ax. Results are presented in a medium coarse mesh with grid size Ax of 2.0 mm. The same grid size is used in all coordinate directions.

The solution is computed on a 2D Cartesian mesh to make the problem tractable with converged meshes. This restriction makes only qualitative comparisons with experiments possible. Full 3D computations are planned in future works after a massively parallel implementation of the acoustic model is completed.

„.C^-.«. F—

step and grid size in each case are chosen such that the Courant number is always less than 0.2. Below a Courant number of 0.2, the computed pressures at the validation point are identical: all results are therefore presented with a Courant number of 0.1, corresponding to time steps of the order of ^s.

Another simulation is then run for the case of liquid aluminium in a crucible as depicted in Figure 3, corresponding to the setup available at Brunel University London [50]. The crucible walls are fully reflective and a 180° shift occurs upon reflection from the free surface. The liquid height is 17.5 cm, the radius of the cylindrical base is 6 cm corresponding to a charge of 5.2 kg of commercially pure aluminium at 700 °C. The transducer operates at 17.7 kHz and 3.5 kW input power, the sonotrode tip (20 mm in diameter) is immersed 20 mm below the free surface. The displacement of the horn is calculated from the operating power [51].




Figure 3: Schematic of aluminium treatment setup [50]. The origin is the axis of the vibrating surface of the sonotrode. Clear dots represent (numerical) probe positions.

4. Results and discussion

! origin is rical) prob

This section describes the comparison of predicted pressures in water with experimental data and the profile of the predicted bubble cloud below the horn. Following the comparison with water, the aluminium sonication case is presented with qualitative comparison with experimental data.

s using th

4.1. Water

Simulations using the high-order acoustic model are compared with experimental pressure measurements from [14]. Figure 4 shows the predicted and measured pressure evolution at hydrophone position 4 cm below the sonotrode.

the hy

03 fin






Campos-Pozuelo et al. [14] Numerical prediction

S 2.00E+05






Figure 4: Comparison between numerical results and pressure measurements 4 cm below the sonotrode.

The maximum predicted press

ure of 570

kPa is not significantly far from the maximum

recorded pressure of 610 kPa. A pressure wave is also emitted at each bubble collapse. Both the peak pressures and the negative pressures are of the same order of magnitude for both numerical predictions and measured values. An exact realization of the experimental data is not possible, since the exact operating conditions (including position of initial nuclei, roughness of surface of vessel, precision in the location of the sonotrode within the mesh resolution ...) cannot be possibly determined: this is also why cavitation pressure

easurements appear chaotic. Nevertheless, the broad features of the cavitation dynamics, namely the peak pressures and intervals, are correctly predicted.

£ 0.40

Numerical Prediction Campos-Pozuelo et al. [14]

0 10 20 30 40 50 60 70 80 90 100

Frequency (kHz)

Figure 5: Comparison between experimental


Srum an

and numerical prediction.

Figure 5 shows the comparison between the numerical and experimental spectra. There is a general good agreement with the predicted frequencies. The strong subharmonic peaks are due to the influence of the bubble cloud below the sonotrode that grows and collapse at a rate of roughly % of the forcing frequency, consistent with the experimental supercavitation observation from [52].

The decay of pressure amplitude with distance from the source is reasonably predicted, as shown in Figure 6. The values for pressure amplitude are obtained by applying a low-pass Butterworth filter to the predicted numerical pressures in the computational cells located at is of the sonotrode.

03 fin

1.6E+5 -| 1.4E+5 -1.2E+5 -1.0E+5 -

• Campos-Pozuelo et al. [14] —Numerical Prediction

S3 8.0E+4 -

^ 6.0E+4 -

4.0E+4 -

2.0E+4 -


0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11

Spatial coordinate (m)

res and r

Figure 6: Comparison between measured pressures and numerical prediction at cavitating conditions. The pressure amplitude from the numerical simulation is obtained by filtering the numerical pressures using a Butterworth filter.

Figure 7: Development of the bubble cloud below the sonotrode in water. Light grey contours represent volume fractions of 0.1 % and darker contours represent volume fractions of 0.5 %.

Figure 7 shows the evolution of the bubble cloud structure below the sonotrode in water. The cavity first grows (50 |js) and then collapses from the outer rim to form a toroidal shape (79 |js).

This is a less dramatic effect of the "acoustic supercavitation" effect reported by [52], due to the lower energy densities involved in this experiment. The large conical bubble structure, that is routinely observed in water experiments [53], is established after a period of initial transient although the structure does not keep a perfectly stable shape (as shown from 174 ^s).

A time averaged volume fraction in the last cycle of the simulation is shown in Figure 8. The numerical method does not predict the correct size of the conical structure, with the cone

ith the cor J experirm

occupying the whole sonotrode surface instead of 2/3 of the area as observed experimentally [18]. With a qualitative 2D comparison, it is not possible to infer the source of this discrepancy since 2D images of the conical bubble structure cover the whole plane below the sonotrode.

Figure 8: Time averaged bubble structure below the sonotrode.

1.1. Sensitivity of results to bubble density and initial radius A parametric study of the effect of the choice of initial radii RQ and bubble density nQ on the pressure predictions has been conducted with Dakota [54] with a Courant number of 0.02 for all cases. The multidimensional study was performed with 3 partitions for each variable and the bounds were 1 ^m < RQ < 10 pm (a factor of 10 lower than bubble size observed during

cavitation) and 1 x 1011 m-3 < nQ < 1 x 1012 m-3, values commonly used in the literature for

water [30-33]. The results of the sensitivity analysis are given in Table 2. The pressure predictions are insensitive to variations of bubble density at low initial volumes and low initial radii. However, large initial radii lead to variations in predictions and should be avoided in numerical simulations: these low values are consistent with the small sizes of stable hydroge bubbles in the melt before cavitation [1].

Table 2: Multidimensional parametric study

R0 (p,m) n0 (m-3) Initial volume (m3) Maximum pressure (MPa) Minimum pressure (MPa)

1 1.00E+11 4.19E-07 1.10 -°.23

1 4.00E+11 1.68E-06 1.17 -0.24

1 7.00E+11 2.93E-06 1.23 -0.23

1 1.00E+12 4.19E-06 1.84 -0.26

4 1.00E+11 2.68E-05 1.39 -0.21

4 4.00E+11 1.07E-04 0.87 -0.13

4 7.00E+11 1.88E-04 0.71 -0.11

4 1.00E+12 2.68E-04 0.61 -0.10

7 1.00E+11 1.44E-04 0.84 -0.11

7 4.00E+11 5.75E-04 0.51 -0.08

7 7.00E+11 1.01E-03 0.33 -0.07

7 1.00E+12 1.44E-03 0.35 -0.06

10 1.00E+11 4.19E-04 0.74 -0.09

10 4.00E+11 1.68E-03 0.25 -0.06

10 7.°°E+H 2.93E-03 0.12 -0.04

10 1.00E+12 4.19E-03 0.10 -0.04

4.2. Aluminium

Figure 9: Development of the bubble cloud below the sonotrode in aluminium. Light grey contours represent volume fractions of 0.1 % and darker contours represent volume fractions of 0.5 %.

number of 0.1, n0 = 1 x 1011 m-3, and R0 = 1 ^m. Figure 9 shows the evolution of the bubble

Results for the aluminium sonication case are presented in this section, using a Courant

cloud structure below the horn. After an initial transient period (< 100 ^s), a stable cone-like structure is developed below the sonotrode: this is the region of intense cavitation activity where nucleation sites are generated during melt treatment.

Figure 10 shows the predicted pressure at 5 selected points in the computational domain: 3 points along the axis of the sonotrode at 2 cm, 6 cm and 10 cm below the sonotrode surface, and 2 off-axis positions all 2 cm below the sonotrode surface - at x = 3 cm (corresponding to the midpoint between the axis and the crucible) and x = 6 cm (at the edge of the crucible). The

ustic pressures predicted decrease dramatically with distance away from the sonotrode and away from the axis, with the point 2 cm below the sonotrode axis registering the highest pressure of 9.5 MPa. Maximum values at the selected points are listed in Table 3: Maximum predicted pressures at selected probe positions in 2 ms to 3 ms range. This pressure dependence on distance from the source is in agreement with experiments reported in [55,

aco and

56]. A zoom on the minimum acoustic pressures are shown in Figure 11, with minimum pressures confined to 2 bars.

9.0E+06 -

7.0E+06 -

5.0E+06 -

3.0E+06 -

1.0E+06 -


"x = 0 cm, y = -2 cm" "x = 0 cm, y = -6 cm" "x = 0 cm, y = -10 cm" "x = 3 cm, y = -2 cm" "x = 6 cm, y = -2cm"

2.0E-03 2.2E-03

Figure 10: Acoustic pressure

re predict



03 2.6E-03 2.8E-03

Time (s)


at selected points in aluminium domain.

O.OE+OO -2.0E+04 -4.0E+04 -6.0E+04

fS -8.0E+04

3 -1.0E+05

Ö -1.2E+05

-1.4E+05 -1.6E+05 -1.8E+05


2.0E-03 2.2E-03 2.4E-03 2.6E-03

Time (s)

"x = 0 cm, y = -2 cm"----"x = 0 cm, y = -6

"x = 3 cm, y = -2 cm" "x = 6 cm, y

pressures prec


"x = Ocm, y = -10 cm"

Figure 11 : Zoom on the minimum acoustic pressures predicted in the aluminium crucible. Table 3: Maximum predicted pressures at selected probe positions in 2 ms to 3 ms range.

Position Distance from sonotrode (cm) Maximum pressure (MPa)

x = 0 cm, y = -2 cm 2.0 \ 9.5

x = 0 cm, y = -6 cm 6.0 1.4

x = 0 cm, y = -10 cm 10.0 0.9

x = 3 cm, y = -2 cm 3.6 2.7

x = 6 cm, y = -2 cm 6.3 1.6

03 fin

Numerical prediction Tzanakis et al. [56] 600

y = 0.0268x-L454< R2 = 0.9564

y = 4.09

1 300 j

0.02 0.04 0.06


stance fr(

Figure 12: Variation of maximum pressure with distance from sonotrode and comparison with cavitation intensity from [56].

Figure 12 shows the variation of the maximum pressure with distance from the sonotrode surface. The pressure-distance relationship obeys a power law with an R2 value of 0.95. The decay of pressure with distance is of the order of 1.45 per metre. This pressure dependence on distance is in agreement with qualitative experimental observations with a high-temperature probe [55, 56] plotted on a separate scale, with a decay of the order of 1.28 per metre. The experimental values are quoted in mV since the maximum pressures cannot vered from the experimental data. This large decay is expected as the efficiency in ustic radiation is proportional to the ratio of horn radius to wavelength. The large wavelength in aluminium and the comparatively small sonotrode makes the pressure decrease with distance pronounced.

be reco aco

5. Conclusions

This paper presents a high-order numerical model for predicting acoustic pressures in bubbly liquids subject to acoustic cavitation. This model is computationally attractive for simulations of ultrasonic melt treatment since it is stable for Courant numbers of 0.2 and couples two separate time scales: a bubble dynamics time scale (in the range of 10-20 s to 10-8 s) and the flow time scale (or the order of 1 ^s). The model has been validated by a water experiment presented in the literature and predicted the reported cavitation characteristics with reasonable agreement between pressure magnitudes, spectra, and decay with distance from the source. When applied to the case of sonication of liquid aluminium, the model predicts a reasonable dependence of acoustic pressure with distance from the ultrasonic source.


Spatial discretization coefficients

Aa = 0.047569386

Temporal discretization coefficients

Ah = 0.08307437

al = 8 + Aa

a2 = -2(è+Aa)

a3 = —A„

3 10 a

fco = 1 + Afe

= 3A*-è



The authors are grateful to the UK Engineering and Physical Sciences Research Council

(EPSRC) for financial assistance for this research in contract numbers EP/K00588X/1 and EP/K005804/1. A representative sample of research data is provided in supplementary data at The underlying raw data is not shared online due to its size.


[1] G.I. Eskin, D.G. Eskin, Ultrasonic Treatment of Light Alloy Melts, Second ed., CRC Press, Boca Raton, FL, USA, 2014.

[2] J. Campbell, Effects of vibration during solidification, Int. Mater. Rev., 26 (1981) 71-108.

[3] O.V. Abramov, Action of high intensity ultrasound on solidifying metal, Ultrasonics, 25 (1987) 73-82.

[4] G.I. Eskin, G.S. Makarov, Y.P. Pimenov, Effect of ultrasonic processing of molten metal on structure formation and improvement of properties of high-strength Al-Zn-Mg-Cu-Zr alloys, Adv. Perform. Mater., 2 (1995) 43-50.

[5] G.I. Eskin, Cavitation mechanism of ultrasonic melt degassing, Ultrason. Sonochem., 2 95) S137-S141.

[6] E. Neppiras, Acoustic cavitation series: part one: Acoustic cavitation: an introduction, Ultrasonics, 22 (1984) 25-28.

[7] S. Komarov, K. Oda, Y. Ishiwata, N. Dezhkunov, Characterization of acoustic cavitation in water and molten aluminum alloy, Ultrason. Sonochem., 20 (2013) 754-761.

[8] L. Nastac, Mathematical Modeling of the Solidification Structure Evolution in the Presence of Ultrasonic Stirring, Metall. Mater. Trans. B, 42 (2011) 1297-1305.

[9] A.K. Singhal, M.M. Athavale, H. Li, Y. Jiang, Mathematical basis and validation of the full cavitation model, J. Fluids Eng., 124 (2002) 617-624.

[10] J.B. Keller, M. Miksis, Bubble oscillations of large amplitude, J. Acous. Soc. Am., 68 (1980) 628-633.

[11] G.S.B. Lebon, K. Pericleous, I. Tzanakis, D.G. Eskin, Application of the "Full Cavitation Model" to the fundamental study of cavitation in liquid metal processing, IOP Conf. Ser. Mater. Sci. Eng., 72 (2015) 052050.

[12] G.S.B. Lebon, K. Pericleous, I. Tzanakis, D.G. Eskin, A model of cavitation for the treatment of a moving liquid metal volume, Int. J. Cast Met. Res., (2016).

[13] G.S. Djambazov, C.H. Lai, K.A. Pericleous, Staggered-mesh computation for aerodynamic sound, AIAA J., 38 (2000) 16-21

[14] C. Campos-Pozuelo, C. Granger, C. Vanhille, A. Moussatov, B. Dubus, Experimental and theoretical investigation of the mean acoustic pressure in the cavitation field, Ultrason. Sonochem., 12 (2005) 79-84.

IOP Conf , A model of cavit;

ihille, A


eva [181

[15] R. Mettin, S. Luther, C.-D. Ohl, W. Lauterborn, Acoustic cavitation structures and simulations by a particle model, Ultrason. Sonochem., 6 (1999) 25-29.

[16] O. Dahlem, J. Reisse, V. Halloin, The radially vibrating horn: A scaling-up possibility for sonochemical reactions, Chem. Eng. Sci., 54 (1999) 2829-2838.

[17] B. Dubus, C. Granger, P. Mosbah, A. Moussatov, H. Schmit, L. Sztor, C. Campos-Pozuelo, Numerical modeling of the ultrasonic cavitation field and experimental

luation of the bubble density, in: Forum Acusticum 2002, 2002.

[18] A. Moussatov, C. Granger, B. Dubus, Cone-like bubble formation in ultrasonic cavitation field, Ultrason. Sonochem., 10 (2003) 191-195.

[19] A. Moussatov, R. Mettin, C. Granger, T. Tervo, B. Dubus, W. Lauterborn, Evolution of acoustic cavitation structures near larger emitting surface, in: WCU 2003, Paris, 2003.

j yaa uuuui

Rations fi &

es in bubbly liqi 1989) 732-

[20] J. Eisenerm, J. Schneider, C. Cairos, R. Mettin, Sound films of acoustic cavitation, in: Proc. the 40th Italian (AiA) Annual Conf. on Acoustics and the 39th German Annual Conf. on Acoustics (DAGA), Merano, 2013., 2013.

P1, A. Zn¡darcic, R. Me«in, M- Dular' Modelin9 in a ^ changing p^e « -

application to a small ultrasonic horn, Ultrason. Sonochem., 22 (2015) 482-492.

[22] L. van Wijngaarden, On the equations of motion for mixtures of liquid and gas bubbles, J. Fluid Mech., 33 (1968) 465-474.

[23] R.E. Caflisch, M.J. Miksis, G.C. Papanicolaou, L. Ting, Effective equations for wave propagation in bubbly liquids, J. Fluid Mech., 153 (1985) 259-273.

[24] K.W. Commander, A. Prosperetti, Linear pressure waves in bubbly liquids: Comparison between theory and experiments, J. Acous. Soc. Am., 85 (1989) 732-746.

[25] O. Louisnard, A simple model of ultrasound propagation in a cavitating liquid. Part I: Theory, nonlinear attenuation and traveling wave generation, Ultrason. Sonochem., 19 (2012) 56-65.

[26] O. Louisnard, A simple model of ultrasound propagation in a cavitating liquid. Part II: Primary Bjerknes force and bubble structures, Ultrason. Sonochem., 19 (2012) 66-76.

[27] S. Dahnke, K.M. Swamy, F.J. Keil, Modeling of three-dimensional pressure fields in sonochemical reactors with an inhomogeneous density distribution of cavitation bubbles: Comparison of theoretical and experimental results, Ultrason. Sonochem., 6 (1999) 31-41.

[28] E.A. Zabolotskaya, S.I. Soluyan, Emission of harmonic and combination-frequency waves by air bubbles, Soviet physics. Acoustics, 18 (1973) 396-396.

. Plesset, The dynamics of cavitation bubbles, J. Appl. Mech., 16 (1949) 277-282. ] C. Vanhille, C. Campos-Pozuelo, Nonlinear ultrasonic propagation in bubbly liquids: A numerical model, Ultrasound Med. Biol., 34 (2008) 792-808. [31] C. Vanhille, C. Campos-Pozuelo, Nonlinear ultrasonic waves in bubbly liquids with nonhomogeneous bubble distribution: Numerical experiments, Ultrason. Sonochem., 16 (2009) 669-685.

[29] M.S [30

ncy reactor,

[32] C. Vanhille, C. Campos-Pozuelo, Nonlinear ultrasonic standing waves: Two-dimensional simulations in bubbly liquids, Ultrason. Sonochem., 18 (2011) 679-682.

[33] C. Vanhille, C. Campos-Pozuelo, Acoustic cavitation mechanism: A nonlinear model, Ultrason. Sonochem., 19 (2012) 217-220.

[34] G. Servant, J.P. Caltagirone, A. Gerard, J.L. Laborde, A. Hita, Numerical simulation of cavitation bubble dynamics induced by ultrasound waves in a high frequency reactor Ultrason. Sonochem., 7 (2000) 217-227.

[35] G. Servant, J.L. Laborde, A. Hita, J.P. Caltagirone, A. Gerard, Spatio-temporal dynamics of cavitation bubble clouds in a low frequency reactor: comparison between theoretical and experimental results, Ultrason. Sonochem., 8 (2001) 163-174.

[36] G. Servant, J.L. Laborde, A. Hita, J.P. Caltagirone, A. Gérard, On the interaction between ultrasound waves and bubble clouds in mono- and dual-frequency sonoreactors, Ultrason. Sonochem., 10 (2003) 347-355.

[37] V.F.K. Bjerknes, Fields of force, Columbia Univ., USA, 1906.

[38] I. Tudela, V. Saez, M.D. Esclapez, M.I. Diez-Garcia, P. Bonete, J. Gonzalez-Garcia, Simulation of the spatial distribution of the acoustic pressure in sonochemical reactors with numerical methods: A review, Ultrason. Sonochem., 21 (2014) 909-919.

[39] G. Djambazov, Numerical techniques for computational aeroacoustics, in, University of

:ov, Nui 8.

Greenwich, 1998.

[40] S. Clerc, Numerical simulation of the homogeneous equilibrium model for two-phase

flows, J. Comput. Phys., 161 (2000) 354-375.

[41] W.W. Xu, I. Tzanakis, P. Srirangam, S. Terzi, W.U. Mirihanage, D.G. Eskin, R.H. lathiesen, A.P. Horsfield, P.D. Lee, In Situ Synchrotron Radiography of Ultrasound

Cavitation in a Molten Al-10Cu Alloy, in: TMS2015 Supplemental Proceedings, John Wiley & Sons, Inc., Hoboken, NJ, USA, 2015, pp. 61-66.

[42] K.W. Morton, D.F. Mayers, Numerical Solution of Partial Differential Equations: An Introduction, Cambridge University Press, New York, NY, USA, 2005.

[43] C.K.W. Tam, J.C. Webb, Dispersion-relation-preserving finite difference schemes for computational acoustics, J. Comput. Phys., 107 (1993) 262-281.

[44] E. Shams, J. Finn, S.V. Apte, A numerical scheme for Euler-Lagrange simulation of bubbly flows in complex systems, International Journal for Numerical Methods in Fluids, 67

t of ultras^ esium alloys, Mate

[52 ultr

(2011) 1865-1898.

[45] J.G. Speight, others, Lange's handbook of chemistry, McGraw-Hill New York, 2005.

[46] I. Waller, I. Ebbsjo, On the calculation of the bulk modulus of a liquid metal using the pressure fluctuations in a microcanonical ensemble and a numerical application to liquid aluminium, J. Phys. C: Solid State Phys., 12 (1979) L705.

[47] J.J. Jasper, The surface tension of pure liquid compounds, J. Phys. Chem. Ref. Data, 1 (1972) 841-1010.

[48] Y.-J. Chen, W.-N. Hsu, J.-R. Shih, The effect of ultrasonic treatment on microstructural and mechanical properties of cast magnesium alloys, Materials Transactions, 50 (2009) 401-408.

[49] G.S.B. Lebon, K. Pericleous, I. Tzanakis, D.G. Eskin, Dynamics of two interacting hydrogen bubbles in liquid aluminum under the influence of a strong acoustic field, Phys. Rev. E, 92 (2015) 043004.

[50] I. Tzanakis, G.S.B. Lebon, D.G. Eskin, K.A. Pericleous, Characterisation of the ultrasonic acoustic spectrum and pressure field in aluminium melt with an advanced cavitometer, Journal of Materials Processing Technology, 229 (2016) 582-586.

[51] M.M. van Iersel, N.E. Benes, J.T.F. Keurentjes, Importance of acoustic shielding in sonochemistry, Ultrason. Sonochem., 15 (2008) 294-300.

] A. Znidarcic, R. Mettin, C. Cairós, M. Dular, Attached cavitation at a small diameter ultrasonic horn tip, Phys. Fluids, 26 (2014) 023304.

[53] B. Dubus, C. Vanhille, C. Campos-Pozuelo, C. Granger, On the physical origin of conical bubble structure under an ultrasonic horn, Ultrason. Sonochem., 17 (2010) 810-818.

[54] B.M. Adams, L.E. Bauman, W.J. Bohnhoff, K.R. Dalbey, M.S. Ebeida, J.P. Eddy, M.S. Eldred, P.D. Hough, K.T. Hu, J.D. Jakeman, J.A. Stephens, L.P. Swiler, D.M. Vigil, T.M.

Wildey, Dakota, A Multilevel Parallel Object-Oriented Framework for Design Optimization, Parameter Estimation, Uncertainty Quantification, and Sensitivity Analysis: Version 6.0 User's Manual, in, Sandia Technical Report SAND2014-4633, July 2014, 2015.

[55] I. Tzanakis, G.S.B. Lebon, D.G. Eskin, K. Pericleous, Effect of input power and temperature on the cavitation intensity during the ultrasonic treatment of molten aluminium, Transactions of the Indian Institute of Metals, 68 (2015) 1023-1026.

[56] I. Tzanakis, G.S.B. Lebon, D.G. Eskin, K. Pericleous, Investigation of the factors influencing cavitation intensity during the ultrasonic treatment of molten aluminium, Materials and Design, 90 (2016) 979-983.

1. A higher order numerical scheme is used in the context of acoustic cavitation for the computation of acoustic pressure propagation in bubbly Newtonian liquids.

2. The model is stable for Courant numbers of less than 0.2 and couples two different time scales for bubble dynamics and fluid flow: these make computations computationally attractive for modelling ultrasonic melt treatment.

3. This model is validated used a water experiment from the literature and corr the trends observed during ultrasonic treatment: pressure peaks, acousti pressure decay with distance from the source.

4. The model is applied to liquid aluminium and predicts the large decay of acoustic pressure with distance from the radiative source that has been observed in aluminium experiments.