Thermal Lattice Boltzmann Simulation of Entropy Generation within a Square Enclosure for Sensible and Latent Heat Transfers

This paper deals with the numerical simulation of heat transfer and entropy generation in a 2D square enclosure for convective melting. A thermal lattice Boltzmann method (TLBM) is used to handle the study, which has been conducted for Prandtl numbers from 0.02 to 70 at Rayleigh numbers of 10 4 and 10 5. The results are presented in terms of the total entropy generation, average Bejan number and average Nusselt number. Within the range considered for the key parameters, the entropy generation is found to be controlled by the heat transfer loss for low Prandtl numbers. However, for the large Prandtl numbers, its variation is dominated by shearing losses. Moreover, the presence of the latent heat state decreases the overall thermodynamic losses while increasing the quantity of heat transferred.


Introduction
Convective flows have been widely investigated for various applications in engineering, such as heat exchangers, drying processes, solar collectors, nanofluids, etc.With the recent development and interest in nano-fields, studies on convection were extended to include electrohydrodynamic and magnetohydrodynamic convection flows [1,2], besides natural convection, which is widely treated in the literature.Despite this, natural convection problems [3,4] are still of considerable interest in many applications, especially in heat transfer systems.
In natural convection problems, hydrodynamic and thermal mechanisms are coupled and strongly affected by the fluid thermo-physical properties and the temperature gradients present in the system.Thermal gradients, shearing, phase changes, chemical reactions, etc., are irreversible processes that are frequently encountered in thermodynamic systems.They are generally provoked by losses emanating from two main causes: the transmission of heat through a limited temperature difference and the pressure drop caused by friction [5].The sum of all of these losses induces what is called the entropy generation within the system.Therefore, design criteria and optimum efficiency for thermodynamic systems can be obtained by analyzing the entropy production.Therefore, this production should be better understood [6].Indeed, the quest for deeper understanding of the entropy generation in enclosures has been underway for a long time and gave rise to fruitful research.However, up to date, studies involving entropy generation during phase changes remain sparse.All of this motivates this paper to perform the current research.Concerning entropy generation in natural convection, one can quote the works of Magherbi et al. [7], Abu-Hijleh et al. [8], Mahmud and Island [9] and Ellahi et al. [10], among others.
The main aim of the present work is to study the entropy generation during melting with natural convection in a cavity via a numerical simulation.Here, we intend to analyze variations of the entropy generation in relation to the relevant parameters of the study, namely the Rayleigh number, the Prandtl number and the irreversibility's coefficient.With this parametric study, we seek to deepen our understanding of these mechanisms to better assess the performance of the considered system.The objective here is to depict the effect of the presence of a phase change on the variation of entropy generation.Note that any achieved decrease in the total entropy production would give a further advantage of the use of latent heat over the sensible heat in mechanical systems.Furthermore, given that the analysis of the entropy production was used as a thermodynamic optimization tool, it seems advisable, hence, to give a new perspective based on the second law of thermodynamics to the heat transfer.
To perform simulations of heat transfer with convection and phase change, we use the single relaxation lattice Boltzmann model (SRT-LBM).This approach is a discrete particle-based method that numerically solves the Boltzmann equation as opposed to conventional methods that are based on the Navier-Stokes (NS) equations.Its principle is to assign a distribution function (DF) that provides the probability to find a pseudo-particle at a position with a given velocity [11].However, its extension to flows with heat transfer is not straightforward due to the numerical instabilities engendered.In the thermal lattice Boltzmann method (TLBM), a separate DF is used to solve for the temperature [12,13].In other words, two sets of distribution functions are defined, one for the velocity field and the other for the temperature.Thereby, such an approach can easily handle arbitrary Prandtl numbers.This double distribution function (DDF) approach has been successfully used to solve thermal problems in two dimensions [3,14].Nevertheless, it only applies when the fluid density depends weakly on the temperature.This is the approach we have adopted here to conduct this work.However, to simulate the phase change process, the partial bounce-back approach was employed, with a modification of the velocity field to mimic the solid, liquid and mushy zones.This approach is considered an enhancement of the existing models based on the full bounce-back approach [15,16].
The layout of this paper is as follows.The problem description and mathematical formulation are introduced in Section 2. This is supplemented by the definition of local and total entropy generation and the Bejan number.Subsequently, the numerical model is outlined in Section 3. In this section, the TLBM method and related issues are briefly introduced and completed by the relevant equations and boundary conditions.Afterwards, in Section 4, a validation of the adopted numerical approach via available results has been performed.Section 5 presents and comments on our investigations for a range of relevant parameters, such as Rayleigh, Prandtl and the irreversibility's distribution.Finally, some conclusions are drawn in Section 6.

