A Piezoelectric Micropump

This tutorial model is kindly provided by Riccardo Vietri, James Ransley, and Andrew Spann at Veryst Engineering, LLC. Piezoelectric micropumps are frequently used in medical applications because of their ability to precisely control the metering of very small volumes of fluids or gases. This example shows how to simulate a simple, non-resonant micropump, suitable for low flow rate applications. The model demonstrates the combination of piezoelectric materials with a fluid-structure interaction, and also illustrates the use of a simple velocity-dependent formula to account for the presence of valves on the inlet and outlet boundaries.

Model Definition

The model geometry consists of an annular piezoelectric actuator on top of the fluid domain which is connected to a flexible membrane (highlighted in Figure 1). Due to the symmetry of the physics, only half of the geometry needs to be included. The fluid flows out when a voltage is applied, and the actuator expands. The circumference of the membrane is fixed as the actuator moves, and it imposes a force on the fluid beneath it, pulling in fluid from the left and pushing fluid out of the channel on the right.

Figure 1: The geometry of the piezoelectric micropump. The membrane is highlighted.

In a real micropump, a stacked piezoelectric actuator would be employed with many separate piezoelectric layers and electrical contacts. For the purposes of this model, neglect the stiffness of the thin metal layers and model the actuator as a monolithic block of piezoelectric. Consequently, apply a high potential difference of 1500 V across the entire piezoelectric. This corresponds to an electric field of 0.2 V/µm. The voltage required in a real device would depend on the thickness of each layer in the stacked actuator.

The inlet and outlet work by employing check valves to ensure a one-way flow. In the model the valve is represented by a simple boundary condition based on K-factor piping losses, where the losses are high when flowing against the valve and low when flowing in the direction of operation. The back pressure resulting from the valve is represented by the following equation:

(1)

where uav is the average velocity of the fluid normal to the boundary, ρ is the fluid density and A is a dimensionless constant that changes magnitude depending on the sign of uav. The back pressure is applied as a normal stress at the end of a short length of pipe - which ensures that the fluid flow in the domain is realistic, despite this approximate boundary condition. This boundary condition can be employed to represent a simple fluid valve or diode. The constants used for the outlet boundary are reversed with respect to those used in the inlet, representing a different orientation of a similar valve. This encourages flow through the pump in the desired direction. To represent a low resistance valve (such as a simple flapper valve), set A to 5,000 for closed condition and to 0.1 for the open condition. Careful tuning of these values and potentially refinement of this model would be desirable in a real application.

To save computation time for the purposes of this example, a relatively coarse mesh is used, as shown in the figure below.

Figure 2: The mesh used in the model. Note the fluid boundary layers next to the highlighted fluid walls.

Results and Discussion

Figure 3 shows the inlet and outlet flow rates and confirms the conservation of fluid volume within the device. The drive voltage is ramped up during the first 3/4 of the actuation period. Afterwards, a consistent time periodic flow is quickly established. The difference in inlet and outlet flow matches the volume of the fluid displaced by the membrane due to the piezoelectric stroke, confirming the volume conservation.

Figure 3: Flow rates and volume conservation.

The time integrated flows into the inlet and out of the outlet are calculated using a global ODE and are shown in Figure 4. The net flow through the pump in 0.05 s is 1.8 µL, corresponding to 36 µL/s or 2.16 ml/min. This is typical for such a nonresonant design.

Figure 4: Net fluid flow through the inlet and the outlet over time, computed using a global ordinary differential equation.

Figure 5 shows the flow through the inlet (left) and outlet (right) when fluid is being drawn into the chamber, whilst Figure 6 shows the case when fluid is being forced out of the chamber.

Figure 5: Velocity (base 10 logarithm in mm/s) when fluid is being drawn through the inlet.

Figure 6: Velocity (base 10 logarithm in mm/s) when fluid is being pushed through the outlet.

Figure 7 and Figure 8 show the fluid streamlines in the corresponding cases.

Figure 7: Flow streamlines when fluid is being drawn through the inlet. The color shows the base 10 logarithm of the velocity in mm/s.

Figure 8: Flow streamlines when fluid is being pushed through the outlet. The color shows the base 10 logarithm of the velocity in mm/s.

MEMS_Module/Fluid-Structure_Interaction/piezoelectric_micropump

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In addition to piezoelectricity and laminar flow, a global equation will be used to keep track of the total flow through the inlet and the outlet.

In the tree, select .

Click .

In the tree, select .

Click .

In the tree, select .

Click .

Click

In the tree, select .

Click

Geometry 1

The Model Wizard starts the COMSOL Desktop at the node. Take the opportunity to set the length unit to mm for convenience.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Enter some global parameters. The expression containing the parameter r_inlet appears red initially and turns to black once r_inlet is defined.

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
w_block	30 [mm]	0.03 m	Width of base
depth_block	30 [mm]	0.03 m	Depth of base
h_block	5 [mm]	0.005 m	Base thickness
h_exit	5*r_inlet	0.005 m	Height of Inlet/Outlet
h_memb	1 [mm]	0.001 m	Height of the Membrane
ID	8[mm]	0.008 m	Disc actuator inner diameter
OD	15[mm]	0.015 m	Disc actuator outer diameter
r_inlet	1 [mm]	0.001 m	Fluid Inlet Radius
r_memb	12 [mm]	0.012 m	Radius of the Membrane
r_outlet	1 [mm]	0.001 m	Fluid Outlet Radius
t0	0.1[mm]	1E-4 m	Piezoelectric layer thickness
n	75	75	Number of layers in actuator
E0	0.2[V/um]	2E5 V/m	Electric field strength
V0	E0t0n	1500 V	Applied voltage
frequency	60[Hz]	60 Hz	Frequency of piston actuation
high_stress	5e3	5000	Boundary Stress (High)
low_stress	1e-1	0.1	Boundary Stress (Low)

Build the geometry. We will model only half of the device using the symmetry of the problem. Use with to allow different meshes in the solid and the fluid domains, while still creating a matching pair of surfaces for the solid-fluid interface.

The fluid chamber is a simple block.

Geometry 1

Block 1 - Fluid Chamber

In the toolbar, click

In the window for , type Block 1 - Fluid Chamber in the text field.

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

In the text field, type depth_block.

In the text field, type h_block-h_memb.

Locate the section. From the list, choose .

In the text field, type -h_memb/2.

The piezoelectric stack is a hollow cylinder, which can be constructed by taking the difference between two cylinders.

Cylinder 1 - Piezo OD

In the toolbar, click

In the window for , type Cylinder 1 - Piezo OD in the text field.

Locate the section. In the text field, type OD/2.

In the text field, type t0*n.

Locate the section. In the text field, type h_block/2.

Cylinder 2 - Piezo ID

Right-click and choose .

In the window for , type Cylinder 2 - Piezo ID in the text field.

Locate the section. In the text field, type ID/2.

Difference 1 - Piezo

In the toolbar, click

and choose .

In the window for , type Difference 1 - Piezo in the text field.

Select the object only.

Locate the section. Find the subsection. Select the

toggle button.

Select the object only.

The actuated membrane is a simple cylinder, so are the inlet and outlet pipes.

Cylinder 3 - Membrane

In the toolbar, click

In the window for , type Cylinder 3 - Membrane in the text field.

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

In the text field, type h_memb.

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

Cylinder 4 - Inlet

In the toolbar, click

In the window for , type Cylinder 4 - Inlet in the text field.

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

In the text field, type h_exit.

Locate the section. In the text field, type -(OD+ID)/4.

In the text field, type -h_block/2-h_exit.

Cylinder 5 - Outlet

Right-click and choose .

In the window for , type Cylinder 5 - Outlet in the text field.

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

Locate the section. In the text field, type (OD+ID)/4.

Use a workplane to divide the geometry to two symmetric parts and remove one of them.

Work Plane 1 - Symmetry Plane

In the toolbar, click

In the window for , type Work Plane 1 - Symmetry Plane in the text field.

Locate the section. From the list, choose .

Partition Objects 1 (par1)

In the toolbar, click

and choose .

Select the object only.

Click in the window and then press Ctrl+A to select all objects.

In the window for , locate the section.

From the list, choose .

Delete Entities 1 (del1)

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

On the object , select Domain 1 only.

Unite the solid and fluid parts of the geometry separately with a selection defined for each, and then finalize the geometry using .

Union 1 - Solid

In the toolbar, click

and choose .

In the window for , type Union 1 - Solid in the text field.

Select the objects and only.

Locate the section. Select the check box.

Union 2 - Fluid

In the toolbar, click

and choose .

In the window for , type Union 2 - Fluid in the text field.

Select the objects , , and only.

Locate the section. Select the check box.

Form Union (fin)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Select the check box.

In the toolbar, click

Create selections for the subsequent settings. First define the piezoelectric stack domain using a box selection.

Definitions

Box 1 - Piezo

In the toolbar, click

In the window for , type Box 1 - Piezo in the text field.

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

Locate the section. From the list, choose .

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

In the text field, type h_block/2-h_memb/2.

The membrane domain is given by the difference between the solid domain and the piezoelectric stack.

Difference 1 - Membrane

In the toolbar, click

In the window for , type Difference 1 - Membrane in the text field.

Locate the section. Under , click

In the dialog box, select in the list.

Click .

In the window for , locate the section.

Under , click

In the dialog box, select in the list.

Click .

Define the inlet and outlet boundaries using box selections.

Box 2 - Inlet

In the toolbar, click

In the window for , type Box 2 - Inlet in the text field.

Locate the section. From the list, choose .

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

In the text field, type -h_block/2-h_exit/2.

Locate the section. From the list, choose .

Use wireframe rendering to see the inlet boundary selection better.

Click the

button in the toolbar.

Box 3 - Outlet

Right-click and choose .

In the window for , type Box 3 - Outlet in the text field.

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

In the text field, type Inf.

Use a box selection to define the symmetry plane.

Box 4 - Symmetry Plane

In the toolbar, click

In the window for , type Box 4 - Symmetry Plane in the text field.

Locate the section. From the list, choose .

Locate the section. In the text field, type min(min(r_inlet,r_outlet),ID)/2.

Locate the section. From the list, choose .

Create a selection for the fluid walls by removing the symmetry plane and the in/outlets from the boundaries adjacent to the fluid. This selection can be used to create boundary layer mesh for the flow problem.

Adjacent 1 - All Fluid Boundaries

In the toolbar, click

In the window for , type Adjacent 1 - All Fluid Boundaries in the text field.

Locate the section. Under , click

In the dialog box, select in the list.

Click .

Difference 2 - Fluid Walls

In the toolbar, click

In the window for , type Difference 2 - Fluid Walls in the text field.

Locate the section. From the list, choose .

Locate the section. Under , click

In the dialog box, select in the list.

Click .

In the window for , locate the section.

Under , click

In the dialog box, select in the list.

Click .

In the window for , locate the section.

Under , click

In the dialog box, select in the list.

Click .

In the window for , locate the section.

Under , click

In the dialog box, select in the list.

Click .

Create a selection for the membrane boundary to be used to keep track of the volume change of the fluid chamber.

Box 5

In the toolbar, click

In the window for , locate the section.

From the list, choose .

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

In the text field, type eps.

In the text field, type -eps.

In the text field, type eps.

In the text field, type h_block/2-h_memb-eps.

Intersection 1 - Fluid Membrane

In the toolbar, click

In the window for , type Intersection 1 - Fluid Membrane in the text field.

Locate the section. From the list, choose .

Locate the section. Under , click

In the dialog box, in the list, choose and .

Click .

Finally select the top surface of the piezoelectric stack and the circumference of the membrane for the fixed constraint in the solid problem.

Explicit 1 - Fixed Boundaries

In the toolbar, click

In the window for , type Explicit 1 - Fixed Boundaries in the text field.

Locate the section. From the list, choose .

Select Boundaries 18, 25, and 30 only.

Add material data to the corresponding selections.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Click in the window toolbar.

In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Materials

Water, liquid (mat1)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Lead Zirconate Titanate (PZT-5H) (mat2)

In the window, click .

In the window for , locate the section.

From the list, choose .

Membrane

In the window, right-click and choose .

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

Locate the section. From the list, choose .

At this point the set of properties shown for the user-defined membrane material is not as intended, because the physics selections haven’t been configured yet.

Set up the physics selections relevant for the material data.

Solid Mechanics (solid)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Piezoelectric Material 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Electrostatics (es)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Laminar Flow (spf)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Now we are ready to enter the material properties for the membrane.

Materials

Membrane (mat3)

In the window, under 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	200[MPa]	Pa	Basic
Poisson’s ratio	nu	0.45	1	Basic
Density	rho	2320[kg/m^3]	kg/m³	Basic

Continue setting up the physics. Add fixed constraint and symmetry boundary condition to solid mechanics, using the boundary selections defined earlier.

Solid Mechanics (solid)

In the window, under click .

Fixed Constraint 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Symmetry 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Add electrical terminals to the top and bottom surfaces of the piezoelectric stack, using a terminal boundary condition on the top and a ground boundary condition at the bottom.

Electrostatics (es)

In the window, under click .

Terminal 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Ground 1

In the toolbar, click

and choose .

Select Boundary 24 only.

The micropump is driven by a sinusoidal voltage at the terminal. Define a ramp function to smoothly ramp up the sinusoidal drive. Then enter the expression for the terminal voltage.

Definitions

Ramp 1 (rm1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type 0.1.

In the text field, type 1.2.

Select the check box.

Click to expand the section. Select the check box.

Select the check box.

Electrostatics (es)

Terminal 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

In the V0 text field, type V0*sin(t*frequency*2*pi)*rm1(t*frequency*4/3).

Use the selections defined earlier to set up the symmetry, inlet, and outlet boundary conditions for the fluid flow.

Laminar Flow (spf)

In the window, under click .

Inlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Outlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Symmetry 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Define average operators to compute the average fluid velocity at the in/outlet and then use it in the analytic formula for the pressure to account for the effect of the valve at the in/outlet.

Definitions

Average 1 - Inlet

In the toolbar, click

and choose .

In the window for , type Average 1 - Inlet in the text field.

In the text field, type av_in.

Locate the section. From the list, choose .

From the list, choose .

Average 2 - Outlet

Right-click and choose .

In the window for , type Average 2 - Outlet in the text field.

In the text field, type av_out.

Locate the section. From the list, choose .

Laminar Flow (spf)

Inlet 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the p0 text field, type if(av_in(w2)>0,-low_stress, high_stress)*(av_in(w2)^2)*av_in(spf.rho).

From the list, choose . Specify the du vector as


0	x
0	y
1	z

Outlet 1

In the window, click .

In the window for , locate the section.

In the p0 text field, type if(av_out(w2)<0,low_stress, -high_stress)*(av_out(w2)^2)*av_out(spf.rho).

Define integration operators to compute the flow rate at the in/outlet and then use it in the global equation to keep track of the accumulated flow volume at the in/outlet.

Definitions

Integration 1 - Inlet

In the toolbar, click

and choose .

In the window for , type Integration 1 - Inlet in the text field.

In the text field, type int_in.

Locate the section. From the list, choose .

From the list, choose .

Integration 2 - Outlet

Right-click and choose .

In the window for , type Integration 2 - Outlet in the text field.

In the text field, type int_out.

Locate the section. From the list, choose .

Global ODEs and DAEs (ge)

Global Equations 1 - Accumulated Flow Volume

In the window, under click .

In the window for , type Global Equations 1 - Accumulated Flow Volume in the text field.

Locate the section. In the table, enter the following settings:


Name	f(u,ut,utt,t) (1)	Initial value (u_0) (1)	Initial value (u_t0) (1/s)	Description
Q_in	Q_int-int_in(w2)	0	0	Accumulated in flow
Q_out	Q_outt-int_out(-w2)	0	0	Accumulated out flow