Passive flow control valve for protein delivery

Abstract

Passive flow control valves are usually intended to deliver or drain a fluid at a constant rate independently of pressure variations. Microfluidic devices made of a stack of two plates are considered here: the first plate comprises a flexible silicon membrane having through holes while the second plate is a rigid substrate with a cavity, an outlet hole and pillars aligned with the through holes of the membrane. The liquid flows through the holes etched in the membrane and through the gap between the membrane and the pillars. Each gap can be considered as a valve which progressively closes as the pressure increases. Numerical modelling of the fluid dynamics inside the device associated with FEM simulations of the membrane distortion have been performed to design a device that exhibits a constant flow rate in a specified range of pressure. To make the design more reliable, the device characteristics have been optimized using a genetic algorithm, the fitness function taking notably into account machining and alignment tolerances. This algorithm has been finally used to design flow control valves for wearable injectors dedicated to the infusion of viscous drug formulations (hyaluronic acid, adalimumab, golimumab …) at high pressure. Prototypes have been characterized using solutions of 12 and 24 cP. It has been demonstrated experimentally that this technology is suitable to passively infuse biological products at flow rates up to 1 mL/min. The numerical model has then been refined further so as to obtain a good correlation with experimental data.

Keywords:

• microfluidic device
• genetic algorithm
• power-law fluid
• flow control valve
• protein delivery

Public Interest Statement

Passive (battery less) flow control valves are usually intended to deliver or drain a fluid at a constant rate independently of pressure variations. A typical application is the intravenous or intradermic protein infusion for immunotherapy. Using syringes prefilled with such biologics can cause injectability related issues. On-body injectors are therefore the preferred solution to infuse either high volume or viscous drug formulations. A flow control valve in series with the pressurized reservoir limits the risk of leakage, the pain associated with drug overflow and allows a better control of the injection duration. Microfluidic devices made of silicon and glass are considered here. The device characteristics have been optimized using a genetic algorithm. It has been demonstrated experimentally that this technology is suitable to passively infuse viscous (24 cP) biological products at constant flow rates up to 1 mL/min considering reservoir pressure variations of 10 bar or more.

1. Introduction

A passive flow control valve or flow regulator for biomedical applications delivers a constant amount of drug over time independently of pressure variations. The electrical equivalent of such microfluidic device is a current source. Typical examples of flow control devices are proposed by Chappel, Dumont-Fillon, and Mefti (2014) and Zhang et al. (2015).

A microfluidic device made of a stack of Silicon-On-Insulator (SOI) and glass wafers is considered here (Figure 1). The SOI wafer (Top wafer) has been chosen for the etch-stop potential of the buried oxide. The handle part has a typical thickness of about 400 μm in order to provide stiffness to the valve whereas the device layer is usually thin (typically 30 μm) to serve as a membrane. The etching of the holes in the membrane (SOI front side) is obtained by DRIE (Deep Reactive Ion Etching) to get vertical sidewalls. The silicon membrane is then opened and released from the backside whereas the remaining SiO₂ is removed by wet etching.

Figure 1. Schematic cross-section of the valve (not to scale).

The substrate wafer (Bottom wafer) is made of Pyrex or Borosilicate. Silicon can also be used as an alternative material but surface quality constraints are much stricter if fusion bonding is targeted to assemble the stack. The cavity defining a gap between the membrane and the pillars after assembly is made either by dry or wet etching. Controlling the depth and uniformity of this gap is critical for the fluidic behavior of the device. The next step consists in an additional etching at the bottom of the cavity in order to create the pillars as well as the fluidic pathway around them to the outlet. Finally, the outlet is created by sandblasting from the backside of the wafer. The final process step consists in an anodic bonding of the two machined wafers.

A typical low pressure application is the derivation of cerebrospinal fluid from the brain ventricles toward the peritoneum. Conversely high pressure applications include, for instance, intravenous protein infusion for immunotherapy. Using prefilled syringes can cause problems with the injectability or syringeability of such biologics. On-body injectors are therefore the preferred solution to infuse either high volume or viscous drug formulations. Depending on the technology used to pressurize the fluid, it could be relevant to add a flow control valve to the device. This additional valve actually limits the risk of leakage induced by overpressure, the pain related to drug overflow and allows a better control of the injection duration.

The dimensioning of the valve to get the expected flow regulation profile is made step by step (Chappel et al., 2014). This iterative process starts with the positioning of a first valve located on the outer edge of the membrane. This valve is intended to determine the high pressure fluidic characteristics of the device. Additional valves located near the center of the membrane are then introduced to gradually define the expected flow rate profile at lower pressures. There are potentially a large number of valve configurations that meets the requirement for the targeted flow rate profile. Therefore, in order to find an optimum configuration, there is a need to establish design rules to sort them. Instead of mimicking the manual iterative design phase described above, we have considered a genetic algorithm that is especially adapted to solve complex problems with a large number of variables.

The key parameters to achieve delivery accuracy within the range ±10% are respectively: the membrane thickness, the gap depth, the hole diameters and the alignment between the holes and the pillars. The two first parameters can be easily controlled by using SOI wafers. The impact of the other machining or alignment tolerances on delivery accuracy shall be reduced to the maximum by design. Instead of trying to reduce further the tolerances, we chose to improve the valve design in order to make it more reliable. But, given the complexity of this optimization problem, a specific program has been developed in order to take into account machining tolerances and flow regulation specifications.

This paper aims to demonstrate that flow control valves based on MEMS technology is suitable for the infusion of viscous drug formulations at high pressure. In addition, modelling tools are provided to adapt the design to the injection requirements.

After a presentation of the numerical modelling, the genetic algorithm is briefly described. Using this original design tool, passive flow control valves for delivery of viscous solutions of 12 and 24 cP have been conceived and fabricated in clean room. After a complete metrology of the samples, a fluidic characterization has been carried out and the experimental results are compared to simulations to refine the numerical modelling.

2. Valve modelling

