A hybrid lattice Boltzmann – random walk method for heat transfer in gas–solids systems

Journal of Computational Physics: X 1 (2019) 100007

Contents lists available at ScienceDirect

Journal of Computational Physics: X www.elsevier.com/locate/jcpx

A hybrid lattice Boltzmann – random walk method for heat transfer in gas–solids systems Aaron M. Lattanzi a , Xiaolong Yin b , Christine M. Hrenya a,∗ a b

University of Colorado at Boulder, Dept. of Chemical and Biological Engineering, Boulder, CO, USA Colorado School of Mines, Petroleum Engineering, Golden, CO, USA

a r t i c l e Article history: Received 4 October Received in revised Accepted 5 January Available online 19 Keywords: DNS Heat transfer LBM Random walk Multiphase

i n f o 2017 form 30 November 2018 2019 February 2019

a b s t r a c t The development of accurate and robust heat transfer correlations for gas–solids flows is integral to the development of efficient multiphase unit operations. Direct numerical simulation (DNS) has been shown to be an effective method for developing such correlations. Specifically, the highly-resolved fields present in DNS may be averaged for use at the macroscopic level. Most commonly, particle-resolved immersed boundary or thermal lattice Boltzmann methods are employed. Here we develop a hybrid DNS framework where the hydrodynamics are resolved by the lattice Boltzmann method and the temperature field by random walk particle tracking (Brownian tracers). The random walk algorithm provides an efficient means for simulating scalar transport and can easily handle complex geometries. However, discontinuous fields pose a significant challenge to the random walk framework – e.g., a particle and fluid with different diffusivities. We derive a technique for handling discontinuities via the diffusivity, arising at a particle–fluid interface, and implement said method within the tracer algorithm. In addition, the heat transfer coefficient in the random walk method is defined and a technique for handling phases with different volumetric heat capacities is also developed. Moreover, the present algorithm is shown to correctly characterize intra-particle temperature gradients present in high Biot number systems. Verification of the code is completed against a host of cases: effective diffusivity of a static gas–solids mixture, hot sphere in unbounded diffusion, cooling sphere in unbounded diffusion, and uniform flow past a hot sphere. Predictions made by the new code are observed to agree with analytical solutions, numerical solutions, empirical correlations, and previous works. © 2019 Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction Fundamentally quantifying the transport phenomena within gas–solids flows often poses a significant challenge, owing to the presence of many different scales within the system (spatial and/or temporal). However, the transfer of thermal energy within a particle-laden flow is a key function of many industrial unit operations: fluidized beds [1–7], packed beds [8–13], spouted bed dryers [14,15], rotary kilns [16–18], and granular heat exchangers [19]; and thus, further theoretical advancements in this area are of great practical significance. Particle-resolved, direct numerical simulation (DNS) may be utilized to obtain high-fidelity predictions for the transport phenomena within a particle-laden flow, as the smallest scales present

*

Corresponding author. E-mail address: [email protected] (C.M. Hrenya).

https://doi.org/10.1016/j.jcpx.2019.100007 2590-0552/© 2019 Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

2

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

within the system are resolved. While computationally demanding, DNS has emerged as an effective tool for simulating non-isothermal gas–solids flows, as well as quantifying system’s heat transfer coefficients (HTCs) [20–26]. The types of DNS frameworks may be broadly grouped according to the type of numerical grid they employ. Namely, the domain may either be resolved by a body-conforming mesh or the system may be superimposed over a non-conforming grid. Boundary-conforming meshes often utilize integral methods for solution since the complex system geometry is resolved by an unstructured grid – i.e., the Navier–Stokes equations are approximated via finite element and finite volume techniques [27,28]. Strict enforcement of boundary conditions may be achieved by conforming the numerical grid to the system boundaries. However, if the boundaries are allowed to move, then the domain must be continuously re-meshed, creating a significant computational overhead. By contrast, this overhead is eliminated in methods that use non-conforming grids. For instance, immersed boundary methods (IBMs) do not conform to domain surfaces but instead introduce forcing terms at the boundary nodes to satisfy the desired boundary condition(s) [29–32]. Alternatively, the domain boundaries may be projected onto a fixed structured grid to generate a grid conforming approximation of the domain boundaries [33–35]. As the boundaries move, the approximation is updated. This update, however, incurs a much lower computational overhead than an unstructured grid. The lattice Boltzmann method (LBM) is a widely used method for solving the Navier–Stokes equations in particle-laden flows. The concepts and solution procedure for the LBM are founded in the kinetic theory, and as a result, they greatly differ from other frameworks such as finite element, finite difference, and finite volume. Specifically, LBM involves discretization of the continuous Boltzmann equation in a special manner that ensures the recovery of the Navier–Stokes hydrodynamics [36,37]. Thus, LBM involves discrete velocity distribution functions and not the hydrodynamic variables in the Navier–Stokes equations. For this reason, boundary conditions in the LBM framework are easily imposed directly upon the distribution functions. The ability of the LBM to efficiently simulate the Navier–Stokes equations (hydrodynamics) is now well documented. The incorporation of scalar transport (thermal energy or concentration) may be achieved in a variety of ways. Moment propagation is one method that may be utilized to simulate the advection-diffusion equation [38–40]. The moment propagation method utilizes an update scheme where portions of the scalar field are translated to adjacent LBM nodes according to the post-collisional velocity distribution. However, a Taylor expansion of the moment propagation method shows that the diffusivity will inherently be anisotropic and depend upon the fluid stress tensor (spatially varying) [38]. To minimize the spurious effects on diffusivity, the Peclet number must be kept small [39] or a modified moment propagation method must be applied [38]. In contrast to moment propagation, a secondary distribution function may be utilized for the scalar transport – i.e., lattice schemes [41,42] and thermal LBM [42–48]. The secondary distribution function is updated in an analogous manner to the velocity distribution function in classic LBM and the diffusivity is directly set through the relaxation parameter(s). In addition, a secondary framework is often coupled to the LBM [49–53]. Random walk particle tracking (RWPT) is one such framework that may be coupled to LBM to resolve scalar transport. In RWPT the position of passive Brownian tracers is tracked as they undergo movements that depend upon the velocity field (advection) and a stochastic term (diffusion). By contrast to the previous methods (grid-based), RWPT does not employ any discretization or numerical grid, and thus it is a meshless method which does not suffer from numerical dispersion in advection-dominated flows [54]. Furthermore, RWPT is capable of handling non-Fickian transport [55], complex geometries [56,57], and is straightforward to parallelize. In general, the scaling of computational costs for stochastic methods (RWPT) is quite favorable when compared to grid-based methods. As discussed in [38], for a spatial resolution of dx, the computational cost for stochastic methods tends to scale as dx−3 or dx−5/2 for the advection and diffusion dominated regimes, respectively. By contrast, grid-based methods scale as dx−4 or dx−5 . However, stochastic methods are subject to inherent statistical noise and the statistical noise scales inversely proportional to the square root of the tracer concentration [58,59]. Therefore, statistical methods, like the RWPT employed here, are computationally efficient for advection dominated flows, complex geometries, and when precision is required of a single value (e.g. heat transfer coefficient here); but to smoothly resolve the entire field, the tracer number required to reduce stochastic fluctuations will serve as a significant computational overhead. Here we expand the hybrid LBM–RWPT framework of [49,50] to allow for the simulation of heat (or mass) transfer in systems with a discontinuous diffusivity field. Previous work in this area has highlighted the significant challenges posed to the RWPT algorithm by discontinuous fields [56,60–62] – e.g., a particle and fluid with different thermal diffusivities. To address the issue of a discontinuous diffusivity field, we derive a method for altering the stochastic (diffusive) movement of tracers within the higher diffusivity medium such that the correct tracer distribution is achieved at equilibrium (homogeneous). Specifically, the modulus of the stochastic step is increased within the high diffusivity medium while only a subset of the tracers within the high diffusivity medium are allowed to undergo a diffusive step. The algorithm developed here is most similar to that of [60], though significant differences occur between the present work and that of [60]. Within the present work and that of [56,60–62], the thermal diffusivities associated with the two mediums (gas and solid-particles) are assumed to be constant (but not necessarily equal) and the generation of thermal energy due to shear work is considered negligible. In addition to addressing a discontinuous diffusivity field, the heat transfer coefficient in RWPT is defined in the present work and a technique for handling phases with different volumetric heat capacities is developed. The RWPT method developed here is shown to capture the heat transfer occurring in both low and high Biot number (Bi = hD p /k s ) systems. Thus, the present method may provide an efficient means for simulating heat transfer in large Biot systems since it does not suffer the computational overhead associated with carrying a highly resolved numerical grid for each particle. The new hybrid LBM–RWPT code is verified against a variety of systems: effective diffusivity of a static gas–solids mixture, hot sphere

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

3

in unbounded diffusion, cooling sphere in unbounded diffusion, and uniform flow past a hot sphere. Predictions made by the new code are found to be in agreement with analytical and numerical solutions, as well as previous works [21,47,63]. 2. Numerical techniques 2.1. Lattice-Boltzmann method (LBM) In this work the LBM is employed to resolve the hydrodynamic fields of fluid only – i.e., the density (ρ ) and velocity (u) are obtained from the LBM. The algorithm utilized within this work matches that of [33,34,64], and thus will only be briefly discussed here. Rather than discretizing the Navier–Stokes equations themselves, LBM involves the discretization of the continuous Boltzmann equation (BE) in time, space, and molecular velocity. While LBM is strongly tied to the BE, it should be noted that, unlike the BE, LBM cannot faithfully reproduce the flow occurring at high Knudsen numbers and thus is not valid for rarefied flows [65,66]. As a result of the LBM discretization, a collection of molecular velocity distributions (population densities) are stored at each node, for every time step. The discrete velocity distributions are updated through a streaming and collision process:





ni (r + ci t , t + t ) ≡ n∗i (r, t ) = ni (r, t ) + i n(r, t )

(1)

where ni is the discrete velocity distribution (population density of molecules at position r with a velocity ci ), t is the LBM time step (scaled to be unity), i is the collision operator, and n∗i is the post-collision distribution function. The collision operator in this work is linearized about the local equilibrium state (neq ) and matches that in [64]. The hydrodynamic quantities may then be found by taking the moments of the discrete distribution functions:

ρ=



j = ρu =

ni

i



=

n i ci

i



n i ci ci

(2)

i

where j is the momentum and  is the fluid stress tensor. Via a Taylor expansion of the finite difference approximation in Eq. (1) and completing a Chapman–Enskog expansion, it can be shown that the LBM recovers the Navier–Stokes equations with the following closures for the shear and bulk viscosities [64]:



η = −ρ c 2s

1

λ

+

1 2



ηb = −

2ρ c 2s 3



1

λb

+

1



2

(3)

