A high-ﬁdelity simulation of the primary breakup within suspension high velocity oxy fuel thermal spray using a coupled volume of ﬂuid and discrete phase model

In the suspension high velocity oxy fuel (SHVOF) thermal spray, a suspension is injected into the com- bustion chamber where the jet undergoes primary and secondary breakup. Current knowledge of the primary breakup within the combustion chamber is very limited as experimental investigations are im- peded due to direct observational inaccessibility. Numerical methods are also limited due to the computational costs associated with resolving the entire range of multiphase structures within SHVOF thermal spray. This paper employs a coupled volume of ﬂuid and discrete phase model, combined with a combustion model, to simulate primary breakup at a fraction of the cost of a fully resolved simulation. A high-ﬁdelity model is employed within this study to model the combustion chamber; the model shows a backﬂow region that will contribute to clogging within the nozzle. This study modiﬁes the injector type for SHVOF thermal spray by introducing a co-ﬂow around the liquid injection to reduce clogging within the combustion chamber. This study shows that introducing a co-ﬂow of gas at a velocity of 200 m/s around the liquid injection reduces the backﬂow region by 40% within the combustion chamber. The addition of a gas co-ﬂow results in a smaller region of backﬂow. Small suspension droplets with insuﬃcient momentum are unable to overcome the backﬂow and will likely deposit themselves onto the wall of the combustion chamber. The deposition of the particles on the walls causes clogging of nozzles often seen in SHVOF thermal spray. The addition of a gas co-ﬂow results in an increase in the velocity of droplets formed during primary breakup. The greater droplet velocity allows for small droplets to overcome the small backﬂow region near the liquid injection. The Sauter mean diameters predicted from the numerical model are compared to experimental measurements available within the literature and shows good agreement.


Introduction
The breakup of liquid jets has been extensively studied both experimentally and computationally due to the wide range of applications that utilise atomization of liquid jets. The breakup of liquid jets plays a fundamental role within diesel engines, suspension thermal spray, gas turbines and rocket engines to name a few ( Xiao et al., 2014 ;Tian et al., 2015 ;Ghiji et al., 2017 ). Breakup investigations within a suspension high velocity oxy fuel (SHVOF) thermal spray combustion chamber through experimental methods are very limited due to the direct observational inaccessibility of the combustion chamber. Additionally, there are several challenges when trying to measure dense regions of flow. The liquid is sur- ( Esther Dongmo et al., 2009 ;Jadidi et al., 2016 ). Traditionally, the blob method is employed to model the injection of suspension into the combustion chamber ( Chadha et al., 2019b ;Chadha et al., 2019a ). For the blob method the jet is simplified to an injection of discrete droplets that can undergo secondary breakup. The suspension is typically injected as a set of discrete droplets using the discrete phase model (DPM) that can undergo secondary breakup, evaporation, heating and acceleration. The DPM framework provides a robust treatment for droplet tracking away from the dense regions of the flow; however, in dense regions of the flow such as near the liquid column; this approach is flawed.
Traditional means of modelling primary breakup in CFD using methods such as interface reconstruction methods are unfeasible to resolve the full scale of droplet structures within SHVOF thermal spray. The most popular interface reconstruction method, the volume of fluid (VOF) model, requires a mesh fine enough to resolve the smallest droplet structures. Droplets formed from the primary breakup of jets occupy a wide range of scales and are of orders of magnitude smaller than the injector diameter ( Tomar et al., 2010 ). Hence, it would be computationally expensive to resolve the breakup phenomenon from large scale structures down to the smallest scale droplets. Within recent years there has been a growing body of literature looking to couple the VOF and DPM models together. Here, the primary breakup is modelled using the VOF model and the droplets formed from primary breakup are then transferred to a DPM framework. The small-scale structures that require the most computational resources within the VOF framework can then be modelled at a lower computational cost using the DPM framework. This allows for large scale structures such as the liquid core, ligaments and large droplets to be resolved whilst allowing the small-scale droplet structures to be modelled without significantly impacting the computational cost of the model.
Early attempts to couple the VOF model with the DPM model had to address the criteria of conversion of droplets from a VOF framework to a DPM framework. Grosshans et al. (2011) utilised a "coupling layer" where VOF droplets are converted into DPM droplets as they pass through a plane known as the coupling layer. One of the challenges in employing a coupled VOF and DPM model for breakup investigations is the different mesh requirements within the VOF and DPM frameworks. With the implementation of a coupling layer suitable meshes can be employed within the respective zones. For the VOF framework droplets must be of an order larger than the mesh size. Whilst for the DPM framework the mesh should be of an order larger than the droplets. The coupling layer approach however is significantly more expensive than later alternative approaches that have been employed to couple the two frameworks. Adeniyi et al. (2017) developed a coupled DPM to VOF framework for bearing chamber applications using the commercial CFD code Ansys Fluent. Within this implementation, droplets were modelled using the DPM framework and the film formation on the bearing chamber walls was modelled using the VOF framework. The criteria for conversion between DPM and VOF was the proximity of DPM droplets to the interface. Tomar et al. (2010) employed a two-way coupled VOF and DPM framework to model jet breakup. The conversion of VOF droplets to DPM droplets was determined by the diameter of the VOF droplets. If the droplets were smaller than the specified diameter then the droplets would convert to DPM. This approach can offer greater savings in computational cost over the coupling layer method as the droplets can be switched from a VOF framework to a DPM framework as soon as the droplets are formed. The conversion of droplets from DPM to VOF was based on the proximity of the droplets to an interface. If the droplets are close enough to an interface the DPM droplets will be deleted and a spherical secondary phase structure is patched in its location. Kim et al. (2014) and Kim et al. (2007) employed a coupled VOF and DPM ap- proach to model the atomization of fuel within a gas turbine injector. Within this study droplets were converted from a VOF framework to a DPM framework if two criteria were met. The first criterion was based on the droplet volume and if the droplets were smaller than a specified volume, they were converted from VOF to DPM. The second criterion was based upon the droplet sphericity and if the droplets were sufficiently spherical the droplets were converted from VOF to DPM. Shinjo and Umemura (2019) and Umemura (2016) extended the application of the coupled VOF and DPM model to include the effects of combustion for a diesel jet. They investigated four cases: diesel jet in cold flow, diesel jet in hot flow, with the inclusion of a single step reaction and finally the inclusion of a multistep reaction. The multistep reaction demonstrated good agreement in predicting the ignition delay with existing literature.
The aim of this paper is to apply the coupled VOF and DPM model to investigate the primary breakup of water within a combustion chamber for SHVOF thermal spray application. The injector design within this study is varied from the standard SHVOF configurations by including a co-flow around the liquid injection. There are several non-dimensional parameters that describe the flow in coaxial jets these include the gas Reynolds number, Re g , the liquid Reynolds number, Re l , the momentum flux ratio, MR, and the Ohenesorge number, Oh. This paper investigates the impact of the Weber number ( Eq. (1) ) by varying the co-flow velocity on the primary breakup within SHVOF thermal spray. The Weber number controls the tendency for breakup ( Xiao et al., 2014 ;Pai et al., 2009 ) and characterises the relative importance of the fluid inertia to its surface tension ( Xiao et al., 2016 ). The terms in Eq. (1) refer to liquid density, ρ l , gas velocity, U g , liquid velocity, U l , diameter of injector, d , and the surface tension, σ . The addition of a co-flow provides a significant design change for thermal spray.
This paper has employed a multiscale modelling approach for a real-world engineering application. Prior numerical models within this field of application ignore the primary breakup which is a key physical process in the combustion chamber. The detailed insight into the flow physics within the combustion chamber for the current application is new and has been previously unknown. This simulation approach utilised here has never before been used for a practical combustion application such as SHVOF. The approach can be employed on modest HPC facilities for design and optimization studies within a combustion chamber and also for a wide array of multiphase combustion applications. An injection of water into the combustion chamber is considered to provide useful and highly detailed information on the physical processes for the current application. This paper uses the Coupled VOF and DPM model in the commercial CFD software Ansys Fluent V19.3 (Pennsylvania, USA).

