A matrix calculation method for multiaxial fatigue life based on finite element results

Abstract

The invention discloses a matrix calculation method for multi-axis fatigue life based on finite element results. Based on the stress / strain tensor corresponding to the global three-dimensional coordinate axis, the stress / strain tensor corresponding to the three-dimensional coordinate axis of the inclined plane is calculated by using coordinate transformation and matrix multiplication, and the normal stress / strain data and shear stress / strain data in the stress / strain tensor are classified; a stress / strain storage array, a spatial discrete coordinate transformation array and a spatial transformation stress / strain array are established; the normal strain data and the shear strain data in the spatial transformation strain array are respectively extracted, and the maximum normal strain amplitude and the maximum shear strain amplitude corresponding to all nodes are calculated; based on the maximum normal strain amplitude and the maximum shear strain amplitude, as well as the normal strain critical surface fatigue life calculation formula and the shear strain critical surface fatigue life calculation formula, the fatigue life corresponding to all nodes is solved. The invention adopts matrix operation to replace the loop, judgment and scalar algebra operation in the traditional algorithm.

Description

Technical Field

[0001] The invention relates to the technical field of fatigue life calculation, and in particular to a multi-axis fatigue life matrix calculation method based on finite element results. Background Art

[0002] In engineering practice, fatigue life assessment usually involves complex geometries under multiaxial stress states, which makes the calculation of fatigue life very complicated. A large number of experiments have shown that cracks always form and grow on certain specific material planes, so fatigue criteria based on critical planes have been widely used. The algorithms of these criteria all involve finding the direction of the critical plane of the material with maximum fatigue damage. However, the computational cost of the critical plane is very high, and existing algorithms need to traverse the entire orientation space based on the stress-strain history of each point.

[0003] For real engineering cases, the load history in the structure is usually calculated by finite element software. For large structures, the finite element model usually has hundreds of thousands or millions of nodes in space and hundreds or thousands of discrete points in time. For this situation, the existing critical plane fatigue assessment method may become too expensive, especially in the early design iteration stage, when a single fatigue life may take tens or even hundreds of hours. Summary of the invention

[0004] The purpose of the present invention is to provide a multi-axis fatigue life matrix calculation method based on finite element results to solve the technical problem that the critical plane fatigue assessment method in the prior art may become too expensive, especially in the early design iteration stage, when a single fatigue life may take dozens or even hundreds of hours.

[0005] In order to solve the above technical problems, the present invention specifically provides the following technical solutions:

[0006] A matrix calculation method for multi-axial fatigue life based on finite element results includes the following steps:

[0007] Step 100: Based on the stress / strain tensor corresponding to the global three-dimensional coordinate axis, the stress / strain tensor corresponding to the three-dimensional coordinate axis of the inclined plane is calculated by using coordinate transformation and matrix multiplication, and the normal stress / strain data and shear stress / strain data in the stress tensor are classified;

[0008] Step 200, setting the storage structure of the entire stress / strain data of the stress / strain tensor of the inclined surface, storing the stress / strain varying with time as a three-dimensional array, and establishing a stress / strain storage array;

[0009] Step 300, storing the spatial discrete matrix of the three-dimensional coordinate axes of the inclined plane as a three-dimensional array, and establishing a spatial discrete coordinate transformation array;

[0010] Step 400, performing matrix multiplication operation on the three-dimensional data of the spatial discrete coordinate transformation array and the three-dimensional data of the stress / strain storage array to obtain a spatial transformation stress array of all inclined surface stresses / strains after spatial discretization;

[0011] Step 500, obtaining a calculation formula for a normal strain critical surface and a calculation formula for a shear strain critical surface according to different critical surface criteria, and extracting normal stress / strain data and shear stress / strain data in the spatial transformation stress / strain array respectively, and calculating the maximum normal strain amplitude and the maximum shear strain amplitude corresponding to all nodes;

[0012] Step 600: Based on the maximum normal strain amplitude and the maximum shear strain amplitude, as well as the SWT fatigue life calculation formula and the Fatemi-Socie fatigue life calculation formula, the fatigue life corresponding to all nodes is solved.

[0013] As a preferred solution of the present invention, in step 100, the stress tensor corresponding to the three-dimensional coordinate axes of the inclined plane has symmetry, and the calculation formula for obtaining the inclined plane stress / strain using the Voigt method is:

[0014] {σ″}=[T]{σ}

[0015] {ε″}=[T]{ε}

[0016] Among them, σ is the stress tensor in the global coordinates xyz of any stress state, and ε is the strain tensor in the global coordinates xyz of any stress state;

[0017] [T] is a 6×6 coordinate transformation matrix;

[0018] σ″ is the stress tensor in the three-dimensional coordinates of the inclined plane, and {σ″} and {σ} are 6×1 column vectors using Voigt notation.

[0019] ε″ is the strain tensor in the three-dimensional coordinates of the inclined plane, and {ε″} and {ε} are 6×1 column vectors using Voigt notation.

[0020] As a preferred solution of the present invention, the Voigt method is used to obtain the inclined plane stress calculation method, the normal stress / strain and shear stress / strain on the inclined plane are calculated according to matrix multiplication, the normal stress and shear stress in {σ″} are classified, and the normal strain and shear strain in {ε″} are classified.

[0021] As a preferred solution of the present invention, in step 200, the first dimension of the stress-strain storage array represents the number of nodes, and the length is i;

[0022] The second dimension represents the discrete number of time, with a length of j;

[0023] The third dimension represents the stress / strain vector, with a length of 6;

[0024] The expression of the stress / strain storage array is (i, j, 6).

[0025] As a preferred solution of the present invention, in step 300, the spatial discrete coordinate transformation array discretizes the coordinate transformation matrix in three-dimensional space, and the first dimension of the spatial discrete coordinate transformation array represents the number of three-dimensional spatial discrete coordinate systems of any stress / strain state, and the length is k;

[0026] The second dimension and the third dimension represent the dimensions of the coordinate transformation matrix, and both have a length of 6;

[0027] The expression of the spatial discrete coordinate transformation array is (k, 6, 6).

[0028] As a preferred solution of the present invention, in step 400, the spatial transformation stress / strain array of all inclined surface stress / strain after spatial discretization is a four-dimensional array, and the size of the spatial transformation stress / strain array is (i, j, k, 6).

[0029] As a preferred solution of the present invention, in step 500, first, the maximum normal strain amplitude is used as the critical surface criterion, the spatial transformation strain array is calculated, the normal strain data is extracted from the spatial transformation strain array, the maximum normal stress amplitude of the inclined surface is calculated as the critical surface and the corresponding coordinate system is recorded. The specific implementation steps are:

[0030] Based on the normal strain classification in ε″, the normal stress part in the spatial transformation strain array (i, j, k, 6) is taken out to establish the normal strain subset (i, j, k, 1);

[0031] The maximum value minus the minimum value is calculated for the second dimension j of the subset to obtain the normal strain amplitude on the inclined plane. The maximum value is obtained in the third dimension k direction to obtain the maximum normal strain amplitude corresponding to all nodes and the inclined plane coordinate system corresponding to the critical surface.

[0032] As a preferred solution of the present invention, in step 600, first, the maximum shear strain amplitude is used as the critical surface criterion, the maximum shear strain amplitude critical surface is calculated in combination with the spatial transformation strain array, and the shear strain data is extracted from the spatial transformation strain array to calculate the inclined surface maximum shear strain critical surface and the corresponding coordinate system. The specific implementation steps are:

[0033] Based on the shear strain classification in ε″, the shear strain part in the spatial transformation strain array (i, j, k, 6) is taken out to establish the shear strain subset (i, j, k, 2);

[0034] Perform the minimum inscribed rectangle operation on the second dimension j and the fourth dimension of the subset to obtain the shear strain amplitude on the inclined plane, and find the maximum value in the third dimension k direction to obtain the maximum shear strain amplitude corresponding to all nodes and the inclined plane coordinate system corresponding to the critical surface.

[0035] As a preferred solution of the present invention, the calculated maximum shear strain amplitude and the maximum normal strain amplitude are respectively substituted into the calculation of fatigue life N f In the two formulas, the nonlinear equations are solved by Newton's method to obtain the fatigue life N corresponding to all nodes. f ; Among them, the calculated fatigue life N f They are:

[0036]

[0037] Where σ n,max represents the maximum normal stress on the critical surface, represents the normal strain amplitude, σ′ f is the fatigue strength coefficient, E is Young's modulus, N f is fatigue life, b is fatigue strength index, ε′ f is the fatigue ductility coefficient, and c is the fatigue ductility index.

[0038]

[0039] Where σ n,max represents the maximum normal stress on the critical surface, k is the material parameter, σ y is the yield stress of the material, represents the shear strain amplitude, τ′ f is the shear fatigue strength coefficient, G is the shear modulus, N f is fatigue life, b γ is the shear fatigue strength index, γ′ f is the shear fatigue ductility coefficient, c γ is the shear fatigue ductility index.

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

[0041] The present invention is applicable to the solution and calculation of various critical surface fatigue models, based on the most general three-dimensional complex load conditions, applicable to any geometric shape, applicable to any complex stress-strain input, and matrix multiplication is very suitable for parallel computing, with a rich and mature code base, and extremely high parallel efficiency. As long as the configuration of computer hardware is improved, the life assessment of complex structures can be achieved in a very short time. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] In order to more clearly illustrate the implementation methods of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for the implementation methods or the description of the prior art. Obviously, the drawings in the following description are only exemplary, and for ordinary technicians in this field, other implementation drawings can be derived from the provided drawings without creative work.

[0043] Figure 1 A three-dimensional coordinate diagram of an arbitrary stress state of a critical surface provided by an embodiment of the present invention;

[0044] Figure 2 A stress diagram of an inclined surface provided by an embodiment of the present invention;

[0045] Figure 3 A schematic diagram of calculating inclined surface stress provided by an embodiment of the present invention;

[0046] Figure 4 A schematic diagram of matrix calculation of normal stress provided by an embodiment of the present invention;

[0047] Figure 5 A schematic diagram of matrix calculation of shear force provided in an embodiment of the present invention;

[0048] Figure 6 A three-dimensional diagram of a stress / strain storage array provided by an embodiment of the present invention;

[0049] Figure 7 A three-dimensional graph of a spatial discrete coordinate transformation array provided by an embodiment of the present invention;

[0050] Figure 8 A three-dimensional graph of a spatially transformed stress array provided by an embodiment of the present invention;

[0051] Fig. 9 A flowchart of matrix calculation of multi-axis fatigue life provided by an embodiment of the present invention. DETAILED DESCRIPTION

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

[0053] like Fig. 9 As shown, the present invention provides a matrix calculation method for multi-axis fatigue life based on finite element results. The traditional critical surface search method needs to traverse the entire orientation space based on the stress-strain history of each point, such as Figure 1 and Figure 2As shown in the figure, the calculation process of the critical surface of the maximum normal stress of the inclined surface is as follows:

[0054] Given an arbitrary stress state:

[0055]

[0056] Then the normal stress σ on the inclined plane Δ is n (t) and shear stress τ n The calculation formula of (t) is

[0057] n x = sin(θ)·cos(φ)

[0058] n y = sin(θ)·sin(φ)

[0059] n z =cos(θ)

[0060] a x = -sin(φ)

[0061] a y =cos(φ)

[0062] a z =0

[0063] b x = -cos(θ)·cos(φ)

[0064] b y = -cos(θ)·sin(φ)

[0065] b z = sin(θ)

[0066] t x (t) = σ x (t)·n x +τ xy (t)·n y +τ xz (t)·n z

[0067] t y (t) = τ xy (t)·n x +σ y (t)·n y +τ yz (t)·n z

[0068] t z (t) = τ xz (t)·n x +τ yz(t)·n y +σ z (t)·n z

[0069] σ n (t) = t x (t)·n x +t y (t)·n y +t z (t)·n z

[0070] τ na (t) = t x (t)·a x +t y (t)·a y +t z (t)·a z

[0071] τ nb (t) = t x (t) b x +t y (t) b y +t z (t) b z

[0072]

[0073] At this point, if we want to retrieve the critical surface of the maximum normal stress, the algorithm is as follows:

[0074] ① Loop through the nodes in the finite element model;

[0075] ② Loop through the theta angle;

[0076] ③ Loop through the phi angle;

[0077] ④ Cycle through the load time t;

[0078] ⑤Calculate σ at the current angle and time n (t);

[0079] ⑥Calculate the current theta, Δσ under the phi surface n ;

[0080] ⑦The largest Δσ under all theta and phi surfaces n ;

[0081] Then the theta, phi surface is the critical surface of the current node. Assuming that the number of nodes in the finite element is n, the number of theta angles is i, the number of phi angles is j, and the number of load times is m, the algorithm contains a total of n*i*j*m cycles. The algorithm is inefficient and not conducive to parallel computing.

[0082] Based on the finite element calculation results, this implementation method designs a matrix calculation method for the critical surface fatigue model, which facilitates the parallel calculation of large structures and greatly improves the efficiency of structural life calculation.

[0083] This embodiment includes the following steps: Fig. 9 As shown:

[0084] Step 100: Based on the stress / strain tensor corresponding to the global three-dimensional coordinate axis, the stress / strain tensor corresponding to the three-dimensional coordinate axis of the inclined plane is calculated by using coordinate transformation and matrix multiplication, and the normal stress / strain data and shear stress / strain data in the stress tensor are classified.

[0085] Combined with Figure 1 and Figure 2 As shown, σ is the stress tensor under the global coordinate xyz, σ″ is the stress tensor under the inclined plane coordinate nba, the horizontal plane is the xy plane of the xyz coordinate system and the normal is z. The shadow plane is the ab plane of the nba coordinate system, and the normal is n. The global coordinate system is an artificially defined reference coordinate system. The calculation formula of the stress tensor under the traditional inclined plane coordinate nba is: σ″=MσM T ;

[0086] in,

[0087]

[0088] Right now

[0089] Then, the normal stress σ on the inclined plane is n (t) = σ″ 11 ;

[0090] Shear stress τ on the inclined plane na (t) = σ″; τ nb (t) = σ″ 13 ;

[0091] In step 100, the stress tensor corresponding to the three-dimensional coordinate axes of the inclined plane is symmetrical, and the calculation formula for obtaining the inclined plane stress / strain using the Voigt method is:

[0092] {σ″}=[T]{σ}

[0093] {ε″}=[T]{ε}

[0094] Among them, σ is the stress tensor in the global coordinates xyz of any stress state, and ε is the strain tensor in the global coordinates xyz of any stress state;

[0095] [T] is a 6×6 coordinate transformation matrix;

[0096] σ″ is the stress tensor in the three-dimensional coordinates of the inclined plane, and {σ″} and {σ} are 6×1 column vectors using Voigt notation.

[0097] ε″ is the strain tensor in the three-dimensional coordinates of the inclined plane, and {ε″} and {ε} are 6×1 column vectors using Voigt notation.

[0098] in,

[0099] The relationship between the [T] matrix and the M matrix is ​​as follows, "**" represents the exponentiation:

[0100] T11=M11**2; T12=M12**2; T13=M13**2; T14=2*M11*M12; T15=2*M11*M13; T16=2*M12*M13;

[0101] T21=M21**2; T22=M22**2; T21=M23**2; T24=2*M21*M22; T25=2*M21*M23; T26=2*M22*M23;

[0102] T31=M31**2; T32=M32**2; T33=M33**2; T34=2*M31*M32; T35=2*M31*M33; T36=2*M32*M33;

[0103] T41=M11*M21; T42=M12*M22; T43=M13*M23; T44=M12*M21+M11*M22; T45=M13*M21+M11*M23; T46=M13*M22+M12*M23;

[0104] T51=M11*M31; T52=M12*M32; T53=M13*M33; T54=M12*M31+M11*M32; T55=M13*M31+M11*M33; T56=M13*M32+M12*M33;

[0105] T61=M21*M31; T62=M22*M32; T63=M23*M33; T64=M22*M31+M21*M32; T65=M23*M31+M21*M33; T66=M23*M32+M22*M33.

[0106] The Voigt method is used to obtain the inclined surface stress calculation method, and the normal stress / strain and shear stress / strain on the inclined surface are calculated according to matrix multiplication. The normal stress and shear stress in {σ″'} are classified, and the normal strain and shear strain in {ε″} are classified. Figure 4 and Figure 5 As shown, based on the definition of normal stress and shear stress. Direction 1 is the normal direction of the plane, so it is normal stress, and directions 2 and 3 are tangent to the plane, so they are shear stress. Similarly, the normal stress σ on the inclined plane n (t) = σ″; shear stress τ on the inclined plane na (t) = σ″ 12 ; τ nb (t) = σ″ 13 .

[0107] Step 200 , setting the storage structure of the entire stress / strain data of the stress / strain tensor of the inclined surface, storing the stress / strain that changes with time as a three-dimensional array, and establishing a stress / strain storage array.

[0108] like Figure 6 As shown, the first dimension of the stress / strain storage array represents the number of nodes, and the length is i;

[0109] The second dimension represents the discrete number of time, with a length of j;

[0110] The third dimension represents the stress / strain vector, with a length of 6;

[0111] The expression of the stress / strain storage array is (i, j, 6).

[0112] Step 300, storing the spatial discrete matrix of the three-dimensional coordinate axes of the inclined plane as a three-dimensional array, and establishing a spatial discrete coordinate transformation array;

[0113] like Figure 7 As shown, in step 300, the spatial discrete coordinate transformation array discretizes the coordinate transformation matrix into three-dimensional space, and the first dimension of the spatial discrete coordinate transformation array represents the number of three-dimensional spatial discrete coordinate systems of any stress / strain state, and the length is k;

[0114] The second dimension and the third dimension represent the dimensions of the coordinate transformation matrix, and both have a length of 6;

[0115] The expression of the spatial discrete coordinate transformation array is (k, 6, 6). If the three-dimensional space is discretized at 180°×180°, the size of the matrix is ​​(32400, 6, 6).

[0116] Step 400: Perform matrix multiplication operation on the three-dimensional data of the spatial discrete coordinate transformation array and the three-dimensional data of the stress / strain storage array to obtain a spatial transformation stress array of all inclined surface stresses / strains after spatial discretization.

[0117] In step 400, the spatial transformation stress / strain array of all inclined surface stress / strain after spatial discretization is a four-dimensional array, and the size of the spatial transformation stress / strain array is (i, j, k, 6)

[0118] like Figure 8 As shown, the spatial transformation stress array of all inclined surface stresses after spatial discretization is a four-dimensional array, and the size of the spatial transformation stress array is (i, j, k, 6). This embodiment realizes all operations of steps ①②③④ in the above-mentioned traditional algorithm through one matrix multiplication, that is, matrix operations are used to replace all loops, judgments and scalar algebraic operations in the traditional algorithm.

[0119] Step 500: obtaining a calculation formula for the normal strain critical surface and a calculation formula for the shear strain critical surface according to different critical surface criteria, and extracting the normal stress / strain data and the shear stress / strain data in the spatial transformation stress / strain array respectively, and calculating the maximum normal strain amplitude and the maximum shear strain amplitude corresponding to all nodes.

[0120] Step 600: Based on the maximum normal strain amplitude and the maximum shear strain amplitude, as well as the SWT fatigue life calculation formula and the Fatemi-Socie fatigue life calculation formula, the fatigue life corresponding to all nodes is solved.

[0121] In step 500, first, the maximum normal strain amplitude is used as the critical surface criterion, the spatial transformation strain array is calculated, the normal strain data is extracted from the spatial transformation strain array, the maximum normal stress amplitude of the inclined surface is calculated as the critical surface and the corresponding coordinate system is recorded. The specific implementation steps are:

[0122] Based on the normal strain classification in ε″, the normal stress part in the spatial transformation strain array (i, j, k, 6) is taken out to establish the normal strain subset (i, j, k, 1);

[0123] The maximum value minus the minimum value is calculated for the second dimension j of the subset to obtain the normal strain amplitude on the inclined plane. The maximum value is obtained in the third dimension k direction to obtain the maximum normal strain amplitude corresponding to all nodes and the inclined plane coordinate system corresponding to the critical surface.

[0124] In step 600, first, the maximum shear strain amplitude is used as the critical surface criterion, and the maximum shear strain amplitude critical surface is calculated in combination with the spatial transformation strain array. The shear strain data is extracted from the spatial transformation strain array to calculate the inclined surface maximum shear strain critical surface and the corresponding coordinate system. The specific implementation steps are:

[0125] Based on the shear strain classification in ε″, the shear strain part in the spatial transformation strain array (i, j, k, 6) is taken out to establish the shear strain subset (i, j, k, 2);

[0126] Perform the minimum inscribed rectangle operation on the second dimension j and the fourth dimension of the subset to obtain the shear strain amplitude on the inclined plane, and find the maximum value in the third dimension k direction to obtain the maximum shear strain amplitude corresponding to all nodes and the inclined plane coordinate system corresponding to the critical surface.

[0127] Substitute the calculated maximum shear strain amplitude and maximum normal strain amplitude into the two formulas for calculating fatigue life Nf, and solve the nonlinear equations by Newton's method to obtain the fatigue life Nf corresponding to all nodes; among which, the calculated fatigue life Nf is:

[0128] Among them, the fatigue life N is calculated f They are:

[0129]

[0130] Where σ n,max represents the maximum normal stress on the critical surface, represents the normal strain amplitude, σ′ f is the fatigue strength coefficient, E is Young's modulus, N f is fatigue life, b is fatigue strength index, ε′ f is the fatigue ductility coefficient, and c is the fatigue ductility index.

[0131]

[0132] Where σ n,max represents the maximum normal stress on the critical surface, k is the material parameter, σ y is the yield stress of the material, represents the shear strain amplitude, τ′ f is the shear fatigue strength coefficient, G is the shear modulus, N f is fatigue life, b γ is the shear fatigue strength index, γ′ f is the shear fatigue ductility coefficient, c γ is the shear fatigue ductility index.

[0133] Substitute the maximum normal strain (i,1) and the maximum shear normal strain (i,1) into the above two formulas (1) and (2), respectively. The fatigue life Nf is the variable to be solved, and the rest are material parameters that can be obtained through experiments or manuals. By solving the nonlinear equations using the Newton method, the fatigue life Nf corresponding to all nodes can be obtained.

[0134] Therefore, this implementation is suitable for solving and calculating various critical surface models, based on the most general three-dimensional complex load conditions, applicable to any geometric shape, applicable to any complex stress and strain input, and matrix multiplication is very suitable for parallel computing, with a rich and mature code library and extremely high parallel efficiency. As long as the configuration of computer hardware is improved, the life assessment of complex structures can be achieved in a very short time.

[0135] The above embodiments are only exemplary embodiments of the present application and are not intended to limit the present application. The protection scope of the present application is defined by the claims. Those skilled in the art may make various modifications or equivalent substitutions to the present application within the essence and protection scope of the present application, and such modifications or equivalent substitutions shall also be deemed to fall within the protection scope of the present application.

Claims

1. A matrix calculation method for multi-axial fatigue life based on finite element results, It is characterized in that The following steps are involved: Step 100: Based on the stress / strain tensor corresponding to the global three-dimensional coordinate axis, the stress / strain tensor corresponding to the three-dimensional coordinate axis of the inclined plane is calculated by using coordinate transformation and matrix multiplication, and the normal stress / strain data and shear stress / strain data in the stress tensor are classified; In step 100, the stress tensor corresponding to the three-dimensional coordinate axes of the inclined plane is symmetrical, and the calculation formula for obtaining the inclined plane stress / strain using the Voigt method is: {σ″}=[T]{σ} {ε″}=[T]{ε} Where σ is the stress tensor in the global coordinates xyz of any stress state, and ε is the strain tensor in the global coordinates xyz of any stress state; [T] is a 6×6 coordinate transformation matrix; σ″ is the stress tensor in the three-dimensional coordinates of the inclined plane, and {σ″} and {σ} are 6×1 column vectors using Voigt notation; ε″ is the strain tensor in the three-dimensional coordinates of the inclined plane, and {ε″} and {ε} are 6×1 column vectors using Voigt notation; Step 200, setting the storage structure of the entire stress / strain data of the stress / strain tensor of the inclined surface, storing the stress / strain varying with time as a three-dimensional array, and establishing a stress / strain storage array; Step 300, storing the spatial discrete matrix of the three-dimensional coordinate axes of the inclined plane as a three-dimensional array, and establishing a spatial discrete coordinate transformation array; Step 400, performing matrix multiplication operation on the three-dimensional data of the spatial discrete coordinate transformation array and the three-dimensional data of the stress / strain storage array to obtain a spatial transformation stress array of all inclined surface stresses / strains after spatial discretization; Step 500, obtaining a calculation formula for a normal strain critical surface and a calculation formula for a shear strain critical surface according to different critical surface criteria, and extracting normal stress / strain data and shear stress / strain data in the spatial transformation stress / strain array respectively, and calculating the maximum normal strain amplitude and the maximum shear strain amplitude corresponding to all nodes; Step 600, based on the maximum normal strain amplitude and the maximum shear strain amplitude, as well as the SWT fatigue life calculation formula and the Fatemi-Socie fatigue life calculation formula, solve the fatigue life corresponding to all nodes; Substitute the calculated maximum shear strain amplitude and maximum normal strain amplitude into the fatigue life N f In the two formulas, the nonlinear equations are solved by Newton's method to obtain the fatigue life N corresponding to all nodes. f ; Among them, the calculated fatigue life N f They are: Where σ n,max represents the maximum normal stress on the critical surface, represents the normal strain amplitude, σ′ f is the fatigue strength coefficient, E is Young's modulus, N f is fatigue life, b is fatigue strength index, ε f ′ is the fatigue ductility coefficient, c is the fatigue ductility index; Where σ n,max represents the maximum normal stress on the critical surface, k is the material parameter, σ y is the yield stress of the material, represents the shear strain amplitude, τ′ f is the shear fatigue strength coefficient, G is the shear modulus, N f is fatigue life, b γ is the shear fatigue strength index, γ f ′ is the shear fatigue ductility coefficient, c γ is the shear fatigue ductility index.

2. According to the multi-axis fatigue life matrix calculation method based on finite element results according to claim 1, Features: The Voigt method is used to obtain the inclined plane stress calculation method, the normal stress / strain and shear stress / strain on the inclined plane are calculated according to matrix multiplication, the normal stress and shear stress in {σ″} are classified, and the normal strain and shear strain in {ε″} are classified.

3. The matrix calculation method of multi-axial fatigue life based on finite element results according to claim 1, Features: In step 200, the first dimension of the stress / strain storage array represents the number of nodes, and the length is i; The second dimension represents the discrete number of time, with a length of j; The third dimension represents the stress / strain vector, with a length of 6; The expression of the stress / strain storage array is (i, j, 6).

4. The matrix calculation method of multi-axial fatigue life based on finite element results according to claim 2, Features: In step 300, the spatial discrete coordinate transformation array discretizes the coordinate transformation matrix into three-dimensional space, and the first dimension of the spatial discrete coordinate transformation array represents the number of three-dimensional spatial discrete coordinate systems of any stress / strain state, and the length is k; The second dimension and the third dimension represent the dimensions of the coordinate transformation matrix, and both have a length of 6; The expression of the spatial discrete coordinate transformation array is (k, 6, 6).

5. The matrix calculation method for multi-axial fatigue life based on finite element results according to claim 4, Features: In step 400, the spatially transformed stress / strain array of all inclined surface stresses / strains after spatial discretization is a four-dimensional array, and the size of the spatially transformed stress / strain array is (i, j, k, 6).

6. The matrix calculation method for multi-axial fatigue life based on finite element results according to claim 5, Features: In step 500, first, the maximum normal strain amplitude is used as the critical surface criterion, the spatial transformation strain array is calculated, the normal strain data is extracted from the spatial transformation strain array, the maximum normal stress amplitude of the inclined surface is calculated as the critical surface and the corresponding coordinate system is recorded. The specific implementation steps are: Based on the normal strain classification in ε″, the normal stress part in the spatial transformation strain array (i, j, k, 6) is taken out to establish the normal strain subset (i, j, k, 1); The maximum value minus the minimum value is calculated for the second dimension j of the subset to obtain the normal strain amplitude on the inclined plane. The maximum value is obtained in the third dimension k direction to obtain the maximum normal strain amplitude corresponding to all nodes and the inclined plane coordinate system corresponding to the critical surface.

7. The matrix calculation method for multi-axial fatigue life based on finite element results according to claim 6, Features: In step 600, first, the maximum shear strain amplitude is used as the critical surface criterion, and the maximum shear strain amplitude critical surface is calculated in combination with the spatial transformation strain array. The shear strain data is extracted from the spatial transformation strain array to calculate the inclined surface maximum shear strain critical surface and the corresponding coordinate system. The specific implementation steps are: Based on the shear strain classification in ε″, the shear strain part in the spatial transformation strain array (i, j, k, 6) is taken out to establish the shear strain subset (i, j, k, 2); Perform the minimum inscribed rectangle operation on the second dimension j and the fourth dimension of the subset to obtain the shear strain amplitude on the inclined plane, and find the maximum value in the third dimension k direction to obtain the maximum shear strain amplitude corresponding to all nodes and the inclined plane coordinate system corresponding to the critical surface.

Images