Komputasi Fisika 1 Simulasi dan Visualisasi Fisika dengan C, C++, Gnuplot, dan OpenGL

Komputasi Fisika 1 Simulasi dan Visualisasi Fisika dengan C, C++, Gnuplot, dan OpenGL

Editor:

Sparisoma Viridi Suprijadi

Departemen Fisika, FMIPA, Institut Teknologi Bandung 2015

Komputasi Fisika 1 Simulasi dan Visualisasi Fisika dengan C, C++, Gnuplot, dan OpenGL

Editor:

Sparisoma Viridi Suprijadi

Departemen Fisika, FMIPA, Institut Teknologi Bandung 2015

Buku ini disusun dengan bantuan dukungan Riset ITB tahun 2015.

Sparisoma Viridi dan Suprijadi (editor), Komputasi Fisika 1: Simulasi dan Visualisasi Fisika dengan C, C++, Gnuplot, dan OpenGL, Departemen Fisika, Institut Teknologi Bandung Lisensi ITB 2015 Penyebaran dalam bentuk aslinya diperbolehkan selama untuk pendidikan dan tak-komersial. Modifikasi dan pembuatan karya turunan memerlukan ijin dari penulis [email protected] | [email protected]

Pengantar Buku Komputasi Fisika 1: Simulasi dan Visualisasi Fisika dengan C, C++, Gnuplot, dan OpenGL ini merupakan kumpulan tugas RBL (research based learning) pengganti U2 (ujian kedua) pada matakuliah FI5005 Komputasi Sistem Fisis yang merupakan layanan dalam Program Magister Fisika, Fakultas Matematika dan Ilmu Pengetahuan Alam, Institut Teknologi Bandung, Indonesia. RBL merupakan suatu metode pembelajaran yang sudah umum diberikan pada program studi-program studi dalam eks-Departemen Fisika, Fakultas Matema-tika dan Ilmu Pengetahuan Alam, Institut Teknologi Bandung, sejak sekitar tahun 2007-an. Saat itu matakuliah-matakuliah Fisika Dasar (IA, IB, IIA, dan IIB) meru-pakan menjadi implementasi pertamanya. Antusiasme dari pengajar dan peserta ajar membuat metode ini kemudian menjadi rutin dan populer digunakan dalam semua jenjang (S1, S2, dan S3). Khusus dalam kuliah terkait komputasi metode pengajaran ini dianggap baik karena peserta ajar dapat menerapkan hasil pemahaman mereka dalam suatu permasalahan nyata yang telah dipilih dan kemudian menyelesaikannya secara numerik. Setiap tugas RBL tersebut disajikan dalam tiap bab, telah dilengkapi dengan deskripsi permasalahan fisisnya, penyelesaiannya sampai pada bentuk (baik secara sengaja ataupun tidak) yang tidak diselesaikan secara analitik melainkan secara numerik, metode numerik yang digunakan, algoritma program, hasil yang diperoleh dan diskusinya, referensi, dan lampiran kode-kodenya (minimal berkas .cpp dan .gps). Pembaca diharapkan dapat memanfaatkan contoh-contoh RBL tersebut dalam mempelajari sistem fisis yang dipilih, metode numerik, dan penerapannya dalam C++ dan visualisasinya menggunakan Gnuplot (dan OpenGL). Salah satu alasan mengapa C++ (lebih tepatnya G++) dan Gnuplot yang dipilih adalah karena keduanya telah tersedia dalam sistem operasi Linux sehingga diharapkan dapat menurunkan hambatan awal bagi seorang mahasiswa fisika untuk belajar simulasi dan visualisasi dengannya, yang tidak perlu terlebih dahulu melakukan instalasinya (yang kadang tidak mudah karena banyaknya kebergantungan pada aplikasi dan pustaka lain). Aplikasi pemrograman yang lebih ramah pengguna (user friendly) seperti Matlab dan Scilab sedapat mungkin dihindari karena terlalu banyak menyembunyikan detil proses komputasi, yang dirasakan dapat memanjakan peserta ajar sehingga tidak terbiasa berpikir kritis dan teknis. Bandung, Desember 2015 Sparisoma Viridi & Suprijadi (editor)

i

Komputasi Fisika 1: Pengantar

ii

Cara merujuk Setiap bab dalam buku Komputasi Fisika 1: Simulasi dan Visualisasi Fisika dengan C, C++, Gnuplot, dan OpenGL ini dapat dirujuk dengan cara Nama_Penulis_1, Nama_Penulis_2, .., "Judul bab" dalam Komputasi Fisika 1: Simulasi dan Visualisasi Fisika dengan C, C++, Gnuplot, dan OpenGL oleh S. Viridi, Suprijadi (editor), 2015, pp. halaman_mulai-halaman_akhir. Atau dengan menggunakan rujukan untuk format referensi berbentuk chapter of the book, misalnya untuk bab 1 N. Faizin, A. M. Ilmah, Hardiyanto, " Pemodelan Aliran Fluida Pipa dengan Penghalang " dalam Komputasi Fisika 1: Simulasi dan Visualisasi Fisika dengan C, C++, Gnuplot, dan OpenGL oleh S. Viridi, Suprijadi (editor), 2015, pp. 1-10.

iii

Komputasi Fisika 1: Cara merujuk

iv

Daftar isi i

Pengantar

iii

Cara merujuk Daftar isi

v

1

Pemodelan Aliran Fluida Pipa dengan Penghalang N. Faizin, A. M. Ilmah, Hardiyanto

1

2

Simulasi Bouncing Ball H. D. Ibrahim, A. Yanitama, H. D. Bhakti

11

3

Model Lintasan Shuttlecock dengan Variasi Gaya Hambat Udara A. Sutiono, R. Syam

25

4

Simulasi Gerakan Air Dangkal 1 Dimensi M. M. Fahmi

35

5

Analisis Titik Luluh Material Dengan Menggunakan Metode Secant F. Nakul, M. S. Lestari, L Handayani, Deni

45

6

Gelombang Tali Fatimah, N. Novita, A. Wiliardy

59

7

Pemodelan Lintasan Benda Titik pada Tong Setan Rustan, E. L. Yani N. S. T., W. Wahyuni, M. Husnah

103

8

Osilasi Teredam Sistem Pegas dengan Massa Beban tak Konstan Kamirul, L. Y. Roodhiyah, H. Syahwanti, Fatriani

113

9

Orientasi Gerak Jatuh Bebas Sistem 3 Partikel J. Pebralia, Y. C. Dewi, D. Dwiputra, Rouf

129

10 Inverted Pendulum Pada Sebuah Kereta H. Hamzah, Firmansyah, M. F. Alfian R. S., M. B. Mutofa

v

141

Komputasi Fisika 1: Daftar isi

vi

|1

Pemodelan Aliran Fluida Pipa dengan Penghalang Nur Faizin | [email protected] Aurista M. Ilmah | [email protected] Hardiyanto | [email protected] Aliran fluida dalam pipa dibagi menjadi dua yaitu aliran laminer dan turbulen. Pemodelan dilakukan untuk fluida laminer. Fluida yang dimodelkan merupakan fluida Bernoulli yaitu = 0, v = konstan dan inkompresibel. Fungsi aliran fluida dan distribusi laju aliran diselesaikan dengan menggunakan metode Finite Different. Dengan menggunakan metode ini akan dapat divisualisasikan dengan komputasi bentuk aliran fluida yang diberi gangguan berupa bola dengan jari-jari tertentu. Dalam aliran suatu fluida bergantung pada tekanan dan kecepatan. Distribusi kecepatan dan tekanan untuk sebuah medan aliran dapat diketahui dari pola garis arus dan dari penerapan persamaan bernoulli. Pola-pola garis arus dapat membentuk suatu jaring-jaring aliran yang berkesusaian, begitu banyak persamaan yang digunakan dalam menjelaskan persamaan aliran fluida salah satunya adalah persamaan dalam menghitung tekanan dan kecepatan, persamaan yang digunakan adalah persamaan bernauli dan persamaan kontinuitas. Pemodelan secara matematis dapat dikelompokkan atas dua bagian yaitu model fisis dan model matematis. Dimana model matematis adalah model sistem persamaan fisika-matematik untuk menjelaskan sifat fisis dari proses yang diteliti dan dianalisa dari persamaan mekanika pergerakan fluida yang dirangkum dalam persamaan kontinuitas dan persamaan bernoulli.

1

2

| 1. PEMODELAN ALIRAN FLUIDA PIPA DENGAN PENGHALANG

1.1

Teori

Di sekitar kita banyak sekali fenomena aliran fluida. Misalnya, aliran pada keran air. Ketika keran dibuka sedikit, air yang keluar dari keran tampak bening, aliran seperti ini disebut aliran laminar. Namun, ketika keran dibuka lebih lebar lagi, fluida akan mengalir seperti semburan berwarna putih. Aliran inilah yang disebut dengan aliran turbulen Contoh lain yang dapat kita ambil adalah aliran sungai dibawah jembatan. Ketika kita berdiri di jembatan sungai, dapat terlihat aliran sungai dibawah jembatan. Misalnya aliran sungai yang melewati kaki jembatan. Aliran yang menabrak kaki jembatan akan terlihat seperti gambar 1.

Gambar 1.1: Aliran dengan Penghalang (bimoauliac.blogspot)

Aliran sungai yang menabrak tiang jembatan akan tampak terbelah menjadi dua aliran yang sejajar dengan permukaan tiang. Namun, pada bagian belakang tiang dimana aliran sungai meninggalkan tiang makan akan terlihat putaran air. Selain itu, aliran sungai pada bagian tengah dan pinggir sungai dapat juga kita amati. Pada bagian tengah aliran sungai, jika ada objek yang terseret aliran (daun, sampah, dll), maka objek tersebut akan bergerak dengan cepat dan tanpa berputar, tetapi objek yang terseret aliran di pinggiran sungai akan bergerak lebih lambat dan berputar. Agar lebih jelas, mari kita perhatikan gambar2. Kemampuan mengalirnya fluida disebabkan adanya perpindahan partikel fluida karena pengaruh tegangan geser, tegangan geser dalam suatu fluida sebanding dengan laju perubahan kecepatan ruang yang normal terhadap aliran. Aliran fluida yang incompresible dalam pipa dalam mempelajari aliran fluida seringkali digunakan asumsi fluida ideal, fluida ideal merupakan fluida yang tidak kental, pada fluida nyata timbul tegangan geser antara partikel - partikel fluida ketika partikel tersebut bergerak pada kecepatan yang berbeda pada fluida ideal yang mengalir melalui suatu tabung lurus, semua partikel bergerak pada garis-garis sejajar dengan kecepatan sama. Pada aliran fluida nyata, kecepatan terdekat dengan dinding akan nol, dan akan bertambah besar pada jarat pendek dari dinding (Orianto, 1989).

3

1.2. ALIRAN FLUIDA INCOMPRESIBLE

Gambar 1.2: Aliran sungai dengan Penghalang (bimoauliac.blogspot)

1.2

Aliran Fluida Incompresible

Garis-garis arus adalah kurva-kurva khayal yang ditarik melalui suatu fluida untuk menunjukkan arah gerakan di berbagai bagian aliran dari sistem fluida garis, untuk aliran steady dan incompresible dengan menggunakan persamaan bernouli.

V2 P + + gZ = Konstan 2 ρ

1.3

(1.1)

Pola Aliran dengan dua Dimensi

Garis-garis arus membagi aliran menjadi beberapa saluran dengan laju aliran yang sama, jika total aliran persatuan kedalaman normal ke bidang garis arus adalah Q, serta jumlah dari saluran adalah n, maka aliran yang melalui panjang tiap-tiap saluran adalah q = Q/n dengan kecepatan rata-rata V ,dalam saluran dengan lebarnb adalah V = q/b dan kecepatan pada suatu titik berbanding terbalik dengan jarak garis arus pada titik. Pada Gambar 3 dimisalkan aliran dua dimensi pada suatu pipa dinyatakan dengan empat garis arus , dua garis arus serupa dengan batas padat, dan aliran total dibagi menjadi tiga saluran. Bagian dari pipa yang memiliki kecepatan V0 dan lebar saluran b0 , sehingga titik P dapat dinyatakan dengan persamaan kontinuitas sebagai berikut:

V 1 b 1 = V 0 b0

(1.2)

4


Gambar 1.3: Pola aliran dua dimensi pada suatu pipa (Ngana,1996)

1.4

Solusi Numerik dan Algoritma

Metode beda hingga merupakan metode yang digunakan untuk memecahkan persamaan diferensial parsial secara numerik, dengan menggunakan deret taylor yang diputuskan pada orde tertentu yang sesuai kebutuhan yang ada. Pada metode beda hingga memperkirakan turunan parsial dalam persamaanpersamaan fisis dengan perbedaan-perbedaan diantara nilai node pada titik suatu jarak terpisah terhingga. Sehingga nilai node digunakan sebagai nilai pengganti persamaan differensial ke persamaan aljabar. Pemakaian metode beda hingga ini dapat digunakan sebagai penghitung harga-harga fungsi arus di titik titik perpotongan kisi siku-siku yang dapat menutup seluruh domain. Garis-garis arus digambarkan dengan menginterpolasikan harga-harga fungsi arus yang ditentukan pada titik-titik tersebut. Pada persamaan laplace untuk aliran fluida sangat cocok untuk teknik-teknik solusi numerik dengan finite difference (beda hingga). Variabel tidak bebas yang sesuai adalah fungsi aliran ψ yang memenuhi persamaan laplace.

∂2ψ ∂2ψ + =0 ∂x2 ∂y 2

(1.3)

Teknik beda hingga membagi medanaliran kedalam node yang sama seperti pada gambar 4. Subskripsi i dan j menunjukkan posisi dari sembarang node ruang yang sama

1.4. SOLUSI NUMERIK DAN ALGORITMA

5

Gambar 1.4: Sketsa Beda Hingga (Nganah,1996) dan ψ( i, j) Merupakan nilai fungsi aliran pada titik node dengan didapat bentuk finnite difference dari persamaan dari pendekatan node adalah :

ψi,j = (ψ(i+1,j) + ψ(i−1,j) + ψ(i,j+1) + ψ(i,j−1) )/4

(1.4)

Dan dari persamaan tersebut dapat diaplikasikan kedalam notasi pemrograman secara langsung kedalam bahasa komputer, jika A(i, j) merupakan variabel fungsi aliran. Pernyataan bahasa pemrograman komputer dapat dituliskan denga persamaan aljabar yaitu:

A(i, j) = (A(i + 1, j) + A(i − 1, j) + A(i, j + 1) + A(, j − 1))/4

(1.5)

Aliran Fuida dengan penghalang. Aliran fluida pada dasarnya menempati ruang yang dia lalui, namun bila ada penghalang disekitar fluida tersebut, hal itu membuat fluida tidak akan mengalir pada tempat yang dilalui penghalang, pada pemodelan ini diberikan penghalang yang berupa elips dan lingkaran.

6


Syarat Batas Dalam pemodelan hanya meninjau kecepatan pada arah x dan y. v = vx + vy

(1.6)

Kecepatan yang ditinjau konstan terhadap waktu dv =0 dt

(1.7)

Algoritma L01. L02. L03. L04. L05. L06. L07. L08.

1.5

vx(awal) , , vy(awal) ? Nx , Ny ? Sawal = 0 Sold = S(i, j) Sgangguan (i, j) Ser = Sgangguan (i, j) − S(i, j) < 10−10 Sold = Sgangguan (i, j) vx(i,j) = S(i + 1, j) − S(i, j), vy(i,j) = S(i, j + 1) − S(i, j)

Hasil dan diskusi

Telah dilakukan pemodelan aliran fluida Bernoulli. Aliran fluida tersebut diberi gangguan berupa lingkaran dan elips. Nilai-nilai parameter yang digunakan dalam pemodelan sebagaimana ditunjukkan pada Tabel 1.1

Tabel 1.1: Nilai-nilai parameter simulasi. Parameter vx(awal) vy(awal) nodax noday Rlingkaran Pelips dan Qelips

Nilai 50 0 147 67 17 17dan9

Pemodelan dilakukan berulang dengan iterasi ke 6017, hasil pemodelan yang diperoleh untuk gangguan berupa lingkaran sebagaimana ditunjukkan pada Gambar5.

1.5. HASIL DAN DISKUSI

7

Gambar 1.5: Kontur aliran fluida dengan gangguan lingkaran

Gambar 5.menunjukkan bahwa laju aliran fluida dalam pipa hampir sesuai dengan teori. Tetapi nilai laju aliran di pojok kiri atas dan kanan atas diperoleh nilai yang sangat besar yaitu sekitar 4, 5x105 - 5x105 . Hal ini kemungkinan disebabkan oleh perhitungan nilai matriks yang diulang-ulang. Profil laju aliran disekitar gangguan terlihat ada anomali, hal ini menunjukkan adanya perubahan laju aliran saat mengenai gangguan. Jika gangguan yang berupa lingkaran dipindahkan pada posisi kiri bawah pipa diperoleh gambar berikut.

Gambar 1.6: Gangguan lingkaran di kiri bawah Sedangkan untuk gangguan berupa elips diperoleh profil laju aliran sebagai berikut. Bentuk gangguan yang terbentuk pada Gambar 7 tidak membentuk gangguan berupa elips. Hal ini berhubungan dengan jumlah matriks yang digunakan dan ketelitian perhitungan pada program.

8


Gambar 1.7: Kontur aliran fluida dengan gangguan berupa elips

1.6

Referensi

Munson-Bruce R, Young-Donald F and Okiishi-Theodore H, 1994.Fundamentals of Fluid Mechanics, 4thedition, 1.6 viscosity, p.18 Ngana , Frederika Rambu (1996) Perhitungan distribusi kecepatan aliran fluida pada aliran potensial dengan menggunakan metoda beda hingga dua dimensi. Undergraduate thesis, FMIPA Undip Orianto, M., Pratikto, W. A.. 1989. ”Mekanika Fluida 1”. yogyakarta: BpFE bimoauliac.blogspot

1.7

Lampiran

Code C #include <stdio.h> #include <math.h> int main(){ FILE *file1, *file2; int maxIter = 10000; int Ny = 67; int Nx = 147; int i,j;\ float vx = 50; float vy = 0.;

1.7. LAMPIRAN

float hawal = 0; float float float float

S[Ny][Nx]; Vx[Ny][Nx],Vy[Ny][Nx]; V[Ny][Nx]; S_old;

// inisialisasi medium bernilai 0 diseluruh posisi for( i = 0; i < Ny ; i++){ for( j = 0 ; j < Nx ; j++){ S[i][j] = 0; } } // Profil kecepatan awal for( i = 0; i < Ny ; i++){ S[i][Nx-1] = 20 * vx* i; //di dinding kanan S[i][0] = 20 * vx* i; //di dinding kiri } for( j = 0; j < Nx ; j++){ S[0][j] = 20 * vy* j; //di dinding bawah S[Ny-1][j] = 20 * vy *j; //di dinding atas } //profil poiseuille flow S_old = hawal; int k = 0; float S_sum; while(k <= maxIter){ for(i = 1; i < Ny-1 ; i++){ for(j = 1; j < Nx-1; j++){ if(((j-130)*(j-130) + (i-0)*(i-0)) <= 289) S[i][j] = S[34][74]; // yg ini mksdnya apa? else S[i][j] = 0.25*(S[i-1][j]+S[i+1][j] + S[i][j-1]+S[i][j+1]); } } S_sum = 0; for(i = 0; i < Ny ; i++){ for(j = 0; j < Nx; j++){ S_sum += S[i][j]; } } S_sum = S_sum/(Nx*Ny); printf("iterasi ke-%d -> %.8f \r\b ",k,S_sum - S_old); fflush(stdout); if (fabs(S_sum - S_old) < 1e-10) break;

9

10


else S_old = S_sum; k++; } printf("\n"); file1 = fopen("dataKecepatan.txt","w"); file2 = fopen("dataVektorKecepatan.txt","w"); for( i = 1; i < Ny-1 ; i++){ for( j = 1; j < Nx-1 ; j++){ Vx[i][j] = S[i+1][j] - S[i][j]; Vy[i][j] = S[i][j+1] - S[i][j]; V[i][j] = sqrt(Vx[i][j]*Vx[i][j] + Vy[i][j]*Vy[i][j]); fprintf(file1,"%.3f\t",V[i][j]); fprintf(file2,"%d\t%d\t%.3f\t%.3f\t\n",i,j,Vx[i][j],Vy[i][j]); } fprintf(file1,"\n"); } fclose(file1); fclose(file2); return 0; }

Kode C di atas dikompilasi dengan cara:

g++ code.cpp -o fluidajadi dan dipanggil dengan cara:

./fluidajadi

Cara plot di Gnuplot Cara mengeplot gambar dalam Gnuplot dapat menulis perintah-perintah berikut :

> set xlabel ’arah x’; > set ylabel ’arah y’; > plot "dataKecepatan.txt" matrix with image;

|2

Simulasi Bouncing Ball Haerul Djusmar Ibrahim | [email protected] Arka Yanitama | [email protected] Henny Dwi Bhakti | [email protected] Simulasi fenomena bouncing ball dilakukan dengan validasi eksperimen dan metode Euler sebagai solusi numerik. Eksperimen dilakukan dengan bola ping pong yang dijatuhkan secara vertikal pada ketinggian 0.2m. Dalam simulasi ini dilakukan variasi perbedaan ketinggian, massa, dan jari-jari bola terhadap waktu berhenti bola.

2.1

Teori

Ketika sebuah benda dijatuhkan dari ketinggian tertentu, maka percepatan yang dialami akan sama dengan percepatan gravitasi yaitu ~a = −~g . Berdasarkan Hukum kedua Newton, maka persamaan geraknya dapat ditulis sebagai berikut

X

F~ = m~a,

(2.1)

dengan m adalah massa benda dan ~a adalah percepatan. Pada kasus ini, jika kondisi sekitar dan dimensi benda tersebut diperhitungkan, maka bekerja drag force yang arahnya berlawanan dengan arah geraknya. 1 F~d = ρv~2 Cd A 2

1

(2.2)

2

| 2. SIMULASI BOUNCING BALL

maka persamaan geraknya menjadi m~a = F~g − F~d 1 m~a = m~g − ρ~v Cd A|~v | 2 ~a = ~g −

1 ρ~v Cd A|~v | 2 m

(2.3)

dimana ρ adalah densitas udara, v adalah kecepatan gerak benda, Cd adalah koefisien drag, A adalah luas permukaan benda, dan a merupakan percepatan gerak benda tiap saat.

Gambar 2.1: skema drag force Pada sistem Bouncing Ball ini terdapat keadaan transisi pada saat bola menumbuk permukaan lantai(y = 0) 3. Karena gaya drag berlawanan arah dengan arah gerak maka pada persamaan 2.2 kecepatannya bernilai negatif. Pada saat kondisi bola mengarah ke atas, gaya drag dan gaya gravitasi sama-sama mengarah ke bawah. Kemudian ketika kondisi bola mengarah kebawah, gaya drag mengarah ke atas sedangkan gaya gravitasi tetap mengarah ke bawah seperti pada gambar 2.1. Maka, pada persamaan umum gaya drag (F = 0.5CdAv 2 ), kita dapat memberi nilai absolut pada salah satu kecepatan dan mengalikannya dengan nilai −1 sehingga persamaannya menjadi (−0.5CdAv|v|)[2]

2.2

Syarat batas

Kecepatan dibatasi hanya memiliki komponen pada arah y karena bola bergerak turun naik pada arah y. Selain itu diberi kondisi tumbukan terhadap permukaan lantai seperti pada gambar 2.2 dan 2.3. Pada kondisi ketika tumbukan δy = y − 0, kondisi tersebut menunjukkan bahwa y bernilai negatif, maka δy = −y. Jika kondisi ini terjadi dianggap bola telah menumbuk permukaan lantai. Kemudian jika kondisi ini terpenuhi maka y = −y

3


Gambar 2.2: tumbukan

2.3

Kondisi ketika

Gambar 2.3: tumbukan

Kondisi setelah

Solusi numerik dan algoritma

Pada simulasi ini digunakan metode yang cukup sederhana yaitu metode Euler. Persamaan (2.3) dapat memberikan kecepatan dan posisi setiap saat untuk arah x, dimana berlaku

1 ρ~v Cd A|~v | 2 m vx (t + ∆t) = vx (t) + ax (t)∆t

(2.5)

rx (t + ∆t) = rx (t) + vx (t)∆t

(2.6)

ax (t) = −

(2.4)

dan untuk arah y berlaku

1 ρ~v Cd A|~v | 2 m vy (t + ∆t) = vy (t) + ay (t)∆t ry (t + ∆t) = ry (t) + vy (t)∆t ay (t) = −mg −

(2.7) (2.8) (2.9)

Dari persamaan-persamaan yang telah didapat, kemudian dibuat algoritma untuk program sebagai berikut: L01. L02. L03. L04. L05. L06. L07. L08. L09.

∆t, t0 = 0, Tmax = 10? G, Cd , e, ρ, R, A, m, x0 , y0 , vx , vy t = t0 ay (t) = G − 0.5ρvy (t)Cd Am|vy (t)|; ax (t) = −0.5ρvx (t)Cd Am|vx (t)| vy (t + ∆t) = vy (t) + ay (t)∆t; vx (t + ∆t) = vx (t) + ax (t)∆t collision? → L07 vy (t + ∆t) = −vy (t) + evy y(t + ∆t) = y(t) + vy (t)∆t Ek (t + ∆t) = 0.5m(vx (t)vx (t) + vy (t)vy (t))

4


L10. StopCriteria = |Ek (t + ∆t) − Ek (t)| L11. t = t + ∆t L12. t ≤ tmax or StopCriteria < 10−8 → L04

2.4

Hasil dan diskusi

Eksperimen yang dilakukan menghasilkan video Bouncing Ball dengan menggunakan bola ping-pong dengan massa 0.0022kg dan jari-jari 0.02m yang dijatuhkan dari ketinggian 0.2m. Dari data posisi dan waktu yang diperoleh dari video tersebut kemudian digunakan untuk validasi hasil simulasi. seperti pada gambar 0.2

Simulation Experiment

Position (m)

0.15

0.1

0.05

0 0

1

2

3

4 Time (s)

5

6

7

8

Gambar 2.4: grafik posisi terhadap waktu untuk validasi hasil eksperimen dan simulasi Hasil validasi tersebut digunakan untuk melakukan simulasi dengan parameter yang tercantum dalam Tabel 2.1. Tabel 2.1: Nilai-nilai parameter simulasi. Parameter ∆t t0 tmax h StopCriteria ··

Nilai 0.001 0 10 0.2 10−8 ··

Parameter m R g Cd ρ e

Nilai 0.0022 0.02 9.8 0.47 1.29 0.05

Dari hasil validasi pada gambar 2.4, menunjukkan bahwa hasil simulasi valid sampai t≈1 s. Pada saat eksperimen, ada beberapa faktor yang menyebabkan

5


hasil yang kurang akurat pada waktu berikutnya, hal tersebut disebabkan karena ketika bola dijatuhkan bola memiliki komponen kecepatan pada arah x. Selain itu, bola yang dijatuhkan mengalami rotasi sehingga ada energi yang berkurang. Pada metode Euler, terdapat kekurangan yaitu pada saat menggunakan ∆t yang lebih besar yang ditunjukkan pada Gambar 2.5. Hal menunjukkan bahwa pada ∆t = 0.01 terlihat ketidakstabilan pada waktu setelah t=4. Sehingga hasil yang didapat menjadi tidak akurat.

0.2 ∆t = 0.001 ∆t = 0.01

Position (m)

0.15

0.1

0.05

0 0

1

2

3

4 Time (s)

5

6

7

8

Gambar 2.5: grafik posisi terhadap waktu untuk variasi ketinggian pada variasi ∆t Simulasi yang dilakukan menggunakan parameter dari hasil validasi kemudian memvariasikan nilai untuk perbedaan ketinggian, massa, dan jari-jari.

0.4 h = 0.4 m h = 0.3 m h = 0.2 m h = 0.1 m

0.35

0.3

Position (m)

0.25

0.2

0.15

0.1

0.05

0 0

1

2

3

4 Time (s)

5

6

7

8

Gambar 2.6: grafik posisi terhadap waktu untuk variasi ketinggian

6


Pada simulasi dengan variasi ketinggian seperti yang ditunjukkan pada Gambar 2.6, dapat dilihat bahwa bola dijatuhkan dari beberapa posisi yaitu 0.1m, 0.2m, 0.3m, dan 0.4m. Hasil dari simulasi dapat diketahui bahwa semakin rendah bola dijatuhkan semakin cepat waktu berhentinya. Kemudian berikutnya dilakukan simulasi dengan variasi massa dan dihasilkan grafik pada Gambar 2.7 0.2 m = 0.22 kg m = 0.022 kg m = 0.0022 kg m = 0.00022 kg

Position (m)

0.15

0.1

0.05

0 0

1

2

3

4 Time (s)

5

6

7

8

Gambar 2.7: grafik posisi terhadap waktu untuk variasi massa Pada simulasi dengan variasi massa, hasil simulasi menunjukkan semakin kecil massa bola yang dijatuhkan mengakibatkan semakin curam penurunan puncak pada pantulan berikutnya. Untuk simulasi yang ketiga, dilakukan dengan memvariasikan jari-jari bola dan didapatkan hasil pada Gambar 2.8 0.2 R = 0.05 m R = 0.03 m R = 0.02 m R = 0.01 m

Position (m)

0.15

0.1

0.05

0 0

1

2

3

4 Time (s)

5

6

7

8

Gambar 2.8: grafik posisi terhadap waktu untuk variasi jari-jari Pada simulasi dengan variasi jari-jari, hasil yang didapat menunjukkan bahwa semakin kecil jari-jari bola yang dijatuhkan mengakibatkan semakin tingginya pantulan berikutnya. Kemudian, semakin besar jari-jari bola akan mengakibatkan lebih curam penurunan puncak pada pantulan berikutnya.

7

2.5. RINGKASAN

h = 0.3 m h = 0.2 m h = 0.1 m

0.4

t = 4.68

t = 6.75

t = 7.8

Velocity (m/s)

0.2

0

-0.2

-0.4

4

4.5

5

5.5

6

6.5

7

7.5

8

Time (s)

Gambar 2.9: grafik kecepatan terhadap waktu Pada Gambar 2.9 menunjukkan perbedaan waktu berhenti dimana untuk bola yang dijatuhkan dari ketinggian 0.1m, waktu berhenti bola menunjukkan pada 4.68s. Sedangkan waktu berhenti bola yang dijatuhkan dari ketinggian 0.2m adalah 6.75s dan bola yang dijatuhkan dari ketinggian 0.3m adalah 7.8s. Dalam grafik tersebut terlihat bahwa setelah kecepatannya ≈ 0 terjadi perubahan kecepatan lagi secara konstan. Hal ini disebabkan karena pada simulasinya, drag force memiliki arah berlawanan dengan arah geraknya sehingga masih memiliki nilai yang mempengaruhi kecepatannya.

2.5

Ringkasan

Dalam simulasi bouncing ball ini telah dilakukan validasi dengan eksperimen menggunakan bola ping pong bermassa 0.0022kg dan jari-jari 0.02m dan dijatuhkan dari ketinggian 0.2m. Hasil validasi untuk simulasi valid hingga t ≈ 1s. Dalam simulasi bouncing ball dilakukan beberapa perlakuan, yaitu variasi massa bola, jari-jari bola, dan ketinggian bola saat dijatuhkan. Pada simulasi dengan variasi massa, hasil simulasi menunjukkan semakin kecil massa bola yang dijatuhkan mengakibatkan semakin curam penurunan puncak pada pantulan berikutnya. Hasil pada simulasi dengan variasi jari-jari bola menjukkan bahwa semakin kecil jari-jari bola yang dijatuhkan mengakibatkan semakin tingginya pantulan berikutnya, sedangkan semakin besar jari-jari bola akan mengakibatkan lebih curam penurunan puncak pada pantulan berikutnya. Dan untuk variasi ketinggian bola saat dijatuhkan menunjukkan bahwa semakin rendah bola dijatuhkan maka semakin cepat waktu berhentinya.

8


2.6

Referensi

1. Faith A, Morrison. 2013. ”Data Correlation for Drag Coefficient for Sphere”. Michigan Technological University, Houghton 2. Shiflet, Angela B. 2006. ”Introduction to Computational Science: modeling and simulation for sciences. Princeton Univercity Press: New Jersey 3. Mathwork. ”Simulation of a Bouncing Ball”. 10 Desember 2015. http://www.mathworks.com/help/simulink/examples/ simulation-of-a-bouncing-ball.html#zmw57dd0e284

2.7

Lampiran

2.7.1

Kode C++

/** * * Program : bounceSim.cpp * Description : Simulation of Bouncing Ball phenomena ... * * - This code use GLUT API to display animation * - This code works fine running in Ubuntu 15.04 using g++ ver4.9.2 * and OpenGL ver2.1 * * Compile : g++ bounceSim.cpp -o bounceSim -lGL -lglut -lGLU -lm * Execute : ./bounceSim * * Haerul Jusmar Ibrahim | [email protected] * */ using namespace std; #include #include #include #include #include #include #define #define #define #define

<sstream> <stdio.h> <string>

PI 4.0*atan(1.0) G -9.8f // Gravity LOSS 0.05f // Restitutiom Coeffient Cd 0.47f // Drag coefficient

// Medium’s Parameter

2.7. LAMPIRAN

9

const double L = 0.5f; // Square Wall : -L < x < L and 0 < y < L const double RHO = 1.29; // Air density 0.001293 gr/cm^3 -> 1.29 kg/m^3 // http://hyperphysics.phy-astr.gsu.edu/hbase/tables/density.html // Ball’s parameter double initPosition = 0.2; // 20cm -> 0.2m const double R = 0.02f; // Radius D:4cm R:2cm -> 0.02m const double MASS = 0.0022f; // Mass of Ball 2.2 g -> 0.0022 kg const double A = PI*R*R; // Area of circle const int N = 1; // Amount of ball const double dt = 0.001; // Time Step double t = 0.; // Time double max_t = 10.; // Maximum time of simulation int step_number = 0; // Number of Step FILE *file = fopen("dataa.txt","w"); class BouncingBall { public: double x, vx, y, vy, r, temp_x, temp_y; BouncingBall(double x, double vx, double y, double vy, double r) { this->x = x; this->vx = vx; this->y = y; this->vy = vy; this->r = r; } void step(double dt) { /// Euler Scheme // Velocity vx += -0.5*RHO*vx*fabs(vx)*Cd*A/MASS * dt; vy += ((G - 0.5*(RHO*vy*fabs(vy)*Cd*A)/MASS))*dt; // Position temp_x = x + vx* dt; temp_y = y + vy* dt; /// Collision Detection (Ball & Wall) // right if (temp_x + R >= L){ x = L - (temp_x + R - L); vx = -vx; }else{ x = temp_x; }

10


// left if (temp_x - R < 0){ x = abs(temp_x - R) + R; vx = -vx; }else{ x = temp_x; } // top if (temp_y + R > L){ y = L - (temp_y + R -L); vy = -vy; }else{ y = temp_y; } // bottom if (temp_y - R < 0){ y = abs(temp_y - R) + R; vy = -vy + LOSS*vy; }else{ y = temp_y; } } }; void title() { cout << "-------------------------------------------------" << endl; cout << " Bouncing Ball Simulation " << endl; cout << "-------------------------------------------------" << endl; cout << " Coef of Restitution \t: " << LOSS << endl; cout << " Coef of Drag \t: " << Cd << endl; cout << " Air Density \t: " << RHO << endl; cout << endl; cout << " dt : " << dt << endl; cout << " y0 : " << initPosition << " m" << endl; cout << " m : " << MASS << " Kg" << endl; cout << " R : " << R << " m" << endl; cout << endl; cout << " >> Click on GLUT window to begin simulation! <<"<< endl << endl; } BouncingBall **balls; void create_balls() { balls = new BouncingBall*[N]; double xR; double yR; // initialize position and velocity of balls for (int i = 0; i < N; i++) {

11

2.7. LAMPIRAN

xR = L*0.5; yR = initPosition + R; balls[i] = new BouncingBall(xR, 0 /*vx*/, yR, 0 /*vy*/, R); } } double EKin = 0., temp_EKin = 1.; // Simulation void Simulation() { cout.flush(); cout << " time : " << t << "

\r";

step_number++; t += dt; for (int i = 0; i < N; i++) { balls[i]->step(dt); } // calculate Energy Kinetic double VX, VY; EKin = 0.; for (int i = 0; i < N; i++) { VX = balls[i]->vx; VY = balls[i]->vy; EKin += 0.5*MASS*(VX*VX + VY*VY); } EKin = EKin/N; // write to file every 2 step if (step_number%2 == 0){ // this code to write parameter for one ball only fprintf(file,"%f %f %f %f\n", t, EKin, balls[0]->y - R,balls[0]->vy); glutPostRedisplay(); // update animation picture } double exitCriteria = temp_EKin = EKin;

fabs(temp_EKin - EKin);

double minPosition = fabs(balls[0]->y - R); // minimum // Terminate Conditions if ((exitCriteria < 1e-8 && minPosition <= 1e-4) || (t >= max_t)) { fclose(file); glutIdleFunc(NULL); cout.flush(); cout << "\n\n >> Press Enter key <<" << endl;

12

| 2. SIMULASI BOUNCING BALL cout.flush(); getchar();

exit(1); } } // Draw text in GLUT Display void drawText(const string &str, double x, double y) { glRasterPos2d(x, y); int len = str.find(’\0’); for (int i = 0; i < len; i++) glutBitmapCharacter(GLUT_BITMAP_HELVETICA_12, str[i]); } // Update Display void display() { glClear(GL_COLOR_BUFFER_BIT); glColor3ub(1, 1, 1); glLineWidth(3.); glBegin(GL_LINE_STRIP); glVertex2d(0, L); glVertex2d(L, L); glVertex2d(L, 0); glVertex2d(0, 0); glVertex2d(0, L); glEnd(); static GLubyte color_table[6][3] = {{255, 0, 0}, {0, 255, 0}, {0, 0, 255}, {0, 255, 255}, {255, 0, 255}, {255, 255, 0}}; // create circle with radius r for (int i = 0; i < N; i++) { glColor3ubv(color_table[i%6]); glBegin(GL_TRIANGLE_FAN); double phi = 2 * PI / 24; for (int j = 0; j < 25; j++) { glVertex2d(balls[i]->x + balls[i]->r * cos(phi*j), balls[i]->y + balls[i]->r * sin(phi*j)); } glEnd(); } glColor3ub(0, 0, 0); ostringstream os; os << "Time \t: " << t

<< ends;

13

2.7. LAMPIRAN drawText(os.str(), L*0.05, L*0.9); os.seekp(0); os << "y \t: " << balls[0]->y - R << ends; drawText(os.str(), L*0.05, L*0.8); os.seekp(0); glutSwapBuffers(); } // Glut callback function to reset the viewing transformation void reshape(int w, int h) { glViewport(0, 0, w, h); glMatrixMode(GL_PROJECTION); glLoadIdentity(); gluOrtho2D(-0.1, L+0.1, -0.1, L+0.1); glMatrixMode(GL_MODELVIEW); glLoadIdentity(); } bool running = false;

// is the animation running

// Glut mouse callback function void mouse(int button, int state, int x, int y) { switch (button) { case GLUT_LEFT_BUTTON: if (state == GLUT_DOWN) { if (running) { glutIdleFunc(NULL); running = false; } else { glutIdleFunc(Simulation); running = true; } } } } int main(int argc, char *argv[]) { system("clear"); title(); create_balls(); fprintf(file,"#t\tEKin\ty\tvy\n"); fprintf(file,"%f %f %f %f\n", t, EKin, balls[0]->y - R,balls[0]->vy); // initialize GLUT glutInit(&argc, argv); glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);

14

| 2. SIMULASI BOUNCING BALL glutInitWindowSize(300, 300); glutInitWindowPosition(100, 100); glutCreateWindow("Bouncing Ball Simulation"); glClearColor(1.0, 1.0, 1.0, 0.0); glShadeModel(GL_FLAT); glutDisplayFunc(display); glutReshapeFunc(reshape); glutMouseFunc(mouse); glutMainLoop(); return 0;

} Kode C++ di atas dikompilasi dengan cara g++ bounceSim.cpp -o bounceSim -lGL -lglut -lGLU -lm dan dipanggil dengan cara ./bounceSim

2.7.2

Script Gnuplot

set terminal postscript eps enhanced color font ’,14’ set output ’validasi.eps’ set grid set xlabel "Time (s)" font ",14" set ylabel "Position (m)" font ",14" set xrange [0:8] plot ’simulation.txt’ u 1:3 pt 6 lt rgb "red" title ’Simulation’, \ ’experiment.txt’ u 1:2 pt 3 lt rgb "blue" title ’Experiment’

|3

Model Lintasan Shuttlecock dengan Variasi Gaya Hambat Udara Azwar Sutiono | [email protected] Rahman Syam | [email protected] Bulutangkis adalah salah satu olahraga yang populer di Indonesia. Banyak aspek fisika yang dapat ditinjau dalam olahraga ini, salah satunya adalah analisis lintasan dari shuttlecock ketika berada di udara. Berbeda dengan lintasan peluru, shuttlecock tidak membentak lintasan parabola yang simetris. LungMing Chen et al [1] telah melakukan analisis lintasan shuttlecock dengan menggunakan metode analitik. Dalam simulasi komputasi, persamaan diferensial gerak shuttlecock diselesaikan secara numerik dengan menggunakan metode Runge-Kutta orde-4. Hasil numerik ini akan dibandingkan dengan hasil analisis secara analitik. Selain itu, dengan metode numerik akan ditambahkan variasi gaya hambat udara yang sulit dilakukan dengan metode analitik.

3.1

Teori

Ketika shuttlecock berada di udara, berdasarkan Hukum Kedua Newton : ΣF~ = m~a

(3.1)

~ + F~H = m~a W

(3.2)

~ , F~H , m, dan a berturut-turut adalah gaya gavitasi, gaya hambat dimana W udara, massa dan percepatan shuttlecock. Besar gaya hambat udara bergantung 1

2| 3. MODEL LINTASAN SHUTTLECOCK DENGAN VARIASI GAYA HAMBAT UDARA pada kecepatan shuttlecock di udara dan arahnya selalu berlawanan dengan arah dari shuttlecock. Secara umum, besar gaya hambat udara dapat dituliskan sebagai berikut : (3.3) F~H = −bv n

dimana v adalah kecepatan shuttlecock relatif terhadap udara, parameter n adalah bilangan riil dan b adalah konstanta berdasarkan sifat udara, bentuk dan dimensi shuttlecock. Nilai konstanta b dapat ditentukan oleh kecepatan terminal shuttlecock. Ketika shuttlecock jatuh bebas di udara, kecepatan dan gaya hambat akan bertambah seiring bertambahnya waktu. Percepatan shuttlecok akan menjadi nol ketika gaya hambat udara sebanding dengan berat shuttlecock. Pada saat itu, shutllecock mendapatkan kecepatan terminal vt dan selanjutnya bergerak dengan percepatan nol dan kecepatan konstan. Konstanta b dapat ditentukan dari persamaan (3.2) dengan a = 0 sehingga : ~ − bv n = m~a W

(3.4)

mg − bv n = m~a

(3.5)

n

mg − bvt = 0

(3.6)

vx0 = v0 cos θ, vy0 = v0 sin θ

(3.9)

n

(3.7) bvt = mg mg (3.8) b= n vt Asumsikan shuttlecock mempunyai kecepatan awal v0 , sehingga kecepatan horizontal dan vertikal dapat dituliskan,

Gaya hambat udara diuraikan dalam komponen horizontal dan vertikal F~H = FHxˆi + FHy ˆj

(3.10)

FHx = −bvx n , FHy = −bvy n

(3.11)

dimana Dari persamaan (3.2) dan (3.11) maka dapat dirumuskan persamaan gerak shuttlecock yaitu : a.) Komponen vertikal m

dvy = −mg − bvy n dt dvy b = −g − vy n dt m

b.) Komponen horizontal m

dvx = −bvx n dt

(3.12)

3

3.2. SYARAT BATAS b dvx = − vx n dt m

(3.13)

yang merupakan percepatan shuttlecock pada arah sumbu-x dan sumbu-y.

3.2

Syarat Batas

Kecepatan shuttlecock dibatasi hanya memiliki komponen pada arah x dan y ~v = vx eˆx + vy eˆy

(3.14)

Begitu pula untuk gaya hambat sebagaimana yang dituliskan pada persamaan (3.10). Kecepatan terminal shuttlecock didapatkan melalui eksperimen yang telah dilakukan oleh Pastreal et al [2] pada tahun 1980 dengan besar vt = 6.8 m/s.

3.3


Dengan menggunakan metode Runge-Kutta orde-4, persamaan (3.12) dan (3.13) dapat memberikan kecepatan dan posisi setiap saat untuk arah x ax (t) = −

b n vx m

(3.15)

k1x = ax

(3.16)

k2x = ax + k1x ∆t/2

(3.17)

k3x = ax + k2x ∆t/2

(3.18)

k4x = ax + k3x ∆t

(3.19)

vx (t + ∆t) = vx (t) + 1/6 (k1x + 2k2x + 2k3x + k4x ) ∆t

(3.20)

x(t + ∆t) = x(t) + vx (t)∆t

(3.21)

dan untuk arah y ay (t) = −g −

b n vy m

(3.22)

k1y = ay

(3.23)

k2y = ay + k1y ∆t/2

(3.24)

k3y = ay + k2y ∆t/2

(3.25)

k4y = ay + k3y ∆t

(3.26)

vy (t + ∆t) = vy (t) + 1/6 (k1y + 2k2y + 2k3y + k4y ) ∆t

(3.27)

y(t + ∆t) = y(t) + vy (t)∆t

(3.28)

4| 3. MODEL LINTASAN SHUTTLECOCK DENGAN VARIASI GAYA HAMBAT UDARA

Dalam implementasinya Persamaan-persamaan (3.15) - (3.28) dituangkan dalam algorima berikut ini: L01. L02. L02. L03. L04. L05. L06. L07. L08. L09. L10. L11. L12. L13. L14. L15. L16. L18. L19. L20. L21.

3.4

∆t, tbeg , tend ? m, g, vterm , v0 , θ, n? vx0 = v0 cos θ, vy0 = v0 sin θ b = mg / (vterm )n t = tbeg b vx n ax (t) = − m b ay (t) = −g − m vy n k1x = ax k2x = ax + k1x ∆t/2 k3x = ax + k2x ∆t/2 k4x = ax + k3x ∆t k1y = ay k2y = ay + k1y ∆t/2 k3y = ay + k2y ∆t/2 k4y = ay + k3y ∆t vx (t + ∆t) = vx (t) + 1/6 (k1x + 2k2x + 2k3x + k4x ) ∆t vy (t + ∆t) = vy (t) + 1/6 (k1y + 2k2y + 2k3y + k4y ) ∆t x(t + ∆t) = x(t) + vx (t)∆t y(t + ∆t) = y(t) + vy (t)∆t t = t + ∆t t ≤ tend → L04.

Hasil dan Diskusi

Nilai-nilai parameter yang digunakan dalam simulasi tercantum dalam Tabel 10.1. Tabel 3.1: Nilai-nilai parameter simulasi. Parameter ∆t tbeg tend m g θ v0 vt n

Nilai 10−2 0 10 0.005 kg 9.81 m/s2 30o 30 m/s 6.8 m/s 0.5, 1.0, 1.5, 2.0, 2.5, 3.0

Gambar 3.1 menunjukkan plot lintasan shuttlecock dengan gaya hambat udara yang bervariasi. Panjang lintasan terjauh shuttlecock untuk n = 0.5, 1.0, 1.5, 2.0, 2.5 dan 3.0 berturut-turut adalah 21.51 m, 16.80 m, 13.12 m, 10.93 m, 9.06

5


m dan 8.26 m. Dari data yang didapatkan, terlihat bahwa semakin besar nilai n (semakin besar gaya hambat udara) maka jarak terjauh yang dapat dicapai shuttlecock semakin pendek. Ketinggian maksimum shuttlecock dari nilai n paling kecil berturut-turut adalah 5.19 m, 4.77 m, 4.37 m, 4.00 m, 3.68 m dan 3.40 m. Hal ini juga menunjukkan bahwa semakin besar nilai n, maka ketinggian yang dapat dicapai oleh shuttlecock semakin rendah. Persamaaan polinomial untuk masing-masing nilai n didapatkan dengan menggunakan Microsoft Excel yang hasilnya dapat dilihat pada tabel 3.2.

6

n=0.5 n=1 n=1.5 n=2 n=2.5 n=3

y (m)

4

2

0

0

2

4

6

8

10

12

14

16

18

20

22

24

x (m) Gambar 3.1: Lintasan shuttlecock dengan variasi gaya hambat (n=0.5, 1.0, 1.5, 2.0, 2.5, 3.0)

Tabel 3.2: Persamaan polinomial untuk masing-masing nilai n. n 0.5 1.0 1.5 2.0 2.5 3.0

y y y y y y

Persamaan Polinomial = −0.0003x4 + 0.0086x3 − 0.1128x2 + 0.9844x − 0.3472 = −0.0008x4 + 0.0211x3 − 0.2057x2 + 1.2112x − 0.4351 = −0.0015x4 + 0.0297x3 − 0.2022x2 + 1.0533x − 0.2392 = −0.0017x4 + 0.0193x3 − 0.0759x2 + 0.7374x − 0.0516 = −0.0040x4 + 0.0467x3 − 0.2015x2 + 1.0434x − 0.1636 = −0.0033x4 + 0.0264x3 − 0.1062x2 + 0.9771x − 0.1021

Gambar 3.2 menunjukkan perbandingan antara plot lintasan shuttlecock hasil analitik dan numerik untuk nilai n=1 dengan persamaan lintasan, y=x

v2 vt + vy0 vt vx0 − t ln vx0 g vt vx0 − gx

(3.29)

Terjadi sedikit perbedaan antara panjang lintasan hasil analitik dan numerik namun dengan beda yang masih dapat ditoleransi. Perbandingan ini dapat dilihat pada tabel 3.3.


6

n=1 komputasi n=1 analitik

y (m)

4

2

0

0

2

4

6

8

10

12

14

16

18

x (m) Gambar 3.2: Perbandingan lintasan hasil analitik dan numerik untuk n=1 Tabel 3.3: Perbandingan jarak horizontal hasil analitik dan numerik n 1.0

3.5

Analitik (m) 17.17

Numerik (m) 16.80

Ringkasan

Telah dilakukan simulasi lintasan shuttlecock di udara menggunakan metode Runge-kutta orde-4 dengan variasi gaya hambat untuk n=0.5, 1.0, 1.5, 2.0, 2.5, 3.0. Hasil simulasi menunjukkan bahwa semakin besar nilai n, maka jarak terjauh yang dapat dicapai oleh shuttlecock semakin pendek dan ketinggian maksimum shuttlecock semakin rendah. Hasil tersebut kemudian dibandingkan dengan perhitungan secara analitik untuk n=1. Kelebihan dari metode numerik adalah kita dapat mengatur gaya hambat dengan n yang besar (n > 1) dan bukan bilangan bulat (n = 0.5, 1.5, 2.5, dst), dimana akan sangat rumit jika dihitung secara analitik.

3.6

Referensi

1. Lung-Ming Chen, Yi-Hsiang Pan, and Yung-Jen Chen. (2009) A Study of Shuttlecock’s Trajectory in Badminton, Journal of Sport Science and Medicine 8, 657-662 2. Peastrel, M., Lynch, R. and Angelo, A. Jr. (1980) Terminal Velocity of Shuttlecock in Vertical Fall, American Journal of Physics 48 (7), 511-513 3. Atam P. Arya. (1997) Introduction to Classical M echanics 2nd Edition, Benjamin Cummings.

7

3.7. LAMPIRAN

3.7

Lampiran

3.7.1

Kode C++

//Program Lintasan Shuttlecock n=0.5 #include #include #include <math.h> using namespace std; int main(int argc, char *argv[]) { //parameter awal double m = 0.005; double g = 9.81; double v0 = 30.0; double vterm = 6.8; double theta = 30;

// massa shuttlecock = 0.005 kg // g = 9.81 m/s^2 // kecepatan awal = 30 m/s // kecepatan terminal // sudut elevasi

// Kondisi awal double vx = v0*cos(theta*M_PI/180.00); // kecepatan sumbu-x double vy = v0*sin(theta*M_PI/180.00); // kecepatan sumbu-y double x = 0; double y = 0; // Parameters related to simulation time double tbeg = 0; // waktu awal double tend = 10; // waktu akhir double dt = 0.01; // time step double t = tbeg; // data output ofstream out("data_xy_n=0.5.txt"); // Show header of data cout << "# t\tx\ty" << endl; // Mengatur nilai konstanta b double b=m*g/sqrt(vterm); // Simulasi while(t <= tend) { // Hasil cout << t << "\t"; cout << x << "\t"; cout << y << endl;

8| 3. MODEL LINTASAN SHUTTLECOCK DENGAN VARIASI GAYA HAMBAT UDARA out << x << "\t"; out << y << endl; // percepatan double ax = - (b/m)*sqrt(vx); double ay = - ( g + (b/m)*sqrt(fabs(vy)) ); // Runge-Kutta 4th order double k1x = ax; double k2x = ax + k1x * dt / 2; double k3x = ax + k2x * dt / 2; double k4x = ax + k3x * dt; double double double double

k1y k2y k3y k4y

= = = =

ay; ay + k1y * dt / 2; ay + k2y * dt / 2; ay + k3y * dt;

// Menghitung kecepatan vx = vx + (k1x + 2*k2x + 2*k3x + k4x) * dt / 6; vy = vy + (k1y + 2*k2y + 2*k3y + k4y) * dt / 6; // Menghitung posisi x = x + vx * dt; y = y + vy * dt; // increase time t += dt; } return 0; }

Kode C++ di atas dikompilasi dengan cara g++ lintasan_n=0.5.cpp -o lintasan-n=0.5 dan dipanggil dengan cara ./lintasan-n=0.5

3.7.2

Script Gnuplot

# Set terminal to eps set term post eps color enhanced

3.7. LAMPIRAN # Set output file set output "gambar_n=0.5.eps" # Set input file input="data_xy_n=0.5.txt" # Set image size set size 1, 0.67 # Set bottom margin -- for some distro it needs correction set bmargin 3.5 # Set x-axis set xlabel "{/Italic x} (m)" font "Times, 24" set xtics 2 font "Times, 22" set xrange [0:24] # Set y-axis set ylabel "{/Italic y} (m)" font "Times, 24" set ytics 2 font "Times, 22" set yrange [0:8] # Set grid set grid xtics ytics # Plot data plot input u 1:2 title "n=0.5" pt 4 ps 2.5 Script Gnuplot di atas dipanggil dengan cara gnuplot plot-n=0.5.gps

9


|4

Simulasi Gerakan Air Dangkal 1 Dimensi Muhammad Mirza Fahmi | [email protected] Persamaan diferensial parsial sering muncul dalam sejumlah permasalahan fisis, diantaranya aliran fluida, perpindahan panas, mekanika benda padat dan proses biologis. Salah satunya adalah permasalahan gerak air dangkal.Persamaan ini muncul dari persamaan dasar mekanika fluida untuk fluida inviscid dan tak termampatkan.

4.1

Teori

Persamaan air dangkal bermanfaat untuk menjelaskan sebuah lapisan tipis fluida dengan kerapatan konstan dalam keseimbangan hidrostatik. Lapisan tipis ini dibatasi dari bawah oleh topografi dasar air dan dari atas oleh permukaan bebas. Karena hanya 1 dimensi, persamaan ini dapat dituliskan sebagai berikut

∂h ∂(uh) + =0 ∂x ∂t 1 ∂ ∂(uh) u2 h + gh = 0 + ∂t ∂x 2

(4.1) (4.2)

h merupakan ketinggian air terhadap permukaan, u adalah kecepatan horizontal fluida, dan g adalah percepatan gravitasi. Dalam dinamika fluida, persamaan dasar fluida terdiri dari kekekalan massa dan kekekalan momentum. Kekekalan massa pada persamaan air dangkal dtuliskan oleh persamaan (4.1), sedangkan kekekalan momentum dituliskan oleh persamaan (4.2). 1

2

| 4. SIMULASI GERAKAN AIR DANGKAL 1 DIMENSI

Dua persamaan (4.1) dan (4.2) dapat dituliskan dalam bentuk matriks sebagai berikut uh h =0 (4.3) + 2 uh t u h + 21 gh2 x Karena persamaan air dangkal merupakan persamaan diferensial parsial, pendekatan secara numerik yang pertama kali harus dilakukan adalah diskritisasi. Diskritisasi adalah proses konversi persamaan diferensial partial dan syarat bantu (syarat awal dan syarat batas) kedalam sebuah sistem aljabar diskrit. Proses diskritisasi dapat diidentifikasi melalui beberapa metode yang saat ini sering digunakan. Dalam hal ini kami akan menggunakan metode beda hingga dengan memanfaatkan skema Lax-Wendroff yang akan dibahas pada bab berikutnya.

4.2

Syarat batas

Dalam menyelesaikan persaman diferensial parsial ini, ada beberapa pilihan syarat batas yang sudah cukup umum dipakai. Dalam hal ini, kami memilih syarat batas free boundary conditions (FBC) yang berlaku untuk ketinggian air h dan kecepatan massa air uh. Syarat batas h dapat ditulis sebagai berikut,

h0 = h1 = hnx−2

(4.4) (4.5)

uh0 = uh1

(4.6)

uhnx−1 = uhnx−2

(4.7)

hnx−1 kemudian syarat batas uh

syarat batas ini memungkinkan untuk gelombang yang terpropagasi dapat dilanjutkan geraknya tanpa dipantulkan kembali.

4.3


Penyelesaian persamaan diferensial parsial salah satunya adalah dengan mendiskritisasi menggunakan metode beda hingga. Karena persamaan air dangkal ini adalah persamaan hiperbolik maka pendekatan yang tepat yaitu memanfaatkan skema Lax-Wendorf yang dibangun berdasarkan metode beda hingga. Nama skema ini diambil dari nama belakang penemunya yaitu Peter Lax dan Burton Wendorf.


3

Tinjau sebuah permasalahan memiliki persamaan sebagai berikut ∂u(x, t) ∂f (u(x, t)) + =0 ∂t ∂x

(4.8)

dimana x dan t adalah variabel bebas dan diberikan keadaan awal u(x, 0). Skema Lax-Wendorf terbagi dua langkah, langkah pertama

n+1/2

1 n ∆t (u + uni ) − (f (uni+1 ) − f (uni )), 2 i+1 2 ∆x ∆t 1 = (uni + uni−1 ) − (f (uni ) − f (uni−1 )). 2 2 ∆x

ui+1/2 = n+1/2

ui−1/2

(4.9) (4.10)

kemudian langkah kedua = uni − un+1 i

i ∆t h n+1/2 n+1/2 f (ui+1/2 ) − f (ui−1/2 ) . ∆x

(4.11)

f (u) merupakan fluks massa. Implementasinya pada persamaan (4.1) dan (4.2) ditulis sebagai berikut, untuk langkah pertama

λ 1 n (hi+1 + hni ) − (uhni+1 − uhni ), 2 2 n 2 2 1 1 λ (uh ) i = (uhni + uhni−1 ) − + · g · hni+1 2 2 hi+1 2 2 1 (uhi ) 2 − · g · (hi ) . − hi 2 n+1/2

hi+1/2 =

n+1/2

uhi−1/2

(4.12)

(4.13)

selanjutnya langkah kedua

uhn+1 i

