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

: 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 num-bers. 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.


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 321 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][9][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][29][30][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][18][19][20][21][22][23][24][25][26][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.

Mathematical Model
The schematic drawing of a two-dimensional lid-driven cavity model within the Cartesian coordinate system is shown in Fig. 1 are the temperature gradients along the fixed wall in x-and y-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. The continuity and governing equations for the two-dimensional in-compressible viscous unsteady Newtonian fluid flow in x-and y-directions are given by [34] ∂u ∂x where t is time and g x = g y = g is 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: left BC : where Δx = Δy = h is the step size, j is the row, i is the column, n is the number of nodes and l s is a slip length. If the slip length, l s is 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 where T is a temperature. Substituting Eqs. (5) and (6) into Eq. (12) will result in The temperature, T at 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 bottom BC : top left corner BC : bottom left corner BC : bottom right corner BC : where λ = kΔt ρch 2 is Fourier number, Δt is a time step and t is the current time step.

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, Δt should be evaluated again so that it can fulfill Eqs. (26) and (27). Thus, the adopted maximum allowable time step, Δt max is given by where min ( ) 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: where n is 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 = t end .
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 −5 is chosen [37].
Step IV: The interior nodes of vorticity, ζ are updated using Eq. (23). The maximum allowable time step, Δt max is 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, t max is utilized in Step IV.
Step VI: The values of temperature at those four walls and four corners namely T t+1 i,j=1 , T t+1 i,j=ny , i=nx,j=ny are calculated using Eq. (14) till Eq. (21).
Step VII: Steps II-VI are repeated until the steady-state condition of the flow is achieved or at any desired time step.

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 of e = 10 −5 ). The applied velocity is U wall = 10 m/s (except U wall = 5 m/s in the 3rd and 4th test cases). The fluid density, ρ is set to 15 kg/m 3 , dynamic viscosity, μ = 2 Pa.s, thermal conductivity, k = 10 W/mK, specific heat, c = 10J/kg·K (except c = 20 J/kg·K in the 5th and 6th cases) and initial temperature, T i of the whole fluid inside the cavity is set to 25 • C. The time step, Δt is 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, Re L can be reckoned via Eq. (33) [38] as where L is height of the cavity. Then, free-slip BC with slip length value, l s = 1 m is assigned in the second, fourth and sixth test cases while the no-slip BC is assigned with the slip length value, l s = 0 m in other test cases. Tab. 1 shows the values of Pr, Re L and α in all these test cases.  Throughout this paper, the 1st and 2nd test cases will serve as the benchmarks (both cases with identical Reynold number, Re L , 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 of u-and v-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 of e = 10 −5 ) in the stream function distribution for the 1st test case (no-slip BC) are shown in Fig. 4. The steady-state criteria of e = 10 −5 [33] and e = 10 −8 are 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 of u-and v-velocity profiles between FDM (MATLAB®) and FEM (Ansys Fluent®) at t = 1.50 s. Fig. 6 shows a good agreement is achieved for u-velocity profile along the centre horizontal line and v-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.  Fig. 7. Besides, the associated contour plots of u-and v-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.   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, Re L than 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, U wall which 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, ψ max and the resultant viscous heating will bring about a reduction in the temperature profile T(ψ 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 of c will 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]:  where δ is thickness of the velocity boundary layer and δ t is thickness of the thermal boundary layer. It is also noticeable that the steady-state period at the criterion of e = 10 −5 for 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, Re L and free-slip BC whereas there is nothing to do with thermal conductivity, k and specific heat, c of 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, T max at the end of the iteration is constant for all the test cases.

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, Re L , 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.

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