One-dimensional numerical calculation for unsteady gas-liquid flow in micro channels with air gaps

A numerical calculation model for gas-liquid unsteady two-phase flow in a micro channel was established and validated. The model focuses on flow in a channel including air gaps. It also reduces calculation cost by minimizing the number of control-volume elements compared with conventional numerical methods for gasliquid two-phase flow. Gas-liquid two-phase flow in channels with diameters of 1 to 2 mm is important in liquid transportation in chemical processing and analysis. As for designing flow channels in chemical processing and analysis, simple numerical models are desirable to evaluate many possible patterns quickly. The model uses moving boundaries that correspond to gas-liquid interfaces, but it represents boundaries between different channels with fixed boundaries. By setting several different configurations of air gaps, the numerical model was validated in regard to position and volume of air gaps in the channel. As an application of the numerical model to the design of flow channels in chemical processing and analysis, the necessary movement conditions of a syringe pump to achieve quick liquid transportation were investigated. By applying the numerical model, the necessary conditions to minimize flow rate oscillation were determined. In simulation-based design of microchannel, the numerical model in this research is an effective tool to determine design parameters quickly.


Introduction
In applications such as chemical processing and analysis, tubes with diameters of 1 to 2 mm are used to handle liquids with volumes of 1-100 microliters. In order to prevent different kinds of liquids mixed, air gaps are often included in liquid transportation during chemical processing and analysis. In this case, air gaps flow as slugs and occupy most of the cross-sectional area of the pipe (Triplett et al., 1999;Chen et al., 2002).
To achieve high reliability and throughput chemical processing and analysis, quick liquid transportation is required. Pressure oscillation in fluid transportation was numerically studied in Nagaoka et al. (2010), and the effects of diameter and viscosity of fluids were investigated. Pressure and flow rate oscillation should be reduced in liquid transportation to achieve high reliability and throughput.
Unsteady flows in pipes have been studied in the context of water hammer caused by opening and/or closing of valves in flow channels. As for unsteady gas-liquid two phase flows, several numerical methods were proposed. For example, a model using distribution of void fraction in pipes was proposed in Ito et al. (2004) and Ozawa et al. (2012). This model requires a large number of elements to resolve gas-liquid interfaces. Models that consider air gaps and liquids as springs connecting concentrated masses were proposed in Matsui et al. (1979), Akagawa et al. (1982), and Ishikawa et al. (2014). These models, however, are not suitable for pipelines with localized air gaps. Models that track gas slugs directly were proposed in Al-Safran et al. (2004) and Rosa et al. (2015). These slug-tracking methods do not consider the unsteady flows with slugs that occupy most of cross-sectional area of pipes. In applications to simulation-based design (Greco et al., 2019, as an example) of flow channels in chemical processing and analysis, flow in pipes should be calculated with a small number of elements to reduce calculation cost and confirm design parameters.
In above research, numerical calculations require a large number of elements to calculate unsteady flows containing air gaps. Especially with small volume of air gaps, volume of the elements must be reduced, which leads to large number of elements. The purpose of this research is to establish a numerical model with a small number of control-volume elements for unsteady flow in channels with diameters of 1-2 mm including air gaps. Volume and position of the air gaps are important factors in chemical processing and analysis. Thus, the numerical models must be able to track the position of gas-liquid interfaces. In this study, a low-cost numerical model for designing flow channels in chemical processing and analysis was established and validated by experiments. In addition, as an application of the established numerical model, the movement of a syringe pump was designed by a simulation based on the model.

Numerical model 2.1 Basic equations
A numerical model based on one-dimensional flow in a pipe is established as follows. In consideration of onedimensional flow, conservation of mass and momentum is described as follows. (1) Here, friction loss in pipes, expansion, and contraction at the connection between pipes, and loss of potential energy due to the gravity are included in pressure loss, Δp. For frictional pressure loss in pipes, both laminar and turbulent flows are considered. In order to consider compressibility of fluids, the following barotropic equation is used to describe the change in volume of fluids in pipes: where δ represents deviation in physical properties. Radial deformation of a pipe is given by the following equation: (1)-(4), the basic equations for one-dimensional flow in a pipe are obtained as follows.

