Method for generating three-dimensional polyhedral unit model of complex morphology material containing holes and cracks

Abstract

The invention relates to a method for generating a three-dimensional polyhedral unit model of a hole and crack-containing complex morphology material, and belongs to the technical field of material science, geotechnical engineering and biomechanical modeling. Comprising the following steps: S1, obtaining a three-dimensional digital outer contour of a to-be-modeled material; s2, analyzing geometric parameters of the three-dimensional digital outer contour; and S3, based on the geometric parameters, setting a statistical distribution rule of types and size parameters of the target polyhedral units, and a target volume ratio of the total volume to the outer contour volume of all the polyhedral units. The method has high flexibility, allows to define the morphology of the filling units according to the actual microstructure of the material, controls the size distribution and the filling compactness of the filling units, realizes reasonable distribution of the filling units in a space containing complex defects, can provide different polyhedron filling schemes for different types of models, and has good application prospects. The filling compactness of the polyhedrons can be controlled, and the geometric morphology of the filled polyhedrons and the arrangement statistical rule of the polyhedron group can be set.

Description

Technical Field

[0001] This invention relates to a method for generating a three-dimensional polyhedral element model of a material with complex morphology containing pores and cracks, belonging to the fields of materials science, geotechnical engineering and biomechanical modeling technology. Background Technology

[0002] In materials science, geotechnical engineering, and biomechanics, many engineering materials and natural media (such as ceramic components, fractured rock masses, composite materials, and biological tissues) commonly possess complex discontinuous structures such as pores and fissures. The macroscopic mechanical properties of these materials, including but not limited to their strength, deformation characteristics, and failure modes, are largely controlled by their internal microstructural features. These features mainly include the geometry, characteristic dimensions and their statistical distribution, spatial arrangement, and spatial distribution of pores and fissures of the infilling elements. Therefore, constructing a high-quality three-dimensional polyhedral element model is a crucial prerequisite for ensuring the reliability of numerical simulations.

[0003] Numerical simulation techniques, especially numerical methods based on discontinuous media theory (such as discontinuous deformation analysis and discrete element method), have become key tools for studying the aforementioned mechanical behaviors of materials. However, the reliability of numerical simulation results heavily depends on the quality of the initial discretization model. A high-quality discretization model must meet the following core requirements: it should be able to accurately characterize the internal structure of the material, that is, the filling elements should have realistic geometric morphology, characteristic size distribution, spatial distribution pattern, and overall filling density, and it should be able to accurately reproduce the geometric information of internal defects (such as pores and cracks).

[0004] Currently, the methods for establishing discrete models of three-dimensional polyhedral elements are mainly divided into two categories: curved polyhedrals represented by ellipses and conventional polyhedrals represented by blocks. The main methods for modeling curved polyhedrals are: (1) Rain method: filling regular geometric regions by free fall, suitable for simple shapes such as cubes and spheres, but unable to handle complex hole and crack structures; (2) Compaction method: increasing filling density by mutual compression between sphere elements, but with low efficiency and poor stability under complex boundary conditions, and difficult to accurately control the geometric relationship between sphere elements and crack surfaces; (3) Voxelization layering method: generating models quickly based on spatial discretization, suitable for large volume calculations, but the distribution of sphere elements is uneven, the boundary transition is not smooth, and it is difficult to realistically reproduce the morphology of holes and cracks; (4) Geometric method: generating sphere elements by wavefront method or tetrahedral mesh filling method, which can build sphere element models in complex boundary and irregular geometric regions, and the position of sphere elements can be accurately controlled and can maintain the geometric relationship with the boundary. However, this method is computationally intensive and inefficient when dealing with large volume models, and the spatial density and contact rationality of spherical elements are still difficult to guarantee in complex areas with dense holes and cracks.

[0005] For conventional polyhedral elements, represented by blocks, the spatial mesh discretization method is often used to construct their discrete models. This method decomposes the target region into several continuous, non-overlapping polyhedral elements through uniform division of regular rectangular meshes, tetrahedral discretization, or other spatial partitioning methods, thereby constructing a well-structured discrete model. However, the models generated by this method are usually single, uniform block systems, making it difficult to flexibly control the size distribution, spatial density, and morphological combination of the elements. In particular, it cannot effectively adapt to situations with irregular boundaries such as complex holes and cracks, and it is even more difficult to meet the requirements for constructing discrete models that are not completely dense or have specific gradation characteristics.

[0006] Currently, existing technologies have significant shortcomings in generating discrete models suitable for materials with complex internal structures that can simultaneously meet multiple requirements mentioned above. Specifically, these shortcomings mainly manifest in: difficulty in flexibly specifying the geometric type of the filling elements (such as polyhedra, ellipsoids, etc.); lack of precise control over the statistical distribution of element sizes; inability to effectively guarantee the preset filling density of the overall model and local areas; and difficulty in achieving adaptive optimization of the spatial distribution of filling elements when dealing with complex internal pores and cracks.

[0007] Based on this, the present invention provides a method for generating a three-dimensional polyhedral unit model of a material with complex morphology containing pores and cracks. Summary of the Invention

[0008] In view of this, the present invention provides a method for generating three-dimensional polyhedral element models of materials with complex morphologies containing pores and cracks. This method is highly flexible and allows for the definition of the morphology of the filling elements, control of their size distribution and filling density, and reasonable division within a space containing complex defects, based on the actual microstructure of the material. It can propose different polyhedral filling schemes for different types of models, control the filling density of polyhedra, the set geometric morphology, and the statistical regularity of polyhedral groups. This method is a way to control the generation of high-quality models by controlling multiple parameters in various numerical computation fields that require discrete models.

[0009] This invention provides a method for generating a three-dimensional polyhedral element model of a material with complex morphology containing pores and cracks. The proposed technical solution includes the following steps:

[0010] S1: Obtain the three-dimensional digital outer contour of the material to be modeled. The three-dimensional digital outer contour defines the closed boundary including the surface of internal holes and cracks. Specific acquisition methods include, but are not limited to, three-dimensional CT scanning, MRI imaging, three-dimensional laser scanning, or digital models directly constructed by three-dimensional modeling software.

[0011] S2: Analyze the geometric parameters of the three-dimensional digital outer contour to provide a data basis for setting the distribution rules of polyhedral units; the analyzed geometric parameters include at least one of the following: the total volume of the outer contour, the volume of the holes, the volume of each closed subspace divided by the crack, and the ratio of the volume of each subspace to the total volume, etc.

[0012] S3: Based on the geometric parameters analyzed in step S2, intelligently set the type of the target polyhedral unit, the statistical distribution law of its size parameters, and the target volume ratio of the total volume of all polyhedral units to the volume of the outer contour; the target polyhedral unit is at least one of a sphere, an ellipsoid, or a polyhedron; preferably, different unit size statistical laws can be set according to the volume of different subspaces to achieve adaptive distribution. When the target polyhedral unit is a sphere, its size parameter is the radius; when it is an ellipsoid, its size parameters are the major axis and the minor axis; when it is a polyhedron, its size parameters are the circumscribed sphere radius and the number of faces; the statistical distribution law is a uniform distribution, a normal distribution, or a probability distribution customized based on the geometric parameters;

[0013] S4: Based on the rules set in step S3, determine the number of polyhedral units and generate an initial number and size of polyhedral units. Scale down each initial polyhedral unit proportionally and place it within the three-dimensional digital outer contour. Ensure that all scaled-down initial polyhedral units do not contact each other using a spatial partitioning algorithm or iterative detection algorithm. Initially place the scaled-down initial polyhedral units within the three-dimensional digital outer contour using random point distribution or a uniform grid layout. Ensure that the scaled-down initial polyhedral units do not contact each other using a spatial partitioning algorithm or iterative detection algorithm. The spatial partitioning algorithm is a spatial grid. The partitioning method includes: dividing the three-dimensional digital outer contour into uniform or non-uniform three-dimensional mesh units, and placing at most one reduced initial polyhedral unit in each mesh unit; the spatial partitioning algorithm is a spatial tetrahedral partitioning method, including: performing tetrahedral meshing on the three-dimensional digital outer contour, and placing at most one reduced initial polyhedral unit in each tetrahedral unit; the iterative detection algorithm abstracts a circumscribed sphere boundary for each polyhedral unit and uses a collision detection algorithm to ensure that the circumscribed sphere boundaries do not overlap, thereby achieving that the reduced initial polyhedral units do not contact each other.

[0014] S5: Using the three-dimensional digital outer contour as a fixed constraint boundary, apply the discontinuous deformation analysis algorithm (DDA) and adjust the deployed polyhedral elements to their target size step by step according to the preset scaling rules; scale up the deployed polyhedral elements step by step according to the preset mathematical sequence until all elements reach their target size; the preset mathematical sequence includes an arithmetic sequence, a geometric sequence, or other incremental sequences customized based on model requirements, and the elements are scaled up step by step according to this sequence until they reach their target size;

[0015] The algorithms for analyzing discontinuous deformation include:

[0016] The overall equilibrium equations of the polyhedral unit system are constructed as follows:

[0017]

[0018] Where M is the system mass matrix, C is the damping matrix, K is the stiffness matrix, D is the displacement vector, and F is the load vector; and at each calculation step, the contact equations are dynamically constructed and solved according to the contact state between the polyhedral elements to update the system forces and displacements.

[0019] S6: After the polyhedral elements are enlarged to the target size, it is determined whether the preset model convergence criterion is met. When the convergence criterion is met, the calculation is terminated and the final polyhedral element discrete model is output. The preset model convergence criterion is: all polyhedral elements reach their target size, and the maximum displacement of the entire polyhedral element system is less than a set threshold. Under this state, multiple consecutive calculation steps are maintained to achieve quasi-static equilibrium. When the convergence criterion is met, the calculation is terminated and the final three-dimensional polyhedral element discrete model is output. The multiple consecutive calculation steps are the multiple calculation steps required to ensure that the system reaches the quasi-static equilibrium state.

[0020] Furthermore, the geometric parameters analyzed in step S2 include the volumes of multiple closed subspaces divided by the crack;

[0021] In step S3, based on the volume of the subspace, different statistical rules for the size of the target polyhedral units are set for different subspaces. The larger the volume of the subspace, the larger the average size of the target polyhedral units set for it. In step S4, according to the statistical rules for the size of the target polyhedral units set for each subspace, polyhedral units are generated and arranged to ensure that the ratio of the total volume of the polyhedral units in each subspace to the volume of the subspace is kept near a set target ratio.

[0022] Step S1 is as follows:

[0023] The acquisition of three-dimensional digital contours is mainly achieved through two approaches: First, for physical entities, non-destructive testing technologies such as 3D laser scanning or CT scanning can be used to obtain their surface point cloud and internal volume data. Then, through point cloud processing, surface reconstruction, and internal structure surface recognition algorithms, a complete "watertight" triangular mesh model representing the geometric features inside and outside the entity can be generated. Second, for existing CAD, BIM, or geological 3D numerical models, their boundary representation data can be extracted through data interfaces, and the geometric engine can be used to discretize them into triangular face models to completely preserve the original model's outer contour and internal discontinuities.

[0024] Step S2 is as follows:

[0025] This paper describes the process of calculating geometric parameters and identifying features of complex 3D contours and various polyhedra. For volume calculation of complex 3D contours, two main methods are employed: The first is the 3D mesh discrete integration method. This method uniformly or adaptively subdivides the circumscribed cuboid region of the contour along three orthogonal coordinate axes, generating a dense cubic mesh (voxels). Then, geometric inclusion judgment algorithms such as ray casting or angle summation are used to determine whether the center point of each voxel is located within the closed space defined by the target contour. Finally, the total volume estimate is obtained by accumulating the volumes of all internal voxels, and the calculation accuracy can be controlled by adjusting the voxel size. The second method is a precise algorithm based on the divergence theorem. This method is suitable for models with existing surface triangular meshes. It calculates the directed volume of the pyramid formed by each triangular facet and the origin, and sums the directed volumes of all facets, taking the absolute value to achieve high-precision volume calculation.

[0026] In calculating the volume of polyhedra, regular polyhedra (such as cuboids, spheres, and ellipsoids) can be directly calculated using their analytical volume formulas. For example, the volume formula for an ellipsoid is: For any polyhedron, the internal structure can be divided into several non-overlapping tetrahedrons using the volume triangulation method, and then the summation can be performed. Alternatively, the surface integral formula can be used based on the divergence theorem. Perform general calculations. Here, S is the surface area of ​​the polyhedron, n is the unit normal vector of the outer surface, and r is the position vector of a point on the surface relative to point O.

[0027] In terms of geometric feature identification, regular polyhedra can directly extract feature elements from their definition parameters, such as the three mutually orthogonal principal axis directions and corresponding semi-axis lengths of an ellipsoid; arbitrary polyhedra, on the other hand, can abstract macroscopic geometric features into the principal axis directions and semi-axis lengths of an equivalent ellipsoid by calculating its minimum volume circumscribed ellipsoid, thus providing a standardized geometric description basis for subsequent contact search and motion analysis.

