Difference between revisions of "Christopher Salendu Erwin"

From ccitonlinewiki
Jump to: navigation, search
(Modelling: Truss Arrangement, Material Properties, Boundary Conditions)
(Modelling: Truss Arrangement, Material Properties, Boundary Conditions)
Line 966: Line 966:
  
 
====== Modelling: Truss Arrangement, Material Properties, Boundary Conditions ======
 
====== Modelling: Truss Arrangement, Material Properties, Boundary Conditions ======
[[File:TugasBesar3DTruss_1.png|600px|center]][[File:TugasBesar3DTruss_2.png|600px|center]]
+
[[File:TugasBesar3DTruss_1.png|600px|center]][[File:TugasBesar3DTruss_2.png|800px|center]]
  
 
====== Solving: Displacement, Reaction Forces, Check Static Sum ======
 
====== Solving: Displacement, Reaction Forces, Check Static Sum ======

Revision as of 14:34, 6 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.

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:

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)

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
Modelling: Truss Arrangement, Material Properties, Boundary Conditions
TugasBesar3DTruss 1.png
TugasBesar3DTruss 2.png
Solving: Displacement, Reaction Forces, Check Static Sum
PostProcessing: Stress & Safety Factor

Week 09 [06/01/21]

Student Data

NPM: 1806201056