摘要：利用ANSYS和ABAQUS软件，对内压作用下的柱壳开孔模型进行了有限元计算，研究了不同的软件、不同的单元类型、不同的单元选项对应力分析结果的影响。结果表明，应用有限元软件对压力容器进行应力分析时，单元类型的选择、单元选项的设置对计算结果的影响是比较大的，对同一模型的最大相对误差能达到9%左右。
9%差别很大吗？作为一个挑剔的软件用户，肯定感觉大。作为一个曾经参与过有限元求解器开发的程序员，这两个求解器在这个问题下的对比结果已经非常不错了。15年前我们自主研发的处女作求解器诞生的时候，跟ANSYS的误差30%多，那时候我们都兴奋的睡不着觉！
目前在国内压力容器界，对压力容器进行应力分析时，一般采用ANSYS或ABAQUS软件，这两种软件均是世界上知名的应力分析软件，均提供了三维实体单元。
在压力容器应力分析中，通常采用实体单元对结构进行网格划分，但不同的软件提供的实体单元有着不同的单元特性选项，这些单元特性选项对应力分析最终结果的影响是较大的，对于不同的受压元件、不同的求解问题应选择合适的单元才可以得到理想的结果。
1 ANSYS单元特性
ANSYS软件提供了8节点三维实体单元Solid185和20节点三维实体单元Solid186，这两种单元的有着不同的单元选项，具体如下：
1.1 Solid185单元（8节点三维实体单元）
solid185单元用于构造三维固体结构。单元通过8个节点来定义,每个节点有3个沿着xyz方向平移的自由度。单元具有超弹性，应力刚化，蠕变，大变形和大应变能力。还可采用混合模式模拟几乎不可压缩弹塑材料和完全不可压缩超弹性材料。
SOLID185的单元选项KEYOPT(2)的设置控制单元在求解时的选项：
0 --使用B方法全积分（默认选项）
1 --带沙漏控制的一致缩减积分
2 --增强应变公式
3 --简单增强应变公式
以上选项具体意义如下所述：
1）B方法全积分(选择缩减积分方法)：
这种方法有助于在几乎不可压缩情况下预防体积锁定，这种方法在单元的高斯积分点上用体积应变代替平均体积应变。然而，此种方法不能预防任何在弯曲问题时的剪切锁定，在这种情况下，应该采用应变增强模式。如果不能十分确定弯曲行为是否发生，推荐使用应变增强模式。
2）一致缩减积分
因为只有一个积分点，所以此种方法也可预防在几乎不可压缩情况下的体积锁定。而且比B方法(选择缩减积分方法)更有效。然而，为了控制沙漏而引入的人工应变能可能会对结果精度产生影响。
当使用这种方法时，通过比较总应变能(ETABLE中的SENE)和控制沙漏引入的人工应变能(ETABLE中的SENE)来检查求解是否正确，如果人工应变能和总应变能的比值小于5%，结果一般可以接受，如果比值大于5%，重新细划网格。可以通过在求解阶段使用命令OUTPR,VENG来监控总应变能和人工应变能。
3）增强应变公式
为了预防弯曲变形时的剪切锁定和几乎不可压缩时的体积锁定，此种方法引入了13个内部自由度（用户不可见），如果混合应变模式和增强应变公式同时使用，仅9个为克服剪切锁定的自由度被引入，所有的内部自由度自动在内部被凝聚。因为额外的内部自由度和静态凝聚，这种方法没有B方法(选择缩减积分方法)和一致缩减积分方法有效。
4）简单增强应变公式
这是一个特殊的增强应变公式，为了预防弯曲变形时的剪切锁定，总是有9个内部自由度被引入（用户不可见）。因为没有内部自由度处理体积锁定，当材料是几乎不可压缩时，此种方法不可用。当和混合u-P模式一起使用时，简单增强应变公式和增强应变公式得到了一样的结果。所有的内部自由度自动在内部被凝聚。因为额外的内部自由度和静态凝聚，这种方法没有方法(选择缩减积分方法)和一致缩减积分方法有效。但由于使用了较少的内部自由度，比增强应变公式有效。
1.2 Solid186单元（20节点三维实体单元）
solid186是一个高阶3维20节点固体结构单元，SOLID186具有二次位移模式，可以更好的模拟不规则的网（例如通过不同的CAD/CAM 系统建立的模型）。 单元通过20个节点来定义，每个节点有3个沿着xyz方向平移的自由度。Solid 186可以具有任意的空间各向异性，单元支持塑性，超弹性，蠕变，应力钢化，大变形和大应变能力。还可采用混合模式模拟几乎不可压缩弹塑材料和完全不可压缩超弹性材料
SOLID186的单元选项KEYOPT(2)的设置控制单元在求解时的选项：
0 -- 一致缩减积分（默认选项）
1 -- 全积分
1）一致缩减积分
在预防几乎不可压缩时的体积锁定是有效的，然而，如果在每个方向上不能保证至少2层单元时，沙漏模式可能在模型中传播。
2）全积分
全积分方法不会引起沙漏模式，但在材料不可压缩时可能导致体积锁定，此方法常用于纯线性分析，或者当模型在每个方向只有一层单元时。
压力容器工程师在使用ANSYS软件进行应力分析时，基本采用的均是这种单元的默认值。
2 ABAQUS单元特性
ABAQUS也同样提供了8节点三维实体单元（C3D8，C3D8R）、20节点三维实体单元（C3D20，C3D20R），这两类单元又按缩减积分、全积分分为两种类型。
缩减积分单元类型为C3D8R、C3D20R，C3D8R单元具有8 个节点，C3D20R单元具有20个节点，R表示该单元为缩减（Reduced）积分单元。通常而言，减缩积分可以满足要求，并且可以避免使用完全积分产生的刚度过大，应力不连续等问题。
全积分单元类型为C3D8、C3D20，C3D8单元具有8 个节点，C3D20单元具有20个节点。
减缩积分单元比完全积分单元在每个方向少用一个积分点。减缩积分的线性单元只在单元的中心有一个积分点。实际上，在ABAQUS中这些一阶单元采用了更精确的均匀应变公式，即计算了单元应变分量的平均值。
ABAQUS在一阶减缩积分单元中引入了一个小量的人工“沙漏刚度”以限制沙漏模式的扩展。在模型中应用的单元越多，这种刚度对沙漏模式的限制越有效，这说明只要合理地采用细划的网格，线性减缩积分单元可以给出可接受的结果。对多数问题而言，采用线性减缩积分单元的细划网格所产生的误是在一个可接受的范围之内。建议当采用这类单元模拟承受弯曲载荷的任何结构时，沿厚度方向上至少应采用四个单元。
3 分析实例
为比较不同的软件、不同的单元类型、不同的单元选项对应力分析结果的影响，在本节中，以压力容器中较为常见的柱壳径向开孔受内压的模型对上述各情况进行计算，例题参数如下：
筒体：Di=2000 mm，厚度δe = 20 mm。
接管：外径Do = 355.6 mm，厚度δt = 14.5 mm。
弹性模量：E = 195000 MPa。
内压：Pc = 1.2 MPa。
沿筒体、接管的壁厚划分3层，在接管附近的单元大小控制为25mm，两种软件的网格划分如下：
图1. ANSYS 软件中的单元划分
图2. ABAQUS软件中的单元划分
3.1 ANSYS软件计算结果
1）Solid 185 单元（8节点实体单元）
a. Solid 185 (KEYOPT = 0，使用B方法全积分（默认选项）。
图3. ANSYS Solid 185 单元（KEYOPT = 0）--应力强度分布云图
b. Solid 185 (KEYOPT=1，带沙漏控制的一致缩减积分)。
图4. ANSYS Solid 185 单元（KEYOPT = 1）--应力强度分布云图
c. Solid 185 (KEYOPT=2，增强应变公式)。
图5. ANSYS Solid 185 单元（KEYOPT = 2）--应力强度分布云图
d. Solid 185 (KEYOPT=3，简单增强应变公式)。
图6. ANSYS Solid 185 单元（KEYOPT = 3）--应力强度分布云图
2）Solid 186 单元（20节点实体单元）
a. Solid 186 (KEYOPT = 0，一致缩减积分（默认选项）。
图7. ANSYS Solid 186 单元（KEYOPT = 0）--应力强度分布云图
b. Solid 186 (KEYOPT = 1，全积分)。
图8. ANSYS Solid 186 单元（KEYOPT = 1）--应力强度分布云图
从以上分析结果可看出，在ANSYS中，采用不同的单元类型或不同的单元选项，在网格划分密度一定的情况下，结构上的最大应力强度值是不同的。
3.2 ABAQUS软件计算结果
1）C3D8R单元（8节点缩减积分实体单元）
图9. ABAQUS C3D8R 单元--应力强度分布云图
2）C3D8单元（8节点全积分实体单元）
图10. ABAUQS C3D8单元--应力强度分布云图
3）C3D20R单元（20节点缩减积分实体单元）
图11. ABAQUS C3D20R单元--应力强度分布云图
4）C3D20单元（20节点全积分实体单元）
图12. ABAQUS C3D20单元--应力强度分布云图
从以上分析结果可看出，在ABAQUS中，采用不同的单元类型或不同的单元选项，在网格划分密度一定的情况下，结构上的最大应力强度值也是不同的。详见表一。
3.3 结果比较
表1. 内压作用下的计算结果
通过对表1所列结果的纵向和横向对比可以发现：
1）在ANSYS中，对8节点单元，B方法全积分和缩减积分的计算结果基本相同；增强应变公式和简单增强应变公式的计算结果基本相同，且后两种单元选项比前两种的计算结果大9%左右；对20节点单元，两种单元选项的计算结果基本相同，与8节点单元下的后两种单元选项的计算结果很接近，略大1%；
2）在ABAQUS中，对8节点单元，两种单元选项的计算结果比较接近；对20节点单元，两种单元选项的计算结果基本相同，比8节点单元下的两种单元选项的计算结果大7%左右；
3）对于8节点的实体单元，ANSYS提供了四种单元选项，当选用默认选项时，求解速度最快；与ABAQUS的C3D8R单元的计算结果很接近；
4）当采用20节点实体单元时，两种软件的计算结果是最为接近的，相对误差在1%以内。
4 结论及后续研究
从以上的实例及分析比较，可以得出以下的结论：
1）压力容器行业中，对受压元件进行应力分析时，采用ANSYS或ABAQUS软件时，单元类型的选择、单元选项的设置均是非常重要的；
2）相比较而言，ANSYS中采用Solid185(8节点实体单元—缩减积分)与ABAQUS中采用C3D8（全积分）结果最为相近。推荐在求解精度与速度折衷的情况下，优先选这两种单元。
由于篇幅限制，本文仅对内压作用下的柱壳径向开孔进行了研究，但压力容器的开孔接管大多还会承受管口载荷的作用，后续论文会对本文算例在管口载荷作用下的计算作进一步研究，还会就ANSYS和ABAQUS软件对网格疏密度的敏感性进行研究。
