Numerical simulation of ligament-growth on a spinning wheel

International Journal of Multiphase Flow 77 (2015) 90–103

Contents lists available at ScienceDirect

International Journal of Multiphase Flow journal homepage: www.elsevier.com/locate/ijmultiphaseflow

Numerical simulation of ligament-growth on a spinning wheel Jure Mencinger a,∗, Benjamin Bizjan b, Brane Širok c a

Faculty of Mechanical Engineering, University of Ljubljana, Aškerˇceva 6, Ljubljana 1000, Slovenia ˇ Izoteh d.o.o., Brnˇciˇceva 15b, Ljubljana-Crnuˇ ce 1231, Slovenia c Faculty of Mechanical Engineering, University of Ljubljana, Aškerˇceva 6, Ljubljana 1000, Slovenia b

a r t i c l e

i n f o

Article history: Received 31 March 2015 Revised 26 July 2015 Accepted 1 August 2015 Available online 20 August 2015 Keywords: Ligament growth Volume of fluid Gerris flow solver Fiberization

a b s t r a c t Liquid ligaments can grow from perturbations in liquid film spread on a spinning wheel due to the centrifugal force acting on the film. Typically, the growth is strongly influenced by the surface tension on the evolving liquid–air interface. This phenomenon, frequently exploited in industry for the production of fibers, was investigated numerically and the Volume-of-Fluid (VOF) methodology was used to model the interface. The freely available Gerris code with adaptive mesh refinement (AMR) was used to achieve the fine resolution of the computational grid required at the evolving liquid–air interface. The results were compared with the experimental data captured by a high-speed camera. The influence of the process operating variables on the ligament growth is also presented.

Introduction A spinning wheel atomizer, where a stream of liquid flows onto the mantle surface of the wheel (Fig. 1), is extensively used in industry for the robust disintegration of highly viscous or nonhomogeneous liquids. The most common application of these kinds of rotary devices is in the production of mineral wool and other fibers (Širok et al., 2008), where liquid ligaments rise from the melt film on the wheel and solidify into fibers. Usually, 2–4 spinning wheels are used together in the so-called spinning machines or spinners, with the melt cascading between the wheels and fiberizing. Compared to other standard atomizing techniques (e.g. using spinning disc and cup atomizers), spinning wheels have several advantages. First of all, due to the large liquid film width and the possibility of using multiple spinning wheels in a single device, relatively high flow rates of liquids can be atomized or fiberized on such devices. Furthermore, rotation around the horizontal rather than vertical axis results in a less complicated set-up when liquid droplets or solidified fibers are transported by the air flow. Also, spinning wheels can be used for the fiberization of very high temperature melts (e.g. rock wool) – up to a temperature of about 1500 ºC – due to the simple external geometry and large surface to volume ratio which allows for efficient cooling. While spinning wheels have a broad range of existing and potential applications, the extent of research and modeling of these par∗ Corresponding author’s present address at: R&D Competence Center Cooling, Gorenje d.d., SI-3503 Velenje, Slovenia. Tel.: +386 40380064. E-mail address: [email protected], [email protected] (J. Mencinger).

http://dx.doi.org/10.1016/j.ijmultiphaseflow.2015.08.002 0301-9322/© 2015 Elsevier Ltd. All rights reserved.

© 2015 Elsevier Ltd. All rights reserved.

ticular types of rotary devices has been scarce and mostly limited to integral time and length scales. Only recently, Širok et al. (2014) performed a high-speed camera visualization of a melt film on a wheel of an industrial spinner, analyzing the velocity and structural dynamics of the early fiberization phase. An experimental study has also been conducted recently by Bizjan et al. (2014a,2014b) on a model spinning wheel with non-solidifying Newtonian liquids being used as working media. Ligament evolution was visualized on this kind of atomization device and analyzed to form multiple regression dynamic models for characteristic ligament properties, such as the ligament number, diameter, length, and strain rate. Some other types of rotary atomizers which, unlike the spinning wheels, incorporate a central liquid or melt influx, have been studied more extensively due to less complex boundary conditions, resulting in an easier experimental and mathematical modeling. Eisenklam (1964) formed a theoretical model for inviscid liquid ligament formation on spinning discs and cups. Kamiya (1972) improved the model to also account for viscosity effects and experimentally verified the ligament number equation on spinning discs. Eisenklam’s and Kamiya’s theory has been tested experimentally by many authors, recently by Liu et al. (2012a,2012b). The fundamental hydrodynamic mechanism of liquid disintegration on spinning wheels is similar to the mechanism on spinning discs and cups as pointed out by Bizjan et al. (2014a). Generally, ligament formation on all of these types of rotary devices is driven by a centrifugal force causing the Rayleigh–Taylor hydrodynamic instability. However, the kinematics and dynamics of ligament formation on spinning wheels are significantly different and more complex than their spinning disk or cup counterparts. Ligament formation on

J. Mencinger et al. / International Journal of Multiphase Flow 77 (2015) 90–103

Fig. 1. Schematic (side) view of spinning wheel atomizer (left) and ligament formation/disintegration (right). [Reprinted from Bizjan et al. (2014b) with permission from Elsevier.]

spinning wheels is asymmetric and unsteady, in addition to being random, as identified by Bizjan et al. (2015), who performed a nonlinear analysis of melt film velocity time series, obtained from an industrial mineral wool spinning machine. For these reasons, spinning wheel atomization and fiberization has so far mostly been modeled experimentally since analytical and numerical modeling is very complex. While experimental modeling of the spinning wheel ligament formation mechanism has given some valuable results, such an approach faces serious limitations. First of all, the possibilities for operating parameter variations are limited by time consuming experiments and data post-processing. Also, due to safety and cost issues, experimental modeling of actual mineral melt fiberization is impractical while experimentation in an industrial environment is dictated by the plant operator and subject to many additional sources of measurement uncertainty (while the cost and safety issues still remain). Furthermore, as fiberization is a solidifying multiphase flow, some of the micro-scale ligament properties such as the local temperature and pressure distribution, though of great importance for an understanding of the fiberization process, cannot be measured by any of the currently known measurement methods. For these reasons, numerical ligament formation modeling in a computational fluid dynamics (CFD) environment is necessary for further improvements in the modeling of the process. At this point it is necessary to mention that the numerical modeling approach is not new in the modeling of ligament dynamics and other properties. For example, Westerlund and Hoikka (1989) and Širok et al. (2008) numerically modeled cooling, kinematics and mechanical tension of ligaments. However, the ligaments were treated as cylindrical solids rigid in the longitudinal direction and ideally flexible in the transverse direction, which is not a realistic approach in the early fiberization phase when solidification has not yet occurred. On the other hand, some authors modeled ligament growth and breakup, but often without taking the heat transfer and solidification in consideration. For example, Olesen (1997) developed a numerical model for ligament breakup under the action of capillary instabilities. Tong and Wang (2007) simulated the evolution of initially still axisymmetrical ligaments in the air; they explained well the end-pinching mechanism responsible for the breakup of drops at each end of the ligament. A simulation of a jet breakup into ligaments and droplets – capturing an impressive level of details – was presented by Shinjo and Umemura (2010) who used massive computational resources with up to 5760 computational cores. The study closest to ligament fiberization was conducted by Purwanto et al. (2005) who formed a mathematical and numerical model for molten slag granulation using a spinning disk atomizer and it was also verified experimentally. In this work, the ligament growth phase is considered and simulated numerically, starting from the imposed embryonic stage.

91

Although the current study does not consider solidification or even heat transfer, it represents a step towards the simulations of fiberization. Already, the influence of the operating variables such as the wheel rotational speed f and film thickness h (connected with the liquid flow rate) on ligament growth can be studied. In principle, the ligament disintegration phase could also be considered. However, as shown later in the paper, an extremely high resolution of the computational mesh is required to capture this phenomena correctly; such resolution requires large computational resources even with the used state-of-the-art adaptive mesh refinement algorithm. This paper is organized as follows. In section “Overview of the ligament-formation mechanism”, an overview of approaches to the melt fiberization mechanism is presented. The physical model consisting of the relevant equations and boundary conditions used in the present work is described in section “Physical model” along with the experimental procedure used for the verification of the simulation results. Section “Reasons to use Gerris code and its numerical solution setup” contains a brief presentation of the Gerris flow simulation code and the presentation of the used numerical solution setup. Numerical results are reported in section “Results” and conclusions are given in section “Concluding remarks”. Overview of the ligament-formation mechanism Fiberization is a process in which a hot liquid (melt) is scattered into ligaments which then solidify into fibers. If the liquid does not solidify, it breaks up into droplets and the process is called atomization. However, both fiberization and atomization mechanisms are similar in the initial phase when the liquid ligaments form from the annular liquid film on the wheel mantle surface. The initial ligament formation process can be described as follows. A liquid stream of flow rate Q falls on the mantle surface of the spinning wheel with radius R. The liquid is drawn in motion by adhesion, viscous forces and surface tension (Fig. 1). A thin film of thickness h and width B is formed across the complete wheel circumference. Because of the rotation, the centrifugal force acts upon the film, pushing it radially into the surrounding air. The established unstable condition, due to the density of the liquid being much higher than the density of the air, presents Rayleigh–Taylor instability (Eisenklam, 1964). The Rayleigh–Taylor instability manifests itself in the development of the initial perturbations, which can originate from liquid feed fluctuations, from the wheel vibrations, and according to Westerlund and Hoikka (1989) also from the shear force between the melt film and the surrounding air (Kelvin–Helmholtz instability). The growth or the decay of the perturbations at the interface of two fluids was first explained by Taylor (1950) in terms of the linear stability analysis. Depending on the wavelength, the spatial Fourier modes, proportional to exp (ikx) (k ≡ 2π /λ), of initial perturbation develop with different exponential growth rates exp (t/τ k ) – each represented by time constant τ k . The fastest growing mode, i.e. the one with the smallest possible positive value of τ k , becomes predominant and, in the considered case, its wavelength λm transforms into ligament spacing s (Fig. 1). Using this approach, Eisenklam (1964) analyzed the instability of a liquid layer of thickness h, adhering to the lip of a rotating cup. He obtained the dispersion relation

1

τk2



= Rω − 2

k2 σ

ρ



k tanh (kh)

(1)

where R, ω, ρ , and σ denote radius at the rim of the cup, cup’s angular velocity, the density of the liquid1 , and the surface tension coefficient, respectively. While the perturbation growth is accelerated by the centrifugal force, it is damped by the surface tension, as can

1 The properties of the surrounding air are neglected due to its considerably lower density.

92

J. Mencinger et al. / International Journal of Multiphase Flow 77 (2015) 90–103



be understood from (1). The modes for which k > ρω2 R/σ = kc (or λ < λc ) are fully suppressed. The wavenumber of the fastest growing mode km is obtained by derivation of (1) with respect to k, which yields the equation



Rω2 − 3

k2 σ





ρ

tanh (kh) + Rω2 −

k2 σ



ρ

2

kh sech

(kh) = 0.

(2)

The above equation can be solved numerically: √ √ however, for . . kh  1 and kh  1 we get km = kc / 2 and km = kc / 3, respectively. As km ∝kc , it follows that

λm ∝

1

ω



σ . ρR

(3)

While Taylor (1950) and Eisenklam (1964) disregarded viscosity of fluids in their analyses, Fermigier et al. (1992) experimentally and theoretically studied the Rayleigh–Taylor instability of a thin layer of viscous liquid spread over a glass plate. They observed the development of two dimensional patterns on the liquid–air interface, subject to gravitational acceleration g, after the plate had been turned over. Their work also provides the dispersion relation, based on the lubrication hypothesis, which can be written as

1

τk

=

h3 (ρ gk2 − σ k4 ) 3η

(4)

where η denotes the√dynamic viscosity of the liquid. From (4), kc = √ ρ g/σ and km = kc / 2 can be easily calculated. To a certain extent, (4) could also be valid for the film on a wheel rotating with angular velocity ω by replacing g with R ω2 , provided (i) that R ω2  g and also (ii) that the impact of the principal curvature of the film 1/R can be neglected (when σ /R  ρ h R ω2 ). If so, then an identical expression for kc and similar values of km are obtained from (4) as from (1). However, it should be noted that (4) was obtained by assuming very long timescales (τm = 12σ η/(h3 ρ 2 R2 ω4 )  h2 ρ /η) and consequently with neglect of the inertial terms in the Navier–Stokes equation, which in the considered case can be hardly justified. Furthermore, the use of the lubrication hypothesis to obtain (4) requires that h2 /λ2m  1, which translates to the condition ρ Rω2 h2  8π 2 σ . The dispersion relation, unconditionally valid for the liquid film on the surface of a rotating wheel, has, to the best of our knowledge, not yet been obtained. An interesting parameter for practical use is N – the number of ligaments formed along the circumference of the rotary device. Provided that the ligaments are generated in a single plane with a constant axial coordinate, N can be estimated as N ∼ 2π R/s. Since the parameter s (or λm ) is in practice not easily accessible, several equations for the number of ligaments were derived mathematically or empirically. For example, based on the energy balance of the disturbance in the liquid film, Kamiya (1972) proposed an implicit relation valid for spinning discs



N

2

Fig. 2. Schematic diagram of the experimental set-up for visualization of the atomizing process. [Reprinted from Bizjan et al. (2014b) with permission from Elsevier.].



2

3 + (8N − 3)Oh



1+

1+



1 2

NOh

− We = 0.

(5)

In (5) We and Oh stand for non-dimensional Weber and Ohnesorge  number, respectively, defined as We ≡ ω2 R3 ρ /σ and Oh ≡ μ/ σ ρ R. The analysis leading to (5) is two-dimensional – this approach is enabled by the fact that ligaments on the spinning disc form in a single plane. This cannot be claimed for the spinning wheel where the axial coordinate of the ligaments is limited only by the width of the film. Based on the multiple regression analysis of the experimental data (Bizjan et al., 2014a) obtained the following relation for the number of ligaments:

N = 0.360 × We0.433 × q0.810



(6)

with q denoting the flow number (Širok et al., 2008) q ≡ Q/B ρ /Rσ . Also, the ligaments were found to be non-uniformly spaced along the

wheel circumference (Bizjan et al., 2014a) and to form randomly. One has to keep in mind that the axial coordinate of the ligament was not measured due to limitations of the experimental method. Nevertheless, the statistical coefficient of determination R2 for (6) was high (≈0.9), indicating a strong correlation between N, We and q. Typically, faster rotation (higher We) and higher flow rate (higher q) cause the number of ligaments to rise. At the same time, the mean ligament diameter dL (Fig. 1) is reduced. While growing, the ligaments can also be considered as jets, although they are not flowing at a constant flow rate through a fixed nozzle. Finally, the ligaments disintegrate into drops to minimize the surface energy. The formation of droplets is attributed to the capillary instability induced by the surface tension and damped by the viscous forces. This breakup mechanism is known as the Plateau– Rayleigh breakup and typically occurs when the effect of aerodynamic forces on the ligament surface stability is low (Lefebvre, 1989). This is the case with the spinning wheel atomization as the liquid flow is typically laminar (Westerlund and Hoikka, 1989). Mean diameter of droplets dD to which the ligaments disintegrate can be estimated by (Liu, 2000)

