Experimental and numerical investigation on conjugate performance of fan and heat exchanger of helicopter oil cooling system

,


Introduction
The combination of axial fan and high-efficiency heat exchanger (HE) is normally applied in the oil cooling system of helicopter. As the core components of the system, the fan and HE with reliable coupling operation performance are crucial to ensuring the flight safety and service life of helicopter. In engineering applications, these factors such as the possible change of the exhaust duct structure of oil cooling system, the effect of helicopter downwash flow [1], the fouling of HE [2,3] and the deviation of the system design can lead to the system operating point deviating from its design point. This deviation may occur surge or stall phenomenon in the axial fan used in the cooling system which could damage the fan and system. It is noted that conjugate performance of fan and HE under different flow rates can be helpful in predicting the performance effect accurately under various system operating conditions. As a consequence, the surge or stall can be avoided and even the system can continue to be used under a small flow rate if it can still meet the heat dissipation requirement.
For the axial fan, a large amount of research on fan design method [4][5]6,7,8] and performance analysis [9][10]11,12] have been conducted. Many significant achievements were obtained over the past years. In addition, with the development of computational fluid dynamics (CFD), the air flow field in axial fan can be accurately calculated. For the rotational motion of the airflow and the rotor, Multiple Reference Frame (MRF) had been proved to be an effective method [13][14][15][16]. For axial fan, the variable air inlet conditions [17][18]19] and ambient pressure [20][21]22] of the axial fan will have a significant impact on its performance. Nowadays, the researches on the axial fans mainly focus on the rotor or the rotor with behind guide vane. For the axial fan with inlet guide vane (IGV), due to its advantages of relatively small size, light weight, and small efficiency change with different working conditions, it is suitable for the helicopter oil cooling system which requires small weight and space. However, the related research on axial fan with IGV needs to be further performed. Shaw et al. [23] studied the effect of variable IGVs on the stability boundary of a transonic fan with distortion inlet airflow. Their results showed that the variable IGVs can improve the performance and stability of fan. The static characteristics of a small axial fan with IGV were investigated for different outlet angles by Liu et al. [24]. It was found that the static characteristics of the fan could be significantly improved after optimizing the IGVs. In addition, the effect of blade tip grooving on fan performance with or without IGVs was studied by Kharati-Koopaee et al. [25]. They stated that the performance parameters of the fan with IGVs were less affected by blade tip grooving changes than that without IGVs.
When the air flow rate was smaller than a certain value, the fan could suffer stall or surge and even came into hazards. Thus, there are many studies published on stall and surge of axial fan and compressors. Greitzer firstly established an original nonlinear model of the surge and rotating stall of axial compressor system and the well-known B-parameter was proposed [26,27]. Afterwards, the Moore-Greitzer nonlinear model was developed by Moore and Greitzer which was widely used to analyze the stall and surge [28,29]. The hysteresis region of compressor performance curve can be explained by their research findings. Zhang et al. [30,31] numerically studied the aerodynamic noise change and the entropy increase after the axial fan stalled. Their results showed that the noise source of the stalling fan mainly comes from the vortex in rotor blades region, and the entropy in the flow field increased after fan stalled. In addition, Zhang et al. [32,33] and Li et al. [34] studied the influence of an abnormal blade on the occurrence and development of rotating stall in axial fan. The abnormal blade would lead to the reduction of the fan performance and stability, that is, the fan would fall into unstable state in advance. Canepa et al. [35] experimentally investigated the unsteady leakage flow in low-speed axial fan by using Laser Doppler Velocimeter (LDV) and particle image velocimetry (PIV). They pointed out that leakage flow noise could be reduced by decreasing its periodicity. Moreover, Mo et al. [36] numerically studied the unsteady flow characteristic of an automotive cooling fan with Large Eddy Simulation (LES) method. Their results showed that there was strong tonal noise near the blades leading edge.
In the aspect of high efficient plate-fin HE, the studies on the design method and performance prediction of HE were carried out sufficiently [37][38]39,40,41,42,43]. In addition, Al-Zahrani et al. [44] numerically investigated the channel flow resistance and pressure drop of plate HEs. The friction factor correlation of HEs were obtained in their study. Baranyuk et al. [45] also simulated the aerodynamic resistance of HEs with different types HE surfaces. The impact mechanism of geometrical parameters on the aerodynamic resistance was explained. Moreover, Hosseini et al. [46] numerically analyzed the effect of particle deposition on the performance of a HE by simplifying the HE model with periodic boundary. Khoshvaght-Aliabadi et al. [47] studied the heat transfer and pressure resistance characteristic of plate-fin HE with different offset-strip channels by using the Al 2 O 3 -water nanofluid. They found that the channel height played significant role in thermal hydraulic performance of HE. Arsenyeva et al. [48] analyzed the effect of water fouling on the performance of plate-fin HE. They found that the fouling mitigation could be performed through optimizing the plate corrugation geometry. For a complex fin structure of HEs, the simplified methods using periodic boundary or porous medium model were normally used to analyze the heat transfer performance [49][50]51,52,53,54,55]. In addition, HE coupling fan was an important and difficult issue for the oil cooling system. In terms of coupling systems, Gao et al. [56,57] designed an oil cooling system and analyzed the flow resistance, which was mainly based on theoretical design and simplified numerical model. Liu et al. [58] conducted numerical research on a heat dissipation system of vehicle. The HE was simplified into a three-dimensional porous medium zone. The average error of heat medium temperature between the simulation and experimental results was about 2.3 %, which showed the effectiveness of their simulation method. Yu et al. [59] numerically studied the three-dimensional indoor distribution of flow field and heat transfer mechanism of fan and HEs. And Mesalhy et al. [60] simulated the coupling heat transfer performance of an axial fan with different HEs. The fan was simplified by fan plane model in Fluent in simulation. Their results showed that the optimal fin number of HEs was related to fan rotating speed and ambient pressure. Alam et al. [61] introduced a method to utilize the heat of proton exchange membrane fuel cell. The heat loss of the fuel cell was discharged into the ambient by using the axial cooling fans. Somwanshi et al. [62] developed a mathematical model for a desert cooler with new design. They also used the cooling fan as the heat discharged source. Moreover, Mesalhy [63] studied a dual-fan system used to cooling the electromechanical actuators. Their results indicated that there would be a better system performance when only operating the upstream fan, while the system performance would significantly had a obvious decrease when only operating the downstream fan due to the flow blockage of upstream fan.
To the best of the authors' knowledge, the related research on the operating performance of the helicopter oil cooling system coupled axial fan with IGV and plat-fin HE is still rare. The main contribution of the present study is to experimentally and numerically investigate the operation performance of an oil cooling system in the hysteresis region of axial fan with IGV. In addition, the physical mechanism that the fan surge and stall does not occur at small flow rate is further revealed. The objective of the current work is to present the total pressure performance of the axial fan with IGV under both two flow rate adjusting modes. The sample pieces of HEs with different fins and air side flow lengths are designed and manufactured. Moreover, the operation performance of the conjugate system including axial fan with IGV and HE was studied both experimentally and numerically in a systematic manner.

Experimental setup
Based on GB1236-2017 industrial fan-performance testing using standardized airways [64], the schematic diagram of the experimental system was shown in Fig. 1. Fig. 1 (a) showed the composition of the experimental system. The axial fan was driven by a motor with a rotating speed of 2882 rpm, and its power was measured by power The empirical constant of k-ε model meters. At the distance of 50 mm from the HE inlet and outlet, two NS-2 pressure sensors with 0.25F.S % accuracy (TM-sensor Company) were allocated at the centre of two adjacent surfaces in the transition air duct. The flow rate and dynamic pressure in the air duct behind the HE were measured by Annubar flowmeter (ANB-45) with an accuracy of ± 1 %, which was located at the point of 600 mm behind the rectifier grille. The static pressure was measured by U-shape water column manometer with an accuracy of ± 0.5 mm water column, which located at the point of 125 mm in front of the rectifier grille. By adjusting the conical throttle valve behind the air duct, the performance curves of fan and the system including fan and HE under different flow rates were obtained. Fig. 1 (b) showed the experimental test section which was composed of one HE and two transition ducts. Since the axial shape of HE was square, two identical circular-square transition ducts were manufactured to connect the HE. The square side and circular side sizes of the transition duct were consistent with the sizes of HE and experimental duct, respectively. The axial length of transition ducts was 303 mm. And two pressure measuring holes were designed at the distance of 50 mm from the square side in the transition duct, where the static pressure at the HE inlet and outlet was measured.
The axial fan with IGV used in helicopter oil cooling system is shown in Fig. 2. The main design parameters of the axial fan are listed in Table 1.
When designing the plate-fin HE, two fin structures of rectangular and offset strip were used. The flow lengths of 30 mm, 60 mm and 90 mm at the air side was determined. A total of six types of plate-fin HE were manufactured. The core size and sample piece of the HE were shown in Fig. 3, and the design parameters and nomenclature of the HEs were listed in Tables 2 and 3, respectively.
In order to determine the flow resistance of the HE, the flow performance of the system without HE (only two transition ducts and the fan) was tested as a reference. Due to the flow resistance of transition duct was small, the performance of the system without HE was considered as the fan performance in the current study. Afterwards, the system performance with axial fan and HE was tested and compared. The total pressure difference between the two systems with and without HE at the same flow rate was the HE flow resistance. When measuring the fan performance, the flow rate was adjusted with two adjusting modes as follows: (1) Adjusting mode1: the flow rate was gradually adjusted from large to small; (2) Adjusting mode2: the flow rate was gradually adjusted from small to large.
In the present study, the accuracy of thermometer was ± 0.3 ℃. The lowest temperature was about 20 ℃, thus the maximum uncertainty of temperature was 1.5 %. With an accuracy of ± 0.5 mm water column, the minimum difference of water column was 7.0 mm and the maximum uncertainty of static pressure was 7.1 % correspondingly. In the same way, the uncertainty of air flow rate and total pressure were 0.6 % and 7.1 %, respectively.

Computational details
To study system performance in detail with CFD method, the flow fields of the system including axial fan and HE-OSF-30 at different air flow rates were simulated. The mechanism of the system did not surge or stall under small flow rate would be explained by simulation results.

Governing equation
When describing the flow of airflow, the general form of the continuity, momentum can be written as follows: whereρ, ϕ and Γ ϕ are represents density, generalized variable and diffusion coefficient, respectively. S ϕ is a source term.
The steady simulation method was used when the air volume flow rate was greater than that at the stall point. Due to the system may be unstable when the air volume flow rate was less than the stall point, the transient simulation method was used when the system operated in this condition.
For both two methods, the airflow turbulence was described by the standard k-ε model:   Fig. 3. The core size and sample piece of the HE.
where, G k and G b are the turbulent kinetic energy caused by the average velocity gradient and lift force, Y M represents the contribution of pulsation expansion in compressible turbulence to the overall turbulent dissipation rate, S k and S ε are the source terms,C 1ε , C 2ε and C 3ε are the empirical constant with the value of 1.44 and 1.92, respectively. σ k and σ ε are the Prandtl number corresponding to k and ε, with values of 1.0 and 1.3, respectively.

Mesh generation and boundary conditions
The three-dimensional model of the system including fan and HE-OSF-30 was established, and the mesh of computational domain was generated by using commercial software ANSYS Fluent 2019R3. The mesh and boundary conditions of the system are shown in Fig. 4.
The governing equations are solved using finite volume method (FVM). As shown in Fig. 4, the inlet and outlet boundary conditions were set to total pressure inlet at atmospheric pressure and mass flow rate outlet. By setting different mass flow rates, the total pressure curve of system was obtained. Using the Multiple Reference Frame (MRF) model, the fan rotor was set as rotating domain with the rotating speed of 2882 rpm. The HE was set as porous medium zone, and its porous medium parameters were obtained according to the theoretical flow resistance curve of pressure drop and airflow velocity of HE. The airflow turbulence was described by k-ε model, and the equations were solved by Simple algorithm. When the residual error was less than 10 − 4 , it was considered to be convergent. The time step of transient simulation was 3.47 × 10 -4 s, which corresponding to 60 • of rotor rotating.
In order to ensure the simulation accuracy, mesh independence analysis was carried out. Under different mesh elements (3.19 × 10 6 , 3.66 × 10 6 , and 4.11 × 10 6 ), the difference of system total pressure was less than 1 % when the mesh elements exceeded 3.66 million. Considering the simulation cost and accuracy, the computational domain was finally meshed with 3.66 million elements.

Simplification of HE model
In the current work, the fin thickness of the HE was only 0.15 mm, which was several orders of magnitude smaller than the maximum size of the HE. If the actual HE model includes the fins was built, it is difficult to discretize the whole domain of the HE model. Therefore, the porous media model was used to simplify the HE model. The parameters of porous medium could be determined by Darcy model. According to the design theory of HEs [65], the relationship between the pressure drop and velocity can be determined by Eq. (4).
where a and b are the coefficients of the fitting equation. ΔP is the pressure drop and v is the velocity of airflow.
Porous media can be expressed as a source term including viscous loss term and inertia loss term. The equation of this source term is defined as follows: where S i is an additional source item. α and C 2 are the permeability and inertia resistance coefficient of porous medium, respectively. μ and ρ are the viscosity coefficient and density of airflow.

Model validation
According to the system pressure difference under the same flow rate with and without HE, the flow resistance characteristic curves of HE-OSF-30 can be obtained in the measured flow rate range in experiment. The comparison of theoretical and experimental flow resistance characteristic curves of HE-OSF-30 is shown in Fig. 5.
As shown in Fig. 5, for the five flow rate points tested in the experiment, the average flow resistance errors between theoretical and experimental of HE-OSF-30 is 8.5 %. It indicates that the theoretical calculation results of the HE resistance were accurate. In addition, the simulation results of the total pressure difference between the inlet and outlet of HE-OSF-30 under the flow rates of 15 m 3 /min, 30 m 3 /min and 45 m 3 /min are 32 Pa, 116 Pa and 231 Pa, respectively. Compared with the flow resistance at the corresponding flow rate on the theoretical flow resistance curve of HE, the errors are 7.8 %, 3.0 % and 8.7 %, respectively. That demonstrates the effectiveness of the simplified simulation method using the porous medium model to simplify the HE.   shown in Fig. 6, it can be found that the performance curves under both air flow rate adjusting modes were almost coincident when the volume flow rate was larger than 42 m 3 /min or smaller than 22 m 3 /min. But the performance curves showed significant differences in the range from 22 m 3 /min to 42 m 3 /min. The total pressure under adjusting mode 1 was obviously larger than that under adjusting mode 2. The total pressure presented a sudden drop from point 2 to point 3 under adjusting mode 1. While the total pressure showed a sudden rise from point 6 to point 7 under adjusting mode 2. The volume flow rate of the system reached the maximum value of 57 m 3 /min at point 1.

Performance of axial fan
For the case of adjusting mode 1, when gradually reducing the volume flow rate by controlling the throttle valve, the total pressure is gradually increased to the maximum value of 582 Pa. Then gradually decreased to 575 Pa and the fan stalled as the volume flow rate decreased to 31 m 3 /min at point 2 along the black line. Due to the rotating stall may cause serious harm to the fan, the motor would be shut off when the rotating stall occurred during the test. Then, keeping the position of throttle valve same as that before the motor shut off and restarting the motor, the operating point of the fan changed from point 2 to point 3. Afterwards, further reducing the volume flow rate, the fan working point would move from point 3 to point 4. And the system can operate normally until the flow rate was as low as about 12.5 m 3 /min.
For adjusting mode 2, when the volume flow rate was increased from 12 m 3 /min to 35.5 m 3 /min, the working point slowly changed from point 4 to point 6 along the red line. The total pressure gradually increased to 361 Pa. However, if a very small air flow rate was increased at point 6, the total pressure and the volume flow rate would increase suddenly to 563 Pa and 42 m 3 /min, respectively. The working point jumped from point 6 to point 7, and then the total pressure increased along with the volume flow rate decreased. The working point gradually moved to point 1.
For the air flow rate range between point 3 and point 7, as illustrated in Fig. 6, two total pressure performance curves under both adjusting modes seemed to form a loop, which was called hysteresis region of the axial fan. The mechanism forming the hysteresis region can be addressed by the aid of the delay effect of the post-stall in axial compression system. When the axial fan operated from the steady state to the unsteady state, the change could be completed within a very short time. But when it changed from unsteady state to steady state, it could take a long time to complete. That was, there would be a certain delay. The delay effect could lead to the hysteresis region of the axial fan.

