A prediction model of wall shear stress for ultra-high-pressure water-jet nozzle based on hybrid BP neural network

ABSTRACT Two hybrid back-propagation neural network (BPNN) models optimized by two heuristic search algorithms, namely genetic algorithm (GA-BP) and particle swarm optimization (PSO-BP), are proposed in this paper to predict radial maximum wall shear stress instead of traditional computational fluid dynamics (CFD) methods. The two proposed models are trained and validated using a database of 150 radial maximum wall-shear-stress values obtained via CFD simulations. The best fit model is identified from various BPNN models based on the Akaike information criterion (AIC) and the Bayesian information criterion (BIC), balancing the trade-off between goodness-of-fit and model complexity. The model performance is evaluated by MAE, RMSE, regression coefficient (R), and Nash-Sutcliffe efficiency (NSE). The best fit model is a three-layered BPNN model consisting of a 4:4:1 topology. In almost all evaluation indicators, the two hybrid BPNN methods outperform three existing algorithms, namely classical BPNN, random forest (RF), and gradient boosting decision tree (GBDT). Both PSO-BP and GA-BP can provide a more precise assessment of radial maximum wall shear stress, with maximum errors being 5.81% and 8.24% respectively. The proposed PSO-BP prediction model is promising and has great feasibility in predicting the radial maximum wall shear stress of UHP water-jet nozzle in engineering applications.


Introduction
With the rapid globalization of economy, the exchange of international trades, and the increasing frequency of human activities, the number of ships has been expanding at an unprecedented rate in recent years (Gu et al., 2020). However, the submerged surfaces of ship hull are highly susceptible to bio-fouling and rusting after longterm travel and berthing. It reveals that fouling organisms corrode the ship hull and increase the roughness of hull surface, thereby reducing the shipping speed (Gu et al., 2015;Schultz, 2007) and increasing fuel consumption and greenhouse gas emissions (Dupraz et al., 2018;Feng et al., 2018). The governments worldwide are increasingly concerning about this environmental and energy consumption issue. It becomes meaningful to clean up the hull attachments and corrosion spots when they grow up to be problems.
As a novel and environmentally friendly surfacecleaning technique, the UHP water jet is widely used to remove hull rust and marine attachments. Since the water jet cleaning technology allows for quick washing jobs in narrow, harsh situations without causing second-hand environmental pollution, it has gained wide acceptance, CONTACT Zheng-Shou Chen aaaczs@163.com including governments and shipbuilding companies. The UHP water-jet nozzle is a crucial unit in implementing the water-jet technology in practice. In this context, enhancing the operational efficiency of water-jet technology by improving the hydrodynamic performance of the UHP water-jet nozzle becomes a topic of intense interest.
In the past few years, most of the studies have concentrated on improving the rust removal efficiency of wallclimbing robots, such as selecting different adsorption modes, movement controls, and drive schemes . However, only a few studies were conducted on hydrodynamic performance analysis of water-jet nozzle. Undoubtedly, the hydrodynamic performance of the water-jet nozzle directly reflects the impinging exfoliation ability of the UHP water-jet device. The impinging exfoliation via a UHP water jet imposed on the target surface generally comes into play to remove the corrosion and other residual coatings on the target surface only when the impinging force (composed of wall shear stress and water wedge effect) of the water jet is more significant than a particular threshold value. Several lines of evidence suggest that these parameters, i.e. standoff distance (L), jet angle (θ ), length to diameter ratio (μ), and contraction angle (α), influence the radial maximum wall shear stress from water impingement (Kim et al., 2022). Thus, assessing the effect of these parameters on radial maximum wall shear stress forms one of the primary focuses of this work. The accurate estimation of radial maximum wall shear stress of impinging water jet remains an essential step in assessing the hydrodynamic performance of nozzle. Wall shear stress generated by an impinging supersonic water jet is notoriously hard to be measured (Kim et al., 2022). Theoretical and experimental investigations on wall shear stress attributed to impinging jet can be found in the literature. Young et al. (2013) measured the wall shear stress created by an impinging supersonic water jet using oil-film interferometry. Their experiments were promising, but the oil-film interferometer's precision is restricted. Tu and Wood (1996) used the Preston tube measurement method to explore wall shear stresses created by subsonic impinging jet. However, their results are constrained by the measuring devices, making extrapolation into compressible jets challenging. Smedley et al. (1999) and Phares et al. (2000) investigated the removal of microspheres via impinging jets and calculated the wall shear stress of target plate using theoretical relationships. Those studies established a link between wall shear stress and particle forces, but they did not account for compressibility or turbulent effects. Loureiro and Freire (2012) showed that laser Doppler anemometry could be used to measure velocity within 50 μm of a wall. On the other hand, the viscous layer thickness crucial to microparticle clearance can be much less than 20 μm (Kottapalli & Novosselov, 2021). According to these previous studies, it is difficult to measure the wall shear stress directly.
Plenty of numerical researches were done to address this issue regarding the supersonic impinging water jet. The Reynolds averaged Navier-Stokes (RANS) models were used in most CFD investigations in this field (Xiao et al., 2020;Xing et al., 2022), as they have lower computational costs and produce good predictions when compared to experimental data. Chin et al. (2013) used turbulence models, such as k-ε and k-ω-SST, to numerically study the compressible impinging at various jet angles and standoff distances. In parallel, high-speed optical technologies such as Schlieren photography and shadowgraph were used to verify the simulation results. Although CFD approaches provide a potential to develop a comprehensive and accurate prediction of radial maximum wall shear stress resulting from UHP impinging water jet in a complex working condition, it is impracticable to evaluate hydrodynamic performance using such a computationally expensive method.
As a result, creating a semi-phenomenological prediction approach is required to attain high accuracy and low processing cost in a numerical solution based on physics. Artificial intelligence has recently attracted considerable attention, and various machine learning (ML) approaches, e.g. artificial neural networks (ANNs), support vector machine (SVM), RF, GBDT and so on, have been extensively implemented to model data in numerous applications. Singh et al. (2022) proposed two hybrid ML-based pedotransfer functions, i.e. two combinations of GA with multilayer perceptron and SVM respectively, to predict hydraulic conductivity. Deng et al. (2022) developed a framework for modeling, analyzing, and forecasting the spatiotemporal variations of water quality indicators based on ML algorithms and environmental Kuznets curves (EKC) analysis. Malik et al. (2022) estimated the pan evaporation prediction based on two innovative techniques: deep learning and gradient boosting machine models. In the water-jet fields, some researchers also made many attempts at ML methods, and their findings suggest that ANNs can be appropriate techniques for the prediction of physical characteristics of water jets (Ficko et al., 2021;Jou, 2000;Onen, 2014;Zhang et al., 2019). The use of ANNs for maximum wall-shear-stress prediction of impinging water jet, especially the prediction based on CFD modeling schemes, remains largely unexplored.
This research aims to filling this gap by implementing hybrid BPNN models to predict the radial maximum wall shear stress of impinging water jet. In Section 2, the efficiency of CFD simulations is validated through the comparison with existing experimental data. It makes certain to use CFD method with confidence in the following text. Next, the input variables (i.e. L, θ, μ and α) and their related value ranges are carefully selected. In order to create the simulation database, 150 simulation jobs are created by Latin hypercube sampling (LHS) method. In Section 3, using simulation data, a thorough analysis of BPNN configuration is carried out to identify an optimal network structure. Furthermore, two powerful heuristic search optimization methods, namely GA and PSO, are utilized to choose best the initial weights and thresholds for the best performing BPNN models. In Section 4, a rigorous statistical and graphical evaluation process demonstrates that the two developed hybrid BPNN models effectively predict the radial maximum wall shear stress generated by UHP impinging water jet. Section 5 compares the proposed hybrid BPNN algorithm, PSO-BP with GA-BP and other existing algorithms, i.e. classical BPNN, RF, and GBDT, to validate its prediction performance. Finally, a summary of key findings and conclusions are presented in Section 6.