Problem Description and Formulation
Figure 1a depicts the cavity and coordinate system along with the boundary conditions considered in this study.It concerns a 2D square cavity differentially heated with a width L and a height H.In the present study, we set L = H.The confined fluid in the cavity is assumed incompressible and Newtonian, and the flow occurring inside is laminar.The Cartesian coordinate system is labeled (x, y) of which the x-axis is horizontal and the y-axis upwardly directed in the direction opposite gravity.During the numerical simulation process, the lower and upper walls are assumed perfectly insulated (adiabatic), and the two vertical walls are isothermal and held at uniform temperatures h T and c T Recall that the temperature difference between the vertical sides produces a temperature gradient in the cavity and a consequent density difference that induces a fluid motion, that is the convection.In addition, all thermo-physical properties of the fluid are assumed constant, except the density, which varies in so far as the Boussinesq's approximation is respected.When dealing with convection melting, Figure 1b shows the solid-liquid phase change problem, also called the melting by convection problem.The left vertical side is kept at h T with h m T T > , m T being the melting temperature.As for the right side, it is kept at ( ) . Meanwhile, the top and bottom sides remain adiabatic.Initially, the medium is in a solid state at m T .

Mathematical Model
To solve the entropy generation problem due to natural convection under the conjectures stated above, the governing equations, in the transient state, are the following: ( ) where ∇ is the nabla operator, / D DT is the material derivative, ( ) is the velocity vector, p is the pressure, T is the temperature, 0 T is the reference temperature ( ( ) / 2) , ρ is the fluid density, ν is the kinematic viscosity, β is the thermal expansion coefficient, α is the thermal diffusivity, g  is the gravity downward and ε ( ) is the latent heat term (source or sink), ε being the liquid fraction.Note that this quantity is negative for melting, positive for solidification and zero in the absence of a phase change.
To cast the above equations in a dimensional form, we employ the following dimensionless variables: Equations ( 1)-( 3), which govern the convective melting problem, then take the following dimensionless form: where the symbol (  ) indicates that the quantity concerned is dimensionless and Y   is the dimensionless vertical direction.Therefore, this problem involves the Prandtl number, the Rayleigh number and the Stefan number, defined by Pr These equations are subjected to the following boundary (BCs) and initial conditions (ICs): The non-slip boundary condition for the hydrodynamic problem is applied at all cavity walls.For the thermal problem, the hot wall was kept with 0.5 T =  and the cold wall with 0.5 = − T  . As for horizontal sides, the Neumann boundary condition (zero heat flux) is applied, i.e., / 0.0 In the cavity, the initial conditions (ICs) are as follows: zero for the velocity field and 0.5 for the temperature field.

