Bingham Fluid Simulation in Porous Media with Lattice Boltzmann Method

Abstract

Generate a lattice Boltzmann model (LBM), which allows to simulate the behavior of a Bingham fluid through a rectangular channel with the D2Q9 model. For this purpose, a relaxation parameter is proposed based on the rheological parameters of the Bingham model. The validation will be carried out with the solution of the movement equation, and velocity profiles will be obtained for three different Bingham numbers (Bn). Other simulations will be made in a rectangular channel in the presence of arbitrarily and randomly generated porous media. The main objective is to propose a method to predict the behavior of non-Newtonian fluids (Bingham fluid) through porous media, which have many applications in the chemical industry avoiding at the same time the expensive experimentation of these systems, with predicting models.

1. Introduction

A continuous medium is characterized by the fact that its atoms or molecules are so close together, in such way that they could be considered macroscopically as a homogeneous mass, whose behavior can be foreseen without considering the movement of each of its elementary particles that compose it. In this sense it is assumed that there are no gaps or separations between the particles.

The movement of fluids can have a wide variety of behaviors for both simple and complex flows (biological and food systems). In addition to this, using Reynolds number, we can know if the flow has turbulent or laminar regime.

Mass conservation is the basic principle of fluid movement, which requires that when the fluid is in motion, it moves so that the mass is preserved. The movement of fluids is governed in general by the continuity equation.

$$\frac{\mathrm{\partial }\mathrm{\rho }}{\mathrm{\partial }\mathrm{t}}+\nabla \cdot \left(\mathrm{\rho }\overrightarrow{\mathrm{v}}\right)=0$$E1

On the other hand, the Navier-Stokes equation, which in general terms corresponds to the application of Newton’s second law of classical mechanics to fluid movement, is described as follows:

$$\mathrm{\rho }\left[\frac{\mathrm{\partial }\overrightarrow{\mathrm{v}}}{\mathrm{\partial }\mathrm{t}}+\overrightarrow{\mathrm{v}}\nabla \cdot \overrightarrow{\mathrm{v}}\right]=-\nabla \mathrm{P}+\mathrm{\mu }{\nabla }^2\overrightarrow{\mathrm{v}}+\mathrm{\rho }\overrightarrow{\mathrm{g}}$$E2

for a fluid with velocity $\overrightarrow{\mathrm{v}}$, density ρ, pressure p, and kinematic shear viscosity μ. The Navier–Stokes equation is a second-order partial differential equation, which can have an analytical solution only for a small number of cases.

There are basically three types of fluids from the point of view of fluid dynamics, which are called Newtonians, non-Newtonians, and viscoelastic.

1.1 Newtonian fluids

In this type of fluid, the shear stress or shear force per unit area is proportional to the viscosity gradient:

$${\mathrm{\tau }}_{\mathrm{yx}}=\mathrm{\mu }\left(\frac{{\mathrm{dv}}_{\mathrm{x}}}{\mathrm{dy}}\right)$$E3

In Eq. (3), τ_yx is the shear stress, (dv_x/dy) is the shear rate, and μ is the viscosity, and its graphic representation can be seen in Figure 1.

Figure 1.

Types of times-independent flow behavior [2].

All gases, liquid water, and liquids of single molecules (ammonia, alcohol, benzene, petroleum, chloroform, butane, etc.) are Newtonian. Many food materials such as milk, apple juice, orange juice, wine, and beer have a Newtonian behavior [1].

1.2 Non-Newtonian fluids

When the relationship between shear stress and shear rate is not linear, the fluid is called non-Newtonian. There are many these types of fluids, and their behaviors are shown in Figure 1.

In the case of these fluids, the viscosity is no longer constant; therefore, the relationship between the shear stress and the shear rate of the fluid is no longer linear. For this reason, a new term is now introduced and is known as apparent viscosity or also known as shear rate-dependent viscosity [3].

Some fluids have a yield stress, from which the fluid begins to move. Below this tension the shear rate would be zero. This relationship is not linear or, if it is, it does not pass through the origin [4]. Complex mixtures are considered non-Newtonian fluids: grouts, pastes, gels, polymer solutions, etc. Most non-Newtonian fluids are mixtures with constituents of very different sizes. For example, toothpaste is composed of solid particles suspended in an aqueous solution of several polymers. Solid particles are much larger than water molecules, and polymer molecules are much larger than water molecules.

Much of the research that is carried out in the field of non-Newtonian fluids has been focused in the measure of its shear stress-shear strain curves and to look for mathematical descriptions of these curves. The study of the behavior of the flow of materials is called rheology (a term that originates from Greek words that give the meaning of “the study of flow”); thus, diagrams such as the one shown in Figure 1 are often called rheograms.

In the case of Bingham fluids, sometimes called Bingham plastics, they resist a small shear force indefinitely, but they flow easily under large shear stresses. In other words, at low stresses the plastic viscosity of Bingham is infinite, and at greater stresses the viscosity decreases with the increase in the velocity gradient. Examples are bread dough, toothpaste, apple sauce, some paints, plastics, mayonnaise, ketchup, aleas, and some grouts [5, 6].

2. Bingham model

Eugene C. Bingham described the paintings with this model in 1919, published in his book Fluidity and Plasticity in 1922. The model was analyzed by Oldroyd (1947), Reiner (1958), and Prager (1961).

The main feature of the Bingham model is the yield stress, necessary for the fluid to deform or flow. Above this minimum yield stress, the fluid begins to move. If this yield stress is not exceeded, the fluid behaves like a rigid or quasi-rigid body, with zero shear rate.

The relationship between the shear stress and the velocity gradient is linear, but it does not go through the origin for a Bingham plastic (Figure 1); its mathematical model is given by

$${\mathrm{\tau }}_{\mathrm{yx}}={\mathrm{\tau }}_0+{\mathrm{\mu }}_{\mathrm{B}}\stackrel{\dot{\mathrm{\gamma }}}{\overbrace{\frac{{\mathrm{dv}}_{\mathrm{x}}}{\mathrm{dy}}}}\text{if}\left|{\mathrm{\tau }}_{\mathrm{yx}}\right|>{\mathrm{\tau }}_0$$E4

$$\frac{{\mathrm{dv}}_{\mathrm{x}}}{\mathrm{dy}}=0;\text{if}\left|{\mathrm{\tau }}_{\mathrm{yx}}\right|<{\mathrm{\tau }}_0$$

$$\mathrm{\mu }\left(\dot{\mathrm{\gamma }}\right)={\mathrm{\mu }}_{\mathrm{B}}+\frac{{\mathrm{\tau }}_0}{\dot{\mathrm{\gamma }}}\text{for}{\mathrm{\tau }}_{\mathrm{yx}}>{\mathrm{\tau }}_0$$E5

$$\dot{\mathrm{\gamma }}=0;\text{for}{\mathrm{\tau }}_{\mathrm{yx}}\le {\mathrm{\tau }}_0$$

where τ₀ is the yield stress, μ_B is the plastic viscosity of Bingham, and $\mathrm{\mu }\left(\dot{\mathrm{\gamma }}\right)$ is the apparent viscosity, which decreases with the increase in the magnitude of the shear rate $\dot{\mathrm{\gamma }}$.

Other examples of Bingham-type fluids in foods are tomato sauce, whipped cream, whipped egg white, margarine, and mustard-type condiments [7, 8].

2.1 Analytical solution of the equation of motion in the case of a Bingham-type fluid

In this section, we will obtain a mathematical expression that shows how to quantify the velocity profile for a non-Newtonian fluid for the Bingham model.

In the case of a Poiseuille flow, the flow of a Bingham-type fluid in the x-direction is considered, between two plates separated by a distance 2H, taking into account the steady-state conditions, constant cross section, absence of gravitational effects and isothermal flow, and being incompressible, such as the one shown in Figure 2.

Figure 2.

Flow of a Bingham fluid between two plates one half view [2].

From the equation of motion, we can obtain the stress profile, as well as the velocity profile:

$$\mathrm{\rho }\left[\frac{\mathrm{\partial }\overrightarrow{\mathrm{v}}}{\mathrm{\partial }\mathrm{t}}+\overrightarrow{\mathrm{v}}\cdot \nabla \overrightarrow{\mathrm{v}}\right]=\nabla \cdot \overline{\overline{\mathrm{\tau }}}-\nabla \mathrm{P}+\mathrm{\rho }\overrightarrow{\mathrm{g}}$$E6

The pressure gradient effect is considered as a favorable driving force for fluid movement. Usually, the pressure decreases at a constant rate from the initial end to the end in the x direction.

$$\nabla \mathrm{P}\cong \frac{\mathrm{\Delta P}}{\mathrm{\Delta z}}=\frac{\left({\mathrm{P}}_{\mathrm{L}}-{\mathrm{P}}_0\right)}{\mathrm{L}}=\frac{\mathrm{\Delta P}}{\mathrm{L}}$$E7

The component of the flow density tensor of the amount of movement of Eq. (6) is τ_yx; therefore, considering the above conditions, we have the following equation to solve:

$$\frac{\mathrm{d}{\mathrm{\tau }}_{\mathrm{yx}}}{\mathrm{dy}}=\left(\frac{\mathrm{\Delta P}}{\mathrm{L}}\right)$$E8

Making a separation of variables, in addition to performing the corresponding integrals in Eq. (8), we obtain

$${\mathrm{\tau }}_{\mathrm{yx}}=\left(\frac{\mathrm{\Delta P}}{\mathrm{L}}\right)\mathrm{y}+{\mathrm{c}}_1$$E9

The integration constant is zero, when τ_yx = 0 at y = 0, i.e., at the center of the plates, the shear stress is minimal. Therefore, the shear stress profile is

$${\mathrm{\tau }}_{\mathrm{yx}}=\left(\frac{\mathrm{\Delta P}}{\mathrm{L}}\right)\mathrm{y}$$E10

Equating Eq. (10) with the Bingham model, we obtain

$$\stackrel{\text{Bingham}\text{model}}{\overbrace{{\mathrm{\mu }}_{\mathrm{B}}\left(\frac{{\mathrm{dv}}_{\mathrm{x}}}{\mathrm{dy}}\right)+{\mathrm{\tau }}_0}}=\left(\frac{\mathrm{\Delta P}}{\mathrm{L}}\right)\mathrm{y}$$E11

Making a second separation of variables, in addition to performing the corresponding integrals, we have

$${\mathrm{v}}_{\mathrm{x}}=\left(\frac{\mathrm{\Delta P}}{\mathrm{L}}\right)\frac{1}{{\mathrm{\mu }}_{\mathrm{B}}}\frac{{\mathrm{y}}^2}{2}-\frac{{\mathrm{\tau }}_0}{{\mathrm{\mu }}_{\mathrm{B}}}\mathrm{y}+{\mathrm{c}}_2$$E12

In Figure 2, it is shown that velocity is zero on the plates; i.e., v_x = 0 at y = ±H. Using this condition, the value of c₂ is obtained. Substituting in Eq. (12) we obtain

$${\mathrm{v}}_{\mathrm{x}}=\frac{{\mathrm{\tau }}_0}{{\mathrm{\mu }}_{\mathrm{B}}}\mathrm{H}\left(1-\frac{\mathrm{y}}{\mathrm{H}}\right)+\left(\frac{\mathrm{\Delta P}}{\mathrm{L}}\right)\frac{{\mathrm{H}}^2}{2{\mathrm{\mu }}_{\mathrm{B}}}\left[{\left(\frac{\mathrm{y}}{\mathrm{H}}\right)}^2-1\right]$$E13

To know the velocity profile in the region of the plug flow, the condition for the yield stress is proposed according to Eq. (10), when y = y₀:

$$\mathrm{y}={\mathrm{y}}_0;{\mathrm{\tau }}_0=\left(\frac{\mathrm{\Delta P}}{\mathrm{L}}\right){\mathrm{y}}_0$$E14

Substituting Eq. (14) into Eq. (13), the velocity in the plug flow region is obtained:

$${\mathrm{v}}_0=\left(\frac{\mathrm{\Delta P}}{\mathrm{L}}\right)\frac{{\mathrm{H}}^2}{2{\mathrm{\mu }}_{\mathrm{B}}}{\left(1-\frac{{\mathrm{y}}_0}{\mathrm{H}}\right)}^2$$E15

Commonly, velocity profiles are usually represented with on the Bingham number, which is defined as

$$\mathrm{Bn}=\frac{{\mathrm{\tau }}_0}{{\mathrm{\mu }}_{\mathrm{B}}}\frac{\mathrm{H}}{\mathrm{v}}$$E16

In Eq. (16), v is a characteristic velocity. Dividing Eq. (13) by this velocity,

$$\frac{{\mathrm{v}}_{\mathrm{x}}}{\mathrm{v}}=\mathrm{Bn}\left(1-\frac{\mathrm{y}}{\mathrm{H}}\right)+\left(\frac{\mathrm{\Delta P}}{\mathrm{L}}\right)\frac{{\mathrm{H}}^2}{2\mathrm{v}{\mathrm{\mu }}_{\mathrm{B}}}\left[{\left(\frac{\mathrm{y}}{\mathrm{H}}\right)}^2-1\right]$$E17

is obtained. Eq. (17) represents the velocity profile for the Poiseuille flow between two parallel plates, in the case of a Bingham-type fluid, and its graphic representation is that shown in Figure 3, for values of Bn = 0.1, 0.2, 0.3, and 0.4, with dimensionless values [4].

Figure 3.

Velocity profile for a Bingham fluid; Bn = 0.1, 0.2, 0.3, and 0.4 [4].

The graph of the shear stress vs. normalized shear rate (rheogram) is shown in Figure 4.

Figure 4.

Shear stress vs. normalized shear rate with analytical solution.

3. Lattice Boltzmann Bhatnagar-Gross-Krook (BGK)

The Lattice Boltzmann Method (LBM) generally consists of a discrete lattice; each site (node) is represented by the distribution function, which is defined by the velocity of a particle and limited to a discrete group of allowed velocities. During each time step, the movement or jump of particle to nearby lattice sites, along its direction of movement, a collision with another particle can occur when they reach the same site. The result of the collisions is determined by means of the solution of the kinetic equation of Boltzmann for the new function of distribution of the particle to that site, and in this way the function of distribution of the particle is updated [9, 10].

There are different lattice models, which are given by DmQn; m indicates the dimension and n the permitted velocity, thus, the D2Q9 model (two-dimensional with nine speed directions), of which four sites correspond to nearby neighbors (right, left, up, and down), four other points of the lattice vectors along the diagonal faces of the following sites, along these diagonals. In this way the particles can travel in eight directions for each lattice site. The circle in the center of the square represents the vector, which has a value of zero and represents particles that have no movement, that is, particles at rest. Then, there are a total of nine real numbers that describe the distribution function of the particle for a lattice site (see Figure 5).

Figure 5.

Allowed directions for particle movement. (a) Model D2Q9, (b) model D2Q15.

The process of the propagation and collision of particles generally occurs in two stages: the first is to denote the advance of the particles to the next lattice site along the directions of movement; this is for each time step Δt. The second stage is to simulate the collisions of the particles, which causes them to propagate in different directions at each lattice site [11, 12]. These stages can be described through the discretized Boltzmann equation on a lattice.

$${\mathrm{f}}_{\mathrm{i}}\left(\overrightarrow{\mathrm{x}}+{\overrightarrow{\mathrm{e}}}_{\mathrm{i}},\mathrm{t}+1\right)-{\mathrm{f}}_{\mathrm{i}}\left(\overrightarrow{\mathrm{x}},\mathrm{t}\right)=-\mathrm{\Omega }$$E18

To simplify Eq. (18), the BGK approximation is usually used; this approximation replaces the term Ω:

$$\mathrm{\Omega }=\frac{{\mathrm{f}}^{\mathrm{eq}}-\mathrm{f}}{\mathrm{\tau }}$$E19

This operator models the effect of the collision as a relaxation of the distribution function towards the Maxwell equilibrium distribution. The parameter τ (relaxation time) has dimensions of time and controls the frequency with which the distribution function relaxes to reach equilibrium, that is, this time determines the rate at which the fluctuations in the system lead to this state of equilibrium [13].

$${\mathrm{f}}_{\mathrm{i}}\left(\overrightarrow{\mathrm{x}}+{\overrightarrow{\mathrm{e}}}_{\mathrm{i}},\mathrm{t}+1\right)-{\mathrm{f}}_{\mathrm{i}}\left(\overrightarrow{\mathrm{x}},\mathrm{t}\right)=-\frac{1}{\mathrm{\tau }}\left[{\mathrm{f}}_{\mathrm{i}}\left(\overrightarrow{\mathrm{x}},\mathrm{t}\right)-{\mathrm{f}}_{\mathrm{i}}^{\mathrm{eq}}\left(\overrightarrow{\mathrm{x}},\mathrm{t}\right)\right]$$E20

The macroscopic variables (mass (ρ) and velocity $\left(\overrightarrow{\mathrm{u}}\right)$) can be calculated directly from the values of the distribution function as

$$\begin{array}{l}\mathrm{\rho }\left(\overrightarrow{\mathrm{x}},\mathrm{t}\right)=\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{f}}_{\mathrm{i}}\left(\overrightarrow{\mathrm{x}},\mathrm{t}\right)\\ \overrightarrow{\mathrm{u}}\left(\overrightarrow{\mathrm{x}},\mathrm{t}\right)=\frac{1}{\mathrm{\rho }\left(\overrightarrow{\mathrm{x}},\mathrm{t}\right)}\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{f}}_{\mathrm{i}}\left(\overrightarrow{\mathrm{x}},\mathrm{t}\right){\overrightarrow{\mathrm{e}}}_{\mathrm{i}}\end{array}$$E21

In the case of a Newtonian fluid, the ratio of kinematic viscosity and relaxation time is given by [13] υ = 1/3 (τ-1/2).

4. Simulation of Bingham fluids with the lattice Boltzmann method

For the simulations of the Bingham fluid with the Lattice Boltzmann Method, a modification was made to the LBGK approach presented by Wang and Ho [14] and Tang et al. [15], for a D2Q9 model, which consists in proposing the relaxation parameter τ based on the apparent viscosity:

$$\mathrm{\tau }=3\left[{\mathrm{\mu }}_{\mathrm{B}}+{\mathrm{\tau }}_0\left(\dot{\mathrm{\gamma }}\right)\right]+\frac{1}{2}$$E22

In Eq. (22), ${\mathrm{\mu }}_{\mathrm{B}}$, ${\mathrm{\tau }}_0$, and $\dot{\mathrm{\gamma }}$ are the Bingham viscosity, yield stress, and shear rate, respectively. Parameter τ was used in the Lattice Boltzmann equation (Eq. (20)). The simulations were carried out on a 64 × 64 lattice, using “bounce back” conditions on the solid walls to ensure that the velocities are zero and periodic boundary conditions at the fluid inlet and outlet, so that the nodes located in the border will have their neighboring nodes on the opposite border. The steady state was reached at 360,000 time steps.

The validation of the proposal in the LBM was performed by comparing the results of the analytical solutions for a Poiseuille flow between two separate plates a distance 2H, shown in Figure 2 using Eq. (17). The used conditions were pressure force = 2.66E-2, yield stress = 2.00E-5, and Bingham viscosity of 0.4 for a Bin = 0.1; pressure force = 5.83E-3, yield stress = 1.10E-5, and Bingham viscosity of 0.08 for a Bin = 0.2; pressure force = 5.19E-3, yield stress = 1.40E-5, and Bingham viscosity of 0.07 for a Bin = 0.3; and pressure force = 1.88E-3, yield stress = 6.50E-6, and Bingham viscosity of 0.025 for a Bin = 0.4.

Simulations were performed in porous media, applying the LBM in the case of deterministic porous medium and random porous media. A modification of a “Box-Muller method,” which is a random number generator, was inserted in random porous media, and blocks were inserted arbitrarily in the lattice for deterministic porous media [16, 17].

An alternative way or method is proposed for obtaining local permeabilities for deterministic and random porous media. This one consists in a modification of Darcy’s law; for this the permeability is calculated based on the apparent viscosity according to the following equation:

$$\mathrm{K}=\left({\mathrm{\mu }}_{\mathrm{B}}+\frac{{\mathrm{\tau }}_0}{\dot{\mathrm{\gamma }}}\right)\mathrm{v}/\left(\raisebox{1ex}{$\mathrm{dP}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dx}$}\right.\right)$$E23

In Eq. (23), K is the permeability, v is the velocity, and dP/dx is a pressure force.

The discrete macroscopic pressure P is given by the state equation that relates the discrete density to the pressure $\mathrm{P}={\mathrm{c}}_{\mathrm{s}}^2\mathrm{\rho }$, where c_s is the speed of sound and ρ is the density that is calculated through Eq. (21).

The relationship is valid for incompressible fluid simulations and is only allowed to fluctuate locally around a fixed value [18]. Hidemitsu Hayashi proposed two LBMs for the flow generated by the pressure gradient (FGPG) and the flow driven by an external force (FDEF), which are consistent with each other [19].

The criterion used by Wang and Ho [14] was taken, for which yielding occurs when the magnitude of the extra shear stress tensor exceeded the yield stress, ${\mathrm{\tau }}_0$, i.e., be yielded when $\left|{\mathrm{\tau }}_{\mathrm{yx}}\right|⟩{\mathrm{\tau }}_{\mathrm{o}}$ and unyielded if $\left|{\mathrm{\tau }}_{\mathrm{yx}}\right|\le {\mathrm{\tau }}_{\mathrm{o}}$.

Figure 6 shows the validation of the velocity profiles with the analytical solution and the simulations with LBM. The error between both solutions was less than 2.0%.

Figure 6.

Comparison of velocity profiles with different Bingham numbers for the analytical solutions and the proposed LBM. Normalized velocity profiles for (a) Bin = 0.1, (b) Bin = 0.2, (c) Bin = 0.3, (d) Bin = 0.4.

For the development of the work, three porous media with a deterministic structure and nine random were proposed; in each of them all the simulations were performed for three Bingham numbers (0.2, 0.3, and 0.4).

In Figures 7–15, the speed patterns for all simulations are shown; for this the values of the yield stress, the pressure forces, and the viscosities were varied.

Figure 7.

Bingham = 0.2, yield stress = 1.1E-5, pressure force = 5.83E-3, and viscosity = 0.08. Porosity 81.68%, deterministic porous medium: (a) vectorized flow, (b) velocity patterns. Porosity 81.62%, random porous medium: (c) vectorized flow, (d) velocity patterns.

Figure 8.

Bingham = 0.2, yield stress = 1.1E-5, pressure force = 5.83E-3, and viscosity = 0.08. Porosity 73.75%, deterministic porous medium: (a) vectorized flow, (b) velocity patterns. Porosity 73.68%, random porous medium: (c) vectorized flow, (d) velocity patterns.

Figure 9.

Bingham = 0.2, yield stress = 1.1E-5, pressure force = 5.83E-3, and viscosity = 0.08. Porosity 65.82%, deterministic porous medium: (a) vectorized flow, (b) velocity patterns. Porosity 65.75%, random porous medium: (c) vectorized flow, (d) velocity patterns.

Figure 10.

Bingham = 0.3, yield stress = 1.4E-5, pressure force = 5.19E-3, and viscosity = 0.07. Porosity 81.68%, deterministic porous medium: (a) vectorized flow, (b) velocity patterns. Porosity 81.62%, random porous medium: (c) vectorized flow, (d) velocity patterns.