Geometric modeling
As shown in Figure 1, the straight cone convergent nozzle is extensively used in firefighting, hydraulic coal mining, and other industrial fields as a high-efficiency nozzle with good jet performance (Li et al., 2010;Wen et al., 2016). Because almost all of the straight cone convergent nozzles have a distinct character, i.e. concentrated velocity distribution, its hydrodynamic characteristics are commonly superior to other types of nozzles. Moreover, this kind of nozzle can be used under ultra-high injection pressure up to 300 MPa, which is suitable for ship rust removal operations.
The principal role of a straight cone convergent nozzle is to translate static pressure supplied by a high-pressure pump device into the dynamic force of UHP impinging water jet. Therefore, a continuous cone water jet can be obtained. As shown in Figure 2, the longitudinal-section (along axis) structure of straight cone convergent nozzle can be separated into three primary sections. The access section (I), the contraction section (II), and the exit section (III). First, the access section (I) connects a high-pressure hose. Second, the contraction section (II) is used to shrink the flow longitudinal section of the nozzle's internal chamber to raise the fluid velocity rapidly. Third, the exit section (III) is used to stabilize the flow state of fluid and enable it to be ejected from the outlet at high velocity, eventually forming a high-pressure water jet with concentrated energy (Sun, 1992). In addition, the structural design of straight cone convergent nozzle must  obey with the Equation (1) below: where D, α, d, l 2 , l 3 are the inlet diameter, contraction angle, outlet diameter, contraction section length, and exit section length, respectively. Liu et al. (2021) simulated the performance of a 3D impinging jet nozzle using CFD methods. The results demonstrate that the axisymmetric 3D model can accurately describe the flow field features of actual situation. The cylindrical symmetric nozzle, found by Chin et al. (2013), may utilize the symmetry of flow problem, requiring just half of the fluid domain to be modeled, and their simulation results match well with experimental results. Therefore, half of the fluid field can enhance calculation efficiency instead of using the entire fluid field for the simulation of a cylindrical symmetric nozzle. Half of the 3D model is chosen as the study objective in this paper, since the straight cone convergent nozzle is also characterized by a cylindrical axisymmetric feature. The computation domain comprises fluid volume inside the nozzle's internal chamber and fluid volume between the water jet outlet and target wall surface. Figure 3 shows a 3D depiction of the computational domain, where R a is the radius of the outer cylindrical computational domain. The static pressure inlet boundary condition is imposed on the cross-section (perpendicular to central axis of water jet) of the access section (I). Water jet from the nozzle impinges the target wall, spreads, reflects and ultimately splashes, leaving the domain via the cylindrical pressure outlet boundary condition. According to the numerical simulation setup performed by Jaramillo et al. (2012), R a should be large enough to illustrate the impinging details after the water jet impacts on the target wall. Under this criterion, the R a was set to be 60d (i.e. R a = 15 mm), which can sufficiently satisfy the requirement of capturing more details related to water jet and splashing droplets.

Multiphase model
The CFD method is the most common way to calculate the wall shear stress, particularly in the cases related to UHP impinging water jet, instead of experiments. The UHP impinging water jet is a turbulent flow in the liquid-vapor-gas three-phase flow field (Xiao et al., 2020). The multiphase volume of fluid (VOF) model available in STAR-CCM+ is chosen to effectively simulate the flows inside the nozzle chamber. The VOF model can predict the interface shape of two or more immiscible fluids (Xiao et al., 2020). The concept of volume fraction related to each phase is introduced in this model and the interface shape can be obtained by calculating the volume fraction of the phases in each control volume. All multiphase components at any given position have the same velocity and pressure. In the computational domain, the equations regulating continuity and momentum are as follows (Cebeci & Bradshaw, 1977): The phase transfer equations of the three phases are listed as follows: whereṘ e andṘ c are the phase transition mass transfer rates between the vapor and liquid phases, respectively. The volume fractions of the liquid, vapor, and air phases are represented by α l , α v , and α i respectively. The Rayleigh-Plesset equations are utilized to establish the transport relationship between the liquid and vapor phases in the current mainstream cavitation model (Zhang et al., 2021). A homogeneous equilibrium model based on barotropic fluids is an option for addressing the phase shift of liquid to the liquid-vapor mixture (Nezamirad et al., 2022). However, due to discrepancies in the density and pressure gradients during the expansion and compression of vapor bubbles, this model makes it challenging to capture the impact of baroclinic torque. In the present work, the transport equation-based model was adopted.
The Schnerr-Sauer cavitation model, which is adequate for complicated orifice flow circumstances, is used to model the source termsṘ e andṘ c given above (Guo et al., 2017;Yu et al., 2017). Condensation and vaporization rates are defined as follows: where the saturation vapor pressure of liquid is defined as P v , and ρ is the density of the mixture defined as (8), which expresses the relationship between the number of bubbles per unit volume N b of pure liquid and the bubble radius R B .
Here, the only required input parameter N b for the numerical solution was assigned as default value of 10 12 , and it was testified that this is an optimal value (Liu et al., 2012). The density of bubble nuclei in the fluid strongly relates to the bubble number density. According to conclusion drawn by Li et al. (2014), it was found the vapor volume appears to be unaffected by the higher bubble density than the default amount. On the other hand, a smaller bubble density is likely to indicate an enormous vapor volume. The overshoot in source intensity fluctuations, proportional to bubble density, might be one plausible explanation for this occurrence. This overshoot appears to be a pseudo-proposition cavitation model with little to do with physics. As a result, in this investigation, the default value of N b = 10 12 is employed.

Turbulence model
In order to choose an appropriate turbulence model for the numerical simulations, the result comparison associated with three different turbulence models, i.e. SST k-ω, RNG k-ε, and Realizable k-ε models, is shown in Figure 4, where the radial distribution of wall shear stress, 2τ w /(ρv 0 ), is illustrated on the y-axis. Although the differences among 2τ w /(ρv 0 ) attributed to different turbulence models are observable, their radial distributions have the same trend from experimental observation when 2 < r/b < 9, where r is the distance from 'O' on the target wall shown in Figure 3, and b is the half-width of the impingement pressure profile (Loureiro & Freire, 2012). By comprehensive comparison with experiments (Loureiro & Freire, 2017;Poreh et al., 1967), the predicted values of 2τ w /(ρv 0 ) through simulation cases via RNG kε model agree well with the experimental results, rather than cases related to the other two turbulence models. In addition, Li et al. (2017) and Celik et al. (2014) verified that the RNG k-ε model can get favorable results in the CFD simulations concerning cavitation flow. Here, the RNG k-ε turbulence model was adopted to evaluate the hydrodynamic performance of UHP water jet nozzle. Yakhot and Orszag (1986) developed RNG k-ε turbulence model, which is capable of simulating flows with high strain rates and huge curvature with acceptable accuracy. The following transport equations are used to calculate the turbulence kinetic energy k and the rate of dissipation ε: where G k is the turbulence kinetic energy generated as mean velocity gradients. The C 1ε and C 2ε coefficients are set to be 1.42 and 1.68 in default respectively. For k and ε, the values α k and α ε are the inverse effective Prandtl numbers. The effective viscosity μ eff in the case of a high Reynolds number is stated as:

Simulation setup and boundary conditions
It is critical to select appropriate boundary conditions for simulation convergence and accuracy. The boundary conditions are marked in Figure 3, together with the dimension details of the computational domain. The crucial boundaries that are employed are stated as follows. The static UHP up to 200 MPa and standard atmosphere pressures are imposed on the inlet and outlet boundaries respectively, as listed in Table 1. The no-slip boundary conditions are used on the inner nozzle wall and the target wall, with the average gradient of scalar function set to be zero. Proper simulation settings can accelerate simulation convergence. A SIMPLEC scheme is used to implement the pressure-velocity coupling scheme. Density, turbulent transports, and Reynolds stresses are interpolated using the first-order upwind approach. In order to discretize the momentum equation, the QUICK scheme was chosen. The compressive scheme is used for discretizing the vapor volume fraction transport equation. For the pressure interpolation, the body force weighted scheme was chosen. The least-squarecell-based approach is used to determine the cell gradients. Finally, a bounded first-order implicit scheme is used to conduct time integration. Table 2 summarizes the fluid thermodynamic properties involved. Using the Tait equation of state, the liquid is supposed to be compressible, but the vapor is assumed to be incompressible (Koukouvinis et al., 2017). The Tait equation, which represents the weak compressibility of liquid water, is shown where B, n, ρ l and P sat are the bulk modulus, material exponent, liquid density and saturation pressure respectively, and their values are provided in Table 2.

Grid independency analysis
The grid density has a significant impact on the accuracy of numerical simulation. Undoubtedly, a refined mesh can improve the simulation accuracy of the impinging process for the UHP water jet. For high-speed multiphase flow simulations, the grid density should simultaneously cater to the accuracy requirement and speed calculation. Therefore, there are strict requirements for mesh quality inside the volumes where the speed of multiphase media is high, or the phases change rapidly, e.g. contraction field. On the other hand, a coarse mesh can be adopted surrounding the refined mesh volumes. The graded hexahedral grid topology is suitable for this kind of computational domain with significant differences in speed and phase, because it can provide more accurate solution to the flow field interested and spend less time to solve the fluid field unconcerned. According to the knowledge of fluid mechanics, a more satisfactory mesh resolution is used in the domain center neighborhood. As the distance from the centerline increases, the mesh resolution becomes coarser, as shown in Figure 5. Grid-independency analysis was carried out firstly to quantify the effect of the proposed CFD model on the numerical simulation results. Based on the same mesh topology scheme, four simulation cases with increasing grid density, namely Mesh-1 to Mesh-4, were performed. The mesh configurations of the proposed CFD model and the obtained maximum pressure are shown in Table 3. As  seen in the table, the maximum pressure value achieved using Mesh-3 is similar to that obtained via Mesh-4, indicating that increasing the number of mesh cells does not affect the calculation results, once a certain number of mesh cells is reached. The total number of elements is estimated to be around 1.3 million, with an average element quality of 0.97. Hence, the mesh topology and grid density utilized in Mesh-3 were employed in the following numerical simulations.

Database establishment
The preliminary experimental design for selecting initial sampling points becomes essential before establishing a surrogate model for some experimental data. In order to obtain an excellent initial sample with the fewest number of sampling points and a good space-filling effect, a sampling method based on LHS was used in this  aper. By using parameter sets that are non-overlapping in any of the dimensions, LHS may effectively sample over multi-dimensional parameter space (Bulut et al., 2021). Thus, the LHS is useful for dispersing the training/testing points more uniformly over the design space. Model parameters used for LHS are given in Table 4, along with a plausible range for each parameter. The LHS algorithm generates the design space for training and testing points (including L, θ, μ, and α values). The MAT-LAB algorithm used to generate the LHS of the design space is shown in Appendix A. Once the training points from the design space are generated, the next step is to perform CFD simulations on the training points (i.e. L, θ, μ, and α values) to calculate radial maximum wall shear stress.

Basic principles of BPNN
BPNN is a typical data processing system that stimulates brain activity and can simulate and forecast data information rapidly and correctly (Belkacem et al., 2017). The classical BPNN consists of the input, hidden, and output layers, as shown in Figure 6. In this study, the input layer comprises 4 input variables, i.e. L, θ, μ, and α. The hidden layer comprises several processing units called neurons that use activation or transfer functions to convert input signals into output signals when they process their inputs. The unit step function, linear function, sigmoid function, and hyperbolic tangent function are commonly used activation functions. The output layer outputs the prediction result (i.e. radial maximum wall shear stress) for the given inputs. The basic BPNN algorithm uses signal forwardpropagation and error back-propagation to map the correct output via related inputs . When forward propagation is active, the information received from the input layer is sent to the output layer after transiting through the hidden layer. Error back-propagation occurs when the actual and intended output values do not match. The output error correction modifies each layer's weight and threshold value in a specific way and propagates back to the input and hidden layers. The input, hidden, and output layers of BPNN autonomic learning process continually alter the weight and threshold values, which is also the process of forwarding information propagation and error back-propagation. The learning process is complete when the output error is smaller than the default, or the number of epochs approaches the predefined value (Xu et al., 2019). The steps in the classical BPNN training process are shown as follows: Step 1: BPNN initialization. First, determine the number of neurons n, l, m of the input layer, hidden layer, and output layer of BPNN based on the input and output characteristics of the system. Next, randomly initialize four parameters, including two connection weights (ω ij between input and hidden layer neurons, and ω jk between hidden and output layers), the threshold a for the hidden layer, and the threshold b for the output layer. Finally, Give the activation functions and the learning rate.
Step 2: Hidden layer output calculation. The output H of the hidden layer is calculated according to the input vector X, the connection weight ω ij , and the threshold a, expressed by Equation (13).
In Equation (14), f (x) is the excitation function of the hidden layer. Commonly used activation functions include the tanh, sigmoid, relu, and softmax. The sigmoid function was adopted as the activation function in this paper, and its mathematical expression is given as follows: Step 3: Output layer input calculation. According to the hidden layer output H, connection weight ω jk , and threshold b, the prediction output O p of the BP neural network is calculated as follows: Step 4: Error calculation. According to the network forecast output O p and the expected output Y, the network forecast error e can be calculated as follows: Step 5: Weight update. Update the network connection weight ω according to the network prediction error e.
Step 7: Determine whether the iteration of algorithm is ended. If not, return to Step 2.

Genetic algorithm
GA is a stochastic non-linear optimization approach based on artificial intelligence that simulates natural selection and genetic mechanisms (Katoch et al., 2021). It is a random and efficient search and optimization approach that aims to identify optimal or suboptimal solutions to problems that are difficult to solve algebraically (Athanasiou & Konkoli, 2020;Sun et al., 2020). In this paper, the GA algorithm is used to optimize the BPNN by modifying the connection initial weights and thresholds to minimize MSE, and its flowchart is depicted in the subplot (A) of Figure 7. The parameters associated with GA algorithm are initialized and the detailed parameter settings used in this paper are summarized in Table 5. All weights and thresholds of one network are encoded into a chromosome X i (i = 1, 2, · · · , N p ), which is contained in a population of randomly initialized individuals {X i } Np i=1 . The fitness value for each chromosome is assessed using the fitness function . Subsequently, the iterative process comes to a halt if the preset stopping condition is met and the individual fitness values are  recorded at the moment. For the optimization algorithm outcome, two scenarios are evaluated, if the results are satisfactory. Otherwise, it is necessary to conduct an additional evaluation of genetic operators, which include three critical operations: selection, crossover, and mutation. The selection operation is performed by the roulette wheel selection approach based on the fitness value for selecting the elite populations for crossover. It starts with the calculation of the fitness value [X i (gen)], the relative fitness value r [X i (gen)] and, the cumulative fitness value c [X i (gen)] for the X i (gen), where gen is the generation counter. The relevant formulae are given as follows: Generate a random number r ∈ [0, 1] and compare it with c [X i (gen)] for i = 1, 2, . . . , N p . If c [X z−1 (gen)] < r < c [X z (gen)], X z gets chosen to be a part of the new population. A new population of N p members is selected by repeating this selection process. Two chromosomes from a population are selected for crossover operation. A random selection test is performed to determine whether these two chromosomes will undergo a crossover, and a random number p (between 0 and 1) is generated for the random selection test. The chromosome will be selected if p < P c . Similarly, another chromosome will be selected. Two crossover points are randomly selected with equal probability on the chromosomes of each individual, and the part of the parent chromosomes between selected points will be swapped with crossover probability P c . The genes after the crossover point will be exchanged to create two new chromosomes. The operations are repeated until all chromosomes have been considered.
The mutation operation also starts with a random selection test for each chromosome. If a generated random number p ∈ [0, 1] for a chromosome is larger than P m , the chromosome undergoes mutation. A random number will be generated for the chosen component with a value within the component limits. The procedures will be repeated until all chromosomes have been considered. After recalculating the fitness value of each chromosome, the best member X B (gen) that has the largest fitness value and the worst member X w (gen) that has the smallest fitness value can be identified. X B (gen) will be compared with the best one in the last generation (i.e. X B (gen-1)). If the fitness value of X B (gen) is smaller than the one of X B (gen-1), the content of X B (gen-1) will replace the content of X B (gen). Afterward, the content of X B (gen-1) will be substituted into X w (gen-1). The individual with the best fitness in the final generation is used to obtain optimal weights and thresholds for the network.

