3.1. Microfluidic chip

The microfluidic chip, schematically depicted in figure 2a is designed to investigate the effect of different temperature gradients on fluid motion. Four Peltier modules are employed to establish the temperature gradient. Each module is connected to a power generation unit and regulated via an Arduino-UNO microcontroller. The microcontroller is interfaced with five negative-temperature-coefficient (NTC) thermistors, four of which are put in contact with each Peltier module to measure their individual temperatures, the fifth being deployed to monitor the ambient temperature. The microcontroller acts on four relay switches that control the Peltier modules by a sequence of ON/OFF cycles to maintain the desired set-point temperatures ( $T_{hot}$ and $T_{cold}$ ), with this approach, the temperature is maintained around the set-point values with a deviation of $\pm$ 0.5 K.

Figure 2.

(a) Exploded view of the microfluidic device used in the experimental campaign; (b) schematic representation of the optical systems used; (c) experimental apparatus used in this work.

The microfluidic chip is composed of a capillary glass tube (VitroTubes, VitroCOM) with a square cross-section of $800\times 800$ $\unicode {x03BC}$ m $^2$ and a length of 50 mm (1). The external sidewalls of the capillary are put in contact with four ad hoc 3D-printed aluminium cases (2) used to facilitate the heat transfer between the Peltier modules (3) and the capillary walls. All the components were assembled together in a 3D-printed polylactic acid case (4). A syringe pump (Harvard PHD 4400) was used to impose the streamwise, pressure-driven Poiseuille flow in the microchannel. The string pump is driven by a stepper motor that can generate small periodic fluctuations in the streamwise flow velocity at low flow rates. The impact of these fluctuations in our analysis will be discussed in detail in the results section. The whole experimental apparatus was mounted on an inverted microscope (Nikon Eclipse TE2000-U) during operation to allow the flow velocity measurements and it is shown in figure 2c.

Figure 3.

(a) Single measurement obtained with GDPT corresponding to a $1400 \times 800 \times 100\mu m^3$ measurement volume. Tracer particle trajectories were collected through 200 image recordings taken at 8.16 Hz. (b) Six measurement volumes were considered across the microchannel height with a spacing of 100 $\mu$ m. The $y$ -axis is plotted with direction from right to left for consistency with the 3-D plots.

3.2. General defocus particle tracking

The 3-D flow field in the microchannel was measured using GDPT (Barnkob et al. Reference Barnkob, Kähler and Rossi2015; Barnkob & Rossi Reference Barnkob and Rossi2020), a well-established single-camera method in which the depth position of the particles (i.e. parallel to the optical axis of the microscope objective) is determined by the defocusing pattern of the corresponding particle images.

The optical and image acquisition system is composed of the epi-fluorescent inverted microscope coupled to a high-sensitivity CCD camera (PCO SensiCam). A sketch of the optical system is presented in figure 2b. A high-power LED (Solis 525C, Thorlabs) emits green light with a wavelength of 525 nm. Thanks to a dichroic mirror, the light is directed to the microfluidic device through a 10 $\times$ objective lens with a numerical aperture NA = 0.25. Here, the fluorescent tracer microparticles absorb the green light and re-emit it at 607 nm. This emitted light passes through the dichroic mirror and reaches the CCD sensor of the camera. An additional cylindrical lens placed in front of the sensor induces astigmatic aberration in the defocused particle images, providing more efficient and unambiguous encoding of defocusing information for depth position determination (Barnkob et al. Reference Barnkob, Cierpka, Chen, Sachs, Mäder and Rossi2021). Polystyrene microparticles with a diameter of 5.0 $\unicode {x03BC}$ m and density of 1050 kg/m $^3$ (Microparticles GmbH) were used as tracer particles in the flow at low concentration. The images were processed using the open-source software DefocusTracker (Barnkob & Rossi Reference Barnkob and Rossi2021).

With this configuration, we achieved a measurement volume of approximately $1400 \times 800 \times 100$ $\unicode {x03BC}$ m $^3$ (length $\times$ width $\times$ height). To cover the entire volume of the channel, six measurement volumes across the microchannel channel height, spaced with a distance of 100 $\unicode {x03BC}$ m, were measured for each experimental case, as reported in figure 3. For each case, 200 images were recorded at a frame rate of 8.16 Hz. The estimated uncertainty in the determination of the position of individual particles was approximately 0.2 $\unicode {x03BC}$ m in the in-plane direction and 3 $\unicode {x03BC}$ m in the out-of-plane direction (Barnkob & Rossi Reference Barnkob and Rossi2020). The particle trajectories were calculated using a nearest-neighbour scheme, and a smoothing step based on the locally weighted scatter plot smooth (loess) method across 5 consecutive data points was applied to reduce the random velocity errors. A further correction procedure based on a reference measurement of a Poiseuille flow was applied to reduce the bias error resulting from the field curvature aberration and refractive index mismatch (Coutinho et al. Reference Coutinho, Moita, Ribeiro, Moreira and Rossi2023). The final uncertainty in the velocity measurement estimated on straight trajectories was approximately 2 and 8 $\unicode {x03BC}$ m/s for the in-plane and out-of-plane directions, respectively.

3.3. Experimental results

The experimental analysis focuses on the fluid flow motion under different temperature gradients and on the case of pure natural convection (without pressure-driven flow) and mixed convection (pressure-driven flow with a constant flow rate of Q = 0.2 ml/h). We perform in total 12 different experiments, considering 6 temperature differences between the two walls of the microchannel, corresponding to $\Delta \,T$ = 5, 10, 15, 20, 25, 30 K. A representative result of the GDPT measurement for $\Delta T$ = 30 K between the sidewalls is presented in figure 4, with and without the super-imposed pressure-driven flow.

Figure 4.

(a) Tracer particle trajectories measured in a swirling flow without the pressure-driven flow (Q = 0 ml/h) and a temperature difference $\Delta T$ = 30 K. (b) Corresponding velocity profile across the depth direction of the horizontal velocity component, $v_y$ . (c) Tracer particle trajectories measured in a swirling flow with the pressure-driven flow (Q = 0.2 ml/h) and a temperature difference $\Delta T$ = 30 K. (d) Corresponding velocity profile across the depth direction of the horizontal velocity component, $v_y$ . The missing data correspond to regions outside the measurement volumes.

The cold and hot walls are drawn in blue and red, respectively. The colour coding of the trajectories represents the value of the transverse velocity $v_y$ . In the case of natural convection (Figure 4a), the trajectories are confined in a plane approximately perpendicular to the channel axis. The projection of the trajectories in the $yz$ -plane shows clearly the convection vortex, with a magnitude of approximately $|v_y|_{max}$ $\simeq$ 90 $\unicode {x03BC}$ m/s. In the case with mixed convection (Figure 4b), the flow in the streamwise direction corresponds to a parabolic Poiseuille flow, with a maximum streamwise velocity of $|v_x|_{max}$ $\simeq$ 180 $\unicode {x03BC}$ m/s, and a buoyancy-driven vortex in the transverse direction, with a maximum magnitude of $|v_y|_{max}$ $\simeq$ 65 $\unicode {x03BC}$ m/s. The corresponding trajectories are, therefore, arranged in a helical motion. A point to consider is that the estimated sedimentation velocity of the 5- $\unicode {x03BC}$ m-diameter polystyrene tracer particles in water, across the temperature ranges observed in the experiments, varies from 0.52 to 1.37 $\unicode {x03BC}$ m/s. These values are small compared with the transverse velocities, but not negligible, especially for the cases with $\Delta T$ = 5 K. To avoid possible bias errors in the evaluation of the magnitude of the transverse velocity, we chose to estimate it in terms of the $y$ -component, $|v_y|_{max}$ , which is perpendicular to the gravitational acceleration and not directly affected by sedimentation.

Several interesting phenomena emerge from the experimental analysis. First, as shown in figure 5, the actual streamwise velocity imposed by the syringe pump oscillates significantly around its nominal value (corresponding to Q = 0.2 ml/h) due to the working principle and limitation of the syringe pump used. These oscillations can reach up to 50 $\unicode {x03BC}$ m/s above the mean value. On the other hand, the transverse velocity, driven by the buoyancy forces and by the thermal gradient steadily controlled by the Peltier elements, remains stable with no appreciable fluctuations, independently of the applied.

Figure 5.

Maximum value of the $y$ and $x$ components of velocity obtained for three different temperature differences ( $\Delta T$ = 5, 15 and 25 K ), plotted as a function of time.

Figure 6.

Top view of the particles trajectories in different cases: (a) $\Delta T$ = 30 K and Q = 0.2 ml/h; (b) $\Delta T$ = 30 K and Q = 0.0 ml/h; (c) $\Delta T$ = 20 K and Q = 0.0 ml/h.

A second observation is that the plane of the buoyancy-driven transverse vortex in the natural convection case is not perfectly perpendicular to the channel streamwise axis, as could be expected from symmetry reasons. This could be explained by a slight tilt of the channel axis with respect to the gravity vector, i.e. the microfluidic device or the microscope stage was not perfectly horizontal. This explanation is supported by the fact that the tilting angle stays approximately constant also for different temperature differences and vortex magnitude, as shown in figure 6b–c from the top view of the trajectories presented for $\Delta T$ = 30 K and $\Delta T$ = 20 K.

Another aspect to consider concerns the magnitude of the transverse vortices which is, as expected, proportional to the imposed $\Delta T$ . However, for the cases without the pressure-driven flow, given the same magnitude of the temperature gradient, the magnitude of the transverse flow is larger with respect to the corresponding cases with mixed convection. This was consistently observed in all experiments, as shown in figure 7. In this plot, the differences in terms of absolute transverse velocity between the two configurations become clear. In all the cases with Q = 0 ml/h (blue line), $|v_y|_{max}$ is always larger than the corresponding case with Q = 0.2 ml/h (orange line). Interestingly, the difference becomes approximately constant for $\Delta T \gt 10$ K. This effect can be related to the improved heat transfer provided by the pressure-driven flow that reduces the effective temperature difference between the walls. Finally, the findings observed in figure 7 for the mixed convection case are summarised in table 1 in terms of the corresponding Reynolds number, Re = $Q / (D_h \nu )$ , Grashof number, Gr = $g \beta \Delta T D_h^3/ \nu ^2$ , Richardson number, Ri = Gr/Re $^2$ and the non-dimensional transverse velocity $v^* = |v_y|_{\textrm {max}}/|v_x|_{\textrm {max}}$ (Hardt et al. Reference Hardt, Pennemann and Schönfeld2006). The parameters $\nu$ and $\beta$ are the kinematic viscosity and the coefficient of volume expansion of the fluid, respectively. It is interesting to notice that, although the large Ri indicates fully natural convection, the ratio between transverse and streamwise velocity shows that forced convection is of the same order as or larger order of magnitude than the natural convection, indicating that a mixed convection regime is established in the system.

Table 1.

Values of Re, Gr, Ri and velocity ratio $v^*$ for the investigated experimental condition. The flow rate is constant for all cases at Q = 0.2 ml/h

Figure 7.

Comparison of the magnitude of the transverse buoyancy-driven vortex, expressed as $|v_y|_{max}$ , as a function of the $\Delta T$ , for the case without and with pressure-driven flow (Q = 0 and 0.2 ml/h).

4.1. Numerical methodology

The simulations were conducted with water as the working fluid, considered as a non-compressible and Newtonian fluid. The flow was assumed to be in a laminar regime under steady-state conditions. The computational analyses were performed utilising STAR-CCM+, a commercially available computational fluid dynamics software. This software employs the finite-volume method to solve the Navier–Stokes–Fourier equations. The set of governing equations employed in this study is reported below

(1) \begin{equation} \nabla \cdot (\rho \mathbf{u}) = 0, \end{equation}

(2) \begin{equation} \nabla \cdot (\rho \mathbf{u} \mathbf{u}) = - \nabla p + \mu \nabla ^2\mathbf{u} + \rho \mathbf{g}, \end{equation}

