METODA ELEMEN HINGGA UNTUK MENYELESAIKAN PERSAMAAN DIFUSI NEUTRON SATU DIMENSI DUA GRUP

Daftar Isi Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(77-86)

METODA ELEMEN HINGGA UNTUK MENYELESAIKAN PERSAMAAN DIFUSI NEUTRON SATU DIMENSI DUA GRUP Utaja* , Topan Setiadipura** , Khairina Ns

ABSTRAK METODA ELEMEN HINGGA UNTUK MENYELESAIKAN PERSAMAAN DIFUSI NEUTRON SATU DIMENSI DUA GRUP. Distribusi neutron di dalam teras reaktor nuklir harus dapat ditentukan dengan teliti, karena menyangkut keselamatan reaktor maupun lingkungan. Penelitian ini menguraikan pemakaian metoda elemen hingga untuk menyelesaikan distribusi neutron berdasarkan teori difusi. Konsep dasar penyelesaian dengan metoda elemen hingga adalah membagi bentuk yang dianalisis menjadi sejumlah besar bentuk kecil yang dinamakan elemen. Untuk ini diambil elemen berbentuk batang. Di dalam elemen diberlakukan persamaan linier untuk menyatakan distribusi neutron. Karena persamaan ini bukan persamaan sebenarnya maka substitusi ke persamaan distribusi neutron akan menimbulkan sisa. Agar harga sisa sekecil-kecilnya, dipakai teori Galerkin. Hasil penyelesaian akhir berupa sejumlah persamaan linier yang disusun dalam persamaan matrik. Penyelesaian persamaan matrik dilakukan dengan komputer, yang program komputernya disusun sendiri. Pada dasarnya penyelesaian difusi neutron adalah penyelesaian persoalan eigen, maka hasil akhir adalah nilai eigen dan vektor eigen. Kata-kata kunci: elemen hingga, difusi neutron, Galerkin, nilai eigen, vektor eigen.

ABSTRACT. FINITE ELEMENT METHOD FOR ONE DIMENS ION TWO GROUPS NEUTRON DIFFUS ION SOLUTION. The neutron distribution in the reactor core should be determined accurately, because it will affect both the reactor and environment safety. This research describes the finite element method application for neutron distribution solution based on the diffusion theory. The fundamental concept solution with finite element method is geometrical discretisation to the finite shape that was named the element. In this research the bar element shape is choosen. Along the element, the linier function will applied for neutron distribution equation. Off the non exact equation, so that the substitution linier function to the neutron distribution equation will give any residuel. To minimized the residual , the Galerkin method should be used. The results are the amount of linier equations which are arranged on a matrix equation. The matrix solution is done by computer, where the software is custom made. The base of neutron diffusion solution is eigen problem, so the final results are eigen value and eigen vector. Keywords: finite element , neutron diffusion, Galerkin, eigen value, eigen vector

*

Pusat Rekayasa dan Perangkat Nuklir - BATAN Pusat Pengembangan Informatika Nuklir - BATAN

**

77

Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(77-86)

PENDAHULUAN Salah satu faktor yang penting di dalam disain reaktor nuklir adalah penentuan distribusi neutron. Distribusi neutron harus dapat ditentukan dengan teliti karena menyangkut pengoperasian reaktor nuklir. Salah satu teori yang dipakai untuk menentukan distribusi neutron di dalam teras reaktor nuklir adalah teori difusi. Model matematik difusi neutron berupa persamaan diferensial biasa orde dua homogen. Umumnya persamaan difusi neutron diselesaikan dengan metoda beda hingga (finite different). Di dalam makalah ini persamaan difusi diselesaikan dengan metoda elemen hingga, dengan harapan dapat diperoleh hasil yang lebih teliti. Hal ini disebabkan pada metoda elemen hingga ditribusi neutron di dalam elemen merupakan suatu fungsi, sedangkan pada metoda beda hingga harga distribusi neutron di dalam sebuah elemen merupakan konstanta. Pada penyelesaian dengan metoda elemen hingga, diperoleh sejumlah persamaan linier (sesuai dengan jumlah nodal yang dipakai) yang disusun dalam bentuk persamaan matrik. Syarat batas yang dikenakan adalah harga distribusi neutron di ujung adalah nol. Karena persamaan matrik yang didapat pada dasarnya merupakan persamaan homogen (di bagian kanan merupakan vektor berharga nol), maka penyelesaiannya merupakan penyelesaian untuk memperoleh nilai eigen dan vektor eigen. Ada dua penyelesaian pokok pada persoalan ini, yaitu penyelesaian invers matrik dan penyelesaian eigen. Invers matrik diselesaikan dengan metoda dekomposisi LU, sedangkan penyelesaian eigen dilakukan dengan metoda Invers iteration. Pemakaian dekomposisi LU sangat menguntungkan karena dapat mempercepat penyelesaian Invers iteration. Seluruh proses perhitungan dilakukan dengan komputer, yang perangkat lunaknya dikembangkan sendiri dengan memakai Visual Basic. Dengan metoda elemen hingga yang diselesaikan dengan komputer diharapkan penyelesaian difusi neutron dapat dilakukan dengan cepat dan hasil yang lebih teliti.

DASAR TEORI Uraian di dalam makalah ini dibatasi pada uraian distribusi neutron satu dimensi dengan dua grup energi. Persamaan distribusi neutron satu dimensi dua grup dinyatakan dengan [1]: G

G

-∇Dg∇Φg + ΣgΦg = χ/keff ∑νΣfg’Φg + ∑ Σsg’→gΦg g’=1

(1)

g’=1

dengan: Dg = konstanta difusi grup, Φg = fluks neutron grup, Σg = konstanta absorbsi neutron grup, ν = kecepatan neutron , keff = faktor perlipatan, Σsg’→g = konstanta moderasi grup g’ ke g, G = jumlah grup, χ=1 78


Persamaan (1) dapat ditulis dengan notasi diferensial lain: G

G

-d{Dg d(Φg)/dr}/dr + ΣgΦg = χ/keff ∑νΣfg’Φg + ∑ Σsg’→gΦg g’=1

(2)

g’=1

Di dalam teori elemen hingga, fluks neutron suatu grup Φg, dapat dinyatakan dengan :

i

x

(l-x)

i,j = adalah node

j

x = posisi antara i dan j

Gambar 1. Elemen garis Fluks neutron di x dinyatakan dengan: [ 2 ] Φg1 = Ni Φi(g1) + NjΦj(g1) Φg2 = Ni Φi(g2) + NjΦj(g2)

(3)

dengan: Ni, Nj = fungsi bentuk atau fungsi parameter (fungsi linier) Φi(g1) = fluks grup 1 di node i , Φj(g1) = fluks grup 1 di node j Φi(g2) = fluks grup 2 di node i , Φj(g2) = fluks grup 2 di node j Ni = (Xj – X)/(Xj-Xi), Nj = (X-Xi)/(Xj-Xk) Persamaan 3) dapat ditulis dengan bentuk : Φ = N ae dengan: N =

Ni 0 0 Ni

(4)

Nj 0 0 Nj

ae = [ Φi(g1) Φi(g2)

Φj(g1)

Φj(g2) ]T

Ruas kiri persamaan (2) merupakan persamaan yang tidak terkopel (Φg1 dan Φg2 tidak terhubung), sedangkan ruas kanan merupakan persamaan yang terkopel. Baik ruas kiri maupun ruas kanan mengandung suku Φg, sehingga ruas kanan dapat dipindahkan ke ruas kiri, sehingga persamaan (2) menjadi persamaan homogen (ruas kanan = (0) Substitusi persamaan (4) ke dalam persamaan (2) akan memberikan sisa R [2]: G

G

d{Dgd(Nae)/dr}/dr - Σg Nae - χ/keff ∑νΣfg’ Nae + ∑ Σsg’→g Nae ≠ 0 = R g’=1

(5)

g’=1

79


Untuk mendapatkan sisa R sekecil-kecilnya (penyelesaian mendekati harga benar), dipakai metoda Galerkin, dimana persamaan (5) dikalikan dengan Ne dan diintegralkan ke seluruh panjang elemen [ 2 ]: G

G

∫Le NT [d{Dgd(Nae)/dr}/dr - Σg Nae - χ/keff ∑νΣfg’ Nae - ∑ Σsg’→g Nae ] dr = 0 . .(6) g’=1

g’=1

Penyelesaian persamaan (6) akan memberikan: (Ke - Me ) Φe = 0 dengan: Ke = Ker + Keab Me = Mefg1 + Mefg2 + Metr Ker = D∫∫le (dNT/dr) (dN/dr) dr Keab = Σa ∫ NTN dr Mefg1 = 1/keff ∫ νΣfg1 NTN dr Mefg2 = 1/keff ∫ νΣfg2 NTN dr Metr = ∫ Σsg’_g NTN dr

(7) (7a)

Penyelesaian integral persamaan (7a) dengan bantuan sitem kordinat Serendipity [2] memberikan: 1 0 -1 0 0 1 0 -1 (7b) Ker = D/L -1 0 1 0 0 -1 0 1

Keab = ΣaL /6

2 0 1 0

0 2 0 1

1 0 2 0

0 1 0 2

+ Σs1->2L /6

2 0 1 0

0 0 0 0

1 0 2 0

0 0 0 0

Perbedaan antara Mefg1 dengan Mefg2 dan Metr terletak pada posisi koefisien matrik, dimana koefisien matrik dari Mefg1 terletak seperti hasil integral persamaan (7a). Koefisien matrik Mefg2 seperti pada matrik Mefg1 tetapi posisinya digeser satu baris ke atas, sedangkan untuk matrik Metr posisi koefisien matrik digeser satu baris ke bawah.

Mefg1 = LνΣfg1/keff

2 0 1 0

0 0 0 0

1 0 2 0

0 0 0 0 80


Mefg2

= LνΣfg2/keff

Metr = LΣsg’_g

0 0 0 0

2 0 1 0

0 0 0 0

1 0 2 0

0 2 0 1

0 0 0 0

0 1 0 2

0 0 0 0

Persamaan (7) adalah persamaan matrik (4x4) yang berlaku pada setiap elemen, dimana setiap elemen memiliki dua buah nodal (simpul). Untuk seluruh panjang yang ditinjau, perlu dilakukan assemblage (penggabungan) seluruh matrik dari setiap elemen. Setelah proses assemblage akan diperoleh persamaan matrik berikut: ( K - M)Φ= 0

(8)

Persamaan (8) akan memberikan harga Φ bukan nol, hanya bila determinan (K – M) nol. Penyelesaian persamaan 8) merupakan penyelesaian untuk mencari nilai eigen dan vektor eigen. Nilai eigen adalah kebalikan dari harga keff, atau dengan istilah lain harga keff adalah kebalikan harga eigen. Penentuan harga eigen Untuk menentukan harga eigen, matrik K dan M persamaan 8) dipisah menjadi persamaan berikut. (9) K Φ = MΦ Persamaan (9) diselesaikan dengan metoda iterasi dengan mengambil harga Φ awal sebagai harga coba-coba. Untuk itu dipakai metoda Invers iteration yang algoritma nya sebagai berikut.[ 3 ] 1. Mula-mula diambil harga awal vektor kolom Φ sebagai harga coba-coba (sebarang . harga). 2. Harga vektor kolom dikalikan dengan matrik M untuk memperoleh harga yk. yk = MΦ Φ 3. Lakukan evaluasi perhitungan berikut untuk harga k = 1, 2, 3 dan seterusnya. KΦ Φ(k+1) = yk Penyelesaian persamaan di atas dilakukan dengan metoda dekomposisi L-U. 4. Bentuk harga yk dengan: Φ(k+1) y(k+1) = MΦ 81


5. Bentuk harga eigen dengan: λ(xk+1) = (xTk+1 yk) / (xTk+1 yk+1) 6. Bentuk harga yk baru dengan yk+1 = yk+1 / (xTk+1 yk+1)1/2 Setelah didapat harga Yk+1, diulang lagi proses 2 sampai dengan 6. Untuk harga k ~ akan didapat : yk+1 ≅ vektor eigen λ(xk+1) ≅ harga eigen Seluruh proses dikerjakan dengan komputer yang programnya dikembangkan sendiri dengan Visual Basic 5.0 [4]

