CN113268874B - LBM-DDA coupling-based accumulation body seepage erosion damage simulation calculation method - Google Patents

LBM-DDA coupling-based accumulation body seepage erosion damage simulation calculation method Download PDF

Info

Publication number
CN113268874B
CN113268874B CN202110578213.2A CN202110578213A CN113268874B CN 113268874 B CN113268874 B CN 113268874B CN 202110578213 A CN202110578213 A CN 202110578213A CN 113268874 B CN113268874 B CN 113268874B
Authority
CN
China
Prior art keywords
dda
lbm
block
blocks
contact
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110578213.2A
Other languages
Chinese (zh)
Other versions
CN113268874A (en
Inventor
黄达
孟秋杰
宋宜祥
岑夺丰
钟助
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hebei University of Technology
Original Assignee
Hebei University of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hebei University of Technology filed Critical Hebei University of Technology
Priority to CN202110578213.2A priority Critical patent/CN113268874B/en
Publication of CN113268874A publication Critical patent/CN113268874A/en
Application granted granted Critical
Publication of CN113268874B publication Critical patent/CN113268874B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention discloses a LBM-DDA coupling-based accumulation body seepage corrosion damage simulation calculation method, which comprises the steps of inputting DDA blocks of an LBM flow field to obtain an LBM-DDA fluid-solid coupling model; calculating the one-way action of the DDA block body on the LBM flow field; calculating LBM flow field momentum; calculating the fluid drag force of the LBM flow field on the DDA block body; contact search of DDA bulk; calculating the large displacement and the large deformation of the DDA block body; and judging whether the maximum calculation time is reached, and if so, obtaining a numerical simulation result of seepage erosion damage. The method takes the LBM method as a fluid calculation part and the DDA method as a solid power calculation part, realizes the processing and solving of the accumulation seepage corrosion damage fluid-solid coupling based on the LBM and DDA methods, and more truly simulates the accumulation seepage damage water-force coupling cooperation mechanism.

Description

LBM-DDA coupling-based accumulation body seepage erosion damage simulation calculation method
Technical Field
The invention belongs to the field of fluid-solid coupling numerical simulation of interaction of fluid and rock-soil mass, and particularly relates to a heap body seepage erosion simulation calculation method based on LBM-DDA coupling.
Background
A gravel-sandwiched land body formed by combining bedrock, an ancient landslide, a modern landslide body, a fourth-system sediment body, and the like is called a heap body. The material structure is disordered and lacks organic matters, and the material structure has the characteristics of poor sorting property, poor inter-particle binding force, strong permeability and the like, and the inner thinner soil and stone bodies are easy to erode under the seepage action of rainfall. Along with the erosion of the fine block bodies, the whole body of the accumulation body deforms and settles until the accumulation body is damaged, so that geological disasters such as landslide, debris flow and collapse are caused, and the life of human beings and the safety of buildings are seriously threatened. However, the bank is widely present. According to statistics, most of large geological disasters are caused by rainfall. Therefore, the research on seepage erosion damage of the accumulation body is of great significance for evaluating the stability of the side slope and researching the instability mechanism of the side slope.
In the process of seepage of rainfall into the soil body along the cracks of the accumulation body, the infiltration rule of the rainfall can be obviously changed. However, the crack shape of the accumulation is different, and the influence of the crack on the side slope infiltration is still unknown, so the infiltration mechanism of rainfall in the accumulation still needs to be further explored. When the erosion damage occurs, the large-size stones in the accumulation body only play a role of a skeleton, and the erosion degree of the fine soil stone blocks is the key for changing the permeability and the stability of the accumulation body, so under the action of rainfall infiltration, the erosion damage of the accumulation body needs to be researched from a more detailed view.
In order to study the migration mechanism of the fine earth and rock mass in the accumulation body under the action of hydraulic erosion and discuss the fluid dynamics interaction mechanism between the fine earth and rock mass and the water saturation pore, researchers often adopt an experiment and numerical simulation calculation method. The experimental research of seepage erosion can provide a large amount of valuable research data, but the shapes and the sizes of coarse blocks in the accumulation body are different, the pore structure of a seepage field is very complex, and the physical modeling difficulty is high. And because the physical experiment is limited by the size, the experimental result is often influenced by the size effect and can not be directly applied to the engineering practice. In contrast, the numerical simulation method has the advantages of short research period, high research efficiency, strong reproducibility and more comprehensive monitoring. The numerical method is introduced into the mechanism research of the seepage erosion damage of the accumulation body, so that the infiltration process of rainwater in the fracture can be researched more intuitively, and the pore water pressure change in the complex fracture can be researched; the inside tiny soil and stone blocks of the accumulation body can be tracked in real time, and the change of the microstructure of the accumulation body caused by the internal erosion and the change of the strength caused by the change of the microstructure of the accumulation body can be researched. Therefore, the numerical simulation for researching seepage erosion damage has great advantages.
As a complex discrete material, a stack has the characteristics of discontinuity, anisotropy, random pore structure and the like on a microscopic structure, and cannot be simply considered as a continuous medium. The DDA (Discontinuous Deformation Analysis) method based on the mechanics of the Discontinuous medium can effectively simulate the mechanical behavior of the Discontinuous medium under the action of complex load, an Analysis unit can be cut into any polygonal block by a contact boundary such as joints, cracks and the like, the rigid body translation, rotation and unit self Deformation of the unit are used as basic unknown quantities, any possible contact form in a block system can be obtained based on a set of efficient contact search method, the interaction among microscopic blocks can be truly reflected, and the DDA method has special advantages for simulating the erosion of the microscopic blocks. The LBM (Lattice Boltzmann Method) is very suitable for complex geometric structures, and can well solve the problem of difficult modeling of complex fractures.
Disclosure of Invention
Aiming at the defects of the prior art, the invention aims to provide a heap seepage erosion damage simulation calculation method based on LBM-DDA coupling.
The technical scheme for solving the technical problem is to provide a LBM-DDA coupling-based accumulation body seepage erosion damage simulation calculation method, which is characterized by comprising the following steps of:
step 1, dividing a flow field into Euler grids, and setting initial parameters in the Euler grids to obtain an LBM flow field; according to the grading of the accumulation body block to be simulated, a block representing the geometric shape of the accumulation body block is generated by programming, a generation position is randomly selected, a grid in a block boundary area is set as a Lagrange grid, contact parameters of the block are defined, and the DDA block is obtained; further obtaining an LBM-DDA fluid-solid coupling model;
step 2, adding a boundary condition equation of an interpolation format into the LBM-DDA fluid-solid coupling model, and dispersing the effect of the boundary of the DDA block on an Euler grid into the whole LBM flow field to obtain the one-way effect of the DDA block on the LBM flow field;
step 3, after the calculation of the one-way effect of the DDA block body on the LBM flow field is completed, collision and migration between Euler grids of the LBM flow field are carried out, and the calculation of the LBM flow field is completed;
step 4, according to the LBM flow field, calculating the fluid drag force of the LBM flow field acting on the DDA block body by applying a momentum exchange method based on Galileo invariance, and then applying the fluid drag force as a load to the corresponding DDA block body;
step 5, contact search of DDA blocks: judging the relative state between DDA blocks; the relative state comprises a contact state and an opening state, and the contact state comprises a locking state and a sliding state;
if the judgment result is the contact state, setting a contact sub-matrix formed by the contact spring, and calculating the cohesive force on the contact surface between the DDA blocks; if the judgment result is in an open state, the contact sub-matrix is not arranged, and the cohesive force on the contact surface does not exist between the DDA blocks; in the calculation of the cohesive force, if an oversliding state or an opening state appears among the DDA blocks, the cohesive force on the contact surface among the DDA blocks disappears, and the cohesive force does not exist among the DDA blocks in the subsequent time step;
step 6, dynamic calculation of the DDA blocks: according to the judgment result of the step 5, when the relative state among the DDA blocks is a contact state, calculating a unit stiffness matrix and an external force matrix of the DDA blocks to obtain a contact force when the DDA blocks generate relative sliding or collision; then synthesizing the contact force and the fluid drag force obtained in the step (4) to obtain resultant force and resultant moment; then, the acceleration of the DDA block is solved according to a Newton second law, and the speed and the displacement of the DDA block are further solved by time integration, so that the DDA block is continuously accumulated, and the large displacement calculation of the DDA block is realized;
step 7, firstly, judging whether the maximum calculation time is reached, if not, returning to the step 2, and if so, firstly, outputting the coordinates of the DDA block to obtain the porosity of the LBM-DDA fluid-solid coupling model after the seepage erosion damage and the mass and diameter distribution of the DDA block lost after the model is damaged by the seepage erosion; and then, obtaining the permeability coefficient of the model after the model is damaged by seepage corrosion through the LBM flow field, and further obtaining the calculation result of the LBM-DDA fluid-solid coupling model after the seepage corrosion damage.
Compared with the prior art, the invention has the beneficial effects that:
(1) in order to research the influence of factors such as pore water pressure, particle grading and the like in a accumulation body on the seepage rule, the LBM method is used as a fluid calculation part, the DDA method is used as a solid power calculation part, the processing and solving of accumulation body seepage erosion damage fluid-solid coupling based on the LBM and DDA methods are realized, and a accumulation body seepage damage water-force coupling cooperation mechanism is simulated more truly.
(2) Aiming at the problem that the conventional seepage model is difficult to accurately simulate the multi-scale pore seepage in a accumulation body and the complicated boundary conditions, the method adopts the collision-migration process of the LBM method to simulate the flow of seepage liquid in the accumulation body, and adopts the momentum exchange method based on Galileo invariance to calculate the fluid drag force.
The method solves the fluid drag force by adopting a momentum exchange method based on Galileo invariance. The method ensures that the movement speed of the DDA block boundary is considered in the process of transferring momentum from the fluid to the DDA block, and the Galileo constant in the DDA block movement process is followed.
(3) Aiming at the problems that the energy exchange process generated by the interaction of particles in the accumulation body is complex, and the dynamic information such as block motion, stress and the like is difficult to obtain through experiments, the DDA method is applied to calculating the particle migration and displacement in the accumulation body for the first time. Setting an initial contact state for the DDA block body and continuously updating the contact state to simulate the relative state of particles generated by the particles under the action of hydraulic erosion; the DDA dynamic calculation method is used for solving a unit rigidity matrix and an external force matrix of a DDA block body in the accumulation body, is used for simulating large displacement and large deformation in the particle migration process, and accords with the actual condition that soil and stone particles in the accumulation body migrate along with seepage.
(4) The LBM method is based on the molecular dynamics theory, and has the advantages of convenience in processing complex boundaries, easiness in programming and parallelism. The method can well simulate microscopic pores, and is convenient for researching the seepage mechanism and influencing factors of the accumulation body. Because a large number of gaps, pores and capillaries are formed in the accumulation body, and the internal structure of the accumulation body is extremely complex, compared with other fluid calculation schemes, simple boundary setting in the LBM method is particularly important for improving the calculation efficiency of the accumulation body seepage. Due to the introduction of the LBM method, the problem that the flow velocity and the direction of the crack are changed in a complex way, so that accurate simulation is difficult is solved. Due to the fact that the LBM method has natural parallelism, the method can remarkably improve the calculation efficiency. Various free interface methods are available for LBM, so that the phenomena of saturated-unsaturated seepage erosion damage and the like in the seepage of a accumulation body can be further researched in the future.
(5) The invention uses DDA method to simulate the displacement and interaction of the particles in the accumulation body, and solves the problem that the large displacement and large deformation are difficult to calculate during the migration of the particles. By iterating the displacement and deformation of each DDA mass, the erosive failure mechanism and ultimate failure mode of the heap under the effect of percolation can be manifested; for the large-size block body which plays a role of a framework, through iterative calculation of a certain time step, the deformation and displacement of the block body tend to be constant, and the final stress and deformation analysis result of the accumulation body can be obtained.
Drawings
FIG. 1 is a D2Q9 discrete velocity model of an embodiment of the invention;
FIG. 2 is a diagram of a numerical model according to an embodiment of the present invention;
FIG. 3 is a diagram of an LBM collision and migration process according to an embodiment of the present invention;
fig. 4 is a GME solution diagram according to an embodiment of the present invention.
Detailed Description
Specific examples of the present invention are given below. The specific examples are only intended to illustrate the invention in further detail and do not limit the scope of protection of the claims of the present application.
The invention provides a LBM-DDA coupling-based accumulation body seepage erosion damage simulation calculation method (method for short), which is characterized by comprising the following steps of:
step 1, dividing a flow field into rectangular Euler grids, and setting a distribution density function, a speed boundary condition and other initial parameters in the Euler grids according to engineering requirements to obtain an LBM flow field; programming to generate a block representing the geometric shape of the accumulation body block according to the grading of the accumulation body block to be simulated, randomly selecting a generation position, setting grids in a block boundary region as Lagrange grids, and defining all contact parameters such as a friction angle, normal rigidity, tangential rigidity, cohesive force and the like of block contact to obtain a DDA block; further obtaining an LBM-DDA fluid-solid coupling model;
preferably, in step 1, the distribution density function in the euler grid is determined by an equilibrium distribution density function, which is expressed by equation (1):
Figure BDA0003085086570000041
in the formula (1), fi eqIs an equilibrium distribution density function;
Figure BDA0003085086570000042
ρ is the density of the Euler grid, fiIs a distribution density function; omegaiIs a weight; e.g. of the typeiIs a discrete velocity vector;
Figure BDA0003085086570000043
v＝(vx,vy) Representing the inlet and outlet flow rates.
Step 2, adding a boundary condition equation of an interpolation format into the LBM-DDA fluid-solid coupling model, and dispersing the effect of the boundary of the DDA block on an Euler grid into the whole LBM flow field to obtain the one-way effect of the DDA block on the LBM flow field;
preferably, in step 2, the boundary condition equation of the interpolation format is as shown in formula (2):
Figure BDA0003085086570000044
in the formula (2), fi'and f'inv(i)Subscripts inv (i) and i indicate that their velocities are in opposite directions for the distribution density function after completion of the collision; q is the linear distance from a certain Euler grid to the DDA block boundary; δ fi' speed correction for DDA bulk surface motion:
Figure BDA0003085086570000045
in the formula (3), csIs the grid velocity uωIs the velocity at the DDA bulk boundary.
Step 3, after the calculation of the one-way effect of the DDA block body on the LBM flow field is completed, collision and migration between Euler grids of the LBM flow field are carried out, and the calculation of the momentum of the LBM flow field is completed;
preferably, in step 3, the collision and migration process is: redistributing the distribution density function in each Euler grid in each time step by collision, and then transferring the redistributed distribution density function to the nearest Euler grid along the direction of discrete speed;
preferably, in step 3, the collision and migration process can be expressed by the evolution equation of collision-migration:
fi'(x+ciΔt,t+Δt)-fi(x,t)＝Ωi (4)
in the formula (4), the evolution equation of collision-migration represents the collision and migration process of the distribution density function in each Euler grid in each time step; f. ofi(x, t) is a distribution density function at coordinate x at time t; f. ofi'(x+ciΔ t, t + Δ t) is fi(x, t) a distribution density function after collision-migration; omegaiIs a collision operator; c. CiThe direction of the discrete velocity.
Step 4, according to the momentum of the LBM flow field, calculating the reaction force of the LBM flow field to the DDA block body by applying a momentum exchange method (GME) based on Galileo invariance, namely the fluid drag force of the LBM flow field acting on the DDA block body, and then applying the fluid drag force as a load to the corresponding DDA block body;
preferably, in step 4, the calculation process of the fluid drag force is as follows:
1) based on the Galileo invariant momentum exchange method, calculating the fluid drag force component at the joint of each Euler grid-Lagrange grid (namely the contact area of the LBM flow field and the DDA block body) at the time t:
Fi(xs)＝(ci-vs)fi'(xf,t)-(cinv(i)-vs)f′inc(i)(xb,t) (5)
in the formula (5), xs，xbAnd xfCoordinates of the Lagrangian mesh, contact area and Euler mesh, v, respectively, of the DDA blockssIs the boundary velocity of the DDA mass;
2) solving the fluid drag force and moment T of the LBM flow field acting on the DDA block body:
Figure BDA0003085086570000051
in formula (6), R is the centroid of the DDA mass.
Step 5, contact search of DDA blocks: judging the relative state of DDA blocks by a molar-coulomb criterion; the relative state comprises a contact state and an opening state, and the contact state comprises a locking state and a sliding state;
if the judgment result is the contact state, setting a contact sub-matrix formed by the contact spring, and calculating the cohesive force on the contact surface between the DDA blocks; if the judgment result is in an open state, the contact sub-matrix is not arranged, and the cohesive force on the contact surface does not exist between the DDA blocks; in the calculation of the cohesive force, if an oversliding state or an opening state appears among the DDA blocks, the cohesive force on the contact surface among the DDA blocks disappears, and the cohesive force does not exist among the DDA blocks in the subsequent time step;
preferably, in step 5, the locking state is: the normal contact stress being compressive stress (k)ndn> 0) and shear force is less than maximum shear strength
Figure BDA0003085086570000052
A sliding state: the normal contact stress being compressive stress (k)ndn> 0) and the shear force is greater than the maximum shear strength
Figure BDA0003085086570000053
An opening state: the normal contact stress being tensile stress (k)ndn≤0)；knAnd ksNormal and tangential spring rates applied, respectively, c joint cohesion, l joint calculated length,
Figure BDA0003085086570000061
to joint the friction angle. dsAnd dnContact strain in the shear direction and normal direction, respectively;
preferably, in step 5, the judgment criterion of the relative state between the DDA blocks is a molar-coulomb criterion;
step 6, dynamic calculation of the DDA blocks: according to the judgment result of the step 5, when the relative state among the DDA blocks is a contact state, calculating a unit stiffness matrix and an external force matrix of the DDA blocks to obtain a contact force when the DDA blocks generate relative sliding or collision; then synthesizing the contact force and the fluid drag force obtained in the step (4) to obtain resultant force and resultant moment; then, the acceleration of the DDA block is solved according to a Newton second law, and the speed and the displacement of the DDA block are further solved by time integration, so that the DDA block is continuously accumulated, and the large displacement calculation of the DDA block is realized;
preferably, step 6 is in particular:
1) the assembly of the global stiffness matrix k and the calculation of the contact force with the DDA mass:
when the DDA blocks inside the DDA block system deform or move under the action of seepage erosion, all the DDA blocks follow the contact rule of no overlapping and no embedding; calculating the contact force between the contact surfaces of the DDA blocks by a molar-coulomb rule;
by assembling the stiffness of all DDA blocks themselves and the contact stiffness between DDA blocks, the global stiffness matrix k of a DDA block system can be expressed as:
Figure BDA0003085086570000062
in the formula (7), KnnIs a cell stiffness matrix of the DDA mass; the DDA block system is a system consisting of n DDA blocks;
2) calculating the displacement vector D of any point of the DDA block body through displacement in the x and y directions, rotation angle, normal strain and shear strainj(ii) a In the DDA block system, each polygonal DDA block has six degrees of freedom for displacement change and elastic deformation; displacement vector D of any point in DDA block jjExpressed as:
Figure BDA0003085086570000063
in the formula (8), u0And v0Respectively DDA block centroid (x)0,y0) Displacement in the x and y directions; r is0Is the rotation angle of the DDA block around the centroid; (εx εy γxy) Is the normal and shear strain of the DDA bulk j;
3) updating acceleration and velocity of DDA mass: and (3) calculating the acceleration and the speed of the DDA block based on a Hamilton principle and a minimum potential energy principle:
Figure BDA0003085086570000071
in the formula (9), the reaction mixture is,
Figure BDA0003085086570000072
and D are acceleration, speed and displacement vector of DDA block body respectively; m is a quality matrix; c is a damping matrix; since the energy consumption in the DDA block system is represented by the mutual friction between the DDA blocks, the damping effect is not considered, so C is 0;
the coordinates of all DDA blocks are correspondingly updated according to the calculated displacement and rotation angle at the end of each time step;
step 7, firstly, judging whether the maximum calculation time is reached, if not, returning to the step 2, and if so, firstly, outputting the coordinates of the DDA block to obtain the porosity of the LBM-DDA fluid-solid coupling model after the seepage erosion damage and the mass and diameter distribution of the DDA block lost after the model is damaged by the seepage erosion; and then, obtaining the permeability coefficient of the model after the model is damaged by seepage corrosion through the LBM flow field, and obtaining the calculation result of the LBM-DDA fluid-solid coupling model after the seepage corrosion damage, namely the numerical simulation result of the seepage corrosion damage.
Examples
The invention is further explained by taking a two-dimensional problem as an example and combining the attached drawings. The embodiment adopts a D2Q9 model, which is only shown in a schematic diagram, not an actual diagram, and is not to be construed as limiting the invention; to better illustrate the embodiments of the present invention, some parts of the drawings may be omitted, enlarged or reduced, and do not represent the size of an actual product; it will be understood by those skilled in the art that certain well-known structures in the drawings and descriptions thereof may be omitted.
Step 1, setting initial parameters of an Euler grid. The remaining part of the euler mesh, except for the entrance and exit boundaries, is set to 0.
In the embodiment, a D2Q9 model is adopted, and as shown in FIG. 1, the value range of i is 0-8; d2 is two-dimensional, Q9 is nine directions of 0-8; upper boundary flow velocity v after the calculation startsN(0, v), lower boundary flow velocity vSThe model left and right boundaries are closed (fig. 2).
Taking the D2Q9 model as an example:
Figure BDA0003085086570000073
Figure BDA0003085086570000074
in the formula (1), the reaction mixture is,
Figure BDA0003085086570000075
in FIG. 1, in the Euler grid points at the upper boundary of the model, the D2Q9 model has 3 distribution density functions f propagating downward4,f7,f8Will be set to:
Figure BDA0003085086570000081
Figure BDA0003085086570000082
Figure BDA0003085086570000083
Figure BDA0003085086570000084
in fig. 1, in the euler grid points located at the lower boundary of the model, the D2Q9 model of which has 3 distribution density functions f propagating upward2,f5,f6Will be set to:
Figure BDA0003085086570000085
Figure BDA0003085086570000086
Figure BDA0003085086570000087
Figure BDA0003085086570000088
step 3, as shown in fig. 3, the LBM flow field is divided into rectangular euler grids with side length Δ x, and the distribution density function fi(x, t) collide at the center of the grid (fig. 3 a-b) and then migrate to the surrounding grid (fig. 3 c).
In step 4, calculating each Euler grid-Lar at the time t based on a Galileo invariant momentum exchange methodThe fluid drag force component at the junction of the grignard meshes (i.e., the region where the LBM flow field contacts the DDA mass) is as shown in fig. 4, the adjacent fluid and solid meshes form a fluid-solid mesh junction, the distribution density function from the fluid meshes is collided by the solid meshes, the direction is turned by 180 °, and two directional distribution density functions are formed at the same junction. (c)i-vs)fi'(xfAnd t) is the distribution density function f of the fluid gridi'(xfT) an increase in momentum of the DDA mass resulting from entering the DDA mass in the direction of i of the flow fastening connection; and f'inv(i)(xbT) the amount of momentum reduction of the blocks caused by the same flow fastening connection flowing out of the DDA blocks in the opposite direction inv (i).
Nothing in this specification is said to apply to the prior art.

