DOI: doi.org/10.21009/SPEKTRA.021.08

## TUMBUKAN TAK-ELASTIK PARTIKEL SEBAGAI MODEL TERBENTUKNYA ASTEROID
Sparisoma Viridi 1,a) , Budi Dermawan 2,b)
1 Prodi Fisika, FMIPA, Institut Teknologi Bandung, Jalan Ganesha 10, Bandung 40132
2 Prodi Astronomi, FMIPA, Institut Teknologi Bandung, Jalan Ganesha 10, Bandung 40132
Email: a) dudung@fi.itb.ac.id, b) budider@as.itb.ac.id

## Abstrak
Suatu asteroid diduga merupakan kumpulan partikel yang bertumbukan tak-elastik sehingga membentuk suatu agregat. Ikatan yang terjadi antar partikel diasumsikan hanyalah disebabkan gaya gravitasi. Gaya normal diperkenalkan untuk mencegah dua partikel berada pada posisi yang sama akibat gaya gravitasi. Parameter kontaktopi diperkenalkan dan teramati bahwa nilainya bertambah dengan semakin padatnya susunan  agregat  dan  pada  akhirnya  akan  mencapai  nilai  saturasi.  Dalam  sistem  yang  disimulasikan dengan jumlah partikel N = 22 diperoleh bahwa kontaktopi mencapai nilai saturasi pada C = 35.
Kata-kata kunci: partikel granular, dinamika molekuler, asteroid, kompaksi.

## Abstract
An  asteroid  is  assumed  to  be  constucted  from  particles  inealastic  collision,  which  later  forms  an agregate  of  particles.  Binding  between  particles  is  caused  only  by  gravitational  force.  Normal  force  is introduced to prevent two particles collapsing into single point as they are colliding and attracted to to each  other  due  to  gravitational  force.  Contactotpy  as  parameter  is  presented  and  it  is  observed  that  its value  increasing  as  the  agregate  becoming  denser,  and  at  the  end  it  achieves  saturation  value.  In  the simulated system with number of particles N = 22 it is obtained that the contactopy reaches saturation value C = 35.
Keywords: granular particles, molecular dynamics, asteroid, compaction.

## 1. PENDAHULUAN
Asteroid  merupakan  suatu  obyek  yang  menarik  apabila  dilihat  dari  bidang  material  granular, terlebih  setelah  tersedianya  foto-foto  resolusi  tinggi  mengenai  permukaan  asteroid  Itokawa  yang diperoleh  dalam  misi  Hayabusa  [1].  Terdapatnya  bongkahan  batuan  berukuran  lebih  besar  di  atas lapisan batuan yang berukuran lebih kecil dikaitkan dengan fenomena efek kacang Brasil (EKB) [2]. Simulasi mengenai material granular dalam sistem planet kemudian dilakukan dengan menggunakan pendekatan bola keras (hard sphere) [3], maupun bola lunak (soft sphere) [4] pada tumbukan antara dua buah partikel baik berbentuk bola ataupun lingkaran. Dilaporkan pula terdapatnya segregasi pada sistem dengan gravitasi sendiri (self-gravitating) suatu agregat material butiran [5], yang agak mirip dengan  EKB.  Eksperimen  dan  simulasi  fenomena  EKB  umumnya  menempatkan  partikel  granular dalam wadah dan semua partikel mengalami gaya gravitasi yang tegak lurus dengan permukaan alas
dalam  wadah.  Hal  ini  kurang  sesuai  dengan  keadaan  di  asteroid  yang  terlalu  kecil  dibandingkan dengan permukaan bumi yang dapat dianggap datar dibandingkan dengan wadah  material granular. Untuk itu suatu alas dalam wadah dengan gravitasi sendiri diusulkan dalam tulisan ini, yang diduga terbentuk dari batuan-batuan yang berinteraksi dengan gaya gravitasi dan gaya normal (untuk saat ini gaya gesek diabaikan).

## 2. SIMULASI
Terdapat  dua  interaksi  antar  partikel,  yaitu  atraktif  dan  repulsif.  Interaksi  atraktif  menjangkau rentang  yang  jauh  (long  range),  sementara  interaksi  repulsif  yang  sebatas  ukuran  partikel  (short range).  Interaksi  pertama  tak-lain  adalah  gaya  gravitasi,  di  mana  gaya  gravitasi  yang  bekerja  pada partikel i akibat adanya partikel j adalah
ij ij j i ij r r m m G G ˆ 2    (1)
dengan
j i ij r r r      , (2.a)
ij ij ij ij r r r r       , (2.b)
ij ij ij r r r   ˆ , (2.c)
dan G adalah  konstanta  universtal  gravitasi.  Lalu,  interaksi  kedua  adalah  gaya  normal  yang menggambil model pegas linier (linear spring-dashpot) [6]
  ij ij ij N ij N ij r dt d k
N ˆ sign              , (3)
dengan
  ij j i ij r R R    , 0 max  , (4.a)
ij ij v dt d    , (4.b)
            , 0 , 1 , 0 , 0 , 0 , 1 sign
x x x x (4.c)
dan R i , R j adalah jari-jari partikel i , j , berturut-turut.
Dengan  menggunakan  hukum  kedua  Newton  tentang  gerak,  percepatan  partikel i pada  saat t akibat berinteraksi dengan N –1 partikel lainnya dalam sistem dapat diperoleh melalui        N j i j ij ij i i N G m a 1 , 1    . (5)
Massa partikel m i diperoleh dari
i i i R m   3 3 4  . (6)
Selanjutnya,  metode  integrasi  numerik  Euler  akan  dapat  memberikan  kecepatan  dan  posisi  baru dari partikel pada waktu t + Δ t , yaitu     t a t v t t v i i i         (7.a)       t t v t r t t r i i i         . (7.b)
Dalam implementasinya, pengulangan Persamaan (5), (6.a), dan (6.b) yang akan mengubah posisi partikel-partikel  –  dikenal  sebagai  metode  dinamika  molekuler  (molecular  dynamics,  MD),  perlu diselipkan  suatu  langkah  untuk  membuat  kecepatan  partikel-partikel  menjadi  nol  sehingga  proses perhitungan menjadi lebih cepat. Suatu periode T digunakan sehingga akan terpenuhi bahwa
  .. , 2 , 1 , 0 , 0   n nT v i  . (8)
Bagaimana perngaruh dari nilai T ini akan diamati dalam simulasi.
Selanjutnya  susunan  partikel  perlu  dikarakterisasi  melalui  parameter  kontaktopi C yang  dapat diperoleh dari energi potensial kontak U C [7], yang dimodifikasi menjadi       N i N j ij C 1 1 sign  . (9)
Telah diamati bahwa pada fenomena EKB nilai C akan semakin besar dan mencapai maksimum saat fenomena ini terhenti [8].

