A soil cracking simulation method based on two-dimensional particle discrete element method

Abstract

This invention discloses a soil cracking simulation method based on two-dimensional particle discrete elements, relating to the field of soil physics. The method includes the following steps: discretizing the region to be solved using particles and assigning control parameters; performing contact searches between particles and between particles and their boundaries to determine the domains of mechanical and water diffusion, establishing force chains and virtual water diffusion channels between them; calculating water diffusion between particles and between particle boundaries and updating the water content of the particles; and performing mechanical calculations between particles and between particle boundaries and updating crack information. This method can simulate soil cracking due to water loss in two dimensions.

Description

Technical Field

[0001] The present invention relates to the field of soil physics, and more particularly to a soil desiccation simulation method based on two-dimensional particle discrete elements. Background Art

[0002] Currently, laboratory experiments are often used to study soil cracking. However, these studies focus on phenomenological observations and lack insight into underlying mechanisms. Furthermore, sample preparation and laboratory or field testing are often expensive and time-consuming. To address these issues, several numerical methods have been developed to simulate desiccation-induced soil cracking. Continuum-based methods, such as the finite element method (FEM), meshless methods, finite volume methods, and numerical manifold methods, are commonly used to study water diffusion in media. However, when considering the coupling of water diffusion and mechanical properties, involving complex contact or cracking problems, continuum-based methods have limited success due to their difficulty in simulating crack initiation and propagation. In contrast, discontinuous methods, such as discrete element methods (DEM), discontinuous deformation analysis, and particle dynamics, can better describe material cracking behavior. However, few researchers have yet to couple water diffusion models with discrete element models for soil cracking, and this work remains flawed. Summary of the Invention

[0003] In view of this, the present invention provides a soil cracking simulation method based on two-dimensional particle discrete elements, which can simulate the soil dehydration and cracking problem in two-dimensional conditions.

[0004] In order to achieve the above object, the present invention adopts the following technical solutions:

[0005] A soil cracking simulation method based on two-dimensional particle discrete element method includes the following steps:

[0006] Use particles to discretize the area to be solved and assign control parameters;

[0007] Conduct contact searches between particles and between particles and boundaries, determine the domains of mechanics and water diffusion, and establish force chains and virtual water diffusion channels between them;

[0008] Calculate the water diffusion between particles and at the boundaries of particles and update the water content of particles;

[0009] Based on the updated particle moisture content, mechanical calculations are performed between particles and between particle boundaries and crack information is updated.

[0010] Optionally, the control parameters include geometric parameters, mechanical parameters and water diffusion parameters; the geometric parameters include particle position coordinate information and particle size; the mechanical parameters include particle density ρ, normal contact stiffness K n, tangential contact stiffness K s , tensile strength f t , bond strength C, internal friction angle and the shrinkage coefficient α s ; The conduction parameters include the water diffusion coefficient k and the boundary flow mass density q w .

[0011] Optionally, the two-dimensional particle model uses a finite number of random or closely arranged particles to discretize the area to be solved, and uses walls to impose constraints on the boundaries of the solution area.

[0012] Optionally, the particles in the two-dimensional particle model are subjected to a normal force F n and tangential force F s , acting on their respective centers, the normal force is along the line between the centers, and the tangential force F n Then perpendicular to the normal force F s , normal force F n With tangential force F s All obey Hooke's law as follows:

[0013] ΔF n =-K n Δd n ;

[0014] ΔF s =-K s Δd s ;

[0015] Where K n , K s are the normal contact stiffness and the tangential contact stiffness respectively; Δd n , Δd s are the changes in the normal and tangential relative displacements of the contact pair within each calculation time step, respectively; the negative sign indicates that the spring force hinders the increase in relative displacement.

[0016] Optionally, in the two-dimensional particle model, in a continuous medium, water diffusion follows the exchange rate per unit time being proportional to the water content gradient, as shown below:

[0017]

[0018] Where q is the exchange rate per unit time, k is the water diffusion coefficient, is the water content gradient; the negative sign indicates that water is transferred from the higher water content end to the lower water content end;

[0019] The total flow Q flows to a point with a given dry mass M, and the change in water content at that point is defined by the following formula:

[0020]

[0021] Optionally, the two-dimensional particle model compensates and adjusts the particle density, and the compensated particle density ρ p As shown below;

[0022]

[0023] Where ρ is the density of the continuous medium, f v is the volume fraction of the particle system.

[0024] Optional water diffusion methods between particles and at particle boundaries are as follows:

[0025] The radius is R j With R i The transmission surface of the virtual channel for water diffusion between particles j and i is a rectangle, and the width of the rectangle is defined as the arithmetic mean of the diameters of the two particles, that is, the sum of the radii of the two particles; the cross-sectional width S of the uncorrected virtual water diffusion channel is ij As shown below:

[0026] S ij =R i +R j ;

[0027] The distance L between the centers of the two particles ij The calculation is as follows

[0028]

[0029] The center coordinates of particle i are (x i ,y i ), the center coordinate of particle j is (x j ,y j );

[0030] Moisture diffusion coefficient k of the virtual moisture diffusion channel ij The harmonic mean of the water diffusion coefficients of the two particles is shown as follows:

[0031]

[0032] Where k i and k j are the water diffusion coefficients of particles i and j, respectively;

[0033] according to The flow rate Q from particle j to particle i per unit time j→i It is expressed by the following formula:

[0034]

[0035] Where Wi and W j are the water contents of particles i and j, respectively; α is a correction coefficient introduced to correct the width of the virtual water diffusion channel, which is calculated by the following formula:

[0036]

[0037] Where the coordination number CN is a measure of the average number of adjacent particles for each particle, f v represents the volume fraction of particle aggregates;

[0038] The flow rate Q from the boundary w to the particle i per unit time w→i as follows:

[0039] Q w→i =βq w S iw ;

[0040] Where, the width of the virtual moisture diffusion channel S is iw And the correction coefficient β is as follows:

[0041] S iw =2R i ;

[0042]

[0043] Where n is the number of particles in contact with the boundary w.

[0044] Optionally, update the moisture content of the pellets as follows:

[0045] The flow rate flowing into particle i is the sum of the flow rates flowing into particle i from all adjacent particles and adjacent boundaries. Particle i has n adjacent particles and m adjacent boundaries. The total flow rate flowing into particle i per unit time is expressed by the following formula:

[0046]

[0047] The total flow rate Q obtained by a particle i per unit time i , combined with Calculate the water content change ΔW within a time step Δt i , and then complete the update of the moisture content of particle i, and calculate:

[0048] M i =ρ p πR i 2 ;

[0049]

[0050] Where M iis the dry mass of particle i, and Δt is the time step.

[0051] Optionally, a linear strain method is applied to reflect the effect of water content changes on mechanics, as follows:

[0052] The initial radius of particle i is R i0 , the water content changes by ΔW i , then its new radius R after water loss and shrinkage i Expressed as:

[0053] R i =R i0 (1+α s ΔW i );

[0054] Where α s is the linear shrinkage coefficient;

[0055] The distance between particles i and j in the normal direction is:

[0056]

[0057] The coordinates of the centers of the two particles remain unchanged. Only the effect of the change in water content is considered. The water content of particle i changes by ΔW within the calculation step. i , the volume of particle i shrinks, causing the distance between the two particles in the normal direction to change Δd n ,have to:

[0058] Δd n =-R i0 α s ΔW i ;

[0059] The increment of the normal force between the two particles ΔF n for:

[0060] ΔF n =K n R i0 α s ΔW i ;

[0061] F t is the maximum allowable contact tension, F c is the bonding contact force, is the internal friction angle, where:

[0062] F t =2R·f t ;

[0063] F c =2R·C.

[0064] Optionally, the flow reduction coefficient method is used to reflect the influence of mechanics on water diffusion, as follows:

[0065] When particles i and j are in a bonded state, the spring breaks, the bond fails, and microcracks are generated. The water diffusion coefficient will decrease, and the water diffusion coefficient k of the virtual channel will be ij Calculated according to the following formula:

[0066]

[0067] Where C r is the flow reduction coefficient, 0≤C r ≤1, where k i and k j are the water diffusion coefficients of particles i and j, respectively.

[0068] It can be seen from the above technical solution that, compared with the prior art, the present invention provides a soil cracking simulation method based on two-dimensional particle discrete elements, which has the following beneficial effects:

[0069] 1. Realize the coupling of water diffusion model and discrete element mechanics model.

[0070] 2. Use the linear strain method to reflect the effect of water content changes on mechanical action, and then realize water loss shrinkage

[0071] 3. Flow reduction coefficient C r , considering the hindering effect of cracks on water diffusion, it can be determined based on experiments or engineering according to different material working conditions, so as to ensure the accuracy of calculation simulation. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are merely embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying any creative work.

[0073] Figure 1a Schematic diagram of the bonding particle model of the present invention;

[0074] Figure 1b Schematic diagram of the bonding particle model of the present invention;

[0075] Figure 2a Schematic diagram of the first diffusion of water between particles and between particle boundaries through virtual channels of the present invention;

[0076] Figure 2bSchematic diagram of the second diffusion of water between particles and between particle boundaries through virtual channels of the present invention;

[0077] Figure 3 This is a schematic diagram of the shrinkage of particles according to the present invention;

[0078] Figure 4a It is a linear schematic diagram of the normal force of the present invention;

[0079] Figure 4b It is a linear schematic diagram of the tangential force of the present invention;

[0080] Figure 5 Schematic diagram of the barrier effect of cracks on water diffusion in the present invention;

[0081] Figure 6 This is a diagram of a soil restrained drying test specimen of the present invention;

[0082] Figure 7 The diagram shows the distribution of water content and crack expansion of the specimen at different times;

[0083] Figure 8 Flowchart of the method of the present invention. DETAILED DESCRIPTION

[0084] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0085] The embodiment of the present invention discloses a soil cracking simulation method based on two-dimensional particle discrete elements, which can simulate the soil dehydration cracking problem in two-dimensional conditions. Figure 8 As shown, the following steps are included:

[0086] Step 1: Use particles to discretize the area to be solved and assign control parameters;

[0087] Step 2: Conduct contact searches between particles and between particles and boundaries to determine the domains of mechanics and water diffusion, and establish force chains and virtual water diffusion channels between them.

[0088] Step 3: Calculate the water diffusion between particles and between particle boundaries and update the water content of the particles;

[0089] Step 4: Perform mechanical calculations between particles and between particle boundaries and update crack information.

[0090] Repeat steps 1 to 4 until the calculation ends.

[0091] Furthermore, in step 1, the control parameters include geometric parameters, mechanical parameters and water diffusion parameters. Geometric parameters include particle position coordinate information and particle size; mechanical parameters include particle density ρ, normal contact stiffness K n , tangential contact stiffness K s , tensile strength f t , bond strength C, internal friction angle and the shrinkage coefficient α s ; The conduction parameters include the water diffusion coefficient k and the boundary flow mass density q w wait.

[0092] Furthermore, the 2D particle model (2DPM) uses a finite number of random or tightly packed particles to discretize the region to be solved and uses walls to impose constraints on the boundaries of the solution region. This model uses the following assumptions:

[0093] 1) In discrete element modeling of continuum media, the domain to be solved is discretized into densely packed particles. Contacting particles are connected using a bonded particle model. In 2DPM, contact forces and water diffusion only occur between adjacent particles and between adjacent particles and their boundaries. Therefore, the particle contact mesh can be considered the domain for mechanical and water diffusion effects.

