A method of meshing an image after reconstruction

Abstract

The application discloses a kind of gridding methods after image reconstruction, comprising the following steps, 1, the three-dimensional reconstruction of object sample geometric structure image is carried out, and scanning data is processed as the object image of different section;2, image pre-processing;3, the object is segmented into multiple phases, and the segmented image slice is three-dimensionally stacked and reconstructed, or further to the core after three-dimensionally stacked and reconstructed rock core generates filler;Digital core is cut into lengthX again×lengthY×lengthZ solid unit small core and stored in the corresponding core file;4, the generated core file is imported in MATLAB grid program, and one-dimensional array is recombined into three-dimensional calculation matrix;5, grid unit numbering and generating node coordinates are carried out to calculation matrix;6, face grid calculation is carried out to grid unit;7, grid file is written, and the grid file is converted into the grid format accepted by application software in Hypermesh.

Description

Technical Field

[0001] This invention belongs to the field of image processing technology, specifically relating to a method for reconstructing and meshing numerical images of objects obtained by means of photography, CT scanning, MRI, electron microscopy, etc., including three main parts: image processing, reconstruction, and meshing. Background Technology

[0002] In physics experiments, challenges exist such as difficulties in preparing laboratory samples, high equipment manufacturing or experimental costs, and difficulties in analyzing microscopic mechanisms. With the rapid development of computer technology, numerical simulation technology has gradually become an experimental analysis method of equal importance to theoretical analysis, laboratory experiments, and field experiments. The construction of geometric meshes is the foundation of numerical simulation. The type, connection method, and accuracy of the mesh determine the iterative transfer of constitutive equations, numerical interpolation between nodes, and computational accuracy in numerical simulation.

[0003] As basic science and engineering science deepen, there is an increasing demand for research on solid mechanics, fluid flow, heat transfer, mass transfer, chemical reactions, and electromagnetic properties of actual physical structures. This is to meet the needs of complex engineering projects, and numerical simulation technology has proven its advantage in revealing microscopic mechanisms. In the past, various mathematical functions were used to construct complex microscopic geometric models, but these still fell far short of real objects. Therefore, emerging technologies in recent years utilize techniques such as imaging, CT scans, MRI, and electron microscopy to obtain tomographic images of real objects, and then use reconstruction techniques to quantify these real objects. This addresses the need for numerical simulations to capture the structure of real objects. However, for use in numerical simulation experiments, the constructed numerical model still requires necessary meshing to facilitate the subsequent construction and iterative calculations of the geometric model in the numerical simulation software.

[0004] Chinese patent document CN 114494637 A, published on May 13, 2022, discloses a method for reconstructing a three-dimensional realistic model of sandstone based on a structure matrix. This method includes: acquiring geometric and topological image data of a standard sandstone specimen using X-ray imaging equipment; converting the image data into a discrete three-dimensional digital matrix; performing two-phase segmentation of the sandstone digital matrix using the Chan-Vese model; establishing a structure matrix pattern library: determining the corresponding central voxel unit pattern and coordinate information based on the structure matrix of the sandstone data; obtaining geometric and topological data and linear interpolation through structure matrix inversion to obtain the point coordinates and direction matrix of the two-phase microstructure of the sandstone, thereby realizing the visualization of the sandstone three-dimensional model reconstruction process and results.

[0005] This patent only focuses on reconstructing a realistic 3D sandstone model to visualize the digital core, without meshing the digital core, thus hindering subsequent scientific analysis. While existing technologies, such as Avizo image processing software, offer meshing for extracting the pore network model (a ball-and-stick model), the resulting meshes are difficult to adapt to other commercial software due to format and accuracy issues. Furthermore, in fields like electrical properties and modulus calculations, open-source software is typically used to directly perform calculations on the digital core. This requires researchers to program not only the mesh nodes but also the numerical simulation components, which is detrimental to multi-field coupled numerical simulations in solid mechanics, fluid flow, and heat transfer. Therefore, there is an urgent need for an image-based meshing method to improve the applicability of numerical models for microscale multi-field coupled numerical simulations. Summary of the Invention

[0006] To address the problems existing in the prior art, the technical problem to be solved by the present invention is to provide a meshing method after image reconstruction, which can solve the problem of meshing numerical models and establish a complete process for image processing, reconstruction, and meshing of physical samples.

[0007] The technical problem to be solved by the present invention is achieved through the following technical solution, which includes:

[0008] Step 1: Perform CT scans, nuclear magnetic resonance, and electron microscopy experiments on the object sample to obtain geometric images of the surface or internal fractures of the object as required by numerical simulation, and process the scan data into object images of different cross sections through three-dimensional reconstruction.

[0009] Step 2: Image preprocessing is performed on the tomographic image of the object, including image segmentation, contrast and brightness adjustment, and filtering and noise reduction. Factors that are not conducive to meshing, such as the outer corners of the object and scanning noise disturbances, are removed in advance.

[0010] Step 3: Perform grayscale analysis on the preprocessed image, divide the object into multiple phases through dual threshold segmentation, including skeletal particles and pore space, and perform three-dimensional stacking reconstruction on the segmented image slices. Then, perform filtering and noise reduction processing again, and save it as the original XXX.raw format three-dimensional image; or further generate filling material on the three-dimensional stacked and reconstructed core.

[0011] The stacked or reconstructed digital cores with filling are then cut into smaller cores of lengthX×lengthY×lengthZ, and the core data are stored in the corresponding core files in ascending order of x, y, z.

[0012] Step 4: Import the core file generated in Step 3 into the MATLAB grid subroutine. The core file is a one-dimensional array arranged in ascending order of x, y, and z. Set the grid pixel size (lengthX, lengthY, lengthZ), recombine the one-dimensional array into a three-dimensional calculation matrix, and judge the correctness of the reconstructed data through three-dimensional matrix slice image visualization processing.

[0013] Step 5: Number the grid cells and generate node coordinates for the computation matrix.

[0014] The meshing subroutine numbers nodes sequentially from position (1, 1, 1) in the 3D computation matrix to (lengthX, lengthY, lengthZ). Each mesh cell generates 8 nodes, resulting in a node matrix of size (lengthX+1, lengthY+1, lengthZ+1). It calculates node coordinates as node matrix coordinates multiplied by resolution and sorts the calculated node coordinates according to the cell number. The node coordinates are stored in a matrix Na(lengthX+1, lengthY+1, lengthZ+1), while the mesh cell numbers are stored in cell1, cell2, and cell3 according to different phases.

[0015] Step 6: The MATLAB mesh subroutine performs surface mesh calculations on mesh elements according to the phase of the skeleton particles, pore spaces, or filling material. The connection format of the six surface meshes of a mesh element is (n1 n2 n3 n4 c1 c2), where ni is the node that makes up the surface mesh and the connection order, and ci is the neighboring mesh element of the surface mesh; the node connection format of the two-dimensional mesh element is (n1 n2 c1c2); the surface meshes are stored in facem1, facem2, or facem3 according to different phases.

[0016] Step 7: After the node coordinates and surface mesh calculations are completed, write the mesh file in the standard Fluent mesh format and save the generated mesh file as XXX.msh. In Hypermesh, convert the XXX.msh file into a mesh format that the application software can accept.

[0017] The technical effects of this invention are:

[0018] This invention can reconstruct and mesh various images, enabling numerical simulation studies of solid mechanics, fluid mechanics, thermodynamics, chemistry, and electromagnetics using the microstructure of real objects in the images. Attached Figure Description

[0019] The accompanying drawings of this invention are described below:

[0020] Figure 1 This is a schematic diagram of the image reconstruction and meshing process;

[0021] (a) Fault section image of CT scan core;

[0022] (b) Image preprocessing slicing;

[0023] (c) Three-dimensional reconstruction diagram;

[0024] (d) Hypermesh mesh visualization;

[0025] (e) Numerical simulation application diagram after image meshing;

[0026] Figure 2 Flowchart for filling pores with hydrates;

[0027] Figure 3 for Figure 2 The flowchart of the subroutine for loop search and hydrate generation in the process;

[0028] Figure 4 The interface diagram for running a 50×50×50 three-dimensional three-phase meshing program. Detailed Implementation

[0029] The present invention will be further described below with reference to the accompanying drawings and embodiments:

[0030] An example includes the following steps:

[0031] Step 1: Perform CT scans, nuclear magnetic resonance, and electron microscopy experiments on the object sample to obtain geometric images of the surface or internal fractures of the object as required by numerical simulation, and process the scan data into different fracture slice images through three-dimensional reconstruction.

[0032] like Figure 1 As shown in (a), the 3D image reconstructed from 479 1000×1000 pixel CT geometric structure images is a CT tomographic slice. See the literature "3D Reconstruction and Numerical Experimental Study of Asphalt Mixture Based on X-ray CT and Finite Element Method", Wan Cheng, South China University of Technology, pp. 20-44, April 2010, which describes a method for reconstructing images from CT scans.

[0033] Step 2: Image preprocessing is performed on the tomographic slice image to remove the non-core parts on the outside of the image, contrast and brightness are adjusted, and filtering and noise reduction are performed to remove the factors that are not conducive to meshing caused by scanning noise in advance.

[0034] like Figure 1As shown in (b), in this embodiment, the tomographic slice image is cut into a slice sequence with a length × width of 331 × 331 using image processing software such as IMAGEJ and Avizo, and the non-core parts on the outside of the sample are removed; then the brightness, contrast and other adjustments are made to the image to improve the resolution between core skeleton particles, pores and filling materials, and at the same time the image is filtered and denoised to facilitate core segmentation.

[0035] Step 3: Perform grayscale analysis on the preprocessed image, divide the object into multiple phases through dual threshold segmentation, including skeletal particles and pore space, and stack the segmented image slices in sequence to reconstruct a three-dimensional digital core. Then, perform filtering and noise reduction again to remove unnecessary pore structures caused by threshold segmentation, and save it as the original XXX.raw format three-dimensional image.

[0036] See the literature "Pore Fractal Characteristics of Hydrate-Bearing Sands and Implications to the Saturated Water Permeability", Zhun Zhang et al, JGR (Journal of Geophysical Research): Solid Earth, 125, pp. 5-6, Mar, 2020, (“Pore Fractal Characteristics of Hydrate-Bearing Sands and Their Influence on Saturated Water Permeability”, Zhang Zhun, Journal of Geophysical Research, Vol. 125, pp. 5-6, March 2020), which describes a dual-threshold segmentation method. This image segmentation method uses dual-threshold segmentation to separate the skeletal particles and pore space. The segmented image is then binarized for subsequent operations, with skeletal particles recorded as "0" and pore space as "1".

[0037] like Figure 1 (c) As shown in the embodiment, dual threshold segmentation is used. By thresholding the grayscale values ​​of the preprocessed image in step 2, the core is divided into two parts: skeletal particles and pores. After segmentation, the CT tomographic slice image sequence is stacked in the slice arrangement order to form a three-dimensional digital core, and the core data is saved as the original image data of three-dimensional A1particle.raw.

[0038] Furthermore, based on the saturation of the pore filler, the pore filler generated after 3D reconstruction is processed by either randomly replacing pixels without considering the generation conditions or manually controlling the replacement of pixels under the generation conditions (the conditions are determined by the specific pore filler generation mechanism).

[0039] Numerical values ​​for generating filled samples with different saturations are generated through pixel replacement calculations. The saturation of the filled phase is calculated as: (Number of pixels in the filled phase) ÷ (Number of pixels in the filled phase + Number of pixels in the remaining pores). This filler generation step overcomes the problem of limited saturation in laboratory-generated multiphase porous media samples, enabling precise control of saturation within a target range.

[0040] In this embodiment, the pore space is segmented using a watershed method, dividing the pores into pore blocks of a certain size according to the throat and pores. The geometric center coordinates and average pore radius of the pore blocks are calculated and stored. In addition, the skeleton space is segmented using a watershed method to obtain the skeleton center coordinates and skeleton particle radii. See the literature "Versatile and efficient pore network extraction method using marker-based watershed segmentation", Jeff T. Gostick, Physical Review E 2017 ("A general and efficient pore network extraction method based on marker-based watershed segmentation", Jeff T. Gostick, Physical Review E, No. 96, pp. 1-13, August 2017), which describes a watershed method for segmenting pore space. Through basic operations such as calculating Euclidean distance, peak picking, peak removal, and merging of neighboring peaks, the pores are finally segmented into pore blocks of a certain size according to throats and pores. The geometric center coordinates and average pore radius of the pore blocks are calculated and stored. In addition, the skeleton space is segmented by watershed to obtain the skeleton center coordinates and skeleton particle radii.

[0041] The following steps are involved in developing a pore-filling program based on MATLAB to generate the filler:

[0042] 1) Import the XXX.raw file, pore information file, and skeleton information file, read the image data into a three-dimensional matrix, and form a calculation matrix a(lenx,leny,lenz) to provide basic skeleton and pore data for generating the infill material;

[0043] 2) Calculate the pore search range, which is a spherical range with a radius of ratio*plr from the pore center outward; and cyclically search for pore pixels with the pore center in mind.

[0044] 3) Generate filling material in the pores within the pore search range: The pore filling material extends from the pore center to the framework particles until the minimum gap between the filling material and the framework particles is not less than one voxel.

[0045] 4) After the filling material of one pore is generated, the next pore search begins. After all pores are searched, the filling material saturation is calculated. If the filling material saturation does not reach the preset saturation, the pore search range and the preset pore radius threshold plr0 are adjusted until the preset filling material saturation is met. When the filling material saturation reaches the preset saturation, the filling material generation is complete.

[0046] like Figure 2 As shown, the hydrate generation procedure follows these steps:

[0047] In S101, import the image data from the A1particle.raw file;

[0048] In S102, import the pore information file and the skeleton information file;

[0049] The pore information includes the geometric center coordinates of the pore block and the average pore radius; the skeleton information includes the skeleton center coordinates and the skeleton particle radius.

[0050] In S103, set the dimensions of the 3D image: lenx, leny, lenz.

[0051] In S104, the image data from S101 is read into a three-dimensional matrix to form a computation matrix a(lenx,leny,lenz);

[0052] In S105, the expected hydrate saturation Sh0 and the pore radius threshold plr0 are set to control the hydrate formation saturation and determine whether pore filler can be formed in pores of different sizes.

[0053] In S106, the pore information is used, and the pore pixels are searched cyclically within a "sphere" with a radius of ratio*plr from the pore center outward.

[0054] In S107, hydrates are generated in the pores, and the pore pixel value is replaced with "7";

[0055] In S108, calculate the hydrate saturation Sh = hydrate phase volume / (hydrate phase volume + remaining pore volume);

[0056] In S109, determine whether the hydrate saturation Sh has reached the preset saturation Sh0. If yes, proceed to S111; otherwise, proceed to S110.

[0057] In S110, adjust the pore search range ratio*plr and the preset pore radius threshold plr0, and return to S105;

[0058] Adjust the hydrate limit saturation to below 40% to meet the actual saturation limit requirements (20% to 40%) for pore-filled natural gas hydrates.

[0059] In S111, pore-filling hydrate is completed.

[0060] The subroutines for cyclical search and hydrate generation in S106 and S107 above are as follows: Figure 3 As shown, this subroutine follows these steps:

[0061] In S201, each pore pixel is searched according to its different pore center coordinates;

[0062] In S202, determine whether the pore radius is greater than the threshold plr0. If so, proceed to S203; otherwise, return to S201.

[0063] In S203, calculate the search radius of the pores, ratio*plr;

[0064] In S204, using the pore center coordinates as the "center of the sphere," search for pore pixels within a sphere with a radius of ratio*plr. This step can significantly reduce the computational load of hydrate formation cycles and judgments.

[0065] In S205, determine whether the pore pixels within the search range are in contact with the skeleton material. If so, proceed to S207; otherwise, proceed to S206.

[0066] At least one pixel's spacing must be maintained between hydrates and skeleton particles.

[0067] In S206, a hydrated material is generated by replacing the pore pixel value with "7", and then proceed to S208.

[0068] In S207, the cycle continues to the next aperture pixel;

[0069] In S208, determine whether all the aperture pixels within the "sphere" range of radius ratio*plr have been searched. If so, proceed to S209; otherwise, return to S204.

[0070] In S209, the search loop for aperture pixels within the radius ratio*plr "sphere" ends, and the next aperture loop 106 begins.

[0071] use Figure 2 The hydrate generation procedure shown, without considering the physical and chemical conditions of generation, searches for pore pixels in the middle of the pores towards the framework particles. When the pore space is not in contact with the framework particles, the pore space pixels are replaced with filler material. This process continues until the filler saturation meets the preset saturation level, and then stops. Finally, numerical cores with hydrate saturation levels of 12.52%, 18.66%, 27.07%, and 31.39% are generated.

[0072] To ensure that the generated mesh meets the memory and CPU requirements of subsequent COMSOL multi-field coupled numerical simulation calculations (the computer used for this example has a 24-core CPU and 128GB of memory), the digital core obtained in this step is further divided into 125,000 small cores of 50×50×50 pixels (due to the limitations of COMSOL's computational scale, it is difficult to divide the cores into smaller ones). The core data is then stored in the corresponding core files (txt text files) in ascending order of x, y, z.

[0073] Step 4: Import the core file generated in Step 3 into the MATLAB grid subroutine. The core file is a one-dimensional array arranged in ascending order of x, y, z (a one-dimensional array is a set of data of the same type that is stored continuously in memory). Set the grid pixel size (lengthX, lengthY, lengthZ). The grid subroutine will recombine the one-dimensional array into a three-dimensional calculation matrix and use the visualization processing of CT scan slice images to determine the correctness of the data reconstructed in Step 3 (the correctness of threshold segmentation).

[0074] Note: The calculation matrix here is similar in form to the "3D calculation matrix" in step S104, but contains different information. The calculation matrix in step S104 mainly contains pore and skeleton information, which is used for filling the pores with filler. The calculation matrix here contains pore / skeleton and generated filler information, which is used for calculating the mesh.

[0075] As an example, the core file data is imported into the MATLAB meshing subroutine. The mesh pixel size [xyz] of the core sample is set to [50 50 50]; the pixel values ​​of the three phases [m1 m2 m3] are [1 7 3], which are the pixel values ​​marking each phase. In this case, the skeleton particles are marked as "1", the filler is marked as "7", and the pore space is "3". Then, the pixel scale is input, which is the actual physical size represented by one pixel. In this example, the pixel scale [s1 s2 s3] is used for mesh calculation: [1 1 1]. The initialization data input for the meshing program is as follows:

[0076] Clear all;

[0077] size1 = input(“Please enter the grid pixel size, format [x,y,z]:”); % Sample pixel size

[0078] lenx=size1(1); leny=size1(2); lenz=size1(3);

[0079] im = load(“hydrate.txt”); % Read CT core data

[0080] mass = input(“Please enter the pixel values ​​of the three-phase material, in the format [m1,m2,m3], fill in 0 for each single or double term:”); %Pixel calibration of different phases

[0081] m1=mass(1); m2=mass(7); m3=mass(3);

[0082] scal = input(“Please enter the scale (how many physical dimensions one pixel represents) [s1,s2,s3], default value 1:”); %Physical size scale

[0083] s1=scal(1); s2=scal(1); s3=scal(1);

[0084] sum1 = sum(im(:) == m1);

[0085] sum2 = sum(im(:) == m2);

[0086] sum3 = sum(im(:) == m3);

[0087] a=reshape(im,lenx,leny,lenz);

[0088] subplot(1,2,1);

[0089] imshow(a(:,:,2),[]);

[0090] The interface for running the gridded program is as follows: Figure 4 As shown, the MATLAB meshing subroutine first reassembles the read one-dimensional (125000) core file into a three-dimensional core calculation matrix a(50,50,50). It then displays the second slice of the calculation matrix a, which shows good agreement between the results and the original slice's porosity and skeletal grains, verifying the correctness of the reconstructed data.

[0091] Step 5: Number the grid cells and generate node coordinates for the computation matrix.

[0092] A grid cell refers to the smallest structure that makes up a grid. In step 3, small rock cores are cut into 50×50×50 sizes, and in step 4, the grid pixel size is set to [50 50 50]. Thus, one small rock core corresponds to one grid cell.

[0093] The meshing subroutine sequentially numbers the eight vertices of the cube mesh element from position (1, 1, 1) to (50, 50, 50) in the 3D calculation matrix, with the numbers increasing in x, y, z. Each mesh element generates 8 nodes (8 nodes form a structured mesh), and the final generated node matrix size is (51, 51, 51). The subroutine calculates the node coordinates as node matrix coordinates × resolution, i.e., node coordinates (x, y, z) = index(X, Y, Z) × scale(X0, Y0, Z0). The calculated node coordinates are then sorted according to the mesh element number.

[0094] The MATLAB mesh subroutine calculates mesh nodes and element numbers by looping through x, y, z as follows:

[0095]

[0096]

[0097] The node coordinates are stored in a matrix Na(51, 51, 51), and the grid cell numbers are stored in cell1, cell2, cell3 (corresponding to the material pixel values ​​[m1 m2 m3]) according to different phases (skeleton, pores, filler), and their size is determined by the actual number of cells of different phases.

[0098] Step 6: The MATLAB mesh subroutine performs surface mesh calculations on the mesh units according to the three phases of skeleton particles, pore space, and filling material. The surface mesh refers to the six faces that make up the three-dimensional mesh unit and the four edges of the two-dimensional mesh unit; making the meshes of different phases independent of each other facilitates subsequent numerical simulation processing.

[0099] The structured mesh with a face mesh has a node connection format of (n1, n2, n3, n4, c1, c2), where ni represents the nodes constituting the face mesh and their connection order, and ci represents the neighboring mesh cells. For example, mesh cell 1 (8 nodes: 1, 2, 4, 5, 7, 8, 10, 11 (sorted in ascending order of coordinates)) and mesh cell 2 (8 nodes: 2, 3, 5, 6, 8, 9, 11, 12) are two adjacent cells. The face (2, 5, 8, 11) and its neighboring cells 1 and 2 have their four nodes connected in the order of the four vertices of a square as (2, 5, 11, 8). This mesh is denoted as (2, 5, 11, 8, 1, 2). For mesh cell 2, the node connection format of this face mesh is (2, 5, 11, 8, 2, 1).

[0100] Four edges are extracted from the two-dimensional mesh element on the surface of the digital core for calculation. The node connection format of the two-dimensional mesh element is to remove n3 and n4.

[0101] The surface meshes are stored in facem1, facem2, and facem3 according to different phases until all surface meshes of the mesh element have been calculated.

[0102] Step 7: After calculating the node coordinates and surface mesh, write the mesh file in the standard Fluent mesh format, including mesh dimension information, node information, surface mesh element information, boundary surface mesh information, and mesh name. Save the generated mesh file as XXX.msh. In Hypermesh, convert the XXX.msh file into a numerical simulation mesh for COMSOL, ABAQUS, ANSYS, etc., for subsequent numerical simulation work.

[0103] Mesh dimensional information includes two-dimensional and three-dimensional information. Node information includes node coordinates and the number of nodes. Surface mesh cell information includes surface mesh connections, classification, and quantity.

[0104] Write the node coordinates and surface mesh information into a .msh file according to the Fluent mesh format (this format is publicly available). The mesh file is written as follows:

[0105] fid = fopen(strcat('outName.msh'), 'w'); % Mesh name declaration

[0106] fprintf(fid,'%s%d%s\n','(',0,'"Created by WANG from CT rawdate")'); %Grid author information declaration

[0107] fprintf(fid,'%s%d%d%s\n','(',2,3,')'); % Mesh geometry information

[0108] fprintf(fid,'%s%d%s\n','(',0,'"Note Section:")'); %% Node information

[0109] fprintf(fid,'%s%d%s%d%2x%2x%d%d%s\n','(',10,'(',0,1,(lenx+1)*(leny+1)*(le nz+1),0,3,'))');

[0110] fprintf(fid,'%s%d%s%d%2x%2x%d%d%s\n','(',10,'(',1,1,(lenx+1)*(leny+1)*(le nz+1),1,3,')');

[0111] fprintf(fid,'%s\n','(');

[0112] fprintf(fid,'%d%d%d\n',(i-1)*s1,(j-1)*s2,(k-1)*s3);

[0113] fprintf(fid,'%s\n','))');%%

[0114] fprintf(fid,'%s%d%s\n','(',0,'"faces:")');%% face mesh declaration

[0115] fprintf(fid,'%s%d%s%2x%2x%2x%2x%2x%s\n','(',13,'(',0,1,lenf1+lenf2+lenf3+lenf4+lenf5+lenf6+lenf7,0,0,'))');%face

[0116] %bar2=waitbar(0,'Writing interior1');

[0117] fprintf(fid,'%s%d%s\n','(',0,'"Interior faces of zone VOID")');

[0118] if lenf1==0

[0119] else

[0120] fprintf(fid,'%s%d%s%2x%2x%2x%2x%2x%s\n','(',13,'(',7,1,lenf1,2,4,')');%interior face

[0121] fprintf(fid,'%s\n','(');

[0122] %for i=1:lenf1

[0123] fprintf(fid,'%x%x%x%x%x%x\n',interior1);

[0124] str1=['Mesh calculation in progress

[0125] ',num2str(((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenf1)*100 / ((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenx*leny*lenz*3+lenx*leny+lenx*lenz+leny*lenz)),'%'];

[0126] waitbar(((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenf1) / ((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenx*leny*lenz*3+lenx*leny+lenx*lenz+leny*lenz),bar1,str1);

[0127] %str2=['interior1 being written...',num2str((i) / lenf1),'%'];

[0128] %waitbar((i) / lenf1,bar2,str2);

[0129] %end

[0130] fprintf(fid,'%s%s\n',')',')');

[0131] end%%

[0132] ……

[0133] %bar2=waitbar(0,'wall being written');%%Boundary grid declaration

[0134] if lenf7==0

[0135] else

[0136] fprintf(fid,'%s%d%s\n','(',0,'"wall")');

[0137] fprintf(fid,'%s%d%s%2x%2x%2x%2x%2x%s\n','(',13,'(',13,lenf1+lenf2+lenf3+lenf4+lenf5+lenf6+1,lenf1+lenf2+lenf3+lenf4+lenf5+lenf6+lenf7,3,4,')');%wall

[0138] fprintf(fid,'%s\n','(');

[0139] %for i=1:lenf7

[0140] fprintf(fid,'%2x%2x%2x%2x%2x%2x\n',wall);

[0141] str1 = ['Grid computing in progress']

[0142] ',num2str(((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenf1+lenf2+lenf3+lenf4+lenf5+lenf6+lenf7)*10 0 / ((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenx*leny*lenz*3+lenx*leny+lenx*lenz+leny*lenz)),'%'];

[0143] waitbar(((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenf1+lenf2+lenf3+lenf4+lenf5+lenf6+lenf7) / ((le nx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenx*leny*lenz*3+lenx*leny+lenx*lenz+leny*lenz),bar1,str1);

[0144] % str2 = ['walls under development...', num2str((i) / lenf7), '%'];

[0145] % waitbar((i) / lenf7,bar2,str2);

[0146] %end

[0147] fprintf(fid,'%s%s\n',')',')');

[0148] end%%

[0149] The final generated mesh file includes mesh declaration, mesh dimension information, node coordinates, face mesh information, mesh structure information, etc., as shown below:

[0150] (0"Created by WANG from CT raw date")%Author's Statement

[0151] (2 3)% grid dimension information

[0152] (0"Note Section:")% Node Declaration

[0153] (10(0 1 2062b 0 3))%% Node Information

[0154] (10(1 1 2062b 1 3) ( 0 0 0 1 0 0 200 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 9 0 0

[0166] ... 46 50 50 47 50 50 48 50 50 49 50 50 50 50 50

[0172] ))%%

[0173] 0"faces:")%%face mesh information

[0174] (13(0 1 5d624 0 0))

[0175] (0"Interior faces of zone VOID")(13(7 1 1802d 2 4)(

[0176] a60 a93 14bc 1489 2 1

[0177] 1488 1489 14bc 14bb 359 1a61 a94 14bd 148a 3 2a93 a94 14bd 14bc 18 2

[0178] 1489 148a 14bd 14bc 35a 2a62 a95 14be 148b 4 3a94 a95 14be 14bd 19 3

[0179] 148a 148b 14be 14bd 35b 3…… ))

[0181] ...

[0182] (0"wall")%% Boundary Mesh (13(d 59b8d 5d624 3 4)(

[0183] 1 34a5d a2a 0af30

[0184] 1 2a2b a2a 0af30

[0185] 1 2 35 34 0af30

[0186] 2 3a2c a2b 0af31

[0187] 2 3 36 35 0af31

[0188] 3 4a2d a2c 0af32

[0189] 3 4 37 36 0af32

[0190] ... ))

[0192] (0"cells:")%%Grid cell declaration(12(0 1 1e848 0 0))

[0193] (12(4 1 99ea 1 4))

[0194] (12(5 99eb af2f 1 4))

[0195] (12(6af30 1e848 1 4))

[0196] (0"Zone Sections:")%% Mesh Phase Information Declaration

[0197] (45(4fluid void)())

[0198] (45(5solid hydrate)())

[0199] (45(6solid solid)())

[0200] (0"faces Sections:")%% Faces Mesh Information Declaration

[0201] (45(7interior VF)())

[0202] (45(8interior HF)())

[0203] (45(9interior SF)())

[0204] (45(a interface VS)())

[0205] (45(b interface VH)())

[0206] (45(d wall WALLS)())%% Boundary Mesh Information Declaration

[0207] Figure 1 (d) is the mesh display diagram of the Hypermesh mesh processing software. In this embodiment, the mesh size was adjusted and the mesh format was converted in Hypermesh, and the converted .nas format mesh file was imported into COMASOL to carry out subsequent numerical simulation work.

[0208] The core of this invention is the meshing of the reconstructed image, and specific applications of this meshing include... Figure 1 As shown in (e), this is the geometric display diagram of the COMSOL multi-field coupling numerical simulation application. The mesh file was imported into COMSOL and a fluid-structure interaction numerical simulation study was carried out.

Claims

1. A meshing method for reconstructed images, characterized in that, Includes the following steps: Step 1: Perform CT scans, nuclear magnetic resonance, and electron microscopy experiments on the object sample. Obtain geometric images of the surface or internal fractures of the object as required by numerical simulation. Process the scan data into object images with different cross sections through three-dimensional reconstruction. Step 2: Image preprocessing is performed on the tomographic image of the object, including image segmentation, contrast and brightness adjustment, and filtering and noise reduction. Factors that are not conducive to meshing, such as the outer corners of the object and scanning noise disturbances, are removed in advance. Step 3: Perform grayscale analysis on the preprocessed image, divide the object into multiple phases through dual threshold segmentation, including skeletal particles and pore space, and perform three-dimensional stacking and reconstruction on the segmented image slices. Then, perform filtering and noise reduction processing again, and save it as the original XXX.raw format three-dimensional image. Alternatively, filling material can be generated from the three-dimensional stacked and reconstructed core. The stacked or reconstructed digital cores with filling are then cut into smaller cores of lengthX×lengthY×lengthZ, and the core data are stored in the corresponding core files in ascending order of x, y, z. Step 4: Import the core file generated in Step 3 into the MATLAB grid subroutine. The core file is a one-dimensional array arranged in ascending order of x, y, and z. Set the grid pixel size (lengthX, lengthY, lengthZ), recombine the one-dimensional array into a three-dimensional calculation matrix, and judge the correctness of the reconstructed data through three-dimensional matrix slice image visualization processing. Step 5: Number the grid cells and generate node coordinates for the computation matrix. The meshing subroutine numbers nodes sequentially from position (1, 1, 1) in the 3D computation matrix to (lengthX, lengthY, lengthZ). Each mesh cell generates 8 nodes, resulting in a node matrix of size (lengthX+1, lengthY+1, lengthZ+1). It calculates node coordinates as node matrix coordinates multiplied by resolution and sorts the calculated node coordinates according to the cell number. The node coordinates are stored in a matrix Na(lengthX+1, lengthY+1, lengthZ+1), while the mesh cell numbers are stored in cell1, cell2, and cell3 according to different phases. Step 6: The MATLAB mesh subroutine performs surface mesh calculations on mesh elements according to the phase of the skeleton particles, pore spaces, or filling material. The connection format of the six surface meshes of a mesh element is (n1 n2 n3 n4 c1 c2), where ni is the node that makes up the surface mesh and the connection order, and ci is the neighboring mesh element of the surface mesh; the node connection format of the two-dimensional mesh element is (n1 n2 c1 c2); the surface meshes are stored in facem1, facem2, or facem3 according to different phases. Step 7: After the node coordinates and surface mesh calculations are completed, write the mesh file in the standard Fluent mesh format and save the generated mesh file as XXX.msh. In Hypermesh, convert the XXX.msh file into a mesh format that the application software can accept.

2. The image reconstructed meshing method according to claim 1, characterized in that, In step 3, the core filling material is generated according to the following steps: 1) Import the XXX.raw file, pore information file, and skeleton information file, read the image data into a three-dimensional matrix, and form a calculation matrix a(lenx,leny,lenz) to provide basic skeleton and pore data for generating the infill material; 2) Calculate the pore search range, which is a spherical range with a radius of ratio*plr from the pore center outward; and cyclically search for pore pixels with the pore center in mind. 3) Generate filling material in the pores within the pore search range: The pore filling material extends from the pore center to the framework particles until the minimum gap between the filling material and the framework particles is not less than one voxel. 4) After the filling material of one pore is generated, the next pore search begins. After all pores are searched, the filling material saturation is calculated. If the filling material saturation does not reach the preset saturation, the pore search range and the preset pore radius threshold plr0 are adjusted until the preset filling material saturation is met. When the filling material saturation reaches the preset saturation, the filling material generation is complete.

3. The image reconstructed meshing method according to claim 1 or 2, characterized in that: The lengthX, lengthY, and lengthZ are all 50.

Images