Numerical Solutions for Heat Transfer of An Unsteady Cavity with Viscous Heating

H.F.Wong,Muhammad Sohail,Z.Siri and N.F.M.Noor,＊

1Institute of Mathematical Sciences,Faculty of Science,University of Malaya,50603 Kuala Lumpur,Malaysia

2Unit of Engineering Design and Product Development,Technology and Engineering Division,Malaysian Rubber Board,47000 Sungai Buloh,Malaysia

3Department of Applied Mathematics&Statistics,Institute of Space Technology,44000 Islamabad,Pakistan

Abstract:The mechanism of viscous heating of a Newtonian fluid filled inside a cavity under the effect of an external applied force on the top lid is evaluated numerically in this exploration.The investigation is carried out by assuming a two-dimensional laminar in-compressible fluid flow subject to Neumann boundary conditions throughout the numerical iterations in a transient analysis.All the walls of the square cavity are perfectly insulated and the top moving lid produces a constant finite heat flux even though the fluid flow attains the steady-state condition.The objective is to examine the effects of viscous heating in the fully insulated lid-driven cavity under no-slip and free-slip Neumann boundary conditions coupled with variations in Reynolds and Prandtl numbers.The partial differential equations of time-dependent vorticity-stream function and thermal energy are discretized and solved using a self-developed finite difference code in MATLAB® environment.Time dependence of fluid thermodynamics is envisaged through contour and image plots.A commercial simulation software,Ansys Fluent®utilizing a finite element code is employed to verify the finite difference results produced.Although the effect of viscous heating is very minimal,Neumann no-slip and free-slip boundary conditions are able to trap the heat inside the fully insulated cavity as the heat flux is constantly supplied at the top lid.A lower Reynolds number and a greater Prandtl number with free-slip effects reduce temperature distribution in the cavity with a faster velocity than in the no-slip condition as the free-slip behaves as a lubricant.

Keywords: Lid-driven; viscous heating; vorticity-stream function; finite difference method; finite element method

Nomenclature

1 Introduction

Viscous heating is generated by a fluid’s viscous friction when the fluid is sheared by an external motion, provided the fluid is not inviscid and sluggish [1].Viscous heating process is important in creating energy sources while Neumann boundary condition (BC) is important in renewable energy applications such as in energy saving or to increase energy efficiency in buildings.Viscous heating is studied in a slipped flow region using a free-slip boundary condition by Mohamed [2].His study shows a momentous effect on the heat transfer and fluid flow peculiarity in microchannel studies even though the fluid has tenuous viscous flow and possesses a tiny value of hydraulic diameter.Koo et al.[3] discovered that viscous dissipation is effectively influenced by Reynolds number and geometrical aspect ratio.The viscous dissipation effect on free-slip BC is also discussed by Rij et al.[4,5] in a 3D micro-channel gas flow where it is dependent on the amount of rarefaction and other factors.

Another important study of viscous dissipation merged with Joule heating, carried out by Kumar et al.[6], shows an upward trend of temperature distribution in a non-Newtonian fluid flow.The heat transfer process of a Newtonian fluid which takes into account the viscous dissipation in a fully-developed region was performed using an explicit numerical method [7].The variation of the heat transfer process presented an unusual degree in the fully-developed region,both hydro-dynamically and thermally, due to the thermal energy balance between the heated walls is combined with adiabatic walls and the viscous heating [8–10].Other applications of combined Neumann BC and constant wall temperature can be discovered in magnetized nanofluids in a porous cylinder [11], buoyancy effects [12], open trapezoidal cavity [13] and in variation of heat source lengths [14].The catastrophe such as the source of heat in storm or hurricane is derived from the viscous dissipation near the surface [15] and it can be estimated via the turbulent kinetic energy [16] where the heat viscous dissipation is a cubic function of the wind speed.

In some cases, such as in magneto-hydrodynamic (MHD) and mixed convection studies [17,18], the electrically conducting fluid motion is retarded by the magnetic field [19], hence it produces an insignificant viscous heating as compared with the amount of heat generated by a buoyancy force.Generalized differential quadrature method was utilized to study the MHD non-Newtonian nanofluid flow and the impacts of various physical properties [20].Other numerical methods such as Gear–Chebyshev–Gauss–Lobatto collocation method [21] and the fourth-fifthorder Runge–Kutta–Fehlberg method [22] were used to solve an MHD nanofluid flow with a thermal radiation and an applied magnetic field respectively.The study of MHD flow was also conducted in a fluid thermal analysis where three modes of convection specifically referred as free, forced and Marangoni convections are considered [23].Besides, Zhou et al.[24] examined Boussinesq approximation [25] in the thermal energy equation when a Newtonian fluid is subject to a mixed convection process with different heat sources in a lid-driven cavity.Rahman et al.[26]revealed that the velocity and temperature profiles of a lid-driven flow with a semi-circular heated wall are independent of the viscous heating when they are subject to magnetic field and Joule heating.Anderson [27] derived the thermal energy equation for a highly viscous flow which includes the viscous terms used to estimate the meteorite surface temperature.Newtonian homogenous and heterogenous nanofluids [28–31] can also be considered in the future studies as the occupying fluids in the cavity flows with viscous heating and slip effects.

The cited literatures [17–27] obviously show that viscous heating is important in fluid dynamics and thermal analyses where the viscous heating is comparable with other effects.The intention of this work is to utilize the lid-driven cavity model combined with Neumann BCs [32].To make our analysis easier, the moving top lid is not taken into account in the heat transfer analysis thus it will be presumed as an insulator.In order to facilitate the calculations, all the governing equations employed in our paper are given in terms of the vorticity-stream function [33] and thermal energy equation [25].Then we use an explicit finite difference method (FDM) to discretize all the governing equations and conduct the fluid flow and heat transfer analyses from an unsteady-state until the steady-state flow is achieved.Those results are then compared with the solutions of finite element method (FEM) using Ansys Fluent®.We also impose free-slip BC in the cavity model to investigate how the viscous heating process influences the temperature distribution inside the whole cavity.

2 Mathematical Model

The schematic drawing of a two-dimensional lid-driven cavity model within the Cartesian coordinate system is shown in Fig.1.In this case, the external velocity,Uwallis imposed on the top lid,uandvare the velocity vectors inxandydirections,ρandμare fluid density and dynamic viscosity,andare the velocity gradients along the fixed wall inx- andy-directions,HandLare height and length of the cavity.Subsequently,kandcare thermal conductivity and specific heat of the fluid,andare the temperature gradients along the fixed wall inx- andy-directions.The fluid flow process is presumed as twodimensional, in-compressible, laminar and the Newtonian fluid properties are not in the function of temperature or they behave as constants.Then the buoyancy force is excluded from the heat transfer analysis as the heat flux (due to viscous heating) is generated from the moving top lid.Therefore, no convection occurs from the bottom fixed lid to the top moving lid while Boussinesq approximation [25] is not included in the thermal energy equation.The top moving lid is presumed fully insulated.

Figure 1:Schematic drawing of a 2D cavity flow with Neumann BC

The continuity and governing equations for the two-dimensional in-compressible viscous unsteady Newtonian fluid flow inx- andy-directions are given by [34]

wheretis time andgx=gy=gis gravity acceleration.The curl [35] of Eqs.(2) and (3) is known as vorticity-stream function scheme or vorticity transport equation and it is given by [33] as

where stream function,ψis given by Eqs.(5) and (6) while vorticity,ζis given by Eq.(7),

The vorticity,ζat the moving top lid can be obtained via Taylor series expansion [36] and it is given by:

The vorticity,ζat the left, right and bottom walls for the free-slip BC are given by:

where Δx=Δy=his the step size,jis the row,iis the column,nis the number of nodes andlsis a slip length.If the slip length,lsis equal to zero, Eq.(9) till Eq.(11) can be used for vorticity,ζat the left, right and bottom walls for the no-slip BC respectively.

Two-dimensional unsteady thermal energy equation [25] which includes the effect of viscous dissipation is given by

whereTis a temperature.Substituting Eqs.(5) and (6) into Eq.(12) will result in

The temperature,Tat each of the walls and at each of the corners can be obtained by introducing Neumann BC [32] in the unsteady-state condition where they are given by

3 Numerical Method

In this section, some matters related to numerical modelling for solution of the twodimensional vorticity transport equation and thermal energy equation will be discussed.These issues are applications of the simplest numerical method, stability issue and iteration algorithm.Since Eq.(4) and Eq.(13) are nonlinear second-order partial differential equations, an explicit FDM is the simplest method among the feasible techniques for the numerical solution.However,it is conditionally stable [33].Eqs.(4) and (13) can be transformed into equivalent finite difference equations (FDE) using forward and central finite-divided-difference formulas.

Rearranging Eq.(22) yields

Similarly, the resultant FDE for the thermal energy equation is given by

Rearranging Eq.(24) yields

From Eq.(23), the criterion is determined by requiring that the coefficient associated with the node of interest at the previous time is greater than or equal to zero [25], such that

and from Eq.(24),

Eqs.(26) and (27) are the conditional stability criteria.However, the maximum allowable time step, Δtshould be evaluated again so that it can fulfill Eqs.(26) and (27).Thus, the adopted maximum allowable time step, Δtmaxis given by

wheremin( ) means the smallest element is taken inside the bracket.Eq.(28) is called Poisson equation [32].Substituting Eqs.(5) and (6) into Eq.(7) yields

Similarly, Eq.(29) can be discretized using the FDE as follows:

wherenis the number of iterations.Eq.(29) can later be solved by using Gauss–Seidel method [35].The time-marching method to solve the forward FDE of the vorticity transport equation and thermal energy equation are carried out in MATLAB®environment following the procedure below:

Step I:Initial values of stream function,ψand vorticity,ζare set to zero at the first time step,t=0.Values of stream function,ψat the top, bottom, left and right boundaries are set to zero along the iterative process from the 1st time step,t=0 till the last time step,t=tend.

Step II:The interior nodes of the stream function,ψare solved using Eq.(29) by Gauss–Seidel method or any other method until all these values converge.In this paper, the convergence criterion of 10−5is chosen [37].

Step III:The values of vorticity at those four boundaries,ζtop,ζbottom,ζleftandζrightare calculated using Eqs.(8)–(11).

Step IV:The interior nodes of vorticity,ζare updated using Eq.(23).The maximum allowable time step, Δtmaxis ensured to be within the limit prescribed by Eq.(28).

Step V:An initial value of temperature throughout the lid-driven is given.Then Eq.(25)is coupled to the function,ψwith the value obtained from Eq.(30).Similarly, the maximum allowable time step,tmaxis utilized inStep IV.

Step VI:The values of temperature at those four walls and four corners namelyare calculated using Eq.(14) till Eq.(21).

Step VII:Steps II–VIare repeated until the steady-state condition of the flow is achieved or at any desired time step.

4 Results and Analysis

In this study, the aspect ratio,A=1 (1 m×1 m) will be used throughout our paper.There are six case studies considered to characterize the flow field and temperature distribution in the two-dimensional lid-driven cavity flow with insulated walls in order to understand the effect of viscous heating in both no-slip and free-slip Neumann BCs.The six test cases with different fluid properties and laminar unsteady flow conditions are performed in an iterative procedure until the analysis period,t=1.5 s is achieved (longer than the steady-state criterion ofe=10−5).The applied velocity isUwall=10 m/s (exceptUwall=5 m/s in the 3rd and 4th test cases).The fluid density,ρis set to 15 kg/m3, dynamic viscosity,μ= 2 Pa.s, thermal conductivity,k=10 W/mK,specific heat,c=10J/kg·K (exceptc=20 J/kg·K in the 5th and 6th cases) and initial temperature,Tiof the whole fluid inside the cavity is set to 25°C.The time step, Δtis set to 0.001 s and the entire computation is performed on uniform grids 40 × 40 (Δx=Δy=h).The Prandtl number,Pr is defined as

whereνis a kinematic viscosity andαis a thermal diffusivity defined by

The Reynolds number,ReLcan be reckoned via Eq.(33) [38] as

whereLis height of the cavity.Then, free-slip BC with slip length value,ls=1 m is assigned in the second, fourth and sixth test cases while the no-slip BC is assigned with the slip length value,ls=0 m in other test cases.Tab.1 shows the values of Pr,ReLandαin all these test cases.

Table 1:Flow conditions of six test cases

Throughout this paper, the 1st and 2nd test cases will serve as the benchmarks (both cases with identical Reynold number, ReL, but subject to no-slip and free-slip BCs) for comparison with other test cases 3 until 6.The seven steps of FDM as discussed in Section 3 are implemented here and the self-developed FDM code is used to analyse the fluid flow and thermal energy results of the square cavity model in contour and image plots.First and foremost, the 1st test case results are compared with those results of FEM using a commercial engineering software Ansys Fluent®.Figs.2 and 3 show the comparison of contour plots ofu-andv-velocity and temperature distributions respectively produced using the MATLAB®code (FDM) and Ansys Fluent®(FEM).Stream function distribution at the end of iteration and at the location of the maximum stream function value,ψmax(designated as point A, evaluated based on the steady-state criterion ofe=10−5) in the stream function distribution for the 1st test case (no-slip BC) are shown in Fig.4.The steady-state criteria ofe=10−5[33] ande=10−8are shown in the stream function profile of the point A as shown in Fig.5.Centre vertical line and centre horizontal line of the lid-driven model as shown in Fig.1 were used to create the comparison ofu- andv-velocity profiles between FDM (MATLAB®) and FEM (Ansys Fluent®) att=1.50 s.Fig.6 shows a good agreement is achieved foru-velocity profile along the centre horizontal line andv-velocity profile along the centre vertical line at the steady-state condition.Once the FDM MATLAB®code is verified using Ansys Fluent®, the 2nd test case (free-slip BC) is then simulated using FDM.

Figure 2:Comparison of contour plots at t=1.50 s:(a) u-velocity using FDM MATLAB® and(b) FEM Ansys Fluent® (c) v-velocity using FDM MATLAB® and (d) FEM Ansys Fluent®

Figure 3:Comparison of image plot at t=1.50 s (a) using FDM MATLAB® and (b) FEM Ansys

Figure 4:Stream function distribution at t=1.50 s

Figure 5:Stream function profile at the point A

Stream function, temperature distribution and location of the maximum temperature,Tmaxfor the 2nd test case are shown in Fig.7.Besides, the associated contour plots ofu- andv-velocity distributions of the cavity flow are shown in Fig.8.Based on Figs.3a and 7b, the temperature distributions are quite similar in both no-slip (1st test case) and free-slip (2nd test case) BCs.The heat is generated from the moving top lid and concentrated at the top right corner first.Subsequently the heat is slowly propagated to the right bottom corner and to the entire cavity.However, the overall temperature distribution is lower when the lid-driven flow is subject to the free-slip BC in Fig.7b (2nd test case).Formation of eddies in the free-slip BC is also restricted in contrast with the no-slip BC (1st test case) as depicted in Fig.7a.Comparison of temperature profiles at the point A for both test cases is shown in Fig.9.

Figure 6:Comparison of results for (a) the u-velocity profile and (b) the v-velocity profile

Figure 7:(a) Stream function and (b) temperature distribution for the 2nd test case (free-slip BC)

Figure 8:(a) U-velocity and (b) V-velocity distributions for the 2nd test case (free-slip BC)

Figure 9:Temperature profile at the point A:(a) 1st test case (no-slip) and (b) 2nd test case (freeslip)

Both slip BCs show similar trend in temperature profile at the point A (location of maximum stream function,ψmaxas shown in Figs.3 and 7) and the temperatures keep increasing, but the no-slip BC gives an upper offset than the free-slip BC as shown in Fig.9.This is due to the free-slip of the fluid at the wall or solid interface, thus no frictional heat can be generated at the wall.These two benchmarks (1st and 2nd test cases) are then compared with the 3rd and 4th test cases (with lower Reynold number,ReLthan those benchmarks) as shown in Fig.10.The 3rd and 4th cases show a reduction in temperature profile due to the lower value of applied velocity,Uwallwhich leads to a smaller stream function distribution in the whole lid-driven flow as prescribed by

Eqs.(5) and (6).Because of the last term on the right-hand side of thermal energy equation in Eq.(12), the viscous heating term is much related to the stream function,ψmaxand the resultant viscous heating will bring about a reduction in the temperature profileT(ψmax)at the point A.The benchmarks (the 1st and 2nd test cases) are then compared with the 5th and 6th test cases(with higher Prandtl number, Pr than those benchmarks) as shown in Fig.11 for the temperature profile at the point A.All the results of the 1st until the 6th test cases are shown in Tab.2.

Both the test cases (no-slip and free-slip BCs) show lower temperature profiles at the point A as compared to the benchmarks.Since Pr is directly proportional to specific heat,c, a higher value ofcwill result in a reduction of thermal boundary layer thickness.Recall that Prandtl number, Pr is a ratio of momentum and energy transport by diffusion in the velocity and thermal boundary layers given by [25]:

Figure 10:Temperature profile at the point A (ReL=100 vs.ReL=50)

Figure 11:Temperature profile at the point A (Pr=1 vs.Pr=2)

whereδis thickness of the velocity boundary layer andδtis thickness of the thermal boundary layer.It is also noticeable that the steady-state period at the criterion ofe=10−5for the 5th and 6th test cases are similar with the benchmarks.Although the fluid flow has attained its steady-state condition as elucidated in Fig.5, no steady-state condition can be observed in the temperature profiles of all the test cases as shown in Figs.10 and 11.Moreover, temperature profile increases in quasi-linear behaviour even though Eq.(13) is an elliptic equation [32] due to the moving top lid provides an additional finite heat flux to the entire cavity.It can also be deduced that the temperature profile at the point A (or even at any point inside the cavity) will keep increasing if the iteration is extended at any period.The 1st test case is executed at the unsteady-state condition for 3 s as shown in Fig.12 and it proves that the temperature profile is increasing quasi-linearly.This temperature rise is due to each of the fluid layers is being sheared relative to each other and this will create viscous heating caused by the fluid viscosity.The heat is being trapped inside the cavity flow by the insulated walls with no way to escape.Thus, temperature profile at every point inside the cavity flow (or temperature distribution of the cavity flow) will keep increasing.Based on Tab.2, location of the maximum stream function,ψmax(point A) at the end of the iteration is merely influenced by Reynold number, ReLand free-slip BC whereas there is nothing to do with thermal conductivity,kand specific heat,cof the fluid.This is because solution of the stream function,ψfrom Eq.(23) can be done without consideration of Eq.(25).Also, the location of the maximum temperature,Tmaxat the end of the iteration is constant for all the test cases.

Table 2:Results of six test cases

Figure 12:Temperature profile at the point A (1st test case at prolonged period)

5 Conclusion

In our study, we assume the two-dimensional unsteady lid-driven cavity model is perfectly insulated and no heat loss from the heat flux (viscous heating) is produced by the moving top lid to the environment.The FDM MATLAB solutions are verified earlier with FEM Ansys Fluent®before the six test cases are examined based on no-slip and free-slip BCs.The main findings are as follows:

• Temperature profile shows an increasing trend in all test cases as the produced constant heat flux is supplied continuously by the moving top lid.If the iteration period is prolonged to any other period, the temperature profiles of these test cases will keep increasing.This is the reason why the surface temperature of a meteorite during its free-fall keeps increasing although the frictional heating or viscous heating due to the atmospheric air contributes a very small portion of heat.

• Neumann BC is able to trap the heat produced inside the cavity model as the temperature profiles keep increasing quasi-linearly.Trapped heat is important for energy efficiency such as for building insulation during winter season.

• The lower value of Reynolds number, ReL, the lower value of thermal diffusivity,αand the free-slip effect of the fluid flow can create lower temperature distributions.

• Greater value of Prandtl number, Pr can reduce the temperature distribution inside the cavity model.

• Free-slip effect will cause the fluid to flow faster but with a lower temperature distribution than that of the no-slip BC as the free-slip effect is assumed to act like a lubricant.

• Although FDM scheme is simple and sufficient to explain the phenomenon of Neumann BC and free-slip effect, this method is limited to simple geometry only.

Funding Statement:The authors acknowledged the funding received from the Ministry of Higher Education, Malaysia and University of Malaya (https://umresearch.um.edu.my/) under the Project No:IIRG006C-19IISS leaded by Z.Siri for this study.

Conficts of Interest:The authors declare no conflicts of interest to report regarding the present study.