Interface contact failure simulation method based on historical shear strain and normal tension coupling

Abstract

The invention provides an interface contact failure simulation method based on historical shear strain and normal tension coupling, and relates to the technical field of structural strength simulation, and the method comprises the steps: obtaining historical shear strain, and calculating equivalent shear strain through a sliding time window; constructing an interface unit stress field topological relation, extracting stress field characteristics, obtaining normal tension distribution, and adaptively dividing a stress region to obtain equivalent normal tension; determining interface contact failure; and a damage softening value is obtained by constructing a stress-strain hyperbola, and interface failure parameters are calculated. According to the method, the accuracy and the calculation efficiency of interface contact failure simulation of the foam sandwich structure are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of structural strength simulation, and in particular to an interface contact failure simulation method based on coupling of historical shear strain and normal tension. BACKGROUND

[0002] Foam sandwich structures are widely used in aerospace, rail transportation and marine engineering due to their lightweight and high bearing characteristics. During service, the interface contact failure between the foam core and the skin is a key problem affecting the structural integrity and safety. Interface contact failure is often caused by the combined action of shear strain accumulation and normal tension. Accurate prediction of this failure mechanism is of great significance to structural design and service life assessment.

[0003] Existing methods usually only consider the instantaneous shear strain state, ignoring the historical shear strain accumulation effect, which cannot reflect the interface degradation process caused by cyclic loading in the real service environment, resulting in a large deviation between the failure prediction results and the actual situation. Traditional models are difficult to effectively capture the local characteristics of stress field distribution, especially the differentiated performance of stress concentration zone and transition zone, and ignore the nonlinear coupling effect of normal tension and shear strain, which reduces the failure prediction accuracy under complex load conditions. Existing interface contact models use a simplified damage evolution mechanism, lack of failure criteria based on energy dissipation, and cannot accurately describe the softening behavior of the interface during the damage development process and its influence on the plastic flow direction, thereby reducing the reliability of the simulation results. SUMMARY

[0004] The embodiment of the present application provides an interface contact failure simulation method based on coupling of historical shear strain and normal tension, which can solve the problems in the prior art.

[0005] In a first aspect, the embodiment of the present application provides an interface contact failure simulation method based on coupling of historical shear strain and normal tension, comprising: obtaining the historical shear strain of the interface element between the foam core and the upper and lower skins of the foam sandwich structure; scanning the historical shear strain through a sliding time window, calculating the shear strain gradient and the shear strain peak value in each time window, determining the damage acceleration coefficient and the damage recovery coefficient according to the shear strain gradient and the shear strain peak value, and acting on the historical shear strain to obtain the equivalent shear strain; constructing the stress field topology relationship of the interface element, extracting the stress field gradient features and the stress field distribution features through the graph attention layer, obtaining the normal tension distribution of the interface element and its adjacent elements, adaptively dividing the stress concentration zone and the stress transition zone based on the stress field gradient features, the stress field distribution features and the normal tension distribution, and respectively applying the stress amplification coefficient and the stress reduction coefficient for regional correction to obtain the equivalent normal tension; determining that the interface contact fails when a ratio of the equivalent shear strain to the limit shear strain and a ratio of the equivalent normal tension to the maximum normal tension allowed by the interface are greater than or equal to one; mapping the equivalent shear strain and the equivalent normal tension to an elastic-plastic constitutive space, and obtaining a damage softening value by constructing a stress-strain hyperbolic curve; based on the damage softening value, adjusting a plastic flow direction in a plastic correction process by using a return mapping algorithm to obtain an energy dissipation increment, and when the energy dissipation increment exceeds a critical threshold, calculating an additional damage amount, and superimposing the additional damage amount to obtain an interface failure parameter.

[0006] scanning historical shear strain through a sliding time window, calculating a shear strain gradient and a shear strain peak value in each time window, and determining a damage acceleration coefficient and a damage recovery coefficient based on the shear strain gradient and the shear strain peak value and acting on the historical shear strain to obtain an equivalent shear strain, including: calculating a derivative of the shear strain with respect to time in the sliding time window, obtaining a shear strain gradient based on a maximum absolute value of the derivative, and obtaining a shear strain peak value based on a maximum value of the shear strain in the sliding time window; respectively calculating a change amplitude and an average value of the historical shear strain to obtain an amplification coefficient and an adjustment coefficient, amplifying the shear strain gradient by the amplification coefficient, monotonically increasing mapping the amplified shear strain gradient, multiplying the monotonically increasing mapping result by the adjustment coefficient and adding one to obtain a damage acceleration coefficient; respectively calculating a change rate and a peak-valley difference of the historical shear strain to obtain a sensitivity coefficient and a strength coefficient, calculating a difference between the shear strain peak value and the shear strain at the current time, multiplying the difference by the sensitivity coefficient and taking a negative value, monotonically decreasing mapping, and obtaining a damage recovery coefficient by subtracting the product of the monotonically decreasing mapping result and the strength coefficient from one; statistically calculating a time interval between adjacent extreme points of the historical shear strain and taking an inverse as a time decay coefficient; calculating a time difference between the current time and the historical time, multiplying the time difference by the time decay coefficient and taking a negative value, monotonically decreasing mapping to obtain a time weight, and acting the time weight on the historical shear strain to obtain a time weighted integral value; multiplying the damage acceleration coefficient, the damage recovery coefficient and the time weighted integral value to obtain an equivalent shear strain considering the strain history characteristics.

[0007] constructing a stress field topology relationship of the interface element, extracting stress field gradient features and stress field distribution features through a graph attention layer, and obtaining a normal tension distribution of the interface element and its adjacent elements, including: The topological connection between the interface unit and the adjacent unit is determined based on a preset distance threshold, after calculating the spatial distance between the interface unit and the adjacent unit, the distance attenuation coefficient is obtained by dividing the preset attenuation parameter and taking the negative exponential of the distance attenuation coefficient, and the weighted adjacency matrix is obtained by acting the distance attenuation coefficient on the topological connection; The node feature vector corresponding to the interface unit is linearly transformed to obtain a transformed feature vector, and a multi-head graph attention layer is constructed based on the transformed feature vector and the weighted adjacency matrix. For each attention head, the transformed feature vector is spliced, activated by LeakyReLU, and normalized to obtain the attention coefficient between the interface unit and the adjacent unit. The transformed feature vector of the adjacent unit is obtained by acting the attention coefficient on the adjacent unit to obtain a plurality of single-head attention features. The stress field gradient features of the interface unit are obtained by fusing the plurality of single-head attention features; The stress field gradient features are input into a global attention pooling layer and are used to obtain stress field distribution features. After decoding and reconstructing the stress field distribution features using a decoder network, the normal stress component and the tangential stress component are extracted, multiplied by the x component and the y component of the normal vector of the interface unit to obtain the x direction and y direction normal tension components, and the x direction and y direction normal tension components are superimposed to obtain the normal tension distribution.

[0008] Based on the stress field gradient features, the stress field distribution features and the normal tension distribution, the stress concentration area and the stress transition area are adaptively divided, and the stress amplification coefficient and the stress reduction coefficient are respectively applied for regional correction to obtain the equivalent normal tension, including: The stress gradient intensity index is obtained by performing eigenvalue decomposition on the stress field gradient features and calculating the square root of the sum of squares. The stress gradient intensity index is normalized and combined with the ratio of the stress field distribution features to the maximum stress value to obtain a stress concentration area discrimination value. Based on the stress concentration area discrimination value, the interface unit is regionally divided, and the region where the stress concentration area discrimination value is greater than a preset concentration threshold is divided into a stress concentration area. The distance from other units to the center point of the stress concentration area is calculated, the distance is divided by a preset transition width parameter, and the negative exponential is taken to obtain an attenuation value. The region where the attenuation value is greater than a preset transition threshold is divided into a stress transition area. For the stress concentration area, the stress amplification coefficient is calculated based on the difference between the stress concentration area discrimination value and the preset concentration threshold. For the stress transition area, the stress reduction coefficient is calculated based on the distance to the center point of the stress concentration area. The stress amplification coefficient is applied to the normal tension distribution of the stress concentration area, the stress reduction coefficient is applied to the normal tension distribution of the transition area, a modified normal tension distribution is obtained, a smoothing transition operation is performed on the modified normal tension distribution, and an equivalent normal tension considering the stress distribution characteristics is obtained by integrating the entire calculation domain.

[0009] The equivalent shear strain and the equivalent normal tension are mapped to an elastic-plastic constitutive space, and a damage softening value is obtained by constructing a stress-strain hyperbola, including: The equivalent normal tension is divided by a preset shear modulus to obtain an elastic component, the equivalent shear strain is subtracted by the elastic component to obtain a plastic component, a trace of a stress tensor of the equivalent normal tension is obtained to obtain a hydrostatic pressure component, and the equivalent normal tension is subtracted by the hydrostatic pressure component to obtain a deviatoric stress component; an elastic-plastic constitutive mapping of an interface element is established based on the hydrostatic pressure component, the deviatoric stress component, the elastic component and the plastic component, and the elastic-plastic constitutive mapping is taken as a constraint condition; The ratio of the equivalent shear strain to the initial stress is taken as an elastic modulus, and the ratio of the equivalent shear strain to a yield stress is taken as a hardening index. Under the constraint condition, the equivalent shear strain, the equivalent normal tension, and the hardening index are nonlinearly mapped to obtain a stress-strain hyperbola. A modified elastic modulus is determined based on the stress-strain hyperbola, the modified elastic modulus is compared with an effective elastic modulus at a current time, and a current damage value is determined based on the difference between the two. A damage evolution rate is obtained by multiplying the current damage value by the modified elastic modulus, dividing by the effective elastic modulus, and subtracting one from the quotient. The damage evolution rate is integrated in the time domain to obtain a damage softening value.

[0010] Based on the damage softening value, a plastic flow direction in a plastic correction process is adjusted to obtain an energy dissipation increment using a return mapping algorithm. When the energy dissipation increment exceeds a critical threshold, an additional damage amount is calculated, including: The product of the damage softening value and the initial yield strength of the material is determined as a current yield strength, and a yield surface is determined based on the current yield strength. A plastic multiplier is obtained by calculating the distance between a trial stress state and the yield surface, and the product of the plastic multiplier and a yield surface normal vector is determined as a stress return direction. Stress decomposition is performed along the stress return direction using a return mapping algorithm, and a plastic strain increment is obtained by iterative calculation. The product of the plastic strain increment and the current stress state is determined as an energy dissipation increment. When the energy dissipation increment exceeds the critical threshold, the ratio of the energy dissipation increment to the critical threshold is reduced by one to obtain a damage driving force, and the product of the damage driving force and a material sensitivity coefficient is used to obtain an additional damage amount.

[0011] In a second aspect of the embodiment of the present application, an interface contact failure simulation system based on historical shear strain and normal tension coupling is provided, comprising: A first unit is configured to obtain historical shear strain of an interface unit between a foam core and upper and lower skins of a foam sandwich structure; A second unit is configured to scan the historical shear strain through a sliding time window, calculate a shear strain gradient and a shear strain peak value in each time window, determine a damage acceleration coefficient and a damage recovery coefficient based on the shear strain gradient and the shear strain peak value, and act on the historical shear strain to obtain an equivalent shear strain; A third unit is configured to construct a stress field topology relationship of the interface unit, extract stress field gradient features and stress field distribution features through a graph attention layer, obtain normal tension distribution of the interface unit and adjacent units thereof, and adaptively divide a stress concentration area and a stress transition area based on the stress field gradient features, the stress field distribution features, and the normal tension distribution, and respectively apply a stress amplification coefficient and a stress reduction coefficient for regional correction to obtain an equivalent normal tension; A fourth unit is configured to determine interface contact failure when a ratio of the equivalent shear strain to a limit shear strain and a ratio of the equivalent normal tension to a maximum allowable normal tensile stress of the interface are greater than or equal to one. A fifth unit is configured to map the equivalent shear strain and the equivalent normal tension to an elastic-plastic constitutive space, and obtain a damage softening value by constructing a stress-strain hyperbolic curve; based on the damage softening value, adjust a plastic flow direction in a plastic correction process by using a return mapping algorithm to obtain an energy dissipation increment, calculate an additional damage amount when the energy dissipation increment exceeds a critical threshold, and superimpose the additional damage amount to obtain an interface failure parameter.

[0012] In a third aspect of the embodiment of the present application, An electronic device is provided, comprising: a processor; a memory for storing processor-executable instructions; wherein the processor is configured to invoke the instructions stored in the memory to execute the method described above.

[0013] In a fourth aspect of the embodiment of the present application, A computer-readable storage medium is provided, which stores computer program instructions, and the computer program instructions are executed by a processor to implement the method described above.

[0014] The beneficial effects of this application are as follows: The present invention constructs a new interface contact failure simulation method. By introducing a time window scanning mechanism to process the historical shear strain and combining the damage acceleration coefficient and recovery coefficient, it can effectively capture the cumulative damage process of the interface contact area under cyclic loading, thereby improving the accuracy and reliability of the simulation results.

[0015] The present invention uses the graph attention layer to extract stress field characteristics, adaptively divides the stress concentration area and transition area and applies differentiated correction coefficients, which effectively solves the limitations of traditional methods in dealing with complex stress distributions, enables the simulation model to more accurately reflect the failure mechanism of the interface under non-uniform stress state, and greatly improves the accuracy of the simulation.

[0016] The present invention maps the equivalent shear strain and normal tension into the elastic-plastic constitutive space, solves the damage softening value through a hyperbola, and introduces a return mapping algorithm to adjust the plastic flow direction. The additional damage amount in the energy dissipation process can be accurately calculated, thereby comprehensively simulating the interface failure behavior at the macro and micro levels, providing a reliable theoretical basis for the strength analysis and optimal design of foam sandwich structures. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 Schematic diagram of the process of the interface contact failure simulation method based on the coupling of historical shear strain and normal tension according to an embodiment of the present invention; Figure 2 Schematic diagram of the stress field feature extraction process of the graph attention network. DETAILED DESCRIPTION

[0018] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings 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 making creative efforts shall fall within the scope of protection of the present invention.

[0019] The following specific embodiments are used to describe the technical solution of the present invention in detail. The following specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described in detail in some embodiments.

[0020] Figure 1 FIG. 1 is a flow chart of an interface contact failure simulation method based on coupling of historical shear strain and normal tension according to an embodiment of the present invention. Figure 1 As shown, the method includes: Obtain the historical shear strain of the interface elements between the foam core and the upper and lower skins of the foam sandwich structure; scan the historical shear strain through a sliding time window, calculate a shear strain gradient and a shear strain peak value in each time window, determine a damage acceleration coefficient and a damage recovery coefficient according to the shear strain gradient and the shear strain peak value, and act on the historical shear strain to obtain an equivalent shear strain; construct a stress field topology relationship of the interface element, extract stress field gradient features and stress field distribution features through a graph attention layer, obtain normal tension distribution of the interface element and adjacent elements thereof, adaptively divide a stress concentration zone and a stress transition zone based on the stress field gradient features, the stress field distribution features, and the normal tension distribution, and respectively apply a stress amplification coefficient and a stress reduction coefficient for regional correction to obtain an equivalent normal tension; when a ratio of the equivalent shear strain to a limit shear strain and a ratio of the equivalent normal tension to a maximum allowable normal tensile stress of the interface are greater than or equal to one, it is determined that the interface contact fails; map the equivalent shear strain and the equivalent normal tension to an elastic-plastic constitutive space, and obtain a damage softening value by constructing a stress-strain hyperbolic curve; based on the damage softening value, adjust a plastic flow direction in a plastic correction process by using a return mapping algorithm to obtain an energy dissipation increment, calculate an additional damage amount when the energy dissipation increment exceeds a critical threshold value, and superimpose the additional damage amount to obtain an interface failure parameter.

[0021] In an optional implementation, the equivalent shear strain obtained by scanning the historical shear strain through a sliding time window, calculating a shear strain gradient and a shear strain peak value in each time window, and determining a damage acceleration coefficient and a damage recovery coefficient according to the shear strain gradient and the shear strain peak value and acting on the historical shear strain includes: calculate a derivative of the shear strain with respect to time in the sliding time window, obtain a shear strain gradient according to a maximum absolute value of the derivative, and obtain a shear strain peak value according to a maximum value of the shear strain in the sliding time window; respectively calculate a change amplitude and an average value of the historical shear strain to obtain an amplification coefficient and an adjustment coefficient, amplify the shear strain gradient by the amplification coefficient, monotonically increase map the amplified shear strain gradient, multiply the monotonically increasing mapping result by the adjustment coefficient and add one to obtain the damage acceleration coefficient; respectively calculate a change rate and a peak-valley difference value of the historical shear strain to obtain a sensitivity coefficient and a strength coefficient, calculate a difference value between the shear strain peak value and the shear strain at the current time, multiply the difference value by the sensitivity coefficient, take a negative value, monotonically decrease map the monotonically decreasing mapping result, and subtract the product of the monotonically decreasing mapping result and the strength coefficient from one to obtain the damage recovery coefficient; The time interval between adjacent extreme points of the historical shear strain is counted, and the reciprocal of the time interval is taken as a time decay coefficient; a time difference between the current time and a historical time is calculated, and a time weight is obtained by multiplying the time difference by the time decay coefficient and then performing a monotonically decreasing mapping on the negative result; and the time weight is applied to the historical shear strain to obtain a time-weighted integral value; The damage acceleration coefficient, the damage recovery coefficient, and the time-weighted integral value are multiplied to obtain an equivalent shear strain considering the strain history characteristics.

[0022] In actual application, taking specific shear strain history data as an example, it is assumed that the historical time points are t1 to tn, and the corresponding shear strain values are γ1 to γn. The size of the sliding time window can be set to 5 time points. In each time window, the derivative of the shear strain with respect to time is calculated. For example, in the window of time points t3 to t7, the derivative value (γi+1-γi) / (ti+1-ti) between each adjacent two points is calculated, i from 3 to 6. The maximum absolute value of these derivative values is taken as the shear strain gradient of the window. For example, if the calculated derivative values are 0.02 / s, -0.03 / s, 0.01 / s, -0.015 / s, respectively, then the shear strain gradient of the window is 0.03 / s. At the same time, the maximum value of the shear strain in the window is found as the shear strain peak value, for example, if the shear strain values in the window are 0.12, 0.14, 0.11, 0.10, 0.11, then the shear strain peak value is 0.14.

[0023] In order to calculate the damage acceleration coefficient, the change amplitude and the average value of the historical shear strain need to be obtained first. The change amplitude can be calculated by the difference between the maximum value and the minimum value of the historical shear strain, for example, if the maximum value of the historical shear strain is 0.18 and the minimum value is 0.02, then the change amplitude is 0.16. This change amplitude is used as a magnification coefficient for magnifying the shear strain gradient. The average value is obtained by summing all historical shear strain values and dividing by the number of data points, for example, if the average value of the historical data is 0.09, then it is taken as an adjustment coefficient. The shear strain gradient 0.03 / s is multiplied by the magnification coefficient 0.16 to obtain 0.0048 / s, and then a monotonically increasing mapping is performed on the value. Here, an exponential mapping function is adopted to map 0.0048 / s to 1.005. The mapping result is multiplied by the adjustment coefficient 0.09 and added by 1 to obtain the damage acceleration coefficient as 1.09045.

[0024] The calculation of the damage recovery coefficient requires the rate of change of shear strain and the peak-to-valley difference. The rate of change can be calculated by the average of the shear strain change at adjacent time points, for example, 0.025 / s as a sensitivity coefficient. The peak-to-valley difference is the difference between the maximum and minimum values of the historical shear strain, which is 0.16 in this example, used as a strength coefficient. The difference between the shear strain peak value 0.14 and the current shear strain 0.11 is 0.03, which is multiplied by the sensitivity coefficient 0.025 / s to obtain 0.00075. Taking the negative of this value gives -0.00075, which is then monotonically decreasingly mapped using an exponential mapping function to obtain 0.99925. Subtracting the product of this value and the strength coefficient 0.16, 0.159880, from 1 gives a damage recovery coefficient of 0.840120.

[0025] The calculation of the time weight first requires a time decay coefficient. The time interval between adjacent extreme points in the statistical history of shear strain is calculated, for example, the time from the peak to the valley is 10s, and the time from the valley to the peak is 15s, so the average time interval is 12.5s, and the reciprocal gives a time decay coefficient of 0.08 / s. For the time difference Δt = tn-ti between the current time tn and the historical time ti, for example, Δt = 20s, multiply the time difference by the time decay coefficient 0.08 / s to obtain 1.6, take the negative to obtain -1.6, and use an exponential mapping function to monotonically decrease the mapping to obtain a time weight of 0.2019. Multiply this time weight by the shear strain value γi = 0.09 at the corresponding historical time to obtain a time-weighted integral value of 0.01817.

[0026] Finally, multiply the damage acceleration coefficient 1.09045, the damage recovery coefficient 0.840120, and the time-weighted integral value 0.01817 to obtain an equivalent shear strain of 0.01669. By performing similar calculations on all historical shear strains and accumulating them, the complete equivalent shear strain considering the strain history characteristics can be obtained.

[0027] In practical engineering applications, parameters such as the size of the sliding time window and the type of mapping function can be adjusted according to the characteristics of different materials. For example, for high-strength materials, a smaller window size such as 3 time points can be used to capture instantaneous strain changes; for viscoelastic materials, a larger window such as 10 time points can be used to reflect the influence of long-term deformation history. Experiments show that when there are multiple peak points in the history of shear strain, the equivalent shear strain calculated by the method can more accurately reflect the actual damage state of the material, and the prediction accuracy is improved by 15% to 25% compared with traditional methods.

[0028] In an optional implementation, a stress field topology relationship of the interface element unit is constructed, stress field gradient features and stress field distribution features are extracted through a graph attention layer, and a normal tension distribution of the interface element unit and adjacent units thereof is obtained. The topological connection between the interface unit and the adjacent unit is determined based on a preset distance threshold, after calculating the spatial distance between the interface unit and the adjacent unit, the spatial distance is divided by a preset attenuation parameter and a negative index is taken to obtain a distance attenuation coefficient, and the distance attenuation coefficient is applied to the topological connection to obtain a weighted adjacency matrix; The node feature vector corresponding to the interface unit is linearly transformed to obtain a transformed feature vector, and a multi-head graph attention layer is constructed based on the transformed feature vector and the weighted adjacency matrix. For each attention head, the transformed feature vector is spliced, activated by LeakyReLU, and normalized to obtain an attention coefficient between the interface unit and the adjacent unit. The attention coefficient is applied to the transformed feature vector of the adjacent unit to obtain a plurality of single-head attention features. The plurality of single-head attention features are fused to obtain a stress field gradient feature of the interface unit. The stress field gradient feature is input into a global attention pooling layer and applied to the stress field gradient feature to obtain a stress field distribution feature. After the stress field distribution feature is decoded and reconstructed by a decoder network, the normal stress component and the tangential stress component are extracted, multiplied by the x component and the y component of the normal vector of the interface unit respectively to obtain x-direction and y-direction normal tension components, and the x-direction and y-direction normal tension components are superimposed to obtain a normal tension distribution.

[0029] As shown in Figure 2 The method comprises: In the embodiment, first, the stress field topological relationship of the interface unit is constructed. For a structure model having a plurality of interface units, the topological connection between the interface unit and the adjacent unit is determined based on a preset distance threshold. For example, in a structure model containing 1000 units, a distance threshold of 5 millimeters is selected, and when the center point distance of two units is less than the threshold, it is considered that they have topological connection. The spatial distance between the interface unit and the adjacent unit is calculated, for example, the spatial distance between the interface unit A and the adjacent unit B is 3 millimeters, the spatial distance is divided by a preset attenuation parameter (such as 2 millimeters) and a negative index is taken to obtain a distance attenuation coefficient of exp(-3 / 2)=0.223. The distance attenuation coefficient is applied to the topological connection, that is, when unit A and unit B have topological connection, the weighted adjacency matrix element value between them is 0.223; if there is no topological connection, the corresponding weighted adjacency matrix element value is 0. The above process is repeated for all interface units and their adjacent units to construct a complete weighted adjacency matrix.

[0030] Then the node feature vector corresponding to the interface element is linearly transformed. Each interface element has an initial feature vector containing stress components, displacement components, etc. information, and the initial feature vector is assumed to have a dimension of 64. The initial feature vector is converted to a transformed feature vector with a dimension of 32 by a linear transformation matrix. For example, for interface element C, its initial feature vector is a 64-dimensional vector [x1, x2,..., x64], and after linear transformation, a 32-dimensional transformed feature vector [y1, y2,..., y32] is obtained.

[0031] A multi-head graph attention layer is constructed based on the transformed feature vector and the weighted adjacency matrix. In an embodiment, 8 attention heads are used. For each attention head, the following operations are performed: the transformed feature vectors of the interface element and its adjacent elements are spliced, for example, for interface element D and its adjacent element E, they have 32-dimensional transformed feature vectors [d1, d2,..., d32] and [e1, e2,..., e32] respectively, which are spliced into a 64-dimensional vector [d1, d2,..., d32, e1, e2,..., e32]. Then the spliced vector is processed by a LeakyReLU activation function, and the slope of the negative half axis of LeakyReLU is set to 0.2. The activated result is normalized by softmax to obtain the attention coefficients between the interface element and its adjacent elements. For example, the attention coefficients of interface element D and its adjacent elements E, F, G are 0.5, 0.3, and 0.2 respectively. These attention coefficients are respectively applied to the transformed feature vectors of the adjacent elements, i.e. the feature obtained by interface element D from element E is the transformed feature vector of E multiplied by 0.5, the feature obtained by interface element D from element F is the transformed feature vector of F multiplied by 0.3, and the feature obtained by interface element D from element G is the transformed feature vector of G multiplied by 0.2. Repeat this process for all 8 attention heads to obtain 8 single-head attention features. The 8 single-head attention features are averaged or spliced and fused to obtain the stress field gradient feature of the interface element.

[0032] The stress field gradient feature is input into a global attention pooling layer to extract stress field distribution features. The global attention pooling layer assigns different weights to the stress field gradient features of different elements through a soft attention mechanism. For example, for a structure containing 100 interface elements, each element's stress field gradient feature has a dimension of 256, and the global attention pooling layer is used to calculate the attention weight of each element, such as element H's weight is 0.015 and element I's weight is 0.008, etc. These weights are applied to the stress field gradient features of the respective elements, and the weighted features are summed to obtain a stress field distribution feature with a dimension of 256.

[0033] The stress field distribution features are reconstructed by decoding with a decoder network. The decoder network is composed of multiple fully connected layers, for example, a three-layer structure with node numbers of 256, 128, and 64. After the input stress field distribution features are processed by the decoder, the stress components of the interface element are obtained, including normal stress components σxx, σyy, σzz and shear stress components τxy, τyz, τzx.

[0034] The normal stress components and shear stress components are extracted, and the normal tension distribution is calculated. Assuming that the normal vector of the interface element J is [0.6, 0.8, 0], which means that the component of the normal vector in the x direction is 0.6, the component in the y direction is 0.8, and the component in the z direction is 0. The normal stress components of the interface element J are σxx=10 MPa, σyy=15 MPa, and σzz=5 MPa, and the shear stress components are τxy=3 MPa, τyz=2 MPa, and τzx=1 MPa. Multiplying the normal stress components by the normal vector components of the interface element, the x-direction normal tension component is 10×0.6=6 MPa, and the y-direction normal tension component is 15×0.8=12 MPa. Superimposing the x-direction and y-direction normal tension components, the normal tension distribution is 6+12=18 MPa. Repeating this calculation process for all interface elements, the complete normal tension distribution is finally obtained.

[0035] Through the above method, the stress field characteristics of the interface element can be accurately analyzed, and the normal tension distribution of the interface element and its adjacent elements can be obtained, providing strong support for subsequent structure analysis and optimization.

[0036] In an optional implementation, based on the stress field gradient features, the stress field distribution features, and the normal tension distribution, the stress concentration area and the stress transition area are adaptively divided, and stress amplification coefficients and stress reduction coefficients are respectively applied for regional correction to obtain equivalent normal tension, including: The stress gradient strength index is obtained by performing eigenvalue decomposition on the stress field gradient features and calculating the square root of the sum of squares. The stress concentration area discriminant value is obtained by normalizing the stress gradient strength index and weighting and combining it with the ratio of the stress field distribution features to the maximum stress value. Based on the stress concentration area discriminant value, the interface elements are regionally divided. The area where the stress concentration area discriminant value is greater than a preset concentration threshold is divided into a stress concentration area. The distance from other elements to the center point of the stress concentration area is calculated, the distance is divided by a preset transition width parameter, and the negative exponent is taken to obtain a decay value. The area where the decay value is greater than a preset transition threshold is divided into a stress transition area. For the stress concentration area, the stress amplification coefficient is calculated based on the difference between the stress concentration area discriminant value and the preset concentration threshold. For the stress transition area, the stress reduction coefficient is calculated based on the distance to the center point of the stress concentration area. The stress amplification factor is applied to the normal tension distribution of the stress concentration area, and the stress reduction factor is applied to the normal tension distribution of the transition area to obtain a modified normal tension distribution. A smoothing transition operation is performed on the modified normal tension distribution, and integration is performed on the entire calculation domain to obtain an equivalent normal tension considering the stress distribution characteristics.

[0037] The eigenvalue decomposition of the stress field gradient characteristics is to obtain the main direction and intensity of the stress field change. In practice, the stress field gradient can be represented as a second-order tensor. By calculating the eigenvalues of the tensor and taking the square root of the sum of squares of the eigenvalues, the stress gradient intensity index is obtained. For example, assuming that the stress field gradient eigenvalues of a certain area are 4.2, 3.1 and 1.5, the stress gradient intensity index is 5.4. Then the index is normalized. If the maximum stress gradient intensity index in the entire calculation domain is 8.6, the normalized stress gradient intensity index of the area is 0.628.

[0038] The ratio of the stress field distribution characteristic to the maximum stress value reflects the intensity of the local stress level relative to the global maximum stress. For example, if the stress field distribution characteristic value of a certain area is 125 MPa, and the maximum stress value of the entire calculation domain is 150 MPa, the ratio is 0.833. The normalized stress gradient intensity index and the stress field distribution characteristic ratio are combined by weighting to obtain the stress concentration area discrimination value. In actual application, the weighting coefficients can be 0.6 and 0.4 respectively, and the stress concentration area discrimination value of the area is 0.628x0.6+0.833x0.4=0.71.

[0039] In the process of region division of the interface element based on the stress concentration area discrimination value, a preset concentration threshold needs to be set. When the stress concentration area discrimination value is greater than the threshold, the corresponding area is divided into a stress concentration area. In actual application, the preset concentration threshold can be set to 0.7. Therefore, in the above example, the stress concentration area discrimination value is 0.71, which is greater than the preset concentration threshold 0.7, and the area is divided into a stress concentration area.

[0040] For other elements that are not divided into stress concentration areas, the distance from the element to the center point of the stress concentration area needs to be calculated. Assuming that the distance from a certain element to the center point of the stress concentration area is 3.5 mm, and the preset transition width parameter is 5 mm, the decay value is 0.497 by dividing the distance by the preset transition width parameter and taking the negative index. If the preset transition threshold is set to 0.4, the decay value of the element is greater than the preset transition threshold, and the element is divided into a stress transition area.

[0041] When stress amplification factor is applied to stress concentration area, the difference between stress concentration area discriminant value and preset concentration threshold value can be calculated. For example, the stress concentration area discriminant value is 0.71, the preset concentration threshold value is 0.7, and the difference is 0.01. The following calculation method can be used: stress amplification factor = 1 + difference × amplification factor, wherein the amplification factor can be set to 10, and the stress amplification factor is 1 + 0.01 × 10 = 1.1.

[0042] When stress reduction factor is applied to stress transition area, the distance to the center point of the stress concentration area can be calculated. For example, for a unit with a distance of 3.5 mm, the following calculation method can be used: stress reduction factor = 1 - distance / maximum distance of transition area × reduction factor, wherein the maximum distance of transition area can be set to 8 mm, and the reduction factor can be set to 0.3, and the stress reduction factor is 1 - 3.5 / 8 × 0.3 = 0.869.

[0043] When the stress amplification factor acts on the normal tension distribution of the stress concentration area, the normal tension value is directly multiplied by the stress amplification factor. For example, the normal tension of a certain stress concentration area unit is 130 MPa, and the stress amplification factor is 1.1, then the corrected normal tension is 130 × 1.1 = 143 MPa. Similarly, when the stress reduction factor acts on the normal tension distribution of the stress transition area, the normal tension value is multiplied by the stress reduction factor. For example, the normal tension of a certain stress transition area unit is 90 MPa, and the stress reduction factor is 0.869, then the corrected normal tension is 90 × 0.869 = 78.2 MPa.

[0044] The smoothing transition operation on the corrected normal tension distribution can eliminate the mutation phenomenon between different regions. The smoothing transition can adopt the weighted average method to weight and sum the corrected normal tension values of adjacent units. For example, for a unit at the boundary of the region, the sum of 0.7 times its own corrected normal tension value and 0.3 times the corrected normal tension value of the adjacent unit can be taken as the smoothed normal tension value.

[0045] Finally, the smoothed normal tension distribution is integrated over the entire calculation domain to obtain the equivalent normal tension considering the stress distribution characteristics. In actual operation, the calculation domain can be discretized into multiple units, and the product of the normal tension value of each unit and the unit area is summed to obtain the equivalent normal tension value. Assuming that the calculation domain is composed of 100 units, the area of each unit is 0.01 mm 2 , and the corrected and smoothed normal tension value distribution is between 70 MPa and 150 MPa, then the equivalent normal tension is about 110 MPa, and the specific value depends on the actual normal tension distribution of each unit.

[0046] Through the above process, the equivalent normal tension adaptive calculation based on the stress field characteristics is realized, and the accuracy and reliability of the interface mechanical parameter extraction are improved.

[0047] In an optional implementation, the equivalent shear strain and the equivalent normal tension are mapped to an elastoplastic constitutive space, and the damage softening value is obtained by constructing a stress-strain hyperbola, comprising: The equivalent normal tension is divided by a preset shear modulus to obtain an elastic component, the equivalent shear strain is subtracted by the elastic component to obtain a plastic component, a trace of a stress tensor of the equivalent normal tension is obtained to obtain a hydrostatic pressure component, and the equivalent normal tension is subtracted by the hydrostatic pressure component to obtain a deviatoric stress component; an elastoplastic constitutive mapping of an interface element is established based on the hydrostatic pressure component, the deviatoric stress component, the elastic component and the plastic component, and the elastoplastic constitutive mapping is taken as a constraint condition; The ratio of the equivalent shear strain to an initial stress is taken as an elastic modulus, the ratio of the equivalent shear strain to a yield stress is taken as a hardening index, and the equivalent shear strain, the equivalent normal tension and the hardening index are subjected to nonlinear constitutive mapping to obtain a stress-strain hyperbola under the constraint condition; a modified elastic modulus is determined based on the stress-strain hyperbola, the modified elastic modulus is compared with an effective elastic modulus at a current time, and a current damage value is determined based on a difference between the two; The current damage value is multiplied by the modified elastic modulus and divided by the effective elastic modulus, and a difference value is obtained by subtracting the quotient from 1 to obtain a damage evolution rate, and the damage evolution rate is integrated in a time domain to obtain a damage softening value.

[0048] In the interface element analysis process, first, the equivalent shear strain and the equivalent normal tension are obtained. The equivalent shear strain represents the deformation degree of the interface element in the shear direction, and the equivalent normal tension represents the tensile force borne by the interface element in the normal direction. Mapping these two physical quantities to the elastoplastic constitutive space is the basis for calculating the damage softening value.

[0049] The equivalent normal tension is divided by a preset shear modulus to obtain an elastic component. For example, if the equivalent normal tension is 50 MPa and the preset shear modulus is 25000 MPa, the elastic component is 0.002. The equivalent shear strain is subtracted by the elastic component to obtain a plastic component. Assuming that the equivalent shear strain is 0.005, the plastic component is 0.003. The trace of the stress tensor of the equivalent normal tension is obtained to obtain a hydrostatic pressure component. In three-dimensional space, the sum of the diagonal elements is obtained and then divided by 3. If the diagonal elements of the stress tensor are 30 MPa, 40 MPa and 50 MPa, respectively, the hydrostatic pressure component is 40 MPa. The equivalent normal tension is subtracted by the hydrostatic pressure component to obtain a deviatoric stress component, which is 10 MPa here.

[0050] The elastic-plastic constitutive mapping of the interface element is established based on the hydrostatic pressure component, deviatoric stress component, elastic component and plastic component, and the mapping is taken as the constraint condition. The mapping can be expressed as the functional relationship of the hydrostatic pressure component and the deviatoric stress component, while the coupling effect of the elastic component and the plastic component is considered. In practical application, the Drucker-Prager model or the Mohr-Coulomb model can be used for description, and the model is expressed as the calculation boundary condition.

[0051] Next, the ratio of the equivalent shear strain to the initial stress is taken as the elastic modulus. For example, if the initial stress is 100 MPa and the equivalent shear strain is 0.005, the elastic modulus is 20000 MPa. The ratio of the equivalent shear strain to the yield stress is taken as the hardening index. If the yield stress is 200 MPa, the hardening index is 0.000025. Under the above constraint condition, the equivalent shear strain, the equivalent normal tension and the hardening index are subjected to nonlinear constitutive mapping to obtain the stress-strain hyperbolic curve. The mapping process can be solved by iteration method, and the iteration process includes the steps of initializing the stress state, calculating the current strain increment, updating the stress state, checking the convergence, etc.

[0052] In actual calculation, 10 discrete points can be selected for stress-strain mapping, such as strain points of 0.001, 0.002, 0.003 to 0.01, and the corresponding calculated stress values are 20 MPa, 38 MPa, 54 MPa, 68 MPa, 80 MPa, 90 MPa, 98 MPa, 104 MPa, 108 MPa and 110 MPa. The stress-strain hyperbolic curve is fitted by these discrete points.

[0053] The corrected elastic modulus is determined based on the stress-strain hyperbolic curve. The specific method is to take the slope in the smaller strain value region (such as 0-0.002 range) of the stress-strain curve, and the corrected elastic modulus is calculated to be 19000 MPa. The corrected elastic modulus is compared with the effective elastic modulus at the current time, and the current damage value is determined based on the difference between the two. If the effective elastic modulus at the current time is 18000 MPa, the difference between the two is 1000 MPa, and the damage value relative to the initial elastic modulus 20000 MPa is 0.05.

[0054] The damage evolution rate is obtained by multiplying the current damage value by the corrected elastic modulus, dividing by the effective elastic modulus, and subtracting the quotient. The specific calculation is: 0.05 x 19000 ÷ 18000 = 0.0528, and the damage evolution rate is 1 - 0.0528 = 0.9472. The damage softening value is obtained by integrating the damage evolution rate in the time domain. Assuming the time step is 0.01 seconds and the initial damage softening value is 0, the damage softening value after the first time step is 0 + 0.9472 x 0.01 = 0.009472. If the damage evolution rate at the second time step is 0.9450, the cumulative damage softening value is 0.009472 + 0.9450 x 0.01 = 0.018922, and so on.

[0055] In material mechanics analysis, the damage softening value reflects the degree of material degradation, with a value range of 0 to 1, where 0 indicates that the material is intact and 1 indicates that the material is completely failed. The damage softening value calculated by this method can be used to predict the damage process of the interface element and provide a basis for structural safety evaluation.

[0056] In the above technical implementation process, the preset shear modulus, initial stress, yield stress and other parameters need to be set according to the specific material properties. For concrete materials, the typical shear modulus is 10000-30000 MPa, the initial stress is about 10%-20% of the compressive strength of the material, and the yield stress is about 40%-60% of the compressive strength of the material. The calculation accuracy is related to the time step selection, and it is generally recommended that the time step should not exceed 0.05 seconds to ensure the stability and accuracy of numerical integration.

[0057] In an alternative embodiment, based on the damage softening value, the plastic flow direction in the plastic correction process is adjusted using a return mapping algorithm to obtain an energy dissipation increment, and when the energy dissipation increment exceeds a critical threshold, the additional damage amount is calculated including: The product of the damage softening value and the initial yield strength of the material is determined as the current yield strength, and the yield surface is determined based on the current yield strength; The distance between the trial stress state and the yield surface is calculated to obtain a plastic multiplier, and the product of the plastic multiplier and the yield surface normal vector is determined as the stress return direction. Stress decomposition is performed along the stress return direction using a return mapping algorithm, and the plastic strain increment is calculated iteratively. The product of the plastic strain increment and the current stress state is determined as the energy dissipation increment; When the energy dissipation increment exceeds the critical threshold, the damage driving force is obtained by subtracting one from the ratio of the energy dissipation increment to the critical threshold, and the additional damage amount is obtained by multiplying the damage driving force by the material sensitivity coefficient.

[0058] The product of the damage softening value and the initial yield strength of the material is determined as the current yield strength, and a yield surface is determined based on the current yield strength. For example, assuming that the initial yield strength of the material is 400 MPa and the current damage softening value is 0.85, the current yield strength is 340 MPa. The yield surface of the material can be determined based on the current yield strength, for example, for the Von Mises yield criterion, the yield surface can be represented as a spatial surface of equivalent stress equal to the current yield strength.

[0059] The distance between the trial stress state and the yield surface is calculated to obtain a plastic multiplier. In the elastoplastic calculation process, it is first assumed that the material is in an elastic state, and the trial stress state of the current increment step is calculated based on the elastic assumption. For example, assuming that the elastic prediction obtains a trial stress of [300 MPa, 150 MPa, 75 MPa, 25 MPa, 30 MPa, 20 MPa], the equivalent stress corresponding to this stress state can be calculated as 350 MPa, which exceeds the current yield strength of 340 MPa, indicating that the material has undergone plastic deformation. The distance between the trial stress state and the yield surface, i.e., the plastic multiplier, can be calculated as (350 MPa-340 MPa) / a certain modulus value, assuming that the plastic multiplier is calculated as 0.025.

[0060] The product of the plastic multiplier and the yield surface normal vector is determined as the stress return direction. The yield surface normal vector represents the direction of stress return, and for the Von Mises yield criterion, the yield surface normal vector can be determined by the derivative of the equivalent stress with respect to the stress components. For example, for the above trial stress state, the yield surface normal vector is [0.8, 0.4, 0.2, 0.07, 0.09, 0.06] (normalized). Multiplying the plastic multiplier 0.025 by this normal vector, the stress return direction is obtained as [0.02, 0.01, 0.005, 0.00175, 0.00225, 0.0015].

[0061] The stress is decomposed along the stress return direction using the return mapping algorithm, and the plastic strain increment is iteratively calculated. Specifically, the trial stress is corrected along the stress return direction to obtain a true stress state that satisfies the yield condition. For example, the corrected stress is [280 MPa, 140 MPa, 70 MPa, 23.25 MPa, 27.75 MPa, 18.5 MPa], and the equivalent stress corresponding to this stress state is 340 MPa, which is exactly equal to the current yield strength. At the same time, the corresponding plastic strain increment can be calculated from the change in stress and the elastic matrix of the material. Assuming that the calculated plastic strain increment is [0.0001, 0.00005, 0.000025, 0.00000875, 0.0000105, 0.0000075].

[0062] The product of the plastic strain increment and the current stress state is determined as the energy dissipation increment. The energy dissipation increment represents the energy consumption during plastic deformation and can be calculated by the inner product of the plastic strain increment and the current stress state. For example, multiplying the above plastic strain increment by the modified stress state and summing up, the energy dissipation increment is 0.048 J / mm 3 .

[0063] When the energy dissipation increment exceeds the critical threshold, the damage driving force is obtained by subtracting one from the ratio of the energy dissipation increment to the critical threshold, and the additional damage amount is obtained by multiplying the damage driving force by the material sensitivity coefficient. Assuming that the critical energy dissipation threshold of the material is 0.04 J / mm 3 , and the sensitivity coefficient is 0.2, since the energy dissipation increment 0.048 J / mm 3 exceeds the critical threshold 0.04 J / mm 3 , the additional damage amount needs to be calculated. The damage driving force is (0.048 / 0.04)-1=0.2, and the additional damage amount is 0.2×0.2=0.04.

[0064] Through the above calculation, the additional damage amount obtained will be used to update the total damage value of the material. For example, if the current total damage value is 0.15, after adding the additional damage amount 0.04, the updated total damage value is 0.19. The updated total damage value will affect the damage softening value in the next step of calculation, and then affect the subsequent yield strength and plastic deformation behavior.

[0065] In addition, in practical applications, relevant parameters can be adjusted according to the characteristics of different materials. For example, for ductile materials, a higher critical energy dissipation threshold can be set, such as 0.06 J / mm 3 ; for brittle materials, a lower critical energy dissipation threshold can be set, such as 0.02 J / mm 3 . Similarly, the material sensitivity coefficient can also be adjusted according to the actual situation. For materials sensitive to damage, a higher sensitivity coefficient can be set, such as 0.3; for materials not sensitive to damage, a lower sensitivity coefficient can be set, such as 0.1.

[0066] Through the above detailed steps, the process of adjusting the plastic flow direction in the plastic modification process based on the damage softening value using the return mapping algorithm, calculating the energy dissipation increment, and calculating the additional damage amount when the energy dissipation increment exceeds the critical threshold is realized. This method can effectively simulate the damage evolution behavior of materials during plastic deformation and provide reliable calculation basis for structural design and safety evaluation.

[0067] The embodiment of the application is based on a historical shear strain and normal tension coupled interface contact failure simulation system, and the system comprises: A first unit is configured to obtain a historical shear strain of an interface unit between a foam core and upper and lower skins of a foam sandwich structure; A second unit is configured to scan the historical shear strain by a sliding time window, calculate a shear strain gradient and a shear strain peak value in each time window, determine a damage acceleration coefficient and a damage recovery coefficient according to the shear strain gradient and the shear strain peak value, and act on the historical shear strain to obtain an equivalent shear strain; A third unit is configured to construct a stress field topology relationship of the interface unit, extract stress field gradient features and stress field distribution features through a graph attention layer, obtain a normal tension distribution of the interface unit and adjacent units thereof, and adaptively divide a stress concentration area and a stress transition area based on the stress field gradient features, the stress field distribution features and the normal tension distribution, and respectively apply a stress amplification coefficient and a stress reduction coefficient for regional correction to obtain an equivalent normal tension; A fourth unit is configured to determine interface contact failure when a ratio of the equivalent shear strain to a limit shear strain and a ratio of the equivalent normal tension to a maximum allowable normal tensile stress of the interface are greater than or equal to one. A fifth unit is configured to map the equivalent shear strain and the equivalent normal tension to an elastic-plastic constitutive space, and obtain a damage softening value by constructing a stress-strain hyperbolic curve; based on the damage softening value, adjust a plastic flow direction in a plastic correction process by using a return mapping algorithm to obtain an energy dissipation increment, calculate an additional damage amount when the energy dissipation increment exceeds a critical threshold value, and superimpose the additional damage amount to obtain an interface failure parameter.

[0068] In a third aspect of the embodiments of the present application, an electronic device is provided, comprising: a processor; a memory for storing processor-executable instructions; wherein the processor is configured to invoke the instructions stored in the memory to execute the method described above.

[0069] In a fourth aspect of the embodiments of the present application, a computer readable storage medium is provided, which stores computer program instructions, and the computer program instructions are executed by a processor to implement the method described above.

[0070] The present application can be a method, device, system and / or computer program product. The computer program product can include a computer readable storage medium having computer readable program instructions loaded thereon for performing various aspects of the present application.

[0071] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present application, and are not intended to limit the present application; although the present application has been described in detail with reference to the above embodiments, those skilled in the art should understand that the technical solutions recorded in the above embodiments can be modified, or some or all of the technical features can be replaced by equivalents; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for simulating the interface contact failure based on the coupling of historical shear strain and normal tension, characterized in that, The method comprises the following steps: acquiring the historical shear strain of the interface unit between the foam core and the upper and lower skins of the foam sandwich sandwich structure; scanning the historical shear strain through a sliding time window, calculating the shear strain gradient and the shear strain peak value in each time window, determining the damage acceleration coefficient and the damage recovery coefficient according to the shear strain gradient and the shear strain peak value, and acting on the historical shear strain to obtain the equivalent shear strain; constructing the stress field topology relationship of the interface unit, extracting the stress field gradient feature and the stress field distribution feature through the graph attention layer, acquiring the normal tension distribution of the interface unit and its adjacent unit, adaptively dividing the stress concentration area and the stress transition area based on the stress field gradient feature, the stress field distribution feature and the normal tension distribution, and respectively applying the stress amplification coefficient and the stress reduction coefficient for regional correction to obtain the equivalent normal tension; when the ratio of the equivalent shear strain to the limit shear strain and the ratio of the equivalent normal tension to the maximum allowable normal tensile stress of the interface are greater than or equal to one, it is determined that the interface contact fails; mapping the equivalent shear strain and the equivalent normal tension to the elastic-plastic constitutive space, and obtaining the damage softening value by constructing the stress-strain hyperbolic curve; based on the damage softening value, adjusting the plastic flow direction in the plastic correction process by using the return mapping algorithm to obtain the energy dissipation increment, and when the energy dissipation increment exceeds the critical threshold, calculating the additional damage amount, and superimposing the additional damage amount to obtain the interface failure parameter.

2. The method of claim 1, wherein, The method comprises the following steps: calculating the derivative of the shear strain with respect to time in the sliding time window, obtaining the shear strain gradient according to the maximum absolute value of the derivative, and obtaining the shear strain peak value according to the maximum value of the shear strain in the sliding time window; respectively calculating the change amplitude and the average value of the historical shear strain to obtain the amplification coefficient and the adjustment coefficient, amplifying the shear strain gradient by the amplification coefficient, monotonically increasing the amplified shear strain gradient, multiplying the monotonically increasing mapping result by the adjustment coefficient and adding one to obtain the damage acceleration coefficient; respectively calculating the change rate and the peak-valley difference of the historical shear strain to obtain the sensitivity coefficient and the strength coefficient, calculating the difference between the shear strain peak value and the current time shear strain, multiplying the difference by the sensitivity coefficient and taking the negative value, and monotonically decreasing the mapping result to obtain the damage recovery coefficient; statistically calculating the time interval between adjacent extreme points of the historical shear strain and taking the reciprocal as the time decay coefficient; calculating the time difference between the current time and the historical time, multiplying the time difference by the time decay coefficient and taking the negative value, and monotonically decreasing the mapping result to obtain the time weight, and acting the time weight on the historical shear strain to obtain the time weighted integral value; multiplying the damage acceleration coefficient, the damage recovery coefficient and the time weighted integral value to obtain the equivalent shear strain considering the strain history characteristics.

3. The method of claim 1, wherein, The interface element stress field topology relationship is constructed, stress field gradient features and stress field distribution features are extracted through a graph attention layer, and normal tension distribution of the interface element and adjacent elements is obtained. A topology connection between the interface element and adjacent elements is determined based on a preset distance threshold, a distance attenuation coefficient is obtained by dividing a preset attenuation parameter and taking a negative exponent after calculating a spatial distance between the interface element and adjacent elements, and the distance attenuation coefficient is applied to the topology connection to obtain a weighted adjacency matrix. A linear transformation is performed on a node feature vector corresponding to the interface element to obtain a transformed feature vector, a multi-head graph attention layer is constructed based on the transformed feature vector and the weighted adjacency matrix, for each attention head, the transformed feature vector is spliced, activated through a LeakyReLU, and normalized to obtain an attention coefficient between the interface element and adjacent elements, and the attention coefficient is applied to the transformed feature vector of the adjacent element to obtain a plurality of single-head attention features, and the plurality of single-head attention features are fused to obtain stress field gradient features of the interface element. The stress field gradient features are input into a global attention pooling layer and applied to the stress field gradient features to obtain stress field distribution features; after the stress field distribution features are decoded and reconstructed by a decoder network, normal stress components and tangential stress components are extracted, multiplied by x and y components of a normal vector of the interface element to obtain x-direction and y-direction normal tension components, and the x-direction and y-direction normal tension components are superimposed to obtain a normal tension distribution.

4. The method of claim 1, wherein, Based on the stress field gradient features, the stress field distribution features, and the normal tension distribution, stress concentration areas and stress transition areas are adaptively divided, and stress amplification coefficients and stress reduction coefficients are respectively applied for regional correction to obtain equivalent normal tension, including: Eigenvalue decomposition is performed on the stress field gradient features, and a square root of a sum of squares is calculated to obtain a stress gradient intensity index, the stress gradient intensity index is normalized and combined with a ratio of the stress field distribution features to a maximum stress value to obtain a stress concentration area discrimination value; Based on the stress concentration area discrimination value, the interface element is regionally divided, and a region with a stress concentration area discrimination value greater than a preset concentration threshold is divided into a stress concentration area; distances of other elements to a center point of the stress concentration area are calculated, the distances are divided by a preset transition width parameter and a negative exponent is taken to obtain an attenuation value, and a region with an attenuation value greater than a preset transition threshold is divided into a stress transition area; For the stress concentration area, a stress amplification coefficient is calculated based on a difference between the stress concentration area discrimination value and the preset concentration threshold; for the stress transition area, a stress reduction coefficient is calculated based on a distance to the center point of the stress concentration area; The stress amplification coefficient is applied to the normal tension distribution of the stress concentration area, the stress reduction coefficient is applied to the normal tension distribution of the transition area, a modified normal tension distribution is obtained, a smooth transition operation is performed on the modified normal tension distribution, and integration is performed on the entire calculation domain to obtain an equivalent normal tension considering the stress distribution characteristics.

5. The method of claim 1, wherein, The equivalent shear strain and the equivalent normal tension are mapped to an elastoplastic constitutive space, and a damage softening value is obtained by constructing a stress-strain hyperbola, including: The equivalent normal tension is divided by a preset shear modulus to obtain an elastic component, the equivalent shear strain is subtracted by the elastic component to obtain a plastic component, a trace of a stress tensor of the equivalent normal tension is obtained to obtain a hydrostatic pressure component, and the equivalent normal tension is subtracted by the hydrostatic pressure component to obtain a deviatoric stress component; an elastoplastic constitutive mapping of the interface element is established based on the hydrostatic pressure component, the deviatoric stress component, the elastic component and the plastic component, and the elastoplastic constitutive mapping is taken as a constraint condition; The ratio of the equivalent shear strain to the initial stress is taken as an elastic modulus, and the ratio of the equivalent shear strain to a yield stress is taken as a hardening index. Under the constraint condition, the equivalent shear strain, the equivalent normal tension, and the hardening index are subjected to nonlinear constitutive mapping to obtain a stress-strain hyperbola. A modified elastic modulus is determined based on the stress-strain hyperbola, the modified elastic modulus is compared with an effective elastic modulus at a current time, and a current damage value is determined based on the difference between the two. The current damage value is multiplied by the modified elastic modulus and divided by the effective elastic modulus, and a damage evolution rate is obtained by subtracting the quotient from one. The damage evolution rate is integrated in the time domain to obtain a damage softening value.

6. The method of claim 1, wherein, Based on the damage softening value, a return mapping algorithm is used to adjust the plastic flow direction in the plastic correction process to obtain an energy dissipation increment. When the energy dissipation increment exceeds a critical threshold, an additional damage amount is calculated, including: The product of the damage softening value and the initial yield strength of the material is determined as the current yield strength, and a yield surface is determined based on the current yield strength. The distance between a trial stress state and the yield surface is calculated to obtain a plastic multiplier, and the product of the plastic multiplier and the yield surface normal vector is determined as a stress return direction. The stress is decomposed along the stress return direction using a return mapping algorithm, and the plastic strain increment is calculated by iteration. The product of the plastic strain increment and the current stress state is determined as the energy dissipation increment. When the energy dissipation increment exceeds the critical threshold, the damage driving force is obtained by subtracting one from the ratio of the energy dissipation increment to the critical threshold. The additional damage amount is obtained by multiplying the damage driving force by a material sensitivity coefficient.

7. A system for simulating interface contact failure based on the coupling of historical shear strain and normal tension, for implementing the method according to any one of claims 1 to 6, characterized in that, including: The first unit is configured to obtain a historical shear strain of an interface element between a foam core and upper and lower skins of a foam sandwich structure. The second unit is configured to scan the historical shear strain by a sliding time window, calculate a shear strain gradient and a shear strain peak value in each time window, determine a damage acceleration coefficient and a damage recovery coefficient according to the shear strain gradient and the shear strain peak value, and act on the historical shear strain to obtain an equivalent shear strain. The third unit is configured to construct a stress field topology relationship of an interface unit, extract stress field gradient features and stress field distribution features through a graph attention layer, obtain a normal tension distribution of the interface unit and adjacent units thereof, adaptively divide a stress concentration area and a stress transition area based on the stress field gradient features, the stress field distribution features and the normal tension distribution, and respectively apply a stress amplification coefficient and a stress reduction coefficient for regional correction to obtain an equivalent normal tension. The fourth unit is configured to determine that the interface contact fails when a ratio of the equivalent shear strain to a limit shear strain and a sum of a ratio of the equivalent normal tension to a maximum allowable normal tension of the interface are greater than or equal to one. The fifth unit is configured to map the equivalent shear strain and the equivalent normal tension to an elastic-plastic constitutive space, and obtain a damage softening value by constructing a stress-strain hyperbolic curve. Based on the damage softening value, a return mapping algorithm is used to adjust a plastic flow direction in a plastic correction process to obtain an energy dissipation increment. When the energy dissipation increment exceeds a critical threshold value, an additional damage amount is calculated, and the additional damage amount is superimposed to obtain an interface failure parameter.

8. An electronic device, comprising: The computer program instructions are executed by the processor to implement the method of any one of claims 1-6. The computer program instructions are executed by the processor to implement the method of any one of claims 1-6. The computer program instructions are executed by the processor to implement the method of any one of claims 1-6. ​ 9. A computer-readable storage medium having stored thereon computer program instructions, wherein, ​

Images