12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19–22, 2022

FRACTAL PROPERTIES OF ATOMIZING TURBULENT TWO-PHASE JETS

Josef Hasslberger, Elias Trautner & Markus Klein
Institute of Applied Mathematics and Scientific Computing
Department of Aerospace Engineering, University of the Bundeswehr Munich
Werner-Heisenberg-Weg 39, 85577 Neubiberg, Germany
josef.hasslberger@unibw.de

INTERFACES IN TURBULENT FLOWS

stretched, eventually rupture and form droplets. Secondary atomization, i.e. the breakup of droplets into smaller droplets, isnaturally included in the simulation but is not at the focus of
this work.
The rectangular computational domain extends 12D inaxial direction and 6D in both lateral directions, where D denotes the diameter of the round nozzle. Spatial discretization
by a uniform Cartesian mesh of grid size ∆ = D/64 yields an
overall number of 768×384×384 ≈ 1.13×108 computational
cells. At the side boundaries, homogeneous Neumann conditions are imposed, including a narrow band where the velocity is filtered to avoid numerical instabilities. At the outflow
boundary, back-flow into the domain is prevented by clipping
negative velocities to zero, also for stability reasons. To imitate
a turbulent nozzle flow at the inlet, a mean velocity profile issuperimposed by a homogeneous isotropic fluctuation level ofu′ /U0 = 5% and an integral turbulent length scale of Lt = D/4
based on the scaling and Gaussian filtering of random fluctuations (Klein et al., 2003). Corresponding to Diesel engine
conditions at a pressure of around 5.2 MPa and a temperatureof around 900 K, the density and dynamic viscosity ratio are
specified as ρl /ρg = 40 and µl /µg = 40, respectively. The liquid viscosity µl and the surface tension σ are adjusted according to the inflow-based Reynolds number Re = ρl U0 D/µl andWeber number We = ρl U02 D/σ in Table 1, respectively. The
same database has been used earlier to study the flow topology in both phases of atomizing liquid jets (Hasslberger et al.,2019). For one of the cases, i.e. Re = 5000 & We = 5000, Fig.
1 (top) provides an impression of the topologically complex
interface.
Based on the prescribed turbulence quantities at the inlet, relevant length scales characterizing this problem are presented in Table 1: Kolmogorov scale

Using a power law ansatz, fractal theory based sub-gridscale modeling is a popular approach in the Large Eddy Simulation (LES) of reactive flows, e.g. combustion. Due to theflame-interface analogy, a similar closure problem exists forinterface-dominated multi-phase flows, as further elaborated
by Klein et al. (2019); Hasslberger et al. (2020). The filteredsurface tension term in multi-phase LES reads σ⃗nκ|∇α| (using the surface tension coefficient σ , interface curvature κ, interface normal ⃗n = ∇α/|∇α|, volume fraction α) but cannot
be evaluated directly in a-posteriori LES.
According to literature, the flame iso-surface in the
flamelet regime (Kerstein, 1988) is characterized by a fractal
dimension D f of 7/3 = 2.333, while the maximum for a flamein the well-mixed regime (Hawkes et al., 2012) as well as apassive scalar iso-surface (Mandelbrot, 1975) has been foundto be 8/3 = 2.666. A theoretical maximum of D f = 3 is givenby the dimensionality of the confining three-dimensional Euclidean space.
Hence, it is the goal of this study to examine the fractal
properties of atomizing multi-phase jets and to indicate a potential modeling strategy for LES of such flows.

DIRECT NUMERICAL SIMULATION DATABASEThe open-source solver PARIS1 (Aniszewski et al., 2021)
is based on the one-fluid formulation of the unsteady incompressible Navier-Stokes equations including capillary forces.
Two immiscible fluids are represented by a jump in densityand viscosity. Propagation of the phase interface is implicitly
calculated by an advection equation for the cell-based volume
fraction α of one of the phases. State-of-the-art interface treatment is based on a geometrical Volume-of-Fluid method including piece-wise linear interface reconstruction and a heightfunction method for interface curvature calculation. The codeuses a second-order Runge-Kutta technique for time integration and second-order spatial discretization is realized by thefinite-volume approach on a regular, cubic staggered grid witha QUICK scheme for momentum advection and central differencing for diffusive fluxes. The projection method invokes
a successive over-relaxation solver for the pressure Poisson
equation.
In the simulation configuration, the liquid is injected
through a round nozzle into an ambient gas and the interface
immediately starts to wrinkle according to the turbulent inflow conditions. Interface corrugations grow in downstream
direction and ligaments are formed. These ligaments are being−3/4

