SPH modeling of natural convection around a heated horizontal cylinder: A comparison with experiments

The equations describing the fluid motion and heat transfer around a horizontal cylindrical heat source are the continuity, the Navier-Stokes, and the energy balance equations. In Lagrangian form these equations are:

where $d/dt$ is the material time derivative, $\rho$ is the density, ${\bf v}$ is the fluid velocity vector, $p$ is the pressure, $\nu=\eta/\rho$ is the kinematic viscosity coefficient and $\eta$ is the dynamic viscosity of the fluid, $T$ is the temperature, $c_{p}$ is the specific heat at constant pressure, and $\kappa$ is the thermal conductivity. Since at low ${\rm Ra}$-values heat generation due to viscous dissipation is negligible [27], only the thermal conduction term is retained in the energy balance equation (3). This is true because the total entropy generation due to thermal dissipation is by far much more intense than that due to viscous dissipation at ${\rm Ra}<10^{3}$ [27]. In the Boussinesq approximation, the body force ${\bf F}$ on the right-hand side of Eq. (2) is given by [37]

where ${\bf g}$ is the gravitational acceleration, $\beta$ is the coefficient of thermal expansion, and $T_{\rm r}$ is a reference temperature. In most applications this is the fluid temperature at room conditions. For natural convection flow in the Boussinesq approximation, the two relevant dimensionless numbers are: the Prandtl number

where $\alpha=\kappa/(\rho c_{p})$ is the thermal diffusivity, and the Rayleigh number

where $D$ is a characteristic length scale and $\Delta T$ is the leading temperature difference between the heat source ($T_{s}=T_{\rm bath}$) and the surrounding fluid. In particular, the value of Ra quantifies the ratio of buoyancy over viscosity in the fluid for a given Prandtl number.

Equations (1)-(3) are closed by a pressure-density relation of the form [38]

where $B=c_{0}^{2}\rho_{\rm r}/\gamma$, $\rho_{\rm r}$ is a reference density, $\gamma=7$ for liquids and $\gamma=1.4$ for gases, and $c_{0}$ is a numerical speed of sound which is chosen to be 10 times higher than the maximum fluid velocity in order to keep the density fluctuations below 1% and satisfy the incompressibility condition [39]. The reference density $\rho_{\rm r}$ in Eq. (7) is allowed to change with the temperature according to the following expression

where $T_{\rm f}$ and $T_{0}$ are, respectively, the final and initial temperatures in a numerical timestep $\Delta t$.

Equations (1)-(3) are solved numerically using a variant of the open-source code DualSPHysics [40], which relies on SPH theory [39, 41]. The temporal variation of the fluid density at the position ${\bf x}_{a}$ of particle $a$ is calculated according to the SPH representation [39]

where ${\bf v}_{ab}={\bf v}_{a}-{\bf v}_{b}$, $W_{ab}=W(|{\bf x}_{ab}|,h)$ is the kernel interpolation, ${\bf x}_{ab}={\bf x}_{a}-{\bf x}_{b}$, $h$ is the width of the kernel (or smoothing length), and the sum is over $n$ neighbours around particle $a$. This form is Galilean invariant and is particularly convenient in the presence of free surfaces or density discontinuities. For the momentum equation (2), we use the following representation

where ${\hat{\rho}}_{ab}=(\rho_{a}+\rho_{b})/2$ and $\epsilon^{2}=0.01h^{2}$ is a small correction factor used to avoid singularities when $|{\bf x}_{ab}|\ll 1$. Following Cleary [22], the SPH form of the energy equation is given by

where $T_{ab}=T_{a}-T_{b}$. This form ensures that the heat flux is automatically continuous across material interfaces. In its original form derived by Cleary [22] with $\kappa\to 2\kappa_{a}\kappa_{b}/(\kappa_{a}+\kappa_{b})$, it allows multiple materials with substantially different conductivities and specific heats to be accurately modelled. In a Lagrangian frame of reference, the SPH particles are moved according to

which must be solved simultaneously with Eqs. (9)-(11).

Since direct evaluation of second-order derivatives of the kernel function is not required, a low-order Wendland C² function is used as the interpolation kernel [42]

for $0\leq q<1$ and zero otherwise, where $q=|{\bf x}-{\bf x}^{\prime}|/h$. Wendland functions have positive Fourier transforms and so they can support arbitrarily large numbers of neighbours without suffering from a close pairing of particles [42, 43, 44]. In addition, Wendland functions are very reluctant to allow for particle motion on a sub-resolution scale, and in contrast to most commonly used kernels they maintain a very regular particle distribution, even in highly dynamical tests as occurs in turbulent regimes [45].

