Numerical simulation method for simulating influence of root-soil complex on slope stability by discrete element software

Abstract

A numerical simulation method for simulating influence of a root-soil complex on slope stability by discrete element software belongs to the field of data processing and comprises the following steps: acquiring slope soil parameter data; establishing a plain soil model by adopting discrete element software according to the field construction data and the high-filled roadbed slope model; utilizing CAD and PFC2D to simplify a plant root system, and obtaining a plant root system simplified model and a root-soil complex slope model; by applying the same gravity load and operating to a system balance state, analyzing the soil particle displacement condition in a plain soil and root soil complex state; analyzing a force-displacement law of the soil particles; and analyzing the motion law of the soil particles. The actual condition of the slope can be effectively restored by simulating and simplifying the plant root system, operation is easy and convenient, calculation efficiency is high, and practical significance and application value are achieved.

Description

Technical Field

[0001] The invention belongs to the technical field of data processing, and in particular relates to a numerical simulation method for simulating the influence of a root-soil complex on slope stability using discrete element software. Background Art

[0002] A large number of artificial slopes and natural disasters can cause serious soil erosion and landslides. Plant roots have a reinforcing effect on unstable slopes, and ecological slope reinforcement projects can reduce the losses caused by related natural risks. In order to gain a deeper understanding of the reinforcing effect of roots on soil, many scholars have continued to deepen their research on the reinforcement mechanism of root-soil complexes, not only through experiments but also from the perspective of numerical simulation to study the influence of roots on soil properties. Cheng Lei (Cheng Lei, Zhang Shengxu, Hao Yanzhou, et al. Experimental study on the shear strength of root-soil complexes on pile slopes [J]. Science, Technology and Engineering, 2017, 17(23): 10-13.) et al. used the GDS multi-field coupled triaxial test system to study the failure mode of root-soil complexes under different vegetation coverage conditions and found that the contribution coefficient of root tensile strength to soil shear strength was between 0.15-0.33. Huang Xiaole's team (Huang Xiaole, Xu Wennian, Xia Zhenyao. Direct shear test of root-soil complexes of two herbaceous plants on vegetated concrete substrate [J]. Research on Soil and Water Conservation, 2010, 17(4): 51-53.) compared the root reinforcement effect of Cynodon dactylon and Medicago sativa and pointed out that the shallow root system (0-30 cm) contributed 62% to the enhancement of slope stability, while the deep root system (>60 cm) mainly played an anchoring role. Jiao Zhen (Jiao Zhen, Wang Daojie, Xie Hong, et al. Experimental analysis of shear strength of undisturbed soil in the Leucaena leucaena forest area in Jiangjiagou, Dongchuan, Yunnan [J]. China Soil and Water Conservation, 2019, 8(11): 42-44.) et al. established a synergistic model of root content (0.5%-3.5%), soil depth (0-2 m) and heterogeneity coefficient based on multiple regression analysis, revealing that the root reinforcement efficiency shows an exponential decay law with the deterioration of soil particle grading. Wang Yuanzhan et al. (Wang Yuanzhan, Liu Xufei, Zhang Zhikai, et al. Experimental study on the effect of root content on the strength of original and reshaped grass root reinforced soil [J]. Chinese Journal of Geotechnical Engineering, 2015, 37(8): 61-63.) conducted comparative tests on reshaped soil and original soil and found that the peak shear strength attenuation rate of grass root reinforced soil could reach 23%-35%. They also constructed a quantitative model of mechanical parameter differences considering the root orientation characteristics. Huang Gang's team (Huang Gang, Peng Jing, Xiao Henglin. Effect of root content on soil shear strength [J]. China Soil and Water Conservation, 2016, 6(5): 42-44.) measured the critical shear strength of different plant roots (0.18 MPa for bermudagrass and 0.25 MPa for alfalfa) using a large direct shear device and established a slope strength prediction formula considering the root diameter effect. Zhang Xiaoming et al. (Zhang Xiaoming, Wang Yujie, Xia Yiping, et al. Study on the shear strength of undisturbed soil and remolded soil of typical vegetation in Jinyun Mountain, Chongqing [J]. Transactions of the Chinese Society of Agricultural Engineering, 2016, 22(011): 6-9.) compared and analyzed the root-soil complex characteristics of five vegetation types in Jinyun Mountain and found that the shear strength of undisturbed soil was 19%-27% higher than that of remolded soil, revealing the inhibitory effect of construction disturbance on the root reinforcement effect.The above studies are all analyzed from a macroscopic perspective, but lack the study of soil and plant roots from a microscopic perspective. By modeling the plant roots and soil particles, the impact of various influencing factors on slope stability can be judged.

