Statistical simulation of molecular diffusion effect on turbulent tetrad dispersion

International Journal of Heat and Mass Transfer 103 (2016) 87–98

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Statistical simulation of molecular diffusion effect on turbulent tetrad dispersion Fei Fei, Zhaohui Liu ⇑, Chuguang Zheng State Key Laboratory of Coal Combustion, School of Energy and Power Engineering, Huazhong University of Science and Technology, 430074 Wuhan, China

a r t i c l e

i n f o

Article history: Received 31 January 2016 Received in revised form 1 June 2016 Accepted 1 July 2016

Keywords: Fokker–Planck model Multi-scale numerical scheme Turbulent dispersion Differential diffusion

a b s t r a c t Molecular diffusion can significantly affect the dispersion and mixing processes at small temporal and spatial scales in turbulent flows, especially when differential diffusion among species is considered. Since species mixing takes place at the molecular level, it is very beneficial to study and model this phenomenon from a kinetic viewpoint. However, most molecular simulation methods, such as the direct simulation Monte Carlo (DSMC) method, are inefficient for high Reynolds or low Knudsen number flows. The diffusive information preservation (D-IP) method, which based on the Fokker–Plank model, can overcome this difficulty. In the current paper, we firstly validate the D-IP method by simulating the decaying homogeneous isotropic turbulence, and then apply it to investigate the effects of molecular diffusion in the turbulent tetrad dispersion. Two typical kinds of particles, i.e., the gas-phase molecules (Schmidt number, Sc = 1) and marked fluid particles (approximately to the soot particle, Sc = 1), have been considered. The effects of differential diffusion are found to be closely related to the tetrad initial size and weakly dependent on the Taylor-scale Reynolds number. Although the discrepancy of the tetrad size and shape between the gas-phase molecules and marked fluid particles is more significant in the turbulent dissipation scale, this influence can extend to larger scales, i.e., the inertial range, due to the correlated fluctuations and intermittency of the turbulent flows. It indicates that the molecular diffusion as well as the differential diffusion cannot be neglected if interested processes belong to small scales, such as in the combustion. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Scalar dispersion and mixing in turbulent flows are closely related to many types of industrial and environmental processes, such as the atmospheric pollutant dispersion and combustion. Generally, turbulent scalar dispersion and mixing processes depend on both turbulent motions in the macro-scale and molecular transport in the micro-scale. For physical or chemical phenomena that take place at small temporal or spatial scales, molecular transport becomes more important compared to turbulent motions. An excellent example is the combustion, in which chemical reactions happen in the molecular collision scale and strongly interact with the molecular and turbulent transport [1]. As well known, the molecular transport, which determines the species mixing in combustions, can influence the flame instability and structure [2], while their differences among species can further affect the dynamic process of combustion and pollutant emissions [3–6]. ⇑ Corresponding author. E-mail address: [email protected] (Z. Liu). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.07.002 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.

Since species mixing as well as chemical reactions take place at the molecular level, it is very beneficial to study and model these phenomena from a kinetic viewpoint [7–9]. In recent years, several efforts have been made to numerically investigate combustion phenomena by using the direct simulation Monte Carlo (DSMC) method, such as for the spontaneous combustion [10,11] and the unstretched diffusion flame [12] of hydrogen oxygen mixture. Detailed reaction mechanism has been confirmed in these simulations, however the transport-chemistry interaction is rarely considered because of the limited computational size and time step of the DSMC method. Therefore, applying the DSMC method to the multi-scale problems such as combustions is still challenging [13]. Simulating molecular motions following the Fokker–Planck model instead of the Boltzmann equation [14–17] can break these limitations in the DSMC method. Since a continuous stochastic process replaces the binary molecular collision compared to the DSMC method, the numerical methods in the framework of the Fokker– Planck model are less restricted by the temporal and spatial scales and more efficient for small or moderate Knudsen number (Kn) flows. Therefore, they are feasible for simulating multi-scale

88

F. Fei et al. / International Journal of Heat and Mass Transfer 103 (2016) 87–98

problems. For example, the Fokker–Planck model has been verified both in the rarefied [15,16] and continuum gas flows [17]. In the continuum regime, the Fokker–Planck model would degenerate to a diffusive form, based on which, the diffusive information preservation (D-IP) method has been proposed [17]. In the previous papers, the D-IP method has shown advantages in simulating small Knudsen number gas flows [17], and has been further applied to the homogeneous isotropic turbulent (HIT) [18] and turbulent channel flows [19] at low Reynolds numbers. Since the D-IP method simulates the turbulent flows and molecular motions from the kinetic viewpoint uniformly, their interactions, which determine the turbulent dispersion and mixing processes, can be directly computed and sampled. In direct numerical simulations (DNS), the turbulent dispersion and mixing processes are often investigated by the Eulerian approach [51,52]. However, as transport issues are addressed naturally from the Lagrangian viewpoint, the Lagrangian approach [53–55] has also been used by many authors dated back to the work of G. I. Taylor (1921) [20]. Afterwards, the molecular diffusion effects on turbulent dispersion processes have been studied starting from the pioneer work of Saffman (1960) [21], such as for the scalar dispersion in the isotropic turbulence [22] and the turbulent channel flows [23–25], the concentration fluctuation of pollutant [26,27] and the chemical reactions [28–30]. Most of the above literatures focus on the Lagrangian behavior of a single particle or particle pairs, which are closely connected to the average expansion and concentration fluctuations of a cloud of passive contaminants, respectively. Nevertheless, the geometric evolution of the particle cluster is also important to understand the physical mechanism of dispersion and mixing processes in turbulent flows, such as their statistical conservation laws [31]. The geometric evolution of a cloud of Lagrangian particles is preferable to be modeled through the multi-particle dispersion [32], where the study of tetrad (four-particle cluster) that is the minimum configuration to define a volume, is of special interest. Although many authors [33–37] have investigated the evolution properties of tetrads expansion and distortion in turbulence through DNS and stochastic models by tracking marked fluid particles, the effects of molecular diffusion on tetrad dispersion have been rarely discussed, which would be significant in the small-scale processes, such as chemical reactions and combustions [2–6]. Thus, in order to study the molecular diffusion effects on the tetrad dispersion by using a kinetic numerical scheme, we would first review and verify the Fokker–Planck model and the D-IP method in Sections 2 and 3 in this paper. Then, the D-IP method is applied to compute the temporal and spatial evolutions of the tetrad dispersion and their Lagrangian statistics in Section 4. As mentioned by Pitsch and Bisetti [38,39], the differential diffusion among soot particles (Schmidt number, Sc  1) and gas-phase species (Sc  1) can cause soot to be transported towards the flame sheet where it is oxidized, which would influence the soot formation and evolution significantly. In analogy to the soot transport in combustions, two typical kinds of particles, i.e. the gas-phase molecules and the marked fluid particles, are considered in Section 4. In addition, their differential diffusion effects are going to be analyzed. Finally, conclusions and remarks of this paper are presented in Section 5.

2. Diffusive information preservation method

processes. When the external force is neglected, the Fokker–Planck equation can be written as,


 @f ðc; x; tÞ @f ðc; x; tÞ @ @ 2 fkB Tðx; tÞ þ ci f ; ¼ ½fðci  U i ðx; tÞÞf  þ @t @xi @ci m @ci @ci ð1Þ where ci is the molecular velocity, Ui(x, t) and T(x, t) are the mean velocity and temperature of the flow fields at location x and time t; kB and m are Boltzmann constant and molecular mass, respectively; f is the friction coefficient which characterizes the time scale of the relaxation process for the velocity distribution function f(c, x, t). In equilibrium, f(c, x, t) relaxes to corresponding Maxwell distribution fM,

   3=2 m m fM ¼ n exp  ðci  U i Þ2 ; 2pkB T 2kB T

ð2Þ

where n is the number density. In numerical simulations [14–17], the Fokker–Planck Eq. (1) is usually solved through the Langevin model by simulated particles (hereafter, referred to as simply ‘‘simulator”), i.e.,

dX i ¼ ci ; dt

ð3aÞ

  dci 2kB T eq dW i ðtÞ ; Þ þ ¼ fðci  ueq i dt dt m

ð3bÞ

where X is the position of the simulator; W(t) is a Wiener process with dWi(t) = Wi(t + dt)  Wi(t), hdWi(t)i = 0 and hdWi(t)dWj(t)i = dtdij, where dij is the Kronecker delta. The symbol h  i denotes an ensemble mean. Considering the location of simulators, we define ðX; tÞ and T eq ðX; tÞ as the mean velocity the Lagrangian quantities ueq k;i k and temperature of simulator k at the location X and time t. They can be considered as the local valuables corresponding to the Eulerian quantities at field position x, such as Ui(x, t) and T(x, t), 2 eq2 /2 + cvTeqi in a computational cell, and U i ¼ hueq i i; U /2 + cvT = hu where cv is the specific heat at constant volume. Additionally, we 00 defined c0i ¼ ci  ueq i and c i ¼ c i  U i as the thermal velocity correeq sponding to ui and Ui, respectively, where hc0i i ¼ hc00i i ¼ 0. In Refs. [14–17] the molecular motions are calculated based on explicit solutions of the Langevin model (3), where the corresponding mean quantities of simulators are assumed to be constant during the molecular motion in a time step Dt, i.e., eq eq eq 0 0 ueq k;i ðX; t 0 þ t Þ ¼ uk;i ðX0 ; t 0 Þ and T k ðX; t 0 þ t Þ ¼ T k ðX0 ; t 0 Þ;

ð4Þ

where t0 e [0, Dt), X0 = X(t0) is the simulator’s location at the initial time t0.Combining Eqs. (3) and (4), the molecular motion of the simulator k can be solved explicitly, D DX k;i ¼ X tk;i0 þDt  X tk;i0 ¼ ueq k;i ðX0 ; t 0 ÞDt þ DX k;i ;

ð5Þ

where DX Dk;i is the thermal motion. Afterwards, by using the integral solutions of the Langevin model [14,17,40], the thermal motion and molecular velocities of simulators can be sampled as

DX Dk;i ¼ c0k;i ðt 0 Þ

1  efDt ~ 1;i ; þ F 1=2 w f

ck;i ðt 0 þ DtÞ 

ueq k;i ðt 0 Þ

¼

c0k;i ðt0 ÞefDt

ð6aÞ 2 H ~ 1;i þ G  H þ F 1=2 w F F

!1=2 ~ 2;i ; w ð6bÞ

2.1. The linear Fokker–Planck model with The Fokker–Planck equation can be considered as an approximation of the Boltzmann equation for gas flows [56], where molecular binary collisions are replaced by continuous random

F¼

kB T eq mf2

ð4efDt  e2fDt þ 2fDt  3Þ;

ð7aÞ

89

F. Fei et al. / International Journal of Heat and Mass Transfer 103 (2016) 87–98

H¼

kB T eq 2 ð1  efDt Þ ; mf

ð7bÞ

G¼

kB T eq ð1  e2fDt Þ: m

ð7cÞ

~ 1;i and w ~ 2;i are independent, normal distributed random variables. w Finally, the mean velocity and temperature of the flow filed at t0 + Dt can be calculated by ensemble averaging from the molecular velocities, i.e., Ui = hcii and T ¼ hc002 i=2cv , respectively. Since the Fokker–Planck equation degenerates to a linear form under the assumption of Eq. (4), we refer to the particle numerical algorithm following Eqs. (4)–(7) as the linear Fokker–Planck (LFP) model. According to the LFP model, the time evolution of the velocity distribution function f of simulators can be divided into two parts: the equilibrium part and the remained non-equilibrium part,

f ðc; X; t 0 þ t 0 Þ ¼ f M ðc; X; t 0 Þ þ f RE ðc; X; t 0 þ t 0 Þ:

ð8Þ

The first term in the right-hand side of the Eq. (8) indicates that the equilibrium part for the velocity distribution function of simulator is unchanged during a time step according to the constant mean quantities as assumed above, while the remained nonequilibrium part would relax to the equilibrium state following the right-hand side of the Fokker–Plank Eq. (1). Furthermore, the conservation equations of density, momentum and energy for the LFP model can be derived from the moment equations of Eq. (1),

@  @ ðnQ Þ þ ðnci Q Þ ¼ D½Q ; @t @xi

ð9Þ

R where Q = m, mc, mc2/2, Q ¼ Qfdc=n, and D[Q] is the ensemble average of Q for the right-hand side of Eq. (1). Since density, momentum and energy are conserved in the Fokker–Planck equation [14], D[Q] = 0. According to Eqs. (1), (8) and (9), these conservation equations can be written as,

@ q @ qU i ¼ 0; þ @t @xi

ð10aÞ

@ @ @ @ ðqU j Þ þ ðqU i U j Þ þ ðqc0i cj jM Þ þ ðqc0i cj jRE Þ ¼ 0; @t @xi @xi @xi

ð10bÞ

@ @ @ @ ð2cv qT þ qU 2 Þ þ ½qU i ð2cv T þ U 2 Þ þ ðqc0i c2 jM Þ þ ðqc0i c2 jRE Þ ¼ 0; @t @xi @xi @xi ð10cÞ

where subscripts M and RE represent the ensemble average of fM and fRE, respectively. The first and second terms in Eq. (10) are time derivative and convective terms, respectively. The last two terms in the left-hand side of Eqs. (10b) and (10c) represent the momentum and energy transport due to the equilibrium and non-equilibrium parts of the velocity distribution function by thermal velocity c0i ; respectively. When the time step in the LFP model is much less than the molecular mean collision time sc, the equilibrium distribution fM(c, X; t0) in Eq. (8) approximates to fM(c, X; t0 + Dt), thus the last two terms in the left-hand side of Eqs. (10b) and (10c) can be rewritten as,

@ @ @ ðqc0i cj jM Þ þ ðqc0i cj jRE Þ ¼ ðqc00i c00j Þ; @xi @xi @xi @ @ @ ðqc0i c2 jM Þ þ ðqc0i c2 jRE Þ ¼ ðqc00i c2 Þ: @xi @xi @xi

ð11aÞ

ð11bÞ

Therefore, Eq. (11) is the same as the governing equation of the FPM method [14], which is appropriate for the gas flows at small temporal scales.

In contrast, when Dt > > sc, the non-equilibrium part in Eqs. (10b) and (10c) vanished, and the equilibrium part can be derived as,

@ @ @ @ ðqc0i cj jM Þ þ ðqc0i cj jRE Þ ¼ ðqc0i ueq ðqc0i c0j jM Þ; j jM Þ þ @xi @xi @xi @xi

ð12aÞ

@ @ @ ðqc0i c2 jM Þ þ ðqc0i c2 jRE Þ ¼ ðqc0i ð2cv T eq þ ueq2 ÞjM Þ @xi @xi @xi @ þ ðqc0i ðc2  2cv T eq þ ueq2 ÞjM Þ: ð12bÞ @xi We note that the thermal transport terms in Eq. (12) are identical to those in the D-IP method [17], which degenerate to diffusive form at large temporal scales. Given all this, the LFP model can both capture the kinetic behavior in the collision limit (Dt < < sc) as the FPM method and the diffusive behavior in the diffusion limit (Dt > > sc) as the D-IP method. Moreover, in the intermediate temporal scale, these two methods can be unified in the framework of the LFP model smoothly, which will be shown in our recent work [57]. All of these indicate that the LFP model has advantages in the simulation of multi-scale flows from the rarefied to continuum regimes. 2.2. Governing equation of the D-IP method In the current paper, one of our aims is to validate the capability of the LFP model to simulate gas flows at large scales, which is the fundament of further applications to study the multi-scale processes, such as the interaction between the transport phenomena and chemical reactions in the turbulent flows. Since the temporal scale in turbulent flows is much larger than the molecular mean collision time, the LFP model is consistent with the D-IP method in this situation as mentioned above. On one hand, the D-IP method tracks the molecular motions following the Langevin model without limitations on the time step and cell size. On the other hand, combined with the basic idea of the IP method [41–45], statistic noise of the D-IP method can be reduced significantly as same as the IP method, which is requisite for simulating the unsteady flows, such as the turbulent flows. Briefly, mean quantities of flow fields are calculated by sampling the simulators’ velocity ueq and temperature Teq (which are identical to the IP velocity and temperature in Ref. [17], i.e., ueq = uIP and Teq = TIP) instead of the molecular velocity itself in the D-IP method. The governing equation of the D-IP method for IP velocity can be derived from Eqs. (10b) and (12a). In Eq. (12a), the nonequilibrium part of the velocity distribution function vanished. Therefore, the last term of Eq. (12a) can be further simplified and directly calculated by averaging from the Maxwell distribution instead of the simulators’ thermal velocities [43]. In the D-IP method, this term represents the effect of pressure, i.e., the third term in the right-hand side of Eq. (13). Other terms in Eqs. (10b) are the time derivative and transport terms related to the IP velocity, respectively. Hence, by applying Eq. (10) to individual simulator as in Ref. [17], the governing equation of the D-IP method for IP velocity can be written as, @  IP  @  IP IP  @  0 IP  quj þ q ui uj þ qci uj @t @xi @xi |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflffl} collectiv e motion

thermal motion

@  0 IP  1 @  0 IP  ¼ ðwg  1Þ qci uj  wg q c i ui  @xi 3 @xj |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} v iscosity modified

@p @xj |{z}

þqaFj ðxk ; tÞ

pressure modified

ð13Þ eq

IP

where u has been replaced by u , and the detailed derivation can be found in Ref. [17]. The last two terms in the left-hand side of

90

F. Fei et al. / International Journal of Heat and Mass Transfer 103 (2016) 87–98

Eq. (13) correspond to the transport of the IP velocity due to molecular collective and thermal motions, respectively. Moreover, in order to satisfy the compressive Navier–Stokes (NS) equations in the continuum regime, additional viscosity-modified terms, which are not contained in the LFP model, should be considered in the D-IP method, i.e., the first two terms in the right-hand side of Eq. (13). The modified coefficient wg = l/(qD⁄), where D⁄ = D[1  (1  e-fDt)/fDt], l and D refers to the gaseous viscosity and diffusion coefficients, respectively [17]. The third term refers to pressure-modified term as mentioned above, while the last term of Eq. (13) represents the forcing acceleration term. The last force term is zero for the decaying HIT flows and non-zero for the forced HIT flows, which will be further discussed in the following section. In addition, to compare with the pseudo-spectral method [46] and reference data in the incompressible condition, the Mach number of the r.m.s. fluctuation velocity for the HIT flows is less than 0.1 in the current paper. Therefore, the energy equation as Eq. (10c) is not considered, and the computational domain of the turbulent flows is assumed to be contacted with a heat-reservoir in following simulations. 2.3. Numerical implementation The D-IP method can be divided into the following six steps to calculate the HIT flows in the current study. 2.3.1. Initialization Each simulator k is initialized by given a position Xk,i (k = 1,2,. . ., N), a thermal velocity c0k;i and IP velocity uIP k;i according to the initial gaseous flow field. In the HIT simulation, the initial flow field is a fully developed HIT flow field obtained from the pseudo-spectral method, whose velocity field is stored in the wave number space. To initialize the IP velocity of simulators efficiently, these simulators are initially placed in uniform nodes in each computational cell. Then, the Fast-Fourier-Transform (FFT) method is used to transform the initial velocity field from wave number space to physical space. The initial IP velocity of each simulator equals to the flow velocity at related nodes. 2.3.2. Molecular motion During a time step, simulator k acquires a new position and thermal velocity according to the LFP model [Eqs. (4)–(7)]. The friction coefficient f = kBT/mD is determined by the diffusion coefficient. 2.3.3. Interact with boundaries In the HIT simulation, the periodic boundary condition is used. 2.3.4. Change IP velocity due to collision effect The IP velocity of simulator k equals to the average value of all simulators at the same location after a time step, i.e., Dt uIP;tþ ðXÞ ¼ huIP;t ðXÞi . k

2.3.5. Update IP velocity After molecular motions, the IP velocity of simulator k is updated according to the governing equation Eq. (13). In practice, the pressure and viscosity modified terms can be calculated by the finite volume scheme. For example, in a control volume, the moment flux due to the thermal motions, such as qc0i uIP j in Eq. (13), could be sampled from the interface flux of computational cells during a time step,

2

qc0i uIPj ¼ Np 4

Nþ i

X

N i

muIP k1;j 

k1¼1

X k2¼1

3,

5 muIP k2;j

½DtS;

ð14Þ

where Np is the number of real molecules represented by a simula tor; and Nþ i ; N i are the numbers of simulators across the cell interface during a time step from the inlet and outlet directions, respectively, and S is the surface area of the cell interface. Additionally, the pressure p = nkBT could be obtained from the state equation. The forcing acceleration as the last term in Eq. (13) should be computed in every time step in the simulation of the forced HIT flows. The stochastic forcing scheme proposed by Eswaran and Pope [47] is used here. In their scheme, the random forcing acceleration aFj ðk; tÞ adds to the low wave numbers of the velocity field in wave number space k to maintain HIT flows. Nevertheless, the D-IP method is solved in physical space through simulators. The forcing acceleration should be transformed to the physical space by the FFT method first, then interpolated into the location of simulators, i.e., aFj ðxk ; tÞ. 2.3.6. Sample and average Finally, the mean velocity of the flow field is obtained by sampling the IP velocity of simulators in each cell, U i ¼ huIP i i. The density is calculated by using the mass conservation Eq. (10a). The turbulent kinetic energy ET and energy dissipation rate e of the HIT flows can be directly obtained from their definitions in physics, which are respectively given by

ET ¼

e¼

N 3 02 1 X 1 IP2 u ¼ u ; 2 N k¼1 2 k

ð15Þ

dET ET ðt0 þ DtÞ  ET ðt 0 Þ ¼ ; dt Dt

ð16Þ

where N is the number of simulators in the computational domain, and u0 is the r.m.s. fluctuating velocity of each component. The other turbulent characters used in the current paper can then be calcupffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi lated from ET and e, such as the Taylor microscale kT ¼ 15mu02 =e; m = l/q; the Taylor-scale Reynolds number Rek ¼ kT u0 =m; the Kolmogorov microscale g = (m3/e)1/4; and the Kolmogorov time-scale tg = (m/e)1/2 . 3. Validation with homogeneous isotropic turbulence In the previous paper [18], the D-IP method has been validated by the decaying and forced HIT flows at low Reynolds number (Rek ¼ 38), where their decaying rate of turbulent kinetic energy and energy dissipation rate, the Kolmogorov-scaled dissipation spectra and the iso-surface of vorticity are in well agreement with the results of the pseudo-spectral method. In the current paper, the D-IP method is used to simulate turbulent flows at higher Reynolds numbers. The flow medium is argon gas in standard condition (1 atm and 273 K), where l = 2.117  105 (Pa s), and q = 1.78 (kg/m3). The Schmidt number (Sc) is assumed to be 1.0, and other computational parameters for the decaying HIT flows are given in Table 1. Npc is the initial number of simulators in each grid, and the grid number in each direction is Ncell. The length of the computation domain is L, and the dimensionless value is 2p, where the character length Lc = L/2p. In addition, the velocity is normalized by the r.m.s. fluctuating velocity u0 . e0 ¼ e0 Lc =u03 is the nondimensional quantity of the initial turbulent energy dissipation

Table 1 Computational parameters for decaying HIT flows. Case

Ncell

Npc

Re0k

t⁄

e0

k0g0

kmaxg0

k0l0

TE0/tg0

1 2 3

64 128 256

512 343 216

38 86 130

0.0170 0.0040 0.0018

0.675 0.506 0.484

0.0500 0.0186 0.0104

1.51 1.12 1.26

1.18 1.00 1.00

7.440 11.25 16.29

91

F. Fei et al. / International Journal of Heat and Mass Transfer 103 (2016) 87–98

time steps to adjust to the initial incompressible flow fields. Therefore, as seen in Fig. 1(b), there are small fluctuations of the energy dissipation rate at the beginning of the simulation of the D-IP method. Since the Mach number is small and the energy equation has not been considered in the current simulations as mentioned in Section 2, the influence of the inconsistent initial condition would be negligible as shown in Fig. 1. However, for compressible turbulent flows, where the compressibility effects are significant, consistent initial conditions [58] should be employed. Additionally, three times are marked for each case in Fig. 1(c), i.e., (t1)–(t3), and their Kolmogorov-scaled dissipation spectra are represented in Fig. 1 (d), respectively. All results of the D-IP method are close to those obtained via the pseudo-spectral method with the same grid number, which demonstrates the validity of the D-IP method for the turbulent simulation. Moreover, to reduce the computational cost, small physical domains are applied for the D-IP method in this study, such as L = 3.25  104 m, 6.5  104 m and 1.3  103 m for cases 1–3, respectively. The cell size and time step in the HIT simulations are 80k and 25sc for all cases, where k and sc are the molecular mean free path and mean collision time. For the most computational cost case—case 3, the D-IP method takes 17.4 h by using 1024 CPU cores of Intel Xeon E5 in Tianhe-2 supercomputer center, for which the computational time is mostly used to track the molecular motions of 3.6  109 simulators. However, for traditional molecule-based method, such as the DSMC method, the computational costs mentioned above would be too huge to be afforded. After validated the D-IP method with the decaying HIT flows, two-particle relative dispersion of marked fluid particles is

rate e0, and Re0k is the corresponding Taylor-scale Reynolds number. t ¼ kT =ðRek Lc Þ is the non-dimensional kinematic viscosity. l0 and g0 are the initial integral scale and Kolmogorov microscale, respecpffiffiffi tively, while k0 = 2p/L and kmax ¼ 2k0 N cell =3 are the lowest and highest wave number. The initial turbulent eddy-turnover time TE0 = l0/u0 . The initial flow field for the decaying HIT is a fully developed HIT flow, which is evaluated by the pseudo-spectral method with stochastic forcing scheme proposed by Eswaran and Pope [47]. Corresponding to Table 1, the computational parameters for the stochastic forcing scheme are presented in Table 2. In analogy to the definition in Ref. [47], KF is the maximum wavenumber of the forced modes; eF and T F are the normalized input parameters of the energy dissipation rate and forcing time-scale, respectively (defined as Eqs. (9) and (12) in Ref. [47]). Decaying from the fully developed HIT flows, Fig. 1(a)–(c) show the temporal evolution of the turbulent kinetic energy, energy dissipation rate and Taylor-scale Reynolds number, respectively. Similar to e0, E0T is defined as the initial turbulent kinetic energy. As naturally a compressive scheme, the D-IP method needs several Table 2 Computational parameters for the stochastic forcing scheme of forced HIT.

1

38

2

86

3

130

KF pffiffiffi 2 2 pffiffiffi 2 2 pffiffiffi 2 2

eF

T F

0.0045

1.21

0.0024

0.75

0.0022

0.77

1.0

1.0

0.8

0.8

0.6

0.6

ε/ε0

Rek

E T/E T0

Case

PS (1) PS (2) PS (3) D-IP (1) D-IP (2) D-IP (3)

0.4

0.2

0.0 -1 10

0.2

PS PS PS D-IP D-IP D-IP

0.0 -1 10

100

0.4

101

100

102

t/tη0

PS (1) PS (2) PS (3) D-IP (1) D-IP (2) D-IP (3)

80

(ηk) E(k)/(εν )

Reλ

100

10

0

5 1/4

120

t3

10

-1

2

60

t2

40

100

101

t/tη0

(c)

PS (1) PS (2) PS (3) D-IP (1) D-IP (2) D-IP (3) 10

t1 10-1

102

(b)

140

0 -2 10

101

t/tη0

(a)

20

(1) (2) (3) (1) (2) (3)

102

103

-2

0.2

ηk

0.4

0.6

0.8 1 1.21.4

(d)

Fig. 1. The comparison of the D-IP method (lines) with the pseudo-spectral method (symbols) on (a) the kinetic energy, (b) the dissipation rate, (c) the Taylor-scale Reynolds number and (d) the Kolmogorov-scaled dissipation spectra at marked times t1–t3 for three cases in Table 1, respectively.

92

F. Fei et al. / International Journal of Heat and Mass Transfer 103 (2016) 87–98

calculated and compared with the results in Ref. [48] at Rek 38 before investigating the tetrad dispersion. Forced HIT flow evaluated following the case 1 in Table 2 is simulated to sample the Lagrangian statistics. Previously, the D-IP method with stochastic forcing scheme has been examined in [18]. As a molecule based numerical method, the Lagrangian statistics of the molecular motions can be directly sampled from the simulators in the D-IP method. Meanwhile, marked fluid particles are extra tracer particles arranged in a flow field without influencing turbulent flow field; their velocities u equal to the IP velocities of simulators at the same location, i.e., u(X) = huIP(X)i. In addition, three initial separations r0 between particle pairs around the Kolmogorov microscale are considered, such as r0/g = 0.25, 1 and 4, respectively. In two-particle relative dispersion, relative differential velocity ðrÞ

Dij

ðrÞ

and displacement covariance F ij

are calculated, which are

defined as [48], ðrÞ

ðrÞ

ðrÞ

^i ðt 1 Þu ^ j ðt 2 Þi; Dij ðt 1 ; t2 Þ ¼ hu

ð17aÞ

ðrÞ ðrÞ ðrÞ F ij ðt1 ; t 2 Þ ¼ h^xi ðt1 Þ^xj ðt 2 Þi;

ð17bÞ

^ iðrÞ ¼ u ^ ð1Þ ^ ð2Þ ^ðrÞ ^ð1Þ  ^xð2Þ where u u i i , xi ¼ xi i , and the superscript notations (1) and (2) represent each particle in the pair. In addition, ^ i ¼ ui  ui;0 , ^xi ¼ X i  ui;0 t  X i;0 ; where ui,0, Xi,0 are the initial mean u ðrÞ

ðrÞ

velocity and location of particle, respectively. Dij and F ij represent the evolution of the relative velocity and displacement between two particles, and correspond to an inertial reference frame moving with the particle initial velocity. Moreover, a cumulative average [35,48] is used to reduce statistical noise arising from temporal variability of the space-averaged energy dissipation rate in the stochastic forcing scheme, which is given by

et1=3 ¼ t1

Z

t

e1=3 ðt0 Þdt0 :

ð18Þ

0 3

In the following part, e(t) would refer to ðe1=3 t Þ , and 1=3 3=2 , t Þ

1=3 t 0 ðtÞ ¼ r 2=3 . tg ðtÞ ¼ ðm1=3 =e 0 =et Fig. 2 shows the mean-square relative differential velocity and relative differential dispersion of Eq. (17), while the data of the structure function (Dii(t, t)/3et) and mean square differential dispersion (Fii(t, t)/et3) versus (t/tg) are also provided for a single particle. For all cases, the results of the D-IP method agree well with the reference data as seen in Fig. 2, which verifies the D-IP method for Lagrangian statistics in the turbulent simulation.

10

4. Molecular diffusion effect on the tetrad dispersion The tetrad dispersion includes more information than twoparticle dispersion in geometry, not only in terms of the tetrad size expansion but also on its shape distortion. In this section, we will analyze the effects of differential diffusion between the gasphase molecules and marked fluid particles on the temporal and spatial evolutions of turbulent tetrad dispersion. Tetrads are uniformly arranged in the computational domain initially, for which three particles are set along the coordinate axes separated from the vertex particle by a distance r0 as same as that in Ref. [35]. The number of tetrads in every case is equal to the number of computational cells. Similar to the two-particle dispersion, forced HIT flows evaluated from the cases 1–3 in Table 2 are simulated, while r0/g = 0.25, 1 and 4, respectively. In order to characterize the evolution of the tetrad size and shape, the moment of inertia tensor I = qqT is introduced [34,35], where q ¼ ½qð1Þ qð2Þ qð3Þ : The three column vectors of q are written as,

pffiffiffi

qð1Þ ¼ ðXð2Þ  Xð1Þ Þ= 2;

ð19aÞ

pffiffiffi ¼ ð2Xð3Þ  Xð2Þ  Xð1Þ Þ= 6;

qð2Þ

ð19bÞ pffiffiffiffiffiffi

qð3Þ ¼ ð3Xð4Þ  Xð3Þ  Xð2Þ  Xð1Þ Þ= 12;

ð19cÞ

(i)

where X represents the location of each particle in tetrad (i = 1,2. . .4). By assigning the eigenvalues of matrix I as gi (g1 > g2 > g3), the gyration radius of the tetrad can be obtained,

R2 ¼

3 X g a ¼ TrðIÞ:

ð20Þ

a¼1

The gyration radius of the initial tetrad is R20 ¼ 9r 20 =4. As the sum of the eigenvalues represents the size of the tetrad, their ratios, i.e.,

Ia ¼ g a =R2 ; ða ¼ 1; 2; 3Þ;

ð21Þ

characterize the shape of tetrad, e.g., I3 = 0 for a sheet-like tetrad and I2 = I3 = 0 for a needle-like tetrad. Another convenient dimensionless parameter to measure the tetrad shape, named as the shape parameter, is defined as

K ¼ V 2=3 =R2 ¼ 32=3 ðI1 I2 I3 Þ1=3 ;

ð22Þ

where V = (1/3)(g1g2g3)1/2 is the volume of the tetrad. The shape parameter K varies from 0 for sheet-like tetrad to 35/3 for regular

1

single particle

10

0

10

-1

10

-2

10

-3

(r)

Fii (t,t)/εt

3

D ii(r)(t,t)/εt

100

single particle

10-1

10

-2

10

-1

10

0

10

t/t0

(a)

1

10

2

10-1

100

t/t0

101

102

(b)

Fig. 2. The comparison of (a) the mean-square relative differential velocity and (b) relative differential dispersion for two-particle dispersion. The lines and symbols are results of the D-IP method and Ref. [48], respectively. The solid lines and corresponding symbols represent the data for particle pairs vs. (t/t0), from bottom to top, r0/g = 0.25, 1 and 4, respectively. The dashed line and rhombus symbols respect the structure function and mean square differential dispersion for the single particle vs. (t/tg).

93

F. Fei et al. / International Journal of Heat and Mass Transfer 103 (2016) 87–98

tetrahedral, and K = 0.1346 for the initial tetrad in the current paper. 4.1. Temporal evolution of the tetrad dispersion First, we consider the temporal evolution of the tetrad size, which is usually represented by the r.m.s. gyration radius of tetrad defined in Eq. (20). Fig. 3(a) shows their evolution scaled by Kolmogorov variables at Rek 38 for three initial separations, i.e., r0/g = 0.25, 1 and 4, respectively. The curves for the gas-phase molecules and marked fluid particles are significantly different at early times (t < tg), and collapse to the turbulent diffusive limit at large times. At early times, turbulent motions expand the gyration radius of tetrad as t scaling. When the tetrad initial size is in the dissipation range, the expression for the mean-squared gyration radius can be written as [35], 2

ðhR i 

hR20 iÞturbulence 2

g



 2 3 r 20 t ¼ : 4 g2 tg

ð23Þ

However, molecular diffusion expands the gyration radius of tetrad as t1/2 scaling, i.e.,

ðhR2 i  hR20 iÞdiffusion

g

2

¼

18Dtg

g

2

  t : tg

ð24Þ

As the Schmidt number Sc = l/qD  1 in gaseous flows, their molecular diffusivity D  g2/tg. Therefore, the time scale of molecular diffusion for the gas-phase molecule is characterized by tg in the current study, wherein the scalar diffusion range coincides with the viscous dissipation range. Furthermore, we note that the ratio between Eqs. (23) and (24) approximates to 24(g/r0)2(t/tg)1, which indicates that the influence of differential diffusion depends on the initial size of tetrad in the dissipation range. As seen in Fig. 3(b), when initial separation r0 is less than the Kolmogorov microscale, the effects of molecular diffusion are more significant than the turbulent dissipation scale motions at early times. Therefore, with molecular diffusion, the curves for initial size r0/g = 0.25 and 1 approximate to t1/2 scaling before tg in Fig. 3(a). Moreover, since molecular diffusion is independent on the separation distance of particles, the influence of the initial configuration will be quickly vanished when molecular diffusion dominates. Hence, curves for two small initial tetrads (r0/g = 0.25 and 1) approach each other at t  tg in Fig. 3(a). In contrast, without molecular diffusion, the

marked fluid particles in small initial tetrads are slowly separated by the turbulent dissipation scale motions at early times, where the separation rate of particle pairs is proportional to r0 in this situation. Thus, as the integral of relative velocity, the r.m.s. gyration radius for marked fluid particles has a long memory of the initial configuration as pointed out by Yeung [49]. Therefore, the discrepancy of the r.m.s. gyration radius between the gas-phase molecules and marked fluid particles at early times remains for a long period due to their differential diffusion as shown in Fig. 3(b). At large times (t TL), the turbulent diffusive limit of the meansquared gyration radius is also presented in Fig. 3(a), whose approximate solution [35] is,

hR2 i

g

2

¼

18Rk 15

1=2

   TL t ; tg tg

ð25Þ 3=4

where TL/tg is estimated by T L =tg ¼ ½4:77 þ ðRk =12:6Þ4=3  [48]; the Lagrangian integral time scale TL are nearly 5.2tg, 8.6tg, and 12tg at Rek 38, 86, and 130, respectively. In the intermediate regime between the dissipation and turbulent diffusive range, both the influences of molecular and turbulent diffusion become weak. Generally, when the dissipation and turbulent diffusive range are widely separated, the expansion of tetrad’s r.m.s. gyration radius would satisfy the inertial range t3/2 scaling as a dashed line of slope 3/2 in Fig. 3(a), i.e., inertial range scaling,

hR2 i

g2

¼

 3 3 t g ; 2 tg

ð26Þ

where g is Richardson constant. However, this scaling cannot be observed in Fig. 3(a) at such low Reynolds number. In Fig. 4, the temporal evolution of the r.m.s. gyration radius at higher Reynolds numbers, such as Rek ¼ 86 and 130, are calculated for initial separations r0/g = 1. In contrast to the dependence of tetrad initial size, the temporal evolution of the r.m.s. gyration radius scaled by Kolmogorov variables is independent on the Tayler scale Reynolds number at early dispersion times according to Eqs. (23) and (24). However, at large times, these curves would collapse to different turbulent diffusive limits of Eq. (25) in the diffusive range, which are functions of the Reynolds number. Moreover, corresponding data in Ref. [35] for r0/g = 1 at Rek 140 are plotted as reference in Fig. 4. They are identical with the results of Rek 38, 86 and 130 calculated via the D-IP method at early times for the independence of the Tayler scale Reynolds number.

103

1

marked fluid particle

2

t1/2

Reλ=38, r0=0.25η r0=η r0=4η

0.8

gas-phase molecule

Reλ=38, r0=0.25η r0=η r0=4η

0.6

ωR

1/2/η

102

Reλ=38, r0=0.25η r0=η r0=4η

101

0.4 3/2

t 10

0.2

0

t

1/2

0 10

-1

10

-1

10

0

t/tη

(a)

10

1

10

2

0

20

40

60

t/tη

80

100

120

(b)

Fig. 3. (a) Comparison the temporal evolution of the r.m.s. gyration radius of tetrad between the gas-phase molecules and marked fluid particles for r0/g = 0.25, 1 and 4 at Rek 38. The dashed line of slope 1.5 provides reference for inference of inertial range scaling. The solid line t1/2 at large times shows the turbulent diffusive limit, which is calculated from Eq. (25), while the solid line t1/2 at early times represents the molecular diffusive scaling. (b) The discrepancy of the r.m.s. gyration radius between the gash i 2 1=2 2 1=2 phase molecules and marked fluid particles at different r0, xR ¼ hR2 i1=2 gas molecule  hR ifluid particle =hR igas molecule .

94

F. Fei et al. / International Journal of Heat and Mass Transfer 103 (2016) 87–98

103 marked fluid particle

1/2/η

102

t1/2

r0=η, Reλ=38 Reλ=86 Reλ=130 gas-phase molecule

r0=η, Reλ=38 Reλ=86 Reλ=130

t3/2

r0=η, Reλ=140,Hackl.

101

100 10

-1

10

0

10

1

10

2

t/tη Fig. 4. Comparison the temporal evolution of the r.m.s. gyration radius of tetrad between the gas-phase molecules and marked fluid particles for r0/g = 1 at Rek 38, 86 and 130. The dashed line of slope 1.5 provides reference for inference of inertial regime scaling. The lines t1/2 at large time show turbulent diffusive limit at Rek 38 (black solid line) and 86 (red dash-dot line), respectively, which are calculated from Eq. (25). Symbols represent data of Ref. [35] at Rek 140. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Additionally, it is noted that there is little difference among the curves of these three Reynolds numbers at the very beginning in Fig. 4. We attribute it to the slight deviation of the initial separation from the requested value in the simulations; however, these errors do not influence the conclusions above. Another interesting phenomenon in tetrad dispersion process is the shape distortion. Fig. 5(a) shows the temporal evolution of the mean value of K for tetrad with initial size r 0 =g ¼ 0:25, 1 and 4 at Rek 38. For marked fluid particles without molecular diffusion, tetrads are strongly distorted from the initial value K = 0.1346 to a minimum value close to zero due to turbulent stretching motions in the dissipation range, which represent a sheet-like configuration. Afterwards, hKi recovers to the diffusive limit value, i.e., hKi = 0.0645, at large times (t TL), where the motions of tetrad

