Eulerian–Lagrangian CFD-microphysics modeling of a near-field contrail from a realistic turbofan

Abstract

Aircraft contrails contribute to climate change through global radiative forcing. As part of the general effort aimed at developing reliable decision-making tools, this paper demonstrates the feasibility of implementing a Lagrangian ice microphysical module in a commercial CFD code to characterize the early development of near-field contrails. While engine jets are highly parameterized in most existing models in a way that neglects the nozzle exit-related aspects, our model accounts for the geometric complexity of modern turbofan exhausts. The modeling strategy is based on three-dimensional URANS simulations of an aircraft nozzle exit involving a bypass and a core jet (Eulerian gas phase). Solid soot and ice particles (dispersed phase) are individually tracked using a Lagrangian approach. The implemented microphysical module accounts for the main process of water-vapor condensation on pre-activated soot particles known as heterogeneous condensation. The predictive capabilities of the proposed model are demonstrated through a comprehensive validation set based on the jet-flow dynamics and turbulence statistics in the case of compressible, turbulent coaxial jets. Simulations of contrail formation from a realistic nozzle-exit geometry of the CFM56-3 engine (short-cowl nozzle delivering a dual stream jet with a bypass rate of 5.3) were also carried out in typical cruise flight conditions ensuring contrail formation. The model provides reliable predictions in terms of the plume dilution and ice-particle properties as compared to available in-flight and numerical data. Such a model can then be used to characterize the impact of nozzle-exit parameters on the optical and microphysical properties of near-field contrails.

Keywords

Contrails CFD microphysics coupling aircraft engine Eulerian–Lagrangian approach ice particle growth

Introduction

The continuous growth in air passengers by about 5% per year has increased scientific concerns about the climate impact of commercial aviation.¹ Cruise flights contribute to global warming through soot and water-vapor emissions from aircraft engines that are responsible for visible white condensation trails (contrails) behind aircraft.^2,3 Persistent contrails composed mainly of supersaturated ice crystals survive in the atmosphere for several hours and even evolve to form large cirrus formations referred to as aircraft-induced clouds (AIC).⁴ The spreading and persistence of this artificial cloudiness in the atmosphere determine their impact on the earth’s radiative balance.⁵ Global assessments⁶ show that the impact of linear contrails is about 3–9·10⁻³ W·m⁻², while the impact of AIC is 10 times higher,⁶ which corresponds to three times the impact of aviation CO₂ alone.⁷ For this reason, the contrail formation process has been a topic of increasing interest in the literature. Owing to both complexity and onerous costs of in-flight measurements, the development of innovative mitigation strategies has to be based on reliable decision-making tools. As such, several modeling approaches with different levels of microphysics/dynamics representations have been used in literature to investigate contrails. The following is a short review of main approaches and findings.

Modeling contrail life cycle involves several scales ranging from local scales in the engine jet phase (t ~ 0–10 s after the exhaust emission) to the global scale in the dispersion phase (few minutes to hours after the exhaust emission; see Figure 1). From the thermodynamic point of view, the revised Schmidt–Appleman criterion⁸ helps to determine atmospheric regions of contrail occurrence based on water liquid saturation conditions. This criterion was used to initiate contrail formation in global models, such as in the Contrail Cirrus Prediction Tool (CoCiP),⁹ and made it possible to assess the global environmental impact of contrails. To allow coverage of the full evolution cycle of a large population of contrails, the physics modeling was relatively simplified, such as for initialization of the species concentration with a Gaussian profile. Even though the results of regional air-traffic simulations^9,10 were quite consistent with in-flight measurements and satellite observations, the microphysical processes of ice production remain misunderstood. A more precise characterization of ice particles was achieved using 0-D microphysical trajectory box models.^11,12 This class of models allows detailed chemical kinetics and microphysics in the engine jet phase. For instance, simulations of young contrails^11,13 highlighted the domination of soot-induced heterogeneous freezing due to the large amount of emitted soot condensation nuclei from kerosene-fueled engines. In this respect, the reduction of soot emissions with the use of alternative jet fuels can help alter the optical properties of contrails.¹⁴ The box model results for plume dilution, however, showed some discrepancies compared to in-flight measurements due to (1) the absence of a detailed turbulence modeling^11,14 and (2) the simplified level of dynamics using averaged bulk exhaust (provided by semiempirical models).^11,14 Some authors remedied the shortcomings through offline-coupling strategies with an external fluid-flow solver.^15,16

Figure 1.

Evolution of a contrail in the wake of an aircraft.

In the last two decades, computational cost-effectiveness has increased interest in 3D CFD-based models to refine the early plume dynamics and associated mixing process. Early plume properties are known to influence later contrail stages (optical depth, radiative forcing) as shown by the studies on the engine jet and vortex phase.^7,17 As such, they are usually used as inputs in mesoscale models.^18,19 Solid particles motion can be tracked in CFD using either an Eulerian²⁰ or a Lagrangian.^21–23 The former can capture the motion of higher particle concentrations²⁴ with either bulk parametrization^25,26 or bins²⁰ for microphysics modeling. The simplest representation of particle-size distribution considers a single bin (monodisperse distribution) as used by Khou et al.^27,28 who investigated the near-field contrails (engine jet phase) behind realistic commercial aircraft geometries using Reynolds-averaged Navier–Stokes (RANS) equations for turbulence modeling. Results revealed that wingtip vortices drained and cooled jets and highlighted the role of ambient relative humidity on ice-crystal properties. Using a single bin could, however, prevent accurate estimates of ice-crystal concentration.²⁸ The most physically accurate treatment of ice particles in Eulerian terms considers size-resolved (binned) microphysics,²⁰ which avoids the problems of numerical diffusion inherent in bulk models^29,30 that could lead to slight discrepancies in the dispersion of ice and water-vapor fields. The computational cost, however, is excessive (see the detailed review³¹).

In contrast, Lagrangian tracking allows for an explicit process-oriented modeling of ice microphysics,³⁰ which accounts for its popularity in contrail research for a while.^15,22,30 Most if not all Lagrangian studies are highly parameterized in a way that neglects nozzle exit-related aspects. For example, Paoli et al.³² modeled free turbulent jets with classical expansion laws to investigate contrail formation for two- and four-engine configurations. Unterstrasser and Görsch³³ also used Gaussian-like plume profiles to study the vortex phase and to characterize the impact of aircraft type on contrail evolution. These simplified representations of the nozzle exhaust jet do not allow for accurately taking into consideration the features of modern turbofan nozzles. As the flow structure and mixing in the early plume depend on engine exhaust parameters, the impact of high bypass rates and complex geometries (e.g. short and long cowls and chevrons) on the near-field microphysics remain unclear. This particularly motivated the development of a CFD-microphysical model that integrates the nozzle-exit geometry and avoids using parameterizations to describe the exhaust plume.

Therefore, the main objective of this study was to implement a Lagrangian microphysical module in a CFD code to simulate ice-particle growth in near-field contrails starting from a realistic turbofan nozzle exit. To the authors’ knowledge, no contrail study has integrated the actual nozzle-exit geometry in a Lagrangian approach. A high-resolution spatial unsteady RANS method was therefore used to demonstrate the feasibility of the present approach. The paper is organized as follows. Section 2 begins with a description of the modeling approach, followed by the mathematical equations of fluid dynamics, particle motion, and ice microphysics. Section 3 describes the model and presents how the microphysical module was coupled to Simcenter Star-CCM+ CFD software. Section 4 discusses the validation test results for jet dynamics, plume dilution, and ice-particle properties.

Model description

Approach for modeling the formation of near-field contrails

Figure 2 shows that the nozzle exit flow from a turbofan engine at cruise flight consists of two jets: a hot core flow and a surrounding cold bypass flow. Both jets are mixed with the freestream moist air. The bypass flow simply contains ambient moist air that has been slightly accelerated by the engine fan, while the core jet drives the fuel-combustion products composed of gaseous emissions laden with soot particles (dispersed or solid phase). Water vapor in the dilution plume condenses onto soot particles, which serve as condensation nuclei when activated with sulfuric acid (heterogeneous nucleation¹¹). Under saturation conditions, ice particles form and grow by taking up more condensable matter – that is, water vapor – from the plume and ambient air, leading to contrail formation.

Figure 2.

2D Scheme of exhaust flow from a turbofan engine with ice-particle formation.

In this study, the Eulerian–Lagrangian approach was used to model the formation of contrail ice particles in the near field of an aircraft engine. The continuous gas phase was modeled with the Eulerian approach, while the dispersed-phase of solid particles (soot and ice particles) were tracked with the Lagrangian approach. Lastly, a one-way microphysical coupling was considered to account for ice-particle growth due to water-vapor condensation.

