冲击载荷下柔性储液罐动态响应数值模拟及规律分析＊

曹 源，金先龙，李 政

（上海交通大学机械系统与振动国家重点实验室，上海 200240）

储液容器广泛应用于储运、航空等领域，随着材料性能的提高，在许多抗冲击场合，如空投、爆炸中，柔性储液容器开始被越来越多的应用。许多学者已对容器内液面晃动、容器冲击响应以及流固耦合等问题进行了长期深入的研究［1－3］。与早期研究相比，柔性容器冲击响应问题的难点是大变形和非线性，使得其在冲击过程中的动态响应规律更加复杂。对于这类问题的研究，目前还是以实验为主。随着计算方法的改进和计算机性能的提高，数值模拟手段已有可能发挥更大的作用。

对于大变形的流固耦合问题，很难用单一的Lagrange方法或Euler方法进行研究。所以ALE方法开始应用于流固耦合分析［4－5］。ZHANG Ai－nian等［6］分别使用附加质量法、Lagrange方法和 ALE方法分析了船舶碰撞时流固耦合作用，并认为与ALE方法相比，附加质量方法和Lagrange方法低估了液体与结构之间相互作用的影响。M.Anghileri等［7］研究了储液箱体与地面碰撞时，箱内液面大幅晃动及箱体变形的问题，并将ALE方法与其他数值模拟方法进行了比较，并且与实验结果吻合很好。

本文中针对跌落或空投中的柔性储液罐动态响应进行实验研究和数值模拟。采用基于Lagrange方法描述下的Belytschko－Tsay（B.T）膜单元建立柔性储液容器模型，流体部分则是基于ALE方法建模，并利用罚函数法实现流固耦合。在实验基础上，使用数值模拟的方法深入分析了空投高度、装水量等因素对容器动态响应的影响规律。

1 流固耦合方法与原理

利用ALE方法建模与罚函数法相结合对冲击过程中的流固耦合现象进行数值模拟。罚函数法的优点在于保证了耦合接触过程中的能量守恒。首先，对于每个结构节点，搜索包含该节点的ALE单元。然后，通过计算穿透速度与时间来计算流体节点的穿透深度。最后，通过正比穿透深度计算基于罚函数的耦合力。耦合力可以看作是节点总力中的一个外部力。在此基础上，通过计算每一个时间步长的总节点力，可以得到流体和结构体耦合界面上每一点的速度、位移和应力等物理量，从而实现流固耦合界面上的相互作用，即耦合计算。ALE方法则综合了Lagrange方法与Euler方法的优点，在材料域与空间域外引入了参考域，并在参考域网格上求解，既解决了Lagrange方法中材料严重扭曲的问题，又解决了Euler方法中移动边界引起的复杂性问题。

基于ALE方法的连续性方程和动量方程为

式中：χ为ALE坐标；vi为水的流动速度；xi和xj为空间坐标；ci和cj均为对流速度；ρ为流体密度；bi为流体体力，σij为应力张量。

对储液容器变形的描述使用弹性体的连续方程

式中：X是Lagrange坐标，ρs为容器密度，fi为体力，u为固体结构位移。

对于水这种弱可压缩流体，引入线性牛顿流体本构方程

式中：p为水的静态压力，μ＝μ（sij）为动力粘性系数，sij为应变率张量，δij为克罗内克常数。

流固耦合需在液体和结构界面上满足几何相容条件和力的平衡条件，即

式中：vi和vs分别为流固耦合界面上水的流动速度和固体材料速度，Ff和Fs分别是流体和固体结构作用在流固耦合界面上的力，这些力可由罚函数方法计算获得

2 实验简介

图1 储液罐结构Fig.1 Stucture of container

为验证本文中方法和模型的正确性，首先对装水量η＝80%的柔性储液罐进行100m高度的空投冲击实验。所使用的储液罐为∅0.5m×0.8m的圆柱体，如图1所示。储液罐壁厚为10mm。罐体采用帘布材料，该帘布采用尼龙织物，基体采用天然橡胶。其顶端由钢盖密封。在100m高度，以零初始速度垂直抛落，整个抛落过程由高速照相机记录。

