Category / Openseespy

by Agus Daud

OpenSees mulai dari yang sederhana

Mempelajari OpenSees memang susah-susah gampang. Bagi user yang belum terbiasa dengan finite element dan bahasa memrograman, pasti akan sedikit tersiksa. Itu pula kesan yang saya rasakan saat menggunakan OpenSees.

Terlebih lagi tutorial yang ada banyak memperlihatkan script pemrograman yang rumit dan menjadi momok bagi user.

Prof. Scott mengatakan bahwa OpenSees dan nonlinear analisis struktur secara umum adalah sederhana, tidak kompleks dan tidak perlu kompleks.

Walaupun memang tujuan OpenSees dibuat adalah untuk simulasi performa struktur dan geoteknik terkait gempa, yang terlihat tidak sederhana.

Mulai dari yang sederhana

Seperti biasa, untuk mempelajari aplikasi baru biasanya saya lebih suka memulai dengan model yang sederhana.

Menurut Minjie Zhu di sini, secara umum script OpenSees dapat dibagi menjadi 2 yaitu bagian domain dan analysis.

  1. Domain terkait pembuatan geometri seperti node, element, kondisi tumpuan dsb.
  2. Analysis berisi pengaturan algorithm, system, numberer, constraint, integrators, dsb.

Sebagai tambahan, ada bagian ke 3 yaitu output dan postprocessing. Terkait merekam (record), menampilkan hasil (output) dan plot dalam grafik.

Baiklah.. kali ini akan dicoba cek lendutan yang terjadi pada balok sederhana dengan data sebagai berikut.

Menggunakan rumus praktis,

lendutan = PL^3 / 48EI = 1,92 mm

Berikut deskripsi tahapan pemodelan balok sederhana dengan OpenSeesPy.

(1) Pertama, tahapan awal adalah meng-import modul OpenSeesPy, dan menulis data penampang, material, serta geometri struktur.

# memasukkan modul opensees
from openseespy.opensees import *
# membersihkan model yang terdahulu
wipe()
# Unit mm, N, MPa
P = 10000
L = 6000
b = 250
h = 500
A = b*h
Iz = (b*h**3)/12
E = 9000

(2) Kedua, mengatur tinjauan dimensi dan derajat kebebasan. pada contoh ini struktur 2D, dengan derajat kebebasan 3 untuk balok.

# model dimensi (-ndm) dan derajat kebebasan (-ndf)
model('basic','-ndm', 2,'-ndf', 3)

(3) Ketiga, memasukkan node dan kondisi tumpuan

# membuat geometri dan tumpuan
# format node(id, x, y)
node(1, 0.0, 0.0) ; fix(1,1,1,0)
node(2, L/2, 0.0)
node(3, L, 0.0)   ; fix(3,0,1,0)

(4) Keempat, memasukkan jenis transformasi dan elemen balok. Bentang balok dibagi menjadi dua elemen seperti gambar soal.

# transformasi
geomTransf('Linear',1)
# memasukkan elemen
element('elasticBeamColumn', 1, 1,2, A,E,Iz,1)
element('elasticBeamColumn', 2, 2,3, A,E,Iz,1)

(5) Kelima, masukkan time series, load pattern, dan beban. Digunakan timeSeries Linear, karena contoh ini menggunakan beban statik yang sama sepanjang waktu.

# beban yang bekerja
timeSeries('Linear',1)
pattern('Plain',1,1)
load(2,0,-P,0)

(6) Keenam, pengaturan analisis dan running analisa

# analisis
system("BandSPD")
numberer("RCM")
constraints("Plain")
integrator("LoadControl", 1.0)
algorithm("Linear")
analysis("Static")
analyze(1)  #running analisa

(7) Ketujuh, melihat hasil respon struktur. Pada contoh ini hasil lendutan di tengah bentang (node 2), arah sumbu y atau 2.

# lendutan yang terjadi
uy = nodeDisp(2,2)
print(uy)

Hasil lendutan = 1,92 mm, sama persis dengan rumus praktis. Dengan ini model dan script OpenSees dianggap sudah benar.

Adapun analisa menggunakan OpenSeesPy secara keseluruhan dilakukan dengan script sebagai berikut. Penjelasan tahapan setiap perintah dapat dilihat di bawah ini.

from openseespy.opensees import *
wipe()
P = 10000
L = 6000
b = 250
h = 500
A = b*h
Iz = (b*h**3)/12
E = 9000
model('basic','-ndm',2,'-ndf',3)
node(1, 0.0, 0.0) ; fix(1,1,1,0)
node(2, L/2, 0.0)
node(3, L, 0.0)   ; fix(3,0,1,0)
geomTransf('Linear',1)
element('elasticBeamColumn',1,1,2,A,E,Iz,1)
element('elasticBeamColumn',2,2,3,A,E,Iz,1)
timeSeries('Linear',1)
pattern('Plain',1,1)
load(2,0,-P,0)
system("BandSPD")
numberer("RCM")
constraints("Plain")
integrator("LoadControl", 1.0)
algorithm("Linear")
analysis("Static")
analyze(1)
uy = nodeDisp(2,2)
print(uy)

Untuk perintah pembuatan model material, elemen, tumpuan, dsb; secara lengkap dapat dilihat pada koleksi perintah model (Model Commands) disini. Sedangkan untuk opsi pengaturan analysis dapat dilihat di sini.

Pada contoh ini, respons struktur yang dicek hanyalah lendutan. Berbagai jenis respons (reaksi, gaya, dsb) pada struktur (elemen dan node) dapat menggunakan perintah ouput (Output Command) di sini.

* Tahapan di atas bukan tahapan baku yang harus diikuti. Tapi merupakan tahapan untuk penyederhanaan agar lebih mudah diikuti.

Demikian semoga bermanfaat.

Advertisement
by Agus Daud

Nyoba Stability Challenge

Kali ini mau nyoba stability challenge yang diadakan oleh blog opensees di sini A Verry Stable Challenge. Ini bukan masalah jawabannya, karena jawabannya bisa dibaca di sini Stability Challenge Results.

Ini lebih kepada bagaimana cara untuk menyelesaikannya, bagaimana script opensees yang sesuai dengan strategi yang diberikan oleh blog tersebut?

I’m just try to understand the strategy written on the blog and try to write an opensees script that seems to fit the problem.

Baik, soalnya sebagai berikut:

Berapakah nilai faktor pengali (Lambda), sehingga terjadi perpindahan horizontal 0,25 inch di puncak struktur?

Strategi pertama,

yaitu dengan menganggap struktur sebagai rangka batang elastik (truss) biasa.

Analisis dilakukan dengan menggunakan openseespy, maka ditulis script dengan bahaya Python sebagai berikut:

from openseespy.opensees import *
import openseespy.preprocessing.DiscretizeMember as opsdm
import openseespy.postprocessing.Get_Rendering as opsplt
import matplotlib.pyplot as plt
import numpy as np
%matplotlib notebook

wipe()
model('basic','-ndm',2,'-ndf',2)

E = 3000.0
Px = 100.0
Py = 50.0

node(1,0.0,0.0)   ; fix(1,1,1,0)
node(2,144.0,0.0) ; fix(2,1,1,0)
node(3,168.0,0.0) ; fix(3,1,1,0)
node(4,72.0,96.0)

uniaxialMaterial("Elastic", 1, E)
element("Truss",1,1,4,10.0,1)
element("Truss",2,2,4,5.0,1)
element("Truss",3,3,4,5.0,1)

timeSeries("Linear",1)
pattern("Plain",1,1)

# input beban
load(4,Px,-Py)

dMax = 0.25
Nsteps = 50

# analisis
system("BandSPD")
numberer("RCM")
constraints("Plain")
integrator("DisplacementControl",4,1,dMax/Nsteps)
algorithm("Linear")
test("NormUnbalance",1e-8,30)
analysis("Static")

disp1 = []
factor1 = []

currentDisp = 0.0
ok = 0

while ok == 0 and currentDisp < dMax:
    ok = analyze(1)
    if ok == 0:
        test("NormDispIncr",1e-8,1000,0)
        algorithm("ModifiedNewton",'-initial')
        ok = analyze(1)
        test("NormDispIncr",1e-8,6,0)
        algorithm("Newton")
    currentDisp = nodeDisp(4,1)
    disp1.append(nodeDisp(4,1))
    factor1.append(getLoadFactor(1))
    
print(disp1[-1],factor1[-1])
plt.xlabel('Perpindahan Horizontal (inch)')
plt.ylabel('Faktor Pengali (Lambda)')
plt.plot(disp1,factor1,Label="Linear Truss")
plt.grid()
plt.legend()
plt.show()

Hasilnya sebagai berikut:

Dengan strategi pertama ini, untuk memperoleh perpindahan horizontal 0,25 inch di puncak struktur; diperoleh faktor pengali Lambda = 0,47.

Nilai ini sama dengan tiga partisipan lain. Berikut hasil dari partisipan lain.

Strategi kedua

Dijelaskan dalam blog tersebut, bahwa tiga partisipan mempertimbangkan non linear geometri, transformasi P-∆ atau Corotational; kemudian merubah elemen rangka menjadi elemen balok (beam). Sedangkan hasilnya akan tetap sama seperti strategi pertama yang hanya menggunakan elemen rangka.

Soal ini adalah masalah P-δ bukan P-∆; terlihat betapa pentingnya ketelitian melihat soal/kondisi struktur yang akan dianalisis.

Efek P-δ disimulasikan dengan membagi elemen balok (beam) menjadi beberapa elemen diskrit. Misalnya setiap balok dibagi menjadi 4 elemen.

Elemen diskrit

Inilah yang dilakukan oleh dua partisipan yang lain memperoleh nilai yang lebih kecil yaitu sebesar Lambda = 0,11.

Berikut script analisis dengan openseespy.

from openseespy.opensees import *
import openseespy.preprocessing.DiscretizeMember as opsdm
import openseespy.postprocessing.Get_Rendering as opsplt
import matplotlib.pyplot as plt
import numpy as np
import math
%matplotlib notebook

wipe()
model('basic','-ndm',2,'-ndf',3)

A = [10.0,5.0]
Iz = []
E = 3000.0
Px = 100.0
Py = 50.0

for val in A:
    b = math.sqrt(val)
    Iz.append(b**4/12)

node(1,0.0,0.0)   ; fix(1,1,1,0)
node(2,144.0,0.0) ; fix(2,1,1,0)
node(3,168.0,0.0) ; fix(3,1,1,0)
node(4,72.0,96.0)

geomTransf("Corotational",1)
section('Elastic',1,E,A[0],Iz[0])
section('Elastic',2,E,A[1],Iz[1])
beamIntegration('Legendre',1,1,10)
beamIntegration('Legendre',2,2,10)

#(ndI, ndJ, numEle, eleType, integrTag, transfTag, nodeTag, eleTag)
nEle = 4
opsdm.DiscretizeMember(1,4,nEle,'dispBeamColumn',1,1,5,1)
opsdm.DiscretizeMember(2,4,nEle,'dispBeamColumn',2,1,8,5)
opsdm.DiscretizeMember(3,4,nEle,'dispBeamColumn',2,1,11,9)

timeSeries("Linear",1)
pattern("Plain",1,1)
load(4,Px,-Py,0.0)

# analisis
Nsteps = 50
dMax = 0.25

integrator("DisplacementControl",4,1,dMax/Nsteps)
system("BandGeneral")
numberer("RCM")
constraints("Plain")

algorithm("Newton")
#test("NormDispIncr",1e-8,30)
analysis("Static")

disp2 = []
factor2 = []

currentDisp = 0.0
ok = 0

while ok == 0 and currentDisp < dMax:
    ok = analyze(1)
    if ok != 0:
        #test("NormDispIncr",1e-8,10000,0)
        #algorithm("ModifiedNewton",'-initial')
        ok = analyze(1)
        test("NormDispIncr",1e-8,100,0)
        algorithm("Newton")
        
    currentDisp = nodeDisp(4,1)
    disp2.append(nodeDisp(4,1))
    factor2.append(getLoadFactor(1))

print(disp2[-1],factor2[-1])
#opsplt.plot_model()

plt.plot(disp1,factor1)
plt.plot(disp2,factor2)
plt.xlabel('Perpindahan Horizontal (inch)')
plt.ylabel('Faktor Pengali (Lambda)')
plt.grid()
plt.show()

Hasilnya.. dengan cara ini diperoleh faktor pengali Lambda = 0,11. Sama dengan dua partisipan lain.

Merubah elemen rangka (truss) menjadi balok (beam), akan menyebabkan terjadinya momen pada join puncak. Kondisi ini disiasati dengan menambahkan node ekstra tambahan di puncak (saling menimpa); kemudian dihubungkan dengan perintah equalDOF.

Sayangnya nilai Lambda = 0,11 lebih besar dari batas bawah yaitu sebesar Lambda = 0,0606.

Nilai ini merupakan rasio kapasitas tekuk elastik (tekuk Euler) batang bagian paling kanan Pcr = 3,35 kip , dengan gaya tekan yang terjadi (pada Lambda = 1) sebesar Pc = 55,3 kip. Atau Lambda = Pcr / Pc = 3,35/55,3 = 0,0606.

Untuk mendapatkan hasil dengan nilai batas bawah Lambda = 0,0606, maka perlu ditambahkan nilai δ awal (ketidak-sempurnaan geometri) di node tengah setiap elemen sebesar L/1000 sampai L/500.

Need to add:
– interior nodes of the corotational mesh with a small out of straightness L/1000 to L/500.

Update 19 Agustus 2021

Di kolom komentar Prof. Scott menyarankan metode yang bisa digunakan untuk me-release momen di node puncak.

Selain itu, untuk memasukkan ketidak-sempurnaan geometri (initial imperfection); saya mencoba cara lain yaitu dengan memasukkan beban notional, Ni (horizontal) di tengah balok sebesar Pc/500 = 55,3/500 = 0,1.

Namun, karena pembebanan menggunakan faktor pengali; maka untuk mendapatkan gaya notional yang sesuai, gaya notional dibagi dengan faktor beban atau Ni = 0,1/0,0606 = 1,7 Kip.

Instead using corotational mesh with a small out of straightness L/1000 to L/500; to simulate initial imperfection, I try to put notional load at mid interior node of compression member by

Notional force, Ni = (Pc/500) / 0,0606 = 1,7 Kip

Berikut script opensees yang digunakan,

from openseespy.opensees import *
import openseespy.preprocessing.DiscretizeMember as opsdm
import openseespy.postprocessing.Get_Rendering as opsplt
import matplotlib.pyplot as plt
import numpy as np
import math
%matplotlib notebook

wipe()
model('basic','-ndm',2,'-ndf',3)

A = [10.0,5.0]
Iz = []
E = 3000.0
Px = 100.0
Py = 50.0

for val in A:
    b = math.sqrt(val)
    Iz.append(b**4/12)

node(1,0.0,0.0)   ; fix(1,1,1,0)
node(2,144.0,0.0) ; fix(2,1,1,0)
node(3,168.0,0.0) ; fix(3,1,1,0)
node(4,72.0,96.0) ; fix(4,0,0,1)
node(144,72.0,96.0); equalDOF(4,144,1,2)
node(244,72.0,96.0); equalDOF(4,244,1,2)
node(344,72.0,96.0); equalDOF(4,344,1,2)

geomTransf("Corotational",1)
section('Elastic',1,E,A[0],Iz[0])
section('Elastic',2,E,A[1],Iz[1])
beamIntegration('Legendre',1,1,10)
beamIntegration('Legendre',2,2,10)

#(ndI, ndJ, numEle, eleType, integrTag, transfTag, nodeTag, eleTag)
nEle = 4
opsdm.DiscretizeMember(1,144,nEle,'dispBeamColumn',1,1,5,1)
opsdm.DiscretizeMember(2,244,nEle,'dispBeamColumn',2,1,8,5)
opsdm.DiscretizeMember(3,344,nEle,'dispBeamColumn',2,1,11,9)

timeSeries("Linear",1)
pattern("Plain",1,1)
load(4,Px,-Py,0.0)
load(9,1.7,0.0,0.0)
load(12,1.7,0.0,0.0)


# analisis
Nsteps = 50
dMax = 0.25

integrator("DisplacementControl",4,1,dMax/Nsteps)
system("BandGeneral")
numberer("RCM")
constraints("Plain")
algorithm("Newton")
analysis("Static")

disp3 = []
factor3 = []

currentDisp = 0.0
ok = 0

while ok == 0 and currentDisp < dMax:
    ok = analyze(1)
    if ok != 0:
        test("NormDispIncr",1e-8,10000,0)
        algorithm("ModifiedNewton",'-initial')
        ok = analyze(1)
        test("NormDispIncr",1e-8,100,0)
        algorithm("Newton")
        
    currentDisp = nodeDisp(4,1)
    disp3.append(nodeDisp(4,1))
    factor3.append(getLoadFactor(1))

print(disp3[-1],factor3[-1])

plt.plot(disp1,factor1,Label="Linear Truss")
plt.plot(disp2,factor2,Label="NL Geometri")
plt.plot(disp3,factor3,Label="Non Geometri, M.Release, & init. imperfect")
plt.xlabel('Perpindahan Horizontal (inch)')
plt.ylabel('Faktor Pengali (Lambda)')
plt.grid()
plt.legend()
plt.show()

Hasilnya.. resulte Lambda = 0,060596

Resulte

Dengan ini diperoleh nilai Lambda = 0,060596 mendekati nilai batas bawah yaitu sebesar Lambda = 0,0606.

by Agus Daud

Opensees: Analisis Elastik Orde-2

Sebelumnya telah dibahas mengenai anaStruct: Tekuk Elastik dan Opensees: Perkenalan. Maka, kali ini saya ingin mencoba eksplor bagaimana sih, hasil analisis elastic buckling menggunakan OpenSeespy.

Sama dengan anaStruct, Openseespy juga merupakan program open source alias gratis dan tentunya non-GUI. Namun, Openseespy dibekali dengan banyak fasilitas hingga mampu memodelkan struktur, geoteknik maupun kedua secara bersamaan. Lebih dikenal dengan soil structure interaction (SSI).

Dengan soal yang sama yaitu kolom baja sederhana 8,5 m, dalam artikel anaStruct: Tekuk Elastik, dilakukan analisis elastik orde-2. Berikut souce code yang ditulis menggunakan bahasa python3 dan package openseespy.

import openseespy.opensees as ops
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import figure
import openseespy.postprocessing.Get_Rendering as opsplt
import openseespy.preprocessing.DiscretizeMember as opsdm
%matplotlib inline

ops.wipe()
ops.model('basic', '-ndm', 2, '-ndf', 3)

cm = 1/100. # m
kg = 0.00981 #kN

# Data Penampang
A = 39.65 *cm**2
I = 563. *cm**4
MPa = 1000. # kPa
w = 31.1 *kg

# Data material
E0 = 200000.*MPa
L = 8.5

Nsteps = 25
Py = -160.
Px = 0.002*160

# membuat node
ops.node(1, 0.0, 0.0)
ops.node(2, 0.0, L)

# tumpuan
ops.fix(1, 1,1,0)
ops.fix(2, 1,0,0)

# elemen
ops.geomTransf('Corotational', 1)
ops.section('Elastic',1,E0,A,I)
ops.beamIntegration('Lobatto',1,1,10)
opsdm.DiscretizeMember(1, 2, 8, 'elasticBeamColumn', 1, 1, 3, 1)

ops.timeSeries("Linear", 1)
ops.pattern("Plain", 1, 1)
ops.load(2, 0., Py, 0.)
ops.load(5, Px, 0., 0.)


# Analisis
ops.system("ProfileSPD") # create SOE
ops.numberer("Plain") # create DOF number
ops.constraints("Plain") # create constraint handler
ops.integrator("LoadControl", 1.0/Nsteps, 10) # create integrator
ops.algorithm("Newton") # create algorithm
ops.analysis("Static") # create analysis object

data = np.zeros((Nsteps+1,2))
for j in range(Nsteps):
    ops.analyze(1)
    data[j+1,0] = ops.nodeDisp(5,1)
    data[j+1,1] = ops.getLoadFactor(1)*Py*-1

#Plot data dalam grafik
figure(num=None, figsize=(5, 5), dpi=80, facecolor='w', edgecolor='k')
plt.plot(data[:,0], data[:,1],'-',label='anaStruct Pcr = 153.60 kN')
plt.plot([0,1.25],[153.815,153.815],'r--',label='Exact = 153.815 kN')
plt.xlabel('Lateral Displacement at Mid Point(m))')
plt.ylabel('Vertical Load at top (kN)')
plt.legend(bbox_to_anchor=(0.5, 0.2), loc='upper center', borderaxespad=0.)
plt.show()
Output

Dari grafik di atas, terlihat perubahan perpindahan lateral yang tiba-tiba dan besar, yang menandakan terjadinya tekuk. Hasil analisa diperoleh kuat tekuk kritis (Pcr) sebesar 153,6 kN. Maka, dapat dibandingkan dari hasil analisa sebelumnya sebagai berikut.

===> geser tabel

AnalisisEulerSAP2000anaStructOpensees
Eksak153,815
Linear Elastic153,820153,820
Elastic Orde-2153,600153,000153,600

Terlihat bahwa hasil yang diperoleh dengan SAP2000 sama persis dengan Openseespy, yaitu sebesar 153,60 kN. Dengan demikian pemodelan dianggap sudah mendekati nilai eksak sebesar 153,815 kN.

*Catatan:

Opensees adalah program yang ditulis dalam bahasa Tcl, sedangkan Openseespy merupakan program yang ditulis dalam bahasa Python3.6. Oleh karena itu untuk membuka opensees harus dipersiapkan terlebih dahulu persyaratannya. Dapat dilihat di halaman berikut:

  1. Opensees: https://opensees.berkeley.edu/
  2. Openseespy: https://openseespydoc.readthedocs.io/en/latest/