d D = dL

 6

4.5π 2 (1 + 3OhL )



(7)

with OhL = μ/ ρσ dL denoting ligament Ohnesorge number. To summarize, the lifetime of a chosen ligament can be roughly divided into three stages: (i) early, i.e. embryonic stage represented by a disturbance which will grow into the ligament, (ii) ligament growth, and (iii) ligament disintegration. In the production of fibers, however, the ligaments should solidify before the disintegration stage begins. As stated in Introduction, the present work considers mainly the ligament growth phase which is simulated numerically, starting from the imposed embryonic stage. Physical model Present numerical simulations are based on a physical model customized to the experimental study by Bizjan et al. (2014a, 2014b). The experimental setup consisted of a spinning wheel, i.e. a plastic cylinder with R = 45 mm radius and 50 mm length. The cylinder was spun by a DC-powered motor around its lateral axis which was oriented horizontally (Fig. 2). The wheel rotational speed f (f ≡ ω/2π ) was varied in the 5–25 Hz range with 5 Hz increments. The working medium – different glycerol–water mixtures – was gravity-fed from a supply tank and poured onto the wheel through a vertically oriented circular nozzle with dN = 3.0 mm diameter at two different flow rates (1.63 mL/s and 3.27 L/s). To limit the scope of the present paper, 3 we consider only the 85% glycerol–water mixture (ρ = 1224 kg/m , μ = 0.107 Pa s, σ = 0.0647 N/m) and use only the experimental data obtained at Q = 1.63 mL/s flow rate. Experimental data was obtained by means of high-speed camera visualization whereby the observed liquid structures were

J. Mencinger et al. / International Journal of Multiphase Flow 77 (2015) 90–103

93

Fig. 3. Experimental data – an image showing the growing ligaments on the bottom half of the cylinder. The circular line presents the (assessed) surface of the cylinder.

illuminated by a diffuse coaxial LED lighting system. Bizjan et al. (2014a,2014b) performed both qualitative and quantitative analyses of the acquired images, while in this paper, experimental images are used to provide the initial and boundary conditions for the numerical simulations as well as for the evaluation of the simulation results. The used experimental setup, however, does not provide the axial coordinate of the evolving ligaments nor the film width B. Had the wheel dried out in one rotation and the film thickness was constant in the axial direction (i.e. function of azimuthal coordinate only) then B could be obtained (assuming h  R) from Q ≈ 2π fRhB. As in practice the wheel does not dry in one rotation, we can assess B ࣡ Q/(2π fRh). In this work, we focus on the ligament growth from the bottom half of the wheel surface shown in Fig. 3. As the figure shows, both the ligaments of variable length as well as ligament embryos exist in this region. The chosen area is therefore suitable for a numerical investigation of the ligament growth. Moreover, by choosing the bottom part, the need for simulating the liquid flow through the nozzle is evaded, thus reducing the computational demands. As stated above, the simulations deal with the ligament growth from the embryonic phase; the causes for the initial appearance of embryos are not studied. The modeling equations Although we are interested primarily in the liquid phase forming of the ligaments, we obviously are dealing with two-phase flow, having a well-defined interface between liquid and gas (air). An established approach to deal with such two-phase flow is the use of the Volume-of-Fluid (VOF) model (Hirt and Nichols, 1981; Scardovelli and Zaleski, 1999). The model is based on the introduction of the volume fraction c(x, t)



c(x, t ) =

1, 0,

x in liquid x in gas

(8)

Fig. 4. Domain defined around the bottom half of the wheel.

with p(x, t) the pressure field, ρ (x, t) the fluid density, μ(x, t) the dynamic viscosity, and f(x, t) the body force field. Using the VOF model, the fluid properties in (10) are written as function of c(x, t)

ρ(c) = cρl + (1 − c)ρg , μ(c) = cμl + (1 − c)μg

(11)

where ρ l , ρ g , μl , and μg denote the density and viscosity of liquid and air, respectively. Introducing (11) in (10) makes the latter equation valid in both phases simultaneously. Additionally, the surface tension which plays an important role in the ligament formation process must be taken into account. In the VOF model, this is achieved by defining the body force field fσ as

ˆ fσ = σ κδs n

(12)

using the Dirac distribution δ s which is nonzero on the interface. ˆ its normal pointing In (12), κ is the curvature of the interface and n (in our case) from gas to liquid phase so that, for example, the curvature of the interface of a spherical liquid drop is positive and that of a spherical air bubble is negative. The total body force field f in (10) therefore equals

f = fσ + ρ g

(13)

with g denoting the gravitational acceleration. which implicitly defines the interface as the discontinuity in c. The volume fraction is passively transported by the fluid as described with the advection equation

∂c + u · ∇c = 0 ∂t

(9)

where u(x, t) ≡ (u, w, w) represents the fluid velocity field. The hydrodynamic description of the liquid film, developing ligaments, and surrounding air is given by the Navier–Stokes equation, which for incompressible fluid (∇ · u = 0) reads



ρ

 ∂u t + u · ∇ u = −∇ p + ∇ · (μ(∇ u + ∇ u )) + f ∂t

(10)

The computational domain The computational domain is defined by firstly aligning the z-axis with the horizontal axis of the rotation so that (hereafter, us ing also the cylindrical coordinates r ≡ x2 + y2 and θ ≡ arctan −y x where appropriate) the surface r = R coincides with the surface of the rotating cylinder. The domain extents are defined with the interW vals x ∈ [− W 2 , 2 ], y ∈ [−H, 0], and z ∈ [0, D] as shown in Fig. 4. We set H = W = 0.1 m based on the experimental evidence – the size of 2 the domain is sufficient to ensure it has the necessary space for the development of ligaments.

94

J. Mencinger et al. / International Journal of Multiphase Flow 77 (2015) 90–103

The size of the domain in z (axial) direction D is chosen by considering the following aspects: (i) D should be larger than the expected thickness of the ligaments (including the evolving droplets) – as no forces act (or they are neglected) in the axial direction the ligaments should not protrude out of the axial boundaries, (ii) smaller D decreases computational burden but it should be large enough to diminish the influence of insufficiently well-defined boundary conditions, (iii) although film width B is not known, accepting it to be wider than the domain enables us to circumvent the modeling of the film contact angle at the wheel surface, (iv) in the used Gerris code the domain is composed of cubical boxes – it is very convenient to have D = H/2m where m is integer number. In the present case D = H/26 = 1.5625 mm is chosen. The description of the boundary conditions for both (9) and (10) follows.

1, 0,

r ≤ Rh r > Rh

(14)

with Rh ≡ R + h. Eq. (14) describes the surface of the supposedly undisturbed film of height h which can thus be described with the equation r = Rh . Expressing the boundary velocity components ub and vb with the tangential and radial component uθ and ur (u = −uθ sin θ + ur cos θ , v = −uθ cos θ − ur sin θ ), respectively, we state

 ωr,

n uθ b = ωRh Rrh ,



c0 (x) =

r ≤ R∗h r > R∗h

1, 0,

with

R∗h ≡ R + h +



r ≤ Rh , r > Rh

urb = 0.

(15)

