A computer numerical simulation method, system and device for a composite material

Abstract

The present invention discloses a computer numerical simulation method, system and device for composite materials, specifically relating to the field of simulation technology, and is used to solve the problem of simulation distortion caused by ignoring the randomness of the microstructure in the prior art; it obtains the microscopic structure image data of the composite material, extracts the fiber arrangement direction and pore distribution parameters, generates the fiber direction probability density function and the pore topological relationship; constructs a non-uniform grid model with local coordinate system constraints according to the fiber direction probability density function; identifies the weak connection regions and marks the nodes based on the pore topological relationship; applies dynamic decoupling constraints to the weak connection regions, including the initial stiffness and the stiffness decay rate; applies an external load to the model and calculates the real-time stress coefficient of the weak connection regions; dynamically adjusts the constraint parameters through threshold judgment and iterative correction until the stress coefficient converges; finally outputs the corrected macroscopic mechanical property results.

Description

Technical Field

[0001] The present invention relates to the field of simulation technology, and more specifically, to a computer numerical simulation method, system and device for composite materials. Background Art

[0002] Composite materials are widely used in fields such as aerospace and automotive manufacturing due to their excellent properties. Their design optimization and performance evaluation highly rely on computer simulation technology. In the prior art, the multi-scale numerical simulation of composite materials usually constructs a model based on the finite element method, and maps the microscopic structural characteristics to the macroscopic mechanical response through the homogenization theory. However, there is natural randomness in the microscopic structure of actual composite materials (such as fiber arrangement, pore distribution, etc.). In order to simplify the calculation model, the existing methods generally use the homogenization hypothesis to describe the material properties, resulting in the difficulty of the simulation process to truly reflect the dynamic influence of microscopic non-uniformity on macroscopic behavior.

[0003] Due to the neglect of the random characteristics of the microscopic structure of composite materials by the existing simulation methods, there are systematic deviations between the predicted results of their macroscopic properties and the actual material behavior. Especially when analyzing the mechanical response of non-uniformly reinforced composite materials, the homogenization hypothesis will cover up key mechanisms such as local defect evolution and interface failure, making the simulation model unable to accurately capture the true failure mode of the material under complex loads, thus reducing the reliability of design optimization. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, the embodiments of the present invention provide a computer numerical simulation method, system and device for composite materials to solve the problems raised in the above background art.

[0005] To achieve the above object, the present invention provides the following technical solutions:

[0006] A computer numerical simulation method for composite materials, comprising the following steps:

[0007] S1. Obtain the microscopic structure image data of the composite material, extract the fiber arrangement direction and pore distribution parameters, and generate the fiber direction probability density function and the pore topological relationship;

[0008] S2. Generate a local coordinate system according to the fiber direction probability density function, and perform element division with direction constraints on the finite element mesh based on the local coordinate system to obtain a non-uniform mesh model;

[0009] S3. Based on the pore topological relationship, identify the weak connection regions between adjacent elements in the non-uniform mesh model and mark the nodes;

[0010] S4. Apply dynamic decoupling constraints to the nodes in the weak connection regions, and the dynamic decoupling constraints include the initial stiffness and the stiffness decay rate;

[0011] S5. Apply an external load to the model after applying the dynamic decoupling constraint, and calculate the real-time stress coefficient in the weak connection area;

[0012] S6. Determine whether the real-time stress coefficient exceeds a preset threshold; if it exceeds, iteratively correct the dynamic decoupling constraint based on the stiffness decay rate;

[0013] S7. After the real-time stress coefficient is lower than the preset threshold, output the corrected macroscopic mechanical property results.

[0014] In a preferred embodiment, S1 includes:

[0015] S1a. Process the microscopic structure image data through the gray gradient analysis method, identify the fiber edge contour and calculate the main axis direction angle of each fiber, and statistically generate the fiber direction probability density function by the distribution frequency of the main axis direction angle;

[0016] S1b. Perform binary segmentation on the same image data, extract the pore geometry and spatial position, and construct an adjacency relationship matrix representing the pore topological characteristics based on the spatial position relationship between the pores and the fibers;

[0017] S1c. Discretely store the fiber direction probability density function in a preset angle interval, and merge it with the adjacency relationship matrix into a statistical model representing the randomness of the microscopic structure.

[0018] In a preferred embodiment, S2 includes:

[0019] S2a. According to the angle interval distribution in the fiber direction probability density function, assign a corresponding local coordinate system to each finite element unit, and the X-axis direction of the local coordinate system is consistent with the median direction of the main axis direction angle of the fiber;

[0020] S2b. Based on the X-axis direction of the local coordinate system, adjust the stiffness matrix direction of each finite element unit through a rotation matrix, so that the main axis direction of the unit stiffness matrix is aligned with the fiber direction;

[0021] S2c. Perform mesh division on the unit after adjusting the stiffness matrix to generate a non-uniform mesh model embedding the randomness of the fiber direction.

[0022] In a preferred embodiment, S3 includes:

[0023] S3a. Based on the adjacency relationship matrix in the pore topological relationship, locate the adjacent unit pairs covering the pore area, and calculate the coincidence degree between the interface connection area of the adjacent unit pairs and the pore projection area;

[0024] S3b. If the coincidence degree exceeds a preset ratio threshold, determine that there is a weak connection area between the adjacent unit pairs, and record the numbers of the unit pairs and the coordinates of the common nodes;

[0025] S3c. Mark the common nodes in the weak connection area in the non-uniform grid model according to the recorded unit pair numbers and node coordinates, and generate a node marking list.

[0026] In a preferred embodiment, S4 includes:

[0027] S4a. Set the initial binding stiffness according to the material properties of the nodes in the weak connection area, and the initial binding stiffness is lower than the standard stiffness value of the nodes in the non-weak connection area;

[0028] S4b. Determine the stiffness decay rate of the nodes in the weak connection area based on a preset stiffness decay rate table, and the stiffness decay rate table is dynamically adjusted according to the difference between the real-time stress coefficient and the preset threshold;

[0029] S4c. Associate the initial binding stiffness with the stiffness decay rate, and apply dynamic decoupling constraints to the nodes in the weak connection area. The dynamic decoupling constraints are realized by iteratively updating the node stiffness values.

[0030] In a preferred embodiment, S5 includes:

[0031] S5a. Define the external load application area and load direction in the model after the dynamic decoupling constraints, and the load direction is consistent with the actual working conditions of the composite material;

[0032] S5b. Calculate the real-time stress distribution of the nodes in the weak connection area based on a finite element solver, and extract the maximum principal stress or Von Mises stress as the real-time stress coefficient;

[0033] S5c. Calculate the stress concentration coefficient evaluation index according to the stress gradient distribution of the nodes in the weak connection area, and the stress concentration coefficient evaluation index is the ratio of the node stress to the average stress of the adjacent non-weak connection area;

[0034] S5d. Screen the nodes in the weak connection area whose stress concentration coefficient evaluation index exceeds the preset threshold, and record their real-time stress coefficients and corresponding coordinate information.

[0035] In a preferred embodiment, S6 includes:

[0036] S6a. Judge whether the real-time stress coefficient exceeds the preset threshold. If it exceeds, calculate the difference interval between the real-time stress coefficient and the preset threshold;

[0037] S6b. Query the stiffness decay rate table according to the difference interval, and match the corresponding stiffness decay rate correction step;

[0038] S6c. Update the stiffness decay rate of the dynamic decoupling constraints based on the correction step, and recalculate the real-time stress coefficients of the nodes in the weak connection area.

[0039] In a preferred embodiment, S7 includes:

[0040] S7a. Verify that the real-time stress coefficient is lower than a preset threshold and the change rate of the real-time stress coefficient between two adjacent iterations is less than a preset convergence accuracy;

[0041] S7b. Extract the final stiffness values and stress distribution data of the nodes in the weak connection area from the modified non-uniform grid model as key parameters;

[0042] S7c. Associate the key parameters with the homogenized constitutive equation of the macroscopic model to generate a macroscopic mechanical property report file containing the modified material properties;

[0043] S7d. Update the composite material model data in the finite element simulation database according to the report file and output the updated macroscopic mechanical property results.

[0044] On the other hand, the present invention provides a computer numerical simulation system for composite materials, including:

[0045] Microstructure acquisition module: Obtain the microscopic structure image data of the composite material, extract the fiber arrangement direction and pore distribution parameters, and generate a fiber direction probability density function and a pore topological relationship;

[0046] Grid orientation module: Generate a local coordinate system according to the fiber direction probability density function, and perform element division with direction constraints on the finite element grid based on the local coordinate system to obtain a non-uniform grid model;

[0047] Weak area identification module: Based on the pore topological relationship, identify the weak connection area between adjacent elements in the non-uniform grid model and mark the nodes;

[0048] Dynamic constraint module: Apply dynamic decoupling constraints to the nodes in the weak connection area, and the dynamic decoupling constraints include an initial stiffness and a stiffness decay rate;

[0049] Load application module: Apply an external load to the model after applying the dynamic decoupling constraints and calculate the real-time stress coefficient of the weak connection area;

[0050] Iterative correction module: Determine whether the real-time stress coefficient exceeds a preset threshold; if it exceeds, iteratively correct the dynamic decoupling constraints based on the stiffness decay rate;

[0051] Result output module: When the real-time stress coefficient is lower than the preset threshold, output the corrected macroscopic mechanical property results.

[0052] On the other hand, the present invention provides a computer numerical simulation device for composite materials, which includes: a processor, a memory, and a program or instruction stored on the memory and executable on the processor. When the program or instruction is executed by the processor, a computer numerical simulation method for composite materials is implemented.

[0053] Compared with the prior art, the present invention has the following beneficial effects:

[0054] 1. By extracting microstructure image data and constructing probability density functions, the randomness of fiber direction and the characteristics of pore space distribution are embedded into the finite element mesh modeling, enabling the simulation model to accurately represent the spatial correlation of micro-defects, thereby effectively capturing the dynamic impact of non-uniformity on the overall performance of the material. Compared with the traditional homogenization method, this method not only solves the problem of mapping distortion between micro-random characteristics and macro-mechanical behaviors, but also provides a modeling basis closer to the actual working conditions for multi-scale simulation of composite materials.

[0055] 2. Through the node stiffness iterative correction and threshold judgment logic in the weak connection area, a dynamic decoupling path based on the real-time stress state is constructed, enabling the simulation process to adaptively adjust local constraint conditions and truly reflect the progressive characteristics of interface failure. The dynamic coupling simulation strategy significantly improves the prediction accuracy of composite material damage evolution. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 is a flowchart of a computer numerical simulation method for composite materials according to the present invention;

[0057] Figure 2 is a schematic structural diagram of a computer numerical simulation system for composite materials according to the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0058] The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0059] Embodiment 1: Figure 1 A computer numerical simulation method for composite materials according to the present invention is given, which includes the following steps:

[0060] S1. Obtain the microstructure image data of the composite material, extract the fiber arrangement direction and pore distribution parameters, and generate the fiber direction probability density function and pore topological relationship;

[0061] S2. Generate a local coordinate system according to the fiber direction probability density function, and perform element division with direction constraints on the finite element mesh based on the local coordinate system to obtain a non-uniform mesh model;

[0062] S3. Based on the pore topological relationship, identify the weak connection regions between adjacent elements in the non-uniform mesh model and mark the nodes;

[0063] S4. Apply dynamic decoupling constraints to the nodes in the weak connection regions, and the dynamic decoupling constraints include the initial stiffness and the stiffness decay rate;

[0064] S5. Apply an external load to the model after applying the dynamic decoupling constraints, and calculate the real-time stress coefficient of the weak connection regions;

[0065] S6. Determine whether the real-time stress coefficient exceeds a preset threshold; if it exceeds, iteratively correct the dynamic decoupling constraints based on the stiffness decay rate;

[0066] S7. After the real-time stress coefficient is lower than the preset threshold, output the corrected macroscopic mechanical property results.

[0067] S1. Obtain the microscopic structure image data of the composite material, extract the fiber arrangement direction and pore distribution parameters, and generate the fiber direction probability density function and the pore topological relationship, including:

[0068] Process the microscopic structure image data of the composite material through the gray gradient analysis method, identify the fiber edge contours and calculate the principal axis direction angles of each fiber, and statistically generate the fiber direction probability density function of the distribution frequency of the fiber principal axis direction angles; the gray gradient analysis method first performs graying processing on the microscopic structure image to convert the color image into a gray image, and the graying processing converts the pixel values of the red, green, and blue channels into a single gray value through the weighted average method.

[0069] Determine the fiber edge contours by calculating the gradient values of each pixel point in the gray image, and the gradient values are calculated by the Sobel operator. The Sobel operator uses a 3×3 convolution kernel to calculate the gradient components of the pixel point in the horizontal and vertical directions respectively, and the gradient amplitude is the square root of the sum of the squares of the horizontal and vertical components; mark the pixel points with gradient amplitudes greater than the preset threshold as fiber edge points, and the preset threshold is dynamically adjusted according to the image contrast, for example, automatically determined by the Otsu algorithm; according to the spatial distribution of the fiber edge points, use the least squares method to fit the principal axis direction angles of each fiber. The least squares method determines the optimal straight line by minimizing the sum of the squares of the perpendicular distances from the fiber edge points to the fitted straight line, and the included angle between the fitted straight line and the horizontal axis of the image is the fiber principal axis direction angle.

[0070] Statistically analyze the distribution frequency of the principal axis direction angles of all fibers to generate a fiber direction probability density function. The statistical method is to divide the angle values into preset intervals (such as 5 degrees or 10 degrees), count the number of occurrences of the principal axis direction angles of the fibers in each interval, and convert the number of occurrences into probability density values through normalization. The normalization process is to divide the number of occurrences in each interval by the total number of fibers.

[0071] Perform binary segmentation on the same microscopic structure image data, extract the geometric shape and spatial position of the pores, and construct an adjacency relation matrix that characterizes the topological features of the pores based on the spatial position relationship between the pores and the fibers.

[0072] Binary segmentation divides the image into a pore region and a non-pore region by setting a gray threshold. The pixel values in the pore region are set to 1, and the pixel values in the non-pore region are set to 0. The gray threshold is automatically calculated by the Otsu algorithm or manually set according to experience; noise points are removed and the pore edges are smoothed through morphological operations (such as erosion and dilation). The erosion operation traverses the image with a 3×3 rectangular structuring element, and the dilation operation restores the pore shape with the same structuring element.

[0073] The geometric shape parameters of the pores extracted include area, perimeter, and aspect ratio. The area is calculated by counting the number of pixels in the pore region. The perimeter is calculated by the chain code method or the edge tracking algorithm to count the number of pixel points on the pore boundary. The aspect ratio is determined by the ratio of the long side to the short side of the minimum bounding rectangle that fits the pore.

[0074] The spatial position is represented by the centroid coordinates of the pore region, and the centroid coordinates are obtained by calculating the average of the horizontal and vertical coordinates of all pixel points in the pore region; an adjacency relation matrix is constructed based on the spatial position relationship between the pores and the fibers. The rows and columns of the adjacency relation matrix represent the numbers of the pores and the fibers respectively, and the matrix elements are boolean values or weight values. The boolean value indicates whether the pore and the fiber are adjacent, and the definition of adjacency is that the minimum distance from the pore edge to the fiber edge is less than a preset threshold (such as 5 pixels). The weight value is calculated according to the overlapping area or the edge distance between the pore and the fiber, such as the proportion of the overlapping area to the total pore area or the normalized Euclidean distance.

[0075] The storage form of the adjacency relation matrix is a two-dimensional array or a sparse matrix. The sparse matrix stores non-zero elements through triples (row number, column number, element value) to save memory space.

[0076] Discretely store the fiber direction probability density function in a preset angle interval and merge it with the adjacency relation matrix into a statistical model that characterizes the randomness of the microscopic structure.

[0077] The fiber direction probability density function is discretized into a sequence of probability values according to a preset angular interval (e.g., 5 degrees or 10 degrees). Each interval corresponds to a probability value. The data structure for storing the discretized data is an array or a list, and the array index corresponds to the starting angle of the angular interval.

[0078] The adjacency relationship matrix is stored in the form of a two-dimensional array or a sparse matrix. The row numbers and column numbers of the matrix correspond to the unique identifiers of the pores and fibers. When combining the fiber direction probability density function and the adjacency relationship matrix into a statistical model, the data structure of the statistical model is a dictionary or a JSON file. The keys of the dictionary are "fiber direction probability density" and "pore adjacency relationship", and the corresponding values are the discretized probability value sequence and the adjacency relationship matrix.

[0079] During the combination process, the data of the fiber direction probability density function and the adjacency relationship matrix are associated through unique identifiers. The unique identifiers are the numbers of the fibers and pores, and the numbers are generated by an image coordinate hashing algorithm. For example, the centroid coordinates are converted into a string and the first 8 bits of the MD5 hash value are taken.

[0080] The statistical model is used to describe the randomness characteristics of the microstructure of the composite material, such as the distribution law of the fiber direction and the spatial adjacency relationship between the pores and the fibers. In the subsequent steps, the finite element mesh generation and mechanical simulation are carried out by reading the data in the statistical model.

