A method for simulating whole process of large deformation of earthquake landslide based on grid isolation technology

Abstract

The application discloses a large deformation whole process simulation method of earthquake landslide based on a grid isolation technology and relates to the technical field of geotechnical engineering.The application constructs a separated background grid and a multi-medium model containing a safety isolation zone based on a second-order convection particle domain interpolation material point method, carries out dynamic analysis on the earthquake landslide by using a strain-driven hardening-softening coupling constitutive mechanism, realizes flexible input of the earthquake wave at the isolation zone and anti-drifting of the base, and can accurately simulate the complete evolution process of the landslide from undrained shear cracking to large deformation friction accumulation.The application effectively eliminates wave reflection interference of the dynamic boundary and numerical drift of long-time integration, realizes real physical reproduction of the whole process of the earthquake landslide from cracking to friction accumulation, makes calculation of large deformation and permanent displacement of the slope under the action of the earthquake more accurate and stable, and lays a foundation for prediction and evaluation of the whole process of the landslide disaster evolution.

Description

Technical Field

[0001] This invention relates to the field of geotechnical engineering technology, specifically to a method for simulating the entire process of large deformation in earthquake landslides based on grid isolation technology. Background Technology

[0002] Existing mesh-based methods (such as the finite element method) often face convergence challenges due to mesh distortion when dealing with large deformation problems in earthquake-induced landslides. While the Material Point Method (MPM) effectively expands the computational capabilities for large deformations, it still faces challenges in refining the dynamic boundary: conventional boundary treatment methods easily introduce non-physical reflection interference between the master model and the free field, while adopting a separate computation strategy increases the complexity of data interaction and time-step coordination. Furthermore, conventional single softening constitutive models have limitations in reproducing the complete evolution of landslides from crack initiation to frictional deposition, easily exhibiting continuous flow-slip characteristics that do not conform to reality, and long-duration dynamic integral calculations are also susceptible to the influence of basement velocity drift. Therefore, proposing an integrated simulation method that can take into account accurate input of dynamic boundaries, suppress numerical drift, and realistically reflect the mechanical behavior of the entire landslide process has both theoretical and practical significance. Summary of the Invention

[0003] This invention provides a simulation method for the entire process of large deformation in earthquake landslides based on grid isolation technology, in order to solve the technical problems of severe free-field boundary wave reflection interference, significant long-duration dynamic integral numerical drift, and difficulty in continuously and accurately simulating the entire process of landslides from brittle initiation to large deformation friction accumulation under the same computational framework when numerical simulation methods deal with earthquake landslide problems.

[0004] According to the first aspect, one embodiment provides a method for simulating the entire process of large deformation in earthquake landslides based on grid isolation technology, the method comprising:

[0005] Based on the slope geometry data, a background grid was established using the second-order convection particle domain interpolation material point method. The computational domain was discretized into free field columns, main slope body and safety isolation zone regions, and different mechanical parameters were assigned to each region to construct a multi-media material point discretization model.

[0006] Based on the multi-medium material point discrete model, quasi-static calculations are carried out by combining gravity loading and local damping to obtain the physical static equilibrium state.

[0007] Seismic acceleration time history records were acquired and integrated, baseline corrected, and stress transformed to construct an equivalent incident shear stress time history suitable for elastic absorbing boundaries.

[0008] Using the static equilibrium state as the initial condition, the explicit dynamic time step loop is entered to calculate particle momentum and nodal forces. Virtual viscous coupling forces are applied on both sides of the isolation zone to achieve flexible connection. A viscous energy-absorbing boundary is applied at the bottom of the model, and drift prevention is implemented at the bottom boundary of the model to complete the seismic dynamic analysis calculation. The equivalent incident shear stress time history is used as the dynamic load condition applied to the energy-absorbing boundary at the bottom of the model.

[0009] Based on the applied boundary conditions, the momentum equation of the background grid nodes is solved to obtain the node acceleration. Then, the node motion is mapped back to the particles, the particle velocity and position are updated, and strain-driven soft-hardening coupled constitutive calculation is performed. The plastic strain increment is calculated according to the stress state and yield criterion of the particles. The key parameters characterizing soil strength, namely soil cohesion and internal friction angle, are updated in real time to simulate the evolution process of landslide from crack initiation to accumulation and obtain the full-time particle position coordinate sequence and particle state variable sequence.

[0010] Based on the obtained full-time particle position coordinate sequence and particle state variable sequence, dynamic video and full-process displacement cloud map containing the spatiotemporal evolution characteristics of landslide are generated, and the shear zone penetration path and final accumulation morphology are analyzed.

[0011] Furthermore, based on the slope geometry data, a background mesh was established using the second-order convection particle domain interpolation material point method. The computational domain was discretized into free-field columns, the main slope body, and the safety isolation zone region. Differentiated mechanical parameters were assigned to each region, and a multi-medium material point discretization model was constructed, specifically including:

[0012] Based on MATLAB software, a regular Euler background mesh is generated in the computational domain. Using the abscissa topological constraint condition, a gap is forced between the free field cylindrical mesh and the main slope mesh to form a physical safety isolation zone.

[0013] The location of the soil-rock interface is determined based on the set thickness of the soil and bedrock layers. Particles with a vertical coordinate below the soil-rock interface are identified using a spatial coordinate mask and defined as bedrock. Particles located above the bedrock and inside the isolation zone are defined as soil. A transition boundary is constructed between the bedrock and soil layers to prevent numerical delamination, so that the bedrock layer only acts as a wave propagation medium in seismic dynamic analysis. Particles outside the isolation zone are identified using a spatial coordinate mask and defined as free-field columns.

[0014] Based on the regional division, different constitutive models and mechanical parameters are assigned to each region.

[0015] Furthermore, based on the regional division, differentiated constitutive models and mechanical parameters are assigned to each region, specifically including:

[0016] The bedrock layer and free-field column are endowed with high elastic modulus and infinite yield strength, and linear elastic constitutive model is performed so that they only act as waveguide medium in seismic dynamic analysis and do not undergo plastic yielding and failure.

[0017] For the soil layer, assign it an initial cohesion peak value c. peak Peak internal friction angle peak With the residual cohesion value c res Residual value of internal friction angle res The soft-hardening coupled elastoplastic constitutive model is implemented to provide physical criteria for simulating the entire evolution of landslides.

[0018] For the transition boundary, a higher cohesion value than that inside the main slope is assigned to enhance the numerical stability of the bottom contact.

[0019] Furthermore, based on a multi-medium material point discrete model, quasi-static calculations are performed using a combination of gravity loading and local damping to obtain the physical static equilibrium state, specifically including:

[0020] Gravitational acceleration is applied to all material points in the field. At the same time, in order to quickly dissipate the unbalanced kinetic energy of the system, inviscid local damping force is applied to the background grid nodes at each time step.

[0021] Continuously perform explicit time integration to monitor the total kinetic energy E of the system in real time. k Compared with historical peak kinetic energy E k,peak The ratio of the system's kinetic energy to its kinetic energy reflects the degree of dissipation of the system's kinetic energy. When the ratio satisfies the preset small-order convergence criterion, i.e., when the ratio is lower than the preset threshold, it indicates that the system's kinetic energy has been basically exhausted, and the system has reached a state of physical static equilibrium.

[0022] After determining that an equilibrium state has been reached, the kinematic state of all particles is reset, that is, the velocity and acceleration are reduced to zero to eliminate the small numerical noise remaining in the quasi-static calculation. At the same time, the current particle position, deformation gradient and Cauchy stress tensor are completely preserved. The corresponding physical quantities constitute the real initial geostress field and will be used as the starting state for seismic dynamic analysis.

[0023] Furthermore, seismic acceleration time history records are acquired and integrated, baseline corrected, and stress transformed to construct an equivalent incident shear stress time history suitable for elastic absorbing boundaries, specifically including:

[0024] Import the target earthquake acceleration time history file from an external seismic motion database, read the externally input earthquake record file, and extract the time series t and acceleration series a(t). Use the trapezoidal numerical integration method to convert the discrete earthquake acceleration series a(t) into the original velocity time history v. raw (t);

[0025] For the original velocity time history v raw (t) Perform linear detrending term processing and cosine window attenuation to eliminate low-frequency drift and force the post-vibration velocity to zero, and obtain the corrected velocity time history v(t);

[0026] Based on the one-dimensional elastic wave propagation theory and the elastic bottom boundary condition, the modified velocity time history v(t) is converted into the equivalent incident shear stress time history τ. in (t).

[0027] Furthermore, using the static equilibrium state as the initial condition, the system enters an explicit dynamic time-step cycle to calculate particle momentum and nodal forces. Virtual viscous coupling forces are applied to both sides of the isolation zone to achieve flexible connections. A viscous energy-absorbing boundary is applied to the bottom of the model, and drift prevention measures are implemented at the bottom boundary. This completes the seismic dynamic analysis calculation. The equivalent incident shear stress time history serves as the dynamic load condition applied to the energy-absorbing boundary at the bottom of the model, specifically including:

[0028] Based on the solution process of the second-order convection particle domain interpolation material point method, the mass and momentum carried by the particles are mapped to the background grid nodes, and the nodal momentum and nodal force are calculated.

[0029] Based on the theory of one-dimensional elastic wave propagation and the principle of wave field superposition, a Lysmer-Kuhlemeyer viscous boundary condition is applied to the bottom of the model, and the calculated incident stress F is superimposed. input , A node The effective force-bearing area represented by the boundary node. The equivalent incident shear stress time history;

[0030] On both sides of the isolation zone, the corresponding free-field particles and main slope particles are identified, their relative velocities Δv are calculated, and a virtual viscous coupling force F is applied. dashpot ; through coupling force F dashpot It transmits the seismic wave response of the free-field column to the main slope body, while absorbing the reflected and scattered waves generated inside the main slope body, and allows the main slope body to slip relative to the free-field column when large deformation occurs, preventing the boundary from being rigidly locked.

[0031] At each time step, the base boundary layer particles with ordinate y less than the threshold δ are dynamically identified, and their vertical velocity and vertical displacement are set to 0 to achieve anti-drift correction.

[0032] Furthermore, based on the applied boundary conditions, the momentum equations of the background mesh nodes are solved to obtain the node accelerations. Then, the node motion is mapped back to the particles, updating the particle velocities and positions. A strain-driven soft-hardening coupled constitutive calculation is performed, calculating the plastic strain increment based on the particle stress state and yield criterion. Key parameters characterizing soil strength, namely soil cohesion and internal friction angle, are updated in real time to simulate the landslide evolution from crack initiation to deposition. The entire time-series particle position coordinate sequence and particle state variable sequence are obtained, specifically including:

[0033] During stress updates, the equivalent plastic strain increment at the current time step is calculated using the backtracking mapping algorithm of the constitutive model. The cumulative equivalent plastic strain of the particles is obtained by time integration: ;

[0034] Based on accumulated equivalent plastic strain p Define the structural damage state of a material and introduce the critical plastic strain. crit As a softening rate control parameter, the damage factor S is calculated using the following formula:

[0035]

[0036] In the formula, S=1.0 represents that the soil structure is intact, and S=0.0 represents that the soil structure is completely destroyed and enters the residual state. p To accumulate equivalent plastic strain, crit The critical plastic strain is set based on the brittle / ductile characteristics of the soil material;

[0037] This simulation of soil strength degradation due to structural loss under strong earthquakes is used to determine whether the soil will yield in the next moment, based on the current soil cohesion c. cur As the destruction factor S decreases, from the initial peak c peak Linear decay to residual value c res c cur The calculation formula is as follows:

[0038]

[0039] In the formula, c peak c is the initial peak value of cohesion. res The set residual cohesion value;

[0040] Establish a soil friction recovery mechanism, including:

[0041] a. Initial stage, i.e., S≈1: Set the initial internal friction angle of the soil layer. peak=0°, simulating a very short earthquake period, saturated soil is in an undrained shear state, initial material parameters: undrained shear strength S u Provided solely by cohesion, the material exhibits undrained shear properties, and the failure envelope degenerates into the Tresca circular yield criterion.

[0042] b. Evolutionary stage, i.e., S→0: As the sliding body undergoes large deformation and failure, the excess pore water pressure dissipates or particles rearrange, and the frictional effect gradually recovers, with the current internal friction angle... cur As S decreases, it dynamically increases to the residual value. res ;

[0043] in, cur The calculation formula is as follows:

[0044]

[0045] In the formula, peak This represents the initial peak value of the internal friction angle. res The set residual value of the internal friction angle;

[0046] By utilizing the high confining pressure of deep soil to generate frictional resistance in the later stages of landslide movement, the kinetic energy dissipation and accumulation evolution process of landslides under natural conditions due to the recovery of frictional intensity can be realistically simulated, thereby accurately predicting the final disaster range.

[0047] Furthermore, based on the obtained full-time particle position coordinate sequence and particle state variable sequence, a dynamic video and full-process displacement cloud map containing the spatiotemporal evolution characteristics of the landslide are generated. The shear zone penetration path and final deposition morphology are analyzed, specifically including:

[0048] Using MATLAB software, the full-time particle position coordinate sequence was read. With the particle state variable sequence, i.e., cumulative plastic strain ;

[0049] Using time t as the axis, the particle positions are plotted frame by frame to generate a dynamic video reflecting the landslide's trajectory; the vector difference between the particle's final position and its initial position is calculated to generate a total displacement cloud map of the entire process.

[0050] Extracting the cumulative plastic strain at the final moment High-value areas are used to identify the penetration path of shear slip zones and to determine the final disaster impact range based on the leading edge coordinates of the accumulation body.

[0051] Furthermore, based on MATLAB software, a regular Eulerian background mesh is generated within the computational domain. Using abscissa topological constraints, a forced gap is created between the free-field cylindrical mesh and the main slope mesh, forming a physical safety isolation zone. Specifically, this includes:

[0052] The horizontal distance between the free-field column mesh and the main slope mesh is set to satisfy:

[0053]

[0054] Among them W gap For the width of the safety barrier, R s The radius of the support domain of the material point shape function is used to physically isolate the interpolation interference of the shape function.

[0055] Furthermore, the method also includes:

[0056] Original velocity time history v raw The formula for calculating (t) is as follows:

[0057]

[0058] In the formula, t i The current time is given, and Δt is the sampling time interval. k Represents the number of moments from the start time to the current time. k Each time step a ( t k ) is the first k The input ground motion acceleration amplitude corresponding to each time step; a ( t k-1 ) is the first k -1 time step (i.e., the previous moment) corresponds to the input ground motion acceleration amplitude;

[0059] The formula for calculating the corrected velocity time history v(t) is as follows:

[0060]

[0061]

[0062] In the formula, A and B are the linear trend coefficients obtained by fitting using the least squares method. A is the slope of the linear trend, representing the linear velocity drift rate caused by constant acceleration error during the integration process when the velocity baseline is corrected; B is the intercept of the linear trend, representing the initial velocity offset constant at the start of the integration. v detrend ( t The time history is the intermediate velocity after removing the linear trend; w(t) is the cosine decay window function, representing the transition at the end of the time history, i.e., t > T.start Defined as For the remaining time periods, w(t) = 1.0. T start This is the starting time of the cosine attenuation window (usually taken at 95% of the total duration of the seismic wave). T end This is the end time of the seismic wave record (i.e., the total duration).

[0063] Equivalent incident shear stress time history τ in The formula for calculating (t) is as follows:

[0064]

[0065] In the formula, ρ base C is the density of bedrock. s,base Let be the shear wave velocity of the bedrock, derived from the formula. Confirmed, G base 1 represents the shear modulus of the bedrock; coefficient 2 is the wave field superposition coefficient of the elastic bottom boundary, used to counteract the attenuation effect of the absorbing boundary on the incident wave.

[0066] This invention provides a method for simulating the entire process of large deformation in earthquake landslides based on grid isolation technology, which has the following advantages: This invention effectively eliminates wave reflection interference from dynamic boundaries and numerical drift from long-duration integrals, and realizes a true physical reproduction of the entire process of earthquake landslides from crack initiation to friction accumulation; it makes the calculation of large deformation and permanent displacement of slopes under seismic action more accurate and stable, laying the foundation for the prediction and evaluation of the entire evolution process of landslide disasters. Attached Figure Description

[0067] Figure 1 A flowchart illustrating a method for simulating the entire process of large deformation in earthquake-induced landslides based on grid isolation technology, as provided in one embodiment of the present invention;

[0068] Figure 2 This is a schematic diagram of the model layout and safety isolation zone in a method for simulating the entire process of large deformation of earthquake landslides based on grid isolation technology, provided in one embodiment of the present invention.

[0069] Figure 3 The image shows the expected effect of simulating plastic strain in a method for simulating the entire process of large deformation in an earthquake landslide based on grid isolation technology, as provided in one embodiment of the present invention. Detailed Implementation

[0070] The present invention will now be described in further detail with reference to specific embodiments and accompanying drawings. Similar elements in different embodiments are referred to by associated similar element reference numerals. In the following embodiments, many details are described to facilitate a better understanding of the invention. However, those skilled in the art will readily recognize that some features may be omitted in different situations, or may be replaced by other elements, materials, or methods. In some cases, certain operations related to the present invention are not shown or described in the specification. This is to avoid obscuring the core parts of the invention with excessive description. For those skilled in the art, detailed description of these related operations is not necessary; they can fully understand the related operations based on the description in the specification and general technical knowledge in the art.

[0071] Furthermore, the features, operations, or characteristics described in the specification can be combined in any suitable manner to form various embodiments. At the same time, the steps or actions in the method description can be rearranged or adjusted in a manner obvious to those skilled in the art. Therefore, the various orders in the specification and drawings are only for the clear description of a particular embodiment and do not imply a necessary order, unless otherwise stated that a particular order must be followed.

[0072] The first embodiment of this invention provides a method for simulating the entire process of large deformation in earthquake landslides based on grid isolation technology. The following is a combination of... Figure 1 Please provide a detailed explanation.

[0073] S1. Constructing and Initializing a Separate Numerical Model: Based on the slope's geometric data, a second-order convection particle domain interpolation material point method is used. A background mesh is established using MATLAB software. The topological relationship between particle coordinates and the preset boundary is determined based on spatial location indexing. The computational domain is discretized into free-field columns and the main slope body (i.e.,...). Figure 2 The model consists of four independent regions: the central master model and the safety isolation zone (Gap). Differentiated mechanical parameters are assigned to each region based on its physical properties. This constructs a multi-medium material point discrete model containing complete material information, which serves as the input model for the initial geostress calculation in step S2.

[0074] Second-order Convected Particle Domain Interpolation (CPDI2) is an advanced interpolation technique used in the Material Point Method (MPM) to address numerical noise and instability issues arising from particle crossings of the background mesh boundary in standard MPM simulations of large deformation problems. The computational process of CPDI2-MPM mainly includes: (a) Particle-corner-mesh mapping: mapping the physical information carried by particles to background mesh nodes using particle-corner characteristic functions and mesh-corner interpolation functions; (b) Mesh node information calculation: solving the momentum equation on the background mesh to update the velocity and internal force state information of the mesh nodes; (c) Physical information mapping back to particles: interpolating and mapping the updated mesh node motion information back to the particles and their corners; (d) Particle and particle domain shape and position update: updating the spatial coordinates of the particles based on the motion information, causing the particle control domain to deform and update its position accordingly.

[0075] Step S1 above specifically includes the following steps:

[0076] S11. Background Mesh Construction and Isolation Zone Definition: A regular Euler background mesh is generated within the computational domain. Using the horizontal coordinate topological constraint, a gap is forced between the free field column mesh and the main slope mesh to form a physical safety isolation zone.

[0077] The horizontal distance between the free-field column mesh and the main slope mesh is set to satisfy:

[0078]

[0079] Among them W gap For the width of the safety barrier, R s The support domain radius of the material point shape function (typically 1.5 to 2 times the mesh size) is used to physically isolate shape function interpolation interference;

[0080] S12. Main Slope Generation and Parameterization of Elastic Bottom Boundary: Based on the established discretized model and preset slope geometric parameters, a spatial coordinate mask is used to automatically identify particles with a ordinate (y) lower than the soil-rock interface, defining them as bedrock layers. The location of the soil-rock interface can be determined according to the user-defined soil and bedrock layer thicknesses. Particles located above the bedrock layer and inside the isolation zone are defined as soil layers. A transition boundary is constructed between the two to prevent numerical delamination, ensuring that the bedrock layer only serves as a wave propagation medium in subsequent dynamic analysis. Figure 2 As shown;

[0081] S13, Free Field Pillar Generation: Using a spatial coordinate mask to identify particles outside the isolation zone, they are defined as free field pillars;

[0082] S14. Assignment of Differentiated Mechanical Parameters: Based on the above regional division, assign differentiated constitutive models and mechanical parameters to each medium:

[0083] a. Assign high elastic modulus and infinite yield strength to the bedrock layer and free field column, and perform linear elastic constitutive modeling so that it only acts as a waveguide medium during subsequent seismic wave input (step S3) and transmission, without plastic yielding and failure;

[0084] b. For the soil layer, assign it an initial peak strength. peak , peak With residual strength c res , res The soft-hardening coupled elastoplastic constitutive model is executed to provide physical criteria for simulating the entire evolution of the landslide in step S5.

[0085] c. Assign a higher cohesion value to the transition boundary than to the main slope body to enhance the numerical stability of the bottom contact.

[0086] S2. Static Equilibrium and Geostress Initialization: On the discrete model established in step S1, quasi-static calculations are performed using a combination of gravity loading and local damping. This step aims to eliminate non-physical kinetic energy oscillations caused by sudden gravity loading, calculate the initial geostress field and deformation state of the slope under its own weight, and use the particle stress tensor and position coordinates in this equilibrium state as the initial conditions for the dynamic analysis in steps S4 and S5.

[0087] Step S2 above specifically includes the following steps:

[0088] S21. Gravity Loading and Damping Application: A gravitational acceleration g is applied to all material points in the field. Simultaneously, to rapidly dissipate the unbalanced kinetic energy of the system, a non-viscous local damping force F is applied to the background mesh nodes at each time step. damp :

[0089]

[0090] In the formula: α is the local damping coefficient, which is usually selected from an empirical range of 0.5 to 0.8 to achieve a balance between computational convergence rate and numerical stability; |F unbal| represents the magnitude of the resultant force of the external load, internal stress, and gravity currently acting on the node, and represents the remaining driving force driving the node's motion; sgn(v) is the velocity sign function, used to extract the direction of the node's instantaneous velocity (takes 1 when v>0, and -1 when v<0); through this formula, the damping force is always applied in the opposite direction of the node's motion, and its magnitude is proportional to the current unbalanced force, thereby efficiently attenuating the system oscillation until it comes to rest.

[0091] S22. Quasi-static convergence criterion: Continuously perform explicit time integration and monitor the total kinetic energy E of the system in real time. k Compared with historical peak kinetic energy E k,peak The ratio of kinetic energy to kinetic energy. This ratio reflects the degree of dissipation of the system's kinetic energy, and it is determined when it meets a preset small-order convergence criterion (E). k / E k,peak <1.0×10 -5 When the system reaches a state of static equilibrium, it indicates that the system's kinetic energy has been largely exhausted, and the system has reached a state of physical static equilibrium.

[0092] S23. Initial dynamic values: After reaching the equilibrium criterion, the kinematic state of all particles is reset (i.e., velocity and acceleration are reduced to zero) to eliminate the small numerical noise remaining from the quasi-static calculation, while fully preserving the current particle positions, deformation gradients, and Cauchy stress tensors. These physical quantities constitute the true initial geostress field, which will be directly passed to steps S3 and S4 as the starting state for the seismic dynamic analysis (t=0).

[0093] S3. Seismic wave conversion and stress time history construction: Import the target seismic acceleration time history record file from the external seismic motion database, read the seismic acceleration time history record, integrate it, correct the baseline and convert the stress to construct the equivalent incident shear stress time history suitable for the elastic absorbing boundary. The time history will be used as the dynamic load condition applied to the energy-absorbing boundary at the bottom of the model in step S4 to realize the input of seismic energy into the computational domain.

[0094] Step S3 above specifically includes the following steps:

[0095] S31. Data Import and Integration Conversion: Read the externally input seismic record file and extract the time series t and acceleration series a(t). Use the trapezoidal numerical integration method to convert the discrete seismic acceleration time history record a(t) into the original velocity time history v. raw (t);

[0096] v raw The formula for calculating (t) is as follows:

[0097]

[0098] In the formula, t i The current time is given, and Δt is the sampling time interval.k Represents the number of moments from the start time to the current time. k Each time step a ( t k ) is the first k The input ground motion acceleration amplitude corresponding to each time step; a ( t k-1 ) is the first k -1 time step (i.e., the previous moment) corresponds to the input ground motion acceleration amplitude;

[0099] S32, Baseline Correction: Adjusting the original velocity time history v raw (t) Perform linear detrend processing and cosine taper attenuation to eliminate low-frequency drift and force the post-oscillation velocity to zero, and obtain the corrected velocity time history v(t);

[0100] The formula for calculating v(t) is as follows:

[0101]

[0102]

[0103] In the formula, A and B are the linear trend coefficients obtained by fitting using the least squares method. A is the slope of the linear trend, representing the linear velocity drift rate caused by constant acceleration error during the integration process when the velocity baseline is corrected; B is the intercept of the linear trend, representing the initial velocity offset constant at the start of the integration. v detrend ( t () represents the intermediate velocity time history after removing the linear trend. w (t) is a cosine decay window function, in the transition segment at the end of the time history (t>T). start ) is defined as:

[0104]

[0105] For the remaining time periods, w(t) = 1.0; T start This is the starting time of the cosine attenuation window (usually taken at 95% of the total duration of the seismic wave). T end This is the end time of the seismic wave record (i.e., the total duration).

[0106] S33. Stress Transformation: Based on the one-dimensional elastic wave propagation theory and elastic bottom boundary conditions, the modified velocity time history v(t) is converted into the equivalent incident shear stress time history τ. in (t);

[0107] τin The formula for calculating (t) is as follows:

[0108]

[0109] In the formula, ρ base C is the density of bedrock. s,base Let be the shear wave velocity of the bedrock, derived from the formula. Confirmed, G base 1 represents the shear modulus of the bedrock; coefficient 2 is the wave field superposition coefficient of the elastic bottom boundary, used to counteract the attenuation effect of the absorbing boundary on the incident wave.

[0110] S4. Dynamic loading and boundary condition application: Using the equilibrium state determined in step S2 as the initial condition, enter the explicit dynamic time step loop, calculate particle momentum and nodal force, apply virtual viscous coupling force on both sides of the isolation zone to achieve flexible connection, apply viscous energy-absorbing boundary at the bottom of the model, and prevent drift at the bottom boundary of the model to achieve lossless input of seismic waves and absorption of reflected waves, and realize dynamic analysis calculation.

[0111] The specific steps of step S4 are as follows:

[0112] S41. Entering the dynamic cycle and momentum mapping: According to the CPDI2-MPM solution process, the mass m carried by the particle is... p With momentum m p v p Mapping to background mesh nodes, calculating nodal momentum and nodal forces provides a basis for subsequent deformation and position updates; firstly, based on the particle deformation gradient F... p Update the spatial coordinates of the four corner points of the particle to construct a quadrilateral particle domain that deforms with the material; based on the CPDI2 shape function N i (x) precisely integrates the physical quantities carried by the particle and maps them to the background mesh nodes:

[0113] a. Nodal momentum mapping (p i ):

[0114]

[0115] In the formula, i is the grid node index, p is the particle index, and m p v p These represent particle mass and velocity, respectively; Ni(x) p ) is a CPDI2 shape function, the value of which depends on the topological relationship between the particle corner point and the grid node.

[0116] b. Calculation of nodal internal forces:

[0117]

[0118] In the formula, V pLet σ be the current volume of the particle. p For the particle Cauchy stress tensor N i This is the gradient of the shape function. In CPDI2, N i Linear interpolation calculations are performed at the four corner points throughout the entire particle domain to more accurately capture the stress gradient distribution during large deformation processes.

[0119] S42. Bottom Dynamic Input and Energy Absorbing Boundary Application: Based on the one-dimensional elastic wave propagation theory and the principle of wave field superposition, a Lysmer-Kuhlemeyer (LK) viscous boundary condition is applied, and the incident stress calculated in step S3 is superimposed. Specifically, a boundary nodal force F is applied at the bottom node. base accomplish:

[0120]

[0121]

[0122]

[0123] In the formula, Density of bedrock For bedrock shear wave velocity, Let Anode be the instantaneous velocity of a node, and let Anode be the effective force-bearing area represented by the boundary node. For the equivalent incident shear stress time history, F viscous The design aims to simulate the radiation boundary effect of a semi-infinite space, absorbing reflected waves propagating downwards from the interior of the model. input Overcoming the damping effect of viscous boundaries enables the effective input of seismic energy into the computational domain.

[0124] S43. Flexible Coupling Application: On both sides of the isolation zone, identify the corresponding free-field particles and main slope particles, calculate their relative velocity Δv, and apply a virtual viscous coupling force F. dashpot ;

[0125] F dashpot The calculation formula is as follows:

[0126]

[0127] In the formula, C side Δv is the coupling damping coefficient, and Δv is the velocity difference between the corresponding particles on both sides of the isolation zone. Through this coupling force, the seismic wave response of the free field column is transmitted to the main slope body, while absorbing the reflected and scattered waves generated inside the main slope body. It also allows the main slope body to slip relative to the free field column when large deformation occurs, preventing the boundary from being rigidly locked.

[0128] S44, Anti-drift correction: At each time step, dynamically identify base boundary layer particles whose ordinate y is less than the threshold δ (typically 0.5 times the background mesh size) and perform constraints;

[0129] The constraint formulas are as follows:

[0130]

[0131] In the formula, v y Let u be the vertical velocity. y This is a vertical displacement, used to eliminate base floating caused by numerical integration errors.

[0132] S5. Mesh Solving and Hybrid Constitutive Evolution: Based on the dynamic boundary conditions applied in step S4, solve the momentum equations of the background mesh nodes to obtain the node accelerations; then map the node motion back to the particles to update the particle velocity and position (v). p , x p The strain-driven soft-hardening coupled constitutive calculation is performed, and the plastic strain increment is calculated based on the stress state of the particles and the yield criterion, updating the key parameters characterizing soil strength (i.e., cohesion c and internal friction angle) in real time. (This simulates the evolution of a landslide from initiation to deposition.)

[0133] Step S5 above specifically includes the following steps:

[0134] S51. Calculation of Cumulative Plastic Strain and Plastic Damage Factor: During stress update, the equivalent plastic strain increment at the current time step is calculated using the backtracking mapping algorithm of the constitutive model. The cumulative equivalent plastic strain of the particles is obtained by time integration. Based on cumulative equivalent plastic strain p Define the structural damage state of a material and introduce the critical plastic strain. crit As a softening rate control parameter, calculate the damage factor S and proceed to step S52;

[0135] The formula for calculating the destructive factor S is as follows:

[0136]

[0137] In the formula, S=1.0 represents that the soil structure is intact, and S=0.0 represents that the soil structure is completely destroyed and enters the residual state. p To accumulate equivalent plastic strain, crit This is the critical plastic strain (a fixed value input by the user based on the brittle / ductile characteristics of the soil material).

[0138] S52. Soil Cohesion Softening Mechanism: Simulates the strength attenuation caused by the loss of soil structure under strong earthquakes, used to determine whether the soil will yield in the next moment. Current soil cohesion c cur As the destruction factor S decreases, from the initial peak c peak Linear decay to residual value c res Proceed to step S53;

[0139] c cur The calculation formula is as follows:

[0140]

[0141] In the formula, c peak c is the initial peak value of cohesion. res This is the residual cohesion value (the input material parameter).

[0142] S53. Soil friction recovery mechanism:

[0143] a. Initial stage (S≈1): Set the initial peak value of the internal friction angle of the soil layer. peak =0°, simulating a very short time during an earthquake, the saturated soil is in an undrained shear state, and the shear strength (i.e., the undrained shear strength S of the initial material parameters) is... u Provided solely by cohesion, the material exhibits undrained shear properties, and the failure envelope degenerates into the Tresca circular arc yield criterion.

[0144] b. Evolutionary Stage (S→0): As the sliding body undergoes large deformation and failure, the excess pore water pressure dissipates or particles rearrange, and the frictional effect gradually recovers. Current internal friction angle cur As S decreases, it dynamically increases to the residual value. res .

[0145] cur The calculation formula is as follows:

[0146]

[0147] In the formula, peak This represents the initial peak value of the internal friction angle. res This is the residual value of the internal friction angle (the input material parameter).

[0148] By using the strain-driven soft-hardening coupled constitutive mechanism described in steps S51, S52, and S53, frictional resistance is generated by the high confining pressure of the deep soil in the later stage of landslide movement. This realistically simulates the kinetic energy dissipation and accumulation evolution process of landslides under natural conditions due to the recovery of frictional intensity, thereby accurately predicting the final disaster range.

[0149] S6. Post-processing and output: The output includes dynamic video and full-process displacement cloud map containing the spatiotemporal evolution characteristics of the landslide, and analyzes the shear zone penetration path and final deposition morphology.

[0150] The specific steps for S6 are as follows:

[0151] S61. Using MATLAB software, read the full-time particle position coordinate sequence output in step S5. With particle state variable sequence (cumulative plastic strain) Using time t as the axis, particle positions are plotted frame by frame to generate a dynamic video reflecting the landslide's trajectory; the vector difference between the particle's final and initial positions is calculated to generate a total displacement cloud map of the entire process; and the accumulated plastic strain at the final moment is extracted. High-value areas are used to identify the penetration path of shear slip zones, and the final hazard impact range is determined based on the leading edge coordinates of the accumulation. (The expected effect of simulated plastic strain can be seen in...) Figure 3 ).

[0152] The above examples illustrate the present invention only to aid in understanding it and are not intended to limit the scope of the invention. Those skilled in the art can make various simple deductions, modifications, or substitutions based on the principles of this invention.

Claims

1. A method for simulating the entire process of large deformation in earthquake-induced landslides based on grid isolation technology, characterized in that, The method includes: Based on the slope geometry data, a background grid was established using the second-order convection particle domain interpolation material point method. The computational domain was discretized into free field columns, main slope body and safety isolation zone regions, and different mechanical parameters were assigned to each region to construct a multi-media material point discretization model. Based on the multi-medium material point discrete model, quasi-static calculations are carried out by combining gravity loading and local damping to obtain the physical static equilibrium state. Seismic acceleration time history records were acquired and integrated, baseline corrected, and stress transformed to construct an equivalent incident shear stress time history suitable for elastic absorbing boundaries. Using the static equilibrium state as the initial condition, the explicit dynamic time step loop is entered to calculate particle momentum and nodal forces. Virtual viscous coupling forces are applied on both sides of the isolation zone to achieve flexible connection. A viscous energy-absorbing boundary is applied at the bottom of the model, and drift prevention is implemented at the bottom boundary of the model to complete the seismic dynamic analysis calculation. The equivalent incident shear stress time history is used as the dynamic load condition applied to the energy-absorbing boundary at the bottom of the model. Based on the applied boundary conditions, the momentum equation of the background grid nodes is solved to obtain the node acceleration. Then, the node motion is mapped back to the particles, the particle velocity and position are updated, and strain-driven soft-hardening coupled constitutive calculation is performed. The plastic strain increment is calculated according to the stress state and yield criterion of the particles. The key parameters characterizing soil strength, namely soil cohesion and internal friction angle, are updated in real time to simulate the evolution process of landslide from crack initiation to accumulation and obtain the full-time particle position coordinate sequence and particle state variable sequence. Based on the obtained full-time particle position coordinate sequence and particle state variable sequence, dynamic video and full-process displacement cloud map containing the spatiotemporal evolution characteristics of landslide are generated, and the shear zone penetration path and final accumulation morphology are analyzed.

