Click to expand the full Transfer Function Example code
.. code-block:: python
# ############################################################################################################
# Title: Transfer Function Calculation
# Description: This script calculates the transfer function of a soil layer and applies it to a given acceleration record.
# Author: Pedro Arduino
# UW Computational Geotechnical Group
# Date: 2024
# All Rights Reserved
# ############################################################################################################
# %%
import numpy as np
import json
import matplotlib.pyplot as plt
from numpy.fft import fft, ifft
from scipy import integrate
from respSpectra import resp_spectra
class TFunctionClass:
def __init__(self, damping, H, Vs):
# Define the variables
self.m_freq = None
self.m_time = None
self.m_acc = None
self.m_absFft = None
self.m_absSoilTF = None
self.m_absIFft = None
self.m_accT = None
self.m_vel = None
self.m_disp = None
self.m_velT = None
self.m_dispT = None
# Define soil layer parameters
self.m_damping = damping / 100.0 # damping from percentage to number
self.m_H = H
self.m_Vs = Vs
def calculateResponse(self):
SoilTF = np.empty_like(self.m_freq, dtype=np.complex_)
absSoilTF = np.empty_like(self.m_freq, dtype=np.float_)
# Compute the Fourier amplitude spectrum
fas = fft(self.m_acc)
# fas = fas[:self.nyquist_index]
absfas = np.abs(fas)
self.m_absFft = absfas
# Compute transfer function of soil layer
SoilTF = self.calcSoilTf()
self.m_absSoilTF = np.abs(SoilTF)
# Compute surface soil response
ifas = fas * SoilTF
absfas2 = np.abs(ifas)
self.m_absIFft = absfas2
accT = ifft(ifas)
self.m_accT = accT.real # Take only the real part
def calcSoilTf(self):
tf = []
if self.m_freq is None:
print("Frequency vector is not defined")
else:
for f in self.m_freq:
"""
* The uniform damped soil on rigid rock transfer function
* 1
* H = -------------------------------------------------
* cos ( 2* PI * freq * H / (Vs(1+ i*damping))
"""
kstar = 2.0 * np.pi * f / self.m_Vs - self.m_damping * 2.0 * np.pi * f / self.m_Vs * 1j
Vsstar = self.m_Vs + self.m_damping * self.m_Vs * 1j
tf.append(1.0 / np.cos(2.0 * np.pi * f * self.m_H / Vsstar))
return tf
def calculate_nat_freq(self):
n_pt = len(self.m_freq)
N_freq = []
N_freqVal = []
dfreq = self.m_freq[-1] / n_pt
TF_tan = 1.0
for i in range(1, len(self.m_freq)):
TF_tan1 = (self.m_absSoilTF[i] - self.m_absSoilTF[i - 1]) / dfreq
if TF_tan1 * TF_tan <= 0 and TF_tan > 0:
N_freq.append(self.m_freq[i])
N_freqVal.append(self.m_absSoilTF[i])
TF_tan = TF_tan1
return N_freq, N_freqVal
def calculate_ratio(self):
grav = 9.81 # m/s2
dT = self.m_time[1] - self.m_time[0]
accAux = [self.m_acc[ii]*grav for ii in range(len(self.m_acc))]
self.m_vel = integrate.cumtrapz(accAux, dx=dT)
# self.mvel = np.insert(self.m_vel, 0, 0.0)
self.m_disp = integrate.cumtrapz(self.m_vel, dx=dT)
# mdisp = np.insert(self.m_disp, 0, 0.0)
self.m_velT = integrate.cumtrapz((self.m_accT*grav), dx=dT)
# self.mvel = np.insert(self.m_vel, 0, 0.0)
self.m_dispT = integrate.cumtrapz(self.m_velT, dx=dT)
# mdisp = np.insert(self.m_disp, 0, 0.0)
ratioA = abs(max(self.m_accT))/abs(max(self.m_acc))
ratioV = abs(max(self.m_velT))/abs(max(self.m_vel))
return ratioA, ratioV
def sin_record(self, f):
n_points = 2000
self.m_dt = 0.02
self.m_acc = [0] * n_points
accel = []
for s in range(n_points):
accel.append(0.4 * np.sin(2 * f * np.pi * s * self.m_dt))
self.m_acc = accel
self.set_time()
self.set_freq()
def sweep_record(self):
n_points = 8000
self.m_dt = 0.002
self.m_acc = [0] * n_points
self.m_time = [0] * n_points
for i in range(len(self.m_time)):
time = i * self.m_dt
self.m_time[i] = time
self.m_acc[i] = np.sin(25.0 * time + 150.0 * (time * time / 2.0) / 16.0)
self.set_freq()
def load_file(self, file_name):
self.m_filename = file_name
try:
with open(file_name, 'r') as file:
# Read file contents into a JSON object
jsonObj = json.load(file)
except FileNotFoundError as e:
print(f"Cannot read file {file_name}: {e}")
return
events = jsonObj.get("Events", [])
if events:
patterns = events[0].get("pattern", [])
timeseries = events[0].get("timeSeries", [])
pattern_type = patterns[0].get("type", "")
tsname = patterns[0].get("timeSeries", "")
units = events[0].get("units", {})
acc_unit = 1.0
acc_type = units.get("acc", "")
if acc_type == "g":
acc_unit = 1.0
elif acc_type == "m/s2":
acc_unit = 1.0 / 9.81
elif acc_type in ["cm/s2", "gal", "Gal"]:
acc_unit = 1.0 / 981.0
timeseries_data = timeseries[0].get("data", [])
dT = timeseries[0].get("dT", 0.0)
self.read_GM(timeseries_data, dT, acc_unit)
def read_GM(self, acc_TH, dT, acc_unit):
n_points = len(acc_TH)
self.m_dt = dT
# self.m_acc = [acc_TH[ii].toDouble() * acc_unit for ii in range(n_points)]
self.m_acc = [acc_TH[ii] * acc_unit for ii in range(n_points)]
if n_points % 2 == 0:
self.m_acc.append(0.0)
self.m_acc = np.array(self.m_acc) # Convert to numpy array
self.set_time()
self.set_freq()
def set_freq(self):
if self.m_dt == 0:
self.m_dt = 0.005
nfreq = 1 / self.m_dt*10
sample_freq = 1.0 / self.m_dt
else:
nfreq = len(self.m_acc)
sample_freq = 1.0 / self.m_dt
# self.m_freq = [0] * (len(self.m_acc) // 2 + 1)
# self.m_freq = [0] * (len(self.m_acc)) # m_freq as a list
self.m_freq = np.zeros(nfreq) # m_freq as a numpy array
sample_freq = 1.0 / self.m_dt
self.nyquist_freq = sample_freq / 2.0
self.nyquist_index = int(len(self.m_freq) / 2)
for i in range(len(self.m_freq)):
self.m_freq[i] = i * sample_freq / len(self.m_acc)
def set_time(self):
# self.m_time = [0] * len(self.m_acc) # m_time as a list
self.m_time = np.zeros(len(self.m_acc)) # m_time as a numpy array
for i in range(len(self.m_time)):
self.m_time[i] = i * self.m_dt
def plot_acc(self):
plt.figure()
plt.plot(self.m_time, self.m_acc, 'b-', label='input')
plt.plot(self.m_time, self.m_accT, 'r-', label='output')
plt.xlabel('Time [sec]')
plt.ylabel('Acc [g]')
plt.legend()
plt.show()
def plot_fft(self):
plt.figure()
plt.plot(self.m_freq[:self.nyquist_index], self.m_absFft[:self.nyquist_index], 'b-', label='input')
plt.plot(self.m_freq[:self.nyquist_index], self.m_absIFft[:self.nyquist_index], 'r-', label='output')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Fourier Amplitude')
plt.legend()
plt.show()
def plot_tf(self):
plt.figure()
plt.plot(self.m_freq[:self.nyquist_index], self.m_absSoilTF[:self.nyquist_index], 'b-')
plt.xlabel('Frequency [Hz]')
plt.ylabel('TF')
plt.show()
def plot_spectra(self):
n_points = len(self.m_acc)
accAux = [self.m_acc[ii]*9.81 for ii in range(n_points)]
accTAux = [self.m_accT[ii]*9.81 for ii in range(n_points)]
periods, psa = resp_spectra(self.m_time, accAux, 0.05)
periodsT, psaT = resp_spectra(self.m_time, accTAux, 0.05)
plt.figure()
plt.plot(periods, psa, 'b-', label='input')
plt.plot(periodsT, psaT, 'r-', label='output')
plt.xlabel('Periods [s]')
plt.ylabel('PSA [cm/s2]')
plt.legend()
plt.show()
def main():
# Define input parameters
damping = 5.0 # damping ratio in %
H = 20.0 # layer height in m
Vs = 200.0 # shear wave velocity in m/s
TF = TFunctionClass(damping, H, Vs)
# Sinusoidal record
f = 0.5 # frequency in Hz
TF.sin_record(f)
# Calculate response
TF.calculateResponse()
# Calculate ratios
ratioA, ratioV = TF.calculate_ratio()
print(f"Acceleration Ratio: {ratioA}")
print(f"Velocity Ratio: {ratioV}")
# Plot acceleration
TF.plot_acc()
# Plot Fourier Transform
TF.plot_fft()
# Plot Transfer Function
TF.plot_tf()
# Plot Spectra
TF.plot_spectra()
if __name__ == "__main__":
main()
.. raw:: html
.. raw:: html
Click to expand the full Response Spectra Python code
.. code-block:: python
#########################################################
#
# Postprocessing python script
#
# Copyright: UW Computational Mechanics Group
# Pedro Arduino
#
# Participants: Alborz Ghofrani
# Long Chen
#
#-------------------------------------------------------
import numpy as np
def resp_spectra(a, time, nstep):
"""
This function builds response spectra from acceleration time history,
a should be a numpy array,T and nStep should be integers.
"""
# add initial zero value to acceleration and change units
a = np.insert(a, 0, 0)
# number of periods at which spectral values are to be computed
nperiod = 100
# define range of considered periods by power of 10
minpower = -3.0
maxpower = 1.0
# create vector of considered periods
p = np.logspace(minpower, maxpower, nperiod)
# incremental circular frequency
dw = 2.0 * np.pi / time
# vector of circular freq
w = np.arange(0, (nstep+1)*dw, dw)
# fast fourier Horm of acceleration
afft = np.fft.fft(a)
# arbitrary stiffness value
k = 1000.0
# damping ratio
damp = 0.05
umax = np.zeros(nperiod)
vmax = np.zeros(nperiod)
amax = np.zeros(nperiod)
# loop to compute spectral values at each period
for j in range(0, nperiod):
# compute mass and dashpot coeff to produce desired periods
m = ((p[j]/(2*np.pi))**2)*k
c = 2*damp*(k*m)**0.5
h = np.zeros(nstep+2, dtype=complex)
# compute transfer function
for l in range(0, int(nstep/2+1)):
h[l] = 1./(-m*w[l]*w[l] + 1j*c*w[l] + k)
# mirror image of Her function
h[nstep+1-l] = np.conj(h[l])
# compute displacement in frequency domain using Her function
qfft = -m*afft
u = np.zeros(nstep+1, dtype=complex)
for l in range(0, nstep+1):
u[l] = h[l]*qfft[l]
# compute displacement in time domain (ignore imaginary part)
utime = np.real(np.fft.ifft(u))
# spectral displacement, velocity, and acceleration
umax[j] = np.max(np.abs(utime))
vmax[j] = (2*np.pi/p[j])*umax[j]
amax[j] = (2*np.pi/p[j])*vmax[j]
return p, umax, vmax, amax
.. raw:: html
.. raw:: html
.. table:: Table 1. Natural Frequencies in the Transfer Function
:widths: auto
+------------+---------------------------------------------+
| Peak | Amplification Factor at Natural Frequencies |
+============+=============================================+
| 1 | 20.49 |
+------------+---------------------------------------------+
| 2 | 7.03 |
+------------+---------------------------------------------+
| 3 | 4.20 |
+------------+---------------------------------------------+
| 4 | 2.98 |
+------------+---------------------------------------------+
.. raw:: html
With the transfer function plotted, the input motion can be transformed using the transfer function to reflect the motion at the surface. The figure below reflects the large amplification that occurred. The value of the highest acceleration increased from ~0.4g in the rock to ~1.25g in the soil. This amplification is also reflected in other frequencies.
.. figure:: ./images/case2_Full_Results_TF.png
:scale: 40 %
:align: center
Fig. 17. Amplification of Ground Motion.
The spectral acceleration spectra can be also determined for each of the motions. These spectra can be used to determine if a structure will be affected by the amplification. A single story structure (~0.1 second period) might be at risk due to this amplification. Any periods with a large amplification ratios should be further analyzed to ensure the safety of the structure and site.
.. figure:: ./images/case2_SpectralAcc_Results_TF.png
:scale: 40 %
:align: center
Fig. 18. Amplification in Spectral Acceleration.
Remarks
-------
I'd like to thank everyone at SimCenter, specifically Sang-ri Yi, Frank McKenna, Jinyan Zhao, Aakash Bangalore Satish, and Barbaros Cetiner, for all of their effort and assistance they provided during the entire quarter. Navigating these tools and creating examples for them would've been a lot more stressful without their help.
Transfer function is one of my favorite topics in geotechnical engineering. I'd really like to continue working with site response and performance based design so being able to create this example along with my class was great.
Finally, I'd like to thank Prof. Arduino who made all of this possible. His determination and motivation was contagious throughout the academic quarter. There are many great professors but there is only one Pedro Arduino.