11institutetext: Department of Physics, P.O. Box 64, FI-00014 University of
Helsinki, Finland
11email<EMAIL_ADDRESS>22institutetext: Max Planck Institute for
Solar System Research, Justus-von-Liebig-Weg 3, D-37077 Göttingen, Germany
33institutetext: ReSoLVE Centre of Excellence, Department of Computer Science,
Aalto University, PO Box 15400, FI-00076 Aalto, Finland 44institutetext:
Leibniz-Institute for Astrophysics Potsdam, An der Sternwarte 16, 14482
Potsdam, Germany 55institutetext: Center of Excellence in Information Systems,
Tennessee State University, Nashville, TN 37209, USA 66institutetext:
Department of Physics and Astronomy, Uppsala University, Box 516, 751 20
Uppsala, Sweden

# Long-term spot monitoring of the young solar analogue V889 Her ††thanks:
Based on observations made with the Nordic Optical Telescope, operated on the
island of La Palma jointly by Denmark, Finland, Iceland, Norway, and Sweden,
in the Spanish Observatorio del Roque de los Muchachos of the Instituto de
Astrofisica de Canarias. Based on observations made with the HARPSpol
instrument on the ESO 3.6 m telescope at La Silla (Chile), under the programme
ID 091.D-0836.††thanks: Photometric data is only available in electronic form
at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via
http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/

T. Willamo 11 T. Hackman 11 J. J. Lehtinen 22 M. J. Käpylä, 2233 I. Ilyin 44
G. W. Henry 55 L. Jetsu 11 O. Kochukhov 66 N. Piskunov 66

(Received 2 November 2018 / Accepted 12 January 2019)

###### Abstract

Context. Starspots are important manifestations of stellar magnetic activity.
By studying their behaviour in young solar analogues, we can unravel the
properties of their magnetic cycles. This gives crucial information of the
underlying dynamo process. Comparisons with the solar cycle enable us to infer
knowledge about how the solar dynamo has evolved during the Sun’s lifetime.

Aims. Here we study the correlation between photometric brightness variations,
spottedness, and mean temperature in V889 Her, a young solar analogue. Our
data covers 18 years of spectroscopic and 25 years of photometric
observations.

Methods. We use Doppler imaging to derive temperature maps from high-
resolution spectra. We use the Continuous Period Search method to retrieve
mean V-magnitudes from photometric data.

Results. Our Doppler imaging maps show a persistent polar spot structure
varying in strength. This structure is centred slightly off the rotational
pole. The mean temperature derived from the maps shows an overall decreasing
trend, as does the photometric mean brightness, until it reaches its minimum
around 2017. The filling factor of cool spots, however, shows only a weak
tendency to anti-correlate with the decreasing mean brightness.

Conclusions. We interpret V889 Her to have entered into a grand maximum in its
activity. The clear relation between the mean temperature of the Doppler
imaging surface maps and the mean magnitude supports the reliability of the
Doppler images. The lack of correlation between the mean magnitude and the
spottedness may indicate that bright features in the Doppler images are real.

###### Key Words.:

Stars: activity – starspots – solar-type – individual: HD 171488

##  1 Introduction

The single, late-type star V889 Her (or HD 171488) is a young solar analogue
of spectral class G2 V, classified as a BY Draconis variable (Kazarovets &
Samus, 1997), with an estimated age of 30-50 Myr (Strassmeier et al., 2003).
Young, solar-type stars are important for the study of stellar activity and
evolution since they provide us a look into what the Sun might have been like
in the past. In this early evolutionary stage the stars are much more
magnetically active than the Sun today. An approximately 10-year activity
cycle has been detected from photometry of V889 Her, with a current trend of
decreasing brightness (Lehtinen et al., 2016). Many Doppler images of it have
been published, with a large high-latitude spot visible in all of them (e.g.
Strassmeier et al., 2003; Järvinen et al., 2008; Frasca et al., 2010).

Marsden et al. (2006) also studied the surface magnetic field of V889 Her
using Zeeman-Doppler imaging (ZDI) with Stokes I&V spectra. They found no
latitude dependence of the radial field, as opposed to the azimuthal field
forming ring-like structures around the pole. A similar result was reported in
the ZDI study by Jeffers & Donati (2008). Neither of these studies revealed
any clear connection between the radial magnetic field and cool spots.

Due to its young age, V889 Her is still a very rapid rotator, with a
rotational period $P_{\mathrm{rot}}\approx 1.33$ days (e.g. Järvinen et al.,
2008). The crucial parameters of stellar dynamos that drive the magnetic
activity are rotation and rotationally affected turbulent convection (see e.g.
the review of Brandenburg & Subramanian, 2005). The dynamo activity and the
resulting magnetic activity thus increase as a function of rotation, at least
to a certain limiting value (e.g. Wright et al., 2011; Wright & Drake, 2016).
In comparison to the Sun, the observations indicate that these stars exhibit
spot structures that have filling factors much larger than the Sun. The
resolution of the observations is not yet accurate enough to reveal whether
these large spotted areas are individual spots or consist of groups of smaller
spots. On rapid rotators, spots are commonly observed at high latitudes. This
has been confirmed even with interferometric observations on $\zeta$ And
(Roettenbacher et al., 2016).

In this paper we present both spectroscopic and photometric observations of
V889 Her, ranging for 18 and 25 years, respectively. From the spectroscopic
observations we calculate surface temperature maps of the star, using the
Doppler imaging method, and study the evolution of the mean temperature, and
the spot coverage and distribution. From the photometry we follow the long-
term brightness changes in the star, thus tracing its activity cycles.

##  2 Observations

###  2.1 Spectroscopy

Spectroscopic observations of V889 Her have been gathered using the SOFIN and
FIES instruments at the Nordic Optical Telescope (NOT) at La Palma and HARPS
at the ESO 3.6m telescope at La Silla. We have analysed 19 sets of spectra,
ranging from 1999 to 2017. The observations are listed in Table 1. Some of the
earlier sets have already been published by Järvinen et al. (2008). To
calculate the rotational phases we used their ephemeris

$\mathrm{HJD}_{0}=2449950.550+1.33697\times E,$ (1)

where HJD0 corresponds to phase zero.

The FIES observations were reduced with the standard FIEStool pipeline
(Telting et al., 2014), and the HARPS observations by the REDUCE package
(Piskunov & Valenti, 2002). The SOFIN observations were reduced with the new
SDS tool, which is described in more detail in Section 2.1.2. The spectral
resolutions of HARPS and FIES are 120 000 and 67 000, respectively.

####  2.1.1 SOFIN observations

Most of the spectroscopic observations were carried out with the SOFIN high-
resolution polarimetric spectrograph (Tuominen et al., 1999; Ilyin, 2000)
mounted at the NOT. We used the medium resolution optical camera with 1 m
focal length, which provides the resolving power of about 80 000 with the
entrance slit size of $0.5^{\prime\prime}$. The spectra were imaged with a
Loral $2048\times 2048\times 15\,\mu\mbox{m}$ CCD which provides a limited
wavelength coverage at the selected spectral range totaling 24 orders having
about 45 Å in length at 5500 Å.

The data in this paper were obtained with different spectral set-ups, which is
why different spectral regions were used for different seasons. A typical
exposure time for V889 Her is about 40 min, achieving a signal-to-noise ratio
of about 270 depending on the seeing and sky conditions, as well as air mass.

####  2.1.2 SOFIN data reduction

The SOFIN data reduction was done on a generic software platform written in
C++ on Linux called the Spectroscopic Data Systems (SDS). SDS can be used to
analyse data from different échelle spectrographs with the same tools. The
image processing pipeline specifically designed to handle SOFIN data
calibration flow and image specific content provides an automated process from
the beginning to the end product without human interaction and relies on
statistical inference and robust statistics. The software numerical toolkit
and graphical interface is designed upon previous software (Ilyin, 2000), and
its complete description is in preparation (Ilyin 2019, in prep.).

The image processing includes the following steps:

  * •

Bias subtraction using the CCD overscan and variance estimation of the pixel
intensity in the source images;

  * •

Master flat-field correction for the CCD spatial pixel-to-pixel noise. The
master flat is obtained as the normalised sum of 100 flat-field spectra with
the slit decker opened enough to overlap spectral orders;

  * •

Scattered light subtraction is achieved by a 2D smoothing spline fit to the
gaps between spectral orders of the science frames. The resistant statistics,
which has a higher breakdown point than robust statistics with regard to how
unaffected it is by outliers, is used to define the gaps;

  * •

The échelle orders are defined from the target image by a 3D fit of the
elongated Gaussian profiles along polynomial paths of the spectral orders;

  * •

The wavelength solution for the ThAr images uses a 3D Chebyshev fit of the
spectral line positions, and order number to their wavelengths. It also uses
resistant statistics to find the initial wavelength in the image. Typically,
all available 1200 ThAr lines are used in the fit which results in about 7 m/s
error of the fit in the image centre. The ThAr image is taken prior to the
science exposure during observations to alleviate bending and flexure of the
Cassegrain mounted spectrograph;

  * •

The optimal extraction of the spectral orders and cosmic spikes elimination in
the target image is done by constructing a spatial profile from the image
which is resampled and smoothed with 2D spline into a regular grid. The flux
in each pixel is obtained with a weighted linear least-squares fit of the
spatial profile at a given pixel to the data. The cosmic spikes elimination is
done prior to the fit with the use of robust statistics;

  * •

The wavelength calibration of the orders in the target image from the ThAr
wavelength solution is followed by the heliocentric correction of the
wavelength scale. We use NOVAS 2 (Kaplan et al., 1989) to perform this
correction as implemented in SDS;

  * •

The extracted images are subjected to a 2D smoothing spline fit to the
continuum. This is achieved by taking an average image of all observations
which follows by its continuum normalisation with the aid of the resistant
statistics to exclude spectral lines from the fit. The ratio of each
individual image and the normalised average, which removes most of the
spectral lines in the ratio, is then used to fit a smoothing 2D robust spline
which constitutes the intermediate continuum for the individual observation;

  * •

The final continuum normalisation is done during each iteration of the Doppler
imaging inversion by comparing the observed spectrum to the modelled spectrum.

####  2.1.3 Quality of the data sets

To estimate the reliability of the temperature inversions produced from our
observations, we have calculated the phase coverage $F$ of each set of spectra
by assuming that one spectrum covers 10% of one rotation, as was done by
Kochukhov et al. (2013). The phase coverage of our data sets vary from very
good (the whole rotation is covered in Jul 2005, i.e. $F=1.0$) to quite poor
(in Oct 1999 $F=0.325$); see Table 1 for a full list of the coverage. The
average signal-to-noise ratio ($S/N$) and the deviation in the final Doppler
imaging inversions ($\sigma$) are also shown there as an indicator of the
quality of each set of spectra.

###  2.2 Photometry

We analysed V-band photometry of V889 Her obtained between the years 1994 and
2018 with the Tennessee State University T3 0.4m Automatic Photoelectric
Telescope (APT) at Fairborn Observatory in southern Arizona.

The observations of V889 Her were made differentially with respect to the
comparison star HD 171286 and the check star HD 170829. The check minus
comparison differential magnitudes show that both stars are constant to better
than 0.005 mag from night to night and to 0.001 mag from year to year. The V
magnitude of HD 171286 (V=6.84; Oja, 1987) was added to the variable minus
comparison star differential magnitudes to derive the V magnitudes of V889
Her. See Lehtinen et al. (2016) for further information on the acquisition and
reduction of V889 Her observations with the T3 APT.

Table 1: Spectroscopic observations of V889 Her. Dates, HJD interval, number of spectra ($N$), esimated phase coverage ($F$), mean signal-to-noise ratio ($\langle S/N\rangle$), deviation between the Doppler imaging solution and the observations ($\sigma$) and instrument. Date | HJD-2400000 | $N$ | $F$ | $\langle S/N\rangle$ | $\sigma[\%]$ | Instrument
---|---|---|---|---|---|---
Jul – Aug 1999 | 51383.482 – 51394.511 | 11 | 0.596 | 252 | 0.455 | SOFIN
Sep 1999 | 51443.401 – 51449.387 | 6 | 0.434 | 199 | 0.743 | SOFIN
Oct 1999 | 51471.359 – 51475.337 | 4 | 0.325 | 285 | 0.400 | SOFIN
Jul – Aug 2001 | 52120.452 – 52131.470 | 12 | 0.595 | 234 | 0.457 | SOFIN
Aug 2002 | 52507.447 – 52515.429 | 8 | 0.552 | 226 | 0.490 | SOFIN
Jul – Aug 2004 | 53215.532 – 53228.485 | 11 | 0.723 | 239 | 0.786 | SOFIN
Jul 2005 | 53567.406 – 53577.562 | 16 | 1.0 | 250 | 0.511 | SOFIN
Aug – Sep 2006 | 53979.443 – 53991.449 | 6 | 0.476 | 225 | 0.651 | SOFIN
Jul 2007 | 54300.520 – 54311.576 | 15 | 0.950 | 235 | 0.500 | SOFIN
Sep 2008 | 54717.368 – 54723.486 | 8 | 0.703 | 228 | 0.680 | SOFIN
Aug – Sep 2009 | 55069.520 – 55081.520 | 14 | 0.606 | 196 | 0.760 | SOFIN
Jul 2010 | 55395.516 – 55408.642 | 10 | 0.707 | 208 | 0.682 | SOFIN
Aug – Sep 2012 | 56166.418 – 56174.413 | 9 | 0.458 | 260 | 0.598 | SOFIN
Aug 2013 | 56521.424 – 56525.446 | 5 | 0.408 | 147 | 0.744 | SOFIN
Sep 2013 | 56545.484 – 56551.581 | 7 | 0.479 | 210 | 0.723 | HARPS
Aug 2014 | 56881.493 – 56889.454 | 10 | 0.710 | 344 | 0.386 | FIES
Jul 2015 | 57204.511 – 57213.563 | 10 | 0.657 | 299 | 0.482 | FIES
Jun 2016 | 57555.553 – 57563.571 | 11 | 0.764 | 259 | 0.558 | FIES
Jun 2017 | 57908.615 – 57918.531 | 11 | 0.496 | 116 | 0.882 | SOFIN
Table 2: Stellar parameters of V889 Her. Parameter | Value | Reference
---|---|---
Spectral class | G2 V | Montes et al. (2001)
Age | 30-50 Myr | Strassmeier et al. (2003)
Period | 1.33697 d | Järvinen et al. (2008)
$v\sin i$ | 38.5 km/s | This work
Inclination | $60^{\circ}$ | Marsden et al. (2006)
Microturbulence | 1.6 km/s | Järvinen et al. (2008)
Macroturbulence | 3.0 km/s | Strassmeier et al. (2003)

##  3 Doppler imaging

Inhomogeneities in the surface temperature of a star, caused mainly by cool
spots, leave traces in photospheric absorption lines, which can be tracked if
the lines are sufficiently broadened by the star’s rotation. These line
profiles can be inverted to the surface temperature distribution which would
cause the observed line profiles. For this purpose we use a modified version
of the Doppler imaging code INVERS7, originally written by Piskunov et al.
(1990). For an optimal solution, INVERS7 uses Tikhonov regularisation to damp
high spatial frequencies of temperature gradients. Recent improvements in the
code include a molecular equilibrium solver (Piskunov & Valenti, 2017).

To reproduce the stellar spectra we use MARCS111http://marcs.astro.uu.se/
model atmospheres for temperatures in the range 3500-6500 K (Gustafsson et
al., 2008), combined with absorption line parameters obtained from the Vienna
Atomic Line Database222http://vald.astro.uu.se/ (VALD) (Piskunov et al., 1995;
Kupka et al., 1999). In the inversion we have used spectral regions including
seven dominant absorption lines: the Fe I line at 6265.1319 Å, the V I line at
6266.3069 Å, the Fe I line at 6411.6476 Å, the Fe I line at 6430.8446 Å, the
Fe II line at 6432.6757 Å, the Ca I line at 6439.0750 Å, and the Ni I line at
6643.6304 Å. Their line parameters are given in Table 3. For the Fe II and V I
lines and one of the Fe I lines we used the standard values for $\log(gf),$ as
found in the database, but for the two other Fe I lines and the Ca I and Ni I
lines we modified the value slightly, as shown in Table 3, in order to get
consistent line strengths in the different regions. We chose to adjust the
$gf$-values, although the inconsistency could be caused by different effects
related to difficulties in modelling stellar line profiles using standard LTE-
calculations. Adjusting the $gf$ -values is equivalent to changing the
strength of individual lines. If the Doppler image is calculated using just
one spectral line, the adjustment changes the mean temperature of the image,
but not the spot configuration. When using multiple lines, some of the
modelled lines may be too strong while others are too weak compared to the
observations. This will prevent reaching a satisfactory Doppler imaging
solution, if no adjustments are done. Since the aim here is not to calculate
absolute stellar parameters or element abundances, adjusting the $gf$-values
is acceptable.

In the SOFIN data the exact range of the available spectrum varies due to
different instrument setups, so the same lines could not be used for all data
sets. When available, we used lines 3-6, as listed in Table 3; when these
lines were not available, we used other combinations. In Aug 2013 line 6 was
not used, although available, due to obvious errors in some of the observed
phases. The lines used in each set are listed along with the results in Table
4. In the FIES and HARPS observations, the full visible spectrum is always
available, so we could always use lines 3-6. Line 6 is the only one available
in all sets for all instruments (although it was not used in Aug 2013), so we
tried to account for the effects of using different line combinations by also
repeating the inversions using only this Ca I line.

In addition to the strongest lines, we included 91 atomic lines and 607
molecular lines (mainly TiO and CN) from four different wavelength regions
(6260-6269 Å, 6405-6417 Å, 6425-6443 Å, and 6639-6648 Å) in the spectral line
synthesis. They are not strong enough to be distinguished individually, but
they can still affect the continuum level.

A telluric absorption line was cut away from the spectra around 6432-6433 Å.
The exact location of the line varies between the spectra because they have
been corrected for the Earth’s motion to the rest frame of V889 Her. This
overlaps heavily with the Fe II line at 6432.6757 Å, but the parts of the line
that remain after the removal of the atmospheric line are still useful.

Important stellar parameters for the inversion process are listed in Table 2.
These include the rotation period, $v\sin i$, the inclination, and micro- and
macroturbulence. Similarly to Järvinen et al. (2008), we have used solar
element abundances for V889 Her. As an initial guess for the temperature, we
used 5800 K.

There have been several studies where significant differential rotation on
V889 Her has been reported (Marsden et al., 2006; Jeffers & Donati, 2008), but
also studies that yielded lower values for it (Huber et al., 2009; Järvinen et
al., 2008) or no differential rotation at all (Kővári et al., 2011).
Observations (Reiners & Schmitt, 2003), theory (Kitchatinov & Rüdiger, 1999),
and global magnetoconvection simulations (Viviani et al., 2018) show that
strong differential rotation should not be common in rapid rotators, such as
V889 Her. Since the studies focusing on differential rotation yield very
inconsistent results, we have excluded it from our study, assuming the star to
rotate as a solid body. This is not necessarily true, but a moderate
differential rotation should not significantly affect the general picture of
spottedness, although it can affect the details of the spot configuration.
Nevertheless, our data can be successfully reproduced without differential
rotation.

Table 3: Parameters for the absorption lines used in the inversion.

