CN101509855A - 一种高温管道损伤及寿命的在线预测方法 - Google Patents

Abstract

本发明是一种高温管道损伤及寿命的在线预测方法，其实施步骤包括：(1)对高温管道进行损伤耦合的有限元模拟分析；(2)根据分析结果，找出重点监测部位，设置传感器，监测其应变；(3)针对不同工况，进行有限元(包含本构方程的分析子程序)分析，建立损伤－应变及剩余寿命－应变的数据库；(4)将在线监测的应变值和载荷工况的数值与数据库的数据进行在线查询对比，得出相应损伤和剩余寿命的评估值；本发明的优点是：在正常生产的同时能对运行中的高温管道进行实时监测，及时反映重要部件、关键部位的变形和损伤，对管道的使用和剩余寿命做出正确的估算，有利于保障生产安全、调整生产负荷、合理维修规划并有效地延长生产设备的寿命。

Description

一种高温管道损伤及寿命的在线预测方法

【技术领域】

本发明涉及高温构件的寿命预测领域，涉及一种基于应变测量的寿命预测技术，可用于实时监测高温管道的损伤情况及剩余寿命。

【背景技术】

电力工业是关系国计民生的重要行业。火电机组是电厂的重要组成部分，它的使用寿命与电厂运行的经济性和安全性密切相关。高温高压管线尤其是主蒸汽管道的设计寿命是普遍关心的重要问题。长期工作在高温高压下的主蒸汽管道材料，在温度和应力作用下，材料的显微结构、力学性能、耐热性能、物理性能都会缓慢地变化。材料在使用过程中出现的老化现象，实质上是蠕变损伤积累和材质劣化，并最终在应力作用下开裂失效的现象。

对于高温高压管线寿命预测技术，特别是热电厂中的主蒸汽管道的寿命设计与预测，国内外的科研机构和生产单位开展了大量的研究工作，积累了一些有意义的研究成果。目前国内外预测主蒸汽管的使用寿命的研究基本都是取样进行实验分析，基于不同的模型或方法离线进行剩余寿命预测。

预测主蒸汽管道寿命的外推方法多种多样，早期的外推方法主要采用等温线外推法，这种方法精度低，试验时间较长。到了20世纪50年代，出现了外推精度较高的时间—温度参数外推法，其中以Larson-Miller公式最为著名。近年来，主蒸汽管道寿命预测方法已由持久强度为主要指标的传统方法转向以蠕变变形量为主要指标的θ法、由θ法发展而来的C射影法以及蠕变曲线逐步外推法。

另外，以持久强度试验为主的综合分析法一直被人们用来预测主蒸汽管道的寿命。基于可靠性的寿命评估方法和蠕变损伤及裂纹扩展的寿命评估方法也是近年来人们提出的预测主蒸汽管道寿命的新方法。

在国内，已有多项专利技术用于研究高温部件寿命预测的问题，如中国专利CN03134314.7“高温部件蠕变寿命的测试方法”、CN200610048003.8“高温构件材料的寿命及寿命消耗的预测方法”、CN200710013519.3“高温炉管剩余寿命评估方法及装置”、CN200710039899.8“一种汽轮机高温部件蠕变寿命的预测方法及系统”等多是利用温度、压力等参数间接地算出应力及应变，进而推算出寿命，其精度难于保证；或者是截取部分材料进行持久性试验数据分析，外推其服役寿命，不能用于在线监测。

【发明内容】

本发明的目的在于弥补现有技术的不足，提出一种直接基于应变测量的方法，可应用于在线监测系统的高温管道损伤及寿命的预测，实现对管道损伤及寿命的在线实时评价和预测。

为实现上述目的，本发明采取的技术方案为：

一种高温管道损伤及在线寿命的预测方法，其特征是，利用测量得到的应变，结合有限元分析，预测高温管道的损伤及在线寿命，其实施步骤包括：

(1)对高温管道进行损伤耦合的有限元模拟分析；

(2)根据有限元分析结果，找出重点监测部位，对重点监测部位设置传感器，监测其应变；

(3)针对不同载荷工况，进行有限元分析，建立损伤-应变及剩余寿命-应变的关系数据库；

(4)将在线监测的应变值和载荷工况的数值，与损伤-应变及剩余寿命-应变数据库的数据进行在线查询对比，得出相应损伤和剩余寿命的评估值。

在步骤(1)中，采用UG建立管道的几何模型，并导入有限元分析软件中，选用管单元分析管线的轴向位移、应变和应力，同时得到管道的周向应力；管单元在保证足够的精度要求的同时，可减少大量的计算时间(相对于常用的实体单元和壳单元)，从而为大规模管线应力的分析提供了实际的操作性。

步骤(3)中所述的有限元分析，包含以基于连续损伤力学的本构方程为准则编制的含损伤计算的有限元分析子程序，通过对高温管道进行大规模的有限元模拟计算，反演出应变—损伤关系和应变—剩余寿命的关系。

在步骤(4)中分析应变、温度、压力等现场数据，在损伤-应变和寿命-应变数据库中检索相应的应变所对应的损伤和剩余寿命，或者根据温度压力的变化分析损伤累积；利用WEB服务器响应终端用户的查询请求，并将应用服务器的处理结果生成相应的表单和参数曲线，通过Internet网页的形式发布；经过授权的用户可通过Internet浏览器实时查询管线的损伤状况和剩余寿命。

步骤(3)中所述的有限元分析中包含的本构方程为：以连续损伤力学理论的本构方程为准则，通过对其模拟求解出相应的应变和损伤结果；包含：采用单状态变量的本构方程：

