模态叠加法的基本原理及其实现方法(Matlab)

  • A+
所属分类:学霸天地

一、系统动力学方程的解耦

坐标变换

将物理坐标x通过坐标变换x=Φq转换到模态空间对系统原振动微分方程进行解耦,转换为模态坐标q下的(用模态坐标表示的)系统振动微分方程。

对于系统的振动微分方程:

模态叠加法的基本原理及其实现方法(Matlab)

这里质量矩阵和刚度矩阵均为实对称矩阵,阻尼矩阵也为实对称矩阵。

由于此常微分方程组中各个微分方程相互耦合,无法单独求解。因此,要想办法通过坐标变换,将此微分方程组解耦——即各个微分方程相互独立,可以单独求解。

x=Φq入上式有:

模态叠加法的基本原理及其实现方法(Matlab)

模态叠加法的基本原理及其实现方法(Matlab)

当方程可以解耦时,应有:

模态叠加法的基本原理及其实现方法(Matlab)

式中MpCpKp均为如下形式的对角阵:

模态叠加法的基本原理及其实现方法(Matlab)

即:

模态叠加法的基本原理及其实现方法(Matlab)

其中:

模态叠加法的基本原理及其实现方法(Matlab)

公式(8)为n阶微分方程组,各个坐标qi(i=1,2,...,n)之间相互独立,可以单独求解,即实现了原方程组的解耦。

问题是,如何求得此Φ?

二、坐标变换矩阵Φ的求解
求解广义特征值问题

将(4)式变为:

模态叠加法的基本原理及其实现方法(Matlab)

入(6)式得:

模态叠加法的基本原理及其实现方法(Matlab)

其中MpKp均为对角阵,从而有

模态叠加法的基本原理及其实现方法(Matlab)

其中对角线上元素为:

模态叠加法的基本原理及其实现方法(Matlab)

于是式(11)变为:

模态叠加法的基本原理及其实现方法(Matlab)

式(14)为矩阵理论中的广义特征值问题。

如果将矩阵Φ按列分开:

模态叠加法的基本原理及其实现方法(Matlab)

其中:

模态叠加法的基本原理及其实现方法(Matlab)

由分块矩阵乘法:

模态叠加法的基本原理及其实现方法(Matlab)

模态叠加法的基本原理及其实现方法(Matlab)

则式(14)就变为如下熟悉的特征值形式:

模态叠加法的基本原理及其实现方法(Matlab)

之所以叫“广义”是因为M的存在。如果矩阵M为单位阵,那么就是常规的特征值问题。

于是,原方程(1)的解耦问题最终归结为求解特征值问题(14),相应的特征值方程为:

模态叠加法的基本原理及其实现方法(Matlab)

一旦求出这个多项式的根λ回式(14),便可以确定变换矩阵Φ

三、归一化问题
关于主质量归一化

习惯上我们让方程的第一项系数为1,即:

模态叠加法的基本原理及其实现方法(Matlab)

问题是,如何找到这样的Φ,使得Mp=E,我们令此时的ΦΦN ,则有:

模态叠加法的基本原理及其实现方法(Matlab)

问题是ΦN如何求?

首先引入振型的正交性

方法一:

由前面的论述

模态叠加法的基本原理及其实现方法(Matlab)

可见,若λiλj,则有:

模态叠加法的基本原理及其实现方法(Matlab)

因为:

模态叠加法的基本原理及其实现方法(Matlab)

可得:

模态叠加法的基本原理及其实现方法(Matlab)

因此有:

模态叠加法的基本原理及其实现方法(Matlab)

由这两个式子可见:任意2个振型之间,既有对M的正交性,又有对K的正交性,它们统称为振型的正交性。

方法二:

设振系的第i个与第j个振型向量分别为ΦiΦj,按照振型方程(19)有:

模态叠加法的基本原理及其实现方法(Matlab)

可得:

模态叠加法的基本原理及其实现方法(Matlab)

模态叠加法的基本原理及其实现方法(Matlab)

因为MK都是对称矩阵,故将上式转置后得:

模态叠加法的基本原理及其实现方法(Matlab)

所以有:

模态叠加法的基本原理及其实现方法(Matlab)

λiλj则有正交关系:

模态叠加法的基本原理及其实现方法(Matlab)

同样有:

模态叠加法的基本原理及其实现方法(Matlab)

由这两个式子可见:任意2个振型之间,既有对M的正交性,又有对K的正交性,它们统称为振型的正交性。因此,由振型的正交性有:

模态叠加法的基本原理及其实现方法(Matlab)

模态叠加法的基本原理及其实现方法(Matlab)

故对于ΦN有:

模态叠加法的基本原理及其实现方法(Matlab)

因此,必须有:

模态叠加法的基本原理及其实现方法(Matlab)

令:

模态叠加法的基本原理及其实现方法(Matlab)

入上式有:

模态叠加法的基本原理及其实现方法(Matlab)

因此:

模态叠加法的基本原理及其实现方法(Matlab)

所以:

模态叠加法的基本原理及其实现方法(Matlab)

此时,由式(14)式得:

因此,对原系统振动微分方程做如下变换:

(1) 令x=ΦNqN代入上式有:

模态叠加法的基本原理及其实现方法(Matlab)

(2) 可得

模态叠加法的基本原理及其实现方法(Matlab)

即:

模态叠加法的基本原理及其实现方法(Matlab)

上式中qN称为正则坐标,fNp为正则激励力。

四、解耦方程的求解

对于动力学方程

模态叠加法的基本原理及其实现方法(Matlab)

当阻尼矩阵Cp可对角化时,其分量形式为:

模态叠加法的基本原理及其实现方法(Matlab)

两边同除以mi有:

模态叠加法的基本原理及其实现方法(Matlab)

其中:

模态叠加法的基本原理及其实现方法(Matlab)

则有:

模态叠加法的基本原理及其实现方法(Matlab)

(之所以这样做是仿效一元二次方程的形式)上式化为:

模态叠加法的基本原理及其实现方法(Matlab)

当采用正则振型解耦时,

模态叠加法的基本原理及其实现方法(Matlab)

模态叠加法的基本原理及其实现方法(Matlab)

可化为:

模态叠加法的基本原理及其实现方法(Matlab)

(注意:由于已经关于质量正则化,故这里的mi=1)

上述方程的解为:

模态叠加法的基本原理及其实现方法(Matlab)

其中:

模态叠加法的基本原理及其实现方法(Matlab)

利用x=ΦNqN可求得物理坐标下的响应。

注意,我们已知的初始条件是物理坐标下的,而解耦方程是模态坐标下的方程,因此,需要将物理坐标下的初始条件转换到模态坐标下去。

如何将初始条件变换到模态坐标(正则坐标)下?

我们知道x=ΦNqN,因此似乎可以直接进行变换。但实际上,由于进行矩阵的求逆运算计算量非常大,而且数值计算中会产生误差,因此,我们不采用此方法进行变换。

将其两端同时左乘则有:

模态叠加法的基本原理及其实现方法(Matlab)

由于:

模态叠加法的基本原理及其实现方法(Matlab)

因此有:

模态叠加法的基本原理及其实现方法(Matlab)

五、模态叠加法步骤

1、求解广义特征值问题,求得特征值λr和特征向量Φr

广义特征值问题为=MΦΛ,其特征方程为:

模态叠加法的基本原理及其实现方法(Matlab)

在MATLAB中利用如下语句求得特征值和特征向量:

[V,D]=eig(K,M);

%特征向量矩阵和特征值矩阵

其中D为特征值矩阵(对角阵),V为对应的特征向量矩阵。

2、对特征值和特征向量进行排序

求得固有频率ωi和振型Φi后,进而得到振型矩阵Φ(习惯上令振型Φi中第一个元素为1)。

我们习惯上对特征值按从小到大的顺序进行排序,这就要求对得到的特征值进行排序,其对应的特征向量也要进行相应的排序。关于排序的算法,有很多,这里采用冒泡法进行排序。

由于:

模态叠加法的基本原理及其实现方法(Matlab)

可求得圆频率,进而可得到固有频率

模态叠加法的基本原理及其实现方法(Matlab)

V1=zeros(n,n);

VN=zeros(n,n);

for i=1:n

V1(:,i)=V(:,i)/V(1,i);

%求振型矩阵,各列除以各列的第一个元素

end

3、振型矩阵Φ关于主质量归一化,得到正则振型矩阵ΦN

由前面论述,可利用:

模态叠加法的基本原理及其实现方法(Matlab)

求得正则振型矩阵ΦN

Mp=V1'*M*V1;

%主质量矩阵

Kp=V1'*K*V1;

%主刚度矩阵

Cp=V1'*C*V1;

for i=1:n

VN(:,i)=V1(:,i)/sqrt(Mp(i,i));

%求正则振型矩阵

end

验证程序:

MNp=VN'*M*VN;

%结果为1 对角阵,与理论相符合

KNp=VN'*K*VN;

%结果为特征值对角阵,与理论相符合

CNp=VN'*C*VN;

%若CNp为对角阵,则原方程组可以解耦

4、将初始条件变换到正则坐标上。

模态叠加法的基本原理及其实现方法(Matlab)

模态叠加法的基本原理及其实现方法(Matlab)

MATLAB中如下实现:

x=zeros(n,1);

%物理坐标系下位移向量

xd=zeros(n,1);

%物理坐标系下速度向量

x0=[0 0]';

%物理坐标系下初位移向量

xd0=[0 0]';

%物理坐标系下初速度向量

qN=zeros(n,1);

%模态坐标系下位移向量

qNd=zeros(n,1);

%模态坐标系下速度向量

%模态坐标系下初位移

qN0=VN'*M*x0;

%物理坐标系下初位移转换到模态坐标系下初位移

%模态坐标系下初速度

qNd0=VN'*M*xd0;

%物理坐标系下初速度转换到模态坐标系下初速度

5、求振系在正则坐标系下的响应。

模态叠加法的基本原理及其实现方法(Matlab)

MATLAB中代码如下:

deltat=0.01;

%正则坐标系下的响应公式中卷积积分t的分段

t0=0;tf=1;

%动力学响应的起止时间

t=t0:deltat:tf;

qNt=zeros(n,length(t));

%用于存储各时刻模态坐标系下的位移

m=length(t);

%%%%%%%%%%%%%%%

f=[-50000*1.6*sin(57.5*t)-10*1.6*57.5*cos(57.5*t);

50000*1.6*sin(57.5*t)+10*1.6*57.5*cos(57.5*t)];

fN=VN'*f;

%fN 为正则激励力

Q=zeros(n,2*m-1);

%动力学响应中的卷积积分部分

%%%%%%%%%%%%%%%

for i=1:n

h=1/Wd(i)*sin(Wd(i)*t).*exp(-ksi(i)*Wn(i)*t);

Q(i,:)=deltat*conv(fN(i,:),h);

qNt(i,:)=exp(-ksi(i)*Wn(i)*t).*(qN0(i)*cos(Wd(i)*t)+(qNd0(i)+ksi(i)*Wn(i)*qN0(i))/Wd(i)*sin(Wd(i)*t)+Q(i,1:1/deltat+1));

end

6、将正则坐标系下的响应再变回到物理坐标下。

模态叠加法的基本原理及其实现方法(Matlab)

MATLAB中如下实现:

xt=VN*qNt;

%求解物理坐标x

六、模态叠加法的修正

前面讨论的模态叠加法是基于系统的特征值没有重根。

利用振型解耦的前提是已经有了全部的振型向量,如果所有特征值互异,那么肯定存在对应的特征向量,它们对刚度矩阵和质量矩阵正交。

重频特征向量的非唯一性:如果特征方程有两个根λ1=λ2,那么数学上可以证明确实存在Φ1Φ2两个特征向量(至少就对称的刚度矩阵和质量矩阵,这是成立的)。但是正交性却无法自动保证。(陈奎孚机械振动基础P214)

可以证明:对应于重根λ有无穷多振型向量。

我记得线性代数中指出什么时候即使有重根也会有n个相互独立的特征向量的。方法是施密特正交化。

如果已经知道了重频所对应的两个独立振型向量,则总可以通过线性代数中的格拉姆——施密特正交化过程找到两个正交向量。例如有两个独立振型向量Φ1Φ2,那么先保留Φ1,而第二个向量假定为:

模态叠加法的基本原理及其实现方法(Matlab)

它肯定也是对应于同一频率的特征向量,现在目标是要调整μ使得上诉变量与已经选定的Φ1也加权正交,也就是:

模态叠加法的基本原理及其实现方法(Matlab)

即:

模态叠加法的基本原理及其实现方法(Matlab)

因为质量矩阵总是正定的,所以上式第二项非系数部分定非0,故可以解出:

模态叠加法的基本原理及其实现方法(Matlab)

这样就保证了加权的正交性。

三个重根的证明:

模态叠加法的基本原理及其实现方法(Matlab)

则应有:

模态叠加法的基本原理及其实现方法(Matlab)

即:

模态叠加法的基本原理及其实现方法(Matlab)

由于

模态叠加法的基本原理及其实现方法(Matlab)

故有:

模态叠加法的基本原理及其实现方法(Matlab)

模态叠加法的基本原理及其实现方法(Matlab)

始发于微信公众号:声振之家

发表评论

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: