2.1. Governing equations

The flow is assumed to be incompressible and isothermal, governed by the Navier–Stokes equations

(2.1)\begin{gather} \frac{\partial u_i^*}{\partial t^*} + u_j^*\,\frac{\partial u_i^*}{\partial x_j^*}={-}\frac{\partial p^*}{\partial x_i^*}+\frac{1}{Re}\,\frac{\partial^2 u_i^*}{{\partial x_j^*}\,{\partial x_j^*}}, \end{gather}

(2.2)\begin{gather} \frac{\partial u_j^*}{\partial x_j^*}=0, \end{gather}

where $u_i^*$ is the Cartesian velocity component in the $i$th direction, and $p^*$ is the pressure. The asterisk denotes non-dimensional variables, defined as

(2.3a–d)\begin{equation} u_i^*=\frac{u_i}{U_o}, \quad p^*=\frac{p}{\rho U_o^2}, \quad x_i^*=\frac{x_i}{h}, \quad t^* = \frac{t}{h/U_o}, \end{equation}

where $U_o$ is the jet velocity, $h$ is the jet width at the inlet, and $\rho$ is the fluid density. The Reynolds number is defined as $Re=U_o h/\nu$, where $\nu$ is the kinematic viscosity.

The nanoparticle dynamics is described by the PBE, which assumes the form (Friedlander Reference Friedlander2000)

(2.4)$$\begin{gather} \frac{\partial n}{\partial t} + u_j\,\frac{\partial n}{\partial x_j}= \frac{\partial }{\partial x_j} \left( D_p\,\frac{\partial n}{\partial x_j} \right) + \frac{1}{2} \int_0^v \beta(w,v-w)\,n(w)\,n(v-w) \,{\rm d}w\nonumber\\ -\int_0^\infty \beta(v,w)\,n(v)\,n(w) \,{\rm d}w, \end{gather}$$

where $v$ and $w$ denote particle volume, $n=n( {\boldsymbol x},t,v)$ is the number density function per unit volume of fluid (such that $n\,{\rm d}v$ is the concentration of particles with volume between $v$ and $v+{\rm d}v$), $D_p$ is the diffusion coefficient of particles of volume $v$ yielded by the Stokes–Einstein equation (Friedlander Reference Friedlander2000), and $\beta (v,w)$ is the coagulation kernel, a function that represents the frequency of collisions that lead to coagulation events and will be discussed further below. Coagulation is described by the two integral terms in (2.4), which denote particle birth and death. The birth term accounts for all possible particle pairs whose coalescence results in particles of volume $v$, while the death term accounts for the disappearance of particles of volume $v$ due to coagulation with particles of any volume. The PBE can include additional mechanisms representing nucleation, growth and breakage, but the focus of the present study is on coagulation.

When the particle diameter is much smaller than the mean free path of the carrier gas, the particles can be assumed to behave like rigid elastic spheres, and an expression for the collision frequency $\beta (v,w)$ can be derived from the kinetic theory of gases. In this case, coagulation takes place in the free-molecule regime and the kernel is

(2.5)\begin{equation} \beta (v,w) = A_1 \left( \frac{1}{v}+\frac{1}{w} \right)^{1/2} \left( v^{1/3} + w^{1/3} \right)^2, \end{equation}

with the factor $A_1$ defined as

(2.6)\begin{equation} A_1 = \left( \frac{3}{4{\rm \pi}} \right)^{1/6} \left( \frac{6 k_b T}{\rho_p} \right)^{1/2}, \end{equation}

where $\rho _p$ is the particle density, $T$ is the temperature, and $k_b$ is the Boltzmann constant.

For future reference, we also present the transport equation for the moments of the distribution $n=n( {\boldsymbol x},t,v)$. The $k$th moment, $M_k$, is defined as

(2.7)\begin{equation} M_k = \int_0^\infty v^k\,n(v) \,{\rm d}v. \end{equation}

