基于四参数棘轮理论的压力容器安定分析

鲍宏琛，沈 鋆，刘应华

(清华大学 航天航空学院，北京 100084)

0 引言

压力容器是过程工业中的关键设备，广泛应用于核能、石油、化工、电力、制药等领域。随着工业装备向着大型化、高参数发展，很多结构单元和组件服役的条件(比如面临高温、高压、腐蚀等环境)越来越苛刻，并且载荷时常是循环反复作用。为确保结构免于交替塑性或增量塑性垮塌，需要确定结构在循环载荷作用下的安定域，从而为工程设计和安全评定提供准确可靠的依据。

经典的Bree图[1]常用来展示理想弹塑性结构在循环热机械载荷下不同的力学响应行为，同时也是ASME规范中承压设备安定性设计准则的理论基础。Bree问题考虑了一个轴对称圆柱壳模型在恒定的一次膜应力和循环的热弯曲应力共同作用下的安定和棘轮行为。它的基本假设是应力作用截面不允许转动，且二次应力中只考虑了径向热应力。KALNINS[2]通过有限元弹塑性分析研究径向或轴向热梯度及内压相互作用的壳体棘轮行为时发现，当热应力中热薄膜应力占主导时，Bree棘轮边界以及3S准则存在不保守的情形。2008年，REINHARDT[3]在安定分析中考虑了恒定机械膜和循环热弯曲、热薄膜应力间的相互作用，通过非循环方法推导建立了三参数棘轮边界，基于他的理论推导及简化建议，ASME Ⅷ-2—2013[4]对热应力棘轮校核准则做了少量修改，增加了对热薄膜应力范围的限制;同时简化弹塑性分析中的“热弯曲应力”也被修改为“热应力”。然而REINHARDT[3]的三参数棘轮边界没有考虑一次弯曲应力与其他三种类型应力间的相互作用。2018年,SHEN等[5]打破了截面平移变形这一基本假设，在棘轮分析中考虑了机械弯曲应力的影响，应用非循环方法,首次推导建立了四参数棘轮边界理论。

从以上研究可以发现，目前安定和棘轮分析中考虑的应力类型逐渐增多，安定和棘轮边界控制方程也越来越复杂，必须进行简化，以适合工程应用的需要。上述推导仅基于一类载荷形式，即机械载荷保持恒定、热载荷不断循环。对于一般载荷工况下，考虑两种应力作用的Bree类修正问题，目前已经报道了很多理论和数值研究工作。如BRADFORD[6-7]理论推导了一次载荷和二次载荷同相循环及异相循环下的修正Bree图;PENG等[8]数值研究了三类载荷相位下的Bree问题。然而对复杂载荷工况下,考虑三种或四种应力作用的Bree类修正问题，目前相关研究还很少见。实际上以应力分类法为基础的压力容器分析设计方法考虑了多种类型应力，因此有必要对同时考虑广义载荷工况和多种类型应力作用的Bree类修正问题开展进一步研究。

本文的目标是基于四参数棘轮理论,开展复杂热机械载荷工况下的压力容器安定分析。由于复杂载荷形式下的理论推导非常困难，且对于实际运行的工业设备，一般不能随意改变载荷条件，试验研究非常昂贵，而数值分析方法可以考虑任意载荷变化，快速构建安定和棘轮边界，因此是一种比较可取的研究方法。本文采用高效准确的线性匹配方法，结合一种新颖灵活的双平面模型开展严格安定分析，考虑机械膜、热薄膜、热弯曲三种类型应力间的相互作用，基于四参数棘轮理论，给出三类工况下的安定域及对应的参数方程。

1 基于线性匹配法的安定分析数值格式

线性匹配法(LMM)用于求解各种复杂结构在任意热机械载荷组合下的非线性响应，该方法采用一系列修正弹性模量的线弹性分析来模拟结构的弹塑性力学行为，是一种快速而直接的方法。线性匹配法的鲁棒性和适用性已经通过有限元逐步非弹性分析得到了验证。通过对一系列经典问题和复杂工程实例的求解，该方法被证明能准确快速地计算得到安定极限和棘轮极限[9-10]。以下简要介绍严格安定分析的数值格式，该算法已在ABAQUS分析软件中使用用户子程序实现。

1.1 循环塑性问题

对于一个给定的结构体，其体积为V，表面积为S,假设其在体积V内承受变化的热载荷为λθθ(x,t) ，在表面St上承受一个变化的机械载荷λPP(x,t)，在剩余表面Su(Su=S-St)上满足零位移边界条件，且热载荷和机械载荷的变化具有相同的周期T。对于一个周期内的时间历程0≤t≤Δt，结构体的线弹性应力解可以表示为：

(1)

假设结构材料满足 Drucker条件，在循环载荷作用下，一个典型时间周期内的应力应变率逐渐达到稳定循环状态，即：

(2)

对于任意的循环载荷历史，其应力解σij(x,t)可表示为:

(3)

(4)

1.2 安定分析的全局最小化过程

本文选用的LMM安定分析格式基于Koiter上限安定定理，安定极限乘子的求解涉及到一个能量的最小化过程，其增量格式可以表示为：

(5)

(6)

(7)

(8)

一系列线性增量关系可以定义为：

(9)

式中，上角标(′)为偏量。

将一个周期内的线性增量表达式相加，可得：

(10)

1.3 安定极限乘子

一系列单调递减的上限安定极限乘子可由式(11)计算得到，最终逼近到真实的安定极限乘子：

(11)

2 双平面模型

为了得到结构在机械膜应力、热弯曲应力、热薄膜应力三种类型应力相互作用下的各类三维弹性安定边界，本文采用文献[5]提出的双平面模型，机械弯曲应力在本文研究中假设为零，其模型如图1所示。

图1 双平面模型示意

双平面模型由两个相同的平面应力模型组成，两个面通过内部靠近的两个边刚性耦合在一起，两个外部边界的Y方向位移进行约束。由于四参数棘轮理论打破了截面平移变形这一假设，因此截面的弯曲应变不被限制。图1的模型中，均匀的一次膜应力可以通过直接将机械载荷施加在耦合边界上实现，由于耦合条件，机械载荷由两个平面共同承受。均匀分布的热薄膜应力和线性分布的热弯曲应力可以通过设计双平面施加的温度场实现。如图1所示，两个平面在Y方向分别施加线性分布的温度梯度，温度场的梯度分布由不同的色块表示。两个平面具有相同的温度梯度变化范围、但具有不同的平均温度。其中热弯曲应力由温度梯度变化范围控制，模型中的热梯度引起的最大的纯弹性热弯曲应力σsb=EαΔT/2(其中,E为弹性模量，α为热膨胀系数，ΔT为模型施加的温度梯度变化范围)。而热薄膜应力则由两个平面的平均温度引起。假设σsm=Rσsb,其中R≥0定义为二次膜弯比，是本文新引入的一个参数。

应用ABAQUS有限元软件进行分析计算，模型有限元网格及几何尺寸如图2所示。理想弹塑性材料属性见表1，假设材料参数在温度循环中是恒定的。

图2 模型有限元网格和几何尺寸

表1 材料属性

材料模型假设为均匀、各向同性，采用小位移假定、von Mises屈服函数和相关的流动法则。模型共包含180个平面应力单元，单元类型为CPS8(ABAQUS),每个单元有3×3个高斯积分点。根据单轴应力状态下三种类型应力的线弹性叠加分析可知，模型中的最大应力会出现在结构边缘处，由于交替塑性边界的准确确定对网格尺寸比较敏感，因此在模型边缘处采用了更加稠密的网格。文中的计算结果表明，采用该网格可以得到足够精确的安定极限乘子。

3 安定分析算例

由于四参数棘轮理论中包含了4个变量，对应于4种类型的应力，即机械膜、机械弯、热薄膜、热弯曲，因此棘轮边界控制方程非常复杂，在工程中难以应用。事实上REINHARDT[3]推导的三参数棘轮边界由5个复杂方程控制，在引入ASME规范中的热应力棘轮评定条例时进行了简化。因此，四参数棘轮理论必须进行简化才能用于工程评估。此外，由于四参数棘轮边界只能在三维坐标系中表示，因此需要任意给定一个变量。为尝试对四参数棘轮边界进行初步简化，本文给定机械弯变量，且假设机械弯曲应力为0。本文研究的3种加载路径如图3所示。

图3 安定分析的三类加载工况

3.1 第一类载荷工况

对于第一类载荷工况，热载荷是循环的，机械载荷保持恒定，如图3(a)所示。注意恒定的机械载荷会引起恒定的机械膜应力，而循环的热弯曲应力和循环的热薄膜应力均由循环的热载荷引起。为定量分析热应力中热薄膜应力和热弯曲应力各自所占比例及区分两者对安定边界的影响程度，本文引入了“二次膜弯比”的概念，表示为R，计算式为：

(12)

式中，Δσsm为热薄膜应力范围；Δσsb为热弯曲应力范围。

二次当量热应力范围可表示为：Δσsmb=Δσsm+Δσsb。

本文选取了R的一系列值进行计算分析，即R=0，1/10，1/4，3/7，2/3，4/5，1，3/2，7/3，4，9，+∞(即1/R=0，表示纯热膜) ,研究不同热应力分布下安定域的变化。计算结果如图4所示，图中σpm为机械薄膜应力。定义X=σpm/σy，Y=Δσsb/σy，Z=Δσsm/σy，Y′=Δσsmb/σy。

从图4中可以看到，随着R值的不断增大，安定区逐渐减小。除纯热膜情况外，弹性安定边界均由两段直线组成，其中一段为交替塑性边界；另一段为棘轮边界。当R=0时，该问题退化为经典Bree问题，安定边界控制方程如下：

(13)

图4 LMM计算得到的第一类工况下的安定极限

对于不同的R值，交替塑性边界始终保持不变，且安定域在XY′坐标系中始终过2个固定点(0,2)和(1,0)。因此，一旦确定交替塑性边界与棘轮边界的交点，即可以快速构建弹性安定域。根据四参数棘轮边界理论[5]，当机械弯曲应力为零时，恒定机械膜应力与循环热应力作用下对应的棘轮边界段可以表示如下。

(1)对于Δσsm+Δσsb≤2σy,且Δσsm≥Δσsb。

(14)

其中:

(15)

(16)

(2)对于Δσsm+Δσsb≤2σy,且Δσsm<Δσsb。

(17)

(18)

假设Y′=2 ,则Δσsm+Δσsb=2σy，由式(12)和式(14)～(16)可以解得安定边界拐点横坐标为：

(19)

通过数学计算和数值拟合，可由式(17)(18)得到安定边界拐点横坐标为：

(20)

因此，完整的弹性安定边界可以用参数方程表示如下。

(21)

(22)

对于R>1:

(23)

图5选取几个典型的R值对比了数值计算结果、参数方程解及四参数棘轮理论解，可以看出三者吻合程度较好，证明了该方法的正确性，同时该参数方程解是安全保守的。

图5 第一类工况下的数值解、参数方程解及理论解之间的对比

图6示出第一类工况下的完整的弹性安定域，借助安定边界和参数方程，可以为工程设计和安全评定提供指导。

图6 第一类工况下的弹性安定域

3.2 第二类载荷工况

对于第二类载荷工况，二次载荷与一次载荷呈比例变化，即同相循环，如图3(b)所示。当R=0时，本节考虑的问题退化为Bradford问题[6]。Bradford安定边界[6]控制方程为：

X′+Y=2 (0≤X′≤1)

