A method, system and storage medium for modeling three-dimensional ground forces

Abstract

The present application relates to the technical field of three-dimensional graphics data processing, and particularly relates to a three-dimensional ground stress modeling method, system and storage medium. The method comprises the following steps: constructing a three-dimensional space grid model of a geological structure, and marking the spatial parameters of stratum, fault and fissure distribution; inputting the spatial parameters into a preset probabilistic mechanics model to generate a multi-scale fissure field; embedding the multi-scale fissure field into an initial tensor field of a preset three-dimensional stress field to generate continuous stress distribution data containing a crack propagation path through fracture evolution simulation; constructing a numerical updating strategy of fluid-solid coupling based on the continuous stress distribution data to solve the fluid pressure field and the solid stress field in a coupled manner and update the three-dimensional stress field. Through the fusion of multi-scale fissure analysis, fracture evolution tracking and fluid-solid coupling solution, the present application realizes the dynamic reconstruction and three-dimensional visualization of the stress field of the geological structure, and significantly improves the accuracy and expression ability of the ground stress distribution under complex geological conditions.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional graphics data processing technology, and in particular to a three-dimensional geostress modeling method, system and storage medium. Background Technology

[0002] Three-dimensional geostress modeling is a core tool in geological engineering and underground energy development, and its development has evolved from empirical estimation to numerical simulation. Early methods mainly relied on borehole testing and seismic inversion data to establish regional stress fields through linear elasticity theory, but they were difficult to characterize the mechanical effects of complex discontinuous structures such as faults and fractures. With the advancement of computational mechanics and three-dimensional modeling techniques, finite element, discrete element, and boundary element methods have been introduced into geostress simulation, which can predict the spatial distribution of stress fields to a certain extent. However, traditional models usually assume that fracture distribution is fixed and static, lacking dynamic characterization of fracture propagation and evolution; at the same time, the coupling between fluid pressure and solid stress often uses simplified assumptions, resulting in a lag in the fluid-solid interaction response, making it difficult to reflect the stress reconstruction process under real geological conditions. In addition, the visualization of existing models is mostly limited to the scalar field level, failing to intuitively show the multidimensional variation characteristics of the stress tensor. Summary of the Invention

[0003] Therefore, it is necessary for the present invention to provide a three-dimensional geostress modeling method, system and storage medium to solve at least one of the above-mentioned technical problems.

[0004] To achieve the above objectives, a three-dimensional geostress modeling method includes the following steps:

[0005] Step S1: Construct a three-dimensional spatial mesh model of the geological structure and label the spatial parameters of the distribution of strata, faults and fractures;

[0006] Step S2: Input the spatial parameters into the preset probabilistic mechanics model to generate a multi-scale fracture field;

[0007] Step S3: Embed the multi-scale fracture field into the initial tensor field of the preset three-dimensional stress field, and generate continuous stress distribution data containing the fracture propagation path through fracture evolution simulation;

[0008] Step S4: Construct a numerical update strategy for fluid-structure interaction based on continuous stress distribution data, couple the fluid pressure field and the solid stress field for solution, and update the three-dimensional stress field;

[0009] Step S5: Extract the stress tensor data from the updated three-dimensional stress field, map the stress tensor data to a three-dimensional spatial voxel mesh, and generate a three-dimensional geostress visualization model.

[0010] Preferably, the present invention also provides a three-dimensional geostress modeling system for performing the above-described three-dimensional geostress modeling method, the three-dimensional geostress modeling system comprising:

[0011] The geological spatial modeling module is used to construct a three-dimensional spatial mesh model of geological structures and to annotate the spatial parameters of strata, faults and fractures.

[0012] The probabilistic fracture generation module is used to input spatial parameters into a preset probabilistic mechanical model to generate a multi-scale fracture field;

[0013] The fracture evolution simulation module is used to embed multi-scale fracture fields into the initial tensor field of a preset three-dimensional stress field, and generate continuous stress distribution data containing fracture propagation paths through fracture evolution simulation.

[0014] The fluid-structure interaction update module is used to construct a numerical update strategy for fluid-structure interaction based on continuous stress distribution data, to couple and solve the fluid pressure field and the solid stress field, and to update the three-dimensional stress field.

[0015] The 3D stress visualization module is used to extract stress tensor data from the updated 3D stress field, map the stress tensor data to a 3D spatial voxel mesh, and generate a 3D geostress visualization model.

[0016] Preferably, the present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed, implements the three-dimensional geostress modeling method.

[0017] This invention comprehensively considers the coupling effects of multiple physical field factors such as stratigraphic structure, fault activity, fracture distribution, and fluid seepage within the same computational framework, thereby significantly improving the spatial resolution and temporal response accuracy of the geostress field in complex geological environments. By introducing probabilistic mechanics to generate a multi-scale fracture field, the spatial distribution and occurrence of fractures possess statistical randomness and geological rationality, avoiding stress field biases caused by traditional static fracture assumptions. By embedding a fracture network in a three-dimensional stress tensor field and tracking fracture propagation paths, the invention achieves a dynamic reflection of fracture evolution on local stress concentration, transfer, and attenuation, enabling the model to realistically reproduce the entire process of fracture initiation, propagation, and interconnection. Furthermore, by constructing a flow... The system couples the control equations and updates the interaction between pore pressure and solid stress at the time-step level, overcoming the response lag problem in previous unidirectional loading or weak coupling modes. It can synchronously reflect the stress redistribution effect caused by pore pressure changes. By introducing a frequency domain pre-filtering mechanism, numerical noise is effectively suppressed and the stability of the stress update process is enhanced. Finally, by mapping the multidimensional stress tensor to a regular voxel mesh and establishing the correspondence between stress scalars and color values, a continuous visualization of the stress field in three-dimensional space is realized. It can intuitively display the principal stress direction, shear concentration zone, and stress gradient distribution, providing an accurate and reliable stress field characterization method for geomechanical analysis, rock mass stability evaluation, and underground engineering design. Attached Figure Description

[0018] Other features, objects, and advantages of the invention will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings:

[0019] Figure 1 This is a schematic diagram of the steps of a three-dimensional geostress modeling method according to the present invention;

[0020] Figure 2 This is a frequency distribution diagram of the crack scale in this invention;

[0021] Figure 3 This is a schematic diagram of the three-dimensional geostress model of the present invention. Detailed Implementation

[0022] The technical method of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0023] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0024] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0025] To achieve the above objectives, please refer to Figures 1 to 3 This invention provides a three-dimensional geostress modeling method, the method comprising the following steps:

[0026] Step S1: Construct a three-dimensional spatial mesh model of the geological structure and label the spatial parameters of the distribution of strata, faults and fractures;

[0027] Step S2: Input the spatial parameters into the preset probabilistic mechanics model to generate a multi-scale fracture field;

[0028] Step S3: Embed the multi-scale fracture field into the initial tensor field of the preset three-dimensional stress field, and generate continuous stress distribution data containing the fracture propagation path through fracture evolution simulation;

[0029] Step S4: Construct a numerical update strategy for fluid-structure interaction based on continuous stress distribution data, couple the fluid pressure field and the solid stress field for solution, and update the three-dimensional stress field;

[0030] Step S5: Extract the stress tensor data from the updated three-dimensional stress field, map the stress tensor data to a three-dimensional spatial voxel mesh, and generate a three-dimensional geostress visualization model.

[0031] Preferably, step S1 includes:

[0032] Collect geological exploration data of the target area and perform preprocessing on the geological exploration data to extract spatial geometric information of stratigraphic interfaces, fault trajectories and fracture groups;

[0033] Unstructured spatial grids are determined based on spatial geometric information in order to construct three-dimensional spatial grid models of geological structures;

[0034] In a three-dimensional spatial grid model, the spatial geometric information of stratigraphic interfaces, fault trajectories, and fracture groups is mapped to the corresponding grid cells, and their spatial location, attitude, and attribute parameters are labeled to form spatial parameters.

[0035] In this embodiment, firstly, geological exploration operations are conducted in the target area using an electromagnetic seismic probe, ground-penetrating radar, and a high-precision borehole sounder. The seismic probe acquires seismic reflected wave signals at sampling intervals of 50 meters, with the reflected wave sampling frequency set to 1000 Hz, to identify the location of stratigraphic interface changes. The ground-penetrating radar uses a 200 MHz center frequency antenna to perform detailed scanning of shallow faults and fracture distribution, with a scanning depth range of 0–60 meters and a resolution of 0.2 meters. The borehole sounder records lithological changes and physical parameters at each depth along the longitudinal direction, including density, wave velocity, and porosity, with a sampling step size of 1 meter. After the above data acquisition is completed, the original geological exploration dataset of the target area is obtained.

[0036] When preprocessing the collected geological exploration data, noise removal is first performed on the seismic reflection wave data. A sliding window smoothing filter is used to smooth the original waveform signal, with a window width set to 5 sampling points, removing high-frequency random interference signals. Then, the depth of the stratigraphic interface is calculated based on the time history difference of the reflected waves using a time-depth conversion formula. This calculation determines the spatial coordinates of the stratigraphic interface. For ground-penetrating radar signals, envelope peak detection is used to extract reflection feature points of faults and fractures. Points with reflection intensities exceeding twice the average value are identified as potential fracture points, and fracture groups are delineated through spatial clustering. For borehole data, lithological abrupt changes are identified by comparing the lithological parameter variation trends at each measuring point, and the results are verified against the stratigraphic interface. After processing, a spatial geometric information set of stratigraphic interfaces, fault trajectories, and fracture groups is formed.

[0037] Based on the aforementioned spatial geometric information, an unstructured spatial mesh for the geological structure was established. Mesh generation employed a tetrahedral-based spatial discretization method. First, node sets were arranged at stratigraphic interfaces, with node spacing controlled within 10 meters. The mesh was then densified in areas containing faults and fractures, with the minimum element side length set to 2 meters to enhance the analytical accuracy of stress distribution. To ensure continuity between fault planes and fracture surfaces within the mesh, interpolation points were set along the fault plane normal direction, with an interpolation step size not exceeding 1 meter. A complete spatial mesh element set was generated by tetrahedral connections. Each element contains a unique number and node coordinates.

[0038] In the established 3D spatial grid, the spatial geometric information of stratigraphic interfaces, fault trajectories, and fracture groups is mapped to the corresponding grid cells. The mapping operation is implemented using a coordinate matching method: for stratigraphic interface points, the grid cell number is determined based on their 3D coordinates, and the cell is marked as a stratigraphic boundary cell; for fault trajectory points, the normal angle of the cell is calculated, and if the angle is less than 10°, the cell is marked as a fault cell; for fracture group points, based on their trace length and aperture, the fractures are projected onto the inner surface of the grid cell, and the fracture attitude angle (strike angle) is labeled. The range is 0° to 180°, and the tilt angle is... The range is 30° to 90°.

[0039] Each labeled unit has additional attribute parameters, including lithological density. (scope ), wave speed (Range 2500~6000 m / s), porosity (Range 0.05~0.25) and spatial location coordinates Simultaneously, the friction coefficient of the fault unit was recorded. (Range 0.3~0.7) and the interfacial spacing of fracture elements (Range 0.1–1.5 m). The labeled set of grid cells forms a complete spatial parameter dataset.

