A basic scouring fluid-solid-soil coupling simulation method based on SPH-DEM coupling and multiphase flow theory

Abstract

The application discloses a foundation scouring flow-solid-soil coupling simulation method based on SPH-DEM coupling and multiphase flow theory, which comprises the following steps: constructing a particle model, setting fluid particles and bed sand particles based on the multiphase flow theory, and setting rigid particles based on the DEM theory; correcting the fluid control equation based on the flow-solid coupling theory, and performing control solving of the fluid particles; performing control solving of the DEM rigid particles based on Newton's second law and introducing flow-solid coupling force; introducing flow-soil and soil-solid coupling force to correct the sediment incipient model based on the Shield criterion of sediment particles, and performing control solving of the foundation scouring; completing one time step and entering the next cycle; through secondary development based on the SPH-DEM coupling and the multiphase flow theory, the application introduces various structure and state models, and can realize fine numerical simulation of the foundation scouring flow-solid-soil coupling.

Description

TECHNICAL FIELD

[0001] The application relates to a scour simulation method based on SPH-DEM coupling and multiphase flow theory and considering fluid-structure coupling, in particular to a secondary development based on an SPH multiphase flow algorithm, introduction of a DEM calculation theory and a plurality of experimental criteria for real-time fine simulation of a whole process of fluid-structure-soil coupling under scour, and belongs to the field of numerical simulation calculation. BACKGROUND

[0002] Bridge foundation scour disease is one of the main reasons for the failure of the function of the bridge structure and the loss of safety performance, which has attracted widespread attention of many scholars. Numerical simulation is still one of the most effective methods for studying bridge scour, and has many advantages such as low cost, high efficiency and short cycle. The mainstream algorithm used at present is to obtain fluid dynamics based on the Euler form solution of N-S fluid control equation, and to obtain the numerical solution of bed elevation change by introducing a numerical bed sand transport model. However, the numerical simulation method based on the Euler form is often difficult to converge when solving the breaking of complex fluid surface and the overtopping of waves and tsunamis, or requires a large amount of solving time and resources, and there is no real soil model, so it is difficult to track the evolution trajectory of the soil. The existing numerical model basically treats the scour foundation as a fixed boundary, and there is almost no real-time analysis of the influence of scour excavation on the stability of the foundation, and there is no case considering the influence of foundation extrusion pressure on the initiation of sediment, and there is a lack of effective tools for simulating the multi-field coupling effect of fluid-structure-soil under the background of foundation scour.

[0003] The smooth particle hydrodynamics (SPH) method is a numerical solution method based on the Lagrangian form. Compared with the Euler method, the particles themselves have mass, and mass conservation can be ensured without additional calculation. When simulating complex free surface flow, it is not necessary to track the fluid boundary and different fluid interfaces. When simulating scour problems, a soil model can be constructed without introducing a numerical sand transport model, which can further improve the efficiency and accuracy. The discrete element (DEM) algorithm also belongs to a numerical solution method based on the Lagrangian form, and is often used to simulate strong nonlinear mechanical behaviors such as large deformation and fracture of structures, and has good compatibility with SPH. However, at present, there are few cases of scour simulation based on SPH-DEM coupling, and there is a lack of design schemes for related algorithms, which need to be further solved. SUMMARY

[0004] The technical problem to be solved by the application is to introduce a plurality of theoretical or experimental criteria for fine simulation of the whole process of fluid-structure-soil coupling under the background of foundation scour, based on SPH-DEM coupling and multiphase flow theory and considering fluid-structure coupling calculation, so as to realize real-time analysis of the dynamic response of the foundation under scour excavation, and consider the influence of foundation extrusion pressure on the initiation of sediment, and have the characteristics of easy programming, high accuracy and strong operability.

[0005] The application adopts the following technical solutions to solve the above technical problems:

[0006] A foundation scouring fluid-solid-soil coupling simulation method based on SPH-DEM coupling and multiphase flow theory, comprising the following steps:

[0007] Step 1, constructing a particle model, setting fluid particles and bed sand particles based on the multiphase flow theory, and setting rigid particles based on the DEM theory; the specific steps are as follows:

[0008] Step 1.1, based on the basic principle of multiphase flow, respectively constructing a fluid particle model and a soil particle model, setting the fluid particle model as Newtonian fluid and the soil particle model as non-Newtonian fluid, determining the soil viscosity μ based on the HBP model HBP ;

[0009]

[0010] Where τ c is the yield stress of the model material, II D is the second invariant of the fluid strain rate tensor, m is the stress exponent growth coefficient, μ is the water viscosity, and n is the power related to the shear stress

[0011] Step 1.2, introducing the DP yield criterion to calculate the specific material yield stress τ y :

[0012] Step 1.3, substituting the specific material yield stress τ y obtained in step 1.2 into the soil viscosity μ HBP calculation formula to replace the model material yield stress τ c , and establish the soil viscosity calculation model

[0013] Step 1.4, setting the rigid particles based on the DEM theory, and establishing the normal contact stiffness K n , normal contact damping γ n , tangential contact stiffness K t , and tangential contact damping γ t between the rigid particles

[0014] Step 1.5, based on the normal contact stiffness K n , normal contact damping γ n , tangential contact stiffness K t , and tangential contact damping γ t between the rigid particles obtained in step 1.4, calculating the normal contact force F n and the tangential contact force F t ;

[0015] Step 2, based on the particle model obtained in step 1, introduce the fluid- solid coupling theory to modify the fluid control equation, and perform the control solving of the fluid particles, the specific steps are as follows:

[0016] Step 2.1, introduce the particle fluid-solid coupling force F fs , modify the fluid momentum conservation equation;

[0017] Step 2.2, based on the kernel function theory of SPH algorithm, give the discrete form of the fluid momentum conservation equation obtained in step 2.1:

[0018] Step 3, based on the particle model obtained in step 1, combined with Newton's second law, considering the particle fluid-solid coupling force F fs , modify the rigid body particle control equation and solve:

[0019] Step 4, based on the particle model obtained in step 1, introduce the Shield criterion of sediment particles, consider the fluid-soil coupling force obtained in step 1 and the solid-soil coupling force obtained in step 2 to modify the sediment incipient model, and perform the control solving of the foundation scour; The specific steps are as follows:

[0020] Step 4.1, based on the fluid particle velocity u obtained in step 2 momentum equation, introduce Einstein logarithmic flow velocity distribution formula, calculate the fluid shear stress τ b acted on the soil particles:

[0021] Step 4.2, based on the soil particle viscosity μ HBP model obtained in step 1.1, introduce the Shield criterion of sediment particles, calculate the sediment incipient critical stress τ cr,0 ;

[0022] Step 4.3, based on the sediment incipient critical stress τ cr,0 obtained in step 4.2, introduce the fluid-solid coupling force F fs of the sediment particles and the rigid body particles obtained in step 2, calculate the modified sediment incipient critical stress τ vr considering the external load and the effect of slope:

[0023] Step 4.4, based on the fluid shear stress τ b obtained in step 4.1 and the modified sediment incipient critical stress τ cr obtained in step 4.3, judge the starting state of the sediment particles, if τ cr ≥ τ b , then τ cr is substituted into the formula in step 1.1 to replace the yield stress τ cr , update the bed sediment particle viscosity μ HBP as the moving load particle, and if the condition is not met, it is processed according to steps 1.2 to 1.3;

[0024] ​Step 4.5, according to the obtained pusher plasmid of step 4.4, the Mastbergen formula is introduced to calculate the critical flow velocity u of the pusher phase transformation suspension phase lift , if the fluid particle flow velocity u satisfies u lift , the suspension phase is transformed into the particle, and the equivalent viscosity μ lift is calculated to replace the bed sand particle viscosity μ HBP , if the condition is not satisfied, it is not processed;

[0025] Step 4.6, according to the obtained suspension phase particle of step 4.5, the Mastbergen formula is introduced to calculate the critical flow velocity u of the suspension phase transformation pusher phase set , if the actual particle flow velocity u satisfies u set , the suspension phase is transformed into the pusher phase, and the step 4.4 is processed, if the condition is not satisfied, it is not processed;

[0026] Step 5, repeat steps 1 to step 4 until the solution is completed.

[0027] In the step 1.2, the material yield stress τ y is calculated as follows:

[0028] |τ y |=αp+β

[0029] Wherein, p is the hydrostatic pressure acting on the saturated sediment particles, and α and β are given by the Mohr-Coulomb yield criterion parameters, and the calculation formula is as follows:

[0030]

[0031] Where: θ is the internal friction angle, and c is the soil cohesion.

[0032] In the step 1.4, the normal contact stiffness K n , the normal contact damping γ n , the tangential contact stiffness K t and the tangential contact damping γ t between the rigid particles are calculated as follows:

[0033]

[0034] Wherein, C n =1*10 -5 ,K n,ij is the normal contact stiffness between the DEM rigid particle i and the rigid particle j, K t,ij , γ n,ij , γ t,ij are the same; E * is the equivalent elastic modulus, R * is the equivalent particle radius, and M *For equivalent mass, the calculation is as follows respectively:

[0035]

[0036] Wherein, E i , E j Respectively, the elastic modulus of the material endowed by DEM rigid particle i, rigid particle j;

[0037] μ i , μ j Respectively, the Poisson ratio of the material endowed by DEM rigid particle i, rigid particle j;

[0038] r i , r j Respectively, the radius of DEM rigid particle i, rigid particle j;

[0039] m i , m j Respectively, the mass of DEM rigid particle i, rigid particle j.

[0040] In step 2.1, the momentum conservation equation of the fluid is modified as:

[0041]

[0042] In the formula, u is the fluid velocity vector, P is the fluid pressure term, v is the fluid viscosity, g is the fluid gravity term, and the rest of the symbols are the same as above.

[0043] In step 2.2, the discrete form of the control equation obtained in step 2.1 is:

[0044]

[0045] Wherein, Is the Hamiltonian operator, subscript a represents the center particle, subscript b represents the neighborhood particle, W ab Is the kernel function with a particle as the search center r a And r b Respectively represent the position vector of the center particle and the neighborhood particle, q is the ratio of the particle spacing to the smoothing length, h is the smoothing length;

[0046] μ0 is the kinematic viscosity;

[0047] τ ij Is the fluid SPS stress vector;

[0048] m f , P f Respectively, the mass and pressure of the fluid particle searched by the rigid particle s within the kernel function radius;

[0049] ms , P s is the equivalent mass and pressure of rigid particle,

[0050] The remaining variables are the same as above.

[0051] In step 3, the rigid particle control equation is modified as:

[0052]

[0053] Where, m a , v a are the mass and velocity vector of rigid particle a, G is the gravity vector of rigid particle a, and the remaining symbols are the same as above.

[0054] In step 4.1, the fluid shear stress τ b is calculated as:

[0055]

[0056] Where, d is the characteristic particle size of the particle; Δu is the velocity difference between the sediment particle and the nearby water particle, κ is the von Karman constant, taken as 0.41, and ρ is the density of the sediment particle.

[0057] In step 4.2, the sediment incipient critical stress τ cr,0 is calculated as:

[0058] τ cr,0 = θ cr ·(ρ s -ρ)gd

[0059] Where, θ cr is the critical Hillz number, only related to sediment parameters, ρ s is the saturated sediment density, ρ is the water density, g is the acceleration of gravity, and d is the particle size of the soil body.

[0060] In step 4.3, the modified sediment incipient critical stress τ cr is calculated as:

[0061]

[0062]

[0063] Where, P is the modulus of F fs , e is the unit vector in the direction of flow velocity, W is the gravity experienced by the sediment particle, α, β, and γ are the angles between the slope normal vector and the X, Y, Z axis directions, w i is the components of in the X, Y, Z axis directions, is the unit vector of gravity in the slope direction, is the internal friction angle of the sediment particle.

[0064] In step 4.5, the bed load transformation critical flow velocity u lift The calculation is as follows:

[0065]

[0066] Where, alpha i is the sediment transport coefficient, n s is the bed normal vector, d * is the sediment particle size coefficient, rho s is the saturated sediment density, rho is the water density, g is the acceleration of gravity, d is the soil particle size, theta b is the actual hillz number of the soil particle, theta cr is the critical hillz number.

[0067] The equivalent viscosity mu lift The calculation is as follows:

[0068]

[0069] Where, mu is the viscosity of the water body, C v is the concentration of the sediment particles in the core function radius;

[0070] In step 4.6, the critical flow velocity u set The calculation formula is as follows:

[0071]

[0072] Where, mu is the viscosity of the water body, d is the soil particle size, d * is the sediment particle size coefficient.

