LS-DYNA爆炸模块计算算法概述

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