[0040] Preferably, step S2 includes:

[0041] Feature extraction is performed on the fault and fracture distribution data in the spatial parameters to generate a fracture feature dataset, which includes spatial density distribution cloud map, dominant orientation distribution map and fracture scale frequency distribution data.

[0042] The fracture feature dataset is input into a pre-defined probabilistic mechanical model, and the spatial location and orientation of the fracture are generated by the stochastic process unit in the pre-defined probabilistic mechanical model.

[0043] Based on the frequency distribution data of fracture scale, a corresponding scale attribute is assigned to each generated fracture to form a multi-scale fracture field.

[0044] In this embodiment, based on spatial parameter data, fault and fracture distribution data are first extracted. Fault data includes the three-dimensional coordinates of the fault plane center point. , directional angle (Range 0°~180°), Inclination (Range 30°~90°), fault plane length (Range 100–2000 m) and inter-surface spacing (Range 2–10 m). Fracture data includes the coordinates of the fracture center point. , directional angle ,inclination , trace length (Range 1–50 m), opening (Range 0.1–2.0 mm).

[0045] In the feature extraction stage, a three-dimensional mesh metaspace is first constructed to... Using a unit cell size, fault and fracture point data are mapped to the center of the volume element. The number of fracture segments within each volume element is then calculated. and according to the formula Determine the fracture density ,in For unit volume ( ).Will The distribution is projected onto a three-dimensional coordinate space to form a fracture spatial density distribution cloud map. Subsequently, based on the strike and dip angle data of the faults and fractures, an angle statistical method is used to count the number of fractures in each direction within a 360° range at 10° intervals, forming a dominant orientation distribution map. For the fracture scale parameter, the trace length of each fracture is statistically analyzed. With opening Based on the trace length intervals [1–5 m], [5–15 m], [15–30 m], and [30–50 m], the scale levels were divided, and the proportion of fractures in each level was calculated to obtain fracture scale frequency distribution data. This resulted in a fracture feature dataset, including spatial density distribution cloud maps, dominant orientation distribution maps, and fracture scale frequency distribution data.

[0046] In the stage of generating the spatial location and attitude of the fractures, the aforementioned fracture feature dataset is imported into the probabilistic mechanics calculation unit. This unit uses a random number generator as its core to generate the coordinates of the fracture center point in three-dimensional space. The seed value of the random number generator is set to a fixed integer (e.g., 20251020) to ensure consistent results in repeated calculations. The location of the fracture center point is determined as follows: within the boundary range of the target area (0–1000 m in the X direction, 0–1000 m in the Y direction, and 0–500 m in the Z direction), a set of uniformly distributed spatial coordinate sequences is generated. Based on the density of the fracture space Determine the number of fractures generated in each region, and the number of fractures. ,in This corresponds to the volume of the volume element.

[0047] For the attitude of each fracture, the orientation of the main fracture group is determined according to the direction interval with the highest frequency in the dominant orientation distribution map. The strike angle is taken as the center value of this interval, and the dip angle is randomly selected within ±10° of the principal dip angle according to a normal distribution. The geometric shape of the fracture is represented by a disk shape, and the radius of the disk is... Corresponding trace length Derived from fracture-scale frequency distribution data. Fracture plane orientation is determined by strike angle. With tilt angle Determine the plane normal vector. This is used to define the spatial orientation of the fracture. The center coordinates of each fracture are calculated. Normal vector and geometric radius This forms a preliminary fracture space dataset.

[0048] Based on the fracture scale frequency distribution data, a corresponding scale attribute is assigned to each generated fracture. The fracture scale attribute includes trace length. With opening During allocation, the scale level is randomly selected based on the proportion in the frequency distribution. For example, if the proportion of a certain scale level is 30%, then 30% of the generated scale levels will be randomly assigned the trace length and aperture parameters of that scale. Scale trace length Range 1–50 m, opening degree Range: 0.1–2.0 mm. To maintain spatial continuity, the distance between adjacent fracture centers must not be less than 0.8 times the fracture trace length; otherwise, the fracture center coordinates must be re-extracted.

[0049] Through the above process, the spatial location, attitude, and scale properties of all fractures have been determined. Each fracture is then... The data is recorded in a specific format and aggregated to form a multi-scale fracture field. This fracture field is represented by a three-dimensional array structure, where each cell stores the fracture geometry parameters and spatial location index.

[0050] Preferably, feature extraction of fault and fracture distribution data in spatial parameters includes:

[0051] Identify the spatial coordinates of fault trajectories and fracture groups in spatial parameters;

[0052] Calculate the spatial density values ​​of fault trajectories and fracture groups within a preset three-dimensional volume element, and generate a spatial density distribution cloud map;

[0053] Statistical analysis of fault trajectory and fracture group attitude data generates a dominant orientation distribution map;

[0054] Based on the trace length and aperture of the fracture group, the fracture scale is divided into levels, and the distribution frequency of fractures at each level is statistically analyzed to generate fracture scale frequency distribution data.

[0055] In this embodiment, based on spatial parameter data, the spatial coordinates of the fault trajectory and fracture group are first identified. The spatial coordinates of the fault trajectory are derived from ground-penetrating radar scans and seismic reflector fitting results. The fault reflector data is then converted into a three-dimensional coordinate point cloud format, with each point containing coordinates. The coordinate accuracy is 0.1 m. The spatial coordinates of the fracture group are derived from borehole imaging data and surface scan data, and the center point of each fracture is determined by point set matching. and endpoint coordinates , All coordinates are uniformly projected to the National Geodetic Coordinate System CGCS2000 to ensure consistent spatial correspondence. Coordinate data is stored in tabular form, with each record corresponding to a fault or fracture unit.

[0056] After obtaining the spatial coordinates of the faults and fractures, the target region is divided into regular three-dimensional voxels for calculating the spatial density. The voxels have a cubic structure, with a side length of 10 m and a volume of [missing information]. Iterate through the coordinates of each fault or fracture point and determine its volume element number. Calculate the number of fractures in each volume element. Then, according to the formula Calculate fracture density ,in Volumetric volume. Fault number density. The calculation method is the same. The density calculation results are in the form of... The density values ​​are mapped to three-dimensional spatial coordinates, and a linear interpolation method is used to transform the discrete density values ​​into a continuous spatial field. Color gradient intervals (low-density areas) are defined based on the density values. Blue, medium density area Green, high-density area (The above is in red), generating a spatial density distribution cloud map. Each volume element in the cloud map corresponds to a density value and a color identifier.

[0057] Subsequently, the attitude data of fault trajectories and fracture groups were statistically analyzed. The attitude of each fault and fracture was determined by its strike angle. and tilt angle Description. Strike angles are measured clockwise from 0° to 180° with north as the reference direction, and dip angles are measured from 0° to 90° with the horizontal plane as the reference. All strike angle data for fractures and faults are grouped at 10° intervals, and the number of fractures within each interval is counted. Calculate the proportion of this interval. The dip angle data were statistically analyzed at 5° intervals to create a dip angle frequency table. The strike and dip angle statistics were used to plot a dominant orientation distribution map. The distribution map is represented in polar coordinates, with the radial angle representing the strike and the radius representing the dip angle frequency. Directions with higher frequencies form prominent peaks, thus showing the main orientation characteristics of fractures and faults.

[0058] After completing the attitude statistics, the trace length and aperture data of the fracture swarm were extracted and scaled. Trace length Calculated from the coordinates of the fracture endpoints, Opening degree Based on the conversion of borehole imaging grayscale reflection distance, the grayscale conversion coefficient was set to 0.005 mm / grayscale level. The trace length of fractures was classified into four scale levels: Level I (1–5 m), Level II (5–15 m), Level III (15–30 m), and Level IV (30–50 m). The aperture of fractures was classified into three categories: Category A (0.1–0.5 mm), Category B (0.5–1.0 mm), and Category C (1.0–2.0 mm). For each fracture, its category was determined based on its trace length and aperture parameters. The number of fractures in each category was counted, and their distribution frequency was calculated. This generates frequency distribution data at the crack scale.

[0059] The spatial density distribution cloud map, dominant orientation distribution map, and fracture-scale frequency distribution data are stored in a unified data table. Each data entry includes the geographic coordinate range, density value, strike angle interval, dip angle interval, and corresponding fracture level frequency.

[0060] Preferably, step S3 includes:

[0061] The multi-scale fracture field is spatially superimposed with the initial tensor field of the preset three-dimensional stress field to form an initial stress field containing a fracture network.

[0062] In the initial stress field containing a fracture network, the initiation and propagation paths of cracks are traced using fracture evolution simulation methods;

[0063] Determine whether the crack tip in the initial stress field of the cracked network meets the crack initiation condition. When the crack initiation condition is met, simulate the crack propagation process and record the spatial coordinates of the crack tip to form a crack propagation path point set.

[0064] A continuous spatial trajectory is determined based on the crack propagation path point set, which serves as the crack propagation path.

[0065] The stress distribution is updated based on the crack propagation path, generating continuous stress distribution data that includes the crack propagation path.

[0066] In this embodiment, based on the multi-scale fracture field and the preset three-dimensional stress field initial tensor field, a spatial superposition operation is first performed. The initial tensor field is based on the three-dimensional stress tensor. The multi-scale fracture field is stored in the form of a grid, with each grid node having six independent component values ​​in MPa. The fracture center point coordinates are used to represent the fracture field. , trace length Opening degree Normal components To achieve superposition, the fracture field is projected onto the same three-dimensional grid coordinate system as the stress tensor field, and the grid cell number where the fracture is located is determined using a spatial matching method.

[0067] During the stacking process, the stress tensor of the mesh containing the crack needs to be corrected based on the presence of the crack. The normal stress of each crack element... With tangential stress Calculated using the following formula:

[0068] ;

[0069] ;

[0070] Will and Friction parameters with cracks and bond strength ( Range 0.3–0.7 Comparing within the range of 1–5 MPa, if If the mesh element is determined to be a stress discontinuity element, a stress correction factor is introduced into the tensor field. (Range 0.1–0.8). The corrected formula is as follows: ,in These are the corrected stress component values. The initial stress field, forming a fractured network, is calculated.

[0071] In the initial stress field containing a fracture network, the path of crack initiation and propagation is traced through fracture evolution. The local stress intensity factor is calculated starting from the tip of each fracture. , , The stress intensity factor is calculated through nodal stress interpolation:

[0072] , , ,in The crack half-length is defined as 0.1–1.0 m. With material fracture toughness Compare, The range of values ​​is ,when At that time, it was considered that the crack met the crack initiation conditions.

[0073] For the crack tip that meets the initiation condition, the crack propagation process is simulated. Propagation direction. Determined by the maximum tangential stress criterion, i.e., moving along the direction of maximum tangential stress. Expanding step size. Take 0.5 m as the reference depth, and update the coordinates of the crack tip with each expansion. Repeat the calculation up to the tips of all cracks. Until then, the coordinates of the tip of each expansion are recorded in the crack propagation path point set. This point set data includes the point number, spatial coordinates, and step size number.

[0074] A continuous spatial trajectory is determined based on the point set of the crack propagation path. Using the crack initiation and termination coordinates as the two endpoints, a continuous curve trajectory is calculated using 3D spline interpolation. The spline order is set to third order, and the interpolation interval is 0.1 m. The interpolation result is output as a crack propagation path data sequence, where each path consists of multiple continuous spatial points, and the path curvature does not exceed 30° / m to ensure spatial continuity.

[0075] The stress distribution is updated based on the crack propagation path. The stress in the elements along the crack path is redistributed, and the new stress values ​​are calculated using the stress equilibrium conditions on both sides of the crack surface.

[0076] ;

[0077] in This represents the average stress around the path. , The elastic modulus of the rock mass (range 20–60 GPa). For the local strain at the crack tip (range of values) After the update, the stress tensor values ​​of the elements near the crack path are rewritten into the tensor field to form continuous stress distribution data containing the crack propagation path.

[0078] Preferably, spatially superimposing the multi-scale fracture field with the initial tensor field of the preset three-dimensional stress field includes:

[0079] Extract the spatial location and geometric parameters of each fracture in the multi-scale fracture field;

[0080] The spatial location and geometric parameters of the crack are embedded as internal boundary conditions into the mesh model corresponding to the initial tensor field;

[0081] Set the stress transfer coefficient and construct the contact relationship of the stress discontinuity surface;

[0082] The initial stress field of the fractured network was calculated based on the contact relationship.

[0083] In this embodiment, the spatial location and geometric parameters of each fracture are first extracted from the multi-scale fracture field. The extracted data fields include the three-dimensional coordinates of the fracture centerline or center point. Coordinates of the two ends of the fracture, normal vector of the fracture plane , trace length (1–50 m), local half-width (0.1~1.0 m), opening (0.1~2.0 mm), coefficient of friction (0.3~0.7) and bond strength (1–5 MPa). All coordinates use a unified reference system (e.g., CGCS2000), with a coordinate accuracy of no less than 0.1 m. Geometric parameters are stored in tabular or array form, with array indices corresponding one-to-one with 3D mesh cell indices.

[0084] After extraction, the spatial location and geometric parameters of the cracks are embedded as internal boundary conditions into a mesh consistent with the preset initial tensor field of the three-dimensional stress field. The mesh uses the unstructured tetrahedral mesh or voxel mesh generated in step S1; the mesh element size is 10 m to 20 m in ordinary regions and 2 m to 5 m in crack-dense regions. The embedding operation is implemented as follows: traverse the coordinates of the endpoints and midpoints of each crack to determine its mesh element number; for cracks that pass through the element volume, the corresponding element is cut into two parts along the crack plane and an interface element is inserted at the cut surface; the interface element is a thin-shell element with a geometric thickness of... Set as fracture aperture Alternatively, a value of 0.001 m can be used to ensure numerical stability, with the interface element center normal aligning with the crack normal. If the crack is located in multiple adjacent elements, it is segmented element by element in the order it passes through the elements, and an interface element is inserted into each segmentation plane. The interface element number maintains a mapping relationship with the original element number to facilitate writing back the results. For irregular cracks, a triangular piece is used to approximate its plane, and the above segmentation and insertion operations are performed at each triangular piece.

[0085] Stress transfer coefficients are set on the interface elements to construct the contact relationship of the stress discontinuity surface. The contact relationship is represented by the interface normal stiffness. interface tangential stiffness The relationship between the two is described as a ratio. Regulation, The range is 0.1 to 1.0. Interface normal stiffness. The range of values ​​is (i.e., the unit stress caused by unit volume displacement), interface tangential stiffness Interfacial friction is determined by the coefficient of friction. With bond strength Control: Interface shear threshold ,in Let represent the normal stress at the interface. To describe the different force transmission capabilities of the interface under open and closed conditions, a normal opening factor is introduced. When the relative displacement in the normal direction jumps When (on), the normal stiffness is adjusted according to the coefficient. Downgraded to , Take 0.01 to 0.2; when (When closed) maintain Shear transfer in the sliding case is calculated using the slip correction factor. Treatment: When tangential shear stress At that time, the tangential stiffness is calculated according to Downgraded to Furthermore, the relative tangential displacement of the interface is treated according to the law of kinetic friction.

[0086] Based on the aforementioned contact relationship, the initial stress field forming the fractured network is solved using a global equilibrium process. The solution employs a linear elastic constitutive model. In each solid element, the constitutive matrix is... (Two-dimensional or three-dimensional isotropic elastic coefficient matrix) From elastic modulus (20–60 GPa) and Poisson's ratio (0.2~0.35) Determined: In the three-dimensional case, The contribution of interface elements is given by the following formula (in matrix form, denoted by symbols here). This contribution is expressed through the interface stiffness matrix. Superimposed on the global stiffness matrix That is, the global equilibrium equation is ,in , This is the global node displacement vector. The applied load vector (including gravity, boundary constraint reactions, and initial far-field stress converted to equivalent volumetric force terms) is used. The numerical solver employs either direct decomposition (LU decomposition) or iterative solution (conjugate gradient), with a solution tolerance set as... The maximum number of iterations is no more than 2000. The node displacements are obtained by solving the problem. Then, through the strain-displacement relationship Calculate the element strain, and then use the stress-strain relationship. Calculate the stress tensor of solid elements. The stress of interface elements is calculated according to the interface constitutive relation: normal direction. Tangential shear force (Pressing friction treatment in a sliding state).

[0087] To ensure the stability and continuity of the numerical implementation, several specific restrictions and processing rules are imposed on the embedding and solution process: the node matching of interface elements and solid elements adopts constraint coupling technology, and the interface element nodes and adjacent solid element nodes are connected through equivalent constraint matrices. Coupling; if element partitioning results in small-sized elements with a volume smaller than the minimum volume threshold. (For example If this is not the case, then neighboring small elements are merged, and interface elements are re-inserted after merging to avoid numerical singularities. The far-field principal stresses of the initial tensor field are applied at the boundary nodes with equivalent displacement constraints. The range of far-field principal stress values ​​is determined according to the formation depth; for example, at a depth of 200–500 m, the vertical stress is... Take 20–40 MPa, horizontal principal stress Take 12-30 MPa, secondary principal stress Take 8-20 MPa.

[0088] Preferably, step S4 includes:

[0089] Based on continuous stress distribution data, fluid-structure interaction control equations are constructed, and a numerical update strategy is developed. The fluid-structure interaction control equations include solving the fluid pressure field control equations and the solid stress field control equations.

[0090] Solve the governing equations of the fluid pressure field to obtain the pore pressure distribution at the current time step;

[0091] Substituting the pore pressure distribution as a volume force term into the governing equation of the solid stress field, the updated solid stress field is obtained by solving the equation.

[0092] Iterate until the convergence condition is met to complete the coupled solution;

[0093] During the coupled solution process, physical pre-filtering based on frequency domain analysis is applied to the input physical quantities;

[0094] The updated solid stress field is superimposed on the initial stress field;

[0095] The updated three-dimensional stress field is determined based on the superposition results, and the data of each component of the stress tensor are stored.

[0096] In this embodiment, based on continuous stress distribution data including crack propagation paths, fluid-structure interaction control equations are first constructed. These control equations consist of both fluid pressure field control equations and solid stress field control equations. For the fluid pressure field portion, the seepage control equation is used to describe the flow behavior of the porous fluid in the porous medium, and its expression is:

[0097] ;

[0098] in, fluid density (take) ), Porosity (taken as 0.05–0.25). For the fluid velocity vector, Source and sink items (unit: Fluid velocity Calculated according to Darcy's Law:

[0099] ;

[0100] In the formula, Permeability (take) ), For fluid viscosity (take) ), This represents the pore pressure (in Pa). Combining mass conservation and boundary conditions, the governing equations of the fluid pressure field can be discretized into time-step equations. Transient form:

[0101] ;

[0102] in The specific storage coefficient (taken as) ), This is the time step number.

[0103] The governing equations for the stress field in the solid are based on the three-dimensional linear elasticity equilibrium equations:

[0104] ;

[0105] In the formula For stress tensor, The density of the rock mass (taken as...) ), The gravitational acceleration vector The relationship between stress and strain is determined by the constitutive relation. It means that, among them Let be the material stiffness matrix, and take the elastic modulus. = 20~60 GPa, Poisson's ratio Isotropic form with a strain of 0.2–0.35. With displacement The gradient satisfies .

[0106] After establishing the aforementioned governing equations, a numerical update strategy for fluid-structure interaction is constructed. Iterative calculations are performed using a time-stepping method: at the... First, the governing equations of the fluid pressure field are solved to obtain the current pore pressure distribution. Then, the pore pressure is added as a volume force term to the solid stress field equation. Specifically, the additional stress is introduced from the pore pressure into the stress equilibrium equation. ,in The effective stress coefficient (taken as 0.8 to 1.0). The unit tensor is used. The updated solid equilibrium equations are:

[0107] ;

[0108] Discretize the equations using the finite volume method or the finite element method to form a global linear system of equations. ,in Here is the stiffness matrix of the solid. This represents the volume force term caused by pore pressure. The system of equations is solved using an iterative approach, updating the pore pressure and solid displacement in each iteration until the convergence condition is met.

[0109] The convergence condition is defined as follows: in two consecutive iterations, the changes in the pore pressure field and stress field satisfy...

[0110] and .

[0111] When both of the above conditions are met, the fluid-structure interaction calculation is considered to have converged.

[0112] During the iteration process, to eliminate the influence of high-frequency noise in the input physical quantities on coupling stability, the input variables (stress tensor, pore pressure, displacement field) undergo physical pre-filtering based on frequency domain analysis. The filtering steps include: performing a Fast Fourier Transform (FFT) on the time series data of each grid node, within the frequency range of 0~ The effective amount is retained. The frequency is set to 10 Hz, and components above this frequency are treated as 0; then, the time-domain signal is recovered through inverse IFFT. After filtering, a smooth sequence of input physical quantities is obtained. This process ensures that the time series of stress, displacement, and pressure are continuous and stationary.

[0113] After the coupling calculation is completed, the updated solid stress field is superimposed on the initial stress field. The superposition operation is performed point-by-point on the tensor components:

[0114] ;

[0115] in The initial stress field tensor obtained in step S3, This updates the tensor for the current step. After stacking, the stress components at each node are calculated. , , , , , .

[0116] Based on the superposition results, the updated 3D stress field is determined and the data of each component is stored. The storage format is a tensor component array corresponding to node indices, with each node recording six component values. The data precision is represented as double-precision floating-point numbers in MPa. A unified stress field data file containing node coordinates is established for the visualization mapping operation in subsequent step S5. Stress components and pore pressure values .

[0117] Preferably, step S5 includes:

[0118] Extract the stress tensor data from the updated three-dimensional stress field and normalize the stress tensor data.

[0119] Establish a regular voxel mesh corresponding to the three-dimensional spatial mesh model of the geological structure;

[0120] The normalized stress tensor data is mapped to the center point of the voxel mesh to generate mapped voxel mesh data.

[0121] Spatial consistency verification is performed on the mapped voxel grid data;

[0122] Based on the mapped voxel mesh data, the correspondence between stress scalars and color values ​​is constructed, and a three-dimensional geostress visualization model is generated.

[0123] In this embodiment, firstly, stress tensor data is extracted from the three-dimensional stress field. The three-dimensional stress field stores six independent components at the node level, including the normal stress component. , , and shear stress components , , Each component is measured in MPa, and the data format is double-precision floating-point. The extraction operation reads each node sequentially according to its index, extracting the three-dimensional coordinates of each node. Write the corresponding six stress component values ​​into the stress tensor matrix. ,in Number the nodes. Number the stress components. After extraction, a three-dimensional array containing all nodal stress tensors is obtained.

[0124] To eliminate numerical imbalances caused by differences in stress component amplitudes across different regions, the stress tensor data needs to be normalized. Normalization employs a linear scaling method. For each component... Calculate the maximum value of the entire field. and minimum value Then follow the formula

[0125] ;

[0126] All stress components are normalized to the 0–1 range. The relative proportions of the six components are preserved during normalization. To prevent numerical truncation errors, when… When, take the denominator as After processing, the normalized stress tensor matrix is ​​obtained. Each of its components is dimensionless and has a fixed range.

[0127] After normalization, a regular voxel mesh corresponding to the three-dimensional spatial mesh model of the geological structure is established. The voxel mesh adopts a cubic voxel structure, with each voxel having a side length of 1 m. The overall spatial extent of the voxel mesh is consistent with the outer boundary of the geological structure, and the boundary extent is determined by the spatial mesh model in step S1. For example, when the geological structure has an extent of 0–1000 m in the X direction, 0–800 m in the Y direction, and 0–500 m in the Z direction, the voxel mesh size is 1000×800×500 elements. Each voxel has a unique index number. The index order increases in the X, Y, and Z directions.

[0128] The normalized stress tensor data is mapped to the center points of the voxel mesh. This mapping operation is achieved through spatial coordinate mapping: for each voxel center point coordinate... Find the coordinates of the node whose spatial location is closest to it in the stress tensor matrix. The distance metric used is Euclidean distance. .like If the average mesh element size is less than 1 m, then the normalized stress tensor components corresponding to the nodes are assigned to the voxel centers; if For values ​​greater than 1 m, interpolation is performed using the inverse distance weighting method, with the weights... The weighted average yields the voxel central stress tensor components:

[0129] ;

[0130] Number of nodes within the interpolation range Eight voxels are selected to ensure spatial continuity. After mapping, a mapped voxel mesh dataset is generated, where each voxel contains six components of the normalized stress tensor and their spatial indices.

[0131] Spatial consistency verification is performed on the mapped voxel mesh data. The verification process is based on the stress gradient between adjacent voxels. For each pair of adjacent voxels... , Calculate the difference in stress components ,like If it is, then it is determined to be an inconsistent unit, where Take 0.05 (normalized units). If the average stress gradient of adjacent voxels in three directions... If the value exceeds a threshold of 0.03, local smoothing is performed in that region. The smoothing process uses a neighborhood averaging method: the target voxel stress value is replaced with a weighted average of its 26 neighboring voxels, with weights... ,in This represents the grid distance from neighboring voxels to the target voxel. After verification and smoothing, the entire voxel data field is ensured to be spatially continuous and consistent.

[0132] Based on the verified mapped voxel mesh data, a correspondence between stress scalars and color values ​​is established. First, the principal stress scalar value for each voxel is calculated. The von Mises equivalent stress formula is used:

[0133] ;

[0134] Will After normalization to the 0-1 range, a color mapping is established. The color correspondence uses a five-level gradient color band: blue corresponds to... Cyan corresponds Green corresponds to Yellow corresponds to Red corresponds to Color values ​​are stored in RGB numerical format, such as blue (0,0,255) and red (255,0,0). A value is calculated for each voxel. The corresponding RGB values ​​are obtained by looking up the table, and a voxel dataset containing coordinates, stress scalars and color ternary information is generated.

[0135] The color mapping results are input into a 3D visualization rendering system. The rendering system uses the center point of each voxel as a node and plots the 3D stress field distribution according to coordinates and color values. The plotting order is from bottom to top in the Z direction to avoid occlusion. The rendering resolution is one voxel per meter in each direction, and the image output format is 24-bit color depth. The final generated 3D geostress visualization model fully displays the stress distribution of each voxel in space, with continuous color gradients and accurate spatial correspondences.

[0136] Therefore, the embodiments should be considered as exemplary and non-limiting in all respects, and the scope of the invention is not limited by the foregoing description. Thus, all changes falling within the meaning and scope of the equivalents of the application are intended to be included within the scope of the invention.

[0137] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A three-dimensional geostress modeling method, characterized in that, Includes the following steps: Step S1: Construct a three-dimensional spatial mesh model of the geological structure and label the spatial parameters of the distribution of strata, faults and fractures; Step S2: Input the spatial parameters into the preset probabilistic mechanics model to generate a multi-scale fracture field; Step S3: Embed the multi-scale fracture field into the initial tensor field of the preset three-dimensional stress field, and generate continuous stress distribution data containing the fracture propagation path through fracture evolution simulation; Step S4: Construct a numerical update strategy for fluid-structure interaction based on continuous stress distribution data, couple the fluid pressure field and the solid stress field for solution, and update the three-dimensional stress field; Step S5: Extract the stress tensor data from the updated three-dimensional stress field, map the stress tensor data to a three-dimensional spatial voxel mesh, and generate a three-dimensional geostress visualization model.

