Ceramic matrix composite turbine blade disc grid reconstruction method based on CT image slice

Abstract

The invention discloses a ceramic matrix composite turbine blade disc grid reconstruction method based on CT image slices, and belongs to the field of composite structure modeling and computer-aided optimization. According to the method, the accurate grid model is generated through the space curvature point cloud extraction and grid division algorithm. According to the invention, through optimized image preprocessing, binarization processing and an efficient 3D curvature point cloud sampling strategy, in combination with a grid reconstruction algorithm, rapid grid reconstruction is realized. According to the method, through the steps of adding protective layer points, removing degraded tetrahedrons, calculating circumscribed balls and the like in the grid division process, the situation that the grids are unstable is reduced. According to the method, the gray level image is converted into the high-fidelity ceramic matrix composite material turbine bladed disc model, the generated grid model is closer to the real structure of the material, and the prediction precision of subsequent turbine bladed disc finite element analysis is remarkably improved. The method is particularly suitable for a complex ceramic-based composite material structure, and three-dimensional mesh models on the surface and inside of the complex ceramic-based composite material structure can be accurately reconstructed.

Description

Technical Field

[0001] The invention belongs to the field of composite material structure modeling and computer-aided optimization, and relates to a grid model fast reconstruction method for CT images of ceramic-based composite material components. Background Art

[0002] With the development of aerospace and energy fields, the demand for high-performance materials continues to increase. Ceramic matrix composites (CMC) have become the material of choice in these fields due to their excellent high-temperature performance, corrosion resistance and mechanical strength, especially in key components such as turbine blades. Despite this, the complex internal structure of CMC materials and defect detection during the manufacturing process remain technical challenges. Traditional detection methods are difficult to fully and accurately obtain detailed information about the interior of the material, limiting its widespread application in engineering. CT scanning technology can provide high-resolution three-dimensional internal structure images and has been widely used in internal defect detection of composite materials. It can accurately reflect defects such as pores and cracks inside the material.

[0003] CT scanning reconstructs the internal structure image of an object through multi-angle X-ray scanning and computer processing. The grayscale value distribution reflects the density, composition and internal structure of the material. However, CT images only provide grayscale information and are difficult to use directly for engineering calculations such as finite element analysis. For this reason, mesh reconstruction of CT images has become an important research direction. By converting geometric and physical information into finite element format, the calculation accuracy and efficiency can be significantly improved. Especially for CMC turbine blades with complex structures, accurate mesh reconstruction methods can effectively reduce simulation errors and provide reliable data support for optimized design and performance prediction. Summary of the invention

[0004] The purpose of the present invention is to provide a ceramic matrix composite turbine blade disk mesh reconstruction method based on CT image slices, which realizes efficient and accurate turbine blade disk mesh reconstruction through CT image binarization, point cloud extraction based on spatial curvature, and surface and volume meshing algorithms. Compared with traditional methods, the mesh model generated by the present invention is closer to the real structure, can improve the accuracy and reliability of finite element analysis of ceramic matrix composite turbine blade disks, and is helpful for the optimization design and fault prediction of turbine blade disks.

[0005] In order to achieve the above object, the technical solution adopted by the present invention is:

[0006] The present invention discloses a method for reconstructing a grid of a ceramic matrix composite turbine blade disk based on CT image slices, comprising the following steps:

[0007] Step 1: Perform CT scanning on the ceramic matrix composite turbine blade disk to obtain its grayscale image slice;

[0008] Step 2: After preprocessing the grayscale image slice obtained in step 1 to remove noise, the image is binarized through image preprocessing, image enhancement, and threshold segmentation algorithm to obtain a binary image of the turbine blade disk; the binary processing is used to distinguish materials and gaps and highlight the internal structural features of the turbine blade disk;

[0009] Step 3: Based on the binary image of the turbine blade disk obtained in step 2, a 3D edge detection algorithm is used to generate a point cloud, the curvature of the point cloud is estimated, and the point cloud is sampled based on the curvature of the point cloud, and finally the 3D point cloud data of the blade disk is obtained;

[0010] Step 3.1: Load 3D image data: load all binary image files in the specified directory and sort the files according to the image acquisition order; then stack all binary images along the Z axis to form a three-dimensional array to obtain a three-dimensional image;

[0011] Step 3.2: 3D edge detection: First, convert the 3D image into floating point type; then, by setting the upper and lower thresholds, variance and maximum error parameters, apply the Canny edge detection filter to extract the edges in the 3D space and obtain the edge image;

[0012] Step 3.3: Generate 3D point cloud: Convert the edge image into a NumPy array, and find the element with a value of 1 to obtain the coordinates of the edge point, which are converted into actual physical space points according to the resolution; the actual physical space point is created as a point cloud object, and each actual physical space point contains its three-dimensional space coordinates to form a complete 3D point cloud structure;

[0013] Step 3.4: Sampling point cloud based on 3D curvature:

[0014] First, the curvature of each actual physical space point is estimated; the normal vector of each actual physical space point is estimated; then the neighbor point set of each actual physical space point is found by KD tree search, and the covariance matrix of the neighbor point set is calculated; the minimum eigenvalue is obtained by eigenvalue decomposition, and the curvature is expressed by the ratio of the minimum eigenvalue to the sum of eigenvalues; the formula is as follows:

[0015]

[0016] Where: min is the minimum eigenvalue of the covariance matrix, λ i is the sum of all eigenvalues;

[0017] The covariance matrix of the neighborhood point set is:

[0018] Where p j is a point in the neighborhood point set, As the centroid of the neighborhood point set, select point pi The neighborhood point set is N(p i ), the point set is obtained by KD tree search; for the covariance matrix C i Perform eigenvalue decomposition to obtain eigenvalue λ;

[0019] Then, the high curvature point cloud and low curvature point cloud areas are divided according to the curvature size, and point cloud sampling is performed on the two areas respectively to obtain the 3D point cloud data of the blade disk;

[0020] Step 4: Use the 3D point cloud data obtained in step 3 as mesh nodes, perform Delaunay triangulation by initializing and expanding the data, establish the connection relationship between tetrahedrons, and reconstruct the surface mesh by calculating the circumscribed sphere and eliminating tetrahedrons that do not conform to the Delaunay property, and finally generate the surface triangular mesh of the turbine blade.

