One of the challenges in numerical modelling of technological processes in granular media is an accurate and efficient handling of random geometry. The issue to be considered is especially challenging when the coupling of physico-chemical processes at different scales occurs in such random media. Thus, the idea of multi-scale modelling of coupled processes takes a considerable place in fluid mechanics and chemical engineering as an approach of solving problems not only in academic-type situations but also in industrial processes. One of the aims is to assure sufficient accuracy of averaged description applied in macro-scale simulations, thus a detailed investigation in smaller scale is highly desirable, for example to obtain closure relations. One of the measurable effects of such investigations are averaged relationships utilized in 1D/2D models provided by the use of conservation laws closed by empirical correlations (Polesek-Karczewska et al. 2015); very often those models were experimentally validated or accordingly fitted to meet the macro-scale needs. An example, taken into consideration in the mentioned scale, is the thermal processing of fossil fuel (like coal). Starting from the micro-scale accounting for a pore level of a single particle, the flow and heat transfer through a solid matrix have to be modelled along with chemical kinetics. At meso-scale as the next level of description, reactive flow thermomechanics in domain with volume of grains dependent on temperature, plastic deformation, and the computation of the resulting stresses in a packed bed of solid fuel particles are considered. Last but not least, as the results at this level are interesting for industry, the description of averaged fields in a macro-scale geometry relevant to the industry has to be accounted for.

The detailed computations at the meso-scale with the use of a proper model of the micro-scale phenomena (at the pore level) can accurately predict macroscopic fields in the geometry of representative element of volume (REV), for regular and random geometry, see (Xu et al. 2017). In (Asinari et al. 2007), the Authors analyse phenomena observed in a solid oxide fuel cell (SOFC) at the level of a single pore in a REV being in fact a randomly created porous medium whose topology is closely related to that of SOFC. Aside of fluid flow in a complex geometry, problems of chemical reactions and species transport are highlighted there. In result, it is argued in (Asinari et al. 2007) that optimization of the fluid element paths would produce an increase in macroscopic performance. An example of modelling of fluid flow in granular media where various kinds of disorder is obtained for a given porosity is presented by Dai et al. (2019); in their work, phenomena of invading fluid injection are modelled with special attention given to various deviation of local porosity. Heat transport in complex geometry was also analyzed by Wang and Pan (2008), focused on the REV domain; discussion is mainly concentrated on the comparison of computational results with theoretical predictions.

Whenever available, theoretical or experimental approximations are often used in case of macro scale modelling to provide the necessary closure relationships such as the effective transport coefficients. To improve the accuracy, those relations are still subject of research. Thus, multiscale studies recently encountered in the literature present detailed investigations, especially for cases of coupled phenomena in complex geometries (see Xu et al. 2017). Another example is (Miura and Silveston 1980) where the Authors present valuable results regarding pore evolution during thermal processing of coal. They account for thermomechanics of flow in pores at micro- and macro-scales based on a comparison with experimental results. The paper concludes with a validated pore development model for coking coals. Numerical analysis at the grain scale was presented by Kayser et al. (2006) who investigated the influence of the size distribution of particles in overall flow of gas through a packed bed of coal grains. The Authors checked how particle size distribution influences the numerically predicted pressure drop and compared their results with empirical correlation based on the Ergun equation. Similar investigation regarding heat transfer phenomenon, including the experimental and numerical setups along with results, are presented in (Polesek-Karczewska 2003). In that work, the effective thermal conductivity of packed bed of spheres is estimated with coupled phenomena taken into account. However, the numerical predicitions underestimated the experimentally determined effective thermal conductivity since the radiation effects were not accounted for in modelling.

In the cases where process conditions prevent from making detailed measurements, the computational modelling has become an important player. Such investigations can be aided by so-called structure models proposed for estimation of effective thermal conductivity (or other transport coefficients) for geometries where the structure is ordered (like packed beds, successively arranged layers, etc.). For such cases, depending on dominant heat transfer direction, perpendicular or parallel with respect to the main flow direction, effective transport coefficients are calculated with the use of known correlations, see (Wang et al. 2019). Recent investigations on effective transport coefficients concern geometrically complex cases such as granular and fibrous porous media, see (Xu et al. 2017; Wang and Pan 2008; Askari et al. 2015; Demuth et al. 2014) for some known models. Wang and Pan (2008) summarised some of available solutions for the effective thermal conductivities obtained for two-component systems, like the effective medium theory (EMT). They also analysed theoretical approximations of the effective conductivity in materials with microstructure (porous, fibrous, granular, etc.). In the work of Gong et al. (2014), aside of a review of models presented in (Wang and Pan 2008), an example of further developments was presented. The authors proposed a novel EMT for components (small spheres) which are uniformly dispersed in a fluid or solid matrix. Towards providing a detailed description of heat transfer in complex geometries, a model with three materials as components taken into account was developed in a work of Thiele et al. (2014). The authors describe change of effective thermal conductivity of granular media in the form of capsules, basing on the heat transfer analysis with the use of finite element method. They check the influence of obstacles layout and core/shell/matrix thermal conductivities. Also, Thiele et al. (2014) analysed the influence of the size distribution of spherical particles on the effective thermal conductivity. The authors presented a model for prediction of the thermal conductivity with numerical uncertainty and identified conditions under which heat transport coefficient remains of satisfactory accuracy. One of the conclusions was that the layout of the second material (solid grains) had irrelevant influence on effective heat transfer (where the case is a non-reactive flow, in contrast to Asinari et al. 2007). The work of Dai et al. (2019) is another example of granular media heat transfer. Their work contains detailed description of the numerical tool where important phenomena are modelled at the grain scale. The important feature is the experimental validation of the constructed model. In the case, where more detailed analysis is required, one can use empirical correlations for physico-chemical phenomena typical of the process. The importance of the meso-scale simulations gives the potential to determine effective transport coefficient appropriate for averaged modelling. In a recent work, Askari et al. (2015) model heat transfer at the pore scale level with a number of relevant phenomena taken into account like the contact resistance and roughness of grains. The authors clearly showed the importance of modelling surface characterizing parameters on the effective thermal conductivity. In the mentioned work by Polesek-Karczewska (2003) indicates underestimated effective thermal conductivities from empirical correlations in a thermal processing of solid grains (balls). In conclusion, the author provides a possible explanation of discrepancies between experiment and used model.

