Multi-material gradient interface microscopic crack path prediction method

Abstract

A multi-material gradient interface microscopic crack path prediction method comprises the following steps: S1, according to material similarity or geometric adjacency, orienting to form an edge, and generating an interface diagram; s2, extracting an energy release rate, a principal stress gradient and a temperature; s3, the three field quantities are mapped back to the graph through KD-Tree nearest neighbor and differentiable Gaussian interpolation, and node and edge embedding features are obtained; s4, the temperature-stress-energy cascade gating message passing network calculates a crack initiation probability and an edge probability, and extracts a crack path according to the maximum cumulative probability; s5, the path is gridded into an indication tensor, and an error gradient for the probability and the energy density is generated; s6, synchronously correcting the network weight and the local energy density according to the gradient, and circulating the steps S2 to S5 until the score is converged; s7, freezing the weight, screening crack initiation nodes, performing spline smoothing on a path, formatting an energy field and a reliability label, and outputting the energy field and the reliability label; according to the method, high-precision failure evaluation is obtained with less calculation amount, and the crack initiation position, the three-dimensional crack curve and the reliability level are output.

Description

Technical Field

[0001] This invention relates to the field of thermo-mechanical coupling failure analysis technology for composite materials, and in particular to a method for predicting microcrack paths at multi-material gradient interfaces. Background Technology

[0002] In the field of predicting crack initiation and propagation under cyclic thermal loading at multi-material gradient interfaces, published literature generally falls into two technical routes: analytical paradigm and data-driven paradigm. The analytical paradigm, centered on the finite element method (FEM), extended finite element method (XFEM), or phase-field method, characterizes cracks by arranging bonded elements or zero-thickness interface elements at the interface, or by introducing damage variables into the phase-field equations. The common practice is to first perform global meshing on the 3D CAD or micro-CT model, then write the material gradient information into the element properties layer by layer or subdomain; subsequently, couple the heat conduction equation and the elastic dynamics equation, and use an incremental-iterative approach to solve the temperature and stress fields. In post-processing, the crack initiation location is determined based on the maximum principal stress, energy release rate threshold, or J-integral criterion. When higher resolution is required at the crack tip, sub-models are enabled locally, the mesh is refined, or secondary crack descriptions are inserted to improve the numerical accuracy of the crack front; some improved schemes also combine GPU parallelism or reduced-order basis function techniques to shorten the solution time. Another empirical paradigm relies on deep learning, using a large number of simulation or experimental samples to train convolutional neural networks or graph neural networks. This transforms material distribution, boundary loads, and environmental parameters into voxel or point cloud inputs, learning the mapping between crack morphology and field quantities. After training, the network can output crack initiation probabilities and propagation paths in seconds or even milliseconds during the inference phase, making it suitable for rapid design iterations or online monitoring scenarios. Some literature also embeds physics terms into the loss function to enhance generalization ability for unseen samples.

[0003] Although analytical paradigms have clear physical mechanisms, the microstructural scale of three-dimensional multi-material gradient interfaces spans a large range. To accurately capture local high gradients, extremely fine meshes or minimal phase field length scales must be arranged at potential crack channels, leading to a sharp increase in unknowns. Simultaneously, high-temperature gradients and asymmetric constraints make the solution process highly nonlinear and have poor convergence, often requiring manual adjustment of damping or time steps to complete the computation. For cracks that jump across layers or propagate in parallel from multiple sources, interface elements and XFEMs must be dynamically inserted into cracks or the mesh reconstructed, making the process complex and prone to numerical noise. While data-driven paradigms offer fast inference speeds, they are highly sensitive to training set coverage and annotation accuracy. When material gradients, loads, or thermal cycling parameters fall within sparse sample regions, the network output often conflicts with energy conservation or mechanical equilibrium, lacking physical interpretability. More importantly, both types of techniques largely separate the simulation and prediction modules, lacking a real-time energy field-crack path bidirectional correction mechanism, making it difficult to balance efficiency and accuracy while ensuring physical consistency.

[0004] This invention aims to provide a technical solution that can achieve high-precision prediction of crack paths at multi-material gradient interfaces within limited computing resources and maintain real-time bidirectional consistency with the thermo-mechanical coupled energy field. The core technical problem to be solved can be described as: how to construct a closed-loop coupling process between the simulation domain and the graph network domain, so that the crack initiation position, propagation direction and local energy release rate are mutually calibrated, thereby improving prediction efficiency and physical reliability under complex thermal cycling conditions. Summary of the Invention

[0005] To address the aforementioned issues, this invention proposes a method for predicting microcrack paths at multi-material gradient interfaces. By using a directional complete interface diagram, a lightweight thermo-mechanical-phase field solver, and a temperature-stress-energy gated graph network coupled with a discriminator, the crack path and energy field are iteratively and synchronously corrected online. The method outputs the crack initiation location, three-dimensional crack curve, and reliability level, thereby achieving high-precision failure assessment with minimal computation.

[0006] A method for predicting microcrack paths at multi-material gradient interfaces includes the following steps:

[0007] S1. Normalized model coordinates, gradient adaptive sampling to form nodes, and edge orientation based on material similarity or geometric adjacency to generate interface graph;

[0008] S2. Lock the local domain with nodes, perform gradient meshing to reduce the order, and only calculate the thermo-mechanical-crack coupling field in the phase field of the fine mesh to extract the energy release rate, principal stress gradient and temperature.

[0009] S3, KD-Tree nearest neighbor and differentiable Gaussian interpolation map the three fields back to the graph to obtain node and edge embedding features;

[0010] S4. The temperature-stress-energy cascaded gated message passing network calculates the crack initiation probability and the crack path along the direction, and extracts the crack path according to the maximum cumulative probability.

[0011] S5. The path is meshed into an indicator tensor, and after being differiated from the energy field, it is convolved by multi-directional irradiation to obtain a consistent score, and an error gradient for probability and energy density is generated.

[0012] S6. Based on the gradient synchronization, adjust the network weights and local energy density, and repeat S2-S5 until the score converges;

[0013] S7. Freeze weights, filter fracture initiation nodes and smooth the path with splines, format the energy field and reliability label, and then output.

[0014] Optionally, step S1 may include the following sub-steps:

[0015] S11. Load the multi-material gradient interface micro-CT model and call the rigid body transformation matrix to precisely align the model's local coordinate system with the origin of the global reference coordinate system, eliminating initial measurement offsets. Perform scale normalization on the model, ensuring the distance between the two farthest nodes is normalized to 1. The resulting normalized coordinate set X = {x} k ∣k=1,…,n0}, where n0 represents the original number of grid nodes.

[0016] S12. Based on the implicit function g(r) = 0 at the material boundary and the centroid c of the particles... p Based on the distribution characteristics, adaptive sampling is implemented on the interface. The sampling density is driven by a gradient-driven function. Controlled in the steep gradient region ( To ensure complete representation of geometrical continuity details, a high ρ value is maintained in the matrix region; however, ρ is reduced in the matrix region to decrease redundant nodes. A node set is then formed after sampling. Where n is significantly less than n0, the efficiency of subsequent graph algorithms is improved.

[0017] S13, For any pair of nodes (v) i ,v j ), calculate the normalized geometric distance d ij Difference from normalized materials m ij And construct a comprehensive metric γ ij .

[0018] Only when γ ij The node pair is added to the candidate edge set ε0 only when the value is ≤1.

[0019]

[0020] Where, x i For node v i The three-dimensional coordinate vector; x j For node v j The three-dimensional coordinate vector; d max p represents the maximum Euclidean distance within the node set, used for distance normalization. i =(E i ,α i ) is node v i , where E is the material property vector. i For elastic modulus, α i p is the coefficient of thermal expansion. j =(E j ,α j ) is node v j Material property vector; m max The maximum Euclidean distance for material property differences within the node set, used for property normalization; d ij For ||x i -x j||2 / d max Normalized geometric distance; m ij For ||p i -p j ||2 / m max Normalized material differences; γ ij It is a distance-material integrated metric used for candidate edge discrimination.

[0021] S14. First, apply the material similarity criterion within the candidate set ε0. When When directly forming an edge, and according to Assign similarity weights to the materials; for the remaining node pairs, apply the geometric adjacency criterion. Perform edge patching and weighting

[0022] S15, Using nodal gradient coordinates g z (x) determines the edge direction based on the material difference. For edges that have already formed (v) i ,v j ), when (E j -E i )+κ(α i -α j When )≥0, record the direction label dir ij =+1, otherwise record dir ij =-1. The coefficient κ is given by the external working conditions, so that all edge directions are uniformly directed from the soft phase or low temperature phase to the hard phase or high temperature phase, providing a symbolic prior for the ordering of crack tip driving forces.

[0023] S16. Write the attribute vector for each node. Write an attribute vector for each edge. Node attributes and edge attributes together form a high-dimensional feature library.

[0024] S17, Viewing the diagram If an reachability search is performed and an isolated node or a non-physical short cycle below a triangular ring is found, S12 to S16 are re-executed in the corresponding local area to make the graph structure satisfy the conditions of global connectivity and local equilibrium.

[0025] S2 includes the following sub-steps:

[0026] S21. Read the coordinates of all nodes in graph S1. Constructing the envelope in the 3D model using the spread factor δ The outer region of the envelope is rigidly frozen, and only Ω′ is subsequently meshed and the field equations are solved.

[0027] S22. Using each graph node as a mandatory mesh vertex, execute a gradient-driven quad- and octagonal hybrid partitioning algorithm within Ω′ to generate a reduced-order mesh from fine to coarse.

[0028] The partition length h(x) depends on the interface gradient modulus. Automatic adjustment, the control function is:

[0029]

[0030] Where h(x) is the target element size at position x; h max h min The maximum and minimum allowable cell sizes; β is the gradient sensitivity coefficient; g z (x) represents the continuous interpolation of the gradient coordinates of nodes in S1 along the normal direction.

[0031] This control function ensures automatic subdivision to h in regions with steep gradients. min Gradually magnify to h in the flat region max This ensures that the phase field crack width can be analyzed by the smallest element, while avoiding overly dense subdivision in the far field.

[0032] S23, the ladder is based on the centroid c of the grid unit. e The gradient coordinates are used to call the material mapping function (E(c e ),α(c e Write spatially variable thermo-elastic parameters for each element; then in The external working condition for synchronous loading is the temperature cycle Θ(t) = Θ0 + ΔΘsin(2πft) and the mechanical load vector P.

[0033] S24. Only when h(x) ≤ ηh max The phase field variable φ is activated in the fine mesh region, where η is the fine mesh threshold coefficient. φ is frozen for coarse mesh elements, and only the thermo-mechanical coupling degree of freedom is retained, reducing the global degree of freedom from N0 to N0(1-ξ), where ξ is the compression ratio. This can reduce the number of degrees of freedom to above 0.7 while maintaining crack resolution accuracy.

[0034] S25. A steady-state solution is obtained by iterating in the order Θ→u→φ on a reduced-order grid. The thermal field equation, elastic equilibrium equation, and phase field evolution equation are decomposed into a sequence system:

[0035]

[0036] Among them, K Θ Here is the temperature conduction stiffness matrix; Θ (k+1) F is the temperature field vector for the (k+1)th iteration; Θ K is the thermal boundary load vector; el (Θ (k+1)) represents the elastic stiffness matrix considering thermal strain; u (k+1) For displacement field; F el M is the mechanical load vector; G is the phase field mass matrix; c φ is the interfacial fracture toughness; l is the phase field length scale; I is the identity matrix; (k+1) For phase field variables; G(u) (k+1) ) represents the crack evolution source term driven by elastic energy.

[0037] Through a thermo-mechanical-phase field serial cycle, only k max The convergence criterion ||φ can be satisfied by the round. (k+1) -φ (k) || / ||φ (k) ||<ε, where ε is 10 -4 .

[0038] S26. After solving, traverse the set of fine mesh elements. Energy release rate G is recorded at the center of gravity of each unit. e Maximum principal stress gradient Λ e With steady-state temperature Θ e .

[0039] Stack the three sets of scalars into sparse vectors G, Λ, Θ according to their cell indices, and preserve the nearest neighbor mapping index table from cells to graph nodes.

[0040] S27. Call the index table Map G, Λ, Θ and write them into the 3D slot vector [G] that has the exact same node order as the S1 graph node. node ,Λ node ,Θ node Missing mapping positions are represented by empty placeholders. Preserve the sparse format. Obtain the first-round physical prior tensor, which is seamlessly aligned with the S1 graph topology, providing input for the next step of cross-scale physical-graph feature mapping.

[0041] S3 includes the following sub-steps:

[0042] S31, Based on the set of node coordinates With S2 fine mesh element centroid set Constructing the nearest neighbor index function using KD-Tree The index radius r is chosen to be twice the phase field length scale l, so that the node can cover its crack tip influence area, while ensuring that the interpolation process is constrained within the local physical correlation range.

[0043] S32, For each node v i In the index set Perform differentiable Gaussian weighting within the inner domain to obtain the nodal field vectors s for energy release rate, maximum principal stress gradient, and temperature, respectively. i =[G i ,Λ i ,Θ i ] T :

[0044]

[0045] Among them, s i For node v i The three-dimensional physical vector; G k ,Λ k ,Θ k σ represents the energy release rate, principal stress gradient, and temperature scalar of element k; σ is the Gaussian kernel width, which is of the same order as l and is used to control the interpolation smoothness. For node v i The set of neighboring cell indices.

[0046] Using equal-weighted exponent kernels in both the numerator and denominator ensures normalized weights and differentiability, which is beneficial for subsequent end-to-end gradient backpropagation.

[0047] S33. For each directed edge e ij ∈ε, calculate the endpoint difference vector d ij =x j -x i . The node physical vector s i ,s j Along d ij By performing sign-preserving decomposition, we obtain the positive and negative components of stress and energy along the edge directions. First find the unit vector of direction. Then, project the difference between the two ends and separate the positive and negative parts.

[0048] S34, Direction maintenance code is e ij Construct a dual-channel vector with positive gain and negative dissipation Amplitude normalization is then performed to ensure that the four-dimensional components fall within [0,1]. Channel partitioning avoids positive and negative cancellation during the aggregation stage of subsequent pooling operations, thus fully preserving the directionality of the crack tip driving force.

[0049] S35. Combine the nodal three-dimensional physical vectors from step S32 with the nodal intrinsic material properties [E]. i ,α i ,g z (x i )] T By concatenating the nodes, we obtain the node embedding features. Vectors carry both intrinsic properties and local field quantities, enabling cross-scale fusion of materials and physics.

[0050] S36. Maintain the direction of step S34 while keeping the vector and the geometric properties of the edge itself. splicing to form direction-sensitive edge features

[0051] S37. According to the original index order of nodes and edges, With {f eij}e ij∈ ε write to buffer The buffer provides a one-time loading method for subsequent S4 graph networks, avoiding repeated interpolation operations during multi-threaded reading. The kernel width σ and neighborhood radius r in the formula will be automatically optimized through gradient descent during the training phase to ensure that the interpolation resolution matches the crack scale.

[0052] S4 includes the following sub-steps:

[0053] S41, Invoking Node Features Initialize trainable hidden state Call edge features Initialize trainable weights Simultaneously, establish a layer counter l=0 and set the maximum number of layers L. max This provides a starting point for subsequent recursion.

[0054] S42. In the l-th layer, along each directed edge e ij Calculate temperature difference Using the trainable amplification factor β T Construct a temperature-gated function, corresponding to the scaled principal stress gradient as follows:

[0055] S43, Adjust the steps in S42 as described above. As the gated input, with trainable coefficients β Λ Continue screening for energy release rate Prioritize temperature, stress, and energy gradients:

[0056]

[0057] in, For the l-th layer node v i v j Temperature component; σ(·) is the Sigmoid activation function; β T β Λ These are trainable gating coefficients; The principal stress gradient of the l-th layer; The energy release rate of the l-th layer edge; The principal stress gradient and energy release rate after two-stage gating.

[0058] S44, will and Combined into weighted messages, and mapped by a linear mapping matrix W e The data is passed to the terminal node, linearly superimposed with the node's self-sustaining vector, and then updated using a nonlinear function φ to update the node's hidden state.

[0059]

[0060] in, The hidden state of the l-th layer node; W v W e This is a linear mapping matrix between nodes and edges; For node v i The set of incoming edge neighborhoods; φ(·) is the Parametric-ReLU with leakage parameters; The weighting coefficient reflects the strength of the driving force after temperature-force coupling.

[0061] S45. Repeat S42 to S44 and set l←l+1. When the condition is met... or l = L max Stop the recursion at the appropriate time to ensure that the temperature-force coupling information is fully propagated within the graph and that the network state is stable.

[0062] S46, in the convergence layer l * Above, the hidden state of a node is mapped to the crack initiation probability. Map the hidden edge weights to the probability of lateral growth. This forms a node-edge joint probability graph.

[0063] S47, Select The largest node v s As the source point, within its neighborhood edge set, according to... Perform a stepwise walk until the cumulative probability falls below the threshold τ. p Or, upon reaching the preset endpoint, the first candidate crack path is generated. (1) Path Γ (1) The probability sequence is sent together with the data to the subsequent verification and encryption simulation module.

[0064] S5 includes the following sub-steps:

[0065] S51, the three-dimensional crack space curve Γ output by S4 (1) Subdivide the sample evenly according to the arc length, ensuring that the sampling interval is no greater than the minimum unit size h of S2. min This avoids mapping omissions. A nearest neighbor mapping method is used to locate each sampling point to the element index e of the S2 reduced-order mesh. Multiple records of a curve traversing the same element are retained only once to reduce redundant markings and maintain consistency at the element-level granularity. The final set of crack-penetrating elements is obtained. The collection accurately depicts the spatial projection of cracks at the current mesh resolution.

[0066] S52. Based on the meshing results, construct a crack indicator tensor χ that is consistent with the energy release rate field G at the same scale. Specifically, at ε c Assigning the value χ to the cell center within e =1, assign the value χ to the remaining units. e =0.0-1 tensors preserve crack geometry information and have a resolution that perfectly matches G, enabling direct arithmetic operations at the element-by-element level and avoiding energy ambiguity caused by traditional interpolation.

[0067] S53. Normalize the amplitude of G, then use the weighting parameter γ to balance the importance of crack geometry and energy amplitude. Next, perform element-wise multiplication in the following formula to obtain the difference tensor Δ. In multi-material gradient interface scenarios, the difference tensor can simultaneously reflect two types of physical mismatches: the crack has passed through but there is still too much energy remaining, and the energy is highly concentrated but the crack has not covered. This provides a sensitive and symmetrical input signal for the downstream convolutional block.

[0068]

[0069] Where Δ is the element-wise energy difference tensor; χ is the crack indicator tensor; G is the energy release rate field; G max γ is the maximum value in the energy release rate field, used for amplitude normalization; γ is the trade-off coefficient between crack geometry and energy amplitude; ⊙ is the Hadamard element-wise multiplication operator.

[0070] S54. Oriented core design is performed based on the radiation characteristics of anisotropic stress fields in composite materials. An internal direction set is configured. And generate an irradiation sector weight kernel K for each direction. p During the scanning process, the convolution kernel performs cross-scale sliding accumulation on the difference tensor Δ and introduces an orientation weighting factor. Sensitivity to principal stress direction modulation. This design can capture the consistency between crack propagation direction and energy gradient direction at the local scale, while simultaneously statistically analyzing residual energy inhomogeneity at the global scale, generating a physically meaningful mismatch scalar.

[0071] S55, the discriminator calculates the global mismatch output by the convolutional block. Converted into an intuitive percentage-based physical consistency score S phys The mapping formula uses a small constant c to avoid zero denominators, ensuring that low mismatches correspond to high scores and high mismatches correspond to low scores, and provides a uniform threshold τ across the score dimension. s The scoring system is used to determine the physical validity of crack paths. In practice, the scoring can quickly indicate whether the crack path satisfies the order of energy release under multi-layer material gradients, providing a quantifiable target for closed-loop optimization.

[0072] S56, when S phys Below the threshold τ s At that time, the automatic differentiation engine backpropagates the convolution computation chain of S54 to obtain the gradient of the difference tensor. Subsequently, based on the following formula, the gradient is split into two categories according to the dependency relationship: one category influences the output of the S4 network along the crack path direction. net Another type of fracture toughness G is influenced by the energy release rate field in the S2 phase field model. c :

[0073]

[0074] in, K represents the total physical mismatch. p The irradiation convolution kernel is in direction p; w p S is the direction weighting factor; phys For physics-consistent percentage fractions; c is a positive constant to prevent the denominator from being zero; o net Output the S4 graph network (joint tensor of node initial probability and edge growth probability); G c This represents the interfacial fracture toughness in the phase-field model.

[0075] S57, will The data is directly fed back to the S4 graph network to update the temperature-force coupling gating coefficients and the linear mapping matrix, making the graph network more biased towards physically consistent crack growth directions in the next iteration; at the same time, it also... Send to S2 and adjust the interface energy density parameters accordingly to match the sensitivity of the phase field model to local material gradients with the actual energy consumption of the crack.

[0076] After completing the two-way feedback, the loop re-enters S2→S3→S4→S5, until S... phys ≥τ s It may reach the upper limit of iteration. The closed-loop mechanism enables crack propagation prediction and energy field evolution to converge to a physically consistent state, significantly improving the failure assessment accuracy of multi-material gradient interfaces under thermal cycling loading.

[0077] S6 includes the following sub-steps:

[0078] S61, Gradient Unpacking

[0079] Read the gradient pairs output by the discriminator And only the component in the direction of interface energy density is retained. Subsequently, based on the cell mapping table Write the gradient tensor into the material parameter database and construct a binary mask m. tgt ∈{0,1} m And select the target unit set Ω with the threshold ζ0.tgt ={e∣m tgt (e)=1}.

[0080] S62, for Ω tgt For each unit e within the array, first query the original interface energy density. With historical offset δ c (e) Then, based on the gradient sign, perform a small adjustment in the same direction as if Then G c (e)↑, conversely G c (e)↓. The amplitude is updated via the scaling factor. Controlled, and subject to global cumulative limits ξ max Constraints, make δ is refreshed immediately after each update. c (e) Avoid non-physical crack tip overspeed propagation caused by continuous iteration.

[0081] S63, Around the target set Ω tgt Expanding one layer to form a local domain Ω δ =ring(Ω) tgt ,1). In Ω δ The phase field degree of freedom φ, displacement degree of freedom u, and temperature degree of freedom Θ are unlocked internally, while a frozen state is applied to the far field region. The incremental solver only assembles the submatrices. M δ and the coupled residual vector R δ The incremental solution Δφ is obtained through a three-step loop: residual evaluation → linearization → iterative convergence. The energy release rate tensor G(e)(e∈Ω) is updated immediately after the solution is completed. δ ), and record the local residual r loc Used for subsequent convergence statistics.

[0082] S64. Utilizing cell hash indexes Ω δ The generated new tensor G new Write back to the global tensor G. The write-back follows the rule of new values ​​overwriting old values, unupdated cells remaining unchanged, and a sparse compression is performed after completion to ensure index continuation and alignment with the graph mapping table. Consistent.

[0083] S65. Return to step S3 with the latest energy field G and recalculate the node physical vector f. v (i) and the edge physical vector f e (ij). Subsequently, the temperature-force coupled message-passing graph network is run in S4 to obtain the updated crack initiation probability. With the probability of growth in the direction S5 calls the discriminator again to calculate the physical consistency score. This completes the closed-loop mapping of G→f→p→S.

[0084] S66, Continuously track the score increase over two rounds And the iteration counter t. When |ΔS phys |<ε S or t≥N max When the system converges, the steady-state result package is exported as the final interface energy density table. Steady-state energy release rate field G * Collaborative crack path Γ ＊ and its probability sequence and convergent physical consistency score

[0085] S7 includes the following sub-steps:

[0086] S71, Receive the network weights and physical field tensors from the final round of S6, and convert the node initial probability vectors. With edge expansion probability vector Perform memory locking and generate a read-only identifier ID. snap .

[0087] S72. Set a single threshold τ according to engineering safety specifications. src ,right Threshold truncation is performed to filter out the candidate node index set Λ cand Perform hierarchical clustering based on Euclidean distance within the physical coordinate system, for Λ cand The resulting cluster set Extract centroid nodes sequentially Each cluster is merged into a single crack initiation location to avoid outputting redundant crack sources to downstream analysis modules.

[0088] S73, Set of high-probability edges along the frozen area Constructing a polyline trajectory Π based on node order k Next, spline interpolation is used in the product's global coordinate system to... k Smoothing is a quadratic continuous curve

[0089] S74. Update the energy release rate tensor G synchronously. ＊ Remapping to the original reduced-order mesh index And save it as a list of fields. <CellID,G ＊ The interface adapter is then invoked to encapsulate the list into a specific format that can be directly read by mainstream engineering CAE software.

[0090] S75, Overall Crack Initiation Probability Path cumulative probability and final state discriminant score Constructing a multidimensional feature vector q k And calculate the overall reliability level ρ within the empirical Bayesian confidence model. k ∈{I,II,III}. This level is used as a note field. Binding.

[0091] S76. Set the coordinates of the crack initiation point. Smooth crack curve Energy release rate field <CellID,G ＊ >and reliability label {ρ k} Summarized into a unified result package The result package is generated automatically. Structure and The geometry file is in dual format and includes a timestamp t. stamp and snapshot identifier ID snap To meet the requirements of version traceability and quality auditing; and to complete the delivery of the algorithm domain to the engineering application end.

[0092] This invention provides a method for predicting microcrack paths at multi-material gradient interfaces, which has beneficial effects. The technical advantages of this invention are as follows:

[0093] This application's technical solution constructs an adaptive interface graph within a geometry-material joint space, replacing the traditional approach of applying a uniform fine mesh across the entire domain. Based on a 3D solid obtained through micro-CT, the material gradient field is normalized while preserving the realistic contour, and the sampling density is driven synchronously based on the rate of change of material properties and the curvature of the topography. The generated nodes not only record geometric coordinates but also simultaneously bind initial values ​​for multiple physical quantities such as temperature, elastic modulus, and coefficient of thermal expansion. The weighted adjacency matrix between nodes is jointly determined by direction vectors, material dissimilarity, and Euclidean distance. To reduce the dimensionality of unknowns in the subsequent thermo-mechanical-phase field solution stage, higher-order shape functions are locally activated around potential crack channels, while order reduction constraints are applied to nodes far from high-gradient regions, keeping the ratio of the phase field length scale to the mesh size within a convergent range. This graph structure is logically separated from the finite element mesh, allowing nodes to be inserted or deleted as needed during iteration, avoiding large-scale reconstruction caused by dynamic crack propagation, and providing a flexible platform for the closed-loop update mechanism.

[0094] Based on the interface graph, a temperature difference-stress-energy multi-gated graph neural network is constructed to achieve edge-by-edge prediction of crack initiation and propagation probabilities. In each round of message passing, the node first calculates the thermal expansion driving force gating coefficient based on the local temperature gradient, then uses this coefficient to adjust the stress gradient information from neighboring nodes, and finally maps the gated stress gradient into the energy release rate increment and writes it into the hidden state of the target node. This two-level gating not only preserves the interpretability of explicit physical quantities during propagation but also limits the network's overfitting tendency in sparse sample regions. The direction vector introduced into the edge features can enhance the network's ability to identify nonlinear geometric effects such as crack tip deflection and bifurcation during the training phase, supporting continuous angle prediction of the three-dimensional crack front without changing the mesh topology. Unlike explicit insertion strategies that rely on critical criteria, the graph network does not require explicit geometric encoding of the crack and can output the crack initiation position and propagation channel through a probability field, providing a differentiable bridge for subsequent physical consistency correction.

[0095] To ensure real-time consistency between the prediction results and the energy field, the scheme constructs an energy-crack differential discriminator after the graph network output. The energy release rate field obtained from the phase field solution is aligned with the crack probability field output by the graph network at the voxel level. The mismatch is accumulated by sliding an irradiation convolution kernel across multiple directions and scales to generate a physical consistency score. When the score falls below a preset threshold, the framework extracts both the network weight gradient and the initial phase field gradient based on automatic differentiation. A gradient routing mechanism is used to limit high-weight updates to local subgraphs within the mismatch set. On the phase field side, the temperature-displacement-phase field triple equations are iterated again only for the corresponding subdomain, thus avoiding global repetitive solutions. The two gradient flows are synchronized in the scheduler before the next iteration, ensuring mutual calibration between the crack path and the local energy release rate until the score converges. Attached Figure Description

[0096] Figure 1 This is a flowchart of an embodiment of the present invention;

[0097] Figure 2 This is a schematic diagram of energy-crack differential discrimination according to an embodiment of the present invention;

[0098] Figure 3 This is a schematic diagram comparing the embodiments of the present invention with the prior art. Detailed Implementation

[0099] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments:

[0100] A method for predicting microcrack paths at multi-material gradient interfaces includes the following steps: (See flowchart for an example). Figure 1 As shown in the diagram, the energy-crack differential discrimination is as follows: Figure 2 As shown in the diagram, the embodiment is compared with the prior art. Figure 3 As shown;

[0101] S1. Normalized model coordinates, gradient adaptive sampling to form nodes, and edge orientation based on material similarity or geometric adjacency to generate interface graph;

[0102] S2. Lock the local domain with nodes, perform gradient meshing to reduce the order, and only calculate the thermo-mechanical-crack coupling field in the phase field of the fine mesh to extract the energy release rate, principal stress gradient and temperature.

[0103] S3, KD-Tree nearest neighbor and differentiable Gaussian interpolation map the three fields back to the graph to obtain node and edge embedding features;

[0104] S4. The temperature-stress-energy cascaded gated message passing network calculates the crack initiation probability and the crack path along the direction, and extracts the crack path according to the maximum cumulative probability.

[0105] S5. The path is meshed into an indicator tensor, and after being differiated from the energy field, it is convolved by multi-directional irradiation to obtain a consistent score, and an error gradient for probability and energy density is generated.

[0106] S6. Based on the gradient synchronization, adjust the network weights and local energy density, and repeat S2-S5 until the score converges;

[0107] S7. Freeze weights, filter fracture initiation nodes and smooth the path with splines, format the energy field and reliability label, and then output.

[0108] Optionally, step S1 may include the following sub-steps:

[0109] S11. Load the multi-material gradient interface micro-CT model and call the rigid body transformation matrix to precisely align the model's local coordinate system with the origin of the global reference coordinate system, eliminating initial measurement offsets. Perform scale normalization on the model, ensuring the distance between the two farthest nodes is normalized to 1. The resulting normalized coordinate set X = {x} k ∣k=1,…,n0}, where n0 represents the original number of grid nodes.

[0110] S12. Based on the implicit function g(r) = 0 at the material boundary and the centroid c of the particles... p Based on the distribution characteristics, adaptive sampling is implemented on the interface. The sampling density is driven by a gradient-driven function. Controlled in the steep gradient region ( To ensure complete representation of geometrical continuity details, a high ρ value is maintained in the matrix region; however, ρ is reduced in the matrix region to decrease redundant nodes. A node set is then formed after sampling. Where n is significantly less than n0, the efficiency of subsequent graph algorithms is improved.

[0111] S13, For any pair of nodes (v) i ,v j), calculate the normalized geometric distance d ij Difference from normalized materials m ij And construct a comprehensive metric γ ij .

[0112] Only when γ ij The node pair is added to the candidate edge set ε0 only when the value is ≤1.

[0113]

[0114] Where, x i For node v i The three-dimensional coordinate vector; x j For node v j The three-dimensional coordinate vector; d max p represents the maximum Euclidean distance within the node set, used for distance normalization. i =(E i ,α i ) is node v i , where E is the material property vector. i For elastic modulus, α i p is the coefficient of thermal expansion. j =(E j ,α j ) is node v j Material property vector; m max The maximum Euclidean distance for material property differences within the node set, used for property normalization; d ij For ||x i -x j ||2 / d max Normalized geometric distance; m ij For ||p i -p j ||2 / m max Normalized material differences; γ ij It is a distance-material integrated metric used for candidate edge discrimination.

[0115] S14. First, apply the material similarity criterion within the candidate set ε0. When When directly forming an edge, and according to Assign similarity weights to the materials; for the remaining node pairs, apply the geometric adjacency criterion. Perform edge patching and weighting

[0116] S15, Using nodal gradient coordinates g z (x) determines the edge direction based on the material difference. For edges that have already formed (v) i ,v j ), when (E j -E i )+κ(αi -α j When )≥0, record the direction label dir ij =+1, otherwise record dir ij =-1. The coefficient κ is given by the external working conditions, so that all edge directions are uniformly directed from the soft phase or low temperature phase to the hard phase or high temperature phase, providing a symbolic prior for the ordering of crack tip driving forces.

[0117] S16. Write the attribute vector for each node. Write an attribute vector for each edge. Node attributes and edge attributes together form a high-dimensional feature library.

[0118] S17, Viewing the diagram If an reachability search is performed and an isolated node or a non-physical short cycle below a triangular ring is found, S12 to S16 are re-executed in the corresponding local area to make the graph structure satisfy the conditions of global connectivity and local equilibrium.

[0119] S2 includes the following sub-steps:

[0120] S21. Read the coordinates of all nodes in graph S1. Constructing the envelope in the 3D model using the spread factor δ The outer region of the envelope is rigidly frozen, and only Ω′ is subsequently meshed and the field equations are solved.

[0121] S22. Using each graph node as a mandatory mesh vertex, execute a gradient-driven quad- and octagonal hybrid partitioning algorithm within Ω′ to generate a reduced-order mesh from fine to coarse.

[0122] The partition length h(x) depends on the interface gradient modulus. Automatic adjustment, the control function is:

[0123]

[0124] Where h(x) is the target element size at position x; h max h min The maximum and minimum allowable cell sizes; β is the gradient sensitivity coefficient; g z (x) represents the continuous interpolation of the gradient coordinates of nodes in S1 along the normal direction.

[0125] This control function ensures automatic subdivision to h in regions with steep gradients. min Gradually magnify to h in the flat region max This ensures that the phase field crack width can be analyzed by the smallest element, while avoiding overly dense subdivision in the far field.

[0126] S23, the ladder is based on the centroid c of the grid unit. eThe gradient coordinates are used to call the material mapping function (E(c e ),α(c e Write spatially variable thermo-elastic parameters for each element; then in The external working condition for synchronous loading is the temperature cycle Θ(t) = Θ0 + ΔΘsin(2πft) and the mechanical load vector P.

[0127] S24. Only when h(x) ≤ ηh max The phase field variable φ is activated in the fine mesh region, where η is the fine mesh threshold coefficient. φ is frozen for coarse mesh elements, and only the thermo-mechanical coupling degree of freedom is retained, reducing the global degree of freedom from N0 to N0(1-ξ), where ξ is the compression ratio. This can reduce the number of degrees of freedom to above 0.7 while maintaining crack resolution accuracy.

[0128] S25. A steady-state solution is obtained by iterating in the order Θ→u→φ on a reduced-order grid. The thermal field equation, elastic equilibrium equation, and phase field evolution equation are decomposed into a sequence system:

[0129]

[0130] Among them, K Θ Here is the temperature conduction stiffness matrix; Θ (k+1) F is the temperature field vector for the (k+1)th iteration; Θ K is the thermal boundary load vector; el (Θ (k+1) ) represents the elastic stiffness matrix considering thermal strain; u (k+1) For displacement field; F el M is the mechanical load vector; G is the phase field mass matrix; c φ is the interfacial fracture toughness; l is the phase field length scale; I is the identity matrix; (k+1) For phase field variables; G(u) (k+1) ) represents the crack evolution source term driven by elastic energy.

[0131] Through a thermo-mechanical-phase field serial cycle, only k max The convergence criterion ||φ can be satisfied by the round. (k+1) -φ (k) || / ||φ (k) ||<ε, where ε is 10 -4 .

[0132] S26. After solving, traverse the set of fine mesh elements. Energy release rate G is recorded at the center of gravity of each unit. e Maximum principal stress gradient Λ e With steady-state temperature Θ e .

[0133] Stack the three sets of scalars into sparse vectors G, Λ, Θ according to their cell indices, and preserve the nearest neighbor mapping index table from cells to graph nodes.

[0134] S27. Call the index table Map G, Λ, Θ and write them into the 3D slot vector [G] that has the exact same node order as the S1 graph node. node ,Λ node ,Θ node Missing mapping positions are represented by empty placeholders. Preserve the sparse format. Obtain the first-round physical prior tensor, which is seamlessly aligned with the S1 graph topology, providing input for the next step of cross-scale physical-graph feature mapping.

[0135] S3 includes the following sub-steps:

[0136] S31, Based on the set of node coordinates With S2 fine mesh element centroid set Constructing the nearest neighbor index function using KD-Tree The index radius r is chosen to be twice the phase field length scale l, so that the node can cover its crack tip influence area, while ensuring that the interpolation process is constrained within the local physical correlation range.

[0137] S32, For each node v i In the index set Perform differentiable Gaussian weighting within the inner domain to obtain the nodal field vectors s for energy release rate, maximum principal stress gradient, and temperature, respectively. i =[G i ,Λ i ,Θ i ] T :

[0138]

[0139] Among them, s i For node v i The three-dimensional physical vector; G k ,Λ k ,Θ k σ represents the energy release rate, principal stress gradient, and temperature scalar of element k; σ is the Gaussian kernel width, which is of the same order as l and is used to control the interpolation smoothness. For node v i The set of neighboring cell indices.

[0140] Using equal-weighted exponent kernels in both the numerator and denominator ensures normalized weights and differentiability, which is beneficial for subsequent end-to-end gradient backpropagation.

[0141] S33. For each directed edge e ij ∈ε, calculate the endpoint difference vector dij =x j -x i . The node physical vector s i ,s j Along d ij By performing sign-preserving decomposition, we obtain the positive and negative components of stress and energy along the edge directions. First find the unit vector of direction. Then, project the difference between the two ends and separate the positive and negative parts.

[0142] S34, Direction maintenance code is e ij Construct a dual-channel vector with positive gain and negative dissipation Amplitude normalization is then performed to ensure that the four-dimensional components fall within [0,1]. Channel partitioning avoids positive and negative cancellation during the aggregation stage of subsequent pooling operations, thus fully preserving the directionality of the crack tip driving force.

[0143] S35. Combine the nodal three-dimensional physical vectors from step S32 with the nodal intrinsic material properties [E]. i ,α i ,g z (x i )] T By concatenating the nodes, we obtain the node embedding features. Vectors carry both intrinsic properties and local field quantities, enabling cross-scale fusion of materials and physics.

[0144] S36. Maintain the direction of step S34 while keeping the vector and the geometric properties of the edge itself. splicing to form direction-sensitive edge features

[0145] S37. According to the original index order of nodes and edges, and Write to buffer The buffer provides a one-time loading method for subsequent S4 graph networks, avoiding repeated interpolation operations during multi-threaded reading. The kernel width σ and neighborhood radius r in the formula will be automatically optimized through gradient descent during the training phase to ensure that the interpolation resolution matches the crack scale.

[0146] S4 includes the following sub-steps:

[0147] S41, Invoking Node Features Initialize trainable hidden state Call edge features Initialize trainable weights Simultaneously, establish a layer counter l=0 and set the maximum number of layers L. max This provides a starting point for subsequent recursion.

[0148] S42. In the l-th layer, along each directed edge e ij Calculate temperature difference Using the trainable amplification factor β T Construct a temperature-gated function, corresponding to the scaled principal stress gradient as follows:

[0149] S43, Adjust the steps in S42 as described above. As the gated input, with trainable coefficients β Λ Continue screening for energy release rate Prioritize temperature, stress, and energy gradients:

[0150]

[0151] in, For the l-th layer node v i v j Temperature component; σ(·) is the Sigmoid activation function; β T β Λ These are trainable gating coefficients; The principal stress gradient of the l-th layer; The energy release rate of the l-th layer edge; The principal stress gradient and energy release rate after two-stage gating.

[0152] S44, will and Combined into weighted messages, and mapped by a linear mapping matrix W e The data is passed to the terminal node, linearly superimposed with the node's self-sustaining vector, and then updated using a nonlinear function φ to update the node's hidden state.

[0153]

[0154] in, The hidden state of the l-th layer node; W v W e This is a linear mapping matrix between nodes and edges; For node v i The set of incoming edge neighborhoods; φ(·) is the Parametric-ReLU with leakage parameters; The weighting coefficient reflects the strength of the driving force after temperature-force coupling.

[0155] S45. Repeat S42 to S44 and set l←l+1. When the condition is met... or l = L max Stop the recursion at the appropriate time to ensure that the temperature-force coupling information is fully propagated within the graph and that the network state is stable.

[0156] S46, in the convergence layer l *Above, the hidden state of a node is mapped to the crack initiation probability. Map the hidden edge weights to the probability of lateral growth. This forms a node-edge joint probability graph.

[0157] S47, Select The largest node v s As the source point, within its neighborhood edge set, according to... Perform a stepwise walk until the cumulative probability falls below the threshold τ. p Or, upon reaching the preset endpoint, the first candidate crack path is generated. (1) Path Γ (1) The probability sequence is sent together with the data to the subsequent verification and encryption simulation module.

[0158] S5 includes the following sub-steps:

[0159] S51, the three-dimensional crack space curve Γ output by S4 (1) Subdivide the sample evenly according to the arc length, ensuring that the sampling interval is no greater than the minimum unit size h of S2. min This avoids mapping omissions. A nearest neighbor mapping method is used to locate each sampling point to the element index e of the S2 reduced-order mesh. Multiple records of a curve traversing the same element are retained only once to reduce redundant markings and maintain consistency at the element-level granularity. The final set of crack-penetrating elements is obtained. The collection accurately depicts the spatial projection of cracks at the current mesh resolution.

[0160] S52. Based on the meshing results, construct a crack indicator tensor χ that is consistent with the energy release rate field G at the same scale. Specifically, at ε c Assigning the value χ to the cell center within e =1, assign the value χ to the remaining units. e =0.0-1 tensors preserve crack geometry information and have a resolution that perfectly matches G, enabling direct arithmetic operations at the element-by-element level and avoiding energy ambiguity caused by traditional interpolation.

[0161] S53. Normalize the amplitude of G, then use the weighting parameter γ to balance the importance of crack geometry and energy amplitude. Next, perform element-wise multiplication in the following formula to obtain the difference tensor Δ. In multi-material gradient interface scenarios, the difference tensor can simultaneously reflect two types of physical mismatches: the crack has passed through but there is still too much energy remaining, and the energy is highly concentrated but the crack has not covered. This provides a sensitive and symmetrical input signal for the downstream convolutional block.

[0162]

[0163] Where Δ is the element-wise energy difference tensor; χ is the crack indicator tensor; G is the energy release rate field; Gmax γ is the maximum value in the energy release rate field, used for amplitude normalization; γ is the trade-off coefficient between crack geometry and energy amplitude; ⊙ is the Hadamard element-wise multiplication operator.

[0164] S54. Oriented core design is performed based on the radiation characteristics of anisotropic stress fields in composite materials. An internal direction set is configured. And generate an irradiation sector weight kernel K for each direction. p During the scanning process, the convolution kernel performs cross-scale sliding accumulation on the difference tensor Δ and introduces an orientation weighting factor. Sensitivity to principal stress direction modulation. This design can capture the consistency between crack propagation direction and energy gradient direction at the local scale, while simultaneously statistically analyzing residual energy inhomogeneity at the global scale, generating a physically meaningful mismatch scalar.