PERANGKAT LUNAK Perangkat lunak yang dikembangkan terdiri dari dua bagian, yaitu perangkat lunak untuk penyiapan data (PRE PROCESSOR), dan perangkat lunak untuk proses penyelesaian matrik dan analisis eigen (PROCESSOR) Pada program PRE PROCESSOR dilakukan diskretisasi geometri, meliputi jumlah nodal (simpul) dan jumlah elemen, kordinat nodal, data elemen, pemberian data material, dan pemberian syarat batas. Pada akhir proses, data yang sudah disusun disimpan dalam file data. Pemberian data material dapat dilakukan dengan mengambil data material yang sudah disiapkan, atau dengan membuat data baru. Program PROCESSOR akan membaca data yang dibuat oleh program PRE PROCESSOR. Mula-mula dibaca jumlah nodal dan jumlah elemen, dilanjutkan pembacaan data kordinat nodal. Selanjutnya dibaca data material, kemudian dilanjutkan dengan pembacaan data elemen. Pada pembacaan data elemen, dilakukan pembentukan matrik Ker, Keab, Mefg1, Mefg2, dan Metr. Pada saat pembentukan matrik, dilakukan optimasi, yaitu mengubah matrik bujur sangkar (n x n) menjadi matrik satu kolom dan sekali gus dilakukan proses assemblage. Metoda penyimpanan matrik Ker dan Keab berbeda dengan metoda penyimpanan matrik Mefg1, Mefg2 dan Metr, dimana matrik Ker dan Keab disimpan dengan metoda skyline, sedangkan matrik Mefg1, Mefg2, dan Metr disimpan dengan metoda banded. Proses selanjutnya adalah proses penentuan harga eigen dan vektor eigen. Pada proses ini invers matrik K dilakukan dengan metoda dekomposisi LU. Metoda ini menguntungkan karena pada proses iterasi untuk mendapatkan harga eigen, hanya bagian substitusi balik saja yang diperlukan. Proses penentuan eigen, dilakukan seperti algoritma di atas, dan dilakukan dengan memakai variabel dobel prosisi. Untuk membatasi proses iterasi, dilakukan pemberian toleransi berikut: TOL = (λk+1 - λk) / λk+1 82


Besarnya toleransi dapat dipilih, dianjurkan kurang dari 0.001. Harga keff merupakan kebalikan dari λ. Seluruh program komputer didulis dengan Visual Basic 5.0 [ 4 ]

HASIL DAN BAHASAN. Analisis dilakukan pada difusi neutron satu dimensi dua grup dengan data material berikut: Tabel 1. Data material No 1

D1 1,4360

Σ1 0,0095042

D2 0,36353

Σ2 0,075058

Σ12 0,017754

νΣf1 0,0058708

νΣf2 0,096067

Keterangan: D1 = koefisien difusi neutron cepat, D2 = koefisien difusi neutron termal Σ1 = tampang lintang absorbsi neutron cepat Σ2 = tampang lintang absorbsi neutron termal Σ12 = koefisien transfer dari neutron cepat ke neutron termal νΣf1 = koefisien fisi neutron cepat νΣf2 = koefisien fisi neutron termal Panjang yang dianalisis 50 cm, dibagi menjadi 20 elemen dan mempunyai 21 nodal. Seluruh panjang dianggap terdiri dari satu macam material seperti Tabel 1. Elemen memiliki panjang yang sama, setiap elemen memiliki dua nodal di ujungnya. Seluruh data disiapkan dengan program penyiapan data (PRE PROCESSOR). Hasil analisis dinyatakan denan Tabel 2. Tabel 2. Hasil distribusi fluks relatif No 1 2 3 4 5 6 7 8 9 10 11

Neutron cepat 0 0,1113 0,2197 0,3226 0,4177 0,5024 0,5748 0,6330 0,6566 0,7017 0,7104

Neutron termal 0 1,4604 E-3 2,8498 E-3 4,1784 E-3 5,4070 E-3 6,5035 E-3 7,4401 E-3 8,1938 E-3 8,7457 E-3 9,0825 E-3 9,1956 E-3

No 12 13 14 15 16 17 18 19 20 21

Neutron cepat 0,7017 0,6756 0,6330 0,5748 0,5024 0,4177 0,3226 0,2197 0,1113 0

Neutron ltermal 9,0825 E-3 8,7457 E-3 8,1937 E-3 7,4400 E-3 6,5032 E-3 5,4063 E-3 4,1762 E-3 2,8432 E-3 1,4396 E-3 0

83


Hasil yang ditampilkan baik neutron cepat maupun kecenderungan sama, yaitu maksimum di bagian tengah. berikut.

pada Tabel 2 adalah harga relatip distribusi neutron, neutron termal. Kedua hasil distribusi memiliki berharga nol pada kedua ujungnya dan berharga Hal ini dapat dilihat pada Gambar 1 dan Gambar 2

Gambar 1. Grafik Fluks cepat vs Posisi

Gambar 2. Gafik Fluks termal vs Posisi Harga koefisien efektip keff sebesar 0.42809, yang dicapai pada iterasi ke 10. Untuk sementara hasil pada Tabel 2 belum divalidasi. Seberapa jauh perbedaan hasil antara metoda elemen hingga dengan metoda beda hingga dapat dilihat pada Tabel 3 berikut untuk perhitungan difusi satu dimensi. Hasil pada Tabel 3 didapat dari perhitungan untuk panjang 20 cm, material seperti Tabel 1, dengan diskretisasi 20. Tabel 3. Perbandingan hasil distribusi fluks relatip antara MEH dan MBH No 1 2 3 4 5 6 7 8 9 10 11

