Difference between revisions of "Muhammad Rayhan Faturrahman"

From ccitonlinewiki
Jump to: navigation, search
(Update Tugas Besar)
(Update Tugas Besar)
Line 780: Line 780:
  
 
== Update Tugas Besar ==
 
== Update Tugas Besar ==
 +
 +
model TugasBesar_Trusses_3d_Safety
 +
//define initial variable
 +
parameter Integer Points=size(P,1); //Number of Points
 +
parameter Integer Trusses=size(C,1); //Number of Trusses
 +
parameter Real Yield=215e6; //Yield Strength (Pa)
 +
parameter Real Area=0.000171;  //Area L Profile (Dimension=0.03, Thickness=0,003) (m2)
 +
parameter Real Elas=193e9;    //Elasticity SS 304  (Pa)
 +
//define connection
 +
parameter Integer C[:,2]=[1,5;
 +
                        2,6;
 +
                        3,7;
 +
                        4,8;
 +
                        5,6;  //1st floor
 +
                        6,7;  //1st floor
 +
                        7,8;  //1st floor
 +
                        5,8;  //1st floor
 +
                        5,9;
 +
                        6,10;
 +
                        7,11;
 +
                        8,12;
 +
                        9,10; //2nd floor
 +
                        10,11;//2nd floor
 +
                        11,12;//2nd floor
 +
                        9,12; //2nd floor
 +
                        9,13;
 +
                        10,14;
 +
                        11,15;
 +
                        12,16;
 +
                        13,14;//3rd floor
 +
                        14,15;//3rd floor
 +
                        15,16;//3rd floor
 +
                      13,16];//3rd floor                                                         
 +
//define coordinates (please put orderly)
 +
parameter Real P[:,6]=[0.3,-0.375,0,1,1,1;    //1
 +
                      -0.3,-0.375,0,1,1,1;    //2
 +
                      -0.3,0.375,0,1,1,1;    //3
 +
                      0.3,0.375,0,1,1,1;      //4
 +
                         
 +
                      0.3,-0.375,0.6,0,0,0;  //5
 +
                      -0.3,-0.375,0.6,0,0,0;  //6
 +
                      -0.3,0.375,0.6,0,0,0;  //7
 +
                      0.3,0.375,0.6,0,0,0;    //8
 +
                         
 +
                      0.3,-0.375,1.2,0,0,0;  //9
 +
                      -0.3,-0.375,1.2,0,0,0;  //10 
 +
                      -0.3,0.375,1.2,0,0,0;  //11
 +
                      0.3,0.375,1.2,0,0,0;    //12
 +
                         
 +
                      0.3,-0.375,1.8,0,0,0;  //13
 +
                      -0.3,-0.375,1.8,0,0,0;  //14
 +
                      -0.3,0.375,1.8,0,0,0;  //15
 +
                      0.3,0.375,1.8,0,0,0];  //16                       
 +
//define external force (please put orderly)
 +
parameter Real F[Points*3]={0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,-500,
 +
                          0,0,-1000,
 +
                          0,0,-1000,
 +
                          0,0,-500};
 +
//solution
 +
Real displacement[N], reaction[N];
 +
Real check[3];
 +
Real stress1[Trusses];
 +
Real safety[Trusses];
 +
Real dis[3];
 +
Real Str[3];
 +
protected
 +
parameter Integer N=3*Points;
 +
Real q1[3], q2[3], g[N,N], G[N,N], G_star[N,N], id[N,N]=identity(N), cx, cy, cz, L, X[3,3];
 +
Real err=10e-10, ers=10e-4;
 +
algorithm
 +
//Creating Global Matrix
 +
G:=id;
 +
for i in 1:Trusses loop
 +
for j in 1:3 loop
 +
q1[j]:=P[C[i,1],j];
 +
q2[j]:=P[C[i,2],j];
 +
end for;   
 +
  //Solving Matrix
 +
  L:=Modelica.Math.Vectors.length(q2-q1);
 +
  cx:=(q2[1]-q1[1])/L;
 +
  cy:=(q2[2]-q1[2])/L;
 +
  cz:=(q2[3]-q1[3])/L;
 +
  X:=(Area*Elas/L)*[cx^2,cx*cy,cx*cz;
 +
                    cy*cx,cy^2,cy*cz;
 +
                    cz*cx,cz*cy,cz^2];
 +
  //Transforming to global matrix
 +
  g:=zeros(N,N);
 +
  for m,n in 1:3 loop
 +
    g[3*(C[i,1]-1)+m,3*(C[i,1]-1)+n]:=X[m,n];
 +
    g[3*(C[i,2]-1)+m,3*(C[i,2]-1)+n]:=X[m,n];
 +
    g[3*(C[i,2]-1)+m,3*(C[i,1]-1)+n]:=-X[m,n];
 +
    g[3*(C[i,1]-1)+m,3*(C[i,2]-1)+n]:=-X[m,n];
 +
  end for; 
 +
G_star:=G+g;
 +
G:=G_star;
 +
end for;
 +
//Implementing boundary
 +
for x in 1:Points loop
 +
if P[x,4] <> 0 then
 +
  for a in 1:Points*3 loop
 +
    G[(x*3)-2,a]:=0;
 +
    G[(x*3)-2,(x*3)-2]:=1;
 +
  end for;
 +
end if;
 +
if P[x,5] <> 0 then
 +
  for a in 1:Points*3 loop
 +
    G[(x*3)-1,a]:=0;
 +
    G[(x*3)-1,(x*3)-1]:=1;
 +
  end for;
 +
end if;
 +
if P[x,6] <> 0 then
 +
  for a in 1:Points*3 loop
 +
    G[x*3,a]:=0;
 +
    G[x*3,x*3]:=1;
 +
  end for;
 +
end if;
 +
end for;
 +
//Solving displacement
 +
displacement:=Modelica.Math.Matrices.solve(G,F);
 +
//Solving reaction
 +
reaction:=(G_star*displacement)-F;
 +
//Eliminating float error
 +
for i in 1:N loop
 +
reaction[i]:=if abs(reaction[i])<=err then 0 else reaction[i];
 +
displacement[i]:=if abs(displacement[i])<=err then 0 else displacement[i];
 +
end for;
 +
//Checking Force
 +
check[1]:=sum({reaction[i] for i in (1:3:(N-2))})+sum({F[i] for i in (1:3:(N-2))});
 +
check[2]:=sum({reaction[i] for i in (2:3:(N-1))})+sum({F[i] for i in (2:3:(N-1))});
 +
check[3]:=sum({reaction[i] for i in (3:3:N)})+sum({F[i] for i in (3:3:N)});
 +
for i in 1:3 loop
 +
check[i] := if abs(check[i])<=ers then 0 else check[i];
 +
end for;
 +
//Calculating stress in each truss
 +
for i in 1:Trusses loop
 +
for j in 1:3 loop
 +
q1[j]:=P[C[i,1],j];
 +
q2[j]:=P[C[i,2],j];
 +
dis[j]:=abs(displacement[3*(C[i,1]-1)+j]-displacement[3*(C[i,2]-1)+j]);
 +
end for;   
 +
  //Solving Matrix
 +
  L:=Modelica.Math.Vectors.length(q2-q1);
 +
  cx:=(q2[1]-q1[1])/L;
 +
  cy:=(q2[2]-q1[2])/L;
 +
  cz:=(q2[3]-q1[3])/L;
 +
  X:=(Elas/L)*[cx^2,cx*cy,cx*cz;
 +
              cy*cx,cy^2,cy*cz;
 +
              cz*cx,cz*cy,cz^2]; 
 +
  Str:=(X*dis);
 +
  stress1[i]:=Modelica.Math.Vectors.length(Str);
 +
end for;
 +
//Safety factor
 +
for i in 1:Trusses loop
 +
if stress1[i]>0 then
 +
  safety[i]:=Yield/stress1[i];
 +
else
 +
  safety[i]:=0;
 +
end if;
 +
end for;
 +
end TugasBesar_Trusses_3d_Safety;
  
  
 
dan didapati hasil berikut dari simulasi tersebut
 
dan didapati hasil berikut dari simulasi tersebut

Revision as of 16:46, 6 January 2021

BIODATA

Nama : Muhammad Rayhan Faturrahman

NPM : 1906301255

Jurusan : Teknik Mesin

Angkatan : 2019

Hobi : Kulineran

Jargon : Selaras Dalam Dinamika


Selamat sore teman-teman, perkenalkan nama saya Fatur, teman-teman saya di Departemen Teknik Mesin biasa memanggil saya bob. Saya memilih jurusan Teknik Mesin karena terisnpirasi dari saudara saya yang seorang process engineer, saya sering memerhatikannya melakukan simulasi desain pabrik dan menurut saya itu merupakan pekerjaan yang seru sekaligus menantang, berangkat dari sana saya memilih untuk menjadi seorang engineer juga tapi pada disiplin yang berbeda dengan harapan nantinya kami dapat membuat perusahaan konsultan sendiri.

Komputasi

Saya memiliki ketertarikan pada dunia komputasi sejak awal memasuki perkuliahan teknik, namun masih muncul banyak pertanyaan "harus mulai belajar dari mana?" dan pertanyaan lain sebagainya, alasan inilah yang semata-mata menjadi pemicu saya untuk mencoba mengambil mata kuliah metode numerik di semester 3 ini.

Perkuliahan Minggu Pertama (Rabu, 11 November 2020)

Pada pertemuan pertama yang diisi oleh Bapak Dr. Ir. Ahmad Indra Siswantara. Mahasiswa dijelaskan mengenai esensi atau tujuan pembelajaran metode numerik, yaitu:

1. Memahami konsep-konsep dan prinsip-prinsip dasar dalam metode numerik seperti:Persamaan algoritma, aljabar, pencocokan kurva, persamaan diferensial parsial, dan lainnya

2. Mengerti dan mampu menerapkan pemahaman atau aplikasi terhadap konsep metode numerik

3. Mampu menerapkan metode numerik dalam hal terkait persoalan keteknikan

4. Mendapat added value (nilai tambah)/adab sehingga mahasiswa menjadi orang yang memiliki budi pekerti yang baik

Selain empat poin yang telah disebutkan, Mahasiswa juga mendapatkan pelajaran terkait penggunaan aplikasi OpenModelica. Dimana Mahasiswa dijelaskan dengan sangat sistematis mengenai penggunaan sederhana OpenModelica.

Untuk memenuhi tugas minggu pertama terkait pengantar aplikasi OpenModelica, saya telah membuat model matematis untuk menyelesaikan problem fisika dengan tema Hukum Kekekalan Momentum

LINK TUGAS : https://youtu.be/JjA5Otn82pY


Perkuliahan Minggu Ke-2 (18 November 2020)

Catatan Penting

Menjadi pribadi yang lebih baik dari hari kemarin, bahwa menjadi manusia haruslah terus berkembang setiap harinya, merupakan pesan moral sekaligus pelajaran hidup pembuka kuliah yang disampaikan Pak Dai pada perkuliahan kali ini.

Beberapa hal penting yang pahami pada perkuliahan kali ini adalah:

- Metode Numerik sangat mengassist perhitungan di bidang engineering, faktor cepat dan akurat memegang peranan penting dalam penggunaannya di lapangan.

- Efisiensi adalah alasan utama mengapa kita sebagai calon engineer semestinya memiliki keahlian metode numerik, terutama pemodelan/

- Berbicara mengenai pemodelan, Pak Dai sempat membahas mengenai mengapa pada perkuliahan metode numerik beliau memilih menggunakan OpenModellica dibanding bahasa pemrograman lain, seperti Fortran dan C, beliau menekankan bahwa Modelica merupakan bahasa pemodelan, hal ini menjadikannya tools yang mudah dioperasikan oleh kami yang berlatar belakang bukan dari Computer Science.

- Selain itu, Modelica adalah aplikasi yang bersifat Open-Source alias gratis, sehingga mahasiswa tidak perlu mengeluarkan biaya sepeserpun.

Tugas Minggu Ke-2

Pak Dai menugaskan kami untuk membuat sebuah class dan fungsi panggil. Disini saya akan mengerjakan sebuah soal Aljabar Linier dan mencari penyelesaiannya menggunakan perintah Modelica.Math.Matrices.solve


Perkuliahan Minggu Ke-3 (25 November 2020)

Pada perkuliahan ini Pak Dai menjelaskan tentang beberapa pengaplikasian metode numerik di dunia engineering, antara lain

- Computational Fluid Dynamic

- Metode Stokastik

- Finite Element Analysis

Sebelum membahas tentang aplikasi metode numerik pada permasalahan engineering, kita mereview tugas minggu ke-2 mengenai class dan function panggil, kemudian ternyata satu kelas masih kurang mencapai target Pak Dai, karena masih menggunakan function yang sudah disediakan oleh OpenModelica dan bukan membuat sendiri. Akhirnya cristo mencoba menjelaskan penyelesaian yang sudah dia buat seperti di bawah ini :

// Gauss-Jordan Algorithm
// Transforms input matrix A into reduced row echelon form matrix B
// Christopher S.E. November 2020
input Real [:,:] A;     // An augmented matrix of m*n
output Real [:,:] B;    // Output matrix in reduced row echelon form
// Local variables
protected
Integer h = 1;          // Initialize pivot row
Integer k = 1;          // Initialize pivot column
Integer m = size(A,1);  // Number of rows in matrix
Integer n = size(A,2);  // Number of columns in matrix
Integer c = 0;          // Index counter
Integer max_row;        // Row index of max number in pivot column
Real [:] pivot_column;  // Vector containing pivot column data
Real [:] pivot_row;     // Vector containing backwards pivot row data
Real [:,:] temp_array;  // Stores matrix data when switching rows
Real r;                 // Ratio value used for row operations
// Limit to handle floating point errors
Real float_error = 10e-10;
algorithm
// Transfer input matrix A into variable B
B := A;
while h <= m and k <= n loop
 
// Dealing with floating point errors
for i in 1:m loop
   for j in 1:n loop
     if abs(B[i,j]) <= float_error then
       B[i,j] := 0;
     end if;
   end for;
 end for;
 // Finding the pivot
   pivot_column := {B[i,h] for i in h:m};
 
   // Get position index of lowest row with greatest pivot number
   c:= h-1;
   for element in pivot_column loop
     c := c+1;
     if abs(element) == max(abs(pivot_column)) then
       max_row := c;
     end if;
   end for;
 
 // No pivot in this column, move on to next column
 if B[max_row, k] == 0 then
   k := k+1;
 
 else
   // Swap Rows h <-> max_row
   temp_array := B;
   temp_array[h] := B[max_row];
   temp_array[max_row] := B[h];
   B:= temp_array;
   
   // Divide pivot row by pivot number
   B[h] := B[h]/B[h,k];
   
   // For all rows below the pivot
   for i in (h+1):m loop
     
     // Store the ratio of the row to the pivot
     r := B[i,k] / B[h,k];
     
     // Set lower part of pivot column to zero
     B[i,k] := 0;
     
     // Operations on the remaining row elements
     for j in (k+1):n loop
         B[i,j] := B[i,j] - B[h,j] * r;
     end for;
     
   end for;
   
   // Move on to next pivot row and column
   h := h+1;
   k := k+1;
   
 end if;
 
end while;
// The matrix is now in row echelon form
// Set values of (h,k) to (m,n)
h := m;
k := n;
while h >= 1 and k >=1 loop
 // Dealing with floating point errors
 for i in 1:m loop
   for j in 1:n loop
     if abs(B[i,j]) <= float_error then
       B[i,j] := 0;
     end if;
   end for;
 end for;
 // Finding the pivot
   pivot_row := {B[h,i] for i in 1:k};
   
   // Get position index k of pivot
   c := 0;
   for element in pivot_row loop
     c := c+1;
     if element <> 0 then
       break;
     end if;
   end for;
   k := c;
 // No pivot in this row, move on to next row
 if B[h, k] == 0 then
Bold text   h := h-1;
   
 else
   
   // Perform row operations
   for i in 1:(h-1) loop
     r := B[i,k];
     B[i] := B[i] - B[h] * r;
   end for;
   
   // Move on to next pivot row and column
   h := h-1;
   k := k-1;
   
 end if;
end while;
// The matrix is now in reduced row echelon form
end GaussJordan; 

Tugas Minggu Ke-3

Soal m3.jpg

Persamaan

class Trusses_Tugas3_MRayhanFatur

parameter Integer N=8; //Global matrice = 2*points connected
parameter Real A=0.001; //Area m2
parameter Real E=200e9; //Pa
Real G[N,N]; //global
Real Ginitial[N,N]; //global
Real Sol[N]; //global dispplacement
Real X[N]={0,0,-1035.2762,-3863.7033,0,0,-1035.2762,-3863.7033};
Real R[N]; //global reaction force
Real SolMat[N,1];
Real XMat[N,1];

//boundary condition
Integer b1=1;
Integer b2=3;

//truss 1
parameter Real X1=0; //degree between truss
Real k1=A*E/1;
Real K1[4,4]; //stiffness matrice
Integer p1a=1;
Integer p1b=2;
Real G1[N,N];

//truss 2
parameter Real X2=0; //degree between truss
Real k2=A*E/1;
Real K2[4,4]; //stiffness matrice
Integer p2a=2;
Integer p2b=3;
Real G2[N,N];

//truss 3
parameter Real X3=90; //degree between truss
Real k3=A*E/1.25;
Real K3[4,4]; //stiffness matrice
Integer p3a=2;
Integer p3b=4;
Real G3[N,N];

//truss 4
parameter Real X4=90+38.6598; //degree between truss
Real k4=A*E/1.6;
Real K4[4,4]; //stiffness matrice
Integer p4a=1;
Integer p4b=4;
Real G4[N,N];

//truss 5
parameter Real X5=90-38.6598; //degree between truss
Real k5=A*E/1.6;
Real K5[4,4]; //stiffness matrice
Integer p5a=3;
Integer p5b=4;
Real G5[N,N];

/*
for each truss, please ensure pXa is lower then pXb (X represents truss element number)
*/

algorithm

//creating global matrice
K1:=Stiffness_Matrices(X1);
G1:=k1*Local_Global(K1,N,p1a,p1b);

K2:=Stiffness_Matrices(X2);
G2:=k2*Local_Global(K2,N,p2a,p2b);

K3:=Stiffness_Matrices(X3);
G3:=k3*Local_Global(K3,N,p3a,p3b);

K4:=Stiffness_Matrices(X4);
G4:=k4*Local_Global(K4,N,p4a,p4b);

K5:=Stiffness_Matrices(X5);
G5:=k5*Local_Global(K5,N,p5a,p5b);

G:=G1+G2+G3+G4+G5;
Ginitial:=G;

//implementing boundary condition
for i in 1:N loop
 G[2*b1-1,i]:=0;
 G[2*b1,i]:=0;
 G[2*b2-1,i]:=0;
 G[2*b2,i]:=0;
end for;

G[2*b1-1,2*b1-1]:=1;
G[2*b1,2*b1]:=1;
G[2*b2-1,2*b2-1]:=1;
G[2*b2,2*b2]:=1;

//solving displacement
Sol:=Gauss_Jordan(N,G,X);

//solving reaction force
SolMat:=matrix(Sol);
XMat:=matrix(X);
R:=Reaction_Trusses(N,Ginitial,SolMat,XMat);

end Trusses_Tugas3_MRayhanFatur;

Fungsi Panggil Yang Digunakan

Matrice Transformation

function Stiffness_Matrices
input Real A;
Real Y;
output Real X[4,4];
Real float_error = 10e-10;

final constant Real pi=2*Modelica.Math.asin(1.0);

algorithm

Y:=A/180*pi;
    
X:=[(Modelica.Math.cos(Y))^2,Modelica.Math.cos(Y)*Modelica.Math.sin(Y),-(Modelica.Math.cos(Y))^2,-Modelica.Math.cos(Y)*Modelica.Math.sin(Y);

Modelica.Math.cos(Y)*Modelica.Math.sin(Y),(Modelica.Math.sin(Y))^2,-Modelica.Math.cos(Y)*Modelica.Math.sin(Y),-(Modelica.Math.sin(Y))^2;

-(Modelica.Math.cos(Y))^2,-Modelica.Math.cos(Y)*Modelica.Math.sin(Y),(Modelica.Math.cos(Y))^2,Modelica.Math.cos(Y)*Modelica.Math.sin(Y);

-Modelica.Math.cos(Y)*Modelica.Math.sin(Y),-(Modelica.Math.sin(Y))^2,Modelica.Math.cos(Y)*Modelica.Math.sin(Y),(Modelica.Math.sin(Y))^2];

for i in 1:4 loop
 for j in 1:4 loop
   if abs(X[i,j]) <= float_error then
     X[i,j] := 0;
   end if;
 end for;
end for;

end Stiffness_Matrices;

Global Element Matrice

function Local_Global
input Real Y[4,4];
input Integer B;
input Integer p1;
input Integer p2;
output Real G[B,B];

algorithm

for i in 1:B loop
 for j in 1:B loop
     G[i,j]:=0;
 end for;
end for;

G[2*p1,2*p1]:=Y[2,2];
G[2*p1-1,2*p1-1]:=Y[1,1];
G[2*p1,2*p1-1]:=Y[2,1];
G[2*p1-1,2*p1]:=Y[1,2];

G[2*p2,2*p2]:=Y[4,4];
G[2*p2-1,2*p2-1]:=Y[3,3];
G[2*p2,2*p2-1]:=Y[4,3];
G[2*p2-1,2*p2]:=Y[3,4];

G[2*p2,2*p1]:=Y[4,2];
G[2*p2-1,2*p1-1]:=Y[3,1];
G[2*p2,2*p1-1]:=Y[4,1];
G[2*p2-1,2*p1]:=Y[3,2];

G[2*p1,2*p2]:=Y[2,4];
G[2*p1-1,2*p2-1]:=Y[1,3];
G[2*p1,2*p2-1]:=Y[2,3];
G[2*p1-1,2*p2]:=Y[1,4];

end Local_Global;

Reaction Matrice Equation

function Reaction_Trusses
input Integer N;
input Real A[N,N];
input Real B[N,1];
input Real C[N,1];
Real X[N,1];
output Real Sol[N];
Real float_error = 10e-10;

algorithm
X:=A*B-C;

for i in 1:N loop
 if abs(X[i,1]) <= float_error then
   X[i,1] := 0;
 end if;
end for;

for i in 1:N loop
 Sol[i]:=X[i,1];
end for;

end Reaction_Trusses;

Gauss Jordan

function Gauss_Jordan input Integer N; input Real A[N,N]; input Real B[N]; output Real X[N]; Real float_error = 10e-10;

algorithm X:=Modelica.Math.Matrices.solve(A,B);

for i in 1:N loop

 if abs(X[i]) <= float_error then
   X[i] := 0;
 end if;

end for;

end Gauss_Jordan;


Grafik

Grafik Displacement
Grafik Reaction Forces


KUIS

Soal m3.jpg

Flowchart

Flowcharttruss.jpg

Minggu Ke-4

Pada minggu ke-4 ini, perkuliahan terhalang oleh libur karena adanya Pilkada, sehingga Pak Dai memberikan tugas untuk mengerjakan soal di bawah ini:


1607498641008.jpg


Untuk memahami alur logika untuk mengerjakan soal di atas dengan metode numerik, saya mencoba membuat Class Diagram dan Flow Chart terkait soal berikut:


DisplacementDisplacement


Dan berikut penjabaran dari flow chart di atas


FBD8.jpg


Berikut adalah pemrograman OpenModellica saya


function StiffnessMatrixElement

 input Real [:,9] inisiasi_mat;
 output Real [size(inisiasi_mat,1),6,6] Ke_mat;
 protected
   Real cos_x;
   Real cos_y;
   Real cos_z;
   Real [6] StiffTrig;
   Real [6,6] StiffTrans;
   Real [size(inisiasi_mat,1)] k_vec;

algorithm

 k_vec := {(inisiasi_mat[i,7] * inisiasi_mat[i,8] / inisiasi_mat[i,9]) for i in 1:size(inisiasi_mat,1)};
 // Finding stiffness matrix of each element member
 for i in 1:size(inisiasi_mat,1) loop
 // Clearing the matrices
 StiffTrig := zeros(6);
 StiffTrans := zeros(6,6);
 
 // Converting degrees to radians
 cos_x := inisiasi_mat[i,4];
 cos_y := inisiasi_mat[i,5];
 cos_z := inisiasi_mat[i,6];
 // {cos^2, sin^2, sincos}
 StiffTrig := {(cos_x)^2,
               (cos_y)^2,
               (cos_z)^2,
               (cos_x*cos_y),
               (cos_x*cos_z),
               (cos_y*cos_z)};
 
 // Construct stiffness transformation matrix
 StiffTrans := [  StiffTrig[1],    StiffTrig[4],    StiffTrig[5], -1*StiffTrig[1], -1*StiffTrig[4], -1*StiffTrig[5];
                  StiffTrig[4],    StiffTrig[2],    StiffTrig[6], -1*StiffTrig[4], -1*StiffTrig[2], -1*StiffTrig[6];
                  StiffTrig[5],    StiffTrig[6],    StiffTrig[3], -1*StiffTrig[5], -1*StiffTrig[6], -1*StiffTrig[3];
               -1*StiffTrig[1], -1*StiffTrig[4], -1*StiffTrig[5],    StiffTrig[1],    StiffTrig[4],    StiffTrig[5];
               -1*StiffTrig[4], -1*StiffTrig[2], -1*StiffTrig[6],    StiffTrig[4],    StiffTrig[2],    StiffTrig[6];
               -1*StiffTrig[5], -1*StiffTrig[6], -1*StiffTrig[3],    StiffTrig[5],    StiffTrig[6],    StiffTrig[3]];
               
 // Multiply in stiffness constant of element, add final stiffness matrix to Ke_mat
 for m in 1:6 loop
   for n in 1:6 loop
     Ke_mat[i,m,n] := k_vec[i] * StiffTrans[m,n];
   end for;
 end for;

