對于一個實(shí)際結(jié)構(gòu)，由有限元法離散化處理后，動力學(xué)方程可寫為：
按照有限單元法的一般規(guī)則，經(jīng)過邊界條件的約束處理，結(jié)構(gòu)在強(qiáng)迫振動時多自由度體系的運(yùn)動平衡方程可以表示為：
其中，M是體系的質(zhì)量矩陣，C是體系的阻尼矩陣，而K則是剛度矩陣，R為外荷載向量。
上述動力平衡方程實(shí)質(zhì)上是與加速度有關(guān)的慣性力和與速度有關(guān)的阻尼力及與位移有關(guān)的彈性力在時刻t與荷載的靜力平衡。
振型疊加法是把多自由度體系的結(jié)構(gòu)的整體振動分解為與振型次數(shù)相對應(yīng)的單自由度體系，求得各個單自由度體系的動力響應(yīng)后，再進(jìn)行疊加得出結(jié)構(gòu)整體響應(yīng)。
振型疊加法原理是利用結(jié)構(gòu)無阻尼自由振動的振型矩陣作為變換矩陣，將結(jié)構(gòu)動力方程式(1)式變換成一組非耦合的微分方程。逐個地求解這些方程后，將解疊加即可得到動力方程的解。
將體系單元節(jié)點(diǎn)的位移向量表示為如下的變換形式：
式中的變換矩陣Φ是由動力方程對應(yīng)的無阻尼自由振動方程解出的前m階振型矩陣。即Φ=[φ1，φ2，...φm]，x(t)是與時間有關(guān)的m階向量，x的各分量稱為廣義位移。
將式(2)代入動力方程(1)并左乘以Φt ，則可得廣義位移為未知數(shù)的方程：
式中
現(xiàn)在進(jìn)一步考察式(4)，考慮到特征向量的正交性，可得：
于是對應(yīng)于振型的廣義位移的平衡方程(3)可改寫為：
其中，Λ為特征值
將式(2)稍加運(yùn)算可得廣義位移用有限元位移表示的形式：
在(6)式中，當(dāng)忽略了阻尼的影響，平衡方程為互不耦合的，可以對每個方程逐個地進(jìn)行時間積分。出于相同的考慮，在對有阻尼的體系進(jìn)行分析時仍然希望采用相同的計(jì)算過程去求解互不耦合的平衡方程式。問題是式(6)中的阻尼陣C通常不能象體系的質(zhì)量陣和剛度陣那樣由單元的剛度陣和質(zhì)量陣裝配而成。但當(dāng)假定阻尼與固有頻率成比例嗎，即假定：
式中，ξ1是振型阻尼參數(shù)；δij是Kronecker符號(當(dāng)i=j時，δij=1。當(dāng)i≠j時，δij=0)。
這時式(6)可簡化為如下形式的若干個方程式：
其中xi(t)的初始條件為下式：
式(10)表示了一個具有單位質(zhì)量，剛度為ω12的自由度體系當(dāng)阻尼比為ξi時的運(yùn)動平衡控制方程。這個平衡方程的求解可通過計(jì)算Duhamel積分求得。
式中
當(dāng)利用式(9)來考慮阻尼的影響時意味著假設(shè)結(jié)構(gòu)的總阻尼是每個振型的阻尼之和，而每個振型上的阻尼是能夠量測的，況且在大多數(shù)情況下結(jié)構(gòu)的阻尼比更易于量測。因而便于用來近似地反映結(jié)構(gòu)體系的阻尼特性。
同時在計(jì)算上也避免計(jì)算阻尼陣而只需計(jì)算剛度陣和質(zhì)量陣。
對以上方程式(10)，考慮某一模態(tài)的振動，并略去下標(biāo)i可寫為：
考慮初始條件為：
此時其定解為：
式中：
將上式對時間求導(dǎo)，并將t，t+2△(其中2△為時間步長)帶入t0，最后通過拋物線法則計(jì)算式中的卷積則可得到：
記
則上述A矩陣元素為：
應(yīng)用上述遞推公式，以前一時刻來求后一時刻的結(jié)果。計(jì)算不重復(fù)。當(dāng)aij求出后，以后在時域中的步進(jìn)求解只是一些簡單的數(shù)組相乘。計(jì)算速度很快。
%%%%%%%%%%%%%%%%%
%單自由度Duhamel積分
%內(nèi)部采用Simpson積分
%初始狀態(tài)為靜止?fàn)顟B(tài)
%可計(jì)算無阻尼和阻尼強(qiáng)迫震動
%輸入可為數(shù)組或函數(shù)荷載
%只能得出位移時程圖
%%%%%%%%%%%%%%%%%
clear all;
%w=30
w=input('輸入自振頻率 ω:');
%n=10
n=input('輸入荷載步數(shù)n ：');
%m=96.6/32.3
m=input('輸入單元質(zhì)量m ：');
%T=0.05
T=input('輸入外荷載持續(xù)時間T：');
%si=0.0
si=input('輸入阻尼系數(shù)ξ：');
%k=2700
k=input('輸入剛度系數(shù)k：');
deta=T/n;
wD=w*(1-si^2)^0.5;
G=deta/(m*wD)/3;
t1=[0:deta:T];
reply = input('輸入荷載數(shù)組直接回車或敲擊Y，輸入函數(shù)荷載輸入N： [Y]:','s');
if isempty(reply)
reply = 'Y';
end
if reply=='Y'
p2=input('輸入荷載數(shù)組并用[]抱住為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('擴(kuò)大四倍后位移隨時間的變化');
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('函數(shù)荷載(關(guān)于未知數(shù)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('荷載作用時間內(nèi)位移變化');
v2
disp('荷載作用時間內(nèi)位移最大值');
max(v2)
disp('荷載作用時間內(nèi)速度變化');
v011
disp('荷載作用時間內(nèi)速度最大值');
max(v011)
disp('荷載作用時間內(nèi)加速度變化');
v022
disp('荷載作用時間內(nèi)加速度最大值');
max(v022)
subplot(2,3,1);
plot(t1,v2),grid on;
title('荷載作用時間內(nèi)位移時程曲線');
subplot(2,3,2);
plot(t1,v011),grid on;
title('荷載作用時間內(nèi)速度時程曲線');
subplot(2,3,3);
plot(t1,v022),grid on;
title('荷載作用時間內(nèi)加速度時程曲線');
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('擴(kuò)大四倍后位移隨時間的變化');
v33
disp('擴(kuò)大四倍后位移最大值');
max(v33)
disp('擴(kuò)大四倍后速度隨時間的變化');
V03
disp('擴(kuò)大四倍后速度最大值');
max(V03)
disp('擴(kuò)大四倍后加速度隨時間的變化');
V003
disp('擴(kuò)大四倍后加速度最大值');
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
