跟着带孔平板学习APDL（上）
ANSYS参数化语言（APDL）是一种用来完成有限元常规分析操作或通过参数化变量方式建立分析模型的脚本语言。APDL扩展了传统有限元分析范围之外的能力，并扩展了更高级的运算，包括灵敏度研究、零件库参数化建模，设计修改和优化设计等。
作为弹性力学中的一个经典问题，带孔平板的平面应力分析一直是各类有限元教程如《xxx有限元分析由入门到精通》、《xxx软件有限元分析初级应用》中的常客，也是各大院校相关课程大小作业的热点。
公众号由Workbench小学生升级为CAE中学生，谨以此文作为我的Workbench小学生涯的毕业论文，有什么不对的欢迎各路大神指正。
1. 带孔平板的基本分析流程
如图一带孔薄板，长200mm，宽100mm，中心开孔直径为10mm，薄板厚度为1mm，左右两侧边线受到100mpa的载荷作用，方向相反。薄板弹性模量为21000mpa，泊松比为0.3。
1.1基本设置
开始一个分析之前，根据实际情况，一般需要结束上一个分析，清除内存数据，设定分析类型和名称等，以此作为一个命令流的开头。
FINISH:表示结束上一个分析；
CLEAR:清除数据；
PREP7：进入前处理器，前处理器的主要的工作包括了设定单元类型，建立模型，设置材料参数，划分网格和设置边界条件；
TITLE,plate：设定分析名称为plate。
1.2单元类型与材料设置
在个人没有配备超级计算机之前，单元类型的选择都还是有限元分析当中值得重视的一个问题，理论上说，选择实体单元可以解决所有的问题，但是，同时采用实体单元的计算规模也是各种不同类型的单元当中最大的。因此，在符合相关的力学假设的条件下，应该优先选用维数低的单元，梁杆>板壳>实体。
ET:定义单元类型，1表示单元类型编号，当采用多种单元时单元类型编号依次增加，PLANE183表示平面单元，183为单元编号，以上命令等价于ET,1，183 。
KEYOPT:设置单元选项，PLANE183单元既可以作为平面单元，也可以作为轴对称单元，设置KEYOPT(3)=3，即为考虑厚度的平面应力模型，参考如下：
R:设置单元实常数，通常用于设置平面厚度，杆单元的截面，当采用PLANE183单元时，仅当KEYOPT(3)=3可以设置。
MP:设定材料参数，EX为弹性模量，PRXY为泊松比，1为材料编号，同样，当有多种材料时编号依次增加，由于ANSYS APDL中没有具体的单位制，所以在输入参数前需要自行规划，所有参数使用统一的单位制，否则计算时可能会出现错误。
1.3建模
在ANSYS APDL中主要有直接创建有限元模型和由几何模型通过网格划分转化为有限元模型两种建模方式，除了杆梁单元，常见的模型通常都是用第二种方式创建模型，对于一些复杂的模型，也可以通过第三方的建模软件进行建模，然后导入到ANSYS中。
BLC4：通过角点创建矩形面（长方体）,以上命令表示以坐标点（-100，-50）为起点，创建长度为200（x方向），高度为100（y方向）的矩形面，在此命令基础上增加一个z方向的参数，即可创建长方体，例如：BLC4，-100，-50，200，100，100。
CYL4：通过圆心和半径创建圆面（圆柱体）。
ASBA:布尔减运算，从一个面中减去另外一个面，1和2分别为面编号，通用命令格式为XSBY，X为被减几何体，其中X、Y可以为L、A和V，分别对应了线、面和体，需要注意的是只能用更高维数的几何体去减低维数的几何体，因此X的维数应该低于或者等于Y的维数，另外，同样可以使用工作平面去切割几何体，需要注意的是切割平面仅限于XY平面，因此需要配合坐标系的平移和旋转使用，命令为LSBW,ASBW和VSBW。
1.4网格划分
ESIZE：定义全局单元大小。
1.5边界施加
SFL：在线上施加分布载荷，2为线编号，通用格式为SFX，SF为在节点组处施加均布载荷，当X为BEAM、E、L和A分别为在梁、单元、线和面上施加均布载荷；
DL：约束线自由度，通用格式为DX，D为约束节点自由度，当X为K、L、A和SYM时，分别为约束关键点、线、面和对称面自由度。
1.6结果输出
PLNSOL：以云图的形式显示节点计算结果，类似地，PLESOL为以云图的形式显示单元计算结果。此处为显示x方向应力结果云图。
1.7小结
以上为带孔平板的一个完整的分析流程，在编写命令流时需要注意的就是线和面的编号不能混淆，否则会出现错误的分析结果，因此，在调试命令流时要经常查看目标的编号。
在边界施加处，注意到施加了两个约束来替代相反分析的作用力，使整个分析可以顺利进行。如果在两边同时施加方向相反的作用力，则有可能会因为网格不对称导致刚体位移。
如果把L4自由度完全约束，L1自由，则应力结果如下所示：
由图可以看到，虽然整体最大应力只是发生了微小变化，但是在结果图示上面可以明显看到应力的不对称，靠近左侧固定端有部分区域明显偏大，这里住着赫赫有名的圣维南老爷子，为了纪念他，该区的名字叫圣维南区。
2. 对称分析
当结构在几何、约束、载荷和材料上都具有对称性的时候，可以对结构的一部分进行分析，减小计算规模。明显，本案例关于x轴和y轴均对称，因此我们可以采用1/4模型对它进行分析。
2.1直接使用对称
类似地，以坐标原点为起点，创建一个100*50的矩形，然后通过布尔运算，减去圆心在原点，半径为5的圆，创建1/4模型。
下面介绍另外一种建模方式，通过创建关键点，然后由关键点连接成线，最后通过选取封闭线的方式完成平面的建模。
2.1.1 设置建模参数
APDL中可以直接在任何的位置定义参数，然后通过引用该参数进行建模或者其他操作。参数类型可以为数值型和字符型，可以通过在程序中赋值或者通过交互结界面输入，定义参数时需要注意以下几点：
以字母开头，长度32个字符以内；
不能用宏命令专用的局部参数名；
不能使用ANSYS的标识和已定义的组件名称，例如temp，ux，all等；
简单易记，不要非主流。
2.1.2 创建点
K:创建关键点命令，第一行命令完整格式为K,1,0,R,0，1代表关键点编号，0，R，0分别为xyz轴坐标值，平面模型，z=0时忽略，R值即为前面定义的半径5，在APDL的命令中同样支持参数的数学运算，如命令“K,3,L/2,0”等价于“K,3,100,0”。
KX(3)：参数直接引用，表示该参数数值为3号关键点的x轴坐标，对于节点，同样适用，将K改为N即可。
2.1.3 创建线
LARC：创建圆弧指令，通过起点（关键点1），终点（关键点2），曲率中心上任意点（关键点6，可以不是圆心），半径（R）创建圆弧线。
L:通过连接关键点创建线，跟坐标系有关，可以为直线或者曲线。
2.1.4 创建面
AL：通过依次连接收尾闭合的一组线创建平面，后面数字为线编号；
KDELE：删除关键点，通用格式为XDELE，当X为K、L、A和V时，分别表示了关键点、线、面和体的删除操作。当删除对象包含有更低维度的元素时，可以选择是否保留低纬度元素。相反地，不能独立删除依附于高维度元素下的低纬度元素，例如，不能在保留线的条件下删除线的端点，同时，也不能删除已划分网格的实体模型。
2.1.5 边界施加
载荷施加命令中的P，即为前面参数设置时“P=-100”,采用参数输入的方法，可以方便快捷地修改一系列具有数学关系的建模尺寸或者载荷，减少参数遗漏等人为失误。
DL:具体详见前面描述，DL,2,1,SYMM具体含义为1号面上的2号线施加对称约束。2.1.6 结果输出
整体网格尺寸设置为2，同样，输出x方向应力如下：
2.2使用约束等效
对于关于线或者是平面的对称的结构，例如本案例，根据生活经验，在相同拉力下，平板中心位置必定保持原位，不会偏向其中一边，所以在数学角度来说，对称轴（面）上的点不会产生沿其法线方向上的移动，因此我们可以根据这一特点对结构进行等效约束。
整体网格尺寸设置为2，同样，输出x方向应力如下：
2.3小结
当其余输入条件一样时两种约束方式的计算结果完全一样，所以第二种方法可以作为对称约束的等效输入方式，输入以下命令，可以将对称模型进行全模型扩展显示：
3. 网格收敛性研究
对比前面两个案例，可以看到利用全模型分析的x分析最大正应力为311.789MPa，对称模型的x方向最大正应力为298.513MPa，而且默认的网格大小均为2，难道全模型和对称模型的分析结果差别是这么大吗？其实，细心对比两个结构的网格，可以看到，即便是默认网格尺寸一样，但是由于结构差异，两者划分的网格还是不一样的，全模型的网格不会自动划分得完全和1/4模型扩展显示一样（手动不知道，没有那个能力），否则案例1施加两边一样的载荷时就不会出现刚体位移了。
结果对比是大部分新手很热衷的一件事，往往会因为一两个结果的差异而在网上与别人争个面红耳赤，充分发扬了传统技术交流的喷倒为止。赢者沾沾自喜，败者怀疑自我，怀疑软件，甚至上升到怀疑有限元理论。其实在我看来，很多的对比都有关公战秦琼之嫌，说到打乒乓球，现在的我还真不怕3岁的刘国梁。
既然要对比，就必须要在同一个起跑线上（边界一样），对比的基准必须是两者的收敛解（结果都通过网格无关性验证），这就是我们下面要讨论的问题。
在没有应力奇异的前提下，通过不断的网格细化可以得到网格无关解，但是盲目地细化整体网格会造成计算量的急剧增加，或许超出了电脑的极限也得不到一个收敛解，因此我们需要有针对性地进行局部网格细化。
3.1全模型
LESIZE：为选择的线设置网格尺寸，5为线编号，0.1为线上单元尺寸，该命令后仍有许多选项可以设置，具体查看帮助文件，尺寸0.1为多次测试后的结果，算得结果如下：
3.2对称模型
LREFINE：局部单元重新划分命令，前面三个1依次为起始线编号，最终线编号，和递增间隔，5为细化指数，默认值为1，最大值为5，大概细化程度如下：
算得收敛值如下：
3.3小结
对比以上结果，当结果收敛后，全模型算的结果与1/4模型计算结果误差为0.3%，这部分差异来源于网格的差异，而并非边界条件和载荷的影响。对于简单的案例，看不出来对称结构分析的优势，当模型变得复杂后，合理利用对称分析可以极大地提高计算效率。
4. 不同单元类型的收敛性比较
前面案例主要讨论了网格数量对计算精度的影响，下面通过四个简单的案例来探讨一下提取不同单元类型时网格的收敛性研究，将前面1/4模型算得的应力值303.44MPa作为参考标准。为减少计算量，以下案例均采用1/4模型进行分析。选取单元类型分别为PLANE183（划分为四边形和三角形）和PLANE182（划分为四边形和三角形）。
PLANE183：
PLANE182:
4.1高阶四边形
高阶四边形单元为PLANE183的默认设置，结构形状不复杂时一般无需特别设置:
MSHPE：网格形状控制命令，0表示划分为四边形（六面体），若为1，则表示划分为三角形（四面体），2D表示单元维数，若为实体模型，则2D变为3D。
多次修改网格尺寸，计算结果如下：
4.2高阶三角形
同样地，改变网格形状，多次修改网格尺寸，计算结果如下：
4.3低阶四边形
直接采用PLANE182单元，单元形状采用默认设置，多次修改网格尺寸，计算结果如下：
4.4低阶三角形
参考高阶三角形的设置命令，多次修改网格尺寸，计算结果如下：
4.5小结
综合以上计算结果，对于本案例，在应力结果收敛性方面，高阶单元要优于低阶单元，即便是低阶四边形单元，在收敛性方面也略次于高阶三角形单元，四边形单元要优于三角形单元，而低阶的三角形单元的计算精度是最差的。本次的分析的结果对其他的静态结构分析在单元选取方面有一定的参考作用。
特别地，对于PLANE183，可以MSHMID,KEY命令对中间节点进行控制，KEY=0时为默认设置，边界区域单元中间节点与边界线或者边界面曲率保持一致；KEY=2时，设置所有中间节点位于单元直边上；KEY=3时，不生成中间节点。有兴趣的朋友可以比较KEY=3时PLANE183单元与PLANE182单元的收敛性。
5. 应力集中系数计算
应力集中是指应力在固体局部区域内显著增高的现象，多出现于尖角、孔洞、缺口、沟槽等区域，反映局部应力增高程度的参数为应力集中系数k，它是峰值应力和不考虑应力集中时的应力（名义应力）的比值，与载荷大小无关。
5.1数组定义与循环
之前第二节中也简单介绍了变量参数的定义，它只能存储一个参数值，而数组参数是按照多个行、列与面的结构存储多个参数值，包含多个元素。
*DIM,PM,,7：定义一个名称为PM的一维数组，包含7个元素，可以通过直接赋值的方式确定数组中各个元素的值，例如：PM(1)=(1，2，3，4，5，6，7)，对于有特定规律的数组，我们可以通过循环语句实现数组元素赋值的操作，例如：
*Do,i,1,7
PM(i)=P*100/(D(i)-100)
*ENDDO
即实现了不同直径时带孔平板的平均应力数组的定义，i为循环变量，1为循环初始值，7为循环最终值，默认递增值为1，使用循环语句时必须要做到有始有终，一个do必定对应一个enddo，以下为实现不同直径建模的循环命令，表示创建了直径为D(1)=(5，10，15，20，25，30，35)系列的带孔平板，并且划分网格（网格无关性方面尚未验证，仅仅按照经验粗略划分）。
5.2结果提取
边界条件施加详见前文描述，不再赘述，进入结果提取。
NSORT,S,X：表示对X方向正应力（节点值）进行排序，同样命令包含了ESORT（单元值）的应用，具体应用详见help文件。
*GET,PK(i),SORT,,MAX：提取排序后的最大正应力，并将结果存储到数组PK当中，关于*GET用法只能说功能强大，几乎无所不能，具体详见help。
ACLEAR,1：清除所有的网格。
ADELE,1,,,1：清除所有的实体
完成特定半径的最大正应力计算并且赋值后，删除所有的网格和实体，结束循环，进入下一个直径的计算。
5.3结果绘制
SMOOTH,D(1),k(1),,3,x(1),y(1),2：结果绘制命令，横坐标为直径，纵坐标为应力集中系数，绘图命令具体含义详见help。
5.4结果的路径映射
限元分析中最重要的是结果的后处理，将分析结果通过通俗易懂的数据、图片等展示给相关人士是分析人员的必备技能之一，大逆不道地说一句，分析可以做不好，但是PPT一定要炫。
在本案例当中，比较关心的是X方向应力沿纵向截面的变化情况，即通过图片显示“应力在固体局部区域内显著增高”的现象。下面以直径为10mm举例说明结果的路径映射的操作方法，具体前处理方面参考前文。
PATH,PATH1,2：定义路径命令，PATH1为路径名称，2为关键点数量；
PPATH,1,,0,5：定义路径第一个端点，完整命令为“PPATH,POINT,NODE,X,Y,Z,CS”，POINT为路径点的编号，NODE为节点编号，若此处为空，则采取后面的坐标进行定义，CS为输入坐标系，默认为当前坐标系，需要注意的是，当需要进行应力线性化操作时则必须要采取按节点的定义方式；
PDEF,XS,S,X：将X方向应力在路径上进行插值计算，并且命名为XS;
PLPAGM,XS,200,'NODE'：沿路径显示X方向应力，放大200倍，同时显示节点。
程序运行完成后结果如下:
5.5小结
本小节主要内容是数组应用和do循环语句的应用，类似地，可以将参考绘制带孔平板不同直径的应力集中系数的方法，绘制不同网格尺寸下的应力变化曲线进行网格无关性研究。