Design and optimization of piezoelectric actuators for microflap-embedded micropumps

In this paper, we propose a new method to use cost-effective multi-sheet off-the-shelf piezoelectric material (e.g. PZT) as an actuator for micropumps. Instead of one customized single PZT sheet that is typically expensive, multiple commercially available PZT sheets are utilized to decrease the cost of fabrication. For this purpose, we have derived analytic equations for expressing the natural frequency and mode shape of the actuator. The FEM simulations are utilized to verify the analytic equations. Thanks to their high accuracy, we can utilize the derived analytic equations as fitness functions of genetic algorithm (GA) for the optimization purpose of PZT physical aspects. Our experimental measurement results show that the GA is capable of optimizing multiple physical parameters of the piezoelectric actuator. Moreover, one-way compliant microflaps are presented for the first time to act as one-way valves for a PZT micropump aided by our proposed multi-sheet PZT actuator. The flow rate of this configuration is compared with a single-sheet PZT actuator in order to demonstrate the effect of the optimized PZT actuators in the practical applications of micropumps.


Introduction
The concept of microfluidics was introduced over 50 years ago by IBM. 1 Later Gravesen et al. 2 published a review paper about micropumping technologies and different actuating principles. Afterward many new fabrication technologies and in turn new micropump techniques have been developed, such as the positive displacement micropump by Cunneen et al., 3 the electro-osmotic micropump by Chen et al., 4 the pneumatic polydimethylsiloxane (PDMS) micropump by Jeong and Konishi, 5 the piezoelectric micropump by Kaviani et al., 6 and the high-pressure peristaltic micropump by Loth and Fo¨rster. 7 These are just some examples of the new fabrication methods in the microfluidics area, where the micropump design is highly dependent on specific applications. supply for better control. Their device was able to reach a flow rate of 4.5 microL/min. A couple of years later, the research group under Alvarez-Bran˜a et al. 9 worked on micropumps through 3D printing technology. They found that the material and design used to print were related to the power output of the pump. The designs they used had a maximum flow rate of 4.1 microL/min. Within the same year of 2021, Qi and Shinshi 10 focused on a bidirectional micropump, which was driven by cylindrical magnets. These magnets exerted a force on three diaphragms, which could drive the fluid at a flow rate of 8.5 microL/min forward and a flow rate of 9 microL/min backwards.
In 2022, at least two more research teams have contributed to this area. Under the leadership of Sravani et al., 11 a design was developed focusing on the use of a stacked style of single layer piezoelectric actuators within the micropump in order to maximize flow rate while minimizing power input. Their device worked at a driving frequency of 100 Hz to achieve an 800 microL/ min flow rate. But no piezoelectric actuator modeling or its performance analysis was reported from this work. Nishikata et al. 12 and his group went to another direction, looking at a new waveform-driven braille actuator for smoother flow. While achieving a 7.4 microL/min flow rate, they also reduced pulsating flow by 79% and backflow in the pump by 63%.
In general, being applied a voltage, a piezoelectric (e.g. PZT) single layer bends inwards as a micropump. 13 This action can push fluid out of the chamber through the outlet valve. In the suction mode, when the voltage is removed, the PZT layer would back up to allow the fluid to enter the chamber. This reciprocating process would eventually cause the pumping action. A recent attempt of improving PZT micropump was done by Hu et al., 14 who added a new passive layer to the previous single layer design to decrease the resonant frequency as an advantage. The authors verified the finite element method (FEM) model with their experimental setup.
Moreover, an analytic modeling for a single layer circular piezoelectric actuator was developed by Mo et al. 15 They investigated the effect of the thickness ratio and the radius of the piezoelectric layer, which is bonded to a metallic layer with the maximum deflection by using their proposed analytic model. In addition, an analytic modeling method based on the thin plate theory and Kelvin-Voigt laws for a single-layer piezoelectric actuator was proposed by Monemian Esfahani and Bahrami. 16 The authors also studied the vibration of an edge clamp for a rectangular PZT actuator in a fluidic environment. Although effective, the derived analytic modeling in those two works above was not used in any follow-up optimization flow and there was no report from them on the performance comparison between single layer and multilayer piezoelectric actuators.
In the past few years, there have been many different trials involving microvalves within various projects. In 2020, Gunda et al. 17 developed a low power, proportionally controlled piezoelectric microvalve. They were able to achieve a device that could work with as high as 1 bar driving pressure with very little leakage (0.8% of the open flow). Around the same time, Bamido et al. 18 were working on thermally actuated microvalves which were meant to aid in optimizing water usage in agriculture. With their apparatus formed through softlithography techniques, they achieved 60% less flow rate. Near the end of 2020, another team led by Amnache et al. 19 worked on a novel idea of selfadaptive microvalves that were sensitive to environmental conditions to help minimize pumping power. While the process was proven effective, the team admitted integration of such a device onto microchips is still a challenge to overcome. Furthermore, Gong et al.'s 20 research group worked on a new on-chip liquid-metalenabled mechanism. Through their research, they achieved a highly flexible, low leak rate device, which can be controlled to form a deformable valve boss to block the flow path.
Above we introduce the recent advancement mainly in the single-sheet PZT actuators (unimorph actuator). Harris et al. 21 presented a multilayer PZT actuator and introduced a technique for replacing the conventional bonded bulk PZT transducer. Following them, Haldkar et al.'s 22 team conducted their work on multilayer strip piezoelectric bimorphs joined with silicone diaphragm, which improved performance compared to the circular PZT actuator. They claimed the new actuator could increase the flow rate by 12 times, with a lower applied voltage. Another work was done by Liu et al. 23 to investigate a bimorph PZT actuator. They conducted both experiment and FEM simulations and observed an enhancement in the deformation. However, they did not offer a systematic design and optimization methodology, which would allow rapid estimation of critical design parameters (e.g. layer thickness) in order to enable faster development of such systems for higher performance without a need of invoking computationally expensive FEM simulations.
Despite the recent advances (e.g. the strip piezo bimorph disk technique), there is still a need for a reliable low-cost micropump. On the top of the various modeling methods achieved in the existing works above, a dedicated optimization process can even boost the performance of actuators and in turn micropumps. Toward this objective, in this paper we are motivated to analytically model and optimize a multilayer PZT actuator with genetic algorithm (GA) for the fluid pumping purpose. The optimized multilayer PZT actuator can be used in micropumps to control the flow rate. Moreover, the micropump valves can be better designed by using the proposed pseudo rigid-body model technique.
This paper is organized as follows. In Section II, the governing analytic equations of the optimized PZT layer with a bonding layer will be presented with their natural frequency and mode shape derived. In Section III our GA-based optimization flow will be discussed, while the actuator design along with FEM simulation and measurement setup will be presented in Section IV. In Section V, we discuss the design of a typical micropump with either a single-PZT-sheet actuator or a double-PZT-sheet actuator. Then the experimental results will be presented in Section VI. Finally, a conclusion is drawn in Section VII.

Analytic modeling
The primary purpose of an optimized piezoelectric actuator is to provide sufficient bending displacement in the transverse direction. As illustrated in Figure 1, a piezoelectric (e.g. PZT) layer, which might be physically implemented by multiple PZT sheets, is glued with the diaphragm layer. These two layers are assumed to connect each other by a massless and linear epoxy bonding layer. The radius of the PZT layer and bonding layer is r O , while their thicknesses are t p and t b , respectively. In this model, since the thickness of the electrodes is less than 0.5 mm, their effect on the deflection of the PZT layer will be ignored. 6 To derive the analytic expression of the natural frequency and mode shape, the dynamic behavior of the actuator has to be studied. By utilizing the Lagrangian method, we can derive the axisymmetric governing equations as follows 6 : where T 1 and T 2 are the applied tension on the piezoelectric layer and diaphragm layer respectively, w 1 and w 2 are the displacement of the piezoelectric layer and diaphragm layer respectively, k is the stiffness coefficient of the epoxy layer, m 1 and m 2 are the mass per unit area of the piezoelectric layer and diaphragm layer, respectively. The piezoelectric layer, which may physically include multiple PZT sheets, is clamped along the edge such that there is no displacement. As a result, the boundary condition equations can be written as follows: Moreover, the PZT layer has no initial velocity: By applying the separation variation method, w 1 and w 2 can be rewritten as follows: Let's assume Z t ð Þ = Z 0 Á e Àivt , then: These equations can be rewritten as: ð7 À 1Þ ð7 À 2Þ With simplification: . By putting rR i into the equations, the following equations will be obtained: By eliminating R 2 from these equations, the following equation can be obtained: where C 1 and C 2 are equal to: The solution of this equation can be represented as the sum of the roots from the Bessel function as follows: where d is the root of the Bessel function. Now by solving this simple algebraic equation, natural frequency v of the system can be derived. We can see that this analytic model is very valuable, since the designers can use it to estimate the natural frequency of a new PZT layer quickly for different dimensions without any computationally expensive FEM simulations. To calculate the mode shape function, (12) has to be solved. The general solution can be achieved as follows: where d i is: By applying the first boundary condition, because of R 1i (0) 6 ¼ ', b i should be equal to zero. In addition, by applying the second boundary condition R 1i r o , t ð Þ= 0, we can get the following equations: or in the matrix form: In order to have a non-zero solution in (18), the following condition must be satisfied: In the end, the characteristic equation will be obtained as: where d n is the root of the Bessel function, and n denotes the n-the root of the Bessel function. Now by solving this simple algebraic equation, natural frequency v of the system can be calculated. In addition, the mode shapes of the piezoelectric layer and diaphragm layer can be expressed respectively as follows: where a 1n is equal to:

Optimization algorithm
The genetic algorithm (GA) is an evolutionary computation algorithm for optimizing sophisticated problems by mimicking biological evolution. 24 As depicted in Figure 2, the GA optimization starts with the initialization block, where the variables are coded in the form of fixed length binary strings. Each variable can be randomly selected with equal probability. Usually, the first operator that performs on a population is the reproduction, which strives to find appropriate individuals in a population and interpolate them into a mating pool. A number of methods for individual selection have been proposed in the literature, although the main idea of this operation is to choose, duplicate, and insert certain preferable individuals from the current population into the pool. The next operator within GA is the crossover, where typically two individuals are selected from the mating pool and a certain quota of these individuals are exchanged in between. In other words, the recombination between these pairs produces new individuals, called offspring. Finally, the mutation operator is performed to change 1 b from 1 to 0 or vice versa. This process is also random with a very low probability (called mutation rate) on the entire population. All the three genetic operators above are performed on the entire population in one GA generation. Thus, the search and optimization aspect of GA is mainly provided by the crossover and mutation operators. The multi-dimensional search capability offered by GA can effectively prevent it from being entrapped by local optima. Therefore, a significant feature of GA in comparison with the conventional optimization approaches is its advantageous access to the global optimum. In this study, our proposed GA-based optimization method is performed to identify the best physical aspects of the piezoelectric layer to increase the displacement of the diaphragm. The coverage of the electrodes is defined to be identical to the size of the piezoelectric film. In the following section, the capability of the GA-based optimization methodology will be further studied.

