## Abstract

In order to improve the fine optimization technology of flow measurement facilities in irrigation areas, the hydraulic characteristics optimization system of an airfoil pillar-shaped flume is constructed. The optimization system integrates three modules of airfoil shape reconstruction by Hicks–Henne function, CFD numerical simulation calculation, and NSGA-II algorithm. It can efficiently find the shape design variable value that makes the hydraulic characteristics optimal. An example of the airfoil pillar-shaped flume in a rectangular channel at 30 L/s of discharge rate is given, and the specific optimization design concept of backwater height and head loss is presented. The results show that the backwater height is reduced by 9.49% compared with the prototype, and the head loss is reduced by 8.10% compared with the prototype, proving the feasibility of the optimization system. It can provide a theoretical reference for the optimization of hydraulic characteristics of similar measuring flumes.

## HIGHLIGHTS

The airfoil pillar-shaped flume itself has the basic advantages of convenient installation and streamlined flow characteristics.

The hydraulic characteristics optimization system integrating the Hicks–Henne function, numerical simulation, and NSGA-II algorithm is designed.

The design concept of the optimization case of the backwater height and head loss of the airfoil pillar-shaped flume is concretely presented.

## INTRODUCTION

Open channel flow measurement is the basis for water resource allocation and water fee collection in irrigation areas and plays a vital role in agricultural irrigation management. Among many flow measurement facilities, a series of measuring flumes developed based on the Venturi principle have been widely used because of their advantages of low construction cost, reliable measurement accuracy, and wide application range. The principle of flow measurement is that when the water flows through the channel section narrowed by the measuring flume, its flow velocity will increase, the water level will drop, and a critical flow state will be formed near its throat. Under this state, there is a stable and single stage–discharge relationship. However, because of the differences in their hydraulic characteristics corresponding to various facility structures, as well as the improvement of the applicable requirements of the measuring facilities in irrigation area management, research on the improvement and optimization of the flow measurement facilities remains an urgent problem to be solved.

The earliest Venturi flume belongs to the short-throat flume. After comparing several different ratios of throat widths to flume lengths, throat lengths, and arrangements of end wings, the results showed that a greater length of converging and diverging section and a rounding of the throat section would result in less head loss and great accuracy in measurement of flow, but the standard was chosen as a compromise between accuracy and cost (Cone 1917). The Parshall flume improved the throat length and slope of the original Venturi flume, which made critical flow easier to occur and improved the measurement accuracy (Parshall 1926). But the overall structure of the Parshall flume is complex, and the construction requirements are high. Based on this problem, Skogerboe & Hyatt (1967) removed the throat section and designed the standard shape of the cutthroat flume, which makes the flume simple and economical. However, the cutthroat flume will cause large head loss under free flow conditions (Ramamurthy *et al.* 1988). In order to find flow measurement facilities with wider adaptability, the simplified NACA airfoil equation profile is introduced into the flow measurement field (Wang & Wen 1990), and it has the advantages of smooth flow, small head loss, high critical submergence, and high flow measurement accuracy (Hong *et al.* 2005; Lu *et al.* 2006; Liu *et al.* 2008).

In consideration of the difficulty of building or installing the measuring flume in different channels and the convenience of flow measurement, a cylindrical measuring flume was designed to be fixed in the center of the channel, which is applicable to all symmetric prismatic channels (Hager 1985). Then, the cylindrical diffuser section was modified into a V-shaped tail, which improved the flow pattern and reduced head loss (Liu *et al.* 2013). Moreover, the principle of optimizing or designing a cylindrical moving flume was discussed by comparing it with other types of the cylindrical moving flume (Li *et al.* 2020). In addition, Peruginelli & Bonacci (1997) also proposed a moving pier-shaped prism structure in the center of a rectangular channel, which was later called the central baffle flume. The influence of different geometric parameters of the central baffle flume was studied through experiments (Kolavani *et al.* 2019) and the optimization design results were verified (Bijankhan & Ferro 2019). According to the central layout principle, Liu *et al.* (2019) conducted comparative experiments and simulations on the application performance of the airfoil pillar-shaped flume in three common channels. The results show that the advantages of smooth flow, small head loss, high flow measurement accuracy, and high submergence are still retained. However, the measuring flume is an additional device relative to the channel. In the actual project, the construction of the channel and the installation of the measuring flume cannot be carried out simultaneously. The installation of the measuring flume will have an impact on the channel flow, causing the backwater upstream and head loss. Especially for the plain irrigation area or the environmental conditions with poor head, the airfoil pillar-shaped flume, like other types of measuring flume, still faces the problem of water blocking.