## 3. HASIL DAN PEMBAHASAN
Nilai  Parameter-parameter  yang  digunakan  dalam  simulasi  MD  ini  diberikan  dalam  Tabel  1, apabila tidak diberikan keterangan lainnya. Akibat adanya gaya gravitasi ij G  maka susunan partikel akan menjadi semakin padat, yang dapat ditunjukkan dengan bertambahnya nilai kontaktopi C seperti dalam  Gambar  1.  Untuk  susunan  dengan N =  22 ini nilai saturasi C tercapai  dengan  nilai 35,  yang mungkin dapat sedikit agar berbeda apabila dilakukan variasi posisi awal partikel-partikel.
Gambar 1. Kontaktopi C sebagai fungsi waktu t dan beberapa ilustrasi susunan partikelnya untuk t = 0, 3, 9, dan 30.
Tabel 1. Nilai parameter-parameter standar dalam simulasi.
No. Parameter Value 1. k N 10 4 2. γ N 5×10 -1 3. T 10 2 4. G 10 -7 5. Δ t 10 -2 6. N 2 7. R i 1 – 2.5 8. ρ i 2000
Dengan menggunakan syarat awal yang sama dilakukan variasi nilai-nilai dari k N , γ N , dan T . Hasil akhir pada t = 30 untuk masing-masing variasi diberikan dalam Gambar 2, 3, dan 4.
0 5 10 15 20 25 30 35 0 5 10 15 20 25 30 C t
(a) (b) (c) (d) (e)
Gambar 2. Susunan akhir agregat untuk berbagai nilai k N : (a) 10 3 , (b) 2×10 3 , (c) 6×10 3 , (d) 10 4 , dan (e) 2×10 4 .
Variasi nilai k N tidak memberikan susunan akhir partikel yang terlalu berbeda sebagaimana pula variasi nilai γ N sebagaimana ditunjukkan dalam Gambar 2 dan 3. Hal ini dapat dipahami karena k N hanya  akan  mengubah  seberapa  besar  irisan  antar  dua  butiran  dan γ N menentukan  energi  yang hilang  dalam  setiap  tumbukan  antara  dua  partikel.  Akan  tetapi  terdapat  hal  yang  berbeda  dengan variasi  nilai T yang  memberikan  hasil  seperti  diberikan  dalam  Gambar  4,  di  mana  susunan  akhir agregat partikel ditentukan oleh nilai T ini.
(a) (b) (c) (d) (e) (a) (b) (c) (d) (e) (f)
Gambar 3. Susunan akhir agregat untuk berbagai nilai γ N :
(a) 0.1, (b) 0.2, (c) 0.3, (d) 0.4, dan (e) 0.6.
Gambar 4. Susunan akhir agregat untuk berbagai nilai T:
(a) 1, (b) 10, (c) 100, (d) 200, (e) 300, dan (f) 500.
Nilai T adalah suatu periode untuk  memaksa kecepatan semua partikel menjadi nol seperti telah diungkapkan  dalam  Persamaan  (8).  Langkah  ini  bertujuan  agar  simulasi  dapat  berjalan  lebih  cepat ketimbang menunggu suku yang mengandung γ N dalam Persamaan (3) berpengaruh. Bagaimana pengaruh dari T terhadap C susunan akhir partikel dapat dilihat dalam  Gambar 5, di mana terlihat bahwa dengan semakin besarnya nilai T akan semakin kecil nilai C . Evolusi kompaksi C sebagai fungsi waktu t dalam Gambar 1 dan bagaimana pengaruh dari parameter-parameter k N , γ N , dan T tidak  terlepas  dari  susunan  awal  partikel  dan  juga  distribusi  ukuran  partikel,  yang  diberikan dalam Gambar 6. Terlihat bahwa partikel berukuran rata-rata R = 1 memiliki jumlah paling banyak yang  kemudian  diikuti  dengan  berukuran  rata-rata  R  =  1.5,  sementara  ukuran  rata-rata  lainnya berjumlah kurang dari 30 %.
15 20 25 30 35 40 0 100 200 300 400 500 C T 0 2 4 6 8 10 0 1 2 3 4 5 6 N 2 R Gambar 5. Kontaktopi akhir C susuanan partikel sebagi fungsi T. Gambar 6. Distribusi ukuran partikel D = 2R yang digunakan dalam simulasi.
Untuk  satu  jenis  ukuran  partikel  nilai  kompaksi  dapat  diperkirakan  dalam  susunan  dua-dimensi berbentuk  HPC  (hexagonal  close  packed)  sebagaimana  diilustrasikan  dalam  Gambar  7  berikut. Diberikan susunan-susunan akhir partikel yang dianggap stabil untuk N = 2, 3, .. 7. Susunan-susunan lain  seperti  linier  dan  zig-zag  tidak  ditampilkan  karena  dianggap  masih  akan  berubah  menunju susunan-susunan seperti dalam gambar.
Gambar 7. Susunan HCP dalam dua-dimensi partikel berukuran sama dengan N = 1, 2, .., 7.
N = 3 N = 4 N = 5 N = 6 N = 2 N = 7
Kontaktopi untuk partikel berukuran sama dengan N = 1, 2, .., 7 diberikan dalam Gambar 8.
C avg = 0.6 N 2 - 3.78 N + 11.49 R² = 0.999 0 4 8 12 16 2 3 4 5 6 7 C avg N Gambar 8. Rata-rata kontaktopi C avg untuk berbagai jumlah partikel berukuran seragam N.
Suatu  perumusan  empirik  dapat  dibuat  untuk  mengaitkan  antara  nilai  rata-rata  kontaktopi C avg sistem  HCP  dua-dimensi  yang  stabil  untuk  partikel  berukuran  seragam  (monodispersed)  dengan jumlah partikel N dalam bentuk
49 . 11 78 . 3 6 . 0 2 mono avg,    N N C , (10)
yang  berlaku  untuk N ≥  4.  Persamaan  (10)  ini  dimulai  saat  susunan  akhir  partikel  pertama  kalinya memiliki  degenerasi,  yaitu  mulai N =  4.  Setelah  itu,  dengan  semakin  besar N semakin  banyak susunan akhir stabil dari partikel yang dapat dibentuk, atau degenerasi akan semakin besar.
Dengan menggunakan Persamaan (10) dalam sistem dengan ukuran partikel beragam yang telah disimulasikan  sebelumnya  diperoleh C avg =  217.73  atau C ≈  218  untuk  satu  susunan  akhir,  yang dalam  hal  ini  jauh  lebih  besar  dari  nilai  C  =  35  yang  diperoleh.  Untuk  sementara  dapat  dituliskan bahwa Persamaan (10) akan memberikan nilai maksimum kontaktopi yang dapat dicapai oleh sistem dengan ukuran partikel beragam (polydispersed)
mono avg, poly C C  , (11)
yang  tentunya  memerlukan  pembuktian  lebih  lanjut,  selain  dari  yang  ditunjukkan  beberapa kemungkin-annya dalam Gambar 9.
Susunan  dalam  Gambar  9  belum  menunjukkan  kemungkinan  yang  memenuhi  Persamaan  (11) karena titik N = 4 merupakan titik pertama dalam Persamaan (10). Diperlukan ilustrasi untuk susunan dengan  N  ≥  5  untuk  menjelaskannya.  Untuk  N  =  5  dapat  dibentuk  berbagai  susunan  akhir  partikel dan nilai kontaktopinya sebagaimana diberikan dalam Gambar 10. Susunan dengan nilai kontaktopi tertentu, misalnya C = 12 hanya dapat diperoleh bila rasio ukuran partikel mengikuti aturan tertentu. Sekilas dengan melihat Gambar 10 susunan akhir yang tidak memenuhi Persamaan (11) diduga tidak terlalu banyak sehingga persamaan tersebut untuk sementara tetap berlaku.
Gambar 9. Beberapa kemungkinan susunan akhir partikel dan nilai kompaktopinya untuk N = 1, 2, 3, 4.
Gambar 10. Beberapa kemungkinan susunan akhir partikel dan nilai kompaktopinya untuk N = 5.

## 4. SIMPULAN
Telah  dapat disimulasikan pembentukan  awal  asteroid  dalam  2-d.  Teramati  bahwa  kontaktopi C bertambah  dengan  semakin  kompaknya  susunan  partikel.  Parameter k N dan γ N tidak  terlalu berpengaruh  ke  susunan  akhir  partikel  yang  terbentuk.  Waktu  membuat v =  0  dengan T tertentu terlihat  berpengaruh,  di  mana  semakin  besar T akan  membuat  nilai C semakin  kecil.  Dugaan  batas atas  nilai C untuk  sistem  partikel  berukuran  beragam  telah  pula  diusulkan  dikaitkan  dengan C avg sistem partikel berukuran seragam.

## UCAPAN TERIMAKASIH
Terima  kasih  kepada  Kemenristekdikti  atas  dukungan  dalam  penelitian  ini  pada  tahun  2016 dengan skema Penelitian Unggulan Perguruan Tinggi (Nomor 583x/I1.C01/PL2016) dan Panitia SNF UNJ  2016  atas  undangan-nya  sebagai  pembicara  undangan  pada  sesi  paralel  dalam  melakukan sosialisasi hasil penelitian ini.

## DAFTAR ACUAN
[1] H.  Miyamoto,  H.  Yano,  D.J.  Scheeres,  S.  Abe,  O.  Barnouin-Jha,  A.F.  Cheng,  H.  Demura, R.W. Gaskell, N. Hirata, M. Ishiguro, T. Michikami, A.M. Nakamura, R. Nakamura, J. Saito, S. Sasaki, Regolith migration and sorting on asteroid Itokawa, Science 316 (2007), p. 1011- 1014.
C = 9 C = 7 C = 12 C = 7 N = 2 N = 3 N = 4 C = 3 C = 6 C = 9 C = 6
[2] G.V. Punyakoti, Brazil nuts help unlock the asteroid’s boulder mystery, Current Science 107 (2014), p. 945-946.
[3] D.C.  Richardson,  K.J.  Walsh,  N.  Murdoch,  P.  Michel,  Numerical  simulations  of  granular dynamics: I. Hard-sphere discrete element method and tests, Icarus 212 (2011), p. 427-437.
[4] S.  Matsumura,  D.C.  Richardson,  P.  Mitchel,  S.R.  Schwartz,  R.-L.  Ballouz,  The  Brazil  nut effect and its application to asteroids, MNRAS 443 (2014), p. 3368-3380.
[5] P. Sánchez, D.J. Scheeres, M.R. Swift, Impact driven size sorting in gelf-gravitating granular aggregates, 41st Lunar and Planetary Science Conference (2010), p. 2634.
[6] J. Schäfer, S. Dippel, D.E. Wolf, Force schemes in simulations of granular materials, J. Phys. I France 6 (1996), p. 5-20.
[7] S.  Viridi,  S.N.  Khotimah,  Widayani,  Novitrian,  Fluctuation  theorem  application  on  2-d granular materials configurations, Advanced Materials Research 1123 (2015), p. 12-15.
[8] S.  Viridi,  S.N.  Khotimah,  Widayani,  Novitrian,  Competition  between  potential  energy minimization  and  number  of  contact  points  in  two-dimension  granular  materials  exhibiting Brazil-nut effect, AIP Conf. Proc. 1677 (2015), p. 070001.