Microvalves are the key elements of MEMS micropumps and other microfluidic systems. Reviews of the different valve technologies are provided by Oh and Ahn (2006) and Au, Lai, Utela, and Folch (2011). A comprehensive mathematical description of fluidics in microvalves is required when the dynamic of the valve opening has an impact on the device behavior. Several analytical formulae related to laminar fluid flow through annular valve are proposed here since such configurations are not fully addressed in classical textbooks of fluid mechanic (see for instance Bird, Stewart, & Lightfoot, 2007; Guyon, Hulin, & Petit, 2001; White, 1991 for detailed descriptions of many solutions to the Newtonian viscous-flow equations). These formulae have been implemented in simulation models for passive flow control valves. Non-idealities like valve tilt and misalignment have been also considered in order to build a model accounting for realistic manufacturing conditions. Finally, to account for the specific geometry of the valves, notably for the transition between the hole and the diffuser, corrective factors have been introduced into the original model (Chappel et al., 2014). The singular losses now include a better description of the bend of the flow between the hole and the diffuser as well as the radial flow expansion inside the diffuser itself. Additionally the flow constriction between the hole and the diffuser has been better simulated simply by using a corrective multiplication factor equal to ε = 0.6 for the inner radius of the diffuser itself. These corrections are based on the experimental analysis of specific valves containing only a single centered hole. This investigation has been conducted using Silicon-On-Insulator wafers having device layers of 30 and 129.5 μm respectively to form membranes (from 3.76 to 5.12 mm in diameter) and Pyrex wafers for the substrate and pillars. Both wafers have then been anodically bonded together after alignment. The target gap was fixed to 10.46 μm. The ratio between pillar and hole radii varies from 1.25 to 2.8. The fluidic chips have been tested using water at Reynold numbers comprised in the range 0.5–35. A detailed description of this experimental plan is provided elsewhere (Musard, Dumont-Fillon, van Lintel, & Chappel, 2018).

2.1. Axisymmetric case

The annular valve geometry considered here consists of a flexible membrane having a hole of radius R_in and a pillar of radius R_out as depicted in Figure 2. Fluid flows through the hole and the small gap between the pillar and the membrane (diffuser). The opening of this diffuser is usually a function of the pressure gradient through the valve. We consider a power law or Ostwald de Waele fluid relationship between the shear stress σ and the shear rate $\dot{\mathit{\gamma }}$ where the apparent dynamic viscosity η can be approximated by:

Figure 2. Cross-section of axisymmetric annular diffuser geometry (the origin of the axis is located along the pillar axis).

(1) $$\mathit{\eta }=\raisebox{1ex}{$\mathit{\sigma }$}\!\left/ \!\raisebox{-1ex}{$\dot{\mathit{\gamma }}$}\right.=m{\left(\dot{\mathit{\gamma }}\right)}^{n-1}$$(1)

where n and m in Equation (1) are two empirical curve-fitting parameters, n = 1 corresponding to Newtonian fluids. The incompressible fluid flow is assumed to be fully developed at the diffuser entrance and the pressures P_in and P_out outside the diffuser are considered as homogenous.

According to the diffuser symmetry, the small value of the gap h and the nature of the fluid, four assumptions are made regarding this flow:

(1)	Steadiness (∂/∂t = 0).
(2)	Axisymmetry (∂/∂θ = 0 and $u_{\mathit{\theta }}=0$).
(3)	The vertical velocity is neglected (uz = 0).
(4)	Small laminar $R{e}_{\mathit{genPL}}$ for Power-Law fluid (Madlener, Frey, & Ciezki, 2009), and therefore the convective terms can be neglected:

(2) $$R{e}_{genPL}=\frac{\mathit{\rho }{\left(h\right)}^n{\left(\frac{\overline{Q}}{2\mathit{\pi }h{R}_{\mathit{in}}}\right)}^{2-\text{n}}}{m{\left(\frac{3n+1}{4n}\right)}^n{8}^{n-1}}$$(2)

The Navier-Stokes equations coupled to the boundary conditions:

(3) $$u_r\left(r,z=\pm h/2\right)=0$$(3)

Yield the fluidic resistance $R_f$ as a function of the flow rate $Q$:

(4) $$R_f={(P_{\mathit{in}}-P_{\mathit{out}})}^{1/n}/Q$$(4)

Following the value of n, two different expressions are derived as shown in Equations (5) and (6) depicting the function $R_f(h,R_{\mathit{in}},R_{\mathit{out}},n,m)$:(5) $$R_f=\frac{6m}{\mathit{\pi }{h}^3}ln\left(\frac{R_{\mathit{out}}}{R_{\mathit{in}}}\right),n=1$$(5)

(6) $$R_f=\frac{\left(2n+1\right)m^{\frac{1}{n}}}{4\mathit{\pi }n{\left(\frac{h}{2}\right)}^{2+\frac{1}{n}}}{\left(\frac{R_{\mathit{out}}^{1-n}-R_{\mathit{in}}^{1-n}}{1-n}\right)}^{\frac{1}{n}},n\ne 1$$(6)

Detailed derivations are provided in the Appendix A. FEM simulations corroborate these calculations with a relative error of less than 0.3% for $0.5\le n\le 2$.

2.2. Non-idealities

Two non-idealities have been introduced into the model:

•	the misalignment between the hole of the membrane and its facing pillar;
•	the tilt of the membrane with respect to the pillar.

The surface of the pillar is considered as the reference plane z = 0. The membrane center is noted O. Two off-axis parameters δ and $\mathrm{\Omega }$ are also introduced to account for machining tolerances and wafer misalignment if the valve is a stack of two wafers:^{Figure 3}

•	δ is equal to the distance between the pillar center Op and the projection $O_h$ of the hole center on the plane z = 0.
•	$\mathrm{\Omega }$ is equal to the angle between OOh and OhOp (see Figure 3).

Figure 3. Geometry (right) and new coordinates (left) for a misaligned diffuser (top view), the angles being conventionally measured counter clockwise.

