Experimental and Numerical Studies of the Temperature Field in a Dielectrophoretic Cell Separation Device Subject to Joule Heating

1. Introduction

In recent years, technologies for the rapid and accurate separation of rare cells from a large number of diverse cell populations have attracted much attention. In particular, cell separation technologies based on dielectrophoresis (DEP) have attracted considerable interest because they are noninvasive to cells [1,2,3]. DEP is a phenomenon in which an electromechanical force (DEP force) is generated on dielectric microparticles placed in a nonuniform electric field [4].

The application of the DEP principle utilizing an A.C. electric field makes it possible to selectively manipulate biological microparticles, such as cells, with high precision according to their inherent electrical properties. This capability contributes to the development of critical fundamental technologies in life sciences, such as cell separation, concentration, purification, and manipulation technologies. Therefore, there is a strong demand for the development of rapid and high-precision cell separation technologies that leverage the advantages of DEP [5,6,7].

Among the many forms of DEP-based microdevices, the electrode-based dielectrophoresis (eDEP) devices are distinct for the large electric field gradient generated within the microdevice by microelectrodes [8,9,10,11]. By applying an A.C. voltage to the eDEP device, a large electric field gradient is generated around the fine electrodes, and this electric field gradient generates a DEP force that acts on cells. Because the electric field gradient determines the magnitude and direction of the DEP force acting on cells, it is necessary to generate within the device an electric field of adequate strength to manipulate target cells with high precision. However, cell separation devices based on dielectrophoresis have the drawback that the cells inside the device are exposed to a high temperature environment due to Joule heating [12,13]. The degree of Joule heating mostly depends on the square of the electric field magnitude [14,15]. Thus, there is a risk that cells inside the device are exposed to a high-temperature environment due to the Joule heat generated by the application of high voltage [16,17].

It has been confirmed that the high-temperature exposure may affect the physiological functionality of the cells [18]. Moritz et al. conducted experiments using human and porcine epidermal cells, and reported that irreversible cell damage occurs when the cells are exposed to a high-temperature environment above 44 $℃$ for approximately 6 h [19]. Harmon et al. reported that apoptosis was significantly promoted by exposing Murine Mastocytoma cells to a high-temperature environment of around 44 $℃$ [20]. Yamada et al. reported that cell division was halted when MKN-1 and SNG-M cells were exposed to a high-temperature environment of 42 $℃$ [21]. Consequently, recent studies have focused on Joule heating and temperature fields within dielectrophoresis-based cell separation devices to evaluate their effects on cells [22,23,24]. In particular, the correlation between solution conductivity and the temperature increase has been extensively investigated. For example, Tay et al. measured the temperature field of planar electrodes using a NaCl solution with ${\sigma }_f=1$ S/m, and reported a temperature rise of up to 100 $℃$ [22]. Sridharan et al. conducted an observation experiment on fluid circulation within a device using a phosphate-buffered solution with ${\sigma }_f=4.7\times {10}^{-2}$ S/m and predicted a temperature change of about 30 $℃$ in their numerical simulations [23]. Nakano et al. measured the temperature field within an insulated DEP (iDEP) device using a phosphate-buffered solution with ${\sigma }_f=1.0\times {10}^{-1}$ S/m, and found that the temperature rose about 70 $℃$ [24]. Among other relevant studies that measured the temperature rise due to Joule heating in microdevices, Williams et al. measured the temperature using the laser-induced fluorescence (LIF) method with varying solution conductivity and applied voltage [25]. Pramod et al. measured the temperature distribution in the microchannel due to Joule heating using the micro-LIF method [26,27].

Additionally, comparisons between experimentally observed and numerically simulated results of temperature field analysis in flow fields have been reported. For example, in one of the studies focused on experimental measurement, Erickson et al. measured the temperature rise due to Joule heating inside a three-dimensional microfluidic device made of poly (dimethylsiloxane) (PDMS) and glass, and reported that the results matched the numerical simulations with an accuracy of $\pm 3$ $℃$ [28]. Similarly, Williams et al. compared experimental results obtained using the LIF method with numerical simulation results, and reported that the temperature rise was proportional to the square of the applied voltage [29].

With respect to the previous studies focused on numerical simulation, Tang et al. proposed a three-dimensional numerical simulation model of thermal fluid flow to investigate the correlation between Joule heating and temperature gradient focusing, and reported that the numerical simulation results matched well with the experimental results [30]. Rasmussen et al. conducted modeling to optimize flow sensors within microchannels and reported that the simulation results and the temperature field measurements using sensors based on the simulations matched with an error of less than 1 % [31]. Many analyses of microflows have been performed using commercial software. Gallo et al. numerically predicted the temperature field due to Joule heating in an insulator-based DEP (iDEP) device [32]. Kwak et al. performed numerical simulations focusing on localized dielectric loss-induced heating in DEP devices and reported that the origin of the heating in DEP devices with electrically insulated electrodes was Debye relaxation [33]. Wang et al. numerically predicted the interactions between electrothermal flow and electroosmotic flow near the constriction region in the microchannel and clarified the temperature rise and flow field near the constriction [34]. Khoshmanesh et al. analyzed the temperature field in the unique DEP device they developed [35].

In this study, with the aim of studying the effects of high-temperature exposure on cells, the temperature rise in the microchannel of a previously developed eDEP cell separation device was evaluated experimentally and numerically [36]. For the experimental evaluation, the micro-LIF method was adopted to measure the temperature in the device. In the numerical approach, the thermal structure of the device was investigated by performing a numerical simulation in which the effect of heat generation due to a nonuniform electric field was incorporated in the thermofluidic model. Furthermore, parametric numerical simulations were performed with changes in the solution conductivity and sample solution flow rate, which affect cell separation performance. Finally, the correlation between the device operating conditions and temperature rise, and their effects on cell viability were explored and discussed.

2. Thermofluidic Numerical Simulation

2.1. Thermofluidic Governing Equations

Figure 1 shows a conceptual diagram of the numerical simulation model. The bottom surface of the microfluidic channel with height $H=500$ µm consisted of a 400 µm thick glass substrate with counter-interdigitated electrodes printed on it, and the top surface consisted of a 1100 µm thick indium tin oxide (ITO) planar electrode. The length of the microfluidic channel was $L_z=3600$ µm. Cell separation is performed by generating a nonuniform electric field while flowing sample solution-suspended cells through the microfluidic channel at a constant volumetric flow rate $Q$.

In the case of analyzing the flow of the electrolyte solution inside the microfluidic channel with the applied A.C. electric field, it was also necessary to evaluate the effects of electroosmotic flow and electrothermal flow near the electrode surfaces on the temperature field, in addition to Joule heating. It is known that the effect of the electroosmotic flow is likely to appear when the electric conductivity of the solution is ${\sigma }_f\le 8.4\times {10}^{-2}$ S/m [37]. Ren et al. reported that the order of magnitude of the electroosmotic flow was less than $1.0\times {10}^{-6}$ m/s under the conditions of an applied voltage $V=30$ V_pp and conductivity ${\sigma }_f=1.0\times {10}^{-3}$ S/m [38]. Electrothermal flows have been reported to be more pronounced when high conductivity solutions (${\sigma }_f~0.2-1.0$ S/m) were used [37]. Akutsu et al. reported that the order of magnitude of electrothermal flow was about $7.5\times {10}^{-5}$ m/s under conditions of an applied voltage $V=20{\mathrm{V}}_{\mathrm{p}\mathrm{p}}$ and conductivity ${\sigma }_f=2.9\times {10}^{-1}$ S/m [39]. In this study, an applied voltage of $V=10{\mathrm{V}}_{\mathrm{p}\mathrm{p}}$ and conductivity of ${\sigma }_f=1.0\times {10}^{-3}~1.0\times {10}^{-1}$ S/m were adopted. Furthermore, taking into account that the average flow velocity of the sample solution was $U_m~0.5\times {10}^{-3}$ m/s, the effects of electroosmotic and electrothermal flows on the thermofluidic field in the microchannel could be considered negligible. Regarding the buoyancy effect due to the temperature distribution in the microfluidic channel, the temperature difference in the microfluidic channel, $∆T$, was only a few degrees Celsius, and the order of the magnitude of the Rayleigh number, $\mathrm{R}\mathrm{a}$, was $\mathrm{R}\mathrm{a}~100$. Therefore, the buoyancy effect could also be considered negligible. The spanwise direction of the microchannel was defined as the $x$-axis, the height direction as the $y$-axis, and the flow direction as the $z$-axis. The sample solution was assumed to be an incompressible homogeneous viscous fluid. A fully developed Poiseuille flow was assumed for the velocity distribution. Thus the velocity vector of the fluid was given as

$$\begin{array}{c}\mathit{v}=\left(0,0,{\displaystyle \frac{6Q}{H^3}}\left(H-y\right)y\right)\end{array}$$

(1)

Here, $H$ is the height of the microfluidic channel and $y$ is the height from the bottom surface ($0\le y\le H$). $Q$ represents the volumetric flow rate of the sample solution, and it was given using the width of the microfluidic channel, $B$, as $Q=U_{m}BH$. The governing equations for liquid and solid phases in the device are given as

$$\begin{array}{c}{\displaystyle \frac{𝜕T}{𝜕t}}+\left(\mathit{v}·\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\right)T={\displaystyle \frac{{\kappa }_f}{{\rho }_f{C}_{pf}}}\left({\displaystyle \frac{{𝜕}^{2}T}{𝜕x^2}}+{\displaystyle \frac{{𝜕}^{2}T}{𝜕y^2}}+{\displaystyle \frac{{𝜕}^{2}T}{𝜕z^2}}\right)+{\displaystyle \frac{q}{{\rho }_f{C}_{pf}}}\end{array}$$

(2)

$$\begin{array}{c}{\displaystyle \frac{𝜕T}{𝜕t}}={\displaystyle \frac{{\kappa }_g}{{\rho }_g{C}_{pg}}}\left({\displaystyle \frac{{𝜕}^{2}T}{𝜕x^2}}+{\displaystyle \frac{{𝜕}^{2}T}{𝜕y^2}}+{\displaystyle \frac{{𝜕}^{2}T}{𝜕z^2}}\right)\end{array}$$

(3)

Here, the values $t$, $\kappa$, $\rho$, and $C_p$ are time, thermal conductivity, density, and isobaric specific heat, respectively. Subscripts $f$ and $g$ represent fluid and solid, respectively. The second term of the right-hand side of Eq. (2) represents the heat generation. The Joule heating $q=q\left(x,y\right)$ in the electrolyte solution can be given as

$$\begin{array}{c}q={\sigma }_f{E}^2\end{array}$$

(4)

where ${\sigma }_f$ is the electric conductivity of the electrolyte solution and $E$ is the effective value of the electric field.

2.2. Thermofluidic Boundary Conditions

The boundary conditions of thermofluid analysis are shown in Figure 2(a). The constant temperature boundary condition $T=T_0$ was given for the upstream end of the device ($z=0$). The freely developed boundary condition of the temperature field was given at the downstream end of the device. The constant wall heat flux condition

$$\begin{array}{c}-{\kappa }_f{\left({\displaystyle \frac{𝜕T}{𝜕y}}\right)}_f=-{\kappa }_g{\left({\displaystyle \frac{𝜕T}{𝜕y}}\right)}_g\end{array}$$

(5)

was applied to the interface between the liquid and the solid phases in the microfluidic channel. At the top surface of the device, the boundary condition

$$\begin{array}{c}-{\kappa }_g{\left({\displaystyle \frac{𝜕T}{𝜕y}}\right)}_g=h_t\left(T_0-T_{wt}\right)\end{array}$$

(6)

was applied. Here, $h_t$ is the convective heat transfer coefficient at the top surface and $T_{wt}$ is the local temperature at the top surface. The $h_t$ was evaluated using the following empirical formula:

$$\begin{array}{c}{h}_t={\displaystyle \frac{{\kappa }_a}{L}}{C}_t{{\mathrm{R}\mathrm{a}}_t}^{\zeta }\end{array}$$

(7)

where $L$ represents the characteristic length of the device and was defined as follows using the perimeter length, $l$, and area, $A$, of the electrode substrate:

$$\begin{array}{c}L={\displaystyle \frac{4A}{l}}\end{array}$$

(8)

The coefficients $\zeta$ and $C_t$ were determined empirically. ${\mathrm{R}\mathrm{a}}_t$ is the Rayleigh number at the top surface defined as

$$\begin{array}{c}{\mathrm{R}\mathrm{a}}_t={\displaystyle \frac{g{\beta }_a\left(T_0-T_{wt}\right)L^3}{{{\nu }_a}^2}}{\mathrm{P}\mathrm{r}}_a\end{array}$$

(9)

where the values $g$, ${\beta }_a$, ${\nu }_a$ and ${\mathrm{P}\mathrm{r}}_a$ are the gravitational acceleration, the thermal expansion coefficient of the air, kinematic viscosity of the air, and Prandtl number of the air, respectively. At the bottom surface of the device, the boundary condition

$$\begin{array}{c}{h}_b\left(T_0-T_{wb}\right)=-{\kappa }_g{\left({\displaystyle \frac{𝜕T}{𝜕y}}\right)}_g\end{array}$$

(10)

was applied. Here, $h_b$ is the convective heat transfer coefficient at the bottom surface and $T_{wb}$ is the local temperature at the bottom surface. In a similar manner as for $h_t$ in Eq. (7), $h_b$ was evaluated using the following empirical formula:

$$\begin{array}{c}{h}_b={\displaystyle \frac{{\kappa }_a}{L}}{C}_b{{\mathrm{R}\mathrm{a}}_b}^{\zeta }\end{array}$$

(11)

Here, ${\mathrm{R}\mathrm{a}}_b$ is the Rayleigh number at the bottom surface defined as

$$\begin{array}{c}{\mathrm{R}\mathrm{a}}_b={\displaystyle \frac{g{\beta }_a\left(T_0-T_{wb}\right)L^3}{{{\nu }_a}^2}}{\mathrm{P}\mathrm{r}}_a\end{array}$$

(12)

The coefficients $\zeta$ and $C_b$ were determined empirically. For both the lateral faces of the computational domain, the symmetry boundary condition

$$\begin{array}{c}{\displaystyle \frac{𝜕T}{𝜕x}}=0\end{array}$$

(13)

was applied.

2.3. Governing Equation for Electric Field

The Joule heating per unit time, $q$, as expressed in Eq. (4), is a function of the effective value of the A.C. electric field, $E$. Therefore, in this study, the electric field in the microchannel was obtained by electrostatic field analysis with $E$ as the dependent variable. Initially, because fingers of the interdigitated electrodes were arranged along the flow direction ($z$-axis), the electric field in the microfluidic channel could be assumed to be uniform along the flow direction. In other words, the electric field only had a distribution within the cross-section of the microfluidic channel. Therefore, the electric field in the microchannel was obtained by solving the Laplace equation for the electrostatic potential, $\psi$, in the $xy$-plane of the microchannel

$$\begin{array}{c}{\displaystyle \frac{{𝜕}^2\psi }{𝜕x^2}}+{\displaystyle \frac{{𝜕}^2\psi }{𝜕y^2}}=0\end{array}$$

(14)

The electrostatic field $\mathit{E}$ ($\left|\mathit{E}\right|=E$) was obtained from the relation

$$\begin{array}{c}E=\left(E_x,E_y\right)=\left(-{\displaystyle \frac{𝜕\psi }{𝜕x}},-{\displaystyle \frac{𝜕\psi }{𝜕y}}\right)\end{array}$$

(15)

using $\psi$ obtained by numerically solving Eq. (14).

2.4. Boundary Conditions for the Electric Field

Figure 2 (b) shows the electric field analysis model. The computational domain was a part of the microfluidic channel cross-section, including half of a pair of interdigitated electrodes. The boundary conditions for $\psi$ were given as

