论文导读：基于PANDA框架。能够分析千万自由度规模的弹塑性静力学问题。非线性求解策略。形成了面向对象有限元并行计算框架PANDA。并行计算，基于PANDA框架的非线性静力学有限元。

　　关键词：PANDA，静力学，非线性，有限元，并行计算

　　1 引言

　　特种武器结构复杂，在整个库存到靶序列（Stockpile to TargetSequence，STS）全寿命周期内要经历复杂严酷的载荷和环境条件，结构响应呈现出高度的材料非线性、边界非线性和几何非线性。为提高特种武器的设计、试验和库存维护水平，对武器结构在各种条件下响应的精细建模和分析至关重要，需要充分考虑结构的几何细节和物理内涵，所建立的有限元模型可达上千万自由度规模乃至更高，而传统的商用有限元程序由于国外对我国的出口限制，非线性有限元模型的分析规模被限制在几百万自由度以下，且计算周期较长，无法快速响应设计和维护的需要。

　　为了提升特种武器的工程数值模拟能力，适应不断提高的武器工程数值模拟需求，迎接和加速由现阶段小规模低效率计算向大规模高效并行计算的转变，2007年中国工程物理研究院启动了院预研重大项目“武器工程大规模并行计算框架研究及基础平台开发”。该项目在已有源码程序的基础上，通过在有限元并行计算方法方面开展研究与软件开发，初步形成了面向对象有限元并行计算框架PANDA，并基于PANDA框架初步开发了可应用于部分静力、振动、冲击和传热武器工程问题求解的大规模有限元并行计算模拟程序。

　　针对特种武器研制中的非线性静力学有限元大规模精细分析需求，充分消化吸收开放源代码的程序设计思想和技巧，基于PANDA框架，开发非线性静力学有限元分析所需的单元类型、材料模型、非线性并行求解策略，集成大规模线性方程组并行求解算法，初步形成了可求解小应变、有限应变线弹性和弹塑性静力学问题的非线性静力学程序。悬臂梁弹塑性有限元分析模型达到了千万自由度规模，并行求解时间低于一小时。本文介绍了基于PANDA框架的单元类型、材料模型、非线性求解策略设计，并初步验证了非线性静力学有限元并行计算程序的计算精度和千万自由度规模分析能力。

　　2 基于PANDA框架的非线性静力学有限元并行计算程序设计

　　通过中国工程物理研究院的预研重大项目，采用面向对象、层次化、组件化的设计思想，对工程结构非结构网格有限元分析程序的基本数据结构、并行通信、求解控制等方面的共性和可重用部分进行抽象和程序实现，并集成了区域分割、解法器等服务组件，形成了面向对象有限元并行计算框架PANDA，提供经过系统规划设计的应用程序开发接口，以提供服务的形式引导应用程序的设计和实现，初步建立了结构分析有限元并行计算应用程序的集成开发环境。科技论文，并行计算。

　　基于PANDA框架，结构分析有限元并行计算应用程序的开发工作变得较为简单和高效，程序开发工作量大为减少。在PANDA框架既设的应用软件架构下，应用程序开发者可以将精力集中到本应用程序独有的个性部分，并充分利用框架中集成的经过充分验证的高效解法器等服务组件，场和网格数据的组织、存储和管理由框架负责，应用程序开发者无需关心其底层数据结构等实现细节。对非线性静力学有限元并行计算程序的开发而言，通过继承PANDA框架中场、节点、单元类型、材料模型、空间积分器、求解控制等共性部分，开发适用于非线性静力学有限元计算的场、节点、单元类型、材料模型和非线性求解策略，使用由框架提供的输入参数解析、并行通信、区域分割、线性方程组解法器等组件，就可较快速形成可求解大规模非线性静力学问题的高性能有限元并行计算程序。

　　2.1 单元类型和材料模型

　　PANDA框架中的ElementBaseT基类抽象出了各种单元类型的共性（属性和操作），在应用程序中，具体单元类型都可由该基类逐级派生并添加自身特有的属性和操作得到，添加新单元类型较为方便。目前，在PANDA非线性静力学有限元并行计算程序中已基于PANDA框架开发了结构分析中常用的八节点六面体小应变单元和更新拉格朗日有限应变单元，后者在每一载荷增量步对几何构形进行更新。

　　基于PANDA框架，应用程序中具体材料模型在各个类继承层次上定义，在程序结构设计时便于添加各种具体的材料模型。目前已实现了适用于小应变单元和更新拉格朗日单元的二维、三维各向同性线弹性和弹塑性材料模型。其中的弹塑性材料模型是线弹性模型与J2各向同性屈服条件的组合，塑性变形阶段的当前屈服应力由下述四种各向同性塑性硬化函数之一进行描述。

　　（1）线性各向同性硬化函数

　　……（1）

　　上式中为初始屈服应力，为硬化模量

　　K＝（EET）/（E-ET）……（2）

　　其中，E为弹性模量，ET为割线模量。科技论文，并行计算。

　　当硬化模量K为0时，采用线性各向同性硬化函数的小应变各向同性弹塑性材料模型就成为小应变理想弹塑性材料模型。

　　（2）linear-exponentialsaturation各向同性硬化函数

　　并行计算…（3）

　　上式中为初始屈服应力，为线性硬化模量。为saturation硬化模量，称为saturation应变，是一特征等效塑性应变量，用于描述等效塑性应变对saturation硬化的影响，的值越大，相同等效塑性应变时的当前屈服应力越小。

　　随着等效塑性应变值增大，linear-exponentialsaturation各向同性硬化函数逐步趋近于由并行计算表示的线性各向同性硬化函数（其初始屈服应力为）。

　　在ANSYS中，Nonlinear Isotropic Hardening （NLISO）材料模型对塑性硬化规律的描述与此相似。

　　（3）general cubic spline（立方样条）各向同性硬化函数

　　并行计算（4）

　　上式中，n为用于定义样条的点数（对应于样条的n-1个子区间和2个端点外区域）。在每个子区间上，样条具有下面的形式：

　　（5）

　　上式中，。从样条点计算系数a（i）的条件为：样条通过这些数据，且样条的一阶导数在子区间之间连续。此外，样条的端条件被规定为：，以满足样条点数据范围之外的线性延伸。科技论文，并行计算。

　　尽管不是必须和强制性的，人们一般选择，以便用定义初始屈服应力。

　　（4）power law（幂函数）各向同性硬化函数

　　……（6）

　　由上式可知，初始屈服应力。

　　2.2 非线性求解策略

　　PANDA非线性静力学有限元并行计算程序实现了标准牛顿法、带线性搜索的非线性牛顿法、带线性搜索的非线性预处理共轭梯度法等三种策略。它将整个载荷划分为多个载荷增量（可根据迭代求解过程的情况进行自动增减），在每个载荷增量内进行反复迭代求解直至满足某个收敛准则，而每一迭代过程中通过调用PANDA框架中集成的线性方程组解法器组件实现线性方程组的求解（线性方程组的求解可选用直接解法或迭代解法）。三种求解策略使用最多17个参数对非线性求解的迭代过程、精度、自动载荷增量步、线性搜索、预处理过程等进行控制，都使用基于残余力范数的两个收敛评估准则（绝对容差和相对容差，后者是两次迭代的残余力范数的比值），只要满足其中一个准则，一次非线性迭代结束。

　　3 测试算例

　　为了初步验证基于PANDA框架开发的非线性静力学有限元并行计算程序的计算精度、并行计算规模和计算速度，我们用该程序分别计算了杆、悬臂梁和及来自于工程实际问题的某离心机的静力学响应，并与理论解或商用有限元程序计算结果进行了对比。对比结果表明，该并行程序具有很好的计算精度和计算速度，能够分析千万自由度规模的弹塑性静力学问题。3.1 杆的均匀拉伸弹塑性加、卸载响应

　　杆尺寸为1m×1m×10m，沿其轴向承受8.5×107N拉力，杆材料依次用小应变、有限应变各向同性线弹性和双线性弹塑性材料模型进行描述（屈服应力为80MPa），结构用与材料模型相匹配的八节点六面体小应变单元或有限应变单元建模。分析了结构分别在加载和卸载状态下的位移和应力。

　　采用八节点六面体小应变单元、小应变线弹性或双线性弹塑性材料模型的计算结果与同等有限元模型规模的ANSYS计算结果完全一致；而采用有限应变单元、有限应变线弹性或双线性弹塑性材料模型的计算结果与同等有限元模型规模的ANSYS计算结果相比，轴向应力相差最大0.046%，轴向位移相差最大0.078%。

　　均匀拉伸载荷作用下杆位移和应力响应的有限元分析测试算例表明，PANDA静力学有限元程序中的线弹性和双线性各向同性硬化弹塑性材料模型具有很高的计算精度。

　　3.2 某离心机的静力学线弹性分析

　　建立了来自于实际工程研制工作中的某型号离心机（图1）的ANSYS和PANDA静力学线弹性分析有限元模型。PANDA模型的自由度为25万，由8节点六面体单元及其退化单元组成。科技论文，并行计算。计算了该离心机在重力作用下结构的变形。PANDA静力学程序计算结果（见图2）与ANSYS计算结果的对比分析表明，两者之间分布相似，最大位移相差4.1%。

　　图1 某型号离心机的PANDA有限元模型（一半结构）

　　Fig.1 PANDA finite element model of a centrifuge（half of the structure）

　　（a）PANDA（b） ANSYS

　　图2 某型号离心机的静力学有限元分析位移结果

　　Fig.2 The displacement of a centrifuge bystatics FEM analysis

　　3.3 悬臂梁的千万自由度双线性弹塑性分析

　　悬臂梁尺寸为1m×1m×10m，约束固定端的全部位移自由度，在自由端面施加7.5MPa的切向面载荷（梁中剪力为7.5×106N，固定端弯矩为7.5×107N·m），材料的弹性模量为210GPa、泊松比为0.3、屈服应力为300MPa、割线模量为6.1GPa。建立的PANDA有限元模型采用八节点六面体有限应变单元、有限应变双线性弹塑性材料模型。共分析了三种规模的有限元模型，其节点数分别为76400、674081和3385900（对应的自由度数分别为228000、2017200和10143000）。相应的ANSYS有限元模型具有12221个节点。

　　表1 不同规模时悬臂梁的最大横向位移及其相对于ANSYS的误差

　　Table 1 Maximum vertical displacement of thebeam and the relative error

　　自由度数228000201720010143000

　　横向位移计算值（mm）0.1632980.1644210.164725

　　相对于ANSYS的误差（%）-0.7928-0.11060.0741

　　（a）横向位移分布云图（m）

　　（a） The vertical displacement （m）

　　（b）von Mises等效应力云图（Pa）

　　（b） The von Mises equivalent stress （Pa）

　　图32017200自由度时悬臂梁的弹塑性静力学响应

　　Fig.3 The elasto-plastic statics response of thebeam with 2017200 dofs

　　PANDA静力学有限元计算结果表明，对同一规模的有限元模型，采用不同并行进程数（1～128个进程）时的计算结果完全相同。不同规模时悬臂梁的最大横向位移计算结果如表1所示。图3为674081个节点（2017200个自由度）时的PANDA弹塑性静力学计算结果。由于我们进行PANDA非线性静力学有限元计算时所取的非线性迭代求解收敛准则值较大，使得同一规模时其计算结果精度比ANSYS略低。各个规模下的计算结果均具有较高的计算精度，与ANSYS相比，相对误差均较小，且随着模型规模增大，计算精度逐步提高。

　　PANDA非线性静力学并行计算在16节点双路四核的曙光小型服务器上进行。3385900个节点（10143000个自由度）的双线性弹塑性有限元静力学并行计算在采用128进程时主进程的运算时间为2468.95秒（不包括区域分割计算的时间）。科技论文，并行计算。

　　4 结语

　　基于PANDA框架，开发非线性静力学有限元分析所需的单元类型、材料模型、非线性并行求解策略，集成大规模线性方程组并行求解算法，初步形成了可求解小应变、有限应变线弹性和弹塑性静力学问题的非线性静力学程序。杆、来自于工程实际问题的某离心机和悬臂梁测试算例验证了该程序的千万自由度规模弹塑性静力学高效并行计算能力及计算精度，表明该程序具有了初步应用于部分实际工程结构的大规模高效非线性静力学有限元分析的能力。科技论文，并行计算。

　　为了进一步提升和丰富PANDA非线性静力学有限元并行计算程序的分析功能和计算效率，基于PANDA框架，我们可进一步丰富结构分析常用的其它单元类型、材料模型和非线性求解策略，比如接触单元、多尺度关联单元、超弹性材料模型、粘弹性材料模型、粘塑性材料模型、新型的高效精细非线性迭代求解策略、高精度高效线性方程组求解算法，从而为武器工程结构的精细建模和分析奠定条件，跨越式地提升我国的武器工程结构以及其它重大工程与装备的数字化设计水平。

　　参考文献

　　[1]王勖成，邵敏。有限单元法基本原理和数值方法。北京：清华大学出版社，1997.3第2版。

　　[2]Sandia National Laboratory. TahoeUser Guide - Input Version 3.4.1. May 21, 2003.