A numerical simulation method of reinforced soil based on morphological evolution and boundary response

Abstract

The present application relates to the field of geotechnical engineering numerical simulation, and discloses a reinforced soil numerical simulation method based on morphology evolution and boundary response, to solve the problems of existing model particle morphology distortion, broken mechanism simplification and large deformation calculation error caused by rigid boundary. The method obtains real granular morphology, establishes the directional mapping relationship between the parent critical broken stress and the real fragment template; based on the characteristic particle size of the particles, a discrete element reinforced model is constructed, and a self-adaptive flexible boundary wrapped around the sample is constructed by using discrete element-finite difference coupling; under the cyclic loading of servo confining pressure, the spatial vector mixed product method is used to calculate the closed topological body deformation in real time and update the node external force; at the same time, the octahedral shear stress suffered by the particles is monitored, and after reaching the standard, the fidelity fragment template is replaced in situ according to the mapping relationship and the skeleton is reconstructed. The present application realizes high-fidelity closed-loop simulation of the mechanical behavior of reinforced granular materials from macro and micro dual dimensions.

Description

Technical Field

[0001] This invention relates to the interdisciplinary field of geotechnical engineering and computer numerical simulation technology, and more specifically, to a numerical simulation method for reinforced soil based on morphological evolution and boundary response. Background Technology

[0002] With the continuous advancement of transportation infrastructure and marine engineering construction, various projects often face the problem of foundation treatment for special soil types (such as calcareous sand). These special soil types typically exhibit highly irregular shapes, high porosity, and extreme fragility, among other physical and mechanical characteristics. To improve the bearing capacity of these special foundations and control uneven settlement, reinforcement techniques such as geogrids are widely used for their reinforcement. Under long-term cyclic dynamic loads such as traffic loads or wave loads, reinforced soils undergo complex macroscopic and microscopic structural evolution. The Discrete Element Method (DEM), due to its ability to reveal the discontinuous, large-deformation, and multiphase contact characteristics of granular materials at the particle scale, has become an important numerical tool for studying the macroscopic and microscopic mechanical behavior of reinforced granular systems.

[0003] However, existing discrete element method (DEM) numerical simulations have significant limitations when dealing with reinforced systems containing irregular and fragile particles. Firstly, in representing particle morphology, most existing numerical models simplify the model by using ideal spherical or simple quasi-spherical particle systems. Real, specialized soil particles often possess strong angularity and complex geometric topology; these morphological characteristics are the physical basis for the structural anisotropy of soil and the interlocking effect between particles. Using simplified spherical particles not only eliminates the real mechanical interlocking force between particles but also fails to reflect the interlocking resistance between particles and geogrid mesh when simulating reinforced interfaces. To compensate for the insufficient macroscopic strength of the model, studies often artificially amplify the contact friction coefficient between particles, leading to significant deviations of microscopic mechanical parameters from their true physical meaning.

[0004] Secondly, in simulating particle breakage behavior, particle breakage under high stress or cyclic shear is the dominant factor causing volume compression and strength evolution in special soil masses. Existing particle breakage simulations typically rely on simple fragment replacement methods or bonded particle models, whose breakage products are mostly aggregates of spherical particles with scaled-down dimensions. This simplified breakage mechanism ignores the initial influence of the parent particle shape on breakage strength, and the generated fragments have a uniform morphology, completely failing to reproduce the flaky or strip-like debris characteristics formed after real particle breakage. The loss of morphological information makes it impossible for existing models to accurately track the continuous weakening mechanism of the overall soil strength caused by particle size distribution evolution, porosity redistribution, and local skeleton instability due to breakage.

[0005] Furthermore, in terms of boundary condition application and macroscopic response monitoring, traditional numerical triaxial tests primarily rely on rigid cylindrical sidewalls to apply confining pressure constraints. In actual indoor triaxial tests, the specimen is encased in a flexible rubber membrane. As shear failure progresses, the specimen's side surface exhibits significant local bulging or non-uniform uneven deformation. The rigid cylindrical boundary imposes excessive lateral constraints, completely limiting this non-linear lateral deformation, causing the development of internal stress paths and shear bands to deviate from reality. More seriously, once the specimen undergoes the aforementioned non-uniform lateral deformation, the traditional total volume calculation formula based on the ideal cylinder assumption becomes completely invalid. Errors in volumetric strain calculation not only severely interfere with the accuracy of the servo control system in maintaining constant confining pressure but also prevent researchers from accurately obtaining the true shear contraction or dilatation evolution of the soil.

[0006] In summary, existing numerical simulation techniques suffer from systematic distortions in reconstructing realistic particle morphology, microphysical fragmentation evolution, and flexible boundary constraints, making it difficult to accurately characterize the progressive failure process of special particle systems. Breakthroughs in numerical experimental methods are urgently needed. Summary of the Invention

[0007] This invention discloses a numerical simulation method for reinforced soil based on morphological evolution and boundary response, the method comprising: Acquire point cloud data of granular surface and reconstruct it into a three-dimensional solid model, classify it into parent particles of various typical morphologies and extract representative feature particle size; Obtain the single-particle crushing strength parameters and fragment morphology, and establish the directional mapping relationship between the critical crushing stress threshold of the parent particle and the rigid block template of the actual fragment; Based on representative characteristic particle size, a bidirectional geogrid model is constructed, and a flexible boundary is constructed to enclose the particle aggregate and the outside of the geogrid model. The total volume change of the flexible boundary envelope is calculated and fed back to the confining pressure servo control module; Under servo confining pressure, cyclic loads are applied, and the octahedral shear stress on each parent particle is monitored. If the critical fracture stress threshold is reached, the corresponding real fragment rigid block template is generated in the original spatial position according to the directional mapping relationship, and iteratively updated.

[0008] The particles were classified into several typical forms of parent particles, and representative characteristic particle sizes were extracted. Specifically, this included: obtaining the minimum bounding cuboid of the reconstructed 3D solid model and extracting the lengths of the major axes of the orthogonal parts of the cuboid. Central axis length and minor axis length Representative characteristic particle size is calculated based on the axis lengths of three orthogonal directions. The formula for calculating the representative characteristic particle size is: The volume calculated in reverse is equal to the actual volume of the particle. Find the equivalent perfect sphere, and calculate the surface area of ​​the equivalent perfect sphere as its surface area. According to the formula Calculate the sphericity of the particle. ,in This represents the actual surface area of ​​the particle.

[0009] The establishment of the directional mapping relationship between the critical fragmentation stress threshold of the parent particles and the rigid block template of the actual fragments includes the following conversion criteria: a conversion reference size is introduced based on the geometric classification of the parent particles. For blocky and spindle-shaped particles, the conversion reference size is as follows: Take representative characteristic particle size For strip-shaped and sheet-shaped particles, take the minor axis length of the smallest circumscribed cuboid; the peak failure load in the single-particle crushing test is... Transformed into the ultimate tensile stress of particles at the moment of failure : Extracting characteristic intensities to characterize the macroscopic intensity of particles of various shapes and Weber modulus This was used to determine the critical fracture stress threshold for parent particles of various morphologies. .

[0010] Based on representative characteristic particle size, a bidirectional geogrid model is constructed. Specifically, the maximum value among the representative characteristic particle sizes of various typical particle morphologies is extracted as the effective characteristic particle size. Set the net mesh width of the bidirectional geogrid model. With effective feature particle size It satisfies the following linear matching model: in, The mesh matching coefficient is set to a value range of 1. .

[0011] When constructing the bidirectional geogrid model, a parallel bonding model is introduced. The parallel bonding model is equivalent to generating an elastically bonded cylinder at the contact surface between two adjacent particles, and the equivalent maximum normal tensile stress acting at the critical point of the cylinder section is calculated in real time. and equivalent maximum shear stress The calculation model is as follows: in, For bonding normal force, For the bonding tangential force, For bond bending moment, For bonding torque; The cross-sectional bonding radius; This represents the cross-sectional area of ​​the bonding section. The moment of inertia of the bonded section; The polar moment of inertia of the bonded section.

