Dependence of the Piezoelectric Micropump Operating Mode on Its Geometry

A mathematical and computer axisymmetric models in which periodic oscillations of ring piezoelectric actuators placed on an elastic tube of circular cross-section leads to radial deformations of the tube is constructed. The result of these displacements is a local compression of the microchannel, which leads to changes in its volume and the corresponding pushing of the fluid. If a non-symmetrical oscillation scheme is used, the average fluid flow over the period will be non-zero and the device can be used as a micropump. The aim of the work is a computer study of the geometric features of an axisymmetric piezoelectric micropump, taking into account the processes of heat transfer by a fluid in the channel. The dependence of the average fluid flow rate on the channel parameters and the frequency of oscillations of the piezoelectric actuators is determined by the method of factorial experiment. The parameters preventing heat backflow have been determinate, which makes it possible to use a device for supplying coolant to a microgripper cooling system.


Introduction
Mathematical and computer simulations of microfluidic systems have recently been actively evolving in connection with the development of new types of microdevices [1]. The designs of new technical devices become possible only due to strict consideration of the peculiarities of fluid flow in microsystems. As the scale of the Reynolds Re and Péclet Pe numbers system decreases, the dynamic properties of the system become small. The small number Re implies the laminar nature of the flow pattern, corresponding to both most technical microdevices and biological systems [2]. Consequently, it becomes possible to develop efficient fluid delivery and dosing systems [3], as well as miniature cooling systems [4]. The small number of Pe means that the contribution of molecular thermal conductivity in heat transfer prevails over convective heat transfer, so that both transfer mechanisms [5] must be considered in micro-devices (microheat exchanger, cooling system, microgripper). The basic characteristics of microdevices largely depend on the operating principles and design of micropumps, as well as hydraulic resistances. The most common are peristaltic pumps in which only a flexible tube is in contact with the working fluid [6], thereby insulating the fluid and preventing the ingress of foreign impurities.
The intensive growth in the production of microelectromechanical systems [7] and the widespread use of microrobots [8] provide a strong incentive for the design of new promising devices for careful micro-components manipulation [9]. One of that devices is a capillary microgripper (MG) for manipulating flat microobjects and membranes with a sufficient area of the flat side, the retention of which is impossible with finger grips [10,11].
The creation of a model of this MG, built based on the miniature Peltier element, involves the use of a cooling system to remove heat from the hot surface of the Peltier element. Afshari studies [12] have shown that with a decrease in the size of the chamber, the efficiency of the fluid cooling system increases in comparison with the air-cooling system. The removal of the heat power from the Peltier element requires the design of a miniature pump, and the development of a mathematical and computer model to ensure the operation of the MG mode.
Existing pumps are being developed for a specific application [13]. It is difficult to find a micropump for a MG cooling system with a small size, a sufficient fluid flow rate and a flexible fluid flow control system. In the design of new technical devices, this requires the formulation of new multiphysical mathematical models which unite the hydrodynamics, elasticity theory, heat transfer, parametric optimization and device control theory [14].
With the development of mathematical and computer models, numerical simulation software using the finite element method (FEM) have appeared [15], which allow simulating devices with almost arbitrary geometry [16]. In the paper [17] an axisymmetric computer model for the interaction of a Newtonian fluid with a hyperelastic incompressible body was developed and FEM modeling to analyze the stability of differential equations in the FreeFem++ software [18] was used.
The authors of this paper propose the idea of a piezoelectric micropump (PEMP) consisting of an elastic tube with a system of ring piezoelectric actuators (PEAs) placed on it. New algorithms have been developed for modeling fluid flow in a channel with dynamically variable geometry [19]. In [20] the deformation of the tube was calculated by the elastic equations on the pipe wall. In the subsequent work [21], a mathematical and computer models of the fluid flow regulator were presented using the dynamic formation of hydraulic resistance by compressing an elastic tube with PEA. A model of a PEMP that creates a fluid flow (flow rate up to 50 μl/s) in a flat channel using a bending PEA immersed in it was developed [22]. And in the paper [23], an axisymmetric computer model of a PEMP was proposed. The model can be used both in the development of the design of a PEMP operating in a closed circuit, and in the synthesis of software for the MG fluid cooling control system.
A complete numerical simulation of the device using a multi-parameter set of input data is required to determine the operating mode of the technical device. In the 2D model of the MG fluid cooling system [24], the relationship between the fluid flow rate and the heat transfer coefficient was analyzed, and the influence of flow pulsations on cooling efficiency was studied by the orthogonal central compositional design (OCCD). The result of a factorial computational experiment is an analytical function that can be used in the software of the device control system in real time.
In this paper, we analyze the geometric features of an axisymmetric PEMP [23] considering the processes of heat and mass transfer.

