Ważna własnośd membran elektrochemicznych polega na tym, że wywołują one powstawanie różnicy potencjałów elektrycznych między dwoma roztworami, zwaną potencjałem membranowym

-1-

Projekt 2: Modelowanie sygnału czujnika elektrochemicznego

Projekt dotyczy wykorzystanie modelu matematycznego i oprogramowania bazującego na tym modelu do symulowania sygnału czujnika elektrochemicznego. Podstawowymi elementami jakiegokolwiek czujnika chemicznego są (i) chemicznie selektywna warstwa receptorowa, (ii) element przetwornikowy. Projekt dotyczy modelowania chemicznej warstwy receptorowej, która w przypadku tzw. elektrod jonoselektywnych jest zbudowana ze specjalnej membrany półprzepuszczalnej.

Ważna własnośd membran elektrochemicznych polega na tym, że wywołują one powstawanie różnicy potencjałów elektrycznych między dwoma roztworami, zwaną potencjałem membranowym

Mem R L

,

  

   gdzie  

,

to potencjały roztworów odpowiednio po prawej i lewej stronie membrany, którą można zmierzyd. Zależnie od struktury rozróżnia się homogeniczne i heterogeniczne membrany elektrochemiczne. W przypadku membran homogenicznych cała membrana uczestniczy w tworzeniu potencjału, natomiast w membranach heterogenicznych zasadnicza częśd procesu potencjałotwórczego ma miejsce na brzegu (np. membrana szklana używana w klasycznym pH-metrze do oznaczania stężenia jonów wodorowych z wysoką selektywnością i dokładnością).

Rys. 1. Schemat modelu membrany znajdującej się pomiędzy roztworami stężeniach i-tego składnika odpowiednio

c

c

Roztwory mogą albo stykad się bezpośrednio (wtedy konieczne jest stosowanie odpowiedniej matrycy uniemożliwiającej mieszanie się roztworów), albo mogą byd rozdzielone przegrodą (membraną) o dostatecznie dużych porach, aby przynajmniej niektóre składniki mogły wędrowad do lub z membrany. Przegrodą taka może byd np. spiek szklisty lub płytka ceramiczna, które są przepuszczalne dla wszystkich składników i tylko mechanicznie blokują mieszanie się składników.

Może to byd też membrana z PCW, w której są rozpuszczone różne składniki takie jak walinomycyna (przykład tzw. jonoforu), plastyfikator i sól lipofilowa. Taka membrana określana jest mianem półciekłej.

Lewy roztwór Membrana

jonoselektywna

Prawy roztwór

c

c

= c

i

(x,t)

 =  (x,t)

c

L M R

0 d

Krzysztof Szyszkiewicz, AGH

-2-

Matematyczny opis modelu

Rozkłady stężeo i potencjału elektrycznego w membranie można uzyskad opierając się na odpowiednich prawach fizycznych. Wykorzystanie tych praw pozwala na sformułowanie właściwych równao obowiązujących w membranie oraz ograniczeo na brzegach układu (warunki brzegowe).

Ponieważ symulacje będą miały charakter dynamiczny (zależny od czasu), więc należy jeszcze pamiętad o warunkach początkowych (stan układu w chwili t=0).

Niewiadome. Model pozwala obliczyd rozkłady stężenia składników w membranie ( , )

i i

c  c x t dla x  [0, ], d t  0 i  1 ,..., r

oraz rozkład natężenia pola elektrycznego oraz potencjału elektrycznego w membranie E  E( , ), x t  ( , ) x t dla x  [0, ], d t  0

jako funkcje położenia i czasu. Znając natężenie pola elektrycznego łatwo uzyskujemy całkowity spadek potencjału wewnątrz membrany jako funkcję czasu

0

( ) ( , ) .

d

M

t E x t dx



    ⁽¹⁾

Prawa fizyczne

Podstawą do opisu rozważanego modelu jest układ równao NernstaPlanckaPoissona (NPP) opisujący ruch jonów w przestrzeni. Szczegółowo model opiera się na odpowiednich prawach fizycznych wyrażonych matematycznie jak poniżej.

(i) Prawo zachowania i-tego składnika (bez reakcji chemicznych wewnątrz membrany) w geometrii jednowymiarowej

0 ( 1, , ),

i i

c J

i r

t x

 

  

  (2)

gdzie c

(mol / m )

to stężenie, a J

(mol / (m s))

 to strumieo i-tego składnika, (ii) Prawu Gaussa (elektrodynamika) w geometrii jednowymiarowej

0

E 1

x



 

 

 (3)

gdzie    ( , ) x t to objętościowa gęstośd ładunku elektrycznego, 

 8,854 10 

N m / C 

to

przenikalnośd dielektryczna próżni (stała elektryczna), 

to względna przenikalnośd elektryczna

ośrodka (np. dla wody w temperaturze 20°C, 

 80, 2). Gęstośd ładunku  można powiązad ze

stężeniami molowymi ( ) c

oraz liczbami ładunkowymi jonów (wartościowościami) następująco

-3-

1

( , ) ( , )

r i i i

x t F z c x t





  ⁽⁴⁾

gdzie z

 ładunek i-tego jonu (liczba ładunkowa, np. dla Ca

: z   2, dla Cl

: z   1), F  stała Faradaya.

Połączenie równao (3) i (4) daje jedno równanie

0 1

E ( , ).

r i i r i

F z c x t x  

 

  ⁽⁵⁾

(iii) Wyrażenie na strumieo i-tego składnika w (strumieo Nernsta  Plancka) zakłada, że on sumą dwóch składników (dyfuzja i migracja):

i

E

i i i i i

dyfuzja migracja

c F

J D D z c

x RT

  

        ⁽⁶⁾

gdzie D

współczynnik dyfuzji i-tego jonu, z

 ładunek i-tego jonu; E E( , )  x t  natężenie pola elektrycznego; R, T, F oznaczają odpowiednio uniwersalną stałą gazową, temperaturę w skali Kelwina oraz stałą Faradaya.

Warunki początkowe i brzegowe

Warunki początkowe określają stan układu w chwili początkowej umownie przyjmowanej jako

„chwila zerowa” t  0. W modelu dopuszcza się dowolne początkowe rozkłady stężeo składników a natężenie pola elektrycznego w chwili początkowej t  0 zerowe

(0, )

( ) ( 1, , ), ( , 0) 0.

i i

c x c x i r

E x

  

 

 ⁽⁷⁾

Warunki brzegowe mogą określad wartośd niewiadomej funkcji (warunki typu Dirichleta), bądź strumienia na brzegu (warunki typu Neumana). Sposób zadania warunków brzegowych jest często podyktowany względami eksperymentalnymi.

Postuluje się następujące warunki brzegowe (tzw. warunki Changa  Jaffé’go)

0 , , ,

, , ,

( ) (0, ),

( ) ( , ),

i i Lf i L i Lb i

id i Rf i R i Rb i

J t k c k c t

J t k c k c d t

   

      

 ⁽⁸⁾

gdzie J t

( )  J

(0, ), t J

( ) t  J d t

( , ) wyrażają odpowiednio strumieo i-tego składnika dla na lewym ( x  0) lub prawym ( x  d ) brzegu, zaś

są heterogenicznymi stałymi wymiany jonowej (czasami są oznaczane tak: k

 k k ,

 k ), c

i c

są stężeniami jonów odpowiednio w po lewej i prawej stronie membrany (Rys. 1).

Ponieważ postad równania na pole elektryczne nie jest wygodna do obliczeo, w których występuje

czas, to w numerycznej implementacji problemu NPP w jednym wymiarze (1D) często wykorzystuje

Krzysztof Szyszkiewicz, AGH

-4-

się równoważne równanie zawierające tzw. prąd przesunięcia

(ang. displacement current)

  

 E /  t  :

0 0 1

1 ( ) .

r i i

r r i

E F

I t z J

t    

  

  ⁽⁹⁾

Wzór Hendersona

Przybliżoną wartośd potencjału membranowego  

w stanie stacjonarnym można obliczyd ze wzoru Hendersona, który daje dobre wartości przede wszystkim dla grubszych membran. Wzór ten wyprowadza się o taką samą jak dla modelu NPP postad strumienia, ale bez uwzględnienia zależności od czasu. Ponadto zakład się: (i) sumaryczny prąd elektryczny równy zero: 

z J

 0, (ii) rozkłady stężeo w membranie są liniowe i nie ma skoku stężenia na brzegach. W oparciu o te założenia można przeprowadzid całkowanie równania na potencjał elektryczny i otrzymad wzór:

2

, , ,

1 1

2 2

, , ,

1 1

( )

ln ,

( )

r r

i i i R i L i i i R

i i

M r r

i i i R i L i i i L

i i

z D c c z D c

RT

F z D c c z D c



 



  



 

  ⁽¹⁰⁾

gdzie wszystkie parametry mają takie samo znaczenie jak w opisanym modelu NPP.

Szczególnym przypadkiem wzoru Hendersona jest wzór Lewisa  Sargenta, który wykorzystuje się do obliczania potencjał membranowy dla elektrolitów 1:1, o takim samym stężeniu po obu stronach membrany i wspólnym kationie (lub anionie). Przykładem może byd układ KCl | NaCl, gdzie oba elektrolity mają takie samo stężenie. Dla przypadku wspólnego anionu X wzór (10) przybiera postad

ln

,

M

X M

D D

RT

F D D





   

 ⁽¹¹⁾

gdzie M

oznacza kation po lewej stronie, a N

kation po prawej stronie membrany. Zauważmy, że ze wzoru LewisaSargenta wynika, że w tym przypadku potencjał nie zależy od stężenia roztworów, oraz że potencjał ten zależy tak naprawdę od stosunku współczynników dyfuzji,

,

,

X X

D D

D D

a nie od ich wartości bezwzględnych.

Program do symulacji modelu NPP

Do symulacji potencjału oraz stężeo w membranach można wykorzystad program komputerowe npp.exe. Jest to program pracujący w tzw. trybie wierszowym, w którym w kolejnych wierszach wpisuje się polecenie (komendy) z ewentualnymi opcjami, a uruchomienie następuje przez naciśnięcie klawisza Enter. Po uruchomieniu programu bez parametrów odczytuje on dane z pliku

1 Chod powszechnie używany jest termin „prąd”, to należy pamiętad że w istocie jest to gęstośd prądu (jednostka A/m²),

-5-

input.txt , który musi znajdowad się w tym samym katalogu co sam program. Gdy program wykona obliczenia, to wyniki zapisywane są do trzech plików: res-c(x).txt, res-E(x), res-mempot(t).txt, które są tworzone przez program i umieszczane w tym samym katalogu. Ponadto program wypisuje pewne dodatkowe informacje i obliczenia na ekranie w trakcie pracy – w szczególności potencjał membranowy w kolejnych chwilach.

Przygotowanie plików danych

Plik danych (standardowa nazwa input.txt, ale można ją zmienid) jest zwykłym plikiem tekstowym ASCII zawierającym dane potrzebne do wykonania obliczeo. Plik ten jest podzielony na sekcje:

 Kolejne wiersze, zawierające nagłówki postaci:

-- Section -- opis parametrów

a po nagłówku jest linia zawierające wartości parametrów. Parametry zawsze mają wartości w podstawowych jednostkach SI.

 Parametry fizykochemiczne dla poszczególnych składników (ang. species, components) rozpoczynają się od nagłówka

-- Component --

i są podawane są w wierszach poniżej tego nagłówka wg schematu:

nazwa parametru/ów: wartość/ci

Przykładowy wygląd pliku danych wraz z komentarzami zamieszczony jest poniżej. Opisuje on układ trójskładnikowy (trzy jony): dwa kationy (I++, J+) oraz jeden anion (R-). Tak więc liczy ładunkowe wynoszą +2, +1, -1.

-- Section -- membrane thickness, d [m]

0.0001

-- Section -- r, number of components 3

-- Section -- T, temperature 298.16

-- Section -- eps_r, relative dielectric constant 10.0

-- Section -- t_end, dt_out: duration of the process, time interval for output 2000 10

-- Section -- Grid type: 0 = uniform; 1 = denser near boundary; 2 = denser in the center; 3 = both 1 and 2 1

-- Section -- n, m, a: number of nodes; number of const length intervals; thickness of boundary layer 1000 300 1.0E-5

-- Section -- rescaling factors: cs, xs, ts 1000 2.0E-9 1.0

-- Section -- Error tolerances (relative, absolute) for ODE integrator: rtol_err, atol_err 1.0E-5 1.0E-6

-- Component 1 (I++) -- D : 1.0E-11 z : +2 cL : 1.0E+3 cML cMR : 0.5 0.5 cR : 1000.0

kLf kLb : 1.0E-3 1.0E-3

Krzysztof Szyszkiewicz, AGH

-6-

-- Component 2 (J+) -- D : 1.0E-11 z : +1 cL : 150.0 cML cMR : 0.0 0.0 cR : 0.0

kLf kLb : 4.472E-7 1.0E-3 kRf kRb : 4.472E-7 1.0E-3 -- Component 3 (R-) -- D : 1.0E-11 z : -1 cL : 0.0 cML cMR : 1.0 1.0 cR : 0.0 kLf kLb : 0 0 kRf kRb : 0 0 -- end of data

Uruchamianie i działanie programu

Program wraz z przykładowym plikiem danych jest spakowany w jednym pliku npp.zip i umieszczony na stronie: http://ceramrtr.ceramika.agh.edu.pl/~szyszkin/kwim/

Po wejściu na stronę klikamy prawym przyciskiem nazwę pliku (npp.zip) i wybieramy „Zapisz element docelowy jako”. Pojawi się okienko pobierania i zapisywania pliku. Nazwy nie zmieniamy, ale możemy ewentualnie wybrad dogodną lokalizację (folder) na plik i zapamiętad ją! Po zapisaniu pliku na dysku przechodzimy do foldera (katalogu), w którym zapisaliśmy plik, klikamy go prawym przyciskiem i wybieramy „Rozpakuj” czy „Wyodrębnij wszystkie …” (trochę to zależy od wersji systemu Windows – ale generalnie chodzi o rozpakowanie archiwum ZIP). Po zatwierdzeniu, w katalogu pojawi się podkatalog npp i w nim dwa pliki: npp.exe (program) oraz input.txt (przykładowy zestaw danych).

Program można uruchomid na dwa sposoby.

I sposób (trudniejszy, ale bardziej elastyczny). Przechodzimy do wiersza poleceo systemu Windows (np. klikając przycisk Start systemu Windows i wpisując w okienku „Wyszukaj programy” polecenie cmd lub wybierając Start Wszystkie programy  Akcesoria  Wiersz polecenia). Następnie używamy polecenia cd, aby przejśd do katalogu, w którym umieściliśmy program npp.exe. Teraz wpisujemy w wierszu polecenia komendę (nazwa programu i nazwa pliku danych) i naciskamy klawisz Enter:

C:\>npp.exe input.txt <Enter>

Jeśli nie zostanie podany parametr, to program załaduje plik danych input.txt (czyli zachowa się tak, jakby został uruchomiony z parametrem input.txt).

II sposób (łatwiejszy). Używając przeglądarki folderów (np. Eksplorator Windows) otwieramy folder (katalog), w którym umieściliśmy program npp.exe i klikamy go.

Jeżeli poprawnie uruchomiliśmy program, należy jeszcze dwa razy nacisnąd Enter, a poniżej opisane kroki będą już wykonywały się automatycznie – bez udziału użytkownika.

Program wczytuje wartości z pliku danych, wyświetla je na ekranie oraz tworzy katalog o nazwie

takiej jak podany parametr. W nim zapisywane są wszystkie pliki wynikowe.

-7-

Program dokonuje przeskalowania zmiennych i wyświetla dla kontroli na ekranie wyliczone wartości (m. in. oszacowanie długości Debye’a).

Następnie generowana jest siatka dyskretyzacji wg schematu pokazanego na Rys. 2. Wartości kolejnych dyskretnych punktów x

są w plikach res-c(x).txt i res-E(x).txt.

Plik res-mempot.txt ma prostą strukturę: dwie kolumny liczb. Pierwsza zawiera czas, a druga wartośd potencjału membranowego (w woltach).

Po każdym kolejnym kroku czasowym dane dotyczące kroku dopisywane są do pliku res-mempot.txt.

Dzięki temu nawet w przypadku przerwania programu przed planowanym koocem obliczeo otrzymamy pliki wynikowe dla ostatniego czasu. Plik res-E(x).txt zawiera rozkład natężenia pola elektrycznego w membranie w chwili koocowej, czyli wartości x oraz E(x). Plik res-c(x).txt wygląda podobnie z tym, że liczba kolumn dla stężeo odpowiada liczbę składników w układzie.

Po osiągnięciu czasu t_end następuje koniec pracy integratora. Wyliczony koocowy potencjał membranowy oraz potencjał obliczony na podstawie wzoru Hendersona wypisywane są na ekranie.

Informacje o siatce (ang. grid)

Układ punktów siatki pokazany jest poniżej. Membrana jest reprezentowana przez odcinek [0, ] [ , d  x x

]. Mamy więc ogólną liczbę węzłów siatki równą n  1. Węzły są zagęszczone w pobliżu brzegu. Na odcinku o długości

jest

równoodległych interwałów wyznaczonych przez węzły x

, , x

. Odległości pomiędzy węzłami x

, x

, x

rosną wg ciągu geometrycznego o ilorazie q  1. Druga połowa jest symetryczna względem środka x



.

Do obliczenia rozkładu węzłów można przyjąd, że d  1, a potem obliczone wartości x

_przez

rzeczywistą grubośd d . Parametry użyte do zbudowania siatki są następujące: n   1 całkowita liczba węzłów (włączając węzły brzegowe: x

 0, x

 1),

długośd odcinka przy brzegu, w którym węzły są równoodległe (długośd ta jest zazwyczaj wybierana jako rzędu kilku długości stałej Debye’a, 

),

m

liczba równych odcinków na długości

q  iloraz ciągu geometrycznego dającego kolejne odległości pomiędzy węzłami h

 x

 x



q

dla k  m , ,

 1.

Wartości ilorazu q  1 ciągu obliczana jest z równania

2 2 1

1

1 ,

2

n m

a a q q q

m





       

  ⁽¹²⁾

co można przepisad następująco

2

1

1 ( 1) 1.

2

n m

q m q

a



          ⁽¹³⁾

x₀ x_n/2 xn

a

q¹ q² q³

1 2 n m

q ^{ } q⁰

xm

…

x_m+1 x_n/2-1

Rys. 2. Siatka węzłów do dyskretyzacji układu NPP. W ogólności musi byd stosowana siatka niejednorodna. W podanym przykładzie pokazana jest siatka typu 1 (wg opisu z pliku danych: Grid type) – zagęszczona w pobliżu brzegów membrany.

Krzysztof Szyszkiewicz, AGH

-8-

Zadanie 1. Symulacja kontaktowego potencjału dyfuzyjnego w zależności od stężenia.

W celu przetestowania modelu oraz oprogramowania zostaną wykonane obliczenia potencjału elektrycznego dla klasycznego układu „liquid junction” (kontakt cieczy). W ogólności potencjał taki powstaje, gdy dwa roztwory elektrolitów są w bezpośrednim kontakcie. Mogą to byd takie same lub różne elektrolity, ale ważne jest aby przynajmniej jeden jon występujący w układzie wykazywał skok stężenia na granicy styku cieczy. Ze względu na różne ruchliwości jonów prowadzi to do powstania obszaru elektrycznie naładowanego w pobliżu kontaktu, co skutkuje pojawieniem się pola elektrycznego, a w konsekwencji niewielkiego potencjału elektrycznego, który można zmierzyd.

Ponieważ istotną rolę w powstawaniu tego potencjału odgrywa zjawisko dyfuzji, to jest on także nazywany potencjałem dyfuzyjnym. Nie jest on duży, zazwyczaj od kilku do kilkudziesięciu miliwoltów. Jednak w dokładnych pomiarach potencjometrycznych może mied znaczenie. Na przykład w pomiarach aktywności jonów H

przy użyciu pehametru mamy 59 mV/pH, więc nawet niewielki potencjał dyfuzyjny 3,1 mV oznacza różnicę 0,05 pH czyli 12% [H

].

Układ kontaktujących się wodnych roztworów elektrolitów MX i NY przedstawimy tak roztwór MX(aq,c

) | roztwór NY(aq,c

)

gdzie c

i c

to stężenia roztworów po umownej lewej i prawej stronie granicy kontaktu. Można wyróżnid kilka kombinacji prowadzących do potencjału dyfuzyjnego. Najprostszym jest układ z jednym elektrolitem, ale o różnych stężeniach: MX(c

)|MX(c

) gdzie c

 c

, na przykład KCl(3,5 M)|KCl(1,0 M). Mechanizm powstawania różnicy potencjałów w takim przypadku jest pokazany na

. Jony Cl

są bardziej ruchliwe niż K

, dlatego dyfundują przez granicę z lewej na prawą stronę szybciej. Po chwili pojawia się nadmiar ładunku ujemnego po prawej stronie kontaktu i jednocześnie nadmiar ładunku dodatniego po lewej stronie kontaktu.

+ +

+

+

+ +

+ +

+ + +

+ + + + +

+

+ +

+ +

+

+ + + + +

+

+ + +

+ +

–

–

–

–

–

– –

– – –

– – – – –

– – –

–

– –

– –

–

– –

– –

– –

– –

+

+ + +

+ +

+ +

+

– –

– –

– – –

t = 0

t > 0

+ +

+

+

+ +

+ +

+ + +

+ +

+

+ +

+

+ +

+ +

+

+ + + + +

+

+ +

+ + +

–

–

–

–

–

– –

– – –

–

– – – –

– – –

–

– –

– –

–

– –

– –

– –

– –

+

+ + +

+ +

+ +

+

– –

–

–

– – –

- - - -

+ +

Rys. 3. Układ KCl(cL)/KCl(cR) gdzie cL > cR i potencjał dyfuzyjny. W chwili t=0 jest nieciągłośd stężenia na granicy styku obu roztworów. Dla t>0 następuje dyfuzyjny transport z lewej strony na prawą. Jony Cl^- są bardziej ruchliwe, więc przechodzą przez granicę szybciej niż jony K⁺. Pojawia się nadmiar ładunku ujemnego po prawej i dodatniego po lewej stronie, co prowadzi wytworzenia się różnicy potencjału elektrycznego.

-9-

C

. Należy wykonad serię symulacji potencjału dyfuzyjnego dla układu KCl(3,5 M)|KCl(x M), dla różnych wartości stężenia x=10

, 10

, 10

, 10

M i różnych szerokości układu. Wynikiem ma byd zestaw liczb z kolumny „potencjał dyfuzyjny” jak w

.

T

. Jak już wyjaśniono klasyczny dyfuzyjny potencjał kontaktowy powstaje na granicy styku dwóch elektrolitów. W abstrakcyjnym modelu matematycznym przyjmuje się, że rozwory rozciągają się od kontaktu do nieskooczoności w obie strony. Oczywiście w obliczeniach numerycznych musimy zawsze przyjąd jakiś skooczony obszar zawierający w sobie tę granicę. Jeżeli jest on dostatecznie szeroki, to dopóki efekty transportowe nie dotrą do granicy nie powinna ta szerokośd wpływad na uzyskany wynik – por. Rys. 3. Na koocach obszaru przyjmiemy ponadto, że wymiana jonowa jest praktycznie bez przeszkód – bardzo duże stałe heterogenicznej wymiany . k

W praktyce k

= 100 m/s oznacza „nieskooczonośd”. Stężenie lewego roztworu (L=left) wynosi 3.5 M, natomiast stężenie prawego roztworu (R=right) będzie się zmieniało od 1 do 0.0001 M. Tabela 1 przedstawia przykładowe wyniki potencjału dyfuzyjnego uzyskane w temperaturze 298 K, dla liczby punktów dyskretyzacji N = 1000 oraz dla membran o grubości d=200, 400, 1000, 2000 μm. Przyjęto następujące dane fizyczne:

1) współczynniki dyfuzji własnej jonów *m

/s]: D

 1.98 10 

, D

 2.01 10 ; 

2) heterogeniczne stałe szybkości wymiany jonowej [m/s]:

dla

Pozostałe parametry są zawarte w przykładowym pliku danych input.txt. Należy pamiętad o właściwych jednostkach. W

są m, mV, M. W pliku z danymi zawsze podajemy w podstawowych jednostkach SI (m, V, mol/m

).

Tabela 1. Potencjały dyfuzyjny układu KCl(3,5 M)/KCl(x M) obliczone dla różnych stężeo w prawym roztworze (R).

Potencjał powstaje w tym przypadku tylko na granicy kontaktu dwóch elektrolitów, ale do obliczeo musimy przyjąd jakąś grubośd całego układu. Kontaktowy potencjał dyfuzyjny nie powinien zależed od tej szerokości! Dla sprawdzenia wybrano trzy: 200, 400, 1000 m. Liczba punktów dyskretyzacji N=1000, temperatura 298 K.

Szerokośd obszaru d *μm+

Stężenie roztworu (L) – 3.5 M KCl Stężenie

roztworu (R) [M] Potencjał dyfuzyjny[mV]

200

10⁰ -0,2200

10^-1 -0,6242

10^-2 -1,0102

10^-3 10^-4

-1,2202 -1,2579

400

10⁰ -0,2200

10^-1 -0,6242

10^-2 -1,010

10^-3 10^-4

-1,2203 -1,2579

1000

10⁰ -0,2200

10^-1 -0,6241

10^-2 -1,0103

10^-3 10^-4

-1,2203 -1,2579

Krzysztof Szyszkiewicz, AGH

-10-

Zadanie 2. Wykonanie symulacji krzywej kalibracyjnej dla elektrody jonoselektywnej.

C

. Celem tego projektu będzie uzyskanie w oparciu o model NPP (opisującego zjawiska w membranie) wykresu o charakterze krzywej kalibracyjnej jak na Rys. 5, a następnie obliczenie współczynnika selektywności dla membrany.

T

. Podstawowym elementem elektrody jonoselektywnej jest membrana selektywna względem określonego jonu. Selektywnośd ta oznacza, że potencjał elektryczny powstający na takiej membranie jest zdeterminowany w głównej mierze przez stężenie (a dokładniej przez aktywnośd termodynamiczną) danego jonu w roztworze kontaktującym się z membraną.

Ponieważ potencjał membrany da się stosunkowo łatwo zmierzyd więc dzięki temu można zbudowad czujnik, który poprzez pomiar potencjału będzie służył do oznaczania stężenia (aktywności) wybranego jonu w badanej próbce. Idealnym czujnikiem byłby oczywiście taki, który reagowałby tylko na wybrany jon (np. Na

) w roztworze, a występowanie innych jonów nie miałoby żadnego wpływu na powstający potencjał. W praktyce jednak membrana jest także czuła na inne jony, a jej przydatności wynika z faktu, że wrażliwośd na jon główny jest co najmniej o rząd większa niż na jony przeszkadzające. Koleją ważną cechą takich czujników jest tzw. granica wykrywalności (detekcji) czyli wartośd stężenia poniżej którego czujnik nie daje sygnału pozwalającego na zmierzenie stężenia jonu.

W praktyce granica wykrywalności jest rządu 10

10

M, ale stosując odpowiednie techniki można ją obniżyd nawet do 10

M.

Elektryczny potencjał równowagowy idealnej elektrody jonoselektywnej dany jest przez równanie Nernsta

ln _i log _i,

E E RT a const S a

  nF  

(14)

gdzie E to potencjał standardowy, a

to aktywnośd oznaczanego jonu (jon główny, ang. primary ion) na powierzchni membrany jonoselektywnej, n w pewnych sytuacjach może byd liczbą elektronów uczestniczących w reakcji redoks, ln oznacza logarytm naturalny (o podstawie e=2,78…), a log to logarytm dziesiętny, który jest wygodniejsze w praktycznych obliczeniach. Jak widad ze wzoru (14) ściśle rzecz biorąc potencjał ten zależy od aktywności termodynamicznej oznaczonego jonu, a nie bezpośrednio od jego stężenia. W naszych obliczeniach posłużymy się jednak uproszczeniem i zastąpimy aktywności stężeniami. Druga forma wzoru (14) jest często stosowana, gdyż podkreśla liniowy charakter zależności mierzonego potencjału od stężenia, a parametry const i S są w praktyce wyznaczane doświadczalnie dla konkretnej elektrody. Na przykład idealna odpowiedź membrany jonoselektywnej na bazie Ag

S (membrana krystaliczna) do oznaczania kationów srebra ma postad

2.303RTlog[Ag ], E const

F

  

(15)

co w temperaturze T=25°C daje zależnośd