Numerical modelling
A SHVOF thermal spray nozzle is comprised of a combustion chamber and a barrel as shown in Fig. 1 . Typically when modelling SHVOF thermal spray, the entire nozzle and the free jet is modelled. However, the combustion chamber alone has been modelled for this investigation, as the primary breakup is confined to the combustion chamber. This allows for a significant reduction in the overall mesh count. Fig. 2 shows the fully structured mesh and the inlet geometry employed to model the combustion chamber; the mesh is comprised of 4.5 million cells before mesh refinement. The initial mesh was created following LES near-wall resolution recommendations of x + < 100, y + < 1 and z + < 30. The structured mesh was generated using a multiblock method within Ansys ICEM. To reduce the computational cost, an adaptive mesh refinement is employed here to increase the mesh density around the interface between the two phases. The refined mesh is coarsened back to the base mesh once the droplet has been converted from VOF to DPM. The base mesh is of the order of 100 μm which is twice as large as the largest DPM droplets and the refined mesh is off the order of 3 μm which is smaller than the smallest VOF droplets. This allows for a finer mesh to be employed for VOF resolved droplets and a coarser mesh to be employed with DPM modelled droplets. Five levels of dynamic mesh refinement have been employed with the gradient of the volume fraction. Premixed oxygen and hydrogen enter the combustion chamber through the 24 circular holes at the back end of the combustion chamber, the hydrogen undergoes combustion. The combustion occurs only in the gas phase. Water is injected into the centre of the combustion chamber. The impact of solid particles, particle size distribution, shape of particles and the mass loading of particles on liquid suspension surface tension properties is not well understood. Further extensive experimental measurements are needed for detailed characterisation before extending the breakup model to a suspension jet. The governing equations of the methods are provided in the following section. The boundary conditions for the case are provided in Table 1 . The operating conditions employed are standard conditions for SHVOF thermal spray ( Chadha et al., 2019b ). The mass continuity equation is given by Eq. (2) . For a density, ρ, time, t and velocity vector u = ( u, v , w ) T , it reads: The momentum conservation equation is given by Eq. (3) where the terms within Eq. (3) refer to the static pressure, p , a dynamic viscosity, μ, the unit tensor, I and a gravitational body force, ρg .
A wall resolved large eddy simulation (LES) approach is employed. The small-scale structures are distinguished from the largescale structures using a filter. The filter is given by Eq. (4) and applied to the Navier-Stokes equations. The wall-adapting local eddy viscosity (WALE) SGS model is used and defined by Eq. (5) . The mixing length for subgrid scales, L s , is given by Eq. (6) where κ represents the von Karman constant, d represents the nearest wall distance, C w represents the WALE constant, and V represents the cell volume. The rate of strain for the resolved scales, S d i j , is given by Eq. (7) .
The VOF method constructs the interface between two phases by solving the volume of fluid, φ, of each phase given by Eq. (9) .
The ideal gas law is given by Eq. (10) where p op , refers to the operating pressure, R , refers to the universal gas constant, M w , refers to the molecular weight of the gas, and T refers to the temperature ( Jafari et al., 2017 ). The energy, E , equation is given by Eq. (11) , where the effective conductivity of the gas is given by k eff , and J i represents the diffusion flux for species i, h i represents the sensible enthalpy for species i, and the heat source term is given by S h . The species equation is given by Eq. (12) , where Y i , refers to the mass fraction of species i , and R i represents the rate of production of species i .
For turbulent premixed combustion the rate of reaction is limited by the turbulent mixing of cold premixed fuel and oxidizer with the hot gaseous products ( E. Dongmo et al., 2009 ). The eddy dissipation model (EDM) has been employed due to its low computational cost. The current numerical modelling approach is intended to be suitable for practical design and optimization studies the computational cost needs to be carefully considered. The EDM was developed to account for the effects of turbulent mixing on the reaction rate. The reaction rate within the eddy dissipation model is given by the smaller of the two terms in Eq. (13a) and (13b) . Here the stoichiometric coefficient for reactant, i, in reaction, r , is given by, ν i,r , the molecular weight of species, i , is given by M W,i , and A and B refer to empirical constants of 4.0 and 0.5 respectively. The quantities k and ε refer to the turbulent kinetic energy per unit mass and the dissipation rate respectively. The quantities Y R and Y p refer to the mass fraction of a reactant, R , and the mass fraction of product species, P . The combustion reaction accounting for the dissociated species is determined from the chem-ical equilibrium package NASA CEA ( Chadha et al., 2019a ) and is given by Eq. (14) .
VOF resolved droplet are converted to DPM modelled droplets based upon their size and sphericity. Values for the radius standard deviation and the radius surface orthogonality of 0.5 and 0.5 were used respectively. The volume equivalent sphere diameter is a user specified range. The lower value was set to zero to ensure that the smallest scales converted to DPM. Care was taken to ensure the upper limit was larger than the largest expected droplet diameter. The DPM model tracks droplets as they are formed from the primary breakup process. The motion of the droplets is given by Newton's second law, Eq. (15) . Here, the mass of the discrete phase is given by, m p , the gas velocity is given by, u g , and the particle velocity is given by, u p . The drag coefficient, c D , is determined by the correlation of Crowe ( Crowe, 1967 ) given by Eq. (16) . The correlation of Crowe has shown to provide good agreement of particle velocities for suspension high velocity oxy fuel thermal spray within the free jet ( Jadidi et al., 2016 ;Chadha et al., 2020 ). The particle temperature, T P , can be determined from Eq. (17) , where the particle surface area is given by, A P and the gas temperature is given by, T g . The heat transfer coefficient, h , is computed using the Nusselt number, Nu , correlation of Ranz and Marshall ( Ansys Inc, 2019 ) given by Eq. (18) . The diameter of the particle is given by, d p , the thermal conductivity of the gas is given by, K 0 , the dimensionless Prandtl number is given by, Pr and the particle Reynolds number is given by Re d .
Similar DPM modelling approaches to that employed within this study have been validated by our recent study ( Chadha et al., 2020 ) and a similar study by ( Jadidi et al., 2016 ). The VOF model has been extensively used simulate breakup and atomization of liquid jets, for example Sekar et al. (2014) , Delteil et al. (2011) and Srinivasan et al. (2011) . This paper brings both modelling approaches together and employs the model to investigate the breakup within SHVOF thermal spray. To convert structures from a VOF framework to a DPM framework, a lump detection algorithm is employed within Ansys Fluent. The detection algorithm scans the fluid domain searching for enclosed droplet structures. The lumps are then assessed on their asphericity based upon the standard deviation of the droplet radius and the radius-surface orthogonality. If the droplets meet the criteria for conversion, they are switched from a VOF framework to a DPM framework. The cells containing the droplet are patched with the primary phase to ensure a continuity of the volume fraction. The mass of the droplet with its equivalent diameter is then injected into the cell using the DPM framework.