Actuator design
Since the analytic equations, which show the relationship among the thickness of the piezoelectric layer, the thickness of the diaphragm layer, the radius of the layers, and the maximum deflection of the layers as discussed in Section II, are able to accurately estimate the thicknesses of the piezoelectric and diaphragm layers in addition to the radius, we have deployed them (mainly (21-2)) as a fitness function of our GA optimization algorithm. We may include multiple physical variables as optimizable parameters such as the piezoelectric thickness (h p ) and diaphragm layer thickness (h d ).
We have implemented the GA-based optimization method in MATLAB by using its genetic algorithm tool box (Version 2017) in order to enhance the vertical displacement amplitude of the layers. The applied fitness function and constraints of the optimization are defined by (22): Maximize : fR 2 g, Subjectto : design rules of the optimizable parameters: For the optimizable physical parameters (e.g. h p or h d ), the upper and lower bounds can be defined as per their design rule constraints (e.g. [1m, 1000 m]). According to the listed resultant data in Table 1, the optimized thickness for the piezoelectric layer is 92.95 mm. In other word, based on the other factors and parameters, such a thickness of the piezoelectric layer can offer the maximum displacement performance. Instead of using any expensive customized one-layer piezoelectric material, this required thickness can be achieved by bonding multiple cost-effective commercially available off-the-shelf piezoelectric sheets. For our case study, we have chosen to use standard commercial piezoelectric sheets (i.e. PZT7BB) with constant radius (i.e. 7.5 mm). By using the closest integer with reference to the optimized piezoelectric layer thickness of 92.95 mm (unimorph actuator)   . 42 mm). The material properties within the system are shown in Table 2.
In order to calculate the natural frequency by using the analytic method, equation (20) must be solved. By solving this simple algebraic Bessel function, natural frequency v of the system can be also derived. To validate our methodology, we have run two groups of FEM simulations. In Group-1, one piezoelectric layer (unimorph actuator) with the optimized thickness was simulated in COMSOL Multiphysics (Version 5.3). Figure 3 shows the natural frequencies and mode shapes of these simulations. In the simulations for Group-2, the equivalent thickness was managed by attaching multiple piezoelectric sheets together (bimorph actuator) along with the epoxy bonding layer. The results of the simulations are shown in Figure 4. To better compare the performance, all the conditions (i.e. radius, initial conditions, and boundary conditions) are the same in both groups of simulations. It is worth mentioning that here we focus on the amounts of the natural frequencies instead of the displacement amplitudes that are limited by the utilized modal analysis method.    Table 3 shows the comparison between the multisheet model and the optimized model. The absolute percentage error (APD) is calculated by using the following equation: where v F i is the natural frequency output from the Group-2 FEM simulation, and v opt i is the natural frequency of Group-1 also calculated by the FEM simulations.
As listed in Table 4, the natural frequency comparison between the analytically computed model and the FEM simulation for the optimized unimorph piezoelectric actuator above shows the acceptable accuracy of our proposed analytic model. We can see that this analytic model is very valuable, since the designers can use it to estimate the natural frequencies of a new system configuration with different geometric dimensions very efficiently and then realize it with the cost-effective offthe-shelf commercial PZT sheets without any computationally expensive FEM simulations. Moreover, the transverse vibration of the optimized PZT layer can be turned into a time-dependent vibration when the input voltage signal changes as a function of time. Thus, the designers may opt to determine the best input signal, which can maximize the efficiency of the system. In order to reach this goal, after calculating the natural frequencies, the first natural frequency will be chosen to apply to the system. A high voltage signal with the same frequency v 1 will be applied to the system and the maximum deflection will be measured and compared between the FEM simulation and experimental setup.
To verify the accuracy of the proposed model, we have built up an experimental measurement setup. A high voltage signal (i.e. 200sin(v 1 t)) with the same frequency v 1 was applied to the system, while the maximum deflection was reported from both the FEM simulation and the experimental measurement setup. The absolute values of the system response are shown in Figure 5. The maximum stress, which occurs at the edge, is shown as a function of time in this figure. Figure 5 also reveals that the maximum deflection of the system with respect to that excitation signal is equal to 121 mm as indicated by the reference scale located on the right.
The experimental measurement setup, as shown in Figure 6, includes a fixture, a high-precision laser measurement sensor (LK-H022, resolution of 0.02 mm), a signal generator, a high-precision high-voltage power amplifier (TEGAM-2350), and a data acquisition system. In order to verify the FEM simulation with reference to the experimental measurement results, an excitation signal was applied to the system. Figure 7 shows the center deflection of the double PZT sheets. As can be seen from this figure, the maximum value of the deflection is equal to 117 mm, which is very close to the FEM result of 121 mm with the error of 3.3% (i.e. (121 2 117)/121 = 3.3%).