Mathematical model
An axisymmetric fluid flow through an elastic microtube with the internal R1 and external R2 radii and the length L is considered (figure 1). The central part of the microtube is located inside a system of ring PEAs with a length of ℓ. The cylindrical coordinate system contains the radial Or and central Oz axes of the tube with the origin located in its geometric center.
The operating principle of the PEMP is as follows. At the first stage, all PEAs are radially compressed. At the next stage, all PEAs are successively returned to their original state, starting from the leftmost one. Therefore, during the period of the compression and recovery phases, an average non-zero flow is generated in the positive direction (Oz).
The fluid flow is described by the Navier-Stokes equations and the incompressibility equation in axisymmetric geometry [16]:  Figure 1. The PEMP geometry: Γ1 and Γ2 are the inlet and outlet openings; Γ3 is the axis of symmetry; Γ4, Γ5 and Γ6 are the inner, side and outer walls; are the contact of the i-th PEA with Γ6. The extended geometry (dotted lines) for the heat problem: Γ7 is the thermostat-tube boundary; Γ0 and Γ8 are the inlet and outlet openings.
where ρ is the density of the fluid; u = (ur, uz) is the velocity vector of the fluid consisting of the radial r and axial z components; t is the time; is the nabla operator; p is the pressure; µ is the coefficient of dynamic viscosity of the fluid.
The deformations of the elastic microchannel are described by the Navier-Cauchy equations [15]: (2) where s = (sr, sz) is the strain vector; λ and η are the Lame coefficients associated with the Young's elasticity modulus E and the Poisson's ratio ν; I is the unit tensor. The investigated model of the PEMP is supposed to be used in the MG cooling system. Therefore, to check the degree of influence of the heated fluid from the MG on the cooling fluid of the microchannel, the heat equation is solved [25]: where T is the temperature of fluid; cp is the isobaric specific heat of the fluid; λT is the coefficient of thermal conductivity.
Equations (1)-(3) were solved numerically by the FEM in the numerical simulation software FreeFem++ [18]. The time discretization for the Navier-Stokes and heat equations were performed by the implicit Euler method of the first order. The initial conditions correspond to the absence of deformations on the outer wall s = 0 and the resting fluid u = 0.
The boundary conditions for the fluid velocity have the following form: where Δp is the pressure difference at the inlet and outlet openings of the tube. The microtube is fixed at the ends (at the boundary Γ5): sr = 0, sz = 0. The Neumann boundary conditions for deformations are determined from the equilibrium condition σrr = -pp,i(t), where the stress tensor component σrr is balanced by the pressure -pp,i induced on the tube by the i-th PEA at the boundary : The dependence of the PEA pressure on the microchannel outer wall is defined in the form pp(t) = max(pp)fp(t), where the signal form fp(t) = (1 -cos(2πft)) / 2 is a normalized harmonic function in the region [0, 1] with the oscillation frequency f.
Boundary conditions for the heat equation are T = T0 on Γ1 and T = T1 on Γ7, where T0 is the ambient temperature; T1 is the temperature of the thermostat TS.

Results
The following values were determined in computer simulations: the ratio of the pipe's outer to inner surfaces deformations in contact with the PEA -; (5) the average for the period steady-state flow rate - where Vj and Vj-1 are the volumes of fluid flowing through the tube at the end of the j-th and (j -1)-th periods of oscillation of the PEAs. The time of establishing the flow τu, for which Qa reaches a constant value, is determined by the formula of the characteristic relaxation time of the velocity: Computational experiments were performed for the following parameters: the length of the tube L = 10 mm; the number of PEAs np = 5; the width of the PEA Lp = 1 mm; the maximum radial compression of the PEA sm = 5 mkm. Parametric analysis was carried out using OCCD [26] for the following parameters: inner radius of the tube R1 = 1, ..., 1.5 mm (X1); the outer radius of the tube R2 = 1.75, ..., 2.25 mm (X2); the oscillation frequency of the PEAs f = 1, ..., 2 kHz (X3). The parameters Xi (i = 1, 2, 3) were normalized to the interval [-1, 1] using the following transformation: where xi is the normalized parameter; Xi,0 is the central point of the interval; Xi,ℓ is the maximum deviation of the parameter from the central point.
The matrix of the computational experiment was constructed and the OCCD response variables were determined in the form of second-degree polynomials. For the ratio (5):   (2) and then the coefficient k was calculated using the formula (5). The analysis of the (8) leads to the following conclusions. The coefficient k increases with an increase in R1 (a1 > 0), a decrease in R2 (a2 < 0), and with a proportional increase in R1 and R2 (a3 > 0). This is due to the fact that as the tube wall thickness increases, the  2  2  2  1 2 3  0  1 1  2 2  3 3  4 1 2  5 1 3  6 2 3  7 1  3  8 2  3  9 3  deformations on the inner surface are redistributed wider along the Oz axis. The dependence of k(R1) is nonlinear -with the growth of R1, the increase in the contribution to k slows down (a4). The coefficient k is practically independent of (a5 ≈ 0). The solution of the related hydrodynamic problem and the elasticity problem was determined numerically using equations (1) and (2). To ensure that the average fluid flow rate Qa (6) reaches the saturation, the time interval for calculation was chosen to be equal to 1.2τu (7). An analysis of the (9) showed that the leading linear contribution to the average flow rate of the fluid Qa is made by the oscillation frequency of the PEAs f (x3, figure 2, lines 1-4 and 5-8) and the inner radius of the tube R1 (x1, figure 2, lines {1, 3, 5, 7} and {2, 4, 6, 8}). With an increase in these parameters, the fluid flow rate increases (b3 > 0 and b1 > 0), their couple interaction is positive (b5 > 0). With an increase in the outer radius of the tube R2 (x2), the fluid flow rate Qa decreases (b2, figure 2, lines 1 and 3, 2 and 4, 5 and 7, 6 and 8) due to a decrease in the ratio k. With a proportional increase in R1 and R2, the coupled contribution is negative (b4 < 0) due to a decrease in the ratio sm/R1 (a decrease in the hydrodynamic resistance). The quadratic terms make a small positive contribution to Qa (b7 > 0, b8 > 0, b9 > 0).
For the possibility of using this PEMP for fluid cooling of the MG, we will estimate the average flow rate of the fluid, which will ensure that the heat propagation from the right opening of the tube is blocked (see figure 1). In the flow absence, the heat equation and the boundary conditions have the form with a solution corresponding to a linear temperature distribution (10) The condition for blocking the heat propagation is a negligible temperature value near the rightmost PEA T(np·Lp/2) ≈ 0, so, there is no flow of the heated fluid in the opposite direction. Figure 3 shows the deviation of the fluid temperature T from the ambient temperature T0 (ΔT = T -T0) on the tube axis r = 0 at the z coordinate zT = 3 mm versus the time t for different values of the average fluid flow rate Qa. The line 1 corresponds to the absence of fluid flow (over time, the temperature reaches a stationary value of T0 + 3.5 °C in accordance with the formula (10)). In the presence of a flow (lines 2-5), with an increase in Qa, the temperature value at the point zT decreases. For Qa > 800 nl/s, the temperature value is almost zero, thus there is no heat backflow. From a set of computational experiments, it can be seen that heat backflow to the z < zT area is absent at R1 = 1 mm and f > 1.5 kHz or at R1 > 1.25 mm and f = 1 kHz.   Figure 4 shows the temperature distributions inside the channel at time tT = 120 s after the start of the process. With an increase in the fluid flow rate, the time to reach the stationary temperature distribution τT decreases. For Q < 300 nl/s, the time tT is not sufficient to reach the stationary mode. For Qa = 300 nl/s, the time is τT ≈ 120 s, for Qa = 500 nl/s -τT ≈ 40 s, for Qa = 1 μl/s -τT ≈ 15 s.

Conclusion
The paper proposes mathematical and computer models of a PEMP. The parametric analysis of the operating mode of a PEMP depending on the internal and external radii of the pipe and the frequency of oscillations of the PEAs was fulfilled in framework of the OCCD. The efficiency of the pump increases with an increase in R1, f and with a decrease in the relative wall thickness of the tube R2/R1. When design a prototype of a PEMP, the OCCD approximation will allow one to determines the parameters of the device to ensure the required average fluid flow. When using a PEMP to supply fluid to the MG cooling chamber, it is necessary to ensure a fluid flow rate Qa > 800 nl/s, which prevents heat backflow. Based on simulation results there is no backflow with a system of 5 PEAs providing a radial compression value of 5 mkm at an oscillation frequency f > 1.5 kHz on a tube with an internal diameter R1 > 1.25 mm.