particles are independent according to the energy containing eddies. Given that, the Reynolds number is not sufficiently high; no tendency of hKi to the inertial range scaling value 0.045 can be captured in Fig. 5(a) for the marked fluid particles, such as an inflection point in Ref. [35]. However, the evolution of hKi is significantly different for the gas-phase molecules when molecular diffusion is considered. Firstly, in the dissipation range, since molecular diffusion separates particles in tetrad independently the same as the turbulent motions in the diffusive range, hKi approaches to the diffusive limit value instead of distorted to zero in the early dispersion. Secondly, as molecular diffusion dominates in the dissipation range as mentioned above, the turbulent fluctuations would not influence the tetrad distortion until the tetrad size expanding to a larger scale than the Kolmogorov microscale. Therefore, it seems that the tetrads are less affected by the turbulent motions in the dissipation range for the gas-phase molecules, whose mean value of K approximates to the inertial range scaling value in the intermediate dispersion time as shown in Fig. 5(a). Moreover, the discrepancy of hKi between the gas-phase molecules and marked fluid particles are plotted in Fig. 5(b). Similar to Fig. 3(b), the discrepancy is more significant for tetrad with small initial size. However, comparing Figs. 3(b) and 5(b), we can find that the discrepancy of the r.m.s. gyration radius vanished after a longer time than that of hKi, which indicates that the temporal evolution of the tetrad size has a longer memory of the initial configuration than that of the tetrad shape. Fig. 6 shows the temporal evolution of hKi at higher Reynolds numbers for r0/g = 1. Similar to the evolution of the r.m.s gyration radius, the dependence of the Tayler scale Reynolds number is not obvious for shape distortion in early dispersion both with and without molecular diffusion. However, as Reynolds number increasing, hKi is delayed to recover to the diffusive limit, which is attributed to the more separated dissipation and diffusive ranges. A more detailed description of the temporal evolution of tetrad shape can be observed from shape factors I1, I2 and I3, respectively. Fig. 7 shows the evolution of their mean values for r0/g = 1 at Rek 38, 86 and 130. In analogy to hKi, either molecular or turbulent diffusion leads shape factors approach to the diffusive limit value at the early and large times, i.e. hI1i = 0.75, hI2i = 0.22, and hI3i = 0.03, respectively. Compared to the marked fluid particles, the curves for the gas-phase molecules have an extreme point closer to the inertial regime scaling value similar to hKi in Figs. 5 and 6, i.e., hI1i = 0.825, hI2i = 0.16, and hI3i = 0.015 [35]. Moreover, Fig. 7 shows that hI1i hI2i > hI3i and hI3i  0 after the 0.05

0.15 marked fluid particle Reλ=38, r0=0.25η r0=η r0=4η gas-phase molecule Reλ=38, r0=0.25η r0=η r0=4η

0.03

ωΛ

<Λ>

0.10

Reλ=38, r0=0.25η r0=η r0=4η

0.04

0.02

0.01

0.05

0.00

0.00

10-2

10-1

100

101

t/t0

(a)

102

-0.01

0

20

40

60

80

100

120

t/tη

(b)

Fig. 5. (a) Comparison the temporal evolution of the mean value of K between the gas-phase molecules and marked fluid particles for r0/g = 0.25, 1 and 4 at Rek 38. The horizontal solid line of 0.0645 represents the turbulent diffusive limit; and the horizontal dashed line of 0.045 marks the inertial range scaling value. (b) The discrepancy of the mean value of K between the gas-phase molecules and marked fluid particles vs. t/tg for different r0, xK ¼ hKigas molecule  hKifluid particle .

95

F. Fei et al. / International Journal of Heat and Mass Transfer 103 (2016) 87–98

0.15

tetrad dispersion for a moment, which implies that tetrad distorted by HIT flows tends to form sheet-like and needle-like configurations regardless of whether molecular diffusion is considered.

marked fluid particle

r0=η,Reλ=38 Reλ=86 Reλ=130

4.2. Spatial evolution of the tetrad dispersion

gas-phase molecule

r0=η,Reλ=38 Reλ=86 Reλ=130

<Λ>

0.1

In addition to the temporal evolution, the molecular diffusion would also affect the spatial evolution of the tetrad dispersion, such as the distribution of the tetrad size, the expansion rate and the shape parameter versus tetrads’ mean gyration radius. The distribution of the tetrad size can be generally represented by the  2 i, where the mean variance of gyration radius, i.e., r2 ¼ hðR  RÞ  ¼ hRi. The mean square deviation r implies the gyration radius R

0.05

0

10

-1

10

0

10

t/t0

1

10

width of the distribution of tetrad gyration radius around its mean value hRi. Fig. 8(a) shows the discrepancy of r between the gas-phase molecules and marked fluid particles versus hRi/g for r0/g = 0.25, 1 and 4 at Rek 38. At the initial value of hRi/g, i.e., hRi  R0, r2 approaches to the right-hand side of Eqs. (23) and (24), which indicates that the variance of the gyration radius for marked fluid particles is less than that for the gas-phase molecules at the very beginning of the curves as seen in Fig. 8. As the tetrad expending, the marked fluid particle pairs in tetrad are separated exponentially in the dissipation range, which is sensitive to the intermittent nature of the turbulent stretching rate [50]. Thus,

2

Fig. 6. Comparison the temporal evolution of the mean value of K between the gasphase molecules and marked fluid particles for r0/g = 1 at Rek 38, 86 and 130, respectively. The horizontal solid line of 0.0645 represents the turbulent diffusive limit; and the horizontal dashed line of 0.045 marks the inertial range scaling value.

0.5

1

marked fluid particle

0.9

r0=η, Reλ=38 Reλ=86 Reλ=130

0.4

gas-phase molecule

0.8

r0=η, Reλ=38 Reλ=86 Reλ=130

0.7





0.3 marked fluid particle

r0=η, Reλ=38 Reλ=86 Reλ=130

0.6

0.2

gas-phase molecule

r0=η, Reλ=38 Reλ=86 Reλ=130

0.5

0.4

10

-1

10

0

10

0.1

1

10

0 -1 10

2

10

0

10

1

10

2

t/t0

t/t0

(a)

(b) 100 marked fluid particle

r0=η, Reλ=38 Reλ=86 Reλ=130 gas-phase molecule

10-1



r0=η, Reλ=38 Reλ=86 Reλ=130

10

-2

10

-3

10

-1

10

0

10

1

10

2

t/t0

(c) Fig. 7. Comparison the temporal evolution of the mean value of (a) I1, (b)I2 and (c) I3 for r0/g = 1 at Rek 38, 86 and 130, respectively. The horizontal solid line represents the turbulent diffusive limit; and the horizontal dashed line marks the inertial range scaling value.

96

F. Fei et al. / International Journal of Heat and Mass Transfer 103 (2016) 87–98 6.0

5.0

r0=η, Reλ=38 Reλ=86 Reλ=130 r0=4η, Reλ=38 Reλ=86 Reλ=130

3.0

Reλ=38, r0=0.25η r0=η r0=4η

4.0

ωσ

ωσ

2.0

3.0

2.0

1.0

1.0

0.0 100

0.0 100

101

102

101

102

/η

/η

(a)

(b)

Fig. 8. The ratio of the mean square deviation of gyration radius between the marked fluid particles and gas-phase molecules vs. hRi/g (a) for tetrad initial size r0/g = 0.25, 1 and 4 at Rek 38; and (b) for tetrad initial size r0/g = 1 and 4 at Rek 38, 86 and 130, respectively, where xr ¼ rfluid particle =rgas molecule .

of the initial configuration. When the tetrad mean gyration radius further expands to the diffusive range, the effects of molecular diffusion and initial configuration are both negligible; therefore, all of the curves in Fig. 8 decrease to 1.0 monotonously for large mean gyration radius. Additionally, we note that the curves for Rek 38 decrease to 1.0 faster than the other two higher Reynolds numbers when hRi > 100g in Fig. 8(b). We infer that it should result from the unseparated dissipation and diffusion range at this small Reynolds number. Moreover, the discrepancy of the variance of the gyration radius weakly depends on the Tayler scale Reynolds number as shown in Fig. 8(b). It indicates that the large-scale turbulent motions would not influence the differential diffusion, which closely related to the discrepancy between the molecular diffusion and turbulent motions in small scales. Furthermore, the spatial evolution of the tetrad size distribution is expected to strongly influence the expansion rate and shape distortion of tetrad. We compare the expansion rate of the mean gyration radius, i.e.,

tetrads in the regime with small stretching rate would remain compact for longer times, while tetrads in the regime with large stretching rate would separate faster in contrast. Therefore, widely distribution of gyration radius would be obtained in the dissipation range for the marked fluid particles. However, for the gas-phase molecules, molecular diffusion dominates molecular motions in the dissipation range, which leads a chi squared distribution of R2 alone. Therefore, we can find r for gas-phase molecules is smaller than that for the marked fluid particles at small mean gyration radius in Fig. 8(a). Moreover, as the initial separation of tetrad increasing, the effects of molecular diffusion become less significant, the discrepancy of the variance of the gyration radius tends to be vanished, such as for r0/g = 4. Fig. 8(b) shows the ratio of the mean square deviation of the gyration radius for r0/g = 1 and 4 at Rek 38, 86 and 130. Although the discrepancy of the distribution of tetrad size caused by their differential diffusion is significant in the dissipation range, it decreases slowly. For example, in Fig. 8(b), when r0/g = 1, though the ratio of r has a maximum value at about hRi = 7g in the dissipation range, it is still obvious for hRi in the inertial range (where the mean gyration radius is about hRi > 60g). We infer that it can mainly be attributed to the left tail of the widely distribution of the gyration radius for marked fluid particles in small scales, which has a long memory

v ¼ d½hRi=g=d½t=tg ;

ð27Þ

between the gas-phase molecules and marked fluid particles at the same mean gyration radius in Fig. 9(a) for r 0 =g ¼ 0:25; 1 and 4 at Rek 38. According to the competition between molecular diffusion 8

8 6

6

Reλ=38, r0=0.25η

r0=η, Reλ=38 Reλ=86 Reλ=130 r0=4η,Reλ=38

r0=η 4

4

r0=4η

ωχ

ωχ

Reλ=86 Reλ=130 2

2

1

1

10

0

10

1

/η

(a)

10

2

10

0

10

1

10

2

/η

(b)

Fig. 9. The ratio of the expansion rate of mean gyration radius between the gas-phase molecules and marked fluid particles vs. hRi/g (a) for tetrad initial size r0/g = 0.25, 1 and 4 at Rek 38; and (b) for tetrad initial size r0/g = 1 and 4 at Rek 38, 86 and 130, respectively, where xv ¼ vgas molecule =vfluid particle .

97

F. Fei et al. / International Journal of Heat and Mass Transfer 103 (2016) 87–98 0.06

0.05

Reλ=38, r0=0.25η r0=η r0=4η

r0=η, Reλ=38 Reλ=86

0.04

Reλ=130 r0=4η, Reλ=38

0.04

Reλ=86 Reλ=130

0.03

ωΛ

ωΛ

0.02

0.02

0.01 0

0

-0.01

-0.02 100

101

102

-0.02 100

/η

101

102

/η

(a)

(b)

Fig. 10. The discrepancy of the mean value of K between the gas-phase molecules and marked fluid particles vs. hRi/g (a) for tetrad initial size r0/g = 0.25, 1 and 4 at Rek 38; and (b) for tetrad initial size r0/g = 1 and 4 at Rek 38, 86 and 130, respectively, xK ¼ hKigas molecule  hKifluid particle .

and turbulent motions, the discrepancy of the expansion rates is more obvious for small mean gyration radius in the dissipation range and vanished for large mean gyration radius in the diffusive range. Because of the turbulent fluctuations are spatial correlated in small scales, in analogy to the discrepancy of the distribution of tetrad size in Fig. 8, the effects of differential diffusion on the tetrad expansion rate also have a long memory. Meanwhile, it is weakly dependent on the Tayler scale Reynolds number as seen in Fig. 9 (b) similar to the mean square deviation of the tetrad size distribution in Fig. 8(b). In addition, the same tendency can be observed from the shape distortion in Fig. 10, where the discrepancy of the mean value of K between the gas-phase molecules and marked fluid particles has been calculated. At the initial value of hRi/g, the molecular diffusion makes the shape parameter approach to the diffusive limit compared to the marked fluid particles as mentioned above, therefore, the discrepancy of hKi is negative. Then, the discrepancy of hKi becomes positive as the tetrads with marked fluid particle are more stretched away from their initial configuration by the turbulent motions, and this discrepancy would be more significant for smaller initial tetrad size. However, as seen in Fig. 10(b), the effects of differential diffusion are also weakly dependent on the Tayler scale Reynolds number for the tetrad distortion. Moreover, this discrepancy of hKi, which is raised from the dissipation range, extends to a large scale of the mean gyration radius as same as that in Figs. 8 and 9. Therefore, we conclude that the differential diffusion between the gas-phase molecules and marked fluid particles can affect the spatial evolution of the tetrad dispersion at a much larger range than the Kolmogorov microscale in gaseous turbulent flows. 5. Conclusion In the present paper, the D-IP method is used to compute the Lagrangian statistics of turbulent tetrad dispersion. Two typical kinds of particles in combustion are simulated in the current study: one respects the gas-phase molecules, Sc = 1; and the other one respects solid particles such as soot, Sc = 1. The effects of molecular diffusion on tetrad dispersion processes have been analyzed through these two species. The result shows that: (1) For the temporal evolution of the tetrad dispersion, compared to the turbulent motions, molecular diffusion significantly affects the turbulent dispersion at the small temporal scale (t < tg), which expands tetrad size faster

and distorts the tetrad shape less, thus makes the shape distortion scaling of tetrad closer to the inertial range value. Similarly, for the spatial evolution, the influence of the molecular diffusion is profound at the small spatial scale. For example, it results in a more narrowed distribution of tetrad size in the dissipation range. (2) The effects of differential diffusion between the gas-phase molecules and marked fluid particles should not be restricted in the dissipation scale. When the dispersion time and mean gyration radius beyond the dissipation range, the discrepancies of the tetrad size and shape as well as its size distribution remain for a long period due to the correlated fluctuations and intermittency of the turbulent flows, especially for tetrads with small initial separation (less than the Kolmogorov micro-scale). (3) Although closely related to the tetrad initial size, the effects of differential diffusion on tetrad dispersion are weakly dependent on the Tayler scale Reynolds number, which means that the differential diffusion is not influenced by the large-scale turbulent motions. Therefore, if interested processes belong to small initial scales, the effects of differential diffusion would be significant and should be considered at a structure larger than the Kolmogorov scale [6], such as for the species mixing and soot transport in combustion [39] and the pollutant released from small source [26]. Therefore, applying the D-IP method to investigate the mixture and dispersion processes with differential diffusion and chemistry reactions from the kinetic level may be an interesting subject in future study, especially when the non-equilibrium behavior of molecular transport is included.

Acknowledgments The supports of National Natural Science Foundation of China under Grant Nos. 51506063 & 51276076 and China Postdoctoral Science Foundation under No. 2014M562019 are gratefully acknowledged. This work is also supported by the Special Program for Applied Research on Super Computation of the NSFCGuangdong Joint Fund (the second phase). In addition, the authors are grateful to Prof. J. Fan, C.X. Xu and Dr. Y. Song for discussions and useful suggestions in the numerical simulation of the turbulent flows.

98

F. Fei et al. / International Journal of Heat and Mass Transfer 103 (2016) 87–98

References [1] Y. Ju, Recent progress and challenges in fundamental combustion research, Adv. Mech. 44 (2014) 201402. [2] A.N. Lipatnikov, J. Chomiak, Molecular transport effects on turbulent flame propagation and structure, Prog. Energy Combust. Sci. 31 (2005) 1–73. [3] R.W. Bilger, Molecular transport effects in turbulent diffusion flames at moderate Reynolds numbers, AIAA J. 20 (7) (1982) 962–970. [4] R. Hilbert, F. Tap, H. El-Rabii, D. Thevenin, Impact of detailed chemistry and transport models on turbulent combustion simulations, Prog. Energy Combust. Sci. 30 (2004) 61–117. [5] M.J. Dunn, R.S. Barlow, Effects of preferential transport and strain in bluff body stabilized lean and rich premixed CH4/air flames, Proc. Combust. Inst. 34 (1) (2013) 1411–1419. [6] C. Bruno, V. Sankaran, H. Kolla, et al., Impact of multi-component diffusion in turbulent combustion using direct numerical simulations, Combust. Flame 162 (11) (2015) 4313–4330. [7] S. Chapman, T.G. Cowling, The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases, Cambridge University Press, 1970. [8] G.A. Bird, Molecular Gas Dynamics and Direct Simulation of Gas Flows, Clarendon Press, Oxford, 1994. [9] D. Livescu, Numerical simulations of two-fluid turbulent mixing at large density ratios and applications to the Rayleigh-Taylor instability, Philos. Trans. R. Soc. Lond. A 371 (2003) (2013) 20120185. [10] G.A. Bird, Chemical reactions in DSMC, in: D.A. Levin, et al. (Eds.), 27th International Symposium on Rarefied Gas Dynamics, AIP Conference Proceedings 1333, 2010, pp. 1195–1202. [11] C. Yang, Q. Sun, Investigation of spontaneous combustion of hydrogen-oxygen mixture using DSMC simulation, in: J. Fan (Ed.), 29th International Symposium on Rarefied Gas Dynamics, AIP Conference Proceedings 1628, 2014, pp. 1261– 1267. [12] I.B. Sebastiao, A.A. Alexeenko, DSMC investigation of nonequilibrium effects in a H2–O2 unstretched diffusion flame, the 45th AIAA Thermophysics Conference, AIAA 2015–3372. [13] J. Fan, Rarefied gas dynamics: advances and applications, Adv. Mech. 43 (2) (2013) 185–201. [14] P. Jenny, M. Torrilhon, S. Heinz, A solution algorithm for the fluid dynamic equations based on a stochastic model for molecular motion, J. Comput. Phys. 229 (4) (2010) 1077–1098. [15] M.H. Gorji, M. Torrilhon, P. Jenny, Fokker–Planck model for computational studies of monatomic rarefied gas flows, J. Fluid Mech. 680 (2011) 574–601. [16] M.H. Gorji, P. Jenny, An efficient particle Fokker–Planck algorithm for rarefied gas flows, J. Comput. Phys. 262 (2014) 325–343. [17] F. Fei, J. Fan, A diffusive information preservation method for small Knudsen number flows, J. Comput. Phys. 243 (2013) 179–193. [18] F. Fei, J. Fan, et al. Statistical simulation of decaying and forced homogeneous isotropic turbulence, in: J. Fan (Ed.), 29th International Symposium on Rarefied Gas Dynamics, AIP Conference Proceedings 1628, 2014, pp. 1356–1362. [19] D.D. Zeng, F.Q. Zhong, J. Fan, Molecular simulation of 3D turbulent channel flow, in: J. Fan (Ed.), 29th International Symposium on Rarefied Gas Dynamics, AIP Conference Proceedings 1628, 2014, pp. 1349–1355. [20] G.I. Taylor, Diffusion by continuous movements, Proc. Lond. Math. Soc. A 20 (1921) 196–212. [21] P.G. Saffman, On the effect of the molecular diffusivity in turbulent diffusion, J. Fluid Mech. 8 (2) (1960) 273–283. [22] K. Kontomaris, T.J. Hanratty, Effect of molecular diffusivity on turbulent diffusion in isotropic turbulence, Int. J. Heat Mass Transfer 36 (5) (1993) 1403– 1412. [23] K. Kontomaris, T.J. Hanratty, Effect of molecular diffusivity on point source diffusion in the center of a numerically simulated turbulent channel flow, Int. J. Heat Mass Transfer 37 (13) (1994) 1817–1828. [24] Y. Mito, T.J. Hanratty, Lagrangian stochastic simulation of turbulent dispersion of heat markers in a channel flow, Int. J. Heat Mass Transfer 46 (6) (2003) 1063–1073. [25] A.K. Karna, D.V. Papavassiliou, Near-wall velocity structures that drive turbulent transport from a line source at the wall, Phys. Fluids 24 (3) (2012) 035102. [26] B.L. Sawford, J.C.R. Hunt, Effects of turbulence structure, molecular diffusion and source size on scalar fluctuations in homogeneous turbulence, J. Fluid Mech. 165 (1986) 373–400.

[27] M.S. Borgas, B.L. Sawford, Molecular diffusion and viscous effects on concentration statistics in grid turbulence, J. Fluid Mech. 324 (1996) 25–54. [28] B.M. Mitrovic, D.V. Papavassiliou, Effects of a first-order chemical reaction on turbulent mass transfer, Int. J. Heat Mass Transfer 47 (2004) 43–61. [29] K.T. Nguyen, D.V. Papavassiliou, Effects of a reacting channel on turbulent mass transfer, Int. J. Heat Mass Transfer 51 (2007) 2940–2949. [30] K.T. Nguyen, D.V. Papavassiliou, Flow effects on the kinetics of a second-order reaction, Chem. Eng. J. 140 (2008) 370–380. [31] G. Falkovich, K. Gawȩdzki, M. Vergassola, Particles and fields in fluid turbulence, Rev. Mod. Phys. 73 (4) (2001) 913–975. [32] B.L. Sawford, S.B. Pope, P.K. Yeung, Gaussian Lagrangian stochastic models for multi-particle dispersion, Phys. Fluids 25 (5) (2013) 055101. [33] A. Pumir, B.I. Shraiman, M. Chertkov, Geometry of Lagrangian dispersion in turbulence, Phys. Rev. Lett. 85 (25) (2000) 5324. [34] L. Biferale, G. Boffetta, A. Celani, B.J. Devenish, A. Lanotte, F. Toschi, Multiparticle dispersion in fully developed turbulence, Phys. Fluids 17 (11) (2005) 111701. [35] J.F. Hackl, P.K. Yeung, B.L. Sawford, Multi-particle and tetrad statistics in numerical simulations of turbulent relative dispersion, Phys. Fluids 23 (6) (2011) 065103. [36] B.J. Devenish, Geometrical properties of turbulent dispersion, Phys. Rev. Lett. 110 (6) (2013) 064504. [37] B.J. Devenish, D.J. Thomson, A Lagrangian stochastic model for tetrad dispersion, J. Turbul. 14 (3) (2013) 107–120. [38] H. Pitsch, E. Riesmeier, N. Peters, Unsteady flamelet modeling of soot formation in turbulent diffusion flames, Combust. Sci. Technol. 158 (1) (2000) 389–406. [39] F. Bisetti, G. Blanquart, M.E. Mueller, H. Pitsch, On the formation and early evolution of soot in turbulent nonpremixed flames, Combust. Flame 159 (1) (2012) 317–335. [40] S. Chandrasekhar, Stochastic problems in physics and astronomy, Rev. Mod. Phys. 15 (1) (1943) 1–89. [41] J. Fan, C. Shen, Statistical simulation of low-speed unidirectional flows in transition regime, in: R. Brun (Ed.), 21th International Symposium on Rarefied Gas Dynamics, Cepadus-Editions, Toulouse, 1999, p. 245. [42] J. Fan, C. Shen, Statistical simulation of low-speed rarefied gas flows, J. Comput. Phys. 167 (2) (2001) 393–412. [43] Q. Sun, I.D. Boyd, Theoretical development of the information preservation method for strongly nonequilibrium gas flows, AIAA Paper 4286 (2005) 2005. [44] N.D. Masters, W. Ye, Octant flux splitting information preservation DSMC method for thermally driven flows, J. Comput. Phys. 226 (2) (2007) 2044–2062. [45] J. Zhang, J. Fan, J. Jiang, Multiple temperature model for the information preservation method and its application to nonequilibrium gas flows, J. Comput. Phys. 230 (19) (2011) 7250–7265. [46] R.J. Rogallo, Experiments in Homogeneous Turbulence, Report TM-81315, NASA-Ames Research Center, 1981. [47] V. Eswaran, S.B. Pope, An examination of forcing in direct numerical simulations of turbulence, Comput. Fluids 16 (3) (1988) 257–278. [48] B.L. Sawford, P.K. Yeung, J.F. Hackl, Reynolds number dependence of relative dispersion statistics in isotropic turbulence, Phys. Fluids 20 (6) (2008) 065111. [49] P.K. Yeung, Direct numerical simulation of two-particle relative diffusion in isotropic turbulence, Phys. Fluids 6 (10) (1994) 3416–3428. [50] R. Scatamacchia, L. Biferale, F. Toschi, Extreme events in the dispersions of two neighboring particles under the influence of fluid turbulence, Phys. Rev. Lett. 109 (14) (2012) 144501. [51] F.A. Williams, Combustion Theory, second ed., Addison-Wesley, Redwood City, 1985. [52] P.K. Yeung, Correlations and conditional statistics in differential diffusion: scalars with uniform mean gradients, Phys. Fluids 10 (10) (1998) 2621–2635. [53] P.K. Yeung, S.B. Pope, Lagrangian statistics from direct numerical simulations of isotropic turbulence, J. Fluid Mech. 207 (1989) 531–586. [54] P.K. Yeung, Lagrangian investigations of turbulence, Annu. Rev. Fluid Mech. 34 (2002) 115–142. [55] F. Toschi, E. Bodenschatz, Lagrangian properties of particles in turbulence, Annu. Rev. Fluid Mech. 41 (2009) 375–404. [56] S. Heinz, Molecular to fluid dynamics: the consequences of stochastic molecular motion, Phys. Rev. E 70 (2004) 036308. [57] F. Fei, Z.H. Liu, C.G. Zheng, Unified discrete Fokker–Planck model for rarefied and continuum gas flows, in submitted. [58] J.R. Ristorcelli, G.A. Blaisdell, Consistent initial conditions for the DNS of compressible turbulence, Phys. Fluids 9 (1) (1997) 4–6.