Micropump design
In this section, the double-PZT-sheet actuator has been modeled onto a micropump to explore the effect of this optimized actuator on the flow rate of the micropump. A typical configuration of the micropump with the double-PZT-sheet actuator is modeled as shown in Figure 8. In order to design the microvalve, the pseudorigid-body model will be introduced. The purpose of the pseudo-rigid-body model is to provide a simple method for analyzing systems, which are under large nonlinear deflections. The pseudo rigid-body model concept is used to model the deflection of flexible segments by using rigid-body segments, which have the same force-deflection characteristics. 25 For each microflap, a pseudo rigid-body model can evaluate the deflection path of the flexible microflap. The motion will be modeled by rigid links attached at pin joints. Figure 9(a) shows the original and deflected position of microflap, while Figure 9(b) illustrates the pseudo rigid-body model. In this model, the spring is added in order to accurately predict the force-deflection relationship of the compliant microflap. The spring constant and deflection, and stress can be calculated as follows: where k is spring constant, E is the Young's module, I is the moment of inertia, and l is the length of microflap. r and k Y are the coefficients, which can be found from Howell. 25 where s max is maximum stress, a and b are the deflected position of the beam end point as shown in Figure 9(b), c is the distance from the neutral axis to the outer surface of the beam, Y is the pseudo rigid-body angle, n is a constant and can be found in Howell, 25 nP is the vertical component of the force, and A is the beam cross section area. The maximum pressure that can be produced by the double-PZT-sheet actuator under excitation signal of 200sin(v 1 t) is equal to 0:00216sin(v 1 t). By applying this pressure to the microflap, the maximum stress will be 1:01e 4 pa while the FEM simulation, as shown in Figure 9(c), reveals 1.22e 4 pa. As one major advantage, the pseudo rigid-body model makes it possible for the designers to analyze the compliant valve and height of channel without any computationally expensive FEM simulations.
For fabrication of microflaps and micropump, we propose using PDMS, which is able to provide unmatched advantages, including optical transparency, bubble permeability, quick and low-cost processing. PDMS, which typically mold and reversibly bond with very good elasticity, can be used for leak-free valving. In contrast, silicon is opaque, bubble and gas impermeable, and known for its expensive fabrication process. Another reason why PDMS is chosen in this study is its ability to form valves without leakage. Although the   silicon nitride or other inorganic materials have been used in the previous studies from the literature, they actually cannot provide reliable seal due to surface roughness and particles. However, the PDMS valves are able to provide more reliable sealing in presence of back pressure. Figure 10 illustrates the fabrication process of microfluid pump with one-way compliant microflaps. Fabrication starts with the first wafer substrate (step 1), followed by the deposition of a photosensitive resist substrate (step 2). The photoresist is photolithographically patterned (step 3) and developed. The PDMS will be first degassed and then is poured on top of the first wafer (step 4), cured, and then will be separated (step 5). The PDMS piece containing the inverse of the molded feature, which includes the micropump chamber and microflaps, is then attached to the second and third wafer to form enclosed chamber (step 6). Finally, the optimized double-PZT-sheet actuator will be placed on the top of the second substrate in step 6.

Experimental results
After designing the microflap, the double-PZT-sheet actuator micropump has been modeled in COMSOL Multiphysics as shown in Figure 11. The fully coupled fluid-structure interaction (FSI) COMSOL Multiphysics module is utilized to combine the fluid flow, which is formulated by using an Eulerian description and a spatial frame, with the solid mechanics, which is formulated by using a Lagrangian description and a material frame. 26 Figure 11(a) shows the micropump chamber in the suction mode, while Figure 11(b) exhibits the chamber of micropump in the pump mode. In those figures, the first and second color bars on the right show the velocity and stress, respectively.
In order to investigate the effect of the double-PZTsheet actuator on the system, the flow rate is calculated by FEM simulations. Figure 12 shows the flow rate of the micropump. As can be seen from this figure, the maximum flow rate reaches up to 22:9 3 10 À8 m 2 =s , which is suitable for pumping biological samples, such as blood. The result reveals that this flow rate is 14.7 times higher than a single layer PZT-sheet actuator. In addition, it is obvious from Figure 13 that with this new compliant microflap, the flow leakage will reach zero. The comparison of multiple conventional micropumps with our proposed device is summarized in Table 5. Our device is featured by its small footprint but a decent level of flow rate.

Conclusion
In this paper, we designed and optimized a PZT layer as an actuator in addition to microvalve for the micropump. The analytic analysis of the PZT-layer-based actuator was first conducted and then verified by using FEM simulations and experimental measurement. The comparison between the analytic computation and COMSOL Multiphysics simulation reveals the high accuracy of our proposed analytic solution, whose error is less than 2% for the first natural frequency. This is valuable since the designers can efficiently estimate the natural frequencies and mode shapes of a new configuration with different dimensions and further design the microflaps by using our proposed analytic solution without depending on any computationally expensive FEM simulations. Moreover, an experimental setup was established to measure the maximum displacement for a comparison with the FEM simulation. The results show that the difference between our experimental measurement and FEM simulation is less than 4%. Furthermore, the optimized actuator was applied to a typical micropump to evaluate its functionality. The comparison between a regular micropump with a single-PZT-sheet actuator and the same micropump with our optimized double-PZT-sheet actuator reveals that the flow rate could be increased up to 14.7       28 PZT Active F7:5 mm 32 mbar 3.5 mL/min Water Ma et al. 29 PZT Active F10 mm -9.1 mL/min Water Le and Hegab 30 PZT Passive F20 mm 2x 18 kpa 32 cm/min Gas Cheng et al. 31 PZT Active -2.35 kpa 0.24 mL/min Air Dong et al. 32 PZT