[0073] The application has the following technical effects compared with the prior art:

[0074] 1. The application is based on SPH-DEM coupling and multiphase flow theory, considering flow-solid coupling calculation, introducing a plurality of theoretical or experimental criteria for real-time fine simulation of the whole process of flow-solid-soil coupling under foundation scour background, so as to realize real-time analysis of foundation dynamic response under scour, and considering the influence of foundation extrusion force on sediment starting, having the characteristics of easy programming implementation, high accuracy, strong operability and the like.

[0075] 2. The application realizes scour numerical simulation based on SPH multiphase flow secondary development, the particles themselves have mass, without additional calculation to ensure mass conservation; when simulating scour problems, the soil model can be constructed, without introducing numerical sediment transport model, which can further speed up the efficiency and improve the accuracy.

[0076] 3. The application truly realizes real-time simulation of local scour calculation and basic dynamic and stability analysis. BRIEF DESCRIPTION OF DRAWINGS

[0077] Figure 1 is a calculation flowchart of a basic scour fluid-solid-soil coupling simulation method based on SPH-DEM coupling and multiphase flow theory of the application;

[0078] Figure 2 is a fluid particle distribution and retrieval mechanism diagram simulated by the application;

[0079] Figure 3 is a DEM rigid particle contact mechanism diagram simulated by the application;

[0080] Figure 4 is a sediment type classification diagram of the scour calculation simulated by the application;

[0081] Figure 5 is an enlarged diagram of the interface between bed load and suspended load.

[0082] 1 is a fluid particle;

[0083] 2 is a center fluid particle to be searched;

[0084] 3 is a kernel function radius h;

[0085] 4 is the center of mass of a DEM rigid particle i;

[0086] 5 is the center of mass of a DEM rigid particle j;

[0087] 6 is a normal spring in the contact part of particle i and particle j, and its stiffness is K n , and the radial contact length is the radial length of the overlapping part;

[0088] 7 is the normal contact damping γ n of the contact part of particle i and particle j;

[0089] 8 is a tangential spring in the contact part of particle i and particle j, and its stiffness is K t , and the tangential contact length is the tangential length of the overlapping part;

[0090] 9 is the tangential contact damping γ t of the contact part of particle i and particle j;

[0091] 10 is ordinary soil;

[0092] 11 is yield soil;

[0093] 12 is bed load;

[0094] 13 is suspended load. DETAILED DESCRIPTION

[0095] Embodiments of the present application are described in detail below with reference to the attached drawings. The embodiments described below are examples of implementations of the present application and are not intended to limit the scope of the present application, as described in the appended claims.

[0096] The present application is a secondary development based on SPH multiphase flow algorithm and DEM calculation theory. Therefore, the default premise of the present application is that the SPH multiphase flow algorithm and the DEM calculation theory are known or open source, and the solution of partial differential equations such as fluid control equations and DEM particle motion control equations is not within the scope of discussion of the present application, and relevant variables can be artificially extracted or obtained.

[0097] The fluid particle distribution and retrieval mechanism based on the present application is shown in Figure 2 , 1 is a fluid particle, 2 is a center fluid particle to be searched, and 3 is a kernel function radius h.

[0098] As shown in Figure 3 , the DEM rigid particle contact model based on the present application, 4 is the center of mass of DEM rigid particle i, 5 is the center of mass of DEM rigid particle j, 6 is the normal spring in the contact part of particle i and particle j, and the stiffness is K n , the radial contact length is the radial length of the overlapping part, 7 is the normal contact damping γ n in the contact part of particle i and particle j, 8 is the tangential spring in the contact part of particle i and particle j, and the stiffness is K t , the tangential contact length is the tangential length of the overlapping part, and 9 is the tangential contact damping γ t .

[0099] As shown in Figure 4 , the distribution of the sediment based on the present application, the riverbed part is 10 ordinary soil, the scour pit below is 11 yield soil, the distribution in the scour pit is 12 bed load and 13 suspended load, 12 bed load is adjacent to 13 yield soil, and 13 suspended load is distributed above 12 bed load.

[0100] Figure 5 is an enlarged view of the interface between the bed load and the suspended load.

[0101] The specific implementation method of the present application takes the SPH open source calculation software Dualsphysics as an example, and the solution of partial differential equations such as fluid control equations and DEM motion control equations is handled by the corresponding algorithm module, and the equation variables can be artificially extracted and obtained.

[0102] In combination with Figure 1The specific calculation process of the basic scouring fluid-solid-soil coupling simulation method based on the SPH-DEM coupling and the multiphase flow theory is as follows:

[0103] The pretreatment module:

[0104] Step 1, construct a particle model, set fluid particles and bed sand particles based on the multiphase flow theory, and set rigid particles based on the DEM theory; the specific steps are as follows:

[0105] Step 1.1, based on the basic principle of multiphase flow, construct a fluid particle model and a soil particle model, respectively, the fluid particle model is set as Newtonian fluid, and the soil particle model is set as non-Newtonian fluid, based on the HBP model, the soil viscosity μ is determined HBP :

[0106]

[0107] Where, τ c is the yield stress of the model material, II D is the second invariant of the fluid strain rate tensor, m is the stress exponent growth coefficient, μ is the viscosity of water, and n is the power related to the shear stress;

[0108] Step 1.2, introduce the DP yield criterion to calculate the specific material yield stress τ y :

[0109] |τ y |=αp+β

[0110] Where, p is the static water pressure acting on the saturated sand particles, and α and β are given by the Mohr-Coulomb yield criterion parameters, and are calculated as follows:

[0111]

[0112] Where: θ is the internal friction angle, and c is the soil cohesion;

[0113] Step 1.3, substitute the specific material yield stress τ y obtained in step 1.2 into the soil viscosity μ HBP calculation formula to replace the model material yield stress τ c , and establish the soil viscosity calculation model;

[0114] Step 1.4, set the rigid particles based on the DEM theory, and establish the normal contact stiffness K n , the normal contact damping γ n , the tangential contact stiffness K t and the tangential contact damping γ t :

[0115]

[0116] Among them, C n =1*10 -5 ,K n,ij Let K be the normal contact stiffness between rigid particle i and rigid particle j in the DEM. t,ij γ n,ij γ t,ij Synonyms; E * For the equivalent elastic modulus, R * M is the equivalent particle radius. * For equivalent mass, the calculations are as follows:

[0117]

[0118] Among them, E i E j These are the elastic moduli of the material imparted by rigid particles i and j in the DEM, respectively, μ. i μ j The Poisson's ratios, r and r, are the ratios assigned to the material by rigid particles i and j in the DEM, respectively. i r j The radii, m, of rigid particles i and j in the DEM are respectively. i m j Let i and j be the masses of rigid body particles i and j in DEM, respectively.

[0119] Step 1.5: Determine the normal contact force F based on the contact stiffness between rigid particles. n and tangential contact force F t :

[0120]

[0121] Where, δ ij , The radial and tangential contact distances between rigid particles i and j in the DEM; e ij Let be the unit vector pointing from the centroid of rigid particle i in the DEM to the centroid of rigid particle j. These represent the radial and tangential deformation rates of rigid particles i and j in the DEM, respectively. Let u be the radial and tangential relative velocities of rigid particles i and j in the DEM. f Let be the friction coefficient between rigid body particle i and rigid body particle j in the DEM.

[0122] Fluid Solver Module:

[0123] Step 2, based on the particle model obtained in step 1, introduce the fluid-solid coupling theory to modify the fluid control equation, and perform the control solving of the fluid particles; the specific steps are as follows:

[0124] Step 2.1, introduce the particle fluid-solid coupling force F fs , the momentum conservation equation of the fluid is modified as:

[0125]

[0126] In the formula, u is the fluid velocity vector, P is the fluid pressure term, υ is the fluid viscosity, g is the fluid gravity term, and the remaining symbols are the same as above;

[0127] Step 2.2, based on the kernel function theory of SPH algorithm, give the discrete form of the control equation obtained in step 2.1:

[0128]

[0129] Wherein, is the Hamiltonian operator, subscript a represents the center particle, subscript b represents the neighborhood particle, W ab is the kernel function with a particle as the search center r a and r b respectively represent the position vectors of the center particle and the neighborhood particle, q is the ratio of the particle spacing to the smoothing length, h is the smoothing length; μ0 is the kinematic viscosity, τ ij is the fluid SPS stress vector, m f , P f are the mass and pressure of the fluid particles searched by the rigid particle s within the kernel function radius, m s , P s are the equivalent mass and pressure of the rigid particle,