Some of the low-order moments have a physical interpretation. For example, the zeroth moment, $M_0$, is the total number concentration of particles, while the first moment, $M_1$, represents the particle volume fraction. Furthermore, when the diameter of the particles is smaller than the wavelength of the incident light, the second moment, $M_2$, is proportional to the total light scattering by particles (Friedlander Reference Friedlander2000). The following transport equation for the evolution of $M_k$ can be derived (Williams & Loyalka Reference Williams and Loyalka1991; Whitby & McMurry Reference Whitby and McMurry1997):

(2.8)$$\begin{gather} \frac{\partial M_k}{\partial t}+u_j\,\frac{\partial M_k}{\partial x_j} = \frac{\partial}{\partial x_j} \left( D_p\,\frac{\partial M_k}{\partial x_j} \right) + \frac{1}{2} \int_0^\infty \int_0^\infty (v+w)^k\,\beta(v,w)\,n(v)\,n(w) \,{\rm d}v \,{\rm d}w \nonumber\\ - \frac{1}{2} \int_0^\infty \int_0^\infty (v^k+w^k)\,\beta(v,w)\,n(v)\,n(w)\,{\rm d}v\,{\rm d}w. \end{gather}$$

In the absence of convection and diffusion (i.e. in a perfectly stirred reactor), it can be shown that (2.4) admits self-similar solutions. It must be noted that $M_1$ is constant under these conditions. If the number density function and particle volume are scaled as

(2.9)\begin{gather} \psi = \frac{M_1}{M_0^2}\,n, \end{gather}

(2.10)\begin{gather} \eta = \frac{M_0}{M_1}\,v, \end{gather}

and substituted in (2.4), then the solution reaches a time-independent distribution $\psi (\eta )$ after a sufficient amount of time. The initial particle distribution affects the time to reach the self-preserving (SP) state but not the distribution itself (Friedlander Reference Friedlander2000). For an initially monodisperse aerosol undergoing coagulation in the free-molecule regime, the time to reach the self-similar distribution is (Landgrebe & Pratsinis Reference Landgrebe and Pratsinis1990; Vemury, Kusters & Pratsinis Reference Vemury, Kusters and Pratsinis1994)

(2.11)\begin{equation} \tau_{SP}= \frac{\tau_f}{ A_1 v_o^{1/6} N_0 }, \end{equation}

where $\tau _f \approx 5$ (Vemury et al. Reference Vemury, Kusters and Pratsinis1994), $v_o$ is the initial particle volume, and $N_0$ is the initial particle concentration.

Using (2.3a–d) and the non-dimensional variables

(2.12a,b)\begin{equation} n^*=\frac{n}{n_o}, \quad v^*=\frac{v}{v_o}, \end{equation}

where $n_o,v_o$ are the number density and volume of particles at the inlet, respectively, and substituting into (2.4), we obtain

(2.13)$$\begin{gather} \frac{\partial n^*}{\partial t^*} \!+\! u_j^*\,\frac{\partial n^*}{\partial x_j^*}= \frac{1}{Sc}\, \frac{1}{Re}\,\frac{\partial^2 n^*}{{\partial x_j^*}\,{\partial x_j^*}} \!+\! Da_{coag} \left(\frac{1}{2} \int_0^{v^*} \beta^*(w^*,v^*\!-\!w^*)\,n^*(w^*)\,n^*(v^*\!-\!w^*) \,{\rm d}w^* \right.\nonumber\\ \left. -\int_0^\infty \beta^*(v^*,w^*)\,n^*(v^*)\,n^*(w^*) \,{\rm d}w^*\right), \end{gather}$$

where $\beta ^*(v^*,w^*) = ( {1}/{v^*}+{1}/{w^*} )^{1/2} ( {v^*}^{1/3} + {w^*}^{1/3} )^2$ is the normalised collision function, and $Sc=\nu /D_p$ is the particle Schmidt number (defined as the ratio of the kinematic viscosity to the particle diffusion coefficient). The factor $Da_{coag}$ is the coagulation Damköhler number, which is equal to the ratio of a convection time scale to a coagulation time scale. The time scales are defined here as in Garrick et al. (Reference Garrick, Lehtinen and Zachariah2006), and for the case of the free-molecule regime, $Da_{coag}$ is given by

(2.14)\begin{equation} Da_{coag}=\frac{\tau_{conv}}{\tau_{coag}}=\frac{A_1 h N_0}{U_o}\,v_o^{1/6}, \end{equation}

where the convection time scale is $\tau _{conv} = h/U_o$, and the coagulation time scale is $\tau _{coag} = (A_1 N_0 v_o^{1/6})^{-1}$. For the case of an initially monodisperse aerosol, (2.11) shows that the coagulation time scale is one-fifth of the time needed to reach the self-preserving distribution. Based on (2.11) and (2.14), the coagulation Damköhler number can be written as

(2.15)\begin{equation} Da_{coag} = \frac{5 \tau_{conv}}{\tau_{SP}}. \end{equation}

The assumptions of the approach adopted in the present paper are listed below.

1. (i) Coagulation in the free-molecule regime is considered, which means that particles should be smaller than the mean free path in the carrier fluid. It must be mentioned that our simulations are carried out in terms of dimensionless variables, but if we assume, for example, that the particles at the inlet of the jet have diameter 5 nm, then the maximum particle mean size in the jet would be around 25 nm, which is less than the mean free path in air at normal temperature and pressure (approximately 65 nm).

2. (ii) The particles have no effect on the fluid motion and do not affect the fluid properties (density, viscosity etc.). If we assume that the particles at the inlet have diameter 4–5 nm, the jet nozzle width is roughly 1 cm and the temperature of the jet is close to $50\,^{\circ }{\rm C}$, conditions not far from those in the experimental work of Junzong et al. (Reference Junzong, Haiying and Jinsheng2013), then the particle volume fraction at the inlet is less than $10^{-6}$. This value is small enough to ensure the validity of the assumption (Elghobashi Reference Elghobashi2006).

3. (iii) The particles are spherical due to instantaneous coagulation. While we are not considering a particular physical system, this is usually true for nanoparticle coagulation in the free-molecule regime.

4. (iv) The flow is assumed to be isothermal.

5. (v) Electrostatic effects (Van der Waals forces) are not taken into account.

6. (vi) The Stokes number $St$ – which is the ratio of a characteristic time scale of the particles, $\tau _p = \rho _p d_p^2 /(18 \mu )$ (Friedlander Reference Friedlander2000), to that of the flow, $\tau _f$ – is much smaller than unity and therefore the particles follow the fluid. This assumption is satisfied easily for nanoparticles, and an example is given here to justify this assertion. If we assume that particles at the inlet have diameter 5 nm, then the largest particles found in the jet (at the right edge of the PSD with very low number density) would have diameter 200 nm (as will be found later in this study). Then $\tau _p$ can be calculated by using $\rho _p = 2200\,{\rm kg}\,{\rm m}^{-3}$ , which is a typical value of density for silica nanoparticles, and $\mu = 1.82 \times 10^{-5}$. The smallest time scale in our problem is the Kolmogorov time scale, $\tau _{\eta } = \sqrt {\nu / \epsilon }$, and this is taken as $\tau _f$ in this example. The smallest non-dimensional value of $\tau _{\eta }$ calculated by the present DNS was $0.38$, while the way the variables are non-dimensionalised has been shown in (2.3a–d). If we assume a jet nozzle width $h \approx 1\,{\rm cm}$, then the corresponding dimensional value of $\tau _{\eta }$ would be $7.59 \times 10^{-4}$. Based on these assumptions, the estimated Stokes number would be $St \approx 4 \times 10^{-4} \ll 1$.

7. (vii) A Schmidt number $Sc=1$ is assumed for all particle classes. In reality, $Sc$ takes a range of values because it depends on the particle diffusion coefficient $D_p$, which in turn is a function of particle size according to the Stokes–Einstein equation (Friedlander Reference Friedlander2000). For temperature 298 K, for example, particles with 20 nm diameter have Schmidt number 1415. The use of $Sc=1$ implies higher diffusive transport, and studies of flows with high $Sc$ have so far been conducted only in the context of 2-D flows, such as the study of Garrick & Khakpour (Reference Garrick and Khakpour2004) on a 2-D mixing layer and the study of Cifuentes et al. (Reference Cifuentes, Sellmann, Wlokas and Kempf2020) on a 2-D flame–vortex interaction. However, a high $Sc$ results in Batchelor scales, $\eta _c=\eta /{\sqrt {Sc}}$, that are much smaller than the Kolmogorov scale, $\eta$ (Davidson Reference Davidson2004). For example, for $Sc=1415$, $\eta _c$ would be roughly 37 times smaller than $\eta$, hence in a 3-D flow, it would require a grid with $37^3$ more cells. For the present simulation, which currently uses $52$ million cells, this would amount to 2.6 trillion cells. Therefore, the study of the effect of high $Sc$ in a 3-D flow simulation has to be left for future work.

