# CFD Simulations AC3-08

**Spray evaporation in turbulent flow**

**Application Challenge 3-08** © copyright ERCOFTAC 2004

**CFD Simulations**

The evaporating spray was considered as a test case for the 6th Workshop on Two-Phase Flow Predictions (Sommerfeld 1993). Several research groups have participated in the test case calculations, mostly using in-house codes. Additionally, calculations of this test case were presented by Sommerfeld (1998) and Kohnen and Sommerfeld (1998). After a review of the numerical method employed, some of the numerical results are presented.

**Overview of CFD simulations and solution strategy**

The numerical calculations of spray droplet motion and evaporation in a turbulent flow were based on the Eulerian/Lagrangian approach (in-house 2D code) for the gas and droplet phase, respectively. The gas phase was computed by solving the axi-symmetric form of the time-averaged conservation equations in connection with the well-known k-ε turbulence model of Launder and Spalding (1974) employing the standard constants. The conservation equations consist of the continuity, momentum, temperature and species equations, which also include source terms resulting from the dispersed phase heat, mass and momentum transfer. Two-way coupling in the k and ε equations was also considered by using the source terms obtained through Reynolds averaging (Kohnen 1997).

The set of gas-phase conservation equations is solved by using a finite-volume discretisation scheme and applying an iterative solution procedure based on the SIMPLE algorithm. The convective terms were discretised by a flux correction method (a combination of upwind and central differencing) and the diffusive terms were discretised by second central differences.

The droplet phase was treated by the Lagrangian approach, where a large number of droplet parcels, representing a number of real droplets with the same properties, were traced through the flow field. The representation of the droplets by parcels was used in order to allow the consideration of the droplet size distribution and to simulate the appropriate liquid mass flow rate at the injection locations. The droplet parcels were traced through the flow field by solving a set of ordinary differential equations for the droplet location, diameter, velocity components, and temperature. For the formulation of the droplets equation of motion only drag and gravity were taken into account. Moreover, the equations of motion were solved in a Cartesian co-ordinate system in order to avoid the singularity problem at the centreline. The instantaneous fluid velocity in the above equations (i.e. along the droplet trajectory) was determined using the Discrete Eddy Concept and a drift correction for the transverse direction as it was described by Sommerfeld et al. (1993). The droplet evaporation was calculated using the infinite conductivity model described by Abramzon and Sirignano (1989).

The source terms in the gas phase equations for mass, momentum, and heat exchange were evaluated for each control volume by averaging a number of droplet trajectories. The calculation of the source terms was based on a modified PSI-CELL approach. In the momentum and energy source terms, the effects of both, evaporation and possible condensation, were taken into account. For the k and ε equations the source terms include effects due to mass and momentum transfer. The Eulerian and Lagrangian part were solved sequentially using an under-relaxation for the source terms (i.e. g = 0.1), until a converged solution for both phases was obtained. This was normally achieved after 50 coupling iterations. More details about the numerical method can be found in Kohnen and Sommerfeld (1998) and Sommerfeld (1998).

**Computational Domain**

The computational domain had a dimension of 0.1 m x 1.2 m in the radial and axial directions, respectively. Four different numerical grids were used with 13 x 22, 24 x 42, 46 x 82 and 90 x 162 cells in the radial and axial direction (i.e. 0.1 m x 1.5 m). A grid-independent result was obtained with 90 x 162 cells. The computations were performed on an equidistant grid in the radial direction and continuously expanded in the axial direction with the finest mesh near the inlet.

**Boundary Conditions**

The stream-wise air velocity profile at the annular inlet (20 < r < 32) was specified according to the measurements, which are available for the single phase flow and the different spray cases. The radial and tangential air mean velocity was assumed to be zero. For the determination of the turbulent kinetic energy at the air inlet, measurements for all the components are available. The dissipation rate (ε) at the inlet was prescribed according to:

where c_{μ} = 0.09 and Δr is the width of the annulus,
i.e. Δr = 12 mm. At the walls, all velocity components were set to zero and the standard law of the wall was applied in the κ-ε turbulence model. For the air temperature a constant value was prescribed at the inlet. Additionally, the mean wall temperature along the test section is available, which was used as a boundary condition for solution of the temperature transport equation. However, the influence of uncertainty in the inlet conditions on the computations was not analysed in detail.

The injection of the droplets was based on the measured properties at several radial positions 3mm downstream of the nozzle holder (Fig. 1). For this purpose profiles of the mean properties, such as velocity, rms values and droplet mass flux are available. Moreover, droplet size distributions and the correlations between mean and rms velocities and the droplet size are available for radial locations 0< r < 9 mm in intervals of 1 mm. Hence, when a droplet is injected at the inlet, its location is randomly sampled. The size is sampled from the measured size distribution. According to this size, the mean velocities are obtained and the fluctuating velocity components are sampled from normal distribution functions with the respective rms values.

The particle size spectrum was resolved by 20 classes having a width of 5 mm. In order to get statistically reliable results 42,000 parcel trajectories were calculated during each of the required 50 coupling iterations (i.e. sequential calculation of the Eulerian and Lagrangian part, see Kohnen 1997).

**Application of Physical Models**

The numerical calculations were based on the Euler/Lagrange approach (in-house 2D code) for the gas and droplet phase, respectively. The gas phase was computed by solving the axi-symmetric form of the time-averaged conservation equations in connection with the well-known k-ε turbulence model of Launder and Spalding (1974) employing the standard constants. Two-way coupling was accounted for in all conservation equations of the gas phase and the species equations.

The droplet phase was treated in a Lagrangian way, where a large number of droplet parcels, representing a number of real droplets with the same properties, were traced through the flow field by accounting for the droplet size distribution. In the formulation of the droplets equation of motion only drag and gravity were taken into account. The instantaneous fluid velocity along the droplet trajectories was determined using the Discrete Eddy Concept and a drift correction for the transverse direction as it was described by Sommerfeld et al. (1993). The droplet evaporation was calculated using the infinite conductivity model described by Abramzon and Sirignano (1989).

**Numerical Accuracy**

The Eulerian and Lagrangian part were solved sequentially using an under-relaxation for all the source terms (i.e. g = 0.1), until a converged solution for both phases was obtained. This was normally achieved after 50 coupling iterations depending on the degree of coupling. A grid-independent result for both phases was obtained with 90 x 162 cells in the axial and radial direction.

© copyright ERCOFTAC 2004

Contributors: Martin Sommerfeld - Martin-Luther-Universitat Halle-Wittenberg

Site Design and Implementation: Atkins and UniS