R.A. Sani *), M. Barmawi, T.Winata, Sukirno Lab. Semikonduktor, KBK Fisika Material Jurusan Fisika FMIPA-ITB

APLIKASI METODA ELEMEN HINGGA DENGAN ALGORITMA CRANKNICOLSON DALAM PENYELESAIAN PERSAMAAN DIFUSI PANAS UNTUK PEMBUATAN SEL SURYA SILIKON AMORF (A-SI) TERPADU1)

R.A. Sani*), M. Barmawi, T.Winata, Sukirno Lab. Semikonduktor, KBK Fisika Material Jurusan Fisika FMIPA-ITB

Abstrak Penyelesaian persamaan difusi panas yang didasarkan pada perhitungan numerik dengam metoda elemen hingga digunakan untuk menentukan sel surya. Program ditulis dalam bahasa Fortran 77 dengan algoritma Crank-Nicolson dan dibantu dengan perangkat lunak dari library IMSL untuk menyimpan elemen matrik kekakuan dalam modus pita dan untuk menyelesaikan persamaan matrik kekakuan struktur. Diskritasi struktur dianalisa dalam elemen-elemen cincin segitiga simetri sumbu. Konfigurasi dari hasil perhitungan diperlihatkan untuk 50, 98, dan 140 buah elemen. Diperlihatkan bahwa hasil perhitungan distribusi temperatur pada saat sinar laser untuk tiga lapisan pertama sel surya a-Si terpadu, masih lebih tinggi dibandingkan dengan hasil perhitungan Kishi dkk [19].

Abstract The solution of thermal diffusion equation computed numerically by finite element method is used to determine optimum conditions of laser irradiations in the solar cell fabrication. Structure discretization is analyzed in axisymmetric triangular ring elements. The program are codes in written using Fortran 77 language. The Crank-Nicolson algorithm and Fortran 77 IMSL library were used in solving stiffnen matrix equations of the structure. The convergence of the results are shown 50, 98, and 140 elements. The temperature distribution of at the center of laser beam in the all three layers of a-Si solar cells is also shown and superior compared with that of Kishi et al [19]. *) 1)

Staf dosen FPMIPA IKIP Medan (alumni S2 Jurusan Fisika ITB angkatan 1993) Diperluas dari Makalah yang dipresentasikan dalam Seminar Metoda Elemen Hingga di Lab. Perancangan Mesin, Jurusan Teknik Mesin, FTI-ITB. 24

JMS Vol.1 No.1, April 1996

25

1. PENDAHULUAN Metode pembuatan pola untuk sel surya terpadu adalah sebuah faktor penting dalam mendapatkan sel surya yang murah. Biasanya sub-modul sel surya terpadu, dibuat dengan menggunakan sebuah topeng dari logam atau dengan teknik fotolitografi. Metode-metode ini mempunyai beberapa permasalahan, seperti : banyak sekali proses pembuatan, adanya pembatasan ukuran substrat, dihasilkan sel surya dengan luas efektif yang kecil, dan sangat tingginya kemungkinan terjadi lubang-lubang kecil pada lapisan-lapisan silikon amorf. Maka dikembangkan suatu teknik pembuatan yang baru dengan menggunakan sinar laser berdaya tinggi untuk membuat pola pada sel surya sehingga masalah-masalah yang disebutkan di atas, tidak muncul. Untuk mengatasi masalah kontrol daya laser, telah diteliti suatu metode baru dalam teknik pola dengan laser untuk pembuatan sel surya terpadu dengan efisiensi tinggi, yaitu metode LWS (Laser Welding and Scribing) yang tidak memerlukan kontrol daya laser yang canggih.

2. ANALISA KEADAAN PANAS DARI SINAR LASER

Sebagai bagian dari analisa panas, kita menentukan sifat panas dari tumpukan lapisan film tipis dengan memecahkan persamaan panas tiga dimensi. Dianggap bahwa sinar laser yang diradiasikan pada permukaan struktur tumpukan lapisan film tipis, berdistribusi Gauss. Rapat daya laser, I(r,t) yang dijatuhkan pada film adalah[19] : I (r , t ) =

⎛ 2r 2 ⎜− exp ⎜ w2 πw 2 ⎝

2 P(t )

⎞ ⎟ ⎟ ⎠

(1)

dengan : r adalah jarak dari pusat sinar, w adalah jari-jari berkas laser, dan P(t) adalah daya laser

Rapat daya laser yang diasorpsi W(r,z,t) dalam tumpukan lapisan (pada lapisan ke-k) bisa dinyatakan sebagai

