Kelompok 1
Kelompok 1 Anggota:
- Andika Ramadhan Gurnida
- M. Pasha Wibisono
- M. Yusuf Abdurrahman
- Rizkyawan Wibawa
Tugas: Menyelesaikan problem set bab 2 dengan program Python
- Agenda untuk tugas ini berdasarkan pemahaman kami atas informasi yang ada adalah pembuatan program Python untuk melakukan eliminasi Gauss, yang kemudian akan digunakan untuk menyelesaikan salah satu soal pada problem set yang terdapat pada Bab 2.
- Kami memilih soal nomor 4 di Problem Set 2.1 dari buku untuk diselesaikan.
- Langkah-langkah yang terdapat dalam program yang kami buat adalah sebagai berikut:
- 1. Mengimpor modul numpy untuk membuat matriks
- 2. Membuat matriks persamaan yang terdapat pada soal untuk diselesaikan
- 3. Menulis algoritma untuk melakukan fase eliminasi
- 4. Menulis algoritma untuk melakukan back substitution
- 5. Mencetak matriks penyelesaian soal
- Program yang telah kami buat secara spesifik untuk menyelesaikan soal tersebut sebagai berikut:
import numpy as np
A = np.array ([[2,-3,-1],[3,2,-5],[2,4,-1]]) b = np.array ([[3],[-9],[-5]])
def gaussElimin(A,b):
n = len(b) for k in range (0,n-1): for i in range(k+1,n): if A[i,k] != 0.0: lam = A[i,k]/A[k,k] A[i,k+1:n] = A[i,k+1:n] - lam*A[k,k+1:n] b[i] = b[i] - lam*b[k]
for k in range(n-1,-1,-1): b[k] = (b[k] - np.dot(A[k,k+1:n],b[k+1:n]))/A[k,k] return b
x = gaussElimin(A,b) print (x)
Tugas: Aplikasi Metode Numerik pada masalah statika struktur
- Agenda untuk tugas ini berdasarkan pemahaman kami atas informasi yang ada adalah menjelaskan aplikasi dari apa yang telah dipelajari pada mata kuliah metode numerik untuk menyelesaikan permasalahan statika struktur satu dimensi.
- Kami memilih menggunakan analisa finite element untuk menyelesaikan contoh soal statika struktur satu dimensi.
- Soal tersebut sebagai berikut:
- Sebuah batang baja dengan panjang 600mm memiliki luas daerah cross-section sebesar 650mm^2 dan 350mm^2 pada kedua ujungnya. Ujung dengan luasan yang lebih besar menempel pada dinding dan diberikan dua gaya aksial sebesar 40kN pada ujung yang kecil dengan arah menjauhi dinding, dan 10kN pada posisi 200mm dari ujung yang kecil dengan arah mendekati dinding. Modelkan batang tersebut dengan tiga finite element dan tentukan nodal displacement yang terjadi.
- Pertama-tama, kami memodelkan batang tersebut dengan tiga finite element terlebih dulu.
- Ketiga elemen tersebut dipisahkan berdasarkan rata-rata luas daerah cross-section pada ujung-ujungnya, dengan panjang masing-masing elemen sebesar 200mm. Elemen yang paling dekat dengan dinding memiliki rata-rata luas daerah cross-section sebesar 600mm^2, elemen di tengah sebesar 500mm^2, dan elemen paling ujung sebesar 400mm^2.
- Titik-titik tempat menyatunya elemen-elemen tersebut diperlakukan sebagai node yang akan dihitung perpindahannya.
- Kita juga menentukan posisi konektivitas antar elemen, yang nantinya akan berguna untuk membuat matriks.
- Node yang dimiliki batang tersebut ada 4 buah (1 pada ujung setiap elemen), maka dapat kita tentukan bahwa elemen paling dekat dengan dinding memiliki konektivitas dengan dinding dan elemen tengah di titik 1 dan 2, elemen tengah terhubung dengan elemen dekat dinding dan elemen ujung di titik 2 dan 3, dan elemen ujung terhubung dengan elemen tengah dan bagian ujung batang di titik 3 dan 4.
- Informasi ini digunakan untuk membuat matriks global stiffness terlebih dahulu, dimana matriks tersebut diisi dengan nilai stiffness dari masing-masing elemen. Nilai tersebut ditentukan oleh rumus k = (E.A) / l, dimana k adalah nilai stiffness, E adalah modulus Young untuk baja (200 GPa = 200 x 10^3 MPa [N/mm^2]), A adalah luas daerah cross-section, dan l adalah panjang setiap elemen.
- Kita juga membuat matriks global load vector, dengan nilai beban sesuai pada beban yang terdapat di setiap node.
- Setelah kita membuat kedua matriks tersebut, maka nodal displacement dapat dicari (dengan catatan node dinding memiliki perpindahan 0 sehingga kita hanya akan mendapatkan tiga nilai perpindahan) dengan program berikut:
# Mengimpor modul untuk membuat matriks import numpy as np
# Membuat matriks stiffness dan matriks gaya pada objek K = np.array ([[5.5,-2.5,0],[-2.5,4.5,-2],[0,-2,2]]) F = np.array ([[0],[-0.05],[0.2]])
# Mendefinisikan eliminasi Gauss def gaussElimin(K,F):
# Algoritma fase eliminasi n = len(F) for k in range (0,n-1): for i in range(k+1,n): if K[i,k] != 0.0: lam = K [i,k]/K[k,k] K[i,k+1:n] = K[i,k+1:n] - lam*K[k,k+1:n] F[i] = F[i] - lam*F[k]
# Algoritma back substitution for k in range(n-1,-1,-1): F[k] = (F[k] - np.dot(K[k,k+1:n],F[k+1:n]))/K[k,k] return F
# Mendefinisikan variabel untuk dicetak pada layar Q = gaussElimin(K,F) print (Q)
- Catatan 1: nilai matriks beban dibagi dengan nilai modulus Young untuk mempermudah pembuatan matriks.
- Catatan 2: hasil matriks Q dalam satuan mm
- Catatan 3: untuk tugas aplikasi metode numerik untuk memodelkan masalah statika struktur, video tidak dapat kami upload ke wiki karena ukurannya yang terlalu besar. Kami mohon pemakluman dari Pak Indra dan Pak Radon mengenai hal ini, dan apabila berkenan, kami bersedia menunjukkan video yang telah kami buat secara langsung.
Tugas: Spring-mass system
- Soal: Diilustrasikan sebuah sistem pegas-massa yang terpasang pada suatu dinding. Pegas tersebut memiliki konstanta kekakuan 75 N/m dan massa benda adalah 2,5 kg. Gaya P diberikan pada sistem ke arah menjauhi dinding dengan besaran 10t N selama 2 detik pertama dan 20 N setelahnya. Persamaan untuk percepatan sistem diberikan sebagai berikut: y" = (P(t)/m) - ((k/m)y). Sistem awalnya berada dalam keadaan diam. Hitung perpindahan maksimum sistem, dan plotkan pada grafik perubahan posisi sistem terhadap waktu!
- Analisa fenomena fisika: Pada sistem pegas-massa akan berlaku hukum II Newton, dimana arah percepatan sistem sama dengan arah total gaya yang bekerja pada sistem. Kita tahu bahwa benda dengan massa memiliki perumusan gaya F = m.a, dimana m adalah massa benda dan a (atau y") adalah percepatan benda, dan pegas memiliki perumusan gaya F = k.y, dimana k adalah konstanta kekakuan pegas dan y adalah perubahan panjang pegas dari panjang awal. Pada sistem pegas-massa dalam soal, gaya yang diberikan pada sistem arahnya menjauhi dinding. Hal ini akan mengakibatkan massa bergerak menjauhi dinding dengan percepatan tertentu. Pegas yang terpasang pada dinding akan memberikan gaya reaksi akibat meregang, yang arahnya berlawanan dengan gaya yang bekerja pada sistem. Oleh karena itu, perumusan gaya pada sistem dapat dinyatakan dengan:
- Gaya yang bekerja pada benda = Gaya yang diberikan pada sistem - gaya reaksi pegas
- m.y" = P(t) - k.y
- Apabila kedua ruas kita bagi dengan m, maka akan kita dapatkan persamaan seperti yang tertulis pada soal
- y" = P(t)/m - (k/m)y
- Kami menggunakan aplikasi pycharm dan kodingan yang kami gunakan adalah:
import numpy as np; from numpy import sin, asarray, zeros, interp from run_kut4 import integrate from printSoln import printSoln import pylab
k = 75 m = 2.5
def P(t):
if t < 2:
return 10.*t
else:
return 20.
def F(x,y):
F = zeros( (2,), dtype=float )
F[0] = y[1]
F[1] = P(x)/m - (k/m) * y[0]
return F
x0 = 0.0 # start of the integration xStop = 6.0 # end of the integration y0 = asarray( [ 0., 0. ], dtype=float ) h = 0.01 freq = 20
X,Y = integrate(F,x0,y0,xStop,h) #printSoln(X,Y,freq)
pylab.plot( X, Y[:,0] )
pylab.xlabel('time (t)');
pylab.ylabel('y(t)');
pylab.axhline(y=0,color='green')
pylab.grid('on')
pylab.title("Kelompok 1")
pylab.show()
.jpg)
.jpg)
pertama kami install modul numpy dan install modul matplotlib setelah itu kami menggunakan modul run_kut4 dan juga modul pprintSoln yang dua-duanya terdapat di github.com
## module run_kut4
X,Y = integrate(F,x,y,xStop,h).
4th-order Runge-Kutta method for solving the
initial value problem {y}' = {F(x,{y})}, where
{y} = {y[0],y[1],...y[n-1]}.
x,y = initial conditions
xStop = terminal value of x
h = increment of x used in integration
F = user-supplied function that returns the
array F(x,y) = {y'[0],y'[1],...,y'[n-1]}.
from numpy import array
def integrate(F,x,y,xStop,h):
def run_kut4(F,x,y,h):
K0 = h*F(x,y)
K1 = h*F(x + h/2.0, y + K0/2.0)
K2 = h*F(x + h/2.0, y + K1/2.0)
K3 = h*F(x + h, y + K2)
return (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0
X = []
Y = []
X.append(x)
Y.append(y)
while x < xStop:
h = min(h,xStop - x)
y = y + run_kut4(F,x,y,h)
x = x + h
X.append(x)
Y.append(y)
return array(X),array(Y)
## module printSoln
printSoln(X,Y,freq).
Prints X and Y returned from the differential
equation solvers using printput frequency 'freq'.
freq = n prints every nth step.
freq = 0 prints initial and final values only.
def printSoln(X,Y,freq):
def printHead(n):
print "\n x ",
for i in range (n):
print " y[",i,"] ",
print
def printLine(x,y,n):
print "%13.4e"% x,
for i in range (n):
print "%13.4e"% y[i],
print
m = len(Y)
try: n = len(Y[0])
except TypeError: n = 1
if freq == 0: freq = m
printHead(n)
for i in range(0,m,freq):
printLine(X[i],Y[i],n)
if i != m - 1: printLine(X[m - 1],Y[m - 1],n)
Setelah semua telah dimasukkan dalam aplikasi pycharm maka kami melakukan run dan kemudian akan langsung muncul grafik.
Tugas optimasi aerodinamika
- Membuat geometri aerofoil di software CAD yang digunakan (kelompok kami menggunakan ukuran airfoil standar NACA dengan nomor 63-412).
- Melakukan simulasi aliran udara untuk geometri airfoil di software CFDSOF untuk mendapatkan nilai lift dan drag dari airfoil pada berbagai angle of attack.
- Menggunakan program Python yang ditujukan untuk optimasi nilai angle of attack, dimana nilai tersebut akan memiliki gaya lift besar dan drag kecil.
- Program Python yang digunakan:
import numpy as np from scipy.optimize import minimize
def calc_drag(x):#drag x1 = x[0] drag = 0.0193*x1**5 - 0.3797*x1**4 + 2.1688*x1**3 - 2.2064*x1**2 - 3.0016*x1 + 4.0659 return drag
def calc_lift(x): #lift x1 = x[0] lift = -0.0179*x1**5 + 0.55*x1**4 - 6.0283*x1**3 + 26.852*x1**2 - 38.461*x1 + 21.939 return lift
def objective(x): #volume yang diminimalkan return calc_drag(x)
def constraint1(x): #variable SUDUT yang meminimalkan persamaan garis drag return calc_drag(x) def constraint2(x): #variable SUDUT yang meminimalkan persamaan garis lift return calc_lift(x)
con1=({'type':'ineq','fun':constraint1})
con2=({'type':'ineq','fun':constraint2})
cons = (con1,con2)
a = (0,90)
batas = [a]
x1_guess = 50
x0 = np.array([x1_guess])
sol = minimize(objective,x0, method='SLSQP',bounds=batas, constraints=cons, options={'disp':True})
xopt = sol.x forceopt = -sol.fun
dragopt = calc_drag(xopt) # drag optimal liftopt = calc_lift(xopt) # lift optimal
print ('sudut optimal = '+str(xopt[0]))
print ('total force optimal = '+str(forceopt))
print ('drag force optimal = '+str(dragopt))
print ('lift force optimal = '+str(liftopt))
# In[10]:
import numpy as np from scipy.optimize import minimize
def calc_drag(x):#drag x1 = x[0] drag = 0.0193*x1**5 - 0.3797*x1**4 + 2.1688*x1**3 - 2.2064*x1**2 - 3.0016*x1 + 4.0659 return drag
def calc_lift(x): #lift x1 = x[0] lift = -0.0179*x1**5 + 0.55*x1**4 - 6.0283*x1**3 + 26.852*x1**2 - 38.461*x1 + 21.939 return lift
def objective(x): #volume yang diminimalkan return calc_lift(x)
def constraint1(x): #variable SUDUT yang meminimalkan persamaan garis drag return calc_drag(x) def constraint2(x): #variable SUDUT yang meminimalkan persamaan garis lift return calc_lift(x)
con1=({'type':'ineq','fun':constraint1})
con2=({'type':'ineq','fun':constraint2})
cons = (con1,con2)
a = (0,90)
batas = [a]
x1_guess = 50
x0 = np.array([x1_guess])
sol = minimize(objective,x0, method='SLSQP',bounds=batas, constraints=cons, options={'disp':True})
xopt = sol.x forceopt = -sol.fun
dragopt = calc_drag(xopt) # drag optimal liftopt = calc_lift(xopt) # lift optimal
print ('sudut optimal = '+str(xopt[0]))
print ('total force optimal = '+str(forceopt))
print ('drag force optimal = '+str(dragopt))
print ('lift force optimal = '+str(liftopt))