Figure 11.

Bingham = 0.3, yield stress = 1.4E-5, pressure force = 5.19E-3, and viscosity = 0.07. Porosity 73.75%, deterministic porous medium: (a) vectorized flow, (b) velocity patterns. Porosity 73.68%, random porous medium: (c) vectorized flow, (d) velocity patterns.

Figure 12.

Bingham = 0.3, yield stress = 1.4E-5, pressure force = 5.19E-3, and viscosity = 0.07. Porosity 65.82%, deterministic porous medium: (a) vectorized flow, (b) velocity patterns. Porosity 65.75%, random porous medium: (c) vectorized flow, (d) velocity patterns.

Figure 13.

Bingham = 0.4, yield stress = 6.5E-6, pressure force = 1.88E-3, and viscosity = 0.025. Porosity 81.68%, deterministic porous medium: (a) vectorized flow, (b) velocity patterns. Porosity 81.62%, random porous medium: (c) vectorized flow, (d) velocity patterns.

Figure 14.

Bingham = 0.4, yield stress = 6.5E-6, pressure force = 1.88E-3, and viscosity = 0.025. Porosity 73.75%, deterministic porous medium: (a) vectorized flow, (b) velocity patterns. Porosity 73.68%, random porous medium: (c) vectorized flow, (d) velocity patterns.

Figure 15.

Bingham = 0.4, yield stress = 6.5E-6, pressure force = 1.88E-3, and viscosity = 0.025. Porosity 65.82%, deterministic porous medium: (a) vectorized flow, (b) velocity patterns. Porosity 65.75%, random porous medium: (c) vectorized flow, (d) velocity patterns.

In Figures 7–15 the decrease in porosity corresponds to a decrease in velocity. It can be noted that for a Bingham number of 0.2 (Figures 7–9), higher velocities are shown in deterministic porous media than random media. But in the latter, there are more places to pass the fluid. These behaviors are presented in the simulations for the three Bingham numbers used (0.2, 0.2, and 0.3). Comparing the Bingham number of 0.2 with that of 0.4, a decrease in the velocity in the deterministic porous media is observed as well as the random ones according to the conditions handled; this is due to the decrease of the initial effort, the pressure force, and the viscosity.

Local permeabilities were simulated based on apparent viscosities for all deterministic and random porous media. In Figures 16–19, only some of the results of local permeabilities for deterministic porous media 81.68 and 65.82% with Bingham numbers of 0.2 and 0.4, respectively, are shown. Likewise, only the result of the simulations for two random porous media 81.62 and 65.75% for Bingham numbers of 0.2 and 0.4, respectively, are shown.

Figure 16.

Local permeabilities for deterministic porous media with porosity of 81.68%, Bingham = 0.2, yield stress = 1.1E-5, pressure force = 5.83E-3 and viscosity = 0.08. (a–c) The rough numerical date of permeability, (d) pattern porous media.

Figure 17.

Local permeabilities for random porous media with porosity of 81.62%, Bingham = 0.2, yield stress = 1.1E-5, pressure force = 5.83E-3 and viscosity = 0.08; (a–c) The rough numerical date of permeability, (d) pattern porous media.

Figure 18.

Local permeabilities for deterministic porous media with porosity of 65.82%, Bingham = 0.4, yield stress = 6.5E-6, pressure force = 1.88-3 and viscosity = 0.025; (a–c) the rough numerical date of permeability, (d) pattern porous media.

Figure 19.

Local permeabilities for random porous media with porosity of 65.75%, Bingham = 0.4, yield stress = 6.5E-6, pressure force = 1.88-3 and viscosity = 0.025; (a–c) The rough numerical date of permeability, (d) pattern porous media.

Figures 16–19 show the zones in blue of the local permeabilities. It is remarkable that the blue areas predominate in random porous media. By comparing the Bingham number of 0.2 for the two deterministic and random porous media according to Figures 16 and 17 with that of 0.4 of Figures 18 and 19, an increase in permeabilities can be seen.

