Internally heated convection of viscoplastic fluids in enclosures using a lattice Boltzmann method

Kefayati, Gholamreza

doi:10.1063/5.0139057

Thermal convection driven by an internal heat source in a two-dimensional enclosure filled with viscoplastic fluids is investigated numerically. Two vertical side walls of the cavity are isotherms with the same low temperatures, while the horizontal walls are adiabatic and insulated. An exact Bingham model is applied in the constitutive equation for the viscoplastic fluid. A lattice Boltzmann method (LBM) is developed to solve the introduced non-dimensional macroscopic equations, and the derivations of the LBM are presented and discussed. The implemented LBM is validated against previous studies of internal natural convection. The effects of the Rayleigh–Roberts number, the Prandtl number, the aspect ratio of the cavity, and the inclined angle of the enclosure on the yielded/unyielded parts are investigated and reported. The maximum (or critical) Bingham (Bn) and yield (Y) numbers for the studied parameters are investigated through the defined Nusselt number. The results are depicted by the contours of isotherms, streamlines, yielded/unyielded zones, vorticity, and horizontal and vertical velocities. In addition, the temperatures and velocities in the middle of the cavity as well as the Nusselt number are shown and discussed. It was revealed that the maximum (or critical) yield number is independent of Rayleigh–Roberts and Prandtl numbers same as external natural convection. The values of the critical yield number decrease gradually as the inclined angle rises counterclockwise. However, the critical yield number enhances with the increase in the aspect ratio although the augmentation is not linear and steady.

I. INTRODUCTION

The study of natural convection due to internal heat sources was considered to reveal the earth-mantle convection process.¹ For the internal heating source (Q_i), a modified Rayleigh number was suggested which is also called “Rayleigh–Roberts number” (RR) after the theoretical study of Roberts,² following the experimental investigation of Tritton and Zarraga.³ The temperature difference due to the applied internal heat source is calculated by $Δ θ = \frac{Q_{i} L^{2}}{k}$ ⁠, where k is the thermal conductivity, L is the characteristic length, and Q_i is a constant uniform volumetric energy source.⁴

In the previously mentioned studies, the fluid was assumed to be Newtonian and incompressible flow where the total stress tensor (T) is defined by

T = - p 1 + σ, σ = η I_{1}, I_{1} = grad v + grad v^{T} .

(1)

I_{1}, v

⁠, p,

σ

⁠, and

1

are first Rivlin–Ericksen tensor, velocity vector, pressure, the extra stress tensor, and unit vector; respectively. However, in some of the studied cases, the considered fluid shows a viscoplastic behavior. For example, in the convection of the earth's mantle, viscoplasticity of the material plays an important role and multifarious numerical methods have been applied to solve this problem.⁵ Viscoplastic fluids, in contrast to Newtonian fluids, are divided into solid (unyielded regions) and liquid (yielded regions) phases where the yield stress (σ_y) and its relationship with the second invariants of the extra stress tensor (⁠

ϒ

⁠) are used as a criterion to distinguish the mentioned regimes as

