# Mass entrainment rate of an ideal momentum turbulent round jet

Fermín Franco Medrano1,2, Yasuhide Fukumoto1, Clara M. Velte3
and Azur Hodžić3
1 Institute of Mathematics for Industry, Kyushu University
West 1 Bldg., 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japan
2 Research Institute of Composite Materials, Fukuoka University
Bldg. 6, 8-19-1 Nanakuma, Jonan-ku, Fukuoka 814-0180, Japan
3 Department of Mechanical Engineering, Technical University of Denmark
Bldg. 404, Nils Koppels Allé, Kgs. Lyngby DK-2800, Denmark Corresponding
author<EMAIL_ADDRESS>
###### Abstract

We propose a two-phase-fluid model for a full-cone turbulent round jet that
describes its dynamics in a simple but comprehensive manner with only the apex
angle of the cone being a disposable parameter. The basic assumptions are that
(i) the jet is statistically stationary and that (ii) it can be approximated
by a mixture of two fluids with their phases in dynamic equilibrium. To derive
the model, we impose conservation of the initial volume and total momentum
fluxes. Our model equations admit analytical solutions for the composite
density and velocity of the two-phase fluid, both as functions of the distance
from the nozzle, from which the dynamic pressure and the mass entrainment rate
are calculated. Assuming a far-field approximation, we theoretically derive a
constant entrainment rate coefficient solely in terms of the cone angle.
Moreover, we carry out experiments for a single-phase turbulent air jet and
show that the predictions of our model compare well with this and other
experimental data of atomizing liquid jets.

##  1 Introduction

Liquid jets appear in a vast range of applications. A very active field of
application is that of fuel jet injection engines, widely used in the
automotive and aerospace industry; other fields of application include medical
apparatuses, so-called “atomizers” used by commercial products in many
industries, flows through hoses and nozzles for various industrial purposes as
well as firefighting. Some of these applications and atomization methods are
described in Jiang et al. (2010).

After the pioneering theoretical work of Rayleigh (1878) many other early
works advanced and improved upon this subject (Weber, 1931; Tomotika, 1935;
Taylor, 1962) (see Gorokhovski and Herrmann, 2008). However, these results
focused on what is termed “primary breakup”, i.e. the transition of the jet
from a cylindrical geometry to the formation of the first detached droplets
and “ligaments”. This breakup has been found to depend on numerous parameters,
and a characterization of the breakup mechanism, based on the dominant
physical forces acting on it has been achieved with some success. This was
summarized in what is called the “breakup regimes” and it can be described
broadly in terms of two parameters: the jet’s Reynolds and the Weber numbers.
Another classification which separated these breakup regimes was based on the
jet’s speed and the “Z length” or the length of the “continuous” part of the
jet. Four main breakup regimes have been identified: (i) the capillary or
Rayleigh regime; (ii) the first wind-induced regime; (iii) the second wind-
induced regime; and (iv) the atomization regime. Some review articles on the
topic include McCarthy and Molloy (1974); Lin and Reitz (1998); Liu (2000);
Birouk and Lekic (2009); Jiang et al. (2010).

The range of applications involving atomizing liquid jets forming two-phase
fluid flows is still large. The complexity of the atomizing process, involving
numerous physical phenomena and many variables, ranging from the conditions
inside the nozzle (or some generating source) to the interaction between the
atomization process and the environment into which the jet is penetrating, all
account for numerous challenges in physical and mathematical modeling.
Notwithstanding, several mathematical models have been attempted to describe
different aspects of the jets in this regime. For example, differential
equations for a fuel jet’s tip penetration distance as a function of time
(Wakuri et al., 1960; Sazhin et al., 2001; Desantes et al., 2005; Pastor et
al., 2008); models for the gas entrainment rate in a full-cone spray (Cossali,
2001); and a one-dimensional model for the induced air velocity in sprays
(Ghosh and Hunt, 1994). None of these models is sufficient by itself as
explained below.

In this study we propose an original 1D mathematical model for the macroscopic
dynamics of a full-cone turbulent round jet ensuing from a circular nozzle
into an stagnant fluid. This kind of jet serves as a basis for many industrial
processes in modern manufacturing industry (Jiang et al., 2010). An advantage
of our model over other analytical 1D models is that this model has a single
experimentally measurable parameter while it maintains reasonable predictive
power and gives theoretical understanding that allows it to analytically
calculate other physical quantities of interest. Moreover, this model can be
extended to accommodate an energy conservation approach with simple turbulence
and energy dissipation models, resulting in an increased accuracy as will be
reported in a subsequent paper.

Figure 1: Diagram of the mass entrainment process and the relevant physical
variables.

Our particular concern lies in the mass entrainment rate, $\hat{m}_{e}$, a
dimensionless quantity that describes the proportion of entrained mass in a
jet. The mass entrainment process is sketched in Fig. 1. This quantity is
important in applications like diesel engine combustion, as it is related to
the rate of mixing of fuel and air (Post et al., 2000). The entrainment rate
has been extensively measured for broad experimental settings, as defined from
early work on single-phase gas jets (Ricou and Spalding, 1961; Hill, 1972) and
later for two-phase liquid-gas jets (Ruff et al., 1989; Hiroyasu and Arai,
1990; Hosoya and Obokata, 1993; Cossali et al., 1996). Formulas derived from
data fitting or dimensional analysis are available in the aforementioned
works. A more recent definition of the mass entrainment rate is given by Post
et al. (2000) as

$\hat{m}_{e}=\frac{\dot{m}_{e}(x)}{\dot{m}_{0}}=\frac{\dot{m}(x)-\dot{m}_{0}}{\dot{m}_{0}}=K_{e}\left(\frac{x}{D_{0}}\right)\left(\frac{\rho_{e}}{\rho_{0}}\right)^{1/2},$
(1)

where $\dot{m}_{e}(x)=\dot{m}(x)-\dot{m}_{0}$, and $\dot{m}_{0}$ and
$\dot{m}(x)$ are respectively the initial and total mass flux rates, defined
as the fluid mass fluxes per unit time through a cross-section at respectively
axial distance 0 and $x$ from the nozzle. Here $\rho_{0}$ and $\rho_{e}$ are
the initial and entrained fluid densities, respectively, $D_{0}$ is the nozzle
diameter and $x$ is the distance from the nozzle’s exit. The dimensionless
number $K_{e}$ is called the entrainment rate coefficient, which taking the
derivative of Eq. (1) can thus be calculated as (Post et al., 2000)

$K_{e}(x)=\frac{d\dot{m_{e}}}{dx}\left(\frac{D_{0}}{\dot{m}_{0}}\right)\left(\frac{\rho_{0}}{\rho_{e}}\right)^{1/2}=\frac{d\hat{m_{e}}}{d\hat{x}}\rho_{*}^{-1/2};$
(2)

where the last equation uses dimensionless units as defined in Sects. 3 and 4.
The definition of Eq. (1) for gaseous jets was motivated by experimental
evidence showing that $K_{e}$ tends to a constant in the far field (Ricou and
Spalding, 1961; Hill, 1972; Abraham, 1996; Post et al., 2000), but it can be
applied to a non-constant entrainment rate coefficient using Eq. (2),
particularly in the near and intermediate fields. Despite being originally
formulated for gaseous jets, the definition in Eq. (2) is also used for the
entrainment rate coefficient of (two-phase) atomized liquid jets (Ruff et al.,
1989; Cossali et al., 1996; Cossali, 2001; Hosoya and Obokata, 1993; Song and
Abraham, 2003; Rabadi et al., 2007). To the the best of the authors’ knowledge
however, there is so far no theoretical derivation of it from first principles
that is of practical use; although there are some models available (Cossali,
2001; Rabadi et al., 2007), the most relevant of which are discussed below.
One of the goals of this paper is to present this derivation. Our obtained
formula generalizes Eq. (1) and reduces to a constant $K_{e}$ at large
distances from the nozzle (far field). An approximation for the near field is
also given.

In the context of a fuel injection engine, Sazhin et al. (2001) derived a
specialized model for the fuel jet’s tip penetration distance as a function of
time. Their model, however, relied on the unrealistic assumption that the mass
ratio of droplets-to-gas ($\alpha_{d}$, in their notation) is constant
throughout the jet’s length, which they then used in order to be able to
integrate the resulting 1D ODE. In the same context of a combustion engine and
fuel injection, Desantes et al. (2005) proposed another 1D model also for the
diesel spray tip penetration. Their model relied only on conservation of
momentum and on fixing the jet’s geometry by means of measuring the cone
angle. They replaced the dynamic analysis of a spray by the analysis of an
assumed analogously-behaved incompressible gas jet (i.e. ignoring the
droplets), by exchanging the original nozzle’s exit diameter for an
“equivalent diameter” and setting the density of the jet equal to that of the
gas. Pastor et al. (2008) also proposed a related model for the more general
case of the transient (time dependent) evolution of the tip penetration
variable. They based their calculations on mass, momentum and enthalpy
conservation equations. Their results also present in an appendix for the
steady state (time independent) as a special case. However, even in this
special case, their solutions are always numerical, except in a constant-
density flow case discussed in the appendix. Cossali (2001) proposed a model
for the gas entrainment in a full cone spray. The model contained newly
proposed physical parameters that are not easy to measure experimentally (like
three different mean droplet sizes, integrals over the radial distributions of
gas and droplets concentrations, velocity, etc.) so it is unclear how the
overall theory may produce any prediction from first principles. Thus,
compared to the present study, previous models provide a less comprehensive
description of the density or liquid fraction of the spray, make unrealistic
assumptions or introduce parameters unavailable experimentally.

There are also widely used numerical models based on turbulence modeling for
jets, like the one by Vallet et al. (2001) used in CFD programs like Star-CD,
KIVA-3V, Ansys Fluent, as well as open source codes like OpenFOAM (Stevenin et
al., 2016). The present work is concerned with developing an analytical model
with the advantage of producing closed expressions.

The article is organized as follows: Sect. 2 describes the experimental
setting. In Sect. 3 we derive from physical conservation laws the basic
mathematical model for the ideal momentum atomizing liquid jet. In Sect. 4 we
present the theoretical derivation of the entrainment rate coefficient, based
on the derived model, and in Sect. 5 we apply the model to some common
limiting and special cases that recover the known solutions (viz. very thin
ambient fluid, (single-phase) submerged jet, and cylindrical jet without
breakup). Next, in Sect. 6 we compare the predictions with laboratory
experiments from the literature and our own from the Technical University of
Denmark. Finally, in Sect. 7 we summarize the main conclusions of the present
work.

##  2 Experimental method

The experiments for a single-phase turbulent air jet were carried out at the
Department of Mechanical Engineering, Technical University of Denmark.
Although this is not an atomizing liquid jet, the special case equations of
Sect. 5.4 can be used to test the present theory for this case and general
insight into how to determine the jet’s geometry is obtained (needed for the
present model).

### Jet facility and enclosure

The jet was the same as that used by Ewing et al. (2007) and the enclosure was
similar as well. The jet was a cubic box of dimensions $58\times 58.5\times
59\,cm^{3}$ fitted with an axisymmetric plexiglas nozzle, tooled into a fifth-
order polynomial contraction from an interior diameter of $32\,mm$ to an exit
diameter of $D_{0}=10\,mm$. The interior of the box was stacked with foam
baffles in order to damp out disturbances from the fan that supplied the
generator with pressurized air. For further details on the generator box, see
Jung et al. (2004); Gamard et al. (2004).

The flow generating box rested rigidly on an aluminum frame. The exit velocity
was monitored via a pressure tap in the nozzle positioned upstream of the
contraction and connected to a digital manometer by a silicon tube. The
ambient pressure was monitored by an independent barometer. The air intake was
located inside the jet enclosure so that both jet and ambient are seeded to
provide as homogeneous a seeding distribution in the measured flow as
possible. The seeding particles were generated from liquid Glycerin using an
atomizing nozzle built in-house, providing tracers with a typical size of
about $2-3\,\mu m$. The enclosure consisted of a large tent of dimension
$2.5\times 3.0\times 4.5\,m^{3}$, yielding a cross-sectional area of
$2.5\times 3.0\,m^{2}$ that resulted in a momentum conservation of $100.0\%$
at $x/D_{0}=30$ and $98.4\%$ at $x/D_{0}=100$, where $x$ is the downstream
distance from the jet exit (Hussein et al., 1994).

Depending on ambient conditions, such as barometric pressure and ambient
temperature, a pressure drop of about $570\,Pa$ was measured across the jet
exit nozzle corresponding to a target jet exit velocity of $30\,m/s$ and
Re=20,000.

### Measurement equipment

The ensemble averaged flow field was measured using planar 2-component
Particle Image Velocimetry (PIV). The scatterers were illuminated using a
DualPower $200-15$ $200\,mJ$ double cavity Nd:YAG $532\,nm$ pulsed PIV laser
with maximum power output of $1200\,mJ$ and a pulse duration of $4\,ns$. In
order to produce a light sheet following the downstream development of the
jet, the laser was positioned below the jet nozzle where an adjustable optical
mirror was placed.

Figure 2: Diagram of the experimental setup and measuring devices.

Acquisition of the flow field was carried out in four different regions to
utilize the spatial resolution capacity of the cameras efficiently. The near
field ($x/D_{0}=0.2-6.8$) and the intermediate field ($x/D_{0}=6.8-28$) were
acquired independently using a HiSense MkII camera with a $1344\times 1024$
pixels CCD sensor and a $60\,mm$ focal length lens. The aperture was set to
$f^{\sharp}2.8$. The far field ($x/D_{0}=27.4-68.6$ and $66.4-106.3$,
respectively) was acquired using two Dantec FlowSenseEO $16\,MPix$ cameras
($4872\times 3248\,Pix$ with a pixel pitch of $7.4\,\mu m$) with $60\,mm$
Nikon lenses with an aperture of $f^{\sharp}2.8$. The different regions
sampled and their respective final Interrogation Area (IA) side lengths are
listed in Table 1. A detailed diagram of the whole setup is shown in Fig. 2.

Table 1: Acquisition parameters for the measured regions. Measured region | extent, $\hat{x}$ | side length, mm
---|---|---
Near field | $1.2-6.8$ | 0.80
Intermediate fld. | $6.8-28.0$ | 2.59
Far field 1 | $27.4-68.6$ | 1.69
Far field 2 | $66.4-106.3$ | 2.66

### Data processing

The images were acquired and processed using Dynamic Studio 2015a (4.15.115).
A correlation scheme computing $32\times 32$ pixels interrogation areas with
$50\%$ overlap in both directions using 2 refinement steps, deforming windows
and subpixel interpolation for enhanced accuracy. Local median validation with
a $3\times 3$ vector neighborhood validation base was implemented to discard
invalid vectors.

Preliminary measurements in the most downstream far field measurement provided
an estimate of the integral time scale of $0.059\,s$ from the local jet-half
width velocity. In order to ensure statistical independence of the acquired
samples, a sampling rate of $1\,Hz$ was chosen. In the near and intermediate
fields, $2000$ realizations of the flow were collected to ensure statistical
convergence of the mean velocity. Averaging over $500$, $1000$ and $1500$
samples did not display significant deviations from the average from $2000$
realizations and the statistics can therefore be considered converged to
sufficient degree for the current purposes. For the far field, 11,000
realizations were collected.

In relevant cases, to remove faint reflections and other adverse effects of
the background, the images were first preprocessed by computing, for each
pixel separately, the minimum intensity of the ensemble of measurement images
and subtracting this from the same ensemble. The average velocity field was
then obtained by only including valid vectors, discarding vectors substituted
by the validation algorithm of the software.

##  3 The ideal momentum atomizing liquid jet model

Consider a stationary full-cone turbulent round jet ensuing from a circular
nozzle (of small diameter) into an ambient fluid, with a constant high gauge
pressure input fluid and a small conical jet angle. We want to calculate the
dynamical variables of the two-phase jet at some distance from the nozzle. The
relevant variables and parameters of the model are depicted in Fig. 1, where
$\theta$ is the spread half-angle of the conical jet, $\rho$ the density of
the fluid, $D$ the diameter of jet, $p$ the fluid’s pressure, $v$ the axial
velocity of a fluid element averaged over a cross section of constant $x$,
where $x$ is the axial distance from the nozzle. The subscript “0” indicates
initial values (at the nozzle’s exit position), e.g. $D_{0}$ is the diameter
of the jet at the nozzle’s exit (equal to the nozzle’s orifice diameter).

###  3.1 Initial momentum flux

The momentum of a flat disc of fluid of infinitesimal width exiting the
circular nozzle is

$d\Pi_{0}=m_{0}v_{0},$ (3)

where the mass of the flat disc is $m_{0}=\frac{1}{4}\rho_{0}\pi
D_{0}^{2}\,dx$. Note that $dx=v_{0}\,dt$. This velocity may be calculated from
the input gauge pressure, $p_{0}$, by Bernoulli’s theorem, neglecting the
dynamic pressure inside the nozzle and assuming the static pressure is totally
converted to the jet’s dynamic pressure just outside the nozzle. Thus
$p_{0}\approx\frac{1}{2}\rho_{0}v_{0}^{2}$ implies

$v_{0}=\sqrt{\frac{2p_{0}}{\rho_{0}}}.$ (4)

Substituting $m_{0}$ and $v_{0}$ into equation (3) we get

$\dot{\Pi}_{0}:=\left.\frac{d\Pi}{dt}\right|_{x=0}=\frac{1}{2}\pi
D_{0}^{2}p_{0},$ (5)

which is the momentum flux per unit time coming out of the nozzle as a result
of the input gauge pressure inside the nozzle.

###  3.2 Conservation of momentum

We assume that the fluid at a distance $x$ is a two-phase fluid. The present
model has been termed “ideal momentum” to distinguish it from “ideal energy”
and “lossy” models including momentum, energy or mass loss parameters,
developed by the present authors and to be reported in detail in a subsequent
article. Accordingly, first assume conservation of momentum, i.e. the momentum
of the two-phase fluid is solely that from the original input pressure. The
droplets transfer momentum to the initially static fluid by drag forces
(Sazhin, 2006; Fuchimoto et al., 2009) and they reach local dynamic
equilibrium in such a way that both phases move at the same speed $v$ inside
the jet (Desantes et al., 2011); this is also the main assumption under the
wide class of “Locally Homogeneous Flows” (LHF) (Faeth, 1983, 1987; Faeth et
al., 1995). We assume that the latter process occurs so fast immediately
outside the nozzle’s exit that we may neglect the non-equilibrium zone near
the nozzle. The latter assumption is reasonable for high-speed pressure
atomized jets, e.g. like the ones used in real life diesel engines (Sazhin et
al., 2001). This effectively allows us to treat the two-phase flow as a single
fluid with a composite density $\rho(x)$ and a single velocity $v(x)$,
depending on the distance from the nozzle, $x$.

Analogous to the calculations in Sect. 3.1, the momentum at some distance $x$
is $d\Pi=mv$, where $m=\frac{1}{4}\rho\pi D^{2}\,dx$ and $dx=v\,dt$, by which

$\frac{d\Pi(x)}{dt}=\frac{1}{4}\pi\rho D^{2}v^{2},$ (6)

which is the momentum flux at a cross-section of the two-phase fluid jet at a
distance $x$ from the nozzle.

Equating Eq. (5) and Eq. (6) by conservation of the momentum contained between
the two planes $x=0$ and at axial distance $x$ in a steady state, using again
Eq. (4) and solving for $v$ we get

$v^{2}=\left(\frac{D_{0}}{D}\right)^{2}\frac{\rho_{0}}{\rho}v_{0}^{2}.$ (7)

This can be written in simpler form by using dimensionless units as
$\hat{v}^{2}=\hat{\rho}^{-1}\hat{D}^{-2}$, where $v$, $D$ and $\rho$ have been
scaled respectively by $v_{0}$, $D_{0}$ and $\rho_{0}$ (the quantities at the
nozzle’s exit) and denoted by a hat symbol.

###  3.3 Volume of entrained mass in the two-phase fluid

The total fluid volume of the thin disc at the target distance is
$dV=dV_{0}+dV_{e}$, where the subscript “e” denotes the quantities related to
the entrained fluid, i.e. the total volume $dV$ of the two-phase fluid is just
the sum of the volume of the original quantity of fluid coming out of the
nozzle $dV_{0}$ plus the added volume of entrained fluid in dynamic
equilibrium, $dV_{e}$. Also, the total volume of the jet at the target
distance is straightforward to calculate from the conical geometry. Then
$dV_{e}=dV-dV_{0}$ from where we can calculate a volumetric flow rate of the
mass entrainment:

$\frac{dV_{e}}{dt}=\frac{1}{4}\pi(D^{2}v-D_{0}^{2}v_{0})$ (8)

Here we can substitute for $v_{0}$ from equation (4) and
$D=D_{0}+2x\tan(\theta)$.

###  3.4 Density of the two-phase fluid

The mean composite density of a two-phase fluid thin-disc element is just the
total mass over the total volume:

$\rho=\frac{dm}{dV}=\frac{dm_{0}+dm_{e}}{dV_{0}+dV_{e}}=\frac{\rho_{0}dV_{0}+\rho_{e}dV_{e}}{dV_{0}+dV_{e}}.$
(9)

After substituting $dV_{e}$ from Eq. (8), and the initial volume element
$dV_{0}$ from analogous calculations, simplifying and solving for $\rho$, we
obtain

$\rho=\rho_{e}+\frac{D_{0}^{2}v_{0}}{D^{2}v}(\rho_{0}-\rho_{e})\quad\textrm{or}\quad\hat{\rho}=\rho_{*}+(1-\rho_{*})/\hat{D}^{2}\hat{v},$
(10)

the density of the two-phase fluid at the distance $x$ from the nozzle.
Fortunately, it depends on $v$ which makes the dependency closed as we can see
from Eq. (7). The second equation in Eq. (10) is the dimensionless form, where
$\rho_{*}=\rho_{e}/\rho_{0}$.

We implicitly assume, by calculating the _mean composite density_ of the two-
phase fluid in Eq. (9), that the droplets distribution throughout the disc
two-phase fluid element does not differ greatly from a uniform distribution.
Notice that we approximate the front of the jet by a planar front of equal
density, i.e. a “top-hat” radial distribution. In reality this is not true,
since the front should be spherical in the first order, and it is then in a
spherical shell within the jet’s cone that we should consider $\rho$ to be
approximately constant, not in a plane. However, for small half-angles
$\theta$ and short distances $x$ a plane should be suffice as a first order
approximation. The same could be said of the front’s velocity $v$. Overall, we
may take the above considerations as utilizing “top-hat” velocity and density
distributions as a first approximation. Note that slicing the spherical jet
front with a $x$-normal plane provides a non-constant $\rho$ density
distribution in this plane. This distribution should, however, be similar to a
two-dimensional Gaussian distribution centered around the $x$-axis, i.e. the
jet’s “centerline”. There are some models (Ghosh and Hunt, 1994; Desantes et
al., 2005; Rabadi et al., 2007; Pastor et al., 2008; Desantes et al., 2011)
which apply self-similar Gaussian velocity distributions as initial
assumptions; however, this calculation will be included in a later work since
we anticipate that it would not lead to a major refinement of the axial
centerline quantities.

###  3.5 Explicit expressions for the axial velocity and composite density of
the spray

From Eqs. (7) and (10) we identify a system of two nonlinear equations with
two unknowns, $v$ and $\rho$. We can eliminate $\rho$ from the system to get
$v$ explicitly in dimensionless form as

$\hat{v}=\frac{\rho_{*}-1+\sqrt{(\rho_{*}-1)^{2}+4\hat{D}^{2}\rho_{*}}}{2\hat{D}^{2}\rho_{*}}.$
(11)

Eq. (11) is equivalent to the one derived independently by Wakuri et al.
(1960) using a similar approach. Analogously, we obtain an explicit
dimensionless form for the density, $\hat{\rho}=\rho/\rho_{0}$, by eliminating
$v$ from the same described system, obtaining:

$\hat{\rho}=\tilde{\rho}+\sqrt{\tilde{\rho}^{2}-\rho_{*}^{2}},$ (12)

where $\tilde{\rho}=\rho_{*}+(1-\rho_{*})^{2}/2\hat{D}^{2}$. From Eqs. (11)
and (12), the dimensionless form of the dynamic pressure,
$\hat{p}=p/p_{0}=\hat{\rho}\hat{v}^{2}$ (where
$p_{0}=\frac{1}{2}\rho_{0}v_{0}^{2}$ is the output dynamic pressure at the
nozzle’s exit), which accounts for the total pressure at the target distance,
may be calculated as

$\hat{p}=\left(\frac{\hat{\rho}}{\hat{D}}\right)^{2}.$ (13)

Notice the implicit single-variable dependence of $\hat{v}(\hat{x})$,
$\hat{\rho}(\hat{x})$ and $\hat{p}(\hat{x})$ on the dimensionless axial
distance $\hat{x}=x/D_{0}$ by substituting $\hat{D}=1+2\hat{x}\tan\theta$ in
Eqs. (11)-(13).

##  4 Mass entrainment rate

We have the volumetric entrainment rate $dV_{e}/dt$ as given by Eq. (8); which
multiplying by $\rho_{e}$ and normalizing by the volumetric flux at the
nozzle, $\dot{m}_{0}=\frac{1}{4}\pi\rho_{0}D_{0}^{2}v_{0}$, we get the mass
entrainment rate. Using the dimensionless notation introduced above, we can
write

$\hat{m}_{e}=\rho_{*}(\hat{D}^{2}\hat{v}-1).$ (14)

as defined in Eq. (1). Notice $\frac{d\hat{V}_{g}}{dt}=\hat{D}^{2}\hat{v}-1$,
by which (14) can be written simply as
$\hat{m}_{e}=\rho_{*}\frac{d\hat{V}_{e}}{dt}$. Substituting $\hat{v}$ from Eq.
(11) we have

$\hat{m_{e}}(\hat{x})=\frac{1}{2}\left[\sqrt{(\rho_{*}-1)^{2}+4\rho_{*}\hat{D}(\hat{x})}-\rho_{*}-1\right].$
(15)

Equation (15) gives us an explicit relationship between the mass entrainment
rate $\hat{m_{e}}$ and the axial distance, $x$, remembering that
$\hat{D}(\hat{x})$ depends solely on this distance and the jet’s spread angle,
$\theta$. It is thus a theoretical derivation of Eq. (1) in a generalized form
and from first principles. It is easy to check that when $x\rightarrow 0$,
i.e. at the nozzle’s exit, $\hat{m}_{e}\rightarrow 0$, as it should be since
there is no entrained gas there yet.

Taking the derivative of Eq. (14) with respect to $\hat{x}$ and substituting
$\hat{v}$ from (11), we get

$\frac{d\hat{m_{e}}}{d\hat{x}}=\frac{4\rho_{*}\hat{D}\tan\theta}{\sqrt{(\rho_{*}-1)^{2}+4\rho_{*}\hat{D}^{2}}}.$
(16)

We can transform Eq. (16) into its dimensional form and substitute it into Eq.
(2) to obtain a new explicit formula for $K_{e}$, which can then be written as

$K_{e}(x)=\frac{4\rho_{*}^{1/2}D(x)\tan\theta}{\sqrt{(\rho_{*}-1)^{2}D_{0}^{2}+4\rho_{*}D(x)^{2}}}.$
(17)

Equation (17) is a new general explicit expression for the mass entrainment
rate coefficient as a function of the axial distance from the nozzle. Notably,
$K_{e}(x)$ includes the half angle $\theta$ of the cone apex as an adjustable
parameter. In Sect. 5.1 below we will theoretically deduce that our formula in
Eq. (17) can be approximated by a constant for the far field, as proposed by
other authors (Ricou and Spalding, 1961; Hill, 1972; Abraham, 1996; Post et
al., 2000) based on experimental results. We also will give a new constant
value for the near field.

##  5 Limiting and special cases

Let us now calculate some special and limit cases of the ideal momentum
atomizing liquid jet model.

###  5.1 The near and far fields

Very near the nozzle, taking $\hat{x}\rightarrow 0$ in Eq. (11) we get
$\hat{D}\rightarrow 1$ and $\hat{v}\rightarrow 1$, correspondingly
$v\rightarrow v_{0}$ in dimensional form, as expected since $v_{0}$ is the
velocity of the liquid exiting the nozzle. From Eq. (12) we get
$\hat{\rho}\rightarrow 1$ thus $\rho\rightarrow\rho_{0}$, also as expected. We
can also check from Eq. (13) that $\hat{p}\rightarrow 1$, thus $p\rightarrow
p_{0}$. All the limiting values agree with the ones inside the nozzle. As for
the entrainment rate, taking the limit of Eq. (17) for $x\rightarrow 0$ we
have the near field entrainment rate constant

$K_{e,\,near}=\frac{4\rho_{*}^{1/2}\tan\theta}{\rho_{*}+1},$ (18)

because $D(x)\rightarrow D_{0}$ as $x\rightarrow 0$. On the other hand, in the
far limit, taking $\hat{x}\rightarrow+\infty$ in (11) we get
$\hat{v}\rightarrow 0$ which also corresponds to $v\rightarrow 0$ in
dimensional form, as expected since the ideal jet expands infinitely,
$\hat{D}\rightarrow+\infty$, and thus eventually tends to diffuse all its
kinetic energy to the increasing amounts of entrained gas. From Eq. (12) we
get $\hat{\rho}\rightarrow\rho_{*}$ and thus $\rho\rightarrow\rho_{e}$, also
as expected, since the initial amount of fluid should diffuse with the ever-
increasing amounts of entrained gas and eventually the composite density tend
to the one of the latter. From Eq. (13) we get $\hat{p}\rightarrow 0$, as does
$p$ then, which follows from the zero velocity. For the entrainment rate
coefficient, taking the approximation of Eq. (17) for $x\gg D_{0}$, i.e.
$\hat{D}\gg 1$, we get the far field entrainment rate constant

$K_{e,\,far}=2\tan\theta.$ (19)

Considering a half-angle of $8^{\circ}\lesssim\theta\lesssim 10^{\circ}$ as
found in the experiments by Hiroyasu and Arai (1990) for high-speed (over 70
m/s) pressure-atomized liquid jets, we get $0.28\lesssim K_{e,far}\lesssim
0.35$. Remarkably, the value 0.32 used by Ricou and Spalding (1961) for
gaseous jets in their experiments falls into the range above, suggesting a
half-angle of $9^{\circ}$.

By comparing Eqs. (18) and (19), we have

$K_{e,near}=\frac{2\rho_{*}^{1/2}}{\rho_{*}+1}K_{e,\,far}.$ (20)

We can check that the coefficient of $K_{e,far}$ in the right-hand-side of Eq.
(20) is always positive and less than one since $\rho_{*}>0$. Therefore we
prove that $K_{e,near}<K_{e,far}$, implying that the entrainment rate has a
tendency of growing with the distance from the nozzle, independently of the
density ratio $\rho_{*}$. We have thus confirmed that Eq. (17) is a
generalization of the empirical one given by Eq. (1), and provides a
theoretical account for it.

###  5.2 Liquid jet without breakup

In this case, the liquid jet does not disintegrate, i.e. it does not atomize,
and the jet consists of pure liquid. This corresponds to a spread half-angle
zero, $\theta=0$. Taking $\theta\rightarrow 0$ in Eq. (11) and thus
$\hat{D}\rightarrow 1$, we immediately see that $\hat{v}\rightarrow 1$, so
that $v\rightarrow v_{0}$, which has the physical interpretation that the
liquid keeps propagating at the same initial speed. As for the density, we see
from Eq. (12) that $\rho\rightarrow\rho_{0}$ since $\hat{\rho}\rightarrow 1$,
which again agrees with the physical meaning as we should only have liquid.
$\hat{p}\rightarrow 1$ so that the dynamic pressure is unaltered along the
jet’s flight. Finally it is immediate to see from Eq. (17) that if $\theta=0$
then $K_{e}=0$ as expected.

###  5.3 The ideal liquid jet in a very thin atmosphere

In its dimensionless form, this condition means $\rho_{*}\approx 0$. With this
assumption, $\rho_{e}\approx 0$ or $\rho_{e}\ll\rho_{0}$, which physically
means a liquid jet inside a very thin ambient gas. Taking the limit in Eq.
(11) we see that $\hat{v}\rightarrow 1$, which corresponds to $v\rightarrow
v_{0}$; i.e. the velocity remains constant at the same initial value, which is
physically acceptable in the sense that liquid droplets would not loose any
momentum to the surrounding thin ambient gas since there would be little drag
forces present. As for the density, taking the same limit in Eq. (12), we see
that $\tilde{\rho}\rightarrow 1/2\hat{D}^{2}$ implies
$\hat{\rho}=\hat{D}^{-2}$, and $\hat{p}=\hat{D}^{-6}$. The entrainment rate in
this case is also zero as can be verified from (17) by taking the limit
$\rho_{*}\rightarrow 0$; this is physically consistent as well.

From the above, in dimensional form, the state equations for an ideal momentum
atomizing liquid jet in a thin atmosphere are then

$\displaystyle
v=v_{0},\quad\rho=\rho_{0}\left(\frac{D_{0}}{D_{0}+2x\tan\theta}\right)^{2},$
(21) $\displaystyle
p=p_{0}\left(\frac{D_{0}}{D_{0}+2x\tan\theta}\right)^{6},\quad\textrm{and
}K_{e}(x)=0;$ (22)

where we can see that as $x\rightarrow 0$, $\rho\rightarrow\rho_{0}$ and
$p\rightarrow p_{0}$, accordingly. As $x\rightarrow\infty$, $\rho\rightarrow
0$ and $p\rightarrow 0$, the approximately null density and dynamic pressure
of the thin gas, respectively.

###  5.4 The submerged jet

This is a single-phase jet, e.g. liquid-liquid or gas-gas jet;
correspondingly, we should have $\rho_{e}=\rho_{0}$. This case corresponds to
$\rho_{*}=1$ in dimensionless form. From Eq. (12) we have that
$\rho=\rho_{0}(\equiv\rho_{e})$, as expected so the density is just that of
the common fluid for all distances. Landau and Lifshitz called this a
“submerged jet” while Batchelor called it a “point source of momentum jet”.

From Eqs. (11), (13) and (17) we get that for this case
$\hat{v}=\hat{D}^{-1}$, $\hat{p}=\hat{D}^{-2}$ and $K_{e}(x)=2\tan\theta$; so
that, in dimensional form, the state equations for the ideal momentum
submerged jet are

$\displaystyle
v=v_{0}\left(\frac{D_{0}}{D_{0}+2x\tan\theta}\right),\quad\rho=\rho_{0},$ (23)
$\displaystyle
p=p_{0}\left(\frac{D_{0}}{D_{0}+2x\tan\theta}\right)^{2},\quad\textrm{and
}K_{e}=2\tan\theta;$ (24)

which for $x\rightarrow 0$, i.e. at the nozzle’s exit, still gives $v=v_{0}$
and $p=p_{0}$, while for $x\rightarrow\infty$, $v=p=0$, as is physically
consistent. Notice that in this case density remains constant and velocity
decays with distance whilst for the case of a thin ambient gas we had the
converse, i.e. the velocity remains constant and the density decays. Also, the
entrainment rate coefficient is constant and $K_{e}=K_{e,far}$ for all
distances, which is the same as the constant coefficient defined by Ricou and
Spalding (1961). This can be interpreted from Eq. (20) where we can see that,
while on the near-field the effect of the density ratio is considerable, in
the far field it is not present. This is because the mass flux through a
cross-section close to the nozzle is highly influenced by the density of the
fluid exiting it, while this effect gradually decreases becoming negligible in
the far-field, as the mass flux of the original fluid is small compared to the
entrained mass flux.

Finally, if we substitute $v_{0}$ from Eq. (4) into Eq. (23) we obtain

$v=2\pi^{-1/2}\left(\frac{\dot{\Pi}_{0}}{\rho}\right)^{1/2}(D_{0}+2x\tan\theta)^{-1}$
(25)

where $\dot{\Pi_{0}}=d\Pi_{0}/dt$ is the momentum flux at the nozzle exit
given in Eq. (5). For $x\gg D_{0}$ we can drop the $D_{0}$ in the last factor
of Eq. (25) obtaining the following proportionality relationship:

$v\propto x^{-1}\sqrt{\dot{\Pi_{0}}/\rho},$ (26)

which is the classical result for the mean velocity decay of a turbulent jet
due to Landau and Lifshitz (1987).

##  6 Comparison with experiments and discussion

We will show the agreement of the ideal atomizing liquid jet model with two
sets of experimental results for the entrainment rate coefficient of two-phase
liquid-gas jets and two more sets for the axial velocity of a single-phase air
jet.

