Random vibration calculation and simulation verification of single degree of freedom system

1、 Reference resources Single degree of freedom vibration calculation

2、 exercises

Suppose that the mean value of the excitation force received by a degree of freedom system is 0、 Self spectrum is S0=1 Ideal white noise stationary process . Find the mean square value and self spectrum of the system response （ Power spectral density function , Represents the distribution density of the mean square value of the random process in the frequency ）.
set up m=1,k=(2*pi)**2, Damping coefficient zeta=0.16, That is, the natural frequency of the structure is 1Hz, Fast decay system .
When the excitation force is random , Similar to the confirmation , We need to use the impulse response function of the system h(t) Or complex frequency response function H(w).
Impulse response function of the system h(t) See the picture below ：

The self spectrum of ideal white noise excitation force is shown in the figure below ：

The self spectrum of displacement response is shown in the figure below ：

3、python Computational code

import matplotlib.pyplot as plt
import numpy as np

from numpy import pi,sqrt,cos,sin,arctan,exp

plt.rcParams['font.sans-serif']=['SimHei'] # Used to display Chinese labels normally 
plt.rcParams['axes.unicode_minus']=False # Used to display negative sign normally 

##########begin: Calculate the vibration of single degree of freedom ###############
k=(2*pi)**2
m=1
wn=sqrt(k/m)
s_=0.16
wd=sqrt(1-s_**2)*wn
c=2*s_*wn*m
delta=s_*wn

# The mean value of the system excitation force is 0 Self spectrum is S0 Ideal white noise stationary process ：
S0=1
# Impulse response function of the system ：
h_t=lambda t:1/(m*wd)*exp(-delta*t)*sin(wd*t)
# Mean square value of response ：
phi2_x=S0/(2*c*k)
# Self spectrum of response ：
Sxx_w=lambda w:S0/((k-m*w**2)**2+c**2*w**2)

ti=np.array(np.linspace(0,4,1000,False))#0:0.002:2;
h_ti=h_t(ti)
wi=np.array(np.linspace(0,100,100,False))#0:0.002:2;
Sxx_wi=Sxx_w(wi) # One side of the bilateral power spectrum 
Gxx_wi=2*Sxx_wi # Unilateral power spectrum 

plt.figure
fig_mpl, ax =plt.subplots(3,1,figsize=(5,3), facecolor='white')
ax[0].plot(ti,h_ti,label=' Impulse response function of the system ')
ax[0].set_title("h(t)")
ax[0].set_xlabel("t")
ax[0].legend()
ax[1].plot(wi,np.ones(np.size(wi))*S0*2,label=' Input force white noise self spectrum （ unilateral ）')
ax[1].set_title("Gff(w)")
ax[1].set_xlabel("w")
ax[1].legend() # Appoint legend
ax[2].plot(wi,Gxx_wi,label=' Output displacement self spectrum （ unilateral ）; Mean square value :%f'%phi2_x)
ax[2].set_title("Gxx(w)")
ax[2].set_xlabel("w")
plt.legend() #plt.legend Mark only the last picture 
plt.subplots_adjust(left=None, bottom=None, right=None,
top=None, wspace=1, hspace=1)
plt.show()

4、 Simulation verification to be continued