Modeling of Three-Dimensional Linear Pressure Fields in

A new model is presented for the numerical calculation of pressure fields in liquids with an inhomogeneous distribution of cavitation bubbles. To calc...
2 downloads 4 Views 510KB Size
848

Ind. Eng. Chem. Res. 1998, 37, 848-864

Modeling of Three-Dimensional Linear Pressure Fields in Sonochemical Reactors with Homogeneous and Inhomogeneous Density Distributions of Cavitation Bubbles† Sascha Da1 hnke and Frerich J. Keil* Department of Chemical Engineering, Technical University of HamburgsHarburg, Eissendorfer Strasse 38, D-21071 Hamburg, Germany

A new model is presented for the numerical calculation of pressure fields in liquids with an inhomogeneous distribution of cavitation bubbles. To calculate the pressure field in a homogeneous single-phase fluid, the Helmholtz integral and the Kirchhoff integral are solved numerically. The Helmholtz integral equation and the Kirchhoff integral are used for the calculation of the acoustic field in a homogeneous fluid for all kinds of transducers of various shapes. The first term of the integral equation embodies a simple superposition of the pressure fields of several point sources, which serves to simulate a harmonic vibrating surface, while the Kirchhoff integral calculates the pressure field which emerges from the boundaries. With a new technique the three-dimensional time-independent pressure field is calculated gradually in the beam direction. With this procedure one is able to combine the Helmholtz integral with a wave propagation in liquids with inhomogeneous distributions of cavitation bubbles. Compared to a single-phase fluid, gas bubbles in a liquid lead to a heavy change of phase velocity and sound attenuation. These changes are determined and considered for every step in the beam direction. With this technique, one should be able to calculate the pressure field in a sonochemical reactor with a sufficient approximation which serves to predict the spatial distribution of cavitation events. These events are related to the yield of chemical reactions. 1. Introduction Ultrasound has been used in chemical synthesis for more than 50 years. Ultrasound provides a form of energy for the modification of chemical reactivity which is different from, for example, heat, light, and pressure. Mason and Lorimer (1988) give a list of possible benefits from the use of power ultrasound in chemistry: (a) A reaction may be accelerated or less forcing conditions may be required if sonication is applied. (b) Induction periods are often significantly reduced as are the exotherms normally associated with such reactions. (c) Sonochemical reactions often make use of cruder reagents than conventional techniques. (d) The number of steps which are normally required in a synthetic route can sometimes be reduced. (e) In some situations a reaction can be directed to an alternative pathway. Although a lot of work was done on the investigation of the chemical effects of power ultrasound, only a few efforts were directed toward engineering aspects, like optimum design and scaleup of sonochemical reactors. The physical aspects may roughly be divided into two subjects: (1) calculation of pressure fields and cavitation phenomena within sonochemical reactors; (2) physical interpretation of sonochemical phenomena of the chemical reactions. For this purpose quantum chemical methods have to be used. * Author to whom correspondence is addressed. Fax: +4940-7718-2145. E-mail: [email protected]. † This paper is dedicated to Prof. K. Roel Westerterp on the occasion of his 70th birthday.

The present paper is related to the first class of problems. The subsequent results reveal that there is a need for standardized sonochemical reactors and other equipment in order to be able to compare the results of yield and selectivity measurements. These results are strongly determined by the pressure fields within the reactors. One of the most important mechanisms for physical, chemical, and biological effects of ultrasound is cavitation. Cavitation may be defined as the formation and activity of vapor or gas bubbles in a liquid under any form of stress. The word “formation” refers both to the creation of a new cavity or to the expansion of a preexisting one to a size where it can be observed. One may distinguish between four different kinds of cavitation according to the effects which are responsible for their occurrence: hydrodynamic, acoustic, optic, and particle cavitation. For sonochemical applications only hydrodynamic and acoustic cavitation are significant. Hydrodynamic cavitation is produced by pressure variations in a flowing liquid due to the geometry of the system, whereas acoustic cavitation is produced by sound waves in a resting fluid due to pressure variations. The modeling of acoustic cavitation in sonochemical reactors is the purpose of this paper. When a liquid is irradiated with ultrasound, it causes a series of compression and rarefaction cycles which create areas of high and low local pressures. During the rarefaction phase of the pressure wave a cavitation bubble will emerge if the local underpressure reaches a specific value (“cavitation threshold”) (Apfel (1984), Trevena (1984)), which depends on the tension and the temperature of the liquid. The emergence of cavitation bubbles is possible due to the existence of so-called

S0888-5885(97)00339-4 CCC: $15.00 © 1998 American Chemical Society Published on Web 01/31/1998

Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 849

“cavitation nuclei”. In liquids, especially in water, there are three types of such nuclei. First, there is evidence for the existence of a large number of spherical gas bubbles which may be stabilized against gaseous diffusion by a skin of varying permeability (Yount et al. (1984), Yount (1982)). The varyingpermeability (VP) model has been proposed, because otherwise the proven existence of gaseous bubbles with the appropriate radii in a liquid could not be explained. Gas phases larger than the order of 10-6 m should rise to the surface of a standing liquid, and smaller ones should dissolve rapidly via the outward diffusion of gas that results from surface tension. Second, there may be solid particles (motes) in the liquid with gas trapped in crevices of these particles (Apfel (1984), Miller and Nyborg (1983), Miller (1982), Apfel (1970)). Bubbles might be stabilized with regard to buoyancy by attachment to a mote. Third, gas can be trapped in tiny crevices in the walls of the vessel containing the liquid. Two distinct types of bubble motion are possible. The first one is a “stable” motion in which the bubbles oscillate for many periods of the sound field, whereas in the second one the bubbles will exist for less than one cycle. On collapse, transient cavitation bubbles produce very high pressures and temperatures, which cause such phenomena as erosion, emulsification, molecular degradation, sonoluminescence, sonochemical, and biological effects. Temperatures approaching 5000 K have been determined, and pressures of several hundred atmospheres have been calculated. A lot of theoretical work has been done concerning the calculation of bubble dynamics and the dependent physical and chemical effects (Leighton (1994), Young (1989), Brennen (1995), Luche (1993), Gaitan et al. (1992)). Despite all this work, the prediction of the behavior of a given sonochemical reactor system has been tried only in a few cases (Naidu et al. (1994), Horst et al. (1996), Martin and Ward (1992)). One reason is that the research and development for sonochemical reactors are often based on empirical reaction data. Another one is the influence of many important effects, which are not understood thoroughly or in which various parts of different models are joined together to give the model for a sonochemical reactor. Our aim is to develop a model that is capable of predicting at least approximately the yield of sonochemical reactions in reactors of different geometric shapes and applications. One will then be in the position to develop ultrasonic reactors for different reactions with optimized characteristics, like maximum yield of reactions, etc. As a first step in this direction pressure fields in homogeneous media in reactors of different geometric shapes with various dimensions were calculated (Keil and Da¨hnke (1996, 1997a,b)). In this work a new model is proposed to calculate the three-dimensional pressure field in a liquid with an inhomogeneous distribution of cavitation bubbles. The acoustic field will be modeled stepwise in the beam direction, and with every step the bubble density distribution will be calculated and the respective parameters for phase velocity and attenuation will be considered. This model can be used for sonochemical reactors with any type of boundary conditions and will be further developed to reduce the approximations. 2. Theory Acoustic field theory is concerned with the investigation of time-dependent pressure fluctuations around the

static pressure in a compressible fluid, such as air or water. To generate a sound field, one must do mechanical work on the medium, e.g., by setting a liquid into motion by the vibration of a submerged plate. Thus, sound represents the transmission of mechanical energy and relies on the elastic and inertial properties of the medium through which it travels. To model the sound field in a fluid, the majority of past effort (Morris (1995), Partridge and Smith (1995), Soenarko (1993), Tjotta and Tjotta (1980), Banaugh and Goldsmith (1962) and others) was based on the solution of idealized problems, e.g., the scattering of acoustic waves on infinite rigid cylinders or ideal spheres. Most investigations are concerned with analytical solutions of the problems by expanding the wave equation into suitable infinite series. This work is based on the numerical solution of the Helmholtz integral equation and the Kirchhoff integral (Ziomek (1995), Morse and Ingard (1968), Junger and Feit (1986), Skudrzyk (1968, 1971), Born and Wolf (1975)). 2.1. The Three-Dimensional Sound Field in an Unbounded Medium of the Homogeneous Phase. A vibrating surface of finite dimensions generates a pressure field p(r,t) that is a function of four coordinates. Therefore, our starting point is the three-dimensional homogeneous wave equation (John (1982))

∆p(r,t) -

2 1 ∂ p(r,t) )0 cliq2 ∂t2

(1)

cliq describes the sound velocity in the liquid medium. For steady-state conditions (with the assumption of a harmonic time dependence of the pressure p(r,t) ) p(r) e-iωt), this can be written in terms of the wavenumber k in the form of a three-dimensional Helmholtz equation

∆p(r) +

ω2 p(r) ) (∆ + k2)p(r) ) 0 2 cliq

(2)

This equation will now be solved for a pulsating sphere, a simple, physically realizable source. In spherical coordinates the Laplace operator, ∆, takes the form

∆)

∂2 1 2 ∂ 1 ∂2 1 ∂ ∂2 + + + + 2 cot θ 2 2 2 2 2 r ∂r ∂θ ∂r r sin θ ∂φ r ∂θ r (3)

In the case of a pulsating sphere the pressure field is of spherical symmetry and the Laplace operator depends only on the first two terms. Therefore, we get the spherical wave equation

(

)

∂2 2 ∂ + k2 p(r) ) 0 + 2 r ∂r ∂r

(4)

Deriving a solution of eq 4, the boundary condition on the radiating surface has to be added. The radial velocity U(t) of the sphere’s surface may be described by

U(t) ) U0e-iωt

(5)

The sphere’s center is located at the origin r ) 0 and has the equilibrium radius a. With the first-order equation of motion in a fluid (Morse and Ingard (1968))

850 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998

∂U ) -∇p(r,t) ∂t

F

(6)

and its application to a radial motion, the boundary condition for the gradient of the pressure at r ) a gives

|

∂p(r,t) ∂r

r)a

) iFωU0eiωt

(7)

where F describes the density of the medium. The general solution of this equation is a superposition of inward and outward traveling waves whose amplitude is decreasing with increasing radius. If the only boundaries for p(r,t) are r ) a and r ) ∞, converging waves do not arise. Therefore, one gets the solution for p(r,t) in terms of the velocity amplitude of the source (Junger and Feit (1986)).

a2 ei(k(r+a)-ωt) p(r,t) ) -iFωU0 1 - ika r

(8)

Introducing the fluid flow q ) 4πU0a2, one obtains

q ei(k(r+a)-ωt) iFω 4π 1 - ika r

p(r,t) ) -

(9)

If the oscillating sphere is not located at the center of the coordinate but at r′, one has to replace r by |r - r′|. The pressure field of a point source can easily be derived with the condition of a vanishing radius, a, equal to zero, and one obtains

p(r,t) ) -

iFω eik|r-r′| -iωt q e 4π |r - r′|

(10)

Calculating the pressure field of a two- or threedimensional distributed source, one uses the idea that each elemental area of the transducer acts as an acoustic point source, emitting spherical waves which propagate outward with an amplitude that decays as r-1. In our case the source is a flat plate (Figure 1) which is attached to an infinite rigid plane (xy) and oscillates harmonically in a piston-like mode in a direction normal to the xy-plane. The infinite plane reflects acoustic emissions radiated by the transducer such that the entire acoustic energy is projected into the half-space in front of the source. Therefore, the fluid flow, q, becomes q ) 2U0 dS′ where U0 is the velocity amplitude of the infinitesimal area dS′. Due to the supposed linear behavior of the pressure waves, the total pressure of the elementary point sources results from a superposition of the pressure fields generated by all single-point sources:

p(r,t) ) -

(∑

iFω 2π

U0(r′)

r′

eik|r-r′| |r - r′|

)

dS′(r′) e-iωt

(11)

ik|r-r′|

∫∫U0(r′) |re - r′| dS′ e-iωt

iFω 2π

a plane transducer in an infinite rigid baffle, one has to calculate two acoustic fields. The first one is the common superposition of elementary point sources, and the second one is the self-reflected field of the transducer’s surface and the reflected field of the surrounding surfaces. 2.2. The Helmholtz Integral Equation: Embedding Boundary Conditions. Before deriving the formulas for the pressure field which arises due to the boundaries based on given values of pressure, the term under the double integral of eq 12 has to be discussed. The fraction g(r,r′) ) exp(-ik|r - r′|)/|r - r′| is a freespace Green’s function and a solution of the special inhomogeneous Helmholtz Equation

∆g(r,r′) + k2g(r,r′) ) -δ(r-r′)

(12)

S′

If the geometry of the system is not the combination of

(13)

where δ(r-r′) ) δ(x-x′) δ(y-y′) δ(z-z′) is the delta function in three dimensions. The delta function represents an omnidirectional point source at r′ with a unit amplitude. The free-space Green’s function is appropriate for outgoing radiation in an unbounded medium. The expression for the sound field in a closed volume will be determined by the contributions of the sources and the boundary terms, which represent the field reflected at the boundaries or enter through the boundaries from outside. Therefore, calculating the acoustic field in homogeneous fluid media in the presence of boundaries, the three-dimensional, inhomogeneous Helmholtz equation has to be considered (John (1982)):

(∆ + k2)p(r,t) ) x(r,t)

The most exact expression for the acoustic field can be found by increasing the number of point sources to infinity and eq 11 transforms into (Leighton (1994), Junger and Feit (1986), Archer-Hall and Gee (1980), Skudrzyk (1971), Zemanek (1971), Stepanishen (1970)).

p(r,t) ) -

Figure 1. Geometry for the calculation of the pressure field at a point r radiated from a plane piston transducer which is located at the xy-plane.

(14)

x(r,t) represents the input acoustic signal to the fluid medium due to the primary sources (transducer) and to the waves which emerge from the boundaries. Green’s function, G(r,r′), which is capable of satisfying these supplementary conditions, includes the free-space Green’s function g(r,r′)

G(r,r′) ) g(r,r′) + f(r′) f(r′) vanishes if no boundaries exist.

(15)

Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 851

Multiplying eq 14 with G(r,r′) and eq 13 with p(r,t) and subtracting the first from the second equation gives

G(r,r′) ∆p(r,t) - p(r,t) ∆G(r,r′) ) x(r,t) G(r,r′) - p(r,t) δ(r-r′) (16) To proceed further, the position vectors r and r′ are interchanged. With the property

G(r′,r) ) G(r,r′)

and

δ(r′-r) ) δ(r-r′) (17)

eq 16 can be rewritten as

G(r,r′) ∆′p(r′,t) - p(r′,t) ∆′G(r,r′) ) x(r′,t) G(r,r′) - p(r′,t) δ(r-r′) (18) where ∆′ is the Laplacian which works on the coordinate r′. This step enables us to perfom an integration over the volume V′ including the concerning pressure values p(r′,t). Integrating both sides over V′ leads to

∫∫∫∇′[G(r,r′) ∇′p(r′,t) - p(r′,t) ∇′G(r,r′)] dV′ ) V′

∫∫∫[x(r′,t) G(r,r′) - p(r′,t) δ(r-r′)] dV′

Figure 2. Subdivision of the reactor’s boundary into small planes.

(19)

V′

with the transformation

∇′[G(r,r′) ∇′p(r′,t) - p(r′,t) ∇′G(r,r′)] ) G(r,r′) ∆′p(r′,t) - p(r′,t) ∆′G(r,r′) (20) By using of the divergence theorem of Gauss, eq 19 can be rewritten to yield the Helmholtz integral equation

p(r,t) )

∫∫∫x(r′,t) G(r,r′) dV′ + V′

