ANSYS/LS-DYNA 如何模拟液体的晃动过程

液体晃动（sloshing）问题的研究在实际工程中有重要的意义，比如在石化工业中广泛应用的大型储罐，一般直径在几十米，甚至上百米。在地震或其他意外条件下液面的波动情况如何，是否存在安全隐患，都需要进行数值模拟研究。下面就 ANSYS/LS-DYNA 软件在这方面的应用。

众所周知，ANSYS/LS-DYNA 在显式计算领域占据主导地位，随着各种新的算法的不断采用，在求解的广度、精度以及效率上，ANSYS/LS-DYNA 具有同类软件所无法比拟的优势。针对液面晃动问题，ANSYS/LS-DYNA 提供以下三种方法：

1. 流固耦合

流固耦合是 ANSYS/LS-DYNA 计算流体和结构间相互作用的最常用的方法，包括单物质 + 空材料和多物质耦合两大类，流体单元有 Euler 和 ALE 两种。其涉及的主要命令如下：

• *control_ale

算法选择有两种 2、3，分别为 Euler 和 ale 实质上此处二者没有区别，只是因为兼容性进行的设置；两种精度供选择 - 单精度、双精度。

• *section_solid_ale

对单物质 + 空材料为 12 号算法，对多物质耦合为 11 号算法。

• *ale_multi-material_group

进行多物质的定义，最多可以定义 20 种材料。可以根据物质间能否混合将各种材料定义在不同的材料组 ID 中。

• * ale_multi-material_system_group

该命令决定流体物质的算法（Euler 或 Ale），或是在运算过程中切换使用两种算法，并可对流体物质进行自由度约束。该命令多与下列三个命令结合使用：

• * ale_multi-material_system_curve

定义 ale 系统的运动曲线。

• * ale_multi-material_system_node

通过一系列节点定义 ale 的运动参考坐标系统。

• * ale_multi-material_system_switch

定义 Euler 和 ale 参考系统的切换。

上述命令是流体物质涉及的关键字，而我们知道，结构采用 Lagrange 单元来离散，二者之间的耦合通过下列命令来实现：

• *constrained_Lagrange_in_solid

耦合算法分为两种： 罚耦合和运动约束 。前者遵循能量守恒，后者遵循动量守恒。一般令结构网格较流体网格密以保证界面不出现渗透，否则可以增大 NQUAD 参数值来增加耦合点，如设置该值为 4 或 5。在 970 中，此命令第三行又增加了一个控制字 ILEAK-0，1 或 2，一般可设置为 1。

最后给出一个典型算例 - 水箱跌落的部分关键字：

*keyword
*TITLE
boxwater2.k: dropping a water box onto a rigid platform
$========================================================================
$ [1] EXECUTION CONTROLS
$========================================================================
*CONTROL_TERMINATION
$ ENDTIM ENDCYC DTMIN ENDENG ENDMAS
0.0500000 0 0.0000000 0 0.0000000
*CONTROL_TIMESTEP
$ DTINIT TSSFAC ISDO TSLIMT DT2MS LCTM ERODE MS1ST
0.0000000 0.2000000 0 0.0000000 0.0000000 0 0 0
*CONTROL_ENERGY
$ HGEN RWEN SLNTEN RYLEN
2 2 2 2
$========================================================================
$ [3] OUTPUT CONTROLS
$========================================================================
*DATABASE_BINARY_d3plot
$ DT CYCL LCDT Beam
0.0005000 0
*DATABASE_GLSTAT
0.0001000
$========================================================================
$ [5] |SECTIONS|PARTS| DEFs
$========================================================================
*PART
water in the box
$ PID SECID MID EOSID HGID GRAV ADPOPT TMID
1 1 1 1 0 0 0 0
*SECTION_SOLID_ALE
$ SECID ELFORM AET
1 12
$ AFAC BFAC CFAC DFAC START END AAFAC
0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
$-------------------------------------------------------------------------------
*MAT_NULL
$ MID RHO PC MU TEROD CEROD YM PR
1 1000.0000 -1.000+10 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
*EOS_LINEAR_POLYNOMIAL
$ EOSID C0 C1 C2 C3 C4 C5 C6
1 0.0000000 1.50000+9 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
$ E0 V0
0.0000000 1.0000000
$========================================================================
*PART
void portion in the box
2 1 1 1 0 0 0 0
*INITIAL_VOID_PART
2
$========================================================================
*PART
rigid box containing water
$ PID SECID MID EOSID HGID GRAV ADPOPT TMID
3 3 3 0 0 0 0 0
*SECTION_SOLID
$ SECID ELFORM AET
3 0
*MAT_RIGID
3 2000.0000 1.00000+8 0.0000000 0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
$========================================================================
*PART
rigid super-heavy platform
$ PID SECID MID EOSID HGID GRAV ADPOPT TMID
4 4 4
*SECTION_SHELL
$ SID ELFORM SHRF NIP PROPT QR/IRID ICOMP
4 0
$ T1 T2 T3 T4 NLOC
0.011 0.011 0.011 0.011
*MAT_ELASTIC
$ MID RHO E PR DA DB K
4 1000000.0 1.0000+14
$========================================================================
$ [8] BC's + IC's + BODY LOADS + FORCE FIELDS
$========================================================================
*INITIAL_VELOCITY
$ NSID NSIDEX BOXID
0
$ VX VY VZ VXR VYR VZR
0.0 -20.0 0.0
$-------------------------------------------------------------------------------
*LOAD_BODY_Y
$ LCID SF LCIDDR XC YC ZC
1 1.00
*DEFINE_CURVE
$ LCID SIDR SFO OFFA OFFO DATTYP
1
$ X=abcissa Y=ordinate
0.0 981.0
1.0 981.0
$========================================================================
$ [9] LAGRANGIAN CONTACTS CONSTRAINTS, ……
$========================================================================
$ SFS = scale fact on dflt SLAVE penal stifns (see CONTROLL_CONTACT)
$ SFM = scale fact on dflt MASTER penal stifns (see CONTROLL_CONTACT)
*CONTACT_AUTOMATIC_NODES_TO_SURFACE
$ SSID MSID SSTYP MSTYP SBOXID MBOXID SPR MPR
3 4 3 3
$ FS FD DC VC VDC PENCHK BT DT
$ SFS SFM SST MST SFST SFMT FSF VSF
100. 100.
$========================================================================
$ [10] EULERIAN & ALE CONTACTS CONSTRAINTS, ……
$========================================================================
*CONTROL_ALE
$ DCT NADV METH AFAC BFAC CFAC DFAC EFAC
2 1 4-1.0000000 0.0000000 0.0000000 0.0000000
$ START END AAFAC VFACT VLIMIT EBC
0.0000000 0.0000000 0.0000000 0.0
*ALE_REFERENCE_SYSTEM_GROUP
$ SID STYPE PRTYP PRID BCTRAN BCEXP BCROT ICOORD
1 0 5 1
$ XC YC ZC EXPLIM
*SET_PART_LIST
$ SID DA1 DA2 DA3 DA4
1
$ PID1 PID2 PID3 PID4 PID5 PID6 PID7 PID8
1 2
*ALE_REFERENCE_SYSTEM_NODE
$ NSID
1
$ N1 N2 N3 N4 N5 N6 N7 N8
5 6 7

2. SPH 算法

SPH 算法作为 DYNA 中第一种无网格（meshfree）算法，在连续体的破碎或分离分析中得到了广泛的关注和应用。在解决极度变形和破坏类型的问题上 SPH 有着其他方法无法比拟的优势，可以说无网格算法正在成为数值分析领域的研究热点，具有很好的发展前景。

我们知道传统的有限单元法中，单元的形状对结果的精度影响很大，如果单元因为变形过大可能造成矩阵奇异，使得精度降低甚至无法计算下去。而 SPH 算法则是把每个粒子作为一个物质的插值点，各个粒子间通过规则的内插函数计算全部质点即可得到整个问题的解。