仿真笔记——详解LS-DYNA爆炸仿真计算的模型与算法
炸药的起爆和爆炸过程是一种快速的化学反应过程，对于该过程的描述，主要存在 CJ 模型和 ZND 模型两种。LS-DYNA 中包含两种炸药爆轰模型：高能燃烧模型和点火生成模型，前者属于 CJ 模型，后者属于 ZND 模型。点火生成模型需要输入炸药反应率方程参数和未反应炸药的 JWL 方程参数，对于大多数炸药来说这些参数缺乏足够的试验数据支撑，而工程中使用 CJ 模型就足够了，因此下面主要介绍高能燃烧模型。
高能燃烧模型根据炸药上各点距起爆点的距离和炸药爆速来确定每点的起爆时间，如某个炸药单元中心离起爆点位置的距离为ri ，炸药爆速为D ，则该单元起爆时间ti=ri/D ，如果存在多个起爆点，各单元起爆时间按照最近起爆点距离计算。该模型定义爆炸产物压力 P：
其中s P 是依据产物状态方程计算得到的压力， F 为燃烧系数：
其中 Lmin为单元最小特征尺寸。当上式计算的F值大于1时，将F值设置为1，使0≤F≤1 。从式(0.2)可以看出，爆轰波到达前F为零，炸药单元压力为零；当爆轰波经过很短时间后，F迅速增大为1，爆炸产物压力几乎瞬间从零增大到 Ps 。
流体运动方程的描述，按照所采用的坐标系可以分为拉格朗日方法和欧拉方法两大类。拉格朗日方法在物质域内求解流体运动方程，坐标系固定在物质上并跟随物质一起运动和变形，因此也被称为物质描述；欧拉方法在空间域内求解流体运动方程，坐标系固定在空间不动，物质在计算网格之间流动，因此也称为空间描述。
LS-DYNA 中采用任意拉格朗日欧拉方法（ALE）来描述流体运动。该方法在拉格朗日坐标和欧拉坐标之外引入一个可以任意运动的参考坐标系，计算域基于参考域，可以在空间中以任意形式运动。采用ALE 算法的网格同时具有欧拉网格和拉格朗日网格的优点，网格可以随物质一起运动，也可以固定在空间不动，甚至可以在一个方向上随物质运动，而另一个方向上固定不动。
任意拉格朗日欧拉方法描述下的物质导数为：
其中f为物理量，
为质点X 的速度，
为参考点ξ 的速度，也即计算网格运动速度。当
时，计算网格在空间固定不动，退化为欧拉描述；当
时，计算网格随同物体一起运动，退化为拉格朗日描述；当
时，计算网格在空间中独立运动，对应于一般的ALE 描述。
由于爆炸产物和空气的粘性很小，而且爆炸流场运动被视为绝热等熵运动，一般都采用无粘性可压缩流体运动方程来描述爆炸流场。通过式(0.3)将欧拉方法描述的无粘性可压缩流体力学方程变换得到ALE 方法描述的控制方程：
式(0.4)－(0.6)结合空气和爆炸产物的状态方程可以构成封闭的控制方程组。
ALE 方法引入了运动网格，通过在移动边界法向上采用拉格朗日方法，可以很简单地描述边界运动，解决了欧拉方法中移动边界描述困难的问题，给计算带来了很大的方便，但计算过程中需要确定网格的位置。
LS-DYNA 程序中提供了简单 平均算法、体积加权算法、等参算法、等势算法以及混合算法等用于ALE 运动网格位置的确定。但由于爆炸流场计算过程中，爆炸产物和空气界面存在很大的压力和密度梯度，采用以上任何一种算法都会产生异常小的界面网格，从而导致计算无法 正常进行。因此爆炸流场计算中一般仅在边界上采用物质描述，使边界节点速度与界面法向运动速度相等，对于除边界节点外的网格要关闭程序中的网格运动算法，使内部网格退化为空间描述。
当需要考虑壳体影响时，壳体和流场边界可通过共用节点联结，壳体为爆炸流场提供运动边界条件，爆炸流场为壳体施加压力载荷条件，在每一个时间步分步求解即可实现爆炸流场和壳体结构的流固耦合；而当采用刚性壁面假设之后，ALE 方法进一步退化为纯粹的欧拉方法。
炸药爆炸后，爆炸容器内存在爆炸气体产物和空气两种物质，这两种气体的流场都可以通过上述无粘性可压缩方程描述，但计算过程中需要区分两种不同介质，并捕捉两种物质的界面。LS-DYNA 中通过定义物质组号来区分不同介质，采用杨氏流体体积法（Yong’s VOF）来捕捉两种物质的运动界面，具体过程是：采用多物质单元划分爆炸流场计算域，一个单元中允许同时存在多种物质；先假设界面沿单元边界，根据与节点相邻的所有单元中存在的物质计算各个节点上的物质体积分数；由同一单元网格各节点的物质体积分数梯度确定界面的法向，并构造该单元内界面；计算每个时间步内通过四周流到相邻单元的流体体积，修改网格单元和相邻单元中的流体体积分数；由各个边界单元内界面组成整个物质界面。
欧拉和ALE 描述的运动方程求解一般有两种方式：全耦合求解和算子分裂法，前者是在整个计算域同时求解运动方程，后者将每个时间步分为拉格朗日阶段和对流阶段依次计算。全耦合求解一个计算单元只能存在一种物质，不适合求解多物质问题。
LS-DYNA 程序中采用算子分裂法求解。拉格朗日阶段采用有限元方法计算由于外力和内部应力产生的速度、压力和内能变化以及现时密度，单元采用单点积分并通过沙漏粘性控制零能模式，引入人工粘性以捕捉冲击波，时间推进采用二阶精度的中心差分法；求得这些参数后再进行界面捕捉，构造多物质内界面。对流阶段采用有限体积法计算通过单元边界的质量、动量和能量通量，通量的计算可以采用一阶精度的迎风算法或者二阶精度的Van Leer 对流算法；该阶段时间步不发生变化，保持与拉格朗日阶段一致。爆炸流场计算一般采用Van Leer对流算法，因为这种算法不仅具有二阶精度，而且具备总变差递减（TVD）性质。
其中P 为空气压力；γ 为多方指数；ρ 为空气现时密度，ρ0为初始密度；E为空气能量密度。理想气体状态方程可以通过设置LS-DYNA 程序中预定义的线性多项式状态方程的相关常数得到。
目前存在的爆炸产物状态方程，按照是否显含产物组分可以分为两类，一类是显含产物组分的状态方程，如BKW 方程、KHT 方程等，另一类是不显含产物组分的状态方程，如多方方程、JWL 方程等，其中JWL 方程是数值计算中应用最为广泛的状态方程。
JWL 方程是1968 年由美国LLNL 的E. L. Lee 在H. Jones 和M. L. Wilkins 的工作基础上提出的半经验状态方程。JWL 等熵方程形式如下：
其中
为相对比容，即现时比容
与初始比容
之比，A 、B 、C 、1 R 、2 R 、ω 为六个参数。
根据热力学第一定律和等熵条件，比内能增量
，将(0.8)式代入并积分可得比内能
。
根据式(0.8)和(0.9)，可消去参数C ，得到
其中
为以初始体积表示的能量密度，初始能量密度
。式 (0.10)为 LS-DYNA、AUTODYN、CTH、MESA 等爆炸动力学计算软件中所普遍采用 JWL 方程形式。该形式只要求输入
五个参数值和初始能量密度
，结合炸药高能燃烧模型要求输入的炸药密度
、爆速 D 和 CJ 爆压
，可以得到爆炸时刻的
等，然后通过增加相对体积微量
，根据式 (0.10) 和原先的
值，计算得到新的
值，新的
值调整为
，如此反复一直计算。
炸药 JWL 方程参数的确定需要通过圆筒试验和二维流体弹塑性数值计算相结合的方法确定：先进行圆筒试验，将待测炸药装入紫铜管，一端起爆，用高速摄像仪记录下铜管外径的运动轨迹；设定一组
值，根据式(2.11)、(2.12)和(2.13)求得A 、B 、C ，利用假设的JWL 方程通过二维流体弹塑性程序数值模拟炸药驱动圆管的外径膨胀轨迹，如果数值计算结果与和试验结果相对误差小于1%，则假设参数即为真实JWL 方程参数，如果不满足相对误差要求，则继续调整系数，直到和试验结果相对误差小于1%为止。
式中
为炸药等容爆热，
为炸药CJ 状态时的比容，
