APLIKASI PAKET PROGRAM MOLDY UNTUK KARAKTERISASI SIFAT BAHAN Fe, Pb, Bi DAN PENDINGIN REAKTOR Pb-Bi

Kembali Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir XVII, Agustus 2006 (119-130)

APLIKASI PAKET PROGRAM MOLDY UNTUK KARAKTERISASI SIFAT BAHAN Fe, Pb, Bi DAN PENDINGIN REAKTOR Pb-Bi Alan Maulana*, Zaki Suud*, Hermawan K.D**, Khairurijal*

ABSTRAK APLIKASI PAKET PROGRAM MOLDY UNTUK KARAKTERISASI SIFAT BAHAN Fe, Pb, Bi DAN PENDINGIN REAKTOR Pb-Bi. Paket program Moldy adalah suatu program open source yang dipakai untuk simulasi dinamika molekul bahan cair atau padat. Program ini dapat digunakan untuk simulasi sistem bahan monoatomik atau poliatomik. Dalam simulasi ini konfigurasi atom-atom bahan Fe, Pb, Bi diperoleh dari data crystal structure yang berupa parameter kisi dengan jenis strukturnya. Interaksi antar atom-atom Fe-Fe, Pb-Pb dan Bi-Bi diasumsikan memenuhi potensial Lennard-Jones. Parameter Lennard-Jones masing-masing atom diambil dengan melalui fitting data yang terdapat dalam literature yang menghitung interaksi antar-atom dengan metoda ab-initio. Hasil yang diperoleh dalam simulasi ini adalah kurva fungsi distribusi radial dan Mean Square Displacement sebagai fungsi dari waktu yang dapat dipakai untuk melihat keadaan fase system sebagai fungsi dari temperatur. Untuk kristal Fe pada temperatur 300 K ,773 K dan 1000 K Fungsi distribusi radial menunjukkan system dalam keadaan padat sedangkan pada temperatur 1809 K system dalam keadaan cair. Kata-kata kunci:

dinamika molekul, potensial Lennard-Jones, ab-initio, Liquid Lead-Bismuth, fungsi distribusi radial, mean square displacement, Fast Breeder Reactor.

ABSTRACT MOLDY APPLICATION FOR Fe, Pb, Bi AND Pb-Bi REACTOR COOLANT CHARACTERIZATION. MOLDY is an open source program which can be used to simulated mono atomic or polyatomic condensed matter. The initial configuration of crystal Fe, Pb and Bi are taken from crystal structure data by means lattice parameters and structure types. Interaction potentials between FeFe, Pb-Pb and Bi-Bi are assumed follow the Lennard-Jones potential. The Lennard-Jones potential parameters had been derived by fitting the data available in the literature which was calculated by abinitio methods with the Lennard-Jones formula. The results of this simulation are radial distribution functions and mean square displacements of the system which can be used to identify whether the system in the solid state or liquid state as a function of temperature. Radial distribution function and mean square displacement of iron crystal at 300 K, 700 K and 1000 K showed the system was in the solid state, while at 1809 K was in the liquid state. Keywords: molecular dynamic, Lennard-Jones potential, ab-initio, Liquid Lead-Bismuth, radial distribution function, mean square displacement, fast breeder reactor

*

Departemen Fisika ITB Departemen Fisika Teknik ITB Jl.Ganesha 10, Bandung 40132, Telp 62-22-250-0834, FAX 62-22-250-6452

**

119

Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir XVII, Agustus 2006 (119-130)

PENDAHULUAN Paket program Moldy adalah suatu program open source yang dipakai untuk simulasi dinamika molekul bahan cair atau padat yang berbasis sistem LINUX maupun Windows. Program ini memakai bahasa C dan dapat digunakan untuk simulasi sistem bahan monoatomik atau poliatomik. Fungsi potensial yang sudah ada dalam program ini adalah Lennard-Jones, Buckingham termasuk Born-Mayer, MCY dll, serta dapat ditambah dengan potensial lain melalui modifikasi program. Dalam tulisan ini akan disajikan hasil simulasi yang dilakukan penulis terhadap bahan-bahan Fe, Pb, Bi dan Pb-Bi dengan paket MOLDY yang merupakan penelitian awal tentang korosi baja dalam Pb-Bi cair. Sebelumnya penulis telah membuat modifikasi program molekular dinamik Pb-Bi dan telah disampaikan pada LKSTN XVI yang tampaknya masih banyak kekurangan terutama untuk sistem lebih dari 2 jenis atom. Sehingga dengan diperolehnya paket program MOLDY diharapkan dapat menutupi kekurangan hasil modifikasi program yang dikembangkan. Timbal-Bismuth cair dijadikan topik penelitian karena menjadi kandidat yang sangat kuat untuk menjadi pendingin Fast Breeder Reactor (Reaktor Pembiak Cepat) atau Fast Reactor (Reaktor cepat). Sifat-sifat yang menguntungkan Pb-Bi cair dalam aplikasi pendingin reaktor nuklir adalah titik lelehnya 123,5 °C dan titik didihnya 1670 °C dengan perubahan volume ketika berubah menjadi padat sebesar 1,5 %. Keuntungan lain penggunaan Pb-Bi cair dalam reaktor nuklir bukan saja inert terhadap air dan udara tetapi juga dapat menghasilkan 30 neutron untuk setiap 1 GeV proton datang. Kelemahan cairan logam ini adalah sangat agresif terhadap besi dan Stainless Steel terutama pada temperature tinggi. Oleh karena itu dalam reaktor nuklir berpendingin Pb-Bi, Stainless Steel akan mengalami korosi. Penelitian fenomena korosi Stainless Steel dalam Pb-Bi cair pada umumnya dilakukan secara eksperimen. Kegiatan ini tampaknya belum memungkinkan di Indonesia karena fasilitas yang belum memadai. Oleh karena itu penulis melakukan penelitian ini dengan simulasi komputer. Penelitian ini ditujukan untuk memahami mekanisme atomik fenomena korosi Stainless Steel dalam Pb-Bi cair sebagai fungsi dari temperatur dan waktu. Sebagai langkah awal akan dilihat sifat pelelehan masing-masing bahan diatas dengan menghitung Radial Distribution Function (RDF) dan Mean Square Displacement (MSD) yang dapat dipakai untuk menghitung konstanta difusi.

DASAR TEORI Simulasi menggunakan metoda dinamika molekul pada dasarnya dimulai dengan menentukan konfigurasi atom-atom bahan yang ditinjau. Masing-masing atom tersebut saling berinteraksi satu sama lain yang dalam hal ini diasumsikan memenuhi 120


potensial Lennard-Jones, kemudian atom-atom tersebut diberikan kecepatan awal secara random. Persamaan gerak sistem diatas dipecahkan secara numerik dengan komputer. Konfigurasi awal bahan Fe, Pb dan Bi ditentukan dari tabel kristalografi sbb: Tabel 1. Parameter Kisi Fe, Pb dan Bi No

Nama

Struktur

Parameter Kisi (Angstrom)

BCC FCC Monoclinic

2,8665 4,9508 6,674

Sudut

Unsur 1 2 3

Fe Pb Bi

2,8665 4,9508 6,117

2,8665 4,9508 3,304

α=β=γ=90° α=β=γ=90° α=90°; β=110,330° γ=90°

Dari data diatas ditentukan koordinat posisi atom untuk masing-masing unsur dalam satu sel satuan. Setelah satu sel satuan diperoleh maka sel satuan tersebut diperbanyak dalam arah XYZ sehingga jumlah sel satuan tersebut menjadi bertambah sesuai dengan yang kita inginkan memakai persamaan:

r = ruc + n x a + n y b + n z c

(1)

ruc adalah posisi sel satuan ; a , b dan c adalah parameter kisi; nx,ny dan nz adalah jumlah sel satuan dalam arah x,y dan z. Kecepatan awal masing-masing atom diberikan secara random dengan nilai antara 0-1 melalui random generator. Potensial interaksi antar Fe-Fe, Pb-Pb, Bi-Bi dan Pb-Bi diambil dari literatur yang telah menghitung secara ab-initio. Data ini difitting menggunakan persamaan potensial Lennard-Jones berikut:

 σ 12  σ  6  u (r ) = 4ε   −     r    r 

(2)

u(r) adalah energi potensial, ε parameter energi dan σ adalah parameter jarak terdekat antar atom. Hasil fitting tampak dalam table 2.

121


Tabel 2. Parameter Lennard-Jones Fe,Pb,Bi dan Pb-Bi No 1 2 3 4

