River channel three-dimensional water flow numerical simulation method

Abstract

The invention discloses a riverway three-dimensional water flow numerical simulation method, and relates to the technical field of intelligent water conservancy simulation, and the method comprises the steps: generating a longitudinal long and narrow hexahedral mesh of a main riverway in a C gradient direction through employing a mucous network growth algorithm, and outputting a hexahedral mesh with a shear layer identifier; turbulence viscosity coefficient distribution is dynamically initialized through a DDPG reinforcement learning model, and turbulence parameters are output; inputting the hexahedral mesh with the shear layer identifier and turbulence parameters into a coupling solver, performing parallel numerical solution, and outputting transient flow field data; a lambda 2 criterion is adopted to identify a vortex core area and calculate a vorticity modulus, a transverse circulation is defined, an on-way distribution cloud picture is generated, an immersive VR scene is constructed through a Unity 3D engine, and a three-dimensional visual water flow simulation scene with a physical feature label is output. Through fusion of lambda2 criterion vortex identification and a Unity 3D visualization technology, a complete solution with both calculation precision and engineering practicability is provided for river flow simulation.

Description

Technical Field

[0001] The invention relates to the technical field of intelligent water conservancy simulation, in particular to a three-dimensional water flow numerical simulation method for a river channel. Background Art

[0002] In recent years, with the development of computational fluid dynamics (CFD) and digital twin technology, three-dimensional numerical simulation of river flow has played an important role in flood control planning, ecological restoration and waterway management. Traditional methods mainly rely on structured grid or unstructured grid division technology, combined with Reynolds-averaged Navier-Stokes (RANS) equations or large eddy simulation (LES) for solution. However, high-precision surface reconstruction and feature preservation of complex river terrain still face challenges, especially in areas with sudden changes in curvature (such as steep slopes and bends), which are prone to grid distortion or physical field distortion.

[0003] In the prior art, river channel modeling based on traditional grid generation methods (such as octrees or T-splines) usually cannot take into account both the topological preservation of the longitudinal narrow structure of the main river channel and the high-resolution identification of the shear layer, which can easily lead to non-physical numerical dissipation in the grid transition zone. Especially in the area of sudden change in curvature, the existing surface reconstruction algorithms (such as radial basis functions or Poisson reconstruction) are not sensitive enough to the local geometric features of discrete point clouds, and may lose key hydraulic features (such as detached vortices or secondary flows). Summary of the invention

[0004] In view of the above existing problems, the present invention is proposed.

[0005] Therefore, the present invention provides a method for numerical simulation of three-dimensional water flow in a river channel to solve the problem of insufficient accuracy of discrete point cloud surface reconstruction in high curvature areas of the river channel in the prior art.

[0006] In order to solve the above technical problems, the present invention provides the following technical solutions:

[0007] In a first aspect, the present invention provides a method for numerical simulation of three-dimensional water flow in a river, which comprises obtaining river terrain point cloud data, performing preprocessing, reconstructing the surface of the discrete point cloud using a moving least squares method and calculating the curvature field, and generating a terrain surface with curvature labels and a high curvature feature point set;

[0008] The slime mold network growth algorithm is used to generate the longitudinal narrow hexahedral mesh of the main river channel along the C gradient direction, and the hexahedral mesh with shear layer mark is output;

[0009] Dynamically initialize the turbulent viscosity coefficient distribution through the DDPG reinforcement learning model and output the turbulence parameters;

[0010] The hexahedral mesh with shear layer identification and turbulence parameters are input into the coupling solver for parallel numerical solution and output of transient flow field data;

[0011] The λ2 criterion is used to identify the vortex core region and calculate the vorticity modulus, the transverse circulation is defined, the cloud map of the distribution along the course is generated, an immersive VR scene is constructed through the Unity3D engine, and a three-dimensional visual water flow simulation scene with physical feature tags is output.

[0012] As a preferred solution of the three-dimensional water flow numerical simulation method for the river channel described in the present invention, wherein: the river channel terrain point cloud data includes three-dimensional coordinates and echo intensity, reflectivity, and RGB color;

[0013] The preprocessing includes removing abnormal points for denoising, downsampling through voxel grid filtering, and performing high-precision alignment using the iterative closest point algorithm.

[0014] As a preferred solution of the three-dimensional water flow numerical simulation method for the river channel described in the present invention, wherein: the moving least squares method is used to reconstruct the surface of the discrete point cloud and calculate the curvature field, generating a terrain surface with curvature labels and a set of high-curvature feature points. The specific steps are as follows.

[0015] Define a local support domain centered on each point in the preprocessed river channel terrain point cloud data, fit a polynomial surface through weighted least squares, and project it onto the fitted surface to output smooth three-dimensional point cloud data and surface normal vectors;

[0016] Calculate the Gaussian curvature and the mean curvature through eigenvalue decomposition of the covariance matrix, generate a terrain surface with curvature labels, extract the curvature extreme points through the principal curvature threshold and non-maximum suppression, and output a set of high-curvature feature points.

[0017] As a preferred solution of the three-dimensional water flow numerical simulation method for the river channel described in the present invention, wherein: the slime network growth algorithm is used to generate a longitudinally narrow hexahedron grid of the main river channel along the C gradient direction. The specific steps are as follows.

[0018] Based on the terrain surface with curvature labels, use the central difference method to calculate the terrain curvature gradient of each grid node

[0019] Extract the bank line feature points from the set of high-curvature feature points, connect the adjacent feature points using Delaunay triangulation, and generate an initial slime tubular network along the connecting edges;

[0020] Normalize the terrain curvature gradient into a unit vector, map it to a dynamic nutrient concentration field through a Gaussian attenuation function, and calculate the growth direction deflection angle based on the terrain curvature gradient through the vector angle formula Calculate the growth direction deflection angle;

[0021] Optimize the spatial distribution and topological connectivity of the initial slime mold tubular network based on the dynamic nutrient concentration field and the growth direction deflection angle to obtain the slime mold tubular network;

[0022] Extract the longitudinal ridge line along the main path of the slime mold tubular network, and use the control points of the Bézier curve to align with the terrain curvature gradient along the normal direction, fit the cross-section, and adjust the aspect ratio of the hexahedral mesh elements according to the terrain curvature gradient to obtain the longitudinally narrow hexahedral mesh of the main river channel.

[0023] As a preferred scheme of the three-dimensional water flow numerical simulation method for the river channel described in the present invention, wherein: the specific steps of outputting the hexahedral mesh with shear layer identification are as follows.

[0024] Based on the longitudinally narrow hexahedral mesh of the main river channel, calculate the geometric shear factor of the unit, and mark the geometric shear factor of the unit exceeding the critical shear threshold as the candidate shear layer unit.

[0025] Trace the streamlines bidirectionally along the candidate shear layer unit at the local mesh size ratio step, and determine the shear layer where the curvature change rate exceeds the curvature change threshold.

[0026] Perform a closing operation on the shear layer to eliminate holes, and filter the discrete topological regions with an area order lower than the mesh discretization accuracy, and output the hexahedral mesh with shear layer identification.

[0027] As a preferred scheme of the three-dimensional water flow numerical simulation method for the river channel described in the present invention, wherein: the specific steps of dynamically initializing the turbulent viscosity coefficient distribution through the DDPG reinforcement learning model and outputting the turbulent parameters are as follows.

[0028] Use a graph convolutional network to perform feature encoding on the hexahedral mesh with shear layer identification, extract the grid topological connection relationship and physical field features, and obtain the water flow feature vector.

[0029] Normalize the historical water flow feature vector to zero mean and unit variance, convert it into a floating-point tensor, input it into the DDPG reinforcement learning model for training, and input the water flow feature vector converted into a floating-point tensor into the trained DDPG reinforcement learning model to obtain the turbulent parameters.

[0030] As a preferred scheme of the three-dimensional water flow numerical simulation method for the river channel described in the present invention, wherein: convert the hexahedral mesh with shear layer identification into the CGNS format, create an independent field for the turbulent parameters in the HDF5 storage structure, arrange them in the order of the unit ID of the hexahedral mesh with shear layer identification, input the processed hexahedral mesh with shear layer identification and the turbulent parameters into the coupled solver, perform parallel numerical solution through the velocity-pressure coupling algorithm, and output the transient three-dimensional velocity field and pressure field to obtain the transient flow field data.​​

[0031] As a preferred embodiment of the three-dimensional river flow numerical simulation method of the present invention, the following steps are included: output a three-dimensional visual water flow simulation scene with physical feature tags, and the specific steps are as follows.

[0032] Calculate the gradient tensor of the transient three-dimensional velocity field using the λ2 criterion, decompose it into a symmetric strain rate tensor and an antisymmetric rotation tensor, and perform eigenvalue solution and evaluation to identify the vortex core region.

[0033] Based on the transient three-dimensional velocity field, calculate the vorticity vector containing three spatial components, and obtain the scalar vorticity modulus characterizing the local vortex intensity by taking the square root of the sum of the squares of each component.

[0034] On the cross-section perpendicular to the mainstream direction of the vortex core region and the scalar vorticity modulus, select a circular path with a radius R and perform discrete integral calculation to output the lateral circulation distribution curve along the river centerline.

[0035] Interpolate the scalar vorticity modulus, lateral circulation distribution curve, and pressure field to a regular grid using radial basis functions, draw contour maps through the Slice filter, and convert the multi-time-step contour map sequence into a GIF animation to output the along-stream distribution contour map.

[0036] Convert the vortex core region, lateral circulation distribution curve, and along-stream distribution contour map into the FBX mesh format, encode the physical feature tags as vertex attributes, and perform dynamic visualization rendering through a customized Shader to output a three-dimensional visual water flow simulation scene with physical feature tags.

[0037] In a second aspect, the present invention provides a computer device, including a memory and a processor, where the memory stores a computer program, and: when the computer program is executed by the processor, any step of the three-dimensional river flow numerical simulation method described in the first aspect of the present invention is implemented.

[0038] In a third aspect, the present invention provides a computer-readable storage medium, on which a computer program is stored, and: when the computer program is executed by the processor, any step of the three-dimensional river flow numerical simulation method described in the first aspect of the present invention is implemented.

[0039] The beneficial effects of the present invention are as follows: By combining moving least squares surface reconstruction with the slime mold network growth algorithm, high-precision hexahedral mesh generation in high-curvature areas of river channels is achieved, solving the deficiencies of traditional methods in maintaining complex terrain features and identifying shear layers; at the same time, the DDPG reinforcement learning model is used to dynamically optimize turbulence parameters, significantly improving the accuracy of vortex structure capture and the computational convergence speed. Finally, through the integration of the λ2 criterion vortex identification and Unity3D visualization technology, an immersive VR scene with physical feature tags is constructed, providing a complete solution with both computational accuracy and engineering practicality for river channel flow simulation. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for description in the embodiments will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0041] Figure 1 It is a flowchart of the three-dimensional flow numerical simulation method for the river channel in Embodiment 1.

[0042] Figure 2 It is a flowchart of surface reconstruction and curvature field calculation in Embodiment 1.

[0043] Figure 3 It is a flowchart of generating hexahedral meshes by the slime mold network growth algorithm in Embodiment 1.

[0044] Figure 4 It is a flowchart of generating a three-dimensional visual flow simulation scene in Embodiment 1. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0045] In order to make the above objects, features, and advantages of the present invention more obvious and understandable, the following detailed description of the specific embodiments of the present invention will be made in conjunction with the accompanying drawings of the specification.

[0046] Many specific details are set forth in the following description in order to provide a thorough understanding of the present invention. However, the present invention may be implemented in other ways different from those described herein. Those skilled in the art can make similar generalizations without departing from the spirit of the present invention. Therefore, the present invention is not limited by the specific embodiments disclosed below.

