单自由度体系传递函数的验证问题

grandpollux

请教一个基础问题：
    已知单自由度粘滞阻尼体系相关参数值：
          m=0.2533,
          k=10,
          zeta=0.05,
         wn=6.283,
         c=zeta*(2*m*wn)
激振力为：
         f(t)=p0*sin(pi/0.6*t) 其中p0=10.0,
初始条件为：
          位移u0=0
          速度v0=0
则可由相关解析计算公式计算此体系的位移响应u(t),
这个解析公式可在一般教科书中找到。
    在Matlab中取dt=0.001,计算总时长T=2.047秒,则
可得到作用力向量f(t)与位称向量u(t). 应用tfe函数，可
知其传递函数Tf及相应的频率向量omiga为：
    [Tf,omiga]=tfe(f(t),u(t))

    另一方面，对于单自由度体系，其传递函数可以解析得到：
    Tf_analytical=1.0/(k+i.*omiga*c-omiga.^2*m)
式中，omiga是tfe函数计算得到的omigat向量,i是虚单位。
    检查Tf与Tf_analytical这向个复向量可以发现并不相等，
请问这是什么原因？

[ 本帖最后由 grandpollux 于 2010-4-4 00:45 编辑 ]

mao

如果你肯定你计算的位移响应u(t)是正确的话，最好传上一张计算的传函数的幅值与相位图与解析值的幅值与相位图进行对比，大家好分析原因

wanyeqing2003

这个图里没看出传递函数。

grandpollux

文档里有三张图片，中间那张图片是关于传递函数（由matlab中tfe函数计算得到）的，分别画了传递函数的实部、虚部和幅值。

mao

看了你的U(t)，问题就出在这，显然你的U(t)是一个受迫振动的响应，你在用tfe计算时，计算总时长T=2.047秒，而tfe是根据输入输出的互谱除以输入的自谱得到的，因此它会认为U(t)是以2.047秒为周期的周期信号，这样明显会与你真实的U(t)是出入的。如何解决这个问题，　我建议不要用强迫振动的结果来计算，只用自由振动的结果来验证。如何方便，给出你自由振动的验证最好！

grandpollux

感谢 mao 的回复。
不过我不太明白的是：如果采用自由振动来验证的话，那么意味着输入为零，此时如何根据输入与输出来计算传递函数呢？

mao

不好意思，昨天没有细看你的模型，回复不周到，刚才仔细看了一下你的模型，你的自由振动同期为6.283s,而强迫激励周期约5/6s，你取的U(t)长度为2.047s,显然时间过短，自由振动引起的瞬态响应还没有大部分消失，象我前面提到的，tfe计算时认为U(t)是以2.047秒为周期的周期信号，所以会出现较大误差。要想得到正确的验证结果，你可以取几个自由振动周期后的响应U(t)（具体取几个自由振动周期取决你模型所用的阻尼大小，阻尼大就可以用较小的周期，如果阻尼小就要多取几个周期后的值才行，反正要保证自由振动引起的瞬态响应至少要大部分消失才行），并且保证截取的长度为强迫激励周期的整数倍，这样就应该不会出现验证误差了。（你的模型中，自由振动的周期过长，建议缩短，即取较小的k/m值）

grandpollux

我的计算中取wn=6.283, 亦即自振周期T＝2*pi/wn=1.0秒。
强迫振动中w=pi/0.6, 周期：2*pi/w=1.2秒。
我刚试过了，采用增长计算时间的办法处理还是有问题。

grandpollux

下面我贴上代码：
% basic known conditions:
m=0.2533;
k=10;
wn=sqrt(k/m);
zeta=0.05;
c=zeta*(2*m*wn);
wd=wn*sqrt( 1-zeta^2 );
dt=0.001;
u0=0;
v0=0;
p0=10;
w=pi/0.6;
T=12;
%----------------------------------------------------------------------------
% the following parameters for solving analytical solutions
CC=p0/k*( 1-(w/wn)^2 )/( (1-(w/wn)^2)^2+(2*zeta*(w/wn))^2 );
DD=p0/k*( -2*zeta*(w/wn))/( (1-(w/wn)^2)^2+(2*zeta*(w/wn))^2 );
AA=u0-DD;
BB=1.0/wd*(v0+zeta*wn*AA-w*CC);
%----------------------------------------------------------------------------
time=0:dt:T;
force=sin(w.*time)*p0;
% analytical solution:
u_analytical=exp(-zeta.*wn.*time).*(AA.*cos( wd.*time)+BB.*sin( wd.*time))+CC.*sin(w.*time)+DD.*cos(w.*time);

