Calculates the random response of a single degree of freedom subject to a PSD

RR_SDOF.py

import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy.interpolate import interp1d
import math

PSD = np.array([[0.001,1.0e-8],
                [5.0, 0.001],
                [25.0, 0.03],
                [30.0, 0.03],
                [89.0, 0.321],
                [90.0, 1.0],
                [260.0, 1.0],
                [261.0, 0.03],
                [359.0, 0.03],
                [360.0, 0.5],
                [520.0, 0.5],
                [540.0, 0.25],
                [780.0, 0.25],
                [781.0, 0.03],
                [2000.0, 0.03],
                [2.0e6,1.0e-8]])

def W_F(freq):
    '''
    A line connecting two points in a log-log plot are exponential
    '''
    logPSD = np.log10(PSD)
    logW_F = interp1d(logPSD[:,0], logPSD[:,1])
    return np.power(10,logW_F(np.log10(freq)))
##    w_f = []
##    for f in freq:
##        index = np.searchsorted(PSD[:,0], f)
##        if index <= 0:
##            w_f.append(PSD[:,1][0])
##        elif index + 1>= PSD.shape[0]:
##            w_f.append(PSD[:,1][-1])
##        else:
##            x0 = PSD[:,0][index-1]
##            F0 = PSD[:,1][index-1]
##            x1 = PSD[:,0][index]
##            F1 = PSD[:,1][index]
##            w_f.append(F0*(f/x0)**(math.log(F1/F0)/math.log(x1/x0)))
##    return np.array(w_f)

def H(r, zeta):
    '''
    Gain function for base excited system
    '''
    return (r**2)/np.sqrt((1-r**2)**2+(2.0*zeta*r)**2)

def W_X(f, f_n, zeta):
    r = f/f_n
    return (H(r,zeta))**2*W_F(f)

if __name__ == '__main__':
    f = np.arange(5.0, 2000.0, 0.1)
    f_n = 153.0
    zeta = 0.10
    y = W_X(f, f_n, zeta)
    y_rms = math.sqrt(np.trapz(y, x = f))

    print("a_rms:%s" % (y_rms,))

    #plt.plot(f,W_X(f, f_n, zeta),'-',f,H(f/f_n,zeta),'--',f,W_F(f),'-.')
    fig = plt.figure()
    ax = fig.add_subplot(3,1,1)
    line, = ax.plot(f, W_X(f, f_n, zeta), color='blue', lw=2)
    ax.set_yscale('log')
    ax.set_xscale('log')
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('W_X')
    ax.grid(True)
    ax = fig.add_subplot(3,1,2)
    line, = ax.plot(f, H(f/f_n,zeta), color='blue', lw=2)
    ax.set_yscale('log')
    ax.set_xscale('log')
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('H')
    ax.grid(True)
    ax = fig.add_subplot(3,1,3)
    line, = ax.plot(f, W_F(f), color='blue', lw=2)
    ax.set_yscale('log')
    ax.set_xscale('log')
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('W_F')
    ax.grid(True)
    plt.show()