To obtain an exhaustive description of technological processes, one has to undertake the interdisciplinary research of the physico-chemical phenomena in porous media. For the purpose of modelling the fluid thermomechanics, the lattice Boltzmann method (LBM) is utilized. The method could be introduced starting from the Boltzmann equation where appropriate discretization is introduced, see (Aidun and Clausen 2010; Succi 2001; Ya-Ling et al. 2019) for details. LBM is especially suitable for modelling the flow of nearly incompressible fluid in a simple (Xu et al. 2017; Guo et al. 2002) as well as complex (Wang et al. 2019; Wang and Ning 2008; Chai et al. 2010) geometries. Additional phenomena are also modelled, like heat transfer (He et al. 1998; Wang et al. 2007a, b) and chemical species evolution (Xu et al. 2017; Yamamoto et al. 2002; Di Rienzo et al. 2012). As the process is modelled in the meso-scale (i.e. at the grain level), the micro-scale approach is utilized for modelling phenomena at the pores of grains. Utilized are simplified models for chemical reactions and transport inside the coal particles. Regarding estimation of effective heat conductivity with the use of LBM in a complex geometry, some examples are known in the literature. So in (Wang and Ning 2008; Wang et al. 2007a, b) LBM is applied to compute effective thermal conductivity in randomly generated media, accordingly simplified for the studied case (no fluid flow was considered). Also in (Demuth et al. 2014) the LBM was used to compute the effective thermal conductivity in a representative geometry of regular granular media. The authors compared results from the LBM with implementations of finite volume method (FVM) in a commercial software as well as an in-house code. In both works, the comparison of LBM results with other estimations (different numerical schemes, references and empirical correlations) substantiate the usefulness of LBM in modelling of heat transfer in complex geometry.

The long term motivation for the present work is to provide detailed description of the coking process. Coking plants are used in the coal-processing industry, to obtain cleaner coal (coke) with other chemical products known as coal gas and tar. The approach of averaged level modelling is the common way for simulating devolatilisation, usually solving the most important phenomena. Some examples of such simulations are Di Blasi (2008) and Polesek-Karczewska et al. (2013). Hence, for the purpose of present analysis, we further develop unified numerical tool utilizing LBM for flow thermomechanics, (Grucelski and Pozorski 2013,2015,2017).

Here, an effective heat conductivity value is estimated for a bed of coal grains (cylinders in 2D case). The work is intended as a continuation to a preliminary investigation presented earlier in Grucelski (2016). Still, description regarding chosen schemes and models of LBM (for the fluid thermomechanics and the evolution chemical species with global chemical model of devolatilisation) is given in Sect. 2. Results of heat transfer in layers with different heat conductivity (chosen as a benchmark) along with comparison with other works and analytical solutions are presented in Sect. 3. The main findings of the paper, i.e. the resulting estimation of effective thermal conductivity along with other results from LBM calculations in granular media are reported in Sect. 4. In this section, one additional test case is presented for non-reactive flow in the same computational domain.

The representative element of volume (REV) is often chosen for modelling of cases where simulation on a larger scale is prohibitively expensive and the influence of the geometry details cannot be further brushed aside (as it is in the averaged modelling). For modelled case, the size of computational domain (here REV) is larger than diameter of coal particle (micro-scale) and smaller than a typical industrial devices (macro-scale).

2.1 Governing equations

For a more intuitive and versatile description of the process, the conservation equations are briefly presented at the outset. During the simulation a low-Mach number reactive flow is modelled, which requires the solution of mass, momentum and energy conservation equations:

where the transport of species (actually a mixture), although modelled, will not be reported in details (Grucelski and Pozorski 2017; Grucelski 2016); chemical specie release directly impacts the average density ρ and the ratio of volatiles left in the grain volume. Briefly, the gaseous products released to the volume of fluid, depend on volatiles left in grain and this aspect of the model is implemented through the boundary scheme locally affecting the momentum. The kinematic viscosity ν is assumed constant in the simulation and the flow is nearly incompressible. The Boussinesq approximation is used to model natural convection; due to linear dependency of density on temperature, the mass force takes the form F = gβ (θ − θ₀), see (Guo et al. 2002), where θ and θ₀ are the local and reference fluid temperatures, respectively; β is the fluid thermal expansion parameter (here \(\beta = 10^{ - 3} \left[ {K^{ - 1} } \right]\)). In Eq. (3), the compression work done by the pressure is neglected and viscous heat dissipation is taken into consideration, see (He et al. 1998) for details; heat transfer coefficient takes the form α = k_eff/ρc_p. The quantity \(Q_{\theta }\) in Eq. (3) pertain to sink of heat, due to gas release; it is detailed in Sect. 2.4.

2.2 Lattice Boltzmann method

Here only the salient features of LBM are presented; an exhaustive description of the method can be found in Aidun and Clausen (2010), Succi (2001) and He et al. (1998). The Lattice Boltzmann equation, discretised in time, space (by lattice) and a set of advection velocities (e_i,i = 0,…,M − 1) on a regular square lattice, describes important physical fields in terms of distribution functions. Let us use a symbol ψ = [f, g, φ] being a vector of distribution functions of fields of interest, here the fluid density, temperature and chemical species concentration. In the presented simulation approach, two dimensional (2D) lattice is considered, so the chosen discretisation of advection velocity in this work is D2Q9 with M = 9 advection velocities of distribution function.

The distribution function evolves as (see Wang et al. 2007a, b):

where, ψ_i^eq stands for the equilibrium parts of distribution function (DF) at (r,t) and \(Q_{i}^{\psi }\) represents a source term. Relaxation time (here in non-dimensional form) τ_ψ depends on physical transport coefficients, characterising modelled phenomena. In this work, the Single Relaxation Time (SRT) as one of the recognisable schemes is implemented in the numerical tool, see Succi (2001) and Ya-Ling et al. (2019). The scheme is one of the well known, exhaustively validated and remains widely utilized in modelling of physicochemical phenomena (see Di Rienzo et al. 2012), very often in complex geometry, like drying in porous media (see Zachariah et al. 2019). Usually, f is used to track the density and velocity with evolution Eq. (4). In a similar fashion, the evolution of temperature and chemical species concentration is tracked with additional density distribution functions. In case of heat transfer g is used in Eq. (4). In that form, g is known as internal energy density distribution function (IEDDF, see He et al. 1998) with appropriate evolution equation. Chemical specie evolution (in fact a mixture of gaseous products of devolatilisation) is described using another distribution function φ (in fluid domain only).

The equilibrium distribution function for every accounted evolution equation has a similar form:

where, c = δx/δt is the lattice speed of advection of distribution functions, δt is a time step and u is the local fluid velocity and Ω_i are the weight coefficients; for D2Q9, \({\Omega }_{i = 0} = 4/9\), \({\Omega }_{i = \{1, \ldots ,4\}} = 1/9\), \({\Omega }_{{i = \left\{ {5, \ldots , 8} \right\}}} = 1/36\). In Eq. (5), coefficients A_i to D_i vary depending on phenomena modelled as well as advection directions and scheme of discretization utilized. The exact form of Eqs. (4)–(5) regarding coefficients for fluid flow and heat transfer LBM modelling are presented by Wang et al. (2007a, b). Regarding the fluid flow, in a generic form considered here, symbol \(\psi_{i}^{eq} = f_{i}^{eq}\) is used and the coefficients (direction independent) are A = 1, B = 3, C = 4.5, D =  − 1.5 together with β = ρ. In the case of heat transfer modelling with use of IEDDF in Eq. (5) β = θ the coefficients take the form (Grucelski 2016):

In contrast, Wang et al. (2007a, b) use a simplified equilibrium distribution Eq. (5) for IEDDF calculations with zero order polynomial and first order polynomial was used by Wang and Ning (2008). In the present work, natural convection is simulated as well. One of the main results is accurate calculation of a second moment of IEDDF, that is the heat flux, so algorithm operates on equilibrium equation in the full form.

The equilibrium density functions for concentration of chemical species (first introduced by Di Rienzo et al. (2012)) for D2Q9 was presented by Grucelski and Pozorski (2015), with:

Modifications proposed by DiRienzo et al. (2012) account for varied density field in the volume of fluid. They introduce a parameter, also in Eq. (5), in the form ψ = ρ∗ /ρ to be understood as a ratio of the minimum density of fluid in the entire domain and the density at a given lattice node.

In the case of fluid flow with viscosity ν, the relaxation time in Eq. (4) takes the form (Aidun and Clausen 2010):

