Preprint
Article

Optimization of water-jet pump based on the coupling of multi-parameter and multi-objective genetic algorithms

Altmetrics

Downloads

150

Views

54

Comments

0

A peer-reviewed article of this preprint also exists.

Submitted:

26 July 2023

Posted:

27 July 2023

You are already at the latest version

Alerts
Abstract
As the core components of water jet propulsion system, water jet pump, directly affect the ship's propulsion performance. 3D reverse design method can effectively suppress the impeller secondary flow to get better performing hydraulic models. Aiming at fast and efficient optimization of the impellers, we parameterize the impellers and runner, and coupled with optimization algorithms and numerical calculations. It is an efficient and feasible way to coupling optimization based on multi-parameter, multi-objective and optimization algorithms, and parameterizing the impellers and guide vane geometry based on the 3D model designed by 3D reverse design method. In this paper, the waterjet propulsion pump is taken as the research object. According to the design method of the waterjet propulsion pump hydraulic model, the traditional axial flow blade design method, the modern ternary design theory of the blades, and the parametric three-dimensional reverse design method are selected to constructing an optimal design optimization strategy for the impeller of water-jet propulsion pump based on multi-parameter and multi-object genetic algorithm. Under the design conditions, the optimized pump efficiency reaches 89.28%, which is 1.24% higher than the original hydraulic model. The optimized pump head is 13.51 m, which is 0.39m more than the original hydraulic model. Under the design condition, the test measured head is 13.26m, the numerical calculation result is 1.88% higher than the test result, the test measured efficiency is 87.1%, the numerical calculation result is 2.5% higher than the test result. From the test results, the optimized impeller has higher hydraulic performance. Experiments have shown that it is feasible to optimize the design strategy.
Keywords: 
Subject: Engineering  -   Marine Engineering

1. introduction

Water-jet pumps are widely used in military ships and weapon launchers. the water-jet pump uses the reaction force of the water-jet from the pump to propel the ship forward. the mixed flow pump is widely used in the field of water-jet. As a core component of the water-jet propulsion system, it directly affects the overall propulsion performance of the ship. Although the water-jet pump is required to operate in a relatively narrow range of working conditions, its design is very challenging. Specifically, because the propulsion efficiency is directly proportional to the efficiency of the pump, the water-jet pump must ensure efficient hydrodynamic performance, and the fluid output of the water-jet pump must be uniform, axial flow, and no circulation to achieve water spray the thrust is maximized [1].Cavitation will cause the hydraulic performance of hydraulic machinery to drop sharply, that is, cavitation fracture, and also lead to the destruction of the material surface [2].
the improvement of water-jet technology is based on the development of high efficiency and high cavitation resistance of water-jet pumps. the design of the water-jet pump is not only limited by the space size of the ship, but also to ensure efficient and low noise operation at high speed.
At present, there are two main approaches to blade design methods for rotating machinery, namely positive design and inverse design. the positive design refers to obtaining the 3D flow field through a given blade, while the inverse design refers to calculating the corresponding blade geometry by specifying the flow field distribution. in recent years, with the development of CFD, great progress has been made in the numerical calculation of the design of rotating machinery blades. It is indeed an effective method to modify the blade shape and apply the CFD method to predict the hydraulic characteristics of rotating machinery. However, there are actually many difficulties in determining the magnitude and direction of blade shape modification. This difficulty mainly reflects that small changes in any position of the blade will affect the flow distribution of other parts of the blade, and the relationship between the geometric parameters of the blade and the flow field is complicated. therefore, it is a more reasonable way to design the rotating machinery blades by the 3D inverse design method.
in recent years, advances in experimental techniques and numerical calculations have facilitated a better understanding of the flow mechanism. in order to suppress the secondary flow in the impeller, Zangeneh [3,4] proposed a compressible full three-dimensional inverse design method for the design of centrifugal and mixed flow rotating machinery in 1991, which defined the impeller. the distribution of the circumferential average velocity circulation on the meridian flow channel, and the corresponding blade shape is solved iteratively. the inverse design method is based on the non-stick slip condition to make the flow depend on the surface of the blade. This method is widely used in the secondary flow suppression of the mixed-flow pump impeller [5,6,7,8,9].
Automatic optimization technology can reduce the development time in the design process and explore the design space in a more systematic way. Automatic optimization technology has been widely used in the design of aerodynamics and rotating machinery, but considering its high computational cost, its application is only limited to simple applications in industry. the automatic optimization technology of rotating blades relies on three main modules: geometric parameterization, performance evaluation and optimization algorithms. A good geometric parameterization method must represent the widest design space with the fewest parameters. the shape of the blade is usually characterized by direct geometric shapes, the most famous of which are simple B-spline curves [10,11,12,13,14], NURBS (non-uniform rational B-spline) surfaces [15] and quasi three-dimensional sections defined by Bezier curves [16]. in order to reduce the computational complexity and avoid invalid blades, it is generally possible to reduce the degree of freedom of design parameters as much as possible, that is, to fix some parameters, such as stacking, bone lines [12], the thickness of the inlet and outlet edges of the blade, the inclination angle, etc. [16]. the optimization algorithm must be able to seek the world's largest multi-dimensional, multi-mode, nonlinear and discontinuous design space, avoiding local maxima or noise. Due to robustness and simple theory, evolutionary algorithms (EA) are the most popular, including evolutionary strategies (ES) [12], genetic algorithms (GA) [10,15] and multi-objective genetic algorithms (MOGA) [11,12,17]. Adaptive simulated annealing algorithm [15] is also an effective method. Another method based on experimental design and agent model coupling is also gradually mature, such as response surface [18] or artificial neural algorithm [19]. Bonaiuti and Zangeneh [18] proposed an optimization strategy, that is, to realize three-dimensional inverse design method, response surface method, multi-objective optimization algorithm and CFD analysis coupling optimization. the parameterization of the blade is realized through a three-dimensional inverse design method, and the dynamic characteristics and performance are controlled by optimizing the parameterized blade, and the optimized blade performance is significantly improved. the coupling of response surface method and inverse three-dimensional design method can also be used to analyze the sensitivity of various parameters to performance, which will provide guidance for design and improvement. Boselli and Zangeneh [20] proposed a 3D multi-objective multidisciplinary design method that combines a 3D inverse design method and a multi-objective genetic algorithm, and is applied to the design of axial-flow turbine blades. the blade is parameterized by the blade load parameters, span load distribution and maximum thickness, and the sample space is constructed from these few parameters. Huang [21] uses a modified NSGA-Ⅱ algorithm coupled with a dynamic crowding distance (DCD) to achieve multi-objective optimization of water-jet propellant mixed-flow pumps. the meridian flow channel of the fixed pump uses ten-dimensional inverse design methods to select ten variables of the rim and hub to control the blade load. the optimization goals are hydraulic efficiency and impeller head. the radial basis neural network (RBNN) is used to approximate the objective function with 82 training samples. Local sensitivity analysis shows the impact of key variables on pump performance. the optimized model impeller reaches 92%, and the high efficiency area is wider, and the suction surface speed gradient is smoother.
Bonaiuti [22] believes that the CFD method and the automatic optimization method are becoming more and more popular, but the traditional design method based on the physical analysis of flow development and the subsequent modification of geometric parameters is still widely used in industrial design because of the characterization of design quality. there are still great difficulties in converting a large number of factors into mathematical parameters. It adopts CFD method for hydrodynamic design of mixed-flow pump, and its purpose is to improve the efficiency and head of the pump in the whole working range. the geometry of the impeller and guide vane is controlled in a parameterized manner, that is, the influence of key parameters on performance is analyzed simultaneously through several parameters and Bezier curves.
in summary, the 3D inverse design method can effectively suppress the secondary flow of the impeller and obtain a hydraulic model with better performance. Parameterizing the blades and flow channels, coupled with optimization algorithms and numerical calculations, can achieve fast and efficient optimization of the impeller. the initial model of the pump is designed using a 3D inverse design method. by parameterizing the geometry of the impeller and guide vane, coupling optimization based on multi-parameters, multi-objectives and optimization algorithms is an efficient and feasible method.
This paper will take the water-jet pump as the research object, and choose the comparison of traditional axial flow blade design method, modern blade ternary design theory, and parameterized 3D inverse design method for the design method of the water-jet pump impeller hydraulic model. An optimized design and optimization strategy for the impeller of water-jet pump based on multi-parameter and multi-object genetic algorithm was constructed, and a set of hydraulic models with better performance was obtained, and the feasibility of the optimized design strategy was verified through experiments. the optimized model will be used for the follow-up visualization test and numerical calculation method of water-jet pump cavitation.
Figure 1. Water-jet pump hydraulic model design process.
Figure 1. Water-jet pump hydraulic model design process.
Preprints 80601 g001