# | Line | $\lambda(Å)$ | $\log(gf)$ | $\log(gf)_{\mathrm{standard}}$
---|---|---|---|---
1 | Fe I | 6265.1319 | -2.850 | -2.550
2 | V I | 6266.3069 | -2.290 | -2.290
3 | Fe I | 6411.6476 | -0.620 | -0.695
4 | Fe I | 6430.8446 | -2.106 | -2.106
5 | Fe II | 6432.6757 | -3.520 | -3.520
6 | Ca I | 6439.0750 | 0.450 | 0.390
7 | Ni I | 6643.6304 | -1.920 | -2.220
333$\log(gf)$ is the adopted value for the line (modified from the standard
value in some cases), while $\log(gf)_{\mathrm{standard}}$ is the standard
value retrieved from VALD.

##  4 Time series analysis

To analyse photometric data, we use the Continuous Period Search method (CPS)
developed by Lehtinen et al. (2011). The method works by dividing the data
into smaller sets and using a sliding window of constant time length. This
enables the light curve parameters to change between each (usually
overlapping) set. To each set, a function of the form

$y(t_{\mathrm{i}})=M+\sum_{k=1}^{K}\Big{(}B_{\mathrm{k}}\cos(k2\pi
ft_{\mathrm{i}})+C_{\mathrm{k}}\sin(k2\pi ft_{\mathrm{i}})\Big{)}$ (2)

is fitted. Here, $M$ is the mean magnitude of the set, $f=1/P_{\mathrm{rot}}$
is the rotational frequency, and $t_{\mathrm{i}}$ the time of the data point.
The best estimates for $M$, $f$, and the coefficients $B_{\mathrm{k}}$ and
$C_{\mathrm{k}}$ are retrieved for each set. The order $K$ of the model is
determined by using the Bayesian information criterion. An upper limit $K=2$
is used, so models with 0, 1, and 2 orders are tested. The length of one data
set is defined by a maximum time span $\Delta t_{\mathrm{max}}$, such that all
data points within the time $\Delta t_{\mathrm{max}}$ = 30 days of the first
data point of the set are included. Only sets with at least 14 data points are
analysed.

Here, we are primarily interested in the mean magnitude $M$. Periodic changes
in $M$ are interpreted as signs of an activity cycle. By using the mean
magnitude instead of individual measurements, the brightness variations caused
by the star’s rotation are diminished.

Jetsu et al. (2017) showed that the CPS method does not necessarily produce
correct rotation period estimates if there are multiple interfering
periodicities. However, this kind of interference will not affect the mean
magnitudes we use here, while it would affect other light curve parameters,
such as the amplitude. We have also extracted the phases of photometric
minima, which could be affected by interfering periodicities.

|
---|---
|
|
|
|
Figure 1: Doppler imaging maps of V889 Her for 1999-2008. Upper panels: Equirectangular projections of the surface temperature distribution. The longitude is given as the rotational phase. The vertical lines represent the phases of the individual spectra. Lower panels: Observed (black) and modeled (red) line profiles. |
---|---
|
|
|
|
Figure 2: Same as Fig. 1, but for 2009-2017.

##  5 Results

###  5.1 Spottedness

The surface temperature maps of V889 Her are presented in Figs. 1 and 2. From
these maps, the spot filling factor $f_{S}$, i.e. the fraction of stellar
surface covered with spots, has been calculated. We have defined all surface
elements with temperature $T\leq 5300$ K to be treated as spots. This chosen
threshold value is certainly higher than the typical temperature of sunspots,
but lower than the temperature of the unspotted surface of V889 Her (in our
Doppler maps the mean temperature is around 5750 K). The resulting filling
factors are shown in Table 4 and Fig. 3. It should be kept in mind, however,
that the choice of the limiting temperature for a spot is arbitrary, and the
resulting filling factors will vary greatly with small changes in this
threshold. A single data point may thus not be very reliable here; hence, we
need a longer time series where systematic trends can be seen, if we want to
follow the evolution of the spottedness. Furthermore, the filling factor may
be influenced by possible hot features.

A stable, large high-latitude spot, centred in most cases a bit off the pole,
is seen in all Doppler maps. The ones with very few rotational phases observed
(e.g. September and October 1999) or those where the phase coverage is poor
because the observed phases are close to each other (e.g. July-August 1999)
cannot be regarded as reliable as those with better phase coverage (see Table
1 for the phase coverage for each set). For instance in the maps from 1999 in
Fig. 1, the surface features such as bright areas are clearly distributed
around the observed phases in the longitudinal direction, indicating that they
are artefacts typical for insufficient phase coverage. It is also notable,
that with better phase coverage, the maps will in general be smoother. Thus,
most of the small-scale variations in the maps with poorer phase coverage are
probably also artefacts. The polar spot, however, is present in every map, and
thus is most probably a real feature. It has also been observed in previous
Doppler maps published of V889 Her (e.g. Strassmeier et al., 2003; Järvinen et
al., 2008; Frasca et al., 2010). The polar spot seems to be a stable feature,
lasting for decades but evolving in size.

####  5.1.1 Comparison with photometry

In the photometry, the clearly cyclic pattern has changed to a monotonic
decreasing trend for the time between about 2007 to 2017. The decreasing trend
was noted by Lehtinen et al. (2016), but the trend continued even more clearly
in the new data analysed here until the brightness reached its minimum and
turned to a rapid increase around 2017. The star seems thus to have reached a
state of higher activity, possibly something similar to the Modern Maximum,
the period from about 1950 to 2000 when the Sun was particularly active. It
will be interesting to see if the star will remain in this state or return to
its earlier activity state.

To compare photometric and spectroscopic results, both the filling factors
$f_{S}$ and mean magnitudes $M$ are plotted in Fig. 3. The hypothesis is that
with decreasing brightness the spottedness should increase since the
luminosity of rapidly rotating young stars should be dominated by dark spots
(Radick et al., 1990); however, this relation is not seen very clearly. In the
next section some possible reasons for this are discussed.

####  5.1.2 Reliability of $f_{S}$

Inconsistencies in the spot filling factors may be caused by the change in
spectral lines used. The level of $f_{s}$ changes between 2007 and 2008 from
around 0.05 to 0.1. This could be explained by the different spectral lines
rather than any real changes in the star since the preferred lines (3-6 in
Table 3) could not be used in the SOFIN sets from 2008 onward. However, in the
data from 2004 and 2006 the Fe I and Fe II lines between 6411 and 6433 Å
(lines 3-5) were not available, and still $f_{S}$ is on a level similar to
those in the other data sets prior to 2008.

The FIES and HARPS observations between 2013 and 2016, where the full visible
spectrum was available, yield $f_{S}$ values comparable to those from the sets
before 2008. The preferred lines 3-6 have been used here, and the results
differ notably from the SOFIN observations between 2008 and 2017, where
different lines had to be used. In Fig. 3, all SOFIN sets where lines 3-6
could be used are highlighted, and most of the other SOFIN sets clearly differ
from these.

As a test of the effect of using different lines in different data sets, we
repeated the temperature inversion for all maps using only the Ca I line at
6439.0750 Å (line 6), which is the only line present in all the data sets. It
is worth mentioning that as this is the strongest absorption line, it is
formed highest up in the photosphere and may be affected by non-LTE effects
more than the other lines, which makes it more difficult to model. Using only
a single line causes more noise, but should reduce systematic differences
related to different spectral regions. In Fig. 4 we compare the filling
factors from Fig. 3 to those where only the Ca I line was used. The data from
the different instruments are shown separately.

Using only the Ca I line, the rise in the filling factors around 2008 seems
smoother. There is also more variance in the data points produced with only Ca
I, as is expected. In the data point from Aug 2013, $f_{S}$ exceeds 0.3, which
is a clear outlier. In this particular set, we excluded the Ca I line from our
combined line list, due to obvious instrumental errors.

The fact that no consistent trend is found from the filling factors could be
due to the effects of using different line combinations, which yields a higher
noise than the true changes in the spottedness. In the three FIES
observations, which should be comparable to each other, however, a clear
rising trend of $f_{S}$ is seen. The HARPS data from Sep 2013 also agree with
this trend. The SOFIN observations, unfortunately, do not seem to be directly
comparable as they clearly show a higher spottedness.

The differences in $f_{S}$ in Fig. 3 could be explained solely from the use of
different lines, but in Fig. 4 there is a notable difference between the SOFIN
observations during and before 2013 and the following HARPS and FIES
observations, even when only the Ca I line was used. In principle there should
be no differences due to instrumentation.

We calculated a Pearson correlation coefficient for the mean magnitudes and
$f_{S}$ values retrieved from the combined lines as $r=0.18$, which implies a
weak correlation. When using only those $f_{S}$ values where the preferred
lines 3-6 could be used , the correlation is much stronger; $r=0.48$. This
indicates that the choice of lines may significantly affect the results. The
correlation coefficient between mean magnitudes and $f_{S}$ derived only with
the Ca I line is $r=0.24$.

To test for the effect of the choice of limiting temperature for a spot, we
also calculated $f_{S}$ by using the limiting temperature 5000 K and 4700 K
(instead of 5300 K) for a spot. This resulted in $r=0.06$ for 5000 K and
$r=0.20$ for 4700 K. These differences (as $r=0.18$ with 5300 K) further imply
that the values of $f_{S}$ may not be that reliable.

To calculate the correlation coefficient, we linearly interpolated the mean
magnitude of the epoch of the spectroscopic data points from the two closest
phtometric data points. Although $f_{S}$ does not show a very clear
correlation to photometry in this case, we point out that there are other
stars for which the correlation seems to be much clearer (e.g. hackman2012;
hackman2018).

###  5.2 Mean temperature

We also calculated the mean temperature of each Doppler map. Large spots
should have a cooling effect on the temperature, and thus also lower the
luminosity of the star. The photometric mean magnitudes are plotted against
the spectroscopic mean temperatures in Fig. 5. A much clearer correlation is
seen here, suggesting that the decreasing luminosity of V889 Her is indeed
caused by decreasing surface temperature, most likely in the form of
increasing spottedness, although this is not clearly seen from the spot
filling factors. The Doppler imaging method is very sensitive to changes in
the mean temperature because this will affect the strength of the disk
integrated line.