[0094] 2) The particles are connected by a set of normal springs and tangential springs. The forces act on the centers of the circles respectively. The normal force is along the line between the centers of the circles, and the tangential force is perpendicular to the normal force (e.g. Figure 1a 、 Figure 1b shown).

[0095] The diffusion of water between particles and at particle boundaries occurs through Figure 2a 、 Figure 2b The virtual water diffusion channels between particles and between particle boundaries are rectangular.

[0096] 3) Normal force F n With tangential force F s They all follow Hooke's law, as shown in Formula 1-2.

[0097] ΔF n =-K n Δd n (1

[0098] ΔF s =-K s Δd s (2

[0099] Where K n , K sare the normal contact stiffness and the tangential contact stiffness respectively; Δd n , Δd s are the changes in the normal and tangential relative displacements of the contact pair within each calculation time step, respectively. The negative sign indicates that the spring force hinders the increase in relative displacement.

[0100] 4) In a continuous medium, water diffusion should follow the basic assumption that the exchange rate per unit time is proportional to the water content gradient, as shown in Formula 3.

[0101]

[0102] Where q is the exchange rate per unit time, k is the water diffusion coefficient, is the moisture content gradient. The negative sign indicates that moisture is transferred from the higher moisture content end to the lower moisture content end.

[0103] Assuming that the total flow Q flows to a point with a given dry mass of M, the change in water content of the point can be defined by formula (4).

[0104]

[0105] 5) Since the area to be solved is discretized by particles, and there are inevitably gaps between the particles, which may cause the total mass of the area to change, it is necessary to compensate and adjust the particle density to ensure that the total mass of the area to be solved remains unchanged. The compensated particle density ρ p It can be shown by formula 5.

[0106]

[0107] Where ρ is the density of the continuous medium, f v is the volume fraction of the particle system.

[0108] Furthermore, it is characterized by the water diffusion mode between particles and at the particle boundaries, as follows:

[0109] The radius is R j With R i Taking the diffusion of water between particles j and i as an example, the transmission surface of the virtual channel for water diffusion is a rectangle, and the width of the rectangle is defined as the arithmetic mean of the diameters of the two particles, that is, the sum of the radii of the two particles. Based on the above assumptions, the cross-sectional width S of the uncorrected virtual water diffusion channel is ij As shown in formula (6).

