第1章机械系统的理论建模

第1章机械系统的理论建模

目前，对机械结构动态性能的研究主要有三种基本方法，即理论建模及分析方法、实验建模及分析方法和二者相结合的方法。

所谓理论建模及分析方法是基于结构动力学原理，根据结构的设计方案、图样、先验知识和资料等建立起模拟机械结构动力特性的动力学模型，而无需依赖于已有的机械设备。通过对该动力学模型的分析计算，即可获得该机械结构各种模拟的动力特性。这不仅可以检验其动力特性是否满足设计目标，是否需要对结构进行修改，还可以通过对理论模型的计算机仿真，预估结构设计及其改进后的动力特性或对其进行动态优化设计。所以，理论建模及分析方法，可以在机械结构设计方案具体实施之前，建立其动力学模型，利用计算机进行模拟仿真，对各种设计方案反复进行比较、修改，使其动态特性逼近设计目标函数的要求。从而可经济、迅速地达到优化设计的目标，把提高机械结构动态性能的问题解决在方案及图样设计阶段。

理论建模及其分析方法的不足，在于建立能够确切模拟机械结构动力特性的动力学模型较为困难，就目前的各种理论建模方法而言，都存在一些难以确定的因素，如难于对机械结构各种工况下的边界条件考虑和处理得与实际工况完全吻合，也难于把机械结构中各种结合部来的误差，影响了所建动力学模型的模拟精度。所以，提高理论模型对机械结构动态性能的模拟精度，使之满足工程实际的需要，是机械结构动态性能理论分析方法必须解决的首要问题。

目前，对一般机械结构的理论建模，除应慎重使用集中参数法或传递矩阵法建模外，最主要的建模方法是有限元法。因此，本章首先讨论有限元及其建模方法，包括静态有限元法和动力学问题官限单元法。然后讨论在机械结构动力分析中常用的另外两种理论建模及分析方法：集中参数法和传递矩阵法。

有限单元法是20世纪50年代初期根据变分原理发展起来的一种强有力的数值近似解法。该方法以计算机为手段，采用分割近似，进而逼近整体的研究思想求解数学物理问题。目前，有限元法已在许多领域成为分析、解决工程和数学物理问题的有力工具。

由于弹性力学问题静态有限元法的理论与方法均已成熟，有着丰富的书籍和资料，因此，本章只讨论静态有限元方程的建立，以其说明弹性力学和变分法等

通用力学、数学原理与有限元法的关系，并简述静态有限元法的分析过程，直接给出常用单元的位移模式及计算公式。若要对此部分内容作详细了解，可参见有关书籍。本章还将讨论动力学问题有限单元建模方法。顺便指出，本章所讨论的变分法及其近似解法，也是后续内容“机械系统动态性能的最优控制”和“动态优化”的理论基础。

1．1 有限单元法的预备知识

1908年，瑞士科学家里兹首先提出用变分法处理弹性薄板问题，开创了弹性力学中直接求解泛函极值问题的近似解法。后经人们的长期努力，形成了变分问题直接解法中最重要的里兹法。里兹法不从微分方程出发，而是根据某个物理问题建立其泛函表达式，并根据某种泛函驻值条件，直接求解泛函极值的近似解，从而使变分解法具有了重要的实用意义。因为人们发现并在数学上得到证明：在连续介质问题中，许多物理、力学问题既可以转化为微分方程的定解问题，也可以归结为变分极值问题，它们的表达形式不同，但却是等价的。例如，图1-1所示等截面梁的横向弯曲平衡方程可用材料力学方法

建立，也可由变分法求得。

设梁的抗弯刚度为EJ ，受分布载荷q (x )作用产生

下弯变形w (x )，两端固定的边界条件为

w (0)=w'(0)=0，w (l )＝w'(l )＝0 (a) 图1-l 梁的横向弯曲简图

由材料力学知，等截面梁挠曲线的近似微分方程为

d 2w (x )/dx 2＝M (x )/(EI ) (b)

式中，M (x )为梁横截面上的弯矩。若梁横截面上的剪力为Q (x )，则有如下关系

dQ (x )/dx ＝q (x ) ； dM (x )/dx ＝Q (x )

或 d 2M (x )/dx 2＝dQ (x )/dx ＝q (x ) (c) 将式(c)代人式(b)，即可得等截面梁的横向弯曲平衡方程

EId 4w (x )/dx 4- q (x )＝0 (d)

用变分原理导出梁的平衡方程如下：

在梁达到平衡时，梁和载荷作为整体的势能达到最小值。梁的势能是梁在弯曲时所取的弯曲能，其值为 ????? ??=l dx dx w d EI U 02

22121 (e)

另外裁荷q (x )的势能在梁弯曲变形w (x )的位移中降低

()()?=l

dx x w x q U 02 (f) 所以，梁和载荷作为整体时的总势能为

()()?????????-???? ??=-=l

dx x w x q dx w d EI U U U 02222121 (g) 平衡条件为总势能达到最小值，即 δU ＝0 (h)

式(g)的变分为 ()()???????-=l

dx x w x q w dx d dx w d EI U 02222δδδ (i) 通过分部积分，并利用固定端点条件d (δw )/dx │x =0,l ＝0，δw │x =0,l ＝0，得

()0)(044=??????-=?l

dx x w x q dx w d EI U δδ (j) 根据变分法基本原理，即得梁的平衡方程为

EId 4w (x )/dx 4- q (x )＝0 (k)

可见，两种方法获得完全相同的结果。这说明满足微分方程及其边界条件的函数将使泛函取极值或驻值，反过来，使泛函取极值或驻值的函数恰是满足微分方程及其边界条件的解答。所以，就某个物理问题所建立的泛函极值问题的解和就该问题所建立的微分方程边值问题的解是一致的。而且，从求近似解的角度出发，求泛函极值的近似解常常要比求微分方程的近似解更为方便[1，2]。这就为我们实际求解工程及数学物理问题提供了一条重要途径：直接从某个物理问题的泛函变分求其近似解，或者把微分方程的定解问题转化成相应的泛函变分求其近似解。里兹法就属于这一类数值近似解法。

需要特别指出，我们在这里反复强调的是求解泛函极值的近似解，这是因为，从变分法出发，求解泛函极值归结为求解欧拉方程，这又回到了微分方程的求解问题。变分法和欧拉方程代表同一个物理问题，从欧拉方程求解和从变分法求解具有相同的效果。但欧拉方程的求解常常是困难的，甚至是不可能的，而从泛函变分求近似解往往并不困难，这也正是变分法备受重视的原因。里兹法正是从泛函变分求极值近似解的有效方法之一。

目前常用的有限单元法均与泛函变分问题直接解法中的里兹法或伽辽金法有着紧密的联系。里兹法和伽辽金法具有可以直接求出结构解析解的近似表达

式，便于作理论分析等优点，但它们都是以结构整体为研究对象的近似计算方法，必须预先给定待求结构一个变形允许的坐标函数(基函数)，而对于边界形状比较复杂的结构、要想找到合适的坐标函数却非常困难，即使找到，也需用到相当高次的多项式或某些函数的特殊结构。因此，这种方法只能对很简单的结构求解，而对于由板梁组合而成且形状复杂的机械结构，是无能为力的[2，3]。

有限元法则与此不同。虽然有限元法是在求解泛函变分问题的里兹法基础上发展起来的数值近似解法，但有限元法是建立在结构离散并使之单元理想化基础上的近似方法，具有结构上的分割近似性质。因此，有限元法的基础是变分原理和剖分插值。即首先把复杂连续结构离散为有限个形状简单的单元，如梁、三角形、四边形单元等，再对每一个单元选择一组简单的坐标函数，便可方便地将里兹法或伽辽金法应用于这些单元；然后，根据变形协调条件把这些单元重新组合，便可求得整体结构的数值解。

有限元法的数学基础是变分原理，而力学基础则是能量原理。根据弹性力学最小势能原理，变形弹性体在受到外力作用而处于平衡状态时，在很多可能的变形允许曲线中，使总势能为最小的那条曲线是真正的变形曲线。而弹性体的总势能是一个以位移为自变量函数，以坐标x，y，z为自变量的泛函表达式。根据变分原理，这条真正的变形曲线就是泛函的极值曲线。有限元法把复杂连续结构离散为有限个形状简单的单元，因此很容易将这些单元的假设变形曲线设置得很好，利用里兹法求解这些单元，便可获得这些单元真正的变形曲线，再根据各单元变形协调条件，把各单元变形曲线装配成整个弹性体的初始变形曲线，并用里兹法求解，便可获得整个弹性体的真实变形曲线。这就是基于弹性力学能量原理，把连续弹性体的变形求解问题转化为泛函的变分极值问题求里兹解的有限元法的形成过程。

1．1．1变分法简介

有限元法的数学基础是变分法。变分法是研究泛函极值的一种经典方法[4，5]。有限元法是在弹性力学能量原理的基础上发展起来的一种数值计算方法，而弹性体的能量表达式是一个泛函，因此，求解弹性体的能量表达式问题可以归结为求泛函极值的问题。

1．函数与泛函

我们知道，对于变量x 的某一变化域中的每一个x 值，若y 都有一值与之对应，则称y 是x 的函数，记作y =f (x )。

如果对于某一类函数{y (x )}中的每一个函数y (x )，因变量J 都有一个确定的值与之相对应，则称因变量J 是函数y (x )的泛函数，简称泛函，记为J ＝J [y (x )]。因此可以说，函数是变量与变量的关系，泛函是变量与函数的关系。

例如，图1-2所示的平面中有两固定点A (x a ，y a )和B (x b ，y b )，连接这两点的曲线弧长l 是曲线函数y ＝y (x )的泛函。

图l-2 求弧长的变分问题

由弧长的微分知

(dl )2＝(dx )2+(dy )2

dl /dx ＝[1+(dy /dx )2]1/2＝[1+y ′2]1/2

所以 dx y l xb

xa ?+=2'1

A 、

B 两点间的连线可以是各种曲线形式。显然，对不同的曲线y (x )，就有不同的长度l 与之对应。所以，A 、B 间连线的弧长l 是曲线函数y (x )的函数，即泛函。记为

[]dx y L dx y x y J xb

xa xb

xa ??=+=)'('1)(2 式中，L (y')＝(l+y')2。

一般情况下，L 也是x 、y 的函数，因此可写成

dx x y y L J xb

xa ?=),',( (1-1) 两点间的最短弧长是直线y *(x )，见图1-2。即

l min ＝J *＝min J [y (x )]＝J [y *(x )]

我们把式(1-1)这种建立在函数L和变量J之间的关系叫做泛函关系。满足此泛函关系式的任一函数y(x)均称为泛函J[y(x)]的容许函数。

2．泛函的极值与变分

求泛函的极大值或极小值问题均称为变分问题。

对于泛函J[y(x)]，其自变函数y(x)的增量是指两个不同的y(x)值之差y(x)-y0(x)＝Δy。当这种增量很小时称为变分，用δy(x)或δy表示。即

δy(x)＝y(x)- y0(x)

对于自变函数y(x)的变分δy(x)所引起的泛函的增量可表示分[1，4]

ΔJ[y(x)+δy(x)]- J[y(x)]＝L[y(x)，δy(x)]+R[y(x)，δy(x)]max│δy│(1-2) 式中，L[y(x)，δy(x)]是δy的线性连续泛函；R[y(x)，δy(x)]是δy(x)的高阶无穷小项。当max│δy│→0时，R[y(x)，δy(x)]→0，则L[y(x)，δy(x)]就称为泛函J[y(x)的变分。即

δJ＝L[y(x)，δy(x)] (1-3) 由此可见，泛函的变分是泛函增量的线性主部。若泛函有变分，而且增量ΔJ 可用式(1-2)表达时，则称泛函是可微的。

泛函的变分也可定义为

δJ＝?J[y(x)+αδy(x)]/?α│α=0 (1-4) 式中，α为一正的小量。

根据式(1-4)，利用函数的微分法则可方便地进行泛函变分的计算。

变分规则：

设L1和L2是x、y、y'的函数，则有如下变分规则[6]：

(1)δ(L1+L2)＝δL1+δL2

(2)δ(L1L2)＝L1δL2+δL1L2

(3)δ∫L(x，y，y')dx＝∫δL(x，y，y')dx

(4)δy'＝dδy/dx

如果泛函J[y(x)]在任何一条与y0(x)接近的曲线上所取的值不小于J[y0(x)]，即

ΔJ＝J[y(x)]-J[y0(x)]≥0 (1-5) 则称泛函J[y(x)]在y0(x)曲线上达到极小值。反之，若

ΔJ＝J[y(x)]-J[y0(x)]≤0 (1-6)

则称泛函J[y(x)]在y0(x)曲线上达到极大值。

事实上，所谓泛函极值问题就是在泛函的容许函数集中寻求一个容许函数y*(x )，使泛函J [y (x )]在其上的值达到极大或极小[5，7，8]。

根据泛函变分原理，若可做泛函J [y (x )]在函数y 0(x )上达到极值时，则泛函在y ＝y 0(x )上的一阶变分等于零，即

δJ ＝0 (1-7)

考虑(1-1)

?=xb

xa dx x y y L J ),',( 在自变量的取值区间x a ＝a ≤x ≤

x b ＝b 内确定一个函数y (x )，使其

满足边界条件x ＝a ，y ＝α；x ＝b ，

y ＝β并使泛函取极值。 图1-3端点固定的极值

设曲线y (x )如图1-3所示。y*(x )为极值曲线。并使泛函取极值。

设曲线y (x )如图1-3所示。y*(x )为极值曲线。在它的附近有一容许曲线y*(x )+η(x )，其中η(x )＝δy (x )是任意给定的连续可微函数。则

y (x )＝y*(x )+εη(x ) 0≤ε≤1

代表了在y*(x )与y*(x )+η(x )之间所有可能的曲线。式(1-1)即为此问题数学模型的一般形式。因为对变量没有约束，故称为无约束变分问题。

取y (x )的变分为δy (x )，由式(1-4)得泛函J (y )的变分为 ????? ????+??=xb xa dx y y L y y L y J '')(δδδ (1-8) 已知[?L (δy ′)/δy ′]dx ＝{?L [d (d δy )/dx ]/δy ′}dx ，利用分部积分公式对式(1-8)的第二项积分，得

?????+????????-??=??-??+??=xb xa xb xa xb xa xb xa xb

xa y y L ydx y L dx d y L dx y L dx d y y y L ydx y L J δδδδδδ'''' (1-9)

两端状态固定时，y (x a )＝α，y (x b )＝β，因此

δy (x a )＝0，δy (x b )＝0，?L (δy ′)/δy ′│xa xb ＝0

于是式(1-9)变为

????

??

???-??=xb xa ydx y L dx d y L J δδ' (1-10) 当y ＝y*(x )时，J 为最小值，δJ (y*)＝0。

在x a 到x b 区间内，δy (x )是任意的。因此y*(x )应满足下述二阶微分方程：

0'

=??-??y L dx d y L (1-11) 0'

=??xb xa y y L δ (1-12) 式(1-11)就是无约束泛函极值必要条件，由欧拉和拉格朗日两人独立推导得出，简称欧拉(Euler)方程。它与式(1-12)一起，通常就可以求解出y *(x )。可见，欧拉方程把泛函的极值问题转化成了微分方程的边值问题。式(1-12)为边界条件。若待求函数的边界值是已知的，即边界值是固定不变的，如式(1-12)中的y (x a )＝α，y (x b )＝β，则有δy (x a )＝0，δy (x b )＝0，边界条件自然满足。这种边界条件称之为本质边界条件；若待求函数的边界值是变动的，则δy (x a )、δy (x b )均不为零，边界条件必须由式(1-11)和式(1-12)联立确定。这种边界条件称为自然边界条件。这时边界条件亦称为横截条件。

利用以上对泛函极值问题的讨论，很容易推导出多元泛函的极值条件。

设y ＝[y 1 y 2 … y n ]T 为n 维变量，y (x a )＝α，y (x b )＝β。求下列目标泛函的极值轨线。

()()dx x y y y y y y L y y y J n n xb

xa n ;',,',';,,,,,,212121 ?= (1-13) 式中，L 是y i 及其一阶导数y i ′(i =1，2，…，n )的标量函数。

为寻求使目标泛函J (y 1，y 2，…，y n )取极值的必要条件，可令y 1，y 2，…，y n 中之一，例如y i (1≤i ≤n )，进行变分，其余n -1个变量保持不变，或其变分为零。于是n 元泛函极值问题就变为求n 个一元泛函极值的问题，J 取极值的必要条件当然就是满足欧拉方程

n i y L dx d y L i i ,,2,10'

==??-?? (1-14) 或写成下列方程组的形式

()()????

?????====??-??=??-??n i y x y y x y y L dx d y L y L dx d y L ib b i ia a i n n ,,2,1,,0'0'11 (1-14a) 写成向量形式为

()()βα===??-??b a x y x y y L dx d y L ,0' (1-14b) 其中，y 应有连续的二阶导数，而L 至少应二次连续可微。

对于自由端点情况，边界条件可由下列横截条件确定：

y (x a )＝α，?L/?y ′│xb ＝0=0 终端自由 (1-15)

?L/?y ′│xa ＝0=0 y (x b )＝β 始端自由 (1-16)

由以上对泛函极值曲线y*(x )的讨论可见，求泛函极值问题归结为在给定始端x a 和终端x b 的边界条件下，求解欧拉方程的问题。特别是在两端状态y (x a )和y (x b )自由时，求解欧拉方程所需的边界条件由横截条件来补足。这说明，横截条件是欧拉方程求解时必须满足的端点边界条件。事实上，相对于端点状态y (x a )和y (x b )固定的情况而言，自由端点问题因其一个端点y (x a )或y (x b )，或两个端点自由，因此泛函表达式(1-1)的容许曲线类就被扩展了。所求的极值曲线类y*(x )除满足端点固定的情况外，还必须在自由端点情况下的所有容许曲线中取极值。

对于含有高阶导数的泛函

[]

?=xb xa n dx y y y y x L J )(,,",',; (1-17) 的极值问题，认为函数L 对所有的自变函数而言都是(n +2)阶可微的，并且假定端点上有边界条件

y (x a )＝y a ，y'(x a )＝y'a ，…，y (n -1)(x a )＝y (n -1)a ；

y (x b )＝y b ，y'(x b )＝y'b ，…，y (n -1)(x b )＝y (n -1)b (1-17a)

亦即，在端点上不仅给出函数的值，而且还给出它的直至n -1阶导数的值。再假定，极值在2n 阶可微的曲线y=y (x )上达到。

用与上面相同的求泛函变分δJ 的方法，可以证明

???

?

?????++??+??+??=xb xa n n dx y y L y y L y y L y y L J )()(""''δδδδδ (1-17b) 对式(1-17b)右端第二项分部积分一次，得 ()?????? ????-??=??xb xa xb

xa a b ydx y L dx d x x y y L dx y dx d y L δδδ''' 对式(1-17b)右端第三项分部积分两次，得

()?????

? ????+?????????? ????-???? ????=??xb xa a b xb xa a b ydx y L dx d x x y y L dx d x x y y L dx y dx d y L δδδδ""'""2222

如此连续应用分部积分公式，则最后一项经n 次分部积分后，得 ()?????

? ????-++?????????? ????-???? ????=??--xb xa n n n n a b n n xb xa a b n n n n n ydx y L dx d x x y y L dx d x x y y L dx y dx d y L δδδδ)()2()()1()()()1( 考虑到式(1-17a)的条件，在x ＝x a ，x ＝x b 时，δy ＝δy ′＝δy ′′＝…＝δy (n -1)＝0，

根据变分的极值条件δJ ＝0，则由式(1-17b)得

a b n n a b n n n n a b n n n n xb xa n n n n x x y y L x x y y L dx d y L dx d y L x x y y L dx d y L dx d y L ydx y L dx d y L dx d y L dx d y L J )1()()(22)3()(111)(22

')1('')1(''')1("'------??++??

??????-++???? ????-??+??

??????-++???? ????-??+??

??????-+-???? ????+???? ????-??=?δδδδδ (1-17c) 根据变分基本原理，得

0)1("')(22

=??-+-???

? ????+???? ????-??n n n n y L dx d y L dx d y L dx d y L (1-18) 这个2n 阶微分方程叫作欧拉—泊桑方程，它的积分曲线就是所论变分问题的解——极值曲线。方程所含的2n 个待定常数，由下列2n 个端点条件决定

在x ＝x a ，x ＝x b 处，

y 已知，或 0)1(''')(111=??-++???

? ????-??---n n n n y L dx d y L dx d y L y'已知，或 0)1('')(22)3(=??-++???

? ????-??--n n n n y L dx d y

L dx d y L (1-18a) …

y (n -1)已知，或 0)(=??n y L 以上为仅含有一个自变量的泛函变分问题。对于含有多个自变量的函数的泛函问题，可作如下推导：

研究下列泛函

()[]?????? ?

?????=s dxdy y w x w w y x L y x w J ,,,,, (1-19) 极值问题。函数w (x ，y )区域s 的边界c 上的值已经给出，即在边界c 上w =w 1，w 1为已知值。为把上述问题化为微分方程的边值问题，求J 的一阶变分

dxdy y w w L x w w L w w L dxdy w w L w w L w w L J s y x s y y x x ???????

? ??????+????+??=???? ????+??+??=δδδδδδδ (1-19a) 式中w x =?w/?x ，w y =?w/?y

由于

x x x x w w L w w L x w w L x δδδ??+???? ??????=???? ??????，y y y

y w w L w w L y w w L y δδδ??+???? ??????=???? ??????(a) 因此式(1-19a)可改写为