Design criteria and applications of multi-channel parallel microfluidic module

The flow resistance model in a single module was used for analysis, including the flow resistance among neighboring channel inlets (R_f), the resistance of each channel (R_u), the equivalent resistance of a single array (R_c), the flow resistance of vertical channels among neighboring arrays (R_ν), and the motivation in a fluid gravitational field (R_g), as shown in figure 1. In order to achieve uniform flow distribution in the 3D microfluidic system, R_u should be significantly greater than R_f while R_c should be also much greater than $\left| {{R}_{\nu }}-{{R}_{g}} \right|$. Based on the above rules, it is possible to derive the influencing parameters and the modular scale-up design criteria.

Zoom In Zoom Out Reset image size

The R_c is an important parameter which is a bridge that connects 2D parallel array and 3D modular scale-up. Mathematical and physical methods can be used to calculate R_c.

The design criteria of the 2D arrays in ladder geometry have been proposed and demonstrated for the uniform flow distribution [21–23]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 2N({{R}_{f}}/{{R}_{u}})\leqslant 0.01.\nonumber \end{align} \tag{ 1 }$$

In order to reveal the intrinsic relationship among scale-up influencing parameters, flow uniformity and microchannel structure, the above design criterion was also modified based on the ladder-like distribution. The two parameters δ and β were used to characterize the flow distribution performance. The smaller values of these two parameters are, the more excellent the flow distribution performance is, and the more uniform the particle size is. The equation (1) can be modified as equations (2) and (3) for the 2D array and 3D stack. N is the number of channels in the array and M is the number of arrays in the stack.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle N\left(\frac{{{R}_{f}}}{{{R}_{u}}} \right)=\delta \nonumber \end{align} \tag{ 2 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M\left(\frac{\left| {{R}_{v}}-{{R}_{g}} \right|}{{{R}_{c}}} \right)=\beta .\nonumber \end{align} \tag{ 3 }$$

The coefficient k was used instead of δ/N, and the R_f can be derived from equation (4):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}_{f}}=\frac{\delta {{R}_{u}}}{N}=k{{R}_{u}}.\nonumber \end{align} \tag{ 4 }$$

The mathematical formula was used to describe the 2D resistance model in figure 1, as shown in equation (5), where R_c,N and R_c,N−1 was the equivalent resistances of N channels and N  −  1 channels.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}_{c,N}}=k{{R}_{u}}+\frac{{{R}_{u}}{{R}_{c,N-1}}}{{{R}_{u}}+{{R}_{c,N-1}}} \nonumber \end{align} \tag{ 5 }$$

Solve the above general formula, and finally get mathematic relationship for R_c,N, as shown in equation (6):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}_{c,N}}=\left[ \frac{k}{2}+\left(\frac{{{f}^{N}}\left(k \right)+1}{{{f}^{N}}\left(k \right)-1} \right)\frac{\sqrt{{{k}^{2}}+4k}}{2} \right]{{R}_{u}}\nonumber \end{align} \tag{ 6 }$$

where $f\left(k \right)=\frac{k+2+\sqrt{{{k}^{2}}+4k}}{k+2-\sqrt{{{k}^{2}}+4k}}$.

$\xi =\left[ \frac{k}{2}+(\frac{{{f}^{N}}\left(k \right)+1}{{{f}^{N}}\left(k \right)-1})\frac{\sqrt{{{k}^{2}}+4k}}{2} \right]=f(\frac{\delta }{N},N)$ is the equivalent resistance coefficient of the array, then equation (6) can be simplified as shown in equation (7):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}_{c,N}}=\xi {{R}_{u}}=f(\frac{\delta }{N},N){{R}_{u}}.\nonumber \end{align} \tag{ 7 }$$

The calculation of the flow resistance R, is analogous to the basic principle of the circuit, since the flow of the fluid in the channel is similar to the transfer of the current in the wire. According to the ohmic formula R  =  U/I, the single channel resistance can be get as shown in equation (8), where R_channel is the flow resistance of the channel, dP is the pressure drop between the inlet and outlet, and Q is the flow rate of the channel.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}_{channel}}=\frac{dP}{Q} \nonumber \end{align} \tag{ 8 }$$