[0028] Step S3 is as follows:

[0029] First, select the polyhedron type based on the material properties of the simulated object: curved polyhedra (represented by ellipsoids) are suitable for materials where surface curvature has a significant impact, while bulk polyhedra (represented by polyhedra) are suitable for granular materials exhibiting interlocking effects. Then, set the statistical distribution law for the characteristic geometric elements of the selected polyhedron, such as uniform distribution, normal distribution, log-normal distribution, or negative exponential distribution, to simulate the heterogeneous characteristics of the material. This is done by setting density parameters. The theoretical porosity of the model is controlled by the ratio of the sum of the volumes of all polyhedra to the total volume of the outer contour. Based on the average volume derived from the total volume of the outer contour, the target density, and the statistical distribution of characteristic dimensions, the required total number of polyhedra can be calculated, and the specific size of each element can be determined through random sampling. For complex outer contours containing subspaces such as small holes and cracks, a normal or negative exponential distribution is preferable to generate sufficiently small elements to improve infill density; for outer contours with uniform structures, a uniform distribution can be used to make the element distribution more balanced.

[0030] Step S4 is as follows:

[0031] The choice of modeling method depends on the density and element type of the target model: when set to fully dense ( =1) When the element is a polyhedron, a spatial discretization method (such as Voronoi meshing) can be used to directly mesh the outer contour. For models that are not completely dense or contain ellipsoidal or other curved surface elements, the general generation method proposed in this invention is used. This method first performs basic partitioning of the space where the complex outer contour is located. Adaptive partitioning strategies such as three-dimensional regular meshes, unstructured tetrahedral meshes, or octrees can be used to establish a spatial index. Then, polyhedral seeds are placed in the partitioned elements. The total number of seeds should not exceed the total number of elements to ensure that they do not overlap initially. If the outer contour contains multiple subspaces (such as holes or cracks), the number of seeds should be allocated according to the proportion of each subspace's volume to the total volume, and a minimum seed threshold should be set for small subspaces to ensure statistical representativeness. The characteristic geometric dimensions of each seed (such as semi-axis length and circumscribed sphere radius) are assigned values ​​according to a set statistical distribution law (such as normal distribution or uniform distribution). Finally, the linear dimensions of all polyhedra are uniformly multiplied by a scaling factor less than 1 to appropriately reduce the size, ensuring that all elements do not touch each other in the initial state, providing a reasonable geometric model for subsequent calculations.

[0032] Step S5 is as follows:

[0033] The discontinuous deformation analysis (DDA) method is used to drive polyhedral elements to gradually expand from an initial shrinkage state to the target size. This process uses a 3D digital outer contour as a fixed constraint boundary and establishes an efficient search algorithm applicable to ellipsoids, blocks, and their contact scenarios to accurately identify the contact relationships between elements. Based on this, force-motion equilibrium equations are constructed according to system dynamics principles, and a staged variable magnification control strategy is adopted to ensure that the elements maintain numerical stability throughout the expansion process, ultimately achieving a reasonable reconstruction of the discrete system.

[0034] 1. The contact search between ellipsoids is as follows:

[0035] We simplify the contact search between ellipsoids to the contact search between spheres, taking the longest axis of the ellipsoid as the radius of the sphere. First, we divide the computational space into a regular grid (called "boxes"), and then quickly assign the spheres and fixed surfaces to these boxes. This transforms the global contact search problem into a problem that is only performed within locally adjacent boxes, greatly improving the search efficiency.

[0036] (1) Space division and box definition

[0037] ① Define the primary computational region: Determine a three-dimensional cuboid region that can contain all moving spheres as the primary computational region. This region is defined by the minimum coordinates. and maximum coordinates definition.

[0038] ② Spatial Mesh Generation (Box Creation): The main computational domain is uniformly divided into cubic mesh cells with side length L. Each cube is called a "box". The total number of boxes, N_box, is calculated by the following formula:

[0039]

[0040] in, , , The same calculation applies. Regions outside the main calculation area are represented by special "boundary boxes".

[0041] ③ Create a 3D array of boxes: Create a 3D array Boxes[ in the program.] ][ ][ These boxes are managed using a linked list. Each box object contains at least one linked list to store the index numbers of the spheres belonging to that box.

[0042] (2) Mapping of the sphere and the box

[0043] Iterate through all spheres in the system. For each sphere... (The coordinates of the center of the sphere are () , , ), radius is ), calculate the index of the box to which it belongs based on its center coordinates ( , , ):

[0044]

[0045] (To ensure the index does not go out of bounds, it is necessary to...) , , Perform boundary processing.