Control-volume finite-element method using pipe and fluid elements
To calculate unsteady flow in pipes, Eqs. (5) and (6) are solved with the control-volume finite-element method (CVFEM) using pipe and fluid elements. As an example, a computational model for multiple fluids in multiple pipes is shown in Fig. 1. In this example, five types of pipes, including three kinds of fluids, are connected. To reduce the number of elements, boundaries between elements were set only between different pipes and between different fluids. Moving boundaries that correspond to interfaces between different fluids were used to track movement of the fluids. In the figure, the boundary between elements 2 and 3 follows the movement of fluids A and B, and the boundary between elements 4 and 5 follows the movement of fluids B and C. On the contrary, the boundaries between the different pipes are fixed. As interfaces between different fluids are expressed as moving boundaries, this method can only be applied to slug flow and the interfaces cannot be broken or split during calculation.
where suffix i represents each element, and suffix n represents a time step. In Eqs. (7) and (8), velocity, pressure, and cross-sectional area are set to the upstream and downstream end of each element. On the contrary, one set of parameters is allocated to each element for another fluid and pipe parameters. Parameters used in the calculations are listed in Table  1. Here, velocity and pressure at the same boundary (e.g., 1,D and 2,U in Table 1) are set to keep continuity in the numerical calculation. As for velocity, 1,D and 2,U are set to keep continuity of mass flow rate at the boundary of elements 1 and 2.

Experimental system and input conditions 3.1 Experimental system
A liquid-transportation system used in chemical processing and analysis is focused on hereafter. In such a system, air gaps are often included in the flow channels, and syringe pumps are equipped to control the flow rate and volume. The experimental system used in this study is shown schematically in Fig. 2. This system mimics liquid transportation from a water cup. The flow channel between the syringe pump and water cup is isolated by a solenoid valve during measurements of pressure. The water source fills the flow channel with water and does not affect the pressure measurements.
The flow channel between the syringe pump and the water cup is composed of a nozzle, two tubes, a pressure sensor, and a syringe pump. The syringe pump induced unsteady flows in the flow channel. The nozzle is dipped in the water cup, and water is aspirated into the nozzle. A pressure sensor, PGMC-A-200KP (Kyowa), recorded pressure histories in the channel. The pressure sensor was attached to the cylinder block with an inner diameter of 8.0 mm and the cylinder block was connected between tube 1 and tube 2 directly. The pressure data were recorded on an EDX-200A (Kyowa) data logger. Sampling rate of pressure data was set to 1 kHz.
Sizes and physical properties of each component between the syringe pump and water cup are listed in Table 2. The values listed in Table 2 were also used in the numerical calculations described in Section 2.

Initial conditions of experiments
To mimic air gaps in a liquid-transportation system used in chemical processing and analysis, air gaps were artificially induced into the nozzle as initial conditions. By changing the volume and position of the air gaps, the validity of the numerical calculation for unsteady flows with air gaps can be checked. The nozzle was divided into three sections, as shown schematically in Fig. 3, and the fluids included in each section were changed. Each section, named "sections 1, 2, and 3" from top to bottom, has a length of 30 mm. The bottom of section 3 corresponds to the boundary between the nozzle and the water in the cup. As shown in Table 3, seven conditions were set as the initial configurations of fluids in the nozzle. Either water or air was allocated to each section. b represents the position of gas-liquid boundary closest to the top of the nozzle. The position was measured from the bottom of the nozzle. w represents the total length of water in sections 1, 2, and 3. For example, under condition 3, an air gap is included in section 2.
Condition 3 is shown schematically in Fig. 4. Two gas-liquid boundaries exist in the nozzle: one located between sections 1 and 2, and the other located between sections 2 and 3. The boundary between sections 1 and 2, which is 60.0 mm from the bottom of the nozzle, is closest to the top of the nozzle. That is, b = 60.0 mm. Sections 1 and 3 include water. That is, w = 30.0 × 2 = 60.0 mm.
Under the initial condition, conditions 2, 3, and 4 have different values of b , so they can show the effect of the position of the air gap on the pressure histories. Conditions 4, 5, 6, and 7 have different values of w , so they can show the effect of the total volume of air gaps on the pressure histories.

Input condition of unsteady flow
In the same manner as liquid transportation used in chemical processing and analysis, unsteady flow in the flow channel was generated by the syringe pump. As shown in Fig. 5, the syringe changed the flow rate from 0 to 200 μL/s within 10 ms, then unsteady flow was occurred. In this case, water hammer in the flow channel is similar to slow valve opening (Kodura 2016;Yue et al. 2019). The syringe pump was set to same flow rate setting for all conditions listed in Table 3.
In syringe aspiration process, pressure in the channel first dropped and then oscillated in a damped fashion. The pressure history of the oscillation process is expected to depend on parameters b and w . To validate the numerical model, both experimental and numerical dependencies of b and w on pressure were evaluated.

Results and discussion 4.1 Validation of numerical model
To check the validity of the numerical model described in Section 2, the experimentally measured and numerically calculated pressure histories were compared. First, to overview of pressure histories, the results obtained under condition 1 (Table 3) were compared. Second, to check the effects of air-gap position, the results obtained under conditions 2, 3, and 4 (Table 3) were compared. Third, to check the effects of the volume of the air gap, the results obtained under conditions 4, 5, 6, and 7 (Table 3) were compared.
The pressure histories obtained by experiment and numerical calculations under condition 1, which is the simplest case because no air gaps were included in the flow channel, are shown in Fig. 6. The pressure histories were calculated with numerical model of this research and Amesim hydraulic model (see Giuffrida 2002 andChen 2018 as examples of Amesim use cases). As start of syringe movement was not measured in the experiment, we set time origin (t = 0 s) to the point that pressure drops below -3 kPa. Due to pressure propagation time in the flow channel, syringe actually started moving 8 ms prior to time origin (t = 0 s) determined by pressure drop. Therefore, pressure change before time origin (t = 0 s) were observed. Pressure oscillation caused by water hammer is observed between 0-0.2 s. The oscillation during syringe aspiration is gradually damped and converges to a stable pressure, which is determined by frictional pressure loss in the channel. Figure 6 shows that the numerical calculations can reproduce the profile of pressure oscillations caused by unsteady flow in the micro channel.
Pressure amplitudes of experiment, calculations in this research and by Amesim are 48.9, 33.6, and 37.8 kPa, respectively. Additionally, pressure oscillation periods of those are 49, 50, and 54 ms, respectively. Comparing the experiment and the calculations in this research and by Amesim, pressure amplitude of the experimental were larger than that of the calculations. Because mechanical vibration of syringe is not considered in the numerical calculation, pressure amplitude was underestimated. Numerical methods should be carefully used when considering a cavitation due to deep pressure drop. Differences of pressure amplitude and period between calculations in this research and by Amesim were less than 10%. Considering that Amesim uses some control-volume elements in one pipe and the numerical calculation in this research uses one control-volume element in one pipe, smaller numerical cost can be achieved in this research. Note that Amesim cannot calculate the slug flow in pipes such as liquid and air gap. To verify the effect of air-gap position, the numerically calculated and experimentally measured pressure histories obtained under conditions 2, 3, and 4 were compared. Under conditions 2, 3, and 4, the volume of air gap included in the flow path is constant, but the position of the air gap is changed; that is, w is constant, but b differs.
As for water hammer in two-phase slug flow, the length of the liquid slug closest to the open/close valve has major effect on the pressure surge of water hammer (Akagawa et al., 1982). In this study, the syringe causes water hammer in the channel, and b is related to the length of the liquid slug connected to the syringe. Pressure histories under conditions 2, 3, and 4, are shown in Fig. 7: (a) and (b) show numerical and experimental results, respectively. As previously reported (Akagawa et al., 1982), pressure amplitude decreases when the length of the liquid slug closest to the syringe decreases. That is, the smallest pressure amplitude is observed under condition 4, which has largest b . Compared with condition 2, condition 4 showed 14.7 % smaller pressure amplitude in the numerical calculation and 16.7 % smaller pressure amplitude in the experiment. Since the effect of air gaps on pressure history was reproduced by the numerical calculation, it is concluded that the numerical model is valid in regard to the position of air gaps.
(a) Numerical calculations (b) Experimental measurements Fig. 7 Pressure histories during syringe aspiration. The position of the air gap is changed as shown in Table 3 ( To verify the effect of air-gap volume, the numerically calculated and experimentally measured pressure histories obtained under conditions 4, 5, 6, and 7 were compared. Under conditions 4, 5, 6, and 7, the length of the liquid slug connected to the syringe is constant, but the volume of the air gap differs; that is, b is constant, but w differs. Pressure histories obtained under conditions 4, 5, 6, and 7 are shown in Fig. 8: (a) and (b) show numerical and experimental results, respectively. In all cases, b is constant; that is, the length of liquid slug closest to the syringe is constant. It is clear that as the volume of the air gap decreases, large frictional-pressure loss in the flow path leads to low pressures around 0.03 s. Therefore, the largest pressure amplitude is observed under condition 7 (i.e., smallest w ), and the smallest pressure amplitude is observed under condition 4 (i.e., largest w ). Compared with condition 7, condition 4 showed 17.2 % smaller pressure amplitude in the numerical calculation and 16.7 % smaller pressure amplitude in the experiment. In consideration that the effect of air volume on pressure histories was reproduced by the numerical calculation, the numerical model was validated in regard to volume of air gaps. The above results validate the numerical calculation; that is, the position and volume of the air gaps can be appropriately modeled by the numerical calculation.

Discussion on quick liquid transportation
In a chemical processing and analysis, a set volume of fluid should be transported. To enhance throughput of the system, the fluid transportation must finish in a short time. Accordingly, the numerical model described in Section 2 was used to determine the necessary conditions for shortening the fluid-transportation time.
To investigate the necessary conditions, condition 2 (Table 3) was set as the initial condition, and 10 μL water was aspirated into the nozzle. Flow rate setting of the syringe pump used in the numerical model is shown in Fig. 9. Maximum flow rate was set to Max , and duration of acceleration and deceleration of the syringe pump was set to 10 ms. Max was set to 100-500 μL/s, and total duration of aspiration was adjusted to keep total aspirated volume to 10 μL. As Max increases, duration of the syringe movement decreases. To guarantee a set volume of liquid transportation, oscillation of flow rate should be suppressed until the end of the transportation. Accordingly, the time needed to stabilize flow rate at the tip of the nozzle was evaluated. Stabilized time st is defined as the time needed to stabilize the flow rate. Calculated flow rate at the tip of the nozzle is plotted in Fig.  10, in which st is defined as the minimum time that the flow rate is between -10 and 10 μL/s after st . In the measurement of st , the time origin is set to the beginning of syringe movement. The relation between aspiration flow rate Max and stabilized time st calculated by the numerical model is shown in Fig. 11. It is clear that as Max increases, duration of the syringe movement decreases. However, oscillation in flow rate increases at the same time. This leads to larger stabilized time st at Max = 500 μL/s compared with that at smaller Max . Moreover, st reaches a minimum at Max = 220 μL/s; that is, at Max = 220 μL/s, liquid transportation can be achieved in a short time. Since the relation between st and Max shows periodical local minimum and maximum, it is important to control oscillation in flow rate to achieve quick liquid transportation. Flow rate histories calculated by the numerical model with Max = 500, 220, and 140 μL/s are shown in Fig. 12. According to Fig. 11, at Max = 500 μL/s, duration of syringe movement is shortest. At Max = 220 μL/s, st is shortest, and at Max = 140 μL/s, st reaches a local maximum. The solid lines show flow rate with aspiration volume of 10 μL. As a reference, the dotted lines show flow rate when the syringe does not decelerate and maintains Max . Black lines show flow rate at the tip of the nozzle, blue lines show input condition of the syringe.
In the solid lines, oscillations caused by acceleration and deceleration are superposed after syringe aspiration finishes. In the dotted lines, only oscillations caused by acceleration are observed. Comparing these two lines reveals that the oscillations caused by acceleration and deceleration have the same phase at Max = 500 and 140 μL/s. That To check the superposition of oscillations experimentally, pressure in the flow path at Max = 500, 220, and 140 μL/s was measured under the same initial condition and flow rate setting as used for the numerical model.
Pressure histories at Max = 500, 220, and 140 μL/s are shown in Fig. 13. Start and finish of syringe movements are also shown in the figure. The same as the numerical results, the largest oscillation in pressure was observed at Max = 500 μL/s. As a result of superposition of the oscillation, the pressure peak after the finish of syringe movement reaches a minimum at Max = 220 μL/s. The numerical results shown in Fig. 11 and Fig. 12 indicate that to achieve quick fluid transportation, it is important to control the oscillations caused by acceleration and deceleration. Under the necessary conditions to shorten the fluidtransportation time, the oscillations caused by acceleration and deceleration have opposite phase, and amplitude of oscillation is reduced. The numerical model can easily estimate the necessary flow rate in the target flow path. Fig. 13 Experimentally measured pressure histories for aspiration volume of 10 μL. The dotted, dashed, and solid lines respectively correspond to Max = 500, 220, and 140 μL/s.

Concluding remarks
A model for numerically calculating gas-liquid two-phase flow in micro channels was established and validated by experiments. The numerical model includes air gaps in the channel and reduces calculation cost by minimizing the number of control-volume elements. The model tracks gas-liquid interfaces to take position and volume of the air gaps into account. In the experiments, a syringe pump equipped in the channel generated unsteady flow, and pressure histories in the channel were recorded with a pressure sensor. Several patterns of air gaps were set in the channel, and the numerical model was validated by comparing the numerically calculated and experimentally measured pressure histories. The experimentally measured effects of position and volume of the air gaps on the pressure histories were reproduced by the numerical model. As an application of the numerical model to design of flow channels in chemical processing and analysis, the necessary conditions for shortening the fluid-transportation time were determined by the numerical model. In order to shorten fluid-transportation time, the phase of oscillations caused by acceleration and deceleration of the syringe pump should have opposite phase. With the developed numerical model, the necessary flow rate for the target flow path can be easily estimated. By using the numerical method of this research, design parameters of flow channels in chemical processing and analysis are determined quickly. We believe that our numerical method is an effective tool in simulation-based design of micro channels. However, our method can only be applied for two-phase slug flow without mixing. In the future, we will improve our numerical method to calculate various flow patterns in micro channels.