LSDYNA爆炸模块计算算法概述

LS-D YNA爆炸模块计算算法概述Den gguide for Simwe 2009-3-18LS-DYNA是世界最著名的显式非线性动力学计算程序，广泛应用于冲击侵彻、高速碰撞、金属成形等固体结构非线性动力学分析。近年来，LS-DYNA中加入了爆炸流场计算功能，可以用于计算爆炸容器的内部爆炸流场和壁面爆炸载 荷。1.1.1爆轰模型炸药的起爆和爆炸过程是一种快速的化学反应过程，对于该过程的描述，主 要存在CJ模型和ZND模型两种。LS-DYNA中包含两种炸药爆轰模型：高能燃 烧模型和点火生成模型，前者属于CJ模型，后者属于ZND模型。点火生成模型需要输入炸药反应率方程参数和未反应炸药的JWL方程参数，对于大多数炸药来说这些参数缺乏足够的试验数据支撑，而工程中使用CJ模型就足够了，因此下面主要介绍高能燃烧模型。高能燃烧模型根据炸药上各点距起爆点的距离和炸药爆速来确定每点的起爆 时间，如某个炸药单元中心离起爆点位置的距离为ri，炸药爆速为D，则该单元起爆时间ti =ri /D，如果存在多个起爆点，各单元起爆时间按照最近起爆点距离 计算。该模型定义爆炸产物压力 P166:P= FPS(V,E) (0.1)其中巳是依据产物状态方程计算得到的压力，F为燃烧系数：?0,t 岂?F = ? t-ti,t ti (0.2)?1.5Lmin / D其中Lmin为单元最小特征尺寸。当上式计算的 F值大于1时，将F值设置为1， 使0 WF 1 0从式(0.2)可以看出，爆轰波到达前F为零，炸药单元压力为零；当 爆轰波经过很短时间后，F迅速增大为1，爆炸产物压力几乎瞬间从零增大到 R。1.1.2运动方程流体运动方程的描述，按照所采用的坐标系可以分为拉格朗日方法和欧拉方 法两大类。拉格朗日方法在物质域内求解流体运动方程，坐标系固定在物质上并跟随物质一起运动和变形，因此也被称为物质描述；欧拉方法在空间域内求解流 体运动方程，坐标系固定在空间不动，物质在计算网格之间流动，因此也称为空 间描述。LS-DYNA中采用任意拉格朗日欧拉方法（ALE ）来描述流体运动167。该方 法在拉格朗日坐标和欧拉坐标之外引入一个可以任意运动的参考坐标系， 计算域 基于参考域，可以在空间中以任意形式运动。采用 ALE算法的网格同时具有欧 拉网格和拉格朗日网格的优点，网格可以随物质一起运动，也可以固定在空间不 动，甚至可以在一个方向上随物质运动，而另一个方向上固定不动。1.1.2.1 ALE 描述任意拉格朗日欧拉方法描述下的物质导数为：dnXd） = ?Hm +（u-Wi）皿（0.3）dt?t?X其中f为物理量，u为质点X的速度，Wi为参考点E的速度，也即计算网格运 动速度。当Wj =0时，计算网格在空间固定不动，退化为欧拉描述；当wi = u时, 计算网格随同物体一起运动，退化为拉格朗日描述；当wi工0时，计算网格在空间中独立运动，对应于一般的 ALE描述。由于爆炸产物和空气的粘性很小，而且爆炸流场运动被视为绝热等熵运动， 一般都采用无粘性可压缩流体运动方程来描述爆炸流场。通过式（0.3）将欧拉方法 描述的无粘性可压缩流体力学方程变换得到 ALE方法描述的控制方程：? + ( u - Wi)?t?p +Xi噥=0 (04)? uj? uj?pp1 + p(u-Wi)L=-(0.5)?t?x?Xj?E?E? uP+ P(u - Wi)=- p - (0.6)?t?xi?n 式(0.4)(0.6)结合空气和爆炸产物的状态方程可以构成封闭的控制方程组。1.1.2.2 网格运动ALE 方法引入了运动网格，通过在移动边界法向上采用拉格朗日方法， 可以 很简单地描述边界运动，解决了欧拉方法中移动边界描述困难的问题， 给计算带 来了很大的方便，但计算过程中需要确定网格的位置。LS-DYNA 程序中提供了简单平均算法、体积加权算法、等参算法、等势算 法以及混合算法等用于 ALE 运动网格位置的确定 71 。但由于爆炸流场计算过程 中，爆炸产物和空气界面存在很大的压力和密度梯度， 采用以上任何一种算法都 会产生异常小的界面网格，从而导致计算无法正常进行。因此爆炸流场计算中一 般仅在边界上采用物质描述，使边界节点速度与界面法向运动速度相等， 对于除 边界节点外的网格要关闭程序中的网格运动算法，使内部网格退化为空间描述。当需要考虑壳体影响时，壳体和流场边界可通过共用节点联结，壳体为爆炸 流场提供运动边界条件，爆炸流场为壳体施加压力载荷条件，在每一个时间步分 步求解即可实现爆炸流场和壳体结构的流固耦合；而当采用刚性壁面假设之后， ALE 方法进一步退化为纯粹的欧拉方法。1.1.2.3 界面捕捉炸药爆炸后，爆炸容器内存在爆炸气体产物和空气两种物质，这两种气体的 流场都可以通过上述无粘性可压缩方程描述， 但计算过程中需要区分两种不同介 质，并捕捉两种物质的界面。 LS-DYNA 中通过定义物质组号来区分不同介质， 采用杨氏流体体积法( Yongs VOF)168 来捕捉两种物质的运动界面，具体过程 是：采用多物质单元划分爆炸流场计算域，一个单元中允许同时存在多种物质； 先假设界面沿单元边界，根据与节点相邻的所有单元中存在的物质计算各个节点 上的物质体积分数；由同一单元网格各节点的物质体积分数梯度确定界面的法 向，并构造该单元内界面；计算每个时间步内通过四周流到相邻单元的流体体积， 修改网格单元和相邻单元中的流体体积分数； 由各个边界单元内界面组成整个物 质界面。1.1.2.4 求解方式欧拉和 ALE 描述的运动方程求解一般有两种方式：全耦合求解和算子分裂 法，前者是在整个计算域同时求解运动方程， 后者将每个时间步分为拉格朗日阶 段和对流阶段依次计算。全耦合求解一个计算单元只能存在一种物质， 不适合求 解多物质问题。LS-DYNA 程序中采用算子分裂法求解。拉格朗日阶段采用有限元方法计算 由于外力和内部应力产生的速度、压力和内能变化以及现时密度，单元采用单点 积分并通过沙漏粘性控制零能模式，引入人工粘性以捕捉冲击波，时间推进采用 二阶精度的中心差分法；求得这些参数后再进行界面捕捉，构造多物质内界面。 对流阶段采用有限体积法计算通过单元边界的质量、 动量和能量通量，通量的计 算可以采用一阶精度的迎风算法或者二阶精度的 Van Leer 对流算法；该阶段时 间步不发生变化，保持与拉格朗日阶段一致。爆炸流场计算一般采用Van Leer对流算法，因为这种算法不仅具有二阶精度，而且具备总变差递减( TVD )性质 71 。o1.1.3 状态方程1.1.3.1 空气空气采用理想气体状态方程描述：P =(Y- 1)pE/ p (0.7)其中P为空气压力；丫为多方指数；P为空气现时密度，po为初始密度；E为空 气能量密度。理想气体状态方程可以通过设置 LS-DYNA 程序中预定义的线性多 项式状态方程的相关常数得到。1.1.3.2 爆炸产物目前存在的爆炸产物状态方程，按照是否显含产物组分可以分为两类，一类 是显含产物组分的状态方程，如 BKW 方程、 KHT 方程等，另一类是不显含产 物组分的状态方程，如多方方程、 JWL 方程等，其中 JWL 方程是数值计算中应 用最为广泛的状态方程 169 。JWL 方程是 1968 年由美国 LLNL 的 E. L. Lee 在 H. Jones和 M. L. Wilkins 的 工作基础上提出的半经验状态方程170 o JWL等熵方程形式如下：PS = Ae- RV + Be- R2V + CV- (+1) (0.8)其中V=v/v= p。/P为相对比容，即现时比容v与初始比容V。之比，A、B、C、R1、R2、3为六个参数根据热力学第一定律和等熵条件，比内能增量de=- pdv,将(0.8)式代入并积分可得比内能eses = Ae-RlV + B e-谢 +CV- ”v0 (0.9) RR23根据式(0.8)和 (0.9)，可消去参数C ,得到PS = A(1- 土)e 敵 +B(1-仝)eR2V + 眷扌RVR2V(0.10)其中Es = pes为以初始体积表示的能量密度，初始能量密度E = pe = pQv。式(0.10)为LS-DYNA、AUTODYN、CTH、MESA等爆炸动力学计算软件中所普遍采用JWL方程形式。该形式只要求输入 A、B、R、R2、3五个参数值和初始能量密度E0，结合炸药高能燃烧模型要求输入的炸药密度p、爆速D和CJ爆压PCj，可以得到爆炸时刻的Ps、Es、Vs等，然后通过增加相对体积微量 Z， 根据式(0.10)和原先的Es值，计算得到新的Ps值，新的Es值调整为 Es = Es - Ps AV，如此反复一直计算171。炸药JWL方程参数的确定需要通过圆筒试验和二维流体弹塑性数值计算相 结合的方法确定：先进行圆筒试验，将待测炸药装入紫铜管，一端起爆，用高速 摄像仪记录下铜管外径的运动轨迹；设定一组R1、R2、3值，根据式(2.11)、(2.12) 和(2.13)求得A、B、C，利用假设的JWL方程通过二维流体弹塑性程序数值模 拟炸药驱动圆管的外径膨胀轨迹，如果数值计算结果与和试验结果相对误差小于 1%,则假设参数即为真实JWL方程参数，如果不满足相对误差要求，则继续调 整系数，直到和试验结果相对误差小于1%为止172。Ae- CJ + Be- R2Vcj + CVCJ-(3+1) = FCj (0.11)ABC吕e RVcj +FeR2VCJ +CVcj-3 = pcQv- FCj(Vcj-1)/2 (0.12) RiR23ARe RVcj +BRe 弘 +C(1+ 3)Vcj-(3+2) = pD2 (0.13)式中Qv为炸药等容爆热，Vcj为炸药CJ状态时的比容,PCJVcj =1-HJr (0.14)PoD