where η is the shear viscosity, c 2s = 1/3 is the square of the speed of sound, λ is the eigenvalue of the collision matrix that corresponds to the relaxation of the off-diagonal portion of the non-equilibrium stress tensor, ηb is the bulk viscosity, and λb is the eigenvalue of the collision matrix that corresponds to the relaxation of the diagonal portion of the non-equilibrium stress tensor. Completion of a linear stability analysis shows that the eigenvalues in Eq. (3) are bounded by λ ∈ (−2 0) [64] – i.e., the viscosity is bounded and positive. Within the present work, the governing equations for the hydrodynamics (Navier– Stokes) and thermal energy (advection-diffusion) are essentially one-way coupled. Namely, the hydrodynamics will couple to the energy transport, but the temperature field will not couple back to the hydrodynamics – i.e., Eq. (1) does not contain a buoyancy forcing term from the temperature field (Bouissinesq approximation) and the fluid viscosity (Eq. (3)) will not be a function of the thermal temperature. Coupling of the fluid and solids phase is completed by directly imposing a no-slip boundary condition at the particle surface via the link-bounce-back method [64]. The force and torque applied to a given particle is found by accumulating the momentum transfer between fluid and solid at every boundary node. Inter-particle collisions are modeled via the hard sphere approach and thus are binary and instantaneous. Newton’s equation of motion is then employed to find the new velocity (translational and rotational) and position of each particle. 2.2. Random walk particle tracking (RWPT) 2.2.1. Single phase systems Similar to the LBM, RWPT does not directly involve the advection-diffusion equation. Instead, the positions of many Brownian tracers are monitored as they undergo movement that depends upon the local velocity field and random fluctuations. These tracers can cross the fluid-solid boundary, so the temperature within the fluid and solid particles is tracked in time. From the theory of stochastic differential equations, a connection can be drawn between the distribution of tracer particles and the advection-diffusion equation [67]. Here we utilize the explicit Euler–Maruyama method to evolve the position of each tracer [54,68]:

√

ri (t + t ) = ri (t ) + A(ri , t )t + ξ (t ) · B(ri , t ) t

(4)

where ri is the position of tracer ‘i’, A is the drift vector, ξ is a random vector whose entries are sampled from a Gaussian distribution with zero mean and unit variance, B is the diffusive displacement tensor, and t is the random walk time step.

4

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

Note that t in Eq. (4) may be set independently of t in Eq. (1) but should not exceed it. The diffusion tensor in Eq. (4) is evaluated at the position before the step (ri ), which implies Itô calculus and the evaluation of the stochastic integral at the left-points of the interval. In addition, the Euler–Maruyama method can be shown to exhibit strong convergence in time at a rate of 1/2 and weak convergence at a rate of 1 [69]. If A and B are chosen in the following manner

⎡√ A(ri , t ) = u(ri , t ),

B=⎣

2α 0 0

√0

2α 0

⎤

0 ⎦, 0 √ 2α

(5)

where α is the thermal diffusivity, the resulting tracer distribution obeys the Fokker–Planck equation for the advectiondiffusion equation with constant diffusivity [67]:

 T  ∂f 1 ≡ −∇ · (A f ) + ∇∇ : BB f = −∇ · (u f ) + α ( f ) ∂t 2

(6)

where u(ri , t ) is the velocity at the tracer position (ri ) before the step is completed (found via trilinear interpolation of the LBM velocity field in this work), f is the tracer density distribution (probability of finding a tracer between r and r + dr), T

B is the transpose of the displacement tensor. Note that RWPT may be coupled to any hydrodynamic solver (utilized to obtain the fluid velocity u(ri , t )) and that LBM was chosen for the present work due to its favorable application to gas– solids flows. Since Eq. (6) governs the evolution of the probability distribution function ( f ), the ‘white noise’ inherent in a stochastic method has been averaged out. From Itô’s chain rule, Eq. (4) will recover the advection-diffusion equation plus an additional term given by α ∇( f )dW t /dt, where W t is a Wiener process. The additional term corresponds to fluctuations about Eq. (6) whose mean value is zero. The stochastic fluctuations may be reduced at the expense of more Brownian tracer −1/2 particles and will scale as C t , where C t is the total tracer concentration [58,59]; see Section 3.2. For constant thermal properties and negligible shear work, the general energy balance simplifies to the advection-diffusion equation

∂T + ∇ · (uT ) = α ( T ) ∂t

(7)

where T is the thermal temperature. Since the form of Eq. (7) is consistent with Eq. (6), the stochastic differential (Eq. (4)) solves the advection-diffusion equation (Eq. (7)). The probability distribution function in Eq. (6) becomes analogous to the thermal temperature distribution in Eq. (7) – i.e., the tracer concentration represents the thermal temperature field. Alternatively, for the assumptions here, Eq. (7) may be formulated in terms of the specific enthalpy (by way of the thermodynamic relations). For this case, the tracer concentration would become proportional to thermal energy. The tracer concentration was chosen to represent temperature since it is consistent with the solutions utilized for verification here and easier to visualize conceptually. In RWPT, the position of tracer particles is evolved according to a stochastic differential equation (Eq. (4)). As a result, their spatial distribution will obey the Fokker–Planck equation for advection-diffusion. Similarly, dissipative particle dynamics (DPD) involves the update of ‘lumps’ of particles according to a stochastic differential equation, and thus will also obey a Fokker–Planck equation [70–72]. In contrast to RWPT, the stochastic differential in DPD is Newton’s second law. Therefore, the forces acting upon a lump of particles due to external sources and collision interactions (conservative, dissipative, and stochastic) are modeled in DPD, but not in RWPT. Since the stochastic process in DPD corresponds to random dissipations due to collisions, rather than random fluctuations in displacement (RWPT), the Langevin equation for DPD leads to a more complex form for the Fokker–Planck equation that accounts for shear work (not resolved with RWPT) [70,71]. While still a mesoscopic method, DPD is far more akin to molecular dynamics than RWPT and will thus incur a significant computational overhead in comparison with RWPT. 2.2.2. Development of tracer-based method for multiphase systems For the case of multiphase systems, the diffusivity may be a spatially varying quantity (B = B(r)), since different phases may have different diffusivities. Under such circumstances, if the drift vector (A) and displacement tensor (B) are naively chosen to match those in Eq. (5), the resulting tracer distribution will not correspond to the macroscopic transport equation – i.e., gradients in the diffusivity field will not be accounted for. If the diffusivity field is smooth, the drift vector may be defined in an alternative manner (A ≡ u∗ = u +∇ · α (r)) so that the form of the Fokker–Planck equation (Eq. (6)) is consistent with the advection-diffusion equation (proper macroscopic transport equation is recovered). Due to the dependency upon the gradient of the diffusivity field, the drift velocity (u∗ ) may become infinite in a multiphase system – e.g., the diffusivity is discontinuous across a particle–fluid interface. Here we consider the case of a gas–solid mixture in which the thermal diffusivities associated with each phase are not equal and do not have any spatial or temporal dependence. Therefore, a singularity in the drift velocity (u∗ ) will be present in the systems considered here and special considerations must be given to ensure the correct macroscopic phenomena is preserved. In the RWPT method, it is most commonly assumed that the diffusivities associated with each phase are equal. However, a few techniques have been developed for random walks on discontinuous fields. The techniques can be broadly grouped into the ‘reflective barrier’ and ‘interpolation methods’; see [56,61,62]. The random walk scheme developed here differs from both the reflective barrier and interpolation methods and is most analogous to that of

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

5

Fig. 1. An illustration of the geometry for steady-state diffusion in a composite medium. The thermal temperature (tracer concentration) and volumetric heat capacity of both materials are equal (T 1 = T 2 = T , C v ,1 = C v ,2 = C v ), while the diffusivities (conductivities) are not equal (α1 = α2 , k1 = k2 ), and yield a discontinuity at the interface (dashed line).

[60]. Note that in the present work, as well as [56,60–62], the interface between two phases yields a discontinuity in the diffusivity field, but the diffusivity within a phase is still assumed to be constant. The approach developed by [60] defines two time-scales that are associated with the two phases. The time-scales in [60] are related to one another by the ratio of the diffusivities (also assumed constant within a phase in [60]) and the ratio of the tracer concentrations in each phase at equilibrium (referred to as the partition coefficient). Here we employ an interfacial tracer balance to ensure that the net tracer flux is zero when the tracer concentrations in each phase are equal (partition coefficient of 1). An illustrative example as to why the random walk algorithm requires modification can be understood via the simple case of diffusion over a domain with a discontinuous thermal diffusivity – e.g., a solid particle suspended within a fluid. If the system is initialized with a uniform tracer concentration (thermal equilibrium), the tracer concentration should remain homogeneous since there are no driving forces (tracer concentration gradients). However, the system will evolve to a state where the concentration of tracers within the high diffusivity medium is less than that of the low diffusivity medium [56] – i.e., equilibrium will be defined by a state where the high diffusivity medium is cooler than the low diffusivity medium and not by a state of equal thermal temperatures. Such behavior is physically incorrect and serves as a basis for how the random walk algorithm should be modified. By considering the case of steady-state diffusion within a composite medium (Fig. 1), analytical methods may be employed to determine the appropriate modifications to the random walk framework, such that the proper equilibrium state is restored for the case of unequal diffusivities. The system under consideration is static and comprised of two materials which have equal temperatures (tracer concentration) and volumetric heat capacities but different thermal diffusivities (conductivities); see Fig. 1. For the case considered, the tracer update (Eq. (4)) only retains the stochastic term √ (ri (t + t ) = ri (t ) + ξ (t ) 2αi t). Therefore, the criteria for a tracer to reach the interface during a step takes the following form:

 ξx 2αi t ≥ x i

(8)

where x i is the separation distance between the current tracer position and the interface, and

ξx is the random number sampled in the x-direction. Since the stochastic terms (ξi ) are sampled from a Gaussian distribution, the probability of reaching the interface, for a given separation distance (x i ), can be found by calculating the probability of selecting a value such that Eq. (8) is satisfied:

∞

 

P xi ≡

∞ ξx (τ )dτ =

x

x

 i

1

√

2π

exp

−τ 2 2

 dτ =

1 2

erfc

x i

2

√

αi t

 (9)

 i

2αi t

2αi t

√ ξi (τ ) = 1/ 2π exp(−τ 2 /2) is the Gaussian distribution from which the stochastic terms are sampled. The probability of reaching the interface ( P (x i )) may then be utilized to complete a tracer balance across the interface and over a single random walk time step (t): where P (x i ) is the probability of reaching the interface, and



Nt =

 

C 1 (r) P x 1 dV1 −



 

C 2 (r) P x 2 dV2

(10)

where N t is the net number of tracers crossing the interface, C i denotes the concentration of tracers in medium ‘i’ and dVi is the differential volume occupied by medium ‘i’. Since the tracer concentration is uniform at steady-state (C = C (r)), the integrals in Eq. (10) may be directly evaluated:

6

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

L 1 C 1 Ac

2 0



1

erfc



x 1

2

√

α1 t



 

C 1 P x1 dV 1 −

dx = C 1 A c





∞ erfc(χ )dχ = C 1 A c

α1 t 0

 

Ac

C 2 P x2 dV 2 = √ [C 1

π



α1 t π

(11)



α1 t − C 2 α2 t ]

(12)

where A c√is the surface area of the interface, L 1 is the system length scale in the direction normal to the interface (x 1 ), and χ = x 1 /2 α1 t. Since the system length scales are large with respect to the modulus of the diffusive step, the bounds of integration in the middle equality of Eq. (11) may be simplified – i.e., as x i → L i , χ → ∞. By definition, the net number of tracers crossing the interface should be zero at equilibrium (N t = 0). By imposing the equilibrium condition upon Eq. (12), it can be seen √ that the tracer concentration (temperature) within each material is not the same unless the diffusivities are equal (C 1 /C 2 = α2 /α1 ). Furthermore, the analytical solution is in agreement with the qualitative behavior previously asserted (tracers will be biased to the lower diffusivity phase at equilibrium). To rectify the phase biasing issue and restore the proper definition of equilibrium, a twofold approach is taken to modifying the random walk algorithm. First, the displacement matrix (B) is modified for the higher diffusivity medium and a probability of undergoing a diffusive step ( P move,i ) is introduced. Specifically, the magnitude of the displacement matrix is increased for the higher diffusivity medium while the probability of undergoing a stochastic step within the higher diffusivity medium is reduced; see Eqs. (13) and (14).

B=

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

⎤

2α12

α2



0

⎥ ⎥ ⎥ ⎥ 0 ⎥ ⎥  ⎦ 2α12

0

0

2α12

α2

0 0 α2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎡√ ⎤ ⎪ ⎪ 2α2 √ 0 0 ⎪ ⎪ ⎪ ⎪ ⎣ 0 ⎪ 2α2 √ 0 ⎦ ⎪ ⎩ 0

(13)

for α1 > α2

(14)

2α2

0

α2 , P move,1 = α1

for α1 > α2

P move,2 = 1

Returning to the calculation of the net number of tracers crossing the interface (Eq. (12)), the requirement of zero net tracer flux at equilibrium (uniform concentration) is satisfied.





 

C 1 P move,1 P x1 dV 1 −

 





α2 C 2 P x2 dV 2 = √ C1 α1 π Ac



 α12 t − C 2 α2 t = 0 α2

(15)

While the appropriate definition for equilibrium is restored, it is not obvious if the modifications will alter the dynamics within the higher diffusivity phase, due to the alteration of the displacement matrix and introduction of the probability of moving. To show that the appropriate macroscopic behavior is retained, the mean squared displacement within the higher diffusivity material is considered:

 2 ri = P move,1 1

ξi = √

2π





ξi 2

∞

τ exp −∞

α12 t α2 −τ 2 2

2  + (1 − P move,1 )(0) = 2

 d τ = 0,

 2 1 ξi = √

2π

α2 α12 t  2  ξ = 2α1 t α1 α2 i

∞

τ 2 exp −∞

−τ 2 2

(16)

 dτ = 1

(17)

where · denotes an average over all Brownian tracers, ξi is the first moment of the stochastic term, and ξi2 is the second moment of the stochastic term. The mean squared displacement of tracers within the higher diffusivity material (Eq. (16)) matches the classic single-phase solution [56]. Inclusion of these modifications into the RWPT algorithm yields the following update scheme:



r1 (t + t ) = r1 (t ) + u(ri , t )t + ξ (t ) 2





α12 α2 t 1 − H ζ − α2 α1



r2 (t + t ) = r2 (t ) + u(ri , t )t + ξ (t ) 2α2 t where H is the Heaviside function and ζ is a random number sampled from a uniform distribution on the span [0, 1].

(18)

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

7

The second part of the approach addresses how tracers that reach the interface should be handled. For the system to be in equilibrium, the net tracer flux everywhere in the domain must be zero. Applying the tracer balance in Eq. (15) to a plane located within medium two at a separation distance d from the interface yields the following:

∞ N t (d) = C 1 A c

d

 

P move,1 P x 1 dx 1 + C 2 A c

 

P x 2 dx 2 − C 2 A c

0

d

∞

 

P x 2 dx 2

(19)

0

where N t (d) is the net number of tracers crossing the plane, the first integral is the number of tracers from medium one crossing the plane from left to right, the middle integral is the number of tracers from medium two crossing the plane from left to right, and the third integral is the number of tracers from medium two crossing the plane from right to left. Simplification of Eq. (19) shows that

α2 N t (d) = C 1 A c α2

∞



x 1

α2 dx 1 − C 2 A c erfc √ 2 2 α2 t α2 1

d

∞

1 2

erfc

d

x 2

2

√

α2 t



dx 2 .

(20)

For the case of C 1 = C 2 , Eq. (20) is not zero. Therefore, if no modification is made to the tracers crossing the interface then the zero-flux condition will not be satisfied at the plane – i.e., even though the net flux at the interface is zero the same will not hold true for a plane recessed within either medium. To rectify this issue, tracers that cross the interface are first propagated to the position that their step intersects with the interface and their diffusion matrix is converted to that of the new medium. The following step away from the interface must then be completed in a manner that ensures N t (d) = 0. Utilizing Eq. (15) and the fact that tracers are first moved to the interface, the tracer balance at the plane becomes:

√

∞ x α2 t 1 N t (d) = C 1 A c √ P (d) − C 2 A c erfc √ 2

π

2

dx 2

α2 t

2

d

√



(21)

√

where C 1 A c α2 t / π is the number of type one tracers that reach the interface (Eq. (15)), and P (d) is the probability of the tracer reaching the plane from the interface. By noting that the probability distribution is given by the integration of the displacement distribution (Eq. (9)), Eq. (21) can be reformulated as

C 1 Ac

√ α2 t √

∞

ξ (τ )dτ = C 2 A c

π

∞

1 2

erfc

2

d

 d

2α2 t



x 2

√

α2 t

dx 2

(22)

where ξ (τ ) is the distribution that the stochastic step away from the interface is sampled from. Inspection of Eq. (22) shows that if ξ (τ ) is given by

 ξ (τ ) =

π 2

τ

erfc √

2

 ,

(23)

then the probability of a tracer reaching the plane from the interface will be

 P (d) =

∞

π

τ

erfc √

2  d

2



 dτ =

√



 d π d d2 − exp − , erfc √ √ 4α2 t 2 α2 t 2 α2 t

(24)

2α2 t

and the tracer balance at the plane (Eq. (21)) will be zero for any d when C 1 = C 2 . To sample from Eq. (23) (ξ (τ )), the cumulative distribution function (C (x)) may be utilized. If values (u ) are sampled from a uniform distribution on 0 to 1 (U (0, 1)) then the inverse cumulative distribution function (C −1 ) may be evaluated at these points to recover the original distribution (ξ (τ )) – i.e., for values u sampled from U (0, 1)), the evaluation of x = C −1 (u ) will yield a distribution of x values that obey Eq. (23). The cumulative distribution function (C (x)) for Eq. (23) is given by

 C (x) =

π

x

τ

erfc √

2 0

2

 dτ =



π 2

x

xerfc √

2



− exp

−x2 2

 + 1.

(25)

Since Eq. (25) cannot be analytically inverted, a numerical root finding method must be employed. Use of Newton’s method for root finding yields the following update scheme

8

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

Fig. 2. An illustration of the RWPT workflow.

τn+1 = −

1 − u − exp(



−τn2 2

τn π erfc( √ ) 2

)

(26)

2

where τn+1 is the root at the n + 1 iteration, u is the value sampled from U (0, 1), and τn is the root at the n iteration. While Newton’s method displays quadratic convergence, iteration for each tracer crossing the interface may become computationally expensive as the number of particles or tracers becomes large. Alternatively, Eq. (23) may be approximated by the √ right half of a Gaussian distribution (2/ π exp(−τ 2 )). The difference between the Gaussian distribution and Eq. (23) is relatively small and allows the stochastic number generator to be reused (but with a variance of 1/2). Differences between the Gaussian approximation and Eq. (23) will result in errors in the tracer flux across the plane. Since the tracer flux is related √ to the integral of the displacement distribution, the integral of Eq. (23) minus 2/ π exp(−τ 2 ) may be utilized to gauge the magnitude of the errors associated with the Gaussian approximation. Evaluation of said integral shows that a maximum deviation of 0.05 occurs at τ ≈ 1.1 and in general the error is of O (10−2 ). Since Eq. (23) and the Gaussian distribution both decay rapidly and the tracer step is small with respect to the LBM cell length, the errors will manifest themselves through sub-grid errors in the temperature field near the interface. A numerical examination of the errors associated with the Gaussian approximation are given in Section 3.2. Except for the error analysis in Section 3.2, the Gaussian distribution was utilized for the step away from the interface in this work. Overall, tracers that cross the interface are first moved to the position their step intersects with the interface, then the diffusive matrix is set to that of the new medium, and the following step is sampled from the Gaussian distribution with a variance of 1/2 until the final position resides within the new medium. When utilizing the erfc distribution in Section 3.2, a symmetric reflection of Eq. (23) about the τ axis is utilized. Specifically, Eq. (26) is utilized to sample a value from Eq. (23) and then a value is sampled from U (0, 1), if the value is less than 0.5 the root is returned, else the negative root is returned. An outline of the algorithm is given in Fig. 2. Despite the simple nature of the system considered (equilibrium diffusion in a composite medium), the modifications made to the random walk algorithm are applicable to much broader classes of problems. If the tracer step is kept small with respect to the particle diameter (D p ), then the effects of particle curvature may be neglected and the analytical solution for the 1D system holds, due to isotropy of the Gaussian distribution. Furthermore, the velocity field considered here will be continuous and nearly divergence free, due to nearly incompressible flow with no-slip boundary conditions imposed upon particle surfaces, and thus, requires no modification. Therefore, the given modifications may be extended to 3D, gas–solid systems involving flow.

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

9

2.3. Boundary conditions for tracer-based method In addition to the new tracer algorithm, suitable boundary conditions must be defined for the random walk method before simulations may be completed. First, we formalize the definition of thermal temperature utilized in this work. As discussed at the end of Section 2.2.1, the thermal temperature may be taken to be proportional to the absolute tracer concentration. For this interpretation, a ‘hot’ boundary (T s ) would act as a source for tracers while a ‘cold’ boundary (T f ) would act as a sink for tracers. Under these conditions, the number of tracers within the system will change and the tracer arrays must be dynamically reallocated as tracers are destroyed and created at the cold and hot boundaries, respectively. Rather than represent the temperature absolutely by the tracer concentration, we map two tracer types to two different thermal temperatures. In this interpretation, a hot boundary acts as a source for ‘b’ tracers and a cold boundary acts as a source for ‘a’ tracers. By utilizing two tracer types we avoid the need for dynamic reallocation of the tracer array since tracers are not created or destroyed but merely converted in type. Tracers labeled as type b correspond to T s (T hot ) and a type tracers correspond to T f (T cold ). The local temperature and dimensionless temperature then become:

T (r, t ) = T s

θb (r, t ) ≡

C b (r, t ) Ct

+Tf

T (r, t ) − T f Ts − T f

=

C a (r, t )

(27)

Ct C b (r, t )

(28)

Ct

where C b is the concentration of b tracers, C a is the concentration of a tracers, C t is the total tracer concentration, and θb is the dimensionless temperature. 2.3.1. Periodic boundaries The basis for periodic boundaries is that the simulation domain behaves as if it is surrounded by its translational image. To achieve the desired behavior, each of the periodic planes is treated as if they are next to one another. If the regions on either side of a periodic plane are adjacent (similar to Fig. 1), then the flux of tracers across a periodic boundary (due to advection and diffusion) will move to the region adjacent to the other periodic plane. For example, consider a system that is treated as periodic in the x direction with planes at xmin and xmax . If a tracer is initially adjacent to the xmin plane and undergoes a displacement of r (right side of Eq. (18)) that leads to a final position which is less than xmin , then the periodic treatment requires that



r x (t + t ) = xmax + xmin − r x (t ) − r x



(29)

where r x (t + t ) is the final tracer position in the x direction after crossing the periodic plane, r x (t ) is the initial tracer position near the xmin plane, and r x is the displacement in the x direction due to advection and diffusion. The above treatment ensures that the flux of tracers out of the domain is the same as the flux of tracers across the other periodic plane and back into the other side of the domain. 2.3.2. Dirichlet boundaries The maintenance of a constant temperature at a boundary (T = T 0 ), or within a sub volume (particle), is a classic type of boundary condition. From Eqs. (27)–(28), it can be seen that the maintenance of a constant temperature is equivalent to the maintenance of the concentration of tracer types (C a (r, t ), C b (r, t )). In general, the particle being maintained at a constant temperature will not be in thermal equilibrium with the surrounding fluid. Therefore, the flux of tracers into the particle will have a different distribution of tracer types (temperature) than that required to achieve T = T 0 . However, if the types of the tracers entering the particle are modified to be consistent with T 0 , then the particle will remain at the desired temperature. To do this, the dimensionless temperature associated with the boundary is utilized (θb,0 ). The dimensionless temperature characterizes the concentration of type b tracers and spans from 0 to 1. The type of each tracer moving into a particle is set by sampling a value from U (0, 1) and if the value is less than or equal to θb,0 , the tracer type is set to b, else it is set to a. Setting the tracer type in this manner ensures that the temperature of the tracers entering the particle is consistent with the boundary condition and the particle remains at θb,0 . In addition to maintaining a particle at a constant temperature, it is often desirable to maintain a boundary plane at a constant temperature. For the case of a constant temperature boundary plane, the tracers are not resolved outside of the simulation domain and thus the flux of tracers across the boundary plane back into the simulation domain is not directly resolved. Rather than attempt to resolve the boundary region with tracers, we seek to derive a treatment for tracers that cross the boundary plane. In order to only operate on tracers that reach the boundary plane, the net flux of tracers across the boundary plane must be zero. This will only be the case if u · n → 0 at the boundary and the flux across the plane is due solely to diffusion (j = k∇ T ). To mathematically illustrate this fact, we introduce a uniform velocity field (u x ) to the probability in Eq. (9) to yield:

10

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

 

∞

1

P x ≡

√ x ±u t √ x 2α t

2π

exp

−τ 2 2

 dτ =

1 2

erfc

x ± u x t 2

√

α t

 (30)

where the ± corresponds to the two cases where the velocity field points towards the boundary (−) (u · n < 0) and away from the boundary (+) (u · n > 0). Clearly, the probabilities of reaching the boundary plane from either side are not equal and the net number of tracers crossing the boundary plane will not be zero. This behavior is entirely physical and due to the advective contribution to the tracer flux. However, since the net flux of tracers through the boundary plane is not zero for u · n = 0, operations cannot be solely restricted to tracers reaching the boundary plane. The most common case of a constant temperature boundary where u · n = 0 is an inflow plane. The discussion for the inflow plane will be given below in Section 2.3.4 since it is found that the extra tracer operations associated with the inflow plane can be naturally handled with the outflow plane. Therefore, we proceed for the case that u · n = 0 at the boundary plane and the tracer flux is solely due to diffusion. In order to replicate the flux originating from a boundary plane held at θb,0 , we consider a mirror image of the tracer distribution near the boundary plane and set the tracer types in this mirror image according to θb,0 . The tracers in the mirror image undergo displacements that mirror the real tracers near the boundary plane. By doing this, it is ensured that the flux of tracers from the fictitious mirror image into the simulation domain is equal to the flux of tracers from the simulation domain into the mirror image region. The post-step positions of the mirror image tracers become the specular reflection of their real tracer counterpart. Therefore, tracers that cross the boundary should be specularly reflected back into the simulation domain and their type set according to θb,0 . For a boundary plane located at xmin , tracers that undergo a displacement of r that leads to a position which is less than xmin should be handled by first modifying their position to

r x (t + t ) = 2.0xmin − r x (t ) − r x

(31)

and then setting the tracer type according to θb,0 – i.e., the process of sampling from U (0, 1) given above. This process ensures that the net flux of tracers back into the domain is consistent with the diffusive flux that would be obtained from a region held at θb,0 . 2.3.3. Neumann boundaries Boundary conditions of the second type, or specified flux boundaries are also of great practical significance. Here we will handle the case of an adiabatic boundary in which the specified flux is zero. Similar to the Dirichlet boundary in Section 2.3.2, the velocity field at the boundary plane is of great importance. We will again consider the case where u · n = 0 at the boundary plane and the tracer flux is solely due to diffusion. By again considering the mirror image treatment at the end of the Section 2.3.2, it becomes clear that if the tracers within the mirror image have the same distribution of tracer types (temperature) as their real counterparts then the net flux of b tracers will be zero (∇θb = 0). This amounts to the exact same treatment as that proposed for a Dirichlet wall, but without modifying the tracer types after specular reflection. 2.3.4. Inflow/outflow boundaries As prefaced in Sections 2.3.2 and 2.3.3, the effect of the velocity field at the boundary plane requires additional work. It is common to impose a fully developed outflow boundary condition ( ∂∂ un = 0, ∂∂ nT = 0) in conjunction with a uniform inflow boundary condition (u = uo , T = T 0 ). For these cases, the criteria of u · n = 0 will not be satisfied. In Section 2.3.3 a method was proposed for achieving ∂∂ nT = 0 where tracers crossing the boundary are specularly reflected back into the domain. For the case that the velocity field points out of the domain, it becomes clear that such a treatment (impenetrable) would lead to an accumulation of tracers at the boundary since tracers will move in a deterministic manner towards the boundary but never vacate the domain. For the opposite case of a velocity field that points into the domain, there is an advective flux of tracers into the simulation region that must be accounted for – i.e., the imbalance in tracer flux arising from the probabilities in Eq. (30). Thus, some tracers must be allowed to vacate the outflow plane and some tracers must be seeded back at the inflow plane if the desire boundary conditions are to be maintained. Since the velocity field in this work will be approximately divergence free (∇ · u ≈ 0) the Gauss divergence theorem may be applied to the domain





∇ · u dv =

u · n d ≈ 0

(32)

to show that the advective flux out of the domain is approximately equal to the advective flux into the domain. This observation led to a treatment where tracers that cross the outflow plane are randomly re-seeded at the inflow plane and their type is set according to the temperature at the inflow plane. Tracers that cross the inflow plane are specularly reflected back into the domain and their type is set according to the boundary temperature (θb,0 ). This two-fold approach accounts for the advective flux into and out of the domain as well as the diffusive flux at the inflow plane. While this treatment yields acceptable results for the simulations considered here, it is not the most general treatment for a set of inflow-outflow boundary conditions and is only valid for advection dominated flows (large Peclet number). Specifically, the diffusive flux at the outflow plane is not taken into account by this treatment. To see this, consider the case that the

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

11

velocity field asymptotically approaches u · n = 0 and the boundary becomes diffusion dominated. The Neumann boundary treatment derived in Section 2.3.3 is valid for this case and would necessitate that all tracers be specularly reflected back into the domain. Therefore, the open boundary utilized here for the outflow plane is valid for advection dominated flows (large Peclet number) while the reflective barrier is valid for diffusion dominated flows (small Peclet number). A general implementation for outflow boundary conditions that takes into account the advective and diffusive flux, and thus is valid for all Peclet numbers, has been derived by the authors and is presently under review. We omit that derivation here and refer the interested reader to our future work. 3. Verification simulations To ensure that the new hybrid LBM–RWPT framework (implementation of Eq. (18) and Fig. 2) is correct, a variety of systems are simulated. Outputs from the new LBM–RWPT code are compared to analytical solutions, finite difference numerical solutions, empirical correlations, and previous DNS works. Note that simulation parameters for the LBM (u, t , D p , ρ , η ) and RWPT method (α ) must be non-dimensionalized since the time step and grid spacing in LBM–RWPT are scaled to be unity [73]. To prevent ambiguity, quantities in physical space will be given alongside units. 3.1. Case 1 (analytical verification): effective diffusivity of a static gas–solids mixture The first case considered with the new code is conduction within a static gas–solids matrix. The matrix is unbounded and comprised of randomly located, non-overlapping, spheres within a fluid. A thermal impulse is applied to the matrix and the resulting conduction allows an effective diffusivity for the two-phase matrix to be quantified. The system in question is a natural starting point for verifying the new methodology and critically tests three key aspects. First, the simulations assess whether the phase biasing by tracers has been rectified. Namely, if tracers preferentially sample the lower diffusivity medium, then the effective diffusivity of the matrix will be significantly reduced from the true value. Second, the simulations are 3D and employ spherical particles within a cubic domain. Therefore, the assertion that the analytical solution may be extended to more complex geometries is evaluated. Third, the simulations are not at steady-state, and thus show whether the derived solution (Eq. (18)) still holds in the transient state (based on definition of the equilibrium state). Furthermore, an analytical solution [74] as well as rigorous bounds [75] are known for the effective diffusivity of the matrix. The analytical solution of [74] is given by:







αs − α f εs + ψ εs2 α s + 2α f αs ψ = 0.11, 0.208, 1.23, 3.90 for = 0.5, 2.0, 5.0, 50.0 αf αeff = α f 1 − 3

(33)

where αeff is the effective diffusivity of the matrix, α f is the diffusivity of the fluid, αs is the diffusivity of the solids, εs is the solids volume fraction, and ψ is the coefficient of the second order solids volume fraction term; while the bounds determined numerically [75] are:

  1

αeff ,L B =

α

−



− 1

4ε g εs ( α1 − α1 )2 s f 6 1 1 α f + (4ε g + 2ω)( αs − α f )

ε g εs (αs − α f )2 3α f + (ε g + 2ω)(αs − α f ) ω = 0.21068εs − 0.04693εs2 + 0.00247εs3

(34)



αeff ,U B = α −

where αeff , L B is the lower bound for the effective diffusivity, ε g is the fluid volume fraction, αeff ,U B is the upper bound for the effective diffusivity, ω is a structure factor (determined by [75] via Monte Carlo simulation), and X = ε g X 1 + εs X 2 is the volume fraction weighted average. To determine the αeff within the new code, we employ the methods as outlined in [60,61]; where αeff is given by the time rate of change of the mean dispersion of the tracers:



M2 = rr − r 2 lim

t →∞

1 dM2 2 dt



(35)

= α eff

where M2 is the dispersion tensor and · refers to an average over all Brownian tracers. For the simulations here, the tracers are initialized uniformly in a plane located at the vertical centerline of the domain (( y min + y max )/2). Therefore, dispersion only occurs in the vertical y-direction and Eq. (35) may be simplified to

 1 d  r y r y − r y 2 . t →∞ 2 dt

αeff = lim

(36)

12

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

Table 1 Effective diffusivity simulations. Geometry and operating conditions LBM-RW space Number of cells R p,L B M

αf αf αf αf

30 × 30 × 30 3 0.21 2.0 1.12 × 10−2 5.60 × 10−3 5.60 × 10−3 5.60 × 10−4

Physical space Rp

250 μm

αf

2.23 × 10−5 [ ms ]

εs Ct

αs αs αs αs

5.60 × 10−3 1.12 × 10−2 2.80 × 10−2 2.80 × 10−2

2

Fig. 3. (a) The effective diffusivity obtained for each matrix configuration and (b) the resulting confidence interval, for a diffusivity ratio ( ααs ) of 0.5. f

Since the effective diffusivity will depend upon the matrix structure (particle locations), ten random particle configurations are generated for every diffusivity ratio (αs /α f ). The number of particles within the matrix is held constant at 50 (εs = 0.21). A 95% confidence interval is utilized to compare the output obtained from the new code to the closures proposed by [74,75]. All domain walls are simulated with periodic boundary conditions as given in Section 2.3.1 and an overview of the simulation conditions is given in Table 1. As shown in Figs. 3–6 for the various diffusivity ratios examined, outputs from the new code are observed to be in agreement with prior analytical predictions [74] and numerical bounds [75], whereas the original code [49,50], which does not contain any modification for discontinuous diffusivity fields, underpredicts αeff because tracers preferentially sample the lower diffusivity medium (phase biasing), resulting in an αeff that is reduced from the true value. By contrast, the modifications developed here (Eq. (18) and Fig. 2) ensure that the tracers properly sample both phases. Note that the lower bound in [75] is essentially identical to the analytical solution obtained by [74]. This can be traced back to the homogeneous particle distribution assumed by [74], whereas [75] obtains bounds upon αeff which depend upon the arrangement of particles. Ultimately, no statistical difference is found between the αeff obtained from the new code and that of [75]. 3.2. Case 2 (analytical verification): hot sphere in unbounded diffusion For the case of a hot sphere (T s ) suddenly immersed in an unbounded, static, cold fluid (T f ), the spherical heat balance only retains the temporal and radial components. If the sphere is kept at a constant temperature (Dirichlet), the analytical solution for the temperature profile and Nusselt number are given by [63]:

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

13

Fig. 4. (a) The effective diffusivity obtained for each matrix configuration and (b) the resulting confidence interval, for a diffusivity ratio ( ααs ) of 2. f

Fig. 5. (a) The effective diffusivity obtained for each matrix configuration and (b) the resulting confidence interval, for a diffusivity ratio ( ααs ) of 5. f

T (r , t ) − T f Ts − T f Nu (t ) ≡

=

hD p kf

Rp r



1 − er f 2

=2+ √

π

r − Rp 

4α f t

Rp



αf t

 (37) (38)

where R p is the radius of the sphere, h is the heat transfer coefficient, and k f is the thermal conductivity of the fluid. The solution method utilized to obtain Eqs. (37)–(38) (similarity transform) inherently assumes that the fluid is unbounded. For the simulations considered here, the fluid domain will in fact be finite. However, if short time-scales are considered, the

14

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

Fig. 6. (a) The effective diffusivity obtained for each matrix configuration and (b) the resulting confidence interval, for a diffusivity ratio ( ααs ) of 50. f

simulation will behave as an unbounded system since the domain boundary conditions (periodic) will not play a significant role. By contrast, at long time-scales (steady-state), the solutions given in Eqs. (37)–(38) are not valid and are entirely dictated by the particle and domain boundary conditions. To make a comparison between the new code and the analytical Nusselt number (Eq. (38)), a heat transfer coefficient (h) must be defined for the RWPT method. By definition, the heat transfer coefficient relates the heat flux (q ) to the far field thermal temperatures (q = h( T s − T f )). Here, for simplicity, we choose the temperatures such that T f = 0 and T s = 1 and T (r, t ) = θb (r, t ) = C b (r, t )/C t ). In the RWPT method, q will be proportional to the net number of b type tracers crossing the particle interface per time. Therefore, if the net flux of b tracers across the surface of the sphere can be converted to thermal energy, then the heat transfer coefficient may be directly evaluated (q = h( T s − T f ) = h). To accomplish this, the thermal energy corresponding to one b tracer within an LBM–RWPT cell (unit volume) is found via:

E t = V P S (C v , f  T )

= dx3P S

C v, f

1



(39)

Ct

where E t is the energy corresponding to C b = 1, V P S is the volume of a LBM–RWPT cell in physical space, C v , f is the volumetric heat capacity of the fluid (assumed equal to the solids, at this point), and dx3P S is the conversion factor between LBM–RWPT space and physical space. The heat flux and heat transfer coefficient then become:

h≡

q

T

=

Net # b tracers crossing A p ,LBM ∗ t LBM



dx3P S C v , f

1 Ct



1 dx2P S dt P S



(40)

where A p , L B M is the area of the sphere in LBM–RWPT space, t L B M is the number of LBM time steps over which the heat flux (net # of b tracers) is accumulated, and dt P S is the conversion factor between LBM–RWPT time and physical time. Note that Eq. (40) is simply converting the flux of b tracers (left term in parentheses) to thermal energy by way of the energy per tracer (middle term in parentheses) and then converting the sphere surface area and time back to physical space (right term in parentheses). To satisfy the required initial condition and boundary condition for the unbounded diffusion case, the methods given in Section 2.3.1 and Section 2.3.2 are employed. An outline of the simulation conditions is given in Table 2. The Nusselt numbers obtained from simulations are observed to be in strong agreement with the analytical solution (Eq. (38)) over a variety of diffusivity ratios; see Fig. 7. The local, dimensionless concentration of all tracers (θt (r) = (C b (r) + C a (r))/C t ) and b tracers (θb (r)) is illustrated in Fig. 8a and 8b, respectively. A cubic grid (4 × 2 × 4) is employed to decompose the domain and calculate the local concentrations. θt is observed to be homogeneous over the entire domain and equal to one (Fig. 8a). Since the system is initialized at equilibrium, with respect to all the tracers, this result supports the fact that the modifications to the RWPT algorithm (Eq. (18) and Fig. 2) properly maintain equilibrium (Eq. (15) is satisfied). By contrast, θb (r) not only varies spatially but decreases as the radial position from the particle surface is increased (Fig. 8b), which agrees with the analytical solution (Eq. (37)).

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

15

Table 2 Simulations of hot sphere in unbounded diffusion simulations. Geometry and operating conditions LBM space Number of cells R p,L B M Ct

αf αf αf

80 × 10 1.0 1.784 1.784 1.784

Physical space Rp

250 μm

αf kf C v, f

80 × 80

× 10−1 × 10−1 × 10−1

αs αs αs

8.920 × 10−2 1.784 × 10−1 8.920 × 10−1

2

2.230 × 10−5 [ ms ] W 2.624 × 10−2 [ mK ] J 1.176 × 103 [ m3 ]

Fig. 7. The Nu number achieved via simulation compared to the analytical solution, for a diffusivity ratio of ( ααs ) = 0.5, 1, 5. f

Fig. 8. (a) The dimensionless total tracer concentration (θt ) and (b) the dimensionless b tracer concentration (θb ) (temperature field) at a time of 1 × 10−3 [s] and ααs = 5, where the dashed line denotes the outline of the sphere. f

16

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

Fig. 9. The residual between LBM–RWPT and Eq. (38) is observed to exhibit 1/2 order scaling with tracer concentration.

Fig. 10. (Left) The thermal profile extracted with the use of Eq. (23) and the Gaussian approximation. (Right) The relative error in the heat flux at the first three spherical shells obtained from using Eq. (23) and the Gaussian approximation.

The accuracy of RWPT will depend upon the number of tracers utilized within the simulation domain (C t ). As C t increases, the stochastic fluctuations (see fluctuations in Fig. 7) will be attenuated and the Nusselt number predicted by RWPT will better agree with the analytical Nusselt number (Eq. (38)). The residual (E 2 ) between RWPT and Eq. (38) is evaluated at varying C t and shows that the method exhibits 1/2 order convergence with tracer count; see Fig. 9. The scaling −1/2 of the residual with tracer concentration (E 2 ∝ C t ) is consistent with that given in [58,59] (discussed in Section 2.2.1). Additionally, the accuracy of the Gaussian approximation for the step away from the interface is examined here. As discussed at the end of Section 2.2.2, the Gaussian approximation will have small deviations from Eq. (23) that cause errors in the flux at a plane located a distance d from the interface. From the advection-diffusion equation (Eq. (7)) it can be seen that the change in temperature is given by the divergence of the flux. Therefore, errors in the flux from the Gaussian approximation will yield small errors in the thermal temperature field near the interface. Test simulations with Eq. (23) and the Gaussian approximation at C t = 5.0 are conducted to assess the magnitude of these fluctuations. The thermal field and heat flux are extracted on concentric spherical shells that have a width of 0.2dx and are averaged over 3 simulations; see Fig. 10. The results show that the errors due to the Gaussian approximation may be detected in the thermal field (first point at rˆ = 1.01) but are only over length-scales that are small with respect to the LBM grid. To assess the impact of the Gaussian

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

17

approximation on the heat flux, the relative error in the heat flux obtained via the Gaussian approximation and Eq. (23) is considered. It should be noted that a direct comparison between the heat fluxes will inherently include the stochastic noise. While the tracer count has been increased to attenuate the stochastic fluctuations, the errors will be comprised of effects due to the Gaussian approximation as well as the random fluctuations. At very short time-scales, the Gaussian distribution is observed to mildly under-predict the heat flux in the spherical shells (∼1–2% error). However, at slightly longer timescales, the relative error tends to exhibit larger fluctuations that are centered around zero relative error. Examination of Fig. 7 shows that the stochastic fluctuations become more significant at longer time scales when the temperature near the surface of the particle is nearly homogeneous. Therefore, at longer time-scales, it appears that the spurious fluxes due to the Gaussian approximation become blended into the stochastic fluctuations and their exact quantification is challenging. Time averaging of the relative error over the entire simulation yields values of ∼0.45% for all three shells. By again considering the E 2 residual with Eq. (38), it is found that errors in the heat flux at the sphere’s interface, due to the use of the Gaussian approximation, is of O (10−3 ) which is nearly an order of magnitude lower than the E 2 in Fig. 9. This is consistent with the observation that the spurious fluxes are rapidly overshadowed by the stochastic fluctuations. For this reason, we consider that errors from the Gaussian approximation practically negligible in our presented results. Such errors should only be of concern when heat transfer over very short distances and time-scales needs to be accurately addressed. 3.3. Case 3 (numerical verification): cooling sphere in unbounded diffusion Under the conditions given in Section 3.2 (Hot Sphere in Unbounded Diffusion), the coupling of heat transfer between the sphere and the fluid is only one-way coupled. Namely, as the sphere transfers heat to the fluid it does not in turn cool but remains at a constant temperature. Here, we allow the sphere to cool, and thus, develop intra-particle temperature gradients in the radial direction. Due to the added complexity of coupling the two phases, an analytical solution cannot be found for the temperature profile. Instead, the simulation outputs are compared to results from a finite difference (FD) code. The FD algorithm discretizes the spherical heat balance on a staggered grid [76] over the span of [0, 3R p ] with a symmetry boundary condition at r = 0 [77] and an adiabatic boundary condition at r = 3R p [78]:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨

∂T = ⎪ ∂t ⎪ ⎪ ⎪ ⎩

6k1/2 t (T 1 − T 0 ) C v ,0 r 2 2k N −1/2 t ( T N −1 − C v , N r 2 1 [ C v ,i r i2

for i = 0 for i = N

TN)

(r i2+1/2 ki +1/2 ) T i +1 −(r i2−1/2 ki −1/2 +r i2+1/2 ki +1/2 ) T i +(r i2−1/2 ki −1/2 ) T i −1 r 2

(41)

] else

where ‘i’ denotes the node index and N is the total number of nodes (103 ). The temporal derivative in Eq. (41) is approximated by an implicit Euler scheme and the resulting matrix equation is solved by the generalized minimal residual method (GMRES). Since both phases are coupled in the current case, the amount of thermal energy gained by the fluid must be equal to the thermal energy lost by the particle. If an energy balance is performed at the fluid–particle interface, it can be seen that this criterion is only met if the volumetric heat capacities of both phases are equal (C v ,s = C v , f ):

 

1 1 1 + ( N b ,s→ f − N b , f →s ) C v , f =  N b,s C v , f (1 − γ )  E = ( N b , f →s − N b ,s→ f ) C v ,s Ct

 N b,s = ( N b,s→ f − N b, f →s ),

Ct

γ=

Ct

(42)

C v ,s C v, f

where  E is the change in energy between the two phases, N b, f →s is the number of b tracers crossing from fluid to solid, N b,s→ f is the number of b tracers crossing from solid to fluid, (C v ,s C1 ) is the energy per tracer in the solid phase, (C v , f C1 ) t t is the energy per tracer in the fluid phase,  N b,s is the net change in b tracers in the solid phase, and γ is the ratio of volumetric heat capacities. A net energy change when C v ,s = C v , f can be attributed to the fact that a tracer represents a different amount of energy in the solid phase than the fluid phase. Thus, if tracer types are merely exchanged across the interface then the energy balance will not be satisfied, unless C v ,s = C v , f . To address the challenge of unequal volumetric heat capacities, a type tracers entering the high volumetric heat capacity medium are converted to b type. Specifically, the number of tracers which should be converted is found by





1 1 +  N b ,s C v , f = 0 for γ > 1  E = ( N b,convert −  N b,s ) γ C v , f Ct

∴ N b,convert =

Ct

(43)

γ −1  N b ,s γ

  1 1 1  E = (− N b,s ) C v ,s C v ,s + ( N b,s − N b,convert ) = 0 for γ < 1 Ct

γ

Ct

(44)

18

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

Fig. 11. The dimensionless b tracer concentration (θb ) as a function of dimensionless radius (rˆ ) for (a)

γ = 0.1,

Fig. 12. The dimensionless b tracer concentration (θb ) as a function of dimensionless radius (rˆ ) for (a)

γ = 1,

αs α f = 0.5 and (b)

γ = 0.1,

αs α f = 0.5 and (b)

γ = 1,

αs α f = 5.

αs α f = 5.

∴ N b,convert = (1 − γ ) N b,s where N b,convert is the number of a tracers entering the high heat capacity medium that should be converted to b tracers. Since N b,convert is a function of the net number of b tracers crossing the particle interface (∼ q ), it is calculated at the end of the tracer loop and the conversion of the tracers from a to b is completed at the beginning of the following loop (before the tracer steps away from the interface). Simulation inputs for each heat capacity ratio (γ = 0.1, 1, 10) match the αs /α f = 0.5, 5 cases in Table 2 and the domain boundary conditions are identical to those in Section 3.2. The dimensionless temperature profiles (θb (ˆr ), rˆ = r / R p ) obtained from simulations are found to be in agreement with those obtained via the FD approximation; see Figs. 11–13 for different

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

Fig. 13. The dimensionless b tracer concentration (θb ) as a function of dimensionless radius (rˆ ) for (a)

γ = 10,

19

αs α f = 0.5 and (b)

γ = 10,

αs α f = 5.

values of heat capacity and thermal diffusivity ratios. The profiles are extracted from the code by decomposing the domain into 50 concentric spherical shells and calculating θb on each shell. While the shells naturally address the radial dependence of the problem, the shell volumes become small as rˆ → 0 and lead to noise in θb near the center of the particle. The noise observed as rˆ → 0 can again be reduced by adding more tracers, thereby attenuating the stochastic fluctuations and smoothly resolving the temperature field (tracer distribution) at smaller length scales. For some of the cases under consideration here, the particle conductivity (k s ) is small with respect to the fluid conductivity (k f ) (k s /k f < 1). Physically, this condition results in thermal energy being redistributed within the sphere at a slower rate than the fluid. For cases in which k s /k f = γ (αs /α f ) < 1, both the RWPT code and the FD code show substantial temperature gradients near the surface of the sphere; see Figs. 11 and 12a. A common assumption made within DNS or the discrete element method (DEM) is that spheres are isothermal – i.e., the Biot number (Bi = hD p /k s ) is small, and thus, intra-particle thermal gradients are negligible. Recent work has highlighted the challenges associated with resolving intra-particle temperature gradients [79,80]. The RWPT method developed here shows promise for being able to capture the heat transfer occurring in systems with large Biot numbers and does not suffer the computational overhead associated with carrying a highly resolved numerical grid for each particle. 3.4. Case 4 (validation): uniform flow past a hot sphere The final case considered with the new LBM–RWPT code is that of unbounded flow of a cold fluid past a stationary, hot sphere. Here, the ability of the new framework to characterize non-isothermal flows is assessed. Comparisons are made between the new code and experiments [81,82], as well as earlier numerical studies [28,83]. Thus, the ability of the new hybrid code to accurately capture flow and transport (validation) is completed, not just the recovery of a known accurate solution (verification) [84]. Previous works [21,47,63] show strong agreement between numerical predictions and stateof-the-art correlations for the Nusselt number. We compare the Nusselt number obtained via the LBM–RWPT code, at intermediate Reynolds numbers (Re = 10, 20, 30, 40, 50, 100) and a Prandtl number ( P r) of 0.7, to the correlations given by Ranz and Marshall [81], Whitaker [82], Feng and Michaelides [28], and Richter and Nikrityuk [83], respectively:

Nu = 2.0 + 0.6Re1/2 Pr1/3



1/2

Nu = 2.0 + 0.4Re Nu = 0.992 + P e

1/3

for 10 < Re < 104 , Pr > 0.7 0 .4

+ 0.06Re

Pr

1/3

1/3

+ 0.1Re

1/2

Nu = 1.76 + 0.55Re

 2/3

Pr

1/3

Pe

2/3

+ 0.014Re

4

Pr

1/3

for 3.5 < Re < 7.6 × 10 , 0.7 < Pr > 380

(46)

for 0.1 < Re < 4000, 0.2 < P e > 2000

(47)

for 10 < Re < 250, Pr > 0.7

(48)

where P e = RePr is the Peclet number. The hydrodynamic and thermal boundary conditions are as follows:



1. Uniform axial inflow (u = 0

U 0, y



(45)

0 ) and temperature (T = T f ) at the inlet plane ( y = y min );

20

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

Fig. 14. The Nu number as a function of sphere Reynolds number achieved via LBM–RWPT simulations and state-of-the art correlations.

Fig. 15. (a) The dimensionless total tracer concentration (θt ) and (b) the dimensionless b tracer concentration (θb ) (temperature field) at steady-state and R e = 40, where the dashed line denotes the outline of the sphere.

2. Zero gradient in the stream wise direction ( ∂∂ uy = 0, ∂∂ Ty = 0) at the outlet plane ( y = y max ); 3. Periodic in the x and z-direction; and 4. Constant-temperature sphere (T = T s ) with no slip (u = 0) at the surface. Enforcement of the axial inflow and the no slip boundary conditions are accomplished by the modified bounce-back method of Ladd [64] while the outflow is achieved by extrapolation of the incoming distribution functions [47,85,86]. The particle is maintained at a constant temperature in the same manner as discussed in Section 2.3.2. Thermal boundary conditions at the inflow/outflow boundary are achieved by the method discussed in Section 2.3.4. The grid resolution employed here (D p , L B M = 20) has been shown by previous work [21,47,63] to achieve the required accuracy over the aforementioned Reynolds numbers. Furthermore, to minimize boundary effects and achieve an unbounded flow, the domain sizes of [63,87] are utilized (8D p × 16D p × 8D p ). The thermal diffusivity ratio is held constant for these simulations (αs /α f = 5) and set to α f = 1.784 × 10−2 , αs = 8.92 × 10−2 in the RWPT algorithm (α f = 2.23 × 10−5 [m2 /s]). Nusselt numbers found by the new LBM–RWPT code are observed to agree with the correlations of [28,81–83]; see Fig. 14. Verification studies by [21,47,63] show similar results, with the Nusselt numbers found here being most complementary to those of [63]. Furthermore, the θt is again observed to be nearly homogeneous (Fig. 15a) while the θb displays a thin boundary layer of hot fluid around the sphere with a wake of warm fluid following behind it (Fig. 15b) – i.e., the expected qualitative behavior.

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

21

4. Concluding remarks A technique for simulating heat (or mass) transfer in multiphase systems via a hybrid lattice Boltzmann method (LBM) – random walk particle tracking (RWPT) method is presented here. Analytical methods were utilized to modify the RWPT algorithm for use in systems with discontinuous diffusivity fields, which past RWPT works were unable to account for. An interfacial tracer balance was combined with the definition of equilibrium to determine how the stochastic movement of tracers should be altered. By expanding the displacement tensor for the higher diffusivity medium and only allowing a subset of the tracers within the higher diffusivity medium to undergo stochastic movement, the tracer balance across an interface can be satisfied at equilibrium. Furthermore, a means for quantifying the heat transfer coefficient in the RWPT framework, as well as a technique for handling differing volumetric heat capacities, was developed. The energy per tracer is employed to convert the flux of tracers to a heat transfer coefficient as well as develop an energy balance at the interface of the two mediums. The algorithm is verified against four test cases to ensure that it recovers the correct macroscopic transport phenomena. Outputs from the new LBM–RWPT code are observed to be in agreement with past analytical/numerical solutions as well as empirical correlations and previous works. The hybrid method developed in the present work offers a variety of advantages for simulating systems with scalar transport. Namely, the RWPT is a meshless method that does not suffer from numerical dispersion and is straightforward to implement/parallelize, unlike classic finite difference and finite element methods. Furthermore, the present work demonstrates the potential of RWPT to capture intra-particle temperature gradients, and thus, may be a valuable alternative for the simulation of heat transfer in high Biot number systems, since it does not suffer the computational overhead associated with a highly-resolved grid for each particle. While not considered in the present work, the particles may also translate in space. However, further development of the method is required to handle such cases. Specifically, the tracers may change phase due to the discrete particle motion over the tracer field. The creation and destruction of solid or fluid type tracers must be handled in a rigorous fashion to ensure the proper physics are retained. Means for addressing medium motion will be the subject of future work. Conflict of interest statement There are no conflicts of interests to report. Acknowledgement The authors gratefully acknowledge funding support provided by the National Science Foundation under Grant CBET 1512630. References [1] J.A.M. Kuipers, W. Prins, W.P.M. Van Swaaij, Numerical calculation of wall-to-bed heat-transfer coefficients in gas-fluidized beds, AIChE J. 38 (1992) 1079–1091, https://doi.org/10.1002/aic.690380711. [2] A.V. Patil, E.A.J.F. Peters, V.S. Sutkar, N.G. Deen, J.A.M. Kuipers, A study of heat transfer in fluidized beds using an integrated DIA/PIV/IR technique, Chem. Eng. J. 259 (2015) 90–106, https://doi.org/10.1016/j.cej.2014.07.107. [3] D.J. Patil, J. Smit, M. van Sint Annaland, J.A.M. Kuipers, Wall-to-bed heat transfer in gas–solid bubbling fluidized beds, AIChE J. 52 (2006) 58–74, https:// doi.org/10.1002/aic.10590. [4] M.L. de Souza-Santos, Solid Fuels Combustion and Gasification: Modeling, Simulation, and Equipment Operations, second edition, CRC Press, Boca Raton, Florida, 2010. [5] M. Syamlal, D. Gidaspow, Hydrodynamics of fluidization – prediction of wall to bed heat-transfer coefficients, AIChE J. 31 (1985) 127–135. [6] Z. Zhou, Q. Hou, A. Yu, Particle scale simulation of heat transfer in fluid bed reactors, in: Heat Transfer – Mathematical Modelling, Numerical Methods and Information Technology, INTECH Open Access Publisher, Vienna, 2011, pp. 383–393. [7] Z.Y. Zhou, A.B. Yu, P. Zulli, Particle scale study of heat transfer in packed and bubbling fluidized beds, AIChE J. 55 (2009) 868–884, https://doi.org/10. 1002/aic.11823. [8] T.N. Cong, Y. He, H. Chen, Y. Ding, D. Wen, Heat transfer of gas–solid two-phase mixtures flowing through a packed bed under constant wall heat flux conditions, Chem. Eng. J. 130 (2007) 1–10, https://doi.org/10.1016/j.cej.2006.11.006. [9] A. Guardo, M. Coussirat, F. Recasens, M.A. Larrayoz, X. Escaler, CFD study on particle-to-fluid heat transfer in fixed bed reactors: convective heat transfer at low and high pressure, in: The John Bridgwater Symposium: “Shaping the Future of Chemical Engineering”, Chem. Eng. Sci. 61 (2006) 4341–4353, https://doi.org/10.1016/j.ces.2006.02.011. [10] M. Nijemeisland, A.G. Dixon, CFD study of fluid flow and wall heat transfer in a fixed bed of spheres, AIChE J. 50 (2004) 906–921, https://doi.org/10. 1002/aic.10089. [11] W. van Antwerpen, C.G. du Toit, P.G. Rousseau, A review of correlations to model the packing structure and effective thermal conductivity in packed beds of mono-sized spherical particles, Nucl. Eng. Des. 240 (2010) 1803–1818, https://doi.org/10.1016/j.nucengdes.2010.03.009. [12] W.L. Vargas, J.J. McCarthy, Conductivity of granular media with stagnant interstitial fluids via thermal particle dynamics simulation, Int. J. Heat Mass Transf. 45 (2002) 4847–4856, https://doi.org/10.1016/S0017-9310(02)00175-8. [13] W.L. Vargas, J.J. McCarthy, Heat conduction in granular materials, AIChE J. 47 (2001) 1052–1059, https://doi.org/10.1002/aic.690470511. [14] A.S. Mujumdar, Handbook of Industrial Drying, 2nd ed., CRC Press, Boca Raton, Florida, 1995. [15] R.G. Szafran, A. Kmiec, CFD modeling of heat and mass transfer in a spouted bed dryer, Ind. Eng. Chem. Res. 43 (2004) 1113–1124, https://doi.org/10. 1021/ie0305824. [16] B. Chaudhuri, F.J. Muzzio, M.S. Tomassone, Modeling of heat transfer in granular flow in rotating vessels, Chem. Eng. Sci. 61 (2006) 6348–6360, https:// doi.org/10.1016/j.ces.2006.05.034.

22

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

[17] H.N. Emady, K.V. Anderson, W.G. Borghard, F.J. Muzzio, B.J. Glasser, A. Cuitino, Prediction of conductive heating time scales of particles in a rotary drum, Chem. Eng. Sci. 152 (2016) 45–54, https://doi.org/10.1016/j.ces.2016.05.022. [18] D. Shi, W.L. Vargas, J.J. McCarthy, Heat transfer in rotary kilns with interstitial gases, Chem. Eng. Sci. 63 (2008) 4506–4516, https://doi.org/10.1016/j. ces.2008.06.006. [19] A.B. Morris, Z. Ma, S. Pannala, C.M. Hrenya, Simulations of heat transfer to solid particles flowing through an array of heated tubes, Sol. Energy 130 (2016) 101–115, https://doi.org/10.1016/j.solener.2016.01.033. [20] C. Dan, A. Wachs, Direct numerical simulation of particulate flow with heat transfer, in: 7th World Conference on Experimental Heat Transfer, Fluid Mechanics and Thermodynamics (ExHFT-7), Krakow and the Conference on Modelling Fluid Flow, CMFF ’09, Budapest, Int. J. Heat Fluid Flow 31 (2010) 1050–1057, https://doi.org/10.1016/j.ijheatfluidflow.2010.07.007. [21] N.G. Deen, S.H.L. Kriebitzsch, M.A. van der Hoef, J.A.M. Kuipers, Direct numerical simulation of flow and heat transfer in dense fluid–particle systems, Chem. Eng. Sci. 81 (2012) 329–344, https://doi.org/10.1016/j.ces.2012.06.055. [22] Z.-G. Feng, E.E. Michaelides, Heat transfer in particulate flows with direct numerical simulation (DNS), Int. J. Heat Mass Transf. 52 (2009) 777–786, https://doi.org/10.1016/j.ijheatmasstransfer.2008.07.023. [23] Z.-G. Feng, E.E. Michaelides, Inclusion of heat transfer computations for particle laden flows, Phys. Fluids 20 (2008) 040604, https://doi.org/10.1063/1. 2911022. [24] B. Kravets, H. Kruggel-Emden, Investigation of local heat transfer in random particle packings by a fully resolved LBM-approach, Powder Technol. 318 (2017) 293–305, https://doi.org/10.1016/j.powtec.2017.05.039. [25] B. Sun, S. Tenneti, S. Subramaniam, Modeling average gas–solid heat transfer using particle-resolved direct numerical simulation, Int. J. Heat Mass Transf. 86 (2015) 898–913, https://doi.org/10.1016/j.ijheatmasstransfer.2015.03.046. [26] H. Tavassoli, E.A.J.F. Peters, J.A.M. Kuipers, Direct numerical simulation of fluid–particle heat transfer in fixed random arrays of non-spherical particles, Chem. Eng. Sci. 129 (2015) 42–48, https://doi.org/10.1016/j.ces.2015.02.024. [27] J. Feng, H.H. Hu, D.D. Joseph, Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid Part 1. Sedimentation, J. Fluid Mech. 261 (1994) 95–134, https://doi.org/10.1017/S0022112094000285. [28] Z.-G. Feng, E.E. Michaelides, A numerical study on the transient heat transfer from a sphere at high Reynolds and Peclet numbers, Int. J. Heat Mass Transf. 43 (2000) 219–229, https://doi.org/10.1016/S0017-9310(99)00133-7. [29] N.G. Deen, E.A.J.F. Peters, J.T. Padding, J.A.M. Kuipers, Review of direct numerical simulation of fluid–particle mass, momentum and heat transfer in dense gas–solid flows, Chem. Eng. Sci. 116 (2014) 710–724, https://doi.org/10.1016/j.ces.2014.05.039. [30] E.A. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yusof, Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations, J. Comput. Phys. 161 (2000) 35–60, https://doi.org/10.1006/jcph.2000.6484. [31] S. Tenneti, S. Subramaniam, Particle-resolved direct numerical simulation for gas–solid flow model development, Annu. Rev. Fluid Mech. 46 (2014) 199–230, https://doi.org/10.1146/annurev-fluid-010313-141344. [32] M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, J. Comput. Phys. 209 (2005) 448–476, https:// doi.org/10.1016/j.jcp.2005.03.017. [33] A.J.C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation, J. Fluid Mech. 271 (1994) 285–309, https://doi.org/10.1017/S0022112094001771. [34] A.J.C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 2. Numerical results, J. Fluid Mech. 271 (1994) 311–339, https://doi.org/10.1017/S0022112094001783. [35] B. Wachmann, W. Kalthoff, S. Schwarzer, H.J. Herrmann, Collective drag and sedimentation: comparison of simulation and experiment in two and three dimensions, Granul. Matter 1 (1998) 75–82, https://doi.org/10.1007/s100350050012. [36] X. He, L.-S. Luo, Lattice Boltzmann model for the incompressible Navier–Stokes equation, J. Stat. Phys. 88 (1997) 927–944, https://doi.org/10.1023/B: JOSS.0000015179.12689.e4. [37] X. He, L.-S. Luo, Theory of the lattice Boltzmann method: from the Boltzmann equation to the lattice Boltzmann equation, Phys. Rev. E 56 (1997) 6811–6817, https://doi.org/10.1103/PhysRevE.56.6811. [38] D. Yu, S.S. Girimaji, A.J.C. Ladd, Revised moment propagation method for scalar transport, Phys. Rev. E 78 (2008) 056706, https://doi.org/10.1103/ PhysRevE.78.056706. [39] R.M.H. Merks, A.G. Hoekstra, P.M.A. Sloot, The moment propagation method for advection–diffusion in the lattice Boltzmann method: validation and Péclet number limits, J. Comput. Phys. 183 (2002) 563–576, https://doi.org/10.1006/jcph.2002.7209. [40] C.P. Lowe, D. Frenkel, The super long-time decay of velocity fluctuations in a two-dimensional fluid, Phys. A, Stat. Mech. Appl. 220 (1995) 251–260, https://doi.org/10.1016/0378-4371(95)00208-O. [41] L. Li, R. Mei, J.F. Klausner, Lattice Boltzmann models for the convection-diffusion equation: D2Q5 vs D2Q9, Int. J. Heat Mass Transf. 108 (2017) 41–62, https://doi.org/10.1016/j.ijheatmasstransfer.2016.11.092. [42] H. Yoshida, M. Nagaoka, Multiple-relaxation-time lattice Boltzmann model for the convection and anisotropic diffusion equation, J. Comput. Phys. 229 (2010) 7774–7795, https://doi.org/10.1016/j.jcp.2010.06.037. [43] D. Contrino, P. Lallemand, P. Asinari, L.-S. Luo, Lattice-Boltzmann simulations of the thermally driven 2D square cavity at high Rayleigh numbers, J. Comput. Phys. 275 (2014) 257–272, https://doi.org/10.1016/j.jcp.2014.06.047. [44] Z. Hashemi, O. Abouali, R. Kamali, Thermal three-dimensional Lattice Boltzmann simulations of suspended solid particles in microchannels, Int. J. Heat Mass Transf. 65 (2013) 235–243, https://doi.org/10.1016/j.ijheatmasstransfer.2013.06.005. [45] Z. Guo, C. Zheng, B. Shi, T.S. Zhao, Thermal lattice Boltzmann equation for low Mach number flows: decoupling model, Phys. Rev. E 75 (2007) 036704, https://doi.org/10.1103/PhysRevE.75.036704. [46] X. He, S. Chen, G.D. Doolen, A novel thermal model for the lattice Boltzmann method in incompressible limit, J. Comput. Phys. 146 (1998) 282–300, https://doi.org/10.1006/jcph.1998.6057. [47] H. Kruggel-Emden, B. Kravets, M.K. Suryanarayana, R. Jasevicius, Direct numerical simulation of coupled fluid flow and heat transfer for single particles and particle packings by a LBM-approach, Powder Technol. 294 (2016) 236–251. [48] Y. Peng, C. Shu, Y.T. Chew, Simplified thermal lattice Boltzmann model for incompressible thermal flows, Phys. Rev. E 68 (2003) 026701, https:// doi.org/10.1103/PhysRevE.68.026701. [49] L. Wang, D. Koch, X. Yin, C. Cohen, Hydrodynamic diffusion and mass transfer across a sheared suspension of neutrally buoyant spheres, Phys. Fluids 21 (2009) 033303, https://doi.org/10.1063/1.3098446. [50] B. Metzger, O. Rahli, X. Yin, Heat transfer across sheared suspensions: role of the shear-induced diffusion, J. Fluid Mech. 724 (2013) 527–552, https:// doi.org/10.1017/jfm.2013.173. [51] D. Yu, A.J.C. Ladd, A numerical simulation method for dissolution in porous and fractured media, J. Comput. Phys. 229 (2010) 6450–6465, https:// doi.org/10.1016/j.jcp.2010.05.005. [52] P. Lallemand, L.-S. Luo, Hybrid finite-difference thermal lattice Boltzmann equation, Int. J. Mod. Phys. B 17 (2003) 41–47, https://doi.org/10.1142/ S0217979203017060. [53] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: acoustic and thermal properties in two and three dimensions, Phys. Rev. E 68 (2003) 036706, https://doi.org/10.1103/PhysRevE.68.036706.

A.M. Lattanzi et al. / Journal of Computational Physics: X 1 (2019) 100007

23

[54] P. Salamon, D. Fernàndez-Garcia, J.J. Gómez-Hernández, A review and numerical assessment of the random walk particle tracking method, J. Contam. Hydrol. 87 (2006) 277–305, https://doi.org/10.1016/j.jconhyd.2006.05.005. [55] G. Srinivasan, D.M. Tartakovsky, M. Dentz, H. Viswanathan, B. Berkowitz, B.A. Robinson, Random walk particle tracking simulations of non-Fickian transport in heterogeneous media, J. Comput. Phys. 229 (2010) 4304–4314, https://doi.org/10.1016/j.jcp.2010.02.014. [56] H. Hoteit, R. Mose, A. Younes, F. Lehmann, P. Ackerer, Three-dimensional modeling of mass transfer in porous media using the mixed hybrid finite elements and the random-walk methods, Math. Geol. 34 (2002) 435–456, https://doi.org/10.1023/A:1015083111971. [57] Y. Sun, C.-H. Park, G. Pichot, J. Taron, Random walk particle tracking, in: Thermo-Hydro-Mechanical-Chemical Processes in Fractured Porous Media: Modelling and Benchmarking: Closed-Form Solutions, Springer, Switzerland, 2014, pp. 153–158. [58] P. Szymczak, A.J.C. Ladd, Boundary conditions for stochastic solutions of the convection-diffusion equation, Phys. Rev. E 68 (2003) 036704, https:// doi.org/10.1103/PhysRevE.68.036704. [59] W. Kinzelbach, The random walk method in pollutant transport simulation, in: Groundwater Flow and Quality Modelling, in: NATO ASI Ser., Springer, Dordrecht, 1988, pp. 227–245. [60] S. Bekri, P.M. Adler, Dispersion in multiphase flow through porous media, Int. J. Multiph. Flow 28 (2002) 665–697, https://doi.org/10.1016/S03019322(01)00089-1. [61] E.M. LaBolle, G.E. Fogg, A.F.B. Tompson, Random-walk simulation of transport in heterogeneous porous media: local mass-conservation problem and implementation methods, Water Resour. Res. 32 (1996) 583–593, https://doi.org/10.1029/95WR03528. [62] M. Tomadakis, S. Sotirchos, Transport through random arrays of conductive cylinders dispersed in a conductive matrix, J. Chem. Phys. 104 (1996) 6893–6900, https://doi.org/10.1063/1.471356. [63] H. Tavassoli, S.H.L. Kriebitzsch, M.A. van der Hoef, E.A.J.F. Peters, J.A.M. Kuipers, Direct numerical simulation of particulate flow with heat transfer, Int. J. Multiph. Flow 57 (2013) 29–37, https://doi.org/10.1016/j.ijmultiphaseflow.2013.06.009. [64] A.J.C. Ladd, R. Verberg, Lattice-Boltzmann simulations of particle–fluid suspensions, J. Stat. Phys. 104 (2001) 1191–1251, https://doi.org/10.1023/A: 1010414013942. [65] Z. Guo, T.S. Zhao, Y. Shi, Physical symmetry, spatial accuracy, and relaxation time of the lattice Boltzmann equation for microgas flows, J. Appl. Phys. 99 (2006) 074903, https://doi.org/10.1063/1.2185839. [66] X. Shan, X. He, Discretization of the velocity space in the solution of the Boltzmann equation, Phys. Rev. Lett. 80 (1998) 65–68, https://doi.org/10.1103/ PhysRevLett.80.65. [67] K. Itô, On stochastic differential equations, in: Memoirs of the American Mathematical Society, vol. 4, 1951, pp. 289–302. [68] C.W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry, and the Natural Sciences, 2nd ed., Springer-Verlag, New York, 1986. [69] P.E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Stochastic Modelling and Applied Probability, Springer-Verlag, Berlin, Heidelberg, 1992. [70] J.B. Avalos, A.D. Mackie, Dissipative particle dynamics with energy conservation, Europhys. Lett. 40 (1997) 141, https://doi.org/10.1209/epl/i1997-004366. [71] A.D. Mackie, J.B. Avalos, V. Navas, Dissipative particle dynamics with energy conservation: modelling of heat flow, Phys. Chem. Chem. Phys. 1 (1999) 2039–2049, https://doi.org/10.1039/A809502G. [72] M. Ripoll, P. Español, M.H. Ernst, Dissipative particle dynamics with energy conservation: heat conduction, Int. J. Mod. Phys. C 09 (1998) 1329–1338, https://doi.org/10.1142/S0129183198001205. [73] C.K. Aidun, J.R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (2010) 439–472, https://doi.org/10.1146/annurevfluid-121108-145519. [74] D.J. Jeffrey, Conduction through a random suspension of spheres, Proc. R. Soc., Math. Phys. Eng. Sci. 335 (1973) 355–367, https://doi.org/10.1098/rspa. 1973.0130. [75] C. Miller, S. Torquato, Effective conductivity of hard-sphere dispersions, J. Appl. Phys. 68 (1990) 5486–5493, https://doi.org/10.1063/1.347007. [76] S. Patankar, Numerical Heat Transfer and Fluid Flow, 1st ed., CRC Press, Boca Raton, Florida, 1980. [77] J. Thibault, S. Bergeron, H.W. Bonin, On finite-difference solutions of the heat equation in spherical coordinates, Numer. Heat Transf. 12 (1987) 457–474, https://doi.org/10.1080/10407788708913597. [78] A.N. Ford Versypt, R.D. Braatz, Analysis of finite difference discretization schemes for diffusion in spheres with variable diffusivity, Comput. Chem. Eng. 71 (2014) 241–252, https://doi.org/10.1016/j.compchemeng.2014.05.022. [79] T. Forgber, B. Mohan, C. Kloss, S. Radl, Heat transfer rates in sheared beds of inertial particles at high Biot numbers, Granul. Matter 19 (2017) 14, https://doi.org/10.1007/s10035-016-0695-0. [80] T. Oschmann, M. Schiemann, H. Kruggel-Emden, Development and verification of a resolved 3D inner particle heat transfer model for the discrete element method (DEM), Powder Technol. 291 (2016) 392–407, https://doi.org/10.1016/j.powtec.2015.12.008. [81] W.E. Ranz, W.R. Marshall, Evaporation from drops, Chem. Eng. Prog. 48 (1952) 141. [82] S. Whitaker, Forced convection heat transfer correlations for flow in pipes, past flat plates, single cylinders, single spheres, and for flow in packed beds and tube bundles, AIChE J. 18 (1972) 361–371, https://doi.org/10.1002/aic.690180219. [83] A. Richter, P.A. Nikrityuk, Drag forces and heat transfer coefficients for spherical, cuboidal and ellipsoidal particles in cross flow at sub-critical Reynolds numbers, Int. J. Heat Mass Transf. 55 (2012) 1343–1354, https://doi.org/10.1016/j.ijheatmasstransfer.2011.09.005. [84] J.R. Grace, F. Taghipour, Verification and validation of CFD models and dynamic similarity for fluidized beds, Powder Technol. 139 (2004) 99–110, https://doi.org/10.1016/j.powtec.2003.10.006. [85] R. Maier, R. Bernard, D. Grunau, Boundary conditions for the lattice Boltzmann method, Phys. Fluids 8 (1996) 1788–1801, https://doi.org/10.1063/1. 868961. [86] Z. Yang, Lattice Boltzmann outflow treatments: convective conditions and others, in: Special Issue on Mesoscopic Methods in Engineering and Science, ICMMES-2010, Edmonton, Canada, Comput. Math. Appl. 65 (2013) 160–171, https://doi.org/10.1016/j.camwa.2012.11.012. [87] M. Nijemeisland, Influences of Catalyst Particle Geometry on Fixed Bed Reactor Near-Wall Heat Transfer Using CFD, Worcester Polytechnic Institute, Worcester, MA, USA, 2000.