P-方法的ANSYS实现

在有限元分析中，当对模型的整体或局部计算精度有更高要求的时候，总是绕不过两个基本方法：H-方法和P-方法。

其中，H-方法的思想，是通过增加节点、增加单元数量的方法，来对模型进行更精细的模拟，特点是形函数不变，单元数量增加；P-方法则相反，从单元自身做起，采用比较复杂的形函数，增加单元形函数的阶数 (Order Polynomial)而不增加单元数量和节点数量。H-方法适用于任何分析类型，但是P-方法只适用于线性结构精力分析。

P-方法最大的优点就是在精细化分析的过程中，不需要精细的网格，多次分析的时候不需要调整网格密度，减少工作量，提高效率，同时又不损伤计算精度。采用的形函数阶数越高，计算结果的精确程度也就越高，但对计算机的要求也就越高了。因此，一般情况下多项式函数的最高阶次P<9可以接受。P方法的收敛准则可以包括计算模型节点上的位移、应力和应变或整体的应变能。

在建模过程中的实现方式：

 (1)定义P-方法单元

ANSYS支持P-方法的单元包括：

2-D Quadrilateral (PLANE145)，2-D Triangle (PLANE146)，3-D Brick (SOLID147)，3-D Tetrahedron (SOLID148)，3-D Shell (SHELL150)

(2)控制P-方法单元阶数

上文说过，P-方法单元的一个重要特点就是可以控制形函数的阶数。ANSYS允许用户在模型的整体或局部范围上分别进行控制，例如：

PPRANGE,2,8

ET,1,PLANE145,3,5

前面一句，采用命令PPRANGE从整体定义P-单元阶数为2~8；后面一句，在用ET命令定义单元的时候，通过KEYOPT(1)和KEYOPT(2)对145单元进行设置——即模型中，145单元的阶数为3~5，其他单元阶数为2~8。如果不进行设定， ANSYS 默认的 p 值范围为 2~8，SHELL150单元默认的p值范围为3~8。

(3)P-单元的网格与模型设置

建议采用从上向下的建模顺序，避免通过直接生成节点的方式建模，因为P-单元需要的中间节点难以手动指定。可用/EFACET命令设置高阶大单元的细致形状。在边界等区域上，可以通过建立耦合方程进行位移的相对控制。尽量避免对单一节点施加位移约束或荷载，这样容易产生应力集中、奇异，若一定要采用，在需要在收敛性判定中排除这些节点。不要在 p 型单元的边上的中间节点上施加力或力矩。

(4)P-方法的收敛与求解

一般建议采用条件共轭梯度(PCG)算法来作为P-方法的求解器，原因是对于3D问题或大规模2D问题(通常达到 40,000 个自由度 以上)，PCG 算法比 JCG 及波前法要快得多，对于 4 阶或更高阶次的 p 型单元，采用 PCG 的快速计算具有非常重要的意义。对于几何形状不太好的单元或材料性质不连续的问题，采用 PCG 方法需要较多的迭代次数，可以通过命令<EQSLV,PCG>来进行设置。

在 ANSYS 中，默认约束边或面上的所有节点均采用整体坐标系，若将其都旋转到同一个节点坐标系中(采用命令NROTAT)，在节点施加位移约束的同时，也就对内部自由度进行了约束。

若只对模型的某一个局部区域的受力状态特别关注，可以使用局部收敛准则来对该区域进行计算精度的控制(采用命令PCONV)；对于具有奇异性或沿材料分界的节点，本身由于奇异性就很难收敛，一般不能指定为收敛控制的对象。对于这些区域，应该在收敛控制时进行区域的排除处理(采用命令PEXCLUDE)。

可以通过命令ESEL,,PINC来选择或指定 p 方法的计算模型中将包含的单元，使用命令ESEL,,PEXC来选择或指定 p 方法的计算模型中将排除的单元；通过命令PINCLUDE,STAT以及PEXCLUDE,STAT来显示采用 p 方法的计算模型中所包含单元以及所排除单元的状态。

(5)P-方法后处理

在后处理中，可以在 p 型单元内部划分较小的面片来进行显示(采用命令/EFACET来设置每一边的划分)，同时可以采用结果查询的方式来获取面片子网格上的计算结果；使用命令AVRES来对公共节点的结果进行平均化处理；采用命令PRCONV或PLCONV来显示所设置的收敛值。

在采用P-方法的注意事项：

(1)在同一个模型中，p 型单元和 h 单元不能同时被激活。

(2)避免单点载荷或约束，否则会导致应力奇异。

(3)不能在单元的中间节点施加力和矩。

(4)不要在有奇异的节点或沿着材料分界面的节点处规定应力或应变的收敛准则，这些位置由于奇异性较难收敛。

(5)所有的 p 型单元要求中间节点要包括在几何定义中。

(6)对于p型单元，在中间节点和角节点同属于一个耦合对时，只有在角节点上才可以定义主自由度。

实例

/PREP7                                 !进入前处理

/TITLE, p-Method

ANTYPE,STATIC                !静力结构分析

ET,1,PLANE145,3,6,3              !设置 1 号单元(p 型单元)

MP,EX,1,2e5      !设置 材料常数

MP,PRXY,1,0.2

R,1,0.5

RECTNG,0,100,0,60          !建立模型

wpoffs,50,30

PCIRC,10,0,0,360

ASBA,1,2

numcmp,all

SMRTSIZ,4

AMESH,1

NSEL,S,LOC,X,0                 !施加约束

D,ALL,UX,0

NSEL,S,LOC,Y,0

D,ALL,UY,0

NSEL,S,LOC,X,100             !施加均布荷载

SF,ALL,PRES,-200

ALLSEL

FINISH

/SOLU

P1=NODE(50,30+10,0)        !获取圆孔顶部节点编号赋给参数 P1

P2=NODE(50,30-10,0)            !获取圆孔底部节点编号赋给参数 P2

PMOPTS                              !设置 p 方法中固定单元 p 阶次的误差，默认 5%

PCONV,0.1,S,X,P1               !针对 P1 节点，采用 p 方法，设置应力 SX 的误差在 0.1%以内

PCONV,0.1,S,X,P2               !针对 P2 节点，采用 p 方法，设置应力 SX 的误差在 0.1%以内

SOLVE                                 !进行求解

 

/POST1                                 !进入后处理

SET,1                                    !调出结果的 1 号数据集

PLNSOL,S,X                        !图形显示 x 方向的应力计算结果

PLDISP,1                              !图形显示受力后的模型结构

*GET,P1_SX,NODE,P1,S,X !获取节点 P1 处的应力 SX，赋给参数 P1_SX

*GET,P2_SX,NODE,P2,S,X !获取节点 P2 处的应力 SX，赋给参数 P2_SX

PLCONV                                !图形显示 p 方法的收敛曲线

PPLOT                                    !图形显示实际所用单元 p 阶次

*STATUS                                  !列显参数的内容

FINISH                                 !结束

参考资料：曾攀《有限元分析基础教程》