Free convection inside a porous square cavity with convective boundary condition using spline functions

In this paper, we consider the problem of free convection in a square cavity filled with a porous medium with convective boundary condition on the left wall of the cavity. We first transform the governing equations transformed in terms of dimensionless variables and then solve them numerically using a cubic spline colocation method. We discuss the effects of two very important parameters, the Biot (Bi) and Rayleigh (Ra) numbers. We perform a comparison of the average Nusselt number Nu at the hot wall with results from the open literature. We can notice that the comparison is very good, which gives us strong confidence that this cubic spline collocation method works very efficiently for such problems. We also state that the present problem has not been considered before by any researcher.

ample, in heat exchange where the conduction in the solid tube wall is greatly influenced by the convection in the fluid flowing over it (see Chaudhary and Jain [14]).
We also mention some very good papers on efficient methods for such problems: Abdelwahed et al. [15] on reconstruction of Tesla micro-valve using topological sensitivity analysis, Abdelwahed and Chorfi [16] on the convergence analysis of a time-dependent elliptic equations with discontinuous coefficients, and Chorfi et al. [17] on the analysis of a geometrically selective turbulence model.
The objective of the present analysis is investigating the influence of the thermal boundary condition in a square cavity filled with a porous medium using a spline colocation method. We determine numerical solutions for the streamlines, isotherms, and local and average Nusselt numbers.

Mathematical formulation
We consider the free convection in a square cavity filled with a fluid-saturated porous medium, as it is shown in Fig. 1, where the Cartesian coordinates are x and y, and L is the wall length. The temperatures at the left and right sides of the walls are considered constant, and the bottom and top walls are thermally insulated. We also assume that the outside left side surface of the cavity is heated by convection from a hot fluid at constant temperature T f , which provides a heat transfer coefficient h f (see Aziz [18]), namely The homogeneity and the thermal equilibrium of the porous medium are assumed. Then the steady conservation of mass, Darcy's law with the Oberbeck-Boussinesq approximation employed, and the energy equations are (see Nield and Bejan [10]) In these equations, v = (u, v) is the velocity vector, T is the fluid temperature, p is the fluid pressure, k is the thermal conductivity, ∇ 2 is the Laplace operator, and the physical meanings of the other quantities are mentioned in the Abbreviations list.
where α m = k m /(ρc) f is the effective thermal diffusivity. Introducing the dimensionless variables and substituting (2.8) into equations (2.6)-(2.7), we obtain: where Ra = gKρ 0 β T/(α m μ) is the Rayleigh number, and T = T h -T c . Then we complete the mathematical model with the following corresponding boundary conditions: where Bi = Lh f /k f is the Biot number. Note that for large values of Bi( 1), from (2.11) it follows that θ (0) = 1 (isothermal left wall of the cavity). The physical quantities of interest are the local and average Nusselt numbers Nu l and Nu l . The local Nusselt number at the left wall is given by 12) and the mean Nusselt number is

Numerical method
The governing equations (2.9) and (2.10) coupled with the boundary conditions (2.11) are approximated numerically using the SADI (spline-alternating-direction-implicit) method, a combination of spline collocation and standard ADI (alternating-direction-implicit) methods. Each unknown is approximated by a cubic spline function. This is an efficient method, having good accuracy and low cost in implementation, as it leads to a governing linear tridiagonal system. Using the fictitious time technique, we write equations (2.9) and (2.10) in the form in the x-direction, with m ,ij and M ,ij denoting the spline approximations of the firstand second-order partial derivatives with respect to x. Then, in the y-direction, we write where l ,ij and L ,ij are the spline approximations of ( y ) ij and ( yy ) ij . In these equations, i, j = 0, N refer to the mesh nodes, n denotes the false time step, and ij are the cubic spline approximations of the functions θ or ψ at the grid nodes. The coefficients F ij , R ij , and Q ij at the two steps are given in Tables 1 and 2. By the usual cubic spline collocation relations (see Pu and Kahawita [19]) we transform equations (3.1) and (3.2) to and respectively, where ϕ can be the function or one of its derivatives (m or M for equation (3.3) and l or L for equation (3.4)). Obtaining two more equations from the boundary conditions, we get a closed tridiagonal, diagonally dominant linear system, which can be solved efficiently by the Thomas algorithm. For details on the mathematical computations, see Rubin and Graves [20] or Pu and Kahawita [19]. The computations are iterated until a steady-state solution is reached, that is, where n is the number of iterations. When a steady-state solution is reached, the mean Nusselt number is approximated by the trapezoidal rule. The procedure is described in more detail by Micula and Pop [21]. At the first step, in the x-direction, from time n t to (n + 1/2) t, discretizing equations (2.6)-(2.7), for every i = 1, . . . , N -1 and fixed j ∈ {0, 1, . . . , N}, we write and with the coefficients given in Table 1. These equations are then put in the form and Also, combining the boundary conditions (2.11) with relations between cubic spline functions and their derivatives, we get one more equation The derivation of coefficients A θ,0j , B θ,0j , and D θ,0j is given in the Appendix. Equations For the second step, in the y-direction, from time (n + 1/2) t to (n + 1) t, we proceed similarly. For any fixed i ∈ {0, 1, . . . , N} and for every j = 1, . . . , N -1, using the coefficients from Table 2, we get equations and Again, from the boundary conditions (2.11) we get two more equations for l n+1 θ,ij , l n+1 θ,i0 = 0, l n+1 θ,iN = 0, (3.16) and two for ψ n+1 ij , Then equations (3.14) and (3.16) form a closed linear system for the unknowns {l n+1 θ,ij } j=0,N (again, from these the values of {θ n+1 ij } i=0,N and {L n+1 θ,ij } j=0,N can also be found). Similarly, equations (3.15) and (3.17) give an (N + 1) × (N + 1) tridiagonal linear system for the unknowns {ψ n+1 ij } j=0,N .

Numerical results and discussion
We apply the SADI method to our problem taking the values Ra = 10, 100, 1000 and Bi = 0.5, 5, 50, 100. The method is validated by comparing the values of Nu we obtained with other well-known results in the literature (see Table 3). We use uniform meshes of N × N nodes for N = 100, 200, 300, and 400. A sensitivity to the grid size test is performed, and the results are seen in Table 4. As the values in   For the maximum values of the parameters (Ra = 1000 and Bi = 100), we also performed calculations using nonuniform grids. We used a nonuniform mesh obtained by keeping the spacing ratio from the wall to the center h i+1 /h i a constant (strictly greater than 1), which depends on the desired number of nodes. In this case, we produced a symmetric grid with h min ≈ 5 × 10 -5 (near the walls) and h max ≈ 0.02 (at the center). The variation of Nu with the number of nodes for a nonuniform grid is shown in Table 6. As it can be seen, there are insignificant differences from smaller to larger numbers of nodes. Note that the values obtained when a nonuniform mesh was used are slightly higher than those obtained by a uniform grid at the same number of nodes. That was also evident in the simple case of a hot wall, θ = 1, on x = 0 (see Table 3).

Conclusions
We used a cubic spline collocation scheme to obtain the numerical solution of the problem of free convection inside a porous square cavity with convective boundary condition. This is an efficient method as it leads to a sparse (tridiagonal) linear system that can be solved faster using the Thomas algorithm. The SADI method is especially useful when  Table 3 shows. As the numerical results and graphs show, when the values of Bi increase ( 1), the boundary condition on the left wall approaches the case θ = 1 (isothermal left wall of the cavity).
The Biot number enhances the release of heat energy, substantially cooling the system. This analysis may certainly have practical impact on obtaining the conditions needed to design an integrated system to achieve optimal functionality.