The above equation describes the rigid body rotation of the liquid film (r ≤ Rh ) and the decrease of the tangential velocity in the air (r > Rh ) proportionally with r−n . The discussion about the suitable value of the exponent n follows later in the paper. The approximation of the rigid body rotation of the film solves the Navier–Stokes equation for the hypothetical stationary state – when the pouring of the liquid is stopped and no perturbations are present – if the gravitational acceleration g is neglected. The ratio g/(ω2 R) equals 2.21 × 10−1 , 2.45 × 10−2 , and 8.84 × 10−3 at f = 5.0 Hz, f = 15.0 Hz, and f = 25.0 Hz, respectively. Thus, the approximation is better for higher rotational frequencies. By using (14) and (15) we accept that the width of the film is larger than twice the domain depth (B > 2D); this may even not be true in practice, but we assume that the ligament growth is not influenced by the axial coordinate as long as the ligament is sufficiently distant from the axial boundary of the film. For this reason – and also to reduce the CPU usage – the embryos are centered in the plane z = 0 assuming their corresponding symmetry. The boundary conditions at the boundary y = 0 for x > 0 representing (azimuthal) inflow are also defined with (14) and (15) – here, the inflow velocity is normal to the domain boundary. The outflow boundary conditions are assumed at the remaining domain boundaries: x = ± W 2 , y = −H, and y = 0 for x < 0. Finally, on the wheel surface, the no-slip boundary condition is specified for (10) and the Dirichlet boundary condition c = 1 for (9).

(16)

gi (θ , z).

(17)

i

In (17), the summation encompass all the considered embryos and the function gi (θ , z) determines the shape of the i-th embryo prescribed with

gi (θ , z)=

In order to virtually double the domain without increasing the computational burden, the symmetrical boundary conditions are imposed at the plane z = 0; all the derivatives of the considered fields in the normal direction to the plane equal zero. On the other side, at the plane z = D, Dirichlet boundary conditions are imposed by defining c(x, y, z = D) ≡ cb (x, y) and u(x, y, z = D) ≡ ub (x, y) ≡ (ub , vb , 0). The function cb (x, y) equals



Definitions (14) and (15) could also be used to impose the initial conditions u(x, t = 0) = u0 (x) and c(x, t = 0) = c0 (x) if the introduction of the ligament embryos were not necessary. At any rate, the ligaments’ growth is initiated by defining c0 (x) as

⎧  ⎪ ⎨

Boundary conditions

cb =

Initial conditions



Ai cos

⎪ ⎩

si  π R2h (θ −θi )2 +z2 z2 , (θ − θi )2 + 2 ≤

θi R h Rh

0,

1 2

θ i

2

otherwise. (18)

The constants Ai , θ i , and θ i define the height of the embryo (i.e. initial length of the growing ligament), its width at the base (foot), and the position of its maximum, respectively. The exponent si additionally adjusts the shape of the embryo. As stated before, the embryos are centered in the plane z = 0. The initial velocity field u0 (x) ≡ (u0 , v0 , 0) is, as above, more simply written with tangential and radial components uθ 0 and ur0



uθ 0 =

ωr,   n ∗ ωR∗h Rrh ,

r ≤ R∗h r>

R∗h



,

ur0 =

u˜ri , 0,

Rh ≤ r ≤ R∗h r < Rh ∨ r > R∗h (19)

where u˜ri denotes the pseudo velocity in the radial direction inside i-th embryo, intended to approximate the initial velocity of embryos. Namely, the velocity field prescribed by (19) is neither continuous across the boundaries of the embryos nor solenoidal. Although this disagrees with the presumed incompressibility of the fluid, it does not present a problem as the numerical solver of (10) makes the velocity field divergence free in the first time step by using the projection method (Popinet, 2003), i.e. by projecting the velocity field onto the space of solenoidal vector fields. Reasons to use Gerris code and its numerical solution setup Known accuracy issues with VOF simulations The system of Eqs. (9) and (10) can be solved numerically with many existing commercial and non-commercial finite volume (and perhaps finite elements) codes since the VOF model is wellestablished in the CFD community. However, when surface tension plays an important role in the behavior of the considered system, the accuracy of the solution of (9) and the calculation of the forces (12) becomes essential and thus the selection of the solution code is critical. Hereafter, we briefly describe the main issues concerning accuracy when using the VOF model. The solution of (9) principally presents implicit interface tracking: in the context of the finite volume discretization, the interface is found in the so called mixed cells (i.e. finite volumes) with the (average) value of the volume fraction c inside the interval [0, 1] (0 < c < 1). Conversely, the cells with the average value c = 0 or c = 1 contain only single phase: in our case, gas or liquid, respectively. Regardless of the simplicity of (9), its numerical solution requires special discretization schemes to prevent excessive numerical diffusion, which results in the “smearing of the interface”. Furthermore,

J. Mencinger et al. / International Journal of Multiphase Flow 77 (2015) 90–103

the occurrences of unphysical values of c: c < 0 or c > 0 should also be prevented. The way to calculate fluxes on faces between neighboring cells roughly divides the schemes to solve (9) in two groups: algebraic and geometric. In the former group, the face flux is obtained from the chosen (algebraic) relations containing the values of c in the neighboring cells. An example of a well-used algebraic scheme is CICSAM (Ubbink and Issa, 1999). Geometrical schemes, on the other hand, typically contain two steps: the interface reconstruction and the flux computation, by taking into account the geometrical information. The first step is commonly done with the Piecewise Linear Interface Construction (PLIC) (Youngs, 1982), i.e. by representing the interface in each mixed cell by a line/plane (in 2D/3D). The geometrical flux computation consists of (i) the determination of the so-called donating region (DR) – the area occupied by (both) fluid(s) crossing the considered cell-face in the actual time-step, and (ii) the calculation of the area inside the DR occupied by each fluid taking into account the previously constructed interface. Whereas algebraic schemes are simpler to implement, especially on general (i.e. non-orthogonal and/or unstructured) grids, geometric schemes prove to be less diffusive and therefore more accurate. As stated above, besides the accuracy of the solution of (9), another important aspect is the accuracy of the calculation of the surface forces – represented in the VOF model as body force (12). As shown by Brackbill et al. (1992), (12) can be consistently replaced by

fσ = σ κ∇ c

(20)

which is the commonly used Continuum Surface Force (CSF) form. Obviously, the accuracy of the calculated fσ depends on the accuracy of the calculated curvature κ and also ∇ c. The most accurate among existing methods to calculate κ appears to be the method using the so-called Height Functions (HF) (Cummins and Francois, 2005) and the Parabolic Reconstruction of Surface Tension (PROST) (Renardy and Renardy, 2002), both achieving second order accuracy for sufficiently resolved interfaces. The HF method is practically applicable only on Cartesian grids, whereas PROST requires multivariable minimization which can be rather time consuming and/or fragile. Key features of Gerris flow solver Considering the described accuracy issues, the Gerris finite volume flow solver – created as open source code by Popinet (2003,2009) – appears to be the most appropriate solver for the problem at hand. Gerris solves the time-dependent incompressible variable-density Navier–Stokes equations using second-order time and space discretization on Cartesian grids. The code contains a geometrical PLIC–VOF advection scheme for solving (9) with the direction splitting (thus obtaining rectangular donating regions). Local curvature of the interface is calculated with the generalized HF method. Gerris controls the time-step size t; it follows the numerical stability requirement (Popinet, 2009) that the time-step should be smaller than the period of the shortest capillary wave, i.e.



t ≤

ρ 3 πσ

(21)

with denoting the cell size. Additionally, classical CFL condition

t ≤ /|u| must be fulfilled to ensure meaningful use of the VOF advection algorithm. Which of the two criteria is more stringent depends on the considered case. One of the strongest features of Gerris is its adaptive mesh refinement (AMR) capability, based on the octree (quadtree in 2D) division of cells. The resolution is adapted automatically and dynamically to the features of the flow, based on the user prescribed criteria. This enables extreme grid refinement – up to the prescribed level

95

Fig. 5. Left: the simulation domain composed of 6 948 boxes presenting zeroth level grid; visible cylinder walls intersect the boxes. Right: grid in plane z = 0 at initial stage near four introduced embryos; the liquid–gas interface presented with PLIC–VOF planes, colored by the calculated curvature.

lmax – in the vicinity of the interface, thus allowing for its more accurate definition. To account for solid boundaries inside the domain (the cylinder walls in the present case), Gerris uses the so-called embedded boundary technique. The cells which are cut by the solid boundaries can be refined automatically using the octree structure. Surfaces can be defined implicitly; for example, the cylindrical surface of the wheel is described with equation x2 + y2 − R2 = 0. Gerris code runs in parallel, employing the MPI library and using a domain decomposition approach: the computational domain is partitioned utilizing chosen existing box boundaries. Due to AMR, the dynamic load-balancing is implemented by repartitioning the domain after the prescribed number of time-steps to achieve optimal use of the available computational resources. Simulation setup The computational domain is composed of interconnected cubic (square in 2D) boxes. The edges (surfaces) of each box therefore present either the domain boundary or a connection to the adjacent box. All the boxes are split evenly and contain 23l0 (22l0 in 2D) cells (each) with l0 representing the initial level of refinement. For the actual calculation presented here, the domain is composed of Nbox = 6 948 boxes of size D = 1.5625 mm, interconnected with 13 676 internal links as shown in Fig. 5 (left).2 This rather large number of boxes is used to obtain the previously described thin computational domain – in our case determined with box width D. A large number of boxes also facilitates load balancing – distributing the computational burden between the used CPU cores. For example, Fig. 6 shows the domain decomposition for the simulation using 24 cores; at first glance the burden seems to be unevenly balanced, however, one has to take into account that the density of cells varies significantly throughout the domain. Results Grid resolution test In order to estimate the accuracy of the calculations, the case of the growing ligament, defined by the parameters given in Table 1, was

2 The described structure of boxes was obtained by starting with a 2D domain (x · y) composed of two square boxes of size H = 0.1 m; the boxes were split 6 times by using 2D splitting tool which is accompanied with Gerris code; the obtained 2D setup was then transferred manually to 3D.

96

J. Mencinger et al. / International Journal of Multiphase Flow 77 (2015) 90–103

Fig. 6. Example of domain decomposition to 24 parts indicated by colors; grid (showing the z = 0 plane) is refined near the liquid–air interface. Table 1 Parameters defining the case used for the grid resolution test. f0

h

θ1

θ 1

A1

s1

u˜r1

15.0 Hz

0.45 mm

0.305

0.069

0.78 mm

2

0.42 m s−1

Table 2 Grid resolution ( min ), number of cells (Ncells ), average time-step size ( ¯ t), and passed wall clock time using 24 computing cores (twall ) at t = 0.015 s. Grid Coarse Medium Fine

lmax 4 5 6

min

Ncells −5

9.76 × 10 m 4.88 × 10−5 m 2.44 × 10−5 m

∼ 94, 000 ∼ 227, 000 ∼ 760, 000

t

twall −6

5.51 × 10 s 1.48 × 10−6 s 2.88 × 10−7 s

4.7 h 35.5 h 321.7 h

run using three different levels (lmax ) of refinement near the interface. The maximal grid resolution min = D/2lmax equals approximately 0.1 mm, 0.05 mm, and 0.025 mm on the tested coarse (lmax = 4), medium (lmax = 5), and fine (lmax = 6) grids, respectively. Table 2 reveals the number of cells (Ncells ) comprised in each of the used grids at approximately half time during the simulations (at t = 0.015 s). Ncells varies due to the AMR – it typically rises in time, which is consistent with the increasing interfacial area. Naturally, Ncells affects the total computational time in two ways: firstly, the solution of the discretized equations on larger grids is more time consuming at each time-step, and secondly, smaller time-step sizes are required on finer grids to satisfy stability criteria (CFL condition and (21)) thus resulting in an increased number of necessary time-steps. As Gerris code adapts time-step size during the calculation, only the average time-step size t is written in Table 2. The wall clock time (twall ) required to simulate ligament growth up to t = 0.015 s using 24 cores varied from ∼4.5 h on the coarse grid to ∼13.5 days on the fine grid. For the illustration of the effectiveness of AMR: if equivalent uniform grids were used they would contain Nbox × 23lmax cells which equals ∼ 28 million, ∼ 227 million, and ∼ 1820 million cells for the coarse, the medium and the fine grid, respectively. Such large grids would be prohibitively large for the actual time-dependent simulations. Fig. 7 shows the calculated interface positions on the z = 0 plane at t = 0.0(0.005)0.025 s, obtained on the three used grids: good qualitative agreement (i.e. the shape and inclination of the ligament) between the results on all three grids are obtained until t = 0.010 s. At t = 0.015 s, however, good qualitative agreement is only achieved between the results on the medium and the fine grid, as the shape of the ligament acquired on the coarse grid differs more notably from the shapes obtained on the two finer grids. The correspondence between the simulated ligament growth on the three grids further decreases with increasing t as the ligament becomes thinner.

Fig. 7. Grid resolution test for a growing ligament at f = 15 Hz with three degrees of near-interface refinement. Superimposed curves c = 0.5 indicate the interface position on z = 0 at t = 0.0(0.005)0.025 s.

Fig. 8. Pressure field and cell-boundaries at the ligament tip area on plane z = 0 at t = 0.015 s, obtained with three degrees of near-interface refinement.

The discussed results indicate that the simulations can be considered as trustworthy until the minimal ligament thickness dmin decreases down to the size of four cells. Namely, at t = 0.010 s,dmin  0.4 mm which equals approximately four times the minimal cell size on the coarse grid Cmin . The ligament thins further so that at t = 0.015 s dmin  0.3 mm, can be acquired from the medium and fine grid calculations. On the coarse grid, however, at the same instant dmin  0.2 mm  2 Cmin is obtained. The most apparent cause for the deviation starting to increase with insufficient grid resolution is indicated in Fig. 8 showing the pressure field (in the plane z = 0) near the ligament tip at t = 0.015 s calculated on the three used grids. The lateral variation of the pressure field inside the ligament on the coarse grid (Fig. 8, left) obviously does not correspond well with the transverse curvature of the ligament surface. This is not the case on the medium and the fine grid which seem sufficiently dense to capture important details at the considered time (t = 0.015 s). Obviously, when dmin < 4 min the pressure inside the ligament due to the surface tension cannot be reliably calculated. Although the ligaments can be presented even when dmin  2 min the results on such under-resolved grids become mostly unusable. As demonstrated in Fig. 9, such a transverse resolution causes a notable deterioration of the curvature calculation accuracy

J. Mencinger et al. / International Journal of Multiphase Flow 77 (2015) 90–103

97

Table 4 Parameters defining four embryos for the comparison with the experimental data.

Fig. 9. Close-up of the ligament and the pinched-off bubble representation on the coarse grid at t = 0.020 s. Left: ambiguity of the isoline c = 0.5; middle and right: presentation of the interface with planes (PLIC), colored by the curvature. Uncertainty of the calculated interface orientation and curvature clearly depends on the relative grid resolution.

i

θi

θ i

Ai [mm]

si

u˜ri [m/s]

1 2 3 4

0.737 0.652 0.491 0.305

0.084 0.072 0.081 0.152

0.596 1.421 1.091 0.882

4 4 4 6

0.72 0.72 0.72 0.72

ever, this is not unforeseen as Popinet (2009) noted that the curvature calculation method “only shows first-order convergence” at intermediate resolutions ( 12 d  4δmin ). A possible explanation for an even bigger discrepancy is not obvious: it could originate from missing the asymptotic grid convergence range to possible inaccuracies in the code implementation. Nevertheless, the grid convergence is obtained (albeit not at very high rate) and the average error of rmax on the medium grid during the ligament growth phase can be estimated to be around 1.0 mm. Due to the rather large computational cost on the fine grid, the medium size grid (lmax = 5) is used for all other simulations presented in this paper. Comparison with the experiment

