Entropy generation and irreversibility analysis due to steady mixed convection flow in a vertical porous channel

The present problem considers steady mixed convection flow of viscous incompressible fluid in a vertical channel formed by two infinite parallel porous plates. The $x^{\prime}$ axis is taken vertically parallel to one of the porous plates of the channel and normal to the  $y^{\prime}$ axis. The porous plates are stationary and parallel to each other at distance 2$l$ apart as shown in Figure 1. The fluid flow is set up due to the applied pressure gradient along channel walls as well as density change caused by the asymmetric heating of the channel boundary porous plates and under the action of gravitational force, hence, the present situation describes a mixed convection flow in a vertical porous channel.

1.png

Figure 1. Schematic diagram of the flow

In addition, heat transfer in the system is due to the isothermal heating of one of the porous plates as well as viscous dissipation within the channel and under the action of heat generation/absorption. Also, we assume the flow to be steady and fully developed; hence, the temperature and velocity fields are functions of $y^{\prime}$ alone.

The governing equations for the steady flow of viscous incompressible fluid between two heated parallel plates are the conservation of mass

$\nabla \vec{V}=0$     (1)

conservation of momentum

$(\vec{V} . \nabla) \vec{V}=v \nabla^{2} \vec{V}-\frac{1}{\rho} \nabla P+g \beta\left(T^{\prime}-T_{0}\right)$       (2)

and the conservation of energy

$(\vec{V} . \nabla) T=\frac{k}{\rho C_{P}} \nabla^{2} T+\frac{\mu}{\rho C_{P}}(\nabla u)^{2}+\frac{Q_{0}\left(T^{\prime}-T_{0}\right)}{\rho C_{P}}$       (3)

where $\nabla=\frac{\partial}{\partial x} i+\frac{\partial}{\partial y} \underline{j}+\frac{\partial}{\partial z} \underline{k}$ . The physical quantities, $\rho$ ，$\vec{V}$ ， $\mathbf{v}$ ， $k$ and P are defined in the nomenclature.

Considering a two dimensional velocity so that $\vec{V}=\left(u^{\prime}, v^{\prime}, 0\right)$ , where $\mathcal{u}^{\prime}$and $\boldsymbol{v}^{\prime}$ are the vertical and horizontal (suction/injection) components of the velocity respectively. Since the flow is assumed to be fully developed, the conservation of mass equation (1) is reduced to

$\frac{d v^{\prime}}{d y^{\prime}}=0$      (4)

Upon integration of eq. (4), we obtain the horizontal velocity as a constant (say, $v^{\prime}=-v_{0}$ ) which is the velocity of suction/injection.

The equations of motion for steady mixed convection flow of a viscous incompressible heat generating fluid with viscous dissipation are given as,

$v \frac{d^{2} u^{\prime}}{d y^{2}}+v_{0} \frac{d u^{\prime}}{d y^{\prime}}-\frac{1}{\rho} \frac{\partial P}{\partial x^{\prime}}+g \beta\left(T^{\prime}-T_{0}\right)=0$     (5)

$\frac{k}{\rho C_{p}} \frac{d^{2} T^{\prime}}{d y^{2}}+v_{0} \frac{d T^{\prime}}{d y^{\prime}}+\frac{\mu}{\rho C_{P}}\left(\frac{d u^{\prime}}{d y^{\prime}}\right)^{2}+\frac{Q_{0}\left(T^{\prime}-T_{0}\right)}{\rho C_{P}}=0$  (6)

while the boundary conditions are

$u^{\prime}=0, T^{\prime}=T_{0}$ at $y^{\prime}=-l$

$u^{\prime}=0, T^{\prime}=T_{w}$ at $y^{\prime}=l$$\}$    (7)

The following non-dimensional quantities are used to transform equations (5) – (7) to dimensionless form.

$y=\frac{y^{\prime}}{l}, \theta=\frac{T^{\prime}-T_{0}}{T_{w}-T_{0}}, S=\frac{v_{0} l}{v}, \operatorname{Pr}=\frac{\mu C_{P}}{k}, E c=\frac{U^{2}}{C_{P}\left(T_{w}-T_{0}\right)}$

$u=\frac{u^{\prime}}{u_{0}^{\prime}}, \quad G r=\frac{g \beta l^{3}\left(T_{w}-T_{0}\right)}{v^{2}}, \operatorname{Re}=\frac{l u_{0}}{v}, \delta=\frac{Q_{0} l^{2}}{k}$$\}$    (8)

Substituting equation (8) into equations (5) and (6), the momentum and energy equations are rendered in dimensionless form as

$\frac{d^{2} u}{d y^{2}}+S \frac{d u}{d y}-\gamma+\frac{G r}{\operatorname{Re}} \theta=0$      (9)

$\frac{d^{2} \theta}{d y^{2}}+S \operatorname{Pr} \frac{d \theta}{d y}+E c \operatorname{Pr}\left(\frac{d u}{d y}\right)^{2}+\delta \theta=0$    (10)

and the boundary conditions are given as

$u=0, \theta=0$ at $y=-1$

$u=0, \theta=1$ at $y=1$$\}$    (11)

where $\gamma=\frac{l^{2}}{\mu u_{0}} \frac{\partial P}{\partial x^{\prime}}$ 

For brevity, we shall write the mixed convection parameter $G r / \mathrm{Re}$ as Gre . This parameter measures the contribution of either the forced convection (small value of Gre) or natural convection (large value of Gre ). For mixed convection, a moderate value of Gre is required. $P r$ is the Prandtl number which is inversely proportional to the thermal diffusivity of the fluid. S is the dimensionless suction/injection parameter, where positive values denote suction at the porous plate $\left(y^{\prime}=-l\right)$ with a corresponding injection on the plate $y^{\prime}=l$ while negative values denote injection at the porous plate $y^{\prime}=-l$ with a corresponding suction on the other plate. Eckert number ( $E c$ ) is the measure of viscous dissipation in the system while  $\delta$ is the temperature dependent heat source/sink parameter, positive values of $\delta$  represent heat source and negative values represent heat sink. All the physical quantities used in the dimensionless analysis are defined in the nomenclature.

In order to solve equations (9) and (10) subject to the boundary conditions (11), we use the Homotopy perturbation method (HPM) to obtain approximate analytical solutions to the problem. HPM is a powerful solution tool to solve linear as well as nonlinear equations. The method has been shown to perform excellently well in solving this type of problems [28] and is advantageous over the regular perturbation technique because it was able to overcome the small parameter restriction which characterizes the regular perturbation technique.

2.1 Homotopy perturbation method

In order to illustrate the basic ideas of the HPM, we consider the following nonlinear differential equation

$A(u)-f(r)=0, r \in \Omega$    (12)

$B\left(u, \frac{\partial u}{\partial r}\right)=0, r \in \Gamma$      (13)

where $A$ is a general differential operator, $B$ is a boundary operator, $f(r)$ is a known analytic function, $\Gamma$ is the boundary of the domain $\Omega$ . The operator $A$ can be divided into two parts $L$ and $N$ , where $L$ is linear and $N$ is nonlinear. Hence, eq. (12) can be rewritten as

$L(u)+N(u)-f(r)=0$        (14)

Using the Homotopy perturbation technique, we construct a Homotopy $v(r, P) : \Omega \times[0,1] \rightarrow R$ which satisfies

$H(v, P)=(1-P)\left[L(v)-L\left(u_{0}\right)\right]+P[A(v)-f(r)]=0, P \in[0,1], r \in \Omega$       (15)

or

$\mathrm{H} \quad(v, P)=L(v)-L\left(u_{0}\right)+P L\left(u_{0}\right)+P[N(u)-f(r)]=0$     (16)

where $P \in[0,1]$ is called the Homotopy parameter, $u_{0}$ is the initial approximation for the solution of eq. (12) which satisfies the boundary conditions. Clearly, from eqs. (15) and (16) we have;

$\mathrm{H} \quad(v, 0)=L(v)-L\left(u_{0}\right)=0$       (17)

$\mathrm{H} \quad(v, 1)=A(v)-f(r)=0$        (18)

Assuming that the solution of eqs. (15) and (16) are expressed as a power series in $P$ , that is

$v=v_{0}+P v_{1}+P^{2} v_{2}+\ldots$        (19)

Setting $P=1$ gives the approximate solution of eq. (12) as

