Design method of bottom gradual transition lock cut-off dike based on multi-objective navigation flow condition optimization

Abstract

The application discloses a kind of based on the design method of bottom gradual transition waterway separation dike of lock based on multi-objective navigation flow condition optimization, comprising: constructing river reach three-dimensional sublevel mathematical model to obtain the velocity distribution outside dike of proposed location;Accordingly determine the necessity of layout and preliminarily determine the length of structure and average air permeability constraint range;Determine the key design parameter set of waterway separation dike, and define the navigation flow condition evaluation index set containing constraint and target;Based on the above model, numerical simulation is carried out under the constraint of key parameters and the evaluation index value is extracted;Simulation data is used to construct the proxy model connected with parameters and indexes, and a multi-objective optimization algorithm is called to optimize and solve, and the final design parameters are output.The application can overcome the limitations of traditional experience trial-and-error design to some extent, and realize the multi-objective collaborative optimization and fine design of navigation channel and mouth area flow condition under the premise that flow state meets the hard red line requirements of navigation.

Description

Technical Field

[0001] This invention relates to the field of water transport engineering technology, and in particular to a design method for a bottom-gradient permeable lock dike based on multi-objective navigation flow condition optimization. Background Technology

[0002] In hydraulic engineering projects, spillway locks and ship locks are often arranged side-by-side due to topographical constraints. The water flow discharged from spillway locks is typically characterized by large volume, high velocity, and concentrated energy, making its lateral diffusion highly susceptible to encroaching on the approach channels and entrance areas of adjacent ship locks, thus evolving into a complex, unsteady flow field. To ensure the safety of navigation for vessels entering and exiting the locks and the stability of mooring lines for vessels waiting in the approach channels, dikes are usually installed between the ship lock approach channels and the spillway locks. This structure aims to control the lateral flow exchange in a localized area through a physical barrier, reducing the momentum interference of the main waterway on the navigable still water zone, thereby maintaining smooth navigable flow conditions within the waterway.

[0003] Existing dikes mostly employ solid gravity structures, which completely block lateral water interaction, forcing the energy of the downstream flow to concentrate at the dike head for high-speed discharge. This easily generates a strong shear layer near the dike head and leads to a sharp deterioration of the flow pattern in the inlet area. To improve this defect, existing technologies have introduced uniformly permeable dike forms such as pier-type or louver-type dikes. While these permeable solutions alleviate the impact of flow bypass caused by solid dike heads to some extent, specific engineering design often relies on scaled-down physical model tests for empirical trial and error. When dealing with high-velocity impacts and large water level differences, existing permeable structural forms are limited by their single geometric permeability properties, often inducing high-speed jets locally in the pores, and are unable to effectively suppress the large-scale backflow vortices generated at the interface between dynamic and static water due to intense momentum exchange.

[0004] Overall, existing hydraulic diversion structure designs still fall short in matching the spatial structural morphology with fluid energy dissipation mechanisms when dealing with complex hydrodynamic environments involving multiple variables. Especially when facing increasingly stringent navigation safety regulations for pilot channels, traditional trial-and-error models struggle to simultaneously address multiple and physically conflicting requirements, such as lateral velocity control and backflow intensity suppression. Therefore, it is necessary to investigate a systematic approach that can improve the flow adaptability and structural design optimization efficiency of complex hydraulic structures. Summary of the Invention

[0005] The purpose of this invention is to provide a design method for a bottom-gradient permeable lock dike based on multi-objective navigation flow condition optimization, in order to solve the above-mentioned problems in the prior art.

[0006] Technical solution: A design method for a bottom-gradient permeable lock dike based on multi-objective navigation flow condition optimization, comprising: A three-dimensional hierarchical mathematical model of the river section, including the lock, was constructed to obtain the flow velocity distribution outside the dike at the proposed location of the dike. The necessity of setting up the dike is determined based on the flow velocity distribution outside the dike, and the structural length and average permeability constraint range of the dike are initially determined. Determine the key design parameter set for the gradually permeable dike, and define the evaluation index set for assessing the navigation flow conditions of the approach channel and entrance area; Based on the three-dimensional hierarchical mathematical model, numerical simulation is performed under the constraints of the key design parameter set to extract the evaluation index values ​​of navigation flow conditions under each working condition. Based on the evaluation index values ​​extracted from the numerical simulation, a proxy model is constructed to connect the set of key design parameters and the set of evaluation indexes. Based on the aforementioned proxy model, a multi-objective optimization algorithm is invoked to find the optimal solution, and the final design parameters of the bottom-gradient permeable dike that meet the hard constraints of navigation regulations and optimize the navigation flow conditions are output.

[0007] According to one aspect of this application, when the spatial coordinate continuous field function assignment method is used, it specifically includes: Obtain the pre-calibrated analytical expressions for permeability and inertial drag coefficient, and substitute the longitudinal permeability gradient function in the set of key design parameters into the analytical expressions for permeability and inertial drag coefficient to transform them into continuous functions with respect to spatial coordinates. The continuous function relating to spatial coordinates is assigned to the centroids of each computational grid through a custom function, resulting in a spatially continuous and gradual change. A transition layer is introduced in the region where the ratio of drag coefficients between adjacent computational grids exceeds a preset ratio threshold, and a hyperbolic tangent function is used to perform a spatial smooth transition in the transition layer.

[0008] According to one aspect of this application, the finalized design parameters of the bottom-gradient permeable dike that satisfy the hard constraints of navigation regulations and optimize navigable flow conditions specifically include: Obtain the target values ​​of each evaluation index in the Pareto optimal solution set; The target value is then normalized and mapped. Calculate the normalized Euclidean distance from each solution point in the Pareto optimal solution set to the ideal reference point; The design parameters corresponding to the solution set point with the smallest normalized Euclidean distance are selected as the final design parameters.

[0009] According to one aspect of this application, after obtaining the physically graded transparency distribution, the method further includes: Calculate the average permeability value corresponding to the physical gradient permeability distribution; Determine whether the calculated average penetration rate falls within the average penetration rate constraint range; If it does not fall into the trap, the preset maximum permissible lateral inflow velocity is dynamically adjusted according to the deviation direction, and the process returns to the step of using the preset maximum permissible lateral inflow velocity, combined with the water level difference inside and outside the dike and the outflow pattern of the orifice, to reverse-calculate the maximum permeability at each position along the longitudinal direction of the dike, until the average permeability calculation value meets the average permeability constraint range.

[0010] According to one aspect of this application, the extraction method of the reflux zone area and reflux intensity includes: In the surface flow field of the pilot channel and gate area, the connected regions where the longitudinal velocity direction is opposite to the mainstream direction and the absolute value is greater than the preset reverse velocity threshold are identified as the backflow region. The area of ​​the reflow region is obtained by summing the areas of all computational grids within the reflow region. The maximum absolute value of the lateral flow velocity perpendicular to the ship's route direction within the backflow zone is extracted as the backflow intensity.

[0011] According to one aspect of this application, after the output of the finalized design parameters, the method further includes: Based on the established design parameters, a three-dimensional engineering structural model of the bottom-gradient permeable dike is generated, and the corresponding structural design drawings are output.

[0012] The bottom-gradient permeable dike described in this invention refers to a dike structure in which permeable pores or openings are longitudinally arranged in the lower part of the dike structure, and the permeability of the pores is spatially gradually distributed along the longitudinal axis of the dike. This structure differs from traditional solid dikes and uniformly permeable dikes.

[0013] Beneficial effects: This invention reconstructs the design paradigm of hydraulic structures by deeply integrating physical mechanisms, numerical calculations and operations research optimization.

[0014] To address the technical challenge of mismatch between structural morphology and water flow energy dissipation, this solution utilizes the law of lateral momentum balance to reverse-engineer the boundary of permeability distribution along the dike, abandoning subjective and rigid geometric and mathematical assumptions. This allows the infiltration intensity of penetrating water flow to accurately match the spatial attenuation gradient of the dynamic pressure difference outside the dike, thus curbing the generation of local microjets and strong bypass flow around the dike head.

[0015] To address the numerical artifacts that can be caused by complex heterogeneous boundaries in fluid simulation, the drag coefficient is transformed into a spatially continuous field function and a nonlinear smoothing filter is integrated for transition intervention. This eliminates the non-physical pressure step at the interface of the discrete drag grid, improving the fidelity of the numerical flow field from a step approximation to a smooth continuum approximation.

[0016] In addition, in order to overcome the design blind spots caused by the mutual interference of single indicators, this solution constructs a constraint-based hierarchical evolutionary architecture, taking the insurmountable navigation flow velocity red line as a hard criterion for dimensionality reduction and elimination. Within this framework, the Pareto equilibrium point between the scale of the backflow vortex and the lateral disturbance thrust is automatically searched, transforming the experience-based manual trial and error into a global mathematical exact solution that complies with the bottom line of engineering specifications. Attached Figure Description

[0017] Figure 1 This is an overall flowchart of the bottom-gradient permeable dike design method based on multi-objective navigation flow condition optimization in the embodiments of this application.

[0018] Figure 2 This is a graph showing the empirical relationship between the average permeability of the lock dike and the flow velocity outside the dike.

[0019] Figure 3 A flowchart for initially determining the structural length of the dike.

[0020] Figure 4 This is a schematic diagram of the gradually permeable flow barrier at the bottom of the lock. Detailed Implementation

[0021] Example 1: This embodiment provides a design method for a bottom-gradient permeable lock dike based on multi-objective navigation flow condition optimization. The method includes the following steps: Step 101: Construct a three-dimensional layered mathematical model of the river section including the lock to obtain the flow velocity distribution outside the dike at the proposed dike location.

[0022] For water conservancy projects with complex flow patterns and parallel layouts of spillway and lock, a precise three-dimensional physical environment reflecting the riverbed topography and structure boundaries needs to be established in the computer system before conducting navigation flow field optimization design. Specifically, digital topographic elevation data and structure design drawings, including upstream and downstream channels, lock chambers, upstream and downstream approach channels, spillway and gate areas, are acquired. The computational domain is spatially discretized using unstructured hexahedral meshes, and the mesh is locally refined in areas with drastic velocity gradient changes, such as the proposed dike area, approach channel inlet, and gate area, to capture fine turbulent structures. Based on this, the Navier-Stokes equations for incompressible fluids are used as the governing equations, and the RNGk-epsilon turbulence model is employed to close the Reynolds time-averaged equations. The fluid volume function method is used to simulate the water-air two-phase flow interface with a free surface. After the model is established, boundary conditions such as upstream given flow rate boundary and downstream given water level boundary are set, and the hydrodynamic numerical simulation program is run. After the flow field calculation converges and stabilizes, the depth-average velocity vector field of the outer region of the proposed dike is extracted from the simulation results. This vector field characterizes the velocity distribution outside the dike generated by the water discharge from the sluice gate on the approach channel side of the lock, providing dynamic input conditions for the subsequent design of the dike.

[0023] Step 102: Determine the necessity of setting up the dike based on the flow velocity distribution outside the dike, and preliminarily determine the structural length and average permeability constraint range of the dike.

[0024] After obtaining the flow velocity distribution data outside the dike, it is necessary to establish preliminary engineering defense boundaries and physical limits. Specifically, based on engineering experience and standard guidelines, a series of flow velocity thresholds are set. By comparing the relationship between the maximum flow velocity outside the dike and the thresholds through programs or manual comparison, it is determined whether physical isolation facilities are needed, and the extension length of the required facilities. At the same time, since the effect of the dike on water flow is essentially resistance release and energy exchange, it is necessary to set an average range that reflects the overall permeability performance of the dike.

[0025] The average permeability constraint range is set based on the upstream and downstream spatial location of the dike and the corresponding flow velocity outside the dike, and is defined by matching pre-summarized empirical curves of positive and negative correlation in engineering hydrodynamics. Through this step, the solution space can be significantly reduced and the initial structural morphology parameters that meet basic navigation safety can be locked in before high-dimensional parameter optimization.

[0026] For example, based on the obtained flow velocity outside the dike... The following criteria are used for discrimination: (1) If the maximum flow velocity outside the dike is less than 0.5 m / s, then a dike is not required; (2) If the maximum flow velocity outside the dike is between 0.5m / s and 1.0m / s, a short dike is required, and the length of the dike must cover the end of the approach channel mooring section. (3) If the maximum flow velocity outside the dike is greater than 1.0 m / s, a long dike is required, extending to the end of the approach channel braking section.

[0027] In this embodiment, the maximum flow velocity outside the dike is 1.5 m / s > 1.0 m / s, therefore it is determined that a flow-blocking dike is needed. The preliminary plan is to designate the length of the flow-blocking dike from the shielding section to the braking section. It is 360m.

[0028] Since the upstream and downstream dikes have different effects on improving navigation flow conditions, the appropriate permeability ratio needs to be determined based on the location of the dike and the flow velocity outside the dike. Upstream outside the dike flow velocity and optimal permeability The empirical relationship is as follows: when the flow velocity outside the dike is greater than 1.5 m / s, the optimal average permeability is 0.50; when the flow velocity is 1.0 to 1.5 m / s, the optimal average permeability is 0.50 to 0.30; when the flow velocity is 0.5 to 1.0 m / s, the optimal average permeability is 0.30 to 0.15.

[0029] Downstream outside the dike flow velocity and optimal permeability The empirical relationship is as follows: when the flow velocity outside the dike is greater than 1.5 m / s, the optimal average permeability is 0.15; when the flow velocity is 1.0 to 1.5 m / s, the optimal average permeability is 0.25 to 0.15; when the flow velocity is 0.5 to 1.0 m / s, the optimal average permeability is 0.50 to 0.25.

[0030] In some embodiments, the dike is located downstream of the lock, with a maximum flow velocity of 1.5 m / s outside the dike and an initial proposed average permeability of 0.15.

[0031] Step 103: Under the constraints of the structural length and the average permeability range, determine the key design parameter set of the gradually permeable dike, and define the evaluation index set for assessing the navigation flow conditions of the approach channel and the entrance area.

[0032] To achieve a refined and gradually permeable design for the dike, a parametric model was established, abstracting the dike structure into a set of key design parameters. This parameter set mainly includes the length of the permeable section, which determines the start and end positions of the gradual transition section, and the longitudinal permeability gradient function, which determines the longitudinal water flow infiltration intensity.

[0033] Furthermore, to measure the actual improvement effect of the dike on navigation flow, a set of mathematical evaluation indicators needs to be established. This indicator set extracts key parameters affecting the stability of ship navigation and berthing from the flow field, such as the maximum lateral velocity at the entrance area, which characterizes the strength of the crossflow; the area and intensity of the backflow zone, which characterize the extent of unfavorable flow conditions; and the maximum longitudinal velocity, which is strictly limited by national navigation design specifications. These indicators will serve as mathematical criteria for subsequent multi-objective optimization.

[0034] For example, the process is as follows: (1) Extract the key design parameter set of the gradually permeable dike, including: Length of open section The length from the head of the dike to the starting position of the permeable section is typically 100 to 300 meters.

[0035] Longitudinal permeability gradient function: a linear function or an exponential function is used, which will be described in detail later.

[0036] (2) Define the evaluation index set for navigable water flow conditions: The evaluation indicators for navigation flow conditions in the upstream and downstream approach channels and entrance areas of the lock include: maximum longitudinal and transverse flow velocities in the entrance area, maximum longitudinal and transverse flow velocities in the approach channels, and the intensity and extent of backflow.

[0037] Maximum lateral velocity in the entrance region Within the entrance area, the maximum transverse velocity should be calculated. Upstream, it should meet the specification requirement of ≤0.3~0.5m / s, and downstream, it should meet the specification requirement of ≤0.5m / s.

[0038] Maximum longitudinal velocity in the entrance region The maximum longitudinal velocity within the entrance area should meet the specification requirement of ≤2.0 m / s.

[0039] Maximum longitudinal velocity of the pilot channel The maximum longitudinal flow velocity in the pilot channel (berthing section and braking section) should be calculated. The upstream should meet the standard requirement of ≤0.5m / s, and the downstream should meet the standard requirement of ≤1.0m / s.

[0040] Maximum lateral velocity of the pilot channel The maximum lateral velocity in the pilot channel (berthing section and navigation section) should be ≤0.15m / s upstream and ≤0.20m / s downstream.

[0041] Reflux area : Refers to the length and width of the backflow formed in the entrance area or pilot channel, in meters.

[0042] Reflux intensity The flow velocity should be perpendicular to the direction of the ship's route, and should be less than 0.25 m / s in the entrance area and less than 0.10 m / s in the pilot channel.

[0043] Step 104: Based on the three-dimensional hierarchical mathematical model, numerical simulation is performed under the constraints of the key design parameter set to extract the evaluation index values ​​of navigation flow conditions under each working condition.

[0044] After determining the design variables and evaluation indicators, an adaptive design of experiments method is used to generate a series of discrete working condition sample points within the parameter space defined by the key design parameter set. For each generated parameter combination sample, its corresponding geometric and drag source parameters are substituted into the three-dimensional hierarchical mathematical model established in the preceding steps.

[0045] A drag source term characterizing the momentum loss of the porous medium is loaded into the corresponding permeable region of the computational grid, and the fluid dynamics numerical simulation is run again. For each steady flow simulation case, after the flow field reaches steady state, the post-processing module is called to calculate and extract the corresponding navigation flow condition evaluation index values ​​of the surface of the pilot channel and the gate area, forming a one-to-one mapping sample pair between input parameters and output flow field effects.

[0046] Step 105: Based on the evaluation index values ​​extracted from the numerical simulation, construct a proxy model connecting the set of key design parameters and the set of evaluation indexes.

[0047] Because the computational cost of directly using the CFD three-dimensional hierarchical mathematical model for iterative optimization of evolutionary algorithms is too high, surrogate model techniques such as response surface modeling are used to replace high-fidelity fluid numerical simulation.

[0048] The surrogate model can be a second-order polynomial response surface model, a radial basis function model, or a Kriging model, etc. In this embodiment, a second-order polynomial response surface model is used for illustration.

[0049] The response surface model takes the following form: (1) In the formula: Let be the objective function. For the proposed design parameters, in this embodiment, the design parameter is the length of the open section. and initial permeability : , , , The regression coefficients to be estimated are: This is random error; This represents the total number of design parameters. and A numbered index for different design parameters.

[0050] Specifically, the parameter input samples and their corresponding evaluation index output samples generated under the above working conditions are collected. Regression analysis methods such as second-order polynomial response surface fitting are used, and the regression coefficients are estimated using the least squares method, thereby establishing a proxy model function relationship between the input independent variables and the output dependent variables.

[0051] In this process, methods such as cross-validation can be used to verify the prediction accuracy of the surrogate model, ensuring that it can accurately approximate the real flow field response surface while significantly reducing computation time.

[0052] Step 106: Based on the surrogate model, call the multi-objective optimization algorithm to find the optimal solution and output the final design parameters of the bottom gradually permeable dike that meet the hard constraints of navigation specifications and optimize the navigation flow conditions.

[0053] After the surrogate model is constructed, it is embedded into a multi-objective optimization framework, such as a non-dominated sorting genetic algorithm, for high-frequency iteration. During the algorithm optimization process, on the one hand, an evaluation system with constraint handling mechanisms eliminates infeasible solutions that violate navigation regulations and red-line indicators such as longitudinal flow velocity; on the other hand, a Pareto optimal solution set is sought among various mutually constraining soft objective indicators. After the algorithm converges, a normalized multidimensional spatial distance metric is applied to automatically extract a design parameter combination from the Pareto front surface that achieves the best trade-off among the various evaluation indicators. This set of parameters serves as the final design parameters for the bottom-gradient permeable dike.

[0054] The multi-objective optimization algorithm can be NSGA-II, MOEA / D, or SPEA2, etc. In this embodiment, NSGA-II with constraint processing mechanism is used as an example.

[0055] Furthermore, in an optional implementation, after outputting the finalized design parameters, a subsequent step of generating an engineering deliverable is included. Specifically, based on the structural length, open section length, and specific gradual distribution law of the openness ratio contained in the finalized design parameters, these are imported into computer-aided design software or a building information modeling system to generate a three-dimensional engineering structural model of a bottom-gradient openwork dike with specific pier dimensions, opening arrangement, and foundation shape. Based on this model, two-dimensional or three-dimensional structural design drawings are output for construction and manufacturing reference, thereby completing the technological transformation from theoretical algorithm design to physical engineering implementation.

[0056] For example, a porous medium model is used to simulate the gradually changing permeable section: When the pores are dense or the mesh resolution is limited, the permeable section is treated as a porous medium region, and the flow resistance is simulated by adding a momentum source term. The relationship between the permeability and the drag coefficient is calibrated through independent small-scale element flow simulations.

[0057] Select the length of the open section and initial permeability The design parameters are variables. In this example, the length of the permeable section is 120~300m, the initial permeability is 0.07~0.11, a total of 4 permeable section length comparison schemes and 3 initial permeability comparison schemes are set, and a total of 12 test conditions are combined.

[0058] Length of open section : 120m, 180m, 240m, 300m; Initial permeability : 0.07, 0.09, 0.11; Steady flow simulation was performed for each test condition, and the evaluation index of the stabilized flow field was extracted.

[0059] Example 2: Based on the above embodiment one, this embodiment further details the specific steps of determining the layout of the dike, the initial calculation of structural parameters, and the digital extraction of navigation evaluation indicators.

[0060] Step 201: Extract the maximum flow velocity outside the dike from the flow velocity distribution outside the dike.

[0061] After obtaining the velocity distribution outside the proposed dike location, the velocity field data of grid nodes or unit centroids in the computational domain of the three-dimensional layered mathematical model are traversed to calculate and extract the peak velocity to obtain the maximum velocity outside the dike. By obtaining this peak value, the lateral kinetic energy disturbance amplitude of the water discharge from the spillway to the navigation area of ​​the lock can be quantified. After extracting the maximum velocity outside the dike, a comparison is performed with the pre-configured velocity control benchmark to determine the necessity of the physical isolation structure's engineering layout. It should be noted that the structural length mentioned here is specifically limited to the physical sum of the length of the solid section and the length of the open section.

[0062] Under the first working condition, if the maximum flow velocity outside the dike is less than the first preset flow velocity threshold, it is determined that no dike needs to be set up, and the determination result that no dike needs to be set up is output.

[0063] Specifically, the first preset flow velocity threshold can be set to 0.5 m / s. When the peak lateral flow velocity is lower than this set value, the lateral momentum transfer energy of the water flow is low and insufficient to interfere with the ship's entry and exit from the lock or its mooring attitude. The water space of the pilot channel itself has sufficient water flow energy dissipation capacity, so a judgment command is output that no additional physical isolation structure is required.

[0064] Under the second operating condition, if the maximum flow velocity outside the dike is between the first preset flow velocity threshold and the second preset flow velocity threshold, it is determined that a short flow barrier should be set, and the structural length is set to the length extending from the tail of the guideway mooring section; specifically, the preset second preset flow velocity threshold can be set to 1.0 m / s. Within this range, the lateral thrust generated by the downstream crossflow reaches the action threshold, requiring the generation of parameter commands to construct a physical barrier. The extended structural length output value is bound to the coordinates of the tail of the guideway mooring section to cut off the lateral water flow exchange in this section.

[0065] Under the third operating condition, if the maximum flow velocity outside the dike is greater than the second preset flow velocity threshold, it is determined that a long dike should be set, and the structural length is set to the length extending from the tail of the pilot channel braking section. Specifically, when the peak lateral flow velocity exceeds the second preset flow velocity threshold, the high-kinetic-energy water flow will diffuse outward and intrude into the ship braking area, reducing the ship's heading control performance during deceleration and lock entry. At this time, the determination is updated, the length of the physical defense line is increased, and the output structural length endpoint coordinates are extended to cover the tail of the braking section, generating a longer breakwater boundary.

[0066] Step 202: Identify the spatial location of the proposed dike upstream or downstream of the lock.

[0067] After initially determining the structural length, the model coordinate parameters with geographical attributes in the building design drawings are read to identify the specific relative position of the proposed structure from the hydraulic intake or outlet side of the lock head, thus completing the digital assignment of the upstream and downstream spatial location attributes. The motivation for distinguishing this location attribute lies in the systematic differences in the hydrodynamic boundary conditions of the upstream and downstream approach channels, which exhibit diametrically opposed physical mechanisms in their physical requirements for permeability.

[0068] Step 203: Based on the spatial location and the flow velocity distribution outside the dike, match the pre-stored empirical relationship between flow velocity and permeability to determine the average permeability constraint range.

[0069] Before invoking the pre-stored empirical relationship between flow velocity and permeability, this empirical relationship is pre-constructed based on historical engineering data. Specifically, measured flow field data from different lock engineering projects are collected, and a series of discrete outflow velocity variables are input into the standardized generalized model. By permuting and combining these variables, the flow regime under different uniform permeability is scanned. Data points of permeability when the crossflow in the outlet gate area is below the specification limit and the backflow intensity reaches its minimum are selected. A continuous functional relationship curve is fitted and pre-stored in the database. Based on the location labels obtained in step 202, corresponding branch matching is performed.

[0070] The pre-stored empirical relationship between flow velocity and permeability is specifically represented as a continuous function curve with the maximum flow velocity outside the dike as the independent variable and the recommended average permeability as the dependent variable.

[0071] In one alternative implementation, the empirical relationship can be fitted using a first-order linear function or a power function. The upper and lower bounds of the average permeability constraint range are generated by adjusting the recommended permeability curve values ​​upward and downward by a certain percentage (e.g., ±10% to ±20%), respectively, to allow for parameter margins in subsequent optimization searches.

[0072] Those skilled in the art can select the specific function form and fluctuation bandwidth based on the amount of available historical engineering data and the fitting accuracy.

[0073] Under the first spatial location condition, when the location is upstream of the lock, the constraint boundary of the average permeability constraint range is set based on the empirical relationship that the outward flow velocity is positively correlated with the optimal permeability. Specifically, when the ship is in the upstream approach channel entering the lock against the current, the lateral impact kinetic energy of the main flow is highly correlated with the water level difference between the inside and outside of the lock. When a high flow velocity parameter is input, a larger permeability constraint parameter is extracted based on the positive correlation curve mapping, thereby allowing the high-velocity water flow to gradually diffuse and reduce the water level difference between the inside and outside of the lock, thus suppressing the phenomenon of crossflow concentration at the lock head.

[0074] Under the second spatial location condition, when the location is downstream of the lock, the constraint boundary of the average permeability constraint range is set using the empirical relationship that the outflow velocity outside the dike is negatively correlated with the optimal permeability. Specifically, when the ship is downstream and exiting the lock in a downstream manner, the main source of disturbance becomes the transverse backflow vortex generated at the interface between dynamic and static water. When high flow velocity parameters are input, a smaller permeability constraint parameter is extracted through the negative correlation curve to enhance the isolation and guiding effect on the high-speed discharge flow, thereby avoiding the entrainment effect of high-speed fluid on the still water zone of the pilot channel.

[0075] Step 204: Divide the evaluation index set into hard constraint indexes and target indexes.

[0076] The evaluation index set established in the early stage is classified and parsed. Based on the different numerical judgment logic, an insurmountable rule set and a floating evaluation rule set are established in memory respectively, so as to form an input set that supports different penalty operation logics of subsequent hierarchical algorithms.

[0077] Step 205, the hard constraint indicators include the maximum longitudinal flow velocity in the entrance area and the maximum longitudinal flow velocity in the pilot channel, which have clearly defined limits in the navigation regulations.

[0078] The hydraulic limits with non-negotiable characteristics are categorized into a hard constraint array. The longitudinal velocity component matrix within a specified region of the computational domain is extracted to assess whether it exceeds the flow velocity threshold required by navigation regulations. In subsequent surrogate model calculations, variable combinations that violate these indicators will be directly assigned a penal value with exclusionary properties to forcibly eliminate them.

[0079] Step 206, the target indexes include the maximum transverse velocity in the gate region, the backflow intensity, and the backflow area, which characterize the degree of flow field deterioration and need to be minimized.

[0080] These indicators also originate from vector or scalar extraction of the three-dimensional flow field, but they possess complex nonlinear constraint mapping relationships. Suppressing one parameter often causes a rebound increase in another parameter. Incorporating these variables into the target indicator sequence serves primarily as a fitness benchmark function for calculating individuals in the multi-objective algorithm's optimization evolution, driving the optimization operator to converge in the feasible solution space towards the direction where the values ​​of each indicator tend to the lower bound.

[0081] According to one aspect of this application, the extraction method of the reflux zone area and reflux intensity includes the following steps: Step 207: In the surface flow field of the pilot channel and gate area, identify the connected regions where the longitudinal flow velocity direction is opposite to the mainstream direction and the absolute value is greater than the preset reverse flow velocity threshold as the backflow zone.

[0082] To address the issue of extracting backflow field features from the aforementioned target indicators, a surface mesh flow field traversal is performed. The mainstream forward direction is set as the positive vector coordinate, and a preset reverse velocity threshold can be set to 0.02 m / s. When the longitudinal velocity scalar of a mesh centroid meets the criterion of a velocity less than -0.02 m / s, that mesh is pushed into the target set. Further, spatially adjacent meshes within the target set are clustered and merged to form a complete, independent connected domain structure, thereby completing the digital identification of the backflow region.

[0083] Step 208: Sum the areas of all computational grids within the reflow region to obtain the area of ​​the reflow region.

[0084] After defining the independent connected regions, the grid geometric topology information attribute table is retrieved. Using the discrete area integral processing rules, the scalar area values ​​of each element belonging to the recirculation zone set on the projected plane are accumulated and summarized, outputting single-value data representing the geometric span of the recirculation deterioration region. Specifically, this is calculated using the following formula: (2) in, Let be the area of ​​the recirculation region, and ∑ be the area of ​​the th connection within the set of connected components. Instructions for summing elements For the first The horizontal projected area parameter of a fluid grid cell that satisfies a preset reverse flow velocity condition.

[0085] Step 209: Extract the maximum absolute value of the lateral flow velocity perpendicular to the ship's route direction within the backflow zone, and use it as the backflow intensity.

[0086] To quantify the magnitude of the destructive kinetic energy carried by the backflow, an extreme value search is performed on the components of the transverse dimension within the extraction matrix, based on the locked node indices of the connected regions. The magnitude of the transverse component is compared for each discrete data point in the set, and the node data with the highest absolute value is selected for output. This is used to measure the strongest rotational or lateral deflection thrust that the flow field may exert on the ship's sidewalls. Specifically, it is extracted using the following formula: (3) in, For the reflux intensity, max is the function operator that searches for and extracts the maximum value within the target set. To eliminate absolute value function operators that retain only their magnitude and direction of sign, This is the lateral velocity parameter extracted from the j-th grid point within the set, perpendicular to the main flow direction.

[0087] Example 3: There are multiple definite paths for the longitudinal permeability gradient function. This example provides a basic alternative based on the mathematical curve preset method.

[0088] Step 301: The longitudinal permeability gradient function is determined by a mathematical curve preset method or a physical driving back-calculation method based on the fluid lateral momentum balance.

[0089] To address the spatial distribution of water permeability intensity within the permeable section, a functional expression describing the variation of the permeability ratio along the axis of the dike is required. This expression is constructed using two parallel processing logics. In the first logic, a pre-configured monotonically increasing mathematical curve equation is invoked for spatial allocation of geometric parameters. In the second logic, the dynamic water pressure difference generated by the water flow around the dike is introduced as a mechanical boundary condition for analytical derivation. This embodiment focuses on a detailed explanation of the technical actions of the first execution path using the mathematical curve pre-set method.

[0090] Step 302: When the mathematical curve preset method is used, the longitudinal permeability gradient function adopts a linear increasing function or an exponential increasing function, and after the permeability segment is longitudinally divided equally, the corresponding discrete permeability average value is assigned.

[0091] Without delving into the underlying fluid dynamics, a one-dimensional spatial mapping domain is established in memory with the origin of the open section as the coordinate origin. Based on the call command from the operation input, a basic increasing function model is extracted from the candidate curve library to characterize the variation of the open section rate along the path.

[0092] Under one operating condition, a linearly increasing function is used to construct a mapping relationship. This mapping relationship assumes that the permeability ratio increases proportionally with increasing distance at a fixed slope, specifically calculated using the following formula: ； (4) in, Distance from the starting point of the transparent section Local permeability at that location The initial penetration rate is set at the starting point of the open section. The spatial linear growth rate of the permeability The coordinates are the spatial coordinates of the distance from the starting point of the open section along the longitudinal axis of the dike.

[0093] To meet the overall permeability requirements initially determined, a definite integral was performed on the constructed continuous linear equation over the length of the permeable section, and the result was divided by the total length to obtain the expected average value of the function. By matching this expected average value with the obtained average permeability constraint range, the specific range of values ​​for the spatial linear growth rate was obtained through inverse constraint and solution, thus completing the full parameter determination of the linear function.

[0094] Under another operating condition, an exponentially increasing function is used to construct a mapping relationship. This mapping relationship assumes that the permeability ratio increases relatively slowly near the initial solid section, while it accelerates and amplifies near the outer opening area. Specifically, this is calculated using the following formula: (5) in, Distance from the starting point of the transparent section Local permeability at that location The initial penetration rate is set at the starting point of the open section. It is a natural constant. This is the spatial exponential growth coefficient of the permeability, and its dimension is the reciprocal of length. The coordinates are the spatial coordinates of the distance from the starting point of the open section along the longitudinal axis of the dike.

[0095] Similarly, the mean of the definite integral over a specified interval is calculated for the given continuous exponential equation, and this mean is used as the basis to calculate the numerical range of the spatial exponential growth coefficient.

[0096] After constructing the expression and solving the coefficients of the continuous longitudinal permeability gradient function, dimensionality reduction discretization is performed to adapt to the spatial discretization characteristics of the fluid mesh in subsequent numerical simulations. The entire length interval is divided into a predetermined number of uniform spatial segments along the longitudinal extension axis of the permeable segment. For each independent uniform spatial segment, the local integral average value of the segment's start and end coordinates on the continuously increasing function is calculated. This value is extracted as the discrete permeability average value of the corresponding segment and used as a constant geometric attribute label, assigned and bound to all computational grid cells falling within the spatial boundary of that segment.

[0097] Example 4 Unlike Embodiment 3 above, to address the problem of localized jetting of water flow caused by the lack of fluid dynamics basis in the preset mathematical function, this embodiment provides a physical-driven reverse-engineering alternative based on the fluid's lateral momentum balance. This embodiment couples the generation process of the permeability distribution with the fluid dynamics boundary conditions for derivation, and combines a closed-loop iterative mechanism to obtain parameter outputs that meet engineering constraints.

[0098] Step 401: Extract the transverse velocity distribution along the longitudinal direction of the dike in the velocity distribution outside the dike.

[0099] After obtaining the complete external velocity vector field at the proposed dike location, a spatial decomposition operation is performed on this vector field. The lateral velocity scalar component perpendicular to the predetermined axis of the dike is extracted along its longitudinal extension trajectory. This lateral velocity component reflects the physical driving intensity of the lateral momentum transfer generated by the water flow from the spillway towards the lock's approach channel. The attenuation or fluctuation pattern of its value along the path constitutes the basic input conditions for subsequent hydrodynamic back-calculation.

[0100] Step 402: The dynamic pressure difference generated by the transverse velocity distribution is equivalent to the water level difference between the inside and outside of the dike.

[0101] When high-speed water flow impacts a dike structure with physical shielding properties, the kinetic energy of the water is partially converted into potential energy. Based on the fluid dynamics energy equation, it is assumed that the dynamic pressure difference generated by the velocity field on the outside of the structure is the main physical driving force for the water flow to penetrate the dike's pores and form the difference in water level inside and outside. Based on this, an equivalent transformation relationship between the velocity field and the water level boundary is established, specifically calculated using the following formula: (6) in, Distance from the starting point of the transparent section The equivalent value of the water level difference between the inside and outside of the dike at that location. Distance from the starting point of the transparent section The transverse velocity component extracted outside the dike, This is the acceleration due to gravity.

[0102] Step 403: Using the preset maximum permissible lateral inflow velocity, combined with the water level difference inside and outside the dike and the outflow pattern of the orifice, the maximum permeability at each position along the longitudinal direction of the dike is calculated.

[0103] After obtaining the equivalent value of the water level difference, the physical process of water flowing through the pores of the dike is abstracted into an orifice outflow model. According to the orifice outflow dynamics, the inflow velocity penetrating the dike and entering the inner side of the navigation channel depends on the external driving pressure difference and the structural permeability.

[0104] According to the classical orifice outflow dynamics formula, the outflow velocity of water passing through the perforated pores of the dike can be expressed as: (7) in, Distance from the starting point of the transparent section The outflow velocity at the orifice, For flow coefficient, It is the acceleration due to gravity. This refers to the equivalent value of the water level difference between the inside and outside of the dike calculated in step 402. The equivalent relationship established in step 402 is then applied. Substituting into the above equation, we get: (8) Furthermore, the water flows through a region with high porosity. After the structure is completed, the apparent inflow lateral velocity inside the pilot channel is the product of the orifice outflow velocity and the permeability: (9) in, The lateral velocity of the inflow on the inner side of the dike. For flow coefficient, The structural permeability at the current location. C represents the lateral velocity component outside the dike. d The specific geometric shape of the perforated opening of the dike can be determined through hydraulic tests or by consulting hydraulic manuals. Those skilled in the art can select appropriate values ​​based on the actual shape of the opening.

[0105] A physical limit threshold is set for the inflow velocity. The initial value of the maximum permissible lateral inflow velocity can be set according to the limit requirements for lateral flow velocity in the navigation regulations of the pilot channel. Those skilled in the art can select an appropriate initial value based on the navigation level and ship size of the actual project.

[0106] In one optional implementation, the preset maximum permissible lateral inflow velocity can be determined based on the control requirements for lateral flow in the navigation channel to ensure safe navigation. Its specific value can be set with reference to the lateral flow velocity limits in the relevant navigation regulations. For example, the initial value of the maximum permissible lateral inflow velocity can be set to 0.20 m / s. The lateral inflow velocity v on the inner side of the dike in the formula... in (x) is replaced with the preset maximum permissible lateral inflow velocity. By rearranging and solving the above equation, we obtain the following inverse formula: (10) in, To calculate the maximum permissible penetration rate, The preset maximum permissible lateral inflow velocity, For flow coefficient, This represents the lateral velocity component outside the dike.

[0107] In one alternative implementation, The value is 0.62.

[0108] Step 404: The maximum allowable permeability is truncated using a preset upper limit for structural safety permeability and a preset lower limit for minimum achievable permeability to obtain a physically gradient permeability distribution.

[0109] After performing the mathematical back-calculation, the output permeability value may exceed the feasible range of the actual engineering structure. Excessive permeability can lead to insufficient strength of the dike pier structure, while insufficient permeability may be limited by the precision of underwater construction operations or result in high-speed microjet flow generated through extremely narrow pores. Therefore, the upper limit of the structurally safe permeability is determined by the load-bearing capacity calculation of the dike pier structure, and the lower limit of the minimum achievable permeability is determined by the minimum opening size limitation of the underwater construction process, such as based on the precision of underwater construction operations and the minimum clear width of the orifice. Both can be determined by those skilled in the art based on the specific structural form and construction conditions.

[0110] In one optional implementation, the upper limit of the structural safety permeability is set to 0.60, and the lower limit of the minimum achievable permeability is set to 0.05.

[0111] A numerical truncation operation is performed by introducing hardware constraint boundaries, specifically handled by the following formula: (11) in, This is the physically gradient transparency distribution output after truncation. For the minimum value function operator, For the maximum value function operator, For the pre-configured upper limit parameter of structural safety permeability, This is the pre-configured minimum achievable permeability lower limit parameter. This is the maximum allowable permeability obtained by working backward from the previous steps.

[0112] This operation forces the removal of values ​​exceeding the limit, outputting an engineered distribution curve that conforms to both physical laws and manufacturing process constraints.

[0113] Step 405: Perform polynomial fitting on the physical gradient permeability distribution, and extract the polynomial coefficients after fitting as design variables characterizing the longitudinal permeability gradient function.

[0114] To transform the continuous or piecewise nonlinear physical gradient transparency distribution into a dimensionality-reduced data structure that can be efficiently processed by multi-objective optimization algorithms, numerical regression tools such as the least squares method are used to perform a second-order mathematical mapping on the truncated curve. Specifically, this is calculated using the following formula: (12) in, To fit the output longitudinal permeability gradient function, The coefficient of the constant term, The coefficient of the linear term, The coefficient of the quadratic term, These are the spatial coordinates along the longitudinal axis of the dike. After fitting, extract the coordinates containing... , , The multivariate coefficient matrix, including the original spatial distribution sequence, is used as the independent variable for subsequent optimization calculations and iterative updates.

[0115] Assuming the transverse velocity outside the dike decreases along the longitudinal direction of the dike, with a higher velocity near the spillway outlet and decreasing to a lower value further away, after the equivalent transformation in step 402, the water level difference between the inside and outside of the dike shows a similar decreasing trend along the flow path. After back-calculation in step 403, the maximum allowable permeability is lower in the region with higher transverse velocity, indicating a more closed structure to block high-speed flow; conversely, the maximum allowable permeability is higher in the region with lower transverse velocity, indicating a more open structure to allow low-speed flow to infiltrate and diffuse. This flow path distribution characteristic allows the spatial configuration of permeability to adaptively match the attenuation gradient of the external dynamic pressure difference.

[0116] Step 406: Calculate the average permeability value corresponding to the physical gradient permeability distribution.

[0117] After generating the permeability distribution curve with physical driving characteristics, its overall permeability performance needs to be macroscopically verified. A definite integral summation is performed on the obtained gradient function within the domain of the total length of the permeable section. The integral result is then divided by the total length of the permeable section, and the expected mathematical value of this distribution over the entire interval is output. This value serves as a macroscopic representation of the overall permeability.

[0118] Step 407: Determine whether the calculated average permeability falls within the average permeability constraint range.

[0119] The average permeability constraint range generated in the previous steps based on the matching of empirical curves of upstream and downstream spatial location and flow velocity is retrieved. The calculated average permeability value is then compared with the upper and lower boundary values ​​of this constraint range. This comparison is used to verify whether the physical curve derived from the local lateral momentum balance is compatible with the global wave protection empirical criteria.

[0120] Step 408: If it does not fall in, dynamically adjust the preset maximum allowable lateral inflow velocity according to the deviation direction, and return to execute the step of using the preset maximum allowable lateral inflow velocity, combined with the water level difference inside and outside the dike and the outflow law of the orifice, to back-calculate the maximum allowable permeability at each position along the longitudinal direction of the dike, until the average permeability calculation value meets the average permeability constraint range.

[0121] If the comparison and judgment output is negative, a closed-loop adjustment mechanism is triggered. If the calculated average permeability is higher than the upper limit of the constraint range, it indicates that the overall porosity is too large, and the preset maximum allowable lateral inflow velocity is decreased downward by a set numerical step size; conversely, if the calculated average permeability is lower than the lower limit of the constraint range, it indicates that the overall structure is too closed, and the preset maximum allowable lateral inflow velocity is increased upward by the same or adaptive step size. After the parameter update is completed, the program instruction jumps to step 403 to reload the data link for back-calculation, truncation, and integration calculations. This iterative process continues to run in the background until the updated calculated average permeability is within the constraint closed interval, at which point the iteration terminates and the output polynomial coefficients at this time are locked, completing the closed-loop generation of design parameters.

[0122] Example 5 After obtaining the design parameter set output by the aforementioned embodiments, this embodiment provides a basic discrete segmented simulation scheme for the generation of sample points for numerical simulation working conditions and the resistance loading mechanism in Embodiment 1, and details the offline calibration process of the resistance parameters of porous media.

[0123] Step 501: Obtain the dimension of the decision variables of the key design parameter set.

[0124] Before initiating batch numerical simulations of fluid mechanics, the incoming key design parameter set data structure is read, and the total number of variables that serve as independent optimization objects is calculated and extracted. This total number of variables represents the dimensionality of the parameter space that the subsequent optimization algorithm needs to explore, serving as the underlying basis for the subsequent adaptive selection of experimental sampling strategies.

[0125] In the first execution path, in step 502, if the dimension of the decision variable is less than or equal to a preset dimension threshold, then the full factorial design method is used to generate the numerical simulation working condition sample points.

[0126] The obtained dimension of the decision variables is compared with a preset dimension threshold. The preset dimension threshold can be set to 3. When the dimension of the decision variables is in the low-dimensional range, the combinatorial explosion effect of the parameter space is controllable. The full factorial design program is invoked to perform a complete orthogonal combination for each discrete level of each design variable. By traversing all possible parameter nodes, a set of numerical simulation case sample points arranged in a regular pattern is generated, thereby achieving complete coverage of the low-dimensional parameter space.

[0127] In the second execution path, in step 503, if the dimension of the decision variable is greater than the preset dimension threshold, then the Latin hypercube sampling method is used to generate the numerical simulation working condition sample points.

[0128] When the dimensionality of the decision variables exceeds a set threshold, using a full factorial design leads to an exponential increase in computational complexity. In this case, the sampling engine is automatically switched to Latin hypercube sampling. This method divides the value range of each dimension into multiple sub-intervals with equal probability and performs random sampling within each sub-interval, ensuring that the generated working condition sample points have good uniform filling characteristics in the multidimensional high-dimensional space. In this way, highly representative numerical simulation working condition sample points are extracted while maintaining limited computational resource consumption.

[0129] Step 504: A porous medium model is used to simulate the bottom gradually transparent structure, wherein the momentum drag source term is applied to the computational grid by a piecewise uniform assignment method or a spatial coordinate continuous field function assignment method.

[0130] When constructing a three-dimensional layered mathematical model of a dike containing complex perforations, the mesh resolution makes it difficult to directly characterize each physical opening through geometric modeling. A porous media physical model is introduced, abstracting the perforated structure as a computational domain that exerts directional resistance on the fluid. The momentum loss process of water flowing through the permeable structure is simulated by adding momentum resistance source terms to the fluid momentum equation. The loading strategies for momentum resistance source terms are divided into two categories: piecewise uniform assignment and spatial coordinate continuous field function assignment.

[0131] Before applying momentum drag source terms to the computational grid, the drag characteristic baseline parameters of the corresponding structure must be obtained through offline experiments. A cell calibration model is constructed based on the Darcy-Forchheimer law, establishing a local porous medium channel containing inlet and outlet regions. Multiple sets of permeability boundaries with different gradients are set, and multiple gradients of uniform inflow velocities are assigned to each set of permeability to drive multiple steady-flow simulations.

[0132] Pressure drop data for the upstream and downstream of the permeable structure under each simulation condition are extracted, and the pressure drop gradient per unit length is calculated. Then, with apparent velocity as the independent variable and pressure drop gradient per unit length as the dependent variable, a quadratic polynomial fitting operation is performed using the numerical analysis module. The fitting process is described by the following formula: (13) in, For the pressure drop upstream and downstream of the perforated structure, The longitudinal thickness of the porous medium region. The dynamic viscosity of the fluid. For penetration rate, The apparent velocity through the porous medium, The inertial drag coefficient, The fluid density is given.

[0133] The coefficients of the first and second terms in the above formula are obtained using a regression algorithm, and then the discrete values ​​of permeability and inertial drag coefficient at various permeability rates are extracted. Furthermore, these discrete values ​​are mapped and fitted to the structural permeability variable to establish a continuous analytical expression. Specifically, the constant coefficients are extracted using the following formula: (14) (15) in, For penetration rate, The first structure constant is obtained by fitting the calibration results. For structural permeability, The inertial drag coefficient, A1 represents the second structure constant obtained by fitting the calibration results. It should be noted that A1 and A2 are used simultaneously for the dimensional coefficients of absorption.

[0134] Through the above offline calibration process, a basic calculation model is obtained that converts the permeability into the momentum drag coefficient of porous media.

[0135] After obtaining the basic calculation model, this embodiment specifically adopts the piecewise uniform assignment method to apply the source terms.

[0136] Step 505: Divide the open section longitudinally into several sub-regions.

[0137] Retrieve the coordinate axis parameters of the start and end points of the permeable section of the dike. Based on the set division spacing or quantity control command, along the longitudinal dimension parallel to the dike axis, truncate and divide the continuously distributed permeable section fluid computation domain into multiple independent sub-regions of equal physical length.

[0138] Step 506: Assign constant permeability and inertial drag coefficients to the computational grid within each sub-region.

[0139] For each independent sub-region formed by segmentation, calculate the corresponding local average permeability. Use the analytical expressions for permeability and inertial drag coefficient obtained from the calibration in the previous steps, and substitute them into the local average permeability for calculation. Traverse all hexahedral mesh elements within the boundary range of a specific sub-region, and batch-write the calculated permeability and inertial drag coefficient values ​​as constant attribute data into the momentum drag source term setting list of the corresponding mesh elements, overwriting them. This completes the discrete segmented drag loading of the porous media model, causing the water flow to experience stepped-varying drag as it passes through different sub-regions.

[0140] In one embodiment, the specific steps are as follows: calculate the local average permeability of each sub-region according to the longitudinal permeability gradient function; determine the permeability and inertial drag coefficient of each sub-region based on the pre-calibrated mapping relationship between permeability and drag coefficient; and assign the permeability and inertial drag coefficient to the computational grid inside each sub-region.

[0141] Example 6 differs from Example 5 in that it uses a piecewise uniform assignment method. In order to eliminate non-physical flow field jumps and oscillations at the interface of adjacent grids, this example provides a preferred method for assigning spatial coordinate continuous field functions.

[0142] Step 601: Obtain the pre-calibrated analytical expressions for permeability and inertial drag coefficient, substitute the longitudinal permeability gradient function in the key design parameter set into the analytical expressions for permeability and inertial drag coefficient, and transform the analytical expressions for permeability and inertial drag coefficient into continuous functions with respect to spatial coordinates. The analytical expressions for the basic fluid resistance model generated from the offline unit calibration test are retrieved from memory. These analytical expressions establish a purely algebraic mapping between permeability, inertial drag coefficient, and air permeability.

[0143] Furthermore, the previously established longitudinal permeability gradient function is extracted, and this function structure, which characterizes the distribution of permeability along the flow path, is substituted into the above analytical expression as an independent variable.

[0144] By performing algebraic composition operations, a direct relationship is established between the resistance parameters of porous media and the one-dimensional axis coordinates of the three-dimensional model space. Specifically, the continuous field distribution is established through the following formula: (16) (17) in, For spatial coordinates The local permeability function, The first structure constant is obtained by fitting the calibration results. Distance from the starting point of the transparent section The local permeability function at that location, For spatial coordinates The local inertial drag coefficient function, The second structure constant is obtained by fitting the calibration results. The coordinates are the spatial coordinates of the distance from the starting point of the open section along the longitudinal axis of the dike.

[0145] Step 602: Assign the continuous function about spatial coordinates to each computational grid centroid using a custom function to perform a spatially continuous and gradual assignment.

[0146] After establishing continuous analytical mapping relationships, the traditional region-level constant loading method is abandoned, and instead, the user-defined function interface in the computational fluid dynamics solver is invoked. During the initial compilation or data loading phase of the fluid numerical simulation, the user-defined function program traverses all fluid computational mesh elements within the geometric boundaries of the transparent segment.

[0147] For each target grid cell retrieved, the longitudinal coordinate value of the cell's spatial centroid is extracted, and this value is substituted into a continuous function of spatial coordinates for solution. The calculated independent permeability data and inertial drag coefficient data are used as local physical properties and directly written into the memory address of the momentum drag source term of the corresponding grid cell.

[0148] Through unit-level dynamic calculation and input strategies, a non-homogeneous field with smooth resistance transition is generated within the computational domain, enabling the numerical model to accurately fit the essential characteristics of the gradual transparency of the physical structure.

[0149] Step 603: In the region where the ratio of drag coefficients of adjacent computational grids exceeds a preset ratio threshold, a transition layer is introduced, and a hyperbolic tangent function is used to perform spatial smooth transition processing on the transition layer.

[0150] Despite the implementation of continuous coordinate mapping assignment, at specific spatial nodes with large curvature of the permeability distribution function, there may still be significant differences in drag gradients between adjacent grid cells, which may induce divergence or non-physical artifacts in the local momentum conservation calculation.

[0151] To mitigate this numerical defect, drag gradient monitoring is performed during mesh traversal loading. The inertial drag coefficients written to two adjacent fluid mesh cells are extracted and divided to obtain their ratio. A preset ratio threshold can be set to 1.5. When a calculated ratio higher than the preset threshold is detected, a high gradient abruptness risk is identified in the region, triggering an over-intervention command.

[0152] For example, a virtual numerical transition layer is inserted between the two grids that trigger the instruction, and the original local function assignment results are overwritten within this layer. A nonlinear hyperbolic tangent function is called as a smoothing filter to redistribute the drag weights. Specifically, spatial smoothing is performed using the following formula: (18) in, To calculate the local inertial drag coefficient of the target mesh after spatial smoothing, This represents the original inertial drag coefficient of the reference grid corresponding to the left boundary of the transition layer. This represents the original inertial drag coefficient of the reference grid corresponding to the right boundary of the transition layer. For hyperbolic tangent function operators, The centroid longitudinal coordinate of the current target mesh cell. The longitudinal coordinates of the center point of the transition layer set for the system. The smoothing scale parameter is used to control the span of the smoothing effect.

[0153] In the above processing logic, the specific value of the smoothing scale parameter is bound to the spatial discrete characteristics of the grid. The average grid side length of the corresponding region is read, and the smoothing scale parameter is configured as a set multiple of the grid side length.

[0154] By applying a transfer function with physical smoothing properties, the numerical step at the high-gradient interface is corrected, ensuring the stability of the iterative updates of the pressure and velocity fields during the solution of the Navier-Stokes equations.

[0155] In some embodiments, it can also be: the left boundary coordinates of the transition layer. and right boundary coordinates Take the longitudinal coordinates of the centroids of two adjacent grids whose trigger drag coefficient ratio exceeds the limit, and the coordinates of the center point of the transition layer. Pick and The arithmetic mean.

[0156] Example 7: This example is a detailed explanation of the steps involved in Example 1, such as constructing a proxy model and calling a multi-objective optimization algorithm to find the optimal solution.

[0157] Step 701: A non-dominated sorting genetic algorithm with constraint processing mechanism is adopted as the multi-objective optimization algorithm.

[0158] A multi-objective optimizer is initialized in memory, and a genetic algorithm module with non-dominated sorting properties is loaded. Basic operating parameters are set for the dike design problem, configuring constraints such as population size, maximum number of generations, crossover probability, and mutation probability. By establishing the above algorithmic infrastructure, a computational engine is provided for efficient global search in the multi-dimensional design parameter space.

[0159] Algorithm parameters such as population size, crossover probability, mutation probability, and maximum number of generations can be determined by those skilled in the art based on the dimensionality of decision variables and the surrogate model, through conventional parameter sensitivity tests.

[0160] Step 702: Based on the surrogate model, fitness evaluation is performed, and Pareto optimal solution set is obtained through population initialization, crossover and mutation iteration.

[0161] An initial parent population is randomly generated within the physical boundaries of the decision variables. When calculating the objective function value and constraint boundary state for each individual, time-consuming 3D fluid numerical calculations are not triggered; instead, a rapid data prediction is performed using the response surface surrogate model constructed in the previous embodiment. Subsequently, a offspring population is generated by performing simulated binary crossover and polynomial mutation operations on the population. The parent and offspring arrays are merged, and environmental selection based on survival of the fittest is repeatedly performed until the maximum number of generations is reached, outputting a set of non-dominated solution data converged on the multi-objective frontier.

[0162] Step 703, the constraint processing mechanism specifically uses the constraint dominance ranking principle to rate individuals in the population, including the following rating criteria.

[0163] An improvement to the conventional non-dominated sorting module is made by embedding a specialized constraint-dominated sorting operator in the environment selection stage to address the engineering characteristics of uncompromising physical constraints in the navigation channel flow regime. This operator handles the dominance relationship between feasible and infeasible solutions by distributing different decision logics.

[0164] Under the first working condition, it is determined that any feasible solution dominates any infeasible solution. Specifically, when one individual in one of the two populations participating in the comparison does not touch any hard constraint index red line in the surrogate model prediction, while the other individual has a record of exceeding the norm limit, regardless of the numerical superiority or inferiority of the two in other target indexes, the non-violation individual is directly given the retention priority to enter the next generation population.

[0165] Under the second operating condition, when comparing two infeasible solutions, the degree of constraint violation is calculated based on the extent of violation of the hard constraints of the navigation regulations, and the individual with the smaller constraint violation degree is determined to have dominance priority. Specifically, when both individuals violate physical constraints, their respective out-of-bounds limits are extracted and quantitatively compared. A mathematical penalty function is introduced to establish the calculation logic for the overall constraint violation degree, specifically calculated using the following formula: (19) in, For target individuals The total constraint violation rate is denoted by ∑, which represents the cumulative summation operation indicating the summation of all hard constraint indicators. For the maximum value function operator, To compare with the benchmark setting, For individuals In the The boundary function value is calculated on the hard constraints. A comparison operation is performed, and the individual with the smaller total constraint violation is determined as the winning individual, driving the algorithm to always point the search gradient in the infeasible region to the boundary of the feasible region.

[0166] In the calculation of constraint violation degree, the out-of-bounds extent of each constraint should be normalized to eliminate the dimensional differences between different physical quantities. The relevant normalization method can use existing technology, which will not be described in detail here.

[0167] Optionally, taking the maximum longitudinal velocity in the inlet region as an example, the corresponding constraint boundary calculation function is: When the maximum longitudinal velocity in the gate area predicted by the surrogate model exceeds the specification limit, the function value is positive, indicating that the constraint has been violated; where The maximum longitudinal flow velocity in the gate area corresponding to the design variable X under the surrogate model prediction; The longitudinal velocity limit for the entrance area as specified in the navigation regulations.

[0168] Under the third operating condition, when comparing two feasible solutions, they are ranked according to their non-dominated level. Specifically, when all compared individuals meet all hard flow limits, the non-dominated frontier levels are determined based on their numerical values ​​for target-type indicators such as crossflow and backflow. Among individuals in the same level, the crowding distance parameter is calculated for each individual, and individuals with larger crowding distance values ​​are retained to prevent the algorithm from getting stuck in local convergence.

[0169] Step 704: Obtain the target values ​​of each evaluation index in the Pareto optimal solution set.

[0170] After the optimization engine outputs the converged Pareto front set, the data parsing instruction is executed to traverse the non-dominated solution set matrix, extract the numerical vectors of target-type indicators such as the maximum lateral velocity of the gate region, the backflow intensity, and the backflow area corresponding to each feasible solution, and convert them from multi-dimensional array format into a list of target attributes for subsequent use by the distance measurement module.

[0171] Step 705: Perform normalization mapping processing on the target value.

[0172] Considering the order-of-magnitude differences in the dimensions and absolute values ​​of various evaluation indicators, to prevent the distance metric function from biasing towards the indicator dimensions with larger absolute values, extreme value parameters of each dimension component in the target attribute list are extracted. Each original target value is then transformed to a dimensionless standard interval of 0 to 1 through a linear mapping to eliminate scale interference caused by physical units.

[0173] Step 706: Calculate the normalized Euclidean distance from each solution point in the Pareto optimal solution set to the ideal reference point.

[0174] After achieving dimensionless scaling, a virtual ideal reference point is constructed in the multidimensional target space, where all evaluation indicators are minimized. This virtual ideal reference point corresponds to the origin, or ideal zero point, in the normalized space. For each non-dominated solution in the frontier set, its geometric spatial distance relative to this virtual ideal reference point in the multidimensional coordinate system is calculated using the following formula: (20) in, For the first The normalized Euclidean distances corresponding to the solution set points For the square root operator, This is an instruction to sum up all the target-type indicator dimensions currently being examined. For the first The solution is at the th solution. The original values ​​extracted from each target dimension. For the Pareto optimal solution set, the th The minimum value to search for each objective. For the Pareto optimal solution set, the th The maximum search value for each target.

[0175] Step 707: Select the design parameters corresponding to the solution set point with the smallest normalized Euclidean distance as the final design parameters.

[0176] The system receives the distance scalar calculation results from each solution point and calls the minimum value search and sorting function. The solution point ranked first with the smallest spatial distance scalar is identified as a compromise solution in the conflicting multi-objective index system. The data structure of this solution point is analyzed, and its bound transparent segment length parameter and polynomial coefficient matrix parameter are extracted as the finalized data for subsequent engineering model construction.

[0177] In a certain scenario, the spillway and lock of a water conservancy project are arranged side by side on the same side of the river. When the spillway is running at its design flow rate, the water flow from the spillway generates a strong lateral diffusion at the outlet section. The lateral velocity is the largest near the outlet of the spillway and gradually decreases along the longitudinal direction of the navigation channel.

[0178] Step 1: Construct a three-dimensional layered mathematical model of the river section including the layout of the hub and run the simulation to extract the velocity distribution outside the proposed dike. It was found that the maximum lateral velocity outside the dike exceeds the second preset velocity threshold.

[0179] Step Two: Based on this, it is determined that a long dike is needed, with its structural length extending to the tail of the pilot channel braking section. Simultaneously, considering that this location is downstream of the lock, an empirical negative correlation curve is used to determine the average permeability constraint range.

[0180] Step 3: Using a physical-driven reverse calculation method, based on the transverse velocity distribution along the dike, and through dynamic pressure equivalence, orifice outflow reverse calculation, and truncation processing, a gradually changing permeability distribution along the longitudinal direction of the dike is obtained. The permeability is lower at the outlet side of the spillway where the flow velocity is higher (more closed structure), and higher at the far end where the flow velocity is lower (more open structure). Coefficients are extracted through polynomial fitting, and iterative verification is used to ensure the average permeability falls within the constraint range.

[0181] Step 4: Load the gradual resistance source term in the 3D model using the spatial coordinate continuous field function method, apply tanh smooth transition in the high gradient region, and run batch simulations to extract evaluation indicators for each working condition.

[0182] Steps 5 and 6: Construct a response surface surrogate model and embed the constrained dominance NSGA-II optimizer. After iterative evolution, output the Pareto front and select the solution with the smallest normalized distance as the final design parameters.

[0183] Through the above process, the resulting gradually permeable dike can improve the lateral velocity and backflow index in the inlet area compared with the uniform permeable scheme, while meeting the longitudinal velocity specification limits of the inlet area and the approach channel under the working conditions.

[0184] The specific steps for NSGA-II are as follows: Set basic parameters, including population size N, maximum number of generations Gmax, crossover probability, mutation probability, crossover distribution index, and mutation distribution index. Set system constraints, including the permeable length of the dike and the initial permeability of the dike.

[0185] (1) Initialization: Randomly generate 80 individuals within the decision variable boundary as the initial parent population P0; (2) Evaluation: The objective function value of each individual was calculated using the response surface surrogate model. and ; (3) Non-dominated sorting: Performing a fast non-dominated sort on the population to divide individuals into different levels; (4) Crowding degree calculation: calculate the crowding degree distance for individuals in the same front layer to maintain the diversity of solutions; (5) Selection: A binary tournament selection operator is used to select parent individuals based on non-dominance level (prioritizing lower levels) and crowding level (prioritizing higher crowding levels); (6) Crossover and mutation: Perform simulated binary crossover and polynomial mutation on the selected parent individuals to generate the offspring population Q. t The size is N; (7) Merging and Elite Preservation: Merging the parent generation P t With offspring Q t Merge into R t , for R t Perform non-dominated sorting and crowding calculation, and select the top N individuals as the parent population P of the new generation. t+1 .

[0186] If the current generation is greater than the maximum number of generations or the Pareto front change is less than the threshold for 20 consecutive generations, the iteration is terminated; otherwise, return to step (3).

[0187] Experimental results show that, compared with the traditional dike design using a uniform permeability, the bottom-gradient permeability dike output by this invention reduces the maximum transverse velocity, backflow intensity, and backflow area in the inlet area to varying degrees. Furthermore, the maximum longitudinal velocity in both the inlet area and the approach channel meets the navigation code limits. This invention's method can achieve coordinated optimization of multi-objective flow conditions while meeting strict navigation constraints.

Claims

1. A design method for a bottom-gradient permeable lock dike based on multi-objective navigation flow condition optimization, characterized in that, include: A three-dimensional hierarchical mathematical model of the river section, including the lock, was constructed to obtain the flow velocity distribution outside the dike at the proposed location of the dike. The necessity of setting up the dike is determined based on the flow velocity distribution outside the dike, and the structural length and average permeability constraint range of the dike are initially determined. Under the constraints of the structural length and the average permeability range, the key design parameter set of the gradually permeable dike is determined, and the evaluation index set for evaluating the navigation flow conditions of the approach channel and the entrance area is defined. Based on the three-dimensional hierarchical mathematical model, numerical simulation is performed under the constraints of the key design parameter set to extract the evaluation index values ​​of navigation flow conditions under each working condition. Based on the evaluation index values ​​extracted from the numerical simulation, a proxy model is constructed to connect the set of key design parameters and the set of evaluation indexes. Based on the aforementioned proxy model, a multi-objective optimization algorithm is invoked to find the optimal solution, and the final design parameters of the bottom-gradient permeable dike that meet the hard constraints of navigation regulations and optimize the navigation flow conditions are output.

2. The method according to claim 1, characterized in that, The determination of the necessity of constructing the dike based on the flow velocity distribution outside the dike, and the preliminary determination of the structural length of the dike, specifically includes: Extract the maximum flow velocity outside the dike from the flow velocity distribution outside the dike; If the maximum flow velocity outside the dike is less than the first preset flow velocity threshold, it is determined that no dike needs to be set up, and the determination result that no dike needs to be set up is output. If the maximum flow velocity outside the dike is between the first preset flow velocity threshold and the second preset flow velocity threshold, it is determined that a short dike should be set, and the length of the structure is set to the length of the shield to the tail of the pilot channel mooring section. If the maximum flow velocity outside the dike is greater than the second preset flow velocity threshold, it is determined that a long dike should be set up, and the length of the structure is set to the length of the shield to the tail of the pilot channel braking section.

3. The method according to claim 1, characterized in that, The preliminary determination of the average permeability constraint range specifically includes: Identify the spatial location of the proposed diversion dam, whether it is upstream or downstream of the lock. Based on the spatial location and the flow velocity distribution outside the dike, the pre-stored empirical relationship between flow velocity and permeability is matched to determine the average permeability constraint range; Specifically, when the spatial location is upstream of the lock, the constraint boundary of the average permeability constraint range is set by the empirical relationship that the flow velocity outside the dike is positively correlated with the optimal permeability; when the spatial location is downstream of the lock, the constraint boundary of the average permeability constraint range is set by the empirical relationship that the flow velocity outside the dike is negatively correlated with the optimal permeability.

4. The method according to claim 1, characterized in that, The key design parameter set of the gradually permeable dike includes the length of the permeable section and the longitudinal permeability gradient function; The longitudinal permeability gradient function is determined by a mathematical curve preset method or a physical driving back-calculation method based on the fluid lateral momentum balance.

5. The method according to claim 4, characterized in that, When the mathematical curve preset method is used, the longitudinal permeability gradient function adopts a linear increasing function or an exponential increasing function, and after the permeability segment is longitudinally divided equally, a corresponding discrete permeability average value is assigned.

6. The method according to claim 4, characterized in that, When employing the physical-driven reverse propagation method based on fluid transverse momentum balance, the specific methods include: Extract the transverse velocity distribution along the longitudinal direction of the dike from the velocity distribution outside the dike; The dynamic pressure difference generated by the transverse velocity distribution is equivalent to the water level difference between the inside and outside of the dike. By using the preset maximum permissible lateral inflow velocity, combined with the water level difference inside and outside the dike and the outflow pattern of the orifice, the maximum permeability at each position along the longitudinal direction of the dike can be calculated. The maximum allowable permeability is truncated by using a preset upper limit for structural safety permeability and a preset lower limit for minimum achievable permeability to obtain a physically gradient permeability distribution. A polynomial fit is performed on the physical gradient permeability distribution, and the coefficients of the fitted polynomial are extracted as design variables to characterize the longitudinal permeability gradient function.

7. The method according to claim 1, characterized in that, The defined set of evaluation indicators is used to assess the navigation flow conditions of the pilotway and entrance area, specifically including: The evaluation index set is divided into hard constraint indexes and target indexes; The hard constraint indicators include the maximum longitudinal flow velocity in the entrance area and the maximum longitudinal flow velocity in the pilot channel, which have clearly defined limits in navigation regulations. The target indicators include the maximum transverse velocity in the inlet region, the backflow intensity, and the backflow area, which characterize the degree of flow field deterioration and need to be minimized.

8. The method according to claim 1, characterized in that, During the numerical simulation, the required simulation case sample points are generated using an adaptive experimental design method, specifically including: Obtain the dimension of the decision variables for the set of key design parameters; If the dimension of the decision variable is less than or equal to a preset dimension threshold, then the full factorial design method is used to generate the numerical simulation working condition sample points. If the dimension of the decision variable is greater than the preset dimension threshold, then the Latin hypercube sampling method is used to generate the numerical simulation working condition sample points.

9. The method according to claim 4, characterized in that, In the numerical simulation, a porous medium model is used to simulate the bottom gradually transparent structure, wherein the momentum drag source term is applied to the computational grid by a piecewise uniform assignment method or a spatial coordinate continuous field function assignment method. When the segmented uniform assignment method is used, the specific steps are as follows: The open section is divided into several sub-regions along the longitudinal direction; The local average permeability of each sub-region is calculated based on the longitudinal permeability gradient function. Based on the pre-calibrated mapping relationship between permeability and drag coefficient, the permeability and inertial drag coefficient of each sub-region are determined, and the permeability and inertial drag coefficient are assigned to the computational grid inside each sub-region.

10. The method according to claim 1, characterized in that, The step of invoking a multi-objective optimization algorithm based on the surrogate model to find the optimal solution specifically includes: A non-dominated sorting genetic algorithm with constraint handling mechanism is adopted as the multi-objective optimization algorithm. Fitness evaluation is performed based on the surrogate model, and Pareto optimal solution set is obtained through population initialization, crossover and mutation iteration. The constraint-based processing mechanism specifically employs the constraint-dominant ranking principle to rate individuals in the population, including the following rating criteria: It is determined that any feasible solution governs any infeasible solution; When comparing two infeasible solutions, the degree of constraint violation is calculated based on the degree of violation of the hard constraints of the navigation regulations, and the individual with the smaller degree of constraint violation is determined to have the priority of domination. When comparing two feasible solutions, they are ranked according to their non-dominated level.

Images