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.

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

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