Erosion riprap grouting combined simulation modeling method based on random convex polyhedral stone blocks

Abstract

This invention relates to the field of bridge scour protection technology, and discloses a joint simulation modeling method for scour-induced grouting based on random convex polyhedral stones. The method includes the following steps: determining the vertices of the polyhedron, distinguishing the inner and outer surfaces of the polyhedron, modeling the polyhedron solid, analyzing the shear stress on the grouting bed, performing grouting simulation analysis, comparing different stone placement methods, and recommending a joint protection design scheme. This invention uses Python and AutoCAD to interactively model the stones, starting from improving the critical shear stress for grouting and considering the flow velocity and scour conditions around the pile foundation to achieve multi-dimensional simulation analysis. This method is easy to program, highly accurate, and easy to operate. Compared with traditional protection methods, the joint protection method improves efficiency and has greater advantages in practical applications. The influence of water depth and flow velocity on the stone protection efficiency is also analyzed, providing a basis for the design and optimization of protection schemes.

Description

Technical Field

[0001] This invention relates to the field of bridge scour protection technology, specifically to a joint simulation modeling method for scour-induced riprap grouting based on random convex polyhedral stones. Background Technology

[0002] In the field of bridge engineering today, the failure of bridge structural functions and the loss of safety performance are critical issues that urgently need to be addressed. Hydrological factors are one of the main causes of these problems, a phenomenon that has attracted widespread attention from numerous scholars. During bridge operation, the scouring effect of water flow on the area around the piers and pile foundations can severely affect the stability and safety of the bridge. Scouring can cause sediment to be carried away from the piers and pile foundations, reducing the depth of the pile foundations and decreasing the bearing capacity, thereby threatening the safety of the entire bridge structure.

[0003] Currently, bridge scour protection methods are mainly divided into two categories: active protection and passive protection.

[0004] The core idea of ​​active protection is to alter the flow field structure that causes scouring by setting specific structures or features at different locations on the bridge pier foundations, thereby reducing the erosive effect of water flow on sediment. Among various active protection methods, retaining rings are particularly effective. Studies have shown that the protective effect is optimal when the retaining ring is placed flush with the general scouring surface of the riverbed. However, active protection methods may be limited by various factors in practical applications.

[0005] Passive protection primarily focuses on increasing the initial shear stress of the bed. Common passive protection measures include riprap placement, gabions, concrete blocks, and other riprap alternatives, as well as expanded foundation protection. In practical engineering, riprap placement is one of the most commonly used passive protection measures. However, this method has some drawbacks. Directly placed riprap has relatively weak stability and integrity, and is prone to displacement and dispersion under long-term water flow. Therefore, secondary protection is usually required, which increases engineering costs and construction difficulty.

[0006] As research deepens, researchers both domestically and internationally have begun exploring combined protection methods based on traditional single scour protection approaches. The main idea behind this approach is to combine two or more protection methods, hoping to achieve a superior protective effect compared to a single method. These studies have provided new ideas and directions for bridge scour protection, but they also face numerous challenges.

[0007] Based on the current research status of bridge scour protection methods, there is an urgent need for a simulation modeling method for scour protection using a combination of riprap grouting and random convex polyhedral stones. Summary of the Invention

[0008] The purpose of this invention is to provide a joint simulation modeling method for scour-induced rock grouting based on random convex polyhedral stones, which solves the technical problems mentioned in the background art.

[0009] The objective of this invention can be achieved through the following technical solutions:

[0010] The method for joint simulation modeling of scour-filled grouting based on random convex polyhedral stones includes the following steps:

[0011] Step A1: Determining the vertices of the polyhedron

[0012] Random convex polyhedral stone modeling based on Python and AutoCAD interaction

[0013] Step A1.1: Determine the vertices of the polyhedron using spherical coordinates, as shown in the following formula:

[0014]

[0015] Obtain the vertex coordinates (x, y, z) of the polyhedron;

[0016] Where θ and φ are the key parameters for determining the vertex coordinates, r represents the distance from the vertex to the origin, θ is the elevation angle, which ranges from 0 to π, and φ is the azimuth angle, which ranges from 0 to 2π.

[0017] Step A1.2: Determine whether interference occurs between the vertices by judging the distance between the vertices of the polyhedron;

[0018] The formula is as follows:

[0019]

[0020] In the formula, (x1, y1, z1) and (x2, y2, z2) are the coordinates of the two vertices, and L is the distance between the two vertices;

[0021] Compare L with the preset interference distance threshold Lmin:

[0022] If L > Lmin, it means that there is no interference between the two vertices; otherwise, it means that there is interference between the two vertices.

[0023] Step A2: Distinguishing the inner and outer surfaces of the polyhedron

[0024] Distinguish the outer and inner surfaces of a polyhedron from the generated random points;

[0025] First, let plane abc, formed by any three points in space, define the outward normal vector of this plane as...

[0026] Its outward normal direction vector The calculation formula is as follows:

[0027]

[0028] in,

[0029]

[0030] Iterate through all vertices of the polyhedron, and define the direction vector connecting each vertex to any point in plane abc as follows:

[0031] The vertex points to the direction vector connecting to any point in plane abc. The calculation formula is as follows:

[0032]

[0033] If it exists If the value is less than 0, that is, the angle between the two vectors is greater than 90°, it means that all points in the space are not on the same side of the plane, and the plane does not meet the convexity condition and is not the outer surface of the polyhedron, so the plane is discarded.

[0034] If any If the values ​​of the two vectors are all greater than 0, that is, the angle between the two vectors is less than 90°, it means that all points in the space are on the same side of the plane. The plane meets the convexity condition and is defined as the outer surface of the polyhedron.

[0035] Step A3: Polyhedral Solid Modeling

[0036] After determining the outer surface of the polyhedron, multiple sets of data files containing three points are obtained;

[0037] Determine the output order of the three points on the outer surface, and output the coordinates of the three points in a clockwise direction;

[0038] The purpose is to ensure that when a sectioning command is executed on a sphere, the inner part containing the polyhedron is preserved, while the outer part of the polyhedron is discarded;

[0039] The specific steps to determine the clockwise output order of the three points are as follows:

[0040] Step A3.1: Select three points a, b, and c on the outer surface of the polyhedron and calculate...

[0041] Step A3.2, Calculation and The vector product, denoted as

[0042] Step A3.3, Calculation and The dot product;

[0043] If their dot product is less than 0, it indicates that a, b, c are in clockwise order;

[0044] Conversely, if a, b, c are in a counter-clockwise order, then output them again in the order of a, c, b.

[0045] After determining the output order of the three points on the outer surface, the clockwise output data of the three points on the outer surface of the polyhedron is obtained. The sphere is then cut using the three-point sectioning command in CAD.