With the development of demand and the improvement of research technology, numerical simulation technology has been widely used in the field of alternative model computing, which brings more possibilities for innovation. More importantly, the optimization concept has been demonstrated in other fields such as water resource allocation optimization, aerodynamic performance optimization, hydrodynamic performance optimization, and material structure optimization by integrating optimization algorithms and numerical simulation technology (Sadeghi-Tabas *et al.* 2017; Yang *et al.* 2018; Yan *et al.* 2019; Nozari *et al.* 2021; Yin *et al.* 2021; Wang *et al.* 2022). Taking optimization algorithm as the intelligent guide of numerical simulation technology improves the efficiency, fineness of optimization, and rational utilization of computing resources. However, from the above research, it can be seen that this concept is rarely applied to the improvement of flow measurement facilities. At present, the optimization of hydraulic characteristics of measuring flumes is mostly the result of comparative analysis by presetting several limited test groups, which cannot be determined to be the relative optimal or global optimal in a specific area. The reason is that in the traditional method, when designing the test groups for comparison, the design variable values to control the geometric structure of the facility are selected at equal intervals, and then the best group is selected by a unified comparison of the results. However, the optimization algorithm first randomly selects the design variable values and can compare each group of results with the previous results in time. It makes the selection direction of the next design variables close to the area with good performance, and more computing resources can be concentrated in the effective area, so as to improve the efficiency of optimization and the fineness of the geometric structure adjustment of the facility. Therefore, research on intelligent and fine optimization of flow measurement facilities will become an inevitable development trend in irrigation areas.

To summarize, the difficulty of construction or installation of the measuring flume, as well as the impact on the flow of the original channel following the installation of the measuring flume, must be considered concurrently. Based on the airfoil pillar-shaped flume with the advantage of streamlining, the optimization design system of hydraulic characteristics is established by integrating the parameterization method of the airfoil, computational fluid dynamics (CFD) numerical simulation technology, and optimization algorithm in this paper. It can improve hydraulic performance and expand the application scope as much as possible under limited hydraulic conditions. The design concept can provide a certain theoretical reference for the optimization research of similar flumes.

## MATERIALS AND METHODS

### The airfoil shape reconstruction

#### The structure of an airfoil pillar-shaped flume

The airfoil structure is the main factor that affects the flow characteristics of the airfoil pillar-shaped flume, so the optimization of hydraulic characteristics is essentially to find the optimal design scheme for the airfoil shape. The link relationship between airfoil shape and hydraulic characteristics must be transformed from physical property problems into mathematical model problems before continuing research. Therefore, choosing an appropriate airfoil expression method is a prerequisite for optimization.

*C*and the maximum airfoil thickness

*P*, and the maximum airfoil thickness is located at 0.3

*C*of the upstream section of the building. The base shape function is as follows:where

*C*is the chord length (m) and

*P*is the maximum airfoil thickness (m).

#### The airfoil parameterization

The parameterization of the airfoil is to express the airfoil to be optimized through a finite number of functions, and then adjust the airfoil shape by controlling the variable parameter values of the function. The optimization process is to find the optimal solution of the variable parameters. On the basis of ensuring the smoothness and practicability of the airfoil, this paper selects the most widely used Hicks–Henne function method (Hicks & Henne 1978) in the linear superposition method of analytic functions to express the airfoil. The airfoil shape is controlled by parameterizing and superimposing the variation of the airfoil curvature and thickness, which not only changes the key points of the airfoil profile curve, but also ensures the smoothness of the airfoil surface. The Hicks–Henne function expresses the geometry of the airfoil in three parts: the basic airfoil , the shape function and the design variable . The shape function is also called the disturbance function, and the design variable is used as the amplification factor of the shape function, and its role is to constrain the disturbance range.

*n*is the number of shape functions . In general, the more the selected shape functions, the finer the shape adjustment of each key point of the reference airfoil. However, the more design variables corresponding to the shape function, the more difficult it is to find the optimal solution. In this paper, seven shape functions are used to disturb the basic airfoil.

*e*(

*k*) = lg0.5/lg

*x*, 0 ≤

_{k}*x*≤ 1 (assuming the airfoil chord length is 1), when

*k*= 2, 3, 4, 5, 6, 7, takes 0.15, 0.30, 0.45, 0.60, 0.75, and 0.90, respectively, correspond to the abscissa positions of the peak points of the shape function, as shown in Figure 2. When

*x*in the shape function takes 0 and 1, the value of is 0, which can ensure that the position of the front and rear trailing edges of the airfoil remains unchanged during the disturbance.

The MATLAB software is applied to realize the airfoil shape reconstruction module in this study. First, the co-ordinate data corresponding to each point of the basic airfoil is extracted. Then the value of the design variable of each key point is brought into Equation (3) to calculate the disturbance data corresponding to the ordinate, and brought into Equation (2) to realize the airfoil adjustment. Finally, the new airfoil co-ordinate points are converted into a file format recognized by the modeling software.

### Automated CFD numerical simulation

#### Model building and meshing

In this study, the ICEM software is used for modeling and meshing the fluid domain. The ICEM has the function of recording scripts, which can record all relevant commands involved in the normal operation process. After executing the script recording function, call the new airfoil co-ordinate data generated by the airfoil shape reconstruction module to establish the model, divide the fluid domain, generate the mesh block, establish the mapping relationship, define the meshing method and mesh size, and output mesh file, finally generate a script file. When the airfoil is adjusted within a small range, it does not affect the original mapping relationship between the fluid domain model and the block. Just run script file again to get the corresponding model output mesh file automatically.

#### Selection of numerical solution method

^{3});

*t*is the time (s);

*p*is pressure (Pa); is the fluid dynamic viscosity (Ns/m

^{2}); is the stress tensor (Pa); is the component of the fluctuation velocity (m/s).

*k*–

*ε*model, which has been verified in previous study to have a good simulation effect (Sun

*et al.*2021). The transport equations of the turbulent kinetic energy

*k*and dissipation rate

*ε*corresponding to this turbulence model are:where is the generation of turbulence kinetic energy due to the mean velocity gradients; is the generation of turbulence kinetic energy due to buoyancy, is the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate; is the turbulent viscosity, , ; , , and are constants, = 1.44, = 1.92; and are the turbulent Prandtl numbers for

*k*and , , .

In addition, for the capture method of free liquid surface, this study adopts the VOF model. The model achieves the tracking of the interphase interface in the computational domain by introducing the variable phase volume fraction, which is very suitable for calculating fluid flows such as air and water that cannot be mixed with each other. The sum of the phase volume fractions of a unit grid is 1, and the free water surface is formed by finding the combination of grids with two phase volume fractions of 0.5.

The FLUENT software also has a secondary development function, it can read two types of commands such as GUI and TUI. The GUI commands are similar to the script recording function of ICEM software. Script recording can be realized by executing the Start Journal function, which is automatically written in Scheme language. The TUI commands are written in C programming language but cannot be recorded and written according to the user interface. The difference between the two is that the GUI command file can be automatically generated, and the familiar operation is simple, but when reading the Journal file again and executing the corresponding GUI command, it is easy to be disturbed by third-party operations and cause confusion in the settings. The TUI commands are relatively concise and without redundant operations, it is not affected by other operations when executing the commands, but users need to be familiar with the relevant setting command lines in order to write related commands. Considering the stability of settings, this module uses the TUI commands to realize automatic settings such as reading mesh files, adding materials, selecting turbulence models, defining boundary conditions, initializing calculations and file output.

### Optimization algorithm

#### Selection of optimization algorithm

As a means to find the optimal solution, optimization algorithm is a key part of the optimization system. It can constantly adjust the search direction in the design space until it finds the shape scheme that can make the hydraulic performance of the airfoil pillar-shaped flume reach the best, which greatly affects the optimization efficiency and the quality of the optimization results. Moreover, the optimization problem is generally a multi-objective problem in practical engineering. The NSGA-II algorithm is selected as the multi-objective problem solving method in this paper, which belongs to the genetic algorithm (GA).

