In [1]:
%matplotlib inline

import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (12.0, 3.2)
import numpy as np

from IPython.display import display, HTML

Vibration Insulation and Numerical Integration

Data

In [2]:
M = 15.2E3 # kg
P0 = 800.0 # N
w0 = 2*np.pi*8.0 # rad/s
t0 = 5.0 # s
TR = 0.25

Design of Suspension System

For a specific value of $\zeta$ we want to determine the maximum value of $K$ for which the transmissibility ratio is not greater than TR.

Starting from the definition of the transmissibility ratio.

$\text{TR} \ge \frac{\sqrt{1+(2\zeta\beta)^2}}{\sqrt{(1-\beta^2)^2+(2\zeta\beta)^2}} \quad\rightarrow\quad \text{TR}^2((\beta^2-1)^2 +(2\zeta\beta)^2)) - 1 - (2\zeta\beta)^2 \ge 0$

Solving for $\beta^2$, observing that $\beta^2=M\omega^2_0/K$ we can find directly the maximum stiffness needed to not exceed the maximum steady-state transmitted force.

In [3]:
def max_K(M, z, w0, TR):
    T =TR**2
    fz2 = 4*z*z
    fz4 = fz2**2
    K = (2*T + fz2 - T*fz2 - np.sqrt((1-2*T+T*T)*fz4 + 4*fz2*(T-T*T) + 4*T)) / (2*T-2)
    return M*w0*w0*K

K01 = max_K(M, 0.01, w0, TR)
K10 = max_K(M, 0.10, w0, TR)

The damping coefficient is $C=2\zeta\sqrt{KM}$

In [4]:
C01 = 2*np.sqrt(M*K01)*0.01
C10 = 2*np.sqrt(M*K10)*0.10

The required output

In [5]:
print('The suspension system for z=0.01')
print('\tK =%E N/m,\tC =%E Ns/m.'%(K01, C01))
print('The suspension system for z=0.10')
print('\tK =%E N/m,\tC =%E Ns/m.'%(K10, C10))

The suspension system for z=0.01
	K =7.675163E+06 N/m,	C =6.831178E+03 Ns/m.
The suspension system for z=0.10
	K =7.131325E+06 N/m,	C =6.584714E+04 Ns/m.

Loading during transient and after

In [6]:
factor = 4
t = np.linspace(0, 6, 600*factor+1)
a = t/t0
f = w0*t0*np.where(t<t0, a*a*(1-a/3), a-1/3)
plt.plot(t,f)
plt.xlabel('$t/s$')
plt.ylabel('Phase/rad, $\\phi(t)$')
plt.grid();
In [7]:
w = w0*np.where(t<t0, a*(2-a), 1)
plt.plot(t, w)
plt.xlabel('$t/s$')
plt.ylabel('Angular velocity/(rad/s), $\\dot\\phi(t)=\\omega(t)$')
plt.grid();
In [8]:
angacc = w0*np.where(t<t0, 2*(1-a), 0)/t0
plt.plot(t, angacc)
plt.xlabel('$t/s$')
plt.ylabel('Angular acceleration, $\\ddot\\phi(t)$')
plt.grid();
In [9]:
p = P0 * (w*w*np.sin(f)-angacc*np.cos(f)) / w0**2
plt.plot(t, p)
plt.xlabel('$t/s$')
plt.ylabel('$p(t)/{}$N')
plt.grid();
In [10]:
plt.plot(t, w*w*np.sin(f)*P0/w0/w0, label='$\\omega^2\\sin\\phi$')
plt.plot(t, angacc*np.cos(f)*P0/w0/w0, label='$\\dot\\omega\\cos\\phi$')
plt.legend()
plt.xlim((0, 2))
plt.ylim((-300, 300))
plt.title('Relative importance of different contributions to unbalanced force');

Numerical Integration

Straightforward implementation, we print also the peak value of the transmitted force

In [11]:
h = 0.01/factor
for m, c, k, zeta in ((M, C01, K01, '0.01'), (M, C10, K10, '0.10')):
    k_ = k + 3*c/h + 6*m/h/h
    c_ = 3*c + 6*m/h
    m_ = 3*m + c*h/2
    x0, v0 = 0, 0
    x, v = [x0], [v0]
    for p0, p1 in zip(p[:-1], p[1:]):
        a0 = (p0-c*v0-k*x0)/m
        dp = p1 - p0 + c_*v0 + m_*a0
        dx = dp/k_
        x0 = x0 + dx
        v0 = 3*dx/h - 2*v0 - a0*h/2
        x.append(x0), v.append(v0)
    x, v = np.array(x), np.array(v)
    f = k*x+c*v ; f0 = np.max(np.abs(f))
    display(HTML('''
       <h3>$\\zeta=%s$</h3>
       Peak value of transmitted force $f_0=%f\,{}$N.'''%(zeta, f0)))
    plt.plot(t, f)
    plt.xlabel('$t/{}$s')
    plt.ylabel('$f(t)/{}$N')
    plt.grid()
    plt.show();

$\zeta=0.01$

Peak value of transmitted force $f_0=1359.497110\,{}$N.

$\zeta=0.10$

Peak value of transmitted force $f_0=659.983619\,{}$N.

Alternative

If you don;t want to solve a quadratic equation, you can use a numerical root finder, that can be initialized using the value $\beta^2=1+1/\text{TR}$, valid for an undamped system.

In [12]:
from scipy.optimize import newton
z = 0.10; 
f = lambda b2: (1+4*z*z*b2)/((1-b2)**2+4*z*z*b2)-TR**2
b2 = newton(f, (TR+1)/TR)
print('K|z=0.10 =', M*w0*w0/b2,'N/m')

K|z=0.10 = 7131324.519944637 N/m