@[TOC]D3plot文件数据提取及处理
LS-Reader是一套用于提供读取LS-DYNA计算结果接口的库，包括D3plot Reader、Binout Reader、LSDA Reader/Writer等模块，支持上千种仿真计算结果的提取和C、C++以及Python语言，下载地址为https://ftp.lstc.com/anonymous/outgoing/lsprepost/LS-Reader。
借助LS-Reader，对D3plot类型文件中的部分仿真结果进行了提取，包括节点坐标（D3P_NODE_COORDINATES）、节点速度（D3P_NODE_VELOCITIES）、单元节点索引（D3P_SHELL_CONNECTIVITY_MAT）、流体密度（D3P_2D_ALE_DENSITY）、单元内能密度（D3P_SHELL_INTERNAL_ENERGY_DENSITY）、物质体积分数（D3P_2D_ALE_VOLUME_FRACTION）和单元应力（D3P_TSHELL_STRESS）等。基于这些仿真结果，进一步处理后可获取冲击波波阵面、气泡体积、气泡内能和气泡边界处速度等，基本过程如下：
（1）冲击波波阵面
①根据单元节点索引获取单元对应节点的坐标，并将节点坐标(xn,yn)的平均值作为单元中心坐标(xe,ye)，即
式（1）
{
x
e
=
1
4
∑
k
=
1
4
x
n
k
y
e
=
1
4
∑
k
=
1
4
y
n
k
begin{cases} x_{mathrm{e}}=frac{1}{4}sum_{k=1}^4{x_{mathrm{n}k}}\ y_{mathrm{e}}=frac{1}{4}sum_{k=1}^4{y_{mathrm{n}k}}\end{cases}
{xe=41∑k=14xnkye=41∑k=14ynk
②D3plot Reader未提供直接读取流场压力的功能，而是提供了读取流场单元应力的接口，因此需要根据流场应力与压力的关系（式（2）），将应力转化为压力。结果如图1(a)所示，为验证结果是否正确，图1(b)给出了后处理软件LS-PrePost给出的结果。
式（2）
p
=
−
1
3
(
σ
x
x
+
σ
y
y
+
σ
z
z
)
p=-frac{1}{3}left( sigma _{xx}+sigma _{yy}+sigma _{zz} right)
p=−31(σxx+σyy+σzz)
式中，σxx、σyy和σzz为x轴、y轴和z轴方向的法向应力。
(a) 本文处理结果
（b）LS-PrePost处理结果
图1 压力云图
③由于单元空间排布规律未知，为方便处理数据，在极坐标系(θ,d)下，根据单元中心坐标和蝴蝶形网格空间分布特性为单元分配索引角度索引i和爆距索引j，计算公式如下：
式（3）
{
i
=
C
E
I
L
(
θ
n
a
n
g
)
j
=
C
E
I
L
(
log
δ
(
L
W
−
d
+
d
δ
n
a
n
g
L
W
)
)
begin{cases} i=mathrm{CEIL}left( frac{theta}{n_{mathrm{ang}}} right)\ j=mathrm{CEIL}left( log _{delta}left( frac{L_{mathrm{W}}-d+ddelta ^{n_{mathrm{ang}}}}{L_{mathrm{W}}} right) right)\end{cases}
⎩⎨⎧i=CEIL(nangθ)j=CEIL(logδ(LWLW−d+dδnang))
式中，CEIL为向上取整函数。
图2 θ=90°时不同索引处的压力
④对于每一个角度索引i，寻找波峰对应索引的集合{jp}以及计算波峰处相较于流体静压力的差值，将差值大于设定阈值所对应的索引中的最大值作为波阵面索引。图2给出了θ=90°时不同索引处的压力，同时标注了满足要求的索引。
⑤最后根据波阵面索引确定波阵面坐标以及压力值，结果如图3所示。
图3 冲击波波阵面形状及压力提取结果
（2）气泡体积
①根据单元对应四个节点的坐标，计算单元中心坐标以及该单元的面积。其中，单元中心坐标由式（1）给出，单元面积Se的计算公式如下：
式（4）
S
e
=
1
2
∣
∑
k
=
1
3
(
x
n
k
y
n
k
+
1
−
x
n
k
+
1
y
n
k
)
−
(
x
n
1
y
n
4
−
x
n
4
y
n
1
)
∣
S_{mathrm{e}}=frac{1}{2}left| sum_{k=1}^3{left( x_{mathrm{n}k}y_{mathrm{n}k+1}-x_{mathrm{n}k+1}y_{mathrm{n}k} right)}-left( x_{mathrm{n}1}y_{mathrm{n}4}-x_{mathrm{n}4}y_{mathrm{n}1} right) right|
Se=21∣∣∣∣∣k=1∑3(xnkynk+1−xnk+1ynk)−(xn1yn4−xn4yn1)∣∣∣∣∣
式中，节点按顺时针或逆时针排列。对于节点围成的凸四边形，可在极点平移到单元中心的极坐标系中对各节点按照极角的大小进行排序。
②根据单元面积以及单元中心坐标估算单元对应旋转体的体积，公式如下：
式（5）
V
e
≈
2
π
x
e
S
e
V_{mathrm{e}}approx 2pi x_{mathrm{e}}S_{mathrm{e}}
Ve≈2πxeSe
③以爆轰产物体积分数F为权重，将各单元对应的体积进行加权求和得到爆炸气泡的体积，即
式（6）
V
B
=
∑
k
=
1
N
(
V
e
k
F
k
)
V_{mathrm{B}}=sum_{k=1}^N{left( V_{mathrm{e}k}F_k right)}
VB=k=1∑N(VekFk)
式中，N为网格数。
为验证上述方法是否正确，图4中同时给出了借助观测点（by tracer）和根据气泡体积（by volume）获取的球形装药水下爆炸气泡半径-时间曲线。可以看出，除第三次脉动外，其他时刻根据气泡体积获取的气泡半径均与观测点获得的结果重合，即上述方法精度较高。实际上，由于误差的积累，气泡在第三次脉动开始时已经偏离球形（见图5），此时观测点获取的数据已不再可靠，因此两类数据不再重合。
图4 气泡半径-时间曲线计算结果对比
图5 t=0.198s时的气泡形态（x≥0一侧）
（3）气泡内能
①按照式（4）和式（5）估算单元旋转后对应的体积；
②D3plot Reader不能直接读取单元对应的内能，只能获取内能密度ei，对于以y轴为对称轴的二维轴对称模型，单元处单位体积内能eV由式（7）给出：
式（7）
e
V
=
e
i
x
e
e_V=frac{e_{mathrm{i}}}{x_{mathrm{e}}}
eV=xeei
因此，第k个单元对应的内能Eik为
式（8）
E
i
k
=
V
e
k
⋅
e
V
k
E_{mathrm{i}k}=V_{mathrm{e}k}cdot e_{Vk}
Eik=Vek⋅eVk
③以爆轰产物体积分数作为权重，将各单元对应内能进行求和得到气泡对应的内能，即
式（9）
E
i
B
=
∑
k
=
1
N
(
F
k
⋅
E
i
k
)
E_{mathrm{iB}}=sum_{k=1}^N{left( F_kcdot E_{mathrm{i}k} right)}
EiB=k=1∑N(Fk⋅Eik)
为验证计算结果是否正确，可将计算结果与初始时刻气泡内能进行对比。以球形装药为例，根据JWL状态方程参数可知1kg TNT爆炸气泡初始内能为EiB0=4.29448MJ，式（9）给出的计算结果为EiB0´=4.29446MJ，相对误差为10-6量级。
（4）气泡边界及速度
①计算单元中心坐标；
②寻找爆轰产物体积分数介于设定阈值范围内的单元的索引作为边界处单元索引index，例如：
式（10）
i
n
d
e
x
=
{
0
<
F
<
0.95
}
index=left{ 0
该方法较简单，而且能够提取复杂边界。不过，需要注意的是，在某些特殊时刻，气泡边界处单元的体积分数将恰好为1或者0，此时该方法将无法捕捉该单元，即该方法存在一定的缺陷。在实际操作中，该缺陷在球形装药中表现得比较明显，为避免通过少数点判断气泡外轮廓，在数据处理中建议剔除边界点较少时对应时刻的数据。
图6 气泡边界及速度
