# Second Law Analysis for a Variable Viscosity Reactive Couette Flow under Arrhenius Kinetics

Authors: N. S. Kobo; O. D. Makinde
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

This study investigates the inherent irreversibility associated with the Couette flow of a reacting variable viscosity combustible material under Arrhenius kinetics. The nonlinear equations of momentum and energy governing the flow system are solved both analytically using a perturbation method and numerically using the standard Newton Raphson shooting method along with a fourth-order Runge Kutta integration algorithm to obtain the velocity and temperature distributions which essentially expedite to obtain expressions for volumetric entropy generation numbers, irreversibility distribution ratio, and the Bejan number in the flow field.


---

## Body

## 1. Introduction

In fluid dynamics, Couette flow refers to the laminar flow of a viscous fluid in the space between two parallel plates, one of which is moving relative to the other. This type of flow is named in honor of Maurice Marie Alfred Couette, a Professor of Physics at the French University of Angers in the late 19th century [1, 2]. Couette flow occurs in fluid machinery involving moving parts and is especially important for hydrodynamic lubrication [3].  If the surfaces are smooth and flat with constant fluid properties, the solution is the simple linear velocity distribution, with a drag proportional to the relative velocity and inversely proportional to the gap width [1, 4]. It is an important classical example of exact solutions for Navier-Stokes equation. However, most fluids used in engineering and industrial systems like coal slurries, polymer solutions or melts, drilling mud, hydrocarbon oils, grease, and so forth are  chemically reactive and can be subjected to extreme conditions, such as high temperature, pressure, and shear rate during processing [5, 6]. In  fact,  viscous  heating  produced  due  to  friction between  the  fluid  and  the  surrounding  walls  coupled  with  high  shear rates  and Arrhenius kinetics can  lead  to  high temperature being generated within the fluid [7]. This may have a significant effect on the fluid properties. It  is  well  known  that  the  most  sensitive  fluid  property  to  temperature  rise  is  the viscosity [3, 5].  For  instance, the  viscosity of  various  lubricants  used  in  engineering systems  like polymer  solutions,  mineral  oils  with  polymeric  additives,  and so forth  varies  with  temperature. This variation in the fluid viscosity due to temperature may affect the flow characteristics as well as the efficient operation of industrial machinery where lubrication is important [1, 7].  Hence, it is necessary to ensure that the viscosity of such lubricants is at all times maintained at optimum levels.From the application point of view, the determination of thermal criticality in a flow system is extremely important. For example, special attention must be paid to the heating of lubricant by the frictional force and Arrhenius kinetics since viscosity is temperature dependent. Thermal criticality occurs when the rate of heat generation within the  flow  system  exceeds  the  heat  dissipation  to  the  surroundings  [2, 6].  This  condition  is incipient  thermal  runaway  or  ignition  in  the  flow  system  [7, 8].  A  primary  objective  of  thermal criticality  analysis  is  the  prediction  of  the  critical  or  unsafe  flow  conditions  in  order  to  avoid them [9]. One method of accomplishing this is to cycle the lubricant through a cooling reservoir in order to maintain the desired viscosity of the fluid. Another way of handling the excessive heat generation problem is to use commercially available additives to decrease the viscosity’s temperature dependence.Meanwhile, efficiency calculation of heat exchange systems has been very much restricted to the first law of thermodynamics. Calculations using the second law of thermodynamics, which is related to entropy generation, are more reliable than first law-based calculations. Therefore, the second law of thermodynamics can be applied to investigate the entropy generation rate in the flow system due to fluid friction and heat transfer. The determination of entropy generation is also important in upgrading the system performance because the entropy generation is the measure of the destruction of the available work of the system [4]. As entropy generation takes place, the quality of energy decreases. It is important to study the distribution of the entropy generation within the fluid volume for preserving the quality of energy in fluid flow processes or reducing the entropy generation. Entropy generation method as a measure of system performance was first introduced by Batchelor [1] and Bejan [10]. Thereafter, considerable research works were carried out by several authors on the application of the second law of thermodynamics to various aspects of fluid flow and heat transfer problems [4, 10–17]. Moreover, to the best of our knowledge, no study has focused on the combined effects of temperature-dependent fluid viscosity and Arrhenius kinetics on thermal stability and entropy generation in flow systems which are of great importance in many engineering fields such as heat exchangers, cooling of nuclear reactors, automobile lubrication, energy storage systems, and cooling of electronic devices.In  this  study,  the  variable  viscosity  reactive  Couette  flow  is  considered  and  the  inherent  irreversibility together with thermal criticality in the flow system is investigated. The plan of this paper is as follows: in Section2 we describe the theoretical analysis of the problem with respect to the fluid velocity and temperature fields. In Sections 3–5 we introduce and apply some rudiments of perturbation technique coupled with Hermite-Padé approximation procedure and standard Newton Raphson shooting method along with a fourth-order Runge Kutta integration algorithm in order to obtain the fluid velocity temperature profiles as well as criticality conditions in the system. Section 6 describes the volumetric entropy generation rate, irreversibility distribution ratio, and the Bejan number. The results are presented graphically and discussed quantitatively with appropriate physical explanations in Section 7. Finally, the conclusion is outlined in Section 8.