The time integration of Eqs. (9)-(12) is performed using the Verlet algorithm provided by DualSPHysics, where the density, velocity, position, and temperature of particle $a$ are advanced from time $t^{n}$ to $t^{n+1}=t^{n}+\Delta t$ using the following steps

Numerical coupling of the discrete SPH equations is ensured during the evolution by alternating the above steps by the form

The time step $\Delta t=t^{n+1}-t^{n}$ to ensure stability of the above scheme is calculated as the minimum between the following timesteps:

where $c_{a}$ is the speed of sound of the fluid.

Equations (1)-(3) are numerically integrated in a two-dimensional domain. A horizontal cylinder of diameter $D=6$ mm is centered in a square plane of the same size as the square cross-section of the experimental enclosure of Fig. 3. No-slip boundary conditions are applied at the walls of the square region (i.e., ${\bf v}_{w}={\bf 0}$) and on the surface of the circular source pipe (i.e., ${\bf v}_{s}={\bf 0}$) as shown in Fig. 1. Isothermal boundary conditions are implemented at the walls of the cavity. This is a more accurate representation of the experimental reality, because it takes into account the finite flux of heat through the walls. Initially, the square region is filled with 508308 uniformly spaced particles at rest, where the fluid is represented by particles that carry the properties of glycerin at $20^{\circ}$C, i.e., $\rho=1264$ kg m^-3, $c_{p}=2368$ J kg^-1 K^-1, $\kappa=0.286$ W m^-1 K^-1, $\eta=1.519$ kg m^-1 s^-1, and $\beta=5.97\times 10^{-4}$ K^-1. These quantities yield the Prandtl number value ${\rm Pr}=12671.05$. With this spatial resolution the inter-particle separation distance is 0.7 mm in both the $x$- and $y$-directions, while the smoothing length is set to $\approx 0.99$ mm. Incompressibility of the fluid is guaranteed by setting $c_{0}=50$ m s^-1 in Eq. (7).

No-slip boundary conditions are implemented at the walls of the vessel using the method of dynamic boundary particles developed by Crespo et al. [46]. In this method, a linear distribution of uniformly-spaced particles is placed at the walls of the enclosure, with separations of $\approx h/1.42$. A second line of uniformly-distributed particles is placed outside the computational domain at distances $\Delta x/2$ and $\Delta y/2$ from the wall particles so that they are arranged in a staggered grid. This external particles are used to cope with the problem of kernel deficiency outside the computational domain. The wall particles are updated using the same loop as the inner fluid particles and so they are forced to satisfy Eqs. (7) and (9). However, they are not allowed to move according to Eqs. (10) and (12) so that their initial positions and velocities (${\bf v}_{w}={\bf 0}$) remain unchanged in time. In this way, the presence of the wall is modelled by means of a repulsive force, which is derived from the source term of the momentum Eq. (10) and includes the effects of compressional, viscous, and gravitational forces. This force is exerted by the wall particles on the approaching fluid particles only when the latter get closer than a distance $d=2h$ from the wall. Across the wall, a particle $a^{\prime}$ is assigned a mass $m_{a^{\prime}}=m_{w}$, and a density $\rho_{a^{\prime}}=\rho_{w}$, while in order to assure Neumann boundary conditions for the pressure $p_{a^{\prime}}=p_{w}$ (where the subscript $w$ indicates wall properties). With a number of 508308 fluid particles, 2917 linearly distributed particles are needed outside the computational domain, giving a total number of 511225 SPH particles. Isothermal wall boundary conditions are specified by setting $T_{a^{\prime}}=2T_{w}-T_{a}$ [25], where the wall temperature, $T_{w}$, is initially assumed to be in thermal equilibrium with the glycerin at $20^{\circ}$C and $T_{a}$ is the temperature of fluid particle $a$ when it is at a distance $<2h$ from the wall.

The circle at the center of the square cavity in Fig. 1, representing the surface of the cylindrical pipe, is also made of fixed boundary particles carrying a a temperature $T_{s}$. Six model calculations are considered where the value of $T_{s}$ is varied above and below the initial fluid temperature $T_{0}=20^{\circ}$C (293.15 K).

A slender acrylic enclosure of height $H=0.5$ m, length $L=0.5$ m, and width $W=0.025$ m, filled with glycerin (Pr$=12671.05$ at $20^{\circ}$C) as the working fluid, was employed for the experiments. The walls of the enclosure thus formed were 6 mm thick to prevent possible optical distortions (Figs. 2 and 3). A horizontal pipe with a circular cross section of diameter $D=0.006$ m spanned the width of the rectangular cavity. Its symmetry axis passed through the points ($H/2$,$L/2$) located on the front and back walls of the cavity, as shown in Fig. 3. The pipe was part of a closed circuit which allowed the circulation of water. The temperature of water was tightly controlled by means of a thermal bath within an accuracy interval of $\pm 0.001^{\circ}$C. Such an arrangement enabled an approximately uniform and steady temperature difference ($\Delta T$) between the horizontal pipe and the glycerin within the cavity. The temperature difference was measured every four seconds with five type $T$ thermocouples (Tc), where four of them were installed on the back wall with a perpendicular orientation with respect to the pipe, while the other one was attached directly to the pipe’s wall (see Fig. 3). The measurements from the thermocouples were acquired by means of a Digi Sense Scanning Thermometer Cole Parmer (Mod. 92000-00) temperature scanner.

The hydrodynamical characteristics of the flow were determined with a Particle Image Velocimetry (PIV) technique, which enabled an accurate spatial and temporal resolution of the velocity field. All measurements were carried out on a plane, perpendicular to the axis of symmetry of the pipe, generated with a laser light sheet. A 1200 mJ pulsating Nd: YAG Laser, with a pulse duration of 4 ms and a wavelength varying from 1064 nm to 582 nm, was used in this study. Images were recorded with a Hisense MKll camera (with an image resolution of $1260\times 1024$ pixels) equipped with a 60 mm Nikon lens and a 532 nm filter. The glycerin was seeded with polyamide (spherical) tracer particles, with a mean diameter of 20 $\mu$m and relative density $\rho_{\rm glycerin}/\rho_{\rm water}\approx 1$. The Dantec Dynamics Studio software was employed to process the acquired images.

The cell was filled with a prepared mixture of glycerin and tracer particles. The seeding concentration was optimized to yield a high statistical correlation for the displacement of particles in the domain. To achieve a detection probability of at least 90%, the mean number of tracers per interrogation window was set to be greater than 15 [47]. The laser intensity and exposure were adjusted to generate sufficient light scattering by the tracers and guarantee a good signal quality and flow traceability. The cameras were focused and image calibration was attained using a ruler placed in the center plane illuminated by the laser source. In the experiments, the temperature difference was defined as $\Delta T=T_{\rm bath}-T_{\rm glycerin}$. As it is common practice in most natural convection experiments, we first considered heating the horizontal pipe with $\Delta T$ values between 0.1 and 10 K. We then considered the opposite case, where the pipe was cooled with temperature differences in the interval ($-10$ K,$-0.1$ K). Each experiment was repeated three times to check its reproducibility.

Initially, a two-valve bypass prevented the flow of water inside the pipe. The water was thus forced to circulate through the bath, until it acquired a state of thermal equilibrium at $T_{\rm bath}$. At this point, both valves were opened to establish the flow inside the pipe and initiate the heat transfer process with the glycerin. All experiments were performed at a room temperature of $20^{\circ}$C.

The experimental reproducibility was assessed by repeating all experiments three times. Figure 4 shows steady-state profiles of the velocity projected along the line contained in the $z=0$ plane that passes through the center of the heat source. In particular, the plot shows the velocity as a function of height for the case when $\Delta T=10$ K and Ra $=242$ for three different measurements starting from identical conditions, namely Test 1, Test 2, and Test 3. The vertical velocity component is very sensitive to small changes in the initial conditions and, therefore it is a convenient parameter to measure the reproducibility of the experiments. The form of the profiles is consistent with a numerical prediction by Lauriat and Desrayaud [3] for laminar buoyancy-induced flows above a horizontal thin cylinder immersed in an air-filled vessel. A non-zero velocity exists close to the cylindrical surface as expected when the pipe behaves as a heat source. In terms of the vertical velocity component the differences between the three experimental results are small throughout the entire profiles. The maximum velocity for Test 1 is $v_{\rm max}\approx 5.24\times 10^{-4}$ m s^-1 and occurs at a height $y\approx 0.344$ m above the source. For comparison, the maximum velocities for Tests 2 and 3 are, respectively, $v_{\rm max}\approx 5.22\times 10^{-4}$ m s^-1 (at height $y\approx 0.362$ m) and $\approx 5.25\times 10^{-4}$ m s^-1 (at height $y\approx 0.350$ m).

