WO2024098534A1 - 一种基于sph-dem耦合与多相流理论的基础冲刷流-固-土耦合仿真方法 - Google Patents

一种基于sph-dem耦合与多相流理论的基础冲刷流-固-土耦合仿真方法 Download PDF

Info

Publication number
WO2024098534A1
WO2024098534A1 PCT/CN2022/142693 CN2022142693W WO2024098534A1 WO 2024098534 A1 WO2024098534 A1 WO 2024098534A1 CN 2022142693 W CN2022142693 W CN 2022142693W WO 2024098534 A1 WO2024098534 A1 WO 2024098534A1
Authority
WO
WIPO (PCT)
Prior art keywords
particle
fluid
particles
soil
sediment
Prior art date
Application number
PCT/CN2022/142693
Other languages
English (en)
French (fr)
Inventor
熊文
张嵘钊
Original Assignee
东南大学
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 东南大学 filed Critical 东南大学
Publication of WO2024098534A1 publication Critical patent/WO2024098534A1/zh

Links

Images

Definitions

  • the present invention relates to a scour simulation method based on SPH-DEM coupling and multiphase flow theory and considering fluid-solid coupling, and specifically relates to a method based on secondary development of SPH multiphase flow algorithm, introducing DEM calculation theory and multiple experimental criteria to perform real-time refined simulation of the whole process of flow-solid-soil coupling under scour, and belongs to the field of numerical simulation calculation.
  • Bridge foundation scour is one of the main reasons for the failure of bridge structure function and loss of its safety performance today, which has attracted widespread attention from many scholars.
  • Numerical simulation is still one of the most effective methods for studying bridge scour, with many advantages such as low cost, high efficiency and short cycle.
  • the mainstream algorithm currently used is to solve the N-S fluid control equation based on the Euler form to obtain fluid dynamics, and to obtain the numerical solution of bed elevation change by introducing a numerical bed sand transport model.
  • the numerical simulation method based on the Euler form is often difficult to converge when solving problems such as the breakage of complex fluid surfaces and wave and tsunami overflow, or requires a lot of solution time and resources, and there is no real soil model, making it difficult to track the evolution trajectory of the soil;
  • the existing numerical models basically treat the scour foundation as a fixed boundary, and there is almost no real-time analysis of the impact of scour hollowing on the stability of the foundation, nor is there a case that considers the influence of foundation extrusion pressure on sediment initiation.
  • the smoothed 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 guaranteed without additional calculations. When simulating complex free surface flows, there is no need to track fluid boundaries and interfaces between different fluids. When simulating scour problems, a soil model can be constructed without introducing a numerical sand transport model, which can further speed up efficiency and improve accuracy.
  • the discrete element (DEM) algorithm is also a numerical solution method based on the Lagrangian form. It is often used to simulate strong nonlinear mechanical behaviors such as large deformation and crushing of structures, and has good compatibility with SPH. However, at this stage, there are few scour simulation studies and application cases based on SPH-DEM coupling, and there is a lack of design solutions for related algorithms, which urgently need to be further resolved.
  • the technical problem to be solved by the present invention is: based on SPH-DEM coupling and multiphase flow theory, considering fluid-solid coupling calculation, introducing multiple theoretical or experimental criteria to carry out refined simulation of the whole process of fluid-solid-soil coupling under the background of foundation scour, so as to realize real-time analysis of the dynamic response of the foundation under scour and hollowing, and consider the influence of foundation extrusion pressure on sediment initiation, with the characteristics of easy program implementation, high accuracy, strong operability, etc.
  • a basic scour flow-solid-soil coupling simulation method based on SPH-DEM coupling and multiphase flow theory includes the following steps:
  • Step 1 construct a particle model, set the fluid particles and bed sand particles based on multiphase flow theory, and set the rigid body particles based on DEM theory; the specific steps are as follows:
  • 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 a Newtonian fluid, and the soil particle model is set as a non-Newtonian fluid, and the soil viscosity ⁇ HBP is determined based on the HBP model;
  • ⁇ 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 exponential growth coefficient
  • ⁇ is the viscosity of the water
  • n is the power related to the shear stress
  • Step 1.2 Introduce the DP yield criterion and calculate the specific material yield stress ⁇ y :
  • Step 1.3 substitute the specific material yield stress ⁇ y obtained in step 1.2 into the soil viscosity ⁇ HBP calculation formula described in step 1.1, replace the model material yield stress ⁇ c , and establish the soil viscosity calculation model;
  • Step 1.4 based on the DEM theory, rigid particles are set, and the normal contact stiffness K n , normal contact damping ⁇ n , tangential contact stiffness K t and tangential contact damping ⁇ t between rigid particles are established;
  • Step 1.5 based on the normal contact stiffness Kn , normal contact damping ⁇ n , tangential contact stiffness Kt and tangential contact damping ⁇ t between rigid particles obtained in step 1.4, calculate the normal contact force Fn and the tangential contact force Ft ;
  • Step 2 Based on the particle model obtained in step 1, the fluid-solid coupling theory is introduced to correct the fluid control equation and solve the control of the fluid particles.
  • the specific steps are as follows:
  • Step 2.1 introduce the particle fluid-solid coupling force F fs to modify the fluid momentum conservation equation
  • Step 2.2 based on the kernel function theory of the SPH algorithm, the discrete form of the fluid momentum conservation equation obtained in step 2.1 is given:
  • 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 obtained in step 2, modify the rigid body particle control equation and solve it:
  • Step 4 Based on the particle model obtained in step 1, the Shield criterion for sediment particles is introduced, and the sediment initiation model is corrected by considering the flow-soil coupling force obtained in step 1 and the solid-soil coupling force obtained in step 2 to solve the foundation scour control.
  • the specific steps are as follows:
  • Step 4.1 Based on the momentum equation in step 2, obtain the fluid particle velocity u, introduce the Einstein logarithmic velocity distribution formula, and calculate the fluid shear stress ⁇ b on the soil particles:
  • Step 4.2 Based on the soil particle viscosity ⁇ HBP model obtained in step 1.1, the Shield criterion for sediment particles is introduced to calculate the sediment critical starting stress ⁇ cr,0 ;
  • Step 4.3 Based on the sediment critical starting stress ⁇ cr,0 obtained in step 4.2, the flow-solid coupling force F fs between sediment particles and rigid particles obtained in step 2 is introduced to calculate the modified sediment critical starting stress ⁇ cr considering the external load and slope effect;
  • Step 4.4 Determine the starting state of sediment particles based on the fluid shear stress ⁇ b obtained in step 4.1 and the modified sediment starting critical stress ⁇ cr obtained in step 4.3. If ⁇ cr ⁇ ⁇ b is satisfied, substitute ⁇ cr into the formula described in step 1.1, replace the yield stress ⁇ cr , and update the bed sand particle viscosity ⁇ HBP as the bed load particle. If the conditions are not satisfied, proceed to step 1.2 to step 1.3;
  • Step 4.5 Based on the bed load particles obtained in step 4.4, the Mastbergen formula is introduced to calculate the critical velocity u lift of bed load conversion into suspended mass. If the fluid particle velocity u ⁇ u lift is satisfied, it is converted into suspended mass particles, and the equivalent viscosity ⁇ lift is calculated to replace the bed sand particle viscosity ⁇ HBP . If the condition is not satisfied, it will not be processed;
  • Step 4.6 Based on the suspended mass particles obtained in step 4.5, the Mastbergen formula is introduced to calculate the critical velocity u set of suspended mass to bed mass. If the actual velocity of the particles u ⁇ u set is satisfied, the suspended mass is converted into bed mass particles and processed according to step 4.4. If the condition is not satisfied, the suspended mass is not processed.
  • Step 5 Repeat steps 1 to 4 until the solution is completed.
  • step 1.2 the material yield stress ⁇ y is calculated as follows:
  • ⁇ is the internal friction angle
  • c is the soil cohesion
  • step 1.4 the normal contact stiffness K n , normal contact damping ⁇ n , tangential contact stiffness K t and tangential contact damping ⁇ t between rigid particles are calculated as follows:
  • C n 1*10 -5
  • K n,ij is the normal contact stiffness between DEM rigid particle i and rigid particle j
  • K t,ij , ⁇ n,ij , ⁇ t,ij are synonymous
  • E * is the equivalent elastic modulus
  • R * is the equivalent particle radius
  • M * is the equivalent mass
  • E i and E j are the elastic moduli of the materials assigned to the DEM rigid particles i and j, respectively;
  • ⁇ i and ⁇ j are the Poisson’s ratios of the materials assigned to the DEM rigid particles i and j, respectively;
  • ri and rj are the radii of DEM rigid particle i and rigid particle j respectively;
  • m i and m j are the masses of DEM rigid particle i and rigid particle j respectively.
  • step 2.1 the momentum conservation equation of the fluid is modified to:
  • u is the fluid velocity vector
  • P is the fluid pressure term
  • ⁇ is the fluid viscosity
  • g is the fluid gravity term
  • step 2.2 the discrete form of the control equation obtained in step 2.1 is:
  • ra and rb represent the position vectors of the central particle and the neighboring particles, respectively, q is the ratio of the interparticle distance to the smoothing length, and 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;
  • step 3 the rigid body particle control equation is corrected to:
  • ma and va are the mass and velocity vector of rigid particle a respectively
  • G is the gravity vector of rigid particle a
  • the other symbols are the same as above.
  • step 4.1 the fluid shear stress ⁇ b is calculated as:
  • d is the characteristic particle size
  • ⁇ u is the velocity difference between the sediment particles and the nearby water particles
  • ⁇ is the von Karman constant, which is taken as 0.41
  • ⁇ is the density of the sediment particles.
  • step 4.2 the calculation formula for the sediment critical starting stress ⁇ cr,0 is:
  • ⁇ cr,0 ⁇ cr ⁇ ( ⁇ s - ⁇ )gd
  • ⁇ cr is the critical Shields number, which is only related to the sediment parameters
  • ⁇ S is the saturated sediment density
  • ⁇ is the water density
  • g is the gravitational acceleration
  • d is the soil particle size.
  • step 4.3 the modified sediment critical stress ⁇ cr is calculated as:
  • w i The components in the X, Y, and Z axis directions, is the unit vector of gravity in the slope direction, is the internal friction angle of sediment particles.
  • step 4.5 the critical flow rate u lift for converting bed load into suspended mass is calculated as:
  • ⁇ i is the sediment transport coefficient
  • ns is the bed surface normal vector
  • d * is the sediment particle size coefficient
  • ⁇ s is the saturated sediment density
  • ⁇ is the water density
  • g is the gravitational acceleration
  • d is the soil particle size
  • ⁇ b is the actual Shields number of the soil particles
  • ⁇ cr is the critical Shields number.
  • ⁇ is the viscosity of water
  • C v is the concentration of sediment particles within the radius of the kernel function
  • step 4.6 the critical flow rate u set is calculated as:
  • ⁇ is the viscosity of water
  • d is the particle size of soil particles
  • d * is the sediment particle size coefficient
  • the present invention is based on SPH-DEM coupling and multiphase flow theory, considers fluid-solid coupling calculation, introduces multiple theoretical or experimental criteria to perform real-time and refined simulation of the entire process of fluid-solid-soil coupling under the background of foundation scour, thereby realizing real-time analysis of the foundation dynamic response under scour, and considering the influence of foundation extrusion pressure on sediment initiation, with the characteristics of easy program implementation, high accuracy, strong operability, etc.
  • the present invention realizes numerical simulation of scour based on secondary development of SPH multiphase flow.
  • the particles themselves have mass, and mass conservation can be guaranteed without additional calculation.
  • simulating scour problems a soil model can be constructed without introducing a numerical sand transport model, which can further speed up efficiency and improve accuracy.
  • the present invention truly realizes the real-time simulation of local scour calculation and foundation dynamics and stability analysis.
  • FIG1 is a calculation flow chart of a basic scour flow-solid-soil coupling simulation method based on SPH-DEM coupling and multiphase flow theory of the present invention
  • FIG2 is a diagram of the fluid particle distribution and retrieval mechanism simulated by the present invention.
  • FIG3 is a diagram of the DEM rigid particle contact mechanism simulated by the present invention.
  • FIG4 is a diagram of the classification of sediment types for the scour calculation simulated by the present invention.
  • Figure 5 is an enlarged view of the interface between bed load and suspended load.
  • 1 is the fluid particle
  • 3 is the kernel function radius h
  • the present invention is a secondary development based on the SPH multiphase flow algorithm and DEM calculation theory. Therefore, the default premise of the present invention is that the SPH multiphase flow algorithm and 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 the present invention, and the relevant variables can be manually extracted or obtained.
  • the fluid particle distribution and retrieval mechanism based on the present invention is shown in FIG2 , where 1 is a fluid particle, 2 is a central fluid particle to be searched, and 3 is a kernel function radius h.
  • the DEM rigid particle contact model on which the present invention is based 4 is the center of mass of the DEM rigid particle i, 5 is the center of mass of the DEM rigid particle j, 6 is the normal spring in the contact portion between the particles i and j, whose stiffness is K n , and the radial contact length is the radial length of the overlapping portion, 7 is the normal contact damping ⁇ n in the contact portion between the particles i and j, 8 is the tangential spring in the contact portion between the particles i and j, whose stiffness is K t , and the tangential contact length is the tangential length of the overlapping portion, and 9 is the tangential contact damping ⁇ t in the contact portion between the particles i and j.
  • the sediment distribution on which the present invention is based is that the riverbed portion is 10 ordinary soil, the scour pit is 11 yield soil, and the inside of the scour pit is distributed with 12 bedload and 13 suspended sediment, the 12 bedload is adjacent to the 13 yield soil, and the 13 suspended sediment is distributed above the 12 bedload.
  • Figure 5 is an enlarged view of the interface between bed load and suspended load.
  • the specific implementation method of the present invention takes the SPH open source computing software Dualsphysics as an example.
  • the solution of partial differential equations such as fluid control equations and DEM motion control equations is processed by corresponding algorithm modules, and the equation variables can be artificially extracted.
  • Step 1 construct a particle model, set the fluid particles and bed sand particles based on multiphase flow theory, and set the rigid body particles based on DEM theory; the specific steps are as follows:
  • Step 1.1 Based on the basic principle of multiphase flow, construct the fluid particle model and the soil particle model respectively.
  • the fluid particle model is set as Newtonian fluid
  • the soil particle model is set as non-Newtonian fluid.
  • Based on the HBP model determine the soil viscosity ⁇ HBP :
  • ⁇ 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 exponential growth coefficient
  • ⁇ is the viscosity of the water
  • n is the power related to the shear stress
  • Step 1.2 Introduce the DP yield criterion and calculate the specific material yield stress ⁇ y :
  • ⁇ is the internal friction angle
  • c is the soil cohesion
  • Step 1.3 substitute the specific material yield stress ⁇ y obtained in step 1.2 into the soil viscosity ⁇ HBP calculation formula described in step 1.1, replace the model material yield stress ⁇ c , and establish the soil viscosity calculation model;
  • Step 1.4 Set rigid particles based on DEM theory, and establish the normal contact stiffness K n , normal contact damping ⁇ n , tangential contact stiffness K t and tangential contact damping ⁇ t between rigid particles:
  • C n 1*10 -5
  • K n,ij is the normal contact stiffness between DEM rigid particle i and rigid particle j
  • K t,ij , ⁇ n,ij , ⁇ t,ij are synonymous
  • E * is the equivalent elastic modulus
  • R * is the equivalent particle radius
  • M * is the equivalent mass
  • E i and E j are the elastic moduli of the materials assigned to the DEM rigid particles i and j
  • ⁇ i and ⁇ j are the Poisson's ratios of the materials assigned to the DEM rigid particles i and j
  • ri and r j are the radii of the DEM rigid particles i and j
  • mi and m j are the masses of the DEM rigid particles i and j, respectively.
  • Step 1.5 Determine the normal contact force Fn and the tangential contact force Ft according to the contact stiffness between rigid particles:
  • ⁇ ij is the radial and tangential contact distance between DEM rigid particle i and rigid particle j;
  • e ij is the unit vector from the mass center of DEM rigid particle i to the mass center of rigid particle j; are the radial and tangential deformation rates of DEM rigid particle i and rigid particle j, respectively.
  • v ij ⁇ is the radial and tangential relative velocity of DEM rigid particle i and rigid particle j
  • u f is the friction coefficient between DEM rigid particle i and rigid particle j;
  • Step 2 Based on the particle model obtained in step 1, the fluid-solid coupling theory is introduced to correct the fluid control equation and solve the control of the fluid particles.
  • the specific steps are:
  • Step 2.1 introduce the particle fluid-solid coupling force F fs , and the momentum conservation equation of the fluid is modified to:
  • u is the fluid velocity vector
  • P is the fluid pressure term
  • ⁇ is the fluid viscosity
  • g is the fluid gravity term
  • Step 2.2 based on the kernel function theory of the SPH algorithm, the discrete form of the control equation obtained in step 2.1 is given:
  • ra and rb represent the position vectors of the central particle and the neighboring particle, respectively;
  • q is the ratio of the particle spacing to the smooth length;
  • h is the smooth length;
  • ⁇ 0 is the kinematic viscosity;
  • ⁇ ij is the fluid SPS stress vector;
  • mf and Pf are the mass and pressure of the fluid particles searched by the rigid particle s within the kernel function radius;
  • ms and Ps are the equivalent mass and pressure of the rigid particle, respectively.
  • Step 3 Based on the particle model obtained in step 1, combined with Newton's second law and introducing the particle fluid-solid coupling force F fs obtained in step 2, the rigid body particle control equation is corrected and solved:
  • ma and va are the mass and velocity vectors of rigid particle a respectively
  • G is the gravity vector of rigid particle a, and the other symbols are the same as above;
  • Step 4 Based on the particle model obtained in step 1, the Shield criterion for sediment particles is introduced, and the sediment initiation model is corrected by considering the flow-soil coupling force obtained in step 1 and the solid-soil coupling force obtained in step 2 to solve the foundation scour control.
  • the specific steps are as follows:
  • Step 4.1 Based on the fluid particle velocity u obtained in step 2, the Einstein logarithmic velocity distribution formula is introduced to calculate the fluid shear stress ⁇ b on the soil particles:
  • d is the characteristic particle size
  • ⁇ u is the velocity difference between the sediment particle and the nearby water particles
  • ⁇ is the von Karman constant, which is 0.41
  • ⁇ is the density of the sediment particle.
  • Step 4.2 Based on the soil particle viscosity ⁇ HBP model obtained in step 1.1, the Shield criterion of sediment particles is introduced to determine the starting state of sediment particles and calculate the sediment starting critical stress ⁇ cr,0 ;
  • ⁇ cr is the critical Shields number, which is only related to the sediment parameters, ⁇ s is the saturated sediment density, ⁇ is the water density, g is the gravitational acceleration, and d is the soil particle size;
  • Step 4.3 Based on the sediment critical starting stress ⁇ cr,0 obtained in step 4.2, the flow-solid coupling force F fs between sediment particles and rigid particles obtained in step 2 is introduced to calculate the modified sediment critical starting stress ⁇ cr considering the external load and slope effect;
  • w i The components in the X, Y, and Z axis directions, is the unit vector of gravity in the slope direction, is the internal friction angle of sediment particles;
  • Step 4.4 Based on the fluid shear stress ⁇ b obtained in step 4.1 and the modified sediment starting critical stress ⁇ cr obtained in step 4.3, the starting state of the sediment particles is judged. If ⁇ cr ⁇ ⁇ b is satisfied, ⁇ cr is substituted into the formula obtained in step 1.1 to replace the yield stress ⁇ cr , and the bed sand particle viscosity ⁇ HBP is updated to be used as the bed load particle. If the conditions are not satisfied, the process is carried out according to steps 1.2 to 1.3;
  • Step 4.5 Based on the bed load particles obtained in step 4.4, the Mastbergen formula is introduced to calculate the critical velocity u lift of bed load conversion into suspended mass:
  • ⁇ i is the sediment transport coefficient
  • ns is the bed surface normal vector
  • d * is the sediment particle size coefficient
  • ⁇ s is the saturated sediment density
  • ⁇ is the water density
  • g is the gravitational acceleration
  • d is the soil particle size
  • ⁇ is the water viscosity
  • ⁇ b is the actual Shields number of the soil particles
  • ⁇ cr is the critical Shields number.
  • ⁇ is the viscosity of water
  • Cv is the concentration of sediment particles within the kernel function radius.
  • Step 4.6 Based on the suspended mass particles obtained in step 4.5, the Mastbergen formula is introduced to calculate the critical flow rate u set of suspended mass conversion bedload:
  • ⁇ is the viscosity of water
  • d is the particle size of soil particles
  • d * is the sediment particle size coefficient
  • step 4.4 If the actual particle velocity u ⁇ u set is satisfied, it is converted into a bed load particle and processed according to step 4.4. If the condition is not satisfied, it is not processed.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Geometry (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Computer Hardware Design (AREA)
  • Computing Systems (AREA)
  • General Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computational Linguistics (AREA)
  • Data Mining & Analysis (AREA)
  • Biomedical Technology (AREA)
  • Molecular Biology (AREA)
  • Artificial Intelligence (AREA)
  • Biophysics (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Health & Medical Sciences (AREA)

Abstract

本发明公开了一种基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真算法，包括：构建粒子模型，基于多相流理论设定流体粒子及床沙粒子，基于DEM理论设定刚体粒子；基于流-固耦合理论修正流体控制方程，进行流体粒子的控制求解；基于牛顿第二定律并引入流-固耦合力，进行DEM刚体粒子的控制求解；基于泥沙颗粒Shield准则，引入流-土、土-固耦合力修正泥沙起动模型，进行基础冲刷控制求解；完成一个时间步长，进入下一次循环；本发明通过基于SPH-DEM耦合及多相流理论的二次开发，引入多种结构、状态模型，可实现基础冲刷流-固-土耦合的精细化数值模拟。

Description

一种基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真方法 技术领域
本发明涉及一种基于SPH-DEM耦合与多相流理论、考虑流固耦合的冲刷仿真方法，具体涉及一种基于SPH多相流算法二次开发，引入DEM计算理论及多种实验准则进行冲刷下流-固-土耦合的全过程实时精细化模拟，属于数值模拟计算领域。
背景技术
桥梁基础冲刷病害是当今桥梁结构功能失效、丧失其安全性能的最主要原因之一，引起了众多学者的广泛关注。数值模拟仍然是研究桥梁冲刷最有效的方法之一，兼具成本低、效率高、周期短等多项优势。现今主流采用的算法为基于欧拉形式求解N-S流体控制方程获得流体动力，并通过引入数值化床沙输运模型，得到床面高程变化的数值解。然而基于欧拉形式的数值模拟方法在解决复杂流体表面的破碎以及波浪、海啸越流等问题时往往难以收敛，或需要大量求解时间与资源，且不存在真实的土体模型，难以追踪土体演变轨迹；现有的数值模型基本将冲刷基础作为固定边界处理，几乎未见实时分析冲刷掏空对基础稳定性的研究，亦未见考虑基础挤压力对泥沙起动影响的案例，缺少研究基础冲刷背景下，模拟流-固-土多场耦合效应的有效工具。
光滑粒子流体动力学(SPH)方法，是基于拉格朗日形式的数值求解法，与欧拉法相比，粒子本身具有质量，无需额外计算就能保证质量守恒；模拟复杂自由表面流时，无需追踪流体边界及不同流体交界面。在模拟冲刷问题时，可构建土体模型，无须引入数值化的输沙模型，可进一步加快效率并提高精度；离散元(DEM)算法同属于基于拉格朗日形式的数值求解法，常用来进行结构物大变形、破碎等强非线性力学行为的模拟，与SPH有较好的兼容性。然而，现阶段基于SPH-DEM耦合的冲刷模拟研究、应用案例较少，缺少相关算法的设计方案，亟待进一步解决。
发明内容
本发明所要解决的技术问题是：基于SPH-DEM耦合与多相流理论、考虑流-固耦合计算，引入多项理论或实验准则进行基础冲刷背景下流-固-土耦合全过程的精细化模 拟，从而实现实时分析冲刷掏空下基础动力响应，并考虑基础挤压力对泥沙起动影响，具有易程序化实现，精确度高，可操作性强等特点。
本发明为解决上述技术问题采用以下技术方案：
一种基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真方法，包括如下步骤：
步骤1，构建粒子模型，基于多相流理论设定流体粒子及床沙粒子，基于DEM理论设定刚体粒子；具体步骤如下：
步骤1.1、基于多相流基本原理，分别构建流体粒子模型与土体粒子模型，流体粒子模型设定为牛顿流体，土体粒子模型设定为非牛顿流体，基于HBP模型，确定土体粘度μ HBP；
Figure PCTCN2022142693-appb-000001
其中，τ c为模型材料屈服应力，II D为流体应变率张量的第二不变量，m为应力指数增长系数，μ为水体粘度，n为与剪切应力相关的幂；
步骤1.2、引入DP屈服准则，计算具体材料屈服应力τ y：
步骤1.3、将步骤1.2所得具体材料屈服应力τ y代入步骤1.1所述土体粘度μ HBP计算公式，替换模型材料屈服应力τ c，确立土体粘度计算模型；
步骤1.4、基于DEM理论设定刚体粒子，确立刚体粒子间法向接触刚度K n、法向接触阻尼γ n、切向接触刚度K t以及切向接触阻尼γ t；
步骤1.5、基于步骤1.4所得刚体粒子间法向接触刚度K n、法向接触阻尼γ n、切向接触刚度K t以及切向接触阻尼γ t，计算法向接触力F n及切向接触力F t；
步骤2，基于步骤1所得粒子模型，引入流-固耦合理论修正流体控制方程，进行流体粒子的控制求解，具体步骤为：
步骤2.1，引入粒子流固耦合力F fs，修正流体动量守恒方程；
步骤2.2，基于SPH算法的核函数理论，给出步骤2.1所得流体动量守恒方程的离散形式：
步骤3，基于步骤1所得粒子模型，结合牛顿第二定律，考虑步骤2所得粒子流固耦合力F fs,，修正刚体粒子控制方程并求解：
步骤4，基于步骤1所得粒子模型，引入泥沙颗粒Shield准则，考虑步骤1所得流-土耦合力、步骤2所得固-土耦合力修正泥沙起动模型，进行基础冲刷控制求解；具体步骤如下：
步骤4.1、基于步骤2动量方程获得流体粒子速度u，引入Einstein对数流速分布公式，计算土体粒子所受流体切应力τ b：
步骤4.2、基于步骤1.1所得土体粒子粘度μ HBP模型，引入泥沙颗粒Shield准则，计算泥沙起动临界应力τ cr,0；
步骤4.3、基于步骤4.2所得泥沙起动临界应力τ cr,0，引入步骤2所得泥沙粒子与刚体粒子流-固耦合力F fs,计算考虑外加荷载及边坡效应的修正泥沙起动临界应力τ cr；
步骤4.4、基于步骤4.1所得流体切应力τ b及步骤4.3所得修正泥沙起动临界应力τ cr判断泥沙粒子起动状态，若满足τ cr≥τ b，则将τ cr代入步骤1.1所述公式，替换屈服应力τ cr，更新床沙粒子粘度μ HBP以作为推移质粒子，不满足条件则按步骤1.2至步骤1.3处理；
步骤4.5、根据步骤4.4所得推移质粒子，引入Mastbergen公式计算推移质转化悬疑质临界流速u lift，若满足流体粒子流速u≥u lift,则转化为悬疑质粒子，并计算等效粘度μ lift替换床沙粒子粘度μ HBP，若不满足条件，则不予处理；
步骤4.6、根据步骤4.5所得悬疑质粒子，引入Mastbergen公式，计算悬疑质转化推移质临界流速u set，若满足粒子实际流速u≤u set，则转化为推移质粒子，并按步骤4.4处理，若不满足条件，则不予处理；
步骤5，重复步骤1至步骤4，直至求解完成。
所述步骤1.2中，材料屈服应力τ y计算如下：
|τ y|＝αp+β
其中，p为饱和泥沙粒子上作用的静水压力，α与β均由Mohr-Coulomb屈服准则参数给出，计算公式如下：
Figure PCTCN2022142693-appb-000002
其中：θ为内摩擦角，c为土体粘聚力。
所述步骤1.4中，刚体粒子间法向接触刚度K n、法向接触阻尼γ n、切向接触刚度K t 以及切向接触阻尼γ t计算如下：
Figure PCTCN2022142693-appb-000003
其中，C n＝1*10 -5,K n,ij为DEM刚体粒子i与刚体粒子j之间的法向接触刚度，K t,ij、γ n,ij、γ t,ij同义；E *为等效弹性模量，R *为等效粒子半径，M *为等效质量，计算分别如下：
Figure PCTCN2022142693-appb-000004
其中，E i、E j分别为DEM刚体粒子i、刚体粒子j所赋予材料的弹性模量；
μ i、μ j分别为DEM刚体粒子i、刚体粒子j所赋予材料的泊松比；
r i、r j分别为DEM刚体粒子i、刚体粒子j的半径；
m i、m j分别为DEM刚体粒子i、刚体粒子j的质量。
所述步骤2.1中，流体的动量守恒方程修正为：
Figure PCTCN2022142693-appb-000005
式中，u为流体速度矢量，P为流体压力项，υ为流体粘度，g为流体重力项，其余符号同上。
所述步骤2.2中，步骤2.1所得控制方程的离散形式为：
Figure PCTCN2022142693-appb-000006
其中，
Figure PCTCN2022142693-appb-000007
为哈密顿算子，下标a表示中心粒子，下标b表示邻域粒子，W ab为以a粒子为检索中心的核函数
Figure PCTCN2022142693-appb-000008
r a和r b分别表示中心粒子和邻域粒子的位置矢量,q为粒子间距与平滑长度的比值,h为光滑长度；
μ 0为运动粘度；
τ ij为流体SPS应力矢量；
m f、P f分别为刚体粒子s在核函数半径内搜索的流体粒子的质量与压力；
m s、P s为刚体粒子等效质量和压力，
Figure PCTCN2022142693-appb-000009
其余变量同上。
所述步骤3中，刚体粒子控制方程修正为：
Figure PCTCN2022142693-appb-000010
其中，m a、v a分别为刚体粒子a的质量和速度矢量，G为刚体粒子a的重力矢量，其余符号同上。
所述步骤4.1中，流体切应力τ b计算为：
Figure PCTCN2022142693-appb-000011
其中，d为粒子特征粒径；Δu为泥沙粒子与附近水粒子的速度差，κ为冯·卡曼常数，取0.41，ρ为泥沙粒子密度。
所述步骤4.2中，泥沙起动临界应力τ cr,0计算公式为：
τ cr,0＝θ cr·(ρ s-ρ)gd
其中，θ cr为临界希尔兹数,仅与泥沙参数有关，ρ S为饱和泥沙密度，ρ为水体密度，g为重力加速度，d为土体粒子粒径。
所述步骤4.3中，修正泥沙起动临界应力τ cr计算为：
Figure PCTCN2022142693-appb-000012
Figure PCTCN2022142693-appb-000013
其中，P为F fs的模，η＝0.7，
Figure PCTCN2022142693-appb-000014
为流速方向单位向量，W为泥沙粒子所受重力，α,β,和γ分别为坡面法向量与X,Y,Z轴方向的夹角，
Figure PCTCN2022142693-appb-000015
w i为
Figure PCTCN2022142693-appb-000016
在X,Y,Z轴方向的分量，
Figure PCTCN2022142693-appb-000017
为重力在坡面方向的分项单位向量，
Figure PCTCN2022142693-appb-000018
为泥沙粒子内摩擦角。
所述步骤4.5中，推移质转化悬疑质临界流速u lift计算为：
Figure PCTCN2022142693-appb-000019
其中，α i为泥沙输运系数，n s为床面法向量，d *为泥沙粒径系数，
Figure PCTCN2022142693-appb-000020
Figure PCTCN2022142693-appb-000021
ρ s为饱和泥沙密度，ρ为水体密度，g为重力加速度，d为土体粒子粒径，θ b为土体粒子实际希尔兹数，θ cr为临界希尔兹数。
等效粘度μ lift计算为：
Figure PCTCN2022142693-appb-000022
其中，μ为水体粘度，C v为核函数半径内泥沙粒子浓度；
所述步骤4.6中，临界流速u set计算公式为：
Figure PCTCN2022142693-appb-000023
其中，μ为水体粘度，d为土体粒子粒径，d *为泥沙粒径系数。
发明采用以上技术方案与现有技术相比，具有以下技术效果：
1.本发明基于SPH-DEM耦合与多相流理论、考虑流-固耦合计算，引入多项理论或实验准则进行基础冲刷背景下流-固-土耦合全过程的实时精细化模拟，从而实现实时分析冲刷下基础动力响应，并考虑基础挤压力对泥沙起动影响，具有易程序化实现，精确度高，可操作性强等特点。
2、本发明基于SPH多相流二次开发实现冲刷数值模拟，粒子本身具有质量，无需额外计算就能保证质量守恒；在模拟冲刷问题时，可构建土体模型，无须引入数值化的输沙模型，可进一步加快效率并提高精度。
3.本发明真正意义上实现了局部冲刷计算及基础动力和稳定性分析的实时模拟。
附图说明
图1是本发明一种基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真 方法的计算流程图；
图2是本发明所模拟的流体粒子分布及检索机制图解；
图3是本发明所模拟的DEM刚体粒子接触机制图解；
图4是本发明所模拟的冲刷计算的泥沙类型分类图解；
图5是推移质与悬移质交界面放大图。
其中，1为流体粒子；
2为待搜索的中心流体粒子；
3为核函数半径h；
4为DEM刚体粒子i的质心；
5为DEM刚体粒子j的质心；
6为粒子i与粒子j接触部分中的法向弹簧，其刚度为K n，径向接触长度即为重合部分径向长度；
7为粒子i与粒子j接触部分中的法向接触阻尼γ n；
8为粒子i与粒子j接触部分中的切向弹簧，其刚度为K t，切尔向接触长度即为重合部分切向长度；
9为粒子i与粒子j接触部分中的切向接触阻尼γ t；
10为普通土体；
11为屈服土体；
12为推移质；
13为悬移质。
具体实施方式
下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。
本发明是以SPH多相流算法、DEM计算理论为基础的二次开发。因此，本发明默认的前提条件为，SPH多相流算法、DEM计算理论已知或开源，流体控制方程及DEM粒子运动控制方程等偏微分方程的求解不在本发明的讨论范围之内，相关变量可人工提取或获得。
本发明所基于的流体粒子分布及检索机制如图2所示，1为流体粒子，2为待搜索 的中心流体粒子，3为核函数半径h。
如图3所示，本发明所基于的DEM刚体粒子接触模型，4为DEM刚体粒子i的质心，5为DEM刚体粒子j的质心，6为粒子i与粒子j接触部分中的法向弹簧，其刚度为K n，径向接触长度即为重合部分径向长度，7为粒子i与粒子j接触部分中的法向接触阻尼γ n，8为粒子i与粒子j接触部分中的切向弹簧，其刚度为K t，切尔向接触长度即为重合部分切向长度，9为粒子i与粒子j接触部分中的切向接触阻尼γ t。
如图4所示，本发明所基于的泥沙分布情况，河床部分为10普通土体，冲刷坑下为11屈服土体，冲刷坑内部分布有12推移质和13悬移质，12推移质紧邻13屈服土体，13悬移质分布在12推移质之上。
图5是推移质与悬移质交界面放大图。
本发明具体实施方法以SPH开源计算软件Dualsphysics为例，流体控制方程及DEM运动控制方程等偏微分方程的求解由相应算法模块处理，方程变量可人为提取获得。
结合图1，本发明一种基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真方法具体计算流程如下：
前处理模块：
步骤1，构建粒子模型，基于多相流理论设定流体粒子及床沙粒子，基于DEM理论设定刚体粒子；具体步骤如下：
步骤1.1、基于多相流基本原理，分别构建流体粒子模型与土体粒子模型，流体粒子模型设定为牛顿流体，土体粒子模型设定为非牛顿流体，基于HBP模型，确定土体粘度μ HBP：
Figure PCTCN2022142693-appb-000024
其中，τ c为模型材料屈服应力，II D为流体应变率张量的第二不变量，m为应力指数增长系数，μ为水体粘度，n为与剪切应力相关的幂；
步骤1.2、引入DP屈服准则，计算具体材料屈服应力τ y：
|τ y|＝αp+β
其中，p为饱和泥沙粒子上作用的静水压力，α与β均由Mohr-Coulomb屈服准则参数给出，计算如下：
Figure PCTCN2022142693-appb-000025
其中：θ为内摩擦角，c为土体粘聚力；
步骤1.3、将步骤1.2所得具体材料屈服应力τ y代入步骤1.1所述土体粘度μ HBP计算公式，替换模型材料屈服应力τ c，确立土体粘度计算模型；
步骤1.4、基于DEM理论设定刚体粒子，确立刚体粒子间法向接触刚度K n、法向接触阻尼γ n、切向接触刚度K t以及切向接触阻尼γ t:
Figure PCTCN2022142693-appb-000026
其中，C n＝1*10 -5,K n,ij为DEM刚体粒子i与刚体粒子j之间的法向接触刚度，K t,ij、γ n,ij、γ t,ij同义；E *为等效弹性模量，R *为等效粒子半径，M *为等效质量，计算分别如下：
Figure PCTCN2022142693-appb-000027
其中，E i、E j分别为DEM刚体粒子i、刚体粒子j所赋予材料的弹性模量，μ i、μ j分别为DEM刚体粒子i、刚体粒子j所赋予材料的泊松比，r i、r j分别为DEM刚体粒子i、刚体粒子j的半径，m i、m j分别为DEM刚体粒子i、刚体粒子j的质量；
步骤1.5、根据刚体粒子间的接触刚度确定法向接触力F n及切向接触力F t:
Figure PCTCN2022142693-appb-000028
其中，δ ij、
Figure PCTCN2022142693-appb-000029
为DEM刚体粒子i和刚体粒子j粒子径向及切向接触距离；e ij为DEM刚体粒子i质心指向刚体粒子j粒子质心的方向单位向量；
Figure PCTCN2022142693-appb-000030
分别为DEM刚体粒子i和刚体粒子j粒子径向及切向的变形速率，
Figure PCTCN2022142693-appb-000031
v ij、
Figure PCTCN2022142693-appb-000032
为DEM刚体粒子i和刚体粒子j粒子径向及切向相对速度，u f为DEM刚体粒子i和刚体粒子j间的摩擦系数；
流体求解模块：
步骤2，基于步骤1所得粒子模型，引入流-固耦合理论修正流体控制方程，进行流体粒子的控制求解；具体步骤为：
步骤2.1，引入粒子流固耦合力F fs，流体的动量守恒方程修正为：
Figure PCTCN2022142693-appb-000033
式中，u为流体速度矢量，P为流体压力项，υ为流体粘度，g为流体重力项，其余符号同上；
步骤2.2，基于SPH算法的核函数理论，给出步骤2.1所得控制方程的离散形式：
Figure PCTCN2022142693-appb-000034
其中，
Figure PCTCN2022142693-appb-000035
为哈密顿算子，下标a表示中心粒子，下标b表示邻域粒子，W ab为以a粒子为检索中心的核函数
Figure PCTCN2022142693-appb-000036
r a和r b分别表示中心粒子和邻域粒子的位置矢量,q为粒子间距与平滑长度的比值,h为光滑长度；μ 0为运动粘度，τ ij为流体SPS应力矢量，m f、P f分别为刚体粒子s在核函数半径内搜索的流体粒子的质量与压力，m s、P s为刚体粒子等效质量和压力，
Figure PCTCN2022142693-appb-000037
Figure PCTCN2022142693-appb-000038
DEM控制求解模块：
步骤3，基于步骤1所得粒子模型，结合牛顿第二定律并引入步骤2所得粒子流固耦合力F fs,，修正刚体粒子控制方程并进项求解：
Figure PCTCN2022142693-appb-000039
其中，m a、v a分别为刚体粒子a的质量和速度矢量，G为刚体粒子a的重力矢量，其余符号同上；
泥沙求解模块：
步骤4，基于步骤1所得粒子模型，引入泥沙颗粒Shield准则，考虑步骤1所得流 -土耦合力、步骤2所得固-土耦合力修正泥沙起动模型，进行基础冲刷控制求解；具体步骤如下：
步骤4.1、基于步骤2所得流体粒子速度u，引入Einstein对数流速分布公式，计算土体粒子所受流体切应力τ b：
Figure PCTCN2022142693-appb-000040
其中，d为粒子特征粒径；Δu为泥沙粒子与附近水粒子的速度差，κ为冯·卡曼常数，取0.41，ρ为泥沙粒子密度；
步骤4.2、基于步骤1.1所得土体粒子粘度μ HBP模型，引入泥沙颗粒Shield准则判断泥沙颗粒起动状态，计算泥沙起动临界应力τ cr,0；
τ bcr＝θ cr·(ρ s-ρ)gd
其中，θ cr为临界希尔兹数,仅与泥沙参数有关，ρ s为饱和泥沙密度，ρ为水体密度，g为重力加速度，d为土体粒子粒径；
步骤4.3、基于步骤4.2所得泥沙起动临界应力τ cr,0，引入步骤2所得泥沙粒子与刚体粒子流-固耦合力F fs,计算考虑外加荷载及边坡效应的修正泥沙起动临界应力τ cr；
Figure PCTCN2022142693-appb-000041
Figure PCTCN2022142693-appb-000042
其中，P为F fs的模，η＝0.7，
Figure PCTCN2022142693-appb-000043
为流速方向单位向量，W为泥沙粒子所受重力，α,β,和γ分别为坡面法向量与X,Y,Z轴方向的夹角，
Figure PCTCN2022142693-appb-000044
w i为
Figure PCTCN2022142693-appb-000045
在X,Y,Z轴方向的分量，
Figure PCTCN2022142693-appb-000046
为重力在坡面方向的分项单位向量，
Figure PCTCN2022142693-appb-000047
为泥沙粒子内摩擦角；
步骤4.4、基于步骤4.1所得流体切应力τ b及步骤4.3所得修正泥沙起动临界应力τ cr判断泥沙粒子起动状态，若满足τ cr≥τ b，则将τ cr代入步骤1.1所得公式替换屈服应力τ cr，更新床沙粒子粘度μ HBP以作为推移质粒子，不满足条件则按步骤1.2至步骤1.3处理；
步骤4.5、根据步骤4.4所得推移质粒子，引入Mastbergen公式计算推移质转化悬 疑质临界流速u lift：
Figure PCTCN2022142693-appb-000048
其中，α i为泥沙输运系数，n s为床面法向量，d *为泥沙粒径系数，
Figure PCTCN2022142693-appb-000049
Figure PCTCN2022142693-appb-000050
ρ s为饱和泥沙密度，ρ为水体密度，g为重力加速度，d为土体粒子粒径，μ为水体粘度，θ b为土体粒子实际希尔兹数，θ cr为临界希尔兹数。
若满足流体粒子实际流速u≥u lift,则转化为悬疑质粒子，并计算其等效粘度μ lift替换μ HBP：
Figure PCTCN2022142693-appb-000051
其中，μ为水体粘度，C v为核函数半径没泥沙粒子浓度。
若不满足条件，则不予处理；
步骤4.6、根据步骤4.5所得悬疑质粒子，引入Mastbergen公式，计算悬疑质转化推移质临界流速u set：
Figure PCTCN2022142693-appb-000052
其中，μ为水体粘度，d为土体粒子粒径，d *为泥沙粒径系数。
若满足粒子实际流速u≤u set,则转化为推移质粒子，并按步骤4.4处理，若不满足条件，则不予处理；
以上具体步骤只是该冲刷模块运行一个时间步长的具体阐述，实际操作中应是此循环的往复进行，直至到达求解终止时间。
以上方法步骤及基本公式原理也可装载到其他开源的SPH计算软件上，使得在计算机或其他可编程设备上执行一系列操作步骤以实现冲刷数值模拟，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程中指定的功能或步骤。
以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (10)

  1. 一种基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真方法，其特征在于，包括如下步骤：
    步骤1，构建粒子模型，基于多相流理论设定流体粒子及床沙粒子，基于DEM理论设定刚体粒子；具体步骤如下：
    步骤1.1、基于多相流基本原理，分别构建流体粒子模型与土体粒子模型，流体粒子模型设定为牛顿流体，土体粒子模型设定为非牛顿流体，基于HBP模型，确定土体粘度μ HBP；
    Figure PCTCN2022142693-appb-100001
    其中，τ c为模型材料屈服应力，II D为流体应变率张量的第二不变量，m为应力指数增长系数，μ为水体粘度，n为与剪切应力相关的幂；
    步骤1.2、引入DP屈服准则，计算具体材料屈服应力τ y：
    步骤1.3、将步骤1.2所得具体材料屈服应力τ y代入步骤1.1所述土体粘度μ HBP计算公式，替换模型材料屈服应力τ c，确立土体粘度计算模型；
    步骤1.4、基于DEM理论设定刚体粒子，确立刚体粒子间法向接触刚度K n、法向接触阻尼γ n、切向接触刚度K t以及切向接触阻尼γ t；
    步骤1.5、基于步骤1.4所得刚体粒子间法向接触刚度K n、法向接触阻尼γ n、切向接触刚度K t以及切向接触阻尼γ t，计算法向接触力F n及切向接触力F t；
    步骤2，基于步骤1所得粒子模型，引入流-固耦合理论修正流体控制方程，进行流体粒子的控制求解，具体步骤为：
    步骤2.1，引入粒子流固耦合力F fs，修正流体动量守恒方程；
    步骤2.2，基于SPH算法的核函数理论，给出步骤2.1所得流体动量守恒方程的离散形式：
    步骤3，基于步骤1所得粒子模型，结合牛顿第二定律，考虑步骤2所得粒子流固耦合力F fs,，修正刚体粒子控制方程并求解：
    步骤4，基于步骤1所得粒子模型，引入泥沙颗粒Shield准则，考虑步骤1所得流-土耦合力、步骤2所得固-土耦合力修正泥沙起动模型，进行基础冲刷控制求解；具体 步骤如下：
    步骤4.1、基于步骤2动量方程获得流体粒子速度u，引入Einstein对数流速分布公式，计算土体粒子所受流体切应力τ b：
    步骤4.2、基于步骤1.1所得土体粒子粘度μ HBP模型，引入泥沙颗粒Shield准则，计算泥沙起动临界应力τ cr,0；
    步骤4.3、基于步骤4.2所得泥沙起动临界应力τ cr,0，引入步骤2所得泥沙粒子与刚体粒子流-固耦合力F fs,计算考虑外加荷载及边坡效应的修正泥沙起动临界应力τ cr；
    步骤4.4、基于步骤4.1所得流体切应力τ b及步骤4.3所得修正泥沙起动临界应力τ cr判断泥沙粒子起动状态，若满足τ cr≥τ b，则将τ cr代入步骤1.1所述公式，替换屈服应力τ cr，更新床沙粒子粘度μ HBP以作为推移质粒子，不满足条件则按步骤1.2至步骤1.3处理；
    步骤4.5、根据步骤4.4所得推移质粒子，引入Mastbergen公式计算推移质转化悬疑质临界流速u lift，若满足流体粒子流速u≥u lift,则转化为悬疑质粒子，并计算等效粘度u lift替换床沙粒子粘度μ HBP，若不满足条件，则不予处理；
    步骤4.6、根据步骤4.5所得悬疑质粒子，引入Mastbergen公式，计算悬疑质转化推移质临界流速u set，若满足粒子实际流速u≤u set，则转化为推移质粒子，并按步骤4.4处理，若不满足条件，则不予处理；
    步骤5，重复步骤1至步骤4，直至求解完成。
  2. 根据权利要求1所述的基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真方法，其特征在于，所述步骤1.2中，材料屈服应力τ y计算如下：
    |τ y|＝αp+β
    其中，p为饱和泥沙粒子上作用的静水压力，α与β均由Mohr-Coulomb屈服准则参数给出，计算公式如下：
    Figure PCTCN2022142693-appb-100002
    其中：θ为内摩擦角，c为土体粘聚力。
  3. 根据权利要求1所述的基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真方法，其特征在于，所述步骤1.4中，刚体粒子间法向接触刚度K n、法向接触阻 尼γ n、切向接触刚度K t以及切向接触阻尼γ t计算如下：
    Figure PCTCN2022142693-appb-100003
    其中，C n＝1*10 -5,K n,ij为DEM刚体粒子i与刚体粒子j之间的法向接触刚度，K t,ij、γ n,ij、γ t,ij同义；E *为等效弹性模量，R *为等效粒子半径，M *为等效质量，计算分别如下：
    Figure PCTCN2022142693-appb-100004
    其中，E i、E j分别为DEM刚体粒子i、刚体粒子j所赋予材料的弹性模量；
    μ i、μ j分别为DEM刚体粒子i、刚体粒子j所赋予材料的泊松比；
    r i、r j分别为DEM刚体粒子i、刚体粒子j的半径；
    m i、m j分别为DEM刚体粒子i、刚体粒子j的质量。
  4. 根据权利要求1所述的一种基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真方法，其特征在于，所述步骤2.1中，流体的动量守恒方程修正为：
    Figure PCTCN2022142693-appb-100005
    式中，u为流体速度矢量，P为流体压力项，υ为流体粘度，g为流体重力项，其余符号同上。
  5. 根据权利要求1所述的基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真方法，其特征在于，所述步骤2.2中，步骤2.1所得控制方程的离散形式为：
    Figure PCTCN2022142693-appb-100006
    其中，
    Figure PCTCN2022142693-appb-100007
    为哈密顿算子，下标a表示中心粒子，下标b表示邻域粒子，W ab为以a粒子为检索中心的核函数
    Figure PCTCN2022142693-appb-100008
    r a和r b分别表示中心粒子和邻域粒子的位置矢量,q为粒子间距与平滑长度的比值,h为光滑长度；
    μ 0为运动粘度；
    τ ij为流体SPS应力矢量；
    m f、P f分别为刚体粒子s在核函数半径内搜索的流体粒子的质量与压力；
    m s、P s为刚体粒子等效质量和压力，
    Figure PCTCN2022142693-appb-100009
    其余变量同上。
  6. 根据权利要求1所述的基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真方法，其特征在于，所述步骤3中，刚体粒子控制方程修正为：
    Figure PCTCN2022142693-appb-100010
    其中，m a、v a分别为刚体粒子a的质量和速度矢量，G为刚体粒子a的重力矢量，其余符号同上。
  7. 根据权利要求1所述的基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真方法，其特征在于，所述步骤4.1中，流体切应力τ b计算为：
    Figure PCTCN2022142693-appb-100011
    其中，d为粒子特征粒径；Δu为泥沙粒子与附近水粒子的速度差，κ为冯·卡曼常数，取0.41，ρ为泥沙粒子密度。
  8. 根据权利要求1所述的基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真方法，其特征在于，所述步骤4.2中，泥沙起动临界应力τ cr,0计算公式为：
    τ cr,0＝θ cr·(ρ s-ρ)gd
    其中，θ cr为临界希尔兹数,仅与泥沙参数有关，ρ s为饱和泥沙密度，ρ为水体密度，g为重力加速度，d为土体粒子粒径。
  9. 根据权利要求1所述的基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真方法，其特征在于，所述步骤4.3中，修正泥沙起动临界应力τ cr计算为：
    Figure PCTCN2022142693-appb-100012
    Figure PCTCN2022142693-appb-100013
    其中，P为F fs的模，η＝0.7，
    Figure PCTCN2022142693-appb-100014
    为流速方向单位向量，W为泥沙粒子所受重力，α,β,和γ分别为坡面法向量与X,Y,Z轴方向的夹角，
    Figure PCTCN2022142693-appb-100015
    w i为
    Figure PCTCN2022142693-appb-100016
    在X,Y,Z轴方向的分量，
    Figure PCTCN2022142693-appb-100017
    为重力在坡面方向的分项单位向量，
    Figure PCTCN2022142693-appb-100018
    为泥沙粒子内摩擦角。
  10. 根据权利要求1所述的基于SPH-DEM耦合与多相流理论的基础冲刷流-固-土耦合仿真方法，其特征在于，所述步骤4.5中，推移质转化悬疑质临界流速u lift计算为：
    Figure PCTCN2022142693-appb-100019
    其中，α i为泥沙输运系数，n s为床面法向量，d *为泥沙粒径系数，
    Figure PCTCN2022142693-appb-100020
    Figure PCTCN2022142693-appb-100021
    ρ s为饱和泥沙密度，ρ为水体密度，g为重力加速度，d为土体粒子粒径，θ b为土体粒子实际希尔兹数，θ cr为临界希尔兹数。
    等效粘度μ lift计算为：
    Figure PCTCN2022142693-appb-100022
    其中，μ为水体粘度，C v为核函数半径内泥沙粒子浓度；
    所述步骤4.6中，临界流速u set计算公式为：
    Figure PCTCN2022142693-appb-100023
    其中，μ为水体粘度，d为土体粒子粒径，d *为泥沙粒径系数。
PCT/CN2022/142693 2022-11-11 2022-12-28 一种基于sph-dem耦合与多相流理论的基础冲刷流-固-土耦合仿真方法 WO2024098534A1 (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN202211417631.4A CN115795986A (zh) 2022-11-11 2022-11-11 一种基于sph-dem耦合与多相流理论的基础冲刷流-固-土耦合仿真方法
CN202211417631.4 2022-11-11

Publications (1)

Publication Number Publication Date
WO2024098534A1 true WO2024098534A1 (zh) 2024-05-16

Family

ID=85437232

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2022/142693 WO2024098534A1 (zh) 2022-11-11 2022-12-28 一种基于sph-dem耦合与多相流理论的基础冲刷流-固-土耦合仿真方法

Country Status (2)

Country Link
CN (1) CN115795986A (zh)
WO (1) WO2024098534A1 (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106055851A (zh) * 2016-07-19 2016-10-26 山西省交通科学研究院 一种土体大变形流动的计算模拟方法
US20180239848A1 (en) * 2017-02-21 2018-08-23 Livermore Software Technology Corporation Numerical Blast Simulation Methods and Systems Thereof
CN112131633A (zh) * 2020-09-04 2020-12-25 山东大学 一种基于粗粒化计算理论的流固耦合仿真方法及系统
CN114781136A (zh) * 2022-04-07 2022-07-22 河海大学 一种基于dem-cfd-sph建模的滑坡涌浪模拟方法
CN114970399A (zh) * 2022-06-09 2022-08-30 东南大学 一种基于sph多相流的桥梁冲刷模块化仿真方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106055851A (zh) * 2016-07-19 2016-10-26 山西省交通科学研究院 一种土体大变形流动的计算模拟方法
US20180239848A1 (en) * 2017-02-21 2018-08-23 Livermore Software Technology Corporation Numerical Blast Simulation Methods and Systems Thereof
CN112131633A (zh) * 2020-09-04 2020-12-25 山东大学 一种基于粗粒化计算理论的流固耦合仿真方法及系统
CN114781136A (zh) * 2022-04-07 2022-07-22 河海大学 一种基于dem-cfd-sph建模的滑坡涌浪模拟方法
CN114970399A (zh) * 2022-06-09 2022-08-30 东南大学 一种基于sph多相流的桥梁冲刷模块化仿真方法

Also Published As

Publication number Publication date
CN115795986A (zh) 2023-03-14

Similar Documents

Publication Publication Date Title
Mukherjee et al. Crown behavior in drop impact on wet walls
Xue et al. Application of a support vector machine for prediction of slope stability
Vaz Jr et al. Considerations on parameter identification and material response for Gurson-type and Lemaitre-type constitutive models
Vijayaraghavan et al. Sustainable manufacturing models for mass finishing process
Martakis et al. Fusing damage-sensitive features and domain adaptation towards robust damage classification in real buildings
CN114970399A (zh) 一种基于sph多相流的桥梁冲刷模块化仿真方法
Zeng et al. An adaptive peridynamics material point method for dynamic fracture problem
WO2024098534A1 (zh) 一种基于sph-dem耦合与多相流理论的基础冲刷流-固-土耦合仿真方法
Mohanty et al. An adaptive neuro fuzzy inference system model for studying free in plane and out of plane vibration behavior of curved beams
CN113486543B (zh) 基于扩展有限元法的坝-基系统地震响应及破坏分析方法
Ning et al. Physics-informed neural network frameworks for crack simulation based on minimized peridynamic potential energy
Albaijan et al. Optimal machine learning-based method for gauging compressive strength of nanosilica-reinforced concrete
Kim et al. Geometric modification for the enhancement of an airfoil performance using deep CNN
Maalawi et al. Design optimization of mechanical elements and structures: A review with application
Cheng et al. A multi-grid sampling multi-scale method for crack initiation and propagation
Qi et al. Study on failure of an ecological tunnel gate caused by jet-flow from air vents
Pogudalina et al. Coupled calculation of the airflow interaction with elastic rod of square cross section
CN109635326A (zh) 机械结构及航空液压管路震动失效灵敏度分析方法
Hidayat et al. Finite element and generalized regression neural network modelling of multiple cracks growth under the influence of multiple crack parameters
CN106529011B (zh) 一种用于sph算法的并行分区实现方法
Wang et al. Improved back propagation neural network based on the enrichment for the crack propagation
Cheng et al. New discrete element models for elastoplastic problems
Zhang et al. An efficient fully Lagrangian solver for modeling wave interaction with oscillating wave energy converter
CN115495854B (zh) 计算机辅助工程模型的参数标定方法、装置、设备及介质
Yang et al. A direct numerical simulation method for solid-solid collision and coupling with fluid