Design and optimization of bionic Nautilus volute for a hydrodynamic retarder

Abstract

To reduce the energy loss of volute and improve the emergency braking performance of hydrodynamic retarder under design condition, a bionic volute optimization design method is proposed based on the Nautilus shell. First, the accuracy of the numerical method is verified by grid convergence analysis and experiments. Thereafter, the bionic parametric expression is performed based on the cross-section shape of the Nautilus shell. Subsequently, sample points are created using the Optimal Latin Hypercube sampling method and numerically simulated. Radial basis function neural network and multi-objective genetic algorithm are used for optimization. Finally, the flow field inside the volute and the braking performance of the hydrodynamic retarder are compared by numerical calculation and experiment. Results show that the internal flow field of the bionic volute is more uniform, the energy loss is reduced by 78.95%, and the oil flow speed is increased by 40.61%. In addition, the stable flow field can be established faster during braking, the peak torque can be increased by 6.15%, and the onset time can be advanced by 6.58%. The analysis of energy characteristics shows that the improved performance of the bionic volute can be attributed to its improved internal oil flow conditions, thus reducing the energy loss.

KEYWORDS:

• Hydrodynamic retarder
• bionic design
• neural network
• multi-objective optimization
• entropy production theory

1. Introduction

When heavy vehicles are driving on closely spaced ramps or long and steep downhill roads, traditional mechanical braking devices often experience excessive braking power and thermal load, making the braking devices prone to thermoelastic instability and torque decline, and can easily lead to major accidents (Chunbao et al., 2015; Le Gigan et al., 2015). As an auxiliary braking device, a hydrodynamic retarder converts input kinetic energy into liquid internal energy through the interaction of blades and hydrodynamic transmission oil, and dissipates the generated thermal energy to achieve the function of retarding braking. Compared with other auxiliary braking devices, hydrodynamic retarders are widely used in heavy vehicles due to their advantages of high power density, strong heat dissipation capacity and small volume and mass (Mu et al., 2016; Mu et al., 2017; Wang, Ju, et al., 2021).

The onset time of the emergency braking torque under the design condition is one of the important indices that measures the braking performance of a hydrodynamic retarder. Previous studies of the dynamic braking performance of hydrodynamic retarders have primarily focused on the control valve system which fills and drains liquid into the working chamber. Wei et al. (2020) established hydrodynamic retarder calculation models with or without filling and drain valves and compared torque changes in the braking process of the two models. The filling valve exhibits a marked hysteresis effect on the transient braking torque. Mu et al. (2019) established an integrated model of a hydrodynamic retarder open wheel chamber and hydraulic control system, and studied the influence of hydraulic control system parameters on the braking response characteristics of the hydrodynamic retarder. The control pressure of the proportional pressure-reducing valve at the outlet of the hydrodynamic retarder and the length of the tubing have a strong influence on braking performance. Although existing researchers have conducted a large number of detailed studies on the volute of centrifugal fans, centrifugal pumps and other equipment, there are still few studies on the inlet volute of the hydrodynamic retarder. Because the working principle of hydrodynamic retarder is different from the above research objects, the flow direction of the working medium in hydrodynamic retarder and centrifugal pump volute is completely opposite. Therefore, the existing conclusions on volute design and optimization are not directly applicable to guide the design of hydrodynamic retarder inlet volute. As a portion of the hydrodynamic retarder filling system, the inlet volute also has an important influence on its dynamic braking performance. In the actual use of the existing volute, there are some problems, such as oil accumulation and high internal energy loss, which are not conducive to the rapid flow of oil into the working chamber, thus, it is necessary to optimize the volute structure design.

Traditional volute design typically uses a one-dimensional design method based on equal circulation theory, which considers that the flow parameters in volutes are uniformly distributed along the inlet circumference (Bian et al., 2021). Yu et al. (2019) proposed a volute profile based on this issue. The volute has a good diversion effect and can provide a uniform flow field for the turbine air classifier, improving its classification performance. However, experimental results showed that the flow field parameters at the inlet of the volute are different from the assumption of a uniform distribution (Ayder et al., 1993; Elholm et al., 1992). To provide a more accurate design basis for volute profile design, many researchers have conducted more in-depth studies. Chen (Chen 1996) designed a volute in two dimensions by solving a set of extended two-dimensional potential flow equations and streamline equations. Owarish et al. (1992) divided the volute into several segments and designed the shape of the two-dimensional volute by considering the influence of the cross-section shape on the performance of the volute.

With the development of computational fluid dynamics (CFD), a three-dimensional design method that can accurately predict the flow in volutes has gradually become commonplace. Simpson et al. (2009) verified the feasibility of CFD technology in volute research by comparing simulation results with experimental data. Zhang et al. (2021) studied the timing effect of outlet RGV on the hydrodynamic characteristics of centrifugal pump, determined the best installation angle between RGV and volute tongue, and obtained the mechanism of the influence of RGV timing position on pump hydrodynamic performance, which has valuable guiding significance for selecting the installation position of outlet RGV. Because the design of experiment (DOE) method can use a small amount of function evaluation in optimization problems to obtain relatively good solutions, many studies have completed performance optimizations of rotating machinery parts based on this technology (Gan et al., 2022). Du et al. (2019) used the orthogonal design test method to optimize the pear-shaped volute of the supercritical carbon dioxide (SCO₂) centripetal compressor and verified that the centrifugal compressor met the design requirements through numerical simulation. Selvaraj et al. (2020) used a Taguchi L9 orthogonal array to design the blade inlet angle, impeller diameter and impeller width, and determined the optimal parameter combination according to the results of grey relational analysis (GRA), effectively optimizing a radial tip centrifugal blower. Yun et al. (2016) optimized the diffuser parameters of a reactor coolant pump by an L18 (3⁷) orthogonal table, which effectively improved the efficiency and head of the pump under design conditions. Yang et al. (2021) combined numerical simulation and the Taguchi method, optimized the linear parameters of the hub and shield in the meridian of the diffuser of an electric submersible pump (ESP), and improved the head and shaft power under the design flow condition.

However, these methods still have the problems of high computation and time costs, which limits their improvement of optimization efficiency. With the wide application of artificial intelligence technology, the combination of CFD technology and machine learning has become a new development trend of parameter optimization methods (Zhao et al., 2022). First, the sample library is extracted from the design domain according to the DOE method, and then, the surrogate model is constructed according to the numerical simulation results of sample points. Finally, parameter optimization is performed on the surrogate model using the optimization algorithm. Currently, this optimization method using a surrogate model has been widely used to accelerate the design process of volutes (Shim et al., 2016). Bian et al. (2021) analysed the sensitivity of the geometric parameters of volute cross sections based on a Gaussian process surrogate model, and realized an effective optimization process through Bayesian optimization (BO), which made the flow field inside the volute more uniform, significantly reduced the flow loss in the volute and reduced the flow distortion at the outlet. Zhao et al. (2022) optimized nine design parameters of impeller and volute by using flow loss visualization technology based on entropy generation theory and machine learning method, which broadened the effective working area of multistage double-suction centrifugal pump and improved its energy efficiency. Acarer et al. (2020) parameterized an existing radial wind turbine (RWT) with 14 parameters, constructed an artificial neural network (ANN) surrogate model using CFD calculation results, and achieved structural optimization by particle swarm optimization (PSO). Compared with the existing RWT, peak power coefficient increased by 103%. Heinrich (Heinrich et al., 2016) used a cubic B-spline curve to express the cross-sectional shape of volutes, and built a genetic optimization numerical model for centrifugal compressor volutes of light commercial vehicles with control points as design variables, which expanded the working range and improved the total pressure ratio and isentropic efficiency. To improve the efficiency of rectangular cross-section centrifugal compressor volutes throughout the whole working area, Huang et al. (2015) selected eight parameters of volute cross-section shape and radial position as design variables, constructed a response surface approximate model (RSM), and optimized the parameters. After optimization, the efficiency of the volute increased by 7%, and the overall efficiency of the centrifugal compressor also increased by 0.4%. In addition, in order to solve the dimension curse problem in the optimization process, A. Serani (Serani et al., 2023) proposed a parameter model embedding method, which effectively improved the effectiveness and efficiency of optimization by reducing the dimension of design space, and verified its significant industrial significance through cases. Although these studies have discussed the optimization design of hydraulic components to a certain extent, most of these designs are limited by the experience of designers, and only some structural parameters in the existing configuration are selected for optimization. As shown in Figure 1, no innovative hydraulic structure is proposed.

Figure 1. Volute structures in the optimization design process of this study.

In recent years, bionic technology has been widely developed and used in hydraulic machinery design through the simulation of biological structures and evolutionary processes to encourage pioneering ideas and inspiration (Wu et al., 2018). Gu et al. (2018) established a vortex pump model with a bionic groove at the volute tongue according to the drag reduction and vibration reduction mechanism of the groove structure on the abdomen surface of mantis shrimp, which effectively reduced the pressure pulsation around the outlet volute tongue. Inspired by the aerodynamic characteristics of an owl wing's leading edge, Wang et al. (2019, 2021) proposed a bionic design of the volute tongue and optimized it by using a multi-objective efficient global optimization (EGO) algorithm, which reduced the flow loss in the volute and alleviated the pressure abrupt change near the volute tongue. Dong and Wu also performed bionic designs of volute tongues based on the owl wing leading edge (Dong & Dou, 2021; Wu et al., 2020). The former improved the aerodynamic performance of the centrifugal fan under different flow conditions, while the latter effectively reduced the aerodynamic noise of the centrifugal fan. Throughout the evolutionary history of Nautilus over hundreds of millions of years, it has formed a spiral structure that can effectively reduce flow resistance. Li et al. (2022) designed a new type of vertical axis wind turbine with a blade structure that is similar to the shape of the Nautilus shell, achieving a strong wind collection capacity and wind energy utilization rate. Gao et al. (2020) explored the influence of the Nautilus blade number on the torque characteristics of wind turbines and further improved their average torque coefficient. Li et al. (2023) proposed a Nautilus bionic flow channel for proton exchange membrane fuel cells (PEMFC). Compared with other types of flow channels, the reactant distribution in the Nautilus bionic flow channel is more uniform, the concentration range loss is lower, and the peak current density and peak power density can be markedly improved. Lv et al. (2022) proposed a type of Nautilus volute for a multi-blade centrifugal fan, which can improve its aerodynamic performance and effectively reduce noise. In addition, bionic technology is also used in the optimization design process of blades and other hydraulic machinery parts (Dai et al., 2021; Gu et al., 2016; Huang et al., 2021; Li et al., 2019; Ma et al., 2022).

Based on the above research methods, the bionic Nautilus volute of a hydrodynamic retarder is studied and optimized in this paper. First, a bionic volute scheme is proposed and parameterized according to its structural characteristics. Then, using the optimal Latin hypercube sampling and CFD simulation results, the RBF neural network surrogate model is established, and the parameter correlation is analysed. Based on the surrogate model, the NSGA-II algorithm is used to complete parameter optimization. Finally, the internal flow characteristics of the volute under design conditions and the transient torque characteristics of the hydrodynamic retarder during emergency braking are compared. In addition, flow loss visualization technology based on entropy production theory is used to analyse the energy characteristics of the two volutes in detail, and the optimization mechanism is explored.

2. Numerical methodology and verification

2.1. Description of baseline model

The baseline model of this study is a single-cycle circular parallel hydrodynamic retarder equipped with a commercial heavy truck, whose inlet flow is 200 L/min and rotational speed is 2000 r/min under design conditions. The hydrodynamic retarder is introduced using the horizontal division structure shown in Figure 2. The hydrodynamic retarder is primarily composed of a shell, a rotor and a stator. The inlet volute and the outlet ring are processed on the shell.

Figure 2. Cross-sectional view of the hydrodynamic retarder.

The inlet volute and the outer ring are connected with the filling valve and the drain valve, respectively, forming the filling and discharging channels of the hydrodynamic retarder. The stator and the shell are fixed to the car body, and the rotor is connected with the output shaft of the gearbox through gears. When the vehicle is running normally, both the filling and drain valves are closed. At this time, the hydrodynamic retarder is filled with gas and does not provide braking torque. After the driver depresses the brake pedal, the filling valve is fully opened quickly, and the hydrodynamic transmission oil enters the working chamber composed of rotor and stator through the inlet volute. The output shaft of the gearbox drives the rotor to rotate, which drives the oil to circulate in the working chamber. The oil impinges the stator blades to generate reaction force, thus resulting in braking torque, and then, the oil flows back to the oil tank through the outer ring and the drain valve. By adjusting the opening of the drain valve, the filling rate in the working chamber can be changed to realize the adjustment of braking torque.

2.2. Numerical setup

The original flow channel model of the hydrodynamic retarder with an inlet and outlet established by the three-dimensional modelling software is shown in Figure 3. The calculation region is divided into four parts: rotor, stator, inlet volute and outer ring.

Figure 3. Flow channel model of the hydrodynamic retarder with an inlet and outlet.

The CFD calculation in this study is carried out on the platform of Ansys Fluent. The pressure-based solver is used to solve the Navier-Stokes equation. The implicit coupling algorithm is selected to improve the convergence speed of the simulation, and the relaxation factors keep the default values unchanged. At the same time, in order to make the simulation results more accurate, the second-order upwind scheme is used to ensure sufficient calculation accuracy. When the residual error is less than 10e-5, and the monitored physical quantity reaches a constant value, the simulation is considered to converge, and the maximum number of iteration steps is 2000. For the turbulence model, this study uses the k-omega SST model, which uses the k-omega model in the near-wall region and the k-epsilon model in the far-wall region. The smooth transition between the two models is achieved by the mixing function, which combines the advantages of the two models well (Le Hocine et al., 2021).

The internal flow of the hydrodynamic retarder is basically gas–liquid two-phase flow, thus, the VOF model is selected to simulate air and oil, which are two insoluble fluids. The location of the free surface is determined according to the volume fraction of each phase fluid in the grid, and the location of the free surface is tracked. Because the influence of temperature is not considered in this study, heat transfer and temperature rise of fluid are neglected, thus, only the continuity equation (Equation 1) and momentum equation (Equation 2) are involved: (1) $\begin{array}{rl}& \frac{\mathrm{\partial }{\rho }_m}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\rho }_m{v}_m\right)=0\end{array}$(1) (2) $\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\rho }_m{v}_m\right)+\mathrm{\nabla }\cdot \left({\rho }_m{v}_m{v}_m\right)\\ & =-\mathrm{\nabla }p+\mathrm{\nabla }\cdot \left[{\mu }_m\left(\mathrm{\nabla }{v}_m+\mathrm{\nabla }{v_m}^T\right)\right]+{\rho }_{m}g\\ & +F+\mathrm{\nabla }\cdot \left(\sum _{k=1}^n{\alpha }_k{\rho }_k{v}_{dr,k}{v}_{dr,k}\right)\end{array}$(2) where ρ_m is the mixing density; v_m is the average velocity; t is the time; p is the pressure; μ_m is the mixing viscosity; g is the gravitational acceleration; F is the volume force; α_k is the volume fraction of the kth phase; ρ_k is the density of the kth phase; and v_dr,k is the drift velocity of the kth phase.

Regarding boundary conditions, the inlet is set as the velocity inlet, and the velocity of the inlet oil is determined according to the flow conditions and geometric parameters. Because the hydrodynamic retarder outlet is generally directly connected to the radiator and oil tank, it can be considered as a pressure outlet with a gauge pressure of 0. Smooth and nonslip conditions are applied to the remaining surfaces (i.e. the fluid velocity at the wall is 0). Specific simulation parameters are shown in Table 1.

Table 1. Simulation parameters.

Name	Symbol	Value (unit)
Inlet oil flow	q	200 (L/min)
Inter velocity	v	0.82 (m/s)
Outlet pressure	p	0 (bar)
Rotor speed	n	2000 (rpm)
Density	ρ	820 (kg/m^3)
Dynamic viscosity	μ	0.28 (Pa·s)
Temperature	T	80 (°C)

Because the hydrodynamic retarder is a rotating machine, it is necessary to consider the influence of the rotor on the internal working medium during this study. According to the principle of the hydrodynamic retarder, the inlet volute, stator and outer ring are set in the static domain, and the rotor is set in the rotating domain. For steady-state simulation, the frozen rotor method with low calculation cost and sufficient accuracy is used to solve the governing equations in the static coordinate system and rotating coordinate system of the static domain and rotating domain, respectively, which is one of the most commonly used methods in steady-state simulation of rotating machinery currently (Ale Martos et al., 2021; Arani et al., 2019; Tao et al., 2022). For transient simulation, because the relative position of the rotor changes with time, the dynamic grid is used to simulate the rotating domain. The force between the oil and the rotor blades is calculated according to the dynamic fluid body interaction (DFBI), to control the motion of the rotating domain grid, and the braking process simulation of the rotor with an initial speed of 2000 r/min is realized. The information of flow field parameters between calculation domains is transmitted through the interface between each other.

Because the steady-state simulation is independent of time, only the transient simulation needs to consider the influence of time. To consider the convergence and efficiency of transient simulation, it is necessary to select an appropriate time step. In this study, the time step is limited so that the rotating angle of the rotor in each time step shall not exceed the included angle between the two blades of the rotor. Thus, the expression of the time step is: (3) $\Delta t\le \frac{60}{Nn}$(3) where N is the number of rotor blades and n is the rotating speed of the rotor. Finally, the time step is set equal to Δt = 0.001 s.

Some assumptions must be made before numerical simulation, which are applicable to both steady-state and transient simulation. The influence of temperature on braking torque of hydrodynamic retarder is mainly reflected in the change of oil properties. The oil density and dynamic viscosity decrease linearly and exponentially with the increase of temperature, respectively (Chen et al., 2022; Wang et al., 2023). Studies have shown that the decrease of dynamic viscosity will lead to the increase of oil shear rate at the same rotor speed, which will aggravate the impact of oil on blades, and finally show the increase of braking torque (Bu et al., 2018). However, when the oil temperature exceeds 80°C, the dynamic viscosity basically does not change with the change of temperature. Therefore, under this working condition, the further increase of temperature has little effect on braking torque. In order to avoid the influence of oil temperature change on braking torque, the oil will be heated to 80°C before the test. On the other hand, the liquid filling speed in the rapid liquid filling test is fast, and the whole test process only lasts about 5 s (Wei et al., 2020), so the influence of short-time emergency braking test on temperature rise can be basically ignored. Wang, Ju, et al. (2021) measured by experiments that the temperature rise rate of the hydrodynamic retarder is only 0.31°C/s under the braking conditions of rotor speed 1000 r/min and braking torque 2500 N m. Zheng et al. (2016) obtained through simulation calculation that when the vehicle is braked urgently at the initial speed of 60 km/h, the maximum temperature rise in the whole braking process is 10°C. Therefore, the CFD simulation process can be simplified as an isothermal process with an oil temperature of 80°C, and the oil is considered as a viscous incompressible continuous fluid with constant density and dynamic viscosity. This assumption is also applied to other studies on braking torque characteristics of hydrodynamic retarder (Kong et al., 2020; Wang et al., 2022; Chunbao et al., 2015). It is also assumed that there is no oil leakage when working, and each part is an absolute rigid body whose surface will not be deformed.

2.3. Mesh generation

An unstructured grid is more suitable for this study because it can better describe the complex structure of the rotor and the stator of the hydrodynamic retarder. To reduce the calculation time, a polyhedral grid with higher computational efficiency is selected for grid division. In addition, mesh encryption is also performed on the complex structure area and blade surface to capture the flow field information at the corresponding position more accurately. To simulate the boundary layer and flow separation more accurately, a 6-layer prismatic grid with a total height of 1 mm and a growth rate of 1.2 is generated on each surface to ensure that the y⁺ value is less than 1. Therefore, the model has sufficient wall resolution to meet the requirements of the k-omega SST model. The interface between each region is imprinted with a grid, which ensures that the grid at the junction uses common nodes, thus transmitting the information between the two computational domains more accurately.

2.4. Grid convergence analysis

Generally, the density of the grid will affect both the simulation accuracy and calculation cost. With increasing grid density, the CFD calculation results gradually approach the real solution, but the corresponding simulation time will be longer. Therefore, it is necessary to conduct grid convergence analysis and determine the appropriate grid density by considering the simulation accuracy and time cost. Grid density is primarily controlled by grid size. To describe the grid convergence of the model, five grids with different densities are divided according to grid size (see Table 2).

Table 2. Number of meshes with different sizes.

	Number of grids ×10^5
Size	Inlet volute	Rotor	Stator	Outer ring
1 mm	17.93	73.93	58.19	16.15
2 mm	3.85	34.67	23.04	3.85
3 mm	1.83	28.30	18.33	1.93
4 mm	1.15	23.64	15.86	1.34
5 mm	0.90	21.21	13.24	1.10

As one of the important evaluation indices of the hydrodynamic retarder, braking torque is also the key parameter to measure the accuracy of numerical simulation results. The steady-state simulation of the full liquid-filled state is performed for grids with different densities under the design condition to obtain the corresponding braking torque and time cost. The braking torque change rate of grids with different densities is defined as follows: (4) $ϵ\left(T\right)=\frac{T\left(n\right)-T\left(n-1\right)}{T\left(n\right)}\times 100\mathrm{\%}$(4) where n is the grid size and T(n) is the braking torque corresponding to the current grid size. Generally, if the braking torque rate of change across two adjacent grids with different densities is within 1%, the grid density is considered to have little influence on the calculation results.

The result of grid convergence analysis is shown in Figure 4. The braking torque gradually converges with decreasing grid size, and the braking torque change rate is less than 1% under a grid density of n = 3 mm. Further reduction of the grid size will not cause marked fluctuations in the braking torque. At this time, the calculation cost is low (only 30% of the maximum calculation cost, n = 1 mm). To satisfy the calculation accuracy and efficiency requirements, the grid size is set to 3 mm, and the final grid model is shown in Figure 5.

Figure 4. Results of grid convergence analysis.

Figure 5. Hydrodynamic retarder mesh model (a) overall, (b) local.

Table 3 shows the distribution of y⁺ values under the current grid density. Although there are some regions with high y⁺ values in each computational domain, the overall y⁺ values are approximately 1, which meets the requirements of the k-omega SST model. In this case, the grid can be used to predict the braking performance of the hydrodynamic retarder and to describe the internal flow situation. A detailed distribution of y⁺ values in some regions is shown in Figure 6.

Figure 6. Distribution of y+ values near the wall of inlet volute and rotor blades.

Table 3. Distribution of y+ value.

	Inlet volute	Rotor	Stator	Outer ring
Max y^+	4.669	9.852	4.404	2.794
Ave y^+	0.172	1.121	1.385	0.453

2.5. Experimental validation

2.5.1. Test rig set up

To verify the accuracy of the numerical simulation, experiments with the baseline model are performed using the test rig shown in Figure 7. The test rig is primarily composed of a motor, inertia device, torque sensor and hydrodynamic retarder. After keeping the motor speed unchanged, and the hydrodynamic retarder running stably during the test, the motor speed signal and torque sensor signal are transmitted to the computer for analysis using a data acquisition card, and the sampling frequency is 50 Hz. The test is repeated three times for each test condition, and the average is reported as the final data.

Figure 7. Introduction of the test rig.

2.5.2 Performance comparison

The braking torque under full-filled condition refers to the maximum braking torque produced when the working chamber of the hydrodynamic retarder is filled with oil, which is the most intuitive parameter to characterize the braking ability of the hydrodynamic retarder. Restricted by the test rig, the full-filled braking characteristics of the hydrodynamic retarder in the range of 200–700 r/min are obtained in this test. The comparison between the simulation results and test data is shown in Figure 8.

Figure 8. Comparison between simulation results and test data.

The braking torque of the hydrodynamic retarder is affected by the working medium parameters, working condition parameters and design parameters. The specific expression is as follows: (5) $T=\lambda \rho g{n}^2{D}^5$(5) where λ is the torque coefficient; ρ is the density of the working medium; g is the acceleration of gravity; n is the speed of the rotor; and D is the effective diameter of the circulation circle.

Figure 8 shows that the calculated braking torque and the test data show a high degree of consistency throughout the speed range and have a nearly quadratic relationship with the rotating speed of the rotor. Considering that the working chamber of the hydrodynamic retarder cannot be filled with oil during the test, which will affect the torque coefficient, and the working parameters of oil will also change marginally due to the influence of temperature, there is a certain deviation between the simulation results and the test data. Error analysis shows that the average error between them is 7.94%, which is acceptable, and the numerical model is thus considered to be able to be used to study the characteristics of hydrodynamic retarder.

3. Entropy production theory

Due to its own viscosity and Reynolds stress in the process of movement, part of the mechanical energy of the oil flowing in volute is converted into internal energy, which leads to irreversible energy loss. Concurrently, the volute is affected by the dynamic effect of the rotor, the internal flow state is relatively complex, and the mutual impact between the oil will also increase the irreversible loss. The pressure difference between the inlet and outlet, which is commonly used to measure the energy loss of static hydraulic components, can only be analysed numerically. To further explore the energy loss mechanism inside the volute and provide guidance for subsequent optimization, entropy production theory based on the second law of thermodynamics can be used to visualize the specific distribution of flow loss in volutes.

Entropy is a state variable of the system, which is used to characterize the degree of chaos inside the system. The dissipation of mechanical energy consumption caused by irreversible factors will lead to an increase in entropy, which is called entropy production. The entropy production transport equation is as follows (Ji et al., 2020): (6) $\begin{array}{rl}& \rho \left(\frac{\mathrm{\partial }s}{\mathrm{\partial }t}+u\frac{\mathrm{\partial }s}{\mathrm{\partial }x}+v\frac{\mathrm{\partial }s}{\mathrm{\partial }y}+w\frac{\mathrm{\partial }s}{\mathrm{\partial }z}\right)\\ & =\text{div}\left(\frac{\overrightarrow{q}}{T}\right)+\frac{\mathrm{\Phi }}{T}+\frac{{\mathrm{\Phi }}_{\theta }}{T^2}\end{array}$(6) where ρ is density; s is the specific entropy; u, v and w are the velocity components in the x, y and z directions, respectively; T is the temperature; $\text{div}\left(\overrightarrow{q}/T\right)$ is the reversible heat transfer; $\mathrm{\Phi }/T$ is the entropy produced by viscous dissipation; and ${\mathrm{\Phi }}_{\theta }/T^2$ is the entropy produced by heat transfer. Based on the assumption of constant temperature in this study of a hydrodynamic retarder, $\text{div}\left(\overrightarrow{q}/T\right)$ and ${\mathrm{\Phi }}_{\theta }/T^2$ are set to 0, and $\mathrm{\Phi }/T$ is the only source of entropy production in the volute. Due to the Reynolds time averaging method, the time averaging expression of Equation 6(6) $\begin{array}{rl}& \rho \left(\frac{\mathrm{\partial }s}{\mathrm{\partial }t}+u\frac{\mathrm{\partial }s}{\mathrm{\partial }x}+v\frac{\mathrm{\partial }s}{\mathrm{\partial }y}+w\frac{\mathrm{\partial }s}{\mathrm{\partial }z}\right)\\ & =\text{div}\left(\frac{\overrightarrow{q}}{T}\right)+\frac{\mathrm{\Phi }}{T}+\frac{{\mathrm{\Phi }}_{\theta }}{T^2}\end{array}$(6) is: (7) $\rho \left(\frac{\mathrm{\partial }\overline{s}}{\mathrm{\partial }t}+\overline{u}\frac{\mathrm{\partial }\overline{s}}{\mathrm{\partial }x}+\overline{v}\frac{\mathrm{\partial }\overline{s}}{\mathrm{\partial }y}+\overline{w}\frac{\mathrm{\partial }\overline{s}}{\mathrm{\partial }z}\right)=\frac{\overline{\mathrm{\Phi }}}{T}=S$(7) where S is the local entropy production rate, which consists of the direct entropy production rate S_VD caused by the time-averaged flow field and the turbulent entropy production rate S_TD caused by velocity fluctuation: (8) $\begin{array}{rl}{S}_{VD}& =\frac{\mu }{T}\{2\left[{\left(\frac{\mathrm{\partial }\overline{u}}{\mathrm{\partial }x}\right)}^2+{\left(\frac{\mathrm{\partial }\overline{v}}{\mathrm{\partial }y}\right)}^2+{\left(\frac{\mathrm{\partial }\overline{w}}{\mathrm{\partial }z}\right)}^2\right]\\ & +{\left(\frac{\mathrm{\partial }\overline{u}}{\mathrm{\partial }y}+\frac{\mathrm{\partial }\overline{v}}{\mathrm{\partial }x}\right)}^2\\ & +{\left(\frac{\mathrm{\partial }\overline{u}}{\mathrm{\partial }z}+\frac{\mathrm{\partial }\overline{w}}{\mathrm{\partial }x}\right)}^2+{\left(\frac{\mathrm{\partial }\overline{v}}{\mathrm{\partial }z}+\frac{\mathrm{\partial }\overline{w}}{\mathrm{\partial }y}\right)}^2\}\end{array}$(8) (9) $\begin{array}{rl}{S}_{TD}& =\frac{\mu }{T}\{2\left[\overline{{\left(\frac{\mathrm{\partial }{u}^{\mathrm{\prime }}}{\mathrm{\partial }x}\right)}^2}+\overline{{\left(\frac{\mathrm{\partial }{v}^{\mathrm{\prime }}}{\mathrm{\partial }y}\right)}^2}+\overline{{\left(\frac{\mathrm{\partial }{w}^{\mathrm{\prime }}}{\mathrm{\partial }z}\right)}^2}\right]\\ & +\overline{{\left(\frac{\mathrm{\partial }{u}^{\mathrm{\prime }}}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{v}^{\mathrm{\prime }}}{\mathrm{\partial }x}\right)}^2}\\ & +\overline{{\left(\frac{\mathrm{\partial }{u}^{\mathrm{\prime }}}{\mathrm{\partial }z}+\frac{\mathrm{\partial }{w}^{\mathrm{\prime }}}{\mathrm{\partial }x}\right)}^2}+\overline{{\left(\frac{\mathrm{\partial }{v}^{\mathrm{\prime }}}{\mathrm{\partial }z}+\frac{\mathrm{\partial }{w}^{\mathrm{\prime }}}{\mathrm{\partial }y}\right)}^2}\}\end{array}$(9) where $\overline{u}$, $\overline{v}$ and $\overline{w}$ are the time-averaged velocity components in the x, y and z directions, respectively; $u^{\mathrm{\prime }}$, $v^{\mathrm{\prime }}$ and $w^{\mathrm{\prime }}$ are the fluctuating velocity components in the x, y and z directions, respectively; and μ is the dynamic viscosity. Because the velocity fluctuations cannot be obtained by solving the Reynolds time average equation, Knock proposed calculating the turbulent entropy production rate in the k-omega SST model using Equation 10 (Kock & Herwig, 2004, 2005): (10) $S_{TD}=0.09\frac{\rho k\omega }{T}$(10) where k is the turbulent kinetic energy and ω is the specific dissipation rate.

Finally, the total entropy production ΔS in the volute can be obtained by integrating the local entropy production rate with the volume of the computational domain (Zhou et al., 2022): (11) $\mathrm{\Delta }S={\int }_{V}SdV={\int }_V{S}_{VD}dV+{\int }_V{S}_{TD}dV$(11)

4. Optimization methodology

4.1. Parametric design of bionic volute

The existing volute region is divided as shown in Figure 9(a). Regions A and C are the transition regions connected with the filling valve and rotor, respectively. The shapes of these two regions cannot be changed due to design restrictions. Thus, region B is the primary region of this bionic design. The shape of region B can be represented by the typical cross-sectional profile shown in Figure 9(b), from which the existing volute profile has a bilateral symmetrical structure.

Figure 9. (a) Volute region division, (b) Typical cross-section oil distribution.

According to the aforementioned introduction to the working principle of the hydrodynamic retarder, the volute should have a good oil diversion effect to ensure the rapid braking of the hydrodynamic retarder. However, in the actual use process, due to the dynamic effect caused by clockwise rotation of the rotor, most of the oil enters the volute and moves clockwise on the right side, while only a small part of the oil flows in from the left side. Two oil streams with opposite moving directions impact the upper left part of the volute to form a vortex, which hinders the movement of the oil. Therefore, it is necessary to further optimize the design of the volute profile so that the oil can quickly and smoothly enter the working chamber of the hydrodynamic retarder.

A Nautilus is a cephalopod primarily distributed in coral reef waters of the tropical Indian Ocean-western Pacific Ocean, which is called ‘living fossil in the ocean.' The inside of the Nautilus shell is divided into many chambers by diaphragms, and each chamber is connected by a siphuncle. When a Nautilus needs to sink or float, it regulates its buoyancy by sucking seawater into or out of the chamber. When a Nautilus needs to move quickly, the seawater in the chamber is quickly ejected in the opposite direction through the nozzle at the head to predate or avoid danger. The special spiral structure and motion mode of Nautilus have made it known as the ‘ocean propeller' (Liang et al., 2021; Tinello et al., 2016). Due to the frequent intake and discharge of sea water into and out of the Nautilus shell, after hundreds of millions of years of evolution, its shell cross section has a shape that is similar to a logarithmic spiral, as shown in Figure 10 (MacQuitty, 2008). Such a structure can effectively reduce the resistance of sea water flowing in the Nautilus shell.

Figure 10. Nautilus cross-sectional shape.

Because a Nautilus shell has the advantage of low flow resistance, the bionic Nautilus volute of the hydrodynamic retarder with a cross-section profile, as shown in Figure 11, is designed via feature extraction and appropriate modification of its cross section, combined with the oil flow situation in the existing inlet volute. The resulting parametric expression is used to explore the potential of bionic design.

Figure 11. Cross-section profile of bionic nautilus hydrodynamic retarder volute.

The cross-sectional profile of the bionic volute has a right-handed logarithmic spiral shape, which uses the influence of the dynamic effect of the rotor on the oil flow in the volute when the hydrodynamic retarder is operating. The entry section is designed as a circular arc, which can quickly and smoothly change the flow direction of the oil, so that it enters the subsequent volute flow channel according to the expected flow state and concurrently can reduce the local energy loss at the entrance. The diversion section mimics the Nautilus shell and uses its horizontal eccentric design. Considering that the right eccentric structure of the Nautilus may be due to its own volume being much larger than that of the chamber, both left eccentric and right eccentric designs will be discussed in the subsequent studies. The return section is designed as a circular arc that is tangent to the diversion section. Such a design can ensure that part of the oil that does not flow into the hydrodynamic retarder will not impact the inlet oil and cause energy loss. In addition, this recirculation oil can reconverge with the inlet oil at a small angle. The rest of the profile remains unchanged. The design of the above three sections of the volute profile can be expressed by the parameters in Table 4. The range of these four parameters in the design process can be determined based on the structural constraints.

Table 4. Design variables and ranges of the volute.

Name	Symbol	Range (mm)
Entry radius	r	[125, 250]
The horizontal distance of the return section	a1	[20, 50]
The vertical distance of return section	a2	[0, 30]
Eccentricity	e	[−10, 10]

Although the volume of each Nautilus chamber increases at a ratio of approximately 1: 1.618, the thickness of the volute remains the same size as the original design due to the constraint of the installation space.

4.2. Optimization objective and constraint

Considering the problems with using the volute as designed, three optimization objectives are set in this optimization: the pressure difference between the inlet and outlet of the volute Δp, the nonuniformity coefficient of volute outlet flow γ, and the maximum velocity of the oil in the volute v_max. The formulas of the optimization objective are as follows: (12) $\begin{array}{rl}\Delta p& =p_{in}-p_{out}\end{array}$(12) (13) $\begin{array}{rl}\gamma & =\frac{\sum _i|q_i-\overline{q}|A_i}{2|\overline{q}|\sum _i{A}_i}\end{array}$(13) where p_in and p_out are the pressures at the inlet and outlet of the volute, respectively. For the outlet surface of the volute, A_i is the area of each grid element, q_i is the flow rate of each grid element, and $\overline{q}$ is the average flow rate.

The inlet volute should effectively guide the oil, so that it can flow into the working chamber of the hydrodynamic retarder quickly and smoothly, and avoid energy loss caused by oil accumulation and impact inside the volute. The pressure difference between the inlet and outlet of the volute Δp is an important index to characterize the internal energy consumption of the volute, which should be minimized in the optimization process. When the energy loss inside the volute decreases, the oil flow rate in the volute increases theoretically, thus, the maximum velocity of the oil in the volute v_max should also increase after optimization. In addition, the flow direction of oil at the outlet of the volute suddenly changes from the radial direction to the axial direction along the hydrodynamic retarder. This change will lead to an uneven distribution of oil at the outlet of the volute, which makes the torque of the hydrodynamic retarder rotor blades unequal, thus affecting the stability of the hydrodynamic retarder in operation. Therefore, the nonuniformity coefficient of volute outlet flow γ should be reduced as much as possible. The mathematical expression is as follows: (14) $\begin{array}{rl}& \text{min}\left(\Delta p,\gamma \right)\mathrm{\&}\text{max}\left(v_{max}\right)\\ & \text{s.t.}{x}_j^{\text{LB}}\le x_j\le x_j^{\text{UB}}\end{array}$(14) where x_j is the design variable, x_j^LB and x_j^UB are the upper and lower limits of the design variable, respectively.

4.3. Optimization procedure

Figure 12 shows the flow chart of optimization, which can be expressed as four steps. First, the optimal Latin hypercube design method is used to sample uniformly in the design space. The sampling points are divided into two parts: one part is used to construct the surrogate model, and the other part is used to evaluate the fitting accuracy of the surrogate model. Then, the RBF neural network surrogate model is constructed according to the sampling results calculated by CFD. If the fitting accuracy meets the requirements, then the next stage is started; otherwise, more sampling points are added to further train the surrogate model, and the above process is repeated several times until the accuracy requirements are met. Furthermore, a correlation analysis of the parameters is performed. Correlation analysis can describe the influences of different design variables and interaction items on optimization objectives, and can provide a certain degree of reference for optimization design. Finally, the NSGA-II algorithm is used to optimize the parameters of the surrogate model, and the optimization results are compared and verified. If the performance of the optimization result is worse than the original volute, the optimization process will be performed again; otherwise, the optimization result will be output to complete the optimization.

Figure 12. Bionic volute optimization design process.

4.3.1. Design of experiment

To improve the fitting accuracy of the surrogate model in the overall design space and avoid local overfitting or underfitting, it is necessary to make the sampling points fill the design space evenly. Concurrently, the number of sampling points should be reduced as much as possible to reduce the amount of calculation. Compared with the traditional Monte Carlo method, the optimal Latin hypercube design method can obtain the same statistical accuracy with fewer sample points, thus, this method is chosen for sample point selection.

One hundred sample points are selected in the design space, of which 70 are used as the training set and 30 as the test set. The corresponding three-dimensional model is constructed for each sample point, and CFD simulation is performed to obtain the response value, which is the three optimization objectives. The volute cross-sectional profiles corresponding to the sample points in the training set are shown in Figure 13. The red and blue profiles are the two boundaries of the design space, and all sample points should be included in these two profiles. The figure shows that the training set samples show a wide diversity, and thus, the selection of sampling points can meet the requirements of uniform distribution in the design space.

Figure 13. Training set volute cross-section profile.

4.3.2. Surrogate model

DOE only represents the responses of a limited number of sample points. To find the optimal design in the entire design space, it is necessary to build a surrogate model. The surrogate model builds a mathematical model based on the existing sample points, which can approximate the relationship between input variables and output variables. Compared with the complex and time-consuming CFD calculation in the optimization design process, the surrogate model has the advantage of using fewer resources to explore a larger design space.

Currently, the commonly used surrogate models are the RBF neural network model, Kriging model and Gaussian process (Dong et al., 2018). The Kriging model essentially uses the interpolation method, but its disadvantage is that it is sensitive to noisy data points and has poor generalization ability. The Gaussian process is developed based on the Kriging model, which also has some shortcomings. Considering that the RBF neural network model can effectively fit the nonlinear functional relationship between input and output, the RBF neural network is an ideal choice for this study.

The RBF neural network is a kind of three-layer feedforward neural network, that is primarily composed of an input layer, hidden layer and output layer. The input layer is the design variable vector of each sample point, and the hidden layer is obtained by nonlinear transformation of the input layer. There are several neurons in the hidden layer, and the transformation function of neurons is the radial basis function. Generally, the Gaussian function is chosen as the radial basis function, and its expression is as follows: (15) ${\phi }_j\left(\mathit{x}\right)=exp\left(-\frac{||\mathit{x}-{\mathit{c}}_j||}{2{{\delta }_j}^2}\right)j=1,2,\dots ,m$(15) where x is the design variable vector of the input sample point; ${\mathit{c}}_j$ is the centre of the radial basis function; $||.||$ is the Euclidean distance between the input vector and the centre vector; ${\delta }_j$ is the width of the radial basis function; and m is the number of neurons in the hidden layer.

The output layer is a linear combination of the radial basis function and weight coefficient, and its expression is: (16) $y_i=\sum _{j=1}^m{\omega }_{ij}{\phi }_j\left(\mathit{x}\right)i=1,2,\dots ,n$(16) where $y_i$ is each response value; ${\omega }_{ij}$ is the weight coefficient of each radial basis function; and n is the number of response values. The structure of the RBF neural network constructed by this research is shown in Figure 14.

Figure 14. RBF neural network structure.

4.3.3. Optimization algorithm

The mathematical problem in this study can be considered to be a multi-objective optimization problem in a high-dimensional design space. As one of the classical multi-objective optimization algorithms, the NSGA-II genetic algorithm is widely used in the optimization process of rotating machinery impellers and volutes due to its high running efficiency and good distribution of the Pareto solution set. The detailed theory of this algorithm is reported in that paper (Deb et al., 2002). In this study, the NSGA-II genetic algorithm is used to optimize the parameters of the abovementioned multi-objective optimization problem. Detailed settings are shown in Table 5.

Table 5. Settings for NSGA-Ⅱ.

Parameters	Value
Population size	100
Maximum number of generations	200
Crossover fraction	0.9
Mutation fraction	0.1

5. Results and discussions

5.1. Correlation analysis

Generally, for multi-objective optimization problems, the objective functions often conflict with each other, and it is impossible to optimize all objectives concurrently. Therefore, we can only choose between the optimization objectives and obtain a relatively optimal solution set that accounts for most of the optimization requirements, this set is referred to as the Pareto optimal solution set. To select solutions that meet the optimization requirements in the Pareto optimal solution set, parameter correlation analysis is required to evaluate the influence degree of design variables on optimization objectives and provide guidance for the selection of optimum points.

In this study, the Pearson correlation coefficient is used to measure the correlation between the design variables and optimization objectives. The expression is as follows: (17) $p_{xy}=\frac{\sum _{i=1}^n\left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)}{\sqrt{\sum _{i=1}^n{\left(X_i-\overline{X}\right)}^2}\sqrt{\sum _{i=1}^n{\left(Y_i-\overline{Y}\right)}^2}}$(17) where X = [x₁,x₂, … ,x_n] and Y = [y₁,y₂, … ,y_n] are any two sets of data; $\overline{X}$ and $\overline{Y}$ are the average values of X and Y, respectively; and n is the number of samples.

The range of the Pearson correlation coefficient is [−1, 1]. p_xy greater than zero indicates a positive correlation of variables, while a value less than zero indicates a negative correlation of variables. The closer the absolute value is to 1, the stronger the correlation is.

Figure 15 shows the results of the correlation analysis between the optimization objectives and design variables. Results show that r and a₂ play a vital role in the optimization of the three objectives. a₁ primarily affects v_max and γ, but has a negligible effect on ΔP. e has a strong influence on ΔP and γ, but a weak influence on v_max. According to the optimization direction of the three objectives and the sign of the Pearson correlation coefficient of the design variables, the optimization trend of each design variable is described by Table 6.

Figure 15. Correlation between optimization objectives and design variables.

Table 6. Optimization trend of design variables.

	ΔP	vmax	γ
r	Min	Min	Max
a1	Max	Max	Min
a2	Min	Min	Max
e	Max	Max	Min

Table 6 shows that a₁ and e must be as large as possible and r and a₂ as small as possible to reduce ΔP and increase v_max. However, the trend of each design variable is the opposite of the optimization of γ. a₁ and e should be reduced as much as possible, while the values of r and a₂ should be increased, which can also indicate that the optimization objectives in the multi-objective optimization problem often conflict with each other. Considering that the primary goal of this study is to reduce energy loss inside the volute and improve the torque response characteristics of the hydrodynamic retarder, ΔP and v_max should have higher weights in the optimization process, thus, a₁ and e should be as large as possible, and r and a₂ should be as small as possible when selecting the optimal parameters.