[0046] sphere Add the index number to the Boxes[ ][ ][ The spheres are in a linked list of sphere indices. Each sphere belongs to only one box at any given time.

[0047] (3) The relationship between the fixed surface and the box

[0048] Iterate through all fixed surfaces in the system (usually represented by the equation of an infinite plane). For each fixed surface, find all boxes that may come into contact with that fixed surface in the next time step.

[0049] Judgment criteria: ① For bounding boxes, they are assumed to be associated with all fixed faces. ② For inner boxes (within the main calculation region), the distance from the center point of the box to the fixed faces is calculated. .like <= / 2+ + (in It is the maximum radius of the sphere in the system. If the estimated maximum displacement at the next time step is given, then the box is considered to be associated with this fixed surface. Each fixed surface is linked to all its associated boxes (e.g., by storing pointers or indices of associated boxes in a linked list).

[0050] (4) Box-based contact pair search

[0051] ① Ball-to-ball contact search: Traverse each Boxes[ ][ ][ For each sphere A within the box, its potential contact objects are limited to: the same boxes [Boxes]. ][ ][ The other spheres within; and its 26 adjacent boxes.

[0052] ② Ball-fixed-face contact search: Traverse each fixed face. For each fixed face, simply traverse all boxes associated with it (the result of step (3)). For each ball in each associated box, determine whether the ball may be in contact with the fixed face. This avoids traversing all balls and all fixed faces.

[0053] 2. The contact search between blocks is as follows:

[0054] The contact detection between blocks employs a contact judgment method specifically for three-dimensional polyhedra. This method systematically completes the entire process from rapid screening and accurate judgment to determining the contact surface by establishing criteria such as "unit bounding box," "global direct contact search," and "final intrusion point determination."

[0055] It mainly consists of three steps:

[0056] (1) Fast filtering based on unit bounding boxes

[0057] This step aims to quickly identify potentially contacting block pairs and their boundary geometric elements (points, edges, faces). ① Define bounding boxes: Create an axial bounding box for each surface of the block, serving as the reference bounding element for that surface. These elements can also approximately cover the edges and vertices contained in the surface. ② Global coarse screening (finding neighboring blocks): Divide the computational domain into a uniform grid, marking the grid occupied by the global bounding box of each block. If the marked grids of two blocks overlap, they are initially identified as "neighboring block pairs". Local coarse screening (finding potentially contacting geometric elements): For each pair of "neighboring block pairs", check the overlap between all their bounding boxes. ③ Judgment criteria: If there are at least three overlapping bounding boxes between two blocks, they are considered to be in contact. By analyzing the overlapping bounding box pairs, the specific types of potentially contacting geometric elements can be further identified, such as face-to-face, edge-to-edge, point-to-face, etc.

[0058] (2) Global direct contact search and precise judgment

[0059] This step performs precise intrusion analysis and contact information calculation on the selected geometric element pairs. ① Search for potential contact points: Avoiding pre-assuming contact types or locality, a global search is conducted to calculate all possible contact point types, mainly including "vertex-face" contact and "intersecting edge-edge" contact. For each candidate contact point, its initial distance, relative approach velocity, and non-intrusion conditions are calculated. ② Determine the actual contact location: Among all valid candidate contact points, the location of the first contact and its corresponding contact surface are determined according to the "final intrusion point" criterion. This criterion states that the location where two blocks first intrude corresponds to the contact surface defined by the last strictly valid contact point that moves from "non-contact" to "intrusion state". By comparing the "intrusion time" of each contact point and combining it with the relative motion trend of the blocks, the contact surface is accurately determined.

[0060] (3) Contact force calculation

[0061] Based on the specific contact points, contact surfaces and intrusion information determined in step (2), the normal contact force and tangential friction force between the blocks are calculated using the penalty function method or the Lagrange multiplier method, and these forces are substituted into the control equations to solve the motion of the block system.

[0062] 3. The contact search between the ellipsoid and the block is as follows:

[0063] The core of contact search between an ellipsoid and a block lies in the reasonable simplification of the geometrically complex block to utilize an efficient spherical contact search framework. This simplification is achieved by constructing a minimum circumsphere for the block: first, the geometric center of the block is determined; then, the distances from this geometric center to each vertex on the block's surface are calculated, and the maximum value is taken as the radius of the circumsphere. The constructed circumsphere can completely enclose the block.

[0064] 4. Magnification setting:

[0065] The variable magnification is set in stages: starting from the reduced state, it is gradually magnified in the sequence of 3 / 20, 5 / 20, 7 / 20...15 / 20 of the target size until the target size is reached.

[0066] 5. The DDA-based drive and contact force matrix are (changed to three degrees of freedom, rigid body):

[0067] The method for analyzing discontinuous deformation based on polyhedral elements calculates the displacement increment of polyhedral elements, specifically including the following steps:

[0068]

[0069] in, It is a 3×3 stiffness submatrix. With polyhedral units The attribute correlation matrix With the ball and It relates to the interaction between them; Polyhedral unit The fundamental unknown quantum matrix, denoted as ,in Polyhedral unit Displacement in the x-direction, Polyhedral unit Displacement in the y direction, Polyhedral unit Displacement in the z-direction; Polyhedral unit The generalized force matrix, denoted as ,in Polyhedral unit Force in the x direction, Polyhedral unit Force in the y-direction, Polyhedral unit Force in the z-direction;

[0070] Using the principle of minimum potential energy, the contact matrix between polyhedral elements can be obtained, expressed by the following formula:

[0071]

[0072] in, Polyhedral unit With polyhedral units The normal spring stiffness; Polyhedral unit The unit vector matrix between the geometric centers of the polyhedral elements; Polyhedral unit With polyhedral units Normal distance between;

[0073] Cracks, holes, and the outer contour of the model that come into contact with the polyhedral elements are simulated using triangular faces. The contact sub-matrix between the polyhedral elements and the triangular faces is then constructed using the principle of minimum potential energy, and expressed by the following formula:

[0074]

[0075] in, Polyhedral unit With plane The normal spring stiffness; Polyhedral unit With plane The direction vector matrix between contact points; Polyhedral unit With plane Normal distance between contact points;

[0076] The obtained sub-matrices (sub-matrices of contact between polyhedral elements and sub-matrices of contact between polyhedral elements and planes) are superimposed onto the overall equilibrium equation, and the overall equilibrium equation is solved. By solving the overall equilibrium equation, the displacement component of each polyhedral element can be obtained, which is used as the displacement increment of the polyhedral element.

[0077] Furthermore, the discontinuous deformation analysis method based on polyhedral elements described above, which calculates the displacement increment of polyhedral elements, can also construct and solve equations for polyhedral elements with four, six, and seven degrees of freedom. The principle is consistent with the equation construction and solution method for the three degrees of freedom mentioned above; the three degrees of freedom are translations in the X, Y, and Z directions, and the four degrees of freedom are translations in the X, Y, and Z directions and radius. The six degrees of freedom are translation in the X, Y, and Z directions, rotation in the XY plane, rotation in the XZ plane, and rotation in the YZ plane. The seven degrees of freedom are translation in the X, Y, and Z directions, rotation in the XY plane, rotation in the XZ plane, rotation in the YZ plane, and radius. .

[0078] Step S6 is as follows:

[0079] The convergence criteria include the following three items: First, the displacement increment of all polyhedra within a single calculation step must not exceed the set maximum displacement increment threshold; second, the polyhedral system must meet the above displacement increment threshold requirement in multiple consecutive calculation steps; finally, all polyhedra must have been enlarged to the preset original size. When the model simultaneously meets the above criteria, the model generation is considered complete.

[0080] This invention discloses a method for rapidly and efficiently establishing a discrete model of a complex three-dimensional outer contour containing pores and cracks composed of specified polyhedra. It possesses high flexibility, allowing the definition of the morphology of the filling units, control of their size distribution and filling density based on the actual microstructure of the material, and achieving a reasonable distribution of these units within a space containing complex defects. It can propose different polyhedral filling schemes for different types of models, control the filling density of the polyhedra, the set geometric morphology, and the statistical regularity of the polyhedral group. This method is applicable to various numerical computation fields that require discrete models, enabling the control of multiple parameters to generate high-quality models. Attached Figure Description

[0081] Figure 1 This is a flowchart illustrating the method of the present invention.

[0082] Figure 2 This is a geometric outline diagram of a three-dimensional complex morphology material model drawn in an embodiment of the present invention without holes.

[0083] Figure 3 This is a diagram of the initial three-dimensional spherical unit model generated in the hole-free embodiment of the present invention.

[0084] Figure 4 This is a three-dimensional spherical unit model diagram after processing the magnified spherical unit using the spherical unit DDA method in the hole-free embodiment of the present invention.

[0085] Figure 5 This is a geometric outline diagram of a three-dimensional complex morphology material model drawn in an embodiment of the present invention with holes.

[0086] Figure 6 This is a diagram of the initial three-dimensional spherical unit model generated in the embodiment with holes of the present invention.

[0087] Figure 7This is a three-dimensional spherical unit model diagram after processing the magnified spherical unit using the spherical unit DDA method in the embodiment of the present invention with holes.

[0088] Figure 8 This is a schematic diagram of a three-dimensional spherical unit model obtained using the method of this invention in a practical application.

[0089] Figure 9 Charts showing common distribution types and their characteristics. Detailed Implementation

[0090] The preferred embodiments of the present invention will now be described in detail.

[0091] This invention provides a method for generating a three-dimensional polyhedral element model of a material with complex morphology containing pores and cracks, the flowchart of which is shown below. Figure 1 As shown, it includes the following steps:

[0092] S1. Acquisition of 3D Digital Outer Contour: Acquire the 3D digital outer contour of the material to be modeled. The 3D digital outer contour defines the closed boundary including the surface of internal holes and cracks. Specific acquisition methods include, but are not limited to, 3D CT scanning, MRI imaging, 3D laser scanning, or digital models directly constructed by 3D modeling software.

[0093] Step S1 is as follows:

[0094] (1) 3D reconstruction of physical entities: First, non-destructive testing technologies such as 3D laser scanning or computed tomography are used to obtain high-precision surface point cloud data of the entity; then, point cloud registration, denoising, and surface reconstruction algorithms are used to generate a triangular mesh model representing the outer contour of the entity. For discontinuous structural surfaces such as cracks and stratification inside the entity, internal recognition algorithms such as image segmentation and edge detection are used to extract and identify them from its volume data (such as CT sequences), and they are also discretized into triangular patches. Finally, the triangular mesh of the outer contour and the internal structural surface are fused and topologically optimized to form a complete "watertight" triangular patch model representing the geometric features inside and outside the entity;

[0095] (2) Format conversion for existing numerical models: For 3D numerical models already existing in computer-aided design, building information modeling, or geological modeling software, the boundary representation data of the model is extracted through its application programming interface or standard data exchange format. This process requires complete acquisition of the model's outer contour geometry and the various discontinuous interfaces defined within it (such as fault planes, joint surfaces, material interfaces, etc.). Subsequently, these B-Rep data are discretized into triangular patches with high quality through a geometry engine, ultimately obtaining a standardized numerical contour model with triangular patches as basic elements;

[0096] S2. Geometric Parameter Analysis: Analyze the geometric parameters of the three-dimensional digital outer contour to provide a data foundation for subsequently setting the distribution rules of polyhedral elements. The analyzed parameters include the total volume of the outer contour, the volume of the holes, the volume of each closed subspace divided by the cracks, and the ratio of the volume of each subspace to the total volume.

[0097] Step S2 is as follows:

[0098] (1) Volume calculation method for complex three-dimensional outer contours:

[0099] ① Three-dimensional mesh discrete integration method: This is the core algorithm. Its basic idea is to uniformly or adaptively subdivide the circumscribed cuboid space of the target model along three orthogonal coordinate axes, generating a dense cubic mesh (also called voxels). Then, an efficient geometric inclusion judgment algorithm (such as the ray casting method, angle summation method, etc.) is used to determine whether the center point of each voxel is located inside the closed space defined by the target's outer contour. Finally, by accumulating the volumes of all voxels determined to be "inside," the total volume estimate of the complex outer contour is obtained. This method directly controls the computational accuracy by adjusting the mesh density (i.e., voxel size): the denser the mesh (the smaller the voxel size), the closer the volume estimate is to the true value, but the corresponding computational resource consumption is also greater.

[0100] ② This efficient and accurate method can be used for entities whose surface triangular mesh models have been obtained through 3D scanning or CAD modeling. Its principle is based on the divergence theorem (Gauss's formula), transforming volume calculation into an integral operation over the closed triangular mesh surface. Specifically, the directed volume of the pyramid formed by each triangular facet and a fixed point in space (such as the origin) is calculated, and the absolute value of the sum of the directed volumes of all triangular facets is taken to obtain the accurate volume of the entity. This method does not require internal mesh generation, directly utilizing the boundary representation information of the model, resulting in high computational efficiency. Its accuracy depends entirely on how closely the surface triangular mesh approximates the geometry of the entity.

[0101] (2) Methods for calculating the volume of various polygons:

[0102] For calculating the volume of different types of polyhedra, this invention provides the following accurate calculation methods:

[0103] ① Calculation of the volume of a regular polyhedron

[0104] For polyhedra with standard geometric shapes, such as cuboids, spheres, and ellipsoids, their analytical volume formulas can be directly applied for calculation. Taking an ellipsoid as an example, if the lengths of its three mutually perpendicular semi-axis are respectively... , , Then its volume From the formula: Provided.

[0105] This type of method is the most efficient in computation and has perfect accuracy, but it is only applicable to specific regular shapes.

[0106] ② Calculation of the volume of any polyhedron

[0107] For any polyhedron enclosed by several polygonal faces, a general method based on spatial decomposition and vector analysis is employed. The preferred method is volume triangulation: the interior of the polyhedron is divided into several non-overlapping tetrahedron sets, and the total volume is obtained by calculating the volume of each tetrahedron and summing them. The volume of a single tetrahedron can be calculated from the coordinates of its four vertices. Another efficient method is the application of the divergence theorem (Gauss's formula): the volume integral is transformed into a surface integral. Choosing a fixed point in space (such as the origin O), the volume V can be accurately calculated using the following formula:

[0108] in, The surface area of ​​the polyhedron. Let be the unit normal vector of the outer surface. Let be the position vector of a point on the surface relative to point O. This formula holds for any polyhedron. When the surface of the polyhedron has been discretized into triangular facets, the integral can be transformed into a summation operation over all triangular facets, which is convenient for programming implementation.

[0109] (3) Methods for identifying geometric feature elements of various polygons:

[0110] To effectively represent and classify polyhedra geometrically, it is necessary to identify their key geometric feature elements. The specific method is as follows:

[0111] Identification of characteristic elements of regular polyhedra:

[0112] For regular polyhedra such as ellipsoids with standard analytical expressions, their characteristic elements can be directly extracted from their definition parameters. Taking an ellipsoid as an example, its characteristic elements are its three mutually orthogonal principal axes and the lengths of the semi-axes corresponding to each principal axis. Typically, the longest axis is defined as the major axis, and its length is denoted as... The shortest axis is defined as the minor axis, and its length is denoted as . The central axis is defined as the central axis, and its length is denoted as . The ratio of the major axis to the minor axis ( ) is a key shape factor that describes its flatness or elongation.

[0113] ② Identification of characteristic elements of any polyhedron:

[0114] For a polyhedron of arbitrary shape, its characteristic elements are determined based on its smallest circumscribed ellipsoid. This ellipsoid is the only ellipsoid that uniquely encloses a given polyhedron and has the smallest volume, optimally approximating the polyhedron's spatial distribution and orientation. By calculating the three principal axes and semi-axis lengths of this smallest circumscribed ellipsoid, the macroscopic geometric features of any polyhedron (including spatial orientation, approximate size, and shape proportions) can be abstracted and simplified into an equivalent ellipsoid. This abstraction provides a standardized geometric description basis for subsequent contact search, kinematic analysis, and other processes.

[0115] S3. Setting rules for generating polyhedral elements: Based on the geometric parameters analyzed in step S2, intelligently set the statistical distribution rules (e.g., uniform distribution, normal distribution, or custom distribution based on geometric parameters) of the target polyhedral element type (e.g., sphere, ellipsoid, or polyhedron), size parameters (e.g., radius, major and minor axes, circumscribed sphere radius, and number of faces), and the ratio of the total volume of all polyhedral elements to the target volume of the outer contour. Preferably, different element size statistical rules can be set for different subspaces to achieve adaptive distribution.

[0116] Step S3 is as follows:

[0117] (1) Selection of polyhedron type

[0118] Based on the geometric characteristics of the simulated object, two main types of polyhedra can be selected as the basic units:

[0119] ① Curved polyhedra: Represented by ellipsoids, their surfaces are smooth and continuous, making them suitable for simulating materials with rheological properties or whose contact behavior is significantly affected by curvature.

[0120] ② Block-type polyhedra: Represented by polyhedra (such as cubes, chamfered octahedrons, etc.), their surfaces are composed of planes and have edges and corners. They are suitable for simulating granular materials or fractured rock masses with obvious interlocking and locking effects.

[0121] (2) To simulate the heterogeneity of natural or artificial materials, statistical distribution rules need to be set for the characteristic geometric elements of the selected polyhedron (such as the semi-axis length of the ellipsoid, the circumscribed radius of the polyhedron, etc.). Common distribution types and their characteristics are as follows: Figure 9 As shown in the table:

[0122] (3) Setting the model compaction parameters

[0123] Set the compaction parameters of the model It is defined as the ratio of the sum of the volumes of all polyhedra to the total volume of the complex outer contour. By adjusting The value (usually less than 1 to reserve pore space) can control the theoretical porosity of the generated model, thereby indirectly reflecting the initial compaction state or degree of cementation of the material.

[0124] (4) Calculation of the total number of polyhedra and individual dimensions

[0125] Based on the above settings, the system parameters are determined through the following steps: Based on the total volume of the complex outer contour... Target density And the statistical distribution of the selected feature geometric elements (from which the average volume of a single polyhedron can be derived). ), calculate the total number of polyhedra required. Based on the total number N and the established statistical distribution rules, the specific geometric dimensions of each polyhedron are generated through random sampling, thereby uniquely determining its geometric shape.

[0126] (5) Intelligent selection strategy

[0127] For complex outer contours, the following intelligent selection strategies can be used to optimize the model: ① For complex outer contours with "subspaces" such as small holes and cracks: it is advisable to select feature geometric elements with a normal or negative exponential distribution. Such distributions can generate a sufficient number of small-sized polyhedra, thereby effectively filling small spaces and improving the model's spatial filling accuracy and geometric approximation accuracy. ② For outer contours with uniform internal structures and no significant subspaces: it is advisable to select feature geometric elements with a uniform distribution. This distribution can make the distribution of polyhedra of different sizes more uniform, avoid local over-concentration, and better reflect the macroscopic homogeneous characteristics of the material.

[0128] S4. Initial Polyhedral Element Layout: Based on the rules set in step S3, determine the number, size, and shape of the polyhedral elements, generating a corresponding number of initial polyhedral elements. After scaling down all initial elements proportionally, place them within the three-dimensional digital outer contour using methods such as random point scattering or uniform mesh layout. During this process, algorithms such as spatial mesh generation, spatial tetrahedral generation, or iterative collision detection based on circumscribed spheres are used to ensure that all scaled-down initial elements do not contact each other in space.

[0129] Step S4 is as follows:

[0130] Criteria for selecting modeling methods

[0131] The choice of modeling method depends on the density of the target model and the type of basic unit:

[0132] ① High-density regular block model: When the model is set to fully dense (i.e., high density) =1) When the basic unit is a polyhedral block, the complex outer contour can be directly subdivided by spatial discretization method (such as Voronoi subdivision) to quickly obtain a dense discrete block system without initial pores.

[0133] ② Incompletely compact or special element models: When the set compactness is less than 1 (e.g., When the porosity is 0.7 (used to simulate loose media such as gravel piles), or when the basic unit is a curved surface such as an ellipsoid, or a hybrid model of ellipsoids and polyhedra, the spatial discretization method is no longer applicable. Therefore, this invention proposes the following general generation method to achieve model construction with controllable porosity and random distribution characteristics.

[0134] (2) The general generation method of the present invention

[0135] The core of this method lies in decoupling seed placement from geometric entity generation, achieving a controllable follower mechanism model through steps such as spatial partitioning, probabilistic seed placement, and size control. Its main steps include:

[0136] (1) Spatial basic division

[0137] First, the space containing the complex outer contour needs to be fundamentally divided to lay the foundation for precise seed location and statistical scoring. Optional division strategies include:

[0138] ① 3D Regular Mesh Generation: The circumscribed cubic space of the outer contour is uniformly subdivided along the three directions of the Cartesian coordinate system to generate a regular cubic mesh. This method is logically simple and computationally efficient.

[0139] ② 3D unstructured mesh generation: Discretize the internal region of a complex outer contour into a series of non-overlapping tetrahedral elements. This method can better fit complex boundaries, but the meshing algorithm is relatively complex.

[0140] ③ Other adaptive partitioning methods: such as spatial indexing structures like octrees, which can achieve adaptive partitioning of different regions of the model, balancing efficiency and boundary fitting accuracy.

[0141] Polyhedral seed placement and distribution control

[0142] After completing the spatial division, "seeds" (i.e., their positioning points) of polyhedra are placed into the resulting units.

[0143] The following key control principles must be followed during the deployment process:

[0144] ① Total quantity control: The preset total number of polyhedral seeds The number of cells must be less than or equal to the total number of cells generated by the basic partitioning to ensure that each seed can be placed in a unique spatial cell and to avoid initial overlap.

[0145] ② Subspace Volume Ratio: If there are multiple connected or closed subspaces (such as holes or a network of fissures) inside a complex outer contour, the total number of seeds needs to be adjusted according to the proportion of the volume of each subspace to the total volume of the outer contour. Distribute proportionally to each subspace. That is, the first... Number of seeds for subspace allocation ,in Vtotal is the volume of the subspace, and Vtotal is the total volume of the outer contour.

[0146] ③ Minimum number guarantee: To avoid insufficient statistical representativeness due to too few seeds in small subspaces, a minimum seed number threshold should be set for each subspace. When the number of seeds calculated according to the volume ratio is lower than this threshold, the seed number of that subspace should be forcibly set to the threshold.

[0147] (3) Geometric dimension setting and initial non-contact guarantee

[0148] ① Size setting: Assign each seed a characteristic geometric dimension of a polyhedron (such as the semi-axis length of an ellipsoid, the radius of the circumscribed sphere of the polyhedron, etc.), the value of which can be randomly generated according to the set statistical distribution law (such as normal distribution, uniform distribution, etc.).

[0149] ② Initial Non-Contact Processing: To ensure that the generated computational model is free from geometric intrusion in its initial state, all polyhedra, after generation, need to have their linear dimensions (such as radius) uniformly multiplied by a scaling factor less than 1 (e.g., 0.95~0.99) to appropriately reduce their size. This operation effectively ensures that the randomly constructed polyhedra near the seed placement location maintain a small initial gap, achieving an initial non-contact state of the model and providing reasonable initial conditions for subsequent mechanical calculations.

[0150] S5. Using the three-dimensional digital outer contour as a fixed constraint boundary, apply a discontinuous deformation analysis algorithm to gradually enlarge the deployed polyhedral elements to their target size.

[0151] Step S5 is as follows:

[0152] Using the three-dimensional digital outer contour as a fixed constraint boundary, the discontinuous deformation analysis (DDA) algorithm is applied to gradually adjust the deployed polyhedral elements to their target size according to a preset scaling rule. The core of this process is to construct and solve the dynamic equations of the system, and its overall equilibrium equations are in the form of:

[0153]

[0154] Where M is the system mass matrix, C is the damping matrix, K is the stiffness matrix, D is the displacement vector, and F is the load vector. At each calculation step, the contact equations are dynamically constructed and solved based on the contact state between elements, updating the system's forces and displacements.

[0155] S6. Model Convergence Judgment and Output: During the amplification process in step S5, it is determined in real time whether the preset model convergence criteria are met. These criteria include: all polyhedral elements reach their target size, the maximum displacement of the entire system is less than a set threshold, and multiple consecutive calculation steps are maintained in this state to achieve quasi-static equilibrium. When the convergence criteria are met, the calculation is terminated, and the final three-dimensional polyhedral element discrete model is output.

[0156] Step S6 is as follows:

[0157] To ensure that the generated polyhedral system reaches a stable initial state that meets the preset geometric requirements, a clear convergence criterion needs to be established. The convergence criterion used in this method is as follows:

[0158] (1) Criteria for meeting geometric morphology standards

[0159] This criterion is used to verify whether the polyhedron has reached the preset geometric dimensions. During the simulation, the polyhedron gradually recovers from its initial shrunk state (to ensure no contact) to the set original dimensions. When the linear dimension magnification of all polyhedra reaches the preset target value (e.g., 1.0, i.e., completely recovered to the original size), the geometry is considered to meet the standard.

[0160] (2) Criterion for motion stability

[0161] This criterion is used to determine whether the macroscopic motion of a polyhedral system tends to come to rest. ① Maximum displacement increment threshold: Set a very small displacement threshold. During the simulation, it is required that the maximum displacement increment of all polyhedra in the system within a single calculation step must not exceed this threshold. ② Continuous stability step requirement: In addition to meeting the maximum displacement increment threshold, the polyhedral system must further be able to maintain continuous stability steps. Within a computation time step (e.g.) (e.g., values ​​between 50 and 100) must all meet this threshold condition. This requirement aims to eliminate instantaneous pseudo-static states of the system and ensure the continuity of stable motion.

[0162] (3) Comprehensive convergence judgment

[0163] The numerical model can be considered successfully generated and can be used as the initial state for subsequent mechanical analysis only if the polyhedral system simultaneously satisfies the above-mentioned motion stability criterion and geometric shape compliance criterion.

[0164] The method of the present invention will be further illustrated below with examples of a complex morphology without holes and a complex morphology with holes:

[0165] 1. Hole-free example:

[0166] Importing triangular facet data from a non-porous 3D complex-morphology material model to form the model outline, such as... Figure 2 As shown;

[0167] The density of the spherical elements was determined to be 0.92, and the final volume occupied by the spherical elements was 92% of the total model volume. The spherical elements were uniformly distributed, with radii ranging from 0.05m to 0.1m. Ultimately, 351,452 initial spherical elements were generated within the model, with radii ranging from 0.02m to 0.04m. Figure 3 As shown;

[0168] Increase the radius of the initial spherical elements, ensuring overlap between spherical elements and between spherical elements and the model's geometric contour. Calculate the spherical element positions based on the spherical element dynamics analysis (DDA). Move the spherical elements to eliminate overlap. Repeat this process until the spherical element radius reaches 0.05m~0.1m. Then, using the spherical element DDA method, calculate the displacement increment of the spherical elements and update their positions. When the calculated displacement increment of the spherical elements within a single time step is lower than a preset unique increment, the calculation converges, the calculation stops, and the final model is output, as shown below. Figure 4 As shown.

[0169] 2. Example with holes

[0170] Import the triangular facet data of a 3D material model with a complex morphology and holes to form the model outline, such as... Figure 5 As shown;

[0171] The density of the spherical elements was determined to be 0.92, and the final volume occupied by the spherical elements was 92% of the total model volume. The spherical elements were uniformly distributed, with radii ranging from 3m to 4m. Ultimately, 21,918 initial spherical elements were generated within the model, with radii ranging from 1.92m to 2.56m. Figure 6 As shown;

[0172] Increase the radius of the initial spherical elements, ensuring overlap between spherical elements and between spherical elements and the model's geometric contour. Calculate the spherical element positions based on spherical element dynamics analysis (DDA), and move the spherical elements to eliminate overlap. Repeat this process until the spherical element radius reaches 3m~4m. Then, using the spherical element DDA method, calculate the displacement increment of the spherical elements and update their positions. When the calculated displacement increment of the spherical elements within a single time step is lower than a preset unique increment, the calculation converges, the calculation stops, and the final model is output, as shown below. Figure 7 As shown.

[0173] The present invention and its embodiments have been described above. This description is not restrictive, and the accompanying drawings are only one embodiment of the present invention; the actual structure is not limited thereto. In conclusion, if those skilled in the art are inspired by this description and design similar structures and embodiments without departing from the spirit of the invention, such designs should fall within the protection scope of the present invention.

Claims

1. A method for generating a three-dimensional polyhedral element model of a material with complex morphology containing pores and cracks, characterized in that, Includes the following steps: S1: Obtain the three-dimensional digital outer contour of the material to be modeled, wherein the three-dimensional digital outer contour defines a closed boundary including the internal pores and crack surfaces; S2: Analyze the geometric parameters of the three-dimensional digital outer contour; S3: Based on the geometric parameters, set the type of the target polyhedral unit, the statistical distribution law of the size parameters, and the target volume ratio of the total volume of all polyhedral units to the outer contour volume; S4: According to the rules set in S3, determine the number of polyhedral units and generate an initial number and size of polyhedral units. After each initial polyhedral unit is scaled down proportionally, it is placed in the three-dimensional digital outer contour. The spatial partitioning algorithm or iterative detection algorithm is used to ensure that all scaled-down initial polyhedral units do not contact each other. S5: Using the three-dimensional digital outer contour as a fixed constraint boundary, apply a discontinuous deformation analysis algorithm and adjust the deployed polyhedral units to their target size step by step according to the preset amplification rules. S6: After the polyhedral element is enlarged to the target size, it is determined whether the preset model convergence criterion is met; when the convergence criterion is met, the calculation is terminated and the final polyhedral element discrete model is output.

2. The method for generating a three-dimensional polyhedral element model of a material with complex morphology containing pores and cracks according to claim 1, characterized in that: In step S1, the three-dimensional digital outer contour of the material to be modeled is obtained by any one of the following methods: three-dimensional CT scanning, MRI imaging, three-dimensional laser scanning, or a digital model constructed by three-dimensional modeling software.

3. The method for generating a three-dimensional polyhedral element model of a material with complex morphology containing pores and cracks according to claim 1, characterized in that: In step S2, the analyzed geometric parameters include at least one of the following: the total volume of the outer contour, the volume of the hole, the volume of each closed subspace divided by the crack, and the ratio of the volume of each subspace to the total volume.

4. The method for generating a three-dimensional polyhedral element model of a material with complex morphology containing pores and cracks according to claim 1, characterized in that: In step S3, the target polyhedral unit is at least one of a sphere, an ellipsoid, or a polyhedron; when the target polyhedral unit is a sphere, its size parameter is the radius; when it is an ellipsoid, its size parameters are the major axis and the minor axis; when it is a polyhedron, its size parameters are the circumscribed sphere radius and the number of faces.

5. The method for generating a three-dimensional polyhedral element model of a material with complex morphology containing pores and cracks according to claim 1, characterized in that: The statistical distribution pattern is a uniform distribution, a normal distribution, or a probability distribution defined based on the geometric parameters.

6. The method for generating a three-dimensional polyhedral element model of a material with complex morphology containing pores and cracks according to claim 1, characterized in that: In step S4, the proportionally reduced initial polyhedral units are initially placed within the three-dimensional digital outer contour using a random point-scattering or uniform grid layout method; in step S4, a spatial partitioning algorithm or an iterative detection algorithm is used to ensure that the reduced initial polyhedral units do not contact each other. The spatial partitioning algorithm is a spatial mesh partitioning method, which includes: dividing the three-dimensional digital outer contour into uniform or non-uniform three-dimensional mesh units, and placing at most one of the reduced initial polyhedral units in each mesh unit; The spatial partitioning algorithm is a spatial tetrahedral partitioning method, which includes: dividing the three-dimensional digital outer contour into tetrahedral meshes, and placing at most one of the reduced initial polyhedral units in each tetrahedral unit; The iterative detection algorithm abstracts a circumscribed sphere boundary for each polyhedral unit and employs a collision detection algorithm to ensure that the circumscribed sphere boundaries do not overlap, thereby achieving the goal of preventing the reduced initial polyhedral units from contacting each other.

7. The method for generating a three-dimensional polyhedral element model of a material with complex morphology containing pores and cracks according to claim 1, characterized in that: In step S5, the deployed polyhedral units are enlarged step by step according to a preset mathematical sequence until all units reach their target size. The preset mathematical sequence includes an arithmetic sequence, a geometric sequence, or other incremental sequences customized based on model requirements. The units are enlarged step by step according to this sequence until they reach their target size. The multiple consecutive calculation steps are the multiple calculation steps required to ensure that the system reaches a quasi-static equilibrium state.

8. The method for generating a three-dimensional polyhedral element model of a material with complex morphology containing pores and cracks according to claim 3, characterized in that: The geometric parameters analyzed in step S2 include the volume of multiple closed subspaces divided by the crack; In step S3, based on the volume of the subspace, different statistical rules for the size of the target polyhedral unit are set for different subspaces. The larger the volume of the subspace, the larger the average size of the target polyhedral unit set for it. In step S4, polyhedral units are generated and deployed according to the statistical law of the target polyhedral unit size set for each subspace, so as to ensure that the ratio of the total volume of the polyhedral units in each subspace to the volume of the subspace is kept near a set target ratio.

9. The method for generating a three-dimensional polyhedral element model of a material with complex morphology containing pores and cracks according to claim 7, characterized in that: In step S5, the application of the discontinuous deformation analysis algorithm includes: The overall equilibrium equations of the polyhedral unit system are constructed as follows: ； Where M is the system mass matrix, C is the damping matrix, K is the stiffness matrix, D is the displacement vector, and F is the load vector; and at each calculation step, the contact equations are dynamically constructed and solved according to the contact state between the polyhedral elements to update the system forces and displacements.

10. The method for generating a three-dimensional polyhedral element model of a material with complex morphology containing pores and cracks according to claim 1, characterized in that: The preset model convergence criterion in step S6 is: all polyhedral elements reach their target size, and the maximum displacement of the entire polyhedral element system is less than a set threshold, and multiple consecutive calculation steps are maintained under this state.

Images