Fatigue crack propagation simulation and life prediction method considering crack closing effect

Abstract

The invention relates to a fatigue crack propagation simulation and life prediction method considering a crack closing effect. The method comprises the following steps: constructing a finite element model; assembling crack sizes in the finite element model according to user requirements, refining regional grids of the cracks, expanding the cracks in the finite element model, extracting a displacement field of the regional grids at the tips of the cracks, performing nonlinear calculation considering a crack closing effect on the basis of the displacement field, and calculating the crack closing effect. Obtaining an I-type stress intensity factor and an II-type stress intensity factor at the tip of the crack; calculating a crack propagation angle based on the I type stress intensity factor and the II type stress intensity factor and based on a maximum circumferential stress criterion, and calculating a new crack tip coordinate based on the crack tip coordinate, the crack propagation angle and a user-defined crack propagation increment; and updating historical crack tip coordinates based on the new crack tip coordinates to update the cracks, and assembling the finite element model until the I-type stress intensity factor is greater than the fracture toughness, and outputting a result in a table form.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of numerical simulation of material fatigue damage, in particular to a fatigue crack propagation simulation and life prediction method considering crack closure effect. BACKGROUND

[0002] Traditional fatigue experiments have the problems of long time consumption and high cost, and need artificial continuous monitoring of crack propagation, resulting in low efficiency. Therefore, when studying fatigue life, more and more people tend to predict structural failure through numerical simulation. At present, the mainstream finite element software is roughly divided into two categories: crack propagation direction prediction software and fatigue life prediction software. The existing software usually studies crack propagation direction and fatigue life analysis according to specific requirements respectively, and the two cannot be processed at the same time, resulting in fragmentation of function modules and a substantial increase in research time and cost. Although the ABAQUS finite element software provides functions of crack propagation path and fatigue life prediction, it is relatively complex to simultaneously perform the two calculations in the same model and the accuracy is difficult to guarantee. Meanwhile, the XFEM crack propagation module in the software omits the crack tip enhancement function, which improves the calculation efficiency but reduces the calculation accuracy, especially when dealing with mixed mode cracks, the accuracy problem is particularly prominent.

[0003] When using numerical simulation to analyze fatigue crack propagation, the crack closure effect has an important influence on fatigue life, but there are few studies on life analysis considering crack closure effect in commercial software. The existing commercial software cannot dynamically couple the crack closure effect, resulting in insufficient accuracy of the predicted crack fatigue life. At the same time, the existing commercial finite element software cannot accurately predict the crack tip position and propagation path when simulating crack propagation, which further affects the efficiency of engineering decision-making. SUMMARY

[0004] Therefore, it is necessary to provide a fatigue crack propagation simulation and life prediction method considering crack closure effect, comprising: S1: automatically inputting material parameters, boundary conditions, mesh division, pre-crack and crack opening load ratio stored in a CSV file into a finite element software, and constructing a finite element model; S2: assembling crack size in the finite element model according to user requirements and refining the area mesh of the crack, using the XFEM crack propagation module to expand the crack in the finite element model, extracting the displacement field of the area mesh at the crack tip from the ODB file, and performing nonlinear calculation considering crack closure effect based on the displacement field to obtain the mode I stress intensity factor and mode II stress intensity factor at the crack tip; S3: Calculates the crack propagation angle based on the mode I and mode II stress intensity factors and the maximum circumferential stress criterion, and calculates the new crack tip coordinates based on the crack tip coordinates, crack propagation angle, and user-defined crack propagation increment; S4: updating the historical crack tip coordinates based on the new crack tip coordinates to update the crack and assemble it into the finite element model; S5: Repeat steps S2-S4 until the mode I stress intensity factor is greater than the fracture toughness, ending the crack growth cycle and outputting the results in a table format. The results include the mode I stress intensity factor, the mode II stress intensity factor, the crack growth path coordinates, the crack growth length, and the fitted crack fatigue life.

[0005] Preferably, the material parameters include Young's modulus, Poisson's ratio, and shear modulus.

[0006] Preferably, in S2, a nonlinear calculation taking into account the crack closure effect is performed based on the displacement field to obtain the mode I stress intensity factor and the mode II stress intensity factor at the crack tip, including: The displacement field of the regional mesh at the crack tip is expressed as: ; ; in, Indicates the number of cells in the crack tip area. Nodes in x Directional displacement; Indicates the number of cells in the crack tip area. Nodes in y Directional displacement; represents the shear modulus of the material; represents the order of Williams polynomial expansion; represents the least squares fitting Undetermined coefficient of order; Indicates the first grid in the area at the crack tip The length of the node; Indicates the first grid in the area at the crack tip The angle of each node; represents the Kolosov constant; express x Rigid body translation in the direction express y Rigid body translation in the direction Indicates the rotation of the rigid body; Undetermined coefficients of least squares fitting in expressions based on displacement fields and , the stress intensity factor of type I and the stress intensity factor of type II are obtained, and the calculation formula is: ; ; wherein, represents the stress intensity factor of type I at the crack tip; represents the stress intensity factor of type II at the crack tip; represents the first-order undetermined coefficient of the least square fitting; represents the second-order undetermined coefficient of the least square fitting.

[0007] Preferably, when the model is in a plane strain state, the calculation formula of the Kolosov constant is: ; When the model is in a plane stress state, the calculation formula of the Kolosov constant is: ; wherein, represents the Kolosov constant; represents the Poisson's ratio.

[0008] Preferably, it further comprises calculating the effective stress intensity factor ratio based on the crack opening load ratio, and the calculation formula is: ; ; ; ; wherein, represents the effective stress intensity factor ratio, which is used to reflect the proportional relationship between the effective stress intensity factor and the load; represents the crack opening load ratio; represents the load ratio; represents the stress intensity factor corresponding to the maximum load; represents the stress intensity factor corresponding to the minimum load; represents the crack opening stress intensity factor; represents the stress intensity factor of type I at the crack tip; represents the theoretical stress intensity factor; represents the applied load; represents the material thickness; represents the material width; represents the ratio of the crack length to the material width.

[0009] Preferably, S4 further comprises: saving the crack opening displacement of the crack along the distance from the crack tip in the crack propagation iteration process according to the user demand.

[0010] Preferably, the process of obtaining the crack opening displacement comprises: In the crack propagation iteration process, obtaining the vertical coordinates of the corresponding two points on the upper and lower surfaces of the crack on the same vertical line, and calculating the difference between the vertical coordinates as the crack opening displacement of the corresponding position; repeating the process of obtaining the vertical coordinates of the other corresponding two points on the upper and lower surfaces of the crack and calculating the crack opening displacement of the corresponding position to obtain the crack opening displacement of each part of the crack in the corresponding crack iteration process.

[0011] Preferably, the fracture toughness is 0.85 times the fracture toughness of the material.

[0012] Beneficial effects: the method extracts the displacement field at the crack tip, and performs nonlinear calculation considering the crack closure effect based on the displacement field, to obtain the I-type stress intensity factor and the II-type stress intensity factor at the crack tip, thereby enhancing the reliability of the fatigue life prediction; at the same time, the maximum circumferential force principle is used to realize the crack path tracking, successfully solving the technical problem of synchronous processing of crack path prediction and fatigue life calculation, and realizing more accurate and efficient crack propagation and life prediction. BRIEF DESCRIPTION OF DRAWINGS

[0013] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed to be used in the embodiments or prior art description will be briefly introduced as follows. Obviously, the drawings in the following description only constitute some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained without creative labor.

[0014] Figure 1 The flow chart of the fatigue crack propagation simulation and life prediction method considering the crack closure effect in the embodiments of the present application.

[0015] Figure 2 The crack tip expansion diagram in the embodiments of the present application.

[0016] Figure 3a The crack opening displacement curve with crack opening load change in the embodiments of the present application.

[0017] Figure 3b The crack opening rate diagram representing the crack closure effect in the embodiments of the present application.

[0018] Figure 3c The average crack opening rate curve with crack opening load change in the embodiments of the present application.

[0019] Figure 3d Figure 1 is a graph of the stress intensity factor characterizing the crack closure effect in the embodiments of the present application.

[0020] Figure 4a Figure 2 is a graph of the crack opening load as a function of the crack length in the embodiments of the present application.

[0021] Figure 4b Figure 3 is a graph of the crack opening load ratio as a function of the cycle number in the embodiments of the present application.

[0022] Figure 5a Figure 4 is a graph of the comparison between the simulated and experimental life as a function of the crack length in the embodiments of the present application.

[0023] Figure 5b Figure 5 is a graph of the comparison between the simulated and experimental error as a function of the crack length in the embodiments of the present application. DETAILED DESCRIPTION

[0024] In order to make the above objectives, features and advantages of the present application more apparent, specific embodiments of the present application will be described in detail below with reference to the accompanying drawings. In the following description, numerous specific details are set forth in order to provide a thorough understanding of the present application. However, the present application can be practiced in a number of different ways beyond the specific embodiments described and it is to be understood that the present application is not limited in any way only to the embodiments disclosed herein.

[0025] In addition, the terms "first", "second", etc. are used only for the purpose of description and should not be construed as indicating or implying relative importance or implying the number of the technical features indicated. Therefore, the features defined with "first", "second", etc. can explicitly or implicitly include at least one of the features. In the description of the present application, the meaning of "a plurality of" is at least two, such as two, three, etc., unless otherwise specifically limited.

[0026] The embodiments of the present application simulate crack propagation using the extended finite element method, covering the entire process from geometry creation, job submission to post-processing, and fully implementing each step of fatigue crack propagation to fracture. The crack propagation process is continuous, and in the simulation analysis, not only the elastic-plastic behavior of the material needs to be considered, but also the change of the crack opening load ratio during the crack propagation process, which is of great significance for accurately predicting the life of the crack considering the crack closure effect. In view of the complexity of the dynamic crack propagation behavior and the limitations of the solving algorithm, the embodiments of the present application use the quasi-static simulation method of ABAQUS to perform cycle calculation, and the application programming interface (API) of the pre- and post-processing is developed using Python language, thereby effectively solving the above problems.

[0027] As Figure 1As shown, this embodiment provides a fatigue crack growth simulation and life prediction method considering the crack closure effect, including: S1: Automatically input the material parameters, boundary conditions, meshing, pre-cracks, and crack opening load ratio stored in the CSV file into the finite element software and construct the finite element model.

[0028] In this embodiment, the material is nickel-based alloy GH4169 (which can be replaced with a material studied by the user according to actual needs), and the material parameters include Young's modulus, Poisson's ratio, and shear modulus.

[0029] S2: Assemble the crack size in the finite element model according to user requirements and refine the regional mesh of the crack. Use the XFEM crack propagation module to propagate the crack in the finite element model, and extract the displacement field of the regional mesh at the crack tip from the ODB file. Based on the displacement field, perform nonlinear calculations considering the crack closure effect to obtain the mode I and mode II stress intensity factors at the crack tip.

[0030] Specifically, a nonlinear calculation considering the crack closure effect is performed based on the displacement field to obtain the mode I stress intensity factor and mode II stress intensity factor at the crack tip, including: The displacement field of the regional mesh at the crack tip is expressed as: ; ; in, Indicates the number of cells in the crack tip area. Nodes in x Directional displacement; Indicates the number of cells in the crack tip area. Nodes in y Directional displacement; represents the shear modulus of the material; represents the order of Williams polynomial expansion; represents the least squares fitting The coefficient of unknown order is related to the stress intensity factor; Indicates the first grid in the area at the crack tip The length of the node; Indicates the first grid in the area at the crack tip The angle of each node; represents the Kolosov constant; express x Rigid body translation in the direction express y Rigid body translation in the direction Indicates the rotation of the rigid body; Kolosov constant is a material constant determined by Poisson's ratio . Its value depends on whether the model is in plane strain or plane stress state, and the parameter is a key parameter in the calculation formula of the displacement field of the region mesh at the crack tip and , and directly affects the accuracy of the results. When the model is in plane strain state, the calculation formula of the Kolosov constant is: ; when the model is in plane stress state, the calculation formula of the Kolosov constant is: , wherein represents the Poisson's ratio; Based on the least square fitting of the undetermined coefficients in the expression of the displacement field and , the mode I stress intensity factor and the mode II stress intensity factor are obtained, and the calculation formula is: ; ; , wherein represents the mode I stress intensity factor at the crack tip; represents the mode II stress intensity factor at the crack tip; represents the first-order undetermined coefficient of the least square fitting; represents the second-order undetermined coefficient of the least square fitting.

[0031] Before running the script, the user can modify the mesh size and crack propagation increment to effectively control the crack closure effect. The model building process is integrated into the plug-in in this embodiment, so that the user can analyze the model multiple times, thereby reducing repetitive operations and improving work efficiency.

[0032] Further, the effective stress intensity factor ratio based on the crack opening load ratio is calculated, and the calculation formula is: ; ; ; ; , wherein represents the effective stress intensity factor ratio, which is used to reflect the proportional relationship between the effective stress intensity factor and the load; represents the crack opening load ratio; represents the load ratio; represents the maximum load corresponding to the stress intensity factor; represents the minimum load corresponding to the stress intensity factor; represents the crack opening stress intensity factor; represents the mode I stress intensity factor at the crack tip; represents the effective stress intensity factor amplitude; represents the stress intensity factor amplitude; represents the theoretical stress intensity factor; represents the applied load; Indicates material thickness; Indicates the material width; It represents the ratio of crack length to material width, , Indicates the crack length; represents the theoretical stress intensity factor, used for comparison and verification with the stress intensity factor calculated using the displacement field proposed in this embodiment. Different specimen configurations require different theoretical formulas. This embodiment uses a standard compact tensor (CT) specimen designed according to GB / T 6398-2000, "Test Method for Fatigue Crack Growth Rate of Metallic Materials." This embodiment calculates the stress intensity factor based on the displacement field, eliminating the need to modify the theoretical formula based on the specimen configuration. This is an advantage of the method proposed in this embodiment.

[0033] In this embodiment, the Paris formula expresses the fatigue crack growth rate da / dN as ΔK function of , so when the crack extension increment Δa When it is small, the number of cycles can be calculated ΔN , that is, crack extension Δa Required number of cycles ΔN (i.e. fatigue life), the calculation formula is:

[0034] in, represents the crack growth increment under cyclic loading; represents the fatigue life corresponding to the crack growth increment under cyclic loading; Indicates in Fatigue life corresponding to crack growth increment under sub-cyclic loading: Indicates in The stress intensity factor amplitude corresponding to the crack growth increment under the sub-cycle loading. Before the crack fails, the fatigue life corresponding to the crack growth increment under each cyclic loading is accumulated, and the sum obtained is the total fatigue life of the crack.

[0035] S3: Calculate the crack propagation angle based on the mode I stress intensity factor and the mode II stress intensity factor, and based on the maximum circumferential stress criterion. Calculate the new crack tip coordinates based on the crack tip coordinates, the crack propagation angle, and the user-defined crack propagation increment.

[0036] Specifically, the maximum circumferential stress criterion assumes that the crack propagation occurs from the crack tip in the direction of the maximum circumferential tensile stress, and the crack propagation angle is calculated as: ; wherein, denotes the inverse cosine function; denotes the mode I stress intensity factor at the crack tip; denotes the mode II stress intensity factor at the crack tip. The crack propagation angle is obtained from the mode I stress intensity factor and the mode II stress intensity factor After the crack propagation angle is obtained, the crack tip coordinates of the next propagation (x , ) are obtained according to the trigonometric function, and the crack propagation continues in this way until the mode I stress intensity factor reaches the critical stress intensity factor (input by the user according to the material being studied) or 0.85 times the fracture toughness (input by the user according to the material being studied), the crack fails, the specimen is determined to be broken, the program stops operating and outputs the mechanical parameters.

[0037] In the above, Figure 2 the calculation of the new crack tip coordinates includes: ; ; wherein, denotes the horizontal coordinate of the crack tip of the next crack propagation; denotes the vertical coordinate of the crack tip of the next crack propagation; denotes the horizontal coordinate of the crack tip of the current crack propagation; denotes the vertical coordinate of the crack tip of the current crack propagation; denotes the user-defined crack propagation increment; denotes the crack propagation angle.

[0038] S4: Update the historical crack tip coordinates based on the new crack tip coordinates to update the crack, and assemble into the finite element model. The historical crack tip coordinates are encapsulated as a set, and the crack is updated by drawing the crack path, and finally the crack propagation until the fracture is completed.

[0039] In the embodiment, the method further comprises: saving the crack opening displacement of the crack at each position along the distance from the crack tip in each iteration of the crack propagation.

[0040] Further, the process of obtaining the crack opening displacement comprises: In each iteration of the crack propagation, obtaining the vertical coordinates of the corresponding two points on the upper and lower surfaces of the crack on the same vertical line, and calculating the difference between the vertical coordinates as the crack opening displacement at the corresponding position; repeating the process of obtaining the vertical coordinates of the other corresponding two points on the upper and lower surfaces of the crack and calculating the crack opening displacement at the corresponding position to obtain the crack opening displacement (COD) at each position in the corresponding iteration of the crack.

[0041] In the embodiment, the crack opening ratio at the corresponding position can be calculated based on the crack opening displacement at each position, comprising: In each iteration of the crack propagation, dividing the crack opening from the beginning to the maximum opening into M states, and arranging pairs of virtual displacement extensometers (VDEs) on the upper and lower surfaces of the crack on the same vertical line, and calculating the difference between the vertical coordinates of the pairs of VDEs as the virtual opening displacement. According to the principle of similar triangles, the M states in the crack opening process can be obtained: ; wherein, represents the crack opening displacement at the first pair of VEDs position in the Mth state; represents the crack opening displacement at the second pair of VEDs position in the Mth state; denote the crack opening displacement at the first pair of VEDs position in the Mth state when the crack is fully opened; , then: ; wherein, represents the crack opening displacement at the first pair of VEDs position in the Mth state when the crack is fully opened; represents the crack opening displacement at the second pair of VEDs position in the Mth state when the crack is fully opened; represents the crack opening displacement at the first pair of VEDs position in the Mth state; represents the crack opening displacement at the second pair of VEDs position in the Mth state; In the M states, there are: ; If v pairs of VDEs are arranged behind the crack tip, then in the crack opening state 1: ; wherein, represents the crack opening displacement of the vth pair of VEDs position in the 1st state when the crack is fully opened; If there are M states of crack opening, then in each state: ; wherein, represents the crack opening displacement of the vth pair of VEDs position in the Mth state when the crack is fully opened; represents the maximum crack opening displacement of the vth pair of VEDs position; It can be concluded that when M and v are large enough, at any moment t when the crack is just fully opened to any phase of the maximum crack opening, the ratio of the crack opening displacement at any position v behind the crack tip to the maximum crack opening displacement at this position is a constant value, i.e. the crack opening ratio at the corresponding position, and the calculation formula is: ; wherein, represents the crack opening ratio of the vth pair of VEDs position at t, The value range of is [0, 1], which represents the opening degree of the crack tip in the opening-closing process of the crack in a single cycle; represents the crack opening displacement of the vth pair of VEDs position at t.

[0042] Further, at the same moment of crack opening and closing, the crack opening ratios of different VEDs positions are a constant value, therefore, the average crack opening ratio at t is: ; wherein, represents the average crack opening ratio at t; represents the crack opening ratio of the vth pair of VEDs position at t; A represents the total number of VDEs.

[0043] Through the developed finite element product, according to the demand, the crack tip displacement field at a certain cycle in the crack propagation process can be extracted, and the crack opening displacement, crack opening ratio, average crack opening ratio and stress intensity factor under different loads along the crack tip distance can be obtained. The crack closure effect is characterized by the four methods, thereby improving the accuracy and reliability of the crack propagation path prediction and fatigue life analysis. As shown in Figure 3a , Figure 3b , Figure 3c , Figure 3d , the four methods characterize the crack opening load.

[0044] S5: repeat steps S2-S4 until the mode I stress intensity factor is greater than the fracture toughness, end the crack propagation cycle, and output the results in table form, including the mode I stress intensity factor, the mode II stress intensity factor, the crack propagation path coordinates, the crack propagation length, and the fitted crack fatigue life.

[0045] In this embodiment, the fracture toughness is 0.85 times the fracture toughness of the material. When the stress intensity factor calculated by the crack tip displacement field is greater than or equal to 0.85 times the fracture toughness of the material, it indicates that the crack has failed, the crack propagation is stopped, and the data is output. (Supplementary explanation: The fracture toughness is a material property that indicates the material's ability to resist crack propagation. The fracture toughness of different materials is different, and can be determined by consulting the "Stress Intensity Factor Handbook". The user needs to input the corresponding value according to the material studied.) In this embodiment, the output results are shown in Table 1; Table 1 Output results table ; In Table 1, job represents the index of finite element analysis; KI and KII represent the mode I stress intensity factor and the mode II stress intensity factor at the crack tip, respectively; K_theo represents the theoretical stress intensity factor; represents the crack propagation angle; CrackTipX and CrackTipY represent the x coordinate and y coordinate of the crack tip, respectively; crackLength represents the crack length; Δa represents the crack propagation length; ΔN represents the fitted crack fatigue life corresponding to the propagation Δa; N represents the total fitted crack fatigue life of the crack propagation to the current length; K_max and K_min represent the stress intensity factors corresponding to the maximum load and the minimum load , respectively; K_open represents the crack opening stress intensity factor; ΔK_eff represents the effective stress intensity factor amplitude; R_open represents the crack opening load ratio.

[0046] Based on the results output by the finite element analysis considering the crack closure effect calculation method, the simulation values of the crack opening load and the crack life are compared with the experimental data, as shown in Figure 4a 、 Figure 4b the comparison chart of simulated and experimental crack opening loads and Figure 5a 、 Figure 5b the comparison chart of simulated and experimental lives, which shows the high consistency between the results of the product and the actual experimental results, and verifies the accuracy of the method in crack life prediction.

[0047] The fatigue crack propagation simulation and life prediction method considering the crack closure effect provided in this embodiment has the following beneficial effects: 1. The method aims to overcome the shortcomings of existing commercial software in the synchronous implementation of crack path tracking and life calculation, based on ABAQUS API for deep customization development, and builds an embedded fatigue crack propagation analysis method. The method first realizes the integrated simulation and prediction of crack path and its expansion life under the consideration of crack closure effect. The method covers the whole process of modeling-computation-postprocessing, supports automatic tracking of crack path and adaptive updating of mesh through dynamic coupling calculation of crack closure effect, realizes real-time output of fatigue life data, and significantly improves the automation and real-time response capability in the simulation process. It effectively overcomes the problems of low calculation efficiency and dependence on manual intervention of traditional tools, meets the dual requirements of high precision and high efficiency in crack behavior prediction of modern engineering simulation, breaks through the limitations of existing commercial software functions, and innovatively solves the industry technical problems that traditional tools cannot realize synchronous processing of crack path tracking and life calculation.

[0048] 2. By developing a fatigue model plug-in considering crack closure effect, the problem that the existing ABAQUS platform-based two-dimensional crack propagation model based on XFEM cannot calculate the stress intensity factor is solved. The product covers all stages from pre-processing to post-processing, can automatically create a finite element model, perform fatigue crack propagation analysis and output fatigue data, thereby improving the automation and integrity of the analysis process and effectively expanding the existing ABAQUS platform.

[0049] 3. The existing numerical simulation fatigue crack method ignores the significant impact of crack closure effect on fatigue life. The method innovatively combines crack closure effect through numerical simulation method to accurately capture the changes of crack propagation path and life, making up for the shortcomings of existing simulation methods in this regard. Key technical points include that during product operation, crack tip related mechanical data can be extracted in real time according to user demand. The crack closure effect is characterized by COD-load curve, COR-load curve, and stress intensity factor-load curve, which is the first time to characterize the crack closure effect through a finite element model.

[0050] 4. The existing finite element software cannot simultaneously perform crack propagation path prediction and fatigue life analysis in the same model, which leads to low calculation efficiency and cannot guarantee the accuracy of the results. The present invention enhances the reliability of fatigue life prediction by extracting the crack tip displacement field and combining the Williams model and the least squares method to accurately calculate the stress intensity factor. At the same time, the maximum circumferential force principle is used to realize crack path tracking, successfully solving the technical problem of synchronous processing of crack path prediction and fatigue life calculation, and realizing more accurate and efficient crack propagation and life prediction.

[0051] 5、Through the comparison and analysis with the fatigue test results, the accuracy and reliability of the plug-in are verified, and the high-precision calculation ability of the plug-in in the fatigue crack propagation simulation is proved. This provides a theoretical basis and technical support for the application of the plug-in in practical engineering.

[0052] The technical features of the above-described embodiments can be combined in any manner. To make the description concise, all possible combinations of the technical features in the above-described embodiments are not described, but as long as the combinations of the technical features do not exist contradictory, they should be considered within the scope of the present disclosure.

[0053] The above-described embodiments only express several implementation manners of the present application, and the description is relatively specific and detailed, but it should not be understood as a limitation on the patent scope of the application. It should be pointed out that for ordinary skilled persons in the art, several modifications and improvements can be made without departing from the concept of the present application, and these all belong to the protection scope of the present application. Therefore, the patent protection scope of the present application should be subject to the appended claims.

Claims

1. A method of fatigue crack growth simulation and life prediction considering crack closure effect, characterized in that, The method comprises the following steps: S1: automatically inputting material parameters, boundary conditions, mesh division, pre-crack and crack opening load ratio stored in a CSV file into a finite element software, and constructing a finite element model; S2: assembling crack size in the finite element model and refining the area mesh of the crack according to user requirements, expanding the crack in the finite element model by using an XFEM crack expansion module, extracting a displacement field of the area mesh at the crack tip from an ODB file, and performing nonlinear calculation considering crack closure effect based on the displacement field to obtain mode I stress intensity factor and mode II stress intensity factor at the crack tip; S3: calculating a crack expansion angle based on the mode I stress intensity factor and the mode II stress intensity factor and based on a maximum circumferential stress criterion, and calculating new crack tip coordinates based on the crack tip coordinates, the crack expansion angle and a user-defined crack expansion increment; S4: updating the historical crack tip coordinates based on the new crack tip coordinates to update the crack, and assembling the crack into the finite element model; S5: repeatedly performing steps S2-S4 until the mode I stress intensity factor is greater than the fracture toughness, ending the crack expansion cycle, and outputting the results in a table form, wherein the results include the mode I stress intensity factor, the mode II stress intensity factor, crack expansion path coordinates, crack expansion length and fitted crack fatigue life.

2. The fatigue crack growth simulation and life prediction method considering crack closure effect according to claim 1, characterized in that, The material parameters include Young's modulus, Poisson's ratio and shear modulus.

3. The fatigue crack growth simulation and life prediction method considering crack closure effect according to claim 1, characterized in that, In S2, the nonlinear calculation considering crack closure effect based on the displacement field is performed to obtain the mode I stress intensity factor and the mode II stress intensity factor at the crack tip, which comprises the following steps: The displacement field of the area mesh at the crack tip is expressed as: ； ； in, Indicates the number of cells in the crack tip area. Nodes in x Directional displacement; Indicates the number of cells in the crack tip area. Nodes in y Directional displacement; represents the shear modulus of the material; represents the order of Williams polynomial expansion; represents the least squares fitting Undetermined coefficient of order; Indicates the first grid in the area at the crack tip The length of the node; Indicates the first grid in the area at the crack tip The angle of each node; represents the Kolosov constant; express x Rigid body translation in the direction express y Rigid body translation in the direction Indicates the rotation of the rigid body; Undetermined coefficients of least square fitting in expression based on displacement field and The stress intensity factors of type I and type II are obtained, and the calculation formula is as follows: ； ； wherein, K1represents the mode I stress intensity factor at the crack tip; K2represents the mode II stress intensity factor at the crack tip; a1represents the first order undetermined coefficient of the least square fitting; a2represents the second order undetermined coefficient of the least square fitting.

4. The fatigue crack growth simulation and life prediction method considering crack closure effect according to Claim 3, wherein When the model is in a plane strain state, the calculation formula of the Kolosov constant is: ； When the model is in a plane stress state, the calculation formula of the Kolosov constant is: ； wherein represents the Kolosov constant; represents the Poisson's ratio.

5. The method of simulating fatigue crack growth and predicting life with consideration of crack closure effect according to claim 1, wherein, The method further comprises calculating an effective stress intensity factor ratio based on the crack opening load ratio, and the calculation formula is: ； ； ； ； wherein, represents an effective stress intensity factor ratio, used to embody the proportional relationship between the effective stress intensity factor and the load; represents a crack opening load ratio; represents a load ratio; represents a stress intensity factor corresponding to the maximum load; represents a stress intensity factor corresponding to the minimum load; represents a crack opening stress intensity factor; represents a mode I stress intensity factor at the crack tip; represents a theoretical stress intensity factor; represents an applied load; represents a material thickness; represents a material width; represents a ratio of crack length to material width.

6. The method of simulating fatigue crack growth and predicting life with consideration of crack closure effect according to claim 1, wherein, S4 further comprises: saving crack opening displacement of the crack along the crack tip distance in the crack expansion iteration process according to user requirements.

7. The method of simulating fatigue crack growth and predicting life with consideration of crack closure effect according to claim 6, wherein The process of obtaining the crack opening displacement comprises the following steps: In the crack expansion iteration process, the vertical coordinates of the corresponding two points on the upper and lower surfaces of the crack located on the same vertical line are obtained, and the difference between the vertical coordinates is calculated as the crack opening displacement at the corresponding position; the vertical coordinates of other corresponding two points on the upper and lower surfaces of the crack are repeatedly obtained, and the crack opening displacement at the corresponding position is calculated to obtain the crack opening displacement at each position of the crack in the corresponding crack iteration process.

8. The method of simulating fatigue crack growth and predicting life with consideration of crack closure effect according to claim 1, wherein, The fracture toughness is 0.85 times the fracture toughness of the material.

Images