Solusi Numerik Persamaan Difusi Neutron

http://syeilendrapramuditya.wordpress.com

Solusi Numerik Persamaan Difusi Neutron III.1 Pendahuluan

Meto Metode de atau atau tekn teknik ik peme pemeca caha han n pers persam amaan aan mate matema mati tiss terb terbag agii menj menjad adii dua dua golongan besar, yaitu metode analitik dan numerik. Solusi yang dihasilkan dengan metode analitik adalah solusi yang sesungguhnya, sebenarnya, dan juga eksak (exact ), ), sedangkan solusi numerik adalah aproksimasi atau pendekatan dari solusi sebena sebenarny rnya, a, dengan dengan orde orde error error tertent tertentu. u. Beberap Beberapaa persam persamaan aan matemat matematis is pada pada persoalan-persoalan fisika pada kenyatannya relatif sangat sulit untuk dipecahkan secara secara analiti analitik, k, karena karena itulah itulah dikemb dikembang angkan kan metode metode numeri numerik k untuk untuk mencari mencari solusinya.

Metode numerik yang pertama-tama akan dibahas disini adalah teknik pemecahan persamaan persamaan diferensial diferensial dengan dengan menggunaka menggunakan n aproksimas aproksimasii untuk untuk fungsi fungsi turunan turunan pertama dan turunan kedua, karena sebagaimana telah diketahui, sebagian besar persamaan-persamaan dalam fisika adalah berupa persamaan diferensial.

Bila terdapat suatu fungsi sembarang sembarang f ( x ) yang akan dicari turunannya, turunannya, yaitu f ' ( x ) dan f ' ' ( x) , maka pertama-tama kita akan menuliskan ekspansi deret

Taylor untuk f ( x + h) dan f ( x − h) sebagai berikut :

f ( x

+ h) = f ( x) + h f ' ( x) +

h

2

2!

f ' ' ( x )

+

h

−

h

3

3!

f ' ' ' ( x )

+

h

+

h

4

f ( x )

+

f iv ( x )

+

iv

4!

(III.1) f ( x

−h) = f ( x ) − h f ' ( x) +

h

2

2!

f ' ' ( x )

3

3!

f ' ' ' ( x)

4

4!

(III.2)

Kemudi Kemudian an untuk untuk mendap mendapatk atkan an

f ' ( x ) ,

pers persam amaan aan (III. (III.1) 1) diku dikura rang ngii oleh oleh

persamaan (III.2) : f ' ( x) =

f ( x + h) − f ( x − h) 2h

−

h2 3!

f ' ' ' ( x) −

(III.3)

1

h4 5!

f ( x) + v


Persamaan (III.3) diatas bukanlah aproksimasi, tetapi e kspresi eksak dalam bentuk deret Taylor dari turunan pertama. Bila kita mengabaikan semua suku selain suku pertama, maka kita akan memperoleh persamaan berikut :

f ' ( x) ≡

d dx

f ( x ) ≅

f ( x + h) − f ( x − h) 2h

(III.4)

Persam Persamaan aan (III.4 (III.4)) diatas diatas adalah adalah aproks aproksima imasi si numerik numerik untuk untuk fungs fungsii turuna turunan n pertama, dan suku lain yang diabaikan dianggap sebagai suku error.

Untuk mendapatkan f ' ' ( x) , persamaan (III.1) ditambahkan dengan persamaan (III.2) :

f ' ' ( x)

=

f ( x + h)

− 2 f ( x ) + f ( x − h)

h

2

−

h2 12

iv

f ( x)

−

h4 360

f vi ( x )

+

(III.5)

Persamaan (III.5) diatas bukanlah aproksimasi, tetapi e kspresi eksak dalam bentuk deret Taylor dari turunan kedua. Bila kita mengabaikan semua suku selain suku pertama, maka kita akan memperoleh persamaan berikut :

f ' ' ( x) ≡

d 2 dx

2

f ( x) ≅

f ( x + h) − 2 f ( x ) + f ( x − h) h2

(III.6)

Persamaan (III.6) diatas adalah aproksimasi numerik untuk fungsi turunan kedua, dan suku lain yang diabaikan dianggap sebagai suku error. Persamaan (III.4) dan (III.6) diatas dikenal sebagai aproksimasi beda hingga [12] ( finite finite difference).

III.2 Pemecahan Numerik Persamaan-Persamaan Neutronik

Pada perhitungan neutronik teras reaktor, terdapat dua proses utama yang harus dikerjakan, yaitu :

2


1. Pertama Pertama adalah perhitu perhitungan ngan distrib distribusi usi fluks fluks neutron neutron di dalam dalam teras teras reaktor, reaktor, yaitu yaitu dengan dengan cara cara memeca memecahka hkan n persam persamaan aan (II.52) (II.52),, persam persamaan aan difusi difusi multig multigrup rup.. Setela Setelah h distri distribus busii fluks fluks neutro neutron n diketa diketahui hui,, maka maka besaran besaran- besaran lain seperti distribusi kerapatan daya dan suku sumber juga dapat dihitung. Proses ini dikenal sebagai iterasi dalam ( inner iteration ). 2. Prose rosess ked kedua adal adalah ah perh erhitun itunga gan n

kriti ritika kali lita tass

tera terass reak reakto tor, r, yait yaitu u

perhi perhitun tungan gan nilai nilai faktor faktor multip multiplik likasi asi teras. teras. Perhit Perhitung ungan an ini dilaku dilakukan kan dengan menggunakan beberapa persamaan neutronik. Proses ini dikenal sebagai iterasi luar (outer iteration ).

III.2.1 Perhitungan distribusi fluks neutron : iterasi dalam Seperti Seperti telah telah dijelas dijelaskan kan sebelu sebelumny mnya, a, perhit perhitung ungan an distri distribus busii fluks fluks dikerja dikerjakan kan dengan menggunakan persamaan difusi multigrup (II.52) :

r ∂φ g (r , t ) r − ∇ ⋅ D (g vr) ∇φ ∂t v g

1

v

v

v

(g r, t) + Σ t(g r)φ (g r, t) =

G

∑Σ

v

v

(r )φ g' (r , t ) +

s'g g

g '=1

χ g

k eff

G

∑v

(II.52)

Σ

g'

v

v g'

( r )φ (r , t )

fg'

g ' =1

Bila kita meninjau teras pada keadaan tunak ( steady state ), maka variabel variabel waktu waktu dapat diabaikan, dan dengan definisi bahwa material pada setiap region teras adalah homogen, maka persamaan (II.52) akan berbentuk :

v

− D g( r)∇ φ 2

v

v

v

g( r) + Σ tg( r)φ g( r ) =

G

∑Σ

g '=1

v

S (r ) =

G

∑v Σ

g '=1

g'

v

v

v

s'g g( r)φ 'g( r ) +

χ g

k eff

v

S( r )

v

(r )φ g' ( r )

(III.7)

(III.8)

fg' g

Persamaan (III.7) diatas memiliki syarat batas φ g(r S)

= S (r S) = 0 , yaitu fluks dan

suku sumber pada permukaan teras reaktor harus bernilai nol.

3


Teras reaktor yang ditinjau memiliki geometri silinder silinder dua dimensi R-Z, dimana dimana geomet geometri ri ini selanj selanjutn utnya ya dibuat dibuat menjad menjadii diskri diskritt dengan dengan cara cara dibagi dibagi menjad menjadii beberapa partisi radial

∆r

dan aksial

∆ z .

Dengan Dengan demikian demikian,, nilai nilai fluks fluks yang yang

didapa didapatka tkan n nanti nanti tidakl tidaklah ah kontin kontinyu yu di setiap setiap bagian bagian teras, teras, melain melainkan kan berupa berupa distribusi diskrit di titik-titik tertentu.

Radial Partition

Axial Partition

∆ z

∆r Gambar III.1 Partisi geometri silinder teras reaktor

Pada sistem diskrit, vektor posisi dinyatakan dengan cara berikut : r = r , z ≡ i, j

(III.9)

Operator Laplacian Laplacian (persamaan III.7) pada geometri silinder adalah :

∇2 =

1 ∂

2 ∂ ∂2 1 ∂ + 2 + 2 r ∂r ∂r ∂ z 2 r ∂φ

r

(III.10)

Fluks tidaklah bergantung pada sudut azimut, maka dengan menggunakan prinsip simetri, persamaan (III.10) menjadi lebih sederhana : ∂2 ∂2 1 ∂ ∇ = + + r ∂r ∂r 2 ∂ z 2 2

(III.11)

Substitusi persamaan (III.11) ke persamaan (III.7) :

4


