## PEMODELAN DAN SIMULASI DISTRIBUSI MUATAN PEMBAWA IONIK PADA ZAT PADAT ANISOTROPIC FIBER DUA DIMENSI

## MODELLING AND SIMULATION OF IONIC CHARGE CARRIER DISTRIBUTION IN TWO DIMENSIONAL ANISOTROPIC FIBER SOLID STATE MATERIAL

## Susanto Sigit Rahardi
Balai Besar Bahan dan Barang Teknik, Jl. Sangkuriang No. 14 Bandung
Email: susantosr@yahoo.com Diterima: 22 September 2015 Direvisi: 24 November 2015 Disetujui: 22 Desember 2015

## ABSTRAK
Salah satu aspek penting dalam rekayasa material adalah pengaruh morfologi terhadap berbagai sifat bahan. Dalam teknologi baterai lithium-ion, morfologi bahan aktif pada elektroda mempengaruhi perambatan ion-ion pembawa muatan di dalamnya. Pada tulisan ini, ditinjau zat padat yang di dalamnya tersusun dari granul-granul anisotropic fiber sebagai tempat perambatan ion-ion  pembawa muatan. Mekanisme perambatan ion pembawa muatan pada granul-granul tersebut adalah proses difusi dan medan listrik luar. Telah dibuat model matematis dan simulasi untuk mendeskripsikan pengaruh keacakan granul-granul fiber terhadap perubahan distribusi ion pembawa muatan. Hasil simulasi mengindikasikan bahwa berlikunya (tortuosity) jalur migrasi ion-ion   mengurangi perambatan ion-ion pembawa muatan.
Kata kunci: distribusi muatan, difusi ionik, migrasi ionik, keacakan, tortuosity

## ABSTRACT
One of important aspects in material engineering is the effect of morphology on the properties of material. In lithium-ion batteries technology, morphology of active material at electrode influences the propagation of ionic charge carrier inside the batteries.   In this paper, the properties of a solid state material consisting of an anisotropic fiber granules where ions propagate were studied. The propagation mechanism is difusion process and affected by an external electric field on the material. Mathematical modelling and simulation was made to describe the influence of the randomness of anisotropic fiber granules on the charge distribution. The result of simulation indicates that the tortuosity decreases ionic migration.
Keywords: charge distribution, ionic difusion, ionic migration, randomness, tortuosity.

## PENDAHULUAN
Perilaku pembawa muatan listrik yaitu ion yang menjalar pada insterstisi kisi merupakan bahasan yang menarik. Hal ini mengkaitkan antara sifat listrik dengan mikrostrukturnya. Berdasarkan eksperimen yang dilakukan oleh Schroeder  et  al.  didapati bahwa untuk bahan dengan granul yang menyambung menunjukkan konduktivitas ion yang lebih tinggi dari pada bahan yang sama dengan granul berbentuk bulat [1].
Pengembangan selanjutnya adalah pada eksperimen oleh Lee et al.  yang mensintesis bahan lithium-ion dengan bentuk nanorod. Bahan ini menunjukkan sifat listrik yang lebih baik dari
pada berbentuk bulat dan meningkatkan unjuk kerja batere lithium-ion [2].
Untuk itu, perlu dibahas dan dimodelkan pengaruh morfologi yang acak (terdapat tortuosity ) dan tak acak terhadap sifat listrik suatu bahan anisotropic fiber dengan pembawa muatan ion.
Dengan demikian, terdapat potensi untuk dapat meningkatkan sifat listrik dengan memanipulasi mikrostruktur tanpa mengubah komposisi bahan. Hubungan ini dapat meng- kaitkan sifat kelistrikan dengan morfologi bahan tempat mengalirnya ion-ion listrik, bagaimana pengaruh bentuk granul anisotropic fiber dengan tingkat keacakan orientasi, terhadap distribusi muatan pembawa ionik pada bahan isolator
dengan pembawa muatan ionik dengan fenomena difusi interstisial ketika diberi muatan eksternal.

## METODE
Model dikembangkan dari persoalan difusi materi zat padat yang disebabkan oleh perbedaan konsentrasi. Perbedaan konsentrasi ini menjadi driving force/activation difusi materi. Faktor lain yang perlu ditinjau adalah ada medan listrik dari luar yang membuat ion-ion pada daerah interstisial menjadi terpengaruh untuk bergerak sehingga migrasi terjadi. Difusi berkaitan dengan efek thermal sedangkan migrasi berkaitan dengan perbedaan potensial. Penggabungan kedua faktor ini mempengaruhi perubahan distribusi muatan ionik pada struktur anisotropik fiber. Faktor penentu geometrik adalah tingkat keacakan yang memberikan berbagai bias jalur perambatan ion- ion (tortuosity) [3]. Bentuk geometri kompleks ini lebih mudah diimplementasikan pada simulasi komputer.
Simulasi dirumuskan dari model fisis yang telah didapatkan menggunakan metode beda hingga. Untuk setiap node merepresentasikan granul anisotropic fiber dengan tingkat keacakan tertentu. Dengan menggunakan perbedaan nilai konsentrasi di tiap node akan didapatkan nilai konsentrasi pada node berikutnya berdasarkan persamaan baru yang merepresentasikan dinamika sistem fisis.
Implementasi algoritma beda hingga diterapkan sepenuhnya pada sistem MATLAB ® dengan kode-kode yang dibuat berdasarkan persamaan baru akibat diskretisasi. Penulisan kode program dilakukan manual berdasarkan algoritma tanpa solver build-in. Simulasi dimulai dengan bagaimana keadaan representasi sistem fisis tanpa beda potensial eksternal, dan kemudian keadaan dengan adanya beda potensial eksternal. Setiap keadaan tersebut, disimulasikan berbagai urutan waktu dan perbedaan keacakan anisotropic fiber.

## Representasi Fisis
Perambatan materi ionik pada zat disebabkan karena adanya driving force perbedaan konsentrasi dan migrasi yang disebabkan oleh perbedaan potensial eksternal dipahami sebagai energy dissipator . Secara makroskopik disipator ini merupakan tahanan internal. Secara mikroskopik, disipator ini dapat dipandang sebagai mekanisme tumbukan dan getaran akibat thermal yang menghambat
pergerakan spesies pembawa muatan ini. Untuk itu, terdapat faktor difusivitas yang menentukan besar difusi. Faktor ini dipengaruhi oleh suhu, jenis pembawa muatan, dan struktur zat padat.
Adanya partikel penghambat yang menyebabkan tumbukan dan getaran menyebab- kan ion pembawa muatan menjalani jalur yang tidak lurus langsung ke elektroda tujuan. Perbedaan jalur berliku dari jalur lurus ini disebut sebagai tortuosity [4,5]. Semakin besar tortuosity semakin besar hambatan perambatan ion.
Pada penelitian ini, dicermati efek struktural (keacakan dan  tortuosity) terhadap perubahan konsentrasi pembawa muatan. Ion-ion pembawa muatan lebih mudah berdifusi dan bermigrasi pada permukaan sepanjang arah geraknya  dibandingkan  pada bagian dalam zat padat. Hal ini dapat dikarenakan lebih banyak interaksi tumbukan dan getaran dengan atom- atom “tetangga” pada bagian dalam dari pada permukaan. Namun demikian, adanya potensial Schottky yang muncul pada grain boundary atau antar permukaan menyebabkan difusi dan migrasi pada permukaan yang tegak lurus dengan arah geraknya jauh lebih sulit daripada pada bagian dalam zat padat. Pada kasus ini, potensial Schottky menjadi penghalang pergerakan ionik.
Untuk dapat memanfaatkan peran permukaan untuk meningkatkan difusi dan migrasi ionik, diadakan suatu granul dengan bentuk anisotropic  fiber . Untuk permukaan anisotropic fiber yang searah dengan pergerakan ion akan membantu pergerakan ion lebih cepat. Namun, permukaan yang tegak lurus dengan arah gerak ion cenderung menghambat karena potensial Schottky. Walaupun demikian, secara keseluruhan bagian anisotropic fiber lebih banyak yang searah dengan arah gerak ion. Kondisi acak ini merupakan faktor struktur pada zat padat ini.
Pada bahan isolator, perpindahan suatu jenis material di dalam agregat merupakan peristiwa difusi. Untuk kasus ini, dikembangkan sistem kisi dengan ion bergerak pada ruang antara kisi, yaitu interstisial. Gerakan ion ini adalah acak dalam pengaruh termal. Selanjutnya, terdapat medan listrik eksternal yang di- aplikasikan pada bahan itu [6].
Terdapat perbedaan mengenai bagaimana ion itu merambat. Di dalam bahan, ion merambat atau berdifusi pada ruang interstisial, sedangkan pada permukaan bahan, ion dapat merambat dengan lebih baik. Hal ini berkaitan dengan tetangga ion yang ada pada satu sisi. Sedangkan
di dalam  kisi, ion-ion tetangga berada di semua sisi. Hal ini juga diindikasikan oleh energi aktivasi yang berbeda untuk melompat dari ruang satu ke ruang yang lain. Selain itu, fitur yang lain adalah bahwa untuk perambatan ion pada permukaan, gerakan translasi ion terkendala (constraint) pada permukaan.
Pada representasi fisis ini, dipilih sistem kontinu makroskopik dengan determinan adalah struktur anisotropik fiber. Faktor  ini merupakan nilai random tingkat keacakan anisotropic fiber (1,0-0,0). Keacakan ini merupakan nilai yang mengikuti distribusi Binomial dengan tingkat 1,0 merupakan anisotropik fiber- anisotropic fiber yang tersusun sangat rapi sampai nilai 0,0 dengan anisotropic fiber yang acak. Dengan susunan anisotropic fiber ini, diterapkan mekanisme difusi dan migrasi sepanjang anisotropik fiber.
Sistem yang ditinjau adalah zat padat isolator bermorfologi granul anisotropic fiber yang mempunyai tingkat keacakan dimana terjadi difusi dan migrasi ionik, susunan geometrik ini menimbulkan bias jalur rambatan ion-ion yang seharusnya lurus menjadi berliku-liku (tortuosity) [7]. Elemen dari sistem ini adalah anisotropic fiber yang mempunyai struktur tertentu. Selain itu, ada beda potensial eksternal yang me- nyebabkan migrasi ionik. Besaran yang diamati pada sistem ini adalah distribusi ion pembawaa muatan yang mengalami perubahan akibat perbedaan konsentrasi dan adanya beda potensial eksternal [8].
Pada pemodelan dan simulasi fisika ini, semua parameter adalah asumsi untuk ke- pentingan visualisasi dan tujuan dari pemodelan ini untuk menganalisis rasio antara bentuk tak acak dengan bentuk yang acak.

## Model Matematik
Terdapat perbedaan mobilitas ion antara kisi dan permukaan. Peristiwa transport pada bahan adalah melalui proses difusi. Tinjau silinder identik Gambar 1 berikut:
Gambar 1. Geometri  Granul Anisotropic fiber Tidak Acak Searah dengan Medan Listrik Luar Tanpa adanya medan listrik eksternal, proses difusi terjadi disebabkan oleh perbedaan
konsentrasi materi pada zat padat. Secara empiris, difusi ini dinyatakan oleh Fick sebagai bergeraknya materi pada zat padat karena perbedaan konsentrasi:
J=-D ∇ c
Dengan J adalah flux muatan ionik (mol /cm 2 detik), c adalah konsentrasi muatan yang dapat dinyatakan tanpa satuan (dimensionless), D adalah difussivity (cm 2 /detik).
Dengan asumsi bahwa dalam sistem tidak terjadi reaksi kimia,  adanya beda potensial listrik eksternal mengakibatkan fenomena lain lagi, yaitu migrasi ionik. Hal ini secara makroskopik dapat dinyatakan dengan hubungan berikut:
𝐽 = 𝐷𝑐𝑧𝐹 𝑅𝑇 ∇𝜑
Dengan J adalah flux ionik (mol/cm 2 detik),  c adalah konsentrasi muatan, yang dapat dinyatakan tanpa satuan (dimensionless), D adalah difussivity (cm 2 /detik),  z adalah valensi ion yang bermigrasi, untuk ion lithium z = +1, F adalah konstanta Faraday 96487 C/mol, R adalah konstanta gas universal 8,314 J/mol K, T adalah suhu dalam Kelvin, suhu ruang adalah 300 K, φ adalah beda potensial listrik eksternal dalam volt, dengan nilai negatif karena muatan yang ditinjau adalah ion positif yang bergerak dari tegangan lebih tinggi ke tegangan lebih rendah.
Dengan demikian, gabungan antara difusi dan migrasi menghasilkan flux ionik:

## 𝐽 = 𝐷 ∇𝑐 − 𝐷𝑐𝑧𝐹 𝑅𝑇 ∇𝜑
Transport massa menetapkan untuk asumsi aliran yang tak  terkompresi, terdapat ketentuan mengenai jumlah massa input maupun output yang berkaitan dengan kontinuitas massa. Hal ini dinyatakan dengan:
𝜕𝑐 𝜕𝑡 + ∇ . 𝐽 = 0
Yang menyatakan bahwa perbedaan massa yang terdapat dalam sistem sama dengan selisih massa yang masuk ke dalam sistem dan yang keluar, atau divergensi dari flux massa.
Dinamika konsentrasi ion pada zat padat ini dapat dinyatakan sebagai berikut:
Medan listrik E
L n silinder R = L/2n
𝜕𝑐 𝜕𝑡 + ∇ . 𝐽 = 0 𝜕𝑐 𝜕𝑡 + ∇ . ( −𝐷 ∇𝑐 − 𝐷𝑐𝑧𝐹 𝑅𝑇 ∇𝜑 ) = 0
𝜕𝑐 𝜕𝑡 = ∇ . �𝐷 ∇𝑐 + 𝐷𝑐𝑧𝐹 𝑅𝑇 ∇𝜑� (1)
Persamaan terakhir ini merupakan representasi sistem fisis yang ditinjau.
Persamaan  (1)  merupakan persamaan diferensial parsial orde dua dan orde satu, dan juga persamaan diferensial parsial temporal. Beda potensial merupakan nilai yang tetap karena hanya bergantung terhadap nilai ekstern sehingga tidak masuk dalam operasi pada operator diferensial parsial. Suku pertama merupakan diferensial parsial temporal konsentrasi. Suku ke dua, setelah tanda sama dengan, adalah diferensial parsial spasial orde dua. Suku ketiga adalah diferensial parsial spasial orde pertama. Persamaan inilah yang akan dibuat persamaan baru dalam rangka diskretisasi dalam metode beda hingga.
Difussivity D (cm 2 /detik) muncul dua kali dalam persamaan tersebut dan dalam kaitan dengan sistem fisis yang kita tinjau ini merupakan faktor penentu karena efek struktur terdapat pada faktor ini. Berdasarkan referensi [9,10], difusivity bulk material lithium manganat berkisar antara 10 -9 sampai 10 -11 . Untuk itu, dalam pembahasan ini, digunakan difusivity bulk 10 -9 cm 2 /detik. Pada referensi yang lain, difusivity bisa sampai 8 x 10 -5 . Diasumsikan bahwa difusivity untuk permukaan lebih besar dibandingkan dengan bagian dalam (bulk) sampai dua kali lipatnya. Untuk itu, digunakan difusivity efektif sebagai berikut:
Gambar 2. Diskretisasi Material dengan
Geometri Anisotropik Fiber
Nilai node u i,j adalah nilai konsentrasi ion. Sumbu mendatar adalah kedudukan spasial granul (node spasial), sedangkan sumbu vertikal adalah diskretisasi waktu evolusi untuk simulasi (node temporal).
Pada  Gambar 2, diilustrasikan terdapat anisotropic fiber dalam arah sumbu x. Anisotropic fiber ini didiskritkan oleh node yang menggambarkan nilai konsentrasi ion (c) pada posisi itu, dinyatakan oleh u i,j = c(i,j). Difusivity yang digunakan adalah difusivity efektif (D eff ) yang menggabungkan nilai difusivity bulk dan permukaan.
Nilai node pada waktu selanjutnya ditentukan oleh nilai node pada waktu se- belumnya, itu beda hingga ke depan (forward) Untuk itu, diperlukan suatu persamaan yang dapat menggambarkan situasi dinamika pada metode beda hingga ini berdasarkan persamaan representasi, yaitu persamaan (1). Persamaan (1) diubah menjadi persamaan dinamika yang didiskritkan sebagai berikut:
𝑐 ( 𝑖 , 𝑗 + 1) − 𝑐 ( 𝑖 , 𝑗 ) ∆𝑡 = 𝐷 � [ 𝑐 ( 𝑖 + 1, 𝑗 ) − 2 𝑐 ( 𝑖 , 𝑗 ) + 𝑐 ( 𝑖 − 1, 𝑗 )] ( ∆𝑥 ) 2 � + 𝐷𝑧𝐹 𝑅𝑇 𝛻𝜑 � [ 𝑐 ( 𝑖 + 1, 𝑗 ) − 𝑐 ( 𝑖 , 𝑗 )] ( ∆𝑥 ) � (2)
Persamaan (2) ini merupakan salah satu pendekatan numerik untuk membahas persamaan (1). Nilai pada node c(i,j+1) ditentukan dari nilai node c(i+1,j) , c(i,j) dan c(i-1,j).    Nilai dengan indeks i menyatakan node spasial, sedangkan indeks j menyatakan node  temporal. Selain itu, diskretisasi ini juga ditandai dengan nilai ∆𝑡 dan ∆𝑥 yang dipilih dengan kisaran nilai tertentu guna stabilitas nilai hasil komputasi .
Asumsi penting lain yang ditetapkan adalah adanya tumpang tindih jalur tidak serta merta menambah konsentrasi ion pada node bersangkutan. Hal ini diambil berdasarkan ketentuan bahwa ion-ion merambat pada jalur dengan mekanisme antrian satu-satu pada daerah dua dimensi.

## HASIL DAN PEMBAHASAN
Kode dan Hasil Simulasi dengan MATLAB ®
Dengan menggunakan persamaan (2), dibuat algoritma pemrograman pada MATLAB ® . Langkah pertama adalah diskretisasi pada kisi- kisi dengan node yang ditentukan nilai awal yang tetap untuk menggambarkan bahwa node-node itu mewakili anisotropic fiber . Dalam MATLAB ® , dibuat fungsi untuk membangkitkan profil pertama dengan kode seperti berikut:
Fungsi  funcProfile (Lx,Ly,n,m,p) mem- punyai masukan berupa panjang material Lx dan Ly, banyak partisi pada x (n), banyak partisi pada y (m) dan banyak time sequence (p). Fungsi ini akan dipanggil pada  program  utama. Sebagai contoh:
Gambar 3. Diskretisasi Material dengan Geometri Anisotropic fiber pada Matriks yang Dibangkitkan oleh Fungsi funcProfile
Nilai node 1 berarti anisotropic fiber melewati node tersebut, sedangkan nilai 0 berarti tidak ada anisotropic fiber yang melewati node tersebut. Untuk anisotropic fiber yang  tumpang
tindih pada satu node dianggap merupakan node yang sama, yaitu bernilai 1.
Ketentuan lain pada pembangkitan profile ini adalah bahwa tedapat fungsi random dengan Distribusi Poisson. Distribusi ini dipilih karena untuk mendapatkan angka yang mengkategorikan bahwa anisotropic fiber itu lurus, belok ke kanan atau  belok ke kiri. Ketiga kategori ini membuat nilai pada profil berubah seperti pada contoh di atas. Untuk nilai keacakan 0,85 akan menghasilkan susunan anisotropic fiber yang lebih rapih bila dibandingkan dengan nilai yang lebih rendah. Untuk nilai keacakan 1,0 atau 0,0 menghasilkan susunan anisotropic fiber yang rapih dan lurus semua searah dengan arah gerak migrasi ion. Fungsi random ini mempunyai nilai setangkup antara parameter 1,0 sampai 0,5 dan 0,5 sampai 0,0. Oleh karena itu, dapat digunakan salah satu setangkup parameter tersebut. Parameter fungsi random Binomial yang menunjukan tingkat keacakan paling tinggi adalah 0,5.
Syarat awal dari komputasi ini adalah menggunakan nilai tetap untuk seluruh node. Pada node-node di pinggir mempunyai nilai misalnya  c=100,0. Sedangkan node yang lain mempunyai nilai c=0. Kondisi ini digambarkan pada Gambar 4 berikut:
Gambar 4. Nilai Awal
Nilai batas menggunakan nilai awal pada posisi pinggirnya, yaitu pada posisi node x=0 nilai c=100, dan pada sisi lainnya menggunakan nilai c=0. Nilai ini dipertahankan sepanjang kompotasi berlangsung sebagai syarat batas Dirichlet.
Untuk kondisi tanpa beda potensial listrik eksternal, didapati evolusi konsentrasi seperti pada gambar berikut:
Gambar 5. Profil Distribusi pada t = 100 dan Nilai Keacakan 0,85
0 0.5 1 0 0.5 1 1.5 2 0 20 40 60 80 100 function w = funcProfile(Lx,Ly,n,m,p) dx = Lx/n; dy = Ly/m; for pt=1:p %loop untuk sepanjang waktu for i=1:n % loop untuk x for j=1:m %loop untuk y u(i,j,p) = 0 ; %buat 0, belum ada anisotropik fiber w(i,j,p) = 0; y(j) = (j-1)*dy; end x(i)= (i-1)*dx; end end % anisotropic fiber for i=2:n-1 for j=1:m k = random( 'bino' ,2,0.85)-1; if u(i+k,j,1)>0 u(i+k,j,1) = 1; else u(i+k,j,1) = 0; % anisotropic fiber adalah 1, tidak anisotropic fiber adalah 0. end end end %MULTIPLIKASI KE SEQUENCE LAINNYA for b=1:p for i=1 : n for j=1 : m w(i,j,b) = u (i,j,1); end end end
Pada Gambar 5 pada t=100, nilai keacakan 0,85 mulai terdapat node-node dengan konsentrasi rendah akibat tidak adanya anisotropic fiber pada posisi tersebut.
Gambar 6. Profil Distribusi pada t = 100 dan Nilai Keacakan 0,10
Pada Gambar 5. Pada t = 100, nilai keacakan 0,1, terdapat lebih banyak node-node dengan konsentrasi rendah akibat tidak adanya anisotropic fiber pada posisi tersebut. Tapi bila dibandingkan dengan Gambar 5 masih mirip.
Kode program untuk simulasi difusi tanpa melibatkan beda potensial eksternal dengan hasil Gambar 5 dan Gambar 6 adalah sebagai berikut:
Untuk simulasi dengan menambahkan migrasi ionik akibat beda potensial eksternal, ditampilkan pada Gambar 7 berikut:
Gambar 7. Profil Distribusi pada t = 100 dan Nilai Keacakan 0,1
Gambar 8. Profil Distribusi pada t = 100 dan Nilai Keacakan 0,85
Gambar 7 dan Gambar 8 menunjukkan bahwa sebaran konsentrasi semakin jauh dibandingkan dengan gambar 6. Hal ini menunjukkan bahwa beda potensial eksternal mempercepat perpindahan muatan ionik.
Hal yang menarik dari gambar 7 dan Gambar 8 adalah bahwa tingkat keacakan yang jauh berbeda antara 0,85 dan 0,1 tidak ber- pengaruh signifikan terhadap kurva.
Kode program untuk simulasi ini adalah sebagai berikut:
0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 0 20 40 60 80 100 0 0.5 1 0 0.5 1 1.5 2 0 20 40 60 80 100 % persamaan ficks bentuk kedua Dg = 0.000008; % Koefisien pada persamaan diferensial Ficks Dgr = 2 * Dg; % grain boundary Deff = Dg + Dgr; %initial value of the function u(amplitide of the wave) for k =1 :tm for j=1:jmax for i=1:imax
c(i,j,k) =0;
end end end for k=1 :tm for j =1 : (imax) for i=1: (jmax) if i<2 c(i,j,k) = 100; else c(i,j,k) = 0; end x(i) = (i-1)*dx; end y(j) = (j-1)*dy; end end %Adanya profile oleh anisotropik fiber w = funcProfile(Lx,Ly,imax,jmax,tm); %implementasi numerik untuk solusi c persamaan diferensial for k=1: (tm-1) %time loop for j=2: (jmax-1) %loop untuk y for i=2: (imax-1) %loop untuk x c(i,j,k+1)= c(i,j,k)+ dt* (Deff/(dx*dx))* (c(i+1,j,k)-2*c(i,j,k)+c(i- 1,j,k)); end end end %gabungan dengan profile for i=1 : imax for j=1 : jmax for k=1 : tm % C(i,j,k) = w(i,j,k) * c(i,j,k); end end end C = w .* c;
%Gabungan Ficks dan Migrasi %Time evolution for boundary %implementasi numerik untuk solusi c persamaan diferensial for k=1: (tm-1) for j=2: (jmax-1) for i=2: (imax-1) c(i,j,k+1)= c(i,j,k)+ dt* (Deff/(dx*dx))* (c(i+1,j,k)- 2*c(i,j,k)+c(i-1,j,k)) + dt* V*(Deff/(dy))* (c(i,j+1,k)-c(i,j,k)); end
end end

