PRO SIDING SEMINAR PENELITIAN DAN PENGELOLAAN PERANGKAT NUKLIR Pusat Teknologi Akselerator don Proses Bahan Yogyakarta, 28 Agustus 2008
ANALISIS DISTRIBUSI TEMPERA TUR PADA BATANG PANAS BAGIAN un "HEATING-Ol" Ainur Rosidi, M. Juarsa Pusat Teknologi Reaktor dan Keselamatan Nuk/ir. BATAN Mohamad Riyadi Program Studi Matematika Faku/tas Sains dan Teknologi Universitas Islam Negeri Syarif Hidayatullah Jakarta
ABSTRAK ANALISIS DISTRIBUSI TEMPERATUR UNTUK BATANG PANAS BAGIAN UJI "HEA TiNG-01". Analisis distribusi temperatur untuk batang panas berbentuk silinder pada bagian uji Heating-01 telah dilakukan. Analisis dilakukan untuk mengetahui distribusi temperatur serta membandingkan hasil eksperimen dan hasil perhitungan numerik. Eksperimen dilakukan dengan memanaskan batang panas hingga mencapai temperatur awal 85(fC, kemudian didinginkan secara radiasi hingga mencapai temperatur 9(fC. Sedangkan tiga metode numerik digunakan untuk memperkirakan penurunan temperatur pada batang panas selama pendinginan radiasi. Hasil simulasi diperoleh distribusi temperatur (berupa tabel dan grafik) menggunakan metode FTCS, BTCS, dan Dufort Frankel dengan waktu proses pendinginan mencapai temperatur Iingkungan secara radiasi yang sangat cepat, yaitu kurang dari 22 detik (maksimal 21,6 detik). Dari grafik menunjukkan bahwa pola penurunan temperatur terjadi secaf'3 eksponensial, baik untuk hasil eksperimen dan perhitungan numerik. Kata kunci : distribusi temperatur, anulus, Metode Finite Difference
ABSTRACT TEMPERA TURE DISTRIBUTION ANAL YSIS FOR HEA TED ROD IN "HEA TlNG-01" TEST SECTION. The analysis of temperature distribution for heated rod in cylindrical from in heating-01 test section has been done. Analyzing done to know distribution of temperature and compares result of experiment and result of calculation numerik. Eksperiment was conducted with heat-up the heated rod until initial temperature 85(fC and cooling down by radiation until temperature 9(fC. While, three methods of numeric was used to predicted a temperature decrease in heated rod during cooling radiation. Result of simulation is obtained by distribution of temperature (in the form of tables and graph) applies method FTCS, BTCS, and Dufort Frankel with time processed cooling to reach environmental temperature in a real radiation quickly, that is less than 22 seconds (maximum 21,6 seconds). From graph of comparison shows that the pattern of temperature decreasing was occurred in exponential/ly, both for experiment result and numerical result. Keyword: heat distribution, annulus, Finite Difference Method
PENDAHULUAN
Panas energididefinisikan antara dua sebagai sistem bentuk (atau
perpindahan sistem dan lingkungannya) berdasarkan perbedaan temperaturnya. Panas merupakan energi dalam keadaan transisi, yang dikenal hanya ketika panas melintasi batas sebuah sistem. Dalam
Ainur Rosidi, dkk.
termodinamika, keadaan ini dikenal sebagai perpindahan panas. Perpindahan panas dapat dibagi menjadi tiga cara perpindahannya, yaitu konduksi, konveksi, dan radiasi. Semua cara perpindahan panas membutuhkan keberadaan panas yang berbeda dan terjadinya perbedaan dari temperatur tinggi ke temperatur rendah (pendinginan), atau sebaliknya
ISSN 1410 - 8178
119
PENELITIAN
PRO SIDING SEMINAR DAN PENGELOLAAN PERANGKAT
Pusa' Teknologi
Akselera'or
Yogyakarta,
28 Agustus 2008
yaitu terjadinya perubahan dari temperatur rendah ke temperatur tinggi (pemanasan).(I] Berdasarkan analisisnya, perpindahan panas dapat dibedakan menjadi perpindahan panas satu dimensi, dua dimensi dan tiga dimensi. Proses perpindahan panas jika dikaitkan dengan waktu [I] . dapat dikategorikan sebagai perpindahan sebagai perpindahan panas transient (perubahan besaran berdasarkan perubahan waktu). Sedangkan perpindahan panas tunak (steady) adalah proses perpindahan yang tidak tergantung kepada waktu. Salah satu perpindahan panas transient terjadi pada peristiwa yang timbul dalam kecelakaan Pembangkit Listrik Tenaga Nuklir (PLTN) tipe reaktor air bertekanan ringan atau Pressurerized Water Reactor (PWR). Kecelakaan yang dimaksud adalah kecelakaan kehilangan pendingin atau Lost of Coolant Accident (LOCA). Salah satu hasil kepatutan analisis evaluasi desain keselamatan PWR selama pasca LOCA (postLOCA) yang disebabkan oleh kejadian pecahnya pipa utama adalah temperatur kelongsong bahan bakar (berbentuk silinder berongga) mencapai harga maksimum sekitar 930°C dan berada di bawah temperatur kritisnya 1200 °C [2J • Berdasarkan latar belakang tersebut di atas maka permasalahan yang akan dibahas dalam penelitian ini adalah bagaimana distribusi temperatur pada batang dipanaskan berbentuk silinder berongga menggunakan metode finite difference. Pada penelitian ini, perhitungan distribusi temperatur dilakukan berdasarkan proses pendinginan pada bagian uji Heating-OJ. Perhitungan solusi numerik dalam kasus ini dibatasi hanya pada persamaan panas 1 dimensi untuk koordinat tabung. Tiga solusi numerik akan digunakan dengan metode aproksimasi beda hingga (finite difference method), yaitu metode FTCS, BTCS, dan Dufort Frankel.
Proses pendinginan pada peristiwa konduksi diindikasikan dengan terjadinya perpindahan panas dari daerah temperatur tinggi ke daerah temperatur rendah, dengan gradien temperatur tertentu dalam kurva perubahan bergantung terhadap waktu. Dalam peristiwa tersebut, enegi yang dipindahkan secara konduksi dengan angka perpindahan panas persatuan luas berbanding lurus terhadap gradien temperatur normal, yaitu aT Ox
Saat konstanta maka: 120
or
q
=-kA-
(1)
Ox
Dengan q adalah angka perpindahan panas dengan satuan kJ per satuan dan DT/ox adalah gradien temperatur dalam arah aliran panas. Konstanta k positif disebut konduktivitas termal benda, dengan satuan W/moK. Tanda kurang disisipkan agar prinsip termodinamika kedua terpenuhi, yaitu panas harus mengalir turun pada skala temperatur [3]. Pada kasus umum yang temperatumya bergantung pada waktu dan terdapat sumber panas pada benda. Untuk elemen ketebalan dx dapat dibuat keseimbangan energi sebagai berikut : Perubahan da1am energi
Energi terkonduksi kedalam
internal
Permukaan kiri
I =
+ Panas yang dibangkitkan
+ Energi terkonduksi
ke luar
Permukaan kanan
Dalam elemen
Secara matematis ditulis sebagai : " or qx + qAdx = pcA-dx + dX+dx
(2)
m
Karena q
= _ kA aT ax
(3)
Dan
qx+dx=-kA :t+dX (4)
maka persamaan (2.2) dapat ditulis sebagai
dx dx +q=pc~u n ( k or)" or
TEORI
~_ a
NUKLIR
don Proses Bahan
k
yang sebanding disisipkan,
a2T q I or --+-=-dx k am
(5)
Dengan : p = kerapatan (kg/m3) c = panas spesifik (kJ / kgK) a = klpc = difusivitas termal bahan (m2/s) Bahan berbentuk hal/ow cylinder. maka sistem koordinat tabung lebih sesuai digunakan. Berdasarkan[4], asumsi dan idealisasi yang diterapkan adalah I. Bahan uji HeaTing-O I bersifat homogen sehingga konstanta difusivitas termal a tidak bergantung terhadap kedudukan r. 2. Perambatan panas hanya pada tebal pipa, sehingga temperatur T hanya bergantung pada
ISSN 1410 - 8178
Ainur Rosidi, dkk
PROSIDING SEMINAR PENELITIAN DAN PENGELOLAAN PERANGKAT NUKLIR Pusot Teknologi Akselerotor don Proses Bahan Yogyakarta,
posisi dan waktu, yaitu r dan I, yang secara matematis dapat ditulis sebagai T = T(r,t). 3. Tebal pipa tetap sekalipun terjadi perubahan temperatur. 4. Kedua sisi r = a dan r = b hanya dihadapkan pad a masalah radiasi. 5. Tidak ada heal generalion. Oleh karena itu, persamaan diferensial konduksi panas 1 dimensi adalah
--+--=-ar2 r ar a a2T
1 aT
1 aT iTt
dalama
0 T(r,t) = To. pada t = 0 TATA KERJA Tata ketja yang dilakukan adalah mengamati fenomena fisika distribusi panas pada batang pemanas yang dipanaskan secara radiasi hingga temperatur tertinggi mencapai 850°C melalui DAS yang terukur dari 14 tennokopel yang terpasang di batang pemanas. Serta menyederhanakan pennasalahan yang kemudian membuat persamaan matematika dari fenomena tersebut, yaitu membuat persamaan diferensial konduksi panas satu dimensi dan syarat batasnya. Kemudian menampilkan hasilnya dalam bentuk tabel dan grafik menggunakan Matlab 6.5 dan software grafik. Peralatan Eksperimen Gambar Iengkap desain alat eksperimen yang diberi nama: bagian uji HeaTiNG-OI dapat dilihat pada Gambar 2. Komponen Utama (Bagian Uji HeaTiNG-I) : I. Plenum atas (tempat menampung air)
28 Agustus 2008
2. Tabung gelas kuarsa (p=IOOOmm, OD=45mm, 10=41
rom)
3. Batang pemanas yang merupakan simulasi debris untuk geometri anulus, material yang digunakan adalah SS3I6 dengan panjang 1100 mm (healed lenglh = 800 mm). Kemudian 14 buah termokopel tipe K dipasang pada pennukaan bagian luar batang pemanas yang digunakan untuk mengukur perubahan temperatur pennukaan batang pemanas selama pendidihan berlangsung. Gambar 4 menyajikan posisi tennokopel yang telah telah dipasang. Kompenen lainnya adalah flange-flange dan material pengikat antara tabung gelas kuarsa dengan batang pemanas. Komponen Pendukung I. Komponen Listrik Bagian uji Heating-I dipanaskan secara radiasi oleh 2 pasang semi-silinder keramik heater dengan daya total 20.000 watt. Selain itu heater pemanas untuk air dipasang pada plenum atas untuk air pendingin yang akan dimasukkan ke dalam celah sempit. Slide regulator voitage dengan daya maksimal 25.000 watt digunakan untuk mengatur masukan tegangan selama pemanasan berlangsung. Gradual kenaikan daya diperlukan agar distribusi panas dapat merata bagian uji dan dapat menghindari thermal shock pada tabung gelas kuarsa. 2. Komponen Instrumentasi Dalam eksperimen ini yang digunakan adalah DAS yang dimiliki Laboratorium Tennohidrolika (Dataq Instruments, USA), dengan laju perekaman data 5 dataldetik untuk setiap kanal dari 24 kanal yang diterpasang. .- --+ -t'J4i_
,
,- •••no --.
0-1111"""-
---+
(p...-.C8Ic)
1 termokopel
3te.mokDpet (•• dill/)"" ----..
1 termohDpei
--
UermokDpei 1 termolmpel
~.f"A
...o,.M_
1 termokDpef
---...
1 termolcDpd
3 posisi radial
f"'''- --3te.mokDpel Ue.mokDpei 1""(•• diall i
./
/
•
I
l
Gambar I. Foto bagian uji HeaTiNG-OI
Ainur Rosidi, dkk.
bill
.•....• 1 termokopel (punc.k)
;~
Gambar 2. Posisi 14 tennokopel pada batang pemanas
ISSN 1410 - 8178
\2\
PROSIDING SEMINAR PENELITIAN DAN PENGELOLAAN PERANGKA T NUKLIR Pusat Teknologi Akselerator dan Proses Bahan Yogyakalta, 28 Agu~us 2008 HASIL DAN PEMBAHASAN Perhitungan
Oari gambar 3 diperoleh bahwa saat t = 19,8 detik temperatur bahan mendekati temperatur lingkungan. Hasil perhitungan menunjukkan bahwa pada waktu tersebut temperatur pada node x = 2 mm, x = 4 mm, x = 6 mm, dan x = 8 mm sebesar 28,0909 DC, 28,155 DC, 28,1641 °c dan 28,1083 DC.
Numerik
Berdasarkan
asumsi
(6)
pada persamaan
menjadi aT -+--r aTar =-a at a2T
1
1
ar2
Metode BTCS
dalam 0,0105 < r < 0,0205
(7)
Oengan : T(O,t) = T(O,O 1,t) = 28° C pada t > 0 T(r,t) = 850° C pada t = 0 Untuk kesamaan antara metode yang digunakan, bahan dibagi menjadi 5 bagian dengan 0i:; U u 6 node dengan 8 = 0,002 m, maka akan terbentuk e
Untuk metode BTCS, L1t tidak bergantung pada kriteria stabilitas. Maka dengan 8 = 0,002 mm, L1t= 0,45 detik dan nt = 60 seperti pada metode FTCS, distribusi perpindahan panas transien ditunjukkan oleh gambar 4.
~
B
besar temperatur T), T2, T3, T4, Ts dan T6 pada waktu ke-n. Pembagian ini dimaksudkan agar tebal bahan antara node tidak menjadi semakin tipis. Juga ingin melihat selisih temperatur pada node yang simetris. Sedangkan untuk pembagian waktu disesuaikan dengan metode yang digunakan.
-0-
::00 f-o 100 900
5<>0 1000 >00 700 •• 0 ••0
-¢--Nodt x. 8 mm -~-Nodt: x-IO_ -6-Nodt x· 4_x· 6 mm -V-Node
Node x. 1 Ibm.
-a-Node
x - 0 mm
".~
Metode FTCS Untuk metode FTCS, L1tditentukan kriteria stabilitas
oleh •
6
10
1::
I~
16
1e
::0
::
:t
::0:.
Waktu t[ detik]
..,+{r] <
Gambar 4. Grafik distribusi temperatur dengan metode BTCS
2(2xlO-3mY
- 4'295XIO-6m2/{4+C'05Xl~~2m
2xI0
Metode Dufort Frankel m
Untuk metode Dufort FrankEl, skema akan konsisten jika L1t < 8. Jika 8 = 0,002 mm, maka L1t < 0,002 detik. Untuk .M = 0,0018 dan nt 15000, distribusi perpindahan panas ditunjukkan oleh gambar 5.
+5r)
:<>0,4615
maka
kriteria
stabilitas
untuk
L1t harus
memenuhi L1t:S 0,461 detik Untuk L1t= 0,45 detik nt = 60, distribusi dan increment waktu 0U oleh perpindahan panas ditunjukkan gambar 3. 500 1000
0
5:; u ~,.., !'! ~
"""
000 000
'00
~:J\_~
i
i:; 100 20~ 90.
-v-Nodlt -Q-Nodit =8 mm. --tNode Xx-6mm = 10xnm -A-Nodtx-4mm
'00 1000 JO. "'09 ,"0
-O-Nodl!: x=2mm
f-o
-a-Node
X
=
0 mm
~
-a-Node -o-Nod! -.6.-No& -v-Node
x· 0 mm x= ~ n:m X· 4 mm x= 6mm
-<>-No1h
x= 8 mm
"1
--4--Nod. x= lOmm
111111111111111111111
Waktu
W.ktu t [detik]
Gambar 3. Grafik distribusi temperatur dengan metode FTCS
20
t
[detik]
Gambar 5. Grafik distribusi temperatur dengan metode Dufort FrankEl
20
I
21 ,15 122
20
ISSN 1410 - 8178
Oari gambar 5 diperoleh bahwa saat = detik temperatur bahan mendekati Ainur Rosidi, dkk
PROSIDING SEMINAR PENELITIAN DAN PENGELOLAAN PERANGKAT NUKLIR Pusat Teknologi Akselerator don Proses Bahan Yogyakarta,
28 Agustus
temperatur lingkungan. Hasil perhitungan menunjukkan bahwa pada waktu terse but temperatur pada node x = 2 mm, x = 4 mm, x = 6 mm, dan x = 8 mm sebesar 28,0997 DC, 28,1698 DC, 28,1799 DC dan 28,1183 Dc. Pada kasus ini, terlihat dari gambar 3, gambar 4 dan gambar 5 evaluasi temperatur pada node x = 0,0 mm dan x = 10,0 mm tidak mengalami pola penurunan secara eksponensial, hal tersebut menunjukkan bahwa dalam simulasi ini bagian luar dianggap adiabatik (tidak ada perpindahan panas). Namun, pada node x = 2,0 mm, x = 4,0 mm, x = 6,0 mm dan x = 8,0 mm evaluasi temperatur terjadi sebagaimana mestinya. Trend penurunan temperatur secara transien mewakili bent uk eksponensial negatif.
2008
diperoleh dengan memfitting kurva pada Gambar 6, sehingga diperoleh korelasi sebagai berikut: T(l) = 41,4
TC-1
+ 42l2e-oI"'u
+ 385.7e-"'·~·
. T(l) = 41,2 + 425,2,,-<'19Jl2 + 385.1,,-
TC-3 TC-5
: T(l) = 41,9 + 448,?.-""·~'
TC-6 : T(l) = 38,4 + 463,2,,-"29<2.82 TC-7 : T(l) = 8,3'" TC-8
: T(l) = 36,8
339,6e-
+ 380,3e-""'5L4
+ 365.1e-,m.~
+ 355.8,,-'114'., + 327.3<:>-<113:11 + 119,2e-
Perbandingan Trend Hasil dengan Hasil Eksperimen
Numerik
Perbandingan trend antara ketiga metode numerik dengan eksperimen pada bagian uji HeaTing-OI diperlihatkan pada gambar 7. 1000
Hasil Eksperimen Pengujian pemanasan dilakukan untuk menentukan kemampuan pemanasan semisilinder keramik dan kemampuan struktur bagian uji HeaTiNG-OJ. Uji pemanasan dilakukan dengan memanaskan batang pemanas secara radiasi hingga temperatur tertinggi pada batang pemanas mencapai 850°C
900
"
'L E
~ iJ
IIh
0 E" 2110 100 0
700 600 800 100 300 500 Metode
FTCS, node x ~ 6 mm
Data Ekspt:rimtn
ODD
U L 800
Im_
-9000 dt:lik
E-.
Iij"
~ CCI
1'(1)= 41.2+ 4Z~.R·fI'w.a.;I-+ :5&,,1
Iij
.g
500
c::
400
co
l(·-3 T(,$
T(l)-" 41.9~ 442 ..'V~.-.ecu .•. 3~:$.~·"AU
T(>6
T(t) •••)8,4
T(t} '" 41.0
,.
+
43~.w~$"'l.:.
'"
TH)
=
36.3
M5.!lo·lrnt,)
463.a·v•u.lJ
T(I) ••3.3+3J!).~""~
.:9 t:Q
T.,.,._OODC
T(f) = 414. 42:1.~-(';Jto)u + 3t~.74~'~
T(';.1
700
-t.
),n.3t"';ru;
+119.2•.••
Ao<+'
.•. J&'U.Jc:o-flWH4
•
11
1'! 200 V '00
E-<
1000
2000
3000
4000
Waktu,
5000
0000
7000
8000
gOOD
t [detik]
Gambar 6. Kurva penurunan temperatur batang dipanaskan terhadap waktu Gambar 6 menunjukkan kurva penurunan temperatur batang yang dipanaskan terhadap waktu kemudian didinginkan secara radiasi (tanpa pendinginan air). Pada proses pendinginan radiasi ini, pemahaman terhadap karakteristik pendinginan radiasi sangat diperlukan untuk mengetahui berapa lama pendinginan radiasi berlangsung dan bagaiamana bentuk kurva penutunan temperatur transiennya. Berdasarkan Gambar 6, meskipun temperatur awal titik termokopel yang terpasang sepanjang batang dipanaskan memilki temperatur yang berbeda, namun pada detik ke 9000 (2 jam, 30 menit) temperatur pad a setiap titik termokopel hampir sarna (sekitar 50°C). Korelasi yang bisa mendekati kurva penurunan temperatumya Ainur Rosidi, dkk.
18
21
21
27
Waltut[detik]
300
!=!
It v
15
••
28t.5t ..•. '*.1
Gambar 7. Grafik perbandingan distribusi temperatur 20 antara hasil eksperimen dan metode numerik bagian uji HeaTing-01 Perbandingan antara hasil perhitungan numerik dengan eksperimen tidak dimaksudkan untuk membandingkan waktu realistik atau perbandingan time-to-time. Perlu membagi hasil dengan bilangan tertentu agar rentang waktu realistik sarna dengan rentang waktu dalam perhitungan di dalam kurva. Perbandingan yang dimaksudkan untuk melihat trend perubahan temperatur terhadap waktu. Pada bagian uji HeaTing-OI, trend hasil numerik mengikuti trend hasil eksperimen. ketiga metode memperlihatkan perbedaan yang kecil. Metode BTCS dan Dufort Frankel dianggap paling mendekati trend hasil eksperimen. KESIMPULAN Eksperimen perpindahan panas pada celah sempit anulus yang dilakukan dalam mempelajari dan memahami karakteristik distribusi panas pada batang bagian uji HeaTing-
ISSN 1410 - 8178
123
PRO SIDING SEl\1INAR PENELITIAN DAN PENGELOLAAN PERANGKAT NUKLIR Pusat Teknologi Akselerator don Proses Bahan Yogyakarta, 28 Agustus 2009 01 telah dilakukan melalui analisis kurva perubahan kenaikan temperatur batang dipanaskan terhadap waktu. Hasil simulasi diperoleh distribusi temperatur (berupa tabel dan grafik) menggunakan metode FTCS, BTCS, dan Dufort Frankel dengan waktu proses pendinginan mencapai temperatur lingkungan secara radiasi yang sangat cepat, yaitu kurang dari 22 detik (maksimal 21,6 detik). Sebagai evaluasi perbandingan, waktu proses pendinginan hasil eksperimen mencapai 9000 detik. Trend distribusi temperatur proses pendinginan secara radiasi dari hasil numerik mengikuti trend hasil eksperimen. Dari ketiga metode yang digunakan, trend metode Dufort Frankel dan BTCS yang dianggap lebih mendekati trend hasil eksperimen. Trend penurunan temperatur secara transien mewakili bentuk eksponensial negatif. UCAPAN TERIMAKASIH Ucapan terima kasih kami sampaikan kepada ternan-ternan di sub bidang termohidrolika BOFa PTRKN atas bantuannya selama eksperimen. Kepada Ka. BOFa atas bimbingannya. Penulis mengucapakan terimakasih atas dukungan dana melalui DIPA KNRT tahun anggaran 2007 (SK. Menristek No. 126/M/Kp/XI/2006 tanggal 17 Nopember 2006, perihal Program InsetifRiset Dasar KNRT 2007). DAFT AR PUST AKA 1. CENGEL,YUNUS A., and R. H. TURNER, Fundamendals of Thermal-Fluid Sciences, McGraw-Hili Companies,lnc.,1997 2. JUARSA,MULYA and ANHAR R. ANT ARIKSA W AN, Studi Perpindahan Panas Selama Rewetting pada simulasi Pendinginan Pasca Loca, Yogyakarta Seminar Pertemuan dan Presentasi I1miah
124
Penelitian I1miah Penelitian Dasar fPTEK Nuklir, 2006. 3. HOLMAN, J.P,. Heat Transfer Seventh Edition in Sf Units, McGraw-Hili Companies, Inc., 1992 4. PONIDI, dkk, Pendahuluan Persamaan Diferensial Parsial dan Syarat Batas, Diktat Kuliah PDP & SB Dept. Matematika FMIPA UI,2001
TANYA JAWAB Agus Dwiatmaja ~ Dari ketiga metode analisis FTGS, BTCS dan Dufort Frankel, menurut anda metode apa yang mudah digunakan tetapi mempunyai hasil yang akurat? Ainur Rosidi ~ Metode yang mempunyai hasi/ yang akurat adalah metode Dufort Frankel karena trend perubahan temperatur mendekati trend dari penurunan temperature eksperimen. Kussigit S. ~ Apakah ada perbedaan dari ketiga metode tersebut? ~ Pendekatan dengan Deret Taylor, sampai berapajauh orde deret yang dipakai? Ainur Rosidi ~ Perbedaan dari ketiga metode Finite Difference adalah dengan melihat dari struktur penyebaran panas dan bentuk grafik trend penurunan distribusi panas eksponensial negatif. ~ Pendekatan solusi numerik untuk differensial parsial diselesaikan dengan Finite Difference deret Taylor orde 2.
ISSN 1410 - 8178
Ainur Rosidi, dkk