[0021] Step 4.1: Point cloud data reading and initialization: Store the 3D point cloud data obtained in step 3 as a vertex matrix V = [v1, v2, ..., v n ], where each vertex v i Represents a 3D coordinate (x i ,y i ,z i );

[0022] Step 4.2: Add protection layer points: Calculate the minimum and maximum values ​​of the 3D point cloud data, expand a certain proportion to obtain new boundary coordinates, and evenly generate protection layer points on the six faces to obtain protection layer point cloud data; set the number of protection layer points to N shield , the expansion step is Delta; the maximum and minimum values ​​of the 3D point cloud data are:

[0023] coords max =max(V)

[0024] coords min =min(V)

[0025] Extension step:

[0026] Delta=(coords max -coords min )×offset

[0027] Where: offset is a scaling factor used to control the expansion of the point cloud data boundary. Through this scaling factor, an extended bounding box is generated to surround the original point cloud data;

[0028] Step 4.3: Delaunay triangulation: triangulate the 3D point cloud data and the protection layer point cloud data, and use the Delaunay algorithm to generate tetrahedral units; to ensure the correctness of the triangulation, remove the degenerate tetrahedrons through volume judgment;

[0029] Determination of the volume of a tetrahedron:

[0030]

[0031] Where: (z1, y1, z1), (z2, y2, z2), (z3, y3, z3) and (z4, y4, z4) are the coordinates of the four vertices; when the volume is less than the threshold, the tetrahedron is removed from the segmentation result;

[0032] Step 4.4: Calculation of connection relationship of tetrahedrons: Calculate the connection relationship between the generated tetrahedral units. The connection between each tetrahedron and its neighboring tetrahedrons is determined by the shared triangular faces. Establish the connection matrix between tetrahedrons for subsequent surface extraction and manifold construction.

[0033] Step 4.5: Calculation of the circumscribed sphere: For each generated tetrahedral unit, calculate the center and radius of its circumscribed sphere;

[0034] The four vertices of the tetrahedron are v1 = (x1, y1, z1), v2 = (x2, y2, z2), v3 = (x3, y3, z3) and v4 = (x4, y4, z4), and the circumcenter C c =(x c ,y c ,z c ) is calculated using the determinant formula:

[0035]

[0036] Where: D1, D2, D3 are vectors calculated based on the positions of the four vertices, and their specific definitions are as follows:

[0037]

[0038] det(V 23 )、det(V 13 )、det(V 12 )、det(V 41 ) represent the determinants of the following sub-matrices respectively:

[0039] det(V 23 ) is a matrix consisting of the coordinates of vertices v2, v3, and v1, expressed as:

[0040]

[0041] det(V 13 ) is a matrix consisting of the coordinates of vertices v1, v3, and v2, expressed as:

[0042]

[0043] det(V 12 ) is a matrix consisting of the coordinates of vertices v1, v2, and v3, expressed as:

[0044]

[0045] det(V 41 ) is a matrix consisting of the coordinates of vertices v4, v1, v2, and v3, i.e., the volume of the tetrahedron, expressed as:

[0046]

[0047] Calculate the circumscribed sphere center C c The coordinates of , and further calculate the radius of the circumscribed sphere:

[0048] r i =||C c -v i ||

[0049] Where: v i is a vertex in the tetrahedron, i = 1, 2, 3, 4, ||·|| represents the Euclidean distance, that is, the straight-line distance from the center of the sphere to the vertex;

[0050] Step 4.6: Determine whether the tetrahedron conforms to the Delaunay property by calculating the circumscribed sphere of the tetrahedron;

[0051] The intersection factor is used to determine whether two tetrahedrons intersect, and tetrahedrons that do not meet the conditions are retained or eliminated; the circumscribed sphere centers of two adjacent tetrahedrons are C d and C e , the radii are r d and r e The distance between the two sphere centers is d de =||C d -C e ||; then the cross factor If de It is expressed as:

[0052]

[0053] When de When >1, the two circumscribed spheres are separated, and there is no intersection between the two tetrahedrons, which conforms to the Delaunay property, and these tetrahedrons are retained;

[0054] Whende =1, the two circumscribed spheres are exactly tangent, which means that the two tetrahedrons may touch but not intersect, which conforms to the Delaunay property and retains these tetrahedrons;

[0055] When de When <1, the two circumscribed spheres intersect, which means that the two tetrahedrons may intersect, which does not conform to the Delaunay property, and these tetrahedrons are eliminated.

[0056] Step 4.7: Surface triangle mesh extraction and manifold construction:

[0057] Traverse the tetrahedron and extract its faces to build a triangle list, then remove the internally shared duplicate triangles, retain only the surface triangles, and finally generate a surface mesh;

[0058] Step 5: Based on the surface mesh obtained in step 4, a closed surface mesh is generated through hole filling and mesh optimization algorithm;

[0059] Step 5.1: Hole identification and filling: The edges not shared by two triangles are free edges. The triangles connected to the free edges are deleted to corrode the holes, and the holes are filled by the fill_holes function. After the filling is completed, the holes are detected and corroded again. If there are still holes, they are filled by the hole filling iterative algorithm.

[0060] Step 5.2: Optimize and smooth the mesh filled in step 5.1 by using the Laplacian operator to obtain a closed surface mesh;

[0061] Step 6: Using the closed surface mesh obtained in step 5 as input, insert new nodes inside the mesh, generate the tetrahedral mesh of the turbine blade disk through the frontier advancing algorithm, and complete the mesh reconstruction of the turbine blade disk.

[0062] It also includes step seven: performing a quality assessment on the tetrahedral mesh generated in step six to ensure that the geometry of the mesh meets the requirements of finite element analysis, performing simulation prediction of the ceramic matrix composite turbine blade disk based on the mesh of the turbine blade disk selected through the quality assessment, and solving engineering problems related to the application of ceramic matrix composite turbine blade disks.

[0063] Engineering issues related to the application of ceramic matrix composite turbine blades include thermal stress distribution and fatigue life prediction in high temperature environments, micro crack propagation and damage tolerance evaluation of materials, and structural stability and dynamic response analysis under complex loads.

[0064] Step 1 is implemented as follows:

[0065] Step 1.1: Preparation: Fix the CMC turbine blade disk on the sample stage of the CT scanner to ensure that it is stable during the scanning process;

[0066] Step 1.2: Parameter setting: According to the material properties of the turbine blade and the required resolution, set the voltage, power, exposure time and scanning time during CT scanning;

[0067] Step 1.3: Scanning process: Start the CT scanner and perform continuous tomographic scanning to obtain grayscale image slice data of the turbine blade disk layer by layer.

[0068] Step 2 is implemented as follows:

[0069] Step 2.1: Image preprocessing: Use the median filter denoising algorithm to preprocess the original CT image to remove noise and artifacts in the image; median filter formula:

[0070] I filtered (x,y)=median{I(x+i,y+j)|(i,j)∈W}

[0071] Where: I filtered (x, y) is the pixel value of the filtered image, I(x+i, y+j) is the pixel value of the original image at (x+i, y+j), W is the filter window, and i and j are the offsets of the pixel positions;

[0072] Step 2.2: Image enhancement: Apply histogram equalization technology to enhance the image contrast, making the boundaries between materials and gaps more obvious; histogram equalization formula:

[0073]

[0074] Where: I equalized (x, y) is the pixel value of the image after equalization, L is the grayscale level of the image, MN is the total number of pixels in the image, h(i) is the number of pixels with grayscale value i; I(x, y) is the pixel value of the original image;

[0075] Step 2.3: Threshold segmentation: Use Otsu adaptive threshold segmentation algorithm to binarize the preprocessed image and convert the image into a binary image;

[0076]

[0077] Where: is the inter-class variance, t is the threshold, ω1(t) and ω2(t) are the pixel probabilities of the foreground and background respectively, μ1(t) and μ2(t) are the average grayscale values ​​of the foreground and background respectively; select t that maximizes the inter-class variance as the optimal threshold.

[0078] Step 7 is implemented as follows:

[0079] Step 7.2: Calculate the aspect ratio A:

[0080]

[0081] Where: E is the side length of the tetrahedron, and the calculation formula for the side length E is: Where (x a ,y a ,z a ) and (x b ,y b ,z b ) are the coordinates of the vertices at both ends of the connecting edge E, max(E) and min(E) are the longest and shortest edges respectively; the tetrahedral mesh with A<2 is a good mesh, the tetrahedral mesh with 2≤As<3 is a better mesh, and the tetrahedral mesh with A>3 is a poor mesh.

[0082] The core of the Laplacian optimization described in step 5.2 is the encoding and decoding of the local detail features of the mesh model. In the encoding stage, the Euclidean space coordinates of the mesh vertices are converted to Laplacian coordinates, which can accurately capture the local detail features of the mesh. Therefore, the Laplacian mesh deformation algorithm can retain these details when processing the mesh. In the decoding stage, the Laplacian coordinates are reversed back to the Euclidean space coordinates by solving the linear system, thereby restoring and optimizing the original mesh. w The Laplacian coordinates are defined as:

[0083] P w =v w +λ·∑ j∈N(w) v j -v w

[0084] Where: v w is the original position of vertex w, P w is the new position of vertex w after smoothing, N(w) represents the set of vertices adjacent to vertex w, λ is the relaxation factor, which controls the smoothing strength, generally 0<λ<1, v j is the position of vertex j, j is the vertex adjacent to w.

[0085] Preferably, the meshing software is Gmsh Python API or Hypermesh.

[0086] The implementation of the method requires ceramic matrix composite turbine blades, industrial CT testing equipment, and software tools for data processing and code writing.

[0087] The ceramic-based composite material used in the turbine blade disk is widely used in the aerospace field, such as aircraft wings, fuselages, satellite antennas and their supporting structures, launch vehicle casings, engine casings, etc.; as well as the automotive industry, such as automobile bodies, load-bearing components, transmission shafts, engine mounts, etc.

[0088] Beneficial effects:

[0089] 1. The ceramic matrix composite turbine blade disk mesh reconstruction method based on CT images disclosed in the present invention generates an accurate mesh model through spatial curvature point cloud extraction and meshing algorithm, thereby improving the accuracy of finite element analysis.

[0090] 2. The ceramic-based composite turbine blade disk mesh reconstruction method based on CT image slicing disclosed in the present invention realizes rapid mesh reconstruction through optimized image preprocessing, binarization processing and efficient 3D curvature point cloud sampling strategy, combined with a mesh reconstruction algorithm, thereby greatly shortening the finite element calculation time.

[0091] 3. The ceramic-based composite turbine blade disk mesh reconstruction method based on CT image slicing disclosed in the present invention can reduce mesh instability, improve the stability of the generated mesh, and ensure the reliability of the turbine blade disk mesh model used for simulation analysis by adding protection layer points, eliminating degenerate tetrahedrons, and calculating circumscribed spheres during the mesh division process.

[0092] 4. Compared with the existing traditional modeling methods, the ceramic-based composite turbine blade disk mesh reconstruction method based on CT image slicing disclosed in the present invention converts the grayscale image into a high-fidelity ceramic-based composite turbine blade disk model. The generated mesh model is closer to the actual structure of the material, which significantly improves the prediction accuracy of subsequent turbine blade disk finite element analysis and reduces the error between the calculation results and the experimental data.

[0093] 5. The ceramic-based composite turbine blade disk mesh reconstruction method based on CT image slices disclosed in the present invention is particularly suitable for complex ceramic-based composite structures. It can accurately reconstruct the three-dimensional mesh model of its surface and interior, provide more realistic and high-fidelity engineering calculation input data, and help improve the accuracy and efficiency of performance optimization and design of ceramic-based composite turbine blade disks. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] Figure 1 It is a flow chart of the method for reconstructing the grid of a ceramic matrix composite turbine blade disk based on CT image slices disclosed in the present invention;

[0095] Figure 2 The images before and after the CT images are binarized based on the threshold segmentation algorithm of the present invention (a is the grayscale image of the turbine blade disk after CT testing of the present invention; b is the binary image of the turbine blade disk after threshold segmentation);

[0096] Figure 3 Schematic diagram of the 3D point cloud of a turbine blade disk obtained based on 3D edge detection and 3D curvature point cloud sampling algorithm (a is the 3D point cloud structure of the turbine blade disk; b is the 3D point cloud data after 3D point cloud curvature sampling);

[0097] Figure 4 Schematic diagram of reconstructing the mesh model based on the 3D point cloud data of the turbine blade (a is the surface mesh of the turbine blade; b is the body mesh of the turbine blade). DETAILED DESCRIPTION

[0098] In order to better illustrate the purpose and advantages of the present invention, the invention is further described below with reference to the accompanying drawings and examples.

[0099] Embodiment 1:

[0100] like Figure 1 As shown, the ceramic matrix composite turbine blade disk mesh reconstruction method based on CT image slices disclosed in this embodiment has the following specific implementation steps:

[0101] Step 1: Perform CT scanning on the ceramic matrix composite turbine blade disk to obtain its grayscale image slices such as Figure 2 (a) shows a grayscale image of the turbine blade disk after CT testing, including 488 CT images.

[0102] Step 1.1: Preparation: Fix the CMC turbine blade disk on the sample stage of the CT scanner to ensure that it is stable and motionless during the scanning process.

[0103] Step 1.2: Parameter setting: According to the material properties of the turbine blade disk, the resolution is set to 34 μm, the voltage during CT scanning is set to 200 kV, the power is set to 15 W, the exposure time is set to 500 ms, and the scanning time is set to 1 h;

[0104] Step 1.3: Scanning process: Start the CT scanner, perform continuous tomographic scanning, and obtain grayscale image slice data of the turbine blade disk layer by layer. A total of 488 CT images are obtained, and the image size is 2060×1274;

[0105] Step 2: After preprocessing the grayscale image slice obtained in step 1 to remove noise, the image is binarized through image preprocessing, image enhancement, and threshold segmentation algorithm to obtain a binary image of the turbine blade disk; the binarization process can distinguish materials and gaps and highlight the internal structural features of the turbine blade disk;

[0106] Step 2 is implemented as follows:

[0107] Step 2.1: Image preprocessing: First, the median filter algorithm is used to perform denoising on the acquired original CT grayscale image slices. The median filter process is to move a window. In this embodiment, a 3×3 median filter window is used to denoise the CT image. Traverse each pixel of the image, sort the pixel values ​​in the window by size, and use the middle value as the new pixel value of the pixel. This can effectively remove random noise in the image while retaining edge details to obtain a CT image after image preprocessing. The median filter formula is shown in formula (1):

[0108] I fltered (x,y)=median{I(x+i,y+j)|(i,j)∈W} (1)

[0109] Where: I filtered (x, y) is the pixel value of the filtered image, I(x+i, y+j) is the pixel value of the original image at (x+i, y+j), W is the filter window, and i and j are the offsets of the pixel positions;

[0110] Step 2.2: Image enhancement: For the preprocessed CT image obtained in step 2.1, the image contrast is further improved, materials and voids are distinguished, and the image is enhanced using histogram equalization technology. By calculating the cumulative histogram of the image, the grayscale distribution of the image is stretched to the entire grayscale range, thereby enhancing the image contrast and obtaining an enhanced CT image. The histogram equalization formula is shown in formula (2):

[0111]

[0112] Where: I equalized (x, y) is the pixel value of the image after equalization, L is the grayscale level of the image, in this example L = 256, MN is the total number of pixels of the image, in this example MN = 1274 × 2060, h(i) is the number of pixels with grayscale value i;

[0113] Step 2.3: Threshold segmentation: For the enhanced CT image obtained in 2.2, an adaptive threshold segmentation algorithm is used to perform binarization on the image. The Otsu method automatically selects the optimal threshold by maximizing the inter-class variance, and classifies the part of the image with a grayscale value higher than the threshold as the foreground (i.e., material), and the part with a grayscale value lower than the threshold as the background (i.e., gap). The inter-class variance formula of the Otsu method is shown in formula (3):

[0114]

[0115] Where: is the inter-class variance, t is the threshold, ω1(t) and ω2(t) are the pixel probabilities of the foreground and background, μ1(t) and μ2(t) are the average grayscale values ​​of the foreground and background. The t that maximizes the inter-class variance is selected as the optimal threshold. In this example, the threshold t = 133. The image after threshold segmentation is as follows: Figure 2 (b)

[0116] Step 3: According to step 2, 488 binary images of turbine blades are obtained, and a point cloud is generated using a 3D edge detection algorithm. The curvature of the point cloud is estimated, and the point cloud is sampled based on the curvature of the point cloud, and finally the 3D point cloud data of the blade is obtained;

[0117] Step 3.1: Load 3D image data: Load 488 binary image files in the specified directory and sort the files according to the image acquisition order. Then stack all binary images along the Z axis to form a three-dimensional array to obtain a three-dimensional image. This process ensures the spatial relationship and sequence continuity of the image, laying the foundation for subsequent processing based on spatial analysis;

[0118] Step 3.2: 3D edge detection: First, convert the three-dimensional image into a floating point type to ensure the accuracy of the processing. On this basis, the Canny edge detection algorithm is used to perform edge detection on the three-dimensional image to obtain an edge image. The parameters in the Canny edge detection algorithm are set as follows: Upper Threshold = 1, used to control the high threshold of edge detection; set Lower Threshold = 0, used to control the low threshold of edge detection; set Variance = 10, used to control the Gaussian blur degree of the filter; set Maximum Error = 0.5, used to control the allowable error in the filter calculation;

[0119] Step 3.3: Generate 3D point cloud: Convert the edge image into a NumPy array and find the element with value 1 to obtain the coordinates of the edge points. After converting these edge points into coordinates in the actual physical space, use the Open3D library to generate a point cloud object. Each point cloud point contains its x, y, and z coordinates in three-dimensional space. These coordinates are formatted and saved in the point cloud structure to form a complete three-dimensional point cloud structure, such as Figure 3 (a) shows the three-dimensional point cloud structure of the turbine blade disk;

[0120] Step 3.4: Sampling point cloud based on 3D curvature:

[0121] In this step, the curvature of each actual physical space point is first estimated. Using the Open3D library, the normal vector of each point is first estimated; then the KD tree search algorithm is used to find the neighbor point set of each actual physical space point and calculate the covariance matrix of the neighbor point set. The minimum eigenvalue is obtained through eigenvalue decomposition, and the ratio of it to the sum of all eigenvalues ​​is used to express the curvature. The specific formula is as follows:

[0122]

[0123] Where: min is the minimum eigenvalue of the covariance matrix, λ i is the sum of all eigenvalues;

[0124] The covariance matrix of the neighborhood point set is:

[0125] Where p j is a point in the neighborhood point set, As the centroid of the neighborhood point set, select point p i The neighborhood point set is N(p i ), which is found by KD tree search. i Perform eigenvalue decomposition to obtain eigenvalue λ;

[0126] The curvature value of each point is calculated by the above formula, and then these points are classified. The high curvature threshold curvature_threshold is set to control the sparseness of the point cloud area sampling. In this example, curvature_threshold = 75, that is, points above this threshold are considered high curvature points. The high curvature point cloud and the low curvature point cloud are downsampled by voxels, respectively. The high curvature point voxel high_curvature_voxel_size = 3, and the low curvature point voxel low_curvature_voxel_size = 6. The 3D point cloud data obtained after 3D curvature sampling is as follows Figure 3 (b)

[0127] Step 4: Use the 3D point cloud data obtained in step 3 as mesh nodes, perform Delaunay triangulation by initializing and expanding the data, establish the connection relationship between tetrahedrons, and reconstruct the surface mesh by calculating the circumscribed sphere and eliminating tetrahedrons that do not conform to the Delaunay property, and finally generate the surface triangular mesh of the turbine blade disk;

[0128] Step 4.1: Point cloud data reading and initialization: Store the 3D point cloud data obtained in step 3 as a vertex matrix V = [v1, v2, ..., v n ], where each vertex v iRepresents a 3D coordinate (x i ,y i ,z i );

[0129] Step 4.2: Add protection layer points: To prevent boundary effects, it is necessary to add protection layer points outside the original point cloud data. First, calculate the minimum and maximum values ​​of the 3D point cloud data, expand a certain proportion to obtain new boundary coordinates, and evenly generate protection layer points on the six faces to obtain protection layer point cloud data. Set the number of protection layer points to N. shield In this example, Nshield = 5 and the expansion step is Delta. The maximum and minimum values ​​of the original point cloud are calculated as follows:

[0130] coords max =max(V) (5)

[0131] coords min =min(V) (6)

[0132] Extension step:

[0133] Delta=(coords max -coords min )×offset (7)

[0134] Where: offset is a scaling factor used to control the expansion of the point cloud data boundary. Through this scaling factor, an extended bounding box can be generated to surround the original point cloud data. In this example, offset = 1. This expansion can avoid mesh instability or inaccuracy caused by the influence of boundary points in the subsequent triangulation process;

[0135] Step 4.3: Delaunay triangulation: triangulate the 3D point cloud data and the protection layer point cloud data, and use the Delaunay algorithm to generate tetrahedral units; to ensure the correctness of the triangulation, remove the degenerate tetrahedrons through volume judgment;

[0136] Tetrahedron volume determination:

[0137]

[0138] Where: (x1, y1, z1), (x2, y2, z2), (x3, y3, z3) and (x4, y4, z4) are the coordinates of the four vertices; the volume Volume is compared with the threshold T to determine whether to remove the tetrahedron. In this example, T = 1e-15, if the Volume is less than the threshold T, the tetrahedron is removed from the segmentation result;

[0139] Step 4.4: Calculation of connection relationship of tetrahedrons: Calculate the connection relationship between the generated tetrahedral units. The connection between each tetrahedron and its neighboring tetrahedrons is determined by the shared triangular faces. Establish a connection matrix between tetrahedrons, record the shared face information of each tetrahedron and its neighboring tetrahedrons, and provide data support for subsequent surface extraction and manifold construction.

[0140] Step 4.5: Calculation of circumscribed sphere: For each tetrahedral unit generated, calculate the center and radius of its circumscribed sphere. By calculating the circumscribed sphere of a tetrahedron, you can determine whether the tetrahedron conforms to the Delaunay property. An important property of Delaunay triangulation is that for any tetrahedron, its circumscribed sphere should not contain any other points. If a point is within the circumscribed sphere, the tetrahedron violates the Delaunay property and needs to be adjusted;

[0141] The four vertices of the tetrahedron are v1 = (x1, y1, z1), v2 = (x2, y2, z2), v3 = (x3, y3, z3) and v4 = (x4, y4, z4), and the circumcenter C c =(x c ,y c ,z c ) is calculated using the determinant formula:

[0142]

[0143] Where: D1, D2, D3 are vectors calculated based on the positions of the four vertices, and their specific definitions are as follows:

[0144]

[0145] det(v 23 )、det(v 13 )、det(v 12 )、det(v 41 ) represent the determinants of the following sub-matrices respectively:

[0146] det(V 23 ) is a matrix consisting of the coordinates of vertices v2, v3, and v1, expressed as:

[0147]

[0148] det(V 13 ) is a matrix consisting of the coordinates of vertices v1, v3, and v2, expressed as:

[0149]

[0150] det(V 12) is a matrix consisting of the coordinates of vertices v1, v2, and v3, expressed as:

[0151]

[0152] det(V 41 ) is a matrix consisting of the coordinates of vertices v4, v1, v2, and v3, i.e., the volume of the tetrahedron, expressed as:

[0153]

[0154] Calculate the circumscribed sphere center C c The coordinates of the sphere can be calculated by formula (10):

[0155] r i =||C c -v i || (10)

[0156] Where: v i is a vertex in the tetrahedron, i = 1, 2, 3, 4, ||·|| represents the Euclidean distance, that is, the straight-line distance from the center of the sphere to the vertex;

[0157] Step 4.6: Determine whether the tetrahedron meets the Delaunay property by calculating the circumscribed sphere of the tetrahedron; calculate the intersection factor between each tetrahedron to determine whether two tetrahedrons intersect, and retain or eliminate tetrahedrons that do not meet the conditions. For two adjacent tetrahedrons, their circumscribed sphere centers are C d and C e , the radii are r d and r e The distance between the two sphere centers is d de =||C d -C e ||, then the cross factor If de It can be expressed as:

[0158]

[0159] When de When >1, the two circumscribed spheres are separated, and there is no intersection between the two tetrahedrons, which conforms to the Delaunay property, and these tetrahedrons are retained;

[0160] When de =1, the two circumscribed spheres are exactly tangent, which means that the two tetrahedrons may touch but not intersect, which conforms to the Delaunay property and retains these tetrahedrons;

[0161] When deWhen <1, the two circumscribed spheres intersect, which means that the two tetrahedrons may intersect, which does not conform to the Delaunay property, so these tetrahedrons are eliminated;

[0162] Step 4.7: Surface triangle mesh extraction and manifold construction:

[0163] By traversing the tetrahedron connectivity and based on the aforementioned cross factor and geometric consistency, surface triangles are extracted from the three-dimensional tetrahedral mesh, and it is ensured that the extracted surface constitutes a good manifold. First, by traversing each tetrahedron, all its faces are extracted and a triangle list is constructed. Next, duplicate triangles shared within the mesh are removed, leaving only surface triangles. Then, the adjacency relationship between triangles is checked to ensure that each edge is shared by exactly two triangles, thereby ensuring that the surface constitutes a manifold. If there is damage or discontinuity, it is repaired to ultimately generate a complete and seamless surface mesh;

[0164] Step 5: Based on the surface mesh obtained in step 4, a closed surface mesh is generated through hole filling and mesh optimization algorithm;

[0165] Step 5.1: Hole identification and filling: The edges that are not shared by two triangles are free edges. Delete the triangles connected to the free edges to corrode the holes, and fill the holes with the fill_holes function. After filling, perform hole detection and corrosion again. If holes still exist, fill them with the hole filling iterative algorithm. The implementation process is as follows: check whether there are still unclosed boundaries in the filled mesh; if there are still holes, find the triangles connected to the hole boundaries and remove them from the mesh; the above process will be repeated for a specified number of iterations, and each iteration will try to fill new holes and remove possible bad units. In this example, iterations = 3; after all iterations are completed, perform the last filling operation to ensure that there are no remaining holes. The turbine blade surface mesh established is as follows: Figure 4 (a)

[0166] Step 5.2: Use the Laplacian operator to optimize and smooth the mesh filled in step 5.1. The key to Laplacian optimization is the encoding and decoding of local details of the mesh model. During the encoding process, the Euclidean coordinates of the mesh vertices are converted to Laplacian coordinates to accurately capture the local detail features, thereby retaining these details when the mesh is deformed. In the decoding stage, the Laplacian coordinates are converted back to Euclidean coordinates by solving the linear system to restore and optimize the mesh. w The Laplacian coordinates are defined as:

[0167] P w =vw +λ·∑ j∈N(w) v j -v w (12)

[0168] Where: v w is the original position of vertex w, P w is the new position of vertex w after smoothing, N(w) represents the set of vertices adjacent to vertex w, λ is the relaxation factor, which controls the smoothing strength, generally 0<λ<1, in this example λ=0.4, v j is the position of vertex j, j is the vertex adjacent to w;

[0169] Step 6: Using the turbine blade surface triangular mesh obtained in step 5 as input, insert new nodes inside the mesh, and generate the turbine blade tetrahedron mesh through the frontier advancing algorithm. The resulting turbine blade body mesh is as follows: Figure 4 (b) as shown.

[0170] Step 7: Evaluate the quality of the generated tetrahedral mesh to ensure that the mesh geometry meets the requirements of finite element analysis. The specific method is to calculate the aspect ratio and evaluate the mesh quality based on the calculation results;

[0171] Step 7 is implemented as follows:

[0172] Step 7.1: Calculate the aspect ratio: The aspect ratio calculation evaluates the shape quality of the tetrahedron. The ideal tetrahedron is close to a regular tetrahedron, with little difference in side length. A too large aspect ratio will cause the shape to deviate, become flat or elongated, and thus reduce the accuracy and stability of the finite element analysis. The aspect ratio calculation formula is shown in formula (13):

[0173]

[0174] Where: E is the side length of the tetrahedron, and the calculation formula for the side length E is: , where (x a ,y a ,z a ) and (x b ,y b ,z b ) are the coordinates of the vertices at both ends of the connecting edge E, max(E) and min(E) are the longest and shortest edges, respectively. In this example, the calculation shows that the proportion of tetrahedrons with A<2 is 95.2%, the proportion of tetrahedrons with 2≤A<3 is 4.8%, and the proportion of tetrahedrons with A>3 is 0%, which indicates that the overall quality of the generated mesh is good.

[0175] Step 8: Based on the mesh of the turbine blade selected by the quality assessment in step 7, perform simulation prediction of the ceramic matrix composite turbine blade to solve engineering problems related to the application of ceramic matrix composite turbine blades.

[0176] Engineering issues related to the application of ceramic matrix composite turbine blades include thermal stress distribution and fatigue life prediction in high temperature environments, micro crack propagation and damage tolerance evaluation of materials, and structural stability and dynamic response analysis under complex loads.

[0177] Based on the grid of the turbine blade selected by the quality assessment in step 7, the thermal stress field of the turbine blade under high temperature and complex loads is simulated, and the thermal stress concentration area of ​​the key parts of the blade is calculated. In the stress concentration area, the extension path and speed of the micro cracks are further evaluated, and the remaining service life of the blade under the set working conditions is predicted in combination with the fatigue life model. Through the simulation results, the material composition and structural design of the blade are optimized to improve its fatigue life under extreme conditions and meet the requirements of long-term stable operation of aircraft engines.

[0178] The specific description above further illustrates the purpose, technical solutions and beneficial effects of the invention in detail. It should be understood that the above is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A ceramic matrix composite turbine blade disk grid reconstruction method based on CT image slices, characterized by: The following steps are included: Step 1: Perform CT scanning on the ceramic matrix composite turbine blade disk to obtain its grayscale image slice; Step 2: After preprocessing the grayscale image slice obtained in step 1 to remove noise, the image is binarized through image preprocessing, image enhancement, and threshold segmentation algorithm to obtain a binary image of the turbine blade disk; the binary processing is used to distinguish materials and gaps and highlight the internal structural features of the turbine blade disk; Step 3: Based on the binary image of the turbine blade disk obtained in step 2, a 3D edge detection algorithm is used to generate a point cloud, the curvature of the point cloud is estimated, and the point cloud is sampled based on the curvature of the point cloud to obtain the 3D point cloud data of the blade disk; Step 3.1: Load 3D image data: load all binary image files in the specified directory and sort the files according to the image acquisition order; then stack all binary images along the Z axis to form a three-dimensional array to obtain a three-dimensional image; Step 3.2: 3D edge detection: First, convert the 3D image into floating point type; then, by setting the upper and lower thresholds, variance and maximum error parameters, apply the Canny edge detection filter to extract the edges in the 3D space and obtain the edge image; Step 3.3: Generate 3D point cloud: Convert the edge image into a NumPy array, and find the element with a value of 1 to obtain the coordinates of the edge point, which are converted into actual physical space points according to the resolution; the actual physical space point is created as a point cloud object, and each actual physical space point contains its three-dimensional space coordinates to form a complete 3D point cloud structure; Step 3.4: Sampling point cloud based on 3D curvature: First, the curvature of each actual physical space point is estimated; the normal vector of each actual physical space point is estimated; Then use KD tree search to find the neighbor point set of each actual physical space point, and calculate the covariance matrix of the neighbor point set; obtain the minimum eigenvalue through eigenvalue decomposition, and use the ratio of the minimum eigenvalue to the sum of eigenvalues ​​to express the curvature; the formula is as follows: Where: min is the minimum eigenvalue of the covariance matrix, λ i is the sum of all eigenvalues; The covariance matrix of the neighborhood point set is: Where p j is a point in the neighborhood point set, As the centroid of the neighborhood point set, select point p i The neighborhood point set is N(p i ), the point set is obtained by KD tree search; The covariance matrix C i Perform eigenvalue decomposition to obtain eigenvalue λ; The point cloud is divided into high curvature point cloud and low curvature point cloud areas according to the curvature size, and the point cloud sampling is performed on the two areas respectively to obtain the 3D point cloud data of the blade disk; Step 4: Use the 3D point cloud data obtained in step 3 as mesh nodes, perform Delaunay triangulation by initializing and expanding the data, establish the connection relationship between tetrahedrons, and reconstruct the surface mesh by calculating the circumscribed sphere and eliminating tetrahedrons that do not conform to the Delaunay property, thus generating the surface triangular mesh of the turbine blade disk; Step 4.1: Point cloud data reading and initialization: Store the 3D point cloud data obtained in step 3 as a vertex matrix V = [v1, v2, ..., v n ], where each vertex v i Represents a 3D coordinate (x i ,y i ,z i ); Step 4.2: Add protection layer points: Calculate the minimum and maximum values ​​of the 3D point cloud data, expand a certain proportion to obtain new boundary coordinates, and evenly generate protection layer points on the six faces to obtain protection layer point cloud data; set the number of protection layer points to N shield , the expansion step is Delta; the maximum and minimum values ​​of the 3D point cloud data are: coords max =max(V) coords min =min(V) Extension step: Delta=(coords max -coords min )×offset Where: offset is a scaling factor used to control the expansion of the point cloud data boundary. Through this scaling factor, an extended bounding box is generated to surround the original point cloud data; Step 4.3: Delaunay triangulation: triangulate the 3D point cloud data and the protection layer point cloud data, and use the Delaunay algorithm to generate tetrahedral units; to ensure the correctness of the triangulation, remove the degenerate tetrahedrons through volume judgment; Determination of the volume of a tetrahedron: Where: (x1, y1, z1), (x2, y2, z2), (x3, y3, z3) and (x4, y4, z4) are the coordinates of the four vertices; when the volume is less than the threshold, the tetrahedron is removed from the segmentation result; Step 4.4: Calculation of connection relationship of tetrahedrons: Calculate the connection relationship between the generated tetrahedral units. The connection between each tetrahedron and its neighboring tetrahedrons is determined by the shared triangular faces. Establish the connection matrix between tetrahedrons for subsequent surface extraction and manifold construction. Step 4.5: Calculation of the circumscribed sphere: For each generated tetrahedral unit, calculate the center and radius of its circumscribed sphere; The four vertices of the tetrahedron are v1 = (x1, y1, z1), v2 = (x2, y2, z2), v3 = (x3, y3, z3) and v4 = (x4, y4, z4), and the circumcenter C c =(x c ,y c ,z c ) is calculated using the determinant formula: Where: D1, D2, D3 are vectors calculated based on the positions of the four vertices, and their specific definitions are as follows: det(V 23 )、det(V 13 )、det(V 12 )、det(V 41 ) represent the determinants of the following sub-matrices respectively: det(V 23 ) is a matrix consisting of the coordinates of vertices v2, v3, and v1, expressed as: det(V 13 ) is a matrix consisting of the coordinates of vertices v1, v3, and v2, expressed as: det(V 12 ) is a matrix consisting of the coordinates of vertices v1, v2, and v3, expressed as: det(V 41 ) is a matrix consisting of the coordinates of vertices v4, v1, v2, and v3, i.e., the volume of the tetrahedron, expressed as: Calculate the circumscribed sphere center C c The coordinates of , and further calculate the radius of the circumscribed sphere: r i =||C c -v i || Where: v i is a vertex in the tetrahedron, i = 1, 2, 3, 4, ||·|| represents the Euclidean distance, that is, the straight-line distance from the center of the sphere to the vertex; Step 4.6: Determine whether the tetrahedron conforms to the Delaunay property by calculating the circumscribed sphere of the tetrahedron: The intersection factor is used to determine whether two tetrahedrons intersect, and tetrahedrons that do not meet the conditions are retained or eliminated; the circumscribed sphere centers of two adjacent tetrahedrons are C d and C e , the radii are r d and r e The distance between the two sphere centers is d de =||C d -C e ||; then the cross factor If de It is expressed as: When If de When >1, the two circumscribed spheres are separated, and there is no intersection between the two tetrahedrons, which conforms to the Delaunay property, and these tetrahedrons are retained; When If de =1, the two circumscribed spheres are exactly tangent, which means that the two tetrahedrons may touch but not intersect, which conforms to the Delaunay property and retains these tetrahedrons; When If de When <1, the two circumscribed spheres intersect, indicating that the two tetrahedrons may intersect, which does not conform to the Delaunay property, and these tetrahedrons are eliminated; Step 4.7: Surface triangle mesh extraction and manifold construction: Traverse the tetrahedron and extract its faces to build a triangle list, then remove the internally shared duplicate triangles, retain only the surface triangles, and generate a surface mesh; Step 5: Based on the surface mesh obtained in step 4, a closed surface mesh is generated through hole filling and mesh optimization algorithm; Step 5.1: Hole identification and filling: The edges not shared by two triangles are free edges. The triangles connected to the free edges are deleted to corrode the holes, and the holes are filled by the fill_holes function. After the filling is completed, the holes are detected and corroded again. If there are still holes, they are filled by the hole filling iterative algorithm. Step 5.2: Optimize and smooth the mesh filled in step 5.1 by using the Laplacian operator to obtain a closed surface mesh; Step 6: Using the closed surface mesh obtained in step 5 as input, insert new nodes inside the mesh, generate the tetrahedral mesh of the turbine blade disk through the frontier advancing algorithm, and complete the mesh reconstruction of the turbine blade disk.

2. The method for reconstructing the mesh of a ceramic matrix composite turbine blade disk based on CT image slices according to claim 1, characterized in that: It also includes step seven, performing a quality assessment on the tetrahedral mesh generated in step six to ensure that the geometry of the mesh meets the requirements of finite element analysis, and performing a simulation prediction of the ceramic matrix composite turbine blade based on the mesh of the turbine blade selected through the quality assessment to solve engineering problems related to the application of ceramic matrix composite turbine blades.

3. The method for reconstructing the grid of a ceramic matrix composite turbine blade disk based on CT image slices according to claim 1, characterized in that: The implementation method of step one is: Step 1.1: Preparation: Fix the CMC turbine blade disk on the sample stage of the CT scanner to ensure that it is stable during the scanning process; Step 1.2: Parameter setting: According to the material properties of the turbine blade and the required resolution, set the voltage, power, exposure time and scanning time during CT scanning; Step 1.3: Scanning process: Start the CT scanner and perform continuous tomographic scanning to obtain grayscale image slice data of the turbine blade disk layer by layer.

4. The method for reconstructing the grid of a ceramic matrix composite turbine blade disk based on CT image slices according to claim 1, characterized in that: The implementation method of step 2 is: Step 2.1: Image preprocessing: Use the median filter denoising algorithm to preprocess the original CT image to remove noise and artifacts in the image; median filter formula: I filtered (x,y)=median{I(x+i,y+j)|(i,j)∈W} Where: I filtered (x, y) is the pixel value of the filtered image, I(x+i, y+j) is the pixel value of the original image at (x+i, y+j), W is the filter window, and i and j are the offsets of the pixel positions; Step 2.2: Image enhancement: Apply histogram equalization technology to enhance the image contrast, making the boundaries between materials and gaps more obvious; histogram equalization formula: Where: I equalized (x, y) is the pixel value of the image after equalization, L is the grayscale level of the image, and MN is the total pixel value of the image number, h(i) is the number of pixels with gray value i; I(x,y) is the pixel value of the original image; Step 2.3: Threshold segmentation: Use Otsu adaptive threshold segmentation algorithm to binarize the preprocessed image and convert the image into a binary image; Where: is the inter-class variance, t is the threshold, ω1(t) and ω2(t) are the pixel probabilities of the foreground and background respectively, μ1(t) and μ2(t) are the average grayscale values ​​of the foreground and background respectively; select t that maximizes the inter-class variance as the optimal threshold.

5. The method for reconstructing the mesh of a ceramic matrix composite turbine blade disk based on CT image slices according to claim 1, characterized in that: The implementation method of step seven is: Step 7.1: Calculate the aspect ratio A: Where: E is the side length of the tetrahedron, and the calculation formula for the side length E is: Where (x a ,y a ,z a ) and (x b ,y b ,z b ) are the coordinates of the vertices at both ends of the connecting edge E, max(E) and min(E) are the longest and shortest edges respectively; the tetrahedral mesh with A<2 is a good mesh, the tetrahedral mesh with 2≤A<3 is a better mesh, and the tetrahedral mesh with A>3 is a worse mesh.

6. The method for reconstructing the mesh of a ceramic matrix composite turbine blade disk based on CT image slices according to claim 1, characterized in that: The core of the Laplacian optimization described in step 5.2 is the encoding and decoding of the local detail features of the mesh model. In the encoding stage, the Euclidean space coordinates of the mesh vertices are converted to Laplacian coordinates, which can accurately capture the local detail features of the mesh. Therefore, the Laplacian mesh deformation algorithm can retain these details when processing the mesh. In the decoding stage, the Laplacian coordinates are reversed back to the Euclidean space coordinates by solving the linear system, thereby restoring and optimizing the original mesh. w The Laplacian coordinates are defined as: P w =v w +λ·∑ j∈N(w) v j -v w Where: v w is the original position of vertex w, P w is the new position of vertex w after smoothing, N(w) represents the set of vertices adjacent to vertex w, λ is the relaxation factor that controls the smoothing strength, 0<λ<1, v j is the position of vertex j, j is the vertex adjacent to w.

7. The method for reconstructing the mesh of a ceramic matrix composite turbine blade disk based on CT image slices according to claim 1, characterized in that: The meshing software is Gmsh Python API or Hypermesh.

Images