[0110] S ij =R i +R j (6

[0111] The distance L between the centers of the two particlesij The calculation is as follows

[0112]

[0113] The center coordinates of particle i are (x i ,y i ), the center coordinate of particle j is (x j ,y j ).

[0114] Define the moisture diffusion coefficient k of the virtual moisture diffusion channel ij is the harmonic mean of the water diffusion coefficients of the two particles

[0115]

[0116] Where k i and k j are the water diffusion coefficients of particles i and j, respectively.

[0117] According to formula 3, the flow rate Q from particle j to particle i per unit time is j→i It can be expressed by the following formula:

[0118]

[0119] Where W i and W j are the water contents of particles i and j, respectively; α is a correction coefficient introduced to correct the width of the virtual water diffusion channel. An appropriate α can correctly reproduce the macroscopic water diffusion coefficient, and its value can be calculated by the following formula:

[0120]

[0121] Where the coordination number CN is a measure of the average number of adjacent particles for each particle, f v represents the volume fraction of particle aggregates.

[0122] The diffusion of water between particle boundaries is very similar to the diffusion of water between particles. Figure 1a As shown, the boundary w is simplified to a line segment in the two-dimensional case, assuming that the particle i is in contact with the boundary w. In order to simulate the evaporation process of water, a flow mass density q is imposed on the boundary w. w , its value is usually negative, indicating the loss of water per unit time. Then the flow rate Q from the boundary w to the particle i per unit time is w→i You can write:

[0123] Q w→i =βq w S iw (11

[0124] Where, the width of the virtual moisture diffusion channel S is iw And the correction coefficient β is as follows:

[0125] S iw =2R i (12

[0126]

[0127] Where n is the number of particles in contact with the boundary w.

[0128] Furthermore, the iterative method for updating the water content during water diffusion is as follows:

[0129] The flow rate flowing into particle i is the sum of the flow rates flowing into particle i from all adjacent particles and adjacent boundaries. Assuming that particle i has n adjacent particles and m adjacent boundaries, the total flow rate flowing into particle i per unit time can be expressed as follows:

[0130]

[0131] Once the total flow rate Q obtained by a particle i per unit time is obtained i , we can combine formula (4) to calculate the water content change ΔW within a time step Δt i , and then complete the update of the moisture content of particle i, and calculate:

[0132] M i =ρ p πR i 2 (15

[0133]

[0134] Where M i is the dry mass of particle i, and Δt is the time step.

[0135] Furthermore, the linear strain method is applied to reflect the influence of water content change on mechanics, as follows:

[0136] Assume that the initial radius of particle i is R i0 , within one calculation step, the water content changes by ΔW i , then its new radius R after water loss and shrinkage i It can be expressed as:

[0137] R i =R i0 (1+α s ΔW i ) (18

[0138] Where α s is the linear shrinkage coefficient.

[0139] Considering that the distance between particles i and j in the normal direction is:

[0140]

[0141] like Figure 3 As shown, it is assumed that the coordinates of the centers of the two particles remain unchanged and only the effect of the change in water content is considered. The water content of particle i changes by ΔW in this calculation step. i , the volume of particle i shrinks, causing the distance between the two particles in the normal direction to change Δd n , can be combined with formula (18) (19) to get

[0142] Δd n =-R i0 α s ΔW i (20

[0143] In 2DPM, according to formula (1), it can be known that the normal force increment ΔF between the two particles within this calculation step is n for

[0144] ΔF n =K n R i0 α s ΔW i (twenty one

[0145] Currently, only the mechanical interaction between one set of contact pairs is analyzed, focusing on the normal force F n If we generalize this to a wider range of situations, since the water content of each particle may be different, a single particle may be subjected to uneven shrinkage stress from the surrounding adjacent particles, which in turn causes particle movement. This movement will cause the coordinates of the particle center to change and cause a tangential force F s dynamic changes.

[0146] After obtaining the stress conditions of each contact pair, the crack information needs to be updated. 2DPM uses the Moore-Coulomb criterion as the constitutive model and uses the spring fracture to simulate the initiation and propagation of cracks. Figure 4a and Figure 4b As shown in the figure, F t is the maximum allowable contact tension, F c is the bonding contact force, is the internal friction angle, where:

[0147] F t =2R·f t (twenty two

[0148] F c =2R·C (23

[0149] Once the normal force F n F exceeding the normal spring t , the tensile strength of the spring decreases to 0 and a tensile crack is recorded; once the tangential force exceeds the allowable value, the adhesive contact force is lost and a shear crack is recorded.

[0150] Furthermore, the flow reduction coefficient method is used to reflect the influence of mechanics on water diffusion, as follows:

[0151] like Figure 5 As shown in Figure 1, after the continuum is divided by cracks, the cracks will hinder the water diffusion process. In order to characterize this effect, the 2DPM model introduces a flow reduction coefficient. Assuming that a contact pair consists of particles i and j, when the two are in a bonded state and the connecting spring is not broken, the water diffusion coefficient of the virtual channel between the particles can be calculated by formula (8). Once the spring breaks, the bond fails and microcracks are generated. At this time, the water diffusion coefficient will decrease, and the water diffusion coefficient of the virtual channel will decrease. It can be calculated according to the following formula:

[0152]

[0153] Where C r is the flow reduction coefficient, 0≤C r ≤1, formula 24 shows that the coefficient C r The smaller the crack is, the more significant the barrier effect of the crack on water diffusion is. r =1 means that water can be conducted freely without any crack obstruction, C r = 0 means the crack is completely watertight and prevents water from spreading. In practical applications, the correction factor C r Usually determined based on engineering experience or experimental data.

[0154] The following is a soil restrained drying test to verify the effectiveness of the method. Figure 6 As shown in Figure 1, the test specimen is 294 mm long and 16 mm high. The upper portion of the specimen is soil, and the lower portion is a metal base with a 2 mm groove. The grooves in the metal base constrain the soil deformation in the x-direction. A unit area mass flow boundary is applied to the upper boundary of the soil sample, with a flow density of -5.6 × 10 -6 kg / (m 2 ·s) to represent the evaporation of water in the soil. The other boundaries are impermeable and the initial moisture content of the soil is 49%. For simplicity, the shrinkage coefficient α of the model is s The moisture diffusion coefficient k remains fixed during the simulation process. The macroscopic and microscopic parameters used in 2DPM are shown in Table 1.

[0155] Table 1 Macroscopic and microscopic parameters of 2DPM

[0156]

[0157] The distribution of water content and crack expansion of specimens at different times, such as Figure 7 As shown in the figure, cracks initiate from the soil surface and propagate inward. Furthermore, during the simulation, the distribution of soil moisture content on both sides of the crack is discontinuous, demonstrating that this technical solution fully considers the crack's barrier to moisture diffusion. The simulation results are highly consistent with the experimental results, demonstrating the effectiveness of the soil desiccation and cracking simulation method based on 2D granular discrete element methods.

[0158] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. Reference can be made to the common and similar parts between the various embodiments. For the devices disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple, and the relevant parts can be referred to the method description.

[0159] The above description of the disclosed embodiments is intended to enable one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be readily apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention is not limited to the embodiments shown herein but is intended to conform to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A soil cracking simulation method based on two-dimensional particle discrete element method, characterized in that: The following steps are involved: Use particles to discretize the area to be solved and assign control parameters; Conduct contact searches between particles and between particles and boundaries, determine the domains of mechanics and water diffusion, and establish force chains and virtual water diffusion channels between them; Calculate the water diffusion between particles and at the boundaries of particles and update the water content of particles; Perform mechanical calculations between particles and between particle boundaries based on the updated particle moisture content and update crack information; The water diffusion between particles and at particle boundaries is as follows: The radius is R j With R i The transmission surface of the virtual channel for water diffusion between particles j and i is a rectangle, and the width of the rectangle is defined as the arithmetic mean of the diameters of the two particles, that is, the sum of the radii of the two particles; the cross-sectional width S of the uncorrected virtual water diffusion channel is ij As shown below: S ij =R i +R j ; The distance L between the centers of the two particles ij The calculation is as follows The center coordinates of particle i are (x i ,y i ), the center coordinate of particle j is (x j ,y j ); Moisture diffusion coefficient k of the virtual moisture diffusion channel ij The harmonic mean of the water diffusion coefficients of the two particles is shown as follows: Where k i and k j are the water diffusion coefficients of particles i and j, respectively; According to q = -k▽W, the flow rate Q from particle j to particle i per unit time j→i It is expressed by the following formula: Where W i and W j are the water contents of particles i and j, respectively; α is a correction coefficient introduced to correct the width of the virtual water diffusion channel, which is calculated by the following formula: Where the coordination number CN is a measure of the average number of adjacent particles for each particle, f v represents the volume fraction of particle aggregates; The flow rate Q from the boundary w to the particle i per unit time w→i as follows: Q w→i =βq w S iw ; Where, the width of the virtual moisture diffusion channel S iw And the correction coefficient β is as follows: S iw =2R i ; Where n is the number of particles in contact with the boundary w.

2. The soil cracking simulation method based on two-dimensional particle discrete element according to claim 1 is characterized in that: The control parameters include geometric parameters, mechanical parameters and water diffusion parameters; the geometric parameters include particle position coordinate information and particle size; the mechanical parameters include particle density ρ, normal contact stiffness K n , tangential contact stiffness K s , tensile strength f t , bond strength C, internal friction angle and the shrinkage coefficient α s ; The conduction parameters include the water diffusion coefficient k and the boundary flow mass density q w .

3. The soil cracking simulation method based on two-dimensional particle discrete element according to claim 1 is characterized in that: The two-dimensional particle model uses a finite number of random or tightly packed particles to discretize the region to be solved, and uses walls to impose constraints on the boundaries of the solution region.

4. The soil cracking simulation method based on two-dimensional particle discrete element method according to claim 1 is characterized in that: In the two-dimensional particle model, the particles are subject to a normal force F n and tangential force F s , acting on their respective centers, the normal force is along the line between the centers, and the tangential force F n Then perpendicular to the normal force F s , normal force F n With tangential force F s All obey Hooke's law as follows: ΔF n =-K n Δd n ; ΔF s =-K s Δd s ; Where K n , K s are the normal contact stiffness and the tangential contact stiffness respectively; Δd n , Δd s are the changes in the normal and tangential relative displacements of the contact pair within each calculation time step, respectively; the negative sign indicates that the spring force hinders the increase in relative displacement.

5. The soil cracking simulation method based on two-dimensional particle discrete element method according to claim 1 is characterized in that: In a two-dimensional particle model, in a continuous medium, water diffusion follows the exchange rate per unit time that is proportional to the water content gradient, as shown below: Where q is the exchange rate per unit time, k is the water diffusion coefficient, is the water content gradient; the negative sign indicates that water is transferred from the higher water content end to the lower water content end; The total flow Q flows to a given number of points of mass M, and the change in water content at the point is defined by the following formula:

6. The soil cracking simulation method based on two-dimensional particle discrete element method according to claim 1 is characterized in that: The two-dimensional particle model compensates for the particle density, and the compensated particle density ρ p As shown below; Where ρ is the density of the continuous medium, f v is the volume fraction of the particle system.

7. The soil cracking simulation method based on two-dimensional particle discrete element method according to claim 1 is characterized in that: Update the moisture content of the pellets as follows: The flow rate flowing into particle i is the sum of the flow rates flowing into particle i from all adjacent particles and adjacent boundaries. Particle i has n adjacent particles and m adjacent boundaries. The total flow rate flowing into particle i per unit time is expressed by the following formula: The total flow rate Q obtained by a particle i per unit time i , combined with Calculate the water content change ΔW within a time step Δt i , and then complete the update of the moisture content of particle i, and calculate: M i =ρ p ·πR i 2 ; Where M i is the dry mass of particle i, and Δt is the time step.

8. The soil cracking simulation method based on two-dimensional particle discrete element method according to claim 1 is characterized in that: The linear strain method is used to reflect the influence of water content change on mechanics, as follows: The initial radius of particle i is R i0 , the water content changes by ΔW i , then its new radius R after water loss and shrinkage i Expressed as: R i =R i0 (1+a s ΔW i ); Where α s is the linear shrinkage coefficient; The distance between particles i and j in the normal direction is: The coordinates of the centers of the two particles remain unchanged. Only the effect of the change in water content is considered. The water content of particle i changes by ΔW within the calculation step. i , the volume of particle i shrinks, causing the distance between the two particles in the normal direction to change Δd n ,have to: Δd n =-R i0 a s ΔW i ; The increment of the normal force between the two particles ΔF n for: ΔF n =K n R i0 a s ΔW i ; F t is the maximum allowable contact tension, F c is the bonding contact force, is the internal friction angle, where: F t =2R·f t ; F c =2R·C。 9. The soil cracking simulation method based on two-dimensional particle discrete element method according to claim 1 is characterized in that: The flow reduction coefficient method is used to reflect the influence of mechanics on water diffusion, as follows: When particles i and j are in a bonded state, the spring breaks, the bond fails, and microcracks are generated. The water diffusion coefficient will decrease, and the water diffusion coefficient of the virtual channel will decrease. Calculated according to the following formula: Where C r is the flow reduction coefficient, 0≤C r ≤1, where k i and k j are the water diffusion coefficients of particles i and j, respectively.

Images