1.1 爆炸毁伤数值仿真基础

1.1.1 数值仿真技术概况

数值仿真方法主要包括有限差分方法、有限元方法、有限体积法、SPH方法等，如图1-1所示。

1．有限差分方法

有限差分方法是一种直接将微分问题变为代数问题的近似数值解法。这种方法首先将求解域划分为网格，然后通过泰勒级数展开等方式，将控制方程中的导数用网格节点上函数值的差商代替，进行离散操作，从而建立以网格节点上的值为未知数的代数方程组，通过解算这些方程组获得问题的近似解。当采用较多的网格节点时，有限差分方法所得近似解的精度可以得到保证。该方法具有数学概念直观、表达简单等优点，是发展较早且比较成熟的数值方法。

2．有限元方法

有限元方法的基础是变分原理和加权余量法，其求解思想是把计算域离散为一组有限个，且按一定方式相互连接在一起的单元的组合体，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达值，借助变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。由于单元能按不同的连接方式进行组合，且单元本身又可以有不同形状，因此可以模型化几何形状复杂的求解域。

3．有限体积法

有限体积法又称为控制体积法、有限容积法。其基本思路是：将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解的微分方程对每一个控制体积积分，便得出一组离散方程。其中的未知数是网格点上的因变量的数值。为了求出控制体积的积分，必须假定值在网格点之间的变化规律。有限体积法的基本思路易于理解，适于解决复杂的工程问题，且具有良好的网格适应性。

4．SPH方法

SPH（Smoothed Particle Hydrodynamics）是光滑粒子流体动力学方法的缩写，是近年来兴起并逐渐得到广泛应用的一种数值模拟方法，属于无网格方法。该方法的基本思想是将连续的流体（或固体）用相互作用的质点组来描述，各个物质点上承载各种物理量，包括质量、速度等，通过求解质点组的动力学方程和跟踪每个质点的运动轨道，求得整个系统的力学行为。由于SPH方法中质点之间不存在网格关系，避免了极度大变形时因网格扭曲造成的精度下降，因此，在处理结构大变形、冲击破坏等方面具有很大优势，特别适合在弹药毁伤仿真领域的应用。

5．数值仿真实例

数值仿真方法是求取复杂微分方程近似解的非常有效的工具，是现代数字化科技的一种重要基础性原理。在科学研究中，数值仿真方法可称为探究物质客观规律的先进手段。将它用于工程技术中，可成为工程设计和分析的可靠工具。严格来说，数值仿真分析必须包含三个方面：①数值仿真方法的基本数学原理；②基于原理所形成的实用软件；③仿真时的计算机硬件。本书的重点是通过一些典型的弹药毁伤效应实例，运用ANSYS分析平台来系统地阐述数值仿真分析的基本原理，展示具体应用数值仿真方法的建模过程。

基于功能完善的数值仿真分析软件和高性能的计算机硬件对设计的结构进行详细的力学分析，以获得尽可能真实的结构受力信息，就可以在设计阶段对可能出现的各种问题进行安全评判和设计参数修改。据有关资料报道，一个新产品的问题有60%以上可以在设计阶段消除，甚至有些结构的施工过程中也需要进行精细的设计，要做到这一点，就需要数值仿真这样的分析手段。

以下是数值仿真方法在工程实例中的应用：

空客A350后机身第19框的设计与有限元分析过程如图1-2所示。

北京奥运场馆鸟巢的实物和有限元模型对比如图1-3所示。

图1-4所示为军用车辆的底部防地雷模块的设计过程，以及采用数值仿真分析的情况。

图1-5所示为高能炸药爆炸时的试验与数值仿真结果，可通过两种结果进行分析比较，研究炸药爆炸的整个过程。

1.1.2 数值模拟基本过程

针对具有任意复杂几何形状的变形体，完整获取在复杂外力作用下其内部的准确力学信息，即求取该变形体的三类力学信息（位移、应变、应力）。在准确进行力学分析的基础上，设计师就可以对所设计对象进行强度（strength）、刚度（stiffness）等方面的评判，以便对不合理的设计参数进行修改，以得到较优化的设计方案；然后，再次进行方案修改后的数值仿真分析，以获得最后的力学评判和校核，确定出最终的设计方案。以有限元分析为例，其工作流程如图1-6所示。

有限元方法是基于“离散逼近（discretized approximation）”的基本策略，可以采用较多数量的简单函数的组合来“近似”代替非常复杂的原函数。因此，采用有限元方法可以针对具有任意复杂几何形状的结构进行分析，并能够得到准确的结果。一个复杂的函数，可以通过一系列的基底函数（based function）的组合来“近似”，也就是函数逼近，其中有两种典型的方法：①基于全域的展开（如采用傅里叶级数展开）; ②基于子域（sub-domain）的分段函数（pieces function）组合（如采用分段线性函数的连接）。

采用有限元方法分析问题的基本步骤如下所示：

（1）建立积分方程：根据变分原理或方程余量与权函数正交化原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。

（2）区域单元剖分：根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点的位置坐标，同时，还需要列出自然边界和本质边界的节点序号和相应的边界值。

（3）确定单元基函数：根据单元中节点数目及对近似解精度的要求，选择满足一定插值条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的，由于各单元具有规则的几何形状，在选取基函数时，可遵循一定的法则。

（4）单元分析：将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再将近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数（即单元中各节点的参数值）的代数方程组，称为单元有限元方程。

（5）总体合成：在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则进行累加，形成总体有限元方程。

（6）边界条件的处理：一般边界条件有三种形式，即本质边界条件（狄里克雷边界条件）、自然边界条件（黎曼边界条件）、混合边界条件（柯西边界条件）。对于自然边界条件，一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件，需按一定法则对总体有限元方程进行修正满足。

（7）解有限元方程：根据边界条件修正的总体有限元方程组，是含所有待定未知量的封闭方程组，采用适当的数值计算方法求解，可求得各节点的函数值。

有限元方程是一个线性代数方程组，一般有两大类解法：一是直接解法，二是迭代法。直接解法有高斯消元法和三角分解法，如果方程规模比较大，可用分块解法和波前解法。迭代法有雅可比迭代法、高斯-赛德尔迭代法和超松弛迭代法等。

通过选用合适的求解法求解经过位移边界条件处理的公式后，得到整体节点位移列阵，然后根据单元节点位移由几何矩阵和应力矩阵得到单元节点的应变和应力。对于非节点处的位移，通过形函数插值得到，再由几何矩阵和应力矩阵求得相应的应变和应力。

应变要通过位移求导得到，精度一般要比位移差一些，尤其对于一次单元，应变和应力在整个单元内是常数，应变和应力的误差会比较大，特别是当单元数比较少时，误差更大，因此，对于应力和应变，要进行平均化处理：

（1）绕节点平均法，即依次把围绕节点的所有单元的应力加起来平均，以此平均应力作为该节点的应力。

（2）二单元平均法，即把相邻的两单元的应力加以平均，并以此作为公共边的节点处的应力。

整理并对有限元法计算结果进行后处理，一是要得到结构中关键位置力学量的数值（如最大位移、最大主应力和主应变、等效应力等）；二是得到整个结构的力学量的分布（根据计算结果直接绘制位移分布图、应力分布图等）；三是后处理要得到输入量和输出量之间的响应关系。

1.1.3 网格生成技术

网格生成技术是指对不规则物理区域进行离散，以生成规则计算区域网格的方法。网格是CFD模型的几何表达形式，也是模拟与分析的载体。对于复杂的CFD问题，网格的生成极为耗时，并且极易出错，生成网格所需的时间常常大于实际CFD计算时间。因此，有必要对网格生成技术进行研究。