Particle swarm optimization
PSO is a heuristic search algorithm initially presented by Eberhart and Kennedy (1995) and inspired by the behavior of swarms and the dynamic movement of animal species such as fishes. It attracts much attention for tackling engineering challenges involving intricate optimization and searching (Poli, 2008;Wang et al., 2018). Similar with GA, the PSO algorithm is used to find the best optimal weights and thresholds for BPNN instead of the gradient descent method and it can enhance generalization performance. The special PSO algorithm which combines the BPNN algorithm is named PSO-BP. The PSO algorithm starts with a random population generation that includes particle velocities and locations. Then, each particle is evaluated using the fitness function, followed by the stopping conditions. The optimal weights and thresholds are found if the particle fitness meets the requirement. Otherwise, the velocity and position of every particle must be updated. The new velocity and location of particles are calculated using Equations (23) and (24), respectively.
The velocity and position at ith and (i + 1)th iterations are represented by V i , V i+1 , X i and X i+1 , respectively. 'rand()' is a number between 0 and 1 that is chosen at random. C 1 and C 2 are the learning factors of PSO algorithm. ω is the weighting factor for speeding up convergence, with its formula expressed by Equation (25), (25) where the terms ω min and ω max are the smallest and largest weighting factors respectively, referring to the iteration number. The reevaluation of the particles is the final phase (Le et al., 2019). The structure of the proposed PSO to determine the optimal BPNN parameters is depicted in the subplot (B) of Figure 7. The setting parameters of PSO used in this paper are summarized in Table 6.

Hyper-parameter tuning of BPNN
In this subsection, the validation of an experiment to choose crucial hyper-parameters in BPNN model (such as the number of hidden layers and the training functions) was carried out. Training the network by altering hyper-parameter values and selecting the one that delivers the best result on the validation set, are required for each hyper-parameter. This method can be used to experimentally tune the hyper-parameters of BPNN model, in order to improve its performance. It is worth noting that only the training and validation sets are used in this procedure, while the test set is always kept unseen.

Number of layers and neurons
Aforementioned, previous studies have demonstrated that L, θ , μ, and α are four key influencing factors of radial maximum wall shear stress concerning straight cone convergent nozzle. Thus, the input layer has four neurons in total. Aiming at predicting the radial maximum wall shear stress in a nonlinear system, the number of neurons in output layer is set to be one. BPNN can feature multiple hidden layers between input and output layers (Campbell et al., 2017) . Although too many hidden layers have more vital generalization ability, it makes the algorithm convergence difficult. Numerous studies proved that a three-layer neural network might mimic practically any nonlinear mapping function, if there are enough hidden neurons. Thus, the prediction model adopted a three-layered neural network in this paper.
An important step in modeling is to determine the number of hidden layer neurons. During the selection of hidden layer neurons, although numerous neurons may minimize training errors, they take longer time to train and may cause overfitting issues. Hence, the number of hidden layer neurons must be optimized while the model complexity and accuracy are considered. The most commonly used traditional model selection criteria are AIC (Akaike, 1998) and BIC (Schwarz, 1978). They are based on information theory and are driven by the need to balance the goodness-of-fit and model complexity (Marasinghe, 2011). A lower AIC and BIC values indicate that the model fits the data better. Only a brief introduction is given here because extensive theoretical results regarding the derivation and performance of AIC and BIC have already been published (Burnham & Anderson, 2004). This study is conducted based on Qi and Zhang (2001) proposed form. The mathematical expression of AIC is given in Equation (26) below: where N is the number of observations, and p is the number of parameters in the model.σ 2 MLE is the maximum likelihood estimation, expressed by Equation (27): where SSE is the sum of square error, y i andŷ i represent the calculated and predicted results of radial maximum wall shear stress, respectively. The typical expression of BIC is similar to that of AIC and is defined by Equation (28). Because the BIC method imposes a more substantial penalty on model complexity, inherent in the method is a preference for simpler models (Burnham & Anderson, 2004).
For complex non-linear models, p is typically defined as a natural measure of model complexity (Burnham & Anderson, 2004). The number l (i.e. the degrees of freedom) of hidden layer neurons is a critical uncertain parameter in the three-layer neural network structure. Therefore, the parameter l is employed to measure the model complexity of a three-layered neural network. In this study, l takes values from 2 to 11 (Shen et al., 2008;Zhu et al., 2017). The AIC and BIC values of the testing dataset corresponding to the different numbers of hidden layer neurons are calculated according to Equations (26) and (28). As shown in Figure 8, it can be observed that the performance measures firstly decrease and subsequently increase with the number of hidden layer neurons increasing. According to the AIC and BIC metrics, the proposed BPNN model presents the best performance (i.e. the smallest AIC and BIC) when l equals to 4. Thus, the best fit model is a three-layered BPNN model consisting of a 4:4:1 topology.

Training function
As mentioned in subsection 3.1, the learning and training rules of the BPNN are to use error back-propagation to continuously update the weights and biases of individual neurons, in order to minimize the network RMSE. A variety of training functions are provided by the BPNN, including trainlm, traingd, traingda, traingdx, trainrp, traincgb (Liang et al., 2020). In order to determine the optimal training function, the RMSE values of the test datasets are compared under various training functions.
As it is shown in Figure 9, in all the training functions, the RMSE obtained by 'trainlm' function is better than the others. Therefore, it was chosen as the training function. The most popular Levenberg-Marquardt (L-M) algorithm was used in this prediction model, combining with the best features of the Gauss-Newton and gradient descent algorithms (Chen et al., 2022). The L-M algorithm can correct the defects such as the absence of inverse matrix Gauss-Newton algorithm. Its significant merits are characterized by quick convergence rate, difficult convergence to local minimum values, and fewer requirements for training sample features. Thus, the L-M algorithm is ideal for this moderate and small data scale, ensuring accuracy without compromising training speed.

Performance evaluation
The prediction performance of the hybrid BPNN models are assessed using different statistical measures. The MSE, RMSE, MAE, NSE, Willmott's index of agreement (WIA), and confidence index (CI) are some of these indicators. The following are the formulas for those indicators (Keshtegar & Seghier, 2018;Seghier et al., 2020): where N s denotes the number of data sets, and y i , y i and y i denote the calculated, predicted, and average values of the radial maximum wall shear stress, respectively. It is worth noting that the prediction model with the lowest MSE, RMSE, and MAE values are thought to have the best accuracy. However, the prediction model with the greatest CI values has a more robust and efficient performance.

Results and discussion
In the existing problem, initially the dataset consisting of 150 samples having four variables (namely, L, θ, u, and a) and the corresponding radial maximum wall shear stress are considered, as shown in Table 7. The optimal sample is No.52 nozzle (where L = 17.38 mm, θ = 12.74°, μ = 1.52 mm, α = 54.23°) with the corresponding radial maximum wall shear stress 1.3899 MPa. The 3D spatial distributions of 150 samples are shown in Figure 10. Among the 150 samples, 100 randomly chosen samples are used to train purpose and the remaining 50 individual solutions are used to test the hybrid BPNN models. The hybrid BPNN models are implemented in MATLAB using the 'newff ' function. In the subsubsection 3.4.1, it is shown that the hidden layer counts 4 neurons, which is the size giving the best prediction performance. The detailed parameter settings related to hybrid BPNN structure are shown in Table 8. The L-M algorithm is used for training purpose, and the MSE criteria are used to test performance of predicted data. · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 51 45.82 9.02 5. · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·   The primary purpose of this paper is to evaluate the feasibility of two hybrid BPNN models, GA-BP and PSO-BP, and then to select the better one for predicting maximum radial wall shear stress, in order to overcome the limits and disadvantages attributed to CFD methods. Because the BPNN is trained through the training set and its hyper-parameters are tuned via the validation set, it implies that the BPNN models have 'seen' these data in advance. Results based upon these sets may be overly optimistic. An unseen dataset is required to evaluate their performance, providing greater confidence in using them. Therefore, the randomly produced test sets are used to assess the performance of these two hybrid BPNN models. The performance and detailed error analysis of several critical factors in model prediction are presented in this section.
A comparison between the calculated and the predicted values of radial maximum wall shear stress related to GA-BP and PSO-BP models are shown in subplots (A) and (B) of Figure 11. The predicted results of these two hybrid BPNN models are in good agreement with the calculated values. The two hybrid BPNN models differ in the testing phase, and it shows a more precise fitting of the data via the PSO-BP model than that via the GA-BP model.  As shown in Figure 12, subplots (A) and (B) are used to examine the performance of the GA-BP and PSO-BP models, respectively. The training, validation, testing, and overall phases are presented individually, with the regression coefficient R value in each subplot. The R measures the correlation between the inputs (predicted output) and the targets (calculated output), with a value nearer to one indicating a strong correlation. The solid line shows the linear regression, whereas the dashed line represents a near-perfect forecast. Thus, it is clear that the best fit between the computed and predicted outputs are provided by the BPNN model combined with the PSO algorithm. R results in the PSO-BP model are highest in the training, validation, testing, and overall phases, at 0.99031, 0.98889, 0.97565, and 0.98626, respectively, with the testing phase being lower than the validation phase because the missing information leads to  an imperfect match. Training results are satisfactory for GA-BP and PSO-BP models, and regression coefficients are greater than 0.9 for training, validation, testing, and overall phases.
In Figure 13, the subplots (A) and (B) show the relative error between predicted and calculated values for the GA-BP and PSO-BP, respectively. It becomes clear that the relative error of the radial maximum wall-shear-stress values of the two-hybrid BPNN models did not exceed ±0.1. The absolute values of the maximum relative error of GA-BP model in the training and testing phases are 8.24% and 4.63%, respectively; the corresponding values of PSO-BP model are 5.81% and 4.25%, respectively, much smaller than those of GA-BP model. Thus, the prediction accuracy of PSO-BP model is higher than GA-BP model.
The stability of the proposed hybrid BPNN models and the results of CFD calculation correlation during prediction process of radial maximum wall shear stress for the input variables L, θ, μ, and α are investigated.  The ratios of the calculated to predicted values of radial maximum wall shear stress for GA-BP and PSO-BP models are plotted in the subplots (A) and (B) of Figure 14, respectively. In the two subplots, y andŷ represent the calculated and predicted results of radial maximum wall shear stress, respectively. The upper stabilization limit is illustrated by a red dashed line. The stable values of parameters L, θ, μ, and α indicate that the smaller the oscillation between y andŷ is, the closer the value of y/ŷ gets to one. The PSO-BP model provides better prediction results with no notable trends regarding the variables, rather than the GA-BP model.
As an important efficiency-evaluation index, the laboratory-time costs and machine overheads were also taken into account. The running of a given simulation case typically takes about 24 h on a desktop computer with an 8-core 3.7 GHz processor. As a result, it will take approximately 150 days to predict the radial maximum wall shear stress for 150 cycles, making hydrodynamic performance assessment via CFD method quite timeconsuming. For the abovementioned task, the proposed hybrid BPNN models can predict the radial maximum wall shear stress within only a few seconds. Although the hybrid BPNN methods have higher prediction accuracy and faster calculation speed, they cannot interpret predictions of radial maximum wall shear stress, also known as the black-box nature. In fact, the unexplainable blackbox models created by BPNN or other ML techniques are widely recognized. Furthermore, the details on the datasets, modeling process, parameter tuning, and evaluation results are not adequately explained and provided in the vast published literature. This lack of transparency (also known as black-box behavior) dramatically hampers the interpretability of the prediction results and even becomes a primary problem for trusting such methods. One of the main aspects on which future efforts should be focused is the reproducibility of the results.

Validation
In the first place, aiming at validating the prediction ability of the proposed PSO-BP algorithm, the results from it are compared with those associated with the GA-BP, the classical BPNN and another two powerful algorithms, namely RF and GBDT. In this regard, a brief information about the RF and GBDT algorithms is presented.

Random forest algorithm
The RF algorithm is a machine learning approach based on classification and regression trees, and it is another powerful algorithm for constructing a predictive model (Breiman, 2001). The RF can handle a lot of noise and is less prone to over-fitting issues. An ensemble of decision tree models is fitted to a data set using RF regression models. In order to increase the predictability of the response variable, the data is recursively split for each tree into more homogeneous units, known as nodes. Split points are determined by the values of the predictor variables (Everingham et al., 2016). As a result, the variables utilized to split the data are regarded as crucial explanatory variables. A predetermined number of bootstrapped datasets are used in the RF algorithm to fit various decision trees. The mean fitted response from all the individual trees produced from each bootstrapped sample is the predicted value of the continuous response. In this paper, a random forestbased regressor is built via the 'randomForest' package (https://code.google.com/archive/p/randomforest-matla b/) based on the algorithm of Breiman (1996).

Gradient boosting decision tree algorithm
The GBDT combines boosting and decision tree methods and is one of the finest ML techniques for fitting actual distributions. The decision tree is a fundamental classification and regression tree method (CART), and it has the benefits of quick classification and model visualization (Zhai et al., 2020). In order to improve performance, boosting learns multiple classifiers by linearly combining them after varying the weights of the training samples. The main idea is that the gradient direction of the model loss function is established each time when the model is built, causing the loss function to decrease in the gradient direction. The GBDT algorithm trains multiple learners for complex tasks using a gradient descent approach. It then combines the outcomes of many learners to do classification, better than the way it could have done with just a single learner. This Python package 'sklearn' was adopted in this paper to implement the GBDT regression model with default values (Pedregosa et al., 2011).

Validation results
In order to verify the performance of the proposed hybrid BPNN methods, the statistics performance comparison among the PSO-BP, GA-BP, BPNN, RF and GBDT is made. Tables 10 and 11 show the details of essential hyper-parameters for the RF and GBDT algorithms. Moreover, to have a better view of the models, RMSE, R, NSE, and MAE values of all models are given in Table 12. Based on the results presented in Table 12, it can be found that the proposed hybrid BPNN method, PSO-BP,  is superior to the other four algorithms in almost all evaluation indicators. All the validation results adequately demonstrate the accuracy and reliability of PSO-BP.

Conclusions
The contribution of this study consists of proposing two new frameworks (namely GA-BP and PSO-BP) based on the hybrid BPNN model for accurate prediction of the radial maximum wall shear stress, which is required for the hydrodynamic performance estimation of UHP water-jet nozzle. It is verified that the PSO-BP is more promising than GA-BP. Consistent with previous literature findings, it is found that the standoff distance L, jet angle θ , length to diameter ratio μ, and contraction angle α are essential factors influencing the radial maximum wall shear stress. For both hybrid BPNN models, the LHS method is applied to obtain the design space (including L, θ , μ, and α values), and a CFD method is utilized to generate training and testing data. The results based on the AIC and BIC demonstrate that the best fit model is a three-layered BPNN model consisting of a 4:4:1 topology. In almost all evaluation indicators, the PSO-BP method is superior to the other four methods, namely GA-BP, classical BPNN, RF, and GBDT. It can provide a more precise assessment of radial maximum wall shear stress, with maximum errors being 5.81%. In addition to verifying the traditional predictive accuracy, stability analysis and regression analysis are important ways to validate the effectiveness of the novel model. The prediction PSO-BP model can be continuously updated to provide more precise estimation of the radial maximum wall shear stress if the database is continuously supplemented with data from a broader range.
In summary, the attempt to combine the CFD method and the hybrid BPNN algorithm can significantly reduce computation time during the assessment of hydrodynamic performance to ensure prediction accuracy. Like other prediction models, this study has some limitations, because we only investigated one specific type of nozzle with a relatively simple geometric structure. For future study, we will focus on expanding this framework to more nozzle types, which can dramatically inspire the design of higher-efficiency nozzles.

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