## Materials and Methods  
### Mathematical Formulation  
A pair of thin staggered vertical plates of length \(H\) and infinite width are suspended in free air at a distance \(W\) to form a duct, as sketched in Figure 1, where the reference Cartesian coordinate system is also indicated. The edges of the plates are misaligned so that \(L\) is the portion of one plate facing the other, which means that when \(L=H\) the plates are located face to face and when \(L=0\) the plates are completely staggered. The front side of each plate is heated at a uniform temperature \(t_{h}\) and the surrounding undisturbed fluid reservoir is maintained at a lower uniform temperature \(t_{\infty}\). The rear side of each plate is assumed to be perfectly insulated. The resulting buoyancy-induced flow is considered to be two-dimensional, laminar and incompressible, with constant physical properties. Additionally, the buoyancy effects on the momentum transfer are taken into account through the customary Boussinesq approximation, whereas viscous dissipation and pressure work are neglected. Figure 1: Sketch of the geometry, coordinate system and computational domain. Upon incorporating these hypotheses into the equations of continuity, momentum and energy, the following set of governing equations expressed in dimensionless form is obtained:  
\[\
abla\cdot\mathbf{V}=0 \tag{1}\] \[\frac{\partial\mathbf{V}}{\partial\tau}+(\mathbf{V}\cdot\
abla) \mathbf{V}=-\
abla P+\
abla^{2}\mathbf{V}-\frac{Ra}{Pr}T\frac{\mathbf{g}}{g}\] (2) \[\frac{\partial T}{\partial\tau}+(\mathbf{V}\cdot\
abla)T=\frac{1 }{Pr}\
abla^{2}T \tag{3}\]  
where \(\tau\) is the dimensionless time normalized by \(H^{2}/\
u\), \(\mathbf{V}\) is the dimensionless velocity vector possessing horizontal and vertical components \(U\) and \(V\) normalized by \(\
u/H\), \(T\) is the dimensionless temperature excess over the uniform temperature of the undisturbed fluid reservoir normalized by the temperature difference \((t_{h}-t_{\infty})\), \(P\) is the dimensionless sum of the thermodynamic and hydrostatic pressures normalized by \(\rho\
u^{2}/H^{2}\), \(\mathbf{g}\) is the gravity vector, \(Pr=\
u/\alpha\) is the Prandtl number and \(Ra_{H}\) is the Rayleigh number defined as the following:  
\[Ra_{\mathrm{H}}=\frac{g\beta(t_{h}-t_{\infty})H^{3}}{\alpha\
u} \tag{4}\]  
in which \(\
u\) is the kinematic viscosity, \(\rho\) is the mass density, \(\alpha\) is the thermal diffusivity and \(\beta\) is the coefficient of volumetric thermal expansion of the fluid. The related boundary conditions are \(\mathbf{V}=0\) and \(T=1\) at the heated side of any plate and \(\mathbf{V}=0\) and \(\partial T/\partial n=0\) at the thermally insulated rear side of any plate. \(n\) denotes the normal direction relative to the plate surface and \(\mathbf{V}=0\) and \(T=0\) at very large distance from the plates. The integration domain is taken as a two-dimensional square ABCD containing the plates and extending sufficiently far from them. Such an integration domain is filled with a non-uniform Cartesian grid possessing a higher concentration of grid lines near the plates and a lower concentration of grid lines far away from the plates. Notice that due to the infinitesimal thickness assumed for the plates, each grid node situated on any plate actually consists of a pair of coincident yet distinct grid nodes. One belongs to one side of the plate and the other belongs to the opposite side of the plate such that the same spatial location is permitted to possess simultaneously two different temperatures, with one for each side of the plate. The boundary conditions required for the numerical solution of the governing Equations (1)-(3) have to be specified at each of the boundary lines that enclose the two-dimensional integration domain defined above. Once these lines are set sufficiently far away from the plates, the motion of the fluid which enters or leaves the computational domain can be reasonably assumed to normally occur relative to them.