Crack stress intensity factor calculation method based on three-dimensional weight function method

Abstract

The invention discloses a crack stress intensity factor calculation method based on a three-dimensional weight function method, and belongs to the technical field of reliability. Based on a complex variable function Taylor series expansion method, a complex variable finite element program is operated to obtain a weight function of given crack geometry, and then a stress intensity factor of the crack is obtained through stress on a crack surface and weight function integration. Specifically, the method comprises the following steps: establishing a crack-containing finite element model and dividing grids; carrying out perturbation on a crack tip unit node on an imaginary axis; running a complex variable finite element program to obtain an imaginary part of a node crack opening direction on the crack surface so as to obtain a weight function on the crack surface; and solving the stress intensity factor of the corner of the connecting pipe by integrating the weight function of the crack surface and the stress of the crack surface. Compared with a two-dimensional weight function, the method is expanded to three-dimensional elliptical cracks, and the application range is wider; by adopting the method, the calculation result of solving the stress intensity factor of the corner of the connecting pipe under the action of complex load is accurate, and the use is simple and effective.

Description

Technical Field

[0001] The invention belongs to the technical field of reliability, and in particular relates to a method for calculating a crack stress intensity factor based on a three-dimensional weight function method. Background Art

[0002] In the field of engineering, the evaluation of stress intensity factors for cracked structures is crucial. Crack generation is a common and important problem in the preparation of materials and the processing and manufacturing of parts. It is worth noting that the load when fracture occurs is usually much lower than the traditional static strength design load level. In this context, fracture mechanics has become an important discipline used to quantitatively answer important questions related to cracks. As the driving force of cracks, the stress intensity factor K plays an important role in fracture mechanics and is of great significance for evaluating crack propagation and the fatigue life of structures. Therefore, the accurate calculation of the K value has become the core task of fracture mechanics.

[0003] For a long time, the international fracture science field has developed a variety of theories and methods to solve the calculation problem of stress intensity factor. However, due to the singularity of the stress field at the crack tip, pure theoretical solutions are usually only applicable to a very small number of ideal crack geometries and load conditions, such as the central crack of an infinite plate or periodic collinear cracks, embedded elliptical cracks, and semi-infinite cracks. For the vast majority of finite body crack geometry stress intensity factor studies in engineering practice, most of the focus is on approximate analytical solutions, experimental analysis, and numerical calculation analysis. There are currently two methods used: weight function method (WFM) and numerical methods (finite element method (FEM), boundary element method (BEM)). Although numerical methods have shown strong analytical capabilities in dealing with complex crack geometries, numerical methods such as finite element method are not technically difficult to solve crack problems, and for some special field crack problems, this is an attractive option. However, when it comes to complex structures containing a large number of critical crack locations, various numerical methods face huge challenges in modeling and calculation due to the singularity of the crack tip field and the continuous change of crack size.

[0004] The weight function method is a powerful tool for solving stress intensity factors. The weight function considers the load conditions and geometric conditions separately, only includes the geometric characteristics of the crack body, and is not affected by the load. Once the weight function is determined, it can be used without restriction to calculate K under any load condition.

[0005] On this basis, the complex variable Taylor series expansion weight function (WCTSE) proposed by Wanger and Millwater for solving the two-dimensional crack geometric weight function is combined with the three-dimensional crack weight function method proposed by Rice and extended to three-dimensional elliptical cracks.

[0006] Compared with the original weight function method, the three-dimensional complex function Taylor series expansion weight function method has a wider application range in that the applicable crack form has increased from two-dimensional cracks to three-dimensional cracks. The stress distribution form has increased from a single crack depth direction to two dimensions of crack depth and length. The three-dimensional complex function Taylor series expansion weight function method is applicable to cracks with complex stress gradients, thus providing a certain reference for the analysis of fracture mechanics parameters under the framework of structural integrity assessment. Summary of the invention

[0007] In view of this, the present invention proposes a method for calculating crack stress intensity factor based on a three-dimensional weight function method to improve the calculation efficiency of the crack stress intensity factor. The specific technical solution of the present invention is as follows:

[0008] A method for calculating crack stress intensity factor based on a three-dimensional weight function method comprises the following steps:

[0009] S1: For the structure with cracks, a finite element model is established and meshed, and the unit and node information of the corresponding geometric model is derived;

[0010] S2: Read the unit and node information of the geometric model in step S1, and perturb the unit near the required solution point at the crack tip along the crack opening direction;

[0011] S3: Run the complex variable finite element program to calculate the node displacement U r , read the imaginary part of the crack surface normal displacement solution Im[U r (P'; x, y)], x, y represent the abscissa and ordinate values ​​of any point in the crack surface, respectively, and P' is the crack front point;

[0012] S4: Combine the imaginary part of the normal displacement solution on the crack surface obtained in step S3, Im[U r (P'; x, y)], and the normal stress σ of the crack surface under the reference load r (x, y) is calculated to obtain the reference stress intensity factor K r (P'), the reference load is the crack surface subjected to a uniformly distributed stress of 1 MPa;

[0013] S5: obtaining an imaginary crack surface normal stress distribution σ(x, y) at a position corresponding to the crack surface in the cracked finite element model constructed in step S1 on the crack-free finite element model, wherein the load on the model is the actual working condition load;

[0014] S6: Combine the imaginary part Im[U r (P'; x, y)], the reference stress intensity factor K obtained in step S4 r (P') and the normal stress distribution σ of the crack surface obtained in step S5 r(x, y), and calculate the stress intensity factor K(P') under the required load.

[0015] Furthermore, the specific process of step S1 is:

[0016] Finite element software ABAQUS was used to establish the crack model, divide the mesh and export the inp file.

[0017] Furthermore, the specific process of step S2 is:

[0018] MATLAB software is used to read the unit and node information of the geometric model in S1, and the unit nodes near the desired crack tip are perturbed along the crack opening direction so that the coordinates of these nodes carry complex variable information. The perturbation step h' is orthogonally decomposed into three directions: x, y, and z. The original node coordinates are (x0, y0, z0), and the node coordinates after perturbation are (x0+ix i ,y0+iy i ,z0+iz i ), i is the imaginary unit.

[0019] Furthermore, the specific process of step S3 is:

[0020] Run the finite element program to calculate the unit stiffness matrix of all units in the model in step S1, and then combine the unit stiffness matrix into an overall stiffness matrix; introduce boundary conditions and loading conditions, and the loading condition is to apply a reference load σ on the crack surface r (x, y); calculate the node displacement according to the equilibrium equation P = KU, and read the imaginary part of the solution of the normal displacement on the crack surface Im[U r (P'; x, y)].

[0021] Furthermore, the specific process of step S4 is:

[0022] S4.1. Theoretical formula for three-dimensional crack weight function:

[0023]

[0024] The Taylor series of complex function is used to expand the weight function, and at the same time, substitution is used. Let x = racosθ, y = rcsinθ, r is the polar diameter in polar coordinates, θ is the polar angle in polar coordinates, a and c are the crack depth and crack half-length respectively, then:

[0025]

[0026] Perturbation value h = w·h', w is the perturbation unit width;

[0027] S4.2, using Gaussian integral calculation, then:

[0028]

[0029] Where u i ,v j They are the integration points in the polar radial direction and polar angular direction, within the integration limit [-1, +1]; several integration points have corresponding integration weights A i , B j ; Stress σ at the integration point of the reference load r,ij and the imaginary value Im of the normal displacement of the crack surface at the integration point of the crack model ij [U r (P'; x, y)] is located on the crack surface ((u i +1) / 2,π(v j +1) / 4) places;

[0030] S4.3. Based on the self-consistency principle, the reference stress intensity factor K r (P') is equal to:

[0031]

[0032] Furthermore, the specific process of step S5 is:

[0033] The finite element software ABAQUS was used to establish a crack-free model, and the finite element method was used to obtain the normal stress distribution σ(x, y) of the crack surface at the imaginary position of the crack surface under the required load.

[0034] Furthermore, the specific process of step S6 is:

[0035] Combined with the imaginary part Im[U r (P'; x, y)], the reference stress intensity factor K obtained in step S4 r (P') and the normal stress distribution σ(x, y) of the crack surface obtained in step S5 are expanded using the complex variable Taylor series weight function, and at the same time, substitution is used, setting x = ra cosθ, y = rc sinθ, and substituting them into formula (4) to calculate and obtain the stress intensity factor K(P') of the required point under the required load.

[0036] In summary, this application includes the following beneficial technical effects:

[0037] 1. The method of the present invention adopts the complex variable Taylor series expansion method to calculate the weight function on the crack surface, and by integrating the stress on the crack surface and the weight function, the stress intensity factor of the crack under complex loads is calculated at high speed, thereby providing a basis for linear elastic fracture mechanics analysis in damage tolerance assessment.

[0038] 2. The present invention is based on a two-dimensional weight function and extends it to a three-dimensional elliptical crack, which has a wider range of applications.

[0039] 3. The crack stress intensity factor calculation method based on the three-dimensional weight function method of the present invention is more accurate in solving the calculation results of the stress intensity factor under complex loads compared to the finite element method, and the calculation method is simple and effective. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 is a flow chart of the method of the present invention;

[0041] Figure 2 Schematic diagram of geometric parameters of a nozzle corner crack in an embodiment of the present invention;

[0042] Figure 3 is a coordinate schematic diagram of a crack at a nozzle corner in an embodiment of the present invention;

[0043] Figure 4 Schematic diagram of perturbation nodes of a pipe corner crack grid in an embodiment of the present invention;

[0044] Figure 5 It is a comparison of the results of calculating the stress intensity factor by the method of this embodiment and the finite element method respectively when the corner crack of the pipe is subjected to internal pressure load in the embodiment of the present invention; wherein 5a, 5b and 5c are comparisons of the calculation results at the crack front point A, point B and point C respectively. DETAILED DESCRIPTION

[0045] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments of the specification. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0046] On the contrary, the present invention covers any substitution, modification, equivalent method and scheme made on the essence and scope of the present invention as defined by the claims. Further, in order to make the public have a better understanding of the present invention, some specific details are described in detail in the detailed description of the present invention below. Those skilled in the art can fully understand the present invention without the description of these details.

[0047] See also Figures 1 to 5 Taking the pipe corner structure with inner corner crack as an example, this example provides a method for calculating the crack stress intensity factor based on the three-dimensional weight function method, which includes the following steps:

[0048] The crack finite element model is established based on the geometric dimensions of the nozzle corner structure with inner corner cracks. The structural geometric dimension parameters are as follows: Figure 1 As shown, R i=500mm, other geometric parameters are calculated according to the dimensionless parameters in Table 1. Elastic modulus E = 204GPa, Poisson's ratio μ = 0.3. The crack surface is subjected to a uniformly distributed stress of 1MPa. The crack geometry model is established using ABAQUS finite element software.

[0049] Table 1 Geometric parameters

[0050]

[0051] The cracked cylinder surface point A, the deepest point B and the nozzle surface point C are perturbed respectively, as follows: Figure 3 As shown. Figure 4 As shown, w is the crack width, and h' is 10 times the crack depth a. -10 . Run the complex variable finite element code. Obtain the imaginary part of the displacement at the integration point. Calculate the reference stress intensity factor solution K under the reference load r (P').

[0052] A crack-free model was established and subjected to an internal pressure load of 1 MPa. The normal stress distribution on the crack surface was obtained by finite element analysis, and the normal stress at the integral point was extracted. The stress intensity factors of the crack front cylinder surface point A, the deepest point B, and the nozzle surface point C of the nozzle corner crack were obtained by Gaussian integral calculation when the nozzle corner crack was subjected to an internal pressure of 1 MPa.

[0053] The stress intensity factor at the crack tip is calculated using the method of the present invention and the traditional finite element method based on J integral. When the model mesh density is the same, the results calculated by the two methods are as follows: Figure 5 As shown, it can be seen that the method of the present invention has a high calculation accuracy, and its relative error with the finite element solution is within 10% under the action of internal pressure load.

[0054] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for calculating crack stress intensity factor based on a three-dimensional weight function method, characterized in that: The following steps are involved: S1: For the structure with cracks, a finite element model is established and meshed, and the unit and node information of the corresponding geometric model is derived; S2: Read the unit and node information of the geometric model in step S1, and perturb the unit near the required solution point at the crack tip along the crack opening direction; S3: Run the complex variable finite element program to calculate the node displacement U r , read the imaginary part of the crack surface normal displacement solution Im[U r (P'; x, y)], x, y represent the abscissa and ordinate values ​​of any point in the crack surface, respectively, and P' is the crack front point; S4: Combine the imaginary part of the normal displacement solution on the crack surface obtained in step S3, Im[U r (P'; x, y)], and the normal stress σ of the crack surface under the reference load r (x, y) is calculated to obtain the reference stress intensity factor K r (P'), the reference load is the crack surface subjected to a uniformly distributed stress of 1 MPa; S5: Establishing a crack-free finite element model, obtaining an imaginary crack surface normal stress distribution σ(x, y) at a position on the crack-free finite element model corresponding to the crack surface in the cracked finite element model constructed in step S1, wherein the load on the model is the actual working condition load; S6: Combine the imaginary part Im[U r (P'; x, y)], the reference stress intensity factor K obtained in step S4 r (P') and the normal stress distribution σ of the crack surface obtained in step S5 r (x, y), and calculate the stress intensity factor K(P') under the required load.

2. A method for calculating crack stress intensity factor based on three-dimensional weight function method according to claim 1, characterized in that: The specific process of step S2 is: Read the unit and node information of the geometric model in S1, and perturb the unit nodes near the desired crack tip along the crack opening direction, so that the coordinates of these nodes have complex variable information. The perturbation step h' is orthogonally decomposed into three directions: x, y, and z. The original node coordinates are (x0, y0, z0), and the perturbed node coordinates are (x0+ix i ,y0+iy i ,z0+iz i ), i is the imaginary unit.

3. A method for calculating crack stress intensity factor based on three-dimensional weight function method according to claim 1, characterized in that: The specific process of step S3 is: Finite element analysis software is used to calculate the unit stiffness matrix of all units in the model in step S1, and then the unit stiffness matrix is ​​combined into an overall stiffness matrix; boundary conditions and loading conditions are introduced, and the loading condition is to apply a reference load σ on the crack surface r (x, y); calculate the node displacement according to the equilibrium equation P = KU, and read the imaginary part of the solution of the normal displacement on the crack surface Im[U r (P'; x, y)].

4. A method for calculating crack stress intensity factor based on three-dimensional weight function method according to claim 1, characterized in that: The specific process of step S4 is: S4.

1. Theoretical formula for three-dimensional crack weight function: The Taylor series of complex function is used to expand the weight function, and at the same time, substitution is used. Let x = racosθ, y = rcsinθ, r is the polar diameter in polar coordinates, θ is the polar angle in polar coordinates, a and c are the crack depth and crack half-length respectively, then: Perturbation value h = w·h', w is the perturbation unit width; S4.2, using Gaussian integral calculation, then: Where u i ,v j They are the integration points in the polar radial direction and polar angular direction, within the integration limit [-1, +1]; several integration points have corresponding integration weights A i ,B j ; Stress σ at the integration point of the reference load r,ij and the imaginary value Im of the normal displacement of the crack surface at the integration point of the crack model ij [U r (P'; x, y)] is located on the crack surface ((u i +1) / 2,π(v j +1) / 4) places; S4.

3. Based on the self-consistency principle, the reference stress intensity factor K is obtained. r (P') is equal to:

5. A method for calculating crack stress intensity factor based on three-dimensional weight function method according to claim 1, characterized in that: The specific process of step S5 is: Finite element analysis software was used to establish a crack-free model, and the finite element method was used to obtain the normal stress distribution σ(x, y) of the crack surface at the imaginary position of the crack surface under the required load.

6. A method for calculating crack stress intensity factor based on a three-dimensional weight function method according to any one of claims 2 to 5, characterized in that: The specific process of step S6 is: Combined with the imaginary part Im[U r (P'; x, y)], the reference stress intensity factor K obtained in step S4 r (P') and the normal stress distribution σ(x, y) of the crack surface obtained in step S5 are expanded using the complex variable Taylor series weight function, and at the same time, substitution is used, setting x = racosθ, y = rcsinθ, and substituting them into formula (4) to calculate and obtain the stress intensity factor K(P') of the required point under the required load.

Images