Comparison of secondary breakup models for droplet-laden compressor flows

Comparison of secondary breakup models for droplet-laden compressor flows

International Journal of Multiphase Flow 116 (2019) 125–136 Contents lists available at ScienceDirect International Journal of Multiphase Flow journ...

4MB Sizes 0 Downloads 33 Views

International Journal of Multiphase Flow 116 (2019) 125–136

Contents lists available at ScienceDirect

International Journal of Multiphase Flow journal homepage:

Comparison of secondary breakup models for droplet-laden compressor flows C. Storm a,b,∗, F. Joos a,b a b

Laboratory of Turbomachinery, Helmut-Schmidt-University, Germany University of the Federal Armed Forces Hamburg, Germany

a r t i c l e

i n f o

Article history: Received 25 September 2018 Revised 3 April 2019 Accepted 6 April 2019 Available online 15 April 2019 Keywords: Two-phase flow Wet compression Compressor flow Droplet-laden flow Droplet deformation Droplet fragmentation Breakup criterion Numerical modelling Euler-Lagrange approach DPF code Atomization Secondary breakup TAB NLTAB3

a b s t r a c t Different secondary breakup models are assessed to simulate gas-droplet flows that involve large velocity gradients and yield relative velocities between both phases (e.g. wet compression flows). For the simulation of such two-phase flows, the Euler-Lagrange approach is employed, using a commercial flow solver and an in-house FORTRAN source code. The droplet deformation models: Taylor Analogy Breakup (TAB) and a simplified implementation of the Non-linear TAB (NLTAB3) are each tested with two different breakup criteria to evaluate the adequacy for considered droplet-laden compressor flows. A detailed atomization model is applied to account for a temporal fragmentation process. By contrast with experimental results good agreements are found for Weber number distributions in the inter blade flow. Though, different sensitivities and locations for secondary breakup occurrence are observed in dependence of applied droplet deformation model and breakup criterion. In spite of incisive model simplifications to reduce its computational effort, the non-linear droplet deformation model (NLTAB3) in combination with a physically substantiated breakup criterion showed good agreements with qualitative experimental data.

1. Introduction The occurrence of secondary breakup is a prevailing issue within liquid atomization processes and hence of importance in many technical applications that use sprayers, atomizers or nozzles for spray coating, spray combustion, spray drying or evaporative cooling. Principally the atomization is provoked in gas-liquid flow interactions due to high relative velocities between the phases. Whereas primary breakup considers the first disintegration of the injected bulk liquid into filaments and drops, secondary breakup is regarded as the proceeding fragmentation of solitary droplets into smaller ones. Secondary breakup occurs, whenever aerodynamic forces acting on a droplet provoke its deformation to a certain degree that it cannot return to its original form and disintegrates. The deformation alters the droplet drag coefficient and in consequence, the droplet trajectory as well as the breakup position and

Corresponding author. E-mail addresses: [email protected] (C. Storm), [email protected] (F. Joos). 0301-9322/© 2019 Elsevier Ltd. All rights reserved.

© 2019 Elsevier Ltd. All rights reserved.

mode, ultimately affecting the droplet size distribution within the two-phase flow. At the laboratory of turbo machinery of the Helmut-SchmidtUniversity, University of the Federal Armed Forces Hamburg, fundamental research concerning various effects of wet compression applied for power augmentation of gas turbines is undertaken (cf. Ulrichs and Joos, 2006; Eisfeld and Joos, 2009 and Neupert et al., 2014). Wind tunnel experiments are performed to identify dropletladen flow phenomena and to investigate their influence on the aerodynamic performance of a planar compressor cascade. Qualitative results revealed secondary breakup occurrence due to strong velocity gradients in the transonic regime on the blade suction side, as well as in front of the leading edge and in the wake. Supplementary to aforementioned experiments, numerical simulations are conducted employing the in-house Disperse Phase Fortran (DPF) source code introduced in Matysiak (2007) and Storm and Joos (2011). Utilizing the Euler-Lagrange method the DPF code predicts flow behavior of the disperse phase (i.e. water droplets) while concurrently communicating with the ANSYS CFX 13.0 flow solver, which simulates the continuous phase (i.e. airflow). Implemented droplet-specific models consider turbulent


C. Storm and F. Joos / International Journal of Multiphase Flow 116 (2019) 125–136

dispersion, flow induced secondary breakup and droplet-wall interaction, while evaporation remains neglected. The influence of droplets on the compressor flow characteristics is simulated in two-way coupled computations. Generally good agreement to experimental results is yielded (cf. Matysiak, 2007 and Storm and Joos, 2011). Yet, numerical models - and in this case the secondary breakup modeling - are continuously optimized for more accurate simulation results. Secondary breakup modeling involves the computation of deformation, breakup criterion as well as fragmentation for the definition of secondary droplet properties. Various empirical and mechanistic models are given in literature, but none alone capable to predict the entire process accurately. Although mathematically simple and numerically efficient, empirical models are omitted, as these are tied to original experimental data and constraints. Guildenbecher et al. (2009) and Ashgriz (2011) reviewed different mechanistic models and recommend the Taylor Analogy Breakup (TAB) and Droplet Deformation Breakup (DDB) modeling in regard of objective simulation application (i.e. without spray atomization and at prevailing moderate Weber number regimes). However, unrealistic predictions are observed by Schmehl (2003) and Park and Lee (2004) for the DDB model introduced by Ibrahim et al. (1993). Park et al. (2002) and Lee et al. (2012) concluded that the TAB model still outperforms the DDB model, so that the implementation of the TAB model is processed. Further, the NLTAB3 model proposed by Schmehl (2003) is implemented, which is considered more sophisticated as it accounts for nonlinearities within deformation calculations. Both deformation models are contrasted in a case study at shock loading conditions by Matysiak (2007), verifying good agreement with experiments for both models, but observing slightly larger oscillation frequency for the TAB model. Yet, the necessity of a time step size small enough to solve non-linear differential equations of the NLTAB3 model limits its practicability for technical flow applications. In this regard, a simplified implementation is chosen that considers constant flow conditions within each time step of motion. The original TAB model introduced by O’Rourke and Amsden (1987) applies a constant maximum deformation as breakup criterion. This is considered unrealistic as Krzeczkowski (1980) detected larger maximum deformation for increasing Weber numbers. Matysiak (2007) proposed an empirical criterion dependent on the Weber number, but based on shock tube results, which show a more abrupt aerodynamic loading than observed within considered technical application of compressor flows, questioning its applicability. Instead, a more physically substantiated breakup criterion proposed by Park et al. (2002) show good prospects for use in flows characterized by rather gradual than shock loading conditions. The implementation of secondary breakup modeling is subdivided in separate stages, so that different models regarding deformation, breakup criterion and fragmentation can be combined in a new way. The objective of this paper is to test various model combinations in regard of considered technical application of twophase flow (e.g. droplet-laden compressor flow) and to compare to recent experimental results presented in Neupert (2017). Due to limited accessibility within such experiments, results are contrasted to overall Weber number distributions and shadowgraphs of located secondary breakup occurrences. For the latter, the temporal evolution of the droplet fragmentation process needs to be considered, so that a detailed atomization model according to Bartz et al. (2010) is implemented. 2. Secondary breakup modeling In the following, the droplet deformation models (i.e. TAB and NLTAB3), the different breakup criteria as well as the detailed at-