2. The method for simulating the entire process of large deformation in earthquake landslides based on grid isolation technology as described in claim 1, characterized in that, Based on the slope geometry data, a background mesh was established using the second-order convection particle domain interpolation material point method. The computational domain was discretized into free-field columns, the main slope body, and the safety isolation zone region. Differentiated mechanical parameters were assigned to each region to construct a multi-medium material point discretization model, specifically including: Based on MATLAB software, a regular Euler background mesh is generated in the computational domain. Using the abscissa topological constraint condition, a gap is forced between the free field cylindrical mesh and the main slope mesh to form a physical safety isolation zone. The location of the soil-rock interface is determined based on the set thickness of the soil and bedrock layers. Particles with a vertical coordinate below the soil-rock interface are identified using a spatial coordinate mask and defined as bedrock. Particles located above the bedrock and inside the isolation zone are defined as soil. A transition boundary is constructed between the bedrock and soil layers to prevent numerical delamination, so that the bedrock layer only acts as a wave propagation medium in seismic dynamic analysis. Particles outside the isolation zone are identified using a spatial coordinate mask and defined as free-field columns. Based on the regional division, different constitutive models and mechanical parameters are assigned to each region.

3. The method for simulating the entire process of large deformation in earthquake landslides based on grid isolation technology as described in claim 2, characterized in that, Based on the regional division, differentiated constitutive models and mechanical parameters are assigned to each region, specifically including: For the bedrock layer and free-field column, a high elastic modulus exceeding the preset threshold and an infinite yield strength are given, and a linear elastic constitutive model is performed so that it only acts as a waveguide medium in seismic dynamic analysis and does not undergo plastic yielding and failure. For the soil layer, assign it an initial cohesion peak value c. peak Peak internal friction angle peak With the residual cohesion value c res Residual value of internal friction angle res The soft-hardening coupled elastoplastic constitutive model is implemented to provide physical criteria for simulating the entire evolution of landslides. The transition boundary is given a higher cohesion value than the interior of the main slope to enhance the numerical stability of the bottom contact.

4. The method for simulating the entire process of large deformation in earthquake landslides based on grid isolation technology as described in claim 1, characterized in that, Based on a multi-medium material point discrete model, quasi-static calculations are performed using a combination of gravity loading and local damping to obtain the physical static equilibrium state, specifically including: Gravitational acceleration is applied to all material points in the field, while non-viscous local damping forces are applied to the background mesh nodes at each time step. Continuously perform explicit time integration to monitor the total kinetic energy E of the system in real time. k Compared with historical peak kinetic energy E k,peak The ratio of the system's kinetic energy to its kinetic energy reflects the degree of dissipation of the system's kinetic energy. When the ratio meets the preset small-scale convergence criterion, i.e., the ratio is lower than the preset threshold, it indicates that the system's kinetic energy has been exhausted, and the system has reached a state of physical static equilibrium. After determining that a static equilibrium state has been reached, the kinematic state of all particles is reset, that is, the velocity and acceleration are reduced to zero to eliminate the numerical noise remaining from the quasi-static calculation. At the same time, the current particle position, deformation gradient and Cauchy stress tensor are fully preserved. The corresponding physical quantities constitute the real initial geostress field and will serve as the starting state for seismic dynamic analysis.

5. The method for simulating the entire process of large deformation in earthquake landslides based on grid isolation technology as described in claim 1, characterized in that, Seismic acceleration time history records were acquired and integrated, baseline corrected, and stress transformed to construct an equivalent incident shear stress time history suitable for elastically absorbing boundaries. Specifically, this included: Import the target earthquake acceleration time history file from an external seismic motion database. Read the externally input earthquake acceleration time history file and extract the time series t and the earthquake acceleration sequence a(t). Use the trapezoidal numerical integration method to convert the discrete earthquake acceleration sequence a(t) into the original velocity time history v. raw (t); For the original velocity time history v raw (t) Perform linear detrending term processing and cosine window decay to obtain the corrected velocity time history v(t); Based on the one-dimensional elastic wave propagation theory and the elastic bottom boundary condition, the modified velocity time history v(t) is converted into the equivalent incident shear stress time history τ. in (t).

6. The method for simulating the entire process of large deformation in earthquake landslides based on grid isolation technology as described in claim 1, characterized in that, Using static equilibrium as the initial condition, the system enters an explicit dynamic time-step cycle to calculate particle momentum and nodal forces. Virtual viscous coupling forces are applied to both sides of the isolation zone to achieve flexible connections. A viscous energy-absorbing boundary is applied to the bottom of the model, and drift prevention measures are implemented at the bottom boundary. This completes the seismic dynamic analysis calculation. The equivalent incident shear stress time history serves as the dynamic load condition applied to the energy-absorbing boundary at the bottom of the model, specifically including: Based on the solution process of the second-order convection particle domain interpolation material point method, the mass and momentum carried by the particles are mapped to the background grid nodes, and the nodal momentum and nodal force are calculated. Based on the theory of one-dimensional elastic wave propagation and the principle of wave field superposition, a Lysmer-Kuhlemeyer viscous boundary condition is applied to the bottom of the model, and an incident stress F is superimposed. input ,in A node The effective force-bearing area represented by the boundary node. The equivalent incident shear stress time history; On both sides of the isolation zone, the corresponding free-field particles and main slope particles are identified, their relative velocities Δv are calculated, and a virtual viscous coupling force F is applied. dashpot ; through coupling force F dashpot It transmits the seismic wave response of the free-field column to the main slope body, while absorbing the reflected and scattered waves generated inside the main slope body, and allows the main slope body to slip relative to the free-field column when large deformation occurs, preventing the boundary from being rigidly locked. At each time step, the base boundary layer particles with ordinate y less than the threshold δ are dynamically identified, and their vertical velocity and vertical displacement are set to 0 to achieve anti-drift correction.

7. The method for simulating the entire process of large deformation in earthquake landslides based on grid isolation technology as described in claim 1, characterized in that, Based on the applied boundary conditions, the momentum equations of the background mesh nodes are solved to obtain the node accelerations. The node motion is then mapped back to the particles, updating their velocities and positions. A strain-driven soft-hardening coupled constitutive calculation is performed, calculating the plastic strain increment based on the particle stress state and yield criterion. Key parameters characterizing soil strength, namely soil cohesion and internal friction angle, are updated in real time. The evolution of a landslide from crack initiation to deposition is simulated, obtaining the full-time particle position coordinate sequence and particle state variable sequence, specifically including: During stress updates, the equivalent plastic strain increment at the current time step is calculated using the backtracking mapping algorithm of the constitutive model. The cumulative equivalent plastic strain of the particles is obtained by time integration: ; Based on accumulated equivalent plastic strain p Define the structural damage state of a material and introduce the critical plastic strain. crit As a softening rate control parameter, the damage factor S is calculated using the following formula: In the formula, S=1.0 represents that the soil structure is intact, and S=0.0 represents that the soil structure is completely destroyed and enters the residual state. p To accumulate equivalent plastic strain, crit The critical plastic strain is set based on the brittle / ductile characteristics of the soil material; This simulation of soil strength degradation due to structural loss under strong earthquakes is used to determine whether the soil will yield in the next moment, based on the current soil cohesion c. cur As the destruction factor S decreases, from the initial peak c peak Linear decay to residual value c res c cur The calculation formula is as follows: In the formula, c peak c is the initial peak value of cohesion. res The set residual cohesion value; Establish a soil friction recovery mechanism, including: a. Initial stage: Set the peak value of the initial soil layer's internal friction angle. peak =0°, simulating a very short earthquake period, saturated soil is in an undrained shear state, initial material parameters: undrained shear strength S u Provided solely by cohesion, the material exhibits undrained shear properties, and the failure envelope degenerates into the Tresca circular yield criterion. b. Evolutionary Stage: As the sliding body undergoes large deformation and failure, the excess pore water pressure dissipates or particles rearrange, and the frictional effect gradually recovers, with the current internal friction angle... cur As S decreases, it dynamically increases to the residual value. res ; in, cur The calculation formula is as follows: In the formula, peak This represents the initial peak value of the internal friction angle. res The set residual value of the internal friction angle; By utilizing the high confining pressure of deep soil to generate frictional resistance in the later stages of landslide movement, the kinetic energy dissipation and accumulation evolution process of landslides under natural conditions due to the recovery of frictional intensity can be realistically simulated, thereby accurately predicting the final disaster range.

8. The method for simulating the entire process of large deformation in earthquake landslides based on grid isolation technology as described in claim 7, characterized in that, Based on the obtained full-time particle position coordinate sequence and particle state variable sequence, a dynamic video and full-process displacement cloud map containing the spatiotemporal evolution characteristics of the landslide are generated. The shear zone penetration path and final deposition morphology are analyzed, specifically including: Using MATLAB software, the full-time particle position coordinate sequence was read. With the particle state variable sequence, i.e., cumulative plastic strain ; Using time t as the axis, the particle positions are plotted frame by frame to generate a dynamic video reflecting the landslide's trajectory; the vector difference between the particle's final position and its initial position is calculated to generate a total displacement cloud map of the entire process. Extracting the cumulative plastic strain at the final moment High-value areas exceeding preset thresholds are used to identify the penetration path of shear slip zones, and the final disaster impact range is determined based on the leading edge coordinates of the accumulation body.

9. The method for simulating the entire process of large deformation in earthquake landslides based on grid isolation technology as described in claim 2, characterized in that, Based on MATLAB software, a regular Eulerian background mesh is generated within the computational domain. Using abscissa topological constraints, a forced gap is created between the free-field cylindrical mesh and the main slope mesh, forming a physical safety isolation zone. Specifically, this includes: The horizontal distance between the free-field column mesh and the main slope mesh is set to satisfy: Among them W gap For the width of the safety barrier, R s The radius of the support domain of the material point shape function is used to physically isolate the interpolation interference of the shape function.

10. The method for simulating the entire process of large deformation in earthquake landslides based on grid isolation technology as described in claim 5, characterized in that, The method further includes: Original velocity time history v raw The formula for calculating (t) is as follows: In the formula, t i The current time is given, and Δt is the sampling time interval. k Represents the number of moments from the start time to the current time. k Each time step a ( t k) For the first k The input ground motion acceleration amplitude corresponding to each time step; a ( t k-1 ) is the first k -1 time step, which is the input ground motion acceleration amplitude corresponding to the previous moment; The formula for calculating the corrected velocity time history v(t) is as follows: In the formula, A and B are the linear trend coefficients obtained by fitting using the least squares method. A is the slope of the linear trend, representing the linear velocity drift rate caused by constant acceleration error during the integration process when the velocity baseline is corrected; B is the intercept of the linear trend, representing the initial velocity offset constant at the start of the integration. v detrend ( t The time history is the intermediate velocity after removing the linear trend; w(t) is the cosine decay window function, which is applied during the transition period at the end of the time history, i.e., t > T. start When, defined as For the remaining time periods, w(t) = 1.0, where, T start This is the start time of the cosine decay window; T end This is the end time of the seismic wave record; Equivalent incident shear stress time history τ in The formula for calculating (t) is as follows: In the formula, ρ base C is the density of bedrock. s,base 1 represents the shear wave velocity of the bedrock; coefficient 2 is the wave field superposition coefficient of the elastic bottom boundary, used to counteract the attenuation effect of the absorbing boundary on the incident wave.

Images