Difference between revisions of "Metnum03 - Elita Kabayeva"

From ccitonlinewiki
Jump to: navigation, search
(→‎PENYELESAIAN MANUAL KUIS)
 
(6 intermediate revisions by the same user not shown)
Line 20: Line 20:
  
 
Pada open methods diajarkan menggunakan Newton Rhapson, Simple Fix Point dan Secant Method. Sedangkan pada metode bracketing diajarkan menggunakan False methods dan Bisections.
 
Pada open methods diajarkan menggunakan Newton Rhapson, Simple Fix Point dan Secant Method. Sedangkan pada metode bracketing diajarkan menggunakan False methods dan Bisections.
 +
 +
NEWTON RHAPSON
 +
 +
Metode newton-raphson dapat diwakili dengan formula berikut :
 +
 +
[[File:NR_formula_ElitaK.png|400px|center]]
 +
 +
 +
BISECTION METHOD
 +
 +
Metode bisection merupakan salah satu metode incremental search yang mana interval dari dua titik x dibagi dua sehingga mendapatkan nilai x lagi. Berikut adalah langkah-langkah metode bisection :
 +
 +
• Lakukan tembakan pada dua titik terendah (xi) dan titik tertinggi (xu). • Kemudian carilah xr dengan menjumlahkan antara xi dan xu kemudian hasil dari penjumlahan tersebut dibagi dua. • Lakukan evaluasi sebagai berikut :
 +
 +
Jika f(xl) f(xr) < 0, akar persamaan terletak di sub interval terendah. Jadi, atur xu = xr kemudian kembali ke langkah 2.
 +
 +
Jika f(xl) f(xr) > 0, akar persamaan terletak di sub interval tertinggi. Jadi, atur xi = xr kemudian kembali ke langkah 2.
 +
 +
Jika f(xl) f(xr) = 0, akar persamaan sama dengan xr. Sehingga hentikan perhitungan.
 +
 +
  
 
3. Penurunan Numeric Pada penurunan ini kita diajarkan penurunan secara backward, forward dan center.
 
3. Penurunan Numeric Pada penurunan ini kita diajarkan penurunan secara backward, forward dan center.
Line 571: Line 592:
  
 
|}
 
|}
 +
 +
 +
 +
== TUGAS PERTEMUAN V ==
 +
 +
Untuk pertemuan ini kami ditugaskan mengaplikasikan coding 3D Truss milik Ahmad Mohammad Fahmi ke permasalahan 3D Truss yang lainnya.
 +
 +
[[File:Tugas_3D_ElitaK.jpeg|600px|thumb|right]]
 +
 +
{| class="wikitable"
 +
|-
 +
| style='border-style: none  none  solid  solid;' |
 +
 +
''Calculate Length and Cosine''
 +
class Length_Cos
 +
 +
parameter Integer Points=4;
 +
parameter Integer Trusses=6;
 +
parameter Integer C[Trusses,2]=[1,2;
 +
                                1,3;
 +
                                1,4;
 +
                                2,3;
 +
                                3,4;
 +
                                2,4];
 +
parameter Integer P[:,3]=[0,0,3;
 +
                          6,0,0;
 +
                          0,0,-3;
 +
                          0,6,0];
 +
 +
//Element 1,4,5
 +
Real L1 [size(P,1)-1];
 +
  Real cosx1 [size(P,1)-1];
 +
  Real cosy1 [size(P,1)-1];
 +
    Real cosz1 [size(P,1)-1];
 +
 +
//Element 2 6
 +
Real L2 [size(P,1)-2];
 +
  Real cosx2 [size(P,1)-2];
 +
  Real cosy2 [size(P,1)-2];
 +
    Real cosz2 [size(P,1)-2];
 +
 +
//Element 3
 +
Real L3 [size(P,1)-3];
 +
  Real cosx3 [size(P,1)-3];
 +
  Real cosy3 [size(P,1)-3];
 +
    Real cosz3 [size(P,1)-3];
 +
   
 +
equation
 +
  //Element 1 4 5
 +
  for i in 1:1:(size(P,1)-1) loop
 +
    L1[i] = sqrt((P[i+1,1]-P[i,1]^2+(P[i+1,2]-P[i,2])^2+(P[i+1,3]-P[i,3])^2));
 +
  end for;
 +
  for i in 1:(size(P,1)-1) loop
 +
    cosx1 [i] = (P[i+1,1]-P[i,1])/L1[i];
 +
    cosy1 [i] = (P[i+1,2]-P[i,2])/L1[i];
 +
    cosz1 [i] = (P[i+1,3]-P[i,3])/L1[i];
 +
  end for;
 +
 
 +
  //Element 2 6
 +
  for i in 1:1:(size(P,1)-2) loop
 +
    L2 [i] = sqrt((P[i+2,1]-P[i,1]^2+(P[i+2,2]-P[i,2])^2+(P[i+2,3]-P[i,3])^2));
 +
  end for;
 +
  for i in i:(size(P,1)-2) loop
 +
    cosx2 [i] = (P[i+2,1]-P[i,1])/L2[i];
 +
    cosy2 [i] = (P[i+2,2]-P[i,2])/L2[i];
 +
    cosz2 [i] = (P[i+2,3]-P[i,3])/L2[i];
 +
  end for;
 +
 
 +
  //Element 3
 +
  for i in 1:1:(size(P,1)-3) loop
 +
    L3 [i] = sqrt((P[i+3,1]-P[i,1]^2+(P[i+3,2]-P[i,2])^2+(P[i+3,3]-P[i,3])^2));
 +
  end for;
 +
  for i in i:(size(P,1)-3) loop
 +
    cosx3 [i] = (P[i+3,1]-P[i,1])/L3[i];
 +
    cosy3 [i] = (P[i+3,2]-P[i,2])/L3[i];
 +
    cosz3 [i] = (P[i+3,3]-P[i,3])/L3[i];
 +
  end for;
 +
end Length_Cos;
 +
 +
|}
 +
 +
 +
'''MASTER CODE (CLASS)'''
 +
 +
{| class="wikitable"
 +
|-
 +
| style='border-style: none  none  solid  solid;'
 +
class Soal5
 +
  //inisiasi = [ elemen#, dX, dY, dZ, A, E]
 +
  parameter Real [:,6] inisiasi = [1,  0.8944,  0.0000, -0.4472, 1.56, 10.6e6; //isi sesuai data
 +
                                  2,  0.0000,  0.0000, -1.0000, 1.56, 10.6e6;
 +
                                  3,  0.0000,  0.8944, -0.4472, 1.56, 10.6e6;
 +
                                  4, -0.8944,  0.0000, -0.4472, 1.56, 10.6e6;
 +
                                  5, -0.7071,  0.7071,  0.0000, 1.56, 10.6e6;
 +
                                  6,  0.0000,  0.8944,  0.4472, 1.56, 10.6e6];
 +
 +
  //node = [ i, j]                               
 +
  parameter Integer [size(inisiasi,1),2] node = [1, 2; //isi sesuai data
 +
                                                1, 3;
 +
                                                1, 4;
 +
                                                2, 3;
 +
                                                2, 4;
 +
                                                3, 4];
 +
 +
  //jumlah node
 +
  parameter Integer n = 4; //isi sesuai data
 +
 +
  //titik node boundary xyz
 +
  parameter Integer [:] Boundary_xyz = {1,3,4}; //isi sesuai data
 +
 +
  //titik node boundary xy
 +
  parameter Integer [:] Boundary_xy = {3}; //isi sesuai data
 +
 +
  //titik node boundary xz
 +
  parameter Integer [:] Boundary_xz = {0}; //isi sesuai data
 +
 +
  //titik node boundary yz
 +
  parameter Integer [:] Boundary_yz = {0}; //isi sesuai data
 +
 +
  //titik node boundary x
 +
  parameter Integer [:] Boundary_x = {2}; //isi sesuai data
 +
 +
  //titik node boundary y
 +
  parameter Integer [:] Boundary_y = {4}; //isi sesuai data
 +
 +
  //titik node boundary z
 +
  parameter Integer [:] Boundary_z = {1}; //isi sesuai data
 +
                           
 +
  //load = [ F1x, F1y, F1z,..., Fnx, Fny, Fnz]
 +
  parameter Real [3*n] load = {0,    0, 0,  //isi sesuai data
 +
                              0, -200, 0,
 +
                              0,    0, 0,
 +
                              0,    0, 0};
 +
 +
  Real [size(inisiasi,1)] L;
 +
  Real [size(inisiasi,1)] k;
 +
  Real [size(inisiasi,1),6,6] Ke;
 +
  Real [size(inisiasi,1),3*n,3*n] Kg;
 +
  Real [3*n,3*n] KgTot;
 +
  Real [3*n,3*n] KgB;
 +
  Real [3*n] U;
 +
  Real [3*n] R;
 +
 +
  //check force
 +
  Real [3] F;
 +
 +
equation
 +
L = {(sqrt(inisiasi[i,2]^2 + inisiasi[i,3]^2 + inisiasi[i,4]^2)) for i in 1:size(inisiasi,1)};
 +
 +
k = {(inisiasi[i,5] * inisiasi[i,6] / L[i]) for i in 1:size(inisiasi,1)};
 +
 +
Ke = StiffnessMatrixElement(inisiasi);
 +
 +
Kg = StiffnessMatrixGlobal(n, node, Ke);
 +
 +
KgTot = SumStiffnessMatrixGlobal(Kg);
 +
 +
KgB = BoundaryStiffnessMatrixGlobal(KgTot, Boundary_xyz, Boundary_xy, Boundary_xz, Boundary_yz, Boundary_x, Boundary_y, Boundary_z);
 +
 +
U = GaussJordan(KgB, load);
 +
 +
R = ReactionForce(KgTot, U, load);
 +
 +
F = CheckForce(load,R);
 +
 +
 +
end Soal5;
 +
 +
|}
 +
 +
Saya mengikuti langkah-langkah coding yang diajarkan oleh Ahmad Mohammad Fahmi dan saat telah melakukan verifikasi code, seluruh code dinyatakan OK pada OpenModelica. Namun pada saat simulate, dalam proses runningnya, coding mengalami crash sehingga tidak dapat di plotting. Saya masih meneliti coding untuk melakukan revisi agar tidak terjadi crash.
 +
 +
 +
== OPTIMASI ==
 +
 +
Pada kesempatan ini kami mempelajari mengenai optimasi dan aplikasi metoda numeriknya. Berikut adalah coding pada OpenModelica yang kami pelajari;
 +
 +
 +
'''FUNCTION'''
 +
 +
{| class="wikitable"
 +
|-
 +
| style='border-style: none  none  solid  solid;' |
 +
 +
function Func_Optimization
 +
input Real x;
 +
import Modelica.Math;
 +
output Real y;
 +
 +
algorithm
 +
y:=2*sin(x)-x^2/10;
 +
 +
end Func_Optimization;
 +
 +
|}
 +
 +
 +
'''OPTIMIZATION GOLDEN SECTION'''
 +
 +
{| class="wikitable"
 +
|-
 +
| style='border-style: none  none  solid  solid;' |
 +
 +
model Opt_Gold
 +
 +
parameter Real xlo=0;
 +
parameter Real xhi=4;
 +
parameter Integer N=8; // maximum iteration
 +
parameter Real es=0.0001; // maximum error
 +
 +
Real f1[N], f2[N], x1[N], x2[N], ea[N];
 +
Real xopt,  fx;
 +
protected
 +
Real d, xl, xu, xint, R=(5^(1/2)-1)/2;
 +
 +
algorithm
 +
xl := xlo;
 +
xu := xhi;
 +
 +
for i in 1:N loop
 +
d:= R*(xu-xl);
 +
x1[i]:=xl+d;
 +
x2[i]:=xu-d;
 +
f1[i]:=Func_Optimization(x1[i]);
 +
f2[i]:=Func_Optimization(x2[i]);
 +
xint:=xu-xl;
 +
 +
if f1[i]>f2[i] then
 +
  xl:=x2[i];
 +
  xopt:=x1[i];
 +
  fx:=f1[i];
 +
  else
 +
    xu:=x1[i];
 +
    xopt:=x2[i];
 +
    fx:=f2[i];
 +
end if;
 +
 +
ea[i]:=(1-R)*abs((xint)/xopt);
 +
if ea[i]<es then
 +
  break;
 +
end if;
 +
end for;
 +
 +
end Opt_Gold;
 +
 +
|}
 +
 +
 +
'''NEWTON OPTIMIZATION'''
 +
 +
{| class="wikitable"
 +
|-
 +
| style='border-style: none  none  solid  solid;' |
 +
 +
model Opt_Newton
 +
 +
parameter Real g=2.5; //initial guess
 +
parameter Integer N=8; // maximum iteration
 +
parameter Real es=0.0001; // maximum error
 +
 +
Real X[N];
 +
Real xopt, fx, ea[N];
 +
 +
algorithm
 +
X[1]:=g;
 +
ea[1]:=1;
 +
 +
for i in 2:N loop
 +
X[i]:=X[i-1]-Func_Optimization_Der(X[i-1])/Func_Optimization_Der_Der(X[i-1]);
 +
xopt:=X[i];
 +
 +
ea[i]:=abs(1-X[i-1]/X[i]);
 +
if ea[i]<es then
 +
  break;
 +
end if;
 +
end for;
 +
fx:=Func_Optimization(xopt);
 +
 +
end Opt_Newton;
 +
 +
|}
 +
 +
 +
 +
