ganggouren
Duhamel积分、Fourier变换法、中心差分方法、Newmark法、功率谱法……傻傻分不清楚，相信有不少小伙伴像笔者一样好奇：这些方法精度如何，结果有什么特点？以下笔者将以单自由度弹性体系为例，对比这几种方法的精度。为了便于感兴趣的小伙伴重现，同时附上了Matlab代码。
算例为《结构动力学》（刘晶波主编）习题5.2
Duhamel积分是一种时域分析方法，它将荷载分解为一系列脉冲，获得每一个脉冲作用下结构的反应，然后叠加每一脉冲作用下的反应得到结构总的反应。Duhamel积分方法以积分的方式给出了体系运动的解析表达式（忽略离散采样带来的误差的情况下），但从实际应用上看，当采用数值积分时，其计算效率不高，因为对于计算任一个时间点t的反应，积分都要从0积到t，而实际要计算一时间点系列，可能要几百到几千个点，计算量很大。因为使用了叠加原理，仅适用于线弹性分析。
Matlab程序如下：
function y=duhamel(dt)
%duhamel
m=17.5; %质量
k=875.5; %刚度
zeta=0.14138;%阻尼比
%dt=0.01; %时间步长
t0=0; %起始时间
t2=6.4; %结束时间
w0=sqrt(k/m);
w1=w0*sqrt(1-zeta^2);
t=t0:dt:t2;
y=t;
for i=1:(length(t))
x=linspace(t(1),t(i));
px=(100*x).*(x>=0&x<=0.4)+(80-100*x).*(x>0.4&x<0.8);
a=px.*exp(zeta*w0*x).*cos(w1*x);
A=trapz(x,a);
b=px.*exp(zeta*w0*x).*sin(w1*x);
B=trapz(x,b);
y(i)=exp(-zeta*w0*t(i))*(A*sin(w1*t(i))-B*cos(w1*t(i)))/(m*w1);
end
ymax=max(y)
figure
plot(t,y);
Fourier变换法是一种频域分析方法，其基本计算步骤是：
1）对外荷载做Fourier变换，得到外荷载的Fourier谱；
2）利用复频反应函数得到反应的频域解；
3）应用Fourier逆变换，得到反应的时域解。
在用频域法分析中涉及到两次Fourier变换，均为无穷域积分，特别是Fourier逆变换，被积函数是复数，有时涉及复杂的围道积分。当外荷载是复杂的时间函数（如地震动）时，用解析型的Fourier变换几乎是不可能的，实际计算中大量采用的是离散Fourier变换。因为使用了叠加原理，仅适用于线弹性分析。
Matlab程序如下：
function [un,tf]=qiaoFFT(dt,N)
% Fourier
%% 执行FFT点数为64
% 构建原信号
% dt=0.1;
% N=64;
t=[0:N-1]*dt; % 时间序列
xn=(100*t).*(t>=0&t<=0.4)+(80-100*t).*(t>0.4&t<0.8);
subplot(2,2,1)
plot(t,xn) % 绘出原始信号
xlabel(‘时间/s’),title(‘原始信号’)
axis([0 6.4 0 50]) % 调整坐标范围
% FFT分析
NN=N; % 执行64点FFT
XN=fft(xn,NN); %
f0=1/(dt*NN); % 基频
f=[0:NN-1]*f0; % 频率序列
A=real(XN); % 幅值序列
subplot(2,2,2);
stem(f,A),xlabel(‘频率/Hz’) % 绘制频谱
axis([0 10 -200 200]) % 调整坐标范围
title(‘荷载傅里叶谱’);
%% H(iw)
k=875.5;
zeta=0.14138;
fn=1.126288;
HN=ones(1,NN)./(1-(f/fn).*(f/fn)+2i*zeta*(f/fn))/k;
for i=1:NN/2
HN(NN+1-i)=conj(HN(i));
end
UN=HN.*XN;
subplot(2,2,3);
stem(f,abs(HN)),xlabel(‘频率/Hz’) % 绘制频谱
axis([0 10 -0.01 0.01]) % 调整坐标范围
title(‘复频反应函数’);
un=ifft(UN,NN);
subplot(2,2,4)
plot(t,un) %
xlabel(‘时间/s’),title(‘ut’)
axis([0 6.4 -0.1 0.1]) % 调整坐标范围
tf=t
ymax=abs(max(un))
荷载Fourier谱、复频反应函数及时域位移反应如下图所示。
中心差分方法用有限差分代替位移对时间的求导（即速度和加速度），在计算i+1时刻的运动时，需要已知i和i-1两个时刻的运动，属于两步法，它具有2阶精度，是有条件稳定的，稳定条件为, 是显式积分方法，不需要对刚度矩阵求逆，具有较高的计算效率。
Matlab程序如下：
function u=central(dt)
% central
m=17.5; %质量
k=875.5; %刚度
c=35;%阻尼比
%dt=0.1; %时间步长
t0=0; %起始时间
t2=6.4; %结束时间
t=t0:dt:t2;
u=t;u(1)=0;u(2)=0;
k1=m/dt/dt+c/2/dt;
b=m/dt/dt-c/2/dt;
c=2*m/dt/dt;
for i=2:(length(t)-1)
x=t(i);
pi=(100*x)*(x>=0&x<=0.4)+(80-100*x)*(x>0.4&x<0.8);
fs=k*u(i);
pi1=pi-fs+c*u(i)-b*u(i-1)
u(i+1)=pi1/k1;
end
ymax=max(u)
Newmark方法同样将时间离散化，运动方程仅要求在离散的时间点上满足。与中心差分法不同的是，它不是用差分对i时刻的运动方程展开，得到外推计算i+1时刻位移的公式，而是通过对加速度的假设，以i时刻的运动量为初始值，通过积分得到计算i+1时刻的运动公式，计算过程中需要对刚度矩阵求逆，是隐式方法。当γ=1/2, β=1/4时，是无条件稳定的，就是常加速度法；当γ=1/2, β=1/6时，就是线性加速度法。
Matlab程序如下：
function u=newmark(dt,beta)
% newmark
m=17.5; %质量
k=875.5; %刚度
c=35;%阻尼比
%dt=0.1; %时间步长
t0=0; %起始时间
t2=6.4; %结束时间
t=t0:dt:t2;
u=t;u(1)=0;
v=t;v(1)=0;
a=t;a(1)=0;
gama=0.5;
a0=1/beta/dt/dt;
a1=gama/beta/dt;
a2=1/beta/dt;
a3=1/2/beta-1;
a4=gama/beta-1;
a5=dt/2*(gama/beta-2);
a6=dt*(1-gama);
a7=gama*dt;
k1=k+a0*m+a1*c;
for i=2:(length(t))
x=t(i);
pi=(100*x)*(x>=0&x<=0.4)+(80-100*x)*(x>0.4&x<0.8);
pi1=pi+m*(a0*u(i-1)+a2*v(i-1)+a3*a(i-1))+c*(a1*u(i-1)+a4*v(i-1)+a5*a(i-1));
u(i)=pi1/k1;
a(i)=a0*(u(i)-u(i-1))-a2*v(i-1)-a3*a(i-1);
v(i)=v(i-1)+a6*a(i-1)+a7*a(i);
end
ymax=max(u)
将几种方法得到的位移时程曲线绘制到一个图里，Matlab程序如下：
dt=0.1; %时间步长
dtdu=0.01; %duhamel时间步长
t0=0; %起始时间
t2=6.4; %结束时间
t=t0:dt:t2;
tdu=t0:dtdu:t2;
y1=central(dt);
y2=newmark(dt,0.25);
y3=Li(dt);
y4=duhamel(dtdu);
y5=newmark(dt,1/6);
figure
plot(t,y1,’r’,t,y2,’g’,t,y3,’m’,tdu,y4,’k’,t,y5,’b’);
xlabel(‘时间/s’)
axis([0 6.4 -0.1 0.1]) % 调整坐标范围
legend(‘central’,’newmark-常加速度’,’Li’,’duhamel’,’duhamel-线性加速度’);
取dt=0.01s
计算方法
最大位移反应(m)
误差(%)
中心差分
0.0578
0.00
Newmark常加速度
0.0578
0.00
Newmark线性加速度
0.0578
0.00
Duhamel
0.0578
0.00
Fourier
0.0576
-0.346
可认为精确解为0.0578m。
取dt=0.1s
计算方法
最大位移反应(m)
误差(%)
中心差分
0.0610
5.54
Newmark-常加速度
0.0569
-1.56
Newmark-线性加速度
0.0583
0.87
Duhamel
0.0577
-0.17
Fourier
0.0566
-2.076
局部放大如下图：
可见：
（1） Duhamel算法精度最高，Newmark线性加速度法精度次之，Newmark常加速度法与中心差分法精度相当。
（2） FFT算法采样间隔取0.1s时，Nyquist频率=1/(2*dt)=5HZ，满足精度要求的上限频率为2/3*5=3.3HZ，误差较大；
（3） 中心差分法位移反应周期比精确解小，Newmark线性加速度法周期比精确解大，Newmark常加速度法更大，这不难从直观上理解，假设质点运动到接近位移峰值处，中心差分法采用两步外推高估了峰值位移，导致恢复力变大，从而用更少的时间恢复到平衡位置；Newmark常加速度法假设加速度在i和i+1时刻之间为常值，Newmark线性加速度法假设加速度在i和i+1时刻之间线性变化，而事实上加速度时按正弦规律变化，所以Newmark常加速度法和Newmark线性速度法都低估了加速度，从而需要更长的时间恢复到平衡位置。
功率谱方法常用来估计随机反应的均值和均方差，其计算步骤为：
1）确定系统输入的功率谱密度函数Sx(w)；
2）确定结构的复频反应函数H(iw)；
3）计算结构反应的功率谱密度函数Sy(w)；
4）由反应的功率谱密度函数计算自相关函数Ry；
5）计算结构反应的方差
假设本文输入荷载为一随机过程，统计荷载均值为2.5KN，计算位移反应均值、均方差，Matlab程序如下：
dt=0.1;
N=64;
t=[0:N-1]*dt; % 时间序列
xn=(100*t).*(t>=0&t<=0.4)+(80-100*t).*(t>0.4&t<0.8);
subplot(2,2,1)
plot(t,xn) % 绘出原始信号
xlabel(‘时间/s’),title(‘原始信号’)
axis([0 6.4 0 50]) % 调整坐标范围
NN=N; % 执行64点FFT
XN=fft(xn,NN); %
f0=1/(dt*NN); % 基频
f=[0:NN-1]*f0; % 频率序列
A=abs(XN).*abs(XN)/NN; % 幅值序列
subplot(2,2,2);
stem(f,A),xlabel(‘频率/Hz’) % 绘制频谱
axis([0 10 0 500]) % 调整坐标范围
title(‘荷载功率谱’);
%% H(iw)
k=875.5;
m=17.5;
c=35;
zeta=c/m/2/sqrt(k/m);
fn=sqrt(k/m)/2/pi;
HN=ones(1,NN)./(1-(f/fn).*(f/fn)+2i*zeta*(f/fn))/k;
for i=1:NN/2
HN(NN+1-i)=conj(HN(i));
end
UN=HN.*conj(HN).*A;
subplot(2,2,3);
stem(f,UN),xlabel(‘频率/Hz’) %
axis([0 10 0 0.003]) % 调整坐标范围
title(‘位移功率谱’);
un=ifft(UN,NN);
subplot(2,2,4)
plot(t,un) %
xlabel(‘时间/s’),title(‘位移自相关函数’)
axis([0 6.4 -0.0003 0.0003]) % 调整坐标范围
ry=un(1) %Ry(0)
mux=mean(xn) %输入均值
muy=mux*HN(1) %输出均值
sigmay=sqrt(ry-muy*muy) %输出均方差
荷载功率谱、位移功率谱及位移自相关函数如下图所示：
计算得到位移均方差为0.0144m，均值为0.0029m，统计前面Duhamel积分结果（积分间隔0.01s，近似认为是精确解）均方差为0.0147m，均值为0.0029m，可见误差很小。
Duhamel方法和Fourier变换法均基于叠加原理，要求结构体系是线弹性的，当外荷载较大时，结构反应可能进入物理非线性或几何非线性，这时叠加原理将不再适用，此时需要采用时域逐步积分法求解运动微分方程。Newmark方法，特别是β=1/4的无条件稳定格式得到了广泛应用。中心差分法，虽然稳定性略差，但因其所具有的简单、高效的特点得到一系列的应用。对于一些特殊的问题，计算精度的要求有时严于或等于稳定性条件，此时，中心差分法将具有更大的优势。
求解非线性反应时，采用中心差分法无需对计算格式和软件做大的变化，仅是对计算抗力的公式进行改动，其余的与线性反应分析的相同，程序编写方便，便于实现并行计算。SAUSAGE软件即是采用中心差分方法进行非线性动力分析，同时采用了GPU并行技术，大幅提高了计算效率。
在SAUSAGE软件中采用隔震支座来模拟本文单自由度体系，采用瑞丽阻尼，α=2，β=0，中心差分方法，积分步长取0.01s，分析得弹性、弹塑性位移反应如下图所示。
作为对照，在Matlab中，修改弹性分析中心差分法代码，考虑弹塑性， Matlab代码如下：
function [a,u]=centralEP(dt)
% central
m=17.5; %质量
k=875.5; %刚度
c=35;%阻尼比
%dt=0.1; %时间步长
t0=0; %起始时间
t2=6.4; %结束时间
t=t0:dt:t2;
u=t;u(1)=0;u(2)=0;
a=t;
fs=t;
a(1)=0;
fs(1)=0;
k1=m/dt/dt+c/2/dt;
b=m/dt/dt-c/2/dt;
c=2*m/dt/dt;
for i=2:(length(t)-1)
x=t(i);
pi=(100*x)*(x>=0&x<=0.4)+(80-100*x)*(x>0.4&x<0.8);
fs(i)=fs(i-1)+k*(u(i)-u(i-1));
if fs(i)>26.7
fs(i)=26.7;
end
if fs(i)<-26.7
fs(i)=-26.7;
end
pi1=pi-fs(i)+c*u(i)-b*u(i-1);
u(i+1)=pi1/k1;
a(i)=(u(i+1)-2*u(i)+u(i-1))/dt/dt;
end
ymax=max(u)
可见，SAUSAGE分析得到的弹性最大位移为0.0578m，弹塑性最大位移为0.0842m，与Matlab结果完全一致。
1） Duhamel算法精度最高，Newmark线性加速度法精度次之，Newmark常加速度法与中心差分法精度相当。
2） 本文算例时间间隔取0.01s时，几种算法均能取得较为准确的结果；但时间间隔取0.1s时，FFT算法采样误差较大，对于频率大于3.3HZ荷载成分难以准确反应，结果误差较大。
3） 中心差分法位移反应周期比精确解小，Newmark线性加速度法周期比精确解大，Newmark常加速度法更大。
4）功率谱方法常用来估计随机反应的均值和均方差，在已知系统输入的功率谱密度函数时，可以得到随机反应过程的时域强度特征和频域特征，可用于车辆行驶振动分析、风振分析、地震分析等随机过程反应分析中。
5） Duhamel方法和Fourier变换法基于叠加原理，仅适用于线弹性体系。对于非线性体系，需要采用时域逐步积分法求解运动微分方程。Newmark方法，特别是β=1/4的无条件稳定格式得到了广泛应用。中心差分法，虽然稳定性略差，但因其所具有的简单、高效、便于实现并行的特点，得到广泛的应用，对于一些特殊的问题，计算精度的要求有时严于或等于稳定性条件，此时，中心差分法将具有更大的优势。
