如何在ABAQUS中构建奇异单元
1. 奇异单元简介
在利用ABAQUS计算裂纹体的应力强度因子时，为了满足裂尖奇异性的要求，通常会在裂纹尖端附近的区域采用奇异单元进行划分。对于平面单元，通过将四边形单元的其中一条单元边折叠，并调节二阶四边形单元中间节点的位置，可以获得满足不同奇异性要求的奇异单元。
（1）将二阶四边形单元的其中一条单元边折叠，同时将折叠边的三个节点绑定到一起，同时将中间节点移动到单元边的1/4位置处，则可以在裂纹尖端获得r-1/2的奇异性：
为了将折叠边的三个节点绑定到一起，可以对三个节点使用线性约束方程或者多点约束，一种更为简便的方法为对折叠边的三个节点使用相同的节点号（需要通过修改inp文件实现），如图1.1所示。
图1.1 奇异单元的构建（r-1/2奇异性）
（2）如果将四边形单元的其中一条单元边折叠，但折叠边的三个节点可以自由移动，并且中间节点仍然保持在单元边的中间位置，如图1.2所示，则可在裂纹尖端获得r-1的奇异性：
图1.2 奇异单元的构建（r-1奇异性）
（3）如果将四边形单元的其中一条单元边折叠，但压缩边的三个节点可以自由移动，并且将中间节点移动到1/4的位置，如图1.3所示，则可在裂纹尖端获得上述奇异性的组合：
图1.3 奇异单元的构建（组合奇异性）
由弹性理论可知，裂纹尖端附近的应力场可统一表示为如下形式：
其中σij为裂尖附近某点处的应力分量，K为应力强度因子，r为该点到裂纹尖端的距离，θ为该点与裂尖构成的连线与裂纹面所呈的夹角，可以看到裂纹尖端的应力分布具有r-1的奇异性，因此预计通过第一种方法构造的奇异单元计算得到的应力强度因子相比于其它两种方法是最为精确的。
在二维中心裂纹板的应力强度因子分析中，笔者已经详细介绍了具有不同奇异性的奇异单元的构建方法。由于在前文中，奇异单元是利用ABAQUS GUI中构建的，因此尽管裂纹尖端采用三角形单元划分，但ABAQUS会自动通过调整单元节点完成奇异单元的构建。而如果裂纹网格是来自外部导入的孤立网格（Orphan mesh），则ABAQUS会将裂尖单元识别为三角形单元，而非奇异单元（单元边被折叠的四边形单元）。尽管在裂纹尖端采用三角形单元也能完成裂纹应力场的计算，但是由于围线积分区域内不能使用三角形单元，因此无法采用围线积分法完成后续应力强度因子的评估。
为此，本文以一个含有环向穿透裂纹的圆柱体为例，详细介绍如何利用Hypermesh完成裂纹体网格划分以及裂尖奇异单元的构建，并评估裂纹尖端的应力强度因子。
2. 计算实例
假设有一长度为800mm，直径为250mm的圆柱体，在距离圆柱体端部100mm的位置处存在一条环形斜裂纹，该裂纹深度为4mm，裂纹面与圆柱体表面的倾角为45°，如图2.1所示。假定圆柱体远端承受有100Mpa的均匀拉伸应力，试评估裂纹尖端的应力强度因子。
图2.1 裂纹体示意图
2.1 裂纹体网格的构建
考虑到裂纹体的几何和边界条件均具有对称性，因此本文首先将其简化为轴对称模型，并计算应力强度因子，随后，在该模型的基础上，讨论如何在三维模型中构建奇异单元并计算应力强度因子。
对于二维轴对称模型，裂纹尖端的奇异单元为其中一条边被折叠的四边形单元，因此可以在Hypermesh中在裂纹尖端切分出一块圆形区域，并在该区域采用三角形单元进行划分，为了获得精确的应力强度因子，裂纹尖端的奇异单元不能少于16个（单元边夹角22.5°）。随后，在三角形单元外围采用多层环状四边形单元进行过渡，如图2.2所示。
图2.2 裂尖附近区域的网格划分
图2.2中的环状区域为围线积分区域，ABAQUS将以裂纹尖端的三角形单元的单元边作为第一层围道，后面的围道依次类推。需要注意的是，在围线积分区域不能出现三角形单元（裂尖的奇异单元为被折叠的四边形单元），围线积分区域的范围取决于用户定义的围道数量，例如，在图2.2中左上角存在一个三角形单元，则定义围道数量时，应保证围道积分区域不会将该三角形包含在内。
对于远离裂纹尖端的区域，可逐渐采用粗糙的网格进行过渡，最终建立的裂纹体网格如图2.3所示。
图2.3 裂纹体网格
采用Hypermesh划分裂纹体网格时，裂纹面两边的节点是相互连接的，因此裂纹面将无法 正常张开。ABAQUS提供了Seam功能，用于将裂纹面上某一侧的单元的节点赋予新的节点号，来实现裂纹面的分离，而在Hypermesh中，则需要用户通过Detach功能手动将裂纹面上的网格分离（注意不能将裂纹尖端的节点分离），最终分离完成后的裂纹网格如图2.4所示。
图2.4 裂纹面节点分离
Hypermesh的Edge功能可以用于检查网格模型中存在的自由边，如图2.4中的红色边线所示，可以看到裂纹面分离之后应形成自由边。在Hypermesh中完成网格划分后即可导出inp文件。
2.2 裂尖奇异单元的构建
通过Hypermesh划分得到的裂尖单元为三角形单元，并非单元边被折叠的奇异单元，因此无法直接导入到ABAQUS中用于应力强度因子计算，为了构建奇异单元，需要对裂纹单元的节点进行重新编号。打开导出的inp文件，将其定位到单元定义的所在行，如图2.5所示。
图2.5 单元节点编号定义
从图2.5中可以看出，该inp文件中定义了两种单元，一种为CPS3单元，即3节点的平面应力单元，另一种为CPS4R单元，即带有减缩积分的4节点平面应力单元，通过在Hypermesh中查询裂尖三角形单元的编号可知，裂尖单元恰好为图2.5中1到16的单元，这些单元具有共同的节点号2，显然该节点为裂纹尖端的节点。参考节1所述方法，为了获得单元边被折叠的四边形单元，只需将这些单元的节点号定义中复 制一个节点号2，随后将其移动到CPS4R单元的定义中即可。例如，对于裂尖单元1，其节点编号为：
2, 3, 1
则对应的奇异单元的节点编号应为：
2, 2, 3, 1
需要注意的是，由于Hypemesh默认导出的单元为一阶单元，因此不存在中间节点，如果导出的inp文件中的单元为二阶单元，则应注意单元编号顺序。在ABAQUS中，单元节点是沿逆时针方向进行编号的，并且先编号边角上的节点，随后编号中间节点，例如，对于图1.1中未折叠的四边形单元，其节点编号应为：
1,2,3,4,5,6,7,8
因此与之对应的奇异单元的节点编号应为：
1,2,3,3,5,6,3,8
如果没有对单元节点进行正确编号，则导入ABAQUS中将得到畸形的单元。对于图2.5中所示的单元，其对应的奇异单元编号如图2.6所示。
图2.6 奇异单元节点编号
将修改之后的inp文件导入ABAQUS，利用Tool→Query→Element查询裂尖节点编号，查询结果如图2.7所示。
图2.7 裂尖单元查询结果
从图2.8可以看出尽管裂纹单元为三角形，但查询节点显示这些单元为4节点的四边形单元，且具有重复的节点号2，证明裂尖单元的一条单元边被折叠。
为了在裂尖获得r-1/2奇异性，需要使用二阶四边形单元，并将中间节点移动到1/4位置。在Mesh模块中将单元类型更改为本次分析需采用的二阶轴对称单元（CAX8R），注意，对于非围线积分区域存在的少量三角形单元，为避免单元协调问题，应使用CAX6单元而非CAX6M单元。随后，利用Mesh→ Edit→ Node→ Adjust midside功能调整奇异单元的中间节点，选择奇异单元与裂尖连接的单元边，将Mideside node parameter更改为0.25，使得裂尖单元的中间节点移动到1/4位置处，如图2.8所示。
图2.8 移动中间节点
移动完成后利用Mesh→ Verify→ Analysis Check检查是否存在错误单元，如果有些非奇异单元的中间节点也发生了移动并且形成错误单元，可以利用Mesh→ Edit→ Node→ Drag手动将这些节点调整到中间位置。通过上述步骤，即可获得裂尖具有r-1/2奇异性的轴对称单元。
完成上述步骤后即可在Interaction模块中创建裂纹，并在Step模块的时间历程输出中设置输出围线积分，具体步骤可参见二维中心裂纹板的应力强度因子分析，这里不再赘述。图2.9给出了计算得到的裂纹尖端的Von Mises应力分布。
图2.9 裂纹尖端Von Mises应力分布
从图2.9可以看出裂纹尖端的应力云图呈现出蝴蝶形状，这与通过弹性理论得到的裂纹应力分布是比较吻合的。
分别采用三种单元：一阶单元（仅单元边折叠）、二阶单元（折叠单元边，移动中间节点和不移动中间节点）计算裂纹尖端的应力强度因子，利用下式对应力强度因子进行归一化处理：
其中K为应力强度因子，σ为圆柱体远端拉伸应力，a为裂纹长度，归一化后的I型和II型应力强度因子随围道的变化规律如图2.10和2.11所示。
图2.10 I型应力强度因子随围道变化（轴对称单元）
图2.11 II型应力强度因子随围道变化（轴对称单元）
从图2.10和图2.11中可以看出，随着围道数量的增加，三种单元计算得到的应力强度因子均呈现收敛的趋势。在一阶单元和不移动中间节点的二阶单元中，由于在第一个围道上围线积分的路径无关性无法得到严格满足，因此差距较大。而移动中间节点的二阶单元由于具有r-1/2奇异性，因此计算结果精度最高。当然，由于在裂纹尖端采用了比较细密的网格，因此三种单元都给出了非常精确的应力强度因子。并且，即使采用较为粗糙的网格，只要围线积分的路径无关性得到满足，也能获得较为准确的应力强度因子。
2.3 三维实体模型的应力强度因子计算
与二维奇异单元的构建方法类似，三维奇异单元同样可以通过修改inp文件进行创建。但是，当裂纹尖端的单元数量较多时，通过修改inp文件来创建奇异单元是比较繁琐的。对于本文计算分析的环向裂纹，一种比较简便的方法为将前面二维轴对称分析生成的inp文件（已经完成奇异单元的创建以及中间节点节点而的移动）重新导入Hypermesh，随后在Hypermesh通过扫略或拉伸来创建实体网格，经笔者验证发现，这样生成的实体网格可以保证裂纹尖端为三维奇异单元，无需在inp文件中修改节点编号。
利用上述方法，创建单元类型为C3D20R的奇异单元，计算得到裂纹尖端的Von Mises应力分布如图2.12所示。
图2.12 三维实体模型裂纹尖端Von Mises应力分布
图2.13和图2.14给出了轴对称模型和实体模型计算得到的应力强度因子对比。
图2.13 I型应力强度因子随围道变化（轴对称和实体单元对比）
图2.14 II型应力强度因子随围道变化（轴对称和实体单元对比）
从图2.13和图2.14可以看出，两种模型计算得到的应力强度因子非常吻合，这证明对于此类可以简化为二维问题的模型，即使采用二维奇异单元进行计算也能获得非常精确的应力强度因子。
3. 总结
本文简单介绍了奇异单元中不同奇异性的实现方法，并针对ABAQUS中孤立网格无法直接创建奇异单元的问题给出了相应的解决方法。随后，在二维奇异单元的基础上，提出了一种穿透裂纹中三维奇异单元的创建方法，利用这些方法可以精确计算裂纹尖端的应力强度因子。
