Two-Mass Oscillator: Plotting Amplitudes Over Frequencies in Hertz

[tehsportsmaen]

TL;DR Summary

I try to simulate a two-mass oscillator for a speaker-system and i don't know if my approach makes sense because the physiscal and code-wise details are very complex for me right now.

Idea:
Given a system of two coupled oscillators in which 2 masses are connected to a spring in the middle. Each of the two masses is coupled to another spring on the left and right, which have fixed ends but are not connected to each other. So we have 3 springs, two masses and the springs also have their own damping parameters in addition to the spring constants. Mass m1 is driven by a constant sinusoidal force. So parameters should be: m1, m2, c1, c2, c3, k1, k2, k3, F1. I try to plot the amplitudes over the frequncies in hertz (not omega)

• Define equations and parameters
• Calculate Amplitudes of x1,x2, then for-loop with fast fourier tranformation-calculate in a frequency-range (np.logspace) from 10-1000Hz in 81 steps
• Plot the amplitudes of x1,x2 over the f-range
• Change parameters to see how the system reacts

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from scipy.integrate import odeint
def equations(Y, t, m1, m2, c1, c2, c3, k1, k2, k3, F1, omega):

    x1, v1, x2, v2 = Y

    dx1dt = v1

    dv1dt = (-k1*x1 + k2*(x2 - x1) - c1*v1 + F1*np.sin(omega*t)) / m1

    dx2dt = v2

    dv2dt = (-k3*x2 - k2*(x2 - x1) - ((c2+c3)*v2))/ m2

    return [dx1dt, dv1dt, dx2dt, dv2dt]#SYSTEM

m1, m2 = 0.001, 0.090

c1, c2, c3 = 2.85, 0.5, 1

k1, k2, k3 = 1872.5, 1777.3, 130

F1 = 1.63#TIME AND FREQUENCIES

t = np.linspace(0, 10, 1000)

frequencies = np.logspace(1, 3, 81)def integrate_system(frequency, m1, m2, c1, c2, c3, k1, k2, k3, F1, t):

    omega = 2 * np.pi * frequency

    Y0 = [0, 0, 0, 0]  # initial conditions: [x1, v1, x2, v2]

    solution = odeint(equations, Y0, t, args=(m1, m2, c1, c2, c3, k1, k2, k3, F1, omega))

    return solution[:, 0], solution[:, 2]  #return x1 x2amplitudes_x1 = []

amplitudes_x2 = []for frequency in frequencies:

    x1, x2 = integrate_system(frequency, m1, m2, c1, c2, c3, k1, k2, k3, F1, t)

    #FFT

    fft_x1 = np.fft.fft(x1)

    fft_x2 = np.fft.fft(x2)

    #compute amplitudes

    amplitude_spectrum_x1 = np.abs(fft_x1) / len(fft_x1)

    amplitude_spectrum_x2 = np.abs(fft_x2) / len(fft_x2)

    #index of fft frequencies

    frequency_index = np.argmin(np.abs(frequency - np.fft.fftfreq(len(t), d=t[1]-t[0])))

    #amplitudes of driving frequencies

    amplitudes_x1.append(amplitude_spectrum_x1[frequency_index])

    amplitudes_x2.append(amplitude_spectrum_x2[frequency_index])pd.set_option('display.max_rows', None)

df = pd.DataFrame({'Frequency_Hz': frequencies,'Amplitude_m1': amplitudes_x1,'Amplitude_m2': amplitudes_x2})

print(df)plt.figure(figsize=(12, 6))

plt.loglog(frequencies, amplitudes_x1, label='Amplitude of m1', color='red')

plt.loglog(frequencies, amplitudes_x2, label='Amplitude of m2', color='blue')

plt.xlabel('Frequency (Hz)')

plt.ylabel('Amplitude')

plt.title('LS FREQUENCIES')

plt.legend()

plt.grid(True)

plt.show()

The result looks strange but maybe it makes sense.Would could be potentially difficult with my approach?=> (timearray, integration, FFT, calculation, initial conditions, parameter interpretation)I am happy to get advices, comments and remarks on this.

 

[scottdave]

I have not used odeint in python before, so I'd have to read up on that documentation.

When you put your equations into code, it looks like you account for 2 damping coefficients for m₂ but only 1 damping coefficient for m₁. Try drawing free body diagrams for each mass.

 

[scottdave]

Also I think when you have the formula for

dv2dt = (-k3*x2 - k2*(x2 - x1) - ((c2+c3)*v2))/ m2

Is the bolded portion solely dependent on v2, or is it the relative speed? The free body diagrams should help with this.

 

[Haborix]

It might be easier to first try just doing oscillators with no damping and no driving. That can be solved analytically and you can compare to numerical results. This way you can check the barebones of your code is working properly.

 

[tehsportsmaen]

Also I think when you have the formula for

dv2dt = (-k3*x2 - k2*(x2 - x1) - ((c2+c3)*v2))/ m2

Is the bolded portion solely dependent on v2, or is it the relative speed? The free body diagrams should help with this.

def equations(Y, t, m1, m2, c1, c2, c3, k1, k2, k3, F1, omega):
    x1, v1, x2, v2 = Y
    dx1dt = v1
    dx1dt = v1
    dv1dt = (-k1*x1 + k2*(x2 - x1) - c1*v1 - c2*(v1 - v2) + F1*np.sin(omega*t)) / m1

    dx2dt = v2
    dv2dt = (-k3*x2 - k2*(x2 - x1) - c3*v2 - c2*(v2 - v1)) / m2
    return [dx1dt, dv1dt, dx2dt, dv2dt]

Is this better?

 

[tehsportsmaen]

This is the Output now