CN1614294A - Fatigue life safety predicting method for pressure container - Google Patents

Abstract

A method for predicting the endurance life of pressure container includes such steps as inputting initial defect information and the 3D address information of node, calculating the area increment and expansion increment of crack, adding said expansion increment to original crack depth to obtain current crack depth, comparing it with the threshold stored in computer, and controlling the pressure container.

Description

Fatigue life safety predicting method for pressure container

Technical field

The present invention relates to a kind of safety detecting method of pressurized container, particularly fatigue life safety predicting method for pressure container.

Background technique

Pressurized container is the key equipment in the development of the national economy, be again a kind of extraordinary bearing device simultaneously with explosion risk, it is bearing the pressure of high temperature, low temperature, inflammable, explosive, severe toxicity or aggressive medium, will cause irremediable catastrophic failure in case take place to destroy and leak.In the pressurized container malicious event, have greatly because the tired expansion of crackle causes, thereby can learn the propagation of crackle in advance, just can take appropriate measures to pressurized container, avoid the generation of accidents such as exploding.The fatigue life prediction of pressurized container is the emphasis problem of domestic and international breaking power educational circles research always.The method that generally adopts is a measurement method at present, but this method is subjected to the restriction of various objective condition, can't carry out on-the-spot test to inservice pressure vessel, and manpower and materials consumption is big, validity can not guarantee.A kind of numerical simulation technology of the pressurized container fatigue life prediction based on stress intensity factor width of cloth Δ K also has reported in literature at present, and has necessarily been used in actual engineering.But Δ K only is confined to the breakage problem of linear elasticity and small scale yielding, can not be suitable for for a large amount of large scale yielding that exists in the engineering reality and complete yield breakage problem, should adopt in theory comparatively perfect cyclic J integration Δ J as fracture parameter.

In the practice, people wish by computer pressurized container to be carried out being detecting discontinuously, make immediately closing or keeping uses to wait operational order, and again by instructing the pilot pressure container, being perfectly safe of assurance container.

Summary of the invention

At the technical problem that exists in the prior art, the invention provides and a kind ofly come the pilot pressure container by computer, can prevent the sudden malicious event of inservice pressure vessel generation, predict its fatigue surplus life, guarantee the fatigue life safety predicting method for pressure container that pressure vessel safety uses.

The present invention is to realize by such technological scheme: a kind of fatigue life safety predicting method for pressure container is provided, may further comprise the steps for reaching above purpose:

1), is used to import the input step of initial imperfection information;

2), be used to import the input step of the three-dimensional address information of node;

3), be used to obtain the step of crackle area increment, adopt

$A_n=2\sqrt{\mathrm{\Δ}{x}_n^2+\mathrm{\Δ}{y}_n^2+\mathrm{\Δ}{z}_n^2}\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{h}_n\left(s_i\right)J\left(s_i\right)\mathrm{\α}\left(s_i\right)---n=\mathrm{1,2,3};$

Δ x wherein _n, Δ y _n, Δ z _nThe displacement of expression node n on X, Y, Z direction, h _n(s _i) the expression shape function, J (s _i) the expression Jacobian matrix, the s at node 1,2,3 places _iThe expression value is respectively-1,0 ,+1.

(seeing accompanying drawing 1)

4), be used to obtain the step of elastoplasticity cyclic J integration Δ J, adopt

$\mathrm{\ΔJ}=\frac{1}{A_n}\∫\∫{\∫}_V\left[(\mathrm{\Δ}{\mathrm{\σ}}_{\mathrm{ij}}\frac{\∂\mathrm{\Δ}{u}_j}{\∂x_k}-\mathrm{\Δw}{\mathrm{\δ}}_{\mathrm{ik}})\frac{\∂\mathrm{\Δ}{x}_k}{\∂x_i}\right]\mathrm{dV}$

A wherein _nThe increment of the crackle area that the empty expansion of expression crackle causes, σ _IjThe expression stress tensor, u ₁The expression displacement vector, W represents strain energy density, δ represents the Kronecker tensor, Δ x _iThe mapping function of representing transformation of coordinates between two independent coordinate systems, V represent to be subjected in the cracks in body the empty volumetric region that influences of expanding of crackle.

5), be used to obtain the step of crack propagation increment, adopt

$N_{i+1}-N_i=\frac{\mathrm{\Δa}}{c{\left(\mathrm{\Δ}{J}_i\right)}^n}---\mathrm{\Δa}=(a_{i+1}-a_i)/m$

N wherein _i, N _I+1, a _i, a _I+1Number of times and crack length when the expression crack propagation arrives i, i+1, c, n represent material constant, m represents to calculate segments.

6), be used to obtain the step of crack depth, described crack depth is the original crack degree of depth and above-mentioned crack propagation increment sum;

7), step that above-mentioned crack depth and the tolerance limit size that is stored in computer-internal are compared, if crack depth is less than the tolerance limit size, output signal then, pressurized container keeps using; If crack depth is greater than the tolerance limit size, output signal then, closing pressure container.

The crack propagation increment is to be control signal with elastoplasticity cyclic J integration Δ J, gets by the computer conversion, thus the use of pilot pressure container.Crackle area increment can be reduced to

$A_n=\frac{1}{6}L\sqrt{\mathrm{\Δ}{x}_n^2+\mathrm{\Δ}{y}_n^2+\mathrm{\Δ}{z}_n^2}---n=\mathrm{1,3}.$

The present invention controls in real time pressurized container by computer and (by computer pressurized container is detected incessantly, immediately make and close or keep operational orders such as use, come the pilot pressure container by instruction again), for the generation of the sudden malicious event of prevention inservice pressure vessel, predict its fatigue surplus life, guarantee that the safe handling of pressurized container has great application value.The present invention with D elastic-plastic cyclic J integration Δ J as fracture parameter, adopt the Paris extended model, utilize finite element technique simulation tracing crackle tired expansion overall process and predict its fatigue surplus life.This is a kind of numerical simulation new technology of the pressurized container fatigue life prediction based on cyclic J integration Δ J, this technology has explicit physical meaning, intuitive is strong, calculation accuracy is high, with characteristics such as actual degree of agreement is good.

Description of drawings

Fig. 1 is the local grid floor map of the three-dimensional isoparametric element of 20 nodes;

Fig. 2 is the crack propagation schematic representation.

Embodiment

At first structure, material and the force-bearing situation according to pressurized container obtains the tolerance limit size that this pressure is easy to, and deposits this tolerance limit size in Computer Database, and be stand-by.

Qualitative by nondestructive detecting method (for example ray or ultrasound examination), quantitatively obtain defect type, position and the size of pressurized container, with this type of initial imperfection information input computer; To also import computer about the three-dimensional address information of node.

From the Paris-Erdogan extended model as can be known, can be represented by the formula for the Stable Crack Growth stage:

$\frac{\mathrm{da}}{\mathrm{dN}}=c{\left(\mathrm{\ΔK}\right)}^n---\left(1\right)$

For the tired spreading rate of the crackle of a large amount of large scale yielding that exists in the engineering reality and complete yield breakage problem, replace Δ K with Δ J, therefore (1) formula becomes:

$\frac{\mathrm{da}}{\mathrm{dN}}=c{\left(\mathrm{\ΔJ}\right)}^n---\left(2\right)$

Wherein the representation based on the three-dimensional Δ J of tensor representation method is as follows:

$\mathrm{\ΔJ}=\frac{1}{A_n}\∫\∫{\∫}_V\left[(\mathrm{\Δ}{\mathrm{\σ}}_{\mathrm{ij}}\frac{\∂\mathrm{\Δ}{u}_j}{\∂x_k}-\mathrm{\Δw}{\mathrm{\δ}}_{\mathrm{ik}})\frac{\∂\mathrm{\Δ}{x}_k}{\∂x_i}\right]\mathrm{dV}---\left(3\right)$

In the formula: A _nBe the empty increment of expanding the crackle area that causes of crackle; σ _IjBe stress tensor; u _iBe displacement vector; W is a strain energy density; δ is the Kronecker tensor; Δ x _iIt is the mapping function of transformation of coordinates between two independent coordinate systems; V is the volumetric region that influenced by the empty expansion of crackle

The matrix form that is suitable for FEM (finite element) calculation is:

$\mathrm{\ΔJ}=\frac{1}{A_n}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\underset{i=1}{\overset{3}{\mathrm{\Σ}}}\underset{j=1}{\overset{3}{\mathrm{\Σ}}}\underset{k=1}{\overset{3}{\mathrm{\Σ}}}\left\{\mathrm{trace}\right[\left(\right\{\mathrm{\Δ\σ}\left\}\right\{\frac{\∂\mathrm{\ΔU}}{\∂X}\}-\mathrm{\Δw}\{I\left\}\right)\frac{\∂\mathrm{\ΔX}}{\∂X}\left]{\mathrm{\α}}_i{\mathrm{\α}}_j{\mathrm{\α}}_k\right|J\left|\right\}---\left(4\right)$

In the formula: N is the unit that influenced by the node virtual displacement; { J} is a Jacobian matrix; α is corresponding Weighting factor;

All the other cotypes (3).

Crackle area increment A _nCalculating very crucial, directly have influence on J integral Calculation precision, its concrete calculating is described below:

For 20 node isoparametric elements, because shape function has secondary, nodal displacement can cause distortion of the mesh phenomenon as shown in Figure 1.

According to finite element theory, any coordinate of crackle forward position can be expressed as among Fig. 1:

Corresponding shape function is:

In the formula: the s value at node 1,2,3 places is respectively-1,0,1.When node 3 moves to new coordinate (x ₃, y ₃, z ₃) time:

The coordinate of any point can be expressed as on the curve of being determined by node 1,2 and node 3 new coordinates:

So obtained the virtual displacement representation of any point on the new curve be:

Therefore obtain among Fig. 1 that the crackle area increment of dash area is between two curves:

$A_3=\∫\left|\mathrm{\ΔX}\right|\mathrm{dl}=\sqrt{\mathrm{\Δ}{x}_3^2+\mathrm{\Δ}{y}_3^2+\mathrm{\Δ}{z}_3^2}{\∫}_{-1}^{+1}{h}_3\left(s\right)\mathrm{Jds}---\left(10\right)$

In the formula (10): dl is the arc length increment;

$J=\left\{{[s(x_1+x_3-2x_2)+\frac{1}{2}(x_3-x_1)]}^2\right\}+[s(y_1+y_3-2y_2)+\frac{1}{2}(y_3-y_1){]}^2+[s(z_1+z_3-2z_2)+\frac{1}{2}(z_3-z_1){]}^2{\}}^{1/2}$

Can certainly obtain the crackle area increment that is caused by node 1 or node 2 place's virtual displacement, its formula unification is expressed as:

$A_n=\sqrt{\mathrm{\Δ}{x}_n^2+\mathrm{\Δ}{y}_n^2+\mathrm{\Δ}{z}_n^2}{\∫}_{-1}^{+1}{h}_n\left(s\right)\mathrm{Jds}----n=\mathrm{1,2,3}---\left(11\right)$

Formula (11) is carried out Gauss integration, obtains:

$A_n=\sqrt{\mathrm{\Δ}{x}_n^2+\mathrm{\Δ}{y}_n^2+\mathrm{\Δ}{z}_n^2}\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{h}_n\left(s_i\right)J\left(s_i\right)\mathrm{\α}\left(s_i\right)---n=\mathrm{1,2,3}---\left(12\right)$

In the formula (12): I is the Gauss integration point.

By above-mentioned derivation, obtained the formula of crackle area increment on 1 unit that causes by 1 node virtual displacement.It is pointed out that and when producing virtual displacement, as shown in Figure 1, can cause the variation of 2 unit grids in next door, so A by 1 cell node _nShould be the crackle area increment of 2 unit.

In the Practical Calculation process, need some disposal skill of utilization.When being positioned at node 1,3 on line centers as node 2, A _nCalculating just fairly simple.Because x ₁+ x ₃-2x ₂=0, y ₁+ y ₃-2y ₂=0, z ₁+ z ₃-2z ₂=0, then Jacobi's representation J just can be reduced to:

$J=\frac{1}{2}\sqrt{{(x_3-x_1)}^2+{(y_3-y_1)}^2+{(z_3-z_1)}^2}=\frac{1}{2}L---\left(13\right)$

In the formula (13): L represents wire length between the node 1,3.

Crackle area increment A correspondingly _nFormula can be reduced to:

$A_n=\frac{1}{6}L\sqrt{\mathrm{\Δ}{x}_n^2+\mathrm{\Δ}{y}_n^2+\mathrm{\Δ}{z}_n^2}---n=\mathrm{1,3}---\left(14\right)$

The path of J integration should avoid the 1st layer of unit, crackle forward position to cause J integration situation bigger than normal because of splitting sharp Stress singularity as far as possible in addition, also should consider to reduce amount of calculation simultaneously, and chooses along the route that splits the 2nd layer of unit of point.

According to said method, utilize finite element technique just can calculate the cyclic J integration Δ J value of crack propagation on computers to the crackle forward position each point in each stage.

At a certain specific crackle, as shown in Figure 2, if crackle forward position curve defines with several discrete points, the change procedure of then getting 1 i on the crackle is as research object, after crack propagation, the i point becomes the i+1 point on the new leading surface of crackle, at this moment, for research point i, can be with formula (2) both sides integrating conversion:

${\∫}_{a_j}^{a_{i+1}}\mathrm{da}={\∫}_{N_j}^{N_{i+1}}c{\left(\mathrm{\ΔJ}\right)}^n\mathrm{dN}---\left(15\right)$

In the formula: N _i, N _I+1, a _i, a _I+1Number of times and crack length when arriving i, i+1 for crack propagation

In the actual engineering, we can suppose that crackle forward position curve is the arbitrary shape of precognition, and the D elastic-plastic cyclic J integration Δ J value in crackle forward position can be obtained by FEM (finite element) calculation, then replaces (15) formula integration with the Eular integral equation, can get:

$N_{i+1}-N_i=\frac{\mathrm{\Δa}}{c{\left(\mathrm{\Δ}{J}_i\right)}^n}---\mathrm{\Δa}=(a_{i+1}-a_i)/m---\left(16\right)$

As long as control Δ a within the specific limits, crackle is from a _iExpand to a _I+1The time, available numerical method obtains the N that separates in the specified accuracy scope _I+1, recycling formula (16) just can calculate the whole expansion process of this point.

Concerning whole crackle forward position, suppose that the J integration in the initial forward position of crackle is tried to achieve, then obviously can use the Paris-Erdogan model at 2 for j, m:

Δa| _j＝c(ΔJ) ⁿΔN| _j (17)

Δa| _m＝c(ΔJ) ⁿΔN| _m (18)

In formula (17) and the formula (18), Δ N| _j=Δ N| _m, two formulas are divided by and can be got:

$\frac{\mathrm{\Δ}{a}_j}{\mathrm{\Δ}{a}_m}={\left(\frac{\mathrm{\Δ}{J}_j}{\mathrm{\Δ}{J}_m}\right)}^n----\left(19\right)$

According to formula (19) as can be known, as long as obtain on the cyclic J integration Δ J value of each point on the crackle forward position and the crackle certain any expansion increment, just can obtain the crack propagation increment of crackle forward position each point.Can at first calculate the crackle deepest point in the specific implementation process and calculate other each point as a reference, and to define this crackle increment be Δ a _Max, just can calculate thus and the expansion situation of change of simulating crack pattern in whole fatigue process, and can obtain the fatigue surplus life of crackle thus.

Last circulation gained crack depth added the crack propagation increment of this circulation gained be the crack depth (the prediction crack depth value when this crack depth is worked for next) that this circulation can get.When if crack depth less than the tolerance limit size in the database, output signal then, pressurized container keeps using; If the crack propagation degree of depth is greater than the tolerance limit size, output signal then, closing pressure container.

The above is a specific embodiment of the present invention only, should be pointed out that for the person of ordinary skill of the art, can also make many modification and improvement, and all modification or improvement all should be considered as protection scope of the present invention.

Claims (3)

1, a kind of fatigue life safety predicting method for pressure container is characterized in that may further comprise the steps:

1), is used to import the input step of initial imperfection information

2), be used to import the input step of the three-dimensional address information of node;

3), be used to obtain the step of crackle area increment, adopt

$A_n=2\sqrt{{\mathrm{\Δx}}_n^2+{\mathrm{\Δy}}_n^2+{\mathrm{\Δz}}_n^2}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{h}_n\left(s_i\right)J\left(s_i\right)\mathrm{\α}\left(s_i\right)---n=\mathrm{1,2,3};$

4), be used to obtain the step of elastoplasticity cyclic J integration Δ J, adopt

$\mathrm{\ΔJ}=\frac{1}{A_n}\∫\∫\∫\left[({\mathrm{\Δ\σ}}_{\mathrm{ij}}\frac{{\∂\mathrm{\Δu}}_j}{{\∂x}_k}-{\mathrm{\Δw\δ}}_{\mathrm{ik}})\frac{{\∂\mathrm{\Δx}}_k}{{\∂x}_i}\right]\mathrm{dV}$

5), be used to obtain the step of crack propagation increment, adopt

$N_{i+1}-N_i=\frac{\mathrm{\Δa}}{{c\left({\mathrm{\ΔJ}}_i\right)}^n}\mathrm{\Δa}=(a_{i+1}-a_i)/m$

6), be used to obtain the step of crack depth, described crack depth is the original crack degree of depth and above-mentioned crack propagation increment sum;

7), step that above-mentioned crack depth and the tolerance limit size that is stored in computer-internal are compared, if crack depth is less than the tolerance limit size, output signal then, pressurized container keeps using; If crack depth is greater than the tolerance limit size, output signal then, closing pressure container.

2, fatigue life safety predicting method for pressure container according to claim 1 is characterized in that: described crack propagation increment is to be control signal with elastoplasticity cyclic J integration Δ J, gets by the computer conversion.

3, fatigue life safety predicting method for pressure container according to claim 2 is characterized in that: described crackle area increment can be reduced to $A_n=\frac{1}{6}L\sqrt{{\mathrm{\Δx}}_n^2+{\mathrm{\Δy}}_n^2+{\mathrm{\Δz}}_n^2}---n=\mathrm{1,3}.$