High-specific-speed rotating blade multi-source constraint design method based on improved grey wolf algorithm and Poisson disc sampling

Abstract

The invention discloses a high-specific-speed rotating blade multi-source constraint design method based on an improved grey wolf algorithm and Poisson disc sampling, and belongs to the technical field of high-specific-speed rotating blade design. The method comprises the following steps: firstly, constructing a multi-source constraint region, giving an initial model and characteristic simulation data thereof, and setting meridians and other geometric characteristics of the front edge and rear edge contours of a blade; calculating the inlet angle of the initial model blade and the vortex velocity difference of the front and rear edges of the blade, screening design parameters, and calculating and determining the design variable range through the dimensionless vortex velocity. Then, uniformly sampling in the limited area by using a Poisson disk sampling algorithm, ensuring that sampling points are not mutually overlapped, and dynamically updating a sampling set; and finally, applying a multi-objective optimization equation, constructing an agent model, obtaining a Pareto solution set through an improved grey wolf optimization algorithm, and selecting an optimal parameter combination from the Pareto solution set to improve the design precision and performance of the rotating blade. The method has the characteristics of no deviation, minimum distance and maximum sampling area coverage.

Description

Technical Field

[0001] The invention relates to the field of high specific speed rotary blade design, in particular to a multi-objective optimization design method for high specific speed rotary blades. Background Art

[0002] Pump-jet propulsion is an advanced propulsion method that uses jet water to propel underwater vehicles forward. Due to its efficient propulsion capability, anti-cavitation, low vibration and low underwater noise, it is widely used in fields such as petroleum engineering, transportation and shipping, and amphibious vehicle design. As the core component of the pump-jet propulsion turbine, the blades directly affect its performance. With the increasing requirements for the cruising speed of underwater vehicles, the demand for new blades with high specific speed is increasing. At present, many researchers optimize the design of rotating blades by adjusting geometric variables, using various intelligent algorithms and inverse design methods. However, the existing design methods still need to be improved, especially in the control of the number and range of design variables. At high specific speeds (exceeding the critical specific speed n sc ), existing methods have limited effect in performance optimization, especially when facing cavitation problems in high-speed working environments. Summary of the invention

[0003] In order to solve the problems existing in the prior art, the present invention provides a multi-source constraint design method for high specific speed rotating blades based on an improved Grey Wolf algorithm and Poisson disk sampling. The method combines computational fluid dynamics (CFD) analysis, multi-source constraint screening and multi-objective algorithm to increase the number of control sections and reduce design variables to improve the optimization efficiency. The propulsion efficiency can be effectively improved during the optimization process, which is of great significance for promoting the development of high-speed pump-jet propulsion power equipment.

[0004] The technical solution adopted by the present invention is: a multi-source constraint design method for high specific speed rotating blades based on improved grey wolf algorithm and Poisson disk sampling, comprising the following steps:

[0005] Step 1. Construct a multi-source constraint region to preliminarily screen the range of blade design parameters, including the following sub-steps:

[0006] Step 1.1: Given high specific speed rotating blade design parameters; the high specific speed rotating blade design parameters include given high specific speed rotating blade design speed n, design flow rate Q d 、Design head H d , inlet diameter D, calculate the rotating blade specific speed ns, determine the blade flow path parameters; set the meridian of the leading and trailing edge contours of the blade, the number of blades, the blade thickness, and the superposition law, select the wrap angle θ as the only factor to adjust the shape of the rotating blade, and control the rotating blade shape by changing the wrap angle distribution along the streamline.

[0007] Step 1.2: Determine the initial model and calculate the inlet vortex velocity r1V in the initial model simulation data θ1 and the vortex velocity difference ΔrV between the leading and trailing edges θ (s), the formula is defined as follows:

[0008]

[0009]

[0010] Among them, V m1 is the inlet meridian velocity; r1 is the inlet meridian streamline radius; ω is the blade design speed; β b1 is the blade inlet angle; u1 is the circumferential velocity of the inlet meridional streamline; δ is the correction value of the angle; g is the gravitational acceleration; H d is the head under design conditions; η h is the hydraulic efficiency, which is approximately equal to 1; γ is the correction factor, which is approximately equal to 0.

[0011] Step 1.3: Introduce and define a quartic polynomial to represent the dimensionless eddy velocity ξ(s), construct a multi-source constraint region, and convert the calculation formula of the dimensionless eddy velocity into the value of two unknown variables (Φ, Ψ) in the multi-source region under the constraint. The calculation formula is as follows:

[0012]

[0013] Step 1.4: Based on the initial model, calculate the design variables (Φ h ,Ψ h ),(Φ m ,Ψ m ),(Φ s ,Ψ s ) and specify the above-mentioned limited sampling area Ω. The subscripts h, m, and s are the values ​​at the hub, middle section, and wheel cover, respectively.

[0014] Step 1.5: Calculate the vortex velocity rV based on the above parameters θ (s), the calculation formula is as follows,

[0015] rV θ (s) = r1V θ1 +ξ(s)ΔrV θ

[0016] Among them, r1V θ1 and ΔrV θ Calculated from the original model, ξ(s) is the dimensionless eddy velocity.

[0017] Step 1.6: Calculate the wrap angle distribution along the streamline. In order to construct the geometry of the rotating blade, three sections including the hub, the middle section and the wheel cover are used. Each section is controlled by multiple scattered points. Through joint simplification, the wrap angle of a point on each section is directly calculated, and the results are as follows:

[0018]

[0019] Among them, (r1V θ1 ) h ,(r1V θ1 ) m ,(r1V θ1 ) s are the inlet vortex velocities at the hub, middle and cover respectively; (ΔrV θ ) h ,(ΔrV θ ) m ,(ΔrV θ ) s are the difference between the inlet and outlet vortex velocities; (Φ h ,Ψ h ),(Φ m ,Ψ m ),(Φ s ,Ψ s ) are the variables of constraint area control; θ h,j ,θ m,j and θ s,j Represent the wrap angle of the j-th control point of the hub, middle and wheel cover respectively.

[0020] Step 1.7: Take the inlet vortex velocity (r1V θ1 ) h ,(r1V θ1 ) m ,(r1V θ1 ) s , the difference in vortex speed between inlet and outlet (ΔrV θ ) h ,(ΔrV θ ) m ,(ΔrV θ ) s , dimensionless eddy current velocity control variable (Φ h ,Ψ h ),(Φ m ,Ψ m ),(Φ s ,Ψ s ) is the optimization variable; based on the initial model simulation calculation results, the inlet and outlet vortex speeds and the inlet and outlet vortex speed difference variable ranges are determined; based on the multi-source constraint area, the dimensionless vortex speed control variable range is determined.

[0021] Step 2. Given the number of sampling points N and the sampling area Ω, use the Poisson disk sampling algorithm combined with distribution density to accurately control the number of sampling points while maximizing the Poisson disk sampling to uniformly sample the design variables in the limited area, including the following sub-steps.

[0022] Step 2.1: Randomly sample N points in the region Ω to obtain the initial sampling point set

[0023] Step 2.2: Perform Delaumay triangulation on the gaps of the sampling point set X. Specifically, based on the Voronoi diagram theory, locate the gap triangle as a triangle t that satisfies Π(t)>0. The formula is as follows:

[0024] Π(t)=dis(p,c t )-r>0

[0025] Where p is a vertex of triangle t, r is the sampling radius, c t is the Voronoi center of triangle t, dis(p,c t ) for p to c t The Euclidean distance d is as follows:

[0026]

[0027] This step is mainly to analyze whether there are gap triangles in the current point set and their specific locations. If there is no gap, it is considered that the sampling point set has been maximized.

[0028] Step 2.3: Find the shortest edge of all triangles in the sampling point set X and delete an endpoint v on the edge. Regarding the removal strategy of the shortest edge endpoint, the maximum and minimum neighborhood average edge length removal strategy is used, where the neighborhood average edge length of the endpoint is defined as:

[0029]

[0030] Where N(v) is the number of neighboring edges of endpoint v (edges with one endpoint v), l i is the length of the neighborhood edge. Step 2.4: Repeat step 2.2 to analyze whether there is a gap triangle in the point set to detect the gap area. If so, extract the primitives, determine the specific location, calculate the gap area, and continue to step 2.5; if not, it is considered that the sampling area is completely covered by the Poisson disk, terminate the loop, and jump to step 2.6.

[0031] Step 2.5: Perform primitive extraction operations on the gap area to extract the gap area not covered by the sampling disk; generate a sampling point according to the Poisson disk sampling criterion and insert it into the point set X, and dynamically update the Delaumay triangulation. Jump to step 2.2 for iteration.

[0032] Step 2.6: It is assumed that the sampling radius of the point set X can completely cover the sampling area, so the new sampling point is inserted into the center of the largest circumscribed circle, the iteration is terminated, and the final sampling result is generated.

[0033] Step 3. Determine the multi-objective optimization equation; construct the WAS agent model; carry out multi-objective optimization work with the help of the improved multi-objective constrained grey wolf optimization algorithm; obtain the Pareto solution set and select the required optimal parameter combination from it, including the following steps.

[0034] Step 3.1: Determine the input variables, output functions, sample libraries and constraints involved in the optimization, and determine the basic settings of the algorithm such as population size and maximum number of iterations;

[0035] Step 3.2: Initialize algorithm parameters such as population size, maximum number of iterations, and population position; set the algorithm to have three different situations: search strategy flag_1, encirclement strategy flag_2, and attack strategy flag_3, and initialize flag_1=1, flag_2=0, flag_3=0.

[0036] Step 3.3: Based on the sample library, construct a WAS proxy model to generate fitted values ​​for the population individuals. Specifically, for all individual proxy models f i (x) The predicted values ​​are weighted averaged as shown below:

[0037]

[0038] Among them, w i is the proxy model f i The weight of (x) and satisfies ∑ i wi=1 and w i ≥ 0. Weight w i The determination of is based on the proxy performance evaluation of the root mean square error, and the corresponding weight allocation strategy is shown as follows:

[0039]

[0040] Among them, ε i The root mean square error of the data points in the test data set. At the same time, the goodness of fit R of the surrogate model is calculated 2 , to evaluate the accuracy of the surrogate model, the formula is as follows:

[0041]

[0042] Among them, Y i is the true value, is the predicted value of the surrogate model, is the average of the true values. 2 ≥0.9, it is considered that the accuracy of the constructed proxy model meets the requirements.

[0043] If the offensive strategy is entered, that is, flag_3 = 1, the cross-validation error is calculated, and the formula is as follows:

[0044]

