国外从20世纪60年代起就开始用有限元法进行车身强度和刚度的计算。1970年美国国家航空航天局的NASTRAN 结构分析程序被引入汽车结构分析中，对车身结构进行了静强度有限元分析，使车身减小了质量，为车身轻量化开了先河。1974年Nagy L.I.用子结构法对车身进行了静态分析。1995年德国保时捷（Porsche）公司的美国分公司完成了ULSAB 钢质车身轻量化设计的项目，改善了占车身本体净质量25%的零部件，使轿车车身的总质量由271kg降至205kg。目前，国外各大汽车公司利用有限元软件进行车身结构静态分析、模态分析的技术已日臻成熟。

美国福特汽车公司早在20世纪70年代就已使用有限元软件，用板梁单元进行车身的静态分析，找出高应力区，并改进应力分布。日本五十菱汽车公司在80年代末已将CAE 应用到车身设计的各个阶段，从最初设计阶段的粗略模型到设计中后期的细化模型，分析的范围包括强度、刚度、振动、疲劳、碰撞及形状和质量的优化，进入90年代，有限元分析得到了更为广泛的应用。美国通用汽车公司在通用有限元程序的基础上自主开发了后处理程序，将发动机和道路激励载荷集成到数据库中，进行汽车对发动机和道路激励的响应分析和改进，极大地简化了分析过程。日本尼桑汽车公司利用有限元分析仿真来驱动整个设计过程，减少了设计时间，在分析中使用的模型已经包括悬架、发动机、轮胎和转向机构，使花费、质量和NVH（Noise、Vibration、Harshness）性能得到优化。美国福特汽车公司也利用CAE 在新车开发中提高其NVH 性能，并取得良好效果。

国外使用有限元对武器系统进行分析研究的公开文献不是太多，Hopkins D.A.等采用三维梁单元建立了火炮动力学模型，并进行了动态响应分析；Wilkerson S.A.和Kaste R.D.等一直致力于有限元技术在火炮动力学建模分析中的应用研究，并对M256 120mm 坦克炮的研究给出了梁单元和三维实体单元模型与试验测试结果的对比分析；Mcgrath S.V.等用有限元程序对炮管在横向载荷作用下的动力响应进行了分析，取得了一定的成果。

国内用有限元法分析车身结构始于20世纪70年代后期。浙江省交通科学研究所应用有限元方法在西门子7739 计算机上对大客车车身进行了强度计算，长春汽车研究所、吉林工业大学等单位也应用有限元方法对客车进行了静态分析。进入80年代，在汽车结构分析中，有限元分析方法逐步开始推广应用。冯国胜等对客车进行了静强度分析，表明模拟计算可以提供足够准确的车身刚度特性以及整车结构应力分布的大致规律，并能够对高应力区进行改进计算。吉林工业大学的王裴和沈阳轿车厂的刘听等探讨了车身车架模型建立和计算中的一些技术问题，并对SY-622B 客车进行了有限元计算分析。

国内从20世纪80年代后期到90年代开展了有限元法应用到武器设计的一些研究。当时的有限元模型以梁单元为主，单元数目少，模型比较简单，主要仿真振动特性，没有分析动应力、应变和热-应力耦合场；随着有限元理论和大型软件不断的完善，有限元法已经具备解决大规模非线性问题、随机问题、多物理场耦合问题的能力，用它来分析装甲车辆及各种武器在各种工况过程中的动力学特性成为必然。

现代装甲车辆车体结构的分析方法包括数值模拟和实验分析法。现代数值模拟分析方法主要是有限元分析方法。这种方法是依据实物等技术资料建立车体及行驶装置等的有限元分析模型，应用通用有限元软件计算和分析整车结构的静态、动态等特性指标，甚至进行优化设计分析和试验仿真研究。现代数值模拟分析方法可以在整车结构开发初期就预测和优化整车的静态和动态特性指标。从而在产品生产或试制之前就尽可能避免相关设计缺陷，提高产品成功率，缩短产品开发周期。有限元分析方法的精度取决于模型、工况、分析方法、对有限元基本知识的理解、软件应用和工程经验。现代实验方法主要是电测法，应用传感器、测量和分析仪器，对车体实物零部件或小部件模型进行支撑、加载和测试。试验分析方法虽然信息没有数值模拟分析方法充分，但是它可以为数值分析提供对比和模拟验证。

对于坦克装甲车辆工程而言，有限元结构分析的主要研究目的有以下三点：

（1）在设计阶段，检查车辆零部件结构强度是否存在关键性和一般性的潜在隐患，若存在，则与设计人员一起，根据分析结果进行结构改进，防患于未然。

（2）对已发生的车辆零部件结构损坏事故，进行有限元结构分析，找出原因，查明真相，提出结构改进方案和措施。

（3）对车辆零部件或整车进行减重与结构优化。与传统的车辆设计研制方法相比可以大幅缩短研制周期，节省大量的人力和物力。

结构分析与疲劳仿真的通用流程（结合不同软件）描述，分析步骤描述等如下。

有限元法的基本思想是将连续的求解区域离散为一组有限个数并且按一定方式相互连接在一起的单元的组合体。由于单元能按照不同的连接方式进行组合，且单元本身又可以有不同形状，因此可以模型化几何形状复杂的求解域。它利用每一个单元内假设的近似函数来分片地表示全求解域上待求的未知函数。从而将无限自由度问题离散为有限自由度问题，得到近似解。

有限元法的基本步骤是：①问题及其定义域的定义；②定义域的离散化；③各种状态变量的确定；④问题的公式表示；⑤建立坐标系；⑥构造单元的近似函数；⑦求单元矩阵和方程；⑧坐标变换；⑨单元方程的组装；⑩边界条件的引入；○1最终联立方程组求解；○12结果的解释。

有限元问题可按照静力或动力来区分，也可以按照线性或非线性来区分。整车结构通过障碍路面问题为典型的非线性动力学问题。动力学有限元的直接积分有隐式算法和显式算法两种。隐式算法的解是稳定的，但隐式算法中必须形成整体刚度矩阵，从而需要求解大型的联立非线性方程组，因此计算量大且要求计算机具有很大的存储空间。特别对大型结构计算问题，由于节点多、信息量大，因此计算时间很长，同时对计算机的内存要求也高。此外，接触碰撞过程中的高度非线性问题使得迭代计算的收敛性很难保证，这也在一定程度上限制了隐式算法有限元法的应用。显式动力学有限元方法利用中心差分离散时间域，无须构造刚度矩阵即可求解节点的运动方程，有效地回避了因非线性问题引起的收敛性问题。由于算法简单，能够实现高度向量化和并行化，在整车结构过障仿真等领域中得到比较普遍的应用。显式算法的缺点在于解的稳定性是有条件的，即积分时间步长很小，必须满足Courant 准则。下面具体介绍显式算法与隐式算法的区别。

考虑一个非线性结构的动力响应问题，矩阵动力平衡方程的形式是：

式中，u（t）， 分别为位移、速度、加速度的n 维与时间有关的矢量；M 为结构的质量矩阵，对称和正定； 为结构中的内部抵抗力，可以依赖于位移和速度；q（u，t）为外力，一般情况随着时间变化，但是也可以依赖于位移；u0 ，为位移和速度的初值。

在有限元求解时，假设在时刻tn+1 的位移和速度用线性差分公式近似，则运用于式（2-1）的时间积分算子总可写成如下形式：

式中，h 是时间步长；α ，β 是差分公式的特定系数。

1.显式算法

如果两个系数α = β = 0 ，则方程（2-2）称为显式，可从时刻tn ，tn+1 …的已知量推算位移un+1 和速度。

这些推算变成

并能根据以前时刻的解直接计算出来，然后解平衡方程以确定加速度。

显然，显式积分本质上是一个非迭代技术，它涉及两个主要计算：计算外力和恢复力之间的平衡：

解一个相应于质量矩阵的线性系统，如果质量矩阵是对角形式，只涉及平常的运算。

显式方法在时间行进过程中所允许的步长受到公式稳定区的强烈限制，对自由度之间的代数约束也不易处理。

2.隐式算法

在隐式情况下，两个系数α 和β 都不等于零，加速度和速度变成un+1 的函数，