Results and discussion
Effect of co-flow on the gas velocity and temperature within the combustion chamber Traditionally in SHVOF thermal spray modelling the inlets within the combustion chamber are simplified to an annular inlet to simplify the geometry, this allow for a reduction in the cell count ( Chadha et al., 2019b ;Jadidi et al., 2018 ). Within this study the exact inlet geometry comprising of two sets of twelve circular holes spaced at two different radial distances from the centre of the combustion chamber have been meshed. Additionally, it can be seen in Figs. 3 and 4 that recirculation zones form near the walls of the combustion chamber. The recirculation zones aid in the mixing of the hot products of combustion with the cold inlet gases, which is essential for a premixed combustion reaction to progress effectively. Within premixed combustion the reaction takes place within the zone that separates the reactants and products.
Clogging is a serious issue for suspensions with a high particle loading within SHVOF thermal spray. Clogging wastes feedstock material, cost's time to restore and maintain the SHVOF thermal spray nozzles. Clogging occurs when the feedstock material impacts the wall with a velocity greater than a critical velocity ( Bémer et al., 2013 ). Particles that impact the wall with a velocity greater than the critical velocity will bond onto the walls of the nozzle. Particles within small droplets are particularly susceptible to causing clogging as these droplets have a low mass inertia. Small droplets are unable to overcome the backflow region near the liquid column shown in Fig. 6 due to insufficient momentum. Such droplets can deposit themselves onto the walls of the combustion chamber due to backflow illustrated within the velocity vector field in Fig. 5 . Optimizing the suspension injection will allow for a significant reduction in the number of particles that impact with the walls of the combustion chamber. This study has employed a co-flow around the liquid injection to force the flow field around the liquid column in the direction of nozzle exit, as shown in Fig. 2 . With the addition of a co-flow the region of backflow becomes smaller hence droplets are less likely to deposit upon the wall of the combustion chamber and thus reduces the risk of the nozzle clogging.
The coupled VOF and DPM framework employed within this study is able to capture the recirculation zones that form around the liquid column due to the viscous flow dynamics as seen in Fig. 4 . Recirculation zones form due to the velocity gradient between the low velocity liquid injection and the high velocity fuel injection. Prior SHVOF thermal spray modelling has been unable to capture the recirculation zones with the standalone DPM model. This is because the pure DPM model only accounts for momentum transfer between the discrete phase and the continuous phase through a momentum source term. The VOF model includes a source term within the momentum equation to account for the effect of surface tension forces. Further discussion on the bene- fits of this approach in comparison to the DPM model employed currently to better capture the breakup physics is later presented. Fig. 6 shows instantaneous static temperature contours for the Weber numbers of 100, 180 and 415 in Figures A, B and C respectively. It can be seen from Fig. 6 that within the liquid injection region there is a very low temperature of ∼ 300 K. With a standalone DPM model, currently employed within the literature for this application ( Chadha et al., 2020 ;Gerhardter et al., 2019 ), a temperature of 2500 K is predicted within this region. The DPM model only accounts for the effect of the secondary phase on the temperature field through a source term for droplet heat transfer. The VOF model on the other hand directly accounts for the effect of the secondary phase through the Navier Stokes equations. The DPM has been shown to give an over-prediction of gas temperature and the velocity within the liquid injection region ( Chadha et al., 2019a ) as the dilute approximation is not valid within the regions near the injector which will result in evaporation initiating earlier and locally higher droplet velocities. Hence employing the VOF in these regions is able to capture more of the flow physics and better predict the temperature near the liquid injection.

Mesh requirements for coupled VOF and DPM breakup modelling
The VOF model and the DPM model have very different mesh requirements. The VOF model requires the mesh to be much smaller than the droplet diameter, enough cells are required to resolve the curvature of the droplet. On the other hand, the DPM model requires the mesh to be large than the droplet diameter. Recent investigations have looked into modifications to the DPM model to make it suitable for fine meshes. To account for the effect of the discrete phase onto the continuous phase source terms are accounted for within the Navier-Stokes equations. The DPM source term is applied to the cell which contains the DPM point entity. However, once the DPM droplet exceeds the cell size the droplet will affect the more than just the cell that contains the droplet. Due to the point representation of the droplet the source term in the original DPM model does not account for the effect of the droplet on neighbouring cells.
There are several approaches that have been developed within the literature to resolve this flaw within the DPM model ( Zhang et al., 2020 ). The cube averaging method (CAM) developed   by ( Link et al., 2005 ) creates a cubic region a factor larger than the droplet diameter. The droplets are treated as porous cubes within the cubic region where the source terms are calculated and distributed over the cubic region. The source terms are then converted back from the cubic region and mapped to the original grid. The two-grid method (TGM) developed by ( Farzaneh et al., 2011 )  proposed to use a coarse grid for the discrete phase and a fine grid for the continuous phase. The source terms are mapped from the coarse grid to the fine grid and the source terms are weighted by the volume of the fine grid to the coarse grid. The diffusion based method (DBM) developed by ( Capecelatro and Desjardins, 2013 ) distributes source terms from the discrete phase to the continu- ous phase with the use of a weighted function also referred to as a statistical kernel function. Source terms can be distributed over kernels such as a Gaussian kernel which employs a Gaussian distribution for the source term within the cell the DPM droplet is contained and the neighbouring cells.
Within this study an aggressive adaptive mesh refinement algorithm is employed to account for the different mesh resolutions required for VOF resolved droplets and DPM modelled droplets. Fig. 7 shows the mesh density around the VOF interface. As can be seen the mesh is suitably refined to accurately capture the interface. The refined mesh is immediately coarsened back to the base mesh once the droplet has been converted from VOF to DPM. The base mesh is of the order of 100 μm which is twice as large as the largest DPM droplets and the refined mesh is off the order of 3 μm which is three times smaller than the smallest VOF droplets.
Allowing for a finer mesh to be employed to resolve the interface between phases. The adaptive mesh refinement also allows for a significant reduction in the computational cost of the numerical model as a much coarser base mesh can be employed. A fine mesh is only employed within the cells where the interface is found.

Effect of the co-flow on the primary breakup
Prior studies within the SHVOF thermal spray modelling have solely employed the DPM framework to model the injection of suspension within the combustion chamber ( Chadha et al., 2020 ;Jadidi et al., 2016 ;Gozali et al., 2014 ). The DPM model has typically been employed due to its low computational cost, robust treatment of droplets through sub models such as the evaporation and secondary breakup sub models. However, the DPM model produces an inadequate representation of the flow near the injection location where a liquid column and ligaments are the main structures present. The assumptions within the DPM model are that the multiphase structures are in the form of droplets and that the region of the flow is predominately occupied by the primary phase. Both of these assumptions near the liquid injector are not valid. The VOF model is a more suitable model to capture the flow near the liquid injector as it can capture all multiphase structures that form and is suitable for the dense regions. However, resolving the entire range of multiphase structures using a pure VOF model would be too computationally expensive for SHVOF thermal spray applications. Droplets and particles within SHVOF thermal spray can occupy a diameter as small as 10 μm ( Chadha et al., 2020 ). Resolving droplets and particles of this size would require a mesh around the order of one fifth of the diameter of the finest droplets. Employing a DPM framework to model droplets within the combustion chamber allows for a coarser mesh to be employed and a reduction in the computational cost in comparison to a pure VOF model. The current approach uses a 4.5 million base mesh and through mesh adaptation reaches a mesh count of circa 8 million cells. Running a full VOF simulation would require a mesh resolution down to the 3 μm scale currently used but in a much larger volume. To resolve droplets to this extent would require a mesh of the order of 20 billion cells which is currently infeasible. This could be alleviated to some extent through the use of adaptive mesh refinement, however, given the number of droplets involved this is unlikely to be efficient. Figs. 8 shows isosurfaces of the VOF field at the Weber numbers of 100, 185 and 415; the isosurface is coloured by velocity magnitude. The coupled VOF and DPM model has allowed for the primary breakup to be modelled and droplets to be tracked with current computing capabilities available. Therefore, the coupled VOF and DPM approach employed within this study is able to accurately capture the flow field and temperature field within and near the liquid column region compared to the DPM approach currently employed. Alternative models are available within the literature to treat the liquid column region such as the Eulerian-Lagrangian Spray Atomization (EL SA) method. The EL SA method solves transport equations for the gas and liquid phases treating them as separate species rather than separate phases ( Hoyas et al., 2013 ). The model cannot resolve breakup structures such as ligaments nor can the model resolve the interface between phases like interface resolving methods such as the VOF method employed for large scale structures within this model. The ELSA method has a lower fidelity than the approach employed within this investigation however it is considerably less computationally expensive.
It has been shown that the addition of a co-flow results in a smaller backflow region and hence suggested that this will result in a reduction in the clogging within the combustion chamber. It must be further evaluated if any further effects of the co-flow are present on the breakup and droplet dynamics. Fig. 8 shows the instantaneous isosurfaces of the liquid jet injection into the combustion chamber (We  Fig. 8 that the addition of a co-axial gas flow to the liquid injection delays the formation of large ligaments from the surface of the liquid jet. This can be seen in Fig. 8 by comparing the isosurface plots for a Weber number of 100 and a Weber number of 180 close to the injector inlet. As can be seen the liquid jet core is destroyed in a shorter distance for the lower Weber number. This is the opposite trend to that normally observed for liquid jets; however, it should be noted that the gas within the combustion chamber is hot and highly turbulent. It is thought that its interaction with the liquid jet excites the initial interface perturbation followed by large scale interface distortion. The recirculation zones set up within the combustion chamber also mean that the surrounding gas flows in an opposite direction to the liquid jet. The addition of the gas co-flow rights this direction and is of sufficiently lower turbulence intensity and temperature. Hence, there is a delayed break-down of the jet core. Comparisons of the jet length at different Weber numbers can be used to benchmark the numerical accuracy. These have not been included as experimental measurements are taken at conditions quite different to those found within an SHVOF thermal spray combustion chamber. It is not possible to provide a meaningful comparison of the intact core jet length from the numerical predictions with measurements from the literature due to the interaction of the jet and the hot and highly turbulent combustion chamber gases. The addition of the co-axial gas flow results in an increase in the jet length and a higher velocity of the liquid jet. The liquid jet in Fig. 8 (a) and (b) shows a velocity not exceeding 20 m/s, while Figs. 8 (c), (d), (e) and (f) show the liquid jet velocity well exceeding 30 m/s. The higher jet velocity allows for the liquid jet, ligaments and the droplets formed to overcome the smaller region of back flow within the combustion chamber more readily. The mechanism behind the primary breakup of liquid jets in co-axial gas flow has been well documented within the primary breakup literature. Surface perturbations are formed on the interface between the liquid and gas by primary shear instability. Surface perturba- tions are stretched into ligaments due to the relative velocity difference between the two fluids. As the ligaments are accelerated, support from the bulk liquid is diminished and the ligaments elongate. The surface of the ligaments are continually subject to strong accelerating forces from the gas jet leading toacknnnnnn Rayleigh-Taylor instabilities. As the Rayleigh-Taylor waves amplify the ligaments, they break off the main jet forming a droplet. Fig. 9 (a), (b) and (c) show droplet diameter distributions for droplets that form from the primary breakup of the liquid jet at the Weber numbers of 100, 185 and 415 respectively. It can be seen that the droplet diameter distributions at the three Weber numbers follow a log normal distribution. It is reported by ( Kazuya et al., 2004 ) that log-normal distributions are typically seen within primary breakup. There is a peak frequency at a diameter and as the droplet diameter decreases from the peak value the frequency of the droplets decreases. A cut-off diameter is seen at a diameter smaller than the peak value as surface tension forces prevent the formation of smaller droplets. Table 2 shows the mean and standard deviation of the log-normal distributions from the histograms in Fig. 9 . It can be seen from table 2 that at the Weber numbers considered that there is little variation in the mean and standard deviation of the droplet diameters.
Experimental observations and measurements within the combustion chamber are very limited due to lack of optical access within the combustion chamber. Therefore, experimental measurements taken at similar Weber numbers, liquid injector diameter and co-flow velocity can allow for an understanding if the numerical values are in a range similar to that expected. Fig. 10 compares the Sauter mean diameter (SMD) for the droplets from primary breakup compared with the experimental measurements of ( Varga et al., 2003 ). It can be seen from Fig. 10 that the droplet diameters predicted within this study match up very well to those measured by Varga et al. through phase-doppler anemometry (PDA) at Weber numbers of 180 and 415; however, there is a small underprediction in the droplet SMD at the Weber number of 100. The measurements taken by Varga et al. considered co-axial injection of liquid and gas flow where both phases travel in the same direction. For the case of where the liquid jet Weber corresponds to 100 there is no co-flow with the liquid injection. The gas flow is in the direction opposing the liquid injection. Hence, this may have some effect on the mean diameter which leads to the small under-prediction witnessed.
Droplets with a higher velocity will be able to more readily overcome the backflow region in the combustion chamber due to their greater momentum. Fig. 11 shows the relationship between the droplet diameter and the droplet velocity from droplets formed from the primary breakup of the liquid column at Weber numbers of 100, 180 and 415. It can be seen that as the Weber number increases the droplet velocity increases. With a higher co-flow velocity greater momentum is transferred to the liquid from the gas allowing for the droplets to obtain a higher velocity. It can also be seen in Fig. 11 that as the droplet diameter increases the velocity of the droplets decrease. The decay takes an exponential form for all the Weber numbers considered. The distributions shown in Fig. 11 can be represented by Eq. (19) , where a, b and c are coefficients that are fitted to the data for the three Weber numbers. The equation is based on an exponential curve using MATLAB's curve fitting toolbox. The three curves at the different Weber numbers take a similar exponential form. As can be seen, the decay constant 'b' is similar for all three cases taking a value of around 2.0. This would indicate that it is independent of the Weber number. The vertical shift in curves appears a function of Weber number and is controlled by parameter 'c'. Therefore, the addition of a co-flow gas injection provides a benefit of a smaller backflow region and a higher velocity for droplets. Both of which will reduce clogging within the combustion chamber. Fig. 12 shows the liquid column and the droplets within the combustion chamber at Weber numbers of 100, 180 and 415. It can be seen from Fig. 12 that the higher droplet velocity results in a reduction in the number of droplets that are trapped within the recirculation zones around the liquid column. This will in turn will aide in the prevention of clogging of suspension within the combustion chamber as fewer droplets are trapped within recirculation zones within the combustion chamber and therefore, solids within the suspension are less likely to deposit onto the combustion chamber walls.
The droplet Weber number can determine the mode of secondary breakup of the droplets and hence the droplet Weber number subsequent to primary breakup has been evaluated. Fig. 13 shows the effect of co-flow velocity on the Weber number of droplets produced from primary breakup. The droplet diameter distribution in Fig. 9 is used to determine the Weber numbers. It can be seen from Fig. 13 (a), (b) and (c) that as the co-flow velocity increases the spread of droplet Weber numbers increases. In addition to this the distribution of Weber numbers shifts slightly towards the right where droplets occupy a higher Weber number at a higher co-flow velocity. The mode of secondary breakup will be different as the co-flow velocities increases. With no co-flow the droplets have a very low Weber number and hence the mode of secondary breakup will be dominated by vibrational breakup and bag breakup. While for a 200 m/s co-flow velocity the secondary breakup will be dominated by bag and stamen breakup. Finally, at a co-flow velocity of 300 m/s the mode of secondary breakup will be dominated by bag and stamen breakup and stripping breakup.
The coupled VOF and DPM has been employed to model the primary breakup of liquid within the combustion chamber within SHVOF thermal spray. This hybrid approach provides many advantages over the traditional standalone multiphase models. The hybrid VOF and DPM approach allows for significant reduction in the computational cost over a standalone VOF model. Additionally, this approach provides considerably higher fidelity than the DPM model. Within this study the effect of introducing a co-annular injection on the primary breakup and the flow field has been investigated. The SMD from primary breakup is compared to experimental measurements at the same Weber numbers. The numerical model is shown to be in good agreement with the experimental measurements available within the literature.

Conclusion
In summary, this investigation is the first study to employ a coupled VOF and DPM approach to model the primary breakup with combustion for a real world engineering application. This study has investigated the effect of using a gas co-flow injection to inject the liquid into the combustion chamber and has shown: • The addition of a gas co-flow results in a smaller region of backflow. Small suspension droplets with insufficient momentum are unable to overcome the backflow and deposit themselves onto the wall of the combustion chamber. The deposition of the particles on the walls is likely to be a cause of clogging of nozzles often seen in SHVOF thermal spray.
• The addition of a gas co-flow results in a significantly higher velocity for the droplets formed during primary breakup. The greater droplet velocity allows for small droplets to overcome the small backflow region near the liquid injection. • The droplet diameter distributions from the primary breakup have been obtained which follow a log normal distribution. • The Sauter mean diameter (SMD) is compared to experimental measurements available within the literature for the same Weber numbers and it is seen that the SMD predicted from the numerical model match closely to experimental data within the literature that is obtained using phase-doppler anemometry (PDA). • A correlation for the relationship between the droplet diameter and the velocity during the primary breakup of the liquid jet has been established. • The work presented in this paper uses high-fidelity scale resolving LES and employs the VOF model for primary breakup coupled with DPM for the smallest scales. This is combined with a combustion model for the first time. It provides a significant enhancement over standalone DPM that has been traditionally used for this application.

Declaration of Competing Interest
None.