Seepage erosion simulation method based on coupling of cfd-dem and dynamic mesh

Abstract

This invention relates to a seepage erosion simulation method based on CFD-DEM coupled with a dynamic mesh. The steps include: S1. Determining initial parameters and initializing the simulation platform model; S2. The DEM module calculates the particle sample boundary deformation information for the current time step; S3. Determining whether the DEM computational domain exceeds the CFD computational domain boundary. If it does, the fluid mesh is reconstructed, and the flow field is reconstructed through the solver. The particle volume fraction within the fluid mesh is calculated on the current mesh or the reconstructed mesh, completing the CFD-DEM momentum exchange; S4. Solving the current CFD governing equations to obtain flow field information; S5. Calculating the coupling force information of the fluid acting on each particle and transferring it back to the DEM module in S2 for calculation in the next time step; S6. Repeating S2-S5 until the total simulation time is reached or the preset termination condition is met. This method provides an efficient and refined numerical simulation tool capable of completely reproducing the entire macroscopic and microscopic damage process.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation in geotechnical engineering, and in particular to a method for simulating seepage erosion based on CFD-DEM coupled with a dynamic mesh. Background Technology

[0002] Under the combined action of dynamic loads such as earthquakes, waves, or traffic, and seepage, fine particles are prone to migration and loss within the soil of major engineering projects such as earth-rock dams, levees, and building foundations. This can induce internal erosion and seepage damage, seriously threatening the safety of the engineering structure. Therefore, the demand for simulations that can reproduce the process of seepage erosion under dynamic loads, have explainable mechanisms, and be scale-compatible continues to increase.

[0003] Currently, research methods for soil seepage erosion mainly fall into two categories: numerical simulation and laboratory experiments. In numerical simulation, the use of computational fluid dynamics-discrete element method (CFD-DEM) for studying seepage erosion under static or isotropic conditions is relatively mature. However, for the more complex problem of seepage erosion under dynamic loads, related research remains primarily at the stage of laboratory physical experiments. Such experiments have inherent limitations, including high economic and time costs, difficulty in controlling experimental conditions, and a lack of ability to reveal the evolution mechanism of micro-particles.

[0004] There is currently no mature and feasible solution for applying the CFD-DEM method to dynamic load seepage erosion problems. This is mainly because existing CFD-DEM algorithms have the following insurmountable technical defects when dealing with dynamic load seepage erosion problems: When simulating large deformation processes such as triaxial compression, the fixed mesh used in the fluid part cannot adapt to particle boundary movement, which easily leads to particle boundary crossing and inconsistency between the local fluid-solid phase, resulting in decreased accuracy of flow field calculation or even calculation failure; existing algorithms themselves only support unidirectional coupling (particle → fluid or fluid → particle), ignoring the real-time feedback effect of fine particle migration on the porosity and permeability of porous media, and cannot accurately reflect the bidirectional coupling mechanism of seepage-erosion; traditional calculation methods use decoupling strategies, which cannot achieve synchronous simulation of the entire process of "external load-seepage erosion", resulting in low computational efficiency and a lack of optimized algorithms for CPU parallel computing.

[0005] Therefore, in view of the above-mentioned shortcomings of existing technologies, it is urgent to develop a new method to fill the gap in numerical simulation technology in the field of seepage erosion damage under dynamic loads. Summary of the Invention

[0006] In view of the shortcomings of the prior art, the technical problem to be solved by the present invention is to provide a seepage erosion simulation method based on CFD-DEM coupled with dynamic mesh. This method can achieve mesh adaptation as the CFD computational domain changes with the DEM computational domain, and provides a general particle calculation scheme with optimal computational efficiency through CPU parallel computing, which is applicable to various seepage erosion scenarios such as soil, rock and soil, and granular materials.

[0007] This invention is achieved through the following technical solution: a seepage erosion simulation method based on CFD-DEM coupled with dynamic mesh, comprising the following steps:

[0008] S1. Determine initial parameters and initialize the model;

[0009] S2. DEM module: Determines particle contact and overlap, solves DEM control equations, obtains motion and position information of each particle, and obtains particle sample boundary deformation information at the current time step through servo wall position;

[0010] S3. Input the particle sample boundary deformation information from S2 into the CFD solver to determine whether the DEM computational domain exceeds the CFD computational domain boundary. If the DEM computational domain exceeds the CFD computational domain, reconstruct the fluid mesh and reconstruct the flow field using the dynamicMotionSolverFvMesh solver to ensure that the particle sample is completely within the flow field range. If the DEM computational domain does not exceed the CFD computational domain, keep the current fluid mesh unchanged. Calculate the particle volume fraction within the fluid mesh on the current mesh or the reconstructed mesh to complete the CFD-DEM momentum exchange.

[0011] S4. CFD module, solves the current CFD control equations to obtain flow field information;

[0012] S5. Calculate the coupling force information of the fluid acting on each particle based on the results of S4 and transmit it back to the DEM module in S2 for calculation in the next time step;