对显式求解，准确度和稳定性要求支配步长的选择，检查数值不稳定性的一个方法是在整个动力响应过程中监测能量平衡。能量平衡的丧失是由于不稳定计算产生伪能量，这对非线性情况特别有效。若以UN，TN 和WN 表示动能的离散值、内能的改变和外力做功的改变，则要求δ =0.01 ，E 为总能量。一般可以取δ = 0.01 ，δ 要取0.05 才能满足上式时可以认为响应是不稳定的。

在装甲车辆通过不同障碍路面时，结构非线性、材料非线性与接触非线性三种类型的非线性问题都存在于工况之中：负重轮与路面的刚性接触进入大变形状态，负重轮外圈以及弹性支座等材料均为弹塑性材料，负重轮与路面以及平衡肘与限制器接触过程中存在的接触力的变化等。可以说，整车结构通过路面障碍属于典型的接触-碰撞弹塑性大变形动力学计算分析问题。

3.接触算法

在ANSYS/LS-DYNA 程序中处理不同结构界面的接触碰撞和相对滑动是程序非常重要和独特有效的功能。有二十多种不同的接触类型可供选择。主要是变形体和变形体的接触、离散点与变形体的接触、变形体本身不同部分的单面接触、变形体与刚体的接触、变形结构固连以及根据失效准则解除固连，模拟钢筋在混凝土中固连和失效滑动的一维滑动线等。

ANSYS/LS-DYNA 程序处理接触- 碰撞界面主要采用三种不同的算法，即：节点约束法、对称罚函数法和分配参数法。第一种算法现在仅用于固连界面。第三种算法仅用于滑动界面，例如炸药起爆燃烧的气体对结构的爆轰压力作用。炸药燃烧气体与被解除的结构之间只有相对滑动而没有分离。第二种算法是最常用的算法，这里对第一和第三种算法做概括介绍，重点介绍第二种算法。

不同结构可能相互接触的两个表面分别称为主表面（其中的单元表面称为主面、节点称为主节点）和从表面（其中的单元表面称为从面，节点称为从节点）。

对称罚函数法是一种新的算法，1982年8月开始用于DYNA2D 程序。其原理比较简单：每一时步先检查各从节点是否穿透主表面，没有穿透则对该从节点不做任何处理。如果穿透，则在该从节点与被穿透主表面之间引入一个较大的界面接触力，其大小与穿透深度、主面刚度成正比，称为罚函数值。它的物理意义相当于在从节点和被穿透主表面之间放置一个法向弹簧，以限制从节点对主表面的穿透。对称罚函数法同时对各主节点处理一遍，其算法与从节点一样。对称罚函数法编程简单，很少激起网格沙漏效应，没有噪声，这是由于算法具有对称性、动量守恒准确，不需要碰撞和释放条件。罚函数值大小受到稳定性限制。若计算中发生明显穿透，可以放大罚函数值或者缩小时步长来调节。它在每一时步分别对从节点和主节点循环处理一遍，算法相同，这里介绍从节点的处理方法。

对任一从节点ns 的计算步骤如下：

（1）从节点ns 搜索，确定与它最靠近的主节点ms。

（2）检查与主节点ms 有关的所有主片，确定从节点ns 穿透主表面时可能接触的主片s1 ，s2 ，…，如图2-1所示。

图2-1　从节点与主表面的接触

若主节点ms 与从节点ns 不重合，当满足式（2-6）时，从节点ns 与主单元面si 接触。

式中，Ci 和Ci+1 是主面单元上在ms 点的两条边矢量；矢量S 是矢量g 在单元上的投影；g 为主节点ms 指向从节点ns 的矢量。

如果ns 接近或者位于两个单元面交线上，上述不等式可能不确定S =

（3）确定从节点ns 在主单元面上的接触点c 位置。主单元面上任一点位置矢量可表示为

式中，是单元第j 节点的xi 坐标值：i1，i2，i3 是x1，x2，x3 坐标轴的单位矢量。接触点c（ξc，ηc）位置为下式的解：

（4）检查从节点是否穿透主面。

若l = ni[t- r（ξc，ηc）]＜0 ，则表示从节点ns 穿透含有接触点c（ξc，ηc）的主单元面。其中，ni 是接触点处主单元面的外法线单元矢量：(www.daowen.com)

如果l ≥0 ，则表示从节点ni 没有穿透主单元面，也即两物体没有发生接触-碰撞，不做任何处理，从节点ni 处理结束，开始搜索下一个从节点ni+1 。从节点与主单元面的关系如图2-2所示。

（5）若从节点穿透主面，则在从节点ns 和接触点c 之间施加法向接触力

式中，ki 为主单元面的刚度因子，有

图2-2　从节点与主面单元的关系

式中，Ki 为接触单元的体模量；Ai 为主单元面面积；Vi 为主单元体积；Li 为板壳单元最大对角线长度；f 为接触刚度比例因子，默认值为0.10。在ANSYS/LS-DYNA 计算过程中发现穿透量过大时，可以修改控制参数放大罚函数因子，但是如果取f ＞0.4 时可能会造成计算不稳定，除非减小时间步长。

在节点ns 上附加接触力矢量fs ，根据牛顿第三定律，在主单元面的接触点c 作用一个反方向的作用力- fs ，按照式（2-12）将c 点的接触力等效分配到主单元的节点上：

φj（ξc，ηc）为主单元面上的二维形函数，且在接触点c，有。

（6）计算切向接触力（摩擦力）。

若从节点ns 的法向接触力为fs ，则它的最大摩擦力为：

式中，μ 为摩擦系数。设在上一时刻（tn）从节点ns 的摩擦力为Fn ，则现在时刻（tn+1）可能的摩擦力（试探摩擦力）为F*（ F* = Fn- kΔα），k 为界面刚度，。现时刻的摩擦力由下式确定：

按照作用力与反作用力原理，计算主单元面上各节点的摩擦力。若静摩擦力系数为μd ，用指数差值函数平滑：

式中，V =Δe/Δt，Δt 为时间步长；c 为衰减因子。

由库伦摩擦造成界面的剪应力，在某些情况下，可能非常大，以致超过材料承受最大剪应力的能力、程序采用某种限制措施，令：

式中，是考虑库伦摩擦摩擦计算的tn+1 时刻的摩擦力；Ai 是主面的表面积；k是黏性系数。

（7）将接触力矢量fs ，fjm 和摩擦力矢量作为已知向量，组装到总体载荷矢量阵{P} 中，进行动力学分析。

对称罚函数方法是将上述算法对从节点和主节点分别循环处理，如果仅仅对从节点循环处理，称之为“分离和摩擦滑动一次算法”，主要应用于接触主体近似刚体的情况，可以节省计算时间。

近年来，有限元程序对表面和表面的接触算法又做了不少改进。主要是为了使薄板模压成型问题计算更加精确，例如接触搜索，采用刚体近似冲压模具时有的单元长宽比很不好。搜索与从节点最靠近的主节点时会造成困难，为了防止这类问题的发生，采用搜索最接近的主面位置来代替搜索最接近的主节点。又例如在接触计算中考虑壳单元的厚度。在薄板模压成型时，壳单元的厚度变化对接触表面摩擦力影响很大。计算中接触表面位置考虑到单元厚度。另外在接触运算中增加黏性接触阻尼项，以模拟薄板模压成型计算过程中垂直于接触表面的振荡。

4.疲劳寿命估算方法

机械零件的疲劳失效与静强度失效有本质区别。静强度失效是由于零件的危险截面上的应力大于其抗拉强度导致断裂失效，或大于屈服极限产生过大的残余变形导致失效；疲劳失效是由于零件局部应力最大处在循环应力作用下形成微裂纹，然后逐步扩展为宏观裂纹，宏观裂纹再继续扩展而最终导致断裂[39]。

在常温下工作的结构和机械的疲劳破坏取决于外载荷的大小。在循环应力水平较低时，弹性应变起主导作用，此时疲劳寿命较长，称为应力疲劳或高周疲劳；在循环应力水平较高，接近甚至超过其屈服极限时，塑性应变起主导作用，此时疲劳寿命较短，成为应变疲劳或者低周疲劳，其疲劳寿命一般低于5 ×104 次。根据不同的疲劳形式，形成两种主要的疲劳研究方法：名义应力法和局部应力应变法。

