The simulation of MEMS (Micro-Electro-Mechanical-System) containing fluid field could not be well performed by conventional numerical analysis methods. The micro flow field characteristics can be simulated by using macromodel including a nonlinear analysis. This paper set up the macromodel of the micromixer of the microfluidic chip using Krylov subspace projection method. The system functions were assembled through finite element analysis using COMSOL. We took the flow field-concentration field analysis for micromixer finite element model. The finite element functions order is reduced by second-order Krylov subspace projection method based on Lanczos algorithm. It can be shown that the simulation results obtained by using the macromodel are highly consistent with the results of finite element analysis. The calculation using the macromodel is two orders of magnitude faster than the calculation performed by the finite element analysis method. This macromodel should facilitate the design of microfluidic devices with sophisticated channel networks

Hongxiang Wang

Jingliang Li

Xueye Chen

English

HRČAK - Portal of Croatian Scientific and Professional Journals

X. Chen, H. Zeng, H. Wang, J. Li: Nonlinear macromodel based on Krylov subspace for micromixer of the microfluidic chipENGINEERING MODELLING 25 (2012) 1-4, 45-49 45SUMMARYThe simulation of MEMS (Micro-Electro-Mechanical-System) containing fluid field could not be well performedby conventional numerical analysis methods. The micro flow field characteristics can be simulated by usingmacromodel including a nonlinear analysis. This paper set up the macromodel of the micromixer of the microfluidicchip using Krylov subspace projection method. The system functions were assembled through finite element analysisusing COMSOL. We took the flow field-concentration field analysis for micromixer finite element model. The finiteelement functions order is reduced by second-order Krylov subspace projection method based on Lanczos algorithm.It can be shown that the simulation results obtained by using the macromodel are highly consistent with the results offinite element analysis. The calculation using the macromodel is two orders of magnitude faster than the calculationperformed by the finite element analysis method. This macromodel should facilitate the design of microfluidicdevices with sophisticated channel networks.Key words: macromodel, Krylov subspace, micromixer, MEMS.UDC 517.9:532.5:519.61Original scientific paperReceived: 12.12.2012.Nonlinear macromodel based onKrylov subspace for micromixer ofthe microfluidic chipXueye Chen(1), Hong Zeng(1), Hongxiang Wang(1) and Jingliang Li(2)(1)Faculty of Mechanical Engineering and Automation, Liaoning University of Technology, Jinzhou 121001, CHINAe-mail: xueye_chen@126.com(2)The Key Laboratory for Micro/Nano Technology and System of Liaoning Province,Dalian University of Technology, Dalian 116024, CHINA1. INTRODUCTIONConsiderable effort has been put to the developmentof microfluidic systems for performing biological andbiochemical mixing [1, 2]. Micromixers, an importantbranch of MEMS (Micro-Electro-Mechanical-System),have as complex a structure and physics field as otherMEMS. The development of technologies in the fieldof Microsystems has generated a wideswpread interestin exploring fast and accurate design methods tosimulate entire systems. The typical MEMS cancomplete many tasks under certain conditions. Multi-physics fields in MEMS make it difficult to designstructure especially if these contain a flow field. InMEMS, the coupled effect commonly takes placeamong mechanical, electrical, fluidic, magnetic andthermal fields [3-6]. The strongly nonlinear effect isusually intractable [7]. In order to design MEMS, finiteelement and finite volume methods have beensuccessfully applied. These numerical methods arehighly effective and accurate but laborious and timeconsuming, i.e. a solution of three-dimensions Jouleheating dispersion by using the finite element methodcan cost four days and four GB of physical memory[8]. Such computational cost is prohibitive for system-level design of complex microfluidic chips system.Furthermore, when the designed system containsconversion of different signals, such as of non-electricsignals to electric signals or other signal processingunits, the finite element or finite volume methods arepowerless. In order to solve these mentioned problems,the macromodel technologies have been widelyinvestigated in the past few years [9, 10].The macromodel, also called the reduced ordermodel, is a low-order behavioural representation of adevice. The macromodel can be described withX. Chen, H. Zeng, H. Wang, J. Li: Nonlinear macromodel based on Krylov subspace for micromixer of the microfluidic chip46 ENGINEERING MODELLING 25 (2012) 1-4, 45-49hardware language and directly applied to EDA(Electronic Design Automation) environment.Analytical solution [11] and numerical solution [12] aretwo methods for creating the macromodel. The matrixsubspace projection based on matrix transformation isa common numerical method for extractingmacromodel from the finite element analysis.In this work, Lanczos algorithm is applied in settingup the second-order Krylov subspace. We have directlyprojected a large-scale sparse matrix to the subspace.The algorithm has been tested on the model of themicromixer of the microfluidic chip. The resultsobtained by the macromodel, when compared with theresults obtained by applying the finite element methodto the micromixer of the microfluidic chip, shows thatthe model is completly accurate.2. THEORETICALY BACKGROUNDThe static and dynamic behaviours of MEMS ispresented by partial differential equations (PDEs) andcorresponding boundary conditions. Discretizationusing finite element method, finite volume method orfinite difference method can transfer PDEs to ordinarydifferential equations (ODEs), which can be written inthe following state-space form:( t ) ( t ) ( t ):( t ) ( t )= ⋅ + ⋅⎧⎨ = ⋅⎩x A x B uΩy C x (1)where, A∈Rm×m is the system matrix, B∈Rm×n is inputmatrix, C∈Rq×m is output matrix, u∈Rn is inputvariable. y∈Rq is output variable and x∈Rm is statevariable. We replace x=V⋅z, V∈Rm×r, z∈Rr. Thedimensionless column vectors of transition matrix V,whose number is r, can compose a basis subspace,called projection subspace. If r<<m, then the r-dimensional vector z can be considered as thereduced order state variable. The reduced order statefunction is:ˆ ˆ( t ) ( t ) ( t )ˆˆ ( t ) ( t )⎧ = ⋅ + ⋅⎪⎨= ⋅⎪⎩z A z B uy C z (2)where, -1ˆ =A V AV , -1ˆ =Β V B , ˆ =C CV .Performing a Laplace transform for Eq. (1), we canobtain its frequency domain transfer function H(s) by:1 1 1( s ) ( s ) ( s )− − −= ⋅ − ⋅ = − ⋅ −H C I B C A I A BΑ (3)Finding a suitable transition matrix V is the key ofmatrix subspace projection method for the reductionof the original system order. In order to apply Eq. (2)to reduce the order of Eq. (1), ˆmin −y y  andˆmin ( s ) ( s )−H H  in time and frequency domain,respectively, must be satisfied.An r-dimension Krylov subspace ϕr is defined asfollows:{ }r 1r 1 1 1 1 1 1 1( , ) span , , −=φ A b b A b A b (4)where, A1∈Rm×m, b1∈Rm is called starting vector.Lanczos algorithm [13] can convert large linearequations to tri-diagonal equations and create theorthonormal vectors, makig it possible to construct theunit basic Krylov subspace. The Lanczos process isrepresented with the chart given in Figure 1.Fig. 1  The flow chart of Lanczos algorithmThe tri-diagonal process is in fact a similaritytransformation, and can be expressed by:-1n n n=V AV T (5)We can simplify Eq. (3) to:1 11 -1n n nT T 11 n 1( s ) ( s )( ) ( s )i ( s ) iΙ− −−−= − ⋅ − == ⋅ − ⋅ − == ⋅ ⋅ ⋅ − ⋅H C A A BC A I V T V BC A B I T (6)where i∈Rn is the first standard orthogonal basis inorthogonal basis set Ti jw v 0, i j⋅ = ≠ . So, we canobtain the moment of the transition function of initialsystem, 1 T ii 1 n 1m i i−= − ⋅ ⋅A B T . The transition functionof the reduced system becomes:T T 1r 1 r 1( s ) i ( s ) i− −= − ⋅H A B I T (7)where Tr, is the first r×r sub-matrix of Tn, and i1∈Rris the first standard orthogonal basis.If the transition matrix V, obtained from Lanczosprocess is a unit basic Kyle subspace [14], and thefirst r moments between H(s) and ˆ ( s )H  are matchingthen we can obtain the minimum of s sˆ( ) ( )−H H .Therefore, the original system function matrix can beprojected to Krylov subspace and the order can bereduced.X. Chen, H. Zeng, H. Wang, J. Li: Nonlinear macromodel based on Krylov subspace for micromixer of the microfluidic chipENGINEERING MODELLING 25 (2012) 1-4, 45-49 47In many engineering problems, second orderdifferential functions are used for describing thedynamical behaviour of a system, as below:T( t ) ( t ) ( t ) ( t )( t ) ( t )+ + =⎧⎪⎨=⎪⎩Mx Ex Kx Buy C x (8)where, M, E, K∈Rm×m, B∈Rm×n, C∈Rq×m, u∈Rn ,y∈Rq and x∈Rm. The second order Krylov subspacecan be defined by:{ }r 1 2 1 0 1 r 1( , , ) span , , −=φ A A b k k k (9)where, k0=b1, k1=A1b1, ki=A1ki-1+A2ki-2, andA1A2∈Rm×m, b1∈Rm. b1 is called the starting vectorand ki is called the basic vector. The preconditionedLanczos process [15] can be used for creating thesecond order Krylov subspace, which is, in turn, usedto project the dynamic behaviour of a system. As forthe projection x=Gxr, G∈Rm×r, r<<m, if oneP∈Rm×r=G, the first r moments of the original secondorder system and the reduced order model can matchand we obtain the following reduced system:T T T Tr r rTr( t ) ( t ) ( t ) ( t )⎧ + =⎪⎨⎪ =⎩P MGx P EGx P KGx P Buy C Gx(10)3. MICROMIXER FINITE ELEMENTMODELMicromixer is one of the important components inMicro Total Analysis System and Lab on a Chip. Rapidmixing is essential in many of the microfluidic systemsused in biochemical analysis, drug delivery andsequencing or synthesis of nucleic acids. Theinvestigation of micromixers is fundamental tounderstanding transport phenomena on the microscale.The micromixer fabricated in our laboratory is shown inFigure 2. The outside measurements of the chip are 50mm (lenght) × 25 mm (width). The width of the micro-channel is 300 microns and the depth is 30 microns. Themicro-channel was fabricated from a PMMA(polymethyl methacrylate) substrate by hot embossing,and sealed with an other PMMA cover plate enclosingthe microchannel. The micromixer can perform mixingby utilizing the “split and recombine” principle. Sincethe thickness of the microchanel is much smaller thanlength and width, we consider it as a plane problem andmesh it with 2D element.The finite element micromixer model is coupledwith mechanical, fluidic and concentration fields,therefore, it is a strongly nonlinear problem. The finiteelement analysis contains three steps. Firstly, theincompressible Navier-Stokes equation, Eq. (11), andcontinuity equation, Eq. (12), are solved in order to getthe velocity field:1ftµ∆ρ∂+ ⋅∇ = − ∇ +∂v v v P v (11)0∇ ⋅ =v (12)where, ν is the velocity vector, f is the body force, ρ isthe density of the fluid, P is the pressure and µ is thedynamic viscosity of the fluid. Secondly, theconcentration distribution can be calculated fromconvection–diffusion equation, Eq. (13), after obtainingthe velocity field:2( ) Dtφ φ φ∂ + ⋅∇ = ∇∂v (13)where, φ is the concentration of species, D is thediffusion coefficient and ν is the fluid velocity. Then,the concentration can be used for investigating thedegree of fluid mixing by equation:W0W00- dy( x ) 1 - 100%- dyφ φγφ φ∞∞= ×∫∫(14)where, W is the width of the micro-channel, φ is theconcentration of species across the width of the mixingchannel, φ0 and φ∞ are the sample concentrationprofiles with completely unmixed (0 or 1) andcompletely mixed states (0.5), respectively.The finite element model of the micromixer is setup using COMSOL software. The mesh model isshown in Figure 3. The number of elements and nodesare 3936 and 5095 respectively.Fig. 2  The structure of the micromixerFig. 3  The finite element model of the micromixerWe have solved the coupling problem by usingCOMSOL's sequential coupling. The coupling in thematrix equation is shown as:U1U200 0 0φφ φ φ⎡ ⎤⎡ ⎤ ⎡ ⎤⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤+ + =⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦⎣ ⎦ ⎣ ⎦D 0 FM u u K K uFφ φ φD D K(15)where, u and ϕ are velocity and concentration degreesof freedom, respectively. The system matrices M, K,Dϕ and Kϕ are assembled following a standard finiteelement methodology. Matrices DUϕ and KUϕ representthe coupled effect attribute. COMSOL has solved thesenonlinear equations using CG solver based on AMG.X. Chen, H. Zeng, H. Wang, J. Li: Nonlinear macromodel based on Krylov subspace for micromixer of the microfluidic chip48 ENGINEERING MODELLING 25 (2012) 1-4, 45-494. MACROMODEL EXTRACTION FORMICROMIXERThe COMSOL software was used to conduct finiteelement simulation for the micromixer and somecharacteristic parameters are set up and presented inTable 1.Table 1. Finite element simulation parameters5. CONCLUSIONSA macromodel is presented and successfullycalculated for the microfluidic chip micromixer. It hasbeen proved that the macromodel based on Kylesubspace is an effective method for the MEMSsimulation. The strongly coupled nonlinear microfluidicproblems can be solved by the use of the macromodelmethod in combination with the finite element analysis.In the future, the macromodel method based on Kylesubspace will have a great potential in the design andsimulation of microfluidic analysis systems.6. ACKNOWLEDGMENTSThe authors are grateful to Liaoning Province KeyLaboratory Funded Projects (No: LS2010080,LKF2011020), and National Science TechnologySupport Plan Projects (No: 2012BAF12B08-5).7. REFERENCES[1] J. Cao, P. Cheng and F.J. Hong, A numerical studyof an electrothermal vortex enhanced micromixer,Microfluidics and Nanofluidics, Vol. 5, No. 1, pp.13-21, 2008.[2] C.A. Cortes-Quiroz, M. Zangeneh and A. Goto,On multi-objective optimization of geometry ofstaggered herringbone micromixer, Microfluidicsand Nanofluidics,  Vol. 7, No. 1, pp. 29-43, 2009.[3] W. Dai and Y. P. Zhao, The nonlinear phenomenaof thin polydimethylsiloxane (PDMS) films inelectrowetting, International Journal ofNonlinear Sciences and Numerical SimulationVol. 8, No. 4, pp. 519-526, 2007.[4] O.-T. Son, J.W. Roh, S.-H. Song, J.-S. Park, W.Lee and H.-I. Jung, Continuous micro-magnetophoretic separation using a dipolemagnetic field, BioChip Journal, Vol. 2, No. 3,pp. 186-191, 2008.[5] Y.-Q. Wan, Q. Guo and N. Pan, Thermo-electro-hydrodynamic model for electrospinning process,International Journal of Nonlinear Sciences andNumerical Simulation, Vol. 5, No. 1, pp. 5-8,2004.[6] Y.C. Zeng, Y. Wu, Z.G. Pei and C.W. Yu,Numerical approach to electrospinning,International Journal of Nonlinear Sciences andNumerical Simulation,  Vol. 7, No. 4, pp. 385-388, 2006.[7] D.H. Shou and J.H. He, Application of parameter-expanding method to strongly nonlinearoscillators, International Journal of NonlinearSciences and Numerical Simulation,  Vol. 8, No.1, pp. 121-124, 2007.Simulation parameter Value Density of fluid Dynamic viscosity o f the fluid Diffusion coefficient Channel width Wall boundary 1 input concentration 2 input concentration 103 (kg/m3) 1.002×10-3 (Pa·s) 10-10 (m2/s) 3×10-4 (µm) No slip 0.01 (mol/m3) 0.1 (mol/m3) The command “assemble” can extract stiffness matrixand other system matrices from COMSOL finite elementmodel analysis outputs, and the command code is:[ , , , ] assemble( fem );spy( )=K L M N K (16)Using the system matrices extracted from the finiteelement model, Lanczos algorithm is executed inMATLAB language for the steady-state analysis of themicromixer. We reduced the original model to 40 and 20order in MATLAB language and we also calculated themacromodel output value. The calculation times ofCOMSOL, 40 order macromodel and 20 ordermacromodel are 2.2 s, 0.027 s and 0.021 s, respectively.The macromodel is two orders of magnitude faster thanthe finite element model. The micromixer outputconcentration results from COMSOL, 40 and 20 ordermacromodel are shown in Figure 4. As it can be seenfrom Figure 4, the calculation error from 40 order modelis less 5% less than the one obtained by COMSOL,therefore, it can be concluded that the calculation errordecreases along with the order number increment.Furthermore, it means that the compromise betweencalculation precision and speed should be taken intoconsideration. 0.0 0.2 0.4 0.6 0.8 1.00.3400.3420.3440.3460.3480.3500.352ConcentrationChannel coordinate Coordinates Coordinates CoordinatesFig. 4  COMSOL model, 40 order macromodel and 20 ordermacromodel output concentration resultsX. Chen, H. Zeng, H. Wang, J. Li: Nonlinear macromodel based on Krylov subspace for micromixer of the microfluidic chipENGINEERING MODELLING 25 (2012) 1-4, 45-49 49[8] Y. Wang, Q. Lin and T. Mukherjee, A model forJoule heating-induced dispersion in microchipelectrophoresis, Lab on a Chip,  Vol. 4, No. 6,pp. 625-631, 2004.[9] D. Vasilyev, M. Rewienski and J. White,Macromodel generation for bioMEMScomponents using a stabilized balanced truncationplus trajectory piecewise-linear approach, IEEETransactions on Computer-Aided Design ofIntegrated Circuits and Systems, Vol. 25, No. 2,pp. 285-293, 2006.[10] H.M. Park and J.Y. Lim, A reduced-order modelof the low-voltage cascade electroosmoticmicropump, Microfluidics and Nanofluidics,Vol. 6, No. 4, pp. 509-520, 2009.[11] R. Qiao and N.R. Aluru, A compact model forelectroosmotic flows in microfluidic devices,Journal of Micromechanics and Micro-engineering,  Vol. 12, No. 5, pp. 625-635, 2002.[12] X. He and S. Hauan, Microfluidic modeling anddesign for continuous flow in electrokineticmixing-reaction channels, Aiche Journal, Vol. 52,No. 11, pp. 3842-3851, 2006.[13] M. Troyer, H. Tsunetsugu and D. Wurtz, Thermo-dynamics and spin gap of the Heisenberg laddercalculated by the look-ahead Lanczos-algorithm,Physical Review B,  Vol. 50, No. 18, pp. 13515-13527, 1994.[14] S. Sundar and B.K. Bhagavan, Comparison ofKrylov subspace methods with preconditioningtechniques for solving boundary value problems,Computers & Mathematics with Applications,Vo. 38, No.11-12, pp. 197-206, 1999.[15] T.J. Frankcombe and S.C. Smith, Fast, scalablemaster equation solution algorithms, IV. Lanczositeration with diffusion approximationpreconditioned iterative inversion, Journal ofChemical Physics, Vol. 119, No. 24, pp. 12741-12748, 2003.PRIKAZ NELINEARNOG MAKROMODELA TEMELJENOG NA KRYLOVLJEVOMPOTPROSTORU PRIKLADNOG ZA SIMULACIJU RADA MIKROMIKSERAMIKROFLUIDNOG PROCESORASA@ETAKU ovome se radu polazi od pretpostavke da konvencionalne metode numeri~ke analize nisu dostatne za simulacijumikro-elektro-mehani~kog-sustava (MEMS) koji sadr`i polje fluida. Me|utim, korištenjem makromodela s primjenomnelinearne analize mogu se simulirati karakteristike mikro-polja toka. U ovome su radu prikazane postavkemakromodela mikromiksera mikrofluidnog procesora pri ~emu je korištena Krylovljeva metoda projiciranjapotprostora. Funkcije sustava dobivene su metodom kona~nih elemenata uz korištenje softvera COMSOL. Zadiskretizirani model mikromiksera uzeta je analiza polja koncentracije. Red funkcije kona~nih elemenata reduciranje pomo}u Krylovljeve metode projekcije potprostora drugog reda bazirane na Lanczosovom algoritmu. Pokazano jeda se rezultati dobiveni korištenjem makromodela podudaraju s rezultatima dobivenim analizom s primjenom kona~nihelemenata, s tim da je vrijeme trajanja simulacije pri korištenju makromodela stotinu puta kra}e od vremena trajanjasimulacije korištenjem metode kona~nih elemenata. U prakti~noj primjeni, ovaj bi makromodel trebao olakšatidizajniranje mikrofluidnih naprava sa sofisticiranom mre`om kanala.Klju~ne rije~i: makromodel, Krylovljev potprostor, mikromikser, mikro-elektro-mehani~ki sustavi (MEMS).

