信息处理 - 如何获得我的质量弹簧阻尼曲线的包络？ - 吾爱随笔录

如何获得我的质量弹簧阻尼曲线的包络？

信息处理信号分析

2022-02-19 13:53:40

我有我用欧拉方法制作的质量弹簧系统程序......我无法获得蓝色阻尼曲线的包络线。你能帮助我吗？

#****************************************************************************************************************************************************************************************
# Damped spring-mass oscillator
# ****************************************************************************************************************************************************************************************
from pylab import *
from scipy.signal import hilbert,chirp

g =  0          #grawitacja [m/s2] bez grawitacji 0
g = -9.80665       #grawitacja [m/s2] bez grawitacji 0
m = 0.4532      #masa ciężaru [kg]
k = 875.60      #sztywność [N/m]
# c = (2*m)*sqrt(k/m)
c = 0
c_kr = 2*sqrt(k*m)
c = 2*sqrt(k*m)
c = 5.0         #tlumienie [Ns/m]

omega=sqrt(k/m)
f=omega/(2*pi)
T=1/f
T=2*pi/omega
f1=1/T
gamma=c/c_kr
p=c/(2*m)

print('dekrement tłumienia',gamma,'  [s]')
print('tłumienie', c)
print('tłumienie krytyczne',c_kr)
print('omega czestosc własna',omega,'  [rad]')
print('czestotliwosc',f,'  [Hz=1/s]')
print('czestotliwosc',f1,'  [Hz=1/s]')
print('okres',T,'  [s]')
# print('e',e,'  [s]')

u_pocz=-0.0254  #ugięcie, przemieszczenie wstępne sprężyny [m]
u_pocz= 0.0     #ugięcie, przemieszczenie wstępne sprężyny [m]
v_pocz= 0.0     #prędkość w [m/s]
t_pocz= 0.0     #czas w [s]
env_pocz= 0.0   #czas w [s]

dt=0.001        #przyrost czasu [s]
t_kon=0.15      #czas końca obliczeń [s]
t_kon=1.00      #czas końca obliczeń [s]

print('czas',t_kon,'  [s]')

u=u_pocz        #przemieszczenie [m]
v=v_pocz        #prędkość w [m/s]
t=t_pocz        #czas w [s]
env=env_pocz    #czas w [s]

#To mogę policzyć bo wiem, ugięcie, wydłużenie sprężyny w [m] - a dlaczego tak ?
#bo k*u=G (siła ciężkości) czyli k*u=m*g dlatego u=m*g/k

# u_ugiecie_wlasne=m*g/k-u_pocz                 #ugięcie przemieszczenie własne w [m]
# a=-g-(k/m)*(u-u_ugiecie_wlasne)-c*v*abs(v)/m  #przyspieszenie pocztkowe w [m/s2] ale potrzebne to bedzie dopierow później
u_ugiecie_wlasne=m*(g)/k                        #ugięcie przemieszczenie własne w [m] bez LOAD_BODY
print('u_ugiecie_wlasne',u_ugiecie_wlasne,'  [s]')

przechowalnia_u=[]
przechowalnia_v=[]
przechowalnia_a=[]
przechowalnia_F=[]
przechowalnia_t=[]
przechowalnia_env=[]


# ****************************************************************************************************************************************************************************************
#                                                           Obliczenia zasadnicze
# ****************************************************************************************************************************************************************************************

while (t<t_kon):

    a=-g-(k*(u-u_ugiecie_wlasne)/m)-(c*v/m)+g   # dv/dt=a=(-g-k*(u+u_ugiecie_wlasne)/m-c*v/m) z LOAD_BODY czyli przyspieszenie ziemskie *g

    v=v+a*dt                                    # v=v+dt*dv/dt
    u=u+v*dt                                    # u=u+dt*du/dt
    F=(-g-(k*(u-u_ugiecie_wlasne)/m)-(c*v/m))*m                                     
    t=t+dt

    przechowalnia_u.append(u)
    przechowalnia_v.append(v)
    przechowalnia_a.append(a)
    przechowalnia_F.append(F)
    przechowalnia_t.append(t)
    przechowalnia_env.append(env)


# ****************************************************************************************************************************************************************************************
#                                                           Obliczenia zasadnicze
# ****************************************************************************************************************************************************************************************



