CN1329689C - 压力容器疲劳寿命安全预测方法 - Google Patents

Abstract

本发明公开了一种压力容器疲劳寿命安全预测方法，包括通过输入初始缺陷信息、节点三维地址信息获得裂纹面积增量，进而获得裂纹扩展增量，将上述裂纹扩展增量与原有裂纹深度相加，即可获得裂纹深度。将上述裂纹深度与存储于计算机内部的容限尺寸进行比较，进而控制压力容器。本发明通过计算机对压力容器进行实时控制，对于预防在役压力容器的突发性破坏事故的发生、预测其疲劳剩余寿命、保证压力容器的安全使用具有重大的应用价值。

Description

压力容器疲劳寿命安全预测方法

技术领域

本发明涉及一种压力容器的安全测试方法，特别是压力容器疲劳寿命安全预测方法。

背景技术

压力容器是国民经济建设中的关键设备，同时又是一种具有***危险的特种承压设备，它承受着高温、低温、易燃、易爆、剧毒或腐蚀介质的压力，一旦发生破坏和泄漏将导致不可挽回的灾难性事故。在压力容器破坏事故中，有很大一部分是由于裂纹疲劳扩展引起的，因而能预先得知裂纹的扩展量，就能对压力容器采取相应的措施，避免***等事故的发生。压力容器的疲劳寿命预测一直是国内外断裂力学界研究的重点课题。目前普遍采用的方法是实测法，但该方法受各种客观条件的限制，无法对在役压力容器进行现场测试，且人力物力消耗大、有效性不能保证。一种基于应力强度因子幅ΔK的压力容器疲劳寿命预测的数值模拟技术目前也有文献报导，并已在实际工程中得到一定应用。但ΔK仅局限于线弹性和小范围屈服的断裂问题，对于工程实际中存在的大量的大范围屈服和全面屈服断裂问题已不能适用，宜采用理论上较为完善的循环J积分ΔJ作为断裂参量。

实践中，人们希望通过计算机对压力容器进行不间断地检测，即时作出关闭或保持使用等操作指令，再通过指令来控制压力容器，保证容器的绝对安全。

发明内容

针对现有技术中存在的技术问题，本发明提供一种通过计算机来控制压力容器，能够预防在役压力容器的突发性破坏事故的发生、预测其疲劳剩余寿命、保证压力容器安全使用的压力容器疲劳寿命安全预测方法。

本发明为达到以上目的，是通过这样的技术方案来实现的：提供一种压力容器疲劳寿命安全预测方法，包括以下步骤：

1)、用于输入初始缺陷信息的输入步骤；

2)、用于输入节点三维地址信息的输入步骤；

3)、用于获得裂纹面积增量的步骤，采用

$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_n、Δy_n、Δz_n表示节点n在X、Y、Z方向上的位移，h_n(s_i)表示形函数，J(s_i)表示雅可比矩阵，节点1、2、3处的s_i表示值分别为-1、0、+1。(见附图1)

4)、用于获得弹塑性循环J积分ΔJ的步骤，采用

$\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}$

其中A_n表示裂纹虚扩展导致的裂纹面积的增量，σ_ij表示应力张量，u_i表示位移矢量，W表示应变能密度，δ表示Kronecker张量，Δx_i表示两个独立坐标系间坐标变换的映射函数，V表示裂纹体中受裂纹虚扩展影响的体积区域。

5)、用于获得裂纹扩展增量的步骤，采用

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

其中N_i、N_i+1、a_i、a_i+1表示裂纹扩展到i、i+1时次数和裂纹长度，c、n表示材料常数，m表示计算分段数。

6)、用于获得裂纹深度的步骤，所述裂纹深度为原有裂纹深度与上述裂纹扩展增量之和；

7)、将上述裂纹深度与存储于计算机内部的容限尺寸进行比较的步骤，若裂纹深度小于容限尺寸，则输出信号，压力容器保持使用；若裂纹深度大于容限尺寸，则输出信号，关闭压力容器。

裂纹扩展增量是以弹塑性循环J积分ΔJ为控制信号，通过计算机转换而得，从而控制压力容器的使用。裂纹面积增量可简化为

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

本发明通过计算机对压力容器进行实时控制(通过计算机对压力容器进行不间断地检测，即时作出关闭或保持使用等操作指令，再通过指令来控制压力容器)，对于预防在役压力容器的突发性破坏事故的发生、预测其疲劳剩余寿命、保证压力容器的安全使用具有重大的应用价值。本发明以三维弹塑性循环J积分ΔJ作为断裂参量、采用Paris扩展模型、利用有限元技术模拟跟踪裂纹疲劳扩展全过程并预测其疲劳剩余寿命。这是一种基于循环J积分ΔJ的压力容器疲劳寿命预测的数值模拟新技术，该技术具有物理意义明确、直观性强、计算精度高、与实际吻合程度好等特点。

附图说明

图1是20节点三维等参元局部网格平面示意图；

图2是裂纹扩展示意图。

具体实施方式

首先根据压力容器的结构、材料和受力状况获得此压力容易的的容限尺寸，将此容限尺寸存入计算机数据库，待用。

通过无损检测方法(例如射线或超声波检测)定性、定量获得压力容器的缺陷类型、位置及尺寸，将此类初始缺陷信息输入计算机；将关于节点三维地址信息也输入计算机。

从Paris-Erdogan扩展模型可知，对于裂纹稳定扩展阶段可用下式表示：

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

对于工程实际中存在的大量的大范围屈服和全面屈服断裂问题的裂纹疲劳扩展速率，用ΔJ来代替ΔK，因此(1)式变为：

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

其中基于张量表示法的三维ΔJ的表达式如下：

$\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}---\left(3\right)$

式中：A_n为裂纹虚扩展导致的裂纹面积的增量；σ_ij为应力张量；u_i为位移矢量；w为应变能密度；δ为Kronecker张量；Δx_i为两个独立坐标系间坐标变换的映射函数；v为裂纹体中受裂纹虚扩展影响的体积区域

适合于有限元计算的矩阵形式为：

$\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)$

式中：N为受节点虚位移影响的单元；{J}为雅可比矩阵；α为相应的加权系数；

$\left\{\frac{\∂\mathrm{\ΔU}}{\∂X}\right\}=\left|\begin{array}{ccc}\frac{\∂\mathrm{\Δu}}{\∂x}& \frac{\∂\mathrm{\Δu}}{\∂y}& \frac{\∂\mathrm{\Δu}}{\∂z}\\ \frac{\∂\mathrm{\Δv}}{\∂x}& \frac{\∂\mathrm{\Δv}}{\∂y}& \frac{\∂\mathrm{\Δv}}{\∂z}\\ \frac{\∂\mathrm{\Δw}}{\∂x}& \frac{\∂\mathrm{\Δw}}{\∂y}& \frac{\∂\mathrm{\Δw}}{\∂z}\end{array}\right|=\left\{\mathrm{\Δ}{U}_0\right\}\left\{P\right\}{\left\{J\right\}}^{-1};$ $\left\{\frac{\∂\mathrm{\ΔX}}{\∂X}\right\}=\left|\begin{array}{ccc}\frac{\∂\mathrm{\Δx}}{\∂x}& \frac{\∂\mathrm{\Δx}}{\∂y}& \frac{\∂\mathrm{\Δx}}{\∂z}\\ \frac{\∂\mathrm{\Δy}}{\∂x}& \frac{\∂\mathrm{\Δy}}{\∂y}& \frac{\∂\mathrm{\Δy}}{\∂z}\\ \frac{\∂\mathrm{\Δz}}{\∂x}& \frac{\∂\mathrm{\Δz}}{\∂y}& \frac{\∂\mathrm{\Δz}}{\∂z}\end{array}\right|=\left\{\mathrm{\Δ}{X}_0\right\}\left\{P\right\}{\left\{J\right\}}^{-1};$

其余同式(3)。

裂纹面积增量A_n的计算非常关键，直接影响到J积分的计算精度，其具体计算介绍如下：

对于20节点等参元，由于形函数具有二次性，节点位移会引起如图1所示的网格变形现象。

根据有限元理论，图1中裂纹前沿某一点的坐标可以表示为：

$\left|\begin{array}{c}x\\ y\\ z\end{array}\right|=\left|\begin{array}{ccc}{x}_1& x_2& x_3\\ y_1& y_2& y_3\\ z_1& z_2& z_3\end{array}\right|\·\left|\begin{array}{c}{h}_1\left(s\right)\\ h_2\left(s\right)\\ h_3\left(s\right)\end{array}\right|---\left(5\right)$

相应的形函数为：

$\left|\begin{array}{c}{h}_1\left(s\right)\\ h_2\left(s\right)\\ h_3\left(s\right)\end{array}\right|=\left|\begin{array}{c}\frac{1}{2}s(s-1)\\ 1-s^2\\ \frac{1}{2}s(s+1)\end{array}\right|---\left(6\right)$

式中：节点1、2、3处的s值分别为-1、0、1。当节点3移至新的坐标 时：

$\left|\begin{array}{c}{\stackrel{\‾}{x}}_3\\ {\stackrel{\‾}{y}}_3\\ {\stackrel{\‾}{z}}_3\end{array}\right|=\left|\begin{array}{c}{x}_3\\ y_3\\ z_3\end{array}\right|+\left|\begin{array}{c}{\mathrm{\Δx}}_3\\ {\mathrm{\Δy}}_3\\ {\mathrm{\Δz}}_3\end{array}\right|---\left(7\right)$

由节点1、2及节点3新坐标确定的曲线上任一点的坐标可以表示为：

$\left|\begin{array}{c}\stackrel{\‾}{x}\\ \stackrel{\‾}{y}\\ \stackrel{\‾}{z}\end{array}\right|=\left|\begin{array}{ccc}{x}_1& x_2& x_3\\ y_1& y_2& y_3\\ z_1& z_2& z_3\end{array}\right|\·\left|\begin{array}{c}{h}_1\left(s\right)\\ h_2\left(s\right)\\ h_3\left(s\right)\end{array}\right|+h_3\left(s\right)\left|\begin{array}{c}{\mathrm{\Δx}}_3\\ {\mathrm{\Δy}}_3\\ {\mathrm{\Δz}}_3\end{array}\right|---\left(8\right)$

于是得到了新曲线上任一点的虚位移表达式为：

$\left\{\mathrm{\ΔX}\right\}=\left|\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ \mathrm{\Δz}\end{array}\right|=h_3\left(s\right)\left|\begin{array}{c}{\mathrm{\Δx}}_3\\ {\mathrm{\Δy}}_3\\ {\mathrm{\Δz}}_3\end{array}\right|---\left(9\right)$

因此求出图1中两曲线间阴影部分的裂纹面积增量为：

$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)$

式(10)中：dl为弧长增量；

${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}$

当然也可以求出由节点1或节点2处虚位移引起的裂纹面积增量，其计算公式统一表示为：

$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)$

对式(11)进行高斯积分，得到：

$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)$

式(12)中：I为高斯积分点。

通过上述推导，得到了由1个节点虚位移引起的在1个单元上裂纹面积增量的计算公式。需要指出的是，当由1个单元节点产生虚位移时，如图1所示，会引起旁边2个单元网格的变化，故A_n应为2个单元的裂纹面积增量。

在实际计算过程中，需要运用一些处理技巧。如节点2位于节点1、3联线中心时，A_n的计算就比较简单。由于x₁+x₃-2x₂＝0，y₁+y₃-2y₂＝0，z₁+z₃-2z₂＝0，则雅可比表达式J就可以简化为：

$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)$

式(13)中：L表示节点1、3之间连线长度。

相应地裂纹面积增量A_n的计算公式可简化为：

$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)$

另外J积分的路径应尽量避免裂纹前沿第1层单元因裂尖应力奇异性造成J积分偏大的情况，同时也应考虑减少计算量，而选取沿裂尖第2层单元的路线。

根据上述方法，在计算机上利用有限元技术就可以计算出裂纹扩展到每个阶段的裂纹前沿各点的循环J积分ΔJ值。

针对某一特定的裂纹，如图2所示，如果裂纹前沿曲线用几个离散点定义，则取裂纹上一点i的变化过程作为研究对象，当裂纹扩展后，i点成为裂纹新前沿面上的i+1点，此时，对于研究点i，可将式(2)两边积分转换为：

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

式中：N_i、N_i+1、a_i、a_i+1为裂纹扩展到i、i+1时次数和裂纹长度

实际工程中，我们可假定裂纹前沿曲线是预知的任意形状，裂纹前沿的三维弹塑性循环J积分ΔJ值可由有限元计算得到，则用Eular积分方程代替(15)式积分，可得：

$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)$

只要控制Δa在一定范围内，裂纹从a_i扩展到a_i+1时，可用数值方法得到规定精度范围内的解N_i+1，循环利用式(16)，便可计算出该点的整个扩展过程。

对整个裂纹前沿来说，假定裂纹初始前沿的J积分已求得，则对于j、m两点显然可运用Paris-Erdogan模型：

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

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

式(17)和式(18)中，ΔN|_j＝ΔN|_m，两式相除可得：

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

据式(19)可知，只要求出裂纹前沿上各点的循环J积分ΔJ值和裂纹上某一点的扩展增量，便可求出裂纹前沿各点的裂纹扩展增量。具体实施过程中可首先计算出裂纹最深点作为参考来计算其它各点，并定义该点裂纹增量为Δa_max，由此便可计算和模拟裂纹形貌在整个疲劳过程中的扩展变化情况，并由此可以得到裂纹的疲劳剩余寿命。

将前一循环所得裂纹深度加上本次循环所得的裂纹扩展增量，即为本次循环可得的裂纹深度(此裂纹深度为下一工作时的预测裂纹深度值)。当若裂纹深度小于数据库内的容限尺寸，则输出信号，压力容器保持使用；若裂纹扩展深度大于容限尺寸，则输出信号，关闭压力容器。

以上所述仅为本发明的一个具体实施方式，应当指出，对于本领域的普通技术人员来说，还可以作出许多变型和改进，所有的变型或改进均应视为本发明的保护范围。

Claims (3)

1、一种压力容器疲劳寿命安全预测方法，其特征是包括以下步骤：

1)、用于输入初始缺陷信息的输入步骤

2)、用于输入节点三维地址信息的输入步骤；；

3)、用于获得裂纹面积增量的步骤，采用

$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=1,\mathrm{2,3};$

4)、用于获得弹塑性循环J积分ΔJ的步骤，采用

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

5)、用于获得裂纹扩展增量的步骤，采用

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

6)、用于获得裂纹深度的步骤，所述裂纹深度为原有裂纹深度与上述裂纹扩展增量之和；

7)、将上述裂纹深度与存储于计算机内部的容限尺寸进行比较的步骤，若裂纹深度小于容限尺寸，则输出信号，压力容器保持使用；若裂纹深度大于容限尺寸，则输出信号，关闭压力容器。

2、根据权利要求1所述的压力容器疲劳寿命安全预测方法，其特征是：所述裂纹扩展增量是以弹塑性循环J积分ΔJ为控制信号，通过计算机转换而得。

3、根据权利要求2所述的压力容器疲劳寿命安全预测方法，其特征是：所述裂纹面积增量可简化为 $A_n=\frac{1}{6}L\sqrt{\mathrm{\Δ}{x}_n^2+\mathrm{\Δ}{y}_n^2+\mathrm{\Δ}{z}_n^2}n=\mathrm{1,3}\·$