[0012] Constructing a flexible boundary around the particle assembly and the grid model includes an equivalent mapping stage from contact forces to nodal forces: particles within the discrete element computation system collide with the surface elements corresponding to the shell elements in the finite difference system. For the shell elements... The corresponding wall surface element exists The first particle contact point will be the first particle contact point. The contact force vector exerted by the particles on the boundary at each contact point is denoted as . For shell unit nodes The equivalent contact node force vector it receives The calculation formula is: in, For the first One contact point in the housing unit Area coordinate weights within the area.

[0013] Performing dynamic updates of external load application and boundary deformation displacement in a finite difference system: servo confining pressure Converted to apply to each node Equivalent confining pressure nodal force The formula for calculating the equivalent confining pressure nodal force is: in, For nodes The subordinate load area, whose value is shared by this node. Sum the current surface areas of all shell elements by adding up one-third of each element. For nodes The current unit outward normal vector of the surface; the equivalent confining pressure nodal forces With equivalent contact node force The external load vector of the node is obtained by summing the vectors, and the unbalanced force is obtained by subtracting the internal force of the shell element node.

[0014] The specific process for calculating the total volume change of the flexible boundary envelope is as follows: the geometric surface enclosing the particle assembly is considered as a closed set of discrete surface element meshes. The set contains Triangular face elements All face elements The surface elements are all oriented with their normals pointing outwards from the sample. The instantaneous spatial coordinate vectors of the three vertices are respectively , and ; Set the global coordinate origin and the surface element The three vertices are connected to form a tiny triangular pyramid. The directed volume contribution of this tiny triangular pyramid is... Determined by the product of the three position vectors: Traversing the grid set All of them The total real-time volume after irregular deformation is obtained by algebraically summing the individual facets. The calculation formula is: .

[0015] Feedback is sent to the confining pressure servo control module, including: Calculate the real-time volumetric strain that reflects the macroscopic shear response of the soil. : in, This represents the total volume of the initial closed mesh. The servo control module retrieves the preset constant target confining pressure. Based on the reconstructed topology, dynamic calculations are performed and each node is updated. The equivalent confining pressure nodal force required at the location The calculation formula is: in, To update the node's topology value according to the current time step Subordinate load-bearing area, For the nodes after updating the topology Unit outward normal vector.

[0016] Monitor the octahedral shear stress on each parent particle. If the critical fracture stress threshold is reached, replace and generate the corresponding real fragment rigid block template in the original spatial position according to the directional mapping relationship. Specifically, iterate through any block of parent particles. Extract the maximum value of the contact force modulus at all contact points. Based on the current conversion reference size of this particle Calculate the octahedral shear stress inside the particle. : If subjected to octahedral shear stress Reaching the critical fracture stress threshold corresponding to the parent body shape category Then, retrieve the corresponding subset of real fragment rigid block templates; let the sum of the original volumes of each fragment template within the retrieved template subset be... Apply a uniform linear scaling factor to all fragments within this subset. This makes the total volume of the scaled fragments equal to the parent particle. actual volume linear scaling factor The formula for calculation is: Deleting native particles in the numerical space The scaled fragment template set is filled in situ with the original centroid coordinates as the reference point to generate a new pure friction contact pair.

[0017] The beneficial effects of this invention are as follows: By reconstructing realistic irregular particles and coupling them with a discrete element reinforced network based on particle size matching, this invention realistically reproduces the mechanical interlocking and interface interlocking effects caused by the strong angularity of special soils (such as calcareous sand) in numerical space. Based on this, this invention abandons the traditional spherical equivalent fracture mechanism and establishes a directional mapping relationship between the parent particle morphology and a real physical fragment template library. When the high stress concentration caused by interlocking triggers particle failure, in-situ fragment replacement with accurate morphology is performed. This not only strictly ensures the conservation of mass and volume before and after fracture but also transfers the real geometric and topological characteristics of the fragments to subsequent calculations, effectively overcoming the problem of skeletal resistance distortion caused by the loss of fracture morphology in traditional numerical models.

[0018] To address the shortcomings of traditional rigid boundaries in reflecting non-uniform deformation of specimens, this invention constructs an adaptive flexible boundary. This boundary can sense in real time the attenuation of local support forces caused by skeletal rearrangement and sudden particle breakage within the specimen, adaptively generating highly irregular concave-convex surface deformation, truly equivalent to the flexibility constraints in physical experiments. Combined with a spatial vector hybrid product algorithm, this invention can ignore the degree of deformation irregularity and accurately integrate the closed surface topology to obtain the true real-time volume. This not only eliminates the calculation error of volumetric strain indices under large deformations but also provides accurate geometric input for the servo control module to dynamically correct the external forces at boundary nodes.

[0019] In summary, this invention constructs an analytical method that deeply coordinates morphological evolution and boundary response. The microscopic fragmentation of real particles and the force chain reconstruction of the reinforced interface directly drive localized non-uniform deformation of the flexible boundary; while the dynamic update of the external boundary topology, through the recalculation of the actual volume and loaded area, feeds back to the servo system, thereby accurately correcting the external stress conditions acting on the soil. This two-way feedback mechanism ensures a high degree of mechanical consistency in the progressive failure evolution process of reinforced granular materials under complex dynamic loads, providing a high-fidelity numerical simulation method for the stability assessment of special soil foundations in marine and transportation engineering. Attached Figure Description

[0020] Figure 1 This is a flowchart of the numerical simulation of reinforced soil based on morphological evolution and boundary response according to the present invention; Figure 2 These are fragment images of the parent particle before and after it broke apart; Figure 3 The diagram shows the calculation results of local bulging / depression of the flexible boundary, where (a) is the initial state (10 cycles), (b) is the calculation result after 10 cycles, and (c) is the result after 20 cycles. Figure 4 The following diagrams show the evolution results of the particle-grid-flexible boundary: (a) shows the initial interlocking state, (b) shows the particle breakage and fragment re-construction results, and (c) shows the coupling evolution results. Detailed Implementation

[0021] Example 1 This embodiment provides a numerical simulation method for reinforced soil based on morphological evolution and boundary response. The application scenario of this application is for special soils with high porosity, irregular morphology, and fragile characteristics, preferably the analysis of the stress-deformation and fragmentation characteristics of calcareous sand granules under reinforced conditions in marine engineering. The specific experimental equipment layout and initial parameter settings of the numerical model are as follows: During the data acquisition phase, the system involves the arrangement of two core experimental setups. The first setup is a three-dimensional topography acquisition device, employing a high-precision three-dimensional laser scanner to non-destructively capture surface point cloud data of actual soil particles and output a three-dimensional solid model file. The second setup is a single-particle crushing test device, which consists of a top pressure plate, an axial force-displacement sensor, and a bottom pressure plate. In the specific arrangement and operation, the bottom pressure plate is fixedly connected to the base to maintain a static and absolutely horizontal position. The top pressure plate is connected to the drive mechanism and the axial force-displacement sensor to apply a constant vertical load to the particle at a rate preferably one millimeter per minute until the tested granular particle fractures, thereby obtaining the stress intensity parameters at the time of single-particle crushing.

[0022] In selecting the target granular sample and setting classification parameters, this embodiment uses calcareous sand as the target soil. Through morphological observation and clustering of actual samples, the basic particle samples are categorized into four typical three-dimensional geometric shapes: spindle-shaped, blocky, strip-shaped, and sheet-like particles. After acquiring the surface point cloud data of each particle using a high-precision three-dimensional laser scanner and reconstructing it into a three-dimensional solid model, the model file is imported into a data processing terminal for three-dimensional geometric parameter calculation. The size reference of each particle shape is characterized by the characteristic particle diameter. The calculation logic of the characteristic particle diameter is to obtain the lengths of the principal axes in three orthogonal directions of the particle and calculate the geometric average of these three principal axis lengths. The degree of irregularity of the shape is characterized by sphericity, which is defined as the ratio of the surface area of ​​an equivalent sphere with the same volume as the measured particle to the actual surface area of ​​the particle. The sphericity parameter is a dimensionless parameter, with a value greater than zero and less than or equal to one. Specifically, in this embodiment, the characteristic particle sizes of the spindle-shaped, blocky, strip-shaped, and sheet-shaped particles are defined as 4.2 mm, 4.7 mm, 5.1 mm, and 3.3 mm, respectively, with corresponding sphericity values ​​of 0.78, 0.75, 0.66, and 0.45. Furthermore, at least twenty-five independent particles of each shape were obtained for volume and surface area statistical analysis. This confirmed that the volume and characteristic particle size satisfy a specific first power-law relationship, and the surface area and characteristic particle size satisfy a specific second power-law relationship. Distribution parameters were obtained by fitting the measured data using the least squares method. Finally, the median volume and median surface area of ​​each shape were extracted as representative geometric parameters for each particle shape, used for space filling constraints in subsequent numerical models.

[0023] Regarding the construction and setup of the realistic fragment rigid block template library, crushing tests were conducted on twenty-five particles of each typical shape using the aforementioned single-particle crushing test apparatus. By collecting the critical destructive force of each test, the strength data were statistically fitted using the Weber distribution function to extract the characteristic strengths and Weber moduli of spindle-shaped, blocky, strip-shaped, and sheet-like particles. Furthermore, a high-precision 3D laser scanner was used to scan the morphology of each fragment particle after crushing. Based on the scanned images, rigid block particles of the fragments were constructed in a numerical model, forming a fragment rigid block template library for later use. Under these constraints, the sphericity distribution of the fragment rigid block particles ranged from 0.51 to 0.78.

[0024] The apparatus and structure of the numerical triaxial test basic model were used to establish a virtual triaxial test environment within computer discrete element simulation software. The numerically reinforced specimen was constructed as a cylindrical structure filled with rigid block particles, reinforced in the middle, and surrounded by a flexible boundary. The initial design diameter of the cylindrical structure was limited to 100 mm, and the height was limited to 200 mm. The total number of particles generated inside was limited to 6815, with the particle size controlled within the range of 2 mm to 10 mm, the average particle size being 4.5 mm, and the initial porosity of the specimen was kept stable at 0.52. The basic micromechanical parameters of the numerical specimen were set as follows: a linear contact model was used between particles, the particle density was set to 2775 kg / m³, the coefficient of friction between particles was set to 0.55, and the critical damping ratio was set to 0.8.

[0025] Regarding the structural and microscopic assembly constraints of the reinforcement device, a bidirectional deformable geogrid model was installed inside a cylindrical numerical specimen. The geogrid model was set to 60 mm in both length and width, containing multiple matrix-distributed mesh openings, with each mesh opening measuring 30 mm by 30 mm. The reinforcement was not a single continuous medium, but rather composed of two types of discrete element clusters bonded together using a parallel bonding model. The reinforcing nodes of the geogrid consisted of a large particle with a diameter of 4 mm and twenty smaller particles with a diameter of 2 mm evenly distributed around it. The rib structure of the geogrid used three overlapping particle clusters with a diameter of 2 mm as basic components, and the ribs at each mesh edge contained fifteen of these basic particle clusters. The microscopic bonding parameters were pre-tested using a virtual test on the numerical geogrid model, with one end fixed and the other end subjected to a strain rate of 1 / 100th of a minute, and then compared and calibrated with the stress-strain curve obtained from laboratory physical tensile tests.

[0026] In the parameter settings for numerical loading boundary and external excitation conditions, upper and lower rigid loading walls with rectangular cross-sections are respectively arranged at the top and bottom of the reference profile of the cylindrical numerical specimen. The upper and lower rigid walls are configured with contact mechanical property parameters matching those of the rigid block particles, and are limited to vertical displacement along the axial direction during numerical calculations to transmit vertical loading force, without generating lateral friction or displacement constraint effects. After completing the initial consolidation calculations, a constant radial confining pressure of 100 kPa is applied to the specimen through a lateral flexible boundary. During the loading test, axial cyclic loads are applied to the upper and lower rigid rectangular walls, and the applied loading waveform is limited to a half-sine wave to simulate traffic loads in actual engineering. The loading frequency is set to 2 Hz, and the total number of effective cyclic loading cycles is set to twenty.

[0027] Based on the hardware device foundation and numerical parameter limitations set in Embodiment 1 above, Embodiment 2 mainly provides a general description of the overall technical solution of this application. This embodiment provides a numerical simulation method for reinforced soil based on morphological evolution and boundary response, and its overall execution process includes the following steps: S1. Obtain the surface point cloud data of the target granules and reconstruct it into a three-dimensional solid model. Classify the three-dimensional solid model into parent particles of various typical morphologies and extract the representative features of particle size and sphericity of each typical parent particle morphology.

[0028] The purpose of obtaining the true three-dimensional geometric topology of granules is to replace the traditional ideal sphericity assumption in numerical simulations. By extracting geometric feature parameters, especially the representative feature particle size that characterizes size and the sphericity that characterizes the degree of irregularity, we not only provide basic data support for generating equivalent scaled rigid blocks in the numerical computation space, but also establish prerequisites for determining the size matching relationship of the stiffening network apertures.

[0029] S2, obtain the crushing strength parameters of various typical single particles under monotonic loading and the actual fragment morphology after crushing, so as to establish the directional mapping relationship between the critical crushing stress threshold of the parent particles of each morphology and the rigid block template of the actual fragments.

[0030] The shear resistance and fragmentation resistance of granular particles are significantly correlated with their initial geometry. Single-particle crushing tests were conducted beforehand, and fragments were scanned and reconstructed to digitally archive the actual crushing failure patterns and debris morphology evolution. This mapping and correlation mechanism allows for direct replacement of a matching real fragment morphology in subsequent numerical simulations of complex loads once the stress on a parent particle reaches its strength limit. This ensures the continuity of mass and equivalent volume transfer in the reinforced soil during the fracture evolution process, avoiding the severe morphological distortion caused by traditional proportional volume replacement methods.

[0031] S3. Based on the size constraint of the representative characteristic particle size, a bidirectional geogrid model containing mesh, reinforcing nodes and ribs is constructed using discrete element cluster units.

[0032] The mechanical enhancement response of reinforced structures in soil media is highly dependent on the matching state between their mesh size and the surrounding particle morphology. By using the actual particle characteristic particle size obtained in step S1 as boundary conditions to constrain the pore size and structural node size of the reinforced model, the mechanical interlocking effect and interfacial friction resistance of the irregular granular community on the reinforced mesh structure are realistically reproduced from the geometric contact level in three-dimensional space.

[0033] S4 couples the rigid wall elements in the discrete element calculation system with the shell elements in the finite difference calculation system in terms of displacement and mechanical parameters to construct a flexible boundary that wraps around the particle aggregate and the bidirectional geogrid model.

[0034] Traditional rigid cylindrical wall calculations forcibly constrain the non-uniform lateral deformation of the specimen, leading to the development of shear bands within the computational system that deviates from reality. By introducing shell elements without lateral stiffness failure limitations for coupled constraint, the constructed external wall can adaptively generate surface stretching and local concave-convex deformation in real time, following the slippage, rearrangement, and extrusion displacement of the internal particle community. This allows for accurate characterization of the flexibility constraint response behavior of the film boundary under shear conditions in actual indoor experiments.

[0035] S5, the spatial volume enclosed by the flexible boundary is accumulated and calculated in real time using the spatial vector hybrid product method to obtain the total volume change of the sample, and the total volume change is fed back to the confining pressure servo control module.

[0036] Because flexible boundaries can expand or contract locally under internal stress, conventional volumetric algorithms based on standard cylinders will produce significant errors. A method that constructs a spatial vector model and performs cross-product operations can accurately capture any nonlinear deformation of the boundary surface. This calculation process not only objectively reflects the volumetric compression or dilatation evolution of reinforced soil under stress, but also provides accurate volumetric strain feedback parameters for the servo control system, ensuring that numerical samples can accurately maintain the preset lateral constant stress environment even under complex large deformation conditions.