Unsur Fe-Fe Pb-Pb Bi-Bi Pb-Bi

ε (eV) -0,62 -0,125 -0,075 -0,12

σ(Å) 2,26 2,75 2,8 2,75

Secara grafis parameter Lennard-Jones tampak pada gambar 1:

Gambar 1. Potensial Lennard-Jones Pb, Bi dan Pb-Bi Dari tabel 2 tampak kedalaman sumur potensial Fe paling tinggi diantara yang lain yang secara kualitatif berarti ikatan antar Fe-Fe lebih kuat dari yang lain. Dari energi potensial diatas dapat diperoleh gaya antar atom sbb: 13 7 ε  σ  σ   du(r ) = 24 2  −    F =− σ   r   r   dr

(3)

Gaya yang dialami oleh atom i akibat berinteraksi dengan N atom adalah : N Fi = ∑ Fij , (i,j =1,…,N)

(4)

j ≠i

Fij adalah gaya yang dialami oleh atom i akibat atom j yang mana atom i≠j. Untuk mengetahui posisi dan kecepatan atom i setelah mengalami gaya dari N atom lain maka diperlukan persamaan gerak Newton : 122


∂vi ∂ri = mi Fi = mi ∂t ∂t

;

(i = 1,…,N)

(5)

mi massa atom partikel i , vi dan ri adalah kecepatan dan posisi atom i. Persamaan diferensial ini dalam paket program MOLDY dipecahkan secara numerik melalui algoritma Beeman. Dari simulasi diatas dapat diperoleh posisi dan kecepatan dari masing-masing atom setiap saat. Dari mekanika statistik dapat dihubungkan kuantitas mikroskopik ini dengan kuantitas makroskopik seperti temperatur, tekanan, RDF, MSD dll. Hubungan antara temperatur dengan kuantitas mikroskopik diatas adalah sbb:

T=

1 2 mi v i ∑ k B (3 N − N C ) i

(6)

N adalah jumlah atom dalam simulasi, Nc jumlah constraint , kB adalah konstanta Boltzmann dan vi adalah kecepatan atom i. Temperatur sistem dapat diset sesuai sistem yang ditinjau, dalam MOLDY terdapat beberapa metoda untuk menentukan temperatur tersebut diantaranya Scaling dan Noose-Hoover thermostat. Persamaan Fungsi distribusi radial (RDF) adalah sebagai berikut: g (r ) =

N ( r , ∆r )

(7)

1 N ρV ( r , ∆ r ) 2

Di mana ρ adalah rapat atom , V(r,∆r) volume kulit bola pada jarak r±∆r . Fungsi distribusi radial g(r ) merupakan ukuran untuk melihat sejauh mana atom-atom mengatur posisinya pada temperatur dan waktu tertentu. Sehingga dapat dibedakan secara kualitatif apakah suatu sistem dalam keadaan padat atau cair. Kuantitas makroskopik lain yang penting dalam dinamika molekul adalah Mean Square Displacement (MSD). Pada temperatur tinggi atom-atom dalam sistem bergerak setiap saat. Dalam molekular dinamik ini berarti merupakan iterasi dari pemecahan persamaan diferensial dengan jumlah timestep tertentu. Dengan demikian perpindahan kuadrat dari atom-atom setiap saat dapat dirata-ratakan. Kuantitas ini dapat dikaitkan dengan perhitungan konstanta difusi. Persamaan untuk Mean square displacement (MSD) adalah sbb:

MSD= ∆r 2 (t ) = r (t ) − r (0)

2

(8)

123


Konstanta difusi diri dapat diperoleh dari hubungan:

1 6 Nt

D=

N

i

i

∑ r (t ) −r (0) i =1

(9)

Dengan demikian koefisien difusi diri berbanding lurus dengan kemiringan kurva MSD vs timestep.

HASIL DAN ANALISIS Hasil yang diperoleh tampak pada gambar 2 dibawah yang berupa fungsi distribusi radial kristal Fe. Pada suhu 300 K puncak fungsi distribusi radial lebih banyak yang muncul dengan intensitas yang tinggi menandakan sistem dalam keadaan padat. Demikian juga untuk temperatur 773 K dan 1000 K puncaknya banyak tetapi intensitasnya berkurang dengan naiknya temperatur. Sedangkan pada temperatur 1809 K puncak kurva fungsi distribusi radial menjadi berkurang dan intensitasnya semakin menurun. Hal ini terjadi karena posisi atom-atom Fe pada temperatur ini semakin acak yang berarti sistem dalam keadaan cair. Hal ini sesuai dengan hasil eksperimen yang mana titik leleh Fe adalah 1809 K. Hasil simulasi yang berupa mean square displacement juga memperlihatkan bahwa pada suhu 1809 K kemiringannya tinggi yang berarti konstanta difusi pada temperatur tersebut lebih tinggi daripada temperatur yang lain. Semakin tinggi kemiringan kurva mean square displacement semakin acak susunan atomnya. Radial Distribution Function of Fe 6

5

4

3

rdf

Fe T=300 T=773 K 2 T=1000K T=1809 K

1

0 0

2

4

6

8

10

12

-1

r (angstrom)

Gambar 2. Kurva fungsi distribusi radial pada beberapa temperatur

124


Pada gambar 4 tampak kurva fungsi distribusi radial kristal Pb (Timbal) pada temperatur 200 K , 300K dan 600 K. Kurva untuk temperatur 600 K menunjukkan sistem dalam keadaan cair yang faktanya pada temperatur tersebut merupakan titik leleh Pb. Demikian juga temperatur 300 K sistem dalam keadaan cair tetapi intensitasnya lebih tinggi. Sedangkan pada temperatur 200 K pada puncak ke 2 sudah akan muncul puncak baru yang menunjukkan bila sistem akan menuju keadaan padat bila temperatur diturunkan. MSD kristal Fe

Mean Square displacement

12

10

8

T= 300 K T=773 K T= 1000 K T=1809

6

4

2

0 0

2

4

6

8

10

12

14

16

t( x 0,5 pSec)

Gambar 3. Mean Square Displacement pada beberapa temperatur Fungsi Distribusi Radial Pb 5

RDF Pb T=600 K RDF Pb T=300 K RDF Pb T=200 K

4

g(r)

3

2

1

0 0

2

4

6

8

10

12

-1

r (Angstrom)

Gambar 4. Fungsi distribusi radial Kristal Pb pada beberapa temperatur

125


Kurva fungsi distribusi radial Bismuth pada temperatur 100 K, 544 K dan 773 K tampak pada Gambar 5. Temperatur leleh Bismuth adalah 544 K sehingga kurva fungsi distribusi radial intensitasnya rendah. Pada temperatur rendah (100 K) tidak muncul puncak baru tetapi intensitasnya jauh lebih tinggi. Gambar 6 adalah kurva untuk sistem Pb-Bi yang hasilnya tidak jauh berbeda dengan hasil yang telah penulis lakukan dengan program yang dikembangkan sendiri[6]. Fungsi Distribusi Radial Bismuth 4.5

4

3.5

T=100 K

3

T=544 K

g(r)

2.5

T=773 K

2

1.5

1

0.5

0 0

2

4

6

8

10

12

-0.5

r (Angstrom)

Gambar 5. Fungsi distribusi radial kristal Bismuth

Fungsi Distribusi Radial Pb-Bi 2.5 rdf pb-bi T=396,5 K rdf pb-bi T= 773 K rdf pb-bi T=273 K

2

rdf

1.5 1 0.5 0 0

2

4

6

8

10

12

-0.5

r ( angstrom )

Gambar 6. Fungsi distribusi radial Pb-Bi

126


KESIMPULAN • Paket program MOLDY dapat dipakai untuk simulasi bahan yang diantaranya dapat menentukan struktur bahan pada temperatur yang diinginkan melalui fungsi distribusi radial dan mean square displacement. • Kemiringan kurva mean square displacement berbanding lurus dengan konstanta difusi diri dari bahan. • Hasil simulasi dinamika molekul kristal Fe menunjukkan bahwa pada 1809 K Kristal Fe dalam fasa cair sedangkan pada temperatur dibawahnya dalam fasa padat.

UCAPAN TERIMAKASIH Penelitian ini dibiayai oleh Kementrian Negara Riset dan Teknologi melalui program Doktor PUSPIPTEK.