Claims (10)

1. A heap seepage erosion damage simulation calculation method based on LBM-DDA coupling is characterized by comprising the following steps:
step 1, dividing a flow field into Euler grids, and setting initial parameters in the Euler grids to obtain an LBM flow field; according to the grading of the accumulation body block to be simulated, a block representing the geometric shape of the accumulation body block is generated by programming, a generation position is randomly selected, a grid in a block boundary area is set as a Lagrange grid, contact parameters of the block are defined, and the DDA block is obtained; further obtaining an LBM-DDA fluid-solid coupling model;
step 2, adding a boundary condition equation of an interpolation format into the LBM-DDA fluid-solid coupling model, and dispersing the effect of the boundary of the DDA block on an Euler grid into the whole LBM flow field to obtain the one-way effect of the DDA block on the LBM flow field;
step 3, after the calculation of the one-way effect of the DDA block body on the LBM flow field is completed, collision and migration between Euler grids of the LBM flow field are carried out, and the calculation of the LBM flow field is completed;
step 4, according to the LBM flow field, calculating the fluid drag force of the LBM flow field acting on the DDA block body by applying a momentum exchange method based on Galileo invariance, and then applying the fluid drag force as a load to the corresponding DDA block body;
step 5, contact search of DDA blocks: judging the relative state between DDA blocks; the relative state comprises a contact state and an opening state, and the contact state comprises a locking state and a sliding state;
if the judgment result is the contact state, setting a contact sub-matrix formed by the contact spring, and calculating the cohesive force on the contact surface between the DDA blocks; if the judgment result is in an open state, the contact sub-matrix is not arranged, and the cohesive force on the contact surface does not exist between the DDA blocks; in the calculation of the cohesive force, if an oversliding state or an opening state appears among the DDA blocks, the cohesive force on the contact surface among the DDA blocks disappears, and the cohesive force does not exist among the DDA blocks in the subsequent time step;
step 6, dynamic calculation of the DDA blocks: according to the judgment result of the step 5, when the relative state among the DDA blocks is a contact state, calculating a unit stiffness matrix and an external force matrix of the DDA blocks to obtain a contact force when the DDA blocks generate relative sliding or collision; then synthesizing the contact force and the fluid drag force obtained in the step (4) to obtain resultant force and resultant moment; then, the acceleration of the DDA block is solved according to a Newton second law, and the speed and the displacement of the DDA block are further solved by time integration, so that the DDA block is continuously accumulated, and the large displacement calculation of the DDA block is realized;
step 7, firstly, judging whether the maximum calculation time is reached, if not, returning to the step 2, and if so, firstly, outputting the coordinates of the DDA block to obtain the porosity of the LBM-DDA fluid-solid coupling model after the seepage erosion damage and the mass and diameter distribution of the DDA block lost after the model is damaged by the seepage erosion; and then, obtaining the permeability coefficient of the model after the model is damaged by seepage corrosion through the LBM flow field, and further obtaining the calculation result of the LBM-DDA fluid-solid coupling model after the seepage corrosion damage.
2. The method for simulating and calculating the erosive damage to the heap body based on LBM-DDA coupling as claimed in claim 1, wherein in step 1, the initial parameters of Euler's grid include distribution density function and velocity boundary condition.
3. The LBM-DDA coupling-based heap seepage erosion simulation calculation method of claim 2, wherein in step 1, the distribution density function in the Euler grid is determined by an equilibrium distribution density function, and the equilibrium distribution density function is represented by formula (1):
Figure FDA0003500131720000021
in the formula (1), fi eqIs an equilibrium distribution density function;
Figure FDA0003500131720000022
ρ is the density of the Euler grid, fiIs a distribution density function; omegaiIs a weight; e.g. of the typeiIs a discrete velocity vector;
Figure FDA0003500131720000023
v＝(vx,vy) Representing the inlet and outlet flow rates.
4. The LBM-DDA coupling-based heap seepage erosion damage simulation calculation method of claim 1, wherein in the step 2, the boundary condition equation of the interpolation format is as shown in formula (2):
Figure FDA0003500131720000024
in the formula (2), fi'and f'inv(i)Subscripts inv (i) and i indicate that their velocities are in opposite directions for the distribution density function after completion of the collision; q is the linear distance from a certain Euler grid to the DDA block boundary; δ fi' speed correction for DDA bulk surface motion:
Figure FDA0003500131720000025
in the formula (3), csIs the grid velocity uωIs the velocity at the DDA bulk boundary.
5. The LBM-DDA coupling-based heap seepage erosion simulation calculation method according to claim 1, wherein in step 3, the collision and migration process is: and redistributing the distribution density function in each Euler grid in each time step by collision, and then transferring the redistributed distribution density function to the nearest Euler grid along the direction of the discrete speed.
6. The method for simulating and calculating the seepage erosion of the heap body based on LBM-DDA coupling as claimed in claim 1, wherein in the step 3, the collision and migration process can be expressed by the evolution equation of collision-migration:
fi'(x+ciΔt,t+Δt)-fi(x,t)＝Ωi (4)
in the formula (4), fi(x, t) is a distribution density function at coordinate x at time t; f. ofi'(x+ciΔ t, t + Δ t) is fi(x, t) a distribution density function after collision-migration; omegaiIs a collision operator; c. CiThe direction of the discrete velocity.
7. The LBM-DDA coupling-based method for simulating calculation of corrosion damage to the seepage of a heap body according to claim 1, wherein in the step 4, the calculation process of the fluid drag force is as follows:
1) calculating the fluid drag force component of each Euler grid-Lagrange grid connection at the time t based on a Galileo invariant momentum exchange method:
Fi(xs)＝(ci-vs)fi'(xf,t)-(cinv(i)-vs)f′inc(i)(xb,t) (5)
in the formula (5), xs，xbAnd xfCoordinates of the Lagrange grid, the contact area and the Euler grid, respectively, of the DDA blocks，vsIs the boundary velocity of the DDA mass;
2) solving the fluid drag force and moment T of the LBM flow field acting on the DDA block body:
Figure FDA0003500131720000031
in formula (6), R is the centroid of the DDA mass.
8. The LBM-DDA coupling-based heap seepage erosion simulation calculation method of claim 1, wherein in step 5, the locking state is: the normal contact stress being compressive stress, i.e. kndn> 0 and the shear force is less than the maximum shear strength, i.e.
Figure FDA0003500131720000032
A sliding state: the normal contact stress being compressive stress, i.e. kndn> 0 and the shear force is greater than the maximum shear strength, i.e.
Figure FDA0003500131720000033
An opening state: the normal contact stress being tensile stress, kndn≤0；knAnd ksApplied normal and tangential spring rates, respectively; c is joint cohesion, l is joint calculated length,
Figure FDA0003500131720000034
a joint friction angle; dsAnd dnThe strain of the contact in the shear direction and the normal direction, respectively.
9. The LBM-DDA coupling-based stack seepage erosion simulation calculation method according to claim 1, wherein in the step 5, the relative state judgment criterion between DDA blocks is a Moore-Coulomb criterion.
10. The LBM-DDA coupling-based heap seepage erosion damage simulation calculation method according to claim 1, wherein the step 6 specifically comprises:
1) the assembly of the global stiffness matrix k and the calculation of the contact force with the DDA mass:
when the DDA blocks inside the DDA block system deform or move under the action of seepage erosion, all the DDA blocks follow the contact rule of no overlapping and no embedding; calculating the contact force between the contact surfaces of the DDA blocks by a molar-coulomb rule;
by assembling the stiffness of all DDA blocks themselves and the contact stiffness between DDA blocks, the global stiffness matrix k of a DDA block system can be expressed as:
Figure FDA0003500131720000035
in the formula (7), KnnIs a cell stiffness matrix of the DDA mass; the DDA block system is a system consisting of n DDA blocks;
2) calculating the displacement vector D of any point of the DDA block body through displacement in the x and y directions, rotation angle, normal strain and shear strainj(ii) a In the DDA block system, each polygonal DDA block has six degrees of freedom for displacement change and elastic deformation; displacement vector D of any point in DDA block jjExpressed as:
Figure FDA0003500131720000041
in the formula (8), u0And v0Respectively DDA block centroid (x)0,y0) Displacement in the x and y directions; r is0Is the rotation angle of the DDA block around the centroid; (εx εy γxy) Is the normal and shear strain of the DDA bulk j;
3) updating acceleration and velocity of DDA mass: and (3) calculating the acceleration and the speed of the DDA block based on a Hamilton principle and a minimum potential energy principle:
Figure FDA0003500131720000042
in the formula (9), the reaction mixture is,
Figure FDA0003500131720000043
and D are acceleration, speed and displacement vector of DDA block body respectively; m is a quality matrix; c is a damping matrix; since the energy consumption in the DDA block system is represented by the mutual friction between the DDA blocks, the damping effect is not considered, so C is 0;
the coordinates of all DDA blocks are updated accordingly at the end of each time step, based on the calculated displacement and rotation angles.
CN202110578213.2A 2021-05-26 2021-05-26 LBM-DDA coupling-based accumulation body seepage erosion damage simulation calculation method Active CN113268874B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110578213.2A CN113268874B (en) 2021-05-26 2021-05-26 LBM-DDA coupling-based accumulation body seepage erosion damage simulation calculation method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110578213.2A CN113268874B (en) 2021-05-26 2021-05-26 LBM-DDA coupling-based accumulation body seepage erosion damage simulation calculation method

Publications (2)

Publication Number Publication Date
CN113268874A CN113268874A (en) 2021-08-17
CN113268874B true CN113268874B (en) 2022-03-22

Family

ID=77232970

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110578213.2A Active CN113268874B (en) 2021-05-26 2021-05-26 LBM-DDA coupling-based accumulation body seepage erosion damage simulation calculation method

Country Status (1)

Country Link
CN (1) CN113268874B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114239352B (en) * 2021-12-14 2024-04-30 西南交通大学 Fluid-solid coupling method for depth integral fluid model and block system

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110261573A (en) * 2019-05-16 2019-09-20 同济大学 A kind of high position rock landslip stability dynamic value evaluation method
CN111368487A (en) * 2020-03-17 2020-07-03 广西师范大学 Flow field processing method for simulating periodic movement of particles based on lattice Boltzmann model

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110261573A (en) * 2019-05-16 2019-09-20 同济大学 A kind of high position rock landslip stability dynamic value evaluation method
CN111368487A (en) * 2020-03-17 2020-07-03 广西师范大学 Flow field processing method for simulating periodic movement of particles based on lattice Boltzmann model

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Dynamic process analysis of the Niumiangou landslide triggered by the Wenchuan earthquake using the finite difference method and a modified discontinuous deformation analysis;Huang D;《 Journal of Mountain Science》;20210412;全文 *
基于LBM-DEM法的二元颗粒体侵蚀破坏机制研究;马祺瑞等;《中国力学大会（CCTAM 2019）》;20190825;全文 *