# ****************************************************************************************************************************************************************************************
#                                                                Grafika
# ****************************************************************************************************************************************************************************************

# coding: utf-8
from matplotlib.font_manager    import FontProperties
# okno = get_current_fig_manager()

def quit_figure(event):
    if event.key == 'escape':
       close(event.canvas.figure)
    if event.key == 'f10':
       savefig('0.png',dpi=150)

# rcParams.update({'font.size': 24})
rcParams['font.size'] = 24                      #set the value globally
rcParams['axes.linewidth'] = 2                  #set the value globally
rcParams['toolbar'] = 'None'
font = {'family':'ISOCPEUR','weight':'normal','color':'black'}
# font = {'family':'ISOCPEUR','weight':'normal','color':'black','size':10}
fig=figure(num=None,  frameon='False', figsize=(16, 7), facecolor='w')
# title('Identification of Johnson-Cook constitutive equations in terms of FEM simulation\n$\mathrm{Y=}$',font)
ylabel('Przemieszczenie [m]')
xlabel('Czas [s]')
plot(przechowalnia_t, przechowalnia_u, linewidth=4, color='b',label='u - przemieszczenie, [m]')
# plot(analytical_signal, linewidth=2, color='b',label='u - przemieszczenie, [m]')
# plot(amplitude_envelope, linewidth=4, color='r',label='u - przemieszczenie, [m]')

# plot(przechowalnia_t, przechowalnia_a, linewidth=4, color='r',label='a - przyspieszenie, [m/s$^2$]'')
# plot(przechowalnia_t, przechowalnia_v, linewidth=4, color='r',label='v - prędkość, [m/s2]')
# plot(przechowalnia_t, przechowalnia_F, linewidth=4, color='r',label='F - siła, [N]')
axvline(0,linestyle='--',linewidth=2,color='r')
axhline(0,linestyle='--',linewidth=2,color='r')
hlines(u_ugiecie_wlasne,t_pocz,t_kon,linestyle='-',linewidth=4,color='g',label='ugięcie własne w [m]')
legend(loc='upper right')

quit = gcf().canvas.mpl_connect('key_press_event', quit_figure)

show()

# ****************************************************************************************************************************************************************************************
#                                                                Grafika
# ****************************************************************************************************************************************************************************************

我做到了，但我担心我不习惯其他东西 - 我想要对数阻尼增量，我必须弄清楚如何从 elvelope 中提取有关阻尼水平的信息。尽管如此 - 谢谢。

1个回答

因此，您已经对二阶微分方程进行了数值求解并绘制了其结果。现在它似乎是类型的阻尼正弦响应

y (t) = K e^{- α t} \cos (ω_{0} t + ϕ)

$y(t) = K e^{-\alpha t} \cos(\omega_0 t + \phi)$

对于一些常数，和。你想计算阻尼参数。 $\alpha$ $\omega_0$ $K$ $\alpha$

可以通过以下方式获得对的粗略近似。考虑曲线上的两个点，一个在任意时间，另一个之后的时间段。你不知道周期？它可以通过观察两个连续的振荡峰值来近似推断。请注意，由于被指数加权，这些峰略有错位；尽管如此，粗略的近似还是可以的。那么这两个点的值将是和；那么可以看出 $\alpha$ $t_0$ $t_0 + T_0$ $A_1$ $A_2$

\frac{| A_{1} |}{| A_{2} |} = \frac{K e^{- α t_{0}}}{K e^{- α (t_{0} + T_{0})}} = e^{α T_{0}}

$\frac{ |A_1| }{|A_2| } = \frac{ K e^{-\alpha t_0} }{K e^{-\alpha (t_0 + T_0) }} = e^{\alpha T_0}$

请注意，您已经在第一步中因此你可以找到阻尼参数为 $T_0$ $\alpha$

α = \frac{1}{T_{0}} \ln (| A_{1} | / | A_{2} |)

$\alpha = \frac{1}{T_0} \ln ( |A_1| / |A_2| )$

如果您想更准确地估计它，您可以考虑使用不同的估计器。也可以设计卡尔曼滤波器来估计它，但这要复杂得多。

其它你可能感兴趣的问题

上一篇平行形式的过滤器设计下一篇应用带通滤波器后信号开始处的大尖峰