2. The three-dimensional geostress modeling method according to claim 1, characterized in that, Step S1 includes: Collect geological exploration data of the target area and perform preprocessing on the geological exploration data to extract spatial geometric information of stratigraphic interfaces, fault trajectories and fracture groups; Unstructured spatial grids are determined based on spatial geometric information in order to construct three-dimensional spatial grid models of geological structures; In a three-dimensional spatial grid model, the spatial geometric information of stratigraphic interfaces, fault trajectories, and fracture groups is mapped to the corresponding grid cells, and their spatial location, attitude, and attribute parameters are labeled to form spatial parameters.

3. The three-dimensional geostress modeling method according to claim 1, characterized in that, Step S2 includes: Feature extraction is performed on the fault and fracture distribution data in the spatial parameters to generate a fracture feature dataset, which includes spatial density distribution cloud map, dominant orientation distribution map and fracture scale frequency distribution data. The fracture feature dataset is input into a pre-defined probabilistic mechanical model, and the spatial location and orientation of the fracture are generated by the stochastic process unit in the pre-defined probabilistic mechanical model. Based on the frequency distribution data of fracture scale, a corresponding scale attribute is assigned to each generated fracture to form a multi-scale fracture field.

4. The three-dimensional geostress modeling method according to claim 3, characterized in that, Feature extraction of fault and fracture distribution data in spatial parameters includes: Identify the spatial coordinates of fault trajectories and fracture groups in spatial parameters; Calculate the spatial density values ​​of fault trajectories and fracture groups within a preset three-dimensional volume element, and generate a spatial density distribution cloud map; Statistical analysis of fault trajectory and fracture group attitude data generates a dominant orientation distribution map; Based on the trace length and aperture of the fracture group, the fracture scale is divided into levels, and the distribution frequency of fractures at each level is statistically analyzed to generate fracture scale frequency distribution data.

5. The three-dimensional geostress modeling method according to claim 1, characterized in that, Step S3 includes: The multi-scale fracture field is spatially superimposed with the initial tensor field of the preset three-dimensional stress field to form an initial stress field containing a fracture network. In the initial stress field containing a fracture network, the initiation and propagation paths of cracks are traced using fracture evolution simulation methods; Determine whether the crack tip in the initial stress field of the cracked network meets the crack initiation condition. When the crack initiation condition is met, simulate the crack propagation process and record the spatial coordinates of the crack tip to form a crack propagation path point set. A continuous spatial trajectory is determined based on the crack propagation path point set, which serves as the crack propagation path. The stress distribution is updated based on the crack propagation path, generating continuous stress distribution data that includes the crack propagation path.

6. The three-dimensional geostress modeling method according to claim 5, characterized in that, Spatial superposition of the multi-scale fracture field with the initial tensor field of the preset three-dimensional stress field includes: Extract the spatial location and geometric parameters of each fracture in the multi-scale fracture field; The spatial location and geometric parameters of the crack are embedded as internal boundary conditions into the mesh model corresponding to the initial tensor field; Set the stress transfer coefficient and construct the contact relationship of the stress discontinuity surface; The initial stress field of the fractured network was calculated based on the contact relationship.

7. The three-dimensional geostress modeling method according to claim 1, characterized in that, Step S4 includes: Based on continuous stress distribution data, fluid-structure interaction control equations are constructed, and a numerical update strategy is developed. The fluid-structure interaction control equations include solving the fluid pressure field control equations and the solid stress field control equations. Solve the governing equations of the fluid pressure field to obtain the pore pressure distribution at the current time step; Substituting the pore pressure distribution as a volume force term into the governing equation of the solid stress field, the updated solid stress field is obtained by solving the equation. Iterate until the convergence condition is met to complete the coupled solution; During the coupled solution process, physical pre-filtering based on frequency domain analysis is applied to the input physical quantities; The updated solid stress field is superimposed on the initial stress field; The updated three-dimensional stress field is determined based on the superposition results, and the data of each component of the stress tensor are stored.

8. The three-dimensional geostress modeling method according to claim 1, characterized in that, Step S5 includes: Extract the stress tensor data from the updated three-dimensional stress field and normalize the stress tensor data. Establish a regular voxel mesh corresponding to the three-dimensional spatial mesh model of the geological structure; The normalized stress tensor data is mapped to the center point of the voxel mesh to generate mapped voxel mesh data. Spatial consistency verification is performed on the mapped voxel grid data; Based on the mapped voxel mesh data, the correspondence between stress scalars and color values ​​is constructed, and a three-dimensional geostress visualization model is generated.

9. A three-dimensional geostress modeling system, characterized in that, For performing the three-dimensional geostress modeling method as described in claim 1, the three-dimensional geostress modeling system comprises: The geological spatial modeling module is used to construct a three-dimensional spatial mesh model of geological structures and to annotate the spatial parameters of strata, faults and fractures. The probabilistic fracture generation module is used to input spatial parameters into a preset probabilistic mechanical model to generate a multi-scale fracture field; The fracture evolution simulation module is used to embed multi-scale fracture fields into the initial tensor field of a preset three-dimensional stress field, and generate continuous stress distribution data containing fracture propagation paths through fracture evolution simulation. The fluid-structure interaction update module is used to construct a numerical update strategy for fluid-structure interaction based on continuous stress distribution data, to couple and solve the fluid pressure field and the solid stress field, and to update the three-dimensional stress field. The 3D stress visualization module is used to extract stress tensor data from the updated 3D stress field, map the stress tensor data to a 3D spatial voxel mesh, and generate a 3D geostress visualization model.

10. A computer-readable storage medium, characterized in that, It stores a computer program, characterized in that, when the computer program is executed, it implements the three-dimensional geostress modeling method as described in any one of claims 1-8.

Images