**.**The algorithm is based on the simple GA and carries out an additional operation before the selection operation, that is, the non-dominated hierarchical method. According to the dominant and non-dominant relationships among individuals, the hierarchical virtual fitness is normalized by fitness sharing strategy. It maintains the advantage of good individual fitness while also maintaining the diversity of the population, but requires decision makers to specify a sharing radius. Therefore, the Fast Elitist Non-dominated sorting genetic algorithm (NSGA-II) was proposed on the basis of NSGA (Deb

*et al.*2002). First, it adopts a fast non-dominated sorting algorithm, which reduces the computational complexity compared to NSGA. Furthermore, the crowding degree and the crowding degree comparison operator are used to replace the fitness sharing strategy that needs to specify the sharing radius. At this time, the virtual fitness value of the generated subgroup is consistent with the non-dominated sorting result. Finally, the elite strategy is introduced, which combines the parent and child populations to compete together, expands the sampling space, prevents the loss of the best individual, and improves the computational speed and convergence of the algorithm. Therefore, considering the convergence speed and accuracy of the optimization algorithm, the NSGA-II is selected as the optimization algorithm for this study.

#### Optimization algorithm settings

The content of optimization algorithm settings includes objective fitness function design, constraint condition setting and algorithm parameter setting. First, designing a good fitness function according to the objective function is an important part of the optimization system, and it is also the basis for judging the pros and cons of the design scheme. For the airfoil pillar-shaped flume, the optimization objectives of its hydraulic characteristics can mainly be selected from the quantitative parameters of upstream backwater height, Froude number, submerged degree, and head loss. The fitness function can be reasonably designed by the penalty function method according to the physical parameters or flow characteristics of the airfoil pillar-shaped flume and channel. The penalty function method refers to adding a penalty constant to the original objective function to obtain an augmented objective function. Its function is to give a maximum value to the non-feasible point or the point trying to cross the boundary and escape from the feasible region, that is, to convert the constrained optimization problem into the unconstrained optimization problem. After the fitness function is determined, it is necessary to impose constraints on the range of independent variables or dependent variables of each module according to the design requirements to control the degree of change. The parameter settings of the optimization algorithm are factors such as the size of the design population, selection, crossover, mutation characteristics, and convergence criteria, which affect the computational efficiency of finding the optimal solution.

*c*

_{1}–

*c*

_{7}) are initialized randomly by the NSGA-II algorithm, and then they are brought into the corresponding equations in the airfoil reconstruction module to obtain the new airfoil shape, and the co-ordinate points corresponding to the shape constraints are output for judgement. If the shape constraints are not satisfied, the objective fitness function will be affected by the penalty function method, resulting in a bad result, and the direction of the design variables update by the NSGA-II algorithm will be far away from the surrounding area of the last selected values when entering the next iteration. If the shape constraints are satisfied, enter the numerical simulation solution module to output the hydraulic parameters and judge the corresponding constraint conditions. If the constraints of hydraulic parameters are not satisfied, a result with bad fitness value is also obtained to guide the direction of design variables update in the next iteration. If the constraints of hydraulic parameters are satisfied, it will not be affected by the penalty function method to get a better result, and then the next step is to judge whether the number of iterations meets the requirements. If the requirements are not met, continue to the next iteration. If the iteration requirements are met, the airfoil shape corresponding to the best fitness value is output, that is, the optimal airfoil shape.

## OPTIMIZATION CASE DESIGN

### Optimization object function

In this paper, the airfoil pillar-shaped flume in the rectangular channel is selected as the research object, aiming at the problem of water blocking in the flow measurement in plain irrigation area, the airfoil shape is optimized with the reduction of the backwater height and the head loss as the hydraulic characteristics optimization objectives. According to the actual situation, the design discharge rate in the optimization process is taken as 30 L/s, and the relevant design parameters of the rectangular channel and the airfoil pillar-shaped flume are shown in Table 1.

Rectangular channel . | Basic airfoil pillar-shaped flume . | |||
---|---|---|---|---|

Width B (m)
. | Depth H (m)
. | Length L (m)
. | Chord length C (m)
. | Maximum airfoil thickness P (m)
. |

0.6 | 0.3 | 12 | 1 | 0.15 |

Rectangular channel . | Basic airfoil pillar-shaped flume . | |||
---|---|---|---|---|

Width B (m)
. | Depth H (m)
. | Length L (m)
. | Chord length C (m)
. | Maximum airfoil thickness P (m)
. |

0.6 | 0.3 | 12 | 1 | 0.15 |

*h*is the backwater height (m), m; is the water depth of the upstream stable section before the installation of the flume (m), and the corresponding under the 30 L/s discharge rate in this paper is 0.0957 m; is the water depth of the upstream stable sections after the installation of the flume (m).

### Numerical simulation

The fluid domain of the airfoil pillar-shaped flume is symmetrical, so half of the original fluid domain is modeled in ICEM, and the fluid domain is meshed by the Cartesian co-ordinate method. For the two hydraulic characteristics of upstream backwater height and head loss in this study, the key is to extract the water depth data of the upstream stable section and the downstream stable section. The Cartesian co-ordinate meshing method is not disturbed by the change of the airfoil shape, and it can effectively ensure the uniformity of mesh quality in these two key regions. Set the side lengths to 0.025, 0.02, 0.015, and 0.01 m, respectively, for the mesh independent analysis of the cube mesh size, the results are shown in Table 2. As the mesh size decreases, the relative error between meshes becomes smaller and smaller. Considering the requirements of computational accuracy and computational efficiency, this study finally decided to use a mesh size of 0.015 m.

Mesh size (m) . | Number of meshes . | Water depth at 2.5 m (m) . | Relative error (%) . | Water depth at 8.0 m (m) . | Relative error (%) . |
---|---|---|---|---|---|

0.025 | 123,060 | 0.1543 | – | 0.0786 | – |

0.02 | 169,932 | 0.1519 | 1.55 | 0.0771 | 1.91 |

0.015 | 914,254 | 0.1505 | 0.92 | 0.0762 | 1.17 |

0.01 | 1,235,058 | 0.1496 | 0.60 | 0.0756 | 0.79 |

Mesh size (m) . | Number of meshes . | Water depth at 2.5 m (m) . | Relative error (%) . | Water depth at 8.0 m (m) . | Relative error (%) . |
---|---|---|---|---|---|

0.025 | 123,060 | 0.1543 | – | 0.0786 | – |

0.02 | 169,932 | 0.1519 | 1.55 | 0.0771 | 1.91 |

0.015 | 914,254 | 0.1505 | 0.92 | 0.0762 | 1.17 |

0.01 | 1,235,058 | 0.1496 | 0.60 | 0.0756 | 0.79 |

*k*–

*ε*turbulence model and VOF multiphase model, and configure gravity and boundary condition. The channel inlet is set to velocity inlet, the outlet is set to pressure outlet, the air inlet is set to pressure inlet, the central axis plane is set to symmetry, and the rest of the boundaries are set to non-slip solid walls.

### Optimization algorithm settings

In this study, the NSGA-II algorithm is used as the guide for discriminative optimization in the operation framework. The content of optimization algorithm setting includes three aspects: objective fitness function design, constraint condition setting, and algorithm parameter setting.

#### Objective fitness function design

According to the optimization objective and hydraulic parameter characteristics, the fitness function is designed as follows:

#### Constraint condition setting

The value ranges of the design variables are shown in Table 3. The value standard is to ensure that each co-ordinate point of the basis airfoil can move between the minimum and maximum values, that is, the active area of the corresponding airfoil curve is a rectangle (1 m × 0.1 m).

Design variables . | c_{1}
. | c_{2}
. | c_{3}
. | c_{4}
. | c_{5}
. | c_{6}
. | c_{7}
. |
---|---|---|---|---|---|---|---|

Upper limit | 0.4190 | 0.0164 | 0 | 0.0105 | 0.0359 | 0.0710 | 0.1138 |

Lower limit | −0.1662 | −0.1336 | −0.1500 | −0.1395 | −0.1141 | −0.0790 | −0.0362 |

Design variables . | c_{1}
. | c_{2}
. | c_{3}
. | c_{4}
. | c_{5}
. | c_{6}
. | c_{7}
. |
---|---|---|---|---|---|---|---|

Upper limit | 0.4190 | 0.0164 | 0 | 0.0105 | 0.0359 | 0.0710 | 0.1138 |

Lower limit | −0.1662 | −0.1336 | −0.1500 | −0.1395 | −0.1141 | −0.0790 | −0.0362 |

#### Algorithm parameter setting

The parameter values of the NSGA-II algorithm are set according to the optimization requirements and experience, as shown in Table 4.