[0003] Although many geotechnical engineering researchers have studied the mechanical effects of underground roots on soil reinforcement through laboratory and field experiments, numerical simulations, and theoretical analysis, the study of the soil reinforcement effect of plant roots remains subject to numerous uncertainties and difficulties. Plant roots are buried in the soil, making them difficult to observe directly. Furthermore, their distribution is random and subject to environmental influences such as climate and soil mass. Root morphology also changes over time, making research challenging. Existing studies often overlook the dynamic growth of roots, resulting in a lack of temporal understanding of the mechanical properties of the root-soil complex. This oversimplifies root morphology and leads to unreliable results. These factors have led to a significant discrepancy between the understanding of the mechanical mechanisms of root reinforcement and the actual situation. This mechanism of root reinforcement is the mechanical foundation of ecological slope protection and a crucial theoretical basis for selecting slope protection plants and designing ecological slope protection. Therefore, it is crucial to further investigate the root reinforcement mechanisms, taking into account the morphology and hierarchical structure of the root system.

[0004] Therefore, it is necessary to establish an accurate and reasonable root-soil composite slope model. In the process of model establishment, the simulation method of the interaction between the root and the soil is particularly important. Summary of the Invention

[0005] The purpose of the present invention is to provide a numerical simulation method for simulating the influence of root-soil complex on slope stability by using discrete element software, which is used for slope stability analysis.

[0006] The technical solutions adopted by the present invention to solve the technical problems are as follows:

[0007] The present invention provides a numerical simulation method for simulating the influence of root-soil complex on slope stability using discrete element software, comprising the following steps:

[0008] (1) Obtain slope soil parameter data;

[0009] (2) Based on the on-site construction data and the high fill roadbed slope model, a soil model was established using discrete element software;

[0010] (3) Use CAD and PFC2D to simplify the plant root system and obtain a simplified plant root system model and a root-soil complex slope model;

[0011] (4) By applying the same gravity load and running the system to the equilibrium state, the displacement of soil particles in the state of bare soil and root-soil complex is analyzed;

[0012] (5) Force-displacement law analysis of soil particles;

[0013] (6) Analysis of the motion laws of soil particles.

[0014] Furthermore, in step (2), the slope ratio of the bare soil model is 1:1.5.

[0015] Furthermore, in step (2), the bare soil model adopts fixed boundary conditions, that is, no constraints are imposed on the upper boundary, horizontal constraints are imposed on the lower boundary while vertical constraints are imposed, and only horizontal constraints are imposed on the left, right, front and rear boundaries.

[0016] Furthermore, in step (3), the required shape diagram is first drawn in the CAD software, and then imported into the particle flow analysis program PFC2D to form a strip unit, and then the strip unit is converted into a super particle cluster, and the super particle cluster is imported into the plain soil model to generate a root-soil complex slope model.

[0017] Furthermore, in step (4), in the bare soil state, when the roadbed slope moves to the system equilibrium state under the action of gravity load, some soil particles are displaced. The displacement of soil particles is mainly concentrated at the roadbed slope, and the maximum displacement of soil particles is 0.53m.

[0018] Furthermore, in step (4), in the root-soil complex state, when the roadbed slope moves to the system equilibrium state under the action of gravity load, the displaced soil particles are mainly concentrated at the bottom of the roadbed slope, and the maximum displacement of the soil particles is 0.354m.

[0019] Furthermore, in step (5), the contact force between the soil particles in contact with each other is related to the displacement, and the contact force F is described according to the position vector of the contact point between the soil particles. i Decomposed into normal contact forces along the normal line and the tangential contact force along the tangent line

[0020]

[0021] Furthermore, in step (5), the normal contact force and the tangential contact force are calculated respectively; the calculation formula of the normal contact force is:

[0022]

[0023] in, is the normal stiffness of the contact point; is the overlap amount; n i is the line vector between adjacent soil particles;