IInS′(r′) [G(r,r′) ∇′p(r′,t) S′

p(r′,t) ∇′G(r,r′)] dS′ (21) nS′(r′) is the unit normal vector of the vibrating body’s or boundary’s surface. Equation 21 enables one to calculate the acoustic pressure field within and on the surface forming the boundary of the medium. It is composed of a volume integral of the source function x(r′,t) over the bounded medium and a surface integral of the boundary values of p(r,t) and its outward-pointed normal gradient over the surface of the boundary. The difference between g(r,r′) and G(r,r′) lies in the fact that the first one describes the field of a point source in an unbounded medium whereas the second one describes the field of a point source in the presence of boundaries. Equation 21 is Fredholm’s integral equation of the second kind and has to be solved with suitable methods (Pipkin (1991)). In this work we apply approximations which enable us to simplify the calculation of the Helmholtz integral equation. The pressure field will be divided into the primary field, which results from the vibration of the transducer’s surface (see the volume integral in eq 21) into the unbounded medium and the reflected fields emerging from the boundaries. In the present paper only firstorder reflections are taken into consideration, as an inclusion of second-order reflections revealed only minor changes in the pressure field. In principle, the algorithm can take into account any number of reflections. To calculate the primary field, we use the free-space Green’s function g(r,r′) in eq 21. It will be calculated for the entire reactor volume dependent on the transducer’s position in space. In this case, p(r,t) consists

solely of outward-going waves and the second integral in eq 21 vanishes (Morse and Ingard (1968)). The calculation is equivalent to the term described in eq 12. The boundaries of the reactor, the bottom, the water surface, and the outer curved wall are approximated by small rectangular planes (outer wall) or cuttings of a circle (bottom and water surface). The principle of the boundary’s subdivision is illustrated in Figure 2. Due to the approximation that only first-order reflections occur at the reactor’s boundary, we again chose the free-space Green’s function to calculate the reflected field. Therefore, each part of the boundary seems to reflect the primary field into the unbounded medium. The calculation of the reflected field is done with the aid of the Kirchoff integral (Skudrzyk (1968))

p(r,t) ) -

[

1 eik|r-r′| nS′(r′) p(r′,t) ∇′ 4π S′ |r - r′|

II

]

eik|r-r′| ∇′p(r′,t) dS′ (22) |r - r′| which expresses the solution of the homogeneous wave equation in terms of the values of the solution and its first derivatives at all points on an arbitrary closed surface surrounding r. Each boundary, the bottom, the outer wall, and the water surface, is treated as a sole boundary, which radiates the pressure field into an unbounded medium. This condition will be realized by increasing the distance of the remaining elements of the respective closed surface to infinity, where the field vanishes. Therefore, Kirchhoff’s integral only has to be integrated over the surface of the respective boundary where the pressure values and their derivatives with respect to the unit normal vectors at the surface are determined by using the spatial boundary conditions (Born and Wolf (1975)). The resultant field arises from a superposition of the primary and the reflected pressure fields with their individual phases at each point. The Kirchhoff integral thus effectively replaces the radiating boundary by a distribution of point sources and of forces weighted according to the surface acceleration and surface pressure. This procedure is only feasible with the assumption of a linear wave propagation in a homogeneous medium. To calculate the field of the second-order reflection, the pressure values

852 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998

determined by the field of the first-order reflections are used as the boundary conditions used in the Kirchhoff integral for calculating the field emerging from the boundaries. Thus, if there are no boundaries at all and the source of the pressure field consists of a vibrating surface in a rigid baffle, the Kirchhoff integral will be reduced to calculate the part of self-reflection at this surface. By considering only a few point sources, the Kirchhoff integral will vanish as there will be no surface from which reflected waves can originate. In this work, the approach to calculate the acoustic pressure field in a reactor with an arbitrary shape is as follows: First, the pressure field generated in the reactor fluid by the transducer is modeled and calculated numerically with the aid of eq 12 in the range of the reaction vessel. The pressure values at the boundaries are included. As the pressure values of the primary field pP(r′) at the reactor walls are known by this sum, one can calculate the reflected pressure field pR(r) by taking the following boundary conditions into consideration: (a) For a rigid wall a doubling of the pressure values due to reflection at the wall is supposed (pR(r′) ) pP(r′)). Furthermore, the derivative of the total field pT(r′) ) pP(r′) + pR(r′) at the boundary normal to the rigid wall must be zero (n′(r′) [∇′ pT(r′)] ) 0). This simplifies the calculation of the Kirchhoff integral in eq 22 because the second term will vanish. (b) The free surface of the fluid nearly fulfills the condition of an ideal “pressure-release” boundary. Therefore, the total pressure pT(r′) must vanish at this surface (pR(r′) ) -pP(r′)) and the first term of the eq 22 will vanish. Calculating the reflected field generated at this surface, one needs to know the normal derivative with respect to the unit normal vector of the total field pT(r′) at the surface, as requested by the second term on the right-hand side of eq 22. This value cannot be calculated until the reflected field is known. Therefore, one calculates the pressure field reflected at a rigid wall with the aid of the first term of eq 22 and then converts it by a phase shift of π (reversed sign of the reflected pressure field). With this procedure the approximated pressure field of an image source which vibrates π out of phase is calculated. This method is consistent with the theory of “The Pressure-Release Boundary” as outlined by Junger and Feit (1986). One result is a vanishing total pressure field at the liquid’s surface (pT(r′) ) 0). The reflected pressure field arising from a boundary with any acoustic reflection parameters (transmission and reflection coefficients, phase shifts) can be modeled by appropriate mixing of the boundary conditions for an ideal rigid and pressure-release boundary. In this paper only first-order reflections are used to determine the entire acoustic field. This procedure is chosen in order to reduce the calculation time and because of the small values of the reflected pressure fields pR(r) at their corresponding opposite boundary sites compared to the primary field pP(r) at this location. In principle, higher order reflections can be employed in a straightforward manner. 2.3. Step-by-Step Development of the Sound Field. The Helmholtz integral equation is applicable only to wave propagation in homogeneous fluids. In the present work, the important parameters for modeling the acoustic field (phase velocity, damping constants) undergo heavy spatial changes which create important

Figure 3. Image source in acoustics and electrodynamics.

changes in wave propagation. Considering these spatial changes, the pressure field will be calculated step by step in the direction of the sound beam. Therefore, a new procedure is developed, which is based on the theory of image sources in acoustics (Morse and Ingard (1968), Junger and Feit (1986)) and electrodynamics (Jackson (1962)). Consider a point source in front of an infinite rigid plane with a distance d from that plane. Due to the influence of the plane, the resulting acoustic field of the point source and the infinite plane can be described by two sources with distance 2d from each other (Figure 3). The second source, called the image source, has its position on the opposite (imaginary) side of the rigid plane. Therefore, the main effect of the infinite rigid plane is to reverse the direction of the propagating wave normal to that plane. By using the Helmholtz integral equation (eq 21), one is able to calculate approximately the reflected field of the infinite rigid wall, if the pressure field on the wall’s surface is known. Due to the condition of an infinite rigid wall, the appropriate Green’s function G(r,r′) in eq 21 can easily be determined. In the case of a source in contact with the half-space, it will be twice the term of the free-space Green’s function g(r,r′) (Skudrzyk (1968)). On the wall, the normal derivative of the total pressure is zero. The remaining parts of the closed boundary required in the surface integral of the Helmholtz integral equation will vanish with the condition of radiation into an unbounded medium (Born and Wolf (1975)). The integral equation will be only an approximation because one is not able to calculate the pressure field arising from an infinite rigid wall due to the relevant numerical infinite calculation time for the integration over an infinite area. One has to limit the integration area to those sections, where significant pressure values for wave propagation can be found. It is found that in the present model a scaleup of the reactor dimensions by a factor of 2 is sufficient to simulate the condition of an infinite reflecting rigid wall. Consider a beam propagating in the positive zdirection in an unbounded medium (Figure 4). If one inserts an infinite rigid wall perpendicular to the beam

Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 853

Figure 4. Geometry of the reflection of an arbitrary ultrasonic beam at an infinite rigid wall.

Figure 5. Geometry of the step-by-step development of the threedimensional, time-independent pressure field.

