**3. Numerical Procedure**

The dimensionless governing Equations (6)–(9) together with the boundary conditions (10) and (11) have been solved numerically using the commercial software tool, which is known as COMSOL Multiphysics (Version 5.5, COMSOL Inc., Stockholm, Sweden). The software employs the Galerkin finite element method, which enforces the orthogonality of residuals to all basis functions in a basis. In Galerkin formulation, weighting functions are chosen to become identical to basis functions [36]. In this paper, we have employed a segregated and parallel direct (Pardiso) solver to solve those equations. As convergence criteria, 10−<sup>6</sup> has been chosen for all dependent variables.
