Difference between revisions of "Christopher Salendu Erwin"

From ccitonlinewiki
Jump to: navigation, search
(→‎Numerical Methods)
(→‎Complete Code)
 
(130 intermediate revisions by the same user not shown)
Line 3: Line 3:
 
== About Me ==
 
== About Me ==
  
−
[[File:christopherseprofile.jpg|200px|thumb|right]]
+
[[File:christopherseprofile.jpg|200px|thumb|right|Freut mich!]]
  
−
My name is Christo, and I am a student at the University of Indonesia pursuing a Bachelor's in Mechanical Engineering. I am fascinated with the fields of aerospace, robotics, and systems design.
+
I am a student at the University of Indonesia pursuing a Bachelor's in Mechanical Engineering. I am fascinated with the fields of aerospace, robotics, and systems design.
  
−
On-campus I participate in two extracurricular organizations: English debating and robotics. As a debater, I enjoy exploring multiple perspectives on issues and deconstructing policies with critical argumentation. Although I consider myself still relatively new to debating, I am thrilled to say that I have earned a Novice Champion title in the JOVED 2019 debating tournament. Besides debating, I recently joined the robotics club's rocketry division. Since the pandemic happened right after, I haven't done much with robotics yet, but I am looking forward to learning more about it!
+
On-campus I participate in two extracurricular organizations: English debating and robotics. As a debater, I enjoy exploring multiple perspectives on relevant issues and deconstructing policies through critical argumentation. Although I consider myself still far from an expert in debating, I am thrilled to say that I have achieved the Novice Champion title in the JOVED 2019 debating tournament. Besides debating, I recently joined the robotics club's rocketry division. Unfortunately, since the pandemic happened right after, I haven't done much with robotics yet, but I am looking forward to learning more about it!
−
 
 
−
Some my interests include programming and foreign languages. For programming, I'd say I'm currently grappling with Python, and I know a little bit of C++. Another hobby of mine is foreign languages, and I'm able to speak some German, Japanese, Spanish, and French while studying several other languages when I'm feeling up to it.
 
  
 +
Some of my interests include programming and foreign languages. For programming, I have some experience with microcontrollers, Arduino, and Raspberry Pi. The programming language that I'm currently grappling with is Python, and I know a little bit of C++, Java, and HTML. Another hobby of mine is learning foreign languages, and I'm able to speak some German, Japanese, Spanish, and French. Sometimes I study several other languages when I'm feeling up to it.
  
 
== Numerical Methods ==
 
== Numerical Methods ==
Line 25: Line 24:
 
* Interpolation
 
* Interpolation
 
* Numerical Differentiation
 
* Numerical Differentiation
 +
 +
=== Week 01 [11/11/20] ===
 +
 +
==== Tujuan Belajar Metode Numerik ====
 +
 +
# Memahami konsep dan prinsip dasar metode numerik.
 +
# Mengerti aplikasi dan penerapan metode numerik.
 +
# Mampu menerapkan pengetahuan metode numerik di persoalan teknik.
 +
# Mendapat nilai tambah sendiri agar menjadi orang yang lebih beradab.
 +
 +
==== Belajar OpenModelica ====
 +
 +
[[File:LearnOpenModelicaChristopherSE.mkv|720px|center]]
 +
 +
Di video ini, saya berusaha untuk menjelaskan cara basic menggunakan OpenModelica. Saya mencoba beberapa fungsi dasar seperti modelling di text view, modelling di diagram view, dan melakukan simulasi dan plotting. Contoh simulasi sistem yang digunakan adalah:
 +
* Pemodelan sistem fluida dalam tangki berbentuk bola dengan menulis code
 +
* Pemodelan sistem osilasi massa spring-damper dengan merancang komponen visual
 +
 +
=== Week 02 [18/11/20] ===
 +
 +
==== OpenModelica Intro ====
 +
 +
OpenModelica adalah aplikasi yang digunakan untuk melakukan pemodelan. OpenModelica memudahkan penyelesaian permasalahan permodelan sistem yang cukup rumit. Coding yang dilakukan dalam OMEdit ditulis berdasarkan bahasa C++.
 +
 +
OpenModelica digunakan untuk kelas Metode Numerik untuk beberapa alasan:
 +
* OpenModelica memampukan kita untuk melakukan pemodelan sitem teknik dengan mudah tanpa terlalu mendalami pengetahuan pemograman dan ilmu komputer.
 +
* OpenModelica adalah software yang Open Source, dengan akses bebas dan terbuka. Karena itu, terdapat banyak pengguna dan software dapat digunakan gratis dengan bebas secara legal.
 +
 +
===== Class & Function =====
 +
 +
[[File:ModelicaClassSpecial.png|Pilihan ''Specialization'' ketika membuat ''Modelica Class'' yang baru.|thumb|upright]]
 +
 +
