**2.4 Treatment of pressure in incompressible Navier-Stokes equations**

The spatial oscillations occur when central differencing is applied to both the continuity equation and the pressure gradient term in the momentum equations. The momentum equations at the even-numbered nodes depend only on pressures at odd-numbered nodes, and vice versa. The same holds for the continuity equation. This situation permits two different pressure fields to co-exist, which is known as checkerboard pressure field.

Nowadays, the staggered grid arrangement is no more necessary. The remarkable turn toward the colocated grid arrangement was the study presented by Rhie and Chow [22]. In the colocated grid arrangement all variables use the same grid and interpolation is needed. As it was mentioned, the colocated grid arrangement causes problems when central

Fig. 2. Ghost cells around the domain.

Fig. 3. The notations of the ghost cells.

**2.4 Treatment of pressure in incompressible Navier-Stokes equations** 

The spatial oscillations occur when central differencing is applied to both the continuity equation and the pressure gradient term in the momentum equations. The momentum equations at the even-numbered nodes depend only on pressures at odd-numbered nodes, and vice versa. The same holds for the continuity equation. This situation permits two

Nowadays, the staggered grid arrangement is no more necessary. The remarkable turn toward the colocated grid arrangement was the study presented by Rhie and Chow [22]. In the colocated grid arrangement all variables use the same grid and interpolation is needed. As it was mentioned, the colocated grid arrangement causes problems when central

different pressure fields to co-exist, which is known as checkerboard pressure field.

differencing is applied to cell-face velocities, Rhie and Chow [21] presented a method for avoiding the usage of the staggered grid arrangement. In this method central differencing has been applied to the pressure gradient and cell-face pressure, while the Rhie & Chow interpolation has been applied to the cell-face velocity as follow,

$$\begin{split} \mu\_{i-1/2,j,k} &= \frac{1}{2} (\mu\_{i-1,j,k} - \mu\_{i,j,k}) - \frac{1}{2} \left( \frac{V\_{i-1,j,k}}{a\_{i-1,j,k}} + \frac{V\_{i,j,k}}{a\_{i,j,k}} \right) \\ &\times \frac{\partial p}{\partial \mathbf{x}} \Big|\_{i-1/2,j,k} + \frac{1}{2} \left( \frac{V\_{i-1,j,k}}{a\_{i-1,j,k}} \frac{\partial p}{\partial \mathbf{x}} \Big|\_{i-1,j,k} + \frac{V\_{i,j,k}}{a\_{i,j,k}} \frac{\partial p}{\partial \mathbf{x}} \Big|\_{i,j,k} \right) \end{split} \tag{27}$$

#### **2.5 Pressure correction equation and multigrid technique**

The current research indicates that the AC-MG acceleration technique is highly efficient, reliable and robust, which makes it feasible for CPU-intensive computations, such as pressure Correction equations. When compared to the discretized momentum equations, the pressure Poisson equations tend to be very stiff and ill-conditioned, i.e. *ap≡*∑*nb anb* because of these reasons, solving the pressure Poisson equation is usually the CPU bottle-neck for the incompressible N–S equation system and AC-MG technique is required. With this acceleration technique the residuals of the large-scale algebraic equation system are guaranteed to be continuously driven down to the level of the computer machine round-off error and warrants strong conservations of mass and momentum satisfied over all the control volumes. In this cell centered multigrid algorithm both restriction and prolongation operators are based on piecewise constant interpolation.

Fig. 4. Schematic of a cell-centered two-level multigrid configuration.

Microchannel Simulation 199

, , ,2 2,2 2 ,2 2,2 1

*<sup>c</sup> f f ijk ij k ij k f f ij k ij k*

*r r*

, , , , 1, , , , , ,

*a a*

*a a*

*f f ff f ijk pijk eijk ijk i jk*

*ra a*

*f f ff nijk ij k tijk ijk f f ff wijk i jk sijk ij k*

ˆ ˆ

, , , 1, , , ,, 1

 

 

 

ˆ ˆ

1, , , 1, , , , ,

ˆ ˆ

,, 1 ,, , ,

*ff f bijk ijk ijk*

A typical two-level multigrid iterative algorithm consists of restriction, relaxation on the coarse grid and prolongation. After a number of relaxation sweeps, such as TDMA sweeps, on the fine grid level, the residuals are calculated using Equation (37) and are restricted to the coarse grid using Equations (29-36). The restricted residuals are then used as the source terms in Equation (28) and relaxation sweeps are used to solve Equation (28) on the coarse level. The solutions to Equation (28) are then utilized as the corrections to be prolongated back to the fine grid using the following relations to update the current iterative solution of

,2 2,2 2 ,2 2,2 2 , ,

*f f c ij k ij k ijk f f c ij k ij k ijk f f c ij k ij k ijk f f c ij k ij k ijk*

ˆ ˆ

ˆ ˆ

ˆ ˆ

ˆ ˆ

,2 2,2 1 ,2 2,2 1 , ,

,2 1,2 2 ,2 1,2 2 , ,

,2 1,2 1 ,2 1,2 1 , ,

Equation (38) provides the prolongation formulation in the AC multigrid calculation. Obviously, the restriction (Equations (29-36)) and prolongation (Equation (38)) are secondorder accuracy in space and no extra interpolations are needed. The AC-MG solution procedure contains two parts. The first part consists of a subroutine of TDMA sweeps, which is controlled by the residual convergence rate. A flowchart of the TDMA solver is depicted in Figure 5. If the total residual before iteration *n* on the current mesh is *Rn=*∑*i,j,k rni,j,k* and the residual after iteration *n* is *R n+1*, then another TDMA sweeping iteration is performed if the residual convergence rate satisfies *R n+1* ≥*f. R n* where the value for *f* is usually set to 0.5 [23]. If the convergence rate is lower than 0.5, i.e. *R n+1/R n* > *f,* a correction on the coarser grid is required, which invokes the second part of the AC-MG solution procedure. A schematic figure of the three-level AC-MG solution procedure will be shown

*a Sc*

*Sc r r*

ˆ

the residuals on the fine grid level *rf*

using the following relation:

, , ˆ *f i j k* 

, , ˆ*f i j k* ,

in Figure 6.

,2 1,2 2 ,2 1,2 1

(36)

(37)

(38)

*i,j,k* are calculated from the current iterative values of

The additive-correction multigrid scheme described in [23] is used for the two crossstreamwise directions (*y* and *z*). The cell-centred two-level multigrid configuration is sketched in Figure 4, in which the grid point (*i j k*) on the coarse level is surrounded by four grid points on the fine level in the *y* and *z* directions, namely (*i 2j-2 2k-2*), (*i 2j-2 2k-1*), (*i 2j-1 2k-2*) and (*i 2j-1 2k-1*). The variables on the coarse and fine levels are denoted by superscripts *c* and *f*, respectively, in the following equations. Based on the idea introduced in [23], the following algebraic equation system can be used to determine the correction for the fine grid level:

$$\begin{split} a\_{p(\{i,j,k\}}^{\mathcal{C}} \phi\_{l,j,k}^{\mathcal{C}} &= a\_{\mathfrak{e}(i,j,k)}^{\mathcal{C}} \phi\_{l+1,j,k}^{\mathcal{C}} + a\_{w(\{i,j,k\}}^{\mathcal{C}} \phi\_{l-1,j,k}^{\mathcal{C}} \\ &+ a\_{n(\{i,j,k\}}^{\mathcal{C}} \phi\_{l,j+1,k}^{\mathcal{C}} + a\_{s(\{i,j,k\}}^{\mathcal{C}} \phi\_{l,j-1,k}^{\mathcal{C}} \\ &+ a\_{t(\{i,j,k\}}^{\mathcal{C}} \phi\_{l,j,k+1}^{\mathcal{C}} + a\_{b(\{i,j,k\}}^{\mathcal{C}} \phi\_{l,j,k-1}^{\mathcal{C}} + Sc\_{i,j,k}^{\mathcal{C}} \end{split} \tag{28}$$

where the coefficients on the coarse grid are calculated from the following restriction formulae:

$$\begin{aligned} a\_{\epsilon(i,j,k)}^{\epsilon} &= a\_{\epsilon(i,2j-2,2k-2)}^{f} + a\_{\epsilon(i,2j-2,2k-1)}^{f} \\ &+ a\_{\epsilon(i,2j-1,2k-2)}^{f} + a\_{\epsilon(i,2j-1,2k-1)}^{f} \end{aligned} \tag{29}$$

$$\begin{aligned} a\_{w(i,j,k)}^c &= a\_{w(i,2j-2,2k-2)}^f + a\_{w(i,2j-2,2k-1)}^f \\ &+ a\_{w(i,2j-1,2k-2)}^f + a\_{w(i,2j-1,2k-1)}^f \end{aligned} \tag{30}$$

$$a\_{n\{i,j,k\}}^{c} = a\_{n\{i,2j-1,2k-2\}}^{f} + a\_{n\{i,2j-1,2k-1\}}^{f} \tag{31}$$

$$a\_{\mathbf{s}(i,j,k)}^{c} = a\_{\mathbf{s}(i,2j-2,2k-2)}^{f} + a\_{\mathbf{s}(i,2j-2,2k-1)}^{f} \tag{32}$$

$$a\_{t(i,j,k)}^c = a\_{t(i,2j-2,2k-1)}^f + a\_{t(i,2j-1,2k-1)}^f \tag{33}$$

$$a\_{b\left(i,j,k\right)}^{c} = a\_{b\left(i,2,j-2,2k-2\right)}^{f} + a\_{b\left(i,2,j-1,2k-2\right)}^{f} \tag{34}$$

$$\begin{aligned} a\_p^f(i,j,k) &= a\_{p(i,2j-2,2k-2)}^f + a\_{p(i,2j-2,2k-1)}^f \\ &+ a\_{p(i,2j-1,2k-2)}^f + a\_{p(i,2j-1,2k-1)}^f \\ &- a\_{s(i,2j-1,2k-2)}^f + a\_{s(i,2j-1,2k-1)}^f \\ &- a\_{n(i,2j-2,2k-2)}^f + a\_{n(i,2j-2,2k-1)}^f \\ &- a\_{b(i,2j-2,2k-1)}^f + a\_{b(i,2j-1,2k-1)}^f \\ &- a\_{t(i,2j-2,2k-2)}^f + a\_{t(i,2j-1,2k-2)}^f \\ &- a\_{t(i,2j-2,2k-2)}^f + a\_{t(i,2j-1,2k-2)}^f \end{aligned} \tag{35}$$

The additive-correction multigrid scheme described in [23] is used for the two crossstreamwise directions (*y* and *z*). The cell-centred two-level multigrid configuration is sketched in Figure 4, in which the grid point (*i j k*) on the coarse level is surrounded by four grid points on the fine level in the *y* and *z* directions, namely (*i 2j-2 2k-2*), (*i 2j-2 2k-1*), (*i 2j-1 2k-2*) and (*i 2j-1 2k-1*). The variables on the coarse and fine levels are denoted by superscripts *c* and *f*, respectively, in the following equations. Based on the idea introduced in [23], the following algebraic equation system can be used to determine the correction for the fine grid level:

*a a*

*c f f eijk ei j k ei j k f f ei j k ei j k*

*c f f*

*c f f*

*c f f*

*c f f*

*c f f*

*c f f pijk pi j k pi j k f f pi j k pi j k*

*aa a*

*aa a*

*aa a*

*aa a*

*c ccc c c pijk eijk ijk i jk wijk i jk c c cc nijk ij k sijk ij k cc cc c tijk ijk bijk ijk ijk*

, , , 1, , , , 1,

*a a Sc*

, , ,, 1 , , ,, 1 ,,

,2 1,2 2 ,2 1,2 1

,2 1,2 2 ,2 1,2 1

*nijk ni j k ni j k aa a* (31)

*sijk si j k si j k aa a* (32)

*tijk ti j k ti j k aa a* (33)

*bijk bi j k bi j k aa a* (34)

 

   

(28)

(29)

(30)

(35)

, , , , , , 1, , , , 1, ,

where the coefficients on the coarse grid are calculated from the following restriction formulae:

*a a*

*wijk wi j k wi j k f f wi j k wi j k*

*a a*

, , ,2 2,2 2 ,2 2,2 1

, , ,2 1,2 2 ,2 1,2 1

, , ,2 2,2 2 ,2 2,2 1

, , ,2 2,2 1 ,2 1,2 1

, , ,2 2,2 2 ,2 1,2 2

*f f si j k si j k f f ni j k ni j k*

*a a*

*a a*

*a a*

*a a*

*a a*

*f f bi j k bi j k*

*f f ti j k ti j k*

, , ,2 2,2 2 ,2 2,2 1

,2 1,2 2 ,2 1,2 1

,2 1,2 2 ,2 1,2 1

,2 2,2 2 ,2 2,2 1

,2 2,2 1 ,2 1,2 1

,2 2,2 2 ,2 1,2 2

, , ,2 2,2 2 ,2 2,2 1

$$\begin{aligned} Sc\_{i,j,k}^{\mathbb{C}} &= r\_{i,2j-2,2k-2}^{f} + r\_{i,2j-2,2k-1}^{f} \\ &+ r\_{i,2j-1,2k-2}^{f} + r\_{i,2j-1,2k-1}^{f} \end{aligned} \tag{36}$$

the residuals on the fine grid level *rf i,j,k* are calculated from the current iterative values of , , ˆ *f i j k* using the following relation:

$$\begin{split} r^{f}\_{i,j,k} &= -a^{f}\_{p(i,j,k)} \hat{\phi}^{f}\_{i,j,k} + a^{f}\_{e(i,j,k)} \hat{\phi}^{f}\_{i+1,j,k} \\ &+ a^{f}\_{\text{in}\{i,j,k\}} \hat{\phi}^{f}\_{i,j+1,k} + a^{f}\_{t(i,j,k)} \hat{\phi}^{f}\_{i,j,k+1} \\ &+ a^{f}\_{w(i,j,k)} \hat{\phi}^{f}\_{i-1,j,k} + a^{f}\_{s(i,j,k)} \hat{\phi}^{f}\_{i,j-1,k} \\ &+ a^{f}\_{b(i,j,k)} \hat{\phi}^{f}\_{i,j,k-1} + \text{Sc}^{f}\_{i,j,k} \end{split} \tag{37}$$

A typical two-level multigrid iterative algorithm consists of restriction, relaxation on the coarse grid and prolongation. After a number of relaxation sweeps, such as TDMA sweeps, on the fine grid level, the residuals are calculated using Equation (37) and are restricted to the coarse grid using Equations (29-36). The restricted residuals are then used as the source terms in Equation (28) and relaxation sweeps are used to solve Equation (28) on the coarse level. The solutions to Equation (28) are then utilized as the corrections to be prolongated back to the fine grid using the following relations to update the current iterative solution of , , ˆ*f i j k* ,

$$\begin{aligned} \hat{\phi}\_{i,2j-2,2k-2}^{f} &= \hat{\phi}\_{i,2j-2,2k-2}^{f} + \phi\_{i,j,k}^{c} \\ \hat{\phi}\_{i,2j-2,2k-1}^{f} &= \hat{\phi}\_{i,2j-2,2k-1}^{f} + \phi\_{i,j,k}^{c} \\ \hat{\phi}\_{i,2j-1,2k-2}^{f} &= \hat{\phi}\_{i,2j-1,2k-2}^{f} + \phi\_{i,j,k}^{c} \\ \hat{\phi}\_{i,2j-1,2k-1}^{f} &= \hat{\phi}\_{i,2j-1,2k-1}^{f} + \phi\_{i,j,k}^{c} \end{aligned} \tag{38}$$

Equation (38) provides the prolongation formulation in the AC multigrid calculation. Obviously, the restriction (Equations (29-36)) and prolongation (Equation (38)) are secondorder accuracy in space and no extra interpolations are needed. The AC-MG solution procedure contains two parts. The first part consists of a subroutine of TDMA sweeps, which is controlled by the residual convergence rate. A flowchart of the TDMA solver is depicted in Figure 5. If the total residual before iteration *n* on the current mesh is *Rn=*∑*i,j,k rni,j,k* and the residual after iteration *n* is *R n+1*, then another TDMA sweeping iteration is performed if the residual convergence rate satisfies *R n+1* ≥*f. R n* where the value for *f* is usually set to 0.5 [23]. If the convergence rate is lower than 0.5, i.e. *R n+1/R n* > *f,* a correction on the coarser grid is required, which invokes the second part of the AC-MG solution procedure. A schematic figure of the three-level AC-MG solution procedure will be shown in Figure 6.

Microchannel Simulation 201

The total grid number is 73,800 (*Nx×Ny×Nz* is 30 × 82 ×30) for the domain. This type of a fine grid mesh for the *y* and *z* directions was chosen in order to properly resolve the velocity and viscous shear layers, and to more accurately define the conjugate heat transfer at the surface of the channel, thereby improving the temperature resolution. Furthermore, comparison with standard theoretical or numerical results indicates that the finer the mesh size the higher the numerical accuracy. The reasons for the comparative coarse discretization for the *x*-direction are: (*i*) with the exception of the inlet region, the temperature gradients are small compared to the gradients occurring in other directions; and (*ii*) The CPU time as well as the memory storage required increases dramatically as the number of grid nodes is increased.

The velocity field can be determined analytically using a more direct approach. As illustrated in Figure 7, the following relations represent the reasonable results for this type

> 

cosh / 2 1 cos cosh / 2 2

 

y

*n*

3 3 1,3,... <sup>16</sup> <sup>1</sup> <sup>1</sup>

2

*c a*

*u*

2 1

max

*u*

Fig. 7. A schematic of the rectangular channel.

where *m* and *n* are derived from below relations,

following form for the aspect ratio \*

*b a n*

*n*

(40)

(41)

(39)

1 1 /2

*ny a n z nb a a*

*n*

5 5 1,3,...

 

This velocity profile is in excellent agreement with the experimental results. Since Equation (39) involves considerable computational complexity, a simple approximation in the

2a

2b z

0.5 is proposed [24],

1 1 *<sup>n</sup> <sup>m</sup> u z y u ba*

<sup>192</sup> <sup>1</sup> 1 tanh 3 2 *<sup>m</sup> n*

*ca a n b*

**3. Validation of the code** 

of problem by Shah and London [24],

**3.1 Velocity field** 

Fig. 5. Flowchart of the TDMA solver controlled by the residual convergence rate.

In order to check the sensitivity of the numerical results to mesh size, three different grid systems were tested. They consisted of 18×42×16, 30×82×30 and 50×162×58 nodes in the *x*, *y*, and *z* directions, respectively. The results from the last two grid systems were very close to each other and local temperature differences were less than 0.1%. Since less computational time and computer memory were needed for the second grid system, it was employed in the final simulation.

Fig. 6. A schematic of the three-level AC-MG solution procedure.

Calculate the total residual: *ressum1* 

Perform sweeps in three spatial directions

Calculate the total residual: *ressum* 

*ressum* < *ε0 ?* Is

Yes

Fig. 5. Flowchart of the TDMA solver controlled by the residual convergence rate.

Return

Is *decrat* < *ε1*?

Calculate the convergence rate, *decrat = ressum/ressum1* 

No

No

Yes

Fig. 6. A schematic of the three-level AC-MG solution procedure.

employed in the final simulation.

In order to check the sensitivity of the numerical results to mesh size, three different grid systems were tested. They consisted of 18×42×16, 30×82×30 and 50×162×58 nodes in the *x*, *y*, and *z* directions, respectively. The results from the last two grid systems were very close to each other and local temperature differences were less than 0.1%. Since less computational time and computer memory were needed for the second grid system, it was

*ressum1=ressum*

The total grid number is 73,800 (*Nx×Ny×Nz* is 30 × 82 ×30) for the domain. This type of a fine grid mesh for the *y* and *z* directions was chosen in order to properly resolve the velocity and viscous shear layers, and to more accurately define the conjugate heat transfer at the surface of the channel, thereby improving the temperature resolution. Furthermore, comparison with standard theoretical or numerical results indicates that the finer the mesh size the higher the numerical accuracy. The reasons for the comparative coarse discretization for the *x*-direction are: (*i*) with the exception of the inlet region, the temperature gradients are small compared to the gradients occurring in other directions; and (*ii*) The CPU time as well as the memory storage required increases dramatically as the number of grid nodes is increased.

### **3. Validation of the code**

#### **3.1 Velocity field**

The velocity field can be determined analytically using a more direct approach. As illustrated in Figure 7, the following relations represent the reasonable results for this type of problem by Shah and London [24],

$$\begin{aligned} u &= -\frac{16c\_1a^2}{\pi^3} \sum\_{n=1,3,\dots}^{\infty} \frac{1}{n^3} (-1)^{(n-1)/2} \\ &\times \left[ 1 - \frac{\cosh\left(n\pi y / / 2a\right)}{\cosh\left(n\pi b / 2a\right)} \right] \cos\left(\frac{n\pi z}{2a}\right) \end{aligned} \tag{39}$$

$$u\_m = -\frac{c\_1a^2}{3} \left[ 1 - \frac{192}{\pi^5} \left( \frac{a}{b} \right) \sum\_{n=1,3,\dots}^{\infty} \frac{1}{n^5} \tanh\left(\frac{n\pi b}{2a}\right) \right] \tag{40}$$

Fig. 7. A schematic of the rectangular channel.

This velocity profile is in excellent agreement with the experimental results. Since Equation (39) involves considerable computational complexity, a simple approximation in the following form for the aspect ratio \* 0.5 is proposed [24],

$$\frac{\mu}{\mu\_{\text{max}}} = \left[1 - \left(\frac{y}{b}\right)^{n}\right] \left[1 - \left(\frac{z}{a}\right)^{m}\right] \tag{41}$$

where *m* and *n* are derived from below relations,

Microchannel Simulation 203

Fig. 8. Velocity field in channel from the approximate analytical expression Equation (44),

Fig. 9. Velocity field in channel from the numerical calculation, *Δp* =50 *kPa*, *Re*=92.383,

Figure 10 compares the analytical friction coefficient as determined from Equation (46) with the numerical results obtained using the following procedure, should be noted that, *fDarcy=4f*.

*<sup>f</sup>* . constant, , *ui jk y z* (49)

*Δp*=50 *kPa*, *Re*=92.68, *Treference*=20o*C*, *um*=1.0779 *m/s*.

*Treference*=20o*C*, *um*=1.1032*m/s*, *umax*=1.997 *m/s*.

and then, the mean velocity is obtained as

The mass flow rate is calculated from the velocity as,

*<sup>m</sup>*

$$m = 1.7 + 0.5 \left( a^\* \right)^{-1.4} \tag{42}$$

$$m = \begin{cases} 2 & \text{for } a^\* \le 1/3 \\ 2 + 0.3\left(a^\* - 1/3\right) & \text{for } a^\* \ge 1/3 \end{cases} \tag{43}$$

The integration of Equation (41) over the duct cross section yields,

$$\frac{u}{u\_{\rm m}} = \left(\frac{m+1}{m}\right) \left(\frac{n+1}{n}\right) \left[1 - \left(\frac{y}{b}\right)^{n}\right] \left[1 - \left(\frac{z}{a}\right)^{m}\right] \tag{44}$$

$$\frac{u\_{\text{max}}}{u\_{\text{m}}} = \left(\frac{m+1}{m}\right) \left(\frac{n+1}{n}\right) \tag{45}$$

With *fRe* of Equation (46), *um* for the rectangular ducts can be expressed in a closed form as,

$$\begin{aligned} f \operatorname{Re}\_{D\_h} &= 24(1 - 1.3553a^\ast + 1.9467a^\ast^2) \\ &- 1.7012a^{\ast 3} + 0.9564a^{\ast 4} - 0.2537a^{\ast 5}) \end{aligned} \tag{46}$$

$$\mu\_m = -\frac{8c\_1 a^2}{f \operatorname{Re}\_{D\_h} \left[ 1 + \left( a/b \right)^2 \right]} \tag{47}$$

where

$$c\_1 = \frac{dp/d\mathbf{x}}{\mu} = \frac{\Delta p/L\_\chi}{\mu} \tag{48}$$

Substituting Equations (46) and (48) into Equation (47), the mean velocity for a given pressure drop, *um* can be obtained. Then, using the resulting value for *um* and Equations (42- 44), the approximate analytical velocity distribution in the microchannel as shown in Figure 8 can be obtained. The numerically determined velocity profile developed here is illustrated in Figure 9. Comparison of the analytical and numerical results indicates that while the numerical code exactly represents the general trend of the results, there is some disparity between the analytical and numerical results. The small difference of the velocity profile between Figures 8 and 9 is due to the approximations in the analytical solution described in Equations (41) and (44). But clearly, as evidenced by the magnitude of the mean velocities and the Reynolds numbers obtained from the different methods, the agreement between the two methods is quite good and provides sufficient evidence for validation of the numerical method. In this comparison, the thermophysical properties of water were chosen at a temperature of 293 K (20oC). Because the thermophysical properties are temperature dependent, especially the liquid viscosity, the velocity and the Reynolds numbers are different under the same pressure drop conditions. This issue will be discussed in more detail later.

 1.4 \* *<sup>m</sup>* 1.7 0.5 

\* \* 2 1 3 2 0.3 1 3 1 3

 

1 1 1 1 *<sup>n</sup> <sup>m</sup> um n y z u mn b a* 

> *u m n* 1 1 *u mn*

With *fRe* of Equation (46), *um* for the rectangular ducts can be expressed in a closed form as,

Re 24(1 1.3553 1.9467

1.7012 0.9564 0.2537 )

8

*f ab* 

*c a*

Re 1 *<sup>h</sup>*

<sup>1</sup> *<sup>x</sup> dp dx p L <sup>c</sup>* 

Substituting Equations (46) and (48) into Equation (47), the mean velocity for a given pressure drop, *um* can be obtained. Then, using the resulting value for *um* and Equations (42- 44), the approximate analytical velocity distribution in the microchannel as shown in Figure 8 can be obtained. The numerically determined velocity profile developed here is illustrated in Figure 9. Comparison of the analytical and numerical results indicates that while the numerical code exactly represents the general trend of the results, there is some disparity between the analytical and numerical results. The small difference of the velocity profile between Figures 8 and 9 is due to the approximations in the analytical solution described in Equations (41) and (44). But clearly, as evidenced by the magnitude of the mean velocities and the Reynolds numbers obtained from the different methods, the agreement between the two methods is quite good and provides sufficient evidence for validation of the numerical method. In this comparison, the thermophysical properties of water were chosen at a temperature of 293 K (20oC). Because the thermophysical properties are temperature dependent, especially the liquid viscosity, the velocity and the Reynolds numbers are different under the same pressure drop conditions. This issue will be discussed in more

*D*

\*3 \*4 \*5

2 1

2

max m

*n*

m

The integration of Equation (41) over the duct cross section yields,

*Dh f*

where

detail later.

*m*

*u*

(42)

(45)

(48)

(43)

(44)

(47)

\*

\* \*2

(46)

*for*

*for*

 

Fig. 8. Velocity field in channel from the approximate analytical expression Equation (44), *Δp*=50 *kPa*, *Re*=92.68, *Treference*=20o*C*, *um*=1.0779 *m/s*.

Fig. 9. Velocity field in channel from the numerical calculation, *Δp* =50 *kPa*, *Re*=92.383, *Treference*=20o*C*, *um*=1.1032*m/s*, *umax*=1.997 *m/s*.

Figure 10 compares the analytical friction coefficient as determined from Equation (46) with the numerical results obtained using the following procedure, should be noted that, *fDarcy=4f*. The mass flow rate is calculated from the velocity as,

$$\dot{m} = \sum\_{} \sum\_{} \rho\_f \, \iota \iota \{i = \text{constant}, j, k\} \Delta y \Delta z \tag{49}$$

and then, the mean velocity is obtained as

Microchannel Simulation 205

shorter entrance lengths occurring for square-edged entrances than for rounded ones. Thus, because the flow entrance length may be less than 5% of the total length for a rectangular channel heat sink, the assumption of fully developed laminar flow over the entire length of the microchannel is acceptable for the heat transfer analysis, particularly in cases such as

The code was first validated for one dimensional heat conduction by comparing the results with a 1D analytical solution of heat conduction with a specified boundary condition [25]. The agreement was quite good and indicated very good correlation between the numerical results and the *1D* analytical solution, Secondly, using conservation of energy, it can be shown that the maximum possible temperature rise between the fluid inlet and outlet can be

, *s s f con*

*q A <sup>T</sup> m c*

In addition, the average temperature rise between the inlet and outlet of the channel can be

.. ..

.. ..

.

*uc T y z*

*m c*

*uc T y z*

*p*

.

*m c*

*p*

*f ave f x f*

*T TxL Tx*

In this work, three different cases (*qs*=90 *W/cm2*, *Δp*=50, 15 and 6 *kPa*) were investigated. Comparison of the results in Table 2, indicates that the difference between *ΔTf,ave* and *ΔTf,con* is small. This issue also is illustrated in Figure 11. Differences of this magnitude can be attributed to (*i*) Equation (54) is the maximum possible temperature rise in the bulk liquid from the energy balance; (*ii*) the mesh size is not as fine as required (infinitesimal), hence the

*Δp(kPa)* 50 15 6

*ΔTf,con (°C)* 14.62 36.82 78.57 *ΔTf,ave (°C)* 12.79 35.10 77.16

Table 2. Comparison between *ΔTf,ave* and *ΔTf,con* for *Δp*=50, 15 and 6 *kPa, qw*=90 *W/cm2.* 

*Re* 162.68 85.60 47.32

*p*

0

*f p i outlet j k*

*f p i inlet j k*

(54)

, ,

(55)

, ,

this where the Reynolds number is less than 200 (or low mass flow rates).

determined from the numerical analysis as follow,

,

accuracy of the statistical result from Equation (55) is limited.

**3.2 Heat transfer** 

expressed as:

$$
\mu\_m = \frac{\dot{m}}{\rho\_f A} = \frac{\dot{m}}{\rho\_f H W} \tag{50}
$$

then using Equation (51), the friction factor can be determined as Equation (52),

$$
\Delta p = f \frac{4L\_\chi}{D\_h} \frac{\rho\_f u\_m^2}{2} \tag{51}
$$

$$f\_{\rm Dark} \,\mathrm{Re}\_{D\_h} = 2 \frac{\Delta p D\_h^2}{\mu\_m L\_\chi \mu\_f} \tag{52}$$

In References [25-27], the friction coefficient, *fDarcyRe*, is determined numerically for different duct cross-sections. For the rectangular channel with an aspect ratio (height to width) of 3–1 (*H:W*), which approximates the geometry used here, 180µm×57µm, the Darcy friction factor– Reynolds number product, *fDarcyRe*, is 69. The agreement between the numerical calculations here and the calculations obtained by others [25-27] represents that the numerical code developed here is quite accurate for the fully developed laminar flow.

Fig. 10. Comparison among the numerical calculations, the analytical and the experimental data for the friction coefficient.

The length required for the formation of a fully developed laminar profile in a microchannel can be estimated by the following analytical relation [26] that is developed for a round tube,

$$\frac{L\_{\varepsilon}}{D\_{h}} = 0.0575 \,\mathrm{Re}\_{D\_{h}}\tag{53}$$

For a hydraulic diameter of *Dh* = 86.58 µm and for *Re* = 160 the entrance length is 796.5 µm. In References [25–27], it is noted that the shape of the entrance is very important, with much

*f f m m*

 

<sup>2</sup> <sup>4</sup> 2 *x f m*

*L u*

2

*h*

*mx f*

*h*

Re 2 *<sup>h</sup>*

*pD <sup>f</sup> u L*

In References [25-27], the friction coefficient, *fDarcyRe*, is determined numerically for different duct cross-sections. For the rectangular channel with an aspect ratio (height to width) of 3–1 (*H:W*), which approximates the geometry used here, 180µm×57µm, the Darcy friction factor– Reynolds number product, *fDarcyRe*, is 69. The agreement between the numerical calculations here and the calculations obtained by others [25-27] represents that the numerical code

Fig. 10. Comparison among the numerical calculations, the analytical and the experimental

The length required for the formation of a fully developed laminar profile in a microchannel can be estimated by the following analytical relation [26] that is developed for a round tube,

For a hydraulic diameter of *Dh* = 86.58 µm and for *Re* = 160 the entrance length is 796.5 µm. In References [25–27], it is noted that the shape of the entrance is very important, with much

*h L*

0.0575Re *<sup>h</sup> <sup>e</sup> <sup>D</sup>*

*<sup>D</sup>* (53)

*A HW* (50)

(51)

(52)

*m*

then using Equation (51), the friction factor can be determined as Equation (52),

*Darcy D*

developed here is quite accurate for the fully developed laminar flow.

data for the friction coefficient.

*p f <sup>D</sup>*

*u*

shorter entrance lengths occurring for square-edged entrances than for rounded ones. Thus, because the flow entrance length may be less than 5% of the total length for a rectangular channel heat sink, the assumption of fully developed laminar flow over the entire length of the microchannel is acceptable for the heat transfer analysis, particularly in cases such as this where the Reynolds number is less than 200 (or low mass flow rates).

#### **3.2 Heat transfer**

The code was first validated for one dimensional heat conduction by comparing the results with a 1D analytical solution of heat conduction with a specified boundary condition [25]. The agreement was quite good and indicated very good correlation between the numerical results and the *1D* analytical solution, Secondly, using conservation of energy, it can be shown that the maximum possible temperature rise between the fluid inlet and outlet can be expressed as:

$$
\Delta T\_{f,con} = \frac{q\_s \cdot A\_s}{\dot{m} \cdot \varepsilon\_p} \tag{54}
$$

In addition, the average temperature rise between the inlet and outlet of the channel can be determined from the numerical analysis as follow,

$$\begin{aligned} \Delta T\_{f, \text{ave}} &= \overline{T}\_f \left( \mathbf{x} = \mathbf{L}\_x \right) - \overline{T}\_f \left( \mathbf{x} = \mathbf{0} \right) \\ &= \frac{\left\{ \sum \sum \rho\_f \, \mathbf{u} \, \mathbf{c}\_p, T \, \Delta \mathbf{y} \, \Delta z \right\}\_{i=\text{outlet}, j,\mathbf{k}}}{\dot{\mathbf{m}} \, \mathbf{c}\_p} \\ &- \frac{\left\{ \sum \sum \rho\_f \, \mathbf{u} \, \mathbf{c}\_p, T \, \Delta \mathbf{y} \, \Delta z \right\}\_{i=\text{inlet}, j,\mathbf{k}}}{\dot{\mathbf{m}} \, \mathbf{c}\_p} \end{aligned} \tag{55}$$

In this work, three different cases (*qs*=90 *W/cm2*, *Δp*=50, 15 and 6 *kPa*) were investigated. Comparison of the results in Table 2, indicates that the difference between *ΔTf,ave* and *ΔTf,con* is small. This issue also is illustrated in Figure 11. Differences of this magnitude can be attributed to (*i*) Equation (54) is the maximum possible temperature rise in the bulk liquid from the energy balance; (*ii*) the mesh size is not as fine as required (infinitesimal), hence the accuracy of the statistical result from Equation (55) is limited.


Table 2. Comparison between *ΔTf,ave* and *ΔTf,con* for *Δp*=50, 15 and 6 *kPa, qw*=90 *W/cm2.* 

Microchannel Simulation 207

As shown, a variation in the reference temperature, *Treference* from 20 to 32°C, changes the mean velocity from 1.1032 to 1.44 *m/s*, and results in a corresponding change in the Reynolds number from 95.38 to 162.68. The numerical results for the temperature distribution in the heat sinks are shown in Figures 13, 14, 15, 16 and 17 for different locations along the channel. Figures 13, 14 and 15 show the local cross-sectional temperature distribution in the *y*–*z* plane at *x=0*, *x=Lx/2* and *x=Lx*, respectively. As shown in Figure 13, the temperature of the liquid at the inlet is initially uniform (at 20°C). The temperature profiles shown in Figures 14 and 15 are identical in shape due to the assumption of hydrodynamic fully developed flow, but the magnitudes of the

Fig. 13. Local temperature distribution in *y–z* plane at *x=0*, (*Δp*=50 *kPa*, *Re*=162.68,

Fig. 14. Local temperature distribution in *y–z* plane at *x=Lx/2*, (*Δp*=50 *kPa*, *Re*=162.68,

temperature are different.

*Treference*=32°C, *um*=1.44 *m/s*).

*Treference*=32°C, *um*=1.44 *m/s*).

Fig. 11. Comparison between the numerical and analytical results for temperature differences upon the Reynolds number.
