arkottke / pyrotd

Python scripts to compute the rotated response spectrum.
MIT License
21 stars 14 forks source link

What is the unit of 'sd'? #6

Open farshadrasuli opened 1 year ago

farshadrasuli commented 1 year ago

Hello everyone,

I found that the (pseudo-)velocity values (psv, or sv) calculated by the pyrotd.calc_spec_accels(time_step, accel_ts, osc_freqs, osc_damping=0.05, max_freq_ratio=5, osc_type='psv') function must be multiplied by $g$ (the acceleration of gravity). In other words, the psv values are in units of g-s.

I obtained the spectral displacement (sd) of a record by the pyrotd.calc_spec_accels(time_step, accel_ts, osc_freqs, osc_damping=0.05, max_freq_ratio=5, osc_type='sd') function in the following script, but I can't figure out what the unit of these values is. I calculated the spectral displacement (sd) with pseudo-spectral acceleration $$Sd = pSa \cdot g \cdot \frac{T^2}{4\pi^2}$$ and there is a huge error between the calculated values by formula and the computed values by pyrotd library. There is a plot to illustrate of this huge error.

import numpy as np
import matplotlib.pyplot as plt
import pyrotd

#%% system of units
m = 1.0
s = 1.0
kg = 1.0
N = kg*m/s/s

g = 9.80665*m/s/s

#%% read record
record = 'RSN8883_14383980_13849360.AT2'

record_content = open(record, mode='r').readlines()

# Flag indicating dt is found and that ground motion
# values should be read -- ASSUMES dt is on last line
# of header!!!
flag = 0

record_samples = []
for line in record_content:
    if line == '\n':
        # Blank line --> do nothing
        continue
    elif flag == 1:
        words = line.split()
        for word in words:
            record_samples.append(float(word))
    else:
        words = line.split()
        if words[0] == 'NPTS=':
            record_NPTS = words[1]
            record_DT = float(words[3])
            flag = 1

#%% compute spectral parameteres
f = np.logspace(-1, 2, 91)
T = 1/f
xi = 0.05
pyrotd_pSa = pyrotd.calc_spec_accels(record_DT, record_samples, f, xi, max_freq_ratio=5, osc_type='psa')

pyrotd_Sd = pyrotd.calc_spec_accels(record_DT, record_samples, f, xi, max_freq_ratio=5, osc_type='sd')

formulated_Sd = []
for index, psa in enumerate(pyrotd_pSa.spec_accel):
    formulated_Sd.append( psa*g * (T[index]**2) / (2 / np.pi)**2 )
formulated_Sd = np.asarray(formulated_Sd)

#%% plot
fig, axs= plt.subplots(layout='constrained')
fig.suptitle('Displacement spectra of ' + record)
axs.plot(T, pyrotd_Sd.spec_accel, label='Sd')
axs.plot(T, formulated_Sd, label=r'Sd = $pSa \cdot g \cdot \frac{T^2}{4\pi^2}$')
axs.set_xlabel('Period, s')
axs.set_ylabel('Spectral displacement')
axs.legend()
plt.show()

RSN8883_14383980_13849360 AT2 displacement spectra

I figured out that by multiplying sd values calculated by pyrotdlibrary by $g^3$, the error will be reduced significantly.

#%% plot
fig, axs= plt.subplots(layout='constrained')
fig.suptitle('Displacement spectra of ' + record)
axs.plot(T, pyrotd_Sd.spec_accel, label='$Sd$')
axs.plot(T, pyrotd_Sd.spec_accel*(g**3), label='$Sd \cdot g^3$')
axs.plot(T, formulated_Sd, label=r'$Sd = pSa \cdot g \cdot \frac{T^2}{4\pi^2}$')
axs.set_xlabel('Period, s')
axs.set_ylabel('Spectral displacement')
axs.legend()
plt.show()

RSN8883_14383980_13849360 AT2 displacement spectra

arkottke commented 1 year ago

That's correct, there is no unit conversion in the calculation of the response spectra.