The dimensionless parameter used to evaluate the misalignment error is $\mathit{\epsilon }=\frac{\mathit{\delta }}{R_{\mathit{out}}}$. The term $\mathrm{\Omega }$ indicates the direction of the error: for $-{90}^{\circ }\le \mathrm{\Omega }\le {90}^{\circ }$ the hole is closer to the centre of the membrane, for ${90}^{\circ }\le \mathrm{\Omega }\le {270}^{\circ }$ the hole gets closer to the periphery of the membrane.

Note that if $\mathit{\psi }=0$, $\mathrm{\Omega }$ is useless as the direction of misalignment does not influence the local height.

Tilt shall be considered when the pillar is not centered with respect to the flexible membrane or in case of particle contamination preventing a correct closure of the valve. This induces a non-homogeneous distance between the membrane and the pillar h(r, θ) as seen in Figure 4.

Figure 4. Schematic cross-section of a tilted radial diffuser.

As the deformation of the membrane is likely to be planar above the pillar (membrane diameter >> pillar diameter), a constant tilt angle ψ will be considered in the following study, which results, in considering both kinds of non-idealities, in a distance membrane-pillar of:

(7) $$h\left(r,\mathit{\theta }\right)=h_0+rcos\left(\mathit{\theta }\right)tan\left(\mathit{\psi }\right)-\mathit{\delta }cos\left(\mathrm{\Omega }\right)tan\left(\mathit{\psi }\right)$$(7)

where h₀ is the distance between the membrane without hole and the centre of the pillar. An interesting non-dimensional parameter accounting for the influence of this tilt is $\mathit{\xi }=R_{\mathit{out}}tan\left(\mathit{\psi }\right)/h_0$. If ξ ≥ 1 there is contact at some point on the pillar, this will be studied later on.

Another system of coordinates is also introduced for further description of the problem in order to keep the origin in the center of the hole. Thus, the link between the old $\overrightarrow{r^{\prime }}$and the new $\overrightarrow{r}$ can be written as:

(8) $$\left\{\begin{array}{c}\overrightarrow{r^{\prime }}=\left(\begin{array}{c}rcos(\mathit{\theta })-\mathit{\delta }cos(\mathrm{\Omega })\\ rsin(\mathit{\theta })-\mathit{\delta }sin(\mathrm{\Omega })\end{array}\right)\hfill \\ r^{\prime }=\Vert \overrightarrow{r^{\prime }}\Vert \hfill \end{array}\right.$$(8)

The total fluidic resistance is estimated through a parallelization of infinitesimal fluidic resistances around the diffuser as depicted in Figure 5.

Figure 5. Electrical analogy.

The discrete and continuous versions of such a calculus would be

(9) $$R_{f_{\mathit{tot}}}={\left({\sum }^{}\frac{1}{R_{f_i}}\right)}^{-1}$$(9)

And

(10) $$R_{f_{\mathit{tot}}}={\left(\underset{-\mathit{\pi }}{\overset{\mathit{\pi }}{\int }}\frac{d\mathit{\theta }}{R_{f,axi}\left(\mathit{\theta }\right)}\right)}^{-1}$$(10)

where $R_{f,axi}\left(\mathit{\theta }\right)$ represents the infinitesimal fluidic resistance between the angles θ and $\mathit{\theta }+d\mathit{\theta }$. It can be computed as the resistance per angle for a diffuser with an axisymmetric height. This latter value is obtained by summing numerically the resistance:

(11) $$R_{f,axi}\left(\mathit{\theta }\right)=\frac{1}{2\mathit{\pi }}{\left(\underset{R_{\mathit{in}}}{\overset{R_{\mathit{max}}\left(\mathit{\theta }\right)}{\int }}\frac{\mathit{\partial }{R}_f{\left(h\left(r^{{}^{\prime }},{\mathit{\theta }}^{{}^{\prime }}\right),r=R_{\mathit{out}}\right)}^n}{\mathit{\partial }r}dr\right)}^{1/n}$$(11)

where r′ and θ′ are defined in function of r and θ in Equation (8) and $\mathit{\partial }{R}_f^n\left(h\left(r^{{}^{\prime }},{\mathit{\theta }}^{{}^{\prime }}\right),r=R_{\mathit{out}}\right)/\mathit{\partial }r$ is the partial derivative of Equations (5) and (6) evaluated in h(r′, θ′) and $r=R_{\mathit{out}}$ (see the Appendix A for the explanation related to the different exponents used in Equation (11), for an axisymmetric variable height). Finally, R_max(θ) is the exit radius “in front of” the angle θ on Figure 3 and can be evaluated as:

(12) $$R_{\mathit{max}}\left(\mathit{\theta }\right)=R_{\mathit{out}}\left(\mathit{ϵ}cos\left(\mathit{\theta }-\mathrm{\Omega }\right)+\sqrt{1-{\mathit{ϵ}}^{2}si{n}^2\left(\mathit{\theta }-\mathrm{\Omega }\right)}\right)$$(12)

Note also that Equation (11) can be calculated analytically in the case of Newtonian fluid (n = 1).

Finally, in case of contact between the pillar and the membrane (ξ ≥ 1), the integration range in Equation (10) is reduced. For small misalignments ($\mathit{ϵ}\le 0.2$) and closing angles (ξ ≤ 1.2), correlations between the numerical modelling and FEM simulations exhibit deviations better than 15% with the estimated fluidic resistances, in the same range 0.5 < n < 2.

3. Genetic algorithm

Conventional optimization algorithms, like gradient-descent, become very complex for multivariable and non-linear systems. Some bio-inspired algorithms are being increasingly used to solve complex problems (Floreano & Mattiussi, 2008). Based on Darwinian principles, a genetic (or evolutionary) algorithm relies on random mutation, selection and reproduction of the best individuals (being valves in our case) over generations.

In the algorithm, each individual (valve) is described by a genotype (set of parameters): number of holes and pillars, hole position and diameter, pillar position and diameter. These parameters vary from one valve to the other and may change during mutations. Other characteristics like the global geometry of the device or the membrane thickness are kept constant. These parameters are initially defined by simple consideration on the regulation pressure range, burst pressure and fluidic interconnection issues.

A population of N = 20 valves is initially generated randomly. Two models have been developed: a finite element model to compute the membrane deflection as a function of pressure, and an analytical fluidic model, to derive the flow rate from a given set of deflection and pressure. By combining these two models, it is possible to determine the flow going through a valve for a given pressure applied on the membrane. The fitness function is a measure of the valve performance, like a grade. It is used to determine which are the best individuals to be selected and reproduced for the next generation (Dumont-Fillon, Hannebelle, van Lintel, & Chappel, 2016).

The fitness function is set to measure two main parameters:

•	Proper flow control: its flow curve must be as close as possible to a target value in a given range of pressure.
•	Sensitivity to micro-machining tolerances.

The sensitivity to machining tolerances is obtained by performing a fluidic simulation of the design not with the ideal but altered dimensions, taking into account the misalignment and worst combination of tolerances at 3σ for all the critical dimensions that impact the flow rate, including the membrane thickness and radius, the hole and pillar diameters as well as the gap value. A typical example of fluidic characteristics after optimization is provided in Figure 6, considering a maximum misalignment of 5 μm and tolerance of ±0.5 μm at 3σ for the critical dimensions listed here above except for the gap ($\pm 0.3\mathit{\mu }m$). It is shown that an error on the flow rate of about ±8% is expected at 3σ.

Figure 6. Simulation of the impact of process tolerances at 3 sigma on the flow profile of a valve designed to deliver a constant flow rate of 0.5 mL/min at 24 cP in the range 7–21 bar.

Since the validity of this model determines the validity of the design generated by the algorithm, it has been validated on fabricated valves (Chappel et al., 2014). Another validation for valves exhibiting very small gaps even at low pressure has been conducted recently to tune the fluidic model. It has been inferred that FEM simulations tend to overestimate the membrane stiffness after contact with the pillars. The estimation of the gap between the hole and the pillar as a function of the pressure includes therefore some corrections based on experimental data too. Microfluidic diffusers as described in (Musard et al., 2018) have been machined in a similar way, using exactly the same dimensions except for the radial position of the diffuser along the membrane radius. By using the fluidic model described in (Musard et al., 2018) and comparing the fluidic characteristics of out-of-center and centered diffusers, it has been possible to refine the estimation of the gap for different radial positions as a function of the applied pressure, using a fluidic characterization of the diffusers.

4. Design description; microfabrication and metrology

Three different designs have been made for infusing newtonian viscous fluids, in the pressure range 5–20 bar, according to the target flow rates provided in Table 1. For the two solutions at 12 and 24 cP, obtained by mixing glycerol and water (see Table 2 for the glycerol fraction mass), densities of 1.145 and 1.170 have been respectively introduced in the model according to (Segur & Oberstar, 1951).

Table 1. Target flow rates at 20°C for the designs A, B and C

Design	Target flow rate at 20°C
A	0.5 mL/min at 12 cP
B	0.5 mL/min at 24 cP
C	0.135 mL/min at 24 cP

Table 2. Physical characteristics of the water/glycerol mixtures of 12 and 24 cP according to Segur and Oberstar (1951)

Fluid (cP)	Mass fraction of glycerol (%)	Density
12	61.407	1.145
24	70.453	1.170

The silicon membranes have been made using SOI wafers having a device layer 130 μm thick (129.5 μm after oxidation) in order to sustain pressures up to 35 bar. The gap has been fixed at 10.46 μm. All membrane diameters are equal to 5.12 mm. The height of the pillars is 50 μm and the outlet hole has a diameter of 500 μm. Overall chip dimensions are 10 × 10 mm.

The process flow is similar to the one described in (Chappel et al., 2014). Microfabrication has been carried out at CMi (EPFL, Switzerland). A photo of a 4′′ wafer after dicing is provided in Figure 7.

Figure 7. Photo of the 4′′ wafer on its dicing tape.

The hole and pillar dimensions have been carefully characterized using the microscope Nikon Eclipe LV150 associated to the software Kappa Metreo. The targeted and measured hole and pillar dimensions are provided for each design in Tables 3 and 4 respectively. The values indicated by the software have been reported here with all digits even if they are not really significant because of the measurement repeatability estimated at 0.5 μm. The goal of the Tables 3 and 4 is just to show that both holes and pillars are within the specifications at ±1 and ±5 μm respectively.

Table 3. Comparison between theoretical and measured hole diameters for each design

Design	Inlet side (μm)	Outlet side (μm)	Theory (μm)	ΔØ (μm)
A	99.30	97.06	96.70	0.36
	122.65	120.71	120.32	0.39
	118.31	116.31	115.79	0.52
	122.35	120.43	119.81	0.62
B	120.62	118.89	118.26	0.63
	135.81	134.03	133.39	0.64
	151.80	149.94	149.37	0.57
	137.50	135.06	134.95	0.11
C	43.71	40.94	40.08	0.86
	94.24	92.27	91.47	0.80
	87.65	85.67	84.80	0.87
			Max (μm)	0.64
			Target (μm)	± 1

Table 4. Comparison between theoretical and measured pillar diameters for each design

Design	Outlet side (μm)	Theory (μm)	ΔØ (μm)
A	303.21	305.40	−2.19
	201.77	203.00	−1.24
	276.73	278.76	−2.03
	210.31	212.40	−2.10
B	230.04	227.00	3.04
	217.48	214.00	3.48
	202.52	199.35	3.16
	214.14	211.00	3.14
C	301.71	301.20	0.51
	253.12	253.00	0.12
	272.14	273.00	−0.86
		Max (μm)	3.48
		Target (μm)	± 5

5. Experimental setup and test methods

5.1. Test setup

The layout of the test setup, as provided in Figure 8, comprises a 20L pressurized bottle (4.5 nitrogen at 200 bar) associated with a first pressure regulator (Gloor) that maintains at 30 bar the inlet of the pressure controller GE Pace 5000 (regulation at ±1 mbar in the range −1 to 21 bar). The outlet pressure, which is monitored by a second sensor GE Unik 5000 (range −1 to +30 bar, resolution 1 mbar), is applied directly onto the 5 liter stainless steel (SS) tank (Festo CRVZS-5) partially filled with the glycerol/water mixture. The viscous fluid flows through two filters toward the sample holder, the flow meter Sensirion SLQ-QT500 (accuracy ± 5%) and finally a beaker placed onto a scale Sartorius MC1 620 (range 0–620 g, resolution 0.001 g). Pneumatic couplings are made by Festo while high pressure fluidic tubing, frits and fittings are products of Upchurch Scientific. The filtering system, which consists of a 2 μm prefilter in PEEK (frit OC-802) associated with a 2 μm SS filter (C-407), is necessary to prevent the contamination of the chips since the complicated fluidic line and the big SS tank being very difficult to purge properly. A computer is finally connected to the pressure controller, the flow sensor, the pressure sensor and the scale for data acquisition.

Figure 8. Layout of the test setup comprising the pressurization system, the viscous fluid tank, the sample and measurements means (scale, pressure sensor and scale).

5.2. Sample holder

A specific sample holder has been machined for these high pressure tests. A picture is provided in Figure 9. Fittings from Upchuch Scientific are used to connect the sample holder to the fluidic line. The tightness is secured by two orings in MQV, shore 50. Overall sample holder dimensions are 40 mm x 40 mm x 24 mm. Since the pressure is applied onto the top of the chip (silicon membrane side), the sample holder has a flat reference surface with only one outlet hole on the glass wafer side to prevent any bending or breakage of the chip. The fluidic resistance of the sample holder itself is more than 3 orders of magnitude smaller that the chip itself.

Figure 9. Photo of the two parts that compose the sample holder with its oring grooves, sample cavity and wing nuts screws.

5.3. Test preparation

Fluids of 12 and 24 cP are obtained by mixing pure deionized water filtered at 0.2 μm and glycerol (99.5% from Reactolab, density 1.26) at 20°C. Mixture preparation procedure: an initial mass of glycerol is measured using the scale Sartorius MC1 620, and the mass of water needed to get the right viscosity is calculated and mixed with the glycerol using a magnetic stirrer. Typical mixtures used during the tests are described in Table 5, including an estimation of the final density.

Table 5. Description of the measured mass fractions for 12 and 24 cP water/glycerol mixtures

Viscosity Target (cP)	Glycerol mass (g)	Water mass (g)	Total mass (g)	Total volume (mL)	Density (g/mL)
12	613.68	385.69	999.37	872.738	1.1451
24	991.374	415.77	1,407.144	1,202.575	1.1701

The flowmeter, which is originally calibrated for water and ethanol, has been calibrated using the scale for 12 and 24 cP respectively. A corrective factor has been derived to obtain the exact flowrate in μL/s. A systematic check of this corrective factor is made for each test using the scale and the estimated density of the mixture. This test with the scale is made at the highest flowrate to minimize measurement errors. The viscosity of the mixtures has been estimated experimentally using a reference PEEK capillary tubing (length 126 mm, ID = 0.4 mm). The flow rate through this tubing has been measured at 20°C using iteratively water and the different mixtures. Test pressure has been varied so as to reach the typical regulation flow rate of the devices. The applied pressure versus flow rate characteristic is fitted linearly to determine the fluidic resistance of the tubing for each fluid. The ratio of the fluidic resistance obtained at nominally 12 cP and 24 cP against that obtained at 1 cP indicates the real viscosities:

•	for 12 cP experimental mixture: Rf (12 cP) / Rf (1 cP) = 12.04
•	for 24 cP experimental mixture: Rf (24 cP) / Rf (1 cP) = 23.71

The characterization of the different elements of the fluidic line has been made using a similar protocol using the flowmeter:

•	all fluidic restrictions of the fluidic line are removed except the element to be tested.
•	the reservoir pressure is gradually increased in order to get 6 different flow rates from 0 to 2 mL/min.
•	the pressure versus flowrate characteristic ∆(Q) is fitted linearly to derive the fluidic resistance using the formula ΔP = Rf × Q.

The results obtained using pure water are provided in Table 6.

Table 6. Measured fluidic resistances with pure water of the main elements of the fluidic line

Fluidic line element	Rf (Pa.s/m^3)
Reference PEEK tubing (for viscosity cal.)	2.38E+11
2 μm PEEK Prefilter	4.80E+11
2 μm SS filter	1.36E+12
Others (sample holder, flow meter …)	<<1E+10

5.4. Test method

The tests have been performed in a class ISO 7 cleanroom with temperature regulation at 20 ± 2°C. Three samples per design are tested. The test method proposed here is valid for both mixtures (12 and 24 cP):

•	Pressure sweep from 0 to 20 bar using step of 1 bar at low pressure (<7 bar) and steps of 2 bar at higher pressure.
•	Measurement of the flowrate with the flowmeter Sensirion*.
•	Final visual inspection.

*The fluid is infused into a beaker placed on the Sartorius scale for an additional measurement of the flowrate (at least 2 verifications per test). The weight in g measured during 1 minute is converted in mL/min using the estimated density of the fluid (see Table 5).

6. Results

Experimental flow rate versus applied pressure characteristics for 3 different devices of design A, B and C are provided in Figures 10^{Figure 11}–12 respectively. At low pressure, all curves exhibit a linear regime with a slope that depends on the fluidic resistances of the fluidic line and the diffusers themselves. The flow rate reaches a maximum value before decreasing and showing a sudden change of slope at the contact between the membrane and the substrate at about 4, 5 and 7 bar for the designs A, B and C respectively. At higher pressure, the flow rate is basically constant for every design in the range 7–21 bar, as expected. Based on the metrology data provided in paragraph §4, theoretical curves are also shown for the sake of comparison. It is important to note that the model does not take into account the additional fluidic resistances of the fluidic line, which result in a shift of the whole fluidic characteristics towards high pressure. Setup limitations do not allow testing at higher pressures but flow regulation up to 35 bar is theoretically expected. This pressure corresponds to the maximum pressure that such a membrane can sustain. This limit has been evaluated on three different chips by applying pneumatic pressure directly on the chip (by shunting the pressure controller that is limited to 21 bar). This experimental burst pressure matches theoretical expectation corresponding, in FEM simulations, to a maximum stress of 0.8 GPa on the membrane edges. The flatness of the fluidic characteristics confirms the pertinence of the membrane deformation modelling despite its highly non-linear behavior. The flow rate, for the design A tested at 12 cP, is slightly larger than expected (about 15%) while designs B and C tested at 24 cP are in the pressure range of interest neatly aligned with the flow rate targets, within the range of measurement accuracy (approximately 5%). An important result is the low dispersion of the fluidic characteristics despite the measurement errors: it could be attributed to the robustness of the design because the machining tolerances have been considered in the genetic algorithm to select the best candidates.

Figure 10. Theoretical fluidic characteristic for design A at 12 cP compared to experimental data obtained using 3 different samples F6, H4 and G5.

Figure 11. Theoretical fluidic characteristic for design B at 24 cP compared to experimental data obtained using 3 different samples C7, D7 and A4.

Figure 12. Theoretical fluidic characteristic for design C at 24 cP compared to experimental data obtained using 3 different samples E0, F0 and E1.

These experimental data show that silicon microfluidic valves comprising flexible membranes having through holes are suitable for the passive flow regulation of viscous fluids (24 cP) in the range 7–21 bar. It has been confirmed that the utilization of a genetic algorithm is especially useful in order to conceive sturdier valve designs with respect to machining tolerances. Finally, the good match between simulations and experimental data suggests that the modelling is suitable to design more exotic fluidic profiles.

7. Conclusion

The design of passive flow regulators dedicated to infusion of viscous fluids has been optimized using a genetic algorithm. Machining and alignment tolerances for estimating the fitness function have been considered. Prototypes have been micromachined and characterized experimentally using 12 and 24 cP fluids. It has been demonstrated experimentally that this technology is suitable to passively infuse, with very high accuracy, viscous drug formulations (hyaluronic acid, adalimumab, golimumab …) at flow rates up to 1 mL/min independently from reservoir pressure variations that can reach 20 bar or more. The numerical modelling of the fluid dynamics has been adapted to account for the specific geometry of the valves in order to fit correctly the experimental data. The genetic algorithm arises as a powerful tool for fast and reliable microfluidic chip design optimization.

Nomenclature
Rf	=	Fluidic resistance
ΔP	=	Pressure gradient between the inlet and the outlet of the valve
Q	=	Flow rate
u	=	Fluid velocity
η	=	Dynamic viscosity
$\dot{\mathit{\gamma }}$	=	Shear stress
n, m	=	Empirical curve-fitting parameters
Rin	=	Hole radius
Rout	=	Pillar radius
h	=	Gap between the membrane and the pillar
ρ	=	Volumetric mass density
$\mathit{\delta },\mathrm{\Omega }$	=	Off-axis parameters
ε, ξ	=	Dimensionless parameters
θ	=	Polar coordinate angle
r	=	Radial coordinate
z	=	Vertical axis
ψ	=	Tilt angle
$R{e}_{\mathit{genPL}}$	=	Generalized Reynolds number for power-law fluids

Funding

This work was funded by Debiotech S.A.

Acknowledgement

Authors wish to thank the staff of Center of MicroNanoTechnology at EPFL for support.

Additional information

Notes on contributors

E. Chappel

Eric Chappel was born in Chambéry (France) and works as R&D project manager at Debiotech S.A. He received a master degree in physics in 1996 (Université Grenoble Alpes) and accomplished his PhD thesis at French National High Magnetic Field Laboratory (Grenoble, 2000). His research interests are innovative medical devices, including insulin patch micropumps, hydrocephalus shunts, external ventricular drains, implantable pumps and on-body injectors.

Appendix A

We consider first that the valve opening h is small compared to R_in and R_out. We assume a steady, radial and laminar flow through the valve. According to the symmetry of the system described in Figure 1 and the former hypothesis, we have:

(13) $$u_r\gg u_z\gg u_{\mathit{\theta }}$$(13)

where u is the flow velocity.

We write now the Navier-Stokes equation in the cylindrical coordinate system:

(14) $$\left\{\begin{array}{c}\frac{\mathit{\partial }P}{\mathit{\partial }r}=\mathit{\eta }\frac{{\mathit{\partial }}^2{u}_r}{\mathit{\partial }{z}^2}\\ \\ \frac{\mathit{\partial }P}{\mathit{\partial }z}=0\\ \\ \frac{\mathit{\partial }r{u}_r}{\mathit{\partial }r}=0\end{array}\right.$$(14)

where η is the dynamic viscosity of the fluid and P the pressure in the incompressible fluid.

Ideal valve with Newtonian fluid

We assume that the tilt angle ψ between the membrane and the pillar is equal to zero. The integration of the term of Equation (2) leads to:

(15) $$u_r=\frac{1}{2\mathit{\eta }}\frac{\mathit{\partial }P}{\mathit{\partial }r}z\left(z-h\right)$$(15)

The integration of the 3rd part of Equation (14) gives, using the following boundary conditions: P = P_in for $r<R_{\mathit{in}}$ and P = P_out for $r>R_{\mathit{out}}$:

(16) $$\frac{\mathit{\partial }P}{\mathit{\partial }r}=\frac{1}{r}\left(\frac{P_{\mathit{in}}-P_{\mathit{out}}}{ln\left(\frac{R_{\mathit{in}}}{R_{\mathit{out}}}\right)}\right)$$(16)

Finally we derive the expression of the flow rate Q by an additional integration:

(17) $$Q=\underset{0}{\overset{h}{\int }}2\mathit{\pi }r{u}_{r}dz=\frac{\mathit{\pi }{h}^3}{6\mathit{\eta }}\frac{P_{\mathit{in}}-P_{\mathit{out}}}{ln\left(\frac{R_{\mathit{out}}}{R_{\mathit{in}}}\right)}$$(17)

Valve tilted by a single rigid particle

A tilted valve is typically a normally closed valve which has trapped a single rigid particle of diameter h₀ onto its valve seat. Because it results in a permanent opening of the valve, it is interesting to have an estimation of the leakage induced by such defects. We suppose that the width of the valve seat R_out − R_in is small compared to R_in and R_out.

The Equation (15) then becomes:

(18) $$u_r=\frac{1}{2\mathit{\eta }}\frac{\mathit{\partial }P}{\mathit{\partial }r}z\left(z-h(\mathit{\theta })\right)$$(18)

where

(19) $$h\left(\mathit{\theta }\right)=h_{0}si{n}^2\left(\frac{\mathit{\theta }}{2}\right)$$(19)

Finally we obtain the expression of the flow rate Q:

(20) $$Q=\underset{0}{\overset{2\mathit{\pi }}{\int }}d\mathit{\theta }\underset{0}{\overset{h\left(\mathit{\theta }\right)}{\int }}r{u}_{r}dz=\frac{P_{\mathit{in}}-P_{\mathit{out}}}{2\mathit{\eta }ln\left(\frac{R_{\mathit{out}}}{R_{\mathit{in}}}\right)}\underset{0}{\overset{2\mathit{\pi }}{\int }}\frac{h{\left(\mathit{\theta }\right)}^3}{6}d\mathit{\theta }=\frac{h_0^3}{12\mathit{\eta }}\frac{P_{\mathit{in}}-P_{\mathit{out}}}{ln\left(\frac{R_{\mathit{out}}}{R_{\mathit{in}}}\right)}\underset{0}{\overset{2\mathit{\pi }}{\int }}si{n}^6\left(\frac{\mathit{\theta }}{2}\right)d\mathit{\theta }=\frac{5}{16}\frac{\mathit{\pi }{h}_0^3}{6\mathit{\eta }}\frac{P_{\mathit{in}}-P_{\mathit{out}}}{ln\left(\frac{R_{\mathit{out}}}{R_{\mathit{in}}}\right)}$$(20)

We deduce that the flow rate, for a valve showing a homogeneous opening of h, is 3.2 times higher than the same tilted valve having trapped a single rigid particle of diameter h onto its valve seat.

General case

We consider now the general case according to Figure 2, wherein the membrane exhibits a constant tilt angle ψ but not in contact with the valve seat. This case corresponds typically to the passive valves described in (Chappel et al., 2014), in particular in the high pressure regime where the distortion of the membrane is large. The expression (19(19) $$h\left(\mathit{\theta }\right)=h_{0}si{n}^2\left(\frac{\mathit{\theta }}{2}\right)$$(19) ) becomes:

(21) $$h\left(\mathit{\theta }\right)=h-\mathit{\psi }{R}_{\mathit{out}}cos\mathit{\theta }$$(21)

The expression of the flow rate through the valve is therefore:

(22) $$Q=\frac{\mathit{\pi }{h}^3}{6\mathit{\eta }}\frac{P_{\mathit{in}}-P_{\mathit{out}}}{ln\left(\frac{R_{\mathit{out}}}{R_{\mathit{in}}}\right)}\left(1+\frac{3t^2}{2h^2}\right)$$(22)

Analytical solutions for laminar flow of non-newtonian fluids

It is well known that protein solutions submitted to high shear, for instance during a bolus injection, exhibit a non-Newtonian behavior. This behavior is usually due to the stretching, orientation or deformation of the proteins in the direction of the flow, so the viscosity is a function of the shear rate. We focus here on steady laminar flows and time-independent fluid behavior.

Power law fluid

We consider power law or Ostwald de Waele fluid (Chhabra & Richardson, 2008), wherein the relationship between the shear stress σ and the shear rate $\dot{\mathit{\gamma }}$ can be approximated by:

(23) $$\mathit{\sigma }=m{\left(\dot{\mathit{\gamma }}\right)}^n$$(23)

See also Equation (1) for the expression in term of apparent viscosity.

Shear-thinning behavior fluids are characterized by $0<n<1$ while newtonian fluid corresponds to n = 1. The expressions of the flow rate through circular conduits or between two infinite parallel plates can be found in many fluid rheology textbooks, e.g. (Chhabra & Richardson, 2008). We derive here the expression for the annular valve described previously (with cylindrical symmetry, see (Na & Hansen, 1967) for the derivation of the general formula using the Sisko’s model).

We write the conservation of the momentum of an annular fluid element comprised between r and r + dr at ± z. The top side of the valve seat is located at the coordinate $z=-\frac{h}{2}$ while the bottom side of the membrane is located at $z=+\frac{h}{2}$. We assume that the change of the shear stress along r is small. According to the valve geometry, we have:

(24) $$\left\{\begin{array}{c}{u}_r\gg u_z\gg u_{\mathit{\theta }}\\ \\ \mathit{and}\\ \\ \frac{\mathit{\partial }{u}_r}{\mathit{\partial }z}\gg \frac{\mathit{\partial }{u}_r}{\mathit{\partial }r}\end{array}\right.$$(24)

We get:

(25) $$2\mathit{\pi }r\times 2z\times P-2\mathit{\pi }\left(r+dr\right)\times 2z\times \left(p+dp\right)=2\mathit{\sigma }\left(\mathit{\pi }{\left(r+dr\right)}^2-r^2\right)$$(25)

Or, after simplification:

(26) $$\mathit{\sigma }=-z\frac{\mathit{dp}}{\mathit{dr}}$$(26)

Using the symmetry of the system we only solve the Equation (26) in the mid-plane $0<z<\frac{h}{2}$. Since $\frac{d{u}_r}{\mathit{dz}}$ is negative in this region, the shear stress for a power-law fluid is given by:

(27) $$\mathit{\sigma }={\left(-\frac{d{u}_r}{\mathit{dz}}\right)}^n$$(27)

Combining Equations (26), (27) and integrating with respect to z yields for the velocity distribution:

(28) $$u_r=-\frac{n}{1+n}{\left(-\frac{1}{m}\frac{\mathit{dp}}{\mathit{dr}}\right)}^{\frac{1}{n}}{z}^{\frac{n+1}{n}}+Cte$$(28)

To satisfy the no-slip condition at the walls, we have u_r = 0 at $z=+\frac{h}{2}$. Therefore:

(29) $$u_r=\frac{n}{n+1}{\left(\frac{h}{2m}\left(-\frac{\mathit{dp}}{\mathit{dr}}\right)\right)}^{\frac{1}{n}}\frac{h}{2}\left(1-{\left(\frac{2z}{h}\right)}^{\left(\frac{n+1}{n}\right)}\right)$$(29)

To derive the expression of the flow rate, the pressure profile as a function of r shall be determined using the flow continuity equation:

(30) $$\frac{d}{\mathit{dr}}\left(r{u}_r\right)=0$$(30)

Or

(31) $$r{\left(\frac{\mathit{dp}}{\mathit{dr}}\right)}^{\frac{1}{n}}=Cte$$(31)

We get:

(32) $$p=\mathit{\alpha }\frac{r^{1-n}}{1-n}+\mathit{\beta }$$(32)

Using the boundary conditions defined previously, the fluid pressure P in the region $z=\pm \frac{h}{2}$ and R_in < r < R_out is:

(33) $$P=P_{\mathit{in}}+\frac{P_{\mathit{in}}-P_{\mathit{out}}}{R_{\mathit{out}}^{1-n}-R_{\mathit{in}}^{1-n}}\left(R_{\mathit{in}}^{1-n}-r^{1-n}\right)$$(33)

We derive:

(34) $$-\frac{\mathit{dp}}{\mathit{dr}}=\left(1-n\right)\left(\frac{P_{\mathit{in}}-P_{\mathit{out}}}{R_{\mathit{out}}^{1-n}-R_{\mathit{in}}^{1-n}}\right)\frac{1}{r^n}$$(34)

The volumetric flow rate Q of liquid is obtained by integrating the velocity over a cross section of the valve:

(35) $$Q=2\underset{0}{\overset{\raisebox{1ex}{$h$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\int }}2\mathit{\pi }r{u}_{r}dz$$(35)

Combining (29), (34) and solving (35) yields:

(36) $$\begin{array}{c}Q=4\mathit{\pi }\frac{n}{n+1}{\left(\frac{h}{2m}\right)}^{\frac{1}{n}}\frac{h}{2}\\ \times {\left(1-n\right)}^{1/n}{\left(\frac{P_{\mathit{in}}-P_{\mathit{out}}}{R_{\mathit{out}}^{1-n}-R_{\mathit{in}}^{1-n}}\right)}^{\frac{1}{n}}\underset{0}{\overset{\frac{h}{2}}{\int }}\left[1-{\left(\frac{2z}{h}\right)}^{\left(\frac{n+1}{n}\right)}\right]dz\end{array}$$(36)

We obtain:

$$Q=4\mathit{\pi }\frac{n}{n+1}{\left(\frac{h}{2m}\right)}^{\frac{1}{n}}\frac{h}{2}$$

(37) $$\begin{array}{c}\times {\left(1-n\right)}^{\frac{1}{n}}{\left(\frac{P_{\mathit{in}}-P_{\mathit{out}}}{R_{\mathit{out}}^{1-n}-R_{\mathit{in}}^{1-n}}\right)}^{\frac{1}{n}}\left[\frac{h}{2}-\frac{h}{2}\frac{1}{\frac{n+1}{n}+1}{\left(\frac{2z}{h}\right)}^{\left(\frac{n+1}{n}\right)}\right]\end{array}$$(37)

Finally, the volumetric flow rate of power-law fluid in an annular valve can be expressed by:

(38) $$Q=4\mathit{\pi }\left(\frac{n}{2n+1}\right){\left(\frac{1-n}{m}\right)}^{\frac{1}{n}}{\left(\frac{P_{\mathit{in}}-P_{\mathit{out}}}{R_{\mathit{out}}^{1-n}-R_{\mathit{in}}^{1-n}}\right)}^{\frac{1}{n}}{\left(\frac{h}{2}\right)}^{2+\frac{1}{n}}$$(38)

Application to passive flow control valve design

The fluidic pathway of passive flow control valves as described in (Chappel et al., 2014) exhibits channels of circular cross-section and other types of flow restrictions.

The flow rate of power-law fluid in such cylindrical channel is, according to (Chhabra & Richardson, 2008):

(39) $$Q=\frac{n\mathit{\pi }}{3n+1}{\left(\frac{-\mathrm{\Delta }p}{2mL}\right)}^{\frac{1}{n}}{R}^{\frac{(3n+1)}{n}}$$(39)

where R is the channel radius. It can also be written in the following form:

(40) $$Q=\frac{\mathrm{\Delta }{p}^{1/n}}{R_{fcN}}$$(40)

where the term $R_{\mathit{fcN}}$ is not a function of the pressure:

(41) $$R_{fcN}=\frac{3n+1}{n\mathit{\pi }}\frac{{\left(2mL\right)}^{\frac{1}{n}}}{R^{(3+\frac{1}{n})}}$$(41)

The pressure drop induced by the viscosity in the whole device, considering several fluidic pathway i, is therefore:

(42) $$\mathrm{\Delta }{p}_i(viscosity)=Q_i^n\times \sum _j{\left(R_{fij}\right)}^n$$(42)

where j represents the different flow restrictions in each pathway i.

We introduce also singular head losses for the element j of each pathway i:

(43) $$\mathrm{\Delta }{p}_i(singular)=Q_i^2\times \sum _j{\mathit{\alpha }}_{\mathit{ij}}$$(43)

These losses are due to, for instance, sudden changes of channel section.

The total pressure drop through the device is therefore:

(44) $$\mathrm{\Delta }\text{P}=Q_i^2\times \sum _j{\mathit{\alpha }}_{\mathit{ij}}+Q_i^n\times \sum _j{\left(R_{fij}\right)}^n$$(44)

Since this equation has no general analytical solution for any n ≠ 1, numerical computations are used to estimate the terms Q_i by iterative method.