第三章有限单元法

第3章 有限单元法

在工程技术领域内，工程师常常运用数学和力学的知识将实际问题抽象成它们应遵循的基本方程（常微分方程或偏微分方程）和相应的边界条件。对于大多数的工程技术问题，由于物体的几何形状和载荷作用方式是很复杂的，除了方程性质比较简单且几何边界相当规则的少数问题之外，试图按经典的弹性力学和塑性力学方法获得解析解是十分困难的，甚至是不可能的。为了克服这种困难，有两条解决途径：一是引入简化假设，将方程和边界条件简化为能够处理的问题，从而得到它在简化状态下的解答。这种方法只在有限的情况下可行，因为过多的简化将可能导致不正确的甚至错误的答案。另一条解决途径就是数值解法，如有限差分法、边界元法、有限单元法和离散元法等。对于非线性问题，有限单元法更为有效，且已经出现了许多通用程序。

有限单元法的主要优点是：①建立于严格理论基础上的可靠性。因为用于建立有限元方程的变分原理或加权余量法在数学上已被证明是微分方程和边界条件的等效积分形式。只要原问题的数学模型是正确的，同时用来求解有限元方程的算法是稳定、可靠的，如果单元满足收敛准则，则近似解最后收敛于原数学模型的精确解；②适应性强，应用范围广，不仅能成功地分析具有复杂边界条件、非线性、非均质材料、动力学等难题，而且还可以推广到解答数学方程中的其它边值问题，如热传导、电磁场、流体力学等问题；③适合计算机实现的高效性。由于有限元分析的各个步骤可以表达成规范化的矩阵形式，最后导致求解方程可以统一为标准的矩阵代数问题，特别适合计算机的编程和执行。已经出现了许多大型结构分析通用程序，如：NASTRAN、ASKA、ADINA、ANSYS、ABAQUS等，可以直接应用。这些优点使有限单元法得到了广泛的应用和发展。

3.1有限单元法分析的基本步骤

在工程或物理问题的数学模型（基本变量、基本方程、求解域和边界条件等）确定以后，有限单元法作为对其进行分析的数值计算方法的基本步骤如下：

(1) 离散化

一个复杂的弹性体可以看成是由无限个质点组成的连续体，它具有无限个自由度。将一个受外力作用的连续弹性体离散成一定数量的有限个小单元的集合体，单元之间只在节点上互相联系，亦即只有节点才能传递力。因此，该集合体只具有有限个自由度，这就为解算提供了可能。由无限个质点的连续体转化为有限个单元的集合体的过程，称为离散化。在数学意义上说，就是把微分方程的连续形式转化为代数方程组，以便于进行数值求解。

将求解域离散为有限单元，根据基本场变量与坐标的关系而决定采用一维、二维和三维单元。一维单元用线段表示，二维单元可为三角形元或四边形元，三维单元常用四面体或六面体元。单元划分越密，计算精度越高，但计算工作量也越大。

(2) 单元分析

51

根据弹性力学的基本方程和变分原理建立单元节点力和节点位移之间的关系，形成单元有限元方程。

1）确定插值函数（形函数）

有限单元法将整个求解域离散为一系列仅靠公共节点联接的单元，而每一单元本身却视为光滑的连续体。单元内任一点的场变量（如位移），可根据其在单元中的假定分布规律（插值函数）Ne 由本单元的节点值插值求得。

N称为单元的形函数矩阵，它与单元节点坐标、节点数目（即单元形状）及插值形式有

关。形函数矩阵分量的数目应与单元节点自由度数相等。

所谓“形函数”是一种视单元节点量为“已知量”的插值函数，由单元节点的坐标可直接求出形函数的具体形式。

2）建立单元方程

当问题比较简单时，可以直接根据问题的物理概念建立单元方程。不过，在一般情况下，特别是二维和三维单元，直接法会显得过于繁杂而难以应用。为此，需要采用更为一般的数学方法，如变分法、加权余量法或虚功原理。

(3) 计入边界条件，求解有限元方程

组集后的总体特征矩阵（或称为总刚度矩阵）是奇异的，必须计入边界条件才能求得唯一解。计入边界条件的方法有三种：

1）直接代入法

将自由度的已知量（边界条件）从总体方程组中消去，得到一组阶数降低了的修正方程，其原理是按节点位移已知和待定重新组合方程。由于这种方法改变了方程组的阶数，使程序编制复杂化，故程序中一般不采用。

2）对角线元素置1法

当给定位移值是零位移时，例如无移动的铰支座，可以将系数矩阵K中与零节点位移相对应行列的主对角元素改为1，其它元素改为0；再将载荷列阵中与零节点位移相对应的元素改为0即可。

这种计入边界条件的方法简单，不改变原方程组的阶数和未知量顺序，但只适用于边界条件为零值的情况。

3）对角元素乘大数法

当有节点位移为给定值ajj时，第j个方程作如下修改：对角元素Kjj乘以大数（可取1010左右量级），并将Pj用Kjjj取代，即：

52

K11K21Kj1Kn1

K12K22K

j2





K

jj



Kn2



K1na1p1

K2na2p2

 （3-1） KjnajKjjj

Knnanpn

经过修改后的第j个方程为：

Kj1a1K

j2

a2KjjajKjnanKjjj （3-2）

对于多个给定位移则按序将每个给定位移都作上述修正，得到全部进行修正后的K和P，然后解方程则可得到包括给定位移在内的全部节点位移值。这个方法使用简单，对任何给定位移（零值或非零值）都适用。采用这种方法引入强制边界条件时方程阶数不变，节点位移顺序不变，编制程序十分方便，因此在有限单元法中经常采用。

(4) 后处理计算

根据解方程组求得的节点基本场变量（如位移等）计算其它相关量，如应变、应力等，视具体问题而定。

3.2一维问题的有限单元法

下面用一个例子来介绍有限单元法解一维问题的过程。

例3.1 用有限单元法求图3-1a所示受拉阶梯杆的位移和应力。已知杆截面面积A（1）=2×10-4m2，A（2)= 1×10-4m2，各段杆长L(1)=L（2)=0.1m，材料弹性模量E(1)=E(2)= 2×10-5Mpa,作用于杆端的拉力F3=100N。

(a)示意图 (b)有限元模型 (c)单元图

图3-1 受拉阶梯杆

解题过程如下： (1) 单元划分

根据材料力学的平面假设，等截面受拉杆的同一横截面的不同点可认为具有相同的位移和应力，即位移只与截面的轴向坐标有关，所以可将阶梯杆看作由两个“一维单元”组成，同一单元内截面积及材料特性不变。每一单元有两个节点，它们分别位于单元两端，相邻两单元靠公共节点联接，用线段表示。图3-1a中所示的受拉阶梯杆简化为由两个一维单元和三个节点构成的有限元模型，如图3-1b所示，图中①②为单元号，1、2、3是节点号。取节点

53

