Scientific paper
Modeling and Finite Difference Numerical Analysis of Reaction-Diffusion Dynamics in a Microreactor
Igor Plazl1* and Mitja Lakner2
1 Department of Chemical Engineering, Faculty of Chemistry and Chemical Technology, University of Ljubljana, A{ker~eva 5, SI-1000 Ljubljana, Slovenia
2 Civil and Geodetic Faculty, University of Ljubljana, Jamova 2, SI-1000 Ljubljana, Slovenia
* Corresponding author: E-mail: igor.plazl@fkkt.uni-lj.si Phone: +3861 2419 512; Fax: +3861 2419 530
Received: 16-04-2009
Dedicated to the memory of the late Prof. Dr. Valentin Koloini
Abstract
A theoretical description with numerical experiments and analysis of the reaction-diffusion processes of homogeneous and non-homogeneous reactions in a microreactor is presented considering the velocity profile for laminar flows of miscible and immiscible fluids in a microchannel at steady-state conditions. A Mathematical model in dimensionless form, containing convection, diffusion, and reaction terms are developed to analyze and to forecast the reactor performance. To examine the performance of different types of reactors, the outlet concentrations for the plug-flow reactor (PFR), and the continuous stirred-tank reactor (CSTR) are also calculated for the case of an irreversible homogeneous reaction of two components. The comparison of efficiency between ideal conventional macroscale reactors and the microreactor is presented for a wide range of operating conditions, expressed as different Pe numbers (0.01 < Pe < 10). The numerical procedure of complex non-linear systems based on an implicit finite-difference method improved by non-equidistant differences is proposed.
Keywords: Microfluidics; reaction-diffusion dynamics; microreactor; numerical analysis; non-equidistant finite differences
1. Introduction
Microtechnology has uncovered new scientific solutions and challenges in a broad range of areas, from electronics, medical technology, and fuel production and processing to biotechnology, chemical industry, environmental protection, and process safety. Microscale reactors are devices whose operations depend on precisely controlled design features with characteristic dimensions from submillimeter to submicrometer. Because of the small amounts of chemicals needed and the high rate of heat and mass transfer, microscale systems are especially suited for reactions with highly toxic, flammable, and explosive reactants.1-4 In the last decade, Microreactor Technology (MRT), accepted as a new concept in chemical engineering, has impressively demonstrated the advantages of mi-crostructured devices for chemical and biochemical reactions. Numerous reactions, including many notable and in-
dustrially relevant reactions, have been tried out successfully in microreactors and several hundred publications with some research review papers have appeared in peer-reviewed journals.5-10 The small length scale of microreac-tors reduces transport limitations, giving nearly gradient-less conditions desirable for the determination of reaction kinetics. The reactor miniaturization allows us to carry out reactions under more precisely controlled conditions than with conventional macroscale reactors, leading to a possibility of improved yield and selectivity of the desired products. 11-13 In the area of catalytic chemistry, microreactors are an extremely efficient tool for rapid catalyst screening and for combinatorial chemistry.14-18 In order to simulate chemical and biochemical processes in a continuous flow microreactor a coupled system of convection-diffusion-reaction in combination with hydrodynamics has been developed and many studies on microfluidic systems are available.19-26'27'28 Characterization of micro phenomena,
like frictional pressure drop and viscous dissipation effects in microchannels are subject of continuous debate. The discrepancies among the work of many researchers have been summarized in review papers.29-31 The use of the Na-vier-Stokes equations appears to be appropriate for microchannel flows of liquids as long as the hydraulic diameter (Dh) of system is greater than 100 pm for conduits filled with Newtonian fluids such as water. Non-Newtonian fluid effects are expected to be important for polymeric liquids and particle suspensions flows. Wall slip effects are negligible for liquid flows in microconduits and viscous dissipation effects on the friction factor are negligible in smooth microchannels, especially for conduits with Dh > 100 pm.29 For the continuum assumption to hold in gaseous media the characteristic length should be large enough and most researchers agree on classifying the flow to be in continuum for Knudsen number (Kn) < 10-3. However, when Kn > 10, molecular level modeling is required to describe the behavior of the fluid.32,33
The present work is a theoretical and numerical outcome of an ongoing research on characterizing microreac-tor performance for different homogeneous, non-homogeneous, and heterogeneously catalyzed reactions with nonlinear reaction kinetics. The parallel flow of two immiscible fluids with the enzyme catalyzed reaction at the interphase surface in Y-shaped glass microchannel was investigated at the conditions when the continuum assumption can be applied.29 It is an attempt to analyze and to forecast the behavior of reaction-diffusion dynamics in a continuous flow microchannel with the application of a relatively simple finite difference technique. An implicit finite-difference method was used to solve the systems of 2D and 3D partial differential equations of a second order. Special attention was devoted to non-equidistant differences in order to improve the stability and accuracy of the numerical solutions. It is shown that the accuracy of the approximations can be enhanced by mesh refinement in the vicinity of the domain parts where the variations in field unknowns are expected to be significant. The numerical experiments are analyzed and examined for accuracy, stability, and theoretical consistency. To investigate the effectiveness of microreactor performance, we compared the di-mensionless outlet concentrations of a general second-order homogeneous reaction in a microchannel with that in the ideal plug-flow and continuous stirred-tank reactor for a broad range of operating conditions.
1. 1. Theoretical Background
Let us consider first an irreversible homogeneous reaction taking place in the microchannel within the flow driven at a parabolic velocity profile, developed in the smallest y-dimension. Reactive components A and B enter the channel in two parallel flows as schematically presented in Figure 1. Regarding the operating conditions typical of microchannels, the laminar fluid flow in the x-axis di-
rection is considered (Fig. 1). Flow in microchannel is predominantly laminar.1 Molecular effects also become more significant in microchannel when the characteristic length decreases to the point that the continuum assumption becomes invalid.33 Jahnisch and co-authors3 observed the laminar flow for Re numbers up to 2000. However, many researchers realized turbulent flows for much lower Re numbers. Turbulence effects become very important for Reynolds numbers above 1000, mainly due to (upstream) geometric non-uniformities.29 It was shown that the channel aspect ratio and the angle of merging of two inlet channels substantially influence the critical Reynolds number. The results, obtained for smooth glass channels with relative roughness around 1%, shown a fast decrease of the critical Reynolds number (from around 2000 at aspect ratio of 2), until the aspect ratio reaches a value of 6 (critical Reynolds number is around 410).34
At laminar flow conditions, the velocity profile fully developed in the direction of the least dimension ([W,-W]) can be described as a function of y position only (Eq. 1):
(1)
The reaction formula and the rate equation of an irreversible homogeneous reaction proceeding in the micro-reactor are as follows:
a\ + pR -> p?, {-rA)= kc"Ac*B
(2)
where A and B are the reactants and P product; (-rA) and k are the reaction rate (kmol m-3s-1), and the rate constant (m3(a+p-i k mol-(a+3-1)s-1), respectively; cA and cB are the molar concentrations of the components A and B (kmol m-3s-1); and a and ¡3 are the reaction orders in the kinetics term.
Figure 1. Scheme of the microchannel with a parabolic velocity profile developed in the smallest dimension.
The mass conservation equations, containing convection, diffusion, and the reaction terms for components A and B at steady-state conditions in a dimensionless
form with the associated boundary conditions for 2D partial differential equations are given by:
Component A
V 'dc PeA
d2X d-X +
dw3
-Da,X"Y''
(3)
b.c.
, , 17, 0 <[//<}
L J J
[0, — 1 <y/ <0
|x-{m)=o> m^
ci//	W
(4)
V	PeE
-l<if/< I
8JY S'Y + ■
dtp1
-DasXaYp (5)
b.c.
(6)
—(i,±/) = 0; 0<4<-dy	W
—f—= ^
where the dimensionless concentrations and coordinates are defined as
Cjrt Csii if W
(7)
The dimensionless concentration profiles of both reactants in the microchannel depend on the following di-mensionless numbers

Wkc",
W k cAOcB
(8)
In the Eqs. (3-7) DA and DB are the diffusion coefficients of the components A and B (m2 s"1), vmax is the maximum velocity of the parabolic velocity profile (m s"1), and W is the half width of the channel (m), cA 0 and cB 0 are the feed concentrations of the components A and B (kmol m3). The dimensionless number, defined for each of the components j, Daj represents the ratio of the characteristic residence time or fluid motion time scale to the characteristic reaction time, and is called the Damkohler number.35 The Peclet number, Pe, provides an indication of the relative importance of diffusion (random thermal motion of molecules within their surrounding environment) and convection (the transport as a result of bulk motion of a fluid).21,36 The ratio of the molecular diffusivity to the product of the average velocity of the flow and the characteristic length of the system perpendicular to the direction of the flow is defined in the literature as well as the inverse value of the product of Reynolds and Schmidt number, where r is density (kg m~3), and n is dynamic viscosity (Pa s).37,38
In order to analyze the reaction-diffusion dynamics of more realistic and industrially relevant processes in mi-crostructured devices, the enzyme reaction process between two immiscible fluids is considered and simulated at steady-state conditions. As stated in several literature reports,1,6 the fluid flow in microstructured devices is typically laminar and enables even two miscible fluids to flow parallel next to each other.27,39 Beside the reaction at the defined interface with short residence times, the efficient phase separation with the aqueous-organic interface formed in the middle of all regions of the Y-shaped microchannel is an important feature that allows the reuse of aqueous flow with enzyme. The simulation of the heterogene-ously catalyzed reaction, as schematically presented in Fig. 2, demands the velocity profile of the two-phase fluid within the microchannel. Considering the presumptions of
15 mm
332 mm 220 n»t W
a)
organic --wier
phase phase
organic phase + product + residuaIA
Figure 2. Graphical presentation of the microfluidic device used in simulations: a) scheme of the main channel; b) line drawing of the whole mi-croreactor with indicated inflow and outflow constituents (dimensions of microchannel used in calculations: width=220 pm, deep = 50 pm, length = 33.2 cm).
the position of the interface area in the middle of the horizontal microchannel and by neglecting the compressibility and gravitational force, dimensionless x-momentum equations for a fully developed Poiseuille-type flow of water and organic phase have to be solved.39
As shown in Fig. 2, model set-up comprised a separate inflow of water phase with a dissolved enzyme, and organic phase containing reactants A and B. In the model description the diffusional transfer of the reactant B across the water-organic interphase was considered, while transport of reactant A into the water phase was neglected. Arising from this, dimensionless 3D partial differential equations for steady-state conditions in the single pass mi-croreactor system with the associated boundary conditions are as follows27:
8X,lh D.lk( d'XAlh 82X,,h
34
oç
W
3f
cy
s2xAlh
ddr
w
8:Y„
c'Y
w J a
d¿¡~ Ôi//: dor
vAyr.wy

84 W \ dÇ* cy- 8(0-
(9)
(10)
(11)

dy/	WW
da	W
b.c. (Eq. 12):
w diff	w
M M^fi)	L	ft
W	dm	W
8X

8ctt
= 0, 0<£ <—, -/<(✓<«; w
(12)
0<t <—, f<£»<—; W	w
where variables YB/w and YB/h are dimensionless concentrations of reactant B (cB/cB 0) in water (w) and organic phase (h), respectively, and XA^h is the dimensionless concentration of reactant A (cA/cA 0) in the organic phase (/); DA/h, DB/h and DB/w are molecular diffusion coefficients for reactants in both phases (m2/s) and vr is rate of enzyme reaction (mol/m3s); KP is the partitioning coefficient for reac-tant B in organic-water system.
For the selected enzyme model reaction (i.e., enzyme catalyzed isoamyl acetate synthesis), the non-linear kinetic expression for reaction rate assuming the inhibition by reactant B was used:40
dXjJÇA®) r-WÂiji . , L , ^ H dV	cA,oDA,b	W	W
W dy/	W
Figure: 3. 2D and ■ 3D
non-equidistant partition of grid points.
(13)
where dimensionless concentrations of reactant B (cB,0 is inlet concentration) in water YB/w and reactant A (cA,0 is inlet concentration) in organic phase XA/h (/) were considered. vrmax is maximum reaction rate (mol/m3s), KB and KA (mol/m3) are binding constants for B and A, respectively, while Ki B (mol/m3) is the B inhibition constant.
1. 2. Numerical Experiments and Analysis
The finite differences were used to replace the partial derivates in the model equations (Eqs. 3-12). The discretization was done by the finite differences on a 2D and 3D Cartesian grid which demands the implicit approach of solution. The result of discretization of the partial derivatives is the system of non-linear algebraic difference equations for the dependent variables at each grid point. The formulation of the mathematical model including all reaction and flow situations in the microchannel can only be solved by rigorous numerical procedures, especially for the highly irregular and complex geometry of the domain, and various attempts are being made towards investigating the complex and non-linear interactions between convection-diffusion and reaction.20 For example, Micro-nit Microfluidics BV from the Netherlands offers glass microreactor chips with a standard channel width/deep 150pm/150pm and channel length 332 mm or 676 mm, meaning that the ratio of width to length in the grid domain is approximately 1 to 2000 or 4500. In two, and definitely in three spatial dimensions many simulations based on numerical solutions of such problems can hardly be performed with a sufficient spatial resolution on a static equidistant computational grid, and require dynamic, adaptive grids. Many authors have recognized that mesh adaptation can be an effective tool for simulating sharp fronts or moving interface problems, reducing numerical dispersion and oscillation, as well as computational costs and data storages, without reducing the overall level of accuracy. It has been demonstrated that significant improvements in accuracy and efficiency can be made by adapting the mesh nodes, so that they remain concentrated in regions of sharp fronts or interfaces.41-44
In order to improve the resolution, in particular at certain locations of the computational domain by considering the steady-state behavior of the presented problems, a non-equidistant partition of grid points were developed in this work. The non-equidistant partitions of grid points for 2D and 3D domains of microchannel are schematically presented in Figure 3.
Simple mathematical manipulation is needed to transform the static equidistant of central difference form with error of O(A<f}
2A$
(14)

A?
to the non-equidistant finite difference (Eq. 15), as is demonstrated for the central finite difference on 2D domain
+ ¿4,
f
'(•? 1S5 [ Xj+l	^i Xj-l I	j
(15)
where the 2D non-equidistant grid is constructed by the formulas:
i = 0,....n-1 and y, - sign(jl— 1m
; j'^'-m.....in
The improved efficiency of the applied non-equidistant finite difference can be graphically observed at the same partition of grid points in Figure 4.
The stability and the accuracy of the improved numerical procedure were examined with numerical experiments for cases without chemical reaction (the reaction rate constant in the Eqs. 3-8 was set to zero, k = 0). While the component A continuously enters at one half of the microchannel and component B at the other half of the microchannel, the average concentrations of both components at the outlet of the microchannel with a high ratio length/width are known as the result of the convection-diffusion processes only. In Figure 5.a the dimensionless concentration profile for component A is presented for cases when the dilutions of the components take place in a single-pass flow through the microchannel. The predicted di-mensionless outlet concentration profiles are in excellent agreement with the reasonably expected outlet concentrations of both components, cAou/cA0 = 0.5 (Fig. 5.a,b).
The result of the discretization of the system of 2D and particularly 3D partial differential equations inclu-
Figure 4. The example of graphical results of the dimensionless concentration profiles for the component A (Eqs. 3-6, and 8) for the very first part of the microchannel (width/length = 2/6) at the static equidistant and non-equidistant partition of grid points with the same overall mesh density (nxm = 60 X 80).
ding convection-diffusion-reaction terms combined with the microfluidics description (Eqs. 3-8 with the arbitrary reaction orders different from one (a, fi ^ 1), and Eqs. 912 for the 3D description of the heterogeneously catalyzed reaction), is the system of non-linear algebraic equations for the dependent variables at each grid point. The root-finding algorithm based on Newton's method was applied to develop an iterative numerical procedure for solving such complex non-linear systems.
However, in case of 3D description of the heterogeneous reaction-diffusion process (Eqs. 9-12) in the form
Mx + R(x) = 0, where M is matrix and R(x) is rational function, at least 10 times bigger systems have to be solved and the application of the Newton's method was found less efficient in comparison with 2D problems. Therefore, the iterative procedure in the form
(1ó)
was applied which enables fast converging to the solution. Mathematica codes were developed to solve the complex systems of model equations.
a)
c)
b)
d)
Figure 5. a) - Steady-state dimensionless concentration profile of the component A in the microchannel without a chemical reaction (vav=6.6 x 10- ms-1; Da=2 x 10-10 m2s-1; DB = 5 x 10-10 m2s-1) with the outlet profiles of the same component A at different non-equidistant mesh resolutions (b); the influence of the density of non-equidistant grid points on the dimensionless profiles of the component B at higher velocity (vav = 13.2 x 10-3 ms-1) at the 3/4 of the microchannel length (c), and at the outlet (d).
3. Results and Discussion
Consider first that the reaction follows general 2nd order kinetics (-rA) = kcAcB and therefore the reaction orders in the reaction terms (Eqs. 3-6, and 8) are equal to one (a = fi = 1). In that simplified case the non-linear problem can be solved by the iterations of the linear systems. In addition, the solutions of the mass balance equations in terms of relative concentration (cA ou/cA 0) for the ideal mixed-flow reactor and continuous plug-flow reactor are required to compare the performance of different reactor types with the microreactor.
Considering a second-order irreversible homogeneous reaction, a dimensionless concentration of component A at steady-state conditions is given by
(17)
- / + 71 + 2 Da + Per (- / + M)' + Da(- l+M)
2 Da M
-e°"(2-M)^ ^- e0"11 (2- Mp^i M
(18)
for the ideal mixed-flow reactor, XCSTR, and for the plug-flow reactor, XPFR. For M = cB(/cA0 = 1, the above equations (Eqs. 17, 18) are simplified to
(19)
(20)
where Da = k cB0 L/vav. Using a wide range of theoretical operating conditions (0.01 < Pe < 10), dimensionless outlet concentrations as a function of Da were calculated for both types of ideal macro-reactors (Eqs. 19, 20) and the microchannel (Eqs. 3-8, where a = ¡3 = 1) and graphically compared in Figure 6. Generally, an increase in the flow rate results in an increase in Pe, and consequently, our model for the microreactor approaches that for a PFR. At 0.5 < Pe < 10, the model gives close-to-ideal PFR performance. As expected, a further increase in Pe number (Pe > 10) results in the shifting away of the model simulations with the PFR curve. Namely, higher values of Pe mean higher velocities and smaller diffusion coefficients, which is reflected in shorter residence and reactions times. Only in the case of very low velocities, expressed as Pe values below 0.1, the microreactor does shift away from the CSTR behavior (Figure 6). This phenomenon might be explained by the fact that at such low velocities and higher values of diffusion coefficients (e.g., Pe = 0.01 at vmax = 0.0001 ms-1 and DA= 6.6 x 10-7 m2s-1) the diffusion of the components in the flow direction (x-direction) becomes more significant.
The numerical results of the convection-diffusion-reaction dynamics in the microchannel for the selected geometry (width/length = 200^m/1cm) based on non-linear reaction kinetics are presented in Figure 7 for a wide range of arbitrary reaction rates orders a and ¡ .

Figure 6. The performances of the ideal PFR and CSTR classical macro-reactors in comparison with the microreactor for a wide range of theoretical operating conditions, expressed as different Pe numbers (0.01 < Pe < 10).
reactant B
Figure 7. Averaged dimensionless outlet concentrations of components A and B as the function of the variety of the reaction rates orders a and 3 (-rA = k X c^; k = 0.1 m3(a+3-1) kmol-(a+3-1)s-1; vmx=0.002 ms-1; DA = 1 X 10-9 m2s-1; DB = 2 x 0-9 m2s-1).
In order to assure the stability of the described iterative numerical procedure for solving non-linear systems, simple mathematical manipulation is recommended. Namely, the root-finding algorithm based on Newton's method requires the differentiable function fon the Euclidean «-space. Therefore, we define the power exponents in the reaction terms DajXaYp (Eqs. 3 and 5) as Da/X2)a/2(Y2)^/2. Because the power exponent is usually defined as ea lnX, the proposed mathematical manipulation assures the stability of the numerical iteration in cases of feasible appearances of very small but still negative values of dependent variables X and Y.
Some typical simulated graphical results of the enzyme reaction process between two immiscible fluids at steady state conditions (Eqs. 9-12) with the velocity profile of the two-phase fluid within the microreactor are presented in Fig. 8. The results of the numerical simulation of a parallel fluid flow in the microchannel with the position of the interface area in the middle of the channel (Fig. 8.a) revealed that at steady-state conditions a fully developed profile takes place very shortly after the beginning of the microchannel, which was shown also in our previous work on steady-state two-phase extraction in the micro-reactor.39 The real values of model parameters for the selected enzyme reaction-diffusion process were used in the
simulations. Molecular diffusion coefficients of acetic acid (B) and isoamyl alcohol (A) in water (w) and in n-he-xane (h) were estimated by the Scheibel empirical correlation (DA^h=5.33 x 10-9m2/s, DB/h = 9.85 x 1Q-9m2/s, DB/w = 1.29 x 10-9m2/s), and the partitioning coefficient KP for acetic acid in ternary system with water and n-hexane was found to be 0.0167.45 The kinetic model parameters were estimated according to the literature and are summarized in Table 1.40
Table 1. Chosen binding and inhibition constants and kinetic parameters used in simulations.
vr,max'103 (mol/m3 s)	kb (mol/m3)	ka (mol/m3)	Ki,B (mol/m3)
0.32	0.26-10-3	0.96-10-3	2.00-10-3
The results of numerical simulation of reactant B dimensionless concentration profile in water phase YBW and organic phase YB/h at defined conditions with equimolar starting reactants concentration 0.5 x 10-3 mol/m3 are presented in Fig. 8b and c, respectively. As expected due to the low Kp value of B in organic-water system, the component B readily diffused from the organic phase into the
b) d)
Figure 8. (a) A 3D presentation of a flow pattern of both phases in a microchannel at fw = 10 pl/min and fh = 34 pl/min using Navier-Stokes equations. Graphical presentation of results of a numerical simulation for YB concentration profile along the microchannel (b) in organic and (c) in aqueous phase at z = 0.5H/W. (d) 3D graphical isosurface presentation of concentration of reactant B YB/w in organic phase.
aqueous phase. As seen from Fig. 8b, it was almost completely removed from the organic phase and consequently, the concentration of B in the aqueous phase increased at the beginning of the microchannel and reached 3-times higher values than in organic phase due to approximately 3-fold slower flow rate of aqueous phase regarding the organic phase, while further decrease in acid concentration was due to the reaction at the two-liquid boundary. 3D graphical isosurface presentation of dimensionless concentration of reactant B in organic phase YB/h is presented in Figure 8.d where the evident influence of the velocity distribution on the concentration profile can be observed. The process optimization and the determination of real kinetic parameters are feasible by the developed model and experimental work.27
4.	Conclusions
The theoretical descriptions and analysis of the convection-diffusion-reaction processes for the irreversible homogeneous and non-homogeneous catalytic reactions in a microreactor were investigated. Numerical experiments of a 2D model, developed for the specific case without a chemical reaction taking place in the microchannel, accurately predict the reasonably expected outlet concentrations of both components. The microreactor model simulations were compared with the ideal plug-flow reactor and continuous stirred-tank reactor predictions for a broad range of operating conditions in order to asses under what conditions (microchannel geometry, diffusion, convection, chemical kinetics) a given microreactor is more efficient or productive than a classical macro-reactor. The comparison of the model simulations with the ideal reactor predictions effectively shows close to PFR behavior at 0.5 < Pe < 10. The proposed numerical procedure based on a relatively simple implicit finite difference technique which is improved by a non-equidistant partition of grid points assures accurate, stable, and theoretically consistent solutions of the non-linear systems. The developed and presented 3D mathematical model comprising of flow distribution, transport phenomena and enzyme reaction kinetics, could easily be applied to different reaction types in homogeneous and heterogeneous systems. Further mi-croreactor design and process optimization by means of the use of the optimized microchannel geometry, reac-tants/enzyme ratio, and the determination of real kinetic parameters are feasible by the developed model and additional experimental work.
5.	References
1.	K. F. Jensen, Chem. Eng. Sci., 2001, 56, 293-303.
2.	O. Worz, K. P. Jackel, T. Richter and A. Wolf, Chem. Eng.
Sci., 2001, 56, 1029-1033.
3.	K. Jähnisch, V. Hessel, H. Löwe and M. Baerns, Angew Chem Int Ed., 2004, 43, 406-446.
4.	D. M. Roberge, L. Ducry, N. Bieler, P. Cretton, B. Zimmermann, Chem. Eng. Techn., 2005, 28 (3), 318-323.
5 V. Hessel, H. Löwe, F. Schönfeld, Chem. Eng. Sci., 2005a, 60, 2479-2501.
6.	V. Hessel and H. Löwe, Chem. Eng. Techn., 2005, 28 (3), 267-284.
7.	P. Watts and S. J. Haswell, Chem. Eng. Techn., 2005, 28 (3), 290-301.
8.	P. L. Urban, D. M. Goodall and N. C. Bruce, Biotechnology Advances, 2006, 24, 42-57.
9.	K. Mae, Chem. Eng. Sci., 2007, 62, 4842-4851.
10.	A. Pohar and I. Plazl, Chem. Biochem. Eng. Q., 2009, 23 (4) 537-544.
11.	S. K. Ajmera, M. W. Losey, K. F. Jensen, M. A. Schmidt, AIChE J, 2001, 47, 1639-1647.
12.	K. S. Ajmera, C. Delattre, A. M. Schmidt, K. F. Jensen, Stud. Surf. Sci. Catal., 2003, 145, 97-102.
13.	J. Yoshida, A. Nagaki, T. Iwasaki, S. Suga, Chem. Eng. Techn., 2005, 28, 259-266.
14.	G. Wiessmeier and D. Honicke, J Micromech. Microeng., 1996, 6, 285-289.
15.	I. M. Hsing, R. Srinivasan, M. P. Harold, K. Jensen, M. A. Schmidt, Chem. Eng. Sci., 2000, 55, 3-18.
16.	J. C. Schouten, E. V. Rebrov, M. H. J. de Croon, Chimia, 2002, 56, 627-635.
17.	G. N. Jovanovic, P. Znidarsic-Plazl, P. Sakrittichai, K. Al-Khaldi, Ind. Eng. Chem. Res., 2005, 44,44, 5099-5106.
18.	S. Maehara, M. Taneda and K. Kusakabe, Chem. Eng. Res. Des., 2008, 86, 410-415.
19.	C. N. Baroud, F. Okkels, L. Ménétrier, P. Tabeling, Phys. Rev. E, 2003, 67, 60104.
20.	A. N. Waghode, N. S. Hanspal, I. M. T. A. Shigidi, V. Nasse-hi, K. Hellgardt, Chem. Eng. Sci., 2005, 60, 5862-5878.
21.	J. Atencia and D. J. Beebe, Nature, 2005, 437, 648-655.
22.	N. Aoki, S. Hasebe and K. Mae, AIChE J., 2006, 52, 1502-1515.
23.	D. Bothea, C. Stemichb and H. J. Warneckeb, Chem. Eng. Sci., 2006, 61, 2479-2501.
24.	M. V. Kothare, Computers & Chemical Engineering, 2006, 30, 1725-1734.
25.	R. Kieffer, C. Charcosset, F. Puel, D. Mangin, Computers & Chemical Engineering, 2008, 32, 1325-1333.
26.	M. Tisma, B. Zelic, D. Vasic-Racki, P. Znidarsic-Plazl and I. Plazl, Chem. Eng. J, 2009, 149, 383-388.
27.	P. Znidarsic-Plazl and I. Plazl, Proc. Biochem., 2009, 44, 1115-1121.
28.	A. Pohar, P. Znidarsic-Plazl and I. Plazl, Lab Chip, 2009, 9, 23, 3385-3390.
29.	J. Koo and C. Kleinstreuer, Journal of Micromechanics and Microengineering, 2003, 13, 568-579.
30.	S. V. Gokhale, R. K. Tayal, V. K. Jayaraman, B. D. Kulkarni, International Journal of Chemical Reactor Engineering, 2005, 3, R2.
31.	T. Bayraktar and S. B. Pidugu, International Journal of Heat
and Mass Transfer, 2006, 49, 815-824.
32.	D. Jie, X. Diao, K. B. Cheong, L. K. Yong, Journal of Micro-mechanics and Microengineering, 2000, 3, 372-279.
33.	K. A. Alfadhel and M. V. Kothare, Chem. Eng. Sci., 2005, 60, 2911-2926.
34.	A. Pohar and I. Plazl, Ind. eng. chem. res., 2008, 47, 74477455.
35.	E. L. Paul, V. A. Atiemo-Obeng and S. M. Kresta (Ed.): Handbook of Industrial Mixing: Science and Practice, John Wiley and Sons, Hoboken, New Jersey, 2004, pp. 1187-1246.
36.	R. B. Bird, W. E. Stewart and E. N. Lightfoot (Eds.): Transport Phenomena, 2nd ed., John Wiley & Sons, New York, 2001.
37.	R. C. Aguirre and H. J. Catrakis, Physics of fluids, 2005, 17, 038103.
38.	J. P. G. Villaluenga and Y. Cohen, Journal of Membrane Science, 2005, 260, 119-130.
39.	P. Znidarsic-Plazl and I. Plazl, Lab Chip, 2007, 7, 883-889.
40.	M. D. Romero, L. Calvo, C. Alba, A. Daneshfar, Journal of Biotechnology, 2007, 127, 269-277.
41.	M. J. Berger and P. Collela, J. Comput. Phy., 1989, 82, 64-84.
42.	P. B. Bochev, G. Liao, G. C. de la Pena, Numerical Methods for Partial Differential Equations, 1996, 12, 489-506.
43.	F. Liu, S. Ji and G. Liao, SIAM Journal on Scientific Computing, 1998, 20, 811-825.
44.	X. X. Cai, D. Fleitas, B. Jiang, G. Liao, Computers and Mathematics with Applications, 2004, 48, 1077-1085.
45.	R. H. Perry, D. W. Green and J. O. Maloney (Eds.): Perry's Chemical Engineer's Handbook, 7th ed., McGraw-Hill Companies Inc., New York, 1999.
Povzetek
Delo predstavlja teoretičen opis z numeričnimi eksperimenti in analizo reakcijsko difuzijskih procesov homogenih in nehomogenih reakcij v mikroreaktorju. Razvit matematični model v brezdimenzijski obliki vključuje tudi hitrostni profil za laminarni tok mešljivih in nemešljivih tekočin v mikrokanalu pri stacionarnih pogojih. Rezultati napovedi obnašanja mikroreaktorja so primerjani z izračuni v klasičnem kontinuirnem mešalnem in pretočnem reaktorju za primer nepovratne reakcije dveh vstopnih komponent. Primerjava učinkovitosti različnih tipov reaktorjev je podana za široko območje obratovalnih pogojev, izraženih z Pe številom (0.01 < Pe < 10). Za reševanje kompleksnih nelinearnih sistemov je predlagana metoda končnih razlik, izboljšana z uvedbo neekvidistantnih razlik.