Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
PENGGUNAAN METODE PN DAN SN UNTUK PERPINDAHAN KALOR RADIASI PADA REAKTOR NUKLIR TERFLUIDISASI (Studi Kasus untuk Pengambilan Kalor Peluruhan Secara Pasif) Alexander Agung Jurusan Teknik Fisika, Universitas Gadjah Mada Jl. Grafika 2, Yogyakarta 55281 Corresponding author: Tel. 0274 902120, Fax. 0274 902210, email:
[email protected]
Danny Lathouwers, Tim van der Hagen, Hugo van Dam Department of Radiation, Radionuclides and Reactors, Delft University of Technology, Mekelweg 15, 2629 JB Delft, The Netherlands ABSTRAK Telah dilakukan simulasi numerik untuk mempelajari kemungkinan pengambilan kalor peluruhan secara pasif pada reaktor nuklir terfluidisasi. Perpindahan kalor secara konduksi, konveksi dan radiasi diperhitungkan dan model turbulen k-ε bilangan Reynolds tinggi digunakan untuk perhitungan aliran fluida. Metode pendekatan P1 dan SN untuk media berpartisipasi digunakan untuk model radiasi. Reaktor dimodelkan sebagai silinder aksisimetris dua dimensi. Perhitungan dilakukan untuk daya operasi awal sebesar 60 MW. Hasil perhitungan menunjukkan bahwa metode P1 cukup baik untuk diterapkan pada bed partikel tetapi gagal pada bagian freeboard. Metode SN orde tinggi dapat digunakan baik pada bagian bed partikel maupun freeboard. Kekurangan dari metode SN orde tinggi ini adalah lamanya waktu komputasi. Ketidakvalidan metode P1 menyebabkan hasil perhitungan suhu bahan bakar maksimum melebihi batas yang diizinkan. Hal yang demikian tidak dijumpai pada metode SN. Kata Kunci: kalor peluruhan, reaktor nuklir terfluidisasi, pengambilan kalor secara pasif, simulasi numerik, metode PN, metode SN. ABSTRACT Numerical simulations have been performed to investigate the possibility of passive decay heat removal in a fluidized bed nuclear reactor. The conduction, convection and radiative heat transfer were included and the high Reynolds number k-ε turbulence model was applied for the flow calculations. The P1 and SN approximation methods for participating media were used for the radiation model. The reactor was modeled as a 2D axisymmetric cylinder. Calculations were performed for an initial operating power of 60 MW. The results show that P1 method is quite good to be applied in the particle bed but it fails in the freeboard. The high order SN method can be applied extremely good both in the particle bed and the freeboard. The drawback of the SN method is its long computational time. The invalidity of P1 method leads to a high fuel particle temperature, exceeding its allowable maximum value. Such situation, however, does not occur in the SN method. Keywords: decay heat, fluidized bed nuclear reactor, passive heat removal, numerical simulation, PN methods, SN methods.
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
705
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
1. PENDAHULUAN FLUBER (Fluidized Bed Thermal Fission Nuclear Reactor) adalah desain reaktor inovatif yang menggunakan konsep fluidisasi gas-padatan. Reaktor ini tersusun atas partikel bahan bakar TRISO yang diisikan ke dalam silinder berdinding grafit. Helium digunakan sebagai pendingin dan sebagai medium fluidisasi. Ketika tidak ada aliran helium ke dalam teras atau ketika besarnya aliran ke atas masih kecil, partikel bahan bakar terkumpul di bagian bawah teras dan reaktor menjadi sangat subkritis dikarenakan kurangnya moderasi neutron. Ketika aliran bertambah, partikel bahan bakar menjadi terfluidisasi, teras menjadi mengembang dan reaktivitas meningkat dikarenakan bertambahnya proses moderasi pada reflektor grafit.
Spesifikasi umum desain
FLUBER ditampilkan pada Tabel 1 dan tampilan skematis dari reaktor pada keadaan packed dan terfluidisasi dapat dilihat pada Gambar 1. Ketika aliran helium ke dalam teras berkurang baik disengaja maupun akibat kegagalan pompa, susunan partikel bahan bakar menjadi packed bed dan reaktor kembali menjadi subkritis. Akan tetapi kalor peluruhan masih dihasilkan di teras dan harus dipindahkan. Pada skenario kondisi terburuk di mana tidak ada pompa yang bekerja, kalor peluruhan harus diambil secara pasif. Karena daya output ketika reaktor masih beroperasi relatif tinggi, pengambilan kalor peluruhan setelah reaktor shutdown menjadi sangat penting. Ketika reaktor beroperasi pada ekspansi bed yang maksimal, daya total maksimum dapat mencapai 120 MW dengan suhu bahan bakar maksimum sebesar 1165 K. Setelah reaktor mati (dengan cara tidak mengalirkan helium ke dalam reaktor), packed bed terbentuk dengan volume sebesar 2,5 m3. Daya peluruhan pada awal transien sebesar 7% dari daya total sehingga densitas daya peluruhan saat itu mencapai sekitar 3,4 MW/m3. Tujuan utama dari penelitian ini adalah untuk mengetahui apakah kalor peluruhan yang dihasilkan dari pengoperasian FLUBER dapat diambil secara pasif tanpa
melampaui
batas
keselamatan.
Untuk
CP
(coated
particles),
IAEA
merekomendasikan nilai sebesar 1873 K sebagai batas konservatif suhu bahan bakar maksimum pada kondisi kecelakaan [1]. Oleh karena itu diperlukan perhitungan yang bersifat best estimate dengan menyertakan model untuk semua mode perpindahan kalor (konduksi, konveksi dan radiasi). Pada paper sebelumnya [2] telah dilaporkan hasil perhitungan yang menggunakan metode PN untuk model perpindahan radiasi termal. Paper ini merupakan pengembangan lebih lanjut dari kegiatan yang telah dilaporkan
706
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
pada paper tersebut yaitu dengan menggunakan metode SN untuk model perpindahan radiasi termal dan membandingkan hasilnya dengan metode PN. Tabel 1. Spesifikasi FLUBER Jari-jari teras [cm]
79,8
Tinggi teras [cm]
600
Tinggi seluruh reaktor [cm]
800
Tebal reflektor aksial dan radial [cm]
100
Tebal dan tinggi penyerap radial [cm]
50
Konsentrasi penyerap radial [ppm]
50
Tinggi packed bed (pada 40% of porositas) [cm] Massa uranium [kg]
122,36 220
Pengkayaan [% berat]
16,76
Tekanan helium [bar]
60
Daya maksimum (ketika ekspansi penuh) [MW]
120
Suhu bahan bakar maksimum [K]
1165
Suhu gas maksimum [K]
1160
Reflektor grafit Helium cavity
Partikel bahan bakar Penyerap samping Gambar 1. Skema FLUBER pada keadaan packed (kiri) dan terfluidisasi (kanan). Susunan dari paper ini adalah sebagai berikut. Pada Bagian 2 akan diuraikan tentang model matematis yang menggambarkan fenomena yang terjadi pada FLUBER setelah shutdown. Pada bagian ini pula metode PN dan SN diuraikan dengan ringkas. Selanjutnya implementasi numerik terhadap model matematis pada FLUBER dijelaskan
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
707
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
pada Bagian 3 dan diikuti dengan hasil perhitungan yang dibahas pada Bagian 4. Selanjutnya Bagian 5 mengakhiri paper ini dengan kesimpulan. 2. MODEL MATEMATIS Pengambilan kalor peluruhan pada FLUBER dapat dideskripsikan dengan gabungan perpindahan kalor konduksi, konveksi dan radiasi pada medium berpori. Skenario kecelakaan yang diperhitungkan adalah bahwa sistem tetap pada kondisi operasi 60 bar dan bahwa pada sistem tidak terdapat konveksi paksa. 2.1 Persamaan kontinuitas Persamaan kontinuitas untuk fase gas dapat dituliskan sebagai
∇ ⋅ (α g ρ g u ) = 0
(1)
dengan α adalah fraksi volume, ρ adalah densitas, u adalah vektor kecepatan dan indeks g menyatakan fase gas.
2.2 Persamaan momentum Persamaan momentum untuk fase gas dapat dinyatakan sebagai [3,4] ∂ (α g ρ g u ) ∂t
+ ∇ ⋅ (α g ρ g uu ) = −
α pρp u + ∇ ⋅ α g τ g − α g ∇P + α g ρ g g τ 12
(2)
dengan τ 12 adalah skala waktu interaksi gas-partikel yang dihitung berdasarkan relasi Ergun, τ adalah tensor stress, P adalah tekanan, g adalah percepatan gravitasi dan indeks p menyatakan fase padatan. Fase gas diasumsikan memenuhi hukum gas ideal. Untuk menyelesaikan persamaan aliran, digunakan pula model turbulen k-ε untuk bilangan Reynolds tinggi. 2.3 Persamaan energi
Persamaan energi untuk fase gas dan padatan adalah
α g ρ g C p,g
α p ρ p C p. p
∂T p ∂t
D g Tg Dt
]
[
[
]
[
]
]
= ∇ ⋅ λ p ∇T p − h pg a pg T p − Tg + Pd''' − ∇ ⋅ q rad
ρ refl C p ,refl
708
[
= ∇ ⋅ λ g ∇Tg + h pg a pg T p − Tg
∂Trefl ∂t
[
= ∇ ⋅ λ refl ∇Trefl
]
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
(3) (4) (5)
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
dengan Cp adalah kapasitas kalor, T adalah suhu, λ adalah konduktivitas termal, hpg adalah koefisien perpindahan kalor konveksi antara fase gas dan padatan yang diperoleh dari korelasi Ranz, apg adalah luasan antarmuka per satuan volume partikel padatan, Pd''' adalah sumber volumetrik kalor peluruhan, qrad adalah fluks kalor radiasi dan indeks refl menyatakan bagian reflektor. D/Dt pada Pers. 3 adalah turunan substansial yang terhubung dengan turunan parsial terhadap waktu dan koordinat ruang. Sifat-sifat fase gas dan padatan diperoleh dari data-data yang tersedia [5 – 7]. 2.4 Kalor peluruhan
Reaktor diasumsikan telah beroperasi cukup lama sebelum shutdown sehingga daya peluruhan telah mencapai keadaan tunak. Kalor peluruhan, Pd, setelah shutdown dapat dinyatakan dalam bentuk sebagai berikut Pd =
Ptot
Nd
γ Qf + ∑ n n =1 λ n Nd
γn
∑λ n =1
e − λn t
(6)
n
dengan Ptot adalah daya total yang dibangkitkan sebelum shutdown, Qf adalah energi serentak yang dapat diambil per fisi, γ n dan λn adalah yield dan konstanta kalor peluruhan dan Nd adalah banyaknya kelompok kalor peluruhan. Pada penelitian ini digunakan 23 kelompok kalor peluruhan dengan data yang diperoleh dari Ref. [8]. Sumber kalor volumetrik pada Pers. 6 diberikan sebagai Pd''' = Pd / Vcore
(7)
dengan Vcore adalah volume teras aktif. Sumber kalor peluruhan dibangkitkan secara homogen pada packed bed dikarenakan burn-up bahan bakar yang homogen selama proses fluidisasi. Hal ini sebagai akibat proses pengadukan bahan bakar yang sangat bagus yang merupakan salah satu karakteristik dari fluidisasi.
2.5 Model radiasi Perpindahan kalor secara radiasi dimodelkan menggunakan persamaan integrodiferensial untuk media berpartisipasi (dapat menyerap dan menghamburkan sinar) [9] sebagai berikut
Ω ⋅ ∇I + Σ t I = Σ a I b +
Σs 4π
∫ I (Ω')Φ(Ω'⋅Ω)dΩ'
(8)
4π
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
709
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
dengan I adalah intensitas radiasi, Ib adalah intensitas benda hitam ( = σT 4 / π ), Σ t adalah koefisien extinction, Σ s adalah koefisien hamburan, Σ a adalah koefisien serapan, σ adalah konstanta Stefan-Boltzman, Ω adalah vektor arah dan Φ adalah fungsi fase hamburan yang menggambarkan probabilitas sinar dengan arah Ω' akan dihamburkan ke arah Ω tertentu. Pers. 8 di atas biasa disebut sebagai Radiative Transfer Equation (RTE). Fungsi fase hamburan dapat dinyatakan dalam N
Φ(Ω'⋅Ω ) = ∑ (2i + 1)Ai Pi (Ω'⋅Ω )
(9)
i =0
Pada perhitungan ini diasumsikan hamburan yang isotropik. Sifat radiasi dari partikel diperoleh dari teori tak gayut (independent) yang artinya bahwa karakteristik serapan dan hamburan sebuah partikel tidak dipengaruhi oleh partikel-partikel tetangganya. Untuk partikel yang besarnya melebihi panjang gelombang termal berlaku hubungan sebagai berikut
Σa = Σs = Σt = 3
3αp 2 dp
(10)
αp
(11) dp dengan dp adalah diameter partikel. Helium yang terdapat pada freeboard bersifat transparan terhadap radiasi termal sehingga tidak ikut dalam proses radiasi. 2.5.1 Metode PN
Metode PN atau disebut juga pendekatan spherical harmonics mula-mula dikembangkan untuk studi perpindahan radiasi pada atmosfer, kemudian dimodifikasi untuk penyelesaian permasalahan transport neutron dan kemudian digunakan untuk permasalahan perpindahan radiasi termal [10]. Pada metode ini, intensitas radiasi dinyatakan sebagai urutan (series) spherical harmonics dan dapat dituliskan sebagai N
l
I (r, Ω ) = ∑ ∑ Alm (r )Yl m (Ω )
(12)
l =0 m = −l
dengan spherical harmonics ditulis sebagai
dan Pl
710
m
( ) (l + m )!
(m + m ) / 2 l − m !
Yl (Ω ) = (− 1) m
1/ 2
e jmφ Pl
m
(cos θ )
adalah polinomial Legendre terasosiasi.
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
(13)
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
Pada Pers. 12 batas atas N untuk indeks l dikenal dengan istilah orde pendekatan. Penyelesaian eksak dari persamaan transfer radiasi dapat dicapai jika N adalah tak hingga, akan tetapi untuk keperluan perhitungan praktis digunakan nilai N yang berhingga. N=1 menghasilkan pendekatan P1 dan N=3 menghasilkan pendekatan
P3. Biasanya digunakan orde spherical harmonics yang ganjil karena untuk menghindari singularitas matematis dari intensitas radiasi. Untuk pendekatan P1, setelah mengalikan intensitas dengan momen ke-nol dan pertama dan mengintegralkannya terhadap sudut ruang, diperoleh
I (r , Ω ) =
1 [g (r ) + 3q(r ) ⋅ Ω] 4π
(14)
dengan G adalah incident radiation (= ∫ I (Ω )dΩ ) dan q adalah fluks kalor radiasi. 4π
Selanjutnya dengan mensubstitusi nilai intensitas ke RTE (Pers. 8) dan kemudian mengalikannya dengan momen ke-nol dan pertama, akan dihasilkan persamaan difusi
∇ ⋅ D∇G (r ) + Σ a G (r ) = 4Σ aσT 4
(15)
dengan D adalah koefisien difusi (= − (3Σ t − A1Σ s ) ). −1
Untuk metode P1 ini dapat digunakan syarat batas Marshak sebagai berikut:
∫ I (r
w
, Ω )nˆ ⋅ ΩdΩ =
nˆ ⋅Ω > 0
∫ I (Ω)nˆ ⋅ ΩdΩ
(16)
w
nˆ ⋅Ω > 0
Persamaan transfer radiatif (Pers. 8) menjadi terkopel dengan persamaan energi (Pers. 4) melalui sumber radiasi,
[
∇ ⋅ q rad = Σ a 4σT p4 − G
]
(17)
2.5.2 Metode SN Pendekatan
SN atau
pendekatan
discrete
ordinates
diperoleh
dengan
mendiskretisasi seluruh sudut ruang ( Ω = 4π ) menggunakan sejumlah arah ordinat yang terbatas dan faktor pembobot yang sesuai. RTE kemudian dituliskan untuk setiap ordinat dan suku integral digantikan dengan jumlahan kuadratur masing-masing ordinat. Pada metode SN, terdapat nilai diskret N dari kosinus arah ξ n ,η n , µ n yang selalu memenuhi identitas ξ n2 + η n2 + µ n2 = 1 [11,12]. Orde N pada metode SN berupa bilangan genap dan banyaknya ordinat untuk kasus dua dimensi dapat dinyatakan dengan N ( N + 2 ) / 2 . Dengan demikian untuk S2 terdapat 4 ordinat, S4 terdapat 12 ordinat, S6
terdapat 24 ordinat dan S8 terdapat 40 ordinat.
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
711
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
RTE kemudian diselesaikan untuk setiap ordinat dan nilai incident radiation dengan mudah dapat dituliskan sebagai N
G (r ) ≈ ∑ wm I m (r )
(18)
m =1
dengan wm adalah nilai bobot dari titik kuadratur. Berdasarkan pada ekspresi umum untuk permukaan tak tembus cahaya yang memancarkan dan menghamburkan secara menyebar, syarat batas untuk metode SN dapat dituliskan sebagai berikut
I m (rw ) = ε (rw )I b (rw ) +
ρ (rw ) ∑ w j I j (rw ) nˆ ⋅ Ω , nˆ ⋅ Ω m > 0 π nˆ ⋅Ω <0
(19)
j
3. IMPLEMENTASI PADA FLUBER 3.1. Metodologi penyelesaian Persamaan diferensial diselesaikan menggunakan metode volume hingga (finite volume) dengan staggered grid untuk silinder aksisimetri dua dimensi. Semua fluks konvektif didekati dengan skema Total Variation Diminishing (TVD) orde dua. Diskretisasi waktu berdasarkan skema Backward Euler yang dikombinasikan dengan teknik koreksi tekanan. Semua sistem linear yang muncul dari diskretisasi sistem konveksi-difusi diselesaikan dengan metode Conjugate Gradient (CG) untuk persamaan tekanan dan Generalized Minimum RESidual (GMRES) untuk persamaan transport yang lain. Pengkondisian awal dari solver linear ini didasarkan pada faktorisasi ILU. 3.2. Diskretisasi ruang dan sudut ruang Geometri reaktor didiskretisasi dengan koordinat silinder aksisimetri 2D dengan mesh yang tidak homogen. Teras didiskretisasi menjadi 20 (arah radial) kali 120 (arah aksial), sementara tinggi dari packed bed didiskretisasi menjadi 26 volume kontrol. Reflektor dibagi menjadi 20 volume kontrol tambahan untuk arah radial dan aksial. Untuk metode SN, sudut ruang didiskretisasi berdasarkan orde dari SN. Dengan demikian pada masing-masing titik mesh terdapat N(N+2)/2 ordinat yang bersesuaian.
712
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
3.3. Time-stepping
Semua perhitungan dilakukan menggunakan time step 0,25 detik dengan kemampuan menyimpan hasil sementara dan restart perhitungan setiap 900 detik real time. Studi awal mengenai sensitivitas time step menunjukkan tidak adanya instabilitas hasil yang signifikan bila dibandingkan dengan nilai time step yang lebih ketat, terutama untuk perhitungan aliran dan model turbulen. 3.4. Syarat awal dan syarat batas Suhu awal partikel dan helium adalah uniform di seluruh bed dan freeboard sesuai dengan pencampuran yang bagus seperti telah disebutkan di atas. Nilai awal temperatur tersebut diperoleh dari perhitungan yang dilakukan pada langkah sebelumnya menggunakan code dinamika titik. Pada daya 60 MW digunakan suhu awal bahan bakar sebesar 922 K dan suhu awal gas sebesar 920 K. Nilai awal suhu reflektor diperoleh dari perhitungan 2D kondisi tunak dengan suhu batas dinding luar sebesar 298 K dan batas dinding dalam dipasang sebesar 922 K. Kecepatan awal gas diset nol. Syarat batas suhu pada dinding luar samping dan atas diperoleh dari konveksi alami dan radiasi benda abu-abu ke atmosfer. Suhu dinding reflektor bawah diasumsikan konstan pada 500 K. Suhu atmosfer dalam hal ini adalah 298 K. Pada batas internal antara reflektor dan teras, syarat batas diperoleh dari ketiga mode perpindahan kalor. Berkaitan dengan fluks radiasi pada reflektor, dinding dianggap
sebagai
permukaan
tak
tembus
cahaya
yang
memancarkan
dan
menghamburkan secara menyebar. Berkaitan dengan mode konveksi-konduksi, syarat batas diperoleh dari fungsi dinding (wall function) untuk kecepatan dan suhu gas. Fraksi volume bahan bakar pada bed adalah konstan sebesar 0,6.
Pada
freeboard nilai fraksi volume ini bervariasi dan akan ditampilkan pada Bagian 4. 3.5. Code Penyelesaian dari persamaan-persamaan diferensial di atas diimplementasikan pada code DECFLUB. Simulasi dilakukan pada AMD64 2800+ untuk waktu 3 jam real time. CPU time untuk perhitungan semacam ini berbeda-beda untuk masing-masing model radiasi yang digunakan dan akan dibahas pada Bagian 4.
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
713
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
4. HASIL DAN PEMBAHASAN Gambar 2 menunjukkan besarnya kalor yang dipindahkan ke permukaan reflektor sebelah bawah di mana reflektor bersinggungan secara langsung dengan partikel bahan bakar dengan fraksi volume sebesar 0,6. Dengan demikian partikel bahan bakar bersifat participating terhadap radiasi termal. Hasil perhitungan menunjukkan perbedaan yang besarnya kurang dari 1% sehingga hal ini berarti bahwa tidak ada perbedaan yang signifikan antara kedua metode. Fraksi volume partikel = 0,6 -134.00 -133.75
kW
-133.50 -133.25 -133.00 -132.75 -132.50 P1
S2
S4
S6
S8
Gambar 2. Kalor yang dipindahkan pada permukaan reflektor bagian bawah. Perhitungan menggunakan metode P1 dan SN dengan fraksi volume partikel pada bed sebesar 0,6. Besarnya kalor yang dipindahkan ke permukaan reflektor sebelah atas ditampilkan pada Gambar 3. Gambar sebelah kiri menyatakan kondisi ketika fraksi volume partikel (αp) pada freeboard sebesar 10-3 dan gambar sebelah kanan ketika αp sebesar 10-12. Pada masing-masing gambar ditampilkan hasil untuk perhitungan menggunakan metode P1, S2, S4, S6 dan S8. Sebagai pembanding digunakan metode View Factor (VF) yang lebih akurat untuk digunakan pada medium transparan. Gambar 3 (kiri) dengan jelas menunjukkan bahwa semua metode pendekatan yang digunakan (P1 dan SN orde rendah maupun tinggi) memberikan
perbedaan yang signikan
dibandingkan dengan hasil perhitungan VF. Hal ini menunjukkan bahwa αp sebesar 103
belum cukup transparan sehingga medium pada freeboard menjadi bersifat
participating.
714
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005 Fraksi volume partikel = 1e-3
Fraksi volume partikel = 1e-12
-60
-20
-50
0
-40
20
VF
kW
kW
-30
S2
S4
S6
S8
40
-20
60
-10
80
0
P1
100
VF
P1
S2
S4
S6
S8
Gambar 3. Kalor yang dipindahkan pada permukaan reflektor bagian atas. Perhitungan menggunakan metode P1 , SN dan VF (sebagai pembanding) dengan fraksi volume partikel pada freeboard sebesar 10-3 (kiri) dan 10-12 (kanan).
Gambar 4. Distribusi incident radiation pada teras dengan besar fraksi volume partikel sebesar 0,6 pada bed dan 10-12 pada freeboard. Hasil diperoleh dari perhitungan menggunakan metode P1 (kiri), S2 (tengah) dan S8 (kanan). Gambar 3 (kanan) dengan jelas menunjukkan bahwa arah aliran perpindahan kalor pada permukaan untuk metode P1 adalah tidak benar. Pada kondisi ini, kalor mengalir dari permukaan reflektor yang bersuhu lebih rendah menuju ke permukaan tumpukan partikel bahan bakar yang bersuhu lebih tinggi. Hal ini menunjukkan bahwa metode P1 menjadi tidak valid untuk αp yang sangat kecil. Nilai koefisien extinction Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
715
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
menjadi sangat kecil yang selanjutnya membuat nilai koefisien difusi menjadi sangat besar. Metode SN memberikan hasil yang berbeda dengan metode P1. Arah perpindahan kalor untuk metode ini sudah benar dan besarnya dalam rentang yang sesuai dengan metode VF. Semakin tinggi orde SN, semakin kecil selisih hasil dengan metode VF. Metode S2 memberikan hasil yang overpredict dikarenakan ordinat yang terlalu sedikit. Analogi dengan transport neutron, hal ini menyebabkan adanya efek pancaran (ray effects). Gambar 4 menunjukkan distribusi incident radiation pada teras dengan fraksi volume partikel pada freeboard sebesar 10-12. Seperti telah disebutkan di atas, metode P1 menjadi tidak valid untuk αp yang kecil sehingga incident radiation pada freeboard bernilai sangat kecil. Hal semacam ini tidak dijumpai pada metode SN dan hal ini menunjukkan pula bahwa metode ini dapat digunakan meskipun pada αp yang sangat kecil. Perhatikan pula adanya gradasi incident radiation pada freeboard yang berdekatan dengan permukaan atas bed. Suhu maksimum bahan bakar sebagai fungsi waktu untuk berbagai metode perhitungan perpindahan kalor radiasi ditampilkan pada Gambar 5.
Secara umum
perilaku transien ini menunjukkan pola yang hampir sama. Transien termal pada teras aktif menunjukkan adanya peningkatan suhu bahan bakar dalam waktu kurang dari dua jam setelah shutdown sampai mencapai suhu puncak yang kemudian diikuti oleh pendinginan secara perlahan. Teras pada awalnya menjadi panas karena pembangkitan kalor peluruhan melebihi laju pengambilan kalor. Kemudian laju pengambilan kalor pada teras aktif melebihi pembangkitan kalor peluruhan dan semenjak itu suhu teras rata-rata menjadi turun karena kalor ditransfer ke reflektor dan atmosfer. Tampak pada Gambar 5 bahwa penggunaan metode P1 menyebabkan suhu bahan bakar maksimum menjadi sangat tinggi bahkan melampaui nilai batas yang diizinkan (1873 K).
Nilai fraksi volume partikel yang digunakan pada freeboard
sebesar 10-12 sehingga arah perpindahan kalor pada permukaan reflektor yang bersinggungan dengan freeboard
menjadi terbalik. Hal ini berarti bahwa alih-alih
perpindahan kalor dari bahan bakar ke reflektor, yang terjadi adalah pemanasan bahan bakar akibat radiasi termal dari reflektor.
716
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005 2000
Suhu bahan bakar maksimum yang diizinkan
1900 1800
S4, S6, S8
1700 P1
1600 1500
S2
1400 1300 1200 1100 1000 900 0
0.5
1
1.5
2
2.5
3
Waktu (jam)
Gambar 5. Suhu bahan bakar maksimum sebagai fungsi waktu pada daya operasi reaktor sebelum shutdown sebesar 60 MW. Simulasi dilakukan dengan menggunakan semua mode perpindahan kalor. Metode
S2
memberikan
perpindahan
kalor
radiasi
yang
overpredict
dibandingkan dengan metode SN berorde lebih tinggi. Dengan demikian partikel bahan bakar menjadi lebih cepat dingin. Hal ini dapat dilihat pada Gambar 5 bahwa suhu bahan bakar maksimum untuk metode S2 lebih rendah daripada metode SN orde lebih tinggi. Meskipun demikian perbedaan hasil antara metode S2 dengan metode SN orde lebih tinggi tidak lebih dari 30 K. Metode S4, S6 dan S8 memberikan hasil yang tidak berbeda secara signifikan (kurang dari 1 K). Akan tetapi seperti tampak pada Gambar 3 (kanan) hasil perhitungan laju kalor menggunakan metode S4 masih berbeda secara signifikan terhadap hasil perhitungan menggunakan metode VF. Untuk metode S8, hasil yang diperoleh sudah cukup mendekati sehingga penggunaan orde tinggi ini dapat dianjurkan. Meskipun demikian terdapat kendala yaitu mahalnya perhitungan menggunakan orde tinggi. Tabel 2 menunjukkan perbandingan CPU time untuk perhitungan real time sebesar 3 jam dengan menggunakan berbagai metode. Metode S8 ternyata 5,6 kali lebih lambat daripada metode S4. Lamanya waktu perhitungan pada orde tinggi ini disebabkan adanya proses penyisiran (sweeping) untuk setiap ordinat pada masingmasing titik mesh hingga tercapai konvergensi nilai intensitas radiasi. Dengan time disadvantage ratio (yaitu perbandingan antara CPU time dengan real time) yang tinggi pada metode S8 (yaitu sekitar 8,9), perhitungan transien menggunakan metode ini untuk Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
717
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
mensimulasikan kondisi reaktor lama setelah shutdown ( > 1 hari) menjadi sangat tidak ekonomis. Tabel 2. Perbandingan CPU time untuk simulasi real time sebesar 3 jam. Metode CPU Time (jam)
P1
S2
S4
S6
S8
1,53
2,29
4,78
11,56
26,69
5. KESIMPULAN DAN PENELITIAN LANJUTAN Simulasi numerik untuk pengambilan kalor secara pasif pada FLUBER telah dilakukan dengan menggunakan beberapa metode untuk penyelesaian perpindahan kalor radiasi. Dari hasil perhitungan didapatkan kesimpulan sebagai berikut: 1. Metode P1 dapat digunakan pada medium berpartisipasi dengan ketebalan optis (yang direpresentasikan oleh nilai fraksi volume partikel) yang besar tanpa ada perbedaan yang signifikan dengan metode SN. Hal ini tampak pada bagian bed partikel. Akan tetapi pada bagian freeboard di mana ketebalan optisnya sangat kecil metode P1 menjadi tidak valid. 2. Metode S2 memberikan hasil yang overpredict dikarenakan sedikitnya jumlah ordinat yang digunakan. Pada metode SN orde tinggi tidak dijumpai adanya permasalahan bahkan untuk ketebalan optis yang sangat kecil. Semakin tinggi orde SN, semakin akurat hasil perhitungannya. 3. Sebagai konsekuensi dari kegagalan metode P1 pada bagian freeboard, pada daya operasional sebesar 60 MW, suhu bahan bakar maksimum partikel bahan bakar pasca shutdown menjadi terlalu tinggi bahkan melampaui batas yang diizinkan. Hal yang demikian tidak dijumpai pada metode SN. Permasalahan utama yang dijumpai apabila menggunakan metode SN orde tinggi adalah lamanya waktu komputasi. Untuk mengatasi hal ini perlu dilakukan penelitian lanjutan yang bertujuan untuk mempercepat waktu komputasi misalnya dengan menambah algoritma akselerasi maupun dengan menggunakan metode gabungan Full Approximation Scheme dan Newton-Krylov [13].
718
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
DAFTAR PUSTAKA 1. IAEA, “Fuel Performance and Fission Product Behaviour in Gas-Cooled Reactors”, IAEA-TECDOC-978, Vienna, 1997. 2. A. AGUNG, D. LATHOUWERS, T.H.J.J VAN DER HAGEN, H. VAN DAM, C.C. PAIN, C.R.E. DE OLIVEIRA, A.J.H. GODDARD, M.D. EATON, J.L.M.A. GOMES, B. MILES, “Passive Decay Heat Removal in a Fluidized Bed Nuclear Reactor”, Proc. PHYSOR 2004 (on CD-ROM), Chicago, Illinois, 25-29 April 2004. 3. M. KAVIANY, Principles of Heat Transfer in Porous Media, Springer, New York, 1991. 4. D. LATHOUWERS, J. BELLAN, “Modeling of Dense Gas-Solid Reactive Mixtures Applied to Biomass Pyrolysis in a Fluidized Bed”, International Journal of Multiphase Flow, 27, 2155, 2001. 5. E.W. LEMMON, M.O. MCLINDEN, D.G. FRIEND, “Thermophysical Properties of Fluid Systems” in NIST Chemistry WebBook, NIST Standard Reference Database Number 69, Eds. P.J. Linstrom and W.G. Mallard, National Institute of Standards and Technology, Gaithersburg MD, 20899, 2003. 6. J.K. FINK, “Thermophysical Properties of Uranium Dioxide”,
Journal of
Nuclear Material, 279, 1, 2000. 7. A.T. DINSDALE, “SGTE Data for Pure Elements”, CALPHAD, 15, 317, 1991. 8. Berechnung
der
Nachzerfallsleitung
der
Kernbrennstoffe
von
Hochtemperaturreaktoren mit Kügelförmigen Brennelementen, Technical Report DIN 25 485, 1990. 9. R. SIEGEL, J.R. HOWELL, Thermal Radiation Heat Transfer, Taylor and Francis, Washington, 2001. 10. R. VISKANTA, M.P. MENGÜÇ, “Radiation Heat Transfer in Combustion Systems”, Progress in Energy Combustion Science, 13, 97, 1987. 11. E.E. LEWIS, W.F. MILLER, Computational Methods of Neutron Transport, John Wiley & Sons, New York, 1984. 12. N.Z. CHO, Recent Developments in Reactor Physics and Kinetics, Lecture Notes, The 2003 Frédéric Joliot/Otto Hahn Summer School, Forschungszentrum Karlsruhe, Germany, August 20 – 29, 2003.
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam
719
Prosiding Seminar Nasional ke-11 Teknologi dan Keselamatan PLTN Serta Fasilitas Nuklir ISSN : 0854 - 2910 Malang, 15 September 2005
13. D. BALSHARA, “Fast and Accurate Discrete Ordinates Methods for Multidimensional Radiative Transfer. Part I, Basic Methods”, Journal of Quantitative Spectroscopy & Radiative Transfer, 69, 671, 2001.
720
Alexander Agung, Danny Lathouwers, Tim van der Hagen, Hugo van Dam