omization model are briefly outlined. For further information, it is referred to referenced literature. 2.1. Droplet deformation modeling Whereas aerodynamic forces trigger deformation, surface tension and liquid-viscous forces counteracts. Due to the dependence on relative velocity, droplet diameter and fluid properties, the nondimensional Weber number We and Ohnesorge number Oh are applied to characterize droplet deformation and breakup mechanisms (cf. Pilch and Erdman, 1987; Hsiang and Faeth, 1995 and Schmehl, 2003). The subscripts d and f consistently denote droplet and continuous phase specific variables respectively. 

(urel )2 dd ρ f ( u f − ud )2 dd ρ f We = = σ σ Oh =

ηd ρd d d σ


√ We ≈ Red


The classical TAB model proposed by O’Rourke and Amsden (1987) is based on Taylor’s analogy relating droplet deformation behavior to an oscillating mass-spring-damper system (cf. Eq. (3)). Herein, xd refers to the drop equator displacement and md to the droplet mass. The inertia term is given by m x¨d , the damping or rather friction term cx˙ d corresponds to the droplet viscosity force, the restoring spring term k xd to the surface tension force and the external force F to the aerodynamic loading.

md x¨d + cx˙ d + kxd = F


It is proposed by Schmehl (2003) to describe the equatorial disd placement as a dimensionless deformation in terms of yd = d d (cf. d ,0

Fig. 2 (left)) and furthermore to normalize the timescale in terms 3 ρd dd, 0

of Tσ = τtσ with τσ2 = σ defined as the characteristic time of free oscillations. In this regard the non-homogeneous second order linear differential equation for free particle oscillation given in O’Rourke and Amsden (1987) reformulates to 2

d yd dy + Cc Oh d + Ck (yd − 1 ) = 4CF W e d Tσ d Tσ2


The dimensionless model coefficients are determined by O’Rourke and Amsden (1987) from theoretical and experimental results, but these are only partially adopted. Schmehl (2003) claims a mistake in the choice of the damping coefficient Cc , so that ultimately the coefficients are defined with Cc = 40, Ck = 64 and CF = 13 . It is argued by Matysiak (2007) that the temporal resolution of all deformation oscillations is not reasonable, as this would necessitate a reduction of the computational time step for droplet motion, which is evaluated from the characteristic droplet response time, by factor 10−3 . The author therefore introduced a mean deformation Wem yielding following analytical solution for droplet deformation.

yd (Tσ ) = W em + e +





y˙ d,0 +

yd,0 − W em cos (ωTσ )

yd,0 − W em

4CF W e + Ck , Ck 1 1 = Cc Oh, τμ 2 1 ω2 = Ck − 2


 sin (ωTσ )


with: W em =



C. Storm and F. Joos / International Journal of Multiphase Flow 116 (2019) 125–136

A general shortcoming of the TAB model is that acting forces (i.e. restoring, damping and external forces) in Eq. (4) are related to the undeformed droplet. Yet, this does not apply to the drag force, which is determined according to the dynamic drag model proposed by Liu et al. (1993). This computes the drag coefficient Cdrag from a linear interpolation between that of a spherical and a disc shaped droplet. Nonetheless, with progressive deformation the interdependencies of acting forces and droplet shape becomes significant. Accounting for these nonlinearities, Schmehl (2003) developed three nonlinear TAB so-called NLTAB models, for which the approximation of an ellipsoid to the deforming droplet is applied to a certain number of terms. The NLTAB3 model considers for nonlinearities within each term: the inertia, the friction, the restoring and the aerodynamic term. The inertia term and aerodynamic term are considered proportional to the sum of a radial and a longitudinal contribution. For the friction term a velocity potential is assumed constant due to symmetry throughout the volume, yielding a y−2 nonlinearity. To describe the potential energy d of the droplet surface for the restoring term, analytical geometry is utilized to determine the ellipsoidal surface. Ultimately, the model equation is given (in accordance to Eq. (4)) as a function of the dimensionless deformation yd and the dimensionless time Tσ with

π2 +

16 y6d


d yd − π 2 + 16 d Tσ2 + 20

48 1 π 2 + 16 y7d

d yd d Tσ


+ 40Oh

1 d yd y2d d Tσ

1 dS 2C2 = We S0 d yd yd


2.2. Breakup criterion Due to unknown prospective aerodynamic loading, a breakup criterion solely based on the instantaneous Weber number is not expedient. Instead, the maximum deformation yd,max has established as breakup criterion, which defines the deformation of the disc-shaped droplet just before the onset of bulging or disintegration. Under assumption of constant flow conditions within a given time step t, breakup occurs within the first amplitude of oscillation calculation. Referring to Matysiak (2007) an analytical solution for the time of maximum amplitude can be derived from the first derivation for the TAB model, yielding Eq. (8). In the case that tmax < t and yd (tmax ) ≥ yd,max and in the case of t < tmax and yd (t) > yd,max breakup occurs.

tmax = ω1 tan−1

( )2 + π (yd,0 −Wem )(1+ω2 (τσ τμ )2 )+y˙ d,0 (τσ τμ ) y˙ d,0 ω τσ τμ

ceeded until the end of the computational time step t. Determined final deformation yd (t + t ) and its velocity y˙ d (t + t ) are stored as initial values for the next time step computation. A constant maximum deformation yd,max = 1.5 is specified by O’Rourke and Amsden (1987) for proposed TAB model. This has been confirmed by Schmelz (2002) in experiments with accelerated gas flows for the critical Weber number W ecrit = 12, above which secondary breakup occurs in the first place. Schmehl (2003) reports of shock tube experiments that identified yd,max = 1.8 for W ecrit = 13. Similar to Krzeczkowski (1980), also Schmehl (2003) observed increased maximum deformations for larger Weber numbers, e.g. yd,max = 2.1 for We > 20. Based on these results, Matysiak (2007) proposes an adaptive maximum deformation determination with following equations.

W e ≤ 13

yd,max = 1.771

W e > 13