MEH 0 4,9381 E-4 9,7507 E-4 1,4323 E-3 8,8542 E-3 2,2306 E-3 2,5519 E-3 2,8104 E-3 2,9999 E-3 3,1154 E-3 3,1541 E-3

MBH 0 4,9353 E-4 9,7360 E-4 1,4271 E-3 8,8471 E-3 2,2060 E-3 2,5101 E-3 2,7457 E-3 2,9881 E-3 2,9881 E-3 2,9066 E-3

MEH = metoda elemen hingga,

No 12 13 14 15 16 17 18 19 20 21

MEH 3,1154 E-3 2,9999 E-3 2,8105 E-3 2,5519 E-3 2,2306 E-3 1,8543 E-3 1,4323 E-3 9,7507 E-4 4,9383 E-4 0

MBH 2,9066 E-3 2,7458 E-3 2,5101 E-3 2,2060 E-3 1,8471 E-3 1,4271 E-3 9,7360 E-4 4,9353 E-4 0

MBH = metoda beda hingga 84


Dari Tabel 3 dapat dilihat adanya perbedaan hasil dari MEH dengan MBH, di daerah maksimum (nomer 11) sebesar 8%. Hal ini disebabkan adanya perbedaan pendekatan penyelesaian seperti yang sudah diuraikan di atas. Perbandingan lain antara hasil dengan MEH dan MBH dapat dilihat pada harga eigen untuk beberapa diskretisasi, seperti ditampilkan pada Tabel 4. Tabel 4. Harga eigen dari tiga macam diskretisasi Jumlah diskretisasi 10 20 50

MEH 0,12980 0,13043 0,13061

MBH 0,1111 0,12061 0,12650

Dari Tabel 4 dapat dilihat konvergensi (∆ eigen/∆diskretisasi) pada MEH ratarata 0,1 dari MBH, yang berarti MEH lebih cepat konvergen dibanding MBH.

KESIMPULAN Metoda elemen hingga dapat dipakai untuk menghitung distribusi neutronik satu dimensi dua grup, berbasis teori difusi. Program yang dikembangkan masih perlu diverifikasi dan divalidasi dengan program sejenis yang sudah teruji, misal dengan bench marking. Untuk pemakaian dengan jumlah elemen yang tinggi, dan bermacam material yaitu lebih dari 2000 elemen dengan berbagai material masih perlu penelitian lebih lanjut.

UCAPAN TERIMA KASIH Kami sampaikan terima kasih pada KPTF PRPN-BATAN yang telah membantu perbaikan dan editing pada makalah kami.

DAFTAR PUSTAKA 1. WESTON M.STACEY, “Nuclear Reactor Physics”, John Wiley & Sons, INC, New York, USA, 2001. 2. FRANK L.STASA, ”Applied Finite Element Analysis for Engineers”, CBS College Publishing, New York, USA, 1985. 85


3. KLAUS-JURGEN BATHE, ”Finite International, New York, USA, 1996.

Elemen

Procedure”,

Prentice

Hall.

4. EVANGELOS PETROUTSOS, ”Visual Basic 5.0”, Sybex, San Fransisco, USA, 1997.

DAFTAR RIWAYAT HIDUP

1. 2. 3. 4. 5. 6.

7.

Nama : Ir Utaja Tempat/Tanggal Lahir : Wates, 28 Nopember 1946 Instansi : PRPN-BATAN Pekerjaan / Jabatan : Peneliti Utama Riwayat Pendidikan : • S1 Teknik Mesin UGM (1983) Pengalaman Kerja : • Ka. Sub bidang disain, 1990 – 2000 • Ka. Sub bidang komponen Nuklir, 2000 – 2002 • Peneliti Madya – 2007 • Peneliti Utama (2007) - sekarang Publikasi Ilmiah : • Pemakaian koordinat luasan pada FEM • Analisis lintasan elektron dalam medan elektrostatik • MEH berbasis elemen beam untuk analisis defleksi poros turbin

Daftar Isi 86

METODA ELEMEN HINGGA UNTUK MENYELESAIKAN PERSAMAAN DIFUSI NEUTRON SATU DIMENSI DUA GRUP

Recommend Documents