(24)

其中,X′=Δσpm/σy,Δσpm为一次薄膜应力范围。

本节考虑热弯曲应力、热薄膜应力和机械膜应力三者均同相循环下的弹性安定行为，即R≥0。数值研究发现，当R=0 时，计算得到的安定极限与Bradford理论解[6]吻合得比较好。然而与第一类问题不同的是，当任意改变R值的大小，安定边界始终保持不变。因此,可以推广得到第二类工况下的弹性安定边界控制方程如下。

对于R≥0 ：

X′+Y′=2 (0≤X′≤1)

(25)

注意式(25)中纵坐标变量Y′的变化。计算结果与Bradford安定边界的对比如图7所示。完整的弹性安定域如图8所示。

图7 第二类工况下的安定极限与Bradford安定边界的对比

图8 第二类工况下的三维弹性安定域

3.3 第三类载荷工况

对于第三类载荷工况，热载荷和机械载荷可以各自独立变化，如图3(c)所示。LMM计算得到的数值结果如图9所示。

图9 LMM计算得到的第三类工况下的安定极限

从图9可以看出，随R值不断变大，安定域逐渐减小。除去纯热膜的情形，安定边界均由两段组成，其中交替塑性边界与第二类工况下的相同，棘轮边界与第一类工况相应的边界段重合。实际上汇总对比三类安定问题可以发现，对于任意选取的R值，第三类载荷工况下的安定边界始终是第一类和第二类载荷工况下的安定边界的下包络线。通过前两类安定问题的参数方程，可以推导得到第三类工况下弹性安定边界控制方程如下。

(26)

(27)

对于R>1:

(28)

其中[X′]由以下公式确定：

(29)

(30)

图10选取了典型的R值对比了三类工况下的计算结果及第三类工况的参数方程解。可以看出，计算结果与参数方程解吻合得比较好。对于第二类载荷工况，安定边界与ASME等规范中用于弹性安定分析的3S准则重合，3S准则用于弹性安定评定是保守的。对于第一类和第三类载荷工况，3S准则仅与部分安定边界(交替塑性边界)重合，此时3S准则是安定评定的必要不充分条件，需要按实际安定边界公式进行补充评定。因此3S准则用于弹性安定分析需要谨慎使用，且需注意载荷工况类型。图11示出第三类工况下完整的三维弹性安定域。

图10 三类安定问题安定域对比(R=0.8)

图11 第三类工况下的三维弹性安定域

4 结论

本文将线性匹配法与双平面模型相结合，开展了三类载荷工况下压力容器的严格安定分析研究，基于四参数棘轮理论，考虑了三种常见类型应力的相互作用,即一次膜应力、二次热弯曲和热薄膜应力。通过引入一个新参数“二次膜弯比”R，揭示了三类安定问题之间的联系，对四参数棘轮理论进行了简化和推广，首次给出了三类安定问题完整的弹性安定域及参数方程，可直接用于工程设计和安全评定。

通过研究发现，对于第一类和第三类加载工况，随热应力中热薄膜比例增大，安定区均逐渐减小。这一点与考虑截面平移假设时的情形(如Reinhardt的工作[3])不同，这说明对于本文研究的算例，热薄膜应力要比热弯曲应力危险。因此，在工程应用中要注意区分二次应力的影响，同时注意安全校核准则的适用范围。对于第二类加载工况，安定边界在不同热应力分布下保持不变。对任意的二次膜弯比，第三类工况下的安定边界是第一类和第二类载荷工况下相应安定边界的下包络线。

除机械载荷和热载荷同相循环外，ASME等规范中给出的3S准则用于弹性安定分析和评定存在不保守性，需要谨慎使用。基于安定边界参数方程，可确定3S准则的适用范围和补充评定条件。本文提出的分析方法和简化思想可用于复杂热机械载荷工况下考虑多类型应力作用的压力容器安定和棘轮问题研究。