Method for simulating local scour of structure under unsteady flow

Abstract

The application discloses a kind of local scour simulation methods of structure under the action of unsteady flow, it is related to fluid and silt transport simulation technical field, comprising the following steps: step 1: generating calculation area grid;Step 2: build limited area projection grid;Step 3: the calculation area grid is divided into several sub-regions;Step 4: start parallel solution, calculate interpolation weight;Step 5: carry out the solution of flow field;Step 6: calculate the displacement amount of the center point at bed surface grid due to silt transport at each time step;Step 7: calculate the displacement amount of the grid point at bottom boundary, and update grid point coordinates;Step 8: if there is the slope of grid face is greater than the angle of repose, use silt sliding algorithm to correct bed surface grid point coordinates;Step 9: update calculation grid.The method of the application can be parallel computing, and compared with prior art, improve simulation precision and solution stability.

Description

Technical Field

[0001] This invention relates to the field of fluid and sediment transport simulation technology, and in particular to a method for simulating local scour of structures under unsteady flow. Background Technology

[0002] Localized scour of structures such as pile foundations and pipelines in marine engineering often has a significant impact on the safety of these structures. In recent years, with the development of computer technology, many scholars have conducted numerical simulation studies on localized scour of marine structures. Among them, the single-phase flow model combining the RANS method and the ALE dynamic mesh method has been widely used in the simulation research of the above-mentioned problems due to its high computational efficiency and ability to accurately capture the morphological changes of the seabed.

[0003] Different scholars have constructed local scour models of structures using the finite element method, finite difference method, and finite volume method (FVM). These different models share similar solution algorithm flows, as detailed below: Step 1: Treat the water flow as an incompressible Newtonian fluid and use the Navier-Stokes equations (NS equations) to solve for the flow field characteristics (velocity, pressure, etc.): ; ; In the formula, the subscript f represents the water flow phase, and the superscript... In RANS simulations, t represents the time-averaged component, and t represents the time term. Represents the fluid velocity vector. Represents fluid pressure. and These represent the density and kinematic viscosity of the water flow, respectively. To represent the turbulent stress tensor, a suitable turbulence model needs to be introduced to address the turbulence terms in the equations. To close the loop, the following method is usually selected. The SST model performs turbulent closure because it can better predict boundary layer separation flow under the influence of the reverse pressure gradient.

[0004] Step 2: After the flow field solution is completed, calculate the bed shear stress based on the near-bottom velocity gradient of the sediment bed. Based on the bed shear stress and sand particle size Calculate the Shields number of the bed surface : ; In the formula, This represents the density of sediment.

[0005] Step 3: Based on the Shields number of the bed surface Shields number for sediment initiation The relationship between the two determines whether the mud and sand on the bed surface have been disturbed. The bed sediment transport rate can then be calculated using an appropriate bed sediment transport rate formula. (such as Van Rijn's formula, Meyer-Peter's formula, etc.); if Then the bed surface sand transport rate will be... It is 0.

[0006] Step 4: Solve for suspended sediment transport within the computational region using the convection-diffusion equation: ; In the formula, Represents suspended sediment concentration (or volume fraction). Represents turbulent kinematic viscosity. Represents the turbulent Schmidt number. This represents the settling velocity of sediment.

[0007] Step 5: Based on the bedload and suspended sediment transport rates calculated in Steps 3 and 4, solve for the elevation change at the bottom of the bed using the continuity equation (Exner equation) for bed sediment transport: ; In the formula, h represents the change in bed height over a certain period of time. Represents the porosity of the bed. and The sediment erosion rate and sediment deposition rate, representing the bed surface, can be obtained from the change in suspended sediment concentration at the bottom boundary: ; .

[0008] Step 6: Based on the calculated changes in the bed surface topography, update the positions of grid points on the bed surface and within the calculation area using the dynamic mesh algorithm.

[0009] Step 7: Repeat steps 1-6 within each time step until the flushing reaches a stable state or the simulation time is greater than the given end time.

[0010] In recent years, although different scholars have developed single-phase flow sediment transport model solvers for simulating local scour processes of structures in OpenFOAM based on the above algorithms, the existing technology still has the following technical problems: (1) Through code testing and analysis, it was found that the currently open-source scour models / solvers usually cannot perform parallel computing, so the computational efficiency for simulating long-term scour processes is low and cannot meet the requirements for solution efficiency in practical engineering applications. (2) At present, there is still room for further optimization of the scour model based on OpenFOAM. The setting of boundary conditions or solution algorithms in the model can easily lead to insufficient simulation accuracy, which in turn leads to poor solution stability.

[0011] In view of this, this invention is hereby proposed. Summary of the Invention

[0012] The purpose of this invention is to overcome the shortcomings of existing technologies by proposing a method for simulating local scour of structures under unsteady flow conditions. This method allows for parallel computation and improves simulation accuracy and solution stability compared to existing technologies.

[0013] To achieve the above objectives, the present invention adopts the following technical solution: A method for simulating local scour of a structure under unsteady flow includes the following steps: Step 1: Generate the computational domain mesh and set the boundary conditions; Step 2: Construct a finite-area projected mesh with the same planar topology as the bottom boundary mesh; Step 3: Divide the computational domain into several sub-regions based on the number of parallel solution cores; Step 4: Start parallel solving and calculate the interpolation weights of different mesh surfaces based on the projected mesh generated in Step 2; Step 5: Solve the flow field. Use rough wall boundary conditions to calculate and assign values ​​to the flow velocity and turbulence variables at the bottom boundary. Use unsteady inflow boundary conditions to calculate and assign values ​​to the flow velocity and turbulence variables at the inflow boundary at different times. Step 6: Solve for bedload transport and suspended sediment transport, and calculate the displacement at the center point of the bed grid due to sediment transport at each time step; Step 7: Using parallel solving, calculate the displacement at the grid points on the bottom boundary through interpolation, and update the grid point coordinates; Step 8: Determine the slope of each grid surface. If there is a grid surface with a slope greater than the angle of repose, use the sediment slippage algorithm to correct the coordinates of the grid points on the bed surface until the slope of all grid surfaces is less than or equal to the angle of repose. Step 9: Update the computational grid and use the new grid point coordinates as the bottom boundary positions for the flow field solution at the next time step.

[0014] Furthermore, in step 4, based on the projected mesh generated in step 2, the interpolation weights of different mesh surfaces are calculated, including the following steps: Step 4.1: Calculate the sum of the reciprocals of the distances from each grid point to the center of all its adjacent faces. .

[0015] Step 4.2: Calculate the interpolation weights on different grid surfaces using the following formula. : ; In the formula, Represents the center of the i-th mesh face With grid points The reciprocal of the distance between them.

[0016] Furthermore, in step 4.1, for each grid point inside the processor, the reciprocal of the distance from the grid point to the adjacent face center point is calculated and summed to obtain the result. For points belonging to two processors, calculate the sum of the reciprocals of the distances between the grid points and their respective processors. Then, use the Pstream class to implement data exchange between the coupling edges of adjacent processors, and sum the corresponding reciprocals of the distances to obtain the result. For points belonging to three or more processors, the globalData class in the finite area method is used to implement data transmission and coupling between multiple processors, thereby obtaining all grid points. And sum to obtain the final .

[0017] Furthermore, in step 5, the flow velocity and turbulence variables at the bottom boundary are calculated and assigned values ​​using rough wall boundary conditions, including the following steps: Step 5.1: Adjust the distance from the bed surface This serves as the boundary of the bottom rough surface; Step 5.2: Calculate the flow velocity parameters at the bottom boundary using the following formula: ; in, ; ; ; In the formula, The tangential velocity representing the water flow at the bottom boundary. Represents the Karman constant. Represents the roughness height of the bed surface. This represents the height at the center point c of the first-layer grid, obtained from the governing equations of the hydrodynamic model in the single-phase flow model. The normal vector representing the mesh surface. This represents the flow velocity at point c. Represents the velocity of the geometric center point of the mesh surface. This represents the velocity of a grid point in the direction perpendicular to the surface normal caused by its movement. Step 5.3: Calculate the parameters of the turbulence variables at the bottom boundary using the following formula: ; ; In the formula, The value representing the turbulent kinetic energy at the bottom boundary. Represents a constant. This represents the dissipation rate of turbulent kinetic energy at the bottom boundary.

[0018] Furthermore, in step 5, the flow velocity and turbulence variables at the inflow boundary at different times are calculated and assigned using non-steady inflow boundary conditions, including the following steps: Step 5.4: Calculate the bed surface friction velocity at different times using the following formula: ; In the formula, The velocity of the water flow at the upstream reference position of the pipeline at different times; Step 5.5: Calculate the parameters of the inflow boundary velocity at different times using the following formula: ; In the formula, The horizontal velocity representing the water flow at the inflow boundary at different times. The maximum value in the horizontal velocity profile representing the flow motion at the inflow boundary at different times; Step 5.5: Calculate the parameters of the turbulence variables at the inflow boundary at different times using the following formula: ; ; in, ; In the formula, The distribution of turbulent kinetic energy at the inflow boundary. Represents the inflow boundary layer thickness. The distribution representing the dissipation rate of turbulent kinetic energy at the inflow boundary.

[0019] Furthermore, in step 6, when solving for suspended sediment transport, the suspended sediment concentration at the bottom boundary is calculated using the following formula: ; In the formula, Represents the Shields number of the bed surface. This represents the critical Shield number for sediment initiation.

[0020] Furthermore, in step 7, the displacement at the grid points on the bottom boundary is calculated using the following formula: ; In the formula, This represents the displacement at the center point of the bed grid due to sediment transport.

[0021] Furthermore, in step 8, the coordinates of the grid points on the bed surface are corrected, and the change in bed surface elevation is calculated using the following formula: ; In the formula, Represents the bed surface elevation. The sediment slip diffusion coefficient; in, ; ; In the formula, Let C represent the angle of repose, and C represent a constant that can range from 1.5 to 2.4. This represents the maximum bed sediment transport rate at that moment.

[0022] Compared with the prior art, the beneficial effects of this invention are as follows: Parallel computation is possible through the projection in step 2 and the interpolation in step 4. The simulation accuracy is improved through the calculation of the flow velocity and turbulence variables at the bottom boundary and the inflow boundary in step 5, and the suspended sediment concentration at the bottom boundary in step 6. In summary, the stability of the solution is ensured through parallel computation and the improvement of simulation accuracy. Attached Figure Description

[0023] Figure 1 This is a flowchart of the method for simulating local scour of a structure in Example 1; Figure 2 This is a schematic diagram of the projection bed surface in Embodiment 1; Figure 3 This is a schematic diagram of the interpolation process from the face center to the grid point between adjacent processors in the parallel computing of Example 1; Figure 4 This is a comparison diagram of the original boundary condition locations and the boundary condition locations in Example 1; Figure 5 This is a schematic diagram of the constant flow scour calculation region in Example 2; Figure 6Diagram showing the changes in the near-bottom region grid topology during scouring using existing technologies; Figure 7 Example 1: Changes in the mesh topology of the near-bottom region during the scouring process; Figure 8 This is a comparison diagram of the scour terrain profile results simulated by the method in Example 1 and those in references 2, 3, and 4; Figure 9 This is a schematic diagram of the constant flow scour calculation region in Example 3; Figure 10 This is a comparison diagram of the simulated inflow velocity at position 8D upstream of the pipeline center in Example 1 and the measured velocity in Reference 5. Figure 11 This is a comparison chart of the simulated value of the dimensionless scour depth G / D below the pipeline center as a function of scour time in Example 1 and the experimental value in Reference 5. Detailed Implementation

[0024] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0025] Example 1: A method for simulating local scour of structures under unsteady flow, such as Figure 1 As shown, it includes the following steps: Step 1: Generate the computational domain mesh and set the boundary conditions.

[0026] In this embodiment, simulation is performed using OpenFOAM, which allows meshing tools such as pointwise or blockMesh to generate the computational domain mesh. Boundary conditions include inflow, outflow, top boundary, bottom boundary, rough wall, and smooth wall.

[0027] Step 2: Construct a finite area projected mesh with the same planar topology as the bottom boundary mesh.

[0028] In the existing technology, the model is solved using a bottom boundary mesh that is completely consistent with the shape of the bed surface. However, due to the curvature of the bed surface, difference errors are inevitably introduced.

[0029] Therefore, in this embodiment, as Figure 2 As shown, a projected mesh with the same planar topology as the bottom boundary mesh is constructed (in specific implementation, it can be constructed based on the finite area algorithm of OpenFOAM), and a parallel interpolation algorithm is constructed on the projected mesh to solve the sediment slip model.

[0030] Step 3: Divide the computational domain grid into several sub-regions based on the number of parallel solution cores.

[0031] In this embodiment, the number of parallel solution cores can be specified according to computer performance.

[0032] In this embodiment, the decomposePar command can be used to divide the computational domain mesh into several sub-regions based on the number of parallel solution cores.

[0033] Step 4: Start parallel solving and calculate the interpolation weights of different mesh surfaces based on the projected mesh generated in Step 2.

[0034] In this embodiment, the mpirun command can be used to start parallel solving.

[0035] In the numerical simulation of the scouring process of a moving bed, the computational workload is much greater than that of the numerical simulation of the hydrodynamic process of a fixed bed flow field due to the involvement of sediment transport and grid motion simulation. Parallel solution can improve simulation efficiency.

[0036] OpenFOAM's parallel architecture uses the Pstream class as the MPI communication protocol to achieve cross-processor information transfer, and utilizes processor boundary conditions to assign values ​​to the boundaries of adjacent computational domains between different processors. OpenFOAM's parallel architecture can be applied to the parallel solution of hydrodynamics, suspended sediment transport, bedload transport, and bed change in scour models. However, existing parallel algorithms in OpenFOAM cannot be directly applied to dynamic mesh updates and sediment slippage calculations for the following reasons: In solving dynamic mesh models or sediment slip models, it is necessary to obtain the displacement or displacement velocity of each mesh point on the bottom boundary. The displacement at the center of each mesh surface on the bottom boundary can be calculated using the Exner equation. Therefore, an interpolation algorithm is needed to interpolate the displacement at the face center to the grid points. OpenFOAM can perform interpolation from the face center to the grid points on a boundary surface using the primitive class (primitivePatchInterpolation class). While the above interpolation algorithm can perform interpolation operations between the face center and grid points on a boundary surface, the primitive class does not consider the special handling of grid points at the boundary in a parallel computing domain. Therefore, existing parallel algorithms in OpenFOAM cannot be directly applied.

[0037] In this embodiment, step 4, based on the projected mesh generated in step 2, calculates the interpolation weights for different mesh surfaces, including the following steps: Step 4.1: Calculate the sum of the reciprocals of the distances from each grid point to the center of all its adjacent faces. .

[0038] Step 4.2: Calculate the interpolation weights on different grid surfaces using the following formula. : ; In the formula, Represents the center of the i-th mesh face With grid points The reciprocal of the distance between them.

[0039] like Figure 3 As shown in the figure, the red dots represent the face centers, and the black, blue, and green dots are grid points at different locations where interpolation needs to be calculated. The three colors of grid points represent three different interpolation cases. Figure 3 Different Processor numbers represent adjacent computational domains located in different processors, and interpolation calculations are independent within each computational domain. Existing technologies use interpolation algorithms in the `primitive` class to calculate the displacement of grid points (black dots) within each processor. However, for grid points at boundaries, their displacement is only affected by the displacement of the center of the adjacent face within the processor they belong to (as shown in the figure, the blue dot belongs to two processors, and the green dot belongs to four processors). This results in the same grid point having different displacements in different processors, further causing computational divergence. Therefore, in parallel computing, it is not possible to achieve face-to-grid interpolation solely using the `primitive` class in OpenFOAM.

[0040] Based on this, in this embodiment, in step 4.1, for a grid point inside the processor, the reciprocal of the distance from the grid point to the adjacent face center point is calculated and summed to obtain the result. For points belonging to two processors, calculate the sum of the reciprocals of the distances between the grid points and their respective processors. Then, use the Pstream class to implement data exchange between the coupling edges of adjacent processors, and sum the corresponding reciprocals of the distances to obtain the result. For points belonging to three or more processors, the globalData class in the finite area method is used to implement data transmission and coupling between multiple processors, thereby obtaining all grid points. And sum to obtain the final .

[0041] Step 5: Solve the flow field. Use rough wall boundary conditions to calculate and assign values ​​to the flow velocity and turbulence variables at the bottom boundary. Use unsteady inflow boundary conditions to calculate and assign values ​​to the flow velocity and turbulence variables at the inflow boundary at different times.

[0042] In this embodiment, step 5 involves calculating and assigning values ​​to the flow velocity and turbulence variables at the bottom boundary using rough wall boundary conditions, including the following steps: Step 5.1: Adjust the distance from the bed surface This serves as the boundary of the bottom rough surface.

[0043] like Figure 4 As shown, for reference Figure 4 In the original OpenFOAM definition, 'a' means the bottom face of the computational mesh is at the same position as the wall, i.e., a no-slip boundary condition is applied at point 'b', the center of the bottom boundary face of the mesh. Although OpenFOAM also includes a built-in boundary condition for rough walls, it requires the center point 'c' of the first layer of mesh points to be located within the logarithmic region. Therefore, the first layer of mesh near the bottom boundary needs a larger height, with its center height... At least equal to the roughness height of the bed surface However, in simulations of processes such as pipeline scouring, to ensure the topology of the mesh is not destroyed during scouring, a small initial scouring pit is typically set between the lower surface of the pipeline and the bed boundary. If the mesh size near the bottom boundary is large, this results in fewer meshes between the pipeline and the bed surface within the initial scouring pit, leading to uneven mesh distribution and poor resolution. Furthermore, during scouring simulations, mesh stretching within the scouring pit further reduces the accuracy of the simulation of the interstitial flow region below the pipeline. On the other hand, using the original boundary conditions in OpenFOAM presents challenges in assigning values ​​to the sediment concentration at the boundary. The near-bottom concentration of suspended sediment can be obtained using empirical formulas, but the calculated concentration is at a reference height above the sediment bed surface. place ( Figure 4 The concentration value at point b in the dashed line position of 'a'. This reference height is... Figure 4 The position of the bottom boundary of the computational domain shown in 'a' (point c) is not consistent, and there is a certain gap between them, so it is impossible to accurately assign a value to the bottom boundary.

[0044] Therefore, in this embodiment, the bottom boundary of the computational grid is no longer located at the surface of the sediment bed, but is instead located at a distance from the bed surface. Using this point as the boundary, the bottom boundary position coincides with the reference height of suspended sediment concentration, and the governing equations for the flow field and suspended sediment transport can be solved within the same grid. In practical implementation, twice the sediment particle size can be used as... .

[0045] Step 5.2: Calculate the flow velocity parameters at the bottom boundary using the following formula: ; in, ; ; ; In the formula, The tangential velocity representing the water flow at the bottom boundary. This represents the Karman constant, which can be 0.4. Represents the roughness height of the bed surface. This represents the height at point c, the center of the first-layer grid, obtained from the governing equations of the hydrodynamic model in the single-phase flow model. In practical implementation, this can be achieved by... , represent, This represents the flow velocity at point c. The velocity representing the movement of the geometric center point of the mesh surface can be obtained from the difference in coordinates of the geometric center point before and after the mesh movement. This represents the velocity perpendicular to the surface normal caused by the movement of a grid point. Its value is equal to the velocity flux caused by the movement of the grid point divided by the area of ​​the grid surface.

[0046] Step 5.3: Calculate the parameters of the turbulence variables at the bottom boundary using the following formula: ; ; In the formula, The value representing the turbulent kinetic energy at the bottom boundary. This represents a constant and can be 0.09. This represents the dissipation rate of turbulent kinetic energy at the bottom boundary.

[0047] In this embodiment, step 5 involves calculating and assigning values ​​to the flow velocity and turbulence variables at the inflow boundary at different times using non-constant inflow boundary conditions, including the following steps: Step 5.4: Calculate the bed surface friction velocity at different times using the following formula: ; In the formula, This represents the flow velocity at a reference location upstream of the pipeline at different times, such as the measured velocity at a certain location upstream of the pipeline or the cross-sectional average velocity. The changes over time can be determined by specifying a corresponding change function based on the simulated unsteady flow process, or by reading the velocity change sequence from a time series file.

[0048] Step 5.5: Calculate the parameters of the inflow boundary velocity at different times using the following formula: ; In the formula, The horizontal velocity representing the water flow at the inflow boundary at different times. The maximum value in the horizontal velocity profile representing the flow motion at the inflow boundary at different times.

[0049] Step 5.6: Calculate the parameters of the turbulence variables at the inflow boundary at different times using the following formula: ; ; in, ; In the formula, The distribution of turbulent kinetic energy at the inflow boundary. Represents the inflow boundary layer thickness. The distribution representing the dissipation rate of turbulent kinetic energy at the inflow boundary.

[0050] In this embodiment, to achieve numerical simulation of continuous velocity changes during unsteady flow scouring, a boundary condition setting method based on logarithmic velocity distribution is adopted at the inflow boundary to specify the vertical distribution of velocity and related turbulence variables, exhibiting good adaptability and physical consistency. Furthermore, considering the evolution characteristics of velocity over time under unsteady flow, the velocity profile at the inflow boundary is... It is constructed by coupling the classical logarithmic distribution form with a time function.

[0051] Step 6: Solve for bedload transport and suspended sediment transport, and calculate the displacement at the center point of the bed grid due to sediment transport at each time step.

[0052] In this embodiment, in step 6, the solutions for bedload transport and suspended sediment transport, as well as the calculation of the displacement at the center point of the bed grid due to sediment transport, can be obtained using the existing technologies mentioned in the background art, which will not be elaborated here.

[0053] In this embodiment, in step 6, when solving for suspended sediment transport, the suspended sediment concentration at the bottom boundary is calculated using the following formula: ; In the formula, Represents the Shields number of the bed surface. This represents the critical Shield number for sediment initiation.

[0054] Step 7: Using parallel solving, calculate the displacement at the grid points on the bottom boundary through interpolation, and update the grid point coordinates.

[0055] In this embodiment, step 7 is calculated based on the calculation results of steps 2 and 6.

[0056] Specifically, in step 7, the displacement at the grid points on the bottom boundary is calculated using the following formula: ; In the formula, This represents the displacement at the center point of the bed grid due to sediment transport.

[0057] Step 8: Determine the slope of each grid surface. If there is a grid surface with a slope greater than the angle of repose, use the sediment slippage algorithm to correct the coordinates of the grid points on the bed surface until the slope of all grid surfaces is less than or equal to the angle of repose.

[0058] During scouring, if the local slope of the bed is too steep, exceeding a certain critical slope, the sediment on the bed will collapse under gravity, causing a decrease in the local slope. This critical slope is generally the underwater angle of repose of the sediment. However, if no special restrictions are imposed on the bed surface slope during numerical calculations, the local slope of the bed surface may gradually increase. When the bed surface slope angle is greater than... At that time, the bed will face the boundary Shields number. It will approach 0, while the bed surface sediment transport rate The slope will approach infinity, which will cause the calculation to diverge. Therefore, numerical simulations need to use a sediment slippage model to increase the constraint on the local bed slope.

[0059] A more common model is the geometry-based slip model proposed by Niemann et al. (2010). After solving for the coordinates of the bed surface grid points at a new moment, it is necessary to check the slope of each grid face on the bottom boundary. If a face with a slope greater than the sediment repose angle is found, the vertical position of the grid points on that face needs to be slip corrected. However, it should be noted that after slip correction of a face, it may cause the slope angle of its adjacent faces to exceed the sediment repose angle. Therefore, after correcting the grid points on a face, it is necessary to continue traversing all faces on the bottom boundary to check if there are still faces with a slope angle greater than the repose angle, and correct them again. The traversal stops when the slope angle of all faces on the bottom boundary meets the slope stability requirements. In summary, the existing method has significant limitations in sediment slip calculation: its solution process must concentrate all bottom boundary elements on the same processor in order to construct the relative slope relationship between elements. This centralized data approach is only suitable for two-dimensional problems with a limited number of bottom boundary meshes; once the mesh size is large or three-dimensional scouring needs to be dealt with, this traversal method of correcting only one face at a time will be very time-consuming.

[0060] Based on this, in step 8 of this embodiment, the coordinates of the grid points on the bed surface are corrected, and the change in bed surface elevation is calculated using the following formula: ; In the formula, Represents the bed surface elevation. The sediment slip diffusion coefficient; in, ; .

[0061] In the formula, Let C represent the angle of repose, and C represent a constant that can range from 1.5 to 2.4. This represents the maximum bed sediment transport rate at that moment.

[0062] Step 9: Update the computational grid and use the new grid point coordinates as the bottom boundary positions for the flow field solution at the next time step.

[0063] Step 10: Continue the cycle according to the time step until the set solution time is reached or the flushing reaches a stable state.

[0064] This embodiment of a method for simulating local scour of structures under unsteady flow conditions allows for parallel computation through projection in step 2 and interpolation calculation in step 4. Improved calculations of flow velocity and turbulence variables at the bottom and inflow boundaries in step 5, and improved suspended sediment concentration at the bottom boundary in step 6, enhance simulation accuracy. In summary, the stability of the solution is ensured through parallel computation and improved simulation accuracy.

[0065] Example 2: This embodiment is a verification embodiment to verify the effectiveness of the method in Embodiment 1.

[0066] Numerical simulation verification of local scouring in a single pipeline under steady flow was conducted based on experimental data from Reference 1 (Mao Y. The interaction between a pipeline and an erodible bed [D]. Copenhagen: Technical University of Denmark, 1986.). Reference 1 conducted experiments on a moving bed and clear water scouring of a single pipeline under steady flow in a water tank measuring 23m long, 2m wide, and 0.5m deep. The water depth H in the experimental section was 0.35m, the pipeline diameter D was 0.1m, and the median particle size of the sediment on the bed surface was [not specified]. It is 0.36mm. Figure 5 The diagram illustrates the computational domain layout in this simulation. The computational domain is 60D in length and 10D in height, with the pipeline center 30D from the inflow boundary. An initial scour pit with a depth of 0.1D is placed below the pipeline to prevent computational divergence caused by changes in the mesh topology during the scour simulation.

[0067] Figure 6 The changes in the mesh topology of the near-bottom region during scouring using existing techniques (referring to existing OpenFOAM algorithms) are presented. Figure 7The changes in the mesh topology of the near-bottom region during the scouring process of the method in Example 1 are shown. It can be found that the bottom mesh distribution of the method in Example 1 is more uniform, and there are no problems such as negative volume and erroneous distribution due to the near-bottom mesh, which improves the stability of the model solution.

[0068] Figure 8 The simulation results of Example 1 at t=30min are compared with those of Reference 2 (Liang D, Cheng L, Li F. Numerical modeling of flow and scour below a pipeline in currents: Part II. Scour simulation. Coastal Engineering, 2005, 52(1): 43-62.), Reference 3 (Zhao M, Cheng L. Numerical modeling of local scour below apiggyback pipeline incurrents. Journal of Hydraulic Engineering, 2008, 134(10): 1452-1463.), and Reference 4 (Larsen BE, Furman DR, Sumer BM. Simulation of wave-plus-current scourbeneath submarine pipelines. Journal of Waterway, Port, Coastal, and Ocean Engineering, 2016, 142(5): 04016003.). It can be observed that the simulation results of Example 1 are closer to the experimental results of Reference 1, especially the simulation accuracy of the downstream sand dunes has been significantly improved.

[0069] Example 3: This embodiment is a verification embodiment to verify the effectiveness of the method in Embodiment 1.

[0070] Based on the local scour experiment data of a single pipeline under unsteady flow conditions conducted in a small annular flume at the University of Western Australia (Zhang Q. Some aspects of local scour mechanics around subsea pipelines [D]. Perth: University of Western Australia, 2015), a simulation study was conducted on the unsteady flow scour process with gradually increasing inflow velocity. In the experiment, the pipeline diameter was D = 0.05 m, and the median particle size of the bed sediment was... Reference 5's experiment observed the change in the depth of the scour pit below the pipeline center during a 1-hour local scour process. Simultaneously, ADV was used to record the change in inflow velocity at a position 8D upstream of the pipeline center. Within the first 0.5 hours, the upstream inflow velocity was approximately... acceleration from Gradually increase to However, it remains unchanged within the next 0.5 hours. A constant flow rate.

[0071] Figure 9 The computational domain layout in this simulation is shown. The distance between the pipeline center and the outflow boundary is set to 40D to ensure sufficient development of the downstream wake. Since the flow velocity is still relatively low in the initial stage of scouring, the scouring pit below the pipeline develops slowly. Therefore, a smaller initial scouring pit is needed for better comparison with experimental results; the depth of the initial scouring pit below the pipeline is set to 0.05D.

[0072] Table 1 shows a comparison of the efficiency of the method in Example 1 during parallel computing and single-core computing (i.e., the method in Example 1 does not perform parallel computing). This represents the average simulation time required to simulate a 1-minute flushing process. The table shows that the simulation time is significantly reduced in parallel computing compared to serial computing, exhibiting a near-linear acceleration characteristic.

[0073] Table 1 Comparison of computational efficiency of scour models

[0074] Figure 10A comparison of the simulated and measured flow velocities at a location 8D upstream of the pipeline center is presented. The figure shows that the inflow boundary conditions set in Example 1 can accurately simulate the gradual increase and subsequent stabilization of the inflow velocity. During the flow acceleration phase, the simulated and measured values ​​match well within the time interval t ∈ [0.1 h, 0.3 h], while the simulated velocity is slightly lower than the measured velocity within the time interval t ∈ (0.3 h, 0.5 h]. This is because the acceleration of the measured velocity change increases slightly during this time interval, while the numerical simulation maintains a constant acceleration. Within the time interval t ∈ (0.5 h, 1.0 h), the velocity remains constant. At this point, the simulated velocity is slightly lower than the average measured velocity, but the measured velocity fluctuates around the simulated value. Overall, the simulated and measured velocities exhibit a consistent trend with minimal difference in value.

[0075] Figure 11 The paper presents a comparison between the dimensionless scour depth G / D below the pipeline center and the measured values, showing a trend over scour time. It can be observed that within the time interval t ∈ [0.2 h, 0.6 h], the simulated scour depth variation agrees well with the experimental values, accurately reflecting the variation law of the scour pit depth below the pipeline center during the flow acceleration stage, and both reach a quasi-equilibrium state around t = 0.6 h. Within the time interval t ∈ (0.6 h, 1.0 h], the measured value of G / D shows a slight backfilling process, i.e., G / D decreases slightly with increasing scour time. However, in the simulation results, the G / D value remains constant and is slightly larger than the experimental value. This difference may be because the suspended sediment transport process in the annular flume reaches an equilibrium state during long-term scour experiments, while the influence of suspended sediment circulation is ignored in the numerical simulation.

[0076] It is worth noting that during the period of accelerated flow, the pipeline scouring will inevitably undergo a transition from a clear water scouring process to a moving bed scouring process. According to measured data, the initial flow velocity of sediment at a position 0.5D above the horizontal bed surface is approximately... Therefore, the critical point between clean water scouring and moving bed scouring is around t ≈ 0.28 h. Figure 11 The variation of G / D over time reveals a brief decrease in the scouring rate below the pipeline near t = 0.28 h (the location indicated by the black dashed line in the figure), followed by a gradual increase in the scouring rate after the pipeline is fully in a state of dynamic bed scouring. This phenomenon is also observed in the experimental data of reference x, demonstrating the accuracy of the simulation performed by the method in Example 1. The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for simulating local scour of a structure under unsteady flow, characterized in that, Includes the following steps: Step 1: Generate the computational domain mesh and set the boundary conditions; Step 2: Construct a finite-area projected mesh with the same planar topology as the bottom boundary mesh; Step 3: Divide the computational domain into several sub-regions based on the number of parallel solution cores; Step 4: Begin parallel solution and calculate the interpolation weights for different mesh surfaces based on the projected mesh generated in Step 2; including the following steps: Step 4.1: Calculate the sum of the reciprocals of the distances from each grid point to the center of all its adjacent faces. ; For a grid point inside the processor, calculate the reciprocal of the distance from the grid point to the adjacent face center point and sum them to obtain the result. ; For points belonging to two processors, calculate the sum of the reciprocals of the distances between the grid points and their respective processors. Then, use the Pstream class to implement data exchange between the coupling edges of adjacent processors, and sum the corresponding reciprocals of the distances to obtain the result. For points belonging to three or more processors, the globalData class in the finite area method is used to implement data transmission and coupling between multiple processors, thereby obtaining all grid points. And sum to obtain the final ; Step 4.2: Calculate the interpolation weights on different grid surfaces using the following formula. : ； In the formula, Represents the center of the i-th mesh face With grid points The reciprocal of the distance between them; Step 5: Solve the flow field. Use rough wall boundary conditions to calculate and assign values ​​to the flow velocity and turbulence variables at the bottom boundary. Use unsteady inflow boundary conditions to calculate and assign values ​​to the flow velocity and turbulence variables at the inflow boundary at different times. Step 6: Solve for bedload transport and suspended sediment transport, and calculate the displacement at the center point of the bed grid due to sediment transport at each time step; Step 7: Using parallel solving, calculate the displacement at the grid points on the bottom boundary through interpolation, and update the grid point coordinates; Step 8: Determine the slope of each grid surface. If there is a grid surface with a slope greater than the angle of repose, use the sediment slippage algorithm to correct the coordinates of the grid points on the bed surface until the slope of all grid surfaces is less than or equal to the angle of repose. Step 9: Update the computational grid and use the new grid point coordinates as the bottom boundary positions for the flow field solution at the next time step.

2. The method for simulating local scour of a structure under unsteady flow according to claim 1, characterized in that, In step 5, the flow velocity and turbulence variables at the bottom boundary are calculated and assigned values ​​using rough wall boundary conditions, including the following steps: Step 5.1: Adjust the distance from the bed surface This serves as the boundary of the bottom rough surface; Step 5.2: Calculate the flow velocity parameters at the bottom boundary using the following formula: ； in, ； ； ； In the formula, The tangential velocity representing the water flow at the bottom boundary. Represents the Karman constant. Represents the roughness height of the bed surface. This represents the height at the center point c of the first-layer grid, obtained from the governing equations of the hydrodynamic model in the single-phase flow model. The unit normal vector representing the mesh surface. This represents the flow velocity at point c. Represents the velocity of the geometric center point of the mesh surface. This represents the velocity of a grid point in the direction perpendicular to the surface normal caused by its movement. Step 5.3: Calculate the parameters of the turbulence variables at the bottom boundary using the following formula: ； ； In the formula, The value representing the turbulent kinetic energy at the bottom boundary. Represents a constant. This represents the dissipation rate of turbulent kinetic energy at the bottom boundary.

3. The method for simulating local scour of a structure under unsteady flow according to claim 2, characterized in that, In step 5, the inflow velocity and turbulence variables at different times are calculated and assigned using non-steady inflow boundary conditions, including the following steps: Step 5.4: Calculate the bed surface friction velocity at different times using the following formula: ； In the formula, The velocity of the water flow at the upstream reference position of the pipeline at different times; Step 5.5: Calculate the parameters of the inflow boundary velocity at different times using the following formula: ； In the formula, The horizontal velocity representing the water flow at the inflow boundary at different times. The maximum value in the horizontal velocity profile representing the flow motion at the inflow boundary at different times; Step 5.5: Calculate the parameters of the turbulence variables at the inflow boundary at different times using the following formula: ； ； in, ； In the formula, The distribution of turbulent kinetic energy at the inflow boundary. Represents the inflow boundary layer thickness. The distribution representing the dissipation rate of turbulent kinetic energy at the inflow boundary.

4. The method for simulating local scouring of a structure under unsteady flow according to claim 1, characterized in that, In step 6, when solving for suspended sediment transport, the suspended sediment concentration at the bottom boundary is calculated using the following formula: ； In the formula, Represents the Shields number of the bed surface. This represents the critical Shield number for sediment initiation.

5. The method for simulating local scour of a structure under unsteady flow according to claim 1, characterized in that, In step 7, the displacement at the grid points on the bottom boundary is calculated using the following formula: ； In the formula, This represents the displacement at the center point of the bed grid due to sediment transport.

6. The method for simulating local scour of a structure under unsteady flow according to claim 1, characterized in that, In step 8, the coordinates of the grid points on the bed surface are corrected, and the change in bed surface elevation is calculated using the following formula: ； In the formula, Represents the bed surface elevation. The sediment slip diffusion coefficient; in, ； ； In the formula, Let C represent the angle of repose, and C represent a constant that can range from 1.5 to 2.

4. This represents the maximum bed sediment transport rate at that moment.

Images