Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 50, 1, pp. 193-214, Warsaw 2012 50th anniversary of JTAM LATTICE BOLTZMANN SIMULATION OF FLUID FLOW IN POROUS MEDIA OF TEMPERATURE-AFFECTED GEOMETRY1 Arkadiusz Grucelski Institute for Chemical Processing of Coal, Zabrze, Poland and Institute of Fluid-Flow Machinery, Polish Academy of Sciences, Gdańsk, Poland e-mail: agrucelski@imp.gda.pl Jacek Pozorski Institute of Fluid-Flow Machinery, Polish Academy of Sciences, Gdańsk, Poland e-mail: jp@imp.gda.pl The Lattice Boltzmann method (LBM) has been applied for flow and heat transfer computations. The simulations have been performedwith the single- relaxation time model and an advanced formulation of boundary conditions for LBM. For non-isothermal cases, a second distribution function has been used. First, validation tests are reported for heated flowpast a single obstacle as well as over a set of regularly and randomly arranged obstacles (grains) thatmake up a simplifiedmodel of a porousmedium.TheNusselt number for heat transfer in flowpast a single obstacle has been computed. Next, novel si- mulations of non-isothermal flow in a porousmedium of temperature-affected geometry have been undertaken. For the purpose, the thermal dilatation of grainshas been accounted for.Results are presented for the pressurehead loss and time-varying temperature profiles in the medium. Qualitative computa- tions accomplished to date constitute an encouraging first step to proceed further towards the impact of temperature-affected geometry in such flows, in particular for the coking process. Key words: Lattice BoltzmannMethod, porousmedia flow and heat transfer, variable geometry 1. Introduction The coking plants are widely used to provide chemically cleaner coal (coke), coke-oven gas and other products. A detailed analysis of the coking process 1 Presented at theXIXPolish FluidMechanics Conference, Poznań, September 5-9, 2010. 194 A. Grucelski, J. Pozorski has gained a renewed interest nowadays (cf. Nomura and Arima, 2000; Guo andTang, 2005), since ecological trends impose regulations to build industrial objects withBest Available Technology. Apart from fluid flow and heat trans- fer, changing shape and volume of coal grains is observed during the process; at the same time, production of the coking gas takes place as a consequence of chemical reactions occurring with growing temperature. The complexity of physico-chemical phenomena and complication of geometry imply that the use ofmore traditional tools and software of computational fluid dynamics (CFD) becomes prohibitively expensive. Hence the idea of a multiscale approach for detailed simulation purposes. In that approach, a microscopic (single-pore le- vel) computation of the representative element of volume (REV) is followed by amacroscopic CFD (system-level) computation (unsteady 1D/2D). Numerical simulations of the coking process in themacroscale need, as an input, several closure relationships that will provide variables describing geo- metrical properties of the domain that change during simulation.Only limited experimental data exist, along with a few examples of CFD investigations, cf. Guo and Tang (2005). We find it useful to provide the necessary input from detailed simulations in the microscale. As a numerical approach at the single-pore level, we have chosen the lattice Boltzmann method (LBM). TheLBM, based on theBoltzmann equation (He andLuo, 1997), has been developed in early 1980s as a method designed originally to be an extension of cellular automaton (Rothman and Zaleski, 1997) for eliminating large nu- merical noise. The LBM has proven to be suitable for simulations of viscous and nearly incompressible fluid flows (Succi, 2001; Chen andDoolen, 1998) as well as heat transfer (He et al., 1998) in simple (Yu and Girimaji, 2005) and complex geometries (Pan et al., 2006). The present authors have successfully applied the LBM to flowpast obstacles and also to thermal problems (Grucel- ski andPozorski, 2009, 2011). Thephenomenon of changing volume and shape of grains in functionof temperaturecanbedealtwithmoreadvancedboundary schemes for LBmethod. In our approach, the internal energy density distribu- tion function (IEDDF) solves for a temperature field, separately fromvelocity and pressure fields (Wang et al., 2007); the on-site interpolation-free (OSIF) boundary scheme accounts for the change of grains shape with temperature (Kao and Yang, 1998). LBM simulations in variable geometry are still under development. Some disadvantage of such calculation is a higher numerical cost caused by appli- cation of more complicated boundary schemes. An interesting case of flow simulation withmoving boundaries is reported byKrafczyk et al. (2001) whe- re authors present LBM results of fluid flow through a heart valve. Yet, we Lattice Boltzmann simulation of fluid flow in porous media... 195 are unaware of any publication where the heat transfer problem has been implemented for flowmodelling in variable geometry. In thepresentpaper,webuildonourexperience todateand furtherdevelop theLBMapproach.Webriefly recall themain idea of themethod anddescribe our current implementation of LBM to compute heat transfer coefficients in a flowpast a single obstacle. Then,we address the issue of flowandheat transfer in a generic (computer-generated) porousmediumandpresent some results for the pressure loss and temperature fields. Themain novel point of the present work is the accounting for variable geometry effects related to displacements of solid-fluid interfaces due to temperature increase in the porousmedium. 2. Lattice Boltzmann Method 2.1. Governing equations The numerical tool used in our simulations stands between molecular dy- namics (where we track position of every molecule and its velocity, like in simulations of fracture propagation) and computational fluid dynamics (whe- re we concentrate on macroscopic conservation laws, i.e., the Navier-Stokes equations). LBM is based on the Boltzmann equation (Succi, 2001) governing the evolution of density distribution function f that describes the probability density of finding particles with velocity u, at some infinitesimal volumewith the centre at r [∂t+u ·∇r +F ·∇u]f(r,u, t)= ∫ σ(∆u12,Ψ)(f1′2′ −f12) dΨ (2.1) where σ = σ(∆u12,Ψ) in the RHS collision term describes the number of particles with relative speed ∆u12 in a solid angle Ψ; symbols f12 and f1′2′ stand for two-particle distribution functions before and after collision, respec- tively. The LHS of Eq. (2.1) describes the Newtonian mechanics of a set of molecules; the RHS stands for the interaction term describing the number of particles lost from and coming into an element of the phase space as a result of collisions. The crucial feature of LBM (Rothman and Zaleski, 1997; Succi, 2001) is discretisation of the microscopic velocity vectors, both in direction and magnitude, cf. Fig.1 for some schemes used in 2D or 3D flows (9 or 15 velocities, respectively). A lot of algebra is needed to discretise the above equ- ation (He and Luo, 1997), also we have to use the H-theorem (Succi, 2001) to simplify its RHS. 196 A. Grucelski, J. Pozorski Fig. 1. Examples of LBM discretisation of velocity space in 2D and 3D Subsequently, the RHS of Eq. (2.1) is modelled by a relaxation term (cf. below). The equation is discretised in time (by the time step δt), in space (by the lattice size δx), and in velocity space by distinction of admissible directions ei (i subscript, cf. Fig.1) of the particles velocity on a regular square lattice. After these steps, the Lattice Boltzmann equation (suitably non-dimensionalised with δt as the time scale and δx as the length scale) takes the form (cf. Chen andDoolen, 1998) fi(r+eiδt,t+ δt)= fi(r, t)− τ −1 ν (fi−f eq i ) (2.2) where f eq i describes the equilibrium state of density fi(r, t) that for fluid flow problems has the following form f eq i = ρΩi [ 1+3eiv+ 9 2 (eiv) 2− 3 2 v 2 ] (2.3) In Eq. (2.3), ρ and v are the local fluid density and velocity, respectively, the weight coefficients Ωi depend only on the discretisation model of the veloci- ty space. For two presented sets of lattice velocity discretisation (D2Q9 and D3Q15), the coefficients of velocity expansion are constant. In Eq. (2.2), the relaxation time τν has the following non-dimensional form τν = 1 2 +3ν+ (2.4) In the above equation, ν+ = νδt/δx2 is the non-dimensional kinematic visco- sity of the fluid. Historically, two closely related schemes for LBMhave been proposed.The first is the single relaxation time(SRT)approachusedhere,Eq. (2.2), knownas Bhatnagar-Gross-Krook (BGK,Bhatnagar et al., 1954), where themicroscopic relaxation time is related to themacroscopic variable describingfluidviscosity, Eq. (2.4). Another approach developed by d’Humières et al. (2002) uses the Lattice Boltzmann simulation of fluid flow in porous media... 197 matrix of relaxation times; some additional cost and complexity are balanced by better stability properties and precision. The single relaxation time (SRT) approximation in the presented form is not appropriate for simulating heat transfer phenomena (He et al., 1998), inte- resting from the practical point of view. Themain reason is constant Prandtl number Pr= 1 and occurrence of instabilities in the course of computations. A possibility to overcome these difficulties while dealing with heat transfer problems in the BGK approximations is to extend the formulation with two distribution functions (He et al., 1998), namely the mass density fi and an additional, internal energy density distribution function (IEDDF) gi with its own relaxation time. Thus the BGK LBM with an additional equation for solving heat transfer becomes a double relaxation time approach. The scheme appears to be very stablewhen compared to other work for solving heat trans- fer problems with LBM (Peng et al., 2003b). Yet, the method reveals to be complex, mainly because of gradient terms in the evolution equation (Wang et al., 2007). In thepresentwork,weapplya simplifiedvariant of the IEDDF, elaborated byPenget al. (2003b). It iswidelyusedanddevelopedalso to implement source terms due, e.g., to chemical reactions (Wang et al., 2007). The thermal lattice Boltzmann equation can be written as gi(r+eiδt,t+ δt)= gi(r, t)− τ −1 θ (gi−g eq i ) (2.5) where g eq i stands for the internal energy density equilibriumdistribution func- tion in (r, t) and has the following form g eq i = θΩi[ai+ bieiv+ ci(eiv) 2+div 2] (2.6) where Ωi is the weight coefficient, θ is the local temperature (cf. below); the coefficients ai through di dependon thechosenmodel of velocity discretisation andalso on the direction i, unlike those inEq. (2.3).A full description of given equations, with exact arrays of coefficients forEq. (2.6), can be found inWang et al. (2007) for D2Q9 model and Peng et al. (2003a) for D3Q15 and D3Q19 models. The thermal relaxation time in phase m has the following form τθ = 1 2 + 3 2 λm ρcpm δt δx2 = 1 2 + 3 2 α+m (2.7) written both for the fluid f (gases) and solid s (coal grains) in the computa- tional domain, so m ∈ {s,f}. The heat diffusivity αm is expressed in terms of the heat conductivity λm and the specific heat cpm. The relaxation times 198 A. Grucelski, J. Pozorski τν and τθ depend on flow parameters (viscosity and thermal diffusivity). We note that both must be chosen within the range where the solution remains stable, for example τν ∈ (0.5,2). On the basis of the mesoscopic formulation, we use the following summa- tion expressions (discretised integrals) to obtain local macroscopic flow varia- bles (Wang et al., 2007), i.e., the fluid density, velocity and temperature ρ= ∑ i fi v= 1 ρ ∑ i fiei θ= ∑ i gi (2.8) Equations (2.2) and (2.5) with appropriate equilibrium distributions are presented in a non-dimensional form and are valid with the assumption of a low Mach number (M = vx/c ≪ 1) where vx is a flow velocity scale and c= δx/δt is the lattice sound speed.This assumption is controlled in twoways in the performed computation. First, themaximumflowvelocity is checked to be a fraction (less than 10%, say) of the sound speed. Second, the maximum variation of the density field is also checked to be smaller than several per- cent. Our experience has shown that larger variations in density may lead to numerical instabilities. In that case, computation is repeated on a finer mesh with a shorter time step. The LBM simulation is divided into three following steps for every lattice node at each time instant. First is a simple propagation of the distribution functions in discrete directions, cf. the LHS of Eqs. (2.2) and (2.5). In the second step we apply boundary conditions and collision (relaxation) terms for density and internal energy distribution functions. In the last step we perform gathering and calculating of macroscopic data from Eqs. (2.8). 2.2. Domain construction In our first LBM simulations for a single obstacle we considered the case of heat transfer in a hot fluid flow over a cold circular cylinder (cf. Sec. 3.2), originally used to calculate the Strouhal number of unsteady vortex shedding, cf. Grucelski and Pozorski (2009, 2011). As far as the construction of porous media is concerned, the domain con- sists of several REVs. A single REV is created by randomization of radii with uniform distribution within the prescribed range (rmin,rmax); centres of the obstacles are uniformlydistributed in the computational area.Thediameter of obstacles ismuch smaller than the domain size. In our simulation, the obstac- les can penetrate each other and also intersect the boundary of the domain; the algorithm next will periodically shift them (Fig.2). In the modelling of Lattice Boltzmann simulation of fluid flow in porous media... 199 flow and heat transfer in temperature-affected geometry, we start with a re- gular array of cylinders in a staggered arrangement (with centres located at vertices of equilateral triangles), followed by simulations in a random array of obstacles. Fig. 2. Construction of a simple porous medium for flow simulation with zoom of a single REV (gray scale represents the fluid velocity magnitude) In the case of flow in fixed geometry, the obstacles are represented as follows: every LBM node inside the circular cylinder with a given radius is given a “solid” identifier; outside nodes lying close (0.3δx, say) to the cylinder surface are also identified as solid. Next, every solid node which has at least one neighbouringfluid node is labelled as “interface”. For the case of changing shape, the nodes nearest to the real surface of the obstacle are next processed for more complicated boundary schemes. 2.3. Inlet/outlet boundary condtion The inlet condition is applied by simple calculation of fi in unknowndirec- tions from known distribution functions within the flow domain, the imposed inlet velocity and density. The following relationships are used for appropriate indices i∈{0,1,2,4,5,8}, cf. Fig.1a f neq i = f neq −i or fi = f−i− (f eq −i−f eq i ) (2.9) where f −i corresponds to the direction e−i =−ei; the non-equilibriumdistri- bution function is defined as f neq i = fi−f eq i . The imposed (macroscopic) fluid velocity and density at the flow inlet enter the distribution function f eq i thro- ugh Eq. (2.3). In the case of heat transfer computations, we additionally im- pose the inlet temparature. Again, it is used in the equilibriumdistribution of 200 A. Grucelski, J. Pozorski IEDDF, Eq. (2.6), at the flow inlet. In terms of the distribution function g neq i , we apply the following equation (cf.Wang et al., 2007) g neq i −e 2 ifi =−g neq −i +e 2 −if−i with g neq i = gi−g eq i . To implement the boundary condition at the outlet (where neither velo- city nor pressure is imposed), we use the extrapolation of known distribution functions for i∈{3,6,7} fi(r, t)= 1 2 [fi(r−e−iδt,t)+fi(r−2e−iδt,t)] (2.10) To fix the Reynolds number as a given parameter in the flow cases consi- dered, we set the inlet velocity vx, v = [vx,0]. Then, Re = d +v+x/ν + where d+ is the number of lattice nodes per typical obstacle diameter. 2.4. Wall boundary schemes in LBM For flow and heat transfer simulations we apply periodic conditions on the longer side of the domain (parallel to the main flow direction), cf. Fig.5. At the inlet we impose a uniform velocity profile on every fluid node with con- stant density and temperature. For the outlet nodes we use the extrapolation condition, based on nodes inside the domain. For fluid flow past fixed geometry we use the non-equilibrium (NEq) bounce-back scheme for boundary conditions at the solid-fluid interface. The scheme shows second-order accuracy and better stability than the original bounce-back scheme (Zou and He, 1997), acceptable implementation comple- xity and computational cost, with the following condition for unknown values of distribution function x+eiδt=S ⇒ f neq −i (x, t+2δt) = f neq i (x, t) (2.11) where S stands for a solid node. These boundary conditions are a developed formof the no-slip bounce-backmethod (Derksen, 2006) and showgood accor- dancewith results obtainedwithmore advanced boundary schemes (Grucelski and Pozorski, 2011). Basically, phenomena of flow in variable geometry cannot be simulated in LBM with the NEq method; for such computations, instabilities occur after changing the identifier of a single node (closer than 0.1δx from the interface, say) from fluid to solid one. To simulate phenomena of the shape change in function of growing temperature, we have to implement a more complex Lattice Boltzmann simulation of fluid flow in porous media... 201 boundary condition, allowing for reconstruction of the real interface. One of such methods is the on-site interpolation-free scheme (OSIF), described in details by Kao and Yang (1998). Aside of possibility of simulating moving boundaries, this method is free from interpolation that may cause a non-zero mass flux through the solid-fluid interface. TheOSIF scheme is based on weighting the relaxation time at nodes near the interface.As shown inFig.3,wedefine theweight coefficient q (0 20) are still possible in regular porous media where local magnitude of velocity is much smaller than the lattice speed of sound. With usage of the bounce-back scheme for single-relaxation-time (SRT) variant of LBM, it is known that the permeability shows some dependence on viscosity,more important than for themultiple relaxation time (MRT) variant (Pan et al., 2006). We want to point out that Pan et al. (2006) compare the simple bounce-back scheme with a few schemes for SRT LBM. More work is 204 A. Grucelski, J. Pozorski still needed inmatter of testing other boundary conditions for SRTLBM, also for the NEq and OSIF schemes used here for both fixed and temperature- dependent geometry. 3.2. Heat transfer To validate the LBM for non-isothermal flows, we considered the classical problem of cooled flow past a single obstacle, cf. Fig.5. The Nusselt number, Nu=Nu(Re,Pr), was determined for a couple of Re and Pr= ν/α. The local and global Nusselt numbers are obtained from Nu(φ)= d θw(φ)−θin ∂θ ∂n ∣ ∣ ∣ ∣ φ Nu= 1 2π ∫ φ Nu(φ) dφ where n is the unit vector normal to the solid-fluid interface, d is the cylinder diameter, θw is the wall temperature, θin is the flow inlet temperature, and φ is the polar angle measured w.r.t. to the inflow direction. The integral is calculated over the solid-fluid interface. Fig. 5. An example of flow past a circular cylinder at Re=150; (a) stream lines with velocity vectors and the temperature field coded with a scale of gray; (b) isolines of temperature (also with flow velocity vectors) superposed on the grayscalemap of velocity magnitude Theangular distribution of theNusselt number Nu(φ) for a single obstacle is presented in Fig.6. Because of coarse and square mesh (not a body-fitted one), the observed values of Nu show unphysical, local oscillations. They are Lattice Boltzmann simulation of fluid flow in porous media... 205 caused by the stair-shape character of the interface. In the next step, the collected values are averaged (typically, over δ = π/9) to obtain a smoother result according to Nu(φo)= 1 2δ φo+δ ∫ φo−δ Nu(φ) dφ (3.3) For such averaged valueswe calculate best-fitted, low-level polynomials, cf. Fig.6. There, the LBM results with the OSIF boundary scheme are obtained on the 500×250 lattice and comparedwith those reportedby Jie andHuiming (2006). In the approximating polynomial we omit the first-order term because of the symmetry constraints: ∂φNu(φ)|φ=0,π = 0 (i.e., at the upstream and downstream stagnation points). The resulting distributions (points fromLBM simulation) still show some variation, not present in the results of Jie and Huiming (2006). Those authors solved the flow with the LBM, and a more traditional CFD approach was used to solve the energy equation on the O- type body-fitted gridwhat greatly improved smoothness of their results. From the present simulation, the best-fitted curve is in a qualitative agreement with the one reported by Jie and Huiming. For a quantitative comparison, the experimental data of Acrivos et al. (1965) have been added in Fig.6. The present results seem to be quite satisfactory. Fig. 6. Local Nusselt numbers Nu in function of the polar angle φ for a few Re. Lines represent the best nonlinear fit to LBM simulation points: Nu(φ)= a0+a1cosφ+a2cos 2φ We want to point here that the OSIF scheme assumes weighting of the mass density distribution function only, while the IEDDF still sees a stair- shaped solid-fluid interface. An improved boundary scheme for non-isothermal 206 A. Grucelski, J. Pozorski flows remains an open issue. This brings us to conclude that for a simple rectangular lattice used in the present paper, the curve fitting (smoothing) is necessary. With growing Re we observe that local Nu grow for angles close to the downstream separation point (φ= π) because of vortex shedding. To the best knowledge of the authors, the results for angular distribution of the Nusselt number, where heat transfer is solved by the LBM, are reported for the first time here. In LBM computation, we first simulate isothermal fluid flow at the inlet temperature θin past the obstacle (initially, of temperature θo).When thevor- tex shedding becomes regular, the global Nusselt number has been calculated (Grucelski and Pozorski, 2009). The obtained values are presented in Fig.7. As we can see, the results show good conformity with empirical laws, like for example Daloglu and Unal (2000). Here, for approximation (lines) we use the relationship for the Nusselt number basing on the Ranz-Marshall correlation Nu=2+0.572Re0.5Pr0.3, where we fit its coefficients. For flow past a circular cylinder another empirical correlation is Nu=0.3609Re0.5749, cf. Daloglu and Unal (2000). In the case of Re→ 0, forwhich LBMsimulations loose stability, we add the point Nu(Re=0,Pr)=2. Fig. 7. Global Nusselt number Nu as a function of Re for a few Pr numbers. Solid lines represent the best nonlinear fit to simulation points with Nu= a+ bRecPrd Then, we have performed LBM simulations of convective heat transfer in a porous medium of fixed geometry. Figure 8 illustrates the heating of the medium.The temperature evolution in two selected cross-sections is shown in Fig.9. Currently, work is underway to retrieve a macroscopic law describing heat transfer in the porous medium. Lattice Boltzmann simulation of fluid flow in porous media... 207 Fig. 8. Evolution of the temperature field in a computer-created porousmedium (the main flow direction is upwards).White circles represent the boundary of the obstacles; (a) tvx/d=0.45, (b) tvx/d=2.25, (c) tvx/d=5.25 Fig. 9. Averaged temperature profile over lattice nodes at different downstream stations of REV 4. Variable geometry As mentioned in the Introduction, an essential part of the coking process is the thermal expansion and plastic deformation of coal grains. Therefore, we develop LBM simulations for the case of temperature-affected geometry of the porous medium. Here we assume that there is a maximum radius of each obstacle, and we simplify the model by imposing no mechanical interactions of the fluid-solid and solid-solid type (growing grainswould overlap each other eventually). For simulation of variable geometry with the NEq scheme, one will en- counter stability problems that may occur even for a change of only a few fluid nodes into solid ones. With the OSIF scheme, simulations of flow in temperature-affected geometry remain stable. In the case of overlapping inter- 208 A. Grucelski, J. Pozorski faces of the obstacles, one will be faced with stability issues connected with the thin layer (of one node width) of the fluid. To cope with those problems, in the case where there are two obstacles with interfaces in distance of a single lattice step, for the fluid node we choose the obstacle whose interface is clo- ser; for that particular value of q, cf. Eq. (2.12), we weight every distribution function pointing at the obstacle interface. Fig. 10. LBM simulations of the volume change with growing temperature of obstacles in regular geometry; four time instants. Lines represent constant temperature and gray-scalemap codes the magnitude of velocity (the main flow direction is upwards); (a) t=500, (b) t=1500, (c) t=3000, (d) t=4500 First, we considered a regular arrangement of N obstacles (cylinders)with a finite thermal conductivity and a temperature-dependent size, according to dn(t)= do[1+β(θn(t)−θo)] (4.1) where θo is the initial temperature of the set of obstacles (simulated poro- us medium) and dn(t), θn(t) are the (time-varying) diameter and average temperature of obstacle n (n=1, . . . ,N), respectively. Figure 10 shows first simulations with changing size of the obstacles (temperature-affected geome- try).Analogously to results discussed in Secs. 2 and3,wehave computed some integral characteristics of a simulated, variable-geometry porousmedium. Fi- gures 11 and 12 show the time evolution of porosity, the total pressure loss in the system and its average temperature, respectively. The two simulation variants are compared: fixed geometry and changing size of the obstacle accor- ding toEq. (4.1). It is seen that, with the increasing temperature, the pressure Lattice Boltzmann simulation of fluid flow in porous media... 209 loss exhibits a rapid growth,mainly due to the first upstream row of obstacles (already heated) that expand. Consequently, the viscous resistance will con- siderably increase. During flow evolution, where the porosity asymptotically attains a constant value, we also perceive much smaller change of pressure loss than the one observed at the start of flow. Time records of the averaged temperature profile are also different from those observed for flow past a fixed geometry. Figure 12 seems to agree with intuition that for flow with a larger pressure loss, the heat transfer through the array of obstacles is much slower, due to a smaller mass flow rate. This also causes the fluid to cool down faster. Fig. 11. Time evolution of: (a) porosity and (b) pressure difference between the inlet and outlet for flow through computer-created, variable porousmedia compared with results for flow in fixed geometry Fig. 12. Time record of the averaged (over fluid and solid nodes) temperature profile at two different downstream stations of the porous medium. Results for flow in variable geometry are compared with results for fixed geometry Next achievement for the variable-geometry case is an implementation of LBM for a simulation with obstacles that can cross each others interface. For 210 A. Grucelski, J. Pozorski the time being, however, we do not consider a true mechanical interaction of individual obstacles. Solution of the thermal deformation (which represents a problem in itself) is left for a subsequent work, along with the implementa- tion of plastic interactions between obstacles and next (if reasonable) fluid- structure interaction. An example of such a simulation is presented in Fig.13. There, plots (a) show a part of the domain where cold fluid flow is present (the gray maps represent the velocitymagnitude and temperature). Plots (b) present the same flow region (partially with changed geometry) with a partially heated fluid. In plots (c), we observe stopping of flow due to overlapping of growing obstacles. We note that shortly before the flow becomes blocked due to neighbouring obstacles getting close to each other, the local fluid velocity in a narrow gap between them can become quite high, occasionally causing numerical stability problems. Fig. 13. LBM simulations of volume change with growing temperature and intersections of obstacles in random arrangement. Upper plots: color-coded flow velocity magnitude (upwardsmean flow); lower plots: temperature map; (a) t=2000, (b) t=10000, (c) t=18000 5. Discussion and conclusions In the present work, we have developed LBM into simulation of fluid flow and heat transfer inporousmedia to retrievemacroscopic variables describingboth phenomena at theREV level.We have shown a comparison of the results from Lattice Boltzmann simulation of fluid flow in porous media... 211 simulation of fluid flow and heat transfer with empirical laws. The pressure loss in function of the inlet velocity allows us to calculate the permeability fromLBMsimulation data.Thepermeability obtained in thatway shows good conformitywith theDarcy-Forchheimer lawwhere the second-order coefficient has been fitted to data points. As we show, the pressure loss on a single REV for different lengths of a domain may vary, but can still be described by the Darcy equation. Using the Forchheimer term (for the inertial effects occurring in faster flow past a porous medium) we can describe the change of pressure loss also for faster (less viscous) flows. During our work, also a change of the permeabilitywith Rewas observed.Awiderdiscussion canbe foundelsewhere (Pan et al., 2006; Grucelski and Pozorski, 2011). Heat transfer inflowover a single obstacle hasbeen checked alsoby compa- risonwith empirical laws. Local values of Nu(φ) show some oscillations due to a coarse interface, occurring in the heat transfer problem for both body-fitted (OSIF) and mesh-fitted (NEq) schemes. Yet, the best-fit polynomials show a good agreement with the experimental data; also, they have been compared to other results connected with LB method (Jie and Huiming, 2006). Global Nusselt numbers for such a flow can also be well correlated by the best-fit curves of the presumed type. Still more work is needed to obtain a better agreement of this type of results, as the local Nusselt numbers, with other numerical schemes (reference data are available). Because Eq. (2.2) does not have any source term, the distribution function does not depend on temperature. (Later we are going to add sources due to chemical reactions.) Consequently, the applied model is not yet ready to simulate the coking process. At the moment, we can compute heat and fluid flow phenomena in temperature-affected, complicated geometry (of regular or random placement of obstacles). The presented results (change of pressure losswith changingporosity) are consistentwith thephysical expectation about flow past a porous media; basing on the mentioned results, we also conclude that the changes of averaged temperature for fixed and variable geometry are qualitatively correct. Still, there is a need for quantitative results from LBM simulations of flow and heat transfer phenomena in variable geometry. Next, we plan to focus on two subsequent developments. First is an ap- propriate treatment of the situation where the topology changes due to, e.g., nearby solid grains coming in contact (solid-solid interaction). Secondly,we in- tend to implement the OSIF boundary scheme for the heat transfer problem. As the ultimate objective, we intend to develop an efficient and accurate me- thod to model the coking of coal. We will use the obtained coefficients as closures in a physically-sound, one-dimensional simulation of the process. 212 A. Grucelski, J. Pozorski Acknowledgements The investigations presented in this paper have been partially financed within the research project POIG.01.01.02-24-017/08 “Intelligent coking plant meeting the requirements of best available technology” („Inteligentna koksownia spełniająca wy- magania najlepszej dostępnej techniki”). References 1. Acrivos A., Snowden D.D., Grove A.S., Petersen E.E., 1965, The ste- ady separatedflowpast a circular cylinder at largeReynolds numbers, J. Fluid. Mech., 21, 737-760 2. Bhatnagar P.L., Gross E.P., Krook M., 1954, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one- component systems,Phys. Rev. A, 94, 511 3. Chen S., Doolen G.D., 1998, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech., 30, 329-364 4. Grucelski A., Pozorski J., 2009, Lattice Boltzmann simulation of flow and heat transfer in porousmedia,Conference onComputationalMethods for Ther- mal Problems, Naples, Italy 5. GrucelskiA.,Pozorski J., 2011,LatticeBoltzmann simulations of flowpast an obstacle and in simple porous media, submitted toComputers and Fluids 6. Daloglu A., Unal A., 2000, Heat transfer from a cylinder in the wake flow, Int. Comm. Heat Mass Transfer, 27, 569-580 7. Derksen J.J., 2006,The latticeBoltzmannmethod for (multiphase) fluid flow simulations, Lecture Notes, CISM School, Udine 8. d’Humières D., Ginzburg I., Krafczyk M., Lallemand P., Luo L., 2002, Multiple-relaxation-time lattice Boltzmann models in three dimensions, Phil. Trans. R. Soc. Lond. A, 360, 437-451 9. Guo Z., Tang H., 2005,Numerical simulation for a process analysis of a coke oven,China Particuology, 3, 373-378 10. He X., Chen S., DoolenG.D., 1998, A novel thermalmodel for the Lattice BoltzmannMethod in incompressible limit, J. Comput. Phys., 146, 282-300 11. He X., Luo L.S., 1997, Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation, Phys. Rev. E, 56, 6811-6817 Lattice Boltzmann simulation of fluid flow in porous media... 213 12. Jie W., Huiming Z., 2006, Laminar flow and heat transfer from a circular cy- linder in cross-flowusing theLattice-Boltzmannmethod,TermPaper,National University of Singapore 13. Kao P.H., Yang R.J., 1998, An investigation into curved and moving bo- undary treatments in the lattice Boltzmannmethod, Annu. Rev. Fluid Mech., 146, 282-300 14. Kim S.M., Ghiaasian S.M., 2009, Numerical modeling of laminar pulsating flow in porousmedia, J. Fluids Engng., 131, 141203-9 15. Krafczyk M., Tolke J., Rank E., Schulz M., 2001, Two-dimensional si- mulation of fluid-structure interaction using lattice-Boltzmannmethods,Com- put. Struct., 79, 2031-2037 16. Nomura S., Arima T., 2000, Coke shrinkage and coking pressure during carbonization in a coke oven,Fuel, 79, 1603-1610 17. Pan C., Luo L.S., Miller C.T., 2006, An evaluation of lattice Boltzmann schemes for porous media flow simulation,Comput. Fluids, 35, 898-909 18. Peng Y., Shu C., Chew Y.T., 2003, A 3D incompressible thermal lattice Boltzmannmodel and its application to simulate natural convection in a cubic cavity, J. Comput. Phys., 193, 260-274 19. Peng Y., Shu C., Chew Y.T., 2003, Simplified thermal lattice Boltzmann model for incompressible thermal flows,Phys. Rev. E, 68, 026701 20. RothmanD.H.,Zaleski S., 1997,Lattice-GasCellularAutomata,Cambridge University Press 21. Scheidegger A.E., 1957, The Physics of Flow Through Porous Media, Uni- versity of Toronto Press 22. Succi S., 2001,The Lattice Boltzmann Method for Fluid Dynamics and Bey- ond, Clarendon Press, Oxford 23. Vafai K., Amiri A., 1998, Non-Darcian effects in confined forced convective flows,Transport Phenomena in Porous Media, 1, 313-329 24. WangJ.,WangM., Li Z., 2007,A latticeBoltzmannalgorithm forfluid-solid conjugate heat transfer, Int. J. Therm. Sci., 46, 228-234 25. Yu H., Girimaji S.S., 2005, Near-field turbulent simulations of rectangular jets using lattice Boltzmannmethod,Phys. Fluids, 17, 125106 26. Zou Q., He X., 1997, On pressure and velocity boundary conditions for the lattice Boltzmann BGKmodel,Phys. Fluids, 9, 1591-1598 214 A. Grucelski, J. Pozorski Symulacje przepływu metodą siatkową Boltzmanna w ośrodku porowatym o geometrii zależnej od temperatury Streszczenie Metoda siatkowa Boltzmanna (LBM) została zastosowana do obliczeń przepły- wu i wymiany ciepła. Symulacje zostały przeprowadzone dla modelu z pojedynczym czasem relaksacji oraz zaawansowanego schematuwarunkówbrzegowychmetody LB. Dla przepływównieizotermicznych, użyto dwóch funkcji rozkładu. Przedstawiono ob- liczenia testowe nieizotermicznego opływu pojedynczej przeszkody, jak również opły- wu regularnie oraz losowo rozmieszczonych przeszkód (ziaren), które tworzą uprosz- czony model ośrodka porowatego.Wyznaczono liczbę Nusselta dla przepływu ciepła w opływie pojedynczej przeszkody.Podjęto symulacje połączonych zjawisk (przepływ i wymiana ciepła) dla ośrodka porowatego o zmiennej geometrii przeszkód. W tym przypadkuuwzględniono rozszerzalność termiczną ziaren. Zaprezentowanowyniki dla straty ciśnienia oraz zmienne w czasie profile temperatury dla przepływu przez ośro- dek. Uzyskane dotąd jakościowe wyniki stanowią pierwszy krok do dalszych badań wpływu zależnej od temperatury geometrii ziaren ośrodka na przebieg takich prze- pływów, a w szczególności na proces koksowania. Manuscript received December 1, 2010; accepted for print May 13, 2011