Skip to content

What is the unit of 'sd'? #6

@farshadrasuli

Description

@farshadrasuli

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

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions