蒸气云爆炸作用下钢储罐动力响应的双向流固耦合分析

胡 可， 赵 阳， 王 钢， 郭艳军

(1. 浙江大学 空间结构研究中心， 杭州 310058； 2. 四川电力设计咨询有限责任公司， 成都 610094)

钢储罐内部爆炸通常是由罐内所储存的石化物资挥发产生的可燃气体与空气混合后形成的蒸气云遇明火或电火花引发的爆炸事故，其过程主要包括：蒸气云爆炸、爆炸冲击波传播、冲击波与储罐结构相互作用等。相对应的研究包括爆炸冲击波参数计算、储罐内表面爆炸荷载确定以及储罐结构动力响应分析等[1]。其中，结构动力响应分析是钢储罐安全设计与抗爆优化的重要基础。目前，密闭空间内部蒸气云爆炸的CFD (Computational Fluid Dynamics)分析大多仅限于流场的模拟计算，而将密闭容器假定为刚体，未考虑流场与结构的相互作用(即流固耦合效应)[2]。然而，内部爆炸流场的作用会使储罐结构产生运动与变形，这些运动和变形又会反过来作用于内部爆炸流场，进而影响储罐内表面上爆炸荷载的大小与分布[3]。因此，求解钢储罐在内部爆炸作用下的结构动力响应实际上是一个复杂的双向流固耦合问题。

从对控制方程的解法上来区分，流固耦合问题可分为强耦合与弱耦合。强耦合方法把流场、固体结构及耦合作用全部构造在一个控制方程里，在单一时间步内对所有变量同时进行求解，具有物理概念清晰、计算精度高等优点，但由于流固控制方程性质不同，求解难度较大。而弱耦合方法是在每一个时间步内分别对固体结构方程和流体方程进行迭代求解，固体域和流体域的计算结果需要通过搭建中间数据交换平台进行彼此信息的映射来实现耦合求解。弱耦合方法能够使用现有的商业软件或程序，并可以根据实际情况的需要进行固体域和流体域信息的传递，比强耦合方法更加容易实现，近年来得到了广泛的使用[4]。Baum等[5]将弱耦合方法应用到了武器爆炸领域，对爆炸以及残片的运动过程进行了模拟，过程中在流固耦合界面采用了快速插值方法和守恒技术进行数据传递。Chun[6]使用ANSYS和FLUENT模拟了爆炸过程中冲击波对柔性防护结构的影响，研究结果表明，弱耦合方法是分析爆炸荷载下刚体或者柔性结构响应的很好工具。Haupt等[7]采用弱耦合方法对有热效应的流固耦合问题进行了计算，采用的CFD求解器是软件TAU，而CSD求解器则是软件ANSYS，数据之间的传递则采用软件MpCCI。王彬等[8]采用MpCCI平台将有限元软件ANSYS与流体计算软件STAR-CD结合起来，采用弱耦合方法对风与膜结构的耦合作用进行了数值模拟。武岳等[9]采用弱耦合方法针对大跨度张拉结构发展了风振分析的数值风洞方法，在研究中考虑了几何非线性因素，编制了相应的有限元分析程序MFS。薛祖杰[10]采用ANSYS 14.0的Workbench数值仿真平台对超高层建筑的流固耦合风振响应进行了数值模拟，得到了建筑表面风压分布与动向响应。

本文采用弱耦合方法，分别通过软件FLUENT与ANSYS对内部蒸气云爆炸流场和钢储罐结构进行独立求解，并通过ANSYS Workbench仿真平台提供的System Coupling模块进行耦合界面的数据传递，实现蒸气云爆炸作用下钢储罐结构动力响应的双向流固耦合分析。同时进行解耦计算，并将耦合与解耦两种方法进行对比分析。

1 双向流固耦合的基本理论与求解方法

1.1 流固耦合的基本控制方程

流固耦合问题需要求解的耦合方程的定义域包含流体域与固体域。对于钢储罐内部蒸气云爆炸流场来说，其需要遵循质量守恒、动量守恒、能量守恒和化学组分平衡方程。质量守恒即连续性方程，其微分形式为[11]

(1)

式中：ux，uy，uz分别为x，y，z三个方向的速度分量；t为时间；ρ为流体质量密度。

动量守恒的本质是满足牛顿第二定律，式(2)～式(4)为x，y，z三个方向的动量方程

(2)

(3)

(4)

能量守恒的本质是热力学第一定律，表达式为

(5)

组分质量守恒可表示为系统内某化学组分质量对时间的变化率等于通过系统界面净扩散量与化学反应产生的该组分生产率之和，组分方程为[12]

(6)

式中：cm为组分m的体积浓度；ρcm为该组分的质量浓度；Dm为该组分的扩散系数；Sm为生成率。

流体诱发固体振动、位移的控制方程为

(7)

式中：Ms为固体质量矩阵；Cs为阻尼矩阵；Ks为刚度矩阵；r为固体位移；τs为固体受到的应力。

耦合控制方程为

(8)

式中：q为热流量；T为温度；下标f为流体的特性；下标s为固体的特性。

1.2 湍流燃烧模型与数值方法

湍流模型采用在关于湍动能k方程的基础上引入一个关于湍流耗散率ε方程而形成的k-ε两方程模型[13]。该湍流模型的准确性与收敛性都十分符合工程模拟计算的要求，目前被广泛使用。

对于蒸气云的燃烧爆炸过程，软件FLUNET提供了预混燃烧模型、通用有限速率模型以及组分概率密度运输燃烧模型等多种模拟反应的方法。其中，组分概率密度运输燃烧模型通过求解概率密度函数，模拟计算时对反应速率的处理未进行任何假设，严格考虑了化学反应与湍流流动的交互过程[14]。尽管采用此模型的计算量较大、计算时间较长以及对网格质量的要求交高，但从模拟精确性上来讲是高于其他燃烧模型的，因此本文采用组分概率密度运输燃烧模型对钢储罐内部蒸气云爆炸过程进行模拟。

软件FLUENT采用有限体积法，通过求解描述每种组成物质的流动、扩散和反应源的守恒方程来模拟燃烧爆炸。对于压力与速度耦合处理，采用半隐方法来求解压力耦合方程的SIMPLE算法[15]。另外，爆炸流场的湍流流动伴随的化学反应导致模型中质量/动量和物质运输方程强烈耦合，在反应流中获得收敛解比较困难，本文采用亚松弛方法来处理这类耦合问题。

1.3 流固耦合界面的网格映射与数据传递

流固耦合界面是流场与结构的交界面，由于弱耦合算法中的流固耦合作用只发生在此交界面上，所以耦合界面上的数据传递是双向耦合计算的关键。在模拟过程中，结构在流体荷载作用下发生变形或运动，结构需将变形或运动通过耦合界面由CSD (Computational Structural Dynamics)网格传递给CFD网格，然后流场分析软件在新的条件下求解流场，再把新的流体荷载通过耦合界面传递给CSD网格[16]。由耦合控制方程式(8)可以看出，流场和结构在边界上对应的质点之间应该满足作用力平衡条件与位移协调条件。然而，CSD和CFD对各自计算区域的网格划分要求不同(对于钢储罐内部蒸气云爆炸来说，流体网格比结构网格要密一些)。这就会导致耦合界面上网格的不匹配，在多数情况下节点并不重合，无法直接传递交换数据，因此需要采用合理有效的映射方法进行数据传递[17]。在双向流固耦合分析中，耦合界面流场部分网格上的节点必须映射到耦合界面固体网格的节点上，以传递力参数；耦合界面固体部分网格上的节点必须映射到耦合界面流体网格的节点上，以传递位移参数。因此，流固耦合分析中一次完整的数据交换包含两次在耦合界面处的映射[18]。

耦合界面上不同类型网格之间的数据映射可以看成是数据在不同网格之间的双向插值。ANSYS Workbench中的System Coupling耦合求解模块提供了两种插值方式，即Profile Preserving插值法和Profile Conservative插值法[19]。对于位移等参数的传递，一般采用Profile Preserving插值法；而对于力、动量等参数的传递，一般采用Profile Conservative插值法[20]。

2 钢储罐动力响应的双向流固耦合分析

2.1 流体域模型与固体域模型

本文研究对象为工程上较为常用的5 000 m3拱顶钢储罐，直径20 m、罐壁高17.82 m，具体几何尺寸由中国石油化工集团提供[21]，图1为储罐结构示意图。

图1 储罐结构示意图

图1中：D为钢储罐内径；R0为拱顶半径(R0/D=1.2)；Hw为罐壁高度；Hr为顶盖矢高；H为钢储罐结构总高度，表1为钢储罐罐壁与顶盖厚度。

表1 储罐罐壁与顶盖厚度

钢储罐内部爆炸的双向流固耦合分析模型由两部分组成，即流体域(内部爆炸流场)和固体域(钢储罐结构)。双向耦合分析需要同时在ANSYS Workbench中的Transient Structural模块(ANSYS Mechanical)中创建固体域模型和在Fluid Flow模块(FLUENT)中创建流体域模型。图2为流体域耦合模型，即钢储罐内部密闭空间。为了描述密闭空间内可燃气体爆炸的主要特征并便于求解，对模型做出以下简化假设：① 储罐内部无紧密排列的障碍物(如设备、隔板等)，也无储液(LNG (Liquefied Natural Gas)、LPG (Liquefied Petroleum Gas)等)；② 点火前储罐内部是均匀的可燃气体与空气的混合物，处于常温静止状态；③ 可燃气体、燃烧产物以及空气均满足理想气体状态方程；④ 燃烧气体混合物的比热容随温度变化，满足混合规则；⑤ 可燃气体燃烧反应过程是单步不可逆的；⑥ 爆炸过程中不考虑储罐内部空间与外界的热量传递辐射，容器壁面绝热。起爆前，钢储罐内部充满处于常温静止状态的蒸气云(即均匀混合的乙炔与空气)，乙炔体积分数为10%，罐内初始温度为300 K，初始压力为0.1 MPa。起爆时刻，点火区域处于储罐内部中心位置，点火温度为2 500 K，如图3所示。流体域耦合模型采用k-ε湍流模型、组分概率密度运输燃烧模型、SIMPLE压力/速度耦合算法、亚松弛等数值模型与算法进行模拟。为了减少计算量，模型采用二分之一实体建模并在对称面(B)施加对称边界条件；流体域顶面(C)和壁面(D)分别与储罐顶盖和罐壁接触，是能够随储罐运动变形的流体边界而不再被考虑为刚性壁，故将其设置为流体耦合面；流体域底面(A)与固定的储罐底板接触不能产生变形，设置为固定流体边界。

图4为固体域耦合模型，模型同样采用二分之一壳体建模并在储罐对称面(B)施加对称边界条件；储罐顶盖(C)和罐壁(D)分别与流体域顶面和壁面接触，可在内部爆炸流场的作用下产生变形，故设置为固体耦合面；而储罐底板(A)置为固定边界。

图2 流体域耦合模型

图3 储罐内部初始温度场

图4 固体域耦合模型

2.2 动网格模型

流固耦合分析中流场计算包含十分复杂的湍流流动与化学反应，其计算收敛性与模拟结果的精确性均非常依赖于流场的网格质量。由于组分概率密度运输燃烧模型的计算成本较高，流场网格数量不宜过多，因此需要一种计算效率与精度均可以得到保证的网格划分方式。本文利用ANSYS ICEM CFD网格划分工具，通过创建初始块、关联曲线以及采用O-grid方式将钢储罐内部呈半圆柱形且带拱顶的流场进行O型剖分后再进行网格划分，可以得到网格质量的评价系数均在0.8以上的高质量六面体结构化网格(质量评价系数越接近1，网格质量越高；系数越接近0，网格质量越低)。另外，网格尺寸对流场模拟计算也有一定影响。本文作者通过比较不同网格间距(0.1～0.5 m)模型的计算结果选择了计算效率与精度均能够得到保证的网格间距，即网格单元最大尺寸为0.3 m，网格单元总数为155 356。

另外，当固体耦合面的位移参数传递到流体耦合面后，流体网格需要进行相应的调整。软件FLUENT提供的动网格模型能够有效处理流场边界运动变形引起的计算域改变的流动问题。动网格算法主要包括铺层法、弹性光顺法以及局部重构法。本文采用适用于四边形、六面体或三棱柱网格的铺层法，该算法可以根据流体域的扩张或收缩来相应地生成(拉伸)网格或合并(压缩)网格。储罐结构部分采用四边形网格划分，单元网格尺寸为0.4 m，单元网格数为6 003。

2.3 材料模型与参数

相较于静力作用，爆炸冲击荷载作用下的储罐结构动力响应同时存在着塑性应力硬化与高应变率强化效应。随着应变率的提高，钢材的屈服极限和强度极限均有一定程度的提高而延伸率降低[22]，因而分析结构在爆炸荷载下的弹塑性动力响应时，高应变率效应不可忽略。Cowper-Symonds模型通过在双线性弹塑性模式基础上加入幂指数形式的应变率因子来缩放屈服应力[23]，即同时考虑了塑性应力硬化和应变率强化效应，材料本构方程为

(9)

表2 C-S材料模型参数

值得注意的是，材料本构与参数可以在ANSYS Workbench中的Engineering Data模块进行设置定义，如图5所示。其中，A为静屈服压力σ0，B为塑性硬化模量Ep，n为塑性硬化指数，D为应变率参数常量C，q为应变率参数常量p。

2.4 双向流固耦合的求解流程与解算设置

双向流固耦合需将每一时间步的流场计算结果加载到结构引起变形后，又反过来重新计算流场，如此反复迭代下去，如图6所示。在ANSYS Workbench提供的双向耦合分析模块System Coupling中可以设置流固耦合面的数据交换(力、位移)、计算步长、流场与固体的解算顺序、求解总时长等，如图7所示。由于每一时间步的迭代均需要流体域、固体域以及数据传递的计算收敛，因此双向耦合分析过程相当复杂，计算量十分庞大，对计算设备性能要求较高，计算耗时长，模拟时需要综合考虑收敛性与计算成本的问题。

图5 ANSYS Workbench的Engineering Data模块中的钢储罐材料参数

图6 双向流固耦合求解流程

(a) 流体域力传递给固体域

(b) 固体域位移传递给流体域

2.5 双向流固耦合分析结果

图8为起爆后钢储罐内部压力场的变化情况。可以看出，在爆炸冲击波未到达储罐结构时，压力前波阵面以爆源为中心呈褶皱的球形向外扩张，流体耦合面没有运动变形；在t=35 ms左右，冲击波到达储罐并发生碰撞反射，此时罐内乙炔/空气混合物基本燃烧完全，在冲击波作用下顶盖与罐壁开始向外膨胀变形，由于耦合效应，爆炸流场边界也随之向外膨胀；在t=45 ms左右，由于储罐结构塑性应变到达失效应变而耦合计算终止。在整个过程中，动网格根据流场边界的扩张来相应地拉伸(或生成)网格，以适应流体计算区域的变化，确保流固交界面始终耦合且能够有效传递交换数据。

(a) t=10 ms

(b) t=30 ms

(c) t=35 ms

(d) t=40 ms

(e) t=43 ms

(f) t=45 ms

图9与图10分别给出了储罐结构在不同时刻下的变形云图与等效应变云图。可以看出当爆炸冲击波还未到达储罐时(35 ms之前)，结构会受到邻近流体单元压力上升而出现数量级较小的变形与应变；当爆炸冲击波到达储罐时，罐壁中间区域首先出现相对较大的变形与应变；随着爆炸冲击波的继续扩张和反射，变形与塑性应变区域由罐壁中央逐渐向储罐顶盖与底板发展；在43 ms时，储罐结构最大变形出现在顶盖中心区域 (约0.7 m)，而最大等效应变出现在顶盖与罐壁的环向连接处(约0.206，十分接近失效应变0.21)；爆炸流场继续通过流固耦合作用于储罐，在45 ms时，结构因在顶盖与罐壁环向连接处的最大等效应变已超过失效应变而发生破坏，此时最大变形仍位于顶盖中心区域 (约1.03 m)；底板因与地面固接且建模时采用刚性体模型，因此在爆炸流场作用过程中无变形和应变。

(a) t=10 ms

(b) t=30 ms

(c) t=35 ms

(d) t=40 ms

(e) t=43 ms

(f) t=45 ms

(a) t=10 ms

(b) t=30 ms

(c) t=35 ms

(d) t=40 ms

(e) t=43 ms

(f) t=45 ms

图11、图12分别给出了不同时刻储罐顶盖与罐壁两个区域中最大位移时程曲线与最大变形速度时程曲线，可以看出顶盖结构相比于罐壁结构在爆炸冲击波的作用下产生了较大的动力响应；顶盖区域最大变形(1.03 m左右)约为罐壁区域最大变形(0.56 m左右)1.8倍；顶盖区域最大变形速度(230 m/s左右)约为罐壁区域最大变形速度(140 m/s左右)1.6倍。这是因为，一方面实际工程中储罐罐壁厚度要大于顶盖；另一方面较平坦的顶盖形式对于爆炸冲击荷载更为敏感，相比之下呈圆柱形的罐壁结构有着更好的抗爆性能。

图11 储罐顶盖与罐壁最大位移时程曲线

图12 储罐顶盖与罐壁最大变形速度时程曲线

3 钢储罐动力响应的解耦分析

3.1 爆炸荷载的分区加载

对于内部爆炸作用下储罐动力响应的研究，按是否考虑爆炸流场和储罐结构间的流固耦合效应，可分为解耦方法和耦合方法。双向流固耦合分析相较于解耦分析更符合储罐内爆的实际情况且所得计算结果更为精确，但其建模与计算过程也有如下一些不足：① 双向耦合分析模型同时包含流体域模型、固体域模型以及耦合面数据传递模块，建模过程十分繁琐，对模型的调试也非常耗时；② 模型同时包含内部爆炸流场网格与储罐结构网格，其网格数量庞大，对计算机硬件有较高的要求且计算耗时长；③ 对于每一时间步的计算，固体域、流体域以及耦合面的数据传递都需要分别迭代计算收敛才能推进到下一时间步，而内部爆炸流场涉及到复杂的化学反应并伴随湍流流动，储罐结构在爆炸冲击波作用下的动力响应也涉及到几何与材料的非线性，因此同时收敛相对困难，计算成本非常高；④ 由于耦合界面的存在，储罐结构的大变形或者失效破坏会引起流体耦合面的畸变，形成流体单元的负体积而导致数据传递失败、模拟计算终止。此外，也有一些学者的研究结果表明壳壁变形对内部爆炸流场影响小于1%，多数情况下解耦方法可得到近似程度较高的爆炸流场。作者在文献[21]中将拱顶钢储罐结构考虑为刚性壁，利用软件FLUENT建立了能够模拟内部蒸气云爆炸过程的CFD模型(所采用的湍流燃烧模型与数值方法与本文一致)，获得了罐壁与顶盖上任意位置处完整的爆炸荷载时程曲线。因此，本节将基于不考虑流固耦合效应的内部爆炸荷载对钢储罐结构的动力响应进行解耦分析。

由于中心起爆的爆炸荷载在储罐内表面对称分布，可以将储罐顶盖与罐壁区域分别呈圆环状划分，图13给出了四分之一储罐结构内表面荷载分区情况，共18个荷载区域。在软件FLUENT的计算过程中选择每个区域内中点处布置压力观测点(共18个)，用于提取对应位置的爆炸荷载时程曲线，然后均匀施加于该区域表面。图14给出了其中四条典型的荷载分区曲线。在确定了各个分区的爆炸荷载时程曲线后，需要将其加载于钢储罐有限元模型。利用有限元前后处理软件LS-PREPOST分别创建18个壳单元集合(关键字*SET_SHELL_LIST)与18条荷载分区曲线(关键字*DEFINE_CURVE)，并且一一对应加载(关键字*LOAD_SHELL_SET)。图15给出了在储罐顶盖与罐壁上分区加载18条内部爆炸荷载时程曲线的过程，图中选择加载的区域为黑色。

图13 储罐结构内表面荷载分区

(a) 荷载分区曲线R1

(b) 荷载分区曲线R8

(c) 荷载分区曲线S5

(d) 荷载分区曲线S10

(a) 加载R1～R4

(b) 加载R1～R6

(c) 加载R1～R8

(a) 加载S1～S2

(b) 加载S1～S6

(c) 加载S1～S10

3.2 钢储罐动力响应的解耦分析结果

在有限元前后处理软件LS-PREPOST中进行了内部爆炸荷载的分区加载后，保存关键字文件并交付LS-DYNA求解器进行解耦计算。图16和图17分别给出了解耦分析时储罐的变形云图与等效塑性应变云图。与图9和图10对比可知，采用解耦分析所获得的储罐变形情况和塑性区域发展过程与耦合分析结果基本一致；随着冲击波的作用，变形与塑性应变区域首先出现在罐壁中央并逐渐向顶盖与底板发展；在43 ms时，结构最大变形在顶盖中心区域(约0.9 m)，而最大等效应变出现在顶盖与罐壁的环向连接处(约0.208，接近于失效应变)；在45 ms时，顶盖与罐壁的环向连接处的最大等效应变已超过失效应变，储罐发生破坏，此时最大变形仍位于顶盖中心区域(约1.34 m)。

由于解耦分析是将分区荷载曲线施加在储罐单元上，结构响应的计算并不会像双向耦合分析一样在储罐发生破坏时终止，而是删除相应的失效单元后继续求解计算直到所设置的分析终结时间(即荷载时程曲线的终止时间90 ms)。但是所施加的分区荷载时程曲线后半段(45 ms以后)与实际情况中的爆炸荷载时程曲线会有较大差别，这是因为一方面储罐破坏后罐内的高压力会因破坏泄漏而逐渐下降，另一方面储罐从初始状态发展到结构破坏的过程中已出现了较大的变形(耦合效应作用下会引起爆炸荷载分布的改变)。因此，继续求解计算所得的结构响应不能很好地描述储罐破坏后的实际情况，故本节只给出了解耦分析中储罐从初始状态到出现失效破坏过程的计算结果。

(a) t=40 ms

(b) t=43 ms

(c) t=45 ms

(a) t=40 ms

(b) t=43 ms

(c) t=45 ms

3.3 解耦与耦合结果的对比分析

为了更加清楚地了解两种分析方法所得结果的差异，选取了计算过程中储罐结构的最大位移和最大等效塑性应变时程曲线进行对比分析，如图18所示。从图中可以看出，两种分析方法所得到的最大变形和最大等效塑性应变时程曲线变化趋势基本相同；当爆炸冲击波到达并作用于储罐结构后，解耦分析所得的储罐最大位移与最大等效塑性应变均超出双向耦合分析结果。这是因为，解耦分析中没有考虑爆炸冲击波与储罐结构的耦合效应，而储罐结构的变形会降低冲击波的反射超压峰值，故解耦计算中加载于储罐内表面的爆炸荷载比实际情况的大，从而也导致该分析结果中储罐结构的变形与塑性应变偏大。

4 结 论

本文利用ANSYS Workbench仿真平台分别建立了基于软件FLUENT计算的CFD模型和基于ANSYS计算的固体结构力学模型，实现了内部蒸气体爆炸作用下钢储罐结构动力响应过程的双向流固耦合分析，同时进行了解耦计算，并对两种方法进行了对比分析，得到以下主要结论：

(1) 利用ANSYS Workbench仿真平台中的System Coupling模块能够有效地耦合流固交界面并在每一时间步准确传递交换数据(位移、力)。

(2) 在内部爆炸流场作用下，储罐顶盖与罐壁均会产生向外膨胀变形并在应力较大位置出现塑性区；顶盖结构在爆炸流场作用下会产生更大的动力响应，相较之下呈圆柱形的罐壁结构有着更好的抗爆性能。

(a) 最大变形

( b) 最大等效塑性应变

(3) 储罐结构最大塑性应变一般出现在顶盖与罐壁的环向连接处，该区域会因塑性应变最先超过失效应变而发生破坏。

(4) 采用解耦分析所获得的储罐结构变形情况和塑性区域发展过程与双向流固耦合分析所得结果类似，最大变形与失效破坏区域也与耦合分析结果一致；由于未考虑流固耦合效应而被高估的爆炸荷载以及分区加载带来的误差，解耦分析求得的储罐结构变形位与塑性应变等动力响应参数偏大。

(5) 解耦分析过程相对简单，收敛难度较小，计算成本较低，能够较为方便地进行内部爆炸作用下的钢储罐结构动力响应与薄弱危险区域初步预测，在抗爆设计计算时也偏于安全保守；双向耦合分析模型包含流体域、固体域模型以及耦合面数据传递模块，建模过程繁琐、模型调试费时、网格数量庞大、迭代收敛缓慢，计算成本较高，但相比于解耦分析更符合钢储罐内部蒸气云爆炸的实际情况且所得计算结果更为精确。