A conversion mapping method from a three-dimensional geological grid model to a numerical analysis model

Abstract

The present application relates to a kind of three-dimensional geological grid model to numerical analysis model conversion and attribute data mapping method, to improve the data conversion accuracy and efficiency between geological modeling and numerical analysis.The method includes the following steps: first, from the geological modeling software export includes grid geometry information and physical attribute first grid data file;By depth analysis, topological rearrangement and unit type matching and the method for accurately mapping physical attribute data generation second format file in line with the requirement of target software;Finally import second format file in numerical analysis software, automatically generate numerical analysis model has been endowed with initial attribute, complete from geological model to numerical analysis model seamless conversion.The present application is accurately mapped to the attribute data of target grid topological sequence, strictly aligns physical attribute data and geometric data, solves the attribute loss problem caused by topological rearrangement and geometry conversion in traditional method, so as to improve the accuracy and reliability of numerical analysis model.

Description

Technical Field

[0001] This invention relates to the fields of geomechanics and oil and gas exploration, and in particular to a method for converting and mapping a three-dimensional geological grid model to a numerical analysis model. Background Technology

[0002] In fields such as geotechnical engineering, geomechanics, and oil and gas exploration, numerical simulation has become an indispensable analytical tool. The typical workflow consists of two stages: First, a three-dimensional geological model reflecting the real stratigraphic structure is constructed using professional geological modeling software. This model includes detailed mesh geometry and key physical properties (such as lithology, porosity, elastic parameters, etc.). Second, the geological model is imported into numerical analysis software for mechanical, fluid, or coupled process simulations.

[0003] However, fundamental differences in the underlying grid data structures, file formats, and definition methods used by geological modeling software and numerical analysis software create significant obstacles to model conversion. Currently, commonly used conversion methods suffer from the following drawbacks: 1. Manual Conversion: This relies on engineers manually reconstructing the geometric mesh in the numerical analysis software based on the geological model's outline. This process is extremely tedious, time-consuming, and struggles to guarantee accurate geometric reproduction. 2. Loss of Attribute Information: Even when geometric mesh conversion is achieved using neutral formats (such as .STL, .DXF), the rich attribute parameters in the geological model cannot be transferred. These parameters must be manually reassigned in the numerical analysis software, which is not only labor-intensive but also highly susceptible to human error. 3. Workflow Interruption: The entire modeling workflow cannot be automated, severely restricting the efficiency of numerical analysis and the repeatability of the model. This problem is particularly prominent in scenarios requiring repeated calculations of multiple schemes and parameters.

[0004] Therefore, there is an urgent need in this field for an automated, high-fidelity conversion method from geological modeling grids to numerical analysis grids, especially a technical solution that can fully inherit geological attribute parameters. Summary of the Invention

[0005] To address the aforementioned problems, this invention provides a method for converting and mapping a three-dimensional geological grid model to a numerical analysis model, thereby solving the problems of cumbersome manual conversion from geological modeling software to numerical analysis software and inaccurate mapping between physical properties and geometric data in existing technologies.

[0006] The technical solution of this invention is:

[0007] A method for converting and mapping a three-dimensional geological grid model to a numerical analysis model includes the following steps:

[0008] S1. Export a first grid data file containing grid geometry information and at least one physical property parameter from geological modeling software. The first grid data file is in .GRDECL format.

[0009] S2. Perform deep analysis on the first grid data file to extract the column coordinates and corner coordinates, and obtain the original dataset of the grid geometry.

[0010] S3. Using the original dataset, the corner grid data is converted into hexahedral cell grid data by topology rearrangement and cell type matching methods, and the target grid topology sequence is obtained.

[0011] S4. Accurately map the attribute data of the target mesh topology sequence, precisely align the physical attribute data with the geometric data, and generate a second format file, the second format file being in the .f3dat format;

[0012] S5. Import the second format file into the target numerical analysis software to automatically generate a numerical analysis calculation grid model containing complete physical properties and output the calculation results for numerical simulation and analysis.

[0013] Furthermore, the specific steps for deep parsing the first grid data file include:

[0014] S21. Read the COORD array from the first grid data file, divide the data in the COORD array into groups of 6 floating-point numbers, and extract the three-dimensional coordinates of the top and bottom of the column one by one;

[0015] S22. By traversing the grid dimensions, construct the parameterized equation for each bar based on the bar coordinates in the COORD array, which is used to describe the trajectory of the bar in three-dimensional space.

[0016] S23. Read the ZCORN array from the first grid file, and extract the depth values ​​of the 8 corner points of each grid cell from the ZCORN array using an indexing algorithm;

[0017] S24. Read physical property data from the first grid data file and align this property data with the grid data;

[0018] S25. Using the cylindrical trajectory in the COORD array, calculate the three-dimensional coordinates of each corner point through linear interpolation.

[0019] S26. Obtain the original dataset including column coordinates, corner coordinates, and physical property data.

[0020] Furthermore, the specific steps for topology rearrangement and cell type matching are as follows:

[0021] S31. Read the ACTNUM array, and based on the array's labels, remove invalid cells, including those with zero volume and those with overlapping corners, and filter out valid cells.

[0022] S32. Based on the corner coordinates obtained in S2, rearrange the corner order of each effective unit using the corner order replacement method;

[0023] S33. According to the right-handed spiral order specific to the numerical analysis software, the corner index is replaced to ensure that the corner order of each grid cell conforms to the topology requirements of the target software.

[0024] S34. Generate a valid cell ID mapping table based on the correspondence between the original mesh and the target mesh;

[0025] S35. Based on the effective cell ID mapping table and the grid data after corner point order adjustment, and combined with the topology format required by the numerical analysis software, construct the target grid topology sequence that meets the requirements of the target numerical analysis software using the mapping table and corner point order arrangement results.

[0026] Furthermore, the specific steps for accurately mapping the attribute data of the target mesh topology sequence are as follows:

[0027] S41. Initialize the counter ZoneID, and based on the valid cell IDs in the target mesh topology sequence, perform a geometric information traversal of the spatial cells and calculate the attribute data index of each mesh cell.

[0028] S43. Verify the validity of each grid cell by reading the ACTNUM array;

[0029] S44. During the mapping process, a corresponding FLAC3D command string is generated for each valid mesh element;

[0030] S45. Generate a second format file whose physical properties perfectly match the geometric information of each mesh cell.

[0031] Furthermore, the formula for calculating the attribute data index of each grid cell in S41 is as follows: ,

[0032] Where i, j, k are the indices of the grid cells in the X, Y, Z directions; NX, NY, NZ are the dimensions of the grid.

[0033] Furthermore, the method for determining whether each grid cell is valid by reading the ACTNUM array is as follows:

[0034] When ACTNUM[i,j,k]=1, the cell is a valid cell and can be used for attribute mapping;

[0035] When ACTNUM[i,j,k]=0 and other tags are used, the cell is an invalid cell and is removed.

[0036] Furthermore, the method is accomplished using a conversion script implemented in the Python programming language.

[0037] Furthermore, the physical property parameters in S1 include one or more of pore pressure, porosity, permeability, Young's modulus, and Poisson's ratio.

[0038] Furthermore, the geological modeling software includes Petrel, Roxar RMS, SKUA-GOCAD, and ECLIPSE.

[0039] Furthermore, the numerical analysis software includes FLAC3D, 3DEC, and PFC3D.

[0040] The beneficial effects of this invention are:

[0041] 1. Improve data conversion accuracy

[0042] This invention precisely maps the attribute data of the target mesh topology sequence in step S4, strictly aligning the physical attribute data with the geometric data. This solves the attribute loss problem caused by topological rearrangement and geometric transformation in traditional methods, thereby improving the accuracy and reliability of the numerical analysis model.

[0043] 2. Avoid human error

[0044] The automated script conversion process of steps S1 to S5 of this invention avoids the tediousness and errors of manual operation in traditional methods, introduces a direct writing method based on index alignment, reduces human error in data conversion, and ensures data consistency and accuracy. Attached Figure Description

[0045] Figure 1 This is a flowchart of the method of the present invention;

[0046] Figure 2 It is the original porosity model within the petrel;

[0047] Figure 3 It is the porosity model numerically simulated by FLAC3D software after conversion. Detailed Implementation

[0048] The following will combine Figures 1-3The technical solutions in the embodiments of the present invention will be clearly and completely described herein. The described embodiments are merely some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments in this application without creative effort are within the scope of protection of this application.

[0049] Example 1:

[0050] Geological modeling software is typically used to construct detailed 3D geological mesh models, while numerical analysis software is used to simulate mechanical, fluid, or coupled processes based on these models. However, due to fundamental differences in mesh data structure, file format, and definition methods, the conversion process from geological modeling software to numerical analysis software suffers from problems such as loss of accuracy and incorrect attribute mapping. Traditional conversion methods usually rely on manual operation, which is not only cumbersome but also prone to human error, severely limiting the efficiency and accuracy of numerical simulations.

[0051] Therefore, this application provides a method for converting and mapping a three-dimensional geological grid model to a numerical analysis model, including the following steps:

[0052] S1. Export a first grid data file containing grid geometry information and at least one physical property parameter from geological modeling software. The first grid data file is in .GRDECL format.

[0053] Specifically, this implementation uses Petrel as the geological modeling software. First, the software is launched, and the previously created 3D geological model is opened and loaded, such as... Figure 2As shown. In Petrel, export the mesh data through the "Export" option under the "File" menu. In the export settings, select .GRDECL as the export file format. .GRDECL is a standard mesh data format suitable for storing mesh geometry and attribute data. In the export settings, check the mesh's geometric information and physical properties. Geometric information includes the COORD array (cylindrical coordinates) and ZCORN array (corner coordinates), and physical properties include porosity, permeability, Young's modulus, etc. Determine the export area; this could be the entire geological model or only a specific layer or region. You can export the entire mesh or select only data within a specific depth range, depending on your needs. After confirming all settings are correct, click the "Export" button to start generating the .GRDECL format file. This file will contain: column coordinate data (COORD array) describing the top and bottom coordinates of each column; corner coordinate data (ZCORN array) containing depth information for the eight corner points of each grid cell, which will be interpolated to obtain the three-dimensional coordinates later; and physical property data such as porosity, permeability, and Young's modulus, with corresponding attribute values ​​for each grid cell. After exporting the file, save the generated .GRDECL file to the specified directory to ensure that this file can be used for subsequent mesh analysis and attribute mapping.

[0054] S2. Deeply analyze the first grid data file to extract the column coordinates, corner coordinates, and physical property data to obtain the original dataset of the grid geometry.

[0055] Furthermore, the specific steps for deep parsing the first grid data file include:

[0056] S21. Read the COORD array from the first grid data file, divide the data in the COORD array into groups of 6 floating-point numbers, and extract the three-dimensional coordinates of the top and bottom of the column one by one;

[0057] Specifically, this process is automated by a conversion script written in Python. The first step involves reading and parsing the .GRDECL file and performing efficient data processing using the NumPy library. First, the specified .GRDECL file is opened, and data is read and parsed from it. The data in the file is stored in text format and contains different parts of data (such as grid dimensions, COORD data, ZCORN data, etc.). Then, data containing multiple floating-point numbers is read and stored row-wise in a two-dimensional array. More specifically, the COORD array portion of the file is read to obtain the coordinate data describing the grid cylinders. The script first parses the file header information to extract the grid dimensions (NX, NY, NZ), which further helps determine the grid size and the distribution of individual cells. Subsequently, the script locates the main data block COORD storing the cylinder coordinates, which contains the top and bottom coordinates of each cylinder. During the reading process, firstly, Python's string manipulation methods are used to replace 'D' with 'E' in scientific notation. Then, the numpy.array function is used to convert the cleaned list of floating-point numbers into a high-precision NumPy floating-point array, and immediately the numpy.reshape function is called to reshape its dimensions to ((N... x +1)×(N y The array is in matrix format (+1), 6), where each row corresponds precisely to the top and bottom coordinates of a pillar. Each row in this array contains 6 floating-point numbers, representing the top coordinate (X, Y, 6) of each pillar. top ,Y top Z top ) and bottom coordinates (X) bottom ,Y bottom Z bottom Each row in the COORD array represents a bar, containing six floating-point numbers. The script splits the data into groups of six floating-point numbers to ensure that each group represents the coordinates of a complete bar. Furthermore, the first three floating-point numbers (X...) top ,Y top Z top ) represents the top coordinate of the bar, and the last three floating-point numbers (X) represent the top coordinates of the bar. bottom ,Y bottom Z bottomThe top and bottom coordinates of each bar are represented by the number of floating-point numbers (N). Using NumPy array slicing, the script extracts every six floating-point numbers sequentially and assigns them to the top and bottom coordinates. The top and bottom coordinates of each bar are stored in two separate arrays for subsequent mesh parameterization. The extracted bar coordinates are stored in NumPy arrays, which will serve as input for subsequent steps, including bar parameterization equation construction, mesh rearrangement, and attribute mapping. The top and bottom coordinates of each bar are neatly stored in arrays to ensure data consistency and facilitate subsequent operations. It's worth noting that after reading the data stream from the COORD array, the script immediately performs a dimension consistency check. The system parses the mesh dimensions (N) based on the SPECGRID keyword. x N y N z ), calculate the theoretically required total number of coordinate data (Count). expected =(N x +1)×(N y +1)×6. The script will multiply the total number of floating-point numbers actually read by Count. expected The two are compared. If they are not equal, it indicates that the source file may have data truncation or format corruption. The script will throw an exception signal of 'COORD data size mismatch' and terminate the conversion process to prevent erroneous data from entering subsequent calculation stages.

[0058] S22. By traversing the grid dimensions, construct the parameterized equation for each bar based on the bar coordinates in the COORD array, which is used to describe the trajectory of the bar in three-dimensional space.

[0059] Specifically, based on the top and bottom coordinates of the COORD array extracted and stored in the preceding steps, each dimension (X, Y, Z) of the grid is traversed using nested loops. For each grid cell, the script iterates through the corresponding bar. Each bar has a top coordinate (X...). top ,Y top Z top ) and bottom coordinates (X) bottom ,Y bottom Z bottom These are the two endpoints of the cylinder in three-dimensional space. For each cylinder, the script constructs a parametric equation for the cylinder based on its top and bottom coordinates. Assuming the cylinder extends along the z-axis, the form of the equation describes the change in the cylinder's coordinates from top to bottom. For each z-value of each cylinder (from Z... top To Z bottomBy calculating the corresponding X and Y coordinates using the parametric equations described above (within the interval between these values), the points at different heights (Z values) of the cylinder in three-dimensional space can be obtained, thus plotting the cylinder's trajectory. All calculated cylinder trajectory points (X, Y, Z) are stored in a data structure for subsequent mesh topology rearrangement and attribute mapping operations. This trajectory data provides a detailed description of the cylinder's geometry, ensuring the accuracy of the mesh's geometric model and subsequent analysis.

[0060] Furthermore, let the coordinate of the top of a certain column be P. top (x t ,y t ,z t The bottom coordinate is P. bottom (x b ,y b ,z b ).

[0061] For a given corner depth z, its corresponding planar coordinates P(x,y) are obtained through the following linear interpolation formula (i.e., the discrete solution of the parameterized equation), which can be expressed as:

[0062]

[0063]

[0064] Where t is the normalized interpolation coefficient; x is the x-coordinate of the current corner point; x t x is the x-coordinate of the center of the top of the column; b The x-coordinate of the bottom center of the column; the y-coordinate of the current corner point; y t The y-coordinate of the center of the top of the column; y b The y-coordinate is the center of the bottom of the column. It should also be noted that, to prevent geometric distortion, the t value was truncated (e.g., limited to the range [-0.5, 1.5]) to handle extremely distorted meshes.

[0065] S23. Read the ZCORN array from the first grid file, and extract the depth values ​​of the 8 corner points of each grid cell from the ZCORN array using an indexing algorithm;

[0066] Specifically, the ZCORN array stores the corner depth information of grid cells. The file content is in text format, with each data item separated by spaces or other delimiters, and the data arranged line by line. The script locates the ZCORN array section in the file. The script ensures correct location of the ZCORN array by reading the file header or other marker information. Each line of data in this array corresponds to the 8 corner depths of one grid cell. Each grid cell has 8 corners, and each corner has a depth value. The ZCORN array in the file is read using NumPy and converted into a NumPy array. Each line contains 8 floating-point numbers, representing the depth values ​​of the 8 corners of one grid cell. The data for each grid cell is stored sequentially, containing the depth values ​​of the 8 corners. An indexing algorithm is used to extract the depth values ​​of the 8 corners corresponding to each grid cell from the ZCORN array. The script iterates through each line of the ZCORN array, extracting the 8 floating-point numbers from each line into an array or list, each representing the depth value of the 8 corners of that grid cell. The extracted eight corner depth values ​​will be stored in a data structure, typically a two-dimensional array or list, ensuring that the depth data for each grid cell is accurately recorded. The corner depth data for each grid cell is stored along with its index (such as the grid cell ID or other identifier) ​​for subsequent topology reordering and attribute mapping operations. It's worth noting that, similar to COORD processing, after extracting the ZCORN data, the script performs a total count check. The total theoretical depth points (Count) are calculated. zcorn =N x ×N y ×N z ×8. If the amount of extracted data does not meet the formula, it is determined that the grid corner data is missing. The system will automatically stop processing and prompt an error to ensure that each grid cell can be assigned a complete set of 8 depth values.

[0067] Furthermore, the specific execution steps of the indexing algorithm are as follows: traverse each row in the ZCORN array; each row contains 8 floating-point numbers, representing the depth values ​​of the 8 corner points of the grid cell; store the extracted 8 depth values ​​in an appropriate data structure (such as a list, array, or dictionary) for later use.

[0068] S24. Read physical property data from the first grid data file and align this property data with the grid data;

[0069] Specifically, the exported .GRDECL file is used, which contains grid geometry information and physical property data (such as porosity PORO, permeability PERMX, etc.). The script locates the keyword section containing the physical property data and extracts the corresponding physical property array. The specific execution flow is as follows: The script traverses the first grid data file and builds a "core keyword list" (such as SPECGRID, COORD, ZCORN). When the read keyword is not in this list and belongs to a predefined physical property alias library (for example, a configuration table that maps common synonym keywords to standard attribute names, such as uniformly recognizing PORO and POROSITY as poro, and uniformly recognizing PERMX and PERM_X as permx), it is determined to be an attribute data block. The numerical stream following the keyword is read, and comments and newline characters in the file are automatically filtered. The extracted linear numerical stream (List) is converted into a NumPy array and then divided according to the grid dimension (N x N y N z Perform a reshape operation to generate a 3D attribute matrix Mprop[k][j][i] according to the priority order of Z, Y, X (i.e., k, j, i). Check if the total amount of extracted data is equal to N. x ×N y ×N z First, ensure that each spatial grid cell has a unique corresponding attribute value, completing logical alignment. Next, parse the attribute data of each grid cell and align it with the grid's geometric data (such as cylindrical coordinates, corner coordinates, etc.). The physical attribute data of each grid cell is matched according to its grid cell ID or index, ensuring a one-to-one correspondence between physical attributes and geometric data.

[0070] It should be noted that in this embodiment, the physical attribute alias library is predefined in the configuration module of the conversion script as a "key-value dictionary". This library aims to eliminate the barrier of inconsistent naming of the same physical attribute by different geological modeling software (such as Petrel, Eclipse, CMG, etc.). Key: Covers common attribute keyword variations in the geological industry. For example, for permeability, the key set includes 'PERMX', 'KX', 'PERM_X', 'PERMEABILITY_X', etc.; for porosity, the key set includes 'PORO', 'POROSITY', 'PHI', etc. Value: The corresponding standard identifier uniformly used internally by the target numerical analysis software (FLAC3D). For example, the above permeability variations are uniformly mapped to 'permx', and porosity variations are uniformly mapped to 'poro'. The script follows these standard processing logics during parsing: Standardized preprocessing: When the script reads any keyword from a .GRDECL file, it first converts it to all uppercase and removes leading and trailing spaces to eliminate recognition interference caused by case differences. Fuzzy matching mechanism: It traverses the keys in the alias library. The script supports not only exact matching (i.e., the read keyword exactly matches a key in the library) but also prefix matching. For example, if the read keyword is 'POROSITY_EFF' (effective porosity), and the library defines the base key 'PORO' and allows prefix matching, the system can automatically identify and map it as a porosity attribute. This greatly improves compatibility with user-defined attribute names. Priority and conflict handling: The traversal order of the alias library is pre-sorted according to the importance or frequency of the attribute. Once a match is successful, the corresponding standard identifier is immediately returned and subsequent searches stop to prevent ambiguity. Fallback: If no match is found after traversing the entire alias library, the script will automatically classify it as a "non-target attribute" or "custom comment field", skip reading the data block, and log an "Unknown Keyword" warning without interrupting the entire conversion process.

[0071] It should also be noted that, in this embodiment, the matching rule for the alias library is strictly defined as "standard key leader matching". The specific rule is as follows: Let the attribute keyword to be identified be S. input (For example, 'POROSITY_EFF' read from a file), the standard key in the alias library is K. std (For example, predefined 'PORO' or 'POROSITY'). The script will iterate through the alias library if and only if S input With K stdWhen the key is the beginning of the name, a successful match is determined. Example 1: Keyword read: 'POROSITY_EFF' (Effective Porosity) Standard key defined in the library: 'POROSITY' Determination: Because 'POROSITY_EFF' contains the prefix 'POROSITY', a successful match is determined. The system maps it to the standard attribute 'poro'. Example 2 (for abbreviations): Keyword read: 'PERMX_REAL' Standard key defined in the library: 'PERMX' Determination: A successful match is determined, mapped to the standard attribute 'permx'. The advantage of this rule is that it can automatically accommodate derived attribute names with suffixes (such as _EFF, _LOG, _RAW) without exhaustively listing all possible suffix combinations in the library. If multiple keys can match (such as 'P' and 'PORO'), the script prioritizes matching the longer key according to a preset key length priority principle to ensure accuracy.

[0072] S25. Using the cylindrical trajectory in the COORD array, calculate the three-dimensional coordinates of each corner point through linear interpolation.

[0073] Specifically, in the previous steps, the top and bottom coordinates of each bar have been extracted from the COORD array. The coordinates of each bar consist of six floating-point numbers: X... top ,Y top Z top (Top coordinates) and X bottom ,Y bottom Z bottom (Bottom coordinates). To estimate the spatial position of the cylinder at different Z values ​​between known top and bottom coordinates, linear interpolation is also chosen here. This method is computationally efficient and simple, suitable for interpolating cylinder coordinates in a grid. Linear interpolation calculates the coordinates of the intermediate point using two known points (the top and bottom coordinates of the cylinder). For each grid cell, its eight corner points are traversed. Each corner point has a coordinate on the Z-axis (from Z...). top To Z bottom We calculate the corresponding X and Y coordinates for each corner point using the linear interpolation formula described above. The specific steps are as follows: For corner point 1, its Z coordinate is Z1, and we calculate its X1 and Y1 using the formula. For corner point 2, its Z coordinate is Z2, and we calculate its X2 and Y2 using the formula. Similarly, for the remaining corner points 3 to 8, we calculate the X and Y coordinates sequentially using the same interpolation method. Storing the 3D coordinates of each corner point: The calculated 3D coordinates (X, Y, Z) of each corner point will be stored in a data structure. This data structure ensures that the coordinates of each corner point in each grid cell can be accurately recorded and indexed.

[0074] Furthermore, calculating the x and y coordinates of a point at a given depth z can be expressed as:

[0075] ,

[0076] ,

[0077] Where X(z) and Y(z) are the X and Y coordinates of the corner point; z is the height of the corner point on the Z-axis; X top Y top Z top These are the x, y, and z coordinates of the center of the top of the column; X bottom Y bottom Z bottom These are the x, y, and z coordinates of the center of the bottom of the column. Using the Z value and the coordinates of the top and bottom of the column, the X and Y coordinates of the corner point are calculated.

[0078] S26. Obtain the original dataset including column coordinates, corner coordinates, and physical property data.

[0079] Specifically, the script extracts and calculates the column coordinates, corner coordinates, and physical property data through the aforementioned steps, then organizes and aligns this data by grid cell. Finally, this data is stored in a new data structure, providing complete data support for subsequent grid topology rearrangement, property mapping, and numerical simulation.

[0080] S3. Using the original dataset, the corner grid data is converted into hexahedral cell grid data by topology rearrangement and cell type matching methods, and the target grid topology sequence is obtained.

[0081] Furthermore, the specific steps for topology rearrangement and cell type matching are as follows:

[0082] S31. Read the ACTNUM array, and based on the array's labels, remove invalid cells, including those with zero volume and those with overlapping corners, and filter out valid cells.

[0083] Specifically, in this process, the script first reads the ACTNUM array, which marks whether each mesh cell is valid. Specifically, each cell in the ACTNUM array with a value of 1 is valid, while cells with a value of 0 or other markers are invalid. The script identifies and excludes invalid cells by parsing this array. In addition, the script checks the volume of each mesh cell; if the volume is zero, the cell is considered invalid. Similarly, if the corners of a mesh cell coincide, causing an abnormal geometry, it is also considered invalid. Through these checks, the script filters out valid mesh cells and generates a mapping table containing valid cell IDs for subsequent steps such as topology rearrangement and attribute mapping.

[0084] When constructing hexahedral elements, the script introduces a volume detection mechanism. The geometric volume V of each reconstructed element is calculated using the vector mixture product formula. If V≈0, it is identified as a pinch-out or degenerate element, and its removal is determined based on the ACTNUM flag. If V<0, it is identified as a node order disorder (negative volume), and the script automatically activates the node flip logic algorithm to adjust the node index order to a positive order conforming to the right-hand screw rule, ensuring the topological integrity of the numerical analysis model.

[0085] S32. Based on the corner coordinates obtained in S2, rearrange the corner order of each effective unit using the corner order replacement method;

[0086] Specifically, the script first obtains the corner coordinates of each valid cell, calculated based on the parameterized equations of the cylinders in step S2. Then, according to the corner order requirements of the numerical analysis software, the script applies a specific permutation algorithm to adjust the corner order of each grid cell. The core of the corner order permutation method is to rearrange the corners of the grid cells according to the order specified by the target software, ensuring that the corner order conforms to the topological format required by the software and avoiding incorrect processing of grid data due to incorrect corner order. More specifically, based on the corner coordinates obtained from S2, initial logical units are generated, and a corner index mapping table is established. For each grid unit (i,j,k), eight initial corners are first generated according to the logical offset, and the following "logical index-geometric position" correspondence table (index 0-7) is established in memory: Index 0: (i,j,k) — corresponds to the local coordinate origin; Index 1: (i+1,j,k) — corresponds to the positive offset point of the X-axis; Index 2: (i,j+1,k) — corresponds to the positive offset point of the Y-axis; Index 3: (i,j,k+1) — corresponds to the positive offset point of the Z-axis; Index 4: (i+1,j+1,k) — corresponds to the diagonal point of the XY plane; Index 5: (i,j+1,k+1) — corresponds to the diagonal point of the YZ plane; Index 6: (i+1,j,k+1) — corresponds to the diagonal point of the XZ plane; Index 7: (i+1,j+1,k+1) — corresponds to the far diagonal point. The initial sequence is denoted as Seq. init =[0,1,2,3,4,5,6,7].

[0087] S33. According to the right-handed spiral order specific to the numerical analysis software, the corner index is replaced to ensure that the corner order of each grid cell conforms to the topology requirements of the target software.

[0088] Specifically, the right-hand spiral order is a standard corner point arrangement, especially in hexahedral elements, where the corner point order determines the geometry and physical property distribution of the mesh element. Taking FLAC3D as an example, the corner point order requirements for hexahedral elements are as follows: NW t (Northwest Peak), NE t(Northeast Peak), SW t (Southwest Peak), SE t (Southeast Peak), NW b (Northwest bottom), NE b (Northeast Bottom), SW b (Southwest bottom), SE b (Southeast bottom). These corner points are numbered and arranged in a right-hand spiral order to ensure that the corner point order conforms to the topology requirements of the target software. Based on the corner point order requirements of the target numerical analysis software, an index permutation matrix is ​​created to record the mapping relationship between the original mesh corner point order and the order required by the target software. For example, corner point 1 of the original mesh is mapped to NW of the target mesh. t Corner point 2 is mapped to NE t And so on. The corner points of each valid grid cell are reordered using an index permutation matrix to ensure they meet the requirements of the target numerical analysis software. For example, corner points 1, 2, 3, and 4 in the original grid are adjusted to NW sequentially. t ,NE t ,SW t SE t Adjust corner points 5, 6, 7, and 8 to NW. b ,NE b ,SW b SE b After adjusting the corner order, the updated corner order is stored in the grid data structure for use in subsequent processing.

[0089] More specifically, numerical analysis software (such as FLAC3D) requires elements to satisfy the right-hand screw rule, meaning the element volume must be positive. The script calculates the volume V of the hexahedron formed by the initial sequence using the vector mixed product formula. Based on the calculation result, the following permutation rule is executed: Positive volume judgment (V>0): If the volume is positive, it means the initial sequence has satisfied the right-hand screw rule, and the Seq is maintained. init The output sequence remains unchanged, i.e., [0,1,2,3,4,5,6,7]. Negative volume correction (V<0): If the volume is negative, it indicates that the normal vector is involuted due to coordinate system definition or mesh distortion. In this case, "Flip Mapping" is activated to correct the normal vector by swapping the indices of adjacent nodes. The specific permutation mapping table is as follows: swap index 1 with index 2; swap index 5 with index 6; keep the remaining indices unchanged. The corrected sequence is: Seq flip =[0,2,1,3,4,6,5,7]. This explicit algorithm ensures that all elements imported into the numerical analysis software are valid positive volume elements, eliminating "negative volume" errors.

[0090] The volume V of a hexahedron is calculated by decomposing it into the sum of the volumes of five tetrahedrons. Let the eight corner points of the hexahedron be P0, P1, P2, P3, P4, P5, P6, and P7. The formula for the mixed product of vectors is:

[0091]

[0092] In the formula, V is the sum of the volumes of the tetrahedrons; the first term is the tetrahedron P0-P1P2P4 with P0 as the vertex; the second term is the tetrahedron P2-P3P0P6 with P2 as the vertex; the third term is the tetrahedron P4-P5P6P1 with P4 as the vertex; the fourth term is the tetrahedron P6-P7P4P3 with P6 as the vertex; and the fifth term is the inner core tetrahedron P2-P4P6P0. The solution is obtained by decomposing the solution into the sum of the volumes of the tetrahedrons.

[0093] S34. Generate a valid cell ID mapping table based on the correspondence between the original mesh and the target mesh;

[0094] Specifically, the script first reads the ID of each valid grid cell in the original grid. In the preceding step (e.g., S31), invalid cells have already been removed, so the original grid now contains only valid cells. Each valid cell has a unique ID, which will be used as input to the mapping table. For each valid grid cell, the script needs to determine its corresponding position in the target grid based on geometric information (e.g., corner coordinates, volume, topology, etc.). This process depends on the corner coordinates and order obtained in the previous steps, as well as the topological requirements of the numerical analysis software for the grid cells. By replacing the corner order and adjusting the geometry, the corresponding position of the original grid cell in the target grid can be accurately found. The mapping table uses a three-dimensional integer array (i.e., ZoneMap[N]) with the same dimensions as the geological grid. z ][N y ][N xThe array is stored as follows: The index (k,j,i) of the array represents the spatial location of the original geological grid, and the value of the array represents the unique number (Zone ID) of the cell in the numerical analysis model (FLAC3D). The specific steps for generating the mapping table are as follows: First, initialize a three-dimensional array ZoneMap with all zeros and a global counter CurrentID=0; then, traverse all grid cells according to the data storage order of the geological modeling software (usually a nested loop of K→J→I); during the traversal, check the ACTNUM array (validity marker). If ACTNUM[k,j,i]==0: mark the cell as invalid and skip it. If ACTNUM[k,j,i]==1: determine it as a valid cell. Execute CurrentID=CurrentID+1; assign CurrentID to the mapping table: ZoneMap[k][j][i]=CurrentID. The resulting ZoneMap establishes a bidirectional mapping from the original geological grid cells (with spatial indices i, j, k, stored in the array [k][j][i]) to the target numerical model cell IDs. Subsequent attribute writing can be quickly located by consulting this table, eliminating the need for repeated traversal and geometric matching. Based on the aforementioned geometric correspondence and the generated ZoneMap array, the script generates a mapping entry for each valid grid cell. The entry includes: Original Grid Cell ID: A unique identifier for each valid cell in the original grid. Target Grid Cell ID: A unique identifier for that cell in the target grid. Each row in the mapping table represents the mapping relationship for a valid cell, i.e., the correspondence between the cell IDs in the original grid and the cell IDs in the target grid.

[0095] S35. Based on the effective cell ID mapping table and the grid data after corner point order adjustment, and combined with the topology format required by the numerical analysis software, construct the target grid topology sequence that meets the requirements of the target numerical analysis software using the mapping table and corner point order arrangement results.

[0096] Specifically, the script first reads the valid element ID mapping table generated in step S34. This table contains the ID correspondence between the original and target meshes for each valid mesh element. Each entry in the mapping table contains the ID of the original mesh element and the ID of the target mesh element, ensuring that the geometric and attribute data of each valid element can be correctly mapped to the target mesh. Then, the mesh data adjusted in steps S32 and S33 is used, where the corner order has been replaced according to the requirements of the target numerical analysis software, ensuring that the corner order conforms to the software's topology requirements (e.g., right-hand spiral order). The geometric information (corner coordinates) and physical properties (e.g., porosity, permeability) of each mesh element have been adjusted in the previous steps. In this step, the script arranges the data of each mesh element into the target format according to the topology format requirements of the target numerical analysis software (e.g., the specific requirements of FLAC3D or other software). Specifically, the script organizes the data according to the element type required by the target software (e.g., hexahedral elements, tetrahedral elements, etc.), ensuring that the corner order of the elements is consistent with the target format. The script uses a valid cell ID mapping table to map the ID of each valid grid cell from the original grid to the target grid. Through this mapping table, the script correctly copies the geometric data (corner coordinates) and physical properties (such as porosity and permeability) of each valid cell in the original grid to the target grid cell. Following these steps, the script combines the corner coordinates, physical properties, and ID of each valid grid cell to form a new target grid topology sequence. This target grid topology sequence includes the cell order, corner order, and corresponding physical property data required by the target numerical analysis software. At this point, the data for each grid cell fully meets the requirements of the target numerical analysis software.

[0097] Through the above steps, based on the effective cell ID mapping table, the mesh data with adjusted corner order, and the topology format required by the target numerical analysis software, a target mesh topology sequence conforming to the software's requirements was constructed. This target mesh topology sequence contains the geometric information, physical properties, and correct corner order of each effective mesh cell, ensuring that the mesh data can be successfully imported into the target numerical analysis software and used for subsequent numerical simulations and calculations.

[0098] S4. Accurately map the attribute data of the target mesh topology sequence, precisely aligning the physical attribute data with the geometric data, and generate a second format file (.f3dat). Before performing physical attribute mapping, the script performs spatial dimension matching verification on each attribute array (such as porosity and permeability). It also checks whether the total number of read attribute values ​​is strictly equal to the total number of mesh cells N. total =N x ×N y ×N z .

[0099] In addition, numerical range cleaning is performed: for physically impossible values ​​(such as negative porosity, infinite permeability, or non-numeric NaN), the script resets them to the default background value or marks a warning in the log to prevent abnormal parameters from causing the numerical simulation to fail to converge.

[0100] Furthermore, the specific steps for accurately mapping the attribute data of the target mesh topology sequence are as follows:

[0101] S41. Initialize the counter ZoneID, and based on the valid cell IDs in the target mesh topology sequence, perform a geometric information traversal of the spatial cells and calculate the attribute data index of each mesh cell.

[0102] Specifically, first, a counter, ZoneID, is initialized to uniquely identify each grid cell. In practice, ZoneID typically starts from 1 and increments sequentially to ensure each grid cell has a unique identifier. The process iterates through the valid cell IDs in the target grid topology sequence. This topology sequence has been constructed according to the format required by the target numerical analysis software and includes the geometric information (such as corner coordinates, volume, etc.) and other relevant data for each grid cell. At this point, the grid has been adjusted for corner order to ensure each grid cell conforms to the topology requirements of the target software. For each valid cell, the script calculates its corresponding attribute data index based on its position (i,j,k) in the grid. Attribute data is typically stored as a one-dimensional array, with the array order being i→j→k in a cyclic manner. This means that the data in the attribute array is arranged according to the three-dimensional coordinates of the grid, and the index position corresponds to the physical attribute value (such as porosity, permeability, etc.) in the attribute data array. At this point, the geometric data and attribute data of each grid cell achieve a one-to-one correspondence in spatial and physical attribute data through this index.

[0103] More specifically, the formula for the index position is as follows:

[0104] ,

[0105] Where i, j, k are the indices of the grid cells in the X, Y, and Z directions; NX, NY, NZ are the dimensions of the grid. This formula is based on the storage convention of attribute data in the GRDECL format, i.e., linearly arranged in the order of i→j→k. Using this formula, the script can accurately calculate the position of the attribute data for each grid cell in the attribute array.

[0106] S42. Verify the validity of each grid cell by reading the ACTNUM array;

[0107] Specifically, the script verifies the validity of a mesh cell by reading the flags in the ACTNUM array. The ACTNUM array records the validity flag for each cell, where a valid cell is marked with 1, indicating that the cell needs to be processed in the physical calculation; an invalid cell is marked with 0, indicating that the cell should be discarded. If the cell is valid (ACTNUM is marked with 1), the script reads the corresponding physical attribute value from the attribute array based on the attribute data index. For example, attribute values ​​such as porosity (PORO) or permeability (PERMX) will be read and used for subsequent attribute mapping.

[0108] S43. During the mapping process, a corresponding FLAC3D command string is generated for each valid mesh element;

[0109] Specifically, the script first obtains the ZoneID for each valid cell and the physical property values ​​already calculated in the preceding steps, such as porosity or permeability. Then, according to FLAC3D's input format, each physical property is organized into a corresponding command string, containing the property name, property value, and the cell region identifier (i.e., ZoneID) to which the property applies. For example, when a cell's porosity is 0.25, the command format "zone property porosity 0.25 rangeid 1" will be generated. If a cell contains multiple physical properties, the script will generate multiple corresponding command strings.

[0110] S44. Generate a second format file whose physical properties perfectly match the geometric information of each mesh cell.

[0111] Furthermore, the method for verifying the validity of each grid cell by reading the ACTNUM array is as follows:

[0112] When ACTNUM[i,j,k]=1, the cell is a valid cell and can be used for attribute mapping;

[0113] When ACTNUM[i,j,k]=0 and other tags are used, the cell is an invalid cell and is removed.

[0114] S5. Import the second format file into the target numerical analysis software to automatically generate a numerical analysis calculation grid model containing complete physical properties and output the calculation results for numerical simulation and analysis.

[0115] Specifically, once the file is successfully imported, the target numerical analysis software will automatically generate a computational mesh model based on the geometric and physical property data in the file. In software such as FLAC3D, the computational mesh model consists of a series of mesh elements (e.g., hexahedral elements, tetrahedral elements, etc.), and the geometric information and physical properties of each element have been aligned through the aforementioned steps. The software determines the position and shape of the mesh elements based on the geometric information (such as corner coordinates), and further generates the mesh model structure for numerical calculation through the target software's solver.

[0116] Furthermore, the method is accomplished using a conversion script implemented in the Python programming language.

[0117] Furthermore, the physical property parameters in S1 include one or more of pore pressure, porosity, permeability, Young's modulus, and Poisson's ratio.

[0118] Furthermore, the geological modeling software includes Petrel, Roxar RMS, SKUA-GOCAD, and ECLIPSE.

[0119] Furthermore, the numerical analysis software includes FLAC3D, 3DEC, and PFC3D.

[0120] The embodiments described above are merely illustrative of specific implementations of the present invention, and while the descriptions are detailed, they should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the scope of protection of the present invention.

Claims

1. A method for converting and mapping a three-dimensional geological grid model to a numerical analysis model, characterized in that, Includes the following steps: S1. Export a first grid data file containing grid geometry information and at least one physical property parameter from geological modeling software. The first grid data file is in .GRDECL format. S2. Perform deep analysis on the first grid data file to extract the column coordinates and corner coordinates, and obtain the original dataset of the grid geometry. S3. Using the original dataset, the corner grid data is converted into hexahedral cell grid data by topology rearrangement and cell type matching methods, and the target grid topology sequence is obtained. S4. Accurately map the attribute data of the target mesh topology sequence, precisely align the physical attribute data with the geometric data, and generate a second format file, the second format file being in the .f3dat format; S5. Import the second format file into the target numerical analysis software to automatically generate a numerical analysis calculation grid model containing complete physical properties and output the calculation results for numerical simulation and analysis.

2. The method according to claim 1, characterized in that, The specific steps for deep parsing the first grid data file include: S21. Read the COORD array from the first grid data file, divide the data in the COORD array into groups of 6 floating-point numbers, and extract the three-dimensional coordinates of the top and bottom of the column one by one; S22. By traversing the grid dimensions, construct the parameterized equation for each bar based on the bar coordinates in the COORD array, which is used to describe the trajectory of the bar in three-dimensional space. S23. Read the ZCORN array from the first grid file, and extract the depth values ​​of the 8 corner points of each grid cell from the ZCORN array using an indexing algorithm; S24. Read physical property data from the first grid data file and align this property data with the grid data; S25. Using the cylindrical trajectory in the COORD array, calculate the three-dimensional coordinates of each corner point through linear interpolation. S26. Obtain the original dataset including column coordinates, corner coordinates, and physical property data.

3. The method according to claim 1, characterized in that, The specific steps for topology rearrangement and cell type matching are as follows: S31. Read the ACTNUM array, and based on the array's labels, remove invalid cells, including those with zero volume and those with overlapping corners, and filter out valid cells. S32. Based on the corner coordinates obtained in S2, rearrange the corner order of each effective unit using the corner order replacement method; S33. According to the right-handed spiral order specific to the numerical analysis software, the corner index is replaced to ensure that the corner order of each grid cell conforms to the topology requirements of the target software. S34. Generate a valid cell ID mapping table based on the correspondence between the original mesh and the target mesh; S35. Based on the effective cell ID mapping table and the grid data after corner point order adjustment, and combined with the topology format required by the numerical analysis software, construct the target grid topology sequence that meets the requirements of the target numerical analysis software using the mapping table and corner point order arrangement results.

4. The method according to claim 1, characterized in that, The specific steps for accurately mapping the attribute data of the target grid topology sequence are as follows: S41. Initialize the counter ZoneID, and based on the valid cell IDs in the target mesh topology sequence, perform a geometric information traversal of the spatial cells and calculate the attribute data index of each mesh cell. S42. Verify the validity of each grid cell by reading the ACTNUM array; S43. During the mapping process, a corresponding FLAC3D command string is generated for each valid mesh element; S44. Generate a second format file whose physical properties perfectly match the geometric information of each mesh cell.

5. The method according to claim 4, characterized in that, The formula for calculating the attribute data index of each grid cell in S41 is as follows: ， Where i, j, k are the indices of the grid cells in the X, Y, Z directions; NX, NY, NZ are the dimensions of the grid.

6. The method according to claim 4, characterized in that, The method for verifying the validity of each grid cell by reading the ACTNUM array is as follows: When ACTNUM[i,j,k]=1, the cell is a valid cell and can be used for attribute mapping; When ACTNUM[i,j,k]=0 and other tags are used, the cell is an invalid cell and is removed.

7. The method according to claim 1, characterized in that, The method described uses the Python programming language to implement script conversion.

8. The method according to claim 1, characterized in that, The physical property parameters in S1 include pore pressure, porosity, permeability, Young's modulus, and Poisson's ratio.

9. The method according to claim 1, characterized in that, The geological modeling software includes Petrel, RoxarRMS, SKUA-GOCAD, and ECLIPSE.

10. The method according to claim 1, characterized in that, The numerical analysis software includes FLAC3D, 3DEC, and PFC3D.

Images