[0165] S55, the discriminator calculates the global mismatch output by the convolutional block. Converted into an intuitive percentage-based physical consistency score S phys The mapping formula uses a small constant c to avoid zero denominators, ensuring that low mismatches correspond to high scores and high mismatches correspond to low scores, and provides a uniform threshold τ across the score dimension. s The scoring system is used to determine the physical validity of crack paths. In practice, the scoring can quickly indicate whether the crack path satisfies the order of energy release under multi-layer material gradients, providing a quantifiable target for closed-loop optimization.

[0166] S56, when S phys Below the threshold τ s At that time, the automatic differentiation engine backpropagates the convolution computation chain of S54 to obtain the gradient of the difference tensor. Subsequently, based on the following formula, the gradient is split into two categories according to the dependency relationship: one category influences the output of the S4 network along the crack path direction. net Another type of fracture toughness G is influenced by the energy release rate field in the S2 phase field model. c :

[0167]

[0168] in, K represents the total physical mismatch. p The irradiation convolution kernel is in direction p; w p S is the direction weighting factor; phys For physics-consistent percentage fractions; c is a positive constant to prevent the denominator from being zero; o net Output the S4 graph network (joint tensor of node initial probability and edge growth probability); G c This represents the interfacial fracture toughness in the phase-field model.

[0169] S57, will The data is directly fed back to the S4 graph network to update the temperature-force coupling gating coefficients and the linear mapping matrix, making the graph network more biased towards physically consistent crack growth directions in the next iteration; at the same time, it also... Send to S2 and adjust the interface energy density parameters accordingly to match the sensitivity of the phase field model to local material gradients with the actual energy consumption of the crack.

[0170] After completing the two-way feedback, the loop re-enters S2→S3→S4→S5, until S... phys ≥τ s It may reach the upper limit of iteration. The closed-loop mechanism enables crack propagation prediction and energy field evolution to converge to a physically consistent state, significantly improving the failure assessment accuracy of multi-material gradient interfaces under thermal cycling loading.

[0171] S6 includes the following sub-steps:

[0172] S61, Gradient Unpacking

[0173] Read the gradient pairs output by the discriminator And only the component in the direction of interface energy density is retained. Subsequently, based on the cell mapping table Write the gradient tensor into the material parameter database and construct a binary mask m. tgt ∈{0,1} m And select the target unit set Ω with the threshold ζ0. tgt ={e∣m tgt (e)=1}.

[0174] S62, for Ω tgt For each unit e within the array, first query the original interface energy density. With historical offset δ c (e) Then, based on the gradient sign, perform a small adjustment in the same direction as if Then G c (e)↑, conversely G c (e)↓. The amplitude is updated via the scaling factor η. Gc Controlled, and subject to global cumulative limits ξ max Constraints, make δ is refreshed immediately after each update. c (e) Avoid non-physical crack tip overspeed propagation caused by continuous iteration.

[0175] S63, Around the target set Ω tgt Expanding one layer to form a local domain Ω δ =ring(Ω) tgt ,1). In Ω δ The phase field degree of freedom φ, displacement degree of freedom u, and temperature degree of freedom Θ are unlocked internally, while a frozen state is applied to the far field region. The incremental solver only assembles the submatrices. M δ and the coupled residual vector R δ The incremental solution Δφ is obtained through a three-step loop: residual evaluation → linearization → iterative convergence. The energy release rate tensor G(e)(e∈Ω) is updated immediately after the solution is completed. δ ), and record the local residual r loc Used for subsequent convergence statistics.

[0176] S64. Utilizing cell hash indexes Ω δ The generated new tensor G new Write back to the global tensor G. The write-back follows the rule of new values ​​overwriting old values, unupdated cells remaining unchanged, and a sparse compression is performed after completion to ensure index continuation and alignment with the graph mapping table. Consistent.

[0177] S65. Return to step S3 with the latest energy field G and recalculate the node physical vector f. v (i) and the edge physical vector f e (ij). Subsequently, the temperature-force coupled message-passing graph network is run in S4 to obtain the updated crack initiation probability. With the probability of growth in the direction S5 calls the discriminator again to calculate the physical consistency score. This completes the closed-loop mapping of G→f→p→S.

[0178] S66, Continuously track the score increase over two rounds And the iteration counter t. When |ΔS phys |<ε S or t≥N max When the system converges, the steady-state result package is exported as the final interface energy density table. Steady-state energy release rate field G ＊ Collaborative crack path Γ ＊ and its probability sequence and convergent physical consistency score

[0179] S7 includes the following sub-steps:

[0180] S71, Receive the network weights and physical field tensors from the final round of S6, and convert the node initial probability vectors. With edge expansion probability vector Perform memory locking and generate a read-only identifier ID. snap .

[0181] S72. Set a single threshold τ according to engineering safety specifications. src ,right Threshold truncation is performed to filter out the candidate node index set Λ cand Perform hierarchical clustering based on Euclidean distance within the physical coordinate system, for Λ cand The resulting cluster set Extract centroid nodes sequentially Each cluster is merged into a single crack initiation location to avoid outputting redundant crack sources to downstream analysis modules.

[0182] S73, Set of high-probability edges along the frozen area Constructing a polyline trajectory Π based on node order k Next, spline interpolation is used in the product's global coordinate system to... k Smoothing is a quadratic continuous curve

[0183] S74. Update the energy release rate tensor G synchronously. ＊ Remapping to the original reduced-order mesh index And save it as a list of fields. <CellID,G ＊ The interface adapter is then invoked to encapsulate the list into a specific format that can be directly read by mainstream engineering CAE software.

[0184] S75, Overall Crack Initiation Probability Path cumulative probability and final state discriminant score Constructing a multidimensional feature vector q k And calculate the overall reliability level ρ within the empirical Bayesian confidence model. k ∈{I,II,III}. This level is used as a note field. Binding.

[0185] S76. Set the coordinates of the crack initiation point. Smooth crack curve Energy release rate field <CellID,G ＊ >and reliability label {ρ k} Summarized into a unified result package The result package automatically generates a JSON structure and a STEP geometry file in dual format, and includes a timestamp t. stamp and snapshot identifier ID snap To meet the requirements of version traceability and quality auditing; and to complete the delivery of the algorithm domain to the engineering application end.

[0186] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any other way. Any modifications or equivalent changes made based on the technical essence of the present invention shall still fall within the scope of protection claimed by the present invention.

Claims

1. A method for predicting microcrack paths at multi-material gradient interfaces, comprising the following specific steps, characterized in that: S1. Normalized model coordinates, gradient adaptive sampling to form nodes, and edge orientation based on material similarity or geometric adjacency to generate interface graph; S2. Lock the local domain with nodes, perform gradient meshing to reduce the order, and only calculate the thermo-mechanical-crack coupling field in the phase field of the fine mesh to extract the energy release rate, principal stress gradient and temperature. S3, KD-Tree nearest neighbor and differentiable Gaussian interpolation map the three fields back to the graph to obtain node and edge embedding features; S4. The temperature-stress-energy cascaded gated message passing network calculates the crack initiation probability and the crack path along the direction, and extracts the crack path according to the maximum cumulative probability. S5. The path is meshed into an indicator tensor, and after being differiated from the energy field, it is convolved by multi-directional irradiation to obtain a consistent score, and an error gradient for probability and energy density is generated. S6. Based on the gradient synchronization, adjust the network weights and local energy density, and repeat S2-S5 until the score converges; S7. Freeze weights, filter fracture initiation nodes and smooth the path with splines, format the energy field and reliability label, and then output.

2. The method for predicting microcrack paths at multi-material gradient interfaces according to claim 1, characterized in that: Step S1 includes the following sub-steps: S11. Load the multi-material gradient interface micro-CT model and call the rigid body transformation matrix to accurately align the model's local coordinate system with the origin of the global reference coordinate system, eliminating initial measurement offsets. Perform scale normalization on the model to ensure the distance between the two farthest nodes is normalized to 1. After processing, the normalized coordinate set X = {x k |k=1,…,n0}, where n0 represents the original number of mesh nodes; S12. Based on the implicit function g(r) = 0 at the material boundary and the centroid c of the particles... p Based on the distribution characteristics, adaptive sampling is implemented on the interface, and the sampling density is driven by the gradient-driven function. Controlled in the steep gradient region ( Maintaining a high ρ value ensures that geometrically continuous details are fully expressed; reducing ρ in the matrix region reduces redundant nodes, and a node set is formed after sampling. Where n is significantly less than n0, the efficiency of subsequent graph algorithms is improved; S13, For any pair of nodes (v) i ,v j ), calculate the normalized geometric distance d ij Difference from normalized materials m ij And construct a comprehensive metric γ ij ; Only when γ ij The node pair is added to the candidate edge set ε0 only when the value is ≤1. Where, x i For node v i The three-dimensional coordinate vector; x j For node v j The three-dimensional coordinate vector; d max p represents the maximum Euclidean distance within the node set, used for distance normalization. i =(E i ,α i ) is node v i , where E is the material property vector. i For elastic modulus, α i p is the coefficient of thermal expansion. j =(E j ,α j ) is node v j Material property vector; m max The maximum Euclidean distance for material property differences within the node set, used for property normalization; d ij For ||x i -x j ||2 / d max Normalized geometric distance; m ij For ||p i -p j ||2 / m max Normalized material differences; γ ij It is a distance-material integrated metric used for candidate edge discrimination; S14. First, apply the material similarity criterion within the candidate set ε0, when... When directly forming an edge, and according to Assign similarity weights to the materials; for the remaining node pairs, apply the geometric adjacency criterion. Perform edge patching and weighting S15, Using nodal gradient coordinates g z (x) and material difference determine the edge direction, for edges that have already formed (v) i ,v j ), when (E j -E i )+κ(α i -α j When )≥0, record the direction label dir ij =+1, otherwise record dir ij =-1, the coefficient κ is given by the external working conditions, so that all edge directions are uniform from the soft phase or low temperature phase to the hard phase or high temperature phase, providing a symbolic prior for the ordering of crack tip driving forces; S16. Write the attribute vector for each node. Write an attribute vector for each edge. Node attributes and edge attributes together form a high-dimensional feature library. S17, Viewing the diagram If an reachability search is performed and an isolated node or a non-physical short cycle below a triangular ring is found, S12 to S16 are re-executed in the corresponding local area to make the graph structure satisfy the conditions of global connectivity and local equilibrium.

3. The method for predicting microcrack paths at multi-material gradient interfaces according to claim 1, characterized in that: Step S2 includes the following sub-steps: S21. Read the coordinates of all nodes in graph S1. Constructing the envelope in the 3D model using the spread factor δ Rigid freezing is applied to the outer region of the envelope, and subsequent mesh implantation and field equation solving are performed only on Ω′. S22. Using each graph node as a mandatory mesh vertex, execute a gradient-driven quad- and octagonal hybrid partitioning algorithm within Ω′ to generate a reduced-order mesh from fine to coarse. The partition length h(x) depends on the interface gradient modulus. Automatic adjustment, the control function is: Where h(x) is the target element size at position x; h max h min The maximum and minimum allowable cell sizes; β is the gradient sensitivity coefficient; g z (x) represents the continuous interpolation of the gradient coordinates of nodes in S1 along the normal direction. This control function ensures automatic subdivision to h in regions with steep gradients. min Gradually magnify to h in the flat region max This ensures that the phase field crack width can be analyzed by the smallest element, while avoiding overly dense subdivision in the far field; S23, the ladder is based on the centroid c of the grid unit. e The gradient coordinates are used to call the material mapping function (E(c e ),α(c e Write spatially variable thermo-elastic parameters for each element; then in The external working condition for synchronous loading is the temperature cycle Θ(t) = Θ0 + ΔΘsin(2πft) and the mechanical load vector P; S24. Only when h(x) ≤ ηh max The phase field variable φ is activated in the fine mesh region, where η is the fine mesh threshold coefficient. φ is frozen for coarse mesh elements and only the thermo-mechanical coupling degree of freedom is retained, reducing the number of global degrees of freedom from N0 to N0(1-ξ), where ξ is the compression ratio, which can reduce it to more than 0.7 while maintaining crack resolution accuracy. S25. A steady-state solution is obtained by iterating in the order Θ→u→φ on a reduced-order grid. The thermal field equation, elastic equilibrium equation, and phase field evolution equation are decomposed into a sequence system: Among them, K Θ Here is the temperature conduction stiffness matrix; Θ (k+1) F is the temperature field vector for the (k+1)th iteration; Θ K is the thermal boundary load vector; el (Θ (k+1) ) represents the elastic stiffness matrix considering thermal strain; u (k+1) For displacement field; F el M is the mechanical load vector; G is the phase field mass matrix; c φ is the interfacial fracture toughness; l is the phase field length scale; I is the identity matrix; (k+1) For phase field variables; G(u) (k +1) ) represents the crack evolution source term driven by elastic energy; Through a thermo-mechanical-phase field serial cycle, only k max The convergence criterion ||φ can be satisfied by the round. (k+1) -φ (k) || / ||φ (k) ||<ε, where ε is 10 -4 ; S26. After solving, traverse the set of fine mesh elements. Energy release rate G is recorded at the center of gravity of each unit. e Maximum principal stress gradient Λ e With steady-state temperature Θ e ; Stack the three sets of scalars into sparse vectors G, Λ, Θ according to their cell indices, and preserve the nearest neighbor mapping index table from cells to graph nodes. S27. Call the index table Map G, Λ, Θ and write them into the 3D slot vector [G] that has the exact same node order as the S1 graph node. node ,Λ node ,Θ node The missing mapping position uses an empty placeholder. By maintaining the sparse format, we obtain the first-round physical prior tensor, which is seamlessly aligned with the topology of the S1 graph, providing input for the next step of cross-scale physical-graph feature mapping.

4. The method for predicting microcrack paths at multi-material gradient interfaces according to claim 1, characterized in that: Step S3 includes the following sub-steps: S31, Based on the set of node coordinates With S2 fine mesh element centroid set Constructing the nearest neighbor index function using KD-Tree The index radius r is selected as twice the phase field length scale l, so that the node can cover its crack tip influence area, while ensuring that the constraints during the interpolation process are within the local physical correlation range; S32, For each node v i In the index set Perform differentiable Gaussian weighted sampling to obtain the nodal field vectors for energy release rate, maximum principal stress gradient, and temperature, respectively. Among them, s i For node v i The three-dimensional physical vector; G k ,Λ k ,Θ k σ represents the energy release rate, principal stress gradient, and temperature scalar of element k; σ is the Gaussian kernel width, which is of the same order as l and is used to control the interpolation smoothness. For node v i The set of neighboring cell indices; Using equal-weighted exponent kernels in both the numerator and denominator ensures normalized weights and differentiability, which is beneficial for subsequent end-to-end gradient backpropagation. S33. For each directed edge e ij ∈ε, calculate the endpoint difference vector d ij =x j -x i . The node physical vector s i ,s j Along d ij By performing sign-preserving decomposition, we obtain the positive and negative components of stress and energy along the edge directions. First find the unit vector of direction. Then project the difference between the two ends and separate the positive and negative parts; S34, Direction maintenance code is e ij Construct a dual-channel vector with positive gain and negative dissipation Then, amplitude normalization is performed to ensure that the four-dimensional components fall into [0,1]. Channel partitioning avoids positive and negative cancellation during the aggregation stage of subsequent pooling operations, thus fully preserving the directionality of the crack tip driving force. S35. Combine the nodal three-dimensional physical vectors from step S32 with the nodal intrinsic material properties [E]. i ,α i ,g z (x i )] T By concatenating the nodes, we obtain the node embedding features. Vectors carry both intrinsic properties and local field quantities, enabling cross-scale fusion of materials and physics; S36. Maintain the direction of step S34 while keeping the vector and the geometric properties of the edge itself. splicing to form direction-sensitive edge features S37. According to the original index order of nodes and edges, and Write to buffer The buffer provides a one-time loading method for the subsequent S4 graph network, avoiding repeated interpolation operations during multi-threaded reading. The kernel width σ and neighborhood radius r in the formula will be automatically optimized through gradient descent during the training phase to ensure that the interpolation resolution matches the crack scale.

5. The method for predicting microcrack paths at multi-material gradient interfaces according to claim 1, characterized in that: Step S4 includes the following sub-steps: S41, Invoking Node Features Initialize trainable hidden state Call edge features Initialize trainable weights Simultaneously, establish a layer counter l=0 and set the maximum number of layers L. max This provides a starting point for subsequent recursion; S42. In the l-th layer, along each directed edge e ij Calculate temperature difference Using the trainable amplification factor β T Construct a temperature-gated function, corresponding to the scaled principal stress gradient as follows: S43, Adjust the steps in S42 as described above. As the gated input, with trainable coefficients β Λ Continue screening for energy release rate Prioritize temperature, stress, and energy gradients: in, For the l-th layer node v i v j Temperature component; σ(·) is the Sigmoid activation function; β T β Λ These are trainable gating coefficients; The principal stress gradient of the l-th layer; The energy release rate of the l-th layer edge; The principal stress gradient and energy release rate after two-stage gating; S44, will and Combined into weighted messages, and mapped by a linear mapping matrix W e The data is passed to the terminal node, linearly superimposed with the node's self-sustaining vector, and then updated using a nonlinear function φ to update the node's hidden state. in, The hidden state of the l-th layer node; W v W e This is a linear mapping matrix between nodes and edges; For node v i The set of incoming edge neighborhoods; φ(·) is the Parametric-ReLU with leakage parameters; The weighting coefficient reflects the strength of the driving force after temperature-force coupling; S45. Repeat S42 to S44 and let l ← l + 1, until the condition is met. or l = L max Stop the recursion at that time to ensure that the temperature-force coupling information is fully propagated within the graph and the network state is stable; S46, in the convergence layer l * Above, the hidden state of a node is mapped to the crack initiation probability. Map the hidden edge weights to the probability of lateral growth. This forms a node-edge joint probability graph. S47, Select The largest node v s As the source point, within its neighborhood edge set, according to... Perform a stepwise walk until the cumulative probability falls below the threshold τ. p Or, upon reaching the preset endpoint, the first candidate crack path is generated. (1) , path Γ (1) The probability sequence is sent together with the data to the subsequent verification and encryption simulation module.

6. The method for predicting microcrack paths at multi-material gradient interfaces according to claim 1, characterized in that: Step S5 includes the following sub-steps: S51, the three-dimensional crack space curve Γ output by S4 (1) Subdivide the sample evenly according to the arc length, ensuring that the sampling interval is no greater than the minimum unit size h of S2. min To avoid mapping omissions, a nearest neighbor mapping method is used to locate each sampling point to the element index e of the S2 reduced-order mesh. For multiple records of a curve crossing the same element, only one record is retained to reduce redundant markings and maintain consistency at the element level, ultimately resulting in a set of crack-penetrating elements. The collection accurately depicts the spatial projection of cracks at the current mesh resolution; S52. Based on the meshing results, construct a crack indicator tensor χ that is consistent with the energy release rate field G at the same scale. Specifically, at ε c Assigning the value χ to the cell center within e =1, assign the value χ to the remaining units. e =0,0-1 tensors preserve crack geometry information and have a resolution that perfectly matches G, enabling direct arithmetic operations at the element-by-element level and avoiding energy ambiguity caused by traditional interpolation; S53. Normalize the amplitude of G, then use the weighted parameter γ to balance the importance of crack geometry and energy amplitude. Subsequently, perform element-wise multiplication in the following formula to obtain the difference tensor Δ. In the scenario of multi-material gradient interfaces, the difference tensor can simultaneously reflect two types of physical mismatches: the crack has passed through but there is too much energy left, and the energy is highly concentrated but the crack has not covered. This provides a sensitive and symmetrical input signal for the downstream convolutional block. Where Δ is the element-wise energy difference tensor; χ is the crack indicator tensor; G is the energy release rate field; G max is the maximum value in the energy release rate field, used for amplitude normalization; γ is the trade-off coefficient between crack geometry and energy amplitude; ⊙ is the Hadamard element-wise multiplication operator; S54. A directional core design is performed to address the radiation characteristics of anisotropic stress fields in composite materials, with an internal direction set. And generate an irradiation sector weight kernel K for each direction. p During the scanning process, the convolution kernel performs cross-scale sliding accumulation on the difference tensor Δ and introduces an orientation weighting factor. By modulating the principal stress direction sensitivity, the consistency between the crack propagation direction and the energy gradient direction can be captured at the local scale, while the residual energy inhomogeneity can be statistically analyzed at the global scale to generate a mismatch scalar with clear physical meaning. S55, the discriminator calculates the global mismatch output by the convolutional block. Converted into an intuitive percentage-based physical consistency score S phys The mapping formula uses a small constant c to avoid zero denominators, ensuring that low mismatch corresponds to high scores and high mismatch corresponds to low scores, and provides a uniform threshold τ in the score dimension. s To determine the physical rationality of crack paths, in terms of scenarios, the scoring can quickly indicate whether the crack paths satisfy the order of energy release under multi-layer material gradients, providing quantifiable targets for closed-loop optimization; S56, when S phys Below the threshold τ s At that time, the automatic differentiation engine backpropagates the convolution computation chain of S54 to obtain the gradient of the difference tensor. Subsequently, based on the following formula, the gradient is split into two categories according to the dependency relationship: one category influences the output of the S4 network along the crack path direction. net Another type of fracture toughness G is influenced by the energy release rate field in the S2 phase field model. c : in, K represents the total physical mismatch. p The irradiation convolution kernel is in direction p; w p S is the direction weighting factor; phys For physics-consistent percentage fractions; c is a positive constant to prevent the denominator from being zero; o net For the output of the S4 graph network, the joint tensor of node initial probability and edge growth probability; G c For phase-field model interface fracture toughness; S57, will The data is directly fed back to the S4 graph network to update the temperature-force coupling gating coefficients and the linear mapping matrix, making the graph network more biased towards physically consistent crack growth directions in the next iteration; at the same time, it also... Send to S2 and adjust the interface energy density parameters accordingly to match the sensitivity of the phase field model to local material gradients with the actual crack energy consumption. After completing the two-way feedback, the loop re-enters S2→S3→S4→S5, until S... phys ≥τ s Or, upon reaching the upper limit of iteration, the closed-loop mechanism enables crack propagation prediction and energy field evolution to converge to a physically consistent state, significantly improving the failure assessment accuracy of multi-material gradient interfaces under thermal cycling loading.

7. The method for predicting microcrack paths at multi-material gradient interfaces according to claim 1, characterized in that: Step S6 includes the following sub-steps: S61, Gradient unpacking; Read the gradient pairs output by the discriminator And only the component in the direction of interface energy density is retained. Subsequently, based on the cell mapping table Write the gradient tensor into the material parameter database and construct a binary mask m. tgt ∈{0,1} m And select the target unit set Ω with the threshold ζ0. tgt ={e∣m tgt (e)=1}; S62, for Ω tgt For each unit e within the array, first query the original interface energy density. With historical offset δ c (e) Then, based on the gradient sign, perform a small adjustment in the same direction as if Then G c (e)↑, conversely G c (e)↓, update the amplitude using the scaling factor Controlled, and subject to global cumulative limits ξ max Constraints, make δ is refreshed immediately after each update. c (e) Avoiding non-physical crack tip overspeed propagation caused by continuous iteration; S63, Around the target set Ω tgt Expanding one layer to form a local domain Ω δ =ring(Ω) tgt ,1), in Ω δ The phase field degree of freedom φ, displacement degree of freedom u, and temperature degree of freedom Θ are unlocked internally, while a frozen state is applied to the far field region. The incremental solver only assembles the submatrices. M δ and the coupled residual vector R δ The incremental solution Δφ is obtained through a three-step loop of residual evaluation → linearization → iterative convergence. After the solution is completed, the energy release rate tensor G(e), e∈Ω, is updated immediately. δ And record the local residual r loc Used for subsequent convergence statistics; S64. Utilizing cell hash indexes Ω δ The generated new tensor G new Write back to the global tensor G, following the rule that new values ​​overwrite old values, and that unupdated cells remain unchanged. Perform a sparse compression operation after completion to ensure that indices are continuous and consistent with the graph mapping table. Consistent; S65. Return to step S3 with the latest energy field G and recalculate the node physical vector f. v (i) and the edge physical vector f e (ij), and then the temperature-force coupled message-passing graph network is run in S4 to obtain the updated crack initiation probability. With the probability of growth in the direction S5 calls the discriminator again to calculate the physical consistency score. This completes the closed-loop mapping of G→f→p→S; S66, Continuously track the score increase over two rounds And the iteration counter t, when |ΔS phys |<ε S or t≥N max When the system converges, the steady-state result package is exported as the final interface energy density table. Steady-state energy release rate field G ＊ Collaborative crack path Γ ＊ and its probability sequence and convergent physical consistency score 8. The method for predicting microcrack paths at multi-material gradient interfaces according to claim 1, characterized in that: Step S7 includes the following sub-steps: S71, Receive the network weights and physical field tensors from the final round of S6, and convert the node initial probability vectors. With edge expansion probability vector Perform memory locking and generate a read-only identifier ID. snap ; S72. Set a single threshold τ according to engineering safety specifications. src ,right Threshold truncation is performed to filter out the candidate node index set Λ cand Perform hierarchical clustering based on Euclidean distance within the physical coordinate system for Λ cand The resulting cluster set Extract centroid nodes sequentially Each cluster is merged into a single crack initiation location to avoid outputting redundant crack sources to downstream analysis modules; S73, Set of high-probability edges along the frozen area Constructing a polyline trajectory Π based on node order k Next, spline interpolation is used in the product's global coordinate system to transform Π k Smoothing is a quadratic continuous curve S74. Update the energy release rate tensor G synchronously. ＊ Remapping to the original reduced-order mesh index And save it as a list of fields. <CellID,G ＊ Then, the interface adapter is called to encapsulate the list into specific format entries that can be directly read by mainstream engineering CAE software; S75, Overall Crack Initiation Probability Path cumulative probability and final state discriminant score Constructing a multidimensional feature vector q k And calculate the overall reliability level ρ within the empirical Bayesian confidence model. k ∈{I,II,III}, this level is used as a note field and Binding; S76. Set the coordinates of the crack initiation point. Smooth crack curve Energy release rate field <CellID,G ＊ >and reliability label {ρ k } Summarized into a unified result package The result package automatically generates a JSON structure and a STEP geometry file in dual format, and includes a timestamp t. stamp and snapshot identifier ID snap To meet the requirements of version traceability and quality auditing; and to complete the delivery of the algorithm domain to the engineering application end.

Images