[0130] DEM control solving module:

[0131] Step 3, based on the particle model obtained in step 1, combine Newton's second law and introduce the particle fluid-solid coupling force F fs , modify the rigid particle control equation and perform solving:

[0132]

[0133] Wherein, m a , v a are the mass and velocity vector of the rigid particle a, G is the gravity vector of the rigid particle a, and the remaining symbols are the same as above;

[0134] Sediment solving module:

[0135] Step 4, based on the particle model obtained in step 1, introduce the Shield criterion of sediment particles, consider the flow-soil coupling force obtained in step 1 and the solid-soil coupling force obtained in step 2 to modify the sediment incipient model, and solve the foundation scour control; the specific steps are as follows:

[0136] Step 4.1, based on the fluid particle velocity u obtained in step 2, introduce the Einstein logarithmic flow velocity distribution formula to calculate the fluid shear stress τ experienced by the soil particles b :

[0137]

[0138] Where d is the characteristic particle size of the particle; Δu is the velocity difference between the sediment particle and the nearby water particle, κ is the von Karman constant, taken as 0.41, and ρ is the density of the sediment particle;

[0139] Step 4.2, based on the soil particle viscosity μ HBP obtained in step 1.1, introduce the Shield criterion of sediment particles to judge the starting state of the sediment particles, and calculate the sediment incipient critical stress τ cr,0 ;

[0140] τ bcr = θ cr ·(ρ s -ρ)gd

[0141] Where θ cr is the critical Hillz number, which is only related to the sediment parameters, ρ s is the saturated sediment density, ρ is the water density, g is the acceleration of gravity, and d is the particle size of the soil particles;

[0142] Step 4.3, based on the sediment incipient critical stress τ cr,0 obtained in step 4.2, introduce the flow-solid coupling force F fs of the sediment particles and the rigid particles obtained in step 2, and calculate the modified sediment incipient critical stress τ cr considering the external load and the slope effect;

[0143]

[0144]

[0145] Where P is the module of F fs , η = 0.7, is the unit vector in the flow direction, W is the gravity experienced by the sediment particles, α, β, and γ are the angles between the slope normal vector and the X, Y, and Z axis directions, w i is the component of in the X, Y, and Z axis directions, is the unit vector of gravity in the direction of the slope surface, is the internal friction angle of the sediment particle;

[0146] Step 4.4, based on the fluid shear stress τ b and the modified sediment incipient critical stress τ cr of step 4.3, determine the sediment particle incipient state, if τ cr ≥ τ b , then τ cr is substituted into the formula obtained in step 1.1 to replace the yield stress τ cr , and the bed sediment particle viscosity μ HBP is updated as a bed load particle, and if the condition is not met, it is processed according to steps 1.2 to 1.3;

[0147] Step 4.5, according to the bed load particle obtained in step 4.4, introduce the Mastbergen formula to calculate the critical flow velocity u lift of the bed load particle converted into suspended load:

[0148]

[0149] where α i is the sediment transport coefficient, n s is the normal vector of the bed surface, d * is the sediment particle size coefficient, ρ s is the saturated sediment density, ρ is the water density, g is the acceleration of gravity, d is the soil particle size, μ is the water viscosity, θ b is the actual Hillz number of the soil particle, θ cr is the critical Hillz number.

[0150] If the actual flow velocity of the fluid particle u ≥ u lift , it is converted into a suspended load particle, and its equivalent viscosity μ lift is calculated to replace μ HBP :

[0151]

[0152] where μ is the water viscosity, C v is the radius of the core function, and the concentration of the sediment particle.

[0153] If the condition is not met, it is not processed;

[0154] Step 4.6, according to the suspended load particle obtained in step 4.5, introduce the Mastbergen formula to calculate the critical flow velocity u set of the suspended load particle converted into bed load:

[0155]

[0156] wherein, μ is the viscosity of water body, d is the particle size of soil body, d * is the coefficient of sediment particle size.

[0157] If the actual flow velocity of the particle u satisfies u set , the particle is converted into a bed load particle and is processed according to step 4.4, and if the condition is not satisfied, the particle is not processed;