[0037] S6. Apply cyclic load under servo confining pressure constraint and monitor the octahedral shear stress state of each parent particle in the particle assembly in real time. If the octahedral shear stress reaches the critical fracture stress threshold, the corresponding real fragment rigid block template is generated in the original spatial position according to the mapping relationship, and the stress and deformation field of the overall skeleton are iteratively updated.

[0038] Under cyclic stress, particles reaching their strength limit are immediately identified and removed, and replaced in situ with a group of real fragments possessing corresponding physical properties. This micro-fracture and replacement process directly triggers local pore structure collapse and skeletal rearrangement. This redistribution of local stress is further transmitted to the surrounding bidirectional geogrid through interface interlocking and diffuses outward, causing changes in the morphology of the flexible boundary, ultimately triggering a recalculation of the total volume and dynamic correction of the servo boundary pressure. Thus, this scheme constructs a closed-loop numerical analysis of the entire process from micro-particle fracture affecting macro-boundary deformation to boundary response reacting on internal stress evolution.

[0039] Example 3 This embodiment details the execution steps of the initial state construction and feature extraction stage of parent particles in the numerical simulation of reinforced soil. Step 1.1: Use a 3D high-precision laser scanner to acquire the real surface point cloud data of the target special soil, reconstruct the 3D solid model, and classify it into typical shapes according to geometric features.

[0040] When numerically modeling granular materials such as calcareous sand, which have high porosity, highly irregular surfaces, and are rich in primitive biological skeleton features, using standard geometries (such as spheres or ellipsoids) can cause the particle swarm to lose its true initial pore distribution structure during accumulation. Therefore, this embodiment directly acquires the three-dimensional morphology of the particle samples.

[0041] The sieved and cleaned single-grain calcareous sand sample is fixed on the stage of a three-dimensional high-precision laser scanner. The scanner is activated, and the spatial coordinates are calculated by the time difference between laser emission and reflection, obtaining a high-density set of three-dimensional spatial coordinate points covering the surface of the grain, i.e., the real surface point cloud data.

[0042] The collected real surface point cloud data is imported into a 3D reconstruction terminal for meshing. Due to the numerous natural grooves and micro-protrusions on the calcareous sand surface, spatial filtering is performed on the point cloud data to remove environmental noise. Subsequently, a spatial triangulation algorithm is used to connect the discrete point cloud data to form a continuous surface triangular mesh. The mesh model undergoes closure checks and topology repair, filling in mesh holes caused by scanning dead angles, generating a 3D solid model with a defined spatial volume and a closed surface.

[0043] After reconstructing the 3D solid models, morphological classification was performed based on the proportions of the 3D outline dimensions of the particles. The maximum length-width-thickness ratio of the model in space was obtained by calculating the minimum circumscribed cuboid of the reconstructed model. Based on this proportional relationship, all 3D solid models were forcibly classified into four typical morphologies: spindle-shaped, massive, strip-shaped, and sheet-like. These four morphologies encompass the mainstream geometric and topological characteristics of calcareous sands formed during marine sedimentation, and different morphological categories directly determine the initial anisotropic fabric of granular communities during deposition.

[0044] Step 1.2: Extract the three-dimensional geometric parameters of particles of each shape, and calculate the representative characteristic particle size of each shape. and sphericity .

[0045] After completing the reconstruction and classification of the digital 3D solid model, it is necessary to calculate the 3D geometric parameters of the particles to provide a numerical benchmark for subsequent reinforcement network size matching and resistance calculation.

[0046] For each individual 3D solid model, the actual volume of the 3D solid model is calculated by traversing all triangular mesh elements on the model surface using a spatial geometric integral algorithm. and actual surface area Simultaneously, based on the constructed minimum circumscribed cuboid, the side lengths of the three mutually orthogonal directions of the cuboid are extracted and defined as the lengths of the major axes. Central axis length and minor axis length .

[0047] Representative characteristic particle size was calculated based on the above three orthogonal axis lengths. The extracted major axis length Central axis length and minor axis length Perform a series of multiplications, and then take the cube root of the product, i.e., calculate the geometric mean of the lengths of the three orthogonal axes. Use this geometric mean as the representative characteristic particle size. It effectively smooths out the deviation in overall particle size assessment caused by single-dimensional size distortion of extremely flat or elongated particles, thus reflecting the spatial occupancy scale of irregular particles more objectively.

[0048] Calculate the sphericity of this 3D solid model. Due to the numerous surface wrinkles present in marine calcareous sand, its actual surface area is... The material is much larger than a standard sphere of the same volume, a characteristic that is the main reason for the high frictional interlocking force at the reinforced interface. To calculate the extent of this deviation from the ideal sphere, based on the obtained actual volume... Based on the reverse calculation logic of the volume of a standard sphere, a volume equal to the actual volume is calculated in reverse. Find the equivalent perfect sphere, calculate the surface area of ​​the equivalent perfect sphere, and define it as the surface area of ​​the equivalent sphere. The calculated equivalent sphere surface area Divide by the actual surface area of ​​the particle The resulting ratio is the sphericity of the particle. .

[0049] Example 4 For fragile soils such as calcareous sand, traditional numerical simulations typically treat particle breakage as equivalent to the parent material fracturing into several smaller standard spheres. This approach severs the morphological inheritance characteristics of real particles before and after breakage; that is, flaky particles are likely to remain flaky fragments after breakage, causing the virtual skeleton after breakage to lose its realistic mechanical interlocking force. To address this issue, this embodiment constructs an experimentally based real fragment library and mapping mechanism, with the specific steps as follows: Step 2.1: Conduct single-particle crushing tests on real particles of various typical shapes to obtain the characteristic strength and Weber modulus when particles of each shape break, and set the critical crushing stress threshold for each type of particle.

[0050] Due to their marine biogenic origin, calcareous sand particles contain a large number of randomly distributed micropores and natural defects, resulting in significant dispersion and shape dependence in their fracture strength. Therefore, it is necessary to conduct independent mechanical tests for different morphologies.