位移为基本未知量，应力由求得的节点位移算出。

(2) 单元分析

1)确定单元插值函数（形函数） 本例中，每单元有两个节点，可采用线性插值方式。图3-1c为一典型单元图，两节点号分别为i和j，水平方向坐标为X，设单元中坐标为x处节点的位移为

e

x，根据线性插值关系分别记为：

n



式中 n——单元节点数；

Ni——形函数，N1

12

e

x

i1

e

NiiN

1,

2l

N2

12

1；

——单元内节点的自然坐标，

l——单元长度；

x

xc；

xc——单元中心点的总体坐标，xc

x1xn

2

。

2)建立单元方程 由等截面杆变形与拉力的关系（胡克定律）得到：

AeEeeeeeijPil

ee

AEeePe

jijel

式中，Pi和Pj分别为作用于单元e的节点i和节点j的节点力，写成矩阵形式为：

ee

AE1

e

l1

ee

1iPi

ee

P1jj

ee

简记为：

KαP

e

e

e

式中 Ke——单元特性矩阵，在力学问题中常称为单元刚度矩阵；

e

P——单元节点力列阵，PPi

e

e

Pj

e



T

。

e

到目前为止，单元方程尚不能求解，因为节点力列阵Pe尚属未知。其分量Pi和Pj是相邻单元作用于单元e的节点i和j的力，即属于单元之间的作用力，只有将具有公共节点的单元“组集”在一起才能确定上述节点力和节点外载荷之间的关系。

（3） 整体分析

建立总体方程组，为获得总体方程组，必须先将单元方程按照局部自由度（i，j）

54

e

和总体自由度（1、2、3）的对应关系进行扩展。具体来说，单元1的扩展方程为：

110

110

1

01P1



102P2 

003

AEl

1

11

式中，各项上角码表示单元序号，下角码表示自由度总体序号。 单元2的扩展方程为：

A

2

E

2

l

2

0

00

011

010

212P2 

213P3

由于相邻单元公共节点上的位移相同，所以可将扩展后的各单元方程相加。将上两式相加得到：

A1E1

1

l

A1E11

l0

AEl

11

AElA

21

11

2

1

AE

2

lE

2

2



A

2

E

2

lA

2

2

2



E

2

ll

2



P112P2 P33

上述组集过程可记为：

N

N

e

K

e1

α

P

e1

e

式中 N——有限元模型的单元数。

组集后的结果简记为：

KαP

式中 K——总体刚度矩阵；

P——总体节点载荷列阵。

需要指出的是，对单元的一个公共节点而言，除了有相邻单元作用于该节点的力之外，还可能有作用于该节点的外载荷。若一节点上无外载荷作用，如本例中的节点2，则说明各相邻单元作用于该节点的力是平衡的，即该节点的节点合力为零。若某节点上有外载荷作用，如本例中的节点3，则各单元作用于该节点的内力和与该节点的外载荷P3相平衡，即：

22

A2E2AEP30 2322ll

这就是说，列阵P各分量的含义是作用于相应自由度（节点位移）上的节点外载荷。将

55

相应数据代入得：

46

104

0

462

01P1

220 

123

上式即为本题的总体线性方程组，但不能获得唯一解，因为上式中的矩阵是奇异的。这种奇异性是因为总体方程组式只考虑了力平衡条件，而只根据力平衡不能唯一地确定系统的位移，因为系统在任意刚性位移的情况下仍可处于力平衡状态。为获得各节点位移的唯一解，必须消除可能产生的刚体位移，即必须计入位移边界条件。

（4） 计入边界条件，解方程组

本题的位移边界条件为10，那么，式中只剩下两个待求的自由度2和3。也就是说，可从式中消去一个方程。比如，舍去第一方程并将10代入后得到：

66

10

2

220 231

解得2=0.25×10-6m，3=0.75×106 m，这与材料力学法求得的结果相同。

（5） 计算单元应变和应力

由材料力学得知，单元任一点的应变ε

ε

e

e

e

x与位移αx的关系为：

x

dα

e

x

dx

Bα

e

式中 B——应变矩阵。

应力—应变关系为:

ζ

e

x

Eε

ee

x

EBα

ee

对于单元1



对于单元2



1

E

1

12

l

1



21000.2510

5

6

0.1

0.5MPa

2

E

2

21

l

2



2100.250.7510

5

6

0.1

1MPa

56

3.3平面问题的有限单元法

本节以弹性力学中的平面问题为例，介绍有限单元法如何解决平面问题。根据弹性体的几何形状和受力状态，平面问题包括平面应变问题和平面应力问题。

平面应变问题的特点是：物体沿某一坐标方向（如z方向）的尺寸远大于沿其它两个坐标轴的尺寸；垂直于Z轴各截面的形状和尺寸均相同；所有外力与Z轴垂直，且不随z坐标变化；物体的约束条件不随z坐标而改变。在这种情况下，物体远离两端的各截面，可认为没有z向位移，而沿x和y方向的位移对各截面均相同，与z坐标无关，各截面内将产生平面应变。z方向的应变z0，三维应力状态下，只包含有x,y,xy,zxy, yzxz0,而且应力分量均为x，y的函数。

平面应力问题的特点是：物体沿某一坐标方向（如z方向）的尺寸远小于其它两个坐标方向的尺寸；外力沿周边作用且与xy平面平行，且体积力也垂直于z轴；由于物体在z方向厚度很小，故外载的表面力和体积力都可看成是沿z向不变化的。约束条件在xy平面内。在这种情况下，所有应力都发生在xy平面内，沿z方向无应力，即z0,yzzx0。故平面应力状态的应力分量为x 、y 、xy，且都是x，y的函数。

上述两种平面状态的平衡方程、几何方程、边界条件和变形连续方程都是相同的，唯有物理方程有差别。如果把平面应力问题中的E换成E12，换成问题的求解就相同了。

平面问题的有限元计算过程同一维问题： (1) 结构的离散化

在弹性力学平面问题中用以进行离散化的单元形式有很多种。如三角形三节点单元、四边形（矩形）四节点单元、三角形六节点单元、四边形八节点单元等，如图

3-2。

1，则两种

图3-2 典型的二维单元

57

(2) 单元分析

1)确定形函数 以二维问题的三节点三角形单元为例，如图3-3所示，每个节点有2个位移分量。

i

uivi

i,j,m （3-1）

每个单元有6个节点位移即6个节点自由度： ie

ju

i

m

viujvjum

vm （3-2）

T

3节点三角形单元位移模式选取一次多项式：

u12x3yv45x6y

（3-3）

单元内的位移是坐标x，y的线性函数。1～6是待定系数，称之为广义坐标。6个广义坐标可由单元的6个节点位移来表示。在（3-3）式的1式中，代入节点i的坐标xi，yi可得到节点i在x方向的位移ui，同理可得：

ui12xi3yi



uj12xj3yj （3-4） um12xm3ym

解（3-4）式可以得到广义坐标由节点位移表示的表达式。系数行列式为：

Dxixjxm

yi

yj2A （3-5） ym

A——三角形单元的面积。 广义坐标1～6为

58

1

1D

uiujum2

1D

3

1D



1

aiuiajujamumxjyj2A

xmym



uiyi

1

biuibjujbmumujyj （3-6） 2Aumym

xiui

1

ciuicjujcmumxjuj2A

xmum

xi

yi

同理，利用3个节点y方向的位移，即（3-3）式的2式可求得

456

在（3-6）和（3-7）中，

ai

xjxmavaviijjmm

2A



1

bivibjvjamvm （3-7） 2A1

civicjvjcmvm2A

1

av

yjymyjymxjxm

xjymxmyj

bi

yjym

cixjxm





i,j,m （3-8）



式（3-8）i,j,m表示下标轮换，即ij，jm，mi。写成矩阵形式为

uiuj0um



Nmvi

vjvm

uNi

u

v0

0Ni

N0

j

0N

j

Nm0

Ni

N

j

ee

（3-9） NmαNα

将求得的广义坐标1～6代入（3-3）式，将位移函数表示为节点位移的函数：

59

uNiuiNjujNmum

 （3-10）

vNiviNjvjNmum

其中，Ni为当i节点有单位位移且j、m节点固定时所引起的单元内部的位移分布函数

Ni

12A

aibixciy

i,j,m （3-11）

单元面积A为

A

12D

12

a

i

ajam

12

bc

i

j

bjci （3-12）

2)建立单元有限元方程 确定了单元位移后，可以利用几何方程和物理方程求得单元的应变和应力。

物体在平面内的状态有两种描述方式：一是给出各点在x方向的位移u和在y方向的位移v，二是给出各点微小矩形单元的应变x，y，xy。其中x和y分别表示沿x方向和y方向的线应变，xy表示切应变（角应变）。

经过弹性物体内任意一点P，沿x轴和y轴的方向取两个微小长度的线段PAdx和

PBdy，如图3-4，假定弹性物体受力以后，P、A、B三点分别移动到P、A、B。

'

'

'

首先求线段PA和PB的正应变，即x和y，用位移分量来表示。设P点在x方向的位移是u，则A点在x方向的位移，由于x坐标的改变，将是u变是：

udxuux

dx

ux

dx，可见线段PA的正应

x

（3-13）

在这里，由于位移微小，y方向的位移v所引起PA的伸缩量是高一阶微量，因此略去不计。同样可见，线段PB的正应变是

vy

y

（3-14）

然后求出线段PA和PB之间直角的改变，也就是剪应变xy，用位移

分量来表示。由图可见，这个剪应变是由两部分组成的：一部分是由y方向的位移v引起的，即x方向的线段PA的转角；另一部分是由x方向的位移u引起的，即y方向的线段PB

的

60

转角。

设P点在y方向的位移分量是v，则A点在y方向的位移分量将是v线段PA的转角是

v

vdxv

vx

vx

dx。因此，



同样可得线段PB的转角是

dx

（3-15）



uy

（3-16）

于是可见，PA和PB之间直角的改变（以减小时为正）也就是剪应变xy

xy

vxvy



uy

（3-17）

综合式（3-13）、（3-14）、（3-17），得到平面问题由位移求应变的几何方程：

x

ux

,y

,

xy



vx



uy

（3-18）

由上式可见，当物体的位移分量完全确定时，应变分量可完全确定。反之，当应变分量完全确定时，位移分量却不能完全确定。因为当物体发生一定的应变时，由于约束条件的不同，它可能具有不同的刚体位移，因而它的位移并不是完全确定的。在几何方程中，位移用（3-10）代入，得到单元应变：

x

e

εyNαNi

xy

Bi

Bj

N

j

e

Nmα

ee

BmαBα （3-19）

式中 B——应变矩阵，分块矩阵为



xBi0

y



0Ni

y0x

N

i

x00Ni

Niy



0Ni

y

Ni

x

i,j,m （3-20）

求导可得到：

Nix



12A

bi，

61

Niv



12A

ci （3-21）

代入（3-20）式得到：

bi1Bi0

2Aci

0

ci

bi

i,j,m （3-22）

3节点单元的应变矩阵是：

bi1

Bm2A0

ci

0cibi

bj0cj

0cjbj

bm0cm

0

cm （3-23） bm

BBi

Bj

式（3-23）中bi，bj，bm，ci，cj，cm由式（3-8）确定，它们是单元形状的参数。当单元的节点坐标确定后，这些参数都是常量（与坐标变量x，y无关），因此B是常量阵。

单元应力可以根据物理方程求得。在完全弹性的各向同性物体内，应变分量与应力分量之间的关系根据胡克定律：

x

1E1



x



y







yxy （3-25） E21xyxy

E

此外，

z



E



x



y

 （3-26）

可由应变求应力，其物理方程为：



x



E1E1

E

22



x



y



E

2

yyx

xy

21



xy

1





 （3-27） 1

xy

2

为了便于以后的应用，上式可写成矩阵形式：

x

ζy,

xy

62

x

εy （3-28）

xy

1E

D2

1

0



10

0



0 （3-29） 1

2

D为弹性矩阵。它是一个对称矩阵，其元素只与弹性常数E和有关。

在物理方程中代入（3-27）式可以得到：

xee

ζyDεDBαSα （3-30）

xy

式中 S——应力矩阵，SDBDBi

Bj

BmSi

v0ci

ci



10

bi2

Sj

Sm ，分块矩阵

SiDBi

E0

21

2



bi

0biA

10

ci

2

i,j,m （3-31）

式中E0，0——材料常数，对于平面应力问题，E0E,0，对于平面应变问题，E0

E1

2

,0

1

。

与应变矩阵B相同，应力矩阵S也是常量阵，即3节点单元中各点的应力是相同的，在很多情况下，不单独定义应力矩阵S，而直接用DB进行应力计算。

单元刚度矩阵Ke

KiiT

BDBtAKji

Kmi

KijK

jj

K

e

Kmj

Kim



Kjm （3-32） Kmm

式中，t是二维体厚度。代入弹性矩阵D和应变矩阵B后，它的任一分块矩阵可表示成：

KrsBrDBstA

T

E0t41v

2



K1AK2K3K4

r,si,j,m （3-33）

其中：

63

1v0

Kbbcrcs1rs21v0K2v0crbsbrcs2 （3-34） 



Kb1v0

3v0rcs2crbs





K4crcs1v02brbs

对于离散模型，系统位能是各单元位能的和，离散模型的总位能为：

1

p



αeT2BTDBtdxdyαee









α

e



eT





N

T

ftdxdy

eT

NTTtdS

α



e

seij



式中，f是作用在二维体内的体积力；T是作用在二维体边界上的面积力。

T

BDBtdxdyK

e

T



N

ftdxdyPe