## 2. Problem Formulation

The configuration of the problem studied in this paper is depicted in Figure1. The fluid is assumed to be viscous, incompressible, reactive, and flowing steadily in the x̅-direction between two parallel plates of width H and length L.The upper plate is moving with constant velocity U while the lower plate is kept stationary. Following [3, 8], the temperature-dependent viscosity (μ̅) and the chemical reaction kinetic (G) functions can be expressed in Arrhenius type [18] as(2.1)μ̅=μ0eE/RT̅,G=QC0Ae-E/RT̅,
where E is the activation energy, R is the universal gas constant, Q is the heat of reaction, A is the rate constant, C0 is the initial concentration of the reactant species, and μ0 is the fluid reference dynamic viscosity at a very large temperature (i.e., as T̅→∞).Figure 1
Schematic diagram of the problem.Under these conditions the continuity, momentum, and energy equations governing the problem in dimensionless form may be written in Cartesian coordinate(x,y) as [1, 2, 5];(2.2)∂u∂x+∂v∂y=0,ε2Re(u∂u∂x+v∂udy)=-∂p∂x+2ε2∂∂x(μ∂u∂x)+∂∂y[μ(∂u∂y+ε2∂v∂x)],ε4Re(u∂v∂x+v∂v∂y)=-∂p∂y+2ε2∂∂y(μ∂v∂y)+ε2∂∂x[μ(∂u∂y+ε2∂v∂x)],ε2Pe(u∂T∂x+v∂T∂y)=ε2∂2T∂x2+∂2T∂y2+λe(T/(1+βT))+μΦ,
where(2.3)Φ=Br[2ε2(∂u∂x)2+2ε2(∂v∂y)2+(∂u∂y+ε2∂v∂x)2].
We have employed the following nondimensional quantities in (2.2):(2.4)y=y̅εL,x=x̅L,u=u̅U,v=v̅εU,ε=HL,μ=μ̅μ0e-E/RTw,T=E(T̅-Tw)RTw2,P=ε2Lp̅μ0U,β=RTwE,Br=μ0EU2kRTw2e-E/RTw,Pe=ρcpLUk,λ=QC0AEH2kRTw2e-E/RTw,Re=ρULμ0,
where μ,ρ,κ are the dynamic viscosity, fluid density, and thermal conductivity, respectively, T is the fluid temperature, Tw is the plate surface temperature, uis the axial velocity,vis  the normal velocity, cp is the specific heat at constant pressure, p is the pressure,(x,y) are distances measured in streamwise and normal directions, respectively, U is the velocity scale, Pe is the Peclet number, β is the activation energy parameter, Br is the Brinkman number, λ is the Frank-Kamenetski parameter, and Re is the Reynolds number. Since the channel aspect ratio is small (0<ε≪1), the lubrication approximation based on an asymptotic simplification of the governing equation (2.2) is invoked. For Couette flow, the axial pressure gradient is zero (i.e., ∂p/∂x=0) and the flow is solely driven by the uniform motion of the upper plate. Equations (2.2) then become(2.5)0=∂∂y(μ∂u∂y)+O(ε2),(2.6)0=∂p∂y+O(ε2),(2.7)0=∂2T∂y2+μBr(∂u∂y)2+λe(T/(1+βT))+O(ε2),
where μ=e-T/(1+βT). The appropriate boundary conditions in dimensionless form are given as follows: (2.8)u=0,T=0,(foralowerfixedimpermeableplate)aty=0,(2.9)u=1,T=0,(theupperplateissubjectedtoauniformmotion)aty=1.
The boundary conditions, (2.8) and (2.7), can be easily combined to give(2.10)dudy=meT/(1+βT),d2Tdy2+γeT/(1+βT)=0,
where γ=m2Br+λ and m is a constant to be determined. In the following sections, (2.10) is solved both analytically using a perturbation method and numerically using the standard Newton Raphson shooting method along with a fourth-order Runge Kutta integration algorithm [19].