lK ≈ Lt Ret= Lt

 ′ −3/4
u Lt
,
νl

(1)

Taylor microscale

lT ≈

√
−1/2
10Lt Ret(2)

and Hinze scalelH ≈ CH

1 http://www.ida.upmc.fr/∼zaleski/paris/

1

σ
ρl

3/5

ε −2/5(3)

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19–22, 2022

Figure 1. DNS database: Slices of the volume fraction contour for the Re = 5000 & We = 5000 case. Top: core+droplets; bottom:
jet core only.

volume averaging, indicated by ⟨·⟩, and taking the logarithm,
one obtains

with CH = 0.5 according to Deike et al. (2016). The turbulencedissipation is calculated as ε = k3/2 /Lt = (3/2)3/2 (u′ )3 /Lt .

Table 1. DNS database: Estimated Kolmogorov scalelK , Taylor microscale lT and Hinze scale lH , all nondimensionalized by nozzle diameter D, depending on the
Reynolds and Weber numbers.

Re-We

lK /D

lT /D

lH /D

2000-2000

0.022360.15810.0857

2000-5000

0.022360.15810.04955000-2000

0.011250.1

0.0857

5000-5000

0.011250.1

0.0495log

Following fractal theory, the wrinkling of the interface Ξ
can be expressed by the power law ansatz
|∇α|
=
|∇α|



ηo
ηi

D f −2,!
= (D f − 2) log(∆) − (D f − 2) log(ηi ),(5)which can be interpreted as a straight line equation with a slopeof D f − 2 in the double-logarithmic diagram. The intersectionof the straight line fit with the line given by Ξ = 1 yields theinner cut-off scale ηi . In agreement with expectations, Fig. 2
demonstrates that the self-similar fractal behavior can only beobserved for large normalized filter widths ∆/D. Accordingly,
the linear fit is based on the last few data points in this regime.
The procedure is described in more detail by Bambauer et al.
(2021), for instance.
As summarized in Table 2, the fractal dimension D f isclose to the asymptotic reference value of 8/3 = 2.666 in two
out of four cases. Smaller values of D f are likely due to insufficient atomization, hence wrinkling of the interface, under
the least violent conditions (Re = 2000 & We = 2000). Under
the most violent conditions considered here (Re = 5000 & We
= 5000), a significant number of disconnected droplets can beobserved. In order to isolate the wrinkling of the connected inner core region, a flood fill algorithm is applied to remove thedroplets and examine only the jet core in an additional artificial case. A visual comparison of those cases is shown in Fig.
1. Note that thin ligaments can be connected in the third dimension which is not evident from this two-dimensional slice.
Now, after removal of separated droplets, the passive scalarreference value of 8/3 is approached in all three cases undersufficiently violent conditions. This is also consistent withthe fact that increasing the Weber number decreases the multi-

A-PRIORI FRACTAL ANALYSIS

Ξ≡⟨|∇α|⟩⟨|∇α|⟩(4)

where the outer cut-off scale ηo is given by the filter width∆ in the context of explicit filtering (a-priori analysis). Upon2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19–22, 2022

Figure 2. A-priori analysis: Averaged interface wrinkling Ξ depending on the normalized filter width ∆/D for all four Re-We cases.
Last row: the Re = 5000 & We = 5000 case with droplets removed (i.e. jet core only). Blue marker symbols: measured values; redline: linear fit in the double-logarithmic diagram. Note the different scales of the ordinate axes.
Table 2. A-priori analysis: Fractal dimension D f and normalized inner cut-off scale ηi /D, depending on the Reynolds
and Weber numbers.

Re-We

Df

ηi /D

2000-2000