Nonlinear macromodel based on Krylov subspace for micromixer of the microfluidic chip

Hrčak - Portal of scientific journals of Croatia

X. Chen, H. Zeng, H. Wang, J. Li: Nonlinear macromodel based on Krylov subspace for micromixer of the microfluidic chipENGINEERING MODELLING 25 (2012) 1-4, 45-49 45SUMMARYThe simulation of MEMS (Micro-Electro-Mechanical-System) containing fluid field could not be well performedby conventional numerical analysis methods. The micro flow field characteristics can be simulated by usingmacromodel including a nonlinear analysis. This paper set up the macromodel of the micromixer of the microfluidicchip using Krylov subspace projection method. The system functions were assembled through finite element analysisusing COMSOL. We took the flow field-concentration field analysis for micromixer finite element model. The finiteelement functions order is reduced by second-order Krylov subspace projection method based on Lanczos algorithm.It can be shown that the simulation results obtained by using the macromodel are highly consistent with the results offinite element analysis. The calculation using the macromodel is two orders of magnitude faster than the calculationperformed by the finite element analysis method. This macromodel should facilitate the design of microfluidicdevices with sophisticated channel networks.Key words: macromodel, Krylov subspace, micromixer, MEMS.UDC 517.9:532.5:519.61Original scientific paperReceived: 12.12.2012.Nonlinear macromodel based onKrylov subspace for micromixer ofthe microfluidic chipXueye Chen(1), Hong Zeng(1), Hongxiang Wang(1) and Jingliang Li(2)(1)Faculty of Mechanical Engineering and Automation, Liaoning University of Technology, Jinzhou 121001, CHINAe-mail: xueye_chen@126.com(2)The Key Laboratory for Micro/Nano Technology and System of Liaoning Province,Dalian University of Technology, Dalian 116024, CHINA1. INTRODUCTIONConsiderable effort has been put to the developmentof microfluidic systems for performing biological andbiochemical mixing [1, 2]. Micromixers, an importantbranch of MEMS (Micro-Electro-Mechanical-System),have as complex a structure and physics field as otherMEMS. The development of technologies in the fieldof Microsystems has generated a wideswpread interestin exploring fast and accurate design methods tosimulate entire systems. The typical MEMS cancomplete many tasks under certain conditions. Multi-physics fields in MEMS make it difficult to designstructure especially if these contain a flow field. InMEMS, the coupled effect commonly takes placeamong mechanical, electrical, fluidic, magnetic andthermal fields [3-6]. The strongly nonlinear effect isusually intractable [7]. In order to design MEMS, finiteelement and finite volume methods have beensuccessfully applied. These numerical methods arehighly effective and accurate but laborious and timeconsuming, i.e. a solution of three-dimensions Jouleheating dispersion by using the finite element methodcan cost four days and four GB of physical memory[8]. Such computational cost is prohibitive for system-level design of complex microfluidic chips system.Furthermore, when the designed system containsconversion of different signals, such as of non-electricsignals to electric signals or other signal processingunits, the finite element or finite volume methods arepowerless. In order to solve these mentioned problems,the macromodel technologies have been widelyinvestigated in the past few years [9, 10].The macromodel, also called the reduced ordermodel, is a low-order behavioural representation of adevice. The macromodel can be described withX. Chen, H. Zeng, H. Wang, J. Li: Nonlinear macromodel based on Krylov subspace for micromixer of the microfluidic chip46 ENGINEERING MODELLING 25 (2012) 1-4, 45-49hardware language and directly applied to EDA(Electronic Design Automation) environment.Analytical solution [11] and numerical solution [12] aretwo methods for creating the macromodel. The matrixsubspace projection based on matrix transformation isa common numerical method for extractingmacromodel from the finite element analysis.In this work, Lanczos algorithm is applied in settingup the second-order Krylov subspace. We have directlyprojected a large-scale sparse matrix to the subspace.The algorithm has been tested on the model of themicromixer of the microfluidic chip. The resultsobtained by the macromodel, when compared with theresults obtained by applying the finite element methodto the micromixer of the microfluidic chip, shows thatthe model is completly accurate.2. THEORETICALY BACKGROUNDThe static and dynamic behaviours of MEMS ispresented by partial differential equations (PDEs) andcorresponding boundary conditions. Discretizationusing finite element method, finite volume method orfinite difference method can transfer PDEs to ordinarydifferential equations (ODEs), which can be written inthe following state-space form:( t ) ( t ) ( t ):( t ) ( t )= ⋅ + ⋅⎧⎨ = ⋅⎩x A x B uΩy C x(1)where, A∈Rm×m is the system matrix, B∈Rm×n is inputmatrix, C∈Rq×m is output matrix, u∈Rn is inputvariable. y∈Rq is output variable and x∈Rm is statevariable. We replace x=V⋅z, V∈Rm×r, z∈Rr. Thedimensionless column vectors of transition matrix V,whose number is r, can compose a basis subspace,called projection subspace. If r<<m, then the r-dimensional vector z can be considered as thereduced order state variable. The reduced order statefunction is:ˆ ˆ( t ) ( t ) ( t )ˆˆ ( t ) ( t )⎧ = ⋅ + ⋅⎪⎨ = ⋅⎪⎩z A z B uy C z(2)where, -1ˆ =A V AV , -1ˆ =Β V B , ˆ =C CV .Performing a Laplace transform for Eq. (1), we canobtain its frequency domain transfer function H(s) by:1 1 1( s ) ( s ) ( s )− − −= ⋅ − ⋅ = − ⋅ −H C I B C A I A BΑ (3)Finding a suitable transition matrix V is the key ofmatrix subspace projection method for the reductionof the original system order. In order to apply Eq. (2)to reduce the order of Eq. (1), ˆmin −y y  andˆmin ( s ) ( s )−H H  in time and frequency domain,respectively, must be satisfied.An r-dimension Krylov subspace ϕr is defined asfollows:{ }r 1r 1 1 1 1 1 1 1( , ) span , , −=φ A b b A b A b" (4)where, A1∈Rm×m, b1∈Rm is called starting vector.Lanczos algorithm [13] can convert large linearequations to tri-diagonal equations and create theorthonormal vectors, makig it possible to construct theunit basic Krylov subspace. The Lanczos process isrepresented with the chart given in Figure 1.Fig. 1  The flow chart of Lanczos algorithmThe tri-diagonal process is in fact a similaritytransformation, and can be expressed by:-1n n n=V AV T (5)We can simplify Eq. (3) to:1 11 -1n n nT T 11 n 1( s ) ( s )( ) ( s )i ( s ) iΙ− −−−= − ⋅ − == ⋅ − ⋅ − == ⋅ ⋅ ⋅ − ⋅H C A A BC A I V T V BC A B I T (6)where i∈Rn is the first standard orthogonal basis inorthogonal basis set Ti jw v 0, i j⋅ = ≠ . So, we canobtain the moment of the transition function of initialsystem, 1 T ii 1 n 1m i i−= − ⋅ ⋅A B T . The transition functionof the reduced system becomes:T T 1r 1 r 1( s ) i ( s ) i− −= − ⋅H A B I T (7)where Tr, is the first r×r sub-matrix of Tn, and i1∈Rris the first standard orthogonal basis.If the transition matrix V, obtained from Lanczosprocess is a unit basic Kyle subspace [14], and thefirst r moments between H(s) and ˆ ( s )H  are matchingthen we can obtain the minimum of s sˆ( ) ( )−H H .Therefore, the original system function matrix can beprojected to Krylov subspace and the order can bereduced.X. Chen, H. Zeng, H. Wang, J. Li: Nonlinear macromodel based on Krylov subspace for micromixer of the microfluidic chipENGINEERING MODELLING 25 (2012) 1-4, 45-49 47In many engineering problems, second orderdifferential functions are used for describing thedynamical behaviour of a system, as below:T( t ) ( t ) ( t ) ( t )( t ) ( t )+ + =⎧⎪⎨ =⎪⎩Mx Ex Kx Buy C x (8)where, M, E, K∈Rm×m, B∈Rm×n, C∈Rq×m, u∈Rn ,y∈Rq and x∈Rm. The second order Krylov subspacecan be defined by:{ }r 1 2 1 0 1 r 1( , , ) span , , −=φ A A b k k k" (9)where, k0=b1, k1=A1b1, ki=A1ki-1+A2ki-2, andA1A2∈Rm×m, b1∈Rm. b1 is called the starting vectorand ki is called the basic vector. The preconditionedLanczos process [15] can be used for creating thesecond order Krylov subspace, which is, in turn, usedto project the dynamic behaviour of a system. As forthe projection x=Gxr, G∈Rm×r, r<<m, if oneP∈Rm×r=G, the first r moments of the original secondorder system and the reduced order model can matchand we obtain the following reduced system:T T T Tr r rTr( t ) ( t ) ( t ) ( t )⎧ + =⎪⎨⎪ =⎩P MGx P EGx P KGx P Buy C Gx (10)3. MICROMIXER FINITE ELEMENTMODELMicromixer is one of the important components inMicro Total Analysis System and Lab on a Chip. Rapidmixing is essential in many of the microfluidic systemsused in biochemical analysis, drug delivery andsequencing or synthesis of nucleic acids. Theinvestigation of micromixers is fundamental tounderstanding transport phenomena on the microscale.The micromixer fabricated in our laboratory is shown inFigure 2. The outside measurements of the chip are 50mm (lenght) × 25 mm (width). The width of the micro-channel is 300 microns and the depth is 30 microns. Themicro-channel was fabricated from a PMMA(polymethyl methacrylate) substrate by hot embossing,and sealed with an other PMMA cover plate enclosingthe microchannel. The micromixer can perform mixingby utilizing the “split and recombine” principle. Sincethe thickness of the microchanel is much smaller thanlength and width, we consider it as a plane problem andmesh it with 2D element.The finite element micromixer model is coupledwith mechanical, fluidic and concentration fields,therefore, it is a strongly nonlinear problem. The finiteelement analysis contains three steps. Firstly, theincompressible Navier-Stokes equation, Eq. (11), andcontinuity equation, Eq. (12), are solved in order to getthe velocity field:1ftµ∆ρ∂ + ⋅∇ = − ∇ +∂v v v P v (11)0∇ ⋅ =v (12)where, ν is the velocity vector, f is the body force, ρ isthe density of the fluid, P is the pressure and µ is thedynamic viscosity of the fluid. Secondly, theconcentration distribution can be calculated fromconvection–diffusion equation, Eq. (13), after obtainingthe velocity field:2( ) Dtφ φ φ∂ + ⋅∇ = ∇∂ v (13)where, φ is the concentration of species, D is thediffusion coefficient and ν is the fluid velocity. Then,the concentration can be used for investigating thedegree of fluid mixing by equation:W0W00- dy( x ) 1 - 100%- dyφ φγφ φ∞∞= ×∫∫ (14)where, W is the width of the micro-channel, φ is theconcentration of species across the width of the mixingchannel, φ0 and φ∞ are the sample concentrationprofiles with completely unmixed (0 or 1) andcompletely mixed states (0.5), respectively.The finite element model of the micromixer is setup using COMSOL software. The mesh model isshown in Figure 3. The number of elements and nodesare 3936 and 5095 respectively.Fig. 2  The structure of the micromixerFig. 3  The finite element model of the micromixerWe have solved the coupling problem by usingCOMSOL's sequential coupling. The coupling in thematrix equation is shown as:U1U200 0 0φφ φ φ⎡ ⎤⎡ ⎤ ⎡ ⎤⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤+ + =⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦⎣ ⎦ ⎣ ⎦D 0 FM u u K K uFφ φ φD D K  (15)where, u and ϕ are velocity and concentration degreesof freedom, respectively. The system matrices M, K,Dϕ and Kϕ are assembled following a standard finiteelement methodology. Matrices DUϕ and KUϕ representthe coupled effect attribute. COMSOL has solved thesenonlinear equations using CG solver based on AMG.X. Chen, H. Zeng, H. Wang, J. Li: Nonlinear macromodel based on Krylov subspace for micromixer of the microfluidic chip48 ENGINEERING MODELLING 25 (2012) 1-4, 45-494. MACROMODEL EXTRACTION FORMICROMIXERThe COMSOL software was used to conduct finiteelement simulation for the micromixer and somecharacteristic parameters are set up and presented inTable 1.Table 1. Finite element simulation parameters5. CONCLUSIONSA macromodel is presented and successfullycalculated for the microfluidic chip micromixer. It hasbeen proved that the macromodel based on Kylesubspace is an effective method for the MEMSsimulation. The strongly coupled nonlinear microfluidicproblems can be solved by the use of the macromodelmethod in combination with the finite element analysis.In the future, the macromodel method based on Kylesubspace will have a great potential in the design andsimulation of microfluidic analysis systems.6. ACKNOWLEDGMENTSThe authors are grateful to Liaoning Province KeyLaboratory Funded Projects (No: LS2010080,LKF2011020), and National Science TechnologySupport Plan Projects (No: 2012BAF12B08-5).7. REFERENCES[1] J. Cao, P. Cheng and F.J. Hong, A numerical studyof an electrothermal vortex enhanced micromixer,Microfluidics and Nanofluidics, Vol. 5, No. 1, pp.13-21, 2008.[2] C.A. Cortes-Quiroz, M. Zangeneh and A. Goto,On multi-objective optimization of geometry ofstaggered herringbone micromixer, Microfluidicsand Nanofluidics,  Vol. 7, No. 1, pp. 29-43, 2009.[3] W. Dai and Y. P. Zhao, The nonlinear phenomenaof thin polydimethylsiloxane (PDMS) films inelectrowetting, International Journal ofNonlinear Sciences and Numerical SimulationVol. 8, No. 4, pp. 519-526, 2007.[4] O.-T. Son, J.W. Roh, S.-H. Song, J.-S. Park, W.Lee and H.-I. Jung, Continuous micro-magnetophoretic separation using a dipolemagnetic field, BioChip Journal, Vol. 2, No. 3,pp. 186-191, 2008.[5] Y.-Q. Wan, Q. Guo and N. Pan, Thermo-electro-hydrodynamic model for electrospinning process,International Journal of Nonlinear Sciences andNumerical Simulation, Vol. 5, No. 1, pp. 5-8,2004.[6] Y.C. Zeng, Y. Wu, Z.G. Pei and C.W. Yu,Numerical approach to electrospinning,International Journal of Nonlinear Sciences andNumerical Simulation,  Vol. 7, No. 4, pp. 385-388, 2006.[7] D.H. Shou and J.H. He, Application of parameter-expanding method to strongly nonlinearoscillators, International Journal of NonlinearSciences and Numerical Simulation,  Vol. 8, No.1, pp. 121-124, 2007.Simulation parameter Value Density of fluid Dynamic viscosity o f the fluid Diffusion coefficient Channel width Wall boundary 1 input concentration 2 input concentration 103 (kg/m3) 1.002×10-3 (Pa·s) 10-10 (m2/s) 3×10-4 (µm) No slip 0.01 (mol/m3) 0.1 (mol/m3) The command “assemble” can extract stiffness matrixand other system matrices from COMSOL finite elementmodel analysis outputs, and the command code is:[ , , , ] assemble( fem );spy( )=K L M N K (16)Using the system matrices extracted from the finiteelement model, Lanczos algorithm is executed inMATLAB language for the steady-state analysis of themicromixer. We reduced the original model to 40 and 20order in MATLAB language and we also calculated themacromodel output value. The calculation times ofCOMSOL, 40 order macromodel and 20 ordermacromodel are 2.2 s, 0.027 s and 0.021 s, respectively.The macromodel is two orders of magnitude faster thanthe finite element model. The micromixer outputconcentration results from COMSOL, 40 and 20 ordermacromodel are shown in Figure 4. As it can be seenfrom Figure 4, the calculation error from 40 order modelis less 5% less than the one obtained by COMSOL,therefore, it can be concluded that the calculation errordecreases along with the order number increment.Furthermore, it means that the compromise betweencalculation precision and speed should be taken intoconsideration. 0.0 0.2 0.4 0.6 0.8 1.00.3400.3420.3440.3460.3480.3500.352ConcentrationChannel coordinate Coordinates Coordinates CoordinatesFig. 4  COMSOL model, 40 order macromodel and 20 ordermacromodel output concentration resultsX. Chen, H. Zeng, H. Wang, J. Li: Nonlinear macromodel based on Krylov subspace for micromixer of the microfluidic chipENGINEERING MODELLING 25 (2012) 1-4, 45-49 49[8] Y. Wang, Q. Lin and T. Mukherjee, A model forJoule heating-induced dispersion in microchipelectrophoresis, Lab on a Chip,  Vol. 4, No. 6,pp. 625-631, 2004.[9] D. Vasilyev, M. Rewienski and J. White,Macromodel generation for bioMEMScomponents using a stabilized balanced truncationplus trajectory piecewise-linear approach, IEEETransactions on Computer-Aided Design ofIntegrated Circuits and Systems, Vol. 25, No. 2,pp. 285-293, 2006.[10] H.M. Park and J.Y. Lim, A reduced-order modelof the low-voltage cascade electroosmoticmicropump, Microfluidics and Nanofluidics,Vol. 6, No. 4, pp. 509-520, 2009.[11] R. Qiao and N.R. Aluru, A compact model forelectroosmotic flows in microfluidic devices,Journal of Micromechanics and Micro-engineering,  Vol. 12, No. 5, pp. 625-635, 2002.[12] X. He and S. Hauan, Microfluidic modeling anddesign for continuous flow in electrokineticmixing-reaction channels, Aiche Journal, Vol. 52,No. 11, pp. 3842-3851, 2006.[13] M. Troyer, H. Tsunetsugu and D. Wurtz, Thermo-dynamics and spin gap of the Heisenberg laddercalculated by the look-ahead Lanczos-algorithm,Physical Review B,  Vol. 50, No. 18, pp. 13515-13527, 1994.[14] S. Sundar and B.K. Bhagavan, Comparison ofKrylov subspace methods with preconditioningtechniques for solving boundary value problems,Computers & Mathematics with Applications,Vo. 38, No.11-12, pp. 197-206, 1999.[15] T.J. Frankcombe and S.C. Smith, Fast, scalablemaster equation solution algorithms, IV. Lanczositeration with diffusion approximationpreconditioned iterative inversion, Journal ofChemical Physics, Vol. 119, No. 24, pp. 12741-12748, 2003.PRIKAZ NELINEARNOG MAKROMODELA TEMELJENOG NA KRYLOVLJEVOMPOTPROSTORU PRIKLADNOG ZA SIMULACIJU RADA MIKROMIKSERAMIKROFLUIDNOG PROCESORASA@ETAKU ovome se radu polazi od pretpostavke da konvencionalne metode numeri~ke analize nisu dostatne za simulacijumikro-elektro-mehani~kog-sustava (MEMS) koji sadr`i polje fluida. Me|utim, korištenjem makromodela s primjenomnelinearne analize mogu se simulirati karakteristike mikro-polja toka. U ovome su radu prikazane postavkemakromodela mikromiksera mikrofluidnog procesora pri ~emu je korištena Krylovljeva metoda projiciranjapotprostora. Funkcije sustava dobivene su metodom kona~nih elemenata uz korištenje softvera COMSOL. Zadiskretizirani model mikromiksera uzeta je analiza polja koncentracije. Red funkcije kona~nih elemenata reduciranje pomo}u Krylovljeve metode projekcije potprostora drugog reda bazirane na Lanczosovom algoritmu. Pokazano jeda se rezultati dobiveni korištenjem makromodela podudaraju s rezultatima dobivenim analizom s primjenom kona~nihelemenata, s tim da je vrijeme trajanja simulacije pri korištenju makromodela stotinu puta kra}e od vremena trajanjasimulacije korištenjem metode kona~nih elemenata. U prakti~noj primjeni, ovaj bi makromodel trebao olakšatidizajniranje mikrofluidnih naprava sa sofisticiranom mre`om kanala.Klju~ne rije~i: makromodel, Krylovljev potprostor, mikromikser, mikro-elektro-mehani~ki sustavi (MEMS).

Nonlinear macromodel based on Krylov subspace for micromixer of the microfluidic chip

Abstract

Similar works

Full text

Available Versions

HRČAK - Portal of Croatian Scientific and Professional Journals

Hrčak - Portal of scientific journals of Croatia