南京理工大学研究生 有限元方法理论及应用考试 个人答案

目录

1

1.1

等参单元及其应用 ........................................................................................................... 1

概述 ............................................................................................................................. 1 1.1.1 1.1.2 1.2

1.2.1 1.2.2 1.2.3 1.3 1.4 2

2.1 2.2 2.3 2.4 2.5 3

3.1

1.3.1

等参单元的概念、原理 ..................................................................................... 1 等参单元对有限元法工程应用的意义 ............................................................. 1 等参单元刚度矩阵的数值积分方法 ................................................................. 1 确定积分阶的原理 ............................................................................................. 2 全积分单元与减缩积分单元讨论和评价 ......................................................... 3 全积分、减缩积分线性等参单元有关问题的分析讨论 ................................. 3

等参单元的数值积分方法 ......................................................................................... 1

线性等参单元 ............................................................................................................. 3 等参单元的应用 ......................................................................................................... 5 四节点平面等参单元的收敛协调性 ......................................................................... 6 八节点平面等参单元 ................................................................................................. 8 3节点平面三角形单元 .............................................................................................. 9 20节点六面体等参单元 .......................................................................................... 10 20节点六面体等参单元 .......................................................................................... 11 实验一 ....................................................................................................................... 15 3.1.1 3.1.2 3.1.3 3.1.4 3.1.5

实验题目 ........................................................................................................... 15 实验目的 ........................................................................................................... 15 建模概述 ........................................................................................................... 15 计算结果分析与结论 ....................................................................................... 16 实验体会与总结 ............................................................................................... 32 实验题目 ........................................................................................................... 33 实验目的 ........................................................................................................... 33 建模概述 ........................................................................................................... 33 计算结果分析与讨论 ....................................................................................... 34 实验体会与总结 ............................................................................................... 36 实验题目 ........................................................................................................... 36 实验目的 ........................................................................................................... 36 建模概述 ........................................................................................................... 37 计算结果分析与结论 ....................................................................................... 37 实验体会与总结 ............................................................................................... 44

分析与计算 ......................................................................................................................... 6

上机实验 ........................................................................................................................... 15

3.2 实验二 ....................................................................................................................... 33 3.2.1 3.2.2 3.2.3 3.2.4 3.2.5

3.3 实验三 ....................................................................................................................... 36 3.3.1 3.3.2 3.3.3 3.3.4 3.3.5

1 等参单元及其应用

1.1 概述

1.1.1 等参单元的概念、原理

普通单元受到两个方面的：（1）单元的精度。单元的节点数越多，单元精度越高；（2）单元几何上的。普通矩形和六面体单元都不能模拟任意形状几何体，所有几种普通单元都是直线边界，处理曲边界几何体误差较大。

为了解决上述矛盾，方法就是突破矩形单元和六面体单元几何方面的，使其成为任意四边形和任意六面体单元，这类单元位移模式和形函数的构造和单元列式的导出不能沿用构造简单单元的方法，必须引入等参变换，采用相同的插值函数对单元的节点坐标和节点位移在单元上进行插值。这种单元叫做等参单元。所引入的等参变换意义为在把母单元映射为实际单元的几何变换过程中采用了与位移模式相同的参数（插值函数）。 1.1.2 等参单元对有限元法工程应用的意义

等参单元的原理对有限元法工程应用的意义有以下几点：（1）等参单元形状、方位任意，容易构造高阶单元，适应性好，精度高；（2）等参单元列式具有统一的形式，规律性强，采用数值积分计算，程序处理方便；（3）由于等参单元涉及单元几何形状的变换，对实际单元的形态有一定要求。单元形态好坏影响计算结果的精度。单元形态应满足：A单元各方向的尺寸尽量接近；B单元边界不能过于曲折，不能有拐点和折点，尽量接近直线或抛物线；C边之间夹角接近直角。（4）高阶等参元精度高，描述复杂边界和形状的能力强，所需单元少，在结构应力分析中应用最广泛；（5）借助于等参元可以对一般的任意几何形状的工程问题和物理问题方便地进行有限元离散。因此，等参元的提出为有限元法成为现代工程实际领域最有效的数值分析方法迈出了重要的一步；（6）由于等参变换的采用使实际问题物理坐标系内的等参单元的单元刚度、质量、阻尼、载荷等特性矩阵的计算仍能在局部的自然坐标表示的规则域内进行，实现了矩阵积分式内被积函数中的导数、体积微元、面积微元、线段微元的变换及积分限的置换。因此不管各个积分式的被积函数如何复杂，都可以方便地采用标准化的数值积分方法进行计算，从而使各类不同工程实际问题的有限元分析纳入了统一的通用化程序。 1.2 等参单元的数值积分方法

1.2.1 等参单元刚度矩阵的数值积分方法

计算单元特性矩阵的方法有Newton-Cotes积分、高斯积分、Irons积分、Hammer积分等。为了减少积分点的数目和便于程序编制，一般采用高斯数值积分方法。

一维高斯数值积分： I1f(x)dxf(x)Hii11ni

结论：n阶高斯积分公式对2n-1次多项式被积函数可求得精确积分。

1

11nn二位高斯数值积分：I11F(,)ddHHii1j1jF(i,i)

积分公式对ξ，η方向最高方次为2n-1的多项式可求得精确值。 1.2.2 确定积分阶的原理

1 保证刚度矩阵积分精度的积分阶选择 （1）完全积分方案

选择积分阶的基本考虑是保证被积函数所有项次精确积分，这种积分方案称为完全积分方案。

以二维单元刚度矩阵高斯数值积分为例： 单元刚度矩阵为Ke1111BTDBtJdd

根据这个原理，具有规则形状的单元（|J|=常数）完全积分方案如下： 二维4节点、三维8节点等参元分别是2×2、2×2×2积分； 二维8节点、三维20节点等参元分别是3×3、3×3×3积分。 对于|J| ≠常数，需要增加积分阶数。 （2）减缩积分方案

实际应用中选取的积分阶往往可以低于被积函数所有项次精确积分所需要的阶数，这种积分方案称为减缩积分。

对二、三维连续体单元，通常按形函数中完全多项式阶数 p 来选取积分阶，即取 n=p。

根据这个原理，具有规则形状的单元（|J|=常数）减缩积分方案如下： 二维4节点、三维8节点等参元分别是1×1、1×1×1积分； 二维8节点、三维20节点等参元分别是2×2、2×2×2积分。

实际计算表明，采用减缩积分往往可以取得较完全精确积分更好的精度。 2 保证结构总刚度矩阵非奇异的积分阶选择 （1）数值积分条件下刚度矩阵的秩

单元刚度矩阵数值计算公式为：KHiBiTDBiJi根据矩阵秩的基本规则有

ei1ngKengd，ng是高斯积分点数，秩KMngd，为保证正确求解，数值积分K的秩不能小于N，即KN，N是系统的自由度数，也就是系数矩阵K的阶数。因此系数矩阵K非奇异的必要条件是MngdN。 （2）减缩积分情况下刚度矩阵的奇异性

减缩积分得到的Ke往往缺秩，即Ke是奇异的，即存在单元发生非刚体位移模式情况下，单元应变能为零的情况——零能模式（沙漏模式）。因此要注意检查Ke奇异性的必要条件，即MngdN。

2

1.2.3 全积分单元与减缩积分单元讨论和评价

虽然减缩积分会出现零能模式，但在实际计算中减缩积分应用的十分广泛。实际计算表明减缩积分往往可以取得较完全精确积分更好的精度。这是主要由于：

(1) 完全积分是由被积函数中对应形函数高方次非完全部分的精确积分所要求。但这些非完全的最高方次项往往并不能提高精度，反而可能带来不好的影响。取较低阶的高斯积分，使积分精度正好保证完全多项式方次的要求，而不包括更高次的非完全多项式的要求，其实质是相当于用一种新的插值函数替代原来的插值函数，从而一定情况下改善了单元的精度。

(2) 基于最小位能原理的有限元位移法中，位移具有下限性质。有限元的计算模型比实际结构刚度偏大。减缩积分方案可以使有限元计算模型的刚度有所降低，因此往往有助于提高计算精度。

(3) 对于势能泛函中包含罚函数的情况，必须采用减缩积分以保证与罚函数对应的刚度矩阵是奇异的。 由于减缩积分单元的特性，因此在使用时须注意一些问题，在网格扭曲严重的情况下，优先使用网格细化的线性减缩积分单元，对于接触问题，采用细网格的线性减缩积分单元，对于规则单元等选用完全积分单元。 1.3 线性等参单元

1.3.1 全积分、减缩积分线性等参单元有关问题的分析讨论

全积分是由被积函数中对应形函数高方次非完全项部分的精确积分所要求。但位移模式中决定计算精度的是完全多项式的方次，高方次非完全多项式部分往往不能提高精度，并且可能带来负面影响！

如采用较低积分阶的减缩积分，可使积分精度正好保证完全多项式方次对应的被积函数精确积分的要求，而不包括更高次非完全多项式的要求，相当于对形函数作了某种修正，从而一定情况下改善了单元精度！这种情形的减缩积分称为“优化积分方案”。

基于最小势能原理的有限元位移法中，位移解具有下限性质。有限元计算模型比实际结构刚度偏大。减缩积分可以使模型刚度有所降低，因此往往有助于提高计算精度。

对于势能泛函中包含罚函数的情况，必须采用减缩积分以保证与罚函数对应的刚度矩阵是奇异的。 从物理角度看，采用减缩积分计算单元刚度矩阵时，存在单元发生非刚体位移模式情况下，单元应变为零的情况——零能模式（沙漏模式）。因此，采用减 缩积分时，需要检查总刚度矩阵非奇异性的必要条件：

MngdN3

应用减缩积分线性单元时，即使满足上面条件，也不能保证总刚度矩阵非奇异，实际计算表明，一个有限元模型中，特定的节点位移约束和载荷施加可能导致沙漏模式扩展，一般对这类单元应使用足够密的单元网格。

以下采用等截面悬臂梁端部受压，求端部位移的平面应力问题来具体说明全积分、减缩积分线性等参单元在计算精度、剪力自锁、零能模式等方面的差异。

悬臂梁左端固定，右端受一大小为F＝10 N的压力。梁的长度l=0.5m,横截面高度h=0.03m，宽度b=0.01m。材料为钢，弹性模量E=200GPa，泊松比u=0.3，材料力学理论解：

Fl3100.534v0.925910m118 3EIz32102.2510实验结果如下：

（410）：

全积分 减缩积分

(620)：

全积分 减缩积分

4

(16100)：

全积分 减缩积分

从计算结果发现：在同样网格划分密度情况下，全积分结果偏小，减缩积分结果偏大，随着网格密度增大，计算结果皆越来越接近理论值。 这是由于全积分的线性单元出现了剪力自锁，剪力自锁引起单元在弯曲时过于刚硬，造成应变偏小，同时随网格的加密，剪力自锁效应可以得到改善，但与实际结果偏差仍较大。

而减缩积分单元存在沙漏问题（零能模式），使得结构过于柔软，变形偏大，甚至在在粗划网格情况下，产生无意义的结果，但在加密网格后，结果能够改善，但依然偏大。

1.4 等参单元的应用

等参单元在有限元的发展中占有重要的位置，由于它能使自然坐标系内形状规则的单元变换成总体坐标系内的形状为扭曲的单元，从而为求解域是任意形状的实际问题的求解提供了有效的单元形式。两种坐标系内的变换通常采用和位移函数的插值形式，依据坐标变化插值点数和位置插值点数的比较，分别称之为等参元、超参元和次参元。通常应用最多的是两者插值点数相同的等参元。 等参单元通常用来模拟二维、三维实体单元。应用时可选择线性位移模式或二次位移模式，完全积分或减缩积分。二维四边形和三维六面体线性单元有非协调位移模式可选择。

(1) 平面问题中的单元使用

一般应力分析，优先使用8节点四边形二次单元（减缩积分、完全积分），其次使用四节点非协调元、6节点三角形二次单元，但非协调元要避免扭曲的单元形状。在网格扭曲严重的情况下，优先使用网格细化的线性减缩积分单元。 有弯曲变形情况下，避免使用4节点双线性位移模式完全积分单元。因为该单元有剪力自锁，会引起单元弯曲时过于刚硬。对于接触问题，采用细网格的线性减缩积分单元或者非协调模式单元。

5

(2) 三维问题中的单元选择

对于中小规模的一般结构分析，优先采用20节点六面体减缩积分、完全积分单元，应力集中区域，应采用完全积分单元。对大规模分析或接触问题，优先采用8节点六面体非协调模式单元。网格扭曲严重的情况下，应使用细网格8节点六面体减缩积分单元。要求快速建模情况下，可考虑采用10节点四面体二次单元，但得到的模型规模（节点数）大于等效网格的六面体单元，否则结果不精确。

2 分析与计算

2.1 四节点平面等参单元的收敛协调性

证明满足收敛的协调性准则即证明单元位移模式在单元交界面必须具有连续性。

4节点等参单元：

在母单元的位移模式为：

uN1u1N2u2N3u3N4u4 vN1v1N2v2N3v3N4v4

1其中，Ni(1i)(1,4i) (i1,2,3

4以母单元的13边（1边）为例证明： 4个节点的坐标分别为：

1(1,1N)1:1；(1)(1

41 2(1,1):N2(1)(1)；

41 3(1,1):N3(1)(1)；

414(1,1):N4(1)(1)；

4)

6

11在1边：N1(1)；N3(1)；N2N40；

221111 u(1)u(1u)u(u)u(u)131331

22221111 v(1)v(1v)v(v)vv)13133(1

2222可得：

uu10，(u3u1)=常数； 2vv10，(v3v1)常数； 2所以对于13边，位移沿单元边界线性变化，能保证单元的协调性。

其他边采用相同的原理： 对于12边：=-1

1111u(1)u1(1)u2(u1u2)(u2u1)

22221111v(1)v1(1)v2(v1v2)(v2v1)

2222可得：

uu10，(u2u1)=常数； 2vv10，(v2v1)常数； 2对于24边：=1

1111u(1)u2(1)u4(u2u4)(u4u2)

22221111v(1)v2(1)v4(v2v4)(v4v2)

2222可得：

uu10，(u4u2)=常数； 2vv10，(v4v2)常数； 2对于34边：=1

1111u(1)u3(1)u4(u3u4)(u4u3)

22221111v(1)v3(1)v4(v3v4)(v4v3)

22227

可得：

uu10，(u4u3)=常数； 2vv10，(v4v3)常数； 2综上，可证明4节点等参单元满足收敛的协调性准则。

2.2 八节点平面等参单元

建立如图2.2.1所示直角坐标系

图2.2.1直角坐标系下的平面等参单元

图2.2.2 局部坐标系下的等参单元

该单元在母单元中的位移模式为包含完全二次式的不完全三次多项式。插值基函数可以用形函数性质直接构造。由于只有2-6-3边受均布载荷作用，故只有2、6、3这三个节点有等效节点力，其对应的形函数为：

1N2(1)(1)(1)41 N6(12)(1)21N3(1)(1)(1)4又因为2-6-3边对应ξ=-1，进一步化简得

1y2y2N2(1)22bb4y22N6112

b1y2y2N3(1)22bb8

作用在边界上的面积力为T=q，单元等效节点力计算公式为：

PeeNTTtdS

s将N和T代入上式得

P2xN2qtdssb2b2y2y21qtdyqbt26bb4y221qtdyqbt2b3y2y21qtdyqbt26bbP2y0P6xN6qtdssb2b2

P6y0P3xN3qtdssb2b2P3y02.3 3节点平面三角形单元 单元应力：

Sae

其中，S称为应力矩阵：DBi,Bj,Bm；

biE0biDBi2(102)A10ci2e；

10bi2ciT0ciaui,vi,uj,vj,um,vm；

当发生刚体位移时，uiujumu； vivjvmv；

21S(aeae)SaeSaeDBu,v,u,v,u,v(bibjbm)u0(cicjcm)vE(cicjcm)v0(bibjbm)u22(10)A1100(cicjcm)u(bibjbm)v22

9

T

其中：

bi1yj1ymyjym；bj1ym1yiymyi；bm1yi1yjyiyj

bibjbm0；

ci1xj1xmxmxj；cj1xm1xixixm；cm1xi1xjxjxi

cicjcm0；

21S(aeae)SaeSae0

所以，三节点三角形单元发生刚体位移时，单元中不产生应力。 2.4 20节点六面体等参单元 三维高斯积分：

IF(,,)dddHiHjHmF(i,j,m)111m1j1i1111nnni,j,m1nHijmF(i,j,m)

如果F(,,)aijm,2n1，且i,jm则上式给出精确的积分结果。 ijm，

以单元刚度矩阵为例进行证明：KeBTDBdV

Ve

xx在等参单元变换中：dVd(dd)x111yzyzdddJddd

yz代入上式得到：Ke111TBDBJddd

可得：F(,,)BTDBJ

20节点六面体单元位移多项式包含的项分别为：

2222222222221,x,y,z,xy,xz,yz,x,y,z,xy,xy,xz,xz,yz,yz,xyz,xyz,xyz,xyz对应在等参元中：

2222222222221,,,,,,,,,,,,,,,,,,,

10

Nxbijmijm x取0到20，Nx在,,三个方向的阶次最高为2；

BxLNx：Nx在,,三个方向的阶次最高仍然为2；

刚度矩阵的被积函数F(,,)BTDBJ，则在,,三个方向的阶次最高仍然为224，即max(i)max(j)max(m)4；i,j,m2n1n3

所以，精确积分方案为333阶。

当J不是常数的时候，F(,,)BTDBJ在,,上的阶次不能确定，但是一定大于4，所以需要提高阶次。 2.5 20节点六面体等参单元

图2.5.1 三节点杆单元示意

把每个单元拿出来考虑，在单元内假设位移场是二次的，根据节点位移进行插值。

对单元1得到：

2s2s4s2uL2Lu11L22s2su2L2Lu3，矩阵形式为uNd 2s2s 2LL2s2s4s2插值函数矩阵为：N212LLL单元1的节点位移列阵为：du1u2u3

u4u5，N同上。

TT对单元2也有同样的位移场：uNd，du3显然，整个杆上，由各单元上假设的位移场拼接而成的位移试探函数是连续的，只要我们记住u10，得到的就是全域可能位移场。这样的位移场已经把

节点位移自由度作为广义坐标。

下面在上述分片位移试探函数的基础上进一步实施里兹法求解。

11

2LE2xAdxqudx

02在两个单元上分片进行总势能计算：p首先计算应变：x单元内有：x2L0u,xu,s

4s1 2LLd4s18sNdBd，B2dsL2LL所以一个单元内的应变能为：

UL2L2LL1212T1T2TExAdsLxAExdsdLBAEBdsd 222228173331TAE8168或者Udkd，k

2L333187333载荷：qcx。则载荷在两个单元内分别表达为：

Lqcs

23， qcLs

20LcLdTNTcdTsds2

6212带入外力功积分式，对各个单元分别计算外力功。第一个单元外力功为：

L2L2qudsuTqdsL2L2L2L21cldT6第二个单元外力功为： 622总势能是两个单元总应变能减去两个单元外力功。为了能够以矩阵形式相加，将单元势能矩阵表达式中的单元节点位移列阵d用结构整体位移向量

Du1u2u3u4u5代替，单元刚度矩阵k等与单元有关的矩（列）阵扩

T展成结构整体规模（5×5），则相加后系统总势能表达为：

12

8100700000300033008168007810031AE33AET333pD1872LL0081680033333300000187000000333000202cL2TcLDD11660602 上述表达形式上的变换，不改变总势能的大小！

上式简写为：ppD1DTKDDTR 20，得到节点平衡方程KDR，即：

应用驻值条件：

81700333816800u1333u2AE181481cL2u3L333336u8168400u33351870033302 262考虑到u10，并划去第一个方程，解出其余四个方程得到：

u10u472cL3u388， u48AE1174128u5由此得到的位移场在一般位置上均为近似值（小于精确解）。 单元应力由公式：xExEBd 得到。

cL11cL2LLs第一个单元应力为： s

22A6A23cL5cL2LLs第二个单元应力为： s

22A6A2

13

而实际上，该问题的位移精确解是一个三次函数：

cu(12L2xx3)

6AE应力精确解是一个二次函数：

c(4L2x2)

2Au有限元解和精确解σ有限元解精确解o2Lxo2Lx

图2.5.2弹性杆里兹解结果比较

图2.5.2将上述两个有限元解答与精确解作了比较，由图中可以看出，位移的有限元解与精确解一致度很高，应力的有限元解略差，也在可接受范围内。由于在求得结构的节点位移后，要通过导数的运算来求单元的应变和应力，导致了精度的下降，因此位移的近似程度比应力的近似程度更好。

14

3 上机实验

3.1 实验一 3.1.1 实验题目

图示一个简支梁平面应力模型。梁截面为矩形，高度h=160mm，长度L=1000mm，厚度t=10mm。上边承受均布压力q =1N/mm2，材料E=206GPa，μ=0.29。X方向正应力弹性力学理论解为：

6qL2yy232x3(x)yq(42)

hh5h4分别应用3节点三角形单元、4节点线性等参元（完全积分、减缩积分、非

协调模式）、8节点二次等参元完全积分进行下列数值实验：1）用较粗单元网格求解梁中部应力分量x的最大值和上下边法向应力分量，并对各单元计算精度进行比较分析；2）对粗网格下梁中部最大位移进行对比和分析；3）通过网格加密对比试验3节点三角形单元和8节点二次等参元的收敛速度。总结出研究结论，撰写实验报告。（10分）

3.1.2 实验目的

通过实验了解单元网格形状的选取以及网格的粗细对计算精度的影响，通过比较各种网格的计算结果与精确解考察有限元解的收敛性，并了解简支梁受均布载荷时的应力状况。 3.1.3 建模概述

（1） 启动ABAQUS/CAE。 （2） 建模

由于问题具有对称性，建模时只需以其中的一个平面来建模。建立长度为1m，高度为0.16m的平面四边形区域。 （3） 设定材料属性

设定材料参数E=206GPa，0.29。

（4） 在Assembly模块中定义装配件为Independent。 （5） 设置分析步。

（6） 定义边界条件和施加载荷

① 先将梁的左右边界从中点分为两段； ② 在上边界施加沿y方向向下的均布载荷；

15

③ 定义简支梁的左边界条件，约束x和y方向位移；

④ 定义简支梁的右边界条件，约束y方向位移。 （7） 划分网格 粗网格为6x10，细网格为10x100.分别建立3节点三角形单元（Tri，linear），4节点线性等参元（完全积分Quad，Linear；减缩积分Quad，linear，Reduced integration；非协调模式Quad，Linear，Incompatible modes）和8节点二次等参元（Quad，Quadratic）。 （8） 提交计算 （9） 查看结果

3.1.4 计算结果分析与结论 （1）理论解

X方向正应力由下式计算：

6qL2yy232x3(x)yq(42)

hh5h4已知q=1N/mm2,h=160mm,L=1000mm，ymaxh代入上式得 2xmax61061113600.0810429.497MPa 0.1634245

（2）粗网格（610）计算结果如下： 3节点三角形单元：

图3.1. 1粗网格三角形单元应力云图

图3.1.2 梁中部应力分量x变化曲线

16

图3.1.3 相应数值

图3.1.4三角形单元上边法向应力

图3.1.5 相应数值

图3.1.6 三角形单元下边法向应力

17

图3.1.7 相应数据

4节点线性等参元： 完全积分：

图3.1.8 4节点全积分单元应力云图

图3.1.9 4节点全积分梁中部应力分量x变化曲线

18

图3.1.10 相应数值

图3.1.11 4节点全积分单元上边法向应力

图3.1.12 相应数值

19

图3.1.13 4节点全积分单元下边法向应力

图3.1.14 相应数值

减缩积分：

图3.1.15 4节点减缩积分单元应力云图

20

图3.1.16 4节点减缩积分梁中部应力分量x变化曲线

图3.1.17 相应数值

图3.1.18 4节点减缩积分单元上边法向应力

21

图3.1.19 相应数值

图3.1.20 四边形减缩积分单元下边法向应力

图3.1.21 相应数值

22

非协调模式：

图3.1.22 非协调元应力云图

图3.1.23 非协调元轴中部应力分量x变化曲线

图3.1.24 相应数值

23

图3.1.25 非协调元上边法向应力

图3.1.26 相应数值

图3.1.27 非协调元下边法向应力

24

图3.1.28 相应数值

8节点二次等参积分元完全积分：

图3.1.29 8节点二次等参元完全积分应力图

图3.1.30 8节点轴中部应力分量x变化曲线

25

图3.1.31 相应数据

图3.1.32 8节点上边法线应力变化曲线

图3.1.33 相应数据

26

图3.1.34 8节点下边法线应力变化曲线

图3.1.35 相应数据

将上述计算结果制作成表格如下：

表3-1-1 梁中部的最大正向应力 3节点三 角形单元4节点线性等参元 8节点二次等参元理论值 完全积分 完全积分 减缩积分 非协调模式 完全积分 梁中部应力分量x的最大值（Mpa） 17.0 25.4 24.8 28.9 29.7 29.5

从表中可以看出，3节点三角形单元完全积分精确性最差，产生的结果没有意义。

4节点线性等参单元完全积分的线性单元出现了剪力自锁，剪力自锁引起单元在弯曲时过于刚硬，造成应变偏小，实际结果偏差仍较大。

27

4节点线性等参单元完减缩积分单元存在沙漏问题（零能模式），使得结构过于柔软，变形偏大，甚至在在粗划网格情况下，产生无意义的结果。与精确解相比，仍有一定的误差。

8节点二次等参元完全积分的精度最高。

下面讨论上下边法向应力分量:

图3.1.36三角形单元梁上边法向应力分量 图3.1.37三角形单元梁下边法向应力分量

图3.1.38 4节点完全积分梁上边法向应力分量 图3.1.39 4节点完全积分梁下边法向应力分量

图3.1.40 4节点减缩积分梁上边法向应力分量 图3.1.41 4节点减缩积分梁下边法向应力分量

28

图3.1.42 4节点非协调元梁上边法向应力分量 图3.1.43 4节点非协调元梁下边法向应力分量

法向应力分量 法向应力分量

图3.1.44 8节点二次等参元完全积分梁上边 图3.1.45 8节点二次等参元完全积分梁下边

比较以上有限元结果，在较粗单元网格划分情况下，以上几种求解方法得到

的结果基本一致，梁的上边界法向应力分量分布呈凹形，梁的下边界法向应力分量分布呈凸形，但是与弹性力学求得的理论解还是有差距，理论解表明，梁的上边界法向应力分量等于-1，梁的下边界法向应力分量为0。所以，在较粗单元网格划分情况下的误差还是很大的。

粗网格下梁中部最大位移：

图3.1.46 3节点梁中部最大位移0.130386 图3.1.47 4节点全积分梁中部最大位移0.168795

29

为0.2263 为0.193362

图3.1.48 4节点减缩积分梁中部最大位移 图3.1.49 4节点非协调元梁中部最大位移

图3.1.50 8节点二次完全积分梁中部最大位移

为0.197421

将上述计算结果制作成表格如下：

表3-1-2 梁中部最大位移 3节点三角形单元完全积分 4节点线性等参元 完全积分 减缩积分 8节点二次等参元非协调模式 完全积分 梁中部最大位移（mm）

0.130386 0.168795 0.2263 0.193362 0.197421 由以上数据可知8节点二次等参元完全积分和4节点线性等参单元非协调元的计算结果相当，而4节点线性等参单元减缩积分存在严重的沙漏，变形严重失真，完全不能反映实际情况。其次还可发现4节点线性等参单元完全积分计算精度比三角形单元完全积分高。

（3）以下通过网格加密对比试验3节点三角形单元和8节点二次等参元的收敛速度。

密网格（10100）情况下应力分析结果如下：

30

图3.1.51 密网格三角形单元应力云图

图3.1.52 密网格三角形单元计算时间是7s

图3.1.53密网格8节点二次等参单元应力云图

图3.1.53 密网格8节点二次等参单元计算时间为7s

31

粗网格（610）情况下应力分析结果如下：

图3.1.54 粗网格3节点三角形单元计算时间为5s

图3.1.55 粗网格8节点二次等参元计算时间为5s

密网格三角形单元计算时间与密网格8节点二次等参单元计算时间同为7s，粗网格三角形单元计算时间与粗网格8节点二次等参单元计算时间同为5s，由此可知：相同网格密度情况下，不同单元的收敛速度基本一致，网格加密后收敛速度变慢。

3.1.5 实验体会与总结

通过此次实验，对三角形单元、四边形完全积分单元、四边形减缩积分单元、四边形非协调元和8节点二次等参元完全积分有了深刻认识。

(1)三节点三角形网格下划分的单元要比4节点单元精度要差； (2)粗网格要比细网格精度要差；

(3) 4节点非协调单元和8节点二次积分等参单元下，粗细网格的精度基本上接近；

(4)4节点线性等参单元减缩积分出现零能模式，求解的位移值失真，总刚奇异；4节点双线性位移模式完全积分单元，因为该单元有剪力自锁，会引起单元弯曲时过于刚硬。

因此，不同的积分单元、积分方法、积分阶次、网格密度会导致计算结果不同。如果单元算法选择不当，除影响到计算精度外，还有可能导致计算的失败，

32

因此在使用时要注意适当的选择合适的计算单元。一般应力分析，优先使用8节点二次等参元和4节点非协调元；避免使用3节点三角形单元，在有弯曲变形的情况下，避免使用4节点双线性位移模式完全积分单元，因为该单元有剪力自锁，会引起单元弯曲时过于刚硬。

3.2 实验二 3.2.1 实验题目

图示一管接头，内壁受均匀压力。自行建立其三维几何模型，运用二次六面体单元对其建模并求解。要求利用对称性，自行设置载荷大小和位移约束条件,并撰写实验报告。（10分）

3.2.2 实验目的

通过实验了解轴对称二次六面体单元，学习建立三维模型，通过比较不同模型的计算结果考察了解空心球受内压时的应力状况。 3.2.3 建模概述

（1）建立图示轴对称模型； （2）创建材料和截面属性；

Steel弹性模量为200000，泊松比0.3。

（3）在Assembly功能模块中创建一个非实体，定义为Independent； （4）在Step功能模块中定义静力分析步； （5）定义边界条件和载荷；

除左边界和上边界未给约束外，其余给予固定约束。内壁施加均匀压力10MPa。

（6）Mesh划分网格；

在Mesh功能模块中，采用Structure分网，划分二次六面体单元。 （7）在Job功能模块中，创建job,然后提交job,监控运行状态，创建提交分析； （8）在Visualization功能模块中进行后处理操作，查看结果并分析。

33

3.2.4 计算结果分析与讨论

图3.2.1 六面体二次单元网格图

图3.2.2 六面体二次单元约束图

34

图3.2.3 六面体二次单元应力图

图3.2.4 相交圆柱体变形图

35

图3.2.5 相交圆柱体最大应力图

结论：在ABAQUS软件下对三维模型运用六面体网格划分，可以较准确的反应实际受力变形的情况。对复杂六面体用二次六面体单元划分网格要多加练习，才能熟练。

3.2.5 实验体会与总结

通过这道题，对于轴对称单元有了深刻的理解，理解了如何简化模型，了解到了对称简化在有限元中的应用，特别是对于三维轴对称问题，在网格划分比较密的情况下，计算机求解时间较长，如果我们巧妙地运用轴对称，将大大简化计算量和计算时间，同时熟悉如何在ABAQUS软件画三维模型，深刻理解了三维二次六面体单元，掌握了其中的许多ABAQUS软件的实用知识。

3.3 实验三 3.3.1 实验题目

一个矩形平板，长1000mm,宽350mm，厚度8mm。材料的E=200GPa，0.3， 7.8210-6Kg/mm3。板的一对短边简支。在较粗单元网格下，分别用8节点六面体完全积分等参元、8节点六面体非协调单元和20节点六面体完全积分等参元计算其前六阶自由振动频率和振型，对计算结果列表作对比，并进行分析。撰写实验报告。（10分）

3.3.2 实验目的

（1）通过对矩形平板的模态分析和瞬态响应分析，了解有限单元法中有关动力学问题分析和求解的过程和方法，并学习动力学问题的知识；

（2）学习如何利用有限元分析软件ABAQUS进行结构的模态分析和瞬态动力学分析；

36

（3）学习如何利用ABAQUS软件进行有限元的动态载荷施加等。掌握结构的模态分析方法，会从计算结果中提取频率和振型。掌握结构的动力学分析方法，会定义输出步；

（4）理解线性和二次六面体全积分等参元的模态计算的影响，以及分别使用隐式直接积分和振型叠加法对瞬态响应计算的异同。 3.3.3 建模概述 （1）建模：

矩形平板，长1000mm,宽350mm，厚度8mm；

（2）设置材料属性，在Property功能模块中定义材料的属性：

E200GP,0.3，7.820109tmm3；

（3）定义装配件为Independent； （4）设置分析步；

（5）定义边界条件和载荷，在Load功能模块中，定义实体简支约束； （6）Mesh划分网格；

分别划分网格为620，别建立8节点六面体完全积分等参元、8节点六面体非协调单元和20节点六面体完全积分等参元；

（7）在Job功能模块中，创建job,然后提交job,监控运行状态； （8）在Visualization功能模块中进行后处理操作； （9）查看结果并分析。

3.3.4 计算结果分析与结论

（1）8节点六面体完全积分等参元：

图3.3.1 8节点六面体完全积分等参元一阶振型和振频

37

图3.3.2 8节点六面体完全积分等参元二阶振型和振频

图3.3.3 8节点六面体完全积分等参元三阶振型和振频

图3.3.4 8节点六面体完全积分等参元四阶振型和振频

图3.3.3 8节点六面体完全积分等参元三阶振型和振频

38

图3.3.5 8节点六面体完全积分等参元五阶振型和振频

图3.3.6 8节点六面体完全积分等参元六阶振型和振频

（2）8节点六面体非协调单元：

图3.3.7 8节点六面体非协调单元一阶振型和振频

39

图3.3.8 8节点六面体非协调单元二阶振型和振频

图3.3.9 8节点六面体非协调单元三阶振型和振频

图3.3.10 8节点六面体非协调单元四阶振型和振频

40

图3.3.11 8节点六面体非协调单元五阶振型和振频

图3.3.12 8节点六面体非协调单元六阶振型和振频

（3）20节点六面体完全积分等参元：

图3.3.13 20节点六面体完全积分等参元一阶振型和振频

41

图3.3.14 20节点六面体完全积分等参元二阶振型和振频

图3.3.15 20节点六面体完全积分等参元三阶振型和振频

图3.3.16 20节点六面体完全积分等参元四阶振型和振频

42

图3.3.17 20节点六面体完全积分等参元五阶振型和振频

图3.3.18 20节点六面体完全积分等参元六阶振型和振频

表3-3-1 自由振动频率

振型 8节点六面体完全积分等参元频率（HZ） 8节点六面体非协调单元频率（HZ） 20节点六面体完全积分等参元频率（HZ）

43

一阶 二阶 三阶 四阶 五阶 六阶 72.285 101.32 288. 316.12 8.85 665.07 18.499 75.105 75.377 168.45 171.22 291.81 18.469 74.794 76.368 166.59 170.14 281.36

由以上计算结果可知， 8节点六面体非协调单元和20节点六面体完全积分等参元计算其前六阶自由振动频率和振型的精度较高，误差较小；8节点六面体完全积分等参元计算其前六阶自由振动频率和振型的精度则较低，误差较大。 3.3.5 实验体会与总结

通过本实验比较了不同类型单元对矩形平板自由振动频率和振型的影响。通过计算可以知道，高阶单元比低阶单元计算精度稍高一些。对物体的自由振动频率和振型有了深刻认识，掌握了用有限元软件进行模态分析和动态响应分析的方法，尤其是对材料的动态和静态响应特性有了深刻的体会，材料的动态力学特性和静态力学特性有着很大差异，因而在考查冲击过载下材料的强度时，不能简单的用静态强度来代替。

44