Gas-phase equations

The compressible flow of the gas phase (“g” index) is considered to be a mixture of two miscible species: air (“a” index) and water vapor (“w” index). Assuming an ideal gas, the mixture density ${\bar{ρ}}_{g}$ , temperature ${\tilde{T}}_{g}$ , and pressure ${\bar{p}}_{g}$ are linked with the equation of state:

{\bar{p}}_{g} = {\bar{ρ}}_{g} r_{g} {\tilde{T}}_{g}

(1)

where $r_{g}$ is the specific gas constant. The overbar sign – denotes the Reynolds (ensemble) averaging, while ~ refers to the Favre (density-weighted) averaging defined as $\bar{ρ} \tilde{Φ} = \bar{ρ Φ}$ and $Φ = \tilde{Φ} + Φ ″$ . The latter is considered to account for compressibility effects. As the growth of ice particles is a time-dependent process, the flow is assumed to be unsteady. Therefore, the unsteady Favre-averaged Navier–Stokes (FANS) equations expressing the conservation of the gas-phase global mass, momentum, and energy were solved. In Cartesian coordinates ( $i$ , $j$ , $k$ ), these equations read as follow with some modeling assumptions.³⁴

\frac{\partial {\bar{ρ}}_{g}}{\partial t} + \frac{\partial}{\partial x_{i}} ({\bar{ρ}}_{g} {\tilde{u}}_{i, g}) = 0

(2)

\begin{matrix} \frac{\partial}{\partial t} ({\bar{ρ}}_{g} {\tilde{u}}_{i, g}) + \frac{\partial}{\partial x_{j}} ({\bar{ρ}}_{g} {\tilde{u}}_{i, g} {\tilde{u}}_{j, g}) \\ = \frac{\partial}{\partial x_{j}} (- {\bar{P}}_{g} δ_{ij} + {\tilde{τ}}_{ij}^{visc} + {\tilde{τ}}_{ij}^{turb}) \end{matrix}

(3)

\begin{matrix} \frac{\partial}{\partial t} ({\bar{ρ}}_{g} {\tilde{e}}_{g}^{tot}) + \frac{\partial}{\partial x_{j}} ({\bar{ρ}}_{g} {\tilde{u}}_{j, g} {\tilde{e}}_{g}^{tot} + {\tilde{u}}_{j, g} {\bar{p}}_{g}) \\ = \frac{\partial}{\partial x_{j}} ({\tilde{q}}_{j, g}^{visc} + {\tilde{q}}_{j, g}^{turb}) - \frac{\partial}{\partial x_{j}} ({\tilde{u}}_{i, g} ({\tilde{τ}}_{ij}^{visc} + {\tilde{τ}}_{ij}^{turb})) \end{matrix}

(4)

Solving equations (2)–(4) provides the primary variables of the mean flow ${\bar{ρ}}_{g}$ , the velocity field ${\tilde{u}}_{g} = ({\tilde{u}}_{g}, {\tilde{v}}_{g}, {\tilde{w}}_{g})$ , and the total energy ${\tilde{e}}_{g}^{tot}$ , such as $c_{p, g}$ denotes the specific heat $μ_{g}$ the dynamic viscosity, $μ_{t, g}$ the eddy viscosity, $P r_{g}$ the Prandtl number, and $P r_{t, g}$ the turbulent Reynolds number of the gas phase. The latter physical properties are all functions of ${\tilde{T}}_{g}$ , which was derived from the total energy using equation (5). The averaging of nonlinear terms on the right-hand side of equations (3) and (4) yields the Reynolds stress ${\tilde{τ}}_{ij}^{turb}$ and the turbulent heat flux ${\tilde{q}}_{j, g}^{turb}$ , in addition to their laminar (or viscous) counterparts ${\tilde{τ}}_{ij}^{visc}$ .and ${\tilde{q}}_{j, g}^{visc}$ , given by equations (6)–(9), respectively, for a Newtonian fluid.

{\tilde{e}}_{g}^{tot} = c_{v, g} {\tilde{T}}_{g} + \frac{{\tilde{u}}_{k, g} {\tilde{u}}_{k, g}}{2} + \frac{\tilde{u_{k, g}^{″} u_{k, g}^{″}}}{2}

(5)

{\tilde{τ}}_{ij}^{visc} = μ (\frac{\partial {\tilde{u}}_{i, g}}{\partial x_{j}} + \frac{\partial {\tilde{u}}_{j, g}}{\partial x_{i}} - \frac{2}{3} \frac{\partial {\tilde{u}}_{k, g}}{\partial x_{k}} δ_{ij})

(6)

\begin{matrix} τ_{ij}^{turb} = - \bar{ρ_{g} u_{i, g}^{″} u_{j, g}^{″}} = μ_{t, g} \\ (\frac{\partial {\tilde{u}}_{i, g}}{\partial x_{j}} + \frac{\partial {\tilde{u}}_{j, g}}{\partial x_{i}} - \frac{2}{3} \frac{\partial {\tilde{u}}_{k, g}}{\partial x_{k}} δ_{ij}) - \frac{1}{3} {\bar{ρ}}_{g} \tilde{u_{k, g}^{″} u_{k, g}^{″}} δ_{ij} \end{matrix}

(7)

{\tilde{q}}_{j, g}^{visc} = c_{p, g} \frac{μ_{g}}{P r_{g}} \frac{\partial {\tilde{T}}_{g}}{\partial x_{j}}

(8)

{\tilde{q}}_{j, g}^{turb} = c_{p, g} \frac{μ_{t, g}}{P r_{t, g}} \frac{\partial {\tilde{T}}_{g}}{\partial x_{j}}

(9)

To compute the mixture composition in the gas phase ( $Φ_{g} = {\tilde{Y}}_{a} Φ_{a} + {\tilde{Y}}_{w} Φ_{w}$ ), the mass fractions of air and water vapor ( ${\tilde{Y}}_{a}$ and ${\tilde{Y}}_{w}$ , respectively) were derived with equations (10) and (11) defined as follow:

\begin{matrix} \frac{\partial}{\partial t} ({\bar{ρ}}_{g} {\tilde{Y}}_{w}) + \frac{\partial}{\partial x_{i}} ({\bar{ρ}}_{g} {\tilde{Y}}_{w} {\tilde{u}}_{i, g}) \\ - \frac{\partial}{\partial x_{i}} ((\frac{μ_{g}}{σ_{g}} + \frac{μ_{t, g}}{σ_{t, g}}) \frac{\partial {\tilde{Y}}_{w}}{\partial x_{i}}) = - {\bar{ω}}_{ice} \end{matrix}

(10)

\tilde{Y_{a}} = 1 - {\tilde{Y}}_{w}

(11)

Equation (10) translates the mass conservation of water vapor where the source term ${\bar{ω}}_{ice}$ on the right-hand side is the rate of condensed matter representing the mass transfer due to ice-particle growth. This coupling term is presented in greater detail in the next section. The variables $σ_{g}$ and $σ_{t, g}$ denote the Schmidt number and the turbulent Schmidt number, respectively taken as being equal to 1³⁵ and 0.9.^27,28

Lastly, the eddy viscosity $μ_{t, g}$ term introduced by the Boussinesq approximation (equation (7)) was modeled using the Realizable $k - ε$ model³⁶ to close the set of equations. This model is known to provide improved predictions in terms of the dissipation rate for round jets.³⁷ It was chosen after comparing three eddy-viscosity models to the available experimental data in the case of compressible coaxial jets. Validation tests are presented below in section 5.1. The $k - ε$ Realizable formulation uses two transport equations: one for the turbulent kinetic energy k and one for the dissipation rate $ε$ , such as:

S_{ij} = \frac{1}{2} (\frac{\partial {\tilde{u}}_{i, g}}{\partial x_{j}} + \frac{\partial {\tilde{u}}_{j, g}}{\partial x_{i}})

(12)

μ_{t, g} = {\bar{ρ}}_{g} C_{μ} \frac{k^{2}}{ε}

(13)

\begin{matrix} \frac{\partial}{\partial t} ({\bar{ρ}}_{g} k) + \frac{\partial}{\partial x_{i}} ({\bar{ρ}}_{g} k {\tilde{u}}_{i, g}) = - \frac{\partial}{\partial x_{i}} (\frac{μ_{t, g}}{σ_{k}} \frac{\partial k}{\partial x_{i}}) \\ + 2 μ_{t, g} S_{ij} \cdot S_{ij} - {\bar{ρ}}_{g} ε \end{matrix}

(14)

