吴济航 苏文强 顾亚楠 缪 好 陈志平

（浙江大学化工过程机械研究所）

大型石油储罐为典型的薄壁圆柱壳结构，受沉降与风载荷联合作用而破坏的事例屡见不鲜［1～3］，对此学者们开展了大量分析研究。 陈志平等针对储罐基础沉降，提出将基础沉降量作为边界条件的建模方法，可大幅提高罐体应力计算的精度［4］。范海贵等基于实测沉降数据，运用理论推导与数值模拟相结合的方法研究了浮顶储罐的变形，同时提出一种泊松曲线模型来预测未来某时刻储罐的沉降［5，6］。 在沉降预测方面，程健等提出一种指数曲线预测方法［7］。 石磊等比较分析国内外储罐沉降的标准， 认为标准要求都较为保守，应利用有限元方法充分考虑储罐对地基沉降的结构响应，同时还分析了地基沉降与动水压力对储罐屈曲强度的影响［8，9］。 Nassernia S等应用有限元方法对储罐沉降的实验和评判方法进行验证，发现有较好的一致性［10，11］。 李玉坤等利用ANSYS软件对拱顶储罐在静载、动载和风载作用下的网壳顶结构强度进行分析，并提出改进措施以加强其承载能力［12］。 Zhao Y等研究了不同储液高度对圆柱形储罐抗风性能的影响，发现随着储液高度的增加，储罐的临界屈曲风载荷随之增加［13］。 Chiang Y C和Guzey S考虑储罐缺陷的影响，对承受动风载荷的敞口空罐进行几何非线性显式动态分析，研究其振动特性［14］。 苏文强等对不均匀沉降与静风载荷联合作用下的储罐进行屈曲分析，发现沉降会削弱储罐抵抗静风屈曲的能力［15］。

现有的研究主要集中于储罐的静风屈曲分析，即风载荷与时间无关，而实际风载荷的风压是随时间变化的。 虽然储罐静风屈曲分析采用的是某一地区的极限风压值，但由于未考虑风载荷作用时间对储罐屈曲的影响，故会导致静风屈曲分析结果表明不会发生失效的储罐，在实际风载荷作用下也有可能发生破坏，因而有必要对储罐开展动风屈曲研究。 在此，笔者提出一种不均匀沉降下储罐动风屈曲的有限元分析方法，并利用有限元软件ABAQUS对一在役储罐进行案例分析。

1 不均匀沉降下储罐的动风屈曲分析方法

在不均匀沉降下储罐的动风屈曲分析中，储罐受到不均匀沉降和动风载荷的作用，属于轴向静载与径向动载的联合作用。 在该联合载荷作用下，储罐的屈曲分析是一个复杂的非线性动力学分析过程，为此在ABAQUS/Explicit模块下应用显式动力学方法， 设置两个分析步对其进行分析。由于该模块下无法设置静力学分析步，因此在第1个分析步中应用ABAQUS/Explicit准静态方法对储罐在沉降作用下的结构响应进行分析，在此基础上第2个分析步应用显式动力学方法分析储罐的动风屈曲。

1.1 ABAQUS/Explicit准静态方法

ABAQUS/Explicit准静态方法实质上是一种动态求解过程，其求解原理是通过人为控制加载速率使系统惯性力被忽略，从而得到与静力分析相近的结果。 储罐沉降的加载速率由加载时间控制，加载时间越长，加载越慢，分析也越接近静力分析。 加载时间的选取既要考虑计算效率，又要满足能够忽略系统惯性力的要求。 文献［16］认为，当加载时间取一阶固有周期的5倍左右时，惯性力的影响能够降低到可接受范围内。 笔者据此初步估算出一个沉降加载时间值，在该值附近选择3个加载时间分别进行储罐的沉降模拟， 比较其计算效率、惯性力影响程度及其与静力学的计算误差，然后选择合适的加载时间。 通常，若分析过程符合以下两条准则，则可以认为系统惯性力影响不显著［16］：

a. 外力做功EW几乎全部转化为系统内能EI，即两者的绝对值几乎相等；

b. 变形材料动能EK波动不大且值较小，不超过内能EI的5%～10%。

1.2 沉降罐动风屈曲分析方法

在第2个分析步中， 将动风载荷加载到储罐外壁面上，应用显式动力学方法进行计算，目标是得到储罐径向位移云图，提取罐壁上径向内位移最大点的径向位移时程曲线，曲线上位移突变的时刻即为屈曲发生时刻。 另外，将沉降罐的动风屈曲结果与静风屈曲结果作比较，以研究风载荷作用时间对储罐屈曲的影响。 最后，为了分析不均匀沉降分布形式对储罐动风屈曲的影响，可改变储罐的迎风受压面，每隔一定角度将动风模型加载到储罐模型上，然后比较动风载荷以不同角度作用于储罐上时储罐发生屈曲的时间。

2 案例分析

为方便与文献［15］中沉降罐的静风屈曲作比较，笔者选择该文献中舟山地区的20 000 m3固定顶储罐进行不均匀沉降下的动风屈曲分析。

2.1 储罐模型

储罐内径42 m，高17 m，罐壁为变壁厚形式，共分10层，自下向上厚度逐渐减小。 罐壁每层的高度h和壁厚δ见表1。 罐壁材料选用Q345R，弹性模量E=206 GPa，泊松比ν=0.3，密度ρ=7850 kg/m3，材料非线性参数［1］见表2。

表1 罐壁每层的高度h和壁厚δ

表2 材料非线性参数

同时，为了便于分析，对储罐模型做如下简化：

a. 储罐内介质可提高储罐的抗风稳定性，故仅考虑空罐工况，并作为最危险工况，而不考虑储罐内部液体静压力的影响；

b. 固定顶约束罐壁顶部的径向位移，储罐底板约束罐壁底部的径向位移和周向位移，故直接在模型中约束这两处的位移，而不再建立储罐的固定顶模型和底板模型；

c. 不考虑其他附属部件（如接管）的影响。

在ABAQUS软件中按照储罐的实际尺寸结构建立几何模型，如图1所示。

图1 储罐几何模型

采用S4R壳单元对模型进行网格划分， 既能降低计算成本，又能使计算精度得到保证。 网格单元数为165 432，节点数为160 026，经验证满足网格无关性要求。

2.2 边界条件与载荷模型

2.2.1 储罐沉降边界条件

表3列出了储罐罐壁底部的12个沉降实测值S0～S11。

表3 储罐罐壁底部沉降实测值mm

利用文献［17］中提出的沉降函数拟合方法对其进行拟合，提取不均匀沉降表达式：

其中，SU（θ）表示周向角θ处的不均匀沉降值。

2.2.2 储罐动风模型

2.2.2.1 风载荷分布模型

储罐沿周向的风载荷计算式为：

其中，风压幅值λ可根据GB 50009—2012《建筑结构荷载规范》由基本风压w0计算得到。 查得舟山地区基本风压w0=850 Pa，忽略风压幅值沿储罐高度方向的变化，取罐顶处的风压幅值保守计算得到λ=1985.6 Pa。

表4 风载荷体型系数中的傅里叶系数

综上可得风载荷从0°位置作用于储罐时罐壁表面的周向风压分布为：

2.2.2.2 风载荷时程变化模型

在工程核算过程中，风载荷时程模型通常采用各种简化形式，在此笔者选用方波形式来简化动风载荷模型［19］，如图2所示。 由图2可以看出，风载荷在0.05 s内完全加载，持续作用至6 s后结束［20］。

图2 风载荷与时间的关系曲线

2.3 储罐模态分析

模态分析是动力学分析的基础，可获得结构固有频率以避免共振，同时也为后续准静态分析做准备。 表5列出储罐前6阶固有频率。

表5 储罐前6阶固有频率

2.4 储罐沉降结构响应的准静态分析

根据储罐模态分析中的一阶固有频率得到一阶固有周期为0.257 3 s， 初步估算加载时间为1.28 s。 现在1.28 s附近选取3个加载时间1.00、1.25、1.50 s， 分别对沉降罐模型的结构响应做显式准静态分析，绘制动能EK、内能EI和外力做功EW的时程变化曲线（图3）。

图3 不同加载时间下准静态分析的能量和动能时程变化曲线

由3条能量时程变化曲线可以看出， 内能与外力做功曲线基本重合，动能曲线与水平轴基本重合，说明外力做功基本转化为内能，动能与内能相比极小；由3条动能时程变化曲线可以看出，动能均只有两个波峰，没有反复波动，且加载时间越长，最大动能越小。 由此说明3个加载时间下系统惯性力影响均较小。

再对储罐进行不均匀沉降作用下的静力分析，与上述3个准静态分析作比较，分别计算罐壁上最大与最小径向位移的相对误差eΔ1、eΔ2后，与动能、内能之比（EK/EI）一并列入表6。 另外，为了分析3个准静态分析的计算效率，以CPU相对耗时（无量纲）来比较它们的相对计算时间，即当加载时间1.50 s下的实际计算时间为1时，其他情况下的实际计算时间则会相对小于1，将CPU相对耗时也列入表6。 比较表6中的数据可以发现，加载时间为1.00 s时， 最小径向位移的相对误差eΔ2最小仅为0.20%， 最大径向位移的相对误差eΔ1最大但不超过5.00%，EK/EI为1.01%也不超过5.00%，此时CPU相对耗时仅为0.52 s， 约为加载时间1.50 s的一半，极大地提高了计算效率，故最终选定准静态分析的加载时间为1.00 s。

表6 准静态分析各项指标比较

2.5 沉降罐动风屈曲分析

在上述分析的基础上，每隔15°在储罐上加载动风模型进行动风屈曲分析。 首先，比较不同风向角α下储罐发生屈曲的时间长短， 以分析不均匀沉降对储罐动风屈曲的影响。 另外，再将对应的静风屈曲结果作对比，以分析风载荷作用时间对沉降罐屈曲的影响。

2.5.1 沉降罐动风屈曲结果

以风向角α分别为0、120、240、330°为例，给出储罐在动风载荷作用6 s后的径向位移云图（图4，迎风受压面朝外）。

图4 沉降罐在动风载荷作用6 s后的径向位移云图

由图4a、b、d可以看出，罐壁的迎风受压面均产生了明显的凹坑，说明储罐在此方向上受到正向风压作用后容易产生大变形，造成局部的屈曲失稳。例如图4a中，动风载荷以0°风向角作用于储罐时，6 s 后罐壁上向内的最大径向位移有334.0 mm。 但在图4c中，当风向角α为240°时罐壁的迎风受压面并未出现明显凹坑，其径向变形主要是由不均匀沉降产生的，说明从该方向加载动风载荷时对罐壁的影响较小。

为了得到储罐发生屈曲的时间t′， 可提取6 s时罐壁上径向向内位移最大的点，绘制其径向位移r的时程变化曲线。 图5为α分别为0、15、120、330°时的r-t′曲线，以α=120°为例，径向位移r先在10.0 mm位置上下波动， 在C点位置迅速增大至241.6 mm，随后在该位置上下波动，故可将C点判定为储罐发生屈曲失稳的点，从而得到发生屈曲的时间t′=1.18 s。

图5 沉降罐在不同α的动风载荷作用下径向位移最大点的位移时程变化曲线

2.5.2 不均匀沉降对储罐动风屈曲的影响

不均匀沉降对储罐动风屈曲的影响可根据储罐在动风载荷作用下发生屈曲的时间t′来判断，t′越小表示储罐越容易发生屈曲， 该位置处的不均匀沉降对其抗风稳定性的削弱也越大。表7列出了沉降罐在不同风向角α的动风载荷作用下发生屈曲的时间t′， 若未标注屈曲时间，则表示从该位置作用动风载荷时， 储罐不会发生屈曲。 从表7中可以看出，储罐在300～360°区域向上隆起， 该区域储罐发生屈曲所用时间都较短，在0～60°区域向下沉降，该区域的屈曲时间普遍较长，而0～15°发生屈曲时间也较短，可能是因为该处还在提离影响区域的范围内， 因此可以判定沉降罐基础向上隆起区域比向下沉降区域的抗风稳定性弱。 而其他位置在6 s内并未发生屈曲， 或即使在120～135°的向上提离区发生屈曲所用时间较长， 但其原因可能是沉降幅值较小， 并未使罐壁产生过大的径向变形，由此说明沉降幅值越大， 储罐的抗风稳定性越差。

表7 沉降罐在不同风向角的动风载荷作用下发生屈曲的时间和静风屈曲临界载荷因子

2.5.3 沉降罐动风屈曲与静风屈曲的比较

以0°位置受静风作用的储罐为例，LPFwc=0.323表示当风压值达到基本风压下风压幅值的0.323倍，即641.3 Pa时，储罐发生屈曲失稳。表7中某些角度下未列出LPFwc的值，表示该方向的风载荷以基本风压作用于储罐时没有发生屈曲。

分析表7中数据可知， 当在30～60°、120～135°范围内时， 储罐在静风载荷作用下未发生屈曲，但在动风载荷作用下却发生了屈曲，由此说明不均匀沉降下储罐抵抗动风载荷的能力弱于静风载荷，因此仅对沉降罐开展静风屈曲分析不能完全保证储罐的安全性，还需考虑风载荷作用的时间效应，对沉降罐开展动风屈曲分析。

3 结论

3.1 应用ABAQUS/Explicit准静态方法能有效模拟储罐在不均匀沉降作用下的静力响应，在此基础上可顺利进行动风屈曲分析。

3.2 储罐基础向上隆起区域比向下沉降区域对储罐抗风稳定性的削弱更大， 且沉降幅值越大，削弱作用越大。

3.3 仅对沉降罐开展静风屈曲分析不能完全保证储罐的安全性，还需考虑风压值随时间变化的影响，对沉降罐开展动风屈曲分析。