Design of the micropump and mass-transfer compartment of a microfluidic system for regular nonenzymatic glucose measurement.

The aim of this paper is to design and numerically simulate the mass-transfer compartment and piezoelectric micropump of an implantable integrated microfluidic device for regular microdialysis-based nonenzymatic measurement of glucose level in diabetic patients. The device function is based on the process that the piezoelectric micropump pumps the dialysis fluid into the mass-transfer compartment microchannels, where the interstitial fluid (ISF) glucose diffusion into this dialysis fluid gives it a glucose content, then detected and measured in the sensor section. This diffusion takes place through the semipermeable membranes located in the microchannels at the base of the hollow microneedles entering the body skin painlessly. The value of dialysis fluid flow rate (1 μL/min) was chosen so that the best achievable recovery factor can be obtained while the size and time delay of system were being kept at the best minimum possible. In the mass-transfer compartment, the number of microneedles, the dimensions of microchannels and the thickness of membranes were selected so as to achieve the best appropriate recovery factor, minimum possible size as well as considering the fabrication feasibility. Furthermore, in the different parts of micropump, the materials and dimensions were chosen so as to provide the needed flow rate with the best minimum voltage, sufficiently small size and fabrication feasibility.

Introduction

In diabetic patients, the body capability of glucose metabolization abnormally decreases. More precisely, these patients suffer from excessively increased levels of blood sugar (hyperglycemia), which, if not treated, will develop defective vision, brain strokes and cardiovascular complications, to name a few. On the other side, using insulin injection or pills that induce insulin secretion are likely to result in an excessive decrease in the blood sugar level (hypoglycemia) in diabetic patients. As the normal performance of the brain heavily relies on glucose, hypoglycemia proves a real threat, as it may lead to dysfunction of the patient's brain or even put them into coma, in extreme cases [1], [2], [3], [4], [5]. Hence, diabetic patients need a regular measurement system to check their glucose level, which can be painlessly implanted.

To continuously detect glucose, two major types of sensors can be utilized, namely needle-type and microdialysis-using sensors [6,7]. In the microdialysis method, while the sensor is placed on the exterior surface of the body, the dialysis probe that contains the semipermeable membrane is put inside the body. When the glucose-free dialysis fluid is pumped into the probe a process starts in the following order: the glucose molecules inside the interstitial fluid (ISF) diffuse through the semipermeable membrane, and enter the dialysis fluid, which becomes glucose-containing before being sent into the sensor [8]. While this system needs a greater space, because the sensor is put outside the body, it is biologically safer and diagnostically more accurate than the needle-type sensor [7,9].

From the chemical viewpoint, glucose can be electrochemically detected either with an enzymatic or with a nonenzymatic method. While the enzymatic sensors suffer from difficulty of enzyme immobilization, low chemical and thermal stability and activity changes following the changes in environmental conditions, nonenzymatic sensors, which can be more easily manufactured, show higher chemical and thermal stability [10,11].

The commercial microdialysis-using glucose sensors, which are mainly made by Menarini Company, are of enzymatic type and suffer from some major weak points, including their inappropriately large size and deep penetration of the dialysis probe. Besides, this probe should be put under the patient's skin by an experienced medical expert [12,13].

Our main goal is to design an integrated implantable microfluidic system according to the microdialysis method in order to regularly and painlessly measure glucose in diabetic patients in a nonenzymatic way. In our previous work [10], we designed and fabricated a promising novel microfluidic nonenzymatic glucose sensor. In following the very work, here we address the other parts of the entire system. Fig. 1 schematically gives an overall view of how the device will work, while Fig. 2A and B show from two different views the main compartments of the integrated system, which are as follows: The piezoelectric micropump, the compartment for mass transfer which contains the array of hollow microneedles, and the electrochemical amperometric glucose sensor. By inserting the microneedles into the skin, they are filled with the ISF due to surface tension. This ISF is separated from the dialysis fluid by the semipermeable membranes located at the microneedles base. Passing over the membrane surfaces, the pumped dialysis fluid takes in, via the diffusion process, the ISF glucose, which is then measured within the sensor section. Beside the sensor section, the other two vital parts of this device are the mass-transfer compartment and the piezoelectric micropump. Therefore, their designing and numerical simulation will be the focus of the present paper.

Fig. 1

The schematic view of the flow diagram of our intended device.

Fig. 2

Two different views of our integrated device for regular glucose measurement [10].

In our system, the glucose measurement is done by means of ISF, which has proven itself as a reliable source of biomarkers which can substitute blood plasma [14], [15], [16], [17], to which it has shown high similarity, based on analyses [14,15]. In the recent years, much interest has been put in microneedles for being used for minimally invasive drug delivery, ISF sampling and diagnostics [18], [19], [20], [21], [22], [23]. For painless use of the microneedles, they should penetrate the skin to 150–200 μm depth [24,25]. A duly widely used material in microneedle fabrication is silicon [26], whose flexible nature makes it possible to produce different shapes and sizes of the needles [27]. Silicon substrates are great candidates for precise manufacture and batch production [27]. Despite some concerns about silicon biocompatibility, silicon devices, including the microneedles, have been announced clear by the Food and Drug Administration (FDA) [26]. In their research, Mukerjee et al. fabricating silicon microneedles of different profiles, concluded that the microneedles with the "snake fang" profile are the best for extracting the ISF, as they can provide the necessary mechanical strength for penetrating the skin without the bore hole plugging problems [28].

As previously explained, in the microdialysis method, the dialysis fluid is pumped into the system. The diaphragm-pump, in which the diaphragm makes a reciprocating movement like a piston, is one of the most commonly used mechanical pumps at microscale [29]. The diaphragm can be driven by different actuation mechanisms including the electrostatic, piezoelectric, electromagnetic, thermo-pneumatic, bimetallic and shape memory alloy [29], [30], [31], [32], [33], [34], [35], [36], [37], [38], [39], [40]; among them the piezoelectric mechanism is widely used because of the piezoelectric micropump fast response and high displacement capacity [29,30]. Zinc oxide and PZT (lead zirconate titanate) are the most widely used piezoelectric materials [41]. On the other hand, the micropumps can be further categorized into valveless and check-valve micropumps [42]. Although the valveless micropums can be easily fabricated and there is not any risk of valves fatigue, they encounter serious problems of flow control difficulty, fluid reflux and energy loss [35,42]. The check-valve micropumps can be categorized based on having active or passive valves. Since the passive valves are far simpler and, unlike the active valves, do not need any starter and input signal, they are more preferable.

In the present paper, Sections 2 and 3 deal with the design and numerical simulation of the mass-transfer compartment for an implantable integrated microfluidic device for regular glucose measurement, while Sections 4 and 5 address the design and numerical simulation of the piezoelectrically-activated diaphragm-type micropump with flapper valves for the very microfluidic device. The factors considered in the designs include the recovery factor, size, time delay, and fabrication feasibility.

Materials and methods (for mass-transfer compartment)

Design and simulation of the mass-transfer compartment

The schematic performance of the mass-transfer compartment is given by Fig. 3A, showing the cross-section of a microneedle above which the dialysis fluid is passing through the microchannel. Following the microneedle entrance into the skin, the surface tension makes the ISF move up the internal walls of the microneedle to the semipermeable membrane at the base of the microneedle. This membrane, which separates the ISF from the dialysis fluid, allows small molecules such as glucose ones diffuse into the dialysis fluid, while it does not let in large molecules such as those of proteins, which can contaminate the working electrode surface of the sensor.

Fig. 3

A: The schematic performance of the mass-transfer compartment. B: The considered profile for each microneedle.

Since there is only a significantly small amount of glucose diffusion through one microneedle, there should be used an array of microneedles each of which has a "snake fang" profile as is shown in Fig. 3B. Each microneedle is 150 μm high [26], which is shorter than the distance of the nerve tips from the skin surface, which means the entrance of the microneedles causes no pain for the patients.

The "snake fang" profile of the microneedle prevents the ISF tissue from being compressed at the entrance of the microneedle, which makes this profile the best one for extracting the ISF [28]. Therefore, when a patient inserts those microneedles into her/his skin, the ISF rises inside the microneedles to reach the semipermeable membranes. The center-to-center distance between two adjacent microneedles is 300 μm [26]. There was selected an array of 24 × 12 microneedles considering the size, recovery factor and the time delay. Fig. 4A shows the microchannel passage system for the dialysis fluid flow from the micropump (not shown here) to the mass-transfer compartment and then to the sensor. In practice, these microchannels can be created by means of the photolithography method on a plate made of pyrex or silicon.

Fig. 4

A: The microchannel passage system for the dialysis fluid flow. B: The zoomed-in configuration of the main microchannel leading to the smaller microchannels.

As can be seen in Fig. 4A, the dialysis outflow from the micropump is divided into two parts and then, by means of these two main microchannels, connected to the smaller microchannels passing over the microneedles. It is important to keep the fluid flow rates in the smaller microchannels rather the same as far as possible. Dividing the outflow of the micropump into two main microchannels instead of one before connecting to the smaller microchannels makes it easier to keep the flow rates in the smaller microchannels the same. The number of smaller microchannel rows was chosen even for being equally divided between the two main microchannels. Since the mass-transfer compartment is connected to the end of the micropump, the width of the microchannels-containing plate was selected the same as that of the micropump (1 cm). Here, considering 24 rows of smaller microchannels can be a good choice because of the appropriate use of the total plate width and approximately 1.5 mm distance between the plate edges and each of the first and final rows, which is convenient.

Fig. 4B shows the zoomed-in configuration of the main microchannel leading to the smaller microchannels. As can be seen, the widths of these input channels were selected different so that there can be the same fluid flow rate in the rows. Table 1 shows these width values.

Table 1

The width of the input microchannels from left to right.

Microchannel number	1	2	3	4	5	6	7	8	9	10	11	12
input width (μm)	137.30	106.50	76.05	51.10	35.50	22.70	12.50	7.00	5.05	6.20	15.30	149.80

The height of the channels was chosen 60 μm, since the height of the sensor to which these channels carry the dialysis fluid is 60 μm [10], so that there would be no height difference. Since the inner diameter of each microneedle is 30 μm [26], there was chosen a 50 μm width for each microchannel located above the microneedles in order to provide enough margin for covering the inner diameter of microneedles. It is worth mentioning, the closer to 30 μm the widths of the microchannels are, the higher the recovery factor will be. Therefore, it is not recommended to consider the widths over 50 μm.

The numerical simulation was performed using the finite element method (FEM). The related differential equations were solved by means of a two-step segregated solver: the first step for velocity and pressure and the second step for concentration. The linearization of the equations was done by the Newton method and were solved by the use of Generalized Minimal Residual (GMRES) algorithm. For the time-dependent case, the temporal discretization was done using the implicit second order backward differential formula (BDF). First order elements were used for the discretization of pressure, velocity and concentration.

The governing equations

The flow of the dialysis fluid is a 3-D incompressible one, which is governed by the continuity and Navier-Stokes equations as given in Eq. (1) and Eq. (2), respectively.where , , and are the velocity, fluid density, time, pressure, and dynamic viscosity, respectively.

By numerical solving of Eq. (3) one can obtain the concentration distribution, with the assumption that the diffusion is the only process for the mass transfer through the membranes.where and are the concentration and diffusion coefficient, respectively. According to Zahn et al. [43], equals in the dialysis fluid and in the semipermeable membranes.

The glucose mass flux is the same on the semipermeable membrane-dialysis fluid boundary surfaces, where the concentration is yielded by Eq. (4).in which is the partition coefficient considered unity here (Additional information on the partition coefficient is provided in the Supplementary).

The recovery factor (RF) is obtained from Eq. (5).in which the numerator is the glucose concentration in the dialysis fluid that leaves the mass-transfer compartment and the denominator is the glucose concentration in ISF.

Boundary conditions (BCs)

Since the flow rates of the fluid pumped into each of the identical smaller microchannels are the same, it is sufficient to take one row of these channels for the numerical simulation.

For the inlet surface of the microchannel, the constant flow rate (which is the total flow rate divided by 24) is defined as the boundary condition. For all the channel walls and the membranes surfaces, which are in contact with the fluid flow, the no-slip boundary condition is considered and for the outlet surface of the microchannel, the outflow rate boundary condition is defined. In order to solve the concentration equation, the boundary condition of specified ISF glucose concentration is applied for those semipermeable membrane surfaces in contact with the ISF. On the boundaries between the semipermeable membranes and the dialysis fluid, the glucose mass flux is the same. The glucose concentration in the dialysis fluid and its gradient are zero respectively for the inlet and outlet surfaces of the microchannel. On the other boundary surfaces the glucose mass flux equals zero.

Results and discussion (for mass-transfer compartment)

To give an overall picture of the microchannel meshing system, Fig. 5 shows a part of the meshed microchannel. The microchannel and the semipermeable membranes are respectively meshed with tetrahedral and planar triangular elements. The semipermeable membranes are those parts which are meshed with finer elements.

Fig. 5

A part of the meshed microchannel.

Effects of the total dialysis fluid flow rate and channel height

One of the most important issues here is the selection of an appropriate total dialysis fluid flow rate, for which it is necessary to investigate the effect of this flow rate on the recovery factor as one of the most important elements. In order to conduct this investigation, the steady-state average glucose concentrations on the specified cross-sections of the microchannel (shown in Fig. 6A, cross-sections after every 3 consecutive membranes), were calculated for different flow rates.

Fig. 6

A: The cross-sections after every 3 consecutive membranes. B: The cross-sections after each membrane.

As can be seen in Fig. 7A, the ratio of the point-wise average glucose concentrations (identified on the specified cross-sections shown in Fig. 6A) to the constant ISF glucose concentration are plotted versus the number of the microneedles having been passed by the fluid, for 15 nm-thick semipermeable membranes.

Fig. 7

A: The ratio of the point-wise average glucose concentrations (identified on the specified cross-sections shown in Fig. 6A) to the ISF glucose concentration versus the number of the microneedles having been passed by the fluid, for different flow rates. B: The ratio of the point-wise average glucose concentrations (identified on the specified cross-sections shown in Fig. 6A) to the ISF glucose concentration versus the number of the microneedles having been passed by the fluid, for different channel heights. C: The ratio of the point-wise average glucose concentrations (identified on the specified cross-sections shown in Fig. 6B) to the ISF glucose concentration versus the number of the microneedles having been passed by the fluid. D: The ratio of the point-wise average glucose concentrations (identified on the specified cross-sections shown in Fig. 6A) to the ISF glucose concentration versus the number of the microneedles having been passed by the fluid, for different membrane thicknesses.

Fig. 7A shows that the recovery factor increases when the dialysis fluid flow rate decreases. For any specified channel height, when the fluid flow rate increases, although the velocity will increase, the recovery factor decreases, because a greater volume of fluid per unit time flows into the channel. On the other hand, for any specified fluid flow rate, whatever the decrease is in the microchannel height, the velocity increase will lead to an increase in the recovery factor. In Fig. 7B, for the total fluid flow rate of 1 μL/min and channel heights of 60, 80 and 100 μm, there are shown the ratio of the average steady-state glucose concentrations (identified on the specified cross-sections shown in Fig. 6A) to the constant ISF glucose concentration versus the number of the microneedles having been passed by the fluid (for 15 nm-thick semipermeable membranes).

As can be seen in Fig. 7B, the 60 μm height of microchannel results in a recovery factor higher than those for 80 μm and 100 μm. The reason is that, at a specified fluid flow rate, when the microchannel height decreases, the velocity increases, as a result of which the diffused glucose molecules inside the microchannel are removed faster from semipermeable membrane surfaces, which means, a higher rate of glucose diffusion into the microchannel through the semipermeable membranes.

Here, 1 μL/min was chosen as a highly appropriate flow rate. Obviously, the greater the recovery factor, the more accurate the system is. In fact, depending on the geometry of the system and the thickness of the semipermeable membrane, this recovery factor can take a wide range of values, even as low as 10% [44]. Of course, if the flow rate is selected less than 1 μL/min, although the recovery factor is increased, there is also an increase in the time delay. Furthermore, in lower glucose concentrations, depending on the number of active sites on the sensor working electrode surface, there usually is a higher electrochemical sensor sensitivity [45], which is another reason why the fluid flow rate was not chosen lower than 1 μL/min. On the other side, it is preferable for the fluid flow rate not to exceed 1 μL/min, because for example when it is increased to 2 μL/min, it leads to a 21% decrease in the recovery factor, the compensation for which requires an increase in the mass-transfer compartment size. Another problem due to a greater flow rate pumping is the necessity of increased power and, therefore, increased micropump size.

Since 12 microneedles in each microchannel yield a desirable recovery factor (72.2%), there is no need to increase this number, which will undesirably increase the system size and time delay.

As can be seen in Fig. 7C, the ratio of the point-wise average glucose concentrations (identified on the specified cross-sections shown in Fig. 6B, the cross-sections after each membrane) to the constant ISF glucose concentration are plotted versus the number of the microneedles having been passed by the fluid (for the flow rate of 1 μL/min and 15 nm-thick semipermeable membranes).

As shown in Fig. 7C, the longer the path along the microchannel paved by the fluid, the lower the rate of increase in the glucose concentration. The reason is that the closer is the channel outlet, the higher is the glucose concentration in the dialysis fluid, which means the concentration between this fluid and the ISF decreases, which in turn results in a decrease in the glucose concentration flux into the microchannel. If the number of microneedles in each microchannel approaches a sufficiently large number, the recovery factor tends to one.

Effect of the semipermeable membrane thickness

As can be seen in Fig. 7D, the ratio of the point-wise average glucose concentrations (identified on the specified cross-sections shown in Fig. 6A) to the constant ISF glucose concentration are plotted versus the number of the microneedles having been passed by the fluid, for different membrane thicknesses.

As can be seen, an increase in the membrane thickness leads to a higher resistivity against the mass transfer, which decreases the recovery factor. Since it is feasible to fabricate a 15 nm-thick polysilicon semipermeable membrane [46], this value was chosen as the membrane thickness.

The unsteady-state solution

Fig. 8B shows the average glucose concentration versus time on the outlet surface of the microchannel for an assumed constant ISF glucose concentration of 5 mM (Fig. 8A). As can be seen, after approximately 17 s, the glucose concentration reaches its steady-state value. Fig. 8D shows the average glucose concentration versus time on the outlet surface of the microchannel for an ISF concentration step change shown in Fig. 8C. According to Fig. 8C, the ISF glucose concentration is 7 mM for the first 26.5 s, which goes up to 8 mM within 26.5–27.5 s time interval to remain constant hereafter. Fig. 8D shows it takes about 17 s for the outlet glucose concentration to become 72% (the RF) of the 7 mM, remaining constant until 26.5 s, when the ISF glucose concentration changes. Following the next approximately 17 s, the outlet concentration gets to 72% of the 8 mM.

Fig. 8

A: ISF glucose concentration versus time. B: The average glucose concentration versus time on the outlet surface of the microchannel for the ISF glucose concentration of A. C: ISF glucose concentration versus time. D: The average glucose concentration versus time on the outlet surface of the microchannel for the ISF glucose concentration of C.

It is worth mentioning that the outflow rate of the micropump is fluctuating. However, since the mass-transfer compartment simulation with the fluctuating flow rate (at the frequency of 2 Hz and average of 1 μL/min) showed the same recovery factor and response time as the case of simulation with constant 1 μL/min flow rate, all the simulations in this section were done with constant flow rate, for the sake of decreasing run time.

Verification and grid independence study

In order to validate the simulation model, we numerically simulated the hollow fiber membrane contacting system of the experimental work by Ghobadi et al. [47], since the underlying physical principles of mass transfer in their work are highly similar to ours. We checked our results with their experimental ones. They have investigated the deletion of CO2 and NO2 entering the contactor from the shell-side. NaOH absorbent enters counter currently from the tube side. Through the membrane, CO2 and NO2 diffuse into the tube, where they react with NaOH. In the tube and shell, mass transfer mechanisms are advection and diffusion, while in the membrane it is only diffusion. The simulation parameters are shown in Table 2.

Table 2

The simulation parameters for validation.