[0051] From the spindle-shaped, blocky, strip-shaped, and sheet-shaped particle samples classified in Example 3, a sufficient number of independent single-particle samples were randomly selected for each typical morphology. Each single sample was placed on the base of a single-particle crushing test apparatus, and the top pressure plate was driven to apply monotonic axial compression to the particle at a constant displacement rate. Simultaneously, load data was recorded in real time using an axial force-displacement sensor until the particle underwent brittle fracture, and the peak failure load at this point was extracted. A conversion reference size is introduced based on particle geometry classification. For both blocky and spindle-shaped particles, a representative characteristic particle size is taken. For elongated or flattened strip-shaped and sheet-shaped particles, take the length of the minor axis of the smallest circumscribed cuboid, and apply the formula... Peak failure load (Newton, The stress is transformed into the ultimate tensile stress of the particle at the moment of failure. (Pascal, ).

[0052] For all sample test data of each typical morphology, a conventional Weber statistical distribution model from rock mechanics was introduced for processing. By logarithmic fitting of the ultimate tensile stress data and the particle survival probability, characteristic strengths representing the overall macroscopic strength of particles of this shape were extracted. This refers to the stress level at which the survival probability decreases to a specific baseline value, and the Weber modulus, which characterizes the uniformity of defect distribution within particles of this shape. .

[0053] Finally, based on the extracted feature strength and Weber modulus Based on the stress calculation criteria set in the numerical simulation, specific critical breakage stress thresholds were calibrated and assigned to the four types of parent particles: spindle-shaped, blocky, strip-shaped, and sheet-shaped. .

[0054] Step 2.2: Use a three-dimensional high-precision laser scanner to scan each of the real fragments produced after the single-particle crushing test, and reconstruct a template library of real fragment blocks that are consistent with the actual crushing morphology.

[0055] After completing the monotonic loading in step 2.1 and causing the particles to completely fracture, carefully collect all the original fragments remaining on the base of the experimental setup. Since the initial geometry of the particles controls the principal stress transmission path within them, parent particles of different shapes will exhibit drastically different fragment topologies after fracture.

[0056] The collected real fragments were then subjected to high-density surface point cloud scanning one by one using a 3D high-precision laser scanner. Consistent with the reconstruction logic of the parent particle in Example 3, the point cloud data of the fragments was meshed, denoised, and closed to generate a 3D solid model of the fragments.

[0057] All generated 3D fragment models are imported into a database to create a realistic fragment template library. Each fragment template in this library retains the sharp edges, cleavage surfaces, and specific length-width-thickness ratios produced during the actual crushing process.

[0058] Step 2.3: Establish the directional mapping relationship between the parent particle and the fragments, that is, establish an associated database of the parent particle shape category, the specific fragmentation threshold, and the corresponding real fragment template.

[0059] After obtaining the mechanical thresholds of the parent body in various forms and the real fragment templates, it is necessary to establish data retrieval logic in the background of the numerical simulation system to ensure that the fracture evolution in the numerical calculation process conforms to objective laws.

[0060] A relational database is built in the background, and the data call paths are bound one-way to form a directional mapping relationship. Specifically: the first index level of the system is set to the parent shape category (spindle-shaped, block-shaped, strip-shaped, sheet-shaped) as defined in Example 3; the first index level is mapped down to the second level, which is the critical fracture stress threshold specific to this category as marked in step 2.1. The second level is further absolutely mapped down to the third level, which is the specific real fragment template subset corresponding to the physical crushing of this type of specific shape in sub-step 2.2.

[0061] By establishing a relational database with the aforementioned directional mapping relationships, the system possesses an architecture capable of autonomous judgment and execution in numerical multiphysics cyclic calculations. When stress concentration occurs in any rigid block particle within the numerical sample, its original parent body shape category is identified, and the corresponding exclusive fracture threshold is automatically retrieved for resistance determination. Once failure is determined, the system strictly follows the index path of the database to call its corresponding real fragment template for spatial replacement.

[0062] Next, this embodiment further supplements the experimental data and established mapping relationships obtained in the foregoing embodiments. This embodiment discloses a table of statistical parameters obtained from experiments conducted on a specific target soil (calcareous sand), and the dictionary mapping relationships constructed therefrom. These data are directly derived from high-precision scanning and crushing tests on calcareous sand samples from an island reef in the South China Sea.

[0063] Independent samples were randomly selected from the washed and sieved granular material, and surface point clouds were acquired and solid models reconstructed using a 3D high-precision laser scanner. Based on the aspect ratio of the principal axes of the smallest circumscribed cuboid, the samples were strictly divided into four categories: spindle-shaped, block-shaped, strip-shaped, and sheet-shaped. Representative feature particle sizes for each shape were extracted and statistically analyzed using a spatial integration algorithm. and sphericity The specific physical measurements and statistical results are shown in Table 1.

[0064] Table 1: Statistical Table of Initial Geometric Characteristic Parameters of Parent Particles in Special Soils (Calcic Sand) For the samples whose morphology records have been completed in Table 1, they were sequentially fed into a single-particle fragmentation test apparatus and subjected to monotonic axial compression at a constant displacement rate. The peak failure load of each sample was recorded. The stress was then converted to the ultimate tensile stress. The characteristic strengths of each morphological system were calculated and obtained. and Weber modulus The mechanical property data obtained from the experiment are shown in Table 2.

[0065] Table 2: Statistical Table of Physical and Mechanical Characteristic Parameters of Parent Particles Crushing under Monotonic Loading Blocky and spindle-shaped particles have high characteristic strength. And Weber modulus The relatively high characteristic strength indicates a relatively uniform distribution of internal defects; however, plate-like and strip-like particles are highly susceptible to brittle fracture due to localized bending, resulting in low characteristic strength. Significantly low, and the Weber modulus The value is extremely low, reflecting the huge dispersion of fracture.

[0066] The debris and residue from the single-particle crushing experiments were collected and subjected to high-precision scanning and digital reconstruction. A three-level association database containing parent category, threshold conditions, and fragment subsets was established in the numerical control system to achieve precise binding of morphological inheritance throughout the particle crushing process. The specific directional mapping relationships are shown in Table 3.

[0067] Table 3: Database of morphological evolution of parent particles and orientation mapping of rigid block templates of real fragments Example 5 This embodiment details the implementation steps for the reinforcement modeling stage in the numerical simulation of reinforced foundations for special soil types (such as calcareous sand). In reinforced calcareous sand systems, due to the highly irregular geometry of the calcareous sand, its interface interaction with the geogrid is highly dependent on the mechanical interlocking effect of particle edges embedded in the geogrid mesh. Traditional numerical simulations often use rigid, thicknessless structural elements as equivalent geogrids, which cannot reflect the flexible tension and node failure of the geogrid itself, and the smooth boundaries lose the basis for three-dimensional spatial embedding and interlocking with irregular sand particles.

[0068] To ensure that the reinforced body in the numerical space can achieve a deep mechanical engagement with the shaped particles constructed in Example 3, the mesh size of the reinforced body must be matched according to the actual space occupied by the particles.

[0069] Retrieve the representative characteristic particle size data of the target granular particles obtained in Example 3. Since there are four typical morphologies, to ensure that even the coarsest particles can participate in local penetration and interlocking, this step extracts the maximum value among the characteristic particle sizes of the four morphologies as the calibration benchmark, denoted as the effective characteristic particle size. Based on the interlocking principle in foundation treatment engineering, the net mesh width of the bidirectional geogrid is set. The net width of the mesh With effective feature particle size It satisfies the following linear matching model: in, The mesh matching coefficient is set to a value range of 1. . The number of penetrating particles that can be accommodated within a single mesh is insufficient, making it difficult to form a multi-particle cooperative interlocking structure, and the interfacial interlocking degenerates into surface friction as the primary mechanism; if If the mesh is too large, the entire particle group will leak out through the mesh, making it impossible to form an effective local skeletal arch effect.

[0070] A reinforced network entity with finite spatial thickness is assembled using discrete rigid microspheres. For the reinforced node region of the geogrid, a sphere with a radius of... The core main particle, and several smaller particles with radii tightly surrounding it in three-dimensional space. Secondary particles form reinforcing node clusters; for the rib region, a series of partially overlapping particle clusters are generated along the rib extension direction, with the radius of the overlapping particles being... The spatial arrangement of multiple particle clusters gives the assembled geogrid ribs a periodic textured surface, which geometrically replicates the mechanical interlocking.

[0071] To reflect the tensile deformation and tearing failure of local nodes of the geogrid under cyclic loading, a parallel bonding model is introduced between all adjacent micro-cluster particles that make up the geogrid.

[0072] This model is equivalent to generating an elastically bonded cylinder with a specific circular cross-sectional area at the contact surface of adjacent particles. Let the bonding radius of this bonded cylinder be... Its value selection rules are set as follows ,in , The radius of the two adjacent particles that come into contact. As the bonding radius multiplier, it is preferred to set When the reinforced body undergoes relative deformation due to internal sand and soil compression, the bonded model will generate a corresponding bond normal force. Adhesive tangential force Bonding moment and bonding torque .

[0073] Specify the adhesive normal force Tension is taken as positive and compression as negative. The equivalent maximum normal tensile stress acting at the critical point of the bonded cylinder section is calculated in real time. and equivalent maximum shear stress Its calculation model is as follows: in, The cross-sectional area of ​​the bond is calculated using the following formula: The unit is square meters ( ); The moment of inertia of the bonded section is calculated using the following formula: The unit is 4th power meters ( ); The polar moment of inertia of the bonded section is calculated using the following formula: The unit is 4th power meters ( The first term on the right-hand side of the above formula represents the uniform stress on the cross section caused by the normal force and the tangential force, respectively. The final calculated... and The units are all Pascals (Pa) ).

[0074] Set a normal tensile strength threshold for this parallel bond model. and tangential shear strength threshold In each calculation step, if it is detected that... or If the microscopic bonding bond is broken, the corresponding local cluster particles degenerate into pure frictional contact, thus reproducing the failure behavior of the geogrid.

[0075] Based on the aperture matching relationship set in step 3.1, the irregular parent particles constructed in Example 3 will be embedded in the mesh. The internal or external textured surface of the ribs forms a rigid contact constraint. In this state, the parallel bonding bonds defined in step 3.2 have a strong deformation constraint capability, restricting the free slippage of surrounding particles. The combined effect of the load and the strong constraint of the mesh will induce extremely high local contact forces in the grid ribs and their surrounding areas. This contact force is converted into the equivalent principal tensile stress induced inside the particles; due to the stress concentration effect in this area, its equivalent principal tensile stress will quickly reach the specific critical breakage stress threshold set for various parent particles in Example 4. .

[0076] Example 7 In triaxial tests conducted on foundation soil, the actual specimens were wrapped with a flexible rubber membrane. Under cyclic loading, the lateral displacement distribution inside the specimens was uneven due to particle folding, local confinement of the reinforced grid, and skeleton collapse caused by sudden breakage of fragile particles, resulting in complex local bulging or depressions on the side surfaces.

[0077] To enable the numerical boundary to have the geometric degrees of freedom to undergo arbitrary local surface deformation, the cylindrical envelope surface surrounding the particle aggregate and the bidirectional geogrid is discretized.

[0078] The initial cylindrical side surface is divided into several continuously joined triangular planar elements. Each triangular planar element is defined as a shell element in the finite difference computation system. For any single shell element... Its spatial form is determined by the coordinates of three difference nodes. Let these three nodes be nodes 1, 2, 3, and 4 respectively. ,node ,node Their instantaneous spatial position vectors in the global coordinate system are denoted as follows: , , .

[0079] A series of triangular facets are generated synchronously in the discrete element module. Each facet is strictly bound to a unique shell element, and the three vertices of both share the same spatial coordinates at any computation time step (i.e., ...). , , In a discrete element system, block particles only come into contact and collide with the surface element.

[0080] The material constitutive model of all shell elements is set as an isotropic linear elastic model with no strength failure threshold.

[0081] The shell unit is assigned three core material parameters: shell thickness Elastic modulus and Poisson's ratio Preferably, to maintain consistency with the physical properties of the latex film in conventional triaxial testing standards and ensure the feasibility of the numerical model, the parameter setting range in this embodiment is as follows: millimeters ( ), Megapascal (MPa) ), (Unitless). Ensure that the virtual boundary exhibits the same in-plane tensile stiffness response as the real latex film when resisting the compression of internal particles.

[0082] Establish the spatial binding association between DEM surface elements and FDM shell elements, and execute the following bidirectional mapping solution loop in each computation time step: The first stage is the equivalent mapping of contact forces to nodal forces (DEM to FDM). Within any computation time step, block particles collide with rigid wall elements, generating a series of microscopic contact forces. Consider a shell element... The corresponding wall surface element has The particle contact point, the first The contact force vector (including normal and tangential friction components) exerted by the particles on the boundary at each contact point is denoted as... (Unit: Newton) The system, based on the principle of area coordinate interpolation, distributes the discrete contact force to the three vertices of the shell element. For the nodes of this element... ( The equivalent contact node force vector it receives The calculation formula is: in, For the first Each contact point in the unit The area coordinates within (i.e., the three-node linear interpolation weights calculated based on the projection position of the contact point in the element plane, and satisfying...) ).

[0083] The second stage involves the application of external loads and the dynamic updating of boundary deformation displacements (FDM to DEM). This includes servo confining pressure... (Unit: Pascal) Converted into equivalent confining pressure nodal forces acting on each node. .

[0084] For nodes The system iterates through all shell elements sharing the node, and accumulates one-third of the current surface area of ​​each element as the subordinate load-bearing area of ​​the node. and along the unit outward normal vector of the surface at that instant. When pressure is applied in the opposite direction (i.e., inward), the equivalent confining pressure nodal force is calculated as follows: .

[0085] The finite difference module will convert the equivalent confining pressure nodal forces Contact equivalent node force The vectors are summed to form the external load vector of the node. The difference between this external load vector and the nodal internal force corresponding to the current stress state of the shell element is calculated to obtain the actual unbalanced force of the node. Based on the unbalanced force and the preset shell material stiffness matrix, the dynamic equations are solved to calculate the spatial displacement increment vector of each node at the current time step. Simultaneously refresh the coordinates of each shell element vertex and the bound discrete element wall surface vertex in the global coordinate system. ), to complete the reconstruction of the geometric surface.

[0086] The aforementioned two-way coupling mechanism of mechanical perception-force transformation-surface reconstruction directly corresponds to the evolution process of reinforced special soil. During cyclic loading, when a parent particle close to the boundary breaks due to local stress concentration, the particle disintegrates in a very short time, resulting in increased contact force between this local area and the external wall surface element. A sudden drop occurs. At this point, the algorithm will capture the decay of this local contact force in real time. The externally acting servo confining pressure... Driven by this mechanism, the finite difference module calculates the incremental displacement of the nodes pointing towards the center of the specimen, causing the shell boundary corresponding to the location of particle breakage to locally indent inward; while in areas such as the reinforced grid nodes that have not broken and are subjected to strong compression from the internal particles, the shell boundary undergoes tensile bulging outward. This mechanism is directly driven by the redistribution of the microscopic internal stress field, and equivalently reproduces the flexible boundary constraint state of the triaxial test in the numerical computation space.

[0087] Example 8 In the numerical simulation of reinforced granular materials, the strong interlocking of the reinforced network nodes with the particle group and the local rearrangement caused by the shear fracture of special parent particles lead to uneven collapse or expansion of the internal skeleton of the sample. This evolution is reflected in the flexible boundary described in Example 7 as highly irregular curved surface deformation. Once the flexible membrane undergoes local concave-convex deformation, the total volume calculation formula based on the cylindrical shape assumption is no longer applicable, thus causing volumetric strain, i.e., the calculation of the index characterizing the shear dilatation or shear contraction characteristics of the soil is distorted; at the same time, the load-bearing area estimated based on the initial regular shape will also cause the applied equivalent external force to deviate from the set confining pressure.

[0088] To capture the true volume of a specimen under arbitrary irregular deformation, the geometric surface enclosing the specimen is considered as a closed set of discrete surface element meshes. The set consists of two parts: the lateral boundary (i.e., the set of all triangular shell elements that undergo real-time deformation in Example 7) and the upper and lower end faces (i.e., the set of closed virtual triangular face elements attached to the top and bottom plate loading walls).

[0089] To ensure To ensure strict closure throughout the entire process of large deformation under load, the following motion binding constraints are set: The axial degrees of freedom of each node on the top and bottom edge rings of the lateral shell mesh are forcibly bound to the axial displacement of the corresponding rigid loading wall, while the radial and circumferential degrees of freedom remain free, in accordance with the clamping constraint of the membrane at the sample cap and base in the equivalent test; the top sealing surface element set is generated by the fan-shaped triangulation of a central node located on the axis of the upper loading wall and the top edge node of the lateral mesh. The central node and the sealing surface element are refreshed synchronously with the axial displacement of the upper rigid wall, and the follow-up refresh mechanism of the bottom sealing surface element is the same.

[0090] Let this closed mesh set Include There are three triangular facets. The vertex sequence of all facets is arranged according to the right-hand screw rule, ensuring that the normals of the resulting surfaces all point outwards from the sample. For each facet... ( Let the instantaneous spatial coordinate vectors of its three vertices be respectively , and .

[0091] The system will set the global coordinate origin. With face element The three vertices are connected to form a tiny triangular pyramid. The directed volume contribution of this tiny triangular pyramid is... Determined by the product of the three position vectors, the calculation model is as follows: Traversal All of them The directional volume contribution of each surface element is algebraically summed to obtain the real-time total volume of the sample after irregular deformation. : Further calculations were performed to determine the real-time volumetric strain reflecting the soil's shear response. : in, Characterizing volume compression, Characterizes volume expansion; At the moment the consolidation calculation is completed, the initial closed mesh is obtained by summing using the spatial vector hybrid product method, which is exactly the same as in the sub-step.

[0092] The reconstructed closed mesh topology parameters (including the real-time total volume) The current surface area and nodal normal vector of each shell unit are fed back to the confining pressure servo control module. The real-time total volume... With volumetric strain As a monitoring output characterizing the shear contraction / dilatation evolution, it is simultaneously used to verify the volume evolution law of the specimen on the deviatoric stress path under servo feedback conditions; the remaining topological parameters are used to dynamically correct the external forces at the boundary nodes.

[0093] The servo feedback described in this application specifically refers to a follow-up and gradual update mechanism for the equivalent confining pressure nodal forces based on real-time topology reconstruction, so as to ensure that the normal stress of the flexible boundary is always equal to the target confining pressure. .

[0094] In the specific feedback operation, the update node coordinates of the lateral shell unit in Example 7 are extracted. For any node on the lateral boundary... Based on the reconstructed topology, adjacent shell units sharing a connection with this node are re-determined. The calculation rule for the subordinate load-bearing area is consistent with that described in Example 7, with superscript... This represents the updated topology value based on the current time step, where one-third of the total actual surface area of ​​adjacent shell elements is set as the node. subordinate load area .

[0095] At the same time, take the node The unit outward normal vector at a node is obtained by weighting and summing the outward normal vectors of adjacent shell elements according to their element surface areas, and then normalizing the sum. .

[0096] The servo control module retrieves the preset constant target confining pressure. (Unit: Pascal) The equivalent confining pressure nodal force required at the node is dynamically calculated and updated. : Updated It will be returned to the mechanical solution stage of Example 7 as the external load component of the current time step of that node.

[0097] Microparticle breakage and mesh engagement induce local contact force redistribution, subsequently driving non-uniform displacement updates at the flexible boundary. Based on this, the closed topological volume, loaded area, and unit normal vector are re-evaluated using the vector hybrid product method. Then, the servo control module updates the nodal external forces accordingly. This mechanism effectively prevents boundary force distortion caused by large deformations, ensuring the mechanical validity of the nonlinear response simulation of reinforced special soils.

[0098] Example 9 In actual service operations of marine engineering foundations containing fragile granular materials such as calcareous sand, wave cyclic loading leads to stress concentration within the soil, resulting in progressive particle breakage. Particle fracturing directly causes local pore collapse and skeleton rearrangement, subsequently triggering abrupt stress changes in the reinforced grid and deformation of its external boundaries. Traditional numerical simulations often employ simple sphere substitution algorithms, losing morphological inheritance characteristics and exhibiting abrupt changes in mass and equivalent morphology before and after fracturing. To effectively address this issue, this embodiment provides a dynamic fracturing and skeleton reconstruction cyclic method based on a collaborative mechanism. The specific implementation steps are as follows: Step 6.1: Apply cyclic load and calculate the equivalent tensile stress of particles in real time.

[0099] Under the constant servo confining pressure constraint of Example 8, axial cyclic loads were applied to the reinforced specimen by loading walls from both above and below. To evaluate the stress state of irregular rigid block particles under complex contact conditions, all parent particles inside the specimen were traversed in real time at each calculation time step, and the octahedral shear stress and mean stress tensor of the particles were calculated.

[0100] For any block of parent particles Let the actual volume of the space it currently occupies be... The surface of the particle contains The first contact point, the... At the point of contact, external forces act on the particles The contact force vector component is denoted as From the center of mass of the particle to the first The branch vector of each contact point is denoted as ( (These are spatial coordinate components). Calculate the average stress tensor of this particle. : The calculated unit average stress tensor The unit is Pascal ( Since the tensor has an antisymmetric component during the dynamic process, the system takes the symmetric part of the tensor. Used for subsequent stress state monitoring.

[0101] Simultaneously, to determine particle breakage, it is necessary to extract characteristic stresses that originate from the same mechanical mechanism as the single-particle crushing test in Example 4. This involves traversing the particle... Extract the maximum value of the contact force modulus at all contact points. (Newton, ), combined with the current conversion reference size of the particle (According to the size reference corresponding to the parent body shape category in Example 4), the equivalent principal tensile stress induced inside the particle is calculated using the same conversion formula as in the calibration process of Example 4. : Step 6.2: Comparison of critical fracture stress thresholds based on morphological dictionary.

[0102] The system extracts the current particles The initial parent body shape categories (i.e., spindle-shaped, block-shaped, strip-shaped, or sheet-shaped) were classified in Example 3.

[0103] Access the established directional mapping association database, use the parent shape category as the first-level index, and retrieve the specific critical fracture stress threshold bound to that shape. .

[0104] Within the current time step, the system compares the equivalent principal tensile stress in real time. With threshold .like If the particle is determined to be in a stable bearing state and the skeleton remains intact, the particle is retained and the calculation proceeds to the next step; if If the microscopic defects of the particle are unstable and expand, a breakage trigger signal is immediately issued, and the process proceeds to step 6.3.

[0105] Step 6.3: Replace the equivalent fragment in situ.

[0106] When receiving particles Upon receiving the breakage trigger signal, morphological evolution is initiated. The particle to be broken is then extracted. Current centroid space coordinates 3D spatial attitude, translational velocity vector and rotational angular velocity vector .

[0107] Based on the third-level mapping instructions of the relational database, a subset of real fragment rigid block templates corresponding to the shape of the parent body is retrieved. Let the sum of the original volumes of all fragment templates within the retrieved template subset be... Apply a uniform linear scaling factor to all fragments within this subset. This makes the total volume of the scaled fragments equal to... In order to maintain the conservation of mass and volume.

[0108] Deleting native particles in the numerical space With the original centroid coordinates Using the original orientation as the reference point, the scaled fragment template set is filled into the spatial region in situ. The slight geometric overlap that occurs during the refilling instant will be eliminated through contact force relaxation in subsequent computational steps. The original translational velocity vector... and rotational angular velocity vector Each newly generated fragment particle is assigned a contact pair, and contact detection is performed again between fragments and between fragments and surrounding unbroken particles to generate new pure frictional contact pairs. In this embodiment, the replaced fragment particles will participate in subsequent contact calculations as stable skeletal units and will no longer be included in the breakage determination, thus simplifying the simulation of a breakage evolution. This step replaces the traditional sphere equivalence method, transferring the genetic morphology of real calcareous sand debris to the next calculation.

[0109] Step 6.4: Skeleton reconstruction.

[0110] After particle breakage, the irregular fragments generated by the aforementioned in-situ replacement are smaller in size and have altered geometric edges compared to the original parent particles. This microscopic morphological evolution will trigger a series of multiphysics coupling responses outward: Internal force chain reconstruction and reinforced interface response: Fragments fall and rearrange into the surrounding pores, and the loss of local support causes the surrounding soil to be squeezed towards the fracture zone. If the fracture zone is located near the geogrid described in Example 5, the soil displacement will pass through the particles in the geogrid mesh, and the microscopic parallel bond model will bear additional tangential shear and tensile forces; if it exceeds its bond threshold, the reinforced network will experience local fracture, reflecting the progressive failure of the reinforced nodes induced by particle breakage.

[0111] Flexible boundary response: As the local collapse caused by rupture is transmitted to the edge of the specimen, the contact force vector borne by the surface element defined in Example 7. A sudden drop occurs. The equivalent contact node force of the outer shell unit decreases accordingly. Under the continuous action of the external confining pressure, the flexible boundary will produce a real inward concave deformation in the lateral region where particle breakage occurs.

[0112] Macroscopic volumetric deformation update and servo feedback: After the flexible boundary undergoes concave deformation, the spatial vector hybrid product algorithm in Example 8 senses the topological change in real time, recalculates the total volume of the sample, and calculates the positive volumetric strain. Simultaneously, update the subordinate load area of ​​the nodes in the concave region. and the unit outward normal vector, which will be used to correct the equivalent confining pressure nodal force. Feedback is sent to the boundary system to maintain consistency in external force updates.

[0113] The updated skeleton coordinates and stress field will be carried into the next loading time step, and the process will iterate until the set number of half-sine wave cycles is completed. .

[0114] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments of the above methods. Any references to memory, storage, databases, or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory. Non-volatile memory can include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in various forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), Rambus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM), etc.

[0115] In this specification, the same or similar parts between the various embodiments can be referred to mutually. Each embodiment focuses on describing the differences from other embodiments. In particular, the descriptions of the embodiments described later are relatively simple, and relevant parts can be referred to the descriptions of the foregoing embodiments.

[0116] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A numerical simulation method for reinforced soil based on morphological evolution and boundary response, characterized in that, The method includes: Acquire point cloud data of granular surface and reconstruct it into a three-dimensional solid model, classify it into parent particles of various typical morphologies and extract representative feature particle size; Obtain the single-particle crushing strength parameters and fragment morphology, and establish the directional mapping relationship between the critical crushing stress threshold of the parent particle and the rigid block template of the actual fragment; A bidirectional geogrid model was established based on representative characteristic particle size, and a flexible boundary was constructed to enclose the particle aggregate and the outside of the geogrid model. The total volume change of the flexible boundary envelope is calculated and fed back to the confining pressure servo control module; Under servo confining pressure, cyclic loads are applied, and the octahedral shear stress on each parent particle is monitored. If the critical fracture stress threshold is reached, the corresponding real fragment rigid block template is generated in the original spatial position according to the directional mapping relationship, and iteratively updated.

2. The numerical simulation method for reinforced soil based on morphological evolution and boundary response as described in claim 1, characterized in that, The particles were classified into several typical forms of parent particles, and representative characteristic particle sizes were extracted. Specifically, this included: obtaining the minimum bounding cuboid of the reconstructed 3D solid model and extracting the lengths of the major axes of the orthogonal parts of the cuboid. Central axis length and minor axis length Representative characteristic particle size is calculated based on the axis lengths of three orthogonal directions. The formula for calculating the representative characteristic particle size is: ； The volume calculated in reverse is equal to the actual volume of the particle. Find the equivalent perfect sphere, and calculate the surface area of ​​the equivalent perfect sphere as its surface area. According to the formula Calculate the sphericity of the particle. ,in This represents the actual surface area of ​​the particle.

3. The numerical simulation method for reinforced soil based on morphological evolution and boundary response as described in claim 2, characterized in that, The establishment of the directional mapping relationship between the critical fragmentation stress threshold of the parent particles and the rigid block template of the actual fragments includes the following conversion criteria: a conversion reference size is introduced based on the geometric classification of the parent particles. For blocky and spindle-shaped particles, the conversion reference size is as follows: Take representative characteristic particle size For strip-shaped and sheet-shaped particles, take the minor axis length of the smallest circumscribed cuboid; the peak failure load in the single-particle crushing test is... Transformed into the ultimate tensile stress of particles at the moment of failure : ； Extracting characteristic intensities to characterize the macroscopic intensity of particles of various shapes and Weber modulus This was used to determine the critical fracture stress threshold for parent particles of various morphologies. .

4. The numerical simulation method for reinforced soil based on morphological evolution and boundary response as described in claim 1, characterized in that, Based on representative characteristic particle size, a bidirectional geogrid model is constructed. Specifically, the maximum value among the representative characteristic particle sizes of various typical particle morphologies is extracted as the effective characteristic particle size. Set the net mesh width of the bidirectional geogrid model. With effective feature particle size It satisfies the following linear matching model: ; in, The mesh matching coefficient is set to a value range of 1. .

5. The numerical simulation method for reinforced soil based on morphological evolution and boundary response as described in claim 1, characterized in that, When constructing the bidirectional geogrid model, a parallel bonding model is introduced; The parallel bonding model is equivalent to the formation of an elastically bonded cylinder at the contact surface of two adjacent particles, and the equivalent maximum normal tensile stress acting at the critical point of the cylinder's cross-section is calculated in real time. and equivalent maximum shear stress The calculation model is as follows: ; ； in, For bonding normal force, For the bonding tangential force, For bond bending moment, For bonding torque; The cross-sectional bonding radius; This represents the cross-sectional area of ​​the bonding section. The moment of inertia of the bonded section; The polar moment of inertia of the bonded section.

6. The numerical simulation method for reinforced soil based on morphological evolution and boundary response as described in claim 1, characterized in that, Constructing a flexible boundary around the particle assembly and the grid model includes an equivalent mapping stage from contact forces to nodal forces: particles within the discrete element computation system collide with the surface elements corresponding to the shell elements in the finite difference system. For the shell elements... The corresponding wall surface element exists The first particle contact point will be the first particle contact point. The contact force vector exerted by the particles on the boundary at each contact point is denoted as . For shell unit nodes The equivalent contact node force vector it receives The calculation formula is: ; in, For the first One contact point in the housing unit Area coordinate weights within the area.

7. The numerical simulation method for reinforced soil based on morphological evolution and boundary response as described in claim 6, characterized in that, Performing dynamic updates of external load application and boundary deformation displacement in a finite difference system: servo confining pressure Converted to apply to each node Equivalent confining pressure nodal force The formula for calculating the equivalent confining pressure nodal force is: ； in, For nodes The subordinate load area, whose value is shared by this node. Sum the current surface areas of all shell elements by adding up one-third of each element. For nodes The current unit outward normal vector of the surface; the equivalent confining pressure nodal forces With equivalent contact node force The external load vector of the node is obtained by summing the vectors, and the unbalanced force is obtained by subtracting the internal force of the shell element node.

8. The numerical simulation method for reinforced soil based on morphological evolution and boundary response as described in claim 1, characterized in that, The specific process for calculating the total volume change of the flexible boundary envelope is as follows: the geometric surface enclosing the particle assembly is considered as a closed set of discrete surface element meshes. The set contains Triangular face elements All face elements The surface elements are all oriented with their normals pointing outwards from the sample. The instantaneous spatial coordinate vectors of the three vertices are respectively , and ; Set the global coordinate origin and the surface element The three vertices are connected to form a tiny triangular pyramid. The directed volume contribution of this tiny triangular pyramid is... Determined by the product of the three position vectors: ; Traversing the grid set All of them The total real-time volume after irregular deformation is obtained by algebraically summing the individual facets. The calculation formula is: .

9. The numerical simulation method for reinforced soil based on morphological evolution and boundary response as described in claim 8, characterized in that, Feedback is sent to the confining pressure servo control module, including: Calculate the real-time volumetric strain that reflects the macroscopic shear response of the soil. : ; in, This represents the total volume of the initial closed mesh. The servo control module retrieves the preset constant target confining pressure. Based on the reconstructed topology, dynamic calculations are performed and each node is updated. The equivalent confining pressure nodal force required at the location The calculation formula is: ； in, To update the node's topology value according to the current time step Subordinate load-bearing area, For the nodes after updating the topology Unit outward normal vector.

10. The numerical simulation method for reinforced soil based on morphological evolution and boundary response as described in claim 3, characterized in that, Monitor the octahedral shear stress on each parent particle. If the critical fracture stress threshold is reached, replace and generate the corresponding real fragment rigid block template in the original spatial position according to the directional mapping relationship. Specifically, iterate through any block of parent particles. Extract the maximum value of the contact force modulus at all contact points. Based on the current conversion reference size of this particle Calculate the octahedral shear stress inside the particle. : ； If subjected to octahedral shear stress Reaching the critical fracture stress threshold corresponding to the parent body shape category Then, retrieve the corresponding subset of real fragment rigid block templates; let the sum of the original volumes of each fragment template within the retrieved template subset be... Apply a uniform linear scaling factor to all fragments within this subset. This makes the total volume of the scaled fragments equal to the parent particle. actual volume linear scaling factor The formula for calculation is: ； Deleting native particles in the numerical space The scaled fragment template set is filled in situ with the original centroid coordinates as the reference point to generate a new pure friction contact pair.

Images