2.43030.1515

2000-5000

2.66260.0816

5000-2000

2.62380.0931

5000-5000

2.78980.05505000-5000 (core only)

2.69840.0564

A-POSTERIORI LES ANALYSISThe above analysis suggests a model for the total filteredsurface tension term, comprising the resolved and the sub-gridscale contribution, along the lines ofsσ⃗nκ|∇α| = σ⃗nκ |∇α| ≈ D f −2
 8/3−2ηo2∆σ⃗n κ|∇α|≈ σ⃗n κ|∇α|
,
ηi
ηi

(6)

where
surface-weighted
filteringisdefinedasss
φ = φ |∇α|/|∇α| and the surface-filtered product ⃗nκ
is approximated by the product of the corresponding resolved
quantities ⃗n and κ. Hence, the separate sub-grid scale
contribution reads"phase character of the flow with respect to the singular surface
tension force at the interface.
Comparing the measured inner cut-off scale ηi in Table 2
with the estimated characteristic length scales in Table 1, it canbe stated that the Taylor and Hinze scales are both on the sameorder as the inner cut-off scale. The fractal behavior might bedominated by the smaller or larger of both. The Kolmogorov
scale is about one order of magnitude smaller. In the context of breakup (and coalescence), the Taylor microscale andthe Hinze scale are also discussed in the literature (Saeedipour
et al., 2016; Ling et al., 2019; Rivière et al., 2021).

σ⃗n κ|∇α|2∆ηi

8/3−2

#
−1 .

(7)

For the outer cut-off scale, ηo = 2∆ is selected in accordancewith Nyquist’s criterion. For the inner cut-off scale, threeglobalglobalmodel variants are tested, namely ηi = lT
, ηi = lH
andlocal . While the global calculation of the length scalesηi = lH
introduced earlier (Eq. 2 and Eq. 3) is based on the prescribedturbulence characteristics at the inlet, the local calculation of
the Hinze scale accounts for local variations of the resolved

312th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19–22, 2022
and sub-grid scale dissipation, i.e.

consistently in the field statistics. Although both the global andlocal calculation of the Hinze scale lead to robust and reasonable results for this jet atomization case, it is expected that the
local calculation is superior in the general case. The phase inversion benchmark (Estivalezes et al., 2022) might be anotherinteresting test case for the fractal theory based modeling. Finally, it can be stated that the consequences of the model on theflow field can be global, meaning that they are visible in the velocity statistics even at the radial positions where ⟨α⟩ ≈ 1 or⟨α⟩ ≈ 0 in Fig. 5. Albeit the fractal model acts only at theinterface, which is most likely at the radial positions where⟨α⟩ ≈ 0.5, this behavior is explained by the active nature of
the volume fraction field feeding back to the momentum equation via the density and viscosity field as well as the surface
tension force. For example, a more compact liquid jet, i.e. lessfragmented, yields the non-local effect of decreased jet width
as measured by the velocity field.

ε = ε res + ε sgs
sgs= 2νSi j Si j + τi j Si j2
= 2νSi j Si j + 2CSmago
∆2

(8)
q
2Skl Skl Si j Si j ,where the eddy viscosity model by Smagorinsky (1963) has
been assumed.
In this a-posteriori analysis, the density and dynamic viscosity ratio are identical to the DNS database of the a-priori
analysis, but the jet Reynolds and Weber numbers are increased to somewhat more challenging conditions of Re =8000 & We = 5000. The computational domain is a rectangular box extending 16D in the axial (x) and 6D in the lateral(y, z) directions. An equidistant cubic grid with a resolution of∆ = D/128 is used for the reference DNS, leading to a total ofaround 1.2 × 109 cells, and a resolution of ∆ = D/16 is used
for all LES, leading to a total of around 2.4 × 106 cells. Insteadof including a nozzle geometry in the computational domain,a mean axial inflow velocity profile