[0047] Secondly, the so-called "one embodiment" or "embodiment" herein refers to specific features, structures, or characteristics that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a separate or alternative embodiment that excludes other embodiments.

[0048] Example 1, referring to Figures 1 to 4 , this example provides a three-dimensional water flow numerical simulation method for a river channel, including the following steps:

[0049] S1. The river channel terrain point cloud data includes three-dimensional coordinates and echo intensity, reflectivity, and RGB color.

[0050] S1.1. The preprocessing includes removing abnormal points for denoising, downsampling through voxel grid filtering, and high-precision alignment using the iterative closest point algorithm.

[0051] It should be noted that:

[0052] Removing abnormal points for denoising is specifically: calculating the k-nearest neighbors of each point (e.g., k = 50), statistically calculating the average distance μ and standard deviation σ of the points within the neighborhood; setting a threshold (e.g., μ ± 3σ), traversing all points and removing the points whose average distance from the neighborhood exceeds the threshold, and retaining the effective terrain points;

[0053] Downsampling through voxel grid filtering is specifically: dividing the point cloud space into cubic grids of a fixed size (e.g., side length 0.1m), calculating the centroid of the internal points of each non-empty voxel or randomly selecting a representative point, and replacing all the points within the original voxel with the representative point to achieve data compression;

[0054] High-precision alignment using the iterative closest point algorithm is specifically: performing an initial rough registration on two point clouds (e.g., manual or PCA alignment), iteratively performing the closest point search (source point → target point) and solving the optimal rigid body transformation (rotation matrix R + translation vector t) through SVD decomposition until the transformation converges (ΔR < 0.001) or reaches the maximum number of iterations (e.g., 100 times).

[0055] S2. Using the moving least squares method to perform surface reconstruction on the discrete point cloud and calculate the curvature field, generating a terrain surface with curvature labels and a set of high-curvature feature points.

[0056] S2.1. Defining a local support domain centered on each point in the preprocessed river channel terrain point cloud data, fitting a polynomial surface through weighted least squares, and projecting it onto the fitted surface, outputting smooth three-dimensional point cloud data and surface normal vectors.

[0057] It should be noted that, based on the KD tree data structure for spatial indexing, then using radius search (setting the search radius r) or K-nearest neighbor search (setting the number of neighborhood points k) to obtain the neighborhood point set, assigning weights to each neighborhood point, and obtaining a weighted local support domain;

[0058] Under the local support domain, a spatial grid hash index is used to quickly partition the neighborhood points in space. The elevation values of the cube vertices are calculated by trilinear interpolation, and then the isosurface patches are extracted based on the Marching Cubes rule. Finally, the local terrain surface is approximately expressed by triangular meshes. The weights are calculated according to the distances from the neighborhood points to the target point, and the closer the distance, the greater the weight. The coefficients of the surface equation are solved by weighted least squares. For any point (x, y) in the target area, the coordinates of the point are substituted into the surface equation to calculate the corresponding z value, so as to determine the exact position of the point coordinates in the target area on the fitted surface. Finally, all the calculated surface points are connected to construct a complete optimal fitted surface. The target point is vertically projected onto the optimal fitted surface to obtain the smoothed new coordinates;

[0059] Calculate the normal vector direction of the surface at the smoothed new coordinates, and normalize it to obtain the unit normal vector. Repeat this process for all points, and finally output the smoothed three-dimensional point cloud data and the surface normal vector.

[0060] S2.2. Calculate the Gaussian curvature and the mean curvature by eigenvalue decomposition of the covariance matrix, generate the terrain surface with curvature labels, extract the curvature extreme points through the principal curvature threshold and non-maximum suppression, and output the high-curvature feature point set.

[0061] It should be noted that based on the smoothed three-dimensional point cloud data, for each point, search for the neighborhood point set within k nearest neighbors or within a radius r, calculate the mean of the three-dimensional coordinates of these neighborhood points as the center point, use the smoothed surface normal vector to verify the neighborhood consistency (remove the abnormal points with the normal vector angle > 30°), construct a 3×3 covariance matrix, and perform eigenvalue decomposition on the 3×3 covariance matrix to obtain three eigenvalues where is the variance in the first principal direction (the direction with the largest curvature change), is the variance in the secondary principal direction (within the tangent plane orthogonal to я1), is the variance in the normal direction (since the terrain is locally approximated as a plane, tends to 0); Calculate the Gaussian curvature (the product of the eigenvalues) and the mean curvature (half of the sum of the eigenvalues) according to the eigenvalues, and label the curvature value for each point to generate the terrain surface with labels;

[0062] Set the high quantile value (such as the top 10% quantile) of the curvature statistical distribution in the smoothed three-dimensional point cloud data as the principal curvature threshold. Based on the smoothed three-dimensional point cloud data, screen all the points in the smoothed three-dimensional point cloud data through the principal curvature threshold, select the points with significantly high curvature as the candidate set, and then use the non-maximum suppression algorithm (compare the curvature values of the points in each candidate set with the points in the neighborhood, and only retain the local maximum points), and finally output the set of high-curvature feature points that meet the conditions.

[0063] S3. Generate a longitudinally narrow hexahedron grid for the main river channel along the C gradient direction using the slime mold network growth algorithm.

[0064] S3.1. Based on the terrain surface with curvature labels, use the central difference method to calculate the terrain curvature gradient of each grid node.

[0065] It should be noted that the expression for calculating the terrain curvature gradient of each grid node is:

[0066]

[0067] where is the terrain curvature gradient, is the first-order partial derivative of the terrain surface C with curvature labels along the x-axis direction, is the first-order partial derivative of the terrain surface C with curvature labels along the y-axis direction, is the first-order partial derivative of the terrain surface C with curvature labels along the z-axis direction.

[0068] S3.2. Extract the riverbank line feature points from the high-curvature feature point set, use Delaunay triangulation to connect adjacent feature points, and generate an initial slime mold tubular network along the connecting edges.

[0069] It should be noted that density clustering (such as DBSCAN) is performed on the high-curvature feature point set to remove discrete noise points and retain the continuously distributed riverbank line candidate points; constrained Delaunay triangulation is performed based on the two-dimensional projection coordinates (ignoring elevation) of the riverbank line candidate points to generate a triangular network that only connects adjacent feature points, extract all the edges in the triangular network as the centerline of the initial slime mold tubular network, and calculate the initial radius of the pipeline according to the elevation difference between the two endpoints of each edge (the greater the elevation difference, the smaller the initial radius), and perform fusion processing on the crossed or too-close pipelines to generate an initial tubular network with a spatial topological structure;

[0070] The elevation difference is the difference in the original Z coordinates of the riverbank line candidate points in the high-curvature feature point set (i.e., the height difference between two points in the vertical direction).

[0071] S3.3. Normalize the terrain curvature gradient into a unit vector, map it to a dynamic nutrient concentration field through the Gaussian attenuation function, and calculate the growth direction deflection angle based on the terrain curvature gradient using the vector angle formula.

[0072] It should be noted that:

[0073] The mapping to the dynamic nutrient concentration field through the Gaussian attenuation function has the expression:

[0074]

[0075] Among them, F(x, y, z) represents the dynamic nutrient concentration field function to quantify the growth excitation intensity of the slime mold tubular network in the three-dimensional space (x, y, z). represents the modulus of the terrain curvature gradient, and d is the shortest distance from the node of the slime mold tubular network to the riverbank line;

[0076] Through the vector angle formula, the growth direction deflection angle is calculated based on the terrain curvature gradient, and the expression is:

[0077]

[0078] Among them, θ is the growth direction deflection angle, is the vertical direction vector of the current slime mold tubular branch (orthogonal to the main trunk flow direction).

[0079] S3.4. Optimize the spatial distribution and topological connectivity of the initial slime mold tubular network based on the dynamic nutrient concentration field and the growth direction deflection angle to obtain the slime mold tubular network.

[0080] It should be noted that the nutrient concentration F and the terrain curvature gradient at the node are calculated to determine the optimal growth direction, and it preferentially grows towards the neighborhood with high F and aligned with the terrain curvature gradient (θ is small); the inefficient branches with F values lower than the threshold (0.3F) are removed through the competition mechanism, and the adjacent pipe segments with a spacing < 0.1L and a direction angle < 15° are fused; the short closed loops with too low F values are detected and broken, and the iteration is optimized until the network is stable (reaching the maximum number of iterations, such as 100 times), and finally the slime mold tubular network with reasonable spatial distribution and optimized topological connection is output.

[0081] S3.5. Extract the longitudinal ridge line along the main trunk path of the slime mold tubular network, and use the control points of the Bézier curve to align with the terrain curvature gradient along the normal direction to fit the cross-section, and adjust the aspect ratio of the hexahedral mesh elements according to the terrain curvature gradient to obtain the longitudinally narrow hexahedral mesh of the main channel.

[0082] It should be noted that the topological analysis of the slime mold tubular network is carried out, and the main trunk path with the highest nutrient concentration is extracted as the initial ridge line, and the continuity of the ridge line is optimized by using the node sampling algorithm based on the curvature change (such as sampling a point every 0.5 times the characteristic length L, and densifying the sampling at the curvature mutation point);

[0083] A local coordinate system is established at each sampling point (the X-axis is along the tangent direction of the ridge line, the Y-axis is along the direction, and the Z-axis is determined according to the right-hand rule), and the control points of the Bézier curve are arranged in the XY plane, where the offset of the control point in the Y direction (k is the scaling factor, usually taken as 0.1 - 0.3), to make the cross-sectional shape adapt to the change of terrain curvature;

[0084] Use a 3rd-order Bézier curve to generate a closed cross-sectional profile in the XY plane, and the number of control points is dynamically adjusted according to the local value (4 as the base, and 1 additional control point is added for every 0.1 times the maximum curvature gradient); construct an initial hexahedron mesh based on the cross-sectional profile, and adjust the aspect ratio r = 1 + 4α2 of the element according to the standardized curvature gradient intensity at the center point of the element (ensuring r ∈ [1, 5]);

[0085] Keep the transition smooth through a physics-based mesh deformation algorithm (such as As-Rigid-As-Possible), and finally output a longitudinally narrow hexahedron mesh that highly coincides with the river channel terrain, with the long axis strictly aligned with the ridge line direction and the short axis precisely matching the curvature gradient field distribution.

[0086] S4. Output a hexahedron mesh with shear layer identification.

[0087] S4.1. Based on the longitudinally narrow hexahedron mesh of the main river channel, calculate the geometric shear factor of each element, and mark the elements with geometric shear factors exceeding the critical shear threshold as candidate shear layer elements.

[0088] It should be noted that:

[0089] Calculate the center point curvature gradient of each element based on the vertex coordinates of the hexahedron mesh elements and solve the rate of change of the normal vector of each side surface of the element through the finite difference method. Combine the angle between the flow direction and the curvature gradient within the element to calculate the geometric shear factor of the element, and the expression is:

[0090]

[0091] where S is the geometric shear factor of the element, K is the angle between the flow direction and the curvature gradient direction, is the spatial change amount in the normal direction;

[0092] Compare the calculated geometric shear factor S of the element with the preset critical shear threshold E (usually taken as the 90% quantile of the distribution of the geometric shear factor S of the element), and screen out all elements with S ≥ E critical shear threshold and mark them as candidate shear layer elements.

[0093] S4.2. Trace the streamlines bidirectionally along the candidate shear layer elements with a local mesh size ratio step, and determine the curvature change rate exceeding the curvature change threshold as the shear layer.

[0094] It should be noted that starting from the center of the candidate shear layer unit, streamlines are traced in both forward and reverse flow directions with an adaptive step size (initially 0.2 times the unit size). The RK4 method is used to integrate the velocity field and the step size is adjusted in real time (0.1 - 0.5 times the unit size). During the tracing process, the change rate of the curvature gradient at the current point is calculated at each step. When κ exceeds the dynamic threshold (such as 2 times the regional median) for three consecutive steps, the streamline segment is marked as a candidate shear layer region. After all streamlines are traced, the discrete candidate regions are topologically connected (merging adjacent regions with a spacing < 1.5 times the unit size) and morphological processing (closing operation to fill holes) is performed. Finally, a binary shear layer identification field and a complete shear layer streamline network are output, and the geometric characteristic parameters (length, average κ value, spatial range) of each shear layer are recorded.

[0095] S4.3. Perform a closing operation on the shear layer to eliminate holes, and filter out discrete topological regions with an area magnitude lower than the grid discretization accuracy, and output a hexahedral grid with shear layer identification.

[0096] It should be noted that for the morphological closing operation on the binary shear layer identification field, a 3×3×3 cubic structuring element is used to perform the dilation and then erosion operations in three-dimensional space to eliminate small holes and smooth the shear layer boundary.

[0097] Based on connected component analysis, all shear layer regions are extracted, the number of voxels (corresponding to the physical area) in each region is calculated, and discrete regions with the number of voxels less than the grid discretization accuracy threshold (such as 5 times the minimum unit volume) are filtered out.

[0098] The processed shear layer identification (1 represents an effective shear layer, 0 represents a non-shear layer) is written into the hexahedral grid data structure in the form of cell attributes, and the curvature gradient at the shear layer boundary is recorded at the grid nodes. value to generate a hexahedral grid with shear layer identification having a complete physical identification.

[0099] S5. Dynamically initialize the turbulent viscosity coefficient distribution through the DDPG reinforcement learning model and output turbulent parameters.

[0100] S5.1. Use a graph convolutional network to perform feature encoding on the hexahedral grid with shear layer identification, extract the grid topological connection relationship and physical field features, and obtain the water flow feature vector.

[0101] It should be noted that the hexahedral grid with shear layer identification is converted into graph structure data, where the grid cells are used as nodes, and the node features include physical attributes (such as shear layer identification, curvature gradient cell aspect ratio) and geometric attributes (such as cell volume, center coordinates), and the connection relationship between adjacent cells is used as edges.

[0102] Feature encoding is performed using a multi-layer graph convolutional network (GCN). In each layer, the node features are updated by aggregating the information of neighboring nodes (such as mean pooling or attention weighting), and at the same time, the topological correlation features between units are extracted through EdgeConv. The feature vectors of all nodes are aggregated through a graph readout function (such as global mean pooling or attention aggregation) to generate a global feature vector with a fixed dimension. Mean pooling is performed on the node features in the shearing layer identification area to extract local high-dimensional embedding features. The global feature vector and the local high-dimensional embedding features are concatenated along the feature dimension to form a combined vector, which is then reduced in dimension by a fully connected layer and output as a water flow feature vector of a unified length.

[0103] S5.2. Normalize the historical water flow feature vector to zero mean and unit variance, convert it into a floating-point tensor, and input it into the DDPG reinforcement learning model for training. Input the water flow feature vector converted into a floating-point tensor into the trained DDPG reinforcement learning model to obtain the turbulence parameters.

[0104] It should be noted that the historical water flow feature vectors are preprocessed by standardization. The mean and standard deviation of each dimension of the historical water flow feature vectors are calculated, and z-score normalization is performed to obtain the historical water flow feature vectors with zero mean and unit variance. The normalized historical water flow feature vectors are converted into 32-bit floating-point tensors, and an experience replay buffer (usually with a capacity of 1e5 - 1e6 entries) is constructed after adding the batch dimension (batch_size×feature_dim).

[0105] In the DDPG training stage, a small batch of data (batch_size = 64 - 256) is sampled from the buffer in each round and input into an Actor-Critic framework composed of a 4-layer fully connected network (with hidden layer dimensions of 256 - 512 and LeakyReLU activation). The Actor network outputs a turbulence parameterization strategy (constrained within a range by Tanh activation), and the Critic network evaluates the state-action Q value. The policy gradient is calculated through the temporal difference error (TD-error), and the network parameters are updated using an Adam optimizer (Actor learning rate of 1e-4, Critic learning rate of 1e-3), and the target network is synchronized periodically (soft update coefficient τ = 0.01). The training termination condition is that the Critic loss function converges (Δloss < 1e-5) or the maximum number of rounds is reached (usually 1e5 steps), and finally the trained DDPG model is output.

[0106] Input the water flow feature vector converted into a floating-point tensor into the trained DDPG reinforcement learning model. Through the forward propagation of the Actor network (3 fully connected layers, 512-dimensional hidden layer, LeakyReLU activation), high-order features are extracted layer by layer. Finally, the output layer activated by Tanh generates a normalized turbulent parameter recommendation value. Subsequently, physical constraint processing is performed on the normalized turbulent parameter recommendation value (for example, ensuring that the turbulent viscosity coefficient is non-negative through the Softplus function and multiplying by a preset upper limit value to restore the actual dimension). At the same time, the Critic network calculates the matching degree Q value between the turbulent parameter recommendation value and the water flow feature vector, filters out high-confidence parameters through a threshold (such as Q value > 0.8), replaces low-confidence parameters with the historical mean of the sliding window, and combines high-confidence parameters and the substituted low-confidence parameters to form a turbulent parameter tensor corresponding to the input batch.

[0107] S6. Convert the hexahedral mesh with shear layer identification to the CGNS format, create an independent field for the turbulent parameters in the HDF5 storage structure, arrange them in the order of the cell IDs of the hexahedral mesh with shear layer identification, input the processed hexahedral mesh with shear layer identification and the turbulent parameters into the coupled solver, perform parallel numerical solution through the velocity-pressure coupling algorithm, and output the transient three-dimensional velocity field and pressure field to obtain transient flow field data.

[0108] Specifically as follows:

[0109] Extract the cell topology (node connection relationship), physical properties (such as cell ID, shear layer identification), and turbulent parameters (such as turbulent viscosity, turbulent kinetic energy) of the hexahedral mesh with shear layer identification, organize the data according to the CGNS standard structure, write the cell data to the cell topology nodes, write the node coordinates to the grid coordinate nodes, and store the initial flow field variables (such as initial velocity, pressure) at the flow field solution nodes. At the same time, create a turbulent parameter group in the HDF5 storage structure, store the turbulent parameters in the order of cell ID, ensuring one-to-one correspondence with the grid cells. Finally, output an HDF5 file that conforms to the CGNS standard, containing complete grid, physical property, and turbulent parameter data;

[0110] Read the CGNS file, parse the HDF5 structure, load the grid topology, node coordinates, initial flow field, and turbulent parameters into memory, match the turbulent parameters according to the cell ID, and map them to the corresponding computational cells. Initialize the solver computational domain, allocate MPI processes for parallel domain decomposition (such as METIS partitioning), ensure that each process only processes local grid data, and establish an inter-process communication mechanism to synchronize boundary data.

[0111] Finally, complete the solver initialization;

[0112] Based on the pressure field at the current time step, the predicted velocity field is obtained by solving the momentum equation (momentum prediction step). The pressure field is corrected by the pressure Poisson equation to satisfy the mass conservation condition (pressure correction step). The velocity field is corrected using the updated pressure gradient to ensure the continuity equation holds (velocity correction step);

[0113] Combined with the turbulence parameters (such as turbulent viscosity and turbulent kinetic energy) output by the DDPG model, the turbulence transport equation (such as k-ε or k-ω model) is solved, the turbulence parameters are updated and fed back into the momentum equation to achieve the dynamic coupling of turbulence characteristics;

[0114] Throughout the process, each MPI process is responsible for the calculation of the local grid area. Through non-blocking communication (MPI_Isend / MPI_Irecv), the velocity, pressure, and turbulence parameter data of the boundary cells are exchanged in real time to ensure the consistency of the global flow field. The iteration of each time step continues until the residuals of the momentum equation and the pressure correction equation drop to the set residual convergence threshold (such as 1e-5) or reach the maximum number of iterations, and finally the updated velocity field, pressure field, and turbulence parameter field that satisfy the physical conservation law are output;

[0115] At the end of each time step, the local velocity field and pressure field calculated by each process are collected, merged into the global flow field data, and written into the transient HDF5 file in CGNS format to obtain the transient flow field data.

[0116] S7. Identify the vortex core region using the λ2 criterion, calculate the vorticity modulus, define the transverse circulation, generate the cloud map of the distribution along the path, construct an immersive VR scene through the Unity3D engine, and output the three-dimensional visual water flow simulation scene with physical feature labels.

[0117] S7.1. Calculate the gradient tensor of the transient three-dimensional velocity field using the λ2 criterion, decompose it into the symmetric strain rate tensor and the anti-symmetric rotation tensor, and perform eigenvalue solving and evaluation to identify the vortex core region.

[0118] It should be noted that the transient velocity field is numerically discretized in three-dimensional space, and the finite difference or finite element method is used to calculate the gradient tensor at each grid point It is decomposed into the symmetric part (strain rate tensor ) and the anti-symmetric part (rotation tensor ) through eigenvalue decomposition; then the vorticity discrimination matrix L = G2 + Ω2 is constructed, the eigenvalues λ1, λ2, λ3 (sorted in ascending order) are solved, the potential vortex core cells are marked based on the λ2 criterion (the region where λ2 < 0 is the vortex core), the morphological closing operation (3×3×3 cube kernel) is performed on the potential vortex core cells to eliminate noise, and the spatially continuous vortex core region is extracted through connected component analysis.

[0119] S7.2. Calculate the vorticity vector containing three spatial components based on the three-dimensional transient velocity field. By taking the square root of the sum of the squares of each component, the scalar vorticity magnitude characterizing the local vortex intensity is obtained.

[0120] It should be noted that based on the three-dimensional transient velocity field, the central difference method is used to calculate the partial derivatives of the velocity field in three directions to obtain the rate of change of the velocity components with respect to each coordinate axis. Then, according to the definition of vorticity, the three components of the vorticity vector are calculated using these partial derivatives, which respectively describe the characteristics of local fluid motion rotating around the x, y, and z axes. Square each component of the vorticity vector and sum them to characterize the overall intensity of the local vortex. Finally, take the square root of the obtained sum to get the scalar vorticity magnitude.

[0121] S7.3. Select a circular path with a radius R on the cross-section perpendicular to the mainstream direction of the vortex core region and the scalar vorticity magnitude, and perform discrete integral calculations to output the lateral circulation distribution curve along the centerline of the river channel.

[0122] It should be noted that a series of cross-sections perpendicular to the mainstream direction are generated along the centerline of the river channel. On each cross-section, a circular integration path is defined with the center of the vortex core region as the origin and a radius R (usually taking 2 - 3 times the characteristic vortex core size). The circular integration path is discretized into N equally angularly spaced points (N ≥ 36 to ensure accuracy). Interpolate the scalar vorticity magnitude and velocity components at the discrete points to calculate the tangential velocity. At each discrete integration point, multiply the tangential velocity by the scalar vorticity magnitude at the corresponding position as a weighting factor to calculate the local circulation contribution. Subsequently, numerically integrate along the closed path (such as the trapezoidal rule) to obtain the equivalent circulation of the cross-section, and record the corresponding centerline mileage coordinates to output the lateral circulation distribution curve along the centerline of the river channel.

[0123] S7.4. Interpolate the scalar vorticity magnitude, the lateral circulation distribution curve, and the pressure field to a regular grid using radial basis functions, draw contour maps through the Slice filter, and convert the multi-time-step contour map sequence into a GIF animation to output the contour maps of the longitudinal distribution.