[0024] The tangential contact force is calculated in the form of increments. When the contact between soil particles is formed, the total tangential contact force is initialized to zero. The elastic tangential contact force caused by the subsequent relative tangential vector displacement is added to the present value. The calculation formula for this process is:

[0025]

[0026] in, Represents the tangential force component The new component obtained after a specific rotation operation (here marked as "rot1"), is the tangential force component in the j direction, δ ij is the second-order counting symbol, e ijk is the third-order counter, e kmn is the third-order counting symbol, is the unit normal vector of the contact surface at the last iteration, n n are the components of the unit vector, Represents the tangential force component The new component obtained after a specific rotation operation (labeled as "rot2"), ω R is the average rotational angular velocity of soil particles about the normal direction, <ω R > is the average rotation angular velocity vector, Δt is the calculation time interval, and are the angular velocity components of two different soil particles in the reference frame, n j is the direction cosine component, is the relative tangential velocity component, V I is the absolute velocity of soil particles, V j is the component of the velocity vector of the soil particles, To express the increment of the tangential force component in the i direction, K s is the tangential stiffness, is the increment of the tangent vector displacement in the i direction, is the bending moment at the contact point, M is the initial angular momentum, is the contact point position, is the translational displacement of the contact point, F k is the kth component of the force vector.

[0027] Furthermore, in step (6), the motion state of the soil particles is determined by the force, bending moment and torque acting on the soil particles; the calculation formulas for the various parameters involved in the motion process of the soil particles are:

[0028] F i =m(x i -g i )(10)

[0029] M1=I1ω°1+(I3-I2)ω3ω2(11)

[0030] M2=I2ω°2+(I1-I3)ω1ω3(12)

[0031] M3=I33+(I2-I1)ω2ω1(13)

[0032] Among them, F i is the contact force, m is the mass of the soil particles, x i is the acceleration of soil particles, g i is the acceleration of gravity, M1, M2, and M3 are the bending moment components of soil particles; I1, I2, and I3 are the rotational inertia of soil particles; ω1, ω2, and ω3 are the angular displacements of soil particles, and ω°1, ω°2, and ω°3 are the angular velocities of soil particles.

[0033] Furthermore, in step (6), the motion state of soil particles is an iterative process; through cyclic calculation until the equilibrium state, the cyclic calculation formula is:

[0034]

[0035]

[0036] in, is the approximate value of the second-order derivative of the soil particle position X with respect to time t at time t, that is, the acceleration; and are the approximate values ​​of the soil particle position X at time (t+Δt / 2) and (t1-Δt / 2), respectively, Δt is the calculation time step, is the angular acceleration at time t, and are the angular velocity components at time (t+Δt / 2) and (t-Δt / 2), is the contact force F i The component at time t, is the angular velocity component at time (t-Δt / 2), is the component of the moment M at time t, is the component of the soil particle position X at time t+Δt, is the component of the soil particle position X at time t, is the component of the soil particle position X at time (t+Δt / 2), and I is the moment of inertia.

[0037] The beneficial effects of the present invention are:

[0038] The present invention provides a numerical simulation method for simulating the influence of root-soil complex on slope stability using discrete element software, which is mainly used to realize the simulation of plain soil slope and the whole process of root-soil complex slope. The present invention can effectively restore the actual condition of the slope by simulating and simplifying the plant root system. Based on the plain soil slope, the displacement of soil particles inside the slope after the plain soil slope is naturally excavated can be accurately calculated. Based on the root-soil complex slope, the displacement of soil particles inside the slope after the slope is excavated under the action of the plant root system can be calculated. The slope-stabilizing effect of the root-soil complex is analyzed by analyzing the displacement of soil particles under the two states. The present invention restores the actual rooted slope condition, is easy to operate, has high calculation efficiency, and has practical significance and application value in the field of geotechnical engineering. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 It is a plain soil model.

[0040] Figure 2 Simplified model of a plant root system.

[0041] Figure 3 It is a root-soil complex slope model.

[0042] Figure 4 This is the displacement cloud diagram of slope soil particles in the bare soil state.

[0043] Figure 5 This is the displacement cloud diagram of slope soil particles in the root-soil complex state.

[0044] Figure 6 It is the maximum displacement of soil particles in the state of bare soil and root-soil complex. DETAILED DESCRIPTION

[0045] The present invention is further described in detail below with reference to the accompanying drawings.

[0046] The present invention provides a numerical simulation method for simulating the influence of root-soil complex on slope stability using discrete element software, which mainly includes the simulation process of plant roots and the analysis of root-soil complex slope stability. The specific implementation process is as follows:

[0047] Step 1: Simulation process of plant root system;

[0048] S1.1: Slope soil parameters are obtained through field construction data and indoor tests, including various soil performance parameters and rooted soil related parameters, see Table 1.

[0049] Table 1

[0050]

[0051] S1.2: Based on the on-site construction data and the high fill roadbed slope model (such as the existing roadbed slope design drawings), use discrete element software (PFC) to establish a soil model, and set the slope ratio of the soil model to 1:1.5. Figure 1 As shown in the figure, the soil model adopts fixed boundary conditions, that is, no constraints are imposed on the upper boundary, horizontal constraints are imposed on the lower boundary while vertical constraints are imposed, and only horizontal constraints are imposed on the left, right, front and back boundaries.

[0052] S1.3: Simplify plant root systems;

[0053] First, draw the required shape in CAD software, then import it into the particle flow analysis program PFC2D to form a strip unit (geometry, a unit name in PFC2D), and then convert the geometry into a super particle cluster; finally, import the generated super particle cluster into the plain soil model to generate a root-soil complex slope model. The simplified model of the plant root system is as follows: Figure 2 The generated root-soil composite slope model is shown in Figure 3 shown.

[0054] Step 2: Root-soil complex slope stability analysis;

[0055] S2.1: For the bare soil model and the root-soil composite slope model, gravity loads were applied under different conditions (different conditions refer to different slope states, one is bare soil without roots, the other is root-soil composite slope, and the gravity load range is 15-20 kPa);

[0056] S2.2: Analysis of soil particle displacement in bare soil state;

[0057] like Figure 4 As shown in Figure 2, when the roadbed slope moves to the system equilibrium state under the action of gravity load, some soil particles are displaced, and the displacement of soil particles is mainly concentrated on the roadbed slope.

[0058] S2.3: Analysis of soil particle displacement in the root-soil complex state;

[0059] like Figure 5 As shown in the figure, in the root-soil composite slope model, when the same gravity load is applied and the system is operated to a balanced state, the displacement of most soil particles is effectively alleviated compared with the plain soil state, and the displaced soil particles are mainly concentrated at the bottom of the roadbed slope.

[0060] S2.4: Analysis of maximum displacement of soil particles in bare soil and root-soil composite states;

[0061] The bare soil model and the root-soil composite slope model are run to the system equilibrium state, and the displacement of the roadbed slope soil particles of the two models are recorded at the same interval. Figure 6As can be seen, when the calculation time step reaches 100,000, the maximum displacement of the roadbed slope soil particles in the bare soil state and the root-soil complex state approach each other. When the system reaches equilibrium, the maximum displacement of the roadbed slope soil particles in the bare soil state reaches 0.53m, and the maximum displacement of the roadbed slope soil particles in the root-soil complex state reaches 0.354m. The displacement suppression effect of soil particles reaches 33.2%, indicating that the plant roots effectively strengthen the slope soil.

[0062] S2.5: Force-displacement law;

[0063] The contact force between soil particles in contact with each other is related to the displacement and described according to the position vector of the contact point between the soil particles. The contact force F i Decomposed into normal contact forces along the normal line (normal force component) and the tangential contact force along the tangent line (tangential force component), that is:

[0064]

[0065] The specific calculation formula of normal contact force is as follows:

[0066]

[0067] in, is the normal stiffness of the contact point, (E is the elastic modulus, r is the radius of the spherical particle), the unit is N / m 3 ; is the overlap amount; n i is the line vector between adjacent soil particles.

[0068] The tangential contact force is calculated in increments. When contact is established between soil particles, the total tangential contact force is initialized to zero. The elastic tangential contact force caused by the subsequent relative tangential vector displacement is added to the current value. The specific calculation formulas involved in this process are as follows: Formula (3) to Formula (9):

[0069]

[0070] in, Represents the tangential force component The new component obtained after a specific rotation operation (here marked as "rot1"), is the tangential force component in the j direction, δ ij is the second-order counting symbol, e ijk is the third-order counter, e kmn is the third-order counting symbol, is the unit normal vector of the contact surface at the last iteration, n nare the components of the unit vector, Represents the tangential force component The new component obtained after a specific rotation operation (labeled as "rot2"), ω R is the average rotational angular velocity of soil particles about the normal direction, <ω R > is the average rotation angular velocity vector, Δt is the calculation time interval, and are the angular velocity components of two different soil particles in the reference frame, n j is the direction cosine component, is the relative tangential velocity component, V I is the absolute velocity of soil particles, V j is the component of the velocity vector of the soil particles, To express the increment of the tangential force component in the i direction, K s is the tangential stiffness, is the increment of the tangent vector displacement in the i direction, is the bending moment at the contact point, M is the initial angular momentum, is the contact point position, is the translational displacement of the contact point, F k is the kth component of the force vector.

[0071] S2.6: Laws of motion;

[0072] The motion state of soil particles is determined by the forces, bending moments, and torques acting on them. This is an iterative process. The specific calculation formulas for the various parameters involved in this process are as follows:

[0073] F i =m(x i -g i )(10)

[0074] M1=I1ω°1+(I3-I2)ω3ω2(11)

[0075] M2=I2ω°2+(I1-I3)ω1ω3(12)

[0076] M3=I33+(I2-I1)ω2ω1(13)

[0077] Among them, F i is the contact force, m is the mass of the soil particles, x i is the acceleration of soil particles, g i is the acceleration of gravity, M1, M2, and M3 are the bending moment components of soil particles; I1, I2, and I3 are the rotational inertia of soil particles; ω1, ω2, and ω3 are the angular displacements of soil particles, and ω°1, ω°2, and ω°3 are the angular velocities of soil particles.

[0078] The cycle calculation can be performed by iterating according to formula (14)-formula (18) until the equilibrium state is reached:

[0079]

[0080] in, is the approximate value of the second-order derivative of the soil particle position X with respect to time t at time t, that is, the acceleration; and are the approximate values ​​of the soil particle position X at time (t+Δt / 2) and (t1-Δt / 2), respectively, Δt is the calculation time step, is the angular acceleration at time t, and are the angular velocity components at time (t+Δt / 2) and (t-Δt / 2), is the contact force F i The component at time t, is the angular velocity component at time (t-Δt / 2), is the component of the moment M at time t, is the component of the soil particle position X at time t+Δt, is the component of the soil particle position X at time t, is the component of the soil particle position X at time (t+Δt / 2), and I is the moment of inertia.

[0081] Therefore, by According to formula (16) and formula (17), it can be calculated that Then we can calculate it by formula (18): Then the next iteration is performed until the equilibrium state is reached.

[0082] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications should also be regarded as within the scope of protection of the present invention.

Claims

1. A numerical simulation method for simulating the influence of root-soil complex on slope stability using discrete element software, characterized in that: The following steps are involved: (1) Obtain slope soil parameter data; (2) Based on the on-site construction data and the high fill roadbed slope model, a soil model was established using discrete element software; (3) Use CAD and PFC2D to simplify the plant root system and obtain a simplified plant root system model and a root-soil complex slope model; (4) By applying the same gravity load and running the system to the equilibrium state, the displacement of soil particles in the state of bare soil and root-soil complex is analyzed; (5) Force-displacement law analysis of soil particles; (6) Analysis of the motion laws of soil particles.

2. The numerical simulation method of the influence of root-soil complex on slope stability using discrete element software according to claim 1 is characterized in that: In step (2), the slope ratio of the bare soil model is 1:1.

5.

3. The numerical simulation method of the influence of root-soil complex on slope stability using discrete element software according to claim 1 is characterized in that: In step (2), the bare soil model adopts fixed boundary conditions, that is, no constraints are imposed on the upper boundary, horizontal constraints are imposed on the lower boundary while vertical constraints are imposed, and only horizontal constraints are imposed on the left, right, front and rear boundaries.

4. The numerical simulation method of the influence of root-soil complex on slope stability using discrete element software according to claim 1 is characterized in that: In step (3), the required shape diagram is first drawn in the CAD software, and then imported into the particle flow analysis program PFC2D to form a strip unit. The strip unit is then converted into a super particle cluster, and the super particle cluster is imported into the plain soil model to generate a root-soil complex slope model.

