双原子分子气体的定常隐式全流域多尺度算法.pdf

巨大的计算资源需求极大地阻碍了统一气体动理学格式的应用。采用宏观预估技术，基于 模型方程发展全流域适用的保守恒定常隐式算法，协同求解宏观方程和微观方程以加速收敛。在单元界面，通过模型方程特征差分解构造简单高效的多尺度数值通量，并结合非均匀非结构速度空间和速度空间自适应技术进一步降低计算需求、提升计算效率。超声速和高超声速平板绕流和圆球绕流的数值结果验证了算法的准确性与高效性。结果表明，算法能够准确求解二维和三维双原子气体多尺度流动问题，听语音聊科研与作者互动且相比于显式离散统一气体动理学格式可加速一个量级。关键词：多尺度算法；隐式算法；稀薄气体；双原子气体中图分类号：文献标志码：开放科学（资源服务）标识码（）：文章编号：（），（，；，；，）：，：；稀薄气体流动在工程领域普遍存在，如微机电系统 、临近空间高超声速飞行器和亚轨道、轨道航天器 等。如返回舱再入涉及跨流域流动问题，临近空间飞行器流场可同时存在连续与稀薄流动，其局部努森数 可相差几个量级 。对于涉及稀薄气体流动的多尺度流动问题，基于连续介质假设的纳维 斯托克斯（，）方程不再适用。稀薄流动模拟方法，如直接模拟蒙特卡罗 （，）方法，传统求解 方程及其模型方程的离散速度法 （，），其单元大小和时间步长受限于分子平均收稿日期：基金项目：国家自然科学基金资助项目（，）；国家“计划”资助项目（）作者简介：钟诚文（），男，重庆璧山人，教授，博士，博士生导师，：；卓丛山（通信作者），男，湖南桑植人，副教授，博士，硕士生导师，：第 期钟诚文，等：双原子分子气体的定常隐式全流域多尺度算法自由程和碰撞时间，难以高效准确计算连续流动。因此，混合算法 是最直接的计算连续 稀薄流动的算法，但其在计算区域划分、复杂分区界面处理、信息传递等方面依然存在挑战。近年来，徐昆等 在有限体积离散速度框架下发展了全流域适用的统一气体动理学格式（，），实现了从连续流域到自由分子流域统一求解。耦合求解分布函数控制方程和宏观量控制方程，其通量的计算真实地描述气体分子自由迁移与碰撞的物理过程，保证连续地捕捉任意 下的流场，并且其物理网格大小和时间步长不再受限于分子平均自由程和碰撞时间。随后，郭照立等 基于 的思想提出了离散统一气体动理学格式（，），其通过特征差分解构造的多尺度数值通量相对于 的当地积分解析解更简洁，计算效率相对较高。杨鲤铭等 发展的改进型 （，），袁瑞峰等 提出的多尺度 也具备全流域计算能力。尽管 和 提供了一个统一的框架求解全流域流动问题，但 模型方程需要在时间、物理空间和粒子速度空间中对分布函数演化方程进行离散求解，计算量较大。因此，发展加速算法是必由之路。由于宏观量和分布函数的耦合关系，碰撞项的全隐式离散求解存在困难。杨照彦等 直接采用 时间步的平衡态分布函数近似 时间步的平衡态分布，导致隐式算法收敛效率较低 。根据宏观量和分布函数的对应关系将平衡态分布函数进行线性化近似以保证控制方程的全隐式离散。但其映射关系为离散速度空间的超大型非稀疏矩阵，显著地增加了计算复杂性。毛枚良等 采用 和 的方法发展了半隐式 方法，等 基于毛枚良的方法开发了三维半隐式并行 求解器。半隐式 方法对碰撞项中的分布函数作隐式处理，平衡态分布函数作显式处理，收敛效率较低。考虑到 方法耦合求解气体分布函数和宏观量的思路，朱亚军等 通过求解隐式宏观控制方程来预估 时间步的宏观量，提出了全隐式离散的 方法。基于预估的宏观变量求解 时间步的预估平衡态分布函数，从而实现全隐式地求解气体分布函数控制方程（微观方程），宏观预估技术的使用极大提高了隐式 的效率。随后，袁瑞峰等 将隐式 推广到双原子分子气体流动模拟，张瑞等 将其进一步扩展到三维双原子分子气体高超声速流动模拟。杨鲤铭等 基于宏观预估技术发展了系列隐式 ，其宏观通量由 模型方程当地积分解近似求解，微观通量依然采用传统 的方式求解，其宏观微观不一致的问题导致收敛效率降低。苏微等 提出了数十步内收敛的广义合成迭代格式（，），可实现全流域快速收敛。宏观 求解器的权重设计不足导致了其结果存在不光滑区域，其宏观 求解器与微观常规迭代格式（，）的耦合方法是经验性的，在极端情况下存在一些不稳定现象。袁瑞峰等 和潘东鑫等 分别发展了单原子分子气体和双原子分子气体的多重预估隐式算法，对宏观预估的隐式 进行改进。方法采用具有二阶精度的 通量求解器来计算隐式宏观控制方程中的宏观量通量变化（隐式 采用 通量求解器计算），计算效率在隐式 基础上提升一个量级。此外，通量求解器在稀薄流求解时黏性和热流计算存在较大偏差，这导致方法在稀薄流求解时有一定的稳定性问题。离散速度类方法提高计算效率的另一个途径是在保证计算精度的同时尽可能减少离散分布函数（离 散 速 度 空 间 网 格）的 数 量。等 基于笛卡尔速度空间网格提出了速度空间网格自适应技术。等 提出了基于当地速度空间的积分方法，在不同的物理空间区域采用不同的离散速度空间，与传统在全物理空间采用一个全局速度空间相比可以有效提高计算效率。陈松泽等 在 中首次引入速度空间自适应技术，有效提高了 的计算效率。杨鲤铭等 采用参数化降阶建模方法优化离散速度空间，有效降低了离散速度空间数量。袁瑞峰等 、陈健锋等 引入非均匀非结构离散速度空间（，），根据流动基本特征灵活生成对应的离散速度空间网格，极大地降低了离散速度空间网格数量。陈健锋等 发 展 了 全 局 自 适 应 离 散 速 度 空 间（，）技术，离散速度空间数量进一步降低，并对自适应参数的选取进行了详细分析。空气最主要的成分是氧气和氮气，实际多尺度流动问题往往是双原子分子或多原子分子气国 防 科 技 大 学 学 报第 卷体。模型 是考虑了分子平动和转动自由度的动理学方程，能适用于较大范围的双原子分子气体流动问题研究。基于 模型方程发展了宏观预估的保守恒定常隐式算法，并采用简化多尺度数值通量 、非均匀非结构速度空间和速度空间自适应技术提高计算效率。双原子分子 模型方程对于双原子分子气体，仅考虑分子平动和转动自由度时，气体分子的数密度分布函数为（，），其中 是空间坐标，是 维空间中的分子平动速度，是分子速度在三维空间中的其他分量，是分子转动能，是时间。宏观量通过对分布函数在相空间 积分得到：（，）（）（，）（）()（，）（）（，）（）（，）（）（）（，）（）（，）（）（，）（）其中：是分子质量，是密度，是流体速度，是分子热运动速度，、和 分别是气体总能量、平动能量和转动能量，和 分别是平动自由度和转动自由度上的热流输运，是应力张量。分布函数在相空间的演化采用 模型方程 描述。（）和 分别是平动和转动平衡态分布函数：()()()（）（）()()()（）（）其中：是分子数密度，和 分别是平动温度和转动温度，是平衡态温度，是气体常数，和 是采用 正交多项式构造的热流调节项，目的是获得尽可能准确的热流松弛率。()（）()（）()（）()（）其中：是松弛时间，考虑到分子的碰撞频率主要受到分子平动速度的影响，松弛时间依赖于平动温度。（）（）（）其中：是黏性系数，是压力。是转动碰撞数，代表弹性碰撞与非弹性碰撞松弛速率之比，通常根据理论公式结合实验数据拟合，如 公式 。和 中的系数、沿用徐昆 和刘沙 在 方法中的取值，对氮气常取 ，。为了避免计算时对速度分量 和分子转动能离散，采用约化形式的模型方程后可减小计算量和存储量，引入三个约化分布函数 ：（，）（，）（，）（，）（，）（，）（）对应的 模型方程为：（）式中：，；，。平衡态分布函数为：（，）（）（，）（，）（，）（）（，）（，）（，）（）（，）（，）（）第 期钟诚文，等：双原子分子气体的定常隐式全流域多尺度算法其中：（，）（）()（）（，）（）()（）（，）（）（）()（）（，）（）（）()（）（，）（）（）（）（，）（）（）（）其中，是 平衡态分布函数。（）()()（）根据式（）（）可知，平衡态分布函数完全由宏观物理量确定。同时，宏观物理量 （，），平动温度 ，热流 和 ，应力张量 可表示为：（）（）（）（）（）（）（）（）（）（）多尺度定常隐式算法 算法总体框架时间方向采用欧拉向后离散，可以得到离散的隐式微观控制方程：（）（），珘 （）式中，为单元 的体积，代表单元 的相邻单元，（）代表 单元的所有相邻单元的集合，代表单元 和单元 之间的界面，为界面 上由单元 指向单元 的外法线单位向量，为界面 的面积，符号“”表示下一迭代步的预估量，表示数值时间步长。将式（）整理为 增量形式，得：珘 ()（）（）式中：，残差为，珘 （）（）残差中的预估平衡态分布函数，和松弛时间 珘 将由预估宏观量珦 计算，预估宏观量珦 可由隐式微观控制方程求解，将在后文阐述。界面 上的分布函数 的计算方法决定格式是否具有多尺度特性、能否用于从自由分子流到连续连的全流域模拟，其构造方式将在 节详细讨论。式（）中界面分布函数增量 采用简单的一阶迎风格式近似。，（），（）（）其中，（）为单元 的相邻单元中满足 条件的单元集合，而（）为满足 条件的相邻单元集合。最终式（）可整理为如下形式：珘 （）()（）（）采用点松弛对称高斯赛德尔（，）迭代法 求解式（）以获得分布函数增量 ，从而得到新迭代步上的分布函数 。考虑到式（）的 迭代需要在所有离散速度点上进行，计算非常耗时（相比于宏观控制方程），因此微观控制方程的 迭代步数设置为。国 防 科 技 大 学 学 报第 卷为高效求解式（），其残差中的预估平衡态分布函数，和松弛时间 珓 将通过求解隐式宏观控制方程获得。对分布函数控制方程在速度空间求矩，可以得到对应的宏观方程：（）其中：（，（）（）。同理，对式（）在物理空间采用有限体积离散，时间方向采用欧拉向后离散，可以得到离散的隐式宏观控制方程：（珦 ）（）珦 珘 （）将其整理为 增量形式，得：珦 （）珦 （）珦（）通过求解式（），即可获得在宏观量预估中所需计算的预估宏观量珦 。式（）中的宏观量通量增量珦 采用如式（）所示形式计算。珦 珚 （）其中：（）（）（）（）式中的 （）为欧拉通量，（）（）（）（）是欧拉通量雅可比矩阵的谱半径加上运动黏性项以保证计算稳定。（）其中：，槡。将式（）和式（）代入式（）可得：（）()珦 （）珦 （）式中宏观残差项为：（）（）（珦 ）（）（）式（）中的宏观界面通量 可以通过界面分布函数 在速度空间积分求解。（）（）式中界面分布函数 的计算将在下一节详细讨论。式（）中对于守恒量 、而言，源项 为零，关于守恒量 、的增量形成一组与转动能 增量无关的封闭方程组，因此可以先计算守恒量增量，再计算转动能增量。根据式（）和式（），守恒量 、的隐式宏观方程整理为：（）()珟 （）珟 （），（）（珦 ）（）（）式中，代表守恒量。对于转动能量，对源项 中转动能的部分进行线性化处理，可得转动能量增量的隐式控制方程：（）珘 ()珟 ，（）珟 ，珘 （），（）（珦 ）（）（）式中，代表转动能量。相比于隐式微观控制方程，隐式宏观控制方程计算量较小，并且预估的宏观量珦 越准确，隐式算法收敛越快，特别是连续流域。因此，隐式宏观控制方程迭代次数通常设置为 或 。多尺度数值通量在单元界面，采用袁瑞峰、张瑞等提出的简化多尺度数值通量 ，在当地物理时间步长 内沿着以单元界面 为终点的特征线积分分布函数演化方程（），碰撞算子的积分采用隐式处理，可以得到：第 期钟诚文，等：双原子分子气体的定常隐式全流域多尺度算法（，）（，）（，）（，）（，）（）其中：当地物理时间步长 （，），和 根据 条件确定。记 （，），（，），式（）整理为：（，）（，）（）式（）在求解单元界面的分布函数时耦合粒子迁移和碰撞过程。在自由分子流域，粒子以自由迁移为主，其松弛时间 远大于当地物理时间步长 ，式（）中代表粒子自由迁移的分布函数（，）将占据主导地位，从而精确描述自由分子流动状态下粒子自由输运的物理过程。在连续流域，松弛时间 远小于当地物理时间步长 。由于强烈的粒子碰撞，系统趋于平衡状态。因此，单元界面处的平衡态分布（，）将自动占据主导地位，从而准确捕捉连续流域的物理过程。式（）中，分布函数 （，）通过线性重构求解。（，）（）（?）?（）（?）?（）其 中：（?）是 非 结 构 网 格 下 常 用 的 限制器 ，梯度采用最小二乘法求解。式（）中，单元界面处的平衡态分布函数（，）根据界面宏观量求解。对式（）在速度空间积分，可得界面宏观量为：（，）（）（，）（）（，）（，）（，）（）（）（，）（）（）()（，）（，）（）（）（）()（，）（）至此，式（）中分布函数 （，）和平衡态分布函数（，）均已求解，则界面分布函数 重构完成。综上，隐式算法的计算步骤可以归纳为：）初始化流场，分布函数初始化为平衡态分布；）计算界面分布函数（，），求解微观通量和对应的宏观通量；）计算隐式宏观方程的残差，迭代求解隐式宏观方程得到预估宏观量珦 ，计算预估平衡态分布函数，；）计算隐式微观方程的残差，迭代求解隐式微观方程，更新分布函数 ；）对分布函数 在速度空间积分得到宏观量 ，数 值 积 分 误 差 补 偿 方 法 参 见 文献 。离散速度空间宏观量通过对连续速度空间中的分布函数积分得到，实际计算时，需要在离散的速度空间进行数值积分。传统的笛卡尔速度空间在计算稀薄高超声速流动时网格量较大，本文将采用非结构速度空间和自适应速度空间技术减少速度空间网格量。非结构速度空间的数值积分（以计算密度为例）表示为：（）其中：下标 表示离散速度编号；表示积分权重，即速度空间网格单元面积（二维）或体积（三维）。自适应速度空间在计算中根据流场特征自动优化速度空间网格。本文提出的全局自适应速度空间无须在物理空间进行插值求解。另一方面，在自适应得到的速度网格上重构分布函数时，通过“亲缘”关系实现新速度空间上的分布函数的重构，如图 所示，当新速度空间上国 防 科 技 大 学 学 报第 卷的网格是由旧速度空间网格分裂得到时，新速度网格上此处的分布函数等于旧速度网格的分布函数。当新速度空间上的网格是由旧速度空间网格合并得到时，新速度网格上此处的分布函数等于旧速度网格的分布函数的平均值。为了保证算法的守恒性，可利用新旧速度网格上的平衡态偏差来修正通过“亲缘”关系直接得到的分布函数。，珚，（），（珡）（）其中：珚，是通过“亲缘”关系得到的分布函数，珡是对珚，在新速度空间上积分得到宏观量，是旧速度空间上的分布函数积分得到的宏观量。图 速度空间网格树形结构 数值验证 稀薄超声速平板绕流稀薄超声速平板绕流算例源自 和 的风洞实验 。风洞实验给出了受到高温热非平衡效应影响的转动温度数据。的实验条件是：实验气体为氮气，喷嘴出口马赫数为 ，喷嘴出口总压为 ，总温为 ，平板温度恒定为 。计算采用变径硬球模型，转动碰撞数设为定值 。计算选取的计算域和非结构物理空间网格如图 所示，网格单元数量为 。非结构速度空间网格如图 所示，网格单元数量为 ，速度空间网格在对应来流速度和壁面速度附近的区域进行了加密。图 给出了自适应之后的速度空间网格，相比于非结构速度空间，网格数量减少 左右。图 稀薄超声速平板绕流计算所用的物理空间网格（单元）（）图 稀薄超声速平板绕流计算所用的非结构离散速度空间网格（单元）（）图 稀薄超声速平板自适应离散速度空间网格（单元）（）第 期钟诚文，等：双原子分子气体的定常隐式全流域多尺度算法图 是 等于 、处竖直截线上的温度分布以及上平板 横截线处的温度分布。目前的模拟结果和采用 计算的结果吻合较好，并且转动温度和实验结果吻合很好。图 给出了壁面压力、热流影响的对比结果（横坐标 是贴体坐标，起点是平板下表面尾部的点），当前方法结合非结构速度空间和自适应速度空间网格的计算结果基本一致，和 的计算结果吻合很好。表 对比了隐式算法相对于显式 的加速比（计算所用 型号为 （）（），速度空间采用 核 并行计算）。（）（）（）（）图 稀薄超声速平板绕流竖直截线和横截线的温度分布结果与 以及实验结果的对比 （）压力（）（）热流（）图 稀薄超声速平板绕流壁面压力和热流分布 当前隐式算法相比于显式 可加速 倍以上，采用自适应速度空间后离散速度空间数量大国 防 科 技 大 学 学 报第 卷表 隐式算法相比于显式 的加速比（核）（）计算方法计算时间 计算步数加速比 大降低，加速比可达到 倍以上。表 给出了采用不同核数并行计算时自适应速度空间相比非结构速度空间的加速比。随着单个核心离散速度空间数量逐渐增加，自适应速度空间相比非结构速度空间的加速比逐渐接近速度空间总数量之比。速度 空 间 的 消 息 传 递 接 口（，）并行计算单个核心的离散速度数量不能太少，否则并行效率明显下降。表 非结构速度空间和自适应速度空间不同并行核数下的加速比 并行核数 计算时间 计算步数计算时间 计算步数加速比 绕钝头前缘平板的高超声速流动绕钝头前缘平板的高超声速流动算例取自 等的实验 。数值结果将与实验结果 和 结果比较。实验条件为：实验气体为氮气，来流密度为 ，来流速度为 ，来流温度为 ，为 。平板板长为 、宽为 、板的厚度为 ，壁面温度恒定为 ，迎角 为 和 ，以板长为参考长度的 为 。采用变径硬球模型，转动碰撞数设为定值 。计算选取的计算域和非结构物理空间网格如图 所示，网格单元数量为 。非结构速度空间网格如图 所示（，对于 迎角只需将其按照来流角度旋转 即可），网格单元数量为 ，速度空间网格在对应于来流速度和壁面速度附近的区域进行了加密。图 给出了自适应之后的速度空间网格，相比于非结构速度空间，网格数量可减少 左右。图 高超声速平板绕流计算所用的物理空间网格（单元）（）图 高超声速平板绕流计算所用的非结构离散速度空间网格（单元）（）（）（单元）（）（）第 期钟诚文，等：双原子分子气体的定常隐式全流域多尺度算法（）（单元）（）（）图 高超声速平板绕流计算的自适应离散速度空间网格 图 、图 、图 分别给出了平板壁面压力系数、热流系数和剪应力系数分布。采用非结构速度空间和自适应速度空间的计算结果基本一致。壁面压力系数、热流系数和 计算结果 （等 采 用 图 平板壁面压力系数分布 开发的 代码计算）吻合较好。与 等的实验结果 对比，数值结果和实验结果的差别主要来源于边界条件处理，当前方法和 均采用完全热适应的漫反射边界条件，与实验条件存在明显差异。壁面剪应力系数与 等 采用 计算的结果吻合较好。表 和表 对比了隐式算法相对于显式 的加速比（计算所用 型号为 （）（），速度空间采用 核 并行计算），当前隐式算法相比于显式 可加速 倍；若使用自适应速度空间，加速比可达到 倍以上。表 和表 给出了采用不同核数计算时自适应速度空间相比非结构速度空间的加速比。同样地，随着单个核心离散速度空间数量逐渐增加，自适应速度空间相比非结构速度空间的加速比逐渐接近速度空间总数量之比。图 平板壁面热流系数分布 图 平板壁面剪应力系数分布 表 隐式算法相比于显式 算法的加速比（，核）（，）计算方法计算时间 计算步数加速比 国 防 科 技 大 学 学 报第 卷表 隐式算法相比于显式 算法的加速比（，核）（，）计算方法计算时间 计算步数加速比 表 非结构速度空间和自适应速度空间不同并行核数下的加速比（）（）并行核数 计算时间 计算步数计算时间 计算步数加速比 表 非结构速度空间和自适应速度空间不同并行核数下的加速比（）（）并行核数 计算时间 计算步数计算时间 计算步数加速比 超声速和高超声速圆球绕流 进行了非连续流区圆球阻力的测量实验，实验气体为空气。本算例对不同 下的超音速（，）和高超音速（，）圆球绕流进行模拟，计算采用单组分氮气。球体直径 ，自由来流总温为 ，壁面温度分别设置为 （）和 （）。黏性系数采用 公式 计算，转动碰撞数设置为 。图 展示了物理空间网格的截面图，网格量为 个六面体单元，壁面第一层网格高度为 。计算分别采用非结构速度空间网格和自适应速度空间网格，图 展示了 条件的非结构速度空间网格（个单元）和自适应速度空间网格（个单元）。图 高超声速圆球绕流计算所用的物理空间网格 （）非结构速度空间网格（）表 给出了隐式算法结合非结构速度空间和自适应速度空间计算的阻力系数和实验值的对比结果，相对误差均小于 。考虑到实验数据的均方根误差为 ，当前计算结果吻合很好。图 、图 、图 分别给出了圆球（）表面压力系数、热流系数和剪应力系数分布，自适应速度空间网格相比于非结构速度空间网格在驻点热流的计算上更加接近 的结果，压力系数和剪应力系数方面与 的结果基本吻合。图 、图 、图 分别给出了圆球（）表面压力系数、热流系数和剪应力系数分布，压力系数和剪应力系数方面与 的结果第 期钟诚文，等：双原子分子气体的定常隐式全流域多尺度算法（）自适应速度空间网格（）图 高超声速圆球计算所用的非结构离散速度空间网格（单元）和自适应速度空间网格（单元）（）（）表 圆球绕流阻力系数对比 实验 计算值 相对误差 计算值相对误差 图 超声速圆球绕流壁面压力系数分布（）（）图 超声速圆球绕流壁面热流系数分布（）（）图 超声速圆球绕流壁面剪应力系数分布（）（）图 高超声速圆球绕流壁面压力系数分布（）（）吻合，自适应速度空间网格计算的热流系数比非结构速度空间计算的结果更加接近 的结果。自适应速度空间能更好地捕捉流动特征。表 和表 给出了当前隐式算法相比于显式 算法的加速比（计算所用 为 （）（），其中非结构速度空间采用 核 并行计算，自适应速度空间采国 防 科 技 大 学 学 报第 卷图 高超声速圆球绕流壁面热流系数分布（）（）图 高超声速圆球绕流壁面剪应力系数分布（）（）表 隐式算法相比于显式 算法的加速比（）（）计算方法计算时间 计算步数加速比 表 隐式算法相比于显式 算法的加速比（）（）计算方法计算时间 计算步数加速比 用 核 并行计算，两者单核心的离散速度空间数量基本一致（保证两者单核心的离散速度空间数量一致，自适应速度空间应该采用大概 核并行计算，但节点内存受限，因此采用 核并行计算）。隐式算法相比于 可加速 倍左右，结合自适应速度空间后可加速 倍左右，可显著提高计算效率。结论基于 模型方程发展了适用于全流域双原子分子气体流动模拟的保守恒定常隐式算法。全隐式离散、简化多尺度数值通量、非结构和自适应速度空间的应用确保算法相比于 加速一到两个量级以上。超声速和高超声速平板绕流和圆球绕流算例的数值结果和 以及实验结果吻合较好。圆球绕流算例的气动力与实验结果相对误差小于 ，壁面压力分布和热流分布与 结果一致，充分验证了本算法的准确性与高效性，方法初步具备成本可控的三维全流域复杂流动模拟能力。参考文献（），：，：，?，：江定武基于模型方程解析解的气体动理学算法研究 绵阳：中国空气动力研究与发展中心，：，（）：，（）：，（）：第 期钟诚文，等：双原子分子气体的定常隐式全流域多尺度算法 ，（）：，：，（）：，（）：，（），（）：，：，（）：，（）：，（）：毛枚良，江定武，李锦，等气体动理学统一算法的隐式方法研究 力学学报，（）：，（）：（），（）：，：，（）：，：，：，：，（）：，：，（）：，？，：，：，：，：，（）：，（）：，：，（）：，：，（）：，国 防 科 技 大 学 学 报第 卷 ，（）：，（）：，：，（）：，（）：，（）：，（）：?，：，：，：，