|
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<> 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<> 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谱、复频反应函数及时域位移反应如下图所示。
当外荷载较大时，结构反应可能进入物理非线性(弹塑性)，或结构位移较大时，结构可能进入几何非线性，这时叠加原理将不再适用，此时可以采用时域逐步积分法求解运动微分方程。结构动力反应分析的时域直接数值计算方法有：分段解析法、中心差分法、Newmark法、Wilson法、Houbolt法、广义alpha法等，本文仅以最为常用的中心差分法和Newmark方法为例介绍。
中心差分方法用有限差分代替位移对时间的求导（即速度和加速度），在计算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<> 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<> 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<> 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<> 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的无条件稳定格式得到了广泛应用。中心差分法，虽然稳定性略差，但因其所具有的简单、高效、便于实现并行的特点，得到广泛的应用，对于一些特殊的问题，计算精度的要求有时严于或等于稳定性条件，此时，中心差分法将具有更大的优势。
参考文献： 刘晶波，杜修力. 结构动力学[M]. 机械工业出版. 2005年1月 特别感谢： 感谢清华大学刘晶波老师的指导！
|
