Variant slip effects, aspect ratios and steady-state criteria in unsteady cavities with direct PDE simulation of eddies

The importance of fluid flow study can be visualized in the influence of flow by sea buckthorn before the dam [1]. Particle Image Velocimetry (PIV) is a sophisticated technique in flow visualization in many engineering and environmental applications. Another application of fluid flow in nuclear engineering includes coolant leaking detection in flow pattern that is important in the design of heat exchangers. The study of free-slip boundary condition (BC) is crucial at certain occasions in order to imitate or to access its impact on the whole fluid flow pattern especially at a very low pressure [2] or when immersion or penetration of bubble or fluid particle into solid interface happens and results as surface contamination. Thus, the large viscosity ratio of the fluid near the solid boundary happens and causes in the rise of tangential velocity near the solid boundary [3]. In some cases, no-slip BC is well-posed in modeling of viscous flow without penetration condition but it is amiss when analysis of highly inviscid flows is considered. Free-slip is important to simulate free surface applications in geophysical occurrences such as in river and glacier, spilling dynamics, ship hull designs, technologies of coating and fuel spraying or injection. Free-slip is measured by a so-called slip length physical parameter [4] and it can be measured by specific equipment. In this study, free-slip BC is applied on the classical problem of lid-driven cavity flow with different aspect ratio (AR) to see the impacting results. This model is a two-dimensional unsteady cavity with a Newtonian fluid is flowing inside it and subject to a translation of top lid at a constant velocity [5].

Swirling motion or eddy can happen in an unsteady flow at low velocity under certain geometries. Eddies in ocean caused nutrients found in deep water to rise to the flow surface while eddies in enclosures create more dynamics patterns of the fluid inside that entire enclosures. Eddy formation helps one to understand specific behavior of fluid dynamics better especially to determine the steady-state condition of the flow. According to Moffatt [6], two corner eddies at the bottom merge to form a secondary main cell at critical AR of 1.629 in no-slip BC besides a primary main cell which is formed at the beginning of the unsteady-state simulation. Validating our present work as compared to Moffat’s, we found that the primary main cell is rotating in clockwise direction while the two corner eddies at the bottom are rotating in counter-clockwise direction with infinitesimal negative stream function value even after they have formed a secondary main cell (negative stream function value means rotating in counter-clockwise direction). However, no study was carried out to see if the above study is subject to a significant free-slip BC. Venturing further, we found that a secondary main cell is formed when iteration is terminated after a steady-state condition with criterion of e ≤ 10⁻⁸ is achieved in no-slip BC.

Normally, for a Newtonian fluid analysis, there are two important dimension less parameters which are often used often namely, Reynolds number, Re=ρULμ {\mathop{\rm Re}\nolimits} \; = \left( {{{\rho UL} \over \mu }} \right) and the aspect ratio, A=HL {\rm{A}} = \left( {{H \over L}} \right) of the cavity [7]. The model contemplated here is known as “singular driven cavity” as there are two discontinuities in the boundary condition at the top lid [8]. In order to avoid the troubles connected with primitive variable approach (velocity-pressure), vorticity-stream function and vorticity-velocity formulations of the Navier-Stokes Equations (NSEs) are broadly employed [9]. Thus, the numerical scheme which employs vorticity-stream function or vorticity transport equation is used to solve the two-dimensional NSEs inside the lid-driven cavity flows. To our best knowledge, our work is the first attempt to carry out the analysis for both no-slip and free-slip effects which are never considered before by the previous authors as discussed in [6, 9,10,11]. Free-slip BC is important for subsequent analysis such as in heat transfer applications because it will influence the temperature distribution of the fluid. The application of free-slip effect can be incorporated with nanofluids, electro-magnetohydrodynamics (EMHD) and magnetohydrodynamics (MHD) to study their flow patterns where these can be compared with the no-slip effect. In addition, free-slip effect should be included in hydrodynamics analysis as long as the temperature profile along the interaction between fluid and solid surface can be influenced. This influence may change the heat transfer process along the boundary layer through convection and radiation especially for thermal dissipated viscous fluids.

This paper outlines the importance of free-slip BC and demonstrates its effect on lid-driven cavity flow model with AR=1.629. In Section 2, the mathematical model of the problem using vorticity-stream functions is applied to both no-slip and free-slip BCs. Then, the finite difference method (FDM) is chosen and employed to solve the vorticity-stream function using a self-developed Matlab_® algorithm which will be verified for the square model (AR = 1) with no-slip BC in Section 3. Solution for both mentioned BCs will also be discussed. These steady-state results are then compared with finite element method (FEM) using Ansys Fluent_®. No-slip BC and different steady-state criteria are later applied on the lid-driven cavity flow model with AR = 1.629 and in free-slip BC with their flow patterns are compared intensively in Section 4. The main novelty of this work is observation on the eddies formation under both no-slip and free-slip BCs with various steady state criteria and aspect ratios for both square and rectangular cavities which altogether cannot be found in the existing studies. FDM method was directly employed to solve the unsteady partial differential equations without conversion into nonlinear ordinary differential equations and some alternate solutions using FEM in Ansys Fluent_® are provided. It is also worth to mention that the previous works only focused on solving such cavity problems but none of them developed their own GUI to help users in solving and analyzing the unsteady cavity lid-driven problem for both no-slip and free-slip BCs. Besides, our user-interactive GUI includes notification of the maximum allowable step time based on the stability criterion which is bounded by the explicit FDM. In Section 5, results and discussion are presented. Finally, Section 6 introduces the paper by marking main contribution of this paper.

The lid-driven model with no-slip condition is shown in Figure 1, where U is the applied velocity, u and v are the velocity vector in x and y directions, ρ and μ are fluid density and dynamic viscosity, ∂v∂x {{\partial v} \over {\partial x}} and ∂u∂y {{\partial u} \over {\partial y}} are the velocity gradient along the no-slip BC in x and y directions, H and L are the height and width of the cavity. The governing equations of an in-compressible laminar viscous unsteady-state Newtonian fluid flowing through a Cartesian plane constituted by a continuity equation and Conservation of Momentum in two-dimensional problem are given by: (1) continuity equation: ∂u∂x+∂u∂y=0, {\rm{continuity}}\;{\rm{equation}}:\,\,\,{{\partial u} \over {\partial x}} + {{\partial u} \over {\partial y}} = 0, (2) x-momentum: ρ∂u∂t+u∂u∂x+v∂u∂y=−∂p∂x+ρgx+μ∂2u∂x2+∂2u∂y2=0, {\text {x-momentum}}:\,\,\,\rho \left[ {{{\partial u} \over {\partial t}} + u{{\partial u} \over {\partial x}} + v{{\partial u} \over {\partial y}}} \right] = - {{\partial p} \over {\partial x}} + \rho {g_x} + \mu \left( {{{{\partial ^2}u} \over {\partial {x^2}}} + {{{\partial ^2}u} \over {\partial {y^2}}}} \right) = 0, (3) y-momentum: ρ∂v∂t+u∂v∂x+v∂v∂y=−∂p∂x+ρgy+μ∂2v∂x2+∂2v∂y2=0, {\text {y-momentum}}:\,\,\,\rho \left[ {{{\partial v} \over {\partial t}} + u{{\partial v} \over {\partial x}} + v{{\partial v} \over {\partial y}}} \right] = - {{\partial p} \over {\partial x}} + \rho {g_y} + \mu \left( {{{{\partial ^2}v} \over {\partial {x^2}}} + {{{\partial ^2}v} \over {\partial {y^2}}}} \right) = 0, where t is time and g_x = g_y = g is gravity acceleration.

Fig. 1

Vorticity-stream function scheme or vorticity transport equation [9] is broadly used rather than velocity-pressure scheme (primitive variables) as it has a clear advantage such that the pressure variable is unnecessary to be calculated in detailed manner [12]. In two-dimensional flow, every fluid particle will rotate because of its viscosity [13] as shown in Figure 2. Viscous fluid will support shear stress and shear stress will create a torque, τ. The coupling of torque will make the fluid particle to rotate in z-axis only. Vorticity transport equation can easily be obtained by calculating the curl [14] of Eq. (2) and Eq. (3). By incorporating stream function ψ, we get (4) u=∂ψ∂y, u = {{\partial \psi } \over {\partial y}}, (5) v=−∂ψ∂x, v = - {{\partial \psi } \over {\partial x}}, (6) ζ=∂ψ∂x−∂u∂y. \zeta = {{\partial \psi } \over {\partial x}} - {{\partial u} \over {\partial y}}. Then the vorticity-stream function equation can be written as, (7) ∂ζ∂t−∂ψ∂x∂ζ∂y+∂ψ∂y∂ζ∂x=μρ∂2ζ∂x2+∂2ζ∂y2. {{\partial \zeta } \over {\partial t}} - {{\partial \psi } \over {\partial x}}{{\partial \zeta } \over {\partial y}} + {{\partial \psi } \over {\partial y}}{{\partial \zeta } \over {\partial x}} = {\mu \over \rho }\left( {{{{\partial ^2}\zeta } \over {\partial {x^2}}} + {{{\partial ^2}\zeta } \over {\partial {y^2}}}} \right).

The vorticity, ζ at each of the walls for no-slip BC as shown in Figure 1 can be obtained via Taylor series expansion [15].

Fig. 2

The vorticity, ζ at left, right and bottom walls for free-slip BC as shown in Figure 2 can be obtained via Taylor series expansion as well. For top wall, its vorticity equation is the same for both no-slip and free-slip BCs. Vorticity, ζ at each wall is given by: (8) bottom BC: ζbottom=−(ψi,j=ny−1−ψi,j=ny)lsΔy+Δy22, {\rm{bottom}}\;{\rm{BC}}:\,\,\,{\zeta _{bottom}} = {{ - ({\psi _{i,j = ny - 1}} - {\psi _{i,j = ny}})} \over {{l_s}\Delta y + {{\Delta {y^2}} \over 2}}}, (9) left BC: ζleft=−(ψi=2,j−ψi=1,j)lsΔx+Δx22, {\rm{left}}\;{\rm{BC}}:\,\,\,{\zeta _{left}} = {{ - ({\psi _{i = 2,j}} - {\psi _{i = 1,j}})} \over {{l_s}\Delta x + {{\Delta {x^2}} \over 2}}}, (10) right BC: ζrigh=−(ψi=nx−1,j−ψi=nx,j)lsΔx+Δx22, {\rm{right}}\;{\rm{BC}}:\,\,\,{\zeta _{righ}} = {{ - ({\psi _{i = nx - 1,j}} - {\psi _{i = nx,j}})} \over {{l_s}\Delta x + {{\Delta {x^2}} \over 2}}}, (11) top BC: ζtop=−2(ψi,j=2−ψi,j=1)Δy2+2UwallΔy, {\rm{top}}\;{\rm{BC}}:\,\,\,{\zeta _{top}} = {{ - 2({\psi _{i,j = 2}} - {\psi _{i,j = 1}})} \over {\Delta {y^2}}} + {{2{U_{wall}}} \over {\Delta y}}, where Δx = Δy is step size, j is row, i is column, n is number of nodes and U_wall is an applied velocity. If slip length, l_s is equal to zero, Eq.(8) till Eq.(10) can be applied for no-slip BC.

In this section, few issues related to numerical modeling for solution of the two-dimensional vorticity transport equation will be analyzed. These issues are the application of the suitable numerical method, stability issue and iteration algorithm. Due to the nonlinear nature of second-order partial differential equation of the vorticity transport equation, there are few techniques accessible for the numerical solution such as finite element method, finite volume method, boundary element method and FDM [13]. The FDM is chosen among these available numerical methods. Despite being the simplest computational fluid dynamics method, it is conditionally stable [16]. In fact, this explicit FDM is sufficient enough to differentiate the shape of streamline distribution in no-slip and free-slip BCs of the unsteady square and rectangular cavity models with enhanced graphic user interface (GUI) to capture the rare formation of eddies. Figure 3 shows the explicit method with nodes that constitute the spatial and temporal approximations.

Fig. 3

Eq.(7) can be transformed into a finite difference equation (FDE) using forward and central finite-divided-difference formulas: (12) ζi,jt+1−ζi,jtΔt=−ψi,j+1t−ψi,j−1t2hζi+1,jt−ζi−1,jt2h+ψi+1,jt−ψi−1,jt2hζi,j+1t−ζi,j−1t2h+ μρζi+1,jt+ζi−1,jt+ζi,j−1t−4ζi,jth2. \matrix{ {{{\zeta _{i,j}^{t + 1} - \zeta _{i,j}^t} \over {\Delta t}}} \hfill & { = - \left( {{{\psi _{i,j + 1}^t - \psi _{i,j - 1}^t} \over {2h}}} \right)\left( {{{\zeta _{i + 1,j}^t - \zeta _{i - 1,j}^t} \over {2h}}} \right) + \left( {{{\psi _{i + 1,j}^t - \psi _{i - 1,j}^t} \over {2h}}} \right)\left( {{{\zeta _{i,j + 1}^t - \zeta _{i,j - 1}^t} \over {2h}}} \right)} \hfill \cr {} \hfill & { + \;{\mu \over \rho }\left( {{{\zeta _{i + 1,j}^t + \zeta _{i - 1,j}^t + \zeta _{i,j - 1}^t - 4\zeta _{i,j}^t} \over {{h^2}}}} \right).} \hfill \cr } Rearranging Eq. (12) yields: (13) ζi,jt+1=−Δtψi,j+1t−ψi,j−1t2hζi+1,jt−ζi−1,jt2h+Δtψi+1,jt−ψi−1,jt2hζi,j+1t−ζi,j−1t2h+ Δtμρζi+1,jt+ζi−1,jt+ζi,j+1t+ζi,j−1th2+1−4μΔtρh2ζi,jt , \matrix{ {\zeta _{i,j}^{t + 1} = - \Delta t\left( {{{\psi _{i,j + 1}^t - \psi _{i,j - 1}^t} \over {2h}}} \right)\left( {{{\zeta _{i + 1,j}^t - \zeta _{i - 1,j}^t} \over {2h}}} \right) + \Delta t\left( {{{\psi _{i + 1,j}^t - \psi _{i - 1,j}^t} \over {2h}}} \right)\left( {{{\zeta _{i,j + 1}^t - \zeta _{i,j - 1}^t} \over {2h}}} \right)} \hfill \cr { + \;\Delta t{\mu \over \rho }\left( {{{\zeta _{i + 1,j}^t + \zeta _{i - 1,j}^t + \zeta _{i,j + 1}^t + \zeta _{i,j - 1}^t} \over {{h^2}}}} \right) + \left( {1 - {{4\mu \Delta t} \over {\rho {h^2}}}} \right)\zeta _{i,j}^t ,} \hfill \cr } where Δt is step time, t is time, Δx = Δy = h is step size. 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 [17], such that (14) 1−4μΔtρh2≥0 or Δt≤ρh24μ. \left( {1 - {{4\mu \Delta t} \over {\rho {h^2}}}} \right) \ge 0\,\,\, {\rm or} \,\,\, \Delta t \le {{\rho {h^2}} \over {4\mu }}.

Eq. (14) is the conditional stability criterion. Substituting Eq. (4) and Eq. (5) into Eq. (6) yields: (15) −ζ=∂2ψ∂x2+∂2ψ∂y2. - \zeta = {{{\partial ^2}\psi } \over {\partial {x^2}}} + {{{\partial ^2}\psi } \over {\partial {y^2}}}.

Eq. (15) is called Poisson equation [14]. Similarly, Eq. (15) can be discretized using FDM: (16) ψi,jn+1=14(ψi+1,jn+ψi−1,jn+ψi,j+1n+ψi,j−1n+h2ζi,jn), \psi _{i,j}^{n + 1} = {1 \over 4}(\psi _{i + 1,j}^n + \psi _{i - 1,j}^n + \psi _{i,j + 1}^n + \psi _{i,j - 1}^n + {h^2}\zeta _{i,j}^n), where n is the iteration step. Eq.(16) can later be solved by using Gauss-Seidel method [14]. The time-marching method to solve the forward finite difference equation of the vorticity transport equation is carried out as below (in Matlab_®):

1. Initial value of stream function, ψ and vorticity, ζ are set to zero at the 1^st step time, t = 0. Values of stream function, ψ at top, bottom, left and right boundaries are set to zero along the iterative process from the 1^st step time, t = 0 till the last step time, t = t_end.

The present finite difference Matlab_® code used in this research is self-developed and it is embedded into a GUI solver that is user-friendly and user-interactive. The GUI was built with multiple of buttons and columns to let the user to key-in the data. It facilitates the user to study the fluid flow distribution throughout the lid-driven cavity and to obtain the readings from the two-dimensional fluid flow profile. This GUI is consisted of three panels, four axes, seven buttons, ten columns and one radio button. The ten columns are the user input panel i.e. dynamic viscosity (μ), density (ρ), velocity (u), length of cavity (l), height of cavity (H), slip length (l_s), number of nodes (n), time increment (Δt), total time (t) and steady-state criterion (e). Once the end user has assigned a value to each of these inputs, the simulation can be executed by clicking the Plot button. Then results such as stream function distribution, u-velocity distribution, v-velocity distribution and velocity vector distribution can be well displayed in animation style.

Their associated finite difference equations are: (17) ui,j=ψi,j+1−ψi,j−12h, {u_{i,j}} = {{{\psi _{i,j + 1}} - {\psi _{i,j - 1}}} \over {2h}}, (18) vi,j=−ψi+1,j−ψi−1,j2h. {v_{i,j}} = - {{{\psi _{i + 1,j}} - {\psi _{i - 1,j}}} \over {2h}}. Velocity vector is given by: (19) V=ui,j2+vi,j2. V = \sqrt {u_{i,j}^2 + v_{i,j}^2} . Direction of velocity vector, V is given by: (20) θ=tan−1vi,jui,j. \theta = {\tan ^{ - 1}}\left( {{{{v_{i,j}}} \over {{u_{i,j}}}}} \right).

In this work, the model used in Ref. [10] will be studied first and served as a benchmark in this paper to verify the self-developed Matlab code and to investigate the developing numerical error if any. Two dimensional cavity flow (with no-slip BC) were carried out based on the numerical method as discussed in Section 3 in such a way that the flow is laminar and unsteady. The iterative procedure of the numerical process is performed until the flow is in a steady-state condition when criterion of e ≤ 10⁻⁹ is satisfied.

The steady state criterion is given by: (21) e=|max(max(sfj,i,n+1))−max(max(sfj,i,n))||max(max(sfj,i,n+1))|, e = {{|max (max (s{f_{j,i,n + 1}})) - max (max (s{f_{j,i,n}}))|} \over {|max (max (s{f_{j,i,n + 1}}))|}}, where n is iterative time, s f_j,i,n+1 is the stream function at current time and s f_j,i,n is the stream function at previous time.

The test case problem of the two-dimensional lid-driven laminar cavity flow is as shown in Figure 4. The fluid density, ρ is set to 1 kg/m³, dynamic viscosity, μ = 0.001Pa.s, applied velocity, U_wall = 0.2m/s. Step time, Δt is set to 0.1 s and the computation is performed on 40 × 40 grids (uniform grid and Δx = Δy). The Reynolds number, Re can be reckoned as [19]: (22) Re=UwallρHμ, {\mathop{ Re}\nolimits} = {{{U_{wall}}\rho H} \over \mu }, where H is the height of the cavity.

Fig. 4

5.1

First model [20] and validation

Computation is performed on side length of the square cavity with 1 m in no-slip BC. The numerical solution of u− and v− velocity distributions were obtained for cavity flow with Re = 200. A time series of contour plots that shows the growing streamlines under the top moving lid in the u− and v−velocities in Figure 5. An area of swirling motion was created under the top lid until it disappears as it reaches the bottom surface of the cavity. Then computational job is continued until the steady-state criterion of e ≤ 10⁻⁹ is satisfied where the results are shown in Figure 6(a) and Figure 6(c). The steady state results are then compared and verified by using the commercial simulation software Ansys Fluent_® utilizing FEM as shown in Figure 6(b) and Figure 6(d). Additionally, the contour plots of the streamlines of u− and v−velocities obtained from Zhu [10] are compared with those results as shown in Figure 5 and Figure 6. Besides, the contour plots of the stream function distribution along the iteration in transient analysis are shown in Figure 7. Two small eddies (which rotate in opposite direction with respect to primary main cell) are formed at both of the corners at the bottom, whereas the right bottom eddy seems bigger as shown in Figure 7(e). The growing/formation of eddies has become more vivid through Figure 5–7.

Fig. 5

Fig. 6

Fig. 7

5.2

Second model

An interesting phenomenon is observed when the first model is subject to a free-slip BC at l_s = 1 m where the shape of the primary cell is being ‘stretched’ to the corners at the bottom. It means the flow dynamics has reached the entire cavity sufficiently that there is no room for corner eddy to form as observed in Figure 8 while contour plots of the streamlines of u− and v− velocities are as shown in Figure 9. This free-slip BC model however takes a longer period in order to achieve a better steady-state condition at e ≤ 10⁻⁹. The results in Figure 8–9 shows the significant role of free-slip to eliminate eddies formation as compared to no-slip case in the first model.

Fig. 8

Fig. 9

5.3

Third model (with Moffat’s aspect ratio)

The AR of the first model is modified from AR=1 to AR=1.629 as suggested by Moffatt [6] in order to investigate how the two corner eddies at the bottom merge to form a secondary main cell at this critical AR (in no-slip BC). A time series of contour evolution of the streamline, u− and v−velocity distributions at different steady-state criteria (from e ≤ 10⁻² to e ≤ 10⁻¹⁰) and the growing formation of secondary eddies are shown in Figure 10 and Figure 11 respectively. It is observed that under no-slip condition, the growing of eddies becomes more significant as compared to observation under free-slip BC. This is because the acceleration of the fluid at the wall is assumed to be zero in no-slip BC while there is a relative velocity between the wall and the fluid in free-slip BC. The value of maximum stream function, ψ_max in each step time increases along the iteration, then it achieves plateau before approaching the highest point and drops again in a tiny quantity until the end of the iteration as shown in Figure 12. Based on Figure 12, the streamline distribution achieves the highest maximum value at t = 65.6 s (at about steady-state criterion of e ≤ 10⁻⁵). The location of maximum stream function, ψ_max is fixed in position and does not move anywhere after this steady-state criterion. Then the value of ψ_max begins to drop a little bit after t = 65.6 s. Secondary eddies are formed at t = 51.6 s and they grown in size further. Then these two secondary eddies merged at t = 152.9 s (at about steady state criterion of e ≤ 10⁻⁸) and the shape and size of the whole streamline distribution do not change at any further iteration as shown in Figure 10(g).

Fig. 10

Fig. 11

Fig. 12

According to suggestion by Moshikin and Poochinapan [9], the steady flow is achieved when the criterion, e ≤ 10⁻⁸ is achieved. However, Abd-El-Hafiz, Ismail and Matit [18] proposed another criterion in which the iteration is repeated until the velocity profile (velocity and stream function are interrelated through Eq. (4) and Eq. (5)) converges when the criterion of e ≤ 10⁻⁵ is reached. A secondary main cell is not able to form if iteration stops as the steady-state criterion e ≤ 10⁻⁵ is achieved. In addition to this, by performing an iteration beyond the steady-state criterion e e ≤ 10⁻⁸, it does not contribute to any useful result due to iteration and round off errors [19]. It is observed that the value of streamline distribution drops after achieving the steady-state criterion at e ≤ 10⁻¹⁰. If more iterations are needed in order to achieve the smaller steady-state criterion, then the resultant steady-state period will be longer. Hence it will cause a longer execution time.

Based on Table 1, the highest value of maximum stream function, ψ_max for the steady state criterion e ≤ 10⁻⁵ can be found occurs at t(ψ_max) = 65.6 s by using the code as shown below: (23) t(ψmax)=fing(max(max(sf(:,:,:)))=max(max(max(sf(:,:,:))))), t({\psi _{max }}) = fing(max (max (sf(:,:,:))) = max (max (max (sf(:,:,:))))), where s f (:, :, :) is a three-dimensional matrix of stream function along the simulation period. The value and location of the highest value of maximum stream function, ψ_max as shown in Figure 10 (d) can be detected via the code shown as below: (24) ψmax=max(max(sf(:,:,t(ψmax))))), {\psi _{max }} = max (max (sf(:,:,t({\psi _{max }}))))), (25) [y,x]=fing(max(max(sf(:,:,t(ψmax))))=max(max(max(sf(:,:,t(ψmax)))))). [y,x] = fing(max (max (sf(:,:,t({\psi _{max }})))) = max (max (max (sf(:,:,t({\psi _{max }})))))). Similarly, the value and location of minimal stream function, ψ_min can be obtained via similar method as shown in Eq. (24) and Eq. (25) by changing the keyword from ‘max’ to ‘min’. Steady-state period can be obtained when an iterative process is aborted at certain time after the desired criterion, e is satisfied.

Table 1

Nomenclatures.

Symbol	Meanings	Symbol	Meanings
ρ	density [kg/m3]	ψ	stream function [m2/s]
U	applied velocity [m/s]	ψmax	maximum stream function [m2/s]
L	length of cavity [m]	Uwall	applied velocity at top wall [m/s]
μ	dynamic viscosity [Pa.s]	n	non-negative integer, number of nodes
Re	Reynolds number	i	column
A	aspect ratio	j	row
H	height of cavity [m]	Δx	step size [m]
u	velocity in x–direction [m/s]	Δy	step size [m]
v	velocity in y–direction [m/s]	us	slip velocity [m/s]
g	gravity acceleration [m/s2]	ls	slip length [m]
u	applied velocity [m/s]	gx	gravity acceleration in x–direction [m/s2]
θ	angle [deg]	gy	gravity acceleration in y–direction [m/s2]
t	time [s]	Δt	step time [s]
∇	Gradient operator	h	step size [m]
ζtop	vorticity at the top [m/s2]	s fj,i,n	stream function at previous time [m2/s]
ζleft	vorticity at the top [m/s2]	s fj,i,n+1	stream function at previous time [m2/s]
V	velocity vector [m/s]	ζright	vorticity at the right [m/s2]
ζ	vorticity [m/s2]	ζbottom	vorticity at the bottom [m/s2]
e	steady-state criterion

5.4

Fourth model

The third model is analyzed again until the steady-state criterion, e ≤ 10⁻⁸ is achieved but in free-slip BC by gradually increasing the slip length, l_s value. It seems that when slip length, l_s passes through the threshold value of 0.07m, the secondary main cell of eddy is not able to be formed as shown in Figure 13(a). Threshold free-slip means the value of slip length that is just nice to eliminate the formation of secondary main cell of eddy while significant free-slip means an amplified slip length to thoroughly change the shape of streamline. The associated contour plots of the streamline, u− and v− velocity distributions are shown in Figure 14. By further increasing the slip length, l_s value to an extremely high free-slip BC at l_s = 1 m, the shape of the primary cell is being ‘stretched’ to the corners of the bottom cavity wall as shown in Figure 13(b) (similar with Figure 8). Thus, there is no room for eddy formation under this free-slip condition. Meanwhile, its associated contour plots of the streamline, u− and v− velocity distributions are shown in Figure 15 (similar with Figure 9). Further reference on slip length measurement is given by Maali and Bhushan [21].

Fig. 13

Fig. 14

Fig. 15

This study is a unique attempt to consider both no-slip and free-slip effects in both square and rectangular cavity models under different aspect ratios and steady-state criteria which were not considered before in the existing literature. All the FDM codes are self-developed in Matlab_® environment and they are embedded into a self-designed GUI that is user-interactive for beginners in this field by the authors. The first model is executed in order to verify the presently self-developed Matlab_® FDM code with Ansys Fluent_® which employs finite element method. Although this scheme seems simple, it is conditionally stable and sufficient to differentiate the flow pattern in both no-slip and free-slip BCs quantitatively by observing the streamline distribution. The changed shape of streamline distribution in free-slip effect is important for subsequent heat transfer analysis as it will amend temperature distribution of the lid-driven model as well. Thus, it further influences the convection heat transfer process of lid-driven model when subject to ambient air or radiation heat transfer. Also, we are able to show the formation of secondary main cell when lid-driven is subject to critical AR as mentioned by Moffatt [6] in the third model. It is noticed that no significant change to the results and contour plots of streamline and velocity distributions are found by executing the iteration smaller than the steady-state criterion of e ≤ 10⁻⁸. Hence, it would only bring about numerical error, longer computational effort and hence it results in futile. The merging of secondary main cell can be refrained when the slip length (free-slip BC) achieves the threshold value. Stretching effect with free-slip BC in the fourth model regulates the fluid dynamics to reach the entire cavity sufficiently with no room for eddy formation by increasing the slip length to an intense value. With help of the presently self-developed GUI, the rare observations on eddies formation in the square and rectangular lid-driven cavity models under influence of both no-slip and free-slip effects which are never considered in a single existing literature can be captured sufficiently and pointed out clearly to promote uniquely extended findings on the subject of interest.

The parameters given by Table 2 are used to explain κ₁ =Steady state criteria e, κ₂ =Period t to achieve designated e, κ₃ = ψ_max(m²/s) at designed e. κ₃ =Location of ψ_max(m) at designated e, respectively.