[0158] The above specific steps are only the specific description of the operation of the scour module for one time step, and in actual operation, the loop should be reciprocated until the solution termination time is reached.

[0159] The above method steps and basic formula principles can also be loaded on other open source SPH calculation software, so that a series of operation steps are executed on a computer or other programmable device to realize the scour numerical simulation, and the instructions executed on the computer or other programmable device provide the functions or steps specified in the flow Figure 1 of the flow or the plurality of flows.

[0160] The above examples only illustrate the technical idea of the present application, and cannot limit the protection scope of the present application, and any modification made according to the technical idea of the present application on the basis of the technical scheme falls within the protection scope of the present application.

Claims

1. A basic scouring fluid-solid-soil coupling simulation method based on SPH-DEM coupling and multiphase flow theory, characterized in that, Comprise the following steps: Step 1, construct particle model, set fluid particle and bed sand particle based on multiphase flow theory, set rigid particle based on DEM theory; The specific steps are as follows: Step 1.1, based on the basic principle of multiphase flow, respectively construct fluid particle model and soil particle model, set the fluid particle model as Newton fluid, set the soil particle model as non-Newton fluid, based on HBP model, determine the soil viscosity ; ， wherein, is the yield stress of the model material, is the second invariant of the fluid strain rate tensor, m is the stress exponent growth coefficient, is the viscosity of the water body, n is the power related to the shear stress; Step 1.

2. Introduce the DP yield criterion to calculate the yield stress of the specific material ; Step 1.3, substitute the soil viscosity calculated in step 1.1 into the formula of step 1.2 Step 1.3, substitute the soil viscosity calculated in step 1.1 into the formula of step 1.2 Step 1.3, substitute the soil viscosity calculated in step 1.1 into the formula of step 1.2 Step 1.3, substitute the soil viscosity calculated in step 1.1 into the formula of step 1.2 Step 1.4, set rigid particles based on DEM theory, establish normal contact stiffness between rigid particles , normal contact damping γ n , tangential contact stiffness , and tangential contact damping γ t ; Step 1.5: Based on the normal contact stiffness between rigid particles obtained in Step 1.4 Normal contact damping γ n Tangential contact stiffness and tangential contact damping γ t Calculate the normal contact force F n and tangential contact force F t ; Step 2, based on the particle model obtained in step 1, introduce the fluid-solid coupling theory to modify the fluid control equation, and carry out the control solution of the fluid particle, the specific steps are as follows: Step 2.1, introducing particle fluid coupling force F fs , modifying fluid momentum conservation equation; Step 2.2, based on the kernel function theory of SPH algorithm, the discrete form of the fluid momentum conservation equation obtained in step 2.1 is given: ， where is the Hamiltonian, the subscript denotes the central particle, the subscript denotes the neighboring particle, is the kernel function centered at the a particle , and denote the position vectors of the central and neighboring particles, respectively, is the ratio of the inter-particle distance to the smoothing length, is the smoothing length; Kinematic viscosity; and SPS is the stress vector for the fluid; respectively the mass and pressure of the fluid particles searched by the rigid particle s within the core function radius; equivalent mass and pressure for rigid particles, ; other variables as above; Step 3, based on the particle model obtained in step 1, combined with Newton's second law, considering the fluid-solid coupling force F obtained in step 2, the rigid particle control equation is modified and solved: the rigid particle control equation is modified as: fs , the rigid particle control equation is modified as: ， where m a , v a are the mass and velocity vector of rigid body particle a, respectively, G is the gravitational vector of rigid body particle a, and the remaining symbols are as above. Step 4, based on the particle model obtained in step 1, introduce the Shield criterion of sediment particles, consider the fluid-soil coupling force obtained in step 1 and the solid-soil coupling force obtained in step 2 to modify the sediment incipient model, and carry out the control solution of the foundation scour; The specific steps are as follows: Step 4.1, based on the momentum equation of step 2, obtain the fluid particle velocity u, introduce the Einstein logarithmic velocity distribution formula, calculate the fluid shear stress on the soil particles ; Step 4.2, Viscosity of soil particles based on the soil particles obtained in step 1.1 Model, introducing the Shield criterion for sediment particles, to calculate the incipient critical stress of the sediment ; Step 4.3, calculating the modified incipient motion critical stress of the sediment based on the sediment obtained in step 4.2 , introducing the fluid-solid coupling force F of the particle obtained in step 2 fs , calculating the modified incipient motion critical stress of the sediment considering the applied load and the effect of the slope ; Step 4.

4. Fluid shear stress based on the fluid obtained from Step 4.1 and the modified incipient sediment transport stress obtained from Step 4.3 Determine the incipient state of the sediment particles, if the condition is satisfied, then is substituted into the formula described in Step 1.1, replacing the yield stress and the bed sediment particle viscosity is updated to be the sediment transport particle, and if the condition is not satisfied, then the process is continued according to Steps 1.2 to 1.

3. Step 4.

5. Calculate the threshold velocity of the transformed suspension from the transformed granules obtained in step 4.4 according to the Mastbergen formula If the fluid velocity u , the suspension is transformed into granules and the equivalent viscosity is calculated and replaces the bed sand viscosity If the condition is not met, the sample is not processed; Step 4.

6. Calculate the threshold velocity of the transformed particles from the suspension particles obtained in step 4.5 according to the formula of Mastbergen If the actual velocity of the particles u is met, the particles are transformed and treated according to step 4.4, otherwise they are not treated. Step 5, repeat steps 1 to 4 until the solution is completed.

2. The simulation method according to claim 1, wherein, In step 1.2, the material yield stress The calculation is as follows: ， where, is the hydrostatic pressure acting on the saturated sediment particles, and are given by the Mohr-Coulomb yield criterion parameters, and the calculation formula is as follows: ， wherein: is the internal friction angle, is the soil cohesion.

3. The method of claim 1, wherein, The normal contact stiffness between rigid particles in step 1.4 , the normal contact damping γ n , the tangential contact stiffness , and the tangential contact damping γ t is calculated as follows: ， Among them, C n =1*10 -5 Let be the normal contact stiffness between rigid particle i and rigid particle j in the DEM. γ n,ij γ t,ij Synonyms; E * For the equivalent elastic modulus, R * M is the equivalent particle radius. * For equivalent mass, the calculations are as follows: ， wherein E i、 E j respectively the elastic modulus of the material assigned to DEM rigid particle i, rigid particle j. 、 respectively the Poisson's ratio of the material assigned to DEM rigid particle i, j. r i、 r j respectively the radius of DEM rigid particle i, rigid particle j; m i、 m j mi, mjare the mass of DEM rigid particle i, rigid particle j, respectively.

4. The method of claim 1, wherein, In the step 2.1, the momentum conservation equation of the fluid is modified as follows: ， where u is the fluid velocity vector, P is the fluid pressure term, is the fluid viscosity, g is the fluid gravity term, and the remaining symbols are as above.

5. The method of claim 1, wherein, In step 4.1, the fluid shear stress Calculated as: ， wherein, is the particle characteristic diameter; is the velocity difference between the sediment particle and the nearby water particle, is the von Karman constant, taken as 0.41, p d is the sediment particle density.

6. The method of claim 1, wherein, In step 4.2, the sediment incipient critical stress The calculation formula is: ， where, is the critical Shields number, which only depends on the sediment parameters, is the saturated sediment density, p is the water density, g is the gravitational acceleration, and d is the soil particle size.

7. The method of claim 1, wherein, In step 4.3, the threshold stress for incipient motion of sediment is corrected The calculation is: ， wherein , where P is F fs of the module, , is the flow velocity direction unit vector, W is the gravity experienced by the sediment particles, are the angles between the slope normal vector and the X, Y, Z axis directions, respectively, , , w i is the components in the X, Y, Z axis directions, is the unit vector of the gravity in the slope direction, is the internal friction angle of the sediment particles.

8. The method of claim 1, wherein, In step 4.5, the transition of the solid phase to the suspension phase occurs at the critical flow rate The calculation is: ， wherein, is the sediment transport coefficient, is the bed normal vector, is the sediment particle size coefficient, , is the saturated sediment density, p is the water density, g is the gravitational acceleration, d is the soil particle size, is the actual Hillz number of the soil particle, is the critical Hillz number; Equivalent viscosity Calculated as: ， wherein, viscosity of the water body, is the concentration of sediment particles within the kernel radius; In step 4.6, the critical flow rate The calculation formula is: ， wherein, is the water body viscosity, d is the soil particle size, is the sediment particle size coefficient.

Images