5.2. Optimum results

5.2.1. Regression analysis for surrogate model

After constructing the surrogate model, its fitting accuracy must be checked. Generally, the deterministic coefficient R², root mean square error RMSE, maximum absolute error MAE and maximum relative error MRE are used to measure the fitting accuracy of the surrogate model. The R² evaluation index is used in this study. When the R² of the response value is greater than 0.9, the surrogate model is considered to have sufficient fitting accuracy for the response value. The expression of R² is as follows: (18) $R^2=1-\frac{\sum _{i=1}^n{\left(f_i\left(x\right)-\widehat{f_i}\left(x\right)\right)}^2}{\sum _{i=1}^n{\left(f_i\left(x\right)-\overline{f}\left(x\right)\right)}^2}$(18) where $f_i\left(x\right)$ is the actual value of the response; ${\hat{f}}_i\left(x\right)$ is the approximate value calculated by the surrogate model; ${\overline{f}}_i\left(x\right)$ is the average value of $f_i\left(x\right)$; and n is the number of test points.

Analysis results are shown in Figure 16, and it can be seen that the predicted results are in good agreement with the CFD calculation results. Table 7 shows the R² values of the RBF neural network surrogate model for the three objective functions. R² are all greater than 0.9, which meets the standards of engineering applications. Therefore, it is considered that parameter optimization based on this surrogate model has sufficient reliability.

Figure 16. Regression analysis for the surrogate model.

Table 7. R² analysis results.

Objective function	R^2
Δp	0.94922
vmax	0.93836
γ	0.90963

5.2.2. Pareto frontiers

The NSGA-II method is used to optimize the parameters of the bionic volute, and 20,000 solutions are obtained. The solutions in the final solution set that conform to the optimization trend are selected as shown in Figure 17. The figure shows that a stable Pareto boundary appears after iteration, which indicates the trade-off between conflicting objective functions.

Figure 17. Computational Pareto frontier.

Obviously, the optimization of one objective function will inevitably lead to degraded performance in another objective function. The optimization results show that most Pareto solutions have a lower pressure difference between the inlet and outlet of the volute Δp and a higher maximum velocity of the oil in the volute v_max under given design conditions. However, the cost is that the oil distribution at the outlet of these volutes will be more chaotic, thus, it is necessary to choose the most suitable scheme according to the actual engineering needs. In this study, the optimization degree of the objective function must be at least 25%, and the performance sacrifice degree of some objective functions must not exceed 25%.

5.2.3. Comparison of geometry and performance

After averaging and rounding several points in the Pareto solution set that meet the above requirements, the final result of this parameter optimization is obtained, and the optimized bionic volute is verified by CFD simulation. Unfortunately, according to the requirements of our partners, we cannot show the specific parameters. The comparison of the objective functions and profile between the original prototype and the bionic volute are shown in Table 8 and Figure 18, respectively.

Figure 18. Comparison of profiles (a) prototype (b) bionic volute.

Table 8. Comparison of the objective function between the prototype and bionic volute.

	Δp (Pa)	γ	vmax (m/s)
Prototype	1766.6	1.134	4.8937
Bionic volute	371.95	1.377	6.8809

Table 8 shows that the pressure difference between the inlet and outlet of the bionic volute is reduced by 78.95% compared with the original prototype, which shows that the oil in the bionic volute has less energy loss in the flow process and can flow into the working chamber of the hydrodynamic retarder more smoothly, thus, the maximum velocity of the oil in the bionic volute is also increased by 40.61%. After optimization, the nonuniformity coefficient of the volute outlet flow is increased due to the change in the oil flow form in the volute, but it is still within the allowable range, thus, the optimization goal has been achieved.

5.3. Optimization mechanism study

5.3.1. Internal flow field analysis

To explore the mechanism of reducing flow energy loss in bionic volutes, the velocity field and streamline diagram of a typical section of the volute under steady-state condition are compared and analysed.

As shown in Figure 19, the upper left part of the volute of the original prototype forms a large vortex area due to the impact of inlet oil and return oil. The deposition of oil in this study reduces the effective flow area of the volute, which is not conducive to the smooth flow of oil into the volute. In addition, there are also extensive flow separation phenomena in the lower left part, and these separation vortices will also affect its working performance. Compared with the prototype volute, the streamline distribution in the bionic volute is markedly improved, all the oil moves uniformly along the volute profile, and there is no marked flow separation in the volute. These internal flow changes lead to better hydraulic performance of the bionic volute. Correspondingly, the flow rate of oil inside the volute increases significantly, thus, the bionic volute has higher liquid filling efficiency.

Figure 19. Typical section velocity field and streamline diagram comparison.

Furthermore, the Omega method is used to identify vortices in the original prototype and the bionic volute. As the latest generation vortex identification method, the Omega method is widely used in the field of hydraulic machinery research because it is not affected by subjective threshold selection and can display both strong and weak vortices concurrently. The mathematical expression of the Omega method is as follows: (19) $\mathrm{\Omega }=\frac{||B|{|}_F^2}{||A|{|}_F^2+||B|{|}_F^2+ϵ}$(19) (20) $ϵ=1/1000\ast max\left(||B|{|}_F^2-||A|{|}_F^2\right)$(20) where matrices A and B are the symmetric tensor and antisymmetric tensor of the velocity gradient respectively, and $||\cdot |{|}_F^2$ is the square of the matrix norm.

Figure 20 compares the vortex distribution inside the two volutes, and the corresponding velocity cloud map is also attached to the surface of the vortex. Figure 20(a) shows that there are many low-speed vortices in the corresponding vortex area of the original prototype volute in Figure 19(a), which nearly occupies the entire flow area of the volute. However, in the bionic volute shown in Figure 20(b), except for the separation vortex at the entrance and the wall vortex in part of the return section, there are basically no vortices in the flow channel of the entire volute. Therefore, the bionic volute can effectively inhibit vortex generation by improving the internal oil flow.

Figure 20. Flow field vortex identification comparison.

When comparing the performance of the volute before and after optimization, the nonuniformity coefficient of the flow field at the outlet of the bionic volute is improved. Figure 21 compares the flow distribution at the outlet of the two volutes. The red region represents the smooth flow of oil into the working chamber, and the blue region represents the return of oil to the volute. The flow at the outlet of the two volutes is shown to not be evenly distributed along the circumferential direction but rather exhibits a marked regional distribution. The oil inside the volute of the original prototype primarily flows into the working chamber through the upper left of the outlet, which is called the oil-filled region. However, there is a certain reflux phenomenon in the lower right, which is called the reflux region. Because the internal oil flow state of the bionic volute is changed, the flow distribution at the outlet is also different from that of the original prototype. The oil-filled region is concentrated directly above the outlet, and the reflux region is located at the lower left of the outlet. In addition, the flow intensity at the outlet of the two volutes is also different, and the intensity of the oil-filled region and the reflux region of the bionic volute is much greater.

Figure 21. Comparison of outlet flow distribution.

According to the coordinate axis shown in Figure 21, the circumferential velocity and pressure field at the outlet of the volute are analysed with the upper part as the starting point and clockwise as the positive direction. Results are shown in Figure 22. The axial velocity is used to replace the distribution of outlet flow to show the direction of oil flow more clearly. The axial velocity is positive due to the oil flowing into the volute, and the axial velocity is negative dur to the oil flowing into the working chamber.

Figure 22. Comparison of pressure and velocity at the outlet of the volute.

The analysis of the velocity field shows that there is a high-velocity region between θ = 90° and θ = 270° in the original prototype volute, while the bionic volute has a higher velocity region between θ = 130° and θ = 320° due to less internal energy loss. According to the relationship between velocity and pressure in the flow field, the pressure shows a marked down-wards trend at the corresponding high-speed region of the two volutes, which leads to the uneven circumferential distribution of pressure at the outlet of the volute.

Thus, the pressure at the outlet presents two trends. For the prototype volute, θ = 0° to θ = 80° and θ = 190° to θ = 360° are in the forwards pressure gradient region, and θ = 80° to θ = 190° are in the inverse pressure gradient region. For the bionic volute, θ = 0° to θ = 120° and θ = 260° to θ = 360° are in the forwards pressure gradient region, and θ = 120° to θ = 260° are in the inverse pressure gradient region. The oil in the volute is affected by the rotor dynamic effect, and its flow state is primarily clockwise, that is, from θ = 0° to θ = 360°. According to the analysis of the pressure field and axial velocity, the oil can flow smoothly into the working chamber when the direction of movement is along the pressure gradient. When the direction of oil movement is the inverse pressure gradient, the pressure gradient will make oil flow back to the volute. Because the oil flow rate inside the bionic volute is higher, the pressure in the low-pressure region is lower, resulting in stronger oil reflux. This is also the reason for the large nonuniformity coefficient of its outlet flow.

5.3.2. Energy characteristic analysis

To explore the optimization mechanism in more detail, entropy production theory is used to compare the energy loss distribution in the flow field of two volutes. Figure 23 shows that the entropy inside the volute is primarily produced by two regions: the near-wall region, where the viscous force acts in the boundary layer; and the high Reynolds number region caused by unstable flows such as vortex and reflux.

Figure 23. Comparison of flow energy loss distribution in the volute.

There are three regions with high entropy production in the prototype volute, two of which correspond to the vortex in Figure 19(a). Although vortex 1 is large in volume, its entropy production is limited due to its low velocity. The high-speed and high-intensity vortex 2 produced by flow separation at the wall is the primary source of entropy production in the flow field. In the lower right corner of the volute there is also a certain degree of entropy production due to reflux, but the entropy production is low. Therefore, the unreasonable structure of the prototype volute leads to the phenomenon of oil impact and flow separation inside it, resulting in a large energy loss.

Due to the reduction in the vortex structure, the primary source of entropy production in the bionic volute is near the wall. Although the reflux phenomenon is more serious than that in the prototype volute, the entropy production in this study will not have a strong adverse effect on the overall energy loss of the volute. In summary, the bionic volute can restrain the generation of vortices by improving the internal oil flow, thus reducing the energy loss.

5.4. Transient braking process analysis

The optimization of the volute is shown in its own improved performance and in the overall braking performance of the hydrodynamic retarder. Figure 24 shows the comparison of the emergency braking torque characteristics of two hydrodynamic retarders with different volutes under design conditions. The braking characteristic indices of both are shown in Table 9.

Figure 24. Emergency braking characteristics of the hydrodynamic retarder.

Table 9. Braking characteristic index of the hydrodynamic retarder.

	Onset time	Max braking torque
Prototype	0.76 s	5960.15 N·m
Bionic volute	0.71 s	6326.49 N·m

Figure 24 shows that the bionic volute hydrodynamic retarder has a larger slope in the torque rising stage, which means that it can achieve a larger braking torque concurrently, effectively improving the braking performance of the hydrodynamic retarder and driving safety. The data in Table 9 also show that the peak braking torque of the bionic volute hydrodynamic retarder is increased by 6.15%, and the onset time is advanced by 6.58%.

These changes in braking torque characteristics are explained in more detail in Figure 25. Generally, the oil passes through the filling valve and the inlet volute, one part enters the working chamber of the hydrodynamic retarder, and the other part stays in the inlet volute until it is full of the volute. This stage is called the initial filling stage, which corresponds to 0-0.3 s in Figure 25. After the volute is basically filled with oil, its internal flow state tends to be stable, and the subsequent oil can smoothly flow into the working chamber of the hydrodynamic retarder. This stage is called the torque effective stage, which corresponds to 0.3-0.9 s in Figure 25. The figure shows that the rising rate and value of the filling rate of the bionic volute are larger than those of the prototype volute, which means that the volute can establish an internal stable flow field more quickly, ensuring the rapid rising of the filling rate of the working chamber of the hydrodynamic retarder, and thus improve its braking performance.

Figure 25. Comparison of filling rate changes.

Figure 26 shows the gas–liquid two-phase distribution diagram in two volutes during liquid filling. After the initial stage of liquid filling, although the rate of the prototype volute has been high, a stable flow field has not been established inside it. Due to the vortex produced by oil impact, there is still a large cavity region in the upper left part, and this cavity region continues until the end of the torque effective stage. Although there are some cavities in the bionic volute at 0.3 s, the overall oil flow tends to be stable, and the same flow state is basically maintained in the subsequent torque effective stage.

Figure 26. Gas-liquid two-phase distribution during liquid filling.

To sum up, the bionic Nautilus volute effectively inhibits the generation of vortex and improves the oil flow speed by improving the oil flow state inside it, which can make the oil flow into the working chamber of hydrodynamic retarder quickly and smoothly, thus producing better braking effect and improving the driving safety of heavy vehicles. At the same time, this bionic design method can also guide the design and optimization process of volute of centrifugal compressor, centrifugal pump and other equipment.

6. Conclusion

In this paper, a design method for the bionic Nautilus volute profile of hydrodynamic retarder is proposed. Four bionic parameters are selected as design variables, 100 sample cases are generated by the OLHS method, and CFD simulation calculations are carried out. Then, the correlation between the design variables and the optimization objective is analysed, and the NSGA-II algorithm based on the RBF neural network surrogate model is used to optimize it. Finally, the internal energy loss characteristics of the two volutes are analysed in detail, and the reasons to improve the performance of the bionic volute are explored. In addition, the braking performances of the hydrodynamic retarder equipped with the two volutes are compared. The following conclusions are reached:

1. Regression analysis show that the R² of the three optimization objectives are all greater than 0.9, which shows that the relationship between the design variables fitted by the RBF neural network surrogate model and the optimization objectives has high accuracy, and it is feasible to use the optimization objective model to optimize the parameters.

2. The pressure difference between the inlet and outlet of the bionic volute decreased by 78.95%, and the internal oil flow velocity increased by 40.61%. Although the nonuniformity coefficient of the outlet flow increased, it is still within the allowable range, and thus, the optimization goal has been achieved.

3. By analysing the flow and energy characteristics of the two types of volutes, the bionic volute is shown to effectively mitigate vortex generation by improving the internal oil flow, thus reducing the energy loss of the internal flow and improving the performance.

4. The emergency braking peak torque of the hydrodynamic retarder with the bionic volute is increased by 6.15%, the onset time is advanced by 6.58%, and the braking torque response characteristics are improved.

The bionic nautilus volute of the hydrodynamic retarder studied in this paper has shown a good effect in engineering application. However, due to the restriction of time cost in the optimization process, only four main parameters are used to express the cross-sectional profile of the volute, which makes the shape of the bionic volute still have a certain gap with the actual nautilus volute. In the follow-up research process, more design parameters will be used to express the shape of the volute, and the dimension curse problem in the optimization process will be solved based on the parameter model embedding method, so as to further explore the potential of bionic nautilus volute in the rapid effect of braking torque of hydrodynamic retarder.

Disclosure statement

No potential conflict of interest was reported by the author(s).