CONCLUSIONS AND OUTLOOK
Using a DNS database of atomizing two-phase jets, the interface wrinkling and its fractal properties have been analyzed.
For sufficient atomization, i.e. sufficiently high Reynolds andWeber numbers, the fractal dimension is very similar to theasymptotic value of 8/3 corresponding to a passive scalar in
turbulent flows. The Taylor microscale and Hinze scale appear
as potential candidates for the inner cut-off scale.
Preliminary a-posteriori LES analysis revealed that allmodel candidates based on fractal theory show a nonnegligible effect on the statistical behaviour of velocity andvolume fraction fields. Both a global and local evaluation of
the inner cut-off scale result in stable simulations. Although
the first model tests are promising (especially in the far-fieldof the nozzle where sufficient interface fragmentation occurs),
a final assessment is yet to be presented because of the multiple
models and numerical aspects influencing the overall simulation result.
Future work will therefore be directed towards the interaction of sub-grid surface tension models with different numerical schemes and other recent sub-grid scale models forthe closure of momentum advection (Hasslberger et al., 2021).
In addition, it is worthwhile to examine a parameterization or
dynamic modeling of the fractal dimension.

q(y − yc )2 + (z − zc )2 − D2
U0 U0
 (9)
⟨ux ⟩ =−
· tanh 222θ

is prescribed (Stanley et al., 2002) at the round inlet area withdiameter D. Here, (yc , zc ) denotes the center of the nozzle,
and the momentum thickness is set to θ = D/20. This inlet
velocity profile is superposed by turbulent fluctuations withan integral turbulent length scale of Lt = D/8 and a uniform
isotropic fluctuation level of 7.5%. This configuration allows
for so-called entrainment, which describes the suction of air
into the domain across the lateral boundaries. It is noted thatthe (sub-grid) behavior of the jet core is of primary interest in
the present study since separated droplets, including secondary
atomization, are conveniently treated by Lagrangian particles
in practical LES. Since it is essential to minimize numericaldissipation in LES, momentum advection has been discretizedby central differences and the corresponding closure of subgrid scale stresses has been achieved by the state-of-the-art
eddy viscosity model of Nicoud et al. (2011).
The zeroth and first order statistics of the axial velocity
field, radial velocity field and volume fraction field are presented in Figs. 3, 4 and 5, respectively. Three axial distancesfrom the nozzle exit are considered, ranging from the nearfield (x/D = 4) to the far-field (x/D = 12), to obtain a complete
picture of the model performance. For sufficient distance fromthe nozzle exit, i.e. x/D = 8 and x/D = 12, the sub-grid surface
tension models lead to an improvement with respect to mean
fields and the corresponding fluctuation statistics. The several assumptions in Eq. 6, including also the fractal dimension
value of D f = 8/3, seem to be better justified in this region.
The situation in the near-field (x/D = 4) is less unambiguous. Some of the modeling assumptions (fully developed turbulence, sufficient interface fragmentation) are better justified
in the far-field than in the near-field. It is also unclear how well
the dominant shear layers close to the nozzle exit can be represented on the LES mesh and what the relative effects of the
sub-grid stress and sub-grid surface tension models are. In thisa-posteriori configuration, the Taylor microscale model variantshows a weaker influence than the Hinze scale model variant
globalglobal/D = 0.044 > lH
/D = 0.022 and this appears
since lT

ACKNOWLEDGEMENTS
Support by the German Research Foundation (Deutsche
Forschungsgemeinschaft - DFG, GS: KL1456/1-3) is gratefully acknowledged.

REFERENCES
Aniszewski, W., Arrufat, T., Crialesi-Esposito, M., Dabiri, S.,
Fuster, D., Ling, Y., Lu, J., Malan, L., Pal, S., Scardovelli,
R. et al. 2021 Parallel, robust, interface simulator (PARIS).
Computer Physics Communications 263, 107849.
Bambauer, M., Chakraborty, N., Klein, M. & Hasslberger, J.
2021 Vortex dynamics and fractal structures in reactive andnonreactive Richtmyer-Meshkov instability. Physics of Fluids 33, 044114.
Deike, L., Melville, W.K. & Popinet, S. 2016 Air entrainmentand bubble statistics in breaking waves. Journal of Fluid
Mechanics 801, 91–129.
Estivalezes, J.-L., Aniszewski, W., Auguste, F., Ling, Y., Osmar, L., Caltagirone, J.-P., Chirco, L., Pedrono, A., Popinet,4

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19–22, 2022

1/2 for varying
Figure 3. A-posteriori analysis: Radial profiles of the mean axial velocity ⟨ux ⟩ and r.m.s. velocity fluctuation (⟨u′2
x ⟩)axial distance x/D from the nozzle exit.

1/2 for varying
Figure 4. A-posteriori analysis: Radial profiles of the mean radial velocity ⟨ur ⟩ and r.m.s. velocity fluctuation (⟨u′2
r ⟩)axial distance x/D from the nozzle exit.

S., Berlemont, A. et al. 2022 A phase inversion benchmark
for multiscale multiphase flows. Journal of Computational
Physics 450, 110810.
Hasslberger, J., Engelmann, L., Kempf, A. & Klein, M.
2021 Robust dynamic adaptation of the Smagorinsky model
based on a sub-grid activity sensor. Physics of Fluids 33 (1),015117.
Hasslberger, J., Ketterl, S. & Klein, M. 2020 A-priori assessment of interfacial sub-grid scale closures in the twophase flow LES context. Flow, Turbulence and Combustion
105 (2), 359–375.
Hasslberger, J., Ketterl, S., Klein, M. & Chakraborty, N. 2019
Flow topologies in primary atomization of liquid jets: adirect numerical simulation analysis. Journal of Fluid Mechanics 859, 819–838.
Hawkes, E.R., Chatakonda, O., Kolla, H., Kerstein, A.R. &Chen, J.H. 2012 A petascale direct numerical simulation
study of the modelling of flame wrinkling for large-eddy

simulations in intense turbulence. Combustion and Flame
159 (8), 2690–2703.
Kerstein, A.R. 1988 Fractal dimension of turbulent premixed
flames. Combustion Science and Technology 60 (4-6), 441–
445.
Klein, M., Ketterl, S. & Hasslberger, J. 2019 Large eddy simulation of multiphase flows using the volume of fluid method:
Part 1 – governing equations and a-priori analysis. Experimental and Computational Multiphase Flow 1 (2), 130–
144.
Klein, M., Sadiki, A. & Janicka, J. 2003 A digital filter based
generation of inflow data for spatially developing direct
numerical or large eddy simulations. Journal of Computational Physics 186 (2), 652–665.
Ling, Y., Fuster, D., Tryggvason, G. & Zaleski, S. 2019 A twophase mixing layer between parallel gas and liquid streams:
multiphase turbulence statistics and influence of interfacial
instability. Journal of Fluid Mechanics 859, 268–307.

512th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19–22, 2022

Figure 5. A-posteriori analysis: Radial profiles of the mean volume fraction ⟨α⟩ and r.m.s. volume fraction fluctuation (⟨α ′2 ⟩)1/2 forvarying axial distance x/D from the nozzle exit.

Mandelbrot, B.B. 1975 On the geometry of homogeneous turbulence, with stress on the fractal dimension of the isosurfaces of scalars. Journal of Fluid Mechanics 72 (3), 401–
416.
Nicoud, F., Toda, H.B., Cabrit, O., Bose, S. & Lee, J. 2011 Using singular values to build a subgrid-scale model for large
eddy simulations. Physics of Fluids 23 (8), 085106.
Rivière, A., Mostert, W., Perrard, S. & Deike, L. 2021 Subhinze scale bubble production in turbulent bubble break-up.
Journal of Fluid Mechanics 917.
Saeedipour, M., Pirker, S., Bozorgi, S. & Schneiderbauer, S.

2016 An Eulerian–Lagrangian hybrid model for the coarsegrid simulation of turbulent liquid jet breakup. International
Journal of Multiphase Flow 82, 17–26.
Smagorinsky, J. 1963 General circulation experiments withthe primitive equations: I. the basic experiment. Monthly
Weather Review 91 (3), 99–164.
Stanley, S.A., Sarkar, S. & Mellado, J.P. 2002 A study of theflow-field evolution and mixing in a planar turbulent jet using direct numerical simulation. Journal of Fluid Mechanics
450, 377–407.

6