## 3 Numerical Method and Validation  
The dimensionless governing equations Equations (17)-(20) including the selected boundary conditions Equations (21)-(24) are determined via the Galerkin weighted residual finite-element technique. Let \(\
abla=\left(\frac{\partial}{\partial X},\frac{\partial}{\partial Y},\frac{ \partial}{\partial Z}\right)\), and \(\mathbf{V}=\left(\text{U},\,\text{V},\,\text{W}\right)\), the finite element analysis concerning the momentum equation Equation (18) in vectorial form is shown by the following steps. The weak formulation regarding Equation (18) through multiplying this equation with an internal region of (\(\Phi\)) and integrating that across the computational area which remains discretized toward small triangular components as revealed within Figure 2. The resulting weak formulation is reached as  
\[\int_{\Omega}\Phi_{i}\mathbf{V}^{k}\left(\frac{\partial\mathbf{V }^{k}}{\partial X}+\frac{\partial\mathbf{V}^{k}}{\partial Y}+\frac{\partial \mathbf{V}^{k}}{\partial Z}\right)\mathrm{d}\mathbf{X}\mathrm{d}\mathbf{Y} \mathrm{d}Z=-\int_{\Omega}\Phi_{i}P^{k}\left(\frac{\partial}{\partial X}+\frac {\partial}{\partial Y}+\frac{\partial}{\partial Z}\right)\mathrm{d}\mathbf{X} \mathrm{d}\mathbf{Y}\mathrm{d}Z\] \[+\frac{\rho_{f}}{\rho_{nf}}\frac{\mu_{nf}}{\mu_{f}}\int_{\Omega} \Phi_{i}\mathbf{V}^{k}\left(\frac{\partial^{2}}{\partial X^{2}}+\frac{\partial ^{2}}{\partial Y^{2}}+\frac{\partial^{2}}{\partial Z^{2}}\right)\mathrm{d} \mathbf{X}\mathrm{d}\mathbf{Y}\mathrm{d}Z+\frac{(\rho\theta)_{nf}}{\rho_{nf} \beta_{f}}\frac{Ra}{\text{Pr}}\int_{\Omega}\Phi_{i}\theta^{k}\mathrm{d} \mathbf{X}\mathrm{d}\mathbf{Y}\mathrm{d}Z.\]  
Selection of the interpolation functions for providing an approximation for the velocity distribution, pressure, and temperature distribution as  
\[\mathbf{V}\approx\sum_{j=1}^{m}\mathbf{V}_{j}\Phi_{j}(X,Y,Z),\quad P\approx \sum_{j=1}^{m}P_{j}\[PERSON]_{j}(X,Y,Z),\quad\theta\approx\sum_{j=1}^{m}\theta_{j }\[PERSON]_{j}(X,Y,Z).\]  
The nonlinear residual equation for Equation (18) that is obtained from the [PERSON] weighted residual finite-element method is  
\[R_{i} =\sum_{j=1}^{m}\mathbf{V}_{j}\int_{\Omega}\left[\left(\sum_{j=1}^ {m}\mathbf{V}_{j}\[PERSON]_{j}\right)\frac{\partial\[PERSON]_{j}}{\partial X}+\left( \sum_{j=1}^{m}\mathbf{V}_{j}\[PERSON]_{j}\right)\frac{\partial\[PERSON]_{j}}{\partial Y }+\left(\sum_{j=1}^{m}\mathbf{V}_{j}\[PERSON]_{j}\right)\frac{\partial\[PERSON]_{j}}{ \partial Z}\right]\[PERSON]_{i}\mathrm{d}\mathbf{X}\mathrm{d}\mathbf{Y}\mathrm{d}Z\] \[+\frac{\rho_{f}}{\rho_{nf}}\frac{\mu_{nf}}{\mu_{f}}\sum_{j=1}^{m} \mathbf{V}_{j}\int_{\Omega}\left[\frac{\partial\[PERSON]_{i}}{\partial X}\frac{ \partial\[PERSON]_{j}}{\partial X}+\frac{\partial\[PERSON]_{j}}{\partial Y}\frac{ \partial\[PERSON]_{j}}{\partial Y}+\frac{\partial\[PERSON]_{i}}{\partial Z}\frac{ \partial\[PERSON]_{j}}{\partial Z}\right]\mathrm{d}\mathbf{X}\mathrm{d}\mathbf{Y} \mathrm{d}Z\] \[+\frac{(\rho\theta)_{nf}}{\rho_{nf}\beta_{f}}\frac{Ra}{\text{Pr }}\int_{\Omega}\left(\sum_{j=1}^{m}\theta_{j}\Phi_{i}\right)\Phi_{i}\mathrm{d} \mathbf{X}\mathrm{d}\mathbf{Y}\mathrm{d}Z,\]  
where the superscripts \(k\), \(i\), \(j\), and \(m\) are approximate index, residual number, node number, and iteration number, respectively. For the intention of simplifying these nonlinear expressions into the momentum equations, any Newton-Raphson iteration algorithm is employed. The current convergence of the solution exists when the relative error toward variables remains according to the following convergence criteria,  
\[\left|\frac{\Gamma^{m+1}-\Gamma^{m}}{\Gamma^{m+1}}\right|\leq 10^{-5}.\]  
Toward the aim of approving the existing numerical data, the outcomes of the current study are compared to the outcomes described by [PERSON] et al. [37] concerning the mode of free convection inside a square hollow cavity of heated of sides, as designated in Figure 3. Besides, Figure 4a shows a comparison between the current outcomes and the experimental arrangements of [PERSON] et al.