2.2. Numerical method

The simulations were performed with our computational fluid dynamics (CFD) code, PANTARHEI, which is coupled with the in-house code for population balance modelling, CPMOD. PANTARHEI has been used to simulate a number of transitional and turbulent flows (Xiao & Papadakis Reference Xiao and Papadakis2019; Alves Portela, Papadakis & Vassilicos Reference Alves Portela, Papadakis and Vassilicos2020). The incompressible Navier–Stokes equations are discretised by applying the finite volume method in unstructured grids, and all unknown variables are stored at the cell centroids (collocated variable arrangement). A second-order central discretisation scheme is employed for the convection and viscous terms, while a third-order backward differentiation formula (BDF) is used for the transient term. The orthogonal viscous terms are treated implicitly, while the convection term and the non-orthogonal diffusion terms are extrapolated from the previous three time steps. The pressure–velocity coupling is dealt with by the fractional step method. The code is parallelised with the aid of the PETSc and HYPRE libraries (Balay et al. Reference Balay2021).

A method developed recently by Liu & Rigopoulos (Reference Liu and Rigopoulos2019) is employed for the discretisation of the PBE. The method (a) provides an accurate prediction of the distribution with a small number of sections, (b) conserves the first moment (or any other single moment) in coagulation, (c) is applicable to any non-uniform grid (even an adaptive one), and (d) is efficient and robust when coupled with CFD. It is based on the framework of finite-volume methods applied to the particle volume domain, $v$, truncated to $(v_o,v_{max})$. The latter is discretised in $(v_o,\ldots,v_k)$, the intervals $\Delta v_i=v_i-v_{i-1}$, $i=1,\ldots, k$ are defined, and each interval corresponds to a number density $n_i$. One of the main challenges in the numerical solution of the PBE is the conservation of the first moment in coagulation, due to the presence of the integral terms. In our method, a transformation is employed to decompose the birth term into fluxes of coagulating particles, and a map is constructed to record them. An algorithm for constructing this map in an arbitrary non-uniform grid has been developed, and the aggregation birth and death terms are then calculated via these fluxes, resulting in mass conservation without introducing correction terms. The map is constructed once at the start of the simulation and accessed throughout the run, making the method computationally efficient and easy to couple with CFD. Test cases, including analytical solutions and a co-flow laminar flame (Liu & Rigopoulos Reference Liu and Rigopoulos2019), have shown that the method produces very accurate solutions that conserve the first moment, even on coarse grids.

After discretisation in the volume domain, the PBE is converted into a coupled system of partial differential equations that govern the temporal and spatial evolution of the PSD in terms of discretised number densities, $n_i$. The transport equation for $n_i$ is

(2.16)\begin{equation} \frac{\partial n_i}{\partial t} + \frac{\partial (u_j n_i)}{\partial x_j} = \frac{\partial}{\partial x_j} \left( D_p\,\frac{\partial n_i}{\partial x_j} \right) + {S_i}(n_1,\ldots,n_k), \end{equation}

where the source term, ${S_i}(n_1,\ldots,n_k)$, represents the effect of coagulation for particles in the interval $i$. The numerical scheme used for the solution of these transport equations is slightly different from that described previously for the flow field. More specifically, the convection and viscous terms are treated explicitly using third-order extrapolation from the previous three time steps. The values of $n_1,\ldots,n_k$ are also extrapolated similarly and used to calculate ${S_i}$. In order to ensure boundedness for $n_i$, the Gamma differencing scheme (Jasak, Weller & Gosman Reference Jasak, Weller and Gosman1999) was employed for the discretisation of the convection term, and the value 0.1 was used for the scheme parameter $\beta _m$.

2.3. Flow configuration and computational setup

Direct numerical simulations of a 3-D spatially developing planar turbulent jet were performed. The jet is laden with monodisperse nanoparticles and flows into a particle-free co-flow stream. There is a wealth of data for validating the flow field (Rajaratnam Reference Rajaratnam1976; Le Ribault, Sarkar & Stanley Reference Le Ribault, Sarkar and Stanley1999; Pope Reference Pope2000; Stanley, Sarkar & Mellado González Reference Stanley, Sarkar and Mellado González2002; Klein, Sadiki & Janicka Reference Klein, Sadiki and Janicka2003a,Reference Klein, Sadiki and Janickab; Suresh et al. Reference Suresh, Srinivasan, Sundararajan and Das2008; Watanabe et al. Reference Watanabe, Sakai, Nagata, Ito and Hayase2014).

The flow configuration is similar to that of Stanley et al. (Reference Stanley, Sarkar and Mellado González2002). The Reynolds number is $Re= {U_o h}/{\nu } = 3000$, and the ratio of co-flow velocity, $U_\infty$, to the jet velocity, $U_o$, is ${U_\infty }/{U_o}=0.2$. The $x$, $y$ and $z$ coordinates (corresponding to indices $i=1,2,3$ in the notation of (2.1)) are defined in the streamwise, cross-stream and spanwise directions, respectively. Turbulent properties are statistically homogeneous in the $z$ direction. The size of the computational domain $L_x \times L_y \times L_z$ is $25h \times 26h \times 5h$, similar to that of Stanley et al. (Reference Stanley, Sarkar and Mellado González2002) and Das & Garrick (Reference Das and Garrick2010). The computational domain is long enough in the streamwise direction for the flow to reach fully developed turbulence. In addition, the domain is extended in the cross-stream direction in order to minimise the effect of the top and bottom boundaries.

A convective boundary condition was applied at the outlet boundary in the streamwise direction; this boundary condition allows vortices to propagate through the exit plane with minimal reflection. Periodic boundary conditions were used in the spanwise direction. Finally, a symmetry boundary condition was applied in the top and bottom boundaries. This condition does not allow for mass entrainment, and this causes a small reduction in the local co-flow velocity in the streamwise direction. To mitigate this effect, the computational domain was significantly extended in the cross-stream direction. Good agreement between the DNS and reference data from the literature, to be presented later, in § 3.1, suggests that the effect of the boundaries is very small.

At inflow, a top-hat profile with smooth edges is prescribed for the mean (i.e. time-averaged) velocity. The profile is given by (Klein et al. Reference Klein, Sadiki and Janicka2003b; Stanley et al. Reference Stanley, Sarkar and Mellado González2002)

(2.17)\begin{equation} \overline{U}(y)= \frac{U_o + U_{\infty}}{2} + \frac{U_o - U_{\infty}}{2} \tanh{\frac{-{\left|y\right|}+h/2}{2 \theta_o}}, \end{equation}

where $\theta _o = h/20$ (Le Ribault et al. Reference Le Ribault, Sarkar and Stanley1999) is the shear layer momentum thickness. The mean cross-stream and spanwise velocities were set to zero. The mean profile of the number density for the first interval $n_1$ is given by the tangential profile

(2.18)\begin{equation} \overline{n_1}(y)= {n_1}_o + {n_1}_o \tanh{\frac{-{\left|y\right|}+h/2}{2 \theta_o}}, \end{equation}

where ${n_1}_o$ is the time-averaged value of $n_1$ at the centreline.

The method of Klein et al. (Reference Klein, Sadiki and Janicka2003a) was employed to generate artificial turbulent fluctuations at the inlet. The method requires a mean velocity profile, a profile for the Reynolds stresses, turbulence length and time scales, a filter width, and the number of slices (snapshots) of data to be generated. The output of the algorithm is a 3-D matrix that represents the velocity components at the inlet plane for different time snapshots. In the present study, the mean velocity profile is given by (2.17). Strictly speaking, the r.m.s. profile in the $x$ direction, $u_{rms}$, should be zero at $y={\pm }h/2$. However, it grows quickly due to high shear, and for this reason, $u_{rms}/U_o$ is set to $1.25\,\%$ at the jet centreline, and increases parabolically to $5\,\%$ at $y={\pm }h/2$. The turbulence intensity of the co-flow is set to $1\,\%$. This $u_{rms}$ profile is similar to that employed by Stanley et al. (Reference Stanley, Sarkar and Mellado González2002). The velocity r.m.s. profiles in the $y$ and $z$ directions were set equal to $u_{rms}$. The turbulence length scale was $l_{Turb}=0.3h$, and the turbulence time scale was $t_{Turb}=10h/U_o$, where it was assumed that $t_{Turb}=l_{Turb}/U_{rms}$, with $U_{rms}=3\,\% U_o$. A discussion on the value and the effect of the turbulence length scale can be found in Klein et al. (Reference Klein, Sadiki and Janicka2003a). The number of snapshots was equal to 100, which is large enough to generate data for one flow-through time (defined using the length of the computational domain and the co-flow velocity, $T=L_x/U_\infty$). The number of snapshots was calculated based on the prescribed turbulence length and time scales mentioned above.

The simulations were performed on a non-uniform computational grid with approximately 52 million cells, distributed as $900 \times 452 \times 128$ in the $x,y,z$ directions, respectively. The grid spacing was uniform in the $x$ and $z$ directions, and equal to $\Delta x=0.028 h$ and $\Delta z=0.039 h$. In the $y$ direction, it was very fine close to the centreline, where $\Delta y =0.018 h$, and increased exponentially towards the boundaries, where $\Delta y =0.124 h$. The time step was $\Delta t=0.005h/U_o$, resulting in a maximum Courant–Friedrichs–Lewy (CFL) number of 0.24. Convergence studies with a perfectly stirred reactor showed that this time step is adequate for the coagulation dynamics for the conditions considered.

Each simulation was run first for one flow-through time $T$ to reach a state that is independent of the initial conditions (recall that $T$ is defined using the co-flow velocity $U_\infty$, which is five times smaller than $U_o$). Following this, the simulations were run for 80 000 time steps (corresponding to $3.2T$), and statistics were computed by averaging in time and in the homogeneous spanwise direction. Finally, where possible, exploitation of the jet symmetry with respect to the centreline was used to double the sample size. The number of cores was 1008, and 84 672 Cray XC30 CPU hours were required for each case.

The maximum value of the ratio of cell size, defined as $\Delta _{cell} = ( \Delta x\,\Delta y\,\Delta z )^{1/3}$, to the Kolmogorov length scale $\eta = (\nu ^3/\epsilon )^{1/4}$ in the domain was $2.8$, while the mean value in the region of fully developed turbulence was $1.95$. The grid resolution was assessed by checking the balance of the turbulent kinetic energy equation; results are shown in Appendix A, figure 17. The terms of the equation balance well, and the maximum error was less than $5.5\,\%$ of the peak absolute value of convection due to mean flow, indicating that the grid resolution is fine enough. The one-dimensional spectra of the velocity fluctuations for a selected point in the fully turbulent region are also presented in Appendix A, figure 18. The spectra exhibit a $-5/3$ power-law slope over about a decade. There is also a dissipative range with faster decay rate at high frequencies, indicating that the turbulent motions are well resolved. The maximum Reynolds number $Re_\lambda = {\sqrt {\overline {{\left |{\boldsymbol u'}\right |}^2}}\lambda }/{\nu }$, based on the Taylor micro-scale $\lambda = \sqrt {{\nu {\overline {{\left |{\boldsymbol u'}\right |}^2}}}/{\epsilon }}$, was found to be 160.

Simulations for two different coagulation Damköhler numbers, $Da_{coag}=1$ and $Da_{coag}=1/3$, were performed. The values of $Da_{coag}$ are close to those examined in previous studies of nanoparticle coagulation (Garrick et al. Reference Garrick, Lehtinen and Zachariah2006; Garrick Reference Garrick2011). The two cases correspond to different times for the PSD to reach the self-preserving distribution. According to (2.15), in a perfectly stirred reactor, these times would be $\tau _{SP}=5$ and $\tau _{SP}=15$, respectively. As will be shown later, the jet breaks up at $x/h \approx 7.5$. For the chosen values of $Da_{coag}$, the PSD becomes self-similar upstream and downstream of the break-up location for the first and second cases, respectively.

The PBE was discretised with 35 intervals in the particle volume domain. An exponential grid was employed for all cases; this is well suited for fast-growing particles. The range $(v_o,\ldots,v_{max})$ covered was different for each $Da_{coag}$ case; $v_o$ remained the same but the maximum value was $v_{max}=11\,250 v_o$ for $Da_{coag}=1$, and $v_{max}=2800 v_o$ for $Da_{coag}=1/3$. Convergence studies with a perfectly stirred reactor showed that this grid is sufficiently accurate and ensures that the size distribution does not grow beyond the maximum value during the simulation time. The computational parameters of the simulations are also summarised in table 1.

Table 1.

Computational parameters of the simulations.

The simulations were carried out in terms of dimensionless variables. The non-dimensionalisation of the equations is described in § 2.1, where an asterisk was introduced to denote dimensionless variables; refer to (2.3a–d) and (2.12a,b). In the sections that follow, only dimensionless variables are presented, and for the sake of simplicity, the asterisk has been removed.

3.1. Flow field

In this section, the DNS data for the planar jet flow are presented and compared with experimental and other reference simulation data from the literature. Additional data for the balance of the terms in the turbulent kinetic energy equation and the one-dimensional spectra of turbulent fluctuations are given in Appendix A. The turbulent kinetic energy, $k$, in the ($x$–$y$) plane is depicted in figure 1(a) (because of symmetry, only the top half is shown). The jet breaks up at $x \approx 7.5$, and downstream of this location the turbulent kinetic energy rapidly grows and decays again; the maximum of $k$ is located at $(x,y)=(11,0.7)$. Contours of instantaneous spanwise vorticity and passive scalar fields, $\omega _z$ and $\phi$, are shown in figures 1(b) and 1(c), respectively. A downstream region with intense turbulent mixing and a wide range of turbulence length scales is evident.

Figure 1.

Flow field data: (a) turbulent kinetic energy (TKE) $k$; (b) 3-D view of the instantaneous spanwise vorticity field $\omega _z$; (c) instantaneous contours of the passive scalar on ($x$–$y$) plane $z=0$.

The growth of the jet half-width based on the streamwise velocity, $\delta _U$, is shown in figure 2(a). The jet half-width $\delta _U$ is defined as the distance from the symmetry plane (centreline) to the point in the cross-stream direction where the mean streamwise velocity excess with respect to the co-flow velocity, $U_e(x,y) = U(x,y)-U_\infty (x)$, is half that with respect to the centreline, $\Delta U_c(x) = U_c(x) - U_\infty (x)$. The ratio of the velocity excess at the jet nozzle exit, $\Delta U_o = U_o - U_\infty$, to $\Delta U_c(x)$ is shown in figure 2(b).

Figure 2.

Growth of the jet half-width based on the velocity (a) and centreline mean velocity excess decay (b). Blue line, present results; $\Box$, Thomas & Chu (Reference Thomas and Chu1989); $\diamond$, Browne et al. (Reference Browne, Antonia, Rajagopalan and Chambers1983); dotted line, DNS (Stanley et al. Reference Stanley, Sarkar and Mellado González2002).

The growth of the jet half-width based on the passive scalar, $\delta _\phi$, is shown in figure 3(a). The half-width $\delta _\phi$ is defined (similarly to $\delta _U$) as the distance from the symmetry plane (centreline) to the point where the passive scalar concentration is half of that on the centreline. The mean scalar decay, defined as $(\phi _o/\phi _c)^2$, where $\phi _c(x)$ is the value at the centreline and $\phi _o$ is the value at the jet nozzle exit, is shown in figure 3(b).

Figure 3.

Downstream growth of the jet half-width based on the passive scalar (a) and decay of the mean scalar on the jet centreline (b). Blue line, present results; $\diamond$, Browne et al. (Reference Browne, Antonia, Rajagopalan and Chambers1983); $\Box$, Davies, Keffer & Baines (Reference Davies, Keffer and Baines1975); dotted line, DNS (Stanley et al. Reference Stanley, Sarkar and Mellado González2002).

It can be seen that there is good agreement with the experimental and other reference data. Theoretical analysis (Rajaratnam Reference Rajaratnam1976; Pope Reference Pope2000) suggests that $\delta _U \propto x$ and $U_c \propto 1/\sqrt {x}$ for planar jet flow, and the same scaling is also valid for the scalar concentration. It has been documented (Stanley & Sarkar Reference Stanley and Sarkar2000; Stanley et al. Reference Stanley, Sarkar and Mellado González2002; da Silva, Lopes & Raman Reference da Silva, Lopes and Raman2015) that the jet break-up location is highly affected by the inlet conditions at the nozzle. The jet characteristics are also affected by the value of the Reynolds number in cases where $Re<10\,000$ (Klein et al. Reference Klein, Sadiki and Janicka2003b; Deo, Mi & Nathan Reference Deo, Mi and Nathan2008). Therefore, independently of the exact location of the jet break-up, the slope of the curves shown in figures 2 and 3 for the present DNS exhibits good agreement with that of reference data.

The self-similar cross-stream profiles for velocity and scalar are depicted in figures 4(a) and 4(b), respectively. The velocity profiles collapse to a single curve, and there is good agreement with reference (experimental or other DNS) data. The velocity field reaches self-similarity at $x \approx 12$, while the passive scalar is slightly slower to develop and the self-similar region is reached at $x \approx 17$. The same behaviour was also observed by Stanley et al. (Reference Stanley, Sarkar and Mellado González2002).

Figure 4.

Self-similar profiles of mean velocity and mean passive scalar for three axial positions. (a) Mean velocity: solid lines, present results; $\Box$, Gutmark & Wygnanski (Reference Gutmark and Wygnanski1976); $\circ$, Ramaprian & Chandrasekhara (Reference Ramaprian and Chandrasekhara1985); dotted line, DNS (Stanley et al. Reference Stanley, Sarkar and Mellado González2002). (b) Mean scalar: solid lines, present results; $\Box$, Davies et al. (Reference Davies, Keffer and Baines1975); $\circ$, Jenkins & Goldschmidt (Reference Jenkins and Goldschmidt1973); dotted line, DNS (Stanley et al. Reference Stanley, Sarkar and Mellado González2002).

3.2. Particle size distribution

The DNS–PBE approach allows for the prediction of the PSD at every point in the domain. As mentioned in § 2.1, in a perfectly stirred reactor, the PSD reaches a self-preserving distribution. Figure 5 shows the time-averaged PSDs at five different probe points on the centreline for both $Da_{coag}$ number cases. The computed distributions are scaled according to (2.9) and (2.10), and compared with the self-preserving PSD from Vemury & Pratsinis (Reference Vemury and Pratsinis1995). It can be seen that as the PSDs develop, they approach the self-preserving distribution and there is good agreement with the reference data, especially close to the exit of the domain. Turbulence may affect the time to reach the self-preserving state (Friedlander Reference Friedlander2000). As mentioned earlier, the values of $Da_{coag}$ were chosen in such a way that the PSD becomes self-preserving before and after the jet break-up location for the two cases. In particular, for the $Da_{coag}=1/3$ case, the PSD becomes self-preserving slightly downstream of the expected location ($x \approx 15$). However, the difference observed is small. The fact that the PSDs in the turbulent flow field do not deviate a lot from the self-preserving distribution can be attributed to the lack of recirculating flow structures that would bring into contact particles with different distributions and would make the PSD deviate from the self-preserving state. Finally, the good matching with the reference data indicates that the number and distribution of volume intervals are sufficient to resolve the PSD accurately.

Figure 5.

Normalised time-averaged PSDs at several probe points along the centreline. The results are compared with reference data for the self-preserving distribution (Vemury & Pratsinis Reference Vemury and Pratsinis1995). A logarithmic scale was used for the horizontal axes. (a) $Da_{coag}=1$; (b) $Da_{coag}=1/3$.