1、 有限元方法理论及其应用 上机实验报告 2016年 11 月 20日 1 目录 题目 1. 2 1 实验题目 2 2 实验目的 2 3 建模概述 3 4 计算结果分析讨论与结论 3 5 实验体会与小结 16 题目 2. 17 1 实验题目 17 2 实验目的 17 3 建模概述 17 4 计算结果分析讨论与结论 18 5 实验体会与小结 20 题目 3. 21 1 实验题目 21 2 实验目的 21 3 建模概述 22 4 计算结果分析讨论与结论 23 5 实验体会与小结 29 2 题目 1 1 实验题目 图示一个简支梁平面应力问题模型。梁截面为矩形，高度 h=160mm，长度L=1000mm
2、，厚度 t=10mm。上边承受均布压力 q =1N/mm2，材料 E=206GPa， =0.29。 X方向正应力弹性力学理论解为： )534()4(6 22223 hyhyqyxLh qx 分别应用 3节点三角形单元、 4节点线性等参元（完全积分、减缩积分、非协调模式）、 8节点二次等参元完全积分进行下列各项数值实验： 1）用粗网格求解梁中部应力分量 x 最大值和上下边法向应力分量，并通过精确解对采用不同单元的 x 计算精度进行对比分析； 2）对粗网格下梁中部铅直（ y向）位移进行对比分析； 3）通过多次网格加密，对比试验 3节点三角形单元和 8节点二次等参元的收敛速度。总结出研究结论，撰写实
3、验报告。 2 实验目的 通过实验了解单元网格形状的选取以及网格的粗细对计算精度的影响，通过比较各种网格的计算结果与精确解考察有限元解的收敛性，并了解简支梁受均布载荷时的应力状况。 3 3 建模概述 1 启动 Abaqus/CAE软件。 2 选择 part模块， 点击 Creat Part创建名为 beam的 2D Planar部件， 并按题目要求 设置 平板尺寸 。 3 在 Property模块中创建材料和截面属性。 在 Material中创建材料为 Steel弹性模量为 206000Mpa，泊松比 0.29。 用 Creat Section创建界面属性。 用 Assign Section给部
4、件赋予截面属性。 4 在 Assembly模块中定义装配件为 Independent。 5 设置分析步，定义边界条件以及施加载荷。 6 划分网格，粗网格为 6x10，分别建立 3节点三角形单元（ Tri， linear）， 4节点线性等参元（完全积分 Quad， Linear；减缩积分 Quad， linear， Reduced integration；非协调模式 Quad， Linear， Incompatible modes）和 8节点二次等参元（ Quad， Quadratic）。 7 创建并提交分析。 8 查看结果并分析。 4 计算结果分析讨论与结论 4.1 粗网格 下梁中部应力分量和
5、上下边法向应力对比 1 理论解： X方向正应力由下式计算： 已知 q=1N/mm2 ， h=160mm， L=1000mm，max 2hy 代入上式得： 4 6 6m a x 36 1 0 1 1 1 30 0 . 0 8 1 0 4 2 9 . 4 9 70 . 1 6 4 2 4 5x M P a 2 3节点三角形单元计算结果： x的应力云图 梁中部应力分量x变化曲线 上边法向应力分量 5 下边法向应力分量 梁中部应力分量x最大值为 17.03Mpa。 梁上边法向应力分量最大值为 -1.3428Mpa 梁下边法向应力分量最大值为 0.3428Mpa 3 4节点线性等参单元完全积分： x的应
6、力云图 梁中部应力分量x变化曲线 6 上边法向应力分量 下边法向应力分量 梁中部应力分量x最大值为 25.421Mpa。 梁上边法向应力分量最大值为 -2.40466Mpa 梁下边法向应力分量最大值为 1.40466Mpa 4 4节点线性等参单元减缩积分： x的应力云图 7 梁中部应力分量x变化曲线 上边法向应力分量 下边法向应力分量 梁中部应力分量x最大值为 24.7759Mpa。 梁上边法向应力分量最大值为 -1.22287Mpa 梁下边法向应力分量最大值为 0.222869Mpa 8 5 4节点线性等参单元非协调模式： x的应力云图 梁中部应力分量x变化曲线 上边法向应力分量 9 下边法
7、向应力分量 梁中部应力分量x最大值为 28.8806Mpa。 梁上边法向应力分量最大值为 -2.09545Mpa 梁下边法向应力分量最大值为 1.09545Mpa 6 8节点二次等参单元： x的应力云图 梁中部应力分量x变化曲线 10 上边法向应力分量 下边法向应力分量 梁中部应力分量 x 最大值为 29.670Mpa。 梁上边法向应力分量最大值为 -1.05957Mpa 梁下边法向应力分量最大值为 -0.0727158Mpa 将上述计算结果制作成表格如下： 表 3-1-1 梁中部的最大正向应力 3节点三角形单元完全积分 4节点线性等参元 8节点二次等参元完全积分 理论值 完全积分 减缩积分
8、非协调模式 梁中部应力分量 x 的最大值（ Mpa） 17.03 25.42 24.78 28.8 29.67 29.50 11 从表中可以看出， 3节点三角形单元完全积分精确性最差，产生的结果没有意义。 4节点线性等参单元完全积分的线性单元出现了剪力自锁，剪力自锁引起单元在弯曲时过于刚硬，造成应变偏小，实际结果偏差仍较大。 4节点线性等参单元完减缩积分单元存在沙漏问题（零能模式），使得结构过于柔软，变形偏大，在粗网格情况下，产生无意义的结果。与精确解相比，仍有一定的误差。 8节点二次等参元完全积分的精度最高。其次是四节点非协调元单元。 下面讨论上下边法向应力分量 : 3节点三角形单元完全积分
9、 4节点线性等参元 8节点二次等参元完全积分 完全积分 减缩积分 非协调模式 理论解 梁上边法向应力分量 的最大值（ Mpa） -1.34 -2.40 -1.22 -2.10 -1.06 -1 梁下边法向应力分量 的最大值（ Mpa） 0.34 1.40 0.22 1.10 -0.07 0 比较以上有限元结果，在较粗单元网格划分情况下，以上几种求解方法得到的结果基本一致，梁的上边界法向应力分量分布呈凹形，梁的下边界法向应力分量分布呈凸形，但是与弹性力学求得的理论解还是有差距，理论解表明，梁的上边界法向应力分量等于 -1，梁的下边界法向应力分量为 0。所以，在较粗单元网格划分情况下的误差还是很大
10、的。其中 8 节点二次单元精确度最高。 12 4.2 粗网格下梁中部最大位移 粗网格下梁中部最大位移： 3节点梁中部最大位移 为 -0.130665 4节点全积分梁中部最大位移 为 -0.168976 4节点减缩积分 梁中部最大位移 -0.226525 13 4节点非协调元梁中部最大位移 -0.193409 8节点二次完全积分梁中部最大位移 -0.197471 将上述计算结果制作成表格如下： 表 3-1-2 梁中部最大位移 3节点三角形单元完全积分 4节点线性等参元 8节点二次等参元完全积分 完全积分 减缩积分 非协调模式 梁中部最大位移（ mm） -0.130665 -0.168976 -0
11、.226525 -0.193409 -0.197471 由以上数据可知 8节点二次等参元完全积分和 4节点线性等参单元非协调元的计算结果相当，而 4节点线性等参单元减缩积分存在严重的沙漏，变形严重失真，完全不能反映实际情况。其次还可发现 4节点线性等参单元完全积分计算精度比三角形单元完全积分高。 14 4.3 通过网格加密对比 不同单元 的收敛速度 1 3节点三角形单元求解时间（ 106 ）为 11s 3节点三角形单元求解时间（ 30 50 ）为 11s 3节点三角形单元求解时间（ 60100 ）为 11s 15 2 8节点二次等参单元 求解时间（ 610 ）为 11s 8节点二次等参单元求解
12、时间（ 30 50 ）为 11s 8节点二次等参单元求解时间（ 60100 ）为 13s 由此可知：相同网格密度情况下，不同单元的收敛速度基本一致，网格加密后收敛速度 略有 变慢。 3节点三角形单元在加密后没有变化，可能原因是因为模型本来就足够简单，对于复杂模型而言 网格越密收敛速度肯定越慢。 16 5 实验体会与小结 通过此次实验，对三角形单元、四边形完全积分单元、四边形减缩积分单元、四边形非协调元和 8节点二次等参元完全积分有了深刻认识。 (1)3 节点三角形网格下划分的单元要比 4 节点单元精度要差； (2)粗网格要比细网格精度要差； (3)4 节点非协调单元和 8 节点二次积分等参单元
13、下，粗细网格的精度基本上接近； (4)4 节点线性等参单元减缩积分出现零能模式，求解的位移值失真，总刚奇异。 因此， 不同的积分单元、积分方法、积分阶次、网格密度会导致计算结果不同。如果单元算法选择不当，除影响到计算精度外，还有可能导致计算的失败，因此在使用时要注意适当的选择合适的计算单元。 17 题目 2 1 实验题目 设计一个采用减缩积分线性四边形等参元的有限元模型，通过数值试验来研究网格密度、位移约束条件与总刚度矩阵奇异性、沙漏扩展、求解精度的关系，并验证采用减缩积分时保证总刚度矩阵非奇异的必要条件。总结出你的研究结论，撰写实验报告。 2 实验目的 研究网格密度、位移约束条件与总刚度矩阵
14、奇异性、沙漏扩展、求 解精度的关系，并验证采用减缩积分时保证总刚度矩阵非奇异的必要条件。 3 建模概述 悬臂梁左端固定，右端受一大小为 M 6 N m的弯矩。梁的长度 l=0.3m,横截面高度 h=0.03m，宽度 b=0.01m。材料为钢，弹性模量 E=210GPa，泊松比u=0.3， 求梁的端部位移。 1 启动 Abaqus，在 part模块建立模型； 2 进入 property，设置材料属性：弹性模量 E=210GPa，泊松比 u=0.3； 3 进入 assembly，组装模型； 4 进入 step，设置静力分析步； 5 进入 interaction，在梁右上角添加参考点，并使之与梁固结
15、； 6 进入 load模块，在参考点施加大小为 6 N m的弯矩，并在做端添加固定约束。 在所有计算完成之后，删除约束的某一自由度； 7 划分网格，依次选择 1 10、 2 20、 4 20、 8 40四种网格密度，设置单元类型为 4节点减缩积分单元 ； 8 提交 job，分析并查看结果对比分析。 18 4 计算结果分析讨论与结论 1) 由材料力学求 解： 22 41 1 86 0 . 3 0 . 5 7 1 1 0 m2 2 2 . 1 1 0 2 . 2 5 1 0zMlv EI 1 网格为 1 10时： 梁端部位移 u=0.159 m 2 网格为 2 20时： 梁端部位移 -4u=0.7
16、56 10 m 19 3 网格为 4 20时： 梁端部位移 -4u=0.614 10 m 4 网格为 8 40时： 梁端部位移 -4u=0.588 10 m 5 分别删除梁左端 X方向的自由度、 Y方向的自由度，以及整个左端的约束： 提交 job之后分析中断， Abaqus提示分析失败，产生了刚性位移，结果不收敛无法继续计算。 20 2) 结果分析： 1 对于同一单元而言，网格越密计算结果的精确度越高，越接近于理论解。 当网格很稀疏时，缩减积分单元计算的端部位移值与理论解差距巨大，说明在这种情况下存在沙漏扩展现象，使得结构异常柔软，变形比真实大。 当网格密度增大之后，这种情况得到很大程度的缓解
17、，但相较于理论解仍然存在偏大的现象。同时也能说明，采用减缩积分单元时，刚度矩阵 K非奇异的必要条件不一定能满足， 尤其是在网格密度比较低的时候。会导致计算结果偏差很大。 2 当边界条件不足以使模型有足够的刚度时，比如缺少某个方向的自由度限制，就会产生刚性位移，使得计算结果不收敛，总刚度矩阵奇异，求解精度很差甚至无法求解计算。 3 总刚度矩阵非奇异的必要条件为 M ng d N ， 只有这个条件满足时才能使用缩减积分单元。网格密度为 1 10时，上述条件不满足，刚度矩阵奇异，计算结果误差很大。采用更高阶次单元和更精细的网格能使结果更精确。 5 实验体会与小结 对于不同的模型要综合考量各种不同单元
18、的效益。并且要在此之前进行 充分的考量，积累经验。比如缩减积分单元的好处是能简化计算，也能适应较大的变形，但是在使用之前就必须考虑总刚度矩阵的奇异性，避免出现沙漏现象，使得计算结果非常离谱。 21 题目 3 1 实验题目 一个矩形平板，长 1000mm,宽 100mm，厚度 5mm。材料的 E=200GPa，0.3 ， 3-6 /107 .8 2 mmKg 。板的一对短边简支。进行下列计算分析并撰写实验报告。 1) 在相同的 粗网格（厚度方向 1层单元）下， 分别用 8节点线性六面体全积分等参元、 8节点六面体非协调元、 20节点二次六面体全积分等参元计算其前六阶自由振动频率（ Hz），对计算
19、结果用表格作汇总对比和分析（不要附图）。 2) 采用适当网格密度的 8节点六面体非协调元，用隐式直接积分法对该平板进行瞬态响应计算实验。板的上表面受对称三角形时间脉冲均布压力，最大值0.05MPa 。要求：（ 1）通过计算获得平板下表面中心位置沿长边方向的正应力分量响应最大值随载荷脉宽（ 0.2毫秒、 1.0毫秒、 10毫秒、 100毫秒、静载）变化的曲线；（ 2）板的厚度改为 20mm，研究上述规律的变化，并进行归纳和对比分析讨论。 2 实验目的 通过实验掌 握线性六面体完全积分等参单元、二次六面体完全积分等参单元和六面体非协调元的应用，了解三种单元在有限元动力学求解中的不同及优缺点。 通过
20、实验掌握使用 8节点六面体非协调单元在隐式直接积分法进行动态响应的分析，了解有限元动力学问题求解和静力学问题的差别。 22 3 建模概述 3.1 1 启动 ABAQUS/CAE软件 2 按题目要求建立 Part模型 3 设置密度为 7.82 10-9 t/mm3，弹性模量为 200000Mpa，泊松比为 0.3。 4 设置截面属性，并将其赋予模型 5 进入 Assembly模块装配实体。然后进入 Step设置分析步 ，选择频率。并设置阶数为 6阶。 6 设置边界条件，两条短边简支。 7 划分网格，粗网格选择 4*40，分别用 8节点线性六面体全积分等参元、 8节点六面体非协调元、 20节点二次
21、六面体全积分等参元划分网格。 8 提交 Job，查看并分析结果。 3.2 1 模型尺寸材料设置如 3.1； 2 进入 Step模块，插入 Dynamic， Implicit分析步。并分别设置 time period为0.0002、 0.001、 0.01、 0.1； 3 设置三角形载荷，并按照载荷脉宽依次改变载荷形状。 4 划分网格，选择六面体 8节点非协调元； 5 提交 job，得出结果并列表对比分析。 6 复制 model，并改变模型厚度重复上述步骤。 23 4 计算结果分析讨论与结论 4.1 前六阶自由振动频率。 1 单元选用 8节点线性六面体全积分等参元。 2 单元选用 8节点六面体非
22、协调元。 24 3 单元选用 20节点二次六面体全积分等参元。 4 将所有数据导入 Excel对比分析。 可以看到，采用 8节点线性六面体单元时计算出来的不同阶数的自由振动频率误差很大；而采用 8节点六面体非协调元时误差很小，并且与 20节点二次六面体单元的值很接近，误差很小，与理论值更接近。 25 4.2 瞬态响应计算 4.2.1 平板厚度为 5mm时 1 载荷脉宽为 0.2ms时，下表面中心沿长边方向应力最大值为 -0.1999 Mpa； 脉宽为 0.2ms时平板 X轴方向的应力云图 2 载荷脉宽为 1ms时，下表面中心沿长边方向应力最大值为 -0.243 Mpa； 脉宽为 1ms时平板
23、X轴方向的应力云图 3 载荷脉宽为 10ms时，下表面中心沿长边方向应力最大值为 204.317 Mpa； 26 脉宽为 10ms时平板 X轴方向的应力云图 4 载荷脉宽为 100ms时，下表面中心沿长边方向应力最大值为 1256.28 Mpa； 脉宽为 100ms时平板 X轴方向的应力云图 5 静载荷为 0.05Mpa时 下表面中心沿长边方向应力最大值为 1498.91 Mpa； 静载荷时平板 X轴方向的应力云图 27 4.2.2 平板厚度为 20mm时 1 载荷脉宽为 0.2ms时，下表面中心沿长边方向应力最大值为 -0.08 Mpa； 脉宽为 0.2ms时平板 X轴方向的应力云图 2 载
24、荷脉宽为 1ms时，下表面中心沿长边方向应力最大值为 -1.85 Mpa； 脉宽为 1ms时平板 X轴方向的应力云图 3 载荷脉宽为 10ms时，下表面中心沿长边方向应力最大值为 115.566 Mpa； 脉宽为 10ms时平板 X轴方向的应力云图 28 4 载荷脉宽为 100ms时，下表面中心沿长边方向应力最大值为 95.0211 Mpa； 脉宽为 100ms时平板 X轴方向的应力云图 5 静载荷为 0.05Mpa时 下表面中心沿长边方向应力最大值为 93.73 Mpa； 静载荷时平板 X轴方向的应力云图 29 4.2.3 数据对比分析 由以上计算结果可以得出， 8节点六面体非协调单元和 2
25、0节点六面体完全积分等参元计算其前六阶自由振动频率和振型的精度较高，误差较小； 8节点六面体完全积分等参元计算其前六阶自由振动频率和振型的精度则较低，误差较大。 当载荷脉宽很小的时候，平板中间下表面最大应力值也很小。随着载荷脉宽的增大，最大应力值也在增大，并且逐渐趋近于静载荷时的最大应力值。并且可以看到 ，当板的厚度增加后，中间下表面最大应力 减小 了。 5 实验体会与小结 通过 本实验比较了不同类型单元对矩形平板自由振动频率和振型的影响。通过计算可以知道，高阶单元比低阶单元计算精度稍高一些。对物体的自由振动频率和振型有了深刻认识，掌握了用有限元软件进行模态分析和动态响应分析的方法，尤其是对材料的动态和静态响应特性有了深刻的体会，材料的动态力学特性和静态力学特性有着很大差异， 因而在考查冲击过载下材料的强度时，不能简单的用静态强度来代替。