Entropy Generation
It should be noted that the entropy production in a system is a key thermodynamic parameter that can improve the performance of such a system.In other words, the determination of this parameter is used to analyze the total losses and, hence, to propose an optimization procedure.In the natural convection process, the entropy generation is associated with the heat transfer and with the fluid friction.According to Bejan [17], the local entropy generation ( lo S ) can be expressed as: ( ) where k is the thermal conductivity of the fluid, μ is the dynamic viscosity and ij S is the rate of the deformation tensor.In this expression, the first term (right-hand) is due to the heat transfer ( , while the second term is due to the viscous effects of the fluid ( , l d S ).A dimensionless form of their expressions can be obtained as follows: with: ( ) S  can be expressed as: Here, 0 T is ( ) / 2 Therefore, Be 0.5 >> indicates the dominance of the heat transfer irreversibility; Be 0.5 << reflects the dominance of the fluid viscous effects; and Be 0.5 = means that the two phenomena do not compete.
The dimensionless total entropy generation ( T S  ) is obtained via the integration of the local entropy generation ( lo S  ) throughout the entire domain, namely: For the sake of clarity, henceforth, we will omit any symbol indicating dimensionless quantities.

Thermal Lattice Boltzmann Equations
The TLBM consists of simulating the statistical behavior of a set of particles on a lattice with finite velocities.It stems from the discrete Boltzmann equation and allows providing macroscopic fluid properties, such as density, velocity, pressure, etc., through weighted averages, or moments, of the particle distribution for all discrete lattice velocities.The SRT-LBM (also called the Bhatnagar-Gross-Krook (BGK) model) for incompressible thermal flows builds on two distribution functions (DFs), i f and i g , and their corresponding evolution equations [18,19].As such, these are given here (Equations ( 13) and ( 14)) in the form of two steps, which are collision and streaming (advection) processes: g x e t g x t g x t g x t S where i e is the microscopic particle velocity in the i -direction, f τ and h τ are the dimensionless relaxation times and eq i f and eq i g are local equilibrium distributions functions that can be computed from: and: 2 . 1 (16) where: are, respectively, the weight coefficient and the velocity vector of the D2Q9 model; ( ) is the macroscopic velocity, with u and v representing velocities in the x -and y -directions, respectively.
Note that the relaxation times f τ and h τ can be determined via . By this treatment, there is no need to add a force term to the collision operator.On the other hand, in the g distribution function, the source term is treated as per the method proposed by Luo [21].Hence, the resulting force in the LBM frame will be: ∂ε ∂ being the source (or sink) term that handles the phase change.

h c T h h h h h h L h h h h
where h is the local enthalpy defined, at a time step j, by: ( ) It is useful to recall that in the current model, the solid and liquid phases are defined according to the liquid fraction value.Therefore, a mushy zone state is assigned when the value of ε is between zero and one.In this case, the velocity field is partially bounced back and the macroscopic velocity is modified as in Equation ( 19) [22].The procedure of the implementation of the partial bounce-back approach is described in [23].Note that, to implement BCs in the current method, we are led to convert them, at the mesoscopic level, in terms of the distribution function.

Boundary Conditions
In the LBM framework, implementing boundary conditions is a delicate task, because of the necessity of imposing conditions in terms of particle distributions i f and i g .
For the velocity field, the non-slip boundary conditions are used for all four walls of the cavity.They are realized by the on-grid bounce-back (BB) boundary conditions.The procedure for this rule is to reverse the distribution function of the particle as ( , 1) ( , ) , where w x is a fluid node adjacent to the wall, and i and j represent two opposite lattice directions on the boundary site.Note that the BB conditions apply to the DF in non-parallel directions at a solid wall.
To specify a constant temperature at the left and right walls, we use the method proposed by Inamuro et al. [24].Its principle is to substitute unknown DFs for a boundary point with local equilibrium values using an adjusted temperature to set the defined temperature at that point.
Specifically, the adjusted temperature on the left side can be expressed as v being the computed near-wall velocity, and k g represents a known distribution function.Hence, the unknown DFs p g are computed by 2 ' 1 ./ . As for the adiabatic BCs, the Neumann BCs are achieved using the BB boundary conditions for the distribution i g , as prescribed for i f .

Macroscopic Quantities
Finally, the basic thermo-hydrodynamic properties, such as density, ρ , momentum density, u ρ , and temperature, T , are defined as moments of the DFs, i f and i g , as follows, 8 0 The velocity field is then modified in the mushy zone as, This modified velocity is in accordance with the requirements of the partial bounce-back approach [23].Accordingly, the flow in the mushy zone will be treated as flow in a porous medium and, thus, will be governed by the Darcy law, as proven by the derived analytical solutions in [23].

Model Validation
The current model has been implemented to carry out the 2D natural convection and the melting with natural convection.The obtained results have been compared with selected models that have adopted the same LBM scheme: Huber et al. [15] and Jourabian et al. [16].

Comparison with Similar LBM Schemes
This subject was extensively treated in the literature with sound benchmarks presented in the works of Jany and Bejan [25] and Bertrand et al. [26], to name a few.However, for an LBM implementation, one can refer to the works of Miller et al. [27], Semma et al. [28] and, recently, Su and Davidson [29].Note that our scheme is close to the bounce-back approach proposed by Huber et al. [15] and recently used by Jourabian et al. [16] to define the solid regime; however, our enhancement is to adopt a partial bounce-back approach rather than a full bounce-back.Hence, to check out our model, we resumed the convection melting problem handled in these references, keeping the same settings (see Table 1).Here, the key parameters governing this problem are the Stefan number, the Prandtl number, the Fourier number and the dimensionless time, whose definitions are: Ste ( ) / Table 1.The melting benchmark parameters as in [15,16].x being the position of the melting front at each value of Y .From this figure, it is obvious that our current model reproduces quite closely the results already reported [15,16].In addition, the shape of melting fronts confirms the expected scaling theory proposed by Jany and Bejan [25].
As time proceeds, the upper portion recedes faster due to convection, while the lower part remains upright while being dominated by the conduction ( 0.05 ζ = of Figure 3).However, thereafter, the advection part expands to govern the entire height of the front ( 0.2 ζ = and 0.3 of Figure 3).At the end, a shrinking solid regime is attained as the front reaches the cold wall.

Comparison with Experimental Results and the Finite Volume Method
To further check the model considered here, we conducted a comparison with the experimental results of Gau and Viskanta available for Gallium [30].Moreover, we compared our results for the same problem (melting of a pure Gallium) with the finite volume method (FVM), based on the enthalpy-porosity approach, proposed by Brent et al. [31].The considered problem seeks to examine the 2D melting of pure Gallium in a rectangular cavity heated from one side, while its other sides are adiabatic.The used parameters are gathered in Table 2.In their work, Gau and Viskanta [30] presented detailed traces of the morphology of the melt front at various times.For the propagation of the melting front using this model, our obtained findings via the current model are compared to the results already published (see [30,31]).This comparison is exhibited in Figure 4.As inferred, our results are close to the experimental data, especially in the upper portion of the cavity.However, more discrepancy is exhibited in the lower portion.Thereby, this comparison demonstrates the effectiveness of the approach used here for the melting by convection problems.

Entropy Generation
Unless otherwise stated, the entropy generation is computed by integrating the local entropy production over the whole simulation domain.To achieve this, we have considered the case already handled by Magherbi et al. [6].Here, too, we have adopted the same parameters as those of these authors (see Table 3).Figures 5 and 6 depict temporal evolutions of the total entropy generation and the Bejan number, respectively.It should be noted that this latter parameter remains a constant for each simulation.As can be seen, our results match those of [6].Note that these two parameters exhibit oscillations before reaching the steady state.To sum up, we can state that such a comparison verifies and validates the present model.Table 3.The entropy generation parameters as in [6].).

Results and Analysis
It is useful to recall that the numerical prediction of heat transfer and entropy generation in a 2D square cavity is the main concern of this work.This has been performed for different Pr   In what follows, we focus on the transient entropy generation with phase change.It is well known that during melting, the material undergoes an isothermal phase with a release of latent heat.Here, we ask the question about the effect of this phenomenon on the evolution of the generation of entropy through time.For this purpose, two models with different Pr have been investigated.Following this simulation, a comparison of the losses between the transfer of sensible and latent heat is carried out.

Single-Phase and Solid-Liquid Phase Change with High Prandtl Number
To investigate the effect of phase change on the entropy generation, we conduct the simulations on octadecane, whose parameters are summarized in Table 4. Octadecane is known for having a relatively high convective coefficient.This increases the heat transfer and the fluid velocities.The generation of latent heat during the melting process increases the energy stored in the system.However, we will investigate whether this is provoked by an increase or decrease in the total entropy generation.1a).Here, these quantities further exhibit oscillations due to the high value of Ra before reaching their asymptotic values.Such behavior is due to boundary conditions [3].At first, the isotherms are almost vertical and parallel to the walls since the heat transfer is dictated by the conduction.As time passes, the advection mode becomes more dominant and distorts isotherms, inducing vertical temperature gradients.This will result in a decrease in horizontal temperature gradients, which eventually become locally negative.Consequently, the central streamline elongates to give rise to a second eddy.Thus, for smaller Rayleigh numbers, the steady state is sufficiently close to the equilibrium state, allowing the system to return towards the steady state.Thereby, Prigogine's theorem of minimum entropy production is satisfied.However, for larger Ra 's, the equilibrium state being relatively far, the system oscillates increasingly.Note that this finding has already been reported by Magherbi [6].We also observe that the steady state is reached sooner at larger Ra (Figure 10), whereas at lower Ra (Figure 9), oscillations of transient entropy generation are both apparent and smooth.This is due to the effect of the heat transfer coefficient.For more information on this, we plot in Figures 11 and 12 the Nusselt number (Nu) versus the dimensionless time for the cases of 4 Ra 10 = and 5  10 .The rate of variation of Nu affects the losses due to heat transfer.However, the increase of Nu values at steady state has a significant effect on the fluid velocities and, hence, fluid friction losses.Figure 13 shows the Nusselt number (Nu) versus the different values of (Ra).We deduce that somewhere after 4 Ra 10 = , the heat transfer coefficient at steady state increases dramatically.This leads to higher fluid velocities and, therefore, higher fluid friction losses, as revealed in Figure 10.
Figure 14 shows the temporal variation of the entropy generation in the case of melting with convection.Note that the entropy generation oscillates at the beginning with about a 25% decrease compared to the case of heat transfer with no phase change.Nevertheless, at the steady state, the total entropy generation decreases by more than 50% with the presence of phase change and latent heat dissipation, as illustrated in Figure 14.This is also justified by tracing the Nu versus time in      It also can be stated that, as the melting interface proceeds, the total entropy generation keeps decreasing until reaching a stable state.This is likely due to the relatively low diffusivity, where the system exhibits higher losses due to the small distance between temperature gradients at the beginning of the melting, as illustrated in Figure 16.

