# Modelling And Simulation Of Microchannel Reactor Engineering Essay

Published:

This essay has been submitted by a student. This is not an example of the work written by our professional essay writers.

## INTRODUCTION

With the growing progress in the field of miniaturization a new field known as microelctro-mechanical systems (MEMS) came into picture. The fabrication of MEMS devices called for a new discipline employing fluid flows operating under unusual and unexplored conditions namely micro-fluidics. Micro-fluidics is an area of science and engineering in which fluid behaviour differs from conventional flow theory primarily due to non-continuum effects, surface dominated effects, and low Reynolds number effects induced by the small length scale of the microflow systems.

The classification proposed by Mehendale et al .[1] and Kandlikar and Grande[2] categorized the range from 1 to 100 Î¼m as micro-channels. In micro-fluidics, theoretical knowledge for gas flows is currently more advanced than that for liquid flows. Concerning the gas flow through micro-channels, the issues are actually more clearly identified; the main micro-effect that results from shrinking down of device size is 'rarefaction'. In the continuum fluid transport theory, governed by the Navier-Stokes equations, it is assumed that the state variables do not vary appreciably over the length and time scales compared with the molecular mean free path and molecular relaxation time. The local density oscillation near the solid-liquid interface of the micro-channels, significant deviation of liquid viscosity compared with the bulk value, may not necessarily mean the breakdown of continuum theory; at the same time, it is important to understand how the continuum theory works in a micro-flow. There were no evidences that continuum assumptions were violated for the micro-channels tested, most of which had hydraulic diameters of 50 Î¼m or more. There is a clear need for additional systematic studies that carefully consider each parameter influencing transport in micro-channels.

Before going into detail on the transport processes it would be wise to discuss the differences in a micro to macro reactor w.r.t. the optimum process efficiency. Micro-channel reactors are increasingly used in many fields of industry due to the capabilities exceeding those of traditional macro-scale reactors. The high heat and mass transfer rates in micro-channel reactors allow the reaction to be performed under more aggressive conditions with higher yields (Jensen, 2001). The surface to volume ratio can be higher up to 10,000-50,000 m2 m-3 (Minsker and Renken, 2005) providing drastically higher heat and mass transfer rates than the traditional chemical reactors. Thus, micro-reactors can remove heat much more efficiently than traditional chemical reactors and can perform safely for highly exothermic or endothermic reaction. Therefore, high reaction temperature is possible with micro-channel reactors and leads to reduced reactor volumes. Moreover, less amount of catalyst used improves the energy efficiency and reduces the operational costs.

Another benefit of micro-channel reactors is that if the system fails the amount of accidentally released chemicals is rather small and it could be easily controlled. The integrated sensor and control units could allow the failed reactor to be isolated and replaced while other parallel units continued production (Jensen, 2001). Besides the benefits of the micro-channel reactors, some problems still remain in the system. Clogging is one of the main problems occurring in the particle containing processes, for example, catalytic reactions. Clogging has been identified as the biggest problem, however, the particle containing processes also cause rather high pressure drop through the flow channel. These two main factors have to be considered while operating with micro-channel reactor.

The following study is an attempt on describing incompressible flow through micro-channels by analytical method and hence comparing the results with that obtained from a commercial software package for simulation (FLUENT). For liquids and gases the fluid particle need not be as large as 10 ïm, but rather be on the order of 10 nm and a few hundreds of nanometers, respectively, to be considered as macroscopic. Therefore for ordinary liquids (fluids), the current micro fluidic devices are subjected to the rules of classical fluid mechanics.

## 1.1 A brief introduction to FLUENT:Â

Fluent is a general-purpose CFD code based on the finite volume method based on a collocated grid. FLUENT technology offers a wide variety of physical models that can be applied to a wide array of industries. Some of its features include the following:

Dynamic and Moving Mesh: The user simply sets up the initial mesh and prescribes the motion, while FLUENT software automatically changes the mesh to follow the motion prescribed. This is useful for modelling flow conditions in and around moving objects.

Heat Transfer, Phase Change, and Radiation: FLUENT software contains many options for modellingÂ convection, conduction and radiation.

Multiphase: It is possible to model several different fluids in a single domain with FLUENT.

Turbulence: A large number ofÂ turbulenceÂ models are used to approximate the effects of turbulence in a wide array of flow regimes.

Acoustics: TheÂ acousticsÂ model lets users perform "on-the-fly" sound calculations.

Reacting Flows: FLUENT technology has the ability to model combustion as well as finite rate chemistry and accurate modeling of surface chemistry.

Post-processing: Users can post-process their data in FLUENT software, creating - among other things - contours, pathlines, and vectors to display the data.

Before simulating the problem we first need to model the problem with software named GAMBIT. Speaking of GAMBIT, it can be well said as a geometric modelling and grid generation tool i.e. often shipped with FLUENT technology. It allows the users to create their own geometry or import geometry from most CADÂ packages. It can automatically mesh surfaces and volumes while allowing the user to control the mesh through the use of sizing functions and boundary layer meshing.

## Transport Process in Micro channel reactor

Apart from wetting, adsorption, and electro kinetic effects, the slip phenomena play an important role in modelling liquid flow through micro-channel. Helmloltz and Von Piotrowski found evidence of slip between a solid surface and a liquid and later Brodman verified their results. Navier (1823) was the first to model partial slip at the wall for liquids well before Maxwell's slip model for gases (1879). Analytic slip-flow studies have traditionally been based on the continuum form of Navier-Stokes equations and energy equation with the slip-flow effects concentrated in the additional terms in the tangential velocity and temperature boundary conditions [3]. These new boundary conditions represent velocity slip and temperature jump conditions at the gas-surface interface. Slip flow condition is well described by introducing the concept of slip length. Slip length is the distance behind the solid-liquid interface at which the velocity extrapolates to zero. Liquid flow in a micro-channel becomes fully developed after a short entrance length so that it can be modelled as a two-dimensional flow. In view of this, present study aims for two-dimensional micro-channel flow simulation with non-continuum (slip velocity) boundary conditions.

The slip velocity at the wall as given by Navier is

wall = wall

Where Ls is the slip length; n is the normal coordinate pointing inward from the channel wall.

The above proposition is valid both for gas and liquid flow.

## 2.1 Modelling of the governing equation for the one dimensional incompressible flow through the micro-channel with a two-dimensional velocity field and Temperature distribution:

Assumptions:

(1) The fluid is Newtonian in nature with constant viscosity and is incompressible.

(2) The flow is laminar and the shear rates involved are small.

(3) Gravity effects are neglected.

(4) We also assume that u = u(x, y) ; v = w = 0 and T= T(x,y).

Let us consider no slip flow condition and model the governing equations very close to the wall surface.

By boundary layer theory the flow near a solid surface with no slip condition is given by: u = [

where, u: velocity in x-direction (function of both x & y)

Î´= Î´(x): boundary layer thickness

The above equation is derived considering the profile to be a polynomial function of the distance from the wall (y) satisfying the specified boundary conditions.

Let the best suited curve be given by

u=c1 + c2 y +c3y2 +c4 y3 (1)

The boundary conditions are:

B1: At y=0; u=0 : No slip condition;

B2: At y=b/2; u=umax = : Bulk velocity;

B3: At y=0; : Constant shear at the wall (at a particular x);

B4: At y= Î´; : Shear at ;

Substituting B1 in (1); c1=0

Using B3; c3=0

Using B2; c2 (b/2) + c4 (b/2)3= umax =

Using B4; c2+3c4 Î´2 = 0

Solving for c2 and c4; we get,

c4= -umax /2 Î´3 and c2= 3umax /2 Î´

Substituting in equation (1); we get,

u = [

Now considering slip flow condition, the slip velocity, uw, defined by Eqn , on the channel walls is based on the following assumptions:

(1) At a distance of Ls from the present channel wall, on both the sides, there exists an imaginary boundary wall (IBW) - rigid and solid.

(2) At the position where presently there exists a wall, it is assumed that no solid wall is present and there exists a velocity for the fluid moving over the IBW.

The velocity is equal to the slip velocity in the present case.

Considering the above steps, and solving for the wall velocity, we get

Above equation is in general considering an IBW.

Now at the real boundary, Ls distance from the IBW slip velocity prevails. Hence, replacing y with Ls (Ls = b/100) would yield us the slip velocity as:

The problem now is to find the free-stream velocity, Uâˆž, and the boundary layer thickness, Î´. As we are dealing with micro-channels, it can be expected that the velocity of the fluid on the IBW shall be small enough to consider the flow as laminar. This consideration will yield the expression for the boundary layer thickness, Î´, as:

The free-stream velocity (Uâˆž) is nothing but the incident velocity. Thus, combining the above expressions, we get an expression for the slip velocity as a function of Î´, that is, a function of x. mathematically, it can be derived that the slip velocity is actually an inverse function of axial dimension. The radial velocity component is assumed to be zero on the channel walls.

At the entrance of the channel, the radial velocity component is set to zero, and a uniform velocity Uâˆž for the axial velocity component is assumed.

Because of the elliptic nature of the flow fields investigated, the outlet boundary condition will have some influence on the development of the upstream flow. The normal gradient of all dependent variables at the outlet are constants. This is equivalent to the requirement of the second derivative of any dependent variable becoming zero. Using this condition, the numerical procedure is stable and convergent for those flows, which eventually reaches a steady state condition at the outlet.

The analytic equation obtained, forms the base for solving the flow equations. In order to predict the phenomena, the steady state equations are solved by using the commercial software package FLUENT which uses a numerical iterative method for the simulation.

The above equation explains the flow phenomena in the boundary layer region. Unlike flow in macro channels, there exists a slip flow at the boundary wall for micro channels; hence the boundary conditions to the above derivation can be modified to:

B1': At y=0; u= uwall;

B2': At y=b/2; u= umax â‰ uâˆž

Here umax is function of axial distance(x) from the entrance. Substituting these along with B3 and B4 in equation (1) the velocity profile comes out to be

(2)

Similarly the temperature profile in the micro channel can be reasoned to be of the form

(3)

Where, Twall: wall temperature;

Tâˆž: bulk temperature (temperature of the entering fluid)

Î´t: thermal boundary layer thickness

The thermal boundary layer thickness is again a function of momentum boundary layer which can be determined by considering energy balance around a control volume element in the flow boundary.

Apart from the empirical approach, the velocity profile can also be determined by solving the continuity and the momentum balance equation. The postulates for the flow phenomena are

vx=u(x,y); vy0 :negligible radial velocity ; P=P(x) : Pressure only function of x

with boundary conditions :

B1: At y=-Ls; u=0 : slip length(Ls) distance away from wall

B2: At y=Î´ ; : shear at Î´

Hence, continuity equation:

and momentum equation:

Ï (

equation reduces to (assuming vy negligible)

=0

Solving for vx using B1 and B2 we get,

(4)

Now consider flow in a micro channel as shown

Both the equations (2) and (4) equally describe the flow very near to the wall but at yâ†’b/2, radial velocity component is induced and hence flattening of the velocity profile occurs. This results in the decreasing magnitude of the axial velocity component after a peak maxima at about xo (where b/2=Î´).

Expression for thermal boundary layer thickness (Î´t):

For flow over a heated boundary surface the following assumptions are made:

T=T(x,y); : considering linear variation with x

Hence the energy balance equation reduces to

(5)

Where : Dissipation function =

Let the velocity profile be given by equation (2) and temperature by equation (4), hence by substituting and integrating w.r.t. y , we get

]0y +

For y=Î´t ;

And y=0; wall flux

This is the integral equation of the boundary layer for constant properties and constant free stream temperature Tâˆž. Since laminar flow is considered, the viscous dissipation term is very negligible and hence is neglected unless the velocity in the flow field becomes very large.

Where, p= ; ;

Neglecting higher order terms of p,

Hydrodynamic boundary layer thickness (is obtained by substituting the velocity profile into the von Karman integral balance to get

Hence we get,

At x=0; p=0 => C=0

(6)

Similarly by considering the velocity profile as in equation (4), the thermal boundary layer thickness is found to be:

(7)

Therefore substituting either of these equations in equation (3) we get the required temperature profile for slip flow over a heated boundary surface.

Defining the Problem Statement:

In the current simulation we aim at studying the flow phenomena and the temperature distribution occurring in a micro-channel. With reference to the simulations done in [1], we choose the same channel dimensions so as to have a discussion on the results obtained.

The channel considered is viewed as two parallel plates placed one over the other extending into the z-direction, and the flow is in the x-direction. It has been assumed that the z-directional dimension is equal to unity. The channel clearance length (b) is 1mm. As far as the bulk properties of the fluids are concerned, 1mm is still a scale governed by the classical laws of simple liquids and ideal gases under normal pressure. Therefore, for ordinary liquids (fluids), the current micro fluidic device is subjected to the rules of classical fluid mechanics. The channel length (L) considered is equal to 10cm. Again in reference to [1] the inlet Reynolds number is taken as 10.

Approach to the simulation:

The above 2D problem geometry is created (to the given dimensions) by using GAMBIT software. Then the geometry is meshed (100*100, let) into smaller cells. Higher the no. of cells the greater is the accuracy with more rigorous calculations and lower convergence limit.

This then is exported as a 2D mesh file to be solved in FLUENT.

Now specifying the boundary conditions and choosing a suitable solver we initialize the problem and then run the iteration for the solution.

## FLUENT SIMULATION

As discussed above the mesh file is imported into FLUENT and the boundary values are specified. The working fluid considered is liquid water with outlet pressure of 1 atm.

In specifying the boundary conditions,

The inlet velocity corresponding to a Reynolds no. of 10 is calculated as 0.01m/s.

At the wall the slip velocity is accounted by specifying the shear stress in x direction.

The solver specified uses a pressure correction based iterative SIMPLE algorithm with the discretization done with 1st order upwind scheme.

Theoritically the shear stress at the wall is calculated by the following expression:

Since the wall velocity is a function of the position x, so as a first approximation the shear value specified is a constant averaged over the total length of the channel.

Alternately, the function generated for shear stress in terms of the position variable, x, can directly be coupled by defining a user defined function (udf). The udf coded is as follows:

#include "udf.h"

DEFINE_PROFILE(wall_shear_x,thread,i)

{ double x[ND_ND];

face_t f;

double y;

begin_f_loop(f,thread)

y=x[1];

F_PROFILE(f,thread,i)= -0.0000000001* (0.15* sqrt(10000* (y+0.0000001))-0.006/sqrt(10000*(y+0.0000001)))-0.00001*(0.0237* sqrt(10000*(y+0.0000001))-0.000000003795/sqrt(10000*(y+0.0000001)));

end_f_loop(f,thread)

## }

The values hence obtained are in close resemblance to that calculated theoretically. Other methods for shear stress evaluation involve defining a line/rake at the required surface and hence obtaining the wall fluxes (in XY plot option). For this we modify the considered micro-channel by increasing the channel width by twice the slip length (Ls) (i.e. the imaginary boundary layer). Doing so we solve for no slip flow condition and subsequently after simulation we evaluate the wall fluxes at y = Ls.

The wall shear stress obtained by the above methods are summarised as follows:

Interpretation:

It is found that there exists no y-shear at the wall which is physically not true. This is because unlike flow through macro-channels, the velocity here is also a function of x even in the fully developed regime. This induces a velocity gradient in y direction though practically very negligible in comparison to x directional shear.

Existence of dominant y-shear stress can be shown by introducing radial velocity component into the micro-channel, therefore consider an extended micro-channel as shown:

On simulation with inlet velocity of 0.01 m/s, the wall fluxes obtained are:

Hence it is concluded that the computed x-shear by the software is a summation of both shear stresses acting on plane parallel to x axis and y-axis in x-direction.

The velocity profile hence obtained shown below.

## Results and discussion:

Comparing with the simulation result obtained in reference [1].

They refer to a phenomenon called flattening in the velocity profile when 1D flow is concerned. For the same characteristics dimension of the micro-channel with slip boundary condition and for Re= 10, the non-dimensional axial velocity develops to a value of 1.503 at x=1.21 and the flattening in the velocity profile begins onwards. Fig.1 shows the graph between axial velocity and channel clearance at different x values is shown as (from the reference):

But from the above fluent simulation we get the plot (fig.2).

In comparing both the plots we find the reference plot seems to be physically incorrect with their curve starting from zero velocity and regarding the velocity magnitudes, the maximum comes out to be equal to 1.038*10^-2 m/s.

Fig.1

Fig.2

## Fluent Simulation:

Consider a micro-channel of width (b) = 1mm and length (L) =10cm. This geometry is then modelled in gambit and exported for fluent simulation. Prior to the simulation the operating parameters and the boundary conditions are specified. The working fluid considered is liquid water with outlet pressure of 1 atm.

As discussed above the mesh file is imported into FLUENT and the boundary values are specified. The working fluid considered is liquid water with outlet pressure of 1 atm.

In specifying the boundary conditions,

The inlet velocity corresponding to a Reynolds no. of 10 is calculated as 0.01m/s

At the wall the slip velocity is accounted by specifying the shear stress in both the directions.

The solver specified uses a pressure correction based iterative SIMPLE algorithm with the discretization done with 1st order upwind scheme.

## Temperature profile in a microchannel:

By boundary layer theory the flow near a solid surface with no slip condition is given by:

u = [

where, u: velocity in x-direction (function of both x & y)

Î´= Î´(x): boundary layer thickness

The above equation is derived considering the profile to be a polynomial function of the distance from the wall (y) satisfying the specified boundary conditions.

Let the best suited curve be given by

u=c1 + c2 y +c3y2 +c4 y3 (1)

The boundary conditions are:

B1: At y=0; u=0 : No slip condition;

B2: At y=b/2; u=umax = : Bulk velocity;

B3: At y=0; : Constant shear at the wall (at a particular x);

B4: At y= Î´; : Shear at ;

Substituting B1 in (1); c1=0

Using B3; c3=0

Using B2; c2 (b/2) + c4 (b/2)3= umax =

Using B4; c2+3c4 Î´2 = 0

Solving for c2 and c4; we get,

c4= -umax /2 Î´3 and c2= 3umax /2 Î´

Substituting in equation (1); we get,

u = [

The above equation explains the flow phenomena in the boundary layer region. Unlike flow in macro channels, there exists a slip flow at the boundary wall for micro channels; hence the boundary conditions to the above derivation can be modified to:

B1': At y=0; u= uwall;

B2': At y=b/2; u= umax â‰ uâˆž

Here umax is function of axial distance(x) from the entrance. Substituting these along with B3 and B4 in equation (1) the velocity profile comes out to be

(2)

Similarly the temperature profile in the micro channel can be reasoned to be of the form

(3)

Where, Twall: wall temperature;

Tâˆž: bulk temperature (temperature of the entering fluid)

Î´t: thermal boundary layer thickness

The thermal boundary layer thickness is again a function of momentum boundary layer which can be determined by considering energy balance around a control volume element in the flow boundary.

Apart from the empirical approach, the velocity profile can also be determined by solving the continuity and the momentum balance equation. The postulates for the flow phenomena are

vx=u(x,y); vy0 :negligible radial velocity ; P=P(x) : Pressure only function of x

with boundary conditions :

B1: At y=-Ls; u=0 : slip length(Ls) distance away from wall

B2: At y=Î´ ; : shear at Î´

Hence, continuity equation:

and momentum equation:

Ï (

equation reduces to (assuming vy negligible)

=0

Solving for vx using B1 and B2 we get,

(4)

Now consider flow in a micro channel as shown

Both the equations (2) and (4) equally describe the flow very near to the wall but at yâ†’b/2, radial velocity component is induced and hence flattening of the velocity profile occurs. This results in the decreasing magnitude of the axial velocity component after a peak maxima at about xo (where b/2=Î´).

Expression for thermal boundary layer thickness (Î´t):

For flow over a heated boundary surface the following assumptions are made:

T=T(x,y); : considering linear variation with x

Hence the energy balance equation reduces to

(5)

Where : Dissipation function =

Let the velocity profile be given by equation (2) and temperature by equation (4), hence by substituting and integrating w.r.t. y , we get

]0y +

For y=Î´t ;

And y=0; wall flux

This is the integral equation of the boundary layer for constant properties and constant free stream temperature Tâˆž. Since laminar flow is considered, the viscous dissipation term is very negligible and hence is neglected unless the velocity in the flow field becomes very large.

Where, p= ; ;

Neglecting higher order terms of p,

Hydrodynamic boundary layer thickness (is obtained by substituting the velocity profile into the von Karman integral balance to get

Hence we get,

At x=0; p=0 => C=0

(6)

Similarly by considering the velocity profile as in equation (4), the thermal boundary layer thickness is found to be:

(7)

Therefore substituting either of these equations in equation (3) we get the required temperature profile for slip flow over a heated boundary surface.

Fluent Simulation

Consider a micro-channel of width (channel clearance) b= 1mm and length (L) = 10cm. The

VR = and

TR=

Then,

=Î¸âˆž(uâˆž-uwall)]

Where, Î¸âˆž= ( ; L=

[L+VR-L (TR)-(VR) (TR)]

Assuming that Î´t < Î´

=[L+VR-L(TR)-(VR)(TR)]

= [[Ly+ (3y2/4 Î´)-(y4/8 Î´3)] +L [(3y2/4 Î´t)-(y4/8 Î´3)]-[(9y3/12 Î´Î´t) + (y7/28 Î´ 3Î´3)-(3y5/20 Î´3 Î´t)-(3y5/20 Î´t3 Î´)]]

Let Î´t / Î´ = p

= [pÎ´L+ (3p2 Î´/4)-(p4 Î´/8)]-L [(3p2 Î´/4)-(p4 Î´/8)]-[(3p2 Î´/4) + (p4 Î´/28)-( 3p4 Î´/20)]

=L Î´ [p-(3p2/4) + (p4/8)] + Î´ [(3p2/20)-(3p4/280)]

As in the case of a micro channel Î´/ Î´t = p is very small

We can neglect higher order terms, so the integrand is L Î´p

Î¸âˆž (uâˆž-uwall) [] =

Î¸âˆž (uâˆž-ub)=3Î¸âˆžÎ¬/2 Î´t

=3Î¬/2 pÎ´ (uâˆž-uwall)

L [Î´ (dp/dx) +p (dÎ´/dx)] = 3Î¬/2 pÎ´ (uâˆž-uwall)

PÎ´2(dp/dx) +p2Î´ (dÎ´/dx) = 3Î¬/2ub

We know that

Î´dÎ´ =140Ï…dx/(13(uâˆž-ub))

Î´2 =280Ï…x/(13(uâˆž-ub))

PX (dp/dx) +p2/2 = 3Î¬/2ub[(13(uâˆž-ub))/280Ï…]

=39/(560 Pr L)

Homogeneous equation is

PX (dp/dx) +p2/2 =0

=> pdp/p2 + dx/2x =0

ln p= -0.5ln x + ln c

=> p=cx-0.5

Solution of the original differential equation is

p= cx -0.5+ 39/ (560 Pr L)

At x = 0; p= 0

=>p=39/ (560 Pr L)