here, \(c_{s} = {\raise0.7ex\hbox{$c$} \!\mathord{\left/ {\vphantom {c {\sqrt 3 }}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${\sqrt 3 }$}}\) is the speed of sound in the simulation. In case of the IEDDF the relaxation time τ_ψ = τ_k,m written both for fluid (m = f, gas) and solid (m = s, coal grains) in the computational domain (m ∈ {s,f}) takes a form:

where C_p,_m and λ_m stands for heat capacity and heat conductivity of m. Concerning the transport of species, in presented approach, a single specie (actually, a mixture) is tracked with distribution function φ (in fluid nodes) and non-dimensional relaxation time:

with D_ϕ being a diffusivity coefficient of the mixture.

Interesting macroscopic variables, that is density ρ, velocity u, temperature θ, heat flux q and chemical specie concentration Y, in relevant nodes are obtained by:

that is by integration of respective distribution functions (Wang and Ning 2008; Di Rienzo et al. 2012).

2.3 Natural convection

To account for fluid density variations due to local temperature difference in the gravity field with acceleration g the scheme of Guo et al. (2002) is implemented. The scheme is chosen as widely recognisable in the literature and validated by the researchers. Moreover, in Ya-Ling et al. (2019) the Authors argue that Guo and others scheme of the body force remains basically the same. Equations describing transport of the fluid and heat are coupled by the following source term (natural convection phenomenon):

Here β is the thermal expansion coefficient of air, θ and \(\langle\theta\rangle\) are local and averaged temperatures, respectively (cf. Grucelski 2016).

2.4 Pyrolysis gas release, source terms consideration

In a recent work (see Grucelski and Pozorski 2017) detailed description of chemical species along with simplified tar release is presented; here, an evolution of a mixture of gases (actually single compound) is modelled. For detailed analysis one has to consider the evolution of every specie by modelling the evolution of appropriate DF in the LBM approach. This stands in opposition to the model where the mixture is considered.

At the moment, a detailed analysis is not meant to be considered here; for the purpose of 2D meso-scale simulation a model is simplified (single mixture approach, still widely used for engineering purposes). An example of such model available in the literature is described in details by Postrzednik (1994). Here only a short presentation of empirical correlations is given. For a grain G the released gaseous product is described by:

here ε is fluid volume fraction (porosity), V_G and \(\varepsilon_{G}^{vol}\) are the volume and volume fraction of volatiles of grain G, respectively; \(\rho_{ }\) is the averaged density of fluid in the domain.

The term \(f_{V} \left( {\varepsilon_{G}^{vol} , \theta } \right)\), namely the release rate in Eq. (13), takes the form of Arrhenius equation:

Some modifications are introduced regarding activation energy and pre-exponential factor, \(b\left( {\varepsilon_{G}^{vol} } \right)\) and \(a\left( {\varepsilon_{G}^{vol} } \right)\), respectively; for further details see Grucelski (2016).

Equation (13) has form of empirical function and has been proposed by Postrzednik (1994):

The actual values of coefficients \(a_{Y} - g_{Y}\) and the characteristic temperature \(\theta_{Y}\) depend on the variety of coal, see Postrzednik (1994) for sample values. The first phase of chemical specie release occurs at about θ = 370 K and the second is considered till about θ = 700 K. In this work we use the value for Budryk coal (a coal mine in the Silesia region, Poland) is used, here θ_Y = 573 K. In brief,\(g_{V} \left( \theta \right)\) corresponds to change of density due to heating and coupled phenomena like gas release.

The reaction rate \(\omega_{V}\), resulting from Eq. (13), is next used to calculate the source terms. In case of chemical species transport, the yield of gaseous products (source term) is uniformly distributed on the surface of corresponding obstacle. Thus for a grain G, whose surface is discretised by a set of lattice points GS, the source term in Eq. (4) takes the form:

where \(Q_{{S_{r} }}\) is a source of chemical specie at a surface node of solid particle G utilized by the boundary scheme for density distribution function at the surface of grain. After summing over the surface nodes of the source term in Eq. (13) one gets, see Yamamoto et al. (2002):

here, (u₀)_LBM is taken in lattice units as a maximal velocity in the domain, next dimensionalised by the speed of sound in air through u₀ = c(u₀)_LBM; the molar mass use in the numerical modeling M = 150 kg/kmol. The mixture is composed of 10% of CO₂, 9.4% of H₂O, 10.3% of CO, 11% of hydrocarbons, 44% of nonvolatile C, 16% of tar and other small amounts of products (for further details see Serio et al. 1981). To accurately simulate heat transfer in reacting granular media, the source of heat is also calculated in proposed description. Here the heat source term results from similar correlation as for mass: \(\left( {Q/c_{p} \rho } \right)\omega_{LBM} = \mathop \sum \nolimits_{r \in G} \mathop \sum \limits_{i} \Omega_{i} Q_{G}\) whereas Q_G is uniformly distributed among all solid nodes of grain G.

2.5 Some other details of numerical procedure

The details regarding the computational domain are presented by Grucelski (2016) and here some differences and the basic information are summarized. The domain (discretized by 200 × 400 nodes) is constructed as a REV element with a random layout of obstacles (grains mean diameter is d = 4.8 mm). The REV size is of about 80 obstacles of diameter d per diagonal of the domain. In the algorithm, the positions and radii of subsequent obstacles are randomly set until an input porosity (solid volume fraction) ε = 0.51 is reached. Although the position is randomly generated from the uniform distribution, radius of every grain are chosen from the Gauss distribution. For every generated obstacle (fossil fuel grain), the position is checked to set the minimal distance between neighbours is preserved on a given lattice. Alternative way is presented by Wang et al. (2019), where Monte Carlo method is utilized to obtain disorder index. In our case much higher porosity is needed. In this work, grains are allowed to cross the boundaries (eventually). At the top and bottom boundaries we use a simple periodic condition (for velocity and temperature). We have to point out that the pressure difference, due to free convection in the REV geometry, is not taken into account for the sake of simplicity. As the pressure difference between top/bottom horizontal walls (from free convection of fluid) is unknown, the flow driving force in vertical direction (buoyancy) is neglected. In fact, the REV element is assumed to be located closer to the bottom part of macro-scale geometry. In consequence the periodic condition is used for fluid flow at the bottom and top (horizontal) boundaries (in preliminary study a mixed condition was used, see Grucelski 2016). This condition type is easy to implement; in short, the DDF leaving the computational domain at one of the boundaries are simply transported onto the opposite boundary. In the case of heat transfer, zero flux condition is imposed at horizontal boundaries. Similarly, in the case of fluid flow at the right boundary of the domain, a symmetry plane is used. In LBM, this boundary scheme is realised through assigning to the incoming distribution functions (that are unknown) a value of DDF at boundary node pointing outside the numerical domain, with the symmetry taken into account (locally).

The fluid–solid interface is treated with no-slip schemes: the bounce-back scheme (see Grucelski and Pozorski 2015) for fluid flow and the scheme of He et al. see (1998) for heat transfer. To account for the jump of thermal conductivity at the solid–fluid interface the scheme presented in Li et al. (2014) for the case of interface crossing the lattice at the mid-point between the nodes is introduced in simulation. In the case of the fluid flow, additional source term at the boundary is implemented to account for fluid mass increase due to gas release. As the initial condition zero concentration of the chemical species in the volume of the fluid is imposed. The release is modelled with 0D model, i.e. the gaseous product release and internal energy change are uniformly distributed among nodes (surface and grain volume).

2.6 Formula for effective thermal conductivity

For the case of benchmark the steady state is acquired from LBM simulations of heat conduction (as well as in case of fluid flow).

To calculate an effective value of thermal conductivity, we use (Wang and Ning 2008):

where q_x stands for heat flux in x direction and A is slice area; ∆θ the temperature difference between boundaries where the heat flux is calculated. Although presented results refer to the steady state of simulations, calculated heat conductivity shows negligible deviations and one can say that it is nearly constant. At this stage, the main difficulty is to distinguish areas where for accurate description of heat transfer the thermal dilatation of grains has to be additionally implemented. Nevertheless, for biomass pyrolysis where the change of grains volume has a negligible effect, the phenomena like convection, chemical reactions, etc., play a crucial role in heat transfer and thus in the estimation of the effective heat conductivity. For the benchmark case (detailed in Sect. 3) the relative error is estimated as err = 100 (k_eff − k_ref)/k_ref, where k_ref is calculated from analytical relationship, see Eqs. (18)–(19). The results are shown in Table 1.

For the benchmark case heat conduction in a medium built of two layers is considered. The material is placed in a series or parallel layout. This can be done through implementation of two different sets of boundary conditions in a numerical domain presented in Fig. 1. The numerical tool was previously verified in a number of numerical tests regarding (but not limited to) fluid flow (see Grucelski and Pozorski 2013), heat transfer (see Grucelski and Pozorski 2015) and evolution of chemical compounds (see Grucelski and Pozorski 2017). Presented benchmark is chosen mainly for validation of conduction phenomenon.

Sketch of the domain with boundary conditions for benchmark purpose: a parallel layout, b series layout. Hatching is introduced to distinguish the two homogeneous layers

The modelled materials have different thermal conductivities, k₁ and k₂; different conditions at boundaries determine proper test variants, i.e. parallel and series arrangement, see Fig. 1. After a steady state is reached in the computations, the effective thermal conductivity can be calculated for the domain (cf. Wang et al. 2007a, b). To describe the problem, usually one of the recognised structural models for two-component materials is used; those models were studied theoretically and analytically, see Dul’nev and Sigalova (1967). Although well known, this geometrically simple case is still often chosen to validate numerical calculations of heat transfer. Solutions for cases used in the present paper as well as for other, more complicated models are presented in the literature, see Wang et al. (2007a, b), Dul’nev and Sigalova (1967), Wang et al. (2006). Carson et al. (2005) and references therein. For the parallel layout, see Fig. 1a, the analytical formula for the effective thermal conductivity is:

and for the series layout, Fig. 1b:

Although in Wang et al. (2007a, b) the Authors reported satisfactory accuracy of calculations with the use of simplified form of thermal LBM equations (heat conduction only), here the formulation presented in Sect. 2 is used. As the present paper treats about heat transfer in granular media with account of modelling gaseous species release and its transport in the domain, the equilibrium equations in full form are utilized here.

For the case the mesh of 600 × 600 nodes is used. The relaxation parameter is set to τ_ν = 0.51 for the mass DF and c = 10. These two parameters along with thermal conductivity ratio are enough to calculate the relaxation parameters for heat transfer (see discussion for Fig. 3). The overall results are presented in Table 1. Although the convergence is obtained for almost all cases considered (see further discussion for k₁/k₂ = 10⁴), the overall accuracy is not so good as in other studies [for example, see (Polesek-Karczewska et al. 2015)]. The observed discrepancy between numerical and analytical results increases with the ratio of thermal conductivities, especially for series arrangement: for example if k₁/k₂ = 10⁴ then the error reaches about 22%. The accuracy issues have their source in the selection of the appropriate value of the relaxation time for BGK LBM, where one can show that \(\tau_{k,1}^{ - 1} = f\left( {\tau_{\nu } , \nu , k_{1} } \right)\). In case of used parameters, for the fluid as well as heat transfer \(\tau_{k,1}^{ - 1} = 1.96,\) which is close to the limit of stability of numerical calculations in LBM. In that case, the heat transfer for solid has to be modelled with \(\tau_{k,1}^{ - 1} \sim 0.01\) when k₁/k₂ = 10⁴. Theoretically, in the case of LBM with the SRT approximation of the collision term, it is known that results obtained with \(\tau_{*}^{ - 1} = 0.5 - \delta \tau\) are burdened with a numerical error, increasing inversely proportionally to δτ. In case of the parallel layout of layers, the error increases much slower with increasing k₁/k₂; LBM results compared with empirical relation value, proportional dependency the thermal conductivity of every of the layers is present. For the second geometrical case, that is the series layout Eq. (19), the inversely proportional dependence is used to approximate the value of effective thermal conductivity. This causes a significant contribution of the layers with large heat conductivity; in that case very high k₁/k₂ ratio prevents proper selection of modelling parameters to assure suitable relaxation times in each layer. Finally one will observe rather higher errors in results of LBM numerical simulation.

Figure 2 presents the evolution in time of results from LBM calculations regarding k_eff. Estimated k_eff(t) for k₁/k₂ = 10³ reaches a constant value with error close to 1%; at a steady state for smaller ratios of k₁/k₂ similar error is noticed. For the case of k₁/k₂ = 10⁴ after a long simulation time, the plotted values do not converge to 1; yet they become constant. When the time evolution of results in Fig. 2 are compared for k₁/k₂ = 10⁴ and 10³, the constant values are obtained faster for higher ratios of thermal conductivity coefficients although plotted characteristics do not converge to a single value. The variation of resulting quantities is characterized by slower dynamics of changes than in the case of smaller k₁/k₂ ratio.

Time evolution of LBM calculated k_eff(t) normalised to k_eff, series layout; the results from LBM calculations are for k₁/k₂ = 10³ and 10⁴ on left and right plot, respectively. The “max” and “min” lines refer to averaged values in the stream wise direction in the whole domain. LBM results are normed to the value of k_eff—saved after residuals get the specified value

For smaller k₁/k₂ ratios, the resulting error (not shown) is acceptable. For the parallel layout of layers, the error grows until some limit and seems to be constant with the growing ratio of thermal conductivities; with increasing resolution of the lattice, the error decreases. In case of series layout of layers the error decreases; for large ratios of thermal conductivity the error is difficult to determine as the computations take a very long time and do not reach a steady state, see Fig. 2. For comparison, the plot presents how the heat flux changes in function of simulation time until a steady state is reached for k₁/k₂ = 10³, whereas for k₁/k₂ = 10⁴ after a very long simulation time steady state in not reached. In case of the effective thermal conductivities for smaller k₁/k₂, the convergence of the results to the averaged value is clearly visible. This cannot be said about k₁/k₂ = 10⁴, where quartiles hold the constant value not close to the averaged value. Also (not presented on the plot) the maximal and minimal values of k_eff seem to be constant. This brings the conclusion that the Authors of the references (what is not specified) might have used a tuned relaxation parameter for heat transfer in their calculations at a given mesh resolution. Still the analysis shows applicability of LBM in most practical cases. In the example presented, one can assume that thermal conductivities for the coal (biomass) are of the order k ≈ 0.2–3 W/(m K) and for gas \(k = 10^{ - 2} {\text{W}}/({\text{m K}})\); this makes the maximum ratio of approximately 5 × 10³.

In a way to obtain much higher accuracy for large ratios of thermal conductivities k₁/k₂, \(\tau_{k}^{ - 1} = 0.5 + \delta \tau\) with δτ → 0 have to be used. Whatever the case, the stability has to be additionally supported by the increase of the spatial resolution. This will additionally increase the time of achieving the steady state for modelled heat transfer phenomenon. Figure 3 presents resulting error in function of input τ_ν which has influence on τ_λ (after some algebra τ_λ = 0.5 + 3(τ_ν − 0.5)α/ν, where α = λ/(ρc_p) for solid or fluid) parameters for IEDDF in case of a solid and fluid, see Eqs. (8)–(10). It is clearly visible that τ_ν parameter has a non-negligible effect on the accuracy of resulting k_eff as a consequence of (for fluid and solid) dependency on other simulation parameters (or directly on τ_ν). It could be argued that SRT LBM has stability issues. Nevertheless (see Wang et al. 2007a, b) it is not the case even for very high thermal conductivities ratios. Moreover, recent works are still reported where multiphysical phenomena in complex geometry (cf. Zachariah et al. 2019) are successfully described with SRT LBM approach. Research described by Zachariah et al. (2019) is based on a rather sophisticated model supplemented with additional relations to simulate drying in porous media. Similarly in present work, although the ratios of between thermal conductivities between solid and fluid are relatively high, criterial numbers describing the process are small (for details see beginning of Sect. 4).

Error estimation for the two layer system in a function of thermal conductivity ratio. On the X axis and on part of Y axis the logarithmic scale is used. The boundary conditions are imposed according to Fig. 1b, which corresponds to the series layout of layers

In the case of often used relaxation time value of τν = 0.75 calculations of heat transfer give results with unacceptable error even for k₁/k₂ = 50. For τν < 0.57 the error level is kept below 5% for k₁/k₂ = 10² whereas for τν = 0.501 the error is equal to 8.81% also for k₁/k₂ = 10⁴ (for k₁/k₂ = 10³ the error is at the level of 1%). This level of error is larger than the one reported by Wang et al. (2007a, b); nevertheless in discussion presented above it is clearly stated that further decrease of the relaxation parameter would noticeably improve the accuracy of estimated k_eff for higher ratios of k₁/k₂.

The study in the granular media has to account for fluid flow so Eqs. (1)–(3) will be solved by an in-house LBM solver. Free convection and gaseous products emission to the bulk at the surface of coal grains are two sources of fluid flow in considered geometry. Heat transfer is characterized by specified non-dimensional parameters, here for example the Prandtl, Grashof and Rayleigh numbers. The Reynolds number (that is Re = Ud/ν = 5.7 × 10 ⁻¹) is based on diameter of grains (d), maximal velocity (U) and viscosity of the fluid (ν); the Prandtl number for the fluid is Pr = ν/α = 0.71. The Grashof number Gr = gβ(θ_h − θ₀)/ν² = 0.82 is based on the acceleration due to gravity, thermal expansion coefficient and the difference between maximal and minimal temperatures at the domain boundary, respectively. As Gr characterizes the flow, the resulting value in the paper shows that the flow, mainly caused by free convection, is fully laminar. Very low value of Ra = GrPr shows that heat transfer is mainly caused by conduction. Although the low values of Re and Gr indicate small contribution of convection in heat transfer for the case of non-reactive granular media, in this work the directional release of gaseous products is implemented, as detailed in Sect. 4.2. At the grain surface, this may impair the heat transfer if the chemical reactions are introduced.

4.1 Non-reactive granular media case

The main purpose of the paper is to calculate effective heat conductivity for geometry of packed grains being a model of coal accounting for phenomena occurring during the coking process (for fluid flow see Grucelski and Pozorski 2013). The numerical domain (similarly as in Grucelski and Pozorski 2013,2015,2017) is created in a simple manner as described in Sect 2.5. The boundary conditions are set as follows, see Fig. 1b:

In a physically sound simulation in the meso-scale (grain scale), one should utilize fully periodic conditions with shifting of the pressure and temperature (on different pairs of boundaries) in order to represent the heating and cooling, see Grucelski and Pozorski (2016). In the best opinion of the authors, for the case described in their paper the appropriate boundary conditions for REV should be:

A similar condition can be implemented in spanwise direction, as the gaseous products release along with free convection introduce a pressure difference. Nevertheless, this work is not meant to address that problem in the exact manner (where the pressure difference depends on the height), thus only periodic condition is implemented.

An estimated effective thermal conductivity from LBM calculations for granular media is presented in Fig. 4 and compared to reference models known in the literature (see Wang et al. 2007a, b). Good accuracy is obtained within LBM calculations for heat transfer in computer generated granular media. When compared with LBM resulting effective thermal conductivity with reference values (from relationships), the calculations are close to the Maxwell–Eucken model and the Hashin–Shtrikman upper bound (referred as HS +), see Wang and Pan (2008). A good agreement is also obtained for the correlation used for the parallel layout of layers, see Wang and Pan (2008). An over prediction is observed in the case of the Self-Consistent model (SC) and for the series layout of layers (not presented in Fig. 4). Visibly higher deviations from the reference values are easy to observe close to domain boundaries where higher velocity of fluid flow are noticeable. As the model introduces buoyancy force (caused by density gradient in non-isothermal flow), the flow past solid grains causes deviation in calculated heat flux. Very good agreement can be observed in regions where fluid flow velocities are negligible (that is in the centre of the computational domain).

The heat flux (top plot) and local velocity (bottom plot) in the granular medium. The heat flux is calculated locally with temperature difference taken from the whole domain in steady state; chemical reactions are not considered

4.2 Reactive granular media case

With changing temperature the buoyancy force causes the cold fluid to penetrate the deposit of coal grains, which impacts heat transfer by scattering the flow on grains. Additionally, the gas release from the grain volume has a non-trivial influence on the heat transfer. Although the release rate is calculated with the use of solid averaged temperature, the flux of emitted gases at the surface corresponds to temperature gradient to imitate local conditions inside of the solid. In this model of gaseous products release from the solid grain, their flux will be non-uniform on the surface.

The model is implemented as follows. Basing on the layout of surface nodes for a given grain, the mass of gaseous products is normalized to unity over the surface. The used relationship gives a maximum value at a point closest to the heated boundary. Although the model reflects the physics, its influence on the results is only minor in the parameter range considered. Nevertheless, the model of non-uniform gaseous product release is accounted for in the computations.

A non-stationary temperature field in an exemplary geometry is presented in Fig. 5 (left panels). Close to the front of high temperature, sink of heat (negative heat source) is visible in the whole volume of grain. At the beginning of the simulation it is difficult to determine the change of heat flux caused by the occurrence of chemical reactions, due to high temperature variability. Possibly, the temperature condition at the left boundary of numerical domain had been chosen contrary to physical intuition: possibly, an increase of temperature at the heated boundary should be introduced as a function of time. At the moment (to the best of our knowledge) there is no relationship for boundary condition for REV in such a simulation. Figure 5 (right panels) presents curves for temperature and heat flux averaged over the domain for a later instant t = t₁. There are two points where the decrease of heat flux is clearly visible. A decrease of temperature in the solid volume caused by the secondary reactions (occurring at a higher temperature) take place in a much longer time range and the decrease of temperature is not so clearly visible on the colour map. Similarly no clear relationship is visible between pressure or velocity profiles (top right plot) and the magnitude of heat flux. Although the evolution of velocity profile corresponds with the temperature field in the domain, it is mainly connected with buoyancy force acting on the fluid. At the beginning of the simulation, due to a high temperature gradient, the free convection occurs close to the heated wall for which endothermic reactions take place. Further evolution makes both phenomena to proceed separately. Free convection has less significant effect on the effective thermal conductivity than chemical reactions and thermal dilatation of grains.

Temperature map (in greyscale) with the surface of grains and temperature isolines (left plot). The top and bottom plots correspond to t = t₀ and t = t₁, respectively. Right plot: fluid velocity and heat flux (top/bottom respectively) averaged in the spanwise direction; vertical dashed lines represent x-coordinates where the release of chemical compounds occurs (purple/green lines at t = t₀ and t = t₁, respectively)

Figure 6 presents the change of heat flux along with the release rate (both are averaged) in y-direction. The change of heat flux is mostly caused by chemical reactions causing the sink of heat. Other phenomena have much smaller effect on heat transfer. A resulting effective heat transport coefficient is presented in Fig. 7 together with the corresponding plot for non-reactive medium and heat transfer coefficients for solid and fluid. Local extrema correspond to local velocity peaks which also influence the heat flux. The increase of heat transfer coefficient value (blue line, Fig. 7) is caused by free convection. The lower value in case of simulation of reactive medium (purple line, Fig. 7) corresponds to the case where heat transfer is inhibited. At the moment, the model does not account for thermal expansion of grains (neither shrinking, see an earlier attempt Grucelski and Pozorski 2012), which can influence the heat transport in granular medium.

Release rate presented along with heat flux in the domain. The maximum of released products covers the spontaneous change of heat flux due to heat consumption required by chemical reactions

Comparison of effective heat transfer coefficient for granular media (reactive and non-reactive) together with reference value for fluid and solid

The chemical reaction (devolatilization) is treated as the endothermic process which reduces the temperature (uniformly in the grain volume) through a negative source term in the evolution equation for the IEDDF, see Eqs. (3)–(4). Taking into account the discussion regarding Fig. 5 (right bottom plot), the lower and higher values of the heat flux are considered in unsteady state. Thus, basing on the LBM results, the point x = x₀ is found where the averaged temperature becomes \(\langle \theta \rangle_{x_{0}}\) ≃590 K which is the threshold temperature for chemical reactions to occur (which is slightly less than recognised threshold temperature for coal decomposition). During the simulation, the heat flux is averaged over the domain in the following way. For high q on the left side of point x = x₀ in Fig. 5 range x ≤ x₀ it is

For low value of q on the right side of x = x₀ in Fig. 5 where chemical reactions do not occur, it is

Figure 8 presents the evolution of x₀ on the time axis along with a best-fit linear function. From the presented LBM results an effect of endothermic chemical reactions is visible during the process; clearly, the chemical reactions significantly affect the heating of the medium. Although the chosen REV is large enough for a reliable quantitative description of the process, for the smoother evolution curve a larger domain (in the spanwise direction) should be considered. In the chosen domain, after the heating of first obstacles, the thermal front propagates with a speed close to 0.265 (lattice nodes per time step, see the best fit line in Fig. 8). The initial discrepancy is caused by the inlet boundary condition: constant temperature seems to be unsuitable for meso-scale simulation in REV which causes faster propagation of thermal front. The right plot in Fig. 8 presents the evolution of q_max and q_min; the local extrema correspond to specific arrangement of grains in a given geometrical layout.

LBM simulation of gas release in heated granular medium: evolution of the thermal front (left plot) and the effective heat transfer coefficient (right plot) in function of time

The propagation of thermal front is mostly controlled by occurring heat sinks; after the temperature of a given grain has become high enough to trigger the gas release, the negative source (sink) becomes active in the internal energy equation.

This hinders the heating of given grain but not of the surrounding fluid; when chemical reactions will be terminated, the process of heating of the grain becomes faster again.

Figure 8 present estimation of effective thermal conductivity coefficient for q_max and q_min. The chosen REV allows one to calculate k_eff (the spatial averaging is used); here results are presented with the use of smooth curve (raw LBM results are presented with dots). The temperature field at a given x₀ is affected mostly by free convection and gas release. The latter greatly inhibits heat transfer, caused by sinks of heat flux In the considered domain, after the initial phase of calculations, the estimated k_eff → 0.48 whereas q_min is approximately 0.035. An important issue to notice is possible comparison of LBM results with experimental data from thermal processing of fossil fuel. The model in its current form does not account for another physically-relevant phenomenon: the thermal dilatation of grains along with mechanical interactions between grains.

4.3 Boundary scheme considerations for REV

Two different sets of boundary conditions were tested during the numerical tests. In the first approach, an element of volume is considered the element of volume close to the solid heated wall with the symmetry plane on the right boundary of the domain. The second case presented in this subsection assumes absence of the wall, whereas heating is done by infinitesimal wire placed at the left boundary of the domain and no changes are introduced at the right boundary.

Figure 9 presents the pressure and velocity for both cases. The resulting temperature difference is negligible. In the case where solid wall is not present in the domain, in the beginning of the simulation the main flow occurs into both directions.

Profiles of pressure and velocity averaged in stream-wise direction in function of the main-stream coordinate. The vertical lines show the area where devolatilisation occurs. The results are gathered in a moment where release of gaseous components takes place in the same region for the considered cases. The symbol q is a heat flux at west (at the left hand) wall; W, E determine condition at the domain boundaries: solid wall at west and zero velocity at the east boundary respectively

In the next example, at the left side in the domain (because of the solid wall), the fluid flow is blocked, thereby causing growth of the pressure. Due to modelled devolatilisation process, the rapid changes in the thermodynamic state occur in the region indicated by the grey lines in Fig. 9. Close to grains actually releasing the gaseous products one can observe change of velocity (averaged in vertical directions) and occurrence of points of interest on a line of pressure (like maximum value for green and purple lines in Fig. 9). So the main conclusion is that change of boundary conditions in the REV has noticeable influence on fluid flow in the domain. The future work direction would consist in considering various conditions on the REV boundary.

In the work, results from LBM modelling along with detailed investigation of heat transfer in granular medium are presented. The model of the Author has been further developed for the benchmark to estimate the effective heat transfer coefficient in such medium. Although very good agreement with references is obtained, the model has to be refined (in the meaning of chosen simulation coefficients) according to the cases considered. This requires adjusting the relaxation time for fluid (τ_ν) in LBM. In the benchmark case of two material layers, the error of LBM calculations of the effective thermal conductivity coefficient is acceptably small (up to a few per cent), except only for very high ratios of thermal conductivities of the layers in the series layout. Benchmark of LBM calculations for heat transfer in granular media also gives reasonable accuracy when compared with empirical correlations.

For the case of heat transfer phenomena in granular media with implemented model of devolatilisation, the calculations show that free convection along with conduction are greatly influenced by the gas release from grains. The locally occurring reactions cease the heating process by introducing a sink for the heat flux. After devolatilisation is accomplished for a given grain, the increase of temperature of the grain becomes faster. Due to model of chemistry used, the release of gaseous products is solved in two steps. The first step reactions clearly influence the heat transfer, the second stage is latent. The reason is the boundary scheme of the heated walls of REV; the high heat flux at the beginning of the simulation triggers in short time the first and second reactions. Afterwards, the release of chemical species does not create any clearly visible jumps in the value of heat flux. Boundary scheme was chosen to fit the volume to be part of the domain at the heated wall. This condition is not suitable for meso-scale simulation in REV. In case of outlet boundary, the shift-periodic boundary scheme (Grucelski and Pozorski 2016) was chosen to properly set the value of unknown DF. For this purpose as an input data (distribution functions), the model uses DF from the half length of the computational domain. The coefficients were set to keep the temperature at the outlet close to the initial temperature. This brings the conclusion that the condition proposed by Grucelski and Pozorski (2016) is suitable to provide scheme for meso-scale computations in REV with periodic conditions even at the heated boundary.

Described investigation, to the best knowledge of the Author, is a first approach to meso-scale modelling of effective thermal conductivity in reactive granular media of REV geometry. At the moment, thermal dilatation of grains is being implemented to check its influence in heat transfer for this modelling approach. In the next step, quantitative comparison of results from LBM modelling with known from literature results of simulations will be made.