Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

What is the unit of 'sd'? #6

Open
farshadrasuli opened this issue Jun 20, 2023 · 1 comment
Open

What is the unit of 'sd'? #6

farshadrasuli opened this issue Jun 20, 2023 · 1 comment

Comments

@farshadrasuli
Copy link

farshadrasuli commented Jun 20, 2023

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
Copy link
Owner

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants