Nemxyz-part1

Tim Neutronik 2009 Bidang Komputasi PPIN-BATAN Metoda Ekspansi Nodal untuk Solusi Persamaan Difusi Neutron 3-D X-Y-Z Usaha komunitas komputasi sains dan teknologi nuklir untuk melakukan perhitungan distribusi daya dan parameter reaktor lain telah begitu intensif dilakukan, banyak metodametoda komputasi unggul lahir dalam rangka perhitungan ini. Hal ini karena kelakuan dari reaktor nuklir harus diprediksi sebelum betul-betul dioperasikan, terutama pada tahap desain perhitungan ini sangat penting diantaranya untuk memahami kelakukan subkritis atau superkritis dari reaktor.1 Software computer untuk simulasi teras reaktor nuklir terdiri dari dua bagian besar yang saling terkait satu sama lain, namun dalam tahap pengembangan program keduanya dapat dipisahkan. Bagian pertama adalah bagian neutronik untuk menghitung distribusi fluks neutron, daya yang dihasilkan, factor kritikalitasm dan lainnya. Bagian kedua adalah termohidraulik terkait dengan transfer atau distribusi panas pada teras serta distribusi temperaturnya.1 Materi yang diberikan dalam tulisan ini terkait dengan bagian pertama yaitu bagian neutronik. Pada tulisan ini akan dipaparkan mengenai formulasi metoda nodal pada geometri kartesian 3-D x-y-z2. Metoda nodal merupakan salah satu metoda untuk dapat memecahkan persamaan difusi neutron, metoda ini lebih unggul bila dibandingkan dengan metoda beda hingga yang lebih awal digunakan. Dengan metoda nodal maka diskritisasi dapat dilakukan dengan mesh yang lebih besar. Untuk mendapatkan hasil dengan akurasi yang sama, metoda beda memerlukan lebar mesh untuk tiap dimensi ~1cm, sedangkan metoda nodal ~20cm. Hal ini tentunya akan menghasilkan efisiensi perhitungan analisis reaktor nuklir yang lebih besar. Persamaan Kesetimbangan Nodal Berawal dari persamaan difusi neutron tunak beberapa grup berikut 1  ∇ ⋅ J g ( r ) + Σ tg ( r )φ g ( r ) = ∑  χ gυΣ fg ' ( r ) + Σ gg ' ( r ) φ g ' ( r ) + Q g ( r ) λ  J g ( r ) = − D g ∇φ g ( r ) g = 1,2,3,...G

[1] [2]

dimana J g (r)

φg ( r )

= Arus neutron pada grup energi g (cm-2s-1) = Fluks neutron skalar pada grup energi g (cm-2s-1)

Σ tg ( r )

= Penampang lintang total makroskopik untuk grup energi g (cm-1)

χg υ Σ fg ' ( r )

= Spektrum fisi untuk grup energi g = Jumlah rerata neutron yang dihasilkan dari satu reaksi fisi = Penampang lintang fisi makroskopik untuk grup energi g (cm-1) 1/9 Revisi Terakhir : 7/1/2009

Tim Neutronik 2009 Bidang Komputasi PPIN-BATAN Σ gg ' ( r ) Qg ( r )

= Penampang lintang transfer makroskopik dari grup energi g’ ke g (cm-1) = Sumber neutron luar untuk grup energi g (cm-3s-1)

Dg G

= Koefisien Difusi untuk grup energi g (cm) = Jumlah grup energi total

Persamaan.1 diatas merupakan hasil integrasi dari persamaan difusi terhadap grup energi diskrit. Parameter neutronik dari tiap grup seperti penampang lintang dan koefisien difusi diperoleh dari rerata bobot spektrum untuk suatu grup energi (a spectrum-weighted average over the energy group). Untuk penurunan formula lebih jauh pada geometri tiga dimensi x-y-z digunakan koordinat umum u,v,w . Ruang dibagi kedalam partisi geometri balok yang berisi material homogen dengan indeks partisi ul,vm,wn dimana  i = 1,2,..., I  l , m, n =  j = 1,2,..., J dimana I,J,K adalah jumlah partisi pada tiap dimensi x,y,z secara  k = 1,2,...K berurutan. Partisi (nodal) dengan indeks (i,j,k) didefinisikan oleh  x y

Є Є

[ xi , xi +1 ]

z

Є

[ zk , zk +1 ]

[y , y ] j

j +1

Dengan lebar tiap nodal adalah hul ≡u l +1 −u l dimana u=x,y,z sehingga volume dari tiap nodal adalah V i , j ,k ≡ hxi h yj hzk Maka dalam geometri kartesian 3-D pers.1 dinyatakan sebagai berikut ∂ ∂ ∂ J gx ( x, y , z ) + J gy ( x, y , z ) + J gz ( x, y , z ) + Σ tg ( x, y, z )φ g ( x, y, z ) = ∂x ∂y ∂z G 1  = ∑  ⋅ χ g ⋅ υ ⋅ Σ fg ' ( x, y , z ) + Σ gg ' ( x, y, z ) φ g ' ( x, y, z ) + Q g ( x, y, z )  g '=1  λ

serta J gu ( x, y , z ) = − D g ( x, y, z )

∂ φ g ( x, y , z ) ∂u

; u = x,y,z

[3]

[4]

2/9 Revisi Terakhir : 7/1/2009

Tim Neutronik 2009 Bidang Komputasi PPIN-BATAN Langkah yang dilakukan untuk memperoleh persamaan kesetimbangan nodal adalah dengan mengintegrasikan pers.3 terhadap seluruh volume nodal (i,j,k) lalu membaginya dengan volume nodal,Vi,j,k. Untuk integrasi suku kebocoran sebelumnya perlu diterapkan Hukum Gauss berikut ui +1

∫∇⋅ J

g

ui

⋅ dV i ,k , j = ∫ J g ⋅ dAi ,k , j

[5]

sehingga integrasi terhadap volume berubah menjadi integrasi terhadap seluruh permukaan yang membatasi nodal (i,j,k). Secara eksplisit untuk arah x dapat digambarkan sebagai berikut

J gxjk ( xi +1 )

J gxjk ( xi )

Xi

Xi+1

Gambar neutrondua arah x Sehingga integrasi arus neutron arah1.x Arus melibatkan permukaan pada xi+1 dan xi sebagai berikut z k +1 y j +1

∫ ∫ J ( x, y, z ) ⋅ dA

[

zk

]

= J gxjk ( xi +1 ) − J gxjk ( xi ) h yj hzk

jk

gx

[6]

yj

lalu hasil ini dibagi dengan volume nodal Vi,j,k sehingga untuk arah x diperoleh z k +1 y j +1

∫ ∫

zk

J gx ( x, y , z ) ⋅ dA jk =

yj

[

1 jk J gx ( xi +1 ) − J gxjk ( xi ) i hx

]

[7]

Dengan integrasi yang serupa untuk arah y dan z maka persamaan kesetimbangan nodal diperoleh sebagai berikut

[

]

[

]

[

]

1 jk 1 ik ( y j +1 ) − J gyik ( y j ) + 1k J gzij ( z j +1 ) − J gzij ( z j ) + Σ ijktg φ gijk J gx ( xi +1 ) − J gxjk ( xi ) + j J gy i hx hy hz G 1  ijk ijk = ∑  ⋅ χ g ⋅ υ ⋅ Σ ijkfg ' + Σ ijk gg ' φ g ' + Q g  g '=1  λ

[8]

dimana

φ

ijk g

Q

ijk g

1 ≡ ijk V

xi +1

y j +1

xi

yj

1 ≡ ijk V

xi +1

y j +1

xi

yj

∫ dx ∫ ∫ dx

∫

z k +1

dy ∫ dz ⋅ φ g ( x, y, z )

[9]

Zk

z k +1

dy ∫ dz ⋅ Q g ( x, y, z )

[10]

Zk

3/9 Revisi Terakhir : 7/1/2009

Tim Neutronik 2009 Bidang Komputasi PPIN-BATAN J

mn gu

1 (u ) ≡ m n hv hw

v m +1

wn +1

vm

wn

∫ dv ∫ dw ⋅ J ( u, v, w) gu

dengan u = x,y,z dan u ≠ v , w ≠ u,v

[11]

Parameter dengan indeks atas (i,j,k) menunjukkan nilai rerata untuk nodal (i,j,k), seperti fluks rerata nodal, sumber rerata nodal, dan penampang lintang rerata nodal. Sedangkan parameter dengan indeks atas (mn) menunjukkan nilai rerata bidang (mn), seperti arus neutron rerata bidang. Persamaan kesetimbangan neutron diatas memiliki tujuh variable yang tidak diketahui yaitu 6 arus neutron rerata-permukaan (face-average neutron current) pada tiap permukaan nodal, dan satu fluks neutron rerata-nodal (node-average neutron flux). Untuk itu diperlukan 6 persamaan lain yang mengaitkan antara arus neutron rerata-permukaan dengan fluks neutron rerata-nodal. Persamaan tambahan ini biasa disebut dengan persamaan kopling (coupling equations). Untuk mendapatkan persamaan kopling ini akan dilakukan dua cara yaitu metoda beda-hingga terkoreksi dan metoda nodal polinomial.

A. Pendekatan Beda Hingga Terkoreksi Salah satu cara untuk memperoleh persamaan kopling yang diperlukan adalah dengan pendekatan beda hingga. Dimulai dengan mengintegrasikan persamaan yang menghubungkan antara arus neutron dan fluks neutron (Hk. Ficks) dalam nodal dan membaginya dengan volume nodal mn J gu (u ) ≡ −

D glmn d vm +1 wn +1 dv dw ⋅ φ g ( u , v, w) hvm hwn du v∫m w∫n

[12] dengan u = x,y,z ; u Є [ul,ul+1] lalu turunan pada pers.12 diatas dapat didekati dengan persamaan beda (difference equation) sederhana berikut J

mn gu

(u l ) ≈ − D

φ

mn gu

1 (u ) ≡ − m n hv hw

lmn g

φ glmn − φ gmn (u l+ )

[13] hlu 2 dimana fluks rerata-bidang yang muncul pada pers.13 didefinisikan sebagai berikut wn +1 U Ul l-1 ( ) dv dw ⋅ φ u , v , w g ∫ ∫

v m +1

vm

Ul+1

[14]

wn

+

dan ul menunjukkan sisi positif dari +permukaan pada indeks l sebagaimana terlihat pada - perm. l-1nodal l gambar.2 4/9 Konvensi arahRevisi pada Terakhir : 7/1/2009 bidang batas

Gambar 2. Konvensi indeks pada permukaan dan nodal

Tim Neutronik 2009 Bidang Komputasi PPIN-BATAN

Perlu diperhatikan bahwa pembentukan persamaan beda sebagaimana pada pers.13 berarti mengasumsikan bahwa fluks pada nodal berubah secara linier dari permukaan lmn + nodal (dimana fluks bernilai φ g (ul ) ) hingga ke tengah nodal (dimana fluks bernilai

φ glmn ) . Cara yang sama dapat dilakukan pula untuk nodal dengan indek l-1 J

mn gu

(u l ) ≈ − D

lmn g

φ gmn (u l− ) − φ glmn [15]

hlu −1 2

Namun pers.13 dan pers.15 hanya akurat untuk ukuran nodal yang kecil dan dapat memberikan eror yang besar ketika satu perangkat bahan bakar (~20cm) dijadikan satu nodal/partisi. Pada kasus reaktor berpendingin air ringan (PWR atau BWR) metoda beda hingga memerlukan lebar partisi nodal sebesar satu pin bahan bakar (~1cm) untuk menghasilkan hasil yang konvergen. Kekurangan ini diatasi dengan menambahkan faktor koreksi sehingga memaksakan pers.13 dan pers.15 bersifat eksak∴ secara formal. Hal ini dicapai dengan mengalikan pers.13 dan pers.15 dengan faktor koreksi untuk memberikan fluks rerata-permukaan yang sebenarnya pada permukaan tersebut

φ gumn (u l ) = f gulmn− ⋅ φ gumn (u l+ ) = f gul −+1,mn ⋅ φ gumn (u l− )

[16]

karena kedua faktor koreksi diatas umumnya memiliki nilai yang berbeda maka fluks rerata yang terdapat pada pers.13 dan pers.15 haruslah bersifat diskontinu. Untuk alasan inilah faktor koreksi diatas disebut faktor diskontinu. (discontinuity factors). Dengan memasukkan faktor diskontinu ke pers.13 dan pers.15 maka diperoleh

φ mn J gu (u l ) ≡ − D glmn

lmn g

−

φ gmn (u l ) f l u

h 2

lmn gu −

φ gumn (u l ) ≡ − D gl −1,mn

f

l −1, mn gu +

− φ gl −1,mn

hul −1 2

[17]

cek http://www.sosmath.com/diffeq/first/exact/exact.html untuk definisi persamaan eksak atau tidak eksak.

5/9 Revisi Terakhir : 7/1/2009

Tim Neutronik 2009 Bidang Komputasi PPIN-BATAN Faktor diskontinu yang dikenalkan diatas adalah untuk mengkoreksi kesalahan formulasi beda hingga spasial (spatial difference errors). Dalam konsep General Equivalence Theory yang diberikan oleh K.Smith dalam disertasinya tahun 1980 adalah untuk mengkoreksi kesalahan yang muncul karena memperlakukan suatu daerah yang heterogen seakan memiliki komposisi yang homogen (proses homogensisasi). Lebih jauh lagi, faktor diskontinu dapat pula digunakan untuk koreksi terhadap kesalahan aproksimasi teori difusi (sebagai pendekatan terhadap transport neutron) juga kesalahan dalam memperoleh koefisien difusi. Sehingga faktor diskontinu pada akhirnya bermanfaat untuk koreksi dari sisi spasial, homogensisasi, dan juga difusi. Dengan apa yang telah diperoleh diatas, sekarang kita coba untuk memperoleh persamaan kopling yang dicari. Menggunakan pers.16 yang menggambarkan kontinuitas pada permukaan (setelah ada faktor diskontinu) kita dapa eliminasi fluks rerata-bidang dari pers.17 sehingga diperoleh hubungan antara arus rerata-permukaan dengan fluks reratanodal berikut J

mn gu

 hl f gulmn− hul −1  (u l ) ≡ −  ulmn l −1,mn +  2 D gl −1,mn   2 D g f gu +

−1

 f gulmn−   φ glmn − φ gl −1,mn  l − 1 , mn f   gu + 

[18]

Dengan melakukan manipulasi yang sama dapat diperoleh pula arus neutron pada batas nodal ul+1 berikut J

mn gu

 hul f gulmn+ hul +1  (u l +1 ) ≡ −  lmn l +1,mn +  2 D gl +1,mn   2 D g f gu −

−1

 f gulmn+  lmn l +1, mn   φ − φ g g  f l +1,mn   gu − 

[19]

dengan melakukan substitusi pers.18 dan pers.19 untuk arah x,y,z ke dalam persamaan kesetimbangan nodal pada pers.8 diperoleh persamaan nodal dalam bentuk beda-hingga berikut 1 hxi

 hxi f gxijk− hxi −1   ijk i −1, jk +  2 D gi −1, jk   2 D g f gx +

−1

 f gxijk− ijk  i −1, jk   φ − φ g  f i −1, jk g   gx +  −1

1 + i hx

 hxi f gxijk+ hxi +1   ijk i +1, jk +  2 D gi +1, jk   2 D g f gx −

1 + j hy

 h yj f gyijk− h yj −1   ijk i , j −1,k +  2 D gi , j −1,k   2 D g f gy +

−1

 f gyijk−  ijk i , j −1, k   φ − φ g  f i , j −1,k g   gy + 

1 + j hy

 h yj f gyijk+ h yj +1   ijk i , j +1,k +  2 D gi , j +1,k   2 D g f gy +

−1

 f gyijk+  ijk i , j +1, k   φ − φ g  f i , j +1,k g   gy + 

1 + k hz

 hzk f gzijk− hzk −1   ijk ij ,k −1 +  2 D gij ,k −1   2 D g f gz +

−1

 f gxijk+ ijk  i +1, jk   φ − φ g  f i +1, jk g   gx − 

[20]

 f gzijk− ijk  ij , k −1   φ − φ g  f ij ,k −1 g   gz +  6/9 Revisi Terakhir : 7/1/2009

Tim Neutronik 2009 Bidang Komputasi PPIN-BATAN −1

ijk  hzk  f gzijk+ hzk +1   f gz + ijk ij , k +1  φ − φ  ijk ij ,k +1 +  g g  2 D gij ,k +1   f gzij ,+k +1  2 D g f gz +  G 1  ijk ijk + Σ ijk = ∑  χ g Σ ijkfg ' + Σ ijk + Q gijk tg φ g gg '  ⋅ φ g  g '=1  λ

1 + k hz

Pers.20 dapat dituliskan dalam notasi matrik sebagai berikut N gφ g =

[

]

[

]

G 1 G Fgg ' ⋅ φ g ' + ∑ Σ gg ' ⋅ φ g ' + Q gijk ∑ λ g '=1 g '=1

[21]

g '≠ g

dimana Ng

φg Σ gg '

= Matriks pentadiagonal N x N yang mengandung suku kopling untuk grup g, penampang lintang total, dan tumbukan dalam-grup. = Vektor kolom orde N yang mengandung fluks neutron skalar pada grup energi g = Matriks diagonal N x N yang mengandung nilai {Σ gg ' }

Fgg '

Qg ( r )

= Matriks diagonal N x N yang mengandung nilai {χ g ⋅ υ ⋅ Σ lmn fg ' } = Vektor kolom berorde N yang mengandung suku sumber luar

N

= Jumlah nodal total = I x J x K

Bentuk matriks lainnya dari pers.20 dapat pula sebagai berikut Aφ = dimana A

1 Mφ + Q λ

[22]

M

= Matriks NG x NG yang mengandung { N g δ gg ' − Σ gg ' } = Vektor orde N yang mengandung fluks neutron skalar pada grup energi g. Kolom { φ g } = Matriks diagonal NG x NG yang mengandung nilai {Fgg ' }

Q

= Vektor orde NG, representasi dari sumber luar. Kolom {Qg}

φ

Syarat Batas Syarat batas yang diterapkan adalah sebagai berikut mn ( u l ) lˆ ⋅ nˆ φ gumn ( u l ) = Γgumn± J gu

[23]

dimana 7/9 Revisi Terakhir : 7/1/2009

Tim Neutronik 2009 Bidang Komputasi PPIN-BATAN

φ gumn ( u l )

= Rerata fluks-permukaan pada batas

mn ( ul ) J gu

= Arus rerata-permukaan pada batas

ul lˆ nˆ mn Γgu±

= = = =

Batas eksternal Vektor sauan pada arah positif sumbu koordinat Vaktor satuan dari batas eksternal Faktor syarat batas dengan nilai sebagai berikut :

Γgumn± = 0

 fluks nol

Γgumn± = 2

 arus masuk nol

Γgumn± = ∞

 arus nol

Γgumn± = 2 +

4

4 mn  albedo dimana Γgu ± = 2 + α −1 − 1

α −1 g Nilai dari arus neutron pada permukaan terluar (eksternal) diperoleh dari pers.17 dan pers.23 dengan mengeliminasi fluks rerata-permukaan sehingga diperoleh nilai arus neutron pada permukaan terluar bawah (indeks lebih kecil) adalah −1 g

−1

J

mn gu

 Γgulmn− hul  (u l ) ≡ −  lmn + φ glmn lmn   f gu − 2 D g 

[24]

sedangkan untuk permukaan terluar atas (indeks lebih besar) adalah −1

 Γgulmn+ hul  mn J gu (u l ) ≡ −  lmn + φ glmn lmn  f 2 D g  gu + 

[25]

Perhitungan Faktor Diskontinu Dari pers.18 bisa didapatkan hubungan untuk faktor diskontinu berikut

f gul −+1,mn f gulmn−

φ glmn + =

φ

l −1, mn g

hul mn J gu (u l ) 2 D glmn

hul −1 mn − J gu (u l ) 2 D gl −1,mn

[26]

Terkait dengan pendefinisian faktor diskontinu maka ada hal menarik pada syarat batas, untuk itu kita coba ubah pers.24 dan pers.25 ke dalam bentuk lain berikut lmn Γgu −

f gulmn−

 φ glmn hl = − mn + ulmn  J (u ) 2 D g  gu l

   

[27] 8/9 Revisi Terakhir : 7/1/2009

Tim Neutronik 2009 Bidang Komputasi PPIN-BATAN Γgulmn+ f gulmn+

 φ glmn hl =  mn + ulmn  J (u ) 2 D g  gu l

   

[28]

mn Terlihat bahwa apabila ditetapkan syarat batas fluks nol, maka Γg bernilai nol namun

perbandingan

Γgumn−

f gulmn−

tidak harus bernilai nol. Dari sini terlihat bahwa dengan

memperkenalkan faktor diskontinu kita telah memberikan syarat batas yang berbeda sehingga memperhitungkan kebocoran neutron pada batas dengan baik. Referensi : 1. Christensen B, 1985, “Three Dimesional Static and Dynamic Reactor Calculations by the Nodal Expansion Method”, Riso National Laboratory RISOR-496. 2. Gehin JC, 1992, ”Quasi Static Polynomial Nodal Method for Nuclear Reactor Analysis”, MIT Dissertation.

9/9 Revisi Terakhir : 7/1/2009