[0013] S6. Repeat S2-S5 until the total simulation time is reached or the preset termination condition is met.

[0014] Furthermore, the specific steps of reconstructing the fluid mesh and reconstructing the flow field using the dynamicMotionSolverFvMesh solver include:

[0015] 1) The solver reads the boundary deformation data transmitted by the DEM module in S2 to obtain information on the displacement and boundary geometric changes of the particle specimen;

[0016] 2) Based on the boundary displacement information, calculate the new positions of the nodes inside the CFD mesh by solving the Laplace equation;

[0017] 3) Mesh deformation is achieved by stretching nodes, thereby maintaining the integrity of the mesh topology while ensuring smooth and uniform displacement of mesh points, maintaining the quality of the CFD mesh, and avoiding mesh distortion;

[0018] 4) Set mesh quality control indicators and automatically correct when mesh distortion exceeds the limit;

[0019] 5) To balance computational efficiency and numerical stability, the data exchange frequency between CFD and DEM is adjusted to ensure that the sample deformation is completely contained within the flow field boundary, while maintaining the maximum Coulomb number strictly less than 1.

[0020] Furthermore, in step 3), the new positions of the nodes inside the CFD mesh are calculated by solving the Laplace equation. The mathematical expression is:

[0021] (4)

[0022] Where r is the grid node position vector. The diffusion coefficient is calculated using the inverse of the node-to-boundary distance or a weight based on the element volume to reduce local distortions in high-deformation environments.

[0023] Furthermore, the modified DEM governing equations solve for the motion of each particle and calculate the contact force, thereby updating the particle position. Its translation and rotation are controlled by the following equations:

[0024] (2)

[0025] (3)

[0026] Particles quality Particles Moment of inertia, and Particles Translational velocity and angular velocity, and Particles With particles The normal and tangential contact forces between them For fluid to particles The force, The torque generated by contact.

[0027] Furthermore, the modified CFD governing equations are solved using the locally averaged Navier-Stokes equations, and their continuity and momentum equations are as follows:

[0028] (9)

[0029] (10)

[0030] In the formula, , , These represent the fluid's volume fraction, density, and velocity, respectively. For fluid pressure; This refers to fluid viscous stress; This represents the force exerted by the particles on the fluid.

[0031] Furthermore, the CFD solver in S3 calculates the volume fraction and momentum exchange of particles in the fluid mesh, and updates the porosity and permeability distribution. The specific steps are as follows:

[0032] 1) For each CFD mesh cell, particle-mesh mapping technology is used to identify the particle set within the mesh, and porosity is calculated through spatial discretization, thus achieving dynamic updating of porosity;

[0033] 2) Based on the porosity of the current time step and the porosity of the previous time step, calculate the porosity time gradient to accurately capture the spatiotemporal evolution of porosity;

[0034] 3) The modified Kozeny-Carman equation is used to update the permeability of the grid cells in real time based on the current porosity;

[0035] 4) Perform DEM to CFD feedback. Based on the updated porosity and permeability, calculate the momentum source term generated by fluid-particle interaction and introduce it into the CFD governing equation to solve the flow field, thus realizing complete DEM to CFD feedback.

[0036] Furthermore, the specific steps in S5 for calculating the fluid forces acting on each particle and transferring them to the DEM module include:

[0037] 1) For each DEM particle, interpolation is used to obtain flow field information such as fluid velocity and pressure gradient at the particle's location;

[0038] 2) Calculated fluid-structure interaction force Update total force, drag force The calculation was performed using the Di Felice model:

[0039] (11)

[0040] Among them, the drag coefficient and experience correction index It is a function of the particle Reynolds number, and its calculation formula is:

[0041] (12)

[0042] (13)

[0043] In the formula and These represent the velocities of the particles and the fluid, respectively.

[0044] 3) The calculated fluid-structure interaction force is transferred to the DEM module for particle motion calculation in the next time step.

[0045] Furthermore, during the S1 initialization process, the particle parameters need to be systematically configured in the DEM module to establish complete particle property parameters. The database includes particle property parameters and particle size distribution; flow field parameters in the CFD module; CFD-DEM coupling parameters and sample geometric parameters, as detailed below:

[0046] 1) The physical properties of the DEM particles include Young's modulus, Poisson's ratio, coefficient of restitution, and coefficient of friction;

[0047] 2) The CFD flow field parameters include fluid density, dynamic viscosity, inlet velocity or pressure boundary conditions, and outlet boundary conditions;

[0048] 3) The CFD-DEM coupling parameters include the CFD-DEM data exchange time step, the flow field reconstruction judgment threshold, and the maximum allowable Coulomb number;

[0049] 4) The geometric parameters of the sample include the initial size of the particle sample, the initial position of the particles, and the initial boundary of the CFD calculation domain.

[0050] Furthermore, the specific steps for initializing the S1 model are as follows:

