KOMPUTASI ALIRAN FLUIDA DINAMIK DENGAN CITRA DIGITAL DAN PIV (PARTICLE IMAGE VELOCIMETRY), KHUSUSNYA DALAM APLIKASI NUKLIR Muhammad Arifin Sanusi*
ABSTRAK KOMPUTASI ALIRAN FLUIDA DINAMIK DENGAN CITRA DIGITAL DAN PIV (PARTICLE IMAGE VELOCIMETRY) KHUSUS NYA DALAM APLIKASI NUKLIR. Aliran fluida yang ditambahkan pembiak (seed), diamati dengan mengambil 2 citra secara berurut dalam selang ∆t . Kedua citra ini dihitung dengan metode FFT, korelasi-silang dan oto-korelasi, kemudian diinversi kembali menjadi fungsi waktu. Dengan mengetahui ∆t maka peta vektor kecepatan partikel aliran fluida dapat dihitung. Hasil fenomena aliran kecepatan fluida ini dapat diaplikasi lanjut pada penggunaan sangat luas di bidang teknik. Khususnya dalam aplikasi nuklir, aliran fluida, turbulensi dalam reaktor, dapat diketahui tanpa memasang dan menyentuh, karena penggunaan laser, kamera CCD, serta perangkat komputasi kendali jauh.
ABSTRACT COMPUTATION OF DYNAMIC FLUID FLOW USING DIGITAL IMAGE AND piv (PARTICLE IMAGE VELOCIMETRY) PARTICULARLY IN NUCLEAR APPLICATION. Two images captured sequentially in ∆t period, by CCD camera and double-pulse laser in cross-section of flow with seeding particle, were reviewed. The digitalized images were calculated by FFT method, crosscorrelation and auto-correlation as an impulse function equation, therefore inversed in to time domain resulting velocity vector maps with known ∆t . The computation result can be applied in a large technical application, especially in nuclear application, turbulence and the other phenomena fluid dynamic in reactor chamber which can be calculated providing Laser, CCD camera and remotely computational equipment are available.
PENDAHULUAN Citra video digital dapat dilihat sebagai medan sinyal analog dua dimensi di mana dapat dijadikan bentuk digital deret satu dimensi terhadap waktu t. Teknik pengolahan *
BPP Teknologi
sinyal satu-dimensi dapat secara langsung dikembangkan ke dua-dimensi, seperti yang telah dibahas oleh Rosenfeld & Kak 1976; dan Pratt tahun 1984. Pada metode Analisa Citra Digital (ACD) - PIV yang dipakai dalam komputasi ini prinsipnya adalah: 2 citra digital berurut (sekuensiel) dicacah pada satu daerah melalui introgasi jendela (window) - Gambar 1. Dalam cacah citra ini suatu pergeseran rata-rata partikel dapat diamati dari satu cacah dalam citra itu terhadap pasangannya pada citra lainnya. Pergeseran spasial itu dapat dijelaskan sederhana dengan suatu model pengolahan sinyal digital yang linier, seperti dalam Gambar 2. Satu daerah cacah f(m,n) dapat dipakai sebagai masukan ke suatu sistem yang mempunyai keluaran g(m,n) yang sesuai dengan daerah cacah dari citra lainnya yang diambil citranya pada ∆t kemudian. Sistem itu sendiri terdiri atas dua komponen, fungsi pergeseran spasial s(m,n) yang juga dikenal sebagai sistem impuls jawab, dan sinyal derau tambahan d(m,n).
f(m,n) citra 1
g(m,n) citra 2 t0 + ∆t
s’(m,n) Fungsi geser dari medan yang di estimasi
t0
Gambar 1. Konsep susunan bingkai-ke-bingkai pencacahan ACD –PIV Sinyal derau ini disebabkan oleh perpindahan partikel daerah cacah partikel yang tidak nampak pada gerak 3 dimensi di lembar laser, dan jumlah total partikel yang ada dalam jendela. Tentu saja cacah awal f(m,n) dan g(m,n) juga berderau. Masukan citra 1 f (m,n) F(u,v)
Keluaran citra 2 s (m,n) S(u,v)
g’(m,n) G’(u,v)
+
d(m,n) D(u,v) Sinyal derau /noise tambahan
g(m,n) G(u,v)
Gambar 2. Model pengolahan sinyal dengan uraian hubungan fungsional antara 2 bingkai suksesif yang memuat partikel, fungsi F(u,v), S(u,v), G’(u,v), G(u,v) & D(u,v) sebagai transformasi Fourier dari fungsi huruf kecilnya yang menyatakan domain frekuensi spasial. Apabila kita abaikan sinyal derau tambahan yang mempengaruhi keluaran g(m,n) maka perhitungan akan lebih mudah dinyatakan sebagai fungsi dirac. Tentunya pendekatan ini diambil untuk menyederhanakan perhitungan yang demikian kompleks, karena kondisi ideal fungsi masukan f(m,n) dan g(m,n) juga akan mengandung sinyal derau sehingga perhitungan akan semakin kompleks. Pendekatan perlu diambil dalam menyelesaikan dengan metode konvolusi - transformasi Fourier.
PENERAPAN KOMPUTASI ACD - PIV. Gambar 3 menjelaskan metoda ACD -PIV yang dipakai dalam analisa ini, di mana pemakaian persamaan konvolusi - transformasi Fourier untuk mempercepat proses perhitungan korelasi-silang dan auto-korelasi, dibanding dengan metode lainnya yang membutuhkan waktu perhitungan yang demikian lama serta ketelitian yang kurang, dengan tranformasi Fourier operasi dapat disederhanakan dengan perkalian konjugasi bilangan kompleks setiap pasangan koefesien Fourier yang bersesuaian. Kumpulan kofesien baru ini selanjutnya dengan transformasi inversi diperoleh φfg . Oleh karena periodisitas FFT dalam ruang tidak diperlukan normalisasi φfg seperti pada kasus linear. Pada kenyataannya, meskipun pergeseran ini sering terlalu besar, komputasi ini bekerja baik, karena nisbah derau pada korelasi-silang akan menurun dengan naiknya pergeseran spasial. Itu berarti angka pasangan citra partikel menurun di daerah pencacahan dan akan lebih banyak partikel citra yang tak berpasangan. Bila diberikan sisi jendelaF(u,v N, kita memperoleh 1/3 sisi ini (N/3) pada batas yang dapat f(m,n Korelasi-silang diperoleh pegeseran vektornya. Ini secara langsung sesuai dengan kriteria PIV FFT ) Φ'( u, v) φ'( m,n ) dx,dy
g(m, n)
FFT
FFT *1
G(u,v
Φ'(u,v) = F(u,v)G*(u,v)
Gambar. 3 Metoda ACD-PIV
vx(i,j) vy(i,j)
Puncak korelasi-silang mulanya ditentukan dengan mendapatkan harga tertinggi dalam nilai korelasi matriks 2 dimensi itu. Untuk ketelitian sub-piksel, sekitar elemen ini suatu kurva parabolik atau eksponensial yang memenuhi arah horisontal dan vertikal dari lokasi mendekati, yang mempunyai puncak korelasi. Dibanding dengan metode konvensional dengan pusat-massa sentroid, tiga kurva yang sesuai yang dipilih. Beberapa urutan pengujian disampaikan dalam makalah ini juga dimonitor dengan teknik sentroid, dan ini memperlihatkan bahwa kurva eksponensial 3 titik sesuai, dan dapat di tentukan puncak korelasi silang dengan kesalahan yang kecil dibanding dengan teknik sentroid pusat-massa. Gejala ini dapat memberikan secara kasar bentuk Gaussian dari puncak korelasi-silang itu sendiri. Besaran tertinggi kecepatan partikel yang dapat di deteksi citra digital sekarang ini pada kecepatan 30 HZ, dengan memberikan ukuran jendela pencacahan N x M. Pergeseran maksimum partikel yang dapat diukur dalam aliran dapat ditentukan dengan mengalikan pergeseran piksel dengan besaran faktor antara bidang citra dan bidang objek. Dengan membagi pergeseran dengan konstante waktu δ t antara citra tangkapan (1/30 detik), maka kecepatan maksimum aliran yang diamati dapat ditentukan. Nisbah (ratio) pencacahan jendela dapat dirubah untuk memenuhi suatu aliran yang skala besar.
ALGORITMA PEMROGRAMAN Citra video yang telah direkam pada kaset video VHS dari percobaan meja optik dengan laser dan sirkuit fluida air yang telah diberikan partikel pembiak, dipindahkan ke komputer. Dengan menggunakan video player, TV dan frame grabber Matrox Pulsar yang dihubungkan dengan komputer serta dengan menjalankan software Matrox Lite (Intelcam) memungkinkan citra itu dapat direkam gambar per gambar sebagai file bitmap. Pada percobaan di meja optik telah diketahui debit air, faktor pembesaran piksel dan laju frekuensi pencacahan. Untuk komputasi secara on line pada sirkuit fluida, digunakan kamera CCD, Laser yang terintegrasi dengan komputer sehingga rekaman gambar dua citra berurut dipindakan (scan) dan hasil komputasi vektor aliran dapat secara langsung diperoleh. Tetapi untuk sementara hanya dilakukan secara batch processing seperti yang dijelaskan di atas. Berikut ini digambarkan diagram alir secara umum.program perhitungan ACPIV.
ALGORITMA MATEMATIS Transformasi fourier diskret dua dimensi Apabila kita nyatakan suatu citra yang ditangkap oleh CCD dan ditampilkan pada komputer dengan representasi gambar sederhana berikut: n2
n1
Gambar 4. piksel urutan priodik Deret Fourier diskret dua dimensi dan koefisien deret Fourier dapat dinyatakan secara langsung (Ionis Pitas, Digital Image Processing Algoritms) sebagai berikut:
X ( k1 , k 2 ) =
N1 −1 N 2 −1
∑ ∑ x (n , n ) exp( −i 1
n1
x ( n1 , n 2 ) =
1
− i 2πNn21k 2 )
(1)
n2
N1 −1 N 2 −1 1 N 1N 2
2π n1k1 N1
∑ ∑ X (k , k ) exp(i 1
n1
n2
1
2πn1k1 N1
+ i 2πNn21k 2 )
(2)
Persamaan (1) dan (2) di atas mendefinisikan transformasi Fourier dua dimensi (2-D DFT), di mana transformasi diskret ini yang digunakan sebagai algoritma pengolahan dan analisa sinyal diskret dua dimensi di komputer. Untuk lebih menyederhanakan persamaan di bawah ini dituliskan kembali pasangan DFT di atas dalam bentuk:
X ( k1 , k 2 ) =
N1 −1 N 2 −1
x( n1 , n1 )W Nn k W Nn k ∑ ∑ n n 1
x (n 1 , n 2 ) =
1 1
2 2
1
2
(4)
2
N1 −1 N 2 −1 1 N1N 2
X (k1 , k1 )W N−n k W N−n k ∑ ∑ n n 1 1
1
1
2 2
2
(5)
2
di mana : W N = exp( − i 2Nπ ) j
j
j = 1,2
Oleh karena dalam menghitung korelasi silang dan auto-korelasi akan melibatkan dua matriks dari dua citra yang diberikan sesuai yang dijelaskan pada bagian di atas, maka dalam menyele saikan komputasi dengan konvolusi linier.
y( n1 , n2 ) = x (n1 , n2 ) ⊗ ⊗ h(n1 , n2 ) ⇔ Y(k1 , k 2 ) = X (k1 , k 2 ) H ( k1 , k 2 ) y( n1 , n2 ) = IDFT [DFT [x( n1 , n 2 )]DFT [h( n1 , n 2 )]] di mana: IDFT DFT
= = x ( n1 , n 2 ) = h( n1 , n2 ) =
(6) (7)
Inversi Discrete Fourier Transform Discrete Fourier Transform matriks pada daerah citra 1 matriks pada daerah citra 2
Sehingga perhitungan konvolusi dari pers 5 di atas dapat dilakukan sbb. suatu daerah dipilih sebagai RN1 N 2 di mana N1 ≥ L1 , N 2 ≥ L2 , urutan x ( n1 , n 2 ) ,
h( n1 , n2 ) diberikan nilai awal nol pada urutan x p (n1 , n2 ) , h p (n1 , n2 ) pada daerah RN1 N 2 : x(n1 , n2 ) xp (n1 , n2 ) = 0 h(n1 , n2 ) hp (n1 , n2 ) = 0
(n1 , n2 ) ∈ RP1P2 (n1 , n2 ) ∈ RP1 P2 (n1 , n2 ) ∈ RP1P2 (n1 , n2 ) ∈ RP1P2
(8)
(9)
Dengan mudah dibuktikan bahwa konvulasi :
yp (n1, n2 ) = xp (n1 , n2 ) ⊗⊗hp (n1, n2 )
(10)
memberikan hasil yang sama dengan pers. (6) konvulasi linier, hasil konvulasi linier diberikan oleh:
y(n1 , n2 ) = y p (n1 , n2 )
(n1 , n2 ) ∈ RL1L2
(11)
Algoritma perhitungan konvulasi dapat diringkaskan sbb. Pilih N1 dan N2 pada daerah citra tampilan di monitor, di mana Ni = 64,32, 16,8, i =1,2 2. Pemberian nilai awal matriks x ( n1 , n2 ), h( n1 , n2 ) dengan nol 3. Hitung DFT dari x p ( n1 , n 2 ), h p (n1 , n2 ) 1.
4.
Hitung DFT Yp ( k1 , k 2 ) sebagai hasil kali X p ( k 1 , k 2 ) * H p ( k1 , k 2 )
5.
Hitung y p (n1 , n2 ) dengan menggunakan inversi DFT
PERHITUNGAN KORELASI SILANG DAN AUTOKORELASI Apabila kita menuliskan kembali transformasi fourier yang dipakai untuk menghitung korelasi dua dimensi: N1 −1N2 −1
Rxy(m1,m2 ) = ∑∑x(n1, n2) y(n1 + m1 ,n2 +m2 )
(12)
n1 =0 n1 =0
N1 −1N2 −1
Rxx(m1,m2 ) = ∑∑x(n1, n2)x(n1 + m1,n2 + m2 ) n1 =0 n1 =0
(13)
Rxy ( m1 , m2 ) dipakai untuk perhitungan korelasi silang sedang Rxx (m1 , m2 ) pada auto-korelasinya. Citra x dan y pada daerah yang masing-masing memenuhi R P1 P2 , RQ1Q2 , selanjutnya metode yang dipakai untuk perhitungan korelasi 2-d adalah: Rxy (n1, n2 ) = IDFT[DFT[x(n1, n2 )]* DFT[ y(n1, n2 )]]
Appendiks: Beberapa contoh hasil komputasi
Contoh 1. hasil komputasi ACD - PIV
(14)
Contoh 2. dari hasil perekaman terakhir
KESIMPULAN Dari hasil algoritma matematis yang telah diuraikan, dan dideklarasikan dalam pemrograman di komputer, komputasi vektor kecepatan aliran fluida yang telah dilakukan belum dapat dicapai secara langsung (on-line) di sirkuit fluida yang ada. Kendala adalah pengaturan kecepatan aliran fluida yang ditambahkan pembiak itu mengalami keterbatasan sirkuit, peralatan dan instrumentasi, sehingga gambar lembar aliran dengan sinar laser direkam dari CCD, dengan TV dan video recorder vhs, setelah itu dipindahkan ke komputer yang selanjutnya gambar bitmap diolah dengan program tersebut. Hasil komputasi masih perlu validasi lebih lanjut meskipun perhitungan telah mendekati kebenaran.
DAFTAR PUSTAKA 1. STEVEN C CHAPRA, RAYMOND P CANALE. “Numerical Methodes for Engineer”, Mc. GRAW – HILL, International Edition. 2. MARKUS RAFFEL, CHRISTIAN E. WILLERT, JURGEN KOMPENHANS “Particle Image Velocimetry: A Practical Guide (Experimental Fluid Mechanics)“, Springer. 3. IONNIS PITAS, “Digital image processing algorithms”, Prentice-Hall, Inc. Upper Saddle River, NJ, USA, 1993.
DAFTAR RIWAYAT HIDUP
1. Nama
: Muhammad Arifin Sanusi
2. Tempat/Tanggal Lahir
: Bone, 21 Desember 1955
3. Instansi
: Puspiptek
4. Pekerjaan / Jabatan
: Pegawai/ Staf Bidang Perencanaan Puspiptek
5. Riwayat Pendidikan
:
• 1976-1984, Jurusan Elektroteknik, Fakultas Teknik Universitas Hasanuddin • 1987, (S2),USTL-Montpellier-DESS Informatika • 1992, Teknik Komputer, Ecole d’Engineur de Marseille-Perancis 6. Pengalaman Kerja
:
• 1985, LAPAN • 1985-Sekarang, BPPTeknologi • 1988-sekarang, diperbantukan di ASDEP Puspiptek 7. Organisasi Profesional
:
• 1985-Sekarang, Anggota PPI 8. Makalah yang pernah disajikan : • Dalam LKSTN tahun 1995. • Dalam Seminar LTMP-BPPT • Dalam Lokakarya LTMP-BPPT
Back
