基于PANDA框架的非线性静力学有限元

基于PANDA框架的非线性静力学有限元
论文导读：基于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框架开发的非线性静力学有限元并行计算
程序的计算精度、并行计算规模和计算速度，我们用该程序分别计算了杆、悬臂梁和及来自于工程实际问题的某离心机的静力学响应，并与理论解或商用有限元程序计算结果进行了对比。

对比结果表明，该并行程序具有很好的计算精度和计算速度，能够分析千万自由度规模的弹塑性静力学问题。