2. Physical Model

the hydraulic components of the water-jet pump are composed of the following four components: suction, impeller, diffuser and outpipe, as shown in Figure 2, where the impeller is a mixed flow impeller and the diffuser is a constricting guide vane.

2.1. Pump Design Parameters

the main design parameters of the water-jet pump are listed in Table 1, including flow rate, head, rotating speed.

2.2. Requirements for the Design Geometric Parameters of the Pump

Considering the space constraints of the ship, the flow channels and structural dimensions of the hydraulic components of the water-jet pump are also restricted. the number of blades of the impeller is not less than 5, and the number of guide vanes is not less than 9. the inlet and outlet diameter of the runner is Dj = 270 mm, Dd = 250 mm, respectively. the maximum diameter of the hub is not less than 250 mm, and the diameter of the hub at the inlet side of the impeller blade DhLE is not less than 120 mm. the maximum diameter of the rim Ds is not more than 370 mm. the total length L of the pump is not more than 400 mm. the layout of flow channel size parameters is shown in Figure 3. the flow components include the impeller, guide vane, inlet section, and outlet section. Add 1080 mm length straight pipe (4Dj) at pump inlet, add 1000 mm length straight pipe (4Dd) at pump outlet. the pressure measurement cross section is taken in the middle of the pipe.

3. CFD calculation Model Requirements

3.1. Establishment of the initial Model

After determining the hydraulic design parameters and geometric flow channel of the water-jet pump, the method of determining the shape of the water-jet pump impeller and diffuser blades need to be considered.
At present, lift method, streamline method or quasi-three-component method design are all based on satisfying the hydrodynamic performance, and based on the flow of the cylindrical cascade. these methods are impossible to effectively control the radial flow direction of the blade and pulsation velocity component in the blade channel, and cannot constrain the blade load distribution, pulsation pressure distribution and cavitation performance [23]. therefore, in order to quickly obtain an initial model with better performance, the design of the initial model adopts the blade parametric 3d inverse design method.
initial conditions given:
(1) the axial projection shape of the blade, including hub, shroud, leading edge, trailing edge.
(2) Shaft speed n and blade number B, guide vane blade shaft speed is 0;
(3) Fluid density ρ, blade inlet reference pressure p0;
(4) Design flow rate Q;
(5) the distribution law of the annular volume Vθr of Leading edge and Trailing edge of the blade in the radial direction;
(6) the distribution law of the circulation Vθr along the axial streamline;
(7) the stacking angle distribution of the blades.
According to the above conditions, the initial hydraulic model of the water-jet pump is designed using the 3D inverse design method. the cross-sectional shape of the impeller is shown in Figure 4, where p-1, s-1,…p17, s-17 represent the geometric profiles of pressure and suction surfaces at different spanwise.

3.2. Parametric Modeling

Before the multi-parameter optimization, the initial hydraulic model needs to be parameterized. the blade parameter optimization design carried out in this study is based on the Numeca software AutoBlade module. the parameterization process is shown in Figure 5.

3.2.1. Profile Control of Rim Hub End Wall Surface

Firstly, it is assumed that the wall surfaces of the rim and the hub are axisymmetric, and the corresponding end wall surfaces can be defined by two curves of the (Z, R) coordinate system. Where Z represents the axial direction and R represents the radial direction. the curve type is a B-spline curve with 5 control points. the reason for choosing the B-spline curve is that the curve passes all points on the curve and the second derivative of the curve is continuous.
Figure 6. B-spline curve.
Figure 6. B-spline curve.
Preprints 80601 g006

3.2.2. Definition of Flow Surface

Select the form of flow surface with linear interpolation from Hub to Shroud in the span direction as shown in Figure 7a. these flow surfaces are constructed with interpolation points from the beginning to the end of the streamline. the rule of this linear interpolation point in the (Z, R) plane satisfies the equation:
p = s p a n × p 1 + ( 1 s p a n ) × p 0
A 5-layer cross section is defined in the span direction, which are 0, 0.25, 0.5, 0.75, and 1, respectively. This definition method is the fastest way to adapt to the end wall curve update, and also applies to the definition of each type of flow surface.

3.2.3. Definition of Parameterization of Blade Stacking Line

in the field of rotating machinery, the position of the two-dimensional blade interface is usually defined on the flow surface according to the reference position (stacking point). the stacking point is first defined on the two-dimensional blade, and then positioned to the corresponding flow surface. It is described by the meridian plane and the circumferential direction based on two independent accumulation rules.
(1) Stacking point: 5 types of stacking point positions are provided in the AutoBlade meridian and tangential position modules, namely the leading-edge point (LE), the trailing edge point (TE), the center of gravity (CG), the point on the camber in the specified position (CA), and the point on the chord at the specified position (CH), as shown in Figure 7b. in the parameterization of the impeller, the stacking point is selected as the leading-edge point.
(2) Meridian position: Sweep law curve is used to define the blade section in the meridian direction. This method defines a curve on the meridian plane from the hub to the rim, which can be a straight line (Figure 8a), Bezier curve (Figure 8b), Bezier -line-Bezier combination curve (Figure 8c) and Bezier or B-Spline curve controlled by n points. in the parameterization of the impeller, the sweep mode is a straight line, and the corresponding control parameter is the sweep angle βs (SWEEP_BETA).
(3) Circumferential position: This method defines a line passing through the intersection of cross-sectional areas in the circumferential direction to describe the “bent” characteristics of the blade to determine the circumferential position of each section. Circumferential position can be defined by straight line, Bezier curve, Bezier-line-Bezier combination curve and Bezier or B-Spline curve controlled by n points. the definition above are similar to the Sweep law curve. in the parameterization of the impeller, a straight line is used to determine the circumferential position, and the control parameter is the lean angle β (LEAN_BETA).

3.2.4. Definition of the Main Blade Profile

Any blade cross-section is constructed from the mid-arc curve on the specified two-dimensional structure surface, so the structure of the blade type includes three processes, as shown in Figure 9. Firstly, the reference length and the configuration coordinate system are selected according to the pump type. Secondly, define the mid-arc. Finally, define the blade profile.
Various types of curves can be used in the definition of mid-arc, such as Simple Bezier curve, Bezier curve controlled by n points, B-spline curve controlled by n points, etc. in the parameterization of the pump impeller, Simple Bezier curve is selected. the first and last control points are associated with the leading edge and the trailing edge, respectively. the relative position between the leading and trailing edges is controlled by the stagger angle γ. the third control point is constructed by the angles β1 of the leading edge and β2 of the trailing edge. Corresponding control parameters are γ, β1 and β2.
there are two ways to define the profile of the airfoil. One is to construct the airfoil by high-order pressure and suction surface curves. the high-order Bezier curves are used for the pressure and suction surface profiles. the other is to define the thickness distribution rule by means of thickness + middle arc. in the parameters of the pump impeller, the first method is used to construct the blade profile.

3.2.5. Parametric Fitting

Fit the above parameters, including the rim, hub, meridian area stacking law, circumferential position stacking law, main blade mid-arc, main blade pressure edge and suction edge. the important features for qualified fitting results are: the shape of the runner and the shape of the front and rear edges, and the shape of the blade of the B2B section.

4. Optimization of the Design Process

in the traditional optimization design process, in order to obtain the optimal scheme, the numerical simulation of multiple schemes is usually carried out based on empirical modification of design parameter variables, which consumes a lot of computing resources and time. This article is based on the Numeca-Design3D optimization design platform, combined with CFD numerical simulation and optimization algorithms, to carry out optimization design to save design cycle and computing resources. Design of Experiments (DOE) uses Latin hypercube to construct the optimized parameter sample space, conduct numerical calculation on the constructed samples, take the performance of the design point as the optimization goal, and use artificial neural network to establish the approximate function model between the optimization goal and the optimization parameters. NSGA-Ⅱ genetic algorithm is used to optimize the approximate function model, and finally the optimized blade geometric parameters are obtained. the optimized design process is shown in Figure 10.

4.1. Optimization Parameters

in the parametric modeling in section 3.2, under normal circumstances, the flow path of the overflow component is determined. the parameters of the rim and the end wall line of the hub are unchanged. in order to control the stacking laws of the cross-sections of the blades, straight lines are used to determine the circumferential position, and the control parameter lean angle β is variable. in the parameterized optimization, the shape of the blade is mainly controlled to achieve the flow control to meet the performance requirements, so the cross-sectional control parameters of each layer of the blade are the stagger angle γ (SN (i) _CAMBER_GAMMA), the leading edge angle β1 (S (I ) _CAMBER_BETA1), the trailing edge angle β2(S(i)_CAMBER_BETA2) and the lean angle β (LEAN_BETA) have a total of 4 parameters and 16 variables. the variable name, initial value and control range of the above control parameters are shown in Table 2. Considering the influence of too large selection range of each parameter on deformed blades during blade modeling, the optimized parameter range Δγ = [-2, 2], Δβ1 = [ -3,3], Δβ2 = [-3,3], Δβ = [-1,3].

4.2. Optimization Goals and Objective Functions

the optimization goal is to improve the efficiency of the water-jet pump and improve the cavitation performance of the pump as much as possible on the premise of ensuring the pump head. Since this optimization is only for the impeller, the head of the impeller is defined as the total pressure drop at the inlet and outlet of the impeller. the cavitation performance is measured based on the lowest pressure point (P_min) [24]. Although this method has been adopted in the public literature, its accuracy and feasibility have not been widely recognized in the industry. in the multi-parameter automatic optimization of single-phase flow, there is only the lowest pressure point in the public research. therefore, in multi-parameter automatic optimization, the cavitation performance is measured based on the point with the lowest pressure, and this method will be demonstrated later in the article.
the purpose of optimal design Is to find the appropriate geometry to make the performance closest to the expected value. in multi-objective optimization, minimization is achieved by transforming objective constraints into penalty functions. the penalty function is constructed as follows:
P = W Q i m p Q s Q r e f k
in the formula, Qimp is the given value of the optimization target; Qs is the calculated value of each sample; Qref is the reference value used for dimensionless; k is the index of the penalty function; W is a weighting factor, used to increase or suppress the influence of a penalty function.
the types of penalty functions are divided into less than, equal to, and greater than. This optimization goal is applied to efficiency η, head H and the lowest pressure point P-min, see Figure 11. in terms of efficiency η, given a value of 1, a penalty function of less than the type is used. When the calculated value is less than 1, a penalty value is added to the objective function, so that the optimization parameters are far from the current value. in terms of head H, the given value is the total pressure difference between the inlet and outlet of 145000 Pa, and a penalty function equal to the type is used. When the calculated value is more than or less than the given value of 145000 Pa, the penalty value will be added to the objective function, thereby keeping the optimization parameters away from the current value. in terms of the lowest pressure point P_min, the given value is -118000 Pa, and a penalty function greater than the type is used. When the calculated value is more than the given value -118000 Pa, the penalty value will be added to the objective function, thereby keeping the optimization parameters away from the current value. the given values used above are based on the calculation results of the initial model, and correspondingly improved according to the desired target.

4.3. Samples

the design principle of Latin hypercube experiment [25,26] is to devise a scheme with uniform distribution of each parameter in space for each optimization parameter. Taking 3 parameters as an example, a sample construction is shown in Figure 12. in the construction of the sample space, the final number of samples is determined by the number of discrete factors. the number of samples is 1000.

4.4. Artificial Neural Network

An artificial neural network [27,28,29,30] is used to construct an approximate function to describe the relationship between the objective function and the optimization parameters. the neural network is generally composed of multiple layers. the first layer is the input layer, one or more hidden layers are in the middle, and the last layer is the output layer, as shown in Figure 13.
the most commonly used learning algorithm is the B-P algorithm which return the difference in reverses. the output values are calculated sequentially, and then the error (the difference between the calculated value and the correct value) is returned in reverse order to adjust the weights in the calculation model. the B-P algorithm is a relatively simple method for calculating the network performance change value by calculating a single weight change.

4.5. Genetic Algorithm

Genetic algorithm is a method to find the optimal solution by simulating the process of biological evolution. It can automatically accumulate information related to the search space during the search process and obtain the optimal relationship by adaptively controlling the search process. Non-dominated sorting algorithm NSGA-Ⅱ [31] has become one of the most popular multi-objective genetic algorithms because of its better performance in global optimization.

5. Verification of Optimization Results

According to the optimization parameters in Table 2, Latin hypercube is used to construct the optimization parameter sample space. Numerical calculations are carried on the constructed samples to obtain the optimization target head, efficiency and minimum pressure. the artificial neural network is used to establish an approximation function model between the optimization target and the optimization parameters. Topological structure is 2 layers. NSGA-Ⅱ genetic algorithm is used to optimize the approximate function model. After 10,000 iterative calculations, the optimized blade geometry parameters and performance are finally obtained.

5.1. Comparison of Blade Geometry and Performance between Original Model and Optimized Model

the geometric shapes of the blades of the original model and the optimized model in different directions are shown in Figure 14. From the figure, it can be found that compared with the original model, the optimized blade shape changes mainly at the trailing edge angle β2, the stagger angle γ and the stacking laws.

5.1.1. Numerical Simulation Method

the mesh of the impeller is shown in Figure 15. the tip clearance of the grid is set to Nomal Distance = 0.3 mm. the mesh of the tip clearance is crucial for accurately predicting the cavitation in the tip clearance. A 25 layers grid is used to divide the tip clearance region to obtain a better flow field. the number of grids of each hydraulic component is: 1.27 million grids in the impeller, 2.04 million grids in the diffuser, 450,000 grids in suction pipe, and 370,000 grids in outpipe. the total grid number is about 4.13 million. the mesh meets with the mesh-independence requirement.
the commercial software ANSYS CFX was used to calculate the internal flow the mixed flow in the pump[32,33,34]. the liquid phase was 25 °C water with a density of 997 kg/m3 and a kinematic viscosity of 8.899×10−4 kgm−1s−1. When the cavitation occurs, the vapor phase uses 25 °C water vapor, the density is 0.02308 kg/m3, and the dynamic viscosity is 9.8626×10−6 kgm−1s−1. the convergence residual is set to 10-5. the inlet boundary condition is pressure inlet and the import turbulence is set to Medium (intensity=5%). the outlet boundary condition is mass outflow. the non-slip wall is used for the solid wall. the impeller domain rotates with a speed of 1450 r/min. the blade and hub are arranged to rotate, and the shroud wall speed is set to Counter rotating wall. the interface between the rotating part and the stationary part is set to the Frozen Rotor interface. High Resolution is selected for the Advection Scheme Option, and the preset iteration step is 3000.

5.1.2. Performance Comparison by CFD

the hydraulic performances of original model and optimized model are shown in Table 3. the efficiency of the pump automatically optimized by artificial neural network and genetic algorithm reaches 89.28%, which is 1.24% higher than the efficiency of the original hydraulic model. the head of the optimized pump is 13.51 m, which is 0.39 m higher than the original hydraulic model. After the optimization, the efficiency of the impeller reached 96.33%, which is 0.35% higher than the original hydraulic model. the optimized impeller head reached 14.54 m, which is 0.272 m higher than the original hydraulic model. the optimized pump power reached 67.91 kW, an increase of 1.03 kW compared to the original hydraulic model. After optimization, the lowest pressure point on the blade reaches -291232 Pa.

5.2. Comparison of Original Model and Optimized Model Flow Field

