http://www.paper.edu.cn
杆系结构自由振动的精确有限元法与 1动力刚度法的等价性 袁 驷，叶康生
清华大学土木工程系（100084）
Email: yuans@tsinghua.edu.cn, yeks@tsinghua.edu.cn
摘 要：本文论述了动力刚度法和精确有限元法的等价性，证明了等价性定理、梁氏定理及本文导出的微分等式定理，并为新近提出的导护型Newton法求解振型和频率提出了一个改进的外推公式。文中给出了算例用以表明改进公式的良好效果。
关键词：自由振动，精确有限元法，动力刚度法，杆系结构 1．引言
对于杆系结构的静力分析计算，传统的刚度法（矩阵位移法）和现代的有限元法在某些情况下（如等截面杆件）殊路同归，人们常将二者混为一谈。其实，两种方法的做法是有区别的：
刚度法：求得满足控制微分方程的解，由其解通过取导数直接计算杆端刚度系数，汇
集后得到单元刚度矩阵。
有限元法：构造单元（杆件）形函数，用能量法（虚功原理、最小势能原理）通过积
分生成单元刚度矩阵。
这里，之所以可以得到相同的单元刚度矩阵，其要点是有限元中采用了精确的形函数，即满足控制微分方程的形函数（如等截面梁单元采用三次多项式为形函数）。换言之，两种方法的共同点是都采用了满足控制微分方程的解，刚度法利用其解直接计算刚度系数，而有限元法则利用其解构成形函数并用能量法计算刚度系数。
对于杆系结构的动力分析计算（自由振动问题），也相应地存在两种方法：动力刚度法（Dynamic Stiffness Method, 简称DSM）和精确有限元法（FEM）。但是，多年来，基本上只有动力刚度法得到发展，而精确有限元法则很少为人提及。本文首先介绍动力刚度法，然后讨论精确有限元法的列式，在此基础上进一步讨论两种方法的等价性，由此等价性可导出梁氏定理及另一个等式定理，并为新近提出导护型Newton法求解振型和频率提出了一个改进的外推公式。 本课题由高等学校博士学科点专项科研基金（20020003045）和教育部科技创新工程重大项目培育资金项目（704003）联合资助。 - 1 -1
http://www.paper.edu.cn
为简明起见，以杆件轴向振动作为模型问题来举例说明，其控制微分方程可表示为[1][8]
?d?du(x)?2?EA??ω(x)=0, 0 其中，u(x)为轴向振幅函数，为分布质量，EA为截面抗压刚度，ω为自振频率，l为杆件长度。对于常截面的情况，其控制微分方程可进一步简化为 d2u(ξ)2, ， ， +νu(ξ)=0ν=ωlξ=x0<ξ<12EAdξ 此时，不难得到用结点位移表示的解析解 (1b) u(ξ)=N1(ξ)?1+N2(ξ)?2 (2a) (2b)N1(ξ)=sinν(1?ξ)sinνξ， N2(ξ)= sinνsinν 其中形函数N1(ξ)和N2(ξ)满足控制微分方程(1)，是精确形函数，且是频率ω的函数。记x1=0和x2=l，则若将形函数看作是x的函数N1(x)和N2(x)，有Ni(xj)=δij，这是有限元形函数的一个基本性质。 3．动力刚度法简述 和静力刚度法类似，动力刚度法直接求得满足自由振动控制微分方程的解，由其解直接计算杆端动力刚度系数，汇集后得到单元动力刚度矩阵。以模型问题为例，则杆端动力刚度系数可以用精确形函数直接计算如下： ?dNi??ij=(?1)j?EA??， i,j=1,2 dx??xj 这里，刚度系数上划线表示为局部坐标的量。汇集刚度系数后得单元动力刚度矩阵 e(ω)=(3)νEA lsinν?cosν???1??1?， ?ν=ωlcosν?EA?(4) 这里虽然以模型问题为例，对于横向振动以及其他类型问题，做法大致类似。 将e(ω)按照常规的方法集成可得如下形式的整体动力刚度方程 K*(ω)D=0 (5) 其中，K*是整体动力刚度矩阵；D是振型向量。注意K*(ω)是频率的函数。 - 2 - http://www.paper.edu.cn k [1][2]ωk （以下简称WW计数法）。由WW计数法，在结构的所有频 率中，低于某个给定值ω*的频率的个数J由下式给出： J(ω*)=J0(ω*)+s{K*(ω*)} (6) 其中J0是所有低于ω*的单元固端频率ωF的数目总和，s{K*(ω*)}代表用普通Gauss消元法将K*(ω*)消成上三角阵K*?(ω*)后，主对角线上负元素的个数。WW计数法可以方便地（如配合二分法）得到精确的结构自振频率。 对于振型计算，近几年本文的作者和Williams课题组紧密合作研究之后，提出了一个高效可靠的导护型Newton法[3]，使得多年空缺的振型的精确计算问题得到了很好的解决。导护型Newton法将式(5)关于近似值ωa展成级数，忽略二阶及高阶微量，将问题归结为如下矩阵广义特征值问题 *K* aD=μKa′D 其中 μ=ωa?ωk (7) ***dω在ωa处的取值。用逆幂迭代法得到(μ,D)后，还其中K* a=K(ωa)，且Ka′为dK 可以外推一个精度更高的频率： ωμ=ωa?μ (8) 由此得到的ωμ比ωa更为接近精确解ωk。取ωμ为新的ωa，再用逆幂迭代法解式(7)并得到新的(μ,D)以及新的ωμ，如此迭代直至收敛。 4．精确单元的构造 仍然以模型问题为例，采用式(2b)中的精确形函数，按照常规有限元的做法，单元刚度矩阵和质量矩阵分别为 ij=∫EAx1x2x2dNidNjdx，ij=∫iNjdx x1dxdx(9) 对于等截面单元，由此得 =1νEA 2lsinν?sinνcosν+ν??sym.??(νcosν+sinν)?? sinνcosν+ν?? (10) e=12νsin2ν?ν?sinνcosνsinν?νcosν????? sym.νsinνcosν??? 整体刚度矩阵K和质量矩阵M通过常规的有限元集成手段获得。注意，K和M都是频率ω的函数，即K=K(ω)和M=M(ω)。精确有限元法将问题归结为如下矩阵广义特征 - 3 - http://www.paper.edu.cn KaD=ω2MaD (11) 5．两种方法的等价性 本节将证明如下三个关系 &*=?M， K&=ω2M& K*=K?ω2M， K(12a,b,c)
杆系结构自由振动的精确有限元法与动力刚度法的等价性.
