+ 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: 50 %
+ :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.
diff --git a/source/images/case2_Gaussian_Mixture_RatioV_histogram.png b/source/images/case2_Gaussian_Mixture_RatioV_histogram.png
index 4beaab7..a84a797 100644
Binary files a/source/images/case2_Gaussian_Mixture_RatioV_histogram.png and b/source/images/case2_Gaussian_Mixture_RatioV_histogram.png differ
diff --git a/source/images/case2_Gaussian_Mixture_RatioV_histogramOLD.png b/source/images/case2_Gaussian_Mixture_RatioV_histogramOLD.png
new file mode 100644
index 0000000..4beaab7
Binary files /dev/null and b/source/images/case2_Gaussian_Mixture_RatioV_histogramOLD.png differ