Figure 16 shows the meridian pressure distribution of the original model and the optimized model. It can be seen from the figure that the pressure distribution of the original model and the optimized model on the meridian surface is consistent as a whole. That is, along the flow direction, the pressure gradually increases from the inlet to the outlet. Locally, in the area close to the hub on the meridian surface and before the leading edge of the blade, the area of the relatively higher pressure area existing in the original model is larger than the optimized model. This means the optimized impeller will generate greater suction near the blade inlet, which makes the impeller head higher and more powerful. It is consistent with the performance of the impeller of the original and optimized model shown in Table 3. Figure 17 shows the meridian velocity distribution of the original model and the optimized model. It can be seen from the figure that the velocity distribution on the meridian surface before and after optimization is slightly different, but not obvious.
Figure 18 shows the blade pressure distribution of the original model and the optimized model. Along the flow direction, the pressure on the suction surface of the blade increases gradually from the inlet side to the outlet side. Locally, the suction surface of the blade is close to the low-pressure area at the leading edge and rim of the blade. the area of the low-pressure area of the optimized model is larger than that of the original model, and the low-pressure area extends toward the outlet. Usually cavitation will occur and develop in the low-pressure area. From this perspective, the cavitation performance of the optimized impellerhas not been significantly improved.
the difference simply seen from the cloud is too intuitive, and a more detailed comparison of the pressure distribution in different spans is needed. Figure 19 shows the blade pressure distribution in different directions for both the original model and the optimized model. Figure 19a shows the pressure distribution on the blade at span = 0.1. the abscissa is the dimensionless position in the flow direction of the blade, and the ordinate is the static pressure on the blade. It can be seen that on the pressure surface, the pressure distribution of the optimized model is very close to that of the original model in the streamline direction of 0-0.05. in the streamline direction of 0.05-0.15, the pressure of the optimized model is higher than that of the original model, and the pressure change of the optimized model in this range is relatively slow. in the streamline direction of 0.2-0.4, the pressure of the optimized model is slightly lower than that of the original model, but the difference is not obvious, and presents a very similar upward trend along the flow direction, which keeps close to the streamline direction of 0.6. in the streamline direction of 0.6 ~ 0.9 area, the geometric difference between the optimized model and the original model has a significant impact on the pressure field of the blade. the pressure of the optimized model is significantly smaller than the original model. in the streamline direction of 0.9 ~ 1.0 area, the pressure of the optimized model is obvious. It is significantly higher than the original model. That is, the pressure of the original model drops sharply in this area, forming a larger pressure gradient, while the optimized model leaves a high pressure in this area. therefore, the pressure at the outlet of the blade keeps a higher value, and can make the blade pressure more smooth transition. On the suction surface, the pressure of the optimized model is significantly lower than that of the original model in the 0-0.3 area of the streamline direction, i.e. at the leading edge of the inlet of the blade, which means that in this area where cavitation is most likely to occur, the optimized model is more likely to have cavitation in this area than the original model. in the 0.55-0.9 area of the streamline direction, the pressure distribution of the optimized model suction surface is lower than that of the original model. the overall pressure difference between the suction side and the pressure side of the optimized model is larger than that of the original model. This shows that the optimized blade has better performance and the efficiency of the impeller is higher.
Figure 19b shows the pressure distribution of blade pressure surface and suction surface at span = 0.5 in spanwise direction, and the distribution is similar to that of blade at span = 0.1 in spanwise direction in Figure 19a. the difference is obvious that the optimized model extends backward in the low-pressure region of the suction surface of the blade with span = 0.5. Figure 18c shows the pressure distribution of blade pressure surface and suction surface at span = 0.9 in spanwise direction, which is similar to the distribution law of blade in span = 0.1 in spanwise direction in Figure 19a. the obvious difference is that the pressure distribution of the suction surface of the optimized model is higher than that of the original model in the area of 0.6-0.9 streamline direction. the pressure difference between the suction surface and the pressure surface of the optimized model is smaller than that of the original model. This shows that the work done by the blade in the optimized model in the tip region is smaller than that of the original model, which is beneficial to reduce the speed of the blade tip near the outlet and avoid the flow separation caused by a large speed difference.

5.3. Numerical Calculation and Test Verification of Optimized Model Performance

According to the optimized parameters, a three-dimensional prototype of the impeller is manufacture, as shown in Figure 20. the model pump test was carried out on standard pump performance test bench [32,33,34]
Figure 21 shows the comparison between the numerical results of the flow head curve of the optimized model and the test results. It can be seen from the comparison that the numerical results in the small flow rates are lower than the test results, while the numerical results in the large flow rates are slightly higher than the test results. in the design condition test, the measured head is 13.26 m, and the numerically calculated head is 13.51 m. the head error of numerical calculation and experiment is 1.89%. Figure 22 is the comparison between the numerical results and the experimental results of the flow efficiency of the optimized model. the efficiency measured in the test under design condition is 87.1%, while the efficiency calculated by simulation is 89.28%. the efficiency error of numerical calculation and experiment is 2.5%. From the comparison of experimental results and numerical calculations, it can be known that the optimization method is feasible and a better hydraulic model is obtained.

5.4. Cavitation Performance Verification

the basic equation of pump cavitation can be constructed[33,34] as Equation (3).
z c + p c ρ g + v c 2 2 g z K h c K p V ρ g + p V ρ g p K ρ g = K 1 v 1 2 2 g + K 2 w 1 2 2 g
Where, pv is saturated vapour pressure, Pa; pK is pressure at K point near leading edge, Pa; pc is pressure at suction fluid surface, Pa; vc is velocity at suction fluid surface, m/s; zc is the height at fluid surface, m; zK is the height at K point, m; hc-k is the hydraulic loss of c-k, m; v1 is absolute speed, m/s; w1 is relative speed, m/s.
the items in the square brackets on the left are determined by the device, representing the energy of the liquid at the pump inlet minus the energy that is rich in overcoming the head of the vaporization pressure. This is called the available Net Positive Suction Head (NPSHa). the defined equation is as follows,
Preprints 80601 i001
Preprints 80601 i002
Preprints 80601 i003
According to Equation (4), it is possible to derive the condition that the pump has a critical cavitation state, pK = pVNPSHa = NPSHr. NPSHa can be calculated by the follow Equation (7)
Preprints 80601 i004
in the process of cavitation development, the head decreases with the dropping of NPSHa until the head curve breaking, so NPSHa corresponds to different inlet total pressure and cavitation generation region and state. Usually when the head drops to 3%, the corresponding NPSHa is considered as the pump required NPSHr. This condition is defined as the critical cavitation condition. At this state, NPSHr = NPSHa = NPSHc, NPSHc is defined as critical Net Positive Suction Head.
Figure 23 displays cavitation performance comparison between calculated value and test results. It can be seen from the figure that measured head first presents an approximate horizontal state with the decrease of NPSHa of the device. When the cavitation curve of the pump approaches the first critical cavitation point or approaches the fractured cavitation state, the head rises briefly with the further decrease of NPSHa of the device, When the cavitation curve of the pump approaches the first critical cavitation point or the fractured cavitation state, with the further reduction of the NPSH of the device, the head of the pump rises briefly, and then increases with the NPSH of the device. As the volume decreases, the head begins to drop sharply, which means that the cavitation in the pump begins to have a significant impact on the performance of the pump, and it enters a state of fractured cavitation. This type of pump fracture cavitation characteristic curve belongs to E-type cavitation curve. the definition and explanation of this type are detailed in section 5.52 of pump cavitation foundation published by Pan [2]. in the book, it is pointed out that the increase of head in the first critical cavitation point area is due to the increase of head caused by cavitation in the blade. With the decrease of NPSHa, the cavitation performance curve obtained by numerical calculation keeps the approximate horizontal state of head first, without obvious fluctuation. With the further reduction of NPSHa, the head drops slowly when approaching the state of fracture cavitation, and starts to drop sharply when entering the state of critical cavitation. the fracture cavitation type of this pump belongs to B-type cavitation curve. the critical NPSHr of the pump can be determined according to the 3% drop of the head. From the cavitation performance curve in the figure, the critical NPSH measured by the test is 7.58 m, and the critical NPSH calculated by the numerical value is 7.36 m. the critical NPSHr measured by the experiment is 2.99% higher than that calculated by the numerical method.

5.5. Numerical Calculation of Cavitation Flow in the Evolution of Cavitation

the conditions of cavitation flow structure comparison between numerical calculation and experiment are shown in Figure 24. the selected operating conditions are the key operating points in the development of the cavitation process. Figure 24 shows the cavitation isosurface with a vapor volume fraction of 0.1 during the evolution of cavitation. It can be seen that at the initial stage, when NPSHa = 9.81 m, the inlet pressure of the pump is high and the pressure inside the pump is high. there is no cavitation in the impeller flow passage and the blade surface. With the decrease of the inlet pressure, when NPSHa = 8.08 m, cavitation appears in local region, where is about 0.7~1.0 span at the blade suction surface near the tip, and is about 0.3~0.6 chord at the blade tip. From the top view, there is the elongate cavitation zone adjacent to the suction surface and the tip of the blade. With the inlet pressure is further reduced, when NPSHa = 7.57 m, it can be found that cavitation has developed to a certain level, and it begins to occupy a part of the suction surface of the blade, mainly distributed in about 0.6~1.0 span, and in about 0.25~0.85 chord. Compared with NPSHa = 8.08 m, the cavitation area is larger, cavitation region extends to the leading edge, extends to the trailing edge and extends to the hub. From the view of the axial plane, the cavitation zone in the tip region occupies more chord length region, and the cavitation attached to the blade suction surface is thicker. That is, the cavitation extends toward the blade flow passage. When the cavitation develops to a critical cavitation condition, the pump performance begins to decrease. At this moment, NPSHa = 7.31 m. the area covered by the sheet cavitation has approached the trailing edge of the blade, and approaches the 0.5 span direction. But the sheet cavitation does not extend forward to leading edge region. From the view of the axial plane, the cavitation area in the tip region increases, and the cavitation is further thickened. As the inlet pressure is further reduced, when NPSHa = 7.16 m, the cavitation phenomenon is serious, and the cavitation area has covered the trailing edge of the blade, which seriously affects the work of the blade, which directly leads to a steep drop in the performance of the pump.

6. Conclusions

in order to obtain a set of high-performance water-jet pumps, this paper mainly focuses on the mixed flow impeller design of water-jet pumps based on multi-parameter and multi-objective genetic algorithm. the main work and conclusions are as follows:
(1) Regarding the water-jet pump, an initial model is designed by the parametric 3D inverse design method. then a multi-parameter and multi-objective genetic algorithm optimization design strategy for the impeller water-jet is constructed, and a water-jet pump with better performance is designed.
(2) in the parameterized optimization, the shape of the blade is mainly controlled to achieve the flow control to meet the performance requirements. the cross-section control parameters of each layer of the blade are staggered angle γ, leading edge angle β 1, trailing edge angle β 2 and lean angle β. there are 4 parameters and 16 variables in total. the Latin hypercube is used to construct the sample space of optimization parameters. Numerical calculations are carried out for the constructed samples to obtain the optimization target head, efficiency and minimum pressure point of the design point of the pump. the approximate function model between the optimization target and the optimization parameters is established by the artificial neural network. NSGA - Ⅱ genetic algorithm is used to optimize the approximate function model. Finally, the optimized blade geometry and performance are obtained.
(3) the efficiency of the optimized pump is 89.28%, which is 1.24% higher than that of the original model. the head of the optimized pump is 13.51 m, which is 0.39 m higher than that of the original model. From the comparison of experimental results and numerical calculations, it can be known that the optimization method is feasible and a better hydraulic model is obtained.

Acknowledgements

This work was funded by the China Postdoctoral Science Foundation Funded Project (Grant No.2019M651734,2023M733355), Jiangsu University Youth Talent Development Program (2020), the Chunhui Program Cooperative Scientific Research Project of the Ministry of Education, Research Project of State Key Laboratory of Mechanical System and Vibration (MSV202203), Natural Science Foundation of China (Grant No.51906085, Grant U20A20292), Jiangsu Province innovation and Entrepreneurship Doctor Project (2019), Zhejiang Postdoctoral Project (2019). This work was also financially supported by the China and Germany international Postdoctoral Exchange Fellowship Program by the office of China Postdoctoral Council and the Helmholtz Centre (grant No: 2020030). Thanks to Li Yajie (Colleague, Jiangsu University) for proof reading the English writing of this article. Thanks to Daniela Piccioni Koch (advisor, Karlsruhe insti-tute of Technology) for proof reading the writing of this article.

References

  1. Matsumoto, K.; Tanaka, H.; Ozawa, H. Optimisation of Design Parameters of Water-Jet Propulsion System. Proceedings of the FAST93. 1993.
  2. Pan, Z.; Yuan, S. Fundamentals of Cavitation in Pumps Zhenjiang: Jiangsu University Press; 2013.
  3. Zangeneh-Kazemi, M.; Dawes, W.N.; Hawthorne, W.R. Three Dimensional Flow in Radial-inflow Turbines. ASME 1988 international Gas Turbine and Aeroengine Congress1988. P. V001T01A46.
  4. Zangeneh, M. A Compressible Three-Dimensional Design Method For Radial and Mixed Flow Tuurbomachinery Blades. international Journal For Numerical Methods in Fluids. 1991, 13, 599–624. [Google Scholar] [CrossRef]
  5. Goto, A.; Zangeneh, M. Hydrodynamic Design of Pump Diffuser Using inverse Design Method and CFD. Journal of Fluids Engineering. 2002, 124, 319–28. [Google Scholar] [CrossRef]
  6. Zangeneh, M.; Goto, A.; Harada, H. On the Design Criteria For Suppression of Secondary Flows in Centrifugal and Mixed Flow Impellers. Journal of Turbomachinery. 1998, 120, 723–35. [Google Scholar] [CrossRef]
  7. Zangeneh, M.; Goto, A.; Takemura, T. Suppression of Secondary Flows in A Mixed-Flow Pump Impeller by Application of 3D inverse Design Method: Part 1—Design and Numerical Validation. ASME 1994 international Gas Turbine and Aeroengine Congress and Exposition: American Society of Mechanical Engineers; 1994. P. V001T01A14-VT01A14.
  8. Bonaiuti, D.; Zangeneh, M.; Aartojarvi, R.; Eriksson, J. Parametric Design of A Waterjet Pump by Means of inverse Design, CFD Calculations and Experimental Analyses. Journal of Fluids Engineering. 2010, 132, 031104. [Google Scholar] [CrossRef]
  9. Daneshkah, K.; Zangeneh, M. Parametric Design of A Francis Turbine Runner by Means of A Three-Dimensional inverse Design Method. IOP Conference Series: Earth and Environmental Science. 2010, 12, 012058. [Google Scholar] [CrossRef]
  10. Oyama, A.; Liou, M.S.; Obayashi, S. Transonic Axial-Flow Blade Optimization: Evolutionary Algorithms/Three-Dimensional Navier-Stokes Solver. Journal of Propulsion & Power. 2004, 20, 612–9. [Google Scholar]
  11. Sonoda, T.; Arima, T.; Olhofer, M.; Sendhoff, B.; Kost, F.; Gieß, P.A. A Study of Advanced High-Loaded Transonic Turbine Airfoils. 2006. P. 1275–83.
  12. Sonoda, T.; Yamaguchi, Y.; Olhofer, M.; Sendhoff, B.; Schreiber, H.A. Advanced High Turning Compressor Airfoils For Low Reynolds Number Condition—Part I: Design and Optimization. international Gas Turbine and Aero Engine Congress, Atlanta, Ga, June2004. P. 350-9.
  13. Dennis, B.; Egorov, I.; Han, Z.X.; Dulikravich, G.; Poloni, C. Multi-Objective Optimization of Turbomachinery Cascades For Minimum Loss, Maximum Loading, and Maximum Gap-To-Chord Ratio. international Journal of Turbo & Jet Engines. 2013, 18, 201–10. [Google Scholar]
  14. Kontoleontos, E.A.; Giannakoglou, K.C.; Koubogiannis, D.G. Robust Design of Compressor Cascade Airfoils, Using Evolutionary Algorithms and Surrogate Models. 2005.
  15. Chen, B.; Yuan, X. Advanced Aerodynamic Optimization System For Turbomachinery. Journal of Turbomachinery. 2008, 130, 111–20. [Google Scholar] [CrossRef]
  16. Öksüz Ö, Akmandor İS. Multi Objective Aerodynamic Optimization of Axial Turbine Blades Using A Novel Multi-Level Genetic Algorithm. Journal of Turbomachinery. 2008, 132, 2375–88.
  17. Zangeneh, M.; Daneshkhah, K.; Dacosta, B. A Multi-Objective Automatic Optimization Strategy For Design of Waterjet Pumps. 2008.
  18. Bonaiuti, D.; Zangeneh, M. On the Coupling of inverse Design and Optimization Techniques For the Multiobjective, Multipoint Design of Turbomachinery Blades. Journal of Turbomachinery. 2009, 131, 021014–29. [Google Scholar] [CrossRef]
  19. Papila, N.; Dorney, D.J.; Griffin, L.; Wei, S.; Papila, N.; Dorney, D.J.; et al. Shape Optimization of Supersonic Turbines Using Global Approximation Methods. Journal of Propulsion & Power. 2002, 18, 509–18. [Google Scholar]
  20. Boselli, P.; Zangeneh, M. An inverse Design Based Methodology For Rapid 3D Multi-Objective/Multidisciplinary Optimization of Axial Turbines. ASME 2011 Turbo Expo: Turbine Technical Conference and Exposition2011. P. 1459-68.
  21. Huang, R.; Luo, X.; Ji, B.; Wang, P.; Yu, A.; Zhai, Z.; et al. Multi-Objective Optimization of A Mixed-Flow Pump Impeller Using Modified NSGA-II Algorithm. Science China Technological Sciences. 2015, 58, 2122–30. [Google Scholar] [CrossRef]
  22. Bonaiuti, D.; Arnone, A.; Corradini, U.; Bernacca, M. Aerodynamic Redesign of A Mixed-Flow Pump Stage. 2003.
  23. Yang, Q.; Wang, Y. Principle and Application of Low Noise Pump Jet Design: HUST Press; 2016.
  24. Cao, L.L.; Che, B.X.; Hu, L.J.; Wu, D.Z. Design Method of Water-Jet Pump Towards High Cavitation Performances. 2016, 129, 012067. [Google Scholar] [CrossRef]
  25. Shen, Y.; Jia, M.; Zhu, L. Noise Analysis of A Centrifugal Leaf Vacuum Based On Latin Hypercube Sampling. Journal of Vibration and Shock. 2016, 35, 93–7. [Google Scholar]
  26. Wang, W.; Pei, J.; Yuan, S. Optimization of Impeller Meridional Shape Based On Radial Basis Neural Network. Transactions of the Chinese Society of Agricultural Machinery. 2015, 46, 78–83. [Google Scholar]
  27. Fonseca, C.M.; Fleming, P.J. Genetic Algorithms For Multiobjective Optimization: Formulationdiscussion and Generalization. international Conference On Genetic Algorithms1993. P. 416-23.
  28. Horn, J.; Nafpliotis, N.; Goldberg, D.E. A Niched Pareto Genetic Algorithm For Multiobjective Optimization. Evolutionary Computation, 1994 IEEE World Congress On Computational intelligence, Proceedings of the First IEEE Conference On2002. P. 82-7 Vol.1.
  29. Schaffer, J.D. Multiple Objective Optimization With Vector Evaluated Genetic Algorithms. international Conference On Genetic Algorithms1985. P. 93-100.
  30. Srinivas, N.; Deb, K. Muiltiobjective Optimization Using Nondominated Sorting in Genetic Algorithms: MIT Press; 1994.
  31. Deb, K.; Pratap, A.; Agarwal, S.; Meyarivan, T. A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II. IEEE Transactions On Evolutionary Computation. 2002, 6, 182–97. [Google Scholar] [CrossRef]
  32. Long, Y.; An, C.; Zhu, R.; Chen, J. Research On Hydrodynamics of High Velocity Regions in A Water-Jet Pump Based On Experimental and Numerical Calculations At Different Cavitation Conditions, Physics of Fluids 2021, 33, 045124.
  33. Long; Yun; Zhang, Y.; Chen, J.; Zhu, R.; Wang, D. A Cavitation Performance Prediction Method For Pumps PART2-Sensitivity and Accuracy. Nuclear Engineering and Technology (2021).
  34. Yun, L.; Rongsheng, Z.; Dezhong, W. A Cavitation Performance Prediction Method For Pumps PART1-Proposal and Feasibility. Nuclear Engineering and Technology, 2020.
Figure 2. Water-jet mixed flow pump hydraulic components.
Figure 2. Water-jet mixed flow pump hydraulic components.
Preprints 80601 g002
Figure 3. Parameters of flow passage.
Figure 3. Parameters of flow passage.
Preprints 80601 g003
Figure 4. 3D blade shape of water-jet pump initial hydraulic model.
Figure 4. 3D blade shape of water-jet pump initial hydraulic model.
Preprints 80601 g004
Figure 5. Blade parameterization process.
Figure 5. Blade parameterization process.
Preprints 80601 g005
Figure 7. Stream surface traces definition.
Figure 7. Stream surface traces definition.
Preprints 80601 g007
Figure 8. Available sweep law.
Figure 8. Available sweep law.
Preprints 80601 g008
Figure 9. Definition of the camber curve.
Figure 9. Definition of the camber curve.
Preprints 80601 g009
Figure 10. Technical route for optimized design.
Figure 10. Technical route for optimized design.
Preprints 80601 g010
Figure 11. Penalty objective function.
Figure 11. Penalty objective function.
Preprints 80601 g011
Figure 12. DOE for Latin Hypercube of 3 parameters.
Figure 12. DOE for Latin Hypercube of 3 parameters.
Preprints 80601 g012
Figure 13. Artificial Neural Network.
Figure 13. Artificial Neural Network.
Preprints 80601 g013
Figure 14. Blade geometry of the original impeller and the optimized impeller.
Figure 14. Blade geometry of the original impeller and the optimized impeller.
Preprints 80601 g014
Figure 15. Impeller and diffuser mesh [32].
Figure 15. Impeller and diffuser mesh [32].
Preprints 80601 g015
Figure 16. the pressure distribution of Meridian plane.
Figure 16. the pressure distribution of Meridian plane.
Preprints 80601 g016
Figure 17. the velocity distribution of Meridian plane.
Figure 17. the velocity distribution of Meridian plane.
Preprints 80601 g017
Figure 18. the pressure distribution on the blade.
Figure 18. the pressure distribution on the blade.
Preprints 80601 g018
Figure 19. the pressure distribution on different span.
Figure 19. the pressure distribution on different span.
Preprints 80601 g019
Figure 20. Real optimized impeller.
Figure 20. Real optimized impeller.
Preprints 80601 g020
Figure 21. the flow-head curve of CFD and experiment.
Figure 21. the flow-head curve of CFD and experiment.
Preprints 80601 g021
Figure 22. the flow-efficiency curve of CFD and experiment.
Figure 22. the flow-efficiency curve of CFD and experiment.
Preprints 80601 g022
Figure 23. Cavitation performance curves of CFD and experiment.
Figure 23. Cavitation performance curves of CFD and experiment.
Preprints 80601 g023
Figure 24. Cavitation volume fraction αv = 0.1 isosurface in the process of cavitation.
Figure 24. Cavitation volume fraction αv = 0.1 isosurface in the process of cavitation.
Preprints 80601 g024
Table 1. Design Parameters of Water-jet Pump.
Table 1. Design Parameters of Water-jet Pump.
Parameter Variables
Flow rate Q = 0.46 m3/s
Head H ≥ 13 m
Rotating speed n = 1450 rpm
Table 2. Variables, original value and boundary of blade parameters.
Table 2. Variables, original value and boundary of blade parameters.
Parameters Variables Lower bound Original value Upper bound
Γ S1_CAMBER_GAMMA 54.3856 56.3856 58.3856
S2_CAMBER_GAMMA 57.8627 59.8627 62.8627
S3_CAMBER_GAMMA 60.4048 62.4048 65.4048
S4_CAMBER_GAMMA 62.4439 64.4439 66.4439
S5_CAMBER_GAMMA 64.0463 66.0463 68.0463
β1 S1_CAMBER_BETA1 55.0341 58.0341 61.0341
S2_CAMBER_BETA1 58.8031 61.8031 64.8031
S3_CAMBER_BETA1 62.1567 65.1567 68.1567
S4_CAMBER_BETA1 61.4439 64.4439 67.4439
S5_CAMBER_BETA1 65.7197 68.7197 71.7197
β2 S1_CAMBER_BETA2 45.5651 48.5651 51.5651
S2_CAMBER_BETA2 50.6644 53.6644 56.66436
S3_CAMBER_BETA2 53.3824 56.3824 59.3824
S4_CAMBER_BETA2 57.3875 60.3875 63.3875
S5_CAMBER_BETA2 55.6307 58.6307 61.6307
β LEAN_BETA 1.5408 2.5408 4.5408
Table 3. Comparison of the hydraulic performance between with Original and Optimized model.
Table 3. Comparison of the hydraulic performance between with Original and Optimized model.
Parameter Pump Head H Pump Efficiency η Impeller Head Hi Impeller Efficiency ηi Pump Power P
unit M % m % kW
Original model 13.12 88.04 14.27 95.98 66.88
Optimized model 13.51 89.28 14.54 96.33 67.91
Disclaimer/Publisher’s Note: the statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

© 2024 MDPI (Basel, Switzerland) unless otherwise stated