Conclusions
This study deals with the analysis of entropy generation during natural convection within a differentially-heated cavity.In this process, the entropy generation for sensible and latent heat transfers has been investigated at different Prandtl numbers via a thermal lattice Boltzmann method.The influences of the Rayleigh number, the Prandtl number and the irreversibility distribution on the entropy generation are assessed.Within the range considered for the key problem parameters, it is shown that entropy generation is controlled by the heat transfer losses for low Pr .However, for large for large Prandtl numbers, the variation of T S is dominated by , l d S , because , l h S has reached its asymptotic value.In addition, it turned out that the increase of the Rayleigh number is required to ensure better heat transfer.Yet, the system will be prone to more losses.Likewise, a suitable choice of the Rayleigh number ( Ra ) and the irreversibility distribution ( ϕ ) will permit getting an indication of a possible optimum design of the enclosure.Moreover, the presence of the latent heat state decreases the overall thermodynamic losses while increasing the quantity of heat transferred.

Figure 1 .
Figure 1.Schematic view of the differentially-heated square cavity with initial and boundary conditions (a) Natural convection (no melting); (b) melting with convection.
lattice sound speed.It should be noted that, the lattice viscosity and diffusivity are selected so as to conform to the intended Prandtl number Pr ( / ) additional force term related to Boussinesq force Fb is incorporated in the model by shifting the velocity field by a term of b f F τ ρ , as proposed by Shan and Chen[20], where

Figure 2
Figure2shows the average position evolution of the melting front m X ( ) 0

Figure 2 .
Figure 2. The average position evolution of the melting front versus ζ .Comparison with the results of [15,16].

Figure 3 .
Figure 3. Position and shape of the melting front at different dimensionless times ζ .

SFigure 7 .SS
Figure 7. (a) T S , , l h S , , l d S and Nu vs. Pr ; (b) zoom on T S , , l h S , , l d S and Nu vs. Pr ( = 0.02 till 10).

Figure 9 .
Figure 9. Entropy generations with no phase change (problem in Figure 1a) for Pr 50 = and

Figure 15 .
The values of Nu decrease in comparison with the case of no phase change.The reason for that is that the fluid volume increases with time as the melting front recedes.Here, we note that, in reality, we can define a Rayleigh number at each time step depending on the position of the melting front ( m X ) by

Figure 13 .
Figure 13.Variation of Nu with Ra (problem in Figure 1a) for Pr 50 = .

Figure 15 .
Figure 15.Nu versus dimensionless time (problem in Figure 1b) for Pr 50 = and

Figure 16 .
Figure 16.Entropy generation versus m X during latent heat transfer with Pr 50 = and Ra = 10 5 .

Table 4 .
The simulation parameters of Octadecane.Figures 9 and 10 show the temporal evolution of the entropy generation during heat transfer with no phase change (the problem described in Figure