[0051] 1) Within a three-dimensional calculation area of ​​a preset size, based on the particle size distribution requirements of the actual project, a discontinuous gradation numerical sample containing coarse and fine particles is randomly generated using the Monte Carlo method;

[0052] 2) The isotropic consolidation process with servo boundary control is implemented by using servo wall technology to apply target confining pressure to the numerical sample and adjusting the boundary stress in real time through a PID controller to ensure stress balance;

[0053] 3) Modify the boundary conditions of the servo wall at the top of the model to allow fine particles to pass through while maintaining impermeability to coarse particles. Under this condition, perform secondary consolidation until the confining pressure recovers to the target value.

[0054] 4) Establish a standardized fluid parameter system and construct a corresponding CFD calculation domain based on the geometry of the DEM sample to ensure that the particle motion range is fully included;

[0055] 5) Establish an initialization process for a steady seepage field by applying a pressure difference between the bottom and top boundaries of the model to simulate a bottom-up steady seepage field;

[0056] 6) Set the side boundary to a no-slip wall boundary condition to ensure the integrity of the boundary conditions for fluid calculation.

[0057] Furthermore, a simulation platform for a seepage erosion simulation method based on CFD-DEM coupled with a dynamic mesh includes a DEM module, a CFD module, and a two-way coupling module. The DEM-CFD module is coupled with the dynamic mesh to construct a simulation platform in which the CFD computational domain changes with the DEM computational domain.

[0058] Beneficial effects of this invention:

[0059] By applying a coupled algorithm of dynamic mesh and CFD-DEM to the simulation of large deformation seepage erosion processes, an adaptive fluid-solid interface was achieved. When triaxial shear causes significant deformation of the soil sample, the dynamic mesh can automatically and smoothly reconstruct the CFD flow field mesh based on the actual boundary displacement calculated by the DEM. This effectively avoids the calculation failures caused by the distortion of traditional fixed meshes, thus enabling a complete and continuous simulation of the entire "loading-failure" process.

[0060] (2) It can adapt to the fluid solution domain with large deformation and moving boundary to avoid particles crossing the boundary and maintain the consistency of fluid-solid phase. The particle-fluid bidirectional coupling is adopted to incorporate the feedback of fine particle migration on porosity, permeability and local flow field into the solution in real time. It also realizes the synchronous coupling of dynamic stress loading and seepage process to accurately reproduce the physical scene of "seepage while being stressed". Thus, it reveals the laws of particle migration, channel formation and instability evolution at the microscale, and provides reliable numerical support for engineering safety assessment and disaster prevention and mitigation.

[0061] (3) By using bidirectional coupling, the formation, blockage and evolution of local seepage channels caused by fine particle migration can be dynamically captured, thereby improving the reliability of seepage erosion prediction.

[0062] (4) It provides a new method for efficient and refined numerical simulation, supports a large-scale parallel high-efficiency simulation platform, significantly improves computational efficiency and engineering applicability, and effectively solves the pain point that the research on dynamic load seepage erosion is limited to the high-cost, macroscopic indoor test stage. It provides a powerful numerical tool for the mechanism analysis and prevention of related engineering disasters (such as earthquake liquefaction, dam instability, etc.). Attached Figure Description

[0063] Figure 1 This is a schematic diagram of a seepage erosion simulation method based on CFD-DEM coupled with dynamic mesh;

[0064] Figure 2 This is a schematic diagram visualizing the simulation results of the simulation platform of the present invention;

[0065] Figure 3 This is a schematic diagram of an indoor corrosion test according to an embodiment of the present invention;

[0066] Figure 4 This is a schematic diagram comparing the mass loss of fine particles in a collection bottle during a simulation platform and an indoor test, according to an embodiment of the present invention.

[0067] Figure 5 This is a schematic diagram comparing the stress-strain curves output by the simulation platform and the indoor test in Embodiment 2 of the present invention;

[0068] Figure 6 This is a schematic diagram comparing the amount of fine particles lost in the simulation platform and the indoor test output in Embodiment 2 of the present invention;

[0069] Figure 7 This is a schematic diagram comparing the particle distribution of the simulation platform with an initial fine particle content of 20% and the indoor test in Embodiment 2 of the present invention;

[0070] Figure 8 This is a schematic diagram comparing the particle distribution of the simulation platform with an initial fine particle content of 35% and the indoor test in Embodiment 2 of the present invention;

[0071] Figure 9 This is a schematic diagram illustrating the velocity information of particles in Embodiment 2 of the present invention. Detailed Implementation

[0072] To further illustrate the technical means and effects of the present invention in achieving its intended purpose, the following detailed description of the specific implementation methods, structures, features, and effects of the present invention, in conjunction with the accompanying drawings and preferred embodiments, is provided below.

[0073] Reference Figure 1 As shown, this invention provides a seepage erosion simulation method based on CFD-DEM coupled with a dynamic mesh. This method is based on a simulation platform where the CFD computational domain varies with the DEM computational domain. The steps include:

[0074] S1. Determine initial parameters and initialize the simulation platform model.