Also Published As

Publication number Publication date
CN113268874A (en) 2021-08-17

Similar Documents

Publication Publication Date Title
Wu et al. Post-failure simulations of a large slope failure using 3DEC: The Hsien-du-shan slope
Ye et al. Seismic dynamics of offshore breakwater on liquefiable seabed foundation
Wang et al. Numerical modelling of fluid-induced soil erosion in granular filters using a coupled bonded particle lattice Boltzmann method
Li et al. Experimental and numerical investigations on the shear behavior of a jointed rock mass
Chen et al. Global concurrent cross-scale nonlinear analysis approach of complex CFRD systems considering dynamic impervious panel-rockfill material-foundation interactions
Qiu et al. Numerical simulation of submarine landslide tsunamis using particle based methods
Abe et al. Numerical simulation for runout process of debris flow using depth-averaged material point method
CN113268874B (en) LBM-DDA coupling-based accumulation body seepage erosion damage simulation calculation method
Lueprasert et al. Three dimensional finite element analysis for preliminary establishment of tunnel influence zone subject to pile loading.
Lai et al. Stress analysis of CFG pile composite foundation in consolidating saturated mine tailings dam
Issa et al. Numerical analysis of micromechanical behaviour of granular materials
Di et al. An operator‐split ALE model for large deformation analysis of geomaterials
Svalova Mechanical-mathematical modeling and monitoring for landslide processes
Hu et al. Numerical simulation of the entrainment effect during mass movement in high-speed debris avalanches
Wu Applying discontinuous deformation analysis to assess the constrained area of the unstable Chiu‐fen‐erh‐shan landslide slope
Peng et al. Dynamic simulation of the water inrush process in tunnel construction using a three-dimensional coupled discontinuous deformation analysis and smoothed particle hydrodynamics method
Zhang Challenges of high dam construction to computational mechanics
CN112241603B (en) Numerical simulation method for high-order landslide impact scraping and underlayer converging process
CN111597712A (en) Calculation method for movement of vehicle in discrete meta-model
CN116562189A (en) Optimization method, system and storage medium of plugging particle material for dynamic crack leakage
Cheng et al. Calculation models and stability of composite foundation treated with compaction piles
Li et al. Particle swarm optimization algorithm coupled with finite element limit equilibrium method for geotechnical practices
Zhao et al. Advances in discontinuous numerical methods and applications in geomechanics and geoengineering
Xiong et al. Arching development above active trapdoor: insight from multi-scale analysis using FEM–SPH
Yuan et al. Three-dimensional numerical modeling study on the failure evolution of tailing dam

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant