*Article* **Forced Convection in Wavy Microchannels Porous Media Using TiO<sup>2</sup> and Al2O3–Cu Nanoparticles in Water Base Fluids: Numerical Results**

**Kholoud Maher Elsafy and Mohamad Ziad Saghir \***

Department of Mechanical and Industrial Engineering, Ryerson University, Toronto, ON M5B 2K3, Canada; kelsafy@ryerson.ca

**\*** Correspondence: zsaghir@ryerson.ca

**Abstract:** In the present work, an attempt is made to investigate the performance of three fluids with forced convection in a wavy channel. The fluids are water, a nanofluid of 1% TiO<sup>2</sup> in a water solution and a hybrid fluid which consists of 1% Al2O3–Cu nanoparticles in a water solution. The wavy channel has a porous insert with a permeability of 10 PPI, 20 PPI and 40 PPI, respectively. Since Reynolds number is less than 1000, the flow is assumed laminar, Newtonian and steady state. Results revealed that wavy channel provides a better heat enhancement than a straight channel of the same dimension. Porous material increases heat extraction at the expenses of the pressure drop. The nanofluid of 1% TiO<sup>2</sup> in water provided the highest performance evaluation criteria.

**Keywords:** porous cavity; wavy channels; nanofluids; forced convection; heat enhancement; pressure drop; mesh model

#### **1. Introduction**

Thermal engineering is one of the major fields in which energy conservation and sustainable development demands are increasing requiring more research efforts towards more energy efficient equipment and processes. The petroleum and process industries helped in many ways through the years to enhance the efficiency and find promising solutions [1]. A new era of microelectronic applications and devices has started requiring more thermal management efficiency specially with the obvious increase of the heat flux of chips. Microchannels heat sinks have been investigated firstly by Tuckerman and Pease [2] and a lot after because of the excellent heat transfer efficiency and capacity plus the small dimensions that are needed nowadays in most new technological applications.

Flow and heat transfer from irregular surfaces are often encountered in many engineering applications to enhance heat transfer such as micro-electronic devices, flat-plate solar collectors and flat-plate condensers in refrigerators, geophysical applications, underground cable systems, electric machinery, cooling system of micro-electronic devices, etc. In addition, roughened surfaces could be used in the cooling of electrical and nuclear components where the wall heat flux is known [3]. An extensive amount of research has been performed in the last decades to obtain a better understanding of flow mixing and heat transfer enhancement in channels with geometrical inhomogeneities, such as serpentine channels, asymmetric and symmetric wavy walls, natural convection heat transfer in wavy porous enclosures, wavy microchannels, etc. [3–6]. In recent years, it has been shown that nanofluids can be applied in heat exchangers to enhance the heat transfer, leading to higher heat exchanger efficiency [7].

The advantages of wavy channels are the ease of manufacturing and the significant enhancement of heat transfer as it was operated [8,9]. Direct liquid cooling incorporating microchannels is considered one of the promising solutions to the problem [10,11]. Secondary flow (Dean vortices) is generated when liquid coolant flows through the wavy

**Citation:** Elsafy, K.M.; Saghir, M.Z. Forced Convection in Wavy Microchannels Porous Media Using TiO<sup>2</sup> and Al2O3–Cu Nanoparticles in Water Base Fluids: Numerical Results. *Micromachines* **2021**, *12*, 654. https:// doi.org/10.3390/mi12060654

Academic Editors: Junfeng Zhang and Ruijin Wang

Received: 5 May 2021 Accepted: 31 May 2021 Published: 2 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

microchannels. It is found that along the flow direction, the quantity and the location of the vortices may change leading to chaotic advection, which can greatly enhance the convective fluid mixing.

Choi [12] was the first to mention the term nanofluid which refers to any liquid that contains solid metallic particles in submicron scale/nano scale (e.g., Ag, TiO2, Cu or Al2O3). These nanoparticles then are added to water to form a nanofluid which was found to increase the thermal conductivity and enhance the heat transfer performance of the base working fluids [13]. Widely speaking, it was found from different results that the type of nanoparticles used in the base fluid is the main factor to enhance the performance of heat transfer of traditional work fluids such as Ag nanofluid which gave experimentally the higher performance value due to the inherent properties of this metal oxide that overcomes the thermal diffusion that happens due to the small sizes of its particles [14,15]. Lee et al. [16] used the base fluids such as water and ethylene glycol and added nanoparticles of Al2O<sup>3</sup> and CuO to measure the thermal conductivity. The results obtained by them showed that ethylene glycol based nanofluid achieved a higher value of thermal conductivity than the water-based type. Xuan and Roetzel [17] were also studying the performance of heat transfer using nanofluids and had some interesting results. They were able to derive many equations using single and double-phase flow techniques for the analysis of convective heat transfer in nanofluid. The heat transfer enhancement of nanofluids and the laminar flow were investigated by Heris et al. [18,19] in a circle tube, and the results were typical compared to the other research. Flow and heat transfer characteristic of copper-water nanofluid in a two-dimensional channel was studied by Santra et al. [20] using Cu water nanofluid. The results were notable as with the increase of volume fraction of the solid nanoparticles and Reynolds number, the heat transfer rate was found to be higher.

Saghir research team [21–32] have experimentally and numerically investigated the forced convection of Al2O3–Cu hybrid nanofluid, Al2O3/water nanofluid in porous media at different flow rates and heating conditions in a straight channel. Results revealed that fluid such as water with metallic nanoparticle can enhance the heat extraction by no more than 6% when compared to water circulation. They concentrated their effort in straight channel with identical dimension.

In the present paper, attempt is made to investigate the performance of wavy channel using the same geometrical parameters and heating condition as the straight channel's cases done by Saghir et al. The novelty of this research is to be able to enhance the heat removal by allowing the flow to circulate a longer pathway when compared to rectangular case. Pressure drop is a main concern for engineering application thus the need to evaluate the Performance Enhancement Criteria. Section 2 presents the problem description, followed by the finite element formulation in Section 3. Section 4 is allocated for the case studied and the conclusion is in Section 5.

#### **2. Problem Description**

In this work, we are investigating a numerical approach of forced convection of nanofluids and hybrid fluid in three-wavy porous channels configuration as shown in Figure 1. The objective is to investigate the heat transfer performance of these fluids that could act as coolants. The three fluids used in the current analysis are water, titanium dioxide/water based nanofluid (1% TiO2/water) and aluminum oxide/copper nanoparticle/water-based hybrid fluid (1% Al2O3–Cu/water). For the hybrid fluid, the single nanoparticles are composed of 90% Al2O<sup>3</sup> and 10% copper [27].

**Figure 1.** (**a**) Model setup. (**b**) Channels' block.

u<sup>୧୬</sup> T<sup>୧୬</sup> − − − − − − − The numerical setup model consists of an inlet tube, a mixing chamber, three wavy porous channels insert, an exit chamber and an outlet tube. The described setup is placed over a heated aluminum block that represents the hot surface. The heated surface of the aluminum block is in direct contact with the bottom of the three channels block. The channel dimensions have a width of 0.00535 m, a height of 0.0127 m and a length of 0.0375 m. The heated aluminum block dimensions are 0.0375 m × 0.0375 m × 0.0127 m. The fluid enters the system with a specific velocity uin and temperature Tin. The temperature was calculated numerically along the fluid path 1mm below the interface of the aluminum block at the centre as shown in Figure 1a. Different flow rates were applied corresponding to 0.05 US gallon per minutes (USGPM), 0.1 USGPM, 0.15 USGPM and 0.2 USGPM. This corresponds to a flow rate of 3.15 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>m</sup>3/s, 6.3 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>m</sup>3/s, 9.45 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>m</sup>3/s and 1.26 <sup>×</sup> <sup>10</sup>−<sup>5</sup> <sup>m</sup>3/s, respectively. The wavy channels are assumed porous, and the permeabilities used were 10 pore per inches (10 PPI), 20 PPI and 40 PPI. This corresponds to a permeability of 9.557 <sup>×</sup> <sup>10</sup>−<sup>7</sup> <sup>m</sup><sup>2</sup> , 2.38 <sup>×</sup> <sup>10</sup>−<sup>7</sup> <sup>m</sup><sup>2</sup> and 3.38 <sup>×</sup> <sup>10</sup>−<sup>8</sup> <sup>m</sup><sup>2</sup> , respectively. The porosity is maintained constant at 0.91. The inlet pipe diameter is set equal to 0.01 m, and the heat flux applied to the model, as shown in Figure 1a, has an intensity of 75,000 W/m<sup>2</sup> . Table 1 presents the physical properties of the fluid used in our simulation. The main reason for selecting these fluids is that we have conducted experimental measurement with these fluids in a straight channel configuration, and no sedimentation of the nanoparticles has been observed.


