ANSYS 静力 分析实例 悬臂梁

Blog Links

• DalNur | 博客总目录

• ANSYS 有限元分析 概述

• ANSYS 有限元分析 坐标系/工作平面

• ANSYS 有限元分析 几何建模

• ANSYS 有限元分析 网格划分

• ANSYS 有限元分析 选择与组件

• ANSYS 有限元分析 修改与编辑

• ANSYS 有限元分析 接触分析

• ANSYS 有限元分析 加载/求解/输出

• ANSYS 有限元分析 后处理 General Postproc

• ANSYS 有限元分析 后处理 结点解与单元解

• ANSYS 有限元分析 命令流 实例

• ANSYS 静力 分析实例 悬臂梁

• nCodeDL 高周疲劳 分析实例 悬臂梁

一、前言

  本文以力学中最简单的悬臂梁受载为例，介绍了采用 ANSYS 进行有限元分析的基本操作流程，并结合基础的有限元知识，探讨了有限元近似解的精度水平。

  本文还讨论了单元类型、单元阶次、积分方案、网格尺寸等对计算结果的影响。

  本文主要章节安排如下，第二部分介绍悬臂梁集中荷载作用下位移和应力的弹性力学上的理论解答，第三部分介绍相应的有限元分析流程及解答，第四部分介绍各种因素对有限元计算结果的影响及误差来源，第五、六部分简单说明了后处理中内力的提取方法。第七部分给出一些有限元分析的合理性建议。

二、弹性力学解答

  悬臂梁自由端受竖向集中荷载作用时，其挠曲线方程的弹性力学解答为：

$$y=-\frac{F}{6EI}{x}^3+\frac{Fl}{2EI}{x}^2+\frac{F{h}^2}{8GI}x\text{\hspace{1em}(1)}$$

式中，$h$ 为梁的高度；$G$ 为剪切模量，可由弹性模量及泊松比计算得到。

  悬臂梁自由端受竖向集中荷载作用时，任意截面处最大弯曲正应力的材料力学一致为：

