Water diffusion-mechanical coupling calculation method based on particle DDA

Abstract

The invention discloses a moisture diffusion-mechanical coupling calculation method based on particle DDA, and the method specifically comprises the steps: carrying out the discretization of a to-be-solved region through the particle DDA, and endowing a control parameter; contact retrieval is carried out between the particles and between the particles and the boundary, the action ranges of mechanics and moisture diffusion are determined, and a force chain and a virtual moisture diffusion channel are established between the action ranges; calculating the mechanical part of the DDA of the particles, and updating the position information of the particles; assembling a moisture diffusion equation set of the particle DDA, solving the equation set and updating particle moisture content information. According to the moisture diffusion-mechanical coupling calculation method based on the particle DDA, an implicit solution strategy is adopted, the time step is large, and calculation is stable; and the influence of cracks on the moisture diffusion process can be fully considered. Benefited from the improvement, the constructed water diffusion model not only improves the calculation efficiency and simplifies the material parameter assignment process, but also can more accurately describe the water migration behavior in the crack-containing structure.

Description

Technical Field

[0001] This invention belongs to the field of water diffusion research technology, specifically involving a water diffusion-mechanical coupling calculation method based on particle DDA. Background Technology

[0002] Currently, laboratory experiments remain a common method for investigating the diffusion patterns of moisture in soil. However, these experimental studies often focus on phenomenological observations and struggle to reveal the underlying physical mechanisms. Furthermore, both sample preparation and laboratory or field testing are typically costly and time-consuming. To overcome these limitations, researchers have developed various numerical simulation methods to characterize the moisture diffusion process in soil. Numerical methods based on continuous medium theory are widely used when dealing with moisture diffusion problems within a medium, such as the finite element method, meshless method, finite volume method, and numerical manifold method.

[0003] In the coupling problem of soil moisture diffusion and mechanical behavior, there is a strong bidirectional feedback between moisture migration and soil mechanical deformation. Moisture infiltration reduces soil matrix suction, softens the structure, and induces settlement or instability; while soil compression and shear deformation change pore connectivity, inversely affecting the moisture diffusion path and rate. This poses a significant bottleneck to traditional continuous numerical methods (such as the finite element method and the finite difference method). The continuous medium assumption on which they are based makes it difficult to accurately characterize the discrete pore structure and particle contact relationships inside the soil. When capturing microscopic moisture diffusion characteristics such as preferential flow and fissure seepage, as well as nonlinear mechanical responses such as plastic deformation and local failure, the numerical stability and computational accuracy are greatly limited, making it difficult to truly reflect the coupling mechanism.

[0004] In contrast, discontinuous methods such as the Discrete Element Method (DEM), Discontinuous Deformation Analysis (DDA), and particle dynamics (e.g., SPH) can reveal the coupled nature of moisture diffusion and mechanical behavior at the microscale by directly simulating the motion, contact, and interaction of discrete particles or elements. For example, the DEM can simultaneously simulate interparticle mechanical contact and pore moisture migration, discontinuous deformation analysis excels at handling the coupling of large block displacement and crack seepage, and particle dynamics is applicable to the interaction between free surface flow and soil deformation. However, existing research mostly focuses on preliminary explorations of single methods, concentrating on basic simulations under simple boundary conditions, with very few coupled studies on complex engineering scenarios such as unsaturated soils and layered soils. Coupled models do not adequately consider key hydraulic-mechanical properties such as the hysteresis effect of soil-water characteristic curves and suction-dependent shear strength, resulting in limited physical realism and prediction accuracy. At the same time, related work lacks systematic comparative verification with indoor test and field monitoring data, the sensitivity analysis of model parameters and the definition of applicable scope are unclear, and the computational efficiency issues of discontinuous methods themselves further limit their engineering application. The overall depth and completeness of research urgently need to be improved. Summary of the Invention

[0005] The purpose of this invention is to provide a moisture diffusion-mechanical coupling calculation method based on particle DDA, which not only improves calculation efficiency and simplifies the material parameter assignment process, but also more accurately describes the moisture migration behavior inside fractured structures.

[0006] The technical solution adopted in this invention is a moisture diffusion-mechanical coupling calculation method based on particle DDA, which is implemented according to the following steps: Step 1: Discretize the region to be solved using particle DDA and assign control parameters; Step 2: Conduct contact retrieval between particles and between particles and boundaries to determine the domain of action of mechanical and moisture diffusion, and establish force chains and virtual moisture diffusion channels between them. Step 3: Perform mechanical part calculations for particle DDA and update particle position information; Step 4: Assemble the moisture diffusion equations for particle DDA, solve the equations and update the particle moisture content information.

[0007] The invention is further characterized in that, In step 1, the control parameters include geometric parameters, mechanical parameters, and moisture diffusion parameters; the geometric parameters include particle position coordinates and particle size; the mechanical parameters include particle density. Normal contact stiffness Tangential contact stiffness ,tensile strength Bond strength internal friction angle With shrinkage coefficient Moisture diffusion parameters include the moisture diffusion coefficient. With boundary flow mass density .

[0008] Step 2 specifically involves: Step 2.1: Construct a unified particle contact network; Step 2.2: Traverse each pair in the particle contact network and determine whether the center distance of the pair is less than or equal to the sum of the particle diameters of the two particles; if the mechanical contact criterion is satisfied, record the pair and its contact geometry information to form a force chain list; Step 2.3: Traverse each pair in the particle contact network; determine whether the center distance of the pair is less than or equal to the preset hydraulic search radius; if the hydraulic diffusion criterion is satisfied, record the pair and its channel geometry information to form a list of hydraulic channels; for particle units adjacent to the boundary of the simulation domain, apply the mechanical contact criterion and the hydraulic diffusion criterion respectively according to their distance from the boundary to establish force chains or hydraulic channels between particles and the boundary, and add them to the corresponding lists respectively.

[0009] Step 3 specifically involves: Step 3.1, Force initialization and application of volume forces; Iterate through each computational particle in the system, initialize the resultant force and resultant torque acting on it to zero; calculate and apply a preset volume force based on the mass of each computational particle, and add the force to the resultant force of that particle. Step 3.2, Calculation of contact force based on force chain: Iterate through the force chain list generated in Step 2.2; for each force chain data item in the list, calculate the normal force and tangential force, and update the force and moment. Step 3.3, Motion State Update: After completing the calculation of all forces, traverse all computational particle elements; based on the current net force and net torque acting on it, combined with its own mass and moment of inertia, solve for its linear acceleration and angular acceleration; use the time integration algorithm to update the velocity, angular velocity, position and angle information of each computational particle element; apply position constraints to computational particle elements that exceed the boundary of the simulation domain after updating their positions to ensure that they always move within the simulation domain.

[0010] In step 4, the specific mechanisms of moisture diffusion between particles and between particle boundaries are as follows: Let the particles The radius is Let the particles The radius is The transmission surface of the virtual channel for moisture diffusion is a rectangle. The width of the rectangle is defined as the side length being the arithmetic mean of the diameters of the two particles, that is, the sum of the radii of the two particles. Define the water diffusion coefficient of the virtual water diffusion channel. It is the harmonic mean of the water diffusion coefficients of the two particles, as shown in formula (3); (3) in, and Particles With particles The water diffusion coefficient; Assume that water is conducted through a virtual uniform rectangular channel between two particles or between a particle and a wall; as shown in formula (5); (5); in, and Particles With particles Moisture content; A correction coefficient is introduced to correct the width of the virtual water diffusion channel, as shown in formula (6); (6) Among them, coordination number It is a measure of the average number of neighboring particles for each particle. This represents the volume fraction of a particle aggregate. From the boundary per unit time Flow towards particles Traffic As shown in formula (7); (7); In the formula, length It is the distance from the center of the particle to the wall, the width of the virtual moisture diffusion channel. and correction factor As shown in formulas (8) and (9); (8) (9) Let the particles have Adjacent particles and Adjacent boundaries, inflowing particles i The total flow includes contributions from adjacent particles and walls, as shown in Equation (10); (10) By substituting formula (10) into formula (5), we obtain formula (11); (11) in, Represents the particles of the previous time step i moisture content, It is a wall w Temperature, Δ t It is the time step; ; Equation (11) can be further written as equation (12); (12) Suppose the particulate system consists of n Composed of individual particles, equation (12) can be written in matrix form, as shown in formula (13); AW = B (13) in, ; ; ; ; AIt is a coefficient matrix that describes the interaction relationships of water diffusion between particles; a ii : reflects the first i The moisture retention capacity of each particle itself; a ij : Reflecting the first j The particle pairs the first i The water diffusion and transfer effect of individual particles; W It's the water content, for each element. W i Corresponding to the i The current moisture content of each particle; B It is the right-hand term vector, determined by the initial water content of the particles and the boundary / source terms; Formula (13) is the final set of water diffusion equations. The water content of the particles can be updated by solving formula (13). Iterate through all particles and update their water content with the new value.

[0011] The beneficial effects of this invention are: The present invention provides a particle-based DDA-based coupled calculation method for water diffusion and mechanics, employing an implicit solution strategy, a large time step, and computational stability; it also fully considers the influence of cracks on the water diffusion process. Thanks to these improvements, the constructed water diffusion model not only enhances computational efficiency and simplifies the material parameter assignment process, but also more accurately describes the water migration behavior within cracked structures. Attached Figure Description

[0012] Figure 1 This is a schematic diagram of the bonded particle model under the normal spring; Figure 2 This is a schematic diagram of the adhesive particle model under the tangential spring. Figure 3 Schematic diagram of moisture diffusion between particles via virtual channels Figure 4 This is a schematic diagram of moisture diffusion between particle boundaries through virtual channels. Detailed Implementation

[0013] The present invention will now be described in detail with reference to specific embodiments and accompanying drawings.

[0014] Example 1 This invention is based on a moisture diffusion-mechanical coupling calculation method using particle DDA, and is implemented according to the following steps: Step 1: Discretize the region to be solved using particle DDA and assign control parameters; Step 2: Conduct contact retrieval between particles and between particles and boundaries to determine the domain of action of mechanical and moisture diffusion, and establish force chains and virtual moisture diffusion channels between them. Step 3: Perform mechanical part calculations for particle DDA and update particle position information; Step 4: Assemble the moisture diffusion equations for particle DDA, solve the equations and update the particle moisture content information.

[0015] Example 2 This invention is based on a moisture diffusion-mechanical coupling calculation method using particle DDA, and is implemented according to the following steps: Step 1: Discretize the region to be solved using particle DDA and assign control parameters; The particle shape adopts round / spherical particles to simplify contact calculations; the particle size range is designed according to the regional scale and computational efficiency, and a log-normal distribution or uniform distribution is adopted to ensure rational particle-level matching.

[0016] Particle arrangement method: Particles are arranged according to a preset grid to simulate regular structures such as bedding and joints. The particle spacing is controlled to be 1.0 to 1.2 times the particle size to ensure effective contact. Boundary adaptation arrangement: For irregular boundary regions, a boundary particle refinement strategy is adopted, and the particle size near the boundary is reduced to 1 / 2 to 1 / 3 of that of the internal particles. Local densification ensures the accuracy of boundary shape fitting and avoids calculation errors caused by boundary effects.

[0017] Particle coverage: The volume / area ratio occupied by particles in the region to be solved is ≥90%; Contact effectiveness: The average number of contact points per particle is controlled to be 3~6 (two-dimensional / 6~12 (three-dimensional) to ensure the continuity of force transmission; Convergence calculation: The stability of the discretization results is verified through particle size sensitivity analysis; Control parameters include geometric parameters, mechanical parameters, and moisture diffusion parameters; geometric parameters include particle location coordinates and particle size; mechanical parameters include particle density. Normal contact stiffness Tangential contact stiffness ,tensile strength Bond strength internal friction angle With shrinkage coefficient Moisture diffusion parameters include the moisture diffusion coefficient. With boundary flow mass density ; Since the region to be solved is discretized by particles, and gaps inevitably exist between the particles, the total mass of the region changes, necessitating an adjustment of the particle density. Compensation adjustments are made to ensure that the total mass of the region to be solved remains constant, and the compensated particle density... for: ;in, For the density of a continuous medium, This represents the volume fraction of the particulate system.

[0018] Step 2: Conduct contact retrieval between particles and between particles and boundaries to determine the domain of action of mechanical and moisture diffusion, and establish force chains and virtual moisture diffusion channels between them. In particle DDA modeling of continuous media, contact forces and moisture diffusion only occur between adjacent particles and between adjacent particles and boundaries. Therefore, the particle contact mesh can be considered as the domain of action of mechanics and moisture diffusion; specifically: Step 2.1: Construct a unified particle contact network; Acquire the position and size information of all computational particle units at the current moment, as well as the geometric boundary information of the simulation domain. Based on the maximum particle size of all computational particle units in the system, set the element edge length of a background grid, ensuring that this edge length is greater than the maximum particle size. Construct a spatial index grid defined by the element edge length on the simulation domain. Traverse all computational particle units and calculate the grid element index of each unit based on its position coordinates. Map the unique identifier of each computational particle unit to its corresponding grid element, forming a grid-particle index relationship. Traverse every grid element in the spatial index grid. For each target particle unit within the current grid element, perform pairing detection within that unit and all its neighboring units. For each potential pairing, calculate the center distance between them. Determine if this center distance is less than or equal to a preset neighborhood search radius. The neighborhood search radius should not be less than the maximum particle size to ensure coverage of all potential interaction pairs. If the condition is met, record the pairing relationship, forming a particle contact network containing all spatially adjacent particle pairs. This network is the common domain of mechanical and moisture diffusion.

[0019] Step 2.2: Mechanical action domain screening and force chain list generation: Traverse each pair in the particle contact network and determine whether the center distance of the pair is less than or equal to the sum of the particle diameters of the two particles.

[0020] If the mechanical contact criterion is met, the pairing and its contact geometry information (including contact normal vector and normal overlap) are recorded to form a force chain list. Step 2.3: Moisture diffusion domain screening and hydraulic channel list generation: Traverse each pair in the particle contact network; determine whether the center distance of the pair is less than or equal to the preset hydraulic search radius. If the hydraulic diffusion criterion is met, record the pair and its channel geometry information (including the center distance) to form a hydraulic channel list.

[0021] Similarly, for particle units adjacent to the boundary of the simulation domain, force chains or hydraulic channels between the particles and the boundary are established by applying the mechanical contact criterion and the hydraulic diffusion criterion respectively, based on their distance from the boundary, and added to the corresponding lists respectively.

[0022] Step 3: Perform mechanical part calculations for particle DDA and update particle position information; Step 3.1, Force initialization and application of volume forces; Iterate through each computational particle in the system, initializing its net force and net torque to zero. Calculate and apply a preset volume force based on the mass of each computational particle, and add this force to the net force of that particle.

[0023] Step 3.2, Contact force calculation based on force chains: The force chain list generated in Step 2.2 is traversed and calculated; for each force chain data item in the list: Normal force calculation: Based on the preset normal contact stiffness, normal overlap and normal relative velocity, calculate the normal contact force between the pairs.

[0024] Tangential force calculation: Based on the preset tangential contact stiffness, tangential relative displacement increment and internal friction angle, combined with the Coulomb friction criterion, the tangential contact force between the pairs is calculated.

[0025] normal force of particles With tangential force All follow Hooke's Law. ;in, , These are the normal contact stiffness and the tangential contact stiffness, respectively. , The negative sign represents the change in the relative displacement of the contact in the tangential and normal directions over a time step, indicating that the spring force opposes the increase in relative displacement.

[0026] Force and torque update: The calculated normal force and tangential force are combined into a total contact force, and according to Newton's third law, the total contact force and the torque it generates are applied to the resultant force and resultant torque of the calculated particle units at both ends of the force chain.

[0027] Step 3.3, Motion State Update: After completing the calculation of all forces, traverse all computational particle elements. Based on the current net force and net torque acting on it, combined with its own mass and moment of inertia, solve for its linear acceleration and angular acceleration. Use a preset time integration algorithm to update the velocity, angular velocity, position, and angle information of each computational particle element. Apply position constraints to computational particle elements that exceed the boundary of the simulation domain after the position update to ensure that they always move within the simulation domain.

[0028] Step 4: Assemble the moisture diffusion equations for particle DDA, solve the equations and update the particle moisture content information; Example 3 Furthermore, step 4 specifically involves: The particles are connected by a set of normal and tangential springs. Forces act on their respective centers; the normal force runs along the line connecting the centers, while the tangential force is perpendicular to the normal force (e.g., ...). Figure 1 , Figure 2 (As shown).

[0029] Moisture diffusion between particles and between particle boundaries occurs through methods such as Figure 3 , Figure 4 The virtual channels shown are used for diffusion. The virtual moisture diffusion channels between particles and between particle boundaries are rectangular.

[0030] The specific mechanisms of moisture diffusion between particles and between particle boundaries are as follows: Let the particles The radius is Let the particles The radius is The transmission surface of the virtual moisture diffusion channel is a rectangle, and the width of the rectangle is defined as the side length being the arithmetic mean of the diameters of the two particles, i.e., the sum of the radii of the two particles. Based on the above assumptions, the cross-sectional width of the uncorrected virtual moisture diffusion channel is... As shown in formula (1); (1) Distance between the centers of the two particles As shown in formula (2); (2) Among them, particles The center coordinates are granules The center coordinates are ; Define the water diffusion coefficient of the virtual water diffusion channel. It is the harmonic mean of the water diffusion coefficients of the two particles, as shown in formula (3); (3) in, and Particles With particles The water diffusion coefficient.

[0031] Assume that water is conducted through a virtual uniform rectangular channel between two particles or between a particle and a wall. According to formula (4), the amount of water transferred per unit time by the particles can be calculated. Flow towards particles Traffic As shown in formula (5); In a continuous medium, the diffusion of moisture should follow the principle that the flow rate exchanged per unit time is proportional to the water content gradient, as shown in formula (4); (4) in, The amount of traffic exchanged per unit of time. The water diffusion coefficient is... This represents the moisture content gradient. The negative sign clearly indicates that the driving direction of moisture diffusion is from the high moisture content region to the low moisture content region.

[0032] (5); in, and Particles With particles Moisture content; A correction factor is introduced to adjust for the width of the virtual water diffusion channel; an appropriate The macroscopic moisture diffusion coefficient can be accurately reproduced, as shown in formula (6); (6) Among them, coordination number It is a measure of the average number of neighboring particles for each particle. This represents the volume fraction of the particle aggregate.

[0033] Moisture diffusion at particle boundaries is essentially similar to moisture diffusion between particles. For example... Figure 1 As shown, the boundary In the two-dimensional case, this is simplified to a line segment. Once a particle comes into contact with the wall, a virtual, uniform rectangular water diffusion channel is formed between the particle and the wall. Let the particle... With boundary Contact, boundary A flow mass density is applied Its value is usually negative, representing the loss of water per unit time. Therefore, the value of water lost per unit time from the boundary... Flow towards particles Traffic As shown in formula (7); (7); In the formula, length It is the distance from the center of the particle to the wall, the width of the virtual moisture diffusion channel. and correction factor As shown in formulas (8) and (9); (8) (9) In the formula, i =1,2… m , Is with the boundary The number of particles in contact; the center coordinates of a particle in contact with the boundary ( , ), the center coordinates of another particle in contact with the boundary ( , ); Inflow particles The flow rate is the flow rate of all adjacent particles and the inflow particles at the adjacent boundary. The total flow rate. Assume the particle... have Adjacent particles and Adjacent boundaries, inflowing particles i The total flow includes contributions from adjacent particles and walls, as shown in Equation (10); (10) By substituting formula (10) into formula (5), we obtain formula (11); (11) in, Represents the particles of the previous time step i moisture content, It is a wall w Temperature, Δ t It is the time step; ; If only the right side of the equation is... and If the water content of the previous time step is replaced by the water content of the previous time step, then equation (11) can be solved explicitly. However, the time step in the explicit method must be very small, which significantly affects computational efficiency. In addition, implicit solutions are a major feature of DDA compared to DEM. Therefore, the implicit solution form of equation (11) is derived. Considering W i and W j It is the unknown quantity to be solved. Equation (11) can be further written as equation (12). (12) Suppose the particulate system consists of n Composed of individual particles, equation (12) can be written in matrix form, as shown in formula (13); AW = B (13) in, ; ; ; ; A It is the coefficient matrix ( n * n dimension, n (This refers to the total number of particles), describing the interaction between particles in terms of water diffusion. a ii (Diagonal element): embodies the first i The moisture retention capacity of each particle itself includes: an intrinsic term; and a term representing the resistance to moisture diffusion between particles (within time step Δ). t diffusion coefficient k ij Particle spacing / channel parameters S ij Related); the effect of moisture and mechanical coupling (with parameters) Interparticle mechanical interaction parameters , (Related).

[0034] a ij (Non-diagonal element): Reflects the first j The particle pairs the first i The water diffusion and transfer effect of individual particles only occurs when j yes i The value is non-zero when the particles are within the diffusion domain, otherwise it is 0; its value is determined by parameters such as diffusion resistance and time step, reflecting the strength of moisture transfer between particles.

[0035] W It's the water content, for each element. W i Corresponding to the i The current moisture content of each particle; B It is the right-hand term vector, determined by the initial water content of the particles and the boundary / source terms; Formula (13) is the final set of water diffusion equations. The water content of the particles can be updated by solving formula (13). Iterate through all particles and update their water content with the new value.

[0036] Example 4 Furthermore, the particle size can be determined based on the change in particle moisture content; the change in moisture content will cause linear strain, affecting the particle size. i The change in water content is set as granules i The initial radius is Then the granules after drying shrinkage i radiusR i As shown in formula (14); (14) in, The shrinkage coefficient is obtained experimentally. Moisture diffusion can be accounted for using equation (14). Example 5 Furthermore, continuous numerical methods yield little success in addressing the coupling between moisture diffusion and mechanical behavior. In contrast, discontinuous methods such as the discrete element method, discontinuous deformation analysis, and particle dynamics demonstrate superior descriptive capabilities for such problems.

[0037] Example 6 The present invention provides a particle-based discontinuous deformation analysis (DDA)-mechanical coupling calculation method for water diffusion. Within this particle-based DDA analysis framework, it employs an implicit solution strategy, resulting in a large time step and computational stability. Furthermore, it fully considers the influence of cracks on the water diffusion process. Thanks to these improvements, the constructed water diffusion model not only enhances computational efficiency and simplifies the material parameter assignment process, but also more accurately describes the water migration behavior within cracked structures.

Claims

1. A method for calculating water diffusion-mechanical coupling based on particle DDA, characterized in that, The specific steps are as follows: Step 1: Discretize the region to be solved using particle DDA and assign control parameters; Step 2: Conduct contact retrieval between particles and between particles and boundaries to determine the domain of action of mechanical and moisture diffusion, and establish force chains and virtual moisture diffusion channels between them. Step 3: Perform mechanical part calculations for particle DDA and update particle position information; Step 4: Assemble the moisture diffusion equations for particle DDA, solve the equations and update the particle moisture content information.

2. The water diffusion-mechanical coupling calculation method based on particle DDA as described in claim 1, characterized in that, In step 1, the control parameters include geometric parameters, mechanical parameters, and moisture diffusion parameters; the geometric parameters include particle position coordinates and particle size; the mechanical parameters include particle density. Normal contact stiffness Tangential contact stiffness ,tensile strength Bond strength internal friction angle With shrinkage coefficient Moisture diffusion parameters include the moisture diffusion coefficient. With boundary flow mass density .

3. The water diffusion-mechanical coupling calculation method based on particle DDA as described in claim 2, characterized in that, Step 2 specifically involves: Step 2.1: Construct a unified particle contact network; Step 2.2: Traverse each pair in the particle contact network and determine whether the center distance of the pair is less than or equal to the sum of the particle diameters of the two particles; if the mechanical contact criterion is satisfied, record the pair and its contact geometry information to form a force chain list; Step 2.3: Traverse each pair in the particle contact network; Determine whether the center distance of the pair is less than or equal to the preset hydraulic search radius; if the hydraulic diffusion criterion is met, record the pair and its channel geometry information to form a list of hydraulic channels; for particle units adjacent to the boundary of the simulation domain, apply the mechanical contact criterion and the hydraulic diffusion criterion respectively to establish force chains or hydraulic channels between the particles and the boundary, and add them to the corresponding lists respectively.

4. The water diffusion-mechanical coupling calculation method based on particle DDA as described in claim 3, characterized in that, In step 2.1, the specific steps are as follows: obtain the position information, particle size information, and geometric boundary information of all computational particle units at the current moment; set the cell side length of a background mesh according to the maximum particle size of all computational particle units in the system, ensuring that the side length is greater than the maximum particle size; Construct a spatial indexed grid defined by the cell side length on the simulation domain; traverse all computational particle cells and calculate the grid cell index to which each cell belongs based on its position coordinates; map the unique identifier of each computational particle cell to its corresponding grid cell to form a grid-particle index relationship. Traverse every grid cell in the spatial index grid; for each target particle cell in the current grid cell, perform pairing detection in that cell and all its neighboring cells; for each potential pairing, calculate the center distance between them; determine whether the center distance is less than or equal to a preset neighborhood search radius; if the condition is met, record the pairing relationship to form a particle contact network containing all spatially adjacent particle pairs; this network is the common action domain of mechanics and moisture diffusion.

5. The water diffusion-mechanical coupling calculation method based on particle DDA as described in claim 3, characterized in that, Step 3 specifically involves: Step 3.1, Force initialization and application of volume forces; Iterate through each computational particle in the system, initialize the resultant force and resultant torque acting on it to zero; calculate and apply a preset volume force based on the mass of each computational particle, and add the force to the resultant force of that particle. Step 3.2, Calculation of contact force based on force chain: Iterate through the force chain list generated in Step 2.2; for each force chain data item in the list, calculate the normal force and tangential force, and update the force and moment. Step 3.3, Motion State Update: After completing the calculation of all forces, traverse all calculated particle units; Based on the current net force and net torque acting on it, combined with its own mass and moment of inertia, its linear acceleration and angular acceleration are solved; a time integration algorithm is used to update the velocity, angular velocity, position and angle information of each computational particle; position constraints are applied to computational particles that exceed the boundary of the simulation domain after the position is updated to ensure that they always move within the simulation domain.

6. The water diffusion-mechanical coupling calculation method based on particle DDA as described in claim 5, characterized in that, In step 3.2, the normal force is calculated as follows: the normal contact force between the pairs is calculated based on the preset normal contact stiffness, normal overlap, and normal relative velocity; the tangential force is calculated as follows: the tangential contact force between the pairs is calculated based on the preset tangential contact stiffness, tangential relative displacement increment, and internal friction angle, combined with the Coulomb friction criterion. Force and torque update: The calculated normal force and tangential force are combined into a total contact force, and the total contact force and the torque it generates are applied to the resultant force and resultant torque of the calculated particle elements at both ends of the force chain.

7. The water diffusion-mechanical coupling calculation method based on particle DDA as described in claim 5, characterized in that, In step 4, the specific methods of moisture diffusion between particles and between particle boundaries are as follows: Let the particles The radius is Let the particles The radius is The transmission surface of the virtual channel for moisture diffusion is a rectangle. The width of the rectangle is defined as the side length being the arithmetic mean of the diameters of the two particles, that is, the sum of the radii of the two particles. Define the water diffusion coefficient of the virtual water diffusion channel. It is the harmonic mean of the water diffusion coefficients of the two particles, as shown in formula (3); （3） in, and Particles With particles The water diffusion coefficient; Assume that water is conducted through a virtual uniform rectangular channel between two particles or between a particle and a wall; as shown in formula (5); （5）； in, and Particles With particles Moisture content; A correction coefficient is introduced to correct the width of the virtual water diffusion channel, as shown in formula (6); （6） Among them, coordination number It is a measure of the average number of neighboring particles for each particle. This represents the volume fraction of a particle aggregate. From the boundary per unit time Flow towards particles Traffic As shown in formula (7); （7）； In the formula, length It is the distance from the center of the particle to the wall, the width of the virtual moisture diffusion channel. and correction factor As shown in formulas (8) and (9); （8） （9） Let the particles have Adjacent particles and Adjacent boundaries, inflowing particles i The total flow includes contributions from adjacent particles and walls, as shown in Equation (10); （10） By substituting formula (10) into formula (5), we obtain formula (11); （11） in, Represents the particles of the previous time step i moisture content, It is a wall w Temperature, Δ t It is the time step; ; Equation (11) can be further written as equation (12); （12） Assume the particulate system consists of n Composed of individual particles, equation (12) can be written in matrix form, as shown in formula (13); AW = B （13） in, ; ； ； ； A It is a coefficient matrix that describes the interaction relationships of water diffusion between particles; a ii : reflects the first i The moisture retention capacity of each particle itself; a ij : Reflecting the first j The particle pairs the first i The water diffusion and transfer effect of individual particles; W It's the water content, for each element. W i Corresponding to the i The current moisture content of each particle; B It is the right-hand term vector, determined by the initial water content of the particles and the boundary / source term conditions; Formula (13) is the final set of water diffusion equations. The water content of the particles can be updated by solving formula (13). Iterate through all particles and update their water content with the new value.

8. The water diffusion-mechanical coupling calculation method based on particle DDA as described in claim 7, characterized in that, The particle size can be determined based on the change in particle moisture content. i The change in water content is set as granules i The initial radius is Then the granules after drying shrinkage i radius R i As shown in formula (14); （14） in, This is the shrinkage coefficient.

Images