DAFTAR PUSTAKA 1. J.M.HAILE,Molecular Dynamics Simulations, Elementary Methods, John Wiley & sons, INC, 1997. 2. Y.QI and M.TAKAHASHI, Study on Corrosion phenomena of steels in Pb-Bi flows, Proceeding ICONE 11, Tokyo, Japan 2003. 3. D.FRENKEL, B.SMITH,Understandings Molecular Simulation, Acadmic Press, 2002. 4. Y.QI and M.TAKAHASHI, Computer Simulation of Diffusion of Pb-Bi Eutectic in Liquid Sodium By Molecular Dynamic Method, Proceeding ICONE 10, Arlington, USA, 2002. 5. A.MAULANA., Z.SUUD, HERMAWAN K.D., dan Khairurijal, Corrosion Study of Steels In Pb-Bi Using Molecular Dynamics Method, COE-INES-Indonesia Internasional Symposiun 2005.

127


6. A.MAULANA., Z.SUUD, HERMAWAN K.D., dan Khairurijal, Studi Sifat Pb-Bi Dengan Metode Molekular Dinamik, LKSTN XVI, P2TIK-BATAN, Jakarta, 2005. 7. A.MAULANA., Z.SUUD, HERMAWAN K.D., dan Khairurijal, Aplication Of Molecular Dynamic Methods To Study Steels Corrosion Phenomena in Liquid Lead-Bismuth Cooled Reactors. 8. K. REFSON, ‘Moldy: a portable molecular dynamics simulation program for serial and parallel computers, Computer Physics Comunications, 126(3), (2000) 309-328

DISKUSI

FAIZAL RIZA 1. 2. 3.

Fase sistem sebagai fungsi dari temperatur, apa yang dimaksud dengan fase sistem ? Sejauh mana konduktivitas SS dalam Pb Bi cair ? Berapa besar dan lama hingga terjadi korosi ? Untuk kristal Fe pada temperature 300 K, 773 K dan 1000 K. Apakah hanya melihat kondisi kristalnya saja ? (padat, cair dan gas)

ALAN MAULANA 1. 2. 3.

Fase sistem disini adalah keadaan sistem apakah dalam keadaan cair atau padat. Dalam presentasi ini baru dilakukan simulasi untuk atom-atom Fe, Pb, Bi dan Pb-Bi secara individual. Sehingga data-data yang ada belum bisa mengatakan apa-apa tentang korosi. Ya, karena saya sedang mencoba melihat kemampuan MOLDY. Selanjutnya Fe tersebut akan disimulasikan dalam interaksinya dengan Pb atau Pb-Bi cair pada beberapa temperatur.

128


DARMAWAN Mohon dijelaskan lebih lanjut fungsi Mean Square Displacement karena dari hasil yang ditulis di ABSTRAK ini masih belum lengkap ? ALAN MAULANA Mean Square Displacement secara matematis dituliskan sebagai beriku ;

MSD= ∆r 2 (t ) = r (t ) − r (0)

2

Yang mana kemiringan kurvanya bisa dipakai untuk menghitung konstanta difusi serta bisa dipakai untuk melihat apakah sistem dalam keadaan padat atau cair. Hasil simulasi yang saya lakukan terhadap Fe pada 1809 K kemiringannya sangat tinggi. Dibandingkan 300 K, 773 K dan 1000 K yang menunjukkan pada temperature 1809 sistem dalam keadaan cair TOPAN Parameter penting apa yang perlu diamati dalam fenomena agresifnya Pb-Bi terhadap korosif ? ALAN MAULANA Parameter yang penting diamati alam fenomerna korosi ini adalah posisi tom-atom Fe, Pb, Bi sebagai fungsi dari waktu pada temperature tertentu. Sehingga bisa dihitung kecepatan korosinya.

129


DAFTAR RIWAYAT HIDUP 1. Nama

: Alan Maulana

2. Tempat/Tanggal Lahir

: Bandung, 16 Juli 1964

3. Instansi

: BATAN

4. Pekerjaan / Jabatan

:

5. Riwayat Pendidikan

:

• S1 Fisika, ITB • S2 Teknik Material, ITB • Saat ini sedang S3 di Fisika ITB 6. Pengalaman Kerja

:

• BATAN 7. Makalah yang dipublikasikan •

Corrosion Study of SS in Liquid Pb-Bi using molecular dynamic methods, COE-INES, Bandung (2005)

130

APLIKASI PAKET PROGRAM MOLDY UNTUK KARAKTERISASI SIFAT BAHAN Fe, Pb, Bi DAN PENDINGIN REAKTOR Pb-Bi

Recommend Documents