**2.1 Problem description**

The micro-heat sink modeled in this investigation consists of a 10 mm long silicon substrate with a silicon cover. The rectangular microchannels have a width of 57 µm and a depth of 180 µm. The hydraulic diameter of microchannel under development is about 86.58 µm and is expected to decrease to 10 µm. This yields a typical Knudsen number for water as a coolant to be between 3.5×10-5 and 3.5×10-4 which lies in the continuous flow regime (*Kn*<10-3) [20]. Hence, the conservation equations based on the continuum model (Navier-Stokes equations of motion) can still be used to describe the transport processes.

A schematic of the rectangular microchannel heat sink is shown in Figure 1 where a unit of cell consisting of one channel was selected because of the symmetry of the structure. The unit cell of the heat sink can be represented by an asymmetric rectangular channel with the cross-sectional dimensions as shown in Table 1. The channel geometry is similar to that employed in the experimental work of Kawano et al. [14] and in the numerical investigations of Qu and Mudawar [5] and Fedorov and Viskanta [6]. It is assumed that the heat flux input at the bottom of the heat sink is uniform.

#### **2.2 Model equations and boundary conditions**

Consider a steady 3D flow in a silicon microchannel heat sink with heating from below and with adiabatic conditions at the other boundaries, as showed in Figure 1. There are some assumptions in this numerical simulation, e.g. the transport processes are considered to be steady-state and three-dimensional, the flow is incompressible and laminar, thermal radiation is neglected, the thermophysical properties are temperature dependent. Under the stated assumptions, the governing equations and related boundary conditions for a fully developed 3D flow heat transfer are given as [21],

Continuity equation

$$\frac{\partial}{\partial \mathbf{x}}(\rho u) + \frac{\partial}{\partial y}(\rho v) + \frac{\partial}{\partial z}(\rho w) = 0 \tag{1}$$

Fig. 1. Part I

Microchannel Simulation 193

*c uT c vT c wT p pp xyz*

*uvw* 0, 0, 0 (6)

*<sup>s</sup>* (7)

, (1 ), 0, 0 *x f out x L p p atm v w* (8)

(9)

*y* 

> *y*

> > *f in*

*<sup>s</sup> <sup>s</sup>*

*<sup>s</sup> <sup>s</sup>*

 

 

0

0

0

(5)

(10)

(11)

(12)

(13)

 

*TTT kkk xx yy zz*

where *u, v, w, p, ρ, μ, T, k* and *cp* are the velocity in *x*-direction, *y*-direction, *z*-direction, pressure, density, dynamic viscosity, temperature, thermal conductivity and specific heat at constant pressure respectively. The hydrodynamic boundary conditions are as follow:

> 0, , 1 , 0, 0 *f in m*

> > *s s s*

*TTT kkk xx yy zz*

0 , 0,0 *<sup>s</sup> x z ss <sup>T</sup> xLy zL k q*

0 , ,0 <sup>0</sup> *<sup>s</sup> xy z s <sup>T</sup> xLy L zL k*

*L H yL <sup>x</sup> T T W zW W*

 / 2 / 2 <sup>0</sup>

*<sup>T</sup> else <sup>k</sup> <sup>x</sup>*

*yy f*

*L H yL T*

 

*W zW W x*

*<sup>T</sup> else <sup>k</sup> <sup>x</sup>*

/ 2 / 2 <sup>0</sup> /2 /2

/2 /2

*x L k*

*w w*

*x f*

*w w*

*y y*

*x pp u v w*

Energy equation

at the inner wall surface (no slip)

the heat conduction in the solid section is,

the thermal boundary conditions can be stated as:

at the inlet of channel

at the outlet of channel

Fig. 1. Schematic of a rectangle microchannel heat sink and the unit of cell.


*y xxyyzz*

Table 1. Geometric dimensions of the unit cell.

Momentum equations

