开通VIP，畅享免费电子书等14项超值服
首页
好书
留言交流
下载APP
联系客服
2023.07.04湖北
目前，数值计算方法在科学研究和工程技术中得到了广泛应用，与理论研究和实验研究组成了现代科学技术的三大支柱，并具有快捷、安全和低成本的优势。然而，超高速碰撞、冲击侵彻、爆炸、金属冲压成形、裂纹动态扩展、流固耦合、应变局部化等强非线性问题涉及材料的特大变形、破碎甚至熔化和汽化，包含了一系列复杂的物理过程，给数值计算方法带来了巨大的挑战。
根据所采用的运动描述方法，数值计算方法可以划分为拉格朗日法、欧拉法和混合法三大类。
1.1拉格朗日法
在拉格朗日法中，计算网格固连在物体上，随其一起变形，因此材料与网格之间不存在相对运动(即迁移运动，也称对流运动)，控制方程中不存在对流项，大大地简化了控制方程及其求解过程。单元质量在计算过程中始终不变，但其体积由于网格变形而不断发生变化。拉格朗日法具有如下优点。
(2)网格线（面）和物体的外表面及材料界面在求解过程中始终重合，因此拉格朗日法易于处理边界条件和跟踪材料界面。
图：拉格朗日网格
上图给出了一个典型的拉格朗日法的例子。如前所述，拉格朗日法的计算网格随物质一起变形扭曲。当计算网格严重扭曲甚至相互重叠时，将引起较大的误差，甚至使计算失败。这个问题可以通过使用滑移界面(slidelines)和重分网格(rezone)技术在一定程度上
得以克服。在物体之间可能产生相互滑移的问题中，需要使用滑移界面单元，如流体与固体壁之间的相互作用问题、子弹穿透靶体过程和碰撞物体之间的相互接触问题等。滑移界面在离散可能出现非常大的剪切变形和失效区域也是很有效的。大多数滑移界面元都是将加速度和速度分解为与界面相切的切向分量和与界面垂直的法向分量。当物体之间相互接触时，法向运动是连续的，而当物体之间脱离接触后，法向运动是完全独立的。当物体之间相互分离或接触面之间无摩擦时，切向运动是相互独立的，但当物体相互接触后并存在摩擦力时，切向运动需满足一定的关系。两个相互接触的物体在满足一定条件下可以相互分离,而已分离的两个物体可以再次相互碰撞。
而每个旧网格都有不同的历史，因此在变换内变量时也将遇到困难。在得到一个成功的新网格之前，一般需要多次网格重分过程，而且最终结果的精度在很大程度上取决于计算者的经验。
克服网格畸变的另外一种方法是采用侵蚀(erosion)算法，当单元的变形过大，其拉力或等效塑性应变达到一定阈值时就将该单元删除。侵蚀算法在侵彻问题的模拟中取得了一定的成功，但也带来新的问题，如单元侵蚀造成的质量不守恒给系统带来扰动，忽略了被侵蚀单元对系统的作用，进而造成计算误差：另外被侵蚀的单元也无法继续参与描述材料的变形，无法模拟超高速碰撞中碎片云的形成过程及其作用效果。
HEMP程序是美国劳伦斯利福摩尔国家实验室（LawrenceLivermoreNationalLaboratory,LLNL)的Wilkins在20世纪50年代到80年代研发的拉格朗日型显式有限差分程序，早期被用于计算爆炸、冲击问题。EPIC6)(ElasticPlasticImpactComputations).是由Johonson等在20世纪70年代研制的拉格朗日型显式积分有限元程序，采用三角形单元和四面体单元，并引入了重分区(rezoning)和侵蚀(erosion)算法，可以计算高速碰撞和爆炸等问题。EPIC程序的发展略早于DYNA程序，它们的区别在于EPIC主要用三角形单元，而DYNA主要用四边形单元4.EPIC在研究领域一直发挥着作用，如美国陆军研究实验室用EPIC研究材料的强度和局部化行为。
DYNA系列的程序是最为知名的显式动力学有限元分析程序，现在使用的很多程序或软件(比如MSC.DYTRAN)都是在DYNA程序的基础上发展而来的。最初的DYNA3D程序源自20世纪70年代的美国LLNL国家实验室，其发展背景和PRONTO3D程序类似，但其商业化版本LS-DYNA取得了更为广泛的应用。一些商业有限元软件（如ANSYS)对LS-DYNA求解器进行了包装，提供了部分前后处理功能。Aup公司的OASYS也为LS-DYNA提供了前后处理环境。对早期各种程序更详细的总结和介绍可参见Zukas主编的《HighVelocityImpactDynaimcs》的第9章。
1.2欧拉法
对于可能产生严重网格扭曲的问题，或者是在相互分离的材料将发生混合的问题中，需要使用欧拉法。在欧拉法中，计算网格固定在空间中，不随物体运动，而材料相对于网格运动(如下图所示)，因此不存在网格畸变问题。各时刻的速度、压力、密度和温度等物理量是在空间点上计算的，而不是像拉格朗日法那样在物质点上计算，因此质量、动量和能量等物理量将跨越单元边界在单元间输运。各单元的体积在计算过程中保持不变，但单元质量由于输运而不断发生变化。欧拉方法比较适合模拟大变形问题，因此在计算流体力学中多采用欧拉方法，早期用于冲击爆炸问题的流体动力学程序也多采用欧拉方法。
图：欧拉网格
欧拉法只计算质量、动量和能量等物理量在单元间跨越网格边界的输运量，难以精确确定材料界面和自由表面的位置，因此欧拉法施加边界条件比拉格朗日法更困难，精度也更低。除非使用特殊的材料分界面定位技术，否则材料分界面将很快在计算网格中弥散。目前已提出了多种方法，如通过引入无质量的拉格朗日示踪点来计算材料分界面的实际位置。这些示踪点也可用于提供材料历史数据，但将大大地增加计算的复杂程度和计算量。
北京理工大学宁建国研究组开发了使用欧拉网格的有限差分程序MMIC2D和MMIC3D,引入了界面处理技术，并用于聚能射流等问题的模拟，同时，为了配合求解程序的使用，也开发了前处理程序MESH2D以及可视化程序VISC2D和VISC3D.北京计算数学与应用物理研究所进行了大量流体动力学程序的研究开发工作，开发了欧拉网格的有限差分程序MEPH2D和MEPH3D29,3o,并引入自适应加密网格技术和界面处理技术，用于高速碰撞、炸药爆轰、聚能射流等问题的模拟。中科院力学所张德良和北京大学刘凯欣等基于高精度时空守恒元解元算法开发了SUPERCE/SE程序用于高速侵彻和爆轰问题的模拟。杨秀敏院士课题组基于高精度、高分辨率差分格式WENO和贴体坐标编制了空气冲击波流场计算程序EF3D,可以研究不同爆炸方式、不同入射反射背景条件下的各种空气波流场参数变化规律。
1.3混合方法
拉格朗日法和欧拉法都存在严重的缺陷，但也都具有各自的优势，如果能将二者有机地结合，充分吸收各自的优势，克服各自的缺点，则可解决一大批只用拉格朗日法或欧拉法所解决不了的问题。
1.3.1任意拉格朗日-欧拉法
任意拉格朗日-欧拉法(arbitraryLagrangian-Eulerianmethod,ALE)最早是由Noh(1964年)以耦合欧拉-拉格朗日法的术语提出的，并用有限差分法求解带有移动边界的二维流体动力学问题。在Noh的研究工作中，网格点可以随物质点一起运动，但也可以在空间中固定不动，甚至网格点可以在一个方向上固定，而在另一个方向上随物质点一起运动，因此ALE描述也被称为耦合欧拉-拉格朗日描述。例如在液体表面波的传播问题中，网格点在垂向随物质点一起运动，而在水平方向上固定不动，如图所示。这样很容易描述液体表面的运动，而且网格不会发生扭曲。
图：流体表面波传播问题
在ALE法中，计算网格可以在空间中以任意的形式运动，即可以独立于物质坐标系和空间坐标系运动。这样通过规定合适的网格运动形式可以准确地描述物体的移动界面，并维持单元的合理形状。类似于欧拉法，在ALE法的控制方程中也将出现对流项，因此使得系数矩阵不对称，还可能得到振荡解，需要进行相应的数值处理，如使用迎流权函数或彼得洛夫-伽辽金法(Petrov-Galerkin)法。其实质是在原微分方程中增加一扩散平衡项。纯拉格朗日法和欧拉法实际上是ALE法的两个特例，即当网格点的运动速度等于物质点的运动速度时就退化为拉格朗日法，而当网格固定于空间不动时就退化为欧拉法。
ALE法最早是为了解决流体动力学问题而引入的，并且使用有限差分法。由于核反应堆结构安全分析的需要，ALE法被引入有限元法中，用以求解流体与结构相互作用问题，目前已被应用于固体力学领域中求解大变形问题，如接触碰撞、弹性断裂力学和加工成形等。
1.3.2质点网格法(PIC)
为了解决PIC法中的粒子脉动和存储量大的问题，在Harlow的指导下，Gentry、Martin和Daly于1966年进一步提出了流体网格法(fluid-in-cellmethod),简称FLIC法。它和PIC法一样采用欧拉网格，但在第二步中不计算质点的迁移，而是计算连续流体的迁移，即先算出通过网格边界的质量输送量，得出每个网格的新密度，再算出通过网格的质量所携带的动量和能量的输送量，最后得到每个网格的新速度和能量。FLIC法还有一套局部网格单元的计算格式，能计算一些边界形状比较复杂的问题。
为了模拟不可压缩自由表面流动问题，Harlow和Velch在PIC法的基础上进一步建立了标记网格法(markerandcellmethod,MAC)。MAC法仍然采用欧拉矩形网格单元，对纳维-斯托克斯方程则用差分近似，而把压力和速度分量作为基本未知量。此外，MAC法还在网格中布置适量的无质量标记点，并在整个计算中跟踪每个标记点，以判定网格里有哪种流体存在。MAC法可用于计算多种流体和带有自由面的问题，是模拟不可压缩流动问题的一个成功的数值方法。
经典的PIC法不是一种完全的拉格朗日质点方法，它只在质点上存储质量和位置信息，而其他物理量仍然存储在网格上。动量在网格和质点之间的迁移将产生较大的数值耗散，严重损害了PIC的计算精度。有两条途径可以提高PIC的精度：采用更精细的动量迁移算法，或者令质点携带所有物质信息。Brackbill等将流体的所有物质信息均赋予质点，基于有限差分法发展了一种用于求解流体流动问题的低数值耗散的PIC方法——FLIP(fluid-implicit-PIC)方法。FLIP是一种完全的拉格朗日质点格式。
1.3.2物质点法(MPM)
图：物质点法示意图(a)材料区域；(b)质点代表区域；(c)物质点法离散
图：MPM与SPH计算效率比较
图：MPM与SPH求解Taylor杆冲击的比较
在物质点法中，动量方程在计算网格结点上求解，因此需要建立计算网格单元的质量阵，既可采用集中质量阵，也可采用协调质量阵。Tan等在计算裂纹尖端的能量释放率时指出协调质量阵能够获得更好的能量守恒性质：Lov等分析了质量阵形式对能量和动量守恒性的影响，结果表明协调质量阵有较好的能量守恒性。
Bardenhagen等发现，物质点法中形函数导数在计算网格结点处不连续，使得质点穿越计算网格时产生较强的数值噪声，导致应力振荡。为此，他们采用彼得洛夫-伽辽金法对MPM进行了扩展，提出了广义插值物质点法（generalizedinterpolationmaterialpointmethod,GMP),建立了一族具有C1连续性的形函数，有效减弱了数值噪声。在GMP中，质点的影响范围由质点的特征函数表示。在计算过程中，需要实时计算该质点的形状以判断该质点的影响范围。为了避免实时计算质点形状的困难，一般假定质点的形状不变，通过其现时位置和初始形状大小判断其影响范围，该方法称为uGIMP。Zhang和Ma等通过将标准物质点法的形函数导数值与FLIP中基于背景网格结点的形函数导数值加权平均，建立了一种新的形函数导数计算方式，以降低质点穿越网格引起的噪声。该方法不需要跟踪计算质点的形状，并且质点仅影响其所在的单元，因此便于非结构化计算网格的计算。
物质点法采用质点离散，各质点之间无拓扑关系，因此易开展h自适应分析。Tan等建立了质点和计算网格的层级结构化自适应分析算法，通过局部加密质点和计算网格计算了裂纹尖端的能量释放率。针对冲击爆炸中易出现的非物理现象的数值断裂问题，马上和张雄等提出了质点自适应分裂的物质点法，有效避免了数值断裂问题。
在并行算法方面，Parker等基于消息传递模型(MP)研究了MPM的并行算法，开发了Uintah程序。黄鹏和张雄等提出基于信息共享模型(OpenMP)的MPM并行算法，张衍涛和张雄等进一步建立了基于OpenMP的区域信息交错更新的MPM并行算法，精巧地避免了数据竞争问题。Ma和Lu等基于并行程序开发平台SAMRAI2建立了GIMP的并行算法。至此，已经建立了基于两种并行机制的MPM并行算法。
经过不断发展和开发，物质点法已经在超高速碰撞、冲击侵彻、爆炸、裂纹扩展、材料破坏、颗粒材料流动和岩土冲击失效等一系列涉及材料特大变形的问题中体现了优势。
在超高速碰撞问题方面，马上和张雄等以空间碎片防护问题为背景采用MPM研究了弹丸碰撞薄板和厚板等一系列问题，黄鹏和张雄等则进一步采用并行的MPM方法精细模拟了超高速碰撞问题中的碎片云问题。宫伟伟和张雄等基于泡沫铝的Micro-CT扫描图片重构了具有真实微观结构的三维泡沫铝物质点模型，并模拟了两种含有泡沫铝的Whipple防护结构的超高速碰撞问题，计算结果与实验结果吻合。
在冲击侵彻问题方面，Sulsky等模拟了Taylor杆问题、金属成形问题和钢球侵彻铝靶体问题；黄鹏和张雄等进一步采用物质点接触算法模拟了中低速冲击侵彻问题，所得计算结果与实验结果吻合；马志涛和张雄等基于局部多重网格物质点接触算法模拟了同样的弹体侵彻问题；宫伟伟和张雄等将Deshpande-Fleck泡沫模型引入到物质点法中，对泡沫铝材料动态压缩性能进行了三维仿真模拟：王晓军和张雄等将JH2、HJC和RHT本构模型引入到物质点法，模拟了陶瓷和混凝土的冲击问题。
在爆炸问题方面，Hu和Chen83基于物质点法模拟了爆炸对混凝土墙的破坏作用：马上和张雄等将MPM应用于高能炸药爆轰问题；王宇新等采用MPM系统研究了爆炸焊接问题；张忠等基于MPM分析了非均质固体炸药的爆轰过程及其对金属材料的破坏作用。
在裂纹扩展问题方面，Nairn的研究组采用MPM做了大量工作，模拟了二维、三维的裂纹扩展问题。采用多速度场描述裂纹面以避免裂纹面之间的穿透，在物质点法中计算了裂纹尖端各种参数，如动态J积分、I、Ⅱ和I型裂纹的应力强度因子,模拟了脆性材料的断裂问题并计算了材料断裂过程中的能量释放率。此外，Gilabert等采用MPM模拟了陶瓷材料中裂纹萌生和扩展过程；Daphalapurkar等在GMP中实现了内聚力模型(cohesivezonemodel)模拟了延性和脆性材料中IⅡ型裂纹动态扩展的问题；Wang等采用非规则的计算网格计算了二维混合型裂纹扩展问题。
在模拟材料失效问题方面，Chen的研究组做了大量工作，采用MPM模拟了冲击载荷下脆性材料的动态失效问题、受压缩的薄膜脱层问题和材料在局部加热情况下的失效问题;Shen采用MPM模拟了玻璃在冲击作用下的失效破碎问题；Schreyer和Sulsky等在MPM中采用脱聚本构模型(decohesionconstitutivemodel)研究了复合材料中的脱层失效模式；Ionescu等采用MPM模拟了各项异性的软组织失效问题；Li和Pan等基于MPM研究了脆性材料在冲击作用下的失效问题。
在岩土冲击动力学问题方面，Coetzee等采用二维物质点法模拟了岩土中错墩受力问题；Andersen等采用GIMP分析了边坡失效问题；黄鹏采用MPM模拟了半球壳体侵彻岩土和边坡失效问题。
在颗粒物流动问题方面，Wieckowski等采用MPM模拟了储料垛内的颗粒流动问题；Bardenhagen等采用MPM模拟了颗粒物受剪作用下的流动过程；Coetzee等基于MPM模拟了挖斗铲玉米颗粒堆的过程；王津龙等基于MPM模拟了柱状堆石成形过程。
此外，MPM在海冰动力学问题、切削问题、碳纳米管增强型复合材料、薄膜流固耦合问题和多孔介质的动力学响应分析等问题中也有应用。
为综合物质点法和其他方法的优势，一些学者将物质点法与其他方法耦合。Guilkey建立了多物质欧拉方法和MPM的耦合方法，分别采用多物质欧拉算法模拟流体区域、物质点法模拟固体结构。Gilmanov和Acharya将MPM与浸没边界法进行了杂交，其中MPM用于求解固体结构。Guo和Yang采用握手区方法耦合了MD和MPM,模拟了Cu-Cu和Si-Si的群簇高能碰撞问题。Lu等通过逐步细化MPM计算网格间距过渡到分子尺度，在过渡区建立质点和分子的一一对应关系，将MPM与MD进行无缝耦合，分析了硅在纳米尺度下的拉伸过程。Chen等采用耦合的MD/MPM进行了纳米尺度的切削和薄膜生产过程的多尺度分析。
求解格式，GMP算法，接触算法，自适应算法，有限元法，杂交物质点有限元法，耦合物质点有限元法和自适应物质点有限元法，具有多种材料模型、状态方程和失效模型，并分别利用OpenMP和MPI实现了SMP和MPP两种体系下的并行化。2009年，美国的ADI公司发布了ASTE-P1.0软件，将物质点法加入到了FSI、CSD、CFD模块，用来计算气动弹性问题。2010年，英国的JveerAerospaceLtd公司推出了商业化的物质点法计算软件MPMsim。
注：本文摘自清华大学张雄教授的专著——《物质点法》，如有侵权，请联系删除。