Ketika membuat file ''Class'' baru di OpenModelica, ada pilihan ''Specialization'' dimana jenis ''Class'' dapat dipilih. Setiap ''Class'' mempunyai karakteristik dan tujuan masing-masing. Informasi lebih lanjut tentang setiap jenis ''Class'' bisa dibaca dalam dokumentasi [//build.openmodelica.org/Documentation/ModelicaReference.Classes.html disini].
 +
 +
Salah satu jenis ''Class'' yang sudah digunakan adalah ''Model'' yang dapat melakukan pemodelan dan simulasi. Ada jenis ''Class'' lain yaitu ''Function'' dimana code yang diketik di dalamnya akan menjadi suatu fungsi yang dapat dipanggil dalam model. Ketika dibuat, ''Function'' mempunyai bagian ''algorithm'' dimana perhitungan fungsi dimasukkan. Sebelum bagian ''algorithm'', dapat dimasukkan ''input'' dan ''output'' yang akan ditangani oleh fungsinya.
 +
 +
Sebagai contoh, dibikin suatu fungsi ''PlusTen'' yang menambahkan angka 10 ke nilai input x. Fungsi ini terus akan dipanggil dalam model ''TestPlusTen'' yang menerima parameter ''num'' dan menghasilkan variabel ''answer''. Ketika dicompile, hasil penggunaan fungsi di model ini dapat dilihat.
 +
 +
<gallery heights=175px widths=225px>
 +
File:PlusTenFunction.png|Contoh fungsi yang menambahkan 10 ke nilai input.
 +
File:PlusTenModel.png|Fungsi dipanggil dalam model.
 +
PlusTenResult.png|Hasil model yang menggunakan fungsi.
 +
</gallery>
 +
 +
==== Persamaan Aljabar Simultan ====
 +
 +
OpenModelica dapat melakukan perhitungan aljabar linier. Matriks ditulis sebagai ''Array'' yang dapat dilakukan kalkulasi seperti perkalian, transpose, inverse, dll. Informasi lebih detail tentang array dapat dibaca [//www.openmodelica.org/images/docs/Courses/lecture06%20-%20arrays%20algorithms%20and%20functions.pdf di sini] dan [//mbe.modelica.university/behavior/arrays/declarations/ di sini]. Dokumentasi operator dan fungsi matriks dalam OpenModelica dabat dibaca [//build.openmodelica.org/Documentation/ModelicaReference.Operators.html di sini] dan [//www.maplesoft.com/documentation_center/online_manuals/modelica/Modelica_Math_Matrices.html di sini].
 +
 +
Salah satu fungsi di OpenModelica adalah solve(). Fungsi ini dapat menyelesaikan suatu sistem persamaan linier A*x=B dan mencari solusi x dengan metode Gaussian elimination dengan partial pivoting. Dokumentasi tentang fungsi ini dapat dibaca [//build.openmodelica.org/Documentation/Modelica.Math.Matrices.solve.html di sini]. Berikut adalah video yang menjelaskan tentang fungsi ini dan contoh pemakaiannya.
 +
 +
[[File:LinearEquationSystemChristopherSE.mkv|720px|LinearEquationSystem.mkv|center]]
 +
 +
===== Gauss-Jordan =====
 +
 +
Berikut adalah contoh kode fungsi Modelica yang melakukan operasi Gauss-Jordan matriks menjadi bentuk reduced row echelon.
 +
 +
<syntaxhighlight lang="modelica">
 +
function GaussJordan
 +
 +
// 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
 +
    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;
 +
 +
</syntaxhighlight>
 +
 +
=== Week 03 [25/11/20] ===
 +
 +
==== Metode Numerik dan Masalah Teknik ====
 +
 +
[[File:AplikasiMetNumMasalahTeknik.jpg|Aplikasi metode numerik dalam permasalahan teknik.|thumb|upright]]
 +
 +
Masalah teknik dapat diselesaikan dengan menggunakan metode numerik. Beberapa contoh metode numerik yang sering digunakan untuk melakukan analisa masalah teknik adalah [//www.simscale.com/docs/simwiki/cfd-computational-fluid-dynamics/what-is-cfd-computational-fluid-dynamics/ Computation Fluid Dynamics (CFD)], [//www.simscale.com/docs/simwiki/fea-finite-element-analysis/what-is-fea-finite-element-analysis/ Finite Element Analysis (FEA)], dan Metode Stokastik. CFD dan FEA berdasarkan hukum-hukum fisika, sementara metode stokastik berbasis data dan statistik. Salah satu contoh metode stokastik yang sering digunakan adalah metode [//towardsdatascience.com/an-overview-of-monte-carlo-methods-675384eb1694 Monte Carlo].
 +
 +
Untuk menyelesaikan suatu masalah teknik menggunakan metode numerik, perlu dilakukan beberapa tahap. Pertama, masalah teknik perlu dilakukan analisis. Setelah itu, hasil analisis masalah perlu dilakukan pemodelan matematis, yang kemudian diterjemahkan menjadi model numerik yang dapat dihitung oleh sebuah komputer. Dengan melakukan perhitungan di software komputer, dapat diperolehkan solusi dari masalah teknik.
 +
 +
==== Truss: Finite Element Formulation ====
 +
 +
[[File:StaticalDeterminacy.jpg|Contoh masalah ''statically determinate'' dan ''statically indeterminate''.|thumb|upright]]
 +
 +
Truss adalah struktur yang terdiri dari beam lurus yang saling dihubungkan di ujungnya. Masalah teknik yang terkait dengan truss dapat dibagikan menjadi ''statically determinate'' atau [//en.wikipedia.org/wiki/Statically_indeterminate ''statically indeterminate'']. Pada masalah jenis ''statically determinate'' dilakukan asumsi bahwa truss yang dianalisis merupakan ''rigid body'' dan defleksi tidak diperhitungkan. Masalah jenis ini dapat diselesaikan dengan ''methods of joints or sections'' yang pernah dibahaskan di materi teknik yang lebih dasar. Masalah jenis ''statically indeterminate'' tidak dapat dilakukan analis dengan asumsi ''rigid body'', sehingga diperlukan ''finite element method'' untuk menyelesaikannya.
 +
 +
'''Equivalent Stiffness: Average Stress, Average Strain, & Hooke's Law'''
 +
 +
''(under construction)''
 +
 +
'''Displacements & Forces: Global & Local Coordinates'''
 +
 +
''(under construction)''
 +
 +
'''Internal Forces & Stiffness Matrix'''
 +
 +
''(under construction)''
 +
 +
==== Naive Gauss Algorithm ====
 +
Berikut adalah contoh ''pseudocode'' fungsi yang menggunakan algoritma ''Naive Gauss Elimination''. Fungsi ini dapat digunakan untuk menyelesaikan sistem persamaan linear di mana rank matriks sama dengan jumlah baris.
 +
 +
<syntaxhighlight lang="modelica">
 +
function NaiveGauss
 +
 +
/*
 +
// Pseudocode Figure 9.4
 +
 +
//Forward Elimination
 +
DOFOR k = 1, n - 1
 +
  DOFOR i = k + 1, n
 +
    factor = a[i,k] / a[k,k]
 +
    DOFOR j = k + 1 to n
 +
      a[i,j] = a[i,j] - factor * a[k,j]
 +
    END DO
 +
    b[i] = b[i] - factor * b[k]
 +
  END DO
 +
END DO
 +
 +
// Back Substitution
 +
x[n] = b[n] / a[n,n]
 +
DOFOR i = n - 1, 1, -1
 +
  sum = b[i]
 +
  DOFOR j = i + 1, n
 +
    sum = sum = a[i,j] * x[j]
 +
  END DO
 +
  x[i] = sum / a[i,i]
 +
END DO
 +
*/
 +
 +
input Real [:,:] A; // Left-hand Coefficients of the Linear Equation System in array form
 +
input Real [:] B;  // Right-hand Constants of the Linear Equation System in array form
 +
output Real [:,:] y;// Array containing coefficient matrix after NGE operation
 +
output Real [:] x;  // Array containing solved values of x
 +
 +
protected
 +
Real [:,:] a;
 +
Real [:] b;
 +
Integer m = size(A,1);  // Number of rows in matrix
 +
Integer n = size(A,2);  // Number of columns in matrix
 +
Real k = 1;      // Pivot column pointer
 +
Real i = 1;      // Row counter
 +
Real j = 1;      // Row element counter
 +
Real factor = 1;  // Factor value used for forward elimination
 +
Real sum = 1;    // Sum value used for back substitution
 +
 +
algorithm
 +
 +
// Transfer input matrix (A,B) into variables (a,b)
 +
a := A;
 +
b := B;
 +
 +
// Forward Elimination
 +
for k in 1:(n-1) loop
 +
  for i in (k+1):n loop
 +
    factor := a[i,k] / a[k,k];
 +
    for j in (k+1):n loop
 +
      a[i,j] := a[i,j] - (factor * a[k,j]);
 +
    end for;
 +
    b[i] := b[i] - (factor * b[k]);
 +
  end for;
 +
end for;
 +
 +
// Back Substitution
 +
x[n] := b[n] / a[n,n];
 +
for i in (n-1):(-1) loop
 +
  sum := b[i];
 +
  for j in (i+1):n loop
 +
    sum := sum - (a[i,j] * x[j]);
 +
  end for;
 +
  x[i] := sum / a[i,i];
 +
end for;
 +
 +
end NaiveGauss;
 +
 +
</syntaxhighlight>
 +
 +
==== 2D Truss Analysis using OpenModelica ====
 +
 +
[[File:Tugas minggu 3.jpeg|600px|thumb|center|Chapter 3 Problem 4 (''Finite Element Analysis, 4th Edition, Saeed Moaveni'')]]
 +
 +
===== Class: Solving 2D Truss Problem =====
 +
{| class=wikitable
 +
| style="vertical-align:top"|'''Declarations: Parameters and Variables'''
 +
<syntaxhighlight lang="modelica">
 +
// DECLARATIONS
 +
// Data for each member: [element #, node i, node j, theta, area, modulus, length]
 +
// SI Units: [integer, integer, integer, degrees, meters^2, pascals, meters]
 +
// Imperial Units: [integer, integer, integer, degrees, inches^2, lb/inches, inches]
 +
parameter Real [:,7] member = [1, 1, 2,      0, 10e-4, 200e9, 1;
 +
                              2, 2, 3,      0, 10e-4, 200e9, 1;
 +
                              3, 1, 4, 308.66, 10e-4, 200e9, 1.60;
 +
                              4, 2, 4, 270.00, 10e-4, 200e9, 1.25;
 +
                              5, 3, 4, 231.34, 10e-4, 200e9, 1.60];
 +
 +
// External loads and boundary conditions for each node: [node #, FX, FY, UX(1:Fixed), UY(1:Fixed)]
 +
parameter Real [:,5] node_load = [1,        0,        0, 1, 1;
 +
                                  2, -1035.28, -3863.70, 0, 0;
 +
                                  3,        0,        0, 1, 1;
 +
                                  4, -1035.28, -3863.70, 0, 0];
 +
 +
// Vector for equivalent stiffness constant of each member
 +
Real k [size(member,1)];
 +
 +
// Arrays for Stiffness Matrices
 +
Real Ke [size(member,1), 4, 4]; // Element
 +
Real KGe [size(member,1), 2*size(node_load,1), 2*size(node_load,1)]; // Element in Global Form
 +
Real KGt [2*size(node_load,1), 2*size(node_load,1)]; // Total Global
 +
 +
// Displacement Matrix
 +
Real U [2*size(node_load,1)];
 +
 +
// Postprocessing Phase
 +
Real R [2*size(node_load,1)]; // Reaction forces for each node
 +
Real InForce_NormStress [size(member,1),7]; // Internal forces and normal stress for each member
 +
// Format: [member #, uix, uiy, ujx, ujy, internal force, normal stress]
 +
 +
</syntaxhighlight>
 +
 +
| style="vertical-align:top"|'''Equations'''
 +
<syntaxhighlight lang="modelica">
 +
// EQUATIONS
 +
equation
 +
 +
// Calculate stiffness constants for each member element
 +
k = {(member[i,5] * member[i,6] / member[i,7]) for i in 1:size(member,1)};
 +
 +
// Find stiffness matrix for each element
 +
Ke = StiffnessMatrixElement(member);
 +
 +
// Assemble the global stiffness matrix
 +
KGe = StiffnessMatrixGlobal(member, node_load);
 +
KGt = StiffnessMatrixFinal(KGe);
 +
 +
// Solve for the displacement matrix
 +
U = DisplacementMatrix(node_load, KGt);
 +
 +
// Postprocessing Phase
 +
(R, InForce_NormStress) = TrussPostProcess(member, node_load, KGt, U);
 +
</syntaxhighlight>
 +
 +
|}
 +
 +
===== Functions: Solving 2D Truss Problem =====
 +
{| class=wikitable
 +
| style="vertical-align:top"|'''Function: Local Element Stiffness Matrix'''
 +
<syntaxhighlight lang="modelica">
 +
function StiffnessMatrixElement
 +
 +
input Real [:,7] member_mat;
 +
output Real [size(member_mat,1),4,4] Ke_mat;
 +
 +
protected
 +
Real theta;
 +
Real [3] StiffTrig;
 +
Real [4,4] StiffTrans;
 +
Real [size(member_mat,1)] k_vec;
 +
Real float_error = 10e-10;
 +
 +
algorithm
 +
 +
k_vec := {(member_mat[i,5] * member_mat[i,6] / member_mat[i,7]) for i in 1:size(member_mat,1)};
 +
 +
// Finding stiffness matrix of each element member
 +
for i in 1:size(member_mat,1) loop
 +
 +
  // Clearing the matrices
 +
  StiffTrig := zeros(3);
 +
  StiffTrans := zeros(4,4);
 +
 
 +
  // Converting degrees to radians
 +
  theta := Modelica.SIunits.Conversions.from_deg(member_mat[i,4]);
 +
 +
  // {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_mat
 +
  for m in 1:4 loop
 +
    for n in 1:4 loop
 +
      Ke_mat[i,m,n] := k_vec[i] * StiffTrans[m,n];
 +
    end for;
 +
  end for;
 +
 +
end for;
 +
 +
end StiffnessMatrixElement;
 +
 +
</syntaxhighlight>
 +
 +
| style="vertical-align:top"|'''Function: Global Element Stiffness Matrix'''
 +
<syntaxhighlight lang="modelica">
 +
function StiffnessMatrixGlobal
 +
 +
input Real member_mat [:,7];
 +
input Real node_mat [:,5];
 +
output Real kge_mat [size(member_mat,1),2*size(node_mat,1),2*size(node_mat,1)];
 +
 +
protected
 +
Real [size(member_mat,1),4,4] Ke_mat;
 +
Integer [2] nodes;
 +
Integer ke_row = 1;
 +
Integer ke_col = 1;
 +
Integer g_row;
 +
Integer g_col;
 +
 +
algorithm
 +
 +
// initialize element stiffness matrix and its global format
 +
Ke_mat := StiffnessMatrixElement(member_mat);
 +
kge_mat := zeros(size(member_mat,1),2*size(node_mat,1),2*size(node_mat,1));
 +
 +
// for each member element stiffness matrix, find the corresponding global matrix coordinates
 +
for member in 1:size(member_mat,1) loop
 +
  nodes:= integer({member_mat[member,node] for node in 2:3}); // nodes {i,j}
 +
 
 +
  ke_row := 1;
 +
  while ke_row <= size(Ke_mat,2) loop
 +
    g_row := (2 * nodes[div(ke_row + rem(ke_row,2),2)]) - rem(ke_row,2);
 +
   
 +
    ke_col := 1;
 +
    while ke_col <= size(Ke_mat,3) loop
 +
      g_col := (2 * nodes[div(ke_col + rem(ke_col,2),2)]) - rem(ke_col,2);
 +
     
 +
      // Transfers values from the element format into its position in the global format
 +
      kge_mat[member, g_row, g_col] := Ke_mat[member, ke_row, ke_col];
 +
     
 +
      ke_col := ke_col + 1;
 +
    end while;
 +
   
 +
    ke_row := ke_row + 1;
 +
  end while;
 +
 
 +
end for;
 +
 +
end StiffnessMatrixGlobal;
 +
 +
</syntaxhighlight>
 +
 +
|}
 +
 +
{| class=wikitable
 +
| style="vertical-align:top"|'''Function: Global Truss Stiffness Matrix'''
 +
<syntaxhighlight lang="modelica">
 +
function StiffnessMatrixFinal
 +
 +
input Real [:,:,size(kge_mat,2)] kge_mat;
 +
output Real [size(kge_mat,2),size(kge_mat,2)] kgfinal_mat;
 +
 +
algorithm
 +
 +
for row in 1:size(kgfinal_mat,1) loop
 +
  for col in 1:size(kgfinal_mat,2) loop
 +
  kgfinal_mat[row,col] := sum({kge_mat[i,row,col] for i in 1:size(kge_mat,1)});
 +
  end for;
 +
end for;
 +
 +
end StiffnessMatrixFinal;
 +
 +
</syntaxhighlight>
 +
 +
| style="vertical-align:top"|'''Function: Displacement (U) Matrix'''
 +
<syntaxhighlight lang="modelica">
 +
function DisplacementMatrix
 +
 +
input Real node_mat [:,5];
 +
input Real kgfinal_mat [2*size(node_mat,1),2*size(node_mat,1)];
 +
output Real displacement_mat [2*size(node_mat,1)];
 +
 +
protected
 +
Real F [2*size(node_mat,1)];
 +
Real aug_mat [2*size(node_mat,1),2*size(node_mat,1)+1];
 +
 +
// [K]{U} = {F}
 +
algorithm
 +
F := {node_mat[integer((n+rem(n,2))/2),integer(3-rem(n,2))] for n in 1:size(F,1)};
 +
 +
// Construct Augmented Matrix of [kgfinal_mat|F]
 +
for row in 1:size(aug_mat,1) loop
 +
  for col_kg in 1:size(kgfinal_mat,1) loop
 +
    aug_mat[row,col_kg] := kgfinal_mat[row,col_kg];
 +
  end for;
 +
 
 +
  aug_mat[row,end] := F[row];
 +
 
 +
end for;
 +
 +
// Changes rows according to boundary conditions
 +
for row in 1:size(aug_mat,1) loop
 +
  if node_mat[integer((row+rem(row,2))/2),integer(5-rem(row,2))] == 1 then
 +
    for col in 1:size(aug_mat,2)-1 loop
 +
      aug_mat[row,col] := 0;
 +
    end for;
 +
    aug_mat[row,row] := 1;
 +
  end if;
 +
end for;
 +
 +
// Solve into Reduced Row Echelon Form by Gauss Jordan function
 +
aug_mat := GaussJordan(aug_mat);
 +
 +
displacement_mat := {aug_mat[i,end] for i in 1:size(aug_mat,1)};
 +
 +
end DisplacementMatrix;
 +
 +
</syntaxhighlight>
 +
 +
|}
 +
 +
{| class=wikitable
 +
| style="vertical-align:top"|'''Function: Postprocessing Phase'''
 +
<syntaxhighlight lang="modelica">
 +
function TrussPostProcess
 +
 +
input Real member_mat [:,7];
 +
input Real node_mat [:,5];
 +
input Real kgfinal_mat [2*size(node_mat,1),2*size(node_mat,1)];
 +
input Real displacement_mat [2*size(node_mat,1)];
 +
 +
// Reaction Forces Vector
 +
output Real reaction_mat [2*size(node_mat,1)];
 +
 +
// Internal Forces and Normal Stresses Matrix
 +
output Real uf_mat [size(member_mat,1),7];
 +
 +
protected
 +
Real F [2*size(node_mat,1)];
 +
Real k_vec [size(member_mat,1)];
 +
Real theta;
 +
Real trig [2];
 +
Real T [4,4];
 +
Real mem_nodes [2];
 +
Real U_mat [4];
 +
Real local_u [4];
 +
 +
algorithm
 +
// Load Matrix
 +
F := {node_mat[integer((n+rem(n,2))/2),integer(3-rem(n,2))] for n in 1:size(F,1)};
 +
 +
// Stiffness Constant Vector
 +
k_vec := {(member_mat[i,5] * member_mat[i,6] / member_mat[i,7]) for i in 1:size(member_mat,1)};
 +
 +
// Reaction Forces
 +
reaction_mat := kgfinal_mat * displacement_mat - F;
 +
 +
// Cleaning float errors
 +
for force in 1:size(reaction_mat,1) loop
 +
  if abs(reaction_mat[force]) <= 10e-6 then
 +
    reaction_mat[force] := 0;
 +
  end if;
 +
end for;
 +
 +
// Internal Forces and Normal Stresses
 +
for row in 1:size(uf_mat,1) loop
 +
  uf_mat[row,1] := member_mat[row,1]; // Transfer member #
 +
 
 +
  // Construct transformation matrix
 +
  theta := Modelica.SIunits.Conversions.from_deg(member_mat[row,4]);
 +
  trig := {Modelica.Math.sin(theta), Modelica.Math.cos(theta)};
 +
  T := [trig[2], trig[1], 0, 0; -trig[1], trig[2], 0, 0;
 +
        0, 0, trig[2], trig[1]; 0, 0, -trig[1], trig[2]];
 +
 
 +
  // Retrieve U displacement values
 +
  mem_nodes := {member_mat[row,n] for n in 2:3}; // Nodes {i,j}
 +
  U_mat := {displacement_mat[integer(2*mem_nodes[integer((k+rem(k,2))/2)] - rem(k,2))] for k in 1:4};
 +
 
 +
  // Calculate internal forces u
 +
  local_u := T * U_mat;
 +
 
 +
  // Transfer values to uf matrix
 +
  for u_col in 2:5 loop
 +
    uf_mat[row,u_col] := local_u[u_col-1];
 +
  end for;
 +
 
 +
  // Calculate internal forces and normal stresses
 +
  uf_mat[row,6] := k_vec[row] * (uf_mat[row,2] - uf_mat[row,4]);
 +
  uf_mat[row,7] := uf_mat[row,6] / member_mat[row,5];
 +
 
 +
end for;
 +
 +
// Cleaning float errors
 +
for value_col in 1:size(uf_mat,1) loop
 +
  for value_row in 2:size(uf_mat,1) loop
 +
    if abs(uf_mat[value_row,value_col]) <= 10e-6 then
 +
      uf_mat[value_row,value_col] := 0;
 +
    end if;
 +
  end for;
 +
end for;
 +
 +
end TrussPostProcess;
 +
 +
</syntaxhighlight>
 +
 +
| style="vertical-align:top"|'''Function: Gauss Jordan Elimination'''
 +
<syntaxhighlight lang="modelica">
 +
function GaussJordan
 +
 +
// 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 [size(A,1),size(A,2)] 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 rows
 +
Real r;                // Ratio value used for row operations
 +
 +
// Limit to handle floating point errors
 +
Real float_error = 10e-6;
 +
 +
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
 +
    if sum({(B[h,i] - B[max_row,i]) for i in 1:n}) <> 0 then
 +
      for i in 1:n loop
 +
        B[h,i] := B[h,i] - B[max_row,i];
 +
        B[max_row,i] := B[h,i] + B[max_row,i];
 +
        B[h,i] := B[max_row,i] - B[h,i];
 +
      end for;
 +
    end if;
 +
   
 +
    // Divide pivot row by pivot number
 +
    r := B[h,k];
 +
    for i in 1:n loop
 +
      B[h,i] := B[h,i] / r;
 +
    end for;
 +
   
 +
    // 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
 +
    h := h-1;
 +
   
 +
  else
 +
   
 +
    // Perform row operations
 +
    for i in 1:(h-1) loop
 +
      r := B[i,k];
 +
      for g in 1:n loop
 +
      B[i,g] := B[i,g] - B[h,g] * 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 reduced row echelon form
 +
 +
end GaussJordan;
 +
 +
</syntaxhighlight>
 +
 +
|}
 +
 +
=== Week 04 [02/12/20] ===
 +
 +
==== Quiz: Diagarm Class & Flowchart ====
 +
Josiah Enrico membuat [[JosiahEnrico#Finite_Element_Method_for_Trusses_-_Metode_Numerik.2F2_Desember_2020|kode Finite Element Method Truss]] yang dapat digunakan dalam aplikasi Modelica untuk menyelesaikan perhitungan ''displacement'' dan ''reaction force'' pada trusses. Kode ini dapat diuraikan dalam bentuk ''Flowchart'' dan ''Diagram Class''. Berikut adalah contoh saya membuat ''flowchart'' untuk penyelesaian kasus pertama:
 +
<gallery mode="slideshow">
 +
File:TrussProblem1Flowchart.jpg|200px|model Trusses
 +
File:MatrixTransformFlowchart.jpg|200px|function Stiffness_Matrices
 +
</gallery>
 +
 +
==== 3D Truss Analysis using OpenModelica ====
 +
 +
[[File:Soal qius 2.2.jpeg|600px|thumb|center|Chapter 3 Problem 8 (''Finite Element Analysis, 4th Edition, Saeed Moaveni'')]]
 +
 +
===== Class: Solving 2D Truss Problem =====
 +
{| class=wikitable
 +
| style="vertical-align:top"|'''Declarations: Parameters and Variables'''
 +
<syntaxhighlight lang="modelica">
 +
// DECLARATIONS
 +
// Data for each member:  [element #, node i, node j, area, modulus]
 +
// SI Units:              [integer, integer, integer, meters^2, pascals]
 +
// Imperial Units:        [integer, integer, integer, inches^2, lb/inches]
 +
parameter Real [:,5] member = [1, 1, 4, 15e-4, 70e6;
 +
                              2, 1, 3, 15e-4, 70e6;
 +
                              3, 1, 2, 15e-4, 70e6];
 +
 +
// Coordinates, external loads and boundary conditions for each node
 +
// [node #, X, Y, Z, FX, FY, FZ, UX(1=Fixed), UY(1=Fixed), UZ(1=Fixed)]
 +
parameter Real [:,10] node_load = [1, 2,  0,    0, 0, -5000, 0, 0, 0, 0;
 +
                                  2, 0,  0, -1.5, 0,    0, 0, 1, 1, 1;
 +
                                  3, 0,  0,  1.5, 0,    0, 0, 1, 1, 1;
 +
                                  4, 0, 1.5,  0  , 0,    0, 0, 1, 1, 1];
 +
 +
// Arrays for Stiffness Matrices
 +
Real Ke [size(member,1), 6, 6]; // Element
 +
Real KGe [size(member,1), 3*size(node_load,1), 3*size(node_load,1)]; // Element in Global Form
 +
Real KGt [3*size(node_load,1), 3*size(node_load,1)]; // Total Global
 +
 +
// Displacement Matrix
 +
Real U [3*size(node_load,1)];
 +
 +
// Postprocessing Phase
 +
Real R [3*size(node_load,1)]; // Reaction forces for each node
 +
</syntaxhighlight>
 +
 +
| style="vertical-align:top"|'''Equations'''
 +
<syntaxhighlight lang="modelica">
 +
// EQUATIONS
 +
equation
 +
 +
// Find stiffness matrix for each element
 +
Ke = StiffnessMatrixElement(member);
 +
 +
// Assemble the global stiffness matrix
 +
KGe = StiffnessMatrixGlobal(member, node_load);
 +
KGt = StiffnessMatrixFinal(KGe);
 +
 +
// Solve for the displacement matrix
 +
U = DisplacementMatrix(node_load, KGt);
 +
 +
// Postprocessing Phase
 +
R = TrussPostProcess(member, node_load, KGt, U);
 +
 +
</syntaxhighlight>
 +
 +
|}
 +
 +
=== Week 05 [09/12/20] ===
 +
Pembahasan pemahaman algoritma Modelica untuk aplikasi metode numerik dalam kasus statika struktur statis Truss 2D (bidang) dan 3D (space).
 +
 +
==== Solving Trusses: 2D and 3D Flowcharts ====
 +
''(under construction)''
 +
 +
[[File:SolvingTrussesChristopherSE.mkv|720px|SolvingTrusses.mkv|center]]
 +
 +
=== Week 06 [16/12/20] ===
 +
 +
==== Fluid Systems ====
 +
System curve: a graphical representation of the pump head, and other properties of a fluid system.
 +
 +
<youtube>uP1ZiZ4khDM</youtube>
 +
 +
==== Golden Ratio Bracket Optimization ====
 +
Golden-section search: a technique for finding an extremum of a function inside a specified interval. Discovered in 1953 by Jack Kiefer, an American mathematical statistician.
 +
 +
<youtube>wpGN2xus75w</youtube>
 +
 +
=====One-Dimensional Unconstrained Optimization using OpenModelica=====
 +
 +
<syntaxhighlight lang="modelica">
 +
model bracket_optimization3
 +
 +
parameter Integer n=8;
 +
 +
Real x1[n];
 +
Real x2[n];
 +
Real xup;
 +
Real xlow;
 +
Real d;
 +
Real f1[n];
 +
Real f2[n];
 +
Real xopt;
 +
Real yopt;
 +
 +
algorithm
 +
xup :=4; xlow:=0;
 +
 +
for i in (1:n) loop
 +
  d:= (5^(1/2)-1)/2*(xup-xlow);
 +
  x1[i]:= xlow+d;
 +
  x2[i]:= xup-d;
 +
  f1[i]:= f_obj3(x1[i]);
 +
  f2[i]:= f_obj3(x2[i]);
 +
 
 +
  if f1[i]>f2[i] then
 +
    xup:= xup;
 +
    xlow:= x2[i];
 +
    xopt:= xup;
 +
    yopt:= f1[i];
 +
  else xlow:= xlow;
 +
    xup:= x1[i];
 +
    xopt:= xup;
 +
  end if;
 +
end for;
 +
 +
end bracket_optimization3;
 +
</syntaxhighlight>
 +
 +
=== Week 07 [23/12/20] ===
 +
 +
==== Powell's Method ====
 +
 +
Powell's method (Powell's conjugate direction method) is an algorithm proposed by Michael J. D. Powell for finding a local minimum of a function. The function need not be differentiable, and no derivatives are taken. The algorithm iterates an arbitrary number of times until no significant improvement is made. Powell's method can be used for multivariable optimization.
 +
 +
<youtube>1Z_4sBNoZj4</youtube>
 +
 +
=== Week 08 [30/12/20] ===
 +
 +
==== Pemodelan dan Optimisasi Rangka Batang ====
 +
 +
===== 3D Truss Analysis =====
 +
[[File:Tugas Besar Metnum Geometri Jos.jpg|600px|center]]
 +
 +
[[File:Captureaaaaaaaaaaaa.PNG|300px|center]]
 +
 +
====== Modelling: Truss Arrangement, Material Properties, Boundary Conditions ======
 +
[[File:TugasBesar3DTruss_1.png|600px]]
 +
 +
[[File:TugasBesar3DTruss_2.png|800px]]
 +
 +
====== Solving: Displacement, Reaction Forces, Check Static Sum ======
 +
[[File:TugasBesar3DTruss_3.png|800px]]
 +
 +
[[File:TugasBesar3DTruss_4.png|600px]]
 +
 +
[[File:TugasBesar3DTruss_5.png|500px]]
 +
 +
[[File:TugasBesar3DTruss_6.png|800px]]
 +
 +
====== PostProcessing: Stress & Safety Factor ======
 +
[[File:TugasBesar3DTruss_7.png|800px]]
 +
 +
===== Optimization =====
 +
 +
====== Data Handling & Curve Fitting ======
 +
'''Least-Squares Curve-Fitting'''
 +
 +
[[File:TugasBesar3DTruss_8.png|800px]]
 +
 +
'''Data Handling'''
 +
 +
[[File:TugasBesar3DTruss_9.png|800px]]
 +
 +
====== Golden Section Search Optimization ======
 +
'''Golden Section Search: Maximum Optimization'''
 +
 +
[[File:TugasBesar3DTruss_10.png|800px]]
 +
 +
[[File:TugasBesar3DTruss_11.png|800px]]
 +
 +
=== Week 09 [06/01/21] ===
 +
 +
'''Verifikasi Muhasabah'''
 +
 +
1. Value (adab/attitude/moral value seperti tingkat kerajinan (kehadiran, menulis (tugas wiki), ketekunan, kesungguhan dll.)
 +
 +
2. Pemahaman pengetahuan/ knowledge (konsep2 dan prinsip2 metoda numerik)
 +
 +
3. Keterampilan (skill dlm mengaplikasikan konsep dalam penyelesaian kasus (algoritma, pemrograman dll)
 +
 +
=== Week 10 [13/01/21] ===
 +
 +
== Jawaban UAS ==
 +
 +
=== Foto per Soal ===
 +
<gallery mode="slideshow">
 +
File:UAS_ChristopherSE_1.jpg|200px|Soal 1
 +
File:UAS_ChristopherSE_2.jpg|200px|Soal 1, 2
 +
File:UAS_ChristopherSE_3.jpg|200px|Soal 1, 2, 3
 +
File:UAS_ChristopherSE_4.jpg|200px|Soal 1, 2, 3, 4, 5, 6
 +
</gallery>
 +
 +
=== Pengerjaan Soal Final ===
 +
<gallery mode="slideshow">
 +
File:UAS_ChristopherSE_5.jpg|200px|Soal 1
 +
File:UAS_ChristopherSE_6.jpg|200px|Soal 2
 +
File:UAS_ChristopherSE_7.jpg|200px|Soal 3
 +
File:UAS_ChristopherSE_8.jpg|200px|Soal 4
 +
File:UAS_ChristopherSE_9.jpg|200px|Soal 5, 6, 7
 +
</gallery>
 +
 +
=== Code Input & Results ===
 +
<gallery mode="slideshow">
 +
File:UAS_ChristopherSE_10.jpg|200px|Code Input
 +
File:UAS_ChristopherSE_11.jpg|200px|Displacement
 +
File:UAS_ChristopherSE_12.jpg|200px|Reaction
 +
File:UAS_ChristopherSE_13.jpg|200px|Stress
 +
</gallery>
 +
 +
=== Complete Code ===
 +
<gallery mode="slideshow">
 +
File:UAS_ChristopherSE_14.jpg|200px|Line 01-52
 +
File:UAS_ChristopherSE_15.jpg|200px|Line 50-94
 +
File:UAS_ChristopherSE_16.jpg|200px|Line 95-146
 +
File:UAS_ChristopherSE_17.jpg|200px|Line 147-183
 +
</gallery>
  
 
== Student Data ==
 
== Student Data ==
  
 
NPM: 1806201056
 
NPM: 1806201056

Latest revision as of 21:44, 13 January 2021

Greetings! Welcome to my page :)

About Me

Freut mich!

I am a student at the University of Indonesia pursuing a Bachelor's in Mechanical Engineering. I am fascinated with the fields of aerospace, robotics, and systems design.

On-campus I participate in two extracurricular organizations: English debating and robotics. As a debater, I enjoy exploring multiple perspectives on relevant issues and deconstructing policies through critical argumentation. Although I consider myself still far from an expert in debating, I am thrilled to say that I have achieved the Novice Champion title in the JOVED 2019 debating tournament. Besides debating, I recently joined the robotics club's rocketry division. Unfortunately, since the pandemic happened right after, I haven't done much with robotics yet, but I am looking forward to learning more about it!

Some of my interests include programming and foreign languages. For programming, I have some experience with microcontrollers, Arduino, and Raspberry Pi. The programming language that I'm currently grappling with is Python, and I know a little bit of C++, Java, and HTML. Another hobby of mine is learning foreign languages, and I'm able to speak some German, Japanese, Spanish, and French. Sometimes I study several other languages when I'm feeling up to it.

Numerical Methods

This page was made for the Numerical Methods class that I'm taking this semester (Ganjil 2020/2021). So far up to midterms, we have learned the following topics:

  • Taylor Series Approximation
  • Bisection Method
  • Newton-Raphson Method
  • Secant Method
  • Systems of Non-linear Equations
  • Gauss Elimination
  • Least Squares Regression
  • Interpolation
  • Numerical Differentiation

Week 01 [11/11/20]

Tujuan Belajar Metode Numerik

  1. Memahami konsep dan prinsip dasar metode numerik.
  2. Mengerti aplikasi dan penerapan metode numerik.
  3. Mampu menerapkan pengetahuan metode numerik di persoalan teknik.
  4. Mendapat nilai tambah sendiri agar menjadi orang yang lebih beradab.

Belajar OpenModelica

Di video ini, saya berusaha untuk menjelaskan cara basic menggunakan OpenModelica. Saya mencoba beberapa fungsi dasar seperti modelling di text view, modelling di diagram view, dan melakukan simulasi dan plotting. Contoh simulasi sistem yang digunakan adalah:

  • Pemodelan sistem fluida dalam tangki berbentuk bola dengan menulis code
  • Pemodelan sistem osilasi massa spring-damper dengan merancang komponen visual

Week 02 [18/11/20]

OpenModelica Intro

OpenModelica adalah aplikasi yang digunakan untuk melakukan pemodelan. OpenModelica memudahkan penyelesaian permasalahan permodelan sistem yang cukup rumit. Coding yang dilakukan dalam OMEdit ditulis berdasarkan bahasa C++.

OpenModelica digunakan untuk kelas Metode Numerik untuk beberapa alasan:

  • OpenModelica memampukan kita untuk melakukan pemodelan sitem teknik dengan mudah tanpa terlalu mendalami pengetahuan pemograman dan ilmu komputer.
  • OpenModelica adalah software yang Open Source, dengan akses bebas dan terbuka. Karena itu, terdapat banyak pengguna dan software dapat digunakan gratis dengan bebas secara legal.
Class & Function
Pilihan Specialization ketika membuat Modelica Class yang baru.

Ketika membuat file Class baru di OpenModelica, ada pilihan Specialization dimana jenis Class dapat dipilih. Setiap Class mempunyai karakteristik dan tujuan masing-masing. Informasi lebih lanjut tentang setiap jenis Class bisa dibaca dalam dokumentasi disini.

Salah satu jenis Class yang sudah digunakan adalah Model yang dapat melakukan pemodelan dan simulasi. Ada jenis Class lain yaitu Function dimana code yang diketik di dalamnya akan menjadi suatu fungsi yang dapat dipanggil dalam model. Ketika dibuat, Function mempunyai bagian algorithm dimana perhitungan fungsi dimasukkan. Sebelum bagian algorithm, dapat dimasukkan input dan output yang akan ditangani oleh fungsinya.

Sebagai contoh, dibikin suatu fungsi PlusTen yang menambahkan angka 10 ke nilai input x. Fungsi ini terus akan dipanggil dalam model TestPlusTen yang menerima parameter num dan menghasilkan variabel answer. Ketika dicompile, hasil penggunaan fungsi di model ini dapat dilihat.

  • Contoh fungsi yang menambahkan 10 ke nilai input.

  • Fungsi dipanggil dalam model.

  • Hasil model yang menggunakan fungsi.

Persamaan Aljabar Simultan

OpenModelica dapat melakukan perhitungan aljabar linier. Matriks ditulis sebagai Array yang dapat dilakukan kalkulasi seperti perkalian, transpose, inverse, dll. Informasi lebih detail tentang array dapat dibaca di sini dan di sini. Dokumentasi operator dan fungsi matriks dalam OpenModelica dabat dibaca di sini dan di sini.

Salah satu fungsi di OpenModelica adalah solve(). Fungsi ini dapat menyelesaikan suatu sistem persamaan linier A*x=B dan mencari solusi x dengan metode Gaussian elimination dengan partial pivoting. Dokumentasi tentang fungsi ini dapat dibaca di sini. Berikut adalah video yang menjelaskan tentang fungsi ini dan contoh pemakaiannya.

Gauss-Jordan

Berikut adalah contoh kode fungsi Modelica yang melakukan operasi Gauss-Jordan matriks menjadi bentuk reduced row echelon. 

function GaussJordan

// 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
    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;

Week 03 [25/11/20]

Metode Numerik dan Masalah Teknik

Aplikasi metode numerik dalam permasalahan teknik.

Masalah teknik dapat diselesaikan dengan menggunakan metode numerik. Beberapa contoh metode numerik yang sering digunakan untuk melakukan analisa masalah teknik adalah Computation Fluid Dynamics (CFD), Finite Element Analysis (FEA), dan Metode Stokastik. CFD dan FEA berdasarkan hukum-hukum fisika, sementara metode stokastik berbasis data dan statistik. Salah satu contoh metode stokastik yang sering digunakan adalah metode Monte Carlo.

Untuk menyelesaikan suatu masalah teknik menggunakan metode numerik, perlu dilakukan beberapa tahap. Pertama, masalah teknik perlu dilakukan analisis. Setelah itu, hasil analisis masalah perlu dilakukan pemodelan matematis, yang kemudian diterjemahkan menjadi model numerik yang dapat dihitung oleh sebuah komputer. Dengan melakukan perhitungan di software komputer, dapat diperolehkan solusi dari masalah teknik.

Truss: Finite Element Formulation

Contoh masalah statically determinate dan statically indeterminate.

Truss adalah struktur yang terdiri dari beam lurus yang saling dihubungkan di ujungnya. Masalah teknik yang terkait dengan truss dapat dibagikan menjadi statically determinate atau statically indeterminate. Pada masalah jenis statically determinate dilakukan asumsi bahwa truss yang dianalisis merupakan rigid body dan defleksi tidak diperhitungkan. Masalah jenis ini dapat diselesaikan dengan methods of joints or sections yang pernah dibahaskan di materi teknik yang lebih dasar. Masalah jenis statically indeterminate tidak dapat dilakukan analis dengan asumsi rigid body, sehingga diperlukan finite element method untuk menyelesaikannya.

Equivalent Stiffness: Average Stress, Average Strain, & Hooke's Law

(under construction)

Displacements & Forces: Global & Local Coordinates

(under construction)

Internal Forces & Stiffness Matrix

(under construction)

Naive Gauss Algorithm

Berikut adalah contoh pseudocode fungsi yang menggunakan algoritma Naive Gauss Elimination. Fungsi ini dapat digunakan untuk menyelesaikan sistem persamaan linear di mana rank matriks sama dengan jumlah baris. 

function NaiveGauss

/*
// Pseudocode Figure 9.4

//Forward Elimination
DOFOR k = 1, n - 1
  DOFOR i = k + 1, n
    factor = a[i,k] / a[k,k]
    DOFOR j = k + 1 to n
      a[i,j] = a[i,j] - factor * a[k,j]
    END DO
    b[i] = b[i] - factor * b[k]
  END DO
END DO

// Back Substitution
x[n] = b[n] / a[n,n]
DOFOR i = n - 1, 1, -1
  sum = b[i]
  DOFOR j = i + 1, n
    sum = sum = a[i,j] * x[j]
  END DO
  x[i] = sum / a[i,i]
END DO
*/

input Real [:,:] A; // Left-hand Coefficients of the Linear Equation System in array form
input Real [:] B;   // Right-hand Constants of the Linear Equation System in array form
output Real [:,:] y;// Array containing coefficient matrix after NGE operation
output Real [:] x;  // Array containing solved values of x

protected
Real [:,:] a;
Real [:] b;
Integer m = size(A,1);  // Number of rows in matrix
Integer n = size(A,2);  // Number of columns in matrix
Real k = 1;       // Pivot column pointer
Real i = 1;       // Row counter
Real j = 1;       // Row element counter
Real factor = 1;  // Factor value used for forward elimination
Real sum = 1;     // Sum value used for back substitution

algorithm

// Transfer input matrix (A,B) into variables (a,b)
a := A;
b := B;

// Forward Elimination
for k in 1:(n-1) loop
  for i in (k+1):n loop
    factor := a[i,k] / a[k,k];
    for j in (k+1):n loop
      a[i,j] := a[i,j] - (factor * a[k,j]);
    end for;
    b[i] := b[i] - (factor * b[k]);
  end for;
end for;

// Back Substitution
x[n] := b[n] / a[n,n];
for i in (n-1):(-1) loop
  sum := b[i];
  for j in (i+1):n loop
    sum := sum - (a[i,j] * x[j]);
  end for;
  x[i] := sum / a[i,i];
end for;

end NaiveGauss;

2D Truss Analysis using OpenModelica

Chapter 3 Problem 4 (Finite Element Analysis, 4th Edition, Saeed Moaveni)
Class: Solving 2D Truss Problem
Declarations: Parameters and Variables 
// DECLARATIONS
// Data for each member: [element #, node i, node j, theta, area, modulus, length]
// SI Units: [integer, integer, integer, degrees, meters^2, pascals, meters]
// Imperial Units: [integer, integer, integer, degrees, inches^2, lb/inches, inches]
parameter Real [:,7] member = [1, 1, 2,      0, 10e-4, 200e9, 1;
                               2, 2, 3,      0, 10e-4, 200e9, 1;
                               3, 1, 4, 308.66, 10e-4, 200e9, 1.60;
                               4, 2, 4, 270.00, 10e-4, 200e9, 1.25;
                               5, 3, 4, 231.34, 10e-4, 200e9, 1.60];

// External loads and boundary conditions for each node: [node #, FX, FY, UX(1:Fixed), UY(1:Fixed)]
parameter Real [:,5] node_load = [1,        0,        0, 1, 1;
                                  2, -1035.28, -3863.70, 0, 0;
                                  3,        0,        0, 1, 1;
                                  4, -1035.28, -3863.70, 0, 0];

// Vector for equivalent stiffness constant of each member
Real k [size(member,1)];

// Arrays for Stiffness Matrices
Real Ke [size(member,1), 4, 4]; // Element
Real KGe [size(member,1), 2*size(node_load,1), 2*size(node_load,1)]; // Element in Global Form
Real KGt [2*size(node_load,1), 2*size(node_load,1)]; // Total Global

// Displacement Matrix
Real U [2*size(node_load,1)]; 

// Postprocessing Phase
Real R [2*size(node_load,1)]; // Reaction forces for each node
Real InForce_NormStress [size(member,1),7]; // Internal forces and normal stress for each member
// Format: [member #, uix, uiy, ujx, ujy, internal force, normal stress]
 Equations 
// EQUATIONS
equation

// Calculate stiffness constants for each member element
k = {(member[i,5] * member[i,6] / member[i,7]) for i in 1:size(member,1)};

// Find stiffness matrix for each element
Ke = StiffnessMatrixElement(member);

// Assemble the global stiffness matrix
KGe = StiffnessMatrixGlobal(member, node_load);
KGt = StiffnessMatrixFinal(KGe);

// Solve for the displacement matrix
U = DisplacementMatrix(node_load, KGt);

// Postprocessing Phase
(R, InForce_NormStress) = TrussPostProcess(member, node_load, KGt, U);
Functions: Solving 2D Truss Problem
Function: Local Element Stiffness Matrix 
function StiffnessMatrixElement

input Real [:,7] member_mat;
output Real [size(member_mat,1),4,4] Ke_mat;

protected
Real theta;
Real [3] StiffTrig;
Real [4,4] StiffTrans;
Real [size(member_mat,1)] k_vec;
Real float_error = 10e-10;

algorithm

k_vec := {(member_mat[i,5] * member_mat[i,6] / member_mat[i,7]) for i in 1:size(member_mat,1)};

// Finding stiffness matrix of each element member
for i in 1:size(member_mat,1) loop

  // Clearing the matrices
  StiffTrig := zeros(3);
  StiffTrans := zeros(4,4);
  
  // Converting degrees to radians
  theta := Modelica.SIunits.Conversions.from_deg(member_mat[i,4]);

  // {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_mat
  for m in 1:4 loop
    for n in 1:4 loop
      Ke_mat[i,m,n] := k_vec[i] * StiffTrans[m,n];
    end for;
  end for;

end for;

end StiffnessMatrixElement;
 Function: Global Element Stiffness Matrix 
function StiffnessMatrixGlobal

input Real member_mat [:,7];
input Real node_mat [:,5];
output Real kge_mat [size(member_mat,1),2*size(node_mat,1),2*size(node_mat,1)];

protected
Real [size(member_mat,1),4,4] Ke_mat;
Integer [2] nodes;
Integer ke_row = 1;
Integer ke_col = 1;
Integer g_row;
Integer g_col;

algorithm

// initialize element stiffness matrix and its global format
Ke_mat := StiffnessMatrixElement(member_mat);
kge_mat := zeros(size(member_mat,1),2*size(node_mat,1),2*size(node_mat,1));

// for each member element stiffness matrix, find the corresponding global matrix coordinates
for member in 1:size(member_mat,1) loop
  nodes:= integer({member_mat[member,node] for node in 2:3}); // nodes {i,j}
  
  ke_row := 1;
  while ke_row <= size(Ke_mat,2) loop
    g_row := (2 * nodes[div(ke_row + rem(ke_row,2),2)]) - rem(ke_row,2);
    
    ke_col := 1;
    while ke_col <= size(Ke_mat,3) loop
      g_col := (2 * nodes[div(ke_col + rem(ke_col,2),2)]) - rem(ke_col,2);
      
      // Transfers values from the element format into its position in the global format
      kge_mat[member, g_row, g_col] := Ke_mat[member, ke_row, ke_col];
      
      ke_col := ke_col + 1;
    end while;
    
    ke_row := ke_row + 1;
  end while;
  
end for;

end StiffnessMatrixGlobal;
Function: Global Truss Stiffness Matrix 
function StiffnessMatrixFinal

input Real [:,:,size(kge_mat,2)] kge_mat;
output Real [size(kge_mat,2),size(kge_mat,2)] kgfinal_mat;

algorithm

for row in 1:size(kgfinal_mat,1) loop
  for col in 1:size(kgfinal_mat,2) loop
  kgfinal_mat[row,col] := sum({kge_mat[i,row,col] for i in 1:size(kge_mat,1)});
  end for;
end for;

end StiffnessMatrixFinal;
 Function: Displacement (U) Matrix 
function DisplacementMatrix

input Real node_mat [:,5];
input Real kgfinal_mat [2*size(node_mat,1),2*size(node_mat,1)];
output Real displacement_mat [2*size(node_mat,1)];

protected
Real F [2*size(node_mat,1)];
Real aug_mat [2*size(node_mat,1),2*size(node_mat,1)+1];

// [K]{U} = {F}
algorithm
F := {node_mat[integer((n+rem(n,2))/2),integer(3-rem(n,2))] for n in 1:size(F,1)};

// Construct Augmented Matrix of [kgfinal_mat|F]
for row in 1:size(aug_mat,1) loop
  for col_kg in 1:size(kgfinal_mat,1) loop
    aug_mat[row,col_kg] := kgfinal_mat[row,col_kg];
  end for;
  
  aug_mat[row,end] := F[row];
  
end for;

// Changes rows according to boundary conditions
for row in 1:size(aug_mat,1) loop
  if node_mat[integer((row+rem(row,2))/2),integer(5-rem(row,2))] == 1 then
    for col in 1:size(aug_mat,2)-1 loop
      aug_mat[row,col] := 0;
    end for;
    aug_mat[row,row] := 1;
  end if;
end for;

// Solve into Reduced Row Echelon Form by Gauss Jordan function
aug_mat := GaussJordan(aug_mat);

displacement_mat := {aug_mat[i,end] for i in 1:size(aug_mat,1)};

end DisplacementMatrix;
Function: Postprocessing Phase 
function TrussPostProcess

input Real member_mat [:,7];
input Real node_mat [:,5];
input Real kgfinal_mat [2*size(node_mat,1),2*size(node_mat,1)];
input Real displacement_mat [2*size(node_mat,1)];

// Reaction Forces Vector
output Real reaction_mat [2*size(node_mat,1)];

// Internal Forces and Normal Stresses Matrix
output Real uf_mat [size(member_mat,1),7];

protected
Real F [2*size(node_mat,1)];
Real k_vec [size(member_mat,1)];
Real theta;
Real trig [2];
Real T [4,4];
Real mem_nodes [2];
Real U_mat [4];
Real local_u [4];

algorithm
// Load Matrix
F := {node_mat[integer((n+rem(n,2))/2),integer(3-rem(n,2))] for n in 1:size(F,1)};

// Stiffness Constant Vector
k_vec := {(member_mat[i,5] * member_mat[i,6] / member_mat[i,7]) for i in 1:size(member_mat,1)};

// Reaction Forces
reaction_mat := kgfinal_mat * displacement_mat - F;

// Cleaning float errors
for force in 1:size(reaction_mat,1) loop
  if abs(reaction_mat[force]) <= 10e-6 then
    reaction_mat[force] := 0;
  end if;
end for;

// Internal Forces and Normal Stresses
for row in 1:size(uf_mat,1) loop
  uf_mat[row,1] := member_mat[row,1]; // Transfer member #
  
  // Construct transformation matrix
  theta := Modelica.SIunits.Conversions.from_deg(member_mat[row,4]);
  trig := {Modelica.Math.sin(theta), Modelica.Math.cos(theta)};
  T := [trig[2], trig[1], 0, 0; -trig[1], trig[2], 0, 0;
        0, 0, trig[2], trig[1]; 0, 0, -trig[1], trig[2]];
  
  // Retrieve U displacement values
  mem_nodes := {member_mat[row,n] for n in 2:3}; // Nodes {i,j}
  U_mat := {displacement_mat[integer(2*mem_nodes[integer((k+rem(k,2))/2)] - rem(k,2))] for k in 1:4};
  
  // Calculate internal forces u
  local_u := T * U_mat;
  
  // Transfer values to uf matrix
  for u_col in 2:5 loop
    uf_mat[row,u_col] := local_u[u_col-1];
  end for;
  
  // Calculate internal forces and normal stresses
  uf_mat[row,6] := k_vec[row] * (uf_mat[row,2] - uf_mat[row,4]);
  uf_mat[row,7] := uf_mat[row,6] / member_mat[row,5];
  
end for;

// Cleaning float errors
for value_col in 1:size(uf_mat,1) loop
  for value_row in 2:size(uf_mat,1) loop
    if abs(uf_mat[value_row,value_col]) <= 10e-6 then
      uf_mat[value_row,value_col] := 0;
    end if;
  end for;
end for;

end TrussPostProcess;
 Function: Gauss Jordan Elimination 
function GaussJordan

// 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 [size(A,1),size(A,2)] 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 rows
Real r;                 // Ratio value used for row operations

// Limit to handle floating point errors
Real float_error = 10e-6;

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
    if sum({(B[h,i] - B[max_row,i]) for i in 1:n}) <> 0 then
      for i in 1:n loop
        B[h,i] := B[h,i] - B[max_row,i];
        B[max_row,i] := B[h,i] + B[max_row,i];
        B[h,i] := B[max_row,i] - B[h,i];
      end for;
    end if;
    
    // Divide pivot row by pivot number
    r := B[h,k];
    for i in 1:n loop
      B[h,i] := B[h,i] / r;
    end for;
    
    // 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
    h := h-1;
    
  else
    
    // Perform row operations
    for i in 1:(h-1) loop
      r := B[i,k];
      for g in 1:n loop
      B[i,g] := B[i,g] - B[h,g] * 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 reduced row echelon form

end GaussJordan;

Week 04 [02/12/20]

Quiz: Diagarm Class & Flowchart

Josiah Enrico membuat kode Finite Element Method Truss yang dapat digunakan dalam aplikasi Modelica untuk menyelesaikan perhitungan displacement dan reaction force pada trusses. Kode ini dapat diuraikan dalam bentuk Flowchart dan Diagram Class. Berikut adalah contoh saya membuat flowchart untuk penyelesaian kasus pertama:

  • model Trusses

  • function Stiffness_Matrices

3D Truss Analysis using OpenModelica

Chapter 3 Problem 8 (Finite Element Analysis, 4th Edition, Saeed Moaveni)
Class: Solving 2D Truss Problem
Declarations: Parameters and Variables 
// DECLARATIONS
// Data for each member:  [element #, node i, node j, area, modulus]
// SI Units:              [integer, integer, integer, meters^2, pascals]
// Imperial Units:        [integer, integer, integer, inches^2, lb/inches]
parameter Real [:,5] member = [1, 1, 4, 15e-4, 70e6;
                               2, 1, 3, 15e-4, 70e6;
                               3, 1, 2, 15e-4, 70e6];

// Coordinates, external loads and boundary conditions for each node
// [node #, X, Y, Z, FX, FY, FZ, UX(1=Fixed), UY(1=Fixed), UZ(1=Fixed)]
parameter Real [:,10] node_load = [1, 2,   0,    0, 0, -5000, 0, 0, 0, 0;
                                   2, 0,   0, -1.5, 0,     0, 0, 1, 1, 1;
                                   3, 0,   0,  1.5, 0,     0, 0, 1, 1, 1;
                                   4, 0, 1.5,  0  , 0,     0, 0, 1, 1, 1];

// Arrays for Stiffness Matrices
Real Ke [size(member,1), 6, 6]; // Element
Real KGe [size(member,1), 3*size(node_load,1), 3*size(node_load,1)]; // Element in Global Form
Real KGt [3*size(node_load,1), 3*size(node_load,1)]; // Total Global

// Displacement Matrix
Real U [3*size(node_load,1)];

// Postprocessing Phase
Real R [3*size(node_load,1)]; // Reaction forces for each node
 Equations 
// EQUATIONS
equation

// Find stiffness matrix for each element
Ke = StiffnessMatrixElement(member);

// Assemble the global stiffness matrix
KGe = StiffnessMatrixGlobal(member, node_load);
KGt = StiffnessMatrixFinal(KGe);

// Solve for the displacement matrix
U = DisplacementMatrix(node_load, KGt);

// Postprocessing Phase
R = TrussPostProcess(member, node_load, KGt, U);

Week 05 [09/12/20]

Pembahasan pemahaman algoritma Modelica untuk aplikasi metode numerik dalam kasus statika struktur statis Truss 2D (bidang) dan 3D (space).

Solving Trusses: 2D and 3D Flowcharts

(under construction)

SolvingTrusses.mkv

Week 06 [16/12/20]

Fluid Systems

System curve: a graphical representation of the pump head, and other properties of a fluid system.

Golden Ratio Bracket Optimization

Golden-section search: a technique for finding an extremum of a function inside a specified interval. Discovered in 1953 by Jack Kiefer, an American mathematical statistician.

One-Dimensional Unconstrained Optimization using OpenModelica
model bracket_optimization3

parameter Integer n=8;

Real x1[n]; 
Real x2[n]; 
Real xup; 
Real xlow; 
Real d; 
Real f1[n];
Real f2[n];
Real xopt; 
Real yopt;

algorithm
xup :=4; xlow:=0;

for i in (1:n) loop 
  d:= (5^(1/2)-1)/2*(xup-xlow); 
  x1[i]:= xlow+d; 
  x2[i]:= xup-d; 
  f1[i]:= f_obj3(x1[i]);
  f2[i]:= f_obj3(x2[i]);
  
  if f1[i]>f2[i] then
    xup:= xup;
    xlow:= x2[i];
    xopt:= xup; 
    yopt:= f1[i];
  else xlow:= xlow; 
    xup:= x1[i];
    xopt:= xup;
  end if;
end for;

end bracket_optimization3;

Week 07 [23/12/20]

Powell's Method

Powell's method (Powell's conjugate direction method) is an algorithm proposed by Michael J. D. Powell for finding a local minimum of a function. The function need not be differentiable, and no derivatives are taken. The algorithm iterates an arbitrary number of times until no significant improvement is made. Powell's method can be used for multivariable optimization.

Week 08 [30/12/20]

Pemodelan dan Optimisasi Rangka Batang

3D Truss Analysis
Tugas Besar Metnum Geometri Jos.jpg
Captureaaaaaaaaaaaa.PNG
Modelling: Truss Arrangement, Material Properties, Boundary Conditions

TugasBesar3DTruss 1.png

TugasBesar3DTruss 2.png

Solving: Displacement, Reaction Forces, Check Static Sum

TugasBesar3DTruss 3.png

TugasBesar3DTruss 4.png

TugasBesar3DTruss 5.png

TugasBesar3DTruss 6.png

PostProcessing: Stress & Safety Factor

TugasBesar3DTruss 7.png

Optimization
Data Handling & Curve Fitting

Least-Squares Curve-Fitting

TugasBesar3DTruss 8.png

Data Handling

TugasBesar3DTruss 9.png

Golden Section Search Optimization

Golden Section Search: Maximum Optimization

TugasBesar3DTruss 10.png

TugasBesar3DTruss 11.png

Week 09 [06/01/21]

Verifikasi Muhasabah

1. Value (adab/attitude/moral value seperti tingkat kerajinan (kehadiran, menulis (tugas wiki), ketekunan, kesungguhan dll.)

2. Pemahaman pengetahuan/ knowledge (konsep2 dan prinsip2 metoda numerik)

3. Keterampilan (skill dlm mengaplikasikan konsep dalam penyelesaian kasus (algoritma, pemrograman dll)

Week 10 [13/01/21]

Jawaban UAS

Foto per Soal

  • Soal 1

  • Soal 1, 2

  • Soal 1, 2, 3

  • Soal 1, 2, 3, 4, 5, 6

Pengerjaan Soal Final

  • Soal 1

  • Soal 2

  • Soal 3

  • Soal 4

  • Soal 5, 6, 7

Code Input & Results

  • Code Input

  • Displacement

  • Reaction

  • Stress

Complete Code

  • Line 01-52

  • Line 50-94

  • Line 95-146

  • Line 147-183

Student Data

NPM: 1806201056

Retrieved from "http://air.eng.ui.ac.id/index.php?title=Christopher_Salendu_Erwin&oldid=55402"