{\begin{matrix} I_{1} = 0, & ϒ \leq σ_{y}, solid phase, \\ σ = η (\dot{γ}, σ_{y}, η_{p}, θ) I_{1}, & ϒ > σ_{y}, liquid phase, \end{matrix}

(2)

where

\dot{γ} = \frac{1}{\sqrt{2}} \sqrt{I_{1} : I_{1}}

and

ϒ = \frac{1}{\sqrt{2}} \sqrt{σ : σ}

⁠. η_p and θ are the dynamic plastic viscosity and temperature, respectively. The known temperature-independent model of Bingham for

η (\dot{γ}, σ_{y}, η_{p})

is defined by

η (\dot{γ}, σ_{y}, η_{p}) = η_{p} + \frac{σ_{y}}{\dot{γ}} .

(3)

For numerical reasons, the constitutive equation without any regulation can be written as⁶

η (\dot{γ}, σ_{y}, η_{p}, λ) = η_{p} + \frac{\sqrt{2} σ_{y} λ}{I_{1}},

(4)

where

λ

is a symmetric traceless tensor and called the “viscoplasticity constraint tensor,” which reveals yielded and unyielded zones (see Appendix A).

External natural convection of viscoplastic fluids, using various analytical, experimental, and numerical models, has been studied widely, in the past two decades. Some of the significant previous studies based on imposed boundary conditions are Rayleigh–Bénard convection, different isotherm temperatures on the walls of enclosures and channels, and applying heat fluxes on boundaries in cavities.^6–12 However, internal natural convection of viscoplastic fluids, which considers the internal heating (and/or cooling) sources, has not been studied yet, based on our best knowledge.

Here, the lattice Boltzmann method (LBM) is used to simulate and study natural convection due to internal heating in an enclosure filled with viscoplastic fluids. The LBM is a robust mesoscopic technique to study a wide range of fluid flows and heat transfer problems.^13–20 To solve this problem, we developed and applied a modified LBM which we introduced and utilized in several isothermal, thermal, and solutal viscoplastic problems before.^21–25

The structure of the paper is arranged in this manner: in Sec. II, the problem statement, the dimensional and non-dimensional governing equations, and different non-dimensional parameters are indicated and explained. The implemented LBM equations to solve the governing equations including the imposed algorithm are defined and discussed in Sec. III. In Sec. IV, the results are demonstrated and analyzed for various considered non-dimensional numbers and their effects on the fluid flow, heat transfer, and yielded/unyielded sections. The main and new findings are noted and determined in Sec. V. In Appendix B, it is shown how the implemented LBM can recover the non-dimensional continuity and momentum equations. In Appendix C, the Courant–Friedrichs–Lewy (CFL) condition is assessed in the introduced LBM by defining the requirement relations. Appendix D demonstrates how the non-dimensional equations are obtained from the presented LBM energy distribution relations.

II. PROBLEM STATEMENT

Internal natural convection of viscoplastic fluids in an enclosure with the width of W and the height of H is studied (see Fig. 1), using the exact model of Bingham. All lengths have been scaled with W. The aspect ratio of the enclosure is defined by A = H/W. An inclination of $ϕ$ against the horizontal line is applied to the cavity. The only applied external acceleration is gravity, and the Oberbeck–Boussinesq approximation is applied for density change. A volumetric uniform heating rate (Q_i) is induced through the enclosure, which drives the convection. Isotherm low temperature is imposed on the vertical sides, while the horizontal walls are adiabatic.

FIG. 1.

View large Download slide

Schematic of the studied problem.

A. Dimensional macroscopic governing equations and boundary conditions

The dimensional equations (continuity, momentum, and energy) are defined as

div v = 0,

(5)

\partial_{t} v + div [(v \otimes v) - \frac{T}{ρ}] = β (θ_{C} - θ) g (cos ϕ e_{y} + sin ϕ e_{x}),

(6)

\partial_{t} θ + div [v θ - α div θ] = \frac{1}{ρ c} [Q_{i} + σ : grad v],

(7)

where the parameters of ρ,

θ_{C}

⁠, α, and c are density, the cold temperature, thermal diffusivity, and thermal capacity; respectively. e_x and e_y are the unit vectors in x and y directions. The definitions of div and grad are the operators of divergence and gradient. The total stress tensor is defined by

T = - p 1 + η_{p} I_{1} + \sqrt{2} σ_{y} λ .

(8)

The dimensional boundary conditions are defined as

{v |}_{Γ_{i}} = 0, i = 1, 2, 3, 4,

(9a)

{θ |}_{Γ_{1}} = {θ |}_{Γ_{2}} = θ_{C}, {\partial_{y} θ |}_{Γ_{3}} = {\partial_{y} θ |}_{Γ_{4}} = 0 .

(9b)

B. Non-dimensional macroscopic parameters

The non-dimensional parameters with the asterisks are defined here as

x^{*} = \frac{x}{W}, U = \sqrt{g β Δ θ W}, v^{*} = \frac{v}{U},

(10a)

t^{*} = \frac{t U}{W}, p^{*} = \frac{p}{ρ U^{2}}, θ^{*} = \frac{(θ - θ_{C})}{Δ θ},

(10b)

I_{1}^{*} = \nabla v^{*} + \nabla v^{* T}, Δ θ = \frac{Q_{i} W^{2}}{k},

(10c)

T^{*} = \frac{T}{ρ U^{2}} or T^{*} = - p^{*} 1 + \sqrt{\frac{Pr}{RR}} σ^{*},

(10d)

where U is the buoyancy velocity scale and the non-dimensional extra stress tensor (⁠

σ^{*}

⁠), Bingham number (Bn), Prandtl number (Pr), and Rayleigh–Roberts number (RR) are defined as

Bn = \frac{σ_{y} W}{U η_{p}}, σ^{*} = I_{1}^{*} + \sqrt{2} Bn λ,

(11)

Pr = \frac{η_{p} c}{k}, RR = \frac{ρ^{2} c g β Q_{i} W^{5}}{η_{p} k^{2}},

(12)

u^{*}

and

v^{*}

are the non-dimensional velocities in x and y directions, respectively. The non-dimensional boundary conditions are defined as

{v^{*} |}_{Γ_{i}} = 0, i = 1, 2, 3, 4,

(13a)

{θ^{*} |}_{Γ_{1}} = {θ^{*} |}_{Γ_{2}} = 0, {\partial_{y} θ^{*} |}_{Γ_{3}} = {\partial_{y} θ^{*} |}_{Γ_{4}} = 0 .

(13b)

C. Non-dimensional governing equations

The following non-dimensional governing equations by using the defined non-dimensional relations are found as

div v^{*} = 0,

(14)

\partial_{t} v^{*} + div [(v^{*} \otimes v^{*}) - T^{*}] = θ^{*} (cos ϕ e_{y} + sin ϕ e_{x}),

(15)

\partial_{t} θ^{*} + div [v^{*} θ^{*} - Γ div θ^{*}] = Q^{*} + Φ,

(16)

Γ = Q^{*} = \frac{1}{\sqrt{Pr RR}} .

(17)

Φ

represents the non-dimensional viscous dissipation term, which is calculated for the studied cases as

Φ = Ec \sqrt{\frac{Pr}{RR}} σ^{*} : grad v^{*},

(18)

where Ec is the Eckert number and is defined by

Ec = \frac{U^{2}}{c Δ θ} \equiv \frac{g β W}{c} .

(19)

It is observed that the Eckert number is similar to the Gebhart number (⁠

Ge = \frac{g β H}{c}

⁠),¹¹ considering the applied velocity characteristic. It was reported¹¹ that the parameter of (⁠

\frac{g β}{c}

⁠) for various fluids is small (⁠

\approx 10^{- 5}

to

10^{- 9}

⁠). For example, in the earth-mantle convection, we have [

β = 2 \times 10^{- 5} (1 / K), c = 1.2 \times 10^{3} ({kg / m}^{3})

], which generates

\frac{g β}{c} \approx 1.67 \times 10^{- 8}

⁠. So, the viscous dissipation term only can be effective in extensive large-scale problems. Consequently, this term plays an effective role in the earth-mantle convection process due to the large thickness of the convection layer around 100 km (Ref. 1), which generates significant Eckert (Gebhart) number (e.g., Ec = Ge = 0.001). We also consider this parameter here to show the derivations in macroscopic dimensional and non-dimensional equations and demonstrate how to apply this term in the implemented lattice Boltzmann method. Here, we assume Ec = Ge = 0 in most parts of this study and compare the results between Ec = 0 and 0.001 in the section on the viscous dissipation effect.

The maximum temperature (⁠

θ_{max}^{*}

⁠) is found in the centerline (⁠

x_{max}^{*} = 0.5

⁠), and the Nusselt number is calculated by the maximum temperature as⁴

Nu = \frac{1}{8 θ_{max}^{*}} .

(20)

The pure conduction (RR

\to

0–10³) is obtained at

θ_{max}^{*} = \frac{1}{8} = 0.125

⁠, and therefore, the Nusselt number would be equal to one (Nu = 1).

Another parameter [called yield number (Y), which is the ratio of the plastic and the buoyancy forces]^7–10 can be introduced to diminish the non-dimensional parameters in the non-dimensional constitutive equation as

Y = \frac{σ_{y}}{ρ g β Δ θ W} = Bn \sqrt{\frac{Pr}{RR}} .

(21)

The maximum value of Y regardless of the applied Prandtl and Rayleigh numbers generates a full unyielded position in the cavity and the pure conduction (or Nu = 1) in external natural convection as mentioned in Refs. 8–10. So, this observation will be assessed for the internal convection here. The maximum Bingham number (⁠

{Bn}_{max}

⁠) and the maximum yield number (⁠

Y_{max}

⁠) are the values that the unyielded parts occupy the enclosure (⁠

ϒ \leq σ_{y}

⁠). It should be noted that, by the definition, the maximum yield number (⁠

Y_{max}

⁠) for all studied Rayleigh and Prandtl numbers is the same as we have

Y_{max} = {Bn}_{max} \sqrt{\frac{Pr}{RR}} .

(22)

This definition is similar to the observed pattern in previous studies^7–10 as the momentum equation is still similar to the studies and the difference is observed in the modified Rayleigh number.

The stream function (⁠

ψ

⁠) and the vorticity (⁠

ω

⁠) are given by

v^{*} = - \partial_{x^{*}} ψ, u^{*} = \partial_{y^{*}} ψ,

(23)

ω = - \nabla^{2} ψ = \partial_{x^{*}} v^{*} - \partial_{y^{*}} u^{*} .

(24)

III. THE APPLIED NUMERICAL MODEL

To solve the non-dimensional macroscopic equations, a non-dimensional lattice Boltzmann method (LBM) is implemented, following the proposed method of So et al.,²⁶ which we developed in various problems of viscoplastic fluids.^23–25 The dimensional equations and proofs of the applied methods are explained and discussed in detail in our previous study.^21–25 Here, the equations and relations of the LBM are presented, and, in Appendixes A–D, we show how the presented non-dimensional lattice Boltzmann equations can satisfy the derived non-dimensional macroscopic equations.

A. Recovering the non-dimensional continuity and momentum equations using the LBM

To recover the non-dimensional continuity and momentum equations (see Appendix B), a non-dimensional discrete particle distribution function (⁠

f_{α}

⁠) is written as^21–25

ε λ [\partial_{t} f_{α} + ξ_{α} \cdot grad f_{α} - F_{α}] = f_{α}^{e} - f_{α},

(25)

where

ξ_{α}

is the non-dimensional mesoscopic velocity for particles in the α direction, which is defined through a D2Q9 lattice as

ξ_{α} = χ {\begin{array}{l} (0, 0), & α = 0, \\ (cos φ_{α}, sin φ_{α}), & α = 1, 3, 5, 7, \\ \sqrt{2} (cos φ_{α}, sin φ_{α}), & α = 2, 4, 6, 8. \end{array}

(26)

The angles of different particles (⁠

φ_{α}

⁠) in the lattice are found by

φ_{α} = (α - 1) π / 4

⁠. The parameter of

χ

is selected to control stability and the Courant–Friedrichs–Lewy (CFL) (see Appendix C) by the time step (⁠

Δ t^{*}

⁠) and the lattice space (⁠

Δ x^{*}

⁠). The parameter ε represents the Knudsen number (Kn), which is a small value (⁠

Kn < 0.001

⁠) to satisfy the continuum mechanics and no-slip boundary conditions.

λ

is the relaxation time.

f_{α}^{e}

(the non-dimensional equilibrium particle distribution function) is found as^21–25

f_{α}^{e} = A_{α} + ξ_{α} \cdot B_{α} + (ξ_{α} \otimes ξ_{α}) : C_{α},

(27)

where

A_{α}, B_{α}

⁠, and

C_{α}

are found as

A_{0} = 1 + \frac{T_{x x}^{*} + T_{y y}^{*}}{χ^{2}} - \frac{| v^{*} |^{2}}{χ^{2}}, A_{α} = 0, α = 1, \dots, 8

(28)

and

B_{α} = \frac{v^{*}}{2 χ^{2}}, α = 1, 3, 5, 7; B_{α} = 0, α = 2, 4, 6, 8.

(29)

Next,

C_{α} = \frac{1}{2 χ^{4}} [\begin{matrix} u^{* 2} - T_{x x}^{*} & 0 \\ 0 & v^{* 2} - T_{y y}^{*} \end{matrix}], α = 1, 3, 5, 7

(30)

and

C_{α} = \frac{1}{8 χ^{4}} [\begin{matrix} 0 & u^{*} v^{*} - T_{x y}^{*} \\ u^{*} v^{*} - T_{y x}^{*} & 0 \end{matrix}], α = 2, 4, 6, 8 .

(31)

Different elements of the non-dimensional total stress are defined as

T_{x x}^{*} = - p^{*} + 2 \sqrt{\frac{Pr}{RR}} \partial_{x} u,

(32a)

T_{y y}^{*} = - p^{*} + 2 \sqrt{\frac{Pr}{RR}} \partial_{y} v,

(32b)

T_{x y}^{*} = \sqrt{\frac{Pr}{RR}} [(\partial_{x} v + \partial_{y} u) + \sqrt{2} Bn λ_{x y}],

(32c)

T_{y x}^{*} = \sqrt{\frac{Pr}{RR}} [(\partial_{x} v + \partial_{y} u) + \sqrt{2} Bn λ_{y x}] .

(32d)

The force term of

F_{α}

is calculated by

F_{α} = \frac{1}{2 χ^{2}} {\begin{array}{l} 0, & α = 0, 2, 4, 6, 8, \\ θ^{*} (cos ϕ e_{y} + sin ϕ e_{x}) \cdot ξ_{α}, & α = 1, 3, 5, 7. \end{array}

(33)

B. Recovering the non-dimensional energy equation using the LBM

To obtain the non-dimensional energy equation, a non-dimensional energy distribution function

g_{α}

is introduced as (see Appendix D)

ε λ [\partial_{t} g_{α} + ξ_{α} \cdot grad g_{α} - G_{α}] = g_{α}^{e} - g_{α} .

(34)

g_{α}^{e}

(the equilibrium energy distribution function) is given by^23–25

g_{α}^{e} = D_{α} + ξ_{α} \cdot E_{α},

(35)

where

D_{α}

is defined as

D_{α} = θ^{*}, α = 0; D_{α} = 0, α = 1, \dots, 8 .

(36)

The vectors

E_{α}

are expressed as

E_{α} = \frac{1}{2 χ^{2}} [v^{*} θ^{*} - Γ div θ^{*}] - Φ^{*}, α = 1, 3, 5, 7,

(37)

E_{α} = 0, α = 0, 2, 4, 6, 8 .

(38)

The parameter of

Φ^{*}

is calculated by

Φ^{*} = Ec \sqrt{\frac{Pr}{RR}} σ^{*} v^{*} .

(39)

Finally, the parameter

G_{α}

is obtained as

G_{α} = Q^{*}, α = 1, 3, 5, 7,

(40)

G_{α} = 0, α = 0, 2, 4, 6, 8.

(41)

C. The applied algorithm

The lattice Boltzmann equations (LBEs) of $f_{α}$ and $g_{α}$ are broken into two parts (streaming and collision) by the splitting method as

Streaming:

\partial_{t} f_{α} + ξ_{α} \cdot grad f_{α} - F_{α} = 0,

(42a)

\partial_{t} g_{α} + ξ_{α} \cdot grad g_{α} - G_{α} = 0 .

(42b)

Collision:

ε λ \partial_{t} f_{α} = f_{α}^{e} - f_{α},

(43a)

ε λ \partial_{t} g_{α} = g_{α}^{e} - g_{α} .

(43b)

Taking

ε λ \equiv Δ t^{*}

and using the Euler method, the collision parts are altered to

f_{α}^{(n + 1)} = f_{α}^{e n},

(44a)

g_{α}^{(n + 1)} = g_{α}^{e n},

(44b)

where

t^{* n} = n Δ t^{* (n + 1)}

⁠. The numerical process is as follows:

In this current LBM, the boundary conditions can be applied directly by the non-dimensional macroscopic values. The boundaries and initial values are utilized for the first non-dimensional macroscopic values of $(v^{* 0}, p^{* 0}, θ^{* 0})$ ⁠, and then, the initial amounts of the equilibrium distribution functions $f_{α}^{e 0}$ and $g_{α}^{e 0}$ are found.
The initial values of the distribution functions (⁠ $f_{α}^{0}$ and $g_{α}^{0}$ ⁠) are calculated through the collision part.
With $f_{α}^{0}$ and $g_{α}^{0}$ ⁠, the new distribution functions (⁠ $f_{α}^{1}$ and $g_{α}^{1}$ ⁠) are found in the streaming part.
Using $f_{α}^{1}$ and $g_{α}^{1}$ ⁠, the updated macroscopic parameters [ $(v^{* 1}, p^{* 1}, θ^{* 1})$ ] are determined as

$u^{*}^{1} \equiv \sum_{α = 0}^{8} f_{α}^{1} ξ_{α_{i}},$

(45a)

$v^{*}^{1} \equiv \sum_{α = 0}^{8} f_{α}^{1} ξ_{α_{j}},$

(45b)

$p^{* 1} \equiv \frac{1}{2} \sum_{i = 0}^{8} f_{α}^{1} {| ξ_{α} |}^{2} - {| v^{1} |}^{2} + σ_{x x}^{* 0} + σ_{y y}^{* 0},$

(45c)

$θ^{*}^{1} \equiv \sum_{α = 0}^{8} g_{α}^{1} .$

(45d)
With the renewed quantities [ $(v^{* 1}, p^{* 1}, θ^{* 1})$ ], the corresponding new equilibrium distribution functions become updated (⁠ $f_{α}^{e 1}$ and $g_{α}^{e 1}$ ⁠).
The updated distribution functions $f_{α}^{2}$ and $g_{α}^{2}$ are achieved in the collision part.
The convergence criteria (⁠ $ϵ \equiv Δ t^{2}$ ⁠) for different main parameters are written as

$| 1 - \sum_{α = 0}^{8} f_{α} | < ϵ,$

(46a)

$\sqrt{| v^{* n + 1} - v^{* n} |} \leq ϵ, \sqrt{| θ^{* n + 1} - θ^{* n} |} \leq ϵ .$

(46b)

In the present study, there is the opportunity to apply different boundary conditions directly, using the macroscopic variables similar to macroscopic numerical methods. So, we do not need any extra distribution functions for boundaries and the macroscopic boundaries just need to be written in the code straightforwardly. It is discussed in detail in Refs. 21–25.

IV. RESULTS

A. Validation and numerical characteristics

The present code was validated for external differentially heated natural convection in a cavity filled with viscoplastic fluids in our previous studies.^23–25 To evaluate our developed code for adding the external heat source, our results in the format of the non-dimensional temperature are compared with the study of Joshi et al.²⁷ The depicted plot in Fig. 2 showed a good agreement. To study the effects of different meshes and their required running time (T_r), the relative L₂-norm (or 2-norm) errors for the non-dimensional temperature (⁠

E_{θ}

⁠) are defined and calculated by

E_{θ} = \frac{| | Δ H | |}{| | H_{e} | |}

where we have (see Table I)

Δ H = \sqrt{\sum_{i = 1}^{k} {(θ^{*} - θ_{a_{i}}^{*})}^{2}}, H_{e} = \sqrt{\sum_{i = 1}^{k} θ_{a_{i}}^{* 2}},

(47)

where

θ_{a}^{*}

is the value of the analytical velocity in Joshi et al.²⁷ The final (developed) stage was defined at the non-dimensional time

t^{*}

= 30, and to observe the progress of the yielding parts through time progress, the yielded/unyielded zones are depicted at RR = 10⁵, Bn = 0.3, and Pr = 1 for A = 1 and

ϕ = 0 °

in different non-dimensional times (⁠

t^{*}

= 5, 15, 20, and 30) in Fig. 3. The white part is the yielded zone (fluid phase), and the black section represents the unyielded region (solid phase). Table II shows the required number of meshes and the running times for this case, using the relative errors of the Nusselt number (E_n), which is calculated as

E_{n} = | \frac{{Nu}^{n + 1} - {Nu}^{n}}{{Nu}^{n}} | .

(48)

FIG. 2.

Comparison of non-dimensional temperature ( θ ∗) at y ∗ = 0.5 with the study of Joshi et al.27 at RR = 103, Bn = 0, and Pr = 0.7 for A = 1 and ϕ = 0 °.

View large Download slide

Comparison of non-dimensional temperature (⁠ $θ^{*}$ ⁠) at $y^{*} = 0.5$ with the study of Joshi et al.²⁷ at RR = 10³, Bn = 0, and Pr = 0.7 for A = 1 and $ϕ = 0 °$ ⁠.

TABLE I.

Studying errors $(E_{θ})$ and running times $(T_{r})$ for different uniform meshes (⁠ $Δ x^{*} = Δ y^{*}$ ⁠) compared with Joshi et al.²⁷ at RR = 10³ and Pr = 0.7 for Bn = 0, A = 1, and $ϕ = 0 °$ ⁠.

$Δ x^{*}$	$E_{θ} (%)$	$T_{r} (s)$
1/25	12.2	120
1/50	2.94	283
1/75	0.51	480
1/100	0.28	692
1/150	0	1109

FIG. 3.

Comparison of yielded (white)/unyielded (black) parts in different non-dimensional times at RR = 105, Bn = 0.3, Pr = 1, A = 1, and ϕ = 0 ° for (a) t ∗ = 5, (b) t ∗ = 15, (c) t ∗ = 20, (d) t ∗ = 30.

View large Download slide

Comparison of yielded (white)/unyielded (black) parts in different non-dimensional times at RR = 10⁵, Bn = 0.3, Pr = 1, A = 1, and $ϕ = 0 °$ for (a) $t^{*} = 5$ ⁠, (b) $t^{*} = 15$ ⁠, (c) $t^{*} = 20$ ⁠, (d) $t^{*} = 30$ ⁠.

TABLE II.

Relative errors (E_n) of Nusselt number (Nu) in different meshes and running times (T_r) at RR = 10⁵, Bn = 0.4, and Pr = 1 for A = 1 and $ϕ = 0 °$ ⁠.

$Δ x^{*}$	1/50	1/100	1/150	1/200	1/250
E_n	$5.33 \times 10^{- 5}$	$2.97 \times 10^{- 6}$	$6.52 \times 10^{- 7}$	$1.02 \times 10^{- 8}$	$1.02 \times 10^{- 8}$
$T_{r} (s)$	290	720	1200	1950	3022

$Δ x^{*}$	1/50	1/100	1/150	1/200	1/250
E_n	$5.33 \times 10^{- 5}$	$2.97 \times 10^{- 6}$	$6.52 \times 10^{- 7}$	$1.02 \times 10^{- 8}$	$1.02 \times 10^{- 8}$
$T_{r} (s)$	290	720	1200	1950	3022

It confirmed that the grid size (⁠ $Δ x^{*}$ = $Δ y^{*}$ = 1/200) ensures a grid-independent solution as demonstrated. In all calculations, we set the time step of $Δ t^{*} = 10^{- 4}$ and the CPU of Intel(R)-Core(TM)i7-7700 is used.

B. The effect of Rayleigh–Roberts number

Figures 4–6 display the isotherms and yielded/unyielded zones for various Bingham and Rayleigh–Roberts numbers at Pr = 1, A = 1, and $ϕ = 0 °$ ⁠. It is seen the stagnant unyielded areas in the corners of the cavity are generated in small Bingham numbers in all studied RRs. Three main other zones are shaped in the middle of the cavity and in parallel to the vertical walls. The enhancement of Bingham number expands the mentioned primary unyielded parts. The isotherms illustrate the convection process weakens (it could be observed with a constant isotherm in the middle of the enclosure in various RRs) gradually as Bingham number rises. Figure 7 exhibits the patterns of streamlines (ψ), horizontal velocity (⁠ $u^{*}$ ⁠), vertical velocity (⁠ $v^{*}$ ⁠), and vorticity (ω) at Bn = 0.1, A = 1 and $ϕ = 0 °$ ⁠, and Pr = 1 in various RRs. The contours demonstrate the small values and close to zero values of velocities are generated in the corners and the middle of the cavity. In addition, the vertical velocities reveal the reason for the unyielded zones in parallel to the vertical sides, which are formed in the center of the two vertical velocities circulations in the left and right sides.

FIG. 4.

Yielded (white)/unyielded (black) sections (left) and isotherms (right) at RR = 104, Pr = 1, A = 1, and ϕ = 0 ° in different Bingham numbers (a) Bn = 0.05 (Y = 0.0005), (b) Bn = 0.1 (Y = 0.001), and (c) Bn = 0.25 (Y = 0.0025).

View large Download slide

Yielded (white)/unyielded (black) sections (left) and isotherms (right) at RR = 10⁴, Pr = 1, A = 1, and $ϕ = 0 °$ in different Bingham numbers (a) Bn = 0.05 (Y = 0.0005), (b) Bn = 0.1 (Y = 0.001), and (c) Bn = 0.25 (Y = 0.0025).

FIG. 5.

Yielded (white)/unyielded (black) sections (left) and isotherms (right) at RR = 105, Pr = 1, A = 1 and ϕ = 0 ° in different Bingham numbers (a) Bn = 0.1 (Y = 0.000 316), (b) Bn = 0.4 (Y = 0.001 26), and (c) Bn = 0.8 (Y = 0.002 53).

View large Download slide

Yielded (white)/unyielded (black) sections (left) and isotherms (right) at RR = 10⁵, Pr = 1, A = 1 and $ϕ = 0 °$ in different Bingham numbers (a) Bn = 0.1 (Y = 0.000 316), (b) Bn = 0.4 (Y = 0.001 26), and (c) Bn = 0.8 (Y = 0.002 53).

FIG. 6.

Yielded (white)/unyielded (black) sections (left) and isotherms (right) at RR = 106, Pr = 1, A = 1, and ϕ = 0 ° in different Bingham numbers (a) Bn = 0.4 (Y = 0.0004), (b) Bn = 1 (Y = 0.001), and (c) Bn = 2 (Y = 0.002).

View large Download slide

Yielded (white)/unyielded (black) sections (left) and isotherms (right) at RR = 10⁶, Pr = 1, A = 1, and $ϕ = 0 °$ in different Bingham numbers (a) Bn = 0.4 (Y = 0.0004), (b) Bn = 1 (Y = 0.001), and (c) Bn = 2 (Y = 0.002).

FIG. 7.

View large Download slide

Contours of streamlines (ψ), horizontal velocity (⁠ $u^{*}$ ⁠), vertical velocity (⁠ $v^{*}$ ⁠), and vorticity (ω) at Bn = 0.1, A = 1, $ϕ = 0 °$ ⁠, and Pr = 1 in different Rayleigh–Roberts numbers.

Figure 8 shows the Nusselt number drops steadily as Bingham number rises in different RRs. The maximum (or critical) Bingham number, which represents the full unyielded condition in the cavity and the Nusselt number is equal to one with the dominant conduction process, is enhanced by the increase in RR number. The maximum Bingham numbers are ${Bn}_{max}$ = 0.4, 1.3, and 4 for RR = $10^{3}, 10^{4}$ ⁠, and 10⁵, respectively, at Pr = 1, A = 1, and $ϕ = 0 °$ ⁠. It was observed that the slope of the decrease in the Nusselt number due to the rise of Bingham number is steep at $Bn < 0.5 {Bn}_{max}$ and around 90% of the dissipated Nusselt number is obtained in this part. Figure 9 reveals the maximum (or critical) yield number for all considered Rayleigh–Roberts numbers is the same (⁠ $Y_{max} = 0.004$ ⁠), which is a similar trend that was reported in external natural convection in previous studies.^8–10

FIG. 8.

Nusselt number at Pr = 1, A = 1, and ϕ = 0 ° in different Bingham (Bn) and Rayleigh–Roberts (RR) numbers (a) RR = 104, (b) RR = 105, and (c) RR = 106.

View large Download slide

Nusselt number at Pr = 1, A = 1, and $ϕ = 0 °$ in different Bingham (Bn) and Rayleigh–Roberts (RR) numbers (a) RR = 10⁴, (b) RR = 10⁵, and (c) RR = 10⁶.

FIG. 9.

View large Download slide

Nusselt number at Pr = 1, A = 1, and $ϕ = 0 °$ in different yield numbers (Y) and Rayleigh–Roberts (RR) numbers.

C. The effect of Prandtl number

Figure 10 shows how the unyielded zones enlarge as the Prandtl number increases from Pr = 0.1 to 100 in a fixed Bingham number of Bn = 0.1 at RR = 10⁵, A = 1, and $ϕ = 0 °$ ⁠. Figure 11 illustrates the influences of Bingham numbers on the non-dimensional temperature at $x^{*} = 0.5$ and vertical velocity at $y^{*} = 0.5$ for RR = 10⁴, Pr = 0.1, A = 1, and $ϕ = 0 °$ ⁠. It shows the convection process is weakened by the increase in the Bingham number where the curve forms of velocity and temperature in the middle of the cavity become linear steadily.

FIG. 10.

Yielded (white)/unyielded (black) sections at Bn = 0.1, RR = 105, A = 1, and ϕ = 0 ° in different Prandtl numbers.

View large Download slide

Yielded (white)/unyielded (black) sections at Bn = 0.1, RR = 10⁵, A = 1, and $ϕ = 0 °$ in different Prandtl numbers.

FIG. 11.

Non-dimensional vertical velocity ( v ∗) at y ∗ = 0.5 and non-dimensional temperature ( θ ∗) at x ∗ = 0.5 for RR = 104, Pr = 0.1, A = 1, and ϕ = 0 ° in different Bingham numbers.

View large Download slide

Non-dimensional vertical velocity (⁠ $v^{*}$ ⁠) at $y^{*} = 0.5$ and non-dimensional temperature (⁠ $θ^{*}$ ⁠) at $x^{*} = 0.5$ for RR = 10⁴, Pr = 0.1, A = 1, and $ϕ = 0 °$ in different Bingham numbers.

Figure 12 displays the Nusselt number in different Prandtl, Bingham, and yield numbers at RR = 10⁵ and A = 1, and $ϕ = 0 °$ ⁠. It indicates that the rise of Prandtl number enhances the Nusselt number in Newtonian fluid (Bn = 0), but, interestingly, the maximum Bingham number (⁠ ${Bn}_{max}$ ⁠) drops by the increase in Prandtl number. For Pr = 0.1, 1, 10, and 100, the critical Bingham number is ${Bn}_{max} = 4, 1.3, 0.4,$ and 0.15; respectively. This figure proves the maximum yield number is independent of Prandtl number where $Y_{max} = 0.004$ for various studied Prandtl numbers. It shows that the Nusselt number increases against the augmentation of Prandtl number in various studied yield numbers.

FIG. 12.

Nusselt number in different Prandtl (Pr), Bingham (Bn), and yield (Y) numbers at RR = 105, A = 1 and ϕ = 0 °.

View large Download slide

Nusselt number in different Prandtl (Pr), Bingham (Bn), and yield (Y) numbers at RR = 10⁵, A = 1 and $ϕ = 0 °$ ⁠.

D. The effect of Eckert (Gebhart) number

Figure 13 illustrates the effect of Eckert number between Ec = 0 and 0.001 on isotherms, streamlines, and yielded/unyielded zones at Pr = 1, A = 1, RR = 10⁵, and $ϕ = 0 °$ ⁠. The contour of isotherms in the left side represents Ec = 0, and the right one belongs to Ec = 0.001. The isotherm line of $θ^{*} = 0.11$ demonstrates the decrease in the convection process by the presence of the Eckert number as it starts moving from the centerline toward the sides. This phenomenon also is depicted in the streamlines. It is observed the unyielded parts expand slightly when the Eckert number rises from Ec = 0 (the green lines) to Ec = 0.001 (the red lines). To observe the influence of Eckert number on the convection process and its viscous dissipation process, the non-dimensional temperature in the middle of the cavity (⁠ $y^{*} = 0.5$ ⁠) for different Rayleigh–Roberts numbers at Pr = 1, A = 1, and $ϕ = 0 °$ is demonstrated in Fig. 14. It shows the maximum value of the dimensional temperature (θ_max) due to the viscous dissipation increases significantly, which generates lower Nusselt number.

FIG. 13.

Isotherms [Ec = 0 (left), Ec = 0.001 (right)], streamlines [Ec = 0 (solid), Ec = 0.001 (dash)], yielded/unyielded zones [Ec = 0 (green), Ec = 0.001 (red)] at Pr = 1, A = 1 RR = 105, and ϕ = 0 °.

View large Download slide

Isotherms [Ec = 0 (left), Ec = 0.001 (right)], streamlines [Ec = 0 (solid), Ec = 0.001 (dash)], yielded/unyielded zones [Ec = 0 (green), Ec = 0.001 (red)] at Pr = 1, A = 1 RR = 10⁵, and $ϕ = 0 °$ ⁠.

FIG. 14.

The non-dimensional temperature for Ec = 0 (red line) and 0.001 (green line) at Pr = 1, A = 1 and ϕ = 0 ° for (a) RR = 104, (b) RR = 105, and (c) RR = 106.

View large Download slide

The non-dimensional temperature for Ec = 0 (red line) and 0.001 (green line) at Pr = 1, A = 1 and $ϕ = 0^{°}$ for (a) RR = 10⁴, (b) RR = 10⁵, and (c) RR = 10⁶.

E. The effect of the enclosure aspect ratio

To see how the aspect ratio of the enclosure affects the convection process and unyielded parts, two cases of A = 0.5 and A = 2 are investigated to cover both scenarios of A <1 and A >1. Figures 15 and 16 depict the isotherms and yielded/unyielded zones in different Bingham numbers of Bn = 0.2 and 0.4 at Pr = 1, $ϕ = 0 °$ ⁠, and RR = 10⁵ for A = 0.5 and 2, respectively. At A = 0.5, the unyielded parts are formed in the center of the enclosure, but at A = 2, a pattern like A = 1 for unyielded parts is seen where three separate zones are produced, one in the center and two symmetric ones parallel to the side walls. In both aspect ratios, the stagnant unyielded parts are generated in the corners and in the middle of the horizontal walls. The isotherms show the convection process becomes weaker where the isotherms of $θ^{*} = 0.12$ and 0.11 move away from the top horizontal walls at A = 0.5 and 2. Figure 17 displays the Nusselt number at A = 0.5 drops compared to A = 1 generally in all considered Rayleigh–Roberts numbers. The Nusselt number decreases gradually as the Bingham number rises. However, this decline is not linear and similar to A = 1 and the most reduction in the Nusselt number is observed at $Bn < = 0.5 {Bn}_{max}$ ⁠. From the comparison between Figs. 17 and 8, it is clear the maximum (or critical) Bingham number falls due to the decrease in the aspect ratio from A = 1 to 0.5. The percentage of this drop in all RRs are nearly the same, and it is approximately $32.5 %$ ⁠. The yield number in Fig. 18 can approve this trend clearly where it shows the unique maximum yield number for all RRs are dropped by the same percentage of $32.5 %$ ⁠, compared to the value at A = 1, which is presented in Fig. 9. Figure 19 depicts the Nusselt number enhances in all Bingham numbers at A = 2 compared to A = 1 and the maximum Bingham numbers in RRs become nearly double. The unique maximum yield number in Fig. 20 demonstrates this trend where it rises by 205.5% (⁠ $Y_{max} = 8.27 \times 10^{- 3}$ ⁠) compared to the one at A = 1.

FIG. 15.

Isotherms and yielded (white)/unyielded (black) zones in different Bingham numbers at Pr = 1, A = 0.5, ϕ = 0 °, and RR = 105 for (a) Bn = 0.2 and (b) Bn = 0.4.

View large Download slide

Isotherms and yielded (white)/unyielded (black) zones in different Bingham numbers at Pr = 1, A = 0.5, $ϕ = 0 °$ ⁠, and RR = 10⁵ for (a) Bn = 0.2 and (b) Bn = 0.4.

FIG. 16.

Isotherms and yielded (white)/unyielded (black) zones in different Bingham numbers at Pr = 1, A = 2, ϕ = 0 °, and RR = 105 for (a) Bn = 0.2 and (b) Bn = 1.2.

View large Download slide

Isotherms and yielded (white)/unyielded (black) zones in different Bingham numbers at Pr = 1, A = 2, $ϕ = 0 °$ ⁠, and RR = 10⁵ for (a) Bn = 0.2 and (b) Bn = 1.2.

FIG. 17.

Nusselt number at Pr = 1, A = 0.5, and ϕ = 0 ° in different Bingham (Bn) and Rayleigh–Roberts (RR) numbers (a) RR = 104, (b) RR = 105, and (c) RR = 106.

View large Download slide

Nusselt number at Pr = 1, A = 0.5, and $ϕ = 0 °$ in different Bingham (Bn) and Rayleigh–Roberts (RR) numbers (a) RR = 10⁴, (b) RR = 10⁵, and (c) RR = 10⁶.

FIG. 18.

View large Download slide

Nusselt number at Pr = 1, A = 0.5, and $ϕ = 0 °$ in different yield numbers (Y) and Rayleigh–Roberts (RR) numbers.

FIG. 19.

Nusselt number at Pr = 1, A = 2, and ϕ = 0 ° in different Bingham (Bn) and Rayleigh–Roberts (RR) numbers (a) RR = 104, (b) RR = 105, and (c) RR = 106.

View large Download slide

Nusselt number at Pr = 1, A = 2, and $ϕ = 0 °$ in different Bingham (Bn) and Rayleigh–Roberts (RR) numbers (a) RR = 10⁴, (b) RR = 10⁵, and (c) RR = 10⁶.

FIG. 20.

View large Download slide

Nusselt number at Pr = 1, A = 2, and $ϕ = 0 °$ in different yield numbers (Y) and Rayleigh–Roberts (RR) numbers.

F. The effect of the enclosure inclined angle

To observe the effects of the inclined angle on isotherms and yielded/unyielded zones, they are evaluated in two angles of $ϕ = 40 °$ and $80 °$ at Pr = 1, A = 1, and RR = 10⁵ in Figs. 21 and 22. They illustrate that the symmetric shape of unyielded zones disappears with the inclined angle's rise since the symmetric shapes of velocities and streamlines are changed due to the inclination. The size and value of the secondary circulation of streamlines diminish, and the unyielded zone is formed in this part of the cavity. So, the left side of the enclosure is filled with solid material due to the increase in the Bingham number. However, the unyielded zone still grows in the right half of the enclosure and corners the same as the inclined angle of $ϕ = 0 °$ ⁠. Figures 21 and 22 demonstrate the increase in the Bingham number causes the isotherms to move to the side walls, which is a sign of the diminish in the convection process. Figures 23–26 show the Nusselt number drops as the inclined angle increases from $ϕ = 0 °$ to 80° for all studied RRs. Further, the maximum Bingham and yield numbers decrease steadily with the increase in the inclination angle. It was observed there is a unique single value of the maximum yield number for all studied Rayleigh–Roberts numbers for each angle where they are $Y_{max} = 3.2 \times 10^{- 3}$ and $1.6 \times 10^{- 3}$ for $ϕ = 40 °$ and $80 °$ ⁠, respectively. In other words, the maximum yield number declines by 20% and 60% for $ϕ = 40 °$ and 80° against the value at $ϕ = 0 °$ ⁠.

FIG. 21.

Isotherms and yielded (white)/unyielded (black) zones at Pr = 1, A = 1, ϕ = 40 °, and RR = 105 for (a) Bn = 0.3 and (b) Bn = 0.8.

View large Download slide

Isotherms and yielded (white)/unyielded (black) zones at Pr = 1, A = 1, $ϕ = 40 °$ ⁠, and RR = 10⁵ for (a) Bn = 0.3 and (b) Bn = 0.8.

FIG. 22.

Isotherms and yielded (white)/unyielded (black) zones at Pr = 1, A = 1, ϕ = 80 °, and RR = 105 for (a) Bn = 0.1 and (b) Bn = 0.2.

View large Download slide

Isotherms and yielded (white)/unyielded (black) zones at Pr = 1, A = 1, $ϕ = 80 °$ ⁠, and RR = 10⁵ for (a) Bn = 0.1 and (b) Bn = 0.2.

FIG. 23.

Nusselt number at Pr = 1, A = 1, and ϕ = 40 ° in different Bingham (Bn) and Rayleigh–Roberts (RR) numbers (a) RR = 104, (b) RR = 105, and (c) RR = 106.

View large Download slide

Nusselt number at Pr = 1, A = 1, and $ϕ = 40 °$ in different Bingham (Bn) and Rayleigh–Roberts (RR) numbers (a) RR = 10⁴, (b) RR = 10⁵, and (c) RR = 10⁶.

FIG. 24.

View large Download slide

Nusselt number at Pr = 1, A = 1, and $ϕ = 40 °$ in different yield numbers (Y) and Rayleigh–Roberts (RR) numbers.

FIG. 25.

Nusselt number at Pr = 1, A = 1, and ϕ = 80 ° in different Bingham (Bn) and Rayleigh–Roberts (RR) numbers (a) RR = 104, (b) RR = 105, and (c) RR = 106.

View large Download slide

Nusselt number at Pr = 1, A = 1, and $ϕ = 80 °$ in different Bingham (Bn) and Rayleigh–Roberts (RR) numbers (a) RR = 10⁴, (b) RR = 10⁵, and (c) RR = 10⁶.

FIG. 26.

View large Download slide

Nusselt number at Pr = 1, A = 1 and $ϕ = 80 °$ in different yield numbers (Y) and Rayleigh–Roberts (RR) numbers.

V. SUMMARY AND DISCUSSION

Natural convection driven by a uniform internal heat source in a cavity filled with viscoplastic fluids are studied motivated by its application in process of the earth's mantle convection.⁵ The constitutive equation is defined by an approach so-called the exact Bingham model, which recovers the Bingham model without any simplification. To solve the non-dimensional governing equations, a lattice Boltzmann method is developed based on the proposed method of So et al.²⁶ which we introduced and implemented in different non-Newtonian fluids, validated against previous studies into isothermal and thermal viscoplastic fluids.^21–25 The LBM code also was validated by previous studies into internal natural convection in a cavity and demonstrated a good agreement. We can summarize our findings as

At Pr = 1, A = 1, and $ϕ = 0 °$ ⁠, the maximum (or critical) Bingham number was obtained at ${Bn}_{max}$ = 0.4, 0.13, and 4 for RR = 10⁴, 10⁵, and 10⁶, respectively.
As it was reported in external natural convection^8–10 in cavities (which the yield number is constant in all studied Rayleigh numbers), it was observed the internal convection follows the same pattern and the maximum yield number is independent of $Rayleigh - Roberts$ number. The critical yield number for Pr = 1, A = 1, and $ϕ = 0 °$ was found at $Y_{max}$ = $4 \times 10^{- 3}$ ⁠.
At A = 1, the unyielded zones are expanded in the center and close to side walls with the increase in the yield number, while stagnant points are generated in the corners.
The convection process becomes weakened as the yield (or Bingham) number enhances, and this weakness is demonstrated through the isotherms and velocities in the middle of the cavity.
It was displayed by the Nusselt number that the highest part of the decrease in the convection process (around 90%) is achieved in the first half of the maximum Bingham and yield numbers (⁠ $Y_{max} / 2$ and ${Bn}_{max} / 2$ ⁠).
Similar to the observed trend in the external natural convection of viscoplastic fluids,^8–10 the maximum (or critical) yield number is independent to Prandtl number and it is $Y_{max}$ = $4 \times 10^{- 3}$ at A = 1, and $ϕ = 0 °$ ⁠.
The maximum Bingham number drops steadily as Prandtl number rises from Pr = 0.1 to 100. However, the Nusselt number in different yield numbers augments regularly due to the increase in Prandtl numbers.
The convection process diminishes with the addition of significant value of Eckert (Gebhart) number (Ec = 0.001) in the energy equation (which can only be observed in large amount of the length characteristic such as the convection layer in the earth's mantle convection). In addition, the unyielded zones expand slightly due to the rise of Eckert (Gebhart) number and the convection process drops.
The decrease in the aspect ratio from A = 1 to 0.5 causes the Nusselt number to decline generally. The maximum Bingham number and the maximum yield number rise as the aspect ratio enhances. The general trend was also mentioned in the external natural convection in previous studies.^9,10
It was demonstrated there is a unique number of the maximum yield number for any studied aspect ratios. It is $Y_{max}$ = $3.2 \times 10^{- 3}$ and $Y_{max}$ = $8.22 \times 10^{- 3}$ for A = 0.5 and A = 2, respectively. If we compare the $Y_{max}$ of the aspect ratios against the value at A = 1, it drops by 20% at $A = 0.5$ and increases by 205.5% at A = 2, accordingly.
The rise of the inclined angle from $ϕ = 0 °$ to 80° causes the maximum Bingham (and Yield) numbers to decrease and the convection process weakens. It was observed the maximum yield number for all studied Rayleigh–Roberts number is a single specific number in various inclined angles as reported in the external natural convection. For example, at A = 1, the maximum yield number is declined to $Y_{max}$ = $3.2 \times 10^{- 3}$ and $Y_{max}$ = $1.6 \times 10^{- 3}$ for $ϕ = 40 °$ and $ϕ = 80 °$ ⁠, respectively.
The symmetric shape of the unyielded zones is changed, and the locations of the unyielded parts in the cavity are altered due to the rise of the inclined angle counter-clockwise. These alterations become stronger as the Rayleigh–Roberts number and the inclined angle enhances. In other words, the strengthen of the buoyancy term in the x direction due to the inclined angle expands the unyielded zones rigorously in the left side of the cavity in the counter-clockwise rotation.

In future studies, the combination of internal and external natural convection of yield stress fluids for differentially heated walls and Rayleigh–Bénard convection will be studied. Further, the mixed convection due to the internal heat source for viscoplastic fluids will be evaluated and reported.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Gholamreza Kefayati: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

NOMENCLATURE

A: Aspect ratio of the enclosure (⁠ $A = H / W$ ⁠)
Bn: Bingham number (⁠ $Bn = \frac{σ_{y} W}{η_{p} U}$ ⁠)
c: Specific heat capacity $(J / kg K)$
$C_{O}$: Courant–Friedrichs–Lewy (CFL)
e: Unit vector
f: Non-dimensional density distribution function
$f^{e}$: Non-dimensional equilibrium density distribution function
F: Non-dimensional discrete body force term
g: Gravity acceleration vector $(m / s^{2})$
g: Non-dimensional energy distribution function
$g^{e}$: Non-dimensional equilibrium energy distribution function
G: Non-dimensional discrete source term in the energy distribution functions
H: Height of the enclosure (m)
$I_{1}$: First Rivlin–Ericksen tensor $(1 / s)$
k: Thermal conductivity (W/m K)
Kn: Knudsen number
L: Characteristic length (m)
Nu: Nusselt number
p: Pressure (Pa)
Pr: Prandtl number (⁠ $Pr = \frac{η_{p} c}{k}$ ⁠)
RR: Rayleigh–Roberts number (⁠ $RR = \frac{ρ^{2} c g β Q_{i} W^{5}}{η_{p} k^{2}}$ ⁠)
Q_i: Internal volume heat source (⁠ $\frac{W}{m^{3}}$ ⁠)
T: Total stress tensor (Pa)
t: Time (s)
x, y: Cartesian coordinates (m)
u: Velocity in x direction (m/s)
U: The buoyancy velocity scale (m/s)
v: Velocity in y direction (m/s)
W: Width of the enclosure (m)
Y: Yield number $Y = \frac{σ_{y}}{ρ g β Δ θ W}$

Greek letters

$θ$: Temperature (K)
$Δ θ$: Temperature difference (K)
$ζ$: Discrete non-dimensional lattice velocity
β: Thermal expansion coefficient $(1 / K)$
$ϕ$: Inclined angle
$σ$: Extra stress tensor (Pa)
$σ_{y}$: Yield stress (Pa)
$Δ x$: Lattice spacing in x direction (m)
$Δ y$: Lattice spacing in y direction (m)
$Δ t$: Time increment (s)
ρ: Density $(kg / m^{3})$
η: Dynamic viscosity (Pa s)
χ: Parameter for the stability control in LBM
ε: Small parameter represents Knudsen number
$\dot{γ}$: Shear rate $(1 / s)$
$λ$: Non-dimensional relaxation time
$λ$: Viscoplasticity constraint tensor
$ϒ$: Second invariant of the extra stress tensor (Pa)
$ω$: Vorticity
$ψ$: Stream function

Subscripts

i: Internal
x: x Direction
y: y Direction
1: First
θ: Temperature
p: Plastic
r: Running
C: Cold
max: Maximum
α: Direction α in a lattice

NOMENCLATURE

A: Aspect ratio of the enclosure (⁠ $A = H / W$ ⁠)
Bn: Bingham number (⁠ $Bn = \frac{σ_{y} W}{η_{p} U}$ ⁠)
c: Specific heat capacity $(J / kg K)$
$C_{O}$: Courant–Friedrichs–Lewy (CFL)
e: Unit vector
f: Non-dimensional density distribution function
$f^{e}$: Non-dimensional equilibrium density distribution function
F: Non-dimensional discrete body force term
g: Gravity acceleration vector $(m / s^{2})$
g: Non-dimensional energy distribution function
$g^{e}$: Non-dimensional equilibrium energy distribution function
G: Non-dimensional discrete source term in the energy distribution functions
H: Height of the enclosure (m)
$I_{1}$: First Rivlin–Ericksen tensor $(1 / s)$
k: Thermal conductivity (W/m K)
Kn: Knudsen number
L: Characteristic length (m)
Nu: Nusselt number
p: Pressure (Pa)
Pr: Prandtl number (⁠ $Pr = \frac{η_{p} c}{k}$ ⁠)
RR: Rayleigh–Roberts number (⁠ $RR = \frac{ρ^{2} c g β Q_{i} W^{5}}{η_{p} k^{2}}$ ⁠)
Q_i: Internal volume heat source (⁠ $\frac{W}{m^{3}}$ ⁠)
T: Total stress tensor (Pa)
t: Time (s)
x, y: Cartesian coordinates (m)
u: Velocity in x direction (m/s)
U: The buoyancy velocity scale (m/s)
v: Velocity in y direction (m/s)
W: Width of the enclosure (m)
Y: Yield number $Y = \frac{σ_{y}}{ρ g β Δ θ W}$

Greek letters

$θ$: Temperature (K)
$Δ θ$: Temperature difference (K)
$ζ$: Discrete non-dimensional lattice velocity
β: Thermal expansion coefficient $(1 / K)$
$ϕ$: Inclined angle
$σ$: Extra stress tensor (Pa)
$σ_{y}$: Yield stress (Pa)
$Δ x$: Lattice spacing in x direction (m)
$Δ y$: Lattice spacing in y direction (m)
$Δ t$: Time increment (s)
ρ: Density $(kg / m^{3})$
η: Dynamic viscosity (Pa s)
χ: Parameter for the stability control in LBM
ε: Small parameter represents Knudsen number
$\dot{γ}$: Shear rate $(1 / s)$
$λ$: Non-dimensional relaxation time
$λ$: Viscoplasticity constraint tensor
$ϒ$: Second invariant of the extra stress tensor (Pa)
$ω$: Vorticity
$ψ$: Stream function

Subscripts

i: Internal
x: x Direction
y: y Direction
1: First
θ: Temperature
p: Plastic
r: Running
C: Cold
max: Maximum
α: Direction α in a lattice

APPENDIX A: THE VISCOPLASTICITY CONSTRAINT TENSOR

The symmetric traceless tensor

λ

has the following conditions:

{\begin{matrix} 1 : λ = 0, \\ λ : λ < 1 \Rightarrow I_{1} = 0, (unyielded parts) \\ λ : λ = 1 \Rightarrow I_{1} \neq 0, (yielded parts) . \end{matrix}

(A1)

The above-presented conditions provide the required relation for the stress tensor of viscoplastic fluids, viz.,

| σ | \leq σ_{y}

at

I_{1} = 0,

and

σ_{y} < | σ |

for

I_{1} \neq 0 .

The unyielded part (or the rigid section) and the yielded section (or the liquid part) are defined through the tensor

λ

with satisfying Eq. and applying the following steps and process:

The introduced relations in the Eqs. (A1) are imposed over the entire flow domain (including both yielded and unyielded sections).
Finding the solutions of velocity $v$ and the viscoplasticity constraint tensor $λ$ to reveal the yielded/unyielded regions. There are no singularities due to $I_{1} / \dot{γ}$ as $I_{1} \to 0 .$
The two unknown vectors of $v$ and $λ$ should be defined, which requires a connection between these two vectors. It is proved that there is a relation between them using a simple projection operation (⁠ $P R_{M}$ ⁠)^6,23 as

$λ = P R_{M} (λ + r σ_{y} I_{1}), \forall r > 0,$

(A2)

where

$M = {Π | Π = {(Π_{i j})}_{1 \leq i, j \leq 2} \in {(L^{2} (Ω))}^{4}, | | Π | | \leq 1 a . e . on Ω}$

and

$P R_{M} : {(L^{2} (Ω))}^{4} \to M$

(A3)

is the projection operator defined so that $P R_{M} (Π) = Π,$ if $| | Π | | \leq 1,$ and $P_{M} (Π) = Π / | | Π | |$ otherwise. Note that in the context of Eq. (A2), the tensor $Π = λ + r σ_{y} I_{1}$ and it is symmetric. Further, the tensor $Π$ must be dimensionless for $λ$ is also dimensionless. $r > 0$ is a value, which can be defined by non-dimensional parameters based on our studied case and problem. The iterations will continue to reach the convergence point based on our desired accuracy. It should be noted that the boundary between the yielded and unyielded regions is shaped between $| | λ | | < 1$ and $| | λ | | = 1$ ⁠.

APPENDIX B: THE LBM FOR THE NON-DIMENSIONAL CONTINUITY AND MOMENTUM EQUATIONS

The non-dimensional equilibrium density distribution function (⁠

f_{α}^{e}

⁠) should satisfy the following constraints:

\sum_{α = 0}^{8} f_{α}^{e} = 1,

(B1a)

\sum_{α = 0}^{8} f_{α}^{e} ξ_{α} = v^{*},

(B1b)

\sum_{α = 0}^{8} f_{α}^{e} ξ_{α} \otimes ξ_{α} = M,

(B1c)

\sum_{α = 0}^{8} f_{α}^{(m)} = \sum_{α = 0}^{8} f_{α}^{(m)} ξ_{α} = 0, m \geq 1 .

(B1d)

Regarding the first constraint in the Eq. , we fix the density equal to one for

f^{e}

and these criteria control the incompressible condition and convergence in this method. In fact, we apply the following condition as

| 1 - \sum_{α = 0}^{8} f_{α}^{e} | < ϵ ϵ \equiv Δ t^{2} .

(B2)

Here, m is the power of the elements in the Chapman–Enskog expansion through

f_{α}

as

f_{α} = f_{α}^{e} + ε f_{α}^{(1)} + ε^{2} f_{α}^{(2)} + O (ε^{3}) .

(B3)

The matrix M and the vector of the force term are defined by

M = [\begin{matrix} u^{* 2} - T_{x x}^{*} & u^{*} v^{*} - T_{x y}^{*} \\ u^{*} v^{*} - T_{y x}^{*} & v^{* 2} - T_{y y}^{*} \end{matrix}] .

(B4)

For the force term, it is also needed to follow two other constraints as

\sum_{α = 0}^{8} F_{α} = 0, \sum_{α = 0}^{8} F_{α} ξ_{α} = F = θ^{*} (cos ϕ e_{y} + sin ϕ e_{x}) .

(B5)

Substituting and summing over the Chapman–Enskog expansion in , we have

\begin{matrix} \partial_{t} {\overset{↗}{(\sum_{α = 0}^{8} f_{α}^{e})}}^{1} + div {\overset{↗}{(\sum_{α = 0}^{8} f_{α}^{e} ξ_{α})}}^{v^{*}} = {\overset{↗}{\sum_{α = 0}^{8} F_{α}}}^{0} - λ^{- 1} {\overset{↗}{\sum_{α = 0}^{8} f_{α}^{(1)}}}^{0} + O (ε) . \end{matrix}

(B6)

It is reduced to the continuity equation as

0 + div v^{*} = 0 + O (ε) .

(B7)

Having the following condition as:

ξ_{α} \cdot (grad f_{α}^{e} \cdot ξ_{α}) \equiv div (f_{α}^{e} ξ_{α} \otimes ξ_{α}) .

(B8)

We multiply the Eq. by

ξ

and sum over it, considering the above relations, the resulting equation is written as

\begin{matrix} \partial_{t} {\overset{↗}{(\sum_{α = 0}^{8} f_{α}^{e} ξ_{α})}}^{v^{*}} + div {\overset{↗}{\sum_{α = 0}^{8} (f_{α}^{e} ξ_{α} \otimes ξ_{α})}}^{M} \\ = {\overset{↗}{\sum_{α = 0}^{8} F_{i} ξ_{α}}}^{F} - λ^{- 1} {\overset{↗}{\sum_{α = 0}^{8} f_{α}^{(1)} ξ_{α}}}^{0} + O (ε) . \end{matrix}

(B9)

It is written as

\partial_{t} v^{*} + div M = F + 0 + O (ε) .

(B10)

Substituting the tensor M generates the momentum equation as

\begin{matrix} \partial_{t} v^{*} + div [(v^{*} \otimes v^{*}) - T^{*}] = θ^{*} (cos ϕ e_{y} + sin ϕ e_{x}) + 0 + O (ε), \end{matrix}

(B11)

so it demonstrates the introduced momentum equation.

APPENDIX C: CFL CONDITIONS

Recalling the non-dimensional lattice Boltzmann equation as

ε λ [\partial_{t} f_{α} + ξ_{α} \cdot grad f_{α} - F_{α}] = f_{α}^{e} - f_{α} .

(C1)

The Chapman–Enskog expansion has the following relation:

f_{α} = f_{α}^{e} + ε f_{α}^{(1)} + ε^{2} f_{α}^{(2)} + O (ε^{3}) .

(C2)

Summing over directions, multiplying by

{| ξ_{α} |}^{2}

⁠, it is given by

\begin{matrix} \partial_{t} \sum_{α = 0}^{8} f_{α}^{e} {| ξ_{α} |}^{2} + div (\sum_{α = 0}^{8} f_{α}^{e} ξ_{α}) {| ξ_{α} |}^{2} \\ = {\overset{↗}{\sum_{α = 0}^{8} F_{α} {| ξ_{α} |}^{2}}}^{0} + λ^{- 1} {\overset{↗}{f_{α}^{(1)} {| ξ_{α} |}^{2}}}^{0} + O (ε) . \end{matrix}

(C3)

As we know,

f_{α}^{e}

in different directions are found as

f_{0}^{e} = \frac{1}{χ^{2}} [χ^{2} - | v^{*} |^{2} + T_{x x}^{*} + T_{y y}^{*}],

(C4)

f_{1}^{e} = \frac{1}{2 χ} u^{*} + 0.5 χ^{- 2} [u^{*}^{2} - T_{x x}^{*}],

(C5)

f_{2}^{e} = f_{6}^{e} = 0.125 χ^{- 4} [u^{*} v^{*} - T_{x y}^{*}],

(C6)

f_{3}^{e} = \frac{1}{2 χ} v^{*} + 0.5 χ^{- 2} [v^{*}^{2} - T_{y y}^{*}],

(C7)

f_{4}^{e} = f_{8}^{e} = - 0.125 χ^{- 4} [u^{*} v^{*} - T_{x y}^{*}],

(C8)

f_{5}^{e} = - \frac{1}{2 χ} u^{*} + 0.5 χ^{- 2} [u^{*}^{2} - T_{x x}^{*}],

(C9)

f_{7}^{e} = - \frac{1}{2 χ} v^{*} + 0.5 χ^{- 2} [v^{*}^{2} - T_{y y}^{*}] .

(C10)

As a result, the first term in the Eq. is obtained by

\begin{matrix} \sum_{α = 0}^{8} f_{α}^{e} {| ξ_{α} |}^{2} = (u^{*}^{2} + v^{*}^{2}) - T_{x x}^{*} - T_{y y}^{*} = | v^{*} |^{2} - T_{x x}^{*} - T_{y y}^{*} \end{matrix} .

(C11)

From , the second term of is found by

(\sum_{α = 0}^{8} f_{α}^{e} ξ_{α}) {| ξ_{α} |}^{2} = v^{*} χ^{2} .

(C12)

is written as

\partial_{t} [{| v^{*} |}^{2} - T_{x x}^{*} - T_{y y}^{*}] + χ^{2} (div v^{*}) = O (ε) .

(C13)

The Courant–Friedrichs–Lewy (CFL) condition is determined by

C_{O} = \frac{u^{*} Δ t^{*}}{Δ x^{*}} + \frac{v^{*} Δ t^{*}}{Δ y^{*}} \leq 1.

(C14)

is altered to

[{| v^{*} |}^{2} - T_{x x}^{*} - T_{y y}^{*}] + χ^{2} C_{O} = O (ε) .

(C15)

Thus,

χ = \frac{1}{\sqrt{C_{O}}} \sqrt{| T_{x x}^{*} + T_{y y}^{*} - {| v^{*} |}^{2} |} .

(C16)

APPENDIX D: THE LBM FOR THE NON-DIMENSIONAL ENERGY EQUATION

g_{α}^{e}

should hold the following conditions:

\sum_{α = 0}^{8} g_{α}^{(m)} = 0, m \geq 1,

(D1a)

\sum_{α = 0}^{8} g_{α}^{e} = θ^{*},

(D1b)

\sum_{α = 0}^{8} g_{α}^{e} ξ_{α} = v^{*} θ^{*} - Γ div θ^{*} - Φ^{*},

(D1c)

\sum_{α = 0}^{8} G_{α} = Q^{*} .

(D1d)

The Chapman–Enskog expression for

g_{α}

gives

g_{α} = g_{α}^{e} + ε g_{α}^{(1)} + ε^{2} g_{α}^{(2)} + O (ε^{3}) .

(D2)

Then, the above relation for

g_{α}

is put in , summing over all directions as

\begin{matrix} \partial_{t} {\overset{↗}{(\sum_{α = 0}^{8} g_{α}^{e})}}^{θ^{*}} + div {\overset{↗}{(\sum_{α = 0}^{8} g_{α}^{e} ξ_{α})}}^{(v^{*} θ^{*} - Γ div θ^{*} - Φ^{*})} \\ = {\overset{↗}{\sum_{α = 0}^{8} G_{α}}}^{Q^{*}} - \frac{1}{λ} {\overset{↗}{\sum_{α = 0}^{8} g_{α}^{(1)}}}^{0} + O (ε) . \end{matrix}

(D3)

It gives

\partial_{t} (θ^{*}) + div (v^{*} θ^{*} - Γ div θ^{*} - Φ^{*}) = Q^{*} + O (ε) .

(D4)

We have

div (θ^{*} v^{*}) = θ^{*} \cdot {\overset{↗}{(div v^{*})}}^{0} + v^{*} \cdot (grad θ^{*}) \equiv v^{*} \cdot grad θ^{*}

(D5)

and

div (Φ^{*}) \equiv div (σ^{*} v^{*}) = tr (σ^{*} grad v^{*}),

(D6)

so the above equation can be written as

div (Φ^{*}) \equiv σ^{*} : grad v^{*} \equiv Φ .

(D7)

As a result, the Eq. is changed to the presented non-dimensional energy equation as

\partial_{t} θ^{*} + div [v^{*} θ^{*} - Γ div θ^{*}] - Φ = Q^{*} .

(D8)

REFERENCES

1.

D. P.

McKenzie

,

J. M.

Roberts

, and

N. O.

Weiss

, “

Numerical models of convection in the earth's mantle

,”

Tectonophysics

19

,

89

–

103

(

1973

).

https://doi.org/10.1016/0040-1951(73)90034-6

Google Scholar

Crossref

2.

P. H.

Roberts

, “

Convection in horizontal layers with internal heat generation

. Theory,”

J. Fluid Mech.

30

,

33

–

49

(

1967

).

https://doi.org/10.1017/S0022112067001284

Google Scholar

Crossref

3.

D. J.

Tritton

and

M. N.

Zarraga

, “

Convection in horizontal layers with internal heat generation.

Experiments,”

J. Fluid Mech.

30

,

21

–

31

(

1967

).

https://doi.org/10.1017/S0022112067001272

Google Scholar

Crossref

4.

D.

Goluskin

,

Internally Heated Convection and Rayleigh–Bénard Convection

(

Springer International Publishing

,

2016

).

Google Scholar

Crossref

5.

A.

Rozel

,

G. J.

Golabek

,

R.

Näf

, and

P. J.

Tackley

, “

Formation of ridges in a stable lithosphere in mantle convection models with a viscoplastic rheology

,”

Geophys. Res. Lett.

42

,

4770

–

4777

, https://doi.org/10.1002/2015GL063483 (

2015

).

Google Scholar

Crossref

PubMed

6.

R.

Huilgol

and

G.

Kefayati

, “

Natural convection problem in a Bingham fluid using the operator–splitting method

,”

J. Non-Newtonian Fluid Mech.

220

,

22

–

32

(

2015

).

https://doi.org/10.1016/j.jnnfm.2014.06.005

Google Scholar

Crossref

7.

J.

Zhangi

,

D.

Vola

, and

I. A.

Frigaard

, “

Yield stress effects on Rayleigh–Bénard convection

,”

J. Fluid Mech.

566

,

389

–

419

(

2006

).

https://doi.org/10.1017/S002211200600200X

Google Scholar

Crossref

8.

A.

Vikhansky

, “

Thermal convection of a viscoplastic liquid with high Rayleigh and Bingham numbers

,”

Phys. Fluids

21

,

103103

(

2009

).

https://doi.org/10.1063/1.3256166

Google Scholar

Crossref

9.

A.

Vikhansky

, “

On the onset of natural convection of Bingham liquid in rectangular enclosures

,”

J. Non-Newtonian Fluid Mech.

165

,

1713

–

1716

(

2010

).

https://doi.org/10.1016/j.jnnfm.2010.09.003

Google Scholar

Crossref

10.

I.

Karimfazli

,

I.

Frigaard

, and

A.

Wachs

, “

A novel heat transfer switch using the yield stress

,”

J. Fluid Mech.

783

,

526

–

566

(

2015

).

https://doi.org/10.1017/jfm.2015.511

Google Scholar

Crossref

11.

O.

Turan

,

R. J.

Pool

, and

N.

Chakraborty

, “

Influences of boundary conditions on laminar natural convection of Bingham fluids in rectangular enclosures with differentially heated side walls

,”

Heat Transfer Eng.

35

,

822

–

849

(

2014

).

https://doi.org/10.1080/01457632.2014.852870

Google Scholar

Crossref

12.

S.

Yigit

and

N.

Chakraborty

, “

Laminar natural convection of Bingham fluids in square cross-sectioned cylindrical annular cavity with differentially heated vertical walls subjected to constant heat fluxes

,”

Heat Transfer Eng.

38

,

1171

–

1188

(

2017

).

https://doi.org/10.1080/01457632.2016.1239937

Google Scholar

Crossref

13.

H.

Sajjadi

,

A. A.

Delouei

,

M.

Atashafrooz

, and

M.

Sheikholeslami

, “

Double MRT lattice Boltzmann simulation of 3-D MHD natural convection in a cubic cavity with sinusoidal temperature distribution utilizing nanofluid

,”

Int. J. Heat Mass Transfer

126

,

489

–

503

(

2018

).

https://doi.org/10.1016/j.ijheatmasstransfer.2018.05.064

Google Scholar

Crossref

14.

H.

Sajjadi

,

H.

Mohammadifar

, and

A. A.

Delouei

, “

Investigation of the effect of the internal heating system position on heat transfer rate utilizing Cu/water nanofluid

,”

J. Therm. Anal. Calorim.

139

,

2035

–

2054

(

2020

).

https://doi.org/10.1007/s10973-019-08611-z

Google Scholar

Crossref

15.

H.

Sajjadi

,

A. A.

Delouei

,

R.

Mohebbi

,

M.

Izadi

, and

S.

Succi

, “

Natural convection heat transfer in a porous cavity with sinusoidal temperature distribution using Cu/water nanofluid: Double MRT lattice Boltzmann method

,”

Commun. Comput. Phys.

29

,

292

–

318

(

2021

).

https://doi.org/10.4208/cicp.OA-2020-0001

Google Scholar

Crossref

16.

G.

Kefayati

,

A.

Tolooiyan

, and

A. P.

Dyson

, “

A finite difference lattice Boltzmann method for modelling dam break debris flows

,”

Phys. Fluids

35

,

013102

(

2023

).

https://doi.org/10.1063/5.0130947

Google Scholar

Crossref

17.

G.

Kefayati

, “

A macroscopic and mesoscopic model of Newtonian and non-Newtonian nanofluids with a two-energy equation method

,”

Phys. Fluids

34

,

112005

(

2022

).

https://doi.org/10.1063/5.0124292

Google Scholar

Crossref

18.

G.

Kefayati

, “

A two-and three-dimensional mesoscopic method for an updated non-homogeneous model of Newtonian and non-Newtonian nanofluids

,”

Phys. Fluids

34

,

032003

(

2022

).

https://doi.org/10.1063/5.0084694

Google Scholar

Crossref

19.

G.

Kefayati

and

A. P.

Bassom

, “

A lattice Boltzmann method for single-and two-phase models of nanofluids: Newtonian and non-Newtonian nanofluids

,”

Phys. Fluids

33

,

102008

(

2021

).

https://doi.org/10.1063/5.0067744

Google Scholar

Crossref

20.

G.

Kefayati

,

A.

Tolooiyan

,

A. P.

Bassom

, and

K.

Vafai

, “

A mesoscopic model for thermal–solutal problems of power-law fluids through porous media

,”

Phys. Fluids

33

,

033114

(

2021

).

https://doi.org/10.1063/5.0042526

Google Scholar

Crossref

21.

R. R.

Huilgol

and

G. H. R.

Kefayati

, “

From mesoscopic models to continuum mechanics: Newtonian and non-Newtonian fluids

,”

J. Non-Newtonian Fluid Mech.

233

,

146

–

154

(

2016

).

https://doi.org/10.1016/j.jnnfm.2016.03.002

Google Scholar

Crossref

22.

R. R.

Huilgol

and

G. H. R.

Kefayati

, “

A particle distribution function approach to the equations of continuum mechanics in Cartesian, cylindrical and spherical coordinates: Newtonian and non-Newtonian fluids

,”

J. Non-Newtonian Fluid Mech.

251

,

119

(

2018

).

https://doi.org/10.1016/j.jnnfm.2017.10.004

Google Scholar

Crossref

23.

G. H. R.

Kefayati

and

R. R.

Huilgol

, “

Lattice Boltzmann Method for simulation of mixed convection of a Bingham fluid in a lid-driven cavity

,”

Int. J. Heat Mass Transfer

103

,

725

–

743

(

2016

).

https://doi.org/10.1016/j.ijheatmasstransfer.2016.07.102

Google Scholar

Crossref

24.

G. H. R.

Kefayati

, “

Lattice Boltzmann simulation of double diffusive natural convection of viscoplastic fluids in a porous cavity

,”

Phys. Fluids

31

,

013105

(

2019

).

https://doi.org/10.1063/1.5074089

Google Scholar

Crossref

25.

G. H. R.

Kefayati

, “

An immersed boundary-lattice Boltzmann method for thermal and thermo-solutal problems of Newtonian and non-Newtonian fluids

,”

Phys. Fluids

32

,

073103

(

2020

).

https://doi.org/10.1063/5.0013977

Google Scholar

Crossref

26.

R. M. C.

So

,

R. C. K.

Leung

,

E. W. S.

Kam

, and

S. C.

Fu

, “

Progress in the development of a new lattice Boltzmann method

,”

Comput. Fluids

190

,

440

(

2019

).

https://doi.org/10.1016/j.compfluid.2019.04.009

Google Scholar

Crossref

27.

M. V.

Joshi

,

U. N.

Gaitonde

, and

S. K.

Mitra

, “

Analytical study of natural convection in a cavity with volumetric heat generation

,”

J. Heat Transfer

128

,

176

–

182

(

2006

).

https://doi.org/10.1115/1.2137761

Google Scholar

Crossref

Published open access through an agreement with University of Tasmania

Internally heated convection of viscoplastic fluids in enclosures using a lattice Boltzmann method

I. INTRODUCTION

II. PROBLEM STATEMENT

A. Dimensional macroscopic governing equations and boundary conditions

B. Non-dimensional macroscopic parameters

C. Non-dimensional governing equations

III. THE APPLIED NUMERICAL MODEL

A. Recovering the non-dimensional continuity and momentum equations using the LBM

B. Recovering the non-dimensional energy equation using the LBM

C. The applied algorithm

IV. RESULTS

A. Validation and numerical characteristics

B. The effect of Rayleigh–Roberts number

C. The effect of Prandtl number

D. The effect of Eckert (Gebhart) number

E. The effect of the enclosure aspect ratio

F. The effect of the enclosure inclined angle

V. SUMMARY AND DISCUSSION

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

NOMENCLATURE

Greek letters

Subscripts

NOMENCLATURE

Greek letters

Subscripts

APPENDIX A: THE VISCOPLASTICITY CONSTRAINT TENSOR

APPENDIX B: THE LBM FOR THE NON-DIMENSIONAL CONTINUITY AND MOMENTUM EQUATIONS

APPENDIX C: CFL CONDITIONS

APPENDIX D: THE LBM FOR THE NON-DIMENSIONAL ENERGY EQUATION

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Internally heated convection of viscoplastic fluids in enclosures using a lattice Boltzmann method

I. INTRODUCTION

II. PROBLEM STATEMENT

A. Dimensional macroscopic governing equations and boundary conditions

B. Non-dimensional macroscopic parameters

C. Non-dimensional governing equations

III. THE APPLIED NUMERICAL MODEL

A. Recovering the non-dimensional continuity and momentum equations using the LBM

B. Recovering the non-dimensional energy equation using the LBM

C. The applied algorithm

IV. RESULTS

A. Validation and numerical characteristics

B. The effect of Rayleigh–Roberts number

C. The effect of Prandtl number

D. The effect of Eckert (Gebhart) number

E. The effect of the enclosure aspect ratio

F. The effect of the enclosure inclined angle

V. SUMMARY AND DISCUSSION

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

NOMENCLATURE

Greek letters

Subscripts

NOMENCLATURE

Greek letters

Subscripts

APPENDIX A: THE VISCOPLASTICITY CONSTRAINT TENSOR

APPENDIX B: THE LBM FOR THE NON-DIMENSIONAL CONTINUITY AND MOMENTUM EQUATIONS

APPENDIX C: CFL CONDITIONS

APPENDIX D: THE LBM FOR THE NON-DIMENSIONAL ENERGY EQUATION

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only