Finally, pressures for all porous media were simulated. Figures 20–25 show some of the results obtained in the case of deterministic porous media with porosities of 81.68, 73.75, and 65.82% for Bingham numbers 0.2, 0.3, and 0.4, respectively. Similarly, the results are presented for random porous media with porosities of 81.62, 73.68, and 65.75% for Bingham numbers 0.2, 0.3, and 0.4, respectively.

Figure 20.

Pressures for deterministic porous medium with porosity of 81.68%, Bingham = 0.2, yield stress = 1.1E-5, pressure force = 5.83E-3 and Viscosity = 0.08. (a) Pressure distribution, (b) the rough numerical date of pressure.

Figure 21.

Pressures for random porous medium with porosity of 81.62%, Bingham = 0.2, yield stress = 1.1E-5, pressure force = 5.83E-3 and Viscosity = 0.08. (a) Pressure distribution, (b) the rough numerical date of pressure.

Figure 22.

Pressures for deterministic porous medium with porosity of 73.75%, Bingham = 0.3, yield stress = 1.4E-5, pressure force = 5.19E-3 and Viscosity = 0.07. (a) Pressure distribution, (b) the rough numerical date of pressure.

Figure 23.

Pressures for random porous medium with porosity of 73.68%, Bingham = 0.3, yield stress = 1.4E-5, pressure force = 5.19E-3 and Viscosity = 0.07. (a) Pressure distribution, (b) the rough numerical date of pressure.

Figure 24.

Pressures for deterministic porous medium with porosity of 65.82%, Bingham = 0.4, yield stress = 6.5E-6, pressure force = 1.88E-3 and Viscosity = 0.025. (a) Pressure distribution, (b) the rough numerical date of pressure.

Figure 25.

Pressures for random porous medium with porosity of 65.75%, Bingham = 0.4, yield stress = 6.5E-6, pressure force = 1.88E-3 and Viscosity = 0.025. (a) Pressure distribution, (b) the rough numerical date of pressure.

Figures 20–25, you can see in all cases the difference in pressures, higher pressure in the red colors, and less pressure in the green colors, in addition to observing the pressures in the different zones.

5. Conclusions

In the present work, the Lattice Boltzmann Method was applied to a problem of the flow of a non-Newtonian Bingham-type fluid between two plates, in the case of a Poiseuille flow.

This method is an alternative to the conventional ones used in computational fluid mechanics, its programming is not complicated, and today it is applied to many engineering problems.

Validations were carried out with the analytical solution of the velocity profiles for the case of a Poiseuille flow and the simulations with Lattice Boltzmann, for the case of Bingham-type fluids, for values of the Bingham number (Bin) of 0.1, 0.2, 0.3, and 0.4. The results of all the simulations were quite acceptable, since the percentage of error between both results did not exceed 2.0%.

The LBM proves to be kind for simulations with small lattices, such as the one used in the present work 64 × 64. All simulations were performed in a laminar regime.

Three deterministic porous media with porosities of 81.68, 75.75, and 65.82% and nine randomized ones with porosities of 81.62, 73.68, and 65.75% were proposed for three Bingham numbers (0.2, 0.3, and 0.4), to perform all simulations. In them the pressure forces, yield stress, and viscosities were varied.

Profiles of velocities, permeabilities, and local pressures were obtained, in all cases the results and behaviors were acceptable for all porous media, and the three Bingham numbers, although only some of the results obtained, were presented at work.

The LBM with the necessary restrictions allows to perfectly simulate the behavior of fluids, as is the case of the Bingham type; the importance of this is the application of multiple industrial processes, in the displacement of fluids reducing costs and time.

Finally, it would be convenient to perform simulations with turbulent flows to verify the goodness of the method with this type of fluid, in which its description is more complex.

Acknowledgments

The present work was developed under the sponsorship of the Facultad de Estudios Superiores Cuautitlán-Universidad Nacional Autónoma de México.

Conflict of interest

The author declares no conflicts of interest regarding the publication of this paper.