f



NT

TtdSPe

SeS

ij

单元等效节点载荷矩阵Pe为：

PePee

bPs 将单元节点位移αe用结构节点位移α表示：

αe

Gα 其中

αu1

v1

u2

v2



ui

vi



um

vmT

122i12i2m12m2j12j2n001000000000100000G0

00000100

000000010000010000

00





1





0

式中 G——单元节点转换矩阵；

n——结构的节点数。

（3） 整体分析

离散形式的总位能可表示为

64

3-35）

3-36）

3-37）

3-38）

3-39）

3-40） （（（（（ （



p

α

T

1

G2

T

KGαα

e

T

G

T

P

e

 （3-41）

Te

KGKG

（3-42） Te

PGP

式（3-41）可以写作



p



12

αKααP （3-43）

TT

离散形式的总位能p的未知变量是结构的节点位移α，根据变分原理，泛函p取驻值的条件是它的一次变分为零，即

pα

0 （3-44）

得到有限元的求解方程是：

KαP （3-45）

可以看出结构整体刚度矩阵K和结构节点载荷列阵P都是由单元刚度矩阵Ke单元等效节点载荷矩阵Pe集合而成。

例3.2 如图3-5a所示一悬臂梁，梁的宽度为1m，长度为2m，厚度为t，板外端受垂直向下的外载荷如图，设泊松比13，弹性模量为E。用有限单元法求出节点位移和单

元上的应力。 解题步骤如下：

（1） 离散化

这是一个弹性力学平面应力问题，用3节点三角形单元离散该结构。离散化的单元划分及单元编号如图3-5b所示，节点采用两种编号，一是节点整体编号，四个节点统一编号为1、2、3、4，在建立节点的平衡方程时，按照这个整体节点编号顺序建立；二是节点的局部编号，规定每个单元的三个节点按反时针方向顺序，各自编号为i、j、m用来进行单元分析。单元与节点整体编号关系如图c所示。作用在结构上的载荷和支座反力，在结构离散化后，都要作用在节点上，因此，悬臂梁固定端支承面上的节点1和4均可简化为铰支座。

65

（2） 单元分析 1) 求形函数 单元Ⅰ：

x2

N

0

0x2

y0

0y

x1y

2



 x1y

2

单元Ⅱ：

x

1

20

01

x2

1y0

01y

x1y

2



 x1y

2

N

Ⅱ

2) 写单元方程 单元应变矩阵B

1

b0

cm0

bm

0

001b

001a

01a0

1b



01 a1b

bi

1B0

2A

ci

0cibi

bj0cj

0cjbj

bm0cm

01a

1E

D2

1

0



10

0

0 1

2

1bE

SDB2

1b

0

e

T

T

0012b

0012a

a1a0





1b

b2a

1



a

1 a

1



2b



KBDBtABStA

66



abt2



E1

2

1b0001b0

0001a0012b12ab012ab12b

22

1a

01b11bab001a1b

012ab12a012a12ab

22

0012b

0012a

a1a0





1b



b12a

a1



a1

2b





abEt21

2

1

b200ab1b2ab



ab001





1

2

b

12ab12a

2

a

2

ab1a

2

ab

11

22

2ab

1ab2ab









ab

1



2ab122b

1

2a1 ab2ab



1122a2b



将a=1m，b=2m代入式中，得到

1

400EtⅠ

K2

1

2142



 98

181401418

1201214

1

2

32414

1

Ⅰ单元和Ⅱ单元的单元刚度矩阵相等，KⅠ为对称矩阵，且有KⅠKⅡ。

（3） 求整体刚度矩阵K

67

324

14

1

4

1Et4

K2

1

00

12

2

1498



14



14180



0012



001411498



12140





2



240

11401418098



2800

32

198



24



1214



12140



2

3214

144

1140



240



141



28

32

14

11





（4） 引入边界条件

支承条件共有4个，即：

u1v1u4v40

对于节点位移为零的约束，修改主对角线元素改为1，其它元素改为0，即修改K中的1、2、7、8各行和各列，修改后的单元刚度方程如下：

1

0000000

01000000

00K33K43K53K6300

00K34K44K54K6400

00K35K45K55K6500

00K36K46K56K6600

00000010

0u10

00v1

0u20P0v2

2 

0u30P0v320u40

1v40

将整体刚度矩阵K代入上式，得：

32

40Et2

11

2

14

0988



12





24



2

3214

1

1

04

u

2P1v22

 

1u30



P4v3

92

8

68

将μ=1/3 代入，解得：

u21.50

P8.42v2

 u3Et1.88v38.99

结构的节点位移向量为：

δu1

PEt

v1u20

v21.50

u3v3u4

v4 8.99

0

T

T

0

8.421.88

（5）求单元应力

由节点位移求单元应力的公式为：

ζSδ

e

e

e

1bES2

1b

0

0012b

0012a

a1a0





1b

b2a

1



b

1 a

1



2b



对于单元Ⅰ，有：

1

29E1S

860

1

2P9E1ζ

Et86

0

0016

1310

0013

121613

1310



121613

1

31 1

6

0016

0013



11.508.4234.500.84403PP11.500.281 t0168.421.5801060

对于单元Ⅱ：

69

S

Ⅱ



E21

2





bi

biA

1

ci

2

ci

ci12

bi

bjcj

cj

cj

12

bj

bm

bj

12

bm

12

cm

cm

cm

 1bm2



E21

2



aaA0

001223

0012232 13

bbb

a

a

a

12

b

b



1

a2

b



19E1



163

0

0013

0023

113

23

20

ζ

Ⅱ

3P3E1Et16

0

001

002

260

312

0



24.500.844

P1.503P61.500.289 8.4216t2.230.41811.888.99

3.4 空间问题的有限单元法

本章将要讨论的空间问题包括空间杆系结构问题和三维弹性问题。求解空间问题的有限元法，它的基本概念、基本理论、建立计算格式的基本过程以及实施步骤与平面问题完全相

同，勿需重述。这里只着重研究其单元和单元刚度矩阵，并就几种常用单元来阐明空间问题的有限元法。

空间问题的有限元分析中常用单元的单元包括四面体单元、六面体单元，杆单元和梁单元。以四面体单元为例来说明求解空间单元刚度矩阵的步骤，4节点四面体单元如图3-6所示。 (1) 位移函数

从单元的自由度数可知位移函数应取线性

形式：

u12x3y4z

v56x7y8z （3-46） wxyz

9101112

在i、j、m、p四个节点，分别有

ui12xi3yi4ziuxyzj12j3j4ju m12xm3ym4zmup12xp3yp4zp

根据克莱姆法则，由式（3-47）求出1，2，3，4，代入式（3-46），得

u

1

6V

aibixciydizuiajbjxcjydjzuj

ambmxcmydmzumapbpxcpydpzup

其中

xiyiziV

1xjyjzj6x mymz mxp

y

p

zp