[0045] Among them, i is the i-th sample point, n i is the number of samples in the validation set, y j is the true value of the sample point, The proxy model fits the sample point. The sample library is screened to build a WAS proxy model with higher accuracy requirements. The screening rules are as follows:

[0046] D t ={x1,...,x d μ-3σ<x i <μ+3σ,i=1,...,d}

[0047] Among them, μ is the mean value of the current sample library, and σ is the standard deviation of the current sample library.

[0048] Step 3.4: Perform non-dominated sorting on all individuals based on the fitted values. Specifically, considering that the parameter space is a limited region with constraints, the non-dominated partial order sorting under constraints is updated as follows: (1) Any feasible solution has a better non-dominated rank than any infeasible solution; (2) All feasible solutions calculate the clustering distance according to the objective function value, and the larger the clustering distance, the better the rank; (3) For infeasible solutions, the sorting with smaller constraint function violation value is prioritized. Ensure that all individuals are layered according to non-dominated relationships, and the individuals contained in each layer are not dominated by any other individuals in the current layer.

[0049] Step 3.5: Calculate the crowding distance CD of all individuals in the population to maintain the diversity of the population and prevent the algorithm from converging to the local optimal solution too early. The formula is as follows:

[0050]

[0051] Among them, CD im represents the crowding degree of the i-th individual in the m-th objective function, f m represents the mth objective function, x i-1 and x i+1are the previous individual and the next individual of individual i on target m, respectively, x min and x max They represent the individuals with the minimum and maximum values ​​under the m function respectively.

[0052] Step 3.6: Use the tournament selection method to update the new generation of α, β, and δ leaders in the Pareto frontier according to the non-dominated sorting results and the crowding distance; if the encirclement strategy is entered, that is, flag_2 = 1, then on the basis of the above, it is also necessary to satisfy that the α, β, and δ leaders are all in different clusters. The formula is as follows:

[0053]

[0054] Where k is the number of cluster centers, d is the Euclidean distance between two sample points, the formula is the same as the distance calculation formula in step 2.2, θ is the distance threshold for cluster classification, set to the adaptive change type, the formula is as follows:

[0055]

[0056] Where m and n are the initial value and final value of the adaptive factor, respectively. The value must be in the range [0,1] and m>n; g is the current number of iterations; g max is the maximum number of iterations.

[0057] Step 3.7: Update the positions of other gray wolves according to the positions of the leaders α, β, and δ. Given control parameters such as convergence factor, use the position update formula in the multi-objective gray wolf algorithm; consider the position of the leader and the current position of the gray wolf, and simulate the behavior of the gray wolf tracking prey to help the algorithm shift from exploration to development. Step 3.8; Perform this round of multi-objective optimization work in the complete limited area and update the current Pareto frontier; if flag_1 = 1, further narrow the optimization space and limit the area to the optimization space composed of the leaders α, β, and δ. The formula is as follows:

[0058]

[0059] Among them, Upper is the upper limit of the optimization space, and Floor is the lower limit of the optimization space. In addition, if the search strategy is entered, that is, flag_1 = 1, the optimization goal is set to maximize the agent uncertainty, and the formula is as follows:

[0060] in and are the fitting values ​​of the single agent model respectively; if the encirclement strategy or attack strategy is entered, that is, flag_2=1 or flag_3=1, the optimization goal is set to find the minimum value of the objective function, and the formula is as follows:

[0061]

[0062] Where f is the fitted value of the WAS proxy model.

[0063] Step 3.9: Repeat steps 3.3 to 3.6 to build the WAS proxy model; calculate the fitted values ​​of the Pareto frontier; determine the non-dominated profiles and α, β, δ leaders.

[0064] Step 3.10: According to the position of the α leader, determine the strategy selection corresponding to the current situation. If flag_1=1, set flag_1=0, flag_2=1; if flag_2=1, when the new α leader can dominate the original α leader, set flag_2=0, flag_3=1, otherwise keep flag_2=1; if flag_3=1, when the new α leader can dominate the original α leader, set flag_3=0, flag_1=1, otherwise keep flag_3=1.

[0065] Step 3.11: Determine whether the current iteration meets the maximum number of generations or other termination conditions. If not, jump to step 3.7 and continue to repeat the above process; if it meets the conditions, output the current Pareto front solution. Furthermore, the Poisson disk sampling criterion in step 2 specifically needs to meet the following conditions:

[0066] (1) Unbiased sampling property, that is, every uncovered point in the sampling area has the same probability of accepting a new sampling point;

[0067] (2) minimum distance property, that is, the distance between any two sampling points is greater than the given sampling radius;

[0068] (3) Maximization property, that is, the sampling area is completely covered by all sampling disks.

[0069] The beneficial effects of the present invention are:

[0070] The method first constructs a multi-source constraint region, gives the initial model and its characteristic simulation data, and sets the meridians and other geometric features of the leading and trailing edge profiles of the blade. The inlet angle of the initial model blade and the vortex velocity difference between the leading and trailing edges of the blade are calculated, the design parameters are screened, and the range of the design variables is determined by dimensionless vortex velocity calculation. Subsequently, the Poisson disk sampling algorithm is used to uniformly sample in the limited area to ensure that the sampling points do not overlap with each other, and the sampling set is dynamically updated. Finally, the multi-objective optimization equation is applied to construct a proxy model, and the Pareto solution set is obtained through the improved gray wolf optimization algorithm, from which the optimal parameter combination is selected to improve the design accuracy and performance of the rotating blade. This method has the characteristics of no deviation, minimum distance and maximum sampling area coverage, and is suitable for the efficient design of high specific speed rotating blades.

[0071] This method combines computational fluid dynamics (CFD) analysis, multi-source constraint screening and multi-objective algorithm to increase the number of control sections and reduce design variables to improve optimization efficiency. The optimization process can effectively improve propulsion efficiency, which is of great significance for promoting the development of high-speed pump-jet propulsion power equipment. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 It is a flow chart of a multi-source constraint design method for high specific speed rotating blades based on an improved grey wolf algorithm and Poisson disk sampling provided by an embodiment of the present invention.

[0073] Figure 2 This is a multi-source constraint area distribution diagram provided by an embodiment of the present invention.

[0074] Figure 3 It is a schematic diagram of the distribution of the hub, middle section, wheel cover position and pump-jet propeller blade angle control points provided by an embodiment of the present invention.

[0075] Figure 4 It is a schematic diagram of a Poisson disk sampling algorithm combined with distribution density provided in an embodiment of the present invention.

[0076] Figure 5 It is a schematic diagram of the target function distribution of the sample library provided by an embodiment of the present invention.

[0077] Figure 6 It is a mapping diagram between the optimization parameters of the parameter set and the objective function provided by the embodiment of the present invention.

[0078] Figure 7 It is a distribution diagram of goodness of fit results provided by an embodiment of the present invention.

[0079] Figure 8 It is the Pareto solution set obtained in the multi-objective constrained optimization provided by the embodiment of the present invention.

[0080] Fig. 9It is a comparison diagram of geometric differences between the initial model and the optimized model provided by an embodiment of the present invention.

[0081] Fig.10 It is a performance comparison diagram of the initial model and the optimized model provided by the embodiment of the present invention.

[0082] Specific implementation method The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0083] The present invention proposes a multi-source constraint design method for high specific speed rotating blades based on an improved grey wolf algorithm and Poisson disk sampling. The specific process is as follows: Figure 1 The multi-objective optimization design method for high specific speed rotating blades includes:

[0084] Step 1.1: Given high specific speed rotating blade design parameters; the high specific speed rotating blade design parameters include given high specific speed rotating blade design speed n, design flow rate Q d 、Design head H d , inlet diameter D, calculate the rotating blade specific speed ns, determine the blade flow path parameters; set the meridian of the leading and trailing edge contours of the blade, the number of blades, the blade thickness, and the superposition law, select the wrap angle θ as the only factor to adjust the shape of the rotating blade, and control the rotating blade shape by changing the wrap angle distribution along the streamline.

[0085] Step 1.2: Determine the initial model and calculate the inlet vortex velocity r1V in the initial model simulation data θ1 and the vortex velocity difference ΔrV between the leading and trailing edges θ (s), the formula is defined as follows:

[0086]

[0087] Among them, V m1 is the inlet meridian velocity; r1 is the inlet meridian streamline radius; ω is the blade design speed; β b1 is the blade inlet angle; u1 is the circumferential velocity of the inlet meridional streamline; δ is the correction value of the angle; g is the gravitational acceleration; H d is the head under design conditions; η h is the hydraulic efficiency, which is approximately equal to 1; γ is the correction factor, which is approximately equal to 0.

[0088] Step 1.3: Introduce and define a quartic polynomial to represent the dimensionless eddy velocity ξ(s), construct a multi-source constraint region, and convert the calculation formula of the dimensionless eddy velocity into the value of two unknown variables (Φ, Ψ) in the multi-source region under the constraint. The calculation formula is as follows:

[0089]

[0090] Step 1.4: Based on the initial model, calculate the design variables (Φ h ,Ψ h ),(Φ m ,Ψ m ),(Φ s ,Ψ s ) and define the above-mentioned limited sampling area Ω, such as Figure 2 As shown in the figure, the subscripts h, m, and s are the values ​​at the hub, middle section, and wheel cover, respectively.

[0091] Step 1.5: Calculate the vortex velocity rV based on the above parameters θ (s), the calculation formula is as follows,

[0092] rV θ (s) = r1V θ1 +ξ(s)ΔrV θ

[0093] Among them, r1V θ1 and ΔrV θ Calculated from the original model, ξ(s) is the dimensionless eddy velocity.

[0094] Step 1.6: Calculate the wrap angle distribution along the streamline. In order to construct the geometry of the rotating blade, three sections including the hub, the middle section and the wheel cover are used. Each section is controlled by multiple scattered points. Through joint simplification, the wrap angle of a point on each section is directly calculated, and the results are as follows:

[0095]

[0096] Among them, (r1V θ1 ) h ,(r1V θ1 ) m ,(r1V θ1 ) s are the inlet vortex velocities at the hub, middle and cover respectively; (ΔrV θ ) h ,(ΔrV θ ) m ,(ΔrV θ ) s are the difference between the inlet and outlet vortex velocities; (Φ h,Ψ h ),(Φ m ,Ψ m ),(Φ s ,Ψ s ) are the variables of constraint region control; θ h,j ,θ m,j and θ s,j Respectively represent the wrap angle of the j-th control point of the hub, middle and wheel cover, such as Figure 3 shown.

[0097] The limited sampling areas of the above variables are clarified, as shown in Table 1. At the same time, the objective function to be optimized is determined, as shown in Table 2.

[0098] Table 1 Parameter optimization variables and ranges