%  ***********plot ***u(t)*** graph***************
figure
plot(time, u_analytical,'.')
axis tight
xlabel('Time(s)')
title('U(z) history curve by analytical method')
grid on
%  ***********************************************

sz=length(time);
% using tfestimate() fuction to solving transfer function Hn_tfe,
Fc=1.0/dt;
window=hanning(sz);
[Tfn_tfe,f]=tfestimate(force,u_analytical,window,sz*0.5,sz,Fc);   % f's unit is Hz

%  ******************plot ****tranfer function****graph***************
figure
subplot(1,3,1)
plot( f, real(Tfn_tfe),'r.')
title('transfer function--real part')
grid on
xlabel('f(Hz)')
subplot(1,3,2)
plot(f,imag(Tfn_tfe),'b.')
title('transfer function--imag part')
grid on
xlabel('f(Hz)')
subplot(1,3,3)
plot(f,abs(Tfn_tfe),'k.')
title('transfer function--amplitude')
grid on
xlabel('f(Hz)')
%  ********************************************************************

omiga=2*pi*(1/dt)*(0:sz/2-1)/sz;     % omiga's unit is rad/s
Hn=p0./(-omiga.^2.*m+j.*omiga.*c+k);  

%  **************plot ***   Comparing  ****graph**********************
headpos=100;
intval=1;
omiga=omiga/2/pi;   % omiga's unit is Hz
figure
subplot(1,3,1)
plot(omiga(headpos:intval:end),real(Hn(headpos:intval:end)),'r',f(headpos:intval:end),real(Tfn_tfe(headpos:intval:end)),'bo')
grid on
xlabel('f(Hz)')
legend('analytical method','numerical method')
title('Comparing two kind of transfer function value--real part')
subplot(1,3,2)
plot(omiga(headpos:intval:end),imag(Hn(headpos:intval:end)),'r',f(headpos:intval:end),imag(Tfn_tfe(headpos:intval:end)),'bo')
grid on
xlabel('f(Hz)')
legend('analytical method','numerical method')
title('Comparing two kind of transfer function value--imag part')
subplot(1,3,3)
plot(omiga(headpos:intval:end),abs(Hn(headpos:intval:end)),'r',f(headpos:intval:end),abs(Tfn_tfe(headpos:intval:end)),'bo')
grid on
xlabel('f(Hz)')
legend('analytical method','numerical method')
title('Comparing two kind of transfer function value--amplitude')
%  ********************************************************************

mao

你没有看懂我的意思，先把这条语句改成：time=100:dt:（100+T）看看结果如何？

grandpollux

我刚试过了，仍然有问题。理论解与数值解仍有明显不同。
但不论结果如何，非常感谢mao的回复。这个问题我困扰很久了，一直找不到原因。

mao

呵，看来做事不认真还真不行呢，不过原因我想还是我说的那样，只是昨天我没有仔细看你的力的表达方式，这样，你把相应的部分再这样改：
time=0:dt:（100+T）;    % 有改动
force=sin(w.*time)*p0;
% analytical solution:
u_analytical=exp(-zeta.*wn.*time).*(AA.*cos( wd.*time)+BB.*sin( wd.*time))+CC.*sin(w.*time)+DD.*cos(w.*time);
u_analytical=u_analytical((100/dt):((100+T)/dt));      % 增加部分，目的是为了去掉瞬态响应的影响

我在网吧，没有装MATLAB，你去运行看看结果如何，应该没有问题了，如果效果还不太满意的话，把上面所有的100改成更大的值，如1000。

grandpollux

谢谢mao的耐心指点。

LPJ870626

请教高手
用Matlab画单自由度强迫振动（偏心轮引起的强迫振动）的bode图和幅频、相频图。
幅频、相频图用等角度采样实现