Earthquake as a triggering factor of an infinite slope mechanism

Earthquake as a triggering factor of an infinite slope mechanism#

© 2024 Daniel F. Ruiz, Exneyder A. Montoya-Araque y Universidad EAFIT.

This notebook can be interactively run in Google - Colab.

This notebook runs the model pyNermarkDisp developed by Montoya-Araque et al. [2024] based on the classical sliding rigid block method by Newmark [1965].

Required modules and global setup for plots#

import os
import sys
import subprocess

if 'google.colab' in str(get_ipython()):
    print('Running on CoLab. Installing the required modules...')
    # subprocess.run('pip install ipympl', shell=True);
    subprocess.run('pip install pynewmarkdisp', shell=True);
    from google.colab import output, files
    output.enable_custom_widget_manager()
else:
    import tkinter as tk
    from tkinter.filedialog import askopenfilename

import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
from pynewmarkdisp.newmark import direct_newmark, plot_newmark_integration
from pynewmarkdisp.infslope import factor_of_safety, get_ky

# %matplotlib widget
%matplotlib inline
mpl.rcParams.update({
    "font.family": "serif",
    "font.serif": ["Computer Modern Roman", "cmr", "cmr10", "DejaVu Serif"],  # or 
    "mathtext.fontset": "cm",  # Use Computer Modern fonts for math
    "axes.formatter.use_mathtext": True,  # Use mathtext for axis labels
    "axes.unicode_minus": False,   # Use standard minus sign instead of a unicode character
})

# Create a folder where the records will be saved
workdir = os.getcwd()
records_dir = os.path.join(workdir, "records")
os.makedirs(f"{records_dir}", exist_ok=True)

Functions#

def load_file(testing_data):
    if testing_data:
        url = "https://raw.githubusercontent.com/eamontoyaa/data4testing/main/pynewmarkdisp/"
        # record = np.loadtxt(f"{url}earthquake_data_friuli.csv", delimiter=";", skiprows=1)
        # record[:, 1] = - 2.0 * 9.81 * record[:, 1]  # Convert to SI units
        record = np.loadtxt(f"{url}earthquake_data_armenia1999.csv", delimiter=",")
    elif testing_data is False and 'google.colab' in str(get_ipython()):
        file = files.upload()
        # Move the file to the records folder
        os.rename(list(file.keys())[0], os.path.join(records_dir, list(file.keys())[0]))
        record = np.loadtxt(os.path.join(records_dir, list(file.keys())[0]), delimiter=",")
    else:  # GUI for file selection from local machine if not in CoLab
        tk.Tk().withdraw() # part of the import if you are not using other tkinter functions
        file = askopenfilename()
        # Move the file to the records folder
        os.rename(file, os.path.join(records_dir, os.path.basename(file)))
        record = np.loadtxt(os.path.join(records_dir, os.path.basename(file)))
    time, accel = record.T
    return time, accel
    
def get_arias_intensity(time, accel, g):
    return np.trapz(y=accel**2, x=time) * np.pi / (2 * g)

def get_pgv(time, accel, g):
    accel = 9.81 * accel / g  # Earthquake acceleration to SI units
    vel = np.cumsum(accel) * (time[1] - time[0])
    return np.max(np.abs(vel))

def get_pga(accel, g):
    return np.max(np.abs(accel / g))

Running a case#

testing_data = True  # Set to False to use the GUI to load the data from an external file

Inputs (conditioning factors)#

# Geometry
β = 35  # [°] - Slope angle - Plane inclination
d = 3.0  # [m] - Depth of the slip surface - Block height
d_w = 3.0  # [m] - Depth of the water table

# Material parameters
φ = 27  # [°] - Friction angle
c = 15  # [kPa] - Cohesion
γ = 19  # [kN/m³] - Unit weight

Inputs (triggering factor - Earthquake record)#

time, accel = load_file(testing_data)
time, accel

(array([0.0000e+00, 5.0000e-03, 1.0000e-02, ..., 3.4985e+01, 3.4990e+01,
        3.4995e+01]),
 array([0.00246051, 0.00246493, 0.00246935, ..., 0.02309875, 0.05239547,
        0.08169219]))

# Value of gravity acceleration in the same units as in the earthquake record
# g = 1.0   # → It means that acceleration units in the input file are given in [g]
g = 9.81  #  → It means that acceleration units in the input file are given in [m/s²]
# g = 981  #  → It means that acceleration units in the input file are given in [cm/s²]

PGA = get_pga(accel, g)
print(f"PGA: {PGA:.2f} g")

PGA: 0.54 g

Static factor of safety, \(\mathrm{FS}_\mathrm{static}\)#

fs = factor_of_safety(d, d_w, β, φ, c, γ, ks=0)
print(f"Static factor of safety: {fs:.2f}")

Static factor of safety: 1.29

Critical seismic coefficient \(k_\mathrm{y}\)#

ky = get_ky(d, d_w, β, φ, c, γ)
print(f"Critical seismic coefficient, ky: {ky:.2f}")

Critical seismic coefficient, ky: 0.15

Pseudo-static factor of safety, \(\mathrm{FS}_\mathrm{pseudostatic}\), when \(k_\mathrm{s}\) is 70% of \(k_\mathrm{y}\)#

ks = 0.7 * ky
fs_ks = factor_of_safety(d, d_w, β, φ, c, γ, ks)
print(f"Seismic coefficient, ks: {ks:.3f}")
print(f"Pseudostatic factor of safety: {fs_ks:.2f}")

decrease_fs = (fs - fs_ks) / fs * 100
print(f"Decrease in factor of safety: {decrease_fs:.2f} %")

Seismic coefficient, ks: 0.104
Pseudostatic factor of safety: 1.07
Decrease in factor of safety: 16.54 %

Pseudo-static factor of safety, \(\mathrm{FS}_\mathrm{pseudostatic}\), when \(k_\mathrm{s}\) is 50% of PGA#

ks = 0.5 * PGA
fs_ks = factor_of_safety(d, d_w, β, φ, c, γ, ks)
print(f"Seismic coefficient, ks: {ks:.3f}")
print(f"Pseudostatic factor of safety: {fs_ks:.2f}")

decrease_fs = (fs - fs_ks) / fs * 100
print(f"Decrease in factor of safety: {decrease_fs:.2f} %")

Seismic coefficient, ks: 0.270
Pseudostatic factor of safety: 0.83
Decrease in factor of safety: 35.56 %

Calculating and plotting \(u_\mathrm{p}\)#

permanent_disp = direct_newmark(time, accel, ky, g)
# permanent_disp = direct_newmark(time, accel, .135, g)
fig = plot_newmark_integration(permanent_disp)
fig.canvas.header_visible = False 
fig.canvas.toolbar_position = 'bottom' 
plt.show()
../../_images/37e7f1f0de4820d2f3fb6883c471af2aa1cf7a0d35ab6b776a8edb3da01c7aaa.png
Contents