Dynamic Modeling and Flow Distribution of Complex Micron Scale Pipe Network

A fluid simulation calculation method of the microfluidic network is proposed as a means to achieve the flow distribution of the microfluidic network. This paper quantitatively analyzes the influence of flow distribution in microfluidic devices impacted by pressure variation in the pressure source and channel length. The flow distribution in microfluidic devices with three types of channel lengths under three different pressure conditions is studied and shows that the results obtained by the simulation calculation method on the basis of the fluid network are close to those given by the calculation method of the conventional electrical method. The simulation calculation method on the basis of the fluid network studied in this paper has computational reliability and can respond to the influence of microfluidic network length changes to the fluid system, which plays an active role in Lab-on-a-chip design and microchannel flow prediction.


Introduction
Nowadays, the rapid development of science and technology constantly drives the upgrading of scientific and technological products. The portability, integration and intelligentization of scientific and technological products are the key points and difficulties in current science and technology product development. Meanwhile, with the improvement of living standards, health issues have become the focus of people's attention so that people put forward higher standards for the accuracy of test results and detection methods. However, the real-time detection and rapid diagnosis technology in medical and health fields are currently less developed and need to be improved urgently. Therefore, the microfluidic chip technology arises at the right moment.
Microfluidic chip technology refers to building biological or chemical laboratories on chips that are only a few square centimeters in size. Different types of fluid channels, whose diameters are in the dozens to hundreds of micrometers range, are built on the microfluidic chip based on various functions. These channels have many utilities, such as biological or chemical reactions and extraction or detection of reagents. Meanwhile, the microfluidic transmission process can be controlled through the design of the microchannel network [1].
The main functions of the microfluidic system are transferring and mixing, reaction separation, and flow control of microfluidics in the microchannels. Because of the small size of the fluid channel, the fluid flowing in the microchannel is mostly in the laminar flow state, such as particles, droplets or bubbles, which generally belong to the field of low Reynolds number fluid theory in microchannels [2].
In order to calculate and control the flow of the microfluidic network more accurately, the study method that is commonly used by researchers is the electrical equivalent method [3]. By using the electrical equivalent method to liken the micropipe to an electrical device, the microfluidic network can be integrated and analyzed. However, it cannot simulate the issues of the local resistance well in the microfluidic pipe, and the diverting 2 of 12 resistance in the microfluidic network needs to be calculated by mathematic methods to be more accurate [4]. This manner aforementioned is the calculation method of a fluid network which is on the basis of the three conservation equations of basic fluid mechanics, including the conversation of mass, momentum and energy [5]. Pressure value and flow rate can be obtained by simultaneous calculations of mass and momentum and heat is solved by the energy equation. However, this study method is rarely reported at present.
Recently, the flow characteristics of microfluidic networks have been studied by many scholars. However, most of them focus on computational fluid dynamics (CFD) and electrical analogy methods [3,6]. For example, Ali Y. Alharbi's group studied the flow process of fractal-like branching networks in micropipes and the pressure drop characteristics and force mode in branching microchannels using three-dimensional computational fluid dynamics approaches [7]. Dalei Jing and Yongping Chen et al. also studied the heat and mass transfer in branched microfluidic channels, but the equivalent methods of microfluidic resistance are all hypotheses of using circuits [8,9]. Bassiouny's group [10,11] developed analytical equations to predict the fluid flow distribution in rectangular, U-shaped and Z-shaped manifolds. Kim's group [12] analyzed the influence of a Z-manifold on channel flow distribution and studied three different shapes, including rectangle, trapezoid and triangle. Delsman's group [6] performed a numerical analysis of the flow distribution in a multi-channel microfluidic device with Z-shape and inline manifolds. They obtained the results of multiple manifold designs by changing the shape and flow direction of the device's inlet and outlet. Wang's group [13] took friction and momentum effects into account in the manifold and modified the U-shaped manifold model established by Bassiouny and Martine; they found that the effects of friction and momentum could increase and decrease the pressure drop in the manifold, respectively. Therefore, it is observed that a uniform distribution of fluid flow can be achieved by balancing these two opposite effects properly.
In this paper, we use the fluid network calculation method to analyze the flow characteristics of the microfluidic network and the conservation equation of mass and momentum to solve the pressure and flow in the microfluidic network. The method that we used was on the basis of the structural characteristics and flow process of the actual network. The flow distribution in a microfluidic device with five microchannels is analyzed and compared with the conventional electrical method. The fluid network method takes the formula of resistance term in the calculation program into account, and its calculation results are closer to the actual flow process than those using the conventional electrical method.

