Numerical simulation of heat transfer in microchannel using energy conservative dissipative particle dynamics

Abstract

Dissipative particle dynamics is a particle-based mesoscopic simulation method. Classic dissipative particle dynamics cannot be used to simulate heat transfer in fluids since the total energy of the system is not conserved. In this article, two-dimensional unsteady heat conduction is first investigated using dissipative particle dynamics with energy conservation. The energy conservative dissipative particle dynamics results are compared with the FLUENT simulation data, and it demonstrates that they are in good agreement with each other. Then, forced convection heat transfers in microchannel of the same wall temperature and different wall temperatures are simulated, respectively, by using periodic boundary condition of dimensionless temperature. The results show that the velocity, temperature, and dimensionless temperature distributions are consistent with theoretical results. Finally, we give a qualitative analysis about the applicability of the energy conservative dissipative particle dynamics approach in simulating flow and heat transfer in rough microchannel.

Introduction

In recent years, with the rapid development of microelectro mechanical systems (MEMS), the demand for probing the mechanism of the microscale flow and heat transfer problems is increased, in which the Knudsen number is probably larger or equal to unity. This restricts the applicability of continuum mechanics in these types of problems. The most straightforward and accurate approach is probably the molecular dynamics (MD). Sofos et al. ¹ studied how wall properties control diffusion in grooved nanochannel by using MD method and achieved some valuable conclusions. However, MD is too inefficient for mesoscale simulations due to the vast computing resources required. ² A number of “coarse-grained” approaches have been suggested and developed to simplify the underlying microscopic model while retaining the essential physics. Although Brownian dynamics simulations (BDS), lattice gas automata (LGA), and lattice Boltzmann method (LBM) are mesoscale simulation methods, it is difficult for BDS to deal with a complex flow field and for LGA and LBM to cope with complex fluids. ³ Dissipative particle dynamics (DPD) has emerged as a promising new mesoscale simulation technique introduced by Hoogerbrugge and Koelman. ⁴ It is a coarse-grained, particle-based simulation method. ⁵ In the DPD system, the basic unit is a set of discrete momentum carriers called particles, moving in continuous space and in discrete time steps. The momentum carrier is a coarse-grained entity of mass m in a three-dimensional simulation domain. The particle’s motion represents the collective dynamic behavior of a large number of molecules (a fluid “particle”). ⁶ DPD conserves not only the number of particles but also the total momentum of the system and satisfies Galilean invariance and the detailed balance equations. Since its introduction, DPD has been widely used for simulations of various complex fluids and complex flows, such as porous flow, ³ colloidal suspensions,^7–9 polymer suspensions, ¹⁰ or multicomponent flows. ¹¹

However, classic DPD can only be applied to simulate isothermal systems, and modeling of thermal systems was not possible since the total energy of the classic DPD system is not conserved. ¹² In order to solve this energy conservation drawback, Español ¹² and Avalos and Mackie ¹³ introduced an internal energy variable and a temperature for each particle in the DPD system, which enabled the modeling of heat transfer. This DPD model with energy conservation is called energy conservative dissipative particle dynamics (eDPD). Each particle in eDPD system is associated with its internal energy in addition to other quantities such as mass, position, and momentum.

The eDPD has been used to model a few of heat transfer problems. By using eDPD, Ripoll et al. ¹⁴ studied the model of eDPD for the simple case of thermal conduction, and it is shown that the model displays correct equilibrium fluctuations and reproduces Fourier law. Chaudhri and Lukes ¹⁵ presented a multicomponent framework of the eDPD; the framework was validated for the special case of identical components in one and two dimensions. Abu-Nada ¹⁶ used eDPD to investigate conduction heat transfer in two dimensions under steady-state condition. Various types of boundary conditions were implemented to the conduction domain. It was found that DPD appropriately predicts the temperature distribution in the conduction regime. Abu-Nada ¹⁷ also used eDPD to study natural convection via Rayleigh–Benard problem and a differentially heated enclosure problem. The eDPD results were compared with the finite volume solutions, and it was found that the eDPD method predicts the temperature and flow fields throughout the natural convection domains properly. The eDPD model recovered the basic features of natural convection. Abu-Nada ¹⁸ further studied natural convection in liquid domain over a wide range of Rayleigh numbers. By using eDPD, Yamada et al. ¹⁹ simulated forced convection in parallel-plate channels with boundary conditions of constant wall temperature and constant wall heat flux. The periodic boundary condition of the dimensionless temperature in the flow direction was adopted. The results agreed well with the exact solutions to within 2.3%. Yamada et al. ²⁰ also applied eDPD to heat conduction in nanocomposite. The effects of the interfacial thermal resistance and the Brownian motion of nanoparticles were incorporated in the model. The results agreed well with the experimental results. Qiao and He^21,22 simulated heat conduction in nanocomposite materials and nanofluids by using eDPD. Their simulation results indicate that the eDPD method can be a versatile method for studying thermal transport in heterogeneous materials and complex systems. More recently, Cao et al. ²³ used eDPD to simulate natural convection in eccentric annulus. A numerical strategy is presented for dealing with irregular geometries in DPD system. And it was found that the temperature and flow fields for the natural convection in complex geometries are correctly predicted using eDPD. Cao et al. ²⁴ also adopted eDPD to simulate the mixed convection in eccentric annulus. The results show that the forced, natural, and mixed convection flow and heat transfer in complex geometries are correctly predicted using eDPD. Furthermore, Li et al. ²⁵ verified the eDPD model for reproducing correctly the temperature-dependent properties including Schmidt number and Prandtl number and analyzed the relationships between the dynamic properties of eDPD fluid and the parameters in the dissipative and random forces. Moreover, Willemsen et al. ²⁶ used the eDPD method to model phase change with the consideration of energy equation.

In the classic DPD method, particles may penetrate the wall because of the soft potentials used. Several methods have been proposed for the wall boundary condition to reintroduce these particles back. In general, there exist two main approaches. In the first approach, the fluid–wall interactions are represented by certain effective forces with the combination of proper reflections to prevent particle penetration. ²⁷ The second and the most widely used approach is based on the representation between the fluid and the wall particles that combine freezing particles along the walls of the physical domain with different wall reflection laws, for example, simple fluid, ²⁸ colloidal suspension, ⁸ and a polymer solution between bounded walls. ²⁹ No matter which methods they take, the aim is to achieve the no-slip boundary condition. In this work, we implement the no-slip boundary condition to the fluid flow in microchannel.

Although quite a number of investigations have been carried out to study the heat transfer problems by using eDPD, it is still necessary to further examine the potential of the eDPD method. In this study, the eDPD was used to simulate fluid flow and forced convective heat transfer in microchannel. The article is organized as follows. In section “Numerical method,” we give a brief review of eDPD schemes. In section “Results and discussion,” we present the validation of the eDPD model by the simulation of two-dimensional (2D) heat conduction. Numerical results for forced convection in microchannel channels are also presented. We conclude with a brief summary in section “Conclusion.”

Numerical method

eDPD

Similar to the classic DPD method, each eDPD particle is still considered as a coarse-grained particle, which represents a cluster of actual fluid molecules instead of atom/molecule directly. ²⁵ The time-evolution of eDPD particles is governed by conservation of momentum and energy. For a simple eDPD particle $i$ ¹⁶

$$\frac{d{\overrightarrow{r}}_i}{\mathrm{dt}}={\overrightarrow{v}}_i$$

(1)

$$\frac{d{\overrightarrow{v}}_i}{\mathrm{dt}}={\overrightarrow{f}}_i={\displaystyle \sum _{i\ne j}}({\overrightarrow{f}}_{\mathrm{ij}}^C+{\overrightarrow{f}}_{\mathrm{ij}}^D+{\overrightarrow{f}}_{\mathrm{ij}}^R)$$

(2)

$$C_v\frac{{\mathrm{dT}}_i}{\mathrm{dt}}=q_i={\displaystyle \sum _{i\ne j}}(q_{\mathrm{ij}}^{\mathrm{visc}}+q_{\mathrm{ij}}^{\mathrm{cond}}+q_{\mathrm{ij}}^{\mathrm{rand}})$$

(3)

where ${\overrightarrow{r}}_i$ and ${\overrightarrow{v}}_i$ are its position and velocity vector, and ${\overrightarrow{f}}_i$ is the interparticle force on particle $i$ by all of other particles (except itself). $T_i$ is the temperature, $C_v$ is the heat capacity at constant volume for eDPD particle, and $q_i$ is the heat flux between particles. The unit of mass is taken to be the mass of a particle so that the force acting on a particle equals its acceleration. The sum runs over all other particles within a certain cutoff radius $r_C$ . Three interparticle forces act upon the particles, and these are the conservative force ${\overrightarrow{f}}_{\mathrm{ij}}^C$ , dissipative force ${\overrightarrow{f}}_{\mathrm{ij}}^D$ , and random force ${\overrightarrow{f}}_{\mathrm{ij}}^R$ ⁴

$${\overrightarrow{f}}_{\mathrm{ij}}^C=a_{\mathrm{ij}}{\omega }_C\left(r_{\mathrm{ij}}\right){\overrightarrow{e}}_{\mathrm{ij}}$$

(4)

$${\overrightarrow{f}}_{\mathrm{ij}}^D=-{\gamma }_{\mathrm{ij}}{\omega }_D\left(r_{\mathrm{ij}}\right)({\overrightarrow{e}}_{\mathrm{ij}}·{\overrightarrow{v}}_{\mathrm{ij}}){\overrightarrow{e}}_{\mathrm{ij}}$$

(5)

$${\overrightarrow{f}}_{\mathrm{ij}}^R={\sigma }_{\mathrm{ij}}{\omega }_R\left(r_{\mathrm{ij}}\right){\xi }_{\mathrm{ij}}\Delta t^{-1/2}{\overrightarrow{e}}_{\mathrm{ij}}$$

(6)

$a_{\mathrm{ij}}$ is a maximum repulsion between particles i and j. ${\gamma }_{\mathrm{ij}}$ is a coefficient which controls the extent of dissipation in a simulation time step, and ${\sigma }_{\mathrm{ij}}$ is a coefficient characterizing the strength of the random forces. For detailed information, one can refer to Chen et al. ⁶

The energy equation (3) includes three parts, in which viscous heat flux $q_i^{\mathrm{visc}}$ accounts for the viscous heating resulting from the dissipative interactions between the particles, collisional heat flux $q_i^{\mathrm{cond}}$ corresponds to the heat transfer resulting from a temperature difference, and random heat flux $q_i^{\mathrm{rand}}$ takes into account the fluctuations due to random heat fluxes. They are expressed as follows ¹²

$$q_i^{\mathrm{visc}}=\frac{1}{2}{\displaystyle \sum _{i\ne j}}\{{\omega }_D\left(r_{\mathrm{ij}}\right)[{\gamma }_{\mathrm{ij}}{({\overrightarrow{e}}_{\mathrm{ij}}·{\overrightarrow{v}}_{\mathrm{ij}})}^2-\frac{{\left({\sigma }_{\mathrm{ij}}\right)}^2}{m}]-{\sigma }_{\mathrm{ij}}{\omega }_R\left(r_{\mathrm{ij}}\right)({\overrightarrow{e}}_{\mathrm{ij}}·{\overrightarrow{v}}_{\mathrm{ij}})\Delta t^{-\frac{1}{2}}{\xi }_{\mathrm{ij}}\}$$

(7)

$$q_i^{\mathrm{cond}}={\displaystyle \sum _{i\ne j}}{k}_{\mathrm{ij}}{\omega }_{\mathrm{cond}}\left(r_{\mathrm{ij}}\right)(\frac{1}{T_i}-\frac{1}{T_j})$$

(8)

$$q_i^{\mathrm{rand}}={\displaystyle \sum _{i\ne j}}{\alpha }_{\mathrm{ij}}{\omega }_{\mathrm{rand}}\left(r_{\mathrm{ij}}\right)\Delta t^{-\frac{1}{2}}{\xi }_{\mathrm{ij}}^e$$

(9)

where ${\overrightarrow{r}}_{\mathrm{ij}}={\overrightarrow{r}}_i-{\overrightarrow{r}}_j$ , $r_{\mathrm{ij}}=\left|{\overrightarrow{r}}_{\mathrm{ij}}\right|$ , ${\overrightarrow{e}}_{\mathrm{ij}}$ is the unit vector, and ${\overrightarrow{v}}_{\mathrm{ij}}={\overrightarrow{v}}_i-{\overrightarrow{v}}_j$ . ${\omega }_C$ , ${\omega }_D$ , ${\omega }_R$ , ${\omega }_{\mathrm{cond}}$ , and ${\omega }_{\mathrm{rand}}$ are the weight functions. $C_v$ is the heat capacity at constant volume for DPD particle. The heat capacity $C_v$ of DPD particles is usually given as a dimensionless form, $\overline{C_v}=C_v/k_B$ , where $k_B$ is the Boltzmann constant. The random numbers ${\xi }_{\mathrm{ij}}$ and ${\xi }_{\mathrm{ij}}^e$ that appear in equations (6), (7), and (9) have a zero mean and unit variance. Each pair of interaction particles has a corresponding value of ${\xi }_{\mathrm{ij}}$ and ${\xi }_{\mathrm{ij}}^e$ , in which ${\xi }_{\mathrm{ij}}={\xi }_{\mathrm{ji}}$ and ${\xi }_{\mathrm{ij}}^e=-{\xi }_{\mathrm{ji}}^e$ . This can ensure that the momentum and energy of each interacting pair of particles are conserved.

All of the weight functions decrease monotonically with the increase in the particle–particle separation distance and vanish for $r<r_C$ . In this study, linear weighting functions are adopted, which are given as

$${\omega }_C\left(r_{\mathrm{ij}}\right)=\{\begin{array}{cc}\hfill 1-\frac{r_{\mathrm{ij}}}{r_C}\hfill & \hfill r_{\mathrm{ij}}<r_C\hfill \\ \multicolumn{1}{c}{0}& \hfill r_{\mathrm{ij}}>r_C\hfill \end{array}$$

(10)

$${\omega }_D\left(r_{\mathrm{ij}}\right)={\omega }_D{\left(r_{\mathrm{ij}}\right)}^2=\{\begin{array}{cc}\hfill {(1-\frac{r_{\mathrm{ij}}}{r_C})}^s\hfill & \hfill r_{\mathrm{ij}}<r_C\hfill \\ \multicolumn{1}{c}{0}& \hfill r_{\mathrm{ij}}>r_C\hfill \end{array}$$

(11)

In order to improve the viscosity of the fluid, $s=1/2$ is adopted. ⁶

The weight functions ${\omega }_{\mathrm{cond}}(r_{\mathrm{ij}})$ and ${\omega }_{\mathrm{rand}}(r_{\mathrm{ij}})$ in equations (8) and (9) are given as

$${\omega }_{\mathrm{cond}}\left(r_{\mathrm{ij}}\right)={\omega }_{\mathrm{rand}}{\left(r_{\mathrm{ij}}\right)}^2=\{\begin{array}{cc}\hfill {(1-\frac{r_{\mathrm{ij}}}{r_C})}^{s_T}\hfill & \hfill r_{\mathrm{ij}}<r_C\hfill \\ \multicolumn{1}{c}{0}& \hfill r_{\mathrm{ij}}>r_C\hfill \end{array}$$

(12)

where $s_T$ is the exponent of the weight functions. The case of $s_T=2$ corresponds to the typical quadratic weighting function, which is adopted in this work following the choices of Abu-Nada. ¹⁶ The parameter $a_{\mathrm{ij}}$ , ${\gamma }_{\mathrm{ij}}$ , and ${\sigma }_{\mathrm{ij}}$ in equations (4) and (7) are the amplitude of the conservative, dissipative, and random forces, respectively. The $k_{\mathrm{ij}}$ and ${\alpha }_{\mathrm{ij}}$ in equations (8) and (9) control the strength of the collisional and random heat flux, respectively. The relation between the parameters ${\gamma }_{\mathrm{ij}}$ and ${\alpha }_{\mathrm{ij}}$ is governed by the fluctuation-dissipation theorem ¹²

$${\gamma }_{\mathrm{ij}}=\frac{{\sigma }_{\mathrm{ij}}^2(T_i+T_j)}{4k_B{T}_i{T}_j}$$

(13)

$${\alpha }_{\mathrm{ij}}=\sqrt{2k_B{k}_{\mathrm{ij}}}$$

(14)

where $k_B$ is the Boltzmann constant. Parameter $k_{\mathrm{ij}}$ is given as

$$k_{\mathrm{ij}}=\frac{C_v^2{k}_0{(T_i+T_j)}^2}{4k_B}$$

(15)

where $k_0$ is interpreted as heat friction that controls thermal conductivity.^21,30 For simplicity of presentation, the mass of DPD particles and the cutoff radius are set to unity in the following simulations.

Initial and boundary conditions

The simple eDPD particles are initially generated at the face-centered cubic (FCC) lattice cites, and it is expected that the density of the fluid is 4. The initial velocities of the fluid particles are set randomly according to the given temperature, but the wall particles are frozen.

For the nanoscales, since the continuum assumption is not satisfied in confined problems, the no-slip boundary condition is not valid. ³¹ In this work, we used eDPD to model the flow and heat transfer in mesoscale; similar to a few of cases, ³² the no-slip boundary condition for the velocity was imposed. As shown in Figure 1, in order to implement no-slip boundary condition and prevent particles from penetrating the walls, the solid walls are represented by using three layers of frozen particles according to Fan et al., ³² in which a thin layer near the wall is assumed, and particles in this layer always leave the wall. Additionally, we enforce a random velocity distribution in this layer with zero mean corresponding to a given temperature.

Figure 1.

Sketch of boundary condition.

In order to mimic an infinite system, periodic flow boundary conditions are applied on the fluid boundary of the computational domain in the x and y directions. However, for the periodic thermally developed region, the dimensionless temperature field repeats itself periodically. The formula of dimensionless temperature for the case of constant wall temperature boundary condition is shown below ¹⁹

$$\Theta =\frac{T-T_w}{\overline{T}-T_w}$$

(16)

where $T_w$ is the wall temperature. The temperature can be derived using the following equation

$$T=\Theta \times (\overline{T}-T_w)+T_w$$

(17)

where the average temperature $\overline{T}$ is given by

$$\overline{T}=\frac{\int {\mathrm{Tv}}_x\mathrm{dy}}{\int v_x\mathrm{dy}}$$

(18)

For the periodic thermally developed region, the dimensionless temperature satisfies the following relationship

$${\Theta }_{\mathrm{cwt}}\left(x\right)={\Theta }_{\mathrm{cwt}}(x+L)=\dots$$

(19)

When a particle moves out of the simulating domain at one side, the average temperature of this side cross section can be calculated by equation (18). Substituting the average temperature into equation (16), the dimensionless temperature $\Theta$ can be calculated. Meanwhile, the average temperature of cross section of inlet side ${\overline{T}}_{\mathrm{inlet}}$ can be calculated. Then, $\Theta$ and ${\overline{T}}_{\mathrm{inlet}}$ are substituted into equation (17), and the updated temperature of the particle is determined and can be assigned to the mirror particle at inlet side.

Results and discussion

Validation of the methodology

Similar to Abu-Nada, ¹⁶ in this study, the eDPD model is validated by simulating a 2D heat conduction problem, and the results are compared with the FLUENT simulation results. Here, we simulated the 2D unsteady heat conduction, while in Abu-Nada, ¹⁶ the steady heat conduction was studied. The schematic of the solution domain and the corresponding boundary conditions of this problem are shown in Figure 2. Three sides of the rectangular plate are maintained at a constant high temperature $T_H=500\mathrm{DPD}\mathrm{units}$ , while the fourth side is set at a constant low temperature $T_C=300\mathrm{DPD}\mathrm{units}$ . Initially, all of the eDPD particles are assigned an initial temperature of 400 DPD units. The simulation system consists of 6085 particles initially randomly distributed in a box of 40 × 30 (x × y). The time step is kept at $2.5\times {10}^{-6}$ . The statistical averaging of the eDPD parameters (temperature in this case) is carried out in the last 800 time steps. Time-evolution of the temperature profile at the mid-slab height (y = 0) and at the mid-slab width (x = 0) is shown in Figure 3. At 17,000 steps, the heat conduction reaches a steady state. Figure 4 shows the comparison of the steady-state temperature isotherms between the eDPD and FLUENT results, and the results agreed well with each other within 0.34%.

Figure 2.

Schematic of the computational domain for 2D plate conduction.

Figure 3.

(a) Time-evolution of the temperature profile at the mid-slab height (y = 0) and (b) time-evolution of the temperature profile at the mid-slab width (x = 0) (solid lines: FLUENT; symbols: eDPD).

Figure 4.

Steady-state temperature isotherms (solid lines: FLUENT; dashed lines: eDPD).

Forced convection heat transfer in microchannel with smooth surfaces

For the simulation of forced convection heat transfer in microchannel, the simulation parameters used are given in Table 1. Note that $a_{\mathrm{ij}}$ and ${\sigma }_{\mathrm{ij}}$ are constant, while the other DPD parameters ( ${\gamma }_{\mathrm{ij}}$ , $k_{\mathrm{ij}}$ , and ${\alpha }_{\mathrm{ij}}$ ) are calculated based on equations (13)–(15) as the temperature of each DPD particle changes at each time step. ¹⁹

Table 1.

Simulation parameters.

$L_x\times L_y\times L_z$	$50\times 10\times 10$
Density (n)	4
Time step ( $\Delta t$ )	0.002
$\overline{C_v}$	$1\times {10}^5$
$k_0$	0.000126
g	0.02
$a_{\mathrm{ijf}}$	18.75
$a_{\mathrm{ijw}}$	5.0
$a_{\mathrm{ij}}$	3.0
Total number of iterations	30,000
Grid number	$50\times 10$
Number of iteration for averaging	5000

Case I: same wall temperature

In this case, the temperatures at top and bottom walls are kept the same, that is, $T_C=0.8$ , and the initial temperature of the fluid confined between the top and bottom walls is set as $T_0=1.0$ . The schematic of the computational domain is shown in Figure 5. A body force $g=0.02$ is applied on the fluid particles. Time-evolution of the temperature profile is shown in Figure 6. As shown in Figure 7, dimensionless temperature distribution is developed to steady state in the fully developed temperature regime.

Figure 5.

Schematic of the solution domain.

Figure 6.

Time-evolution of the temperature profile.

Figure 7.

Dimensionless temperature profile.

Sofos et al. ³³ examined the effect of channel width and system temperature on diffusion coefficient, shear viscosity, and thermal conductivity in the MD simulations. They mentioned that the rise in temperature results in a significant decrease in viscosity. Therefore, the effect of viscosity and its variation cannot be ignored and is also considered in this work.

The kinematic viscosity of the eDPD fluid can be computed by ³⁴

$$\eta =\frac{315{\mathrm{mk}}_{B}T}{128\pi \gamma r_c^3}+\frac{512\pi \gamma n^2{r}_c^5}{51,975}$$

(20)

Using the above-mentioned parameters (i.e. m = 1,r_c  = 1, n = 4, T = 1.0, and γ = 4.5), the value of the kinematic viscosity of the eDPD fluid is 2.402. Furthermore, for the Poiseuille flow, the velocity distribution in the x direction across the channel can be calculated ¹⁸ by

$$V_x=\frac{{\mathrm{nfL}}_Z^2}{8\eta }[1-{\left(\frac{z}{L_Z/2}\right)}^2]$$

(21)

And the maximal velocity of $V_x$ is

$$V_{\mathrm{xmax}}=\frac{{\mathrm{nfL}}_Z^2}{8\eta }$$

(22)

As shown in Figure 8, the simulated maximal velocity for this case was 0.341, and the difference from the analytical solution calculated by equation (22) was 18%. The difference is larger, and it is mainly due to the variation in the viscosity with the variation in the temperature. In fact, in this case, the fluid temperature approaches the wall temperature infinitely for enough time steps. So, the wall temperature could be adopted as the approximations of the fully developed fluid temperature to calculate the maximal velocity by equation (22). The newly calculated value of the analytical solution is 0.345, and the difference from the simulated value of 0.341 is only about 1.2%. Figure 8 shows the comparison between the eDPD and the analytical solutions, and it is found that they are in good agreement with each other.

Figure 8.

Velocity profile.

The Nusselt number (Nu) is the ratio of convective to conductive heat transfer, and it is given by the following equation

$$\mathrm{Nu}=\frac{{\mathrm{hD}}_h}{k}$$

(23)

where $D_h$ is the hydraulic diameter, $k$ is the thermal conductivity of the fluid, and $h$ is the convective heat transfer coefficient of the fluid and can be calculated from the gradient of the temperature profile as follows

$$h={-\frac{k}{\overline{T}-T_w}·\frac{\partial T}{\partial z}|}_{z=z_0}$$

(24)

where the temperature gradient is given by

$$\frac{\partial T}{\partial z}=\frac{\partial (T-T_w)}{\partial z}·\frac{\overline{T}-T_w}{T-T_w}=\frac{\partial \Theta }{\partial z}(\overline{T}-T_w)$$

(25)

So, the convective heat transfer coefficient of the fluid is expressed as

$$h={-k·\frac{\partial \Theta }{\partial z}|}_{z=z_0}$$

(26)

Substituting equation (26) into equation (23), the Nusselt number (Nu) is expressed as

$$\mathrm{Nu}={D_h(-\frac{\partial \Theta }{\partial z})|}_{z=z_0}$$

(27)

Three layers of fluid particles close to the wall are chosen to calculate the Nusselt number. So, the variation of density near the wall has a great effect on the final results. As shown in Figure 9, particle distribution is relatively uniform when adopting the bounceback boundary condition mentioned in section “Initial and boundary conditions” in which $a_{\mathrm{ijfw}}=10$ . The calculated Nusselt numbers for z = 0 and z = 10 by equation (27) are 7.72 and 7.28, respectively, and the average value is 7.50. The analytical value of fully developed convective heat transfer in grooved tube has been given by Shah and London, ³⁵ that is, 7.54 in this case. The error between the presented calculated value and analytical value is within 0.53%.

Figure 9.

Density profile.

Different wall temperatures

The forced convection in which the temperatures at top and bottom walls are different was also simulated. The schematic of the computational domain is shown in Figure 10. A body force $g=0.02$ is also applied on the fluid confined between a hot wall of $T_H=1.0$ and a cold wall of $T_C=0.8$ . The initial temperature of the fluid confined between the top and bottom walls is set at $T_0=0.9$ . As shown in Figure 11, the temperature profile is developed to steady state (linear profile) for 10,000 time steps. The result agrees well with the theoretic analytic result. The development of velocity profiles is shown in Figure 12. Since the viscosity increases with the decrease in temperature, the peak of the velocity profile shifts to the side corresponding to the high temperature fluid. ²⁵

Figure 10.

Schematics of the solution domain.

Figure 11.

Time-evolution of the temperature profile.

Figure 12.

Velocity profile of forced laminar convective heat transfer.

As shown in Figure 13(a), density distribution is horizontal as similar to that simulated by the classic DPD method. However, the character of fluids such as water is that density distribution depends on the temperature distribution, that is, lower density for the hot fluid and higher density for the cold fluid. Then, we follow the method used by Li et al. ²⁵ to regulate the conservative coefficient $a_{\mathrm{ijf}}=18.7$ into $a_{\mathrm{ijf}}=18.75k_B(T_i+T_j)/2$ . The result after adjustment is shown in Figure 13(b). Similar to Li et al., ²⁵ it reveals the effect of the temperature on the fluid density. We will adopt this adjustment in the later simulation.

Figure 13.

Density profile. (a) a _ijf  = 18.7 and (b) a _ijf  = 18.75k_B (T_i +T_j )/2.

Forced convection heat transfer in grooved microchannel

In this section, the problem of forced convection heat transfer on rough surface is studied using the eDPD approach and compared with FLUENT solution. The fluid in this microchannel is incompressible, Newtonian fluid, as shown in Figure 14. The temperatures of the hot and cold walls are maintained at a constant 1.0 and 1.5, respectively. The initial temperature of the eDPD fluid particles is set at 1.0. Other parameters are the same as before. Introducing the eDPD parameters into the predictions, equation (20) gives $\eta =2.402$ ; hence, approximately $\mathrm{Re}=100$ .

Figure 14.

Schematic of the solution domain.

Comparison between present eDPD simulation and FLUENT solution for the temperature isotherm is shown in Figure 15(a) and (b). Figure 16 shows the streamline of flow around the single obstacle in the channel. There is a clockwise vortex near the front of the lower left corner, and a clockwise recirculation is also developed downstream from the obstacle. It is found that the eDPD simulated results are in agreement with FLUENT data. These results demonstrate that the eDPD is effective in simulating this kind of flow and heat transfer. However, the size of the vortexes between the eDPD and FLUENT results is slightly different. We guess it may be caused by the boundary condition mentioned in section “Initial and boundary conditions,” and we expect it will be our future work.

Figure 15.

Temperature isotherms of flow around the single obstacle in the channel (a) FLUENT and (b) eDPD.

Figure 16.

Streamline of flow around the single obstacle in the channel (a) FLUENT and (b) eDPD.

Conclusion

In this article, forced convection heat transfer problems in microchannel were simulated numerically by eDPD. This concludes as follows.

The 2D unsteady heat conduction problem was simulated by eDPD. The results agree well with FLUENT results.

For the simulation of the forced convection on smooth surface that the temperatures of the upper wall and lower wall are different, the distribution of temperature is linear, which is in a good agreement with the theoretical solution. Since the viscosity changes with the temperature, the velocity distribution is asymmetric. The peak of the velocity profile shifts to the side corresponding to the high temperature fluid. In the situation that the upper wall and lower wall temperatures are same, in order to calculate the viscosity correctly, the variation in the fluid temperature should be taken into account. Furthermore, the simulated maximum velocity agrees well with the analytical value.

Finally, we give a qualitative analysis about the applicability of the eDPD approach in simulating flow and heat transfer in rough microchannel, and a good agreement is observed with each other. Above all, the eDPD model is a powerful tool in simulate flow and heat transfer problems. In the future work, we will extend the application of eDPD to simulate flow and heat transfer in more complex structures.