(3) \begin{equation} \nabla \cdot (\rho c_p T \mathbf{u}) = -\nabla \cdot (k \nabla T) ,\end{equation}

(4) \begin{equation} \nabla \cdot (\rho c \mathbf{u}) = -\nabla \cdot (D \nabla c). \end{equation}

Here, $\mathbf{u}$ is the flow velocity field, $\rho$ is the fluid density, $p$ is the pressure, $T$ is the temperature, $\mu$ is the dynamic viscosity, $k$ is the thermal conductivity coefficient and $c_p$ is the specific thermal capacitance at constant pressure. The last equation represents the convection–diffusion equation, where $c$ is a non-dimensional concentration and the diffusion coefficient $D$ is taken equal to 2.8 $\cdot 10^{-10}$ m $^2$ /s, which corresponds to the diffusivity of a large dye molecule like Rhodamine B in water. The associated Schmidt number $Sc$ is equal to 3180. The concentration $c$ is normalised with respect to the maximum concentration of the Rhodamine, taking the value between 0 (pure water) and 1 (maximum/initial concentration in one half of the channel).

The computational domain was segmented into three distinct regions: the inlet region with a length $L_{in}$ = 15 $D_h$ where the thermal gradient is not applied, a central or active zone with a length $L$ = 50 mm where the temperature gradient is applied and an outlet region with a length $L_{out}$ = 5 mm used to avoid reversed flow in the central domain (see figure 1). The length of the inlet region was chosen to be dependent on the hydraulic diameter to facilitate comparisons between different geometries of the channel since the diffusive mixing is also active in the region without thermal mixing. With this choice, for different geometries and a fixed Reynolds number, the levels of mixing achieved at the entrance of the active zone will be the same.

The origin of the reference frame is located at the beginning of the active region. The mesh utilised in this study consisted of a structured grid composed of uniform hexahedral elements. A sensitivity analysis based on the relative error calculated with respect to the finest grid (6.4 million cells) of the $v_y$ value across a line at $z$ = 0 in the centre of the channel was performed to determine the optimal number of grid elements. The results, reported in figure 8, showed that 1.22 million elements, corresponding to 26 $\times$ 26 elements on the inlet plane, were sufficient to achieve mesh convergence. Regarding the flow boundary conditions, a fully developed Poiseuille velocity profile was prescribed at the inlet ( $x=-L_{in}$ ), while an outflow condition was applied at the outlet ( $x=L+L_{out}$ ). The fluid enters the domain at ambient temperature ( $T_{in}$ = 303.15 K). Concerning the thermal boundary conditions for the channel walls, both in the inlet and outlet passive zones, the walls are considered adiabatic. In the central region, Dirichlet boundary conditions were applied at the walls to produce the necessary temperature gradient. The wall at $y=+D_h/2$ was set to $T_{hot}$ and the wall at $y=-D_h/2$ to $T_{cold}$ . Throughout all simulations, $T_{cold}$ and $T_{hot}$ were kept constant at 293.15 K and 353.15 K, respectively, and with a temperature gradient between the two walls equal to $\Delta T/D_h$ = 50 K/mm. The top ( $z=+D_h/2$ ) and bottom ( $z=-D_h/2$ ) channel walls were treated as adiabatic surfaces. It is crucial to note that the gravity vector points in the negative z-direction, defined as g = (0,0,−9.81) m/s $^2$ . Additionally, all the thermophysical and transport properties of the fluid are modelled as functions of temperature following the International Association for the Properties of Water and Steam (IAPWS) model (Wagner & Kretzschmar Reference Wagner and Kretzschmar2008).

Figure 8.

Sensitivity analysis of the mesh on the $v_y$ . (a) Relative velocity error as a function of the number of elements. The selected mesh with 1.22 million elements (highlighted with a larger marker) exhibits a similar error to the finest grid. (b) Error distribution along the considered line. Also in this case, the grid with 1.22 million elements follows a similar trend to the finest meshes.

4.2. Numerical results

First, the numerical results were compared with the experimental measurements for validation, showing a good agreement in the distribution and magnitude of the transverse vortex. The case with Re = 1 and $\Delta T$ = 20 K is presented in figure 9.

Figure 9.

Comparison between the v $_{y}$ field from the experimental and numerical analyses for the case with Re = 1 and $\Delta T$ = 20 K.

Figure 10.

Concentration maps obtained from numerical simulations for different values of Reynolds number (0.1,1 and 10).

For the mixing analysis, we considered the distribution of the non-dimensional concentration $c$ . At the inlet, we arbitrarily chose to put the mixture of water and Rhodamine at maximum concentration ( $c = 1$ ) in the half-channel region close to the hot wall ( $y\lt 0$ ) and the pure water ( $c=0$ ) in the other region ( $y\gt 0$ ). Figure 10 shows some exemplary concentration fields, taken in a plane in the centre of the channel ( $z/D_h = 0$ ) for different Re values and $D_h$ = 800 $\unicode {x03BC}$ m. The figure provides a qualitative comparison between the swirl-on and swirl-off cases. The swirl-off cases were obtained at the same conditions (temperature, flow rate) but removing the gravitational field. As expected, in the swirl-off cases the mixing mechanism is slow and relies only on diffusion. When the swirl is activated (swirl on), the contact surface between the two fluids increases, thus improving mixing. Since the swirl magnitude is fixed by the , the active mixing (swirl on) is more effective at low Re. At larger Re, where the magnitude of the transverse flow is small compared with the streamwise velocity, the difference between the swirl-on and swirl-off cases is negligible.

Figure 11.

Mixing efficiency $\phi$ against the dimensionless position $x/D_h$ for different Reynolds numbers and different hydraulic diameters ( $D_h$ ); the black dashed line refers to the beginning of the active region, always located at $x/D_h = 0$ .

The mass transport enhancement can be characterised by introducing a mixing efficiency parameter ( $\phi$ ), defined as

(5) \begin{equation} \phi = 1- \frac {\sum _{i=1}^{n_{el}} | c_i-\bar {c}|}{2|\bar {c}| n_{el}} ,\end{equation}

where $c_i$ is the value of the non-dimensional concentration of the tracking fluid evaluated in the ith element of a certain cross-section, $\bar {c}$ is the average value of the concentration in the considered cross-section (here always equal to 0.5 for the mass conservation) and $n_{el}$ is the number of elements of the cross-section. A value of $\phi = 1$ corresponds to a complete mixing. This parameter and the pressure drop $\Delta p$ were analysed for different geometries ( $D_h$ ) and Reynolds numbers (Re). Each result was compared with the respective reference swirl-off case, indicated with the subscript 0.

In figure 11 we present the evolution of the mixing efficiency parameter $\phi$ in Eq. (5), evaluated across the capillary length for different hydraulic diameters and Re values. The streamwise direction was normalised with the hydraulic diameter of the channel ( $x/D_h$ ), with the position $x/D_h=0$ corresponding to the entrance of the active region (where the temperature gradient is activated). The results are presented from $x/D_h = -15$ , which is the beginning of the inlet region, up to $x/D_h = 25$ . Using this scaling, it is possible to compare the performance of different microchannels with different hydraulic diameters. The black dashed line, shown in figure 11, refers to the starting point of the active region. For the smallest diameter $D_h=100$ $\unicode {x03BC}$ m, the swirling effect is almost negligible since viscous forces dominate over buoyancy forces; the mixing is due to pure molecular diffusion and the concentration evolution is the same for the swirl-on and swirl-off cases. For all the other cases, the effect of the swirling is not negligible, and one can appreciate a clear enhancement of the mixing of the fluid for all the Reynolds numbers after the fluid enters the active region. The higher mixing efficiencies at the exit of the channel are observed, as expected, for the lower Reynolds numbers due to the larger residence time of the fluid in the channel. Interestingly, for the cases with larger $D_h$ , the mixing enhancement starts before the active region (before the dashed red line), up to 5 times the hydraulic diameter of the capillary. This is linked to the backward propagation of the swirling motion from the active region, which becomes more effective as the $D_h$ is increased.

Figure 12.

(a) Value of $\phi$ at $x/D_h = 25$ as a function of the hydraulic diameter for the different Reynolds numbers. (b) Mixing efficiency parameter $PEC_{\phi }$ as a function of the hydraulic diameter for the different Reynolds numbers.

A meaningful insight into the overall mixing efficiency is obtained by comparing the swirl-on and swirl-off mixing efficiency ratio ( $\phi /\phi _0$ ) at the reference distance $x/D_h = 25$ , as a function of $D_h$ and parametrised with the Reynolds number. As illustrated in figure 12a, when the swirling is activated the mixing is always enhanced. This is true for all the cases analysed except for the smaller channel case ( $D_H=100 \ \unicode {x03BC}$ m), where the convective motion in the flow is not activated.

Even if the general trend of mass transport enhancement is clear, the mixing efficiency ratio behaviour is nonlinear with respect to Re. For the cases with a low Reynolds number (Re = 0.1), the value of $\phi /\phi _0$ as a function of $D_h$ reaches an asymptotic value for $D_h \ge 800$ $\unicode {x03BC}$ m, showing that the complete mixing has been achieved already before $x/D_h = 25$ . On the other hand, the asymptotic value of $\phi /\phi _0$ is equal to 1.1, showing that the effective mixing enhancement introduced by mixed convection is limited for low Re. At larger Reynolds numbers (Re = 10), the diffusion and the swirl magnitude are small compared with the streamwise advection, and $\phi /\phi _0$ remains low as well. To achieve good mixing, a longer channel length is required, and the mixing efficiency starts to increase effectively only for the larger hydraulic diameters, with negligible enhancement observed for $D_h \lt 1200 \ \unicode {x03BC}$ m. The most significant improvement was observed for the case of Re = 1, where, for $D_h=2000 \ \unicode {x03BC}$ m, the swirling motion has more than 60 % mixing enhancement.

Finally, to have a clearer picture with respect to the advantages of using the thermally induced swirling motion for mass transport enhancement, it is of interest to check the mechanical energy needed to drive the pressure-driven flow (i.e. the pressure drop $\Delta p$ ), which increases as the swirl magnitude is increased. A comparison was performed by analysing the performance efficiency parameter given by

(6) \begin{equation} PEC_\phi = \frac {\phi }{\phi _0}\big / \Big (\frac {\Delta p}{\Delta p_0}\Big )^{1/3}, \end{equation}

where the subscript 0 refers to the swirl-off case (Figure 12b). A value of $PEC_\phi \gt 1$ means that the mass transfer enhancement is predominant with respect to the pressure-drop increase, corresponding to a net overall enhancement of the efficiency of the system. One can readily see that, for the cases where swirling is present and large compared with streamwise advection (low Re), one needs to exert higher pumping power with respect to swirl-off flows, resulting in a low $PEC_\phi$ . Also for the case of large Re, the $PEC_\phi$ is larger than one only for larger channels with $D_h \gt 1600$ $\unicode {x03BC}$ m, therefore moving outside the microfluidic application.

In conclusion, from a fluid dynamic perspective, we can approximate the swirling flow to a superposition of a buoyancy-driven transverse vortex and a pressure-driven flow, that can be regulated independently by setting the temperature gradient and the flow rate. The transverse and streamwise flows remain essentially decoupled, although a larger buoyancy-driven vortex requires a larger pressure drop to sustain the streamwise flow. The mixing efficiency with respect to a molecular solute is promoted when the transverse velocities are large with respect to molecular diffusion and not too small with respect to the streamwise velocities. If these conditions are not verified, the mixing enhancement obtained with this approach is limited. For the fluid and temperature differences considered in this work, an optimal range of operation emerges for Re = 1 and $D_h$ between 400 and 1600 $\unicode {x03BC}$ m. This range is particularly interesting from a technological point of view, since it covers microfluidic devices working at low flow rates, where diffusion is still too slow and active mixing is required.