We may then calculate the differences between the profiles of Test 2 and 3, with respect to Test 1, by utilizing the root-mean-square error (RMSE) as a metric. The approximate errors are $\approx 2.7\times 10^{-5}$ m s^-1 between the profiles 1 and 2 and $\approx 2.6\times 10^{-5}$ m s^-1 between the profiles 1 and 3. These results demonstrate that the experiments are reproducible with good accuracy.

From the numerical point of view, a particle-number independence test was necessary to determine the quality of the numerical solution as the number of SPH particles is increased. The experimental profile corresponding to Test 1 was selected as a validation standard. Figure 5 illustrates the quality of the numerical results as the number of particles is increased. The initial inter-particle separations at four different resolutions are 1.55 mm for $N=104976$ (dotted line), 1 mm for $N=251001$ (dot-dashed line), 0.7 mm for $N=511225$ (dashed line), and 0.5 mm for $N=10002001$ (solid line). The numerical profiles show an asympototic tendency to globally converge as the number of particles is increased. In terms of the RMSE metric the errors between the numerical and the experimental data for Test 1 are approximately $5.63\times 10^{-5}$ m s^-1 (for $N=104976$), $3.25\times 10^{-5}$ m s^-1 (for $N=251001$), $3.01\times 10^{-5}$ m s^-1 (for $N=511225$), and $2.83\times 10^{-5}$ m s^-1 (for $N=10002001$). On the other hand, the maximum velocities and heights are $v_{\rm max}\approx 5.10\times 10^{-4}$ m s^-1, $y\approx 0.350$ m (for $N=104976$), $v_{\rm max}\approx 5.15\times 10^{-4}$ m s^-1, $y\approx 0.355$ m (for $N=251001$), $v_{\rm max}\approx 5.20\times 10^{-4}$ m s^-1, $y\approx 0.365$ m (for $N=251001$), and $v_{\rm max}\approx 5.21\times 10^{-4}$ m s^-1, $y\approx 0.365$ m (for $N=10002001$). The curves for $N=511225$ and $N=10002001$ are seen to overlap with a RMSE deviation of $5.05\times 10^{-6}$ m s^-1. On the basis of these results, a spatial resolution resulting from $N=511225$ particles was adopted for all subsequent simulations conducted in the present study. At this spatial resolution, the numerical errors between the numerical and experimental profiles are observed to fall within the experimental uncertainty (see Fig. 4). The robustness of the method is further confirmed by the fact that these errors decay faster than the characteristic timescale of the heat diffusion process.

The SPH scheme is further validated against the lid-driven cavity test. This test has been used as a benchmark for convection heat transfer problems (see Hopp-Hirschler et al. [29] and references therein). The test consists of a circular, isothermal flow which sets in within a square cavity as a result of a constant velocity on the top wall of the cavity. A schematic of the cavity is shown in Fig. 6. We use the same parameters as in Hopp-Hirschler et al. [29], i.e., the side length of the cavity is $L=1$ mm, while the fluid density and dynamic viscosity are $\rho=1000$ kg m^-3 and $\eta=0.01$ Pa$\cdot$s, respectively. Three different calculations are considered for varying velocity of the top wall of the cavity, i.e., $v_{w}=0.1$ m s^-1, 1 m s^-1, and 10 m s^-1, corresponding to Reynolds numbers of ${\rm Re}=\rho Lv_{w}/\eta=10$, 100, and 1000, respectively. No-slip boundary conditions are applied at the walls of the cavity, while Neumann boundary conditions are employed for the pressure at the walls. At the bottom corners of the cavity the pressure gradients are either very low or even zero. Therefore, in order to avoid the onset of spurious pressures at the bottom corners, the reference pressure there is set to zero [29].

Figure 7 shows the dimensionless $y$-component of the velocity, $v_{y}/v_{w}$, as a function of the dimensionless horizontal coordinate position, $x/L$, at $y/L=0.5$ (left column of plots) and the dimensionless vertical coordinate position, $y/L$, as a function of the dimensionless $x$-component of the velocity, $v_{x}/L$, at $x/L=0.5$ (right column of plots) for ${\rm Re}=10$, 100, and 1000. The SPH profiles at a spatial resolution of $N=240\times 240$ particles are compared with the OpenFOAM simulations of Hopp-Hirschler et al. [29] using a mesh of $240\times 240$ regular cells. The results from these simulations are taken as reference solutions for comparison with the SPH calculations. The velocity profiles are very well reproduced by the present SPH scheme. Using as an indicator of convergence the position of the center of the vortex, the relative errors of the velocity profiles at this position between the $N=240\times 240$ SPH calculations and the reference OpenFOAM solutions are ($\epsilon_{v_{x}}=0.57$%, $\epsilon_{v_{y}}=0.49$%) for ${\rm Re}=10$, ($\epsilon_{v_{x}}=0.25$%, $\epsilon_{v_{y}}=0.30$%) for ${\rm Re}=100$, and ($\epsilon_{v_{x}}=0.59$%, $\epsilon_{v_{y}}=0.34$%) for ${\rm Re}=1000$. For comparison, Hopp-Hirschler et al. reported for this test relative errors of ($\epsilon_{v_{x}}=0.60$%, $\epsilon_{v_{y}}=0.42$%) for ${\rm Re}=10$, ($\epsilon_{v_{x}}=0.24$%, $\epsilon_{v_{y}}=0.31$%) for ${\rm Re}=100$, and ($\epsilon_{v_{x}}=0.60$%, $\epsilon_{v_{y}}=0.32$%) for ${\rm Re}=1000$ between their $N=240\times 240$ SPH results and the OpenFOAM solution.

Flow visualization was carried out by recording the displacement of the tracers within the cell. The first and third rows of Fig. 8 display the experimental velocity maps (first row) and their corresponding streamlines (third row) for $\Delta T=10$ K (${\rm Ra}=242$), $\Delta T=5$ K (${\rm Ra}=60.8$), and $\Delta T=1$ K (${\rm Ra}=15.2$). For comparison, the second and fourth rows show the corresponding fields produced by the SPH simulations for the same temperature differences. For all three temperature differences, an ascending stationary flow forms, which consists of a pair of counter-rotating vortices filling the upper part of the cavity just above the heated cylinder. The size and shape of the vortices depend on the temperature difference and Rayleigh number. From the last two rows of Fig. 8 it is evident that the vortices grow in size with the temperature difference. The red regions in the velocity maps correspond to zones of maximum velocity. In these regions the flow is vertically ascending (i.e., it moves away from the heat source) and the flow streamlines (experiments) and velocity vectors (SPH calculations) are closer together. As was previously found by Cesini et al. [5] using air at atmospheric pressure as the working fluid and Atmane et al. [11] using water, the flow pattern on the upper part of the cavity is characteristic of naturally convected heat transfer from a heated circular cylinder. The circulations expand in the upper part of the cavity, and as the temperature difference increases, both the buoyancy driven flow pattern and convection heat transfer increase. While the convective flow is confined in the upper cavity, heat transfer by conduction prevails at the bottom. Apart from minor details, the SPH calculations reproduce other relevant features of the experimentally obtained flow patterns. For instance, with an increase of the Rayleigh number, the density variation becomes greater, causing a more elongated thermal plume (the yellow and red regions in both the experimental and SPH velocity maps).

Figure 9 depicts the profiles of the vertical velocity component, $v$, at different heights in the upper part of the cavity for $\Delta T=10$ K (left) and $\Delta T=5$ K (right); the symbols correspond to the experimental data and the solid lines to the SPH results. At all heights, the vertical velocity component peaks just above the cylindrical heat source. The maximum velocities occur close to the source. At a height $y/H=0.68$ the velocities are $v\approx 5.62\times 10^{-4}$ m s^-1 (experiment) and $v\approx 5.54\times 10^{-4}$ m s^-1 (SPH) for $\Delta T=10$ K, while at a height $y/H=0.52$ the velocities are $v\approx 4.89\times 10^{-4}$ m s^-1 (experiment) and $v\approx 4.81\times 10^{-4}$ m s^-1 (SPH) for $\Delta T=5$ K. At greater heights, the maximum vertical velocity component decreases. Moreover, at either side of the cylinder the vertical velocity component decreases steeply across the vortices, reaching values as low as $-0.5\times 10^{-4}$ m s^-1, regardless of the height. Because of the circulations, the negative values of the vertical velocity components indicate an ascending flow near the lateral walls of the cavity. Again, the agreement between the experimental and SPH velocity profiles of Fig. 9 is measured in terms of the RMSE metric. These values are listed in Table 1. The lowest errors in both cases always occur when $y/H=0.52$, i.e., close to the heat source. At greater heights, the RMSE errors increase marginally. However, in all cases the errors are always less than $5\times 10^{-5}$ m s^-1.