R_i is the resistance of a circular channel under laminar flow which is determined by channel diameter, length and fluid viscosity. We can use equation (9) to estimate it:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}_{i}}=\frac{128\,\mu {\rm l}}{\pi D_{h}^{4}}=\frac{8\,\mu {\rm l}}{\pi {{r}^{4}}}.\nonumber \end{align} \tag{ 9 }$$

It is found that R_c,N is related to ξ and R_u from the definition of equivalent resistance in the 2D array. Either of these two parameters increase will lead to the increment of R_c,N. The equivalent resistance coefficient ξ, which is an integrated function of N and δ/N, is related to 2D parallelization number N and the parameter of flow distribution performance δ. δ is mainly affected by the structure of the distribution, such as the parallel array or the circular array. Since ξ is a complex physical parameter, it needs to be discussed in detail. The relationship among ξ, N and δ in the 2D array is studied.

N directly affects the production rate. The higher N indicates the higher production rate. δ directly affects the uniformity of distribution. The smaller δ results in the higher uniformity. Therefore, in the 2D array, it is necessary to obtain a higher value of N and a smaller value of δ. To further study their internal relationship, we have got figure 2(created by Matlab) which is the curve family of ξ changing with δ/N and N. In the figure, the X-axis is δ/N while the Y-axis is ξ. For example, if we parallelize 8 channels in 2D array, N will be 8 and a purple curve will be selected. Different δ are selected to obtain different ξ. As can be seen from the figure, when N is small, the correlation among ξ, δ and N is obvious. Higher δ indicates the higher ξ. But increasing N will lead to a smaller ξ. With the increase of N, the curve tends to be gentle, and the influence of δ on ξ becomes lower, mainly depending on the value of N. When N equals 64 or more, these curves basically coincide and achieve convergent. Then ξ converges to 0.0733, and it has the lowest correlation with δ and N.

Zoom In Zoom Out Reset image size

Therefore, in the 2D array, if the value of N is high enough, the ξ is the minimum and the curve is gentle, regardless of N and the δ; if the value of N is limited, N and δ should be considered to determine ξ, and then calculate the value of equivalent resistance.

Substitute equation (3) of design criterion in 3D stack into the equation (7) of equivalent resistance. Equation (10) can be obtained:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\left| {{R}_{\nu }}-{{R}_{g}} \right|/{{R}_{u}}}{\beta /M}=\xi .\nonumber \end{align} \tag{ 10 }$$

Similar to N and δ of 2D array, M is the number of vertically stacked arrays. It directly affects the droplet throughput. The higher M will lead to the higher throughput. β is an indicator of vertical flow distribution performance. The smaller β value results in higher uniform distribution. Therefore, in the 3D stack, it is also necessary to obtain a higher value of M and a smaller value of β.

Equation (10) can be modified as equation (11), where ${{K}_{R}}=\left| {{R}_{\nu }}-{{R}_{g}} \right|/{{R}_{u}}$ is defined as the design parameter of the stack/module. K_R means the ratio of the resistance of vertical flow distribution channel to the resistance of single channel, which is the core in designing the single channels and the vertical flow distribution channels. The smaller K_R leads to the more conservative channel design. On the other word, it is much easier to obtain the higher throughput and better uniformity.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{{K}_{R}}}{\beta /M}=\xi \nonumber \end{align} \tag{ 11 }$$

Equation (11) reveals the relationship between ξ and the 3D stack. Figure 3 is the curve family of ξ changing with β/M and K_R(created by Matlab). It can be seen from the figure that, for the determined K_R, the smaller the ξ leads to the higher β/M. The smaller K_R, which results in the more conservative design, leads to the smaller β/M with the same ξ. Therefore, to obtain the same value of β, the number of stacked array M can be increased appropriately. When the value of K_R is conservative enough to be below 5  ×  10⁻⁵, several curves are basically coincided. Therefore, it is not necessary to get further conservative design. A higher level of processing difficulty is related to the more conservative design. When K_R and ξ is constant, for a fixed β/M value, the higher uniformity means that we should decrease β as well as M, or we can decrease the uniformity(β increases) to get a higher M.

Zoom In Zoom Out Reset image size

Equivalent resistance is an important parameter that relates 2D array and 3D stack, and ξ can also reflect the performance of two-stage scale-up. From the perspective of ξ, in the 2D array, it is desirable to increase N and lower δ to achieve higher output and better uniformity, but the result is to achieve a small ξ. When the ξ is small, the fixed K_R in the 3D stack will result in a higher β/M. It is actually detrimental because we need to get a smaller β/M to accommodate the higher throughput with the higher M and the higher uniformity with the smaller β. It can be seen that the requirements of two-stage scale-up for ξ are in conflict. According to equations (2) and (3), it can also be seen that increasing R_u can get smaller δ and β at the same time, which is good for two-stage scale-up. However, the improvement of R_u means the increased difficulty of mechanical processing and the channel structure design. In practice, combined with the mechanical processing and design conditions, we can determine R_u as much as possible. And then according to the parameter N, ξ can be determined from figure 2. β/M can be determined by certain K_R and ξ in figure 3. The values of β and M should be properly assigned.

The multi-stage module design criterion equation (12) can be got by combining equation (7) with (11):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{{K}_{R}}}{\beta /M}=f(\frac{\delta }{N},N).\nonumber \end{align} \tag{ 12 }$$

In this criterion, there are three categories of parameters, which are the throughput scale-up parameters M and N, the distribution uniformity parameters δ and β, and the microchannel design parameter K_R. The two throughput scale-up parameters are directly related to the whole throughput. Higher M and N lead to higher throughput. The two distribution uniformity parameters characterize the distribution performance of 3D multi-stage scale-up module. The smaller two values are, the more uniform the distribution performance is. And the microchannel design parameter is related to the machining geometric parameters. The smaller value indicates the more conservative design, but it will be subject to processing accuracy, space volume and other factors. We discussed these parameters separately based on the derived multi-stage module design criterion. Since δ and β have the same meaning, for convenience, we made these two values equal.

The relationship between throughput scale-up parameters and microchannel design parameters was discussed. Better uniformity can be achieved while setting the two uniformity parameters δ and β to 0.005. Figure 4 shows the curve family of M changing with different N and K_R while δ and β are fixed at 0.005(created by Matlab). For different K_R (K_R is respectively 2.5  ×  10⁻⁵, 5  ×  10⁻⁵, 7.5  ×  10⁻⁵, 1  ×  10⁻⁴, 2.5  ×  10⁻⁴, 5  ×  10⁻⁴, 7.5  ×  10⁻⁴, 1  ×  10⁻³.), the selected value of M, N is different. The higher K_R indicates that the curve is closer to the coordinate axes, and the values of M and N will be selected more strictly. Conversely, smaller K_R indicates the more conservative design, so throughput scale-up parameters may be appropriate to increase. For either case, the maximum of M can be determined from the figure 4 after selecting parameter N. Or according to the demand for throughput scale-up parameters, K_R can be determined by the graph, and then it is used to guide the microchannel design. It can also be found from the figure that when increasing R_u, considering that δ and β are fixed values, N and M will be further improved as K_R decreases.

Zoom In Zoom Out Reset image size

On the other hand, the relationship between throughput scale-up parameters and uniformity parameters was discussed. As shown in figure 5 (created by Matlab), K_R is set to be 2.5  ×  10⁻⁴. For different uniformity requirements, the values of M and N are also different. Figure 5 is the curve family of M changing with different uniformity parameters while K_R is fixed at 2.5  ×  10⁻⁴ (Both δ and β values are 0.02, 0.015, 0.01, 0.008, 0.006, 0.004, 0.002, and 0.001). As can be seen from the figure, when K_R is fixed, the higher uniformity requirements leads to a lower value of δ and β. Meanwhile, the value of M and N are selected more strictly and the increment of throughput is more limited. The appropriate increment of δ and β (flow distribution performance degradation) can greatly enhance the throughput. Based on the value of M and N, the distribution uniformity can be roughly determined.

Zoom In Zoom Out Reset image size

Figure 6 is production throughput (N  ×  M) changing with different δ, β and K_R. The Y-axis is the product of 2D parallelization number N and the 3D stack number M, which is N  ×  M. N is fixed at 8. The value of N  ×  M is directly related to the production throughput. The dispersed phase flow rate of single unit Q_cs multiplies by the total number of units N  ×  M, and Q_cs  ×  N  ×  M is the production throughput of a single module. For convenience, we set two distribution uniformity parameters δ  =  β since they have the same meaning. It can be quickly seen from figure 6 that the increment of throughput (N  ×  M) directly leads to an increment of the values of δ and β, which refers to the deterioration of the flow distribution performance. To improve the flow distribution performance, the values of δ and β need to be reduced, and then N  ×  M decreases. Therefore, it can be found that the demand for scale-up and the uniformity of flow distribution are in conflict. In addition, while δ and β are kept constant, smaller K_R will lead to the obvious increment of N  ×  M. The result is consistent with the discussion in section 3.2. Smaller K_R reflects the improvement of R_u (resistance of each channel), which is related to the mechanical processing and the channel structure design.

Zoom In Zoom Out Reset image size

The above analysis was based on the ideal condition, and the flow resistance model based on ladder-shaped asymmetric structure was deduced without considering the influence of the machining error. Therefore, it is necessary to modify the above derivation process. Two basic criteria equations (2) and (3) are modified to equations (13) and (14) respectively:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle N\left[ \frac{\left({{R}_{f}}+\Delta {{R}_{u}} \right)}{{{R}_{u}}} \right]=\delta \nonumber \end{align} \tag{ 13 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M\left[ \frac{(\left| {{R}_{v}}-{{R}_{g}} \right|+\Delta {{R}_{u}})}{{{R}_{c}}} \right]=\beta , \nonumber \end{align} \tag{ 14 }$$

where $\Delta {{R}_{u}}$ is the error of microchannel resistance ${{R}_{u}}$ in machining process, which is the main influencing factor of flow distribution. Through resolving equations (13) and (15) can be obtained.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{{R}_{f}}}{{{R}_{u}}}+\frac{\Delta {{R}_{u}}}{{{R}_{u}}}=\frac{\delta }{N} \nonumber \end{align} \tag{ 15 }$$

The left side of equation (15) is composed of two terms, $\frac{{{R}_{f}}}{{{R}_{u}}}$ and $\frac{\Delta {{R}_{u}}}{{{R}_{u}}}$, respectively, revealing two key issues of fluid distribution. $\frac{{{R}_{f}}}{{{R}_{u}}}$ is the structural issue, which refers to the asymmetry factor of the branch structure. And it mainly reflects the pressure difference among neighboring channel inlets. $\frac{\Delta {{R}_{u}}}{{{R}_{u}}}$ is non-structural issue, which refers to machining process error and it is not related to the structure. To meet the uniformity of flow distribution, some conditions should be achieved:

Structural condition: the resistance among neighboring channel inlets R_f should be far less than the resistance of each channel R_u so as to ensure that the inlet pressure of each channel is nearly the same.

Non-structural conditions: the error of the resistance of each branch is negligible compared with R_u so as to ensure that each R_u is nearly the same.

Equation (14) can also be split in a similar way, as shown in equation (16):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\left| {{R}_{v}}-{{R}_{g}} \right|}{\xi {{R}_{u}}}+\frac{\Delta {{R}_{u}}}{\xi {{R}_{u}}}=\frac{\beta }{M}.\nonumber \end{align} \tag{ 16 }$$

When achieving the completely symmetrical ladder-like structure, the distribution channel evolves into a central reservoir. And the pressure of each branch is almost the same, which means reducing R_f in the maximum extent and reaching the limit of structural condition. Therefore, the distribution uniformity mainly depends on the processing error of the microchannel. In order to meet the non-structural conditions, the serpentine microchannel needs to be adopted to improve the resistance of each branch as shown in figure 7, so that the flow uniform distribution can be conveniently achieved. As the length of the serpentine channel increases, the distribution performance will also gradually increase. Our research team has presented a flow distributor based on polymethylmethacrylate (PMMA) and laser microfabrication engraving (Universal VLS2.30) [24], which uses serpentine channel [26]. Serpentine channel can balance the pressure drop at each branch inlet and the non-uniform flow resistance caused by machining process to improve flow distribution performance [27–29].

