COMSOL 6.4 - Micropump Mechanism

Micropump Mechanism¹

Micropumps are key components of microfluidic systems with applications ranging from biological fluid handling to microelectronic cooling. This model simulates the mechanism of a valveless micropump, that is designed to be effective at low Reynolds numbers, overcoming hydrodynamic reversibility. Valveless pumps are often preferred in microfluidic systems because they minimize the risk of clogging and are gentle on biological material. The Fluid-Structure Interaction interface is used to solve for the flow of the fluid and the associated deformation of the structure. In addition, the Global ODEs and DAEs interface is used to demonstrate how to perform a time-resolved integration of the total flow throughout the pumping cycle.

Introduction

Many valveless pump designs are ineffective when the system has a low Reynolds number, and consequently are unsuitable for viscous fluids and applications with small length scales or low flow rates. This is largely because without valves it is difficult to achieve sustained flow in a given direction.

The mechanism simulated in this model overcomes this limitation by converting oscillatory fluid motion, induced by a simple reciprocating pumping mechanism, into a net flow in one direction. It is relatively easy to create an oscillatory pumping mechanism in a microfluidic system, for example, a membrane can be vibrated by a piezo oscillator to periodically vary the volume of a microchamber. In this model, an oscillatory flow is fed into a channel containing bendable microflaps. The deformation of the microflaps in response to the motion of the fluid alters the flow and results in a net flow rate in a consistent direction. The passive nature of the flow regulator allows for directional control of the flow without the use of the complicated synchronized actuation mechanisms that would be required in a valve based system.

In this model the Fluid-Structure Interaction Multiphysics interface is used to specify the input oscillatory flow, along with the mechanical properties of the flaps. The deformation of the flaps, and the flow of the fluid, is calculated as a function of time for two full cycles of the pumping mechanism. This allows the physical mechanism responsible for generating the unidirectional flow to be visualized clearly using an animation. As well as visualizing flow rate and direction as a function of time throughout the pumping cycle, integration coupling components are used in conjunction with the Global ODEs and DAEs interface to calculate the net volume pumped from left-to-right as a function of time. This is an example of how the functionality of one COMSOL interface can be enhanced by using a custom equation specified in a mathematics interface, and demonstrates the ease with which user defined equations can be incorporated into COMSOL models.

Model Definition

The model geometry is shown in Figure 1. It consists of a horizontal channel that is 600 μm in length and 100 μm high. A vertical chamber connects to the channel at the midpoint along its length. Two tilted flaps are attached to the bottom of the channel such that they partially obstruct flow along the channel length. They are spaced to be centered on the midpoint of the channel length, and they are both angled at 45 degrees to the horizontal channel edge. Note that this 2D model represents a cross-section through the midpoint of the channel in the out-of-plane direction. An out-of-plane thickness of 10 μm has been used for the purpose of calculating the volume of fluid pumped as a function of time. However, as no edge effects due to walls that are out-of-plane are included, this is equivalent to modeling a 10 μm deep section of a much thicker channel.

Figure 1: The model geometry consists of a horizontal channel and a vertical chamber. Tilted flaps are positioned within the channel, the response of these flaps to the oscillatory fluid motion induced via the labeled boundary results in a net flow from left-to-right.

The physics required for the model is configured within the Fluid-Structure Interaction multiphysics interface. An Inlet boundary feature is applied to the top of the vertical chamber. This specifies the inflow velocity, via a user input expression, to vary sinusoidally in time with a period of 1 s. An Outlet boundary feature is applied to the left and right boundaries of the channel. Two boundary integration coupling components, named intopL() and intopR() for the left and right outlets respectively, are also applied to these outlet boundaries. These are used to compute the flow rate out of each outlet. This is achieved using some user defined variables in the node within . The flow rate from each outlet is calculated by integrating the depended variable u_fluid, which is the horizontal component of the fluid velocity, and multiplied by the out-of-plane length scale of 10 μm. The net flow rate out of the channel, UoutNet, is then calculated from the difference between the flow from the left and right outlets, such that positive values correspond to a net flow in the left-to-right direction.

A Global ODEs and DAEs interface is added to compute the integrated net flow as a function of time. This is achieved using a Global Equation which integrates UoutNet with respect to time to obtain Vpump. This step is necessary as UoutNet gives the instantaneous net flow rate as a function in time, however a more useful metric for evaluating a pump is the total volume pumped throughout an entire pumping cycle. Note that the timeint() operator can also be used to visualize the time integration of a variable. However, use of timeint() is not as accurate as directly solving for an integrated quantity, as the timeint() operator only uses the time steps which are saved in the solution but solving directly uses every time step taken by the solver.

The mesh is configured to be tightest around the tilted flaps, in order to resolve the stress within the bending flaps. The mesh is shown in Figure 2.

Figure 2: Initial mesh prior to any structural deformation.

A single time dependent study is performed over a duration of 2 seconds, which corresponds to two full oscillations of the inlet velocity. Some minor amendments to the default solver sequence, described in the step-by-step instructions, are required to ensure reliable results. In particular, the Time Stepping sets the maximum step to 0.01 s. This is needed because the feedback between the flap deformation and the fluid flow can lead to abrupt changes in the flow. The default time stepping setting allows the solver to vary the time step it takes so that larger steps can be taken during times when the solution is not rapidly varying. However, because of the potential for sudden changes to the fluid flow in fluid-structure interaction models it can be helpful to impose a maximum time step to ensure that rapid changes are not missed.

It should be noted that the geometry is parameterized so that the channel dimensions and flap angle can be modified by amending the relevant entry in the Global Parameters table. The average flow rate at the inlet and the fluid properties are also parameterized in the same way. In addition, the effective Reynolds number can be easily changed as the viscosity and average inlet velocity are appropriately scaled by a shared coefficient (the parameter coeff) that is computed from the target Reynolds number (the parameter Re). Note that the parameter Re_check is provided to confirm that the Reynolds number does indeed take the specified value. The model is set up with a Reynolds number of 16, but it is straightforward to verify that the pumping action occurs even for Reynolds numbers significantly less than one.

Results and Discussion

The mechanism by which the flow direction is regulated can be observed in a combined Flow and Stress plot. During the downstroke, when fluid is pushed from the vertical chamber into the channel, the right-hand flap is bent down toward the bottom of the channel whilst the left-hand flap is bend away from the channel bottom. This configuration is shown in Figure 3, which shows the solution at a time of 0.26 s which corresponds to when the velocity flowing into the vertical chamber via the inlet is at its maximum. Due to the asymmetric bending of the flaps, fluid can more easily flow out of the right-hand outlet. During the upstroke, when fluid is drawn from the channel into the vertical chamber, the flaps are bent in the opposite directions. This configuration is shown in Figure 4, which shows the solution at a time of 0.74 s. Now the right hand flap restricts the flow more than the left-hand flap, and the majority of the fluid that is drawn into the vertical chamber is from the left-hand outlet.

The result of this behavior is that there is a net flow rate from left-to-right inside the channel. This has many possible applications in microfluidic systems. For example, this device could be used to deliver fluid from a droplet reservoir connected to the left-hand outlet into a microfluidic pathway connected to the right-hand outlet. Alternately, this device could be used to create a circulating system where a fluid is pumped around a continuous loop to cool a microelectronic system.

Figure 3: Velocity magnitude and velocity field, along with the von Mises stress within the flaps during the pumping down-stroke.

Figure 4: Velocity magnitude and velocity field, along with von Mises stress within the flaps, during the pumping upstroke.

The net volume of fluid that is pumped from left-to-right is shown in Figure 5. As expected, the gradient of the curve, which is the net flow rate, varies sinusoidally with a period equal to the inlet velocity. The maximum gradients occur at intervals of odd multiples of 0.5 s, which correspond with the peaks in the magnitude of the inlet velocity.

Figure 5: Net volume pumped from left-to-right as a function of time.

MEMS_Module/Fluid-Structure_Interaction/micropump_mechanism

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Geometry 1

Add some parameters which will be used to control the fluid properties and device geometry.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
Re	16	16	Reynolds number
coeff	4/sqrt(Re)	1	Coefficient to change Reynolds number
dens	1000[kg/m^3]	1000 kg/m³	Fluid density
visc	0.001[Pas]coeff	0.001 Pa·s	Fluid dynamic viscosity
U	16[cm/s]/coeff	0.16 m/s	Average inlet flow speed
H	100[um]	1E-4 m	Channel height
W	10[um]	1E-5 m	Domain width
rp	2[um]	2E-6 m	Pillar radius
hp	70[um]	7E-5 m	Pillar height
L	600[um]	6E-4 m	Length of channel
beta	45[deg]	0.7854 rad	Flap tilt angle
x0	150[um]	1.5E-4 m	Flap center location
Re_check	densUH/(visc)	16	Reynolds number

Next create the geometry for the device.

Geometry 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type L.

In the text field, type H.

Locate the section. In the text field, type -L/2.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2*rp.

In the text field, type 2*hp.

Locate the section. From the list, choose .

In the text field, type x0-hp*sin(beta)/2.

Locate the section. In the text field, type -beta.

Copy 1 (copy1)

In the toolbar, click

and choose .

Select the object only.

Intersection 1 (int1)

In the toolbar, click

and choose .

Select the objects and only.

In the window for , locate the section.

Clear the checkbox.

Fillet 1 (fil1)

In the toolbar, click

On the object , select Points 3 and 4 only.

It might be easier to select the correct points by using the window. To open this window, in the toolbar click and choose . (If you are running the cross-platform desktop, you find in the main menu.)

In the window for , locate the section.

In the text field, type rp.

Copy 2 (copy2)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

In the text field, type -2*x0.

Rectangle 3 (r3)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.8*H.

In the text field, type H.

Locate the section. In the text field, type -0.4*H.

In the text field, type H.

Click

Add two nonlocal integration couplings, one to the boundary which forms the left outlet and one to the boundary which forms the right outlet. These couplings will be used to calculate the flow rate out of each outlet.

Definitions

Integration 1 (intop1)

In the toolbar, click

and choose .

In the window for , type intopL in the text field.

Locate the section. From the list, choose .

Select Boundary 1 only.

Integration 2 (intop2)

In the toolbar, click

and choose .

In the window for , type intopR in the text field.

Locate the section. From the list, choose .

Select Boundary 17 only.

Create some variables to calculate the net flow rate in the channel. The first two variables calculate the average flow rate out of each outlet by performing an average over u_fluid , which is the horizontal component of the fluid flow. Note the sign convention, that assigns positive values to left-right flow and negative values to right-left flow. The third variable calculates the difference between the average flow rate out of each outlet, positive values correspond to a net flow from left to right.

Variables 1

In the window, right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
UoutL	-intopL(u_fluid)*W	m³/s	Flow rate from left outlet
UoutR	intopR(u_fluid)*W	m³/s	Flow rate from right outlet
UoutNet	UoutR-UoutL	m³/s	Net flow rate

In order to integrate the average net flow rate over time to calculate the volume of fluid that is pumped, add a Global ODES and DAEs interface to the model.

Add Physics

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Global ODEs and DAEs (ge)

Global Equations 1 (ODE1)

In the window for , locate the section.

In the table, enter the following settings:


Name	f(u,ut,utt,t) (1)	Initial value (u_0) (1)	Initial value (ut_0) (1/s)	Description
Vpump	Vpumpt-UoutNet	0	0	Net volume pumped left-to-right

Locate the section. Click

In the table, enter the following settings:


Dependent variable quantity	Unit
Custom unit	m^3

Click

In the table, enter the following settings:


Source term quantity	Unit
Custom unit	m^3/s

Configure the interface.

Moving Mesh

Deforming Domain 1

In the window, under > click .

In the window for , locate the section.

In the C2 text field, type 100.

Select Domain 1 only.

Fixed Boundary 1

In the toolbar, click

Select Boundaries 1, 2, 7, 9, 16, and 17 only.

Prescribed Normal Mesh Displacement 1

In the toolbar, click

Select Boundaries 3 and 12 only.

Laminar Flow (spf)

In the window, under click .

Select Domains 1 and 3 only.

Inlet 1

In the toolbar, click

and choose .

Select Boundary 10 only.

In the window for , locate the section.

In the U0 text field, type U*6*edgparnal*(1-edgparnal)*sin(2*pi*t/(1[s])). Here edgparnal is the normalized arc length parameter, a predefined variable in COMSOL Multiphysics.

Outlet 1

In the toolbar, click

and choose .

Select Boundaries 1 and 17 only.

In the window for , locate the section.

Clear the checkbox.

Solid Mechanics (solid)

In the window, under click .

Select Domains 1, 2, and 4 only.

Linear Elastic Material 1

In the window, under > click .

In the window for , locate the section.

From the list, choose .

Fixed Constraint 1

In the toolbar, click

and choose .

Select Boundaries 4 and 13 only.

With the interfaces configured, add materials to the geometry domains. In this case, it was beneficial to wait until all the physics settings were configured before adding materials, as COMSOL automatically adjusts which materials properties are needed depending on the physics assigned to each domain.

Materials

Fluid

In the window, under right-click and choose .

In the window for , type Fluid in the text field.

Solid

Right-click and choose .

In the window for , type Solid in the text field.

Select Domains 2 and 4 only.

Fluid (mat1)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	dens	kg/m³	Basic
Dynamic viscosity	mu	visc	Pa·s	Basic

Solid (mat2)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Young’s modulus	E	3.6e5	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.499	1	Young’s modulus and Poisson’s ratio
Density	rho	970	kg/m³	Basic

Configure the mesh.

Mesh 1

Free Triangular 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Size

In the window, click .

In the window for , locate the section.

Click the button.

Locate the section. In the text field, type 10.

In the text field, type 1.4.

In the text field, type 0.6.

In the text field, type 0.7.

Size 1

In the window, right-click and choose .

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Domains 2 and 4 only.

Locate the section. Click the button.

Locate the section.

Select the checkbox. In the associated text field, type 2.

Select the checkbox. In the associated text field, type 1.5.

In the window, right-click and choose .

Configure the study.

Set the time steps and range for the study.

Study 1

Step 1: Time Dependent

In the window, under click .

In the window for , locate the section.

In the text field, type range(0,0.02,2).

In the toolbar, click

Results

Velocity (spf)

Modify the default plot group to show the von Mises stress within the flexible rods, as well as the velocity magnitude of the fluid within the channel. Add an plot of the fluid velocity. An animation of the data series allows the action of the pump to be visualized.

Study 1/Solution 1 (2) (sol1)

In the window, expand the > node.