direction at z0, determines the reflected field by the method of image sources (Figure 4a), and converts the direction of the reflected field again (Figure 4b), nothing will change. This seems to be a complicated method to determine the acoustic field of a transducer, but for the present model this an important element. Because one is able to calculate approximately the reflected field of this wall with the aid of the Helmholtz integral equation, the pressure field in the beam direction is calculated gradually by inserting an approximated infinite rigid wall for every step in the beam direction zi (Figure 5). Therefore, the step-by-step calculation of the time-independent, three-dimensional pressure field will be done by using the boundary conditions for an infinite rigid wall for every step in the beam direction with the aid of the second term on the right-hand side of eq 21. Another possibility is the use of the Kirchhoff integral (eq 22) by creating a new boundary for every step in the beam direction. The boundary then consists of an infinite plane orthogonal to the beam direction which emits the acoustic waves into the unbounded medium. In this case, the exact boundary conditions for the pressure and its derivative normal to the plane can be used. In our work, we calculate the pressure field step by step with the Helmholtz integral equation, because this integral needs less calculation time because of the smaller number of integrands. Both options were examined and led to the same results.

With this procedure one is able to use various damping constants and phase velocities for every step in the beam direction. This is a first approximation to calculate three-dimensional wave propagation in inhomogeneous media. For further applications it is intended to take into consideration different damping constants and phase velocities for every differential volume element of the reactor. 2.4. Wave Propagation in a Two-Phase Medium. Mixtures of liquids and small gas bubbles occur in many industrial processes. For the case of sonochemistry, much work has been done on formulating equations of motions for a bubbly liquid (Nakoryakov et al. (1993), Prosperetti et al. (1988); Prosperetti and Commander (1989), Gumerov (1994), Sangani (1991), Blake et al. (1994)). In this work, the authors will mainly concentrate on sound propagation in such a two-phase system. Gas bubbles in a liquid affect profoundly the acoustic behavior of the liquid. For example, the speed of sound in the medium depends roughly on the product of the compressibility and density of the medium (Morse and Ingard (1968)). Since the compressibility is primarily determined by the density of gas bubbles in the reactor medium, the resulting sound velocity in the mixture (down to 20 ms-1) is quite different from that of a pure liquid (in pure water approximately 1500 ms-1) or the pure gas (340 ms-1 for air). A number of investigators have studied theoretically and experimentally the changes of phase velocity and damping constants due to the presence of gas bubbles in a liquid (Ye and Ding (1995), Nakoryakov et al. (1993), Prosperetti and Commander (1989), Temkin (1988), among others). The current paper employs the theory of Prosperetti and Commander (1989) for a certain distribution of bubble radii. This theory is based on the model of van Wijngaarden-Papanicolaou (van Wijngaarden (1968)). The two essential equations of the model needed in the theory of Prosperetti and Commander are the continuity and momentum equations belonging to bubbly liquids. The continuity equation of the model by van Wijngaarden takes the form

∂β(t) 1 ∂P(r,t) + ∇u(r,t) ) 2 ∂t ∂t Fcliq

(23)

where F and cliq are the density and velocity of sound of the liquid and P(r,t) and u(r,t) the time- and spatialdependent pressure and velocity of the medium, respectively. The local fraction of volume occupied by the gas, β, is, in the case of the occurrence of bubbles with the same equilibrium radius R0, given by

β)

4π 3N R 3 V

(24)

R is the present radius of the bubbles, and N/V is their density. For a mixture containing bubbles of different sizes, the volume fraction β is determined by

β(r,t) )

∫0∞R3(r,R0,t) f(r,R0) dR0

4π 3

(25)

where R(r,R0,t) defines the instantaneous bubble radius at time t at position r having an equilibrium radius R0 and f(r,R0) dR0 defines the number of bubbles per unit volume with equilibrium radius R0 between R0 and R0 + dR0. The second essential equation of the van Wijngaarden-Papanicolaou model is the momentum

854 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998

equation belonging to bubbly liquids and reads

F

∂u + ∇P(r,t) ) 0 ∂t

(26)

Combining eqs 23 and 26 with the aid of equations of bubble dynamics, thermodynamics, and a linear approximation for ∂β/∂t (Prosperetti and Commander (1989)) yields a three-dimensional wave equation for the propagation of pressure waves in a liquid containing gas bubbles:

∆P(r,t) -

2 1 ∂ P(r,t) ) cliq2 ∂t2

Figure 6. Different types of modeled sonochemical reactors.

2

∂ R(r,R0,t)

∫0∞

-4πF

∂t2

R02f(r,R0) dR0 (27)

Assuming a sinusoidal bubble motion, Prosperetti and Commander (1989) find that

[

R(r,R0,t) ) R0 1 -

Q(r,t) FR02

ω02

]

1 - ω2 + 2ibω

(28)

where ω0 describes the resonance frequency, Q(r,t) is defined by Q(r,t) ) P(r,t) - p∞, where p∞ is the equilibrium pressure in the liquid, and b is a damping constant that arises from viscous, thermal, and acoustic effects. It exhibits a strong dependence on ω. The changes of b near the resonance frequency ω0 are significant. Substitution of eq 28 into eq 27 leads to

∆Q(r,t) +

ω2 Q(r,t) ) c2 R f(r,R )

∫0∞ω 2 -0ω2 + 02ibω dR0

-4πω2Q(r,t)

(29)

0

This equation serves as the final result to calculate the complex wavenumber, km, of the liquid medium containing gas bubbles of different densities and radii, defined by

∆Q(r,t) + km2(r) Q(r,t) ) 0

(30)

Therefore, the expression for the complex wavenumber, km, reads (Prosperetti and Commander (1989))

km(r) )

ω2 + 4πω2 2 cliq

R f(r,R )

∫0∞ω 2 -0ω2 + 02ibω dR0

(31)

0

In the present work a new wavenumber, km, will be calculated for every step in the beam direction depending on the bubble density determined as well. 3. Types of Sonochemical Reactors Modeled Three different types of ultrasonic reactors (Figure 6) were analyzed and their pressure fields calculated as a function of the density distribution of cavitation bubbles. The first reactor has a simple geometry: the ultrasonic transducer is a circular plate with a diameter of 0.06 m which performs a harmonic movement in the z-direction. Acoustic waves are generated and emerge from the bottom of the reaction vessel. The outer shape of the reactor is a circular tube with a height of 0.15 m.

Reactor 2 is a circular tube with open ends in both z-directions. This means that it may represent the active (insonated) zone of a sonochemical loop reactor. The acoustic waves are transmitted by three ultrasonic horns that are placed in an equilateral triangle in the same plane (z ) 0) around the tube. The vibrating cross section of each horn has the same shape as the transducer of the first reactor. The centers of the three transducers’ surfaces have their positions at the coordinates (φ ) 0, z ) 0) and (φ ) (2π/3, z ) 0) with the same distance r from the center as the corresponding reactor walls. Reflected waves originating from a surface orthogonal to the z-direction are neglected. For this reactor the pressure values were calculated for two planes: One is the rz-plane in cylindrical coordinates as used for the other reactor types. Pressure values are calculated up to a height of 0.15 m. The other plane is orthogonal to the rz-plane that meets the transducer surfaces in the center. The axes are described as the x- and y-axes like the Cartesian coordinates of a system. The third reactor that was modeled consisted of a cylindrical tube with the same dimensions as reactor 1. An ultrasonic horn emits the acoustic waves downward to the bottom from a height of 0.1125 m. All the calculations were carried out for transducer frequencies of 25 and 50 kHz and a uniform surface amplitude of 10-5 m. The reaction vessels described have unitary diameters of 0.12 m. For all types of reactors the calculations were made for the pure liquid and for a liquid with a homogeneous density distribution of cavitation bubbles. A procedure to calculate the pressure field step by step in the beam direction in order to calculate the time-independent acoustic field in a medium with an inhomogeneous density distribution is used only for reactor types 1 and 3. The method of step-by-step calculation has to be modified for reactor type 2 and is not yet implemented. 4. Details of the Program Implementation The analytical calculation of eqs 21 and 22 is possible only in a few special cases and with appropriate simplifying assumptions (Harris (1981), Seybert et al. (1986), Tjotta and Tjotta (1980)). As the modeling of pressure fields of different ultrasonic reactors is strived for and in most practical cases the shapes of the transducers’ surfaces or the boundaries are not symmetric, a numerical approach to the solution of eqs 21 and 22 was chosen. The surfaces of the transducers and the boundaries were subdivided into small regular pieces (Figure 2) with a minimum subdivision of 300 × 300 elements. The

Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 855

Figure 7. Pressure field of a circular transducer without the influence of boundaries for (a) 25 and (b) 50 kHz.

integrals have been calculated by a simple Riemann integration where the algorithms were described in FORTRAN 90 in connection with the Message Passing Interface (Gropp et al. (1994)). The Message Passing Interface is a library for parallel programming and makes it possible to reduce the calculation time effectively due to embedding several CPU’s. With this interface every CPU gets a certain amount of coordinate points for calculating the pressure field. For our application the following CPU’s were used: one Pentium Pro 200 MHz, three Double Pentium Pro 200 MHz, one Double Pentium Pro 180 MHz, and three Pentium 133 MHz. Concerning the model of reactor type 1, with a quantity of 100 × 200 pressure values and a surface subdivision into 300 × 300 pieces the computing time needed for the parallel algorithm takes 90 min. By using only one 200 MHz Pentium Pro CPU and a sequential FORTRAN 90 source code, the computing time takes at least 15 h. The calculation time for an IBM Workstation Model RS/6000 Model 3ct (which was used earlier) is nearly the same as that for a Pentium Pro 200 MHz CPU. Considerably longer calculation times were needed for reactor type 2 which does not have a cylindrical symmetry. One reason is that the pressure values at all boundaries must be calculated in addition to the primary field of the transducer that is only represented for one cross-section plane. If one wants to calculate the acoustic field of sonochemical reactors with a cylindrical shape and a symmetrical pressure field, the determination of the primary field for one cross-section plane, the plane that has to be considered, is sufficient to get the pressure values at all boundaries. For example, due to the angle-dependent pressure field in the xy-direction of reactor type 2, one has to calculate additionally the pressure field at the wall for one-third of a circle. Another reason is that the primary field of three transducers must be calculated whereas the other two types are modeled with only a single transducer. Furthermore, two independent acoustic fields for two planes (xy-plane and xz-plane) have to be modeled. For this type, the calculation time needed comes to 11 h. 5. Results In this section some results of the calculations are presented for the three different types of sonochemical reactors of different dimensions. All the figures are a representation of the amplitude of the pressure field |p(r,t)| in the rz-plane or xy-plane of the reactors.

5.1. Pressure Fields in a Fluid with a Homogeneous Phase. In Figure 7a,b the undisturbed pressure fields of a circular transducer that vibrates at a frequency of 25 kHz (Figure 7a) and of 50 kHz (Figure 7b) are represented. To make the influence of boundaries clear, Figure 8 shows the linear pressure field of all three types of reactors for a transducer frequency of 25 kHz. Figure 8a shows the field of reactor type 1 and is similar to Figure 8b, which describes the field of reactor type 3. The pressure amplitude in the region occupied by the ultrasonic horn (Figure 8b) is not defined, and therefore these values are set to zero. One recognizes in Figure 8b that the amplitudes for z greater than 0.1125 m are small compared to the remaining area. This is due to a shielding effect of the transducer’s body against the reflected waves from the boundaries. The pressure field for the maximum zvalue is zero due to the water surface which acts as an ideal pressure release boundary. Figure 8c is a representation of the pressure field of reactor type 2 in the xy-plane. In this figure a remarkable increase of the amplitudes can be observed. This is due to the presence of three transducers and the converging of the pressure waves to the center of the system. In Figure 8c the acoustic field of the same reactor for the plane perpendicular to the xy-plane is shown. The center of the transducers is located at the origin of the coordinates. The boundaries are located at |r| ) 0.06 m, whereas r ) 0.06 m and z ) 0 describe the center of one transducer surface and r ) -0.06 m and z ) 0 a boundary point at the center of the other two transducers. Therefore, the resultant pressure values are not symmmetrically distributed around the z-axis. Due to the concentration of acoustic energy at the center of the system, a pressure maximum at the origin can be observed. In Figure 9 the pressure fields of the same reactor types for a transducer frequency of 50 kHz are presented. Due to the same surface amplitude of 10-5 m as for the 25 kHz simulation, the amplitudes of all types show increased values compared to Figure 8. Another effect is an increased number of maxima and minima due to the reduced wavelength by a factor of 2. 5.2. Calculated Values for Phase Velocity and Damping Coefficients for the Refined Bubble Densities and Radii Distributions. To include the effects of wave propagation in a bubbly liquid, first of all one has to commit oneself to a certain radial distribution and volume fraction of cavitation bubbles.

856 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998

Figure 8. Pressure field for different types of reactors for a transducer frequency of 25 kHz: (a) reactor type 1, (b) reactor type 3, and (c and d) reactor type 2.

Prosperetti and Commander (1989)) assumed a Gaussian distribution

{

-(R-R0) /σ f(R0) ) ce 0 2

2

R1 < R < R2 otherwise

(32)

with a range of bubble radii from R1 ) 5 × 10-6 m to R2 ) 3 × 10-3 m and a standard deviation σ of about 2 × 10-3 m. The parameter c must be chosen to match the gas volume void fraction, β, which ranges from 10-5 to 10-2. In the calculations of Prosperetti and Commander the parameter R0 usually had a value of (R1 + R2)/2. Ye and Ding (1995) and Yount et al. (1984) took somewhat different distributions. They supposed a decreasing number of bubbles with the radius R:

f(R0) ) c(1/R4)

(33)

The range of radii and the volume fraction were limited similar to Prosperetti and Commander. In the present work a Gaussian radii distribution like eq 32 is assumed with a maximum value at 5 × 10-4 m as outlined in Figure 10. The ranges R1 and R2 span the entire area of previously observed bubbles and are the same as mentioned above. The maximum R0 is chosen to be in the vicinity of the lower limit due to the theories of Ye and Ding. Four different volume fractions of gas (β1 ) 10-5, β2 ) 10-3, β3 ) 10-2, and β4 ) 2 × 10-1) were used, and their influence on the acoustic field was analyzed. In Figure 11a the results for the attenuation coefficient R calculated from eq 31 are presented. It is the

imaginary part of the complex wavenumber, km, and determines the exponential decay of the three-dimensional propagating acoustic waves (with the relation p(x) ) p0e-Rx). The curves represent the results for the four different volume fractions βi as a function of the transducer frequency f in the range from 102 to 106 Hz. Considering the reactor’s size in the beam direction in the range of 10-1 m, zero amplitudes just at the transducer surfaces can be expected for β ) β4 for the transducer frequencies of 25 and 50 kHz and a remarkable decrease for β ) β3 at a frequency of 25 kHz. In the latter case the decrease of amplitude due to attenuation will be approximated by the ratio of pressure with attenuation to the pressure without attenuation:

|p(z)0.15m;β)β3)| |p(z)0.15m;β)0)|

≈ e-1.2×10 (1/m) 0.15m ≈ 0.17 1

(34)

This means that the amplitude of a single point source will be reduced to 17% compared to its value in the undamped case. Figure 11b shows the phase velocity cmix calculated from eq 31 for the four different volume fractions βi depending on the frequency f. The phase velocity for the highest volume fraction β4 increases in the range of the used frequencies by a factor of approximately 4. This means that the wavelength will be increased by the same factor compared to its original value. Therefore, a remarkable decrease of pressure maxima can be expected. This holds to a lesser extent for the volume fraction β3 and a frequency f of 25 kHz.

Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 857

Figure 9. Pressure field for different types of reactors for a transducer frequency of 50 kHz: (a) reactor type 1, (b) reactor type 3, (c and d) reactor type 2.

Figure 10. Assumed radial distribution of cavitation bubbles.

5.3. Pressure Fields in a Liquid with a Homogeneous Distribution of Cavitation Bubbles. For the modeling of the acoustic fields in liquids with a homogeneous distribution of gas bubbles, four different bubble volume fractions, βi, are assumed and the related attenuation coefficents, Ri, and phase velocities, cmix,i, calculated. The parameters are included in eq 21 or eq 22, and the pressure fields are calculated as in the case of a single-phase fluid. Solutions are obtained for all types of reactors and for transducer frequencies of 25 and 50 kHz. It is found that the acoustic field for bubble

volume fractions β of 10-5 and 10-3 remains almost unchanged. A very minor decrease of pressure amplitudes is observable for β2 ) 10-3, whereas the number of pressure maxima remains constant. For β3 ) 10-2 the decrease of pressure amplitudes is significant. The results of the calculations for this gas bubble volume fraction with a transducer frequency of 25 kHz are presented in Figure 12a-d. In Figure 12a the pressure field for reactor type 1 is shown. The acoustic field vanishes in the upper third of the reaction vessel. For reactor type 3 (Figure 12b) the acoustic wave amplitude decreases nearly to zero before it reaches the bottom of the vessel. The only reactor type with observable pressure amplitudes in the entire volume of the reaction vessel is the second one due to the concentration effects of the three transducers. The maxima of the amplitudes are located in the vicinity of the transducer surfaces (Figure 12c). Only a flat maximum can be obtained near the center of the system, as can be seen in Figure 12d. The results for a transducer frequency of 50 kHz and the same bubble volume fraction of β3 ) 10-2 are shown in Figure 13. The damping effects of the acoustic waves are similar to the calculated ones with a stimulation frequency of 25 kHz. Different are the pressure amplitudes at the top of the reaction vessel of reactor type 1 and the bottom of the vessel of reactor type 3. This is due to the lower attenuation coefficient calculated by eq 31 and presented in Figure 11a. A remarkable difference can be seen by comparing Figure 12c with Figure 13c. In the case of 50 kHz excitation the acoustic

858 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998

Figure 11. Attenuation coefficient R (a) and phase velocity cmix (b) depending on the transducer frequency f: s, β1; - - -, β2; -‚-, β3; -..-, β4.

Figure 12. Pressure field for different types of reactors for a transducer frequency of 25 kHz and an assumed homogeneous gas bubble void fraction of β3 ) 10-2: (a) reactor type 1, (b) reactor type 3, (c and d) reactor type 2.

field shows a homogeneous amplitude distribution in the xy-plane (this means that no point in this plane can be assigned to show the highest amplitudes), whereas in Figure 12c the amplitudes dominate in the vicinitiy of the transducer surfaces. Regarding the results presented in Figure 13d, the pressure fields reveal a resemblance to those in the undamped case in a homogeneous phase (Figure 9) with the exception of the pressure maximum in the center of the reaction vessel. Regarding both stimulation frequencies, an increasing

number of pressure maxima due to an increase of the wave velocity cmix (due to the cavitation bubbles) cannot be observed. Analyzing eq 31 and comparing it with Figure 11b, this effect will be observable for frequencies in the range of 10 kHz, where the wave velocity can attain the 3-fold value compared to the undamped propagation velocity. For volume fractions larger than β4 ) 2 × 10-1 the acoustic field amplitude is nearly zero in the vicinity of the transducer surfaces. In Figure 12a the acoustic field

Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 859

Figure 13. Pressure field for different types of reactors for a transducer frequency of 50 kHz and an assumed homogeneous gas bubble void fraction of β3 ) 10-2: (a) reactor type 1, (b) reactor type 3, (c and d) reactor type 2.

Figure 14. Pressure field for two types of reactors for a transducer frequency of 50 kHz and an assumed homogeneous gas bubble void fraction of β4 ) 2 × 10-1: (a) reactor type 1, (b) reactor type 2.

of reactor type 1 for a stimulation frequency of 50 kHz is presented and shows a small maximum just at the beam source, whereas in other areas the field vanishes totally. Even for reactor type 2 (Figure 14b), which produces the highest amplitudes, no pressure values can be observed except in the vicinity of the transducer surfaces. 5.4. Pressure Fields with Inhomogeneous Density Distributions of Cavitation Bubbles. The calculation of pressure fields in a medium with an inhomogeneous distribution of cavitation bubbles is

based on the procedure introduced in section 2.3 supported by the theory outlined in section 2.4 (Prosperetti and Commander (1989)). In this paragraph we investigate the influence of different bubble volume fractions on a three-dimensional wave propagation in an enclosed volume. This shows the importance of regarding different acoustic properties in modeling the pressure field in an ultrasonic stimulated fluid and gets one a little closer to the real conditions in the very complex phenomena occurring in a sonochemcial reactor. The range that the beam has to propagate until it

860 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998

Figure 15. Pressure field of reactor type 1 for a transducer frequency of 25 kHz, an assumed cavitation threshold of 106 Pa, and different inhomogeneous gas bubble void fractions (a) β ) β2 ) 10-3, (b) β ) β3 ) 10-2, and (c) β ) β4 ) 2 × 10-1.

reaches the opposite boundary was subdivided into 100 planes with equal distances from each other. The pressure field in each plane is calculated by using eq 21. For z ) 0 the first integral of eq 21 is employed referring to the movements of the transducer’s surface. For z > 0 the second term on the right side of eq 21 is employed referring to the boundary conditions existing on the previous plane as outlined in section 2.3. After calculating the acoustic field the regions which show pressure amplitudes greater than the cavitation threshold ptre (Young (1989)) are counted and the average bubble volume fraction βave was approximated by averaging the assumed bubble volume fraction βi over the relevant plane. Therefore, a homogeneous average bubble volume fraction for the relevant plane is assumed, which changes with every plane in the beam direction. For planes in the vicinity of the beam source the averaging of the volume fraction βi is carried out not in the entire area of the plane but in the area of the beam cone. This procedure is chosen to take only those parts of the plane into consideration which are essential for the beam propagation. The areas left out are not essential for the geometric structure of the beam and could be dropped for an approximate calculation of the propagating beam. This method is chosen to avoid the determination of bubble volume fractions βave in the vicintiy of the beam source which are small due to the nearly zero pressure amplitudes for small z values and radii r larger than 0.03 m (see Figure 7a,b). In addition to the contribution of the fields reflected at the boundaries, a main bubble volume fraction βmain is calculated

as an average value of all approximated volume fractions βave used for the calculation of the propagating beam. The averaging is also done over the volume of the beam cone and used for the whole volume of the reaction vessel. Based on the variety of experimentally determined or theoretically predicted cavitation thresholds ptre (Young (1989), Trevena (1984), Overton et al. (1984), Lauterborn (1970)), two different thresholds and three different bubble volume fractions β2, β3, and β4 are used to calculate the pressure fields in the different types of ultrasonic reactors. The first one is assumed to be ptre,1 ) 106 Pa which is an average value indicated by Trevena (1984), for ordinary (ptre ) 9 × 105 Pa) and boiled tap water (ptre ) 11.5 × 105 Pa). The second one is chosen to be ptre,2 ) 5 × 105 Pa, which is one of the parameters indicated by Young (1989), for a transient cavitation threshold. Figure 15 is a representation of the pressure fields of reactor type 1 for a transducer frequency of 25 kHz and a cavitation threshold of ptre,1 ) 106 Pa. In Figure 15a the assumed bubble volume fraction of β2 ) 10-3 Pa for pressure values greater than 106 Pa seems to have no observable effect on the amplitudes of the acoustic field in comparison with Figure 8a. By modeling the acoustic field with a volume fraction of β3 ) 10-2 (Figure 15b), a decrease of the pressure values with increasing distance from the source is observable which is, however, less than in the case of an assumed homogeneous distribution. Figure 15c shows the result for an assumed volume fraction of β4 ) 2 × 10-1. The damping of the

Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 861

Figure 16. Pressure field of reactor type 1 for a transducer frequency of 50 kHz, an assumed cavitation threshold of 106 Pa, and different inhomogeneous gas bubble void fractions (a) β ) β2 ) 10-3, (b) β ) β3 ) 10-2, and (c) β ) β4 ) 2 × 10-1.

acoustic wave is very strong, but pressure values different from zero can be observed even for distances far from the beam source. This is a considerable difference from the model of a homogeneous bubble distribution where the pressure field tends to zero immediately in the vicinity of the source. In Figure 16 the same results are shown for the case of a 50 kHz stimulation. Like the 25 kHz stimulation for volume fractions of β2 ) 10-3 (Figure 16b), no remarkable difference can be observed. The results for a void fraction of β3 ) 10-2 show a total decrease of nearly all amplitudes of about 20%. The highest pressure values are located at the beam source, as expected. Figure 16c reveals a remarkable difference between acoustic fields in media with homogeneous and inhomogeneous distributions of cavitation bubbles. Compared to Figure 14a, the acoustic field can be observed in the entire plane presented and the geometric structure of the pressure field is kept at the beam source. The results for the same cavitation threshold ptre,1 ) 106 and the bubble volume fractions β2, β3, and β4 for reactor type 3 are represented in Figure 17. The transducer frequency is 25 kHz. The pressure field is similar to that presented in Figure 16, whereas in the case of the highest volume fraction β4 (Figure 17c), the amplitudes at the bottom of the reaction vessel are larger than those shown in Figure 15c due to the lower distance of the boundary. In Figure 18 the results for both reactor types for the lower cavitation threshold of ptre,2 ) 5 × 105 Pa for a transducer frequency of 25 kHz

and an assumed bubble volume fraction of β3 ) 10-2 are represented. Compared to Figures 15b and 17b, no remarkable difference can be observed. This shows that the influence of the cavitation threshold ptre in this range, with the used frequencies and calculated primary amplitudes, is not a decisive parameter for calculating the wave propagation in a medium with an inhomogeneous distribution of cavitation bubbles. 6. Conclusions In this paper the theories of the Helmholtz integral equation (Morse and Ingard (1968)) and the Kirchhoff integral (Born and Wolf (1975)) were combined with the conditions of wave propagation in bubbly liquids (Prosperetti and Commander (1989)). Furthermore, a new procedure was introduced which enables one to calculate the time-independent, three-dimensional propagating wave step by step in the beam direction using different acoustic properties such as attenuation and wave velocity. This procedure makes it possible to model wave propagation in liquids with an inhomogeneous density distribution of cavitation bubbles. This is a common situation in sonochemical reactors. In a first approximation, the linear pressure fields in a pure liquid were calculated, with homogeneous and inhomogeneous density distributions of cavitation bubbles. Three different types of sonochemical reactors and two common stimulation frequencies were taken. Referring to homogeneous bubble distributions, it is found that the influence on the created acoustic field is

862 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998

Figure 17. Pressure field of reactor type 3 for a transducer frequency of 25 kHz, an assumed cavitation threshold of 106 Pa, and different inhomogeneous gas bubble void fractions (a) β ) β2 ) 10-3, (b) β ) β3 ) 10-2, and (c) β ) β4 ) 2 × 10-1.

Figure 18. Pressure field of (a) reactor type 1 and (b) reactor type 3 for a transducer frequency of 25 kHz, an assumed cavitation threshold of 5 × 105 Pa, and an assumed inhomogeneous gas bubble void fraction of β3 ) 10-2.

negligible for bubble volume fractions below β ) 10-3. Only for volume fractions greater than β ) 10-3 could remarkable changes be observed. The attenuation of the acoustic field in the bubbly liquid for bubble volume fractions of β ) 10-2 causes a considerable reduction of the pressure amplitudes. Regarding the reactor types 1 and 3 (Figure 6), the amplitudes vanish near the opposite boundaries to the beam source whereas for the second reactor type a considerable attenuation of the amplitudes could be observed. The pressure field of the latter one does not vanish due to the presence of three

beam sources and the concentrating of the acoustic energy to the center of the reaction vessel. For bubble volume fractions of β ) 2 × 10-1 the damping coefficient of the pressure waves increases to very high values, and the amplitudes of the acoustic field are zero just in the vicinity of the transducer surfaces. Investigating the acoustic field calculated in liquids with an inhomogeneous distribution of cavitation bubbles, remarkable differences to the fields modeled with homogeneous distributions could be observed for high bubble volume fractions such as β ) 10-2 and β ) 2 ×

Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 863

10-1. The damping of the pressure field decreases considerably with increasing the distance from the source. Acoustic field amplitudes could be observed in the entire area of the reaction vessel. In the case of a volume fraction of β ) 10-2 the structure of the field keeps its original form with a nearly homogeneous amplitude decrease of approximately 20%. For a volume fraction of β ) 2 × 10-1 the acoustic field keeps its original structure in the vicinity of the beam source, decreases to a very low value, and shows then an undamped propagation due to the calculated low bubble concentration more distant from the transducer. This is due to the the very high bubble volume fractions calculated in the vicinity of the beam source surface. Only minor effects could be observed by decreasing the value of the cavitation threshold ptre down to 50%. This parameter seems not to be important in this range for the investigated reactor types. Summarizing the results, one can conclude that the influence of inhomogeneous density distributions of cavitation bubbles compared to homogeneous distributions has remarkable effects on the propagation of the acoustic waves and, therefore, on the geometric structure of the ultrasonic field. An important common characteristic is that the number of active cavitation regions is not a simple function of the transducer frequency or surface amplitude: higher frequencies or amplitudes do not need to rase the number of active zones in the reactor (where sound-induced chemical reactions may happen). If the initial pressure field amplitudes reach certain values, the damping effects are dominating and the cavitation volume is restricted to the vicinity of the beam source. To raise the number of active zones to a maximum value depending on the geometry of the reaction vessel, the number and properties of the sound sources and the properties of the fluid medium have to be investigated in more detail. Therefore, our future efforts are concentrated on a refinement of the approximations used in the new proposed model of wave propagation in an inhomogeneous medium and on the use of extended theories for cavitation thresholds and wave parameters in bubbly liquids. The next step will be the inclusion of nonlinear effects. The results obtained so far reveal that the structure of pressure fields in sonochemical reactors has to be known in detail in order to be able to design such reactors. There is a need for standardized sonochemical reactors which enable comparison of yields and selectivities of sonochemical reactions. Acknowledgment The authors are grateful to the Max-Buchner-Forschungsstiftung and the Fonds der Chemischen Industrie for financial support. Notation a ) equilibrium radius of the pulsating sphere b ) damping constant arising from viscous, thermal, and acoustic effects due to the presence of gas bubbles cliq ) sound velocity in the liquid medium cmix ) sound velocity in the mixture of liquid and gas bubbles d ) distance of the point source in front of an infinite rigid plane f ) stimulation frequency of the transducer

f(r,R0) ) bubble density at r for bubbles having an equilibrium radius R0 g(r,r′) ) free-space Green’s function i ) imaginary unit k ) wave vector ()2π/λ) km ) wave vector in the mixture of liquid and gas bubbles n′ ) unit normal vector of the boundary S′ p(r,t) ) pressure value depending on the vector r at a time t ptre ) cavitation threshold pressure pP(r′) ) pressure values of the free space at the boundary at r′ pR(r′) ) pressure values of the field arising from the boundary at the boundary at r′ pT(r′) ) sum of pP(r′) and pR(r′) p∞ ) equilibrium pressure of the liquid far from the cavitation bubbles q ) fluid flow of the oscillating sphere or the point source r ) radial coordinate in the cylindrical or spherical coordinate system r ) coordinate vector r′ ) coordinate vector on a vibrating or reflecting surface t ) time coordinate u(r,t) ) velocity of the liquid medium depending on the vector r at a time t G(r,r′) ) Green’s function in consideration of boundaries N/V ) number of bubbles per unit volume P(r,t) ) pressure value depending on the vector r at a time t Q(r,t) ) difference between the pressure p(r,t) at the point r and time t and the equilibrium pressure p∞ R(r,R0,t) ) instantaneous bubble radius at time t at position r having an equilibrium radius R0 S′ ) area of the boundaries of the reaction vessel U(t) ) radial velocity of the oscillating sphere U0 ) radial velocity amplitude of the oscillating sphere V′ ) volume of the reaction vessel Greek Letters R ) attenuation coefficient due to the presence of cavitation bubbles in the liquid β ) bubble volume fraction βave ) average bubble volume fraction middled at the planes perpendicular to the beam direction βmain ) main bubble fraction middled over the whole volume of the beam cone φ ) axial angle in spherical coordinates θ ) azimuthal angle in spherical coordinates F ) density of the medium σ ) standard deviation of the Gauss distribution ω ) circle frequency ω0 ) resonance frequency of the cavitation bubble ∇ ) Nabla operator ∆ ) Laplace operator

Literature Cited Apfel, R. E. The Role of Impurities in Cavitation-Threshold Determination. J. Acoust. Soc. Am. 1970, 48, 1179-1186. Apfel, R. E. Acoustic cavitation inception. Ultrasonics 1984, 22, 167-173. Archer-Hall, J. A.; Gee, D. A single integral computer method for axisymmetric transducers with various boundary conditions. NDT Int. 1980, 13, 95-101. Banaugh, R. P.; Goldsmith, W. Diffraction of Steady Acoustic waves by surfaces of Arbitrary Shape. J. Acoust. Soc. Am. 1963, 35, 1590-1601. Blake, J. R., Boulton-Stone, J. M., Thomas, N. H., Eds. Bubble dynamics and Interface Phenomena; Kluwer Academic Publishers: London, 1994. Born, M.; Wolf, E. Principles of Optics; Pergamon Press: Oxford, U.K., 1975.

864 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 Brennen, C. E. Cavitation and bubble dynamics; Oxford University Press: New York, 1995. Gaitan, D. F.; Crum, L. A.; Church, C. C.; Roy, R. A. Sonoluminescence and bubble dynamics for a single, stable, cavitation bubble. J. Acoust. Soc. Am. 1992, 91, 3166-3183. Gropp, W.; Lusk, E.; Skjellum, A. Using MPIsPortable Parallel Programming with the Message-Passing- Interface; MIT Press: Cambridge, MA, 1994. Gumerov, N. A. Equations describing the propagation of nonlinear modulation waves in bubbly liquids. In Classical Electrodynamics; Jackson, J. D., Ed.; John Wiley & Sons: New York, 1962. Harris, G. R. Review of transient field theory for a baffled planar piston. J. Acoust. Soc. Am. 1981, 70, 10-20. Horst, C.; Chen, Y.-S.; Kunz, U.; Hoffmann, U. Design, Modeling and Performance of a Novel Sonochemical Reactor for Heterogeneous Reactions. Chem. Eng. Sci. 1996, 51, 1837-1846. Jackson, J. D. Classical Electrodynamics; John Wiley & Sons, Inc.: New York, 1962. John, F. Partial Differential Equations; Springer-Verlag: New York, 1982. Junger, M. C.; Feit, D. Sound, Structures, and Their Interaction; MIT Press: London, 1986. Keil, F.; Da¨hnke, S. Numerical calculation of pressure fields in sonochemical reactors (in german). Chem.-Ing.-Tech. 1996, 68, 419-422. Keil, F.; Da¨hnke, S. Numerical calculation of pressure fields in sonochemical reactorssLinear effects in homogeneous phase. Period. Polytech. Ser. Chem. Eng. 1997a, 41, 41-55. Keil, F.; Da¨hnke, S. Numerical calculation of scale-up effects of pressure fields in sonochemical reactorsshomogeneous phase. Hung. J. Ind. Chem. 1997b, 25, 71-80. Lauterborn, W. Contribution to a Theory of Cavitation Thresholds (Zu einer Theorie der Kavitationsschwellen). Acustica 1970, 22, 48-54. Leighton, T. G. The Acoustic Bubble; Academic Press: London, 1994. Luche, J.-L. Sonochemistrysfrom experiment to theoretical considerations. Adv. Sonochem. 1993, 3, 85-124. Martin, P. D.; Ward, L. D. Reactor Design for sonochemical engineering. Trans. Inst. Chem. Eng. 1992, 70, 296-303. Mason, T. J.; Lorimer, J. P. Sonochemistry: Theory, Applications and Uses of Ultrasound in Chemistry; Ellis Horwood Ltd.: New York, 1988. Miller, D. Experimental investigation of the response of gas-filled micropores to ultrasound. J. Acoust. Soc. Am. 1982, 71, 471476. Miller, D. L.; Nyborg, W. L. Theoretical investigation of gas-filled micropores and cavitation nuclei to ultrasound. J. Acoust. Soc. Am. 1983, 73, 1537-1544. Morris, P. J. The scattering of sound from a spatially distributed axisymmetric cylindrical source by a circular cylinder. J. Acoust. Soc. Am. 1995, 97, 2651-2656. Morse, P. M.; Ingard, K. U. Theoretical Acoustics; McGraw-Hill: New York, 1968. Naidu, D. V. P.; Rajan, R.; Kumar, R.; Gandhi, K. S.; Arakeri, V. H.; Chandrasekaran, S. Modelling of a batch sonochemical reactor. Chem. Eng. Sci. 1994, 49, 877-888. Nakoryakov, V. E.; Pokusaev, B. G.; Shreiber, I. R. Wave Propagation in Gas-Liquid Media; CRC Press: London, 1993.

Overton, G. D. N.; Williams, P. R.; Trevena, D. H. The influence of cavitation history and entrained gas on liquid tensile strength. J. Phys. D 1984, 17, 979-987. Partridge, C.; Smith, E. R. Acoustic scattering from bodies: Range of validity of the deformed cylinder method. J. Acoust. Soc. Am. 1995, 97, 784-795. Pipkin, A. C. A Course on Integral Equations; Springer-Verlag: New York, 1991. Prosperetti, A.; Commander, K. W. Linear pressure waves in bubbly liquids: Comparison between theory and experiment. J. Acoust. Soc. Am. 1989, 85, 732-746. Prosperetti, A.; Crum, L. A.; Commander, K. W. Nonlinear bubble dynamics. J. Acoust. Soc. Am. 1988, 83, 502-514. Sangani, A. S. A pairwise interaction theory for determining the linear acoustic properties of dilute bubbly liquids. J. Fluid Mech. 1991, 232, 221-284. Seybert, A. F.; Soenarko, B.; Rizzo, F. J.; Shippy, D. J. A special integral equation for acoustic radiation and scattering for axisymmetric bodies and boundary conditions. J. Acoust. Soc. Am. 1986, 80, 1241-1247. Skudrzyk, E. Simple and Complex Vibratory Systems; The Pennsylvania State University Press: University Park, PA, 1968. Skudrzyk, E. The Foundations of Acoustics; Springer-Verlag: New York, 1971. Soenarko, B. A boundary element formulation for radiation of acoustic waves from axisymmetric bodies with arbitrary boundary conditions. J. Acoust. Soc. Am. 1993, 93, 631-639. Stepanishen, P. R. Transient Radiation from Pistons in an Infinite Planar Baffle. J. Acoust. Soc. Am. 1970, 49, 1629-1638. Temkin, S. Attenuation and dispersion of sound in bubbly fluids via the Kramers-Kronig relations. J. Fluid Mech. 1990, 211, 61-72. Tjotta, J. N.; Tjotta, S. An analytical model for the nearfield of a baffled piston transducer. J. Acoust. Soc. Am. 1980, 68, 334339. Trevena, D. H. Cavitation and the generation of tension in liquids. J. Phys. D 1984, 17, 2139-2164. van Wijngaarden, L. On Equations of motion for mixtures of liquid and gas bubbles. J. Fluid Mech. 1968, 33, 465-474. Ye, Z.; Ding L. Acoustic dispersion and attenuation relations in bubbly mixture. J. Acoust. Soc. Am. 1995, 98, 1629-1636. Young, F. R. Cavitation; McGraw-Hill: London, 1989. Yount, D. E. On the evolution, generation, and regeneration of gas cavitation nuclei. J. Acoust. Soc. Am. 1982, 71, 1473-1481. Yount, D. E.; Gillary, E. W.; Hoffman, D. C. A microscopic investigation of bubble formation nuclei. J. Acoust. Soc. Am. 1984, 76, 1511-1521. Zemanek, J. Beam Behavior within the Nearfield of a vibrating Piston. J. Acoust. Soc. Am. 1971, 49, 181-191. Ziomek, L. J. Fundamentals of Acoustic Field Theory and SpaceTime Signal Processing; CRC Press: London, 1995.

Received for review May 13, 1997 Revised manuscript received December 14, 1997 Accepted December 16, 1997 IE9703393