Simulation method for interface contact failure based on coupling of historical shear strain and normal tension

Abstract

The application provides an interface contact failure simulation method based on historical shear strain and normal tension coupling, relates to the technical field of structural strength simulation, and comprises the following steps: obtaining historical shear strain and calculating equivalent shear strain through a sliding time window; constructing a stress field topology relationship of an interface unit, extracting a stress field feature, obtaining a normal tension distribution, and obtaining equivalent normal tension by adaptively dividing a stress area; determining interface contact failure; and calculating interface failure parameters by constructing a stress-strain hyperbola to obtain a damage softening value. The application improves the accuracy and calculation efficiency of interface contact failure simulation of a foam sandwich structure.

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 area and transition area, and ignore the nonlinear coupling effect of normal tension and shear strain, which makes the failure prediction accuracy insufficient 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:

[0006] Obtaining the historical shear strain of the interface element between the foam core and the upper and lower skins of the foam sandwich structure;

[0007] Scanning the historical shear strain through a sliding time window, calculating the shear strain gradient and shear strain peak value in each time window, determining the damage acceleration coefficient and 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;

[0008] 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, normal tension distribution of the interface element and adjacent elements is obtained, the stress concentration area and the stress transition area are adaptively divided based on the stress field gradient features, the stress field distribution features and the normal tension distribution, and stress amplification coefficients and stress reduction coefficients are respectively applied for regional correction to obtain equivalent normal tension;

[0009] When the ratio of the equivalent shear strain to the limit shear strain and the ratio of the equivalent normal tension to the maximum normal tension allowed by the interface are greater than or equal to one, it is determined that the interface contact fails.

[0010] 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 hyperbolic curve; 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, and when the energy dissipation increment exceeds a critical threshold, an additional damage amount is calculated, and the interface failure parameter is obtained by superimposing the additional damage amount.

[0011] The historical shear strain is scanned through a sliding time window, the shear strain gradient and the shear strain peak value in each time window are calculated, and the damage acceleration coefficient and the damage recovery coefficient are determined according to the shear strain gradient and the shear strain peak value and are applied to the historical shear strain to obtain the equivalent shear strain, including:

[0012] The derivative of the shear strain with respect to time in the sliding time window is calculated, the shear strain gradient is obtained according to the maximum absolute value of the derivative, and the shear strain peak value is obtained according to the maximum value of the shear strain in the sliding time window;

[0013] The change amplitude and the average value of the historical shear strain are calculated to obtain the amplification coefficient and the adjustment coefficient, the shear strain gradient is amplified through the amplification coefficient, the amplified shear strain gradient is monotonically increasing mapped, the monotonically increasing mapped result is multiplied by the adjustment coefficient and added by one to obtain the damage acceleration coefficient, the change rate and the peak-valley difference of the historical shear strain are calculated to obtain the sensitivity coefficient and the strength coefficient, the difference between the shear strain peak value and the current shear strain is calculated, the difference is multiplied by the sensitivity coefficient and taken negatively, and monotonically decreasing mapping is performed, and the product of one minus the monotonically decreasing mapped result and the strength coefficient is obtained to obtain the damage recovery coefficient;

[0014] The time interval between adjacent extreme points of the historical shear strain is counted and taken as a time decay coefficient; the time difference between the current time and the historical time is calculated, the time difference is multiplied by the time decay coefficient and taken negatively, and monotonically decreasing mapping is performed to obtain a time weight, and the time weight is applied to the historical shear strain to obtain a time weighted integral value;

[0015] The damage acceleration coefficient, the damage recovery coefficient and the time weighted integral value are multiplied to obtain an equivalent shear strain considering strain history characteristics.

[0016] The stress field topology relationship of the interface element is constructed, the stress field gradient features and the stress field distribution features are extracted through a graph attention layer, and the normal tension distribution of the interface element and its adjacent elements is obtained.

[0017] The topological connection between the interface element and the adjacent elements is determined based on a preset distance threshold, the spatial distance between the interface element and the adjacent elements is calculated, then divided by a preset attenuation parameter and taken as a negative exponent to obtain a distance attenuation coefficient, and the distance attenuation coefficient is applied to the topological connection to obtain a weighted adjacency matrix.

[0018] The node feature vector corresponding to the interface element is linearly transformed 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 the adjacent elements, and the attention coefficient is applied to the transformed feature vector of the adjacent elements to obtain a plurality of single-head attention features, and the plurality of single-head attention features are fused to obtain the stress field gradient features of the interface element.

[0019] 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 components and y components of normal vectors of the interface element 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 the normal tension distribution.

[0020] 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:

[0021] The stress field gradient features are subjected to eigenvalue decomposition 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.

[0022] perform regional division on the interface elements based on the stress concentration area discrimination value, divide a region where the stress concentration area discrimination value is greater than a preset concentration threshold value into a stress concentration area; calculate a distance from other elements to a center point of the stress concentration area, divide a region where the distance divided by a preset transition width parameter and taking a negative exponent to obtain an attenuation value is greater than a preset transition threshold value into a stress transition area;

[0023] for the stress concentration area, calculate a stress amplification coefficient based on a difference between the stress concentration area discrimination value and the preset concentration threshold value; for the stress transition area, calculate a stress reduction coefficient based on the distance to the center point of the stress concentration area;

[0024] apply the stress amplification coefficient to a normal tension distribution of the stress concentration area, apply the stress reduction coefficient to a normal tension distribution of the transition area, to obtain a modified normal tension distribution, perform a smooth transition operation on the modified normal tension distribution and integrate on the entire calculation domain to obtain an equivalent normal tension considering stress distribution characteristics.

[0025] map the equivalent shear strain and the equivalent normal tension to an elastoplastic constitutive space, and obtain a damage softening value by constructing a stress-strain hyperbola, including:

[0026] divide the equivalent normal tension by a preset shear modulus to obtain an elastic component, subtract the elastic component from the equivalent shear strain to obtain a plastic component, take a trace of a stress tensor of the equivalent normal tension to obtain a hydrostatic pressure component, subtract the hydrostatic pressure component from the equivalent normal tension to obtain a deviatoric stress component; establish an elastoplastic constitutive mapping of the interface element based on the hydrostatic pressure component, the deviatoric stress component, the elastic component and the plastic component, and take the elastoplastic constitutive mapping as a constraint condition;

[0027] take a ratio of the equivalent shear strain to an initial stress as an elastic modulus, take a ratio of the equivalent shear strain to a yield stress as a hardening index, under the constraint condition, perform a nonlinear constitutive mapping on the equivalent shear strain, the equivalent normal tension, the hardening index to obtain a stress-strain hyperbola; determine a modified elastic modulus based on the stress-strain hyperbola, compare the modified elastic modulus with an effective elastic modulus at a current time, and determine a current damage value based on a difference between the two;

[0028] multiply the current damage value by the modified elastic modulus and divide by the effective elastic modulus, subtract one from the quotient to obtain a damage evolution rate, and integrate the damage evolution rate in the time domain to obtain a damage softening value.

[0029] Based on the damage softening value, an energy dissipation increment is obtained by adjusting the plastic flow direction in the plastic correction process using a return mapping algorithm, and when the energy dissipation increment exceeds a critical threshold, an additional damage amount is calculated, comprising:

[0030] A 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;

[0031] A distance between a trial stress state and the yield surface is calculated to obtain a plastic multiplier, a 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 iteratively calculated to obtain an energy dissipation increment;

[0032] When the energy dissipation increment exceeds the critical threshold, a damage driving force is obtained by subtracting one from the ratio of the energy dissipation increment to the critical threshold, and an additional damage amount is obtained by multiplying the damage driving force by a material sensitivity coefficient.

[0033] The second aspect of the embodiment of the application provides an interface contact failure simulation system based on historical shear strain and normal tension coupling, comprising:

[0034] A first unit is configured to obtain historical shear strain of an interface element between a foam core and upper and lower skins of a foam sandwich structure;

[0035] 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;

[0036] A third unit is configured to 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, 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;

[0037] 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 sum of a ratio of the equivalent normal tension to a maximum allowable normal tension of the interface, is greater than or equal to one.

[0038] 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 plastic flow direction in a plastic correction process is adjusted by using a return mapping algorithm to obtain an energy dissipation increment, and when the energy dissipation increment exceeds a critical threshold value, an additional damage amount is calculated, and the interface failure parameter is obtained by superimposing the additional damage amount.

[0039] A third aspect of the embodiments of the present application,

[0040] An electronic device is provided, comprising:

[0041] a processor;

[0042] a memory for storing processor-executable instructions;

[0043] The processor is configured to invoke the instructions stored in the memory to execute the method described above.

[0044] A fourth aspect of the embodiments of the present application,

[0045] 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.

[0046] The beneficial effects of the present application are as follows:

[0047] The present application constructs a brand-new interface contact failure simulation method, which can effectively capture the cumulative damage process of the interface contact area under cyclic loading by introducing a time window scanning mechanism to process the historical shear strain, and combining a damage acceleration coefficient and a recovery coefficient, thereby improving the accuracy and reliability of the simulation results.

[0048] The present application uses a graph attention layer to extract stress field features, adaptively divides stress concentration zones and transition zones and applies differentiated correction coefficients, effectively solving the limitations of traditional methods in dealing with complex stress distribution, enabling the simulation model to more accurately reflect the failure mechanism of the interface under non-uniform stress state, and greatly improving the simulation accuracy.

[0049] The present application maps the equivalent shear strain and the normal tension to the elastic-plastic constitutive space, solves the damage softening value by the hyperbolic curve, and adjusts the plastic flow direction by introducing the return mapping algorithm, which can accurately calculate the additional damage amount in the energy dissipation process, thereby comprehensively simulating the interface failure behavior at the macro and micro levels, and providing a reliable theoretical basis for the strength analysis and optimization design of the foam sandwich sandwich structure. BRIEF DESCRIPTION OF DRAWINGS

[0050] Figure 1A flowchart of an interface contact failure simulation method based on historical shear strain and normal tension coupling of an embodiment of the present application is shown in FIG.

[0051] Figure 2 A flowchart of a stress field feature extraction process of a graph attention network is shown in FIG. DETAILED DESCRIPTION

[0052] To make the objectives, technical solutions, and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described below in connection with the drawings of the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative work fall within the scope of protection of the present application.

[0053] The technical solutions of the present application will be described in detail below with specific embodiments. The following specific embodiments can be combined with each other, and some embodiments may not be described again for the same or similar concepts or processes.

[0054] Figure 1 A flowchart of an interface contact failure simulation method based on historical shear strain and normal tension coupling of an embodiment of the present application is shown in FIG. Figure 1 The method comprises:

[0055] obtaining historical shear strain of an interface element between a foam core and upper and lower skins of a foam sandwich structure;

[0056] 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, 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 to obtain an equivalent shear strain;

[0057] 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, obtaining normal tension distribution of the interface element and its adjacent elements, adaptively dividing 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 applying a stress amplification coefficient and a stress reduction coefficient for regional correction to obtain an equivalent normal tension;

[0058] 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, is greater than or equal to one, determining that the interface contact failure occurs;

[0059] 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, calculating an additional damage amount when the energy dissipation increment exceeds a critical threshold, and superimposing the additional damage amount to obtain an interface failure parameter.

[0060] In an optional embodiment, historical shear strain is scanned by a sliding time window, a shear strain gradient and a shear strain peak value in each time window are calculated, a damage acceleration coefficient and a damage recovery coefficient are determined according to the shear strain gradient and the shear strain peak value, and the historical shear strain is affected to obtain an equivalent shear strain, including:

[0061] derivative of the shear strain with respect to time in the sliding time window is calculated, a shear strain gradient is obtained according to a maximum absolute value of the derivative, and a shear strain peak value is obtained according to a maximum value of the shear strain in the sliding time window;

[0062] a change amplitude and an average value of the historical shear strain are calculated respectively to obtain an amplification coefficient and an adjustment coefficient, the shear strain gradient is amplified by the amplification coefficient, the amplified shear strain gradient is monotonically increasing mapped, a product of the monotonically increasing mapped result and the adjustment coefficient is added by one to obtain the damage acceleration coefficient, a change rate and a peak-valley difference value of the historical shear strain are calculated respectively to obtain a sensitivity coefficient and a strength coefficient, a difference value between the shear strain peak value and a shear strain at a current time is calculated, a product of the difference value and the sensitivity coefficient is taken negatively and monotonically decreasing mapped, and a product of one minus the monotonically decreasing mapped result and the strength coefficient is obtained to obtain the damage recovery coefficient;

[0063] a time interval between adjacent extreme points of the historical shear strain is counted and taken as a time decay coefficient; a time difference between a current time and a historical time is calculated, a product of the time difference and the time decay coefficient is taken negatively and monotonically decreasing mapped to obtain a time weight, and the time weight is applied to the historical shear strain to obtain a time weighted integral value;

[0064] the damage acceleration coefficient, the damage recovery coefficient and the time weighted integral value are multiplied to obtain an equivalent shear strain considering strain history characteristics.

[0065] In practical application, taking a specific shear strain history data as an example, suppose 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 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) is calculated between each adjacent two points, 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 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.

[0066] 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 and minimum values 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 to magnify the shear strain gradient. The average value is obtained by summing all historical shear strain values and dividing by the number of data points, if the average value of the historical data is 0.09, then it is taken as the adjustment coefficient. Multiply the shear strain gradient 0.03 / s by the magnification coefficient 0.16 to get 0.0048 / s, then perform a monotonic increasing mapping on this value. Here, an exponential mapping function is used to map 0.0048 / s to 1.005. Multiply the mapping result by the adjustment coefficient 0.09 and add 1 to get the damage acceleration coefficient as 1.09045.

[0067] The calculation of the damage recovery coefficient requires the change rate of shear strain and the peak-to-valley difference. The change rate can be calculated by the average value of the shear strain change between adjacent time points, for example, 0.025 / s is obtained as the 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 the strength coefficient. Calculate the difference between the shear strain peak value 0.14 and the current shear strain 0.11 as 0.03, multiply this difference by the sensitivity coefficient 0.025 / s to get 0.00075. Take the negative of this value to get -0.00075, then perform a monotonic decreasing mapping, using an exponential mapping function to get 0.99925. Subtract this value from 1 and multiply by the strength coefficient 0.16 to get 0.159880, and the damage recovery coefficient is 0.840120.

[0068] The calculation of the time weight first needs the time decay coefficient. The time interval between adjacent extreme points in the statistical history shear strain is calculated, for example, the time from the wave crest to the wave trough is 10 s, and the time from the wave trough to the wave crest is 15 s, so the average time interval is 12.5 s, and the reciprocal thereof is 0.08 / s. For the time difference Δt=tn-ti between the current time tn and the historical time ti, for example, Δt=20 s, 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 the exponential mapping function to obtain the time weight 0.2019. Multiply the time weight by the shear strain value γi=0.09 corresponding to the historical time to obtain the time-weighted integral value 0.01817.

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

[0070] In actual engineering applications, the sliding time window size, mapping function type and other parameters 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 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.

[0071] In an optional implementation, a stress field topology relationship of an interface element 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 thereof is obtained, including:

[0072] A topological connection between the interface element and the adjacent elements is determined based on a preset distance threshold, a distance decay coefficient is obtained by dividing a spatial distance between the interface element and the adjacent elements by a preset decay parameter and taking a negative exponential after the spatial distance is calculated, and the distance decay coefficient is applied to the topological connection to obtain a weighted adjacency matrix;

[0073] The node feature vector corresponding to the interface unit is linearly transformed 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 by LeakyReLU and normalized to obtain an attention coefficient between the interface unit and an adjacent unit, and the attention coefficient is applied to the transformed feature vector of the adjacent unit to obtain a plurality of single-head attention features, and the plurality of single-head attention features are fused to obtain a stress field gradient feature of the interface unit.

[0074] 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, a normal stress component and a tangential stress component are extracted, multiplied by x and y components of a normal vector of the interface unit to obtain x and y direction normal tension components, and the x and y direction normal tension components are superimposed to obtain a normal tension distribution.

[0075] As shown in Figure 2 The method comprises:

[0076] 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 the negative exponential 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.

[0077] Then, the node feature vector corresponding to the interface unit is linearly transformed. Each interface unit has an initial feature vector containing stress component, displacement component and other information, and the initial feature vector has a dimension of 64. The initial feature vector is converted into a transformed feature vector with a dimension of 32 by a linear transformation matrix. For example, for the interface unit 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.

[0078] The 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 adopted. 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, which have 32-dimensional transformed feature vectors [d1, d2,..., d32] and [e1, e2,..., e32] respectively, they 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 the 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, that is, 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.

[0079] The stress field gradient feature is input into the global attention pooling layer to extract the stress field distribution feature. The global attention pooling layer gives 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, the stress field gradient feature of each element has a dimension of 256, and the attention weight of each element is calculated by the global attention pooling layer, such as the weight of element H is 0.015 and the weight of element I 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.

[0080] The stress field distribution feature is decoded and reconstructed by 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 respectively. After the input stress field distribution feature is processed by the decoder, the stress components of the interface element are obtained, including the normal stress components σxx, σyy, σzz and the shear stress components τxy, τyz, τzx.

[0081] The normal stress component and the tangential stress component 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 tangential stress components are τxy=3 MPa, τyz=2 MPa, and τzx=1 MPa. Multiplying the normal stress component by the normal vector component 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. Repeat this calculation process for all interface elements to obtain the complete normal tension distribution.

[0082] By 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, which provides strong support for subsequent structure analysis and optimization.

[0083] In an optional implementation, based on the stress field gradient characteristics, the stress field distribution characteristics, and the normal tension distribution, the stress concentration area and the stress transition area are adaptively divided, and a stress amplification coefficient and a stress reduction coefficient are respectively applied for regional correction to obtain an equivalent normal tension, including:

[0084] The stress gradient strength index is obtained by performing eigenvalue decomposition on the stress field gradient characteristics 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 performing weighted combination with the ratio of the stress field distribution characteristics to the maximum stress value, and the stress concentration area discriminant value is obtained by normalizing the stress gradient strength index and performing weighted combination with the ratio of the stress field distribution characteristics to the maximum stress value.

[0085] Based on the stress concentration area discriminant value, the interface element is regionally divided, the region with a stress concentration area discriminant value 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 index is taken to obtain a decay value, and the region with a decay value greater than a preset transition threshold is divided into a stress transition area.

[0086] For the stress concentration area, a stress amplification coefficient is calculated based on the difference between the stress concentration area discriminant value and the preset concentration threshold, and for the stress transition area, a stress reduction coefficient is calculated based on the distance to the center point of the stress concentration area.

[0087] The stress amplification coefficient is applied to the normal tension distribution of the stress concentration area, and the stress reduction coefficient 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.

[0088] 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 eigenvalues of the stress field gradient characteristics of a certain region 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 region is 0.628.

[0089] The ratio of the stress field distribution characteristics 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 region 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 region is 0.628*0.6+0.833*0.4=0.71.

[0090] 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 region 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 region is divided into a stress concentration area.

[0091] 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.

[0092] 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.

[0093] 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.

[0094] 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.

[0095] 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.

[0096] 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.

[0097] 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.

[0098] 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, including:

[0099] 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 stress tensor of the equivalent normal tension is traced 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;

[0100] 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, 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.

[0101] 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 a quotient from 1. The damage evolution rate is integrated in a time domain to obtain a damage softening value.

[0102] In the interface element analysis process, first, an equivalent shear strain and an 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 the two physical quantities to an elastoplastic constitutive space is the basis for calculating the damage softening value.

[0103] The elastic component is obtained by dividing the equivalent normal tension by the preset shear modulus. 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 plastic component is obtained by subtracting the elastic component from the equivalent shear strain. Assuming that the equivalent shear strain is 0.005, the plastic component is 0.003. The hydrostatic pressure component is obtained by tracing the stress tensor of the equivalent normal tension. In three-dimensional space, it can be obtained by summing the diagonal elements and then dividing 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 deviatoric stress component is obtained by subtracting the hydrostatic pressure component from the equivalent normal tension, which is 10 MPa in this case.

[0104] Based on the hydrostatic pressure component, the deviatoric stress component, the elastic component, and the plastic component, an elastoplastic constitutive mapping of the interface element is established, and the mapping is used as a constraint condition. This mapping can be represented as a functional relationship between the hydrostatic pressure component and the deviatoric stress component, while considering the coupling effect of the elastic component and the plastic component. In practical applications, the Drucker-Prager model or the Mohr-Coulomb model can be used for description, which is represented as a calculation boundary condition.

[0105] 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 an iterative method, which includes the steps of initializing the stress state, calculating the current strain increment, updating the stress state, and checking the convergence.

[0106] In actual calculations, 10 discrete points can be selected for stress-strain mapping, such as strain points of 0.001, 0.002, 0.003, and 0.01, corresponding to calculated stress values of 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 through these discrete points.

[0107] The modified elastic modulus is determined based on the stress-strain hyperbolic curve. Specifically, the slope of the stress-strain curve in the region with small strain values (e.g., within the range of 0-0.002) is taken to calculate the modified elastic modulus as 19000 MPa. The modified 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 of 20000 MPa is 0.05.

[0108] The damage evolution rate is obtained by multiplying the current damage value by the modified elastic modulus, dividing by the effective elastic modulus, and subtracting the quotient from 1. Specifically, the calculation is as follows: 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 that 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.

[0109] 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 failure process of the interface element and provide a basis for structural safety evaluation.

[0110] In the above technical implementation process, the parameters such as the preset shear modulus, initial stress, and yield stress 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 material compressive strength, and the yield stress is about 40%-60% of the material compressive strength. The calculation accuracy is related to the selection of time step, and it is generally recommended that the time step should not exceed 0.05 seconds to ensure the stability and accuracy of numerical integration.

[0111] In an alternative embodiment, based on the damage softening value, a return mapping algorithm is used to adjust the plastic flow direction in the plastic modification process to obtain an energy dissipation increment, and when the energy dissipation increment exceeds a critical threshold, the additional damage amount is calculated.

[0112] 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.

[0113] A plastic multiplier is calculated by computing the distance between the trial stress state and the yield surface. A stress return direction is determined by multiplying the plastic multiplier with the yield surface normal vector. A stress return mapping algorithm is used to decompose the stress along the stress return direction. A plastic strain increment is iteratively calculated. An energy dissipation increment is determined by multiplying the plastic strain increment with the current stress state.

[0114] When the energy dissipation increment exceeds the critical threshold, a damage driving force is determined by subtracting one from the ratio of the energy dissipation increment and the critical threshold. An additional damage amount is determined by multiplying the damage driving force with a material sensitivity coefficient.

[0115] A current yield strength is determined by multiplying a damage softening value with an initial yield strength of the material. 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 Von Mises yield criterion, the yield surface can be represented as a spatial surface with an equivalent stress equal to the current yield strength.

[0116] A plastic multiplier is calculated by computing the distance between the trial stress state and the yield surface. In the process of elastoplastic calculation, the material is first assumed to be 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 trial stress obtained by elastic prediction is [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.

[0117] A stress return direction is determined by multiplying the plastic multiplier with the yield surface normal vector. The yield surface normal vector represents the direction of stress return. For 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). By multiplying the plastic multiplier 0.025 with the normal vector, the stress return direction is obtained as [0.02, 0.01, 0.005, 0.00175, 0.00225, 0.0015].

[0118] The stress is decomposed along the stress return direction using the return mapping algorithm, and the plastic strain increment is calculated iteratively. 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. Suppose the calculated plastic strain increment is [0.0001, 0.00005, 0.000025, 0.00000875, 0.0000105, 0.0000075].

[0119] 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 corrected stress state and summing them up, the energy dissipation increment is 0.048 J / mm 3 .

[0120] 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. Suppose 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.2x0.2=0.04.

[0121] 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 calculation, and then affect the subsequent yield strength and plastic deformation behavior.

[0122] 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 3Similarly, the material sensitivity coefficient can also be adjusted according to actual conditions. For materials sensitive to damage, a higher sensitivity coefficient, such as 0.3, can be set. For materials not sensitive to damage, a lower sensitivity coefficient, such as 0.1, can be set.

[0123] Through the above detailed steps, the process of adjusting the plastic flow direction in the plastic correction 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 value is realized. This method can effectively simulate the damage evolution behavior of materials during plastic deformation, providing reliable calculation basis for structural design and safety evaluation.

[0124] The embodiment of the application is based on an interface contact failure simulation system coupling historical shear strain and normal tension, and the system comprises:

[0125] The first unit is used to obtain the historical shear strain of the interface element between the foam core and the upper and lower skins of the foam sandwich structure;

[0126] The second unit is used to scan the historical shear strain through a sliding time window, calculate the shear strain gradient and shear strain peak value in each time window, determine the damage acceleration coefficient and damage recovery coefficient according to the shear strain gradient and the shear strain peak value, and act on the historical shear strain to obtain the equivalent shear strain;

[0127] The third unit is used to construct the stress field topology relationship of the interface element, extract the stress field gradient feature and stress field distribution feature through a graph attention layer, obtain the normal tension distribution of the interface element and its adjacent elements, and adaptively divide 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 apply the stress amplification coefficient and the stress reduction coefficient for regional correction to obtain the equivalent normal tension.

[0128] The fourth unit is used to determine the interface contact failure 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.

[0129] The fifth unit is used to map the equivalent shear strain and the equivalent normal tension to the elastic-plastic constitutive space, and obtain the damage softening value by constructing a stress-strain hyperbolic curve. Based on the damage softening value, the plastic flow direction in the plastic correction process is adjusted using the return mapping algorithm to obtain the energy dissipation increment. When the energy dissipation increment exceeds the critical threshold value, the additional damage amount is calculated, and the interface failure parameter is obtained by superimposing the additional damage amount.

[0130] In a third aspect, the embodiment of the application provides an electronic device, comprising:

[0131] a processor;

[0132] a memory for storing processor-executable instructions;

[0133] The processor is configured to invoke the instructions stored in the memory to perform the method described above.

[0134] In a fourth aspect, the present application provides a computer readable storage medium having stored thereon computer program instructions, which when executed by a processor, implement the method described above.

[0135] The present application can be a method, apparatus, system, and / or computer program product. Computer program products can include computer-readable storage media having computer-readable program instructions loaded thereon for performing various aspects of the present application.

[0136] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part or all of the technical features; 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