$u=\lim _{P \rightarrow 1} v=v_{0}+v_{1}+v_{2}+\cdots$        (20)

Applying the Homotopy perturbation technique to solve the governing equations in the present problem, we construct Homotopy on eqs. (9) and (10) to get

$H(u, P)=(1-P)\left(\frac{d^{2} u}{d y^{2}}+G r e \theta\right)+P\left(\frac{d^{2} u}{d y^{2}}+S \frac{d u}{d y}+G r e \theta-\gamma\right)=0$      (21)

$H(\theta, P)=(1-P)\left(\frac{d^{2} \theta}{d y^{2}}\right)+P\left(\frac{d^{2} \theta}{d y^{2}}+S \operatorname{Pr} \frac{d \theta}{d y}+E c \operatorname{Pr}\left(\frac{d u}{d y}\right)^{2}+\delta \theta\right)=0$        (22)

Simplifying

$\frac{d^{2} u}{d y^{2}}+G r e \theta+P\left(S \frac{d u}{d y}-\gamma\right)=0$        (23)

$\frac{d^{2} \theta}{d y^{2}}+P\left(S \operatorname{Pr} \frac{d \theta}{d y}+E c \operatorname{Pr}\left(\frac{d u}{d y}\right)^{2}+\delta \theta\right)=0$         (24)

Assume the solutions of eqs. (9) and (10) to be written as

$u=u_{0}+P u_{1}+P^{2} u_{2}+\ldots$

$\theta=\theta_{0}+P \theta_{1}+P^{2} \theta_{2}+\ldots$$\}$         (25)

Substituting eq. (25) into eqs. (23) and (24) and simplifying, we have the following sets of boundary value problems which are obtained by comparing the coefficients of powers of $P$ .

Zeroth order:

$\frac{d^{2} \theta_{0}}{d y^{2}}=0$        (26)

$\frac{d^{2} u_{0}}{d y^{2}}+G r \theta_{0}=0$        (27)

$\begin{array}{ll}{\theta_{0}(-1)=0,} & {u_{0}(-1)=0} \\ {\theta_{0}(1)=1,} & {u_{0}(1)=0}\end{array}$        (28)

First order:

$\frac{d^{2} \theta_{1}}{d y^{2}}+S \operatorname{Pr} \frac{d \theta_{0}}{d y}+E c \operatorname{Pr}\left(\frac{d u_{0}}{d y}\right)^{2}+\delta \theta_{0}=0$      (29)

$\frac{d^{2} u_{1}}{d y^{2}}+G r e \theta_{1}+S \frac{d u_{0}}{d y}-\gamma=0$      (30)

$\begin{array}{ll}{\theta_{1}(-1)=0,} & {u_{1}(-1)=0} \\ {\theta_{1}(1)=0,} & {u_{1}(1)=0}\end{array}$       (31)

Second order:

$\frac{d^{2} \theta_{2}}{d y^{2}}+S \operatorname{Pr} \frac{d \theta_{1}}{d y}+2 E c \operatorname{Pr}\left(\frac{d u_{0}}{d y}\right)\left(\frac{d u_{1}}{d y}\right)+\delta \theta_{1}=0$         (32)

$\frac{d^{2} u_{2}}{d y^{2}}+G r e \theta_{2}+S \frac{d u_{1}}{d y}=0$         (33)

$\begin{array}{ll}{\theta_{2}(-1)=0,} & {u_{2}(-1)=0} \\ {\theta_{2}(1)=0,} & {u_{2}(1)=0}\end{array}$         (34)

Solving the sets of Boundary Value Problems and taking limit as $P \rightarrow 1$ , the solution of the momentum and energy equations are

$u(y)=u_{0}(y)+u_{1}(y)+u_{2}(y) \cdots$     (35)

$\theta(y)=\theta_{0}(y)+\theta_{1}(y)+\theta_{2}(y)+\cdots$         (36)

where:

$u_{0}(y)=a_{1} y^{3}+a_{2} y^{2}+c_{3} y+c_{4}$

$u_{1}(y)=b_{1} y^{8}+b_{2} y^{7}+b_{3} y^{6}+b_{4} y^{5}+b_{5} y^{4}+b_{6} y^{3}+b_{7} y^{2}+c_{7} y+c_{8}$