end for; end StiffnessMatrixElement;


function StiffnessMatrixGlobal

 input Integer [:,2] n;
 input Integer x;
 input Integer y;
 input Real [y,6,6] Ke_mat; 
 output Real [y,x,x] Kg_mat;
 

algorithm

 for i in 1:y loop
   for a in 1:x loop
     for b in 1:x loop
       Kg_mat[i,a,b]:=0;
     end for;
   end for;
  end for;
 for i in 1:y loop
   Kg_mat[i,3*n[i,1],3*n[i,1]]:=Ke_mat[i,3,3];
   Kg_mat[i,3*n[i,1],3*n[i,1]-1]:=Ke_mat[i,3,2];
   Kg_mat[i,3*n[i,1],3*n[i,1]-2]:=Ke_mat[i,3,1];
   Kg_mat[i,3*n[i,1]-1,3*n[i,1]]:=Ke_mat[i,2,3];
   Kg_mat[i,3*n[i,1]-1,3*n[i,1]-1]:=Ke_mat[i,2,2];
   Kg_mat[i,3*n[i,1]-1,3*n[i,1]-2]:=Ke_mat[i,2,1];
   Kg_mat[i,3*n[i,1]-2,3*n[i,1]]:=Ke_mat[i,1,3];
   Kg_mat[i,3*n[i,1]-2,3*n[i,1]-1]:=Ke_mat[i,1,2];
   Kg_mat[i,3*n[i,1]-2,3*n[i,1]-2]:=Ke_mat[i,1,1];
   Kg_mat[i,3*n[i,2],3*n[i,2]]:=Ke_mat[i,6,6];
   Kg_mat[i,3*n[i,2],3*n[i,2]-1]:=Ke_mat[i,6,5];
   Kg_mat[i,3*n[i,2],3*n[i,2]-2]:=Ke_mat[i,6,4];
   Kg_mat[i,3*n[i,2]-1,3*n[i,2]]:=Ke_mat[i,5,6];
   Kg_mat[i,3*n[i,2]-1,3*n[i,2]-1]:=Ke_mat[i,5,5];
   Kg_mat[i,3*n[i,2]-1,3*n[i,2]-2]:=Ke_mat[i,5,4];
   Kg_mat[i,3*n[i,2]-2,3*n[i,2]]:=Ke_mat[i,4,6];
   Kg_mat[i,3*n[i,2]-2,3*n[i,2]-1]:=Ke_mat[i,4,5];
   Kg_mat[i,3*n[i,2]-2,3*n[i,2]-2]:=Ke_mat[i,4,4];
   Kg_mat[i,3*n[i,2],3*n[i,1]]:=Ke_mat[i,6,3];
   Kg_mat[i,3*n[i,2],3*n[i,1]-1]:=Ke_mat[i,6,2];
   Kg_mat[i,3*n[i,2],3*n[i,1]-2]:=Ke_mat[i,6,1];
   Kg_mat[i,3*n[i,2]-1,3*n[i,1]]:=Ke_mat[i,5,3];
   Kg_mat[i,3*n[i,2]-1,3*n[i,1]-1]:=Ke_mat[i,5,2];
   Kg_mat[i,3*n[i,2]-1,3*n[i,1]-2]:=Ke_mat[i,5,1];
   Kg_mat[i,3*n[i,2]-2,3*n[i,1]]:=Ke_mat[i,4,3];
   Kg_mat[i,3*n[i,2]-2,3*n[i,1]-1]:=Ke_mat[i,4,2];
   Kg_mat[i,3*n[i,2]-2,3*n[i,1]-2]:=Ke_mat[i,4,1];
   Kg_mat[i,3*n[i,1],3*n[i,2]]:=Ke_mat[i,3,6];
   Kg_mat[i,3*n[i,1],3*n[i,2]-1]:=Ke_mat[i,3,5];
   Kg_mat[i,3*n[i,1],3*n[i,2]-2]:=Ke_mat[i,3,4];
   Kg_mat[i,3*n[i,1]-1,3*n[i,2]]:=Ke_mat[i,2,6];
   Kg_mat[i,3*n[i,1]-1,3*n[i,2]-1]:=Ke_mat[i,2,5];
   Kg_mat[i,3*n[i,1]-1,3*n[i,2]-2]:=Ke_mat[i,2,4];
   Kg_mat[i,3*n[i,1]-2,3*n[i,2]]:=Ke_mat[i,1,6];
   Kg_mat[i,3*n[i,1]-2,3*n[i,2]-1]:=Ke_mat[i,1,5];
   Kg_mat[i,3*n[i,1]-2,3*n[i,2]-2]:=Ke_mat[i,1,4];
 end for;

end StiffnessMatrixGlobal;

function GaussJordan

 input Integer x;
 input Real [x,x] KgB_met;
 input Real [x] load_met;
 output Real [x] U_met;
 
 protected
 Real float_error = 10e-10;

algorithm

 U_met:=Modelica.Math.Matrices.solve(KgB_met,load_met);
 for i in 1:x loop
   if abs(U_met[i]) <= float_error then
    U_met[i] := 0;
   end if;
 end for;