We calculated a Pearson correlation coefficient for the mean magnitude and
mean temperature values as $r=-0.60$. This is a significantly stronger
correlation than between mean magnitudes and $f_{S}$. Thus, the mean
temperature seems to be a more reliable indicator of the evolution of
photospheric activity than the spot coverage, at least in this case. This is
important since the mean temperature is independent of spot models or
definitions of spotted regions, and is affected even by small features, which
go easily undetected in the Doppler images and spot filling factors.
Furthermore, the mean temperature is much less sensitive to Doppler imaging
artefacts since they tend to be alternating cool and hot features which will
eliminate each other in the mean temperature.

Figure 3: Mean magnitudes of each independent set in the CPS analysis (blue
dots) and spot filling factors (red stars). The sets where the preferred lines
at 6411-6439 Å were used are highlighted with a circle.  Figure 4: Spot
filling factors calculated with combined lines, compared to those where only
the Ca I line was used. SOFIN, FIES, and HARPS observations are distinguished
from each other.  Figure 5: Mean magnitudes (blue dots) plotted against the
mean temperature of each spectral set (red stars).

###  5.3 Distribution of spots

In Figs. 6 and 7 we show the averaged surface temperature distribution of V889
Her over each latitude and longitude, respectively, to study changes in the
average spot locations. The latitudinal spot distribution could reveal effects
resembling the butterfly diagram of the Sun, whereas the longitudinal
distribution could reveal active longitudes. However, no clear trends are
revealed from our data. In Fig. 6 the mean spot latitude remains fairly
constant.

In Fig. 7, there is some evolution in the longitudinal spot distribution, but
no clear trends. It should also be kept in mind that possible differential
rotation, which is not taken into account here, can affect the apparent
longitudes of spots.

Figure 7 also shows the phases of the photometric minima, which should
coincide with the dominating spot regions. This can be used as an independent
check of the reliability of the spot longitudes. There is certainly not a
perfect match between them, but in many cases they are located close to each
other, suggesting that the spot longitudes are still fairly reliable.
Photometric minima are also concentrated in most cases close to phase 0,
although there is a great deal of scatter. Jetsu (2018) showed that phases of
photometric minima are accurate only if there is one dominating spot region.

Figure 6: Surface temperature distribution of V889 Her from our Doppler images averaged over all longitudes for each latitude. The time intervals of the observations are extended for clearer plots.  Figure 7: Surface temperature distribution of V889 Her from our Doppler images averaged over all latitudes for each longitude. The time intervals of the observations are extended for clearer plots. Phases of the photometric minima are also shown with their error bars.  Table 4: Filling factors $f_{S}$ and mean temperatures $T_{\mathrm{mean}}$ for each set of spectral observations. The lines used in each set are numbered as in Table 3. Set | Lines | $f_{S}$ |  $T_{\mathrm{mean}}$ [K]
---|---|---|---
Jul – Aug 1999 | 3,4,5,6 | 0.0725 | 5796
Sep 1999 | 3,4,5,6 | 0.0606 | 5771
Oct 1999 | 3,4,5,6 | 0.0318 | 5809
Jul – Aug 2001 | 3,4,5,6 | 0.0489 | 5774
Aug 2002 | 3,4,5,6 | 0.0455 | 5743
Jul – Aug 2004 | 1,2,6 | 0.0523 | 5857
Jul 2005 | 3,4,5,6 | 0.0515 | 5790
Aug – Sep 2006 | 1,2,6 | 0.0693 | 5764
Jul 2007 | 3,4,5,6 | 0.0504 | 5776
Sep 2008 | 1,2,6,7 | 0.1329 | 5707
Aug – Sep 2009 | 1,2,6,7 | 0.1067 | 5778
Jul 2010 | 1,2,6,7 | 0.1448 | 5721
Aug – Sep 2012 | 1,2,6,7 | 0.1228 | 5743
Aug 2013 | 1,2,7 | 0.1018 | 5704
Sep 2013 | 3,4,5,6 | 0.0421 | 5752
Aug 2014 | 3,4,5,6 | 0.0550 | 5725
Jul 2015 | 3,4,5,6 | 0.0590 | 5706
Jun 2016 | 3,4,5,6 | 0.0706 | 5697
Jun 2017 | 1,2,6 | 0.1329 | 5626

###  5.4 Chromospheric activity

As an additional check for the activity level of V889 Her, we studied the
evolution of the H$\alpha$ line at 6562.8 Å. The line is seen as an absorption
line, but magnetic activity will increase H$\alpha$ emission from the
chromosphere, which will weaken the strength of the absorption line. None of
the SOFIN spectra, unfortunately, contain the line nor any of the other
significant chromospheric lines, such as the H$\beta$ or Ca II H&K lines, but
from our FIES and HARPS spectra we get four data points for the average
strength of the line. The equivalent width (EW), characterising the strength
of a spectral line, is defined as

$EW=\int_{\lambda_{1}}^{\lambda_{2}}\big{(}1-\frac{F_{\lambda}}{F_{c}}\big{)}\mathrm{d}\lambda\approx\sum_{\lambda=6560\AA}^{6566\AA}\big{(}1-\frac{F_{\lambda}}{F_{c}}\big{)}\Delta\lambda,$
(3)

where $F_{\lambda}$ is the intensity at wavelength $\lambda$; $F_{c}$ is the
continuum intensity, which is here normalised to 1; and $\Delta\lambda$ is the
interval between wavelength points. The line strength is integrated over the
line, which we have defined as the region between $\lambda_{1}=6560$ Å and
$\lambda_{2}=6566$ Å. We have calculated this for every spectrum in each set,
and taken the average value as the indicator for the chromospheric activity of
the set. The EW values are shown in Table 5.

As we have not been able to use SOFIN spectra to check for chromospheric
emission, we have too few data points to draw any reliable conclusions. Still,
the EW from the FIES and HARPS spectra behave as would be expected: as the
brightness of the star decreases from 2013 to 2016, the EW also decreases. As
photospheric absorption dominates the H$\alpha$ line, the decrease in EW
suggests that chromospheric activity is increasing. This supports the
conclusion that V889 Her has indeed entered a state of enhanced activity.

Table 5: Equivalent widths of the H$\alpha$ line for the FIES and HARPS spectra. Set | EW [Å]
---|---
Sep 2013 | 0.657
Aug 2014 | 0.640
Jul 2015 | 0.574
Jun 2016 | 0.538

##  6 Conclusions

We make the following conclusions from our study:

  1. 1.

The clear cyclic behaviour of V889 Her seems to have been replaced with an
almost monotonic luminosity decrease between 2007 and 2017 when the star
reached its activity maximum, after which its luminosity increased rapidly,
although it is still notably lower than before 2007. We interpret it to be in
a state of enhanced activity, possibly a grand maximum.

  2. 2.

The correlation between brightness and mean temperature is much clearer than
the anti-correlation between brightness and spottedness.

  3. 3.

The polar spot, also reported in previous Doppler imaging studies of V889 Her,
is a fairly stable feature. It does not show any major changes even with the
decreasing brightness and mean temperature.

From our results, it seems clear that the choice of which spectral lines are
used in the Doppler imaging inversion may affect the outcome of the
temperature map. This may cause serious problems when using instruments for
which the spectral region varies between different observing seasons. There
may also be some effects that depend on the instrument.

Our results indicate, that V889 Her has entered a state of decreased
brightness and mean temperature. The cyclic behaviour seems to have continued
after the long luminosity decrease, as the brightness has again increased in
the most recent photometry. The fact that the spottedness did not follow the
same monotonic evolution may indicate that bright surface features also
influence the mean magnitude. Bright spots in Doppler images based on
insufficient data are, of course, not reliable. However, the maps based on
high $S/N$ and good phase coverage also show regions with surface temperatures
slightly higher than expected for the unspotted surface. These hotter surface
areas could be generated by magnetic activity, but could also be of non-
magnetic origin. Käpylä et al. (2011) showed that both cooler and hotter
regions may arise from large-scale vortices generated by rotating turbulent
convection. With rotation rates comparable to those of V889 Her, hotter
anticyclonic vortices would be preferred from this mechanism.

###### Acknowledgements.

The work of T.W. was financed by the Emil Aaltonen Foundation. M.J.K. was
partially supported by the Academy of Finland Centre of Excellence ReSoLVE
272157. The Automated Astronomy Program at Tennessee State University has been
supported by NASA, NSF, TSU, and the State of Tennessee through its Centers of
Excellence programme. O.K. acknowledges support by the Knut and Alice
Wallenberg Foundation, the Swedish Research Council, and the Swedish National
Space Board. This work has made use of the VALD database, operated at Uppsala
University, the Institute of Astronomy RAS in Moscow, and the University of
Vienna. The authors thank the anonymous referee for the good comments and
suggestions which improved the manuscript.

## References

  * Brandenburg & Subramanian (2005) Brandenburg, A. & Subramanian, K. 2005, Phys. Rep, 417, 1
  * Frasca et al. (2010) Frasca, A., Biazzo, K., Kővári, Z., Marilli, E., & Çakırlı, Ö. 2010, A&A, 518, A48
  * Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951
  * Huber et al. (2009) Huber, K. F., Wolter, U., Czesla, S., et al. 2009, A&A, 501, 715
  * Ilyin (2000) Ilyin, I. V. 2000, PhD thesis, Astronomy Division Department of Physical Sciences P.O.Box 3000 FIN-90014 University of Oulu Finland
  * Järvinen et al. (2008) Järvinen, S. P., Korhonen, H., Berdyugina, S. V., et al. 2008, A&A, 488, 1047
  * Jeffers & Donati (2008) Jeffers, S. V. & Donati, J.-F. 2008, MNRAS, 390, 635
  * Jetsu (2018) Jetsu, L. 2018, ArXiv e-prints [arXiv:1808.02221]
  * Jetsu et al. (2017) Jetsu, L., Henry, G. W., & Lehtinen, J. 2017, ApJ, 838, 122
  * Kaplan et al. (1989) Kaplan, G. H., Hughes, J. A., Seidelmann, P. K., Smith, C. A., & Yallop, B. D. 1989, AJ, 97, 1197
  * Käpylä et al. (2011) Käpylä, P. J., Mantere, M. J., & Hackman, T. 2011, ApJ, 742, 34
  * Kazarovets & Samus (1997) Kazarovets, E. V. & Samus, N. N. 1997, Information Bulletin on Variable Stars, 4471
  * Kővári et al. (2011) Kővári, Z., Frasca, A., Biazzo, K., et al. 2011, in IAU Symposium, Vol. 273, Physics of Sun and Star Spots, ed. D. Prasad Choudhary & K. G. Strassmeier, 121–125
  * Kitchatinov & Rüdiger (1999) Kitchatinov, L. L. & Rüdiger, G. 1999, A&A, 344, 911
  * Kochukhov et al. (2013) Kochukhov, O., Mantere, M. J., Hackman, T., & Ilyin, I. 2013, A&A, 550, A84
  * Kupka et al. (1999) Kupka, F., Piskunov, N., Ryabchikova, T. A., Stempels, H. C., & Weiss, W. W. 1999, A&AS, 138, 119
  * Lehtinen et al. (2011) Lehtinen, J., Jetsu, L., Hackman, T., Kajatkari, P., & Henry, G. W. 2011, A&A, 527, A136
  * Lehtinen et al. (2016) Lehtinen, J., Jetsu, L., Hackman, T., Kajatkari, P., & Henry, G. W. 2016, A&A, 588, A38
  * Marsden et al. (2006) Marsden, S. C., Donati, J.-F., Semel, M., Petit, P., & Carter, B. D. 2006, MNRAS, 370, 468
  * Montes et al. (2001) Montes, D., López-Santiago, J., Gálvez, M. C., et al. 2001, MNRAS, 328, 45
  * Oja (1987) Oja, T. 1987, A&AS, 71, 561
  * Piskunov et al. (1990) Piskunov, N., Tuominen, I., & Vilhu, O. 1990, A&A, 230, 363
  * Piskunov & Valenti (2017) Piskunov, N. & Valenti, J. A. 2017, A&A, 597, A16
  * Piskunov et al. (1995) Piskunov, N. E., Kupka, F., Ryabchikova, T. A., Weiss, W. W., & Jeffery, C. S. 1995, A&AS, 112, 525
  * Piskunov & Valenti (2002) Piskunov, N. E. & Valenti, J. A. 2002, A&A, 385, 1095
  * Radick et al. (1990) Radick, R. R., Lockwood, G. W., & Baliunas, S. L. 1990, Science, 247, 39
  * Reiners & Schmitt (2003) Reiners, A. & Schmitt, J. H. M. M. 2003, A&A, 412, 813
  * Roettenbacher et al. (2016) Roettenbacher, R. M., Monnier, J. D., Korhonen, H., et al. 2016, Nature, 533, 217
  * Strassmeier et al. (2003) Strassmeier, K. G., Pichler, T., Weber, M., & Granzer, T. 2003, A&A, 411, 595
  * Telting et al. (2014) Telting, J. H., Avila, G., Buchhave, L., et al. 2014, Astronomische Nachrichten, 335, 41
  * Tuominen et al. (1999) Tuominen, I., Ilyin, I., & Petrov, P. 1999, in Astrophysics with the NOT, ed. H. Karttunen & V. Piirola, 47
  * Viviani et al. (2018) Viviani, M., Warnecke, J., Käpylä, M. J., et al. 2018, A&A, 616, A160
  * Wright & Drake (2016) Wright, N. J. & Drake, J. J. 2016, Nature, 535, 526
  * Wright et al. (2011) Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, ApJ, 743, 48


#  ILASS-Americas 31st Annual Conference on Liquid Atomization and Spray
Systems, May 2021

Improved Methods for Mixing-Limited Spray Modeling

Majid Haghshenas111Corresponding Author<EMAIL_ADDRESS>Mechanical and
Industrial Engineering Department, University of Massachusetts Amherst, MA
Peetak P. Mitra  Mechanical and Industrial Engineering Department, University
of Massachusetts Amherst, MA  Chu Wang  Convergent Science Inc., Madison, WI
Fabien Tagliante  Sandia National Laboratories, Livermore, CA  Lyle Pickett
Sandia National Laboratories, Livermore, CA  David P. Schmidt  Mechanical and
Industrial Engineering Department, University of Massachusetts Amherst, MA

( Abstract
The realization that interfacial features play little role in diesel spray
vaporization and advection has changed the modus operandi for spray modeling.
Lagrangian particle tracking has typically been focused on droplet behavior,
with sub-models for breakup, collision, and interfacially-limited
vaporization. In contrast, the mixing-oriented spray models are constructed so
that gas entrainment is the limiting factor in the evolution of momentum and
energy. In this work, a new spray model, ELMO (Eulerian Lagrangian Mixing-
Oriented), is implemented in a three-dimensional CFD code with two-way
coupling with the gas phase. The model is verified and validated with three
canonical sprays including spray A, H, and G from the ECN database. )

## Introduction

Spray modeling in CFD has long been dominated by drop-oriented Lagrangian-
Eulerian methods, where the stochastic nature of the algorithm provides
efficient sampling of a sub-set of the actual drops [1] without the
frustrations of numerical diffusion [2]. The solution of the spray, which is
statistically represented as parcels, depends on the droplet breakup models,
drag models, collision models, and turbulent dispersion models. The essence of
the spray model was drop-oriented, while there could be partial consideration
for the presence of neighboring drops in dense sprays [3]. This method was
successful under many situations, and for this reason, Lagrangian-Eulerian
modeling, such as Kelvin-Helmholtz Rayleigh Taylor (KH-RT), remains a mainstay
of spray modeling with CFD. One major drawback of these models is that there
are many adjustable constants that need to be tuned, each of which have an
effect on model predictions. For example, $B_{0}$, $B_{1}$, $C_{t}$ and
$C_{RT}$ are the KH droplet size constant, KH breakup time constant, RT
droplet size and RT breakup time constant respectively.

Experimental studies conducted by Siebers et al. [4, 5, 6] contradict the
basis of past spray models such as KH-RT. Siebers [4, 5, 6] concluded ’the
processes of atomization and the ensuing interphase transport of mass and
energy at droplet surfaces are not limiting steps with respect to fuel
vaporization in Direct Injection (DI) diesel sprays but rather they are
limited by the turbulent mixing between phases’. Moreover, a recent numerical
study on spray evolution by Agarwal and Trujillo [7] questioned the validity
of common spray models’ assumed linking linear stability with primary
atomization.

Siebers [4] proposed a stand-alone steady jet one-dimensional simplified model
for predicting liquid length for non-vaporizing sprays. Other studies,
including Pastor et al. [8, 9], Desantes et al. [10, 11] and Musculus et al.
[12] extended the original formulation using the Eulerian control volume
approach to add the physics of vaporization, model reacting flows, and include
the end of injection transients with great success. These models assumed a
radial profile of the spray contrary to the assumptions in the original work
[10, 8]. This is consistent with experimental and high-fidelity simulation
observations for sprays. The typical profile for a spray during its evolution
is shown in Figure 1. The spray starts off as a slug flow near the nozzle
exit, and as it progresses downstream of the nozzle exit, acquires a radial
profile due to aerodynamic interactions. Here $\alpha$ is the shape factor
that allows the profiles to evolve from a uniform shape at the nozzle
($\alpha=\infty$) to a fully developed shape ($\alpha=1.5$) far downstream
[12].

Figure 1: The spray profile behavior from nozzle exit to downstream. $U_{C}$
is centerline velocity and $R$ is radius of the spray profile.

Most of the mixing-oriented models were formulated to be used only as stand-
alone models and not be coupled with a multidimensional CFD solver. There have
been several implementations of reduced-order Mixing-limited models in CFD
codes in recent years [13, 14, 15, 16], however, these promising works assumed
a constant spray angle and constrained the mass and momentum exchange to be
mostly predetermined. The Eulerian Lagrangian Mixing Oriented (ELMO) model
presented here is a fully coupled model that is well-suited to time-varying
spray angle. The inputs of the ELMO model are informed by the internal flow
and experimental observations, including discharge rate and spray angle. The
near-nozzle spray evolution is based on the more accurate mixing-limited
hypothesis inspired by previous efforts in this area [4, 10, 8, 12]. ELMO is
coupled to a CFD solver, and in the dilute region downstream of the nozzle
this model transitions to a standard Lagrangian Eulerian treatment with its
droplet centric formulation. A schematic is shown in Figure 2.

Figure 2: The ELMO model idealization is shown in this figure. In the near
nozzle region, the ELMO operates under the mixing limited assumption. When the
spray becomes sufficiently dilute, the ELMO transitions to standard LE and its
droplet centric approach.

The advantages of this approach include: 1) The model can benefit from
sophisticated transient internal models, informing the inputs for the external
spray modeling paradigm. 2) The near nozzle behavior is consistent with the
modern mixing limited hypothesis, therefore expected to be more accurate. 3)
In the downstream region, the model employs conventional droplet models
(collisions, drag) where they work best [17]. 4) Since the spray angle is an
input to the model, it can be more effectively modeled than in the current
state-of-the-art. Siebers [4] shows the spray angle to be transient during
high/low needle lift. Most of the current modeling approaches ignore this
transient, yet important, behavior.

The Musculus-Kattke model [12] is used as the starting point for the numerical
formulation for the ELMO model. While the previous implementation is in an
Eulerian frame of reference, the ELMO model is in the corresponding Lagrangian
frame of reference. This essentially provides the capability to use the
Lagrangian ’capsules’ as self-contained liquid fuel parcels that entrain
ambient air during the translation process along the spray axis. The capsule
involves a moving Lagrangian parcel over an underlying Eulerian mesh. This
capsule therefore has two phases, the liquid fuel and gas, including ambient
air as well as vaporized fuel due to vaporization caused by the entrainment of
hot ambient air. The in-capsule vaporization is modeled using the Desantes
[10] formulation based on equilibrium evaporative physics. Therefore, the ELMO
model is built on assumptions of inertial and thermal equilibrium, limited
only by the entrainment rate of surrounding ambient air.

In the following sections, first, the numerical methodology of ELMO method is
presented. Next, the model is verified against Musculus-Kattke and Desantes
models for kinematics and thermodynamics respectively. And finally, validation
results are presented based on standard test cases from Engine Combustion
Network [18], including spray A, H, and G.

## Methodology

In this section we will introduce the detailed discussion about the numerical
methods for the formulation of the ELMO model, including some of the
assumptions.

### The capsule formulation

The Eulerian formulation in Musculus Kattke [12], while suitable for a
standalone spray prediction model, is converted into its corresponding
Lagrangian formulation. This is done to use the available capabilities of the
underlying CFD code CONVERGE [19, 20] in handling Lagrangian parcels. The ELMO
model has a few assumptions:

  * •

The mass of the fuel in a capsule remains constant throughout the evolution of
the transient, up to the point of transition. It can be represented
mathematically as:

$\frac{dm_{f}}{dt}=0$ (1)

where $m_{f}$ is the fuel mass

  * •

The momentum of the capsule remains constant. This translates to the capsule
velocity slowing down in proportion to the mass of the ambient air entrained.
Mathematically, at an individual capsule level this can be represented as:

$\frac{dM}{dt}=0$ (2)

where $M$ is the momentum

  * •

Each capsule starts off with its spray angle at birth. This spray angle
dictates the shape and entrainment within the capsule and remains constant
throughout the capsule’s life, even as the spray angle of injection changes
over time.

  * •

Underlying all of the equations (introduced later in this section) are
assumptions of equilibrium. Specifically, the phases all move at the same
velocity and experience equal pressure. All phases are in thermal equilibrium
as well. The capsule is in thermodynamic equilibrium, with an assumption of
locally homogeneous flow. Therefore, the liquid and gas are at the same
temperature.

Applying the Reynolds Transport Theorem (RTT), we convert the Eulerian MK
control volumes [12] into the Lagrangian ELMO capsules.

$\frac{dm_{f}}{dt}=\frac{\partial}{\partial
t}\int_{V}\rho_{f}\bar{X_{f}}dV+\int_{CS}\rho_{f}\bar{X_{f}}(\vec{u}-\vec{u_{r}})dA$
(3)

By design, the left hand side of the equation is set to zero. The volumetric
integral time is also zero due to fuel mass conservation principles, which
reduces the above equation to

$\int_{FS}\bar{X_{f}}\vec{u}.dA-\int_{FS}\bar{X_{f}}\vec{u_{r}}\cdot dA=0$ (4)

which further can be reduced to:

$\vec{u_{r}}=\beta\vec{u}$ (5)

where $m_{f}$ is the mass of fuel, $u$ is the axial velocity, $u_{r}$ is the
control surface velocity, $\beta$ is the shape factor of fuel volume fraction,
$\rho_{f}$ is the density of fuel, $\bar{X_{f}}$ is the fuel volume fraction,
$CS$ is the control surface and $FS$ is the front surface. The term $\beta$
accounts for the shape of the fuel volume fraction and velocity profiles,
ranging from $\beta=1$ for a uniform profile to approximately $\beta=2$ for
fully developed jets and is based on the representation shown in Figure 1 for
$\alpha$. For details about the shape factor term, please refer the original
source of this work [12].

Similarly, applying the RTT, the momentum formulation can be represented as:

$\frac{dM}{dt}=\frac{\partial}{\partial
t}\int_{V}\bar{\bar{\rho}}\vec{u}dV+\int_{CS}\bar{\bar{\rho}}\bar{\bar{u}}(\vec{u}-\vec{u_{r}}.n)dA$
(6)

per the design of the capsule concept and in the absence of a pressure
gradient, the equation above reduces (approximately) to

$\frac{\partial}{\partial t}\bar{\bar{u}}m=0$ (7)

where $M$ is the momentum of the capsule, $\bar{\bar{\rho}}$ is the mixture
density of the capsule, $\vec{u}$ is the axial velocity, $\bar{\bar{u}}$ is
the cross-sectionally averaged axial velocity, $u_{r}$ is the control surface
velocity and $m$ is the total mass of capsule.

Since the fuel mass in a capsule remains constant, it can either exist in a
liquid or vapor (gas) phase. The growth in volume (and mass) of the capsule at
successive time instants is accounted by the entrainment of the ambient air
(schematic shown in Figure 3). Therefore the mass and volume conservation
equations have the following functional form:

$m=m_{l}+m_{v}+m_{a}$ (8)

$V=V_{l}+V_{v}+V_{a}$ (9)

where the subscripts $l$, $v$, and $a$ are liquid fuel, vapor, and non-
condensible gas respectively. The presence of the non-condensible gas in the
capsule necessitates calculating phase interactions with the underlying
Eulerian CFD control volume. This approach will be detailed in the coming
sections.

Figure 3: The schematic of the ambient air entrainment in the capsule.

The increase in capsule mass due to the entrainment of air can be represented
as

$m^{t+1}_{capsule}=m^{t}_{capsule}+\Delta m_{air}$ (10)

In the present implementation, the effect due to entrained momentum is
ignored, similar to the original work [12]. This therefore leads to a constant
momentum for each capsule, with decelerating capsule velocity. Like the
implementation of Musculus and Kattke, the consideration of a radial velocity
and density profile is included here, which has implications for momentum.
This radial profile is manifested by the term $\beta$ such that the velocity
is higher along the spray axis, where most of the mass is concentrated. As the
capsule advances in time, the position is updated by a simple explicit Euler
time integration given by:

$z^{t+1}=z^{t}+\beta\bar{u}\Delta t$ (11)

where $z$ is the axial location of the capsule, and $\bar{u}$ is the capsule
velocity at the location $z$. While the capsule volume increases radially, its
axial extent diminishes as the control surfaces of the capsules are advected
at different velocities based on their location along the spray axis. This
same reduction in velocity represents a deceleration of the center of mass of
each capsule. The derivation of the front side (fs) and back side (bs)
location of each capsule begins with the assumption that these capsule
boundaries lie halfway between the centroid $z$ of a capsule at its current
position and the next position for the front side or halfway between the
current centroid location and previous location for the back side. Expressing
this mathematically:

$fs^{t}=\frac{z^{t}+z^{t+1}}{2}$ (12)

$bs^{t}=\frac{z^{t}+z^{t-1}}{2}$ (13)

In order to make the train of capsules continuous, we constrain the back
surface of the leading capsule to match with the front surface of the trailing
capsule.

Once the locations of the front and back sides are updated, the geometry of
the capsule is fully determined and the kinematics of the advection process is
complete. Once the axial extents of the capsules are known, the other features
such as velocity, volume can be easily calculated.

### Thermodynamic formulation

To model the process of vaporization, we follow the analysis first proposed by
Desantes et al.[10]. Desantes et al. extended a one-dimensional non-
evaporative spray model to take into account the effects of evaporation under
direct-injection diesel engine conditions. Their approach is based on the
assumption that evaporation is limited by the mixing process between fuel and
ambient gas, and makes use of non-ideal gas state relationships for the
description of the whole spray. Their model has been shown to accurately
predict the influence of fuel type and both ambient and injection conditions
on liquid spray penetration. The hot entrained ambient air is the primary
driver of the vaporization process. The local thermodynamic equilibrium in the
capsule depends on the mass fractions. Mathematically, they are represented as
below.

$Y_{f}=Y_{f,l}+Y_{f,v}$ (14)

$Y_{a}=1-Y_{f}$ (15)

where suffixes $f$, $a$, $l$, $v$ correspond to fuel, air, liquid, and vapor
phase respectively. While in the original work, non-ideal state relationships
were considered, in this implementation we assume all gases to be ideal. This
greatly simplifies our calculations without appreciable differences between
predictions. The local mixture density is calculated using an ideal mixture
assumption

$\rho=\frac{1}{\frac{Y_{f,l}}{\rho_{f,l}}+\frac{Y_{f,v}}{\rho_{f,v}}+\frac{Y_{a}}{\rho_{a}}}$
(16)

where $Y_{i}$ is the mass fraction of the mixture $i$ component and $\rho_{i}$
is the density for the pure component $i$ at the mixture temperature $T$ and
total pressure $P$. The local mixture enthalpy $h$ is also calculated through
an ideal mixture assumption as

$h(T)=Y_{f,l}\cdot h_{f,l}(T)+Y_{f,v}\cdot h_{f,v}(T)+Y_{a}\cdot h_{a}(T)$
(17)

where $Y_{i}$ is the mass fraction of the mixture $i$ component and $h_{i}$ is
the enthalpy for the pure component $i$ at the mixture temperature $T$.

Due to the liquid-vapor equilibrium, the composition and temperatures of
phases are coupled i.e., for a given $Y_{f}$, the proportion of liquid or
evaporated fuel is not known, nor is the mixture temperature. It can however
be solved iteratively, using the formulation below