###  6.1 Atomizing diesel jet in air

Cossali et al. (1996); Cossali (2001) reported experimental results for the
air entrainment coefficient of two-phase pressure-atomized diesel jets in air
inside pressure and temperature-controlled chambers, approximating the
conditions inside a fuel-injection combustion engine. For all experiments they
used single-hole small opening nozzles with $D_{0}=0.25$ mm, using standard
diesel fuel of density $\rho_{0}=820\,\mathrm{kg/m^{3}}$. Of these, only one
experiment was performed with a chamber at an ambient pressure of 101.325 kPa
(1 atm) and an chamber temperature equal to an ambient temperature of 25∘C,
denoted as experiment “C”. At this temperature, $\rho_{e}\approx
1.169\,\mathrm{kg/m^{3}}$. The injection pressure is reported therein to have
been 21.4 MPa (214 bar). Using Eq. (4) we can estimate $v_{0}\approx 228.5$
m/s and consequently $\dot{m}_{0}\approx 9.1$ g/s. The results are shown in
Fig. 3. The solid line is Eq. (17), using a parameter $\theta=2.99^{\circ}$
optimized by least squares fitting of the curve. In this case, the trend
agrees with the experimental data for the far field, viz. $\hat{x}\gtrsim 30$
with a correlation coefficient of $R=0.95299$. It must be noted that the best
experimental curve fit available for the same data, as given by the original
authors themselves (Cossali et al., 1996), is not much closer (this curve is
not shown in their article), having a correlation coefficient of $R=0.97705$.
The dotted curve in Fig. 3 labeled “Exp. Fit 8(C)” corresponds to this
experimental power-law fit given in Cossali et al. (1996) as

$K_{e}(x)=B\left(\frac{x+x_{0}}{D_{0}}\right)^{0.5}\left(\frac{\rho_{e}}{\rho_{0}}\right)^{0.31}\left(\frac{T_{g}}{T_{0}}\right),$
(27)

where $B=0.044$ is some scaling parameter, $x_{0}$ is the position of the
nozzle’s exit (defined when $x=0$) and $T_{g}$ and $T_{0}$ are the
temperatures of the gas and liquid, respectively. The exponents of the power
law were calculated using an optimized fit of the parameters in a RMS sense,
adjusting not only for this but for the data of several experiments
altogether. This experimental power-law seems to approximate better in other
cases presented by the same authors.

Finally, we must note that the obtained optimized angle $\theta\approx
3^{\circ}$ is smaller than the angle expected for the whole jet, estimated
between about $8^{\circ}\lesssim\theta\lesssim 10^{\circ}$ as in the
experiments by Hiroyasu and Arai (1990). We must, however, note that there are
great variations in the angle depending on the specific jet characteristics,
most notably the densities ratio $\rho_{*}$ (Lefebvre, 1989; Hiroyasu and
Arai, 1990; Desantes et al., 2005; Xie et al., 2015; Emberson et al., 2016).
The constant-value approximations for the entrainment rate coefficient in the
near and far fields are also shown in Fig. 3 using the same angle in Eqs. (18)
and (19). Here we see that $K_{e,\,far}$ agrees relatively well with the far
field data, while $K_{e,\,near}$ more notably underestimates the experimental
value in this case.

Figure 3: Air entrainment coefficient of a steady pressure-atomized diesel
fuel jet entering stagnant ambient air. Experimental results by Cossali et al.
(1996), far field only ($\hat{x}>30$). Pressure drop was 21.4 MPa and ambient
air temperature was 298${}^{\circ}K$. The optimized jet’s half-angle parameter
is $\theta=2.99\,^{\circ}$ with a correlation coefficient of $R\approxeq
0.953$.

###  6.2 Atomizing water jet in air

Ruff et al. (1989) reported experimental results for a pressure-atomized water
jet entering stagnant ambient air, injected downward using a nozzle with
$D_{0}=9.5$ mm, a fuel density of $\rho_{0}$ = 998 kg/m3, ambient temperature
$T_{g}=298\pm 2^{\circ}$K, reporting for water a Reynolds number (Re = $\rho
vD/\mu$, where $\mu$ is the kinematic viscosity) Re = 5.34$\times 10^{5}$ and
Weber number (We = $\rho v^{2}D/\sigma$, where $\sigma$ is the surface tension
coefficient) of We = 411.5$\times 10^{3}$, while for air We = 492.8. The
results for the air entrainment coefficient using a pressure drop of 2.52 MPa
are shown in Fig. 4 along with our model fit, Eqs. (17), (18) and (19), using
a jet’s half angle parameter of $\theta=1.89^{\circ}$ and obtaining a
correlation coefficient $R=0.98768$. We estimate $v_{0}=71.05$ m/s from Eq.
(4). The power-law experimental fit by Cossali et al. (1996) is also shown in
the same figure getting a correlation coefficient $R=0.98189$, slightly lower
than with our model. It must be noted, however, that in their paper Cossali et
al. gave $B=0.044$, but using this value gives a very bad fit for the
presently considered experiment by Ruff et al. (1989). Therefore, we
recalculated the value of this parameter by least-squares optimization,
obtaining $B=0.0204$. Notably our model seems to adjust very well in the far
field, while still not being so far in the intermediate field. The
experimental power-law follows an opposite trend in this case, giving a better
fit in the near field but differing in the far field. $K_{e,near}$
underestimates the experimental values in the near field but is still a
reasonable approximation. $K_{e,far}$ does seem to overestimate the far field
trend considerably. We must remark that the optimized value of the half-angle,
$\theta\approx 2^{\circ}$, is again smaller than the expected value for the
whole jet geometry (Hosoya and Obokata, 1993). However, which angle to use as
a parameter is not clear, as there is no single universally agreed definition
of $\theta$ (Lefebvre, 1989; Hiroyasu and Arai, 1990; Desantes et al., 2005;
Xie et al., 2015; Emberson et al., 2016). We will give our own original
approach to this issue in Sect. 6.3 below.

Figure 4: Air entrainment coefficient of a steady pressure-atomized water jet
entering stagnant ambient air. Experimental data from Ruff et al. (1989),
intermediate and far fields only ($\hat{x}>10$). $D_{0}=9.5$ mm and pressure
drop of 2.52 MPa. The parameter $\theta=1.89$ is the optimized value,
obtaining a correlation coefficient of $R\approxeq 0.988$.

###  6.3 Single-phase air jet

Although our main focus is for atomizing liquid jets, the availability of the
experimental data for a single-phase air jet made it possible to develop our
own original methods for the analysis, as described below, in particular for
determining the jet’s spread angle, virtual origin and an “equivalent nozzle
diameter”. These original methods can be analogously applied to the atomizing
liquid jet case. The data for our own experiment at DTU was acquired as
described in Sect. 2. For these experiments, the initial velocity was
$v_{0,exp}=30.59$ m/s, as determined experimentally at the axial position of
$\hat{x}=1.16$, remarkably close to the theoretical velocity given by Eq. (4),
$v_{0,bern}=30.51$ m/s. The latter value is the one used to normalize the
other measured velocities. We estimate $Re_{0}\approx 20,000$. Two data sets
were acquired for the single dependent variable of the axial velocity: (i) a
coarse 2D field as a function of the axial and traversal positions as
coordinates; and (ii) a high-resolution sampling 1D line along the jet’s axis.
Both of these data sets are divided in the same near, intermediate and far
fields. A detailed view of the axial velocity in the different fields is shown
in Fig. 5.

Figure 5: Air jet’s axial velocity color map, detailed view of the three
different fields. (a) Near field; (b) intermediate field; and (c) far field.
The velocity is normalized to the Bernoulli theoretical, Eq. (4), and the
distances are normalized by the nozzle’s diameter, $D_{0}$.

From the Fig. 5(a) we can see there is a small jump or acceleration in the
axial velocity of the fluid after exiting the nozzle, reaching a maximum at
about 3 nozzle diameters and decaying thereafter. This initial acceleration
has been reported by several authors in experiments (Rabadi et al., 2007;
Birouk and Lekic, 2009) and, in this case, it is most likely due to the
expansion of the gas as it comes out of the pressurized conditions inside the
nozzle. Our current model does not intend to address this initial acceleration
phenomenon, by which we will hereafter only consider the intermediate and far
fields. In order to do this, we need to calculate an equivalent nozzle
diameter, $d_{eq}$, for a further downstream axial position. This was done by
inspecting the jet’s geometry and determining an appropriate angle. However,
this is not a straightforward task, as there are many possible angle
definitions, depending on the focus of different authors (Lefebvre, 1989; Xie
et al., 2015; Emberson et al., 2016). A common procedure based on the velocity
distribution, as done by Desantes et al. (2011) is to define the “half-
velocity angle”, i.e. the angle for which the velocity has decayed to half the
centerline maximal value. This, however, does not give a full account of the
whole limits of the jet, so we here propose to select an angle that contains
99% of the total momentum of the jet. In order to do this, we have defined and
calculated “Cumulative Momentum Radii” (CMRs), i.e. for each axial position,
we calculate the radius that contains a certain proportion of the jet’s total
momentum in that cross section. These CMRs are shown in Fig. 6. Note that the
near field has not been taken into account to calculate the CMRs, the reason
being that CMRs were calculated using Gaussian curve fits for each velocity
profile at every axial position, and such curves are not a good approximation
to the velocity profile in the near field, as can be seen in Fig. 7. It is
worth noting that the exit velocity at the nozzle is very close to that
predicted by the Bernoulli theorem ($v_{0,\,bern}=30.51$), as can also be seen
in the Fig. 7, as the velocity used for normalization is precisely $v_{0}$
which is calculated from the Bernoulli theorem in Eq. (4). From the same
figure, we observe that the Gaussian profiles are a reasonable approximation
only after about 5 or 6 nozzle diameters; before that, the central portion is
evidently flat and the profile is more consistent with a “top-hat” shape. The
intermediate field dataset starts at $\hat{x}\approx 7$ (as shown in Fig. 5-b)
so the Gaussian profiles are a reasonable approximation thereafter for these
datasets, although they tend to underestimate the velocity at the tails. This,
however, can be neglected for the purpose of calculating the CMRs as they
contain a very small proportion of the total momentum. In the intermediate
field, the Gaussian fits are a good approximation overall. The lines shown in
Fig. 6 are least-squares linear fits (shown in logarithmic scale). The CMR99
(99% CMR) linear fit was calculated by fitting a line with 0.5 as its
intercept (one radius), forcing it to cross at the nozzle’s exit edge. This
defined a virtual origin of the jet inside the nozzle at $\hat{x_{0}}=-2.62$
(and $\hat{y}=0$); this virtual origin was used as an anchor for the fitting
of the lines for subsequent CMRs, so as to produce a consistent set of lines
with the same virtual origin defining the concentric cones geometry of the
jet.

Figure 6: Cumulative Momentum Radii (CMRs) of the air jet and corresponding
angle fits (in degrees). Figure 7: Gaussian curve fits for the axial velocity
profiles of the air jet. The velocity is normalized to the Bernoulli initial
velocity, $v_{0}$. Near field, including actual experimental data.

The initial acceleration of the velocity in the near field, as shown in Fig.
5-a, is a physical phenomenon that was not intended to be dealt with in the
present ideal momentum jet model, by which we cannot apply the theory to this
initial portion of the experimental data. One solution to this is to skip the
near field range and define a “virtual nozzle exit” at some distance
downstream from the real nozzle exit, as shown in Fig. 8. We thus define the
position of the virtual nozzle exit as the first value of the intermediate
field dataset (leftmost value of the horizontal axis in Fig. 5-b), and denote
it by $\hat{x}_{VN}=7.06$. With this virtual nozzle exit position, we can
define the diameter of the virtual nozzle exit using the CMR99 line shown in
Fig. 6, obtaining $\hat{d}_{VN}=1.85$. In the following, we take
$\hat{d}_{VN}$, the corresponding experimental axial velocity, denoted by
$\hat{v}_{VN}=0.88$ (or $v_{VN}=26.88$ m/s), and initial position for the
intermediate field to be the boundary values “at the virtual nozzle exit” for
input in our dynamical models, so as to avoid the initial accelerating region
affecting the calculations. We will thus re-normalize our calculations using
these values, normalizing distances by $\hat{d}_{VN}$, velocities by
$\hat{v}_{VN}$ and shifting axial positions by $\hat{x}_{VN}$.

Figure 8: Diagram of the virtual nozzle exit definition and associated
dynamical parameters.

With the above virtual nozzle setting, the evolution of the jet’s centerline
axial velocity is shown in Fig. 9. The vertical scale is normalized to
$v_{VN}$. All the parameter optimizations were done using only the
intermediate and far fields data; in this case the high resolution centerline
dataset was used. The ideal momentum submerged jet curve given by Eq. (23) is
shown as the dashed curve with the label “ideal-intfar”, which shows excellent
agreement with the experimental data. In the same figure, the solid curve
shows also Eq. (23), but having fixed the initial velocity to $v_{VN}$ and
optimizing the angle only. The other curve shown in a dash-dot line is the
result of a model not described in this paper but also developed by the
present authors; this “lossy energy jet model” and other additional related
models will be described in detail in a sequel article. Suffice it to say now
that an extension of the present model using partial energy conservation by
introducing an energy loss factor, which effectively acts as a simple
turbulence model, and producing a numerical scheme (shown as the dash-dot line
in Fig. 9), gives the best result in terms of accuracy, but the small increase
in accuracy is at the expense of losing analytical solutions.

Figure 9: Axial velocity of the jet vs. axial distance. Results for a pressure-atomized single-phase turbulent air jet. Table 2: Least-squares optimized parameters for the model fittings (except “*” which has been fixed before-hand) and the closes corresponding cumulative momentum radius (CMR). Model | $\theta\,(^{\circ})$ | CMR | $\hat{v_{0}}$ | $\hat{H_{e}}$
---|---|---|---|---
Ideal Momentum | 7.83 | 96 | 1.11 | -
Id.Mo. Bernoulli | 6.75 | 93 | 1.0* | -
Lossy Energy | 4.31 | 80 | 1.04 | 7,053,879

The results of all parameter optimizations by least-squares are shown in Table
2, where $\hat{H}_{e}$ is the dimensionless “energy half-loss” parameter used
in the lossy energy jet model. Overall, the models give similar results for
the initial velocity, where it was left as a free parameter to be optimized,
and the angles are between 4 and 8∘. Focusing on the presently discussed model
of the ideal momentum jet, when $v_{0}$ was left as a free parameter, the
optimization result was still a velocity close to that of the theoretical
Bernoulli velocity (11% greater). Fixing $v_{0}$ or setting it as a free
parameter gave a similar optimized jet half-angle of around 7∘, which is
somewhat lower than the experimental result by Hosoya and Obokata (1993), but
still in the range of other experimental results (Lefebvre, 1989; Xie et al.,
2015; Emberson et al., 2016). In order to determine the best angle for the
model optimization, we investigated the relationship between the cumulative
momentum value, denoted by $\Pi_{c}$, and the jet’s spread half-angle,
$\theta$. Let $\alpha=1-\Pi_{c}$ so that for the 75% CMR, $\alpha=0.25$, for
the CMR99, $\alpha=0.01$, etc. The relationship between $\theta$ and $\alpha$
for the turbulent air jet is shown in Fig. 10, where an experimental curve

$\theta(\alpha)=a(1-\alpha^{b})+c(1-\alpha^{d})$ (28)

optimized by a least squares procedure has produced $a\approx-376474$,
$b\approx-4.8\times 10^{-4}$, $c\approx 112997$, $d\approx-1.6\times 10^{-3}$.
Based on Eq. (28) we can select CMRs closest to a determined half-angle as
shown in Table 2. Thus, as a rough approximation we can recommend the CMR95,
i.e. the cone angle containing approximately 95% of the jet’s total momentum,
as an appropriate angle to be used as a parameter in the current ideal
momentum jet model. For this experiment, the CMR95 corresponds to
$\theta\approx 7.3^{\circ}$.

Figure 10: Half-angle of the jet vs. $\alpha$ tail value of the CMR. Results
for a pressure-atomized single-phase turbulent air jet.

##  7 Conclusions

We have presented a 1D mathematical model applicable to the dynamics of a wide
class of turbulent round jets. The model’s main assumptions include the so-
called Locally Homogeneous Flow (LHF) for a two-phase flow. The model is based
on conservation laws of the momentum and mass, and describes the dynamical
quantities, viz. density, velocity and dynamic pressure, along the jet’s axis.
The main advantages of the model over others in the literature are that the
solutions are analytical, it contains a single free parameter, viz. the jet’s
angle, and the fact that this angle can be approximated from experimental
measurements with an also herein proposed method (Sect. 6.3). In particular,
we theoretically derive a new explicit formula on a sound physical ground,
Eqs. (15) and (17) for the mass entrainment rate. Special cases of the model
give constant approximations of the mass entrainment rate coefficients for the
near and far fields. The comparison with experimental data for atomizing jets
from the literature shows reasonable agreement in the intermediate and far
fields but more data is needed for testing the limits of the theory. Other
special cases considered include a liquid jet in a very thin atmosphere and a
submerged jet. In the latter case, the derived state equations agree with the
classic result for the velocity decay given by Landau and Lifshitz (1987). For
the same special case, we carried out our own experiments with turbulent air
jets, showing excellent agreement with the centerline velocity decay in the
intermediate and far fields. A reasonable approach to finding the jet’s angle
has also been given, introducing the concept of Cumulative Momentum Radii
(CMRs).

The present theory can be extended to include partial conservation of energy
and mass, producing semi-analytical and numerical solutions, improving the
prediction power, as will be described in a subsequent paper.

## Acknowledgements

FFM thanks the support of the Ministry of Education, Sports, Science and
Technology of Japan, the Bank of Mexico’s FIDERH program, the Mexican National
Council of Science and Technology (CONACYT), Kyushu University’s Graduate
School of Mathematics and KUMIAY International Co. This research has also been
partially supported by the Center of Innovation Program from the Japan Science
and Technology Agency, JST. YF was supported in part by a Grant-in-Aid for
Scientific Research from the Japan Society of Promotion of Science (Grant No.
16K05476). AH and CMV wish to gratefully acknowledge the support of DTU MEK
for sponsoring the complete PhD studies program of AH.

## References

  * Abraham [1996] John Abraham.  Entrainment characteristics of transient gas jets.  _Numerical Heat Transfer, Part A: Applications_ , 30(4):347–364, 1996.
  * Batchelor [1967] G. K. Batchelor.  _An introduction to fluid dynamics_.  Cambridge University Press, 1967.  ISBN 978-0-521-66396-0.
  * Birouk and Lekic [2009] Madjid Birouk and Nebojsa Lekic.  Liquid Jet Breakup in Quiescent Atomosphere: A Review.  _Atomization and Sprays_ , 19(6):501–528, 2009\.
  * Cossali [2001] G. E. Cossali.  An integral model for gas entrainment into full cone sprays.  _J. Fluid Mech._ , 439:353–366, 2001.
  * Cossali et al. [1996] G. E. Cossali, A. Gerla, A. Coghe, and G. Brunello.  Effect of Gas Density and Temperature on Air Entrainment in Transient Diesel Spray.  _SAE Tech. Paper_ , 960862, February 1996.  ISSN 0148-7191.
  * Desantes et al. [2005] J. M. Desantes, R. Payri, F. J. Salvador, and A. Gil.  Development and validation of a theoretical model for diesel spray penetration.  _Fuel_ , 85(2006):910–917, 2005.
  * Desantes et al. [2011] J. M. Desantes, F. J. Salvador, J. J. López, and J. De La Morena.  Study of mass and momentum transfer in diesel sprays based on X-ray mass distribution measurements and on a theoretical derivation.  _Exp. Fluids_ , 2011(50):233–246, 2011.
  * Emberson et al. [2016] D. R. Emberson, B. Ihracska, S. Imran, and A. Diez.  Optical characterization of Diesel and water emulsion fuel injection sprays using shadowgraphy.  _Fuel_ , 172:253–262, 2016.  ISSN 0016-2361.
  * Ewing et al. [2007] D. Ewing, B. Frohnapfel, W. K. George, J. M. Pedersen, and J. Westerweel.  Two-point similarity in the round jet.  _J. Fluid Mech._ , 557(309330), 2007.
  * Faeth [1983] G. M. Faeth.  Evaporation and Combustion of Sprays.  _Prog. Energy Combust. Sci._ , 9:1–76, 1983.
  * Faeth [1987] G. M. Faeth.  Mixing, Transport and Combustion in Sprays.  _Prog. Energy Combust. Sci._ , 13:293–345, 1987.
  * Faeth et al. [1995] G. M. Faeth, L.-P. Hsiang, and P.-K. Wu.  Structure and Breakup Properties of Sprays.  _Int. J. Multiphase Flows_ , 21:99–127, 1995.
  * Fuchimoto et al. [2009] Tetsuya Fuchimoto, Shinichiro Yanase, Jiro Mizushima, and Jiro Senda.  Dynamics of vortex rings in the spray from a swirl injector.  _Fluid Dyn. Res._ , 41(4):045503, May 2009.
  * Gamard et al. [2004] S. Gamard, D. Jung, and W.K. George.  Downstream evolution of the most energetic modes in a turbulent axisymmetric jet at high Reynolds number. Part 2: The far-field region.  _J. Fluid Mech._ , 514:205–230, 2004.
  * Ghosh and Hunt [1994] S. Ghosh and J. C. R. Hunt.  Induced air velocity within droplet driven sprays.  _Proc. R. Soc. London, Ser. A_ , 444(1920):105–127, January 1994.
  * Gorokhovski and Herrmann [2008] Mikhael Gorokhovski and Marcus Herrmann.  Modeling Primary Atomization.  _Annu. Rev. Fluid Mech._ , 40(1):343–366, 2008\.  doi: 10.1146/annurev.fluid.40.111406.102200.
  * Hill [1972] B. J. Hill.  Measurement of local entrainment rate in the initial region of axisymmetric turbulent air jets.  _J. Fluid Mech._ , 51(04):773–779, 1972.
  * Hiroyasu and Arai [1990] Hiro Hiroyasu and Masataka Arai.  Structures of Fuel Sprays in Diesel Engines.  _SAE Tech. Paper_ , 900475, 1990.
  * Hosoya and Obokata [1993] Hajime Hosoya and Tomio Obokata.  Efect of Nozzle Configuration on Characteristics of Steady-State Diesel Spray.  _SAE Tech. Paper_ , 930593(3), 1993.
  * Hussein et al. [1994] H. J. Hussein, S. P. Capp, and W. K. George.  Velocity measurements in a high-Reynolds-number, momentum-conserving, axisymmetric, turbulent jet.  _J. Fluid Mech._ , 258(1):31–75, April 1994.
  * Jiang et al. [2010] X. Jiang, G. A. Siamas, K. Jagus, and T. G. Karayiannis.  Physical modelling and advanced simulations of gas-liquid two-phase jet flows in atomization and sprays.  _Prog. Energy Combust. Sci._ , 36:131–167, 2010.
  * Jung et al. [2004] D. Jung, S. Gamard, and W.K. George.  Downstream evolution of the most energetic modes in a turbulent axisymmetric jet at high Reynolds number. Part 1: The near-field region.  _J. Fluid Mech._ , 514:173–204, 2004.
  * Landau and Lifshitz [1987] L. D. Landau and E. M. Lifshitz.  _Fluid Mechanics_ , volume 6 of _Course of Theoretical Physics_.  Butterworth Heinemann, 2nd edition, 1987.
  * Lefebvre [1989] Arthur H. Lefebvre.  Properties of Sprays.  _Part. Part. Syst. Charact._ , 6(1989):176–186, 1989.
  * Lin and Reitz [1998] S. P. Lin and R. D. Reitz.  Drop and spray formation from a liquid jet.  _Annu. Rev. Fluid Mech._ , 30:85–105, 1998.
  * Liu [2000] Huimin Liu.  _Science and Engineering of Droplets_.  Noyes Publications / William Andrew Publishing, 2000.
  * McCarthy and Molloy [1974] M. J. McCarthy and N. A. Molloy.  Review of Stability of Liquid Jets and the Influence of Nozzle Design.  _Chem Eng. J._ , 7:1–20, 1974.
  * Pastor et al. [2008] José V. Pastor, J. Javier López, José M. García, and José M. Pastor.  A 1D model for the description of mixing-controlled inert diesel sprays.  _Fuel_ , 87(13):2871–2885, 2008.
  * Post et al. [2000] Scott Post, Venkatraman Iyer, and John Abraham.  A Study of Near-Field Entrainment in Gas Jets and Sprays Under Diesel Conditions.  _J. Fluids Eng._ , 122:385–395, jun 2000.
  * Rabadi et al. [2007] Said Al Rabadi, Lutz Friedel, and Robert Surma.  Prediction of Droplet Velocities and Rain out in Horizontal Isothermal Free Jet Flows of Air and Viscous Liquid in Stagnant Ambient Air.  _Chem. Eng. Tech._ , 30(11):1546–1563, 2007.
  * Rayleigh [1878] J. S. W. Rayleigh.  On the instability of jets.  _Proc. Lond. Math. Soc._ , 10:4, 1878.
  * Ricou and Spalding [1961] F. P. Ricou and D. B. Spalding.  Measurements of entrainment by axysymmetrical turbulent jets.  _J. Fluid Mech._ , 11(01):21–32, August 1961\.
  * Ruff et al. [1989] G. A. Ruff, A. D. Sagar, and G. M. Faeth.  Structure and Mixing Properties of Pressure-Atomized Sprays.  _AIAA J._ , 27(7):901–908, July 1989.
  * Sazhin [2006] S. Sazhin.  Multiple scales in spray modelling.  _J. Phys. Conf._ , 55:191–202, 2006.
  * Sazhin et al. [2001] S. S. Sazhin, G. Feng, and M. R. Heikal.  A model for fuel spray penetration.  _Fuel_ , 80:2171–2180, 2001.
  * Song and Abraham [2003] Lijun Song and John Abraham.  Entrainment Characteristics of Transient Turbulent Round, Radial and Wall-Impinging Jets: Theoretical Deductions.  _J. Fluids Eng._ , 125(4):605–612, August 2003\.  ISSN 0098-2202.  doi: 10.1115/1.1593707.
  * Stevenin et al. [2016] C. Stevenin, A. Vallet, S. Tomas, M. Amielh, and F. Anselmet.  Eulerian atomization modeling of a pressure-atomized spray for sprinkler irrigation.  _Int. J. Heat and Fluid Flow_ , 57:142–149, 2016.  ISSN 0142-727X.
  * Taylor [1962] G. I. Taylor.  _Generation of ripples by wind blowing over viscous fluids_ , volume 36.  Cambridge University Press, 1962.
  * Tomotika [1935] S. Tomotika.  On the instability of a cylindrical thread of a viscous liquid surrounded by another viscous fluid.  _Proc. R. Soc. London, Ser. A_ , 150:322–337, 1935.
  * Vallet et al. [2001] A. Vallet, A. A. Burluka, and R. Borghi.  Development of an Eulerian model for the atomization of a liquid jet.  _Atomization and Sprays_ , 11(6), 2001.  ISSN 1044-5110.
  * Wakuri et al. [1960] Yutaro Wakuri, Masaru Fujii, Tatsuo Amitani, and Reijiro Tsuneya.  Studies on the penetration of fuel spray in a diesel engine.  _Bull. Japan Soc. Mech. Eng._ , 3(900475):123–130, 1960.
  * Weber [1931] C. Z. Weber.  Zum Zerfall eines Flussigkeitsstrahles.  _Math. Mech._ , 11:136–154, 1931.
  * Xie et al. [2015] Hongzhan Xie, Lanbo Song, Yizhi Xie, Dong Pi, Chunyu Shao, and Qizhao Lin.  An Experimental Study on the Macroscopic Spray Characteristics of Biodiesel and Diesel in a Constant Volume Chamber.  _Energies_ , 8(6):5952–5972, 2015.


# Design of an Externally Driven Current Cloak

A. S. García-Gordillo  Group of Complex Systems and Statistical Physics,
Physics Faculty, University of Havana, 10400 Havana, Cuba  E. Altshuler
<EMAIL_ADDRESS>Group of Complex Systems and Statistical Physics,
Physics Faculty, University of Havana, 10400 Havana, Cuba

###### Abstract

An inhomogeneity into a conductive matrix deforms the flow pattern of an
applied electric current. A usual current cloak can be defined as a permanent
modification of the matrix properties around the inhomogeneity guaranteeing
that the current flow pattern is similar before and after passing by the
modified zone, so it implies the “electrical invisibility” of the
inhomogeneous region. Here we introduce the concept of a current cloak that
can be tuned –switched on and off, for example– by means on an external field.
We demonstrate analytically and using Finite Elements Simulations that a
current cloak can be constructed and manipulated by an external magnetic field
for a concrete system consisting in a magneto-resistive matrix with a
stainless steel inclusion.

###### pacs:

72.80.Tm, 81.05.Ni, 72.15.Gd

##  I Introduction

Lack of homogeneity in materials is generally assumed as deleterious:
uncontrolled pores, second phase inclusions, grain boundaries and many other
defects tend to weaken mechanical, electrical and magnetic properties
Frankel1957 ; Altshuler1999 ; Batista1999 ; Batista2003_resistive ;
Batista2003 ; Mo2013 . In fact, one of the main aims of Materials Science is
to improve fabrication methods in order to eliminate them.

One alternative approach to the problem is not eliminating, but concealing
these defects: this is at the core of the emerging field of cloaks. They can
be defined as materials modified in such a way that an inhomogeneous region
inside them is made invisible (i.e., indistinguishable from its surroundings)
to electromagnetic fields or sound waves, for example. The subject was kicked
off by the pioneering theoretical work of Pendry and Leohardt in optical
cloaks Pendry2006 ; Leonhardt2006 : if a region of space has optical
properties different from the rest of the material, it is possible to engineer
a change of properties in the vicinity of the inhomogeneity so the rays of
light detected from a point far from the region follow the same trajectories
they had before entering the region in question. Pendry and Leohardt coined
the term metamaterial to describe a material engineered in such way. Over the
following years, several types of cloaks were reported in the optical
Valentine2009 ; Gabrielli2009 ; Smolyaninov2009 ; Ergin2010 ; Rama2015  and
microwave Schurig2006 ; Liu2009 ; Tretyakov2009 ; Edwards2009 ; Ma2010
frequency ranges.

Later on, cloaks invaded the scenario of stationary (or near-stationary)
fields. Wood and Pendry were the first to propose a cloak to conceal a defect
in the presence of a low-frequency magnetic field Wood2007 . The idea was
experimentally proven a few years later by a combination of two materials with
“opposite” magnetic behaviours: a superconductor and a ferromagnet
Narayana2012 . The boundary was immediately pushed further to include
stationary fields Gomory2012 . Yungui Ma et al. designed and obtained
experimentally a “bifunctional cloak” able to conceal the lines of electrical
current (current cloak) or heat flux lines (heat cloak) around a cavity in a
metal, or an inclusion with conductivity different from its surroundings
Ma2014 . Importantly, in all the cloaks described above the material
properties must be permanently modified around the inhomogeneity in order to
achieve the desired effect.

Here, we propose a new way to achieve the cloaking effect: instead of
permanently modifying the material, an external field is used to “tune” or
“drive” the material properties around the inhomogeneity. We illustrate the
idea with a current cloak able to suppress the deformation of the current
lines around an inclusion of larger conductivity within a magneto-resistive
matrix. Firstly, the feasibility of achieving the cloaking effect is
demonstrated analytically. Then, we show it by means of Finite Elements
Simulations for a realistic set of materials and magnetic fields. We call the
new metamaterial Externally Driven Current Cloak or, for our particular case,
Magnetically Driven Current Cloak. Finally, we briefly examine the interaction
of cloaks between them, and with a boundary.

##  II The magnetically driven current cloak

###  II.1 General description

Fig. 1(a) shows a sketch of a 2D version of the magnetic field driven current
cloak. The inner circle of radius $R_{c}$ is the inclusion (or core), whose
conductivity $\sigma_{c}$ is larger than that of the magneto-resistive matrix
at zero applied field, $\sigma_{m}$. The circle of radius $R_{B}$ is the area
covered by a magnetic field $B$ which is applied in a direction perpendicular
to the plane of the matrix. As a result, the conductivity increases from
$\sigma_{m}$ to $\sigma_{B}$ within a ring of inner and outer radii $R_{c}$
and $R_{B}$, respectively.

Fig. 1 (b)-(d) illustrates the working principle of the cloak: the inner and
outer circumferences have radii $R_{c}$ and $R_{B}$ respectively, which are
not explicitly indicated for clarity; the black lines correspond to an
electric field (or current) which is uniform far from the core. Fig. 1 (b)
shows the matrix and the core at zero magnetic field: the field lines “bend
towards” the core, due to the fact that $\sigma_{c}>\sigma_{m}$. Fig. 1 (c)
corresponds to the system at a magnetic field $B_{cloak}$ applied within the
region of radius $R_{B}$, without showing the effect of the core: due to the
decrease of the conductivity in the magneto-resistive ring, the current lines
“bend away” from the core –an effect that compensates that described in (b).
Fig. 1 (d) illustrates the effect of $B_{cloak}$ on the matrix-core composite,
so the flow lines are the overlapping between the patterns sketched in (b) and
(c): the result is that the current flow becomes homogeneous immediately
outside the region containing the magnetic field –the current cloak effect has
been achieved. Now, we calculate the magnetic field needed to produce the
cloaking, $B_{cloak}$, as a function of the radii and conductivities.

Figure 1:  A two-dimensional current cloak driven by an external magnetic
field. (a) Parameters used to derive the cloaking condition (see text). (b)
Magneto-resistive matrix with inclusion at zero magnetic field. (c) Matrix
alone at cloaking magnetic field, $B_{cloak}$ (d) Matrix with inclusion at
$B_{cloak}$. In (b-d), darker blue indicates larger conductivity, black lines
represent the electric current (or electric field) lines. In (c and d), the
applied magnetic field is represented by “x”s.

###  II.2 An analytical expression for the cloaking condition

The continuity equation for an electric current density $\vec{J}$, represented
by the black lines in Fig. 1 reads $\vec{\nabla}\cdot\vec{J}=0$. If $\vec{E}$,
$V$ and $\sigma$ are the electric field, the potential and the conductivity,
one can write $\vec{J}=\sigma\vec{E}$ (where $\vec{E}=-\vec{\nabla}V$). By
combining the previous expressions, we get a Laplace equation for the
potential $V$ in cylindrical coordinates.

If we solve it for the system illustrated in Fig. 1(a) and impose the cloaking
condition, i.e., $\vec{E}$ is uniform outside the area where the magnetic
field is applied, we get (see Supporting Online Material):

$\sigma_{m}=\frac{R_{B}^{2}(\sigma_{c}+\sigma_{B})+R_{c}^{2}(\sigma_{c}-\sigma_{B})}{R_{B}^{2}(\sigma_{c}+\sigma_{B})+R_{c}^{2}(\sigma_{B}-\sigma_{c})}\sigma_{B}$
(1)

It is worth noting that the cloaking condition reported by Ma et al. for a
core consisting in a void Ma2014  is a particular case of Eq. 1 corresponding
to $\sigma_{c}=0$. Both conditions stand for any value of applied electric
field, i.e., no matter the value the uniform current injected through the
sample takes, the cloaking parameters will remain the same.

If we follow the thread of previous work in the field, the cloaking condition
implies inserting a ring of conductivity $\sigma_{B}$ and inner and outer
radii $R_{c}$ and $R_{B}$, respectively, so all the parameters fulfill Eq. 1.
Our approach is different: since $\sigma_{B}$ is a function of the magnetic
field, we propose applying a magnetic field Bcloak on the ring-shaped region
in such a way that the resulting conductivity matches the value given by Eq.
1. Naturally, $B_{cloak}$ depends on the magneto-resistive behaviour of the
matrix.

###  II.3 A realistic magnetically driven current cloak

Now, we analyze the case of a composite made of a polycrystalline Bi matrix
shaped as a 25 $\times$ 25 $\times$ 1 mm3 sheet, whose magneto-resistive
behavior is illustrated in Fig. 2 (8.67 $\times$ 105 S/m at zero field
Yang1999 ). A cylindrical inclusion of stainless steel of radius 1.6 mm is
located at the center of the matrix. The conductivity of the inclusion –which
is very weakly dependent on the magnetic field– is taken as 9.93 $\times$ 106
S/m Conductivities . The external field is confined to a region of radius 3.2
mm around the inclusion, which can be easily achieved in practice by using a
couple of cylindrical super-magnets of appropriate radii.

Fig. 3 shows FEM simulations of the electric current density lines along the
sample at different magnetic fields perpendicular to the large face of the
matrix (a current of $5.5$ A is uniformly injected along the left vertical
boundary of the matrix with a uniform density $J_{in}$ $=$ $2.2$ $\times$ 105
A/m2 in the case of no inclusion and zero applied magnetic field). For zero
applied field (Fig. 3 (a)) the lines are deformed towards the inclusion,
corresponding to the case sketched in Fig. 1 (b). For a field $B$ = $1.79$ T
(Fig. 3 (b)), the deformation of the current lines towards the highly
conductive core is compensated by the ring of depleted conductivity due to the
effect of the applied magnetic field, so they are homogeneous outside the zone
where the magnetic field is confined. This situation, analogous to the one
sketched in Fig. 1 (d) implies that $B_{cloak}$ $\approx$ $1.79$ T for our
composite. In fact, this is the field value that satisfies Eq. 1, taking for
the specific geometry and conductivities of the cloak, as well as the curve
shown in Fig. 2. For a magnetic field of $3.28$ T (Fig. 3(c)), the
conductivity of the ring around the inclusion has decreased to a point that it
“over-compensates” the original deformation of the current lines.

Figure 2:  Magnetic field dependence of conductivity for polycrystalline Bi
(data taken from Yang1999 ).  Figure 3:  Magnetically driven cloaking for a
realistic sample. FEM simulations of the current density lines for a 1-mm
thickness polycrystalline Bi matrix with a stainless steel inclusion of
$R_{c}=$ 1.6 mm radius (inner circle), and a magnetic field applied within a
region of radius $R_{B}=$ 3.2 mm (outer circle). (a), (b) and (c) correspond
to magnetic fields of $0$ T, $1.79$ T and $3.28$ T, respectively. Here,
$B_{cloak}$ $\approx$ 1.79 T.

In order to quantify the current density distribution near the cloak, and
properly estimate $B_{cloak}$, a more detailed analysis is needed. Fig. 4
shows how the electric current density changes along a vertical line running
from the top to the bottom of the sample, which is tangent to the region where
the magnetic field is applied (see inset of Fig. 4). The two vertical dotted
lines in the main graph correspond to the upper and lower edges of the region
where the magnetic field is applied, respectively. For $B$ = $0$, the current
density above and below that region decreases from its value when there is no
inclusion nor magnetic field applied, since the current lines “bend towards”
the higher conductivity core. The bending is smaller for $B$ = $0.96$ T and
virtually non-existent for $B$ = $1.79$ T. The lines bend in the opposite
direction for fields of $2.47$ T and $3.28$ T. In the region where the
magnetic field is applied, along the same vertical line, the deformation of
the originally uniform current density field is more evident, even so, for $B$
= $1.79$ T there is no tangible sign of deformation. This magnetic field was
chosen according to Eq. 1 for the cloaking condition to be fulfilled. It
should be pointed out that $B_{cloak}$ does not change if the sample is larger
in the x direction (see reference frame of Fig. 1 (a)); as a matter of fact,
the current density profiles of Fig. 4 also remain the same as expected.

Based on the curves shown in Fig. 4, the deviations of the current density
from its uniform value when there is no inclusion nor magnetic field applied
can be quantitatively evaluated by defining a cloaking parameter:

$P_{cloak}=1-\frac{Max(|J(x)-J_{in}|)}{J_{in}}$ (2)

where $J_{in}$ is the uniform current density for no inclusion or magnetic
field applied and $J(x)$ is the current density dependence along the vertical
line displayed in the inset of Fig. 4 for a certain external magnetic field
and $R_{B}$. This allows us to quantify how far we are from achieving the
cloak. If the applied magnetic field is different from $B_{cloak}$, then
$P_{cloak}<1$, suggesting that the device is not in the cloaking state. On the
other hand, if the applied field coincides with $B_{cloak}$, then
$P_{cloak}=1$ and the device will behave as a genuine current cloak.

Figure 4:  Current density values plotted along a vertical line tangent to the
area of application of the magnetic field (see inset). The straightest curve
($B_{cloak}$ = 1.79 T) fulfills the cloaking condition.

By varying the area of the zone where the external magnetic field is applied
(by modifying its radius $R_{B}$) and then varying the value of the applied
magnetic field in each previous configuration, it is possible to build a
“parameter phase diagram” using the definition of $P_{cloak}$ given by Eq. 2.
Fig. 5 illustrates the behavior of $P_{cloak}$ with $R_{B}/R_{c}$ and the
applied magnetic field $B$. These magnitudes were chosen because they can be
externally modified to achieve the cloaking condition for a given inclusion
into a given matrix. The cloaking and non-cloaking states are visible in Fig.
5 in terms of $P_{cloak}$: the red ridge corresponds to the sets of external
parameters that make possible the cloaking. The smaller the radius $R_{B}$ the
strongest the magnetic field we have to apply to achieve the current cloaking
condition. These sets of parameters, for which $P_{cloak}=1$, match those
calculated using Eq. 1 within negligible uncertainty ($\approx$0.1$\%$).
Although the red ridge follows Eq. 1 for a given magneto-resistive matrix, an
analytical expression for the specific cloaking state given by the rest of the
parameters in Fig. 5 is not easy to obtain.

Fig. 5 also illustrates that the cloaking condition is achievable even for
$R_{B}$ values comparable to the size of the matrix, suggesting that its
boundaries do not affect considerably the cloaking effect. This is shown in
the Supporting Online Material.

Figure 5:  Exploring the cloaking parameter space. Dependance of $P_{cloak}$
with $R_{B}$ and the external applied magnetic field $B$: cloaking and non-
cloaking states.

Up to now we have only considered the modification of the current flow pattern
by the presence of one inhomogeneity in our magneto-resistive matrix. If two
identical inclusions deployed side-by-side along the symmetry axis of the
sample along the y-axis, the cloaking affect is again achieved for the same
value of the magnetic field (see Supporting Online Material). This is
presumably the case for a larger number of inclusions, provided the matrix is
long enough along the y-axis.

In summary, we have proposed the idea that current cloaks can be tuned
externally. We demonstrate that the cloaking condition can be achieved in the
case of a circular inclusion embedded in a magneto-resistive plate, by
applying a perpendicular magnetic field of appropriate intensity. We provide
enough quantitative information to attempt an experimental test of the
concept. It is reasonable to believe that other kinds of externally driven
currents cloaks (based, for example, on temperature gradients) can be
conceived and materialized.

##  III Acknowledgements

We acknoweledge Adrián Enrique for suggesting the possibility of using a
magneto-resistive material for the matrix and Y. Nahmad for letting us use
COMSOL Multiphysics for the simulations. This research was partially supported
by the University of Havana’s Institutional Project “Superconductors and
conductors: from characterization to applications”. E. A. drew inspiration
from the late M. Álvarez-Ponte.

##  Appendix A Supporting Online Material

###  A.1 An analytical expression for the cloaking condition

We consider a cylindrical inclusion (infinite along z) of radius $R_{c}$,
whose conductivity $\sigma_{c}$ is larger than that of the magneto-resistive
matrix $\sigma_{m}$ where it is immersed. In a concentric region of radius
$R_{B}$ an external magnetic field is applied. Due to the presence of this
field the conductivity of a hollow cylinder of inner and outer radii $R_{c}$
and $R_{B}$ changes to $\sigma_{B}$, being $\sigma_{B}$ smaller than
$\sigma_{m}$ for any non-zero magnetic field values. The electric field at
infinity is shown in Fig. 6 with a sketch of the system.

To deal with an electric stationary problem like the current cloak we used the
current density continuity equation:

$\textbf{$\vec{\nabla}$}\cdot\textbf{$\vec{J}$}=0$ (3)

considering that $\vec{J}=\sigma\vec{E}$, $\vec{E}=-\vec{\nabla}V$, and the
divergence of the gradient of a function is equal to the laplacian of itself,
we get:

$\begin{split}\textbf{$\vec{\nabla}$}\cdot(\sigma\textbf{$\vec{E}$})=0\\\
-\textbf{$\vec{\nabla}$}\cdot\sigma\textbf{$\vec{\nabla}$}V=0\\\
\sigma\textbf{$\nabla^{2}$}V=0\\\ \textbf{$\nabla^{2}$}V=0\end{split}$ (4)

For the present system it is suitable to use cylindrical coordinates. The
Laplace equation for the potential in this coordinate system reads:

$\frac{1}{r}\frac{\partial}{\partial r}(r\frac{\partial V}{\partial
r})+\frac{1}{r^{2}}\frac{\partial^{2}V}{\partial\varphi^{2}}+\frac{\partial^{2}V}{\partial
z^{2}}=0$ (5)

If we take into account the geometry of the problem in the way that the lines
of $\vec{E}$ (parallels to the lines of $\vec{J}$) come from $-\infty$ in the
x axis -orthogonal direction to the z axis- (see Fig. 6) it is not difficult
to realize that:

$V(r,\varphi,z)=V(r,\varphi)$ (6)

hence:

$\frac{\partial^{2}V}{\partial z^{2}}=0$ (7)

This independence regarding the z coordinate allows us to treat as equivalent
three qualitative different systems: the one which is infinite in the z
direction, the one with finite z dimensions and the $2D$ system. The case
presented in this paper, due to practical implications, was the system with
finite dimensions consisting in a foil with an inclusion resembling a coin and
a magnetic field affected washer around it.

Figure 6:  Sketch of the magnetically driven current cloak. Parameters used to
derive the cloaking condition (see text).

If we propose the solution $V(r,\varphi)=R(r)\cos\varphi$ and substitutes it
in Eq 5 (without the z related term), the expression reduces to an Euler’s
equation for $R(r)$:

$r^{2}\frac{\partial^{2}R}{\partial r^{2}}+r\frac{\partial R}{\partial r}-R=0$
(8)

Proposing the classic solution to this equation as $R=r^{k}$, we get the
coefficients $k=\pm 1$. Combining the two solution proposals we get the
following expressions for the potential in the three regions of interest:

$V(r,\varphi)=\begin{cases}\big{(}C_{1}r+\frac{C_{2}}{r}\big{)}\cos\varphi\textbf{
}\textbf{ }\textbf{ }\textbf{ }\textbf{ }\textbf{ }r>R_{B}\\\
\big{(}C_{3}r+\frac{C_{4}}{r}\big{)}\cos\varphi\textbf{ }\textbf{ }\textbf{
}\textbf{ }\textbf{ }\textbf{ }R_{c}<r<R_{B}\\\
\big{(}C_{5}r+\frac{C_{6}}{r}\big{)}\cos\varphi\textbf{ }\textbf{ }\textbf{
}\textbf{ }\textbf{ }\textbf{ }r<R_{c}\end{cases}$ (9)

where $C_{1}$, $C_{2}$, $C_{3}$, $C_{4}$, $C_{5}$ and $C_{6}$ are constants to
be calculated from the boundary conditions.

The expression for the potential at the infinity considering an uniform
electric field of intensity $E_{0}$ parallel to the x axis reads:

$V(\infty,\varphi,z)=-E_{0}r\cos\varphi$ (10)

Moreover, the potential at $r=0$ cannot reach infinite, so it must be bounded:

$\mid V(0,\varphi)\mid<\infty$ (11)

From the expressions 10 and 11, and the solution 9 for $V(r,\varphi)$ we get:

$C_{1}=-E_{0}$ (12)

and

$C_{6}=0$ (13)

remaining only four constants to solve.

In order to calculate these constants we need to look at the boundary
conditions at $R_{c}$ and $R_{B}$. In each of these boundaries stands the
continuity of the potential and the continuity of the normal component of the
current density ($\vec{J}$):

$\begin{split}V_{B}(R_{c},\varphi)=V_{c}(R_{c},\varphi)\\\
V_{m}(R_{B},\varphi)=V_{B}(R_{B},\varphi)\\\
J_{Bn}(R_{c},\varphi)=J_{cn}(R_{c},\varphi)\\\
J_{mn}(R_{B},\varphi)=J_{Bn}(R_{B},\varphi)\end{split}$ (14)

Substituting the potentials and the normal components of the current density
($J_{n}$) evaluated at $R_{c}$ y $R_{B}$, and taking into account the gradient
of $V$ in cylindrical coordinates for which:

Figure 7:  Proximity of the boundaries. (a,b) Current density lines for a 1-mm
thickness polycrystalline Bi matrix with a stainless steel inclusion of
$R_{c}=$ 1.6 mm radius (inner circle), and a magnetic field applied within a
region of radius $R_{B}=$ 3.2 mm (outer circle). Upper and lower boundary
positions were changed allowing a separation of $1/7$ its radius from the
device. (a) and (b) correspond to magnetic fields of $0$ T and $1.79$ T,
respectively. (c) Current density values along a vertical line tangent to the
area of application of the magnetic field (see inset).

$\vec{E}=-\frac{\partial V}{\partial r}\hat{r}-\frac{1}{r}\frac{\partial
V}{\partial\varphi}\hat{\varphi}$ (15)

we obtained the following system of equations in four variables:

$\begin{split}C_{3}R_{c}+\frac{C_{4}}{R_{c}}=C_{5}R_{c}\\\
-E_{0}R_{B}+\frac{C_{2}}{R_{B}}=C_{3}R_{B}+\frac{C_{4}}{R_{B}}\\\
(-C_{3}+\frac{C_{4}}{R_{c}^{2}})\sigma_{B}=-C_{5}\sigma_{c}\\\
(E_{0}+\frac{C_{2}}{R_{B}^{2}})\sigma_{m}=(-C_{3}+\frac{C_{4}}{R_{B}^{2}})\sigma_{B}\end{split}$
(16)

which can be solved for $C_{2}$, $C_{3}$, $C_{4}$ and $C_{5}$ as functions of
$\sigma_{c}$, $\sigma_{B}$, $\sigma_{m}$, $R_{c}$ and $R_{B}$. If we impose
that $C_{2}(\sigma_{c},\sigma_{B},\sigma_{m},R_{c},R_{B})=0$, it is not
difficult to see from Eq. 9 that the electric field outside the area where the
magnetic field is applied will be equal in magnitude and direction to the
uniform electric field at infinity. This constitutes the base of the cloaking
condition: the electric field distribution where the magnetic field is applied
will be concealed and this region will be invisible in terms of electric
currents. If we let $C_{2}(\sigma_{c},\sigma_{B},\sigma_{m},R_{c},R_{B})=0$ in
the solution of Eq. 16 for $C_{2}$ we obtain the cloaking condition:

$\sigma_{m}=\frac{R_{B}^{2}(\sigma_{c}+\sigma_{B})+R_{c}^{2}(\sigma_{c}-\sigma_{B})}{R_{B}^{2}(\sigma_{c}+\sigma_{B})+R_{c}^{2}(\sigma_{B}-\sigma_{c})}\sigma_{B}$
(17)

This condition for the externally driven current cloak can be reduced to a
particular case already published Ma2014 .

If we considered $\sigma_{c}=0$, Eq. 17 reduces to:

$\sigma_{m}=\frac{R_{B}^{2}-R_{c}^{2}}{R_{B}^{2}+R_{c}^{2}}\sigma_{B}$ (18)

This condition was obtained by Yungui Ma et al. working with a bifunctional
cloak Ma2014  and corresponds to a core consisting in a void ($\sigma_{c}=0$).
Here $\sigma_{B}>\sigma_{m}$ and all parameters have to be fixed to properly
fulfill the cloaking condition. Our cloak allows us to apply an external field
to “tune” the cloaking state instead of permanently modifying the material.

It is worth noting that in our case $\sigma_{c}$ must be larger than
$\sigma_{B}$ in order to obtain the cloak. If $\sigma_{c}$ is smaller than
$\sigma_{B}$, the $\sigma_{m}/\sigma_{B}$ ratio in Eq. 17 (looking at the
right-hand side of the equation after dividing both sides by $\sigma_{B}$)
would be smaller than $1$ , letting the conductivity of the region where the
magnetic field is applied (excepting the inclusion) be larger than that of the
external material. This scenario is not possible using a magneto-resistive
material and a magnetic field as described above, but it is part of the
solution and a perfect current cloak can be built with such specs.

###  A.2 Proximity of the boundaries

Figure 8:  Cohabiting cloaks. (a,b) Current density lines for a 1-mm thickness
polycrystalline Bi matrix with two stainless steel inclusions of $R_{c}=$ 1.6
mm radii (inner circles), and magnetic fields applied within concentric
regions of radii $R_{B}=$ 3.2 mm (outer circles). (a) and (b) correspond to
magnetic fields of $0$ T and $1.79$ T, respectively. (c) Current density
values plotted along a vertical line tangent to the area of application of the
magnetic field (see inset).

Fig. 7 (a,b) shows a Finite Elements simulation of the current density lines
along a 1-mm thickness polycrystalline Bi matrix with a stainless steel
inclusion of $R_{c}=$ 1.6 mm radius (inner circle), and a magnetic field
applied within a region of radius $R_{B}=$ 3.2 mm (outer circle). The
positions of the upper and lower boundaries were changed allowing a separation
of $1/7$ its radius from the device. The presence of the highly conductive
inclusion in the absence of magnetic field deforms current density lines (Fig.
7 (a)); on the other hand, when we applied $B_{cloak}$ = $1.79$ T in a region
of radius 3.2 mm concentric with the inclusion (Fig. 7 (b)) the
inhomogeneities are concealed and we obtain a perfect current cloak. Fig. 7
(c) allows us to see in detail what happens to the current density along a
vertical line running from the top to the bottom of the sample (red line
presented in the inset).

The cloak is obtained at the same value of $B_{cloak}$, making confined
scenarios suitable to apply Eq. 17 without any correction. This let us
conclude that, as long as the electric field is uniform in the absence of
magnetic field or inclusion, the proximity of the boundaries is not going to
affect the cloaking state.

###  A.3 Cloaks coexisting together

Fig. 8 (a,b) shows the current density lines along a 1-mm thickness
polycrystalline Bi matrix with two stainless steel inclusions of $R_{c}=$ 1.6
mm radii (inner circles), and magnetic fields applied within concentric
regions of radii $R_{B}=$ 3.2 mm (outer circles). The current density field is
deformed by the highly conductive inclusions at zero magnetic field in the
expected way (Fig. 8 (a)); however, when a magnetic field of induction
$B_{cloak}$ = $1.79$ T is applied in a region of radius 3.2 mm around both
inclusions (Fig. 8 (b)) we obtain a perfect concealing of inhomogeneities.
Fig. 8 (c) exhibits the changes of the current density along a vertical line
running from the top to the bottom of the sample (red line presented in the
inset), emphasizing the attainment of the cloaks at $B$ = $1.79$ T.

If the magnetic field in one of the cloaks (or the area where it is applied)
is not correct, both cloaks are going to be out of the cloaking state. The
conditions stated in Eq. 17 have to be individually satisfied so as to achieve
a collective result. We examined several configurations changing position,
number, and radii of the cores; and, in all cases, we obtained a concealing of
inhomogeneities when all cloaks were perfectly tuned.

## References

  * (1) Jacob Porter Frankel. Principles of the Properties of Materials, 1st Ed. (McGraw-Hill Book Company, New York, Toronto, London 1957)
  * (2) E. Altshuler, R. Cobas, A. J. Batista-Leyva, C. Noda, L. E. Flores, C. Martínez, M. T. D. Orlando Phys. Rev. B 60, 3673-3679 (1999)
  * (3) A. J. Batista-Leyva, R. Cobas, M. T. D. Orlando, C. Noda and E Altshuler, Magnetic hysteresis of the zero-resistance critical temperature in YBaCuO, BiSrCaCuO and HgBaCaCuO superconducting polycrystals Physica C: Superconductivity 314 (1-2), 73-80 (1999).
  * (4) A. J. Batista-Leyva, M. T. D. Orlando, L. Rivero, R. Cobas and E. Altshuler, The resistive transition of (Hg0.85Re0.15)(Ba1-ySry)2Ca2Cu3O8+d superconducting polycrystals Physica C 383, 365-373 (2003).
  * (5) A. J. Batista-Leyva, R. Cobas, M. T. D. Orlando, E. Altshuler, Hysteresis and relaxation in $TlBa_{2}Ca_{2}Cu_{3}O_{y}$ superconducting polycrystals Supercond. Sci. Technol. 16, 857-864 (2003)
  * (6) Wenlin Mo, Shanping Lu, Dianzhong Li, Yiyi Li, Effects of filler metal composition on inclusions and inclusion defects for ER NiCrFe-7 weldments. J. Mat. Sci. Technol. 29, 458-466 (2013)
  * (7) J. B. Pendry, D. Schurig and D. R. Smith. Controlling electromagnetic fields. Science 312, 1780–1782 (2006).
  * (8) U. Leonhardt. Optical conformal mapping. Science 312, 1777–1780 (2006).
  * (9) J. Valentine, J. Li, T. Zentgraf, G. Bartal and X. Zhang. An optical cloak made of dielectrics. Nat. Mater. 8, 568–571 (2009).
  * (10) L. H. Gabrielli, J. Cárdenas, C. B. Poitras and M. Lipson. Silicon nanostructure cloak operating at optical frequencies. Nat. Photon. 3, 461–463 (2009).
  * (11) I. I. Smolyaninov, V. N. Smolyaninova, A. V. Kildishev and V. M. Shalaev. Anisotropic metamaterials emulated by tapered waveguides: application to optical cloaking. Phys. Rev. Lett. 102, 213901 (2009).
  * (12) T. Ergin, N. Stenger, P. Brenner, J. B. Pendry and M. Wegener. Three-dimensional invisibility cloak at optical wavelengths. Science 328, 337–339 (2010).
  * (13) Sarath Ramadurgam, Tzu-Ging Lin, Chen Yang, Tailoring optical and plasmon resonances in core-shell and core-multishell nanowires for visible range negative refraction and plasmonic light harversting: a review. J. Mat, Sci. Technol. 31, 533-542 (2015).
  * (14) D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr and D. R. Smith. Metamaterial electromagnetic cloak at microwave frequencies. Science 314, 977–980 (2006).
  * (15) R. Liu, C. Ji, J. J. Mock, J. Y. Chin, T. J. Cui and D. R. Smith. Broadband ground-plane cloak. Science 323, 366–369 (2009).
  * (16) S. Tretyakov, P. Alitalo, O. Luukkonen and C. Simovski. Broadband electromagnetic cloaking of long cylindrical objects. Phys. Rev. Lett. 103, 103905 (2009).
  * (17) B. Edwards, A. Alu, M. G. Silveirinha y N. Engheta. Experimental verification of plasmonic cloaking at microwave frequencies with metamaterials. Phys. Rev. Lett. 103, 153901 (2009).
  * (18) H. F. Ma and T. J. Cui. Three-dimensional broadband ground-plane cloak made of metamaterials. Nat. Commun. 1, 21 (2010).
  * (19) B. Wood and J. B. Pendry. Metamaterials at zero frequency. J. Phys.: Condens. Matter 19, 076208 (2007).
  * (20) S. Narayana y Y. Sato. DC Magnetic Cloak. Adv. Mater. 24, 71–74 (2012).
  * (21) Fedor G$\ddot{o}$m$\ddot{o}$ry, Mykola Solovyov, Ján $\check{S}$ouc, Carles Navau, Jordi Prat-Camps and Álvaro Sanchez. Experimental Realization of a Magnetic Cloak. Science 335, 1466 (2012).
  * (22) Yungui Ma, Yichao Liu, Muhammad Raza, Yudong Wang y Sailing He. Experimental Demonstration of a Multiphysics Cloak: Manipulating Heat Flux and Electric Current Simultaneously. PRL 113, 205501 (2014).
  * (23) F. Y. Yang, K. Liu, K. Hong, D. H. Reich, P. C. Searson and C. L. Chien. Large Magnetoresistance of Electrodeposited Single-Crystal Bismuth Thin Films. Science 284, 1335 (1999).
  * (24) Environmental, Chemistry and Hazardous Materials News, Careers and Resources https://environmentalchemistry.com/
  * (25) Fedor G$\ddot{o}$m$\ddot{o}$ry, Mykola Solovyov, Ján $\check{S}$ouc, Carles Navau, Jordi Prat-Camps and Alvaro Sanchez. Supporting Online Material for Experimental Realization of a Magnetic Cloak. www.sciencemag.org/cgi/content/full/335/6075/1466/DC1