## 3. Perturbation Approach

Due to the nonlinear nature of the velocity and temperature field equations in (2.10), it is convenient to form a power series expansion both in the parameter γ, that is,(3.1)u=∑i=0∞uiγi,T=∑i=0∞Tiγi,m=∑i=0∞miγi.
Substituting the solution series in (3.1) into (2.10) and collecting the coefficients of like powers of γ, we obtained the followings:Order Zero(γ0). One has
(3.2)du0dy=m0eT0/(1+βT0),d2T0dy2=0,
with u0(0)=0, T0(0)=0,  u0(1)=1, and T0(1)=0.Order One(γ1). One has
(3.3)du1dy=(m1+m0T1(1+βT0)2)eT0/(1+βT0),d2T1dy2+eT0/(1+βT0)=0,
with u1(0)=0, T1(0)=0,  u1(1)=0, and T1(1)=0.Order Two(γ2). One has
(3.4)du2dy=(m2+m1T1(1+βT0)2)eT0/(1+βT0)-m0eT0/(1+βT0)2(1+βT0)4(2T12β+2T0T12β2-T12-2T2-4T0T2β-2T2T02β2),d2T2dy2+T1(1+βT0)2eT0/(1+βT0)=0,
with u2(0)=0, T2(0)=0,  u2(1)=0, T2(1)=0,and so on. The above equations for the coefficients of solution series are solved iteratively for the velocity and temperature fields, and we obtain(3.5)T(y)=-12y(y-1)γ+124y(y-1)(y2-y-1)γ2+11440y(y-1)(12βy4-8y4-24βy3+12y3+6βy2+y2-9y+6yβ-9+6β)γ3+O(γ4),(3.6)u(y)=y-112y(2y-1)(y-1)γ-1360y(2y-1)(y-1)(3y2-3y-1)(-2+3β)γ2-160480y(2y-1)(y-1)(-756βy4+540β2y4+204y4-1080y3β2-408y3+1512βy3+63y2-367βy2+378y2β2+141y+162yβ2-378yβ+54β2+47-126β)γ3+O(γ4).
It is noteworthy that, in the limit of γ →0, the fluid velocity profile obtained in (3.6) reduces to u(y) = ywhich corresponds to the classical linear velocity profile for Couette flow with constant fluid viscosity. Using a computer symbolic algebra package (MAPLE) [19], the first few terms of the above solution series in (3.5)-(3.6) are obtained. We are aware that these power series solutions are valid for very small parameter values (of order 10-4). However, using Hermite-Padé approximation  technique,  we  have  extended  the  usability  of  the  solution  series  beyond  small parameter values as illustrated in the following section.


## 4. Hermite-Padé Approximation Technique

From the application point of view, it is extremely important to determine the appearance of criticality or nonexistence of steady-state solution for certain parameter values. In order to achieve this, we first derived a special type of Hermite-Padé approximant. Let(4.1)UN(γ)=∑n=0Nanγn+O(γN+1),asγ→0,
be a given partial sum. It is important to note here that (4.1) can be used to approximate any output of the solution of the problem under investigation (e.g., the series for the wall heat flux parameter in terms of Nusselt number Nu = -dT/dy at y = 1), since everything can be Taylor expanded in the given small parameter. Assume that U(γ) is a  local representation of an algebraic function of γ in the context of nonlinear problems; then we seek an expression of the form(4.2)Fd(γ,U)=∑m=1d∑k=0mfm-k,kγm-kUk,
of degree d≥ 2, such that(4.3)∂Fd∂U(0,0)=1,Fd(γ,UN)=O(γN+1),asγ→0.
The requirement (4.3) reduces the problem to a system of Nlinear equations for the unknown coefficients of Fd. The entries of the underlying matrix depend only on the N given coefficients an and we will take N=(d2+3d-2)/2, so that the number of equations equals the number of unknowns. The polynomial Fd is a special type of Hermite-Padé approximant and is then investigated for bifurcation and criticality conditions using Newton diagram; Vainberg and Trenogin see [20].


## 5. Numerical Approach

The numerical technique chosen for the solution of the coupled ordinary differential (2.10) is the standard Newton Raphson shooting method along with a fourth-order Runge Kutta integration algorithm. Equation (2.10) is transformed into a system of first-order differential equations as follows. Let u=x1, T=x2, and T′=x3, where the prime symbol represents derivative with respect to y. Then, the problem becomes
(5.1)x1′=mex2/(1+βx2),x2′=x3,x3′=-γex2/(1+βx2),
subject to the following initial conditions:(5.2)x1(0)=0,x2(0)=0,x3(0)=s1.
The unspecified initial condition s1 and the undetermined constant m are guessed systematically and (5.1) is then integrated numerically as initial valued problems to the given terminal point at y=1. For each set of parameter values for β and γ, the procedure is repeated until conditions at the y = 1 (i.e., x1(1)=1, x2(1)=0) are satisfied and the desired degree of accuracy (namely, 10-7) of the results obtained is achieved.


## 6. Entropy Analysis

Flow and heat transfer processes between two parallel plates are irreversible. The nonequilibrium conditions arise due to the exchange of energy and momentum within the fluid and at solid boundaries, thus resulting in entropy generation. A part of the entropy production is due to the heat transfer in the direction of finite temperature gradients and the other part of entropy production arises due to the fluid friction. The general equation for the entropy generation per unit volume is given by [2, 10, 12, 16]
(6.1)Sm=kTw2(∇T̅)2+μTwΦ.
The first term in (6.1) is the irreversibility due to heat transfer and the second term is the entropy generation due to viscous dissipation. Using (6.1), we express the entropy generation number in dimensionless form as(6.2)Ns=H2E2SmkR2Tw2=(∂T∂y)2+μBrβ(∂u∂y)2+O(ε2).
In (6.2), the first term can be assigned as N1 and the second term due to viscous dissipation as N2, that is,
(6.3)N1=(∂T∂y)2,N2=μBrβ(∂u∂y)2.
In order to have an idea whether fluid friction dominates over heat transfer irreversibility or vice versa, Bejan [10] defined the irreversibility distribution ratio as Φ=N2/N1.  Heat transfer dominates for 0≤Φ<1 and fluid friction dominates when Φ>1. The contributions of both heat transfer and fluid friction to entropy generation are equal when Φ=1. In many engineering designs and energy optimization problems, the contribution of heat transfer entropy N1 to overall entropy generation rate Ns is needed. As an alternative to irreversibility parameter, the Bejan number (Be) is define mathematically as
(6.4)Be=N1Ns=11+Φ.
Clearly, the Bejan number ranges from 0 to 1. Be=0 is the limit where the irreversibility is dominated by fluid friction effects and Be=1 corresponds to the limit where the irreversibility due to heat transfer by virtue of finite temperature differences dominates. The contributions of both heat transfer and fluid friction to entropy generation are equal when Be=1/2.


## 7. Results and Discussion

We emphasize here that an increase in the parameter value ofβ indicates an increase in the fluid viscosity and a decrease in the fluid activation energy while an increase in the parameter value of γ signifies an increase in the reactive flow Arrhenius kinetics. Table 1 below demonstrates agreement between the results obtained using perturbation technique and purely fourth-order Runge Kutta numerical integration approach coupled with shooting method at small and moderate parameter values. Generally, the difference is of order 10-8.Table 1
Comparison between analytical and numerical results (β=0.1, γ=0.5)
yT(y) Perturbation ResultsT(y) Numerical Results|Tnumer.-Tperturb.|00000.10.023600449430.023600395205.423× 10-80.20.042083800220.042083743845.638× 10-80.30.055355320380.055355252286.810× 10-80.40.063346130330.063346049478.086× 10-80.50.066014408110.066014326978.114× 10-80.60.063346130230.063346049478.076× 10-80.70.055355320370.055355252286.809× 10-80.80.042083800200.042083743845.636× 10-80.90.023600449440.023600395205.424× 10-81.0000The Hermite-Padé approximation procedure in Section4 above was applied to the first few terms of the solution series in Section 2 and we obtained the results as shown in Tables 2 and 3 below.Table 2
Computations showing the criticality procedure rapid convergence (β = 0.1).
dNNuγcN245.08545486714993.9528766995579385.08498312497323.95283121152074135.08498318158073.95283121483905195.08498318156643.95283121483836265.08498318156643.9528312148383Table 3
Computations showing thermal criticality for different parameter values.
β00.10.150.2Nu4.00000000000005.08498318156646.02157317389347.717815638192141γc3.513830719125163.95283121483834.25060382646474.647918009128950The results in Table2 reveal the rapid convergence of Hermite-Padé approximation procedure with gradual increase in the number of series coefficients utilized in the approximants. In Table 3, it is noteworthy that the magnitude of thermal criticality (γc) increases with an increase in the reactive flow activation energy parameter (i.e., β≫0). This invariably will lead to a delay in the development of thermal runaway in the flow system and enhance flow thermal stability.The velocity profiles are reported for increasing values ofγ and β in Figures 2 and 3.  The fluid velocity is zero at the lower stationary plate and increases gradually towards the upper moving upper plate. For  γ=0,  the  fluid  shows  the  standard Couette  linear  velocity  profile  with  maximum  velocity  at  the  moving  upper  plate.  As  the parameter value of γ>0  increases, the Arrhenius kinetic increases, causing the velocity profile  to  increase  nonlinearly  across  the  channel  to  a  maximum  at  y=1.  Furthermore, for increasing value of β, an inflexion point appears in the velocity profile around the centre of the channel as shown in Figure 3.Figure 2
Velocity profile:β=0.3; (——) γ=1; (ooooo) γ=2; (++++) γ=4; (……) γ=5.Figure 3
Velocity profile:γ=5; (——) β=0.3; (ooooo) β=0.5; (++++) β=0.7; (……) β=1.Typical variations of the fluid temperature profiles in the normal direction are shown in Figures4 and 5. The fluid temperature increases with increasing values of γ. This can be attributed to an increase in heat generation within the fluid due to exothermic reaction as illustrated in Figure 4. A decrease in the fluid temperature is observed when the parameter value of β increases; in this case, the fluid activation energy is reduced and its viscosity has increased. Meanwhile, minimum temperature is generally observed at both the lower and the upper plate surfaces while the maximum temperature occurs around the core region of the channel.Figure 4
Temperature profile:β=0.3; (——) γ=1; (ooooo) γ=2; (++++) γ=4; (……) γ=5.Figure 5
Temperature profile:γ=5; (——) β=0.3; (ooooo) β=0.5; (++++) β=0.7; (……) β=1.A slice of the bifurcation diagram for0<β≪1 in the (γ, Nu) plane is shown in Figure 6. It represents the qualitative change in the flow system as the parameter (γ) increases. In particular, for 0≤β≪1 there  is  a  critical  value  γc   (a  turning  point)  such  that,  for  0<γ<γc, there  are  two  solutions (labelled I and II).  The upper and lower solution branches occur due to the temperature-dependent variable viscosity and Arrhenius kinetics in the governing thermal boundary layer equation (2.10). When  γ>γc, the system  has  no  real  solution  and  displays  a  classical  form  indicating thermal runaway. As temperature increases the fluid viscosity decreases exponentially. The velocity gradient  specified by (2.10) increases exponentially with temperature coupled with increasing Arrhenius kinetics and feeds back into the temperature equation, leading to thermal runaway [5, 8, 9].Figure 6
A slice of approximate bifurcation diagram in the (γ, Nu(β=0.1)) plane.Figures7 and 8 display results  for  the  entropy  generation  versus  the  channel  width  for  various  parametric  values. Generally, entropy generation rate is maximum at the plate surfaces and minimum around the core region of the channel. It is interesting to note that the entropy generation rate decreases with increasing value of β and increases with increasing value of γ.Figure 7
Entropy generation rate:γ=2; (——) β=0.5; (ooooo) β=0.6; (++++) β=0.7; (……) β=0.8.Figure 8
Entropy generation rate:β=0.5; (——) γ=2; (ooooo) γ=3; (++++) γ=4; (……) γ=5.Figures9 and 10 display the Bejan (Be) number versus the channel width. It is observed that the fluid friction irreversibility dominates at the channel core region while the heat transfer irreversibility dominates at both the lower and upper plate surfaces. The dominant effect of heat transfer irreversibility at the plate increases with increasing values of β  and γ.Figure 9
Bejan number:γ=2; (——) β=0.1; (ooooo) β=0.3; (++++) β=0.5; (……) β=0.8.Figure 10
Bejan number:β=0.5; (——) γ=2; (ooooo) γ=3; (++++) γ=4; (……) γ=5.


## 8. Conclusion

The evaluation of the entropy production rates for variable viscosity reactive Couette flow was carried out using both analytical and numerical techniques. Solutions are obtained for fluid velocity and temperature profiles. Using a special type of Hermite-Padé approximation technique, we obtain accurately the thermal criticality conditions and the solution branches. The volumetric entropy generation rate and the Bejan number depend on fluid viscosity variation and activation energy parameter (β) and heat generation parameter (γ). Our results reveal that, for all parametric values, fluid friction irreversibility dominates at the channel core region while at both lower fixed and upper moving plate surfaces the heat transfer irreversibility dominates.



---

*Source: 278104-2010-06-22.xml*