[0075] I. Determine the initial parameters.

[0076] In the DEM module, the system configuration of particle parameters is performed to establish a complete database of particle physical properties, including particle physical properties and particle gradation; in the CFD module, the flow field parameters are configured; and the CFD-DEM coupling parameters and sample geometric parameters are established. Details are as follows:

[0077] 1) The physical properties of the DEM particles include Young's modulus, Poisson's ratio, coefficient of restitution, and coefficient of friction;

[0078] 2) The CFD flow field parameters include fluid density, dynamic viscosity, inlet velocity or pressure boundary conditions, and outlet boundary conditions;

[0079] 3) The CFD-DEM coupling parameters include the CFD-DEM data exchange time step, the flow field reconstruction judgment threshold, and the maximum allowable Coulomb number;

[0080] 4) The geometric parameters of the sample include the initial size of the particle sample, the initial position of the particles, and the initial boundary of the CFD calculation domain.

[0081] II. Simulation platform model initialization. The specific steps are as follows:

[0082] 1) Within a three-dimensional calculation area of ​​a preset size, according to the particle size distribution requirements of the actual project, a discontinuous gradation numerical sample containing coarse and fine particles is randomly generated using the Monte Carlo method.

[0083] 2) The isotropic consolidation process with servo boundary control is implemented by using servo wall technology to apply target confining pressure to the numerical sample and adjusting the boundary stress in real time through a PID controller to ensure stress balance.

[0084] 3) Modify the boundary conditions of the servo wall at the top of the model so that it remains impermeable to coarse particles while allowing fine particles to pass through. Under this condition, perform secondary consolidation until the confining pressure returns to the target value.

[0085] 4) Establish a standardized fluid parameter system and construct a corresponding CFD calculation domain based on the geometry of the DEM sample to ensure that the particle motion range is fully included.

[0086] 5) Establish the initialization process for a steady seepage field. A bottom-up steady seepage field is simulated by applying a pressure difference between the bottom and top boundaries of the model. The pressure difference is calculated using the following formula:

[0087] (1)

[0088] In the formula, For the height of the sample, It is the density of water. It is a hydraulic gradient.

[0089] 6) Set the side boundary to a no-slip wall boundary condition to ensure the integrity of the boundary conditions for fluid calculation.

[0090] S2.DEM module.

[0091] The DEM module determines particle contact and overlap, solves the DEM control equations, obtains the motion and position information of each particle, and obtains the particle sample boundary deformation information at the current time step through the servo wall position.

[0092] In the initial state, the DEM solver extracts information from the particle property parameter database and calculates the contact force and overlap between particles based on the particle property parameters. In the next cycle, the DEM solver updates the particle position and velocity based on the fluid-structure interaction force and the contact force between particles received in the previous time step, and re-detects the particle contact and overlap.

[0093] Specifically, the DEM control equations, based on Newton's second law, solve for the motion of each particle and calculate the contact force, thereby updating the particle position. Its translation and rotation are controlled by the following equations:

[0094] (2)

[0095] (3)

[0096] The formula is Particles quality Particles Moment of inertia, and Particles Translational velocity and angular velocity, and Particles With particles The normal and tangential contact forces between them For fluid to particles The force, The torque generated by contact.

[0097] S3.DEM provides complete feedback for CFD.

[0098] S31. Input the particle sample boundary deformation information from S2 into the CFD solver and determine whether the DEM computational domain exceeds the CFD computational domain boundary. When the DEM computational domain exceeds the CFD computational domain, reconstruct the fluid mesh and reconstruct the flow field using the dynamicMotionSolverFvMesh solver to ensure that the particle sample is completely within the flow field range; when the DEM computational domain does not exceed the CFD computational domain, keep the current fluid mesh unchanged.

[0099] S32. Calculate the particle volume fraction within the fluid mesh on the current mesh or the reconstructed mesh to complete the CFD-DEM momentum exchange.

[0100] Specifically, in S31, the fluid mesh is reconstructed and the flow field is reconstructed using the dynamicMotionSolverFvMesh solver, coupling dynamic mesh technology with CFD-DEM to achieve mesh adaptation where the CFD computational domain changes with the DEM computational domain. The specific steps include:

[0101] 1) The solver reads the boundary deformation data transmitted by the DEM module in S2 to obtain information on the displacement and boundary geometric changes of the particle sample.

[0102] 2) The new positions of the nodes inside the CFD mesh are calculated by solving the Laplace equation. The mathematical expression is:

[0103] (4)

[0104] Where r is the grid node position vector. The diffusion coefficient is calculated using the inverse of the node-to-boundary distance or a weight based on the element volume to reduce local distortions in high-deformation environments.

[0105] 3) Mesh deformation is achieved by stretching nodes, thereby maintaining the integrity of the mesh topology while ensuring smooth and uniform displacement of mesh points, maintaining the quality of the CFD mesh, and avoiding mesh distortion.

[0106] 4) Set mesh quality control indicators and automatically correct when mesh distortion exceeds the limit. The mesh quality control indicators include minimum element angle, maximum element deformation rate, etc.

[0107] 5) To balance computational efficiency and numerical stability, the data exchange frequency between CFD and DEM is adjusted to ensure that the sample deformation is completely contained within the flow field boundary, while maintaining the maximum Courant number strictly less than 1.

[0108] Specifically, the CFD solver in S32 calculates the particle volume fraction and momentum exchange in the fluid mesh, and updates the porosity and permeability distribution. The specific steps are as follows:

[0109] 1) For each CFD mesh cell, particle-mesh mapping technology is used to identify the particle set within the mesh. Porosity is calculated through spatial discretization, achieving dynamic updating of porosity. The calculation formula is as follows:

[0110] (5)

[0111] in, The weighting function is Gaussian. This is the smooth length parameter.

[0112] 2) Based on the porosity of the current time step and the porosity of the previous time step, calculate the porosity time gradient to accurately capture the spatiotemporal evolution of porosity. The calculation formula is as follows:

[0113] (6)

[0114] In the formula Let be the porosity of grid cell k at the current time step n; Δt represents the porosity of grid cell k at the previous time step t-1; Δt is the time step size.

[0115] 3) The modified Kozeny-Carman equation is used to update the permeability of the grid cells in real time based on the current porosity. The calculation formula is as follows:

[0116] (7)

[0117] in, .

[0118] 4) Perform DEM-to-CFD feedback: Based on the updated porosity and permeability, calculate the momentum source term generated by fluid-particle interaction and incorporate it into the CFD governing equations to solve the flow field, achieving complete DEM-to-CFD feedback. The calculation formula is as follows:

[0119] (8)

[0120] in, This is the momentum source term within grid cell k, which directly acts on the CFD governing equations. item.

[0121] S4.CFD module.

[0122] The flow field information is obtained by solving the current CFD governing equations. These CFD governing equations are solved using the locally averaged Navier-Stokes equations, and their continuity and momentum equations are as follows:

[0123] (9)

[0124] (10)

[0125] In the formula, , , These represent the fluid's volume fraction, density, and velocity, respectively. For fluid pressure; This refers to fluid viscous stress; This is the reaction force of the particles on the fluid.

[0126] S5. Calculate the fluid forces acting on each particle and transfer them to the DEM module.

[0127] The specific steps include:

[0128] 1) For each DEM particle, interpolation is used to obtain flow field information such as fluid velocity and pressure gradient at the particle's location;

[0129] 2) Calculated fluid-structure interaction force Update total force, drag force The calculation was performed using the Di Felice model:

[0130] (11)

[0131] Among them, the drag coefficient and experience correction index It is a function of the particle Reynolds number, and its calculation formula is:

[0132] (12)

[0133] (13)

[0134] In the formula and These represent the velocities of the particles and the fluid, respectively.

[0135] 3) The calculated fluid-structure interaction force is transferred to the DEM module for particle motion calculation in the next time step.

[0136] S6. Repeat S2-S5 until the total simulation time is reached or the preset termination condition is met.

[0137] This two-way coupling process ensures complete consistency and physical realism between particle motion, pore evolution, and flow field changes. It can adapt to fluid solution domains with large deformations and moving boundaries to prevent particle transgression and maintain fluid-solid phase consistency. The particle-fluid two-way coupling allows real-time incorporation of feedback from fine particle migration on porosity, permeability, and local flow field into the solution. It also achieves synchronous coupling of dynamic stress loading and seepage processes to accurately reproduce the physical scenario of "seepage while under stress," thereby revealing the laws governing particle migration, channel formation, and instability evolution at a microscale, providing reliable numerical support for engineering safety assessment and disaster prevention and mitigation.

[0138] This invention provides a simulation platform for the seepage erosion failure process under dynamic loading, based on the aforementioned CFD-DEM coupled dynamic mesh simulation method. The platform includes a DEM module, a CFD module, and a CFD-DEM bidirectional coupling module. The CFD module is coupled with the dynamic mesh, constructing a simulation platform where the CFD computational domain varies with the DEM computational domain. The platform adopts a modular C++ architecture, integrates the open-source software OpenFOAM (CFD) and LIGGGHTS (DEM), supports multiple operating systems (Windows / Linux), and uses the OpenMP-MPI hybrid parallel framework. It possesses good scalability and engineering applicability. It achieves synchronous simulation of the entire "loading-seepage" process, establishing a multi-physics coupled computational loop. Based on a stable seepage field, the DEM triaxial shear process is synchronously initiated, with the top and bottom servo walls moving towards each other at a constant axial strain rate, and the lateral confining pressure maintained in real time by the servo walls. The simulation platform supports full-process visualization, allowing users to monitor multi-physics information such as particle motion, flow field evolution, porosity, and permeability changes in real time. Figure 2 As shown, the three-dimensional visualization interface of the simulation results intuitively displays the particle distribution, fluid streamlines, and particle loss process.