26


⎛ k −1 ⎞ W (r , z, t ) = (1 − R) I (r , t )exp⎜ − (a1 − a k )d i + a k z ⎟ ⎜ ⎟ ⎝ t =1 ⎠

∑

(2)

dimana R adalah energi reflektansi. Sumber panas dalam Q(r,z,t) dinyatakan sebagai :

Q( q , z , t ) = −

∂W ( r , z , t ) ∂z

(3)

Persamaan difusi panas tiga dimensi bisa dinyatakan sebagai : ρC

∂T K ⎛ ∂ ⎛ ∂T = ⎜ ⎜r ∂t r ⎜⎝ ∂r ⎝ ∂r

⎞ ∂ ⎛ ∂T ⎞ ⎞ ⎟ ⎟⎟−Q(r , z , t ) ⎟ + ⎜r ⎠ ∂z ⎝ ∂z ⎠ ⎠

(4)

dengan ρ adalah kerapatan C adalah panas jenis, K adalah konduktivitas panas, dan T adalah suhu

Bila radiasi dan konveksi panas diabaikan (dianggap kecil sekali), maka syarat batas untuk persamaan (4) adalah : (a) T = T0, pada posisi dengan suhu tetap (b) pada posisi dengan fluks termal, berlaku : ⎛ ∂T ∂T ⎞ K⎜ + ⎟+ g 0 ( r , z , t ) = 0 ⎝ ∂r ∂z ⎠

(5)

Kondisi awal adalah T = φ0 (saat t = 0), dimana φ0 adalah suhu awal dan q0 adalah fluks panas.

3. PERUMUSAN ELEMEN HINGGA

Dalam masalah ini dipilih elemen cincin segitiga simetri sumbu, karena sifat termal bahan dan transfer panas adalah simetri terhadap sumbu z. Dalam tiap elemen segitiga, fungsi medan bisa didekati oleh polinomial orde pertama, yaitu :

φi ( R, Z ) = ai + bi R + ci Z

(6)

Besaran yang tidak diketahui, diwakili oleh sebuah variabel simpul φ yang dapat dinyatakan sebagai [6] :

φ (e) = N1φ1 + N 2 φ 2 + N 3φ 3 Ni adalah koordinat-koordinat alamiah untuk segitiga, yang diberikan oleh persamaan (6):

(7)


Ni =

27

ai + bi R + ci Z , i = 1, 2, 3 2∆

(8)

dengan

1 R1

Z1

2 ∆ = 1 R2 1 R3

Z2 Z3

(9)

Persamaan-persamaan elemen untuk persamaan (4) dan (5) dapat diturunkan dengan metode Galerkin6). Langkah pertama adalah menentukan pendekatan kelakuan temperatur dalam tiap elemen, yaitu : 3

T (e) ( R , Z , t ) = ∑ N i ( R , Z ) Ti ( t ) i =1

Kemudian dengan menerapkan kriteria Galerkin[6], persamaan (4) bisa ditulis :

∫∫D

K ⎡ ∂ ⎛⎜ ∂T ( e) R ⎢ ∂R R ⎣⎢ ∂R ⎜⎝

(e)

⎞ ∂ ⎛ ∂T (e ) ⎟+ ⎜ R ⎟ ∂Z ⎜ ∂Z ⎠ ⎝

(e) ⎤ ⎞ ⎟ + Q − ρc ∂T ⎥ 2πRdRdz = 0, (i = 1,2,3) ⎟ ∂t ⎦⎥ ⎠

(10)

Persamaan (10) menyatakan perata-rataan kesalahan yang diinginkan (atau residu) di dalam batas-batas elemen, tapi tidak memuat pengaruh dari batas. Solusi persamaan (10) dalam bentuk matrik adalah : ⎧ dT ⎫ [C ] = ⎨ ⎬ + [ K ]{T } = {RQ }+{Rq } ⎩ dt ⎭

(11)

dimana,

(R1 + R2 + R3 ) K (b b

k ij =

12∆

i j

+ ci c j

)

(12)

koefisien-koefisien b dan c, dan ∆ dinyatakan dengan persamaan (8) dan (9) ⎡6 R1 + 2 R2 + 2 R3 ρc∆ ⎢ [C ] = 12 ⎢ ⎢⎣simetri (e)

{RQ }

2 R1 + 2 R2 + R3 2 R1 + 6 R2 + 2 R3

2 R1 + R2 + 2 R3 ⎤ R1 + 2 R2 + 2 R3 ⎥⎥ 2 R1 + 2 R2 + 6 R3 ⎥⎦

(13)

⎡2 1 1⎤ ⎧ R1 ⎫ Q∆ ⎢ ⎪ ⎪ = 1 2 1⎥⎥ ⎨R2 ⎬ 12 ⎢ ⎢⎣1 1 2⎥⎦ ⎪⎩ R3 ⎪⎭

(14)

q 0 l ij ⎡2 1 ⎤ ⎧ ri ⎫ ⎢ ⎥⎨ ⎬ 6 ⎣ 1 2 ⎦ ⎩r j ⎭

(15)

{R q ) ( e ) =

28


4. ANALISIS NUMERIK

Prosedur penyelesaian sistem persamaan adalah dengan menurunkan rumus rekrusif yang mengaitkan nilai-nilai {T} pada suatu saat t dengan nilai-nilai {T} pada waktu t + ∆t berikutnya, dimana ∆t adalah selang waktu untuk setiap langkah. Dengan rumus rekrusif itu, perhitungan dimulai dari kondisi awal saat t = 0, lalu dilanjutkan langkah demi langkah sampai dicapai waktu yang diinginkan. Secara umum bisa ditulis tθ = tD + θ∆t, dimana 0 ≤ θ ≤ t, dan n = 0,1,2,...,N. Dengan tθ sebagai waktu, maka persamaan (11) bisa ditulis sebagai : [ C ]{T }θ + [ K ]{T }θ = {R ( tθ )}

(16)

digunakan pendekatan[16] :

[ T ]n +1 − {T }n ∆t

(17)

{T }θ = (1 − θ){T }n + θ{T }n +1

(18)

{R ( tθ )} = (1 − θ ){R}n + θ{R}n +1

(19)

{T }θ =

Karena daya laser yang digunakan selalu konstan untuk satu tembakan, maka {R} bukanlah merupakan fungsi waktu, sehingga {R}n= {R}n+1. Dengan algoritma CrankNicolson diperoleh : 1 ⎡ ⎤ ⎢[ K ] + ∆t [C ]⎥{T }n+1 ⎣ ⎦

2 ⎡ ⎤ =⎢− [ K ] + [C ]⎥{T }n + 2{R}n ∆t ⎦ ⎣

(20)

Untuk setiap langkah pertambahan waktu, nilai {T}n+1 bisa dihitung sebab nilai-nilai matrik [K], [C], {R}n, dan {T}n diketahui sebelumnya, perhatikan bahwa perlu dilakukan pemasukan syarat-syarat batas pada setiap langkah iterasi. Untuk melihat konvergensi hasil komputasi, struktur dibagi dalam 3 macam diskritasi (dari segi jumlah elemen yang mewakilinya), yaitu : 50, 98, dan 140 elemen. Lapisan ke empat (TCO) dan ke lima (gelas) tidak disertakan dalam analisa pemanasan, karena syarat batas alamiah (T = T0) harus diberikan pada lapisan ke tiga, yaitu isolator/konduktor. Matrik-matrik [k], [c], dan {R} dari tiap-tiap elemen hingga, harus dirakit agar bisa (i ) dari elemen ke-i, pada mewakili keseluruhan struktur. Skema penyimpanan elemen matrik k pq

matrik kekakuan [K] dengan nomor baris rp dan nomor kolom cq, adalah :


k (pqi) ⎯ ⎯→ Krp cq , untuk p,q = 1,2,3.

29

(21)

dengan rj = cj = Nij Nij adalah nomor simpul ke-j dari elemen ke-i. Matrik kekakuan yang diperoleh adalah sebuah matrik pita (band matrix) yang bisa disimpan dalam modus penyimpanan pita dengan bantuan SUBROUTINE CRGRB dari perangkat lunak library IMSL. Penyelesaian persamaan matrik yang disimpan dalam modus pita ini, dilakukan dengan bantuan SUBROUTINE LSARB dari perangkat lunak library IMSL, yang menggunakan prosedur eliminasi Gauss.

5. HASIL DAN DISKUSI

Grafik 1 menyatakan kondisi ambang untuk pemotongan lapisan Al dan lapisan a-Si dari struktur Al(5000Å)/a-Si(5000Å)/isolator (27µm)/TCO(5000Å)/gelas, untuk tiga macam diskritisasi. Suhu lapisan Al dan a-Si untuk diskritisasi pertama (50 elemen hingga) lebih tinggi daripada suhu lapisan Al dan a-Si untuk diskritisasi ke dua (98 elemen hingga), dan ke duanya lebih tinggi daripada suhu lapisan Al dan a-Si untuk diskritisasi ke tiga (140 elemen hingga). Terlihat bahwa konvergensi hasil komputasi masih cukup berarti, tetapi karena jumlah maksimum elemen hingga yang bisa digunakan dengan perangkat komputer 386 DX (yang tersedia di Lab. Semikonduktor Fisika Material, Jurusan Fisika FMIPA ITB) adalah 140 buah, maka analisa selanjutnya hanya menggunakan 140 buah elemen. Banyaknya elemen hingga yang dipakai untuk mewakili idealisasi struktur, ternyata sangat mempengaruhi hasil perhitungan untuk suhu simpul dari tiap elemen hingga. Untuk memperoleh hasil yang lebih bisa dipercaya, disarankan untuk melakukan perhitungan ulang dengan menggunakan komputer yang memiliki RAM yang lebih besar agar tidak menemui masalah konvergensi. Cara lain adalah dengan memperbesar jumlah simpul untuk idealisasi struktur pada saat sinar laser, karena sinar laser yang diradiasikan dianggap berdistribusi Gauss. Untuk perbandingan, dilakukan juga perhitungan dengan menggunakan algoritma Galerkin (θ=23). Hasil komputasi dengan algoritma Crank-Nicolson (θ=½) ternyata tidak

30


berbeda dengan hasil komputasi dengan algoritma Galerkin. Jadi tidak ada alasan khusus dalam memilih algoritma dalam penulisan program MEH dalam persoalan ini. Gambar 2 menunjukkan distribusi ruang pada pusat sinar setelah struktur Al(5000Å)/aSi(5000Å)/isolator(27µm)/TCO(2000Å)/gelas disinari selama 145 ns untuk memotong lapisan Al dan a-Si (laser scribing), agar Al dan a-Si terpisah secara listrik. Bila digunakan daya laser yang lebih besar daripada 4,8 x 107W/cm2, Al dan a-Si akan meleleh, tapi TCO belum dipengaruhi karena konduktivitas panas isolator sangat kecil, yaitu 0,0144 J/cm.K.S. Kecilnya nila konduktivitas panas isolator menyebabkan rendahnya perubahan suhu pada permukaan isolator, hal ini bahkan menyebabkan adanya anomali hasil perhitungan suhu, yaitu ada beberapa simpul yang suhunya di bawah suhu awal. Untuk memperoleh lebih banyak data suhu simpul data pada lapisan ke tiga (isolator/konduktor), dilakukan beberapa kali perhitungan dengan mengurangi tebal lapisan ke tiga itu dalam kondisi yang masih memungkinkan. Hal ini bisa dilakukan karena koordinat-koordinat simpul dibuat berkaitan dengan tebal dan jari-jari lapisan. Gambar 3 menunjukkan distribusi suhu ruang pada pusat sinar untuk struktur Al(5000Å)/a-Si(5000Å)/konduktor/(12µm)/TCO(2000Å)/gelas yang disinari selama 20 ns untuk mengkonduksikan Al dengan konduktor (laser welding). Di atas daya laser 1,5 x 108W/cm2, Al, a-Si, dan konduktor akan meleleh. Karena konduktivitas panas konduktor jauh lebih besar daripada konduktivitas panas isolator, yaitu 3,222 J/cm.K.S., maka bagian dalam dari konduktor ikut meleleh. Distribusi suhu ruang yang diperoleh oleh Y.Kishi dkk[9] diberikan dalam grafik 4 untuk laser scribing dan grafik 5 untuk laser welding. Hasil yang diperoleh dalam pekerjaan ini, ternyata sedikit berbeda dengan yang diperoleh oleh Y. Kishi dkk. (lihat grafik 4 dan 5). Perbedaan ini disebabkan oleh beberapa hal, antara lain : (a) Perbedaan pemilihan banyaknya elemen hingga yang digunakan dalam diskritisasi struktur. (b) Perbedaan perhitungan koefisien refleksi energi sinar laser. (c) Perbedaan perkiraan besarnya jari-jari berkas laser


31

6. KESIMPULAN

Hasil perhitungan menunjukkan bahwa kondisi optimum untuk laser scribing adalah penyinaran selama 145 ns dengan rapat daya laser di atas 4,8 x 107 W/cm2. Dan kondisi optimum untuk laser welding adalah penyinaran selama 120 ns dengan rapat daya laser di atas 1,5 x 108 W/cm2. Hasil ini sesuai dengan yang diperoleh oleh Y. Kishi dkk[19].

Ucapan Terima Kasih Penulis mengucapkan terima kasih kepada Departemen Pendidikan dan Kebudayaan yang telah memberikan beasiswa TMPD selama studi S2 di Jurusan Fisika, FMIPA ITB.

DAFTAR PUSTAKA

[1] Baker, A.J., Finite Element Computational FLuid Mechanics, Mc. Graw-Hill, Singapore, 1985. [2] Cook, R.D., Concepts and Applications of Finite Element Analysis, John Wiley & Sons, New York, 1981. [3] Desai, C.S. Terjemahan Wirjosoedirdjo S., Dasar-dasar Metode Elemen Hingga, John Wiley & Sons, New York, 1982. [4] Hamakawa, Y., Amorphous Semiconductor Technologies and Devices, Japan Annual Reviews in Electronics, Computer and Telecomunications (1982), Ohsha Ltd., Tokyo. [5] Heavens, O.S., Optical Properties of Thin Solid Films, Butterworts Scientific Pub. Washington DC, 1955 [6] Huebner, K.H. dan Thornton, E.A., The Finite Element Method for Engineers, John Wiley & Sons, New York, 1982. [7] Incropera, F.P. dan De Witt, D.P., Introduction to Heat Transfer, John WIley & Sons, Canada, 1985.

32


[8] Kant, T., Finite Elements in Computational Mechanics, (Vol. I dan II). Proceedings of the International Conference Held at The Indian Institute of Technology Bombay; 2-6 Dec 85, Pergamon Press, Oxford, 1985. [9] Lewis, R.W., Morgan, K. dan Schrefler, B.A., Numerical Method in Heat Transfer, Vol II, John Wiley & Sons, New York 1983. [10] Microsoft Fortran Ver. 5.1. Quick Reference Guide, Reference, Uiser's Guide (Microsoft Corporation, USA, 1991). [11] Nakamura, S., Computational Methods in Engineering and Science, John Wiley & Sons, New York, 1977. [12] Nakano, S., Matsuoka, T., Kiyama, S., Kawata, H, Nakamura, N., Nakashima, Y., Tsuda, S., Nishiwaki, H., Ohnishi, M., Nagaoka, I., and Kuwano, Y., Jpn.J.Appl.Phys. 25 (1986 1936. [13] Nasution, A., Fortran 77, Penerbit Erlangga, Jakarta, 1990. [14] Ozisik, M.N., Heat Conduction, John Wiley & Sons, Canada, 1980. [15] Pao, Y.C., A First Cource in Finite Elemen Analysis, Allyn and Bacon Inc., Massachusets, 1986. [16] Shur, M., Physics of Semiconductor Devices, Prentice Hall International, New Jersey, 1990. [17] Takahashi, K., dan Konagai, M., Amorphous Silicon Solar Cells, North Oxford Academic Pub., London, 1986. [18] Tsuda, S., Takahama, T., Isomura, M., Tarui, H., Nakashima, Y., Hishikawa, Y., Nakamura, N., Matsuoka, T., Nishiwaki, H., Nakano, S., Ohnishi, M., and Kuwano, Y., Jpn.J.Appl.Phys. 30 (1991) 8.


33

Grafik 1. Distribusi suhu ruang pada pusat sinar dari struktur Al(5000Å)/aSi(5000Å)/isolator(27µm)/TCO(2000Å)/gelas untuk 50 elemen, 98 elemen, dan 140 elemen; yang disinari laser dengan rapat daya 4,8 x 107 W/cm2

34


Grafik 2. Distribusi suhu ruang pada t=145 ns untuk struktur Al(5000Å)/aSi(5000Å)/isolator(27µm)/TCO(2000Å)/gelas pada pusat sinar.


Grafik 3. Distribusi suhu ruang pada t=120 ns untuk struktur Al(5000Å)/aSi(5000Å)/isolator(12µm)/TCO(2000Å)/gelas pada pusat sinar.

35

36


Grafik 4. Perbandingan distribusi suhu ruang dalam laser scribing berdasarkan perhitungan penulis dan perhitungan Y. Kishi dkk.


Grafik 5. Perbandingan distribusi suhu ruang dalam laser welding berdasarkan perhitungan penulis dan perhitungan Y. Kishi dkk.

37

R.A. Sani *), M. Barmawi, T.Winata, Sukirno Lab. Semikonduktor, KBK Fisika Material Jurusan Fisika FMIPA-ITB

Recommend Documents