是四面体的体积，系数ai，bi，ci，di为

xj

yjzjaixm

ymzmi，j，m，p

xpypzpyjzjbiymzmi，j，m，p

ypzpxj

1zjcixm

1zmi，j，m，p

xp

1

zp

71

（3-47）

（3-48）

（3-49）

xj

dixm

xp

yjymyp

111

i，j，m，p （3-50）

式中的i，j，m，p表示下标轮换，即ij，jm，mp，pi。

为了使四面体的体积V不为负值，单元的四个顶点的标号i，j，m，p必须遵循右手螺旋规则，即右手螺旋按顺序转动时，大拇指应指向p。

用同样的方法，可得出其余两个位移分量

v

aibixciydizviajbjxcjydjzvj6V

ambmxcmydmzvmapbpxcpydpzvp

w

aibixciydizwiajbjxcjydjzwj6V

ambmxcmydmzwmapbpxcpydpzwp

11

（3-51）

（3-52）

结合表达式（3-48）、（3-51）、（3-52），可以把位移分量表示成

u =u =INi

v

wNα

e

INjINm

e

INpα （3-53）

其中，I是三阶的单位阵，各插值函数为

Niaibixciydiz6VNjajbjxcjydjz6V

i，m

j，p

（3-54）

(2) 四面体单元的应力矩阵及刚度矩阵

弹性体在一般变形情况下，有3个方向的线应变x、y、z及3对剪应变xyyx，

zy、zxxz。由弹性力学可知，应变和位移间的几何关系用矩阵表示为

yz

72

x0

xy0zε

xy

yyzzx0

z

0y0xz0



00

uz

v （3-55） 0wyx

将四面体单元的位移表达式（3-48）、（3-51）、（3-52）带入（3-55），可以得到单元的应变为

e

εBαBi

BjBm

e

（3-56） Bpα

应变矩阵B的分块矩阵

bi

010Bi

6Vci

0di

0ci0bidi0

0

0di0cibi

i，j，m，p （3-57）

单元的应力用节点位移表示为

ζDεDBα （3-58）

e

令应力矩阵SDB，其分块矩阵为

bi

Ab1iA1biE1

Si

6112VA2ci

0

A2di

A1ciciA1ciA2biA2di0

A1di

A1di

di

0A2ci

A2bi

i，j，m，p （3-59）

其中

A1

e

1

，A2

1221

设节点虚位移为α，代入公式（3-55）得

73

e

εBα

由虚功方程可得

αPα

e

e

eT

B

V

T

DBαdxdydz （3-60）

e

可以得出单元载荷矩阵

eTeee

PBDBVK （3-61）

其中，单元刚度矩阵Ke

kiikji

kmikpi

kijkjjkmjkpj

kimkjmkmmkpm

kip

kjp

 kmp

kpp

K

e

其中

brbsA2crcsdrds



A1crbsA2brcs

36112V

A1drbsA2brds

E1A1brdsA2drbs

A1brcsA2crbs

crcsA2brbsdrdsA1drcsA2crds

krs



A1crdsA2drcs



drdsA2brbscrcs

ri，j，m，p；si，j，m，p

可以得到结构的整体平衡方程

KαP

由平衡方程解出节点位移后，即可根据每个单元上的节点位移αe求得该单元中的应力。 空间问题的有限元分析比较复杂，因为空间问题通常自由度数较多，构成的整体空间计算模型的自由度数和节点数非常大。另外，离散连续体结构时，单元划分也很困难。体单元的空间离散网格，除比较简单的情况外，很难像平面问题那样直观清晰地划分出来，特别是体内单元，往往只能想象。因此，应用有限元法虽然可以建立空间问题的统一数值计算格式，但如果不借助计算机求解，这种计算格式几乎没有实际意义。

3.5 有限元分析计算机程序

有限单元法是通过计算机为工具实现数值求解的，现有软件国际上早在60年代初就投入大量的人力和物力开发具有强大功能的数值模拟分析程序。从那时到现在,世界各地的研究机构和大学发展了一批专用或通用数值模拟分析软件。

NASTRAN是工业标准的FEA原代码程序及国际合作和国际招标中工程分析与校验的首选工具,可以解决各类结构的强度、刚度、屈曲、模态、动力学、热力学、非线性、声学、流体 结构耦合、气动弹性、超单元、惯性释放及结构优化等问题。能求解高度非线性、瞬态动力学、流体及液固耦合等问题，其先进的技术可解决广泛复杂的工程问题，如金属成形、爆炸、

74

碰撞、搁浅、冲击、穿透、安全气囊(带)、液固耦合、晃动、安全防护等。可用于零部件的初始裂纹分析、裂纹扩展分析、应力寿命分析、焊接寿命分析、随机振动寿命分析、整体寿命预估分析、疲劳优化设计等各种分析。同时该软件还拥有丰富的与疲劳断裂有关的材料库、疲劳载荷和时间历程库等，分析的最终结果具有可视化特点。可以处理各种线性和非线性结构分析，包括线性/非线性静力分析、模态分析、简谐响应分析、频谱分析、随机振动分析、动力响应分析、自动的静/动力接触、屈曲/失稳、失效和破坏分析等。可以解决各种高度复杂的结构非线性、动力、耦合场及材料等工程问题，尤其适用于冶金、核能、橡胶等领域。

NASTRAN的前后处理软件PATRAN是世界公认最好的集几何访问、数值模拟建模、分析及数据可视化于一体的新一代框架式软件系统。通过其全新的“并行工程概念”和无可比拟的工程应用模块可将世界著名的CAD/数值模拟/CAM/CAT(测试)软件系统及用户自编程序自然地融为一体。

ANSYS是融结构、流体、电场、磁场、声场分析于一体的大型通用数值模拟分析软件。它能与多数CAD软件实现数据的共享和交换,如I DEAS、UG、Pro/Engineer、NASTRAN等，软件主要包括3个部分:前处理模块、分析计算模块和后处理模块。前处理模块提供了一个强大的实体建模及网格划分工具，用户可以方便地构造数值模拟模型；分析计算模块包括结构分析(可进行线性分析、非线性分析和高度非线性分析)、流体动力学分析、电磁场分析、声场分析、压电分析以及多物理场的耦合分析，可模拟多种物理介质的相互作用，具有灵敏度分析及优化分析能力；后处理模块可将计算结果以彩色等值线显示、梯度显示、矢量显示、粒子流迹显示、立体切片显示、透明及半透明显示(可看到结构内部)等图形方式显示出来，也可将计算结果以图表、曲线形式显示或输出。

ALGOR作为世界著名的大型通用工程仿真软件，广泛应用于各个行业的设计、数值模拟分析、机械运动仿真中。它包括静力、动力、流体、热传导、电磁场、管道工艺流程设计等,能够帮助设计分析人员预测和检验在真实状态下的各种情况，快速、低成本地完成更安全更可靠的设计项目。ALGOR以其分析功能齐全、使用操作简便和对硬件的要求低，在从事设计、分析的科技工作者中享有盛誉。作为中高档数值模拟分析工具的代表之一，ALGOR在汽车、电子、航空航天、医学、日用品生产、军事、电力系统、石油、大型建筑以及微电子机械系统等诸多领域中均有广泛应用。工程师们通过使用ALGOR进行设计，虚拟测试和性能分析，缩短了产品投入市场的时间，并能以更低的成本制造出优质而可靠的产品。

ABAQUS被广泛地认为是功能最强的数值模拟软件之一，可以分析复杂的固体力学和结构力学系统，特别是能够驾驭非常庞大复杂的问题和模拟高度非线性问题。ABAQUS不但可以做单一零件的力学和多物理场的分析，同时还可以做系统级的分析和研究。优秀的分析能力和模拟复杂系统的可靠性使得ABAQUS被各国的工业和研究机构广泛采用，其产品也在大量的高科技产品研究中发挥着巨大的作用。

ABAQUS是一套功能强大的工程模拟数值模拟软件，其解决问题的范围从相对简单的线性分析到许多复杂的非线性问题。ABAQUS包括一个丰富的、可模拟任意几何形状的单元库，并拥有各种类型的材料模型库，可以模拟典型工程材料的性能，其中包括金属、橡胶、高分子材料、复合材料、钢筋混凝土、可压缩超弹性泡沫材料以及土壤和岩石等地质材料。作为通用的模拟工具，ABAQUS除了能解决大量的结构(应力/位移)问题，还可以模拟其它工程领域

75

的许多问题，例如热传导、质量扩散、热电耦合分析、声学分析、岩土力学分析(流体渗透/应力耦合分析)及压电介质分析。

下面以通用有限元软件ANSYS为例，介绍其求解问题的过程。 解决问题的流程一般分为前处理、求解和后处理三大部分：

（1） 在前处理过程中建立有限元模型，输入几何模型的基本资料，如节点坐标、单元长度等；定义单元类型及实常数；设定材料属性，如弹性模量、泊松比等；对几何模型划分单元。

（2） 有限元分析通过求解器完成，其主要步骤包括输入载荷情况、设定边界条件、设定求解选项，然后进行有限元求解计算。

（3） 后处理可以察看计算结果，比如在post1后处理器中，可以查看求解结果的变形情况、振型状态或各种等值线（如等效应力线），还可以绘制求解结果的各种分析曲线。

例3.3 一悬臂梁，长L=1m，梁截面为矩形，截面宽W=0.05m，高H=0.01m，其左端固支，右端受集中载荷P=100N，如图3-7。材料为钢，弹性模量E=2.1×5Mpa，泊松系数μ=0.27。分析其最大位移以及变形后的应力。

解题过程如下： (1) 前处理过程 1)进入ANSYS

由于ANSYS在进行一个分析时，会在工作目录中生成和保存许多文件，为方便管理，最好对每一个分析对象单独使用相应的工作目录和文件名字。

Ansys 8.1 →Interactive →改变工作路径（ working directory）至用户自设的路径 →输入工作文件名字（Initial jobname），如beam（必须是英文）→运行（Run）。

2)设置计算类型

主菜单（ANSYS Main Menu）： 前处理（Preferences） →选择结构分析（Structural）→ OK

3)选择单元类型

悬臂梁定义单元类型为最常使用的2节点三维线性梁单元BEAM188单元。Preprocessor →Element Type→Add/Edit/Delete… →Add… →选择单元类型为2节点梁（Beam 2 node 188）→ OK

4)定义材料参数 需要指定的材料参数为弹性模量Ex和泊松比μ。

76

Preprocessor →Material Props →Material Models →Structural →Linear →Elastic →Isotropic →输入 EX:2.1e11, PRXY:0.27 → OK

5)定义截面 对于梁单元，需要指定其截面形状和尺寸。

Preprocessor →Sections →Beam →Common Sections →定义矩形截面：截面编号ID=1，宽度B=0.05，高度H=0.01 →Apply→OK

6)生成几何模型

首先生成关键点。梁的模型简化为一条直线，需要生成3个关键点，其中两个用来生成梁，另一个为表示截面方向的方向点。

Preprocessor →Modeling →Create →Keypoints →In Active CS →依次输入三个点的坐标：1(0,0,0)，2(1,0,0)，3(0,0.05,0) →OK（点1、2为梁的两个端点，点3是为了生成梁所建立的方向点）。

生成关键点后，再生成梁。

Preprocessor →Modeling →Create →Lines →lines →Straight line →连接两个特征点，1(0,0,0),2(1,0,0) →OK

7)网格划分

首先指定表示梁的直线的属性，将单元类型和材料模型赋予直线，然后指定其网格划分份数为10份，通过网格划分，得到有限元模型。

Preprocessor →Meshing →Mesh Attributes →Picked lines选取上一步所画的线 →OK →Pick Orientation Keypoint(s)：选择YES→选择方向点，拾取点3（0,0.05,0） →OK→网格划分工具（Mesh Tool）→Size Controls中Lines： Set →Pick All →输入划分份数（NDIV）：10→OK (回到网格划分工具 Mesh Tool 窗口) → Mesh →Pick All → Close

这时单元看上去还是一条直线，在菜单PlotCtrls选择Style，打开Size and Shape 对话框，在[/ESHAPE]选项后面的复选框点选on，[/REPLOT]中选择Replot ，用截面方式显示单元的实体形状。

（2） 求解过程

77

经过前处理的悬臂梁有限元模型已经建立，进入求解器进行求解。

1)模型施加约束

首先是最左端节点加约束，悬臂梁的左端为固定约束，所以将节点1的所有自由度全部约束。图中以三角形标志表示约束。

主菜单（ANSYS Main Menu）： Solution →Define Loads →Apply →Structural →Displacement → On Nodes →选取节点1（0,0,0）→ OK → 约束所有自由（ALL DOF）→ OK

然后对最右端节点施加y方向的载荷。

Solution →Define Loads →Apply →Structural →Force/Moment→ On Nodes →选取节点2（1,0,0）→ OK →选择FY方向，在VALUE后输入100→ OK

施加约束及载荷后的有限元模型如图3-8所示。

图3-8 施加载荷和约束的有限元模型

2)分析计算 在求解前应对数据进行保存，此问题是一个简单的线性问题，可以不进行特殊的求解控制设置，直接求解。

Solution →Solve →Current LS →OK （3） 察看结果

进入后处理器post1查看计算结果，分别查看变形、位移和应力分布结果。 1)查看变形情况 由图3-9可以看到受力前后的变形比较

主菜单（ANSYS Main Menu）： General Postproc →Plot Results →Deformed Shape → 选择 Def + undeformed →OK

图3-9 受力后的变形（与变形前比较）

78

2)查看节点位移 由图3-10中下方的标尺条可以看到，最大位移出现在最后端，下垂了0.039923m。