5. The numerical simulation method of the influence of root-soil complex on slope stability using discrete element software according to claim 1 is characterized in that: In step (4), in the bare soil state, when the roadbed slope moves to the system equilibrium state under the action of gravity load, some soil particles are displaced. The displacement of soil particles is mainly concentrated at the roadbed slope, and the maximum displacement of soil particles is 0.53m.

6. The numerical simulation method of the influence of root-soil complex on slope stability using discrete element software according to claim 1 is characterized in that: In step (4), in the root-soil complex state, when the roadbed slope moves to the system equilibrium state under the action of gravity load, the displaced soil particles are mainly concentrated at the bottom of the roadbed slope, and the maximum displacement of the soil particles is 0.354m.

7. The numerical simulation method of the influence of root-soil complex on slope stability using discrete element software according to claim 1 is characterized in that: In step (5), the contact force between the soil particles in contact with each other is related to the displacement, and the contact force F is described according to the position vector of the contact point between the soil particles. i Decomposed into normal contact forces along the normal line and the tangential contact force along the tangent line 8. The numerical simulation method of the influence of root-soil complex on slope stability using discrete element software according to claim 7 is characterized in that: In step (5), the normal contact force and the tangential contact force are calculated respectively; the calculation formula of the normal contact force is: in, is the normal stiffness of the contact point; is the overlap amount; n i is the line vector between adjacent soil particles; The tangential contact force is calculated in the form of increments. When the contact between soil particles is formed, the total tangential contact force is initialized to zero. The elastic tangential contact force caused by the subsequent relative tangential vector displacement is added to the present value. The calculation formula for this process is: in, Represents the tangential force component The new component obtained after a specific rotation operation, is the tangential force component in the j direction, δ ij is the second-order counting symbol, e ijk is the third-order counter, e kmn is the third-order counting symbol, is the unit normal vector of the contact surface at the last iteration, n n are the components of the unit vector, Represents the tangential force component The new component obtained after a specific rotation operation, ω R is the average rotational angular velocity of soil particles about the normal direction, <ω R > is the average rotation angular velocity vector, Δt is the calculation time interval, and are the angular velocity components of two different soil particles in the reference frame, n j is the direction cosine component, is the relative tangential velocity component, V I is the absolute velocity of soil particles, V j is the component of the velocity vector of the soil particles, To express the increment of the tangential force component in the i direction, K s is the tangential stiffness, is the increment of the tangent vector displacement in the i direction, is the bending moment at the contact point, M is the initial angular momentum, is the contact point position, is the translational displacement of the contact point, F k is the kth component of the force vector.

9. The numerical simulation method of the influence of root-soil complex on slope stability using discrete element software according to claim 1 is characterized in that: In step (6), the motion state of soil particles is determined by the force, bending moment and torque acting on the soil particles; the calculation formulas for the various parameters involved in the motion process of soil particles are: F i =m(x i -g i ) (10) M1=I1ω o 1+(I3-I2)ω3ω2 (11) M2=I2ω o 2+(I1-I3)ω1ω3 (12) M3=I33+(I2-I1)ω2ω1 (13) Among them, F i is the contact force, m is the mass of the soil particles, x i is the acceleration of soil particles, g i is the acceleration of gravity, M1, M2, and m3 are all the bending moment components of soil particles; I1, I2, and I3 are all the rotational inertia of soil particles; ω1, ω2, and ω3 are all the angular displacements of soil particles, ω ° 1.ω ° 2.ω ° 3 are the angular velocities of soil particles.

10. The numerical simulation method of the influence of root-soil complex on slope stability using discrete element software according to claim 9 is characterized in that: In step (6), the motion state of soil particles is an iterative process; through cyclic calculation until the equilibrium state, the cyclic calculation formula is: in, is the approximate value of the second-order derivative of the soil particle position X with respect to time t at time t, that is, the acceleration; and are the approximate values ​​of the soil particle position X at time (t+Δt / 2) and (t1-Δt / 2), respectively, Δt is the calculation time step, is the angular acceleration at time t, and are the angular velocity components at time (t+Δt / 2) and (t-Δt / 2), is the contact force F i The component at time t, is the angular velocity component at time (t-Δt / 2), is the component of the moment M at time t, is the component of the soil particle position X at time t+Δt, is the component of the soil particle position X at time t, is the component of the soil particle position X at time (t+Δt / 2), and I is the moment of inertia.

Images