A novel method for forced convection heat transfer in 2D skewed cavities using a non-orthogonal mesh

Abstract

Keywords: Cavity; Convection; Lid-Driven; Non-Orthogonal; Skewed; Tensor Calculus

Introduction

Lid-driven flow in a rectangular cavity is a benchmark problem in the field of computational fluid dynamics (CFD). High resolution results of the velocity distribution for this problem have been reported in various previous studies [1, 2]. An extension of this problem is the lid-driven flow in a skewed cavity [3]. This problem can be used as a benchmark for testing CFD codes that are written for flow simulation in complex, two-dimensional geometries. Many researchers have studied cavity flow and a thorough list of

the early studies can be found in [3]. Benchmark results are also documented in [3] and [4]. Mansour et al. [5] obtained accurate results using the finite difference method for the problem of mixed convection in a lid-driven square cavity that is partially heated from below, using a water-based nano-fluid. They explored the influence of Reynolds number, Nusselt number, and heat source location on streamlines and isotherm contours [5]. Taher et al. [6] investigated the cavity flow problem by considering partial heating applied to the bottom wall, with insulated top wall and with the opposite walls at various uniform temperature values.

*Corresponding author. *E-mail address: saadbin@kfupm.edu.sa This paper was recommended for publication in revised form by Editor-in-Chief Ahmet Selim Dalkılıç Published by Yıldız Technical University Press, İstanbul, Turkey Yıldız Technical University. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).

Abbasian et al. [7] approached this problem by introducing sinusoidal heating on the sidewalls, maintaining insulation on both the top and bottom walls. Kamyar et al. [8] conducted a comprehensive literature review and undertook a collective analysis of diverse forms of heat convection (natural, forced, or mixed) across various geometries that utilized nanofluids. It is to be noted that the two-dimensional, steady-state, lid-driven cavity flow for large Reynolds numbers evolves to a vortex with an inviscid core surrounded by viscous boundary layers [9]. Hence, for low Reynolds numbers, as is the case in this study, the authors have assumed that the boundary layer effect is not important. In the current study, a new approach is introduced to address the problem of a lid-driven flow in a skewed cavity while employing the primitive variables, i.e. the velocity components and the pressure. The main objective of the paper is to derive, present and numerically solve the Navier-Stokes Equations (i.e. the collection comprising the continuity equation, the momentum equations, the energy equation and the pressure Poisson equation) in a non-orthogonal coordinate system. The various desirable features of the formalism employed in the present study and their advantages are listed below: 1) Structured grid

connectivity matrix is not required; results in a sparse, banded matrix.

boundaries are modelled accurately; boundary conditions are implemented easily.

non-zero boundary velocities are not a linear combination of the Cartesian velocity components and are easily implemented.

primitive variables are employed; staggered grids are not required.

As is well known, utilizing structured grids instead of unstructured grids, results in general to a more efficient solution algorithm and this is because the connectivity information is built into a structured gird, whereas this information has to be explicitly supplied in the case of an unstructured grid. For irregular, complex geometries, structured grids can be constructed by means of the procedure of Numerical Grid Generation [10]. This procedure in essentially introduces a new, non-orthogonal, curvilinear coordinate system. However, it is then necessary to derive the balance equations in the new coordinate system. Usually, this implies re-writing of the balance equations in terms of the new curvilinear coordinates (ξ,η). However, in the case of the Navier-Stokes Equations, this also implies that the momentum balance equations be written in terms of the new velocity components as well as the direction of the momentum equation should be along the tangent to the new basis vectors. Hence, three different transformations

have to be carried out. To the best of the knowledge of the authors, the complete set of these three transformations have not been attempted before in a numerical study. Moreover, appropriate mathematical tool for realizing these transformations is that of Tensor Calculus which again appears to be a novel aspect of this study. A familiar problem in Computational Fluid Dynamics is that of the pressure-velocity coupling. Customarily, this is achieved by means of using three (in 2D) staggered grids as well as the SIMPLE algorithm. Another, infrequently used procedure of pressure-velocity coupling is to replace the continuity equation by the Pressure Poisson Equation (PPE). More detail about this equation is given in the next section. The benefit of using this approach is that a single grid is sufficient for all variables and three or more staggered grids are not required. This is the approach utilized in this paper. A drawback of this approach is that the boundary conditions for the PPE are all of the Neumann type and hence an integral constraint on the pressure variable has to be satisfied explicitly. It is hoped that this approach of utilizing non-orthogonal, structured grids in the numerical solution of the Navier-Stokes Equations will be of assistance in the efficient solution of problems involving geometries in which it may be difficult to generate an orthogonal structured grid and therefore one is compelled to resort to using unstructured grids.

Mathematical Formulation

Description of the Problem The problem investigated in the present study is two-dimensional, nevertheless the methodology employed is general and can possibly be extended to unsteady and three-dimensional flows. The schematic of the problem is shown in Figure 1. The fluid in contact with the all the boundaries have zero tangential and normal velocities with the exception of boundary B3 where the fluid has a tangential velocity component of U (m/s). Furthermore, the boundaries B2 and B4 are insulated whereas the temperatures of boundaries B1 and B3 are maintained at uniform but unequal values. Mathematical Modelling A skew cavity is innately represented by an oblique coordinate system and using such a system considerably simplifies the code as well as the implementation of the boundary conditions. The Cartesian system is represented by the (x,y) coordinates while the oblique system is represented by (ξ,η) coordinates. The skew cavity in the (x,y) coordinates (physical domain) transforms in to a rectangle in the (ξ,η) coordinates (computational domain) under the following transformation equations, (1)

Figure 1. Schematic and coordinate system for the skew cavity with inclination angle a.

From many aspects, it is desirable to numerically solve the Navier-Stokes equations (NSE) in the (ξ,η) coordinate system. To this end, it is necessary to write the NSE in terms of the new coordinates. In the current study, this is achieved by means of the mathematical technique of tensor calculus. This technique provides the most efficient path to deriving the NSE in any general curvilinear coordinate system. We note here that since the NSE is a vector, partial differential equation (PDE) (a set of multiple, scalar PDEs) with the dependent variable also a vector, therefore, in writing the NSE in a general curvilinear coordinate system (of which the oblique coordinates are a special case), three separate transformations have to be carried out. These are, 1) Transform the independent coordinates: (x,y) → (ξ,η). 2) Transform the velocity components: {Vx,Vy} → {V(1),V(2)} 3) Transform the set of two momentum equations that are in the {î,Ĵ} direction into a set of two equations that are in the {ê(1), ê(2)} direction. {V(1),V(2)} are the (physical) velocity components, whereas {ê(1), ê(2)}are the (physical) basis vectors in the (ξ,η) coordinate system. Additionally, the pressure-velocity coupling is achieved by means of the Pressure Poisson Equation (PPE). This allows us to avoid using the staggered grid as well as the SIMPLE algorithm. In fact, a single grid is enough for all the variables. Having a single grid instead of three separate ones and avoiding the SIMPLE algorithm results in a more computationally inexpensive and efficient procedure. Moreover, a finite difference technique is employed instead of the finite volume method. In addition, the energy equation is also transformed to the new (ξ,η) coordinate system. Based on equations (1), the first step is to write the position vector r. Various other entities are derived from it and will be described as they appear. The details can be found in [11] and [12]. Note that summation and range conventions are used in some equations. Moreover, in the index

notation, as used in tensor calculus, we denote (x,y) → (y1, y2) and (ξ,η) → (x1, x2) [11]. The terminology used for various variables primarily follows that of [11]. (2) (3)

(8) The velocity vector V of the fluid in the Cartesian and the oblique coordinate systems is now explained. {Vx,Vy} are the velocity components in the Cartesian system whereas {V1,V2}, {V1,V2} and {V(1),V(2)} are respectively the covariant, contravariant and physical components of the velocity

vector in the oblique (ξ,η) coordinate system. In the oblique coordinates, the velocity vector is written in various ways but ultimately, we are interested in the physical components of the velocity vector.

simplification, the two momentum equations in the ê(1) (ξ-direction) and the ê(2) (η-direction) are written as, (17a) (17b)

(9) The physical components of the velocity vector can also be defined directly as [12],

For the determination of the pressure variable, an alternative procedure is to replace the Continuity Equation by the Pressure-Poisson Equation (PPE). This is accomplished by taking the divergence of the Momentum Equations and then utilizing the Continuity Equation for simplification, where ever it is applicable [14]. These equations are presented below. (18)

(10) Let u = V(1) and v = V(2) then we can finally write, (11) In deriving the various equations below, we require the use of the concept of the covariant derivative from tensor calculus. The subscript |i denotes covariant differentiation with respect to the coordinate xi. The Continuity Equation in a general curvilinear coordinate system is written as [13],

(19) The steady-state, incompressible energy equation with negligible viscous dissipation and constant thermal conductivity is written as [15], (20) Or, in component-form, it is written as, (21)

(12) Introducing the physical components of the velocity vector and simplifying we get, (13) The Navier-Stokes Equation (NSE) in a general curvilinear coordinate system, in vector form, is written as [13], (14) Or, in component-form, it is written as, (15) Or,

In equations (15) and (16), ρ, p and µ are the fluid density, pressure and absolute viscosity, respectively. After

After simplification, we write, (22) Equations (17), (19) and (22) require boundary conditions for their solution. The velocity boundary conditions are straightforward. These are, (23) For the PPE we require a boundary condition at each of the four boundaries B1 B2, B3 and B4 (Fig. 1). These boundary conditions are obtained by applying the normal component of the vector momentum equation at the boundary in question. This is termed as the Neumann Boundary Condition for pressure [16]. At a point on the boundaries B1 and B3, the unit normal vectors are and respectively. At a point on the boundaries B2 and B4, the unit normal vectors are and and respectively. So, for example, at a generic point on boundary B2, the normal component of the momentum equation is determined by taking the inner product of the vector momentum equation with the unit normal vector .

Referring to the computational domain in Figure 1, the ξ-axis is discretised in to Nξ number of divisions with the first and last nodes designated as 1 and Nξ+1 and the η-axis is discretised in to Nη number of divisions with the first and last nodes designated as 1 and Nη+1. Equations (17) are non-linear, partial-differential-equations. Replacing the derivatives in these equations by their appropriate finite-difference approximations will result in a system of non-linear, simultaneous, algebraic equations. These can be solved, for e.g., by means of the NewtonRaphson method. Alternately, the technique of iterative update along with under-relaxation can be utilized and this is the approach taken in our investigation. Equations (17) are first ‘linearized’ by replacing the u and v velocity component factors by and in the convective terms as described below,

The side walls are insulated which at the set of points (0,η) and (1,η) results in,

(24) After all necessary simplifications we get, (25) Moreover, it is of utmost importance that the continuity equation be also satisfied at the boundary [16]. At boundary B2 this implies that, (26) For the energy equation, the boundary conditions at the top and bottom wall are,

(28) All the pressure boundary conditions, the simplified continuity equation and the temperature boundary conditions at the four boundaries are presented in Table 1. Finally, we note that for α = 90°, all the relevant equations, velocity components and the independent variables revert to those that are applicable in a Cartesian system. It is now required to numerically solve equations (17), (19) and (22) simultaneously, subject to the boundary conditions described in equations (23), (29) and (30). This task is accomplished by means of the finite difference method. The details are described in the next section. Numerical Approach In the current study, we have selected the finite-difference method for the numerical solution of the relevant equations due to the simplicity of implementation. The finite-volume and the finite element methods can also be used for the task; however, these methods are more helpful in situations where unstructured grids are being used and for cases in which the solution domain is multiply connected.

(31b) When equations (31) are solved iteratively, the and velocity component factors are kept fixed in the current global iteration and are updated immediately with the newly obtained values of the u and v velocity components at the end of the inner iterations (to be discussed later). In equations (31), the diffusion terms are discretized by means of the second-order, central difference formulae. The cross-derivative terms are treated as a source term during the calculations. The pressure derivatives are also discretized by means of the second-order, central difference formulae. In order to obtain stable discretisation for the convective terms, second-order upwind (also called Linear Upwind) discretisation is used. This entails using second-order backward difference formulae for positive flow direction and second-order forward difference formulae for negative flow direction, for the derivatives ∂u⁄∂ξ, ∂u⁄∂η, ∂v⁄∂ξ, ∂v⁄∂η. It was found that for skewed cavities, first-order upwind discretisation is not satisfactory. In equation (19) also, all the derivatives are approximated by means of

Table 1. Pressure/Temperature B.C. and the continuity equation at the four boundaries Boundaries

the second-order, central difference formulae. The various discretised equations are presented below, (32a) (32b) where, (33a)

equations. This system is then solved iteratively by means of the Gauss-Seidel iterative method (constituting the inner iterations). The iterations start by assuming some suitable u- and v- velocity distributions and p-distribution inside the domain and then solving the system of equations to obtain the correct distributions. However, since the coefficients in equations (32) and the source term Sp are themselves functions of the velocity components, therefore the solution obtained is not correct. The coefficients are updated in the next global iteration by means of the new velocity components and the procedure is repeated. The previously described solution procedure is stabilized by employing the technique of under-relaxation. This is achieved by performing only one or a few Gauss-Seidel iterations in the current global iteration and then updating the coefficients in equations (32). However, to get a converged solution, equations (32) are slightly modified to incorporate under-relaxation as described below. It is not required to introduce under-relaxation in equation (36). (38a)

(38b) where ω is the under-relaxation factor and is taken to be equal to 0.3. Application of the velocity boundary conditions is straight-forward, however, the pressure boundary condition implementation needs discussion. Let us focus at a generic point on boundary B2, (ξ,η) = (1,η). Discretising the pressure boundary condition at that point, we get, (39) (40)

It is pertinent to mention that the right-hand-side (RHS) of equation (19) involves derivatives of the velocity components that originate from the convective terms in the momentum equation. Hence, it is also necessary to use upwind differencing for these derivatives. The source term in Sp , equation (37) incorporates such differencing. Applying equations (32) and (36) at each interior node results in a system of (pseudo) linear, simultaneous

In discretizing equation (40) we have utilized the u-velocity values at the fictitious nodes (j,Nξ +2). These values are obtained through satisfying the continuity equation at the boundary, equation (39). In equation (40), we may observe that the pressure at any given boundary node is related to the pressure at its neighboring boundary nodes. Consequently, we need to solve a linear system of equations at every boundary to determine the pressure. Once again, we can use the Gauss-Seidel iterations to calculate the pressure. However, it is not necessary to get a converged solution within every global iteration, instead a few iterations are enough. It is also important here to indicate that because of utilizing the linear (second-order) upwind scheme, at any interior node that is a neighbor of a boundary node, equations (34), (35) and (37) may require values of the velocity variables that are outside the domain.

This situation is handled by utilizing a first-order upwind scheme at all such nodes. Finally, we mention that the boundary conditions for the pressure Poisson equation (19) are of the Neumann type and this implies that an integral constraint has to be satisfied if equation (19) has to have a solution. The mathematical form of the integral constraint in the oblique (ξ,η) coordinate system is described below.

(22) in a fashion similar to that of the momentum equations (17). These equations are given below, (43) where, (44)

(41) (45) Finally, we note that the treatment of the insulated boundary condition of the energy equation for the determination of the unknown temperature at boundaries B2 and B4 is identical to that of the pressure boundary conditions. (42)

In equation (42) above, the left-hand-side (LHS) should be equal to the (RHS), for the solution to exist for equation (19), having boundary conditions given by equations (29) and (30). During the numerical solution process this may not happen and therefore, before the pressure is updated through equation (36), it is required to modify the RHS of equation (19). This is done by adding the quantity described below to the RHS of equation (19) at all interior nodes.

where A is the total area of the skew cavity. In evaluating the second-order derivatives of the velocity components at the boundaries, it is necessary that the continuity equation is satisfied. Once the velocity and pressure fields have been computed, we may then proceed to the determination of the temperature field. To this end we need to discretize equation

Results And Discussion

In the present study, the continuity equation (13), the momentum equations (17), the pressure poisson equation (19) and the energy equation (22) were solved for a skew cavity with equal sides for Re = 100, Pr ≈ 0.69, Ra = 533, and for skew angle a = 60o, 75o, 90o, 105o and 120o. To have negligible free convection effects it was necessary to keep Ra number less than 1000 [17]. The final grid size used had Nξ × Nη = 256 × 256 number of nodes for the case study. Grid independence test showed that a grid size of Nξ × Nη = 128 × 128 was also acceptable but insufficient to resolve finer flow structures. The mathematical formalism was validated by comparing the relevant results with those published in [4] and with those obtained by solving the same problem through the Ansys 2023 R2 commercial software. Considerably good accuracy was obtained for the cases investigated.

Figure 2 shows the grid independence test results. The NSE were solved for a 60o skewed cavity for grids of sizes 32 × 32, 64 × 64, 96 × 96, 128 × 128, 196 × 196 and 256 × 256. Three variable are plotted, these are: the mean of the absolute values of the normalized u-velocity throughout the domain; the mean of the absolute

values of the normalized v-velocity throughout the domain; the maximum of the absolute values of the normalized v-velocity , where U is the lid speed. It is observed that the test variables become effectively independent of grid size for grids of size greater or equal to Nξ × Nη = 128 × 128.

Figure 3. Normalized u-velocity distribution along the line a-b.

Figure 4. Normalized v-velocity distribution along the line c-d.

Figure 3 shows the variation of the normalized u-velocity distribution along the line a-b with the normalized grid numbers. It is seen that the velocity is negative along some portion of the depth and positive along the remaining portion. For all skew angles, the velocity is negative for a larger portion of the depth and is zero at around the 0.75 mark except for the 120o case for which it is zero at around the 0.65 mark. The results are compared with those published earlier in [4] as well as with the simulation performed in Ansys, utilizing the finite volume method (FVM). It is observed that for all relevant skew angles, the results obtained in the present study are in close or acceptable agreement to those obtained in [4] as well as with the simulation results. Figure 4 shows the variation of the normalized v-velocity distribution along the line c-d with the normalized grid numbers. It is observed that the velocity is positive along some portion of the breadth and negative along the remaining portion. For all skew angles, the velocity is negative for a larger portion of the breadth and is zero at around the 0.6 mark except for the 105o case and the 120o case for which it is zero at around the 0.55 and 0.5 marks respectively. The results are compared with those published earlier in [4] as well as with those obtained through the Ansys software. It is observed that for all relevant skew angles the agreement either excellent or acceptable. Figure 5 shows (a coarse version) of the vector plot for the normalized velocity components. The central vortex can be clearly observed whereas some other flow structures are also hinted at. It is observed that the central vortex

maintains the same relative position with respect to the boundaries. Figures 6 and 7 display the variation of the normalized temperature distribution Θ = (T - Tmin)⁄(Tmax-Tmin) along lines a-b and c-d with the normalized grid numbers for the case of the bottom surface at Tmax. In Figure 6, along line a-b, the temperature distribution can be seen to have a point of inflection which corresponds to the simultaneous effect of convection and diffusion of heat. In Figure 7, the slope of the temperature is zero or close to it, at the left and right boundaries which corresponds to the insulated boundary condition. The results are compared with those obtained from the Ansys software and it is observed that for all relevant skew angles the agreement is either close or acceptable. Figures 8 and 9 display the variation of the normalized temperature distribution Θ = (T - Tmin)⁄(Tmax-Tmin) along lines a-b and c-d with the normalized grid numbers for the case of the top surface at Tmax. In Figure 8, along line a-b, the temperature distribution can again be observed to have a point of inflection which corresponds to the simultaneous effect of convection and diffusion of heat. In Figure 9, the slope of the temperature is zero or close to it, at the left and right boundaries which corresponds to the insulated boundary condition. The results are compared with those obtained from the Ansys software and it is observed that for all relevant skew angles the agreement is either close or acceptable.

Figure 5. Velocity vector plots for cavities of various skew angles.

Figure 6. Normalized temperature distribution along the line a-b for bottom heating source.

Figure 7. Normalized temperature distribution along the line c-d for bottom heating source.

Figure 8. Normalized temperature distribution along the line a-b for top heating source.

Figure 9. Normalized temperature distribution along the line c-d for top heating source.

Conclusion

more flexibility than strictly orthogonal coordinates and are

The present study explores the possibility of the numerical solution of Navier-Stokes Equations (NSE) in non-orthogonal coordinate systems. Such coordinate systems offer

believed to be able to successfully describe more complex regions. Moreover, they result in structured grids which are in general computationally less expensive than the unstructured

grids, which are customarily utilized in the finite-volume and the finite-element methods. The correct and complete derivation of the NSE that is applicable in non-orthogonal coordinate systems, involves the change of independent variables, the velocity components and the direction in which the momentum equations are valid. These derivations can be efficiently performed by means of the mathematical tool of Tensor Calculus and our paper describes them in detail. The mathematical formalism so obtained is validated by numerically solving it for a bench mark problem of the skewed lid-driven cavity. The results acquired from these numerical studies are compared with those published earlier and with those obtained through an independent Ansys simulation of a similar problem. It was found that the present formulation produces results that are in good agreement with the earlier studies for relevant skew angles. In general, the approach presented in this paper is more efficient than those based on the finite volume (FV) and finite element (FE) formulations. This is because the customary FV & the FE methods employ unstructured grids that do not have the connectivity information built in to them. Hence, connectivity information is to be provided separately and this contributes to the computational overhead. Moreover, unstructured meshes results in matrices that may be sparse but that are not banded. This means that efficient matrix inversion algorithms, such as the tridiagonal matrix algorithm, cannot be employed. It is surmised that these considerations make the present algorithm more competitive. Finally, it is believed that the algorithm can be extended to three-dimensional as well as unsteady problems without any fundamental difficulties. So, for example, for an unsteady, incompressible flow, equations (13) & (19) do not change. However, a time-derivative term will have to be included on the LHS in equations (17) based on the term ρ ∂V⁄∂t that is included on the LHS of equation (14). Similarly, the term ρcp ∂T⁄∂t is to be included on the LHS of equation (20) and based on it, an appropriate term is to be included in equation (22). For three-dimensional problems we note that equations (12), (16), (18) & (21) are actually not restricted to two-dimensions but are valid in an N dimensional space. Specifically, they can be written for a three-dimensional space.

Nomenclature And Abbreviations

Total area of the skew cavity (m2) Heat capacity (J/kg/K) Physical basis vector (unit vector) in the ξ-direction Physical basis vector (unit vector) in the η-direction Covariant bases vectors Metric tensor Conjugate metric tensor

î Ĵ k Nξ Nh p Pr r Ra Re T U V (V1,V2) (V1,V2) u, Vx v, Vy V(1) V(2) x y α ξ η ρ µ ω Θ

Connection coefficients unit vector in the x-direction unit vector in the y-direction Thermal conductivity (W/m/K) number of divisions along the ξ-axis number of divisions along the η-axis Pressure (Pa) Prandtl number Position vector Rayleigh number Reynolds number Temperature (K) Lid speed (m/s) Velocity vector of a fluid particle (m/s) Covariant components of the velocity vector Contravariant components of the velocity vector x-component of the velocity vector (m/s) y-component of the velocity vector (m/s) ξ-component of the velocity vector (m/s) η-component of the velocity vector (m/s) Cartesian coordinate (m) Cartesian coordinate (m) Skew angle Curvilinear coordinate Curvilinear coordinate Fluid density (kg/m3) Absolute viscosity (Pa.s) Under-relaxation factor Normalized temperature

FVM Finite-volume method NSE Navier-Stokes Equations PDE Partial Differential Equation PPE Pressure Poisson Equation SIMPLE Semi-Implicit Pressure Linked Equations

Acknowledgement

The authors would like to acknowledge the support of King Fahd University of Petroleum and Minerals (KFUPM), Saudi Arabia. Acknowledgement is also extended to the IRC for Sustainable Energy Systems (IRC-SES), King Fahd University of Petroleum and Minerals, Dhahran, Saudi Arabia.

Data Availability Statement

The authors confirm that the data that supports the findings of this study are available within the article. Raw data that support the finding of this study are available from the corresponding author, upon reasonable request.

Conflict Of Interest

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Ethics

There are no ethical issues with the publication of this manuscript.

Share and Cite