$$\begin{array}{c}\left\{\begin{array}{c}\psi =V_{\mathrm{r}\mathrm{m}\mathrm{s}}\left(\mathrm{H}\mathrm{i}\mathrm{g}\mathrm{h}-\mathrm{v}\mathrm{o}\mathrm{l}\mathrm{t}\mathrm{a}\mathrm{g}\mathrm{e}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{d}\mathrm{e}\right)\\ \psi =0\left(\mathrm{G}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{d}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{s}\right)\end{array}\right.\end{array}$$

(16)

The symmetry boundary condition

$$\begin{array}{c}{\displaystyle \frac{𝜕\psi }{𝜕x}}=0\end{array}$$

(17)

was applied to both the lateral faces of the computational domain ($x=0,100$ µm). The boundary condition

$$\begin{array}{c}{\displaystyle \frac{𝜕\psi }{𝜕y}}=0\end{array}$$

(18)

was applied to the bottom surface of the microfluidic channel ($25\mathsf{\mu }\mathrm{m}\le x\le 75\mathsf{\mu }\mathrm{m},y=400$ µm). The boundary condition

$$\begin{array}{c}{E}_x=0\end{array}$$

(19)

was applied to the electrode surfaces, and both the lateral faces of the computational domain. The boundary condition

$$\begin{array}{c}{E}_y=0\end{array}$$

(20)

was applied to the bottom surface of the microfluidic channel ($25\mathsf{\mu }\mathrm{m}\le x\le 75\mathsf{\mu }\mathrm{m},y=400$ µm).

2.5. Numerical Simulation Schemes

The governing equations and boundary conditions were discretized using the first-order accurate finite difference scheme in time and the second-order accurate scheme in space. For example, the discretized Eq. (3) is denoted as

$$\begin{array}{c}{T}_{i,j,k}^{n+1}={\alpha }_0{T}_{i,j,k}^n+∆t\left({\alpha }_1{T}_{i+1,j,k}^{n+1}+{\alpha }_2{T}_{i-1,j,k}^{n+1}+{\alpha }_3{T}_{i,j-1,k}^{n+1}+{\alpha }_4{T}_{i,j+1,k}^{n+1}+{\alpha }_5{T}_{i,j,k-1}^{n+1}+{{\alpha }_{6}T}_{i,j,k+1}^{n+1}\right)\end{array}$$

(20)

$$\left\{\begin{array}{c}{\alpha }_0=1+∆t{\displaystyle \frac{2{\kappa }_g}{{\rho }_g{Cp}_g}}\left({\displaystyle \frac{1}{∆x_i∆x_{i-1}}}+{\displaystyle \frac{1}{∆y_j∆y_{j-1}}}+{\displaystyle \frac{1}{∆z_k∆z_{k-1}}}\right)\\ {\alpha }_1={\displaystyle \frac{2{\kappa }_g}{{\rho }_g{Cp}_g}}{\displaystyle \frac{1}{∆x_i\left(∆x_i+∆x_{i-1}\right)}}\\ {\alpha }_2={\displaystyle \frac{2{\kappa }_g}{{\rho }_g{Cp}_g}}{\displaystyle \frac{1}{∆x_{i-1}\left(∆x_i+∆x_{i-1}\right)}}\\ {\alpha }_3={\displaystyle \frac{2{\kappa }_g}{{\rho }_g{Cp}_g}}{\displaystyle \frac{1}{∆y_j\left(∆y_j+∆y_{j-1}\right)}}\\ {\alpha }_4={\displaystyle \frac{2{\kappa }_g}{{\rho }_g{Cp}_g}}{\displaystyle \frac{1}{∆y_{j-1}\left(∆y_j+∆y_{j-1}\right)}}\\ {\alpha }_5={\displaystyle \frac{2{\kappa }_g}{{\rho }_g{Cp}_g}}{\displaystyle \frac{1}{∆z_k\left(∆z_k+∆z_{k-1}\right)}}\\ {\alpha }_6={\displaystyle \frac{2{\kappa }_g}{{\rho }_g{Cp}_g}}{\displaystyle \frac{1}{∆z_{k-1}\left(∆z_k+∆z_{k-1}\right)}}\end{array}\right.$$

where $T_{i,j,k}^n$ is the temperature at grid point $(i,j,k)$ in the computational domain at time $t=n$. Moreover, $∆x_i$, $∆y_j$ and $∆z_k$ are the distances to the adjacent grid points in the $x$, $y$, and $z$ directions from the grid point $(i,j,k)$, and $∆t$ is the time increment. The convection term in Eq. (2) was discretized using the values from the previous time $t=n$.

The computational grid is shown in Figure 3. Note that the number of grid points in the computational grid shown in the figure has been reduced for visibility. A finer grid system was adopted at the channel inlet, the interfaces between the fluid and glass, and the interfaces between glass and air. The number of computational grid points used was $~\mathrm{1,300,000}$. This number of computational grids was the minimum number of grids that would not change the results of the simulation. For the heat transfer coefficient $h$ of the horizontal plates corresponding to the upper and lower end surfaces of the device, many empirical formulas were proposed in the past, and several of these are listed in Table 1. Regarding the exponent $\zeta$ of the Rayleigh number Ra in Eqs. (7) and (11), for the low Rayleigh numbers ($\mathrm{R}\mathrm{a}~\le {10}^5$), which were relevant to the thermofluidic phenomena targeted in this analysis, a model with $\zeta ~1/4$ was proposed. Therefore, in this analysis, a model of $h~{\mathrm{R}\mathrm{a}}^{1/4}$ was also adopted. For the coefficients $C_t$ and $C_b$ appearing in Eqs. (7) and (11), their magnitudes were generally on the order of ${10}^{-1}$ [40,41,42,43]. Accordingly, the coefficients $C_t$ and $C_b$ for $h_t$ and $h_b$ were determined as

$$\begin{array}{c}\left\{\begin{array}{c}{C}_t=0.1\\ C_b=0.1\end{array}\right.\end{array}$$

(22)

To begin with the numerical simulation for thermofluid analysis, the electric field distribution was solved first using Eqs. (14)–(20). Then, using Eq. (4), the distribution of $q$ in the cross-section of the microfluidic channel was obtained. Next, using Eqs. (2)–(4), the temperature field in the device was obtained. Time integration was performed with a constant time step until the temperature at the bottom surface of the microfluidic channel became almost constant. The analysis conditions are listed in Table 2. The thermophysical properties used are listed in Table 3. The in-house Fortran parallel-computing code managed using OpenMP library was used for the numerical simulations. Computations were performed on the multi-core parallel simulation system of the National Defense Academy.

Table 1. Empirical formulas proposed for the convective heat transfer coefficient $h$.

Range of $\mathbf{R}\mathbf{a}$	Empirical formula	Refs.
${10}^4\le \mathrm{R}\mathrm{a}\le {10}^7$ (Top)	$h=0.54{\mathrm{R}\mathrm{a}}^{1/4}$	McAdams et al. [40]
${10}^5\le \mathrm{R}\mathrm{a}\le 2\times {10}^7$ (Top) ${3\times 10}^5\le \mathrm{R}\mathrm{a}\le 3\times {10}^{10}$ (Bottom)	$h=0.54{\mathrm{R}\mathrm{a}}^{1/4}$ $h=0.27{\mathrm{R}\mathrm{a}}^{1/4}$	Fishenden et al. [41]
$2\times {10}^5\le \mathrm{R}\mathrm{a}\le 4\times {10}^7$ (Top)	$h=0.70{\mathrm{R}\mathrm{a}}^{1/4}$	Al-Arabi et al. [42]
$1\times {10}^2\le \mathrm{R}\mathrm{a}\le 6\times {10}^3$ (Top)	$h=0.297{\mathrm{R}\mathrm{a}}^{1/4}$	Yousef et al. [43]

Table 1. Empirical formulas proposed for the convective heat transfer coefficient $h$.

Range of $\mathbf{R}\mathbf{a}$	Empirical formula	Refs.
${10}^4\le \mathrm{R}\mathrm{a}\le {10}^7$ (Top)	$h=0.54{\mathrm{R}\mathrm{a}}^{1/4}$	McAdams et al. [40]
${10}^5\le \mathrm{R}\mathrm{a}\le 2\times {10}^7$ (Top) ${3\times 10}^5\le \mathrm{R}\mathrm{a}\le 3\times {10}^{10}$ (Bottom)	$h=0.54{\mathrm{R}\mathrm{a}}^{1/4}$ $h=0.27{\mathrm{R}\mathrm{a}}^{1/4}$	Fishenden et al. [41]
$2\times {10}^5\le \mathrm{R}\mathrm{a}\le 4\times {10}^7$ (Top)	$h=0.70{\mathrm{R}\mathrm{a}}^{1/4}$	Al-Arabi et al. [42]
$1\times {10}^2\le \mathrm{R}\mathrm{a}\le 6\times {10}^3$ (Top)	$h=0.297{\mathrm{R}\mathrm{a}}^{1/4}$	Yousef et al. [43]

Table 2. Values of parameters used in the numerical simulation.

Parameters	Values
Number of grid points	$\mathrm{1,300,000}$
Time increment $∆t$ [s]	$1.0\times {10}^{-2}$
Characteristic length $L$ [m]	${6.0\times 10}^{-2}$
Characteristic temperature $T_0$ [K]	$293$
Effective potential $V_{\mathrm{r}\mathrm{m}\mathrm{s}}$ [V]	$5/\sqrt{2}$
Solution conductivity ${\sigma }_f$ [S/m]	$1.0\times {10}^{-3}-1.0\times {10}^{-1}$
Flow rate $Q$ [mL/h]	$5.0-50.0$

Table 2. Values of parameters used in the numerical simulation.

Parameters	Values
Number of grid points	$\mathrm{1,300,000}$
Time increment $∆t$ [s]	$1.0\times {10}^{-2}$
Characteristic length $L$ [m]	${6.0\times 10}^{-2}$
Characteristic temperature $T_0$ [K]	$293$
Effective potential $V_{\mathrm{r}\mathrm{m}\mathrm{s}}$ [V]	$5/\sqrt{2}$
Solution conductivity ${\sigma }_f$ [S/m]	$1.0\times {10}^{-3}-1.0\times {10}^{-1}$
Flow rate $Q$ [mL/h]	$5.0-50.0$

Table 3. Thermophysical properties used in the numerical simulation (20 $℃$, 1 atm).

	$\mathit{\kappa }$ [W/(m$·$K)]	$\mathit{\rho }$ [kg/m3]	${\mathit{C}}_{\mathit{p}}$ [J/(kg$·$K)]	$\mathit{\nu }$ [m2/s]	$\mathit{\beta }$ [1/K]	Pr
Water	${6.03\times 10}^{-1}$	$9.97{\times 10}^2$	$4.19{\times 10}^3$	$1.00\times {10}^{-6}$	$0.20\times {10}^{-3}$	$7.01$
Air	$2.57{\times 10}^{-2}$	$1.17$	${1.01\times 10}^3$	$1.50\times {10}^{-6}$	$3.41\times {10}^{-3}$	$0.71$
Glass	$1.05$	$2.20{\times 10}^3$	${8.40\times 10}^2$	―	―	―

Table 3. Thermophysical properties used in the numerical simulation (20 $℃$, 1 atm).

	$\mathit{\kappa }$ [W/(m$·$K)]	$\mathit{\rho }$ [kg/m3]	${\mathit{C}}_{\mathit{p}}$ [J/(kg$·$K)]	$\mathit{\nu }$ [m2/s]	$\mathit{\beta }$ [1/K]	Pr
Water	${6.03\times 10}^{-1}$	$9.97{\times 10}^2$	$4.19{\times 10}^3$	$1.00\times {10}^{-6}$	$0.20\times {10}^{-3}$	$7.01$
Air	$2.57{\times 10}^{-2}$	$1.17$	${1.01\times 10}^3$	$1.50\times {10}^{-6}$	$3.41\times {10}^{-3}$	$0.71$
Glass	$1.05$	$2.20{\times 10}^3$	${8.40\times 10}^2$	―	―	―

3. Experimental

3.1. Device Fabrication

Figure 4 (a) shows the schematic of the electrode substrate. The counter-interdigitated electrode measured 36 mm long and 50 µm wide. The interdigitated electrode used in the experiment was an aluminum film with a thickness of 300 nm fabricated by standard photolithography. In brief, a glass plate measuring $50\mathrm{m}\mathrm{m}\times 90\mathrm{m}\mathrm{m}$ with a thickness of 0.4 mm was used as the electrode substrate. The plate surface was ultrasonically cleaned using acetone and isopropyl alcohol for 5 min each. A positive photoresist (S1805G; Rohm and Haas Electronic Materials) was spin-coated on the plate surface. The coated plate was baked at 90 $℃$ for 3 min. The resist layer was exposed to UV light through a positive mask image. The exposed photoresist was developed and baked at 130 $℃$ for 3 min. A 300 nm-thick layer of the aluminum was vacuum-deposited over the plate. The uncovered aluminum area was etched with mix acid at 40 $℃$ for 90 s. Finally, the photoresist was removed using an AZ 100 remover (AZ Electronic Materials). Figure 4 (b) shows the components of the cell separation device. The upper surface was a glass substrate with a thickness of 1.1 mm coated with ITO film. It had $\varphi 1.0$ holes for introducing and draining the cell sample solution, which were made using a router (2307396; Sea Force). Polypropylene female luer fittings were adhered to the perforated parts, and lead wires were attached to the electrode terminals using conductive epoxy adhesive. The microfluidic channel of dimensions 0.5 mm ($H$)$\times$36 mm ($L_z$)$\times$10 mm ($B$) was constructed with parallel top and bottom glass plates separated by 0.5 mm with a silicon rubber spacer.

3.2. Sample Solution Preparation

Rhodamine B, which is widely used in the LIF method, was used as the fluorescent dye. A 300 mM mannitol solution, used for cell separation [36], was used as the solvent for the sample solution and the solution conductivity was prepared to ${\sigma }_f=1\times {10}^{-3}$ S/m using cell culture medium, and the final concentration of Rhodamine B (Sigma 83695-250MG; $ex$. 553 nm, $em$. 627 nm) was set at 20 µM.

3.3. Experimental Setup

A schematic diagram of the primary part of the experimental setup is shown in Figure 5. The experimental setup mainly consisted of the cell separation device, a syringe pump (YSP-201; YMC) and a wave generator (AFG3151; Tektronix). The cell separation device has a parallel plate structure with an ITO glass plate electrode on the top face and counter interdigitated electrodes on the bottom face. The device was set on the stage of a confocal laser-scanning microscope (FV-1000; Olympus) and the microscope focus was set on the surface of the counter-interdigitated electrode at the edges where heat generation was expected to be most pronounced. The top electrode plate and one of the paired counter-interdigitated electrodes were connected to the ground side of the wave generator, and the remaining electrode was connected to the high-voltage side of the wave generator using lead wires. The syringe set on the syringe pump was connected to the sample solution introduction port at the microfluidic upstream end via a silicone tube. After the syringe was filled with 5 mL of sample solution, the sample solution flowed into the microfluidic channel at a constant flow rate and the syringe pump was used to degas the microfluidic channel. After the sample solution flow had stabilized, the wave generator was switched on and an A.C. voltage of $V=10$ V_pp was applied to the device. The temperature of the sample solution was measured using a T–type thermocouple. Fluorescence images were acquired using a laser beam with an excitation wavelength of 559 nm and a bandpass filter with a wavelength range of 570–670 nm. To ensure that the acquired images were free from white-out and black-out, the intensity of the laser light was adjusted. By checking the brightness histogram of the fluorescence image, the fluorescence intensity was confirmed to be within the dynamic range of the detector. Before the electric-field application, the initial temperature of the sample solution was measured and a reference image was acquired at the initial temperature with no voltage application (0 V_pp). After acquiring the reference image, the voltage was applied to the device concurrent with the start of fluorescence image acquisition. Time-lapse fluorescence images were acquired every 30 s from the onset of the voltage application ($t=0$ s) until $t=1500$ s elapsed. The acquired fluorescence images were stored in an uncompressed TIFF format with a resolution of 1024 $\times$1024 pixels. The experiment was performed by setting the frequency $f=57.6$ kHz and the sample solution flow rate $Q=5$ mL/h. The experimental conditions are listed in Table 4.

3.4. Image Processing

The procedure for converting the acquired fluorescence images into images of the distribution of the temperature rise, $∆T$, is shown in Figure 6. The NIH (National Institutes of Health) open-source image-processing software Image J (https://imagej.nih.gov/ij/) was used for image analysis. Time-lapse fluorescence images were sequentially read by Image J and converted into images of the fluorescence intensity ratio by dividing each read image by the reference image acquired before the electric field application. The converted images were then transformed into 8-bit (256 shades of grey scale) text images. Subsequently, the intensity values per image pixel were converted to $∆T$ values using the calibration curve. Finally, the distribution of $∆T$ was visualized in pseudo-color using the graphical drawing software Tecplot (https://www.tecplot.com/). The histogram of image pixels with values of $∆T$ was obtained for each image acquired at each time interval. OpenCV library functions were used to generate the histograms. For each image of the fluorescence intensity ratio before temperature transformation, a non-local mean filter was used to remove noise, after which a histogram of brightness was obtained. The average value of $∆T$ was obtained from the distribution of the obtained histograms. These processes were carried out with Python 3 (ver. 3.9.6; https://www.python.org/) code.

3.5. Temperature Calibration

Calibration of the temperature against the fluorescence intensity ratio was performed using a plate reader (DTX-880; Beckman Coulter). 60 µL of the sample solution was dispensed into 10 wells of a microtiter plate and another 60 µL of 300 mM mannitol solution with solution conductivity ${\sigma }_f=1\times {10}^{-3}$ S/m was dispensed into another 10 wells for background acquisition. The wavelength of the excitation light was set to 535 nm and the wavelength of the bandpass filter to 625 nm, and the temperature in the chamber of the plate reader was changed every 5 $℃$ from 25 $℃$ to 45 $℃$ to acquire the fluorescence intensity. For calibration, the background fluorescence intensity was subtracted from the fluorescence intensity of the sample solution acquired at each temperature. The calibration curve obtained is shown in Figure 7 (a). The number of measurements was $n=20$ and the plots represent $\mathrm{M}\mathrm{E}\mathrm{A}\mathrm{N}\pm 2\mathrm{S}\mathrm{E}$. Fluorescence intensities were normalized to that at a temperature of 25 $℃$. The straight line was linearly interpolated by the least-squares method and the coefficient of determination was R²=0.9924. Figure 7 (b) shows the histogram obtained by imaging analysis. The histogram was a smooth curve with only a single maximum at each time point. The histogram curve became broader with time, while the value of $∆T$, which gave the maximum, increased. $∆T$ was obtained as the sum of the multiplication of the temperature rise ${∆T}_p$ with the number of image pixels having the value of the temperature rise $k$(${∆T}_p$) divided by the total number of image pixels. $∆T$ was defined as

$$\begin{array}{c}∆T={\displaystyle \frac{{\sum }_{p=1}^{\mathrm{M}}{∆T}_p·k\left({∆T}_p\right)}{{\sum }_{p=1}^{\mathrm{M}}k\left({∆T}_p\right)}}\end{array}$$

(23)

where $M$ is the largest number of pixels. The resolution of the temperature rise was ${∆T}_{p+1}-{∆T}_p=0.26$ $℃$.

4. Results and Discussion

4.1. Developed Code Verification

4.1.1. Developed Code Verification

Figure 8 (a) shows the distribution of (left) the electrostatic potential, ${\psi }^{*}$, and (right) the electric field, $E^{*}$, in the cross section of the microfluidic channel. The potential and electric field are non-dimensionalized as ${\psi }^{*}=\psi /V_{\mathrm{r}\mathrm{m}\mathrm{s}}$ and $E^{*}=EH/V_{\mathrm{r}\mathrm{m}\mathrm{s}}$, respectively. The color contours represent the results using the developed simulation code, and dashed monochrome lines represent the results using the commercial software FEATool Multiphysics (https://www.featool.com/multiphysics/). The number of grid points used in the commercial software was $2.6\times {10}^6$. The convergence criterion of the relative error of the iterative calculation was defined as $<{1.0\times 10}^{-8}$. As shown in Figure 8 (a), the electric field was extremely strong in the vicinity of the electrode edges. Thus, the strong electric field, which is highly localized in the microfluidic channel, acts as a strong local heat source. Therefore, it is essential to ensure sufficient spatial resolution in this region to accurately evaluate the temperature field. The electric field obtained by the developed simulation code was in good agreement with that obtained by the commercial software, confirming that the developed simulation code has sufficient accuracy. Figure 8 (b) shows the profile of the local Nusselt number for parallel-plates laminar flow of isothermal walls obtained with the developed simulation code, compared with those reported by Shah and London [44]. The boundary conditions are constant wall temperature ($T=23$ $℃$) and constant inflow temperature ($T=T_0$). The dimensionless $z$-coordinate is defined as

$$\begin{array}{c}{z}^{\mathrm{*}}={\displaystyle \frac{z}{2H\mathrm{R}\mathrm{e}·\mathrm{P}\mathrm{r}}}\end{array}$$

(24)

The Reynolds number is defined as

$$\begin{array}{c}Re={\displaystyle \frac{{2HU}_m}{{\upsilon }_f}}\end{array}$$

(25)

and calculations were performed under the condition of $\mathrm{R}\mathrm{e}=400$. The local Nusselt number is defined as

$$\begin{array}{c}{\mathrm{N}\mathrm{u}}_{z^{\mathrm{*}}}={\displaystyle \frac{2H}{T_{\mathrm{w}}\left(z^{\mathrm{*}}\right)-T_{\mathrm{b}}\left(z^{\mathrm{*}}\right)}}{\left(-{\displaystyle \frac{𝜕T\left(z^{\mathrm{*}}\right)}{𝜕y}}\right)}_{y=0}\end{array}$$

(26)

where $T_{\mathrm{w}}$ is the temperature at the bottom surface of the microfluidic channel, and $T_{\mathrm{b}}$ is the bulk mean temperature of the fluid defined as

$$\begin{array}{c}{T}_{\mathrm{b}}={\displaystyle \frac{{\int }_S{\rho }_f{C}_{pf}T\left(x,y,z^{\mathrm{*}}\right)wdS}{{\int }_S{\rho }_f{C}_{pf}wdS}}\end{array}$$

(26)

where $S$ represents the cross-sectional area of the microfluidic channel in the computational domain. The values of ${\mathrm{N}\mathrm{u}}_{z^{*}}$ obtained by the developed simulation code were in good agreement with that obtained by the previous study, confirming that the developed thermofluid simulation code has sufficient accuracy.

4.1.2. Thermal Structure of the Device

Figure 9 shows the transient of the distribution of temperature rise, $\mathsf{\Delta }T$ ($=T-T_0$), in the computational domain. The analysis conditions were set to $\sigma ={4.0\times 10}^{-2}$ S/m and $Q=5.0$ mL/h to match the experimental conditions of the developed cell separation device [36]. Numerical simulation was performed until $t=1440$ s when the cell separation of a 2 mL sample solution was completed in the experiment. The $\mathsf{\Delta }T$ took the maximum value at the bottom of the microfluidic channel, where the generation of Joule heat was the greatest. The temperature distribution had a gradual increase along the flow direction ($z$) across the liquid phase (microfluidic channel) and solid phases (top and bottom glass substrates). The results imply that the heat transport by diffusion is dominant over that by convection in the device. The maximum of the $\mathsf{\Delta }T$ appeared at the downstream end of the microfluidic channel was predicted to reach $~21$ $℃$ at time $t=1440$ s.

Figure 10 shows time variations of distributions of (a) the temperature rise at the bottom surface of the microfluidic channel, $∆T_{\mathrm{w}}$ ($=T_{\mathrm{w}}-T_0$), along the flow direction, and (b) the local Nusselt number, Nu, along the flow direction, up to $t=1500$ s, for $\sigma =1.0\times {10}^{-2}$ S/m and $Q=10.0$ mL/h. $\mathrm{N}\mathrm{u}$ is defined as

$$\begin{array}{c}Nu=-{\displaystyle \frac{2H}{T_{\mathrm{w}}\left(z\right)-T_{\mathrm{b}}\left(z\right)}}{\left(-{\displaystyle \frac{𝜕T\left(z\right)}{𝜕y}}\right)}_{y=0}\end{array}$$

(28)

Here, $T_w$ and $𝜕T\left(z\right)/𝜕y$ are the values of the local temperature and temperature gradient at the bottom of the microfluidic channel, averaged with respect to the span direction

$$\begin{array}{c}{T}_{\mathrm{w}}\left(z\right)={\displaystyle \frac{\int T_{\mathrm{w}}\left(x,0,z\right)dx}{\int dx}},{\left({\displaystyle \frac{𝜕T\left(z\right)}{𝜕y}}\right)}_{y=0}={\displaystyle \frac{\int {\left(𝜕T\left(x,y,z\right)/𝜕y\right)}_{y=0}dx}{\int dx}}\end{array}$$

(29)

Unless otherwise specified, all distributions of thermophysical quantities in the flow direction represent the spanwise averaged distribution. As shown in Figure 10 (a), because the Joule heat was generated almost uniformly along the flow direction over the bottom surface of the microfluidic channel, $\mathsf{\Delta }{T}_{\mathrm{w}}$ began to increase uniformly in the flow direction as soon as the voltage was applied. On the other hand, the upstream end of the microfluidic channel was cooled by the inflow fluid. Therefore, while $\mathsf{\Delta }{T}_{\mathrm{w}}$ increased uniformly with time, a distribution was presented in which the high-temperature plateau gradually narrowed from the upstream side due to the cooling effect of the convection heat transport of the inflow fluid. As shown in Figure 10 (b), the area near to the inlet of the microfluidic channel, $T_{\mathrm{w}}-T_{\mathrm{b}}$, presented the maximum value due to Joule heating, and the Nusselt number, $\mathrm{N}\mathrm{u}$, presented the minimum value. On the other hand, $T_{\mathrm{b}}$ did not increase because the area near the inlet continued to be cooled. Therefore, the minima of $\mathrm{N}\mathrm{u}$ became even smaller as time passed. On the downstream side of the microfluidic channel, $\mathrm{N}\mathrm{u}$ gradually increased uniformly due to Joule heating. Compared to the distribution of $\mathrm{N}\mathrm{u}$ shown in Figure 8 (b), which was obtained under the condition of constant wall temperature, the distribution obtained in this study was different. In the thermofluidic field under study, the fluid in the microfluidic channel was found to be heated spatially nonuniformly by Joule heating due to the presence of the nonuniform electric field. These are the characteristics of the temperature field within the cell separation device.

4.1.3. Influence of Flow Rate and Conductivity on Temperature

Figure 11 shows (top) the distribution of the local temperature rise at the bottom surface of the microfluidic channel, $\mathsf{\Delta }{T}_{\mathrm{w}}$, along the flow direction, and (bottom) the time variation of the mean temperature rise, $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$, for varying $\sigma$ and $Q$. $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$ is defined as

$$\begin{array}{c}\Delta T_{\mathrm{w}\mathrm{m}}\left(t\right)={\displaystyle \frac{\int \mathsf{\Delta }{T}_{\mathrm{w}}\left(t\right)dz}{\int dz}}\end{array}$$

(30)

For $\mathsf{\Delta }{T}_{\mathrm{w}}$, the distribution for the thermal steady state ($t=1500$ s) is shown, and for $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$, the transient variation is shown up to $t=600$ s, when the value changed significantly. In Figure 11 (a), $Q$ was fixed at $Q=10.0$ mL/h, and $\sigma$ was varied in the range $\sigma =1.0\times {10}^{-3}-1.0\times {10}^{-1}$ S/m. In Figure 11 (b), $\sigma$ was fixed at $\sigma =1.0\times {10}^{-2}$ S/m, and $Q$ was varied in the range $Q=5.0-20.0$ mL/h. As shown in the top panel of Figure 11 (a), the distribution of $\mathsf{\Delta }{T}_{\mathrm{w}}$ increased proportionally with $\sigma$. This is because the heat generation term $q$ on the right-hand side of Eq. (4) is proportional to $\sigma$. As shown in the top panel of Figure 11 (b), the effect of convective heat transfer became significant with increasing $Q$, and the shape of the $\mathsf{\Delta }{T}_{\mathrm{w}}$ profile transitioned from a parabolic to a linear form. From the above, it was found that $Q$ is the dominant parameter for the shape of the distribution of $\mathsf{\Delta }{T}_{\mathrm{w}}$ in the microfluidic channel. Therefore, in order to maintain a uniform temperature distribution in the microfluidic channel, which is desirable for cell separation devices, it is essential to minimize the value of $Q$ as much as possible. In this study, the heat exchange by conduction was not considered as a thermal boundary condition at the downstream end of the device. Therefore, under actual cell separation operations, the temperature at the downstream end of the device is expected to be affected by the external thermal environment at the downstream end of the device.

Regarding the $\sigma$ dependence of the transient behavior of $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$, $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$ reached a steady state at$t~100$ s regardless of the value of $\sigma$, as shown in the bottom panel of Figure 11 (a). After reaching a steady state, the level of $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$ increased in proportion to the value of $\sigma$, as did $\mathsf{\Delta }{T}_{\mathrm{w}}$, in the top panel of Figure 11 (a). Regarding the $Q$ dependence of the transient behavior of $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$, the effect of convective heat transport increased with increasing $Q$ as shown in the bottom panel of Figure 11 (b). Therefore, the time until $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$ reaches a steady state became shorter. On the other hand, the level of $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$ decreased with increasing $Q$ due to the cooling effect of the incoming fluid flow. In addition, when $Q$ varied in the range of $Q=5-20$ mL/h, the Peclet number, $\mathrm{P}\mathrm{e}$, varied in the range of approximately $\mathrm{P}\mathrm{e}=2.0-7.8$. Because of the small change and small order of magnitude of $\mathrm{P}\mathrm{e}$, the effect of $Q$ on $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$ was about $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}~1$ $℃$, and thus hardly appeared.

Figure 12 shows a bird’s-eye view of $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$ as a function of both $\sigma$ and $Q$ after the temperature field reached a steady state. Parametric calculations ($~200$ simulations) were performed with $\sigma$ varying in the range of $\sigma =1.0\times {10}^{-3}-1.0\times {10}^{-1}$ S/m, and $Q$ in the range of $Q=5.0-50.0$ mL/h. As shown in Figure 12, the value of $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$ changed rapidly in the range of $\sigma =1.0\times {10}^{-2}-1.0\times {10}^{-1}$ S/m regardless of the value of $Q$. Moreover, at higher values of $\sigma$, it was found that the value of $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$ varied significantly with variations in $Q$. Under these conditions, the effect of convective heat transport on the temperature rise was pronounced due to the high heat generation, resulting in a distinct $Q$-dependence of $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$. The maximum value of $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$ reached approximately 48 ℃. Elengoe et al. reported that the viability of MCF10A and MDA-MB-231 cells have been decreased by exposure to a 41 $℃$ environment for 30 min [45]. Furthermore, it has been reported that exposure to high-temperature environments induces adverse effects on cells, including cell death and reduced viability [19,20,21]. Therefore, the highest value of the predicted $\mathsf{\Delta }{T}_{\mathrm{w}\mathrm{m}}$ is considered to have a non-negligible influence on cells during the cell separation process.

4.1.4. Influence of Flow Rate and Conductivity on the Nusselt Number

Figure 13 shows (top) the distribution of the local Nusselt number at the bottom surface of the microfluidic channel, $\mathrm{N}\mathrm{u}$, along the flow direction, and (bottom) the time variation of the mean Nusselt number, ${\mathrm{N}\mathrm{u}}_{\mathrm{m}}$, for varying $\sigma$ and $Q$. ${\mathrm{N}\mathrm{u}}_{\mathrm{m}}$ is defined as

$$\begin{array}{c}{\mathrm{N}\mathrm{u}}_{\mathrm{m}}\left(t\right)={\displaystyle \frac{\int \mathrm{N}\mathrm{u}\left(t\right)dz}{\int dz}}\end{array}$$

(31)

For $\mathrm{N}\mathrm{u}$, the distribution for the thermal steady state ($t=1500$ s) is shown, and for ${\mathrm{N}\mathrm{u}}_{\mathrm{m}}$, the transient variation is shown up to $t=600$ s, when the value changed significantly. In Figure 13 (a), $Q$ was fixed at $Q=10.0$ mL/h, and $\sigma$ was varied in the range of $\sigma =1.0\times {10}^{-3}-1.0\times {10}^{-1}$ S/m. In Figure 13 (b), $\sigma$ was fixed at $\sigma =1.0\times {10}^{-2}$ S/m, and $Q$ was varied in the range of $Q=5.0~20.0$ mL/h. As shown in the top panel of Figure 13 (a), the distribution of $\mathrm{N}\mathrm{u}$ exhibited the minimum value near the inlet of the microfluidic channel regardless of the value of $\sigma$. Then, due to the heat balance between Joule heating and heat diffusion from the bottom surface of the microfluidic channel, $\mathrm{N}\mathrm{u}$ gradually increased along the flow and eventually asymptotically approached a constant value. The position from the inlet of the microfluidic channel where $\mathrm{N}\mathrm{u}$ asymptotically approached a constant value moved downstream in proportion to the value of $\sigma$, i.e., the amount of heat generated. In other words, the value of $z$ at which the temperature field reaches thermal equilibrium shifted downstream as $\sigma$ increased. In the bottom panel of Figure 13 (a), for a large value of $\sigma$, ${\mathrm{N}\mathrm{u}}_{\mathrm{m}}$ exhibited a minimum value for a short while ($<~20$ s) after the onset of the voltage application. This is because, in cases of large heat generation, the $T_{\mathrm{b}}$ was initially low, resulting in a large $T_{\mathrm{w}}-T_{\mathrm{b}}$ near the inlet. As the temperature rose, heat dissipation from the bottom surface of the microfluidic channel increased, causing ${\mathrm{N}\mathrm{u}}_{\mathrm{m}}$ to gradually increase and eventually reach a constant value when the device reached thermal equilibrium. The top panel of Figure 13 (b) shows the dependence of the Nu distribution on $Q$. As the value of $Q$ increased, the value of $\mathrm{N}\mathrm{u}$ decreased. This is because the temperature rise in the microfluidic channel was suppressed due to the inflow of fluid at room temperature from the inlet. In this study, heat conduction at the downstream end of the microfluidic channel is not considered in the thermal boundary condition. Therefore, ${\mathrm{N}\mathrm{u}}_{\mathrm{m}}$ decreased simply in accordance with the increase in $Q$ as shown in the bottom panel of Figure 13 (b). The convective heat transport due to the increase in $Q$ would suppress the temperature rise in the microfluidic channel.

Figure 14 shows a bird’s-eye view of ${\mathrm{N}\mathrm{u}}_{\mathrm{m}}$ as a function of both $\sigma$ and $Q$ after the temperature field reached a thermal steady state. In the range of $\sigma =1.0\times {10}^{-2.2}-1.0\times {10}^{-1}$ S/m, ${\mathrm{N}\mathrm{u}}_{\mathrm{m}}$ transitioned smoothly with changes in $Q$, while in the range of $\sigma =1.0\times {10}^{-3}-1.0\times {10}^{-2.2}$ S/m, ${\mathrm{N}\mathrm{u}}_{\mathrm{m}}$ remained nearly flat regardless of the value of $Q$. The maximum value of ${\mathrm{N}\mathrm{u}}_{\mathrm{m}}$ appeared when $\sigma$ was the maximum and $Q$ was the minimum, while the minimum value of ${\mathrm{N}\mathrm{u}}_{\mathrm{m}}$ appeared when $\sigma =1.0\times {10}^{-2.2}$ and $Q=50.0$ mL/h. The results suggest that $\sigma ~1.0\times {10}^{-2.2}$ S/m may be the threshold for a significant change in the thermal structure in the microfluidic channel of the present cell separation device.

4.2. Experimental Results and Comparison with Numerical Simulation

Figure 15 shows the transient of the distribution of the temperature rise, $∆T$, at the bottom surface of the microfluidic channel within a $100\mathsf{\mu }\mathrm{m}\times 100\mathsf{\mu }\mathrm{m}$ area in the center of the device, displayed in pseudo-color. The top panel represents the numerical simulation results and the bottom panel represents the experimental results. In the numerical simulation, the $∆T$ increased significantly until $t\cong 120$ s, after which it hardly increased at all. In contrast, the experimental results showed a gradual increase in $∆T$ up to$t\cong 1440$ s. Moreover, the experimental results showed some temperature difference between the electrode and glass substrate, whereas the simulation results showed uniform $∆T$ distributions. The reason for this difference may be attributed to the fact that the numerical simulation model ignored the thickness of the metal electrode (300 nm). The difference in the heat capacity between the metal and glass may be responsible for the difference in the rate of temperature rise. However, since heat transport by conduction was dominant in microfluidic channels, there was actually little temperature difference between the electrode surface and the glass surface, as shown in the experimental results. Figure 16 shows the transient of the averaged $\mathsf{\Delta }T$ with respect to the $100\mathsf{\mu }\mathrm{m}\times 100\mathsf{\mu }\mathrm{m}$ area in the center of the bottom surface of the microfluidic channel, comparing the experimental and numerical simulation results. In the experiment, as can be seen from the transient behavior of the $\mathsf{\Delta }T$ in Figure 16, the temperature field in the device nearly reached a steady state in $t\cong 720$ s, whereas in the numerical simulation, it reached a steady state in $t\cong 240$ s. The difference in the time history of the $\mathsf{\Delta }T$ between the two may be due to the fact that the numerical simulation did not take into account the effect of heat dissipation from both lateral sides of the device (relative to the flow direction). Another factor to consider is the validity of the values for $h_b$ and $h_t$. In this study, the values of $h_b$ and $h_t$ were estimated using reasonable empirical formulas that took into account the length scales under analysis and were roughly of the order of magnitude of 2–3 digits [46]. On the other hand, in previous reports on numerical simulation of Joule heating within the microfluidic channel [23,24], the value of the heat transfer coefficient $h$ was significantly higher than what is typically expected from natural convection of air, ranging from $20-20000$ W/(m²∙K). A value of $h~20000$ W/(m²∙K) is on the same order of magnitude as that observed in boiling heat transport phenomena. The setting of various parameter values used for boundary conditions could potentially have a non-negligible effect on the evaluation of the temperature field. Nakano et al. [24] have pointed out that the results of temperature field analysis through numerical simulation largely depend on the determination of $h$ values. In the modeling of thermofluid analysis, it is believed that physically reasonable assumptions are crucial for constructing models which contribute to the accurate understanding of phenomena.

5. Conclusions

In recent years, dielectrophoresis-based cell separation technologies have attracted attention as a method that is noninvasive to cells. On the other hand, there are concerns about the effects of the high-temperature environment inside the device caused by Joule heating due to the strong electric field on the physiological functions of cells. In this study, we investigated the temperature field inside the microfluidic channels of a previously developed cell separation device through both experimental and numerical approaches. In the experiment, the micro-LIF method was adopted to measure the temperature in the device. In the numerical approach, the thermal structure of the device was investigated by performing a numerical simulation in which the effect of heat generation due to a nonuniform electric field was considered in the model. The experimental results revealed that the temperature rise at the bottom surface of the microfluidic channel was approximately 20 $℃$ after 1440 s of voltage application. The transient of the thermofluidic phenomena in the device was investigated by numerical simulation, and an attempt was made to elucidate the thermal structure of the device. Additionally, a parametric analysis was conducted to investigate the correlation between the flow rate of the cell sample solution or solution conductivity and the temperature field. The results suggested that, under certain operational conditions, the temperature rise from room temperature could exceed 40 $℃$. Our findings suggested that under certain operating conditions, the device could significantly affect the physiological functions of cells. To achieve both an improvement in device performance and a reduction of stress for the cells, detailed exploration of the optimal operating conditions is indispensable.