矩阵的特征值和特征向量

7.3　有限单元法

如上节所述，求解闸门结构动力响应的方法有解耦分析法(振型叠加法)和直接积分法两大类，实际上这两类方法的具体实施均依赖于有限元法。以结构静力的刚度法为基础发展起来的有限元法，现已成为解决工程结构动力问题的一个有力工具。虽然有限元法从命名到现在只不过50多年，但已经相当成熟了。

以节点位移为基本未知量，插值构造单元内部的位移场，从而建立有限元分析的位移协调元，这是有限元形成以来学者们一直研究的主要方面。随着解决问题的日趋大型化、复杂化，采用位移协调元显露出来未知量太多、计算时间过长、费用太高的缺点，同时，位移协调元在分析板壳结构等特殊结构时常常会遇到协调性难以满足的困难，于是，非协调元类型的单元和方法就被提了出来，如混合元、杂交元、加权残数法、样条函数法、半解析半离散法以及耦合法，等等。这些新方法或简洁、或快速、或局部精度更高，为解决结构静力分析和动力分析问题开辟了新的途径，研究者们针对不同的结构予以选择采用，推动了结构动力分析理论的发展。

有限元动力分析遵循一般动力分析的基本原则，在此叙述如下:

(1)自由振动分析——求解闸门结构在特殊初始条件下做自由振动的简谐振动时的自振频率和主振型，最终归结为广义特征值的求解问题或转换为标准特征值的求解问题。

(2)动力响应分析——求解闸门结构由于动力荷载作用下引起强迫振动的振动加速度、动应力和动位移。

有限元法的应用离不开计算机和软件程序，随着有限元法理论的发展和完善，研究者们已开发了许多大型通用有限元程序。如美国开发的国际通用有限元分析程序ANSYS软件等。这些程序容量大、功能强、内容多，大多数工程问题都能被解决，并且还经受了工程实际的检验。

7.3.1　动力分析有限元列式

7.3.1.1　基础知识

动力分析三维有限元——系列方程式的推导，基于三维弹性动力学，在此将三维弹性动力学基本方程予以必要的引述:

以上各式均为张量表达式。式中:

下标“，j”表示求偏导数;

σ_ij、ε_ij为应力张量和应变张量;

u_i、、分别为位移张量、体积力张量和面积力张量;

D_ijkl为弹性系数张量;

n_j为边界外法线的三个方向余弦;

ρ是质量密度;

μ是阻尼系数;

u_i，tt和u_i，t分别是对时间的二次导数和一次导数，即分别表示方向的加速度和速度;

ρu_i，tt和μu_i，t分别代表惯性力和阻尼力(取负值)。

7.3.1.2　连续区域离散化

动力分析问题在引入时间坐标后，事实上处理的是四维(x，y，z，t)问题，不过在有限元分析中都采用部分离散的方法，即对空间域进行离散，这一点与静力分析完全相同。与静力分析类似，动力有限元仍将结构视为仅在节点处连接的有限个单元体的集合体。基本未知量仍为独立的节点位移，当然由于动力响应与时间有关，节点位移是时间的函数。

7.3.1.3　位移模式

设定动力有限元每个单元的动位移分布规律时，所采用的位移模式与静力有限元也相同，即形函数与时间无关，由此构造空间域的插值函数，单元内位移u，v，w的插值表示为:

7.3.1.4　单元集成，建立总体运动方程

集合所有单元的运动方程，建立起结构的总体运动方程。

式(7-8)平衡方程及式(7-12)边界条件等效积分形式的伽辽金提法表示如下:

对上式的第一项∫_Vδu_iσ_ij，j d V进行分部积分，并代入物理方程，可得到

将式(7-14)(此时u₁= u，u₂= v，u₃=w)代入上式，即得到前已提到的结构动力方程式(7-1)，在此重述如下:

7.3.1.5　动力方程与静力方程的比较

(1)动力方程要比静力方程多建立一个质量矩阵和阻尼矩阵。

(2)动力方程是二阶常微分方程，静力方程是线性代数方程。

(3)动力问题要寻求二阶常微分方程组的有效解法，静力问题则要寻求线性代数方程组的有效解法。

7.3.2　自由振动分析

7.3.2.1　特征值问题

1.干扰力为零的无阻尼振动

对于干扰力为零的无阻尼振动，式(7-17)成为

弹性体的自由振动可以分解为一系列的简谐振动的叠加，设这种简谐振动形式的解为

那么二次导数

将(1)、(2)两式代入式(7-18)得

因为{}是非零向量，所以必有

要想获得式(7-18)的解(即(1)式)，只需要确定{}和ω，使之满足(3)式即可。若取λ=ω²，λ称为结构的特征值，(3)式成为前述的广义特征值问题(式7-2)，则

式中:刚度矩阵是对称、正定或半正定的;质量矩阵是对称、正定的;{}则成为对应于特征值λ=ω²的特征向量。

2.式(7-19)的广义特征值问题

由线性代数知，在大多数情况下直接处理它是比较复杂的，一般先将其化为标准特征值问题，然后再进行求解。如式(7-19)就可以化为如下形式的标准特征值问题:

式中:［A］=［M］^－1［K］，虽然质量矩阵和刚度矩阵总是对称矩阵，但经上述的求逆处理后，所得结果未必为对称矩阵。非对称矩阵会带来更多麻烦，所以须事先对［M］矩阵作三角形分解，即将其表述为:

［L］为对角线元素均不为零的下三角阵。在此条件下，式(7-19)才确实化为了形如

的标准特征值问题，上式中有如下关系

7.3.2.2　求取特征值的方法

对于广义特征值问题或标准值问题，求解的要求一般有两种情况:一是求解全部特征值问题，二是求解部分特征值问题。求解的方法有两大类，一类是变换方法，另一类是矩阵迭代方法。

变换方法有雅可比方法、广义雅可比方法、吉文斯-豪斯霍尔德方法、行列式搜索方法等;矩阵迭代方法有幂迭代法、逆迭代法、松弛法、子空间迭代法等。

本文采用子空间迭代法求解结构的部分特征值，下面简述其理论和计算过程。

7.3.3　子空间迭代法的理论及计算方法

子空间迭代法是瑞利—里兹分析法和同时逆迭代法结合的成果。譬如我们若只想求取结构的前p阶特征值(通常结构也只需要前p阶特征值)，那就可以把一个n阶特征值问题化为一个q阶特征值问题(p＜q＜n)，这种方法就是瑞利－里兹分析法。其中基底应用同时逆迭代法，使之接近低阶特征子空间。假如知道n维空间中的一个子空间E_p(p＜n)中的线性无关的向量组{x₁}，{x₂}，…，{x_p}，就可以把求解［K］{u}=λ［M］{u}在子空间E_p的特征值问题化为一个求解［K］{a}=ρ［M］{a}的维数缩小了的特征值部问题，即把n阶特征值问题化为p阶特征值问题，其中ρ为瑞利商。

［K］、［K］分别为［K］、［M］在子空间上的投影。

由此可见，子空间迭代法涉及如下几个问题:

(1)如何寻找线性无关的向量组{x₁}，{x₂}，…，{x_p};

(2)应用瑞利商的极值原理;

(3)如何缩减维数。

下面分别予以讨论。

7.3.3.1　瑞利商的极值原理

将任一向量的瑞利商定义为

对于广义特征值问题

可以根据瑞利商的极值原理把上述问题化为一个等阶的瑞利商的极值问题。

瑞利商的极值原理可表述如下:由式(7-23)可见，瑞利商是一个标量，且是随向量{x}变化的，可以证明，当{x}取为广义特征值问题(式(7-24))的特征向量{φ_i}时，瑞利商ρ({x})达到其极值，这个极值恰恰是和特征向量{φ_i}对应的特征值。

下面证明这一原理。

任一向量{x}可以是特征向量{φ₁}，{φ₂}，…，{φ_n}的任一种线性组合，设

将式(7-25)代入式(7-23)，得

由特征向量正交性，记

将式(7-28)展开，得

而且，当{x}={φ₁}时，我们有a₁= 1，将a₂= a₃=…= a_n= 0代入式(7-29)，有

式(7-31)表明:当{x}取广义特征值问题的一阶特征向量时，ρ({x})取极小值;而当ρ({x})取其i阶特征向量{φ_i}时，ρ({x})取其极值。同时由式(7-31)看出也是最小值也是最大值。

设a₁= a₂=…= a_m－1= 0，则{x}不包含前m－1阶特征向量的分量，即在该基底分量上的投影为0，这就表明{x}和前{m－1}阶特征矢量正交，则由式(7-30)可得

综上可以看出:求解广义特征问题和求瑞利商的极值问题二者是等价的。

7.3.3.2　缩减维数

只要应用基于瑞利商极值原理的瑞利－里兹法就可以把n阶特征值问题降为一个p阶特征值问题。

瑞利－里兹法的实质就是用待定系数和基函数来描述一个任意函数。

现在我们首先把式(7-24)所示特征值问题的前p个特征向量{φ₁}，{φ₂}，…，{φ_p}构成的空间记为E_p空间(p＜n)，这是全部特征向量空间中的一个子空间。

在E_p子空间中任意选择一组线性无关的向量{x₁}，{x₂}，…，{x_p}作为坐标向量，称为子空间E_p的基底。

将任一向量{x}表示为{x₁}，{x₂}，…，{x_p}的线性组合，得

{x}=［X］{a}

式中:［X］=［x₁ x₂…x_p］_n×p;

{a}为p维矢量。

因此，有

式中:［K］=［X］^T［K］［X］

分别为［K］、［M］在E_p子空间上的投影。

式(7-24)所对应广义特征问题的前p阶特征值就是(7-33)所示瑞利商的极值(实质上就是式(7-24)所对应的广义特征值问题在子空间E_p上的特征值)。

而式(7-33)所代表的瑞利商极值总是等价于下述特征值问题:

这样就把式(7-24)所代表的n阶特征值问题化成了式(7-34)所代表的p阶特征值问题。

求解式(7-33)可以求得待定系数{a_i}(i= 1，2，…，p)及特征值ρ_i(i= 1，2，…，p)，此时原问题的前p个特征对为

7.3.3.3初始向量组的选择

在子空间迭代法中要选择一组线性无关的向量组作为E_p子空间的基底(即坐标向量)，这个基底选择得好，则其近似解的精度就高，反之其精度就低。对于复杂结构而言，这个基底很难凭经验选定，这是子空间迭代法的一个缺点。通常先初选一组线性无关的向量组，再利用逆迭代法，使这个初选的基底不断完善，当迭代次数趋于无穷大时，{x^(k+1)}→［φ₁］。也就是说，经过反复迭代后可以使{}→［φ］。_(p)

为了简便，亦可以这样来选择初始向量组，即取{x₁}为［M］矩阵的对角元素，{x₂}，{x₃}，…，{x_p}均为单位向量，且其中“+ 1”放在比值m/k最大的坐标处，即初始向量{x⁽¹⁾}=［M］{x⁽⁰⁾}可以这样构成，其序列号按原问题的对角矩阵diag(m_ii/k_ii)中对角元素的大小排列，而在每一列中把m_ii/k_ii比值最大的那一行置“1”，其余置“0”。

7.3.3.4　计算步骤

首先预选一组线性无关的向量组{x⁽⁰⁾}=［x⁽⁰⁾，x⁽⁰⁾，…，x⁽⁰⁾］，

_12p其方法如前所述，可以直接用［M］{x⁽⁰⁾}表示。如需要求前p阶特征对，通常应使q＞p，以保证所求特征对的精度，但q太多必增加内存容量。在SAP程序中取q=min(2p，p+ 8)。上述{}(i= 1，2，…，q)构成了一个q维子空间，而只是作为q个特征向量{φ_i}(i= 1，2，…，q)构成的子空间E_q的初次近似。通过逆迭代法不断地改善子空间，依次求得，…，，直到→E_q。

1.根据迭代法，作如下迭代:

求出新的向量组{}，形成新的子空间，依次记为，，…，。

2.计算［K］、［M］在子空间上的投影

3.用广义雅可比法求解低阶特征值的问题

由　　　{a}=ρ［］{a}

得到相应的ρ_i、{a_i}等。

4.求出第k+ 1次改进的特征向量矩阵

当k→∞时，［Ω^(k+1)］²→［Q］²

由上述计算步骤可见，在子空间迭代法中，实际上又应用了逆迭代法和广义雅可比法。

7.3.4　动力响应分析

前面在“闸门的动力响应”中已介绍过，求解多自由度体系的强迫振动响应时有两种方法，一是解耦分析法(振型叠加法)，另一是直接积分法。解耦分析法首先采用矩阵迭代法或这里介绍的子空间迭代求解出体系的动态特征参数ω_nj和ξ_i，从而建立主振型矩阵(解耦矩阵)［A_p］;然后用解耦分析法将一个多自由度体系分解为多个单自由度体系的叠加，而单自由度的动力响应问题可以用杜哈美积分求解。这样就可以解决多自由度体系的动力响应问题了。至于用直接积分法求解动力响应问题，就是用计算机求运动微分方程的数值解，亦即用计算机来模拟再现结构随时间变化的过程和现象，所以又称为计算机仿真。

本书着重阐述直接积分法中SAP程序主要采用的方法。

1.欧拉法

直接积分法中最简单的方法是欧拉法。

将运动微分方程式(7-1)写成

可以把整个求解时间区间分为若干个等距的时间步长Δt，根据t时刻的{u(t)}、{u·(t)}及{u··(t)}求出t+Δt时刻的{u(t+Δt)}、{u·(t+Δt)}及{u··(t+Δt)}。

欧拉法的实质是假定求位移时在Δt时间步长内速度为常数，求速度时加速度也为常数，即取位移、速度、加速度为

如果［M］是对角矩阵，则这时完全不需要求解微分方程，这是一种显式解法。

2.线性加速度法

式(7-37)至式(7-39)所示欧拉解法只有当Δt→0时，效果才较好，反之则不稳定。之所以会出现这种问题，是因为实际上在Δt时间内速度、加速度都不一定是常数。

按照泰勒级数，可得

欧拉法实际上只取前两项。为了提高精度，可以从两方面采取措施:

(1)对泰勒级数表达式取更高的项次;

(2)直接用数值积分式(7-40)中的积分值。

从第(2)方面来看，可以说欧拉法事实上是采用了最简单的矩阵形式进行数值积分的。

线性加速度法是综合采用上述两方面的措施来改进欧拉法的一种隐式解法。

这时，位移公式基本上按泰勒展开式的前三项给出:

和泰勒展开式相比较，只是在第三项中用来代替了{u··(t)}，即假定在Δt时间步长内，加速度{u··(t)}是直线变化的。

速度公式中采用梯形公式计算积分{u(t)}d t，即

这里，也是假设加速度为线性变化。

在线性加速度法中，要联立求解式(7-39)、式(7-41)和式(7-42)，可以用迭代法，也可以用直接法。下面介绍一种直接解法(先求{u··(t+Δt)})。

(1)将式(7-41)、式(7-42)代入运动微分方程式，消去{u(t+Δt)}和{u·(t+Δt)}，最后得到

(2)根据式(7-43)求出{^·u^·(t+Δt)}。

其中出现的矩阵［M］+(［C］+)［K］在各阶段求解时(只要Δt为常数)均保持不变，所以仅需作一次三角分解或事先一次求逆。

(3)再根据式(7-41)和式(7-42)求{u(t+Δt)}和{(t+Δt)}。

以上所述直接解法称为第一种直接解法。

也可以先求位移，再求速度和加速度，这种方法称为第二种直接解法。

3.威尔逊－θ法

此法是线性加速度法的一种变形，是该法的进一步扩展。所不同的是，它所取的线性加速度区间不是t～t+Δt，而是t～t+θΔt，通常θ＞1。在这一假定下求得{(t+θΔt)}以后，再作线性内插求得{(t+Δt)}以及{(t+Δt)}和u(t+Δt)}。

其基本关系可以根据线性加速度法的相应关系式获得式(7-37)和式(7-38)，只需要把Δt换成θΔt即可。

然后再和运动微分方程式联立求解即可求得{(t+θΔt)}，再按下式内插求得{(t+Δt)}:

威尔逊－θ法是一种隐式解法。所谓隐式解法是指在求解过程中不断求出新的时刻满足微分方程式的近似解。

隐式解法的稳定性比显式解法好，因为隐式解法中等式右端均有t+Δt时刻的参数(而不像显式解法那样等式右端仅涉及时刻t时的值)，其作用就像自动调节系统中引入了负反馈，可以增强系统的稳定性。当某一阶段计算结果偏差较大时，可以在该步长中得到补偿，从而增强稳定性。

威尔逊－θ法不仅仅是隐式解法，它首先计算时间比时间步长Δt更大的t～θΔt区间内的近似解，又把其后边部分的结果甩掉，而仅取其前面部分(t～t+Δt)的结果作为正式结果，这样就增大了反馈作用，这种巧妙的处理可以提高稳定性，所以威尔逊－θ法是一种非常实用的方法。通常只要取θ＞1.37则不管Δt取多大都是稳定的，即此法是无条件稳定的。但θ也不宜过大，最好取θ≤2，通常取为1.4。

如采用前述第一种直接解法，则其计算步骤如下:

(1)初始计算

①形成闸门结构的刚度矩阵［K］，质量矩阵［M］和阻尼矩阵［C］;

②确定初始值{u₀}，{0}，{0};

③选择时间步长Δt和θ;

④形成等效质量矩阵

(2)依次递推求各时间步长下的终值

①计算t+θΔt时刻的等效荷载项

②求t+θΔt时刻的加速度

③求t+Δt时刻的加速度，速度和位移(按照式(7-46)、式(7-47)、式(7-48)等)。