Figure 10 shows the profiles of the vertical component of the velocity along a vertical line in the $z=0$ plane passing through the center of the source for temperature differences of 1, 5, and 10 K. These profiles correspond to the velocity field inside the red plumes shown in Fig. 8. In all cases, the magnitude of the velocity achieves maximum values in the inner region of the plume close to the source ($0.6<y/H<0.7$). As the temperature difference increases, the region of maximum velocity becomes progressively broader. This observation is consistent with the red plumes of Fig. 8 that become more elongated for higher temperature differences. A non-zero velocity exists at the heat source ($y/H=0.5$). This feature, as well as the form of the profiles shown in Fig. 10, are in good agreement with the predictions of Lauriat and Desrayaud [3], who reported 2D calculations for the time-dependent laminar buoyancy-induced flow above a horizontal line heat source immersed in an air-filled vessel. The SPH profiles match the experimental data with RMSEs of $\approx 1.4\times 10^{-4}$ m s^-1 for $\Delta T=10$ and 5 K, and of $\approx 8.1\times 10^{-5}$ m s^-1 for $\Delta T=1$ K. Further calculations with lower spatial resolution (not reported here) show that the agreement deteriorates particularly around the maximum. Therefore, we expect that by increasing the resolution above $N=511225$ particles the quality of the results will improve. From Figs. 8 and 9 it follows that the steady-state flow patterns are not exactly symmetric with respect to the central vertical axis ($x=0$). Asymptotic flow predictions by Angeli et al. [12] indicate that steady non-symmetric temperature and velocity fields may sometimes appear in double-cell circulation patterns for different ratios of the cylinder’s diameter to the cavity’s side and low ${\rm Ra}$-values.

When $\Delta T<0$, the cylindrical source is colder than the surrounding glycerin, thereby producing a vertically descending flow. A pair of counter-rotating vortices forms again, which expands in the bottom part of the cavity. Figure 11 shows the details of the flow for $\Delta T=-5$ K (${\rm Ra}=93.8$) and $\Delta T=-1$ K (${\rm Ra}=15.7$). The counter-rotating vortices are nearly symmetrical and advect glycerin from the bulk of the fluid located beneath the cylinder. A vertically descending flow with a maximum velocity of $-4.60\times 10^{-4}$ m s^-1 (experiment) and $-4.61\times 10^{-4}$ m s^-1 (SPH) takes place in the red region within the plume. The fine details of the flow are reproduced by the SPH simulations, because the velocity vectors closely follow the experimentally obtained streamlines. When the cooler glycerin reaches the bottom surface of the cavity, it splits up into two nearly symmetric cells occupying the entire lower half of the enclosure. This can be clearly observed in the experimental and SPH velocity maps of Fig. 11.

The profiles of the vertical velocity component at different heights from the bottom surface of the cavity are depicted in Fig. 12 for the descending flows shown in Fig. 11. The SPH solution (solid lines) is compared with the experimental data (symbols). Notice how the descending flow velocity always reaches a maximum just below the source. At either side of the source, the velocity decays sharply and reverts to positive values close to the lateral surfaces of the cavity. In contrast to the cases where $\Delta T>0$, the strength of the flow is more sensitive to the temperature difference. For instance, when $\Delta T=-1$ K, the inverted peaks become much less pronounced and the ascending flow near the lateral walls of the cavity is practically non-existent. Table 1 lists the RMSEs of the SPH solutions along with the experimental data. In all cases, the errors are always less than $6.28\times 10^{-5}$ m s^-1.

Another comparison between the SPH and the measured maximum values of the vertical velocity components is illustrated in Fig. 13, for all positive and negative temperature differences. Positive values of the velocity correspond to the ascending jet produced by a positive temperature difference, while negative values of the velocity are associated with negative temperature differences. The numerically predicted maximum velocities are in reasonably good agreement with the experimentally obtained values. The actual difference between both curves in terms of the RMSE is $8.41\times 10^{-6}$ m s^-1.