“有限元法原理及应用”讲义-2013秋修订(2)

图1-4 弹性杆件

势能泛函为：

L

lE2

,xAdx ∫qudx

02

Πp=

∫

该问题中，在求出近似位移场u=u(x)后，单轴应力σx=Eεx=Eu,x。 最简单的“可能”位移（试探函数）为：u

=a1x，待定参量a1称为广义坐标。

它满足位移边界条件，不满足应力边界条件。将它代入上式势能泛函得到：

AEL2cL3

Πp=a1 a1

23

cL2

由驻值条件 =0，得到 a1=

da13AE

dΠp

cL2cL2

因此， u=。 x 和 σx=

3AE3A

考虑一个更好的试探函数：

u=a1x+a2x2，广义坐标为a1，a2 上式代入势能泛函，并应用驻值条件求解得到：

a1 cL =

a2 12AE

u=

7L

，从而 3

cLcL

(7L 6x) (7Lx 3x2)， σx=12A12AE

图1-5将上述两个解答与精确解作了比较，显然，位移的近似程度比应力的近似程度更好！

图 1-5 弹性杆里兹解结果比较

实际上，该问题的位移精确解是一个三次多项式函数（自行推导）： u=

c

(3L2x x3) 6AE

因此，如果我们的试探函数取到三次项，所得到的里兹解一定是该精确解。为什么？如果试探函数多项式取到四次或更高次，会得到什么解答？

前面研究了经典里兹法解弹性体变形和应力的原理和过程。可以总结出该方法的重要特点：

1) 在求解域整体上假定位移场（试探函数）；

2) 假定的位移场必须是可能位移（或称为许可位移，即满足连续性和边界几何约束条

件）和简单的。

3) 要得到收敛解，试探函数必须是完备的。

4) 里兹解往往过刚，除非位移试探函数包含了精确解。由于假定的位移模式往往给结

构加上了约束，使结构不能按其要求的方式自由变 形，从而刚化了结构。

第二章 弹性力学有限元法基本原理（一）

第一节 里兹法的有限元形式

由于经典里兹法求解时需要在求解域整体上假定位移场，且位移场必须满足连续性和边界位移约束条件（许可位移），因此经典里兹法在解决实际问题时，尤其是复杂几何形状的二、三维问题，具有很大局限性。

解决的办法是在求解区域上分片假设位移场。

下面以一维直杆的分析为例子，研究基于里兹法的有限元位移法基本原理和求解过程。考虑图2-1（a）所示的结构，长度改为3L，把杆分为三个部分。 应用里兹法求解时分三个区域假设位移场：

u=b1+b2x 0≤x≤L u=b3+b4x L≤x≤2L u=b5+b6x 2L≤x≤3L

上述多项式系数bi（i=1,6） 作为待定参量。

为了使得上述假设位移是“可能的”，上述待定参量必须满足一定约束关

系，因此该问题的独立参量（广义坐标）只有3个。

图2-1 （a）截面积A，弹性模量E，轴向受力杆。

（b）杆的有限单元

把上述假设的位移代入势能泛函：

Πp=∫

3L

3LE2

,xAdx ∫qudx

02

进行分段积分，再应用势能驻值条件，立即就求出待定参量，位移场完全确定，进而可求出应力。

上述是原理上的里兹法有限元形式。关键是它采用分片多项式拟合全域上的可能位移场，这就是有限元法的实质。

但是，上述分片形式的假设位移场有下列缺点：

（1） 必须对其表达式进行调整，使其满足连续条件和边界约束条件； （2） 多项式系数广义坐标bi缺乏明显的物理意义。

因此，上述不是通常意义上的标准有限元形式，具有局限性。必须把bi换成节点的未知位移分量，得到标准有限元形式。

图2-1（b）给出了连续结构的三个分区（单元），单元之间的连接点称为节点。这一步骤称为离散化。

下面对每个单元，在单元内假设简单的位移场（位移模式）：根据单元2个节点上的位移未知量进行线性插值得到单元上任意点的位移。过程如下。

对单元1按线性函数插值方法得到：

u=

sL s

u1+u2=N1u1+N2u2，矩阵形式为u= N {d}

LL

L s

=N Ls

=[N1 L

N2] 是插值基函数矩阵，称为“形函数”矩阵；N1，N2

分别是单元两个节点的位移形状函数，简称为“形函数”。

{d}=

u1

。 ，为单元1的节点位移列阵（单元自由度）

u 2

u

u3 u u4

对其它两个单元也可同样获得插值函数形式的假定线性位移场： 单元2：u

= N {d}，{d}= 2 = N {d}，{d}= 3

单元3：u

N 同上。每个单元中，位移都是线性场，但数值不同，取决于单元两端节点位移。显然，整个杆上，由各单元上假设的位移场拼接而成的位移试探函数是连续的，只要满足位移边界条件u1

=0，得到的就是全域可能位移场。这样的位移场已经把节点

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

下面在上述分片位移插值试探函数的基础上进一步实施里兹法求解。 在三个单元上分片进行总势能计算：Πp

=∫

3L

3LE2

xAdx ∫qudx

02

首先计算应变：

εx=u,x=u,s

d 11

单元内有：εx=——单元应变矩阵 N {d}= B {d}， B = ds LL

所以一个单元内的应变能为：

L

U=∫

或者U=

121LT1TLT

εxAds=∫εxAEεxds={d} ∫ B AE B ds {d}

0 2202

1T1 1

{d}[k]{d}，[k]=AE ——杆单元刚度矩阵 2 11L

这里由于三个单元的几何和物理特性完全相同，其单元刚度矩阵相同。 杆上载荷与上节例题相同，q

=cx。则载荷在三个单元内分别表达为：

q=cs

，

q=c(L+s)， q=c(2L+s)

cL2T 1

{d} N csds=d}

6 2

T

T

代入外力功积分式，对三个单元分别计算外力功。第一个单元外力功为：

∫

L

quds=∫uqds=∫

L

T

L

第2、3单元外力功分别为：

7 cL2cL2T 4

{d} 和 {d}T 66 5 8

总势能是三个单元总应变能减去三个单元外力功。为了能够以矩阵形式相加，将单元势能矩阵表达式中的单元节点位移列阵{d}用结构整体位移向量{D}=[u1

u2u3

u4]代

T

替，单元刚度矩阵[k]等与单元有关的矩（列）阵扩展成结构整体规模（4×4），则相加后系统总势能表达为：

1T AEΠp={D}

L2

2

T cL {D}

6

1 1 11

00

00

0000

0 000

01 10 AE +

0 L 0 11

0 000 0 0 7

8

0 0

00 AE +

0 L 0

0 0

000

000 {D}

01 1

0 11

1

2

cL2 + 0 6 0 0

4

cL2 + 5 6 0

上述表达形式上的变换，不改变总势能的大小！ 上式简写为：

Πp=

1

{D}T[K]{D} {D}T{R} 2

应用驻值条件：

Πp D

=0，得到节点平衡方程[K]{D}={R}，即：

0 u1 1 1 10

AE 12 10 u2 cL2 6

=

6 12 L 0 12 1 u3

0 11 8 u4 0

考虑到u1

=0，并划去第一 …… 此处隐藏：2300字，全部文档内容请下载后查看。喜欢就下载吧 ……