Tugas Besar Metode Numerik - Josiah Enrico S

Objective

• Membuat program analisis truss dan optimasi sederhana dengan Modelica

• Mengoptimasi harga pembuatan rangka truss sederhana dengan memvariasi dimensi dan elastisitas material.

Methodology

• Metode Optimasi Golden Section (univariabel)

• Metode Optimasi Golden Section (multivariabel)

Procedure

Geometri dan Load

Force

F1	2000
F2	1000

• 

  Figure 1 Truss geometry and load

Constraint:

- Spesifikasi L (Panjang) dan geometri rangka truss

- Gaya beban terhadap struktur (1000 N dan 2000 N)

Asumsi:

- Beban akan terdistribusi hanya pada point penghubung dan semua gaya yang dipengaruhi moment diamggap tidak ada karena sistem bersifat truss

Koleksi Data

• 

  Table 1 Angle bar standard dimension and area

• 

  Table 2 Cost of material in market

Data Proceessing

Truss Structure Analysis Model

Pada tahap pertama, dilakukanka simulasi menggunakan Modelica untuk mengetahui pengaruh beban gaya yang diberikan terhadap struktur truss. Pengaruh ini direpresantasikan oleh 4 nilai yaitu displacement/perpindahan, gaya reaksi, tekanan di setiap truss dan faktor keamanan di setiap truss. Keempat nilai ini diperoleh dengan menggunakan operasi tranformasi matriks untuk pembebanan truss dalam Finite Element Analysis. Operasi ini diterjemahkan ke dalam bahasa C dalam pemrograman Modelica. Berikut ini adalah hasil simulasi Modelica dan program tertulisnya:

• 

  Figure 2 Modelica Simulation Result of point displacement, reaction force, each truss stress, each truss safety factor

model Truss_3D_Tugas_Besar_Safety // Base Code by Josiah Enrico S. 1906356286 // Editing and Testing by Ahmad M. Fahmi 1806181836 // Commenting and Cleaning by Christopher S. Erwin 1806201056 // === DEFINITION PHASE === // Define the properties of the truss arrangement parameter Integer Nodes = size(node_points,1); // Number of Nodes parameter Integer Members = size(node_pairs,1); // Number of Members // Define the properties of the truss material parameter Real Yield = 215e6; // Yield Strength (Pa) parameter Real Area = 0.000224; // Area L Profile [Dimension=0.03, Thickness=0,004] (m^2) parameter Real Elas = 193e9; // Elasticity of material [SS 304] (Pa) // Define connections; list pair of nodes joined by a member parameter Integer node_pairs[:,2]=[1,5; 2,6; 3,7; 4,8; // Lowest vertical members 5,6; 6,7; 7,8; 5,8; // 1st shelf 5,9; 6,10; 7,11; 8,12; // Middle vertical members 9,10; 10,11; 11,12; 9,12; // 2nd shelf 9,13; 10,14; 11,15; 12,16; // Top vertical members 13,14; 14,15; 15,16; 13,16];// 3rd shelf // Define coordinates and boundaries of each node // [x, y, z, boundary_x, boundary_y, boundary_z] // If the node is bound (stationary) along that axis, then boundary_axis = 1 // otherwise the node is free and boundary_axis = 0 parameter Real node_points[:,6]=[0.3,-0.375,0,1,1,1; // Node 1 -0.3,-0.375,0,1,1,1; // Node 2 -0.3,0.375,0,1,1,1; // Node 3 0.3,0.375,0,1,1,1; // Node 4 0.3,-0.375,0.6,0,0,0; // Node 5 -0.3,-0.375,0.6,0,0,0; // Node 6 -0.3,0.375,0.6,0,0,0; // Node 7 0.3,0.375,0.6,0,0,0; // Node 8 0.3,-0.375,1.2,0,0,0; // Node 9 -0.3,-0.375,1.2,0,0,0; // Node 10 -0.3,0.375,1.2,0,0,0; // Node 11 0.3,0.375,1.2,0,0,0; // Node 12 0.3,-0.375,1.8,0,0,0; // Node 13 -0.3,-0.375,1.8,0,0,0; // Node 14 -0.3,0.375,1.8,0,0,0; // Node 15 0.3,0.375,1.8,0,0,0]; // Node 16 // Define external forces in Newtons on each node // [Fx, Fy, Fz] parameter Real Forces[nodes_xyz]={0,0,0, // Node 1 0,0,0, // Node 2 0,0,0, // Node 3 0,0,0, // Node 4 0,0,0, // Node 5 0,0,0, // Node 6 0,0,0, // Node 7 0,0,0, // Node 8 0,0,0, // Node 9 0,0,0, // Node 10 0,0,0, // Node 11 0,0,0, // Node 12 0,0,-500, // Node 13 0,0,-1000, // Node 14 0,0,-1000, // Node 15 0,0,-500}; // Node 16 // === SOLUTION PHASE === Real displacement[nodes_xyz], reaction[nodes_xyz]; Real check[3]; // To check if forces in xyz directions sum to zero (structure is static) Real safety[Members]; // Safety Factor of each member (based on stress experienced) Real mem_dis[3]; // Displacement in xyz directions experienced by member Real stress_vector[3]; // Stress experienced by member in xyz directions Real stress_mag[Members]; // Magnitude of stress experienced by member protected parameter Integer nodes_xyz = 3 * Nodes; // Size of vector of {node_#_x, node_#_y, node_#_z} containing all nodes Real node_i[3], node_j[3]; // Coordinates in xyz for member's node pairs // Global Stiffness Matrix construction Real G_element[nodes_xyz, nodes_xyz], G_total[nodes_xyz, nodes_xyz], G_total_new[nodes_xyz, nodes_xyz]; // Transformation Matrix construction Real cos_x, cos_y, cos_z, L, T[3,3]; // Limits to eliminate Float Errors Real err_solve = 10e-10, err_check = 10e-4; algorithm //Creating Global Matrix G_total := identity(nodes_xyz); for mem in 1:Members loop // Get coordinates of member's nodes for i in 1:3 loop node_i[i] := node_points[node_pairs[mem,1],i]; node_j[i] := node_points[node_pairs[mem,2],i]; end for; // Finding the Stiffness Matrix [Ke] // {F} = [Ke]{U} // [Ke] = k[T] L := Modelica.Math.Vectors.length(node_j - node_i); cos_x := (node_j[1]-node_i[1])/L; cos_y := (node_j[2]-node_i[2])/L; cos_z := (node_j[3]-node_i[3])/L; T := (Area * Elas / L)*[cos_x^2,cos_x*cos_y,cos_x*cos_z; cos_y*cos_x,cos_y^2,cos_y*cos_z; cos_z*cos_x,cos_z*cos_y,cos_z^2]; // Transforming to Global Matrix G_element := zeros(nodes_xyz,nodes_xyz); for m,n in 1:3 loop G_element[3*(node_pairs[mem,1]-1)+m, 3*(node_pairs[mem,1]-1)+n] := T[m,n]; G_element[3*(node_pairs[mem,2]-1)+m, 3*(node_pairs[mem,2]-1)+n] := T[m,n]; G_element[3*(node_pairs[mem,2]-1)+m, 3*(node_pairs[mem,1]-1)+n] := -T[m,n]; G_element[3*(node_pairs[mem,1]-1)+m, 3*(node_pairs[mem,2]-1)+n] := -T[m,n]; end for; // Adding to total Global Matrix G_total_new := G_total + G_element; G_total := G_total_new; end for; //Implementing boundary conditions to Global Matrix for nod in 1:Nodes loop // Boundary x-axis if node_points[nod,4] <> 0 then for i in 1:Nodes*3 loop G_total[(nod*3)-2,i] := 0; G_total[(nod*3)-2,(nod*3)-2] := 1; end for; end if; // Boundary y-axis if node_points[nod,5] <> 0 then for i in 1:Nodes*3 loop G_total[(nod*3)-1,i] := 0; G_total[(nod*3)-1,(nod*3)-1] := 1; end for; end if; // Boundary z-axis if node_points[nod,6] <> 0 then for i in 1:Nodes*3 loop G_total[nod*3,i] := 0; G_total[nod*3,nod*3] := 1; end for; end if; end for; // Solving displacement (U) // [G_total]{U} = {F} displacement := Modelica.Math.Matrices.solve(G_total, Forces); // Solving reaction (R) // {R} = [G_total]{U} - {F} reaction := (G_total_new * displacement) - Forces; // Eliminating float errors for i in 1:nodes_xyz loop reaction[i] := if abs(reaction[i])<= err_solve then 0 else reaction[i]; displacement[i] := if abs(displacement[i])<= err_solve then 0 else displacement[i]; end for; // Checking Sum of Forces // For static structures, should add up to 0 check[1]:=sum({reaction[i] for i in (1:3:(nodes_xyz-2))}) + sum({Forces[i] for i in (1:3:(nodes_xyz-2))}); check[2]:=sum({reaction[i] for i in (2:3:(nodes_xyz-1))}) + sum({Forces[i] for i in (2:3:(nodes_xyz-1))}); check[3]:=sum({reaction[i] for i in (3:3:(nodes_xyz))}) + sum({Forces[i] for i in (3:3:(nodes_xyz))}); // Eliminating float errors for i in 1:3 loop check[i] := if abs(check[i])<= err_check then 0 else check[i]; end for; // === POSTPROCESSING PHASE === //Calculating stress in each truss for mem in 1:Members loop // Getting node coordinate data and displacement for current member for i in 1:3 loop node_i[i] := node_points[node_pairs[mem,1],i]; node_j[i] := node_points[node_pairs[mem,2],i]; mem_dis[i] := abs(displacement[3*(node_pairs[mem,1]-1)+i] - displacement[3*(node_pairs[mem,2]-1)+i]); end for; // Solving the Displacement Matrix for current member // Stress = F/A = (E/L) * [dL] L:=Modelica.Math.Vectors.length(node_j-node_i); cos_x:=(node_j[1]- node_i[1])/L; cos_y:=(node_j[2]- node_i[2])/L; cos_z:=(node_j[3]- node_i[3])/L; T := (Elas/L)*[cos_x^2,cos_x*cos_y,cos_x*cos_z; cos_y*cos_x,cos_y^2,cos_y*cos_z; cos_z*cos_x,cos_z*cos_y,cos_z^2]; stress_vector:=(T*mem_dis); // The stress in vector components {stress_x, stress_y, stress_z} stress_mag[mem]:=Modelica.Math.Vectors.length(stress_vector); // Magnitude/length of the stress vectors for each member end for; // Safety Factor // Factor of Safety = (Yield Strength of material) / (Stress experienced by Member) for mem in 1:Members loop if stress_mag[mem]>0 then safety[mem] := Yield / stress_mag[mem]; else safety[mem] := 0; end if; end for; end Truss_3D_Tugas_Besar_Safety;

Cost and Cost Efficiency of Area Variation and Constrained Elasticity

Berdasarkan hasil analisis dari struktur truss sebelumnya, pada bagian ini parameter area cross-section dari setiap truss dalam struktur tesebut divariasi dengan menggunakan yann standar ukuran yang tertera pada prosedur. Data perbedaan Area dalam material yang sama tersebut diproses menjadi variasi harga (cost) dan efisiensi harga (safaty factor/cost) agar dapat diproses di akhir dengan optimisasi metode Powell. Dalam Prosesnya, data harga material yang terbatas di pasar akan di-regresi orde 2 untuk mendapatkan nilai harga pada ukuran-ukuran yang tidak memilik harga. Berikut adalah tabel hasil proses data dan grafik yang dihasilkan:

• 

  Table 3 Table of cost and cost efficiency

• 

  Figure 3 Graph of cost and cost efficiency

Cost and Cost Efficiency of Elasticity Variation and Constrained Area

Mirip dengan pengolahan data sebelumnya, bagian ini juga akan menghasilkan nilai variasi harga (cost) dan efisiensi harga (safaty factor/cost) untuk diolah dengan optimisasi metode powell di akhir. Namun, berbeda dengan sebalumnya, pada bagian ini area cross section akan disamakan nilainya dan elastisitas material akan divariasikan. Dalam pengolah ini, metode curve fitting orde 2 juga dilakukan untuk mendapatkan nilai-nilai harga, kekuatan yield, densitas yang kosong. Berikut adalah tabel hasil proses data dan grafik yang dihasilkan:

• 

  Table 4 Table of cost and cost efficiency

• 

  Figure 4 Graph of cost and cost efficiency

Result and Analysis

Hasil dari studi ini adalah nilai cost efficincy optimal yang dipengaruhi oleh 2 faktor yaitu elastisitas dan area. Oleh karena itu jenis optimisasi yang dilakukakn adalah optimasi multivariabel (lebih dari satu variabel) dengan metode powell. Metode ini dipilih karena tergolong cukup sederhana dan cepat mencapai titik konvergen. Secara singkat, agar lebih memahami optimasi multivariabel ini, bisa dibayangkan distribusi cost efficieny digambarkan sebagi grafik 3 dimensi dalam sumbu x, y, z, Cost efficiency akan diwakili skala z, elastisitas sebagai skala y dan area sebagai skala x. Nilai dari cost efficincy(z) pada suatu titik adalah fungsi area(x) dan elastisitas (y). Persamaan kuadrat cost efficiency yang diperolah dari curve fitting tahap sebelumnya pada 2 variasi, yaitu area dan elastisitas, akan di rata-ratakan sehingga menjadi fungsi tunggal. Hasil tersebut akan menghasikan permukaan kontur bidang xy yang memiliki titik optimum. koordinat titik tersebut akan dicari secara linear satu persatu dengan metode Golden. Hasil dari simulasi dan program tertulis Modelica akan ditampilkan sebagai berikut:

Modelica Result

Powell Method model Powell_Method_Opt_Jos_Ratio parameter Integer order=2; //order of regression parameter Integer GoldN=15; //maximum iteration of linear gold optimization parameter Integer N=10; //maximum iteration parameter Real maxerror=1e-50; //maximum error //Assumed to be X function parameter Real RatioArea[size(Area,1)]={3.6534e-5, 3.75348e-5, 3.79639e-5, 3.79507e-05, 3.74821e-05}; parameter Real CostArea[size(Area,1)]={829776,1053955,1308215,1722235,2202651}; parameter Real Area[:]={141e-6, 184e-6, 231e-6, 304e-6, 384e-6}; Real CoeArea[order+1]; Real CoeAreaCost[order+1]; //Assumed to be Y function parameter Real RatioElas[size(Elas,1)]={2.83875e-5, 8.5637e-5, 12.5153e-5, 5.70447e-5}; parameter Real CostElas[size(Elas,1)]={1367994,732878,492691,622939}; parameter Real Elas[:]={192e9, 195e9, 198e9, 201e+9}; Real CoeElas[order+1]; Real CoeElasCost[order+1]; //Guessing Start (1 point, 2 vector) - the second vector will be auto-generated parameter Real StartPoint[2]={171e-6,195e9}; parameter Real VectorPoint[2]={205e-6,200e9}; //Solution Real RatioZ[N]; //Optimum result Real CostZ; //Cost Optimum result Real x3[2]; //Optimum point coordinate (x,y) Real error[N]; //Error protected Real jos, RatioZX, RatioZY, CostZX, CostZY; Real x0tes[2], x0[2], x1[2], x2[2]; Real grad1, grad2, grad3; //Checking approach Real cek0[N], cek1[N], cek2[N], cek3[N], yek0[N], yek1[N], yek2[N], yek3[N]; algorithm //Creating function at both CoeArea:=Curve_Fitting(Area,RatioArea,order); CoeElas:=Curve_Fitting(Elas,RatioElas,order); CoeAreaCost:=Curve_Fitting(Area,CostArea,order); CoeElasCost:=Curve_Fitting(Elas,CostElas,order); //Creating 2 vector defined in 2 gradient (grad1,grad2) x0:=StartPoint; x0tes:=VectorPoint; grad1:=(x0tes[2]-x0[2])/(x0tes[1]-x0[1]); grad2:=-0.5*grad1; //Provide error variable RatioZX:=CoeArea[1]*x0[1]^2+CoeArea[2]*x0[1]+CoeArea[3]; RatioZY:=CoeElas[1]*x0[2]^2+CoeElas[2]*x0[2]+CoeElas[3]; jos:=(RatioZX+RatioZY)/2; //Computing approach for i in 1:N loop cek0[i]:=x0[1]; yek0[i]:=x0[2]; x1:=Gold_Opt_Func_2D_Vector(x0,grad1,GoldN,CoeArea,CoeElas); cek1[i]:=x1[1]; yek1[i]:=x1[2]; x2:=Gold_Opt_Func_2D_Vector(x1,grad2,GoldN,CoeArea,CoeElas); cek2[i]:=x2[1]; yek2[i]:=x2[2]; if (x0[1]==x2[1]) and (x0[2]<>x2[2]) then grad3:=1/0; //divergent gradient, please subtitute vector x3:=Gold_Opt_Func_2D_Vector_1(x2,1,GoldN,CoeArea,CoeElas); else grad3:=(x2[2]-x0[2])/(x2[1]-x0[1]); x3:=Gold_Opt_Func_2D_Vector(x0,grad3,GoldN,CoeArea,CoeElas); end if; cek3[i]:=x3[1]; yek3[i]:=x3[2]; grad1:=grad2; grad2:=grad3; //Function Created RatioZX:=CoeArea[1]*x3[1]^2+CoeArea[2]*x3[1]+CoeArea[3]; RatioZY:=CoeElas[1]*x3[2]^2+CoeElas[2]*x3[2]+CoeElas[3]; RatioZ[i]:=(RatioZX+RatioZY)/2; CostZX:=CoeAreaCost[1]*x3[1]^2+CoeAreaCost[2]*x3[1]+CoeAreaCost[3]; CostZY:=CoeElasCost[1]*x3[2]^2+CoeElasCost[2]*x3[2]+CoeElasCost[3]; CostZ:=(CostZX+CostZY)/2; x0:=x3; error[i]:=abs(1-jos/RatioZ[i]); jos:=RatioZ[i]; if error[i]<maxerror then break; end if; end for; end Powell_Method_Opt_Jos_Ratio;

Gold Optimization Function function Gold_Opt_Func_2D_Vector //only for second order polynomial input Real Point[2]; input Real gradient; //Gradient input Integer N; //Maximum iteration input Real coe1[3]; input Real coe2[3]; output Real optimum[2]; protected Real C=Point[2]-Point[1]*gradient; Real xlo=141e-6; Real xhi=384e-6; Real f1[N], f2[N], x1[N], x2[N], ea[N]; Real fx; Real es=1e-10; // maximum error Real d, xl, xu, xint, R=(5^(1/2)-1)/2; algorithm xl := xlo; xu := xhi; for i in 1:N loop d:= R*(xu-xl); x1[i]:=xl+d; x2[i]:=xu-d; f1[i]:=((coe1[1]*x1[i]^2+coe1[2]*x1[i]+coe1[3])+(coe2[1]*(gradient*x1[i]+C)^2+coe2[2]*(gradient*x1[i]+C)+coe2[3]))/2; f2[i]:=((coe1[1]*x2[i]^2+coe1[2]*x2[i]+coe1[3])+(coe2[1]*(gradient*x2[i]+C)^2+coe2[2]*(gradient*x2[i]+C)+coe2[3]))/2; xint:=xu-xl; if f1[i]>f2[i] then xl:=x2[i]; optimum[1]:=x1[i]; optimum[2]:=(gradient*x1[i]+C); fx:=f1[i]; else xu:=x1[i]; optimum[1]:=x2[i]; optimum[2]:=(gradient*x2[i]+C); fx:=f2[i]; end if; ea[i]:=(1-R)*abs((xint)/optimum[1]); if ea[i]<es then break; end if; end for; end Gold_Opt_Func_2D_Vector; Curve Fitting function Curve_Fitting input Real X[:]; input Real Y[size(X,1)]; input Integer order=2; output Real Coe[order+1]; protected Real Z[size(X,1),order+1]; Real ZTr[order+1,size(X,1)]; Real A[order+1,order+1]; Real B[order+1]; algorithm for i in 1:size(X,1) loop for j in 1:(order+1) loop Z[i,j]:=X[i]^(order+1-j); end for; end for; ZTr:=transpose(Z); A:=ZTr*Z; B:=ZTr*Y; Coe:=Modelica.Math.Matrices.solve(A,B); end Curve_Fitting;