超高速撞击下碎片云的OTM分析

廖祜明; 黎波; 樊江; 焦立新; 于帅超; 林健宇; 裴晓阳

doi:10.11883/bzycj-2021-0275

超高速撞击下碎片云的OTM分析

doi: 10.11883/bzycj-2021-0275

廖祜明^{1, 2,},
黎波³,
樊江^1, ,,
焦立新²,
于帅超²,
林健宇⁴,
裴晓阳⁴

1.
北京航空航天大学能源与动力工程学院，北京 102206
2.
云翼超算(北京)软件科技有限公司，北京 100027
3.
北京大学工学院，北京 100871
4.
中国工程物理研究院流体物理研究所，四川绵阳 621900

详细信息

作者简介:
廖祜明（1987－　），男，博士，liaohuming@buaa.edu.cn

通讯作者:
樊　江（1973－　），男，博士，副教授，fanjiang@buaa.edu.cn

中图分类号: O385
计量
- 文章访问数: 567
- HTML全文浏览量: 333
- PDF下载量: 111
- 被引次数: 1
出版历程
- 收稿日期: 2021-07-01
- 修回日期: 2022-08-26
- 网络出版日期: 2022-09-08
- 刊出日期: 2022-10-31

OTM analysis of debris cloud under hypervelocity impact

LIAO Huming^{1, 2
,},
LI Bo³,
FAN Jiang^{1
, ,},
JIAO Lixin²,
YU Shuaichao²,
LIN Jianyu⁴,
PEI Xiaoyang⁴

1.
School of Energy and Power Engineering, Beihang University, Beijng 102206, China
2.
ESCAAS Co. Ltd., Beijing 100027, China
3.
College of Engineering, Peking University, Beijing 100871, China
4.
Institute of Fluid Physics, China Academy of Engineering Physics, Sichuan Mianyang 621900, China

摘要

摘要: 空间碎片超高速撞击是典型的高温、高压、高应变率的极限力学问题，涉及材料复杂的动态响应，对传统的数值方法提出了巨大挑战。最优运输无网格（OTM）方法通过有机结合最优运输时间积分理论、局部最大熵无网格近似、物质点抽样、基于物理的裂纹扩展算法以及大规模并行计算策略，克服了传统数值方法瓶颈，在理论上保证了不同形式能量耗散的自主耦合分配，为超高速撞击仿真预测提供了高效的解决方案。采用基于OTM方法自主研发的极限力学仿真软件ESCAAS，对不同质量（3、10 g）的铜飞片以不同撞击角度（5.4°、11.7°）和不同撞击速度（5.55、5.12 km/s）撞击铝合金靶板的过程进行数值模拟，获得碎片云的形貌、靶板穿孔孔径等结果，与实验测量数据吻合良好，显示出OTM方法及ESCAAS软件可以作为超高速撞击的有力数值分析手段。
- 超高速撞击 /
- OTM方法 /
- 裂纹扩展 /
- 无网格方法
Abstract: The hypervelocity impact (HVI) of space debris is a typical extreme mechanics problem at high temperature, high pressure and high strain rate. The HVI involves the complex dynamic response of materials. Numerical methods have become very useful and important tools to predict the complex phenomena and look into the details of the entire process. However, the accurate simulation of the HVI is a grand challenge in scientific computing that places exacting demands on physics models, numerical solvers and computing resources. The optimal transportation meshfree (OTM) method is a meshfree updated-Lagrangian methodology for fluid and solid dynamic flows, possibly involving multiple phases, viscosity and general equations of state, general inelastic and history-dependent constitutive relations, arbitrary variable domains and boundary conditions and the interaction between fluid flows and highly deformable structures. The rationale behind the approach is combining concepts from the optimal transportation (OT) theory with material-point sampling, local maximum entropy (LME) approximation, the seizing contact, variational material point failure algorithm, and overcomes the essential difficulties in grid-based numerical methods like Lagrangian and Eulerian finite element method. Owing to those advantages, the OTM method provides an efficient and accurate solution for HVI simulation. In this paper, large scale three-dimensional numerical simulations on the HVI of the copper projectiles with different thicknesses (3.45, 5.13 mm), different masses (3, 10 g), different impact angles (5.4°, 11.7°) and different impact velocities (5.55, 5.12 km/s) impacting the Al6061-T6 plate with different thicknesses (2.87, 4.39 mm) were performed within the software ESCAAS based on the OTM method using a dynamic load balancing MPI/Pthreads parallel implementation. The dynamic response of material including phase transition in the high strain rate, high pressure and high temperature regime expected in this paper was described by the use of a variational thermomechanical coupling constitutive model with the SESAME equation of state, Grüneisen equation of state, rate-dependent J2 plasticity with power law hardening and thermal softening. The simulation results are in good agreement with the experimental measurements, which indicates the capacity of the OTM method and the ESCAAS software for HVI simulation.
- hypervelocity impact /
- OTM method /
- crack propagation /
- meshfree

HTML全文

地球轨道上散布着大量的空间碎片和微流星体，若它们与在轨运行航天器相撞，相对撞击速度可达十几千米每秒，对航天器造成巨大威胁。大型碎片撞击可直接撞毁航天器；中型碎片撞击会造成航天器部件功能失效；而小型碎片撞击可在航天器表面造成累积损伤，导致部件功能下降^[1]。因此，空间碎片超高速撞击防护是航天器结构设计须重点考虑的问题。厘清空间碎片超高速撞击的力学行为是防撞设计的重要基础。碎片超高速撞击的物理过程非常复杂，撞击过程中激波在材料中传播，同时激活各种相互耦合与竞争的能量耗散机制，引起材料的非弹性大变形、相变（液化、气化、凝固等）、裂纹萌生及扩展和层裂、碎裂等复杂的物理过程。采用数值方法模拟超高速撞击面临如何解决速度过大导致的网格畸变，以及如何准确描述超高速撞击过程中各种能量耗散模式之间的竞争与耦合关系等巨大挑战。

目前分析超高速撞击问题的数值计算方法主要包括基于网格的方法和无网格方法两大类。基于网格的数值方法中，传统拉氏有限元法受网格畸变困扰，难以处理超大变形问题；欧氏方法具有网格不变形的固有优势，但如何跟踪历史变量和物质界面仍然是亟待解决的问题，同时在模拟动态裂纹扩展、多体碰撞与应变率相关的材料热力学行为时还需进一步研究；任意拉格朗日-欧拉（arbitrary Lagrangian-Eulerian, ALE）^[2-3]法在一定程度上缓解了网格畸变带来的困难，但ALE继承了拉氏有限元的大部分缺点，同样面临大变形时网格畸变的问题^[4]。无网格法由于不需要进行网格离散及采用高阶插值函数，有利于解决大变形问题，在超高速撞击问题上应用更为广泛。常用的无网格方法有光滑粒子流体动力学法（smoothed particle hydrodynamics, SPH）^[5]、再生核粒子法（reproducing kernel particle method, RKPM）^[6]以及物质点法（material point method, MPM）^[7]等。有大量的研究应用以上方法^[8-9]及其改进方法^[10]对超高速撞击问题进行仿真分析并取得了丰富的研究成果，但在算法的精度、稳定性、健壮性、效率等方面仍然存在许多不足^[11]。

最优运输无网格方法（optimal transportation meshfree, OTM）^[12]是一种专门针对动态冲击问题提出的显式增量更新拉格朗日无网格方法。OTM方法采用局部最大熵无网格插值函数（local maximum entropy, LME）^[13]，在边界上满足Kronecker-Delta属性，解决了传统无网格法在处理位移边界条件时的缺陷；同时在节点的基础上增加物质点，以物质点作为积分点有效避免了计算结果在拉伸载荷下的不稳定性；在时间离散上采用了最优运输理论（optimal transport theory）^[14]，保证了其时间离散形式在理论上严格满足动量守恒，具有严格的收敛性数学证明；在并行计算上结合了共享内存多线程并行与分布式多进程并行实现了大规模多层混合并行^[15]，显著提高了计算效率；在材料失效方面，采用本征侵蚀裂纹扩展算法（Eigenerosion），以材料固有属性能量释放率作为失效判据 ^[16-18]，克服了传统裂纹扩展数值方法无法清晰表征材料真实变形和失效物理机制的缺点，理论证明收敛到Griffith理论，与网格无关。OTM方法的这些特性为超高速撞击模拟预测提供了高效精准的解决方案^[19]。本文采用基于OTM方法的极限力学仿真软件ESCAAS，对空间碎片超高速撞击问题进行数值模拟，并与实验结果进行对比验证。

1. OTM最优运输无网格方法

1.1 时间与空间离散

假设参考构型中的连续介质 ${\varOmega }_{0}\subset {\mathbb{R}}^{3}$ 运动过程由变形映射（位移） $\boldsymbol{\varphi }:{\varOmega }_{0}\times \left[{t}_{0},t\right]\to {\mathbb{R}}^{3}$ 描述，其中 $\left[{t}_{0},t\right]$ 为运动过程的总时间。 $\boldsymbol{X}\in {\varOmega }_{0}$ 表示材料质点在参考构形中的位置， $\boldsymbol{x}=\boldsymbol{\varphi }(\boldsymbol{X},t)$ 表示材料质点在现时构形 ${\varOmega }_{t}=\boldsymbol{\varphi }({\varOmega }_{0},t)$ 的位置，则拉格朗日描述下质量守恒、动量守恒与能量守恒方程为：

$\rho ({\boldsymbol{x}})J = {\rho _0}$

(1)

${\rho _0}\ddot {\boldsymbol{\varphi}} = {\nabla _0} \cdot ({{\boldsymbol{P}}^{\text{e}}} + {{\boldsymbol{P}}^{\text{v}}}) + {\rho _0}{\boldsymbol{B}}$

(2)

$T\dot N = {{\boldsymbol{P}}^{\text{v}}}:{\boldsymbol{F}} + {\boldsymbol{Y}} \cdot \dot {\boldsymbol{Z}} - {\nabla _0} \cdot {\boldsymbol{q}} + {\rho _0}Q$

(3)

式中： ${\rho }_{0}$ 为初始密度， $\rho$ 为当前密度， $\boldsymbol{F}$ 为变形梯度， $J$ 为体积变形率（Jacobian）， ${\boldsymbol{P}}^{\mathrm{e}}$ 为第一Piola-Kirchhoff应力张量平衡力项， ${\boldsymbol{P}}^{\mathrm{v}}$ 为第一Piola-Kirchhoff应力张量黏性力项， $\boldsymbol{B}$ 为体积力， $T$ 为温度、 $N$ 为单位体积的熵， $Q$ 为单位质量热源产生的热量， $\boldsymbol{q}$ 为热通量， $\boldsymbol{Z}$ 为内部变量， $\boldsymbol{Y}$ 为内部变量驱动力。

通过采用变分原理，建立与运动控制方程、初始条件与边界条件等效的变分结构^[20]，然后在变分框架下通过最小化变分结构来求解具有初始条件的边值问题。在连续介质力学中，与速率问题相对应的变分结构 $I\left(\dot{\boldsymbol{\varphi }}\right)$ 可表示为：

$I(\dot {\boldsymbol{\varphi}} ) = \mathop \int \nolimits_{{\varOmega _0}} \left( {\dot K + \dot W\left( {\boldsymbol{F}} \right) - {\rho _0}{\boldsymbol{B}} \cdot \dot {\boldsymbol{\varphi}} } \right){\text{d}}V - \mathop \int \nolimits_{{\varGamma _t}} \overline {\boldsymbol{T}} \cdot \dot {\boldsymbol{\varphi}} \ {\text{d}}S$

(4)

式中： $\dot{K}$ 为动能变化率， $\dot{W}$ 为应变能变化率， $\overline{\boldsymbol{T}}$ 为单位面积上施加的外力， $V$ 为体积， ${\varGamma }_{\mathrm{t}}$ 为Neumann边界， $S$ 为面积。根据最优运输理论，对式(4)进行时间离散后，变分可得：

$\begin{split} {\text{δ}}{I_k} =\int _{{\varOmega _k}} \frac{{2{\rho _k}}}{{{t_{k + 1}} - {t_k}}}\left( {\frac{{{{\boldsymbol{\varphi}} _{k \to k + 1}}\left( {\boldsymbol{x}} \right) - {\boldsymbol{x}}}}{{{t_{k + 1}} - {t_k}}} + \frac{{{\boldsymbol{\varphi}} _{k - 1 \to k}^{ - 1}\left( {\boldsymbol{x}} \right) - {\boldsymbol{x}}}}{{{t_k} - {t_{k-1}}}}} \right) \cdot {\text{δ }}{{\boldsymbol{\varphi}} _{k + 1}}{\text{d}}V +\\ \int _{{\varOmega _k}} {{\boldsymbol{P}}_{k + 1}}:{\text{δ}}{{\boldsymbol{F}}_{k + 1}}{\text{d}}V - \int_{{\varOmega _k}} {\rho _k}{{\boldsymbol{B}}_{k + 1}} \cdot {\text{δ}}{{\boldsymbol{\varphi}} _{k + 1}}{\text{d}}V - \int_{{\Gamma _t}} {{\overline {\boldsymbol{T}}}_{k + 1}} \cdot {\text{δ}}{{\boldsymbol{\varphi}} _{k + 1}}{\text{d}}S \end{split}$

(5)

式中：下标 $k$ 代表时间步长索引。式(5)为时间离散的变分公式，为了便于计算机求解，要对此半离散的变分公式进行空间离散得到时间-空间全离散格式。

在有限元法中，每个单元有特定的形状且各单元间具有固定的连接关系，每个单元携带了材料一定量的质量并具有相应的体积，材料的响应在单元的积分点上进行求解。与有限元法类似，OTM方法中采用物质点 ${\boldsymbol{x}}_{p}$ 与节点 ${\boldsymbol{x}}_{a}$ 相结合的方式对半离散增量变分公式进行无网格空间离散，其中下标 $a$ 表示节点索引， $p$ 表示物质点索引。如图1所示，所有的场变量信息分别由节点 ${\boldsymbol{x}}_{a}$ 和物质点 ${\boldsymbol{x}}_{p}$ 携带。物质点可以理解为有限元中的积分点，而其携带的体积可以看作是有限元中的单元，而不同于有限元的是节点和物质点不再有固定的连接关系，即此“单元”没有固定的形状，其形状将根据物质的变形可以是任意的形状，由于解放了网格的束缚，因此对于大变形问题不存在网格畸变带来的问题^[21]。

图 1 空间离散示意图^[21]

Figure 1. Spatial discrete schematic diagram^[21]

下载: 全尺寸图片幻灯片

材料的动力学信息存储在节点上，包括位移、速度、加速度与温度等。在计算过程中 ${t}_{k+1}$ 时刻节点位置 ${\boldsymbol{x}}_{a,k+1}$ 和节点速度 ${\dot{\boldsymbol{\varphi }}}_{a,k+1}$ 的更新公式为：

${{\boldsymbol{x}}_{a,k + 1}} = {{\boldsymbol{x}}_{a,k}} + {\dot {\boldsymbol{\varphi}} _{a,k}}\Delta t + \frac{1}{2}{\ddot {\boldsymbol{\varphi}} _{a,k}}{{\Delta }}{t^2}$

(6)

${\dot {\boldsymbol{\varphi}} _{a,k + 1}} = {\dot {\boldsymbol{\varphi}} _{a,k}} + \frac{1}{2}\left( {{{\ddot {\boldsymbol{\varphi}} }_{a,k}} + {{\ddot {\boldsymbol{\varphi}} }_{a,k + 1}}} \right)\Delta t$

(7)

从 ${t}_{k}$ 时刻到 ${t}_{k+1}$ 时刻的位移传输（变形）映射 ${\boldsymbol{\varphi }}_{k\to k+1}$ 为：

${{\boldsymbol{\varphi}} _{k \to k + 1}} = \mathop \sum \nolimits_{a \in {N_{\rm{H}}}\left( {{{\boldsymbol{x}}_{p,k}}} \right)} {{\boldsymbol{x}}_{a,k + 1}}{N_a}\left( {{{\boldsymbol{x}}_{p,k}}} \right)$

(8)

式中： ${N}_{a}$ 为LME局部最大熵插值函数， ${N}_{\mathrm{H}}\left({\boldsymbol{x}}_{p,k}\right)$ 代表物质点 ${\boldsymbol{x}}_{p}$ 在 ${t}_{k}$ 时刻下的邻域。

材料的物理信息存储在物质点上，包括质量、体积、密度、变形、应力与材料内部参数等。材料的响应在物质点上求解，物质点携带的质量在计算过程中保持不变，而体积和密度在不断的变化，每个时间步都得到更新。物质点的运动通过节点的动力学信息插值获得^[22]。在计算过程中，物质点在 ${t}_{k+1}$ 的位置更新，通过对邻域内的节点进行插值获得：

${{\boldsymbol{x}}_{p,k + 1}} = \mathop \sum \nolimits_{a \in {N_H}\left( {{{\boldsymbol{x}}_{p,k}}} \right)} {{\boldsymbol{x}}_{a,k + 1}}{N_a}\left( {{{\boldsymbol{x}}_{p,k}}} \right)$

(9)

物质点邻域范围 ${N}_{\mathrm{H}}\left({\boldsymbol{x}}_{p,k}\right)$ 在每个时间步进行更新，同时更新各节点上插值函数值。LME插值函数 ${N}_{a}$ 满足严格的非负性，以及零阶和一阶连续性要求，同时在边界满足Kronecker-Delta属性。此外，在LME插值函数中，通过调整 $\mathrm{\gamma }$ 值，插值函数作用范围可从局部有限元无缝过渡到全局无网格，如图2所示，这种衰减特性建立了与高斯径向基函数的联系，其重要作用是只有少量的节点对待求函数有明显的贡献，大大降低计算成本。

图 2 局部最大熵插值函数^[22]

Figure 2. Local maximum entropy shape function^[22]

下载: 全尺寸图片幻灯片

通过将上述离散的参数代入半离散变分公式中，并实施驻点条件，可以获得完全离散的应力平衡方程：

${m_{a,k + 1}}{\ddot {\boldsymbol{\varphi}} _{a,k + 1}} = {\boldsymbol{f}}_{a,k + 1}^{{\text{int}}} - {\boldsymbol{f}}_{a,k + 1}^{{\text{ext}}}$

(10)

式中： ${m}_{a,k+1}$ 为节点质量， ${\boldsymbol{f}}_{a,k+1}^{\mathrm{i}\mathrm{n}\mathrm{t}}$ 和 ${\boldsymbol{f}}_{a,k+1}^{\mathrm{e}\mathrm{x}\mathrm{t}}$ 分别为节点内力和节点合外力。关于OTM基础理论更为详细的介绍和收敛性证明见文献[12]。

1.2 捕捉式多体动态接触算法

OTM方法采用物质点和节点进行离散，并且动态更新局部邻域，使其具备自动捕获接触（seizing contact）的能力。具体地，物体 ${\varOmega }_{\mathrm{A}}$ 上的节点同时也可以作为物体 ${\varOmega }_{\mathrm{B}}$ 中物质点的邻域，因此对两物体的接触可以不加特殊处理，如图3所示，其中 ${R}_{a}$ 为节点 ${\boldsymbol{x}}_{a}$ 的插值函数的作用范围， ${R}_{b}$ 同理， ${\boldsymbol{v}}_{p,k}$ 代表物质点速度， ${\boldsymbol{v}}_{b,k}$ 代表节点速度。

图 3 捕捉式接触算法及节点邻域示意图^[23]

Figure 3. Schematic of the seizing contact algorithm and the support of nodes ^[23]

下载: 全尺寸图片幻灯片

在捕捉式多体动态接触算法中，当碰撞节点进入物质点邻域后，其线性动量会被自然地抵消，如碰撞节点 ${\boldsymbol{x}}_{a}$ 的线性动量 ${\boldsymbol{l}}_{a,k}$ 为：

${{\boldsymbol{l}}_{a,k}} = \sum\limits_{p = {\text{1}}}^M {{{\boldsymbol{m}}_p}\frac{{{{\boldsymbol{x}}_{p,k}} - {{\boldsymbol{x}}_{p,k - {\text{1}}}}}}{{{{\boldsymbol{t}}_k} - {{\boldsymbol{t}}_{k - {\text{1}}}}}}{N_a}\left( {{{\boldsymbol{x}}_{p,k}}} \right)}$

(11)

式中： $M$ 为节点邻域内的物质点数量，为 ${m}_{p}$ 为物质点质量。节点的动量为其邻域物质点动量的平均值，结合运算中采用的集中质量法，节点的速度本身代表了局部物质点速度的质量加权平均值，平均了速度的所有分量，接触时物体之间的法向速度和切向速度都被用于相互抵消，从而自动确定了碰撞体之间的接触条件。因此，这种捕捉式接触时的摩擦系数可视为无限大，适用于模拟碰撞终端抛射物与目标物混合及伴随的绝热加热过程。

1.3 本征侵蚀裂纹扩展算法

材料在极限条件下的裂纹扩展过程包含了多种能量传播与耗散方式的结合与竞争，准确地描述这些特征需要采用基于物理的裂纹扩展算法。OTM方法中采用能量的方式来描述整个材料的响应过程，提出了本征侵蚀（Eigenerosion）裂纹扩展算法。在该方法中，每个物质点代表一小块物质，通过移除失效物质点的方式来模拟材料中裂纹扩展和碎裂的形成，物质点失效代表该小块物质内部产生裂纹并形成自由表面，而自由表面的产生则必然伴随着能量耗散，裂纹路径是由弹性能量（弹性能量作为能量释放机制促进断裂）、材料特定的断裂能量（与裂纹面积成正比的断裂比能）、单调性和接触约束之间的竞争引起的结果。物质点失效的判断准则为：

${G_p} = \frac{{\lambda \eta }}{{\left\| {{C_\eta }} \right\|}}\mathop \sum \limits_{{{\boldsymbol{x}}_q} \in {C_\eta }({{\boldsymbol{x}}_p})} {W^{\text{e}}}\left( {{{\boldsymbol{x}}_q}} \right){V_q} {\text{≥}} {G_{\rm c}}$

(12)

式中： ${G}_{p}$ 表示物质点 ${\boldsymbol{x}}_{p}$ 的等效能量释放率，下标 $p$ 代表计算域内物质点的索引； ${C}_{\eta }$ 表示局部邻域， $\eta$ 决定了其范围， $\|{C}_{\eta }\|$ 表示局部邻域总体积； $\lambda$ 为几何校正因子； ${V}_{q}$ 为局部邻域内物质点 ${\boldsymbol{x}}_{q}$ 体积，下标 $q$ 代表邻域内物质点的索引。如图4所示，黑色点为一组失效物质点组成的裂纹，圆圈内的红点为位于裂纹尖端的物质点的邻域 $\eta$ ，通过平均该物质点邻域内的弹性能 ${W}^{\mathrm{e}}$ 得到物质点的等效能量释放率 ${G}_{p}$ ，若 ${G}_{p}$ 大于材料本身的临界能量释放率 ${G}_{\mathrm{c}}$ 时，该物质点将失效不再参与计算。

图 4 EigenErosion算法等效能量释放率计算示意图 ^[24]

Figure 4. Schematic of the equivalent energy release rate calculation of the EigenErosion algorithm ^[24]

下载: 全尺寸图片幻灯片

由于临界能量释放率 ${G}_{\mathrm{c}}$ 为材料固有参数，可通过实验测得的断裂韧性（fracture toughness）转换得出。本征侵蚀裂纹扩展算法不需要对裂纹的方向、大小进行显式描述，具备简便的三维几何与拓扑结构的处理方式，通过平均化能量克服了由于网格分布形式而带来的收敛性问题和裂纹路径网格相关的问题，且该算法有严格的数学证明收敛于Griffith解^[16]。此外，由于本征侵蚀裂纹扩展方法从材料失效物理机制出发，考虑了材料内部能量耗散的各种机制，包括塑性变形、裂纹扩展和相变，使得精确预测脆性或者韧性材料在不同温度以及载荷下的裂纹扩展成为可能。

2. 铜飞片超高速撞击铝合金靶板数的值模拟与验证

2.1 数值模拟模型

文献[25]开展了柱形OFHC铜飞片超高速撞击Al6061-T6铝合金靶板实验，并运用高速摄像系统拍摄了碎片云形貌。本节将采用基于OTM方法开发的ESCAAS软件，建立同样的超高速撞击模型（如图5所示），对该超高速撞击过程进行模拟并开展与实验结果的对比验证。

图 5 计算模型示意图

Figure 5. Calculation models schematic diagram

下载: 全尺寸图片幻灯片

工况1：撞击角度为5.4°，质量为3 g的柱形OFHC铜飞片，以5.55 km/s的速度撞击Al6061-T6铝合金靶板，几何模型如图5(a)所示，靶板外侧施加固定边界条件，取铜飞片和靶板的1/2对称模型进行无网格离散，物质点最小特征尺寸为0.14 mm，共计31.3万个物质点，6.26万个节点；

工况2：撞击角度为11.7°，质量为10 g的柱形OFHC铜飞片，以5.12 km/s的速度撞击Al6061-T6铝合金靶板，几何模型如图5(b)所示，靶板外侧施加固定边界条件；取铜飞片和靶板的1/2对称模型进行无网格离散，物质点最小特征尺寸为0.6 mm，共计40.6万个物质点，22.9万个节点，如图5(c)所示。

2.2 材料模型

超高速撞击过程涉及材料应变硬化、应变率硬化和高温软化等效应，J2黏塑性模型（J2-power law viscoplasticity）能较好地模拟铜和铝合金材料塑性响应，其等效屈服应力为：

${\sigma _{\text{y}}}\left( {{\varepsilon ^{\text{p}}},{{\dot \varepsilon }^{\text{p}}},T} \right) = {\sigma _0}{\left( {1 + \frac{{{\varepsilon ^{\text{p}}}}}{{\varepsilon _0^{\text{p}}}}} \right)^n}\left[{1 + {{\left( {\frac{{{{\dot \varepsilon }^{\text{p}}}}}{{\dot \varepsilon _0^{\text{p}}}}} \right)}^m}} \right]{\left( {1 - \frac{{T - {T_0}}}{{{T_{\text{m}}} - {T_0}}}} \right)^l}$

(13)

式中： ${\varepsilon }^{\mathrm{p}}$ 和 ${\dot{\varepsilon }}^{\mathrm{p}}$ 分别为塑性应变与塑性应变率； ${\sigma }_{0}$ 为参考温度下的初始准静态屈服应力； ${\varepsilon }_{0}^{\mathrm{p}}$ 和 ${\dot{\varepsilon }}_{0}^{\mathrm{p}}$ 分别为参考塑性应变和参考塑性应变率； $n$ 为硬化指数， $m$ 为率敏感指数， $l$ 为热软化指数； $T$ 为绝对温度， ${T}_{0}$ 、 ${T}_{\mathrm{m}}$ 分别为材料初始温度与熔化温度； ${\sigma }_{0}$ 、 ${\varepsilon }_{0}^{\mathrm{p}}$ 、 ${\dot{\varepsilon }}_{0}^{\mathrm{p}}$ 、 $m$ 、n和l为待拟合参数。本算例柱形OFHC铜飞片和Al6061-T6铝合金靶板的材料物性参数如表1所示，其J2材料本构模型参数如表2所示。高温高压下Al6061-T6铝合金变形与温度、压力的关系采用SESAME状态数据库处理，OFHC铜采用Grüneisen状态方程进行描述，材料数据引自文献[26]。

表 1 材料物性参数

Table 1. Material parameters

材料	密度/ (kg·m⁻³)	杨氏模量/ GPa	泊松比	临界能量释放率/ (kJ·m⁻²)
OFHC铜	8930.0	129.0	0.35	250.0
Al6061-T6铝合金	2700.0	68.9	0.33	100.0

下载: 导出CSV

| 显示表格

表 2 材料本构模型参数

Table 2. Parameters of material constitutive model

材料	σ₀/MPa	${\varepsilon }_{0}^{\mathrm{p}}$	${\dot{\varepsilon }}_{0}^{\mathrm{p}}$	n	m	l	T_m/K
OFHC铜	120.0	0.0278	1.0	0.45	1.0	1.0	1790.0
Al6061-T6铝合金	270.0	0.0025	1.0	0.148	0.01389	1.0	1200.0

下载: 导出CSV

| 显示表格

2.3 结果分析

工况1：在Linux平台下采用ESCAAS模拟 3 g柱形铜飞片，以5.4°的撞击角度、5.55 km/s的速度撞击铝合金靶板的过程。整个模拟采用16核并行计算，10 μs撞击过程模拟用时约为1.33 h。图6为铜飞片穿透缓冲墙形成碎片的过程，模拟结果很好地呈现出碎片云的形态，包括碎片云呈现出头部椭圆形、内核碎片云的位置、外泡碎片云的形态以及反溅碎片云的形态等特征信息。

图 6 工况1不同时刻下碎片云轮廓

Figure 6. Debris cloud outlines at different moments of 3 g copper impact (condition 1)

下载: 全尺寸图片幻灯片

图7所示为OTM模拟与实验测量碎片云形貌的对比。轴向位置点Point 1为靶板碎片云前端；轴向位置点Point 2为铜飞片碎片云前端；轴向点Point 3为撞击之前靶板的前端；径向位置点Point 4和点Point 5为铜飞片和靶板碎片径向相交距离测量点；径向位置点Point 6为靶板后端反溅碎片云高度。从对比结果中可以看到OTM模拟得到的碎片云形貌特征与实验结果吻合理想。其中，Point 2位置处实验测量的碎片云前端轮廓较为光滑和饱满外凸，而OTM模拟的碎片云前端在Point 2位置处相对平坦和稀疏，其原因是此算例中网格量较少，进一步加密网格可以改善效果。此外，值得注意的是由于实验中碎片云图像是采用激光阴影技术获得，细微的碎片无法在照片中显示，因而实验中的碎片云轮廓较为清晰，而OTM计算得到的碎片云对比图中，所有碎片均有显示，轮廓有一定发散性。

图 7 工况1碎片云形貌对比情况 (6.4 μs)

Figure 7. Debris cloud shape comparison at 6.4 μs (condition 1)

下载: 全尺寸图片幻灯片

除了碎片云轮廓形貌的定性对比外，碎片云特征尺寸定量对比方面：铝合金靶板前端碎片云轮廓轴向最大位置点Point 1距靶板的轴向距离的OTM模拟结果为(50.74±0.66) mm，实验测量结果为51.3 mm，与实验的相对误差为(1.09±1.29)%；铜飞片前端碎片云轮廓轴向最大位置点Point 2距靶板的轴向距离的OTM模拟结果为(31.9±1.4) mm，实验测量结果为33.5 mm，与实验的相对误差为(4.8±4.18)%；铝合金靶板反溅碎片云的高度测量点Point 6距靶板的轴向距离的OTM模拟结果为(7.87±0.74) mm，实验测量结果为7.76 mm，与实验对比相对误差为(1.42±9.54)%；铜飞片和靶板材料碎片径向相交距离（测量点Point 4和Point 5的径向位置）的OTM模拟结果为(27.6±1.3) mm，实验测量结果为27.9 mm，与实验对比相对误差为(1.07±4.66)%。

在超高速撞击过程中，弹丸厚度/板厚(D/H)、打击速度与角度的结合决定了冲击波的分段分区效应。本例中D/H=1.2，属于中厚板情形，其主要失效模式为塑性变形、裂纹扩展、熔化三种能量转换与耗散方式相结合。如图8所示，OTM方法模拟了飞片撞击下，靶板发生塑性变形引起靶板翻边撕裂的过程。测量得到OTM模拟的穿孔直径为(20.9±0.4) mm，实验测试值约为20.1 mm，相对误差为(3.98±1.99)%。从上述分析可见仿真量化结果与实验吻合良好。

图 8 工况1靶板穿孔直径

Figure 8. Perforation diameter of target plate (condition 1)

下载: 全尺寸图片幻灯片

工况2：模拟了10 g柱形OFHC铜飞片，以11.7°的撞击角度，5.12 km/s的速度撞击Al6061-T6铝合金靶板的过程。整个模拟采用16核并行计算，10 μs撞击过程模拟时间约为2.4 h。图9所示为不同时刻下碎片云的轮廓，由于铜飞片撞击角度11.7°较大，铜飞片碎片云与靶板碎片云前端呈现出显著偏转。

图 9 工况2不同时刻下碎片云轮廓

Figure 9. Debris cloud outlines at different moments of 10 g copper impcat (condition 2)

下载: 全尺寸图片幻灯片

图10为7.6 μs时刻下，OTM模拟与实验测量的碎片云形貌对比情况。同样地，轴向位置点Point 1为靶板碎片云前端；轴向位置点Point 2为铜飞片碎片云前端；轴向点Point 3为撞击之前靶板的前端；径向位置点Point 4和点Point 5为铜飞片和靶板碎片径向相交距离测量点；径向位置点Point 6为靶板后端反溅碎片云高度。此外，由于撞击角度11.7°较大，铜飞片碎片云与靶板碎片云前端呈现出显著偏转，增加了碎片云偏转角度与实验测量结果的对比情况，包括碎片云前端偏离中线的角度 $\alpha$ ，靶板上侧反溅碎片云偏离角 $\theta$ ，靶板下侧反溅碎片云偏离角 $\beta$ 。

图 10 工况2碎片云形貌对比情况 (7.6 μs)

Figure 10. Debris cloud shape comparison at 7.6 μs (condition 2)

下载: 全尺寸图片幻灯片

铜飞片速度为5.12 km/s时，铝合金靶板前端碎片云轮廓轴向最大位置点Point 1距靶板轴向距离的OTM模拟结果为(47.2±1.2) mm，实验测量结果为47.5 mm，与实验的相对误差为(0.63±2.53)%；铜飞片前端碎片云轮廓轴向最大位置点Point 2距靶板轴向距离的OTM模拟结果为(27.1±2.1) mm，实验测量结果为26.0 mm，与实验的相对误差为(4.23±8.08)%；铝合金靶板反溅碎片云的高度测量点Point 6距靶板轴向距离的OTM模拟结果为(32.0±2.6) mm，实验测量结果为30.9 mm，与实验对比相对误差为(3.56±8.41)%；碎片云前端偏离中线的角度 $\alpha$ 的OTM模拟结果为(10.55±0.65)°，实验测量结果为11.5°，与实验对比相对误差约为(8.26±5.65)%；靶板上侧反溅碎片云偏离角 $\theta$ 的OTM模拟结果为(109.5±3.5)°，实验测量结果为100.8°，与实验对比相对误差约为(8.63±3.47)%；靶板下侧反溅碎片云偏离角 $\beta$ 的OTM模拟结果为(109.75±3.25)°，实验测量结果为108.0°，与实验对比相对误差约为(1.62±3.00)%，模拟结果与实验结果吻合理想。

图11所示为工况2铜飞片撞击铝合金靶板，本例中D/H=1.17属于中厚板，同样地，在撞击过程中可以观察到铝合金靶板具有显著的塑性变形和翻边效果。OTM模拟的靶板穿孔直径为32.125 mm，实验测试值约为31.6 mm，相对误差大约1.66%，穿孔孔径与实验高度吻合。

图 11 工况2靶板穿孔直径

Figure 11. Perforation diameter of the target plate (condition 2)

下载: 全尺寸图片幻灯片

3. 结　论

超高速撞击涉及材料局部超大变形、裂纹扩展、层裂与碎裂、多体动态接触、固液气相变与多相混合等复杂的材料动态响应过程，对传统基于网格的数值方法提出了巨大挑战。本文引入最优运输无网格方法（OTM）对空间碎片超高速撞击过程进行数值分析，指出了OTM在处理此类问题中有如下优势：

(1) OTM方法基于变分原理，在统一的拉格朗日框架下求解物质运动方程，实现了对不同材料在高温、高压、高应变率等极限条件下多种能量耗散模式结合与竞争的现象进行了完善的考虑，在理论上保证了不同形式能量耗散的自主耦合分配；

(2) OTM方法采用物质点和节点空间离散，在时间与空间上的离散被严格地限制在变分框架内，保证了质量，能量及动量的守恒；结合局部最大熵无网格插值函数，克服了一般无网格法在精准施加位移边界条件上的不足，稳定性和收敛性有完备的数学证明；

(3) OTM方法采用了真实的裂纹扩展算法，考虑了材料内部能量耗散的各种机制，包括塑性变形、裂纹扩展和相变，采用材料的固有属性能量释放率作为失效判据，而非简单采用单元变形量、应变大小，或是基于实验的失效模型来作为材料的失效判据，克服了传统裂纹扩展数值方法无法清晰表征材料真实变形和失效物理机制，以及不收敛和网格相关的缺点；

(4) OTM方法及ESCAAS软件在计算弹孔直径、碎片云形态方面与实验结果吻合理想，验证了其在超高速撞击数值模拟方面的适用性，显示出OTM方法及ESCAAS软件可以作为超高速撞击的有力数值分析手段。

图 1 空间离散示意图^[21]

Figure 1. Spatial discrete schematic diagram^[21]

下载: 全尺寸图片幻灯片

图 2 局部最大熵插值函数^[22]

Figure 2. Local maximum entropy shape function^[22]

下载: 全尺寸图片幻灯片

图 3 捕捉式接触算法及节点邻域示意图^[23]

Figure 3. Schematic of the seizing contact algorithm and the support of nodes ^[23]

下载: 全尺寸图片幻灯片

图 4 EigenErosion算法等效能量释放率计算示意图 ^[24]

Figure 4. Schematic of the equivalent energy release rate calculation of the EigenErosion algorithm ^[24]

下载: 全尺寸图片幻灯片

图 5 计算模型示意图

Figure 5. Calculation models schematic diagram

下载: 全尺寸图片幻灯片

图 6 工况1不同时刻下碎片云轮廓

Figure 6. Debris cloud outlines at different moments of 3 g copper impact (condition 1)

下载: 全尺寸图片幻灯片

图 7 工况1碎片云形貌对比情况 (6.4 μs)

Figure 7. Debris cloud shape comparison at 6.4 μs (condition 1)

下载: 全尺寸图片幻灯片

图 8 工况1靶板穿孔直径

Figure 8. Perforation diameter of target plate (condition 1)

下载: 全尺寸图片幻灯片

图 9 工况2不同时刻下碎片云轮廓

Figure 9. Debris cloud outlines at different moments of 10 g copper impcat (condition 2)

下载: 全尺寸图片幻灯片

图 10 工况2碎片云形貌对比情况 (7.6 μs)

Figure 10. Debris cloud shape comparison at 7.6 μs (condition 2)

下载: 全尺寸图片幻灯片

图 11 工况2靶板穿孔直径

Figure 11. Perforation diameter of the target plate (condition 2)

下载: 全尺寸图片幻灯片

表 1 材料物性参数

Table 1. Material parameters

材料	密度/ (kg·m⁻³)	杨氏模量/ GPa	泊松比	临界能量释放率/ (kJ·m⁻²)
OFHC铜	8930.0	129.0	0.35	250.0
Al6061-T6铝合金	2700.0	68.9	0.33	100.0

下载: 导出CSV

表 2 材料本构模型参数

Table 2. Parameters of material constitutive model

材料	σ₀/MPa	${\varepsilon }_{0}^{\mathrm{p}}$	${\dot{\varepsilon }}_{0}^{\mathrm{p}}$	n	m	l	T_m/K
OFHC铜	120.0	0.0278	1.0	0.45	1.0	1.0	1790.0
Al6061-T6铝合金	270.0	0.0025	1.0	0.148	0.01389	1.0	1200.0

下载: 导出CSV

参考文献(26)

[1]	刘岩, 张雄, 刘平, 等. 空间碎片防护问题的物质点无网格法与软件系统 [J]. 载人航天, 2015, 21(5): 503–509. DOI: 10.16329/j.cnki.zrht.2015.05.013. LIU Y, ZHANG X, LIU P, et al. Meshfree material point method and software system for problems of shielding space debris [J]. Manned Spaceflight, 2015, 21(5): 503–509. DOI: 10.16329/j.cnki.zrht.2015.05.013.
[2]	DONEA J, GIULIANI S, HALLEUX J P. An arbitrary Lagrangian-Eulerian finite element method for transient dynamic fluid-structure interactions [J]. Computer Methods in Applied Mechanics and Engineering, 1982, 33(1): 689–723. DOI: 10.1016/0045-7825(82)90128-1.
[3]	DONEA J, HUERTA A, PONTHOT J P, et al. Arbitrary Lagrangian-Eulerian methods [M]//STEIN E, DE BORST R, HUGHES T J R. Encyclopedia of Computational Mechanics. John Wiley, 2004: 413–437. DOI: 10.1002/0470091355.ecm009.
[4]	QUAN X, BIRNBAUM N K, COWLER M S, et al. Numerical simulation of structural deformation under shock and impact loads using a coupled multi-solver approach [C]// 5th Asia-Pacific Conference on Shock and Impact Loads on Structures. Hunan, 2003.
[5]	GINGOLD R A, MONAGHAN J J. Smoothed particle hydrodynamics: theory and application to non-spherical stars [J]. Monthly Notices of the Royal Astronomical Society, 1997, 181(3): 375–389. DOI: 10.1093/mnras/181.3.375.
[6]	LIU W K, JUN S, ZHANG Y F. Reproducing kernel particle methods [J]. International Journal for Numerical Methods in Fluids, 1995, 20(8/9): 1081–1106. DOI: 10.1002/fld.1650200824.
[7]	ZHANG X, CHEN Z, LIU Y. The material point method [M]//ZHANG X, CHEN Z, LIU Y. The Material Point Method: A Continuum-Based Particle Method for Extreme Loading Cases. Oxford: Academic Press, 2017: 37-101. DOI: 10.1016/B978-0-12-407716-4.00003-X.
[8]	闫晓军, 张玉珠, 聂景旭. 空间碎片超高速碰撞数值模拟的SPH方法 [J]. 北京航空航天大学学报, 2005, 31(3): 351–354. DOI: 10.3969/j.issn.1001-5965.2005.03.019. YAN X J, ZHANG Y Z, NIE J X. Numerical simulation of space debris hypervelocity impact using SPH method [J]. Journal of Beijing University of Aeronautics and Astronautics, 2005, 31(3): 351–354. DOI: 10.3969/j.issn.1001-5965.2005.03.019.
[9]	刘有英, 王海福. 高速碰撞下航天器防护结构效能评价 [J]. 弹箭与制导学报, 2005, 25(4): 359–361. DOI: 10.3969/j.issn.1673-9728.2005.04.117. LIU Y Y, WANG H F. Evaluations of high-velocity impact for spacecraft shields [J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2005, 25(4): 359–361. DOI: 10.3969/j.issn.1673-9728.2005.04.117.
[10]	强洪夫, 范树佳, 陈福振, 等. 基于拟流体模型的SPH新方法及其在弹丸超高速碰撞薄板中的应用 [J]. 爆炸与冲击, 2017, 37(6): 990–1000. DOI: 10.11883/1001-1455(2017)06-0990-11. QIANG H F, FAN S J, CHEN F Z, et al. A new smoothed particle hydrodynamics method based on the pseudo-fluid model and its application in hypervelocity impact of a projectile on a thin plate [J]. Explosion and Shock Waves, 2017, 37(6): 990–1000. DOI: 10.11883/1001-1455(2017)06-0990-11.
[11]	林健宇, 罗斌强, 徐名扬, 等. 铝弹丸超高速撞击防护结构的研究进展 [J]. 高压物理学报, 2019, 33(3): 030112. DOI: 10.11858/gywlxb.20190774. LIN J Y, LUO B Q, XU M Y, et al. Progress of aluminum projectile impacting on plate with hypervelocity [J]. Chinese Journal of High Pressure Physics, 2019, 33(3): 030112. DOI: 10.11858/gywlxb.20190774.
[12]	LI B, HABBAL F, ORTIZ M, et al. Optimal transportation meshfree approximation schemes for fluid and plastic flows [J]. International Journal for Numerical Methods in Engineering, 2010, 83(12): 1541–1579. DOI: 10.1002/nme.2869.
[13]	ARROYO M, ORTIZ M. Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods [J]. International Journal for Numerical Methods in Engineering, 2006, 65(13): 2167–2202. DOI: 10.1002/nme.1534.
[14]	SANTAMBROGIO F. Introduction to optimal transport theory [EB/OL]. arXiv: 1009.3856. (2010-09-20)[2021-07-01]. https://doi.org/10.48550/arXiv.1009.3856.
[15]	LI B, STALZER M, ORTIZ M. A massively parallel implementation of the optimal transportation meshfree method for explicit solid dynamics [J]. International Journal for Numerical Methods in Engineering, 2014, 100(1): 40–61. DOI: 10.1002/nme.4710.
[16]	SCHMIDT B, FRATERNALI F, ORTIZ M. Eigenfracture: an eigendeformation approach to variational fracture [J]. Multiscale Modeling & Simulation, 2009, 7(3): 1237–1266. DOI: 10.1137/080712568.
[17]	LI B, PANDOLFI A, ORTIZ M. Material-point erosion simulation of dynamic fragmentation of metals [J]. Mechanics of Materials, 2015, 80: 288–297. DOI: 10.1016/j.mechmat.2014.03.008.
[18]	PANDOLFI A, LI B, ORTIZ M. Modeling fracture by material-point erosion [J]. International Journal of Fracture, 2013, 184(1/2): 3–16. DOI: 10.1007/s10704-012-9788-x.
[19]	樊江, 袁圆, 廖祜明, 等. 基于最优运输无网格法的Whipple屏超高速撞击数值模拟 [J]. 爆炸与冲击, 2019, 40(7): 074201. DOI: 10.11883/bzycj-2019-0241. FAN J, YUAN Y, LIAO H M, et al. Numerical simulation of Whipple shield hypervelocity impact based on optimal transportation meshfree method [J]. Explosion and Shock Waves, 2019, 40(7): 074201. DOI: 10.11883/bzycj-2019-0241.
[20]	STAINIER L. A variational approach to modeling coupled thermo-mechanical nonlinear dissipative behaviors [J]. Advances in Applied Mechanics, 2013, 46: 69–126. DOI: 10.1016/B978-0-12-396522-6.00002-5.
[21]	廖祜明. 整体拉格日无网格流固耦合计算方法 [D]. 北京: 北京航空航天大学, 2018: 33–34.
[22]	FAN J, LIAO H M, WANG H, et al. Local maximum-entropy based surrogate model and its application to structural reliability analysis [J]. Structural and Multidisciplinary Optimization, 2018, 57(1): 373–392. DOI: 10.1007/s00158-017-1760-y.
[23]	FAN Z Y, WANG H, HUANG Z D, et al. A Lagrangian meshfree mesoscale simulation of powder bed fusion additive manufacturing of metals [J]. International Journal for Numerical Methods in Engineering, 2021, 122(2): 483–514. DOI: 10.1002/nme.6546.
[24]	NAVAS P, YU R C, LI B, et al. Modeling the dynamic fracture in concrete: an eigensoftening meshfree approach [J]. International Journal of Impact Engineering, 2018, 113: 9–20. DOI: 10.1016/j.ijimpeng.2017.11.004.
[25]	PIEKUTOWSKI A J. A simple dynamic model for the formation of debris clouds [J]. International Journal of Impact Engineering, 1990, 10(1): 453–471. DOI: 10.1016/0734-743X(90)90079-B.
[26]	STEINBERG D J. Equation of state and strength properties of selected materials [M]. Livermore: Lawrence Livermore National Laboratory, 1996.

施引文献

期刊类型引用(1)
1. 廖祜明，龚自正，宋光明，樊宗岳，杨宏涛，黎波. OTM/HOTM极限力学仿真在小行星防御中的应用. 空间碎片研究. 2022(03): 5-18 . 百度学术
其他类型引用(0)

资源附件(0)

访问统计

图(11) / 表(2)

计量

文章访问数: 567
HTML全文浏览量: 333
PDF下载量: 111
被引次数: 1

1. OTM最优运输无网格方法
1.1 时间与空间离散
1.2 捕捉式多体动态接触算法
1.3 本征侵蚀裂纹扩展算法
2. 铜飞片超高速撞击铝合金靶板数的值模拟与验证
2.1 数值模拟模型
2.2 材料模型
2.3 结果分析
3. 结　论

超高速撞击下碎片云的OTM分析

doi: 10.11883/bzycj-2021-0275

作者简介: 廖祜明（1987－ ），男，博士，liaohuming@buaa.edu.cn

通讯作者: 樊 江（1973－ ），男，博士，副教授，fanjiang@buaa.edu.cn

计量

出版历程

OTM analysis of debris cloud under hypervelocity impact

1. OTM最优运输无网格方法

1.1 时间与空间离散

1.2 捕捉式多体动态接触算法

1.3 本征侵蚀裂纹扩展算法

2. 铜飞片超高速撞击铝合金靶板数的值模拟与验证

2.1 数值模拟模型

2.2 材料模型

2.3 结果分析

3. 结 论

期刊类型引用(1)

其他类型引用(0)

计量

出版历程

目录

1. OTM最优运输无网格方法

1.1 时间与空间离散

1.2 捕捉式多体动态接触算法

1.3 本征侵蚀裂纹扩展算法

2. 铜飞片超高速撞击铝合金靶板数的值模拟与验证

2.1 数值模拟模型

2.2 材料模型

2.3 结果分析

3. 结 论

作者简介:
廖祜明（1987－　），男，博士，liaohuming@buaa.edu.cn

通讯作者:
樊　江（1973－　），男，博士，副教授，fanjiang@buaa.edu.cn

3. 结　论

3. 结　论