J. Nig. Soc. Phys. Sci. 4 (2022) 231–241 Journal of the Nigerian Society of Physical Sciences Computational investigation of free convection flow inside inclined square cavities Mariya Helen Mercy JK, Prabhakar.V∗ Division of Mathematics, School of Advanced Sciences, Vellore Institute of Technology, Chennai, Tamil Nadu, India. Abstract The temperature and fluid profiles of flow inside tilted square cavities are analysed with two different cases of thermal boundary conditions, (1) Isothermally cold sidewalls of the cavity and the hot bottom wall kept parallel to the insulated top wall, (2) Hot left wall, cold right wall, insulated top and bottom walls. The Galerkin finite element method with penalty parameter is used to solve the nonlinear coupled system of partial differential equations governing the flow and thermal fields. The method is further used to solve the Poisson equation for stream function. The results are presented in terms of isotherms and streamlines. The Gaussian rule with the hybrid function formed from the block-pulse function and Lagrange polynomial is implemented for the evaluation of the definite integrals present in the residual equations. Attempting to affix the hybrid methods in the integration part for solving Finite Element Method (FEM) turned efficacious as the convergence is achieved for streamlines and isotherms with the existing results. The tilted square cavities with inclination angles 0 ◦ ≤ φ ≤ 90 ◦ and Rayleigh number ranging between 103 ≤ Ra ≤ 105 for Pr = 0.71 (air) are investigated. The source code for the finite element analysis is written in Mathematica. The results thus obtained are found to be competent with those of COMSOL, the Software for Multiphysics Simulation. DOI:10.46481/jnsps.2022.637 Keywords: Penalty finite element analysis, inclined square cavities, thermal distribution, fluid distortion, Mathematica, hybrid function, numerical integration, COMSOL Article History : Received: 04 February 2022 Received in revised form: 07 April 2022 Accepted for publication: 07 April 2022 Published: 29 May 2022 c©2022 Journal of the Nigerian Society of Physical Sciences. All rights reserved. Communicated by: T. Latunde 1. Introduction The term ‘convection’ is derived from the Latin convectare, which means to take to a place. During the transmission of heat, convection can be defined as the drift of heat energy from or to a solid because of nearby fluid in motion, in the existence of a temperature gradient. Further classified into two major types depending upon the source of the fluid motion: Forced con- vection, where the fluid movement happens due to an external ∗Corresponding author tel. no: +914439931239 Email address: prabhakar.v@vit.ac.in (Prabhakar.V) influence, such as a fan, extractor, pump etc. It is called as a natural convection, when the fluid movement is solely due to the characteristics of the fluid between two determined points of the process. The fact that the velocity and the temperature equations are solved all together makes the study of natural convection seemingly complex. The natural convection within a cavity plays important roles in many science and engineering applications namely heat exchangers [1], solar stills [2], cooling of electric and electronic components [3], nanoscience [4]-[6], etc. Square cavity with obstacles, lid-driven, heated walls are also in the scope of the researchers [7]-[16]. Navier stokes equations governs the fluid flow, whereas the en- 231 Mariya & Prabhakar / J. Nig. Soc. Phys. Sci. 4 (2022) 231–241 232 ergy balance equations regulate the nature of the temperature. Among the various methods utilized to handle square cavities, FEM, a very competent method in dealing complex geome- tries and unstructured domains, is employed in this work. As the complexity increases, the integration process to calculate the stiffness and mass matrix becomes more complicated. The penalty finite element method [17] is involved in the nonlinear system of dimensionless governing equations of the 2D flow problem to waive the pressure constant. Two non-dimensional numbers, Prandtl number (Pr), signifying how much heat is carried away by how much fluid, transferring from one point to another and Rayleigh Number (Ra), describing the relation- ship between buoyancy and viscosity within the fluid, are ex- plored for the convectional fluid flow. The tilted square cavi- ties with inclination angles 0 ◦ ≤ φ ≤ 90 ◦ and Rayleigh number ranging between 103 ≤ Ra ≤ 105 for Pr = 0.71 (air) are inves- tigated. The Rayleigh number above 103 will display a distinc- tive change in the convective heat transfer rate as the inclination angle is increased. At Ra = 103 and Ra = 105 gives the mini- mum and maximum heat transfer rate as the angle increases [8]. Literature provides with various methods for solving ODEs and PDEs and differential equations are solved using block meth- ods, step points, numerical methods, etc., [18]-[22].In common practice, Gaussian quadrature method is used for evaluating the integrals present in the finite element equations. Among the numerical methods available in the literature, a technique based on the Hybrid function approximation for solving non- linear initial-value problems with applications to Lane-Emden type equations was suggested [23]. They made use of the prop- erties of block-pulse functions and Lagrange interpolating poly- nomials to reduce the nonlinear initial-value problem to a sys- tem of non-algebraic equations to arrive at the accurate results. The same method was also used to provide accurate results for solving nonlinear integro-differential equations such as in Volterrra’s population model [24] and Volterra-Fredholm equa- tion, and also for variational problems. A comparative study of quadrature rules based on Haar wavelet and Hybrid function of block-pulse and Legendre polynomial, for finding the approx- imate value of the definite integrals, was made and the proce- dure was extended to the numerical solution of double and triple integrals with variable limits [25, 26]. It was observed that the hybrid method provides faster convergence when related to Haar Wavelet and that the orders of the block-pulse function and Legendre polynomial can be attuned to attain very precise solution. Hybrid of Block-Pulse function using Lagrange poly- nomial was considered for evaluation by [27] for the evalua- tion of general double and triple integrals with variable limits that shows better accuracy over Haar Wavelet. Solution of vari- ational problems are also derived using the same [28]. Vari- able weights are used for the study in contrast with the constant weights considered in [25]. The square cavities inclined at angles 0 ◦ ≤ φ ≤ 90 ◦ are studied with the fluid (Pr = 0.71) flowing within a two dimensional, steady and laminar flow across the calculation domain with two different cases of thermal boundary conditions (BC), (Case 1) Isothermally cold sidewalls of the cavity and the hot bottom wall kept parallel to the insulated top wall, Figure 1. Square tilted at angle (Case 1) φ (Case 2) Hot left wall, cold right wall, insulated top and bottom walls. are featured in the present work. An endeavour is made to in- corporate this hybrid method in the integration part for solv- ing FEM turned efficacious as the convergence is achieved for streamlines and isotherms with the existing results in the litera- ture [10, 12]. The source code for the finite element analysis is written in Mathematica and the results obtained are compared with COMSOL. Mathematical Procedure involving the Govern- ing Differential Equations (GDEs) and finite element formu- lation is discussed in the first section, Hybrid method for the numerical integration is explained in the next section, and the results are detailed in the last section. 2. Methodology 2.1. Computational technique The physical domain (D) of square cavity (PQRS) with angles of inclination 0 ◦ ≤ φ ≤ 90 ◦ for both thermal bound- ary conditions is considered in the present study. The range of Rayleigh’s Numbers falls between 103 ≤ Ra ≤ 105 for the study. Figure 1, shows a clear picture of the geometry consid- ered, incorporating case 1 boundary conditions. Physical properties are kept constant during the calculation except the density in buoyancy term, where change in den- sity due to temperature variation is estimated using Boussi- nesq approximation. The phenomena of thermal and fluid flow inside the domain are governed by the energy balance and Navier–Stokes equations, respectively. The governing equa- tions for steady natural convection flow using conservation of mass, momentum and energy in dimensionless form are given below: ∂U ∂X + ∂V ∂Y = 0 (1) U ∂U ∂X + V ∂U ∂Y = − ∂P ∂X + Pr ( ∂2U ∂X2 + ∂2U ∂Y 2 ) (2) 232 Mariya & Prabhakar / J. Nig. Soc. Phys. Sci. 4 (2022) 231–241 233 U ∂V ∂X + V ∂V ∂Y = − ∂P ∂Y + Pr ( ∂2V ∂X2 + ∂2V ∂Y 2 ) + Ra Pr θ (3) U ∂θ ∂X + V ∂θ ∂Y = ∂2θ ∂X2 + ∂2θ ∂Y 2 (4) The penalty finite element method [17] is used to eliminate the pressure term by the inclusion of the penalty parameter (γ) in the equations (2) and (3) with the following relationship in- volving the incompressibility. P = −γ ( ∂U ∂X + ∂V ∂Y ) (5) Generally, (γ = 107) is considered for reliable solutions [10]. Applying (5) to (2) and (3) yields U ∂U ∂X + V ∂U ∂Y = −γ ∂ ∂X ( ∂U ∂X + ∂V ∂Y ) + Pr ( ∂2U ∂X2 + ∂2U ∂Y 2 ) (6) U ∂V ∂X + V ∂V ∂Y = −γ ∂ ∂Y ( ∂U ∂X + ∂V ∂Y ) + Pr ( ∂2V ∂X2 + ∂2V ∂Y 2 ) + Ra Pr θ (7) The whole domain is discretized into bi-quadratic elements. Galerkin finite element method is applied to solve the system of governing differential equations (6), (7) and (4). A numerical integration method involving hybrid functions explained in the section 2.3 is used for obtaining the finite el- ement solutions from the residual equations. The thermal (θ) and the velocity components (U&V ) are expanded through ba- sis set {∑9 j=1 N j } given in equation (8). N j s are the shape func- tions. The procedure to obtain N j s are detailed in [29] and those functions considered for this study is provided in the Appendix. U ≈ 9∑ j=1 U j N j(X, Y ), V ≈ 9∑ j=1 V j N j(X, Y ), θ ≈ 9∑ j=1 θ j N j(X, Y ) (8) The residual equations for an element are obtained by using (8) in the governing equations, resulting in the following system of nonlinear residual partial differential equations. R(1)i = 9∑ j=1 U j × ∫ D   9∑ j=1 U j N j  ∂N j∂X +  9∑ j=1 V j N j  ∂N j∂Y  NidXdY + γ  9∑ j=1 U j ∫ D ∂Ni ∂X ∂N j ∂X dXdY + 9∑ j=1 V j ∫ D ∂Ni ∂X ∂N j ∂Y dXdY  + Pr 9∑ j=1 U j ∫ D ( ∂Ni ∂X ∂N j ∂X + ∂Ni ∂Y ∂N j ∂Y ) dXdY (9) R(2)i = 9∑ j=1 V j × ∫ D   9∑ j=1 U j N j  ∂N j∂X +  9∑ j=1 V j N j  ∂N j∂Y  NidXdY + γ  9∑ j=1 U j ∫ D ∂Ni ∂Y ∂N j ∂X dXdY + 9∑ j=1 V j ∫ D ∂Ni ∂Y ∂N j ∂Y dXdY  + Pr 9∑ j=1 V j ∫ D ( ∂Ni ∂X ∂N j ∂X + ∂Ni ∂Y ∂N j ∂Y ) dXdY − Ra Pr 9∑ j=1 θ j ∫ D N jdXdY (10) R(3)i = 9∑ j=1 θ j × ∫ D   9∑ j=1 U j N j  ∂N j∂X +  9∑ j=1 V j N j  ∂N j∂Y  NidXdY + 9∑ j=1 θ j ∫ D ( ∂Ni ∂X ∂N j ∂X + ∂Ni ∂Y ∂N j ∂Y ) dXdY (11) For nine-noded biquadratic element (Figure 2) with three degrees of freedom, residual equations (9)− (11) consists of 27 unknowns in 27 equations. 2.2. Stream function The stream function is used to display the fluid flow and is acquired from velocity components U and V . The relationships between stream function, ψ and velocity components for 2D flows are U = ∂ψ ∂Y and V = − ∂ψ ∂X (12) On differentiating (12), the governing equation for stream function is attained. ∂2ψ ∂X2 + ∂2ψ ∂Y 2 = ∂U ∂Y − ∂V ∂X (13) Expanding the stream function ψ using the basis set ψ ≈ 9∑ j=1 ψ j N j(X, Y ) (14) and the relation for U, V from (8), the residual equations for (13) is, R(ψ)i = 9∑ j=1 ψ j ∫ D ( ∂Ni ∂X ∂N j ∂X + ∂Ni ∂Y ∂N j ∂Y ) dXdY − ∫ Γ Nin.∇ψdΓ + 9∑ j=1 U j ∫ D ∂Ni ∂N j ∂Y dXdY − 9∑ j=1 V j ∫ D ∂Ni ∂N j ∂X dXdY (15) Applying no-slip boundary condition (ψ = 0), ψ’s is ob- tained from the linear system (15). Stream functions (ψ′s) thus obtained might be positive or negative denoting the anti- clockwise and clockwise circulation respectively. The integrands of the definite integrals appearing in (9)-(11) and (15) are functions of the global coordinates X and Y . Fig- ure 2 shows the co-ordinate transformation for the discretized elements from the X − Y plane to s − t the plane [11]. The 233 Mariya & Prabhakar / J. Nig. Soc. Phys. Sci. 4 (2022) 231–241 234 Figure 2. Coordinate transformation of the discretized element form X − Y to s − t plane Table 1. Isotherms of Ra = 103 ; Pr = 0.71; PS & QR-cold; SR-adiabatic; PQ-hot φ Mathematica results COMSOL results 15 ◦ 45 ◦ 75 ◦ integrals appearing in the transformed residual equations are evaluated using hybrid function of block pulse function and La- grange polynomial. The transformed residual equations solved for every node in the domain provides the thermal and velocity components. Finite element procedure to solve this is briefed in the Appendix. It was shown that hybrid function based on block-pulse and Lagrange polynomial gives better accuracy than Haar wavelets and other hybrid functions [26], with comparatively lesser nodal points. Motivated by the accuracy of this method in solv- ing the differential equations (both linear and nonlinear), an in- tegration scheme based on hybrid functions for definite single, double and triple integrals was presented by [27]. An attempt has been made in applying hybrid functions to obtain the finite element solutions. Details regarding the hybrid function and the integration scheme are explained in the following section. Table 2. Isotherms of Ra = 104 ; Pr = 0.71; PS & QR-cold; SR-adiabatic; PQ-hot φ Mathematica results COMSOL results 15 ◦ 45 ◦ 75 ◦ Table 3. Isotherms of Ra = 105 ; Pr = 0.71; PS & QR-cold; SR-adiabatic; PQ-hot φ Mathematica results COMSOL results 15 ◦ 45 ◦ 75 ◦ 2.3. Hybrid Function Definition 1: Block pulse functions: A set of block-pulse function φ j(t), j = 1, 2, ...J defined on the interval [0, 1) are denoted as φ j(t) = 1, t j−1 ≤ t ≤ t j0, otherwise (16) where J represent the number of partitions of [0, 1) or the order 234 Mariya & Prabhakar / J. Nig. Soc. Phys. Sci. 4 (2022) 231–241 235 Table 4. Streamlines of Ra = 103 ; Pr = 0.71; PS & QR-cold; SR-adiabatic; PQ-hot φ Mathematica results COMSOL results 15 ◦ 45 ◦ 75 ◦ Table 5. Streamlines of Ra = 104 ; Pr = 0.71; PS & QR-cold; SR-adiabatic; PQ-hot φ Mathematica results COMSOL results 15 ◦ 45 ◦ 75 ◦ of the block pulse function in [0, 1). Definition 2: Hybrid functions of block pulse and Lagrange polynomial functions: A set of hybrid functions h jk(t), j = 1, 2, ...J, Table 6. Streamlines of Ra = 105 ; Pr = 0.71; PS & QR-cold; SR-adiabatic; PQ-hot φ Mathematica results COMSOL results 15 ◦ 45 ◦ 75 ◦ Table 7. Isotherms of Ra = 103 ; Pr = 0.71; PS-hot; QR - Cold; PQ and SR -adiabatic; φ Mathematica results COMSOL results 0 ◦ 30 ◦ 50 ◦ k = 0, 1, ...K − 1 on the interval [0, 1) are denoted as h jk(t) = Lk(2Jt − 2 j + 1), t ∈ [ j−1 J , j J ] 0, otherwise (17) where Lk(t) denoting the Lagrange polynomial of order K. As the block pulse functions and Lagrange interpolating polynomi- als are complete and orthogonal, { h jk } is a complete orthogonal set in the Hilbert space L2 [23]. 235 Mariya & Prabhakar / J. Nig. Soc. Phys. Sci. 4 (2022) 231–241 236 Table 8. Isotherms of Ra = 104 ; Pr = 0.71; PS-hot; QR - Cold; PQ and SR -adiabatic; φ Mathematica results COMSOL results 0 ◦ 30 ◦ 50 ◦ Table 9. Isotherms of Ra = 105 ; Pr = 0.71; PS-hot; QR - Cold; PQ and SR -adiabatic; φ Mathematica results COMSOL results 0 ◦ 30 ◦ 50 ◦ 2.3.1. Approximation of function of a single variable: A function f (t) ∈ L2 [ 0, 1 ) can be approximated, in terms of the hybrid functions as f (t) = ∞∑ j=1 ∞∑ k=0 c jkh jk(t) (18) For an interval [0, 1) divided into J partitions and K internal Table 10. Streamlines of Ra = 103 ; Pr = 0.71; PS-hot; QR - Cold; PQ and SR -adiabatic; φ Mathematica results COMSOL results 0 ◦ 30 ◦ 50 ◦ Table 11. Streamlines of Ra = 104 ; Pr = 0.71; PS-hot; QR - Cold; PQ and SR -adiabatic; φ Mathematica results COMSOL results 0 ◦ 30 ◦ 50 ◦ nodes in each partition, f (t) is taken in the form f (t) ≈ J∑ j=1 K−1∑ k=0 c jkh jk(t) (19) Considering K nodes in the jth partition [ j−1 J , j J ] , denoted by r jk, it can be shown that c jk = f (r jk) [26]. 236 Mariya & Prabhakar / J. Nig. Soc. Phys. Sci. 4 (2022) 231–241 237 Table 12. Streamlines of Ra = 105 ; Pr = 0.71; PS-hot; QR - Cold; PQ and SR -adiabatic; φ Mathematica results COMSOL results 0 ◦ 30 ◦ 50 ◦ 2.3.2. Approximation of function of two variables: Extending the function approximation of a single variable in terms of hybrid functions, the function f (t, u) defined in [0, 1) × [0, 1) can be approximated as f (t, u) ≈ J′∑ j=1 K′−1∑ k=0 J∑ j=1 K−1∑ k=0 c jklmh jk(t)hlm(u) (20) by considering J partitions and K internal nodes in each partition for [0, 1) in the t-direction, and J′ partitions and K′ internal nodes in each partition for [0, 1) in the u-direction. It can be shown that c jklm = f ( r jk, rlm ) (21) where Gaussian nodes along the t and u directions are and r jk and rlm respectively [25]. 2.3.3. Integration scheme for the integrals in residual equa- tions: Using the function approximation from equations (19) and (20), the integrals (both single and double) of the residual equations can be easily evaluated as follows. ∫ 1 0 f (t)dt = J∑ j=1 K−1∑ k=0 c jkw jk (22) where w jk = ∫ 1 0 h jk(t)dt are the weights. As the hybrid function h jk(t) is a polynomial, its weight can be easily calculated. Similarly, ∫ 1 0 ∫ 1 0 f (t, u)dtdu = J′∑ j=1 K′−1∑ k=0 J∑ j=1 K−1∑ k=0 c jklmw jkwlm (23) where w jk = ∫ 1 0 h jk(t)dt and wlm = ∫ 1 0 hlm(u)du are the weights. For single integral J = 1 and K = 2; for double integrals J = 2, K = 4, J′ = 3 and K′ = 5 are well accommodated for the integrals present in the residual equations. 3. Results and discussion We consider a two dimensional, steady and laminar flow across the inclined square cavity (PQRS) under two different thermal boundary conditions as mentioned in section 1. The physical domain (D) of square cavity with angle of inclination (φ), inclined with the X-axis within the acute angles of 15 ◦ , 45 ◦ and 75 ◦ for case 1 and 0 ◦ , 30 ◦ and 50 ◦ for case 2 thermal boundary conditions is considered for the present study. Pr = 0.7, and Ra varied from 103 to 105 for both the boundary conditions are investigated. In order to assess the accuracy of our numerical procedure, we have tested our algorithm for the results presented in [10] for 28 × 28 elements and found to be exactly matching. In this paper, however, we restrict our study considering 20 × 20 elements (Figure 3) with the grid size of 51 × 51, which were also found to be in good agreement with the isotherms and streamlines reported in [10]. Distributions of isotherms and streamlines of case (1) are portrayed in Tables 1-3 and 4-6 respectively; Similarly Tables 7-9 and 10-12 illustrates the heat and fluid contours of case (2). Observing the isotherms for case (1), when Ra is 103, however inclined the square cavity is, there is not much disruption, indicating heat transfer through conduction as depicted in Table 1, whereas, when Ra is increased, isotherms are concentrated at the hot and cold walls. Both clockwise and anticlockwise circulations seen in all the Rayleigh numbers are pondered at φ = 15 ◦ . When Ra is 103, the secondary cells are comparatively lesser in size but when it is 104 the cells are slightly bigger and when the value of Ra is 105 the two axisymmetric flow exquisitely occupies the entire cavity. As φ increases to 45 ◦ , the strength of anticlockwise circulation cells increases (Table 5). As the angle of inclination increases, there is a significant push of isotherms towards the right wall. At φ = 75 ◦ , the isotherms are found to be qualitatively similar to those of φ = 45 ◦ as the temperature contours are piling towards the right wall (Tables 2 and 3). The overall amount of heat transfer along the right wall increases with inclination angle and decreases along the left wall. In Table 6 at 45 ◦ , strong primary anticlockwise circula- tion cells occupy almost the entire part of cavity except top corner of the wall QR. As inclination angle φ further increases to 75 ◦ , the primary circulation cells span the entire cavity whereas secondary circulation cells completely disappear. It is interesting to observe that the centre of left vortex for the fluid 237 Mariya & Prabhakar / J. Nig. Soc. Phys. Sci. 4 (2022) 231–241 238 circulations is pushed towards the bottom wall and centre of right vortex is pushed towards the top wall due to enhanced tangential components of buoyancy force along the bottom wall which tends to impose an anticlockwise fluid circulation. Exploring case (2), the temperature contours exhibit a gentle migration from left to right side of the calculation domain for all the investigated. In Table 7, dealing with the Ra = 103 , the contours are observed to exhibit an elegant protuberance toward the right wall as there is a surge in the inclination angles. Streamlines in these conditions are not much varied but are laminar in nature as in Tab.12. With increase in Ra, heat transfer due to convection is clearly visible (Table 9). The difference is observed as φ is gradually amplified. The fluid flow is more unruffled than in the case 1. Only when Ra = 105 and φ = 0 ◦ , the centre of the vortex is not circular in shape as in Table 12. The flow in the circulation cells are almost circular in shape near the core and further, they swell and take the shape of the cavity near the walls due to strong convection only increasing the angle of inclination. From what is seen in the isotherms the convection mode of heat transfer dominates over the conduction, as the inclination angle increases. Evident picture of the conduction effect can be seen in the contour for temperature. The movement of fluid backs up the convection, due to the bulk movement of air triggers the convective heat transfer moves the constant temperature line in the direction of its motion. Hence, the circulation cells are strained along these walls, prevailing in the sizeable portion of the central region. The temperature layer in the core breaks down because of the good mixing occurring in the core. The inclined square is the calculation domain, hence the effect of tangential and normal components of buoyancy force relative to hot wall gains importance for flow and thermal phenomena. Stronger anticlockwise fluid circulations are witnessed when increasing the angle of inclination as the buoyancy force along the hot wall also grows. Results indicate that the fluid circulations, isotherms are strongly dependent on the inclination angle of the square cavity. The heat and the flow contours obtained by using the hybrid function are in good agreement with [10] for case (1) and [12] for case (2) where the FEM with traditional Gaussian quadrature were used to carry out the calculation. 3.1. Comparison using COMSOL The COMSOL Model Wizard was used to create a station- ary, 2D Multiphysics model. The computation was carried out in a physical machine with Intel i7 64-bit processor with 8 GB of memory and SSD for storage, loaded with Windows 10 Op- erating system. The time taken for computation for a partic- ular Pr and various Ra is 20 seconds. The 2D geometry is mathematically extended to infinity in both directions along the z-axis, assuming no variation along that z − axis. A square with 1metre on all sides is built, inclined to specific angles for needed cases. Different components with different inclination angles are created in the same model. ‘Parameters’ nodes are present under global definitions, that are used for creating and Figure 3. 20 × 20 meshing of the domain defining the values of the non-dimensional temperature bound- ary conditions for hot, cold and adiabatic walls, Rayleigh Num- ber, Prandtl Number and the reference pressure. Coupling of the fields are at ease even with more than one physics interface. All applicable fields that can be used as in- puts in one physics interface is spontaneously found in the other physics interface’s Settings. Laminar flow (sp f ) takes care of fluid properties, initial conditions, volume force and the pres- sure of the fluid. Heat Transfer in Fluids (ht) handles the nature of the heat transfer and the temperature boundary conditions. The Mesh nodes are generated by meshing the domain, which enables the discretization of the geometry into mesh elements. The parameter Ra varied is listed for the given range using Aux- iliary sweep. The set-up is computed. COMSOL Multiphysics generates plot windows for displaying convergence results and isothermal contours and streamlines. The results from COM- SOL are compared with those of Mathematica for every bound- ary condition and Ra found to be in good agreement. 4. Conclusion In this paper, we considered a two dimensional, steady and laminar flow across the tilted square cavity under two different thermal boundary conditions: case (1) Isothermally cold side- walls of the cavity and the hot bottom wall kept parallel to the insulated top wall and case (2) Hot left wall, cold right wall, in- sulated top and bottom walls. In case 1, it is found that two ax- isymmetric flows dwell in the entire cavity, whereas, as the an- gle of inclination increases, the primary circulation cells occupy the entire cavity and the secondary circulation cells disappear completely. In case 2, only at Ra = 105, the streamlines are cir- cular only on enlarging inclination of the square cavity. In both the cases as inclination increases, buoyancy force along the hot wall gradually increases leading to stronger anticlockwise fluid circulations corresponding to the Rayleigh number. Results in- dicate that the streamlines and isotherms are strongly depen- dent on the inclination angle of the cavity. An attempt has been made to evaluate the definite integrals of the residual equations of the finite element equations using Gaussian quadrature with Hybrid functions of block-pulse function and Lagrange poly- nomial as orthogonal polynomials, and the results are found 238 Mariya & Prabhakar / J. Nig. Soc. Phys. Sci. 4 (2022) 231–241 239 to be in accordance with [10] for case (1) and [12] for case (2). With hybrid functions, desired accuracy can be achieved with the flexibility of choosing the orders of block pulse func- tion and Lagrange polynomial independently. Finite element results presented with Wolfram Mathematica were found to be in good agreement with COMSOL Multiphysics. The model can be applied for a field problem like plastic injection mould flow, porous media flow etc., by varying the variables such as fluid temperature, fluid velocity, mass weight and flow rate. Acknowledgments We thank the referees for the positive enlightening com- ments and suggestions, which have greatly helped us in making improvements to this paper. APPENDIX A.Notations used: X = x L ,Y = y L , (X and Y - distance in dimensionless form) U = uL α , Y = vL α , (U and V - velocity components in dimen- sionless form) θ = T − Tc Th − Tc (θ - temperature component in dimensionless form; Th&Tc-Temperature at hot and cold) P = pL3 ρα2 (P - pressure in the dimensionless form) Pr = ν α (Pr - Prandtl number) Ra = gβ(Th − Tc)L3 Pr ν (Pr - Rayleigh number) α = k ρC p (α - Thermal diffusivity) B.The shape functions of the transformed residuals from X − Y to s − t plane with reference to the local numbering in the Figure 2. N1 = (1 − 3s + 2s2)(1 − 3t + 2t2); N2 = (1 − 3s + 2s2)(4t − 4t2); N3 = (1 − 3s + 2s2)(−t + 2t2); N4 = (4s − 4s2(1 − 3t + 2t2); N5 = (4s − 4s2)(4t − 4t2); N6 = (4s − 4s2)(−t + 2t2); N7 = (−s + 2t2)(1 − 3t + 2t2); N8 = (−s + 2t2)(4t − 4t2); N9 = (−s + 2t2)(−t + 2t2); C.The nonlinear residual equations (9) - (11) are solved by the following procedure Grouping the co-efficient of ∑9 j=1 U j from all three residues s11 = ∫ D   9∑ j=1 U j N j  ∂N j∂X +  9∑ j=1 V j N j  ∂N j∂Y  NidXdY + γ ∫ D ∂Ni ∂X ∂N j ∂X dXdY + Pr ∫ D ( ∂Ni ∂X ∂N j ∂X + ∂Ni ∂Y ∂N j ∂Y ) dXdY (C.1) s12 = γ ∫ D ∂Ni ∂Y ∂N j ∂Y dXdY (C.2) s13 = 0 (C.3) Grouping the co-efficient of ∑9 j=1 V j from all three residues s21 = γ ∫ D ∂Ni ∂Y ∂N j ∂X dXdY (C.4) s22 = ∫ D   9∑ j=1 U j N j  ∂N j∂X +  9∑ j=1 V j N j  ∂N j∂Y  NidXdY + γ ∫ D ∂Ni ∂Y ∂N j ∂Y dXdY + Pr ∫ D ( ∂Ni ∂X ∂N j ∂X + ∂Ni ∂Y ∂N j ∂Y ) dXdY (C.5) s23 = −Ra Pr 9∑ j=1 θ j ∫ D N jdXdY (C.6) Grouping the co-efficient of ∑9 j=1 θ j from all three residues s31 = 0 (C.7) s32 = 0 (C.8) s33 = ∫ D   9∑ j=1 U j N j  ∂N j∂X +  9∑ j=1 V j N j  ∂N j∂Y  NidXdY + ∫ D ( ∂Ni ∂X ∂N j ∂X + ∂Ni ∂Y ∂N j ∂Y ) dXdY (C.9) Rewriting the governing equations, with the above substitu- tions, in matrix form R(1)i R(2)i R(3)i  = 9∑ j=1  si, j11 s i, j 12 s i, j 13 si, j21 s i, j 22 s i, j 23 si, j31 s i, j 32 s i, j 33   U j V j θ j  (C.10) { Ri } = [ ki j ] { U j } , 1 ≤ i ≤ 9 (C.11) The integrands, which leads to the element stiffness matrix, are functions of the global coordinates X and Y . The functions N j can be expressed in terms of the local coordinates s and t. Figure 2 gives an idea of co-ordinate transformation for the dis- cretized elements in the X − Y plane to the s − t plane. X = 9∑ j=1 x j N j and Y = 9∑ j=1 y j N j (C.12) The integrand contains not only functions but also deriva- tives with respect to the global coordinates (x j, y j). There- fore, ∂N j ∂X and ∂N j ∂Y must be related to ∂N j ∂s and ∂N j ∂t . Applying the shape functions to equation (C.10), the transformed finite element equation is obtained. Element stiffness matrices thus obtained are assembled using the element connectivity (Table 13) resulting in a global stiffness matrix. The global stiffness matrix is solved to get the thermal and the velocity components. On applying of the boundary condi- tions to the walls, matrix reduces to 9 nodes (internal), i.e., 27 239 Mariya & Prabhakar / J. Nig. Soc. Phys. Sci. 4 (2022) 231–241 240 Table 13. Element connectivity involving Nodes (N) and Elements (E) NE 1 2 3 4 1 1 3 11 13 2 2 4 12 14 3 3 5 13 15 4 6 8 16 18 5 7 9 17 19 6 8 10 18 20 7 11 13 21 23 8 12 14 22 24 9 13 15 23 25 Figure 4. Isotherms contour of 4 elements inclined φ = 15 ◦ ; Ra = 104; Pr = 0.71 obtained from Mathematica unknowns and equations, for the entire domain discretized to for 2 × 2. The temperature contours of 4 elements for a tilted square cavity with Ra = 104; Pr = 0.71 is displayed in Figure 4. The values of the velocity (U ,V) and thermal components (θ) for the internal nodes of the discretized domain with 2 × 2 (4 elements) are given below. U7 = −0.0000228252, V7 = −0.0000176029,θ7 = 0.414016, U8 = 4.69382E-6, V8 = −1.88969E-6,θ8 = 0.208614, U9 = 0.0000131575, V9 = −1.06333E-6,θ9 = 0.121948, U12 = −0.0000119197, V12 = −0.000053752,θ12 = 0.575621, U13 = 0.0000149714, V13 = −0.0000638262,θ13 = 0.271557, U14 = 0.0000157566, V14 = −0.000030133,θ14 = 0.171982, U17 = 1.58957E-6, V17 = −0.000024564,θ17 = 0.414017, U18 = 1.69972E-6, V18 = −0.0000310296,θ18 = 0.208615, U19 = −1.65325E-6, V19 = −0.0000198665,θ19 = 0.121949. The discretization was gradually increased until the con- vergence was obtained at 20 × 20. References [1] L. Viorel & N. Anisoara, “Numerical modelling of fluid flow and heat transfer in a corrugated channel for heat exchanger applications”, Proce- dia Manufacturing 22 (2018) 634. [2] A. Omri, J. Orfi & S.B. Nasrallah, “Natural convection effects in solar stills”, Desalination 183 (2005) 73. [3] G.D. Mey, M. Wojcik, J. Pilarski, M. Lasota, J. Banaszczyk, B. Vermeer- sch, A. Napieralski & M.D. Paepe, “Chimney effect on natural convection cooling of a transistor mounted on a cooling fin”, Journal of Electronic Packaging 131 (2009) 14501. [4] M. Sabour, M. Ghalambaz, & A Chamkha, “Natural convection of nanofluids in a cavity: criteria for enhancement of nanoflu- ids”,International Journal of Numerical Methods for Heat and Fluid Flow 27 (2017) 1504. [5] K. Wisam, B.Hussama, K. Khanafera, H.J. Salema & Sheard GJ, “Natural convection heat transfer utilizing nanofluid in a cavity with a periodic side-wall temperature in the presence of a magnetic field”,International Communications in Heat and Mass Transfer 104 (2019) 127. [6] N. Nabipour, M. Babanezhad, A.T.Nakhjirin, & S. Shirazian, “Prediction of Nanofluid Temperature Inside the Cavity by Integration of Grid Par- tition Clustering Categorization of a Learning Structure with the Fuzzy System”,International Communications in Heat and Mass Transfer 5 (2020) 3571. [7] A.J. Chamkha, “Hydromagnetic combined convection flow in a vertical lid-driven cavity with internal generation or absorption”,An International Journal of Computation and Methodology 41 (2002) 529. [8] C.S.N. Azwadi, M.Y.M Fairus & Syhrullail.S, “Virtual Study of Natural convection Heat transfer in an Inclined square cavity”,Journal of Applied Sciences 10 (2010) 331. [9] M.K Das & K.S.K Reddy, “Conjugate natural convection heat transfer in an inclined square cavity containing a conducting block”,International Journal of Heat and Mass Transfer 49 (2006) 4987. [10] A.K. Singh, S. Roy, T. Basak, “Analysis of Bejan’s heatlines on vi- sualization of heat flow and thermal mixing in tilted square cavi- ties”,International Journal of Heat and Mass Transfer 55 (2012) 2965. [11] J.K. Mariya Helen Mercy & V. Prabhakar, “Study of fluid flow inside closed cavities using computational numerical methods”,International Journal for Simulation and Multidisciplinary Design Optimization 12 (2012) 1. [12] J. Rasoul & P. Prinos, “Natural convection in an inclined enclo- sure”,International Journal of Numerical Methods for Heat & Fluid Flow 12 (1997) 438. [13] F. Xu, J.C. Patterson, C. Lei, “Natural convection in a differentially Heated cavity with a square obstruction on the sidewall”,Australian Jour- nal of Mechanical Engineering 4 (2015) 77. [14] Bendaraa et al, “Numerical Study of Natural Convection in Square Cav- ity Using Lattice Boltzmann Method: Obstacles Effect”, IOP Conf. Ser: Mater. Sci. Eng 948 (2020) 1. [15] J. Vierendeels, B. Merci & E.Dick, “Benchmark solutions for the natural convective heat transfer problem in a square cavity with large horizontal temperature differences”,International Journal of Numerical Methods for Heat & Fluid Flow 13 (2003) 1057. [16] R.C. Mohapatra, “Study on Natural Convection in a Square Cavity with Wavy right vertical wall Filled with Viscous Fluid”,Journal of Mechanical and Civil Engineering 14 (2017) 32. [17] J.N. Reddy, An Introduction to the Finite Element Method, United States of America, McGraw-Hill, New York (1993). [18] Victor Oboni Atabo, Solomon Ortwer Adee, “A New Special 15-Step Block Method for Solving General Fourth Order Ordinary Differential Equations”,J. Nig. Soc. Phys. Sci.3 (2021) 308. [19] Mark I. Modebei, Olumide O. Olaiya, Ignatius P. Ngwongwo, “Compu- tational study of a 3-step hybrid integrators for third order ordinary dif- ferential equations with shift of three off-step points”,J. Nig. Soc. Phys. Sci.3 (2021) 459. [20] V. J. Shaalini, S. E. Fadugba, “A New Multi-Step Method for Solving Delay Differential Equations using Lagrange Interpolation”,J. Nig. Soc. Phys. Sci.3 (2021) 159. 240 Mariya & Prabhakar / J. Nig. Soc. Phys. Sci. 4 (2022) 231–241 241 [21] Friday Obarhua, Oluwasemire John Adegbor, “An Order Four Continuous Numerical Method for Solving General Second Order Ordinary Differen- tial Equations”,J. Nig. Soc. Phys. Sci.3 (2021) 42. [22] Olumide O. Olaiya, Rasaq A. Azeez, Mark I. Modebei, “Efficient Hy- brid Block Method For The Numerical Solution Of Second-order Partial Differential Problems via the Method of Lines”,J. Nig. Soc. Phys. Sci.3 (2021) 26. [23] H.R. Marzban, H.R. Tabrizidooz & M. Razzaghi, “Hybrid functions for nonlinear initial value problems with applications to Lane-Emden type equations”,Physics Letters A 372 (2008) 5883. [24] H.R. Marzban, S.M. Hoseini & M. Razzaghi, “Solution of Volterra’s pop- ulation model via block-pulse functions and lagrange interpolating poly- nomials”,Math Methods Appl Sci, 32 (2009) 127. [25] Siraj-Ul-Islam, I. Aziz & F. Haq, “A comparative study of numerical in- tegration based on Haar wavelets and hybrid functions”,Computers and Mathematics with Applications, 59 (2009) 2026. [26] I. Aziz, Siraj-Ul-Islam, & W. Khan, “Quadrature rules for numerical in- tegration based on haar wavelets and hybrid functions”,Computers and Mathematics with Applications, 61 (2011) 2770. [27] G. Uma, V. Prabhakar & S. Hariharan, “Numerical integration using hybrid of block-pulse functions and Lagrange polynomial”,International Journal of Pure and Applied Mathematics, 106 (2016) 33. [28] H.R. Marzban, H.R. Tabrizidooz & M.Razzaghi, “Solution of variational problems via hybrid of block-pulse and Lagrange interpolating”,ET Con- trol Theory Appl, 3 (2009) 1363. [29] Larry J.Segerlind, Applied finite element analysis, United States of Amer- ica, John Wiley & Sons, New York (1984). 241