支柱沉降作用下球形储罐的屈曲研究

2022-03-23易天坤欧晨希林冠堂

重庆理工大学学报(自然科学) 2022年2期

黄思，易天坤，欧晨希，林冠堂

(1.华南理工大学机械与汽车工程学院，广州 510641；2.广东省特种设备检测研究院珠海检测院，广东珠海 519002)

球形储罐被广泛应用于石油、化工、冶金、轻纺、城建、市政等工业领域中，主要用于储存液态丙烷、丁烷、丙烯、丁烯及其混合物等[1]。当结构在微小扰动(载荷变化)使结构发生较大的形态变化，则称为结构发生屈曲或失稳，通常结构屈曲的发生表现为变形的突然增大。由于受到基础质量、安装施工、地理环境等因素的影响[2]，很难使球罐支柱保持在同一水平面上，支柱一旦沉降，球罐受力不再均匀，球罐局部受力过大可能会屈曲或破损。

国内外学者相继开展了储罐在沉降下的静力学研究和屈曲研究。静力学研究方面，Dorazio等[3]通过实验对圆柱储罐在不均匀沉降下的变形和应力进行了简单的测量；Yang等[4]采用有限元法计算了支柱不均匀沉降下的球罐应力和支柱位移；高红利等[5-6]研究了基础不均匀沉降对直径为9.2 m的球罐应力的影响；房江祥[7]以1.5 m3球罐为研究对象，运用ABAQUS软件来模拟球罐支柱沉降对其拉杆应力和极限载荷的影响。Chen等[8]通过ANSYS软件建立了圆柱储罐的有限元模型，研究了圆柱储罐在不均匀沉降下的变形行为。屈曲研究方面，Jonaidi等[9-10]采用有限元和实验验证的方法，研究了圆柱储罐在理想谐波沉降下的屈曲响应；Godoy等[11-12]采用实验方法研究了局部沉降下圆柱钢储罐的屈曲响应；范海贵[13]采用现场实测、理论推导和有限元数值模拟相结合的方法，研究了圆柱储罐在沉降和轴向载荷联合作用下的屈曲响应；苏文强[14]建立了大型圆柱储罐在不均匀沉降和风载荷联合作用下的屈曲数值模拟方法。

综上所述，大多数学者已对圆柱储罐在沉降下的静力学和屈曲变形进行了大量的研究。但针对球罐的研究主要集中于沉降下的静力学研究，缺乏对沉降下的球罐的屈曲变形研究。球罐支柱的沉降量达到一定程度时，球罐可能因局部受力不均而发生屈曲变形和破坏。因此，选取1台5 000 m3的在役丙烷球罐为研究对象，对其支柱在3种不均匀沉降形式下进行屈曲响应研究，为球罐设备的管理和维护及稳定运行具有重要工程意义。

1 计算模型

研究对象为1台5 000 m3的在役丙烷球罐，球罐的结构示意图如图1所示，其结构尺寸的基本参数和主要结构的材料特性参数如表1、2所示。

图1 球罐结构示意图

表1 丙烷球形储罐设计参数

表2 球罐主要结构的材料特性参数

2 屈曲响应分析

2.1 网格划分

采用Workbench平台中的Mesh对球罐进行结构网格划分，为避免网格单元数目对计算结果的影响，需要考察计算域在几种不同网格数量下的计算结果。同时考虑到本文仅对球罐进行屈曲变形的研究，因此以单根支柱沉降1 mm为例，图2给出了通过Workbench平台计算得到的球罐的最大总位移xmax随计算网格数的变化情况。

图2 网格无关性检验曲线

由图2可知，当球罐的网格数达到75万之后，球罐最大总位移xmax值基本稳定。考虑到计算的精度和成本，选用767 830为球罐计算域的网格数量，模型网格划分结果如图3所示。

图3 网格划分结果示意图

2.2 载荷与约束条件

考虑球罐充装率φ为90%的工况[16]，该工况下，球罐的载荷有重力G、内壁承受均布的压力p。其中压力载荷p(Pa)随液面高度h(m)呈线性变化，其变化规律为：

p=po+(R-h)ρlg

(1)

式中：po为球罐的工作压力，为1.42 MPa；R为球罐的半径；ρl为液态丙烷密度，为474.04 kg/m3；g为重力加速度，取9.81 m/s2。

图4给出了支柱编号的示意图，对于球罐的支柱不均匀沉降的模型，考虑了3种沉降形式：单根沉降、蝶形沉降和两根沉降。表3给出了3种沉降方式下每根支柱的沉降量(Sd)。在设置支柱的沉降量时，通过对沉降支柱的底板设置沿着轴向向下的位移，从而实现对球罐不均匀沉降的有限元模拟。

图4 支柱编号示意图

表3 3种沉降方式的沉降量(Sd)

2.3 屈曲分析原理

当加载到结构上的载荷达到某一值时，若增加微小的载荷，而结构的位移却发生很大的改变，通常将这种情况叫做结构屈曲，对应的载荷称为屈曲载荷，这时的状态称为临界状态。根据临界状态的特性，屈曲可分为分支屈曲、极值屈曲和跳跃屈曲3种，分别对应图5中的(a)、(b)和(c)所示的载荷-位移曲线[17]。

图5 3种屈曲特性下的载荷-位移曲线

当球罐支柱的不均匀沉降量达到一定水平时，球罐在自重G、液体压力p和支柱位移沉降量Sd等共同作用下，将发生整体或局部弯曲变形，导致结构屈曲。因此，需要对球罐开展屈曲分析，确定其不均匀沉降量的极限值。针对球罐的屈曲分析可采用“弧长法”。“弧长法”最初由Riks[18]提出，后经Crisfield[19]改进并逐渐得到完善，是目前结构屈曲分析中效率较高的方法。

非线性平衡方程组一般可以表示为：

K·x=F

(2)

其中矩阵K为切向刚度矩阵；向量x为节点位移；向量F为外载荷。现在将方程写成增量形式，并引入一个载荷因子λ，平衡方程可以改写成：

K·Δx=ΔλF

(3)

此时平衡方程共有n+1个未知量：Δx1，Δx1，…，Δxn，λ，而方程数量只有n个，因此需要补充一个约束条件：

Δλ2FTF+ΔxTΔx=Δl2

(4)

弧长法的求解流程如图6所示，通过该流程的求解得到完整的沉降-位移关系曲线。弧长法是利用平衡方程和弧长约束条件寻找平衡点的过程，如图7所示，表3中的球罐支柱的沉降量是逐步加载的，在加载过程中，前期OC段是属于小变形的弹性阶段，CD段及以后为大幅的屈曲变形阶段。当沉降量超过C点所对应的沉降值时，球罐局部将会突然丧失原有的几何形状，发生大幅的屈曲变形，因此，可认为C点后球罐发生了屈曲。本文基于弧长法的数值模拟方法，分别研究球罐在3种沉降方式下的屈曲响应过程，以此来分析球罐变形与不均匀沉降之间的关系。

图6 弧长法迭代求解过程框图

图7 弧长法迭代求解过程示意图

3 结果分析

在应用Workbench平台模拟球罐的屈曲响应过程中，表3中的球罐支柱的沉降量Sd是逐步加载的。载荷因子λ表示球罐屈曲响应过程中支柱不均匀沉降的加载比例，如λ=0.5时，球罐支柱加载的不均匀沉降为0.5Sd。当λ从0增大到1时，意味着球罐支柱的不均匀沉降全部加载完成。

图8是通过Workbench求解得到的在单根沉降下的λ-xmax曲线，其中xmax为最大总位移，λ为载荷因子。图9为该情形下的曲线中A、B、C和D点处的变形云图，其中A、B、C和D点对应的λ值分别为0.15、0.3、0.45、0.55。

图8 单根沉降下的λ-xmax曲线

图9 单根沉降下球罐屈曲响应的变形云图(变形缩放系数=20)

对照图8、5可知，该球罐在单根沉降下的屈曲响应符合典型的分支屈曲的特点，在C点前球罐的变形随支柱沉降的增大呈线性变化，且变形量较小，C点之前并未出现位移突变现象。说明在C点之前，球罐处于平衡状态并未发生屈曲。C点时，球罐的变形云图如图9(c)所示，C点对应的最大总位移为13.5 mm(在支柱1的底部)。当支柱沉降量超过C点时(即λ>0.45)，球罐局部发生屈曲，如图9(d)所示，此时最大位移为36.4 mm，屈曲主要发生在支柱1与支柱2之间的拉杆处。

图10是通过Workbench求解得到的在蝶形沉降下的λ-xmax曲线。图11给出了该情形下的曲线中A、B、C和D点处的变形云图，其中A、B、C和D点对应的λ值分别为0.1、0.2、0.3和0.45。

图10 蝶形沉降下的λ-xmax曲线

图11 蝶形沉降下球罐屈曲响应的变形云图(变形缩放系数=20)

对照图10、5可知，该球罐在蝶形沉降下的屈曲响应也符合典型的分支屈曲的特点，在C点前球罐的变形随支柱沉降的增大呈线性变化，且变形量较小，C点之前并未出现位移突变现象。说明在C点之前，球罐处于平衡状态并未发生屈曲。C点时球罐的变形的云图如图11(c)所示，C点对应的最大总位移为7.5 mm(在支柱4的底部)。当支柱沉降超过C点时(即λ>0.3)，球罐局部发生了屈曲，如图11(d)所示，此时最大总位移为55.8 mm，屈曲主要发生在支柱1与支柱12之间的拉杆处。

图12是通过Workbench求解得到的在两根沉降下的λ-xmax曲线。图13给出了该情形下的曲线中A、B、C和D点处的变形云图，其中A、B、C和D点对应的λ值分别为0.15、0.3、0.45和0.6。

图12 两根沉降下的λ-xmax曲线

图13 两根沉降下球罐屈曲响应的变形云图(变形缩放系数=20)

对照图12、5可知，该球罐在两根沉降下的屈曲响应同样也符合典型的分支屈曲的特点，在C点前，球罐的变形随支柱沉降的增大呈线性变化，且变形量较小，C点之前并未出现位移突变现象。说明在C点之前，球罐处于平衡状态并未发生屈曲。C点时球罐的变形的云图如图13(c)所示，C点对应的最大总位移为10.6 mm(在支柱1与球罐连接处)。当支柱沉降超过C点时(λ>0.45)，球罐局部发生了屈曲，如图13(d)所示，此时最大总位移为44.1 mm，屈曲主要发生在支柱11与支柱12之间的拉杆处。