∂ 2φ (gr , z ) 1 ∂φ (gr , z ) ∂2 φ (gr , z ) Σ t(gr , z ) φ (gr , z ) + + − = ∂ r2 ∂r ∂ z2 r D g ( r, z) G  χ g 1 −  ∑ Σ sg' g(r , z)φ g' ( r , z ) + S (r , z )  D( ,r )z  g '=1 g k  eff

(III.12)

Syarat batas untuk persamaan (III.12) diatas adalah :

•

Fluks pada permukaan teras bernilai nol : φ (r S , z ) =φ (r , z S ) = 0

•

Suku

sumber

pada

permukaan

teras

bernilai

nol

:

S (r S , z ) = S (r , z S ) = 0

Persamaan (III.4) dan (III.6) dapat ditulis dalam bentuk berikut : d

f ( x ) ≅

dx

f ( x + h) − f ( x − h) 2h

=

f i +1 − f i −1 2∆ x

(III.13)

d 2

f ( x + h) − 2 f ( x) + f ( x − h)

dx

h2

f ( x) ≅ 2

=

f i +1 − 2 f i + f i −1 ∆ x

2

(III.14)

Dengan menggunakan persamaan (III.13) dan (III.14) diatas, serta definisi vektor posisi pada persamaan (III.9), maka bentuk diskrit dari persamaan (III.12) adalah sebagai berikut :

− 2φg , i , j + φg , i −1, j 1 φ g , i +1, j − φg ,i −1,j φg ,i ,j + 1 − 2φ g ,i ,j + φg ,i ,j − 1 + + − 2 ∆r ∆r 2 i ∆r ∆z 2 (III.15)  Σtg , i , jφ g , i , j χ g 1  G =−  ∑ Σ( sg ' g), i, jφ g ', i , j + S i , j  D g, i, j D g , i, j  g '=1 k eff 

φ g , i +1, j

+ φg , i −1, j 1 φ g , i +1, j − φg , i −1, j φg , i ,j +1 + φg , i ,j − 1 + + − 2 ∆r ∆r 2 i ∆r ∆z 2   Σt g, i , j 2 χ g 2  1  G + 2 + 2  Σ + φ g, i, j = − φ S   ∑ ( sg ' g ), i , j g ', i , j i, j   D ∆ ∆ D r z k  g, i, j  eff  g, i, j   g '=1 φg , i +1, j

5

(III.16)


 Σt g, i, j 2 χ g n  2  n +1 1 G n φ φ + + = Σ +    ∑ ( sg ' g ), i, j g , i , j n S i , j  + 2 2  g, i, j D D ∆ ∆ r z k eff  g, i, j   g, i, j   g '=1 φ gn, i +1, j + φgn, i −1, j φ gn, i +1, j − φgn , i −1, j φgn , i ,j +1 + φgn , i ,j − 1 + + i 2∆r 2 ∆r 2 ∆z 2

(III.17)

Superskrip n pada persamaan (III.17) diatas menunjukan nilai awalnya, sedangkan superskrip n+1 menunjukan nilai barunya setelah dihitung secara iteratif.

φ ng+,1i, j

 1 G  χ g n     ∑ Σ( sg ' g), i, jφ g ',i , j + n S i , j  + k eff  D g, i, j  g '=1    n  n n n n n φ φ φ φ φ φ + − +  g , i+1, j g , i −1, j + g , i +1, j g , i −1, j + g , i ,j +1 g , i ,j −1    ∆ ∆ r2 i 2∆r 2 z2   =  Σ t g, i , j 2 2  + 2 + 2     D g, i, j ∆r ∆z 

(III.18)

Perhat Perhatika ikan n bahwa bahwa suku suku kedua kedua pada pada ruas ruas kiri kiri persam persamaan aan (III.12 (III.12)) mengan mengandun dung g (1 / r ) , maka pada r

= 0 atau i = 0

suku ini akan bermasalah karena akan bernilai

tak hingga, sehingga persamaan (III.18) diatas hanya akan berlaku untuk nilai

r

≠ 0 atau i ≠ 0 .

Untuk r

= 0 atau i = 0 ,

maka maka suku suku yang yang meng mengan andu dung ng (1 / r ) tersebut tersebut harus

ditangani secara khusus, yaitu sebagai berikut [12] :

Berdasarkan teorema limit L’Hospital [13] : 1 ∂φ r ∂r

g r =0

=

∂φ / ∂gr r

= r =0

lim

r → 0

∂ ∂ ∂ ( φ / r ) ∂ 2φ ∂r g = 2 ∂ ∂r (r ) ∂r

g r = 0

(III.19)

Substitusikan persamaan (III.19) ke persamaan (III.12) :

∂ 2φ (gr , z) ∂2φ (gr , z ) Σ (tgr , z ) φ (gr , z ) 2 + − = ∂ r2 ∂ z2 D g ( r, z) G χ g 1 φ ( , ) ( , ) − Σ + r z r z  ∑ sg' g g' D( ,r )z  g '=1 g k

6

 S (r , z )   eff

(III.20)


φ g , 0+1, j − 2φ g , 0, j + φ g , 0−1, j

∆r 2 Σtg , 0, jφ g , 0, j ,D 0g,

1

=−

D

j

, 0g,

+

φ g , 0, j +1 − 2φ g , 0, j

+ φg , 0, j −1

∆z 2

G  χ g Σ + φ S  ∑ ( sg' g), 0, j g', 0, j 0, j  keff  j g '=1

− (III.21)

Berdasarkan simetri sudut azimut pada geometri silinder yang ditinjau :

φ g, i+1, j

= φ g, i−1, j

(III.22)

Dengan menggunakan definisi (III.22), maka persamaan (III.21) akan berbentuk :

2

2φ g, 1,

∆r 2

j

+

φ g, 0,

   φ g , 0, j =   χ g 1  G − + S  ∑ Σ ( sg' g), 0, φ j g',0, j 0, j D , 0, g  g 'j=1 k eff  j+1

+ φ g, 0, j−1  Σ tg, 0, j 4 2 − + +   D g ∆r 2 ∆z2 ∆z 2 

 Σtg , 0, j 4 2  n+1 φ g,0, j = + 2 + 2    D r z ∆ ∆  , 0, g j 

G  χ g  ∑ Σ( sg' g), 0, jφ g',0, j + n S 0n, j  + D, 0, g  g j '=1 k eff  4φ ng, 1, j φ ng, 0, j+1 + φ ng, 0, j−1 + ∆r 2 ∆z 2

(III.23)

1

G χ g n  4φ ng, 1, j φ ng, 0, j+1 + φ ng, 0, j−1 + S 0, j +  ∑ Σ( sg' g), 0, φ j g', 0 , j + n 2 ∆ ∆ 2z D k r g ' 1 =   , 0, g  j eff  n +1 φ g ,0, j =  Σtg , 0, j 4  + 2 + 2 2     D g , 0, j ∆ r ∆ z 

(III.24)

1

(III.25)

Bila dilakukan iterasi terhadap persamaan diskrit (III.18) dan (III.25) diatas, maka pada akhirnya akan tercapai keadaan konvergen numerik dengan akurasi atau orde error tertentu, yaitu :

7


φ ng+,i, j

−

1

φ ng, i, j

φ ng, i , j

< ε ,

untuk untuk selu seluru ruh h g, i, j

(III.25a)

Metode numerik yang digunakan untuk menurunkan persamaan diskrit (III.18) dan (III.25) adalah metode yang disebut Jacobi Iteration Scheme [12], atau metode

Jacobian . Laju konvergensi metode Jacobian sebenarnya tidaklah terlalu tinggi, karena pada meto metode de ini, ini, nila nilaii fluk flukss yang yang baru baru dida didapa patt dari dari hasi hasill perh perhit itun unga gan n deng dengan an menggunakan nilai fluks yang lama seluruhnya, ata u disebut layer-by-layer .

φ k n 21

φ k n 11

φ k n

φ k n 2

φ k n−1

φ k n 21

φ k n 11

+

−

n o i t a r e t i

−

−

−

+

−

−

−

φ k n 11

φ k n 21

φ k

φ k n+1

φ k

φ k n

φ k n 11

φ k n 21

+1

+

+

n

−1

−

+

+

+

n +2

−

+

points points sequence

Gambar III.2 Skema iterasi Jacobian

Gambar Gambar (III.2) (III.2) diatas adalah skema iterasi Jacobian, dimana nilai fluks pada tiap layer dihitung hanya dengan menggunkan nilai fluks pada layer sebelumnya.

Laju konvergensi dapat ditingkatkan dengan menggunakan metode Gauss-Siedel pada meto metode de ini, ini, nila nilaii fluk flukss yang yang baru baru dihi dihitu tung ng deng dengan an Iterati Iteration on Scheme Scheme, pada memanf memanfaatk aatkan an secara secara langsu langsung ng nilai nilai fluks fluks yang yang baru baru saja saja dihitu dihitung ng pada pada titik titik partisi sebelumnya.

8


φ k n 21

φ k n 11

φ k n

φ k n 2

φ k n−1

φ k n 21

φ k n 11

+

+

−

n o i t a r e t i

−

−

−

−

−

−

φ k n 11

φ k n 21

φ k n

φ k n+1

φ k n 2

φ k n

φ k n 11

φ k n 21

+1

+

+

−1

−

+

+

+

+

−

+

points points sequence

Gambar III.3 Skema iterasi Gauss-Siedel Dengan Dengan menggunak menggunakan an skema skema iterasi iterasi Gauss-Sied Gauss-Siedel, el, maka persamaan persamaan (III.18) (III.18) dan (III.25) akan berbentuk : untuk

φ ng+,1i, j

untuk

r ≠ 0 :

 1 G  χ g n  φ S Σ + +    ∑ ( sg ' g ), i, j g ',i , j k n i , j  D  , g, i j g '=1   eff   n  n +1 n n +1 n n +1 φ φ φ φ φ φ + − +  g , i+1, j g , i −1, j + g , i +1,j g , i −1,j + g , i ,j +1 g , i ,j −1    r2 i 2∆r 2 z2  ∆ ∆  =  Σt g, i, j 2 2  + 2 + 2   ∆ r ∆ z  D g , i , j  

(III.26)

r = 0 :

G χ g n  4φ ng, 1, j φ ng, 0, j+1 + φ ng,+10, j−1 + n S 0, j + +  ∑ Σ( sg' g), 0, φ j g', 0 , j 2 ∆ ∆ 2z D k r = g ' 1   , 0, g j e f f   n +1 φ g ,0, j =  Σtg , 0, j 4  + 2 + 2 2     D g , 0, j ∆ r ∆ z  1

(III.27)

Pada persamaan (III.26) dan (III.27) diatas, nilai konstanta grup atau cross section merupakan fungsi posisi, yang dilambangkan dengan subskrip i dan j. Sebenarnya nilai nilai konsta konstanta nta grup grup atau cross cross sectio section n bergan bergantun tung g pada pada jenis jenis dan kompos komposisi isi material di titik tersebut, maka untuk lebih menyederhanakan bentuk persamaan, subskrip ganda i dan j tersebut akan diganti dengan subskrip tunggal m, yang menunjukan jenis material di titik i,j tersebut.

9


Maka persamaan (III.26) dan (III.27) akan berbentuk :

r ≠ 0 :

untuk

φ ng+,1i, j

 1 G  χ g ,m n     ∑ Σ( sg' g), mφ g ', i, j + n S i, j  + k eff   D , g g m'=1   n  n +1 n n +1 n n +1  φ g , i+1, j + φg , i −1, j + φ g , i +1,j − φg , i −1,j + φg , i ,j +1 + φg , i ,j −1    ∆ ∆ r2 i 2∆r 2 z2   =  Σtg , m 2  + 2 + 2 2     D g , m ∆ r ∆ z 

r = 0 :

untuk

1 φ gn +,0,1 j

(III.28)

=

D,

G χ g , m n  4φ ng, 1, j φ ng, 0, j+1 + φ ng,+10, j−1 ∑ Σ( sg' g), mφ g',0, j + n S 0, j + ∆ 2 + ∆ 2z k eff  r g  g m'=1  Σtg , m 4  + 2 + 2 2     D g , m ∆ r ∆ z 

(III.29)

Laju konvergensi dapat ditingkatkan lagi dengan menggunakan metode Succesive

Over Relaxation (SOR). Prinsip dasar metode SOR adalah bahwa hasil konvergen sama dengan hasil iterasi ke n+1 ditambah selisih dari dua hasil iterasi terbaru, yaitu iterasi ke n dan ke n+1, yang dikalikan suatu konstanta.

Secara matematis, persamaan SOR dituliskan sebagai berikut : φ gn,+i1, j

= φ gn,+i1, j + α (φ gn,+i1, j − φ gn,i , j )

(III.30)

Pada persamaan (III.30) diatas, α adalah konstanta akselerasi SOR, yang nilainya bersifat unik untuk setiap kasus/persamaan.

Pada persamaan (III.28) dan (III.29), subskrip g menunjukan index grup, dan subskrip m menunjukan index jenis material, maka sekarang kita telah memiliki satu set lengkap persamaan difusi yang mampu menangani teras reaktor dengan spektrum energi neutron diskrit dan komposisi material heterogen, atau disebut

persamaan difusi multigrup-multiregion .

10


Persamaan Persamaan difusi difusi multigrup multigrup tersebut akan dipecahkan secara numerik numerik dengan dengan menggu menggunak nakan an metode metode SOR, SOR, yaitu yaitu dengan dengan mengg mengguna unakan kan persam persamaan aan (III.28 (III.28), ), (III.29), dan (III.30).

III.2.2 Perhitungan kritikalitas teras reaktor : iterasi luar Pada bagian sebelumnya telah dibahas mengenai teknik yang digunakan untuk memecah memecahkan kan persam persamaan aan difusi difusi multig multigrup rup,, yaitu yaitu untuk untuk menghi menghitun tung g nilai nilai dan distribusi distribusi fluks neutron di dalam teras reaktor. Selanjutnya, Selanjutnya, pada bagian bagian ini akan dibahas perhitungan sumber neutron S ( source source ) dan faktor multiplikasi k . Untuk Untuk memper mempermud mudah ah penuru penurunan nan persam persamaan aan,, pertam pertama-ta a-tama ma akan akan diguna digunakan kan persamaan difusi satu grup (II.48) : 1 ∂φ

v ∂t













−∇⋅ D ( r ) ∇ φ ( r , t ) + Σa (r ) φ (r , t ) = S ( r , t )

(II.48)

Untuk keadaan tunak dan teras homogen :

− D∇ 2φ ( vr) + Σaφ ( vr) =

1

k

v

S( r)

(III.31)

S (r ) = vΣ f φ (r )

(III.32)

Dengan menggunakan operator, persamaan (III.31) berbentuk : 

M φ ( r ) =

1

k

F φ (r )

(III.33)

M ≡ destructio n operator = (− D ∇2 + Σa )

(III.34) F ≡ production

operator

= vΣ f

Bentuk diskrit persamaan (III.31) adalah : 2

− D ∇

n

n

φ i , j + Σaφ i , j =

1 k

n

S ni , j

(III.36)

11

(III.35)

http://syeilendrapramuditya.wordpress.com n +1 Solusi dari persamaan (III.36) diatas adalah nilai fluks yang baru, yaitu φ i , j .

Selanjutnya nilai fluks yang baru ini digunakan untuk menghitung nilai source yang baru :

n +1

n +1

n +1

S i , j = F φ i , j = vΣ f φ i , j

(III.37)

Setelah nilai source yang baru diketahui, selanjutnya nilai k dapat dihitung dengan menggunakan persamaan (III.33) :

n +1

M φ



(r )

=

1 k

n +1

n +1

F φ

(r )

(III.38)

Integralkan persamaan (III.38) terhadap volume teras :

∫ d

3



1



1 r M φ n + ( r ) =

Vcore

k

∫ d

n +1

3



r F φ n + (r ) 1

Vcore

(III.39)

∫ d r F φ ∫ d r M φ 3

k n +1 =

n +1



( r )

Vcore

3



n +1

(III.40)

( r )

Vcore

∫





∫

d 3 r S n +1 ( r )

k

n +1

=

Vcore

1 k

n

∫

3 n





= k

n Vcore 3 n

∫ d r S

d r S ( r )

Vcore



d 3 r S n +1 ( r ) ( r )

Vcore

(III.41)

Operasi Operasi integr integral al

S ( r )

terh terhad adap ap volu volume me teras teras reakt reaktor or sebe sebena narn rnya ya untu untuk k

menghitung populasi source total di dalam teras reaktor. Karena S merupakan fung fungsi si posi posisi si,, maka maka untu untuk k meng menghi hitu tung ngny nya, a, akan akan lebi lebih h muda mudah h bila bila kita kita menggunakan S rata-rata rata-rata yang kemudian dikalikan dengan volume teras reaktor :

12


∫ d

3

r S ( r )

= S AVE V CORE

Vcore

(III.42)

Maka persamaan (III.41) akan berbentuk : n +1

n +1

= k

n

k

S AVE V CORE n S AVE V CORE

= k

n

n +1

S AVE

(III.43)

n S AVE

Atau dalam bentuk persamaan diskrit : S AVE

=

1

N

∑S

N k =1

1

k n +1 = k n

(III.44)

k

nr

nz

∑∑S nr nz i =0 j =0

1

nr

∑∑S

nz

∑∑S i , j

n +1

n i =0 j =0 nr nz

= k

nz

nr nz i =0 j =0

nr

n +1 i , j

∑∑S

n i , j

(III.45) n i , j

i =0 j =0

Untuk kasus difusi multigrup, pertama-tama persamaan (III.7) dipecahkan dengan menggunakan persamaan (III.28), (III.29), dan (III.30). Dengan demikian, kita n +1 akan memiliki nilai dan distribusi fluks untuk setiap grup, yaitu φ g, i, j .

Nilai source dihitung menggunakan persamaan berikut :

Sin, +j 1

G

= ∑ vg 'Σ( fg ' g ) g ,i , jφ gn'+,i1, j

(III.46)

g '=1

Setela Setelah h nilai nilai source diketah diketahui, ui, maka maka selanj selanjutn utnya ya nilai nilai faktor faktor multip multiplik likasi asi k dihitung dihitung dengan dengan menggunak menggunakan an persamaan persamaan (III.45). (III.45). Perhitunga Perhitungan n iterasi iterasi luar ini dilakukan terus-menerus sampai konvergensi numerik tercapai, yaitu :

13


k

n +1

− k n

k

n

< ε 1

dan

1 n S ni ,+ j − S i , j

S ni , j

< ε 2 , untuk seluruh i,j

(III.47)

Dengan demikian, maka kita telah memiliki satu set persamaan lengkap untuk mengerjakan iterasi luar, yaitu menghitung nilai dan distribusi source dan k .

III.2.3 Alur kerja program neutronik Secara garis besar, program neutronik yang penulis kembangkan memiliki alur kerja sebagai berikut : 1. baca baca file file inpu inputt penge pengenda ndali li prog program ram 2. baca baca file file dat dataa cros crosss sect sectio ion n 3. baca aca fil filee pet petaa (core map) komposisi material teras 4. inis inisia iali lisa sasi si nil nilai ai fluk fluks, s, source, dan k 5. defin definis isik ikan an sya syara ratt batas batas 6. kerjak kerjakan an itera iterasi si dalam dalam dan dan iter iterasi asi luar luar 7. tulis hasil-hasil hasil-hasil perhitunga perhitungan n pada pada file-file file-file output output 8. tampilkan tampilkan hasil-hasil hasil-hasil perhitunga perhitungan n pada pada grafik-grafik grafik-grafik

Diagram alir program neutronik akan ditunjukan pada halaman berikutnya.

14


Gambar III.4 Diagram alir program neutronik

15


File input pengendali pengendali program neutronik berada di dalam folder ” input”, dengan nama ”neutronics.txt”, dan berikut ini adalah contohnya :

Tabel III.1 Format file input neutronik Line 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20

neutronics.txt INPUT VARIABLES | VALUES ----------------------------------|-----------GENERAL CONTROL | group |8 core_diameter(cm) | 250 core_height(cm) core_height(cm) | 150 errormax | 1E-12 itermax | 100000 radial_points,max=80 radial_points,max=80 | 20 axial_points(odd),max=81 axial_points(odd),max=81 | 21 flux_guess | 1E+14 source_guess | 1E+14 keff_guess | 1.0 SOR_parameter |0 | CORE DEFINITION | number_of_material_type,max=10 number_of_material_ty pe,max=10 | 3 01 pdsmdl06.c001.Dec04.10.27.33.macs 02 pdsmdl06.c002.Dec04.10.27.52.macs 03 pdsmdl06.c003.Dec04.10.28.06.macs

Keterangan :

•

bar baris is 04 : bany banyak akny nyaa grup grup ener energi gi neut neutro ron n yang yang digu diguna naka kan n dala dalam m perhitungan

•

baris 05 dan 06 : dimensi teras, yaitu diameter dan tingginya

•

baris 07 : tingkat akurasi atau orde error yang akan digunakan dalam perhitungan

•

baris 08 : banyaknya iterasi maksimum maksimum pada tiap loop, untuk menghindari menghindari error stack stack overflow, yaitu bila konvergensi numerik gagal tercapai

•

baris 09 : banyaknya partisi pada arah r (radial)

•

baris 10 : banyaknya partisi pada arah z (aksial)

•

baris 11 : inisialisasi nilai fluks

•

baris 12 : inisialisasi nilai source

•

baris 13 : inisialisasi nilai k

•

baris 14 : parameter akselerasi SOR

•

baris 17 – 20 : banyaknya dan definisi jenis material penyusun teras

16


File peta komposisi material teras berada di dalam folder ”input”, dengan nama ”coremap1.txt”.

Gambar III.5 Sistem pemetaan teras

File coremap1.txt tersebut memberikan informasi jenis material di setiap titk partis partisii teras, teras, baik baik partis partisii radial radial (horis (horisont ontal), al), maupun maupun partis partisii aksial aksial (verti (vertikal kal), ), daerah pemetaan berupa sistem array dua dimensi R-Z, seperti ditunjukan pada gambar III.5 diatas. Berikut ini adalah contohnya :

Tabel III.2 Format file input pemetaan teras Line 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17

coremap1.txt 33333333 33333333 22222222 22222222 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111

3 3 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1

3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

3 3 3 3 3 3 3 3 3 3 3C 3 3 3 3 3 3

17


18 19 20 21

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

2 2 3 3

3 3 3 3

3 3 3 3

Tabel Tabel III.2 III.2 diatas diatas adalah adalah contoh contoh perhit perhitung ungan an dengan dengan menggu menggunak nakan an 20 partis partisii radial dan 21 partisi aksial, dan 3 jenis material. Angka ”1” pada Tabel III.2 menunjukan bahwa material di titik partisi tersebut adalah material jenis ”1” yang data cross section-nya berada di dalam file yang namanya ditunjukan pada baris 18 di tabel III.1. Demikian juga arti angka ”2” dan ”3”. Material di dalam teras biasanya dibedakan berdasarkan level enrichment pada bahan bakarnya.

File File output output hasil hasil perhit perhitung ungan an disimp disimpan an di dalam dalam folder folder ”output”, terdiri terdiri dari dari beberapa file, diantaranya adalah file yang menyimpan data nilai fluks dan source di setiap partisi, dan file output umum yang bernama ” out.neutronics.txt”.

Tabel III.3 Format file output perhitungan neutronik Lin Line 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22

out out.neu .neutr tron onic ics s.tx .txt NEUTRONICS CALCULATION, RESULTS k-eff = 1.01701705819494 Flux Max = Flux[0,0] = 722510203052706 #/cm2.s Flux Mean = 264728516229622 #/cm2.s Flux Peaking Factor = 2.72924962275696 wflux[1] = 0.158018861888309 wflux[2] = 0.226427140655164 wflux[12] = 0.214468792303527 wflux[4] = 0.141663662275637 wflux[5] = 0.0940651397076077 wflux[6] = 0.0752527701029753 wflux[13] = 0.189037055437533 wflux[8] = 1 BG2 = 0.00080882668449286 BG = 0.0284398784190942 Iteration = 36839 Start :Thursday, 28 28 December 2006, 2006, 07:58:42:343 07:58:42:343 Finish :Thursday, 28 December 2006, 07:58:46:750 ---------------------------------------------------------------------------------------------Results are saved to "output/out.neutronics.txt" "output/out.neutronics.txt"

Graf Grafik ik hasi hasill

perh perhit itun unga gan n

seca secara ra

otom otomat atis is lang langsu sung ng dita ditamp mpil ilka kan n

perhitungan selesai dikerjakan. Berikut ini beberapa contohnya :

18

sete setela lah h


Gambar III.6 Contoh grafik distribusi fluks radial 8 grup

Gambar III.7 Contoh grafik distribusi fluks aksial 8 grup

Data cross section untuk setiap jenis material disimpan dalam sebuah file tunggal, yang namanya ditunjukan pada baris 18 – 20 tabel III.1. Pada bagian selanjutnya akan dijelaskan mengenai data cross section tersebut secara lebih rinci.

Keselu Keseluruh ruhan an sistem sistem progra program m komput komputer er yang yang penuli penuliss kemban kembangka gkan n ini dibuat dibuat dengan Borland Delphi 7.0, dan penulis namakan ” Preliminary Nuclear Plant

Analysis Code” atau disingkat PRENPAC, dan bisa di download secara gratis di situs http://syeilendrapramuditya.wordpress.com.

19

Solusi Numerik Persamaan Difusi Neutron

Recommend Documents