[0139] This invention employs OpenMP's shared-memory parallel computing mode, combined with MPI's distributed-memory parallel technology, to construct a hybrid parallel computing framework. Within a single node, OpenMP multi-threaded parallel processing is used for particle contact detection force calculation and CFD mesh calculation. Between multiple nodes, MPI is used for data communication and task coordination, achieving efficient solutions to large-scale problems. A load balancing strategy is employed in the parallel computing, with the particle and fluid domains partitioned spatially in a uniform manner. Structured data packets are used for inter-node communication to ensure data synchronization and efficiency.

[0140] To further clarify the technical solution and advantages of the present invention, the following, in conjunction with the accompanying drawings and two typical indoor test scenarios, details the application process and key details of the simulation method and computing platform for seepage erosion failure under dynamic load based on CFD-DEM and dynamic mesh coupling algorithm.

[0141] Example 1: Based on an indoor one-dimensional seepage erosion test, the application process of the simulation platform of this invention is described in detail below:

[0142] S1. Physical experiment.

[0143] S1-1. This embodiment is based on an indoor one-dimensional seepage erosion test. Sample 1 consists of a mixture of 60 g of sand (coarse particles) and 15 g of spherical glass beads (fine particles), placed in a transparent water tank 2. A constant water head is applied at the bottom of the tank, which is 1.4 times the critical hydraulic gradient for seepage, while the top is free-flowing. During the test, the fine particles are eroded by the seepage flow, and the outflowing liquid is collected in a collection bottle 3, dried, and weighed to obtain the amount of fine particle loss (e.g., ...). Figure 3 (As shown).

[0144] S2. Model building and parameter settings.

[0145] S2-1. Based on the dimensions of the experimental setup, establish a three-dimensional numerical model consistent with the physical experiment, and set the length, width, and height of the computational domain to be consistent with the water tank;

[0146] S2-2. Particle parameter values: Coarse particles and fine particles correspond to 60 g of sand and 15 g of glass beads, respectively. Input the actual particle size distribution, density and sphericity into the platform database.

[0147] S2-3. Numerical simulations of sand using the discrete element method include simulated particle Young's modulus of 70 GPa and density of 2650 kg / m³. 3 Poisson's ratio 0.3, coefficient of restitution 0.2, coefficient of sliding friction 0.5, coefficient of rolling friction 0.1, discrete element method simulation is used with a loading speed of 0.0002 m / s and a time step size of 1e. -7 s;

[0148] S2-4. Fluid parameters: Set the water density to 1000 kg / m³, dynamic viscosity to 0.001 Pa·s, and time step size to 1e. -5 s.

[0149] S3. Boundary and initial condition settings.

[0150] S3-1. The bottom is set as a constant head boundary, with a head value of 1.4 times the critical hydraulic gradient of the burrowing.

[0151] S3-2. The top is a free outflow boundary, allowing fine particles to escape with the water flow, while coarse particles are restricted.

[0152] S3-3. The sidewalls are non-slip solid walls to prevent water and particle leakage;

[0153] S3-4. Initially, all particles are evenly distributed, with fine particles filling the pores of coarse particles;

[0154] S4. Simulation process.

[0155] S4-1. Start the CFD-DEM-dynamic mesh coupled simulation. The platform automatically reads particle and fluid parameters and initializes the mesh.

[0156] S4-2. During the simulation, the platform tracks particle movement, fine particle loss, and changes in porosity and permeability in real time;

[0157] S4-3. The simulation termination condition reaches the preset time;

[0158] S5. Results Collection and Comparison.

[0159] S5-1. The platform outputs data such as fine particle loss, loss rate, and porosity evolution;

[0160] S5-2. Compare the mass of fine particles in the collection bottle from the indoor test (e.g.) Figure 4 As shown in the figure, the accuracy of the simulation platform was verified. To eliminate the dimensional discrepancy between the experimental and simulation results, a dimensionless parameter was used to describe the degree of erosion. The degree of erosion was determined by the cumulative mass of eroded fine particles. Compared with the initial fine particle mass The ratio determines, that is When fine particles are lost in The corresponding time is expressed as dimensionless time, defined as the ratio of simulation time to the time above. ,Right now .

[0161] Example 2: Based on the in-house corrosion test under triaxial loading conditions, the application process of the simulation platform of this invention is described in detail below:

[0162] S1. Triaxial corrosion test in the laboratory.

[0163] S1-1. A specially designed small triaxial testing device was used. The specimen diameter was 20 mm and the height was 40 mm. The top and bottom were perforated aluminum plates. A stepper motor at the top applied a constant displacement rate, and the vertical force and displacement were monitored during the loading process. A constant head system formed a seepage flow from top to bottom. The quality of the outflowing water was continuously monitored by a high-precision balance. After the test, the lost fine particles were collected. The specimens were grouped as follows: Group S1 had an initial fine particle content of 18% and a porosity of 0.32; Group S2 had a fine particle content of 24% and a porosity of 0.27.