实验中，罐体肩部（位置1）破裂，储液罐内的水溅出。所以，罐体肩部被确定为脆弱位置之一。此外，罐体中部（位置2）和底部（位置3）也是易受冲击和易损部位，这些部位也是数值模拟时关注的重点。

3 数值模拟与分析

3.1 数值模型

以实验为基础，建立三维有限元数值模型。该模型包括：地面、罐状柔性储液容器、水和容器内空气。相对于柔性容器，冲击过程中混凝土地面变形较小，模型中将地面材料作为弹塑性体处理。容器的钢盖处理为弹性体。各项参数如表1所示，其中，ρ为密度，E为弹性模量，ν为泊松比，σy为屈服应力，G为剪切模量，β为硬化系数。

容器主体复合材料结构的力学性能是数值模拟的关键，本文中采用 XIA Zhi－hui等［8］提出的周期性边界加载方法计算获得帘布材料的力学性能参数，如表2所示。其中E1、E2和E3分别为材料3个方向的弹性模量；G12、G23和G13分别为材料3个方向的剪切模量；ν12、ν23和ν13分别为材料3个方向的泊松比。由表2可知材料呈横观各向异性。由于容器壁较薄，所以采用基于Lagrange方法描述的B.T膜单元，沿膜厚度方向取3个积分点。

表1 材料参数Table1 Material porperties and parameters

表2 容器材料力学参数Table2 Mechanical properties of the composite

为了节省计算时间，忽略罐体坠落过程中的空气阻力，将抛落高度换算为罐体冲击地面时的速度，作为初始条件。对于内部流体采用基于ALE方法描述的六面体单元。为了确保结果的准确性，所有单元尺寸和形状均仔细控制。整个模型单元数为228 033，节点数为238 484。

数值模拟采用基于显示计算方法的ANSYS／LS－DYNA有限元软件，为确保计算的稳定性，时间步长由模型的最小单元尺寸控制，为约2.0μs。在流固耦合的模拟中，过多的耦合点将导致数值计算不稳定，而耦合点较少则易产生泄漏现象。由于柔性容器为曲面结构，存在较小的夹角，故耦合点数取为4。此外，耦合最小体积参数决定了流体物质达到某个单元体积时，流体与结构发生耦合作用。由于空投物体在触地时属于高速碰撞冲击问题，极易产生泄漏现象，所以计算中将耦合最小体积参数值设为0.1，可以有效地减小流体的泄漏现象。

3.2 数值模拟与实验结果比较

对相同条件下的数值模拟和实验结果进行比较，以验证方法和模型的有效性。定义罐体直径变化量Δd来评估柔性罐体的变形情况。所谓直径变化量，是指冲击过程中不同时刻罐体中间直径d与盛有相应水量的罐体静置于地面时罐体中间直径d0之差

图2 变形量比较Fig.2 The comparison of the deformation

对罐体径向变形量的实验结果和柔性罐状容器的数值模拟结果相比较发现，数值模拟结果相对于物理实验结果要小一些，但仍处于可接受的范围，如图2所示。

不同时刻下柔性容器应力分布和内部流体的变化，如图3所示。从图中可以看出，流体形态与柔性容器变形一致，其内部的空气先是被压缩，随着柔性容器的反弹，空气完全被水包围。当冲击开始时，柔性容器表面最大应力出现在最先与地面接触的底部；反弹时，由于钢盖惯性较大，在顶部形成“弹坑”形状，此时，容器肩部变形最大，应力也最大。

图3 柔性容器表面应力分布及流体形态Fig.3 Stress distribution of the container and configuration of the fluid

柔性容器上3个位置主应力σp变化曲线，如图4所示。从曲线中可以看出，3个位置的应力依次达到最大值，与这些位置受冲击的先后顺序一致。其中，肩部的最大应力几乎是其他2个位置的2倍。这也与实验中罐体反弹过程中，肩部破裂的结果十分吻合。

3.3 响应规律分析及讨论

从抛落高度和装水量等2个方面对容器的冲击响应进行研究。对装水量η＝80%时，100、200和400m等3个抛落高度进行数值模拟。图5显示了这3种高度下，3个位置上最大主应力σm的比较。

图4 关键位置处主应力曲线Fig.4 Principal stresses at key positions

图5 不同高度下最大应力Fig.5 Principal stresses at different heights

柔性容器肩部的应力一直大于另外2个位置的，而随着高度的增加，容器中部的应力增加最快。当高度达到400m时，该位置的最大应力几乎与肩部相同。所以，随着高度的增加，容器中部也成为危险部位。另一方面，容器底部的应力随高度的增加，升高缓慢。特别是当高度低于200m时，容器底部应力变化不大。

图6给出了3个抛落高度下柔性容器直径的变化曲线。容器直径变化量和最大幅度都随高度的增加而变大。其中容器直径变化量的最大值从30.6mm增加到64mm。图中还可以看出，尽管抛落高度不同，但振荡收敛的时间没有大的改变。

图7显示的是空抛落高度100m时，装水量η＝30%，50%，80%和100%等4种情况下容器上3个位置的最大应力。在装水量不同的情况下仍然是容器肩部应力最大，且与装水量呈“V”型关系。当η＝80%时，肩部应力是4种情况下最小的。原因可能是：当装水量少时，钢盖可以上下大幅运动，导致容器肩部产生大变形。随着装水量的增加，肩部变形减小，应力也降低。但是当容器100%充满水时，冲击产生的力和能量可以通过水直接作用在肩部。由于水的可压缩性很小，刚度又大于柔性材料，减震效果小于容器本身，所以导致肩部的应力急剧增加。与肩部应力的变化不同，容器底部和容器中部的应力与装水量的变化的关系不大。还发现容器的变形量随着装水量增加而呈线性增加。

图6 不同高度下柔性容器直径变化Fig.6 Container deformations at different heights

图7 不同装水量下最大应力Fig.7 Principal stresses under different amount of water

4 结 论

在对柔性储液罐进行空投实验的基础上，采用基于ALE方法描述的流固耦合状态和三维有限元模型对柔性储液容器的冲击响应进行研究。实验和数值模拟结果都显示，柔性储液容器的肩部是其薄弱部位。当达到一定的抛落高度时（如400m），容器中部也会变成另外一个薄弱部位。容器的装水量对其抗冲击性有很大影响。其中肩部的应力与装水量呈“V”型关系。而柔性容器的变形则基本与装水量呈线性关系。在其他条件不变的情况下，装水量为80%时，柔性储液容器的抗冲击性最好。

［1］Moiseev N N.Introduction to the theory of oscillations of liquid containing bodies［J］.Advances in Applied Mechanics，1964，8（5）：233－289.

［2］Ma D C，Su T C.Fluid－structure Vibration and Liquid Sloshing［C］∥Proceedings of Pressure Vessels and Piping Conference.New York：American Society of Mechanical Engineers，1987.

［3］Shimizu N.Advances and trends in seismic design of cylindrical liquid storage tanks［J］.Japan Society of Mechanical Engineers International Journal，1990，33（2）：111－124.

［4］Nitikitpaiboon C，Bathe K J.An arbitrary Lagrangian－Eulerian velocity potential formulation for fluid－structure interaction［J］.Computers ＆Structures，1993，47（4／5）：871－891.

［5］Nomura T.ALE finite element computations of fluid－structure interaction problems［J］.Computer Methods in Applied Mechanics and Enginneering，1994，112（1／2／3／4）：291－308.

［6］ZHANG Ai－nian，Suzuki K.A comparative study of numerical simulations for fluid－structure interaction of liquidfilled tank during ship collision［J］.Ocean Engineering，2007，34（5／6）：645－652.

［7］Anghileri M，Luigi M L，Castelletti，Tirelli M.Fluid－structure interaction of water filled tanks during the impact with the ground［J］.International Journal of Impact Engineering，2005，31（3）：235－254.

［8］XIA Zhi－hui，ZHANG Yun－fang，Ellyin F.A unified periodical boundary conditions for representative volume elements of composites and applications［J］.International Journal of Solids and Structure，2003，40（8）：1907－1921.