$\frac{{\mathrm{d\ϵ}}_{\mathrm{ij}}^c}{\mathrm{dt}}=\frac{3}{2}B{\left[\frac{{\mathrm{\σ}}_e}{1-D}\right]}^n\frac{S_{\mathrm{ij}}}{{\mathrm{\σ}}_e}{t}^m---(1-1)$

式中，

为蠕变应变张量，s_ij为应力偏张量，σ₁为最大主应力，σ_e为VonMises应力，D为损伤变量，其值在0～1之间变化；

α(0≤α≤1)是材料常数，用来表征多轴破断准则，B、n、A、υ是与最小蠕变应变速率以及断裂相关的材料常数；

或由多相复合假设得出的本构模型为：

$\frac{{\mathrm{d\ϵ}}_{\mathrm{ij}}^c}{\mathrm{dt}}=\frac{3}{2}B{\mathrm{\σ}}_e^{n-1}{s}_{\mathrm{ij}}[(1-\mathrm{\ρ})+\mathrm{\ρ}{(1-D)}^{-n}]---(1-3)$

$\frac{\mathrm{dD}}{\mathrm{dt}}=g\frac{A}{\mathrm{\φ}+1}\frac{{[\mathrm{\α}{\mathrm{\σ}}_1+(1-\mathrm{\α}){\mathrm{\σ}}_e]}^{\mathrm{\υ}}}{{(1-D)}^{\mathrm{\φ}}}---(1-4)$

$D_{\mathrm{cr}}=1-{(1-g)}^{\frac{1}{\mathrm{\φ}+1}}---(1-5)$

式中，D_cr为临界损伤值，当 $\frac{D}{D_{\mathrm{cr}}}=1$ 时，表明材料达到了其蠕变寿命；

g、φ和ρ是指定非均质损伤的材料常数；

考虑到复杂晶相对蠕变变形和损伤演化的影响，可采用多状态变量本构方程，如两个状态变量的模型形式如下：

$\frac{d{\mathrm{\ϵ}}_{\mathrm{ij}}^c}{\mathrm{dt}}=\frac{3}{2}\frac{G{\mathrm{\σ}}_e^{n-1}{s}_{\mathrm{ij}}}{{(1-\mathrm{\Φ})}^n{(1-D)}^n}{t}^m---(1-6)$

$\frac{\mathrm{d\Φ}}{\mathrm{dt}}=\frac{k_c}{3}{(1-\mathrm{\Φ})}^4---(1-8)$

式中，第二个状态变量Φ描述了碳化物的沉积间隔，G和k_c是材料常数；将服从虎克定理的弹性应变和蠕变应变联立得：

$\frac{{\mathrm{d\ϵ}}_{\mathrm{ij}}}{\mathrm{dt}}=\frac{1+\mathrm{\μ}}{E}[\left(\frac{d{\mathrm{\σ}}_{\mathrm{ij}}}{\mathrm{dt}}\right)-\frac{\mathrm{\μ}}{1+\mathrm{\μ}}(\frac{d{\mathrm{\σ}}_{\mathrm{kk}}}{\mathrm{dt}}\left){\mathrm{\δ}}_{\mathrm{ij}}\right]+\frac{3}{2}B{\mathrm{\σ}}_e^{n-1}{s}_{\mathrm{ij}}[(1-\mathrm{\ρ})+\mathrm{\ρ}{(1-D)}^{-n}]---(1-9)$

首先，将有限元计算分为弹性域中的有限元和含蠕变应变的有限元分析两部分：

(1)弹性域中的有限元分析

有限元法首先将结构离散化成不同的单元，单元与单元间以节点相连；在弹性范围内，应力应变的本构关系可表示为：

{σ}＝[D]{ε}                           (1-10)

式中[D]为弹性矩阵；

若按直接法选定位移模式，则结构位移场与单元节点位移间的关系为：

{u}＝[N]{δ}                           (1-11)

式中{u}为位移场向量，[N]是坐标的函数，反映单元的位移形态，被称作“形函数”；δ为单元节点的位移向量；

又因为应变与位移存在以下关系：

{ε}＝[B]{δ}                            (1-12)

若作用在物体上的外力为P，则由虚功原理可得：

[K]{δ}＝{P}                            (1-13)

其中：

[K]＝∑[k]                             (1-14)

$\left[k\right]={\∫}_V{\left[B\right]}^{-1}\left[D\right]\left[B\right]\mathrm{dV}---(1-15)$

式中[k]为单元刚度矩阵，[K]为总刚度矩阵，[B]为应变矩阵，{P}为节点载荷向量；通过求解代数方程组(1-13)即可求得弹性应力解与应变解；

(2)含蠕变应变的有限元分析

若考虑带有蠕变变形的情况，则：

{ε^e}＝{ε}-{ε^c}                         (1-16)

式中上标e代表相应的弹性分量，c代表相应的蠕变分量；

将(1-16)分别代入(1-10)和(1-11)中，得：

{σ^e}＝[D]({ε}-{ε^c})                    (1-17)

{δ^e}＝[B]^-1({ε}-{ε^c})          (1-18)

这样(1-13)就可以表达成含有蠕变项的平衡方程：

{K]([B]^-1({ε}-{ε^c}))＝{P}      (1-19)；

其次，将非线性方程线性化：

在有限单元法中描述物体离散后的非线性特性，可以用下面的平衡方程式表示：

F^N(u^M)＝0                       (1-20)

式中，F^N是与第N个变量对应的分力，u^M表示第M个变量的值；

对于(1-20)式的非线性方程，假设第i步迭代后可获得的近似解为

令离散方程的精确解与近似解之间的差值为即：

$F^N(u_i^M+c_{i+1}^M)=0---(1-21)$

方程左端在近似解

处按泰勒级数展开得：

$F^N\left(u_i^M\right)+\frac{{\∂F}^N}{\∂u^P}\left(u_i^M\right)c_{i+1}^P+\frac{{\∂}^2{F}^N}{\∂u^P\∂u^Q}\left(u_i^M\right)c_{i+1}^P{c}_{i+1}^Q+.......=0---(1-22);$

忽略二阶以上的偏导数，得到如下的线性方程：

$K_i^{\mathrm{NP}}{c}_{i+1}^P=-F_i^N---(1-23);$

式(1-23)中的 $K_i^{\mathrm{NP}}=\frac{{\∂F}^N}{\∂u^P}\left(u_i^M\right)$ 即为雅可比矩阵，并且 $F_i^N=F^N\left(u_i^M\right);$ 这样第i+1个近似解可表示为：

$u_{i+1}^M=u_i^M+c_{i+1}^M---(1-24)$

上式即为线性化后的迭代公式；

材料的雅可比矩阵定义为：

$J=\frac{\∂\mathrm{\Δ}{\mathrm{\σ}}_{\mathrm{ij}}}{\∂\mathrm{\Δ}{\mathrm{\ϵ}}_{\mathrm{ij}}}---(1-25)$

由前所述，若忽略塑性后总的应变式为：

${\mathrm{\ϵ}}_{\mathrm{ij}}={\mathrm{\ϵ}}_{\mathrm{ij}}^e+{\mathrm{\ϵ}}_{\mathrm{ij}}^c---(1-26)$

程序中每一载荷的施加都是通过两个载荷步来完成的：在第一步，通过定义一个服从虎克定理的材料模型来产生一个弹性响应；第二步蠕变分析中，实现如前所述的蠕变与损伤耦合的材料行为；

由于第一步是与时间无关的弹性响应， ${\mathrm{\ϵ}}_{\mathrm{ij}}={\mathrm{\ϵ}}_{\mathrm{ij}}^e,$ 因此，其应力增量与雅克比矩阵都较为简单，由弹性力学可知，线弹性条件下应力与应变的关系为：

{σ}＝[D]{ε}                   (1-27)

式(1-27)中：

$\left[D\right]=\frac{E(1-\mathrm{\μ})}{(1+\mathrm{\μ})(1-2\mathrm{\μ})}\left[\begin{array}{cccccc}1& \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& 0& 0& 0\\ \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& 1& \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& 0& 0& 0\\ \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& 1& 0& 0& 0\\ 0& 0& 0& \frac{1-2\mathrm{\μ}}{2(1-\mathrm{\μ})}& 0& 0\\ 0& 0& 0& 0& \frac{1-2\mathrm{\μ}}{2(1-\mathrm{\μ})}& 0\\ 0& 0& 0& 0& 0& \frac{1-2\mathrm{\μ}}{2(1-\mathrm{\μ})}\end{array}\right]$

式中的[D]即为第一步中需计算的J；

在第二步中应力和应变都是与时间相关的，不同的时间所对应的J也不相同；由于蠕变损伤方程的表达式是非线性的微分方程，所以无法直接给出相应的雅可比矩阵，需要对蠕变损伤方程进行适当的处理；在子程序中引入一个显式的时间积分法，采用中心差分对目标方程进行离散化处理；采用的中心差分表达式和算符表示如下：

${\stackrel{\·}{f}}_{t+\frac{1}{2}\mathrm{\Δt}}=\frac{\mathrm{\Δf}}{\mathrm{\Δt}}---(1-28)$

$f_{t+\frac{1}{2}\mathrm{\Δt}}{=t}_t+\frac{\mathrm{\Δf}}{2}---(1-29)$

式中f为一任意函数，f_t表示该函数在增量开始时的值，Δf表示在该增量中函数的变化值，Δt是时间增量；

把蠕变损伤连续方程(式(1-9)和式(1-4))按中心差分离散化后得：

$\frac{{\mathrm{\Δ\ϵ}}_{\mathrm{ij}}}{\mathrm{\Δt}}=\frac{1+\mathrm{\μ}}{E}[\left(\frac{\mathrm{\Δ}{\mathrm{\σ}}_{\mathrm{ij}}}{\mathrm{\Δt}}\right)-\frac{\mathrm{\μ}}{1+\mathrm{\μ}}(\frac{\mathrm{\Δ}{\mathrm{\σ}}_{\mathrm{kk}}}{\mathrm{\Δt}}\left){\mathrm{\δ}}_{\mathrm{ij}}\right]+\frac{3}{2}B{({\mathrm{\σ}}_e+\frac{{\mathrm{\Δ\σ}}_e}{2})}^{n-1}(s_{\mathrm{ij}}+\frac{\mathrm{\Δ}{s}_{\mathrm{ij}}}{2})[(1-\mathrm{\ρ})+\mathrm{\ρ}{(1-D)}^{-n}]$

                                                      (1-30)

$\frac{\mathrm{\ΔD}}{\mathrm{\Δt}}=g\frac{A}{\mathrm{\φ}+1}\frac{{[\mathrm{\α}{\mathrm{\σ}}_1+(1-\mathrm{\α}){\mathrm{\σ}}_e]}^{\mathrm{\υ}}}{{[1-(D+\frac{\mathrm{\ΔD}}{2})]}^{\mathrm{\φ}}}---(1-31)$

由于蠕变损伤方程在中心差分后，所得的雅可比矩阵J是一个6×6阶的非对称矩阵；直接给出该方程式的应力增量对应变增量的偏导数存在难度，而应变增量却可以在非线性方程组里以闭合形式给出，所以，首先计算应变增量对应力增量的偏导数：

$J^{-1}=\frac{\∂\mathrm{\Δ}{\mathrm{\ϵ}}_{\mathrm{ij}}}{\∂\mathrm{\Δ}{\mathrm{\σ}}_{\mathrm{ij}}}---(1-32)$

然后通过求逆J^-1得到J；

将以上方法编制成UMAT用户子程序的步骤为：先通过判断模块(JUDGE_HTTD)判断应力状态，再分别进入弹性响应模块(STEP1)和蠕变损伤响应模块(STEP2)进行计算；然后分别转入弹性响应模块和蠕变损伤响应模块的子程序模块进行相应计算：弹性响应模块给出材料在弹性范围内的雅可比矩阵J和相应的应力应变值；蠕变损伤响应模块最后给出材料在增量步下的雅可比矩阵J；

蠕变损伤响应模块中的计算损伤增量模块(DAMSOLVE_HTTD)和其子程序一起计算，得到在当前应力状况下损伤增量(ΔD)的大小，并对方程(1-31)进行变形，得到外部函数计算损伤函数的函数值(DAMFUNC_HTTD)：

$f\left(\mathrm{\ΔD}\right)=g\frac{A}{\mathrm{\φ}+1}\frac{{[\mathrm{\α}{\mathrm{\σ}}_1+(1-\mathrm{\α}){\mathrm{\σ}}_e]}^{\mathrm{\υ}}}{{[1-(D+\frac{\mathrm{\ΔD}}{2})]}^{\mathrm{\φ}}}-\frac{\mathrm{\ΔD}}{\mathrm{\Δt}}---(1-33)$

由于D的值在[0，1]区间单调递增，所以ΔD的变化也不会超过该区间；令ΔD在0至1间以0.05的倍数分别代入(1-31)式进行试根，若f(ΔD₁)×f(ΔD₂)≤0，由零点定理可知公式(1-26)在[ΔD₁，ΔD₂]上有根存在；分别以ΔD₁和ΔD₂作为上下限，代入子程序计算损伤函数的根(ZBRENT_HTTD)，利用布伦特法(Brent’s method)保证所获得的ΔD的误差不超过一定容限；其中最大主应力σ₁是通过调用有限元软件(ABAQUS)内部提供的SPRINC子程序得到；

计算应力增量(NEWT_HTTD)模块包含了外部函数计算范数(FMIN_HTTD)和三个子程序模块：计算雅可比矩阵模块(FJAC_HTTD)、LU分解模块(LUDCMP_HTTD)、解线性方程组模块(LUBKSB_HTTD)；计算雅可比矩阵模块的功能是把损伤增量代入到公式(1-15)中，形成六个联立的非线性方程组，采用牛顿-拉斐森法搜索六个应力增量的大小；若六个应力增量符合如下不等式：

$\underset{j=1}{\overset{3}{\underset{i=1}{\mathrm{\Σ}}}}\{\frac{1+\mathrm{\μ}}{E}[\left(\frac{\mathrm{\Δ}{\mathrm{\σ}}_{\mathrm{ij}}}{\mathrm{\Δt}}\right)-\frac{\mathrm{\μ}}{1+\mathrm{\μ}}\left(\frac{\mathrm{\Δ}{\mathrm{\σ}}_{\mathrm{kk}}}{\mathrm{\Δt}}\right){\mathrm{\δ}}_{\mathrm{ij}}]$

$+\frac{3}{2}B{({\mathrm{\σ}}_e+\frac{\mathrm{\Δ}{\mathrm{\σ}}_e}{2})}^{n-1}{(s_{\mathrm{ij}}+\frac{{\mathrm{\Δs}}_{\mathrm{ij}}}{2})[(1-\mathrm{\ρ})+\mathrm{\ρ}{(1-D)}^{-n}]-\frac{\mathrm{\Δ}{\mathrm{\ϵ}}_{\mathrm{ij}}}{\mathrm{\Δt}}{\}}^2}_{\≤\mathrm{\ξ}}---(1-34)$

式(1-34)中，ξ为允许的误差，则认为此六个应力分量就是所求值。

本发明的技术原理是：通过对超临界发电厂高温管道进行大规模的有限元模拟计算，反演出损伤-应变关系和剩余寿命-应变关系；同时建立损伤-应变和寿命-应变的数据库，用于在线监测系统；应变传感器监测到管道的应变后，与数据库中的应变数据进行对比分析，实现对管道损伤情况和剩余寿命的在线评价和预测。

本发明的积极效果是：能在正常生产的同时，保证对运行中的高温管道系统进行实时监测，及时反映系统中重要部件、关键部位的应力和应变，对管道的使用寿命和剩余寿命做出正确的估算，有利于改进设计、调整生产负荷，合理有效地延长生产设备的寿命。

【附图说明】

附图1是UMAT用户子程序的流程图；

附图2是本发明一种高温管道损伤及在线寿命的预测方法的流程框图；

附图3是T型接管的几何尺寸图；

附图4是T型接管的网格模型图；

附图5是T型接管在25MPa内压下达到临界损伤时的损伤分布云纹图；

附图6是25MPa内压下T型接管最大损伤部位的应变发展曲线图；

附图7是25MPa内压下T型接管最大损伤部位的损伤发展曲线图；

附图8是25MPa内压下T型接管最大损伤部位的损伤-应变曲线图；

附图9是弯管的三维实体模型图；

附图10是弯管的网格模型图；

附图11是弯管在25MPa内压下达到临界损伤时的损伤分布云纹图；

附图12是25MPa内压下弯管最大损伤部位的应变发展曲线图；

附图13是25MPa内压下弯管最大损伤部位的损伤发展曲线图；

附图14是25MPa内压下弯管最大损伤部位的损伤-应变曲线图。

【具体实施方式】

以下结合附图和实施例进一步解释本发明一种高温管道损伤及在线寿命的预测方法，但是，本发明的实施不限于以下的形式。

附图1为UMAT用户子程序的流程图，首先通过判断模块(JUDGE_HTTD)判断应力状态，再分别进入弹性响应(STEP1)和蠕变损伤响应(STEP2)进行计算，其中，判断模块首先判断单元上承受的载荷是否为静水压力；因为当构件在承受静水压力时，材料不会发生蠕变行为，而且一旦单元上承受的是静水压力时，其应力偏张量就将是一个零矩阵，这会导致UMAT在后续的蠕变损伤计算中产生奇异；因此，每次调用UMAT时，程序第一个进入的就是该子程序模块来检测是否有静水压力的存在，一旦有则只对其弹性响应计算，忽略蠕变行为的计算；该模块还将根据载荷子步的编号数(KSTEP值)判断该载荷步是属于第一个弹性响应，还是后续的蠕变损伤响应，然后分别转入各自的子程序模块进行相应计算；弹性响应模块给出材料在弹性范围内的雅可比矩阵J，以及相应的应力应变值；蠕变损伤响应模块通过调用计算损伤增量模块(DAMSOLVE_HTTD)得到每一增量步下的损伤增量，通过调用计算应力增量(NEWT_HTTD)得到相应的应力增量，最后给出材料在该增量步下的雅可比矩阵J。

参见附图2，本发明一种高温管道损伤及在线寿命的预测方法的流程框图。

首先对高温管道进行有限元模拟计算，一方面，采用UG建立管道的几何模型，并导入有限元分析软件中，选用单元分析管线的轴向位移、应变和应力，同时得到管道的周向应力(沿壁厚方向相等，沿周向则可以根据不同的要求分布相应数目的积分点)。管单元在保证足够的精度要求的同时，可以减少大量的计算时间(相对于常用的实体单元和壳单元)，从而为大规模管线应力的分析提供了实际的可操作性。

另一方面，根据载荷数据库和材料数据库的相关参数，利用编制的损伤耦合有限元程序，计算出各种载荷下的应力应变和损伤；找出重点监测部位、对重点监测部位布置高温应变测量传感装置，监测其应变；建立应变-损伤-剩余寿命的关系数据库；根据在线监测的应变值，查询该关系数据库，得出管线剩余寿命的估值，实现对管道系统在线寿命的监测及评价。

以下结合优选实施对本发明一种高温管道损伤及在线寿命的预测方法作进一步的详细说明，以下的实施例仅用于进一步解释本发明，而非用于限定本发明的范围。

实施例1  管道系统T型接管损伤及寿命的在线预测

T型接头是工厂中广泛使用的一种复杂结构，由于在T型接头的连接部分存在着结构不连续，接管附近又有较严重的应力集中现象，因此，T型接头区域往往会先行失效并最终导致整个管道系统的损坏。使用本发明的预测方法对T型接头处的蠕变损伤发展情况进行模拟，通过计算获得相应的应变-损伤及应变-寿命数据库。

材料性能

T型接头的材料为0.5Cr0.5Mo0.25V钢，它在540℃下的蠕变参数如下：

表1  0.5Cr0.5Mo0.25V钢在540℃下的材料参数

 

E/MPa	μ	B	n	g	φ	A	α	v	ρ
										1.6×105	0.28	8.495×10-37	14.45	0.725	1.472	7.963×10-32	0.43	12.65	0.00534

几何尺寸

T型接管的几何尺寸和网格模型分别见附图3和附图4。

使用本发明的预测方法，计算内压为25MPa时该T型接管损伤和应变的发展情况。由附图5可以看出T型接管外壁面的损伤值始终高于内壁面的损伤值，而且随着时间的增加，外壁面处的A点先达到临界损伤值，随后内壁面的损伤量也跟着不断逼近临界损伤值。具体的计算结果见附图6～8；将计算的结果制作成本发明中提到的应变-损伤-寿命数据库，用于在线监测。

由计算可以得到该T型接管的实际寿命。如果监测到的蠕变应变值为0.031，如附图6所示，可以得出该T型接管的寿命分数为0.406。而相应的损伤可由附图8中查出，即D/D_cr＝0.6。

实施例2  管道系统弯管损伤及在线寿命的预测

弯管是主蒸汽管道中最常见的复杂结构，它是一个复杂的空间曲面构件。通过对主蒸汽管道的应力分析发现，弯管是整个管道系统中应力最大的部位。使用本发明的预测方法对某电厂主蒸汽管道的弯管部位进行了模拟，通过计算获得相应的应变-损伤及应变-寿命数据库。

材料性能

弯管的材料同实施例1。

几何尺寸

弯管的平均曲率半径R_m为1400mm，管道尺寸为270×40mm，实物图及有限元网格见附图9和附图10。

使用本发明的预测方法，计算内压为25MPa时弯管损伤和应变的发展情况。由附图11可以看出弯管内弧区的损伤值最大，监测时主要监测此处即可。具体的计算结果见附图12～14。将计算的结果制作成本发明中提到的应变-损伤-寿命数据库，用于在线监测。

由计算得出该弯管的实际寿命。如果监测到的蠕变应变值为0.049，如附图12所示，可以得出该弯管的寿命分数为0.704。而相应的损伤可由附图14中查出，即D/D_cr＝0.7。

Claims (5)

1、一种高温管道损伤及寿命的在线预测方法，其特征在于，利用测量得到的应变，结合有限元分析，实时预测高温管道的损伤及寿命，其实施步骤包括：

(1)对高温管道进行损伤耦合的有限元模拟分析；

(2)根据有限元分析结果，找出重点监测部位，对重点监测部位设置传感器，监测其应变；

(3)针对不同载荷工况，进行有限元分析，建立损伤-应变及剩余寿命-应变的关系数据库；

(4)将在线监测的应变值和载荷工况的数值，与损伤-应变及剩余寿命-应变数据库的数据进行在线查询对比，得出相应损伤和剩余寿命的评估值。

2、根据权利要求1所述的一种高温管道损伤及在线寿命的预测方法，其特征在于，在步骤(1)中，采用UG建立管道的几何模型，并导入有限元分析软件中，选用管单元分析管线的轴向位移、应变和应力，同时得到管道的周向应力；管单元在保证足够的精度要求的同时，可减少大量的计算时间(相对于常用的实体单元和壳单元)，从而为大规模管线应力的分析提供了实际的操作性。

3、根据权利要求1所述的一种高温管道损伤及在线寿命的预测方法，其特征在于，步骤(3)中所述的有限元分析，包含以基于连续损伤力学的本构方程为准则编制的含损伤计算的有限元分析子程序，通过对高温管道进行大规模的有限元模拟计算，反演出应变—损伤关系和应变—剩余寿命的关系。

4、根据权利要求1所述的一种高温管道损伤及在线寿命的预测方法，其特征在于，在步骤(4)中分析应变、温度、压力等现场数据，在损伤-应变和寿命-应变数据库中检索相应的应变所对应的损伤和剩余寿命，或者根据温度压力的变化分析损伤累积；利用WEB服务器响应终端用户的查询请求，并将应用服务器的处理结果生成相应的表单和参数曲线，通过Internet网页的形式发布；经过授权的用户可通过Internet浏览器实时查询管线的损伤状况和剩余寿命。

5、根据权利要求3所述的一种高温管道损伤及寿命的在线预测方法，其特征在于，所述的本构方程为：以连续损伤力学理论的本构方程为准则，通过对其模拟求解出相应的应变和损伤结果，包含：采用单状态变量的本构方程：

$\frac{{\mathrm{d\ϵ}}_{\mathrm{ij}}^c}{\mathrm{dt}}=\frac{3}{2}B{\left[\frac{{\mathrm{\σ}}_e}{1-D}\right]}^n\frac{S_{\mathrm{ij}}}{{\mathrm{\σ}}_e}{t}^m---(1-1)$

式中，

为蠕变应变张量，s_ij为应力偏张量，σ₁为最大主应力，σ_e为VonMises应力，D为损伤变量，其值在0～1之间变化；

α(0≤α≤1)是材料常数，用来表征多轴破断准则，B、n、A、υ是与最小蠕变应变速率以及断裂相关的材料常数；

或由多相复合体假设得出的本构模型：

$\frac{{\mathrm{d\ϵ}}_{\mathrm{ij}}^c}{\mathrm{dt}}=\frac{3}{2}{\mathrm{B\σ}}_e^{n-1}{s}_{\mathrm{ij}}[(1-\mathrm{\ρ})+\mathrm{\ρ}{(1-D)}^{-n}]---(1-3)$

$\frac{\mathrm{dD}}{\mathrm{dt}}=g\frac{A}{\mathrm{\φ}+1}\frac{{[{\mathrm{\α\σ}}_1+(1-\mathrm{\α}){\mathrm{\σ}}_e]}^{\mathrm{\υ}}}{{(1-D)}^{\mathrm{\φ}}}---(1-4)$

$D_{\mathrm{cr}}=1-{(1-g)}^{\frac{1}{\mathrm{\φ}+1}}---(1-5)$

式中，D_cr为临界损伤值，当 $\frac{D}{D_{\mathrm{cr}}}=1$ 时，表明材料达到了其蠕变寿命；

g、φ和ρ是指定非均质损伤的材料常数；

考虑到复杂晶相对蠕变变形和损伤演化的影响，可采用多状态变量本构方程，如两个状态变量的模型形式如下：

$\frac{{\mathrm{d\ϵ}}_{\mathrm{ij}}^c}{\mathrm{dt}}=\frac{3}{2}\frac{{\mathrm{G\σ}}_e^{n-1}{s}_{\mathrm{ij}}}{{(1-\mathrm{\Φ})}^n{(1-D)}^n}{t}^m---(1-6)$

$\frac{\mathrm{d\Φ}}{\mathrm{dt}}=\frac{k_c}{3}{(1-\mathrm{\Φ})}^4---(1-8)$

式中，第二个状态变量Φ描述了碳化物的沉积间隔，G和kc是材料常数；

将服从虎克定理的弹性应变和蠕变应变联立得：

$\frac{{\mathrm{d\ϵ}}_{\mathrm{ij}}}{\mathrm{dt}}=\frac{1+\mathrm{\μ}}{E}[\left(\frac{{\mathrm{d\σ}}_{\mathrm{ij}}}{\mathrm{dt}}\right)-\frac{\mathrm{\μ}}{1+\mathrm{\μ}}\left(\frac{{\mathrm{d\σ}}_{\mathrm{kk}}}{\mathrm{dt}}\right){\mathrm{\δ}}_{\mathrm{ij}}]+\frac{3}{2}{\mathrm{B\σ}}_e^{n-1}{s}_{\mathrm{ij}}[(1-\mathrm{\ρ})+\mathrm{\ρ}{(1-D)}^{-n}]---(1-9)$

首先，将有限元计算分为弹性域中的有限元和含蠕变应变的有限元分析两部分：

(1)弹性域中的有限元分析

有限元法首先将结构离散化成不同的单元，单元与单元间以节点相连；在弹性范围内，应力应变的本构关系可表示为：

{σ}＝[D]{ε}                 (1-10)

式中[D]为弹性矩阵；

若按直接法选定位移模式，则结构位移场与单元节点位移间的关系为：

{u}＝[N]{δ}                 (1-11)

式中{u}为位移场向量，[N]是坐标的函数，反映单元的位移形态，被称作“形函数”；δ为单元节点的位移向量；

又因为应变与位移存在以下关系：

{ε}＝[B]{δ}                  (1-12)

若作用在物体上的外力为P，则由虚功原理可得：

[K]{δ}＝{P}                  (1-13)

其中：

[K]＝∑[k]                   (1-14)

[k]＝∫[B]^-[D][B]dV          (1-15)

式中[k]为单元刚度矩阵，[K]为总刚度矩阵，[B]为应变矩阵，{P}为节点载荷向量；通过求解代数方程组(1-13)即可求得弹性应力解与应变解；

(2)含蠕变应变的有限元分析

若考虑带有蠕变变形的情况，则：

{ε^e}＝{ε}-{ε^c}                (1-16)

式中上标e代表相应的弹性分量，c代表相应的蠕变分量；

将(1-16)分别代入(1-10)和(1-11)中，得：

{σ^e}＝[D]({ε}-{ε^c})           (1-17)

{δ^e}＝[B]^-1({ε}-{ε^c})          (1-18)

这样(1-13)就可以表达成含有蠕变项的平衡方程：

[K]([B]^-1({ε}-{ε^c}))＝{P}    (1-19)；

其次，将非线性方程线性化：

在有限单元法中描述物体离散后的非线性特性，可以用下面的平衡方程式表示：

F^N(u^M)＝0                   (1-20)

式中，F^N是与第N个变量对应的分力，u^M表示第M个变量的值；

对于(1-20)式的非线性方程，假设第i步迭代后可获得的近似解为

令离散方程的精确解与近似解之间的差值为

即：

$F^N(u_i^M+c_{i+1}^M)=0---(1-21)$

方程左端在近似解

处按泰勒级数展开得：

$F^N\left(u_i^M\right)+\frac{{\∂F}^N}{{\∂u}^P}\left(u_i^M\right)c_{i+1}^P+\frac{{\∂}^2{F}^N}{{\∂u}^P{\∂u}^Q}\left(u_i^M\right)c_{i+1}^P{c}_{i+1}^Q+......=0---(1-22)$

忽略二阶以上的偏导数，得到如下的线性方程：

$K_i^{\mathrm{NP}}{c}_{i+1}^P=-F_i^N---(1-23)$

式中的 $K_i^{\mathrm{NP}}=\frac{{\∂F}^N}{{\∂u}^P}\left(u_i^M\right)$ 即为雅可比矩阵，并且 ${F_i}^N=F^N\left(u_i^M\right);$ 这样第i+1个近似解可表示为：

$u_{i+1}^M=u_i^M+c_{i+1}^M---(1-24)$

上式即为线性化后的迭代公式；

材料的雅可比矩阵定义为：

$J=\frac{\∂\mathrm{\Δ}{\mathrm{\σ}}_{\mathrm{ij}}}{\∂\mathrm{\Δ}{\mathrm{\ϵ}}_{\mathrm{ij}}}---(1-25)$

由前所述，若忽略塑性后总的应变式为：

${\mathrm{\ϵ}}_{\mathrm{ij}}={\mathrm{\ϵ}}_{\mathrm{ij}}^e+{\mathrm{\ϵ}}_{\mathrm{ij}}^c---(1-26)$

程序中每一载荷的施加都是通过两个载荷步来完成的：在第一步，通过定义一个服从虎克定理的材料模型来产生一个弹性响应；第二步蠕变分析中，实现如前所述的蠕变与损伤耦合的材料行为；

由于第一步是与时间无关的弹性响应， ${\mathrm{\ϵ}}_{\mathrm{ij}}={\mathrm{\ϵ}}_{\mathrm{ij}}^e,$ 因此，其应力增量与雅克比矩阵都较为简单，由弹性力学可知，线弹性条件下应力与应变的关系为：

{σ}＝[D]{ε}                 (1-27)

式中：

$\left[D\right]=\frac{E(1-\mathrm{\μ})}{(1+\mathrm{\μ})(1-2\mathrm{\μ})}\left[\begin{array}{cccccc}1& \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& 0& 0& 0\\ \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& 1& \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& 0& 0& 0\\ \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& 1& 0& 0& 0\\ 0& 0& 0& \frac{1-2\mathrm{\μ}}{2(1-\mathrm{\μ})}& 0& 0\\ 0& 0& 0& 0& \frac{1-2\mathrm{\μ}}{2(1-\mathrm{\μ})}& 0\\ 0& 0& 0& 0& 0& \frac{1-2\mathrm{\μ}}{2(1-\mathrm{\μ})}\end{array}\right]$

式中的[D]即为第一步中需计算的J；

在第二步中应力和应变都是与时间相关的，不同的时间所对应的J也不相同；由于蠕变损伤方程的表达式是非线性的微分方程，所以无法直接给出相应的雅可比矩阵，需要对蠕变损伤方程进行适当的处理；在子程序中引入一个显式的时间积分法，采用中心差分对目标方程进行离散化处理；采用的中心差分表达式和算符表示如下：

${\stackrel{\·}{f}}_{t+\frac{1}{2}\mathrm{\Δt}}=\frac{\mathrm{\Δf}}{\mathrm{\Δt}}---(1-28)$

$f_{t+\frac{1}{2}\mathrm{\Δt}}=f_t+\frac{\mathrm{\Δf}}{2}---(1-29)$

式中f为一任意函数，f_t表示该函数在增量开始时的值，Δf表示在该增量中函数的变化值，Δt是时间增量；

把蠕变损伤连续方程(式(1-9)和式(1-4))按中心差分离散化后得：

$\frac{{\mathrm{\Δ\ϵ}}_{\mathrm{ij}}}{\mathrm{\Δt}}=\frac{1+\mathrm{\μ}}{E}[\left(\frac{{\mathrm{\Δ\σ}}_{\mathrm{ij}}}{\mathrm{\Δt}}\right)-\frac{\mathrm{\μ}}{1+\mathrm{\μ}}\left(\frac{{\mathrm{\Δ\σ}}_{\mathrm{kk}}}{\mathrm{\Δt}}\right){\mathrm{\δ}}_{\mathrm{ij}}]+\frac{3}{2}B{({\mathrm{\σ}}_e+\frac{{\mathrm{\Δ\σ}}_e}{2})}^{n-1}(s_{\mathrm{ij}}+\frac{{\mathrm{\Δs}}_{\mathrm{ij}}}{2})[(1-\mathrm{\ρ})+\mathrm{\ρ}{(1-D)}^{-n}]$

                                           (1-30)

$\frac{\mathrm{\ΔD}}{\mathrm{\Δt}}=g\frac{A}{\mathrm{\φ}+1}\frac{{[{\mathrm{\α\σ}}_1+(1-\mathrm{\α}){\mathrm{\σ}}_e]}^{\mathrm{\υ}}}{{[1-(D+\frac{\mathrm{\ΔD}}{2})]}^{\mathrm{\φ}}}---(1-31)$

由于蠕变损伤方程在中心差分后，所得的雅可比矩阵J是一个6×6阶的非对称矩阵；直接给出该方程式的应力增量对应变增量的偏导数存在难度，而应变增量却可以在非线性方程组里以闭合形式给出，所以，首先计算应变增量对应力增量的偏导数：

$J^{-1}=\frac{{\∂\mathrm{\Δ\ϵ}}_{\mathrm{ij}}}{{\∂\mathrm{\Δ\σ}}_{\mathrm{ij}}}---(1-32)$

然后通过求逆J^-1得到J；

将以上方法编制成UMAT用户子程序的步骤为：先通过判断模块判断应力状态，再分别进入弹性响应模块和蠕变损伤响应模块进行计算；然后分别转入弹性响应模块和蠕变损伤响应模块的子程序模块进行相应计算：弹性响应模块给出材料在弹性范围内的雅可比矩阵J和相应的应力应变值；蠕变损伤响应模块最后给出材料在增量步下的雅可比矩阵J；

蠕变损伤响应模块中的计算损伤增量模块和其子程序一起计算，得到在当前应力状况下损伤增量(ΔD)的大小，并对方程(1-31)进行变形，得到外部函数计算损伤函数的函数值：

$f\left(\mathrm{\ΔD}\right)=g\frac{A}{\mathrm{\φ}+1}\frac{{[{\mathrm{\α\σ}}_1+(1-\mathrm{\α}){\mathrm{\σ}}_e]}^{\mathrm{\υ}}}{{[1-(D+\frac{\mathrm{\ΔD}}{2})]}^{\mathrm{\φ}}}-\frac{\mathrm{\ΔD}}{\mathrm{\Δt}}---(1-33)$

由于D的值在[0，1]区间单调递增，所以ΔD的变化也不会超过该区间；令ΔD在0至1间以0.05的倍数分别代入(1-31)式进行试根，若f(ΔD₁)×f(ΔD₂)≤0，由零点定理可知公式(1-26)在[ΔD₁，ΔD₂]上有根存在；分别以ΔD₁和ΔD₂作为上下限，代入子程序计算损伤函数的根，利用布伦特法保证所获得的ΔD的误差不超过一定容限；其中最大主应力σ₁是通过调用有限元软件(ABAQUS)内部提供的SPRINC子程序得到；

计算应力增量模块包含了外部函数计算范数和三个子程序模块：计算雅可比矩阵模块、LU分解模块和解线性方程组模块；计算雅可比矩阵模块的功能是把损伤增量代入到公式(1-15)中，形成六个联立的非线性方程组，采用牛顿-拉斐森法搜索六个应力增量的大小；若六个应力增量符合如下不等式：

$\underset{j=1}{\overset{3}{\underset{i=1}{\mathrm{\Σ}}}}\{\frac{1+\mathrm{\μ}}{E}[\left(\frac{{\mathrm{\Δ\σ}}_{\mathrm{ij}}}{\mathrm{\Δt}}\right)-\frac{\mathrm{\μ}}{1+\mathrm{\μ}}\left(\frac{{\mathrm{\Δ\σ}}_{\mathrm{kk}}}{\mathrm{\Δt}}\right){\mathrm{\δ}}_{\mathrm{ij}}]$

$+\frac{3}{2}B{{({\mathrm{\σ}}_e+\frac{{\mathrm{\Δ\σ}}_e}{2})}^{n-1}(s_{\mathrm{ij}}+\frac{{\mathrm{\Δs}}_{\mathrm{ij}}}{2})[(1-\mathrm{\ρ})+\mathrm{\ρ}{(1-D)}^{-n}]-\frac{{\mathrm{\Δ\ϵ}}_{\mathrm{ij}}}{\mathrm{\Δt}}\}}^2\≤\mathrm{\ξ}---(1-34)$

式中，ξ为允许的误差，则认为此六个应力分量就是所求值。

Images