管道稳态裂纹扩展模拟

 管道稳态裂纹扩展模拟

2.1断裂动力学理论简述
断裂力学学科的先导者是A.A.Griffith，他在1920年首先提出将裂纹临界扩展的判据与裂纹的长度定量地联系在一起，建立起脆断理论的基本框架。断裂力学的蓬勃发展则以1948年Irwin和Orowan分别独立建立的工程材料脆性断裂理论为标志。
作为断裂力学的一个重要分支，断裂动力学诞生的标志是1948年英国物理学家N.F.Mott在Griffith理论中考虑动能的影响后所发表的论文。1951年印度女科学家Elizabeth Yoffe最早给出了运动Griffith裂纹的解析解。然而断裂动力学中最重要基本概念的提出，系统分析方法的形成，以及相对成熟的实验研究方法的建立，是20世纪80年代以来的研究成果。本文的主旨是用数值模拟和实验方法解决管道动态断裂评估问题。
2.1.1静止裂纹与动态扩展
断裂力学认为在结构中不可避免有类似裂纹的缺陷存在。在小变形、低能量耗散的情况下，可看作脆性断裂；对于大变形、高能量耗散的情况则按延性断裂来处理。受应力集中影响，除去理想脆性材料外，外加载荷在裂纹尖端附近均伴随非弹性区。若该区域尺寸与其他特征尺寸相比为小量，则可以用线弹性断裂力学处理。裂纹根据加载方式的不同可分为三种模型：张开模式（I型）；滑开模式（Ⅱ型）；撕开模式（Ⅲ型）。本文讨论的范围限于I型裂纹。在断裂力学中，定义引起裂纹产生单位长度的扩展所需要的能量为裂纹驱动力G。在平面问题中，G与应力强度因子K_c的关系可以表达为：

对于平面应力问题，E′=E；对于平面应变问题，E′=E/（1-v²）。其中，E是材料的杨氏模量，v为泊松比。

断裂动力学（FractUre Dynamics）也叫做动态断裂力学（Dynamic Fracture

Mechanics)，其目的是研究那些惯性效应不能忽略的断裂力学问题。这些问题主要归纳为两大类：一是裂纹稳定而外力随时间迅速变化，如振动、冲击、波动等：二是外力恒定或缓漫变化而裂纹发生快速传播。这两类断裂动力学的问题分别称为裂纹动态起始问题和运动裂纹问题。
运动裂纹问题从现象上看可以分为前期加速，扩展轨迹，扩展速度，分叉和止裂几方面问题。本文研究裂纹从快速扩展开始直至止裂的过程。

2.1.2裂纹扩展的极限速度

输气管道上的脆性裂纹扩展速度曾经达到过10³m/s而S的量级。通过改善钢材韧性，降低韧脆转换温度等办法，裂纹扩展的主要形式由脆性断裂转化为延性断裂，裂纹扩展速度也有大幅度的下降。伴随着超高压管道上的裂纹扩展，近年来的实测速度又有上升。
对于本文的裂纹动态模拟和止裂评估，裂纹扩展速度是关键变化量。那么裂纹扩展存在极限速度吗？如果存在，如何量化？本小节旨在探讨这方面的内容。尤其针对输气管道上的应用，文中给出了大致的估算。
1948年Mott认为快速裂纹扩展过程中的动能作用不可忽略。考虑无穷大弹性板中的裂纹扩展，加入动能项的能量平衡方程为：

式中a为半个裂纹长度，U为弹性应变能，K为动能，Γ为表面能。

在推导动能表达式的过程中，Mott引入以下假定:
●围绕裂纹尖端区域的应力场、位移场，可由静态弹性理论确定；

●裂纹扩展速度远小于柱形杆纵波速度C_o= ；
●裂纹扩展阻力不随裂纹速度变化而改变。
在上述假定的前提下，得到无限大弹性板中的动能表达式:

式中k是待定的比例系数，p为密度，v为裂纹扩展速率，σ为无穷远处承受的均匀拉伸应力，E为材料的杨氏模量。

对于无限大板中的裂纹而言，应变能U和表面能厂的表达式同静态清况相比没有形式上的改变。将三种能量的表达式代入能量平衡方程（2-2)，可得：

将单位面积的表面能厂用表示为临界裂纹长度a_c的形式，经过Berry等人的修正，得到裂纹速度的表达式:

Roberts(1954）用数值方法计算了泊松比v=0.25时的系数 ＝0.38,从而得到钢材中的极限裂速v_m≈1929m/s。
在双悬臂梁DCB(double-Cantilever-Beam）试样中测到的结果比上述值略低，约为1500m/s。
Freund(1972）通过分析推导，认为裂纹的极限裂速v_m，应该等于Rayleigh波速C_R。按此方法估算，v=0.25时钢材中的极限裂速v_m=294m/s。

根据Rice(2000）通过奇异裂纹模型做出的最新论断，在典型的远场加载条件下，在I型和Ⅱ型时的极限速度为Rayleigh波速，而在Ⅲ型时的极限速度可达横波波速，对钢材而言约3100m/s。
Rice总结了I型拉伸裂纹的实验观测结果，发现:
●在脆性非晶态材料（如玻璃、PMMA）中，裂纹速度的上限为0.55C_R～0.65C_R;
●裂纹速度v<0.3C_R～0.4C_R时，断裂面呈镜面光滑状。在更高（平均）速度时，裂纹表面变得非常粗糙，并且v开始剧烈振荡；
●存在v接近于C_R的例外情况，如高度各向异性的脆性单晶（钨、云母）,及不完全烧结的固体。
这一结果同钢材断裂从DBF到DDF的发展过程中，实测速度与断口形貌的变化趋势相吻合。简单的说，就是韧性提高可以导致断裂的极限速度降低。

压力管道的裂纹扩展与无限大弹性板不同，Kanninen等（1980）利用类似一维梁的模型提出了输气管道中的裂纹扩展极限速度的估算值。模型中采用弹性基础梁的挠度模拟在对称载荷作用下圆柱壳的变形，并作了以下基本假设：

●以径向变形为主；
●压力沿圆周的变化可以忽略不计；
●裂纹张开位移等于在裂纹区任何截面径向位移沿圆周的积分；

●塑性屈服铰在裂尖后部形成；
●裂纹扩展速度超过流体降压速度，扩展中裂尖后面的压力为零；

●由于裂纹出现导致壳体刚度突变。
得到管道裂纹的极限速度：

式中C₀为柱形杆纵波波速，钢材约5076m/s,h和R为管道壁厚和平均半径。按此式估算的西气东输管道上裂纹扩展的极限速度约为648m/s。

对于裂纹在脆性管道中的快速扩展，在全尺寸实验中可观察到的裂速范围为600～1000m/s，考虑到不同的设计参数，与上述预测比较接近；全尺寸实验在近年来应用的高韧性管线中观测到的裂纹稳态扩展速度一般在150～350m/s,近似可以看作脆性断裂速度的1/3，即：

在Kanninen之后，Emery（1950）提出，尽管裂纹的起裂和小范围扩展可由系统的降压值作为上限而加以估算，但开裂管道中的裂纹扩展却受到流体外流引起的压力降低的强烈影响。在几何尺寸一定的条件下，由管道断裂引起的压力释放波可以由远端反射回来，从而增加了裂尖处的管道压力：同时张开裂纹的自由边，却承受一个显著减少的压力，致使裂纹扩展减速。因此，为了解裂纹的动力学特征，必须考虑管道变形及流体压力。
在Emery的圆柱壳模型下，极限裂纹扩展速度比（2-6）式略低。
2.1.3动态断裂力学参数

断裂动力学理论为结构动态断裂的分析与控制奠定了理论基础，并为其在工程应用方面提供了重要的概念、分析方法、结构参数计算和动态断裂准则。基本概念可以归纳为三个方面：描述动态断裂的特征参数、材料的动态断裂韧性和运动裂纹的止裂判据。接下来对上述三方面问题，尤其是对于和管道问题相关的概念，加以阐述和介绍，以备后文直接引用。
本节的讨论范围限于I型裂纹快速扩展问题，主要目的是给出与断裂力学参数相应的动态应力强度因子、动态能量释放率和裂纹尖端张开角的表达式。

2.1.3.1动态应力强度因子
I型动态应力强度因子  (t）的定义为：

式中t为时间，a_yy为y向正应力分量。
自Yoffe（1951)开始，包括形Rice(1968)、Freund(1990)在内的很多学者推导并不断完善了无限大板等简单模型的以  (t）表示的稳态扩展的裂纹顶端渐进应力场和位移场的解析解，Nishioka等（1996）进而给出了包括瞬态情况下的四阶渐进展开式，在此不做赘述。
这样只要得到  (t）的值，就可以按照上述文献中的公式得到裂纹尖端场的解析解。目前除对于无限介质中的稳态扩展裂纹，通过适当简化可以得到应力强度因子的解析表达式以外，有限介质及构件中运动裂纹的动态应力强度因子的计算主要依靠数值方法，包括有限元法和有限差分法。
本文研究的高韧性钢为弹塑性材料，在裂纹扩展前，往往在裂端区甚至更大范围内有相当大的塑性变形，且伴随着裂端后面的卸载。因此，起裂后必须克服塑性变形才能发生裂纹扩展。此时，作为衡量裂端区应力场强度的力学参量J积分和应力强度因子K并不是严格有效的。Kanninen提出了更先进的J积分-T积分，适合应用于在弹塑性材料中的动态裂纹扩展分析。
2.1.3.2动态能量释放率
对含运动裂纹的一般弹性体，为确定释放到裂尖的能量，考虑围绕裂尖非常小的闭合回路Γ，如图2-1所示，动态能量释放率可以由能量流动定义为：

式中v_n表示裂纹尖端运动速度在回路Γ上的法向分量，v为裂尖速率，T_i为作用于Γ上的应力分量，W为应变能密度，p为物质密度。

这里Γ的形状是任意的，但必须附关于裂尖的运动坐标系。这个表达式对线性与非线性的材料均适用。

在裂纹稳态扩展的条件下， （t）退化为断裂力学中的动态J积分，是与路径无关的量：

对于无限介质而言， （t）和 （t）都可以写成与裂纹速度有关的函数与静态因子乘积的形式：

（t）=k（v）K（0）                   （t）=g（v）G（0）    （2-9）

将用动态应力强度因子表示的无限大平面弹性体应力场和位移场的表达式代入动态能量释放率的定义式，可得二者之间的对应关系：

式中F（v）可以表示为裂尖速率v与材料平面纵波波速C₁及平面横波波速C₂的函数，E是材料的杨氏模量，v为泊松比。

（2-10）式意味着 （t）和 （t）同时达到临界值。也就是说，动态扩展断裂判据既可用临界应力强度因子，也可用临界能量释放率来表示。这是对线弹性材料、恒定裂速扩展的裂纹在无限大弹性体中得出的结论，可以定性地推广到一般问题中使用。

2.1.3.3裂纹尖端张开角
裂纹尖端张开角CTOA（Crack Tip Opening Angle）最早被用来分析核电站管路的稳态裂纹扩展，后来用于输气管线和航空工业。如图2-2，与裂纹张开角COA不同，CTOA对应于一定的断裂扩展状态，不随测量点的位置而改变。CTOA在数学上可以按下式定义：

式中：δ_t为裂纹尖端张开位移，△a为裂纹扩展长度。

2.1.4材料的动态断裂韧性
作为材料抵抗动态断裂能力的度量，引入韧度的概念，称为动态断裂韧性。具体材料的韧性值由规定条件下的实验测定。
表征材料抗力的动态断裂参数有冲击断裂韧性K_Id，动态能量释放率临界值G_d，动态J积分临界值J_d，裂纹尖端张开角临界值（CTOA)c，甚至包括直接从实验机上获得的夏比冲击功和落锤吸收能量等等，均可以作为判断裂纹扩展与止裂的参数。工程上如何选用要视具体问题与实验条件而定。
动态韧性参数除与静态断裂韧性同样与材料性质，环境温度，应力状态等因素有关以外，同时还与裂纹速度有关。这是本文关心的问题。
裂纹发生快速运动时，一个高应力集中区被驱动穿过材料，因此裂纹顶端附近的材料产生很高的应变速率或加载速率。在一般情形下，高加载速率使得材料变脆，因而抗力降低：另一方面，随着裂纹速度的提高，由于裂纹顶端应力场多轴化程度降低，又增强了塑性效应，从而提高了材料的抗力。另外，快速扩展引起的热耗散使得裂纹顶端附近温度上升，也会影响材料的抗断性能。Eftis,Klaft(1965）和Barton,Mall(1960）等的实验结果表明，动态断裂韧性由裂纹起始扩展的材料抗力开始下降，当裂纹速度达到一定值时，动态断裂韧性又随裂纹速度的增长而提高，表明了前面分析的正确性。
在裂纹稳态扩展阶段，CTOA达到临界值而保持恒定不变，称作（CTOA）c。(CTOA）c可在小试样冲击实验中测得，并直接应用到全尺寸管道的断裂上。在描述断裂性能方面，(CTOA）c是衡量延性断裂扩展抗力的指标之一，它被认为可以替代夏比冲击韧性，以解决在高韧性管线上依赖后者而出现的重大偏差。

Kanninen等（1992）认为，对于动态裂纹扩展问题，G_d和（CTOA)c是最适合的断裂韧性描述参数，如表2-1。
表2-1断裂力学基本处理方法

2.1.5运动裂纹的止裂判据

含裂纹结构在动载作用下将发生裂纹的起裂和扩展，导致结构破坏或者发生止裂，产生这些过程的条件称为判据。除了1.2,1节提到过的减压波速判据及在它基础上发展起来的双曲线方法以外，还有应力强度因子K判据，能量释放率G判据和裂纹尖端张开角CTOA判据等。表2-1同样给出了参考性意见。

在止裂条件的建立上，存在两种不同的观点，即动态观点和静态观点。下面用K判据为例介绍一下两种观点的区别和发展，以及本文的选择。动态观点认为，裂纹的快速扩展直至止裂是一个完整的过程，止裂是这种过程的结束，并忽略惯性力的影响。这样，裂纹持续快速扩展的条件为：

=K_ID（v，T）                             （2-12）

即把裂纹动态扩展断裂韧度K_ID表示成为裂纹速率v和温度T的函数。注意这一条件并不适用于起裂时刻，因为K_ID(O，T刀不等于静态韧性K_IC（T）。

按照上面的观点，从止裂发生的t₀时刻一直到裂纹速度降为零的t_a时刻， 加均应小于扩展断裂韧度K_ID（v，T）的最小值，止裂判据表示为：

这样，通过数值方法解出特定问题的K_I(t)，和材料实验定出的K_ID（v，T）相比较，可以得到是否满足止裂条件的结论。
静态观点着眼于止裂瞬间的裂纹尖端条件，而不顾及已经发生的裂纹动态扩展过程和历史。这样，基于静态观点的止裂判据表示为：
≤K_Ia                                                       （2-14）

式中 是对应于止裂长度的静态应力强度因子，K_Ia是材料的静态断裂韧性。两种学术观点的争论长达二十多年，并都得到了特定条件下实验结果的证实。直至1985年方达成一致，解释如下：
引起动态效应的主要原因是由于试样边界对应力波的反射，当边界反射的应力波几乎不返回到裂纹尖端的情况下，采用静态观点相对合适；反之应采用动态止裂分析观点。这一结论得到了实验的证明。
压力管道尤其是高压输气管道上的裂纹扩展是强几何非线性的过程，同时不断伴随应力波的传播与反射，因而本文在分析裂纹扩展与止裂时用了动裂纹约束的概念（动态观点）。根据这一概念，裂纹止裂就是使裂纹不能连续扩展。现将本文用到的以G和CTOA建立的动态判据分述如下：
对于天然气管道，能量释放率判据可以表示为:
G(a,p,D,SDR,E)=G_d(T,v,h)                                    (2-15)

式（2-15）中，当G_max=G_d时，裂纹扩展；若G<G_dmin，裂纹止裂。G表示裂纹驱动力，G_d表示材料的动态断裂韧性，a为裂纹长度，p为流体压力，D为管道直径，SDR(Standard Dimension Ratio）为管道外直径与壁厚的比值，E为材料的杨氏模量，T为温度，v为裂纹扩展速率，h为试样壁厚。
为了防止裂纹在工程管道上扩展事故的发生，G和G_d的相对值为实际应用准则提供了定量的形式。G来源于对结构模型的分析和计算，而G_d要通过对结构材料的实验来分析确定。

类似的，在管道延性裂纹扩展过程中，表征裂纹驱动力的（CTOA）_max与表征管壁断裂韧性（CTOA）c构成止裂判据：

（CTOA）_max（a,p,D,SDR,E）=（CTOA）c（T,v,h）               （2-16）

当（CTOA）_max ≥(CTOA）c时，裂纹扩展；若（CTOA）_max＜(CTOA）c，则裂纹止裂。

2.2断裂动力学的有限元法

管道裂纹动态扩展问题的关键是正确描述管道的实际工作状态。断裂动力学问题的复杂性和目前的研究深度决定了分析解法的求解必须对加载形式、材料参数、裂纹几何和运动状态作很多限制，因而本文以有限元方法为主，在表征韧性的参数确定和计算结果参照中引入实验手段，在试件比拟的塑性区范围确定等处用到了解析方法。
断裂动力学中的有限元法分为动力学分析和断裂分析两方面。在动力学分析中要考虑材料的惯性效应，在动态变化的载荷下求解运动方程；在断裂分析方面，要考虑裂纹尖端场的性质，以及裂纹扩展引起的边界变化等。下面分别从这两方面介绍本文计算中用到的有限元技术及相关处理方法。

2.2.1壳体结构动力学基本方程

壳体有限元计算的一般思路是在单元内描述问题，再根据节点集合。在描述几何方程和运动方程之前，首先说明符号的格式。本节将用到矢量符号、矩阵符号和上下标。矢量符号通过黑体表示，如a，它的分量表示为：a_i，也可以以矩阵｛a｝的形式表示。a_iI是节点I的a矢量的第i个分量。
在本节中，将会用到图2-3所示的三种坐标系：

●空间固定的总体坐标系（x，y，z），或x_i，正交基矢量为 。
●节点坐标系行（  ）习，反映节点的平动与转动；正交基矢量为b_i。

●单元随体坐标系（  ）；正交基矢量为e_i。
对于总体坐标系中的任意矢量a，其在节点坐标系中的分量可以表达为：

其中，b_xl，b_yi，b_zl是节点基矢量b_i在总体坐标系中的三个分量。

类似地，对于单元随体坐标系：

其中[b]^T、[e]^T分别是总体坐标与节点坐标、总体坐标与单元随体坐标之间的转换矩阵。

节点J的初始位置为： （i=1，2，3代表x，y，z三个方向的分量），当发生了位移u_iJ后，新的位置坐标为x_iJ。

x_iJ= +u_iJ                                               （2-19）

2.2.1.1几何方程

基于Kirchhoff假设的板壳单元，壳体中某一点的应变可以由中面应变  和曲率  来表示：

其中[N^P]为壳体面内位移形函数，[N^b]为壳体弯曲位移形函数。[δ^P]、[δ^b]是相应的节点位移。

将方程（2-22）代入（2-21），得到

2.2.1.2运动方程

对于惯性力不可忽略的结构进行分析，要求解运动方程，以获得系统在瞬时状态下的位移和速度、加速度。一般的方法是按照平动和转动的自由度，分别给出运动方程。在有限元中，对于离散化的节点，要相应地给出它的三个平动方程和三个转动方程。

三个平动方程为：

其中i代表节点号，m_i为节点质量，  为平动加速度， 为外力， 为变形引起的内力（其中j=1，2，3，分别表示三个坐标轴方向x，y，z）。

三个转动方程为：

其中I_ij为转动惯量， 分别为角速度和角加速度， 为外力矩， 为变形引起的内力矩。

公式（2-24），（2-25）中的内力分量 和 可以如下求得：

首先根据虚功原理，在一个单元E内，对节点i有：

其中对I 不求和；；j代表各方向，对j求和。 为单元E的应变率，V是单元体积。由上述方程可以得到节点内力和节点内弯矩：

利用叠加原理，将每一个单元内的节点上的内力和内弯矩分别求和，可得到节点上的相应分量：

其中， 是如下定义的相关系灵敏：

=1     当单元E中的节点A与节点i相关；

=0     当单元E中的节点A与节点i无关。

将方程（2-28）代入方程（2-24）和（2-25），可得到各节点的加速度与角加速度：

对于每一时间步长，由运动方程首先解出各节点的平动加速度，然后根据显式中心差分公式解出各节点速度：

利用位移可以根据几何方程和本构方程求得应变和应力。

2.2.1.3本构方程

本文研究钢制管道的裂纹扩展。忽略轧制过程的影响，钢材可看作各向同性材料，则本构关系可以写成分段线性弹性的形式。进一步考虑到本文采用节点力释放法进行裂纹扩展计算，裂尖附近的应力集中影响不大，故可采用线弹性本构关系计算。

对于板壳单元，在单元随体坐标下应变和应力均只有三个分量：

2.2.1.4局部坐标下的边界条件

一般地，可以在边界处约束任意平动与转动分量，即令其为零。如要定义的运动分量不在总体坐标系下，则需在该节点定义局部坐标（X，Y，Z）。计处牙可以通过定义X和Y坐标的方向余弦c_xX，c_yX，c_zX，c_xY，c_yY，c_zY来实现。Z坐标的方向余弦可通过其正交性由下式求得：

将上面的转换矩阵定义为[c]，则边界条件中的平动速度可按（2-36）式转化为局部坐标中的分量，转动角速度按下式定义：

然后将被约束的分量置零，求解。求解后的平动速度和转动速度可分别通过（2-36）和（2-37）的逆变换转化成总体坐标中的分量。因为新的角位移量从方程(2-31）和关于b₁的方程直接解得，边界条件通常反映在更新的位移和节点向量b₁中。
2.2.2断裂动力学基本方程

2.2.2.1生成模式与扩展模式

在管道裂纹扩展问题中，因动态断裂发生的时间很短，以及裂纹扩展速度的非预先确定性，使得测量高阶的物理量如能量分布、瞬时动态能量释放率和动态裂纹尖端应力场变得十分困难。为此针对管道快速断裂问题发展了两种分析计算模式：生成模式（Generation mode）和扩展模式（Propagation mode）。

生成模式就是给定管道的几何尺寸和工作条件，同时必须给出已知的以裂纹扩展长度或者裂纹速度表示的裂纹扩展历史作为补充条件使得问题可解，从而计算裂纹驱动力G（或CTOA和评估断裂韧度G_d（或（CTOA）c）。
本章计算中的裂纹扩展速度取为定值，该值未必是裂纹稳态扩展的真实速度。然而通过指定不同的速度值，可以寻找对应于某一特定稳态裂纹扩展速度的G或CTOA，如果该值小于材料韧性，可以判断发生止裂；否则可能发生扩展。从而得到一些有价值的结论，例如在什么速度下管道失稳断裂等等。这种分析方法是运用了生成模式。
扩展模式是通过实验测量裂纹速度和压力分布，分析G_d（或（CTOA）），找到它与裂纹速度之间的关系，分析确定产生裂纹扩展需要的驱动力，从而得到一定条件下是否止裂的结论。包括第一章提到的双曲线法和后面提到的韧性减速机理等以实验和结合实验的计算为主的分析方法，就是采用了扩展模式。

2.2.2.2节点力释放技术
节点力释放（Node Foroe Release)，顾名思义就是当裂纹顶端通过有限元网格中的某节点时，便解除该节点在网格中本来起到的连接作用，即将其分为两个节点，并释放连接力，如图2-4。

在典型的生成模式计算中，可根据指定的a(t）或v(t）判断裂纹经过的时间和路线，依次释放节点力；另一种运用扩展模式的方法是通过判断裂纹顶端被约束的节点力达到F_c时，这个节点便被松弛。这里的F_c是指定的，它与网格尺寸以及动态断裂韧性K_ID有关。Keergstra(1976）证明，对于给定的网格，裂纹顶端的节点力正比于应力强度因子K_ID。

早期将有限元法应用于动态裂纹扩展时，裂尖运动由不连续的突进进行模拟：在时间增量△t内，裂尖沿裂纹方向从单元的一个节点改变到下一个节点。为了得到比较精确的解答，必须采用较小的时间步长△t，通常为膨胀波在两个最近单元节点间传播所需要的时间。由于裂纹传播速度通常明显低于波速，在时间△t内裂尖实际上通常只能运动到两个相邻节点之间的某一位置。在裂尖从一个节点移动到下一个节点的过程中，裂纹长度的突然增加和位移约束的突然解除，将引起有限元求解中比较严重的高频振荡现象。为了克服这一困难，产生了“回复力”法，Lagrange乘子法，节点力释放率法等缓步释放节点的方法,下面介绍本文采用的节点力释放率法。

在生成模式中，输入信息包括作为时间函数的裂纹尖端位置，每一时间步长中裂纹前进的距离是已知的。裂纹沿着单元扩展的过程中单元节点力逐步释放，此时，能量释放率G，即在裂端区每单位面积上裂纹扩展所引起的能量变化，可以近似地以节点力作功的形式表示:

式中，h是管理壁厚度， 是沿裂纹扩展方向的一个单元长度，△t是裂纹扩展一个单元所用的时间，v_n是垂直裂纹扩展方向的节点位移速度，系数2表示管道对称计算的两个部分。

在公式（2-38）中，F是被释放的约束力，它的大小随着裂纹在单元边界上的位置呈线性变化，表示为：

式中的F_o是被释放的节点约束力，a（t）是裂纹在单元上的扩展长度。指数c曾取过1/2、1、3/2、2等不同取值。本文在与实验对比的基础上取1。

2.2.2.3裂纹驱动力G的能量表达

能量平衡方法计算裂纹驱动力基于Griffiths断裂理论。在本章的计算程序中，外力作功和内部能量的平衡常被用来校核计算结果的数值稳定性，也可将算得的能量值带入Irwin-Orowan能量平衡公式计算裂纹驱动力。该方法与节点释放法更为简便。

Griffiths奠定了断裂问题的能量平衡理论。在裂纹扩展过程中，物体内部能量释放所产生的裂纹驱动力导致了裂纹增长，同时存在阻止形成新裂纹面积的阻力，当裂纹增长da长度时，二者形成平衡。这样裂纹驱动力G可以表示为：

式中h表示裂尖处材料厚度，W表示外力作用于裂纹体的功，U表示应变能。在动裂纹扩展问题中需要考虑动能K的作用，上式扩展为：

早期对流体压力管道裂纹扩展问题的分析中曾应用应变能变分原理，将裂纹动力等效为裂纹前面环向应变能的释放率，Iewin-Conten发展了裂纹驱动力的解析计算：

式中p₀为初始内压，D₀和h为管壁外径和壁厚，E为材料的杨氏模量。

2.2.3数值稳定性分析
在显式时间积分中，除非时间步长选得足够小，否则会出现数值不稳定性，从而使计算结果失去意义。尽管对于完全非线性问题，不能进行稳定性分析，但是经验表明线性分析也会给出有用的时间步长估计。一般说来，将线性稳定的时间步长减少10%到20%，可以克服非线性带来的不稳定性。
对于中心差分，线性分析表明时间步长

其中， 是系统的最高频率， 是最高频率下的临界阻尼系数。
ω是下面问题的特征值：
[K]{X}=ω²[M]{X}                                        （2-44）

其中[K]是系统的刚度矩阵，[M]是质量矩阵。
在实际计算中，我们并不是直接计算系统的ω_max，而是计算单元的 。在所有的单元的 中选择最大值作为ω_max。
对于板单元，有

其中，C_l是平面纵波波速，对四边形板单元来说，L是短边的长度。

2.2.4气体压力模式
1.3.3节提到了解祸算法中气体压力模式的广泛应用。对于含裂纹扩展的管道，其内部管壁受到的气体压力由两种模式组成：裂纹前面气体减压传播和裂纹后部的气体压力衰减。
Battelle的Maxey等人在实验中采用压缩空气、氮气以及甲烷含量高的天然气观察气体的减压行为，发现理想气体规律足以满足这些气体的状态等式。假设气体的膨胀过程是等熵过程，充满气体的管道突然发生开裂，整个状态中气体成分混和均匀，则裂尖附近气体局部压力与裂纹扩展速度有如下关系:

式中p_d为减压后的压力水平（MPa)；p_i为开裂前管内压力（MPa)；v为裂纹扩展速度（m/s)；v_a为初始压力和温度状态下的声速（m/s)； 为起始状态下气体的绝热系数，C_P和C_v分别为定压比热和定容比热。
从（2-46）式可以看出，当裂纹在管道上扩展时，裂纹前面的气体减压传播，其数值低于气体初始压力值p₀，气压减少值主要取决于裂纹的扩展速度，对于稳态扩展的长裂纹，管道内部气体减压趋于稳定的数值。

当裂纹沿管壁快速扩展时，裂纹后面的气体逸出，气压降低，但是与液体压力不同的是，该区域内的气体压力并不是立即降为零值，而是稳定地衰减至大气压力，并有一段明显的衰减区。在这一衰减区内，气体压力仍然作用在已经开裂的管壁上，加速管壁的变形。
根据实验给出的实际管道开裂实验的结果，发展了裂纹后面气体压力衰减模式，将气压的衰减表示成裂纹位置和初始压力的函数。常见的衰减模式有指数衰减的Fourier级数和抛物线衰减函数等。
更为简便的是图2-5所示的气体压力衰减模式是线性衰减函数，图中以裂纹扩展方向z为横轴，压力分布p为纵轴，其计算表达式为

P（z）=P₁（1- ），z<L                                                 （2-47）

式中z是计算截面的位置，L是衰减长度。在计算中，L取1.5倍或2倍的管道外径。
2.3管道稳态裂纹扩展问题的数值模拟
本节着重介绍稳态扩展条件下的开裂管道计算结果，力求从多个角度模拟管道变形的真实状态，为建立止裂判据与后文的减速机理打下基础。

2.3.1计算模型与边界条件
管道上的裂纹一般从起裂点同时向两侧扩展。因管道几何形状和载荷分布具有对称性，本文取一侧管道的一半来分析和计算，即四分之一模型，如图2-6所示。采用四节点板壳单元，沿管道轴向和环向分别划分网格，在管道内壁给出气体压力分布。

在对西气东输管道进行计算时，采用的管道模型长度L=37.5m，直径D=1.016m，厚度h=0.0147m，沿轴向划分为250个单元，环向划分为16个单元，共4000个单元，4267个节点，每个四边形单元轴向长度0.15m，环向长度0.lm。下文如无特殊说明，则网格均按此划分。
在图2-7中，C点表示裂尖位置：裂纹沿LINE4扩展，在A、B点和边界LINEl、LINE2、LINE3和LINE4上施加不同的边界条件。I、Ⅱ和Ⅲ区代表卸载区、减压区和裂前区：其中，点A、B和边界LINEl、LINE2、LINE3是约束

边界条件，C点和边界LINE4为运动边界条件。I、Ⅱ及Ⅲ区随时间而变化。

根据对称性指定边界条件如下：
B、F点是固定点；
C点随时间而变化，AC是自由边；
边界LINE1（包括A点）与LINE3上约束了z向的位移和绕x、y轴的转动；

边界LINE2上约束x、y向的位移和绕x、z轴的转动；
边界LINE4上约束x向的位移和绕y、z轴的转动。
2.3.2表征裂纹驱动力的G与CTOA的计算
本节采用节点力释放法计算裂纹驱动力G。根据（2-38）式，以钢制管道为例，指定设计参数：直径为0.72m，壁厚0.015m，管长取37.5m，压力为0.30MPa,裂纹稳定扩展速度为550m/s。图2-8给出了程序计算的结果，反映了裂纹驱动力与裂纹扩展距离的变化关系。在这一计算中，裂纹后面压力衰减长度，即（2-47)式中的L，取为2倍的直径长度。
如图2-8，随着裂纹起裂后的迅速扩展，裂纹驱动力上升直至达到并保持在稳定的数值，这一平台值1.6KN/m就是裂纹稳态（steady state）扩展的驱动力G。裂纹尖端的移动过程中，G在该值附近有一些小的波动，这是由于裂纹长度的突然增加和位移约束的突然解除，引起的有限元解的高频振荡。振荡的幅度与有限元单元网格的疏密和单元内节点力释放的指定规律有关。

图2-9表示了CTKOA的计算方法。由于在裂纹后方有很大的非线性变形，计算中张开位移用割线近似，采用8个单元的长度来计算CTOA。这一长度经过与细化后的网格结果比较，误差在10%之内。

2.4管道稳态裂纹扩展算例分析

如2.1.2节所述，可在Kanninen（1980）按一维梁模型推导的输气管线裂纹扩展极限速度（2-6a）式的基础上，对高延性管道进行修正。按照西气东输管道选用的设计参数，可得到大致的延生断裂稳态裂纹扩展速度，估算裂纹稳态扩展达到该速度时的裂纹驱动力G和CTOA。另我，通过选用不同的设计参数，可给出稳态裂纹驱动力随内压、壁厚、管径以及径厚比等参数变化的规律。

2.4.1西气东输设计参数下的典型算例

在这一小节，我们研究的管道参数限于：管道长度L=25m，直径D=l.016m,内压10.OMPa，厚度h=0.0147m。网格划分比2.3.1节中的网格两个方向各密一倍，即：沿轴向划分为500个单元，环向划分为32个单元，共16000个单元，16533个节点，每个四边形单元轴向长度0.05m，环向长度0.05m。裂纹稳态扩展速度按（2-6b）式，取为216m/s。轴向压力衰减长度取1.5倍直径，即约1.5m。

首先观察裂尖附近的MISES应力分布，如图2-10。按以针状铁素体和块状铁素体组织为主的X70钢级管线钢的平均屈服强度500～600MPa估算，计算给出的裂尖前部塑性区长度大约在3个单元长度以内，即0.15m左右。

图2-11直观地给出了内压10MPa，裂纹速度216m/s时一维管道气体压力模型（L=1.5D）下的管壁压力沿轴向的分布。每条曲线的衰减起始点即为该时刻的裂尖位置。图中裂尖压力约为6.6MPa。

图2-12展示了管道裂纹稳态扩展中的开裂段位移的分布。图中最大位移点位于裂尖后部第28个单元，即1.4m左右，位移量为0.23m。裂纹后方远端的管道因回弹现象而位移量降低，这与全尺寸实验中观测到的现象是吻合的。

图2-13分别按照三个方向上的位移分布在开裂前管道表面画出云纹图。从图中可以看到，此种设计参数下的位移量主要发生在裂尖附近5m左右的区域，对远端的影响较小。

在裂纹扩展与止裂分析中，代表裂纹驱动力的G和CTOA是最受关注的量。下面分别对本小节稳态扩展情况下的典型算例进行分析。

图2-14反映了节点力释放法计算的动态能量释放率与裂尖位置的关系，其结果因前面提到过的有限元算法的原因而略有抖动，G的稳态值约35KN/m，远高于X70钢材的断裂韧性G_d<10KN/m。

CTOA的计算值大约在14°左右，略高于第四章中测到的（CTOA）c≈12°。

2.4.2不同设计参数对裂纹驱动力的影响

在这一小节，我们分别根据与西气东输有关的不同的管道设计参数，在裂纹稳态扩展的前提下，对G和CTOA的变化作规律性的评测。

目前我国输气管道的最高输气压力为6.4MPa；俄罗斯运营的输气干线最高为7.5MPa；西气东输工程先期设计压力8.4MPa；干线设计压力10MPa；国际上的输送压力在二十世纪末期达到了14MPa。图2-16和2-17给出了上述几种压力下计算的动态能量释放率G和裂纹尖端张开角CTOA的值，管道口径和壁厚均分别为1.016m和14.7mm。从图上看，G和CTOA与P_O呈线性关系。

西气东输工程中，靖边——上海复线管壁厚11.7mm；根据地理结构与危险程度划分的一、二、三、四类地区分别采用的干线臂厚为14.7、17.5、21.0、26.2mm。我国现有输气管道最大直径 660mm（陕京线）；西气东输管线天津——南京支线口径813mm；主要干线口径1016mm；先期设计口径1118mm；预计未来十年间进口管线直径将达1422mm。

根据式（2-15）和（2-16），引入管道无量纲参数SDR，即直径与厚度之比，固定其余参数，分别按上述不同工程参数变化壁厚和直径，总结出G和CTOA随SDR变化的规律。从图2-18和图2-19看，除个别点外，G和CTOA与SDR满足线性关系。但有一点值得注意，这里未考虑SDE对稳态扩展速度的影响。

2.5本章小结

在本章中，我们将断裂动力学中的基本原理（如裂纹极限速度、动态能量释放率、止裂动态观点与静态观点等），基本方法（节点力释放法）和基本止裂判据与壳体动力学有限元方法结合起来，利用气体压力模式确定管壁处的气体压力，形成了求解输气管道上稳态裂纹扩展的数值分析方法。
程序计算各种设计参数下裂纹驱动力的结果同管材实验得到的韧性值处于同一量级，基本反映了目前求解的工程问题所处的状态。通过计算得到的管道稳态扩展过程中的各向应力与位移的具体数值，对分析管道变形量和强度要求有重要的参考意义。
在裂纹稳态扩展的条件下，动态能量释放率G或裂纹尖端张开角CT（刘同管道内压p₀或径厚比SDR之间均呈线性关系。由此关系可便捷地推算任意设计参数下的裂纹驱动力。通过和第四章介绍的小试件实验的动态断裂韧性结果相对照，可以得到裂纹止裂或扩展的结论。
本章的结果还存在一些不足。
一是与计算结果密切相关的稳态裂纹扩展速度难以确定。本文拟定的（2-6b)式仅仅是一个粗略的估算公式，无法保证计算的精度。而测定裂纹扩展速度的全尺寸爆破实验因为其昂贵（一次实验耗费上千万人民币）和测量手段的高要求，目前国内尚未进行；
二是对于止裂过程的描述。高韧性输气管道上已有的全尺寸爆破实验结果表明，动态裂纹的止裂大都不是瞬间完成的，而是有一段明显的减速过程。本章仅能提供稳态裂纹扩展的驱动力，通过和小试件实验测得的动态断裂韧性相比较得到止裂或扩展的结论，却无法提供止裂过程、止裂点以及止裂长度等止裂分析中的关键量；
三是气体压力模式的取值范围与可靠性难以保证。尽管该模式早在三十年前就得到了实验的验证，但仅常见的衰减模式就有指数衰减的Fourie：级数和抛物线衰减函数等一系列表达，本文应用的线性衰减模式中的衰减长度也有不同的取值方法。尤其是高压输送导致欠膨胀超声速射流出现以后，气体压力模式的适用性也在受到挑战。
本文将在下面的三至五章中对上述问题作进一步的研究与分析。