r_1(mm)	inner fiber radius	0.34	D_CO2,s(m^2s)	Diffusivity of CO_2 in shell	1.855×10^−5
r_2(mm)	outer fiber radius	0.6	D_CO2,t(m^2s)	Diffusivity of CO_2 in tube	1.92×10^−9
L (mm)	Module effective length	200	D_CO2,m(m^2s)	Diffusivity of CO_2 in membrane	2.29×10^−6
r_3(mm)	Radius of free surface (based on Happel's free surface method)	1.04	D_NO2,s(m^2s)	Diffusivity of NO_2 in shell	1.54×10^−5
K_r_NO2(m^3mol.s)	NO_2-NaOH reaction rate constant	10^5	D_NO2,t(m^2s)	Diffusivity of NO_2 in tube	1.4×10^−9
K_r_CO2(m^3mol.s)	CO_2-NaOH reaction rate constant	8.37	D_NO2,m(m^2s)	Diffusivity of NO_2 in membrane	1.9×10^−6
H_CO2(kmolkg.atm)	Henry's constant for CO_2	0.034	D_NaOH(m^2s)	Diffusivity of NaOH in tube	9.6×10^−10
H_NO2(kmolkg.atm)	Henry's constant for NO_2	0.007

In their work, Ghobadi et al. experimentally obtained 85% NO2 removal and 73% CO2 removal for the mean flow velocities of 1.05m/s and 0.05m/s in the shell and tube, respectively (for 0.5M NaOH and 2% CO2 and 300 ppm NO2 balanced in N2). Our corresponding numerical simulation results are the removal of 81% for NO2 and 72.3% for CO2. The slight corresponding differences show the satisfactory accuracy of our numerical model. These differences stem from the natural differences between the experimental work and the numerical simulations. For example, these differences can be due to the used appropriate approximate simulation parameters such as the diffusion coefficients in the numerical simulation and the actual ones in the experimental work.

For grid independence study, since the element types used for the meshing of the semipermeable membranes and of the microchannel (Fig. 5) were different, the study was conducted separately, meaning that it was first carried out for the membranes with a specified microchannel element number, the results of which were used for grid independence study of the microchannel. This procedure was repeated again (using the findings of the first studies) and the results are shown in Fig. 9A and B. For the membranes grid independence study, Fig. 9A shows the average glucose concentrations for each of the first five membranes circular surfaces in contact with the dialysis fluid, for different numbers of elements for each membrane (for 1 μL/min total flow rate, 15 nm membrane thickness and 5 mM ISF glucose concentration).

Fig. 9

A: The average glucose concentrations for each of the first five membranes circular surfaces in contact with the dialysis fluid for different numbers of elements for each membrane. B: The point-wise average glucose concentrations, identified on the specified cross-sections shown in Fig. 6A, versus the number of the microneedles having been passed by the fluid for different numbers of elements for the microchannel.

As can be seen, for each membrane, for the numbers of elements equal to or greater than 1808, the results overlap with each other with high accuracy. So, the number of elements equal to 1808 was chosen for meshing of each membrane.

Also, to show the independence of the results from the microchannel numbers of elements, the average glucose concentrations on the specified cross-sections (cross-sections shown in Fig. 6A) versus the number of the microneedles having been passed by the fluid are depicted in Fig. 9B for different microchannel numbers of elements (for 1 μL/min total flow rate, 15 nm membrane thickness and 5 mM ISF glucose concentration).

As can be seen in Fig. 9B, the results for the mentioned numbers of elements overlap with each other accurately. Since the results of 605730 and 995900 microchannel numbers of elements coincide with each other with very high accuracy, 605730 cells were chosen for meshing the microchannel.

Materials and methods (for micropump)

Design and numerical simulation of the micropump

The designed micropump has a main chamber above which on the diaphragm there is attached a piezoelectric actuator. There are also two flapper check valves on both sides of the micropump. Fig. 10A shows the designed micropump. As can be seen, the different layers, made of different materials are shown in different colors. Fig. 10B shows the micropump stack up layers, which from the top to the bottom, are as follows: 1- piezoelectric disk 2- diaphragm 3- silicon layer containing channels and chamber 4- PDMS layer containing flappers 5- silicon bottom layer. Bonded together the layers make up the micropump. The length, width and thickness of the micropump (Fig. 10A) are 10 mm, 10 mm and 0.65 mm, respectively.

Fig. 10

A: The designed micropump. B: The micropump stack up layers, 1- piezoelectric disk 2- diaphragm 3- silicon layer containing channels and chamber 4- PDMS layer containing flappers 5- silicon bottom layer.

Since among the common materials for micropump fabrication are glass and silicon and the latter was chosen for making the array of hollow microneedles, silicon and glass were respectively selected as the materials for the main body and diaphragm of the micropump. The flappers can be made of different materials such as glass, silicon and PDMS. Since among the mentioned materials, PDMS has a greater flexibility and sealing capability, it was chosen for the flappers. Furthermore, PDMS is a highly elastic polymer and therefore highly durable. Table 3 shows the properties of the materials used in different parts of the micropump.

Table 3

The properties of the materials used in different parts of the micropump

The material properties of the part	piezoelectric disk	ρ(kgm^3)=7500, S_Ematrix(1Pa):S_E−11=1.65×10^−11,S_E−12=−4.78×10^−12,S_E−13=−8.45×10^−12,S_E−14=0,S_E−15=0,S_E−16=0,S_E−22=1.65×10^−11,S_E−23=−8.45×10^−12,S_E−24=0,S_E−25=0,S_E−26=0 S_E−26=0, S_E−33=2.07×10^−11,S_E−34=0, S_E−35=0, S_E−36=0, S_E−44=4.35×10^−11,S_E−45=0, S_E−46=0, S_E−55=4.35×10^−11, S_E−56=0, S_E−66=4.26×10^−11 & S_E−ij=S_E−ji
		dmatrix(CN):d_13=−2.74×10^−10,d_33=5.93×10^−10, d_15=7.41×10^−10
	diaphragm	ρ(kgm^3)=2203,E(MPa)=73100, υ=0.17
	flappers	ρ(kgm^3)=970,E(MPa)=0.75, υ=0.49

Fig. 11 shows the geometry used for the simulation. Regarding the full symmetry, it was enough to model half of the micropump. The geometry consists of the fluid domain and the solid parts. These solid parts are the piezoelectric disk, the diaphragm, and the inlet and outlet flappers whose magnified views are also shown in Fig. 11. Table 4 gives the geometric specifications of the solid parts shown in Fig. 11. For a specified piezoelectric disk diameter, its thickness decrease leads to an increase in compression ratio [35]. Considering this fact along with the system general size, commercially available diameter and thickness were chosen for the piezoelectric disc [48].

Fig. 11

The geometry used for the simulation of the micropump.

Table 4

the geometric specifications of the solid parts.

geometric specifications of	piezoelectric disc	diameter: 6.5 mm, thickness: 0.25 mm
	diaphragm	diameter: 8.5 mm, thickness: 0.1 mm
	flappers	height: 50 μm, thickness: 5 μm

In Fig. 11, the fluid domain is made of a main chamber which has inlet and outlet microchannels connected to it. Since the flappers are placed in the inlet and outlet microchannels, the channels should be stepwise. In those sections of the microchannels where the flappers are placed, the height is 60 μm, whereas other sections which are connected to those sections are 20 μm shorter.

The piezoelectric actuator is polarized in the out-of-plane direction. As the piezoelectric actuator deflects frequently and therefore moves the diaphragm to and fro, correspondingly the input and output flappers are opened and closed opposite to each other. The boundary of the solid and fluid domains is the interface through which these interact with each other. In the simulation there were considered a Lagrangian framework for the solid domain and an Eulerian framework for the fluid one. The fluid domain was constituted by a freely moving deformed mesh, whose movement and deformation was obtained by solving a hyperelastic smoothing PDE (Yoeh model). The finite element method (FEM) was used to perform the numerical simulation. For the temporal discretization, the implicit second order backward differential formula (BDF) was employed. All fully coupled, the differential equations were solved. In each time step, the Newton method was used for linearization of the equations, which were then solved by means of MUltifrontal Massively Parallel sparse direct Solver (MUMPS). First order and second order elements were used for the discretization of the pressure and velocity, respectively. The quadratic elements were utilized to discretize the electric potential and quadratic Lagrange elements were used for discretizing the displacement field.

For the solid parts, the second Piola-Kirchhoff stress tensor (S) and the Green-Lagrange strain tensor (∈) are respectively obtained from Eqs. (6) and (7).where ,u, t, and F are respectively the material density, displacement vector, time, volume force vector (considered 0 here) and deformation gradient. F is obtained from the formula, , in which I is the identity matrix.

Assuming linear elasticity, the Piola-Kirchhoff stress tensor is proportional to the Green-Lagrange strain tensor according to Eq. (8).where C, the elasticity tensor, is a function of Young's modulus (E) and Poisson's ratio (υ).

Eq. (9) which relates to the piezoelectric material, is obtained from Eq. (8) via involving the electric field effect in it and rewriting it in the strain-charge form.where , and are the compliance matrix, second Piola-Kirchhoff stress tensor (from Eq. (6)), coupling matrix and electric field, respectively.

The electric displacement is obtained from Eq. (10) .in which and are the vacuum and relative permittivity, respectively.

Eq. (11) shows the charge conservation, while Eq. (12) yields the electric field.in which , V and are the electrical charge concentration, electrical potential and electrical displacement field, respectively.

Eqs. (1) and (2) are the fluid flow governing equations.

The boundary conditions (BCs)

First, it should be explained here that in the simulation process, while the inlet/outlet flapper valve is closed, there is considered a minimum gap between the flapper and the wall of the stepwise part of the microchannel whereas there is no gap in practical reality. The problem lies with the point that in the simulation, the change in meshing topology does not allow the solid to solid contact, when it comes to fluid-solid interaction (FSI).

For the fluid simulation, at stationary walls, the boundary conditions are considered no-slip ( There should be a 90 Pa fluid pressure at the micropump outlet, if the fluid is to flow in the passages shown in Fig. 4A at a 1 μL/min flow rate before entering the waste fluid bag of a zero gauge pressure (for higher and lower flow rates this pressure changes proportionally). Hence the outlet BC is gained by relation (13).in this relation, "if" means that the nonzero pressure outlet BC is exerted only when the fluid is flowing out and otherwise it is zero.

In practice, there should be set a zero gauge pressure for inlet BC. However, as mentioned before, in the simulation it is considered a minimum gap between the inlet/outlet flapper and the micropump wall (while the flapper is closed), which means when the direction of the fluid movement is towards the outlet, there should be found a way to prevent a significant leakage from the inlet due to the nonzero pressure outlet BC. The solution taken here is setting the inlet BC relation so that the inlet pressure will be zero when the fluid flows inwards and when the fluid flows outwards it will have the pressure as the outlet BC.

On the fluid-solid boundaries, the solid velocity is given to the fluid, and in return the solid boundaries are loaded by the fluid flow stress.

An electrical potential wave signal of specific frequency and amplitude, and a constant zero potential were applied to the piezoelectric upper and lower surfaces, respectively. Zero-charge-density BC was applied to the entire lateral surface.

Results and discussion (for micropump)

Fig. 12 shows the meshing system. As can be seen, the meshing is implemented to be smaller around the flappers compared to other parts. The entire geometry consists of 81,022 tetrahedral elements.

Fig. 12

The meshed micropump geometry for the simulation.

In the proposed micropump, with a 2 Hz frequency, giving a 7.5 V input voltage yields an average flow rate of 0.99 μL/min, which is almost identical to the intended value (1 μL/min). Here, to get this voltage conveniently, the commercial batteries, whose voltages are either 1.5 or 3 V, can be used in series. Fig. 13A shows the input voltage variation versus time, while Fig. 13B and C show the corresponding center of diaphragm displacement and total outflow rate, respectively (The graphs have been drawn for the duration of 1 second as they are periodic functions with the period of 0.5 s) .

Fig. 13

A: The input voltage variation versus time. B: The center of diaphragm displacement versus time. C: The outflow rate versus time. D: The entire pumped fluid value in the course of time.

As can be seen in Fig. 13B, the center of the diaphragm is displaced proportional to the input voltage. The diaphragm moves upward when the curve slope is positive and it moves downward when the slope is negative.

Since the initial velocity and gauge pressure are considered zero, at the time zero, the outflow rate is also zero at this time, as can be seen in Fig. 13C. The time passing and the diaphragm starting to move upward, the fluid starts to be sucked into the chamber. Since the outlet flapper valve is not totally closed at this time, the fluid leaks from the outlet to the inlet, which is the reason for a small negative flow rate noticed. Once the diaphragm reaches its maximum positive deflection, it moves downward, causing a positive outflow rate which starts to rise. Then, when the diaphragm reaches the zero position (zero displacement), its velocity becomes maximum, which is when the outflow rate becomes maximum too. With a negative displacement of the diaphragm, while the outflow rate is still positive, its value takes a decreasing trend. When the diaphragm reaches its maximally negative displacement, the outflow rate gets zero, when again with the diaphragm moving upward, there happens the suction of the fluid from the inlet, while at the same time the slight outlet-to-the-inlet leakage gives a minor negative value to the outflow rate.

By integrating the outflow rate versus time, the entire pumped fluid value in the course of time is obtained (Fig. 13D). As can be seen, each time the fluid volume gets to its local maximum, it slightly decreases in the passage of time, due to the outlet-to-the-inlet leakage.

Here, for example, if the frequency increases to 3 and 5 Hz, the respectively required voltages will be 5 V (results in 0.98 μL/min) and 3 V (results in 0.96 μL/min). Of course, while the frequency increase brings about a decrease in the required voltage, which is desirable, it will increase the probability of the decreased lifetime due to the mechanical fatigue.

In order to verify the gained outflow rate from the micropump, its per second value was divided by the frequency and then compared with the stroke volume. Fig. 14A shows the variations of the main chamber volume versus time.

Fig. 14

A: The variations of the main chamber volume versus time. B: The micropump outflow rate versus time for different numbers of elements.

According to Fig. 14A, the stroke volume ( equals 8.98 nL. Dividing the total outflow per second rate by the 2 Hz frequency gives 8.25 nL, which is satisfactorily close to the stroke volume, where the slight difference (0.73 nL) is due to the fluid leakage in the numerical simulation. In each cycle, there is a 0.336 nL leakage from the inlet and 0.394 nL leakage from the outlet. The total sum of these two leakage values is 0.73 nL which is identical to the mentioned difference.

To show the independence of the micropump results from the number of computational elements, the outflow rate of the micropump versus time was depicted for different numbers of elements (Fig. 14B).

As can be seen, there is a highly acceptable level of overlap of the results for 81,022 elements and more. So, the number of elements chosen for the simulation was 81,022.

Conclusion

In this paper, the design and numerical simulation of the mass-transfer compartment and the piezoelectric micropump of our ultimately intended implantable microfluidic device are presented for the microdialysis-based regular glucose measurement in diabetic patients. The mass-transfer compartment has the appropriate recovery factor of 72%, utilizing 24 microchannels, each of which located on 12 hollow microneedles. The height and width of each microchannel are 60 μm and 50 μm, respectively. The penetration depth of the microneedles is 150 μm, so that they would cause no pain to the patient. The thickness of the semipermeable membranes located in each microchannel at the base of the microneedles is 15 nm. The designed piezoelectric micropump uses two passive flapper valves for letting the fluid flow in and out. With a 2 Hz frequency, and using an input voltage of 7.5 V, the micropump can provide the flow rate of 1 μL/min. The materials and the dimensions of the different parts of the micropump were chosen so as to achieve feasible fabrication of the micropump at appropriate dimensions. In summary, the designed mass-transfer compartment and piezoelectric micropump perfectly fit our ultimately intended device, which enjoys the benefits of the microdialysis method and can be conveniently used by diabetic patients painlessly.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.