$u_{2}(y)=g_{1} y^{3}+g_{2} y^{12}+g_{3} y^{11}+g_{4} y^{10}+g_{5} y^{9}+g_{6} y^{8}+g_{7} y^{7}+g_{8} y^{6}+g_{9} y^{5}+$ $g_{10} y^{4}+g_{11} y^{3}+g_{12} y^{2}+c_{11} y+c_{12}$

$\theta_{0}(y)=\frac{1}{2} y+\frac{1}{2}$

$\theta_{1}(y)=a_{3} y^{6}+a_{4} y^{5}+a_{5} y^{4}+a_{6} y^{3}+a_{7} y^{2}+c_{5} y+c_{6}$

$\theta_{2}(y)=d_{1} y^{11}+d_{2} y^{10}+d_{3} y^{9}+d_{4} y^{8}+d_{5} y^{7}+d_{6} y^{6}+d_{7} y^{5}+d_{8} y^{4}+d_{9} y^{3}+d_{10} y^{2}+c_{9} y+c_{10}$

The constants $a_{i}(i=1,2, \cdots 7)$ , $b_{i}(i=1,2, \cdots 7)$ , $c_{i}(i=1,2, \cdots 12)$  $d_{i}(i=1,2, \cdots 10)$ and $g_{i}(i=1,2, \cdots 12)$ are defined in the appendix.

2.2 Entropy generation in the system

For the flow of a Newtonian incompressible fluid due to Fourier law of heat conduction, the volumetric rate of entropy generation rate is given by Bejan² in Cartesian coordinates as

$E_{G}=\frac{k}{T_{0}^{2}}\left(\left(\frac{\partial T}{\partial x}\right)^{2}+\left(\frac{\partial T}{\partial y}\right)^{2}\right)+\frac{\mu}{T_{0}}\left(2\left\{\left(\frac{\partial u}{\partial x}\right)^{2}+\left(\frac{\partial v}{\partial y}\right)^{2}\right\}+\left[\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right]^{2}\right)$      (37)

This form of entropy generation shows that the irreversibility results from two effects; conductive $(k)$ and viscosity $(\mu)$ . Entropy generation rate is finite and positive whenever one of temperature or velocity gradients is present in a medium. Velocity and temperature distribution are often simplified in many fundamental convective heat transfer problems, assuming that the flow is hydro-dynamically developed $\left(\frac{\partial \vec{V}}{\partial x}=0\right)$ and thermally developing $\left(\frac{\partial T}{\partial x} \neq 0\right)$ or thermally developed $\left(\frac{\partial T}{\partial x}=0\right)$ ( White [29]  and Burmeister [30]), then eq. (37) is reduced to the form:

$E_{G}=\frac{k}{T_{0}^{2}}\left(\left(\frac{\partial T}{\partial x}\right)^{2}+\left(\frac{\partial T}{\partial y}\right)^{2}\right)+\frac{\mu}{T_{0}}\left(\frac{\partial u}{\partial y}\right)^{2}$        (38)

By scaling the velocity $u$ with the reference velocity $u_{0}$ , the distance $y$ with $l$, the distance $x$ with $l^{2} u_{0} \alpha^{-1}$ and expressing the dimensionless temperature $\theta$ as $\left(T-T_{0}\right) \Delta T^{-1}$ , the entropy generation is rendered in dimensionless form as

$N_{S}=\frac{1}{P e^{2}}\left(\frac{\partial \theta}{\partial x}\right)^{2}+\left(\frac{\partial \theta}{\partial y}\right)^{2}+\frac{B r}{\Omega}\left(\frac{\partial u}{\partial y}\right)^{2}=N_{c}+N_{y}+N_{f}$   (39)

where $P e=l u_{0} \alpha^{-1}$ is the Peclet number, and $\Omega=\Delta T / T_{0}$ is the dimensionless temperature difference. The first term $\left(N_{c}\right)$ stands for the entropy generation by heat transfer due to axial conduction, the middle term $\left(N_{y}\right)$ represents entropy generation due to heat transfer across different fluid sections within the channel, and the last term $\left(N_{f}\right)$ gives the contribution of viscous dissipation to entropy generation.

 Substituting eqs. (35) and (36) into the dimensionless expression for entropy as presented in eq. (39), we have  

$N_{s}=\lambda_{3}^{2}+B r \Omega^{-1} \lambda_{4}^{2}$       (40)

where:  

$\frac{d \theta}{d y}=\lambda_{3}=11 d_{1} y^{10}+10 d_{2} y^{9}+9 d_{3} y^{8}+8 d_{4} y^{7}+7 d_{5} y^{6}+\left(6 a_{3}+6 d_{6}\right) y^{5}$ $+\left(5 a_{4}+5 d_{7}\right) y^{4}+\left(4 a_{5}+4 d_{8}\right) y^{3}+\left(3 a_{6}+3 d_{9}\right) y^{2}+\left(2 a_{7}+2 d_{10}\right) y+\left(\frac{1}{2}+c_{5}+c_{9}\right)$

$\frac{d u}{d y}=\lambda_{4}=13 g_{1} y^{12}+12 g_{2} y^{11}+11 g_{3} y^{10}+10 g_{4} y^{9}+9 g_{5} y^{8}+\left(8 b_{1}+8 g_{6}\right) y^{7}$ $+\left(7 b_{2}+7 g_{7}\right) y^{6}+\left(6 b_{3}+6 g_{8}\right) y^{5}+\left(5 b_{4}+5 g_{9}\right) y^{4}+\left(4 b_{5}+4 g_{10}\right) y^{3}$ $+\left(3 a_{1}+3 b_{6}+3 g_{11}\right) y^{2}+\left(2 a_{2}+2 b_{7}+2 g_{12}\right) y+\left(c_{3}+c_{7}+c_{11}\right)$

2.3 Irreversibility distribution ratio

Entropy generation rate is accounted for by both heat transfer and fluid friction in a convective problem. Eq. (39) shows the level of distribution of entropy but is unable to point out the relative contributions of each irreversibility to the total entropy generation. To identify whether it is the fluid friction or heat transfer irreversibility that dominates the total entropy, Bejan [1] defined irreversibility distribution ratio $\varphi$ as the ratio of the entropy generation due to fluid friction $\left(N_{f}\right)$ to heat transfer $\left(N_{c}+N_{y}\right)$ . However, Paoletti et al. [31]

showed an alternative irreversibility distribution parameter in terms of the Bejan number (Be) and defined it as the ratio of the entropy generation due to heat transfer $\left(N_{c}+N_{y}\right)$ to the entropy generation $\left(N_{s}\right)$  In the present problem, the flow is assumed to be thermally developed so that heat transfer due to axial conduction becomes negligible so that Bejan number becomes

$B e=\frac{N_{y}}{N_{s}}$        (41)

A critical look at the Bejan number shows that $0 \leq B e \leq 1$ with $B e=1$  indicating that the entropy generation is only by heat transfer irreversibility, whereas, $B e=0$  signify that the total entropy generation is as a result of fluid friction irreversibility only, whenever heat transfer and fluid friction contribute equally to the total entropy generation then $B e=1 / 2$ . In the light of this analysis, heat transfer irreversibility dominates the total entropy generation whenever $B e \rightarrow 1$ and likewise fluid friction irreversibility dominates when $B e \rightarrow 0$      

$B e=\frac{\lambda_{3}^{2}}{\lambda_{3}^{2}+B r \Omega^{-1} \lambda_{4}^{2}}$         (42)

2.4 Pressure gradient

In order to estimate the pressure gradient in the present problem, we assumed that the flow has a constant mass flux so that $\gamma$ is obtained such that the mass flux within the channel is maintained at a constant level.  That is

$\int_{-1}^{1} u(y) d y=2$        (43)

Upon evaluation of eq. (43), we have the following polynomial equation

$n_{4} \gamma^{3}+n_{8} \gamma^{2}+n_{9} \gamma+n_{10}=0$        (44)

where the expressions for $n_{i,} \quad i=0,1, \cdots 11$ are defined in the appendix.

The roots of the cubic equation (44) are obtained numerically using an in-built function on the platform of MATLAB.  

2.5 Validation of results

By setting the heat generation term $\delta$ and the mixed convection term Gre to zero in the present problem, we recover the results of Ajibade et al. [14] The comparison is presented in the table below