\begin{matrix} \frac{\partial}{\partial t} ({\bar{ρ}}_{g} ε) + \frac{\partial}{\partial x_{i}} ({\bar{ρ}}_{g} ε {\tilde{u}}_{i, g}) = \frac{\partial}{\partial x_{i}} (\frac{μ_{t, g}}{σ_{ε}} \frac{\partial ε}{\partial x_{i}}) \\ + 2 C_{1 ε} μ_{t, g} \frac{ε}{k} S_{ij} \cdot S_{ij} - C_{2 ε} {\bar{ρ}}_{g} \frac{ε^{2}}{k + \sqrt{ν_{g} ε}} \end{matrix}

(15)

where $S_{ij} = \frac{1}{2} (\partial {\tilde{u}}_{i, g} / \partial x_{j} + \partial {\tilde{u}}_{j, g} / \partial x_{i})$ represent the deformation rate of the mean flow. $C_{μ}$ , $σ_{k}$ , $σ_{ε}$ , $C_{1 ε}$ , and $C_{2 ε}$ are the model constants (see Refs.^34,38).

Dispersed-phase equations and microphysics model

Particles motion

At cruise conditions, the number density of physical particles is so high that a modeling simplification is required to save simulation time. For example, a typical value for emitted soot particles from an engine nozzle is about 10¹¹ #· m⁻³.³⁹ Therefore, computational particles are considered as spheres of radius $r_{p}$ , with each representing N physical particles.⁴⁰ Numerically, soot and ice particles are treated as passive tracers since their size and relaxation time are relatively small compared to the characteristic scales of the gas-phase flow (for ice particles $τ_{p} ~ 10^{- 5}$ s with the largest $r_{p}$ of a few microns).⁴¹ Thus, they have the same velocity as the carrier phase, but the drag and gravity forces are neglected since they have neither mass nor volume. The equation governing particle movement reads as follows:

\frac{d x_{p}}{dt} = \tilde{u} (x_{p}, t)

(16)

where $x_{p}$ is the particle position vector.

Ice-particle growth

In contrails, ice-particle growth is due to condensation of water vapor onto suitable nucleation sites on pre-activate particles (mostly soot particles).⁸ The ratio between the local water vapor partial pressure ( $P_{w}$ ) and the saturated vapor pressure with respect to liquid water $(P_{w, liq}^{0})$ defines the saturation ratio⁵ noted S:

S = \frac{P_{w} ({\tilde{T}}_{p})}{P_{w, liq}^{0} ({\tilde{T}}_{p})}

(17)

such as the local pressure $P_{w}$ is calculated according to Dalton’s law:

P_{w} ({\tilde{T}}_{p}) = X_{w} P_{tot}

(18)

where $X_{w}$ is the mole fraction of water vapor and $P_{tot}$ is the absolute total pressure. When the saturation ratio is above 1, heterogeneous ice nucleation takes place and spherical-shaped ice particles form and grow in the dilution plume. In this study, the microphysics of ice-crystal growth was modeled, while the processes of soot activation⁴² and homogeneous ice nucleation⁴³ were not taken into account. Therefore, the transition from continuum to molecular regime is supposed as being instantaneous.⁴⁴

The mass change of a single-particle (computational) growth due to condensation/evaporation effects can be expressed by a diffusion law as described by the Fukuta and Walter⁴⁴ model, which describes the change of both particle mass and size as follows:

{\overset{\cdot}{m}}_{p} = \frac{d m_{p}}{dt} = \frac{4 π ρ_{p} r_{p} (S - A)}{C_{T} A + C_{ρ}}

(19)

{\overset{\cdot}{r}}_{p} = \frac{d r_{p}}{dt} = \frac{S - A}{r_{p} (C_{T} A + C_{ρ})}

(20)

where $m_{p}$ denotes the mass of a particle p, $ρ_{p}$ the particle density, and A the Kelvin effect (defined in Appendix A). The factors $C_{T}$ and $C_{ρ}$ control the water-vapor temperature and air density just above the particle surface. The detailed formulations of both factors can be found in the Appendix A. Therefore, we obtain the mass variation of each ice crystal and a fortiori the vapor mass variation. This leads to the coupling term ${\bar{ω}}_{ice}$ in equation (10), which can be derived as follows:

{\bar{ω}}_{ice} = n \sum_{p = 1}^{N_{com}} {\overset{\cdot}{m}}_{p}

(21)

where n is the number density of particles contained in a volume cell V, such that $n = N_{p} / V$ . The total number of computational particles is $N_{com}$ , each of them representing $N_{p}$ physical particles. More details about the model are available in Fukuta and Walter.⁴⁴ The latter model was implemented in a microphysical module and coupled to the CFD commercial code in Simcenter STAR-CCM+ v13.04. A detailed description of the numerical configuration is provided below.

Numerical methodology

The CFD solver is based on a cell-centered finite-volume approach for unstructured grids.²⁰ The SIMPLE algorithm⁴⁵ was used to solve the flow governing equations. A second-order upwind scheme was used for space discretization, while time integration was achieved with a second-order Euler scheme. For the model of ice-particle growth, a fourth-order Runge-Kutta scheme was used to discretize the equation (20) of particle size change. The discretization of the equation (20) is presented in Appendix B. The microphysical module was implemented as a user-defined function in the commercial CFD code in Simcenter STAR-CCM+ v13.04. The script of the microphysical module (coded in C) can be downloaded from Ref.⁴⁶ The methodology for coupling with the CFD solver is detailed in Figure 3. At each time step, the primary flow variables of both gas- and dispersed-phase are derived from the CFD solver. The pressure and temperature at each grid point help to compute the thermodynamic variables ( $P_{w}$ , $P_{w, liq}^{0}$ , $C_{T}$ , and $C_{ρ}$ ), allowing derivation of the local saturation ratio. The residence time (τ_p) of each particle is computed at each time step. Radius of new particles (τ_p = 0) is initialized with a monodisperse distribution using an initial radius of $r_{p, 0} = 20 nm$ , as per numerical studies.^22,27 The number density of soot particles was set to $n = 10^{11}$ #·m⁻³, which corresponds to an emission index $E I_{soot}$ , of $4 \times 10^{14}$ particles per kg of fuel according to in situ measurements.³⁹ The radius of old particles (τ_p > 0) was recovered from a particle-radius table. For the following time steps, the evolution of particle radius depends on the local saturation-ratio values. Nonsaturated particles keep the same radius, while saturated particles grow. The rates of change for their mass ( ${\overset{\cdot}{m}}_{p}$ ) and size ( ${\overset{\cdot}{r}}_{p}$ ) were computed with equations (19) and (20). As such, the particle-radius table was updated and the mass-transfer term ${\bar{ω}}_{ice}$ derived with equation (21). The coupling term was injected into the CFD solver to correct the water-mass fraction in the following time step.

Figure 3.

Schematic of the CFD-microphysics coupling for modeling ice-particle growth.

Validation tests and results

Two validation tests were performed to assess the reliability of the CFD-microphysics model presented above for studying ice-particle growth in the exhaust jet of an aircraft engine. In the first test, the model’s gas-phase resolution (i.e. jet dynamics and turbulence) was assessed with experimental results from the canonical case of compressible coaxial jets. This validation test made it possible to choose the most suitable turbulence model. In the second test, the model’s microphysical resolution in terms of dilution and ice-particle properties (size, distribution, optical depth) were assessed in the exhaust jet of a CFM56-3 engine at realistic cruise flight conditions.

Compressible coaxial turbulent jets

Computational domain and grid

Guitton et al.⁴⁷ experimentally characterized the flow dynamics and turbulence in the case of isothermal compressible coaxial jets of air. This experiment was used as a validation test case for our computational model. As shown in Figure 4, the nozzle had a length $L_{noz} = 0.158 m$ and was composed of two jets: an inner primary jet, whose diameter was $D_{p} = 0.055 m$ , and an outer secondary jet with a diameter of $D_{s} = 0.1 m$ . The mean flow was assumed to be axisymmetric. Indeed, it has been observed experimentally that the vortex structures forming from the rollup of the shear–layers remain axisymmetric over a distance of a few diameters.⁴⁸ Hence, the computational domain was the two-dimensional rectangular box measuring $40 D_{s} \times 15 D_{s}$ based on the half-nozzle model, and the Reynolds-averaged Navier–Stokes (RANS) equations were solved in cylindrical coordinates $(O, \vec{r}, \vec{y})$ . Two velocity-inlet boundaries were used for the engine primary- and secondary-jet flows, with corresponding values of $U_{p} = 170 m . s^{- 1}$ and $U_{s} = 120 m . s^{- 1}$ , respectively (Figure 4). A symmetry condition was applied to the centerline (bottom boundary), while the remaining boundaries were treated as ambient-pressure outlet. A no-slip boundary condition was used on nozzle walls (5 mm thick).

Figure 4.

Scheme of the computational domain with nozzle geometry and boundary conditions.

A grid convergence study was carried out using three grids: coarse (290,110 cells), medium (618,651 cells), and fine (1,227,567 cells), respectively labeled grid 1, 2, and 3. The GCI method⁴⁹ was used to compute the discretization errors. Results of grid errors based on the axial velocity at a point probe located in the mixing layer between the two exhaust jets ( $r = 0.5 D_{p}, y = 10 D_{s}$ ) were sufficiently low: ${GCI}_{coarse}^{12} = 1.45 %$ and ${GCI}_{fine}^{23} = 0.94 %$ . As such, the medium grid was selected for computations. This mesh configuration had a total of 35 prismatic layers next to all wall boundaries, which ensured a $y^{+} \leq 1$ .

Jet dynamics and turbulence statistics

The mixing process and associated turbulence effects have a significant impact on the dilution and ice-particle microphysics in an exhaust jet.³² As such, the model’s gas-phase resolution was first validated by studying the performance of three eddy viscosity turbulence models: $k - ε$ Standard,³⁶ $k - ε$ Realizable,⁵⁰ and $k - ω$ SST.⁵¹ This comparative study covered first- and second-order statistics, that is, mean axial velocity, turbulent kinetic energy (TKE), and the Reynolds stress.

Figure 5 presents the cross-sectional profiles of the exhaust flow against experimental data.³¹ Results of axial velocity profiles (Figure 5(a)) confirm that the predictions of the three tested turbulence models are satisfactory overall. The Standard and SST versions slightly overestimated the axial velocity on the external part of the jet (i.e. $r \geq 0.5 D_{s}$ ). Small discrepancies can also be observed near the primary-jet centerline (i.e. $y \geq 7 D_{s}$ ) for the $k - ε$ Standard, while the Realizable version seems to offer the best fit for experimental measurements.³² The Reynolds-stress profiles (Figure 5(b)) also confirm the latter conclusion in which the Standard and SST models overestimated the experimental profiles for $y \geq 5 D_{s}$ .

Figure 5.

Comparison of exhaust-flow cross-sectional profiles for the three turbulence models tested against experimental results³¹: (a) axial velocity and (b) Reynolds stress.

To better assess the discrepancies observed near the jet axis in Figure 5(a), the performance study on turbulence models was carried out along the primary-jet centerline (see the results in Figure 6). The three models tested yielded similar performance in the potential core region (Figure 6(a)). In contrast, both velocity decay and TKE increase (Figure 6(b)) were best captured by the $k - ε$ Realizable model. Quantitatively, the $k - ε$ Realizable provided the lowest mean relative error with respect to experimental data on the mean axial velocity: 2.2% compared to 3.8% and 4.9% for $k - ω$ SST and $k - ε$ Standard, respectively. Based on this comparative study, the $k - ε$ Realizable was selected.

Figure 6.

Comparison of primary-jet flow results for the three turbulence models tested against experimental results³¹: (a) axial velocity and (b) turbulent kinetic energy.

Contrail formation in the CFM56-3 exhaust jet

Computational domain and grid

Ice-particle formation was modeled in the exhaust plume of the CFM56-3 engine at realistic cruise flight conditions that ensure contrail formation.²⁷ The engine speed was set at $U_{cruise} = 252$ m·s⁻¹ and ambient conditions were taken at $35, 000 ft$ . The static temperature and pressure values were $T_{amb} = 219 K$ and $P_{amb} = 23, 800$ Pa, respectively. The mass fraction of water vapor in ambient air was set at $Y_{amb} = 6.08 \cdot 10^{- 5}$ to account for 60% relative humidity for comparison with numerical results.²⁷

As the mean flow was assumed to be axisymmetric, only the 90° sector based on the engine nozzle was considered for URANS calculations, as shown in Figure 7. The geometry of the CFM56-3 short-cowl nozzle is composed of a core and a bypass flow whose diameters are $D_{c} = 0.646 m$ and $D_{b} = 1.586 m$ . The nozzle configuration is representative of modern commercial turbofan engines,⁵² in which a central plug is included to encourage early mixing between the core and bypass flows (hereafter the inner shear layer). The exhaust conditions used by Garnier et al.⁵³ were used in this study (see Table 1). According to Paoli et al.³² physical soot particles in the range 10⁶–10⁸ preserves the statistical reliability of the simulation results. As such, 10⁷ physical particles were used herein. This ensures 40,000 computational particles, and a total number of emitted soot particles at $n = 10^{11}$ #·m⁻³.

Figure 7.

Computational domain and boundary conditions for the CFM56-3, short-cowl nozzle: (a) isometric view, (b) 2D view, and (c) nozzle geometry.

Table 1.

Exhaust conditions of the CFM56-3 engine.⁵³

Core flow	Mach, $M_{c}$	1
	Velocity, $U_{c}$ (m·s⁻¹)	480
	Static temperature, $T_{c}$ (K)	580
	Number density of soot particles, n (#·m⁻³)	10¹¹
	Number of computational particles, $N_{com}$ (#)	4 10⁴
	Initial soot-particle diameter, $r_{p, 0}$ (nm)	20
	Soot-particle distribution	Monodisperse
	Mass fraction of water vapor, $Y_{c}$	3.5 10⁻²
Bypass flow	Mach, $M_{b}$	1
	Velocity, $U_{b}$ (m·s⁻¹)	311
	Static temperature, $T_{b}$ (K)	242

As shown in Figure 7, the computational domain had a radius $R = 5 D_{b}$ and a length $L = 76 D_{b}$ to ensure ambient pressure at far-field boundaries. Three inlet sections were used upstream: a freestream pressure inlet ( $P_{amb}^{0}$ , $T_{amb}^{0}$ ), a bypass pressure inlet ( $P_{b}$ , $T_{b}$ ), and a core pressure inlet ( $P_{c}$ , $T_{c}$ ). A no-slip adiabatic condition was imposed on engine walls, while a symmetry condition was used on both lateral planes A turbulence intensity of 10% was used for the exhaust jets, as per Garnier et al.,⁵³ while a value of 0.1% was applied on both freestream and far-field boundaries ( $P_{amb}$ , $T_{amb}$ ). Remaining boundaries were treated with a far-field condition.

The geometry and mesh of the computational domain were generated with ANSYS ICEM CFD software.⁵⁴ A grid convergence study was performed on the gaseous phase using three grids: coarse (4.35M cells), medium (9.8M cells), and fine (20.9M cells), respectively labeled grid 1, 2, and 3. The mean errors computed in the inner shear layer (i.e. between the core and the bypass jets) for the static temperature, water-vapor mass fraction, and core axial velocity are below 1% (see Table 2). Furthermore, discretization errors estimated using the GCI method⁴⁹ at a probe point ( $r = 0.7 D_{c}, y = 10 D_{b}$ ) were sufficiently low with the medium mesh. As such, computations were only performed on the medium grid.

Table 2.

Grid study results: mean relative error $ε$ , standard deviation $σ$ , and GCI index⁵³ between grids 1 (coarse), 2 (medium), and 3 (fine) in the inner shear layer ( $r = 0.7 D_{c}$ ).

Flow variable	Static temperature ${\tilde{T}}_{g}$	(H₂O)_v Mass fraction ${\tilde{Y}}_{w}$	Axial velocity ${\tilde{w}}_{g}$
$ε_{coarse}^{12}$ ( $σ_{coarse}^{12}$ ), %	0.2 (10⁻³)	0.6 (3·10⁻³)	0.4 (1.9·10⁻³)
${GCI}_{coarse}^{12}$ , %	0.13	1.47	0.69
$ε_{fine}^{23}$ ( $σ_{fine}^{23}$ ), %	3.8 10⁻³ (4.6·10⁻⁴)	0.4 (3·10⁻³)	0.2 (1.4·10⁻³)
${GCI}_{fine}^{23}$ , %	0.03	1.01	0.22

The medium grid is mainly composed of hexahedral elements, as shown in Figure 8. Three cell numbers $n_{r}$ , $n_{θ}$ , and $n_{z}$ characterize the global mesh in cylindrical coordinates ( $r, θ, z$ ). Refinement blocks were used in the exhaust jet and shear layer, where gradients are expected to be relatively high. Table 3 presents the cell numbers on each boundary condition for the medium mesh. The minimum cell size in the inner shear layer was Δy = 10⁻³ m, Δx = Δz = 5·10⁻⁵ m, compared to Δy = 10⁻³ m, Δz = Δx = 10⁻² m in the outer shear layer (between the freestream and bypass flow). The near-wall refinement ensured a maximum y⁺ of unity, which makes it possible to resolve the boundary layer. The mesh in the far field is relatively coarser with cells about Δy = 0.2 m, Δx = Δz = 1.09 m. A time step Δt = 10⁻³ s was chosen ensuring a maximum Courant–Friedrichs–Lewy (CFL) number below 1 in the entire domain except in the potential core region where the flow speed is equal to core exit velocity $U_{c}$ = 480 m·s⁻¹ ( $M_{c}$ = 1) and the CFL reached a value of 10.

Figure 8.

Computational grid for the medium mesh: (a) axial view, (b) lateral view, (c) close-up view of the nozzle exit, and (d) close-up view of the central plug.

Table 3.

Mesh specifications on boundary conditions for the medium configuration.

Boundary conditions	Number of cells
Boundary conditions	$n_{r}$	$n_{θ}$	$n_{z}$
Freestream inlet	130	20	–
Core inlet	70	20	–
Bypass inlet	60	20	–
Freestream-bypass wall	–	20	20
Bypass-core wall	–	20	50
Short-cowl wall	–	20	25
Symmetry planes	280	–	1895
Top plane	–	20	1895
Far-field outlet	280	20	–

Jet flow and plume dilution

Figure 9(a) and (b) present the stream-wise contours of the axial mean velocity and Mach number, respectively, for the CFM56-3 plume. Since both dual-stream jets were at sonic conditions (M = 1), weak shock waves were observed near the exhaust sections and edges (Figure 9(b)); similar behavior was reported at nozzle-exit flows.^55,56 A flow separation on the central plug can be seen in the core flow due to the over-expansion shock. The velocity contours in Figure 9(a) show that the potential core length was about $5 D_{b}$ ; this distance is comparable to both experimental (5.4D_b)⁵² and numerical (5D_b)⁵⁷ values.

Figure 9.

Stream-wise (xy-plane) contours of: (a) axial mean velocity and (b) Mach number.

The dilution ratio defined by Schumann et al.⁵⁸ was used to characterize the evolution of flow variables (such as temperature and contrail diameter) along the exhaust plume. The dilution ratio based on a passive scalar Z was calculated as

N_{Z} = \frac{N_{exit}}{\tilde{Z}}

(22)

where $N_{exit}$ is the dilution factor at the exhaust section⁵⁸ and $\tilde{Z}$ is the ratio of the amount of fluid from the jet region to that from the freestream.⁵³ The latter was calculated from

\tilde{Z} = \frac{{\tilde{Z}}_{y, plume} - {\tilde{Z}}_{i, atm}}{{\tilde{Z}}_{i, jet} - {\tilde{Z}}_{i, atm}}

(23)

where ${\tilde{Z}}_{y, plume}$ is the passive scalar concentration in the exhaust plume, ${\tilde{Z}}_{i, atm}$ is the initial atmospheric concentration, and ${\tilde{Z}}_{i, jet}$ denotes the initial jet concentration. The concentration measurement was taken along the plume centerline.

Figure 10 shows the evolution of the dilution ratio $N_{Z}$ as function of the plume age compared to the bulk mean dilution interpolated from in-flight measurements.⁵⁹ For instance, the interpolation was based on a large set of data from more than 70 plume encounters with the DLR Falcon research aircraft and the NASA ER-2 aircraft.⁵⁹ It should be noted that the bulk dilution ratio corresponded to the radial average of the local dilutions ratios within the plume at a given distance from the nozzle exit. This explains why the predicted dilution ratio along the core jet as seen in Figure 10 was relatively smaller than the bulk interpolation of the mixed primary and secondary jets (as also highlighted by Schumann et al.⁵⁹), especially near the engine exhaust. Even with that, the variations of the dilution ratio in the simulated plume are of the same order of magnitude as the bulk mean interpolation. The study of ice-crystal properties during contrail formation is presented in the next section.

Figure 10.

Evolution of the predicted dilution ratio along the primary jet centerline against in-flight measurements.⁵⁶

Contrail formation and ice-particle properties

Figure 11(a) presents the spatial distribution of ice-particle radii in the plume exhaust. Up to a distance of 10 $D_{b}$ downstream, the jet temperature was high enough to prevent ice coating of soot particles with an initial radius of $r_{p, 0} = 20 nm$ . The mixing with ambient air downstream cooled the plume, and the first ice crystals can be observed within a thin layer of the plume edge at $y = 10 D_{b}$ . More ice particles formed and grew subsequent to the thickening of the external shear layer by transverse diffusion up to $y ~ 50 D_{b}$ where water-vapor condensation reached the plume core.

Figure 11.

Spatial distribution of particle radii: (a) longitudinal cross-sectional distribution, transverse cross-sectional distributions at (b) $y = 30 D_{b}$ and (c) $y = 72 D_{b}$ .

The rate of condensation along the exhaust plume was fully determined by the local ice saturation ratio S, which, in turn, depends on the local temperature ${\tilde{T}}_{g}$ and water mass fraction ${\tilde{Y}}_{w}$ .³² Figure 11(b) and (c) compare two transverse cross-sectional distributions at $y = 30 D_{b}$ and 72 $D_{b}$ from the core exhaust section, respectively. Both engine core and bypass sections are represented by black circles (solid line). At $y = 30 D_{b}$ , the particles were spatially distributed over a circular area of diameter $2 D_{b}$ . Ice condensation was only observed on the external shear layer ( $r \geq 0.5 D_{b}$ ) while soot particles transported along the plume core remain under-saturated. The transverse dispersion of ice crystals enhanced by mixing with ambient air downstream resulted in a radial expansion of the plume up to a diameter 3 $D_{b}$ at $y = 72 D_{b}$ . At this section, the center of the plume was fully saturated and, due to the relatively high water content in the core jet, ice particles seemed to grow faster along the plume core ( $r \leq D_{b}$ ), up to a radius about 5 $0 r_{p, 0}$ , as compared to the plume edge. At the outer edge of the plume; ice particle radii reached a constant value around 500 nm for $y / D_{b}$ > 50. Indeed, the growth of ice-particles (i.e. condensation rate) is mainly governed by the local value of the saturation ratio S, which is below 1 in the outer edge of the plume. As a consequence, the ice-particles do not grow in the outer edge and reach a constant size of 500 nm.

The evolution of the saturation ratio along the exhaust jet shows that the plume was fully saturated, that is, $S = 1$ , at a distance about $55 D_{b}$ behind the nozzle exit (see Figure 12). This is consistent with the numerical results of Khou⁶⁰ obtained behind a full aircraft geometry with an Eulerian approach for particle tracking. Since the growth rate of ice crystals is directly related to the local saturation-ratio value, the mean particle radius follows the same trend as the bulk mean saturation ratio. It should be noted that that the mean particle radius $< r_{p} (y) >$ corresponds to the arithmetic average of particle radii over the subdomain of length h centered on the position y, defined as³²

< r_{p} (y) > = \frac{1}{h} \int_{y - h / 2}^{y + h / 2} r_{p} d y^{'} = \frac{1}{N} \sum_{i = 1}^{N} r_{p, i}

(24)

where N is the particle number within the considered subdomain and $r_{p, i}$ is the radius of a given particle i. Overall, the results show that our model captures well the global dynamic of ice-particle growth compared to the numerical prediction by Khou.⁶⁰ The mean radius was initially constant near the engine exhaust and started growing around $y = 10 D_{b}$ . At $0.5 s$ of the plume age ( $y = 72 D_{b}$ ), the mean particle radius was in the good range of 500–1000 nm, as confirmed by several numerical studies.^15,16,61

Figure 12.

Evolutions of the mean particle radius and the saturation ratio along the jet plume.

To better assess the particle-size distribution along the exhaust plume, the probability density functions (PDF) were investigated at four distances: $20 D_{b}$ , $30 D_{b}$ , $50 D_{b}$ , and 72 $D_{b}$ behind the engine exhaust (see Figure 13). Results show that the PDF distributions shifted toward large radii following the growth of ice crystals along the exhaust plume. Moreover, Gaussian distributions were obtained from $y = 50 D_{b}$ indicating polydispersed distributions of ice-particle radii, as observed in other numerical studies.^22,32

Figure 13.

Probability density functions (PDF) of particle radius at four axial distances: $20 D_{b}$ , $30 D_{b}$ , $50 D_{b}$ , and 72 $D_{b}$ .

Lastly, the visibility in the predicted contrail was assessed with the optical depth $τ$ . This parameter measures the attenuation of radiation passing through clouds and is often reported at a wavelength in the visible part of the radiation spectrum.⁷ The optical depth $τ$ was evaluated as formulated by Naiman et al.⁶²

τ = \sum_{i = 1}^{N_{com, g}} π r_{p}^{2} Q_{ext} \frac{N_{g}}{Δ V}

(25)

where $N_{com, g}$ denotes the number of computational particles in the grid cell, $N_{g}$ the total number of hysical particles in the grid cell, and $Δ V$ its volume. Equation (26) approximates the Mie extinction efficiency $Q_{ext}$ at a given wavelength $λ_{w}$ , such as the refractive index for ice of $μ_{r} = 1.31$ and the wavelength of visible light is considered as $λ_{w} = 0.55 μ m$ .

Q_{ext} = 2 - \frac{4}{e} [\sin (e) - \frac{1 - \cos (e)}{e}]

(26)

e = \frac{4 π r_{p} (μ_{r} - 1)}{λ_{w}}

(27)

The visibility threshold of jet contrails is defined by the criterion $τ \geq 0.5$ .⁶¹ Figure 14 shows that the predicted contrail from the CFM56-3 engine was visible at a distance about $26 m$ from the engine exhaust. The latter value is good agreement with the in situ observations reporting visible contrails at a distance in the range 10–30 m.^59,63

Figure 14.

Visible contrail predicted in the near-field of the CFM56-3 engine using the visibility criterion $τ \geq 0.5$ .

The comprehensive analysis performed shows that this model’s predictions are consistent with available numerical data. Such model can be used in future investigations to characterize the impact of complex nozzle geometries – such as chevron nozzles – on the optical and microphysical properties of near-field contrails. Research perspectives to be explored next involve the interactions of soot and sulfur species in the early plume. This requires the implementation of a gas-phase chemistry with the microphysics of volatile particles.

Conclusion

This paper demonstrated the feasibility of Eulerian–Lagrangian CFD coupling with ice microphysics to predict the near-field contrails starting from the nozzle-exit geometry of the CFM56-3 engine. The microphysical module implemented accounts for the main process of ice-particle growth based on heterogeneous nucleation. Two validation tests were performed to assess the reliability of the proposed CFD-microphysics model. In the first validation, the canonical case of compressible coaxial jets made it possible to choose the Realizable $k - ε$ turbulence model based on its good performance in terms of the jet dynamics and turbulence properties. In the second test case, the model’s microphysical resolution was assessed in the exhaust jet from the CFM56-3 engine at realistic cruise flight conditions. The predicted plume dilution showed good agreement with available experimental data. Qualitatively, the spatial distribution of ice particles revealed that the first ice crystals were formed at the plume edge at a distance of about $10 D_{b}$ from the engine exhaust. Further downstream, ice particles formed and grew subsequent to the thickening of the external shear layer due to transverse diffusion. Results of the saturation ratio show that the plume core axis was saturated at a distance of about $55 D_{b}$ . Polydispersed distributions of particle radii were obtained in the simulated contrail with the largest ice particle reaching a radius about 1 $μ$ m in the plume core. As the jet evolved downstream, the particle-size distributions shifted toward larger particles. Lastly, the optical-depth analysis shows that the predicted contrail was visible around 26 m behind the engine exhaust, which is consistent with contrail observations.

The model presented herein is intended to serve as the foundation for more advanced numerical simulations to investigate the impact of complex nozzle-exit geometries, such as chevron nozzles. The implementation of gas-phase chemistry and volatile-particle microphysics are planned to extend the model capabilities.

Footnotes

Appendix A: Equations for the microphysical module

This appendix begins with a description of the factors $C_{T}$ and $C_{ρ}$ in the microphysical model.⁴⁴ This is followed by a description of the thermodynamics equations to calculate the saturation pressure, density, surface tension, and latent heat of water.

The factor $C_{T}$ was calculated as:

(28)

C_{T} = \frac{L_{w}^{2} M_{w} ρ_{p}}{k_{g} R {\tilde{T}}_{p}^{2} f_{3 α}}

where $L_{w}$ is the heat latent of water vapor (J·kg⁻¹), $M_{w}$ is the molar fraction of water (M_w = 18.015 10⁻³ kg·mol⁻¹⁾, $ρ_{p}$ is the particle density of water vapor (kg·m⁻³), $k_{g}$ is the conductivity of the gas (W·m⁻¹·K⁻¹), R is the universal gas constant (R = 8.314 J·mol⁻¹·K⁻¹), ${\tilde{T}}_{p}$ the particle temperature (K), and $f_{3 α}$ a correction factor for the temperature difference between the particle’s boundary layer and ambient air.

The conductivity of the gas was calculated from⁶⁴:

(29)

k_{g} = (5.69 \times 10^{- 5} + ({\tilde{T}}_{g} - 273.15) 0, 017 \times 10^{- 5}) 4.18 \times 10^{2}

The factor $f_{3 α}$ was calculated from:

(30)

f_{3 α} = \frac{r_{p}}{r_{p} + \frac{k_{g} \sqrt{2 π M_{a} R {\tilde{T}}_{p}}}{α p_{\infty} (C_{p, a} - R / 2)}}

where $M_{a}$ is the molar fraction of air (M_a = 28.964 10⁻³ kg·mol⁻¹), $α$ is the accommodation coefficient with a consensus value of 0.68,⁶⁵ $p_{\infty}$ is the atmospheric pressure (Pa), and $C_{p, a}$ is the molar heat of air at constant volume (J·K⁻¹·mol⁻¹).

The factor $C_{ρ}$ was calculated as:

(31)

C_{ρ} = \frac{R {\tilde{T}}_{p} ρ_{p}}{P_{w}^{0} D_{w} M_{w} f_{3 β}}

where $P_{w}^{0}$ is the saturation pressure (Pa), $D_{w}$ is the diffusivity of water vapor (m²·s⁻¹), and $f_{3 β}$ is a correction factor for the vapor-density difference between the particle’s boundary layer and the ambient air.

The factor A in equations (19) and (20) was calculated as:

(32)

A = \exp (\frac{2 σ_{p} M_{w}}{r_{p} ρ_{p} r_{g} {\tilde{T}}_{g}})

where $σ_{p}$ denotes the surface tension of water.

The diffusivity of water vapor was calculated from⁶⁴:

(33)

D_{w} = (0.211 {(\frac{{\tilde{T}}_{g}}{273.15})}^{1, 94} \frac{1.013 \times 10^{5}}{p_{g}}) 10^{- 4}

The factor $f_{3 β}$ was derived from:

(34)

f_{3 β} = \frac{r_{p}}{r_{p} + \frac{D_{w}}{β} \sqrt{\frac{2 π M_{w}}{R T_{p}}}}

where $β$ is the condensation coefficient taken as 0.022.⁶⁵

The thermodynamic equation of the model is described. These equations depend on the particle temperature. The model considered the state of water. When the water was liquid ( ${\tilde{T}}_{g}$ (x_p) ≥+5°C), equations (35)–(38) were used. When the water was ice ( ${\tilde{T}}_{g}$ (x_p) < +5°C), equations (39)–(42) were used.

The density of liquid water (kg·m⁻³) was calculated from⁶⁴:

(35)

\begin{matrix} ρ_{liq} ({\tilde{T}}_{g}) = \\ (0.9998 + 0.86 \times 10^{- 4} ({\tilde{T}}_{g} - 273.15) - 0.108 \times 10^{- 4} {({\tilde{T}}_{g} - 273.15)}^{2}) 10^{3} \end{matrix}

The surface tension of liquid water (J·m⁻²) was calculated from⁶⁶:

(36)

σ_{liq} ({\tilde{T}}_{g}) = (76.1 - 0, 155 ({\tilde{T}}_{g} - 273.15)) 10^{- 3}

The condensation latent heat of liquid water (J·kg⁻¹) was calculated from⁶⁴:

(37)

L_{liq} ({\tilde{T}}_{g}) = (597.3 {(\frac{273.15}{{\tilde{T}}_{g}})}^{0.167 + 3.61 \times 10^{- 4} {\tilde{T}}_{g}} 4.18 \times 10^{3})

The saturation pressure of liquid water (Pa) was calculated from⁶⁷:

(38)

\begin{matrix} P_{w, liq}^{0} ({\tilde{T}}_{p}) = 100 \times 10 (10.79574 (1 - \frac{273.16}{{\tilde{T}}_{p}}) - 5.028 \times lo g_{10} (\frac{{\tilde{T}}_{p}^{^}}{273.16}) + 1.50475 \times 10^{- 4} (1 - 10^{(- 8.2969 (\frac{{\tilde{T}}_{p}}{273.16} - 1))})) \\ + 0.42873 \times 10^{- 3} (10^{(- 4.76955 (1 - \frac{{\tilde{T}}_{p}}{273.16}))} - 1) + 0.78614 \end{matrix}

The density of ice was (kg·m⁻³) calculated from⁶⁴:

(39)

\begin{matrix} ρ_{ice} ({\tilde{T}}_{g}) = \\ (0.91676 - 1.75 \times 10^{- 4} ({\tilde{T}}_{g} - 273.15) - 5 \times 10^{- 7} {({\tilde{T}}_{g} - 273.15)}^{2}) 10^{3} \end{matrix}

The surface tension of ice (J·m⁻²) was calculated from⁶⁶:

(40)

σ_{ice} ({\tilde{T}}_{g}) = (104.6 - 0.095 ({\tilde{T}}_{g} - 273.15)) 10^{- 3}

The sublimation latent heat of liquid water (J·kg⁻¹) was calculated from⁶⁴:

(41)

\begin{matrix} L_{ice} ({\tilde{T}}_{g}) = ((597.3 {(\frac{273.15}{{\tilde{T}}_{g}})}^{0.167 + 3.61 \times 10^{- 4} {\tilde{T}}_{g}})) \\ + (79.7 + 0.485 ({\tilde{T}}_{g} - 273.15) \\ - 2.5 10^{- 3} {({\tilde{T}}_{g} - 273.15)}^{2}) 4.18 \times 10^{3} \end{matrix}

The saturation pressure of ice (Pa) was calculated from:

(42)

P_{w, ice}^{0} ({\tilde{T}}_{p}) = 100 \times 10^{(\frac{9.5 ({\tilde{T}}_{p} - 273.15)}{({\tilde{T}}_{p} - 273.15) + 265.5} + 0.7858)}

Appendix B: Discretization of the ice particle growth equation

This appendix describe the discretization of ice particle growth equation (20) with a fourth order Runge-Kutta method. The four constants $k_{1}$ , $k_{2}$ , $k_{3}$ , $k_{4}$ of the Runge-Kutta method are presented below:

(43)

{\begin{matrix} k_{1} = \frac{Δ t}{r_{p} (t_{n - 1})} [\frac{S - \exp (\frac{A^{'}}{r_{p} (t_{n - 1})})}{C_{T} \exp (\frac{A^{'}}{r_{p} (t_{n - 1})}) + C_{ρ}}] \\ k_{2} = \frac{Δ t}{[r_{p} (t_{n - 1}) + \frac{k_{1}}{2}]} [\frac{S - \exp (\frac{A'}{r_{p} (t_{n - 1}) + \frac{k_{1}}{2}})}{C_{T} \exp (\frac{A'}{r_{p} (t_{n - 1}) + \frac{k_{1}}{2}}) + C_{ρ}}] \\ k_{3} = \frac{Δ t}{[r_{p} (t_{n - 1}) + \frac{k_{2}}{2}]} [\frac{S - \exp (\frac{A'}{r_{p} (t_{n - 1}) + \frac{k_{2}}{2}})}{C_{T} \exp (\frac{A'}{r_{p} (t_{n - 1}) + \frac{k_{2}}{2}}) + C_{ρ}}] \\ k_{4} = \frac{Δ t}{[r_{p} (t_{n - 1}) + k_{3}]} [\frac{S - \exp (\frac{A'}{r_{p} (t_{n - 1}) + k_{3}})}{C_{T} \exp (\frac{A'}{r_{p} (t_{n - 1}) + k_{3}}) + C_{ρ}}] \end{matrix}

The factor $A'$ is computed as follows:

(44)

A' = \frac{2 σ M_{w}}{ρ r_{g} {\tilde{T}}_{g}}

where $σ$ denotes the surface tension of water and $ρ$ the particle density.

The new particle radius $r_{p} (t_{n})$ is computed as follows:

(45)

r_{p} (t_{n}) = r_{p} (t_{n - 1}) + (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}) / 6

Then, the condensed matter of ice particles associated to equation (19) is calculated as follows

(46)

m_{p} (t_{n}) = \frac{ρ_{ice} 4 π}{Δ t} (r_{p} (t_{n}) - r_{p} (t_{n - 1})) {(r_{p} (t_{n}))}^{2} n

where n is the number density of particles contained in a volume cell V, such that $n = N_{p} / V$ .

Acknowledgements

The authors thank Charles E. Tinney for providing data measurements on coaxial jets conducted at the Laboratoire d’Etudes Aérodynamiques (LEA) in Poitiers, France. Computations were made on the Beluga supercomputer at the École de technologie supérieure (ÉTS) in Montréal, managed by Calcul Quebec and Compute Canada.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), the Ministère de l’Économie, de la science et de l’innovation du Québec (MESI), 531, and the Fonds de recherche du Québec-Nature et technologies (FRQ-NT). The authors would also like to acknowledge the financial support of Safran Aircraft Engine for this work. This work is part of the Safran Industrial Research Chair on the Development of Sustainable Aero-Propulsion Systems at ÉTS.

ORCID iDs

Sébastien Cantin

Mohamed Chouak

François Morency

References

Lee

Fahey

Forster

, et al. Aviation and global climate change in the 21st century. Atmos Environ 2009; 43: 3520–3537.

Schumann

Heymsfield

. On the life cycle of individual contrails and contrail cirrus. Meteorol Monogr 2017; 58: 3.1–3.24.

Nguyen

Nguyen-Tri

Vancassel

Garnier

. Aero-thermodynamic and chemical process interactions in an axial high-pressure turbine of aircraft engines. Int J Engine Res 2019; 20: 653–669.

Lewellen

Meza

Huebsch

. Persistent contrails and contrail cirrus. Part I: Large-eddy simulations from inception to demise. J Atmos Sci 2014; 71: 4399–4419.

Jensen

Toon

Kinne

, et al. Environmental conditions required for contrail formation and persistence. J Geophys Res: Atmos 1998; 103: 3929–3936.

Solomon

. IPCC (2007): climate change the physical science basis. Agenda 2007; 6(7): 333.

Kärcher

. Formation and radiative forcing of contrail cirrus. Nat Commun 2018; 9: 1824.

Schumann

. On conditions for contrail formation from aircraft exhausts. Meteorol Z 1996; 5(1): 4–23.

Schumann

. A contrail cirrus prediction model. Geosci Model Dev 2012; 5: 543–580.

10.

Schumann

Mayer

Graf

Mannstein

. A parametric radiative forcing model for contrail cirrus. J Appl Meteorol Climatol 2012; 51: 1391–1406.

11.

Kärcher

. Physicochemistry of aircraft-generated liquid aerosols, soot, and ice particles: 1. Model description. J Geophys Res: Atmos 1998; 103: 17111–17128.

12.

Wong

Miake-Lye

. Parametric studies of contrail ice particle formation in jet regime using microphysical parcel modeling. Atmos Chem Phys 2010; 10(7): 3261–3272.

13.

Gierens

. On the transition between heterogeneous and homogeneous freezing. Atmos Chem Phys 2003; 3: 437–446.

14.

Rojo

Vancassel

Mirabel

Ponche

Garnier

. Impact of alternative jet fuels on aircraft induced aerosols. Fuel 2015; 144: 335–341.

15.

Paoli

Vancassel

Garnier

Mirabel

. Large-eddy simulation of a turbulent jet and a vortex sheet interaction: Particle formation and evolution in the near field of an aircraft wake. Meteorol Z 2008; 17(2): 131–144.

16.

Vancassel

Mirabel

Garnier

. Numerical simulation of aerosols in an aircraft wake using a 3D LES solver and a detailed microphysical model. Int J Sustainable Aviation 2014; 1: 139–159.

17.

Unterstrasser

. Properties of young contrails – a parametrisation based on large-eddy simulations. Atmos Chem Phys 2016; 16: 2059–2082.

18.

Burkhardt

Karcher

. Process-based simulation of contrail cirrus in a global climate model. J Geophys Res: Atmos 2009; 114: D16201.

19.

Jacobson

Wilkerson

Naiman

Lele

. The effects of aircraft on climate and pollution. Part I: numerical methods for treating the subgrid evolution of discrete size- and composition-resolved contrails from all commercial flights worldwide. J Comput Phys 2011; 230: 5115–5132.

20.

Lewellen

. A large-Eddy simulation study of contrail ice number formation. J Atmos Sci 2020; 77: 2585–2604.

21.

Unterstrasser

Paoli

Sölch

Kühnlein

Gerz

. Dimension of aircraft exhaust plumes at cruise conditions: effect of wake vortices. Atmos Chem Phys 2014: 14(5): 2713–2733.

22.

Garnier

Maglaras

Morency

Vancassel

. Effect of compressibility on contrail ice particle growth in an engine jet. Int J Turbo Jet Engines 2014; 31: 131–140.

23.

Cantin

Morency

Garnier

. Modeling capabilities on the formation of contrails in a commercial CFD code. In: ISABE 2019 (ed Engines ISfAB), Canberra, Autralia, 22–27 September 2019, p.14.

24.

Zhang

Chen

. Comparison of the Eulerian and Lagrangian methods for predicting particle transport in enclosed spaces. Atmos Environ 2007; 41: 5236–5248.

25.

Lewellen

. Effects of aircraft wake dynamics on measured and simulated NOx and HOx wake chemistry. J Geophys Res 2001; 106: 27661–27672.

26.

Paugam

Paoli

Cariolle

. Influence of vortex dynamics and atmospheric turbulence on the early evolution of a contrail. Atmos Chem Phys 2010; 10: 3933–3952.

27.

Khou

Ghedhaifi

Vancassel

Garnier

. Spatial simulation of contrail formation in near-field of commercial aircraft. J Aircraft 2015; 52(6): 1927–1938.

28.

Khou

Ghedhaifi

Vancassel

Montreuil

Garnier

. CFD simulation of contrail formation in the near field of a commercial aircraft: effect of fuel sulfur content. Meteorol Z 2017; 26(6): 585–596.

29.

Unterstrasser

Sölch

. Optimisation of the simulation particle number in a Lagrangian ice microphysical model. Geosci Model Dev 2014; 7: 695–709.

30.

Unterstrasser

Sölch

. Study of contrail microphysics in the vortex phase with a Lagrangian particle tracking model. Atmos Chem Phys 2010; 10: 10003–10015.

31.

Khain

Beheng

Heymsfield

, et al. Representation of microphysical processes in cloud-resolving models: Spectral (bin) microphysics versus bulk parameterization. Rev Geophysics 2015; 53: 247–322.

32.

Paoli

Nybelen

Picot

Cariolle

. Effects of jet/vortex interaction on contrail formation in supersaturated conditions. Phys Fluids 2013; 25: 053305.

33.

Unterstrasser

Görsch

. Aircraft-type dependency of contrail evolution. J Geophys Res: Atmos 2014; 119: 14015–14027.

34.

CD-Adapco. STAR-CCM+ Documentation Version 13.06.012. 2018.

35.

Kärcher

. Transport of exhaust products in the near trail of a jet engine under atmospheric conditions. J Geophys Res: Atmos 1994; 99: 14509–14517.

36.

Jones

Launder

. The prediction of laminarization with a two-equation model of turbulence. Int J Heat Mass Transfer 1972; 15: 301–314.

37.

Shaheed

Mohammadian

Kheirkhah Gildeh

. A comparison of standard k–ε and realizable k–ε turbulence models in curved and confluent channels. Environ Fluid Mech 2019; 19: 543–568.

38.

Thies

Tam

CKW

. Computation of turbulent axisymmetric and nonaxisymmetric jet flows using the K-epsilon model. AIAA J 1996; 34: 309–316.

39.

Schröder

Brock

Baumann

, et al. In situ studies on volatile jet exhaust particle emissions: impact of fuel sulfur content and environmental conditions on nuclei mode aerosols. J Geophys Res: Atmos 2000; 105: 19941–19954.

40.

Schröder

Karcher

Duroure

, et al. On the transition of contrails into cirrus clouds. J Atmos Sci 2000; 57: 464–480.

41.

Paoli

Hélie

Poinsot

. Contrail formation in aircraft wakes. J Fluid Mech 2004; 502: 361–373.

42.

Koehler

DeMott

Kreidenweis

, et al. Cloud condensation nuclei and ice nucleation activity of hydrophobic and hydrophilic soot particles. Phys Chem Chem Phys 2009; 11: 7906–7920.

43.

Hartmann

Niedermeier

Voigtlaender

, et al. Homogeneous and heterogeneous ice nucleation at LACIS: operating principle and theoretical studies. Atmos Chem Phys 2011; 11(4): 1753–1767.

44.

Fukuta

Walter

. Kinetics of hydrometeor growth from a vapor-spherical model. J Atmos Sci 1970; 27: 1160–1172.

45.

Patankar

Spalding

. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int J Heat Mass Transfer 1972; 15: 1787–1806.

46.

Cantin

. Emissions and Microphysical Modelling for Aircraft Plume (EMMAP), 2020, https://github.com/SebastienCantin/EMMAP (accessed 4 February 2021).

47.

Guitton

Tinney

Jordan

, et al. Measurements in a co-axial subsonic jet. In: 45th AIAA aerospace sciences meeting and exhibit, Reno, VA, USA, American Institute of Aeronautics and Astronautics, 8–11 January 2007.

48.

Salvetti

Lombardi

. Three-dimensional effects on mixing in a coaxia jet configuration. In: Advances in turbulence VII: Proceedings of the seventh European turbulence conference, held in Saint-Jean Cap Ferrat, France, 30 June – 3 July 1998 / Actes de la Septième Conférence Européenne de Turbulence, tenue à Saint-Jean Cap Ferrat, France, 30 Juin – 3 Juillet 1998 (ed Frisch

), Dordrecht, Springer Netherlands, 30 June–3 July 1998, pp.167–170.

49.

Celik

Ghia

Roache

Freitas

. Procedure for estimation and reporting of uncertainty due to discretization in CFD applications. J Fluids Eng 2008; 130: 078001.

50.

Shih

T-H

Liou

Shabbir

Yang

Zhu

. A new k-ϵ eddy viscosity model for high Reynolds number turbulent flows. Comput Fluids 1995; 24: 227–238.

51.

Menter

. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J 1994; 32: 1598–1605.

52.

Timmerman

Skeen

Bryanston-Cross

Graves

. Large-scale time-resolved digital particle image velocimetry (TR-DPIV) for measurement of high subsonic hot coaxial jet exhaust of a gas turbine engine. Meas Sci Technol 2009; 20: 074002.

53.

Garnier

Baudoin

Woods

Louisnard

. Engine emission alteration in the near field of an aircraft. Atmos Environ 1997; 31: 1767–1781.

54.

AnsysInc. ANSYS ICEM CFD user manual. 2019.

55.

Nair

Suryan

Kim

. Computational study on flow through truncated conical plug nozzle with base bleed. Propul Power Res 2019; 8: 108–120.

56.

Goulos

Stankowski

MacManus

Woodrow

Sheaf

. Civil turbofan engine exhaust aerodynamics: Impact of bypass nozzle after-body design. Aerosp Sci Technol 2017; 73: 85–95.

57.

Xiong

Nielsen

Liu

Papamoschou

. Computation of high-speed coaxial jets with fan flow deflection. AIAA J 2010; 48: 2249–2262.

58.

Schumann

Schlager

Arnold

Baumann

Haschberger

Klemm

. Dilution of aircraft exhaust plumes at cruise altitudes. Atmos Environ 1998; 32: 3097–3103.

59.

Schumann

Ström

Busen

, et al. In situ observations of particles in jet aircraft exhausts and contrails for different sulfur-containing fuels. J Geophys Res: Atmos 1996; 101: 6853–6869.

60.

Khou

. Modélisation des traînées de condensation par interaction entre l’aérodynamique, la cinétique chimique et la microphysique. PhD Dissertation, Université Pierre et Marie Curie, ONERA/EMI, Palaiseau, 2016.

61.

Kärcher

Mayer

Gierens

Burkhardt

Mannstein

Chatterjee

. Aerodynamic contrails: microphysics and optical properties. J Atmos Sci 2009; 66: 227–243.

62.

Naiman

Lele

Jacobson

. Large eddy simulations of contrail development: Sensitivity to initial and ambient conditions over first twenty minutes. J Geophys Res: Atmos 2011; 116: 12

63.

Busen

Schumann

. Visible contrail formation from fuels with different sulfur contents. Geophys Res Lett 1995; 22: 1357–1360.

64.

Pruppacher

Klett

Wang

. Microphysics of clouds and precipitation. New York, NY: Springer Publications, 1998, p.381.

65.

Fukuta

Myers

. Simultaneous measurement of condensation and thermal accommodation coefficients for cloud droplet growth in due consideration of a new moving surface-boundary effect. J Atmos Sci 2007; 64: 955–968.

66.

Taleb

D-E

. Etude microphysique et modelisation des sillages d’avions subsoniques et supersoniques. PhD Dissertation, 1997. Université Louis Pasteur, Strasbourg, France.

67.

Goff

. Saturation pressure of water on the new Kelvin temperature scale. Trans Am Soc Heat Vent Eng 1957; 63: 347–354.