$$\begin{aligned} &\frac{\partial}{\partial \mathbf{x}} \Big(\rho u^2\Big) + \frac{\partial}{\partial y} \Big(\rho uv\Big) + \frac{\partial}{\partial z} \Big(\rho uw\Big) \\ &= -\frac{\partial P}{\partial \mathbf{x}} + \frac{\partial}{\partial \mathbf{x}} \Big(\mu \frac{\partial u}{\partial \mathbf{x}}\Big) + \frac{\partial}{\partial y} \Big(\mu \frac{\partial u}{\partial y}\Big) + \frac{\partial}{\partial z} \Big(\mu \frac{\partial u}{\partial z}\Big) \\ &\qquad \frac{\partial}{\partial \mathbf{x}} (\rho uv) + \frac{\partial}{\partial y} \Big(\rho v^2\Big) + \frac{\partial}{\partial z} \Big(\rho ww\Big) \\ &= -\frac{\partial P}{\partial y} + \frac{\partial}{\partial \mathbf{x}} \Big(\mu \frac{\partial v}{\partial \mathbf{x}}\Big) + \frac{\partial}{\partial y} \Big(\mu \frac{\partial v}{\partial y}\Big) + \frac{\partial}{\partial z} \Big(\mu \frac{\partial v}{\partial z}\Big) \end{aligned} \tag{3}$$

$$\begin{aligned} &\frac{\partial}{\partial \mathbf{x}} (\rho wu) + \frac{\partial}{\partial y} (\rho wv) + \frac{\partial}{\partial z} (\rho w^2) \\ &= -\frac{\partial P}{\partial z} + \frac{\partial}{\partial x} \left( \mu \frac{\partial w}{\partial x} \right) + \frac{\partial}{\partial y} \left( \mu \frac{\partial w}{\partial y} \right) + \frac{\partial}{\partial z} \left( \mu \frac{\partial w}{\partial z} \right) \end{aligned} \tag{4}$$

(Part II)

*Ly m*

180 10 100 900 86.58 43 57

*Pu u u xx x y y z z*

*Pv v v y xxyyzz*

*Pw w w zx x y y z z*

 <sup>2</sup> *u uv uw xyz*

 <sup>2</sup> *uv v vw xyz*

 <sup>2</sup> *wu wv w xyz*

 

 

 

*Dh m* *Ww m*

*W m* 

(2)

(3)

(4)

Fig. 1. Schematic of a rectangle microchannel heat sink and the unit of cell.

*Lz m*

*H m* 

Momentum equations

Table 1. Geometric dimensions of the unit cell.

*Lx m*

#### Energy equation

$$\begin{aligned} &\frac{\partial}{\partial \mathbf{x}} \left( c\_p \rho uT \right) + \frac{\partial}{\partial y} \left( c\_p \rho vT \right) + \frac{\partial}{\partial z} \left( c\_p \rho wvT \right) \\ &= \frac{\partial}{\partial \mathbf{x}} \left( k \frac{\partial T}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial y} \left( k \frac{\partial T}{\partial y} \right) + \frac{\partial}{\partial z} \left( k \frac{\partial T}{\partial z} \right) \end{aligned} \tag{5}$$

where *u, v, w, p, ρ, μ, T, k* and *cp* are the velocity in *x*-direction, *y*-direction, *z*-direction, pressure, density, dynamic viscosity, temperature, thermal conductivity and specific heat at constant pressure respectively. The hydrodynamic boundary conditions are as follow:

at the inner wall surface (no slip)

$$u = 0, \quad v = 0, \quad w = 0 \tag{6}$$

at the inlet of channel

$$\mathbf{x} = \mathbf{0}, \quad p\_f = p\_{in\nu} \quad \mu = \mathbf{1} \frac{m}{s}, \quad \upsilon = \mathbf{0}, \quad w = \mathbf{0} \tag{7}$$

at the outlet of channel

$$\mathbf{x} = \mathbf{L}\_{\mathbf{x}'} \quad p\_f = p\_{out} \text{(1} at m), \quad \mathbf{v} = \mathbf{0}, \quad w = \mathbf{0} \tag{8}$$

the heat conduction in the solid section is,

$$\frac{\partial}{\partial \mathbf{x}} \left( k \frac{\partial T}{\partial \mathbf{x}} \right)\_s + \frac{\partial}{\partial y} \left( k \frac{\partial T}{\partial y} \right)\_s + \frac{\partial}{\partial z} \left( k \frac{\partial T}{\partial z} \right)\_s = \mathbf{0} \tag{9}$$

the thermal boundary conditions can be stated as:

$$0 \le x \le L\_{\ge \prime}, y = 0, 0 \le z \le L\_z \implies -k\_s \frac{\partial T\_s}{\partial y} = q\_s \tag{10}$$

$$0.0 \le x \le L\_{\infty}, y = L\_{y, \prime} \\ 0 \le z \le L\_z \implies \\ -k\_s \frac{\partial T\_s}{\partial y} = 0 \tag{11}$$

$$\mathbf{x} = \mathbf{0} \begin{cases} \left( L\_y \; / \; 2 - H \right) \le y \le L\_y \; / \; 2 \\ W\_w \; / \; 2 \le z \le \left( W\_w \; / \; 2 + W \right) \end{cases} \implies T\_f = T\_{in} \tag{12}$$
 
$$\mathbf{x} \in \mathbb{R}^{N \times N} \tag{12}$$

$$\text{else}\qquad\qquad\Rightarrow-k\_s\frac{\partial T\_s}{\partial \chi} = 0$$

$$\begin{aligned} \text{l.x} = \text{L}\_x \begin{cases} \left( \text{L}\_y \ne 2 - H \right) \le y \le \text{L}\_y \ne 2\\ \left( \text{W}\_w \ne 2 \le z \le \left( \text{W}\_w \ne 2 + \text{W} \right) \end{cases} \implies -\text{k}\_f \frac{\partial T\_f}{\partial \mathbf{x}} = \mathbf{0} \end{aligned} \tag{13}$$

$$\begin{aligned} \text{else} \quad \begin{aligned} \text{else} \quad \begin{aligned} \text{l.x} - \text{k}\_s \frac{\partial T\_s}{\partial \mathbf{x}} = \mathbf{0} \end{aligned} \end{aligned} \tag{13}$$

Microchannel Simulation 195

, 1/2, , , 1/2 , , 1/2 , ,

 

where *A* is a surface area of a face of the cell. From the physical point of view, the convection transfers the information downstream. Therefore, the approximation of convection terms must be weighted to the up-stream-side. The simplest stencil is the first-

1/2, , 1, , 1/2, , ; 0

1/2, , , , 1/2, , ; 0

The diffusion terms of Equation (19) are approximated by using the central difference scheme,

where the grid is assumed to be uniform. Also the pressure gradient terms of the momentum equations are approximated by using central differencing. After the integration

*ij k ijk ijk ijk*

(18)

(19)

1/2, ,

*i jk*

*i jk i jk Ox ui jk* (20)

*i jk ijk Ox ui jk* (21)

(23)

1/2, , , , 1, , 0.5 *i jk ijk i jk p p p Ox* (24)

*face* (25)

(26)

(22)

1, , , , 2

1/2, , 1/2, ,

<sup>2</sup>

*i jk ijk O x*

 

*i jk i jk <sup>V</sup>*

To close the partial differential equation system (1), (2), (3), (4) and (5) for a certain problem, the boundary conditions must be specified. In the present solver, the boundary conditions are handled by using ghost cells, which are illustrated in Figures (2) and (3), the principal idea is to use the ghost cell values to give the fixed boundary value at the boundary of the domain, i.e. the ghost cell values are extrapolated from the face and domain values. Thus,

> 2, , 3, , 2 *j k j k*

 

*x x* 

*face*

*d pA pA x*

*x* 

1/2, , 1/2, , , 1/2,

*FFFS* 

*i jk i jk ij k*

*FFF*

1/2, ,

 

> 

1/2 , ,

*x x*

*i jk*

*p*

over the control volume *Vi,j,k* they are obtained in the form

, ,

where central difference schemes like,

are applied.

respectively.

*ijk*

the Dirichlet and Neumann conditions are given in the form

2, , 3, , *jk jk*

 

*F Au*

*i jk*

order upwind (FOU),

$$0 \le x \le L\_{\infty'}, 0 \le y \le L\_{y'}, z = 0 \implies -k\_s \frac{\partial T\_s}{\partial z} = 0 \tag{14}$$

$$0 \le x \le L\_{\chi}, 0 \le y \le L\_{y'}, z = L\_z \implies -k\_s \frac{\partial T\_s}{\partial z} = 0 \tag{15}$$

at the inner wall surface,

$$\begin{aligned} -k\_s \left( \left. \frac{\partial T\_s \left( \mathbf{x}\_s \mathbf{y}\_\tau z \right)}{\partial n} \right|\_\Gamma \right) &= -k\_f \left( \left. \frac{\partial T\_f \left( \mathbf{x}\_s \mathbf{y}\_\tau z \right)}{\partial n} \right|\_\Gamma \right) \\ T\_{s,\Gamma} &= T\_{f,\Gamma} \end{aligned} \tag{16}$$

where, Equation (10) gives the uniform heat flux boundary condition at the bottom wall of the substrate. Equations (11)–(15) assumes no heat loss from the solid to the environment at the boundary except at *x = 0* for the fluid, where *Tf=Tin*. It should be noted that in reality, heat losses from the heat sink to the environment should be considered by conduction and convection at the inlet and outlet and at the top surface of the heat sink. Heat transfer in the unit cell is a conjugate one combining heat conduction in the solid and convection to the cooling fluid. The two heat transfer modes are coupled by continuities of temperature and heat flux at the interface between the solid and fluid, which are expressed by Equation (16). Γ denotes the perimeter of the inner wall of the channel. Equations (1), (2), (3), (4) and (5) form a closed system from which the flow properties *u, v, w, p* and *T* can be solved as a function of space and time. But, in this study only steady-state flows will be calculated.

#### **2.3 Calculation of incompressible flows**

Simultaneous numerical calculation of Equations (1), (2), (3), (4) and (5) is computationally complex. Therefore, the equations are solved one after another, i.e. in a segregeted manner. The basic structure of Equations (2), (3), (4) and (5) is similar to each other, containing an unsteady term, convection, diffusion and possibly source terms, and they are often called convection-diffusion equations. The flow properties *u, v, w* and *T* are solved from Equations (2), (3), (4) and (5) respectively. Therefore, the continuity Equation (1) is to be modified for pressure or pressure-like quantity. The first stage is to derive a convection-diffusion equation to finite form. In this study the control volume approach is utilized. The process is to be studied by the aid of a general convection-diffusion equation for quantity *Φ*,

$$\begin{aligned} &\frac{\partial}{\partial \mathbf{x}} (\rho u \phi) + \frac{\partial}{\partial y} (\rho v \phi) + \frac{\partial}{\partial z} (\rho w \phi) \\ &= \frac{\partial}{\partial \mathbf{x}} \left( \alpha \frac{\partial \phi}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial y} \left( \alpha \frac{\partial \phi}{\partial y} \right) + \frac{\partial}{\partial z} \left( \alpha \frac{\partial \phi}{\partial z} \right) + S \end{aligned} \tag{17}$$

where α is diffusion coefficient and the source term *S* could contain, for instance, pressure gradient and/or body force, etc. Next, Equation (17) is integrated over the control volume *Vi,j,k.* After the rearrangements the integrated equation can be obtained in the form,

0 ,0 , 0 <sup>0</sup> *<sup>s</sup> xy s <sup>T</sup> xL yLz k*

0 ,0 , <sup>0</sup> *<sup>s</sup> x yz s <sup>T</sup> xL yLzL k*

, , , , *<sup>s</sup> <sup>f</sup>*

 

where, Equation (10) gives the uniform heat flux boundary condition at the bottom wall of the substrate. Equations (11)–(15) assumes no heat loss from the solid to the environment at the boundary except at *x = 0* for the fluid, where *Tf=Tin*. It should be noted that in reality, heat losses from the heat sink to the environment should be considered by conduction and convection at the inlet and outlet and at the top surface of the heat sink. Heat transfer in the unit cell is a conjugate one combining heat conduction in the solid and convection to the cooling fluid. The two heat transfer modes are coupled by continuities of temperature and heat flux at the interface between the solid and fluid, which are expressed by Equation (16). Γ denotes the perimeter of the inner wall of the channel. Equations (1), (2), (3), (4) and (5) form a closed system from which the flow properties *u, v, w, p* and *T* can be solved as a function of space and time. But, in this study

Simultaneous numerical calculation of Equations (1), (2), (3), (4) and (5) is computationally complex. Therefore, the equations are solved one after another, i.e. in a segregeted manner. The basic structure of Equations (2), (3), (4) and (5) is similar to each other, containing an unsteady term, convection, diffusion and possibly source terms, and they are often called convection-diffusion equations. The flow properties *u, v, w* and *T* are solved from Equations (2), (3), (4) and (5) respectively. Therefore, the continuity Equation (1) is to be modified for pressure or pressure-like quantity. The first stage is to derive a convection-diffusion equation to finite form. In this study the control volume approach is utilized. The process is

to be studied by the aid of a general convection-diffusion equation for quantity *Φ*,

*uvw xyz*

 

*Vi,j,k.* After the rearrangements the integrated equation can be obtained in the form,

 

*xxyyzz*

where α is diffusion coefficient and the source term *S* could contain, for instance, pressure gradient and/or body force, etc. Next, Equation (17) is integrated over the control volume

 

*S*

(17)

*T xyz T xyz k k n n*

*s f*

, ,

*s f*

*T T*

only steady-state flows will be calculated.

**2.3 Calculation of incompressible flows** 

at the inner wall surface,

*z* 

> *z*

(14)

(15)

(16)

$$\begin{aligned} F\_{\mathbf{i}+1/2,j,k} - F\_{\mathbf{i}-1/2,j,k} + F\_{\mathbf{i},j+1/2,k} - \\ F\_{\mathbf{i},j-1/2,k} + F\_{\mathbf{i},j,k+1/2} - F\_{\mathbf{i},j,k-1/2} = S\_{\mathbf{i},j,k} \end{aligned} \tag{18}$$

$$F\_{i-1/2,j,k} = \left[ A \left( \rho u \phi - \alpha \frac{\partial \phi}{\partial \mathbf{x}} \right) \right]\_{i-1/2,j,k} \tag{19}$$

where *A* is a surface area of a face of the cell. From the physical point of view, the convection transfers the information downstream. Therefore, the approximation of convection terms must be weighted to the up-stream-side. The simplest stencil is the firstorder upwind (FOU),

$$
\phi\_{\mathbf{i}-1/2,j,k} = \phi\_{\mathbf{i}-1,j,k} + O(\Delta \mathbf{x}); \qquad u\_{\mathbf{i}-1/2,j,k} > 0 \tag{20}
$$

$$
\phi\_{i-1/2,j,k} = \phi\_{i,j,k} + O(\Delta x); \qquad \mu\_{i-1/2,j,k} < 0 \tag{21}
$$

The diffusion terms of Equation (19) are approximated by using the central difference scheme,

$$\left[\frac{\partial \phi}{\partial \mathbf{x}}\right]\_{i+1/2,j,k} = \frac{\phi\_{i+1,j,k} - \phi\_{i,j,k}}{\Delta \mathbf{x}} + O\left(\Delta \mathbf{x}^2\right) \tag{22}$$

where the grid is assumed to be uniform. Also the pressure gradient terms of the momentum equations are approximated by using central differencing. After the integration over the control volume *Vi,j,k* they are obtained in the form

$$\int\_{V\_{i,j,k}} \frac{\partial p}{\partial \mathbf{x}} d\Omega = \left( pA \right)\_{i+1/2,j,k} - \left( pA \right)\_{i-1/2,j,k} \tag{23}$$

where central difference schemes like,

$$p\_{i+1/2,j,k} = 0.5(p\_{i,j,k} + p\_{i+1,j,k}) + O\left(\Delta x^2\right) \tag{24}$$

are applied.

To close the partial differential equation system (1), (2), (3), (4) and (5) for a certain problem, the boundary conditions must be specified. In the present solver, the boundary conditions are handled by using ghost cells, which are illustrated in Figures (2) and (3), the principal idea is to use the ghost cell values to give the fixed boundary value at the boundary of the domain, i.e. the ghost cell values are extrapolated from the face and domain values. Thus, the Dirichlet and Neumann conditions are given in the form

$$
\phi\_{\mathbf{z}\_{2,j,k}} = \mathbf{2}\phi\_{\text{face}} - \phi\_{\mathbf{z}\_{3,j,k}} \tag{25}
$$

$$
\phi\_{\mathbf{z}\_{2,j,k}} = \phi\_{\mathbf{z}\_{3,j,k}} - \Delta \mathbf{x} \frac{\partial \phi}{\partial \mathbf{x}}\Big|\_{\text{face}} \tag{26}
$$

respectively.

Microchannel Simulation 197

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

> 1 1 2 2

*p pp V V x a x ax*

 

1, , , ,

1, , , ,

*i jk ijk*

1/2, , 1, , 1, , , , , ,

*i jk i jk i jk ijk ijk*

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

 

1, , , ,

(27)

*i jk ijk*

*i jk ijk*

*V V*

*a a*

interpolation has been applied to the cell-face velocity as follow,

1/2, , 1, , , ,

*i jk i jk ijk*

*u uu*

**2.5 Pressure correction equation and multigrid technique** 

operators are based on piecewise constant interpolation.

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

1 2

Fig. 2. Ghost cells around the domain.

Fig. 3. The notations of the ghost cells.