**Table 1.** Thermo-physical properties of the fluid used in the analysis [27–29].

#### **3. Governing Equation and Boundary Conditions**

In the present work, we attempt to solve the Navier-Stokes equation for the fluid in the entrance and exit chamber combined with the Brinkman formulation for the flow in the porous channels and the energy equation for the fluid in the setup. In addition, the heat conduction equation is solved for the solid surface. The problem is assumed steady state and the flow is in laminar regime. The set of equations used in our model is as follows:

Momentum equations along x, y and z directions, respectively,

$$\rho\_{\rm nf} \left( \mathbf{u} \frac{\partial \mathbf{u}}{\partial \mathbf{x}} + \mathbf{v} \frac{\partial \mathbf{u}}{\partial \mathbf{y}} + \mathbf{w} \frac{\partial \mathbf{u}}{\partial z} \right) = -\frac{\partial \mathbf{p}}{\partial \mathbf{x}} + \mu\_{\rm nf} (\frac{\partial^2 \mathbf{u}}{\partial \mathbf{x}^2} + \frac{\partial^2 \mathbf{u}}{\partial \mathbf{y}^2} + \frac{\partial^2 \mathbf{u}}{\partial z^2}) \tag{1}$$

$$\rho\_{\rm nf} \left( \mathbf{u} \frac{\partial \mathbf{v}}{\partial \mathbf{x}} + \mathbf{v} \frac{\partial \mathbf{v}}{\partial \mathbf{y}} + \mathbf{w} \frac{\partial \mathbf{v}}{\partial z} \right) = -\frac{\partial \mathbf{p}}{\partial \mathbf{y}} + \mu\_{\rm nf} (\frac{\partial^2 \mathbf{v}}{\partial \mathbf{x}^2} + \frac{\partial^2 \mathbf{y}}{\partial \mathbf{y}^2} + \frac{\partial^2 \mathbf{v}}{\partial z^2}) \tag{2}$$

$$\rho\_{\rm nf} \left( \mathbf{u} \frac{\partial \mathbf{w}}{\partial \mathbf{x}} + \mathbf{v} \frac{\partial \mathbf{w}}{\partial \mathbf{y}} + \mathbf{w} \frac{\partial \mathbf{w}}{\partial z} \right) = -\frac{\partial \mathbf{p}}{\partial \mathbf{z}} + \mu\_{\rm nf} \left( \frac{\partial^2 \mathbf{w}}{\partial \mathbf{x}^2} + \frac{\partial^2 \mathbf{w}}{\partial \mathbf{y}^2} + \frac{\partial^2 \mathbf{w}}{\partial z^2} \right) \tag{3}$$

Continuity equation,

$$(\frac{\partial \mathbf{u}}{\partial \mathbf{x}} + \frac{\partial \mathbf{v}}{\partial \mathbf{y}} + \frac{\partial \mathbf{w}}{\partial z}) = \mathbf{0} \tag{4}$$

Energy conservation equation,

$$(\rho \,\mathrm{Cp})\_{\mathrm{nf}} (\mathbf{u} \frac{\partial \mathbf{T}}{\partial \mathbf{x}} + \mathbf{v} \frac{\partial \mathbf{T}}{\partial \mathbf{y}} + \mathbf{w} \frac{\partial \mathbf{T}}{\partial z}) = \,\mathrm{k}\_{\mathrm{nf}} (\mathbf{u} \frac{\partial^2 \mathbf{T}}{\partial \mathbf{x}^2} + \mathbf{v} \frac{\partial^2 \mathbf{T}}{\partial \mathbf{y}^2} + \mathbf{w} \frac{\partial^2 \mathbf{T}}{\partial z^2}) \tag{5}$$

For the porous flow, the following formulation are used. In particular:

$$\frac{\mu\_{\rm nf}}{\kappa}\mathbf{u} = -\frac{\partial \mathbf{p}}{\partial \mathbf{x}} + \mu\_{\rm nf}(\frac{\partial^2 \mathbf{u}}{\partial \mathbf{x}^2} + \frac{\partial^2 \mathbf{u}}{\partial \mathbf{y}^2} + \frac{\partial^2 \mathbf{u}}{\partial z^2}) \tag{6}$$

$$\frac{\mu\_{\rm nf}}{\kappa} \mathbf{v} = -\frac{\partial \mathbf{p}}{\partial \mathbf{y}} + \mu\_{\rm nf} \left( \frac{\partial^2 \mathbf{v}}{\partial \mathbf{x}^2} + \frac{\partial^2 \mathbf{v}}{\partial \mathbf{y}^2} + \frac{\partial^2 \mathbf{v}}{\partial z^2} \right) \tag{7}$$

$$\frac{\mu\_{\rm nf}}{\kappa} \mathbf{w} = -\frac{\partial \mathbf{p}}{\partial \mathbf{z}} + \mu\_{\rm nf} (\frac{\partial^2 \mathbf{w}}{\partial \mathbf{x}^2} + \frac{\partial^2 \mathbf{w}}{\partial \mathbf{y}^2} + \frac{\partial^2 \mathbf{w}}{\partial z^2}) \tag{8}$$

Energy formulation for the porous flow,

$$(\rho\_{\rm nf} \mathbf{C} \mathbf{p}\_{\rm nf})\_{\rm eff} (\mathbf{u} \frac{\partial \mathbf{T}}{\partial \mathbf{x}} + \mathbf{v} \frac{\partial \mathbf{T}}{\partial \mathbf{y}} + \mathbf{w} \frac{\partial \mathbf{T}}{\partial z}) = (\mathbf{k}\_{\rm nf})\_{\rm eff} (\mathbf{u} \frac{\partial^2 \mathbf{T}}{\partial \mathbf{x}^2} + \mathbf{v} \frac{\partial^2 \mathbf{T}}{\partial \mathbf{y}^2} + \mathbf{w} \frac{\partial^2 \mathbf{T}}{\partial z^2}) \tag{9}$$

The effective conductivity and the heat capacity combine the porous material and the flow. For further information, readers should consult the reference by Bayomy and Saghir [29]. For the purpose of studying the performance of the wavy channel, two important parameters have been investigated. The first is the local Nusselt number, and the second is the Performance Enhancement criterion. The Nusselt number is defined as the ratio of the convective heat coefficient multiplied by the inlet pipe diameter over the water conductivity (i.e., hD kw ). The heat convection coefficient is known as the ratio of the heat flux over the temperature θ which is the temperature calculated at 1 mm below the interface minus the inlet temperature Tin (i.e., θ = T − Tin). The blue dots shown in Figure 1a indicate the location where the temperature was calculated.

The performance evaluation criterion is an important parameter which combines the average Nusselt number and the friction factor. Nanofluid and hybrid fluid exhibit a larger pressure drop when compared to water. In our analysis, the performance evaluation criteria is shown in Equation (10).

$$\text{PEC} = \frac{\overline{\text{Nu}}}{\text{(f)}^{0.333}} \tag{10}$$

Here Nu is the Nusselt number defined earlier, and f is the friction coefficient in the channels.

The fanning friction coefficient is known to be represented by Equation (11).

$$\mathbf{f} = \frac{4(\Delta \mathbf{p})}{(\frac{\mathbf{L}}{\mathbf{D}})(\rho\_{\rm nf})(\mathbf{u}\_{\rm in}^2)}\tag{11}$$

The pressure difference shown in Equation (11) is the pressure taken at the middle of the inlet mixing chamber to the one at the middle of the exit chamber. Here L is the channel Length equal to 0.0375 m without waviness.

#### *Boundary Conditions and Solution Approach*

The boundary conditions used in the model consist of applying an inlet velocity uin and an inlet temperature Tin. The heat flux is applied at the bottom of the heated block, and the remaining external surface is insulated. Figure 1a shows the graphical location of the boundary condition. Porous material is used as an insert inside the channel. At the exit of the flow, a free flow boundary condition is applied. Different approaches exist in COMSOL to tackle the convergence criteria. In this particular model, the default solver used was the segregated method. Details about this approach could be found in any finite element's textbook. The convergence criterion is clearly explained in COMSOL manual. In a short summary, the convergence criteria were set as follows: at every iteration, the average relative error of u, v, w, p and T were computed. These were obtained using the following relation:

$$\mathbf{R}\_{\mathbf{C}} = \frac{1}{\mathbf{n} \cdot \mathbf{m}} \sum\_{i=1}^{\mathrm{i}=\mathrm{m}} \sum\_{\mathbf{i}=1}^{\mathrm{i}=\mathrm{m}} \left| \frac{(\mathbf{F}\_{\mathrm{i},\mathrm{j}}^{\mathrm{s}+1} - \mathbf{F}\_{\mathrm{i},\mathrm{j}}^{\mathrm{s}})}{\mathbf{F}\_{\mathrm{i},\mathrm{j}}^{\mathrm{s}+1}} \right| \tag{12}$$

where F represents one of the unknowns, viz., u, v, w, p, or T; s is the iteration number; and (i, j) represents the coordinates on the grid. Convergence is reached if R<sup>c</sup> for all the unknowns is below 1 <sup>×</sup> <sup>10</sup>−<sup>6</sup> in two successive iterations. For further information on detailed solution method, the reader is referred for COMSOL software manual [33].

#### **4. Mesh Sensitivity Analysis, Convergence Criteria and Model Verification**

The mesh sensitivity is examined in purpose of determining the optimal mesh required for the analysis. In the table below, we demonstrated different mesh sensitivity which were investigated following the terminology used by COMSOL software.

The mesh levels that COMSOL supports and the elements number for each mesh level are shown in Table 2. The average Nusselt number was evaluated at 1 mm below the interface in the aluminum block, and the results are represented in the Figure 2a. Here, the heat flux applied is equal to 75,000 W/m<sup>2</sup> , and the conductivity of the water was used to evaluate the Nusselt number. It is evident that a coarse or normal mesh level will be suitable to be used in the COMSOL model. Figure 2b, presents the finite element mesh used in our simulation with normal mesh level.


**Table 2.** Mesh information for different level of meshing [33].

The model has been tested against experimental data to demonstrate the model accuracy. Welsford, Saghir et al. [31], Plant and Saghir [24,27] and Delisle, Saghir et al. [32] conducted experimental measurement of heat enhancement in straight porous channel. They used the same proposed numerical model except the channels were straight channels. Identical boundary conditions are used as well. Results revealed a good agreement between the experimental measurement and the numerical code. Thus, the accuracy of the current numerical model.

**Figure 2.** Mesh sensitivity analysis. (**a**) Different mesh performance. (**b**) Final adopted mesh.

#### **5. Results and Discussion**

− − − − In the present study, an attempt is made to investigate the effectiveness of wavy channel in improving heat enhancement. Different flow rates were applied corresponding to 0.05 US gallon per minutes (USGPM), 0.1 USGPM, 0.15 USGPM and 0.2 USGPM. This corresponds to a flow rate of Q<sup>1</sup> = 3.15 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>m</sup>3/s, Q<sup>2</sup> = 6.3 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>m</sup>3/s <sup>Q</sup><sup>3</sup> = 9.45 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>m</sup>3/s and Q<sup>4</sup> = 1.26 <sup>×</sup> <sup>10</sup>−<sup>5</sup> <sup>m</sup>3/s, respectively. The question the authors raised is whether a wavy channel leads to a better performance enhancement criterion when compared to a straight channel? Different fluids are used in the current analysis with water, then a 1% TiO<sup>2</sup> nanoparticles diluted in 99% water, and finally, a hybrid fluid which consists of 1% nanoparticles containing 90% Al2O<sup>3</sup> and 10% copper diluted in 99% water. The differences between these fluids are their conductivity, density and viscosity. Thermal conductivity may affect the Nusselt number whereas the viscosity, specific heat and density can affect the friction coefficient and thus the pressure drop. Different flow rates will be applied, but the heating condition remains the same with a heat flux of 75,000 W/m<sup>2</sup> . The model was studied for three different permeabilities by maintaining the porosity constant at 0.9.

#### *5.1. Heat Enhancement Using Water as Working Fluid*

Figure 3 presents the temperature distribution 1 mm below the interface for the three porous material types and for four different flow rates as stated earlier. The temperature profile shows that at the beginning of the flow entrance, the boundary layer is very small, thus allowing the heat to pass from the solid block to the fluid. However, as the boundary layer starts developing, it appears that it reduces the amount of heat extracted from block. This temperature profile is identical regardless of the permeability of the material. By carefully examining the temperature magnitude, it is noticeable that as the permeability varies from 10 PPI to 40 PPI, the heat extraction is improving. This is shown as the temperature in the heated block drop in magnitude. It is evident that having porous material helps to absorb more heat. To study further the heat extraction, Figure 4 presents the local Nusselt number variation for the case or a permeability of 40 PPI.

A negative slope for the local Nusselt number variation is an indication that as the boundary layer developed, the heat extraction decreased accordingly. It is evident for this case that as the flow rate increased, the Nusselt number increased. This increase is noticeable around 3% for the highest flow rate of 0.2 USGPM. For a constant flow rate of 0.2 USGPM, one may notice that the average Nusselt number increases by 3% when compared to the case of 10 PPI and by an additional 2% when compared to 20 PPI case. This is another indication that as the permeability increases the heat extraction improves. Similar observations are noted for three remaining flow rates.

#### *5.2. Heat Enhancement Using 1% TiO<sup>2</sup> in Water as Working Fluid*

The previous model is repeated by using Titanites as the working fluid. As indicated earlier, a 1% Titanium oxide in water solution is used. The nanoparticles diameter is found to be around 31 nm. Figure 5 presents the local Nusselt number variation along the flow for the four flow rates used and for a permeability of 20 PPI.

(**c**) Permeability 40 PPI.

**Figure 3.** Temperature distribution with water as working fluid.

Since the thermal conductivity of the Titanite is higher, one expects a heat enhancement improvement when compared to water solution. By carefully examining the average Nusselt number ad for the same permeability, the 1% TiO2/water nanofluid exhibits a higher average Nusselt number when compared to water. This improvement is merely 0.5% for a flow rate of 0.2 USGPM. As the permeability increases to 40 PPI, the additional increase in average Nusselt number is found to be around 0.5% as well. Then, one may conclude that nanofluid increases the amount of heat removal.

**Figure 4.** Local Nusselt number variation for different flow rates.

**Figure 5.** Local Nusselt number variation for a permeability of 20 PPI.

#### *5.3. Heat Enhancement Using 1% (Al2O3–Cu) in Water as Working Fluid*

The model was repeated with a hybrid fluid consisting of a single nanoparticle containing Al2O<sup>3</sup> and Copper. Figure 6 displays the local Nusselt number variation for the case of a 10 PPI permeability and for the four flow rates. The Nusselt number variation profile is identical to the other cases, but its magnitude is different.

**Figure 6.** Local Nusselt number variation for a permeability of 10 PPI. (hybrid fluid is the working fluid).

A similar observation was made regarding the average Nusselt number when compared to water. It is evident that there was not a large change when compared to water. This may be due to the fact that water can provide a good cooling fluid if it is given the chance to circulate more in a hot surface. It is worth noting that using the thermal conductivity of the water in the Nusselt number calculation offers an opportunity to make a correct comparison between the three fluids.

#### *5.4. Performance Evaluation Criteria for All Fluids*

In the previous sections, we have investigated the importance of heat removal based on different flow rates and porous medium permeabilities. However, an important parameter worth investigation is the pressure drop. One may find a suitable fluid for heat removal but at the expense of higher pressure drop. To overcome this issue, it is important to combine the heat effect and the fluid effect by using the performance evaluation criteria. As shown in Equation (10), it is defined as the ratio of the average Nusselt number and the friction factor to the power one third. The friction factor is defined in Equation (11), which combines the physical properties of the fluid, the inlet velocity and the geometrical dimension of the channel. Figure 7 displays the average Nusselt number, the pressure drops and the performance evaluation criteria for all cases when the permeability is set at 40 PPI. Figure 7c displays the PEC and shows that the TiO2/water nanofluid exhibits a slightly better performance than water, and the hybrid is the worst candidate between the three studied fluids. In order to determine the reason, Figure 7a displays the average Nusselt number for the three fluids. Almost all of them have a close to identical heat enhancement which may indicate that the waviness of the channel enhances heat removal regardless the fluid used. However, in Figure 7b, a large pressure drop is found for the hybrid fluid contrary to the other remaining two fluids. Thus, there is justification for a higher PEC for the water and the TiO<sup>2</sup> nanofluid.

As the permeability changes to 20 PPI, the fluid circulates with less obstruction, thus less pressure drop is observed. Figure 8 presents the PEC for all cases at a permeability of 20PPI. It is evident from this figure that the nanofluid slightly outperformed the water and more evident that the hybrid performance is very weak. It is obvious that as the flow rates increase, the PEC increases accordingly.

Finally Figure 9 presents the cases when the permeability is set at 10 PPI. Larger pore provides less pressure drop and thus higher PEC. Nevertheless, the same observation is justified here which indicates that the TiO<sup>2</sup> nanofluid exhibits the highest performance evaluation criteria followed by the water. The hybrid fluid due to its high pressure drop exhibits the lowest PEC.

(**a**) Average Nusselt number for all cases for 40 PPI.

(**b**) Pressure drop for all cases for 40 PPI.

(**c**) PEC for all cases for 40 PPI.

**Figure 7.** Average Nusselt number, pressure drop and PEC for all cases at a permeability of 40 PPI.

**Figure 8.** PEC for all cases when the permeability is set at 20 PPI.

**Figure 9.** PEC for all cases and for a permeability of 10 PPI.

#### **6. Conclusions**

TiO<sup>ଶ</sup> AlଶO<sup>ଷ</sup> TiO<sup>ଶ</sup> AlଶO<sup>ଷ</sup> This paper presented a numerical study of the heat performance of three different fluids mainly a single fluid water, a nanofluid 1% TiO2/water and a hybrid fluid 1% (Al2O3–Cu)/water. Three wavy channels with a porous insert, having the three different permeabilities of 10 PPI, 20 PPI and 40 PPI, were investigated numerically. For each case, four different flow rates were implemented. Results reveal the following:


**Author Contributions:** Conceptualization, K.M.E. and M.Z.S.; methodology, K.M.E.; software, K.M.E.; validation, K.M.E. and M.Z.S.; formal analysis, K.M.E.; investigation, K.M.E.; resources, M.Z.S.; data curation, K.M.E.; writing—original draft preparation, K.M.E.; writing—review and editing, M.Z.S.; visualization, K.M.E.; supervision, M.Z.S.; project administration, M.Z.S.; funding acquisition, M.Z.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Qatar Foundation, grant number NPRP12S-0123-190011.

**Institutional Review Board Statement:** Not Applicable.

**Informed Consent Statement:** Not Applicable.

**Data Availability Statement:** Not Applicable.

**Acknowledgments:** This research was funded by the National Science and Engineering Research Council Canada (NSERC), the Faculty of Engineering and Architectural Science at Ryerson University, and the Qatar Foundation, grant number NPRP12S-0123-190011.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **References**


**Yao Zhao, Kai Zhang \* , Fengbei Guo and Mingyue Yang**

College of Metrology & Measurement Engineering, China Jiliang University, Hangzhou 310018, China; s1802080449@cjlu.edu.cn (Y.Z.); s20020804019@cjlu.edu.cn (F.G.); s1902080445@cjlu.edu.cn (M.Y.) **\*** Correspondence: zkzb3026@cjlu.edu.cn; Tel.: +86-137-57-135-320

**Abstract:** 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.

**Keywords:** microfluidic; flow distributions; fluid network

**Citation:** Zhao, Y.; Zhang, K.; Guo, F.; Yang, M. Dynamic Modeling and Flow Distribution of Complex Micron Scale Pipe Network. *Micromachines* **2021**, *12*, 763. https://doi.org/ 10.3390/mi12070763

Academic Editors: Junfeng Zhang and Ruijin Wang

Received: 28 April 2021 Accepted: 22 June 2021 Published: 28 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. 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

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.

#### **2. 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:


The construction of the microfluidic pipe network model mainly includes the following two parts:


The directed topology model of the microfluidic network is shown in Figure 1. In the microfluidic network, both the micropump, as a power source, and the network outlet

are regarded as boundary nodes. Moreover, the intermediate branches include various resistance components, such as a micropipe and microvalve.

**Figure 1.** Directed topology model of the microfluidic network.

For a stable microfluidic network, the continuity equation (continuous mass) in the flow process can be expressed as follows:

$$V\_i \frac{d\rho\_i}{dt} = \sum\_{j=1}^{N} D\_{ij} G\_{ij} \tag{1}$$

Without considering the influence of heat from *dρ dt* = *∂ρ ∂p dp dt* + *∂ρ ∂H dH dt* ≈ *∂ρ ∂p dp dt* , the following equation can be obtained:

$$\mathbf{C}\_{i}\frac{dp\_{i}}{dt} = \sum\_{j=1}^{N} D\_{ij}\mathbf{G}\_{ij} \tag{2}$$

where, *C<sup>i</sup>* = *V<sup>i</sup> ∂ρ<sup>i</sup> ∂p<sup>i</sup>* 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; *Gij* is the mass flow between nodes *i* and *j*; *V<sup>i</sup>* is the volume of node *i*; *Dij* is the connection mode between nodes *i* and *j* (*i* = 1, 2, . . . , *N*; *j* = 1, 2, . . . , *N*) and the specific meaning is:

*Dij* 1, 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*


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.

**Figure 2.** Schematic diagram of the topological structure of the *i*-th node of the pipeline. *P* is the pressure at node *i*. *Gij* is the flow rate between node *i* and node *j*.

The momentum conservation equation in the microfluidic network can be expressed as follows:

$$\rho\_{i\bar{j}}L\_{i\bar{j}}\frac{d\mathcal{U}\_{i\bar{j}}}{dt} = D\_{i\bar{j}} \times \left(P\_{\bar{j}} - P\_{\bar{i}} + D\_{i\bar{j}} \times H\_{i\bar{j}}\right) - h\_w \tag{3}$$

where, *Hij* is the pressure generated by macro kinetic energy, potential energy and power source between node *i* and *j*; *Uij* is the fluid velocity between node *i* and *j*; *Dij* has the same meaning as *Dij* in Equations (1) and (2).

In the formula, *h<sup>w</sup>* = ∑ h 0.5*λ L <sup>d</sup>* + *ξ* i 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 *Iij* = *Lij Aij* , the friction resistance coefficient on the pipeline *Rij* = *hw* (*Aijuijρij*) 2 , and no power supply in the pipeline, *Hij* = 0 is introduced into the Equation (3):

$$I\_{ij}\frac{dG\_{i\bar{j}}}{dt} = D\_{i\bar{j}} \times \left(P\_{\bar{j}} - P\_{\bar{i}}\right) - R\_{i\bar{j}} \times G\_{i\bar{j}}^2 \tag{4}$$

*Rij* is the resistance characteristic term between node *i* and *j*, which can be expanded as follows:

$$R\_f = \left(\frac{\lambda L}{d} + \tilde{\xi}\right) \frac{1}{2\rho A^2} \tag{5}$$

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:

$$C\_i^{r+1} \frac{P\_i^{r+1} - p\_i^r}{\Delta t} = \sum\_{j=1}^N D\_{ij}^r \sqrt{\frac{D\_{ij}^r \left(p\_j^{r+1} - p\_i^{r+1}\right)}{R\_{ij}}} \tag{6}$$

where, *r* is the discrete-time variable, time is *t* = *t*<sup>0</sup> + *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.

**Figure 3.** Schematic diagram of calculation process using fluid network method.

#### **3. Results and Discussions**

*3.1. 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.** Structure of the microfluidic platform.

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, *RF*<sup>1</sup> is the microchannel resistance of the microvalve output, *RF*<sup>2</sup> is the main channel resistance of fluid flow, and *RF*3, *RF*4, *RF*<sup>5</sup> 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]:

$$R\_f = \frac{12uL}{(1 - 0.63(h/w))h^3w} \tag{7}$$

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.

3 4 5

=

−

1

2

**Figure 5.** Fluidic resistors of the microfluidic circuit.

**Figure 6.** Electrical model of microfluidic platform.

*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:

$$P\_0(t) - P\_1(t) = Q\_1(t)R\_{F1} \tag{8}$$

$$P\_1(t) - P\_2(t) = Q\_2(t)R\_{F2} \tag{9}$$

$$P\_1(t) - P\_4(t) = Q\_4(t)R\_{F4} \tag{10}$$

$$P\_2(t) - P\_3(t) = Q\_3(t)R\_{F3} \tag{11}$$

$$P\_2(t) - P\_5(t) = Q\_5(t)R\_{F5} \tag{12}$$

$$Q\_1(t) = Q\_3(t) + Q\_4(t) + Q\_5(t) \tag{13}$$

$$Q\_2(t) = Q\_3(t) + Q\_5(t) \tag{14}$$

$$P\_0(t) = \frac{P\_{\bar{l}}V\_{\bar{l}}}{(V(t) + V\_{\bar{l}})} \tag{15}$$

$$V(t) = \int q(t)dt\tag{16}$$

#### *3.2. 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.

**Table 1.** Dimensions of the device and values of some parameters of the setup.


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*<sup>2</sup> in the main microchannel is the sum of *Q*<sup>3</sup> and *Q*5. The flow rate of *Q*<sup>1</sup> is the sum of *Q*<sup>2</sup> 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.


**Table 2.** Model Simulation and Simulink calculation results under different pressures.

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.

**Figure 7.** Schematic diagram of different driving pressures.

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.


**Table 3.** Simulation results under different lengths of microchannel 4.

**Figure 8.** Schematic diagram of calculation results under different lengths of microchannel 4.

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, *Rij*, 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.

#### **4. 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.

**Author Contributions:** Conceptualization, Y.Z.; methodology, Y.Z. and K.Z.; validation, Y.Z., K.Z., M.Y. formal analysis, K.Z.; writing—original draft preparation, Y.Z. and F.G.; writing—review and editing, K.Z. and M.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant No. 11472260).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