end GaussJordan;



function ReactionForce

 input Integer x;
 input Real [x,x] KgTot_met;
 input Real [x] U_met;
 input Real [x] load_met;
 output Real [x] R_met;

algorithm

 R_met := (KgTot_met*U_met)-load_met;

end ReactionForce;



class TugasWeek4

 parameter Real [:,9] inisiasi = [1, 1, 2, -0.8,    0, -0.6, 15e-4, 70e9, 2.5;
                                  2, 1, 3, -0.8, -0.6,    0, 15e-4, 70e9, 2.5;
                                  3, 1, 4, -0.8,    0,  0.6, 15e-4, 70e9, 2.5];
                                  
 parameter Integer [:,2] node = [1, 2;
                                 1, 3;
                                 1, 4];
                                 
 parameter Integer y = size(node,1);
 
 parameter Integer x = 3*(size(node_load,1));
 
 parameter Integer z = size(Boundary,1);
 
 parameter Integer [:] Boundary = {2,3,4};
                              
 parameter Real [:,4] node_load = [1, 0, -5000, 0;
                                   2, 0,     0, 0;
                                   3, 0,     0, 0;
                                   4, 0,     0, 0];
                                   
 parameter Real [x] load = {0,-5000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
                                  
 Real [y] k;
 Real [y,6,6] Ke;
 
 Real [y,x,x] Kg;
 
 Real [x,x] KgTot;
 
 Real [x,x] KgB;
 
 Real [x] U;
 
 Real [x] R;
 

equation

k = {(inisiasi[i,7] * inisiasi[i,8] / inisiasi[i,9]) for i in 1:y};
Ke = StiffnessMatrixElement(inisiasi);

Kg = StiffnessMatrixGlobal(node, x, y, Ke);

KgTot = SumStiffnessMatrixGlobal(x, y, Kg);

KgB = BoundaryStiffnessMatrixGlobal(x, z, KgTot, Boundary);

U = GaussJordan(x, KgB, load);

R = ReactionForce(x, KgTot, U, load);

end TugasWeek4;


Kuliah Minggu Ke-7 (23 Desember 2020)

TUJUAN

Kuliah kali ini kita melaksanakan kuis mengenai optimasi 3D Truss menggunakan class algoritma yang sebelumnya telah disusun pada kuliah minggu lalu. 3D Truss yang dimaksud seperti ini.


TubesMetnumRayhan.jpg
Soaltubes1.jpg


Kuis ini juga sekaligus mendukung pengerjaan tugas besar, pada tugas besar Pak Dai meminta untuk menentukan jenis material, harga material, dan luas cross sectional area dari truss. Kami juga diminta untuk melakukan optimasi dan membentuk curve fitting pada variabel harga

PENGERJAAN

Diketahui: panjang= 0,6 m lebar= 0,75 m tinggi= 1,6 m F1= 2000 N F2= 1000 N

kemudian diasumsikan

-Beban akan terdistribusi hanya pada node (karena bersifat trusses).

-Safety factor minimal bernilai 2.

-Batas displacement 0,001 m sebelum buckling (pada truss paling atas).

-Ketinggian trusses pada tiap lantai sama yaitu 0,6 m.


berikutnya: Untuk jenis material yang sama:

-Mencari harga untuk 5 ukuran batang dengan material yang sama yaitu SS304.

-Mengitung nilai safety factor pada 5 ukuran batang dengan coding awal.

-Membuat rasio antara safety factor dengan harga total.

-Membuat persamaan antara rasio (safe/harga) dengan area dengan metode curve fitting.

-Melakukan optimasi menggunakan metode golden section.


Untuk area penampang yang sama:

-Mencari harga untuk 3 jenis material dengan area penampang yang sama yaitu 171 mm^2.

-Mengitung nilai safety factor pada 3 jenis batang dengan coding awal.

-Membuat rasio antara safety factor dengan harga total.

-Membuat persamaan antara rasio (safe/harga) dengan jenis material dengan metode curve fitting.

-Melakukan optimasi menggunakan metode golden section.


Update Tugas Besar

model TugasBesar_Trusses_3d_Safety //define initial variable parameter Integer Points=size(P,1); //Number of Points parameter Integer Trusses=size(C,1); //Number of Trusses parameter Real Yield=215e6; //Yield Strength (Pa) parameter Real Area=0.000171; //Area L Profile (Dimension=0.03, Thickness=0,003) (m2) parameter Real Elas=193e9; //Elasticity SS 304 (Pa) //define connection parameter Integer C[:,2]=[1,5;

                        2,6;
                        3,7;
                        4,8;
                        5,6;  //1st floor
                        6,7;  //1st floor
                        7,8;  //1st floor
                        5,8;  //1st floor
                        5,9;
                       6,10;
                       7,11;
                       8,12;
                       9,10; //2nd floor
                       10,11;//2nd floor 
                       11,12;//2nd floor
                        9,12; //2nd floor
                        9,13;
                       10,14;
                       11,15;
                       12,16;
                       13,14;//3rd floor
                       14,15;//3rd floor
                       15,16;//3rd floor
                      13,16];//3rd floor                                                           

//define coordinates (please put orderly) parameter Real P[:,6]=[0.3,-0.375,0,1,1,1; //1

                     -0.3,-0.375,0,1,1,1;    //2
                     -0.3,0.375,0,1,1,1;     //3
                     0.3,0.375,0,1,1,1;      //4
                         
                     0.3,-0.375,0.6,0,0,0;   //5
                     -0.3,-0.375,0.6,0,0,0;  //6
                     -0.3,0.375,0.6,0,0,0;   //7
                     0.3,0.375,0.6,0,0,0;    //8
                          
                     0.3,-0.375,1.2,0,0,0;   //9
                     -0.3,-0.375,1.2,0,0,0;  //10  
                     -0.3,0.375,1.2,0,0,0;   //11
                     0.3,0.375,1.2,0,0,0;    //12
                          
                     0.3,-0.375,1.8,0,0,0;   //13
                     -0.3,-0.375,1.8,0,0,0;  //14
                     -0.3,0.375,1.8,0,0,0;   //15
                     0.3,0.375,1.8,0,0,0];   //16                         

//define external force (please put orderly) parameter Real F[Points*3]={0,0,0,

                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,-500, 
                          0,0,-1000, 
                          0,0,-1000, 
                          0,0,-500}; 

//solution Real displacement[N], reaction[N]; Real check[3]; Real stress1[Trusses]; Real safety[Trusses]; Real dis[3]; Real Str[3]; protected parameter Integer N=3*Points; Real q1[3], q2[3], g[N,N], G[N,N], G_star[N,N], id[N,N]=identity(N), cx, cy, cz, L, X[3,3]; Real err=10e-10, ers=10e-4; algorithm //Creating Global Matrix G:=id; for i in 1:Trusses loop for j in 1:3 loop

q1[j]:=P[C[i,1],j];
q2[j]:=P[C[i,2],j];

end for;

 //Solving Matrix
 L:=Modelica.Math.Vectors.length(q2-q1);
 cx:=(q2[1]-q1[1])/L;
 cy:=(q2[2]-q1[2])/L;
 cz:=(q2[3]-q1[3])/L; 
 X:=(Area*Elas/L)*[cx^2,cx*cy,cx*cz;
                   cy*cx,cy^2,cy*cz;
                   cz*cx,cz*cy,cz^2];
 //Transforming to global matrix
 g:=zeros(N,N); 
 for m,n in 1:3 loop
   g[3*(C[i,1]-1)+m,3*(C[i,1]-1)+n]:=X[m,n];
   g[3*(C[i,2]-1)+m,3*(C[i,2]-1)+n]:=X[m,n];
   g[3*(C[i,2]-1)+m,3*(C[i,1]-1)+n]:=-X[m,n];
   g[3*(C[i,1]-1)+m,3*(C[i,2]-1)+n]:=-X[m,n];
 end for;  

G_star:=G+g; G:=G_star; end for; //Implementing boundary for x in 1:Points loop if P[x,4] <> 0 then

 for a in 1:Points*3 loop
   G[(x*3)-2,a]:=0;
   G[(x*3)-2,(x*3)-2]:=1;
 end for;

end if; if P[x,5] <> 0 then

 for a in 1:Points*3 loop
   G[(x*3)-1,a]:=0;
   G[(x*3)-1,(x*3)-1]:=1;
 end for;

end if; if P[x,6] <> 0 then

 for a in 1:Points*3 loop
   G[x*3,a]:=0;
   G[x*3,x*3]:=1;
 end for;

end if; end for; //Solving displacement displacement:=Modelica.Math.Matrices.solve(G,F); //Solving reaction reaction:=(G_star*displacement)-F; //Eliminating float error for i in 1:N loop reaction[i]:=if abs(reaction[i])<=err then 0 else reaction[i]; displacement[i]:=if abs(displacement[i])<=err then 0 else displacement[i]; end for; //Checking Force check[1]:=sum({reaction[i] for i in (1:3:(N-2))})+sum({F[i] for i in (1:3:(N-2))}); check[2]:=sum({reaction[i] for i in (2:3:(N-1))})+sum({F[i] for i in (2:3:(N-1))}); check[3]:=sum({reaction[i] for i in (3:3:N)})+sum({F[i] for i in (3:3:N)}); for i in 1:3 loop check[i] := if abs(check[i])<=ers then 0 else check[i]; end for; //Calculating stress in each truss for i in 1:Trusses loop for j in 1:3 loop

q1[j]:=P[C[i,1],j];
q2[j]:=P[C[i,2],j];
dis[j]:=abs(displacement[3*(C[i,1]-1)+j]-displacement[3*(C[i,2]-1)+j]);

end for;

 //Solving Matrix
 L:=Modelica.Math.Vectors.length(q2-q1);
 cx:=(q2[1]-q1[1])/L;
 cy:=(q2[2]-q1[2])/L;
 cz:=(q2[3]-q1[3])/L; 
 X:=(Elas/L)*[cx^2,cx*cy,cx*cz;
              cy*cx,cy^2,cy*cz;
              cz*cx,cz*cy,cz^2];   
 Str:=(X*dis); 
 stress1[i]:=Modelica.Math.Vectors.length(Str);

end for; //Safety factor for i in 1:Trusses loop if stress1[i]>0 then

 safety[i]:=Yield/stress1[i];

else

 safety[i]:=0;

end if; end for; end TugasBesar_Trusses_3d_Safety;


dan didapati hasil berikut dari simulasi tersebut