Model Development
In this study, it is assumed that the fluid flowing in the microfluidic network is water and that there is no phase change during transportation. Therefore, it can be considered a single-phase incompressible fluid. The hypotheses made to establish the model in this study can be expressed as follows:  For a stable microfluidic network, the continuity equation (continuous mass) in the flow process can be expressed as follows: Without considering the influence of heat from dρ dt = ∂ρ ∂p dp dt + ∂ρ ∂H dH dt ≈ ∂ρ ∂p dp dt , the following equation can be obtained: is the compressibility (kg/MPa) of the working fluid; p is the pressure of the working fluid; ρ is the density of the working fluid; t is the time; N is the total number of nodes; G ij is the mass flow between nodes i and j; V i is the volume of node i; D ij is the connection mode between nodes i and j (i = 1, 2, . . . , N; j = 1, 2, . . . , N) and the specific meaning is: there is a connection between nodes i and j, and the flow direction is from j to i 0, there is no connection between nodes i and j -1,there is a connection between nodes i and j, and the flow direction is from i to j The continuity equation is the flow conservation relationship at a node, and the flow state of the node is affected by the nodes connected to it. Figure 2 shows a schematic diagram of the volume elements of all nodes which are connected to i. It is important to note that a node can be connected to N nodes logically. However, a node to connect the number of other nodes is limited in the actual model. The momentum conservation equation in the microfluidic network can be expressed as follows: where, H ij is the pressure generated by macro kinetic energy, potential energy and power source between node i and j; U ij is the fluid velocity between node i and j; D ij has the same meaning as D ij in Equations (1) and (2).
In the formula, h w = ∑ 0.5λ L d + ξ is the on-way and local resistance loss [14]; λ is the on-way resistance coefficient; L is the length of pipe; d is the diameter of the pipe; ξ is the local pressure loss.
According to the pipeline inertia coefficient 2 , and no power supply in the pipeline, H ij = 0 is introduced into the Equation (3): R ij is the resistance characteristic term between node i and j, which can be expanded as follows: where, λ is the equivalent resistance coefficient on the way of the microchannel; ξ is the local resistance coefficient of the microchannel; ρ is the working medium density; A is the cross-sectional area of the microchannel. Equations (2) and (4) are nonlinear differential and algebraic equations. For dynamic simulation, the priority is to ensure real-time performance, which requires high robustness of the calculation process. Using an implicit Euler integration algorithm to calculate, the equation can be obtained by equations (2) and (4) as follows: where, r is the discrete-time variable, time is t = t 0 + r∆t, r = 0, 1, 2, . . . Equation (6) is completely decoupled. The pressure of each node can be obtained by solving the N- order nonlinear algebraic equations composed by equation (6) in the process of dynamic simulation. Under the slight perturbation condition, C t+1 i can be substituted by C t i , which has little influence on precision and can reduce the calculation amount. After solving the pressure of each node, the momentum equation is introduced to obtain the branch flow. The calculation process is shown in Figure 3.

Model Design and Electrical Equivalent
The basic structure of the microfluidic system model includes micropump, micromixer, microvalve, and microchannel with appropriate size. The flow characteristics of the fluid in the chip can be analyzed at the micron or even nanometer level through these unique miniaturized structures and control devices.
In this case, taking the microscopic scale issues of fluid flow into account, a fluid flow model of an incompressible fluid with a low Reynolds number is adopted to study pipeline pressure distribution under microscale conditions and flow distribution characteristics in the process of fluid flow. Therefore, we designed a microfluidic flow distribution model that has a single entry with three exports, such as Figure 4, as a simulation tool. The variation of fluid flow value in different pipe segments under different inlet pressures was observed and compared with the results obtained by using the electrical method. Since there is no mixing of the two fluids, in this case, the micromixer is not considered in the design of the structure, and the position of the microvalve does not affect the final equivalence study of the fluid flow distribution. Therefore, the electrical analogy method can treat microchannels as resistors.  Figure 4 shows the specific device design of the microfluidic model. The first is the drive system, which consists of a pressurized chamber, a microchannel controller, and a microvalve that activates the pulse. According to the actual characteristics of our pipe network, we choose the pressure pump as the driving device, whose main role is to provide driving pressure for the microfluidic network. The size of the pressure chamber can be determined on the basis of the size of the microvalve and microchannel used. In the electrical model, this is a power supply that provides voltage. The pressurized chamber is one of the ways to provide the driving pressure. Through the combination of the pressure chamber and the micro valve as a controllable pressure source, the driving pressure can be adjusted according to the system demand.
In order to verify the accuracy of our algorithm, we have mentioned previously the conventional electrical equivalence method in which the microfluidic network is equivalent to a circuit model and the fluid flow parameters are equivalent to electrical parameters. According to the electrical abstraction of our model, Figure 4, the microfluidic network with circuit characteristics can be obtained, as shown in Figure 5. Where, R F1 is the microchannel resistance of the microvalve output, R F2 is the main channel resistance of fluid flow, and R F3 , R F4 , R F5 is the branch channel resistance of fluid flow, respectively. Corresponding to the resistance of the storage tank where the fluid is located and taking the processing technology of the microfluidic chip into account, the cross-section of the microfluidic tube is generally rectangular, and its resistance is [15]: where, L is the length of the microchannel, h is its height, w is its width, and u is the viscosity of the fluid. Because of the greater sizes of reservoir and outlet, their resistances to remaining resistors are ignored. Excluding the influence of structure and considering only the logical model of the circuit, the equivalent circuit can be obtained, as shown in Figure 6.
cosity of the fluid. Because of the greater sizes of reservo remaining resistors are ignored. Excluding the influence o the logical model of the circuit, the equivalent circuit can 6.   P 0 (t) represents the pressure in the supercharging chamber, that is, namely the pressure value at the inlet of the channel; P 1 (t) and P 2 (t) corresponding to the pressure value at the node of the channel; P 3 (t), P 4 (t), and P 5 (t) corresponding to the pressure value at the outlet of the channel and atmospheric pressure is considered. Q 1 (t) and Q 2 (t) are the fluid flows in the main channel, and Q 3 (t), Q 4 (t), and Q 5 (t) are the fluid flows in each branching channel.
In this paper, we use Simulink to build the model because of using the circuit analogy method to analyze the parameters of the fluid network in the microchannel. Simulink is a graphical modeling tool in MATLAB software, which can build and calculate various circuit models. The circuit analogy built by Simulink satisfies basic electrical laws, such as Ehrhoff's and Ohm's law. According to the structure in Figure 6, the relevant power parameter relationship can be obtained, and the specific expression is as follows:

Model Calculation Results
According to the electrical method, the equivalent resistance and the dimensions of each branching channel in this example are shown in Table 1. The main local drag coefficient in this model is brought about by the micro three-way pipe. For the time being, we do not consider the local resistance coefficient ξ, but only consider the resistance of the micropipe, which is compared with the circuit model of the electrical method. In order to reduce the activation energy required for fluid flow, the microvalve needs to have a higher length-width ratio. A PDMS film deformation valve can be selected, and both the height of the valve and the main channel are designed to be 500 µm, whose main role is to control the inlet pressure of the fluid network. The driving pressure here is the network pressure, so the micro-valve resistance does not need to be considered. Meanwhile, considering the uniformity of the fluid flow in the microchannel, the distance of the fluid movement should be as long as possible, and the cross-sectional area design of the fluid flow channel should be narrow. Therefore, the length percentages of channels 3, 4, and 5 are set to 5:3:3 and the height is set to 350 µm.
Above all, we keep other parameters constant and change the driving pressure, P 0 (t), in the microfluidic network, set as 15 kPa, 18 kPa, and 20 kPa, respectively. In the outlet boundary condition of the microchannel, the pressure values, P 3 (t), P 4 (t) and P 5 (t), are set as 11 kPa, 12 kPa and 11 kPa, respectively. The results of the fluid network algorithm and electrical simulation of Simulink are shown in Table 2. Both the results indicate that the flow in the microfluidic network satisfies the relationship of flow distribution. The flow rate of Q 2 in the main microchannel is the sum of Q 3 and Q 5 . The flow rate of Q 1 is the sum of Q 2 and Q 4 . Moreover, there is a certain nonlinear relationship among the microchannels when the input pressure of the pressure pump increases. The electrical method is very similar to the results of the fluid network method we used, which proves the correctness of the fluid network method and can be used to calculate the microfluidic network. The first example shows that the total flow rate of the system, Q 1 (t), increases with the increase of the input pressure. Microchannel 3 is still the outlet with a smaller flow rate, due to the longer length of microchannel 3, so its equivalent resistance value is greater than that of microchannel 4 and 5. Afterward, the increase of pressure leads to the flow rate rise, which is mainly reflected in the shunting of micro-channel 4 and 5. While the flow rate increase in micro-channel 3 is relatively lower. Both calculation results show this point, as shown in Figure 7, and can prove that the computational results of the fluid network method are closer to those of the electrical method, which has certain computational accuracy. In the second example, in order to study the influence of microchannel structure on flow characteristics, the width and depth of microchannel 4 were kept at 350 µm. The flow distributions at different channel lengths (40 mm, 50 mm and 60 mm) were studied. The boundary condition at the inlet was set as 15 kPa and the three outlet pressures, P 3 (t), P 4 (t), and P 5 (t), are set as 11 kPa, 12 kPa, and 11 kPa, respectively. The calculation results of the two methods are shown in Table 3 and Figure 8. The diagram shows that the flow rate in microchannel 4 decreases with the increase of its length, which is because the increase of the length of the microchannel leads to the resistance value rise. On the other hand, comparing microchannel 3 with microchannel 5, it can be found that the flow rate in these two channels is relatively stable. Under certain pressure conditions, the flow rate change caused by changing the structure of microchannel 4 does not affect the flow rate in microchannel 3 and 5.  In the first example in this section, we have verified the calculation accuracy of the fluid network algorithm and electrical methods. Compared the results, which are similar. However, the resistance items of the fluid network algorithm can be more detailed, such as three or more pipe diversion coefficients of local resistance, and is a good way to add items in resistance, R ij , in the model. Electrical methods can treat micro-channels or micro-valves as resistance items, but it is difficult to deal with local resistance items such as three-way or multi-way, and some approximate formulas are generally adopted to calculate it [15]. According to Figures 7 and 8, we know that the feedback results of the fluid network method and electrical method to the pressure and structure changes in the microfluidic network are similar, which can be used for the simulation of the fluid flow process in the microfluid network and set up test-rig in the follow-up work. According to the actual measurement parameters, the calculation model can be modified to improve its accuracy.

Conclusions
In this paper, a computational method of microfluidic pipeline network on the basis of fluid network is studied and compared with the conventional electrical equivalent method. It is found that the computational results obtained by using the fluid network method are closer to those given by using the electrical method. Therefore, we believe that the computational accuracy of the fluid network method is the same as that of the electrical method. From the research results of the simulation examples established in this paper, it can be concluded that the fluid network method can reflect the fluid flow process under different pressures, and the effects caused by structural changes in the microfluidic network are consistent with the results obtained by the electrical equivalent method. However, the electrical method has limitations in the treatment of local resistance, which is difficult to deal with the network shunt, the gradual shrinkage or expansion of microchannels, etc. that can be solved in the algorithm of fluid network easily. Furthermore, the accuracy of its calculation can be improved by merely modifying the resistance term in the algorithm. The computational method of the microfluidic network that we proposed in this paper can be used as an alternative to the electrical method. The follow-up work is to study the resistance term of the model and how to deal with the various devices presented in microfluidic is an important research direction.