2
$\begingroup$

I have my program for mass spring system made with Euler method...and I can't menage to obtain envelope of the BLUE damping curve. Can you help me?

enter image description here

#****************************************************************************************************************************************************************************************
# 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
# ****************************************************************************************************************************************************************************************

I did it but I am afraid that I wonted something else - I wantef logarithmic damping increment and I have to figure out how to draw from the elvelope the information about level of damping. Nevertheless - thank you.

enter image description here

$\endgroup$
  • 1
    $\begingroup$ so, this looks nice, but what have you tried so far? It's hard to really help you if we don't know where you're stuck. $\endgroup$ – Marcus Müller Nov 26 '19 at 22:00
  • 1
    $\begingroup$ Have you tried the hilbert envelope? $\endgroup$ – Gideon Genadi Kogan Nov 26 '19 at 23:12
  • $\begingroup$ I have tried Hilbert envelope - it doesn't work :( $\endgroup$ – L. Flis Nov 27 '19 at 17:29
  • $\begingroup$ The programs works - ok. I just need the uper and down envelope. I don't know how to do it - ok. Then I looked for some aproximation in Python with "signa envelope" subject etc. In the Internet in most cases there was Hilbert solution answer. Then I implemente to my program but it didn't work...Additionaly I didn't find any example of working Hilbert solution for real data...or if the solution seems good I don't see the data. My curve is just displacement in time function and when I put it to envelope = abs(hilbert(displacement)) and then plot (time, envelope) it doesn't work $\endgroup$ – L. Flis Nov 27 '19 at 17:42
1
$\begingroup$

So, you have numerically solved a second order differential equation and plotted its result. Now it seems to be a damped sinusoidal response of the type

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

for some constants $\alpha$, $\omega_0$ and $K$. And you want to compute the damping parameter $\alpha$.

A crude approximation to $\alpha$ can be obtained by the following. Consider two points, one at an arbitrary time $t_0$ and the other one period later $t_0 + T_0$, on the curve. You don't know the period ? It can be approximately decuded by looking at two consecutive peaks of the oscillations. Note that those peaks are slightly misplaced due to being weigted by the exponential; nevertheless it will be ok for a crude approximation. Then the values at those two points will be $A_1$ and $A_2$; then it can be seen that

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

Note that you have already measured the period $T_0$ in the first step; hence you can find the damping parameter $\alpha$ to be

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

If you want to estimate it more accuretely, you might consider using different estimators for it. A Kalman filter can also be designed to esitmate it, but that's a lot more complex.

$\endgroup$
  • $\begingroup$ finding two peaks is not that difficult, but yes, this is not an automatic approach. may be an approximate solution that could be easily computed can be also found... $\endgroup$ – Fat32 Nov 27 '19 at 17:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.