Population size . | Number of generations . | Crossover probability . | Crossover distribution index . | Mutation distribution index . |
---|---|---|---|---|

60 | 50 | 0.9 | 10 | 20 |

Population size . | Number of generations . | Crossover probability . | Crossover distribution index . | Mutation distribution index . |
---|---|---|---|---|

60 | 50 | 0.9 | 10 | 20 |

## RESULTS AND DISCUSSION

Shape . | Discharge Q (L/s)
. | Upstream (m) . | Downstream (m) . | Backwater height (m) . | Head less (m) . |
---|---|---|---|---|---|

Basic airfoil | 30 | 0.1505 | 0.0762 | 0.0548 | 0.0580 |

Optimized airfoil | 30 | 0.1453 | 0.0759 | 0.0496 | 0.0533 |

Shape . | Discharge Q (L/s)
. | Upstream (m) . | Downstream (m) . | Backwater height (m) . | Head less (m) . |
---|---|---|---|---|---|

Basic airfoil | 30 | 0.1505 | 0.0762 | 0.0548 | 0.0580 |

Optimized airfoil | 30 | 0.1453 | 0.0759 | 0.0496 | 0.0533 |

As can be seen from Figure 6, in terms of shape, compared with the basic airfoil, under the condition of fixed airfoil chord length, the throat position of the optimized airfoil moves backward from the position of 0.3 m to the position of 0.46 m, making the contraction section longer and the curvature of its contour curve smaller, while the corresponding diffusion section is shorter and the curvature of its contour curve larger. It can be seen from Table 5 that the backwater height of the airfoil pillar-shaped flume optimized in terms of hydraulic characteristics is reduced from 5.48 to 4.96 cm, with a relative reduction of 9.49%. The head loss decreased from 5.80 to 5.33 cm, with a relative reduction of 8.10%. According to the analysis, the change trend of the shape is to improve the water contraction effect in time and space as much as possible on the basis of taking into account the diffusion effect of the downstream water flow, which reduces the disturbance and obstruction effect of the flume structure on the upstream water flow, thus retaining more flow velocity head and making the upstream backwater lower than the prototype. To sum up, the analysis proves that the hydraulic characteristics optimization system for the airfoil pillar-shaped flume is reliable.

## CONCLUSION

With the continuous improvement of the requirements for flow measurement facilities in the water transmission and distribution system of the irrigation area, it is necessary to reasonably improve the existing facilities to reduce the impact on the flow and adapt to more severe hydraulic conditions. In order to make the work more systematic, refined, and intelligent, this paper designs an optimization system of hydraulic characteristics based on the airfoil pillar-shaped flume as the research object, which is mainly composed of three modules. Firstly, this paper uses the Hicks–Henne function in the airfoil parameterization method to ensure the smooth transition of the curve while disturbing the adjustment of the basic airfoil curve. Then the flow field of the airfoil pillar-shaped flume is automatically simulated by using the numerical simulation technology and the required hydraulic parameters are extracted. In the optimization process, the NSGA-II algorithm is used to properly constrain the design variables and hydraulic parameters, and the penalty function method is used to process the fitness value to guide the optimization direction in a timely manner.

The optimization design concept of the system is demonstrated concretely by using the example of the airfoil pillar-shaped flume under a discharge rate of 30 L/s. The optimization results show that the backwater height is reduced by 9.49% compared with the prototype, and the head loss is reduced by 8.10% compared with the prototype. Furthermore, the optimized airfoil also performs well in hydraulic characteristics under other discharge conditions. These prove the feasibility of the optimization system. In addition, it is possible to flexibly adjust each module in the system according to the demand and actual situation. For example, other parameterization methods can be used in the shape reconstruction module to represent different flow measurement facilities, different channel fluid domain models can be established in numerical simulation, and other optimization objectives and constraints can be considered in the optimization algorithm.

## ACKNOWLEDGEMENTS

This research was supported by the National Key Research and Development Program of China (No. 2017YFC1501204), the National Natural Science Foundation of China (No. 51909242), the Program for Science and Technology Innovation Talents in Universities of Henan Province (No. 19HASTIT043), and the Outstanding Young Talent Research Fund of Zhengzhou University (1621323001).

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## CONFLICT OF INTEREST

The authors declare there is no conflict.