动力学方程数值解法：振型叠加法与Duhamel积分

对于一个实际结构，由有限元法离散化处理后，动力学方程可写为：

振型叠加法

按照有限单元法的一般规则，经过边界条件的约束处理，结构在强迫振动时多自由度体系的运动平衡方程可以表示为：

其中，M是体系的质量矩阵，C是体系的阻尼矩阵，而K则是刚度矩阵，R为外荷载向量。

上述动力平衡方程实质上是与加速度有关的惯性力和与速度有关的阻尼力及与位移有关的弹性力在时刻t与荷载的静力平衡。

振型叠加法是把多自由度体系的结构的整体振动分解为与振型次数相对应的单自由度体系，求得各个单自由度体系的动力响应后，再进行叠加得出结构整体响应。

振型叠加法原理是利用结构无阻尼自由振动的振型矩阵作为变换矩阵，将结构动力方程式(1)式变换成一组非耦合的微分方程。逐个地求解这些方程后，将解叠加即可得到动力方程的解。

将体系单元节点的位移向量表示为如下的变换形式：

式中的变换矩阵Φ是由动力方程对应的无阻尼自由振动方程解出的前m阶振型矩阵。即Φ=[φ1，φ2，...φm]，x(t)是与时间有关的m阶向量，x的各分量称为广义位移。

将式(2)代入动力方程(1)并左乘以Φt ，则可得广义位移为未知数的方程：

式中

现在进一步考察式(4)，考虑到特征向量的正交性，可得：

于是对应于振型的广义位移的平衡方程(3)可改写为：

其中，Λ为特征值

将式(2)稍加运算可得广义位移用有限元位移表示的形式：

在(6)式中，当忽略了阻尼的影响，平衡方程为互不耦合的，可以对每个方程逐个地进行时间积分。出于相同的考虑，在对有阻尼的体系进行分析时仍然希望采用相同的计算过程去求解互不耦合的平衡方程式。问题是式(6)中的阻尼阵C通常不能象体系的质量阵和刚度阵那样由单元的刚度阵和质量阵装配而成。但当假定阻尼与固有频率成比例吗，即假定：

式中，ξ1是振型阻尼参数；δij是Kronecker符号(当i=j时，δij=1。当i≠j时，δij=0)。

这时式(6)可简化为如下形式的若干个方程式：

其中xi(t)的初始条件为下式：

式(10)表示了一个具有单位质量，刚度为ω12的自由度体系当阻尼比为ξi时的运动平衡控制方程。这个平衡方程的求解可通过计算Duhamel积分求得。

式中

当利用式(9)来考虑阻尼的影响时意味着假设结构的总阻尼是每个振型的阻尼之和，而每个振型上的阻尼是能够量测的，况且在大多数情况下结构的阻尼比更易于量测。因而便于用来近似地反映结构体系的阻尼特性。 

同时在计算上也避免计算阻尼阵而只需计算刚度阵和质量阵。

Duhamel积分递推公式

对以上方程式(10)，考虑某一模态的振动，并略去下标i可写为：

考虑初始条件为：

此时其定解为：

式中：

将上式对时间求导，并将t，t+2△(其中2△为时间步长)带入t0，最后通过抛物线法则计算式中的卷积则可得到：

记

则上述A矩阵元素为：

应用上述递推公式，以前一时刻来求后一时刻的结果。计算不重复。当aij求出后，以后在时域中的步进求解只是一些简单的数组相乘。计算速度很快。

Duhamel积分matlab实例

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

%单自由度Duhamel积分

%内部采用Simpson积分

%初始状态为静止状态

%可计算无阻尼和阻尼强迫震动

%输入可为数组或函数荷载

%只能得出位移时程图

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

clear all;

%w=30

w=input('输入自振频率 ω:');

%n=10

n=input('输入荷载步数n ：');

%m=96.6/32.3

m=input('输入单元质量m ：');

%T=0.05

T=input('输入外荷载持续时间T：');

%si=0.0

si=input('输入阻尼系数ξ：');

%k=2700

k=input('输入刚度系数k：');

deta=T/n;

wD=w*(1-si^2)^0.5;

G=deta/(m*wD)/3;

t1=[0:deta:T];

reply = input('输入荷载数组直接回车或敲击Y，输入函数荷载输入N： [Y]:','s');

if isempty(reply)

reply = 'Y';

end

if reply=='Y'

p2=input('输入荷载数组并用[]抱住为n＋1个元素：例如[019.3200 38.6400 57.9600 77.2800 96.6000 77.2800 57.9600 38.6400 19.32000]:\n');

 

%p2=[0:19.32:96.6,77.28:-19.32:0]

exp1=exp(-2*si*w*deta);

exp2=4*exp(-si*w*deta);

A(1)=p2(1)*cos(wD*0);

for i=2:n/2+1

A(i)=(A(i-1)+p2(2*i-3)*cos(wD*(2*deta*(i-2))))*exp1+p2(2*i-2)*cos(wD*(2*i*deta-3*deta))*exp2+p2(2*i-1)*cos(wD*2*(i-1)*deta);

end

 

B(1)=p2(1)*sin(wD*0);

for i=2:n/2+1

B(i)=(B(i-1)+p2(2*i-3)*sin(wD*(2*deta*(i-2))))*exp1+p2(2*i-2)*sin(wD*(2*i*deta-3*deta))*exp2+p2(2*i-1)*sin(wD*2*(i-1)*deta);

end

 

t1=[0:2*deta:T];

disp('位移随时间的变化');

V=G*(A.*sin(wD*t1)-B.*cos(wD*t1));

disp('弹性力随时间的变化');

f=k*V;

for i=1:n/2+1

p22(i)=p2(2*i-1);

end

disp('加速度随时间的变化');

v11=(p22-k*V)/m;

subplot(1,2,1),plot(t1,V);

title('位移'),grid on,

 

 

t2=[T:deta:T*5];

sw=sin(w*t2);

cw=cos(w*t2);

disp('扩大四倍后位移随时间的变化');

y=G*A(n/2+1)*sw-G*B(n/2+1)*cw;

subplot(1,2,2),

 

plot(t1,V,'r') ,hold on;

plot(t2,y,'r');

title('加载时间增加四倍时的位移位移时程曲线');

grid on;

else

 

syms tao t;

y=input('函数荷载(关于未知数tao)：');

%y=-10000/8*tao+100;

v1=1/m/wD*int(y*exp(-si*w*(t-tao))*sin(wD*(t-tao)),tao,0,t);

v01=diff(v1,t);

v02=diff(v01,t)

t1=[0:T/n:T];

for i=1:n+1

v2(i)=subs(v1,t,t1(i));

v011(i)=subs(v01,t,t1(i));

v022(i)=subs(v02,t,t1(i));

end

disp('荷载作用时间内位移变化');

v2

disp('荷载作用时间内位移最大值');

max(v2)

disp('荷载作用时间内速度变化');

v011

disp('荷载作用时间内速度最大值');

max(v011)

disp('荷载作用时间内加速度变化');

v022

disp('荷载作用时间内加速度最大值');

max(v022)

subplot(2,3,1);

plot(t1,v2),grid on;

title('荷载作用时间内位移时程曲线');

subplot(2,3,2);

plot(t1,v011),grid on;

title('荷载作用时间内速度时程曲线');

subplot(2,3,3);

plot(t1,v022),grid on;

title('荷载作用时间内加速度时程曲线');

A=1/m/wD*int(y*exp(si*w*tao-si*w*t)*cos(wD*tao),tao,0,T);

B=1/m/wD*int(y*exp(si*w*tao-si*w*t)*sin(wD*tao),tao,0,T);

 

A1=subs(A,T);

B1=subs(B,T);

V33=A1*sin(wD*t)-B1*cos(wD*t);

V033=diff(V33,t);

V0033=diff(V033,t);

 

t2=[T:T/n:5*T];

v33=A1*sin(wD*t2)-B1*cos(wD*t2);

for i=1:4*n+1

V03(i)=subs(V033,t,t2(i));

V003(i)=subs(V0033,t,t2(i));

end

disp('扩大四倍后位移随时间的变化');

v33

disp('扩大四倍后位移最大值');

max(v33)

disp('扩大四倍后速度随时间的变化');

V03

disp('扩大四倍后速度最大值');

max(V03)

disp('扩大四倍后加速度随时间的变化');

V003

disp('扩大四倍后加速度最大值');

max(V003)

subplot(2,3,4);

plot(t1,v2),hold on;

title('加载时间增加四倍时的位移时程曲线');

plot(t2,v33),grid on;

subplot(2,3,5);

plot(t1,v011),hold on;

title('加载时间增加四倍时的速度时程曲线');

plot(t2,V03),grid on;

subplot(2,3,6);

plot(t1,v022),hold on;

title('加载时间增加四倍时的加速度时程曲线');

plot(t2,V003),grid on;

end