Month / January 2022

by Agus Daud

Studi kasus: Python untuk Geoteknik

Beberapa bulan terakhir involve dalam proyek yang mempunyai data pengujian geoteknik yang sangat banyak. Baik itu data pengujian in-situ (bor SPT, sondir..) dan data laboratorium (triaxial, analisa saringan, direct shear..)

Pekerjaan pengujian tanah tersebut dilakukan oleh kontraktor geoteknik, kemudian diserahkan ke konsultan untuk analisa lebih lanjut.

Sayangnya, data yang diberikan ada kemiripan/sama satu sama lain. Kontraktor dicurigai melakukan kecurangan.

Singkat kata, Bos memerintahkan saya untuk mengecek data laboratorium tersebut. Data apa saja yang sama/mirip dan teridentifikasi kecurangan.

Sialnya, data yang diberikan hanyalah file pdf. Dan dalam satu file pdf terdapat beberapa format tabel seperti analisa saringan, analisa saringan+hidrometer, geser langsung, triaxial test, berat jenis, dan beberapa uji atterberg.

Sekitar 1020an halaman khusus data lab (total 2000an) yang akan dicek, dan harus diselesaikan kurang dari 2 hari.

Ini bisa saja kita lakukan secara manual, melihat dan membandingkan. Namun, saya tidak begitu yakin apakah bisa diselesaikan dengan waktu kurang dari 2 hari dengan hasil akurat.

Bekerja dengan data yang banyak dan berulang. Disinilah bahasa pemrograman python menjadi senjata ampuh.

Hal yang dilakukan adalah:

  1. Meng-extract file pdf per page ke dalam file text
  2. Mengidentifikasi page termasuk dalam format data apa? (triaxial, direct shear, sieve analysis.. etc)
  3. Memilah tulisan (angka) yang termasuk dalam data, memasukkan angka tersebut dalam dataframe pandas dan meng-export ke file csv untuk diolah.
  4. Membanding data sesuai jenis pengujian/test, dengan cara menghitung selisih setiap data. Jika selisih = nol (0) maka data dikatakan sama persis
  5. Selesai

Tugas ini diselesaikan dalam beberapa jam saja. Karena membutuhkan waktu untuk ngulik google, dan perlu debugging codingan error. Selain itu, ternyata proses ekstrasi text pdf per halaman membutuhkan waktu.

Paling banyak error karena ternyata dalam satu format tabel juga masih beda. Data uji analisa saringan misalnya, ukuran saringan yang digunakan tidak seragam sehingga jumlah baris atau row data berbeda-beda.

Ini hanyalah contoh kasus sederhana, dimana kita dapat menggunakan python untuk keperluan sehari-hari. Tentunya jika diperlukan analisa lebih lanjut akan lebih mudah. Misal koreksi N-SPT untuk setiap titik pengujian dan per 1,5m kedalaman. Koreksi bisa dilakukan dengan cepat dan akurat (tentunya dengan verifikasi).

Demikian semoga bermanfaat.

Advertisement
by Agus Daud

Anastruct: Solusi Beban Terdistribusi Tidak Seragam

Dalam analisa struktur, biasanya ditemukan beban merata yang tidak seragam (non-uniform). Seperti beban akibat tekanan tanah lateral berbentuk segitiga, trapesium atau sembarang bentuk.

Tekanan tanah lateral (Hardiyatmo, C. 2011)

Menjadi masalah ketika program tidak menyediakan opsi input beban tidak seragam, seperti package Anastruct hanya menyediakan beban merata (q) uniform atau seragam.

Artikel ini terinspirasi dari blog Syont di sini : Beban titik kesetaran pada elemen balok

Secara ringkas, beliau menggunakan beban titik sebagai solusi analisa beban tidak seragam.

***

Baiklah..

Artikel kali ini akan membahas solusi analisa struktur beban tidak seragam menggunakan Anastruct/ Python.

Strategi Pemodelan

Strategi yang digunakan yaitu dengan membagi elemen dalam elemen pias/diskrit. Kemudian menggunakan beban q rata-rata (q tengah elemen) sebagai beban merata.

Berikut kira-kira ilustrasinya,

Strategi pemodelan beban tidak seragam

1. Membuat fungsi beban

Mulai dari sini pemodelan akan menggunakan Python dan Jupyter Notebook. Langkah pertama adalah membuat fungsi beban. Fungsi ini akan meng-handle beban merata setiap variasi kedalaman (z) atau jarak.

Lateral pressure due to line load Q:
\sigma_h = \frac{4Q}{\pi} \frac{x^2z}{R^4} 
Lateral active earth pressure Rankine theory:
K_a = \tan (45 - \phi/2)^2 
\sigma_a = \gamma z K_a

Di sini disiapkan dua fungsi. Yaitu fungsi tekanan tanah Rankine dan tekanan lateral akibat beban garis di atas dinding. Berikut script nya:

import math
def Q_pressure(H, q, x, z):
    R = math.sqrt(x**2 + z**2)
    sigma_q = 4*q/math.pi * (x**2)*z / (R**4)
    return(sigma_q)
def Rankine(phi, gamma, z):
    Ka = math.tan(math.radians(45-(phi/2)))**2
    sigma_a = gamma*z*Ka
    return(sigma_a)
Persamaan/fungsi tekanan tanah

2. Menentukan data dimensi, material, dan pengaturan yang digunakan

Sebagai verifikasi awal, maka akan ditinjau

  • Dinding kantilever
  • Tinggi total 6 m
  • Beban = tekanan tanah aktif teori Rankine
  • Jumlah elemen diskrit = 6 elemen

Sehingga panjang setiap elemen menjadi 6/6 = 1 m. Berikut input data keseluruhan:

import numpy as np
import pandas as pd
df = pd.DataFrame()
# Dimensions
H = 6 # wall height (m)
t = H/10 # thickness of the wall (m)
b = 1 # width (m)
A = b*t # area section (m^2)
Ix = (b*t**3)/12 # Momen inertia (m4)
# Materials
## Concrete
fc = 21 # MPa
E = 4700*math.sqrt(fc) * 1000 # concrete modulus elasticity (kPa)
## Soil
gamma = 18 # unit weight of soil
phi = 30 # internal friction angle
# Loads
Q = 112.5 # acting load (kN/m)
x = 2.0 # horizontal distance from wall (m)
# Settings
n_ele = 6 # number of element
dx = H/n_ele # lenght of dicretized element (m)

3. Membuat tabel beban tekanan tanah untuk setiap kedalaman (z)

Pada tahap ini, terlebih dahulu akan dicoba tekanan tanah rankine (beban segitiga).

z_list = [0]
sigma_list = [0]
for n in range(1,n_ele+1,1):
    z = n*dx
    z_list.append(-z)
    
    sigma_list.append(Rankine(phi,gamma,z)) # Rankine only
    #sigma_list.append(Q_pressure(H,Q,x,z)) # Lateral pressure due to line load only
    #sigma_list.append(Q_pressure(H,Q,x,z) + Rankine(phi,gamma,z)) # Rankine + Q lateral pressure
    
df['z (m)'] = z_list
df['sigma (kPa)'] = sigma_list 
df
Tabel tekanan tanah dibuat dalam data frame pandas
Plot data tekanan tanah

4. Analisa struktur dengan Anastruct

Berdasarkan data di atas, kemudian lakukan analisa struktur dengan Anastruct. Perlu diperhatikan, input beban ke dalam elemen dilakukan pada perintah loop “for i in range ()“; dengan terlebih dahulu dilakukan interpolasi pada kedalaman (z) tengah elemen. Ditunjukkan pada garis <<< assign q load each element.

Perlu perhatian khusus pada tumpuan, disesuaikan dengan dinding tinjauan yaitu fixed pada dasar dinding.

from scipy.interpolate import interp1d
from anastruct import SystemElements
ss = SystemElements(EI=E*Ix, EA=E*A)
ss.add_multiple_elements([[0, H], [0, 0]], n_ele)
##------------Support conditions very important!----------##
# first_node == node on the top of the wall
first_node = 1
last_node = n_ele+1
# ss.add_support_roll(first_node,direction='y')
# ss.add_support_hinged(last_node)
ss.add_support_fixed(last_node)
first_ele_z = dx/2
for i in range(n_ele):
    z = (i*dx) + first_ele_z
    interp = interp1d(df['z (m)'], df['sigma (kPa)']) # <<< assign q load each element
    sigma_interp = interp(-z)
    ss.q_load(q=-sigma_interp, element_id = i+1)
# Run Analysis
ss.solve()
# Plot graph
ss.show_structure(verbosity=0)
ss.show_bending_moment(verbosity=1)
ss.show_shear_force(verbosity=1)
# verbosity=0 to show value

Berikut tampilan beban yang telah di input:

Input beban
Output reaksi perletakan

Untuk verifikasi hasil reaksi tumpuan di atas dapat dilakukan dengan cara sederhana sebagai berikut:

Ka = tan(45 – 30/2)2 = 0.333

Pa = 0.5 x gamma x H2 x Ka = 107,89 kN ≈ Ranastruct = 108.0 kN (mendekati)

Mo = Pa x H/3 = 215,78 kNm ≈ Manastruct = 219.0 kNm (cukup mendekati nilai teori)

Terlihat reaksi hasil pemodelan cukup mendekati, perlu di perhatikan ini adalah hasil untuk 6 element diskrit. Nilai akan semakin mendekati jika jumlah elemen di perbanyak; misal 100 elemen.

5. Tahapan post processing

Pada tahap ini, gaya dan perpindahan untuk setiap node akan dirangkum dalam table pandas data frame. Kemudian di plot dalam grafik agar lebih mudah dilihat.

moment_list = ss.get_element_result_range('moment')
shear_list = ss.get_element_result_range('shear')
axial_list = ss.get_element_result_range('axial')
for n in range(0,n_ele+1):
    if n == n_ele:
        df.loc[n,'Moment (kNm)'] = ss.get_node_results_system(node_id=n+1)['Ty']
        df.loc[n,'Shear (kN)'] = -ss.get_node_results_system(node_id=n+1)['Fx']
        df.loc[n,'Axial (kNm)'] = ss.get_node_results_system(node_id=n+1)['Fy']
    else:
        df.loc[n,'Moment (kNm)'] = moment_list[n]
        df.loc[n,'Shear (kN)'] = shear_list[n]
        df.loc[n,'Axial (kNm)'] = axial_list[n]
    df.loc[n,'Displacement (mm)'] = 1000 * ss.get_node_results_system(node_id=n+1)['ux']
    
# Display table
df
Tabel output gaya dan perpindahan

Kemudian di plot menggunakan package plotly:

from plotly.subplots import make_subplots
import plotly.graph_objects as go
fig = make_subplots(rows=1, cols=4)
fig.add_trace(go.Scatter(name="Lateral Pressure (kPa)",x=df['sigma (kPa)'], y=df['z (m)'], mode="lines"), row=1, col=1)
fig.add_trace(go.Scatter(name="Bending Moment (kNm)",x=df['Moment (kNm)'], y=df['z (m)'], mode="lines"), row=1, col=2)
fig.add_trace(go.Scatter(name="Shear (kN)",x=df['Shear (kN)'], y=df['z (m)'], mode="lines"), row=1, col=3)
fig.add_trace(go.Scatter(name="Displacement (mm)", x=df['Displacement (mm)'], y=df['z (m)'], mode="lines"), row=1, col=4)
fig.update_yaxes(title_text="Depth, z (m)", row=1, col=1)
fig.update_xaxes(title_text="Lateral Pressure (kPa)", row=1, col=1)
fig.update_xaxes(title_text="Moment (kNm)", row=1, col=2)
fig.update_xaxes(title_text="Shear (kN)", row=1, col=3)
fig.update_xaxes(title_text="Displacement (mm)", row=1, col=4)
fig.update_layout(height=800, width=1000, title_text="Result",)
fig.show()
Plot kedalaman (z) vs gaya/ perpindahan

***

Pemodelan beban distribusi untuk sembarang bentuk

Eksperimen terakhir, akan dicoba pemodelan beban untuk sembarang bentuk sesuai dua fungsi yang telah dibuat (seperti gambar tekanan tanah pada awal artikel). Berikut deskripsi soalnya:

Data:

  • Tinggi dinding, H = 6 m
  • Kondisi tumpuan = jepit – jepit (fixed)
  • Jumlah element diskrit = 100
  • Data material sama dengan soal di atas

Tinjauan kondisi beban:

  1. Beban tekanan tanah aktif Rankine
  2. Beban tekanan lateral akibat beban garis (Q) di atas dinding
  3. Kombinasi beban 1 dan 2

Berikut hasilnya..

Kondisi 1. Tekanan tanah aktif Rankine
Kondisi 2. Tekanan tanah lateral akibat beban garis Q di atas dinding
Kondisi 3. Kombinasi tekanan tanah aktif Rankine dan beban garis Q di atas dinding

Bagi yang ingin mencoba silahkan download file jupyter notebook berikut:

Lateral Pressure Due to Surcharge

Bila ada yang ingin didiskusikan, silahkan tinggalkan komentar di bawah ini. Demikian semoga bermanfaat.

Catatan:

  • Untuk variasi beban lain dapat dimasukkan dalam tambahan fungsi.
by Agus Daud

LEMSlope: Pseudostatik analisis

Artikel kali ini mau nyoba analisis faktor aman lereng terhadap gempa, menggunakan LEMSlope dan membandingkan hasilnya dengan Geostudio Slope/W.

Metode yang digunakan adalah pseudostatik analisis, dimana beban gempa diberikan dengan cara input koefisien gempa horizontal (kh).

Tujuan utama analisis kali ini adalah untuk mencari koefisien gempa horizontal yang menghasilkan faktor aman lereng FS = 1,0. Atau dikenal dengan k yield (ky).

Nilai ky biasanya digunakan untuk menghitung deformasi permanen akibat gempa; namun, perhitungan deformasi tidak akan dibahas kali ini.

Baiklah..

Berikut lereng yang akan dianalisis:

Lereng Tinjauan

Untuk mendapatkan nilai ky, maka akan dicoba beberapa nilai kh (0 – 0.30g) kemudian akan diplot grafik hubungan antara kh dan faktor aman (FS).

Berdasarkan soal lereng di atas, maka ditulis script LEMSlope sebagai berikut:

#~units: m|kN|kPa|kN/m3|degree

import lem

#~materials
sand = lem.new_material_mohr_coulomb(
	name = 'sand',
	color = '#FFFBA9', 
	unit_weight = 20,
	c = 5, 
	phi = 30)

#~base
base = lem.create_base()
base.set_limits(
	left = 0, 
	right = 40, 
	bottom = 0)
base.add_hline(20)

base.assign_material(sand, (1,1))

#~analysis cases
case1 = lem.new_case('Case 1')
case1.cut_right_slope(
	bottom_elevation = 10,
	slope = '2:1',
	x = 30,
	is_top_x = False)

# deskripsi bidang gelincir entry-exit
	
case1.set_entry_exit_limits(
	left_limits = (0,10), 
	left_points_number = 20, 
	right_limits = (30, 40), 
	right_points_number = 20,
	radius_increment_number = 5, 
	slice_number = 30, 
	left_y_limits = [], 
	right_y_limits = [])
	
# input koefisien gempa kh
case1.apply_seismic_load(
	horizontal_coefficient = 0.2)

Sebagai contoh untuk koefisien gempa kh = 0,2; diperoleh faktor aman lereng FS = 1,077.

kh = 0,20 menghasilkan FS = 1,077.

Adapun nilai faktor aman untuk masing-masing nilai kh, diplot sebagai berikut.

Plot nilai kh vs FS; dipeoleh ky = 0,25g

Dari grafik di atas terlihat faktor aman hasil LEMSlope lebih rendah (merah) dibandingkan dengan hasil Geostudio Slope/W (biru). Namun, masih mendekati.

Untuk kondisi lereng kereng sesuai soal di atas; dari hasil plot nilai FS = 1,0 diperoleh nilai ky = 0,25 g.

Dengan demikian LEMSlope dianggap sudah memenuhi analisa faktor aman lereng akibat gempa; metode pseudostatik analisis.

Demikian semoga bermanfaat.