Fig. 10. Identification of four embryos (e1, e2, e3, and e4) with indicated measurement uncertainty of their height. The experimental image is superimposed by the c = 0.5 isoline in the z = 0 plane of the initial field c0 (x).

(Fig. 9, middle). Also, the insufficient transverse resolution results is the uncertainty of the calculated interface (PLIC) orientation. As Fig. 9 (left) indicates, even the isoline c = 0.5 (although, not directly affecting the calculations) becomes ambiguous when the ligament width approaches 2 min . As the ligament grows further, it becomes thinner and approaches the stage of breaking into drops. However, even the used fine grid (lmax = 6) lacks the resolution necessary to correctly capture the breakup phase; a still finer grid (i.e. lmax > 6) on the current setup should be used. The simulations of this phenomenon with the present approach thus require larger computational resources than used in this work. Quantitatively, we assessed the error of the calculations by comparing the maximal reach rmax of the ligaments in the radial direction calculated on the three grids. Assuming sufficient convergence inside each time-step (i.e. neglecting the iterative convergence error), the error of the calculation of rmax on grid G can be estimated as the disG RE |. Here, r RE denotes Richardson’s − rmax cretization error DE G = |rmax max extrapolation (Roy, 2005)



RE F rmax = rmax +

F M rmax − rmax , p 2 −1

ln p=

M rCmax −rmax M −r F rmax max



(22)

ln (2)

F M , rmax and rCmax representing rmax obtained on the fine, with rmax medium and coarse grid, respectively. The discretization errors at three different instants are given in Table 3; as expected, the errors grow with time and decrease with the grid resolution. The observed order of accuracy p between 0.59 and 0.78 is notably lower than the theoretical value p = 2.0 for the second order discretization. How-

To test how well the simulations correspond with reality in view of the assumptions made, the results of the simulations are compared with the experimental data which is in the form of the recorded images of the evolving ligaments. Specifically, we chose the experimental data obtained at the nominal frequency f = 14.8 Hz. In the experiment, however, a certain slip (i.e. lower angular velocity) of the film surface was detected (Bizjan et al., 2014b). The approximation of the rigid body rotation of the film, imposed as boundary and initial conditions, does not account for slip – therefore a reduced rotational frequency f ∗ = 14.25 Hz is used in the simulations for comparison. This rotational frequency is obtained simply by following the positions of the ligament feet in the chosen sequence of images. The simultaneous growth of four selected ligaments is compared. This requires the identification of four embryos (Fig. 11, left) in the first (initial) image, defining also the instant t = 0.0 s. The selected embryos are described in terms of (17) and (19), i.e. by estimating h and the constants θ i , θ i , Ai , si , and u˜ri which are written in Table 4. The film thickness is estimated at h = 0.35 mm. All the estimates are acquired directly from the images – by counting pixels and comparing the counted number with the number of pixels in the diameter of wheel R. The left columns of Figs. 11 and 12 display the comparison of the simulation against the experimental results: the isocurves c = 0.5 which indicate the interface position (in the plane z = 0) are drawn on the corresponding images from the experiment. The calculated ligament shapes correlate well with the observed ones. To quantify the agreement of the results, the distances between the tips of the calculated and of the observed ligaments are printed on the images. They were measured directly from the images with the accuracy of ∼0.1 mm (1 pixel) and vary from 0 to 1.4 mm, generally increasing with time. If the distances, normalized with the assessed lengths of the corresponding ligaments, are used as a measure of discrepancy between the calculated and the observed results, then this can be

Table 3 Error estimation using Richardson’s extrapolation. t [ms]

rCmax

M rmax

F rmax

p

RE rmax

DEC

DEM

DEF

5.0 10.0 15.0

47.9 mm 54.0 mm 66.1 mm

48.2 mm 54.8 mm 67.3 mm

48.3 mm 55.3 mm 68.0 mm

0.585 0.780 0.778

48.7 mm 55.9 mm 69.0 mm

0.72 mm 1.89 mm 2.88 mm

0.48 mm 1.10 mm 1.68 mm

0.32 mm 0.64 mm 0.98 mm

98

J. Mencinger et al. / International Journal of Multiphase Flow 77 (2015) 90–103

Fig. 11. Left: experimental images superimposed by isoline c = 0.5 in plane z = 0 at t = 1.6(1.6)8.0 ms with indicated distances between the observed and the calculated ligament tips. Right: reconstructed interfaces colored by curvature.

estimated to be less than 5%. In accordance with our previous discussion, the calculated shapes start to deviate more notably from the observed ones at t = 16.0 ms – when the thickness of the ligaments #2 and #3 decreases below the size of four computational cells (∼0.2 mm on the used grid). The ligament thickness can be read from the right columns of Figs. 11 and 12 displaying the reconstructed interface, which is colored depending on the calculated curvature κ : the ligament thickness is – except at the ligament foot and tip – approximately equal to 2/κ . The influence of operating variables on ligament growth Ligament growth depends on the shape and initial velocity of the corresponding embryo and on the operating variables such as the wheel rotational speed f0 and the liquid flow rate Q. As noted in section “Physical model”, both of these affect the liquid film thickness h; the relationship between them can be roughly approximated as Q ≈ 2π hBRf0 supposing that (i) the film thickness does not vary significantly across film width B and (ii) B settles to a constant after a

certain distance from the flow impinging point on the wheel. Except for the assumption B > 2D, B and thus Q do not have a direct influence on the numerical simulations. Nevertheless, the increase of Q at constant f increases h, for which the impact on the ligament growth is studied. To assess the significance of h, its value is varied from 0.15 mm to 0.55 mm in steps of 0.1 mm. The growth of the ligament starts from embryo #4 (i = 4) defined in Table 4 at the wheel rotational frequency f = 15.0 Hz. Fig. 13 reveals that the length of the ligament – not only the distance of ligament tip from the wheel – increases with increasing h. This can be understood by observing Fig. 14 which presents the comparison of two ligaments at h = 0.15 mm and h = 0.55 mm. Clearly, when h is larger the ligament draws more surrounding liquid as indicated by the more concave liquid–air interface in the vicinity of the ligament foot. We also investigated the influence of the wheel rotational frequency on ligament growth. One should keep in mind that in practical application decreasing/increasing f requires also correspondingly decreasing/increasing Q in order to keep h constant. For the present

J. Mencinger et al. / International Journal of Multiphase Flow 77 (2015) 90–103

99

Fig. 12. Left: experimental images superimposed by isoline c = 0.5 in plane z = 0 at t = 9.6(1.6)16.0 ms with indicated distances between the observed and the calculated ligament tips. Right: reconstructed interfaces colored by curvature.

test h is set to 0.35 mm. The ligament growth is again initiated from embryo #4 defined in Table 4. Fig. 15 enables the comparison of the calculated shapes at different points in time, chosen in a way to match the positions of the ligament feet at different rotational frequencies, i.e. the ligaments are displayed at the same values of ft. As expected, since the centrifugal force increases proportionally to f2 , the ligaments grow faster with increasing f. At the lowest tested rotational frequency ( f = 5.0 Hz) the embryo flattens due to the surface tension; the centrifugal force is too weak for the ligament to grow from the considered embryo. Certainly, a critical frequency f∗ exists, under which the embryo flattens instead of growing further into a ligament. f∗ can be assessed by comparing the centrifugal force Fc , which tends to accelerate the embryo in a radial direction, with the force Fσ pulling the embryo in the opposite direction due to the surface tension. A simple estimate can be obtained by approximating the embryo shape as a hemisphere with radius ξ (ξ  R); taking into account also h  R, the centrifugal force . Fc can be expressed as Fc = 23π ξ 3 ρl (2π f )2 R and the opposing force

as Fσ ࣊2π ξ σ . Equating Fc and Fσ yields the estimate of f∗

f∗ 

1



2π ξ

3σ Rρl

(23)

which is inversely proportional with the radius ξ . As in this work the shape (height) of the embryo is approximated with function gi (θ , z), defined in (18), ξ is not directly available. Instead, the equivalent radius ξ i is used; it is obtained by equating the volume of the embryo Vi ≡ ∫gi dS with the volume of a hemisphere 23π ξi3 . Vi can be simply calculated by integrating in local cylindrical coor π /2 dinates which gives Vi = π2 (Rh θi )2 Ai Ii with Ii ≡ 0 x( cos x)si dx. Hence, ξi =

 3

3 A (R θi )2 Ii which gives π2 i h

ξ i ∼ 1.24 mm for the con-

sidered embryo #4; introducing this value in (23) yields f ∗ = 7.59 Hz. This estimate corresponds well with the results of the simulations as 5.0 Hz < f < 10.0 Hz. However, it can be expected that the true values of f∗ are higher than those that are calculated from (23) since the viscous forces were neglected in obtaining the present estimate.

100

J. Mencinger et al. / International Journal of Multiphase Flow 77 (2015) 90–103

Fig. 13. Influence of film thickness h on ligament growth: h = 0.15(0.1)0.55 mm. Note: the calculated ligament for h = 0.15 mm at t = 20.0 ms is drawn although its thickness approaches the computational limiting width of two grid cells.

Fig. 14. Ligament close-up at t = 10.0 ms with curvature scale clipped at the interval [−1000 m−1 , 50 m−1 ]. Left: h = 0.15 mm , right: h = 0.55 mm.

The higher tested values of f elicit ligament growth from the imposed embryo. The relative importance of liquid’s inertia versus viscous and capillary forces, expressed through the ligament Reynolds number ReL = ρl ωRdL /ηl and ligament Weber number WeL = ρl ω2 R2 dL /σ , respectively, increases with f so the shape of the ligament approaches what we can name the kinematic shape. The latter shape is acquired by connecting mutually independent imaginary (liquid) particles repeatedly “released” (in a tangential direction) from the chosen point fixed on the surface of the rotating wheel. Supposing that the liquid particles detach from the wheel at the instants t∗ in time interval [t0 , t] with t0 denoting the starting time of ligament growth at azimuthal coordinate θ 0 , the kinematic shape can be depicted by a curve Kg (x, y, t) written in parametric form (t∗ ∈ [t0 , t]) as





x(t, t ∗ ) = R cos (θ ∗ ) − (t − t ∗ ) ω sin (θ ∗ )





y(t, t ∗ ) = −R sin (θ ∗ ) + (t − t ∗ ) ω cos (θ ∗ ) −

1 g(t − t ∗ )2 2

(24)

where θ ∗ ≡ (t ∗ − t0 ) ω + θ0 denotes the starting azimuthal coordinate of the particle released at t∗ . Without the gravity term

1 ∗ 2 2 g(t − t ) , (24) becomes a parametric equation of an involute of a circle with radius R rotating with angular velocity ω – let us denote this curve by K0 (x, y, t). Both, Kg (x, y, t) and K0 (x, y, t) are considered and also drawn in Fig. 15. While K0 (x, y, θ /ω) is invariant, since the particles’ velocity magnitude equals ωR, the curve Kg (x, y, t) approaches K0 (x, y, t) with increasing f because the importance of gravitational acceleration g decreases with g/(ω2 R). By comparing the actual and kinematic shape of the ligament, the importance of surface, viscous and aerodynamic (i.e. drag) forces can be observed. In the earlier phase of growth, the surface tension tends to straighten the ligament – this is more pronounced at lower values of f (e.g. 10 Hz). As Fig. 15 demonstrates, the length of the ligament (at the same foot position) increases with increasing f and, despite being shorter, the ligament’s shape (longitudinal curvature) closely matches Kg (x, y, t). In the later phase of ligament growth (from ≈ 6 o’clock foot position in the considered case), the shape starts to deviate from Kg (x, y, t) while the ligament tip mostly remains on it, albeit not at its end. The ligament tip (and later the pinched-off drop) moves in an almost straight path – its steepness increases with f approaching the maximum possible steepness defined by the line connecting the

J. Mencinger et al. / International Journal of Multiphase Flow 77 (2015) 90–103

101

Fig. 15. Influence of rotational frequency on film growth: f = 5.0(5.0)30.0 Hz. Gray curves indicate the kinematic ligament shape: thick – without gravity, dashed – with gravity.

ends of Kg (x, y, t). One also notes that the foot of the ligament starts to lag behind the foot of Kg (x, y, t) with the lag increasing in time and with f. Both the deviation of longitudinal curvature and the delay of the foot position, appearing at the later stage of ligament growth, can be attributed to the aerodynamic drag as discussed in the following section. Note on the impact of the velocity field in surrounding air Principally, we could impose boundary conditions for quiescent air at the boundaries of the computational domain to avoid the need of prescribing the velocity profile; however, this would require a much larger domain, impracticable for the present computations. Instead, the velocity field in the air at the domain boundaries is set by (15) stating that the decrease of tangential velocity uθ b in the

air (r ≤ Rh ) is proportional to r−n . If the cylinder is infinite then n = 1 (Landau and Lifshitz, 1987). Because the cylinder is finite and the surrounding air in the axial direction from the cylinder is quiescent, we can expect a more rapid decrease of tangential velocity with r (i.e. n > 1). Fig. 16 shows the influence of n by comparing the growth of the ligaments (from embryo #4 in Table 4, h = 0.35 mm, f = 15.0 Hz) using the values n = 1, n = 3, and n = 5 in (15). The influence of differently set uθ b is practically unnoticeable until the ligament reaches a sufficient length – at t = 12.5 ms for the considered case. With further growth of the ligament, its declination becomes more notably influenced by uθ b . At any rate, n = 3 is used in all previously presented simulations – this value appears appropriate since setting n = 3 results in a good correspondence between the calculated and the observed ligament inclination (Fig. 12).

102

J. Mencinger et al. / International Journal of Multiphase Flow 77 (2015) 90–103

Fig. 16. Influence of the boundary imposed air velocity profile uθ b by (15) with n = 1, n = 3, and n = 5 on ligament growth.

Even though the boundary conditions set through (15) with different value of n present better or worse approximation to reality, Fig. 16 offers a possible explanation of the ligament’s forward (longitudinal) curvature near the tip. As explained in the following paragraph, this shape can be attributed to the aerodynamic drag force FD on the ligaments. The drag is approximately squarely proportional to the difference

u in the velocities of the ligament and of the surrounding air; in the considered case u consists mostly of the azimuthal component

u ࣃ uθ . Fig. 16 indicates that uθ – which increases with n – has a stronger effect on the thin part of the ligament than on the evolving drop at the ligament tip. This can be rationalized as follows: FD on an object moving in gas with relative velocity u is commonly written as FD = 12 ρg ( u)2 CD A where CD represents the drag coefficient and A the cross-sectional area. Approximating the evolving drop with a sphere and the segment of the ligament with a cylinder, their deceleration aD due to the drag equals

⎧ 2 ⎪ ⎨ ρg ( u) d ρl aD = 2 ρ g ( u) ⎪ ⎩ d ρl

3 sph C 4 D 2 cyl CD

π

for drop (sphere) for ligaments segment (cylinder) (25)

sph

cyl

with CD and CD denoting the drag coefficient of the sphere and cylinder, respectively. Acceleration of both the cylindrical ligament’s segment and the spherical drop is inversely proportional to their sph cyl corresponding diameter d. Since 34 CD  π2 CD , the ligament’s elongated body is subject to faster deceleration than the evolving drop with a notably larger diameter. This results in a larger azimuthal velocity of the drop compared to the remaining near part of the

ligament leading to the numerically and experimentally observed shape of the ligament. Concluding remarks The VOF model and its state-of-the-art implementation with the Gerris code can successfully predict the growth phase of the ligaments. The main issue for the accurate simulation of the ligaments’ growth is obtaining sufficient grid density near the liquid–air interface which can be efficiently achieved with the AMR approach. Such simulation closely coincides with the realistic growth of the ligaments as long as the ligaments are covered by at least four cells in the transverse direction. The growth starts from embryos – small bulges rising from the rotating film in a radial direction. Experimental evidence shows various sizes and positions of ligaments indicating that the appearance of embryos is random. In this work we did not study the formation of embryos; instead we started the simulations after they had already formed. This approach enables relatively simple (i.e. without statistical analysis) comparison of the numerical results with the experimental data. Notwithstanding the assumptions made in setting the initial and boundary conditions (rigid body rotation of the film and the reduction of azimuthal velocity as rn ), the simulated development of ligaments closely resembles the experimental one. Because the surface tension opposes the growth, an embryo – initiated either by the liquid feed fluctuations, the wheel vibrations, and/or interaction with the surrounding air – evolves into a ligament only if the wheel rotational frequency f exceeds the critical value f∗ . The simple analytical assessment of f∗ suggests that f∗ is inversely proportional to the size of the embryo, expressed with the equivalent (hemisphere) radius. The results of the simulations agree with the estimated f∗ , indicating that the estimate is appropriate.

J. Mencinger et al. / International Journal of Multiphase Flow 77 (2015) 90–103

The effect of changing operating variables on the ligament’s growth is investigated by varying the film thickness h and f in the simulations. The ligament grows faster and longer with larger h since more liquid is “available” in the vicinity of its foot to be drawn into it. With larger f, the growth is faster – also in terms of nondimensional time ft; this is not surprising since the centrifugal force is proportional to f2 . By increasing f, the ligament grows longer and thinner while its shape approaches the kinematic shape, i.e. the hypothetical limiting shape when capillary, viscous, and drag forces are neglected. However, as the ligament becomes thinner, the importance of the aerodynamic drag increases due to the fact that the deceleration on each ligament’s segment is inversely proportional to its diameter d. The shape of the ligament in its later phase of growth is thus notably influenced by the drag force. The drag is also responsible for the typical forward curvature near the ligament’s tip caused by the relatively faster azimuthal velocity of the drop evolving at the tip. Even with the unsurpassed effectiveness of the AMR approach to obtain a dense grid near the interface, the simulation of the ligament breakup still requires larger computational resources than used in this work. The simulations of the complete life of the ligaments will be pursued in the future. This work also presents a step towards numerical simulation of the industrial fiberization process – the later requires the inclusion of the transport equation for enthalpy, and modeling the material properties of the liquid. Acknowledgments This work was in part supported by Slovenian Research Agency (ARRS), grants P2-0162, P2-0167, and L2-4270. We are also grateful to the prof. S. Popinet and others contributing to the freely available Gerris code. References Bizjan, B., Širok, B., Govekar, E., 2015. Nonlinear analysis of mineral wool fiberization process. ASME. J. Comput. Nonlinear Dynam. 10 (2) 021005-021005-8. ´ A., 2014a. Ligament-type liquid disintegration Bizjan, B., Širok, B., Hoˇcevar, M., Orbanic, on a spinning wheel. Chem. Eng. Sci. 116, 172–182. ´ A., 2014b. Liquid ligament formation dynamBizjan, B., Širok, B., Hoˇcevar, M., Orbanic, ics on a spinning wheel. Chem. Eng. Sci. 119, 187–198.

103

Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling surface tension. J. Comput. Phys. 100, 335–354. Cummins, S.J., Francois, M.M.B.K.D., 2005. Estimating curvature from volume fractions. Comput. Struct. 83, 425–434. Eisenklam, P., 1964. On ligament formation from spinning discs and cups. Chem. Eng. Sci. 19 (9), 693–694. Fermigier, M., Limat, L., Wesfreid, J.E., Boudinet, P., Quilliet, C., 1992. Two-dimensional patterns in Rayleigh–Taylor instability of a thin layer. J. Fluid Mech. 236, 349–383. Hirt, C.W., Nichols, B.D., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 39, 201–225. Kamiya, T., 1972. Analysis of the ligament-type disintegration of thin liquid film at the edge of rotating disk. J. Chem. Eng. Jpn. 5, 391–396. Landau, L.D., Lifshitz, E.M., 1987. Fluid Mechanics, second ed. Pergamon Press. Lefebvre, A., 1989. Atomization and Sprays. Hemisphere, New York. Liu, H., 2000. Science and Engineering of Droplets: Fundamentals and Applications. Noyes Publications, Norwich, New York, U.S.A. Liu, J., Yu, Q., Guo, Q., 2012a. Experimental investigation of liquid disintegration by rotary cups. Chem. Eng. Sci. 73, 44–50. Liu, J., Yu, Q., Li, P., Du, W., 2012b. Cold experiments on ligament formation for blast furnace slag granulation. Appl. Therm. Eng. 40 (1), 351–357. Olesen, M., 1997. Prediction of Drop-size Distributions Based on Ligament Breakup (Ph.D. thesis). Queen’s University at Kingston, Kingston, Ontario. Popinet, S., 2003. Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries. J. Comput. Phys. 190, 572–600. Popinet, S., 2009. An accurate adaptive solver for surface-tension-driven interfacial flows. J. Comput. Phys. 228, 5838–5866. Purwanto, H., Mizuochi, T., Akiyama, T., 2005. Prediction of granulated slag properties produced from spinning disk atomizer by mathematical model. Mater. Trans. 46 (6), 1324–1330. Renardy, Y., Renardy, M., 2002. PROST: A parabolic reconstruction of surface tension for the volume-of-fluid method. J. Comput. Phys. 183, 400–421. Roy, C.J., 2005. Review of code and solution verification procedures for computational simulation. J. Comput. Phys. 205, 131–156. Scardovelli, R., Zaleski, S., 1999. Direct numerical simulation of free-surface and interfacial flow. Annu. Rev. Fluid Mech. 31, 567–603. Shinjo, J., Umemura, A., 2010. Simulation of liquid jet primary breakup: Dynamics of ligament and droplet formation. Int. J. Multiph. Flow 36 (7), 513–532. Taylor, G.I., 1950. The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. Proc. R. Soc. 201 (1065), 192–196. Tong, A.Y., Wang, Z., 2007. Relaxation dynamics of a free elongated liquid ligament. Phys. Fluids 19 (9), 092101. Ubbink, O., Issa, R.I., 1999. A method for capturing sharp fluid interfaces on arbitrary meshes. J. Comput. Phys. 153, 26–50. ´ A., Bajcar, T., 2014. Mineral wool melt fiberization on a Širok, B., Bizjan, B., Orbanic, spinner wheel. transactions of the institution of chemical engineers. Part A Chem. Eng. Res. Des. 92 (1), 80–90. ´ B., Bullen, P., 2008. Mineral Wool: Production and Properties. Širok, B., Blagojevic, Woodhead Publishing Limited, Cambridge. Westerlund, T., Hoikka, T., 1989. On the modeling of mineral fiber formation. Comput. Chem. Eng. 13 (10), 1153–1163. Youngs, D.L., 1982. Time-dependent multi-material flow with large fluid distortion. In: Morton, K.W., Baines, M.J. (Eds.), Numerical Methods for Fluid Dynamics. Academic Press, New York, pp. 273–285.