Table 1. Numerical comparison between the present problem and case 2 of Ajibade et al. [14] 

A slight deviation is observed in the values displayed in the table and this could be attributed to the different pressure gradients used in both problems. For instance, in the work of Ajibade et al. [14], the pressure gradient was assumed to be constant. However, the pressure gradient in the present problem is obtained from the evaluation of the constant mass flux equation (43). This means that the pressure gradient in the present problem is dependent on all other parameters that affect the velocity of the fluid.

In the present study, entropy generation and irreversibility distribution in a steady fully developed flow and heat transfer with viscous dissipation in a vertical channel have been investigated. The governing momentum and energy equations were solved using the Homotopy perturbation method. The impacts of each of the operating parameters are discussed with the aid of graphs. It is found that heat source and the Grashof number exerts significant influence on the temperature, velocity, as well as entropy generation rate and irreversibility distribution within the channel. The following major conclusions have been drawn from the present study:

Velocity and temperature are higher with increasing values of the parameters

Minimum entropy generation is obtained towards the center of the channel

Entropy generation is higher at the porous walls.

The contribution of fluid friction irreversibility dominates that of heat transfer irreversibility near the channel’s porous walls while heat transfer irreversibility dominates fluid friction irreversibility towards the centerline of the channel.

Finally, a comparison made between this work and Ajibade et al.[14] shows that the present work agrees significantly with Ajibade et al.[14]

2.png

Figure 2. Velocity Profiles for different Gre $(P r=0.71, E c=0.3,S=0.5, \delta=0.2)$

3.png

Figure 3. Velocity Profiles for different $\delta(P r=0.71, E c=0.3 .S=0.5,Gre=10)$

4.png

Figure 4. Temperature Profiles for different Gre$(P r=0.71 . E c=0.3 .S=0.5 .\delta=0.2)$

5.png

Figure 5. Temperature Profiles for different $\delta(P r=0.71, E c=0.3 .S=0.5 .Gre=10)$

6.png

Figure 6. Entropy generation number for different Gre$\left(P r=0.71, E c=0.3, S=0.5, \delta=0.2, B r \Omega^{-1}=1\right)$

7.png

Figure 7. Entropy generation number for different  $\delta\left(P r=0.71, E c=0.3, S=0.5, G r e=10, B r \Omega^{-1}=1\right)$

8.png

Figure 8. Entropy generation number for different $S\left(P r=0.71, E c=0.3, \delta=0.2, \text { Gre }=10, B r \Omega^{-1}=1\right)$

9.png

Figure 9. Entropy generation number for different $E c\left(P r=0.71, S=0.5, \delta=0.2, \text { Gre }=10, \mathrm{Br} \Omega^{-1}=1\right)$

10.png

Figure 10. Entropy generation number for different $\operatorname{Pr}\left(E c=0.3, S=0.5, \delta=0.2, \text { Gre }=10, \mathrm{Br} \Omega^{-1}=1\right)$

11.png

Figure 11. Entropy generation number for different $B r \Omega^{-1}(P r=0.71, S=0.5, \delta=0.2, G r e=10, E c=0.3)$

12.png

Figure 12. Bejan Number for different Gre $\left(P r=0.71, E c=0.3, S=0.5, \delta=0.2, B r \Omega^{-1}=1\right)$ 

13.png

Figure 13. Bejan Number for different $\delta\left(P r=0.71, E c=0.3, S=0.5, \text { Gre }=10, B r \Omega^{-1}=1\right)$ 

14.png

Figure 14. Bejan Number for different $S\left(P r=0.71, E c=0.3, \delta=0.2, G r e=10, B r \Omega^{-1}=1\right)$

15.png

Figure 15. Bejan number for different $E c\left(P r=0.71, S=0.5, \delta=0.2, \text { Gre }=10, \mathrm{Br} \Omega^{-1}=1\right)$

16.png

Figure 16. Bejan number for different $\operatorname{Pr}\left(E c=0.3, S=0.5, \delta=0.2, G r e=10, \mathrm{Br} \Omega^{-1}=1\right)$

17.png

Figure 17. Bejan Number for different $B r \Omega^{-1}(P r=0.71, S=0.5, \delta=0.2, \text { Gre }=10, E c=0.3)$