[0164] S2. Model building and parameter settings.

[0165] S2-1. Based on the dimensions of the triaxial test apparatus, establish a three-dimensional numerical model with a diameter of 20 mm and a height of 40 mm;

[0166] S2-2. Particle parameters: coarse particles (2 mm glass beads), fine particles (0.212–0.3 mm glass beads), density 2.52 g / cm³, sphericity 0.97;

[0167] S1-3. Sample grouping: Group S1 initially had a fine particle content of 20% and a porosity of 0.32; Group S2 had a fine particle content of 24% and a porosity of 0.27.

[0168] S2-4. Fluid parameters: Set the water density to 1000 kg / m³, dynamic viscosity to 0.001 Pa·s, and time step size to 1e. -5 s.

[0169] S2. Boundary and initial condition settings.

[0170] S2-1. The bottom is set as a constant head boundary, with the head value being 0.5 times the critical hydraulic gradient of the burrowing.

[0171] S2-2. The top is a free outflow boundary, allowing fine particles to escape with the water flow, while coarse particles are restricted.

[0172] S2-3. The sidewalls are non-slip solid walls to prevent water and particle leakage;

[0173] S2-4. Initially, all particles are evenly distributed, with fine particles filling the pores of coarse particles;

[0174] S3. Simulation process.

[0175] S3-1. Start the CFD-DEM-dynamic mesh coupled simulation. The platform automatically reads particle and fluid parameters and initializes the mesh.

[0176] S3-2. The platform first performs an initial seepage stage simulation, with the hydraulic gradient gradually increasing from 0 to 0.5, and maintaining this gradient until the permeability stabilizes.

[0177] S3-3. During the shearing phase, the platform controls the top boundary to be loaded at a rate of 0.01 mm / min while maintaining a constant hydraulic gradient, and outputs data such as vertical force, displacement, permeability, and fine particle loss in real time. During the simulation, the platform outputs data such as vertical force, displacement, permeability, and fine particle loss in real time.

[0178] S3-4. The dynamic mesh module dynamically adjusts the CFD mesh based on the deformation of the particle domain to ensure the accuracy of the simulation;

[0179] S3-5. The simulation termination condition reaches the preset time;

[0180] S4. Results Collection and Comparison.

[0181] S4-1. The platform outputs porosity, permeability, and fine particle loss at each stage (e.g., ...). Figure 6 As shown), stress-strain curves (as shown) Figure 5 (as shown) etc.;

[0182] S4-2. Post-process the particle information output by the platform, compare the stratified particle distribution, and further verify the accuracy of the platform (e.g., Figure 7 and Figure 8 (as shown)

[0183] S4-3. Display particle velocity information (e.g.) Figure 9 (As shown).

[0184] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A seepage erosion simulation method based on CFD-DEM and dynamic mesh coupling, characterized by the steps of: Comprise: S1. Determine the initial parameters and model initialization; S2. DEM module, determine the particle contact and overlap, solve the DEM control equation, obtain the motion and position information of each particle, and obtain the current time step particle sample boundary deformation information through the servo wall position; S3. The particle sample boundary deformation information in S2 is input into the CFD solver, and it is judged whether the DEM calculation domain exceeds the CFD calculation domain boundary. When the DEM calculation domain exceeds the CFD calculation domain, the fluid mesh is reconstructed and the flow field is reconstructed through the dynamicMotionSolverFvMesh solver to ensure that the particle sample is completely located in the flow field range. When the DEM calculation domain does not exceed the CFD calculation domain, the current fluid mesh is kept unchanged; the particle volume fraction in the fluid mesh is calculated on the current mesh or the reconstructed mesh to complete the CFD-DEM momentum exchange; The specific steps of reconstructing the fluid mesh and reconstructing the flow field through the dynamicMotionSolverFvMesh solver include: 1) The solver reads the boundary deformation data transmitted by the DEM module in S2 to obtain the particle sample displacement and boundary geometry change information; 2) According to the boundary displacement information, the new position of the internal nodes of the CFD grid is calculated by solving the Laplace equation, and the mathematical expression is: (4) where r is the mesh node position vector, is the diffusion coefficient, which employs the inverse of the node-to-boundary distance or a cell volume-based weight to reduce high deformation local distortions; 3) The mesh deformation is realized through node stretching, so as to ensure the smoothness and uniform displacement of the mesh points while maintaining the integrity of the mesh topology, maintain the CFD mesh quality, and avoid mesh distortion; 4) Set the mesh quality control index, and automatically correct when the mesh distortion is out of limit; 5) In order to balance the calculation efficiency and numerical stability, the data exchange frequency between CFD and DEM is controlled to ensure that the sample deformation is completely contained in the flow field boundary, and at the same time maintain that the maximum Courant number is strictly less than 1; S4. CFD module, solving the current CFD control equation to obtain the flow field information; S5. According to the results of S4, the coupling force information of each particle acted by the fluid is calculated and transmitted to the DEM module in S2 for calculation of the next time step; S6. Cycle S2-S5 until the total simulation time or the preset termination condition is reached.

2. The CFD-DEM and dynamic mesh coupled seepage erosion simulation method according to claim 1, characterized in that: The modified DEM control equation is used to solve the motion of each particle and calculate the contact force, so as to update the particle position, and the translation and rotation are controlled by the following equations: （2） （3） where is the mass of the particle , is the moment of inertia of the particle , and are the translational and angular velocities of the particle , and are the normal and tangential contact forces between the particle and the particle , is the force exerted by the fluid on the particle , is the torque generated by the contact.

3. The CFD-DEM and dynamic mesh coupled seepage erosion simulation method according to claim 1 or 2, characterized in that: The modified CFD control equation is solved by the locally averaged Navier-Stokes equation, and its continuity equation and momentum equation are as follows: （9） （10） wherein , , are the volume fraction, density and velocity of the fluid, respectively; is the fluid pressure; is the fluid viscous stress; is the force of the particles on the fluid.

4. The CFD-DEM and dynamic mesh coupled seepage erosion simulation method according to claim 1, characterized in that: The CFD solver in S3 calculates the particle volume fraction and momentum exchange of the fluid mesh, updates the porosity and permeability distribution, and the specific steps are as follows: 1) For each CFD grid unit, the particle-grid mapping technology is used to identify the particle set in the grid, the porosity is calculated through spatial discretization, and the dynamic update of the porosity is realized; 2) Based on the current time step porosity and the previous time step porosity, the porosity time gradient is calculated to accurately capture the spatiotemporal evolution process of the porosity; 3) The modified Kozeny-Carman equation is used to update the permeability of the grid unit in real time based on the current porosity; 4) Feedback from DEM to CFD, based on the updated porosity and permeability, calculate the momentum source term generated by fluid-particle interaction and introduce it into the CFD governing equations for flow field solving, realize the complete feedback of DEM to CFD.

5. The CFD-DEM and dynamic mesh coupled seepage erosion simulation method according to claim 1, characterized in that: The specific steps of S5 include: 1) For each DEM particle, obtain the fluid velocity, pressure gradient and other flow field information at the particle location through interpolation; 2) calculated fluid-structure coupling forces update total forces, drag forces computed using the Di Felice model: （11） where the drag coefficient and the empirical correction index is a function of the particle Reynolds number, calculated as （12） （13） wherein and Vpand Vfare the velocities of the particles and fluid, respectively; 3) Transfer the calculated fluid-structure interaction force to the DEM module for particle motion calculation in the next time step.

6. The CFD-DEM and dynamic mesh coupled seepage erosion simulation method according to claim 1, characterized in that: The S1 initialization process needs to configure the particle parameters in the DEM module, establish complete particle physical parameters, and the database includes particle physical parameters, particle size distribution; CFD module flow field parameters; CFD-DEM coupling parameters and sample geometric parameters, as follows: 1) The DEM particle physical parameters include Young's modulus, Poisson's ratio, recovery coefficient and friction coefficient; 2) The CFD flow field parameters include fluid density, dynamic viscosity, inlet flow rate or pressure boundary condition, outlet boundary condition; 3) The CFD-DEM coupling parameters include CFD-DEM data exchange time step, flow field reconstruction judgment threshold, maximum allowed Courant number; 4) The sample geometric parameters include the initial size of the particle sample, the initial position of the particle, and the initial boundary of the CFD calculation domain.

7. The CFD-DEM coupled with dynamic mesh based seepage erosion simulation method as claimed in claim 1 or 2 or 6, wherein: The specific steps of S1 model initialization are as follows: 1) In a three-dimensional calculation region of a predetermined size, according to the particle size distribution requirements of the actual project, a discontinuous graded numerical sample containing coarse and fine particles is randomly generated by using the Monte Carlo method; 2) Implement the isotropic consolidation process of servo boundary control, apply target confining pressure to the numerical sample by using servo wall technology, and adjust the boundary stress in real time through the PID controller to ensure stress balance; 3) Modify the boundary conditions of the servo wall at the top of the model, which allows fine particles to pass through while keeping coarse particles impermeable, and perform secondary consolidation until the confining pressure returns to the target value; 4) Establish a standardized fluid parameter system, construct the corresponding CFD calculation domain according to the geometric size of the DEM sample, and ensure that it completely contains the particle movement range; 5) Establish an initialization process for a stable seepage field, simulate a stable seepage field from bottom to top by applying a pressure difference between the bottom and top boundaries of the model; 6) Set the side boundary as a no-slip wall boundary condition to ensure the integrity of the fluid calculation boundary conditions.

8. A simulation platform based on the CFD-DEM and dynamic mesh coupling seepage erosion simulation method according to claim 1, characterized in that: The simulation platform is constructed by including DEM module, CFD module and bidirectional coupling module, CFD-DEM module and dynamic mesh coupling.

Images