Zoom In Zoom Out Reset image size

According to the design criteria proposed in this paper, the design and manufacture of scale-up module have been carried out. In the discussion of the previous section, a long serpentine channel was adopted to meet the non-structural condition, and a circular array was adopted to meet the structural condition.

The experiment, taking into account the processing equipment and structure factor, selected 80 mm  ×  80 mm polymethylmethacrylate (PMMA) plate to array 8 channels. Therefore N  =  8 was selected. It has been studied that in 2D scale-up criterion, 2δ should be less than 0.01 so that flow distribution performance is better, hence δ  =  0.005 is selected [21–23]. As discussed in section 3.1 and shown in figure 2, ξ  ≈  0.13 was obtained.

The above parameters belong to 2D scale-up parameters, and then 3D scale-up parameters were discussed. Due to the limitation of PMMA plate size and processing precision, the total length of the serpentine channel was set to be approximately 280 mm, and the channel section was a small triangle with the size of 619 µm which was viewed by an optical microscope with a digital camera (OLYMPUS SZX7). Based on the actual processing size, the R_u value can be calculated, and then K_R can be got. Through the above parameters calculated, it shows that β/M is a certain value by equation (11) and figure 3. From figure 1, the 3D scale-up can also be reduced to the 2D scale-up. Therefore, we set β  <  0.005. When β is less than or equal to 0.005, the performance of flow distribution in the 3D scale-up is better. Figure 8 is the value of β changing with M. When M  =  10 is set, β is just less than 0.005 while β is 0.005 197 with M  =  11, meeting the flow distribution uniformity requirements. According to the above steps, we determined N  =  8, M  =  10.

Zoom In Zoom Out Reset image size

The relevant verification experiments have been carried out by our research team. According to the design criteria, a scale-up module has been manufactured. Taking chitosan system as an example, the chitosan/TiO₂ composite microspheres were prepared and the particle size, surface morphology, CV Values were also analysed. The detailed experiment has been reported by our research team [24]. The chitosan powder and TiO₂ powder were dissolved in acetic acid (Shanghai Lingfeng Chemical Reagents Co., Ltd) and deionized water, forming an aqueous solution as the dispersed phase. The continuous phase was made from a mixture of paraffin liquid and petroleum ether with volume ratio of 7:5 (all from Shanghai Aladdin Bio-Chem Technology Co., LTD), and 5% span 80 was added as stabilization for chitosan droplets. The functional microspheres prepared by chitosan have good adsorption properties and they can be widely used in industrial wastewater treatment, ion adsorption and other fields [30–34].

The experiment schematic and the photograph of experimental equipment are shown in figures 9(a) and (b). The preparation process of chitosan system is shown in figure 9(c). The specific structure can refer to the patent [35]. Before the new test started, the necessary parameters could be adjusted, such as solution composition, experimental temperature and so on. The experimental parameters should be kept constant during the test. Between run to run, the variations that could be controlled were the flow rate of the dispersed phase and continuous phase. The matching of the two-phase flow determines the size of droplets. The scale-up module was run at a continuous phase flow rate Q_c of 80 ml min⁻¹ and a dispersed phase flow rate Q_d of 8 ml min⁻¹ with two-phase flow rate ratio Q_c/Q_d of 10. The production throughput of a single module was approximately 12 l d⁻¹. 7.3  ×  10⁵ droplets could be generated per hour in a single module. Figure 10 shows the particle size of the microspheres prepared by the module. The number on the X-axis means the 2D scale-up array from top to bottom. Due to the vertical stacking of the arrays, Q_c increases while Q_d decreases with the particle diameter of each array generally decreasing from the first array to the tenth array. The average diameter of chitosan microspheres is 539.65 µm and the CV value is about 3.59%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Obviously, many modules can be integrated to build a microfluidic scale-up system. The flow distribution error is related to the processing error of each module and the flow rate error of fluid driving devices.