'''PARABOLIC OPTIMIZATION'''
 +
 +
{| class="wikitable"
 +
|-
 +
| style='border-style: none  none  solid  solid;' |
 +
 +
model Opt_Parabolic
 +
 +
parameter Real g1=0; //initial guess
 +
parameter Real g2=1; //initial guess
 +
parameter Real g3=4; //initial guess
 +
parameter Integer N=5; // maximum iteration
 +
parameter Real es=0.0001; // maximum error
 +
 +
Real x1, x2, x3, xopt, xp[N], ea[N];
 +
//Real xl, xm, xu;
 +
Real fx1, fx2, fx3, fx, A[4], A_star[4];
 +
algorithm
 +
x1:=g1;
 +
x2:=g2;
 +
x3:=g3;
 +
 +
for i in 1:N loop
 +
fx1:=Func_Optimization(x1);
 +
fx2:=Func_Optimization(x2);
 +
fx3:=Func_Optimization(x3);
 +
 +
xp[i]:=(fx1*(x2^2-x3^2)+fx2*(x3^2-x1^2)+fx3*(x1^2-x2^2))/(2*fx1*(x2-x3)+2*fx2*(x3-x1)+2*fx3*(x1-x2));
 +
xopt:=xp[i];
 +
fx:=Func_Optimization(xp[i]);
 +
A:={x1,x2,x3,xp[i]};
 +
A_star:=Modelica.Math.Vectors.sort(A);
 +
 
 +
if xp[i]>x2 then
 +
  x1:=A_star[2];
 +
  x2:=A_star[3];
 +
  x3:=A_star[4];
 +
  else
 +
    x1:=A_star[1];
 +
    x2:=A_star[2];
 +
    x3:=A_star[3];
 +
end if;
 +
end for;
 +
 +
ea[1]:=1;
 +
for i in 2:N loop 
 +
ea[i]:=abs(1-xp[i-1]/xp[i]);
 +
if ea[i]<es then
 +
  break;
 +
end if;
 +
end for;
 +
 +
end Opt_Parabolic;
 +
|}
 +
 +
 +
 +
== TUGAS BESAR METODE NUMERIK ==
 +
 +
 +
'''OBJEKTIF'''
 +
 +
Tujuan dari tugas besar ini adalah>> optimasi harga pembuatan rangka truss sederhana dengan memvariasi dimensi dan elastisitas material.
 +
 +
 +
 +
'''GEOMETRI'''
 +
 +
[[File:Tubes_ElitaKjpg.jpeg|600px|thumb|center|Geometri dan Load Truss]]
 +
 +
 +
 +
'''ASUMSI'''
 +
 +
 +
Dalam tugas besar ini digunakan asumsi-asumsi sebagai berikut:
 +
 +
1. Beban hanya terdistribusi pada node.
 +
 +
2. ''Safety Factor'' minimal 2.
 +
 +
3. Batas ''displacement'' pada truss paling atas sebelum ''buckling'' bernilai 0,001 m.
 +
 +
4. Ketinggian ''trusses'' pada setiap lantai konstan, sebesar 0,6 m.
 +
 +
 +
'''METODOLOGI'''
 +
 +
 +
'''I. ELASTICITY LOCKED'''
 +
 +
1. Mencari data yield strength dan elasticity dari material yang dipilih (disini digunakan SS400/ASTM A36)
 +
 +
2. Mencari harga truss sesuai dengan material.
 +
 +
3. Menghitung nilai safety factor.
 +
 +
4. Membuat rasio (efficiency) dari safety factor dengan total cost.
 +
 +
5. Membuat persamaan (mencari Coe) antara rasio dengan area menggunakan curve-fitting.
 +
 +
6. Melakukan optimasi menggunakan metode golden section.
 +
 +
 +
'''II. AREA LOCKED'''
 +
 +
1. Mendefinisikan luas penampang yang sama untuk seluruh variasi material >> 3,75 x 10^-4 m^2.
 +
 +
2. Mencari harga untuk 3 jenis variasi material >> SS400, SS304, SS316.
 +
 +
3. Menghitung nilai ''safety factor'' pada 3 variasi truss.
 +
 +
4. Membuat rasio antara ''safety factor'' dengan ''total cost''.
 +
 +
5. Membuat persamaan (mencari Coe) antara rasio dengan area menggunakan curve-fitting.
 +
 +
6. Melakukan optimasi menggunakan metode golden section.
 +
 +
 +
'''KALKULASI DISPLACEMENT, REACTION FORCE, & STRESS VIA OPENMODELICA'''
 +
 +
{| class="wikitable"
 +
|-
 +
| style='border-style: none  none  solid  solid;' |
 +
 +
model Trusses_3D_Tugas_Besar_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.000224;  //Area L Profile (Dimension=0.03, Thickness=0,004) (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 Trusses_3D_Tugas_Besar_Safety;
 +
 +
|}
 +
 +
 +
'''KALKULASI COE UNTUK CURVE-FITTING VIA OPENMODELICA'''
 +
 +
{| class="wikitable"
 +
|-
 +
| style='border-style: none  none  solid  solid;' |
 +
 +
function Curve_Fitting
 +
 +
 +
input Real X[:];
 +
 +
input Real Y[size(X,1)];
 +
 +
input Integer order=2;
 +
 +
output Real Coe[order+1];
 +
 +
 +
protected
 +
 +
Real Z[size(X,1),order+1];
 +
 +
Real ZTr[order+1,size(X,1)];
 +
 +
Real A[order+1,order+1];
 +
 +
Real B[order+1];
 +
 +
 +
algorithm
 +
 +
for i in 1:size(X,1) loop
 +
  for j in 1:(order+1) loop
 +
  Z[i,j]:=X[i]^(order+1-j);
 +
  end for;
 +
end for;
 +
ZTr:=transpose(Z);
 +
 +
A:=ZTr*Z;
 +
 +
B:=ZTr*Y;
 +
 +
Coe:=Modelica.Math.Matrices.solve(A,B);
 +
 +
//Coe:=fill(2,size(Coe,1));
 +
 +
end Curve_Fitting;
 +
/*
 +
for i in 1:3 loop
 +
for j in 1:Points loop
 +
  R[j]:=reaction[3*(j-1)+i];
 +
end for;
 +
Sur[i]:=sum(R);
 +
end for;
 +
*/
 +
 +
| style="width: 10cm;"|
 +
 +
 +
model CurveFitting
 +
 +
parameter Real X[6]={141e-6,375e-6,384e-6,575e-6,931e-6,864e-6};//area
 +
 +
parameter Real Y[6]={7702,7833,7664,7708,7927,7805};//Cost/Kg
 +
 +
Real Coe[3];
 +
 +
algorithm
 +
 +
Coe:=Curve_Fitting(X,Y,2);
 +
 +
end CurveFitting;
 +
 +
|}
 +
 +
 +
 +
'''TABULASI DATA DI EXCEL & SNAPSHOT Coe'''
 +
 +
 +
'''I. ELASTICITY LOCKED EXCEL TABULATION'''
 +
 +
 +
[[File:Elasticity_locked_ElitaK.png|1000px|center]]
 +
 +
 +
 +
'''Coe for Curve-Fitting'''
 +
 +
[[File:Coef_elasticity_ElitaK.png|250px|center]]
 +
 +
 +
 +
'''II. AREA LOCKED EXCEL TABULATION'''
 +
 +
 +
[[File:Area_locked_ElitaK.png|1000px|center]]
 +
 +
 +
 +
'''Coe for Curve-Fitting'''
 +
 +
 +
[[File:Coef_material_1CostKg_ElitaK.png|250px|thumb|left|Curve-Fitting for Cost/Kg]]
 +
 +
 +
 +
[[File:Coef_material_2_density_ElitaK.png|250px|thumb|center|Curve-Fitting for Density]]
 +
 +
 +
 +
[[File:Coef_material_3_yield_ElitaK.png|250px|thumb|left|Curve-Fitting for Yield Strength]]
 +
 +
 +
 +
 +
== UAS METODE NUMERIK (13/01/2021) ==
 +
 +
 +
Terlampir di bawah ini adalah jawaban untuk Ujian Akhir Semester Metode Numerik saya pada hari ini Rabu, 13 Januari 2021.
 +
 +
 +
 +
 +
[[File:UAS_1_ElitaK.jpeg|600px|center]]
 +
 +
 +
[[File:UAS_2_ElitaK.jpeg|600px|center]]
 +
 +
 +
[[File:UAS_3_ElitaK.jpeg|600px|center]]
 +
 +
 +
[[File:UAS_45_ElitaK.jpg|600px|center]]
 +
 +
 +
[[File:UAS_67_ElitaK.jpg|600px|center]]
 +
 +
 +
[[File:UAS_ss1_ElitaK.png|600px|center]]
 +
 +
 +
[[File:UAS_ss2_ElitaK.png|600px|center]]

Latest revision as of 21:17, 14 January 2021

BIODATA

Nama : Elita Kabayeva

NPM : 1906435486

Pendidikan Terakhir : D-III


PERTEMUAN I (9/11/2020)

Dalam pertemuan pertama ini, pak Ahmad Indra melakukan review mengenai materi-materi yang telah dipelajari sebelum UTS. Sekaligus memberikan tugas untuk membuat video mengenai pembelajaran metode numerik yang saya lakukan menggunakan software OPENMODELICA. Materi yang diajarkan sebelum UTS adalah sebagai berikut Materi yang diajarkan sebelum UTS adalah sebagai berikut

1. Deret Taylor dan Deret Mclaurin

2. Pencarian Akar (Open Methods dan Bracketing)

Pada open methods diajarkan menggunakan Newton Rhapson, Simple Fix Point dan Secant Method. Sedangkan pada metode bracketing diajarkan menggunakan False methods dan Bisections.

NEWTON RHAPSON

Metode newton-raphson dapat diwakili dengan formula berikut :

NR formula ElitaK.png


BISECTION METHOD

Metode bisection merupakan salah satu metode incremental search yang mana interval dari dua titik x dibagi dua sehingga mendapatkan nilai x lagi. Berikut adalah langkah-langkah metode bisection :

• Lakukan tembakan pada dua titik terendah (xi) dan titik tertinggi (xu). • Kemudian carilah xr dengan menjumlahkan antara xi dan xu kemudian hasil dari penjumlahan tersebut dibagi dua. • Lakukan evaluasi sebagai berikut :

Jika f(xl) f(xr) < 0, akar persamaan terletak di sub interval terendah. Jadi, atur xu = xr kemudian kembali ke langkah 2.

Jika f(xl) f(xr) > 0, akar persamaan terletak di sub interval tertinggi. Jadi, atur xi = xr kemudian kembali ke langkah 2.

Jika f(xl) f(xr) = 0, akar persamaan sama dengan xr. Sehingga hentikan perhitungan.


3. Penurunan Numeric Pada penurunan ini kita diajarkan penurunan secara backward, forward dan center.

4. Regresi Linear dan Interpolasi

5. Gauss-Jordan dan Aljabar

Untuk video pembelajaran tersebut saya upload pada link Youtube berikut.


OPENMODELICA (PENYELESAIAN SISTEM PERSAMAAN PANGKAT TIGA, TIGA VARIABEL

Video diatas adalah perbaikan dari video sebelumnya yang telah saya upload di Youtube. Dikarenakan, video yang sebelumnya terdapat error pada audionya.


PERTEMUAN II (16/11/2020)

Pada pertemuan kedua ini, saya izin tidak mengikuti kelas sinkron melalui Zoom dikarenakan sakit. Namun, saya tetap melakukan pembelajaran mandiri setelahnya dengan membaca rangkuman dari teman-teman sekelas serta mencoba melakukan simulasi pencarian mean dengan menggunakan coding pada OPENMODELICA, sesuai dengan materi yang dilaksanakan pada Zoom di pertemuan kedua.

Coding OPENMODELICA

Diatas adalah coding OPENMODELICA untuk mencari mean. Disini saya menggunakan 25 data angka riil.

Setelah melakukan coding, saya mengecek apakah coding tersebut benar atau tidak.

Hasil Check Coding

Baru kemudian melakukan simulasi.


Proses Simulasi


Berikut adalah hasil dari simulasi yang menyatakan nilai mean dari 25 data yang di input kan.

Hasil Mean


TUGAS PERTEMUAN II

Dalam pertemuan II kemarin, pak Dai memberi tugas untuk menyelesaikan persamaan aljabar simultan dengan menggunakan Gauss-Elimination, Gauss Seidel ataupun metoda lain.

Saya akan menggunakan metoda Gauss-Elimination untuk menyelesaikan persamaan aljabar berikut:

x1 + 2x2 + 3x3 + x4 = 9

3x1 + 5x2 + 7x3 + 4x4 = 12

4x1 + x2 + x3 + 3x4 = 23

6x1 + 7x2 + 5x3 + 2x4 = 0


Berikut adalah hasil input coding saya pada OpenModelica untuk membentuk matriks dari persamaan-persamaan tersebut. Untuk coding, saya menggunakan perintah yang ada di library modelica yaitu "Modelica.Math.Matrices.solve(A,b)" untuk menyelesaikan sistem persamaan linier eliminasi gauss yang ada diatas.

Coding


Kemudian, setelah melakukan input, saya melakukan pengecekan coding.


Coding Check


Setelah dilakukan simulasi, berikut adalah hasil untuk penyelesaian Eliminasi Gauss terhadap persamaan-persamaan diatas.


Hasil


Hasilnya :

X1 = 11.8824

X2 = -17.5294

X3 = 12.9412

X4 = -6.64706


PERTEMUAN III (23/11/2020)

Pada pertemuan hari ini, pak Dai memberi tugas untuk menyelesaikan persoalan pegas yang terdapat pada bab 12 dari buku Numerical Method (Steven Chapras) dan membuktikan secara matematis hasil dari penyelesaian tersebut.

Saya melakukan perhitungan matematis terlebih dahulu dengan pertama-tama membuat Free Body Diagram dari masing-masing spring mass system yang ada di buku dan kemudian menuliskan persamaannya.


Pembuktian Matematis 1


Kemudian dari persamaan tersebut, saya membuat bentuk matriks dan menyelesaikannya dengan metode Gauss-Jordan.

Pembuktian Matematis 2


Dari Gauss-Jordan, saya mendapatkan hasil :

x1 = 10.3575 x2 = 14.55525 x3 = 20.00775


Kemudian setelah melakukan pembuktian secara matematis, saya menggunakan OpenModelica untuk menyelesaikan persoalan tersebut. Berikut adalah coding yang saya lakukan di OpenModelica.


Coding Spring Mass System


Lalu, saya melakukan verifikasi pengecekan coding.


Check Coding Spring OK
'


Setelah itu, saya simulate coding dan mendapatkan hasil sebagai berikut.


Result Spring OpenModelica


Dari simulasi pada OpenModelica, nilai x yang saya dapatkan adalah sebagai berikut :


x1 = 10.3575 x2 = 14.5552 x3 = 20.0077


Bisa dilihat bahwa dari pembuktian secara matematis maupun dari penyelesaian menggunakan simulasi OpenModelica, nilai dari X yang didapatkan bisa dikatakan sama (dengan perbedaan pada ketelitian dan angka dibelakang koma).


TUGAS PERTEMUAN III

Untuk tugas dalam pertemuan III (23/11/2020) lalu, pak Dai memberi tugas untuk menyelesaikan persamaan pada soal berikut dengan menggunakan OpenModelica.

Soal


Untuk menyelesaikan defleksi pada soal tersebut, saya meng-entry coding berikut pada software OpenModelica.

Coding

dan setelah melakukan simulasi, didapatkan hasil sebagai berikut, yang mana hasil ini sama dengan yang tertera pada buku.

Result


Namun, saat saya mengentry coding untuk external force, terdapat error pada saat saya akan melakukan simulasi. Saya masih berusaha melokasikan dimana titik errornya.

berikut adalah coding external force yang dinyatakan error oleh OpenModelica.


Error External Force Coding


QUIZ I METODE NUMERIK (30/11/2020)

Pada pertemuan IV ini, pak Dai memberikan dua soal kuis dan kami diminta untuk menuliskan flowchart pengerjaannya. Berikut adalah dua soal kuis tersebut.

Kuis


Kuis (2)


dan berikut adalah flowchart pengerjaan soal kuis tersebut yang telah saya tulis. Mengingat alur pengerjaan dua soal tersebut sama, maka flowchart berikut mencakup kedua soal tersebut.


Flowchart


PENYELESAIAN MANUAL KUIS

NO 4 Kuis no 4 ElitaK.jpg

NO 8

Kuis no 8 ElitaK.jpg

CODING

NO 4

MASTER CODE (CLASS)

class Quiz_no4 //establish Real parameter in order (element, i, j, theta, Area, E, length) 

 parameter Real[:,7]inisiasi = [1, 1, 0, 10e-4, 200e9;
                                2, 1, 0, 10e-4, 200e9;
                                3, 1, -1.25, 10e-4, 200e9;
                                4, 0, -1.25, 10e-4, 200e9;
                                5, -1, -1.25, 10e-4, 200e9];

//establish nodes 

 parameter Integer[size(inisiasi,1),2] node = [1, 2;
                                2, 3;
                                3, 4;
                                2, 4;
                                1, 4];

//number of nodes 

 parameter Integer n=4;

//xy boundary node points 

 parameter Integer [:] Boundary_xy = {1,3}; //boundary is fixed. no change to overall truss.

//x boundary nodes 

 parameter Integer [:] Boundary_x = {0};

//y boundary nodes 

 parameter Integer [:] Boundary_y = {0};

//load = [F1x, F1y,....,Fnx,Fny] 

 parameter Real [2*n] load = { 0, 0, 
                               -1035.28, -3863.70, 
                                0, 0,  
                                -1035.28, -3863.70};

 Real [size(inisiasi, 1)] L;
 Real [size(inisiasi, 1)] k;
 Real [size(inisiasi, 1), 4,4] Ke;
 Real [size(inisiasi, 1), 2*n, 2*n] Kg;
 Real [2*n, 2*n] KgTot;
 Real [2*n, 2*n] KgB;
 Real [2*n] U;
 Real [2*n] R;

//check force balance 

 Real [2] F;

equation 

 L = {(sqrt(inisiasi[i,2]^2 + inisiasi[i,3]^2)) for i in 1:size(inisiasi,1)};
 
 k = {(inisiasi[i,3] * inisiasi[i,4] / L[i]) for i in 1:size(inisiasi,1)};
 
Ke = StiffnessMatrixElement(inisiasi);
 
Kg = StiffnessMatrixGlobal(n, node, Ke);
 
KgTot = SumStiffnessMatrixGlobal(Kg);

KgB = BoundaryStiffnessMatrixGlobal(KgTot, Boundary_xy, Boundary_x, Boundary_y);

U = GaussJordan(KgB, load);

R = ReactionForce(KgTot, U, load);

F = CheckForce(load,R);

end Quiz_no4;


FUNCTIONS

Stiffness Matrix Element

function StiffnessMatrixElement 

 input Real [:,5] inisiasi;
 output Real [size(inisiasi,1),4,4] Ke;



 protected

   Real theta;

   Real [3] StiffTrig;

   Real [4,4] StiffTrans;

   Real [size(inisiasi,1)] k;

   Real [size(inisiasi,1)] L;

   Real float_error = 10e-10;


algorithm 

 L := {(sqrt(inisiasi[i,2]^2 + inisiasi[i,3]^2)) for i in 1:size(inisiasi,1)};



 k := {(inisiasi[i,3] * inisiasi[i,4] / inisiasi[i,5]) for i in 1:size(inisiasi,1)};



 // Finding stiffness matrix of each element member

 for i in 1:size(inisiasi,1) loop



 // Clearing the matrices

 StiffTrig := zeros(3);

 StiffTrans := zeros(4,4);



 // Converting degrees to radians

 theta := Modelica.SIunits.Conversions.from_deg(inisiasi[i,2]);



 // {cos^2, sin^2, sincos}

 StiffTrig := {(Modelica.Math.cos(theta))^2,

               (Modelica.Math.sin(theta))^2,

               (Modelica.Math.sin(theta)*Modelica.Math.cos(theta))};



 // Handle float error elements in StiffTrig

 for t in 1:size(StiffTrig,1) loop

   if abs(StiffTrig[t]) <= float_error then

     StiffTrig[t] := 0;

   end if;

 end for;



 // Construct stiffness transformation matrix

 StiffTrans := [  StiffTrig[1],    StiffTrig[3], -1*StiffTrig[1], -1*StiffTrig[3];

                  StiffTrig[3],    StiffTrig[2], -1*StiffTrig[3], -1*StiffTrig[2];

               -1*StiffTrig[1], -1*StiffTrig[3],    StiffTrig[1],    StiffTrig[3];

               -1*StiffTrig[3], -1*StiffTrig[2],    StiffTrig[3],    StiffTrig[2]];



 // Multiply in stiffness constant of element, add final stiffness matrix to Ke

 for m in 1:4 loop

   for n in 1:4 loop

     Ke[i,m,n] := k[i] * StiffTrans[m,n];

   end for;

 end for;



 end for;

end StiffnessMatrixElement;


Stiffness Matrix Global


function StiffnessMatrixGlobal 

 input Integer x;

 input Integer [:,2] n;

 input Real [:,4,4] Ke; 

 output Real [size(Ke,1),2*x,2*x] Kg;

algorithm 

 Kg := zeros(size(Ke,1),2*x,2*x);



 for i in 1:size(Ke,1) loop

   Kg[i,2*n[i,1],2*n[i,1]]:=Ke[i,2,2];

   Kg[i,2*n[i,1]-1,2*n[i,1]-1]:=Ke[i,1,1];

   Kg[i,2*n[i,1],2*n[i,1]-1]:=Ke[i,2,1];

   Kg[i,2*n[i,1]-1,2*n[i,1]]:=Ke[i,1,2];



   Kg[i,2*n[i,2],2*n[i,2]]:=Ke[i,4,4];

   Kg[i,2*n[i,2]-1,2*n[i,2]-1]:=Ke[i,3,3];

   Kg[i,2*n[i,2],2*n[i,2]-1]:=Ke[i,4,3];
  
   Kg[i,2*n[i,2]-1,2*n[i,2]]:=Ke[i,3,4];



   Kg[i,2*n[i,2],2*n[i,1]]:=Ke[i,4,2];

   Kg[i,2*n[i,2]-1,2*n[i,1]-1]:=Ke[i,3,1];

   Kg[i,2*n[i,2],2*n[i,1]-1]:=Ke[i,4,1];

   Kg[i,2*n[i,2]-1,2*n[i,1]]:=Ke[i,3,2];



   Kg[i,2*n[i,1],2*n[i,2]]:=Ke[i,2,4];

   Kg[i,2*n[i,1]-1,2*n[i,2]-1]:=Ke[i,1,3];

   Kg[i,2*n[i,1],2*n[i,2]-1]:=Ke[i,2,3];

   Kg[i,2*n[i,1]-1,2*n[i,2]]:=Ke[i,1,4];

 end for;

end StiffnessMatrixGlobal;

Sum Stiffness Matrix Global

function SumStiffnessMatrixGlobal 

 input Real [:,:,:] Kg;
 output Real [size(Kg,2),size(Kg,2)] KgTot;

algorithm 

  for a in 1:size(Kg,2) loop
   for b in 1:size(Kg,2) loop
     KgTot[a,b] := sum(Kg [:,a,b]);
    end for;
   end for;

end SumStiffnessMatrixGlobal;


NO 8

MASTER CODE (CLASS)

class Quiz8 

 //inisiasi = [ elemen#, dX, dY, dZ, A, E]
 parameter Real [:,6] inisiasi = [1,  6,  0, -3, 1.56, 10.6e6; //isi sesuai data
                                  2,  0,  0, -6, 1.56, 10.6e6;
                                  3,  0,  6, -3, 1.56, 10.6e6;
                                  4, -6,  0, -3, 1.56, 10.6e6;
                                  5, -6,  6,  0, 1.56, 10.6e6;
                                  6,  0,  6,  3, 1.56, 10.6e6];

 //node = [ i, j]                                 
 parameter Integer [size(inisiasi,1),2] node = [1, 2; //isi sesuai data
                                                1, 3;
                                                1, 4;
                                                2, 3;
                                                2, 4;
                                                3, 4];

 //jumlah node
 parameter Integer n = 4; //isi sesuai data

 //titik node boundary xyz
 parameter Integer [:] Boundary_xyz = {1}; //isi sesuai data

 //titik node boundary xy
 parameter Integer [:] Boundary_xy = {4}; //isi sesuai data

 //titik node boundary xz
 parameter Integer [:] Boundary_xz = {0}; //isi sesuai data

 //titik node boundary yz
 parameter Integer [:] Boundary_yz = {0}; //isi sesuai data

 //titik node boundary x
 parameter Integer [:] Boundary_x = {3}; //isi sesuai data

 //titik node boundary y
 parameter Integer [:] Boundary_y = {0}; //isi sesuai data

 //titik node boundary z
 parameter Integer [:] Boundary_z = {0}; //isi sesuai data
                            
 //load = [ F1x, F1y, F1z,..., Fnx, Fny, Fnz]
 parameter Real [3*n] load = {0,    0, 0,  //isi sesuai data
                              0, -200, 0, 
                              0,    0, 0, 
                              0,    0, 0}; 

 Real [size(inisiasi,1)] L;
 Real [size(inisiasi,1)] k;
 Real [size(inisiasi,1),6,6] Ke;
 Real [size(inisiasi,1),3*n,3*n] Kg;
 Real [3*n,3*n] KgTot;
 Real [3*n,3*n] KgB;
 Real [3*n] U;
 Real [3*n] R;

 //check force
 Real [3] F;

equation 

L = {(sqrt(inisiasi[i,2]^2 + inisiasi[i,3]^2 + inisiasi[i,4]^2)) for i in 1:size(inisiasi,1)}; 

k = {(inisiasi[i,5] * inisiasi[i,6] / L[i]) for i in 1:size(inisiasi,1)};

Ke = StiffnessMatrixElement(inisiasi);

Kg = StiffnessMatrixGlobal(n, node, Ke);

KgTot = SumStiffnessMatrixGlobal(Kg);

KgB = BoundaryStiffnessMatrixGlobal(KgTot, Boundary_xyz, Boundary_xy, Boundary_xz, Boundary_yz, Boundary_x, Boundary_y, Boundary_z);

U = GaussJordan(KgB, load);

R = ReactionForce(KgTot, U, load);

F = CheckForce(load,R);

end Quiz8;


TUGAS PERTEMUAN V

Untuk pertemuan ini kami ditugaskan mengaplikasikan coding 3D Truss milik Ahmad Mohammad Fahmi ke permasalahan 3D Truss yang lainnya.

Tugas 3D ElitaK.jpeg

Calculate Length and Cosine class Length_Cos

parameter Integer Points=4; parameter Integer Trusses=6; parameter Integer C[Trusses,2]=[1,2; 

                               1,3;
                               1,4;
                               2,3;
                               3,4;
                               2,4];

parameter Integer P[:,3]=[0,0,3; 

                         6,0,0;
                         0,0,-3;
                         0,6,0];

//Element 1,4,5 Real L1 [size(P,1)-1]; 

 Real cosx1 [size(P,1)-1];
  Real cosy1 [size(P,1)-1];
   Real cosz1 [size(P,1)-1];

//Element 2 6 Real L2 [size(P,1)-2]; 

 Real cosx2 [size(P,1)-2];
  Real cosy2 [size(P,1)-2];
   Real cosz2 [size(P,1)-2];

//Element 3 Real L3 [size(P,1)-3]; 

 Real cosx3 [size(P,1)-3];
  Real cosy3 [size(P,1)-3];
   Real cosz3 [size(P,1)-3];

equation 

 //Element 1 4 5
 for i in 1:1:(size(P,1)-1) loop
   L1[i] = sqrt((P[i+1,1]-P[i,1]^2+(P[i+1,2]-P[i,2])^2+(P[i+1,3]-P[i,3])^2));
 end for;
  for i in 1:(size(P,1)-1) loop
    cosx1 [i] = (P[i+1,1]-P[i,1])/L1[i];
    cosy1 [i] = (P[i+1,2]-P[i,2])/L1[i];
    cosz1 [i] = (P[i+1,3]-P[i,3])/L1[i];
  end for;
  
 //Element 2 6
 for i in 1:1:(size(P,1)-2) loop
   L2 [i] = sqrt((P[i+2,1]-P[i,1]^2+(P[i+2,2]-P[i,2])^2+(P[i+2,3]-P[i,3])^2));
 end for;
  for i in i:(size(P,1)-2) loop
   cosx2 [i] = (P[i+2,1]-P[i,1])/L2[i];
   cosy2 [i] = (P[i+2,2]-P[i,2])/L2[i];
   cosz2 [i] = (P[i+2,3]-P[i,3])/L2[i];
  end for;
  
 //Element 3
 for i in 1:1:(size(P,1)-3) loop
   L3 [i] = sqrt((P[i+3,1]-P[i,1]^2+(P[i+3,2]-P[i,2])^2+(P[i+3,3]-P[i,3])^2));
 end for;
  for i in i:(size(P,1)-3) loop
   cosx3 [i] = (P[i+3,1]-P[i,1])/L3[i];
   cosy3 [i] = (P[i+3,2]-P[i,2])/L3[i];
   cosz3 [i] = (P[i+3,3]-P[i,3])/L3[i];
  end for;

end Length_Cos;


MASTER CODE (CLASS)

style='border-style: none none solid solid;'

class Soal5 

 //inisiasi = [ elemen#, dX, dY, dZ, A, E]
 parameter Real [:,6] inisiasi = [1,  0.8944,  0.0000, -0.4472, 1.56, 10.6e6; //isi sesuai data
                                  2,  0.0000,  0.0000, -1.0000, 1.56, 10.6e6;
                                  3,  0.0000,  0.8944, -0.4472, 1.56, 10.6e6;
                                  4, -0.8944,  0.0000, -0.4472, 1.56, 10.6e6;
                                  5, -0.7071,  0.7071,  0.0000, 1.56, 10.6e6;
                                  6,  0.0000,  0.8944,  0.4472, 1.56, 10.6e6];

 //node = [ i, j]                                 
 parameter Integer [size(inisiasi,1),2] node = [1, 2; //isi sesuai data
                                                1, 3;
                                                1, 4;
                                                2, 3;
                                                2, 4;
                                                3, 4];

 //jumlah node
 parameter Integer n = 4; //isi sesuai data

 //titik node boundary xyz
 parameter Integer [:] Boundary_xyz = {1,3,4}; //isi sesuai data

 //titik node boundary xy
 parameter Integer [:] Boundary_xy = {3}; //isi sesuai data

 //titik node boundary xz
 parameter Integer [:] Boundary_xz = {0}; //isi sesuai data

 //titik node boundary yz
 parameter Integer [:] Boundary_yz = {0}; //isi sesuai data

 //titik node boundary x
 parameter Integer [:] Boundary_x = {2}; //isi sesuai data

 //titik node boundary y
 parameter Integer [:] Boundary_y = {4}; //isi sesuai data

 //titik node boundary z
 parameter Integer [:] Boundary_z = {1}; //isi sesuai data
                            
 //load = [ F1x, F1y, F1z,..., Fnx, Fny, Fnz]
 parameter Real [3*n] load = {0,    0, 0,  //isi sesuai data
                              0, -200, 0, 
                              0,    0, 0, 
                              0,    0, 0}; 

 Real [size(inisiasi,1)] L;
 Real [size(inisiasi,1)] k;
 Real [size(inisiasi,1),6,6] Ke;
 Real [size(inisiasi,1),3*n,3*n] Kg;
 Real [3*n,3*n] KgTot;
 Real [3*n,3*n] KgB;
 Real [3*n] U;
 Real [3*n] R;

 //check force
 Real [3] F;

equation 

L = {(sqrt(inisiasi[i,2]^2 + inisiasi[i,3]^2 + inisiasi[i,4]^2)) for i in 1:size(inisiasi,1)}; 

k = {(inisiasi[i,5] * inisiasi[i,6] / L[i]) for i in 1:size(inisiasi,1)};

Ke = StiffnessMatrixElement(inisiasi);

Kg = StiffnessMatrixGlobal(n, node, Ke);

KgTot = SumStiffnessMatrixGlobal(Kg);

KgB = BoundaryStiffnessMatrixGlobal(KgTot, Boundary_xyz, Boundary_xy, Boundary_xz, Boundary_yz, Boundary_x, Boundary_y, Boundary_z);

U = GaussJordan(KgB, load);

R = ReactionForce(KgTot, U, load);

F = CheckForce(load,R);


end Soal5;

Saya mengikuti langkah-langkah coding yang diajarkan oleh Ahmad Mohammad Fahmi dan saat telah melakukan verifikasi code, seluruh code dinyatakan OK pada OpenModelica. Namun pada saat simulate, dalam proses runningnya, coding mengalami crash sehingga tidak dapat di plotting. Saya masih meneliti coding untuk melakukan revisi agar tidak terjadi crash.


OPTIMASI

Pada kesempatan ini kami mempelajari mengenai optimasi dan aplikasi metoda numeriknya. Berikut adalah coding pada OpenModelica yang kami pelajari;


FUNCTION

function Func_Optimization input Real x; import Modelica.Math; output Real y;

algorithm y:=2*sin(x)-x^2/10;

end Func_Optimization;


OPTIMIZATION GOLDEN SECTION

model Opt_Gold

parameter Real xlo=0; parameter Real xhi=4; parameter Integer N=8; // maximum iteration parameter Real es=0.0001; // maximum error

Real f1[N], f2[N], x1[N], x2[N], ea[N]; Real xopt, fx; protected Real d, xl, xu, xint, R=(5^(1/2)-1)/2;

algorithm xl := xlo; xu := xhi;

for i in 1:N loop 

d:= R*(xu-xl);
x1[i]:=xl+d;
x2[i]:=xu-d;
f1[i]:=Func_Optimization(x1[i]);
f2[i]:=Func_Optimization(x2[i]);
xint:=xu-xl;

if f1[i]>f2[i] then
  xl:=x2[i];
  xopt:=x1[i];
  fx:=f1[i];
  else
    xu:=x1[i];
    xopt:=x2[i];
    fx:=f2[i];
end if;

ea[i]:=(1-R)*abs((xint)/xopt);
if ea[i]<es then
  break;
end if;

end for;

end Opt_Gold;


NEWTON OPTIMIZATION

model Opt_Newton

parameter Real g=2.5; //initial guess parameter Integer N=8; // maximum iteration parameter Real es=0.0001; // maximum error

Real X[N]; Real xopt, fx, ea[N];

algorithm X[1]:=g; ea[1]:=1;

for i in 2:N loop 

X[i]:=X[i-1]-Func_Optimization_Der(X[i-1])/Func_Optimization_Der_Der(X[i-1]);
xopt:=X[i];

ea[i]:=abs(1-X[i-1]/X[i]);
if ea[i]<es then
  break;
end if;

end for; fx:=Func_Optimization(xopt);

end Opt_Newton;


PARABOLIC OPTIMIZATION

model Opt_Parabolic

parameter Real g1=0; //initial guess parameter Real g2=1; //initial guess parameter Real g3=4; //initial guess parameter Integer N=5; // maximum iteration parameter Real es=0.0001; // maximum error

Real x1, x2, x3, xopt, xp[N], ea[N]; //Real xl, xm, xu; Real fx1, fx2, fx3, fx, A[4], A_star[4]; algorithm x1:=g1; x2:=g2; x3:=g3;

for i in 1:N loop 

fx1:=Func_Optimization(x1);
fx2:=Func_Optimization(x2);
fx3:=Func_Optimization(x3);

xp[i]:=(fx1*(x2^2-x3^2)+fx2*(x3^2-x1^2)+fx3*(x1^2-x2^2))/(2*fx1*(x2-x3)+2*fx2*(x3-x1)+2*fx3*(x1-x2));
xopt:=xp[i];
fx:=Func_Optimization(xp[i]);
A:={x1,x2,x3,xp[i]};
A_star:=Modelica.Math.Vectors.sort(A);
 
if xp[i]>x2 then
  x1:=A_star[2];
  x2:=A_star[3];
  x3:=A_star[4];
  else
    x1:=A_star[1];
    x2:=A_star[2];
    x3:=A_star[3];
end if;

end for;

ea[1]:=1; for i in 2:N loop 

ea[i]:=abs(1-xp[i-1]/xp[i]);
if ea[i]<es then
  break;
end if;

end for;

end Opt_Parabolic;


TUGAS BESAR METODE NUMERIK

OBJEKTIF

Tujuan dari tugas besar ini adalah>> optimasi harga pembuatan rangka truss sederhana dengan memvariasi dimensi dan elastisitas material.


GEOMETRI

Geometri dan Load Truss


ASUMSI


Dalam tugas besar ini digunakan asumsi-asumsi sebagai berikut:

1. Beban hanya terdistribusi pada node.

2. Safety Factor minimal 2.

3. Batas displacement pada truss paling atas sebelum buckling bernilai 0,001 m.

4. Ketinggian trusses pada setiap lantai konstan, sebesar 0,6 m.


METODOLOGI


I. ELASTICITY LOCKED

1. Mencari data yield strength dan elasticity dari material yang dipilih (disini digunakan SS400/ASTM A36)

2. Mencari harga truss sesuai dengan material.

3. Menghitung nilai safety factor.

4. Membuat rasio (efficiency) dari safety factor dengan total cost.

5. Membuat persamaan (mencari Coe) antara rasio dengan area menggunakan curve-fitting.

6. Melakukan optimasi menggunakan metode golden section.


II. AREA LOCKED

1. Mendefinisikan luas penampang yang sama untuk seluruh variasi material >> 3,75 x 10^-4 m^2.

2. Mencari harga untuk 3 jenis variasi material >> SS400, SS304, SS316.

3. Menghitung nilai safety factor pada 3 variasi truss.

4. Membuat rasio antara safety factor dengan total cost.

5. Membuat persamaan (mencari Coe) antara rasio dengan area menggunakan curve-fitting.

6. Melakukan optimasi menggunakan metode golden section.


KALKULASI DISPLACEMENT, REACTION FORCE, & STRESS VIA OPENMODELICA

model Trusses_3D_Tugas_Besar_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.000224; //Area L Profile (Dimension=0.03, Thickness=0,004) (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 Trusses_3D_Tugas_Besar_Safety;


KALKULASI COE UNTUK CURVE-FITTING VIA OPENMODELICA

function Curve_Fitting


input Real X[:];

input Real Y[size(X,1)];

input Integer order=2;

output Real Coe[order+1];


protected

Real Z[size(X,1),order+1];

Real ZTr[order+1,size(X,1)];

Real A[order+1,order+1];

Real B[order+1];


algorithm

for i in 1:size(X,1) loop 

 for j in 1:(order+1) loop
 Z[i,j]:=X[i]^(order+1-j);
 end for;

end for; ZTr:=transpose(Z);

A:=ZTr*Z;

B:=ZTr*Y;

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

//Coe:=fill(2,size(Coe,1));

end Curve_Fitting; /* for i in 1:3 loop 

for j in 1:Points loop
 R[j]:=reaction[3*(j-1)+i];
end for;
Sur[i]:=sum(R);

end for;

  • /


model CurveFitting

parameter Real X[6]={141e-6,375e-6,384e-6,575e-6,931e-6,864e-6};//area

parameter Real Y[6]={7702,7833,7664,7708,7927,7805};//Cost/Kg

Real Coe[3];

algorithm

Coe:=Curve_Fitting(X,Y,2);

end CurveFitting;


TABULASI DATA DI EXCEL & SNAPSHOT Coe


I. ELASTICITY LOCKED EXCEL TABULATION


Elasticity locked ElitaK.png


Coe for Curve-Fitting

Coef elasticity ElitaK.png


II. AREA LOCKED EXCEL TABULATION


Area locked ElitaK.png


Coe for Curve-Fitting


Curve-Fitting for Cost/Kg


Curve-Fitting for Density


Curve-Fitting for Yield Strength



UAS METODE NUMERIK (13/01/2021)

Terlampir di bawah ini adalah jawaban untuk Ujian Akhir Semester Metode Numerik saya pada hari ini Rabu, 13 Januari 2021.



UAS 1 ElitaK.jpeg


UAS 2 ElitaK.jpeg


UAS 3 ElitaK.jpeg


UAS 45 ElitaK.jpg


UAS 67 ElitaK.jpg


UAS ss1 ElitaK.png


UAS ss2 ElitaK.png
Retrieved from "http://air.eng.ui.ac.id/index.php?title=Metnum03_-_Elita_Kabayeva&oldid=56345"