[0099]

[0100]

[0101] Table 2 Objective function

[0102]

[0103] Step 2.1: Given the number of sampling points N and the sampling area Ω, use the Poisson disk sampling algorithm combined with the distribution density to accurately control the number of sampling points while maximizing the Poisson disk sampling, and uniformly sample the design variables within the limited area. The sampling diagram is shown in Figure 4 As shown. Randomly sample N points in the region Ω to obtain the initial sampling point set

[0104] Step 2.2: Perform Delaumay triangulation on the initial sampling point set X to calculate the gap area. Specifically, based on the Voronoi diagram theory, locate the gap triangle as a triangle t that satisfies Π(t)>0, and the formula is as follows:

[0105] Π(t)=dis(p,c t )-r>0

[0106] Where p is a vertex of triangle t, r is the sampling radius, c t is the Voronoi center of triangle t, dis(p,c t ) for p to c t The Euclidean distance dis is as follows:

[0107]

[0108] This step is mainly to analyze whether there are gap triangles in the current point set and their specific locations. If there is no gap, it is considered that the sampling point set has been maximized.

[0109] Step 2.3: Optimize the distribution density through an iterative method. In each iteration, find the shortest edge of DT and delete an endpoint v on the edge from X. Regarding the removal strategy of the shortest edge endpoint, the maximum and minimum neighborhood average edge length removal strategy is used, where the neighborhood average edge length of the endpoint is defined as:

[0110]

[0111] Where N(v) is the number of neighboring edges of endpoint v (edges with one endpoint v), l i is the length of the neighborhood edge. Step 2.4: Dynamically update triangulation to detect gap areas. If there is no gap triangle, that is, the sampling area is completely covered, then insert the new sampling point into the center of the largest circumscribed circle and terminate the iteration, otherwise continue the iteration.

[0112] Step 2.5: Perform primitive extraction operations on the gap area to extract the gap area not covered by the sampling disk; randomly generate a sampling point in the gap area extracted in the previous step and insert it into X to dynamically update the triangulation. The above process is repeated until there is no gap area.

[0113] Step 2.6: Obtain data of all variables based on the above work, carry out relevant numerical simulations, obtain objective function results including thrust, power, and cavity size, and then improve the sample library. The objective function distribution is as follows: Figure 5 shown.

[0114] Step 3: Determine the multi-objective constrained optimization equation; construct the WAS agent model; carry out multi-objective optimization work with the help of the improved multi-objective constrained grey wolf optimization algorithm; obtain the Pareto solution set and select the required optimal parameter combination from it.

[0115] Step 3.1: Determine the input variables, output functions, sample libraries and constraints involved in the optimization, and determine the multi-objective optimization equation as follows:

[0116]

[0117] Step 3.2: Initialize algorithm parameters such as population size, maximum number of iterations, and population position; set flag_1=1, flag_2=0, and flag_3=0 in the initial state.

[0118] Step 3.3: Based on the sample library, construct a WAS proxy model and fit the nonlinear mapping relationship between parameters and targets, such as Figure 6 shown.

[0119] Specifically, for all single agent models f i (x) The predicted values ​​are weighted averaged as shown below:

[0120]

[0121] Among them, w i is the proxy model f i The weight of (x) and satisfies Σ i w i =1 and w i ≥ 0. Weight w i The determination of is based on the proxy performance evaluation of the root mean square error, and the corresponding weight allocation strategy is shown as follows:

[0122]

[0123] Among them, ε i The root mean square error of the data points in the test data set. At the same time, the goodness of fit R of the surrogate model is calculated 2 , to evaluate the accuracy of the surrogate model, the formula is as follows:

[0124]

[0125] Among them, Y is the true value, is the predicted value of the surrogate model, is the average of the true values. 2 ≥0.9, the constructed proxy model is considered to meet the accuracy requirements. Figure 7 The distribution diagram of the goodness of fit results is shown in the figure. The results show that all the points basically fall on the diagonal line, that is, the predicted values ​​are basically consistent with the numerical calculation results.

[0126] If the offensive strategy is entered, that is, flag_3 = 1, the cross-validation error is calculated, and the formula is as follows:

[0127]

[0128] Among them, i is the i-th sample point, n i is the number of samples in the validation set, y j is the true value of the sample point, The proxy model fits the sample point. The sample library is screened to build a WAS proxy model with higher accuracy requirements. The screening rules are as follows:

[0129] D t ={x1,...,x d μ-3σ<x i <μ+3σ,i=1,...,d}

[0130] Among them, μ is the mean value of the current sample library, and σ is the standard deviation of the current sample library.

[0131] Step 3.4: Perform non-dominated sorting on all individuals based on the fitted values. The non-dominated partial order sorting under constraints is specified as follows: (1) Any feasible solution has a better non-dominated rank than any infeasible solution; (2) All feasible solutions are clustered according to the objective function value, and the larger the clustering distance, the better the rank; (3) For infeasible solutions, the one with a smaller constraint function violation value is prioritized.

[0132] Step 3.5: Calculate the crowding distance CD of all individuals in the population to maintain the diversity of the population and prevent the algorithm from converging to the local optimal solution too early. The formula is as follows:

[0133]

[0134] Among them, CD im represents the crowding degree of the i-th individual in the m-th objective function, f m represents the mth objective function, x i-1 and x i+1 are the previous individual and the next individual of individual i on target m, respectively, x min and x max They represent the individuals with the minimum and maximum values ​​under the m function respectively.

[0135] Step 3.6: Use the tournament selection method to update the new generation of α, β, and δ leaders in the Pareto frontier according to the non-dominated sorting results and the crowding distance; if the encirclement strategy is entered, that is, flag_2 = 1, then on the basis of the above, it is also necessary to satisfy that the α, β, and δ leaders are all in different clusters. The formula is as follows:

[0136]

[0137] Where k is the number of cluster centers, d is the Euclidean distance between two sample points, the formula is the same as the distance calculation formula in step 2.2, θ is the distance threshold for cluster classification, set to the adaptive change type, the formula is as follows:

[0138]

[0139] Where m and n are the initial value and final value of the adaptive factor, respectively. The value must be in the range [0,1] and m>n; g is the current number of iterations; g max is the maximum number of iterations.

[0140] Step 3.7: Update the positions of other gray wolves according to the positions of the leaders α, β, and δ. Given control parameters such as convergence factor, use the position update formula in the multi-objective gray wolf algorithm; consider the position of the leader and the current position of the gray wolf, and simulate the behavior of the gray wolf tracking prey to help the algorithm shift from exploration to development. Step 3.8; Perform this round of multi-objective optimization work in the complete limited area and update the current Pareto frontier; if flag_1 = 1, further narrow the optimization space and limit the area to the optimization space composed of the leaders α, β, and δ. The formula is as follows:

[0141]

[0142] Among them, Upper is the upper limit of the optimization space, and Floor is the lower limit of the optimization space. In addition, if the search strategy is entered, that is, flag_1 = 1, the optimization goal is set to maximize the agent uncertainty, and the formula is as follows:

[0143] in and are the fitting values ​​of the single agent model respectively; if the encirclement strategy or attack strategy is entered, that is, flag_2=1 or flag_3=1, the optimization goal is set to find the minimum value of the objective function, and the formula is as follows:

[0144]

[0145] Where f is the fitted value of the WAS proxy model.

[0146] Step 3.9: Repeat steps 3.3 to 3.6 to build the WAS proxy model; calculate the fitted values ​​of the Pareto frontier; determine the non-dominated profiles and α, β, δ leaders.

[0147] Step 3.10: According to the position of the α leader, determine the strategy selection corresponding to the current situation. If flag_1=1, set flag_1=0, flag_2=1; if flag_2=1, when the new α leader can dominate the original α leader, set flag_2=0, flag_3=1, otherwise keep flag_2=1; if flag_3=1, when the new α leader can dominate the original α leader, set flag_3=0, flag_1=1, otherwise keep flag_3=1.

[0148] Step 3.11: Determine whether the current iteration meets the maximum number of generations or other termination conditions. If not, jump to step 3.7 and continue to repeat the above process. If it meets the conditions, count all the optimization solutions with excellent performance and output a set of 100 Pareto solutions. The optimal parameter combination is determined, such as Figure 8 shown.

[0149] Step 4: Calculate and test the optimization effect. Generate the corresponding impeller optimal model and perform simulation calculation. Fig. 9 As shown; calculate the performance index and analyze the flow field improvement by comparison. This embodiment provides a comparison between two sets of optimized models with better performance and the initial model, such as Fig.10 shown.

[0150] The final results show that in different cruising speed ranges 0.5U d (30kn)~1.3U d (78kn), the optimized pump-jet propulsor blades help reduce the power consumption of underwater vehicles, while reducing the cavity size and increasing thrust. Among them, the maximum thrust is increased by 6.6%, the maximum power is reduced by 3.1%, and the maximum cavitation size is reduced by 2.67%. Therefore, the effectiveness of the optimization design method has been quantitatively verified.

[0151] At the same time, we also noticed that there are relatively large differences between the original model and the optimized model, which provides more options for further application.

[0152] Although the above embodiments describe the technical solution of the present invention in detail, they are only used to illustrate the technical solution of the present invention, rather than to limit it. In the relevant field, line pipe technicians should understand that the specific implementation of the present invention can still be modified or replaced by equivalents, as long as these modifications or equivalents do not deviate from the spirit and scope of the present invention. Therefore, any modification or equivalent replacement should be included in the scope of protection of the claims of the present invention.

Claims

1. A multi-source constraint design method for high specific speed rotating blades based on improved grey wolf algorithm and Poisson disk sampling, characterized in that: Step 1. Construct a multi-source constraint region to preliminarily screen the range of blade design parameters, including the following sub-steps: Step 1.1: Given high specific speed rotating blade design parameters: The parameters include given high specific speed rotating blade design speed n, design flow rate Q d 、Design head H d , inlet diameter D, calculate the specific speed ns of the rotating blade; determine the blade flow path parameters, set the meridian of the leading edge and trailing edge profiles of the blade, the number of blades, the thickness of the blade, the superposition law, select the wrap angle θ as the only factor to adjust the shape of the rotating blade, and control the shape of the rotating blade by changing the wrap angle distribution along the streamline; Step 1.2: Determine the initial model and calculate the inlet vortex velocity r1V in the initial model simulation data θ1 and the vortex velocity difference ΔrV between the leading and trailing edges θ (s), the formula is defined as follows: Among them, V m1 is the inlet meridian velocity; r1 is the inlet meridian streamline radius; ω is the blade design speed; β b1 is the blade inlet angle; u1 is the circumferential velocity of the inlet meridional streamline; δ is the correction value of the angle; g is the gravitational acceleration; H d is the head under design conditions; η h is the hydraulic efficiency; γ is the correction factor; Step 1.3: Introduce and define a polynomial to represent the dimensionless eddy velocity ξ(s), construct a multi-source constraint region, and convert the calculation formula of the dimensionless eddy velocity into the value of two unknown variables (Φ, Ψ) in the multi-source region under the constraint. The calculation formula is as follows: ξ(s)=Ψs 5 +(-2Ψ+Φ-2)s 4 +(Ψ-Φ+3)s 3 Step 1.4: Based on the initial model, calculate the design variables (Φ h ,Ψ h ),(Φ m ,Ψ m ),(Φ s ,Ψ s ) corresponding to the range, and clarify the above-mentioned limited sampling area Ω; Among them, (Φ h ,Ψ h ),(Φ m ,Ψ m ),(Φ s ,Ψ s ) is the dimensionless eddy velocity control variable in the constraint area of ​​the hub, middle section and wheel cover; Step 1.5: Calculate the vortex velocity rV based on the above parameters θ (s), the calculation formula is as follows, rV θ (s)=r1V θ1 +ξ(s)ΔrV θ Among them, r1V θ1 and ΔrV θ Calculated from the original model, ξ(s) is the dimensionless eddy velocity; Step 1.6: Calculate the wrap angle distribution along the streamline: In order to construct the geometry of the rotating blade, three sections including the hub, the middle section and the wheel cover are used; each section is controlled by multiple scattered points; through joint simplification, the wrap angle of a point on each section is directly calculated, and the results are as follows: Among them, (r1V θ1 ) h ,(r1V θ1 ) m ,(r1V θ1 ) s are the inlet vortex velocities at the hub, middle and cover respectively; (ΔrV θ ) h ,(ΔrV θ ) m ,(ΔrV θ ) s are the differences between the inlet and outlet vortex velocities at the hub, middle and cover, respectively; (Φ h ,Ψ h ),(Φ m ,Ψ m ),(Φ s ,Ψ s ) are the dimensionless eddy current velocity control variables in the constraint areas of the hub, middle and wheel cover respectively; θ h,j ,θ m,j and θ s,j denote the wrap angles of the jth control point of the hub, middle, and wheel cover respectively; Step 1.7: Take (r1V θ1 ) h ,(r1V θ1 ) m ,(r1V θ1 ) s 、(ΔrV θ ) h ,(ΔrV θ ) m ,(ΔrV θ ) s 、(Φ h ,Ψ h ),(Φ m ,Ψ m ),(Φ s ,Ψ s ) is an optimization variable; determining the inlet vortex velocity, the difference variable range between the inlet and outlet vortex velocities based on the initial model simulation calculation results, and determining the dimensionless vortex velocity control variable range based on the multi-source constraint area; Step 2. Given the number of sampling points N and the sampling area Ω, use the Poisson disk sampling algorithm combined with distribution density to accurately control the number of sampling points while maximizing the Poisson disk sampling to uniformly sample the design variables in the limited area, specifically including the following sub-steps; Step 2.1: Randomly sample N points in the region Ω to obtain the initial sampling point set Step 2.2: Perform Delaumay triangulation on the gaps of the sampling point set X. Specifically, based on the Voronoi diagram theory, locate the gap triangle as a triangle t that satisfies Π(t)>0. The formula is as follows: Π(t)=dis(p,c t )-r>0 Where p is a vertex of triangle t, r is the sampling radius, c t is the Voronoi center of triangle t, dis(p,c t ) for p to c t The Euclidean distance d is as follows: Analyze whether there are gap triangles in the current point set and their specific locations; if there is no gap, it is considered that the sampling point set has been maximized; Step 2.3: Find the shortest edge of all triangles in the sampling point set X and delete an endpoint v on the edge; Regarding the removal strategy of the shortest edge endpoint, the maximum and minimum neighborhood average edge length removal strategy is used, where the neighborhood average edge length of the endpoint is defined as: Where N(v) is the number of neighboring edges of endpoint v, l i is the length of the neighborhood edge; Step 2.4: Repeat step 2.2 to analyze whether there is a gap triangle in the sampling point set X to detect the gap area; if there is, extract the primitive, determine the specific position to calculate the gap area, and continue to step 2.5; if not, it is considered that the sampling area is completely covered by the Poisson disk, terminate the loop, and jump to step 2.6; Step 2.5: Execute the primitive extraction operation of the gap area to extract the gap area not covered by the sampling disk; generate a sampling point according to the Poisson disk sampling criterion, insert it into the point set X, and dynamically update the Delaumay triangulation; jump back to step 2.2 for iteration; Step 2.6: It is considered that the sampling radius of the point set X can completely cover the sampling area, so the new sampling point is inserted into the center of the largest circumscribed circle, the iteration is terminated, and the final sampling result is generated; Step 3. Determine the multi-objective optimization equation and construct the WAS agent model; carry out multi-objective optimization with the help of the improved multi-objective constrained grey wolf optimization algorithm; obtain the Pareto solution set and select the required optimal parameter combination from it, including the following sub-steps; Step 3.1: Determine the input variables, output functions, sample libraries and constraints involved in the optimization, and determine the basic settings of the algorithm such as population size and maximum number of iterations; Step 3.2: Initialize algorithm parameters and population positions; set the algorithm to have three different situations: search strategy flag_1, encirclement strategy flag_2, and attack strategy flag_3, and initialize flag_1 = 1, flag_2 = 0, flag_3 = 0; Step 3.3: Based on the sample library, construct a WAS proxy model to generate the fitted values ​​of the population individuals; for all individual proxy models f i (x) The predicted values ​​are weighted averaged as shown below: Among them, w i is the proxy model f i (x) and satisfies ∑ i w i =1 and w i ≥0; Weight w i The determination of is based on the proxy performance evaluation of the root mean square error, and the corresponding weight allocation strategy is shown as follows: Among them, ε i The root mean square error of the data points in the test data set; at the same time, the goodness of fit R of the surrogate model is calculated 2 , to evaluate the accuracy of the surrogate model, the formula is as follows: Among them, Y i is the true value, is the predicted value of the surrogate model, is the average value of the true value; if the calculated R 2 ≥0.9, it is considered that the accuracy of the constructed proxy model meets the requirements; If the offensive strategy is entered, that is, flag_3 = 1, the cross-validation error is calculated, and the formula is as follows: Among them, i is the i-th sample point, n i is the number of samples in the validation set, y j is the true value of the sample point, The proxy model fits the sample points; thus, the sample library is screened to build a WAS proxy model with higher accuracy requirements. The screening rules are as follows: D t ={x1,...,x d |μ-3σ<x i <μ+3σ,i=1,...,d} Among them, μ is the average value of the current sample library, and σ is the standard deviation of the current sample library; Step 3.4: Perform non-dominated sorting on all individuals based on the fitted values; ensure that all individuals are stratified according to the non-dominated relationship, and that the individuals contained in each layer are not dominated by any other individuals in the current layer; Step 3.5: Calculate the crowding distance CD of all individuals in the population to maintain the diversity of the population and prevent the algorithm from converging to the local optimal solution too early. The formula is as follows: Among them, CD im represents the crowding degree of the i-th individual in the m-th objective function, f m represents the mth objective function, x i-1 and x i+1 are the previous individual and the next individual of individual i on target m, respectively, x min and x max They represent the individuals with the minimum and maximum values ​​under the m function respectively; Step 3.6: Use the tournament selection method to update the new generation of α, β, and δ leaders in the Pareto frontier according to the non-dominated sorting results and the crowding distance; if the encirclement strategy is entered, that is, flag_2 = 1, then on the basis of the above, it is also necessary to satisfy that the α, β, and δ leaders are all in different clusters. The formula is as follows: Among them, k is the number of cluster centers, d is the Euclidean distance between two sample points, and θ is the distance threshold of cluster classification, which is set to an adaptive change type. The formula is as follows: Where m and n are the initial value and final value of the adaptive factor, respectively, with a range of [0,1] and m>n; g is the current number of iterations; g max is the maximum number of iterations; Step 3.7: Update the positions of other gray wolves according to the positions of the leaders α, β, and δ; given the control parameters, use the position update formula in the multi-objective gray wolf algorithm; Step 3.8: In the complete limited area, perform this round of multi-objective optimization and update the current Pareto frontier; if flag_1 = 1, further narrow the optimization space and limit the area to the optimization space composed of α, β, and δ leaders. The formula is as follows: Among them, Upper is the upper limit of the optimization space, and Floor is the lower limit of the optimization space; in addition, if the search strategy is entered, that is, flag_1 = 1, the optimization goal is set to maximize the agent uncertainty, and the formula is as follows: in and are the fitting values ​​of the single agent model respectively; if the encirclement strategy or attack strategy is entered, that is, flag_2=1 or flag_3=1, the optimization goal is set to find the minimum value of the objective function, and the formula is as follows: Where f is the fitted value of the WAS proxy model; Step 3.9: Repeat steps 3.3 to 3.6 to build the WAS proxy model; calculate the fitted value of the Pareto frontier; determine the non-dominated profile and α, β, δ leaders; Step 3.10: According to the position of the α leader, determine the strategy selection corresponding to the current situation; if flag_1=1, set flag_1=0, flag_2=1; if flag_2=1, when the new α leader can dominate the original α leader, set flag_2=0, flag_3=1; if flag_3=1, when the new α leader can dominate the original α leader, set flag_3=0, flag_1=1; Step 3.11: Determine whether the current iteration meets the termination condition. If not, jump to step 3.7 and continue to repeat the above process; if it meets the condition, output the current Pareto frontier solution.

2. According to the multi-source constraint design method of high specific speed rotating blades based on improved grey wolf algorithm and Poisson disk sampling as described in claim 1, it is characterized in that: The Poisson disk sampling criterion in step 2 specifically needs to meet the following conditions: (1) Each uncovered point in the sampling area has the same probability of accepting a new sampling point; (2) The distance between any two sampling points is greater than the given sampling radius; (3) The sampling area is completely covered by all sampling disks.

3. The multi-source constraint design method for high specific speed rotating blades based on improved grey wolf algorithm and Poisson disk sampling according to claim 1, characterized in that: In step 3.4, the non-dominated partial order sorting under the constraints is stipulated as follows: (1) Any feasible solution has a higher non-dominated level than any infeasible solution; (2) All feasible solutions calculate the clustering distance according to the objective function value. The larger the clustering distance, the higher the non-dominated level. (3) For infeasible solutions, the solutions with smaller constraint function violation values ​​are prioritized.

Images