对于连续介质系统，例如飞行器周围的气体，集中在障碍物上的压力，回路中的电磁场，或者是化学反应器中的浓缩物，都可以用偏微分方程来进行描述的。为了对这些系统进行模拟，需要基于一定数量的时间、空间意义的点对连续性方程进行离散化，并且在这些点上对各种物理量进行计算。离散的方法通常有下列三种：有限差分法、有限体积法、有限元法，都是使用相邻的点来计算所需要的点。

一般来说，通过连接点的方式的不同，可以把网格类型分为两种：结构化网格和非结构化网格。结构化网格是正交的处理点的连线，即意味着每个点都具有相同数目的邻点；而非结构化网格是不规则的连接，每个点周围的点的数目是不同的。图1-7给出了两种网格的例子。

在一些情况中，也有部分网格是结构的，部分网格是非结构的，例如，在黏性流体中，边界一般使用结构网格，其他部分使用非结构网格。

1．离散法类型

离散的主要方法有有限差分法、有限体积法和有限元法，为了说明这些方法，首先来考虑连续性方程。

式中，ρ是密度；S是源项；U是速度，表示各个方向上的质量流的速度。有限差分法是用下面的办法来达到对所需要的点的模拟的。例如，对正交网格，矩形在横轴的长度是h：

有限差分法适用于规则的网格，但对于不规则的网格也可以运用，也可以在特殊的坐标系中对正交网格使用有限差分法（例如在球形极坐标中）。

在有限体积法中，物理空间被分成很多小的体积V，对每一个小的体积运用偏微分方程进行积分：

然后用每个小体积中的每个所求量的平均值来代替所要求的值，用相邻体积中变量的函数来表示流过每个体积表面的流量。运用有限体积法进行离散，适用于结构或者非结构网格。在非结构网格中，每个表面上的流通量依然可以用相邻的变量来进行很好的定义。

有限元方法也是把空间分为很多小的体积，相当于很多小的单元，然后在每个单元里，变量和流通量都用势函数来表示，计算的变量都是这些势函数中的系数。有限元方法在结构网格运用中没有明显的优势，但在非结构网格中被普遍使用。

2．结构化网格

从严格意义上讲，结构化网格是指网格区域内所有的内部点都具有相同的毗邻单元。网格系统中节点排列有序，每个节点与邻点的关系固定不变。结构化网格具有以下优点：

（1）它可以很容易地实现区域的边界拟合，适于流体和表面应力集中等方面的计算。

（2）网格生成的速度快。

（3）网格生成的质量好。

（4）数据结构简单。

（5）对曲面或空间的拟合大多数采用参数化或样条插值的方法得到，区域光滑，与实际的模型更容易接近。

结构化网格最典型的缺点是适用的范围比较窄。尤其随着近几年的计算机和数值方法的快速发展，人们对求解区域的复杂性的要求越来越高，在这种情况下，结构化网格生成技术就显得力不从心了。

结构化网格生成技术主要有：正交曲线坐标系中的常规网格生成法、贴体坐标法和对角直角坐标法。结构化网格生成法结构如图1-8所示。

3．非结构化网格

同结构化网格的定义相对应，非结构化网格是指网格区域内的点不具有相同的毗邻单元。即在这种网格系统中节点的编号命名并无一定规则，甚至是完全随意的，并且每一个节点的邻点个数也不是固定不变的。从定义上可以看出，结构化网格和非结构化网格有相互重叠的部分，即非结构化网格中可能会包含结构化网格的部分。

非结构化网格技术从20世纪60年代开始得到了发展，主要是弥补结构化网格不能够解决任意形状和任意连通区域的网格剖分的缺欠。到90年代时，非结构化网格的文献达到了它的高峰时期。由于非结构化网格的生成技术比较复杂，随着人们对求解区域的复杂性的不断提高，对非结构化网格生成技术的要求也越来越高。从现有的文献情况来看，非结构化网格生成技术中只有平面三角形的自动生成技术比较成熟（边界的恢复问题仍然是一个难题，现在正在广泛讨论），平面四边形网格的生成技术正在走向成熟。而空间任意曲面的三角形、四边形网格的生成技术，三维任意几何形状实体的四面体网格和六面体网格的生成技术还远远没有达到成熟，需要解决的问题还非常多。主要的困难是从二维到三维以后，待剖分网格的空间区非常复杂，除四面体单元以外，很难生成同一种类型的网格。需要各种网格形式之间的过渡，如金字塔形、五面体形等。

对于非结构化网格技术，可以根据应用的领域，分为应用于差分法的网格生成技术（grid generation technology）和应用于有限元方法中的网格生成技术（mesh generation technology）。应用于差分计算领域的网格，除了要满足区域的几何形状要求以外，还要满足某些特殊的性质（如垂直正交、与流线平行正交等），因而从技术实现上来说就更困难一些。基于有限元方法的网格生成技术相对非常自由，对生成的网格只要满足一些形状上的要求就可以了。

非结构化网格生成技术还可以从生成网格的方法来区分，主要有以下一些生成方法：

对平面三角形网格生成方法，比较成熟的是基于Delaunay准则的一类网格剖分方法（如Bowyer-Watson Algorithm和Watson's Algorithm）和波前法（Advancing Front Triangulation）的网格生成方法。另外，还有一种基于梯度网格尺寸的三角形网格生成方法，这一方法现在还在发展当中。基于Delaunay准则的网格生成方法的优点是速度快，网格的尺寸比较容易控制。缺点是对边界的恢复比较困难，很可能造成网格生成的失败，对这个问题的解决方法现在正在研究中。波前法的优点是对区域边界拟合的比较好，所以，在流体力学等对区域边界要求比较高的情况下，常常采用这种方法。它的缺点是对区域内部的网格生成的质量比较差，生成的速度比较慢。

曲面三角形网格生成方法主要有两种：一种是直接在曲面上生成曲面三角形网格；另外一种是采用结构化和非结构化网格技术偶合的方法，即在平面上生成三角形网格以后再投影到空间的曲面上，这种方法会造成曲面三角形网格的扭曲和局部拉长，因此，在平面上必须采用一定的修正技术来保证生成的曲面网格的质量。

平面四边形网格的生成有两类主要方法：一类是间接法，即在区域内部先生成三角形网格，然后分别将两个相邻的三角形合并成一个四边形。生成的四边形的内角很难保证接近直角。所以，再采用一些相应的修正方法（如Smooth）加以修正。这种方法的优点是首先就得到了区域内整体的网格尺寸的信息，对四边形网格尺寸梯度的控制一直是四边形网格生成技术的难点。缺点是生成的网格质量相对比较差，需要多次修正，同时，需要首先生成三角形网格，生成的速度也比较慢，程序的工作量大。

另一类是直接法，二维的情况称为铺砖法（paving method）。采用从区域的边界到区域的内部逐层剖分的方法。这种方法到现在已经逐渐替代间接法而成为四边形网格的主要生成方法。它的优点是生成的四边形的网格质量好，对区域边界的拟合比较好，最适合流体力学的计算。缺点是生成的速度慢，程序设计复杂。空间的四边形网格生成方法到现在还是主要采用结构化与非结构化网格相结合的网格生成方法。

三维实体的四面体和六面体网格生成方法现在还远远没有达到成熟。部分四面体网格生成器虽然已经达到了使用的阶段，但是对任意几何体的剖分仍然没有解决，现在的解决方法就是采用分区处理的办法，将复杂的几何区域划分为若干个简单的几何区域，然后分别剖分再合成，对凹区的处理更是如此。

六面体的网格生成技术主要采用的是间接方法，即以四面体网格剖分作为基础，然后生成六面体。这种方法生成的速度比较快，但是生成的网格很难达到完全的六面体，会剩下部分四面体，四面体和六面体之间需要金字塔形的网格来连接。现在还没有看到比较成熟的直接生成六面体的网格生成方法。

图1-9所示是网格生成技术在具体实例中的应用。