（1）名义应力法。

以名义应力为基本参数的疲劳研究方法称为名义应力法，也称S- N 曲线法或应力寿命法，是最早使用的方法。具体做法是，以材料S- N 曲线为基础，计入有效应力集中系数、零件尺寸系数、表面系数和平均应力系数等影响因素，得到零件的S- N 曲线。并根据零部件的S- N 曲线，按照疲劳累积损伤理论，进行疲劳寿命估计。

名义应力法估算构件或结构的寿命，适用于构件或结构所承受的载荷不大，断面的应力小于材料的屈服极限，应力应变为线性关系，构件及结构的寿命较大，属高周疲劳的情况。因此，疲劳寿命估算的依据是载荷谱或应力谱、S-N 曲线以及累积损伤理论等。

（2）局部应力应变法。

当应力水平较高时，零部件局部最大应力处可能会出现塑性屈服的现象。这时，只用应力参量已不能很好地表述零部件的疲劳特性。以零部件应力集中处的局部应力、应变为基本设计参数的疲劳寿命研究方法，称为局部应力应变法。局部应力应变法应用了低周疲劳的相关理论，用ε- N 曲线代替了S- N曲线，应力与应变之间的关系用循环应力- 应变曲线代替了单调的σ- ε 曲线。经研究，零部件的破坏都是从应力集中部位或应力最高处开始的，应力集中处的塑性变形是疲劳裂纹形成和扩展的主要控制参量，因此局部最大应变决定了零部件的疲劳寿命。一般，低周疲劳是指在102～105 次循环范围内的失效现象。

5.疲劳累积损伤

累积损伤是有限寿命设计的核心问题。当零部件承受高于疲劳极限的应力时，每一循环都使材料内部产生一定量的损伤。在循环载荷的作用下，疲劳损伤会不断累积，当损伤累积到临界值时会发生疲劳破坏。不同研究者根据他们对损伤累积方式的不同假设，提出了不同的疲劳累积损伤理论，归纳起来可分为三大类：

（1）线性疲劳累积损伤理论。材料在各个应力下的疲劳损伤是独立的，总损伤可以线性地累加起来。其中最有代表性的是Palmgren-Miner 理论。

（2）非线性疲劳累积损伤理论。基于假定载荷历程和损伤之间存在着相互干涉作用，即在某应力下产生的损伤与前面应力作用的水平和次数有关。其中最具代表性的是损伤曲线法和Corten-Dolan 理论。

（3）其他累积损伤理论。大都是从实验、观测和分析推导出来的损伤公式，多属于经验和半经验公式，如Levy 理论、Kozin 理论等。

由于Palmgren-Miner 线性累积损伤法则简单实用，在工程上得到了广泛的应用。Palmgren- Miner 理论认为，在疲劳试验中，材料在各个应力下的疲劳损伤是独立进行的，并且总损伤可以线性地累加起来，如图2-3所示。

图2-3　Miner 疲劳累积损伤示意图

（a）实际应力加载曲线；（b）S-N 曲线

Palmgren-Miner 理论有如下假设：

（1）载荷必须是对称循环的，即平均载荷为零；

（2）在受载过程中，每一个载荷循环都损耗一定的有效寿命分量；

（3）疲劳损伤与其所吸收的功成正比，这个功的作用循环次数和在该应力下达到的循环次数成正比；

（4）达到破坏时的总损伤量是一个常数；

（5）损伤与载荷的作用次序无关；

（6）各循环应力产生的所有损伤分量相加为常数时，就发生破坏。

如图2-3所示，应力σ1 作用n1 次，在该应力水平下达到的总循环次数为N1 ，设D 为最终失效时的损伤临界值，根据线性疲劳累积损伤理论，应力σ1每作用一次对材料的总损伤为D/N1 ，经n1 次循环作用后，对材料的损伤为n1D/N1 。同样，应力σ2 作用n2 次，产生的损伤为n2D/N2 ；应力σ3 作用n3次，产生的损伤为n3D/N3 。当各级应力对材料的损伤总和达到临界值D 时，材料即发生破坏，公式如下所示：

即

推广到更普遍的情况：