$\begin{split}&\frac{Y_{f,v}}{1-Y_{f}}=\\\
&=\frac{h_{a}(T_{a,\inf}-h_{a}(T)-\frac{Y_{f}}{1-Y_{f}}[h_{f,l}(T)-h_{f,l}(T_{f,0})]}{\Delta
h_{v}(T)}\end{split}$ (18)

In this equation, if $Y_{f}$ is known the right hand side terms depend only on
local temperature (T), and the left hand side on $Y_{f,v}$. This equation
fully determines the thermodynamic equilibrium state of the capsule. The
equilibrium equation can be written in terms of mass fractions and molecular
weights as

$\frac{Y_{f,v}}{1-Y_{f}}=\frac{MW_{f}}{MW_{a}}\frac{1}{\frac{p_{a}}{p_{v}}-1}$
(19)

where $MW_{f}$ is the molecular weight of the fuel, and $MW_{a}$ is the
molecular weight of the ambient air. This helps to estimate the $Y_{f,v}$ term
and therefore solve the equation 18 using a root finding method to determine
the temperature.

The ideal gas assumption necessitates using corrections (Poynting) in
calculating vapor pressures [10].

$p_{v}=p^{0}_{v}*exp[V_{l}\frac{p_{a}-p^{0}_{v}}{R_{gas}T}]$ (20)

where the uncorrected vapor pressure is $p^{0}_{v}$, $R_{gas}$ is the gas
constant, $p_{a}$ is the ambient pressure. This corrected vapor pressure is
used in equation 19.

A two-way coupling mechanism based on Monte Carlo (MC) integration is used to
sample thermodynamic variables within the capsule geometry. This MC mechanism
is also used to provide gas phase source term, as a result of vaporized fuel,
to the Eulerian CFD mesh. A key feature of this implementation is the presence
of gas in both the capsule (to account for effects of vaporization and
mass/momentum) and the Eulerian mesh. To simplify our calculations, we define
the gas phase gain equal to the loss of the liquid fuel. This evaporated fuel,
now in the gas phase, is used in the source term for the mass. Further, the
energy loss of the capsules at each time step is deposited as the energy
source into the gas solution of the underlying gas cell.

### Coupling formulation

The source terms must be distributed to the underlying Eulerian mesh.
Depending on the capsule spatial and radial extent, and the underlying mesh
resolution each capsule may overlap with numerous gas phase cells. A Monte
Carlo (MC) integration method is used to distribute the source terms to each
gas phase cells. This coupling using MC is a two-way coupling, as the gas
phase properties are obtained from the capsule advection equations. A single
capsule may overlap with multiple mesh cells. This overlap complicates the
redistribution of source terms to the underlying mesh. A simple strategy is
based on the same point finding algorithm used to determine in which cell a
Lagrangian parcel might reside is used. In this algorithm for a given capsule,
for each of these Monte Carlo points, a fraction of $1/N$ of the source term
is deposited into the cell which contains the point. For example, in the
Figure 4, cell # 8 has two randomly generated MC points coinciding with it.
Therefore, $2/N$ of the total mass/momentum/energy in the capsule is added to
the gas phase source term. This stochastic approach obviates the need for
calculating expensive geometric intersections. These interior points are used
only for distributing source terms.

Figure 4: The schematic of a limited number of points distributed randomly
within the capsule geometry. Blue numbers represent cell indices. The
underlying Eulerian mesh computes the overlapping number of points within each
cell and distributes the source terms.

### Transition

As the capsule traverses along the axial direction, it ingests more ambient
gas in the form of entrainment. This causes vaporization as the ambient gas
temperatures are on average higher than the fuel temperature. This combined
effect reduces the liquid fuel in the capsule, making the capsules diluted at
a certain location downstream. At this point the Lagrangian capsule approach
is not further needed, and the capsule is ready to be transitioned to
Lagrangian parcels consisting of very small droplets. This transition has
several advantages including making use of existing frameworks to capture
parcel impacts on cylinder walls. In the current formulation, the capsules
only traverse in the direction of the initial spray, but post transition
parcels are not limited by this over-simplification. The flexibility of the
post-transition parcels will help in the model’s ability to successfully
capture the phenomena of spray bending, as is observed in the case of multi-
injector nozzles.

The criterion for transition in this current implementation is based on three
factors:

  * •

The liquid volume fraction is below 0.005

  * •

The vapor volume fraction of the total fuel volume fraction is above 99%

  * •

The capsule’s average velocity is lower than 10 $m/s$

During transition, the remaining liquid mass in the capsule is proposed to be
divided equally into a predetermined number of parcels at locations generated
randomly within the capsule’s geometrical extent. Since the capsule geometry
is small, of the order of a few Eulerian mesh cells, the random distribution
of parcels is not expected to make a great difference. The droplet size is
based on a critical Weber number of 6, as indicated by experiments [21]. Once
the transition to Lagrangian parcels is complete, the parcels evolve with
typical spray models appropriate for secondary atomization of dilute sprays.

### Overall approach

The kinematic equations are solved first using the same time-step ($\delta t$)
as the CFD solver. The discharge rate curve is interpolated based on the
overall timestep of the solution. The discharge rate determines the mass of
fuel in each capsule upon initialization. The capsules are initialized at the
nozzle exit, with an initial spray angle, and an axial length equal to the
velocity (at birth) times $\delta t$. The capsules then evolve in the domain
using equations 11. The new positions determine the volume of the capsule as
the one-dimensional implementation is projected onto three dimensions using
simplified assumptions of radial profiles. The entrained ambient air mass is
calculated and the new velocity determined by conserving the capsule momentum.
The solution of the thermodynamic state using the relationships discussed,
equations 18 and 19, requires iterative solution due to the coupling of the
variables. Once the new thermodynamic state is found, the capsule calculations
are concluded. At the end of each cycle of calculations, capsules are examined
against the transition criterion and the ones meeting it are transitioned into
smaller Lagrangian parcels.

## Results

The results presented here comprise two categories: verification and
validation. Verification provides evidence that the kinematics and
thermodynamics equations are being solved correctly, and validation provides
evidence that the model represents actual sprays. Verification of the results
begins with a spray penetration comparison to the Musculus-Kattke model [12]
in order to show that the kinematics implementation of the ELMO model is being
solved correctly. The test case is based on the Engine Combustion Network
(ECN) Spray A [18] without vaporization, and the result is shown in Fig. 5.
Mean axial velocity of fuel is plotted versus axial distance at different
times of the solution to check the spray evolution. The prediction of ELMO
model matches with the result of M-K model.

Figure 5:  Fuel jet penetration compared to Musculus-Kattke model [22].

Next, the thermodynamic implementation was tested as a verification measure to
ensure that it was faithful to the thermodynamics of Desantes et al. [23]. The
test case is Spray A, while the source terms are switched off to make the ELMO
code see a uniform ambient condition undisturbed by the spray. The results are
shown in Fig. 6. Mass fraction ratio of vapor to liquid fuel is plotted versus
fuel mass fraction. It shows the vaporization progress versus the gas
entertainment into the fuel capsule. The fuel is completely vaporized at a
fuel mass fraction slightly less than 0.4.

Figure 6: With a static, undisturbed flowfield, the current implementation
reproduces the results of Desantes et al. The abscissa represents the mass
fraction of fuel and the ordinate is the fraction of fuel in the vapor phase.

For validation, three standard spray cases from Engine Combustion Network
(ECN) were considered: Spray A, Spray H, and Spray G [24] [25] [26]. The
validation includes global spray quantities, such as vapor penetration, radial
fuel distribution, and gas motion. The ELMO model was used without any
adjustment to the model parameters in all validation cases. For all the
simulations, the based grid size is 2mm, while embedded grid refinement around
the nozzle and adaptive mesh refinement (AMR) assist the simulation.
Unresolved sub-grid scale features were modeled with RNG $k-\epsilon$
turbulence model. Taylor Analogy Breakup (TAB) [27] and No Time Counter (NTC)
[28] are considered as the breakup and collision models for the transition
parcels. Also, the Frössling [29] method is employed as the vaporization
model.

The first validation case is Spray A and the injection parameters can be found
in the reference [18]. For this case, spray cone angle is 21.5°and the
injection coefficient of area is 0.98. The initial turbulent kinetic energy
and dissipation are 5.02E-4 $m^{2}/s^{2}$ and 0.0187 $m^{2}/s^{3}$
respectively. The predictions of the ELMO model were compared to experimental
data from Pickett et al [30]. Figure 7 shows the penetration versus time (a),
the transverse distribution of fuel mass at 20 $mm$ downstream (b), and
longitudinal distribution of fuel mass from 18 $mm$ to 40 $mm$ downstream of
injection at the end of injection. For comparison, Lagrangian-Eulerian (LE)
simulation results are also included in Fig. 7. The model inputs were
specified as a part of an example case that comes with CONVERGE 2.4. The
penetration is predicted by the ELMO model without any model adjustment. Vapor
penetration result by ELMO is in good agreement with the experiment and
relatively outperformed the LE solution. The transverse and longitudinal
profiles shown in Fig. 7 (b) and (c) permit comparison of the fuel dispersion.
The fuel is entirely vaporized well before these measurement locations. ELMO,
similar to LE, predicts the bell-shaped distribution with the approximate
width and height of the experimental measurements. For the longitudinal
profile, ELMO shows some over-prediction compared with the the experimental
data at locations near the nozzle, while it switches to under-prediction as we
move downstream. There is a similar trend in the LE solution. It worth noting
that fuel capsules evaporate entirely within a few millimeters of injector
exit, and other factors such as turbulence and evaporation models could affect
the downstream results.

a)

b) c)

Figure 7: Spray A simulation compared against experimental result of Pickett
et al [30]. a) vapor penetration, b) transverse and c) longitudinal mixture
fraction distribution.

The next validation case is the ECN Spray H [31]. The injected fuel is
n-heptane, which is more volatile than dodecane in Spray A. Injection
parameters can be found in the reference [ECN]. The value of 24°is considered
as the spray angle [30], while the nozzle coefficient of area is 0.86. The
initial turbulent kinetic energy and dissipation are 5.02E-4 $m^{2}/s^{2}$ and
0.0176 $m^{2}/s^{3}$ respectively. The ELMO results along with the Lagrangian
Eulerian solution are shown in figure 8 including vapor penetration (a),
transverse (b) and longitudinal fuel distribution (c) at 1.13 ms after the
start of injection. ELMO outperformed LE in terms of vapor penetration
prediction. About fuel mixture fraction, ELMO prediction is consistent with
the Lagrangian Eulerian solution. Again, these experimental data are located
far beyond the point where ELMO capsules are terminated, and so other factors
may play a role. Because Spray H is more volatile and the gas temperature so
high, the capsules vaporize quickly and the transition to parcels occurs only
at very early times, near the start of injection. For the large majority of
the simulation, the fuel vaporizes directly and completely from the capsules.

a) b) c)

Figure 8: Spray H simulation compared against experimental result of Sandia
lab. a) vapor penetration, b) transverse, and c) longitudinal mixture fraction
distribution.

The final test for the ELMO model is the gasoline spray G condition [ECN, 32].
This eight-hole gasoline direct injection uses iso-octane under non-flashing
conditions. This test case is fundamentally different than the two previous,
not only because this is a multi-hole injector, but because of the far lower
ambient density. These gasoline direct injection conditions do not correspond
to the diesel-relevant conditions under which Siebers composed the mixing-
limited hypothesis. There is no existing theoretical boundaries that define
when the mixing-limited hypothesis is applicable and when it is not, so this
test case provides empirical illumination of the model applicability.

The two most critical spray input values are the spray cone angle = 20°and
nozzle angle = 34°, taken from the work of Payri et al. [32], and the
coefficient of area, $C_{a}=0.68$, from Moulai et al. [33]. The initial
turbulent kinetic energy and dissipation are 6.4E-3 $m^{2}/s^{2}$ and 0.486
$m^{2}/s^{3}$ respectively. The eight hole injection was reduced to a half
sector with four nozzles. For validation, we use vapor penetration data from
the Engine Combustion Network and measured gas velocity data from Sphicas et
al. [34]. The penetration data are measured in the direction of the injector
axis, not the hole axis, as a function of time. The gas velocity data come
from transient measurements made between plumes, at a point 15 $mm$ below the
injector tip. Figure 9 provides experimental validation for the predicted
penetration (a) and the axial gas velocity at a point on the axis of the
injector (b). The vapor penetration plot shows a good agreement between ELMO
and experimental data. The axial velocity validation is based on the
observation that the downward spray motion induces a temporary upward motion
of air along the injector axis. These velocity data show that shortly after
the start of injection, the ambient gas begins moving towards the injector. In
this test, the ELMO model predicts the right pattern matching the experimental
data, similar to the LE simulation.

a)  b)

Figure 9: Spray G simulation compared against experimental result of Sandia
lab. a) vapor penetration measured in the direction of the injector axis, b)
Axial gas velocity at the centerline located at a downstream position of z =
15 $\mathrm{mm}$. Positive values indicate motion away from the injector and
negative values indicate motion towards the injector.

## Conclusions

A mixing-oriented spray model has been conceived, implemented, and tested.
ELMO is based on ideas of thermodynamic and inertial equilibrium in the dense
spray core, resulting in an approach that targets mixing-limited conditions.
ELMO, as demonstrated by the results generated here, can be applied to diesel
sprays and gasoline sprays without adjustment. Compared to the drop-oriented
models, ELMO has few arbitrary physical constants, providing little temptation
to tune the model. This makes the predictions reproducible and predictive. The
computational cost of the model was roughly equivalent to traditional
Lagrangian/Eulerian simulation. In these non-combusting cases, the spray
represents roughly eight percent or less of the total computational cost, so
the cost of the ELMO model is not of great concern.

This work represents an attempt to make a fully-coupled, mixing-limited, CFD
computation. The comparison with Lagrangian-Eulerian spray modeling is a
particular challenge; given that Lagrangian-Eulerian models have been
optimized and tuned in the last two decades. The fact that the ELMO model
produces comparable accuracy in this initial effort is encouraging.

## Acknowledgements

Funding for the project was provided by the Spray Combustion Consortium of
automotive industry sponsors, including Convergent Science Inc. Cummins Inc.
Ford Motor Co. Hino Motors Ltd. Isuzu Motors Ltd. Groupe Renault, and Toyota
Motor Co. Experiments were performed at Combustion Research Facility, Sandia
National Laboratory with support from the U.S. DOE Office of Vehicle
Technologies. Sandia National Laboratories is a multi-mission laboratory
managed and operated by National Technology and Engineering Solutions for
Sandia LLC, a wholly owned subsidiary of Honeywell International, Inc. for the
U.S. Department of Energy’s National Nuclear Security Administration under
contract DE-NA0003525. The authors also thank Dr. Yajuvendra Shekhawat of
Convergent Science Inc. for his help in using the CONVERGE code’s full
capabilities.

## References

  * [1] David P Schmidt and Frederick Bedford.  International Journal of Multiphase Flow, 102:95–101, 2018.
  * [2] John K Dukowicz.  Journal of Computational Physics, 35(2):229–253, 1980.
  * [3] S. Balachandar, Kai Liu, and Mandar Lakhote.  Journal of Computational Physics, 376:160 – 185, 2019.
  * [4] Dennis L Siebers.  SAE transactions, pp. 703–728, 1999.
  * [5] Jeffrey D Naber and Dennis L Siebers.  SAE transactions, pp. 82–111, 1996.
  * [6] Dennis L Siebers.  SAE transactions, pp. 1205–1227, 1998.
  * [7] Arpit Agarwal and Mario F Trujillo.  International Journal of Multiphase Flow, 109:1–13, 2018.
  * [8] José V Pastor, J Javier López, José M García, and José M Pastor.  Fuel, 87(13-14):2871–2885, 2008.
  * [9] Jose V Pastor, Raul Payri, Jose M Garcia-Oliver, and Francisco J Briceno.  Atomization and Sprays, 21(6):503–520, 2011.
  * [10] Jose M Desantes, J Javier Lopez, Jose M Garcia, and Jose M Pastor.  Atomization and Sprays, 17(3), 2007.
  * [11] JM Desantes, JV Pastor, JM García-Oliver, and JM Pastor.  Combustion and Flame, 156(1):234–249, 2009.
  * [12] Mark PB Musculus and Kyle Kattke.  SAE International Journal of Engines, 2(1):1170–1193, 2009.
  * [13] N Abani, S Kokjohn, SW Park, M Bergin, A Munnannur, W Ning, Y Sun, and Rolf D Reitz.  An improved spray model for reducing numerical parameter dependencies in diesel engine CFD simulations.  Technical report, SAE Technical Paper, 2008.
  * [14] Zongyu Yue and Rolf D Reitz.  International Journal of Engine Research, p. 1468087417744144, 2017\.
  * [15] Zongyu Yue and Rolf D Reitz.  Energy & Fuels, 33(4):3565–3575, 2019.
  * [16] Federico Perini, Stephen Busch, and Rolf Deneys Reitz.  The Journal of Supercritical Fluids, 167:105042, 2021.
  * [17] Madhusudan G Pai and Shankar Subramaniam.  Journal of Fluid Mechanics, 628:181, 2009.
  * [18] Lyle M Pickett.  Engine combustion network, https://ecn.sandia.gov.
  * [19] PK Senecal, E Pomraning, KJ Richards, and S Som.  Journal of Energy Resources Technology, 136(1), 2014.
  * [20] Converge 2.4 manual, 2017.
  * [21] Heinz Pitsch.  Annu. Rev. Fluid Mech., 38:453–482, 2006.
  * [22] Mark PB Musculus and Kyle Kattke.  SAE International Journal of Engines, 2(1):1170–1193, 2009.
  * [23] Jose M Desantes, J Javier Lopez, Jose M Garcia, and Jose M Pastor.  Atomization and Sprays, 17(3), 2007.
  * [24] Alan L Kastengren, F Zak Tilocco, Christopher F Powell, Julien Manin, Lyle M Pickett, Raul Payri, and Tim Bazyn.  Atomization and Sprays, 22(12):1011–1052, 2012.
  * [25] Daniel J Duke, Alan L Kastengren, Katarzyna E Matusik, Andrew B Swantek, Christopher F Powell, Raul Payri, Daniel Vaquerizo, Lama Itani, Gilles Bruneaux, Ronald O Grover Jr, et al.  Experimental Thermal and Fluid Science, 88:608–621, 2017.
  * [26] Anne Kösters and Anders Karlsson.  Atomization and Sprays, 26(8), 2016.
  * [27] Peter J O’Rourke and Anthony A Amsden.  The TAB method for numerical calculation of spray droplet breakup.  Technical report, SAE Technical Paper, 1987.
  * [28] David P Schmidt and CJ Rutland.  Journal of Computational Physics, 164(1):62–80, 2000.
  * [29] N. M. Frössling.  Gerlands Beitr.Geophys., 52:170–216, 1938.
  * [30] Lyle M Pickett, Julien Manin, Caroline L Genzale, Dennis L Siebers, Mark PB Musculus, and Cherian A Idicheria.  SAE International Journal of Engines, 4(1):764–799, 2011.
  * [31] Cherian A Idicheria and Lyle M Pickett.  SAE Transactions, pp. 490–504, 2007.
  * [32] Raul Payri, Jaime Gimeno, Pedro Marti-Aldaravi, and María Martínez.  Nozzle flow simulation of GDi for measuring near-field spray angle and plume direction.  Technical report, SAE Technical Paper, 2019.
  * [33] Maryam Moulai, Ronald Grover, Scott Parrish, and David Schmidt.  Internal and near-nozzle flow in a multi-hole gasoline injector under flashing and non-flashing conditions.  Technical report, SAE Technical Paper, 2015.
  * [34] Panos Sphicas, Lyle M Pickett, Scott Skeen, Jonathan Frank, Tommaso Lucchini, David Sinoir, Gianluca D’Errico, Kaushik Saha, and Sibendu Som.  SAE International Journal of Fuels and Lubricants, 10(1):184–201, 2017.