General Postproc →Plot Results →Contour Plot →Nodal Solu →选择：DOF Solution, →Y-Component of displacement，在Undisplaced shape key 中选择 Deformed shape with undeformed model→

OK

图3-10 梁节点的位移分布图

图3-11 梁等效应力分布图

3)查看等效应力 查看米塞斯等效应力云图，从图3-11中可以看到，最大应力出现在梁的根部，最大等效应力为114Mpa。

General Postproc →Plot Results →Contour Plot →Nodal Solu →选择：查看等效应力Stress →von Mises stress ，在Undisplaced shape key 中选择 Deformed shape with undeformed model→OK

求解的命令流代码： /PREP7 ！进入前处理器

ET，1，BEAM188 ！定义单元类型 MP，EX，1，2.1E11 ！定义杨氏弹性模量

79

MP，PRXY，0.27 ！定义泊松比

SECTYPE，1，BEAM，RECT，，0 ！定义梁的截面 SECOFFSET，CENT SECDATA，0.05，0.01 ！生成关键点 K，1，0，0，0 K，2，1，0，0 K，3，0，0.5，0 L，1，2 ！生成线

LSEL，，，，1 ！定义单元属性及方向点 LATT，1，，，，3，，1 ！生成单元 ESIZE，10 LMESH，1 /ESHAPE，1.0 EPLOT FINISH /SOLU

D，1，ALL ！施加位移约束 F，2，FY，-100 ！施加集中力 SOLVE ！求解 FINISH

/POST1 ！进入后处理器

PLDISP，1 ！绘制结构变形图，并与变形前比较

PLNSOL，U，Y，1 ！用等高线图显示位移结果，并与变形前比较

80

PLNSOL， S，EQV， 1 ！用等高线图显示米塞斯等效应力结果，并与变形前比较 例3.4 空间桁架结构内力计算

空间桁架结构形式、几何尺寸、及所受载荷如图3-12所示，求桁架各节点的位移及各杆件所受力的大小。

图3-12 空间桁架结构

单元类型选用LINK8单元（空间桁架单元），由各点坐标生成节点，并采用直接法建立有限元模型，即由节点直接生成单元。用命令流方式求解，求解结束后，在后处理器中分别察看桁架的结构变形（与受力前作比较）、变形的动态过程，并用绘制位移的等高线图，列表显示各节点的位移及节点反力。

命令流如下：

/PREP7 ！进入前处理器 ET,1,LINK8 ！定义单元类型 R,1,50e-4 ！定义实常数

MP,EX,1,2.1E11 ！定义杨氏弹性模量 MP,PRXY,0.3 ！定义泊松比 ！生成节点 N，1，1，0，1.5 N，2，-1，0，1.5 N，3，0，0，-1.5

81

N，4，0，5，0 ！连接节点生成单元 E，1，4 E，2，4 E，3，4 FINISH

/SOLU ！进入求解器

！在节点1，2，3处施加位移约束 D，1，ALL D，2，ALL D，3，ALL ！施加载荷 F，4，FX，-20E3 F，4，FY，-50E3 F，4，FZ，20E3 SOLVE ！求解 FINISH

/POST1 ！进入通用后处理器

PLDISP，1 ！绘制桁架的变形图，并于变形前作比较 ANDSCL，10，0.5 ！制作结构的变形动画

PLNSOL， U，SUM， 0，1 ！用等高线图显示位移结果 PRNSOL，DOF ！列表显示各节点的位移 PRESOL，FORC ！列表显示各节点的反力

3.6 实例分析

3.6.1球磨机回转体有限元分析

82

球磨机是在选矿、建筑、化学等工业部门中重要的粉磨设备，其主体是一个水平装在两个大型轴承上的低速回转的筒体，某溢流型球磨机结构如图3-13所示，筒体内径为5500mm，筒体长度为8570mm，转速为0.2283r/s，

主电机功率为4500kW。电动机通过减速机及周边大齿轮减速传动，驱动回转部回转。简体内部装有适当的磨矿介质-钢球。磨矿介质在离心力和摩擦力的作用下，被提升到一定的高度，呈抛落或泄落状态。物料由给料口连续地进入简体内部，被运动的磨矿介质所粉碎，并通过溢流和连续给料的力量将物料排出机外，以进行下一段工序处理。

球磨机筒体的强度直接影响到球磨机的使用寿命和传动装置各部件的运动精度。过去一直沿用的经验强度计算方法，是将球磨机筒体简化为简支梁，然后按平面弯曲和扭转组合变形来计算。显然，由于筒体受力十分复杂，这种算法不能准确和全面地反映其应力和变形规律，因此导致磨机设计存在一定程度上的盲目性，而采用有限元方法可以很好的解决此问题。

（1）有限元模型

选择三维实体单元Solid45划分网格。Solid45单元有八个节点，每个节点有x、y、z三个方向的移动自由度，该单元有塑性、蠕变、膨胀、应力刚化、大变形和大应变的性能，有限元模型如图3-14所示。

（2）求解过程 1)工作载荷的确定

大齿轮的输入转矩T＝3.13×10N•mm，由计算得到作用在大齿轮上的切向力Ft＝7.278×105N，轴向力Fa＝6.6876×104N，径向力Fr＝2.6601×104N。

5

9

球磨机满载正常工作时物料重量为3.4×105kg，钢球重量为6×104kg，所以其正常工作时介质总重为4×10kg。但是由于在工作状态时有一些物料及钢球等在离心力和摩擦力的作用下被提升到一定的高度后甩在空中，呈抛落或泄落状态，并不作用在筒体上。因此应计算出直接作用在筒体上的物料及钢球的重量。

由已知的球磨机筒体主要技术参数：筒体有效半径R＝2.6425×10mm，长度L＝1×

634

10mm，转速n0.2283r/s，破碎介质的松散比重4.6510kgmm，可得：

3

83

筒体角速度：2n/601.44rad/s； 最外层脱离角：1cos

1

Rn

2



56.75；

900

1

最小半径脱离角：2cos式中，R1为最小半径，

R1

250n

2

R1n900

2

73.88



1.3310mm。

3

由此可得进行圆周运动的破碎介质的重量等计算所需积分限为：

0



2





19056.7533.25 29073.8816.12









1



2

从而得出圆周运动的破碎介质的重量、离心力及质心坐标分别为

G1

Lg2

4

2

sin22cos2

2

3

01

2.12822610kg

3

5

C

4Lg9

4

3coscos3sin

5

01

9.456610N



sin

4

23

0

xc

2g

sin

1

01

6



2

sin22cos2

23

3

3

1.610mm116

14

0

3

yc

2g

(cossin)sin4

01



1



2

sin22cos2

820mm

