Computer numerical simulation method, system and equipment for composite material

Abstract

The invention discloses a computer numerical simulation method, system and equipment for a composite material, particularly relates to the technical field of analog simulation, and is used for solving the problem of simulation distortion caused by neglecting the randomness of a microstructure in the prior art. The method comprises the following steps: acquiring microstructure image data of a composite material, extracting fiber arrangement direction and pore distribution parameters, and generating a fiber direction probability density function and a pore topological relation; constructing a non-uniform grid model constrained by a local coordinate system according to the fiber direction probability density function; identifying a weak connection area based on a pore topological relation and marking nodes; applying dynamic decoupling constraint to the weak connection area, wherein the dynamic decoupling constraint comprises initial rigidity and rigidity attenuation rate; applying an external load to the model and calculating a real-time stress coefficient of the weak connection area; dynamically adjusting the constraint parameter through threshold judgment and iterative correction until the stress coefficient is converged; and finally, outputting a corrected macroscopic mechanical property result.

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 equipment for composite materials. Background Art

[0002] Composite materials are widely used in aerospace, automobile manufacturing and other fields due to their excellent performance. Their design optimization and performance evaluation are highly dependent on computer simulation technology. In the existing technology, multi-scale numerical simulation of composite materials is usually based on the finite element method to build a model, and the microstructural characteristics are mapped to the macroscopic mechanical response through the homogenization theory. However, the microstructure of actual composite materials has natural randomness (such as fiber arrangement, pore distribution, etc.). In order to simplify the calculation model, the existing method generally adopts the homogenization assumption to describe the material properties, which makes it difficult for the simulation process to truly reflect the dynamic influence of microscopic heterogeneity on macroscopic behavior.

[0003] Since existing simulation methods ignore the random characteristics of the microstructure of composite materials, there are systematic deviations between the predicted results of macroscopic performance and the actual material behavior. Especially when analyzing the mechanical response of non-uniformly reinforced composites, the homogenization assumption will mask key mechanisms such as local defect evolution and interface failure, making it impossible for the simulation model to accurately capture the true failure mode of the material under complex loads, thereby 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-mentioned background technology.

[0005] To achieve the above object, the present invention provides the following technical solutions: A computer numerical simulation method for composite materials comprises the following steps: 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, judging whether the real-time stress coefficient exceeds a preset threshold; if so, iteratively correcting the dynamic decoupling constraint 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.

[0006] In a preferred embodiment, S1 comprises: 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.

[0007] In a preferred embodiment, S2 includes: 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.

[0008] In a preferred embodiment, S3 includes: 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.

[0009] In a preferred embodiment, S4 includes: 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.

[0010] In a preferred embodiment, S5 includes: 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 Von Mises 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.

[0011] In a preferred embodiment, S6 includes: 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.

[0012] In a preferred embodiment, 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.

[0013] In another aspect, the present invention provides a computer numerical simulation system for composite materials, comprising: 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 it exceeds, iteratively corrects the dynamic decoupling constraint 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.

[0014] On the other hand, the present invention provides a computer numerical simulation device for composite materials, the device comprising: 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 composite materials is implemented.

[0015] Compared with the prior art, the present invention has the following beneficial effects: 1. Through the extraction of microstructure image data and the construction of probability density functions, the randomness of fiber direction and the spatial distribution characteristics of pores are embedded in the finite element mesh modeling, so that the simulation model can accurately characterize the spatial correlation of microscopic 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 mapping distortion problem between microscopic random characteristics and macroscopic mechanical behaviors, but also provides a modeling basis for multi-scale simulation of composite materials that is closer to the actual working conditions; 2. Through the iterative correction of node stiffness in weakly connected areas and the threshold judgment logic, a dynamic decoupling path based on the real-time stress state is constructed, so that the simulation process can adaptively adjust local constraints and truly reflect the progressive characteristics of interface failure; the dynamic coupling simulation strategy significantly improves the prediction accuracy of damage evolution of composite materials. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 A flowchart of a computer numerical simulation method for a composite material according to the present invention; Figure 2The present invention is a schematic structural diagram of a computer numerical simulation system for composite materials. DETAILED DESCRIPTION

[0017] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0018] Embodiment 1: Figure 1 A computer numerical simulation method for composite materials is provided, which comprises the following steps: 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, judging whether the real-time stress coefficient exceeds a preset threshold; if so, iteratively correcting the dynamic decoupling constraint 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.

[0019] S1. Obtain the microstructure image data of the composite material, extract the fiber arrangement direction and pore distribution parameters, generate the fiber direction probability density function and pore topology relationship, including: The grayscale gradient analysis method is used to process the microstructure image data of the composite material, identify the fiber edge contour and calculate the main axis direction angle of each fiber, and then the distribution frequency of the fiber main axis direction angle is statistically analyzed to generate the fiber direction probability density function. The grayscale gradient analysis method first grays the microstructure image and converts the color image into a grayscale image. The grayscale processing converts the pixel values ​​of the red, green and blue channels into a single grayscale value through the weighted average method.

[0020] The fiber edge contour is determined by calculating the gradient value of each pixel in the grayscale image. The gradient value is calculated by the Sobel operator. The Sobel operator uses a 3×3 convolution kernel to calculate the gradient components of the pixel in the horizontal and vertical directions respectively. The gradient amplitude is the square root of the sum of the squares of the horizontal and vertical components. Pixels with gradient amplitudes greater than a preset threshold are marked as fiber edge points. 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, the least squares method is used to fit the main axis direction angle of each fiber. The least squares method determines the optimal straight line by minimizing the sum of the squares of the vertical distances from the fiber edge points to the fitting straight line. The angle between the fitting straight line and the horizontal axis of the image is the main axis direction angle of the fiber.

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

[0022] The same microstructure image data is binarized and segmented to extract the geometric shape and spatial position of the pores. Based on the spatial position relationship between the pores and fibers, an adjacency relationship matrix characterizing the pore topological characteristics is constructed.

[0023] Binarization segmentation divides the image into pore area and non-pore area by setting a grayscale threshold. The pixel value of the pore area is set to 1, and the pixel value of the non-pore area is set to 0. The grayscale threshold is automatically calculated by the Otsu algorithm or manually set based on experience. Morphological operations (such as erosion and dilation) are used to remove noise points and smooth the pore edges. The erosion operation traverses the image with a 3×3 rectangular structure element, and the dilation operation uses the same structure element to restore the pore shape.

[0024] The geometric shape parameters of the extracted pores include area, perimeter, and aspect ratio. The area is calculated by counting the number of pixels in the pore area, the perimeter is calculated by the chain code method or edge tracking algorithm to calculate the number of pixels on the pore boundary, and the aspect ratio is determined by fitting the ratio of the long side to the short side of the minimum circumscribed rectangle of the pore.

[0025] The spatial position is represented by the centroid coordinates of the pore area, and the centroid coordinates are obtained by calculating the average values ​​of the horizontal and vertical coordinates of all pixels in the pore area; an adjacency matrix is ​​constructed based on the spatial position relationship between the pores and fibers, and the rows and columns of the adjacency matrix represent the numbers of the pores and fibers, respectively. The matrix elements are Boolean values ​​or weight values. The Boolean value indicates whether the pores and fibers are adjacent. Adjacency is defined as the minimum distance from the edge of the pore to the edge of the fiber is less than a preset threshold (for example, 5 pixels). The weight value is calculated based on the overlapping area or edge distance between the pores and fibers, such as the ratio of the overlapping area to the total pore area or the normalized Euclidean distance.

[0026] The adjacency matrix is ​​stored in the form of a two-dimensional array or a sparse matrix. The sparse matrix stores non-zero elements through triplets (row number, column number, element value) to save memory space.

[0027] The fiber direction probability density function is discretized and stored according to preset angle intervals and combined with the adjacency matrix into a statistical model that characterizes the randomness of the microstructure.

[0028] The fiber direction probability density function is discretized into a probability value sequence according to a preset angle interval (such as 5 degrees or 10 degrees), each interval corresponds to a probability value, and the discretized storage data structure is an array or list, and the array index corresponds to the starting angle of the angle interval.

[0029] The adjacency matrix is ​​stored in the form of a two-dimensional array or a sparse matrix, and the row and column numbers of the matrix correspond to the unique identifiers of the pores and fibers. When the fiber direction probability density function and the adjacency matrix are combined into a statistical model, the data structure of the statistical model is a dictionary or a JSON file, and 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 matrix.

[0030] During the merging process, the fiber direction probability density function and the adjacency matrix data are associated through unique identifiers. The unique identifiers are the fiber and pore numbers, which are generated through an image coordinate hash algorithm, such as converting the centroid coordinates into a string and taking the first 8 bits of the MD5 hash value.

[0031] The statistical model is used to describe the random characteristics of the microstructure of the composite material, such as the distribution law of the fiber direction and the spatial adjacency relationship between pores and fibers. The subsequent steps are to perform finite element meshing and mechanical simulation by reading the data in the statistical model.

[0032] In the process of fiber edge contour recognition, the Sobel operator is used to calculate the gradient amplitude. 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 grayscale image is convolved to obtain the horizontal gradient component and the vertical gradient component respectively. In the process of fiber main axis angle fitting, the least squares method is used to calculate the slope and intercept of the fitting line by solving the linear equation group, and the slope is the tangent value of the fiber main axis angle. The normalization of the fiber direction probability density function is to divide the frequency value of each angle interval by the total number of fibers. The total number of fibers is obtained by the number of independent fiber regions in the image through connected region analysis.

[0033] In the process of extracting the pore geometry, morphological operations are used to remove noise points, the erosion operation is achieved by taking the minimum pixel value in the area covered by the structural element, and the dilation operation is achieved by taking the maximum value.

[0034] The calculation formula of the centroid coordinates is: the centroid abscissa is equal to the sum of the abscissas of all pixels in the pore area divided by the number of pixels, and the same is true for the centroid ordinate.

[0035] In the process of constructing the adjacency relationship matrix, the minimum distance between the pore and the fiber is calculated by traversing the coordinates of the pore edge pixel points and the fiber edge pixel points to calculate the Euclidean distance, and the minimum value of all distances is taken as the distance between the two; the weight value is calculated, for example, as 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 pixel-level logical and operation to calculate the number of pixels where both are 1.

[0036] In the discretely stored fiber orientation probability density function, the angle interval division needs to cover the range of 0 to 180 degrees to avoid periodic repetition of the orientation angle; in the process of merging the adjacency matrix, the generation of the unique identifier needs to ensure that the fiber and pore numbers are unique and can be generated repeatedly, and the hash algorithm is implemented by converting the floating point number of the centroid coordinates into a string and then taking the hash value; the data structure example of the statistical model is in JSON format, including the "fiber_orientation" key corresponding to the discretized probability value array, and the "pore_adjacency" key corresponding to the sparse storage format of the adjacency matrix; in the subsequent steps, by parsing the data structure in the statistical model, the fiber orientation probability density function is used for the material property assignment of the finite element unit, and the adjacency matrix is ​​used to identify weakly connected areas and impose dynamic decoupling constraints.

[0037] In the analysis of the spatial position relationship between pores and fibers, if the minimum distance between the pores and fibers is less than a preset threshold, the two are judged to be 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 to represent 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. Boolean values ​​are suitable for qualitative analysis, and continuous values ​​are suitable for quantitative simulation. The discrete probability value sequence in the statistical model can be converted into a continuous probability density function through interpolation methods, such as linear interpolation or cubic spline interpolation. The choice of interpolation method is determined according to the simulation accuracy requirements.

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

[0039] An example of sparse matrix storage of the adjacency matrix is ​​a triple (row number, column number, element value) that records only non-zero elements. For example, (1,5,0.8) means that the overlapping area between the pore numbered 1 and the fiber numbered 5 is 0.8.

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

[0041] S2. Generate a local coordinate system according to the fiber direction probability density function, and divide the finite element mesh into units with direction constraints based on the local coordinate system to obtain a non-uniform mesh model, including: According to the angle interval distribution in the fiber direction probability density function, a corresponding local coordinate system is assigned to each finite element unit; the angle interval distribution of the fiber direction probability density function is a preset discretized angle interval, such as 0-5 degrees, 5-10 degrees, etc., and each angle interval corresponds to a probability value; the X-axis direction of the local coordinate system is consistent with the median direction of the fiber main axis direction angle in the angle interval, and the median direction is the middle angle value of the angle interval, for example, the median direction of the 0-5 degree interval is 2.5 degrees.

[0042] The local coordinate system of each finite element unit is determined according to the median direction of the angle interval to which it belongs. For example, if the angle interval of a unit is 10-15 degrees, 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.

[0043] Based on the X-axis direction of the local coordinate system, the stiffness matrix direction of each finite element unit is adjusted through the rotation matrix; the rotation matrix is ​​used to transform the unit 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.

[0044] For example, if the angle between the X-axis 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.

[0045] 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 principal axis direction of the fiber.

[0046] The adjustment of the unit stiffness matrix is ​​achieved through matrix multiplication, that is, the stiffness matrix in the local coordinate system is obtained by multiplying the rotation matrix with the original stiffness matrix.

[0047] The elements after adjusting the stiffness matrix are meshed to generate a non-uniform mesh model with embedded randomness in the fiber direction; the meshing is based on a mesh generation module of the finite element software, such as using hexahedral or tetrahedral elements to discretize the composite material geometry model.

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

[0049] The unit direction distribution in the non-uniform grid model matches the statistical law of the fiber direction probability density function. For example, the higher the probability value corresponding to a certain angle interval, the greater the proportion of units in this direction. After the grid division is completed, each unit in the non-uniform grid model contains the adjusted stiffness matrix direction information for subsequent multi-scale mechanical simulation.

[0050] In the process of assigning the local coordinate system, if the discretization angle interval of the fiber direction probability density function is 5 degrees, the median direction of each interval is the starting angle of the interval plus 2.5 degrees.

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

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

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

[0054] 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 unit is 30 degrees, then the main direction of its stiffness matrix and the X-axis direction of the global coordinate system are at an angle of 30 degrees.

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

[0056] When the non-uniform grid model is used for mechanical simulation in subsequent steps, the external load boundary conditions are set according to the actual working conditions. For example, the 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.

[0057] S3. Based on the pore topology, identify the weakly connected areas between adjacent units in the non-uniform grid model and mark the nodes, including: Based on the adjacency relationship matrix in the pore topology, the adjacent unit pairs covering the pore area are located, and the overlap between the interface connection area of ​​the adjacent unit pairs and the pore projection area is calculated.

[0058] The rows and columns of the adjacency matrix represent the numbers of pores and fibers, respectively. The matrix elements are Boolean values ​​or weight values. When the Boolean value is 1, it means that the pores are adjacent to the fibers. The adjacent unit pair covering the pore area refers to two units that share a boundary in the finite element mesh, and one unit belongs to the pore area and the other unit belongs to the fiber area.

[0059] The interface connection area is the geometric area of ​​the shared boundary of 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 area on the shared boundary plane, which is calculated by counting the number of pore pixels in the projection plane and multiplying the actual pixel size; the overlap is the ratio of the interface connection area to the pore projection area. For example, if the interface connection area is 10 square microns and the pore projection area is 15 square microns, the overlap is 66.7%.

[0060] If the overlap exceeds the preset ratio threshold, it is determined that there is a weak connection area between adjacent unit pairs, and the number and common node coordinates of the unit pair are recorded; the preset ratio threshold is set according to the material interface strength experimental data, for example, when the overlap exceeds 50%, it is determined to be a weak connection area.

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

[0062] 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 to generate a node marking list; the marking method is to traverse the recorded unit pair numbers, locate the corresponding units and their common nodes, and add the node identifiers to the node marking list.

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

[0064] After the node marking list is generated, the weakly connected area nodes in the non-uniform mesh model are marked with special attributes, such as adding a "weakly connected" mark in the node attribute field; the node marking list is used to apply dynamic decoupling constraints in subsequent steps, such as adjusting stiffness parameters only for marked nodes.

[0065] In the calculation example of the overlap, if the interface connection area is 8 square microns and the pore projection area is 12 square microns, the overlap is 66.7%; the preset ratio threshold is set according to the material type, for example, brittle material is set to 40%, and ductile material is set to 60%; the example of recording the unit pair number is that unit number 5 and unit number 12 constitute an adjacent unit pair, and their common nodes are node 7, node 8, and node 9; the 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); the example of generating the node mark list is to add nodes 7, nodes 8, and nodes 9 to the list, and mark their coordinates and the unit pair numbers to which they belong.

[0066] In the process of determining the weak connection area, if the overlap is lower than the preset threshold, the adjacent unit pair is ignored. For example, when the overlap is 30%, it is determined to be a normal connection area. The storage format of the node marking list is compatible with the finite element software, such as the node set (Node Set) of Abaqus or the named selection set (Named Selection) of ANSYS. An example of subsequent processing of marked nodes is to reduce their binding stiffness or introduce contact failure criteria in finite element analysis. The node marking information in the non-uniform mesh model is realized by modifying the mesh file, such as adding attribute fields in the node definition part of the INP file.

[0067] The calculation of overlap is based on actual geometric projection, for example, projecting the three-dimensional shape of the pore area onto the shared boundary plane of adjacent unit pairs, 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 unit, for example, the single-layer thickness of the composite laminate is 0.1 mm; an example of statistics of the pore projection area is to use image processing software (such as ImageJ) to measure the number of pixels in the projection area and convert it into actual size; an example of recording unit pair numbers is to store the numbers of adjacent unit pairs as a two-dimensional array, with the first column being the pore unit number and the second column being the fiber unit number.

[0068] An example of the application of the node marking list is to read the list data in the finite element solver and apply dynamic decoupling constraints to the marked nodes; an example of parameter adjustment of the dynamic decoupling constraint is to modify the contact properties of the corresponding nodes in the solver input file according to the node number in the node marking list; the node marking information in the non-uniform mesh model is retained until the simulation is completed and used for subsequent result analysis, such as extracting the stress concentration factor in the weak connection area; the verification method of the node marking list is visual inspection, such as highlighting the marked node in the finite element post-processing software to confirm that it is located at the interface area between the pore and the fiber.

[0069] The basis for setting the overlap threshold includes material interface peel strength test data, such as the critical interface peel force measured by double cantilever beam test and converted into an area ratio threshold; the global uniqueness of the unit pair number is guaranteed by the finite element mesh generation algorithm, such as assigning an increasing integer number to each unit during meshing; the accuracy of the node coordinates is consistent with the meshing resolution, such as when the mesh size is 1 micron, the node coordinates are retained to three decimal places; the efficiency of generating the node mark list is achieved through batch processing, such as traversing all recorded unit pair numbers, extracting all common nodes at one time and removing duplicates.

[0070] An example of marking weakly connected areas is to identify high overlap areas 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 for a node marking list is a CSV file, which contains the node number, X-coordinate, Y-coordinate, Z-coordinate, and the unit pair number to which it belongs. 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. An example of simulation results for a non-uniform mesh model shows that the stress concentration factor in the weakly connected area is significantly higher than the predicted value of the homogenized model, which verifies the effectiveness of the marking logic.

[0071] The selection of projection plane in the overlap calculation is determined according to the fiber direction, for example, projection is performed along the main axis of the fiber to reflect the actual interface contact situation; the calculation example of the interface connection area is the multiplication of the shared boundary length of the rectangular unit by the unit thickness, and the multiplication of the shared boundary arc length of the circular unit by the thickness; the measurement example of the pore projection area is to use the metallographic microscope image and calculate the pore proportion in the projection area through image analysis software; the example of the recording format of the unit pair number is a text file, in which each line records a pair of unit numbers and their overlap values; the example of the tool for generating the node mark list is a Python script, which parses the unit pair records and extracts the node coordinates, and outputs them as a node set file.

[0072] After the nodes in the weakly connected area are marked, the marked nodes in the non-uniform mesh model participate in the dynamic decoupling constraint calculation in the subsequent simulation; an example of 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 maintenance of the node marking list is to synchronously update the node number and coordinates when the mesh is encrypted or coarsened to ensure the consistency of the marking information; an example of simulation output of the non-uniform mesh model includes the stress history data of the marked nodes, which is used to analyze the damage evolution law of the weakly connected area.

[0073] S4. Dynamic decoupling constraints are applied to the nodes in the weakly connected area. The dynamic decoupling constraints include initial stiffness and stiffness decay rate, including: The initial binding stiffness is set according to the material properties of the nodes in the weakly connected area. The material properties include the interface strength parameters between the fiber and the matrix, and the interface strength parameters are obtained through a single fiber pull-out test or a droplet debonding test. 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-weakly connected area, such as the stiffness value corresponding to the elastic modulus of 100 GPa. The initial binding stiffness is set by modifying the node stiffness parameters in the finite element software, such as reducing the elastic modulus field value of the weakly connected node to a preset value.

[0074] The stiffness decay rate of the weakly connected area nodes is determined based on the preset stiffness decay rate table. The stiffness decay rate table is a two-dimensional table and the rows represent the difference interval between the real-time stress coefficient and the preset threshold, and the columns represent the corresponding stiffness decay rate. The difference interval is calibrated by material fatigue test data. For example, the difference of 5MPa corresponds to an interface damage accumulation rate of 0.1. The dynamic adjustment method of the stiffness decay rate is a table lookup method and matches the corresponding rate value according to the stress difference of the current iteration step. For example, when the real-time stress coefficient exceeds the threshold of 5MPa, the stiffness decay rate is 0.1 / iteration step.

[0075] The initial binding stiffness is associated with the stiffness decay rate and a dynamic decoupling constraint is applied to the nodes in the weakly connected 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 implemented by iteratively updating the node stiffness value and multiplying the node stiffness value by (1-stiffness decay rate) in each iteration step; for example, the stiffness value drops to 45 MPa after the first iteration and to 40.5 MPa in the second iteration. 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.

[0076] The stiffness decay rate table is constructed based on the experimental data of material interface damage accumulation. For example, the stiffness decay rate under different stress differences is measured through cyclic load tests. An implementation example of the table lookup method is to write a script program to read the real-time stress difference and match the rate value of the corresponding interval in the rate table. An example of the termination condition of the iterative update is that the node stiffness value drops below a threshold or the number of iteration 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 is associated with the number of iteration steps.

[0077] An example of the logic of applying dynamic decoupling constraints is to read the real-time stress coefficient in each iteration step of finite element analysis, calculate the difference with the preset threshold, look up the table to obtain the stiffness decay rate 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 next iteration step calculation. After the iteration termination condition is met, the final node stiffness value and the corresponding stress distribution result are output. An example of modifying the node stiffness value is to use the USDFLD subroutine in Abaqus to define the field variables and update the material parameters according to the number of iteration steps.

[0078] During the stiffness update process of the nodes in the weakly connected area, if the real-time stress difference decreases to a lower range, 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 is switched 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 elastic modulus of the matrix is ​​30 GPa, the node stiffness value must not be lower than 30 GPa. Examples of the control of the iteration step size are fixed step size, that is, the step size is consistent for each iteration, or 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 the damage evolution.

[0079] The test case of the dynamic decoupling constraint is to simulate the interface delamination of composite laminates and verify the correspondence between the stiffness attenuation of the marked nodes and the delamination expansion. An example of storing the node stiffness value is to record the change of the stiffness value of each iteration 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 simulation of interface performance degradation under cyclic loads through node stiffness iteration. An example of visualizing the node stiffness value is to render the stiffness distribution cloud diagram of different iteration steps in the finite element post-processing software.

[0080] An example of dynamic adjustment of the stiffness decay rate table is to estimate the rate value by linear interpolation in the uncovered interval, for example, when the difference is 7.5MPa, the linear interpolation is 0.15 / step based on the rates of 5MPa (0.1 / step) and 10MPa (0.2 / step). An example of iterative update of node stiffness value is to write an iterative loop in Nastran using DMAP language and update the stiffness parameters and resubmit the solution in each loop. An example of parallel calculation of dynamic decoupling constraints is to distribute the stiffness update tasks of different nodes to multi-core processors.

[0081] After the dynamic decoupling constraints of the nodes in the weakly connected area are applied, the model has time-varying stiffness characteristics and can simulate the progressive failure process of interface damage. The storage and call examples of node stiffness values ​​are to store the stiffness values ​​as time history variables for the post-processing module to generate a stiffness attenuation curve. The convergence judgment example of the dynamic decoupling constraints is to monitor the change rate of the stiffness value and judge convergence if the change rate is less than 1% for three consecutive iterations. The maintenance example of the node mark list is to synchronously update the node number and coordinates when the mesh is encrypted or coarsened.

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

[0083] S5. Apply external loads to the model after dynamic decoupling constraints are applied, and calculate the real-time stress coefficient of the weak connection area, including: In the model after dynamic decoupling constraints, the external load application area and load direction are defined, and the load direction is consistent with the actual working condition of the composite material; the external load application area is determined according to the force position of the actual engineering structure, such as the edge of the composite laminate or the area around the hole; the load direction is set according to the working condition type, such as 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 node or unit in the finite element software and apply force or displacement boundary conditions, such as defining concentrated force or pressure distribution through the Load module in Abaqus.

[0084] The real-time stress distribution of nodes in the weakly connected area is calculated based on the finite element solver, and the maximum principal stress or Von Mises stress is extracted 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 node 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 node stress value, for example, the maximum principal stress of a node is 200MPa or the Von Mises stress is 150MPa.

[0085] The stress concentration factor evaluation index is calculated based on the stress gradient distribution of the nodes in the weakly connected area. The stress gradient distribution is calculated by the stress difference of adjacent nodes. For example, if the stress of a node is 200 MPa, the average stress of the adjacent non-weakly connected area nodes is 120 MPa. The stress concentration factor evaluation index is the ratio of the node stress to the average stress of the adjacent non-weakly connected area, for example, 200 / 120≈1.67. The adjacent non-weakly connected area is defined as a set of nodes that are directly connected to the weakly connected node and are not marked as weakly connected.

[0086] Nodes in weakly connected areas whose stress concentration coefficient evaluation index exceeds the preset threshold are screened, and their real-time stress coefficient and corresponding coordinate information are recorded; the preset threshold is set according to the stress concentration coefficient allowed by the material, for example, the maximum stress concentration coefficient allowed for the composite material interface is 1.5, and nodes exceeding this value are judged to be 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.

[0087] An example of an external load application area is a unilateral tensile condition of a composite laminate, where the load application area is one end of the plate and the load direction is along the length of the plate. An example of a finite element solver calculation is to use the Abaqus / Standard solver to perform a static analysis and output a node stress cloud diagram. An example of extracting the maximum principal stress is to read the stress components from the Abaqus ODB result file and calculate the maximum principal stress value. An example of extracting Von Mises stress is to directly read the Mises stress field data.

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

[0089] The preset threshold is set based on material fatigue test data, such as the critical stress concentration factor of interface failure measured by cyclic load test. 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 to visualize the location of high-risk nodes, such as highlighting nodes near coordinates (10.2, 5.3, 0.0) in finite element post-processing software. An example of the application of node information is to give priority to high-risk nodes in subsequent iterative corrections and adjust their dynamic decoupling constraint parameters in a targeted manner.

[0090] An example of model loading after dynamic decoupling constraints is to apply force loads through the Solution module in ANSYS, and the direction is consistent with the fiber arrangement direction; an example of obtaining real-time stress distribution is to print node stress results through the PRNSOL command in Mechanical APDL; an example of the calculation tool for stress concentration factor evaluation indicators is a Python script that automatically traverses node data and outputs a ratio list; an example of threshold screening is to use the Pandas library to filter records in CSV files with stress concentration factors greater than 1.5.

[0091] An example of stress gradient analysis of nodes in weakly connected areas is that the stress of a node is 180MPa, the average stress of its adjacent non-weakly connected nodes is 100MPa, and the stress concentration factor is 1.8; when the preset threshold is set to 1.5, the node is screened and recorded; an example of recorded coordinate information is the coordinates corresponding to node N205 (15.0, 8.7, 0.0), and the stress coefficient is 1.8; an example of storage of node information is a JSON file, which contains node numbers, stress coefficients, and coordinate key-value pairs; the subsequent steps locate high-risk nodes and adjust the stiffness decay rate by reading the JSON file.

[0092] An example of convergence judgment of the finite element solver is to monitor whether the residual norm is less than the preset tolerance, for example, convergence is determined when the residual is less than 1e-6; an example of updating the real-time stress coefficient is to recalculate the node stress and extract the maximum value after each iteration step; an example of recalculating the stress concentration factor evaluation index is to update the average stress of adjacent nodes and recalculate the ratio after each iteration; an example of dynamic update of recorded node information is to generate a new CSV file for each iteration to record the high-risk node data of the current step.

[0093] An example of verifying 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 visualizing the real-time stress coefficient is to render a stress cloud diagram in ParaView and superimpose the coordinates of the marked high-risk nodes; an engineering case for setting the threshold is that the design of an aviation composite structure requires the stress concentration factor not to exceed 1.8, and the threshold is set to 1.8 in the simulation to screen potential failure areas; the recorded node information is used to generate reports, such as the output PDF file containing the location and stress coefficient of the high-risk nodes.

[0094] During the loading process of the model after dynamic decoupling constraint, 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 real-time stress coefficient requires selecting the corresponding stress component according to the component direction.

[0095] The calculation of the stress concentration factor evaluation index needs to exclude the interference of boundary effects, for example, only the internal nodes that are more than a certain distance away from the boundary are calculated; the average stress calculation example of adjacent non-weakly connected areas is to select the stress values ​​of 5 non-weakly connected nodes around the node and take the arithmetic average; the recorded node information example is the stress coefficient of node N305 is 2.0, the coordinates are (20.1, 10.5, 0.0), and this node is located at the edge of the hole in the laminate, which is consistent with the crack initiation position observed in the actual experiment; the processing example of threshold exceeding the limit is to increase the stiffness decay rate of the node by 50% in subsequent iterations to accelerate the decoupling process.

[0096] S6. Determine whether the real-time stress coefficient exceeds a preset threshold; if so, iteratively correct the dynamic decoupling constraint based on the stiffness decay rate, including: It is determined whether the real-time stress coefficient exceeds a preset threshold value. If it exceeds, the difference interval between the real-time stress coefficient and the preset threshold value is calculated; the real-time stress coefficient is the maximum principal stress or Von Mises stress of the weakly connected region node extracted in step S5, and the preset threshold value 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 value is set to 120 MPa; the difference interval is the difference range between the real-time stress coefficient and the threshold value, for example, the difference of 0-10 MPa is interval one, and 10-20 MPa is interval two, and the interval division is set according to the material performance classification.

[0097] According to the difference interval, the stiffness decay rate table is queried to match the corresponding stiffness decay rate correction step; the stiffness decay rate table is a two-dimensional table, the row represents the difference interval, and the column represents the correction step. For example, the difference interval 0-10MPa corresponds to a step of 0.05 / iteration step, and 10-20MPa 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 15MPa, interval two is matched and a step of 0.1 / iteration step is obtained.

[0098] The stiffness decay rate of the dynamic decoupling constraint is updated based on the revised step size, and the real-time stress coefficient of the nodes in the weakly connected area is recalculated; the stiffness decay rate of the dynamic decoupling constraint is adjusted according to the revised step size, for example, the original rate is 0.05 / iteration step, and the revised rate 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 weakly connected area, for example, the initial stiffness of the node is 100MPa, and it drops to 85MPa after the first iteration at a rate of 0.15 / step; the real-time stress coefficient is recalculated, and static analysis is performed through the finite element solver to output the updated node stress data.

[0099] An example of the division of the difference interval is to preset three intervals: 0-5MPa, 5-10MPa, and above 10MPa, corresponding to step sizes of 0.02, 0.05, and 0.1 respectively; an example of table matching is that the real-time stress coefficient of a node is 130MPa, the threshold is 120MPa, the difference is 10MPa, and the matching interval is two, with a step size of 0.05; an example of updating the stiffness attenuation rate is the original rate of 0.03 / step, which is corrected to 0.03+0.05=0.08 / step, and the node stiffness is reduced from 80MPa to 73.6MPa.

[0100] 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 in 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, convergence is determined when the change rate is less than 1%.

[0101] The stiffness decay rate table is constructed based on the experimental data of material interface damage accumulation. For example, the difference range of 0-5MPa corresponds to a rate of 0.02 / step, and the stiffness decay rate in this range is measured through a cyclic load test. The implementation example of the table lookup method is a Python script that reads the real-time stress difference, matches the rate table in CSV format and returns the corrected step size. The update example of the dynamic decoupling constraint is to modify the node real constants in ANSYS through the APDL command loop, adjust the stiffness value and re-solve.

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

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

[0104] An example of the stiffness update of the dynamic decoupling constraint is to modify the elastic modulus field in the MAT1 card through the DMAP language in Nastran, and reduce the modulus value in each iteration; an example of the adaptive adjustment of the difference interval is to dynamically optimize the interval division according to the historical iteration data; an example of the expansion of the table lookup logic is to introduce the interpolation method to handle non-integer differences, for example, when the difference is 7.5MPa, the correction value is linearly calculated according to the step size of the adjacent interval.

[0105] An engineering case of the iteration termination condition is that a certain aviation composite material structure requires the interface stress not to exceed 100MPa, the threshold is set to 100MPa, and the convergence accuracy is set to 0.5%; after five iterations, the node stress drops from 105MPa to 99.3MPa, with a change rate of 0.4%, which meets the termination condition; an example of outputting the corrected macroscopic mechanical properties results is to generate a report file containing the final stiffness distribution, stress cloud diagram and the number of convergence iterations.

[0106] S7. When the real-time stress coefficient is lower than the preset threshold, the corrected macroscopic mechanical properties results are output, including: Verify that the real-time stress coefficient is lower than the preset threshold and the rate of change of the real-time stress coefficient between two adjacent iterations is less than the preset convergence accuracy; the real-time stress coefficient is the node stress value after the last iteration in step S6, for example, the stress of a node drops from 150MPa to 118MPa after three iterations, and the rate of change between two adjacent iterations is 1.5%. Convergence is determined when the preset convergence accuracy is 2%; the preset threshold is set according to the allowable interface strength of the material, for example, the critical stress of the composite material interface is 120MPa, and the threshold is set to 120MPa; the convergence accuracy is calibrated by engineering experience or experimental data, for example, the stress change rate is required to be less than 1% to ensure simulation stability.

[0107] The final stiffness value and stress distribution data of the weakly connected area nodes are extracted from the corrected non-uniform grid model as key parameters; the non-uniform grid model is the finite element model generated in step S2 and iteratively corrected, 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 weakly connected area nodes after the last iteration, for example, the stiffness of a certain node drops from the initial 100MPa to 75MPa; the stress distribution data includes the node stress components and the Von Mises stress values, which are extracted through the finite element post-processing module, for example, read from the ODB file of Abaqus.

[0108] The key parameters are associated with the homogenized constitutive equation of the macro model to generate a macro mechanical properties report file containing the corrected material properties; the homogenized constitutive equation is a macro material model constructed based on the theory of homogenization of the microstructure of composite materials, such as the Hill yield criterion or the anisotropic elastic model; the association method is to input the final stiffness value of the weakly connected node into the constitutive equation and update the elastic constant or yield parameter in the equation; the report file format is PDF or Excel, including corrected parameters such as elastic modulus, Poisson's ratio and yield strength, for example, the elastic modulus is corrected from 70GPa to 65GPa, and the yield strength is corrected from 500MPa to 480MPa.

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

[0110] An example of the verification process is that the real-time stress coefficient of a node drops from 130MPa to 115MPa, and the adjacent iteration change rate is 0.8%, which is lower than the preset convergence accuracy of 1%, and convergence is determined; an example of extracting key parameters is exporting the node stiffness value and stress component from the Abaqus ODB result file and storing them in CSV format; an example of associating the constitutive equation is inputting the node stiffness mean value of 65GPa into the Hill criterion to generate a corrected anisotropic parameter table; an example of the report file is a PDF document containing a comparison chart of the elastic modulus before and after correction and a data table.

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

[0112] An example of the final stiffness value of the nodes in the weakly connected area is that the elastic modulus of node N1001 is reduced from 80 GPa to 62 GPa. An example of the application of the homogenized constitutive equation is to use the corrected elastic modulus of 62 GPa to calculate the equivalent stiffness matrix of the macro model. An example of the data table in the report file is to list the stiffness correction values ​​of all weakly connected nodes and the corresponding unit numbers.

[0113] An example of composite material model data in the database is an entry containing material name, fiber volume fraction, porosity and corrected mechanical parameters; an example of engineering application of output results is submitting the updated INP file to the Abaqus solver to perform macroscopic structural strength verification; an example of verification of the simulation configuration file is comparing the stress cloud diagrams before and after correction to confirm that the interface stress concentration area has been improved.

[0114] An example of a key parameter extraction tool is a Python script that parses ODB files and extracts stiffness and stress data for specified nodes. An example of an implementation of the associated constitutive equation is a MATLAB script that reads stiffness values ​​from a CSV file and calls a homogenization algorithm to generate macro parameters. An example of automated report file generation is the use of a LaTeX template to integrate data tables and charts into a PDF document.

[0115] An example of permission management for database updates is setting up multi-level user permissions to allow only authorized personnel to modify material parameter fields. An example of visualization of output results is rendering the corrected stress cloud in ParaView and overlaying the original model results to compare the differences. An example of the extensibility of simulation configuration files is supporting the addition of new material parameter fields, such as thermal expansion coefficient or damping ratio.

[0116] Embodiment 2: Figure 2 A structural schematic diagram of a computer numerical simulation system for composite materials of the present invention is given, and a computer numerical simulation system for composite materials comprises: 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 it exceeds, iteratively corrects the dynamic decoupling constraint 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.

[0117] Embodiment 3: A computer numerical simulation device for composite materials, the device comprising: 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 composite materials is implemented.

[0118] The above formulas are all dimensionless and numerical calculations. The formula is a formula that is closest to the actual situation obtained by collecting a large amount of data and performing software simulation. The preset parameters and thresholds in the formula are set by technicians in this field according to actual conditions.

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

[0120] The above embodiments may be implemented in whole or in part by software, hardware, firmware or any other combination thereof. When implemented by software, the above embodiments may 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 process or function described in the embodiment of the present application is generated in whole or in part. The computer may be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions may be stored in a computer-readable storage medium, or may be transmitted from one computer-readable storage medium to another computer-readable storage medium. For example, the computer instructions may be transmitted from one website, computer, server or data center to another website, computer, server or data center by wired (e.g., infrared, wireless, microwave, etc.) means. The computer-readable storage medium may be any available medium that a computer can access or a data storage device such as a server or data center that contains one or more available media sets. The available medium may be a magnetic medium (e.g., a floppy disk, a hard disk, a tape), an optical medium (e.g., a DVD), or a semiconductor medium. The semiconductor medium may be a solid-state hard disk.

[0121] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working processes of the systems, devices and modules described above can refer to the corresponding processes in the aforementioned method embodiments and will not be repeated here.

[0122] In the 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 only schematic. For example, the division of the modules is only a logical function division. There may be other division methods in actual implementation, such as 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 mutual coupling or direct coupling or communication connection shown or discussed can be through some interfaces, indirect coupling or communication connection of devices or modules, which can be electrical, mechanical or other forms.

[0123] 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, and may be located in one place or distributed on multiple network modules. Some or all of the modules may be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0124] In addition, each functional module in each embodiment of the present application may be integrated into one processing module, or each module may exist physically separately, or two or more modules may be integrated into one module.

[0125] If the functions are implemented in the form of software function modules and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present application, or the part that contributes to the prior art or the part of the technical solution, can be embodied in the form of a software product, which is stored in a storage medium and includes several instructions for a computer device (which can be a personal computer, server or network device, etc.) to perform all or part of the steps of the methods described in each embodiment of the present application. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM), random access memory (RAM), disk or optical disk, and other media that can store program codes.

[0126] The above is only a specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any person skilled in the art who is familiar with the present technical field can easily think of changes or substitutions within the technical scope disclosed in the present application, which should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope of the claims.

[0127] Finally: The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should 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