[0081] During the recognition process of the fiber edge contour, the Sobel operator is used to calculate the gradient magnitude. The horizontal convolution kernel of the Sobel operator is [-1,0,1;-2,0,2;-1,0,1], and the vertical convolution kernel is [-1,-2,-1;0,0,0;1,2,1]. The horizontal and vertical gradient components are obtained by performing convolution operations on the grayscale image respectively. During the fitting process of the fiber major axis direction angle, the least squares method calculates the slope and intercept of the fitting line by solving a system of linear equations. The slope is the tangent value of the fiber major axis direction angle. The normalization process of the fiber direction probability density function is to divide the frequency values of each angular interval by the total number of fibers. The total number of fibers is obtained by analyzing the number of independent fiber regions in the image through connected component analysis.

[0082] During the extraction process of the pore geometry, morphological operations are used to remove noise points. The erosion operation is achieved by taking the minimum value of the pixel values within the area covered by the structuring element, and the dilation operation is achieved by taking the maximum value.

[0083] The calculation formula for the centroid coordinates is that the abscissa of the centroid is equal to the sum of the abscissas of all pixel points within the pore region divided by the number of pixels, and the ordinate of the centroid is calculated in the same way.

[0084] During the construction of the adjacency relationship matrix, the minimum distance between pores and fibers is calculated by traversing the coordinates of the edge pixels of pores and the edge pixels of fibers to calculate the Euclidean distance, and the minimum value among all distances is taken as the distance between the two; the calculation of the weight value is, for example, the normalized overlapping area, that is, the overlapping area of the pore and the fiber divided by the total area of the pore, and the overlapping area is counted by the pixel-level logical AND operation for the number of pixels where both are 1.

[0085] In the discretely stored fiber direction probability density function, the division of the angular interval needs to cover the range from 0 degrees to 180 degrees to avoid the periodic repetition of the direction angle; during the merging process of the adjacency relationship matrix, the generation of unique identifiers needs to ensure that the numbers of fibers and pores are unique and can be regenerated repeatedly. The hash algorithm is implemented by converting the floating-point centroid coordinates into a string and then taking the hash value; an example of the data structure of the statistical model is in JSON format, including an array of discretized probability values corresponding to the key "fiber_orientation", and a sparse storage format of the adjacency relationship matrix corresponding to the key "pore_adjacency"; in the subsequent steps, by parsing the data structure in the statistical model, the fiber direction probability density function is used for the material property assignment of finite element cells, and the adjacency relationship matrix is used to identify weak connection regions and apply dynamic decoupling constraints.

[0086] In the analysis of the spatial position relationship between pores and fibers, if the minimum distance between pores and fibers is less than the preset threshold, it is determined that the two are adjacent. The preset threshold is set according to the image resolution. For example, when the resolution is 1 micron / pixel, the threshold is set to 5 pixels, which means 5 microns; the weight value of the adjacency relationship matrix can be selected as a boolean value or a continuous value according to actual needs. The boolean value is suitable for qualitative analysis, and the continuous value is suitable for quantitative simulation; the sequence of discretized probability values in the statistical model can be converted into a continuous probability density function by interpolation methods, such as linear interpolation or cubic spline interpolation, and the choice of interpolation method is determined according to the simulation accuracy requirements.

[0087] An example of the discretized storage of the fiber direction probability density function is that the probability value corresponding to the angular interval 0 - 5 degrees is 0.02, and the probability value corresponding to 5 - 10 degrees is 0.03, and so on. The sum of the probability values of all intervals is 1.

[0088] An example of the sparse matrix storage of the adjacency relationship matrix is to only record the triple (row number, column number, element value) of non-zero elements. For example, (1,5,0.8) means that the overlapping area of pore No. 1 and fiber No. 5 is 0.8.

[0089] An example of the hash algorithm for unique identifiers is to convert the centroid coordinates (123.45,67.89) into the string "123.45_67.89", and then take the first 8 characters "a1b2c3d4" of the MD5 hash value as the number.

[0090] S2. Generate a local coordinate system according to the fiber direction probability density function, and perform element division with direction constraints on the finite element mesh based on the local coordinate system to obtain a non-uniform mesh model, including:

[0091] According to the angular interval distribution in the fiber direction probability density function, assign a corresponding local coordinate system to each finite element cell; the angular interval distribution of the fiber direction probability density function is a preset discretized angular interval, such as 0 - 5 degrees, 5 - 10 degrees, etc., and each angular interval corresponds to a probability value; the X-axis direction of the local coordinate system is consistent with the median direction of the fiber principal axis direction angle within the angular interval, and the median direction is the middle angular value of the angular interval. For example, the median direction of the 0 - 5 degree interval is 2.5 degrees.

[0092] The local coordinate system of each finite element cell is determined according to the median direction of the angular interval it belongs to. For example, if the angular interval a certain cell belongs to is 10 - 15 degrees, then the X-axis direction of its local coordinate system is 12.5 degrees; the Y-axis direction of the local coordinate system is perpendicular to the X-axis and follows the right-hand rule.

[0093] Based on the X-axis direction of the local coordinate system, adjust the stiffness matrix direction of each finite element cell through a rotation matrix; the rotation matrix is used to transform the element stiffness matrix in the global coordinate system to the local coordinate system, and the construction of the rotation matrix is based on the angle between the local coordinate system and the global coordinate system.

[0094] For example, if the angle between the X-axis direction of the local coordinate system and the X-axis of the global coordinate system is θ, then the rotation matrix is a two-dimensional rotation matrix containing cosθ and sinθ elements.

[0095] The principal axis direction of the adjusted element stiffness matrix is aligned with the fiber direction. For example, the principal direction of the element stiffness matrix in the local coordinate system is along the X-axis direction, which is consistent with the fiber principal axis direction.

[0096] The adjustment of the element stiffness matrix is achieved through matrix multiplication, that is, multiplying the rotation matrix by the original stiffness matrix to obtain the stiffness matrix in the local coordinate system.

[0097] Perform mesh division on the elements after adjusting the stiffness matrix to generate a non-uniform mesh model that embeds the randomness of the fiber direction; the mesh division is based on the mesh generation module of the finite element software. For example, use hexahedral or tetrahedral elements to discretize the composite material geometric model.

[0098] During the mesh division process, the stiffness matrix direction of each element is set according to the X-axis direction of its local coordinate system to ensure that the material property direction of the element is consistent with the fiber direction.

[0099] The cell orientation distribution in the non-uniform grid model matches the statistical law of the fiber orientation probability density function. For example, the higher the probability value corresponding to a certain angular interval, the larger the proportion of the number of cells in that direction. After the grid division is completed, each cell in the non-uniform grid model contains the adjusted stiffness matrix direction information for subsequent multi-scale mechanical simulations.

[0100] During the allocation process of the local coordinate system, if the discretization angular interval of the fiber orientation probability density function is 5-degree intervals, the median direction of each interval is the starting angle of the interval plus 2.5 degrees.

[0101] For example, the median direction of the angular interval 15 - 20 degrees is 17.5 degrees, corresponding to the X-axis direction of the local coordinate system. An example of constructing the rotation matrix is a two-dimensional rotation matrix, whose elements consist of cosθ and -sinθ in the first row, and sinθ and cosθ in the second row, where θ is the radian value corresponding to the median direction.

[0102] An example of adjusting the cell stiffness matrix is to convert the isotropic stiffness matrix in the global coordinate system into an anisotropic stiffness matrix in the local coordinate system, with the main direction along the X-axis. During the grid division process, the cell size of the non-uniform grid model is set according to the simulation accuracy requirements. For example, a finer cell division is used in the high stress gradient region.

[0103] The data structure of the non-uniform grid model includes node coordinates, element connection relationships, and the stiffness matrix direction information of each cell. The data structure format is the format supported by finite element software, such as INP files or CDB files.

[0104] The adjusted stiffness matrix direction is aligned with the fiber direction. For example, if the X-axis direction of the local coordinate system of a certain cell is 30 degrees, the main direction of its stiffness matrix forms an angle of 30 degrees with the X-axis direction of the global coordinate system.

[0105] During the generation process of the non-uniform grid model, if the fiber orientation probability density function shows that the probability value of a certain angular interval is 0, no cells are allocated in that direction. For example, when the probability value of the angular interval 170 - 175 degrees is 0, there are no corresponding cells in that direction. The verification of the grid division is achieved by checking the consistency between the cell stiffness matrix direction and the local coordinate system. For example, several cells are randomly selected to verify whether the deviation of the main direction of their stiffness matrices from the X-axis direction of the local coordinate system is less than the preset tolerance.

[0106] When the non-uniform grid model is used for the mechanical simulation in the subsequent steps, the external load boundary conditions are set according to the actual working conditions. For example, a tensile load is applied along the X-axis direction of the global coordinate system, and the simulation results include stress, strain, and damage distribution data.

[0107] S3. Based on the pore topological relationship, identify the weak connection regions between adjacent cells in the non-uniform grid model and mark the nodes, including:

[0108] Based on the adjacency relation matrix in the pore topological relationship, locate the adjacent unit pairs covering the pore region, and calculate the coincidence degree between the interface connection area of the adjacent unit pairs and the pore projection area.

[0109] The rows and columns of the adjacency relation matrix represent the numbers of pores and fibers respectively, and the matrix elements are boolean values or weight values. When the boolean value is 1, it means the pore is adjacent to the fiber; the adjacent unit pairs covering the pore region refer to two units sharing a boundary in the finite element mesh, and one of the units belongs to the pore region and the other belongs to the fiber region.

[0110] The interface connection area is the geometric area of the shared boundary of the adjacent unit pairs, which is obtained by calculating the length of the shared boundary multiplied by the unit thickness; the pore projection area is the projection area of the pore region on the plane of the shared boundary, which is calculated by counting the number of pore pixels in the projection plane multiplied by the actual size of the pixel; the coincidence degree is the ratio of the interface connection area to the pore projection area. For example, if the interface connection area is 10 square micrometers and the pore projection area is 15 square micrometers, the coincidence degree is 66.7%.

[0111] If the coincidence degree exceeds the preset ratio threshold, it is determined that there is a weak connection region between the adjacent unit pairs, and the numbers of the unit pairs and the coordinates of the common nodes are recorded; the preset ratio threshold is set according to the experimental data of the material interface strength. For example, when the coincidence degree exceeds 50%, it is determined as a weak connection region.

[0112] The number of the unit pair is the globally unique identifier of the unit in the finite element mesh. For example, the sequential number of the unit in the mesh file is adopted; the coordinates of the common nodes are the coordinates of the nodes shared by the adjacent unit pairs, and the node coordinates are obtained from the node coordinate list in the finite element mesh model. The node coordinate list stores the X, Y, and Z coordinate values of each node; the recording method is to store the number of the unit pair and its common node coordinates as a list, and the list structure is [unit pair number 1, node coordinates 1; unit pair number 2, node coordinates 2;...].

[0113] According to the recorded numbers of the unit pairs and the node coordinates, mark the common nodes of the weak connection region in the non-uniform mesh model to generate a node marking list; the marking method is to traverse the recorded numbers of the unit pairs, locate the corresponding units and their common nodes, and add the node identifiers to the node marking list.

[0114] The data structure of the node marking list is an array or a table, which contains the node identifiers and their coordinate values, and the node identifier is the unique number of the node in the finite element mesh model.

[0115] After generating the node label list, the nodes in the weak connection area of the non-uniform grid model are labeled with special attributes, such as adding a "weak connection" label in the node attribute field; the node label list is used for the imposition of dynamic decoupling constraints in subsequent steps, such as adjusting the stiffness parameters only for the labeled nodes.

[0116] In the example of calculating the coincidence degree, if the interface connection area is 8 square micrometers and the pore projection area is 12 square micrometers, the coincidence degree is 66.7%; the preset ratio threshold is set according to the material type, for example, 40% for brittle materials and 60% for ductile materials; an example of recording the element pair number is that element number 5 and element number 12 form an adjacent element pair, and their common nodes are node 7, node 8, and node 9; an example of obtaining the node coordinates is the coordinates of node 7 (10.2, 5.3, 0.0) and the coordinates of node 8 (10.5, 5.3, 0.0); an example of generating the node label list is adding node 7, node 8, and node 9 to the list and marking their coordinates and the element pair number they belong to.

[0117] During the determination process of the weak connection area, if the coincidence degree is lower than the preset threshold, the adjacent element pair is ignored, for example, when the coincidence degree is 30%, it is determined as a normal connection area; the storage format of the node label list is compatible with finite element software, such as the node set (Node Set) of Abaqus or the named selection set (Named Selection) of ANSYS; an example of the subsequent processing of the labeled nodes is to reduce their bonding stiffness or introduce a contact failure criterion in finite element analysis; the node label information in the non-uniform grid model is implemented by modifying the grid file, such as adding an attribute field in the node definition part of the INP file.

[0118] The calculation of the coincidence degree is based on the actual geometric projection, for example, projecting the three-dimensional shape of the pore area onto the shared boundary plane of the adjacent element pair, and the projection method is orthogonal projection or inclined projection along the fiber direction; the calculation of the interface connection area takes into account the actual thickness of the element, for example, the single-layer thickness of the composite laminate is 0.1 millimeters; an example of the statistics of the pore projection area is using image processing software (such as ImageJ) to measure the number of pixels in the projection area and convert it to the actual size; an example of recording the element pair number is storing the numbers of the adjacent element pairs as a two-dimensional array, with the first column being the pore element number and the second column being the fiber element number.

[0119] An application example of the node marking list is to read list data in a finite element solver and apply dynamic decoupling constraints to the marked nodes; an example of parameter adjustment for dynamic decoupling constraints is to modify the contact properties of the corresponding nodes in the solver input file according to the node numbers in the node marking list; the node marking information in the non-uniform grid model is retained until the simulation is completed for subsequent result analysis, such as extracting the stress concentration coefficient in the weak connection region; the verification method of the node marking list is visual inspection, such as highlighting the marked nodes in a finite element post-processing software to confirm that they are located in the interface region between pores and fibers.

[0120] The setting basis of the coincidence degree threshold includes material interface peel strength test data, such as measuring the critical peel force of the interface through a double-cantilever beam test and converting it into an area ratio threshold; the global uniqueness of the element pair numbers is guaranteed by the finite element mesh generation algorithm, such as assigning an increasing integer number to each element during mesh division; the accuracy of the node coordinates is consistent with the mesh division resolution, for example, when the mesh size is 1 micron, the node coordinates are retained to three decimal places; the generation efficiency of the node marking list is achieved through batch processing, such as traversing all recorded element pair numbers, extracting all common nodes at once and removing duplicates.

[0121] An example of marking the weak connection region is to identify the high coincidence degree region at the pore-fiber interface of a composite laminate, and the marked nodes are used to simulate interface delamination failure; an example of data storage of the node marking list is a CSV file, which contains node numbers, X coordinates, Y coordinates, Z coordinates and the corresponding element pair numbers; an example of applying dynamic decoupling constraints in subsequent steps is to set the initial stiffness of the marked nodes to 50% of the standard value and gradually reduce it according to the stiffness decay rate; the simulation results of the non-uniform grid model show that the stress concentration coefficient in the weak connection region is significantly higher than the predicted value of the homogenized model, verifying the effectiveness of the marking logic.

[0122] The selection of the projection plane in coincidence degree calculation is determined according to the fiber direction, such as projecting along the main axis direction of the fiber to reflect the actual interface contact situation; an example of calculating the interface connection area is the shared boundary length of rectangular elements multiplied by the element thickness, and the shared boundary arc length of circular elements multiplied by the thickness; an example of measuring the projected area of pores is to use metallographic microscope images and calculate the proportion of pores in the projected area through image analysis software; an example of the recording format of element pair numbers is a text file, with each line recording a pair of element numbers and their coincidence degree values; an example of the generation tool for the node marking list is a Python script, which parses the element pair records, extracts the node coordinates, and outputs them as a node set file.

[0123] After the nodes in the weak connection area are marked, the marked nodes in the non-uniform grid model participate in the dynamic decoupling constraint calculation in subsequent simulations; an example of the iterative correction of the dynamic decoupling constraint is to adjust the stiffness decay rate of the marked nodes according to the real-time stress coefficient until the stress is lower than the threshold; an example of the maintenance of the node marking list is to synchronously update the node numbers and coordinates during grid refinement or coarsening to ensure the consistency of the marking information; an example of the simulation output of the non-uniform grid model includes the stress history data of the marked nodes, which is used to analyze the damage evolution law of the weak connection area.

[0124] S4. Apply dynamic decoupling constraints to the nodes in the weak connection area. The dynamic decoupling constraints include the initial stiffness and the stiffness decay rate, including:

[0125] Set the initial binding stiffness according to the material properties of the nodes in the weak connection area. The material properties include the interface strength parameters of the fiber and the matrix, and the interface strength parameters are obtained through single fiber pull-out tests or micro-droplet debonding tests. The initial binding stiffness is set proportionally according to the interface strength parameters. For example, when the interface strength is 50 MPa, the initial binding stiffness is set to 50% of the standard stiffness value. The standard stiffness value is the stiffness value of the nodes in the non-weak connection area, such as the stiffness value corresponding to an elastic modulus of 100 GPa. The setting method of the initial binding stiffness is to modify the node stiffness parameters in the finite element software, such as reducing the elastic modulus field value of the weak connection node to a preset value.

[0126] Determine the stiffness decay rate of the nodes in the weak connection area based on a preset stiffness decay rate table. The stiffness decay rate table is a two-dimensional table, where the rows represent the difference intervals between the real-time stress coefficient and the preset threshold, and the columns represent the corresponding stiffness decay rates. The difference intervals are calibrated through material fatigue test data. For example, a difference of 5 MPa corresponds to an interface damage accumulation rate of 0.1. The dynamic adjustment method of the stiffness decay rate is the look-up table method, and the corresponding rate value is matched according to the stress difference in the current iteration step. For example, when the real-time stress coefficient exceeds the threshold by 5 MPa, the stiffness decay rate is 0.1 / iteration step.

[0127] Associate the initial binding stiffness with the stiffness decay rate and apply dynamic decoupling constraints to the nodes in the weak connection area. The initial binding stiffness is used as the starting point of the iteration. For example, the initial stiffness value is 50 MPa and the stiffness decay rate is 0.1 / iteration step. The dynamic decoupling constraint is realized by iteratively updating the node stiffness value, and in each iteration step, the node stiffness value is multiplied by (1 - stiffness decay rate); for example, after the first iteration, the stiffness value drops to 45 MPa, and after the second iteration, it drops to 40.5 MPa. The iteration process continues until the real-time stress coefficient is lower than the preset threshold. For example, when the threshold is set to 30 MPa, the iteration termination condition is that the node stiffness value drops below 30 MPa.

[0128] The construction of the stiffness decay rate table is based on experimental data of material interface damage accumulation. For example, the stiffness decay rate under different stress differences is measured through cyclic load tests. An example of the implementation of the look-up table method is to write a script program to read the real-time stress difference and match the rate value in the corresponding interval of the rate table. An example of the termination condition for iterative update is that the node stiffness value drops below the threshold or the number of iterative steps exceeds the maximum allowed number of steps, such as 100 steps. The update of the node stiffness value is achieved through the field variable function of the finite element software. For example, the stiffness value is used as a field variable and associated with the number of iterative steps.

[0129] An example of the application logic of the dynamic decoupling constraint is to read the real-time stress coefficient at each iterative step in the finite element analysis, calculate the difference from the preset threshold, look up the stiffness decay rate in the table, and update the node stiffness value. The update of the node stiffness value triggers the reassembly of the stiffness matrix of the finite element model and enters the calculation of the next iterative step. After the iterative termination condition is met, the final node stiffness value and the corresponding stress distribution results are output. An example of modifying the node stiffness value is to use the USDFLD subroutine in Abaqus to define the field variable and update the material parameters according to the number of iterative steps.

[0130] During the stiffness update process of the nodes in the weak connection area, if the real-time stress difference decreases to a lower interval, the rate is adjusted. For example, when the initial difference is 12 MPa, the rate is 0.2 / step, and when the difference drops to 6 MPa after iteration, it switches to 0.1 / step. The lower limit of the stiffness value is set to the elastic modulus of the material matrix. For example, when the matrix elastic modulus is 30 GPa, the node stiffness value shall not be lower than 30 GPa. An example of controlling the iterative step size is a fixed step size, that is, the step size is the same for each iteration, or an adaptive step size, that is, adjusted according to the convergence speed. An example of verifying the dynamic decoupling constraint is to compare the simulation results with the experimental data to confirm the consistency of damage evolution.

[0131] A test case of the dynamic decoupling constraint is to simulate the interface delamination of a composite laminate and verify the corresponding relationship between the stiffness decay of the marked nodes and the delamination propagation. An example of storing the node stiffness value is to record the change of the stiffness value at each iterative step in the ODB result file. An example of the industrial application of the dynamic decoupling constraint is the fatigue life prediction of aerospace composite structures, and the interface performance degradation under cyclic loads is simulated through node stiffness iteration. An example of visualizing the node stiffness value is to render the stiffness distribution cloud map of different iterative steps in the finite element post-processing software.

[0132] An example of the dynamic adjustment of the stiffness decay rate table is that for the uncovered interval, the linear interpolation method is used to estimate the rate value. For example, when the difference is 7.5 MPa, the rate is linearly interpolated to 0.15 / step based on the rates of 5 MPa (0.1 / step) and 10 MPa (0.2 / step). An example of the iterative update of the nodal stiffness value is to write an iterative loop using DMAP language in Nastran and update the stiffness parameters and resubmit the solution in each loop. An example of the parallel calculation of the dynamic decoupling constraint is to assign the stiffness update tasks of different nodes to a multi-core processor.

[0133] After the dynamic decoupling constraint is applied to the nodes in the weak connection region, the model has time-varying stiffness characteristics and can simulate the progressive failure process of interface damage. An example of the storage and call of the nodal stiffness value is to store the stiffness value as a time history variable for the post-processing module to generate the stiffness decay curve. An example of the convergence judgment of the dynamic decoupling constraint is to monitor the change rate of the stiffness value, and if the change rate is less than 1% for three consecutive iterations, it is determined to converge. An example of the maintenance of the node mark list is to synchronously update the node numbers and coordinates when the mesh is refined or coarsened.

[0134] An example of the stiffness decay rate table is that the rate corresponding to the difference of 0 - 5 MPa is 0.05 / step and the rate corresponding to 5 - 10 MPa is 0.1 / step. An example of the dynamic adjustment is that when the real-time stress difference of a certain node is 8 MPa, the rate of 0.1 / step is selected by looking up the table. An example of the iterative update is that the initial stiffness of 50 MPa drops to 36.45 MPa after three steps of iteration. An example of the modification of the nodal stiffness value is to use the APDL command loop in ANSYS to modify the nodal real constants. An example of the application of the dynamic decoupling constraint is to call the user-defined material model in the solver and dynamically adjust the stiffness parameters according to the number of iterative steps.

[0135] S5. Apply external loads to the model after applying the dynamic decoupling constraint, and calculate the real-time stress coefficient of the weak connection region, including:

[0136] Define the external load application area and load direction in the model after the dynamic decoupling constraint, and the load direction is consistent with the actual working conditions of the composite material; the external load application area is determined according to the stress position of the actual engineering structure, such as the edge or the area around the hole of the composite laminate; the load direction is set according to the type of working conditions. For example, the tensile load is applied along the fiber direction, and the shear load is applied along the interlayer direction; the load application method is to select the target nodes or elements in the finite element software and apply force or displacement boundary conditions. For example, in Abaqus, the concentrated force or pressure distribution is defined through the Load module.

[0137] Calculate the real-time stress distribution of the nodes in the weak connection area based on a finite element solver, and extract the maximum principal stress or Von Mises stress as the real-time stress coefficient; the finite element solver performs static or dynamic analysis, solves the equilibrium equation according to the current stiffness matrix and load conditions, and obtains the nodal stress tensor; the maximum principal stress is used to evaluate the fracture risk of brittle materials, and the Von Mises stress is used to evaluate the yield risk of ductile materials; the real-time stress coefficient is the nodal stress value, for example, the maximum principal stress of a certain node is 200 MPa or the Von Mises stress is 150 MPa.

[0138] Calculate the stress concentration coefficient evaluation index according to the stress gradient distribution of the nodes in the weak connection area; the stress gradient distribution is calculated by the stress difference between adjacent nodes. For example, the stress of a certain node is 200 MPa, and the average stress of its adjacent non-weak connection area nodes is 120 MPa; the stress concentration coefficient evaluation index is the ratio of the nodal stress to the average stress of the adjacent non-weak connection area, for example, 200 / 120≈1.67; the adjacent non-weak connection area is defined as the set of nodes that are directly connected to the weak connection nodes and are not marked as weak connections.

[0139] Screen the nodes in the weak connection area whose stress concentration coefficient evaluation index exceeds the preset threshold, and record their real-time stress coefficients and corresponding coordinate information; the preset threshold is set according to the allowable stress concentration coefficient of the material. For example, the maximum allowable stress concentration coefficient of the composite material interface is 1.5, and nodes exceeding this value are determined as high-risk nodes; the recording method is to store the node number, real-time stress coefficient and coordinate information as a table or list. For example, the node number N1001 corresponds to the stress coefficient 1.67 and the coordinates (10.2, 5.3, 0.0); the recorded node information is used for iterative correction or result analysis in subsequent steps.

[0140] Examples of the external load application area include the uniaxial tension condition of a composite laminate. The load application area is one end of the plate, and the load direction is along the length direction of the plate; an example of the calculation of the finite element solver is to perform a static analysis using the Abaqus / Standard solver and output the nodal stress nephogram; an example of the extraction of the maximum principal stress is to read the stress components from the ODB result file of Abaqus and calculate the maximum principal stress value; an example of the extraction of the Von Mises stress is to directly read the Mises stress field data.

[0141] An example of calculating the stress gradient distribution is to traverse the nodes in the weak connection area and calculate the stress difference between them and the adjacent non-weak connection nodes. The method for determining adjacent nodes is based on the element connection relationship. For example, if node A and node B belong to the same element and node B is not marked as a weak connection, then node B is an adjacent non-weak connection node. An example of the statistical example of the stress concentration factor evaluation index is to calculate the ratio of all weak connection nodes and sort them, and select the top 10% of the high-risk nodes. An example of the recorded node information is a CSV file, which contains fields such as "node number, stress coefficient, X coordinate, Y coordinate, Z coordinate".

[0142] The setting basis of the preset threshold is the material fatigue test data. For example, the critical stress concentration factor of interface failure is measured through cyclic load tests. An example of threshold adjustment is that if the simulation result deviates greatly from the experiment, the threshold is corrected from 1.5 to 1.4. The recorded coordinate information is used for visualizing and locating high-risk nodes. For example, in the finite element post-processing software, the nodes near the coordinates (10.2, 5.3, 0.0) are highlighted. An example of the application of node information is to give priority to processing high-risk nodes in subsequent iterative corrections and targetedly adjust their dynamic decoupling constraint parameters.

[0143] An example of loading the model after dynamic decoupling constraint is to apply a force load through the Solution module in ANSYS, and the direction is consistent with the fiber arrangement direction. An example of obtaining the real-time stress distribution is to print the node stress results through the PRNSOL command in Mechanical APDL. An example of the calculation tool for the stress concentration factor evaluation index is a Python script, which automatically traverses the node data and outputs a list of ratios. An example of threshold screening is to use the Pandas library to filter the records in the CSV file where the stress concentration factor is greater than 1.5.

[0144] An example of the stress gradient analysis of the nodes in the weak connection area is that the stress of a certain node is 180 MPa, the average stress of its adjacent non-weak connection nodes is 100 MPa, and the stress concentration factor is 1.8. When the preset threshold is set to 1.5, this node is screened and recorded. An example of the recorded coordinate information is that the coordinates corresponding to node N205 are (15.0, 8.7, 0.0), and the stress coefficient is 1.8. An example of the storage of node information is a JSON file, which contains key-value pairs of node numbers, stress coefficients, and coordinates. In the subsequent steps, high-risk nodes are located by reading the JSON file and the stiffness decay rate is adjusted.

[0145] An example of the convergence judgment of the finite element solver is to monitor whether the residual norm is less than a preset tolerance. For example, when the residual is less than 1e-6, it is determined to converge. An example of the update of the real-time stress coefficient is to recalculate the nodal stress and extract the maximum value after the end of each iteration step. An example of the recalculation of the stress concentration coefficient evaluation index is to update the average stress of adjacent nodes and recalculate the ratio after each iteration. An example of the dynamic update of the recorded nodal information is to generate a new CSV file at each iteration, recording the data of high-risk nodes at the current step.

[0146] An example of the verification of the application of external loads is to confirm the rationality of the load direction and magnitude by comparing the experimental strain data with the simulation results. An example of the visualization of the real-time stress coefficient is to render the stress cloud diagram in ParaView and overlay the coordinates of the marked high-risk nodes. An engineering case of threshold setting is that the stress concentration coefficient of a certain aviation composite material structure design requirement should not exceed 1.8, and the threshold is set to 1.8 in the simulation to screen potential failure areas. The recorded nodal information is used to generate a report. For example, the output PDF file contains the positions and stress coefficients of high-risk nodes.

[0147] During the model loading process after dynamic decoupling constraints, if there is an angle between the load direction and the fiber direction, the load needs to be decomposed into components along the fiber and perpendicular to the fiber. For example, when the angle between the load direction and the fiber direction is 30 degrees, the component along the fiber direction is the load multiplied by cos30°, and the perpendicular component is the load multiplied by sin30°. The extraction of the real-time stress coefficient needs to select the corresponding stress component according to the component direction.

[0148] The calculation of the stress concentration coefficient evaluation index needs to exclude the interference of boundary effects. For example, only the internal nodes more than a certain distance from the boundary are calculated. An example of the calculation of the average stress of adjacent non-weak connection regions is to take the arithmetic mean of the stress values of 5 non-weak connection nodes around the node. An example of the recorded nodal information is that the stress coefficient of node N305 is 2.0, and the coordinates are (20.1, 10.5, 0.0). This node is located at the edge of the hole in the laminate and is consistent with the crack initiation position observed in the actual experiment. An example of the handling of threshold overrun is to increase the stiffness decay rate of this node by 50% in subsequent iterations to accelerate the decoupling process.

[0149] S6. Judge whether the real-time stress coefficient exceeds the preset threshold; if it exceeds, iteratively correct the dynamic decoupling constraint based on the stiffness decay rate, including:

[0150] Judge whether the real-time stress coefficient exceeds the preset threshold. If it exceeds, calculate the difference interval between the real-time stress coefficient and the preset threshold. The real-time stress coefficient is the maximum principal stress or Von Mises stress of the nodes in the weak connection area extracted in step S5. The preset threshold is set according to the material interface strength or experimental data. For example, when the critical stress of the composite material interface is 150 MPa, the threshold is set to 120 MPa. The difference interval is the difference range between the real-time stress coefficient and the threshold. For example, the difference of 0 - 10 MPa is interval one, and 10 - 20 MPa is interval two. The interval division is set according to the material property grading.

[0151] Query the stiffness decay rate table according to the difference interval and match the corresponding stiffness decay rate correction step. The stiffness decay rate table is a two-dimensional table. The rows represent the difference intervals, and the columns represent the correction steps. For example, the difference interval of 0 - 10 MPa corresponds to a step of 0.05 / iteration step, and 10 - 20 MPa corresponds to a step of 0.1 / iteration step. The table lookup method automatically matches the interval to which the current difference belongs by writing a script program. For example, when the difference is 15 MPa, it matches interval two and obtains a step of 0.1 / iteration step.

[0152] Update the stiffness decay rate of the dynamic decoupling constraint based on the correction step, and recalculate the real-time stress coefficient of the nodes in the weak connection area. The stiffness decay rate of the dynamic decoupling constraint is adjusted according to the correction step. For example, the original rate is 0.05 / iteration step, and after correction, it is 0.05 + 0.1 = 0.15 / iteration step. The updated rate is applied to the iterative decay of the stiffness value of the nodes in the weak connection area. For example, the initial stiffness of the node is 100 MPa, and after the first iteration at a rate of 0.15 / step, it drops to 85 MPa. Recalculate the real-time stress coefficient by performing a static analysis through a finite element solver and output the updated node stress data.

[0153] An example of the division of the difference interval is to preset three intervals: 0 - 5 MPa, 5 - 10 MPa, and above 10 MPa, corresponding to steps of 0.02, 0.05, and 0.1 respectively. An example of table lookup and matching is that the real-time stress coefficient of a certain node is 130 MPa, the threshold is 120 MPa, and the difference of 10 MPa matches interval two with a step of 0.05. An example of the update of the stiffness decay rate is that the original rate is 0.03 / step, and after correction, it is 0.03 + 0.05 = 0.08 / step, and the node stiffness drops from 80 MPa to 73.6 MPa.

[0154] An example of the process of recalculating the real-time stress coefficient is to modify the node material parameters in Abaqus and then submit the Job to extract the maximum principal stress from the ODB result file. An example of the iteration termination condition is that the real-time stress coefficient is lower than the threshold and the stress change rate between two adjacent iterations is less than 5%. The convergence accuracy is set according to engineering requirements. For example, when the change rate is less than 1%, it is determined to converge.

[0155] The construction of the stiffness decay rate table is based on the experimental data of material interface damage accumulation. For example, the rate corresponding to the difference range of 0 - 5 MPa is 0.02 / step, and the stiffness decay rate in this range is measured through cyclic load tests; an example of the implementation of the look-up table method is that a Python script reads the real-time stress difference, matches the rate table in CSV format, and returns the correction step; an example of the update of dynamic decoupling constraints is to modify the nodal real constants through APDL commands in ANSYS, adjust the stiffness value, and re-solve.

[0156] An example of the recalculated stress coefficient is that after a certain node is updated, the stress drops from 140 MPa to 115 MPa, and the difference of 15 MPa triggers further correction; an example of the iterative process is that after three corrections, the stress coefficient drops to 118 MPa, meeting the threshold of 120 MPa and the change rate is less than 1%, and the iteration is terminated; an example of the dynamic adjustment of the difference range is that when the differences of most nodes are concentrated in a certain range, the range is automatically subdivided to improve the correction accuracy.

[0157] An example of the verification of the stiffness decay rate table is to adjust the range division and step assignment by comparing simulation and experimental data; an example of the update of nodal stiffness is that after two iterations, the stiffness of a certain node drops from 90 MPa to 72.9 MPa, and the real-time stress coefficient drops from 135 MPa to 112 MPa; an example of the data saving after convergence is to store the final nodal stiffness, stress coefficient, and iterative parameters as JSON or CSV files for subsequent analysis and call.

[0158] An example of the stiffness update of dynamic decoupling constraints is to modify the elastic modulus field in the MAT1 card through DMAP language in Nastran, and reduce the modulus value in each iteration; an example of the adaptive adjustment of the difference range is to dynamically optimize the range division according to historical iteration data; an example of the extension of the look-up table logic is to introduce interpolation method to handle non-integer differences. For example, when the difference is 7.5 MPa, the correction value is linearly calculated according to the step sizes of adjacent ranges.

[0159] An engineering case of the iteration termination condition is that for a certain aviation composite material structure, the interface stress is required not to exceed 100 MPa, the threshold is set to 100 MPa, and the convergence accuracy is set to 0.5%; after five iterations, the nodal stress drops from 105 MPa to 99.3 MPa, and the change rate is 0.4%, meeting the termination condition; an example of the output of the corrected macroscopic mechanical property results is to generate a report file containing the final stiffness distribution, stress nephogram, and the number of convergence iterations.

[0160] S7. When the real-time stress coefficient is lower than the preset threshold, output the corrected macroscopic mechanical property results, including:

[0161] Verify that the real-time stress coefficient is lower than the preset threshold and the change rate of the real-time stress coefficient between two adjacent iterations is less than the preset convergence accuracy; the real-time stress coefficient is the nodal stress value after the last iteration in step S6. For example, if the stress of a certain node drops from 150 MPa to 118 MPa after three iterations, and the change rate between two adjacent iterations is 1.5%, and the preset convergence accuracy is 2%, it is determined to converge; the preset threshold is set according to the allowable interface strength of the material. For example, if the critical stress of the composite material interface is 120 MPa, the threshold is set to 120 MPa; the convergence accuracy is calibrated through engineering experience or experimental data. For example, it is required that the stress change rate is less than 1% to ensure the simulation stability.

[0162] Extract the final stiffness value and stress distribution data of the nodes in the weak connection area from the modified non-uniform grid model as key parameters; the non-uniform grid model is the finite element model generated in step S2 and iteratively modified, such as the INP file of Abaqus or the CDB file of ANSYS; the final stiffness value is the elastic modulus or stiffness matrix value of the nodes in the weak connection area after the last iteration. For example, the stiffness of a certain node drops from the initial 100 MPa to 75 MPa; the stress distribution data includes nodal stress components and Von Mises stress values, which are extracted through the finite element post-processing module. For example, it is read from the ODB file of Abaqus.

[0163] Associate the key parameters with the homogenized constitutive equation of the macroscopic model to generate a macroscopic mechanical property report file containing the modified material properties; the homogenized constitutive equation is a macroscopic material model constructed based on the homogenization theory of the composite material microstructure, such as the Hill yield criterion or the anisotropic elastic model; the association method is to input the final stiffness value of the weak connection nodes into the constitutive equation to update the elastic constants or yield parameters in the equation; the report file format is PDF or Excel, containing parameters such as the modified elastic modulus, Poisson's ratio, and yield strength. For example, the elastic modulus is modified from 70 GPa to 65 GPa, and the yield strength is modified from 500 MPa to 480 MPa.

[0164] Update the composite material model data in the finite element simulation database according to the report file and output the updated macroscopic mechanical property results; the finite element simulation database is a local or cloud database storing material parameters, such as MySQL or MongoDB; the update operation is to write the modified material properties into the corresponding fields of the database. For example, replace the elastic modulus field value of the original model; the output results include the updated material parameter list and the corresponding simulation configuration file. For example, the JSON file contains the material name, elastic modulus, Poisson's ratio, and yield strength.

[0165] An example of the verification process is that the real-time stress coefficient of a certain node drops from 130 MPa to 115 MPa, and the adjacent iteration change rate is 0.8%, which is lower than the preset convergence accuracy of 1%, so convergence is determined; an example of extracting key parameters is to export the node stiffness value and stress components from the ODB result file of Abaqus and store them in CSV format; an example of associating constitutive equations is to input the average node stiffness of 65 GPa into the Hill criterion to generate a modified anisotropic parameter table; an example of a report file is a PDF document containing a comparison chart and data table of the elastic modulus before and after modification.

[0166] An example of database update is to upload the modified parameters in JSON format to MongoDB and replace the "Composite_Layer_1" entry of the original model; an example of output result is to generate a new INP file, in which the material card is updated to the modified elastic modulus of 65 GPa and yield strength of 480 MPa; an example of a simulation configuration file is a YAML file that defines loads, boundary conditions, and material parameters for subsequent analysis calls.

[0167] An example of the final stiffness value of a node in the weak connection area is that the elastic modulus of node N1001 drops from 80 GPa to 62 GPa; an example of the application of the homogenized constitutive equation is to calculate the equivalent stiffness matrix of the macroscopic model using the modified elastic modulus of 62 GPa; an example of data representation in the report file is to list the stiffness correction values of all weak connection nodes and their corresponding element numbers.

[0168] An example of composite material model data in the database is entries containing material name, fiber volume fraction, porosity, and modified mechanical parameters; an example of the engineering application of the output result is to submit the updated INP file to the Abaqus solver to perform macroscopic structural strength checking; an example of the verification of the simulation configuration file is to compare the stress nephograms before and after modification to confirm that the interface stress concentration area has been improved.

[0169] An example of a key parameter extraction tool is a Python script to parse the ODB file and extract the stiffness and stress data of specified nodes; an example of the implementation of associating constitutive equations is a MATLAB script to read the stiffness values in the CSV file and call the homogenization algorithm to generate macroscopic parameters; an example of the automated generation of a report file is to use a LaTeX template to integrate data tables and charts into a PDF document.

[0170] An example of permission management for database update is to set multi-level user permissions, allowing only authorized personnel to modify the material parameter fields; an example of the visualization of output results is to render the modified stress nephogram in ParaView and overlay the results of the original model for comparison; an example of the scalability of the simulation configuration file is to support adding new material parameter fields, such as the coefficient of thermal expansion or damping ratio.

[0171] Example 2:Figure 2 The structural schematic diagram of a computer numerical simulation system for a composite material according to the present invention is given. A computer numerical simulation system for a composite material includes:

[0172] Microstructure acquisition module: Obtain the microscopic structure image data of the composite material, extract the fiber arrangement direction and pore distribution parameters, and generate the fiber direction probability density function and pore topological relationship;

[0173] Grid orientation module: Generate a local coordinate system according to the fiber direction probability density function, and perform element division with direction constraints on the finite element grid based on the local coordinate system to obtain a non-uniform grid model;

[0174] Weak area identification module: Based on the pore topological relationship, identify the weak connection areas between adjacent elements in the non-uniform grid model and mark the nodes;

[0175] Dynamic constraint module: Apply dynamic decoupling constraints to the nodes in the weak connection areas. The dynamic decoupling constraints include the initial stiffness and the stiffness decay rate;

[0176] Load application module: Apply an external load to the model after applying the dynamic decoupling constraints, and calculate the real-time stress coefficient of the weak connection areas;

[0177] Iterative correction module: Determine whether the real-time stress coefficient exceeds the preset threshold; if it exceeds, iteratively correct the dynamic decoupling constraints based on the stiffness decay rate;

[0178] Result output module: After the real-time stress coefficient is lower than the preset threshold, output the corrected macroscopic mechanical property results.

[0179] Embodiment 3: A computer numerical simulation device for a composite material. The device includes: a processor, a memory, and a program or instruction stored on the memory and executable on the processor. When the program or instruction is executed by the processor, a computer numerical simulation method for a composite material is implemented.

[0180] The above formulas are all dimensionless and take their numerical calculations. The formulas are obtained by collecting a large amount of data for software simulation to obtain a formula closest to the actual situation. The preset parameters and threshold selection in the formulas are set by those skilled in the art according to the actual situation.

[0181] It should be noted that the present invention can be deployed on the device itself to achieve embedded applications, or can also run on a PC or other terminals with a user interface, so as to meet various hardware environments and usage requirements.

[0182] The above embodiments can be implemented in whole or in part by software, hardware, firmware, or any combination thereof. When implemented using software, the above embodiments can be implemented in whole or in part in the form of a computer program product. The computer program product includes one or more computer instructions or computer programs. When the computer instructions or computer programs are loaded or executed on a computer, the processes or functions described in the embodiments of the present application are generated in whole or in part. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable devices. The computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, the computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wired (such as infrared, wireless, microwave, etc.) means. The computer-readable storage medium can be any available medium that can be accessed by a computer or a data storage device such as a server or data center that contains one or more collections of available media. The available media can be magnetic media (such as floppy disks, hard disks, magnetic tapes), optical media (such as DVDs), or semiconductor media. The semiconductor media can be a solid-state drive.

[0183] Those skilled in the art can clearly understand that for the convenience and conciseness of description, the specific working processes of the systems, devices, and modules described above can refer to the corresponding processes in the foregoing method embodiments and will not be elaborated herein.

[0184] In several embodiments provided in the present application, it should be understood that the disclosed systems, devices, and methods can be implemented in other ways. For example, the device embodiments described above are merely illustrative. For example, the division of the modules is only a logical function division, and there can be other division methods in actual implementation. For example, multiple modules or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the displayed or discussed couplings or direct couplings or communication connections to each other can be through some interfaces, and the indirect couplings or communication connections of the devices or modules can be in electrical, mechanical, or other forms.

[0185] The modules described as separate components may or may not be physically separated, and the components shown as modules may or may not be physical modules. They can be located in one place or distributed to multiple network modules. Some or all of the modules can be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0186] In addition, in each embodiment of the present application, each functional module can be integrated into a processing module, can exist physically alone for each module, or two or more modules can be integrated into one module.

[0187] If the above-mentioned function is implemented in the form of a software functional module and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on such an understanding, the technical solution of the present application, in essence, or the part that contributes to the prior art, or a part of this technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to enable a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in each embodiment of the present application. The aforementioned storage medium includes: various media such as USB flash drives, mobile hard disks, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical discs that can store program codes.

[0188] The above is only the specific implementation manner of the present application, but the protection scope of the present application is not limited thereto. Any person skilled in the art can easily think of changes or substitutions within the technical scope disclosed in the present application, and all of them should be covered by the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

[0189] Finally: The above is only the preferred embodiment of the present invention and is not used to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A computer numerical simulation method for composite materials, characterized in that: The steps include: S1. Obtain the microstructure image data of the composite material, extract the fiber arrangement direction and pore distribution parameters, and generate the fiber direction probability density function and pore topological relationship; S2, generating a local coordinate system according to the fiber direction probability density function, performing direction-constrained unit division on the finite element mesh based on the local coordinate system, and obtaining a non-uniform mesh model; S3. Based on the pore topology, identify the weakly connected areas between adjacent units in the non-uniform grid model and mark the nodes; S4. Apply dynamic decoupling constraints to the nodes in the weakly connected area, where the dynamic decoupling constraints include initial stiffness and stiffness decay rate; S5, applying external loads to the model after applying dynamic decoupling constraints, and calculating the real-time stress coefficient of the weak connection area; S6, determining whether the real-time stress coefficient exceeds a preset threshold; If exceeded, the dynamic decoupling constraint is iteratively corrected based on the stiffness decay rate; S7. When the real-time stress coefficient is lower than the preset threshold, the corrected macroscopic mechanical properties results are output.

2. A computer numerical simulation method for composite materials according to claim 1, characterized in that: S1 includes: S1a, processing microstructure image data by grayscale gradient analysis method, identifying fiber edge contours and calculating the main axis direction angle of each fiber, and statistically analyzing the distribution frequency of the main axis direction angle of the fiber to generate a fiber direction probability density function; S1b, perform binary segmentation on the same image data, extract the pore geometry and spatial position, and construct an adjacency relationship matrix that characterizes the pore topological characteristics based on the spatial position relationship between the pores and fibers; S1c, the fiber direction probability density function is discretized and stored according to the preset angle interval, and combined with the adjacency matrix into a statistical model to characterize the randomness of the microstructure.

3. The computer numerical simulation method for composite materials according to claim 1, characterized in that S2 include: S2a, according to the angle interval distribution in the fiber direction probability density function, a corresponding local coordinate system is assigned to each finite element, and the X-axis direction of the local coordinate system is consistent with the median direction of the fiber main axis direction angle; S2b, based on the X-axis direction of the local coordinate system, adjust the direction of the stiffness matrix of each finite element unit through the rotation matrix so that the main axis direction of the unit stiffness matrix is ​​aligned with the fiber direction; S2c, meshing the elements after adjusting the stiffness matrix to generate a non-uniform mesh model with embedded fiber direction randomness.

4. The computer numerical simulation method for composite materials according to claim 1, characterized in that S3 include: S3a, based on the adjacency relationship matrix in the pore topology, locate the adjacent unit pairs covering the pore area, and calculate the overlap between the interface connection area of ​​the adjacent unit pairs and the pore projection area; S3b, if the overlap exceeds the preset ratio threshold, it is determined that there is a weak connection area between the adjacent unit pairs, and the number and common node coordinates of the unit pair are recorded; S3c. According to the recorded unit pair numbers and node coordinates, the common nodes in the weakly connected area are marked in the non-uniform grid model, and a node marking list is generated.

5. The computer numerical simulation method for composite materials according to claim 1, characterized in that S4 include: S4a, setting the initial binding stiffness according to the material properties of the weakly connected area nodes, the initial binding stiffness being lower than the standard stiffness value of the non-weakly connected area nodes; S4b, based on a preset stiffness decay rate table, determining the stiffness decay rate of the nodes in the weakly connected area, the stiffness decay rate table is dynamically adjusted according to the difference between the real-time stress coefficient and the preset threshold value; S4c, associate the initial binding stiffness with the stiffness decay rate, and impose dynamic decoupling constraints on the nodes in the weakly connected area. The dynamic decoupling constraints are achieved by iteratively updating the node stiffness values.

6. The computer numerical simulation method for composite materials according to claim 1, characterized in that S5 include: S5a, define the external load application area and load direction in the model after dynamic decoupling constraints, and the load direction is consistent with the actual working condition of the composite material; S5b, calculating the real-time stress distribution of nodes in the weakly connected area based on the finite element solver, and extracting the maximum principal stress or VonMises stress as the real-time stress coefficient; S5c, calculating the stress concentration factor evaluation index according to the stress gradient distribution of the nodes in the weak connection area, where the stress concentration factor evaluation index is the ratio of the node stress to the average stress of the adjacent non-weak connection area; S5d, filter out weakly connected area nodes whose stress concentration coefficient evaluation index exceeds the preset threshold, and record their real-time stress coefficient and corresponding coordinate information.

7. The computer numerical simulation method for composite materials according to claim 1, characterized in that S6 include: S6a, determining whether the real-time stress coefficient exceeds a preset threshold, and if so, calculating a difference interval between the real-time stress coefficient and the preset threshold; S6b, querying the stiffness attenuation rate table according to the difference interval, and matching the corresponding stiffness attenuation rate correction step length; S6c, based on the modified step size, the stiffness decay rate of the dynamic decoupling constraint is updated, and the real-time stress coefficients of the nodes in the weakly connected area are recalculated.

8. The computer numerical simulation method for composite materials according to claim 1, characterized in that: S7 includes: S7a, verifying that the real-time stress coefficient is lower than a preset threshold and the change rate of the real-time stress coefficient between two adjacent iterations is less than a preset convergence accuracy; S7b, extracting the final stiffness value and stress distribution data of the nodes in the weakly connected area from the corrected non-uniform grid model as key parameters; S7c, associating the key parameters with the homogenized constitutive equation of the macro model, and generating a macro mechanical properties report file containing the corrected material properties; S7d. Update the composite material model data in the finite element simulation database according to the report file, and output the updated macroscopic mechanical performance results.

9. A computer numerical simulation system for composite materials, used to implement a computer numerical simulation method for composite materials according to any one of claims 1 to 8, characterized in that: include: Microstructure acquisition module: obtains microstructure image data of composite materials, extracts fiber arrangement direction and pore distribution parameters, and generates fiber direction probability density function and pore topological relationship; Mesh orientation module: Generate a local coordinate system according to the fiber direction probability density function, divide the finite element mesh into units with direction constraints based on the local coordinate system, and obtain a non-uniform mesh model; Weak area identification module: Based on the pore topology relationship, it identifies the weak connection areas between adjacent units in the non-uniform grid model and marks the nodes; Dynamic constraint module: dynamic decoupling constraints are imposed on nodes in weakly connected areas. The dynamic decoupling constraints include initial stiffness and stiffness decay rate. Load loading module: applies external loads to the model after dynamic decoupling constraints are applied, and calculates the real-time stress coefficient of the weak connection area; Iterative correction module: determines whether the real-time stress coefficient exceeds the preset threshold; If exceeded, the dynamic decoupling constraint is iteratively corrected based on the stiffness decay rate; Result output module: When the real-time stress coefficient is lower than the preset threshold, the corrected macroscopic mechanical properties results are output.

10. A computer numerical simulation device for composite materials, characterized in that: include: A processor, a memory, and a program or instruction stored in the memory and executable on the processor, wherein when the program or instruction is executed by the processor, a computer numerical simulation method for a composite material as described in any one of claims 1 to 8 is implemented.

Images