- Email: [email protected]
Numerical analysis of the pressure drop in porous media flow with lattice Boltzmann (BGK) automata
Computer Physics Communications 129 (2000) 247–255 www.elsevier.nl/locate/cpc
Numerical analysis of the pressure drop in porous media flow with lattice Boltzmann (BGK) automata J. Bernsdorf ∗ , G. Brenner, F. Durst Institute of Fluid Mechanics (LSTM), University of Erlangen–Nuremberg, Cauerstr. 4, D-91058 Erlangen, Germany
Abstract The lattice Boltzmann (LB) method is used for a detailed study on the origins of the pressure drop in porous media flow. In agreement with the experimental results [Durst et al., J. Non-Newtonian Fluid Mech. 22 (1987) 169] it is shown, that the elongation and the contraction of fluid elements is an important factor for the pressure loss in porous media flow, and that a significant error is made, when only shear forces are taken into account. To obtain the geometry information of the porous media for our simulations, we used the 3D computer tomography technique. 2000 Elsevier Science B.V. All rights reserved. Keywords: CFD; Porous media flow; Lattice Boltzmann; Computer tomography
1. Introduction It is a well known fact that the lattice Boltzmann method is a very efficient numerical tool to investigate flows in highly complex geometries, such as porous media [2,3]. It allows a detailed discretization of the porous geometry and thus the exact simulation of the transport of mass and momentum without any of the underlying semi-empirical homogenization models, generally used in engineering applications. Thus, on the one hand, they allow one to investigate the transport phenomena in porous media and to improve the basic understanding, e.g., of the high viscous pressure losses in these flows. On the other hand, valuable information entering in the formulation of homogenization models can be obtained from these data. In that respect, the LB simulation may be considered as a ‘numerical experiment’ with the additional advantage, that more information about local flow properties can be obtained than usually possible in experiments. The present homogenization approaches are based on the assumption of a linear or quadratic relationship between pressure gradients and the mean velocity (Darcy or Forchheimer law). Usually, a proportionality is assumed with a constant permeability. This parameter has to be estimated from experiments or from theoretical considerations, such as the Kozeny–Darcy equation (see, e.g., [4]). In this theory, it is assumed that the porous media can be modeled as a bundle of capillaric tubes. Only shear forces in a laminar Poiseuille like flow are taken into account. Any forces due to elongation and contraction of fluid elements are neglected. In [1] it has been shown, that this approach leads to an underestimation of the momentum loss by a factor of 2.5 and a significant discrepancy ∗ Corresponding author. E-mail: [email protected]. Present address: C&C Research Laboratories, NEC Europe Ltd., Rathausallee 10, D-53757 Sankt Augustin, Germany.
0010-4655/00/$ – see front matter 2000 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 0 0 ) 0 0 1 1 1 - 9
248
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
to experimental results. This discrepancy is usually explained by the introduction of a tortuosity factor considering the effect of the additional length due to the complex structure of the flow path. In the present paper, the results of detailed numerical simulations are used to investigate quantitatively the effect of elongational forces and their contribution to the pressure loss in porous media flows. In the next section, we summarize the basic idea of the capillaric theory in comparison to some experimental results. In the third section, we describe briefly the lattice Boltzmann method applied to this investigation including one validation test case. Then we present two simulations of a porous media flow. In the fourth part, the results of these two simulations are evaluated with regard to shear- and elongational forces and their contribution to the pressure drop.
2. Analytical analysis of the pressure drop In general, the idea behind analytically estimating a viscous porous media flow is to define a relation which describes the pressure drop as a function of the geometry (e.g., porosity, specific surface), fluid parameters (density, viscosity) and flow parameters (velocity), ∂x P = f (geometry, fluid, flow).
(1)
2.1. The Kozeny–Darcy equation The most common models which address this problem can be summarized under the term ’capillaric theories’. There, the porous media flow is being modeled as a flow through a bundle of channels with not strongly changing cross-sections. For each of these channels, the Navier–Stokes equation can be solved, and the relationship of the ex and the pressure drop for a single channel results in: mean flow velocity U ex = − 1 dP Rh2 , (2) U 2µ dx where P is the pressure, µ the fluid viscosity and x the mean flow coordinate. The hydraulic radius Rh , which is defined as the ratio of fluid volume and wetted surface, can be derived from the radius R of the pipe by: Rh =
R πR 2 = . 2 πR 2
(3)
eh and a length L, Eq. (2) can be rewritten Now assuming a bundle of tubes with an average hydraulic radius R as: e2 . ex = − 1 1P R (4) U 2µ 1L h Neglecting the underlying channel-geometry, one can try to use this formula as a general expression for calculating the pressure drop in other types of geometries, assuming that the average hydraulic radius is known eh can be ep . R or can be derived. For example, consider a porous media built up of spheres of an average diameter D e written as a function of Dp and its porosity ε: ep D ε . 6 (1 − ε) Inserting this expression into Eq. (4) yields: eh = R
e2 D ε2 ex = 1 U0 = 1 1P p , U ε 2µ 1L 36 (1 − ε)2 where U0 is the so called ‘effective velocity’ inside the porous media.
(5)
(6)
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
Eq. (6) can be rewritten in so called Ergun coordinates (see, e.g., [4]), ep ep ρ U0 D ε3 1P D . 72 = 1L ρU02 (1 − ε) µ(1 − ε)
249
(7)
Introducing the dimensionless quantities Reynolds number Re and friction factor f as defined by ep ε3 1P D , (8) 2 1L ρU0 (1 − ε) ep ρ U0 D , (9) Re = µ(1 − ε) the advantage of using the Ergun coordinates becomes obvious, because Eq. (7) can be written in the compact form Λth , (10) f= Re with the ‘theoretical’ friction coefficient Λth = 72. f=
2.2. Experimental results The friction coefficient was measured by Durst [1] experimentally for a packed bed of spheres with different diameters at a wide range of Reynolds numbers (see Fig. 1). For the Reynolds numbers below Re = 1, the experimentally determined friction coefficient appears to be a constant with Λexp = 182. This value is about 2.5 times higher than the theoretical one derived within the capillaric theory. To explain this additional pressure loss, it is usually argued, that the capillaric theories do not take into account the complex paths, the fluid normally has to go through a porous media. To take into account this longer flow paths compared to the edge length of a porous media, the so called tortuosity factor is introduced as follows (see, e.g., [4]): τ=
length of flow paths . macroscopic length scale
(11)
It might be doubted that a tortuosity factor of τ = 2–3 is a realistic assumption, because this would imply that the length of the fluid channels is up to three times larger than the length of the porous media. Also, there exists
Fig. 1. Friction coefficient versus Reynolds number (section from [1]).
250
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
no way to measure the tortuosity directly. Some a posteriori derivation based on comparing the measured pressure drop to the one derived theoretically seems to be a highly questionable – although often practiced – method. In the following section, we will show that there exists another, and until now, by the capillaric theories not recognized physical effect causing pressure drop in porous media flow.
3. Numerical simulation of a porous media flow For the analysis of the pressure drop of a flow through a packed bed of spheres, we used the lattice Boltzmann method described shortly below to compute the detailed flow field. This flow field was then analyzed with respect to the elongational stresses, which can produce additional pressure losses. 3.1. Numerical method For the computations, a 3D 19-speed (D3Q19) lattice Boltzmann automata with single time Bhatnager–Gross– proposed by Qian [5] is used: Krook (BGK) relaxation collision operator 4Boltz i , Ni (t∗ + 1, rE∗ + cEi ) = Ni (t∗ , rE∗ ) + 4Boltz i 4Boltz i
eq = ω(Ni
(12)
− Ni ),
(13) eq Ni :
with a local equilibrium distribution function ciα uα uα uβ ciα ciβ eq − δ . Ni = tp % 1 + 2 + αβ cs 2cs2 cs2
(14)
eq
This local equilibrium distribution function Ni is chosen in such a way to recover the time-dependent Navier– Stokes equation [5] in the low mach number limit: ∂t ρ + ∂α (ρuα ) = 0,
∂t (ρuα ) + ∂β (ρuα uβ ) = −∂α p + µ∂β ∂β uα + ∂α uβ .
(15) (16)
3.2. Boundary conditions A parabolic velocity inlet profile and a fixed pressure at the outlet were chosen for all test cases. This was achieved by introducing the equilibrium density distribution at the first and last lattice column computed with an upstream extrapolated pressure and a downstream extrapolated flow velocity for the inlet and outlet respectively. The inlet and outlet region was chosen to be long enough to prohibit any errors introduced by this method from affecting the measured quantities. A slightly improved bounce back boundary condition for the wall boundaries (including one relaxation step per iteration also on the solid domain nodes) ensures an accurate and relaxation parameter independent position of the wall for a range of at least 1.8 6 ω 6 1.95. 3.3. Validation The code used for this investigation was already validated in detail for time dependent flow around simple geometries [8]. To produce quantitatively reliable data when using CFD, it is also necessary to ensure the grid independence of the numerical results. This is usually done by discretizing the same geometry with meshes of increasing sizes and observing the convergence of the results with increasing mesh refinement. Especially a discretization of spherical objects with rectangular elements, as it is common practice within the LB technique, makes it necessary to carefully investigate the discretization error.
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
251
Fig. 2. Packed bed of spheres, surface shaded with the pressure.
Fig. 3. Friction coefficient versus particle diameter.
As a validation test-case, the pressure drop for low Reynolds number flow through an orthorhombic package of spheres was simulated. The domain sizes were chosen to be lx × ly × lz = 120 × 20 × 20 lattice nodes for the 8 × 2 × 2 spheres with a diameter of Dp = 10 for the coarsest resolution and lx × ly × lz = 480 × 80 × 80 lattice nodes with Dp = 40 for the finest resolution. Periodic boundary conditions were applied normal to the main flow direction to make this test case similar to the experimental set up of Durst et al. [1]. A very good convergence of the numerically achieved friction coefficient to the experimental values (see Fig. 1) can be observed in Fig. 3. A sphere diameter of Dp = 20 is sufficient to approach the convergence result within an accuracy of bet ter than 3%. 3.4. Flow simulation through porous media In order to have realistic geometries for the pressure drop studies, two porous media from engineering applications were chosen: One sponge-like SiC matrix, and one catalyst consisting of a tube filled with spheres.
252
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
Fig. 4. Computer tomography data as input for the LB simulation: catalyst.
3.4.1. Geometry preprocessing For both samples, the catalyst as well as the SiC matrix, the geometry was digitized using 3D computer tomography (3D CT). The tomography data can directly be used as geometry input for the lattice Boltzmann simulation, without the need of any further preprocessing [6,7]. 3.4.2. Case 1: Catalyst A cylindrical porous probe with a height of 110 mm and a diameter of 80 mm was scanned using 3D CT with an average resolution of 0.9 mm. This leads to a discretization of lx × ly × lz = 123 × 90 × 90 voxels. With these parameters, the average diameter of the spheres is in a range which allows one to predict correct results according to the convergence study in the foregoing section. The complex geometry data were centered inside an lx × ly × lz = 250 × 99 × 99 sized channel (see Fig. 4), and a flow for a Reynolds number of about Re ≈ 0.1 was simulated using velocity inlet and pressure outlet boundary conditions. The simulation was performed on one processor of a VPP 700 at the Leibniz–Rechenzentrum in Munich; 50,000 iterations were necessary for this set-up, which took about 25,200 CPU seconds, 850 Mbyte of computer memory were necessary for the storage of the ≈ 2.45 × 106 voxels. 3.4.3. Case 2: SiC matrix A cylindrical porous probe with a height of 30 mm and a diameter of 82 mm was scanned using 3D CT with an average resolution of 0.5 mm. This leads to a discretization of lx × ly × lz = 44 × 147 × 147 voxels. The average diameter of the flow channels is again large enough with this resolution to produce lattice-independent results. The complex geometry data were centered inside an lx × ly × lz = 100 × 149 × 149 sized channel, and a flow for a Reynolds number of about Re ≈ 0.1 was simulated using velocity inlet and pressure outlet boundary conditions. The simulation was performed on one processor of a VPP 700 at the Leibniz–Rechenzentrum in Munich; 10,000 iterations were necessary for this set-up, which took about 5760 CPU seconds, 800 Mbyte of computer memory were necessary for the storage of the ≈ 2.2 × 106 voxels.
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
253
Fig. 5. Pressure field for the flow through a catalyst (x–z plane).
Fig. 6. Computer tomography data as input for the LB simulation: SiC matrix (right: section).
4. Analysis of the pressure drop from experimental data As argued above, the tortuosity is obviously not the only reason for the higher pressure drop observed in experiments and numerical simulations when compared to the results derived from the capillaric theories. The total dissipation in the flow when passing through a porous media can be expressed by: Φ = −τij with
∂Uj , ∂xi
∂Uj ∂Ui + −τij = −µ ∂xi ∂xj
(17)
∂Uj 2 , − µδij 3 ∂x | {z i} =0
for incompressible fluid. Eq. (17) can be rewritten as ∂Uj ∂Ui ∂Uj + Φ = −µ ∂xi ∂xj ∂xi
(18)
254
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
Fig. 7. Pressure field for the flow through a porous SiC matrix (left: x-z plane, right: istotache, shaded with the pressure).
z = −µ
Φe (elongation)
∂u1 ∂x1
2
}| { ∂u2 2 ∂u3 2 ∂u3 ∂u2 2 ∂u1 ∂u3 2 ∂u2 ∂u1 2 + + + + + + + . + ∂x2 ∂x3 ∂x ∂x2 ∂x2 ∂x3 ∂x3 ∂x1 {z } | 1
Φs (shear)
(19) The dissipation can thus be expressed as a sum of two parts: the dissipation caused by shear forces Φs and the dissipation caused by elongational strain Φe : Φ = Φs + Φe .
(20)
An evaluation of the detailed flow fields produced by the numerical simulation with the fraction, Φe /Φs , yields: Example 1 (Catalyst): Friction coefficient
Λ = 217.6,
Dissipation (elongation) Dissipation (shear) elongation / shear
Φe = 2.111 ∗ 10−02 , Φs = 2.854 ∗ 10−02 , Φe /Φs = 0.74.
Example 2 (SiC matrix): Friction coefficient
Λ = 342.4,
Dissipation (elongation) Dissipation (shear) elongation / shear
Φe = 3.853 ∗ 10−04 , Φs = 6.486 ∗ 10−04 , Φe /Φs = 0.59.
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
255
In both examples, the friction coefficient Λ is much larger than predicted by the capillaric theory. A considerable amount of the pressure drop is caused by the elongational strain Φe , what can be clearly observed from the relation Φe /Φs . In the derivation of the Kozeny–Darcy law (as well as in similar approaches), only the shear forces are taken into account. The gap for the pressure drop between the theoretical values due to these shear forces and the experimental results is usually explained by the tortuosity of the porous media, while the contribution of the elongational forces is neglected. Although for the above examples tortuosity causes a fraction of the observed pressure drop, it could be shown that this is not the only reason for the difference between the experimental and the theoretical results. So, an a posteriori adjustment of the tortuosity with a much too large tortuosity factor of τ ≈ 2–3, while neglecting the contribution of the elongational forces, may lead to a considerable error.
5. Conclusion With the low Reynolds number flow simulation in regularly packed beds of spheres we could show that the lattice Boltzmann method is able to quantitatively predict pressure drop with high accuracy. Three-dimensional computer tomography was used to digitize the geometry of two porous media: a catalyst and a porous SiC matrix. It could be shown by the a posteriori analysis of the simulated flow fields, that elongational strain gives an important contribution to the pressure drop. This must not be neglected when deriving a theory for predicting the pressure drop in porous media flow. The derivation of a tortuosity factor from pressure drop measurements while neglecting the dissipation due to elongational strain might produce a considerable error.
Acknowledgements The authors would like to thank the Hattinger Prüf- und Entwicklungs-GmbH (HAPEG) for providing the computer tomography data. This project was founded by the Deutsche Forschungsgemeinschaft (Proj. Nr. Br 1864/1). References [1] [2] [3] [4] [5] [6] [7] [8]
F. Durst, R. Haas, W. Interthal, J. Non-Newtonian Fluid Mech. 22 (1987) 169. D.H. Rothman, Geophysics 53 (1988) 509. J. Bernsdorf, M. Schäfer, F. Durst, Int. J. Numer. Meth. Fluids 29 (1999) 251. R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena (New York, 1960). Y.H. Qian, Europhys. Lett. 17 (6) (1992) 479. B. Ferréol, D.H. Rothman, Transport in Porous Media 20 (1995) 3. J. Bernsdorf, O. Günnewig, W. Hamm, O. Münker, GIT Labor-Fachzeitschrift 4/99 (1999) 387. J. Bernsdorf, Th. Zeiser, G. Brenner, F. Durst, Int. J. Mod. Phys. C 9 (8) (1998) 1141.
Numerical analysis of the pressure drop in porous media flow with lattice Boltzmann (BGK) automata J. Bernsdorf ∗ , G. Brenner, F. Durst Institute of Fluid Mechanics (LSTM), University of Erlangen–Nuremberg, Cauerstr. 4, D-91058 Erlangen, Germany
Abstract The lattice Boltzmann (LB) method is used for a detailed study on the origins of the pressure drop in porous media flow. In agreement with the experimental results [Durst et al., J. Non-Newtonian Fluid Mech. 22 (1987) 169] it is shown, that the elongation and the contraction of fluid elements is an important factor for the pressure loss in porous media flow, and that a significant error is made, when only shear forces are taken into account. To obtain the geometry information of the porous media for our simulations, we used the 3D computer tomography technique. 2000 Elsevier Science B.V. All rights reserved. Keywords: CFD; Porous media flow; Lattice Boltzmann; Computer tomography
1. Introduction It is a well known fact that the lattice Boltzmann method is a very efficient numerical tool to investigate flows in highly complex geometries, such as porous media [2,3]. It allows a detailed discretization of the porous geometry and thus the exact simulation of the transport of mass and momentum without any of the underlying semi-empirical homogenization models, generally used in engineering applications. Thus, on the one hand, they allow one to investigate the transport phenomena in porous media and to improve the basic understanding, e.g., of the high viscous pressure losses in these flows. On the other hand, valuable information entering in the formulation of homogenization models can be obtained from these data. In that respect, the LB simulation may be considered as a ‘numerical experiment’ with the additional advantage, that more information about local flow properties can be obtained than usually possible in experiments. The present homogenization approaches are based on the assumption of a linear or quadratic relationship between pressure gradients and the mean velocity (Darcy or Forchheimer law). Usually, a proportionality is assumed with a constant permeability. This parameter has to be estimated from experiments or from theoretical considerations, such as the Kozeny–Darcy equation (see, e.g., [4]). In this theory, it is assumed that the porous media can be modeled as a bundle of capillaric tubes. Only shear forces in a laminar Poiseuille like flow are taken into account. Any forces due to elongation and contraction of fluid elements are neglected. In [1] it has been shown, that this approach leads to an underestimation of the momentum loss by a factor of 2.5 and a significant discrepancy ∗ Corresponding author. E-mail: [email protected]. Present address: C&C Research Laboratories, NEC Europe Ltd., Rathausallee 10, D-53757 Sankt Augustin, Germany.
0010-4655/00/$ – see front matter 2000 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 0 0 ) 0 0 1 1 1 - 9
248
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
to experimental results. This discrepancy is usually explained by the introduction of a tortuosity factor considering the effect of the additional length due to the complex structure of the flow path. In the present paper, the results of detailed numerical simulations are used to investigate quantitatively the effect of elongational forces and their contribution to the pressure loss in porous media flows. In the next section, we summarize the basic idea of the capillaric theory in comparison to some experimental results. In the third section, we describe briefly the lattice Boltzmann method applied to this investigation including one validation test case. Then we present two simulations of a porous media flow. In the fourth part, the results of these two simulations are evaluated with regard to shear- and elongational forces and their contribution to the pressure drop.
2. Analytical analysis of the pressure drop In general, the idea behind analytically estimating a viscous porous media flow is to define a relation which describes the pressure drop as a function of the geometry (e.g., porosity, specific surface), fluid parameters (density, viscosity) and flow parameters (velocity), ∂x P = f (geometry, fluid, flow).
(1)
2.1. The Kozeny–Darcy equation The most common models which address this problem can be summarized under the term ’capillaric theories’. There, the porous media flow is being modeled as a flow through a bundle of channels with not strongly changing cross-sections. For each of these channels, the Navier–Stokes equation can be solved, and the relationship of the ex and the pressure drop for a single channel results in: mean flow velocity U ex = − 1 dP Rh2 , (2) U 2µ dx where P is the pressure, µ the fluid viscosity and x the mean flow coordinate. The hydraulic radius Rh , which is defined as the ratio of fluid volume and wetted surface, can be derived from the radius R of the pipe by: Rh =
R πR 2 = . 2 πR 2
(3)
eh and a length L, Eq. (2) can be rewritten Now assuming a bundle of tubes with an average hydraulic radius R as: e2 . ex = − 1 1P R (4) U 2µ 1L h Neglecting the underlying channel-geometry, one can try to use this formula as a general expression for calculating the pressure drop in other types of geometries, assuming that the average hydraulic radius is known eh can be ep . R or can be derived. For example, consider a porous media built up of spheres of an average diameter D e written as a function of Dp and its porosity ε: ep D ε . 6 (1 − ε) Inserting this expression into Eq. (4) yields: eh = R
e2 D ε2 ex = 1 U0 = 1 1P p , U ε 2µ 1L 36 (1 − ε)2 where U0 is the so called ‘effective velocity’ inside the porous media.
(5)
(6)
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
Eq. (6) can be rewritten in so called Ergun coordinates (see, e.g., [4]), ep ep ρ U0 D ε3 1P D . 72 = 1L ρU02 (1 − ε) µ(1 − ε)
249
(7)
Introducing the dimensionless quantities Reynolds number Re and friction factor f as defined by ep ε3 1P D , (8) 2 1L ρU0 (1 − ε) ep ρ U0 D , (9) Re = µ(1 − ε) the advantage of using the Ergun coordinates becomes obvious, because Eq. (7) can be written in the compact form Λth , (10) f= Re with the ‘theoretical’ friction coefficient Λth = 72. f=
2.2. Experimental results The friction coefficient was measured by Durst [1] experimentally for a packed bed of spheres with different diameters at a wide range of Reynolds numbers (see Fig. 1). For the Reynolds numbers below Re = 1, the experimentally determined friction coefficient appears to be a constant with Λexp = 182. This value is about 2.5 times higher than the theoretical one derived within the capillaric theory. To explain this additional pressure loss, it is usually argued, that the capillaric theories do not take into account the complex paths, the fluid normally has to go through a porous media. To take into account this longer flow paths compared to the edge length of a porous media, the so called tortuosity factor is introduced as follows (see, e.g., [4]): τ=
length of flow paths . macroscopic length scale
(11)
It might be doubted that a tortuosity factor of τ = 2–3 is a realistic assumption, because this would imply that the length of the fluid channels is up to three times larger than the length of the porous media. Also, there exists
Fig. 1. Friction coefficient versus Reynolds number (section from [1]).
250
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
no way to measure the tortuosity directly. Some a posteriori derivation based on comparing the measured pressure drop to the one derived theoretically seems to be a highly questionable – although often practiced – method. In the following section, we will show that there exists another, and until now, by the capillaric theories not recognized physical effect causing pressure drop in porous media flow.
3. Numerical simulation of a porous media flow For the analysis of the pressure drop of a flow through a packed bed of spheres, we used the lattice Boltzmann method described shortly below to compute the detailed flow field. This flow field was then analyzed with respect to the elongational stresses, which can produce additional pressure losses. 3.1. Numerical method For the computations, a 3D 19-speed (D3Q19) lattice Boltzmann automata with single time Bhatnager–Gross– proposed by Qian [5] is used: Krook (BGK) relaxation collision operator 4Boltz i , Ni (t∗ + 1, rE∗ + cEi ) = Ni (t∗ , rE∗ ) + 4Boltz i 4Boltz i
eq = ω(Ni
(12)
− Ni ),
(13) eq Ni :
with a local equilibrium distribution function ciα uα uα uβ ciα ciβ eq − δ . Ni = tp % 1 + 2 + αβ cs 2cs2 cs2
(14)
eq
This local equilibrium distribution function Ni is chosen in such a way to recover the time-dependent Navier– Stokes equation [5] in the low mach number limit: ∂t ρ + ∂α (ρuα ) = 0,
∂t (ρuα ) + ∂β (ρuα uβ ) = −∂α p + µ∂β ∂β uα + ∂α uβ .
(15) (16)
3.2. Boundary conditions A parabolic velocity inlet profile and a fixed pressure at the outlet were chosen for all test cases. This was achieved by introducing the equilibrium density distribution at the first and last lattice column computed with an upstream extrapolated pressure and a downstream extrapolated flow velocity for the inlet and outlet respectively. The inlet and outlet region was chosen to be long enough to prohibit any errors introduced by this method from affecting the measured quantities. A slightly improved bounce back boundary condition for the wall boundaries (including one relaxation step per iteration also on the solid domain nodes) ensures an accurate and relaxation parameter independent position of the wall for a range of at least 1.8 6 ω 6 1.95. 3.3. Validation The code used for this investigation was already validated in detail for time dependent flow around simple geometries [8]. To produce quantitatively reliable data when using CFD, it is also necessary to ensure the grid independence of the numerical results. This is usually done by discretizing the same geometry with meshes of increasing sizes and observing the convergence of the results with increasing mesh refinement. Especially a discretization of spherical objects with rectangular elements, as it is common practice within the LB technique, makes it necessary to carefully investigate the discretization error.
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
251
Fig. 2. Packed bed of spheres, surface shaded with the pressure.
Fig. 3. Friction coefficient versus particle diameter.
As a validation test-case, the pressure drop for low Reynolds number flow through an orthorhombic package of spheres was simulated. The domain sizes were chosen to be lx × ly × lz = 120 × 20 × 20 lattice nodes for the 8 × 2 × 2 spheres with a diameter of Dp = 10 for the coarsest resolution and lx × ly × lz = 480 × 80 × 80 lattice nodes with Dp = 40 for the finest resolution. Periodic boundary conditions were applied normal to the main flow direction to make this test case similar to the experimental set up of Durst et al. [1]. A very good convergence of the numerically achieved friction coefficient to the experimental values (see Fig. 1) can be observed in Fig. 3. A sphere diameter of Dp = 20 is sufficient to approach the convergence result within an accuracy of bet ter than 3%. 3.4. Flow simulation through porous media In order to have realistic geometries for the pressure drop studies, two porous media from engineering applications were chosen: One sponge-like SiC matrix, and one catalyst consisting of a tube filled with spheres.
252
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
Fig. 4. Computer tomography data as input for the LB simulation: catalyst.
3.4.1. Geometry preprocessing For both samples, the catalyst as well as the SiC matrix, the geometry was digitized using 3D computer tomography (3D CT). The tomography data can directly be used as geometry input for the lattice Boltzmann simulation, without the need of any further preprocessing [6,7]. 3.4.2. Case 1: Catalyst A cylindrical porous probe with a height of 110 mm and a diameter of 80 mm was scanned using 3D CT with an average resolution of 0.9 mm. This leads to a discretization of lx × ly × lz = 123 × 90 × 90 voxels. With these parameters, the average diameter of the spheres is in a range which allows one to predict correct results according to the convergence study in the foregoing section. The complex geometry data were centered inside an lx × ly × lz = 250 × 99 × 99 sized channel (see Fig. 4), and a flow for a Reynolds number of about Re ≈ 0.1 was simulated using velocity inlet and pressure outlet boundary conditions. The simulation was performed on one processor of a VPP 700 at the Leibniz–Rechenzentrum in Munich; 50,000 iterations were necessary for this set-up, which took about 25,200 CPU seconds, 850 Mbyte of computer memory were necessary for the storage of the ≈ 2.45 × 106 voxels. 3.4.3. Case 2: SiC matrix A cylindrical porous probe with a height of 30 mm and a diameter of 82 mm was scanned using 3D CT with an average resolution of 0.5 mm. This leads to a discretization of lx × ly × lz = 44 × 147 × 147 voxels. The average diameter of the flow channels is again large enough with this resolution to produce lattice-independent results. The complex geometry data were centered inside an lx × ly × lz = 100 × 149 × 149 sized channel, and a flow for a Reynolds number of about Re ≈ 0.1 was simulated using velocity inlet and pressure outlet boundary conditions. The simulation was performed on one processor of a VPP 700 at the Leibniz–Rechenzentrum in Munich; 10,000 iterations were necessary for this set-up, which took about 5760 CPU seconds, 800 Mbyte of computer memory were necessary for the storage of the ≈ 2.2 × 106 voxels.
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
253
Fig. 5. Pressure field for the flow through a catalyst (x–z plane).
Fig. 6. Computer tomography data as input for the LB simulation: SiC matrix (right: section).
4. Analysis of the pressure drop from experimental data As argued above, the tortuosity is obviously not the only reason for the higher pressure drop observed in experiments and numerical simulations when compared to the results derived from the capillaric theories. The total dissipation in the flow when passing through a porous media can be expressed by: Φ = −τij with
∂Uj , ∂xi
∂Uj ∂Ui + −τij = −µ ∂xi ∂xj
(17)
∂Uj 2 , − µδij 3 ∂x | {z i} =0
for incompressible fluid. Eq. (17) can be rewritten as ∂Uj ∂Ui ∂Uj + Φ = −µ ∂xi ∂xj ∂xi
(18)
254
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
Fig. 7. Pressure field for the flow through a porous SiC matrix (left: x-z plane, right: istotache, shaded with the pressure).
z = −µ
Φe (elongation)
∂u1 ∂x1
2
}| { ∂u2 2 ∂u3 2 ∂u3 ∂u2 2 ∂u1 ∂u3 2 ∂u2 ∂u1 2 + + + + + + + . + ∂x2 ∂x3 ∂x ∂x2 ∂x2 ∂x3 ∂x3 ∂x1 {z } | 1
Φs (shear)
(19) The dissipation can thus be expressed as a sum of two parts: the dissipation caused by shear forces Φs and the dissipation caused by elongational strain Φe : Φ = Φs + Φe .
(20)
An evaluation of the detailed flow fields produced by the numerical simulation with the fraction, Φe /Φs , yields: Example 1 (Catalyst): Friction coefficient
Λ = 217.6,
Dissipation (elongation) Dissipation (shear) elongation / shear
Φe = 2.111 ∗ 10−02 , Φs = 2.854 ∗ 10−02 , Φe /Φs = 0.74.
Example 2 (SiC matrix): Friction coefficient
Λ = 342.4,
Dissipation (elongation) Dissipation (shear) elongation / shear
Φe = 3.853 ∗ 10−04 , Φs = 6.486 ∗ 10−04 , Φe /Φs = 0.59.
J.Bernsdorf et al. / Computer Physics Communications 129 (2000) 247–255
255
In both examples, the friction coefficient Λ is much larger than predicted by the capillaric theory. A considerable amount of the pressure drop is caused by the elongational strain Φe , what can be clearly observed from the relation Φe /Φs . In the derivation of the Kozeny–Darcy law (as well as in similar approaches), only the shear forces are taken into account. The gap for the pressure drop between the theoretical values due to these shear forces and the experimental results is usually explained by the tortuosity of the porous media, while the contribution of the elongational forces is neglected. Although for the above examples tortuosity causes a fraction of the observed pressure drop, it could be shown that this is not the only reason for the difference between the experimental and the theoretical results. So, an a posteriori adjustment of the tortuosity with a much too large tortuosity factor of τ ≈ 2–3, while neglecting the contribution of the elongational forces, may lead to a considerable error.
5. Conclusion With the low Reynolds number flow simulation in regularly packed beds of spheres we could show that the lattice Boltzmann method is able to quantitatively predict pressure drop with high accuracy. Three-dimensional computer tomography was used to digitize the geometry of two porous media: a catalyst and a porous SiC matrix. It could be shown by the a posteriori analysis of the simulated flow fields, that elongational strain gives an important contribution to the pressure drop. This must not be neglected when deriving a theory for predicting the pressure drop in porous media flow. The derivation of a tortuosity factor from pressure drop measurements while neglecting the dissipation due to elongational strain might produce a considerable error.
Acknowledgements The authors would like to thank the Hattinger Prüf- und Entwicklungs-GmbH (HAPEG) for providing the computer tomography data. This project was founded by the Deutsche Forschungsgemeinschaft (Proj. Nr. Br 1864/1). References [1] [2] [3] [4] [5] [6] [7] [8]
F. Durst, R. Haas, W. Interthal, J. Non-Newtonian Fluid Mech. 22 (1987) 169. D.H. Rothman, Geophysics 53 (1988) 509. J. Bernsdorf, M. Schäfer, F. Durst, Int. J. Numer. Meth. Fluids 29 (1999) 251. R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena (New York, 1960). Y.H. Qian, Europhys. Lett. 17 (6) (1992) 479. B. Ferréol, D.H. Rothman, Transport in Porous Media 20 (1995) 3. J. Bernsdorf, O. Günnewig, W. Hamm, O. Münker, GIT Labor-Fachzeitschrift 4/99 (1999) 387. J. Bernsdorf, Th. Zeiser, G. Brenner, F. Durst, Int. J. Mod. Phys. C 9 (8) (1998) 1141.