## Perbandingan Profil Distribusi dengan Tingkat Keacakan
Tingkat konsentrasi ion yang ditunjukkan pada profil distribusi ion dapat dikuantisasi menggunakan fungsi EGP, yaitu fungsi yang merupakan integrasi antara konsentrasi ion terhadap permukaan dua dimensi. Tingkat keacakan adalah parameter fungsi random Pioson. Sedangkan persentase beda EGP adalah modulus persentase perbedaan EGP terhadap tingkat keacakan 1.
Perbandingan profil distribusi ion pada permukaan dua dimensi terhadap tingkat keacakan dapat ditabulasi pada Tabel 1 berikut:
Tabel 1. Perbandingan Tingkat Konsentrasi Ion terhadap Tingkat Keacakan R (tingkat keacakan) EGP Persentase beda EGP (%) 1 72,8039 0 0,8 60,457 16,96 0,6 56,9443 21,78 0,55 55,9467 23,15 0,555 58,5588 19,57 0,5 57,2212 21,40 0,455 58,0166 20,31 0,45 57,7089 20,73 0,4 56,8943 21,85 0,3 57,4617 21,07 0,2 61,4123 15,65 0 72,5328 0,372
Fungsi random Binomial adalah fungsi generik yang disediakan oleh MATLAB ® , dengan ciri adalah fungsi simetris. Hal ini ditandai dengan setangkupnya nilai EGP antara tingkat keacakan 1 sampai 0,5 dan tingkat keacakan 0,5 sampai 0.
Pada  Tabel 1, pada tingkat keacakan 0,55 didapati nilai persentase beda EGP sekitar 23% merupakan nilai tertinggi. Pergerakan ionik yang tidak lurus menandakan jalur gerakan yang berliku-liku, atau tortuosity. Semakin tinggi tingkat keacakannya maka tortuosity  akan meningkat sehingga mempengaruhi difusi ionik, yang dinyatakan oleh nilai persentase beda EGP.

## KESIMPULAN
Nilai presentase beda fungsi EGP 21% dari struktur acak dibandingkan dengan struktur lurus menggambarkan bahwa tingkat  keacakan anisotropic fiber meningkatkan  tortuosity  bahan yang berakibat secara signifikan mempengaruhi distribusi konsentrasi. Terdapat hubungan konsep tortuosity bahan dengan proses migrasi pembawa muatan ionik,  bahwa  semakin tinggi tingkat keacakan akan meningkatkan tortuosity ionik dan menghambat proses migrasi ionik.

## DAFTAR PUSTAKA
[ 1]   M. Schroeder, S. Glatthaar, and J. R.
Binder, 2011, “Influence of Spray
Granulation on the Properties of Wet Chemically  Synthesized Li 1.3 Ti 1.7 Al 0.3 (PO4)3 (LTAP) Powders,” Solid State Ionics , vol. 201 , no. 1, pp. 49–53.
[2] H. Lee, P. Muralidharan, C. M. Mari, R. Ruffo, and D. Kyung, 2011, “Nanorods As a High Power Positive Electrode for Rechargeable  Li-Ion Batteries”, J. Power Sources , vol. 196 , no. 24, pp. 10712– 10716.
[3] Thorat et.al, 2009 ”Quantifying Tortuosity In Porous Li-Ion Battery Materials”, Journal of Power Sources, 188 , 592.
[4] Martin et.al, 2013 “ Tortuosity Anisotropy In Lithium-Ion Battery Electrodes”, Advanced Energy Materials, DOI: 10. 1002/aenm.201301278
[5] Hutzenlaub et.al, “FIB/SEM-Based Calcu- lation of Tortuosity  In A Porous LiCoO 2 Cathode  For a Li-Ion Battery”,  Electro- chemistry Communication 27 , 77-80, 20
[6] Xioojian et.al,  2013,  “ Direct Observation of  Lithium-ion  Transport  Under an Electrical  Field in Li x CoO 2 Nanograins”,
Scientific Reports DOI: 10.1038/srep 01084.
[7] Bae, Erdonmez, Halloran, Chiang,  2013 “Design  of Battery Electrodes with Dual- Scale Porosity to Minimize Tortuosity and Maximize Performance”,  Adv.Mater, 25 , 1254-1258.
[8] Vijayaraghavan et.al, 2012 “An Analytical Method to Determine  Tortuosity  in Rechargeable  Battery  Electrode”, Journal of  the Electrochemical Society, 159 (5) A548-A552.
[9] Myonggu Park et.al,  2010,  "A Review of Conduction Phenomena in Li-Ion Batteries",  Journal of Power Sources, pp26.
[10]  Joachim Maier,  2014,  ”Pushing  to the Limit: Charge Carrier Chemistry in Extremely Small System”, Chemistry Material, 26 , 348-360.