$${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{F(l-x)}{W}\text{\hspace{1em}(2)}$$

式中，$W$ 为抗弯截面模量。

  悬臂梁的最大挠度出现在梁的自由端位置处，由由式 (1) 可计算得悬臂梁的最大挠度为：

$$y_{\max }=\frac{F{l}^3}{3EI}+\frac{Fl{h}^2}{8GI}\text{\hspace{1em}(3)}$$

  式 (3) 等号右侧第一项就是我们熟悉的悬臂梁在竖向集中荷载作用下最大挠度的材料力学解答。显然，第二项是剪切变形对梁挠度的贡献。由此可见，梁任意位置处的挠度由弯曲变形和剪切变形两部分构成，即：

$$挠度=弯曲变形+剪切变形\text{\hspace{1em}(4)}$$

  该梁任意位置处因剪切变形产生的挠度与因弯曲变形产生的挠度比值为：

$$\frac{Fl{h}^2/8GI}{F{l}^3/3EI}=\frac{3h^{2}E}{8l^{2}G}=\frac{3h^2}{8l^2}\cdot 2(1+\nu )=\frac{3h^2}{4l^2}\cdot (1+\nu )=\frac{3(1+\nu )}{4}\cdot \frac{h^2}{l^2}\text{\hspace{1em}(5)}$$

  当 $h\ll l$ 时，上式趋近于 0，这说明当梁足够细长时，因剪切变形产生的挠度可以忽略不计，即梁的挠度主要由弯曲变形引起，可见材料力学的解答具有一定的合理性，事实上，材料力学中所研究的梁就是要满足跨高比大于10。

三、有限元解答

  悬臂梁长度为 360 mm，其横截面尺寸为 H×B = 12 mm × 6 mm。材料为钢材，牌号为 Q235B，其弹性模量为 200 Gpa，泊松比为 0.3。其端部承受集中荷载 F = 500 N，如下图所示：

  本文以如上图所示的悬臂梁为例，分别采用梁单元、壳单元和实体单元建立其有限元模型，通过对比跨中位置处的挠度及固定端截面上的最大弯曲正应力，来说明不同建模方式对计算结果的影响。

  集中荷载施加在梁端部截面的形心位置处，因集中荷载作用处有很大的应力集中，附件的计算结果误差相对较大，根据圣维南原理，跨中位置离开集中荷载作用点已足够远。因此，我们采用跨中位置处的竖向位移为对比对象。

  在材料力学中，推导梁的弯曲正应力时，除了采用平截面假定外，还默认材料的泊松比为 0，即认为横截面两个方向上的变形相互独立互不影响。因此，为了能严格的和理论解答做对比，在有限元建模时，材料的泊松比设置为 0。

  事实上，弹性力学给出的该悬臂梁的位移及应力解答也并不是实际情况下最精确没有误差的解答，因为实际的应力场、位移场相当复杂，我们一旦合理的给出其假设形式，那么在此基础上得到的计算结果就已经偏离了理论上的精确解，只不过，此种解答误差足够小，我们能接受，就将其近似认为是理论上的精确解。

FINISH
/CLEAR
                      ! Units: mm, N
    L    = 360        ! Length
    H    = 12         ! Height
    B    = 6          ! Width
    E    = 200000     ! Young’s modulus
    nu   = 0          ! Poisson’s ratio
    NumL = 10         ! Number of elements in x direction (长度方向单元个数)
    NumB = 2          ! Number of elements in y direction (宽度方向单元个数)
    NumH = 2          ! Number of elements in z direction (高度方向单元个数)
    F    = 500        ! Point load

3.1 梁单元建模

1. 前处理

/PREP7

    K, 1, 0, 0, 0      ! 梁轴线起点
    K, 2, L, 0, 0      ! 梁轴线终点 
    K, 3, 0, 20, 0     ! 方向关键点

    LSTR,1,2           ! 创建直线

    /VIEW,, 1, 1, 1    ! 设置视图
    LPLOT              ! Plot Lines

    ET,1,BEAM188       ! 188一阶一次形函数单元
    KEYOPT,1,3,0

    ET,2,BEAM188       ! 188一阶二次形函数单元
    KEYOPT,2,3,2

    ET,3,BEAM188       ! 188一阶三次形函数单元
    KEYOPT,3,3,3

    ET,4,BEAM189

    MP, EX, 1, E
    MP, NUXY, 1, nu

    SECTYPE,1,BEAM,RECT
    SECOFFSET,CENT
    SECDATA,B,H,NumB,NumH

    LSEL,S,LINE, ,1                ! 选择1号线
    LATT,1,,2,,3,3,1               ! 设置线属性，2为单元类型编号。

    LESIZE,1, , ,NumL,1, , , ,0    ! 布种子
    LMESH,ALL                      ! 执行网格划分

FINISH

2. 加载与求解

/SOLU

    ! Fixed end
    NSEL, S, LOC, X, 0
    D, ALL, ALL, 0

    ALLSEL,ALL                           ! 全显示

    ! Loads at free end
    N1 = NODE(L, 0, 0)
    F, N1, FY, -F

    SOLVE

FINISH

悬臂梁梁单元有限元模型

3. 后处理

/POST1

    /VIEW,, 0, 0, 1
    /ESHAPE,1                            ! 显示梁截面
    PLDISP, 2
    PLNSOL, S, X

    N2 = NODE(L/2, 0, 0)                 ! 获取跨中结点编号
    NSEL,S,NODE, ,N2                     ! 选择结点
    PRNSOL,U,Y                           ! 列出结点位移

FINISH

表 1 实体单元有限元分析

单元名称	单元阶次	Shape functions along the length	跨中竖向位移 (mm)	最大拉伸弯曲正应力 (Mpa)
BEAM 188	一阶单元	Linear	14.029	1187.5
BEAM 188	一阶单元	Quadratic	14.077	1250
BEAM 188	一阶单元	Cubic	14.077	1250
BEAM 189	二阶单元	—	14.077	1250

3.2 壳单元建模

1. 前处理

/PREP7

    BLC4,0,-H/2,L,H 

    /VIEW,, 1, 1, 1    ! 设置视图
    APLOT              ! Plot Areas

    ET,1,SHELL181      ! 181一阶减缩积分单元
    KEYOPT,1,3,0

    ET,2,SHELL181      ! 181一阶完全积分单元
    KEYOPT,2,3,2

    ET,3,SHELL281      ! 281二阶单元

    MP, EX, 1, E
    MP, NUXY, 1, nu

    TYPE, 3                              ! 设置单元类型
    MAT, 1

    SECT,1,SHELL 
    SECDATA, B,1,0.0,3 
    SECOFFSET,MID  

    LSEL,S,LENGTH,,H                     ! 按高度选线
    LESIZE,ALL, , ,NumH, , , , ,0        ! 布种子
    ALLSEL,ALL

    LSEL,S,LENGTH,,L                     ! 按长度选线
    LESIZE,ALL, , ,NumL, , , , ,0        ! 布种子
    ALLSEL,ALL

    MSHAPE,0,2D                           ! 指定划分单元的形状，生成四边形单元。
    MSHKEY,1                              ! 采用映射网格划分方式划分网格
    AMESH,ALL                             ! 将所有面划分网格

    /REPLOT

FINISH

2. 加载与求解

/SOLU

    ! Fixed end
    NSEL, S, LOC, X, 0
    D, ALL, ALL, 0

    ALLSEL,ALL                           ! 全显示

    ! Loads at free end
    N1 = NODE(L, 0, 0)
    F, N1, FY, -F

    SOLVE

FINISH

3. 后处理

/POST1

    /VIEW,, 0, 0, 1
    /ESHAPE,1                            ! 显示壳截面

    PLDISP, 2
    PLNSOL, S, X

    N2 = NODE(L/2, 0, 0)                 ! 获取跨中结点编号
    NSEL,S,NODE, ,N2                     ! 选择结点
    NSEL,A,NODE, ,N1                     ! 选择结点
    PRNSOL,U,Y                           ! 列出结点位移

FINISH

表 2.1 壳中面位于xoy平面内 宽度为壳的厚度 荷载位于壳中面内

单元名称	单元阶次	积分方案	跨中竖向位移 (mm)	最大拉伸弯曲正应力 (Mpa)
SHELL181	一阶单元	减缩积分	18.220	771.898
SHELL181	一阶单元	完全积分	13.925	1179.58
SHELL 281	二阶单元	—	13.023	1159.75

表 2.2 壳中面位于xoz平面内 高度为壳的厚度 荷载垂直于壳中面

单元名称	单元阶次	积分方案	跨中竖向位移 (mm)	最大拉伸弯曲正应力 (Mpa)
SHELL181	一阶单元	减缩积分	14.021	1187.50
SHELL181	一阶单元	完全积分	14.021	1187.50
SHELL 281	二阶单元	—	14.078	1250.03

  当高度为壳的厚度即壳的中面位于 xoz 平面内时，无论网格如何划分，减缩积分与完全积分计算结果居然一样？？？

3.3 实体单元建模

1. 前处理

/PREP7

    BLOCK, 0, L, -H/2, H/2, -B/2, B/2    ! 创建六面体

    /VIEW,, 1, 1, 1                      ! 设置视图
    VPLOT                                ! Plot Volumes

    ET, 1, SOLID185                      ! 185减缩积分单元
    KEYOPT,1,2,1

    ET, 2, SOLID185                      ! 185完全积分单元
    KEYOPT,2,2,0

    ET, 3, SOLID186                      ! 186减缩积分单元
    KEYOPT,3,2,0

    ET, 4, SOLID186                      ! 186完全积分单元
    KEYOPT,4,2,1


    MP, EX, 1, E
    MP, NUXY, 1, nu

    TYPE, 3                              ! 设置单元类型
    MAT, 1

    ALLSEL,ALL                           ! 全显示

    LSEL,S,LENGTH,,H                     ! 按高度选线
    LESIZE,ALL, , ,NumH, , , , ,0        ! 布种子
    ALLSEL,ALL

    LSEL,S,LENGTH,,B                     ! 按宽度选线
    LESIZE,ALL, , ,NumB, , , , ,0        ! 布种子
    ALLSEL,ALL


    LSEL,S,LENGTH,,L                     ! 按长度选线
    LESIZE,ALL, , ,NumL, , , , ,0        ! 布种子
    ALLSEL,ALL

    VSWEEP,ALL                           ! 扫略分网

FINISH

2. 加载与求解

FINISH

/SOLU

    ! Fixed end
    NSEL, S, LOC, X, 0
    D, ALL, ALL, 0

    ALLSEL,ALL                           ! 全显示

    ! Loads at free end
    N1 = NODE(L, 0, 0)
    F, N1, FY, -F

    SOLVE
    
FINISH

悬臂梁实体单元有限元模型

3. 后处理

/POST1

    /VIEW,, 0, 0, 1
    PLDISP, 2
    PLNSOL, S, X

    N2 = NODE(L/2, 0, 0)                 ! 获取跨中结点编号
    NSEL,S,NODE, ,N2                     ! 选择结点
    PRNSOL,U,Y                           ! 列出结点位移

FINISH

表 3 实体单元有限元分析

单元名称	单元阶次	积分方案	跨中竖向位移 (mm)	最大拉伸弯曲正应力 (Mpa)
SOLID 185	一阶单元	减缩积分	18.641	789.692
SOLID 185	一阶单元	完全积分	2.6014	185.164
SOLID 186	二阶单元	减缩积分	14.077	1250.21
SOLID 186	二阶单元	完全积分	14.051	1220.71

四、结果整理

  由公式 (1)，可计算得该悬臂梁受载后跨中位置处的挠度为：

$$y_{\mathrm{m}\mathrm{i}\mathrm{d}}=-\frac{500\times 180^3}{6\times 200000\times 864}+\frac{500\times 360\times 180^2}{2\times 200000\times 864}+\frac{500\times 12^2}{8\times 100000\times 864}$$

$$=-2.8125+16.875+0.0375=14.10\mathrm{m}\mathrm{m}$$

  由此可以看出，因剪切变形产生的挠度仅为总挠度的 0.266 % 。

  由公式 (2)，可计算得该悬臂梁受载后固定端截面上最大弯曲正应力为：

$${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{500\times 360}{144}=1250\mathrm{M}\mathrm{p}\mathrm{a}$$

  现将各有限元计算得到的位移结果与应力结果分别除以相应的弹性力学解答，结果见下表。

表 4 有限元近似解与弹性力学理论解的位移比和应力比

单元名称	单元阶次	积分方案/形函数次数	位移比	应力比
BEAM 188	一阶单元	Linear	0.99496	0.99496
BEAM 188	一阶单元	Quadratic	0.99837	1.0
BEAM 188	一阶单元	Cubic	0.99837	1.0
BEAM 189	二阶单元	—	0.99837	1.0
SHELL181	一阶单元	减缩积分	1.29220	0.61752
SHELL181	一阶单元	完全积分	0.98759	0.94366
SHELL 281	二阶单元	—	0.92362	0.92780
SOLID 185	一阶单元	减缩积分	1.32206	0.63175
SOLID 185	一阶单元	完全积分	0.18450	0.14813
SOLID 186	二阶单元	减缩积分	0.99837	1.00017
SOLID 186	二阶单元	完全积分	0.99652	0.97657

  上表中各比值越接近于 1 表明有限元近似解与弹性力学的理论解越接近，小于 1 表明有限元近似解小于弹性力学理论解，大于 1 表明有限元近似解大于弹性力学理论解。

4.1 位移下限性

  由表 4 的位移比结果可以看出，除了一阶减缩积分单元外，其余单元计算的位移值均小于理论值 (绿色数值)，这也验证了有限元位移解答的下限性。值得注意的是，有限元位移解答的下限性是指有限位移元得到的位移解 总体上 不大于精确解。重要的事情说三遍是 总体上、总体上、总体上 ! 因此，上述各种对比分析，试图只单纯的讨论跨中截面形心点这一点处的竖向位移，来证实有限元位移解答的下限性，显然，不甚很合理。

  由表 4 的位移比结果可以看出，一阶减缩积分单元计算的位移值大于理论值 (绿色数值)，这与有限元位移解的下限性并不冲突，这是由于减缩积分单元 (一阶减缩积分单元更常见) 所独有的沙漏效应引起的，沙漏效应使单元变得更柔软，计算的位移值过分的偏大。

4.2 单元类型对计算结果的影响

  由表 4 可以看出，在如此粗糙的网格划分情况下，采用梁单元建模无论是位移还是应力均得到了很好的解答。而壳单元和实体单元计算所得的结果离散性较大。

  对于梁的模拟，一阶减缩积分壳单元、一阶减缩积分实体单元和一阶完全积分实体单元模拟效果很差，因此分析时不建议采用。

  若采用壳单元模拟梁结构，则建议采用一阶完全积分单元，总体上，位移和应力精度尚可。

  若采用实体单元模拟梁结构，则建议采用二阶减缩积分单元，该单元在计算精度和计算成本上达到了完美的平衡。

  如无特别的必要，对于细长结构还是应尽量采用梁单元来模拟，不但计算成本低，而且计算精度高，尤其是 BEAM 189 单元表现的相当优秀。

4.3 单元阶次对计算结果的影响

  显然，总体上讲，二阶单元 (高次形函数单元) 的计算精度高于一阶单元 (低次形函数单元) 。

4.4 积分方案对计算结果的影响

  由表 4 可以看出，减缩积分单元计算出的位移值总体上偏大，但不一定会小于理论上的精确解答。这是由于，减缩积分方案只保证形函数的完全多项式部分的精确积分，非完全的高次项在积分时精度得不到保证，这相当于对原位移函数做了调整，降低了离散模型的刚度，使位移计算较完全积分偏大，在一定程度上改善了计算精度。

  附带指出，若某一单元它的插值函数中只包括完全多项式的项次，如采用面积坐标边界为直线的三角形单元，那么，不存在完全积分与减缩积分的差别，只要根据被积函数的方次选择精度阶次和其一致的积分方案即可。

  对于一阶减缩积分单元，计算的位移值不是一般的偏大，而是过分的大，由表 4 可以看到，一阶减缩积分实体单元计算的位移值比理论值大了 29.22 %，显然这不主要是由上述原因引起的，而是一阶减缩积分单元所具有的 沙漏效应 引起的。

When using uniform reduced integration, verify the solution accuracy by comparing the total energy (SENE label on the ETABLE command) and the artificial energy (AENE label on the ETABLE) introduced by hourglass control. If the ratio of artificial energy to total energy is less than five percent, the solution is generally acceptable. If the ratio exceeds five percent, refine the mesh. You can also monitor the total energy and artificial energy during the solution phase of the analysis (OUTPR,VENG).

  为了控制能够降低求解精度的沙漏效应，求解器在计算过程中引入 artificial energy (伪能量)，计算完成后，若 artificial energy 与 total energy (应变能) 之比小于 5% ，则认为计算结果可接受。

  以上文提到的一阶减缩积分实体单元为例，分析完成后，运行如下命令流可查看各单元的伪能量与应变能。

/POST1  
    ETABLE, ,SENE                        ! 定义单元SENE表
    ETABLE, ,AENE                        ! 定义单元AENE表
    PRETAB,SENE,AENE                     ! 列出单元AENE表
FINISH

  各单元的伪能量、应变能与能量比如下表所示。

表 5 各单元 AENG、SENG及能量比

  由表 5 可以看到，第37、38、39 和 40 号单元的能量比均远大于 5%，因此，计算结果不可接受。

  事实上，只要保证总的伪能量与总的应变能之比小于 5% 即可，对于该算例总的伪能量即各单元伪能量之和为 38.9772，总的应变能为 14894.3，两者之比为 0.2617%，远小于 5%，满足收敛性要求。但事实上，位移值已为理论值的 132.21%，而最大正应力值仅为理论值的 61.75%，显然，该结果难以接受。

  那么，说了一堆，一阶减缩积分单元以后不要采用了，此单元存在还有存在的意义么？

  当然不是，该算例中虽然从能量比的角度已满足收敛性要求，但从实际上来看，显然结果误差很大，难道上面提到的收敛标准有问题，当然也不是。最根本的原因是本算例特殊，悬臂梁集中荷载作用下产生的变形是以弯曲变形为主导的，而一阶线性减缩积分单元的软肋就是对于此种问题的处理，需要格外的小心才能得到满意的结果，若研究的问题本身不是以弯曲变形为主导即弯曲变形与剪切变形势均力敌，至少剪切变形相较于弯曲变形不可忽略，那么当能量比小于 5% 时，无论是变形还是应力得到的结果肯定都相当精确。

  37、38、39 和 40 号单元的位置如下图所示，即为环绕荷载作用点的 4 个单元，由于这 4 个单元的能量比远大于 5，这表明沙漏效用在这四个单元中已经很显著，由此计算得到的荷载作用点处的位移肯定比理论值更大。 经分析，有限元算得的该点处的竖向位移为 59.733 mm，而由材料力学计算得到的该点的竖向位移为 45 mm，有限元近似解为理论解的 132.74%。

  避免一阶减缩积分单元出现沙漏效应有效的措施： 一是合理细化网格，二是改变荷载作用方式即可将单一结点荷载分散到多个结点上。对于一阶减缩积分实体单元，不同网格划分的有限元分析结果如下表所示。

表 6 一阶减缩积分实体单元不同网格划分的有限元分析

  由上表可以看出，宽度方向即荷载垂直方向上的网格数对计算结果没有任何影响；梁的长度方向即沿梁的轴线方向上的网格数对计算结果影响也不大，甚至是有不利影响。而最有效地控制沙漏效用的措施是加密高度方向即荷载作用方向上的网格数目，由 2 变为 4 后，位移的计算精度明显提高，当然，此方向上的网格数目越多位移的计算结果约精确。高度方向的网格数目为 10 时，位移的精度已经相当精确，但最大应力值为 1128.74 Mpa，为理论值的 90.3%。由此，可以得出如下结论，以弯曲变形为主导的问题不适于采用一阶减缩积分实体单元进行分析。如必须采用，则应注意网格的划分及能量比的查看，采用不同的网格进行对比分析，直至能量比在 5% 以内。但还是建议， 以弯曲变形为主导的问题，就别采用一阶减缩积分实体单元了，即使最终的能量比远远小于 5%，由此计算得到的应力值在一定程度上仍然有很大的误差。

Higher-order elements (actually second-order here) do not have shear locking due to the quadratic shape functions, but they are prone to volumetric locking if the full-integration scheme is used. To avoid this problem, higher-order elements use uniform reduced integration as the default or only option.

  由表 4 可以看出，一阶完全积分实体单元计算的位移值明显偏小，不是一般的小，是相当小，仅为理论值的 18.45%，这是由于一阶完全积分单元的 剪力自锁 效应引起的，使得单元显得过分刚硬。避免一阶完全积分单元出现剪力自锁效用的有效措施是合理细化网格。对于一阶完全积分实体单元，不同网格划分的有限元分析结果如下表所示。

表 7 一阶完全积分实体单元不同网格划分的有限元分析

  由上表可以看出，最有效地控制剪力自锁效用的措施是加密梁长度方向即梁轴线方向上的网格数目。加密高度和宽度方向上的网格数目效果不显著。

  二阶单元存在剪力自锁问题，但如果二阶单元发生扭曲或它的弯曲应力有梯度，也将有可能出现某种程度上的自锁。线性减缩积分单元能够很好的承受扭曲变形，因此，在扭曲变形为主导的问题分析中可采用合理网格细化的线性减缩积分单元。

Even in a linear analysis, lower-order elements are prone to shear locking in bending-dominated problems. Both lower- and higher-order elements experience volumetric locking when materials are incompressible or nearly incompressible, yet most nonlinear materials are either nearly or fully incompressible.

4.5 网格尺寸对计算结果的影响

  显然，有限元网格越细小，计算结果越精确，这是常识。当然，网格越精细，需要消耗的计算资源也就越多。应合理的选择网格尺寸，以平衡计算精度和计算时长。

  事实上，当网格尺寸足够细小时，单元阶次 (一阶单元与二阶单元)、积分方案 (减缩积分与完全积分)、解答类型 (结点解与单元解) 等方面造成的计算结果上的差异基本消除，解答趋近于一致且收敛于理论上的精确解。

  以上文中的悬臂梁为例，采用实体单元进行分析，网格尺寸设置为 0.5mm × 0.5mm × 0.5mm ，即 NumL=720、NumB=12、NumH=24。并选取 82285 号单元上的同一位置处的结点 (一阶单元为20567号结点、二阶单元为62732号结点) 为研究对象，查看其 x 方向正应力的单元解数值和结点解数值 ，各计算结果见下表。

实体单元有限元网格

82285号单元


表 8 0.5mm×0.5mm×0.5mm网格各实体单元有限元分析结果

  显然，由上表可以看到，同一结点处的应力结点解与单元解已经完全一致。与表 3 相比，网格细化并没有使位移精度进一步提高，事实上表 3 计算结果已相当精确，没有什么提高的空间了，反而使应力解答发生了细微振荡，因此，过分的细化网格是没有必要的。在如此细小的网格下，一阶减缩积分单元的最大正应力解答仍然具有一定的误差，可见以弯曲变形为主导的问题，无论如何也要避免采用此种单元进行分析。

五、截面内力提取

  本部分主要介绍后处理中某截面的内力提取，以本文算例中的悬臂梁跨中截面为例，根据材料力学知识可计算得悬臂梁跨中截面上的弯矩和剪力分别为：

$$M_{\mathrm{z}}=+500\times 180=90000\mathrm{N}\cdot \mathrm{m}\mathrm{m}$$

$$F_{\mathrm{s}\mathrm{y}}=+500\mathrm{N}$$

其中，z 轴为截面横轴，y轴为截面竖轴，如下图所示。

材料力学中梁横截面坐标系

5.1 梁单元端部内力的提取

5.1.1 ETABLE

  当悬臂梁采用梁单元建模时，命令流见本文 3.1 节，跨中截面为 5 号单元的终端截面或 6 号单元的始端截面，如下图所示，在进行内力提取时，选择哪一个单元作为提取对象均可，但要注意区分梁的始端和终端。

  若以 5 号单元为提取对象，其终端截面处按如下命令提取内力：

FINISH
/POST1

*DIM,MidSecForce5J,ARRAY,1,6                         ! 创建1行6列名称为MidSecForce5J的数组

elenum = 5                                           ! 单元编号为5

ESEL,S,ELEM, ,elenum                                 ! 选择5号单元

RSYS,0                                               ! 以下内力值采用全局直角坐标系 (对于梁单元无效)

ETABLE,JEndN,SMISC,14                                ! 梁单元终端截面沿局部x轴方向轴力
ETABLE,JEndFsy,SMISC,19                              ! 梁单元终端截面沿局部y轴方向剪力
ETABLE,JEndFsz,SMISC,18                              ! 梁单元终端截面沿局部z轴方向剪力
ETABLE,JEndT,SMISC,17                                ! 梁单元始终端面绕局部x轴扭矩
ETABLE,JEndMy,SMISC,15                               ! 梁单元终端截面绕局部y轴弯矩
ETABLE,JEndMz,SMISC,16                               ! 梁单元始终端面绕局部z轴弯矩


*GET,MidSecForce5J(1,1),ELEM,elenum,ETAB,JEndN       ! 获取梁单元终端截面沿局部x轴方向轴力数值
*GET,MidSecForce5J(1,2),ELEM,elenum,ETAB,JEndFsy     ! 获取梁单元终端截面沿局部y轴方向剪力数值
*GET,MidSecForce5J(1,3),ELEM,elenum,ETAB,JEndFsz     ! 获取梁单元终端截面沿局部z轴方向剪力数值
*GET,MidSecForce5J(1,4),ELEM,elenum,ETAB,JEndT       ! 获取梁单元终端截面绕局部x轴扭矩数值
*GET,MidSecForce5J(1,5),ELEM,elenum,ETAB,JEndMy      ! 获取梁单元终端截面绕局部y轴弯矩数值
*GET,MidSecForce5J(1,6),ELEM,elenum,ETAB,JEndMz      ! 获取梁单元终端截面绕局部z轴弯矩数值

  内力提取结果如下图所示，由下图可以看到，5 号单元终端截面上，JEndN = 0、JEndFsy = 0 、JEndFsz = - 499.999、JEndT = 0 、JEndMy = 89999.999、JEndMz = -7.401，这与理论计算相差无几。此外，也可以看出，即使设置了结果坐标系为整体直角坐标系，ANSYS 仍然采用梁单元的局部坐标系输出各内力值。

  若以 6 号单元为提取对象，其始端截面处按如下命令提取内力：

FINISH
/POST1

*DIM,MidSecForce6I,ARRAY,1,6                         ! 创建1行6列名称为MidSecForce6I的数组

elenum = 6                                           ! 单元编号为6

ESEL,S,ELEM, ,elenum                                 ! 选择6号单元


RSYS,0                                               ! 以下内力值采用全局直角坐标系 (对于梁单元无效)


ETABLE,IEndN,SMISC,1                                 ! 梁单元始端截面沿局部x轴方向轴力
ETABLE,IEndFsy,SMISC,6                               ! 梁单元始端截面沿局部y轴方向剪力
ETABLE,IEndFsz,SMISC,5                               ! 梁单元始端截面沿局部z轴方向剪力
ETABLE,IEndMy,SMISC,2                                ! 梁单元始端截面绕局部y轴弯矩
ETABLE,IEndMz,SMISC,3                                ! 梁单元始端端面绕局部z轴弯矩
ETABLE,IEndT,SMISC,4                                 ! 梁单元始端端面绕局部x轴扭矩


*GET,MidSecForce6I(1,1),ELEM,elenum,ETAB,IEndN       ! 获取梁单元始端截面沿局部x轴方向轴力数值
*GET,MidSecForce6I(1,2),ELEM,elenum,ETAB,IEndFsy     ! 获取梁单元始端截面沿局部y轴方向剪力数值
*GET,MidSecForce6I(1,3),ELEM,elenum,ETAB,IEndFsz     ! 获取梁单元始端截面沿局部z轴方向剪力数值
*GET,MidSecForce6I(1,4),ELEM,elenum,ETAB,IEndT       ! 获取梁单元始端截面绕局部x轴扭矩数值
*GET,MidSecForce6I(1,5),ELEM,elenum,ETAB,IEndMy      ! 获取梁单元始端截面绕局部y轴弯矩数值
*GET,MidSecForce6I(1,6),ELEM,elenum,ETAB,IEndMz      ! 获取梁单元始端截面绕局部z轴弯矩数值

  6 号单元始端截面各内力值，如下图所示：

5.1.2 FSUM

  除了采用 ETABLE 命令外，还可采用 FSUM 命令获取梁端部截面上的内力值，详见本文 5.2 实体单元截面内力的提取。

5.2 实体单元截面内力的提取

  若采用实体单元建模，则某一截面上的内力提起相对复杂，其基本操作步骤如下：

    首先，选取隔离体 (单元)，将要提取内力的截面暴露出来。

    其次，指定 SUMMATION POINT ，该点为提取内力截面的形心，该点是计算弯矩的参考点。

    然后，选取提取内力截面上的全部结点。

    最后，指定结果坐标系，并执行 FSUM 命令。

  SUMMATION POINT 是计算弯矩的参考点，若不指定，系统默认采用整体笛卡尔坐标系的原点，即 (0, 0, 0)，在计算弯矩时，各结点处的集中力对该点取矩，合成后作为该截面的弯矩值，显然，应正确的设置该点的位置，否则得到的截面弯矩值是错误的。

  在有限元计算时，首先进行的是单元分析，即确定单元结点力与单元结点位移的关系，然后进行的是结点分析，即有限元模型中全部结点满足平衡条件，显然，第三步中选择的结点为单元的附属结点，而非有限元模型的结点。

  综上，执行悬臂梁跨中截面上内力提取的命令流如下：

FINISH
/POST1
ALLSEL,ALL

CSYS,0
LOCAL,107,CART,L/2,0,0             ! 创建局部直角坐标系107

SPOINT, ,L/2,0,0                   ! 指定SPoint为(L/2,0,0),该坐标值为全局笛卡尔坐标系下的坐标值

CSYS,0
NSEL, S, LOC, X, L/2-1,L+1         ! 选择跨中截面及其右侧的全部结点(选左侧一样)
ESLN,S,1                           ! 全部结点选中单元才被选中(选择隔离体单元)

NSEL, S, LOC, X, L/2               ! 选择跨中截面上的全部结点

RSYS,107                           ! 指定107号坐标系为当前结果坐标系
FSUM,RSYS                          ! 在当前激活的结果坐标系下求截面合力

  FSUM： Sums the nodal force and moment contributions of elements.

  由上图可以看到，由此得到的内力值 Fy = -500、Mz = -90000，这在数值上与上述按材料力学算得的内力值一致，但符号有出入，下文讨论。当选取跨中截面左侧为隔离体时，得到的跨中截面内力如下图所示。

FINISH
/POST1
ALLSEL,ALL

CSYS,0
LOCAL,107,CART,L/2,0,0             ! 创建局部直角坐标系107

SPOINT, ,L/2,0,0                   ! 指定SPoint为(L/2,0,0),该坐标值为全局笛卡尔坐标系下的坐标值

CSYS,0
NSEL, S, LOC, X, -1,L/2+1          ! 选择跨中截面及其左侧的全部结点
ESLN,S,1                           ! 全部结点选中单元才被选中(选择隔离体单元)

NSEL, S, LOC, X, L/2               ! 选择跨中截面上的全部结点

RSYS,107                           ! 指定107号坐标系为当前结果坐标系
FSUM,RSYS                          ! 在当前激活的结果坐标系下求截面合力

  由上图可以看到，由此得到的内力值 Fy = 500、Mz = 90000，这与取右侧梁作为隔离体时得到的内力值在符号上完全相反。

  综上，可以看到，在 ANSYS 中，由 FSUM 计算得到的某个截面上的内力，在符号上与选取的隔离体有关。这是因为，FSUM 命令得到的各内力值采用结果坐标系标定，各内力数值的正负与所选定的结果坐标系有关，且有此得到的也并不是材料力学上所提及的截面内力，而是截面内力的反作用力。FSUM 命令在帮助文档中的解释为：

Sums the nodal force and moment contributions of elements.

  这说明，FSUM 命令得到的内力为所选择隔离体单元作用在所选择结点上的力，而材料力学所提及的内力，是作用在这个截面上单元的力，因此，FSUM 命令得到的内力取相反数后即为所选隔离体上该截面内的内力，各内力方向按所选结果坐标系标定。于是，当采用实体单元建模时，某一截面上的内力值按如下命令确定：

  这说明，FSUM 命令得到的内力为所选择隔离体单元作用在所选择结点上的力，而材料力学所提及的内力，是作用在这个截面上单元的力，因此，FSUM 命令得到的内力取相反数后即为所选隔离体上该截面内的内力，各内力方向按所选结果坐标系标定。于是，当采用实体单元建模时，某一截面上的内力值按如下命令确定：

*DIM,MidSecForceFromPartRigth,ARRAY,1,6    ! 创建1行6列名称为MidSecForceFromPartRigth的数组

FINISH
/POST1
ALLSEL,ALL

CSYS,0
LOCAL,107,CART,L/2,0,0             ! 创建局部直角坐标系107

SPOINT, ,L/2,0,0                   ! 指定SPoint为(L/2,0,0),该坐标值为全局笛卡尔坐标系下的坐标值

CSYS,0
NSEL, S, LOC, X, L/2-1,L+1         ! 选择跨中截面及其右侧的全部结点(选左侧一样)
ESLN,S,1                           ! 全部结点选中单元才被选中(选择隔离体单元)

NSEL, S, LOC, X, L/2               ! 选择跨中截面上的全部结点

RSYS,107                           ! 指定107号坐标系为当前结果坐标系
FSUM,RSYS                          ! 在当前激活的结果坐标系下求截面合力


*GET,FxVal,FSUM,0,ITEM,FX          ! 提取FX值将其赋予给变量FxVal
*GET,FyVal,FSUM,0,ITEM,FY          ! 提取FY值将其赋予给变量FyVal
*GET,FzVal,FSUM,0,ITEM,FZ          ! 提取FZ值将其赋予给变量FzVal
*GET,MxVal,FSUM,0,ITEM,MX          ! 提取MX值将其赋予给变量MxVal
*GET,MyVal,FSUM,0,ITEM,MY          ! 提取MY值将其赋予给变量MyVal
*GET,MzVal,FSUM,0,ITEM,MZ          ! 提取MZ值将其赋予给变量MzVal


MidSecForceFromPartRigth(1,1) = -FxVal
MidSecForceFromPartRigth(1,2) = -FyVal
MidSecForceFromPartRigth(1,3) = -FzVal
MidSecForceFromPartRigth(1,4) = -MxVal
MidSecForceFromPartRigth(1,5) = -MyVal
MidSecForceFromPartRigth(1,6) = -MzVal

  值得注意的是，此种内力提取方式同样适用于梁单元，但与单元表方式提取内力不同的是，采用此种方式可以指定内力输出的结果坐标系，即内力值方向可采用用户自定义的结果坐标系标定，而采用单元表方式输出的梁单元杆端内力，只能采用梁的局部坐标系。

内力方向采用结果坐标系标定

5.3 壳单元内力的提取

六、支座反力提取

6.1 FSUM

  采用 FSUM 方式也可提前支座反力，此时，无需选择隔离体，直接选取支座处截面上的全部结点后，执行 FSUM 命令，即可得到该支座截面上受到的支座反力，该方式得到的内力值在符号上与实际情况相反，需要取相反数，原因仍然是 FSUM 得到的内力是单元作用在结点上的力，而我们需要的材料力学中所指的内力是结点反作用于单元上的力，各内力方向采用所激活的结果坐标系标定，命令流如下：

FINISH
/POST1
ALLSEL,ALL

NSEL, S, LOC, X, 0            ! 选择跨中截面上的全部结点
SPOINT, ,0,0,0                ! 指定SPoint为(0,0,0),该坐标值为全局笛卡尔坐标系下的坐标值
RSYS,0                        ! 指定全局笛卡尔坐标系为当前结果坐标系
FSUM,RSYS                     ! 在当前激活的结果坐标系下求截面合力
*DIM,Reaction,ARRAY,1,6       ! 创建1行6列名称为Reaction的数组，用于存储各支座反力数值。

*GET,FxVal,FSUM,0,ITEM,FX     ! 提取FX值将其赋予给变量FxVal
*GET,FyVal,FSUM,0,ITEM,FY     ! 提取FY值将其赋予给变量FyVal
*GET,FzVal,FSUM,0,ITEM,FZ     ! 提取FZ值将其赋予给变量FzVal
*GET,MxVal,FSUM,0,ITEM,MX     ! 提取MX值将其赋予给变量MxVal
*GET,MyVal,FSUM,0,ITEM,MY     ! 提取MY值将其赋予给变量MyVal
*GET,MzVal,FSUM,0,ITEM,MZ     ! 提取MZ值将其赋予给变量MzVal

Reaction(1,1) = -FxVal        ! 令数组Reaction第1行第1列元素值为FxVal的相反数
Reaction(1,2) = -FyVal        ! 令数组Reaction第1行第2列元素值为FyVal的相反数
Reaction(1,3) = -FzVal        ! 令数组Reaction第1行第3列元素值为FzVal的相反数
Reaction(1,4) = -MxVal        ! 令数组Reaction第1行第4列元素值为MxVal的相反数
Reaction(1,5) = -MyVal        ! 令数组Reaction第1行第5列元素值为MyVal的相反数
Reaction(1,6) = -MzVal        ! 令数组Reaction第1行第6列元素值为MzVal的相反数

七、结论

  综上，通过各种分析我们可以看到，在进行有限元分析时，应根据所分析结构的几何特性选择合适的单元进行离散，这样才能事半功倍，以最小的代价获得最高的计算精度。正如本文所进行的悬臂梁计算，若采用梁单元分析，即使很粗糙的网格，也能得到相当高的精度；若采用实体单元分析，不但精度一般，计算时间长，还需要特别关注单元本身所具有的缺陷如沙漏效应、剪力自锁效应等，这些缺陷对计算结果影响很大。

  对于细长结构尽量采用梁单元来模拟，优先选用 BEAM 189 单元，对于扁平结构尽量采用壳单元来模拟，优先选用 SHELL 281 单元，对于三个方向上尺寸相当的结构，则采用实体单元模拟。

  若所分析的问题以弯曲变形为主导，即剪切变形相比与弯曲变形可忽略不计时，采用一阶减缩积分单元分析易产生沙漏效用，使单元显得格外柔软，采用一阶完全积分单元分析易产生剪力自锁效应，使单元显得格外刚硬，此时应采用二阶单元，能有效地避免上述缺陷。若能采用二阶减缩积分单元，则在消耗更少机时的前提下，获得较二阶完全积分单元更高的计算精度，这是由于，减缩积分方案只保证形函数的完全多项式部分的精确积分，非完全的高次项在积分时精度得不到保证，这相当于对原位移函数做了调整，降低了离散模型的刚度，使位移计算较完全积分偏大，在一定程度上改善了计算精度。因此，ANSYS 程序默认二阶单元的积分方案为减缩积分。

  若所分析的问题不是以弯曲变形为主导，那么此时有限采用一阶减缩积分单元进行分析，但是为了保证计算精度，需将能量比控制在 5% 以内。

  总之，弯曲变形为主导的问题，采用二阶减缩积分单元分析，非弯曲变形主导的问题，采用一阶减缩积分单元分析且控制能量比在 5% 以内，一阶完全积分和二阶完全积分单元能不采用就不采用，主要是精度没有明显提高，还更耗时。

  网格细分是解决有限元收敛性问题的必杀技能，事实上，网格划分细种程度，有限元分析中的许多问题基本上已迎刃而解，如单元解与结点解差差别较大问题，完全积分与减缩积分对计算结果的影响问题等等。现实计算中，我们必须权衡计算精度与计算时常，过细的网格往往很消耗求解时间，工程师要做的最重要工作就是用差不多的网格得出精度相当高的计算结果。

  最后，值得注意的是以上结论只适用于静力分析，对于其他分析，这些结论是否仍然成立有待确认。

八、尾声

  以上，便是 ANSYS 悬臂梁 静力分析 的简单介绍。

  因篇幅有限，某些内容未做详细介绍，如有疑问，欢迎邮件交流。

  Email： liyang@alu.hit.edu.cn 。

  仅以此文为我 ANSYS 及有限元相关理论的学习做一个总结。

  与此同时，也希望能够为初学者/有需要的人提供多一点参考。

  本文仅用于个人学习，除此之外，无其他任何用途。

  因个人水平有限，文中难免有所疏漏，还请各位大神不吝批评指正。

  胸藏文墨怀若谷，腹有诗书气自华，希望各位都能在知识的 pāo 子里快乐徜徉。

  本文逻辑清楚，内容详实，引例丰富。

  欢迎大家点赞、评论及转载，转载请注明出处！

  为我打call，不如为我打款！

  最后，祝各位攻城狮们，珍爱生命，保护发际线！

---

九、参考文献

[01]. ANSYS工程分析-基础与观念. 李辉煌.

[02]. 均布荷载作用下悬臂梁的弹性力学解. 王小琴. 刘云帅. 冯英杰.

[03]. 弹性理论及其工程应用. 李志鹏.

[04]. 剪力自锁和沙漏. 张强.

[05]. ANSYS 单元公式.

[06]. 剪切锁死、体积锁死、沙漏模式、零能模式.

[07]. ANSYS 中计算实体截面内力的方法.

[08]. Ansys中节点力提取.

[09]. ANSYS Classical 中如何获取实体单元某截面的内力.

[10]. 减缩积分单元与全积分单元.