[0124] It should be noted that based on the radial basis function (RBF) interpolation method, map the scalar vorticity magnitude, the lateral circulation distribution curve, and the pressure field on the irregular grid to a regular grid (such as 0.1 times the characteristic length size), and use a compactly supported Gaussian kernel function for scatter data fitting to ensure smooth transition of the physical field gradient;

[0125] Use a Slice filter to generate a series of equally spaced cross-sections (spacing ≤ 0.2 times the river width) along the centerline of the river on a regular grid, extract the interpolated data of each cross-section and render it as a 2D cloud map (e.g., the vorticity magnitude is rendered as a heat map, and the pressure field is rendered with contour lines superimposed). Sort the cloud map sequence for each time step according to the timestamp, and use the FFmpeg tool to synthesize a GIF animation (frame rate 10fps, color mapping optimized to the Viridis scheme) to dynamically display the vortex evolution process, and output a cloud map of the along-channel distribution containing spatio-temporal evolution characteristics.

[0126] S7.5. Convert the vortex core region, the transverse circulation distribution curve, and the along-channel distribution cloud map into the FBX mesh format, encode the physical feature tags as vertex attributes, and perform dynamic visualization rendering through a customized Shader to output a 3D visual water flow simulation scene with physical feature tags.

[0127] It should be noted that the voxel data of the vortex core region, the transverse circulation distribution curve, and the along-channel distribution cloud map are converted into triangular meshes. Among them, the surface of the vortex core is reconstructed using the Marching Cubes algorithm, and the circulation curve generates a renderable geometry through tubular meshing (radius adaptively to the circulation value).

[0128] When exporting to FBX, encode the physical features (such as vorticity magnitude, circulation value, pressure) as vertex attributes (the RGB channels correspond to different physical quantities, and the Alpha channel stores the normalized intensity). At the same time, retain the time step information as vertex animation frame data. Load the FBX asset in the Unity / Unreal engine and achieve dynamic rendering through a customized Shader: the vertex shader interpolates the animation frames according to the time parameter, and the fragment shader mixes the Phong lighting and the isosurface shading based on the physical tags (e.g., a blue-white gradient for the vortex core region, a red-yellow gradient for the high-circulation region), and add a streamline particle system (driven by the velocity field) to enhance the performance of vortex motion, and output a 3D visual water flow simulation scene with physical feature tags.

[0129] This embodiment also provides a computer device, which is applicable to the case of the three-dimensional water flow numerical simulation method for rivers, including: a memory and a processor; the memory is used to store computer-executable instructions, and the processor is used to execute the computer-executable instructions to implement the three-dimensional water flow numerical simulation method for rivers proposed in the above embodiment.

[0130] The computer device may be a terminal, which includes a processor, a memory, a communication interface, a display screen, and an input device connected via a system bus. Among them, the processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The communication interface of the computer device is used to communicate with external terminals in a wired or wireless manner. The wireless manner can be achieved through WIFI, a carrier network, NFC (Near Field Communication), or other technologies. The display screen of the computer device can be a liquid crystal display screen or an electronic ink display screen. The input device of the computer device can be a touch layer covering the display screen, or buttons, a trackball, or a touchpad provided on the housing of the computer device, or an external keyboard, touchpad, or mouse, etc.

[0131] This embodiment also provides a storage medium, on which a computer program is stored. When the program is executed by a processor, it implements the method for realizing three-dimensional river flow numerical simulation proposed in the above embodiment; the storage medium can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as static random access memory (Static Random Access Memory, abbreviated as SRAM), electrically erasable programmable read-only memory (Electrically Erasable Programmable Read-Only Memory, abbreviated as EEPROM), erasable programmable read-only memory (Erasable Programmable Read Only Memory, abbreviated as EPROM), programmable read-only memory (Programmable Red-Only Memory, abbreviated as PROM), read-only memory (Read-Only Memory, abbreviated as ROM), magnetic memory, flash memory, a magnetic disk, or an optical disc.

[0132] In summary, the present invention combines moving least squares surface reconstruction with the slime mold network growth algorithm to achieve high-precision hexahedral mesh generation in high-curvature areas of the river channel, and solves the deficiencies of traditional methods in maintaining complex terrain features and identifying shear layers; at the same time, the DDPG reinforcement learning model is used to dynamically optimize turbulence parameters, significantly improving the accuracy of vortex structure capture and the computational convergence speed. Finally, through the integration of the λ2 criterion vortex identification and Unity3D visualization technology, an immersive VR scene with physical feature tags is constructed, providing a complete solution with both computational accuracy and engineering practicality for river channel flow simulation.

[0133] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit them. Although the present invention has been described in detail with reference to the preferred embodiments, those of ordinary skill in the art should understand that the technical solutions of the present invention can be modified or equivalently replaced without departing from the spirit and scope of the technical solutions of the present invention, and they should all be covered within the scope of the claims of the present invention.

Claims

1. A three-dimensional numerical simulation method for river flow, characterized in that: including, obtaining the point cloud data of the river channel terrain and performing preprocessing, using the moving least squares method to reconstruct the surface of the discrete point cloud and calculate the curvature field, generating a terrain surface with curvature labels and a high-curvature feature point set; using the slime mold network growth algorithm to generate a long and narrow hexahedral grid for the main river channel longitudinally along the C gradient direction, and outputting a hexahedral grid with shear layer identification; dynamically initializing the distribution of turbulent viscosity coefficients through a DDPG reinforcement learning model and outputting turbulent parameters; inputting the hexahedral grid with shear layer identification and turbulent parameters into a coupled solver for parallel numerical solution and outputting transient flow field data; using the λ2 criterion to identify the vortex core region and calculate the vorticity modulus, defining the transverse circulation, generating a cloud map of the distribution along the path, constructing an immersive VR scene through the Unity3D engine, and outputting a three-dimensional visual water flow simulation scene with physical feature labels.

2. The three-dimensional water flow numerical simulation method for river channels according to claim 1, wherein: The point cloud data of the river channel terrain includes three-dimensional coordinates, echo intensity, reflectivity, and RGB color; The preprocessing includes removing abnormal points for denoising, downsampling through voxel grid filtering, and performing high-precision alignment using the iterative closest point algorithm.

3. The three-dimensional water flow numerical simulation method for river channels according to claim 2, characterized in that: The steps of using the moving least squares method to reconstruct the surface of the discrete point cloud and calculate the curvature field, generating a terrain surface with curvature labels and a high-curvature feature point set are as follows: Defining a local support domain centered on each point in the preprocessed point cloud data of the river channel terrain, fitting a polynomial surface through weighted least squares, and projecting it onto the fitted surface, outputting smooth three-dimensional point cloud data and surface normal vectors; Calculating the Gaussian curvature and mean curvature through eigenvalue decomposition of the covariance matrix, generating a terrain surface with curvature labels, and extracting curvature extreme points through the principal curvature threshold and non-maximum suppression, outputting a high-curvature feature point set.

4. The three-dimensional water flow numerical simulation method for river channels according to claim 3, characterized in that: The steps of using the slime mold network growth algorithm to generate a long and narrow hexahedral grid for the main river channel longitudinally along the C gradient direction are as follows: Based on the terrain surface with curvature tags, the terrain curvature gradient of each grid node is calculated using the central difference method Extracting the bank line feature points from the high-curvature feature point set, using Delaunay triangulation to connect adjacent feature points, and generating an initial slime mold tubular network along the connecting edges; Normalize the topographic curvature gradient into a unit vector, map it into a dynamic nutrient concentration field through a Gaussian attenuation function, and calculate the growth direction deflection angle based on the topographic curvature gradient through the vector included angle formula Calculate the growth direction deflection angle; Optimizing the spatial distribution and topological connectivity of the initial slime mold tubular network based on the dynamic nutrient concentration field and the growth direction deflection angle to obtain the slime mold tubular network; Extract the longitudinal ridge line along the main path of the slime mold tubular network, and use the control points of the Bézier curve to align with the terrain curvature gradient in the normal direction and fit the cross-section. According to the terrain curvature gradient Adjust the aspect ratio of the hexahedron mesh elements to obtain a longitudinally narrow hexahedron mesh for the main river channel according to the terrain curvature gradient and adjust the aspect ratio of the hexahedron mesh elements to obtain a longitudinally narrow hexahedron mesh for the main river channel 5. The three-dimensional water flow numerical simulation method for river channels according to claim 4, wherein: The steps of outputting a hexahedral grid with shear layer identification are as follows: Based on the long and narrow hexahedral grid for the main river channel longitudinally, calculating the element geometric shear factor, and marking the element geometric shear factors exceeding the critical shear threshold as candidate shear layer elements; Tracking the streamlines bidirectionally along the candidate shear layer elements at the local grid size ratio step, and determining the curvature change rate exceeding the curvature change threshold as the shear layer; Performing a closing operation on the shear layer to eliminate holes and filtering discrete topological regions with an area order lower than the grid discretization accuracy, and outputting a hexahedral grid with shear layer identification.

6. The three-dimensional water flow numerical simulation method for river channels according to claim 1, characterized in that: The steps of dynamically initializing the distribution of turbulent viscosity coefficients through a DDPG reinforcement learning model and outputting turbulent parameters are as follows: Using a graph convolutional network to perform feature encoding on the hexahedral grid with shear layer identification, extracting the grid topological connection relationship and physical field features, and obtaining a water flow feature vector; Normalize the historical water flow feature vector to zero mean and unit variance, convert it into a floating-point tensor, and input it into the DDPG reinforcement learning model for training. Input the water flow feature vector converted into a floating-point tensor into the trained DDPG reinforcement learning model to obtain the turbulence parameters.

7. The three-dimensional water flow numerical simulation method for river channels according to claim 1, characterized in that: Convert the hexahedral mesh with shear layer identification into the CGNS format, and create an independent field for the turbulence parameters in the HDF5 storage structure, arranged in the order of the cell IDs of the hexahedral mesh with shear layer identification. Input the processed hexahedral mesh with shear layer identification and the turbulence parameters into the coupled solver, and perform parallel numerical solution through the velocity-pressure coupling algorithm to output the transient three-dimensional velocity field and pressure field, obtaining the transient flow field data.

8. The three-dimensional water flow numerical simulation method for river channels according to claim 1, wherein: Output the three-dimensional visual water flow simulation scene with physical feature labels. The specific steps are as follows. Use the λ2 criterion to calculate the gradient tensor of the transient three-dimensional velocity field, decompose it into a symmetric strain rate tensor and an antisymmetric rotation tensor, and perform eigenvalue solution and evaluation to identify the vortex core region. Based on the transient three-dimensional velocity field, calculate the vorticity vector containing three spatial components, and obtain the scalar vorticity modulus representing the local vortex intensity by taking the square root of the sum of the squares of each component. On the cross-section perpendicular to the mainstream direction of the vortex core region and the scalar vorticity modulus, select a circular path with a radius R and perform discrete integral calculation to output the lateral circulation distribution curve along the river centerline. Use the radial basis function to interpolate the scalar vorticity modulus, the lateral circulation distribution curve, and the pressure field to a regular grid, draw a contour map through the Slice filter, and convert the multi-time-step contour map sequence into a GIF animation to output the along-channel distribution contour map. Convert the vortex core region, the lateral circulation distribution curve, and the along-channel distribution contour map into the FBX mesh format, encode the physical feature labels as vertex attributes, and perform dynamic visualization rendering through a customized Shader to output the three-dimensional visual water flow simulation scene with physical feature labels.

9. A computer device, comprising a memory and a processor, the memory storing a computer program, characterized in that: When the processor executes the computer program, it implements the steps of the three-dimensional water flow numerical simulation method of the river channel according to any one of claims 1 to 8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by the processor, it implements the steps of the three-dimensional water flow numerical simulation method of the river channel according to any one of claims 1 to 8.

Images