[0046] In this process, after each cut, Boolean operations are used to delete the outer surface after the cut, while keeping the inner surface. After traversing all the cutting point sets, the desired polyhedron can be obtained.

[0047] Step A4: Bed Surface Shear Stress Analysis

[0048] When the shear stress in the bed exceeds the preset shear stress threshold, the eddy current will carry the particles away from the riverbed surface and suspend them. The suspended sediment particles are carried by the water flow.

[0049] When suspended sediment particles settle onto the sand bed under the combined action of gravity, buoyancy and friction, sediment deposition occurs.

[0050] Among them, the shear stress threshold is expressed as the critical shear stress;

[0051] Step A4.1: Bed shear stress is the shear stress exerted by the fluid on the surface of the filling bed, which is calculated based on the wall function of three-dimensional turbulence.

[0052] The specific formula is as follows:

[0053]

[0054] In the formula, u τ Shear rate, τ is the bed shear stress; ρ is the bulk density of the water and sediment mixture; Y is the distance between the center point of the element grid and the boundary wall; v is the kinematic viscosity of the fluid; k = 0.4 is the von Kármán constant. s For roughness, k s =c s d 50 d 50 c is the median particle size of the sediment. s In this embodiment, c is a custom coefficient. s The value is 2.5;

[0055] In shallow water flow, the second friction law is used to calculate the shear stress in the riverbed. In the formula ρ is the depth-averaged fluid velocity.f C is the fluid density. D The drag coefficient is determined by the roughness k. s Local determination:

[0056] in,

[0057]

[0058] In the formula, h is the water depth, B = 0.71; z0 = k s / 30;

[0059] Step A4.2, the Shields number is the dimensionless form of the bed shear stress, defined as:

[0060]

[0061] In the formula, θ n ρ is the Shields number, g is the absolute value of gravity; n d represents the mass density of the sediment particles. n The diameter of the particle;

[0062] First, the physical processes related to bed shear stress are explained, then the calculation methods are introduced, including calculations based on the three-dimensional turbulent wall function and the secondary friction law in shallow water flow. Finally, the dimensionless form of the Shields number is given. These are interconnected and jointly describe the characteristics of bed shear stress and its influence on sediment movement.

[0063] Step A5: Grouting Simulation Analysis

[0064] Step A5.1: The grouting section is simulated and analyzed using the critical Shields number for soil improvement.

[0065] The specific method is as follows:

[0066] The critical shear stress of scour for known silt-stabilized soil is extracted, and its value is between 2.5 and 4 Pa. Then, the equivalent critical Shield number of the grouting soil is calculated by back-calculating the critical Shield number formula.

[0067] The scour prevention effect after grouting is then evaluated using weighted shear stress, and the formula is as follows:

[0068]

[0069] Where A wcz The shear stress is the area-weighted value, where τ is the bed surface shear stress. cr The bed surface is subjected to boundary shear stress;

[0070] Step A5.2 Simultaneously, extract the overshear stress on the bed surface, which is used to reflect the shear stress distribution around the pile foundation, and then to predict the scour development trend.

[0071] Among them, the excess shear stress of the bed surface is expressed as the part where the bed surface shear stress exceeds the critical shear stress of sediment;

[0072] The method for extracting the excess shear stress value per unit area is as follows:

[0073]

[0074] In the formula, is the average value of the excess shear stress, which is used to measure the protection effect; τ i is the shear stress value at a certain grid point on the bed surface; A i is the grid area;

[0075] StepA5.3. Conduct relevant analyses based on the fixed-bed scour model for the flow around stones;

[0076] The specific method is as follows:

[0077] Extract the excess shear stress data within the 4D range in front, behind, left, and right of the pile foundation. That is, within the nested grid, and use it as the basis for subsequent excess shear stress analysis;

[0078] StepA6. Comparison of different stone placement methods

[0079] Compare and analyze the protection effects of two stone placement methods: staggered circular arrangement and random arrangement

[0080] StepA6.1. For the staggered circular arrangement, the method is as follows:

[0081] The stones are distributed in three layers. The interval angle between adjacent stones in each layer is 30°, and the two adjacent layers are staggered by a rotation angle of 15°;

[0082] StepA6.2. For the random arrangement, the method is as follows:

[0083] StepA6.2.1. Randomly place polyhedrons within the set protection area. Consider placing them within a circular ring area, and set the center point of placement using polar coordinates:

[0084]

[0085] In the formula, R1 and R2 are the two radius values used to define the placement area in the polar coordinate setting. Among them, R1 is the inner radius and R2 is the outer radius, indicating that when the polyhedron is placed within the circular ring area, this area is the ring range determined by the radii R1 and R2; ρ is the radial variable in the polar coordinate, representing the distance from the pole to a certain point. Among them, the value range of ρ is R1 < ρ < R2, that is, the distance from the placement point of the polyhedron to the center of the placement area should be between R1 and R2, and the pole is represented as the center of the placement area;

[0086] Step A6.2.2: At the same time, the spacing between the centers of any two polyhedra is controlled within a specified range by using a preset lower limit value for spacing. The purpose is to ensure that the polyhedra are distributed as evenly as possible around the cylindrical pile foundation.

[0087] The formula for calculating the spacing is as follows:

[0088]

[0089] In the formula L1 min This is a preset lower limit for the spacing, representing the minimum allowable distance between the center points of two polyhedra, and L1 min =3r1, where r1 is the stone particle size, (xb1, yb1, zb1) and (xb2, yb2, zb2) are the coordinates of the two polyhedra respectively, and L1 is the distance between the center points of the two polyhedra;

[0090] Step A6.2.3: At the center point coordinates (xb) of the generated polyhedron i ,yb i ,zb i After that, a random polyhedron is first generated at the origin (0, 0, 0), and then the random polyhedron is distributed around the cylinder by translation, thus completing the random arrangement of the random polyhedron.

[0091] Step A7: Recommended Design Scheme for Joint Protection

[0092] Perform statistical analysis on the overshear stress mesh data within the computational domain;

[0093] The specific method is as follows:

[0094] The statistical data of the overshear stress cloud map were divided into 12 data intervals: [0, 0.1], [0.1, 0.2], [0.2, 0.3], ..., [1.1, 1.2];

[0095] For each data interval, calculate the corresponding frequency and probability density; the frequency calculation formula is:

[0096]

[0097] Where n i Let N be the number of data points in the i-th data interval, and N be the total number of statistical data points; the probability density formula is:

[0098]

[0099] Where L represents the length of the data interval, and L is set to 0.1 in the data analysis.

[0100] The purpose is to investigate the distribution of overshear stress around the pile foundation under different stone blocks and grouting parameters in order to derive the optimal combined protection recommended design scheme;

[0101] As a further aspect of the present invention: in Step A4, fluid analysis calculations are also performed based on Flow-3D;

[0102] Flow-3D uses the continuity equation and the Navier-Stokes equation for incompressible viscous fluids as the governing equations for fluid motion.

[0103] Considering that the fluid is a Newtonian fluid, and combined with Flow-3D's unique FAVOR mesh technology, the continuity equation and momentum equation include volume and area fraction parameters;

[0104] As a further aspect of this invention: In Step A4, the equations involved in calculating the bed shear stress are solved using the generalized minimum residual GMRES algorithm.

[0105] The solution method is as follows:

[0106] First, an intermediate velocity is introduced. The influence of pressure at the new moment on the velocity field is not considered. Instead, the pressure correction value at the current moment is introduced.

[0107] The intermediate velocity can be obtained by solving the momentum equation, which is the equation describing the change of fluid momentum.

[0108] Then, the relationship between the obtained intermediate velocity and the pressure correction value is substituted into the continuity equation to generate the pressure Poisson equation containing the pressure correction value.

[0109] Finally, the GMRES algorithm is applied to solve the pressure Poisson equation.

[0110] As a further aspect of the present invention: in Step A5, the critical Shield number is used to determine whether sediment will initiate movement from a dimensionless perspective;

[0111] The critical Shield number is used to characterize the critical bed shear stress τ. cr,n The dimensionless form, when the bed shear stress τ≥τ cr,n At that time, the sediment began to move;

[0112] The specific calculation is as follows:

[0113]

[0114] Where θ cr,n Let θ be the critical Shields number. cr,n The calculation equation is determined based on the properties of the sediment and is as follows:

[0115]

[0116] Where, d *,n For dimensionless particle size, d *,n It reflects the influence of factors such as sediment particle diameter, gravitational acceleration, density difference between sediment and fluid, and fluid kinematic viscosity on sediment movement, and is a key intermediate parameter for calculating the critical Shields number.

[0117] d *,n The formula for calculation is:

[0118]

[0119] Where, d n The diameter of the sediment particles is represented by s, g is the acceleration due to gravity, and s is the velocity of the sediment particles. n =ρ n / ρ f , ρ n ρ is the mass density of the sediment particles. f For fluid density, s n V represents the relative density of sediment and fluid. f It is the kinematic viscosity of a fluid, which represents the fluid's ability to resist deformation. It is related to the fluid's viscosity. The smaller the value, the easier the fluid flows, and the force exerted on the sediment particles will change accordingly, thus affecting the movement state of the sediment.

[0120] As a further aspect of the present invention: in Step A5, the net exchange rate between accumulated and suspended sediments is determined by the entrainment velocity and deposition velocity of the sediment.

[0121] Among them, the entrainment velocity and deposition velocity of sediment are regarded as two opposite micro-processes that occur simultaneously and are related to each other.

[0122] The sediment entrainment velocity can be calculated using the following formula:

[0123]

[0124] Where, α n Let n be the entrainment coefficient of the nth type of sediment. b The normal vector of the filling bed surface is outward;

[0125] The sediment settling velocity can be calculated using the following formula:

[0126]

[0127] Where g is the acceleration due to gravity, and u settle,n They are in the same direction.

[0128] As a further aspect of the present invention: In Step A5, the dimensionless form of the sediment transport rate of the nth layer is defined as follows:

[0129]

[0130] Where, Φ n Let q be the sediment transport rate of the nth layer, which is used to comprehensively characterize the intensity and characteristics of sediment transport, facilitating a unified comparison and analysis of sediment transport under different sediment and flow conditions. b,n q represents the unit width volumetric sediment transport rate of the bed surface, which refers to the volume of sediment passing through a unit width of the bed surface per unit time. It directly reflects the amount of sediment transported on the bed surface and is an important physical quantity for measuring the sediment transport process. b,n The calculation formula is as follows:

[0131] q bn =B n (θ n -θ c,n ) 1.5 c b,n

[0132] Among them, B n The sediment transport coefficient of the bed layer;

[0133] Among them, the low-speed transport coefficient is usually taken as 5.0-5.7, about 8.0 for medium transport speed, and up to 13.0 for high transport speed;

[0134] c bn Let n be the volume fractions of the n substances in the bed mass:

[0135]

[0136] As a further aspect of the present invention: in Step A7, the expected value of the overshear stress is also calculated based on the probability density statistical diagram of the overshear stress.

[0137] The expected value is calculated as follows:

[0138]

[0139] Where E(X) is the expected value of the overshear stress, X i Let P(X) be the midpoint of the i-th data interval. i ) represents the probability corresponding to the i-th data interval, and n is the total number of data intervals;

[0140] By integrating the information from the probability density graph of overshear stress, an expected value that can represent the average level of overshear stress is obtained, thus providing a more intuitive understanding of the overall distribution characteristics of overshear stress.

[0141] The beneficial effects of this invention are:

[0142] The proposed simulation modeling method for combined scour protection using riprap and grouting based on random convex polyhedral stones has significant advantages. A modeling method for random convex polyhedral stones based on Python and AutoCAD interaction was developed. This innovative modeling approach is easy to implement programmatically. Complex polyhedral stone models can be efficiently generated through programming, reducing the workload and errors of manual modeling. This provides an accurate model foundation for subsequent scour protection analysis, greatly improving the efficiency and operability of the entire research process.

[0143] This invention achieves comprehensive simulation analysis with high accuracy from multiple key dimensions. Firstly, starting with improving the critical initiation shear stress of the grouting section, it accurately simulates and analyzes the grouting area around the pile foundation based on clean water scouring. This allows for a deeper understanding of the grouting's role in the protection process and its response to different water flow conditions, thus enabling better optimization of grouting parameters. Secondly, it analyzes the flow velocity around the pile foundation, effectively simulating and analyzing the disturbance effect of rocks based on fixed-bed scouring. This method clearly understands the role of rocks in water flow and how to reduce scouring of sediment around the pile foundation by changing the water flow state. Furthermore, starting with the scouring situation around the pile foundation, it achieves comprehensive simulation analysis of the combined protection based on clean water scouring, fully considering the interaction of various factors and providing a reliable basis for evaluating the effectiveness of the combined protection.

[0144] Based on analysis of factors such as flow velocity, scour pit morphology, and scour depth variations, this invention demonstrates significant advantages over traditional protection methods such as sacrificial piles, retaining rings, and riprap placement in its combined rock-block protection system. The combined protection method significantly improves protection efficiency and makes it more competitive in practical applications. This advantage is reflected in the long-term stability protection of bridge structures, effectively reducing maintenance costs and extending bridge service life.

[0145] This invention, through fixed-bed simulation and using the change in the average overshear stress as an indicator, analyzes in depth the impact of water depth and flow velocity on the efficiency of rock protection. The results show that water depth has little effect on the average overshear stress around the pile foundation, providing a reference for protection design under different water depth conditions. Simultaneously, the relationship between flow velocity and rock protection efficiency is clarified: when the flow velocity is low, the rock protection efficiency increases with increasing flow velocity, but once the flow velocity reaches a certain value, the improvement in rock protection efficiency becomes no longer significant. This conclusion helps to rationally select and optimize protection schemes based on actual water flow velocities, further improving the pertinence and effectiveness of protection measures. Attached Figure Description

[0146] The invention will now be further described with reference to the accompanying drawings.

[0147] Figure 1 This is a flowchart illustrating the simulation modeling method for scour-induced rock grouting combined protection based on random convex polyhedral stones according to the present invention.

[0148] Figure 2 This is a schematic diagram of the pre-selected protection range and protection corner coordinates of the multi-pile piers in this invention;

[0149] Figure 3 This is a schematic diagram of the planar layout and top rough surface of the multi-pile pier protection system of the present invention;

[0150] Figure 4 This is a schematic diagram of the three-dimensional convex polyhedron of the present invention;

[0151] Figure 5 This is a schematic diagram illustrating the staggered circular and random placement of stones used in this invention;

[0152] Figure 6 This is a streamline diagram of the area around a single pile foundation under unprotected and combined protection conditions used in this invention;

[0153] Figure 7 This is a schematic diagram of the scour morphology of multi-pile bridge piers with and without combined protection, as used in this invention.

[0154] In the figure, (x1) rect1 ,y1 rect1 (x1) rect1 ,y2 rect1 (x2) rect1 ,y2 rect1 (x2) rect1 ,y1 rect1 (x1) represents the coordinates of the interior corner of the protected area; rect2 ,y1 rect2 (x1) rect2 ,y2 rect2 (x2) rect2 ,y2 rect2 (x2) rect2 ,y1 rect2 1) Coordinates of the outer corner of the protection range; 4) Location of the pile foundation; 5) Location of the pier column; 6) Top rough surface of the riprap protection; 7) Sphere; 8) Partially cut sphere; 9) Polyhedron; 10) Staggered ring placement; 11) Random placement; 12) Morphology and streamline of a single pile without protection scour pit; 13) Morphology and streamline of a single pile with combined protection scour pit; 14) Scour morphology of a multi-pile foundation pier without combined protection; 15) Scour morphology of a multi-pile foundation pier with combined protection.

[0155] Figures 2 to 7 It serves to visually demonstrate the invention's content and verify its effects throughout the text;

[0156] pass Figures 2 to 5It can intuitively demonstrate key technical points such as protection range, structural layout, stone modeling and placement method, and can more clearly understand each link of the simulation modeling method of scour grouting and grouting combined protection based on random convex polyhedral stones.

[0157] Figure 6 and Figure 7 By comparing the scour patterns and depth data under unprotected and combined protection conditions, the effectiveness of the combined protection method is strongly demonstrated. Detailed Implementation

[0158] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0159] exist Figure 2 In the process, based on the protection range and the coordinates of the corner points of the protection, the coordinates of the center point of the stones within the protection range are determined;

[0160] The method for determining it is as follows:

[0161] (x∈(x1 rect1 ,x1 rect2 )∩(x2 rect2 x2 rect1 ),y∈(y1 rect1 ,y2 rect1 ))∪(x∈

[0162] (x1 rect1 x2 rect1 ),y∈(y1 rect1 y1 rect2 )∩(y2 rect2 ,y2 rect1 ))

[0163] exist Figure 6 In the diagram of the scour pit morphology, under the rockfill protection, due to the high flow velocity on both sides of the pile foundation, the water flow carries and transports the rocks, resulting in some scour on both sides of the pile foundation. However, the overall scour depth is lower than that without protection measures. At the end of 600s, the depth of the scour pit on both sides is 1.166cm, which is lower than the scour depth (2.939cm) without protection measures. The maximum scour depth reduction rate is 60.3%.

[0164] exist Figure 7 In the case of combined protection, after 30 minutes of scouring, the maximum scouring depth was 0.63 m; without combined protection, the maximum scouring depth after 30 minutes of scouring was [missing value].

[0165] After applying combined protection, the scour depth on both sides and behind the pile cap was improved, reaching 1.012m.

[0166] The specific implementation method of this invention takes Flow3D as an example, and develops a random convex polyhedral stone block modeling method based on the interaction of Python and AutoCAD, and establishes a joint protection scour simulation model based on random convex polyhedral stones.

[0167] Example 1

[0168] Please see Figure 1 This invention relates to a joint simulation modeling method for scour-induced riprap grouting based on random convex polyhedral stones, including a random convex polyhedral stone model constructed using Python and AutoCAD interaction, comprising the following steps:

[0169] Step A1: Determining the vertices of the polyhedron

[0170] Random convex polyhedral stone block modeling based on Python and AutoCAD interaction;

[0171] Step A1.1: Determine the vertices of the polyhedron using spherical coordinates, as shown in the following formula:

[0172]

[0173] Obtain the vertex coordinates (x, y, z) of the polyhedron;

[0174] Where θ and φ are the key parameters for determining the vertex coordinates, r represents the distance from the vertex to the origin, θ is the elevation angle, which ranges from 0 to π, and φ is the azimuth angle, which ranges from 0 to 2π.

[0175] Step A1.2: Determine whether interference occurs between the vertices by judging the distance between the vertices of the polyhedron;

[0176] The formula is as follows:

[0177]

[0178] In the formula, (x1, y1, z1) and (x2, y2, z2) are the coordinates of the two vertices, and L is the distance between the two vertices;

[0179] Compare L with the preset interference distance threshold Lmin:

[0180] If L > Lmin, it means that there is no interference between the two vertices; otherwise, it means that there is interference between the two vertices.

[0181] Step A2: Distinguishing the inner and outer surfaces of the polyhedron

[0182] Random convex polyhedral stone modeling based on Python and AutoCAD interaction

[0183] Distinguish the outer and inner surfaces of a polyhedron from the generated random points;

[0184] First, let plane abc, formed by any three points in space, define the outward normal vector of this plane as...

[0185] Its outward normal direction vector The calculation formula is as follows:

[0186]

[0187] in,

[0188]

[0189] Iterate through all vertices of the polyhedron, and define the direction vector connecting each vertex to any point in plane abc as follows:

[0190] The vertex points to the direction vector connecting to any point in plane abc. The calculation formula is as follows:

[0191]

[0192] If it exists If the value is less than 0, that is, the angle between the two vectors is greater than 90°, it means that all points in the space are not on the same side of the plane, and the plane does not meet the convexity condition and is not the outer surface of the polyhedron, so the plane is discarded.

[0193] If any If the values ​​of the two vectors are all greater than 0, that is, the angle between the two vectors is less than 90°, it means that all points in the space are on the same side of the plane. The plane meets the convexity condition and is defined as the outer surface of the polyhedron.

[0194] Step A3: Polyhedral Solid Modeling

[0195] Random convex polyhedral stone modeling based on Python and AutoCAD interaction

[0196] After determining the outer surface of the polyhedron, multiple sets of data files containing three points are obtained;

[0197] Determine the output order of the three points on the outer surface, and output the coordinates of the three points in a clockwise direction;

[0198] The purpose is to ensure that when a sectioning command is executed on a sphere, the inner part containing the polyhedron is preserved, while the outer part of the polyhedron is discarded;

[0199] The specific steps to determine the clockwise output order of the three points are as follows:

[0200] Step A3.1: Select three points a, b, and c on the outer surface of the polyhedron and calculate...

[0201] Step A3.2, Calculation and The vector product, denoted as

[0202] Step A3.3, Calculation and The dot product;

[0203] If their dot product is less than 0, it indicates that a, b, c are in clockwise order;

[0204] Conversely, if a, b, c are in a counter-clockwise order, then output them again in the order of a, c, b.

[0205] After determining the output order of the three points on the outer surface, the clockwise output data of the three points on the outer surface of the polyhedron is obtained. The sphere is then cut using the three-point sectioning command in CAD.

[0206] In this process, after each cut, Boolean operations are used to delete the outer surface after the cut, while keeping the inner surface. After traversing all the cutting point sets, the desired polyhedron can be obtained.

[0207] This embodiment realizes the random convex polyhedron stone model based on the interaction of Python and AutoCAD. By determining the vertices of the polyhedron, distinguishing the inner and outer surfaces, and solid modeling, a foundation is laid for the subsequent protection model. The stone model can be accurately constructed for research related to erosion protection.

[0208] Example 2

[0209] As a second embodiment of the present invention, in specific implementation, the technical solution of this embodiment differs from that of embodiment one only in that this embodiment further includes the following steps:

[0210] Step A4: Bed Surface Shear Stress Analysis

[0211] Step A4.1: Bed shear stress is the shear stress exerted by the fluid on the surface of the filling bed, which is calculated based on the wall function of three-dimensional turbulence.

[0212] The specific formula is as follows:

[0213]

[0214] In the formula, u τ Shear rate, τ is the bed shear stress; ρ is the bulk density of the water and sediment mixture; Y is the distance between the center point of the element grid and the boundary wall; v is the kinematic viscosity of the fluid; k = 0.4 is the von Kármán constant. s For roughness, k s =c s d 50 d 50 c is the median particle size of the sediment. s In this embodiment, c is a custom coefficient. s The value is 2.5;

[0215] In shallow water flow, the second friction law is used to calculate the shear stress in the riverbed. In the formula ρ is the depth-averaged fluid velocity. f C is the fluid density. D The drag coefficient is determined by the roughness k. s Local determination:

[0216] in,

[0217]

[0218] In the formula, h is the water depth, B = 0.71; z0 = k s / 30;

[0219] The calculation of bed shear stress involves solving the relevant equations using the generalized minimum residual GMRES algorithm.

[0220] The solution method is as follows:

[0221] First, an intermediate velocity is introduced. The influence of pressure at the new moment on the velocity field is not considered. Instead, the pressure correction value at the current moment is introduced.

[0222] The intermediate velocity can be obtained by solving the momentum equation, which is the equation describing the change of fluid momentum.

[0223] Then, the relationship between the obtained intermediate velocity and the pressure correction value is substituted into the continuity equation to generate the pressure Poisson equation containing the pressure correction value.

[0224] Finally, the GMRES algorithm is applied to solve the pressure Poisson equation;

[0225] Step A4.2, the Shields number is the dimensionless form of the bed shear stress, defined as:

[0226]

[0227] In the formula, θ nρ is the Shields number, g is the absolute value of gravity; n d represents the mass density of the sediment particles. n The diameter of the particle;

[0228] When the shear stress in the bed exceeds the preset shear stress threshold, the eddy current will carry the particles away from the riverbed surface and suspend them. The suspended sediment particles are carried by the water flow.

[0229] When suspended sediment particles settle onto the sand bed under the combined action of gravity, buoyancy and friction, sediment deposition occurs.

[0230] Among them, the shear stress threshold is expressed as the critical shear stress;

[0231] First, the physical processes related to bed shear stress are explained, then the calculation methods are introduced, including calculations based on the three-dimensional turbulent wall function and the secondary friction law in shallow water flow. Finally, the dimensionless form of the Shields number is given. These are interconnected and jointly describe the characteristics of bed shear stress and its influence on sediment movement.

[0232] In this step, fluid analysis calculations are also performed based on Flow-3D;

[0233] Flow-3D uses the continuity equation and the Navier-Stokes equation for incompressible viscous fluids as the governing equations for fluid motion.

[0234] Considering that the fluid is a Newtonian fluid, and combined with Flow-3D's unique FAVOR mesh technology, the continuity equation and momentum equation include volume and area fraction parameters;

[0235] Step A5: Grouting Simulation Analysis

[0236] Step A5.1: The grouting section is simulated and analyzed using the critical Shields number for soil improvement.

[0237] The specific method is as follows:

[0238] The critical shear stress of scour for known silt-stabilized soil is extracted, and its value is between 2.5 and 4 Pa. Then, the equivalent critical Shield number of the grouting soil is calculated by back-calculating the critical Shield number formula.

[0239] The scour prevention effect after grouting is then evaluated using weighted shear stress, and the formula is as follows:

[0240]

[0241] Where A wcz The shear stress is the area-weighted value, where τ is the bed surface shear stress. cr The bed surface is subjected to boundary shear stress;

[0242] Among them, the critical Shield number is used to determine whether sediment will initiate movement from a dimensionless perspective;

[0243] The critical Shield number is used to characterize the critical bed shear stress τ. cr,n The dimensionless form, when the bed shear stress τ≥τ cr,n At that time, the sediment began to move;

[0244] The specific calculation is as follows:

[0245]

[0246] Where θ cr,n Let θ be the critical Shields number. cr,n The calculation equation is determined based on the properties of the sediment and is as follows:

[0247]

[0248] Where, d *,n For dimensionless particle size, d *,n It reflects the influence of factors such as sediment particle diameter, gravitational acceleration, density difference between sediment and fluid, and fluid kinematic viscosity on sediment movement, and is a key intermediate parameter for calculating the critical Shields number.

[0249] d *,n The formula for calculation is:

[0250]

[0251] Where, d n The diameter of the sediment particles is represented by s, g is the acceleration due to gravity, and s is the velocity of the sediment particles. n =ρ n / ρ f , ρ n ρ is the mass density of the sediment particles. f For fluid density, s n V represents the relative density of sediment and fluid. f It is the kinematic viscosity of a fluid, which represents the fluid's ability to resist deformation. It is related to the fluid's viscosity. The smaller the value, the easier the fluid flows, and the force exerted on the sediment particles will change accordingly, thus affecting the movement state of the sediment.

[0252] The dimensionless form of the sediment transport rate in the nth layer is defined as follows:

[0253]

[0254] Where, Φ nLet q be the sediment transport rate of the nth layer, which is used to comprehensively characterize the intensity and characteristics of sediment transport, facilitating a unified comparison and analysis of sediment transport under different sediment and flow conditions. b,n q represents the unit width volumetric sediment transport rate of the bed surface, which refers to the volume of sediment passing through a unit width of the bed surface per unit time. It directly reflects the amount of sediment transported on the bed surface and is an important physical quantity for measuring the sediment transport process. b,n The calculation formula is as follows:

[0255] q bn =B n (θ n -θ c,n ) 1.5 c b,n

[0256] Among them, B n The sediment transport coefficient of the bed layer;

[0257] Among them, the low-speed transport coefficient is usually taken as 5.0-5.7, about 8.0 for medium transport speed, and up to 13.0 for high transport speed;

[0258] c bn Let n be the volume fractions of the n substances in the bed mass:

[0259]

[0260] Step A5.2 Simultaneously, extract the overshear stress on the bed surface, which is used to reflect the shear stress distribution around the pile foundation, and then to predict the scour development trend.

[0261] Among them, the overshear stress of the bed surface is expressed as the portion of the bed surface shear stress that exceeds the critical shear stress of the sediment.

[0262] The method for extracting the overshear stress value per unit area is as follows:

[0263]

[0264] In the formula, τ represents the mean overshear stress, which is used to measure the effectiveness of the protection. i A represents the shear stress value at a certain grid point on the bed surface; i The area is the grid area.

[0265] Step A5.3: Conduct relevant analysis on the fixed-bed scour model based on the analysis of flow around stones;

[0266] The specific method is as follows:

[0267] Extract overshear stress data in a 4D range around the pile foundation, including front, back, left, and right sides. This data is then used as the basis for subsequent overshear stress analysis within the nested mesh.

[0268] In this embodiment, the bed shear stress analysis is added, and its calculation method and related physical processes are clarified, including calculation based on three-dimensional turbulent wall functions and the quadratic friction law of shallow water flow. The Shields number is also introduced, and fluid analysis is carried out in combination with Flow-3D. At the same time, the anti-erosion effect is evaluated in the grouting simulation analysis, which can more comprehensively analyze sediment movement and protection effects.

[0269] Embodiment III

[0270] As Embodiment III of the present invention, in the specific implementation of this application, compared with Embodiment I and Embodiment II, the technical solution of this embodiment is to combine and implement the solutions of the above Embodiment I and Embodiment II. The difference between the technical solution of this embodiment and Embodiment I and Embodiment II is only that this embodiment also includes a comparison of different stone placement methods;

[0271] Compare and analyze the protection effects of two stone placement methods, namely staggered annular arrangement and random arrangement, specifically as follows:

[0272] StepA6.1. For the staggered annular arrangement, the method is as follows:

[0273] The stones are distributed in three layers in total. The interval angle between adjacent stones in each layer is 30°, and the two adjacent layers are staggered by a rotation angle of 15°;

[0274] StepA6.2. For the random arrangement, the method is as follows:

[0275] StepA6.2.1. Randomly place polyhedrons within the set protection area. Considering placing them within a circular ring area, use polar coordinates to set the center point of placement:

[0276]

[0277] In the formula, R1 and R2 are two radius values used to define the placement area in the polar coordinate setting; where, R1 is the inner radius and R2 is the outer radius, indicating that when the polyhedron is placed within the circular ring area, this area is the ring range determined by the radii R1 and R2; ρ is the radial variable in the polar coordinates, representing the distance from the pole to a certain point; where, the value range of ρ is R1 < ρ < R2, that is, the distance from the placement point of the polyhedron to the center of the placement area should be between R1 and R2, and the pole represents the center of the placement area;

[0278] StepA6.2.2. At the same time, control the distance between the centers of any two polyhedrons within a specified range through a preset lower limit value of the spacing, and the purpose is to ensure that the polyhedrons are as evenly dispersed around the cylindrical pile foundation as possible;

[0279] The spacing calculation formula is as follows:

[0280]

[0281] In the formula L1 min This is a preset lower limit for the spacing, representing the minimum allowable distance between the center points of two polyhedra, and L1 min =3r1, where r1 is the stone particle size, (xb1, yb1, zb1) and (xb2, yb2, zb2) are the coordinates of the two polyhedra respectively, and L1 is the distance between the center points of the two polyhedra;

[0282] Step A6.2.3: At the center point coordinates (xb) of the generated polyhedron i ,yb i ,zb i After that, a random polyhedron is first generated at the origin (0, 0, 0), and then the random polyhedron is distributed around the cylinder by translation, thus completing the random arrangement of the random polyhedron.

[0283] This embodiment compares the protective effects of two stone placement methods: staggered circular arrangement and random arrangement. By specifying the arrangement method and parameter settings, a better stone placement scheme can be found to improve the protection efficiency.

[0284] Example 4

[0285] As a fourth embodiment of the present invention, in specific implementation, compared with embodiments one, two and three, the difference between this embodiment and embodiments one, two and three is only that this embodiment also includes the step of: joint protection recommended design scheme;

[0286] The specific method is as follows:

[0287] Perform statistical analysis on the overshear stress mesh data within the computational domain;

[0288] The specific method is as follows:

[0289] The statistical data of the overshear stress cloud map were divided into 12 data intervals: [0, 0.1], [0.1, 0.2], [0.2, 0.3], ..., [1.1, 1.2];

[0290] For each data interval, calculate the corresponding frequency and probability density; the frequency calculation formula is:

[0291]

[0292] Where n i Let N be the number of data points in the i-th data interval, and N be the total number of statistical data points; the probability density formula is:

[0293]

[0294] Where L represents the length of the data interval, and L is set to 0.1 in the data analysis.

[0295] The purpose is to investigate the distribution of overshear stress around the pile foundation under different stone blocks and grouting parameters in order to derive the optimal combined protection recommended design scheme;

[0296] In this step, the expected value of the overshear stress is also calculated based on the probability density graph of the overshear stress.

[0297] The expected value is calculated as follows:

[0298]

[0299] Where E(X) is the expected value of the overshear stress, X i Let P(X) be the midpoint of the i-th data interval. i ) represents the probability corresponding to the i-th data interval, and n is the total number of data intervals;

[0300] By integrating the information from the probability density graph of overshear stress, an expected value that can represent the average level of overshear stress is obtained, thus providing a more intuitive understanding of the overall distribution characteristics of overshear stress.

[0301] This embodiment performs statistical analysis on the overshear stress grid data, calculates the frequency, probability density and expected value, and can explore the distribution of overshear stress around the pile foundation under different parameters, which helps to derive the optimal joint protection recommended design scheme.

[0302] Example 5

[0303] As a fifth embodiment of the present invention, in specific implementation, compared with embodiments one, two, three and four, the technical solution of this embodiment is to combine the solutions of embodiments one, two, three and four.

[0304] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters and thresholds in the formulas are set by those skilled in the art according to the actual situation.

[0305] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A joint simulation modeling method for scour-induced riprap grouting based on random convex polyhedral stones, characterized in that, Includes the following steps: Step A1: Determining the vertices of the polyhedron First, the vertices of the polyhedron are determined using spherical coordinates. Then, the distance L between the vertices is determined to see if interference occurs. This is done by comparing L with a preset interference distance threshold Lmin, and based on the comparison result, it is determined whether interference occurs between two vertices. Step A2: Distinguishing the inner and outer surfaces of the polyhedron First, let plane abc, formed by any three points in space, define the outward normal vector of this plane as... Next, iterate through all the vertices of the polyhedron, and define the direction vector connecting the vertex to any point in plane abc as follows: Subsequently based on The value determines the outer surface of the polyhedron; Step A3: Polyhedral Solid Modeling After determining the outer surface of the polyhedron, multiple data files containing three points are obtained. The output order of the three points on the outer surface is determined, and the coordinates of the three points are output in a clockwise direction. After determining the output order of the three points on the outer surface, the clockwise output data of the three points on the outer surface of the polyhedron is obtained. The sphere is cut using the three-point sectioning command in CAD. After each cut, Boolean operations are used to delete the cut outer surface and keep the inner surface. After traversing all the sectioning point sets, the required polyhedron is obtained. Step A4: Bed Surface Shear Stress Analysis Step A4.1: The bed shear stress is the shear stress exerted by the fluid on the surface of the filling bed, which is calculated based on the wall function of three-dimensional turbulence; then, based on the preset shear stress threshold of the bed shear stress, it is determined whether sediment will be deposited. Step A4.2: Determine the dimensionless form of the bed shear stress: the Shields number; Step A5: Grouting Simulation Analysis Step A5.1: The grouting section is simulated and analyzed using the critical Shields number for soil improvement. Step A5.2 Simultaneously, extract the overshear stress of the bed surface, where the overshear stress of the bed surface is expressed as the portion of the bed surface shear stress that exceeds the critical shear stress of the sediment. Step A5.3: Conduct relevant analysis on the fixed-bed scour model based on the analysis of flow around stones; Step A6: Comparison of different stone placement methods Comparative analysis of the protective effects of two stone placement methods: staggered circular arrangement and random arrangement; Step A7: Recommended Design Scheme for Joint Protection Statistical analysis was performed on the overshear stress mesh data within the computational domain.

2. The method for joint simulation modeling of scour-induced riprap grouting based on random convex polyhedral stones according to claim 1, characterized in that, In StepA1, if L > Lmin, it means that there is no interference between the two vertices; otherwise, it means that there is interference between the two vertices.

3. The method for joint simulation modeling of scour-induced rock grouting based on random convex polyhedral stones according to claim 1, characterized in that, In StepA2, its outer normal direction vector The calculation formula is as follows: in, The vertex points to the direction vector connecting to any point in plane abc. The calculation formula is as follows: In Step A2, if it exists If the value is less than 0, that is, the angle between the two vectors is greater than 90°, it means that all points in the space are not on the same side of the plane, and the plane does not meet the convexity condition and is not the outer surface of the polyhedron, so the plane is discarded. If any If the values ​​of the two vectors are all greater than 0, that is, the angle between the two vectors is less than 90°, it means that all points in the space are on the same side of the plane. The plane meets the convexity condition and is defined as the outer surface of the polyhedron.

4. The method for joint simulation modeling of scour-induced riprap grouting based on random convex polyhedral stones according to claim 1, characterized in that, In Step A3, the specific steps to determine the clockwise output order of the three points are as follows: Step A3.1: Select three points a, b, and c on the outer surface of the polyhedron and calculate... Step A3.2, Calculation and The vector product, denoted as Step A3.3, Calculation and The dot product; If their dot product is less than 0, it indicates that a, b, c are in clockwise order; Conversely, if a, b, c are in a counter-clockwise order, then output them again in the order of a, c, b.

5. The method for joint simulation modeling of scour-induced riprap grouting based on random convex polyhedral stones according to claim 1, characterized in that, In shallow water flow, the second friction law is used to calculate the shear stress in the riverbed. In the formula The depth-averaged fluid velocity; ρ f C is the fluid density. D The drag coefficient is determined by the roughness k. s Local determination: in, In the formula, h is the water depth, B = 0.71; z0 = k s / 30; In StepA4, the bed shear stress is calculated as follows: First, an intermediate speed is introduced, along with a pressure correction value for the current moment; The intermediate velocity is obtained by solving the momentum equation; Then, the relationship between the obtained intermediate velocity and the pressure correction value is substituted into the continuity equation to generate the pressure Poisson equation containing the pressure correction value. Finally, the GMRES algorithm is applied to solve the pressure Poisson equation; In Step A4, when the bed shear stress exceeds the preset shear stress threshold, the eddy current carries the particles away from the riverbed surface and suspends them. The suspended sediment particles are carried by the water flow. Sediment deposition occurs when suspended sediment particles settle onto the sand bed under the combined action of gravity, buoyancy, and friction; Among them, the shear stress threshold is expressed as the critical shear stress; The formula in StepA4.1 is as follows: In the formula, u τ Shear rate, τ is the bed shear stress; ρ is the bulk density of the water and sediment mixture; Y is the distance between the center point of the element grid and the boundary wall; v is the kinematic viscosity of the fluid; k = 0.4 is the von Kármán constant. s For roughness, k s =c s d 50 d 50 c is the median particle size of the sediment. s For a custom coefficient, c s The value is 2.5; StepA4.

2. The dimensionless form of the bed shear stress is defined as: In the formula, θ n ρ is the Shields number, g is the absolute value of gravity; n d represents the mass density of the sediment particles. n The diameter of the particle; In StepA4, fluid analysis and calculation are also performed based on Flow-3D; Flow-3D uses the continuity equation and the Navier-Stokes equation for incompressible viscous fluids as the governing equations of fluid motion.

6. The method for joint simulation modeling of scour-induced riprap grouting based on random convex polyhedral stones according to claim 1, characterized in that, In StepA5, the critical Shields number is used to determine whether sediment will start to move from a dimensionless perspective; The critical Shield number is used to characterize the critical bed shear stress τ. cr,n The dimensionless form, when the bed shear stress τ≥τ cr,n At that time, the sediment began to move; The specific calculation is as follows: Where θ cr,n Let θ be the critical Shields number. cr,n The calculation equation is determined based on the properties of the sediment and is as follows: Where, d *,n For dimensionless particle size, d *,n It reflects the influence of factors such as sediment particle diameter, gravitational acceleration, density difference between sediment and fluid, and fluid kinematic viscosity on sediment movement, and is a key intermediate parameter for calculating the critical Shields number. d *,n The formula for calculation is: Where, d n The diameter of the sediment particles is represented by s, g is the acceleration due to gravity, and s is the velocity of the sediment particles. n =ρ n / ρ f , ρ n ρ is the mass density of the sediment particles. f For fluid density, s n V represents the relative density of sediment and fluid. f It is the kinematic viscosity of a fluid, which represents the fluid's ability to resist deformation. It is related to the fluid's viscosity. The smaller the value, the easier the fluid flows, and the force exerted on the sediment particles will change accordingly, thus affecting the movement state of the sediment. In StepA5, the net exchange rate between deposited and suspended sediments is determined by the entrainment velocity and deposition velocity of sediment; Among them, the entrainment velocity and deposition velocity of sediment are regarded as two opposite micro-processes occurring simultaneously, and these two processes are interrelated and occur simultaneously; The sediment entrainment velocity is calculated according to the following formula: Where, α n Let n be the entrainment coefficient of the nth type of sediment. b The normal vector of the filling bed surface is outward; The sediment settling velocity is calculated according to the following formula: Where g is the acceleration due to gravity, and u settle,n Same direction; In StepA5, the dimensionless form of the sediment transport rate of the nth layer is defined as: Where, Φ n Let q be the sediment transport rate of the nth layer. b,n q represents the unit width volumetric sediment transport rate of the bed surface. b,n The calculation formula is as follows: q bn =B n (i n -θ c,n ) 1.5 c b,n Among them, B n The sediment transport coefficient of the bed layer; Among them, c b,n Let be the volume fraction of n substances in the bed.

7. The method for joint simulation modeling of scour-induced riprap grouting based on random convex polyhedral stones according to claim 1, characterized in that, In StepA5.1, the simulation analysis method is as follows: Extract the critical erosion shear stress of the known silt-solidified soil, whose value is between 2.5 - 4 Pa, and then calculate the equivalent critical Shields number of the grouted soil body according to the critical Shields number formula; Then, the anti-erosion effect after grouting is evaluated using the weighted shear stress, and the formula is as follows: Where A wcz The shear stress is the area-weighted value, where τ is the bed surface shear stress. cr The bed surface is subjected to boundary shear stress; In StepA5.2, the method for extracting the value of the excess shear stress per unit area is as follows: In the formula, τ is the mean value of the overshear stress. i A represents the shear stress value at a certain grid point on the bed surface; i The area is the grid area. In StepA5.3, the relevant analysis method is to extract the excess shear stress data within the 4D range in front, behind, left, and right of the pile foundation, that is, within the nested grid.

8. The method for joint simulation modeling of scour-induced rockfill grouting based on random convex polyhedral stones according to claim 1, characterized in that, For the staggered circular arrangement, the method is as follows: The stones are distributed in three layers, and the interval angle between adjacent stones in each layer is 30°, and the adjacent two layers are staggered by a rotation angle of 15°; For the random arrangement, the method is as follows: StepA6.2.

1. Randomly place polyhedrons within the set protection area, considering placing them within a circular ring area, and use polar coordinates to set the center point of placement: StepA6.2.

2. At the same time, control the distance between the centers of any two polyhedrons within the specified range through a preset lower limit value of the distance; Step A6.2.3: At the center point coordinates (xb) of the generated polyhedron i ,yb i ,zb i After that, a random polyhedron is first generated at the origin (0, 0, 0), and then the random polyhedron is distributed around the cylinder by translation, thus completing the random arrangement of the random polyhedron.

9. The method for joint simulation modeling of scour-induced riprap grouting based on random convex polyhedral stones according to claim 1, characterized in that, In StepA6.2.1, the method for setting the center point of placement using polar coordinates is as follows: In the formula, R1 and R2 are the two radius values used to define the placement area in the polar coordinate setting; among them, R1 is the inner radius and R2 is the outer radius, indicating that when the polyhedron is placed within the circular ring area, this area is the ring range determined by the radii R1 and R2; ρ is the radial variable in polar coordinates, representing the distance from the pole to a certain point; among them, the value range of ρ is R1 < ρ < R2, that is, the distance from the placement point of the polyhedron to the center of the placement area should be between R1 and R2, and the pole represents the center of the placement area; In StepA6.2.2, the distance calculation formula is as follows: In the formula L1 min This is a preset lower limit for the spacing, representing the minimum allowable distance between the center points of two polyhedra, and L1 min =3r1, where r1 is the stone particle size, (xb1, yb1, zb1) and (xb2, yb2, zb2) are the coordinates of the two polyhedra respectively, and L1 is the distance between the center points of the two polyhedra.

10. The method for joint simulation modeling of scour-induced riprap grouting based on random convex polyhedral stones according to claim 1, characterized in that, In StepA7, the statistical analysis method is as follows: The statistical data of the overshear stress cloud map were divided into 12 data intervals: [0, 0.1], [0.1, 0.2], [0.2, 0.3], ..., [1.1, 1.2]; For each data interval, calculate the corresponding frequency and probability density; the frequency calculation formula is: Where n i Let N be the number of data points in the i-th data interval, and N be the total number of statistical data points; the probability density formula is: Where L represents the length of the data interval, and L is taken as 0.1 in data analysis; In Step A7, the expected value of the overshear stress is also calculated based on the probability density graph of the overshear stress. The expected value is calculated as follows: Where E(X) is the expected value of the overshear stress, X i Let P(X) be the midpoint of the i-th data interval. i ) represents the probability corresponding to the i-th data interval, and n is the total number of data intervals.

Images