A method for calculating crack stress intensity factor based on three-dimensional weighting function method

Abstract

This invention discloses a method for calculating crack stress intensity factor based on a three-dimensional weighted function method, belonging to the field of reliability technology. Based on the Taylor series expansion method of complex functions, a complex variable finite element program is run to obtain the weighting function for a given crack geometry. Then, the stress intensity factor of the crack is obtained by integrating the stress on the crack surface with the weighting function. Specifically, a finite element model containing the crack is established and meshed; the nodes of the crack tip are perturbed on the imaginary axis; the imaginary part of the crack opening direction at the nodes on the crack surface is obtained by running the complex variable finite element program, thus obtaining the weighting function on the crack surface; the stress intensity factor at the nozzle corner is solved by integrating the weighting function and the stress on the crack surface. Compared to two-dimensional weighting functions, this method extends to three-dimensional elliptical cracks, broadening its application range; the method provides accurate calculation results for the stress intensity factor at the nozzle corner under complex loads, and is simple and effective to use.

Description

Technical Field

[0001] This invention belongs to the field of reliability technology, specifically relating to a method for calculating crack stress intensity factor based on the three-dimensional weighting function method. Background Technology

[0002] In engineering, assessing the stress intensity factor of cracked structures is crucial. Crack formation is a common and significant problem during material preparation and component manufacturing. Notably, the load at fracture is typically much lower than the conventional static strength design load level. Against this backdrop, fracture mechanics has become an important discipline for quantitatively answering critical questions related to cracks. As the driving force of cracks, the stress intensity factor K plays a vital role in fracture mechanics, significantly impacting crack propagation and the fatigue life of structures. Therefore, accurately calculating the K value is a core task of fracture mechanics.

[0003] For a long time, the international fracture science field has developed various theories and methods to solve the problem of calculating stress intensity factors. However, due to the singularity of the stress field at the crack tip, purely theoretical solutions are usually only applicable to a very small number of ideal crack geometries and load conditions, such as central cracks or periodic collinear cracks in infinitely large plates, embedded elliptical cracks, and semi-infinite cracks. For most finite body crack geometry stress intensity factor studies in engineering practice, most research focuses on approximate analytical solutions, experimental analysis, and numerical calculation analysis. Currently, two methods are used: the weight function method (WFM) and numerical methods (finite element method (FEM) and boundary element method (BEM)). Although numerical methods demonstrate powerful analytical capabilities in handling complex crack geometries, such as the finite element method, solving crack problems is technically not difficult, making it a very attractive option for some special field crack problems. However, when dealing with complex structures containing numerous critical crack locations, the singularity of the crack tip field and the continuous variation in crack size pose significant challenges to the modeling and calculation of various numerical methods.

[0004] The weighted function method is a powerful tool for solving the stress intensity factor. The weighted function considers the load conditions and geometric conditions separately, including only the geometric characteristics of the cracked body and remaining unaffected by the load. Once the weighted function is determined, it can be used without restriction for calculating K under any load condition.

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

[0006] Compared to the original weighted function method, the three-dimensional complex variable function Taylor series expansion weighted function method expands the applicable crack types from two-dimensional to three-dimensional, thus broadening its application range. The stress distribution pattern expands from a single crack depth direction to two dimensions: crack depth and length. The three-dimensional complex variable function Taylor series expansion weighted function method is suitable for cracks with complex stress gradients, thus providing a certain reference for fracture mechanics parameter analysis within 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 the three-dimensional weighting function method, in order to improve the calculation efficiency of 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 the three-dimensional weighting function method includes the following steps:

[0009] S1: For structures with cracks, establish a finite element model and mesh it, and export the element and node information of the corresponding geometric model;

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

[0011] S3: Run the complex variable finite element program to calculate the nodal displacements U. r Read the imaginary part Im[U] of the solution of the normal displacement of the crack surface. r [(P';x,y)], where x and y represent the x-coordinate and y-coordinate of any point within the crack surface, respectively, and P' is the crack leading edge point;

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

[0013] S5: Obtain the hypothetical normal stress distribution σ(x,y) of the crack surface at the corresponding position of the crack surface in the cracked finite element model constructed in step S1 on the crack-free finite element model, where the load on the model is the actual working condition load.

[0014] S6: Combining the imaginary part Im[U] obtained in step S3 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) is used to calculate the stress intensity factor K(P') under the required load.

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

[0016] A crack model was created using the finite element software ABAQUS, and the mesh was generated and exported as an inp file.

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

[0018] MATLAB software was used to read the element and node information of the geometric model in S1. The nodes near the crack tip were perturbed along the crack opening direction to imbue their coordinates with complex variable information. The perturbation step size h' was orthogonally decomposed into three directions: x, y, and z. The original node coordinates were (x0, y0, z0), and the perturbed node coordinates were (x0 + y0, z0). i ,y0+iy i ,z0+iz i ), where i is the imaginary unit.

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

[0020] Run the finite element program to calculate the element stiffness matrix of all elements in the model in step S1, and then assemble the element stiffness matrices into a global stiffness matrix; introduce boundary conditions and loading conditions, with the loading condition being a reference load σ applied to the crack surface. r (x,y); Calculate the nodal displacements according to the equilibrium equation P=KU, and read the imaginary part Im[U] of the solution for the normal displacement of the crack surface. r (P';x,y)].

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

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

[0023]

[0024] Using a Taylor series expansion of the weight function for complex functions, and employing substitution, let x = racosθ, y = rcsinθ, where r is the polar radius in polar coordinates, θ is the polar angle in polar coordinates, and a and c are the crack depth and crack half-length, respectively. Then we have:

[0025]

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

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

[0028]

[0029] In the formula u i ,v j Let A be the integration points along the radial and polar directions, respectively, within the integration limit [-1, +1]. Each integration point has a corresponding integration weight A. i B j The stress σ at the integration point of the reference load r,ij 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)] located on the crack surface ((u i +1) / 2, π(v) j +1) / 4);

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

[0031]

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

[0033] A crack-free model was established using the finite element software ABAQUS, and the normal stress distribution σ(x,y) on the crack surface at the illusory location of the crack surface under the desired load was obtained by the finite element method.

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

[0035] Combined with the imaginary part Im[U] obtained in step S3 r [P';x,y)], the reference stress intensity factor K obtained in step S4 r The stress distribution σ(x,y) of the crack surface obtained in step S5 is expanded using the Taylor series of complex variable functions and the weight function is expanded. At the same time, substitution is used to let x = ra cosθ and y = rc sinθ. Substitute these values ​​into equation (4) to calculate the stress intensity factor K(P') of the point under the load.

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

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

[0038] 2. This invention is based on a two-dimensional weight function and extends it to three-dimensional elliptical cracks, thus broadening its application range.

[0039] 3. The crack stress intensity factor calculation method based on the three-dimensional weight function method of this invention is more accurate in solving the stress intensity factor under complex loads than the finite element method, and the calculation method is simple and effective. Attached Figure Description

[0040] Figure 1 This is a flowchart of the method of the present invention;

[0041] Figure 2 This is a schematic diagram of the geometric parameters of the corner crack in the pipe in an embodiment of the present invention;

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

[0043] Figure 4 This is a schematic diagram of the perturbation node of the corner crack mesh in an embodiment of the present invention;

[0044] Figure 5 This embodiment of the invention compares the results of calculating the stress intensity factor using the method of this embodiment and the finite element method under internal pressure load on a corner crack in a pipe. Among them, 5a, 5b and 5c are the comparison of the calculation results at crack leading edge points A, B and C, respectively. Detailed Implementation

[0045] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0046] Conversely, this invention encompasses any substitutions, modifications, equivalent methods, and solutions made within the spirit and scope of the invention as defined in the claims. Furthermore, to provide a better understanding of the invention, certain specific details are described in detail below. However, those skilled in the art will fully understand the invention even without these detailed descriptions.

[0047] Please see Figures 1-5 Taking a corner structure with an 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] A finite element model of the crack was established for the corner structure of the nozzle containing the inner corner crack. The structural geometric parameters are as follows: Figure 1 As shown, R i=500mm, other geometric parameters were calculated based on the dimensionless parameters in Table 1. Elastic modulus E = 204 GPa, Poisson's ratio μ = 0.3. The crack surface is subjected to a uniformly distributed stress of 1 MPa. The crack geometric model was established using ABAQUS finite element software.

[0049] Table 1 Geometric Parameters

[0050]

[0051] Perturbations were applied to point A on the surface of the cracked cylinder, the deepest point B, and point C on the surface of the nozzle, respectively. Figure 3 As shown. Figure 4 As shown, w is the crack width, and h' is taken as 10 times the crack depth a. -10 Run the complex variable finite element method code. Obtain the imaginary part of the displacement at the integration point. Calculate the solution K of the reference stress intensity factor 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 using finite element method, and the normal stress at the integration point was extracted. The stress intensity factor at point A on the cylinder surface at the crack leading edge, point B at the deepest point, and point C on the nozzle surface were calculated using Gaussian integration when the nozzle corner crack was subjected to an internal pressure of 1 MPa.

[0053] The stress intensity factor at the crack tip was calculated using both the method of this invention and the traditional finite element method based on J-integral. With the same model mesh density, the results obtained by the two methods are as follows: Figure 5 As shown, the method of the present invention has 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 within the protection scope of the present invention.

Claims

1. A method for calculating crack stress intensity factor based on the three-dimensional weighted function method, characterized in that, Includes the following steps: S1: For structures with cracks, establish a finite element model and mesh it, and export the element and node information of the corresponding geometric model; S2: Read the element and node information of the geometric model in step S1, and perturb the elements near the solution point required at the crack tip along the crack opening direction; S3: Run the complex variable finite element program to calculate the nodal displacements U. r Read the imaginary part Im[U] of the solution of the normal displacement of the crack surface. r [(P';x,y)], where x and y represent the x-coordinate and y-coordinate of any point within the crack surface, respectively, and P' is the crack leading edge point; S4: The imaginary part Im[U] of the solution for the normal displacement of the crack surface obtained in step S3. r [P';x,y)], and the normal stress σ on the crack surface under the reference load. r The reference stress intensity factor K is obtained by calculation from (x,y). r (P'), the reference load is a uniformly distributed stress of 1 MPa on the crack surface; S5: Establish a crack-free finite element model and obtain the hypothetical normal stress distribution σ(x,y) at the location corresponding to the crack surface in the cracked finite element model constructed in step S1 on the crack-free finite element model. The load on the model is the actual working condition load. S6: Combining the imaginary part Im[U] obtained in step S3 r [P';x,y)], the reference stress intensity factor K obtained in step S4 r Using the normal stress distribution σ(x,y) on the crack surface obtained in step S5 (P'), the stress intensity factor K(P') under the required load is calculated.

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

3. The method for calculating crack stress intensity factor based on the three-dimensional weighting function method according to claim 1, characterized in that, The specific process of step S3 is as follows: The element stiffness matrix of all elements in the model in step S1 is calculated using finite element analysis software, and then the element stiffness matrices are assembled into a global stiffness matrix. Boundary conditions and loading conditions are introduced, with the loading condition being a reference load σ applied to the crack surface. r (x,y); Calculate the nodal displacements according to the equilibrium equation P=KU, and read the imaginary part Im[U] of the solution for the normal displacement of the crack surface. r (P';x,y)].

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

1. Theoretical formula for three-dimensional crack weight function: ； Using a Taylor series expansion of the weight function for complex functions, and employing substitution, let x = racosθ, y = rcsinθ, where r is the polar radius in polar coordinates, θ is the polar angle in polar coordinates, and a and c are the crack depth and crack half-length, respectively. Then we have: ； The perturbation value h = w∙h', where w is the width of the perturbation unit; S4.2, Calculated using Gaussian integral, then: ； In the formula u i ,v j Let A be the integration points along the radial and polar directions, respectively, within the integration limit [-1, +1]. Each integration point has a corresponding integration weight A. i B j The stress σ at the integration point of the reference load r,ij 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); S4.

3. Based on the principle of self-consistency, the reference stress intensity factor K is obtained. r (P') equals: 。 5. The method for calculating crack stress intensity factor based on the three-dimensional weighting function method according to claim 1, characterized in that, The specific process of step S5 is as follows: A crack-free model was established using finite element analysis software, and the normal stress distribution σ(x,y) on the crack surface at the illusory location of the crack surface under the desired load was obtained by the finite element method.

6. The method for calculating crack stress intensity factor based on the three-dimensional weighting function method according to claim 4, characterized in that, The specific process of step S6 is as follows: Combined with the imaginary part Im[U] obtained in step S3 r [P';x,y)], the reference stress intensity factor K obtained in step S4 r The stress distribution σ(x,y) of the crack surface obtained in step S5 is expanded using the Taylor series of complex variable functions and the weight function is expanded. At the same time, substitution is used to let x=racosθ and y=rcsinθ. Substitute into equation (4) to calculate the stress intensity factor K(P') of the point under the load.

Images