h i n+1/2 n+1/2 hn+1 = hni − λ uhi+1/2 − uhi−1/2 , i 2 uhn+1/2 2 i+1/2 1 n+1/2 =uhni − λ + · g · h i+1/2 n+1/2 2 hi+1/2 2 n+1/2 2 uhi−1/2 1 n+1/2 − − · g · h i−1/2 n+1/2 2 h

(4.14)

(4.15)

i−1/2

selanjutnya implementasi skema dalam algoritma dapat ditulis sebagai berikut: S01. Atur posisi node posisi x dan waktu t S02. dx = xlength /(nx − 1);dt = tlength /(nt); λ = dt/dx

4


S03. Syarat awal: h0i = 2 + cos(2πxi ); uh0i = 0 S04. Syarat Batas (t0 ): h0 = h1 ; hnx−1 = hnx−2 ; uh0 = uh1 ; uhnx−1 = uhnx−2 n+1/2

n+1/2

S05. Langkah Pertama: hi+1/2 ; uhi−1/2

S06. Langkah Kedua: hn+1 ; uhn+1 i i S07. Syarat Batas (tn ): h0 = h1 ; hnx−1 = hnx−2 ; uh0 = uh1 ; uhnx−1 = uhnx−2

4.4

Hasil dan diskusi

Nilai-nilai parameter yang digunakan dalam simulasi tercantum dalam Tabel 10.1. Tabel 4.1: Nilai-nilai parameter simulasi. Parameter nx nt xlength tlength g

Nilai 41 500 1 0.4 9.8

Dari parameter-parameter diatas, akan diperoleh nilai λ sebagai syarat konvergensi dari simulasi ini. Secara teori, syarat ini disebut sebagai syarat CourantFriedrichs-Lewy (CFL). Dalam Kasus satu dimensi, CFL memenuhi C=

uδt δx

(4.16)

dimana u adalah besar kecepatan, sedangkan δt dan δx berturut-turut adalah besar langkah dalam domain waktu dan spasial. Agar komputasi mencapai konvergensi, maka nilai C ≤ Cmax . Nilai Cmax bergantung dari metode yang digunakan (eksplisit atau implisit). Untuk metode yang eksplisit maka Cmax = 1. Salah satu hasil yang diperoleh adalah seperti dalam Gambar 10.1. Simulasi diawali dengan dijatuhkannya sebuah tetesan air dalam bentuk fungsi 2 + cos(2πxi ) sebagai syarat awal. Karena propagasi gelombang tetap dilanjutkan dan tidak dipantulkan kembali maka pada saat t = 0.3 keatas permukaan air sudah kembali tenang. Hasil berikutnya adalah keadaan dari fluks massa selama simulasi berlangsung seperti dalam Gambar 10.2. Kalau dibandingkan dengan Gambar 10.1 maka hasil simulasi memiliki akurasi yang cukup baik meskipun tampak tidak terlalu halus.

4.5

Ringkasan

Permasalahan fisis yang melibatkan persamaan diferensial parsial sangat sering kita temui. Khususnya permodelan air dangkal yang bermanfaat untuk memo-

5

4.6. REFERENSI t=0

t=0.0664

4

4

3

3 Tinggi Air

5

Tinggi Air

5

2

2

1

1

0

0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

Posisi

Posisi

t=0.0856

t=0.1392

4

4

3

3

0.8

1

0.6

0.8

1

Tinggi Air

5

Tinggi Air

5

0.6

2

2

1

1

0

0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

Posisi

Posisi

t=0.168

t=0.2408

4

4

3

3 Tinggi Air

5

Tinggi Air

5

2

2

1

1

0

0 0

0.2

0.4

0.6 Posisi

0.8

1

0

0.2

0.4

0.6

0.8

1

Posisi

Gambar 4.1: Simulasi air dangkal dengan syarat awal 2 + cos(2πxi ) dan nilai parameter: nt = 500, nx = 41, g = 9.8, xlength = 1, dan tlength = 0.4 delkan berbagai macam permasalahan seperti Tsunami dan Tanggul yang jebol. Simulasi persamaan fluida satu dimensi ini akan bergantung dari nilai λ, karena nilai ini akan menentukan tingkat konvergensi dari simulasi yang akan dikerjakan.

4.6

Referensi

1. P. Crowhurst, ”Numerical Solutions of One-Dimensional Shallow Water Equations”, Australian Mathematical Mathematical Sciences Institute (2013). 2. Richard L. Burden and J. Douglas Fairies, ”Numerical Analysis”, Brooks/Cole, 9th (2011) 3. L. Rezolla, ”Numerical Methods for The Solution of Partial Differential Equation”, Lecture Notes, (2011) 4. C. Hirsch, ”Numerical Computation of Internal and External Flows”, Elsevier, 2nd (2007) 5. C. Moler, ”Experiments with Matlab”, Mathworks.

6

| 4. SIMULASI GERAKAN AIR DANGKAL 1 DIMENSI t=0.0664

4

2

2

Tinggi Air

Tinggi Air

t=0

4

0

-2

-2

-4

-4

0.2

0.4

0.6

0.8

1

0

0.2

0.4

Posisi

Posisi

t=0.0856

t=0.1392

4

4

2

2

Tinggi Air

Tinggi Air

0

0

-2

0.6

0.8

1

0.6

0.8

1

0

-2

-4

-4

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

Posisi

Posisi

t=0.168

t=0.2408

4

4

2

2

Tinggi Air

Tinggi Air

0

0

-2

0

-2

-4

-4

0

0.2

0.4

0.6

0.8

Posisi

1

0

0.2

0.4

0.6

0.8

1

Posisi

Gambar 4.2: Keadaan fluks massa uh pada simulasi air dangkal

4.7

Lampiran

4.7.1

Kode C++

#include #include <stdio.h> #include <math.h> #include #include #include #include #define PI 3.141592653589793 using namespace std; int i, j, iter, nx, nt; double lambda, g, delta_x, delta_t, dx, dt; double *h, *h_array, *ht, *uh, *uh_array, *uht, *x, *t; int main() {

4.7. LAMPIRAN

7

FILE *pipe = popen("gnuplot -persist", "w"); /*Parameter Simulasi*/ nx = 41; nt = 500; delta_x = 1; delta_t = 0.4; g = 9.8; /*Alokasi Memory*/ h = new double[nx]; h_array = new double[nx*(nt+1)]; ht = new double[nx-1]; t = new double[nt+1]; uh = new double[nx]; uh_array = new double[nx*(nt+1)]; uht = new double[nx-1]; x = new double[nx]; x = r8vec_linspace_new ( nx, 0, delta_x ); t = r8vec_linspace_new ( nt + 1, 0, delta_t ); /*Time and spasial step*/ dt = delta_t / (double) (nt); dx = delta_x / (double) (nx - 1); lambda = dt / dx; /*Terapkan syarat awal dan syarat batas saat t[0]*/ syarat_awal(nx, nt, x, t[0], h, uh); syarat_batas(nx, nt, x, t[0], h, uh); /*Simpan array pada time step awal*/ for (i = 0; i < nx; i++) { h_array[i + 0 * nx] = h[i]; uh_array[i + 0 * nx] = uh[i]; } /*Proses pemanggilan Gnuplot untuk menampilkan hasil simulasi*/ fprintf(pipe, "set terminal postscript eps enhanced color font \"Helvetica,14\"\n"); /*Mencetak hasil simulasi tiap time step dalam format eps*/ //fprintf(pipe, "set terminal x11 font \"Helvetica,14\"\n");

fprintf(pipe, "set yrange[-5:5] \n"); fprintf(pipe,"set xlabel \"Posisi\"\n"); fprintf(pipe,"set ylabel \"Tinggi Air\"\n"); fprintf(pipe,"unset key \n"); fprintf(pipe,"set grid \n"); for (iter = 0; iter <= nt; iter++) {

8


/*Perhitungan Langkah Pertama */ for (i = 0; i < nx - 1; i++) { ht[i] = 0.5 * (h[i] + h[i + 1]) - 0.5 * lambda * (uh[i + 1] - uh[i]); } for (i = 0; i < nx - 1; i++) { uht[i] = 0.5 * (uh[i] + uh[i + 1]) - 0.5 * lambda * (pow(uh[i + 1], 2) / h[i + 1] + 0.5 * g * pow(h[i + 1], 2) - pow(uh[i], 2) / h[i] - 0.5 * g * pow(h[i], 2)); } /*Perhitungan Langkah Kedua for(i = 1; i < nx - 1; i++){ h[i] = h[i] - lambda * (uht[i] - uht[i-1]); } for (i = 1; i < nx; i++){ uh[i] = uh[i] - lambda * (pow(uht[i],2)/ht[i] + 0.5 * g * pow(ht[i],2) - pow(uht[i-1],2)/ht[i-1] - 0.5 * g * pow(ht[i-1],2)); } /*Terapkan syarat batas untuk time step selanjutnya*/ syarat_batas(nx, nt, x, t[iter], h, uh); /*Salin data dalam array yang besar*/ for (i = 0; i < nx; i++){ h_array[i+iter*nx] = h[i]; uh_array[i+iter*nx] = uh[i]; } /*Mencetak luaran simulasi dalam format eps */ fprintf(pipe,"set output \"fig-%d.eps\"\n",iter); fprintf(pipe,"set title \"t=%g\"\n",t[iter]); fprintf(pipe, "plot ’-’ with lines\n"); for(i = 0; i < nx; i++){ fprintf(pipe, "%f %f\n", x[i], uh[i]); } fprintf(pipe, "%s\n", "e"); fflush(pipe); cout << "\b\riter = " << iter; cout.flush(); } cout << endl; /*Free Memory*/

4.7. LAMPIRAN

9

free (h); free (h_array); free (ht); free (t); free (uh); free (uh_array); free (uht); free (x); waktu(); return 0; } void syarat_awal(int nx, int nt, double x[], double t, double h[], double uh[]) { for (i = 0; i < nx; i++) { h[i] = 2 + cos ( 2 * PI * x[i] ); } for (i = 0; i < nx; i++) { uh[i] = 0; } return; } void syarat_batas(int nx, int nt, double x[], double t, double h[], double uh[]){ h[0] = h[1]; h[nx-1] = h[nx-2]; uh[0] = uh[1]; uh[nx-1] = uh[nx-2]; return; } double *r8vec_linspace_new ( int n, double a_first, double a_last ) { double *a; int i; a = new double[n]; if ( n == 1 ) { a[0] = ( a_first + a_last ) / 2.0; } else { for ( i = 0; i < n; i++ ) {

10


a[i] = ( ( double ) ( n - 1 - i ) * a_first + ( double ) ( i ) * a_last ) / ( double ) ( n - 1 ); } } return a; } void waktu(){ # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct std::tm *tm_ptr; size_t len; std::time_t now; now = std::time ( NULL ); tm_ptr = std::localtime ( &now );

len = std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr ) std::cout << time_buffer << "\n"; return; # undef TIME_SIZE } Kode C++ di atas dikompilasi dengan cara g++ code.cpp -o code dan dipanggil dengan cara ./code

|5

Analisis Titik Luluh Material Dengan Menggunakan Metode Secant Fitriyanti Nakul | [email protected] Mega Silvia Lestari | [email protected] Linda Handayani | [email protected] Deni | [email protected]

5.1

Pendahuluan

Bahan logam maupun non logam memiliki karakteristik mekanik yang menentukan kualitas bahan tersebut. Karakteristik ini berkaitan dengan kemampuan bahan untuk menahan beban luar yang mengenainya. Beban atau gaya luar yang diberikan akan cenderung menimbulkan perubahan bentuk dan ukuran (deformasi) ketika mencapai batasan nilai tertentu [1]. Untuk menghindari masalah yang akan ditimbulkan dari perubahan bentuk tersebut maka perlu dilakukan pengujian karakteristik material melalui pengontrolan dan desain yang tepat, termasuk pentingnya menentukan nilai batasan saat bahan mulai berdeformasi yang direpresentasikan oleh titik luluh material. Pada penelitian ini, kekuatan bahan dideskripsikan sebagai kekuatan ikatan antar partikel dari bahan tersebut. Ikatan antar partikel terdekat didekati dengan pendekatan partikel-pegas bergandeng. Kondisi x ¨=0 dipilih. Penyelesaian numerik dengan metode secant digunakan untuk mendapatkan nilai konstanta dari persamaan Hollomon dan menghubungkan persamaan tersebut dengan per1

2| 5. ANALISIS TITIK LULUH MATERIAL DENGAN MENGGUNAKAN METODE SECA saman Hooke.

5.1.1

Tujuan

Tujuan dari penelitian ini ialah menentukan nilai titik luluh dari beberapa material dengan menggunakan metode secant. Secara spesifik tujuan ini meliputi: 1. Menganalisis hubungan regangan dan tegangan bahan yang muncul akibat gaya yang diberikan. 2. Menganalisis pengaruh variasi konstanta pegas terhadap deformasi pada bahan.

5.1.2

Variasi Parameter

Dalam penelitian ini terdapat dua variasi parameter, antara lain: variasi material yang digunakan yaitu alumunium dan tembaga dan variasi nilai konstanta pegas. Konstanta pegas merupakan besaran yang mempengaruhi besarnya gaya yang diberikan.

5.2

Teori

Titik luluh adalah parameter yang merepresentasikan suatu kondisi saat material tidak lagi berada dalam kondisi elastis, dengan kata lain ketika pembebanan dihilangkan, benda tersebut tidak dapat kembali ke bentuk semula. Titik maksimum ini menjadi batasan antara daerah I (linear/elastis) dan daerah II (plastis), saat mulai terjadi deformasi seperti tampak pada Gambar 5.1.

Gambar 5.1: Hubungan regangan dan tegangan pada bahan. [2]

3

5.2. TEORI

Berdasarkan Gambar 5.1, secara kuantitatif untuk daerah linier (daerah I) berlaku persamaan yang mengikuti Hukum Hooke sedangkan pada daerah plastis (daerah II) berlaku persamaan Hollomon [3] berikut: σ = Kǫn + C

(5.1)

Beberapa parameter penting yang berelasi dengan fokus kajian dari sistem yang akan ditinjau meliputi: Tegangan (σ) adalah besaran pengukuran intensitas gaya (F ) atau reaksi dalam yang timbul persatuan luas (A). Sedangkan regangan (ǫ) merupakan perubahan geometri panjang awal sebagai hasil dari gaya yang menarik atau menekan material. Kemudian, K dan n sebagai konstanta empirik yang diperoleh dari hubungan tegangan dan regangan sebelum terjadi necking, dan C merupakan konstanta berdimensi (N m−2 ). Selain menghindari masalah yang dapat dimunculkan dari adanya deformasi akibat pembebanan yang berlebihan, sifat plastisitas material juga dapat dimanfaatkan untuk mendesain bentuk material sesuai yang diinginkan. Hal ini tentu bergantung pada titik luluh dari material yang digunakan. Setiap material memiliki nilai titik luluh yang berbeda, ini terlihat pada Tabel 5.1 berikut:

Tabel 5.1: Nilai titik luluh untuk material yang berbeda Modulus Young Tegangan luluh Regangan luluh (1010 N/m2 ) (106 N/m2 ) Aluminium 10 28 0.00028 Tembaga 11 69 0.000672 Material

n

K (10 N/m2 ) 6.4 3.15 8

0.2 0.54

[4] Dalam penelitian ini, penentuan titik luluh dilakukan dengan menggunakan model sistem pegas bergandeng (1 partikel, 2 pegas) seperti diilustrasikan pada Gambar 5.2. Dengan menerapkan persamaan Euler-lagrange pada sistem di atas, maka diperoleh persamaan gerak sistem sebagai berikut. m1 x ¨1 + (k1 + k2 )x1 − k2 x2 = 0

(5.2)

Gerak pada sistem ditinjau sebagai gerak statis, maka x ¨1 = 0, sehingga: x2 =

k + k 1 2 x1 k2

(5.3)

4| 5. ANALISIS TITIK LULUH MATERIAL DENGAN MENGGUNAKAN METODE SECA

Gambar 5.2: Model sistem pegas bergandeng

Pada kasus ini, perlu ditentukan fungsi persamaan untuk daerah I dan fungsi persamaan untuk daerah II. Hal ini dilakukan untuk mengetahui di titik mana bahan mulai berdeformasi. (a) Fungsi untuk daerah I 1. Nilai regangan dari sitem diberikan, ǫ=

x1 + x2 l

(5.4)

subtitusi persamaan (5.3) ke (5.4) diperoleh, h k1 + k 2 i x1 ǫ= 1+ k2 l h k + 2k i x 1 2 1 = k2 l

(5.5)

2. Nilai tegangan dari sitem diberikan, σ=

kx1 + kx2 A

(5.6)

subtitusi persamaan (5.3) ke (5.6), maka:

σ=

kx1 + k2

h

k1 +k2 k2

A h 2k + k i 1 2 x1 = A

x1

i

(5.7)

Dari ungkapan tegangan dan regangan yang diperoleh diatas maka didapatkan

5

5.3. SYARAT BATAS fungsi tegangan terhadap regangan yang memenuhi daerah I sebagai berikut, 1 σ h 2k1 + k2 i i x1 h = k1 +2k2 x1 ǫ A k2

l

σ h 2k1 + k2 i lk2 = ǫ A k1 + 2k2 h 2k + k i lk 2 1 2 ǫ σ= A k1 + 2k2 h 2k + k i lk 1 2 2 σ = mǫ, m= A k1 + 2k2

(5.8)

(b) Fungsi untuk daerah II Fungsi tegangan pada daerah II, mengikuti persaaman Holomon yang dinyatakan pada persamaan (5.1). Fungsi ini mengandung suatu konstanta C. Jika dianalisis secara numerik, persamaan tersebut membutuhkan suatu konstanta (C). Konstanta ini perlu ditambahkan agar dua persamaan yang diwakilkan pada dua daerah yang berbeda (daerah linier dan daerah plastis) tersebut dapat menyatu dalam satu titik tertentu tanpa menganggu makna fisis dan matematisnya.

5.3

Syarat Batas

Penentuan titik potong yang menghubungkan Persamaan Hooke dan Hollomon diperoleh dengan menerapkan syarat batas kontinuitas yang harus terpenuhi pada dua daerah yang ditinjau, yaitu σ1 = σ 2 σ2 dσ1 = dǫ1 dǫ2

(5.9)

sehingga diperoleh dσ1 dσ2 = dǫ dǫ m = nKǫn−1 m K = n−1 nǫ 1 m n−1 ǫ= nK σ1 = σ2 mǫ = Kǫn + C

(5.10)

(5.11)


Dengan mensubtitusikan parameter K dan ǫ pada persamaan (5.10) ke persamaan(5.11) diperoleh nilai konstanta C secara analitik sebagai berikut; mǫn−1 + Cǫ−1 nǫn−1 m + Cǫ−1 m= n 1 m 1−n m m= +C n nK 1 m m n−1 C = m− n nK

m=

(5.12)

Selain itu, syarat peralihan nilai regangan dari daerah I (linier) ke daerah II (plastis) ditentukan dengan menerapkan kondisi yang harus dipenuhi sebagai berikut: (a) jika ǫ ≤ ǫintersect −→ (daerah linier) (b) jika ǫ ≥ ǫintersect −→ (daerah plastis)

5.4


Penyelesaian secara numerik diimplementasikan melalui penggunaan metode secant, Adapun metode secant dimanfaatkan untuk menentukan nilai akar-akar dari sebuah fungsi [5]. Secara umum, bentuk metode secant diberikan melalui ungkapan, f (xn )[xn − xn−1 ] (5.13) xn+1 = xn − f (xn ) − f (xn−1 ) Sehingga, pada kasus ini nilai konstanta C juga dapat ditentukan secara numerik dengan menerapkan metode secant pada persamaan (5.14) berikut. F (c) = m −

1 m 1−n m −C n nK

(5.14)

Dalam implementasinya persamaan (5.8),(5.10),(5.12) dan (5.14) dituangkan dalam algoritma berikut: L01. Input parameter k1 , k2 , A, l, K, hn (variasi i nilai k1 dan k2 ) 2 L02. Definisikan parameter m, m = 2k1A+k2 k1lk +2k2 L03. Ambil kondisi syarat awal, x1 = 0, x2 = 0, ∆x = 0,000001, σ = 0, ǫ=0 2 x1 , xtot = x1 + x2 L04. Mulai proses hitung, x1 = x1 + ∆x1 , x2 = k1k+k 2 xtot L05. Hitung, ǫ1 = l , σ1 = mǫ 1 n−1 m L06. Hitung titik luluh, ǫintersect = nK

5.5. HASIL DAN DISKUSI L07. Ca = m −

m n

L08. F (C) = m −

L09. L10. L11. L12.

5.5

m n

m nK

7

1 n−1

−→ (berdimensi N m−2 ), analitik 1 1−n m − C nK , perhitungan numerik untuk C

n )(Cn −Cn−1 ) Gunakan metode secant, Cn+1 = Cn − ff(C (Cn )−f (Cn−1 ) Terapkan syarat batas, jika ǫ ≤ ǫintersect −→ (daerah linier) Terapkan syarat batas, jika ǫ ≥ ǫintersect −→ (daerah plastis) Output, plot grafik hubungan regangan dan tegangan bahan.

Hasil dan Diskusi

Nilai-nilai parameter yang digunakan pada kasus ini tercantum dalam Tabel 5.2

Tabel 5.2: Nilai-nilai parameter Material A(m2 ) L(m) n K(N/m) Aluminium 1.5 × 10−5 2.5 0.2 6.4 × 108 Tembaga 1.6 × 10−5 0.95 0.54 3.15 × 108 Secara spesifik, pada penelitian ini yang dibahas adalah nilai regangan luluh Aluminium dan Tembaga dengan variasi konstanta pegas. Regangan luluh Aluminium dengan konstanta pegas k1 = 1, 65 × 105 N/m dan k2 = 1, 3 × 105 N/m adalah 0, 0043858. Artinya, pada saat Alumunium meregang hingga 0, 0043858. Aluminium mulai mengalami deformasi plastis. Kurva deformasinya dapat dilihat pada Gambar 5.3.

Gambar 5.3: Perbandingan tegangan dan regangan aluminium

8| 5. ANALISIS TITIK LULUH MATERIAL DENGAN MENGGUNAKAN METODE SECA Variasi parameter dilakukan dengan mengubah nilai konstanta pegas untuk melihat posisi tegangan luluhnya. Konstanta pegas yang digunakan secara berurut terdapat pada Tabel 5.3

Variasi ke I II III

Tabel 5.3: Variasi parameter k1 (N/m) k2 (N/m) Regangan luluh 5 1.65 × 10 1.30 × 105 0.00148358 5 1.5 × 10 1.15 × 105 0.00171001 1.35 × 105 1 × 105 0.00200796

Perubahan variasi parameter nilai konstanta pegas memberikan perubahan terhadap titik dimana bahan mulai berdeformasi. Semakin besar konstanta pegas, menyebabkan gaya yang mengenai bahan juga besar. Hal inilah yang kemudian menyebabkan Aluminium dengan konstanta pegas yang lebih besar lebih cepat berdeformasi jika dibandingkan dengan nilai konstanta pegas yang lebih kecil. Pada variasi kedua, bahan mulai berdeformasi saat meregang 0, 00171001, yakni 15, 26% lebih panjang dari variasi pertama. Variasi ketiga, dengan nilai konstanta yang lebih keci lagi, bahan mulai berdeformasi saat meregang 0.00200796. Nilai regangan ini 35% lebih panjang dari variasi pertama, dan 20, 08% lebih panjang dari variasi kedua. Perbedaan antara 3 variasi tersebut, dapat digambarkan pada Gambar 5.4.

Gambar 5.4: Perbandingan Tegangan dan Regangan Aluminium dengan Beberapa Variasi Konstanta Pegas

Selain Alumunium, analisis juga dilakukan pada Tembaga. Pada kasus tembaga ini, nilai konstanta pegas yang digunakan adalah k1 = 2, 89 × 105 N/m dan k2 = 2, 79 × 105 N/m. Titik luluh Tembaga dengan konstanta pegas tersebut adalah 0.0000463538. Artinya, Tembaga mulai mengalami deformasi plastis

9

5.5. HASIL DAN DISKUSI saat Tembaga meregang sejauh 0.0000463538. dilihat pada Gambar 5.5.

Kurva deformasinya dapat

Gambar 5.5: Perbandingan tegangan dan regangan tembaga

Variasi parameter juga dilakukan dengan mengubah nilai konstanta pegas untuk melihat perubahan titik regangan luluhnya. Konstanta pegas yang digunakan secara berurut terdapat pada tabel berikut.

Variasi ke I II III

Tabel 5.4: k1 (N/m) 2.89 × 105 1.89 × 105 0.89 × 105

Variasi parameter

k2 (N/m) 2.79 × 105 1.79 × 105 0.79 × 105

Regangan luluh 0.0000463538 0.000119974 0.000677565

Berdasarkan Tabel 5.4, jelas terlihat bahwa perubahan nilai konstanta pegas, sangat mempengaruhi titik luluh Tembaga. Perubahan variasi parameter nilai konstanta pegas memberikan perubahan terhadap titik dimana bahan mulai berdeformasi. Pada kasus Tembaga, juga terlihat bahwa titik luluh Tembaga semakin lebar saat nilai konstanta pegasnya diperkecil. Pada variasi kedua, bahan mulai berdeformasi saat meregang 0.000119974, 60% lebih panjang dari variasi pertama. Variasi ketiga, dengan nilai konstanta yang lebih kecil, yaitu k1 = 0, 89 × 105 N/m dan k2 = 0, 79 × 105 N/m bahan mulai berdeformasi saat meregang 0.000677565. Nilai regangan ini 93% lebih panjang dari variasi pertama, dan 82% lebih panjang dari variasi kedua. Perbedaan antara 3 variasi tersebut, dapat digambarkan pada Gambar 5.6.


Gambar 5.6: Perbandingan tegangan dan regangan tembaga dengan beberapa variasi konstanta pegas Jika merujuk pada literatur yang telah ada, besar nilai regangan (ǫ) yang dianalisis secara numerik, memiliki nilai yang sama dengan nilai yang dianalisis secara analitik. menggunakan persamaan (5.12). Analisis nilai regangan secara analitik membutuhkan besaran Modulus Young dan Tegangan Luluh. Tabel 5.5 menyajikan parameter Modulus Young (E), Tegangan Luluh (σ), dan hasil analitik regangan (ǫ).

Tabel 5.5: Analisis analitik penentuan regangan Variasi ke σ(N/m2 ) E(N/m2 ) Regangan luluh Aluminium 2.8 × 106 1011 0.00028 6 Tembaga 69 × 10 11 × 1010 0.000672 Berdasarkan Tabel 5.5 Regangan yang diperoleh secara analitik sesuai dengan analisis numerik pada Alumunium saat konstanta pegasnya k1 = 1, 35×105 N/m dan k2 = 1 × 105 N/m, sedangkan pada Tembaga nilai regangan yang diperoleh secara analitik sesuai dengan nilai regangan luluh yang dianalisis secara numerik saat konstanta pegasnya k1 = 0, 89 × 105 N/m dan k2 = 0, 79 × 105 N/m.

5.6

Kesimpulan

Telah diperoleh bahwa adanya perubahan nilai konstanta pegas dan jenis bahan yang berbeda mempengaruhi besarnya nilai regangan luluh suatu material. Pada Aluminium, dengan konstanta pegas k1 = 1, 35 × 105 N/m dan k2 = 1 × 105 N/m, besarnya regangan luluh adalah 0.00200796. Pada Tembaga,

5.7. REFERENSI

11

dengan k1 = 0, 89×105 N/m dan k2 = 0, 79×105 N/m besarnya regangan luluh adalah 0.000677565.

5.7

Referensi

1. Souisa, Metheus. ”Analisis Modulus Elastisitas dan Angka Poison Bahan dengan Uji Tarik”. Jurnal Ilmu Matematika dan Terapan. 2011, 9.14 2. M. Vable ”Mechanics of Materials: Mechanical Properties of Materials”, Chapter 3. 2014. http://www.me.mtu.edu/ mavable/Book/Chapter3.pdf. 3. J.H. Hollomon, Tensile Deformation, Trans. AIME, 162, 268 (1945). 4. Daryus, Asyari. ”Diktat Material Teknik”, Teknik Mesin Universitas Darma Persada, Jakarta. 5. Sparisoma Viridi, ”Catatan Kuliah Komputasi Sistem Fisis”. 2015. Institut Teknologi Bandung.

5.8

Lampiran

5.8.1

Kode C++

#include #include #include #include

using namespace std; double f(double c); double metodeSecant (double cn, double cn_1, double e); //Penentuan Titik Luluh Bahan (Tembaga) dengan Pendekatan Metode Secant int main(int argc, char** argv) { double k1=0.89e5, k2=0.79e5 ; double A=1.6e-5, l = 0.95; double m = (2*k1+k2)/A*(l*k2/(k1+2*k2)); //sigma satu double n=0.54, K=3.15e8; //menulis file output ofstream fout; fout.open("spring.txt"); //definisi nilai awal double x1=0.0, x2=0;

12| 5. ANALISIS TITIK LULUH MATERIAL DENGAN MENGGUNAKAN METODE SECA double eps=0; double sigma=0; double double double double double double

dx1=0.000001; x=0; em = ((2*k1+k2)/A)*(l*k2/(k1+2*k2)); epsIntersect = pow((m/(n*K)), 1/(n-1)); //titik potong c_a = (pow((m/n/K),(1/(n-1)))*(m-m/n)); //nilai C c_n = metodeSecant(25, 80, 0.0001);

cout << "Titik luluh : " << epsIntersect << endl; do{ x1=x1+dx1; x2=x1*(k1+k2)/k2; x=x1+x2; eps=abs(x)/l; //daerah 1 if (eps<=epsIntersect){ sigma=em*eps; //nilai tegangan daerah 1 } //daerah 2 else if(eps>epsIntersect){ sigma=K*pow(eps,n)+c_n; //nilai tegangan daerah dua } fout << eps << "\t"<<sigma<<"\t"<<endl; } while(eps<=0.01); return 0; } double { double double double double return }

f(double c) k1=0.89e5, k2=0.79e5 ; A=1.6e-5, l = 0.95; m = (2*k1+k2)/A*(l*k2/(k1+2*k2)); n=0.54, K=3.15e8; m-(m/n)-c*pow(m/(n*K),(1/(1-n)));

//mencari nilai c dengan metode secant double metodeSecant (double cn, double cn_1, double e) { double d; double error = 1; while (error>= e) { d = (cn_1 - cn) * f(cn_1) / (f(cn_1) - f(cn)); error = abs(f(cn_1)-f(cn)); cn = cn_1; cn_1 = cn - d; } return cn;

5.8. LAMPIRAN } Kode C++ di atas di panggil dengan cara fout.open("spring.txt");

5.8.2 set set set set set set set set

Scrip Gnuplot

term gif size 400,400; output ’spring.jpg’; xrange[1e-4:4.5e-3]; yrange[0.1e6:1.6e7]; xlabel "Regangan" ylabel "Tegangan (N/m2)" title "Grafik Tegangan terhadap Regangan pada Tembaga" key off;

Skrip gnuplot di atas dipanggil dengan cara plot ’./spring.txt’

13


|6

Gelombang Tali Fatimah | [email protected] Nanda Novita | [email protected] Abednego Wiliardy | [email protected] Suatu tali bermassa dengan tegangan tertentu dimodelkan sebagai sejumlah partikel yang saling terhubung dengan pegas. Kondisi khusus dipilih dengan mengambil massa partikel seragam dan konstanta pegas yang juga seragam dengan tujuan untuk memodelkan suatu tali yang homogen. Tali ini terikat di salah satu ujungnya, sedangkan ujung lainnya digetarkan dengan frekuensi tetap. Persamaan gerak ditinjau hanya dengan melibatkan interaksi pegas antarpartikel. Pengaruh gravitasi dan gaya lainnya diabaikan dalam simulasi ini. Pada saat kondisi awal, setiap partikel berada dalam keadaan diam pada posisi terdistribusi merata di sepanjang bidang horizontal. Penyelesaian persamaan diferensial dilakukan secara numerik dengan menggunakan metode Euler. Persamaan diferensial tersebut meliputi persamaan kecepatan dan persamaan posisi setiap partikel. Persamaan diferensial ini didapatkan dengan menyelesaikan persamaan Newton untuk sistem pegas terkopel seperti yang dapat dilihat pada (Halliday, 2008). Dengan menggunakan model ini, diharapkan dapat diperoleh simulasi gerak gelombang berdiri pada tali. Selain itu, akan ditinjau juga kondisi apa saja yang harus dipenuhi agar model ini dapat mensimulasikan gerak tersebut sesuai dengan model teoretik. Model teoretik yang digunakan merupakan model gelombang tali sebagai gelombang analitik yang diperoleh dengan mendapatkan persamaan gelombang seperti yang ditunjukkan oleh (Pain, 2006). Pada bagian pertama laporan ini, akan ditunjukkan penurunan persamaan gelombang pada tali secara analitik. Rumusan inilah yang akan digunakan sebagai acuan model teoretik. Selanjutnya akan dicari persamaan diferensial gerak model tali sebagai partikel-partikel yang berinteraksi dengan pegas. Persamaan diferensial ini akan diselesaikan kemudian secara numerik dengan syarat-syarat batas yang ditentukan. Syarat batas yang digunakan terkait dengan peristiwa gelombang berdiri dengan ujung terikat. Setelah itu, akan ditunjukkan salah 1

2

| 6. GELOMBANG TALI

satu contoh hasil yang didapatkan dari simulasi. Batasan apa saja yang harus dipenuhi model ini ditentukan selanjutnya. Pada akhirnya, akan disimpulkan hasil dari penelitian ini dan semua file yang diperlukan dilampirkan pada bagian akhir dari laporan ini.

6.1

Teori

Pada bagian ini akan dijelaskan teori yang digunakan sebagai ide dasar gelombang berdiri, bagaimana pendekatannya dalam pemodelan yang digunakan, serta rumusan numerik yang akan digunakan dalam pemodelan tersebut.

6.1.1

Persamaan gelombang pada tali

Tinjau suatu segmen infinitesimal dari suatu tali homogen yang bergerak pada arah vertikal. Segmen tersebut akan bergerak dengan memenuhi hukum II Newton tentang gerak di mana gaya-gaya yang bekerja timbul akibat adanya tegangan pada tali. Untuk kasus ini, gaya gravitasi ataupun gaya luar lainnya tidak akan diperhitungkan. Jika kerapatan tali per satuan panjang adalah ρ, maka persamaan Newton tersebut dapat dituliskan sebagai X

F~ = ρ(dS)~a.

(6.1)

di mana dS merupakan elemen panjang dari segmen tali tersebut. Untuk elemen tali yang sangat kecil, dS dapat kita tuliskan sebagai

dS

= =

p

(dx)2 + (dy)2 r ∂y 2 dx 1 + ∂x

(6.2)

Total gaya yang bekerja pada elemen tali pada arah vertikal adalah X

Fy = T (x + dx) sin (θ + dθ) − T (x) sin (θ)

(6.3)

Untuk gelombang tali yang ideal, elemen tali diasumsikan agar bergerak hanya pada arah vertikal. Selain itu, simpangan yang terbentuk pada tali diambil tidak terlalu besar sehingga θ sangat kecil. Karena itu, dapat diambil sin θ ≈ tan θ = ∂y/∂x. Perlu diperhatikan bahwa untuk θ yang kecil, ∂y/∂x akan bernilai kecil, sehingga kuadrat dari nilainya akan jauh lebih kecil dari 1. Akibatnya, pada persamaan 6.2, panjang segmen tali dapat diambil dS ≈ dx. Dengan asumsiasumsi seperti ini, dari persamaan 6.1 dan 6.3 dapat kita peroleh persamaan gerak elemen tali tersebut sebagai berikut.

3

6.1. TEORI

y T(x+dx) +d elemen tali dS

T(x)

x

x + dx

x

Gambar 6.1: Gerak sebuah elemen tali dengan panjang dS dan tegangan tali T (x) dengan sudut θ pada x serta tegangan tali T (x + dx) dengan sudut θ + dθ pada x + dx.

∂y ∂y − T T ∂x x+dx ∂x x ∂ ∂y T ∂x ∂x

= =

∂2y ∂t2 2 ∂ y ρ 2 ∂t

ρdx

(6.4)

Karena simpangan diambil cukup kecil dan tali dianggap tidak mengalami pemelaran yang signifikan, maka tegangan tali T dapat dianggap bernilai tetap untuk sepanjang x, sehingga persamaan 6.4 yang sudah didapatkan sebelumnya dapat ditulis kembali secara lebih sederhana sebagai

T

∂ ∂y ∂x ∂x ∂2y ∂x2

= =

∂2y ∂t2 ρ ∂2y T ∂t2

ρ

(6.5)

Bandingkan persamaan 6.5 ini dengan persamaan gelombang, ∂2y 1 ∂2y = , ∂x2 v 2 ∂t2

(6.6)

di mana v merupakan cepat rambat dari gelombang. Bandingkan persamaan gelombang ini dengan persamaan 6.5 sehingga cepat rambat gelombang pada tali dapat kita peroleh sebagai

4

| 6. GELOMBANG TALI

v=

s

T . ρ

(6.7)

Untuk suatu tali dengan ujung terikat, kondisi resonansi gelombang berdiri akan terjadi ketika

L=n

λ 2

,

(6.8)

di mana L adalah panjang dari tali dan n = 1, 2, 3, . . . menunjukkan orde dari resonansi yang terjadi. Karena itu, dari persamaan 6.7 dan 6.8 dapat diperoleh modus-modus frekuensinya, yaitu

f

= =

v λ s n T . 2L ρ

(6.9)

Ketika osilator dengan frekuensi tersebut digetarkan pada ujung tali lainnya, maka akan terbentuk gelombang berdiri dengan jumlah perut dan simpul yang sesuai dengan orde frekuensinya. Jumlah perut gelombang yang terbentuk akan sama dengan orde resonansinya sedangkan jumlah simpul yang terbentuk akan satu buah lebih banyak dari jumlah perut gelombang yang terbentuk.

6.1.2

Persamaan gerak partikel pada model tali

Ketika tali dimodelkan sebagai partikel-partikel, maka dinamika tali tersebut merupakan dinamika dari keseluruhan partikel-partikel penyusunnya. Untuk suatu partikel ke-i, hukum Newton yang berlaku adalah X

F~i = mi~ai .

(6.10)

Jika sebuah partikel hanya mengalami gaya pegas akibat interaksi dengan dua tetangga terdekatnya, maka X

F~i = −ki,i−1 ∆~xi,i−1 − ki,i+1 ∆~xi,i+1 .

(6.11)

Jika yang dimodelkan adalah seutas tali homogen, maka setiap partikel akan memiliki massa m dan konstanta pegas k yang seragam. Sebagai akibatnya, percepatan yang dialami oleh partikel ke-i dapat diperoleh sebagai

5

6.1. TEORI

y Fi+1 , i+2

partikel ke-i+1 Fi , i+1 Fi , i-1 partikel ke-i

xi

xi+1

x

Gambar 6.2: Suatu tali yang dimodelkan sebagai sejumlah partikel bermassa yang terhubung dengan pegas.

~ai = −

k (∆~xi,i−1 + ∆~xi,i+1 ). m

(6.12)

Jika terdapat N buah partikel yang menyusun tali, maka percepatan setiap partikel dapat ditulis ulang dalam bentuk yang lebih umum sebagai  k  − m (∆~xi,i+1 ), k ~ai = − m (∆~xi,i−1 ),   k − m (∆~xi,i−1 + ∆~xi,i+1 ),

i=1 i=N 1
(6.13)

di mana ∆~xi,j merupakan vektor simpangan pegas oleh partikel ke-i terhadap partikel ke-j yang dapat dirumuskan sebagai berikut.

∆~xi,j = (|~xi − ~xj | − l0 )

~xi − ~xj |~xi − ~xj |

(6.14)

Konstanta l0 dalam persamaan tersebut merupakan panjang pegas dalam kondisi relaksasi. Untuk tali yang homogen, l0 ini bernilai sama untuk seluruh pegas penyusun tali tersebut. Setelah percepatan partikel ~ai didapatkan, maka dapat diperoleh persamaan gerak dari setiap partikel yang memenuhi persamaan diferensial berikut.

d~xi (t) dt d2 ~xi (t) dt2

= ~vi (t) =

d~vi (t) dt

(6.15) =

~ai

(6.16)

6

6.2

| 6. GELOMBANG TALI

Syarat batas

Untuk suatu tali dengan ujung terikat dan ujung lainnya digetarkan secara harmonik, setiap partikel penyusunnya bergerak dengan percepatan yang diberikan pada persamaan 6.13 kecuali untuk partikel yang berada di ujung-ujung tali. Jika ujung tali yang terikat adalah partikel ke-N dan ujung tali yang digetarkan secara harmonik adalah partikel ke-1, maka

~xi (t) ~vi (t)

=

(

=

(

0 eˆx + A sin (ωt) eˆy , −L eˆx + 0 eˆy ,

0 eˆx + Aω cos (ωt) eˆy , ~0,

i=1 i=N i=1 i=N

(6.17) (6.18)

di mana A merupakan amplitudo getaran yang diberikan (amplitudo osilator), ω = 2πf merupakan frekuensi osilator, dan ~0 menunjukkan vektor nol yang bernilai nol untuk setiap komponennya. Frekuensi yang digunakan pada osilator dibatasi hanya untuk frekuensifrekuensi resonansi, sehingga frekuensi osilator yang digunakan akan memenuhi persamaan 6.9 yang dapat ditulis ulang sebagai nπ ω = 2πf = L

s

T , ρ

n = 1, 2, 3, . . .

(6.19)

Kondisi awal dipilih agar setiap partikel terdistribusi secara merata di sepanjang tali dengan kecepatan awal bernilai nol.

~xi (t = 0)

=

−

~vi (t = 0)

= ~0

i−1 Lˆ ex + 0ˆ ey N −1

(6.20) (6.21)

Selain itu, pada mulanya tali sudah diberikan tegangan T tertentu yang cara penerapannya pada model dilakukan dengan mendefinisikan panjang relaksasi pegas sebagai

l0 =

6.3

T L − N −1 k

(6.22)


Model tali sebagai N buah partikel yang dihubungkan dengan N − 1 buah pegas dievaluasi persamaan geraknya dengan menggunakan metode Euler. Setelah

7


mendefinisikan partikel pada kondisi awal yang memenuhi persamaan 6.20 dan 6.21, posisi dan kecepatan partikel setiap saat dapat dihitung dengan menggunakan persamaan 6.17 dan 6.18 untuk partikel yang berada di ujung-ujung tali atau dievaluasi menggunakan metode Euler dengan percepatan yang diberikan pada persamaan 6.12 untuk partikel lainnya. Untuk setiap partikel yang tidak berada di ujung tali, berlaku bahwa

~ai (t + ∆t) ~vi (t + ∆t) ~xi (t + ∆t)

k (∆~xi,i−1 (t) + ∆~xi,i+1 (t)), m = ~vi (t) + ~ai (t + ∆t)∆t, = ~xi (t) + ~vi (t + ∆t)∆t. =

−

(6.23) (6.24) (6.25)

Dapat dilihat bahwa vektor percepatan di setiap waktu ditentukan oleh vektor simpangan pegas di waktu sebelumnya. Vektor simpangan pegas pada persamaan tersebut telah didefinisikan pada persamaan 6.14. Vektor kecepatan pada waktu t + ∆t dihitung dengan metode Euler menggunakan vektor percepatan pada waktu t+∆t. Begitu juga vektor posisi pada waktu t+∆t dihitung dengan metode Euler menggunakan vektor kecepatan pada waktu t + ∆t. Jadi, baik vektor posisi dan kecepatan untuk waktu yang selanjutnya dihitung dengan metode Euler menggunakan nilai turunan pertamanya yang sudah diperbarui. Untuk partikel yang berada di ujung-ujung tali, vektor posisi dan kecepatannya setiap waktu ditentukan secara eksak dengan menggunakan fungsi analitik. Bentuk fungsi ini telah diberikan pada persamaan 6.17 dan 6.18. Posisi dan kecepatan partikel ini ditentukan secara eksak agar syarat batas tali ini dapat tetap dipertahankan. Sebagai catatan, syarat batas yang berlaku adalah ujung tali yang terikat (posisi tetap dan kecepatan nol) dan ujung lainnya yang digetarkan secara harmonik dengan osilator. Seluruh persamaan gerak yang telah diungkapkan pada persamaan 6.23 sampai 6.25 serta persamaan 6.17 dan 6.18 diterapkan dalam program dengan algoritma berikut ini. L01. L02. L03. L04. L05. L06. L07. L08. L09. L10. L11. L12. L13.

Parameter waktu: ∆t, tbeg , tend ? Parameter tali: L, T , N , m, k? Parameter osilator: A, ω? t = tbeg Tave = 0 θave = 0 l0 = (L/(N − 1)) − (T /k) ˆx + 0 eˆy , i = 1 . . . N ~xi (t) = − Ni−1 −1 L e ~vi (t) = 0 eˆx + 0 eˆy , i = 1 . . . N xi+1 (t) ∆~xi,i+1 (t) = (|~xi (t) − ~xi+1 (t)| − l0 ) |~~xxii (t)−~ (t)−~ xi+1 (t)| , ~ xi (t)−~ (t) ∆~xi,i−1 (t) = (|~xi (t) − ~xi−1 (t)| − l0 ) |~xi (t)−~xxi−1 , i−1 (t)| PN k 1 Tave = Tave + N −1 i=2 m |∆~xi,i−1 (t)| PN ∆~ x (t)·ˆ ex θave = θave + N 1−1 i=2 acos( |∆~i,i−1 xi,i−1 (t)| )

k (∆~xi,i−1 (t) + ∆~xi,i+1 (t)), L14. ~ai (t + ∆t) = − m

i = 2 . . . (N − 1)

i = 2 . . . (N − 1)

i = 2 . . . (N − 1)

8 L15. L16. L17. L18. L19. L20. L21. L22. L23.

| 6. GELOMBANG TALI ~vi (t + ∆t) = ~vi (t) + ~ai (t + ∆t)∆t, i = 2 . . . (N − 1) ~xi (t + ∆t) = ~xi (t) + ~vi (t + ∆t)∆t, i = 2 . . . (N − 1) ~v1 (t + ∆t) = 0 eˆx + Aω cos (ω(t + ∆t)) eˆy ~x1 (t + ∆t) = 0 eˆx + A sin (ω(t + ∆t)) eˆy ~vN (t + ∆t) = 0 eˆx + 0 eˆy ~xN (t + ∆t) = −L eˆx + 0 eˆy t = t + ∆t t ≤ tend → L10. Tave = tend∆t −tbeg Tave

L24. θave =

6.4

∆t tend −tbeg θave

Hasil dan diskusi

6.4.1

Simulasi Model Gelombang Tali

Untuk menguji model tali yang sudah dibuat, dilakukan simulasi gelombang berdiri pada tali dengan orde 2 di mana parameter masukan simulasi yang digunakan ditunjukkan pada tabel berikut.

Tabel 6.1: Parameter masukan simulasi gelombang berdiri orde 2. Parameter ∆t tbeg tend L T N m k A ω

Nilai 10−5 s 0s 10 s 1m 10 N 100 1 gram 105 N/m 0.5 mm 20π rad/s

Berikut ini pada gambar 6.3 adalah hasil simulasi gelombang tali pada saat tali mulai digetarkan untuk pertama kali dari keadaan diam. Pada gambar 6.3 tersebut ditunjukkan dinamika partikel-partikel penyusun tali pada saat-saat awal tali digetarkan. Dapat dilihat bahwa gelombang mulai menjalar dari arah kanan ke kiri dengan amplitudo yang sama dengan amplitudo osilator A. Berdasarkan gambar hasil simulasi tersebut, waktu yang diperlukan oleh gelombang untuk bergerak dari ujung tali yang satu ke ujung lainnya adalah Ttempuh = 0.1 s. Jika panjang tali L = 1 m, maka cepat rambat gelombang pada tali adalah

9


Gambar 6.3: Hasil simulasi gelombang pada tali pada saat pertama kali digetarkan dengan nilai parameter sesuai dengan tabel 6.1: (kiri-atas) saat t = 0.025 s, (kanan-atas) saat t = 0.05 s, (kiri-bawah) saat t = 0.075 s, dan (kanan-bawah) saat t = 0.1 s. Sumbu x dan y menunjukkan koordinat ruang dengan satuan meter dan skala yang disesuaikan.

v

= =

L Ttempuh 10 m/s.

=

1m 0.1 s (6.26)

Hasil ini bersesuaian dengan hasil yang didapatkan dengan menggunakan persamaan 6.7, yaitu

v

=

s

T , ρ

=

s

=

10 N (100 × 0.001 kg)/(1 m)

10 m/s.

ρ=

N ×m L

(6.27)

Melihat kesesuaian ini, kita bisa harapkan bahwa simulasi model tali ini dapat menggambarkan model teoretik dengan cukup baik. Meskipun begitu, masih ada aspek-aspek lainnya yang perlu diperhatikan yang akan dibahas selanjutnya. Selanjutnya akan ditunjukkan hasil simulasi gelombang tali pada saat gelombang pantul muncul untuk pertama kali dan mulai berinterferensi dengan gelombang datang seperti yang ditunjukkan pada gambar 6.4. Keadaan ini mulai terjadi sejak detik ke-0.1 s. Pada gambar 6.4 tersebut dapat terlihat saat gelombang pantul muncul untuk pertama kali. Gelombang pantul muncul dari ujung sebelah kiri dan mulai

10

| 6. GELOMBANG TALI

Gambar 6.4: Hasil simulasi gelombang pada tali pada saat gelombang pantul muncul untuk pertama kali dengan nilai parameter sesuai dengan tabel 6.1: (kiri-atas) saat t = 0.125 s, (kanan-atas) saat t = 0.15 s, (kiri-bawah) saat t = 0.175 s, dan (kanan-bawah) saat t = 0.2 s. Sumbu x dan y menunjukkan koordinat ruang dengan satuan meter dan skala yang disesuaikan. berinterferensi dengan gelombang datang. Pada waktu t = 0.125 s, interferensi destruktif terjadi untuk 1/4 bagian tali di sebelah kiri. Hal ini dikarenakan gelombang pantul baru muncul untuk pertama kalinya pada t = 0.1 s dan baru menempuh jarak 25 cm pada t = 0.125 s. Pada t = 0.15 s interferensi yang terjadi bersifat konstruktif dan pada t = 0.175 s terjadi lagi interferensi yang bersifat destruktif yang juga masih terjadi secara sebagian. Gelombang datang dan gelombang pantul mulai berinterferensi secara penuh pada saat t = 0.2 s ketika gelombang pantul mencapai ujung yang satunya. Sejak saat inilah interferensi yang terjadi mulai melahirkan fenomena gelombang berdiri. Dengan mengamati interferensi konstruktif dan destruktif yang terjadi secara terus menerus, kita akan melihat bentuk-bentuk dari perut dan simpul dari gelombang berdiri. Untuk kondisi dengan parameter sesuai dengan tabel 6.1, dapat kita lihat bahwa jumlah perut gelombang berdiri yang terbentuk adalah dua buah, sesuai dengan yang diharapkan secara teoretik. Jika simulasi ini diperhatikan baik-baik, kondisi interferensi konstruktif dan destruktif pada gelombang terjadi secara periodik dengan frekuensi tertentu. Berdasarkan hasil pengamatan, waktu yang diperlukan antara interferensi konstruktif dan destruktif yang berutuan adalah 0.025 s, atau dengan kata lain, waktu antara dua interferensi konstruktif yang terjadi secara berurutan adalah Tinter = 0.05 s. Karena itu, dapat diperoleh frekuensi terbentuknya perut dan simpul gelombang sebagai

ωsw =

2π rad = 40π rad/s = 2ω. Tinter

(6.28)

Sebagai catatan, gelombang yang dimaksud dalam bahasan ini bukanlah gelombang analitik yang sesungguhnya, melainkan merupakan dinamika sejumlah partikel yang bergerak secara bersama-sama membentuk ”gelombang”.

11


Perhatikan gambar 6.3 pada waktu t = 0.1 s dan gambar 6.4 pada waktu t = 0.2 s. Dapat dilihat bahwa pada t = 0.2 s amplitudo gelombang yang terbentuk dua kali lebih besar dibandingkan amplitudo gelombang pada saat t = 0.1 s. Secara analitik hal ini terjadi karena gelombang datang dan gelombang pantul berinterferensi secara konstruktif sehingga amplitudo total menjadi dua kalinya. Akan tetapi, perlu kita ingat bahwa gambar yang ditunjukkan tersebut bukanlah gelombang yang sesungguhnya, melainkan merupakan ”gelombang” partikel, dan simulasi gerak ”gelombang” partikel ini memberikan hasil yang cukup baik dan sesuai dengan hasil yang seharusnya diberikan oleh gelombang analitik. Pada waktu t = 0.2 s, amplitudo gelombang menjadi dua kalinya karena adanya interferensi dari dua buah gelombang. Mengingat getaran dari osilator diberikan secara terus menerus, maka sudah sewajarnya amplitudo gelombang yang terbentuk akan semakin besar. Amplitudo ini akan terus bertambah besar sampai batas tertentu. Hal ini dikarenakan semakin besar amplitudo maka akan semakin besar pula gaya pegas yang melawan simpangannya.

6.4.2

Batasan Model

Berdasarkan rumusan teoretik yang telah diberikan, terdapat beberapa anggapan yang digunakan dalam perumusannya. Batasan model yang dimaksud di sini adalah suatu keadaan di mana model yang disimulasikan mendekati model teoretik tersebut. Model teoretik ini adalah model yang memenuhi persamaan diferensial gelombang seperti yang diungkapkan pada persamaan 6.5. Anggapan yang diambil dalam rumusan teoretik ini adalah 1. simpangan pada tali jauh lebih kecil dibandingkan panjang tali, 2. sudut antara elemen tali dengan bidang horizontalnya sangat kecil, dan 3. tegangan tali dianggap bernilai tetap di sepanjang tali. Poin kedua dan ketiga dari anggapan yang diambil tersebut sesungguhnya merupakan konsekuensi dari diambilnya anggapan yang pertama. Jadi, keaadaan pertama dapat kita katakan terpenuhi apabila keadaan kedua dan ketiga terpenuhi. Sudut antara elemen tali dengan bidang horizontalnya dapat dikatakan sangat kecil apabila berlaku

sin(θ) ≈ tan(θ) =

∂y . ∂x

(6.29)

Secara sederhana, kita bisa anggap persamaan ini berlaku apabila kesalahan yang diberikan tidak lebih dari 1%, sehingga

∀ x,

error =

| tan(θ(x)) − sin(θ(x))| ≤ 0.01, | sin(θ(x))|

−L ≤ x ≤ 0.

(6.30)

12

| 6. GELOMBANG TALI

Kondisi ini dapat terpenuhi ketika ∀ x,

−0.140836 rad ≤ θ(x) ≤ 0.140836 rad,

−L ≤ x ≤ 0.

(6.31)

Tegangan tali dapat dikatakan bernilai tetap di sepanjang tali apabila perubahan tegangan tali yang terjadi tidak lebih dari 1% besar tegangan tali yang diberikan mula-mula, sehingga

∀ x,

error =

|T (x) − Tinitial | ≤ 0.01, Tinitial

−L ≤ x ≤ 0

(6.32)

atau dengan kata lain, tegangan tali tidak dapat dikatakan bernilai tetap apabila

∃ x,

|T (x) − Tinitial | > 0.01, Tinitial

−L ≤ x ≤ 0

(6.33)

Jika persamaan 6.31 dan 6.33 terpenuhi, maka kita bisa katakan bahwa batasan model ini terpenuhi secara ketat. Kondisi ini dapat dianalisa dengan mengamati nilai tegangan maksimum dan sudut maksimum yang terbentuk. Jika tegangan atau sudut maksimum memenuhi kondisi di atas, maka batasan model tersebut tentu terpenuhi secara ketat. Akan tetapi, kondisi terpenuhi secara ketat itu sulit untuk dipenuhi karena cakupannya yang sangat sempit. Karena itu, batasan model ini dapat dikatakan terpenuhi secara longgar apabila besaran yang ditinjau cukup nilai rata-ratanya saja, yaitu

|θ|ave |Tave − Tinitial | Tinitial

≤

0.140836 rad

(6.34)

≤

0.01

(6.35)

Selanjutnya, akan dilakukan variasi beberapa parameter masukan untuk menentukan batasan apa saja yang harus dipenuhi agar kondisi-kodisi di atas tercapai. Parameter yang akan divariasikan adalah amplitudo osilator dan frekuensi osilator. Pada tabel 6.2 ditampilkan nilai parameter yang digunakan. Amplitudo osilator akan divariasikan dari 0.5 mm sampai 10.0 mm untuk empat buah frekuensi yang berbeda yang terkait dengan gelombang berdiri orde 1, 2, 4, dan 8. Jika parameter yang digunakan sesuai dengan tabel 6.2, maka tegangan tali kritis untuk kasus ini adalah Tcrit = 1.01 Tinitial = 10.1 N.

(6.36)

Berikut ini adalah kurva hasil pengamatan parameter tegangan tali dan sudut terhadap perubahan parameter amplitudo osilator untuk beberapa buah frekuensi yang berbeda.

13


Tabel 6.2: Parameter masukan untuk variasi amplitudo dan frekuensi osilator. Parameter ∆t tbeg tend L T N m k

Nilai 10−5 s 0s 5s 1m 10 N 100 1 gram 105 N/m

1.6 ambang orde 1 orde 2 orde 4 orde 8

Sudut Rata-Rata (X 0.140836 rad)

1.4

1.2

1

0.8

0.6

0.4

0.2

0 0

0.002

0.004

0.006

0.008

0.01

Amplitudo Osilator (m)


Pada gambar 6.5, simulasi memenuhi batasan model secara longgar untuk orde 1, 2, dan 4 pada semua rentang amplitudo osilator yang digunakan. Simulasi untuk orde 8 menghasilkan sudut yang lebih besar dibandingkan ambang yang diberikan untuk amplitudo A > 5 mm. Berdasarkan gambar 6.6, batasan model secara ketat tidak dapat dicapai untuk amplitudo A ≥ 1.5 mm bagi semua orde yang disimulasikan. Untuk parameter tegangan, hasil diberikan pada gambar 6.7 dan 6.8. Pada kedua gambar tersebut dapat terlihat bahwa batasan model dapat terpenuhi secara longgar maupun ketat untuk amplitudo A ≤ 1 mm untuk semua orde. Selain itu, dari keempat gambar tersebut, terjadi lonjakan yang cukup besar ketika amplitudo sudah melewati angka 5 mm untuk semua orde. Karena itu, untuk amplitudo osilator yang lebih besar dari itu kemungkinan akan menghasilkan bentuk model yang kurang baik. Jadi, dengan merangkum semua hasil yang didapatkan, bisa dikatakan bahwa model gelombang tali ini dengan parameter sesuai dengan tabel 6.2 memenuhi batasan model secara teoretik untuk amplitudo osilator A ≤ 1 mm.

14

| 6. GELOMBANG TALI 10 ambang orde 1 orde 2 orde 4 orde 8

9

Sudut Maksimum (X 0.140836 rad)

8 7 6 5 4 3 2 1 0 0

0.002

0.004

0.006

0.008

0.01



6.5

Ringkasan

Dari simulasi yang telah dilakukan, dapat diamati terbentuknya fenomena gelombang berdiri dengan menggunakan model tali sebagai partikel-partikel. Telah diperoleh bahwa cepat rambat gelombang pada tali dari hasil simulasi model sesuai dengan prediksi yang diberikan oleh teori. Jumlah perut dan simpul gelombang yang terbentuk dari hasil simulasi yang didapatkan juga sesuai dengan prediksi kondisi resonansi secara teoretik. Selain itu, telah diperoleh juga bahwa frekuensi terbentuknya perut-simpul gelombang berdiri adalah dua kali frekuensi yang diberikan oleh osilator. Telah diamati juga bahwa amplitudo gelombang berdiri ini terus bertambah besar hingga mencapai nilai tertentu, kemudian kembali mengecil menuju nol. Setelah itu, amplitudo ini kembali bertambah besar. Proses ini terjadi terus menerus secara periodik yang menghasilkan fenomena kembang kempis perut gelombang berdiri. Meskipun model ini cukup baik untuk dapat mensimulasikan model teoretik, masih terdapat batasan-batasan yang harus dipenuhi agar model ini dapat memberikan hasil yang baik. Hal ini terkait dengan anggapan-anggapan yang diambil dalam perumusan model gelombang tali secara teoretik. Batasan model ini melingkupi batasan secara ketat dan batasan secara longgar. Batasan-batasan ini diatur dalam persamaan 6.30 sampai 6.35 yang mensyaratkan agar kesalahan model tidak melebihi 1% dari nilai yang diharapkan. Untuk model dengan parameter sesuai dengan tabel 6.2, telah dilakukan simulasi dan berdasarkan gambar 6.5 sampai 6.8 didapatkan bahwa batasan secara ketat dapat tercapai pada orde gelombang berdiri 1, 2, 4, dan 8 untuk amplitudo osilator A ≤ 1 mm. Hal ini menunjukkan bahwa simulasi peristiwa gelombang berdiri sebaiknya menggunakan amplitudo osilator yang sangat kecil agar model dapat menggambarkan model teoretik dengan baik.

15

6.6. PENELITIAN LANJUTAN 6 ambang orde 1 orde 2 orde 4 orde 8

Tegangan Rata-Rata (X 10 N)

5

4

3

2

1

0 0

0.002

0.004

0.006

0.008

0.01



6.6

Penelitian Lanjutan

Masih banyak hal yang perlu ditinjau dalam model tali partikel ini. Dalam penelitian ini, parameter yang ditinjau untuk menentukan batasan model hanya parameter osilator, yaitu amplitudo dan frekuensi. Parameter tali yang digunakan masih belum diperhitungkan dalam menentukan batasan model tersebut. Parameter ini meliputi jumlah partikel, rapat massa tali, dan konstanta pegas. Selain itu, masih terdapat juga parameter keluaran yang masih belum ditinjau, seperti energi kinetik, energi potensial, dan besaran fisis lainnya baik secara total ataupun untuk setiap partikel pada setiap waktu. Karena itu, penelitian selanjutnya dapat dilakukan untuk menganalisa nilai parameter-parameter ini. Penelitian lanjutan juga dapat dilakukan untuk memperbaiki model tali partikel ini. Perbaikan dapat dilakukan dengan mengubah metode numerik yang digunakan untuk menyelesaikan persamaan diferensial. Hal lainnya yang dapat dilakukan yaitu dengan memperkecil time step numerik yang digunakan. Selain itu, perbaikan lainnya juga dapat dilakukan untuk membuat sistem model ini semakin real, yaitu dengan menambah pengaruh redaman pada pegas, pengaruh medan luar seperti gravitasi, pengaruh viskositas medium, dan sebagainya. Pengembangan lebih jauh dapat dilakukan dengan menghubungkan beberapa tali secara paralel membentuk suatu bidang. Hal ini memungkinkan untuk dilakukannya penelitian lebih lanjut mengenai gelombang pada bidang. Gelombang berdiri hanya merupakan salah satu contoh fenomena yang dapat dijumpai pada kasus gelombang tali. Penelitian yang lebih jauh dapat dilakukan dengan menggunakan model tali partikel ini untuk mengamati fenomena lainnya. Fenomena lainnya yang dapat diamati adalah peristiwa gelombang berdiri dengan ujung terikat bebas. Sebagai catatan, pada penelitian ini ujung tali yang dimodelkan terikat mati. Fenomena yang lebih berbeda dapat diamati jika

16

| 6. GELOMBANG TALI 40 ambang orde 1 orde 2 orde 4 orde 8

35

Tegangan Maksimum (X 10 N)

30

25

20

15

10

5

0 0

0.002

0.004

0.006

0.008

0.01


Gambar 6.8: Gerak sebuah elemen tali dengan panjang dS dan tegangan tali T (x) dengan sudut θ pada x serta tegangan tali T (x + dx) dengan sudut θ + dθ pada x + dx. massa partikel ataupun konstanta pegas diambil tidak seragam. Hal ini dapat dilakukan untuk mengamati peristiwa transmitansi dan reflektansi gelombang. Lebih jauh lagi, sumber getaran pada tali juga dapat divariasikan, misalnya dengan menerapkan osilator dengan frekuensi yang berubah secara bertahap selama waktu tertentu. Selain itu, dapat juga digunakan dua buah osilator dengan frekuensi yang berlainan. Letak sumber getaran juga dapat divariasikan di berbagai titik pada tali. Bahkan, dapat juga diamati apabila sumber getaran memiliki frekuensi dengan spektrum tertentu.

6.7

Referensi

1. David Halliday, Robert Resnick, and Jearl Walker, Fundamentals of Physics, John Wiley & Sons (Asia), 8th, Extended, Student, Edition, (2008). 2. H.J. Pain, The Physics of Vibrations and Waves, John Wiley & Sons (Asia), 6th Edition, (2005).

6.8

Lampiran

6.8.1

Kode C++

#include <sys/types.h> #include <sys/stat.h> #include

6.8. LAMPIRAN #include #include #include #include #include #include

17

<sstream> "rwparams.h" "vect2.h"

using namespace std;

// Define oscillator position function and oscillator velocity function vect2 r_oscillator(double t, double A, double w) {vect2 r = A * sin(w * t) * vect2(0, 1) vect2 v_oscillator(double t, double A, double w) {vect2 v = A * w * cos(w * t) * vect2(0 int main(int argc, char *argv[]) { // Get compiled name const char *arg0 = argv[0]; const char *pname = strstr(arg0, "./"); if(pname != NULL) pname += 2; // Check number of argument if(argc < 2) { cout << "Usage: " << pname << " [pfile]"; cout << endl; cout << "pfile\tparameters file, e.g. params.txt"; cout << endl; exit(1); } // Get filename of pfile const char *pfn = argv[1]; // Read parameters from pfile double m = 0.0; readparam(pfn, "MASS", m); double l = 0.0; readparam(pfn, "LENGTH", l); int N = 0; readparam(pfn, "NPAR", N); double k = 0.0; readparam(pfn, "KSPRING", k); double A = 0.0; readparam(pfn, "AMPLIT", A); double W = 0.0; readparam(pfn, "ANGFRE", W); int nA = 0; readparam(pfn, "PARAPP", nA); double Ti = 0.0; readparam(pfn, "TENSION", Ti); int e0fx = 0; readparam(pfn, "E0FIXX", e0fx);

18

| 6. GELOMBANG TALI

int e0fy = 0; readparam(pfn, "E0FIXY", e0fy); int e1fx = 0; readparam(pfn, "E1FIXX", e1fx); int e1fy = 0; readparam(pfn, "E1FIXY", e1fy); int afx = 0; readparam(pfn, "ALFIXX", afx); int afy = 0; readparam(pfn, "ALFIXY", afy); double tbeg = 0.0; readparam(pfn, "TBEG", tbeg); double tend = 0.0; readparam(pfn, "TEND", tend); double dt = 0.0; readparam(pfn, "TSTEP", dt); double Tdata = 0.0; readparam(pfn, "TDATA", Tdata); double x_scale = 1.0; readparam(pfn, "XSCALE", x_scale); double y_scale = 1.0; readparam(pfn, "YSCALE", y_scale); int creani = 1.0; readparam(pfn, "CREANI", creani); string ofnpre = ""; readparam(pfn, "OUTPRE", ofnpre); int digit = 0; readparam(pfn, "DIGIT", digit); string odir = ""; readparam(pfn, "OUTDIR", odir); // Check whether OUTDIR exists struct stat info; const char *odirc = odir.c_str(); bool OUTDIR_NOT_EXIST = stat(odirc, &info); if(OUTDIR_NOT_EXIST) { cout << endl << "Error: "; cout << "Output directory does not exist: " << odirc << endl; cout << "Suggestion: Create the directory manually "; cout << "and run " << pname << " again" << endl; exit(2); } // Calculate granular parameters double mi = m / N; double li = l / (N - 1); // Define initial position and velocity vect2 r[N], v[N]; r[0].x = nA * li;

6.8. LAMPIRAN

19

for(int i = 0; i < N; i++) { double yi = 0; double xi = r[0].x - (i * li); r[i] = vect2(xi, yi); v[i] = vect2(0, 0); } // Define fixed condition vect2 fixed_condition[N]; fixed_condition[0] = vect2(1 - e0fx, 1 - e0fy) * vect2(1 - afx, 1 - afy); fixed_condition[N - 1] = vect2(1 - e1fx, 1 - e1fy) * vect2(1 - afx, 1 - afy); for(int i = 1; i < N - 1; i++) { fixed_condition[i] = vect2(1 - afx, 1 - afy); } // Define normal length double l0[N]; for(int i = 0; i= Ndata) { // Create sequence number for output file stringstream ss; ss << setfill(’0’); ss << setw(digit);

20

| 6. GELOMBANG TALI

ss << j; j++; // Calculating progress int progress = round(100 * (t - tbeg) / (tend - tbeg)); // Construct output filename string ofnt = ofnpre + "-" + ss.str() + ".txt"; string ofns = odir + "/" + ofnt; cout << "Computing physical value " << progress << "% completed cout.flush();

\r";

// Open output file const char *ofn = ofns.c_str(); ofstream fout; fout.open(ofn); // Write header in output file fout << "# t = " << t << endl; fout << "# x\ty\tvx\tvy" << endl; // Write particles position to output file for(int i = 0; i < N; i++) { fout << r[i].x << "\t"; fout << r[i].y << "\t"; fout << v[i].x << "\t"; fout << v[i].y << endl; // change idata back to zero idata = 0; } // Close output file fout.close(); // Open output animation file std::ofstream afout(oafns.c_str(), std::ios::app); // Write plotting syntax in output animation file afout << "set title \"t = " << t << "\"" << endl; afout << "plot \"" << ofnt.c_str() << "\" notitle" << endl; // Close output animation file afout.close(); } // Define dummy variable to measure average tension and average angle double tenavei = 0.0; double angavei = 0.0;

6.8. LAMPIRAN

21

// Calculate particles acceleration vect2 a[N]; for(int i = 0; i < N; i++) { vect2 Fim = vect2(0,0); vect2 Fip = vect2(0,0); // Calculate spring interaction between particle i and i - 1 if(i != 0) { vect2 rim = r[i] - r[i-1]; vect2 nim = rim >> 1; double lim = length(rim); Fim = -k * (lim - l0[i]) * nim; } // Calculate spring interaction between particle i and i + 1 if(i != N - 1) { vect2 rip = r[i] - r[i+1]; vect2 nip = rip >> 1; double lip = length(rip); Fip = -k * (lip - l0[i]) * nip; } a[i] = (Fim + Fip) / mi; // Update max tension value of spring if(tenmax < length(Fim)) {tenmax = length(Fim);} // Calculate sum of tension (to be averaged) tenavei += length(Fim); } // Average tension for a certain time tenavei = tenavei / (N - 1); tenave += tenavei; // Increase idata and t idata++; t += dt; // Update particles velocity and position + look for peak condition for(int i = 0; i < N; i++) { if(i == nA) { v[i] = v_oscillator(t, A, W) * fixed_condition[i]; r[i] = r_oscillator(t, A, W) * fixed_condition[i]; } else { v[i] = v[i] + (a[i] * dt) * fixed_condition[i]; r[i] = r[i] + (v[i] * dt) * fixed_condition[i]; }

22

| 6. GELOMBANG TALI

// checking if asw (amplitude of standing wave) is smaller than particle devia if(asw < abs(r[i].y)) {asw = abs(r[i].y);}

// checking if angmax (maximum angular) is smaller than particle position angl double angi = 0.0; if(i > 0) {angi = angle(r[i] - r[i-1], (r[i].x - r[i-1].x) * vect2(1, 0));} if(angmax < angi) {angmax = angi;} // Calculate sum of angle (to be averaged) angavei += angi; } // Average angle for a certain time angavei = angavei / (N - 1); angave += angavei; } // Calculate average tension & angle all time tenave = tenave / (((tend - tbeg)/dt) + 1); angave = angave / (((tend - tbeg)/dt) + 1); // Open output parameter file const char *opfn = opfns.c_str(); ofstream fout; fout.open(opfn); // Writing ouput parameters into file fout << "# Standing wave amplitude" << endl; fout << "OSCILLATOR\t" << A << "\tAMPLITUDE" << endl; fout << "STANDING-WAVE\t" << asw << "\tAMPLITUDE" << endl; fout << endl; fout << "# Tension data" << endl; fout <<"APPLIED\t" << Ti << "\tTENSION" << endl; fout << "MAXIMUM\t" << tenmax << "\tTENSION" << endl; fout << "AVERAGE\t" << tenave << "\tTENSION" << endl; fout << endl; fout << "# Angle data" << endl; fout << "MAXIMUM\t" << angmax << "\tANGLE" << endl; fout << "AVERAGE\t" << angave << "\tANGLE" << endl; fout << endl; // Close output parameters file fout.close(); // Display output parameters to user cout << endl << endl; cout << "Oscillator Amplitude = " << A << " m" << endl; cout << "Standing Wave Amplitude = " << asw << " m" << endl; cout << "Amplitude Ratio = " << asw/A << endl;

23

6.8. LAMPIRAN cout cout cout cout cout cout cout cout cout cout

<< << << << << << << << << <<

endl; "Initial "Maximum "Average "Maximum "Average endl; "Maximum "Average endl;

Tension Applied Tension Tension Tension Ratio Tension Ratio

= = = = =

" " " " "

<< << << << <<

Ti << " N" << endl; tenmax << " N" << endl; tenave << " N" << endl; tenmax/Ti << endl; tenave/Ti << endl;

Angle Angle

= " << angmax << " rad" << endl; = " << angave << " rad" << endl;

// Creating animation.gif file if(creani == 1) { // Open the pipe FILE *pipe = popen("gnuplot -persist", "w"); fprintf(pipe, "\n"); // Define x range & y range plot area std::ostringstream x_0, x_1, y_0, y_1; x_0 << - x_scale * l; x_1 << x_scale * l; y_0 << - y_scale * asw; y_1 << y_scale * asw; string x_range = "[" + x_0.str() + ":" + x_1.str() + "]"; string y_range = "[" + y_0.str() + ":" + y_1.str() + "]"; // Converting fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe,

output data into gif animation in pipe "cd \’%s\’\n", odir.c_str()); "reset\n"); "set term gif animate\n"); "set output \"animation.gif\"\n"); "set xrange %s\n", x_range.c_str()); "set yrange %s\n", y_range.c_str()); "set xlabel \"x\"\n"); "set ylabel \"y\"\n"); "load \"%s\"\n", oafnt.c_str()); "set output\n"); "\n");

// Close the pipe fclose(pipe); } return 0; }

Kode C++ di atas dikompilasi menggunakan g++ dari MinGW melalui Command Prompt dengan cara

24

| 6. GELOMBANG TALI

g++ tali.cpp -o tali dan dipanggil melalui Command Prompt dengan cara tali params.txt di mana file ’params.txt’ merupakan file yang berisi seluruh parameter masukan program. Pada saat kompilasi, pastikan file ’rwparams.h’ dan ’vect2.h’ berada dalam direktori yang sama dengan file ’tali.cpp’. Pada saat program dipanggil, pastikan folder ’data’, file ’params.txt’ dan ’out-params.txt’ berada dalam direktori yang sama dengan ’tali.exe’. Untuk mencegah duplikasi data, kosongkan folder ’data’ setiap kali pemanggilan program.

6.8.2

File Parameter

Berikut ini adalah isi dari file parameter masukan ’params.txt’: # Particle parameters MASS 0.1 VARIABLE NPAR 100 VARIABLE LENGTH 1 VARIABLE # Interaction parameters KSPRING 1E5 VARIABLE TENSION 10 VARIABLE # Oscillator parameters AMPLIT 0.0005 VARIABLE ANGFRE 62.831853075 VARIABLE PARAPP 0 STATIC # Fix point parameters E0FIXX 0 STATIC E0FIXY 0 STATIC E1FIXX 1 STATIC E1FIXY 1 STATIC ALFIXX 0 STATIC ALFIXY 0 STATIC # Simulation parameters TBEG 0 STATIC TEND 10 STATIC TSTEP 1E-5 STATIC TDATA 0.05 STATIC XSCALE 1.0 STATIC YSCALE 1.1 STATIC

6.8. LAMPIRAN

25

CREANI 1 STATIC # Output file parameters OUTDIR data STATIC OUTPRE out STATIC DIGIT 4 STATIC Berikut ini adalah isi dari file parameter keluaran ’out-params.txt’: # Standing wave amplitude OSCILLATOR 0.0005 AMPLITUDE STANDING-WAVE 0.0089133 AMPLITUDE # Tension data APPLIED 10 TENSION MAXIMUM 10.7414 TENSION AVERAGE 10.1862 TENSION # Angle data MAXIMUM 0.0549925 ANGLE AVERAGE 0.0139773 ANGLE

6.8.3

File Library C++ Tambahan

Berikut ini adalah isi dari library ’rwparams.h’ yang digunakan untuk baca tulis parameter: #include #include <sstream> #include #ifndef rwparams_h #define rwparams_h using namespace std; void void void void void void

readparam(const char*, const char*, int&); readparam(const char*, const char*, double&); writeparam(const char*, const char*, int); writeparam(const char*, const char*, double); writeparam(const char*, const char*, int, double*); writeparam(const char*, const char*);

void readparam(const char*, const char*, string&); void readparam(const char*, const char*, int&, string&);

26

| 6. GELOMBANG TALI

void readparam(const char*, const char*, double&, string&); void readparam(const char*, const char*, string&, string&); void void void void void

lookparam(const findparam(const findparam(const findparam(const findparam(const

char*); char*, int, char*, int, char*, int, char*, int,

string&); string&, int&, string&); string&, double&, string&); string&, string&, string&);

int countparam(const char*); int countparam(const char*, string); void readparam(const char *pf, const char *pttrn, int &x) { ifstream fin; fin.open(pf); while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { stringstream ss; string col1; string col2; ss << line; ss >> col1; ss >> col2; if(col1 == pttrn) { stringstream ss; ss << col2; ss >> x; break; } } } fin.close(); } void readparam(const char *pf, const char *pttrn, double &x) { ifstream fin; fin.open(pf); while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos;

6.8. LAMPIRAN if(!COMMENT && !BLANK && !SINGLE_VALUE) { stringstream ss; string col1; string col2; ss << line; ss >> col1; ss >> col2; if(col1 == pttrn) { stringstream ss; ss << col2; ss >> x; break; } } } fin.close(); } void readparam(const char *pf, const char *pttrn, string &x) { ifstream fin; fin.open(pf); while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { stringstream ss; string col1; string col2; ss << line; ss >> col1; ss >> col2; if(col1 == pttrn) { stringstream ss; ss << col2; ss >> x; break; } } } fin.close(); } void readparam(const char *pf, const char *pttrn, int &x, string &attrib) { ifstream fin; fin.open(pf);

27

28

| 6. GELOMBANG TALI

while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { stringstream ss; string col1; string col2; string col3; ss << line; ss >> col1; ss >> col2; ss >> col3; if(col1 == pttrn) { stringstream ss; ss << col2; ss >> x; ss << col3; ss >> attrib; break; } } } fin.close(); } void readparam(const char *pf, const char *pttrn, double &x, string &attrib) { ifstream fin; fin.open(pf); while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { stringstream ss; string col1; string col2; string col3; ss << line; ss >> col1; ss >> col2; ss >> col3; if(col1 == pttrn) { stringstream ss; ss << col2;

6.8. LAMPIRAN ss >> x; ss << col3; ss >> attrib; break; } } } fin.close(); } void readparam(const char *pf, const char *pttrn, string &x, string &attrib) { ifstream fin; fin.open(pf); while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { stringstream ss; string col1; string col2; string col3; ss << line; ss >> col1; ss >> col2; ss >> col3; if(col1 == pttrn) { stringstream ss; ss << col2; ss >> x; ss << col3; ss >> attrib; break; } } } fin.close(); } void lookparam(const char *pf) { ifstream fin; fin.open(pf); int number = 0; while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos;

29

30

| 6. GELOMBANG TALI

bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { stringstream ss; string col1; string col2; ss << line; ss >> col1; ss >> col2; number ++; cout << number << ".\t" << col1; cout << "\t= " << col2 << endl; } } cout << endl; if(number == 0) {cout << "There is no parameter" << endl;} if(number == 1) {cout << "Only 1 parameter exists" << endl;} if(number > 1) {cout << "There are " << number << " parameters exist" << endl; fin.close(); } void findparam(const char *pf, int parnum, string &pttrn) { ifstream fin; fin.open(pf); int number = 0; while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { number ++; if(number == parnum) { stringstream ss; ss << line; ss >> pttrn; } } } fin.close(); }

void findparam(const char *pf, int parnum, string &pttrn, int &x, string &attr ifstream fin; fin.open(pf); int number = 0; while(fin.good()) { string line; getline(fin, line);

6.8. LAMPIRAN

31

bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { number ++; if(number == parnum) { stringstream ss; ss << line; ss >> pttrn; ss >> x; ss >> attrib; } } } fin.close(); } void findparam(const char *pf, int parnum, string &pttrn, double &x, string &attrib) { ifstream fin; fin.open(pf); int number = 0; while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { number ++; if(number == parnum) { stringstream ss; ss << line; ss >> pttrn; ss >> x; ss >> attrib; } } } fin.close(); } void findparam(const char *pf, int parnum, string &pttrn, string &x, string &attrib) { ifstream fin; fin.open(pf); int number = 0; while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0;

32

| 6. GELOMBANG TALI

bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { number ++; if(number == parnum) { stringstream ss; ss << line; ss >> pttrn; ss >> x; ss >> attrib; } } } fin.close(); } int countparam(const char *pf) { ifstream fin; fin.open(pf); int number = 0; while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) {number ++;} } fin.close(); return number; } int countparam(const char *pf, string attrib) { ifstream fin; fin.open(pf); int number = 0; while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { stringstream ss; string col1; string col2; string col3; ss << line; ss >> col1; ss >> col2; ss >> col3;

6.8. LAMPIRAN

33

if(col3 == attrib) {number ++;} } } fin.close(); return number; } void writeparam(const char *pf, const char *pttrn, int x) { fstream fout; fout.open(pf, fstream::out | fstream::app); fout << pttrn << "\t"; fout << x << endl; fout.close(); } void writeparam(const char *pf, const char *pttrn, double x) { fstream fout; fout.open(pf, fstream::out | fstream::app); fout << pttrn << "\t"; fout << x << endl; fout.close(); } void writeparam(const char *pf, const char *pttrn, int N, double *x) { fstream fout; fout.open(pf, fstream::out | fstream::app); fout << pttrn << "\t"; fout << N << endl; for(int i = 0; i < N; i++) { fout << x[i] << endl; } fout.close(); } void writeparam(const char* pf, const char* x) { fstream fout; fout.open(pf, fstream::out | fstream::app); fout << x << endl; fout.close(); } #endif

Berikut ini adalah isi dari library ’vect2.h’ yang digunakan untuk mendefinisikan jenis variabel vektor 2 dimensi, ’vect2’:

34 #include #include <sstream> #include using namespace std; #ifndef vect2_h #define vect2_h // Structure for 2-d vector struct vect2 { double x; double y; vect2(void); vect2(const vect2&); vect2(double, double); ~vect2(void); string strval(void); }; // Operator overloading vect2 operator+(vect2, vect2); vect2 operator-(vect2, vect2); double operator|(vect2, vect2); double length(vect2); vect2 operator*(double, vect2); vect2 operator*(vect2, double); vect2 operator/(vect2, double); vect2 operator>>(vect2, double); vect2 operator+=(vect2&, vect2); vect2 operator-=(vect2&, vect2); vect2 operator*(vect2, vect2); vect2 operator/(vect2, vect2); // Implementation vect2::vect2(void) { x = 0.0; y = 0.0; } vect2::vect2(const vect2& z) { x = z.x; y = z.y; } vect2::vect2(double zx, double zy) { x = zx; y = zy; }

| 6. GELOMBANG TALI

6.8. LAMPIRAN vect2::~vect2(void) { x = 0.0; y = 0.0; } string vect2::strval(void) { stringstream ss; ss << "(" << x << ", " << y << ")"; return ss.str(); } vect2 operator+(vect2 a, vect2 b) { vect2 c; c.x = a.x + b.x; c.y = a.y + b.y; return c; } vect2 operator-(vect2 a, vect2 b) { vect2 c; c.x = a.x - b.x; c.y = a.y - b.y; return c; } double operator|(vect2 a, vect2 b) { double c = a.x * b.x + a.y * b.y; return c; } double length(vect2 r) { double c = sqrt(r|r); return c; } vect2 operator*(double a, vect2 b) { vect2 c; c.x = a * b.x; c.y = a * b.y; return c; } vect2 operator*(vect2 a, double b) { vect2 c; c.x = a.x * b; c.y = a.y * b; return c; } vect2 operator/(vect2 r, double l) {

35

36

| 6. GELOMBANG TALI

vect2 s; s.x = r.x / l; s.y = r.y / l; return s; } vect2 operator>>(vect2 r, double s) { double l = length(r); vect2 u = (l > 0) ? r / l : vect2(0, 0); return u; } vect2 operator+=(vect2& a, vect2 b) { a = a + b; } vect2 operator-=(vect2& a, vect2 b) { a = a - b; } double double double double return }

angle(vect2 a, vect2 b) { la = length(a); lb = length(b); c = acos((a|b)/(la*lb)); c;

double double double double

angle(vect2 a, vect2 b, string option) { la = length(a); lb = length(b); c = acos((a|b)/(la*lb));

if (option == "DEG") {c = (180/M_PI)*c;} return c; } vect2 operator*(vect2 a, vect2 b) { vect2 c; c.x = a.x * b.x; c.y = a.y * b.y; return c; } vect2 operator/(vect2 a, vect2 b) { vect2 c; c.x = a.x / b.x; c.y = a.y / b.y; return c; } #endif

6.8. LAMPIRAN

6.8.4

37

Program Modifikasi Parameter

Program ini bertujuan untuk melakukan modifikasi parameter masukan dari program tali. Berikut ini adalah code C++ dari program tersebut ’paramsmodifier.cpp’, yaitu: #include <sys/types.h> #include <sys/stat.h> #include #include #include #include #include #include

<sstream> <process.h>

#include "rwparams.h" using namespace std;

void plotresults(string rd, string dn, int nop, string xl, string yl, string *pn, string // Create gnuplot pipe to create graphs FILE *gnuplot; // Construct data file string string dl = rd + "/" + dn + ".txt"; // Check the existance of data file ifstream fin; fin.open(dl.c_str()); if(fin.good()) { // Create graphs of amplitude results gnuplot = popen("gnuplot -persist", "w"); fprintf(gnuplot, "reset\n"); fprintf(gnuplot, "cd \"%s\"\n", rd.c_str()); fprintf(gnuplot, "set term %s\n",ftyp.c_str()); fprintf(gnuplot, "set output \"%s.%s\"\n", dn.c_str(), ftyp.c_str()); // Find maximum y to adjust plot boundary fprintf(gnuplot, "ymax = 0\n"); fprintf(gnuplot, "set yrange [0:*]\n"); fprintf(gnuplot, "do for [i = 2:%i] {\n", nop + 1); fprintf(gnuplot, "plot \"%s.txt\" using 1:i\n", dn.c_str()); fprintf(gnuplot, "if (ymax < GPVAL_Y_MAX) {ymax = GPVAL_Y_MAX}\n"); fprintf(gnuplot, "set output \"%s.%s\"\n}\n", dn.c_str(), ftyp.c_str()); // Set x and y properties fprintf(gnuplot, "set yrange [0:1.5*ymax]\n");

38

| 6. GELOMBANG TALI

fprintf(gnuplot, "set xlabel \"%s\"\n", xl.c_str()); fprintf(gnuplot, "set ylabel \"%s\"\n", yl.c_str()); // Plot the results fprintf(gnuplot, "plot "); for(int i = 0; i < nop; i++) { fprintf(gnuplot, "\"%s.txt\" ", dn.c_str()); fprintf(gnuplot, "using 1:%i ", i + 2); fprintf(gnuplot, "title \"%s %s\"", pn[i].c_str(), pnpre.c_str()); if(i == nop - 1) {fprintf(gnuplot, "\n");} else {fprintf(gnuplot, ",");} } fprintf(gnuplot, "set output\n"); fprintf(gnuplot, "\n"); fclose(gnuplot); // Give information: image is created cout << rd << "/" << dn << "." << ftyp << " is created." << endl << endl; } else { // Give error result: data file does not exist cout << dl << " does not exist." << endl; cout << "Try to evaluate the data to get the results." << endl << endl; } // close result data file fin.close(); } int main(int argc, char *argv[]) { // Get compiled name const char *arg0 = argv[0]; const char *pname = strstr(arg0, "./"); if(pname != NULL) pname += 2; // Check number if(argc < 4) { cout << "Usage: cout << "arg[1] cout << "arg[2] cout << "arg[3] exit(1); }

of argument " << pname << " [pfile]" << endl; -> program_directory" << endl; -> program_file_name" << endl; -> params_file_name" << endl;

// Get directory and name information from argument // Notes : program & params file must be in the same directory string progdir = argv[1]; string prognam = argv[2];

6.8. LAMPIRAN string parfnam = argv[3]; // Define program and parameter file name string string prog = progdir + "/" + prognam; string parf = progdir + "/" + parfnam; // Look parameters from program directory lookparam(parf.c_str()); int nop = countparam(parf.c_str()); // Read output preposition parameter string ofnpre = ""; readparam(parf.c_str(), "OUTPRE", ofnpre); // Define output parameter file name string string opfns = ofnpre + "-temp-" + parfnam; // Show all variable parameters cout << endl << "List of variable parameters" << endl; for(int i = 1; i < nop + 1; i++) { string parnam = ""; string parval = ""; string paratr = ""; findparam(parf.c_str(), i, parnam, parval, paratr); if(paratr == "VARIABLE") { cout << i << ".\t" << parnam << endl; } } // Ask user to select an input parameter to be varied int selpar = 0; cout << "Select a variable input parameter to be varied : "; cin >> selpar; cout << endl; // Get information about selected parameter string selnam = ""; string selval = ""; string selatr = ""; findparam(parf.c_str(), selpar, selnam, selval, selatr); if(selatr != "VARIABLE") { cout << "Error: parameter selected is not variable" << endl; exit(1); } // Define result files name string string resdir = "result"; string ampvsvarn = "result_AMPLITUDE_VS_" + selnam; string tenvsvarn = "result_TENSION_VS_" + selnam; string angvsvarn = "result_ANGLE_VS_" + selnam;

39

40

| 6. GELOMBANG TALI

string ampvsvar = resdir + "/" + ampvsvarn + ".txt"; string tenvsvar = resdir + "/" + tenvsvarn + ".txt"; string angvsvar = resdir + "/" + angvsvarn + ".txt"; // Count output parameters int nrp = countparam(opfns.c_str()); int nampp = countparam(opfns.c_str(), "AMPLITUDE"); int ntenp = countparam(opfns.c_str(), "TENSION"); int nangp = countparam(opfns.c_str(), "ANGLE"); // Construct output parameter title string array string tampp[nampp], ttenp[ntenp], tangp[nangp]; int nampc = 0, ntenc = 0, nangc = 0; for(int i = 1; i < nrp + 1; i++) { string parnam = ""; string parval = ""; string paratr = ""; findparam(opfns.c_str(), i, parnam, parval, paratr); if(paratr == "AMPLITUDE") { tampp[nampc] = parnam; nampc++; } if(paratr == "TENSION") { ttenp[ntenc] = parnam; ntenc++; } if(paratr == "ANGLE") { tangp[nangc] = parnam; nangc++; } } // Repeat the option available to do while(1) { // Choose option to do int option = 0; cout << "Choose task to do" << endl; cout << "1. Evaluate data to get results" << endl; cout << "2. Plot the results" << endl; cout << "3. Exit" << endl; cout << "Choose option no: "; cin >> option; cout << endl; // Executing option no. 1: evaluating the data if(option == 1) { // Ask user to input initial, step, and final value of variable double varinit = 0;

6.8. LAMPIRAN double varstep = 0; double varend = 0; cout << selnam << " initial value = "; cin >> varinit; cout << selnam << " final value = "; cin >> varend; cout << selnam << " step value = "; cin >> varstep; cout << endl; // Compare initial variable value with final variable value if(varinit > varend) { double dum = varinit; varinit = varend; varend = dum; } // Open result files ofstream ampout, tenout, angout; ampout.open(ampvsvar.c_str()); tenout.open(tenvsvar.c_str()); angout.open(angvsvar.c_str()); // Create ampout << tenout << angout << for(int i for(int i for(int i ampout << tenout << angout <<

header for result files "# " << selnam; "# " << selnam; "# " << selnam; = 0; i < nampp; i++) {ampout << "\t" << tampp[i];} = 0; i < ntenp; i++) {tenout << "\t" << ttenp[i];} = 0; i < nangp; i++) {angout << "\t" << tangp[i];} endl; endl; endl;

// Selected variable parameter iteration double varval = varinit; while(varval < varend + (0.5 * varstep)) { // Convert variable value into string stringstream varcon; string varvals = ""; varcon << varval; varcon >> varvals; // Create temporary parameters file string temp = "temp-" + parfnam; ofstream tempar; tempar.open(temp.c_str()); for(int i = 1; i < nop + 1; i++) { string parnam = ""; string parval = ""; string paratr = ""; findparam(parf.c_str(), i, parnam, parval, paratr); if(i == selpar) {parval = varvals;}

41

42

| 6. GELOMBANG TALI

if(parnam == "CREANI") {parval = "0";} tempar << parnam << "\t" << parval << "\t" << paratr << endl; } tempar.close(); // Executing program using temporary parameters file cout << "Progress (" << (varval - varinit + varstep)/varstep << "/"; cout << (varend - varinit + varstep)/varstep << "): Executing "; cout << prognam.c_str() << " for " << selnam << " = " << varvals; spawnl(P_WAIT, prog.c_str(), prognam.c_str(), temp.c_str(), NULL); // Save output parameters into result files ampout << varvals; tenout << varvals; angout << varvals; for(int i = 1; i < nrp + 1; i++) { string parnam = ""; string parval = ""; string paratr = ""; findparam(opfns.c_str(), i, parnam, parval, paratr); if(paratr == "AMPLITUDE") {ampout << "\t" << parval;} if(paratr == "TENSION") {tenout << "\t" << parval;} if(paratr == "ANGLE") {angout << "\t" << parval;} } ampout << endl; tenout << endl; angout << endl; // Increase variable value varval += varstep; } // Close result files ampout.close(); tenout.close(); angout.close(); // Continue to the next option continue; } else if(option == 2) { // Executing option no. 2: plot the results // Plot the resulted data into plotresults(resdir, ampvsvarn, plotresults(resdir, tenvsvarn, plotresults(resdir, angvsvarn, // Continue to the next option continue;

eps image file nampp, selnam, "AMPLITUDE", tampp, "AMPLITUDE", ntenp, selnam, "TENSION", ttenp, "TENSION", "ep nangp, selnam, "ANGLE", tangp, "ANGLE", "eps");

6.8. LAMPIRAN

43

} else {break;} // Last option: break the menu option iteration } // Give information: the program ’paramsmodifier’ has done cout << endl << "end of paramaters modifier" << endl; return 0; } Kode C++ di atas dikompilasi menggunakan g++ dari MinGW melalui Command Prompt dengan cara g++ paramsmodifier.cpp -o paramsmodifier dan dipanggil melalui Command Prompt dengan cara paramsmodifier program tali.exe params.txt di mana ’program’ merupakan folder yang berisi program ’tali.exe’ dan segala file kelengkapannya, ’tali.exe’ adalah nama program yang akan dipanggil, dan ’params.txt’ merupakan nama dari file parameter program. Pada saat kompilasi, pastikan file ’rwparams.h’ berada dalam direktori yang sama dengan file ’paramsmodifier.cpp’. Pada saat program dipanggil, pastikan folder ’data’, folder ’program’, folder ’result’, file ’temp-params.txt’ dan ’out-tempparams.txt’ berada dalam direktori yang sama dengan ’paramsmodifier.exe’. Untuk mencegah duplikasi data, kosongkan folder ’data’ setiap kali pemanggilan program.

6.8.5

Script Gnuplot

Script ini digunakan ungtuk membuat cuplikan gambar dari hasil simulasi ’tali.exe’ dengan sumbu x dan y memiliki skala yang menyesuaikan dengan besar ukuran tali. Script ini digunakan untuk gambar 6.3 dan 6.4. # Set terminal type and output file set term eps set output "fig-03a.eps" # Set y range set yrange [-0.001:0.001] # Set title set title "t = 0.0250 s"

44

| 6. GELOMBANG TALI

# Execute plotting plot "fig-03a.txt" notitle set output Script gnuplot yang berikut ini digunakan untuk memplot kurva parameter keluaran untuk gambar 6.5 sampai 6.8. # Set terminal type and output file set terminal postscript eps enhanced color set output "fig-07.eps" # Set legend position set key left # Set x & y axis set xlabel "Amplitudo Osilator (m)" set ylabel "Tegangan Rata-Rata (X 10 N)" set xrange [0:0.0105] set yrange [0:*]

# Execute plotting plot "fig-07.txt" using 1:2 w l title "ambang", "fig-07.txt" using 1:3 w p tit set output Script Gnuplot di atas dipanggil dengan cara gnuplot script.gps Pastikan file data yang akan diplot dan file ’script.gps’ berada dalam direktori yang sama.

|7

Pemodelan Lintasan Benda Titik pada Tong Setan Rustan | [email protected] Erika Linda Yani NST | [email protected] Wenny Wahyuni | [email protected] Miftahul Husnah | [email protected] Tong Setan (Wall of Death) merupakan istilah untuk salah satu wahana hiburan di pasar malam, berupa atraksi seseorang yang berputar mengelilingi sebuah tong silinder berukuran besar menggunakan motor bahkan mobil. Anehnya, meskipun dinding tong memiliki kemiringan hampir tegak lurus dengan tanah, pengendara dapat berputar-putar tanpa terjatuh. Banyak pendapat mengenai fenomena tersebut, ada berpendapat pemain atraksi menggunakan ilmu magik, sulap, atau memang ada bantuan setan. Fenomena tersebut menarik untuk dikaji dan dijelaskan secara fisika dan dimodelkan untuk meluruskan pendapat-pendapat tersebut. Pada RBL ini, sebuah motor diasumsikan sebagai sebuah benda titik, bergerak dengan kecepatan sudut ω bergerak mengelilingi dinding silinder. Untuk bergerak stabil dan tidak terjatuh, benda titik harus melewati kecepatan sudut ωmin . Persamaan gerak akan diturunkan secara analitik dan numerik menggunakan metode Euler.

7.1

Teori

Suatu benda yang bergerak pada dinding vertikal kasar, maka gaya-gaya yang bekerja antara lain: gaya gravitasi, gaya gesek, dan gaya normal. 1

2

| 7. PEMODELAN LINTASAN BENDA TITIK PADA TONG SETAN

Gambar 7.1: Gaya-gaya yang bekerja pada benda di dinding silinder

Pada sumbu x, v2 R v2 m R

ΣFx

=

N

=

ΣFz W − fs

= =

maz maz

mg − µk N

=

az

=

maz mg − µk N m

m

(7.1)

Pada sumbu z,

(7.2)

Subtitusi persamaan (7.1) ke persamaan (7.2) sehingga:

az

=

az

=

az

=

mg − (µk mv 2 /R) m µk v s 2 g− R g − µk ω 2 R

(7.3)

Saat ω ≥ ωmin , maka: az

=

0

(7.4)

3

7.2. SYARAT BATAS

7.2

Syarat Batas

Agar benda tetap pada lintasan dan tidak jatuh, maka dalam arah sumbu vertikal gaya gesek, fs harus sama dengan gaya berat, mg: fs

=

mg

µs N

=

N

=

mg mg µs

(7.5)

Dalam arah horizontal, gaya yang bekerja adalah gaya sentripetal. Gaya sentripetal adalah resultan gaya yang arahnya menuju/keluar dari pusat silinder. ΣFs = m

v2 R

(7.6)

Pada kasus ini, gaya sentripetal yang bekerja hanya gaya normal, N maka: N =m

v2 R

(7.7)

Subtitusi persamaan (7.5) ke persamaan (7.7) : mg µs

=

v

=

ω

=

v2 m R s gR µs r g µs R

(7.8)

dimana ω adalah kecepatan sudut minimal yang harus dilewati benda titik agar tetap pada lintasan (tidak terjatuh), g adalah percepatan gravitasi, µs adalah koefisien gesek dinding tong, dan R adalah jari-jari silinder.

7.3


Dengan menggunakan metode Euler, maka persamaan (7.3) dapat memberikan kecepatan dan posisi setiap saat dalam arah z untuk ω < ωmin : az (t) vz (t + ∆t)

= =

g − µk ω 2 R

vz (t) + az (t)∆t

(7.9) (7.10)

sedangkan persamaan posisi pada sumbu x, y, dan z yaitu: x(t) y(t)

= =

Rsin(ωt) Rcos(ωt)

(7.11) (7.12)

z(t + ∆t)

=

z(t) + vz (t)∆t

(7.13)

4


Untuk ω ≥ ωmin , persamaan posisi pada sumbu x, y, dan z yaitu: x(t)

=

Rsin(ωt)

(7.14)

y(t) z

= =

Rcos(ωt) z0

(7.15) (7.16)

Dalam implementasinya Persamaan-persamaan (7.9) - (7.16) dituangkan dalam algoritma berikut ini : L01 . t0 , t1 , t2 , ∆t, ω0 , ω1 , ω2 ? L02 . g, R, z0 , vz , µ, θ? L03 . t = t0 −ω0 L04 . ω(t) = ωt11 −t (t − t0 ) + ω0 0 L05 . θ = θ + ω(t)∆t L06 . t = t + ∆t L07 . t ≤ t1 → L04 −ω1 (t − t1 ) + ω1 L08 . ω(t) = ωt22 −t 1 L09 . θ = θ + ω(t)∆t L10 . t = t + ∆t L11 . t ≤ t2 → L08 L12 . ω = ω0 L13 . x = Rsinθ L14 . y = Rcosθ L15 . az = 0 L16 . vz = 0 L17 . z = z0 L18 . ω(t) = ω0 − ∆ω L19 . ω(t) ≥ ωmin → L13 L20 . x = Rsinθ L21 . y = Rcosθ L22 . az = −0.5(g − µω 2 R) L21 . vz (t + ∆t) = vz (t) + az (t)∆t L22 . z(t + ∆t) = z(t) + vz (t)∆t L23 . ω(t + ∆t) = ω(t) + ∆ω L24 . ω(t) < ωmin → L20

7.4

Hasil dan Diskusi

Nilai-nilai parameter yang digunakan dalam simulasi tercantum dalam tabel 1.1. berikut: Hasil simulasi yang diperoleh ditampilkan pada gambar 1.2. Kurva hubungan kecepatan sudut dengan waktu ditampilkan dalam gambar 1.3. Kurva hubungan ketinggian dengan waktu ditampilkan dalam gambar 1.4.

5


Tabel 7.1: Nilai-nilai parameter simulasi. P arameter Jari-jari silinder (R) Ketinggian awal (z0 ) Gaya gravitasi (g) Kec sudut minimum (ωmin ) Kecepatan sudut awal (ωawal ) Kecepatan sudut akhir (ωakhir ) takhir tawal Koefisien gesek (µ)

N ilai 7 10 10 5 10 1 5 0 0.04

Gambar 7.2: Hasil Simulasi, dengan parameter R = 7, z0 = 10, g = 10, wmin = 5, µ = 0.04, tbeg = 0, ∆t = 0.001, tmaks = 5

Gambar 7.3: Kecepatan sudut vs waktu, dengan parameter t1 = 2, t0 = 0, t2 = 5dant1 = 3, t0 = 0, t2 = 5

6


Gambar 7.4: Kecepatan sudut vs waktu, dengan parameter t1 = 2, t0 = 0, t2 = 5dant1 = 3, t0 = 0, t2 = 5

7.5

Ringkasan

Telah diperoleh bahwa lintasan benda titik pada dinding silinder vertikal akan stabil (tidak turun) jika kecepatan sudut benda melebihi atau sama dengan kecepatan sudut minimum, ω ≥ ωmin . Sebaliknya jika kecepatan sudut benda lebih kecil dari kecepatan sudut minimum, ω < ωmin , maka benda akan bergerak turun sebesar z tertentu.

7.6

Referensi

1. Abramowick, M.A. dan Szuszkiewicz, E. 1993. ”The Wall of Death”. American Journal of Physics. http://dx.doi.org/10.1119/1.17349 2. Ahmed, M. 2012. Methods in Computational Physics. North Carolina University. 3. Sparisoma Viridi. ”Visualisasi Data dengan Gnuplot: Modul Praktikum FI2283 Pemrograman dan Simulasi Fisika”. Versi 29 September 2013.

7.7

Lampiran

7.7.1

Kode C++

/* coba3.cpp melihat gerakan dari motor dalam dimensi x, y, z

7.7. LAMPIRAN tugas rbl komfis */ #include #include #include <math.h> #include using namespace std; ofstream fout; double omega (double t){ double t0 = 0; double t1 = 3; double t2 = 5; double w0 = 10; double w1 = 1; double w2 = w0; double w = 0; if ((t0<=t)&&(t<=t1)){ w = (((w1-w0)/(t1-t0))*(t-t0))+w0; } else { w = (((w2-w1)/(t2-t1)*(t-t1))+w1); } return w; } int main(){ const char *ofn = "data2.txt"; fout.open(ofn); fout << "#t\tx\ty\tz" << endl; double t=0; double tmax=5; double dt=0.001; int N=(tmax-t)/dt; //double theta = theta + omega (t) * dt; double wmin = 5; double double double double double double double double

g = 10; r = 7; z = 10; x, y; vz = 0; theta = 0; muk = 0.04; w = 0;

cout << "N = " << N << endl; for (int i = 0; i < N; i++) { w = omega(t); theta = theta + w * dt;

7

8


t=t+dt; if(w>=wmin){ x = r*sin(theta); y = r*cos(theta); double az = 0; vz = 0; z = z + vz * dt; } else { x = r*sin(theta); y = r*cos(theta); double az = -0.5*(g-muk*w*w*r); vz = vz + az * dt; z = z + vz * dt; } //cout << w << "\t" << x << " \t" << y << " \t" << z << endl; cout << t << "\t"; cout << x << "\t"; cout << y << "\t"; cout << z << "\t"; cout << w << "\t"; cout << wmin << endl; fout fout fout fout fout }

<< << << << <<

t x y z w

<< << << << <<

"\t"; "\t"; "\t"; "\t"; endl;

fout.close(); return 0; } Kode C++ di atas dikompilasi dengan cara g++ coba3.cpp -o coba3 dan dipanggil dengan cara ./code > data2.txt

7.7.2

Script Gnuplot

# plot.gps

7.7. LAMPIRAN # Plot data and produce EPS file # rustan # 20151216 # Execute: gnuplot plot.gps # Set terminal to eps set output "gambar_02.eps" # Set label and grid set xlabel "{/Times-Oblique posisi x (meter)}" set ylabel "{/Times-Oblique posisi y (meter)}" set zlabel "{/Times-Oblique ketinggian (meter)} set grid # Execute plotting splot ’data2.txt’ using 2:3:4 Script Gnuplot di atas dipanggil dengan cara gnuplot plot.gps

9

10


|8

Osilasi Teredam Sistem Pegas dengan Massa Beban tak Konstan Kamirul | [email protected] Lisa’ Yihaa Roodhiyah | [email protected] Hezliana Syahwanti | [email protected] Fatriani | [email protected] Peristiwa atau fenomena fisis sederhana sering terjadi disekitar kita dan menarik untuk dikaji. Salah satu contohnya adalah gerak osilasi pegas non-harmonik. Selama ini kita sering mengamati ataupun mensimulasikan gerak pegas harmonik yaitu osilasi tanpa adanya gesekan udara maupun pengurangan massa beban yang digantungkan pada pegas tersebut. Berbanding terbalik dengan apa yang terjadi pada keadaan sebenarnya, justru sangat sulit untuk menemukan bahkan mengkondisikan berosilasi harmonik. Sebelumnya Rodrigues, et.al.(2014) telah mensimulasikan hal serupa namun tidak memasukkan efek gaya gesekan udara pada persamaan tersebut. Oleh karena itu, didalam penelitian ini akan disimulasikan osilasi pegas non-harmonik sebagai akibat adanya faktor gesekan udara. Selain itu, akan ditambahkan pula efek perubahan (pengurangan) massa beban terhadap perilaku gerak osilasi sistem pegas non harmonik. Pada proses simulasi, pengurangan massa beban dikondisikan dengan menggunakan beban berupa wadah berisi air dan dilubangi dibagian dasarnya. Tujuan pokok dari penelitian ini, selain mensimulasikan lintasan gerak osilasi pegas non harmonik adalah mempelajari pengaruh besarnya lubang kebocoran wadah terhadap gerak sistem dan waktu yang dibutuhkan sistem untuk teredam sempurna. Selanjutnya hasil tersebut dikarakterisasi untuk mencari hubungan yang sesuai antara rasio kebocoran dan waktu yang dibutuhkan sistem pegas untuk berhenti berosilasi.

1

2| 8. OSILASI TEREDAM SISTEM PEGAS DENGAN MASSA BEBAN TAK KONSTAN

8.1

Teori

Hukum Newton kedua secara umum dapat dituliskan sebagai berikut (Giancoli, 2009):

F = ma(t)

(8.1)

Berdasarkan persamaan (8.1), a(t) dapat dijabarkan sehingga didapatkan :

F =m

dv dt

(8.2)

dan variasi bentuk lainnya :

F =

d mv dt

(8.3)

Dari persamaan (8.3) diatas, terdapat hubungan antara gaya dengan perubahan momentum linear, sehingga dapat didefinisikan (Giancoli, 2009):

F =

dP dt

(8.4)

Hukum ke-2 Newton bersifat invarian terhadap transformasi Galilean yang didefinisikan oleh x′ = x − ut dan v ′ = v − u. Dengan mengaplikasikan transformasi Galilean pada hukum ke-2 Newton maka akan didapatkan F ′ = F yang menyatakan bahwa semua pengamat harusnya mengukur atau merasakan gaya yang sama tanpa memperhatikan kecepatan relatif mereka. Untuk sistem yang bergerak dengan massa yang tidak konstan, perlunya menganalisa penerapan prinsip konservasi momentum linier pada sistem. Besar perubahan momentum linier sistem terhadap waktu adalah : ∆v ∆m ∆m ∆p =m + ∆v − (w − v) ∆t ∆t ∆t ∆t

(8.5)

Dengan menerapkan limit ∆t → 0, ∆m → 0, dan ∆v → 0, maka persamaannya (8.5) menjadi :

F =

dp dv dm =m − (w − v) dt dt dt

(8.6)

Untuk sistem yang bergerak sambil mengalami pengurangan massa, dm dt < 0, maka didapatkan persamaan gaya sebagai berikut (Rodrigues, et.al, 2014): dv dm =F− (w − v) dt dt d dm F = − (mv) − w dt dt

m

(8.7)

3

8.1. TEORI

Untuk mensimulasikan pengaruh pengurangan massa beban terhadap gerak sistem pegas, digunakan beban berupa wadah berisi air yang dilubangi dibagian dasarnya. Ilustrasi diagram gaya sistem pegas ditampilkan seperti Gambar 8.1. xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

k

Fg = -bv

h0

mw(t) + mb aw

A

w = mg

Gambar 8.1: Ilustrasi diagram gaya sistem pegas. Diasumsikan bahwa arah bocornya air dan gerak osilasi berada pada sumbu yang sama, yaitu sumbu z. Pada kondisi ini gaya-gaya yang bekerja pada sistem adalah gaya berat, gaya pegas, dan gaya gesekan udara. Gerak osilasi sistem diberikan oleh persamaan (8.8) sebagai berikut : (Rodrigues, et. al., 2014): mv˙ =

dm (w − v) − kz(t) − bv(t) − mg dt

(8.8)

Massa m merupakan total massa sistem yang terdiri dari massa wadah(mb ) dan massa air(mw ). Sedangkan z(t) adalah perpindahan pusat massa beban dihitung dari posisi kesetimbangan awal ze . Nilai ze dapat dihitung dengan persamaan : mb g ze = − (8.9) k . Parameter w merupakan kecepatan rata-rata air keluar dari lubang, dan v = z˙ adalah kecepatan osilasi pegas. Parameter b, k, dan g masing-masing adalah konstanta gesekan udara, konstanta pegas, dan percepatan gravitasi bumi 9,81 m/s2 . Selama proses osilasi, massa air mw akan berkurang sehingga merupakan fungsi dari waktu t dan diberikan oleh persamaan sebagai berikut : r g 2 ) (8.10) mw (t) = mw (0)(1 − f t 2h0 Parameter h0 adalah tinggi awal kolom air didalam wadah dan f = aAw adalah perbandingan antara diameter lubang kebocoran (aw ) dengan diameter dasar wadah (A). Asumsikan bahwa kecepatan air mengalir sangat kecil, sehingga nilai dm dt (w − v) dapat diabaikan. Persamaan kecepatan sebagai fungsi t dapat dituliskan sebagai berikut (Rodrigues, et. al., 2014): v˙ =

−k(z − z0 ) bv(t) −g− mb + mw (t) mb + mw (t)

(8.11)

4| 8. OSILASI TEREDAM SISTEM PEGAS DENGAN MASSA BEBAN TAK KONSTAN Ketika massa air berkurang selama osilasi, maka akan tercapai satu keadaan dimana massa air akan habis, sehingga pegas akan berosilasi hanya dipengaruhi oleh massa wadah mb saja. Waktu yang dibutuhkan untuk air keluar seluruhnya dinyatakan dengan τ diberikan oleh persamaan (8.12) berikut (Rodrigues, et. al., 2014): s τ=

1 f

2h0 g

(8.12)

Setelah berosilasi melebihi waktu pengosongan τ , maka pegas akan berosilasi dengan percepatan : bv(t) −k(z − z0 ) v˙ = −g− (8.13) mb mb

8.2

Syarat batas

Pada penelitian ini, arah gerak osilasi pegas hanya diperhitungkan pada arah z saja, yang mana artinya arah kecepatan dan percepatan juga searah dengan arah gerak pegas tersebut. Selain itu, syarat batas lain adalah laju kebocoran air tidak dipengaruhi oleh posisi ketinggian beban (energi potensial). Hal ini mengakibatkan jumlah air yang terdapat pada wadah hanya bergantung pada t saja.

8.3


Persamaan (8.13) merepresentasikan gerak harmonik sederhana dengan gaya gesekan yang memiliki solusi eksak berupa fungsi sinusoidal. Akan tetapi, mendapatkan solusi eksak untuk persamaan (8.11) tidak mudah, sehingga di dalam paper ini simulasi program dilakukan dengan pendekatan numerik. Berdasarkan persamaan (8.11), terlihat bahwa nilai v dan v˙ bergantung pada nilai terakhir dari z(t) dan mw (t). Dengan melakukan diskritisasi terhadap waktu osilasi pegas ti = t0 + ih, maka kita dapat menghitung v(ti+1 ) dari v(ti ). Untuk mendapatkan nilai vi+1 maka dapat digunakan pendekatan v˙ ≈ (vi+1 − vi )/h. Pemilihan nilai selang waktu h yang baik adalah sekecil mungkin, namun harus disesuaikan dengan kemampuan komputer untuk menjalankan simulasi. Dengan cara yang sama pula, kita dapat menghitung posisi z(ti+1 ) dari z(ti ) dengan menggunakan aproksimasi z˙ ≈ (zi+1 − zi )/h. Pada penerapannya, untuk interval simulasi 0 ≤ t ≤ τ , maka nilai v(t) diberikan oleh persamaan berikut : vi+1 = vi − h(

−k(z(ti ) − z0 )) bv(t) +g+ ) mb + mw (ti ) mb + mw (ti )

(8.14)

Setelah pegas berosilasi melebihi waktu pengosongan t ≥ τ , maka persamaan (8.14) menjadi : vi+1 = vi − h(

bv(t) −k(z(ti ) − z0 )) +g+ ) mb mb

(8.15)

5

8.4. HASIL DAN DISKUSI dan untuk menghitung z(t) dapat diselesaikan dengan persamaan : zi+1 = zi − hvi

(8.16)

Sedangkan posisi awal beban z(0) yang digunakan untuk memulai proses simulasi dapat dihitung dengan persamaan: z(0) = −

mb + mw0 g k

(8.17)

Selain itu, posisi beban juga dapat ditentukan dengan hanya mengetahui nilai massa total beban tersebut yang dirumuskan dengan : z¯ = −

mb + m w g k

(8.18)

Karena adanya pengaruh gaya gesekan, maka setelah berosilasi dalam rentang waktu tertentu pegas akan berhenti dan kembali pada posisi kesetimbangan z0 . Untuk mendefinisikan kondisi berhenti pegas digunakan sebuah nilai dz yang dirumuskan sebagai : dz =| zi+1 − zi | (8.19) Apabila nilai dz lebih kecil dari nilai batasan berhenti yang diberikan Cstop , maka simulasi dihentikan. Implementasi program dituangkan dalam bentuk algoritma sebagai berikut: L01. L02. L03. L04.

k, b, g, Cstop ? mb , mw0 , h0 , a, A? z0 = − mkb g w g z = − mb +m qk 2h0 1 L05. τ = − f g L06. t = tb eg q L07. Jika t < τ , maka mw (t) = mw0 (1 − f t 2hg 0 )2

Jika tidak, maka mw (t) = 0 −k(z−z0 ) −g− L08. a(t + ∆t) = m +mw (t) b L09. v(t + ∆t) = v(t) + a(t)∆t L10. z(t + ∆t) = z(t) + z(t)∆t L11. t = t + ∆t L12. t ≤ tend → L07.

8.4

bv(t) mb +mw (t)

Hasil dan diskusi

Tahapan simulasi awal dilakukan menggunakan parameter-parameter yang diberikan pada Tabel 10.1. Posisi awal dan posisi kesetimbangan sistem beban dihitung dengan persamaan 8.17 dan 8.9. Ketika lubang kebocoran masih tertutup, posisi beban dan tempatnya akan berada di z(0) sesuai dengan persamaan 8.17. Pada kondisi

6| 8. OSILASI TEREDAM SISTEM PEGAS DENGAN MASSA BEBAN TAK KONSTAN

Tabel 8.1: Nilai-nilai parameter simulasi. Parameter ∆t tbeg tend b k g mb mw h0 A aw

Nilai 0, 05 0 150 0, 18 100 9.81 1 10 0, 05 1 0, 01

ini sistem tidak mengalami percepatan sehingga akan tetap diam. Percepatan sistem murni hanya dikarenakan adanya perubahan massa sistem yaitu ketika lubang kebocoran mulai dibuka. Hasil simulasi awal ditampilkan pada Gambar 8.2. 0

-0.2

Posisi (m)

-0.4

-0.6

-0.8

-1

Z(t) Z(t)rata-rata Tau

-1.2 0

20

40

60

80

100

120

140

160

Waktu (s)

Gambar 8.2: h0 = 0, 05.

Hasil simulasi dengan parameter beban mb = 1, mw = 10 dan

Grafik dengan warna merah merupakan lintasan osilasi beban sebagai fungsi t, z(t). Terlihat bahwa setelah melewati waktu pengosongan τ (garis hijau tegak), sistem pegas akan berosilasi disekitar titik kesetimbangan ze (garis biru mendatar) dan selanjutnya teredam karena adanya pengaruh gaya gesek udara yang besarnya sebanding dengan −bv. Grafik berwarna biru merupakan beban sistem pegas yang didekati dengan persamaan (8.18) yang merupakan posisi beban sebagai fungsi massa total. Bentuk grafik ini memiliki kemiringan yang berbeda-beda bergantung pada besarnya parameter lubang kebocoran yang diterapkan.

7


0

-0.1

-0.2

-0.3

Posisi (m)

-0.4

-0.5

-0.6

-0.7

-0.8

-0.9 Z(t) rata-rata a=0.001 Z(t) rata-rata a=0.002 Z(t) rata-rata a=0.003 Z(t) rata-rata a=0.004 Z(t) rata-rata a=0.005

-1

-1.1 0

20

40

60 Waktu (s)

80

100

120

Gambar 8.3: Variasi lubang kebocoran berpengaruh pada kemiringan grafik posisi rerata beban. Hasil simulasi pada Gambar 8.3 memperlihatkan bahwa semakin besar lubang kebocoran yang diterapkan semakin kecil kemiringan yang dihasilkan. Terdapat hal yang menarik dari hubungan variasi lubang kebocoran terhadap bentuk lintasan rerata sistem beban. Dengan mendekati setiap lintasan tersebut menggunakan sebuah persamaan eksponensial yang z(t) = z(0)(1 − e)−αt , maka kita dapat menetukan besarnya nilai konstanta α masing-masing lintasan. Nilai α dapat digunakan untuk mengkarakterisasi besarnya lubang kebocoran suatu (aw) wadah yang sedang berosilasi dengan sistem serupa. Nilai α dapat ditentukan dengan mengetahui besarnya waktu yang dibutuhkan sistem pegas untuk berhenti berosilasi (Tstop ). Pada saat t = Tstop , maka posisi beban akan kembali ke posisi kesetimbangan dimana massa air sudah kosong yaitu di z = ze , sehingga berlaku hubungan : α=

ln(1 −

ze z(0) )

Tstop

(8.20)

Kondisi pemberhentian yang diterapkan adalah ketika besar perpindahan absolut dari z(t − ∆t) ke z(t) lebih kecil dari Cstop sesuai dengan persamaan (8.19). Idealnya, semakin besar lubang kebocoran yang diterapkan, semakin singkat waktu dibutuhkan oleh sistem pegas untuk berhenti berosilasi. Hal ini bersesuaian dengan hubungan antara rasio kebocoran terhadap waktu pengosongan yang diberikan pada persamaan 8.12. Hubungan tersebut ditampilkan pada Gambar 8.4. Simulasi selanjutnya dilakukan dengan memvariasikan besarnya nilai f dan mencatat waktu yang dibutuhkan beban untuk berhenti berosilasi (Tstop ). Pada simulasi ini, konstanta pemberhentian Cstop yang diterapkan adalah sebesar 10−10 m. Hasil simulasi sangat dipengaruhi oleh dikritisasi ∆t yang dipilih.

8| 8. OSILASI TEREDAM SISTEM PEGAS DENGAN MASSA BEBAN TAK KONSTAN 110

100

90

Waktu Pengosongan (s)

80

70

60

50

40

30

20

10 0.001

0.002

0.003

0.004

0.005 0.006 f (Rasio Kebocoran)

0.007

0.008

0.009

0.01

Gambar 8.4: Hubungan antara variasi rasio kebocoran (f ) dan waktu pengosongan (τ )(Rodrigues, et.al.(2014).

Hal ini dikarenakan pendefinisian berhentinya osilasi ditentukan oleh selisih perpindahan beban tersebut pada arah sumbu z. Parameter ∆t yang dipilih merupakan kelipatan dari periode terkecil osilasi sistem yaitu ketika massa air sudah kosong (T ). Apabila ∆t yang dipilih terlalu besar, maka akan mengurangi ketelitian serta resolusi posisi beban sebagai fungsi dari waktu t. Hal ini pula akan menyebabkan definisi pemberhentian tidak terdeteksi sehingga sistem pegas tidak berhenti berosilasi pada waktu yang seharusnya. Namun, apabila ∆t yang dipilih terlalu kecil, berarti besar perpindahan beban dari dari waktu t ke t + ∆t juga akan kecil. Terdapat kondisi dimana perubahan yang sangat kecil ini terdeteksi dan memenuhi syarat berhenti padahal sistem belum mencapai titik kesetimbangan ze . Oleh karena itu untuk melakukan proses simulasi ∆t yang dipilih harus tepat. Pemilihan ∆t yang tidak tepat dapat memunculkan kesalahan fatal yang mengakibatkan hasil simulasi tidak fisis. Sebagai contoh, pada kondisi tertentu, nilai f yang lebih besar menghasilkan ∆t yang lebih besar. Hal ini tentunya bertentangan dengan kondisi fisis yang seharusnya yaitu semakin besar nilai f semakin singkat waktu yang dibutuhkan untuk berhenti. Hasil simulasi menunjukkan bahwa dengan menerapkan nilai ∆t yang berkisar antara T/20−T/10 kondisi fisis akan terpenuhi (Gambar 8.5). Pada kondisi ini, nilai grafik nilai Tstop terhadap f mengikuti bentuk grafik hubungan τ terhadap f pada Gambar 8.4. Hal ini mengindikasikan bahwa dengan memilih ∆t pada rentang nilai tersebut hasil simulasi memenuhi kondisi fisis yang seharusnya. Ketika ∆t yang dipilih diluar rentang nilai tersebut, contohnya ∆t = T/5 dan ∆t = T/100, didapatkan nilai Tstop yang tidak fisis. Untuk ∆t = T /5 polanya hampir mirip dengan Gambar 8.4, namun masih terjadi ketidak sesuaian di dua

9

8.5. RINGKASAN 11000

10000

9000

8000

Tstop (s)

7000

6000

5000

4000

3000

2000

1000 0.00001

0.00002

0.00003

0.00004

0.00005 f (Rasio Kebocoran)

0.00006

0.00007

0.00008

0.00009

Gambar 8.5: Grafik hubungan Tstop terhadap f dengan ∆t = T/20.

titik data. Lebih jelasnya disajikan pada Gambar 8.6. Hal yang sama juga terjadi pada ∆t = T /100, pola grafik yang dihasilkan semakin tidak teratur (Gambar 8.7). Hal ini tentunya kembali lagi ke definisi pemberhentian yang sangat bergantung pada ∆t yang dipilih. Setelah mendapatkan nilai ∆t yang sesuai untuk proses simulasi, selanjutnya melakukan karakterisasi rasio kebocoran. Tiap rasio kebocoran menghasilkan Tstop yang berbeda sehingga menghasilkan nilai α yang berbeda pula sesuai dengan persamaan 8.20. Hasil plot nilai α terhadap f disajikan pada Gambar 8.8. Dari gambar 8.8 terlihat bahwa hubungan antara α dan f adalah linier dan memenuhi fungsi persamaan α(f ) = 0, 942f +(2e−8 ). Dengan adanya hubungan ini, maka kita dapat menghitung besarnya lubang keboc

8.5

Ringkasan

Dari hasil simulasi diperoleh bahwa nilai ∆t simulasi yang dipilih berpengaruh besar terhadap waktu berhenti (Tstop ) osilasi sistem pegas. Pemilihan ∆t yang tidak tepat dapat memunculkan makna tidak fisis yaitu lubang kebocoran yang besar mengakibatkan waktu berhenti osilasi pegas lebih lama. Hal ini tentunya tidak sesuai dengan keadaan yang sebenarnya. Dari hasil percobaan secara berulang-ulang, didapatkan bahwa nilai Deltat sesuai dengan kondisi fisis tersebut adalah berkisar antara T /20 − T /10. Pada rentang ∆t tersebut, hubungan nilai α dan rasio kebocoran f adalah linier. Dengan mengetahui informasi ni-

10| 8. OSILASI TEREDAM SISTEM PEGAS DENGAN MASSA BEBAN TAK KONSTAN 11000 10000 9000

Tstop (s)

8000 7000 6000 5000 4000 3000 2000 1000 1

2

3

4

5 f (x10e-6)

6

7

8

9


lai α ini, kita dapat menghitung besar lubang kebocoran yang dimiliki beban pegas hanya dengan menghitung waktu yang dibutuhkan pegas tersebut untuk berosilasi.

8.6

Referensi

1. Plastino A.R. dan Muzzio J.C., ”On the use and Abuse of Newtons Second Law for Variable Mass Problems”, Celest. Mech. Dyn. Astron (1991) 2. Rodrigues H. dan Panda N., ”The Motion of A Leaking Oscillator: A Study for the Physics Class”, Rio de Janeiro, RJ, Brazil (2014)

8.7

Lampiran

8.7.1

Kode C++

Terdapat 2 file kode c++ yang digunakan pada simulasi ini. Pertama, kode c++.cpp digunakan untuk mencatat posisi beban selama berosilasi (z(t)) posisi rata-rata (¯ z ) dan waktu simulasi (t). Data tersebut tercatat ke dalam file txt dengan tujuan untuk keperluan data analisis dan untuk mengetahui posisi beban pegas secara realtime. Berikut kodenya : #include #include <math.h>

11

8.7. LAMPIRAN 11000 10000 9000 8000 Tstop (s)

7000 6000 5000 4000 3000 2000 1000 0 1

2

3

4

5 f (x10e-6)

6

7

8

9


#include using namespace std; int main(int argc, char *argv[]) { //INPUT KONSTANTA double k=100; //konstanta pegas double b=1.8*pow(10,-2); //konstanta gaya gesek udara double g=9.81; //percepatan gravitasi

//PARAMETER BEBAN double mb=1; //massa wadah double mw0=10; //massa cairan double A=1; //luas penampang wadah double aw=0.005; //luas penampang kebocoran double h=0.05; //tinggi wadah double T=2*3.14*(pow((mb/k),(0.5))); //hitung periode osilasi double f=aw/A;

//rasio luas kebocoran terhapad luas penampang wadah

//PARAMETER WAKTU SIMULASI double tau=(1/f)*(pow((2*h/g),0.5)); double tstep=0.05; double tbeg=0; double tend=100;

//hitung waktu pengosongan

12| 8. OSILASI TEREDAM SISTEM PEGAS DENGAN MASSA BEBAN TAK KONSTAN 9e-05 8e-05 7e-05

Alpha (1/s)

6e-05 5e-05 4e-05 3e-05 2e-05 1e-05 0 1

2

3

4

5 f (x10E-5)

6

7

8

9

Gambar 8.8: Grafik hubungan α terhadap f dengan ∆t = T/20.

//INITIAL double v=0; //kecepatan awal; double z0=-(mb*g)/k; //posisi kesetimbangan double z=-((mw0+mb)*g)/k; //posisi awal simulasi double z_bar=0; double EK=0; double EP=0;

int n= double double double double double

(tend-tbeg)/tstep; mw, a; zi[n+1]; zi_bar[n+1]; EKi[n+1]; EPi[n+1];

ofstream file_; //buat file file_.open ("data_out.txt"); for (int i=0; i<=n; i++){ double t=i*tstep; bool con1= (t=tau);

//buka, lalu tulis file

13

8.7. LAMPIRAN if (con1) { mw=mw0*pow((1-(f*t*(pow((g/(2*h)),0.5)))),2); } else if (con2) { mw=0; //setelah air habis } a=-((k*(z-z0))/(mb+mw))-g-((b*v)/(mb+mw)); v=v+(a*tstep); //hitung v(t) z=z+(v*tstep); //hitung z(t) z_bar=-((mw+mb)*g)/k; //hitung z_bar (t)

//ketika masih ada air

//hitung a(t)

//BUAT ARRAY z(t) dan z_bar (t) zi[i]=z-z0; zi_bar[i]=z_bar; //TAMPILKAN t, z(t) dan z_bar (t) cout<< t << "\t"; cout<< zi[i] << "\t"; cout<< zi_bar[i] << "\n"; //SIMPAN file_<< file_<< file_<< file_<< file_<<

t, z(t), dan z_bar (t) t << "\t"; zi[i] << "\t"; zi_bar[i] << "\t"; EKi[i] << "\t"; EPi[i] << "\n";

} cout<<"\n"<<"t(sec)"<< "\t"<<"z(t)"<< "\t"<<"z_bar(t)"<<"\n"; cout<<"\n"<<"tau :"<
<math.h>

using namespace std; int main(int argc, char *argv[]) { //INPUT KONSTANTA double k=100;

14| 8. OSILASI TEREDAM SISTEM PEGAS DENGAN MASSA BEBAN TAK KONSTAN double b=1.8*pow(10,-1); double g=9.81; double Cstop=pow(10,-10);;//pow(10,-1000); //PARAMETER BEBAN double mb=1; double mw0=10; double A=1; double aw=0.0009; double h=0.05; double T=2*3.14*(pow((mb/k),(0.5))); //PARAMETER WAKTU SIMULASI double tstep=T/40; double tbeg=0; double tend=10000000;

double f=aw/A; //rasio luas kebocoran terhapad luas penampang wadah double tau=(1/f)*(pow((2*h/g),0.5)); //hitung waktu pengosongan //INITIAL double v=0; //kecepatan awal; double z0=-(mb*g)/k; //posisi kesetimbangan double z=-((mw0+mb)*g)/k; //posisi awal simulasi

int n= (tend-tbeg)/tstep; double mw, a, DZ;

//INPUT KONSTANTA double Z[2]; double dz[1]; double Tstop[1]; double I[1]; for (int i=0; i<=n; i++){ double t=i*tstep; bool con1= (t=tau); if (con1) { mw=mw0*pow((1-(f*t*(pow((g/(2*h)),0.5)))),2); } else if (con2) { mw=0; //setelah air habis }

//ketika masih ada air

15

8.7. LAMPIRAN a=-((k*(z-z0))/(mb+mw))-g-((b*v)/(mb+mw)); v=v+(a*tstep); //hitung v(t) z=z+(v*tstep); //hitung z(t)

//hitung a(t)

bool con3= (i % 2==0); bool con4= (i % 2==1); if (con3) { Z[2]=z-z0; } else if (con4) { Z[1]=z-z0; } DZ=abs((abs(Z[2]))-(abs(Z[1]))); dz[1]=DZ; Tstop[1]=t; I[1]=i; bool con5=(dz[1] data_out.txt

8.7.2

Script Gnuplot

set terminal postscript eps enhanced color font "Helvetica,20"

16| 8. OSILASI TEREDAM SISTEM PEGAS DENGAN MASSA BEBAN TAK KONSTAN set xlabel ’Waktu (s)’ set ylabel ’Posisi (m)’ set output ’fig03.eps’ set key right bottom plot ’Data3.txt’ u 1:2 w l title ’Z(t)’,’Data3.txt’ u 1:3 w l lt 1 lc 3 title ’Z(t)rata-rata’ Script Gnuplot di atas dipanggil dengan cara gnuplot script3.gps

|9

Orientasi Gerak Jatuh Bebas Sistem 3 Partikel Jesi Pebralia | [email protected] Yunita Citra Dewi | [email protected] Donny Dwiputra | [email protected] Rouf | [email protected] Partikel adalah objek yang memiliki massa, posisi dan kecepatan serta dapat dipengaruhi oleh gaya-gaya. Partikel dibedakan menjadi dua jenis berdasarkan struktur pembentuknya yaitu partikel rigid dan partikel non-rigid. Pada penelitian ini, sistem partikel yang ditinjau adalah sistem partikel non-rigid. Sistem partikel non-rigid adalah sekumpulan partikel-partikel yang membentuk suatu struktur tertentu dimana masing-masing partikel penyusun sistem dihubungkan melalui pegas dengan konstanta pegas tertentu (Kilian, 2005). Adapun struktur partikel dapat dibedakan menjadi dua jenis yaitu regular particles dan irregular particles (Webb, 2002. Regular particles adalah partikel-partikel yang mempunyai susunan kisi-kisi teratur seperti kubus, segitiga, tetrahedral, dan lain sebagainya. Sedangkan, irregular particles adalah partikel-partikel yang mempunyai susunan kisi-kisi yang tidak teratur. Pemodelan pada benda pejal telah sebelumnya dilakukan oleh, contohnya, (Baraff, 1995), (Hahn, 1988), (Baber, 2006) dengan beberapa variasi algoritma yang pada umumnya bertumpu pada perhitungan tensor momen inersia yang cukup kompleks. Pada tulisan ini kami mencoba melakukan simulasi sederhana untuk memodelkan benda tegar (sistem partikel), dengan mengaproksimasi sistem benda pejal menjadi sistem tiga partikel terhubung dengan pegas berkonstanta pegas besar. Dalam penelitian ini, struktur partikel yang digambarkan adalah regular particles dengan sistem partikel berbentuk segitiga. Sistem partikel berbentuk segitiga digambarkan oleh tiga buah partikel yang masing-masing dihubungkan 1

2

| 9. ORIENTASI GERAK JATUH BEBAS SISTEM 3 PARTIKEL

dengan suatu pegas. Sistem partikel dengan masing-masing partikel mempunyai massa m mengalami jatuh bebas dan dipengaruhi oleh gaya hambat F . Orientasi gerak jatuh bebas dimodelkan melalui sistem tiga partikel dengan parameter berupa massa masing-masing partikel dan gaya hambat. Penyelesaian persamaan diferensial dilakukan secara numerik dengan menggunakan Metode Euler.

9.1

Teori

Gerak jatuh bebas adalah gerak sebuah benda yang dijatuhkan tanpa kecepatan awal. Dinamika dari gerak jatuh bebas dijelaskan oleh Hukum Newton kedua F~ =

X

mi~ai ,

(9.1)

i

di mana m adalah massa benda dan a adalah percepatan benda. Persamaan (9.1) menggambarkan resultan gaya-gaya yang bekerja terhadap suatu benda dan hubungannya dengan percepatan benda.

Gambar 9.1: Diagram gaya-gaya yang bekerja pada sistem tiga partikel.

3

9.2. BATASAN MODEL

Adapun gaya-gaya yang bekerja pada sistem partikel meliputi gaya gravitasi, gaya gesek udara dan gaya pegas. Gaya gravitasi adalah gaya tarik bumi yang dinyatakan oleh persamaan F~G = mg yˆ

(9.2)

Gaya gravitasi disebut juga gaya berat, dengan m menyatakan massa partikel dan ~g adalah besar percepatan gravitasi. Gaya lain yang mempengaruhi gerak sistem partikel adalah gaya hambat udara. Gaya hambat udara adalah gaya yang disebabkan oleh adanya gesekan antara benda yang jatuh dengan udara. Gaya hambat udara dinyatakan oleh persamaan berikut F~f i = −bi~vi

(9.3)

dengan bi merupakan koefisien hambatan tiap partikel. Pada sistem tiga partikel, masing-masing partikel mempunyai gaya hambat udara tertentu. Selanjutnya, gaya lain yang mempengaruhi sistem partikel adalah gaya pegas. Gaya pegas dinyatakan oleh persamaan berikut F~si = −ki~ri

(9.4)

dengan i k merupakan konstanta pegas. Adapun sifat dari pegas yang digunakan dalam sistem partikel non-rigid ini adalah pegas yang memiliki konstanta pegas yang sangat besar sehingga akan meminimalkan atau bahkan meniadakan efek tegangan dan regangan yang muncul. Adapun resultan gaya-gaya yang mempengaruhi orientasi gerak dari sistem tiga partikel adalah penjumlahan dari ketiga gaya di atas. Sehingga diperoleh F~i = F~G + F~f i + F~si

(9.5)

Substitusi persamaan (9.5) ke persamaan (9.1) akan memberikan

~ai = ~g −

bi ki ~vi − r~i , mi mi

(9.6)

yang merupakan percepatan partikel setiap saat.

9.2

Batasan model

Orientasi gerak jatuh bebas sistem tiga partikel dimodelkan dalam gerak dua dimensi yaitu dalam sumbu horizontal x dan sumbu vertikal y. Kecepatan sistem partikel hanya memiliki komponen pada arah x dan y

4

| 9. ORIENTASI GERAK JATUH BEBAS SISTEM 3 PARTIKEL ~v = vx eˆx + vy eˆy .

9.3

(9.7)


Dengan menggunakan metode Euler Persamaan (9.6) dapat memberikan kecepatan dan posisi setiap saat untuk arah x, untuk tiap partikel-i,

b k vx (t) − x(t) m m vx (t + ∆t) = vx (t) + ax (t)∆t rx (t + ∆t) = rx (t) + vx (t)∆t

(9.8)

ax (t) = −

(9.9) (9.10)

dan untuk arah y

b k vy (t) − y(t) m m vy (t + ∆t) = vy (t) + ay (t)∆t

(9.12)

ry (t + ∆t) = ry (t) + vy (t)∆t

(9.13)

ay (t) = gy −

(9.11)

Dalam implementasinya Persamaan-persamaan (9.8) - (9.13) dituangkan dalam algorima berikut ini: L01. Inisiasi nilai awal m[3], k[3], l[3], b[3], vx [3], vy [3], x[3], y[3], Fx [3], Fy [3],

a1 , a2 , b1 , b2 , xpm , ypm ? L02. ∆t, tbeg , tend ? L03. t = tbeg L04. Iterasi partikel i = 0, i < 3, i + + b[i] L05. Percepatan sebelum memantul, ax [i](t) = − m[i] vx [i] − g m[i]

b[i] m[i]

k[i] m[i]

x[i](t),

k[i] m[i]

ay [i](t) = − vx [i] − y[i](t) L06. Kecepatan sebelum memantul, vx [i](t + ∆t) = vx [i](t) + ax [i](t)∆t, vy [i](t + ∆t) = vy [i](t) + ax [i](t)∆t L07. Syarat pantul : sumbu x, x[i] < a1 atau x[i] > a2 , sumbu y, y[i] < b1 atau y[i] > b2 Fx [i] b + m[i] vx [i], ay [i](t) = L08. Percepatan setelah memantul, ax [i](t) = − m[i]

b x [i] − m[i] + m[i] vx [i] + gdt L09. Kecepatan setelah memantul, vx [i](t + ∆t) = −vx [i](t) + ax [i](t)∆t, vy [i](t + ∆t) = −vy [i](t) + ax [i](t)∆t L10. x[i](t + ∆t) = x[i](t) + vx [i](t)∆t L11. y[i](t + ∆t) = y[i](t) + vy [i](t)∆t x[i](t)m[i] L12. xp m(t + ∆t) = xp m(t) + M y[i](t)m[i] L13. yp m(t + ∆t) = yp m(t) + M

5

9.4. HASIL DAN DISKUSI L14. L15. L16. L17. L18. L19. L20. L21.

9.4

Energi kinetik Ek(t) = Ek(t) + 0.5m[i]v[i]2 Energi potensial Ep(t) = Ep(t) + m[i]gy[i](t) + 0.5k[i]x[i]2 Energi mekanik total E(t) = Ek(t) + Ep(t) Definisikan dan hitung parameter orientasi yp i[i] = ypm − y[i] Time increment t = t + ∆t t ≤ tend → L03 Output x[3](t), y[3](t), xp m, yp m, Ek , Ep , E, t Export output ke springs.txt, springs-energy.txt, springs-CM-Frame.txt.

Hasil dan diskusi

Simulasi ini menggunakan nilai parameter yang pada dasarnya scale invariant, tetapi masuk akal jika kita interpretasikan satuannya dalam SI. Nilai-nilai parameter yang digunakan dalam simulasi tercantum dalam Tabel 10.1. m[3] k[3] l g b[3] v0x [3], v0y [3] x0 [3] y0 [3]

{100, 100, 100} {1, 1, 1} × 108 , {1, 1.2, 1.3} × 108 1 9.8 {1, α, β} × 0.01225 {0, 0, 0} {2, 3, 2.5} √ {14, 14, 14 + 0.75}

Tabel 9.1: Nilai-nilai parameter simulasi. Dari data pada Tabel 10.1 nampak bahwa susunan partikel adalah m1 dan m2 sejajar dan terletak di bawah m3 . Masing-masing partikel mengalami gaya gesek dengan koefisien gesek yang bervariasi: {1, α, β} × 0.01225 (0.01225 = g/800) berturut turut untuk m1 , m2 , dan m3 . Simulasi ini dilakukan pada dua kasus, kasus pertama untuk konstanta pegas k seragam (108 ) dan kasus kedua untuk k berbeda (c.f. Tabel 10.1). Untuk mengetahui seberapa besar orientasi sistem berubah seiring waktu, definisikan parameter orientasi berupa simpangan vertikal yang terjadi antara posisi pusat massa dan posisi partikel ke-i, ypii ≡ ypm − yi .

(9.14)

Parameter ini menggambarkan seberapa jauh posisi vertikal partikel berubah dari waktu ke waktu. Jika pada awal nilainya positif dan pada suatu keadaan akhir nilainya negatif maka sistem pasti telah mengalami rotasi. Hal ini tidak ditinjau lebih jauh karena dibutuhkan parameter kuantitatif lain untuk mendeskripsikan sifat rotasi ini. Di sini akan ditinjau seberapa besar ”fluktuasi” orientasi partikel yang dialaminya saat jatuh bebas. Pada Tabel 9.4 akan nampak nilai standar deviasi dari parameter ypii , yang dilabelkan σi , untuk nilai k seragam. Sementara Tabel 9.4 menampakkan nilai σi untuk k berbeda-beda. Standar deviasi ini dihitung dari t = 0 sampai t = 10.3, nilai ini diambil karena pada saat ini partikel sudah dalam keadaan stabil.

6


Tabel 9.2: Standar deviasi dari parameter gesekan tiap benda dan k seragam. α β σ1 σ2 1 1 0.000923 0.000922 1 4 0.000466 0.0004666 4 1 0.258379 0.312573 4 4 0.122643 0.253006

orientasi ypi [i] untuk variasi nilai σ3 0.001845 0.000933 0.36744 0.136158

σtotal 0.003690 0.001866 0.938382 0.51181

Tabel 9.3: Standar deviasi dari parameter orientasi ypi [i] untuk variasi nilai gesekan tiap benda dan k = {1, 1.2, 1.3} × 108 . α β σ1 σ2 σ3 σtotal 1 1 0.090436 0.436939 0.426239 0.953615 1 4 0.068643 0.053308 0.01634 0.138288 4 1 0.401586 0.296273 0.309193 1.007054 4 4 0.122057 0.252670 0.136205 0.510932

Dari kedua tabel di atas nampak untuk α 6= 1 nilai standar deviasi total lebih besar dibandingkan dengan yang lain, hal ini terjadi karena nilai α menentukan gerak m1 dan m2 terus sejajar sebelum tumbukan, dengan melihat simetri cermin dari sistem ini. Nilai α yang tidak satu akan membuat m1 dan m2 berfluktuasi pada kondisi jatuh bebas, sehingga σ membesar. Untuk mengamati lebih jauh perilaku perubahan orientasi sistem, amati Gambar 9.4. Gambar ini memberikan informasi tentang dinamika ypii untuk ketiga partikel terhadap waktu.

9.5

Ringkasan

Berdasarkan pemodelan yang telah dibuat, telah diketahui bahwa parameter yang mempengaruhi orientasi jatuh dari sistem tiga partikel adalah koefisien hambatan udara.

9.6

Referensi

1. Baber, R., ”Rigid body simulation”, University of Warwick, (2006). 2. Baraff, D., ”Rigid body simulation”, SIGGRAPH Course Notes 1992, Vol.19, (1995). 3. Hahn, James K., ”Realistic animation of rigid bodies”, ACM SIGGRAPH Computer Graphics, Vol.22.4, (1988), p.299-308. 4. Kilian, A. and Ochsendorf, J., ”Particle-spring systems for structural form

7

9.6. REFERENSI

1

1 Vertical shift particle-1 Vertical shift particle-2 Vertical shift particle-3

Particle-i vertical shift from vertical CM (m)


Vertical shift particle-1 Vertical shift particle-2 Vertical shift particle-3

0.5

0

-0.5

-1

0.5

0

-0.5

-1 0

2

4

6

8

10

0

2

4

Time (s)

6




0.5

0

-0.5

-1

0.5

0

-0.5

-1 0

2

4

6

8

10

0

2

4

Time (s)

6

8

10

Time (s)

1





10


0.5

0

-0.5

-1

0.5

0

-0.5

-1 0

2

4

6

8

10

0

2

4

Time (s)

6

8

10

Time (s)

1





8

Time (s)

1

0.5

0

-0.5

-1

0.5

0

-0.5

-1 0

2

4

6 Time (s)

8

10

0

2

4

6

8

10

Time (s)

Gambar 9.2: Perubahan parameter orientasi ypii (vertical shift) terhadap waktu sampai sistem stabil. Kolom kiri untuk k seragam dan kanan k bervariasi (c.f. Tabel 10.1); baris 1–4: (α, β) = (1, 1), (1, 4), (4, 1), (4, 4).

8

| 9. ORIENTASI GERAK JATUH BEBAS SISTEM 3 PARTIKEL finding”, Journal-International Association for Shell and Spatial Structures, Vol.148, (2005), p.77. 5. Webb, P., ”Interpretation of particle size reported by different analytical techniques”, Homepage Micromeritics, (2002).

9.7 9.7.1

Lampiran Kode C++ (standalone)

#include #include #include <math.h> using namespace std; int main() { double dt=0.0002, t=0; //besar time step (resolusi waktu) int a1=0, a2=10, b1=0, b2=15; //batas-batas ruang (kotak) int m[3]={100,100,100}; //massa benda int k[3]={100000000,120000000,130000000}; //konstanta semua pegas double l[3]={1,1,1}; //panjang awal semua pegas double g=9.8; //percepatan gravitasi double b[3]={0.01225,0.01225,0.01225}; //konstanta gesekan udara //menulis file output ofstream fout; fout.open("3 Particles with springs.txt"); ofstream fout_energy; fout_energy.open("3 Particles with springs energy.txt"); ofstream fout_CM; fout_CM.open("3 Particles with springs CM frame.txt"); ofstream fout_ypi; fout_ypi.open("3 Particles with springs orientation.txt"); double vx[3]={0,0,0}, vy[3]={0,0,0}; //kecepatan (vx,vy) double x[3]={2,3,2.5}, y[3]={12+2,12+2,12+2+sqrt(0.75)};// posisi awal double Fx[3]={0,0,0}, Fy[3]={0,0,0}; //gaya pegas //iterasi selama semenit while(t<=50){ double jarak[3]={sqrt(pow(x[1]-x[0],2)+pow(y[1]-y[0],2)), //panjang pegas tiap waktu sqrt(pow(x[2]-x[0],2)+pow(y[2]-y[0],2)), sqrt(pow(x[2]-x[1],2)+pow(y[2]-y[1],2))}; //gaya pegas yang bekerja pada benda 1 Fx[0]=-k[0]*(jarak[0]-l[0])*(x[0]-x[1])/jarak[0]k[1]*(jarak[1]-l[1])*(x[0]-x[2])/jarak[1]; Fy[0]=-k[0]*(jarak[0]-l[0])*(y[0]-y[1])/jarak[0]k[1]*(jarak[1]-l[1])*(y[0]-y[2])/jarak[1]; //gaya pegas pada benda 2

9.7. LAMPIRAN Fx[1]=-k[0]*(jarak[0]-l[0])*(x[1]-x[0])/jarak[0]k[2]*(jarak[2]-l[2])*(x[1]-x[2])/jarak[2]; Fy[1]=-k[0]*(jarak[0]-l[0])*(y[1]-y[0])/jarak[0]k[2]*(jarak[2]-l[2])*(y[1]-y[2])/jarak[2]; //gaya pegas yang bekerja pada benda 3 Fx[2]=-k[1]*(jarak[1]-l[1])*(x[2]-x[0])/jarak[1]k[2]*(jarak[2]-l[2])*(x[2]-x[1])/jarak[2]; Fy[2]=-k[1]*(jarak[1]-l[1])*(y[2]-y[0])/jarak[1]k[2]*(jarak[2]-l[2])*(y[2]-y[1])/jarak[2]; double M_tot=m[0]+m[1]+m[2]; //massa total double x_pm=0,y_pm=0; // pusat massa (nantinya) //hitung kecepatan terus posisi partikel for(int i=0;i<3;i++){ //syarat pantul if (x[i]a2) { vx[i]=-vx[i]-Fx[i]*dt/m[i] + b[i]*vx[i]/m[i];} if (y[i]b2) { vy[i]=-vy[i]-Fy[i]*dt/m[i]+g*dt + b[i]*vy[i]/m[i];} //waktu gak mantul vx[i]=vx[i]+Fx[i]*dt/m[i]-b[i]*vx[i]/m[i]; vy[i]=vy[i]+Fy[i]*dt/m[i]-g*dt-b[i]*vy[i]/m[i]; x[i]=x[i]+vx[i]*dt; y[i]=y[i]+vy[i]*dt; //hitung pusat massa x_pm=x_pm+m[i]*x[i]; y_pm=y_pm+m[i]*y[i]; //output dan tulis file //cout << jarak[i] << "\t"; //cout << x[i] << "\t" << y[i] << "\t"; fout << x[i] << "\t" << y[i] << "\t"; } t=t+dt; //cout << x_pm/M_tot << "\t" << y_pm/M_tot; //plot pusat massa fout << x_pm/M_tot << "\t" << y_pm/M_tot; //cout << endl; //pindah baris tiap iterasi waktu fout << endl; //plot dalam koordinat pusat massa double x_CMF[3], y_CMF[3]; for(int i=0;i<3;i++){ x_CMF[i]= x[i]-x_pm/M_tot; y_CMF[i]= y[i]-y_pm/M_tot; fout_CM << x_CMF[i] << "\t" << y_CMF[i] << "\t"; } fout_CM << 0 << "\t" << 0 << endl;

9

10

| 9. ORIENTASI GERAK JATUH BEBAS SISTEM 3 PARTIKEL fout_ypi << t << "\t"; //baru bagian waktunya utk file orientasi

double Ek=0, Ep=0, Etot=0; double y_pi[3]; //parameter perubahan orientasi for(int i=0;i<3;i++){ Ek=Ek+0.5*m[i]*(pow(vx[i],2)+pow(vy[i],2)); Ep=Ep+m[i]*g*y[i]+0.5*k[i]*pow(jarak[i]-l[i],2); /*Penting! parameter yang menjelaskan perilaku perubahan orientasi sistem (nebeng iterasi di sini)*/ y_pi[i]=y_pm/M_tot - y[i]; fout_ypi << y_pi[i] << "\t"; } fout_ypi << endl; //menghitung energi potensial, kinetik, total sistem Etot=Ep+Ek; //cout << t << "\t" << Ek << "\t" << Ep << "\t" << Etot << endl; fout_energy << t << "\t" << Ek << "\t" << Ep << "\t" << Etot << endl; } system("pause"); return 0; }

Kode C++ di atas dikompilasi dengan cara g++ springs.cpp -o springs dan dipanggil dengan cara ./springs > springs.txt, springs_energy.txt, springs_CM_Frame.txt

9.7.2 set set set set set q = d =

Script Gnuplot untuk Plot Lintasan 3 Partikel

terminal gif size 400,400 animate delay 0.05; output ’springs.gif’; xrange[-0.1:4.1]; yrange[-0.1:5.1]; key off; 1; 500;

do for [t=0:200]{ plot ’./springs.txt’ using 1:2 every ::d*t::d*t w lp pt 7 ps 3 \ ,’’ using 3:4 every ::d*t::d*t w lp pt 7 ps 2 \

11

9.7. LAMPIRAN ,’’ using 5:6 every ::d*t::d*t w lp pt 7 ps 3.2 \ ,’’ using 7:8 every ::d*t::d*t w lp pt 6 ps 2 \ } Dalam kerangka pusat massa, set set set set set q = d =

terminal gif size 400,400 animate delay 0.05; output springs_CM_frame.gif’; xrange[-1.2:1.2]; yrange[-1.2:1.2]; key off; 1; 500;

do for [t=0:400]{ plot ’./springs_CM_frame.txt’ using ,’’ using 3:4 every ::d*t::d*t w lp ,’’ using 5:6 every ::d*t::d*t w lp ,’’ using 7:8 every ::d*t::d*t w lp } }

1:2 every ::d*t::d*t w lp pt 7 ps 3 \ pt 7 ps 4 \ pt 7 ps 6.4 \ pt 6 ps 4 \

Script di atas dieksekusi dengan menggunakan Gnuplot.

12


| 10

Inverted Pendulum Pada Sebuah Kereta Hardi Hamzah | [email protected] Firmansyah | [email protected] M. Fitrah Alfian R. S. | [email protected] M. Bisri Mutofa | [email protected] Terdapat suatu sistem pendulum terbalik (inverted pendulum) yang diletakkan pada sebuah kereta yang dapat bergerak bolak - balik secara horizontal. Kestabilan dari sistem pendulum ini hanya dapat terjadi ketika tidak ada input gaya yang bekerja pada sistem. Pada penelitian ini kami mengkaji pengaruh massa kereta mc , massa pendulum mp , sudut simpangan awal θ, panjang batang pendulum tak bermassa l , dan gaya penyeimbang sebagai fungsi sudut simpangan F (θ) = kθ terhadap waktu yang dibutuhkan oleh sistem ini hingga berada pada posisi setimbangnya. Solusi gerak diperoleh dengan menggunakan metode Euler.

10.1

Teori

Sistem pendulum terbalik (inverted pendulum) yang diletakkan pada sebuah kereta, di mana batang pendulum dapat bebas berputar / berotasi pada titik poros yang terletak pada kereta yang dapat bergerak bolak - balik secara horizontal, seperti yang ditunjukkan pada Gambar 10.1. dengan l adalah panjang pendulum, θ adalah sudut simpangan awal pendulum, k konstanta dari gaya penyeimbang fungsi sudut, mc adalah massa kereta, dan mp adalah massa pendulum. 1

2

| 10. INVERTED PENDULUM PADA SEBUAH KERETA

Gambar 10.1: Diagram gaya dengan parameter-parameter yang digunakan: l, θ, k, mc , dan mp . Untuk gerak dalam arah harizontal, berlaku hukum Newton kedua, yaitu: ΣF~ = m~a,

(10.1)

F~px = −mp acx + f~ux .

(10.2)

dapat ditulis

Kita asumsikan gaya gesek udara memiliki nilai: dx~p . f~u = b dt dan percepatan keretanya: |acx | =

k θ. mc

(10.3)

(10.4)

di mana θ adalah sudut yang dibentuk pendulum terhadap pusat kesetimbangannya di sumbu y. Kemudian pada arah y, yaitu: F~py = mp~g + f~uy .

10.2

(10.5)

Syarat batas

Gaya luar yang bekerja pada kereta dibuat sebagai fungsi sudut θ dan hanya pada arah x: F~ext = Fx (θ)ˆ ex = kθˆ ex .

(10.6)

3

10.3. SOLUSI NUMERIK DAN ALGORITMA Gaya gesek pada lantai diasumsikan tidak ada: f~g = 0.

(10.7)

Kecepatan pendulum dibatasi hanya memiliki komponen pada arah x dan y v~p = vpx eˆx + vpy eˆy .

(10.8)

Kecepatan kereta dibatasi hanya memiliki komponen pada arah x: v~c = vcx eˆx .

10.3

(10.9)


Sebelumnya kami menggunakan pendekatan Lagrange untuk memperoleh solusi dinamika dari sistem inverted pendulum (Castro, 2012), dimana algoritma yang digunakan merupakan algoritma yang dikembangkan oleh Andreou (2013). Namun, dikarenakan algoritma tersebut terdapat suatu fungsi yang hanya dapat digunakan pada program Matlab yaitu fungsi (acker) sehingga kami beralih menggunakan metode Euler. Dengan menggunakan metode Euler, persamaan (10.2) - (10.5) dapat memberikan kecepatan dan posisi setiap saat untuk kereta, yaitu:

k x(t) , arcsin mc l vcx (t + ∆t) = vcx (t) + acx (t)∆t, xc (t + ∆t) = xc (t) + vcx (t)∆t, acx (t) =

(10.10) (10.11) (10.12)

dan untuk pendulum pada arah x, yaitu

b vpx (t), mp vpx (t + ∆t) = vpx (t) + apx (t)∆t, xp (t + ∆t) = xp (t) + vpx (t)∆t, apx (t) = −acx (t) +

(10.13) (10.14) (10.15)

serta untuk pendulum pada arah y, yaitu

b vpy (t), mp vpy (t + ∆t) = vpy (t) + apy (t)∆t, apy (t) = gy +

yp (t + ∆t) = yp (t) + vpy (t)∆t.

(10.16) (10.17) (10.18)

4


Dalam implementasinya persamaan-persamaan (10.10) - (10.18) dituangkan dalam algoritma berikut ini: L01. L02. L03. L04.

∆t, tbeg , tend ? l, k, mp , mc , b, g? xp0 , yp0 , vpx0 , vpy0 , xc0 , vcx0 ? t = tbeg

L05. acx (t) =

k mc

arcsin

L06. apx (t) = −acx (t) + b mp

vpy (t) vcx (t + ∆t) = vcx (t) + acx (t)∆t vpx (t + ∆t) = vpx (t) + apx (t)∆t vpy (t + ∆t) = vpy (t) + apy (t)∆t xc (t + ∆t) = xc (t) + vcx (t)∆t xp (t + ∆t) = xp (t) + vpx (t)∆t yp (t + ∆t) = yp (t) + vpy (t)∆t

L07. apy (t) = gy +

L08. L09. L10. L11. L12. L13.

L14. Jika arcsin

x(t) l

L15. acx (t) = − mkc

<0 arcsin x(t) l

L16. apx (t) = acx (t) −

b mp

vpx (t)

b mp

vpy (t) vcx (t + ∆t) = −vcx (t) + acx (t)∆t vpx (t + ∆t) = −vpx (t) + apx (t)∆t vpy (t + ∆t) = vpy (t) + apy (t)∆t xc (t + ∆t) = −xc (t) + vcx (t)∆t xp (t + ∆t) = −xp (t) + vpx (t)∆t yp (t + ∆t) = yp (t) + vpy (t)∆t t = t + ∆t Jika t < tend → L05

L17. apy (t) = gy + L18. L19. L20. L21. L22. L23. L24. L25.

x(t) l b mp vpx (t)

10.4

Hasil dan diskusi

Nilai-nilai parameter yang digunakan dalam simulasi tercantum dalam Tabel 10.1. Kemudian setelah itu kami menerapkan variasi parameter-parameter tersebut pada program dengan hasil sebagai berikut: • Keadaan awal acuan

Untuk membandingkan hasil dari variasi parameter, kami menjadikan salah satu grup parameter sebagai acuannya. Acuan parameter tersebut adalah l=1 m, θ=0.1 rad, k=1 kg.m2 /s2 , mc =0.5 kg, mp =0.1 kg, dan b=10−5 kg/s. Kemudian dari parameter-parameter acuan ini kami dapatkan bahwa pendulum menjadi setimbang di sekitar t=16 s. Simulasi gerak dari inverted pendulum ini ditunjukkan pada Gambar dan .

5


Tabel 10.1: Nilai-nilai parameter simulasi. Parameter tbeg tend tstep l θ k mc mp b g

Nilai 0 20 10− 5 0.5 dan 1 0.1, 1, dan 3 1 dan 5 0.1 dan 0.5 0.1 dan 0.5 0 dan 10−5 9.81

Satuan s s s m rad kg.m2 .s− 2 kg kg kg.s− 1 m.s− 2

Gambar 10.2: Simulasi gerak inverted pendulum sebelum setimbang. • Pengaruh massa pendulum dan massa kereta

Jika massa kereta dan pendulum dibuat sama sebesar 0.1 kg, pendulum masih dapat mencapai kesetimbangan di atas 20 s. Namun, jika massa pendulum lebih besar dari massa kereta, dengan massa pendulum 0, 5 kg sedangkan massa kereta 0, 1 kg, pendulum tidak mencapai kesetimbangan, kecuali konstanta gaya diperbesar. Apabila massa kereta yang lebih besar dibandingkan dengan massa pendulum, maka pendulum dapat lebih cepat menuju kesetimbangan.

• Pengaruh nilai konstanta gaya k

Apabila nilai konstanta gaya luar diperbesar dari 1 kg.m2 /s2 menjadi 5 kg.m2 /s2 , pada detik ke-14 pendulum telah mulai terlihat setimbang. Sehingga semakin besar konstanta gayanya, maka pendulum semakin cepat setimbang.

• Pengaruh panjang pendulum

Apabila panjang pendulum diperpendek menjadi 0, 5 m, maka pendulum semakin cepat mencapai kesetimbangan. Sebaliknya, semakin panjang

6


Gambar 10.3: Simulasi gerak inverted pendulum setelah setimbang. pendulum, maka semakin lama mencapai kesetimbangan. Hal ini sesuai dengan prinsip torsi, semakin besar lengan gaya, maka torsi semakin besar yang menyebabkan benda semakin tidak setimbang. • Pearuh besar sudut simpangan

Semakin besar sudut simpangan θ, maka akan semakin lama pendulum setimbang. Namun, ada suatu nilai kritis sudut di mana ketika sudutnya melewati sudut kritis tersebut, maka pendulum pasti akan terjatuh. Kami belum mencari sudut kritis tersebut karena variasi parameter sudut sebesar 1 rad dan 3 rad langsung membuat pendulum terjatuh dan tidak bisa setimbang di posisi atas lagi.

• Pengaruh adanya hambatan udara

Diberikan gaya hambat udara dengan kontanta hambat sebesar 10−5 kg/s dan tanpa adanya hambatan. Dari percobaan tersebut, terlihat bahwa pendulum semakin lama menuju ke kesetimbangan jika ditambahkan hambatan udara karena berlawanan dengan gaya dari kereta.

10.5

Ringkasan

Telah diperoleh bahwa besar nilai setiap parameter (yang divariasikan), yaitu l, θ, k, mp , mc , dan b akan mempengaruhi lama waktu pendulum untuk mencapai kesetimbangan. Bahkan apabila beberapa nilai parameter terlalu besar misalnya θ, maka pendulum akan jatuh dan tidak kembali setimbang di atas kereta. Namun, perlu dikaji kembali pengaruh lain yang tidak kami variasikan, yaitu percepatan gravitasi g dan tstep . Selain itu juga, perlu dikaji kembali beberapa nilai kritis dari sistem inverted pendulum ini misalnya θ, mc , mp , dan lain-lain. Perlu juga adanya percobaan dengan penambahan parameter lain seperti gaya gesek dengan lantai dan kasus batang bermassa.

10.6. REFERENSI

10.6

7

Referensi

1. Andreou, S, ”Revisiting The Inverted Pendulum On A Moving Cart”, European International Journal of Science and Technology, Vol.2(7), (2013). 2. Castro, A, ”Modelling and Dynamics Analysis of A Two-Wheeled Inverted Pendulum”, Thesis of Master of Science at Georgia Institute of Technology: August, (2012). 3. David Halliday, Robert Resnick, and Jearl Walker, Fundamentals of Physics, John Wiley & Sons (Asia), 8th, Extended, Student, Edition, (2008)

10.7

Lampiran

10.7.1

Kode C++ (Termasuk Kode Gnuplot di Dalamnya)

#include <sys/types.h> #include <sys/stat.h> #include #include #include #include #include #include #include

<sstream> "rwparams.h" "vect2.h"

using namespace std; int main(int argc, char *argv[]) { // Get compiled name const char *arg0 = argv[0]; const char *pname = strstr(arg0, "./"); if(pname != NULL) pname += 2; // Check number of argument if(argc < 2) { cout << "Usage: " << pname << " [pfile]"; cout << endl; cout << "pfile\tparameters file, e.g. params.txt"; cout << endl; exit(1); } // Get filename of pfile

8


const char *pfn = argv[1]; // Read parameters from pfile double mp = 0.0; readparam(pfn, "PMASS", mp); double mc = 0.0; readparam(pfn, "CMASS", mc); double l = 0.0; readparam(pfn, "LENGTH", l); double k = 0.0; readparam(pfn, "ACONST", k); double ani = 0.0; readparam(pfn, "ANINIT", ani); double g = 0.0; readparam(pfn, "GRAVACC", g); double eta = 0.0; readparam(pfn, "VISCOSA", eta); double tbeg = 0.0; readparam(pfn, "TBEG", tbeg); double tend = 0.0; readparam(pfn, "TEND", tend); double dt = 0.0; readparam(pfn, "TSTEP", dt); double Tdata = 0.0; readparam(pfn, "TDATA", Tdata); string ofnpre = ""; readparam(pfn, "OUTPRE", ofnpre); int digit = 0; readparam(pfn, "DIGIT", digit); string odir = ""; readparam(pfn, "OUTDIR", odir); // Verbose values cout << "PMASS\t" << mp << endl; cout << "CMASS\t" << mc << endl; cout << "LENGTH\t" << l << endl; cout << "ACONST\t" << k << endl; cout << "ANINIT\t" << ani << endl; cout << "GRAVACC\t" << g << endl; cout << "VICOSA\t" << eta << endl; cout << "TBEG\t" << tbeg << endl; cout << "TEND\t" << tend << endl; cout << "TSTEP\t" << dt << endl; cout << "TDATA\t" << Tdata << endl; cout << "OUTDIR\t" << odir << endl; cout << "OUTPRE\t" << ofnpre << endl; cout << "DIGIT\t" << digit << endl; // Check whether OUTDIR exists struct stat info;

10.7. LAMPIRAN

9

const char *odirc = odir.c_str(); bool OUTDIR_NOT_EXIST = stat(odirc, &info); if(OUTDIR_NOT_EXIST) { cout << endl << "Error: "; cout << "Output directory does not exist: " << odirc; cout << endl; cout << "Suggestion: Create the directory manually "; cout << "and run " << pname << " again" << endl; exit(2); } // Define initial position and velocity int N = 2; vect2 r[N], v[N]; r[0] = vect2(l*sin(ani), l*cos(ani)); r[1] = vect2(0, 0); v[0] = vect2(0, 0); v[1] = vect2(0, 0); // Construct output animation filename string oafnt = "animation.gnuplot"; string oafns = odir + "/" + oafnt; cout << "Writing " << oafns.c_str() << "\n"; // Iterate int Ndata = ceil(Tdata/dt); int idata = 0; double t = tbeg; int j = 0; while(t <= tend) { // Monitor Tdata if(idata >= Ndata) { idata = 0; } // Write output every Tdata periode if(idata == 0) { // Create sequence number for output file stringstream ss; ss << setfill(’0’); ss << setw(digit); ss << j; // Construct output filename string ofntp = ofnpre + "p-" + ss.str() + ".txt"; string ofnsp = odir + "/" + ofntp; string ofntc = ofnpre + "c-" + ss.str() + ".txt"; string ofnsc = odir + "/" + ofntc; cout << "Computing " << 100*j/ceil((tend-tbeg)/Tdata) << "% completed \r"; cout.flush();

10


j++; // Open output file const char *ofnp = ofnsp.c_str(); ofstream foutp; foutp.open(ofnp); const char *ofnc = ofnsc.c_str(); ofstream foutc; foutc.open(ofnc); // Write foutp << foutp << foutc << foutc <<

header in output file "# t = " << t << endl; "# x\ty" << endl; "# t = " << t << endl; "# x\ty" << endl;

// Write foutp << foutp << foutc << foutc <<

particles r[0].x << r[0].y << r[1].x << r[1].y <<

position to output file "\t"; endl; "\t"; endl;

// Close output file foutp.close(); foutc.close(); // Open output animation file std::ofstream afout(oafns.c_str(), std::ios::app); // Write plotting syntax in output animation file afout << "set title \"t = " << t << "\"" << endl; afout << "plot \"" << ofntc.c_str() << "\" w p pt 4 ps 8 notitle, \"" << ofntp.c_str() << "\" w circles notitle" << endl; // Close output animation file afout.close(); } // Calculate particles acceleration vect2 a[N], F[N]; vect2 rpc = r[0] - r[1]; double an = angle(rpc, vect2(0, 1)); vect2 npc = rpc>>1; vect2 ntg; int sign = 1; if(npc.x < 0) {sign = -1;}

10.7. LAMPIRAN

11

ntg.x = sign * npc.y; ntg.y = -sign * npc.x; F[1] = k * an * sign * vect2(1, 0); F[0] = (mp * g * sin(an) - length(F[1]) * cos(an) + length(v[0]) * eta) * ntg; a[0] = F[0] / mp; a[1] = F[1] / mc; // Increase idata and t idata++; t += dt; // Update particles velocity and position for(int i = 0; i < N; i++) { v[i] = v[i] + (a[i] * dt); r[i] = r[i] + (v[i] * dt); } rpc = r[0] - r[1]; npc = rpc>>1; r[0] = r[1] + (l * npc); } // Creating animation.gif file FILE *pipe = popen("gnuplot -persist", "w"); fprintf(pipe, "\n"); // Converting fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe, fprintf(pipe,

output data into gif animation in pipe "cd \’%s\’\n", odir.c_str()); "reset\n"); "set term gif animate\n"); "set output \"animation.gif\"\n"); "set xrange [-1:1]\n"); "set yrange [-0.5:1.5]\n"); "set xlabel \"x\"\n"); "set ylabel \"y\"\n"); "load \"%s\"\n", oafnt.c_str()); "set output\n"); "\n");

// Closing the pipe fclose(pipe); return 0; }

Kode C++ di atas dikompilasi dengan cara

12


g++ Inverted_Pendulum.cpp -o Inverted_Pendulum dan dipanggil dengan cara ./Inverted_Pendulum > params.txt

10.7.2

Script shell Bash

/* rwparams.h Library for reading parameters file Sparisoma Viridi | [email protected] 20141120 Create this class, define first just following two functions void readparam(const char*, const char*, int&) void readparam(const char*, const char*, double&) and test them 20150621 Add function writeparam for integer and double Add function writeparam for array of double Add function writeparam for comment Modify readparam for only two values line, that can differ comment line, single value line, two value line 20150920 Add namespace std, comment header file vect2.h, in order to remove error messages. 20150921 Add readparam for string */ #include #include <sstream> #include //#include "vect2.h" #ifndef rwparams_h #define rwparams_h using namespace std; void readparam(const char*, const char*, int&); void readparam(const char*, const char*, double&);

13

10.7. LAMPIRAN void void void void

writeparam(const writeparam(const writeparam(const writeparam(const

char*, char*, char*, char*,

const const const const

char*, int); char*, double); char*, int, double*); char*);

void readparam(const char*, const char*, string&);

void readparam(const char *pf, const char *pttrn, int &x) { ifstream fin; fin.open(pf); while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { stringstream ss; string col1; string col2; ss << line; ss >> col1; ss >> col2; if(col1 == pttrn) { stringstream ss; ss << col2; ss >> x; break; } } } fin.close(); } void readparam(const char *pf, const char *pttrn, double &x) { ifstream fin; fin.open(pf); while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { stringstream ss; string col1; string col2;

14


ss << line; ss >> col1; ss >> col2; if(col1 == pttrn) { stringstream ss; ss << col2; ss >> x; break; } } } fin.close(); } void readparam(const char *pf, const char *pttrn, string &x) { ifstream fin; fin.open(pf); while(fin.good()) { string line; getline(fin, line); bool COMMENT = line.find("#") != string::npos; bool BLANK = line.length() == 0; bool SINGLE_VALUE = line.find("\t") == string::npos; if(!COMMENT && !BLANK && !SINGLE_VALUE) { stringstream ss; string col1; string col2; ss << line; ss >> col1; ss >> col2; if(col1 == pttrn) { stringstream ss; ss << col2; ss >> x; break; } } } fin.close(); } void writeparam(const char *pf, const char *pttrn, int x) { fstream fout; fout.open(pf, fstream::out | fstream::app); fout << pttrn << "\t"; fout << x << endl; fout.close(); }

10.7. LAMPIRAN

void writeparam(const char *pf, const char *pttrn, double x) { fstream fout; fout.open(pf, fstream::out | fstream::app); fout << pttrn << "\t"; fout << x << endl; fout.close(); } void writeparam(const char *pf, const char *pttrn, int N, double *x) { fstream fout; fout.open(pf, fstream::out | fstream::app); fout << pttrn << "\t"; fout << N << endl; for(int i = 0; i < N; i++) { fout << x[i] << endl; } fout.close(); } void writeparam(const char* pf, const char* x) { fstream fout; fout.open(pf, fstream::out | fstream::app); fout << x << endl; fout.close(); } #endif dan /* vect2.h Library for 2-d vector Sparisoma Viridi | [email protected] 20141112 Create this class 20141119 Add destructor 20141120 Change title of this file from class to library */

#include #include <sstream> #include using namespace std;

15

16


#ifndef vect2_h #define vect2_h // Structure for 2-d vector struct vect2 { double x; double y; vect2(void); vect2(const vect2&); vect2(double, double); ~vect2(void); string strval(void); }; // Operator overloading vect2 operator+(vect2, vect2); vect2 operator-(vect2, vect2); double operator|(vect2, vect2); double length(vect2); vect2 operator*(double, vect2); vect2 operator*(vect2, double); vect2 operator/(vect2, double); vect2 operator>>(vect2, double); vect2 operator+=(vect2&, vect2); vect2 operator-=(vect2&, vect2); vect2 operator*(vect2, vect2); vect2 operator/(vect2, vect2); // Implementation vect2::vect2(void) { x = 0.0; y = 0.0; } vect2::vect2(const vect2& z) { x = z.x; y = z.y; } vect2::vect2(double zx, double zy) { x = zx; y = zy; } vect2::~vect2(void) { x = 0.0; y = 0.0; }

10.7. LAMPIRAN string vect2::strval(void) { stringstream ss; ss << "(" << x << ", " << y << ")"; return ss.str(); } vect2 operator+(vect2 a, vect2 b) { vect2 c; c.x = a.x + b.x; c.y = a.y + b.y; return c; } vect2 operator-(vect2 a, vect2 b) { vect2 c; c.x = a.x - b.x; c.y = a.y - b.y; return c; } double operator|(vect2 a, vect2 b) { double c = a.x * b.x + a.y * b.y; return c; } double length(vect2 r) { double c = sqrt(r|r); return c; } vect2 operator*(double a, vect2 b) { vect2 c; c.x = a * b.x; c.y = a * b.y; return c; } vect2 operator*(vect2 a, double b) { vect2 c; c.x = a.x * b; c.y = a.y * b; return c; } vect2 operator/(vect2 r, double l) { vect2 s; s.x = r.x / l; s.y = r.y / l; return s; }

17

18


vect2 operator>>(vect2 r, double s) { double l = length(r); vect2 u = (l > 0) ? r / l : vect2(0, 0); return u; } vect2 operator+=(vect2& a, vect2 b) { a = a + b; } vect2 operator-=(vect2& a, vect2 b) { a = a - b; } double double double double return }

angle(vect2 a, vect2 b) { la = length(a); lb = length(b); c = acos((a|b)/(la*lb)); c;

double double double double

angle(vect2 a, vect2 b, string option) { la = length(a); lb = length(b); c = acos((a|b)/(la*lb));

if (option == "DEG") {c = (180/m_pI)*c;} return c; } vect2 operator*(vect2 a, vect2 b) { vect2 c; c.x = a.x * b.x; c.y = a.y * b.y; return c; } vect2 operator/(vect2 a, vect2 b) { vect2 c; c.x = a.x / b.x; c.y = a.y / b.y; return c; } #endif Script shell Bash di atas dipanggil dengan cara ./rwparams.h

10.7. LAMPIRAN ./vect2.h

19

Komputasi Fisika 1 Simulasi dan Visualisasi Fisika dengan C, C++, Gnuplot, dan OpenGL

Recommend Documents