LES of turbulent liquid jet primary breakup in turbulent coaxial air flow

A robust two-phase flow Large Eddy Simulation (LES) algorithm has been developed and applied to predict the primary breakup of an axisymmetric water jet injected into a surrounding coaxial air flow. The high liquid/gas density and viscosity ratios are known to represent a significant challenge in numerical modelling of the primary breakup process. In the current LES methodology, an extrapolated liquid velocity field was used to minimise discretisation errors, whilst maintaining sharp treatment of fluid properties across the interface. The proposed numerical approach showed excellent robustness and high accuracy in predicting coaxial liquid jet primary breakup. Since strong turbulence structures will develop inside the injector at high Reynolds numbers and affect the subsequent primary breakup, the Rescaling and Recycling Method (R2M) was implemented to facilitate generation of appropriate unsteady LES inlet conditions for both phases. The influence of inflowing liquid and gas turbulent structures on the initial interface instability was investigated. It is shown that liquid turbulent eddies play the dominant role in the initial development of liquid jet surface disturbance and distortion for the flow conditions considered. When turbulent inflows were specified by the R2M technique, the predicted core breakup lengths at different air/water velocities agreed closely with experimental data.


Introduction
Atomisation of liquid jets in coaxial air flow (air-blast or airassisted atomisation) has been widely used in combustion systems of gas turbines and rocket engines. Rapid liquid fuel atomisation exerts an important influence on fuel/air mixing, and thus affects combustion performance significantly. Study of this atomisation process is fundamentally important, but also very challenging.
In order to describe the atomisation of a single round liquid jet injected into a coaxial annular gas flow, the following characteristic non-dimensional parameters have traditionally been used: the gaseous Weber number We G , liquid and gas Reynolds numbers Re L and Re G , momentum flux ratio M, and momentum ratio MR, defined as: Here, D L is the liquid jet round nozzle diameter; D G is the hydraulic diameter of the annular gas nozzle; A L and A G are cross-sectional areas of round liquid and annular gas nozzles; q L and q G are the densities of liquid and gas; l L and l G are liquid and gas dynamic viscosities; U L is the liquid injection speed; U G is the velocity of the coaxial gas flow; finally r is the liquid surface tension coefficient. Experimental studies of air-assisted atomisation using a coaxial jet configuration have been carried out by many researchers; a recent review by Dumouchel (2008) has provided a useful summary. Faragó and Chigier (1992) classified the air-assisted atomisation into five regimes (axisymmetric Rayleigh breakup, non-axisymmetric Rayleigh breakup, membrane breakup, fibre breakup, and superpulsating breakup) via a map of gaseous Weber number vs. liquid Reynolds number. Lasheras and coworkers (Lasheras et al., 1998;Lasheras and Hopfinger, 2000) carried out their experiments using a different atomiser in term of geometrical dimensions (liquid jet diameter, gas/liquid diameter ratio), and suggested that the momentum flux ratio M is an important and additional parameter to Re L and We G for a universal classification of air-assisted atomisation.
The primary breakup of a liquid jet in a coaxial flow can be divided into two stages: initial jet surface perturbation is triggered near the nozzle exit; this perturbation is then amplified under the influence of aerodynamic forces, resulting in jet breakup. When the liquid flow is turbulent inside the nozzle (as is inevitable if the central and/or annular jet Reynolds numbers are large enough),  and Mayer and Branam (2004) argued that the initial perturbation arose from eddies originating in the liquid jet. For the case that jets are injected from a nozzle under laminar conditions, Marmottant and Villermaux (2004) suggested that the initial destabilisation is caused by a Kelvin-Helmholtz instability; the most unstable wavelength is then proportional to the thickness of the gaseous boundary layer formed in the annular nozzle. In the second stage, the initial surface perturbations grow due to aerodynamic interactions, liquid structures protruding from the liquid surface are accelerated by form drag due to the gas flow, making them subject to the Rayleigh-Taylor instability, and finally ligaments and droplets disintegrate from the liquid jet surface.
The liquid core length (or liquid jet breakup length) L C is the axial downstream location where the continuity of the liquid jet discharged from the nozzle exit is interrupted over the entire jet cross-section, and is considered a fundamental and important parameter for evaluation of atomisation performance. Measurement of L C has been carried out by many authors (see Engelbert et al., 1995;Lasheras et al., 1998;Porcheron et al., 2002;Leroux et al., 2007), and several correlations have been proposed as indicated in Table 1. Although each correlation shows appropriate agreement with the experimental data from which it was deduced, no single correlation is able to predict correct liquid core length for other experiments. Since the characteristics of the flow developed inside the nozzles can considerably influence the primary breakup, the liquid core length will inevitably be strongly dependent on the details of the injector geometry. This is a primary cause of the scatter or discrepancy in currently available empirical correlations. Another cause of inaccuracy is measurement error in the shadowgraph technique, which has commonly been used in experimental studies. Droplets stripped off the periphery of the central liquid jet during the early stages of breakup can obscure observation of the liquid core due to the line-of-sight nature of the technique. A novel optical technique, based on internal illumination of the continuous liquid jet core by Laser Induced Fluorescence (LIF) has been proposed by Charalampous et al. (2009a,b) for conducting measurements of liquid core length. Their data demonstrated that LIF can provide more accurate detection of the liquid jet geometry than the shadowgraph technique.
Numerical modelling of liquid jet atomisation has made significant progress since the 1970s, with the promise of eventually achieving as much success as for single phase flow (Fuster et al., 2009;Gorokhovski and Herrmann, 2008). In the last decade or so, most proposals have adopted the more expensive but more advanced Direct Numerical Simulation (DNS) or Large Eddy Simulation (LES) approach to turbulence modelling in order to capture unsteady effects on the interface dynamics. Unsteady numerical simulations can show many more details than are possible to capture in experiments, providing further insight into atomisation mechanisms. However, such numerical modelling of primary breakup of a liquid jet under the influence of strong aerodynamic and turbulence effects is still very challenging, especially for high liquid/gas density ratio O(1000). Numerical error arising from the high density ratio can be large (sufficient even to cause the simulation to fail), and many published numerical simulations to date have been limited to a liquid/gas density ratio no more than 100 (see Herrmann, 2010, Level Set (LS) interface tracking method; Herrmann et al., 2011, LS;Desjardins et al., 2008, LS;Pai et al., 2008, LS;Kim et al. (2007), LS ;Fuster et al., 2009, Volume of Fluid (VOF); Tomar et al., 2010, VOF;Ménard et al., 2007, Coupled LS and VOF (CLSVOF); Lebas et al., 2009, CLSVOF;Shinjo and Umemura, 2010, CLSVOF). Since the majority of liquid jet atomisation experiments are carried out at atmospheric pressure with high density liquids, quantitative comparison between numerical modelling and experiment is quite rare. A robust method capable of dealing with the high density ratio of for example air-water systems is therefore of great interest (Fuster et al., 2009). In order to deal with this problem, Rudman (1998) proposed to advect the momentum using a density estimated from the interface geometry in cells intersected by the surface, aiming to improve the consistency between interface and momentum transport. This technique as well as a two-velocity Ghost Fluid Method were investigated by Desjardins and Moureau (2010), although they have not yet been demonstrated to work well in simulating liquid jet atomisation. Sussman et al. (2007) proposed an approach using a liquid velocity field extrapolated across the interface into the gas phase region. This approach was applied by Li et al. (2010) to simulate a liquid jet in air cross-flow at a high (650) density ratio. Adaptive Mesh Refinement and the removal of under-resolved small liquid structures were both necessary since the experimental data selected for comparison were far downstream. In spite of the advanced modelling, agreement with measurements was relatively poor (Li et al., 2010). A robust two-phase LES algorithm also making use of a modified extrapolated liquid velocity field has recently been proposed by Xiao (2012), and validated against experiments by simulating droplet and liquid jet primary breakup; it was demonstrated that the proposed method showed high robustness and good accuracy compared with experiments for air-water systems.
The objectives of the current paper are therefore: (i) to simulate a round water jet injected into a coaxial air flow at atmospheric pressure (a high density ratio of 830) and compare the predicted results directly with the experimental data and (ii) to investigate the mechanisms behind the initial jet surface disturbance and the liquid jet primary breakup. Note: the aspects of liquid jet primary breakup which are given particular focus in the present paper are the initial destabilisation of the liquid/air interface and the location of first complete rupture of the jet core. Whilst the subsequent ligament and droplet formation is captured in the simulations shown, the measured data used do not allow quantitative assessment of these aspects, and, in the far downstream region of the solution domain, the mesh density currently used is inadequate for this purpose, so this has been left for a separate study. The experimental tests of Charalampous et al. (2009a,b) are simulated, as the LIF technique used there can provide more accurate measurements of liquid jet core length. Due to the high Reynolds number of both liquid and gas streams, turbulent boundary layers undoubtedly develop on the internal nozzle walls of the injector system, and may influence significantly the primary breakup process. As in any LES prediction, the generation of unsteady 3D correlated inlet conditions is a challenging task. The Rescaling and Recycling Method (R 2 M) developed by Xiao (2012) and Xiao et al. (2010) has shown good performance in single phase flows, and will therefore be implemented here. Since the inflow conditions can be specified as whatever one wishes in numerical modelling, the effect of turbulence on primary breakup is investigated by switching between laminar and turbulent inflow conditions. Since the breakup length is a parameter of great importance in performance assessment of air-assisted atomisation, the dependence Table 1 Correlations for liquid core length.  Lc relationship of the breakup length on the characteristic parameters (M, We G , Re L ) will be carefully examined.
In this paper, we first give a description of the two-phase LES formulation and associated numerical methods; the R 2 M approach for LES inlet condition is then described. Finally, a detailed analysis of the simulation results and both a qualitative and quantitative comparison between numerical predictions and experiments will be given.

Governing equations
In the current study of two-phase flow modelling, both liquid and gas are assumed to be incompressible and immiscible. The derivation of the two-phase flow LES formulation is described in detail in Xiao (2012). The philosophy behind this formulation is as follows: the liquid/gas interface is resolved directly without modelling Sub-Grid Scale (SGS) features; the usual spatially filtered LES formulation is employed in the single-phase flow regions; an appropriate treatment is adopted when discretising the governing equations (when interpolating cell face fluxes) for cells which are intersected by the interface for both phases. Although a sub-grid interface dynamics model has been proposed by Herrmann (2013) for sub-filter surface tension induced interface motion, it requires an interface tracking or capturing method that can provide significant sub-filter resolution, making it difficult to implement in the current code. Similarly, a modelling approach to sub-grid surface tension was recently proposed by Aniszewski et al. (2012) and has produced interesting results by comparing LES results with a-priori DNS results, but this still needs further development and validation before application to the high Re flow which is of prime interest in the current paper. Therefore, these sub-grid interface models have not been implemented in the current LES formulation. Effectively, the present LES formulation is similar to that referred to as quasi-DNS/LES in Gorokhovski and Herrmann (2008), i.e. an under-resolved DNS of interface tracking combined with an LES of the single-phase regions of the flow.
A coupled Level Set (LS) and VOF (Volume of Fluid) method (CLSVOF) is used here to capture and evolve the interface. The Level Set / is a signed distance function from the interface satisfying r/ = 0. The interface is defined by / = 0, with / > 0 representing liquid and / 6 0 representing air. / is evolved by the simple advection equation using the resolved velocity field (note: in what follows a spatially filtered (resolved) quantity is indicated by an overbar) and ignoring the contribution of any Sub-Grid-Scale (SGS) velocity effects: To maintain the signed-distance property, the re-initialisation equation is also solved: where s represents pseudo time, S(u 0 ) is a modified sign function, u 0 ¼ uðx i ; s ¼ 0Þ ¼ /ðx i ; tÞ, d = max(Dx, Dy, Dz). After solving this equation to the steady state in the interface vicinity, / is replaced by u.
The VOF function F is defined as the volume fraction occupied by the liquid in each computational cell. The resolved evolution of the VOF function is governed by: With residual or SGS stress tensor s r ij modelled by a simple Smagorinsky eddy viscosity approach, the governing equations for the resolved velocity field are: Here, g i is gravitational acceleration, D represents the filter width, taken as the cube root of the local cell volume and the value of the Smagorinsky constant C S is set in all the calculations reported below to 0.1. The surface tension term F ST i is computed via: Here, r is the surface tension coefficient, j is the interface curvature, and n i is interface normal vector pointing from the liquid phase to the gas phase. The SGS term arising from filtering the surface tension term is neglected. To maintain the implication of a sharp interface, fluid density and viscosity are in the present approach not considered as spatially filtered quantities, but are set to be the properties of liquid or gas depending on the local value of the resolved Level Set variable /.

Coupled Level Set and VOF method (CLSVOF)
A CLSVOF methodology (Sussman and Puckett, 2000) for interface advection is adopted in the present approach since it offers an optimum combination of the good mass conservation property of the VOF approach and the convenient and accurate capability of LS for evaluation of interface geometrical properties. Coupling of LS and VOF variables is enforced in the interface reconstruction step where both LS and VOF information can be used to good effect: the interface normal vector is computed from LS (a smooth function) rather than VOF (a discontinuous function), since LS gives more accurate information on interface location and shape, the interface position in the cell is constrained by the VOF function, since this gives more accurate information on liquid volume conservation.
Based on such a reconstructed interface, the VOF function is evolved to the next time step, and the LS function is corrected. A flow chart of the CLSVOF method is shown in Fig. 1. The coupling of the LS and VOF methods occurs during the interface reconstruction and LS-redistance processes. The detailed algorithm of the present CLSVOF method is as follows: Initialise the LS and VOF functions at time step n = 0: / n and F n . Reconstruct the interface in cells where 0 < F n < 1. The interface normal vector n i is calculated from the LS function, and the position of the interface within the cell is constrained by the VOF function. Advect the VOF function from F n to F nþ1 based on the reconstructed interface. Advect the LS function from / n to / nþ1;Ã . The operator split method is used for solving both equations (see Xiao, 2012 andSussman andPuckett, 2000 for details). Reconstruct and constrain the interface in cut cells using the new LS function / nþ1;Ã and the VOF function F nþ1 . Perform a re-initialisation step on / nþ1;Ã to obtain the final level set function / nþ1 with a recovered signed distance property.
In the current CLSVOF method, the normal vector is calculated directly by discretising the LS gradient using a finite difference scheme. By appropriately choosing one of three finite difference schemes (central, forward, or backward differencing), it has been demonstrated that thin liquid ligaments can be well resolved see Xiao (2012). Although a high order discretisation scheme (e.g. 5th order WENO) has been found necessary for LS evolution in pure LS methods to reduce mass error, low order LS discretisation schemes (2nd order is used here) can produce accurate results when the LS equation is solved and constrained as indicated above in a CLSVOF method (see Xiao, 2012), since the VOF method maintains 2nd order accuracy. This is a further reason to adopt the CLSVOF method, which has been used for all the following simulations of liquid jet primary breakup.

Discretisation of two-phase flow transport equations
Since the convection and diffusion terms are discontinuous across the interface, a cautious first order forward-Euler projection method was used for temporal discretisation of the two-phase flow governing equations (for more details see Xiao (2012)). First, an intermediate velocity is computed from convection, diffusion and gravitational terms (spatial discretisation is described below) (NB surface tension is treated via the pressure term using the ghostfluid approach (Fedkiw et al., 1999;Kang et al., 2000;Liu et al., 2000)): Second, the intermediate velocity field is updated using a pressure gradient term to obtain the velocity at time step n + 1: Since the velocity field at time step n + 1 must satisfy the continuity equation, a pressure Poisson equation may be derived by taking the divergence of the above equations to allow P nþ1 to be calculated: The variables are arranged in a staggered manner in the current two-phase flow LES formulation: the pressure, LS and VOF are located at the cell centre; velocity components are located at corresponding faces. In general, 2nd order central methods are used to discretise spatial derivatives. Since the gas phase has a much smaller density and viscosity than the liquid phase, the velocity gradient in the gas phase is typically much larger than liquid phase. However, for discretisation of the momentum equation for cells in the vicinity of the interface, it was found important to specify correctly which velocity should be used when calculating cell face fluxes associated with convection and diffusion terms. If the velocity in the adjoining gas phase cell was used to discretise the momentum equation for an interface cut liquid phase cell, large momentum errors were induced in the vicinity of the interface. In order to reduce this momentum error (arising from the density jump across the interface) an extrapolated liquid velocity field U L i was introduced for spatial discretisation of convection and diffusion terms for interface cut cells. The technique of an extended velocity field was used by Nguyen et al. (2001) to simulate incompressible flames. A projection method was used by Tanguy et al. (2007) to impose a divergence-free condition for the extrapolated velocity field across the interface to tackle the problem of parasitic currents in their simulations of two-phase vaporising flows. Though an extrapolation projection approach was used by Sussman (2003) to impose divergence-free condition in simulations of bubble growth and collapse, it was abandoned by Sussman et al. (2007) in their later numerical formulation of incompressible two-phase flows. In the current formulation, a constant extrapolation technique (Fedkiw et al., 1999;Aslam, 2003) was used for liquid velocity extrapolation, and an original method was developed to imposed the divergence free condition for the extrapolated liquid velocity. For a detailed implementation of this technique, readers are referred to Xiao (2012). It was also demonstrated in Xiao (2012) that: (i) a divergence free step for the extrapolated liquid velocity field was necessary to reduce the momentum error when the interface moved rapidly across a fixed grid and (ii) a more accurate capture of interface dynamics was obtained when the extrapolated liquid velocity field was also used in the VOF and LS transport equations. The methodology developed in Xiao (2012) was used for all simulations reported here, including the use of a Box multigrid (Dendy, 1982) preconditioned conjugate gradient method for solution of the Poisson equation. Fig. 2 summaries the final procedure of the developed twophase flow LES. The detailed algorithm for the current two-phase flow LES is as follows:

Algorithm for two-phase flow LES
Based on the interface represented by the LS function / n , discretise the two-phase flow governing equations to solve for the velocity field at the next time step U nþ1 i . Construct the extrapolated liquid velocity field U L;nþ1 i using the extrapolation technique described in Xiao (2012). Ensure continuity for the extrapolated liquid velocity in the gas phase by a divergence free step. Based on U L;nþ1 i advect the LS and VOF functions to the next time step to obtain / nþ1 and F nþ1 using the CLSVOF algorithm in Section 2.2.1 and Fig. 1. Set the velocity in CVs which change from gas to liquid (i.e.

Repeat for further time steps.
The numerical implementation of the methodology described above was carried out in an existing multi-block structured mesh code for Large Eddy Simulation (LES) of single phase constant density turbulent flow. A full description of the numerical techniques adopted in this code (LULES) may be found in Tang et al. (2004). The code solves for contravariant velocity components on a 3D orthogonal mesh, to give the code some complex geometry capability. It has been validated against a range of incompressible turbulent flow problems, for example the unsteady aerodynamics of high swirl fuel injectors (Cheng et al. (2012). For the test cases reported here only Cartesian or cylindrical polar meshes are needed.

Rescaling and Recycling Method (R 2 M)
The Rescaling and Recycling Method (R 2 M) proposed in Xiao et al. (2010) and Xiao (2012) has demonstrated a validated capability for generation of unsteady 3D correlated velocity fields for specification of LES inflow conditions. The R 2 M technique was therefore used here to generate nozzle exit conditions matched as far as possible to the experimental configuration selected for LES of liquid jet primary breakup. In order to generate realistic unsteady inflows for LES using the R 2 M technique, target values for statistical quantities (mean velocity and Reynolds normal stress (or rms intensity) fields) at the ''main simulation'' (MS) solution domain inlet plane need to be prescribed (i.e. for axisymmetric flow inside both liquid and gas nozzles: U target (r), V target (r), W target (r), u 0 target ðrÞ, v 0 target ðrÞ, w 0 target ðrÞ, where r is the radial co-ordinate and U, u 0 (etc) represent mean and rms velocities respectively. How these values are specified and the R 2 M technique applied is described below.
First, extra ''inlet condition'' (IC) solution domains (single block or multi-blocks as necessary) are created upstream of the inlet plane of the MS domain (here the nozzle exit plane). The inflow conditions for the extra IC domain are generated by recycling the velocity field from a selected plane in the downstream region of the IC domain. By rescaling the velocities, and solving the LES equations in the IC domain until a statistically stationary state is achieved, the resulting instantaneous flow field within the IC domain can achieve the target statistical characteristics whilst also possessing self-consistent spatial and temporal correlations. This unsteady velocity field is then fed into the MS domain as inlet conditions.
In the simulation of round liquid jet primary breakup in coaxial gas flow, two cylindrical IC domains for liquid and gas respectively need to be created; in these domains cylindrical polar co-ordinates (x, r, h) and velocity decomposition were used for convenience. The generation of unsteady conditions for the liquid flow in the central round nozzle as well as the gas flow in the surrounding annular nozzle using the R 2 M approach was carried out as follows: 1. Cylindrical IC domains within the central nozzle (and annular gas nozzle) were specified upstream of the MS domain inlet plane (nozzle exit). The size of the IC domains in the radial and azimuthal (r, h) directions were set by the injector nozzle geometry. The size in the streamwise (x) direction was chosen so that the two point axial spatial correlations would fall to zero well within the IC domains, as required by the recycling technique and also to prevent the generated axial turbulent integral scale being constrained by the solution domain size. 2. Use recycling to provide inflow conditions for the IC domains.
The velocity field at a plane a short distance upstream of the IC domain outlets was recycled. In addition, it is important that the mesh in the IC domains should be uniform in the streamwise and azimuthal directions to avoid any varying spatial filtering effects. 3. Initialise the velocity field in the IC domains; the instantaneous velocity field was generated by superimposing white noise with an intensity of u 0 target ðrÞ, v 0 target ðrÞ, w 0 target ðrÞ onto the mean velocity U target (r), V target (r), W target (r). 4. Run the simulation in both IC and MS domains simultaneously.
Rescale the flow field everywhere within the IC domains every k time steps (k = 5 used in the current simulations) in the following way (described here for the axial component only for simplicity): Calculate the mean velocity by spatial averaging in the x and h directions (homogeneous directions) and temporal averaging with a weight that decreases exponentially backward in time (see Lund et al. (1998) Rescale the other two components V and W using the same procedure.

Test case details
In the section, the developed two-phase flow CLSVOF LES formulation is applied to predict the primary breakup of a single cylindrical water jet in a coaxial annular air flow. As stated above, the test cases from Charalampous et al. (2009a,b) are simulated. The coaxial air-blast atomiser used is shown in Fig. 3. The diameter of the liquid nozzle is D L = 2.3 mm; the inner and outer diameters of the annular gas nozzle are 2.95 mm and 14.95 mm respectively. Table 2 lists the flow conditions simulated here (U G and U L are the area-averaged velocities at the exit of the nozzle and M is the gas to liquid momentum flux ratio). The experiment was carried out with water and air at atmospheric pressure, so fluid properties used in the simulations are: q G = 1.205 kg/m 3 , l G = 1.836 Â 10 À5 Pa s, q L = 1000.0 kg/m 3 , l L = 0.848 Â 10 À3 Pa s, and surface tension coefficient r = 0.072 N/m. Note that both Re L and Re G are high enough that turbulent flow exiting the nozzle may be expected.
An initial simulation with uniform and laminar inflows for both water and air was first carried out, with the results shown in Section 4.2. The R 2 M technique was then applied to generate realistic turbulent inflows and this simulation is presented in Section 4.3. The mechanism of initial interface disturbance and subsequent liquid jet primary breakup is investigated in Section 4.4. Liquid jet structures for a range of flow conditions were then simulated and are compared with experimental images in Section 4.5; a comparison of the LES-predicted liquid core length extracted from these simulations against the experimental data is given in Section 4.6. Finally, the property of liquid volume conservation in the current LES predictions of liquid jet primary breakup is examined in Section 4.7.

Simulation with laminar inflows
The simulation domain co-ordinates were [0, 54], [À6, 6], [À6, 6] mm in the x, y and z directions respectively. The axis of the nozzle lies in the x direction, with the Cartesian co-ordinate origin at the central nozzle centre. Fig. 4 shows the Cartesian mesh used in the simulation. In order to provide good resolution of the initial stages of primary breakup, a uniform fine mesh was used in the region [0, 54], [À2.43, 2.43], [À2.43, 2.43] mm with a cell size of 0.09 mm (i.e. $D L /25); in the outer region of the domain, an expanding mesh was used to reduce computational cost. The   number of cells was 600 Â 80 Â 80 in x, y, z directions (4 million in total). The mesh in a part of the inner region is shown in Fig. 4; the red region corresponds to the initial condition for the liquid phase, i.e. the cells where VOF was initialised to 1 and LS was set to be a positive distance to the interface. Uniform and laminar inflows were first specified for both water and air streams at the nozzle exit. Fig. 5 shows an instantaneous liquid jet structure predicted by LES for U G = 47 m/s and U L = 4 m/s; this shows several significant differences in comparison with the experimental shadowgraph image. The predicted liquid jet is continuous without any breakup while ligaments and droplets are observed in the shadowgraph, and there is evidence that core continuity is disrupted towards the end of the experimental image. Some small disturbances are also observed on the interface near the nozzle exit in the experiment. However, in the simulation the predicted interface is very smooth in the first three D L downstream of nozzle exit; only further downstream do surface disturbances appear and grow, but slowly. Since the internal flows inside the nozzles will be turbulent due to the high Reynolds numbers, turbulent eddies emerging from nozzle exit may destabilise the interface. This simulation underlines the importance of accurate representation of nozzle exit conditions for correct simulation of turbulent liquid jet primary breakup.

Generation of turbulent inflows
As detailed above, in order to generate realistic LES turbulent inflow conditions using the R 2 M technique, mean velocity and rms profiles that are representative of experimental conditions at the nozzle exit are required as input. 2D axisymmetric RANS predictions using a low Re Reynolds Stress Transport turbulence model were therefore run using the Fluent CFD code for the internal nozzle flows (central (water) nozzle and outer (air) nozzle separately). Fig. 6 shows the simulation domain for the water flow inside the central nozzle of the atomiser. Note the presence of the optical fibre used for the laser illumination was included (white region in bottom left), since this represents a blockage to the internal water flow. Predicted contours of turbulent kinetic energy (TKE), are shown to provide an indication of the rapid boundary layer growth on the nozzle wall at this Reynolds number. A quad mesh with 180,000 elements was used, with the mesh generated to achieve values of y + $ 0.1 mm near the wall. Fig. 7 shows the predicted radial profiles of mean streamwise velocity U, and the three normal stress intensities u 0 , v 0 , and w 0 at the nozzle exit. These formed the required R 2 M data input for the liquid flow. Fig. 8 shows the TKE contours obtained for the gas flow inside the annular nozzle. Turbulent boundary layers were also observed to develop on the walls after the contraction section. A quad mesh with 170,000 elements was used. Fig. 9 shows the predicted profiles of mean streamwise velocity and rms levels for air flow at the nozzle exit. Again these formed the required target data for the R 2 M technique.
Figs. 10 and 11 show the domain for two-phase flow LES when realistic turbulent inflows were generated using R 2 M at nozzle exit. The main simulation (MS) domain for resolving the liquid jet primary breakup was the same as that used in Section 4.    boundary layer regions. Both IC domains had 58 Â 52 Â 68 cells. In the LES calculations of the flow in the IC domains, the instantaneous velocities from the plane x = À0.15 mm were recycled for use as inflow velocities to the IC domain inlet as explained above; a convective outflow condition was used at the outlet.
Figs. 12 and 13 show that the mean and rms values of the velocity field predicted in the IC domains agree well with the input target profiles. The profiles at locations x = À1, À4, À7 mm collapse together, indicating that the predicted turbulent flow is homogeneous in the streamwise direction.
Simulations in all domains were run simultaneously with the same time step. At every time step, the instantaneous velocities from a selected plane (x = À0.0048 m was used) were mapped from the cylindrical polar IC domain mesh onto the Cartesian MS domain mesh at the inlet plane of the MS domain. Figs. 14 and 15 show contours of instantaneous water and air streamwise velocity at plane x = À0.0048 m and also at the MS domain inlet. It may be observed that the velocity fields at the two planes show the same large eddy structures, indicating a correct mapping process. There is no doubt that the smallest resolved eddy structures are somewhat smoothed in this mapping due to a finer cylindrical mesh being used in the IC domain compared to the Cartesian mesh in the MS domain, but it is not believed this introduced significant error, since these eddies were of course much less energetic than the large scale structures. Fig. 16 presents contours of the instantaneous streamwise velocity in the z = 0 plane for both phases. It is observed (especially evident for the liquid phase) that the turbulent eddies developing downstream of the MS domain inlet have similar structures as those in the IC domain, indicating that the turbulent eddies developed by the R 2 M technique in the nozzles were convected downstream as the liquid and gas were injected from their respective nozzle exits. Fig. 17 compares the LES predicted interface topologies near the nozzle exit with the two different inflow treatments. In contrast to the predicted smooth interface region observed when applying uniform and laminar inlet conditions, the interface is disturbed by turbulent eddies right after the jets exit the nozzles when using turbulent inflow conditions. Fig. 18 shows the overall liquid jet structure predicted by LES when turbulent inflows are specified for both phases using R 2 M. The growth of the surface disturbance in the liquid jet under aerodynamic forces is well reproduced by LES, and the breakup point of the liquid core now agrees well with the experiment. Ligaments and droplets are ejected from the resulting liquid clusters, which is consistent with shadowgraph observations. It is an interesting question as to whether it is liquid or gas flow turbulence (or both) which exerts the larger influence upon the interface, and this will be investigated in the following section.

Mechanism of initial interface disturbance and liquid jet primary breakup
In order to investigate in more detail the mechanism behind the initial interface instability and disturbance, two further simulations were performed: one with turbulent air inflow but uniform and laminar water inflow, and the other with these conditions reversed. Views of the predicted liquid jet structures are shown   in Fig. 19 together with results from the simulations described above using either both laminar or both turbulent inflows. It is evident that the two-phase interface is disturbed immediately after the nozzle exit in the two simulations with turbulent water inflow, while a marked region of smooth interface exists downstream of the nozzle exit in the two simulations with laminar water inflow. The liquid eddies rather than the gas eddies are clearly responsible for the initial interface disturbance, which is indeed rational since the liquid has a much larger inertia than the gas.
When uniform laminar inflow was specified for the liquid phase, the initial interface disturbance is mainly due to the shear force exerted by the gas flow. It is thus perhaps surprising that the interface predicted with turbulent gas inflow had a longer smooth region than with laminar gas inflow. The explanation is provided by Fig. 20. With a turbulent boundary layer, the simulated gas flow had a larger low velocity recirculation (backflow) region behind the nozzle lip separating liquid and gas flows, and the gas velocity in the cell adjacent to the interface reached 10 m/s only by x = 4.4 mm. In contrast, when laminar inflow was used for the gas phase, the gas flow in the cell adjacent to the interface increased to 10 m/s only 2 mm after nozzle exit, and consequently destabilised the interface a shorter distance downstream of nozzle exit.
After the initial interface disturbance, the liquid jet surface perturbations grow under the aerodynamic interactions between gas and liquid flows, resulting in primary breakup. Fig. 21 shows the velocity vector field and pressure contours in a slice through the jet centreline. It is observed that gaseous vortices are predicted by the current LES in the wave troughs and are well resolved by the mesh. For the protruding interface waves, the pressure on the upstream side is higher than on the downstream side, and thus the form drag arising from this pressure imbalance accelerates the protruding liquid. Since the heavier liquid is accelerated by the lighter gas flow, Rayleigh-Taylor instability can develop on the protruding liquid waves, and thus finally ligaments and droplets are created and separate from the liquid jet surface as shown in Fig. 18(b); the liquid jet core then disintegrates into liquid clusters (see Figs. 18 and 19) which undergo further evolution under the action of aerodynamic forces leading to eventual complete primary breakup. As the air velocity increases, the predicted location where drops and ligaments are first seen moves towards the nozzle, in good agreement with experimental observations. Due to the increasing aerodynamic forces, the dimensions of the predicted liquid ligaments and droplets resulting from the primary breakup becomes smaller. If these structures enter a region where mesh resolution is inadequate, then the disintegration of these liquid structures into smaller droplets (secondary breakup) in the furthest downstream region is probably not well resolved for the high speed air velocity cases (U G = 119 m/s and 166 m/s).      current LES but obviously only when the R 2 M approach is applied to generate the turbulent inflow boundary conditions. As the liquid injection velocity grows from 2 m/s to 8 m/s, the liquid jet breakup position moves downstream, in good agreement with the experimental images.

Liquid core length
In the simulations, the liquid core length was determined by examining the cross-sectional area ratio of the liquid jet which was illuminated by the laser beam, in a manner very similar to the experimental LIF technique. Fig. 24 demonstrates the liquid area (A) which is illuminated directly by the laser beam originating    at the nozzle exit at plane x = 19 mm assuming total internal reflection of the light from the interface surface as done in the experiments. A 2D sketch is provided but the illuminated area in 3D is calculated in the same way. The ratio of area illuminated by laser beam to liquid nozzle exit area is calculated by a = A/A 0 . The predicted variation of a in the x direction at t = 24 ms is shown in Fig. 25. A small but non-zero value of a is chosen as the criterion to determine the liquid core length. Fig. 26 shows the instantaneous and mean liquid core length calculated basing on the criterion a = 0.05 for the case U G = 47 m/s, U L = 4 m/s when turbulent inflows were used for both phases. A difference of 5% in the simulated liquid core length is observed when changing a from 0.05 to 0.1. Fig. 27 shows the liquid core length predicted by the current two-phase flow LES for flows 1a-d in Table 2; U L is 4 m/s in all cases, with the air velocity varied to produce different Weber numbers. When laminar inflows were used for both phases, the predicted liquid core length was (as expected from the discussion above) much larger than the experimental measurements of Charalampous et al. (2009b) for lower Weber numbers. When turbulent inflows were specified for both phases, the simulated core length agreed very well with the experimental value for all Weber numbers, confirming that the initial interface perturbations caused by liquid eddies plays an important role in the resulting surface instability development and primary breakup process. For the cases with higher gaseous co-flow (U G = 119 m/s and 166 m/s), strong aerodynamic forces dominate the primary breakup process, and the liquid jet core is destroyed in a short distance and the difference of predicted core length between simulations with different inflow conditions is much smaller. Overall, liquid core lengths are well reproduced by the developed two-phase flow    LES method for all four flows when realistic turbulent inflows are provided for both phases using the R 2 M technique. Fig. 28 shows predicted liquid core lengths for different liquid velocities with the same air co-flow (U G = 70 m/s). For the two lower liquid injection velocities, the core length predicted by the current LES agrees well with the experimental measurements, but the predicted value is considerably larger than the LIF measurement for the case with highest liquid velocity U L = 8 m/s. This marked difference is in all probability due to measurement error. Charalampous et al. (2009a) observed a significant decrease of the LIF signal intensity along the jet length, due to scattering of laser light by refraction at the water-air interface which was much more wavy at the highest liquid velocity. This was observed to lead to undervalued measurements of liquid jet core length. It is noted in Fig. 23(c) that few droplets are stripped off the liquid jet in the test section for this case (U L = 8 m/s). In this instance, the shadowgraph technique is able to provide an accurate measurement, and indicated a core length of $13D L , in good agreement with the current LES. Fig. 29 shows that the predicted liquid core length was a power law function of the momentum flux ratio when turbulent inflows were specified. As the gas velocity increases while the liquid velocity was kept at 4 m/s, the predicted liquid core length decreased with a power law exponent of À0.39. As the liquid velocity increases while the gas velocity was kept at 70 m/s, the predicted core length decreased with a power law exponent of À0.35. The simulations carried out for uniform and laminar inflows produced a predicted power law exponent of À0.5. These predictions are consistent with the power law exponents (from À0.5 to À0.3) reported based on a range of experimental measurements (see Lasheras et al., 1998;Engelbert et al., 1995;Leroux et al., 2007), and confirms the possibility that the difference of exponent values reported by different authors results from different atomisers used. Finally, Fig. 30 indicates that the predicted liquid core length (with turbulent inflows) decreases as a power law function of Weber number when the liquid velocity is kept at 4 m/s, with a power law exponent of À0.39, which is in good agreement with the value of À0.4 reported by .

Liquid volume error in predicted primary breakup
It is very important to conserve fuel mass in simulations of fuel atomisation in an engine, and hence the liquid volume error in any    (Charalampous et al., 2009b) for different U L with U G = 70 m/s. method that claims to be a good candidate for primary breakup simulation is of prime importance. This was one of the main reasons for adopting the CLSVOF methodology. This particular aspect was therefore examined in the present simulations for a test case with U L = 4 m/s and U G = 166 m/s. When it was judged that the simulated liquid jet primary breakup had reached a statistically stable state (e.g. as shown in Fig. 26), the following three quantities were calculated from the LES results (the solution domain considered was the region enclosed by the black line in Fig. 31, an axial distance of 9.4D L , corresponding to 230% of the liquid core length for this case): the liquid volume injected into the domain at inlet during a time period t: V in (t), the liquid volume flowing out of the outlet V out (t), and the liquid volume change in the domain over this time period: DV(t) = V(t) À V(0) (V(t) is the liquid volume in the domain at time (t)). The first two were estimated from surface integrals on the domain boundaries and the third by a volume integral of the VOF function within the domain. If liquid mass is conserved, the three quantities should satisfy the relation: V in (t) = V out (t) + DV(t). Fig. 32 presents the liquid volume budget for the simulation domain shown in Fig. 31. Since the liquid mass influx is constant, V in (t) grows linearly with time. However, the rate of liquid volume exiting the simulation domain varies due to turbulent fluctuations and flapping of the liquid jet, and this also results in a temporal variation of DV(t). It is observed in Fig. 32 that the calculated value of V out (t) + DV(t) collapsed onto with an error of only 0.01%, demonstrating that the liquid volume is conserved well in the simulation of liquid jet primary breakup. In order to examine whether the divergence free step applied to the extrapolated liquid velocity field was influential in the conservation performance of the scheme, a simulation without the divergence free step was run. The result shown in Fig. 33 indicated that a significant increase in liquid volume error was observed in the primary breakup process (to 3.6%), confirming the importance of this part of the numerical algorithm.

Grid resolution
The mesh used in the above simulations is referred to now as Mesh I. A finer mesh, named Mesh II, was generated with a cell size of 0.06 mm in the central uniform-mesh region of the MS domain (a 50% reduction in grid spacing); the mesh in the IC domain was kept the same as that in Mesh I. Mesh II was used to investigate the influence of mesh size, and the results obtained are described below.
Due to the high liquid/gas density ratio of the present flow problem ($800) and the high Re of the flow, the initial interface perturbation of a turbulent liquid jet in a turbulent coaxial gas flow and its early stage development is caused not by any fundamental Plateau-Rayleigh or Rayleigh-Taylor instability, and also not by gaseous turbulent eddies in the gas-flow boundary layer, but by the liquid turbulent eddies convected out of the nozzle, and this is clearly shown in Fig. 19. The large liquid eddies have large turbulent energy and cause large-scale interface disturbance and the restoring surface tension corresponding to large-scale interface disturbances is weak due to small curvature. Small-scale liquid eddies (including those not resolved in the LES) have small turbulent energy and cause small-scale interface disturbance, the restoring surface tension corresponding to small-scale interface disturbances is strong due to large curvature. Therefore, the small-scale disturbances due to small eddies will be smoothed quickly by the strong surface tension near the nozzle exit. The development of the large-scale disturbance due to large liquid eddies (well resolved by the current LES mesh and turbulence inlet condition treatment) thus dominates the interface disturbance which is principally responsible for the subsequent bulk liquid core breakup location. Fig. 34 shows the liquid column surface disturbance due to liquid eddies on Mesh I and II for the case U L = 4 m/s where the gas velocity has also been reduced to 4 m/s to minimise the effects of aerodynamic forces. It is seen that the resulting large-scale interface disturbances are well resolved on both Mesh I and Mesh II. The small-scale disturbances which are under-resolved on Mesh  I but appear on Mesh II have only a marginal effect on the development of the surface disturbance. After the initial interface distortion develops due to the liquid eddies, high and low pressure regions form as the gas flows around the disturbed but still attached interface structures. The form drag due to upstream high pressure on the protruding structures and low pressure on their lee side accelerates the protruding liquid, amplifying the structure growth, which eventually results in breakup of the liquid jet column as well as formation of large liquid clusters (note that the shear drag from the gas boundary layer contributes little at high Re O(10 3 -10 4 )). Therefore the important processes leading to breakup of the bulk liquid column are well resolved on the current mesh, assuring the validity of this numerical scheme for prediction of liquid core length. It is also observed that the pressure on the upstream side of the protruding interface wave is higher than that on the wave crest, causing further protrusion of the liquid structure. For the case with gas velocity of 47 m/s the formation of protruding liquid structures and ligaments is adequately resolved as far as the column breakup process is concerned on the current mesh as the comparison with experiment shown in Figs. 18 and 22 demonstrates. For cases with higher velocities, the breakup of the created liquid ligaments may be under-resolved, but this should not impact negatively on the liquid column breakup process.
Whilst under-resolved or inadequate SGS modelling might influence small scale droplet characteristics, it should not affect the liquid jet core breakup length L C . This was confirmed by a grid resolution study carried out for the case with the highest gas velocity (166 m/s) and the lowest liquid velocity (4 m/s) where the resolution problem is most demanding. Fig. 35 shows the predicted liquid jet breakup for this case on Mesh II, which may be directly compared with the Mesh I results seen in Fig. 22c. While more small drops and ligaments are resolved on Mesh II than on Mesh I in the far downstream region, the liquid core breakup process and the resulting large liquid clusters are the same on both meshes. Although the pinch-off of liquid ligaments and drops from the liquid jet core and the surface tension are probably underresolved, this does not affect significantly the breakup process of the liquid core due to the following. First, the resulting ligaments and droplets move and disperse in the radial direction and thus do not influence the liquid jet core behaviour. Second, the breakup of the liquid jet core into clusters is still dominated by the aerodynamics rather than the surface tension due to the high Weber number based on the dimensions of the liquid clusters. The liquid core length predicted on Mesh II only shows marginal difference (5% longer) from that on Mesh I. Therefore, the breakup processes of the liquid column as a whole is well resolved on the current mesh, assuring the validity of this numerical scheme for prediction of liquid core length.

Conclusions
The primary breakup of a round water jet in coaxial air flow was simulated using a developed two-phase flow CLSVOF LES methodology. The Rescaling/Recycling Method (R 2 M) developed in Xiao et al. (2010) and Xiao (2012) was successfully implemented to generate unsteady turbulent inflows with appropriate statistical properties for both liquid and gas phases in simulations of air-blast atomisation. When these realistic turbulent inflows generated by R 2 M were used, the simulated jet structures agreed qualitatively with experimental shadowgraph images and the predicted liquid jet breakup length agreed quantitatively well with experimental    measurements. The simulation results indicated that the liquid jet breakup length was a power law function of the gas/liquid momentum flux ratio with a power law exponent from À0.39 to À0.35. By specifying different inlet conditions for liquid or gas flows, the effect of turbulent eddies developed inside the injector nozzles could be examined. It was found that liquid turbulent eddies, which will exist under the high Re conditions found in practice, are responsible for the initial interface perturbations. If laminar liquid inflow were to occur, the mean shear stress from the gas flow would then be the primary cause of the initial interface instability. The developed two-phase flow LES methodology showed excellent numerical robustness and accuracy in modelling liquid jet primary breakup at high liquid/gas density ratio.