A method for simulating water diffusion in soil based on particle discrete element method
Through the method based on particle discrete elements, a virtual moisture conduction channel is constructed and the correction coefficient α is introduced, which solves the accuracy and applicability of moisture diffusion simulation in the soil in the prior art, and achieves a more efficient moisture diffusion simulation effect.
Patent Information
- Application Number
- CN202510271806.2
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2025-03-07
- Publication Date
- 2025-08-26
- Estimated Expiration
- 2045-03-07
AI Technical Summary
When simulating moisture diffusion in soil, the experimental method is expensive and difficult to monitor microscopic processes in real time. The numerical simulation method has limitations on applicability when dealing with heterogeneous media and complex boundary conditions, and the coupling of moisture diffusion models and discrete element models is insufficient.
Using a method based on particle discrete elements, a particle moisture diffusion model is established, a particle material parameters are assigned, a virtual moisture conduction channel is constructed, a water conduction between particles and between particle boundaries is calculated, and a correction coefficient α is introduced to adjust the cross-sectional width of the moisture channel is adjusted, and a specific moisture update iteration method is used for iterative calculation.
More precise water diffusion simulation is achieved, simulation accuracy and adaptability are improved, and can be suitable for water conduction simulation of different media, reducing calculation costs and resource waste.
Smart Images

Figure CN120180843B_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to the field of water conduction, and more particularly to a method for simulating water diffusion in soil based on particle discrete elements. Background Art
[0002] Water diffusion in soil is an important research topic in geotechnical engineering. Currently, related research mainly focuses on two aspects: experiments and numerical simulations.
[0003] Experimental methods are traditional means of studying water diffusion in soil, and can intuitively observe the migration patterns of water in the soil. However, experimental studies often focus on phenomenological observations and lack in-depth revelations of the underlying mechanisms. In addition, the experimental process is time-consuming and costly, making it difficult to monitor the microscopic process of water diffusion in real time. For example, Tang Chaosheng et al. studied the evolution of cracks under different temperatures, soil thicknesses, and numbers of dry-wet cycles, and proposed a series of quantitative indicators for the crack network on the soil surface. Although these studies provide rich phenomenological data, they are still insufficient in explaining the microscopic mechanisms of crack formation and development.
[0004] Numerical simulation methods are gaining widespread application in the study of water diffusion in soils. Continuum-based methods, such as the finite element method, meshless method, and finite volume method, have been widely used to simulate water diffusion processes. For example, Zhang Peisen et al. investigated the effects of different unloading paths on the damage characteristics and energy evolution of sandstone under stress-seepage coupling. While these methods can effectively address the physical phenomena involved in water diffusion, their applicability is limited when considering heterogeneous media or complex boundary conditions.
[0005] In recent years, methods based on discontinuous media have gradually gained attention. For example, the discrete element method (DEM) and particle dynamics methods can simulate the interactions between soil particles and their influence on water diffusion at the microscopic level. For example, Lin Zhuyuan et al. studied the boundary effects of soil shrinkage cracking using discrete element simulation. These methods have unique advantages in describing the discontinuity of water diffusion. However, there is still little research on coupling water diffusion models with discrete element models, and related work still has shortcomings in the physical meaning of the models and the selection of parameters. Summary of the Invention
[0006] In view of this, the present invention provides a soil water diffusion simulation method based on particle discrete element method to solve the problems existing in the background technology.
[0007] In order to achieve the above object, the present invention adopts the following technical solutions:
[0008] A method for simulating water diffusion in soil based on particle discrete element method includes the following steps:
[0009] Establish a particle moisture diffusion model and assign particle material parameters and control parameters. The material parameters include geometric parameters, mechanical parameters and moisture diffusion parameters.
[0010] By analyzing the contact relationship between particles and the contact between particles and boundaries, the coordination number of each particle is determined and a virtual water conduction channel is constructed;
[0011] Determine whether there is a moisture content gradient between each contact pair;
[0012] Calculate water conduction between particles and at particle boundaries;
[0013] Based on the water conductivity calculation results, the water content information of the particles is updated and iterative calculation is performed.
[0014] Optionally, the geometric parameters include particle position coordinate information, particle shape and particle size; the mechanical parameters include particle density ρ, normal stiffness K n , tangential stiffness K s , tensile strength σ t , bond strength c and internal friction angle Seepage parameters include water diffusion coefficient k, flow density q w .
[0015] Optionally, the particle water diffusion model is a discretized particle model that establishes polygons based on the actual shape and size of the model, and imports the polygon information into the MATLAB system as the model boundary, and establishes a model based on the polygons consisting of a finite number of randomly or regularly closely arranged discretized particles and walls.
[0016] Optional water transfer methods between particles and between particle boundaries are as follows:
[0017] In the two-dimensional case, water is conducted between particles i and j through a virtual uniform rectangle; the amount of water transferred from particle j to particle i per unit time is expressed as:
[0018]
[0019] Where W i and W j are the moisture contents of particles i and j respectively; S ij is the unmodified virtual water conduction width, which is assumed to be the average diameter of particles i and j; L ij is the distance between the centers of two particles; k ij is the water diffusion rate of the virtual uniform rectangle; S ij , L ij 、k ij The method for determining is as follows:
[0020]
[0021] Where (x i ,y i ) and (x j ,y j ) are the two-dimensional center coordinates of particles i and j, (x i ,y i , z i ) and (x j ,y j , z j ) are the three-dimensional center coordinates of particles i and j respectively; k i , k j are the water diffusion coefficients of particles i and j, respectively;
[0022] For a group of particles, the total flux entering particle i is the sum of all particles in its vicinity. Particle i is adjacent to n particles. The total flux entering particle i per unit time is expressed as:
[0023]
[0024] When moisture conduction occurs in a contact pair consisting of a particle and a boundary, particle i contacts boundary w, and the endpoint coordinates of the two-dimensional boundary are and The endpoint coordinates of the 3D boundary are For a fixed moisture content boundary W w , the flux from the boundary w to the particle i per unit time, denoted as Q w→i , the expression is as follows:
[0025]
[0026] Where S iw Cross-sectional area width, L iw is the vertical distance from the particle center to the boundary, k iw is the particle moisture diffusion coefficient, and β is the correction coefficient. The above parameters are determined by the following formula:
[0027]
[0028] k iw =k i ;
[0029]
[0030] Where n is the number of particles in contact with the boundary w;
[0031] For a fixed flow density q w The flux from the boundary w to the particle i per unit time is recorded as Q w→i, the expression is:
[0032] Q w→i =-βq w S iw ;
[0033] When particle i contacts m boundaries, the total flux flowing into particle i from these boundaries per unit time is expressed as:
[0034]
[0035] Optionally, a specific moisture update iteration method is used to calculate water conduction, as follows:
[0036] Particle i is adjacent to n particles and m boundaries, and the total water flow into particle i is expressed by the following formula:
[0037]
[0038] To calculate the total flux Q i The change of particle moisture content within a time step. The moisture content expression of particle i is as follows:
[0039]
[0040] Optionally, a correction coefficient α is introduced to adjust the cross-sectional width of the moisture channel; in two-dimensional space, the two-dimensional correction coefficient α 2D The expression is:
[0041]
[0042] In three-dimensional space, the three-dimensional correction coefficient α 3D The expression is:
[0043]
[0044] Where CN represents the coordination number, represents the two-dimensional volume fraction, represents the three-dimensional volume fraction.
[0045] It can be seen from the above technical solution that, compared with the prior art, the present invention provides a method for simulating soil water diffusion based on particle discrete elements, which has the following beneficial effects:
[0046] 1. A water content iteration method is proposed, which can calculate the change of water content within a specific time step to simulate water conduction more accurately;
[0047] 2. A correction factor α is introduced for the numerical simulation of the virtual water conduction channel. This correction factor can be adjusted accordingly based on the changes in the cross-section of the virtual water conduction channel, thereby ensuring the accuracy of the simulation. BRIEF DESCRIPTION OF THE DRAWINGS
[0048] 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.
[0049] Figure 1(a) is a schematic diagram of water conduction between particles and between particle boundaries through virtual channels in a two-dimensional case;
[0050] Figure 1(b) is a schematic diagram of water conduction between particles and between particle boundaries through virtual channels in a three-dimensional situation;
[0051] Figure 2(a) is a schematic diagram of the virtual water conduction channel network under dgap=0;
[0052] Figure 2(b) shows dgap = 0.3R min Schematic diagram of the virtual water conduction channel network below;
[0053] Figure 2(c) shows dgap = 0.7R min Schematic diagram of the virtual water conduction channel network below;
[0054] Figure 3 For a particle p represented by a regular pentagonal element 2D Schematic diagram;
[0055] Figure 4 Schematic diagram of particles and CN polygons representing each;
[0056] Figure 5 For particles p 3D With element e 3D Schematic diagram of CN side;
[0057] Figure 6 Schematic diagram of transient moisture diffusion simulation model;
[0058] Figure 7 The distribution diagram of water content in soil strips at different time points;
[0059] Figure 8 It is a schematic diagram of the overall process of the present invention. DETAILED DESCRIPTION
[0060] 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.
[0061] The embodiment of the present invention discloses a method for simulating water diffusion in soil based on particle discrete elements, such as Figure 8 As shown, the following steps are included:
[0062] Step 1: Establish a particle moisture diffusion model (PMDM) and assign particle material parameters and control parameters, where the material parameters include geometric parameters, mechanical parameters, and moisture diffusion parameters. Geometric parameters include particle position coordinate information, particle shape, and particle size; mechanical parameters include particle density ρ, normal stiffness K n , tangential stiffness K s , tensile strength σ t , bond strength c and internal friction angle Seepage parameters include water diffusion coefficient k, flow density q w wait.
[0063] Step 2: By analyzing the contact relationships between particles and the contact between particles and boundaries, the coordination number (CN) of each particle is determined and a virtual water conduction channel is constructed based on this. This step is crucial for studying the interaction mechanism between particles and optimizing the water transfer path.
[0064] Step 3: Determine whether there is a moisture content gradient between each contact pair.
[0065] Step 4: Calculate the water conduction between particles and at the boundaries of particles.
[0066] Step 5: Update the moisture content and other related information of the particles according to the results calculated in step 4 and perform iterative calculations.
[0067] Furthermore, the particle diffusion model (PMDM) establishes a discretized particle model by creating a polygon based on the actual shape and size of the model. This polygon information is imported into the MATLAB system as the model boundary. Based on this polygon, a model consisting of a finite number of randomly or regularly densely arranged discretized particles and walls is established. This model uses the following assumptions:
[0068] 1) In particle-based discontinuous deformation analysis (PDDA), the solution domain is discretized into densely packed particles that mimic the behavior of a continuous medium by bonding together. In the particle water diffusion model (PMDM), water diffusion is assumed to occur within adjacent particles. In addition, water conduction only occurs between contact pairs (particle-particle, particle-boundary), indicating that water diffusion depends on the contact state of the contact pair. For example, if the contact state is open, water transfer is reduced to account for the effects of discontinuities such as cracks on water diffusion.
[0069] 2) Water conduction between particles and at particle boundaries occurs through virtual channels as shown in Figure 1. Virtual channels between particles are square in 2D and cube in 3D; virtual channels between particle boundaries are rectangular in 2D and cuboid in 3D.
[0070] 3) The water transfer between particles is proportional to the moisture content gradient, as shown in formula (1).
[0071]
[0072] Where q is the water flux per unit area, k is the water diffusion coefficient, is the moisture content gradient.
[0073] Based on the total flux into a given mass, the change in water content can be calculated (2) definition.
[0074]
[0075] Where Q is the total water flux, t is time, and M is mass
[0076] Furthermore, the water transfer between particles and between particle boundaries is as follows:
[0077] As shown in Figure 1(a), taking particles i and j as an example, in the two-dimensional case, it is assumed that water is conducted through a virtual uniform rectangle. The amount of water transferred from particle j to particle i per unit time is expressed as .
[0078]
[0079] Where W i and W j are the moisture contents of particles i and j respectively; S ij is the unmodified virtual water conduction width, whose value is assumed to be the average diameter of particles i and j; L ij is the distance between the centers of two particles; k ij is the water diffusion rate of the virtual uniform rectangle. Sij , L ij 、k ij The method for determining is as follows:
[0080]
[0081] Where (x i ,y i ) and (x j ,y j ) are the two-dimensional center coordinates of particles i and j, (x i ,y i , z i ) and (x j ,y j , z j ) are the three-dimensional center coordinates of particles i and j respectively. i , k j are the water diffusion coefficients of particles i and j, respectively.
[0082] For a group of particles, the total flux entering particle i is the sum of all the particles in its vicinity. Assuming that particle i is adjacent to n particles, the total flux entering particle i per unit time can be expressed as:
[0083]
[0084] Moisture conduction can also occur in a contact pair consisting of a particle and a boundary, as shown in Figure 1(b). Assume that particle i is in contact with boundary w. The endpoint coordinates of the two-dimensional boundary are and The endpoint coordinates of the 3D boundary are For a fixed moisture content boundary W w , the flux from the boundary w to the particle i per unit time, denoted as Q w→i , the expression is as follows
[0085]
[0086] Where S iw Cross-sectional area width, L iw is the vertical distance from the center of the particle to the boundary, k iw is the particle moisture diffusion coefficient, and β is the correction coefficient. The above parameters are determined by the following formula
[0087]
[0088] k iw =k i (11)
[0089]
[0090] Where n is the number of particles in contact with the boundary w.
[0091] For a fixed flow density q w The flux from the boundary w to the particle i per unit time is recorded as Q w→i , the expression is:
[0092] Q w→i =-βq w S iw (13)
[0093] If particle i contacts m boundaries, the total flux from these boundaries into particle i per unit time can be expressed as:
[0094]
[0095] Furthermore, a specific moisture update iteration method is used to calculate water conduction, as follows:
[0096] Once the total water flow rate of a particle is calculated, its water content change can be calculated. A particle can be in contact with multiple particles and boundaries simultaneously. Therefore, the total water flow rate flowing into particle i per unit time is the sum of the water flow rates flowing into particle i from all adjacent particles and adjacent boundaries. Assuming that particle i is adjacent to n particles and m boundaries, the total water flow rate flowing into particle i can be expressed as follows:
[0097]
[0098] To calculate the total flux Q i The change of particle moisture content within a time step. The moisture content expression of particle i is as follows:
[0099]
[0100] Furthermore, this embodiment introduces an innovative correction coefficient α, which is as follows:
[0101] In order to satisfy the mass conservation in DEM calculation, the particle density needs to be adjusted. The expression is as follows:
[0102]
[0103] Where: p is the particle density, ρ is the density of the continuous medium, f v is the particle volume fraction.
[0104] Because the transmission surface of the virtual water channel is artificially defined, it cannot represent the actual water flow surface. Calibration is typically used to correct the water conduction surface, but due to the complexity of the calibration process, an analytical expression for α is derived here based on the area equivalence principle. First, a brief introduction to the involved microscopic parameters and coordination numbers is provided.
[0105] α is a geometric correction factor used to adjust the cross-sectional width of the moisture channel to accurately reproduce the macroscopic moisture diffusion coefficient. In order to derive the theoretical formula of α, the concept of coordination number (CN) is introduced. The coordination number (CN) represents the average number of neighboring particles for each particle in a particle system. Under normal circumstances, it is difficult to control because only particles that are in contact with each other are considered. Therefore, dgap is introduced to solve this problem. dgap is defined as the product of the minimum particle size of the particle and a dimensionless number greater than 0. Its function is to regard the particles as forming a contact pair only when the distance between them is less than dgap. When dgap = 0, the coordination number (CN) of particle i in Figure 2 (a) is 4, and when 0.3R min ≤dgap≤0.7R min When , the coordination numbers of particle i in Figure 2(b) and Figure 2(c) are 5 and 6, respectively. As mentioned above, water conduction occurs within the contact pair, so the number of virtual water channels in the particle system depends on the coordination number (CN), and the coordination number (CN) is positively correlated with dgap.
[0106] The introduction of the correction factor α is one of the core innovations of this technology. It can flexibly adjust to different cross-sectional characteristics, thereby significantly improving the adaptability and accuracy of the system. This feature makes α show significant advantages in many aspects:
[0107] 1) By adjusting α, the system can automatically optimize the calculation model according to the geometric characteristics and flow velocity distribution of different cross sections, thereby more accurately reflecting the actual working conditions.
[0108] 2)α can effectively correct the error caused by cross-sectional non-uniformity, and the correction effect is more significant.
[0109] 3)α can optimize the calculation model, reduce resource waste and repeated calculations caused by errors, thereby reducing overall operating costs and improving calculation efficiency.
[0110] The proposed water content iteration method can accurately calculate the water content change within a specific time step, thereby more accurately simulating the water conduction process:
[0111] 1) By iteratively calculating the change in water content, this method can more accurately reflect the dynamic changes in the water conduction process and provide higher simulation accuracy
[0112] 2) The iterative method can be flexibly adjusted according to different time steps and initial conditions, and is suitable for water conduction simulation in a variety of media (such as soil, frozen soil, etc.), with better adaptability.
[0113] like Figure 3 In the two-dimensional case shown, for the 2D and a group of particles with coordination number CN, a typical particle p 2D (The radius is R 2D ) and a circumscribed regular CN polygon element e 2D Under this assumption, the entire continuum is divided into regular CN polygons with the same number of particles. By establishing an equivalent relationship between polygon elements and particles, the characteristic parameters of the particle system can be corrected. It is worth noting that, based on statistical averages, the coordination number (CN) of the CN polygon can be integer or non-integer.
[0114] For polygon element e 2D The cross-sectional width of the water channel between it and the surrounding elements is the side length of the polygon For the particle p 2D The uncorrected width of the water channel between it and the adjacent particles is given by Equation 7 as 2R 2D The equivalent relationship between the section widths is established as follows:
[0115]
[0116] According to the definition of volume fraction
[0117]
[0118] By combining the formula and the formula, the expression of α can be obtained as:
[0119]
[0120] It can be seen from the above formula that the correction coefficient α 2D Depends on the coordination number CN and volume fraction
[0121] Let's verify the validity of this expression for modeling some rules, such as Figure 4 As shown. Taking the hexagonal packing of CN=6 as an example, the side length of the hexagonal element is Volume fraction for α 2D Can be considered as The area of the virtual water conduction surface of the hexagonal element is The corrected particle transport surface area is calculated The equivalence of the areas of the two regions shows that the proposed expression is applicable to hexagonal packing. Figure 4 The other cases shown can be verified similarly.
[0122] In two-dimensional space, any continuous region can be approximated relatively easily by a set of regular polygons (such as triangles, squares, or hexagons). However, filling space with regular polyhedra in three dimensions is more complex; these polyhedra may have a varying number of faces, and each face may have a varying number of edges.
[0123] Let the representative element of the polyhedron be e 3D , it has CN faces, each with a side length of A regular polygon. Figure 5 Shows the particles and their representative elements p 3D A face and element p 3D The face of is an equilateral triangle.
[0124] The correction coefficient α in the three-dimensional case can be derived by referring to the similar method in the two-dimensional case 3D This process involves adjusting the geometric properties of the representative elements (such as side length and area) to more accurately reflect the actual particle characteristics in three-dimensional space. Therefore, even in the face of the difficulty of filling regular polyhedrons in three-dimensional space, it is still possible to find an effective method to describe and derive the correction factor α 3D .
[0125] Let the area of each face be element e 3D The volume is It can be calculated as follows,
[0126]
[0127] The volume fraction of three-dimensional particle accumulation is calculated as follows:
[0128]
[0129] According to formula (4), particle p 3D The average value of the water-conducting surface area of the virtual water-conducting channel without correction is written as,
[0130]
[0131] element e 3D The transmission surface area is CN times the polygon area,
[0132]
[0133] According to the area equivalence criterion of the transmission surface, we have
[0134]
[0135] Where α 3D is the three-dimensional correction coefficient, substitute formula (22) into formula (25), α 3D Can write,
[0136]
[0137] In order to further verify the effectiveness of the present invention, the following experiments were conducted: The numerical experiments aim to verify the effectiveness and accuracy of the proposed PMDM model in simulating transient moisture diffusion problems by comparing the simulation results with the analytical solutions. Figure 6 As shown in the figure, the research object is a rectangular soil strip model with a height of 0.01m and a length of 0.1m. The model contains 29784 particles, the coordination number is 5.2305, the initial moisture content is set to 0%, and the flux density is 0.00001kg / (m 2 s), and the other boundaries are impervious boundaries.
[0138] By selecting specific parameters, the diffusion coefficient k = 1×10 -6 kg / (m·s), ρ0=1333kg / m, the distribution of water content in the model at different time points is as follows Figure 7 shown.
[0139] Figure 7 (Top) shows the initial water content distribution. The darker blue portion on the right indicates lower water content. Over time, the color on the left gradually transitions from blue to red, indicating a gradual increase in water content. Figure 7 The image (bottom) shows water diffusion into the middle and right sides of the soil strip, where the color change is more pronounced, indicating higher water content. These images are arranged in chronological order, from top to bottom, showing the diffusion of water from the left boundary to the right over time. This allows for a visual understanding of how water diffuses through soil conditions over time and how water content varies at different locations. This visualization, consistent with actual conditions, facilitates understanding and analysis of transient water diffusion processes.
[0140] 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.
[0141] 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 method for simulating water diffusion in soil based on particle discrete element method, characterized in that: The following steps are involved: Establish a particle moisture diffusion model and assign particle material parameters and control parameters. The material parameters include geometric parameters, mechanical parameters and moisture diffusion parameters. By analyzing the contact relationship between particles and the contact between particles and boundaries, the coordination number of each particle is determined and a virtual water conduction channel is constructed; Determine whether there is a moisture content gradient between each contact pair; Calculate water conduction between particles and at particle boundaries; Update the information related to the moisture content of the particles based on the water conductivity calculation results and perform iterative calculations; The water transfer between particles and at particle boundaries is as follows: In the two-dimensional case, water is conducted between particles i and j through a virtual uniform rectangle; the amount of water transferred from particle j to particle i per unit time is expressed as: Where W i and W j are the moisture contents of particles i and j respectively; S ij is the unmodified virtual water conduction width, whose value is assumed to be the average diameter of particles i and j; L ij is the distance between the centers of two particles; k ij is the water diffusion rate of the virtual uniform rectangle; S ij 、L ij 、k ij The method for determining is as follows: Where (x i ,y i ) and (x j ,y j ) are the two-dimensional center coordinates of particles i and j, (x i ,y i , z i ) and (x j ,y j , z j ) are the three-dimensional center coordinates of particles i and j respectively; k i , k j are the water diffusion coefficients of particles i and j, respectively; For a group of particles, the total flux entering particle i is the sum of all particles in its vicinity. Particle i is adjacent to n particles. The total flux entering particle i per unit time is expressed as: When moisture conduction occurs in a contact pair consisting of a particle and a boundary, particle i contacts boundary w, and the endpoint coordinates of the two-dimensional boundary are and The endpoint coordinates of the 3D boundary are For a fixed moisture content boundary W w , the flux from the boundary w to the particle i per unit time, denoted as Q w→i , the expression is as follows: Where S iw Cross-sectional area width, L iw is the vertical distance from the particle center to the boundary, k iw is the particle moisture diffusion coefficient, and β is the correction coefficient. The above parameters are determined by the following formula: k iw =k i ; Where n is the number of particles in contact with the boundary w; For a fixed flow density q w The flux from the boundary w to the particle i per unit time is recorded as Q w→i , the expression is: Q w→i =-βq w S iw ; When particle i contacts m boundaries, the total flux flowing into particle i from these boundaries per unit time is expressed as: A specific moisture update iteration method is used to calculate water conduction, as follows: Particle i is adjacent to n particles and m boundaries, and the total water flow into particle i is expressed by the following formula: To calculate the total flux Q i The change of particle moisture content within a time step. The moisture content expression of particle i is as follows:
2. The method for simulating soil water diffusion based on particle discrete element method according to claim 1, characterized in that: Geometric parameters include particle position coordinate information, particle shape and particle size; mechanical parameters include particle density ρ, normal stiffness K n , tangential stiffness K s , tensile strength σ t , bond strength c and internal friction angle Seepage parameters include water diffusion coefficient k, flow density q w .
3. The method for simulating soil water diffusion based on particle discrete element method according to claim 1 is characterized in that: The particle water diffusion model is a discretized particle model that establishes polygons based on the actual shape and size of the model, imports the polygon information into the MATLAB system as the model boundary, and establishes a model composed of a finite number of randomly or regularly closely arranged discretized particles and walls based on the polygons.
4. The method for simulating soil water diffusion based on particle discrete element method according to claim 1, characterized in that: It also includes the introduction of a correction coefficient α to adjust the cross-sectional width of the moisture channel; in two-dimensional space, the two-dimensional correction coefficient α 2D The expression is: In three-dimensional space, the three-dimensional correction coefficient α 3D The expression is: Where CN represents the coordination number, represents the two-dimensional volume fraction, represents the three-dimensional volume fraction.