基于OTM方法的超高速碰撞问题数值模拟

马天宝，粟鑫，郝莉

(1.北京理工大学 爆炸科学与技术国家重点实验室，北京 100081; 2.北京建筑大学 理学院，北京 100044)

在航天技术高速发展的今天，在轨航天器正时刻受到外层空间碎片的撞击威胁.这些空间碎片分布在距地球表面2 000 km的低地轨道上，拥有极大的数量和极高的运行速度，与航天器发生碰撞时其撞击速度可达10 km/s以上，因此现今大多数的航天器多数采用在航天器舱壁外加装Whipple屏[1]的防护形式来防护空间碎片.在对空间碎片防护技术的研究中，由于碰撞速度很高，受限于目前实验器材发展水平，实验研究存在很大的困难，因而数值模拟成为了研究航天器碎片防护技术中超高速碰撞问题的重要方法[2-3].

目前，超高速碰撞问题的数值模拟多采用光滑粒子流体动力学方法(smoothed particle hydrodynamics，SPH)及其相关的衍生方法进行模拟.但SPH方法本身也存在一些问题，如存在拉力不稳定性，施加位移边界条件存在困难等.最优输运无网格方法是一种显式增量更新拉格朗日无网格计算方法，由M.Ortiz等[4]于2010年提出.该方法的理论体系由最优输运理论、物质点采样法以及局部最大熵无网格近似[5]结合而成，与SPH方法相比，OTM方法具有以下优势：① 最优输运理论为OTM方法中的显式时间积分算法提供了可靠的理论基础，保证了时间离散下系统的动量守恒与辛守恒，同时质量守恒自动满足无需额外计算，可以有效地克服现有显式分析中能量不守恒的问题；② OTM在计算边界处具有Kronecker-delta性质，可以更好地引入位移边界条件；③ OTM方法采用物质点与节点结合的方式进行空间离散，由于插值与积分在不同的离散点进行，完全避免了SPH方法中存在的拉力不稳定性.本文采用OTM方法模拟弹丸超高速碰撞单层及双层薄靶板问题，通过与实验结果的对比，验证了OTM方法模拟超高速碰撞问题的有效性.

1 OTM方法

1.1 OTM方法理论

在参考构型中，连续介质运动过程由变形映射φ:Ω0×[t0,t]→3进行描述.假设材料的弹性响应、塑性响应以及比热容是相互独立的，因此系统的内能密度为

W=We,vol+We,dev+Wp+Wh,

(1)

式中：We,vol为弹性响应的体积部分;We,dev为弹性响应的偏离项;Wp为塑性响应;Wh为热储量.得到系统的内能密度之后，可以推出以下变分结构,

(2)

式中：K为动能密度；B为单位质量体积力；T为施加于Neumann边界Γt上的牵引力.该变分结构描述了材料系统的能量变化率.从变分原理出发,对等效势能Φ的变分取极值，再对该变分结构进行时间离散.在给定的时间增量[tn,tn+1]上，若tn时刻的状态变量φn已知，则tn+1时刻下的状态变量φn+1可以通过最小化增量能量近似获得，即

(3)

根据最优运输理论，时间间隔内增量动能的精确最小值在[tn,tn+1]上可以表示为初始质量密度和最终质量密度之间的Wasserstein距离[6]，即

(4)

变分结构(3)中除惯性项以外的其他项则采用后向差分格式进行近似，可以得到以下半离散化(时间离散)的增量变分结构，

(5)

对增量变分结构的时间离散格式(5)求极值，将在时间段[tn,tn+1]内产生动量守恒和能量守恒，使得系统的整个变形过程化为一系列离散的运动过程.

在OTM方法中采用物质点xp与节点xa相结合的方式对半离散增量变分公式进行无网格空间离散.计算过程中，材料的物理信息存储在物质点上,而计算邻域的动力学信息存储在节点上，如图1所示.

计算过程中，物质点的运动通过节点的动力学信息插值获得，采用标准的Ritz-Galerkin法来近似位移场.假设N表示离散域内的总节点数，则位移传输映射表示为

(6)

式中：Na(x)为节点xa的插值函数，插值函数Na(x)将采用局部最大熵无网格插值函数.LME插值函数支持宽度自适应改变，其导数满足连续性要求，在计算边界处具有Kronecker-delta性质，LME的具体理论可详见文献[5].同样地，基于Ritz-Galerkin法，物质点在tn+1时刻的位置通过节点近似获得，

(7)

式中：NH(xp,n)代表物质点xp在tn时刻下的邻域，其在每个时间步中将得到更新.

通过上述无网格近似方案对位移场进行近似，增量变分势能δΦn被完全地离散化，根据变分的极值条件，可得到系统的全离散格式力平衡方程，即

(8)

式中Ma,n+1代表节点xa在tn+1时刻的集中质量，即

(9)

式中mp为物质点xp所携带的质量，节点xa在tn+1时刻的加速度则通过中心差分的形式近似获得，因此节点坐标的更新形式为

xa,n+1=xa,n+(tn+1-tn)×

(10)

式中la,n为节点xa在tn+1时刻的动量，即

(11)

节点xa在tn+1时刻的节点内应力和节点外力分别为

(12)

(13)

通过求解每个节点上的力平衡方程(8)，可以获得整个系统的变形场.

1.2 数值迭代计算流程

基于OTM方法所编制程序的主要算法流程如图2所示.

① 几何模型的初始化.一般地，将d维的问题域(几何模型)Ω0进行网格划分生成初始节点集{xa,0,a=1,2,…,N}，并取每个单元的形心作为一个物质点的初始位置，组成物质点集{xp,0,p=1,2,…,M}，物质点的初始体积可设为该物质点所在单元占据的空间大小.

② 物质点和节点的初始化.在初始时刻n=0时，利用第1步离散得到的节点集合物质点集，初始化节点坐标xa,n、节点插值函数Na(xp,n)及其导数、节点力fa,n、节点动量la,n、节点质量矩阵Ma,n+1，初始化物质点坐标xp,n、物质点邻域NH(xp,n)、物质点体积vp,n、物质点密度ρp,n以及物质点质量mp=ρp,nvp,n.

③ 节点坐标更新.进行tn→tn+1的瞬态分析，求解全离散格式的力平衡方程(8)，或节点坐标更新公式(10)，显示计算得到节点在tn+1时刻的坐标{xa,n+1,a=1,2,…,N}.

④ 物质点的变形更新.根据节点坐标，此时节点就已产生一个位移，同时物质点将产生变形.采用Jaumann应力率模型更新物质点的应力偏张量.

⑤ 根据第3步得到的节点坐标，由式(7)更新物质点的坐标xp,n+1，更新物质点的体积vp,n+1=det(Jp,n→n+1)vp,n与密度ρp,n+1=mp/vp,n+1，其中mp为物质点xp的质量，在计算过程中保持不变，Jp,n→n+1为tn→tn+1时物质点xp的Jacobian矩阵.

⑥ 结合前两步得到的数据，通过本构关系获得物质点的应力σp,n+1，其中应力张量的偏张量部分由Jaumann应力率模型给出，球张量部分由物质点的密度ρp,n+1通过状态方程给出；更新tn+1时刻由于物质点变形引起的节点力fa,n+1，更新节点动量la,n及节点质量矩阵Ma,n+1.

⑦ 邻域和插值函数的更新.根据第5步得到的物质点坐标xp,n+1、体积vp,n+1与密度ρp,n+1，更新tn+1时刻下的邻域NH(xp,n+1)；同时更新节点的插值函数Na(xp,n+1)及其导数.

⑧ 在完成第7步后，即代表完成了物质点和节点在一个时间步内的力学动态分析，判断当前时间步n，如果已计算至最后一个时间步则结束计算，否则转入第3步.

2 弹丸超高速碰撞薄靶问题的数值模拟

2.1 弹丸和靶板的本构模型及状态方程

为描述超高速碰撞中高温高压高应变率对靶板材料屈服强度的影响，采用Johnson-Cook损伤模型[7]，即

σY=(1-D)(A+Bεn)×

(14)

D=∑(Δε/εf)，

(15)

损伤变量D将在每一个时间步上得到累积，其中εf是断裂塑性应变，与材料的应力三轴度、应变率和温度有关，其描述如下，

εf=[D1+D2exp(D3σ*)]×

(16)

式中：D1～D5为材料参数;σ*为应力三轴度;σm为平均正应力;σeq为Mises等效应力.

超高速碰撞问题中常用的状态方程为Mie-Grüneisen状态方程，表达式为

p=pH+Γ0ρ0(e-eH)，

(17)

式中：p和e分别为静水压力和比内能;ρ0为材料的初始密度；Γ0为材料的Grüneisen系数.Mie-Grüneisen状态方程的shock形式定义如下，

(18)

式中：μ=ρ/ρ0-1，pH和eH分别为Hugoniot曲线上静水压力和比内能的参考值；S为冲击波速度U与波后质点速度up之间线性关系的斜率，通常取U=C0+Sup；C0为体积声速.

2.2 弹丸超高速碰撞单层靶板

针对弹丸超高速撞击单层薄靶的过程，采用OTM方法进行模拟，数值工况参照Hiermaier等[8]的实验过程，弹丸和靶板皆采用LY12硬铝合金材料，材料模型采用Johnson-Cook损伤模型和Mie-Grüneisen状态方程，具体材料参数见表1、表2.

表1 本构模型参数

表2 状态方程参数

如图3所示，铝制球形弹丸半径为5 mm，铝靶板的尺寸为50 mm×50 mm×4 mm，铝球的撞击速度设为与实验[8]一致的6.18 km/s.计算过程中，时间步长由CFL条件给出，每个时间步长的量级约为0.1 μs，室温设为300 K.

图4为20 μs时OTM方法计算得到的碎片分布，将其与Hiermaier的实验[8]及Hiermaier数值模拟结果[8]进行对比，可以发现OTM得到的碎片云，整个碎片云的大部分质量中位于外泡碎片云的前端，气泡的大致形貌较为圆润，在形貌上和分布上都与实验结果吻合较好.

表3给出了OTM方法计算得到的碎片云形貌的具体特征参数，并与传统SPH方法[8]、ASPH方法[9]的模拟结果进行对比，其中：d1和d2分别为不包含和包含弹坑边缘粒子的弹坑直径；l为碎片云的膨胀长度；w为碎片云宽度；Δ为l/w与实验结果的相对误差.从表中可以看出，OTM方法得到的碎片云形貌与实验结果更加吻合，证明OTM方法非常适合于模拟超高速碰撞问题.

表3 不同数值方法的碎片云模拟结果对比

Tab.3 Comparison of different algorithms’ simulation result for debris cloud

数值方法D1D2l/mmw/mml/wΔ/%OTM方法模拟结果29.635.6102.478.21.315.8Hiermaier实验结果[2]27.534.51.39Hiermaier模拟结果[2]35.01.1120.1Liu模拟结果[7]28.9105.186.11.2212.2

2.3 弹丸超高速碰撞双层靶板

针对铝弹丸超高速碰撞双层铝板问题，采用OTM方法进行数值模拟，并与实验结果[10]相对比.实验中采用了9 mm直径的铝球以5 km/s的速度正碰撞双层铝靶板，第1层靶板厚1 mm，第2层靶板厚2 mm，两板的尺寸都设为50 mm×50 mm，两板相距17 mm.弹丸和两层靶板皆采用LY12硬铝合金材料，材料模型与材料参数的选取与2.2节相同.

图5给出了弹丸超高速碰撞双层靶板的计算结果，并与实验结果相比较.从图5中的结果对比可以看出：在弹丸撞击第1层板4 μs后，实验中的一次碎片云与后板开始接触，这一点与模拟结果相同；在一次碎片云撞击后板6.8 μs后，模拟得到的二次碎片云长度为15.8 mm，而实验结果显示为15.4 mm，相差2.6%，表明模拟结果与实验结果吻合良好.

3 结 论

本文采用最优输运无网格方法(OTM)，对铝制球形弹丸超高速撞击铝薄板形成碎片云的过程进行了数值模拟，得到以下结论.

① 在超高速撞击单层靶板的算例中，OTM方法计算得到的碎片云形态及内核碎片云的形状和分布均与实验吻合较好.整个碎片云的大部分质量中位于外泡碎片云的前端，气泡的大致形貌较为圆润，碎片云外围粒子分布较为平均.超高速撞击的弹坑直径、碎片云膨胀距离和宽度，与传统SPH方法、ASPH方法的模拟结果进行比较，OTM方法在碎片云的形貌上与实验结果更加接近，说明OTM方法适合模拟超高速碰撞等动态大变形问题.

② 针对球形弹丸超高速碰撞双层靶板的情况，采用OTM方法进行模拟，给出了几个关键时刻超高速碰撞过程的模拟结果.模拟结果显示，一次、二次碎片云的形貌与实验较为吻合，模拟给出的二次碎片云长度仅与实验结果相差2.6%，模拟结果良好.