Coupling performance of fan with different HEs
The total pressure curves of the oil cooling system composed of the axial fan and different HEs are demonstrated in Fig. 7. It can be clearly seen that the total pressure decreased with the volume flow rate increase for each plate-fin HE. Due to the flow resistance of the HE increased with the flow lengths increase for the same fin shape, the flow resistance of the rectangular fin HE with 30 mm, 60 mm and 90 mm flow lengths increased in turn. The increase was also suitable for the offset strip fin HE with 30 mm, 60 mm and 90 mm flow lengths. Since the offset strip fin could cause larger flow resistance than the rectangular fin, the flow resistance of the HE-OSF-30 was larger than that of HE-RF-60. Fig. 7 shows the total pressure of system decreases with the increase of HE flow resistance. When including with small flow resistance HEs (such as with all rectangular fins and 30 mm offset strip fins), the system would occur stall when the flow rate was adjusted with mode 1. And the volume flow rate of rotating stall was basically the same as that without HE. It can be known that, if stall occurred, the flow length and fin structure of HE had little effect on the fan stall boundary. However, when the flow resistance of HE was too large (with 60 mm and 90 mm offset strip fins), the system will not occur stall or surge. The system can still operate normally until the flow rate was far less than the stall flow rate of fan.
The fan performance curve and the theoretical flow resistance curves of offset strip fins HEs with the air side flow lengths of 30 mm, 60 mm and 90 mm are presented in Fig. 8. It can be found that when the HE flow length was 30 mm, the theoretical volume flow rate of the system working point was 54.1 m 3 /min. Compared with the theoretical result, the corresponding experimental result in Fig. 7 is 49.3 m 3 /min, with an error of 8.9 %. The experimental and theoretical results are in good agreement. In addition, increasing the system flow resistance through adjusting the throttle valve, the system of fan coupling with 30 mm offset strip fins HE will occur stall at the flow rate of fan stall. When increasing the air side flow length of HE, there may exist two intersections of the HE flow resistance curve and the fan performance curve. For example, when the air side flow length of the offset strip fins HE is 60 mm, the HE theoretical performance curve and the fan performance curve intersect at two working points where the volume flow rates are 31.8 m 3 /min and 41.0 m 3 / min, respectively. In the current experiment, as shown in Fig. 7, the operating volume flow rate of system is consistent with that of the working point with smaller flow rate, and the system can operate normally. Meanwhile, increasing the system flow resistance through adjusting the throttle valve, there is no surge or stall  phenomenon even flow rate is reduced to close 10 m 3 / min. When the air side flow length of offset strip fins HE is 90 mm, the system performance shows a similar result with that of 60 mm.

Heat transfer performance of the system in hysteretic region
The heat transfer performance of the oil cooling system is primarily characterized by the heat exchange capacity of the HE. For a working condition of the oil cooling system, the inlet temperature and the volume flow rate of the oil at the oil side of the HE was set to 120 • C and 60 L/ min, respectively. The inlet temperature of the cooling air at the air side was 20 • C. When including HE-OFS-60 and HE-OFS-90, the oil cooling system could operate at one of the two working points in the hysteresis region, as shown in Fig. 8 for each HE. According to thermal efficiency and number of transfer unit method [66], the heat transfer capacity of the system at both two working points for each HE could be determined, as illustrated in Fig. 9.
In Fig. 9, the air flow rate at both working points for the HE-OFS-60 was 31. Similarly, the air flow rate at both working points for the HE-OFS-90 was 25.1 m 3 /min and 34.2 m 3 /min. The heat exchange capacity at 25.1 m 3 / min and 34.2 m 3 /min was 23.1 kW and 27.5 kW, respectively. For both HE-OFS-60 and HE-OFS-90, the capacity at small air flow rate was 12.4 % and 19.0 % less than that at large air flow rate. Therefore, when the flow resistance of HE was large, attention should be paid to whether the system had two theoretical working points in fan hysteresis region. If yes, the system heat capacity should not be calculated basing on the working point with large flow rate. That may cause the system did not meet its cooling requirements due to it actually operated at a small flowrate.

The coupling mechanism of the system
The flow fields of the system including fan and HE-OSF-30 were simulated for the cases of 15 m 3 /min (small flow rate region), 30 m 3 / min (hysteresis region) and 45 m 3 /min (normal operation region). The HE was simplified by the porous medium model. The corresponding porous medium parameters were calculated according to the theoretical flow resistance of the HE. The velocity streamline and the total pressure field on the middle section of the system at 45 m 3 /min is shown in Fig. 10.
It can be seen from Fig. 10 (a) that the airflow was accelerated at fan rotor and there was swirl flow when it flowed into the transition duct between fan and HE. The local airflow velocity in the rotor area was more than 54.0 m/s. In addition, the airflow uniformly flowed out from HE (porous media zone) due to the rectification effect of the HE. The main flow velocity was less than 20.3 m/s. As shown in Fig. 10 (b), the total pressure of airflow in the fan rotor increased due to the work of rotor. And when the airflow flowed through the HE, there was an obvious pressure reduction process due to the resistance of the HE. As shown in Fig. 10(c), when the system operated at 45 m 3 /min, the inlet and outlet pressure of the HE were 548 Pa and 317 Pa, respectively. The corresponding pressure drop was 231 Pa. Compared with its theoretical  flow resistance of 253 Pa, the error was 8.7 %. It indicated that the simplified simulation method based on porous media model had good accuracy.
The turbulence kinetic energy in rotor outlet at different flow rates are shown in Fig. 11. As shown in Fig. 11 (a), the volume flow rate of 45 m 3 /min was belong to the normal operation flow rate range of fan. The turbulence kinetic energy in rotor outlet was low and basically the same at every blade passage. The maximum value of turbulence kinetic energy was 45.4 m 2 /s 2 and it was less than 10 m 2 /s 2 in most areas. It indicated that there was no boundary layer separation of airflow at the rotor blade. In Fig. 11(b), the turbulence kinetic energy in rotor outlet was not the same at each blade passage when the system operated at 30 m 3 /min. The maximum value was 105.8 m 2 /s 2 and the turbulence kinetic energy distribution had the fan stall characteristic. There was a certain boundary layer separation near the trailing edge of the rotor blade. It can be inferred that when the system operated in the fan hysteresis region, a certain boundary layer separation at the rotor blade led to an unsteady state flow. Therefore, when the air flow rate was adjusted by mode 1, the fan presented obvious stall phenomenon. But the fan did not stall when the flow rate was adjusted by mode 2. Fig. 11 (c) shows that the turbulence kinetic energy was basically the same at each blade passage at 15 m 3 /min. The turbulence kinetic energy behind all blades had a certain strength which was between turbulence kinetic energy when the system operated at 30 m 3 /min and 45 m 3 /min. It shows that the airflow also appeared separation near the rotor blade. However, the turbulence kinetic energy of airflow was low and basically the same at each blade passage. The separation of airflow is not strong enough to interfere with the airflow main flow field in fan rotor, thus, there was no obvious surge or stall phenomenon.

Conclusions
For the oil cooling system composed of the axial fan with front guide vane and plate-fin HE, the flow and heat transfer performance were investigated both experimentally and numerically. The fan total pressure performance was tested by two adjusting modes of air flow rate. Six types of plate-fin HEs with offset strip fin and rectangular fin were designed and manufactured for the flow lengths of 30 mm, 60 mm and 90 mm at the air side. The coupling characteristics of the axial fan and different HEs are analyzed in a systematic manner. The main conclusions can be achieved as follows: (1) The maximum total pressure and volume flow rate of the axial fan were 582 Pa and 57 m 3 /min. Under both adjusting modes of air flow rate, the hysteresis region formed in the range from 22 m 3 / min to 42 m 3 /min, which was near the stall boundary. (2) When the air flow rate decreased to a certain value, the stall would occur in the system coupling the rectangular or small resistance offset strip HEs. But the stall did not appear for the system coupling the large resistance HEs. The air flow rate at the stall point for the system with different HEs was almost the same as that without HE. The fin structure and the air side flow length of HE have no effect on the fan stall boundary. (3) The system could have two theoretical working points in the fan hysteretic region. The corresponding heat transfer capacity was significantly different at both two working points. For the system including HE-OSF-90, the air flow rate at both working points were 25.1 m 3 /min and 34.2 m 3 /min. The heat transfer capacities were 23.1 kW and 27.5 kW, respectively, with a difference of 19.0 %. In actual operation, the system operated at the working point with smaller flow rate. (4) The boundary layer separation of the airflow occurred at the rotor blade when the system flow rate was far less than the fan stall flow rate. Due to the turbulence kinetic energy of separation was low and basically the same at each blade passage, the system had no obvious surge or stall when operating at small flow rate.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.