由上述计算结果可以看出物料及钢球总计4×105kg，重量中只有重量G1＝2.128×105kg（Fn＝9.2×105N，Ft＝1.92×106N）直接作用在筒体上，其离心力C＝9.46×105N，质心坐标Xc＝1.6×103mm，Yc＝812mm，56.75，73.88，如图3-15所示。





2)质量等效方法

84

在实体建模时没有考虑筒体衬板、中空轴内部的衬套及固定环等，因此在施加重力载荷时，把衬板的重量按等效密度加到筒体及端盖上，把中空轴内部的衬套和固定环等的重量按等效密度加到中空轴上。等效的方法如下：

①筒体（包括对应的衬板）重量为1.47×10kg，体积是9.4704×109mm3，所以等效密度是1.5522×10-9kg/mm3；

②每个端盖（包括对应的衬板）重4.65×10kg，体积是3.377×10mm，故等效密度是1.377×10kg/mm；

③中空轴（包括固定环、进料衬套、出料衬套等）重9.4×10kg，体积是6.528×10mm，故其等效密度是1.44×10-5 kg/mm3。

已知大齿轮的重量为5.2×104kg，所以可以得出球磨机的结构重量为3.86×105kg，它们的重力作用由施加的重力加速度g＝9.8m/s2确定。

3）分析工况

上面简单介绍了与筒体计算关系密切的主要载荷的计算方法。根据筒体的工作特点，选定满载静止、正常工作和最大启动转矩时三个工况进行有限元分析。两端按简支方式施加约束，其载荷及约束形式如图3-16所示。

(3) 察看结果 1)应力

筒体在三种工况下的当量应力分布如3-17所示。

4

9

3

-5

3

9

3

45

85

(a)满载静止工况应力云图 (b)满载静止工况筒体内侧轴向应力曲线

(c)正常工况应力云图 (d)启动工况应力云图

图3-17 回转体应力分布图

从图3-17中可以清楚的看出，满载静止时大部分应力在11MPa左右，只有在筒体与端盖连接处的应力突然变大，且在进料端筒体处的应力要大一些，筒体最大应力为47.8MPa。正常工况，最大应力值为35.3MPa。其位置在大齿轮的边缘载荷加载处，在筒体与右端盖（无齿轮一侧）的连接处的应力值为24.9MPa，与左端盖连接处的应力值为18MPa，其它位置的应力值很小。启动工况最大应力值为106MPa，其位置在大齿轮的边缘载荷加载处附近。

筒体与端盖的材料屈服极限为300MPa，大齿轮的材料屈服极限是540MPa，各工况的最大应力值都小于相应材料的屈服极限，故在各工况下球磨机回转体强度均满足要求。

2)位移分析

(a)满载静止工况位移分布云图 (b)满载静止工况轴向位移分布曲线

86

(c)正常工况位移分布云图 (d)启动工况位移分布云图

图3-18 回转筒体轴向位移云图

从图3-18中可以清楚的看出，满载静止时筒体的最大位移为1.682mm，正常工况下最大位移为1.194mm，启动工况最大位移值为3.103mm，其位置均在载荷作用方向筒体中心处。其余绝大部分变形量都很小，说明各工况下回转体的刚度足够。曲线图中曲线突变是由筒体与端盖连接处的结构突变引起的。

3.6.2 装载机倾翻保护结构有限元分析

装载机的倾翻保护结构，简称ROPS。ROPS的性能要求主要有四项，即侧向承载能力、侧向能量吸收能力、垂直承载能力和纵向承载能力。为掌握装载机ROPS的力学特性，对此类ROPS

分析结果与样机试验结果进行了对比。

(1) ROPS结构

ROPS结构简图如图3-19所示，ROPS室设计成一体，ROPS通过螺栓与车架连接。图3-19中的DLV为的容许挠曲极限量（Deflection-limiting Volume称DLV），表示了驾驶室中的人体极限生存空间，求变形后的ROPS不得侵入DLV。

(2) ROPS性能仿真 1）有限元仿真模型

ROPS性能仿真模型如图3-20。ROPS中的细长构件采用梁单元模拟。加强板筋和底板采用壳单元模拟。连接螺栓用三维弹塑性梁单元模拟。ROPS材料的单调加载屈服极限y=260MPa，车架材料分屈服极限y=345MPa。

87

2）边界条件

图3-20 ROPS有限元模型

根据ISO3471标准要求，仿真时各工况载荷条件见下表3-1

（表中整机质量m=6000kg），载荷作用位置和方向可参见图3-20。

表3-1 载荷条件【ISO3471】

仿真工况

载荷

FL=6m=36.0kN

侧向承载能力 侧向能量吸收 垂直承载能力

指向X轴正向，作用于左上纵梁 E=12500(m/10000)1.25=6601J FV=19.61m=117.7kN

Z轴负向，作用于ROPS左右上纵梁

纵向承载能力

FZ=4.8m=28.8kN

指向Y轴正向，作用于后上横梁中间

约束施加在前后车架的铰接点和后桥中心线处，如图3-20所示的约束符号（架呈固定状态。

(3) 结果对比与分析 1）侧向加载

），使车

侧向加载结束时的仿真与试验变形对比见图3-21a和图3-21b。侧向加载结束后为50mm；侧向模拟地平面LSGP与DLV间的距离为101mm，但都未侵入DLV。

88

(a)试验变形

图3-21 侧向加载结束时ROPS变形对比

侧向载荷与加载点挠度的特性曲线如图3-22所示。由该图可见：①侧向承载力达到36.0kN的标准要求时，挠度为90mm。载荷与挠度基本上呈直线状态，表明ROPS仍处于线弹性变形阶段。此阶段内ROPS吸收的能量为1795J；②当载荷增加到50.5kN、挠度为198mm时，侧向能量吸收达到了6601J的标准要求。此试验

阶段，随着侧向载荷增加，载荷—挠度曲线的斜率逐渐变小，表明ROPS开始屈服。 2) 垂直加载

侧向加载结束后，在侧向变形的基础上，进行ROPS垂直加载试验，垂直载荷与加载点的挠度特性曲线如图3-23所示。由该图可见垂直载荷为117.7kN时，有限元仿真的挠度值为3.5mm，而测试值为4.3mm，但二者都呈直线趋势，即表明ROPS垂直变形在弹性阶段。经分析，该工况下ROPS变形很小，没有任何ROPS构件

侵入DLV，故承载能力和变形量都满足国际标准ISO3471要求。

89

图3-23 垂直载荷－挠度特性曲线 图3-24 纵向载荷－挠度特性曲线

3）纵向加载

在侧向和垂直加载试验的残余变形基础上，进行纵向加载试验，纵向载荷与载荷作用处的挠度特性曲线如图3-24所示。由该图可见仿真曲线与测试曲线吻合的较好，二者都呈斜线上升，表明ROPS在纵向变形也在弹性阶段。纵向载荷满足标准要求的28.8kN时，ROPS纵向挠度为43mm，且ROPS的任何构件未侵入DLV，纵向加载试验满足ISO3471要求。

90