yd,max = 1.771 + 0.329 1 − e−(


Considering the application of the NLTAB3 model, maximum amplitude is reached, when its derivation y˙ d (t ) changes sign. If corresponding deformation exceeds the defined maximum deformation yd,max i.e. yd (t) ≥ yd,max the breakup criterion is met. Otherwise, the calculation of dimensionless deformation yd (t) is pro-


We−13 2




A more physically substantiated breakup criterion is proposed by Park et al. (2002), considering the equilibrium of the external gas pressure, the interfacial surface tension pressure and the internal liquid pressure at the stagnation point and the equator. The authors assume breakup when the gas pressure difference between stagnation point and equator is lower than the surface pressure difference and derive following breakup condition:

2 1+

To solve Eq. (7) for the temporal resolution of deformations a Runge–Kutta algorithm is applied. Schmehl (2003) recommends evaluating the time step from the characteristic time scale of free oscillation τ σ . But for the synchronous solution of the droplet motion this results in an unreasonable expenditure of computation time, e.g. presented simulations afford approximately ten times the CPU time compared to the TAB model. Hence, a simplification is considered whereupon droplet deformation is computed separately from droplet motion. For droplet motion, the computational time step remains evaluated from the characteristic droplet response time. Assuming constant flow conditions within this time step, droplet deformation according to Eq. (7) is computed separately with the time step evaluated from the time scale of free oscillation τ σ as recommended by Schmehl (2003).


1 y 2 d,max


+ 1+

1 y 2 d,max


−4 1+

1 y 2 d,max


> 2.5107W e

(11) Compared to the experimental data from horizontal air streams presented in Krzeczkowski (1980), this model showed good agreements for the bag and shear breakup regime and conforms with the experimental data of Schmelz (2002) i.e. starting with a maximum deformation of yd,max = 1.5 at the critical Weber number W ecrit = 12. As shown in Fig. 1 the breakup criterion by Park et al. (2002) allows lower values of maximum deformation for We < 65 compared to the empirical equations given in Matysiak (2007). 2.3. Fragmentation modeling Referring to Pilch and Erdman (1987), the breakup mode is defined by the Weber number (cf. Fig. 2 (right)). Each mode shows a particular disintegration behavior concerning breakup times as well as number and properties of secondary droplets. Regarding prevalent Weber numbers observed within considered flow region of simulated experiments, solely the bag breakup and bag and plumen (so-called multimode breakup) are of interest within this work. Since vibrational breakup occurs contingently upon aerodynamic forces, when the droplet interacts in such a way that oscillations develop at the natural frequency of the droplet, ensuing the breakup into few large fragments, this breakup mode remains generally unconsidered in secondary breakup studies. Corresponding to defined breakup modes, Schmehl (2003) outlined experimental results (primarily from shock tube experiments) from Dai and Faeth (2001), Krzeczkowski (1980) and Wierzba (1990), illustrating the temporal stages of droplet deformation and atomization over varying Weber number in Fig. 3. Consequently, to consider for definite breakup times and hence locations an explicit empirical modeling of the fragmentation process has been implemented according to the detailed atomization model introduced in Bartz et al. (2010). Incorporating this data from Fig. 3 into the code, the dimensionless times for maximum deformation Tymax as well as starting and ending breakup times Tb,i for different stages (i.e. i = bag, rim,


C. Storm and F. Joos / International Journal of Multiphase Flow 116 (2019) 125–136

Fig. 1. Maximum droplet deformation at breakup as a function of Weber number according to the model by O’Rourke and Amsden (1987), Matysiak (2007) and Park et al. (2002).

Fig. 2. Schematic droplet deformation (left) and breakup mechanisms classified according to Weber numbers (adopted from Pilch and Erdman, 1987) (right).

plume and sheet thinning) are tabulated for distinct Weber numbers W e0 = 12,15, 20, 25, 33, 40, 50, 60, 80 and 100. For intermediate Weber numbers, the dimensionless times are determined by linear interpolation. However, these Weber numbers refer to the initial Weber numbers We0 (i.e. at the start of shock loading), which differ significantly from the Weber number evaluated at breakup condition Weymax . As a consequence of droplet acceleration with accompanied droplet deformation, the relative velocity decreases reducing the instantaneous Weber number until the breakup criterion is met. In order to apply available experimental results, Bartz et al. (2010) proposed a correlation in terms of the Weber number Weymax , Reynolds number Reymax and Ohnesorge number Ohymax at maximum deformation to evaluate the equivalent initial Weber number We0 . For the application of this correlation in computations with gradual aerodynamic loading (instead of shock loading), it is assumed that the Weber number at the maximum deformation Weymax as well as the proceeding fragmentation process are equivalent for both types of loading conditions. Bartz et al. (2010) further propose correlations to evaluate the dimensionless time Tymax to evaluate the fictive onset of loading. Yet, this calculation is omitted, instead the dimensionless time Tmin

denoting the time, when the disc shaped droplet reaches its minimal thickness, is considered as the time of maximum deformation Tymax (cf. Fig. 3). Shaded areas in Fig. 3 determine the time periods Tb_end,i − Tb_init,i , in which the breakup of specified stages are observed to occur. To identify a discrete breakup time Tbu,i for an individual stage (e.g. i = bag, ring, plume, core), this is considered randomly distributed within defined time spans. Exemplified for a bag breakup at W e = 14 from qualitative results of Kulkarni and Sojka (2014) the time bar with dimensionless time definitions is illustrated in Fig. 4. With proceeded tracking of a dummy droplet after met breakup criterion, the locations of secondary droplet formation are determined with the remaining times Tbu,i − Tymax for each individual disintegration stage (i.e. i = bag, ring, plume and core). Corresponding to the disintegrating stage, individual secondary droplet properties are determined from results of a differentiated analysis of experimental data for secondary droplet size, number and velocity presented in Schmehl (2003). As secondary droplet properties are not within the focus of examination and validation of simulated experiments, the description of the modeling is omitted within this work. For further information on the modeling for secondary droplet properties it is referred to the work of

C. Storm and F. Joos / International Journal of Multiphase Flow 116 (2019) 125–136


Fig. 3. Temporal evolution of droplet deformation and atomization (Oh<0,1), adapted from Schmehl (2003) (Experimental data: Dai and Faeth, 2001, Krzeczkowski, 1980 and Wierzba, 1990; red dashed line: Pilch and Erdman, 1987). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. Temporal sequence of deformation and individual stage breakup, (Images adapted from Kulkarni and Sojka, 2014).

Bartz et al. (2010), who applied this data in a detailed atomization model and observed good agreement with experiments from Park et al. (2006). 3. Experimental and computational setup To evaluate different variations of secondary breakup modeling with regard to focused application in droplet-laden compressor flows, the outcomes of numerical simulations are contrasted to results from appropriate experiments published in Neupert (2017). Hence, the experimental and corresponding computational setup are presented in the following. 3.1. Experimental setup Experimental research on the behavior of a two-phase air and water droplet flow in a transonic compressor rotor blade passage was carried out by Neupert (2017). The test facility employs an open loop wind tunnel, which comprises a settling chamber, a transonic planar compressor cascade and a water separator. A radial compressor feeds compressed air into the settling chamber at Ttot = 304K , where sieve plates equalize the flow and a manifold equipped with high-pressure pin jet nozzles inject decalcified

and demineralised water mist. Followed by an acceleration section, the droplet-laden flow is homogenized before entering the test section, which is equipped with a wall boundary layer treatment as reported in Eisfeld and Joos (2009) to eliminate the wall boundary layer effects. Within the test section (cf. Neupert et al., 2014), eight transonic aerofoils that are extruded outer cross sections of the first row of an experimental axial compressor are mounted between two acrylic glass windows for optical measurement accessibility. The geometric properties of the cascade and design flow conditions are presented in Fig. 5 (left) and corresponding geometric definitions and measurement positions are shown in Fig. 5 (right). The cascade blading has been examined at varying incidence angle β1 = 149◦ − 155◦ operated at an inlet Mach number Main = 0.95 measured at the front traverse and a relative outlet static pressure pstat = 0. As illustrated in Fig. 5 (right), measurements are taken between cascade inlet and outlet, respectively front and back measurement traverse. Nonintrusive measurement techniques Laser Doppler and Phase Doppler Anemometry (LDA/PDA) are used to measure droplet velocity and size within the measurement zone. Shadowgraph is applied to visualize inhomogeneity within the flow, such as droplets or shock waves. A water load of χw = 1.0% is applied for the droplet-laden


C. Storm and F. Joos / International Journal of Multiphase Flow 116 (2019) 125–136

Fig. 5. Compressor cascade data (left) and geometric definitions and measurement positions (right) (adopted from Neupert 2017).

Fig. 6. Computational domain within experimental setup (left) with boundary conditions (right).

compressor flow. Corresponding droplet spectrum measured at the inlet of the cascade is depicted in Fig. 7 (left) and spray characteristics by means of representative diameters: arithmetic mean diameter dd,10 , Sauter mean diameter dd,32 and the volume median diameters dd V0.1 , dd V0.5 and dd V0.9 are tabulated in Fig. 7 (right). For more information on experimental setup and procedure it is referred to Neupert (2017). 3.2. Computational setup For numerical in-depth examination of the two-phase flow characteristics within the transonic compressor cascade and alongside curtailed numerical efforts the simulation involves a quasitwo-dimensional domain containing the three center blades (cf. Fig. 6 (left)). The corresponding mesh grid generated with the commercial ANSYS ICEM CFD software consists of 163,0 0 0 to 168,0 0 0 hexahedral cells with blades embedded using O-grid structures. Concerning mesh quality, a minimum determinant of 0.75 and a minimum angle of 22.7◦ are yielded. Within the commercial ANSYS CFX 13.0 environment, boundary conditions are defined according to the experimental operating conditions (cf. Fig. 6 (right)). The computational grids are designed for the use of the k-ε turbulence model with the scalable wall function approach. Exploiting symmetry and periodicity a significant reduction of the computational effort is achieved. Total pressure and total temperature conditions at the inlet and static pressure condition at the outlet are adjusted to match the flow conditions as specified by experiments. The blades are defined with a smooth no slip wall boundary condition. The continuous phase flow is modeled with air ideal gas ap-

plying the total energy model. Additional momentum source terms are used to account for the effect of the disperse phase; yet mass transfer remains unconsidered. The droplet flow is computed separately within the DPF code. The water mass is simulated with 20 0,0 0 0 droplet parcel trajectories starting uniformly distributed from the inlet. The droplet size distribution is modeled with a combination of two Rosin–Rammler approaches given in Eq. (12). The parameters are fitted to the droplet size spectrum observed within experiments at the cascade inlet (cf. Fig. 7) with λ1 = 22, λ2 = 4.4, k1 = 1.42 and k2 = 1.6. Resulting spray characteristics are tabulated in Fig. 7 (right).

1 P ( dd ) = 2dd



k1 e



− λd 1





k 2

k2 e


− λd 2

k 2 


To simulate droplet slip velocities measured within experiments (cf. Neupert 2017), numerical analyses were conducted yielding 

correlations for the evaluation of the droplet inlet velocity u d with 

Eq. (13) in dependence of the mean continuum velocity u f,mean 

and droplet diameter dd . The droplet velocity difference u d is defined applying the Gaussian distribution function. For the oper

ating point of Main = 0.95, the mean velocity difference uμ and 

standard deviation σ are correlated to the droplet diameter dd according to Eqs. (14) and (15) respectively.  ud

= u f,mean + u d


C. Storm and F. Joos / International Journal of Multiphase Flow 116 (2019) 125–136


Fig. 7. Inlet droplet spectrum from experiments (left) and measured and modelled spray characteristics (right), (cf. Neupert 2017; Neupert et al., 2016).

Fig. 8. Weber number distribution at evaluation line of droplet classes over relative x/c-coordinate from experiments (left) and from simulations (right).

 uμ (dd ) = 0.00443dd2 − 0.90221dd − 3.82938 uμ = v μ ( d d ) = 0 w μ ( d d ) = 0

 σu (dd ) = 0.0 0 064dd2 − 0.08697dd + 4.19169  2 σ = σv (dd ) = 0.0 0 021dd − 0.04106dd + 2.52718 σw (dd ) = 0.0 0 035dd2 − 0.05775dd + 2.88276 



With the particle equation of motion, the droplet trajectory is computed within the Lagrangian frame of reference taking drag, gravity, pressure gradient and buoyancy forces into account. Furthermore, particle specific phenomena such as turbulent dispersion and droplet-wall interaction are applied within the disperse phase computation. For more detailed information on phase coupling, droplet trajectory computation and turbulent dispersion and droplet-wall interaction modeling it is referred to Matysiak (2007) and Storm and Joos (2011). 4. Simulation results In order to investigate secondary breakup occurrence within accelerated flows, Neupert (2017) examined the high-speed regime on the suction side within the compressor cascade at the incidence angle β1 = 154◦ and with inlet Mach number Main = 0.95.

Along the evaluation line shown in Fig. 5 (right), the Weber number of individual droplet classes were quantified and plotted over the relative x/c-coordinate in Fig. 8 (left). The experimental results reveal that only droplets of the largest size attain Weber numbers above W ecrit = 12, which Pilch and Erdman (1987) defined as the bag breakup regime. Good agreements are obtained for the Weber number curves of computational results as depicted in Fig. 8 (right). Yet, higher Weber numbers are yielded in the upstream flow particularly for the larger droplet classes, which is accounted to overestimated droplet inlet velocities. Smaller droplet classes are less affected by an overestimation, because these approximate quickly to the air flow velocity. Nonetheless, the minor discrepancy is neglected, in regard of the wellmatched Weber numbers for the most part i.e. for the coordinates x/c ≥ 0.125. Moreover, Neupert (2017) investigate the influence on the Weber number curves considering the largest droplet class dd ≥ 60μm for varying incident angle 149◦ ≤ β 1 ≤ 155◦ . The results shown in Fig. 9 reveal that the curves change accordingly to the high-speed region on the suction side i.e. with decreased angle β 1 , the curve on the suction side become wider, lower and shifts further downstream towards the trailing edge (TE). For lower incident angles β 1 ≤ 152◦ the Weber numbers remain below the critical value W e < W ecrit = 12, so that only vibrational breakup may occur.


C. Storm and F. Joos / International Journal of Multiphase Flow 116 (2019) 125–136

Fig. 9. Measured Weber number of droplet class dd ≥ 60 μm over relative x/c-coordinate at varied incident angles 149◦ ≤ β 1 ≤ 155 (Neupert 2017).

Fig. 10. Simulated Weber number at analysis line over relative x/c-coordinate for droplet sizes dd ≥ 60 μm (left) and corresponding contour plots with marked breakup positions for incident angles β 1 =152◦ and 153◦ (right).

For larger incident angles β 1 ≥ 153◦ an abrupt increase in Weber numbers to W e ≥ W ecrit = 12 is observed, which is decisive for secondary breakup by means of bag breakup. Against the trend, larger maximum Weber numbers are observed for incident angle β1 = 154◦ compared to the incident angle β1 = 155◦ . Neupert (2017) speculates this to be due to strong deformed droplets, which cannot be detected by the PDA measurement system. Within corresponding numerical results depicted in Fig. 10 (left) this discrepancy is not observed. Overall good agreement for the curve progression within the relative coordinate range of x/c = 0.125 − 1.0 is obtained. Fig. 10 (right) show corresponding simulation results for Weber number contour plots with secondary breakup positions (black spheres) computed with implemented TAB model and the applied breakup criterion according to Park et al. (2002). As depicted for the incident angle β1 = 153◦ secondary breakup occurs above the evaluation line (black dashed line) due to Weber numbers We ≥ 12. Compared to the results for the incident angle β1 = 152◦ a larger number of secondary breakup is detected within the high-speed regime below the evaluation line, where higher Weber numbers are obtained.

Above simulation results, depicting Weber number distributions in dependence of different droplet sizes and varying incident angles, show overall good agreement with experiments. Only larger droplets possibly undergo secondary breakup within the transonic regime, provided the incident angle is large enough to induce high aerodynamic forces. In consideration of the droplet spectrum (cf. Fig. 7) introduced to the flow, this consequently applies only to a very small number of droplets. Therefore, the effects of secondary breakup on the droplet size and Weber number distributions remain unverifiable, but the evidence of conformity for flow conditions of examined droplet spectrum and aerodynamic load cases is provided. With regard to above findings, following simulations investigate compressor flows with the incident angle β1 = 155◦ to initiate a high aerodynamic load case. For comparison of applied deformation models the maximum deformation yd,max is evaluated within each time step of motion, without setting the breakup criterion, in order to keep the droplets within the flow. At first, the practicability of utilized simplification for the NLTAB3 model is scrutinized, averaging the maximum deformation exclusively for the droplet class 75μm ≤ dd ≤ 85μm, which are significantly affected.

C. Storm and F. Joos / International Journal of Multiphase Flow 116 (2019) 125–136


Fig. 11. Maximum deformation yd,max for droplet sizes 75μm ≤ dd ≤ 85μm along evaluation line (left) and contour plot (right).

Computations reveal that the simplification made for the implementation of the NLTAB3 model reduces the computational effort significantly for considered simulations i.e. in terms of CPU time by approximately a factor of ten. The results obtained along the evaluation line (shown in Fig. 5) are plotted over the relative x/ccoordinate in Fig. 11 (left). It is shown, that the simplified implementation of the NLTAB3 model renders slightly higher maximum deformations yd,max compared to the implementation proposed by Schmehl (2003). Identifying the highest difference between the maximum deformations predicted of different models with yd,max , the largest deviation between the original and the simplified NLTAB3 model implementation is observed at x/c ≈ 0.47 with yd,max ≈ 0.2. Overestimation of the maximum deformation possibly renders overrated secondary breakup occurrence and due to a higher number of secondary droplets to increased CPU time. Nonetheless, both NLTAB3 model implementations yield a similar mapping of maximum deformations throughout the compressor flow area, as illustrated for the simplified NLTAB3 model in Fig. 11 (right (top)). With regard to the reduction of computational effort, the simplified implementation of the NLTAB3 model appears reasonable and preferable. In comparison to the TAB model implemented as proposed by Matysiak (2007), the simplified NLTAB3 model requires approximately 25–30% more CPU time (in neglect of breakup occurrence and secondary droplets). Moreover, contrasting the results for maximum deformation along the evaluation line (cf. Fig. 11(left)), the TAB model primarily predicts lower values. Compared to the original NLTAB3 model, it shows its largest deviation of the maximum deformation yd,max ≈ 0.4 at x/c ≈ 0.55. Although all deformation models display a local increase of the maximum droplet deformation yd,max within the transonic regime on the suction side, the TAB model predicts an earlier growth. Its peak value is followed by a gradational reduction, yielding a minimum, which is predicted further down to the trailing edge compared to the NLTAB3 models. Corresponding contour plots of the maximum deformation yd,max of the simplified NLTAB3 and the TAB model are presented in Fig. 11 (right). The TAB model predicts high values nearby the leading edge, whereas the NLTAB3 models predicts no deformation within this area. This behavior corresponds to the findings from

Fig. 11 (left) that the TAB model predicts maximum deformations ahead of time respectively instantaneously on present aerodynamic load as also identified in Hwang et al. (1996). In contrast, the delay predicted by the NLTAB3 model is assumed due to considered nonlinearities and progression over finite time. The high aerodynamic load nearby the leading edge rather shows its effect on increased maximum deformation values further downstream at x/c ≈ 0.2. Due to post-pulse oscillations and further aerodynamic loading similar values are regained at x/c ≈ 0.5. Ultimately, the maximum deformation defines in which regions secondary breakup occurs, as the breakup criteria are based on the degree of deformation. In the following, the breakup criteria given by Matysiak (2007) in Eqs. (9) and (10) and by Park et al. (2002) in Eq. (11) are applied in combination with the TAB and the simplified NLTAB3 model. In Fig. 12 the contour plot of the Weber number distribution averaged for the droplet class of dd ≥ 60 μm is displayed, chosen because primarily this droplet class undergoes secondary breakup. The positions, where secondary breakup occurrence is identified, are marked with circles colored according to the primary droplet diameter dd, pri , which specifies the initial droplet diameter before secondary breakup. Secondary breakup is identified when droplets enter the high-speed regime on the suction side, attaining high Weber numbers due to high slip velocities. Corresponding to the results for the maximum deformation yd,max predicted in Fig. 11 for the droplet class 75 μm ≤ dd ≤ 85 μm the breakup positions are located further downstream for the NLTAB3 model. In lieu thereof, droplets of lower diameter show similar breakup positions for both deformations models. Obviously, the delay effect of breakup increases with the droplet diameter, when nonlinearities are considered within the progression of deformation. Moreover, the NLTAB3 model identifies additional breakup occurrences in the second half of the high-speed regime in the near wall flow. This is accounted to the effect of post-pulse oscillations and further aerodynamic loading also observed in Fig. 11. Although the TAB model likewise shows this behavior, it does not reach as high values for the maximum deformations and thus no secondary breakup. Regarding applied breakup criteria, the model by Park et al. (2002) is expected to yield overall more secondary


C. Storm and F. Joos / International Journal of Multiphase Flow 116 (2019) 125–136

Fig. 12. Simulated results of breakup positions for incident angle β 1 =155◦ for TAB and NLTAB3 model with breakup criterion by Matysiak (2007) with Eqs. (9) and (10) and by Park et al. (2002) with Eq. (11).

breakup occurrences, as it predicts breakup at lower maximum deformations (cf. Fig. 1). For the TAB model it seems on first sight that more breakup occurrences are observed for the criterion proposed by Matysiak (2007). This is not the case, as the simulation with the criterion by Park et al. (2002) identifies a few more breakup positions but rather in the same spot. However, with the NLTAB3 model the increase of secondary breakup occurrences by use of the criterion by Park et al. (2002) becomes drastically more pronounced. Comparing the simulations with the breakup criterion by Matysiak (2007), it can be assumed from Fig. 11 that the larger region of high maximum deformation values within the NLTAB3 simulation yields a larger number of secondary breakup occurrences. Yet, an antithetic trend is observed, which is due for

two reasons. While the TAB model shows a smaller region of high maximum deformation values in the vicinity of the leading edge, it predicts higher values than the NLTAB3 model. Secondly, the inter blade flow shows larger Weber numbers in the NLTAB3 simulation severely increasing the breakup criterion given by Eq. (10) (cf. Fig. 1). In turn the larger Weber numbers are yielded, due to present higher droplet diameters, resulting from non-initiated secondary breakup in the inlet flow. Hence, it is shown that loading processes within presented accelerated gas flows do rarely attain the degrees of deformation, which are observed for abrupt shock loading and typified within the criterion developed by Matysiak (2007). Qualitative results from experiments published in Neupert (2017) showing secondary breakup detected in the vicinity of the

C. Storm and F. Joos / International Journal of Multiphase Flow 116 (2019) 125–136


Fig. 13. Shadowgraph of secondary breakup and shock for incident angle β 1 =155◦ with detailed views of (1) droplet deformation, (2) bag breakup, (3) ring breakup and (4) vibrational breakup (Neupert 2017).

Fig. 14. Simulated fragmentation process of secondary breakup nearby shock for incident angle β 1 = 155◦ .

shock that terminates the high-speed region on the suction side are presented in Fig. 13. Supplementary, detailed views of temporal stages of a bag breakup such as droplet deformation (1), bag (2) and ring (3) breakup are depicted. In comparison, with Fig. 12, solely simulations that apply the NLTAB3 model in combination with the breakup criterion by Park et al. (2002) reveal secondary breakup occurrences nearby the shock at x/c ≈ 0.45. In regard of subsequent fragmentation process that develops over several stages as shown in Fig. 2 (right) and Fig. 4, secondary droplets evolve downstream of identified breakup position as illustrated in Fig. 14. Contrasting simulation results to experiments in Fig. 13 good agreement is found for the breakup positions of the individual stages. Solely bag and multimode breakup is computed within the high-speed region according to present Weber number regime (cf. Figs. 2 and 3). The deformation into an oblate form and the detected breakup incidence is followed by the formation of a thin hollow bag anchored to a rim. For the multimode breakup, moreover a stamen (or plume) evolves in the center of the bag. The bag breaks up first, yielding a plenty of small droplets, followed by the rim disintegrating into somewhat larger fragments. The plume breakup results in only few larger droplets. The individual breakup of bag, ring and plume and involved ejection of secondary droplets are depicted in Fig. 14. The droplets are visualized as spheres, colored according to the average diameter of the droplet parcel specified as dd,10 . Secondary droplets of smaller diameter (as particularly observed for bag breakup) approximate quickly the airflow. Neupert (2017) determines no experimental results for secondary

droplet properties, so that a comparison is not facilitated on this. 5. Conclusion Presented examinations considered the TAB and NLTAB3 model in combination with the two breakup criteria proposed by Matysiak (2007) and Park et al. (2002). Simulations are compared and contrasted to experimental quantitative and qualitative results presented in Neupert (2017). As a measure for secondary breakup regime (cf. Fig. 2 (right)), Weber numbers are compared along the evaluation line shown in Fig. 5 (right) to assure similar flow conditions. Moreover, the importance of different droplet classes to the Weber number and the effect of varying incidence angles are well reproduced by simulations. The highest Weber numbers are obtained for the largest incidence angle β 1 and droplet class dd ≥ 60 μm. In consequence, the secondary breakup modeling is investigated within the acceleration on the suction side for the incidence angle β1 = 155◦ . The TAB and NLTAB3 models are considered for droplet deformation modeling. Regarding the NLTAB3 model a simplified implementation is applied, which reduces the computational effort significantly and provides its practicability for technical flow applications. Compared to the originally proposed NLTAB3 model implementation, overall good agreement is found for computational results of the maximum deformation yd,max (cf. Fig. 11). Breakup criteria that account for a dependence on the Weber number are applied. Generally larger numbers of secondary breakup occurrences are obtained with the breakup criterion by


C. Storm and F. Joos / International Journal of Multiphase Flow 116 (2019) 125–136

Park et al. (2002), which defines lower maximum deformations yd,max compared to the criterion by Matysiak (2007) (cf. Fig. 1). In combination with applied deformation models, disparities are observed regarding numbers and positions of secondary breakup. A comparison of the maximum deformation distribution obtained for the individual deformation models (cf. Fig. 11 (right)), exhibit further insight on the cause for these disparities. With the TAB model, secondary breakup is detected solely within the first half of the high-speed region, whereas the NLTAB3 model additionally predicts secondary breakup further downstream towards the shock, which terminates the supersonic regime. The disparities are accounted to nonlinearities considered within the NLTAB3 model that refer to interdependencies i.e. the droplet deformation reaffecting the load distribution on the droplet. In comparison to qualitative results by Neupert (2017), which identified secondary breakup occurrence in the vicinity of the shock, most accurate predictions are obtained with the secondary breakup model that combines the NLTAB3 deformation model with the breakup criterion by Park et al. (2002). Further examination of the detailed fragmentation model according to Bartz et al. (2010) show good agreement for predicted positions of temporal fragmentation of bag, ring and plume observed within bag and multimode breakup (cf. Fig. 2 (right) and Fig. 3). In the lack of appropriate experimental data for secondary droplet properties, the corresponding modeling has not been validated. Acknowledgments The authors kindly acknowledge partial financial support by the German Science Foundation (Deutsche Forschungsgemeinschaft) and thank N. Neupert for the provision of experimental data. References Ashgriz, N.E., 2011. Handbook of atomization and sprays. In: Ashgriz, N. (Ed.), Theory and Applications. Springer Science+Business Media LLC, Boston, MA. Bartz, F.-O., Schmehl, R., Koch, R., Bauer, H.-J., 2010. An extension of dynamic droplet deformation models to secondary atomization. ILASS Europe 23rd Annual conference on Liquid Atomization and Spray Systems, Brno, 23, Brno, Czech Republic. Dai, Z., Faeth, G.M., 2001. Temporal properties of secondary drop breakup in the multimode breakup regime. Int. J. Multiph. Flow 27, 217–236. Eisfeld, T., Joos, F., 2009. Experimental investigation of two-phase flow phenomena in transonic compressor cascades, GT2009-59365. In: ASME Turbo Expo 2009: Power for Land, Sea, and Air, Orlando, Florida, USA, volume 7 of Turbomachinery: Parts A and B, pp. 103–112. Eisfeld, T., Joos, F., 2009. New boundary layer treatment methods for compressor cascades. In: ETC 8th European Conference of Turbomachinery. Guildenbecher, D.R., López-Rivera, C., Sojka, P.E., 2009. Secondary atomization. Exp. Fluids 46, 371–402.

Hsiang, L.P., Faeth, D., 1995. Drop deformation and breakup due to shock wave and steady disturbances. Int. J. Multiph. Flow 21, 545–560. Hwang, S.S., Liu, Z., Reitz, R.D., 1996. Breakup mechanisms and drag coefficients of high-speed vaporizing liquid drops. Atomization Sprays 6, 353–376. Ibrahim, E.A., Yang, H.Q., Przekwas, A.J., 1993. Modeling of spray droplets deformation and breakup. J.Propuls.Power 9, 651–654. Krzeczkowski, S.A., 1980. Measurement of liquid droplet disintegration mechanisms. Int. J. Multiph. Flow 6, 227–239. Kulkarni, V., Sojka, P.E., 2014. Bag breakup of low viscosity drops in the presence of a continuous air jet. Phys. Fluids 26. Lee, M.W., Park, J.J., Farid, M.M., Yoon, S.S., 2012. Comparison and correction of the drop breakup models for stochastic dilute spray flow. Appl. Math. Modell. 36, 4512–4520. Liu, A., Mather, D., Reitz, R.D., 1993. Modeling the Effects of Drop Drag and Breakup on Fuel Sprays, 930072. In: International Congress and Exposition, Detroit, Michigan. SAE Technical Paper Series, Warrendale, PA, p. 930072. Matysiak, 2007. Euler–Lagrange verfahren zur simulation tropfenbeladener strömung in einem verdichtergitter. Helmut-Schmidt-Universität, Universität der Bundeswehr, Hamburg. Neupert, N., 2017. Experimentelle Untersuchung einer tropfenbeladenen Strömung in einer ebenen Verdichterkaskade. Helmut-Schmidt-Universität, Universität der Bundeswehr, Hamburg. Neupert, N., Gomaa, H., Joos, F., Weigand, B., 2016. Investigation and modeling of two phase flow through a compressor stage: analysis of film breakup. ISROMAC International Symposium on Transport Phenomena and Dynamics of Rotating Machinery. Neupert, N., Ober, B., Joos, F., 2014. Experimental investigation on droplet behaviour in a transonic compressor cascade, GT2014-25645. In: ASME Turbo Expo 2014: Turbine Technical Conference and Exposition, Düsseldorf, Germany, volume 2A of Turbomachinery, pp. 1–12. O’Rourke, P.J., Amsden, A.A., 1987. The Tab Method for Numerical Calculation of Spray Droplet Breakup, 872089. SAE International Fall Fuels and Lubricants Meeting and Exhibition. SAE Technical Paper Series, Warrandale, PA. Park, J.H., Yoon, Y., Hwang, S.S., 2002. Improved TAB model for prediction of spray droplet deformation and breakup. Atomization Sprays 12, 387–401. Park, S.W., Kim, S., Lee, C.S., 2006. Breakup and atomization characteristics of mono-dispersed diesel droplets in a cross-flow air stream. Int. J. Multiph. Flow 32, 807–822. Park, S.W., Lee, C.S., 2004. Investigation of atomization and evaporation characteristics of high-pressure injection diesel spray using Kelvin–Helmholz instability/droplet deformation and break-up competition model. Proc. Inst. Mech.Eng. 218, 767–777. Pilch, M., Erdman, C.A., 1987. Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop. Int. J. Multiph. Flow 13, 741–757. Schmehl, R., 2003. Tropfendeformation und Nachzerfall bei der technischen Gemischaufbereitung. university of Karlsruhe. Schmelz, F., 2002. Tropfenzerfall in beschleunigten Gasströmungen. Universität Dortmund, Dortmund. Storm, C., Joos, F., 2011. Euler–Lagrange method for numerical simulation of water droplet-laden compressor flow. In: ISABE 20th International Symposium on Air Breathing Engines, Gothenburg, volume 1. American Institute of Aeronautics and Astronautics, Sweden, pp. 662–672. Ulrichs, E., Joos, F., 2006. Experimental investigations of the influence of water droplets in compressor cascades, GT20 06-90411. In: ASME Turbo Expo 20 06 : Power for Land, Sea and Air, Barcelona, Spain, Vol. 6: Turbomachinery, Parts A and B, pp. 221–230. Wierzba, A., 1990. Deformation and breakup of liquid drops in a gas stream at nearly critical Weber numbers. Exp. Fluids 9, 59–64.