0,05913 log[Ag ],

E  const  

(16)

2 Sokalski et al., „Lowering the Detection Limit of Solvent Polymeric Ion-Selective Electrodes. 1. Modeling the Influence of Steady-State Ion Fluxes”, Anal. Chem. 1999, 71, 12041209.

-11-

gdzie wynik jest wielkością w woltach (V), a stężenie jonów Ag

wyrażamy w mol/dm

.

Wykres potencjału w funkcji mierzonego stężenia dla idealnej elektrody jonoselektywnej będzie miał postad linii prostej w układzie logarytmicznym. Przykład takiego wykresu jest na Rys. 4.

W praktyce wiadomo jednak, że elektrody jonoselektywne reagują na inne jony niż te dla których są dedykowane. Wpływ takich zakłócających jonów (ang. interfering ions) opisuje się pół-empirycznym równaniem Nikolskiego  Eisenmana (NE)

/ /

, ,

ln( ) log( ),

| |

i j i j

z z z z

pot pot

i i j j i i j j

i

E E RT c K c const S c K c

   z F     (17)

gdzie K

jest tzw. potencjometrycznym współczynnikiem selektywności (krótko: współczynnik selektywności; używa się też oznaczeo K

, k

), a z

_i z

to liczby ładunkowe jonów głównego i przeszkadzającego. Równanie to można traktowad jak uogólnienie równania Nernsta (14): dla

,

0

K

 równanie NE przechodzi w równanie Nernsta i elektroda jest idealnie selektywna. Zauważmy tu pewną niezręcznośd językową: im mniejsza jest wartośd współczynnika K

tym lepsza jest selektywnośd elektrody. Zatem gdy współczynnik selektywności jest duży, to selektywności elektrody jest mała. (Z punktu widzenia językowego lepiej byłoby zdefiniowad jako współczynnik selektywności odwrotnośd K

).

Rys. 4. Idealna zależnośd Nernsta odpowiedzi potencjałowej czujnika jonoselektywnego. Wykres przedstawia zależnośd E od log(ci). Zgodnie ze wzorem(14) jest to linii prosta (w układzie logarytmicznym). Różne nachylenia związane są z różnymi wartościami współczynnika S. Tendencja jest taka, że dla jonów o większych liczbach ładunkowych nachylenie będzie mniejsze.

Typowy wygląd krzywej opisującej zależnośd potencjału od logarytmu stężenia na podstawie wzoru NikolskiegoEisenmana (17) jest pokazany na Rys. 5. Widad na nim charakterystyczny fragment

„odpowiedzi nernstowskiej” (fragment prostoliniowy), który przechodzi w obszar „limitu detekcji”

(granicy wykrywalności)  fragment gdzie następuje wygięcie linii i przejście do odcinka poziomego.

Wykres taki nazywamy także krzywą kalibracyjną, gdyż w praktyce aby można było użyd elektrodę jonoselektywną, musimy wyznaczyd współczynniki E°, S, K

doświadczalnie (kalibracja). W tym celu stosuje się wzorcowe roztwory (o znanym składzie) i mierzy powstający potencjał. Tak wyznaczone wartości są punktami w układzie współrzędnych, przez które musi przechodzid krzywa z Rys. 5. W

B

A

Krzysztof Szyszkiewicz, AGH

-12-

teorii wystarczą trzy pomiary dla zwykłego równania NE, ale aby otrzymad większą dokładnośd może byd ich więcej. Numeryczne wyznaczanie E°, S i K

najlepiej przeprowadzid metoda minimalizacji.

Rys. 5. Wykres potencjału od logarytmu stężenia dla przypadku elektrody, która nie jest doskonale selektywna. Wykres ten jest otrzymany w oparciu od wzór NikolskiegoEisenmana (17):

E = E° + Slog(c_i+Kijcj), gdzie przyjęto E°=0,5 V, S=0,077 V, Kij=610^-2, cj=0,1 M.

Membrana jest opisana przez szereg parametrów: grubośd (d), przenikalnośd elektryczna ( 

), współczynniki dyfuzji składników w membranie (D

), liczby ładunkowe jonów (z

), heterogeniczne stałe wymiany jonowej ( k

, k

, k

, k

). Ponadto na działanie membrany, a w szczególności na powstawanie potencjału membranowego, wpływ ma stężenie w roztworach po obu stronach membrany. Jeden z tych roztworów jest elementem składowym czujnika jonoselektywnego (wypełnia przestrzeo pomiędzy membraną a elektrodą odniesienia znajdującą się w obudowie czujnika), dlatego też wchodzi w skład parametrów charakteryzujących membranę. Natomiast drugi z roztworów reprezentuje badana próbkę i w trakcie symulacji, aby uzyskad krzywą kalibracyjną należy wykonad kilka/kilkanaście symulacji z różnymi wartościami tego stężenia (np. w zakresie 10

10

M i zanotowad wartośd potencjału membranowego dla każdej wartości stężenia. W ten sposób otrzymamy punkty, które umieszczone w układzie Elog(c) powinny układad się na krzywej podobnej do tej z Rys. 5. W drugim kroku należy wyznaczyd parametry krzywej kalibracyjnej E°, S, K

metodą optymalizacji z wykorzystaniem arkusz kalkulacyjnego Ms Excel i dodatku Solver. Przykładowy wygląd danych w arkuszy pokazany jest na Rys. 6. W kolumnie log( c

/ M) są logarytmy stężeo, w kolumnie

(sym) / mV

E obliczone potencjały (w miliwoltach) dla tych stężeo. Na podstawie tych wartości budujemy funkcję celu zdefiniowaną jako





Gf ( , ,

)

( ) log(

) ,

k

E S K    E sym    E S c  K c ⁽¹⁸⁾

gdzie sumowanie jest po wszystkich obliczonych punktach (czyli po wierszach w kolumnie; w podanym przykładzie na Rys. 6 jest ich 11). W arkuszu zmienne E S K  , ,

są umieszczone w komórkach B1, B2, B3, a sama funkcja celu korzystająca z funkcji arkusza SUM jest w komórce E1.

Wykorzystuje ona pomocniczo zdefiniowane wartości z kolumny ( (sym) E  E (wzór_EN)) .

Teraz

można uruchomid dodatek Solver, w którym komórki zmieniane to te, gdzie jest funkcja celu (w

naszym przypadku komórka E1). Ustawiamy to w oknie dialogowym Parametry dodatku Solver.

-13-

Należy pamiętad o wyborze opcji Ustaw cel  min. Warto też wprowadzid dodatkowe ograniczenie, aby komórki zmieniane B2 i B3 były nieujemne, gdyż , S K

 0. Zauważmy, że parametr E nie może mied takiego ograniczenia, dlatego że potencjał może byd ujemny.

Rys. 6. Przykładowa organizacja danych w arkuszu kalkulacyjnym do obliczenia parametrów E°, S i Kij występujących we wzorze NE na podstawie potencjału w membranie dla różnych wartości stężenia jonu głównego, cI. Jak widad w kolumnie log(cI/M) są wartości odpowiadające stężeniom (jonu oznaczanego ) 10^-1010⁰ M, a w kolumnie E(sym)/mV) znajdują się odpowiadające tym stężeniom wartości potencjału membranowego.

D

,

. Obliczenia będą wzorowane na publikacji *4+ (w części dotyczącej ISE).

Membrana o grubości d  100  m (0,0001 m) ma byd selektywna względem jon dwudodatniego I

w obecności jedno-dodatniego jonu przeszkadzającego J

. W układzie występuje też jon ujemny (przeciw-jon) X

. W modelu NPP (2)(8) selektywnośd jest zdeterminowana przede wszystkim wartościami heterogenicznych stałych przejścia, a w mniejszym stopniu przez wartości współczynników dyfuzji, liczb ładunkowych czy stałą elektryczną membrany. Jeżeli jon „i” nie może wchodzid do membrany przez lewą granicę (Rys. 7), to k

 0. Jest możliwe, aby ten sam jon mógł opuszczad membranę – wtedy będzie k

 0. W naszym przypadku zakładamy, że jon X

nie może ani wchodzid, ani wychodzid z membrany (przez żadną granicę), więc k

 k

 0,

, ,

0.

X Rf X Rb

k

 k

 Selektywnośd I

względem J

będzie związana głównie z dużo większą wartością stałej

, , J L f

k

w stosunku do

.

I Lf

k

Szczegółowe wartości prezentuje Tabela 2, a także Rys. 7.

Tabela 2. Współczynniki dyfuzji, liczy ładunkowe i stałe przejścia roztwór/membrana.

Jon D

(m s )

z

( )  k

(ms )

k

(ms )

k

(ms )

k

(ms )

I

10

+2 10

10

10

10

J

10

+1 4,47210

10

4,47210

10

X

10

-1 0 0 0 0

Jak widad z danych w tabeli jon I

wchodzi do membrany znacznie łatwiej (10

m/s) niż jon J

(4, 472 10 m / s), 

ale oba jony z taką samą łatwością opuszczają membranę. Jon X

jest całkowicie

blokowany (nie może przechodzid przez granicę membrana/roztwór). Pamiętajmy, że membrana ma

Krzysztof Szyszkiewicz, AGH

-14-

dwie granice ( : L x  0) oraz ( : R x  d ). W modelowanym przypadku zakładamy, że oba interfejsy przepuszczają jony tak samo (symetria).

Początkowe wartości stężeo podane są w kolejnej tabeli – Tabela 3. W lewej części układu (roztwór próbkowany) stężenie jonu I

będzie zmieniane w zakresie 10

 10 mol m ,

stężenie jonu przeszkadzającego J

będzie zawsze takie samo (0,15 M), a stężenie X

można pominąd. Oczywiście w rzeczywistym układzie jakiś przeciw-jon musi występowad, ale w naszym przypadku

obliczenia w membranie „widzą” tylko jony X

, które tam były od początku, a jony X

z roztworów otaczających nie mogą wejśd ani wyjśd. Można więc przyjąd dowolną wartośd stężenia X

w lewym i prawym roztworze (w pliku przyjęto zero) – nie będzie to miało wpływu na obliczenia w tym przypadku. W samej membranie stężenia są na początku stałe i wynoszą: I

(0,0005 M), J

(0), X

(0,001M). Oczywiście stężenia w membranie muszą byd dobrane tak, aby w chwili początkowej był spełniony warunek elektroobojętności. W prawej części układu (roztwór odniesienia, który normalnie nie powinien zmieniad składu, stężenia tam są na ogół duże) mamy: I

(1 M), J

(0 – czyli brak jonu przeszkadzającego w membranie na początku), X

(nie ma znaczenia, bo k

 0).

Tabela 3. Stężenia początkowe poszczególnych jonów w różnych częściach układu. Zauważmy, że stężenia w membranie są tak dobrane, aby był spełniony warunek elektroobojętności. Natomiast stężenie przeciw-jonu X^ w roztworach otaczających można w obliczeniach przyjąd dowolne, gdyż jon ten nie przechodzi przez granice faz.

Jon c

(mol m )

c

(mol m )

c

(mol m )

I

10

10

0,5 10

J

150 0 0

X

dowolne 1 dowolne

W pliku input_ISE.txt omówione dane są już wprowadzone, a także pozostałe dane fizyko- chemiczne. Jeżeli chodzi o parametry numeryczne to liczba węzłów jest 200. Można poeksperymentowad z większą liczbą węzłów, np. 500, 1000. Czas procesu (t_end, duration of the process) jest ustawiony na 100 s, co może byd wartością zbyt krótką, gdyż potencjał w równaniu Nernsta czy Nikolskiego  Eisenmana jest potencjałem równowagowym, lub przynajmniej stacjonarnym (nie zmieniającym się w czasie). Jeżeli po 100 s nie będzie wyraźnej stabilizacji należy wydłużyd czas, np. do 400 s lub więcej.

W

. Do symulacji można użyd programu npp.exe, a dane wejściowe z przygotowanego pliku input_ISE.txt. Z opisu zawartego w sekcji Uruchamianie i działanie programu wiemy, że po uruchomieniu programu szuka on w bieżącym katalogu pliku input.txt i pobiera z niego dane. Dlatego plik z danymi powinien byd w tym samym katalogu co program

Rys. 7. Przykładowe heterogeniczne stałe przejścia w zadaniu modelującym krzywą kalibracji dla elektrody jonoselektywnej.

2

3 ,

7 ,

,

10 4, 472 10 0

I Lf

J Lf

X Lf

k k k













 



2

3 ,

3 ,

,

10 10 0

I Lb

J Lb X Lb

k k k

















2

3 ,

7 ,

,

10 4, 472 10

0

I Rf

J Rf

X Rf

k k k













 



2

3 ,

3 ,

,

10 10

0

I Rb

J Rb X Rb

k k k

















-15-

npp.exe. Ponieważ plik z danymi ma nazwę input_ISE.txt to musimy wykonad jedną z dwóch czynności.

(a) Uruchamiamy z wiersza poleceo, ale z nazwą pliku jako parametrem:

C:\>npp.exe input_ISE.txt <Enter>

Ta metoda wymaga jednak umiejętności otwarcia trybu wiersza poleceo systemu Windows.

(b) Zmieniamy nazwę pliku input_ISE.txt na input.txt i uruchamiamy program przez kliknięcie. Gdy w katalogu był wcześniej jakiś inny plik o nazwie input.txt, to będzie konflikt nazw, dlatego w takim przypadku lepiej zmienid mu wcześniej nazwę np. input1.txt.

Rys. 8. Fragment pliku danych input_ISE.txt, w którym zmieniamy stężenie jonu głównego I²⁺ dla kolejnych uruchomieo programu npp.exe. Typowe wartości to 1E-7, 1E-6, …, 1E3 (czyli 10^-7, 10^-6, …, 10³).

Gdy program zakooczy obliczenia możemy odczytad wartośd potencjału membranowego z pliku mempot.txt (ostatni wiersz, druga kolumna). Jest on też widoczny na ekranie – gdy uruchomiliśmy program wg punktu (a). W ten sposób uzyskujemy parę ( c

, E ). Wykonujemy takich obliczeo kilkanaście, np. dla stężeo 10

, 10

, …, 10

mol/m

i otrzymujemy odpowiadające potencjały E

, E

, itd. Wartości te umieszczamy w arkuszu kalkulacyjnym jak na Rys. 6 (w kolumnie log(c

) i E(sym)).

Ostatnim etapem jest zdefiniowanie funkcji celu danej wzorem (18). Można posłużyd się funkcją arkusza SUMA, przygotowując wcześniej wyrażenia do sumowania, które w przykładzie na z Rys. 6 są w kolumnie ( (sym) E  E (wzór_EN)) .

Szczegóły są po wzorze (18). Teraz uruchamiamy dodatek Solver i powinniśmy uzyskad wynik optymalizacji w komórkach B1, B2, B3. Wartości pokazane na rysunku Rys. 6 nie są wartościami po optymalizacji! Proszę też pamiętad, że algorytm optymalizacji w Solverze zaczyna poszukiwania od wartości które znajdują się w komórkach. Jeżeli optimum jest zbyt daleko to algorytm może nie znaleźd szukanych liczb. Wtedy należy zmienid wartości początkowe i spróbowad kolejny raz. Na ogół wiemy z jakimś przybliżeniem (np. co do rzędu wielkości) jakich wartości optymalnych należy się spodziewad, więc od nich startujemy. Ale w ogólności szukanie optymalnych parametrów może wymagad wielu prób.

Zadania do wykonania w ramach projektu

1) Zapoznanie się z modelem NPP i porównanie z innymi modelami (równanie Nernsta, wzór Nikolskiego  Eisenmana, wzór Hendersona).

W tym miejscu wpisujemy stężenie

jonu I

w próbce. Przykładowy

zapis 1E2 to inaczej 100. Obie

formy zapisu są dopuszczalne.

Krzysztof Szyszkiewicz, AGH

-16-

2) Wykonanie symulacji dla przykładowej pary roztworów wg schematu pracy *1+, tabela 1, str. 2446.

3) Wykonanie symulacji krzywej kalibracyjnej (przykład: praca *1+, str. 2446-2447, Rys. 2, 3) 4) Wyznaczenie współczynnika selektywności ze wzoru Nikolskiego  Eisenmana w oparciu o

wyliczoną krzywą kalibracyjną.

5) Wykonanie symulacji ewolucji potencjału w czasie dla styku cieczy (Liquid junction) (przykład: praca *1+, str. 2446, Rys. 1).

Literatura

[1] T. Sokalski, P. Lingenfelter, and A. Lewenstam, Numerical Solution of the Coupled Nernst- Planck and Poisson Equations for Liquid Junction and Ion Selective Membrane Potentials. J.

Phys. Chem. B (2003), 107, 24432452.

[2] J.J. Jasielec, R. Filipek, K. Szyszkiewicz, M. Danielewski, Computer simulations of electrodiffusion problems based on Nernst–Planck and Poisson equations. Comp. Mat. Sci. 63 (2012), 75-90.

[4] W. Kucza, M. Danielewski, EIS simulations for ion-selective site-based membranes by a numerical solution of the coupled Nernst–Planck–Poisson equations, Electrochem. Comm. 8, (2006), 416420.

[5] A. Kisza, Elektrochemia, cz. II – Elektrodyka, WNT, Warszawa, 2001.