Index of /rozprawy2/10180

Pełen tekst

(1)Akademia Górniczo-Hutnicza im. Stanisława Staszica w Krakowie Wydział Inżynierii Materiałowej i Ceramiki Katedra Fizyko-Chemii Ciała Stałego. Rozprawa doktorska. Metoda elementów skończonych w elektrochemii. Krzysztof Szyszkiewicz-Warzecha. Promotor: prof. dr hab. inż. Marek Danielewski. Kraków, 2009.

(2) Chciałbym wyrazić serdeczne podziękowania dla Panów prof. dr hab. Marka Danielewskiego i dr hab. Roberta Filipka za merytoryczną oraz redakcyjną pomoc w czasie powstawania tej pracy, a także Pani dr Beacie Paczosie-Bator za pomoc w części eksperymentalnej..

(3) Spis treści. 1 Wprowadzenie Wykaz symboli i skrótów . . . . . . . . . . . . . . . . . . . . . . . . 2 Wstęp i cel pracy 2.1 Przeskalowanie i zmienne bezwymiarowe . 2.2 Zakres stosowalności modelu NPP . . . . . 2.3 Przegląd numerycznych metod w NPP . . 2.3.1 Równania NPP niezależne od czasu 2.3.2 Równania NPP zależne od czasu . 2.3.3 Stosowane metody numeryczne . .. . . . . . .. . . . . . .. . . . . . .. . . . . . .. . . . . . .. . . . . . .. . . . . . .. . . . . . .. . . . . . .. . . . . . .. 3 Elektrodyfuzja w stanie stacjonarnym 3.1 Wprowadzenie . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 NPP w stanie równowagi . . . . . . . . . . . . . . . . . . . . 3.3 NPP w stanie stacjonarnym . . . . . . . . . . . . . . . . . . 3.3.1 Planck . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Henderson . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Warunki brzegowe w elektrochemii . . . . . . . . . . 3.3.4 Potencjał membranowy w przypadku elektrolitu binarnego w stanie stacjonarnym . . . . . . . . . . . . . . 3.4 Nowa metoda . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Wyniki obliczeń . . . . . . . . . . . . . . . . . . . . . 4 Metoda elementów skończonych dla układu NPP 4.1 Wprowadzenie do metody elementów skończonych . . 4.2 Słabe sformułowanie ewolucyjnego zagadnienia NPP nym wymiarze . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Prawo zachowania masy . . . . . . . . . . . . 4.2.2 Wyrażenie na strumień składników . . . . . . 4.2.3 Prawo Gaussa . . . . . . . . . . . . . . . . . . 4.2.4 Równanie stanu . . . . . . . . . . . . . . . . . 4.2.5 Warunki brzegowe i początkowe . . . . . . . . 4.2.6 Niewiadome . . . . . . . . . . . . . . . . . . . 4.2.7 Warunki początkowe . . . . . . . . . . . . . . 4.2.8 Przeskalowanie zmiennych i równań . . . . . .. 3. 5 6. . . . . . .. 7 16 19 22 22 26 28. . . . . . .. 34 34 35 41 42 48 50. . 55 . 59 . 67. 74 . . . . . 74 w jed. . . . . 82 . . . . . 82 . . . . . 82 . . . . . 83 . . . . . 84 . . . . . 84 . . . . . 85 . . . . . 85 . . . . . 85.

(4) SPIS TREŚCI. 1. 5 Cześć eksperymentalna 94 5.1 Membrany jonoselektywne . . . . . . . . . . . . . . . . . . . . 94 5.2 Aparatura i odczynniki . . . . . . . . . . . . . . . . . . . . . . 96 5.3 Przygotowanie elektrod . . . . . . . . . . . . . . . . . . . . . . 97 6 Obliczenia numeryczne. 100. 7 Podsumowanie i wnioski. 120. 8 Dalsze kierunki badań. 122. A Dodatek do metody elementów skończonych 124 A.1 Elementy skończone i obliczenia całek . . . . . . . . . . . . . . 124 A.2 Całki do metody elementów skończonych . . . . . . . . . . . . 124 B Różnice skończone, stan stacjonarny. 131. C Bibliografia. 142.

(5) Spis rysunków 1 2 3. 4. Trójwymiarowy model membrany otoczonej rozworami elektrolitów. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Jednowymiarowy model membrany otoczonej rozworami elektrolitów. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Ilustracja do zjawiska nasycenia pola elektrycznego. Potencjał ηW w obszarze W jest dodatni. Nawet, gdy ηW → ∞, to potencjał poza obszarem W zmierza w każdym punkcie x do skończonej wartości granicznej. . . . Nasycenie polem elektrycznym w jednym wymiarze. Bezwymiarowe pole elektryczne u = u(x, η) (rozwiązanie problemu (3.14)) dla różnych wartości potencjału η na brzegu u(0) = η. Nawet, gdy η → ∞, to potencjał poza brzegiem {0} zmierza w każdym punkcie x do skończonej wartości granicznej u ˜(x) (linia przerywana). . . . . . . . . . . . . . . . . . .. 39. 41. Kνz11 Azν22. 5. Membrana otoczona rozworami elektrolitów o koncentracjach cL i cR . . . . . . . . . . . . . . . . . . . . . . . . . . . 56. 6. Wykres dokładnego rozwiązania problemu ²y 0 = − sin(x) cos(x)y, y(0) = 1 dla ² = 10−2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 Wykres rozwiązania problemu ²y 0 = − sin(x) cos(x)y, y(0) = 1 dla ² = 10−2 procedurą numeryczną (niejawną). Przy prawym końcu rozwiązanie całkowicie odbiega od dokładnego (por. Rys. 6). . . . . . . . . . . . . 62 Porównanie profili koncentracji i potencjału elektrycznego dla stanu stacjonarnego otrzymanych dwiema metodami. Obliczenia wykonane dla ¯ 1 = 10, 98, D ¯ 2 = 2, 01, z1 = 1, z2 = −1. . . . . . 68 parametrów: λ = 20, D Koncentracje i potencjał elektryczny obliczone dla stanu stacjonarnego przy pomocy metody II. Obliczenia wykonane dla następujących wartości ¯ 1 = 10, 98, D ¯ 2 = 2, 01, z1 = 1, z2 = −1. . . . . 69 parametrów: λ = 105 , D Porównanie profili koncentracji i potencjału elektrycznego dla stanu stacjonarnego w układzie trójskładnikowym. Obliczenia wykonane dla parametrów: ¯ 1 = 10, 98, D ¯ 2 = 1, 35, D ¯ 3 = 2, 01, z1 = 1, z2 = 1 z3 = −1. . . 69 λ = 20, D Krzywa kalibracji dla PCV DOS ISE - porównanie obliczeń numerycznych z wynikami eksperymentu. (Eksperyment przeprowadzony dzięki pomocy pani dr Beaty Paczosy-Bator, Katedra Chemii Analitycznej, AGH.) ∆φ — potencjał membranowy. . . . . . . . . . . . . . . . . . 69 Obliczone i eksperymentalne ewolucje czasowe potencjału dla PCV DOS ISE w zewnętrznych roztworach elektrolitów NaCl i KCl. . . . . . . . . 70 Profil bezwymiarowej koncentracji c = c(x). Przypadek symetryczny: c(0) = c0 = c1 = c(1). Obliczenia wykonane dla wartości: c0 = 10−2 , 100 , 102 , 103 , 104 . . . . . . . . . . . . . . . . . . . . . . . . . . . 72. 7. 8. 9. 10. 11. 12 13. 2.

(6) SPIS RYSUNKÓW 14. 15 16 17. 18 19. 20. 21. 22. 23. 24. 25. Profil bezwymiarowego pola elektrycznego E = E(x). Przypadek symetryczny: c(0) = c0 = c1 = c(1). Obliczenia wykonane dla wartości: c0 = 10−2 , 100 , 102 , 103 , 104 . . . . . . . . . . . . . . . . . . . . . .. 3. 73. Funkcje poligonalne na odcinku — baza do metody elementów skończonych. Dla funkcji ϕi mamy ϕi (xi ) = 1. . . . . . . . . . 77 Przykład triangulacji obszaru Ω na płaszczyźnie R2 . . . . . . . 79 Krzywe kalibracji dla elektrody jonoselektywnej PCV DOS – porównanie wyników obliczń oraz eksperymentu. φ – potencjał membranowy, pK = − log cK + . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Potencjał membranowy dla roztworu 0,5M CaCl2 w funkcji czasu dla różnych grubości membrany d: 500, 1000 i 2000 µm. . . . . . . . . . . 107 Obliczone potencjały dyfuzyjne dla różnych roztworów o różnym stężeniu dla liczby punktów dyskretyzacji N = 200 w temperaturze 298 K. FEM — metoda elementów skończonych. . . . . . . . . . . . . . . . . . . 109 Obliczone potencjały dyfuzyjne dla różnych roztworów o różnym stężeniu dla liczby punktów dyskretyzacji N = 400 w temperaturze 298 K. FEM — metoda elementów skończonych. . . . . . . . . . . . . . . . . . . 110 Obliczone potencjały dyfuzyjne dla różnych roztworów o różnym stężeniu dla liczby punktów dyskretyzacji N = 800 w temperaturze 298 K. FEM — metoda elementów skończonych. . . . . . . . . . . . . . . . . . . 111 Obliczone potencjały dyfuzyjne dla różnych roztworów o różnym stężeniu dla liczby punktów dyskretyzacji N = 800 w temperaturze 298 K. FEM — metoda elementów skończonych. . . . . . . . . . . . . . . . . . . 112 Ewolucja czasowa potencjału dyfuzyjnego dla układu typu biionic NaCl(0,1 M)|KCl(0,1 M). Wszystkie stałe heterogeniczne takie same i bardzo duże (ki = 103 m/s). Obliczenia wykonane metodą elementów skończonych dla liczby funkcji bazowych N = 1000 (przypadek uwzględniajacy dryft i przypadek bez dryftu). Por. [43] str. 4229, Rysunek 4 – krzywa na dole (K = 157, 08) odpowiada stałej względnej dielektrycznej ²r = 78, 6 (woda). 113 Ewolucja czasowa potencjału dyfuzyjnego dla układu trzech jonów: z1 = +2, z2 = +1, z3 = −1. Współczynniki dyfuzji D1 = 0, 798 · 10−5 , D2 = 1, 98 · 10−5 , D3 = 2, 01 · 10−5 . Koncentracje c1,L = 0 M , c2,L = c3,L = 0, 1 M , c1,R = 0, 5 M , c2,R = 0 M , c3,R = 1 M . Wszystkie stałe heterogeniczne takie same i bardzo duże (ki = 103 m/s). Obliczenia wykonane metodą elementów skończonych (przypadek uwzględniajacy dryft i przypadek bez dryftu). Liczba funkcji bazowych N = 1000 w metodzie FEM. . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 Układ bi-ionic – ewolucja czasowa profilu koncentracji dla kationu M + dla wybranych czasów (w sekundach). Stałe heterogeniczne takie same i bardzo duże (ki = 102 m/s). Obliczenia wykonane metodą elementów skończonych (przypadek bez dryftu). Liczba funkcji bazowych N = 1000 w metodzie FEM. . . . . . . . . . . . . . . . . . . . . . . . . . . 115.

(7) SPIS RYSUNKÓW 26. 27. 28. 29. 30. 4. Układ bi-ionic – ewolucja czasowa profilu koncentracji dla kationu N + dla wybranych czasów (w sekundach). Stałe heterogeniczne takie same i bardzo duże (ki = 102 m/s). Obliczenia wykonane metodą elementów skończonych (przypadek bez dryftu). Liczba funkcji bazowych N = 1000 w metodzie FEM. . . . . . . . . . . . . . . . . . . . . . . . . . . 116 Układ bi-ionic – ewolucja czasowa profilu koncentracji dla kationu R− dla wybranych czasów (w sekundach). Stałe heterogeniczne takie same i bardzo duże (ki = 102 m/s). Obliczenia wykonane metodą elementów skończonych (przypadek bez dryftu). Liczba funkcji bazowych N = 1000 w metodzie FEM. . . . . . . . . . . . . . . . . . . . . . . . . . . 116 Obliczona i ekeperymetalna evolucja potencjału φ(t) membrany jonoselektywnej PCV DOS. Zewnętrzne roztwory to NaCl i KCl. Liczba funkcji bazowych N = 1000 w metodzie FEM. . . . . . . . . . . . . . . . . 117 Ewolucja czasowa potencjału membranowego. Obliczenia wykonane metodą elementów skończonych (przypadek uwzględniajacy dryft i przypadek bez dryftu. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 Ewolucja czasowa potencjału membran. Obliczenia wykonane metodą elementów skończonych (przypadek uwzględniajacy dryft i przypadek bez dryftu). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119.

(8) 5.

(9) ROZDZIAŁ 1. WPROWADZENIE. 6. Rozdział 1 Wprowadzenie. Wykaz symboli i skrótów. N Z R Rn = R × . . . × R int(Ω) ¯ Ω ∂Ω f :X→Y dom(f ) ϕ P x · y = nj=1 xj yj dy , y0 ³ dx ´ ∂f ∂f , . . . , ∇f = ∂x ∂x 1 Pnn div u = ∇ · u = j=1 P 2 ∆f = nj=1 ∂∂xf2 1. H (0, d) µi µoi Bi E V, φ Vo Di Bi ci ai γi Ji zi k F NA R = kNA T τ = kT e M. j. – – – – – – – – – – – –. zbiór liczb naturalnych {1, 2, 3, . . .} zbiór liczb całkowitych zbiór liczb rzeczywistych n-wymiarowa przestrzeń kartezjańska wnętrze zbioru Ω domknięcie zbioru Ω brzeg zbioru Ω odwzorowanie (funkcja) ze zbioru X w Y dziedzina odwzorowania (funkcji) f : X → Y funkcja próbna (do metody elementów skończonych) standardowy iloczyn skalarny wektorów x i y w Rn pochodna funkcji jednej zmiennej y : I → Rn gdzie I ⊂ R. – gradient funkcji f ∂uj ∂xj. – dywergencja funkcji wektorowej u = (u1 , . . . , un ) – laplasjan funkcji f – przestrzeń funkcyjna Sobolewa na odcinku [0, d] – potencjał chemiczny i−tego składnika – potencjał chemiczny i−tego składnika w stanie standardowym – ruchliwość i−tego składnika – wektor natężenia pola elektrycznego – potencjał elektryczny – potencjał standardowy (względem NHE) – współczynnik dyfuzji i−tego składnika – ruchliwość i−tego składnika – koncentracja (stężenie) i−tego składnika – aktywność i−tego składnika – współczynnik aktywności i−tego składnika – strumień i−tego składnika – liczba ładunkowa i−tego składnika – stała Boltzmanna – stała Faradaya – liczba Avogadra – stała gazowa – temperatura w skali Kelvina – temperatura fundamentalna (miara energii termicznej; jednostka J) – ładunek elementarny (protonu) – koncentracja molowa (mol/dm3 ).

(10) Rozdział 2 Wstęp i cel pracy. Celem pracy jest (I) Zastosowanie metody elementów skończonych do numerycznego rozwiązywania równań opisujących transport masy i ewolucję pola elektrycznego w układach membranowych. (II) Zbadanie wpływu prędkości dryftu Darkena na model transportu jonów w oparciu o układ równań Nernsta-Plancka-Poissona (NPP), a w szczególności jego wpływ na potencjał stacjonarny i na sposób dochodzenia do tego stanu. (III) Rozwiązanie stacjonarnego zagadnienia NPP – ze szczególnymi warunkami brzegowymi – zaproponowanego w pracy [77] Układy elektrochemiczne z membranami jonoselektywnymi znajdują bardzo szerokie zastosowania [7], [13]. O końca lat sześćdziesiątych XX wieku elektrody jonoselektywne (ang. ion selective electrodes, ISE ) były i są jednym z najważniejszych narzędzi w chemii analitycznej, a w szczególności analizie klinicznej i środowiskowej [70]. Aby układy takie mogły być w pełni i efektywnie wykorzystane, konieczna jest możliwość obliczenia rozkładu koncentracji składników elektrolitu oraz rozkładu potencjału (lub pola elektrycznego) w takim układzie. Dotyczy to zarówno stanu stacjonarnego jak i ewolucji czasowej układu. W ogólności istnieją dwa rodzaje transportu przez membrany [34]. Jeden występuje wtedy, gdy membrana jest jednorodna, a więc może być traktowana jako pojedyncza faza. Drugi występuje, gdy w membranie znajdują się wyraźne struktury zwane kanałami (typowy promień powyżej 10−9 m), poza którymi membrana jest nieprzepuszczalna lub ogólniej jej przepuszczal7.

(11) ROZDZIAŁ 2. WSTĘP I CEL PRACY. 8. ność poza kanałami jest istotnie mniejsza. Taka struktura membrany oznacza, że jest to medium składające się z dwóch faz, a jony mogą znajdować się w każdej z nich. Granica podziału pomiędzy tymi dwoma przypadkami nie musi być ostra. Jeżeli gęstość kanałów będzie się zwiększać i jednocześnie ich średnica zmniejszać, to opis drugi przechodzi w opis pierwszy. W tej pracy zajmujemy się membranami, w których występuje pierwszy rodzaj transportu. Matematyczny opis takich układów oparty jest na kilku prawach i założeniach. Niektóre z nich mają charakter uniwersalnych praw fizyki inne są świadomymi przybliżeniami spisującymi się dobrze tylko w ramach pewnych założeń. Podstawowe równania rządzące dyfuzją i migracją naładowanych cząstek zostały sformułowane przez Plancka i Nernsta w pracach, które ukazały się w latach 1888-1890. Główne elementy tego opisu stanowią: • Równania wyrażające prawo zachowania masy poszczególnych składników (mogą tu wystąpić człony opisujące reakcje – czyli produkcję lub konsumpcję składników). • Równanie opisujące zjawiska elektrodynamiczne, gdyż układ zawiera nośniki ładunku (w elektrochemii są to przede wszystkim jony). W przybliżeniu „quasistatycznym” jest to prawo Gaussa (często używane w formie równania Poissona, gdy posługujemy się potencjałem elektrycznym zamiast natężeniem pola elektrycznego). • Wyrażenie opisujące postać strumienia poszczególnych składników układu (zwane równaniem konstytutywnym). • W przypadku, gdy w układzie występuje dryft (zwany w elektrochemii konwekcją), należy użyć dodatkowego równania. W najogólniejszym przypadku będzie to hydromechaniczny układ równań Naviera-Stokesa. W pewnym sytuacjach można wykorzystać ideę pochodzącą od Darkena (która pojawiła się w kontekście dyfuzji w ciałach stałych Darken [22] i została rozwinięta dla układów wieloskładnikowych (Holly, Danielewski [49]). Sprowadza się ona do przyjęcia stałości całkowitej koncentracji. Warunek ten wymusza pojawienie się dryftu. W niniejszej pracy warunek ten będzie wykorzystany..

(12) ROZDZIAŁ 2. WSTĘP I CEL PRACY. 9. • Charakter zjawisk opisujących proces transportu na granicy faz. Ostatecznym matematycznym wyrazem tego opisu będą warunki brzegowe narzucone na równania różniczkowe w modelu. Mogą to być proste warunki typu Dirichleta (uwzględniające tylko wartości funkcji lub warunki typu Neumanna uwzględniające także gradienty funkcji). Wymienione powyżej elementy (wraz z warunkami początkowymi) wystarczają w zasadzie do poprawnego matematycznie sformułowania problemu opisującego układ elektrochemiczny. Użyta w poprzednim zdaniu fraza „w zasadzie” ma podkreślić, że z matematycznego punktu widzenia poprawne sformułowanie problemu w teorii równań różniczkowych oznacza, że układ równań ma własność istnienia i jednoznaczności rozwiązania oraz stabilności wzglądem zaburzenia warunków początkowych, brzegowych i ogólnie parametrów występujących w równaniach (przynajmniej w pewnym zakresie zmienności tych wielkości). Mówimy wtedy, że problem jest dobrze postawiony (ang. well posed) [30, 81]. Dla ogólnych równań opisujących elektrodyfuzję z uwzględnieniem prędkości konwekcyjnych (co w ogólnym przypadku wymaga zastosowania równań Naviera-Stokesa) taka pełna teoria aktualnie nie istnieje. W literaturze występuje dużo wyników cząstkowych (np. Biler, Habish, Nadzieja [6], Roubicek [82]). Wyniki zaprezentowane w monografii [82] dotyczą istnienia rozwiązań układu równań typu Naviera-StokesaNernsta-Planka-Poissona dla dowolnej liczby składników w geometrii dwu i trójwymiarowej. Układ ten uwzględnia efekty hydromechaniczne, elektrodyfuzyjne jak i reakcyjne - jest więc bardzo ogólny. Ograniczeniem są warunki brzegowe (tylko typu Dirichleta) oraz równość wszystkich współczynników dyfuzji (w pracy autor używa ruchliwości, co dzięki relacji Einsteina jest równoważne z opisem przy pomocy współczynników dyfuzji). Ponadto dowodzone jest tylko istnienie rozwiązania (bez jednoznaczności). Z drugiej strony elektrochemika interesuje przede wszystkim symulacja numeryczna uzyskanego modelu. W praktyce dobrze sprawdza się zasada, że jeżeli takie symulacje przeprowadzone różnymi metodami, dla różnych dyskretyzacji przestrzennych i czasowych dają rezultaty zgodne, to należy przyjąć, że model jest wiarygodny i poprawnie postawiony matematycznie. Oczywiście jest to pewne założenie metodologiczne – zresztą bardzo powszechne w zastosowaniach – nad słusznością którego można dyskutować. Jednakże sukcesy przyjętej metodologii są silnym potwierdzeniem jej słuszności..

(13) ROZDZIAŁ 2. WSTĘP I CEL PRACY. 10. Dalej zostaną przedstawione dokładnie podstawowe idee fizyczne i pojęcia zarysowane powyżej, a do których będziemy się odwoływać w tej pracy. Pełny przegląd tych zagadnień w kontekście elektrochemicznym można znaleźć w monografiach Bockris, Redy [8] lub Bard, Faulkner [3]. W zagadnieniach elektrodyfuzji, które można sprowadzić do geometrii jednowymiarowej (np. membrany jonoselektywne czy zjawiska elektrodializy) wykorzystujemy wyidealizowany model membrany. Rozwór elektrolitu otacza cienką membranę, która jest utworzona przez dwie płaskie równoległe półprzepuszczalne płaszczyzny graniczne. Słowo „półprzepuszczalny” oznacza, że membrana pozwala przenikać różnym jonom w różnym stopniu. Zależy to od ich znaków, lub ogólniej od ich tożsamości chemicznej. Ponadto zakładamy, że pole elektryczne może być przyłożone tylko w kierunku prostopadłym do płaszczyzn granicy. Nośnikami ładunku są tu wyłącznie jony (przewodnik jonowy) – Rysunek 1.. Rysunek 1: Trójwymiarowy model membrany otoczonej rozworami elektrolitów. Ze względu na zakładaną jednorodność w dwóch prostopadłych kierunkach rozważany model redukuje się do jednego wymiaru. Schematycznie może go zobrazować jak na Rysunku 2.. Rysunek 2: Jednowymiarowy model membrany otoczonej rozworami elektrolitów. Dążenie do stanu równowagi termodynamicznej w układzie, który zawiera składniki jonowe może spowodować powstanie niezerowego potencjału elektrycznego. Początkowa różnica potencjałów chemicznych poszczególnych składników w różnych obszarach będzie siłą napędową prowadzącą do.

(14) 11. ROZDZIAŁ 2. WSTĘP I CEL PRACY. separacji ładunku elektrycznego i uformowania się pewnego rozkładu koncentracji oraz potencjału elektrycznego, przy których potencjał chemiczny będzie stały w całym obszarze. W szczególności podstawowe zagadnienie polega na wyznaczeniu potencjału membranowego ∆φM (t) = φR (t) − φL (t) jako funkcji parametrów opisujących układ – zarówno w stanie stacjonarnym (niezależnym od czasu) jak i niestacjonarnym (dynamicznym, tj. zależnym od czasu).1 Zapiszmy teraz ten słowny opis w terminach równań klasycznej kinetyki fenomenologicznej. Rozważmy więc mieszaninę r składników każdy obdarzony liczbą ładunkową zi , i = 1, . . . , r (są to małe liczby całkowite wyrażające ładunek cząsteczki w jednostkach ładunku elementarnego (protonu); dla obojętnych cząsteczek zi = 0). Koncentracja i-tego składnika opisana jest funkcją ci = ci (x, t) [mol/dm3 ], przy czym zakładamy w tej pracy, że w roztworach otaczających membranę koncentracje te są utrzymywane na stałym poziomie ci,L oraz ci,R . Do opisu transportu materii używa się strumienia i-tego składnika Ji = Ji (x, t) (ilość materii, która przepłynie przez jednostkową powierzchnię prostopadłą do wektora prędkości w jednostkowym czasie). Produkcja lub konsumpcja poszczególnych składników opisana jest funkcją Ri (człon reakcyjny). Zasada zachowania masy i−tego składnika wyrażona w formie różniczkowej ma przy powyższych oznaczeniach postać (2.1). ∂ci + divJi = Ri (c1 , . . . , cr ) ∂t. (i = 1, . . . , r).. Zwróćmy uwagę, że równanie to jest fundamentalnym prawem fizyki, a nie jedynie pewnym przybliżeniem. Oczywiście obowiązuje on w ramach pewnego paradygmatu opisu ośrodka. Głębsze i ogólniejsze podejście do tego 1. Klasyczna elektrochemia była głównie zainteresowana zjawiskami na granicy faz: elektroda (przewodnik elektronowy, np. metal lub grafit) a elektrolit (przewodnik jonowy). Mniejsza waga była przykładana do powstawania potencjału elektrycznego w samej fazie elektrolitu, co związane jest z dość powszechnie przyjmowanym założeniem tzw. elektroobojętności oraz z dużym znaczeniem reakcji redox i ich potencjałotwórczej roli. Jednakże dla membran jonoselektywnych takie podejście – redukujące powstawanie potencjału do zjawisk na granicy faz – jest nieusprawiedliwione, i może czasami prowadzić do poważnych błędów. W szczególności badania prowadzone przez Lewenstama et. al. [89], [90] wykazały, że udział potencjału na granicy faz (ang. boundary potential ) w całym potencjale może stanowić nawet mniej niż 50%! Sytuacja taka występuje np. w elektrodach jonoselektywnych..

(15) 12. ROZDZIAŁ 2. WSTĘP I CEL PRACY. paradygmatu i postaci równania ciągłości można znaleźć w pracach Danielewski, Wierzba [27] i [28], gdzie prezentowany jest nowy model zależności pomiędzy dyfuzją i dryftem poprzez gęstość objętości. Wspomniani autorzy próbują rozwiązać ogólny problem zunifikowanego opisu procesów transportu w kontekście termodynamiki procesów nieodwracalnych łącząc w ramach jednego modelu dryft i dyfuzję.2 Ponieważ w prezentowanej rozprawie nie będą uwzględniane reakcje, więc dalej będzie używane tylko równanie bilansu masy bez członów reakcyjnych (Ri = 0) (2.2). ∂ci = −divJi ∂t. (i = 1, . . . , r).. Kolejnym krokiem niezbędnym do bardziej szczegółowego opisu układu jest użycie jakiegoś prawa konstytutywnego na strumień Ji , które w naszym przypadku powinno wyrażać zależność strumienia od pozostałych wielkości opisujących układ, tj. koncentracji ci = ci (x, t), natężenia pola elektrycznego E = E(x, t) oraz prędkości dryftu mieszaniny υ dryf t = υ dryf t (x, t). Dla rozcieńczonych roztworów, lub ogólnie dla roztworów idealnych (idealnych sensie termodynamicznym, tzn., gdy aktywności składników można zastąpić ich koncentracjami) oraz w warunkach stanu lokalnej równowagi powszechnie jest akceptowane wyrażenie Nernsta-Plancka, które (z uwzględnieniem prędkości unoszenia) ma postać [3] (2.3). Ji = −Di ∇ci +. F zi Di ci E + ci υ dryf t , RT. gdzie Di jest współczynnikiem dyfuzji i−tego składnika (mogącym w ogólności zależeć od koncentracji poszczególnych składników Di = Di (c1 , . . . , cr )), F jest stałą Faradaya, R uniwersalną stałą gazową, a T temperaturą w skali Kelvina. 2. Zasadniczy problem z dyfuzją w mechanice ośrodków ciągłych polega na poprawnej definicja strumienia dyfuzyjnego. W literaturze występuje wiele takich strumieni, które określa się mianem dyfuzyjnych. Przykładowo w monografii Joel L. Plawsky [79] podanych jest sześć – w zależności od tego jakie stężenie i jakie prędkości bierze się za podstawę do definicji strumienia dyfuzyjnego. Ogólnie Strumień=(stężenie)×prędkość. Na pojawiające się w związku z tym trudności i sprzeczności zwracali już uwagę np. Brenner [9] i Öttinger [76]. Problem ten rozwiązali wspomniani w gł. tekście autorzy ([28]) poprzez wykazanie, że właściwym układem odniesienia, w którym należy definiować strumienie dyfuzyjne jest tzw. prędkość objętościowa. W tym celu autorzy wyprowadzają równanie ciągłości objętości..

(16) ROZDZIAŁ 2. WSTĘP I CEL PRACY. 13. Równania (2.2) i (2.3) należy uzupełnić o prawo Gaussa wiążące natężenie pola elektrycznego z rozkładem ładunków poprzez koncentracje ci (gdy w opisie używamy potencjału elektrycznego φ związanego z polem elektrycznym równością definicyjną E = −∇φ, wtedy jest to równanie Poissona) oraz równanie Naviera-Stokesa, gdy założymy, że rozwór jest newtonowskim płynem lepkim. Prawo Gaussa [36] (2.4). div(εr ε0 E) = q,. gdzie εr jest względną przenikalnością elektryczną ośrodka (czasami zwaną stałą dielektryczną ośrodka), ε0 jest stałą dielektryczna próżni, q jest przestrzenną gęstością ładunku, pochodzącą od mogących się poruszać naładowanych jonów oraz ewentualnego stałego ładunku (tła) w membranie opisanego funkcją niezależną od czasu N = N (x) (np. dla stałych membran jonoselektywnych). Tak więc gęstość ta może być zapisana następująco (2.5). q(x, t) = F. r X. zi ci (x, t) + N (x).. i=1. W oparciu o równania (2.2), (2.4), (2.5) można wyprowadzić bardzo ważną zależność na gęstość całkowitego prądu, w której występuje m. in. tzw. prąd przesunięcia (ang. displacement current). Szczególny przypadek tej zależności jest często wykorzystywany w symulacjach numerycznych zamiast wprost równania Gaussa (lub w wersji z potencjałem - równania Poissona) ([16], [11], [89], [90], [61]). Po raz pierwszy została ona użyta w pracy Cohen, Cooley [16] (1965), ale autorzy nie podają żadnego wyprowadzenia uważając prawdopodobnie, że interpretacja fizyczna podanej przez nich równości jest wystarczającym dowodem. Poniżej przedstawiamy jednak pełne wyprowadzenie. W tym celu pomnóżmy każde z równań bilansu masy (2.2) przez F zi i dodajmy je stronami (dla i = 1, . . . , r), aby otrzymać Ã r ! r X ∂ X F zi ci = −div Ji . ∂t i=1 i=1 Teraz zauważmy, że wyrażenie w nawiasie po lewej stronie to gęstość ładunku, którą można zastąpić na podstawie prawa Gaussa (2.4) przez div(εr ε0 E). Tak więc mamy. Ã r ! X ∂ div(εr ε0 E) = −div F zi Ji , ∂t i=1.

(17) 14. ROZDZIAŁ 2. WSTĘP I CEL PRACY. czyli po przegrupowaniu. Ã. r X ∂E div εr ε0 +F z i Ji ∂t i=1. (2.6). ! = 0.. Wykorzystując postać strumienia (2.3) możemy te rachunki podsumować następująco (2.7). divI(t) = 0,. gdzie (2.8) I(t) = −F. r X. Ã zi Di ∇ci +. i=1. ! Ã r ! r X F2 X 2 ∂E z Di ci E + εr ε0 + F zi ci υ, RT i=1 i ∂t i=1. jest gęstością całkowitego prądu). Widać, że gęstość ta jest sumą następujących składników: µ F2 RT. 1. prąd przewodzenia: 2. prąd dyfuzyjny: −F. r P. r P i=1. 3. prąd konwekcyjny:. F. E,. zi Di ∇ci ,. i=1. µ. ¶ zi2 Di ci. r P. ¶ zi ci υ,. i=1. 4. prąd przesunięcia3 : εr ε0 ∂E . ∂t Równość przedstawiona w (2.8) może być interpretowana jako uogólnione prawo Ohma dla układów jonowych. W szczególności czynnik (2.9) 3. r F2 X 2 z Di ci = σ, RT i=1 i. Nazwa „prąd przesunięcia” (ang. displacement current) pochodzi jeszcze od J. C. Maxwella (po raz pierwszy użył jej w pracy “On Physical Lines of Force” (1861-1862)), który wprowadził wyraz ² ∂E ∂t do równania wyrażającego prawo Ampère’a, uzyskując równanie rotB = µJ + µ² ∂E (B – wektor indukcji magnetycznej, J – wektor gęstości prądu). ∂t Wyraz ten jest niezbędny m. in. do tego, aby było spełnione prawo zachowania ładunku. Taką zresztą motywację wprowadzenia tego wyrazu przypisuje się na ogół Maxwellowi co, jak to często bywa w podobnych przypadkach w historii nauki, nie jest prawdą. Oczywiście Maxwell zdawał sobie sprawę, z tego, że prawo zachowania ładunku jest dzięki takiej postaci prawa Ampère’a spełnione, ale jego pierwotna argumentacja opierała się na mechanicznej koncepcji eteru. Kontekst historyczny oraz powody dla których Maxwell wprowadził prąd przesunięcia można znaleźć np. w artykule A. M. Bork’a [4]. Chociaż koncepcja eteru w fizyce została w zasadzie odrzucona, to jednak prowadzone są przez niektórych uczonych badania nad współczesną analogią tej koncepcji, która nosi nazwę kryształu Kleinerta-Plancka..

(18) ROZDZIAŁ 2. WSTĘP I CEL PRACY. 15. jest odpowiednikiem oporu omowego. Równanie (2.6) jest bardzo użyteczne zwłaszcza w układach, które można sprowadzić do geometrii liniowej jednowymiarowej. W tej sytuacji dywergencja jest zwykłą pochodną względem jednej zmiennej przestrzennej, więc równanie ma postać Ã ! r X ∂E ∂ εr ε0 +F zi Ji = 0. ∂x ∂t i=1 Stąd oczywiście wynika, że funkcja pod znakiem pochodnej jest stała, a zatem (2.10). r X ∂E εr ε0 = I(t) − F z i Ji , ∂t i=1. gdzie prąd I = I(t) nie zależy od położenia tylko ewentualnie od czasu. Ostatnie równanie – o którym już wspomnieliśmy wyżej – to równanie opisujące hydromechanikę układu, a więc np. równanie Naviera-Stokesa ¶ µ ∂υ (2.11) ρm + υ∇υ = −∇p + η∆υ + f, ∂t gdzie υ = υ(x, t) jest polem prędkości płynu, ρm = ρm (x, t) jest gęstością masy roztworu, p oznacza ciśnienie, η lepkość kinematyczną, a f = f (x, t) jest gęstością sił objętościowych działających na roztwór (płyn). Najczęściej siły masowe w układach jonowych składają się z części grawitacyjnej (fG ) oraz części elektrycznej (fE ): (2.12). f = fG + fE , fG = ρm g, fE = qE,. gdzie g oznacza stałe przyspieszenie grawitacyjne (w laboratorium, w którym znajduje się nasz układ przestrzenna zmienność pola grawitacyjnego Ziemi jest zaniedbywana). Zauważmy, że posługując się potencjałem elektrycznym oraz wykorzystując równanie Gaussa (2.4) (a dokładniej Poissona) oraz definicję potencjału elektrycznego E = −∇φ możemy gęstość sił elektrycznych fE = qE zapisać następująco (por. [36], str. 116) (2.13). fE =. 1 (∆φ)(∇φ). εr ε0.

(19) ROZDZIAŁ 2. WSTĘP I CEL PRACY. 2.1. 16. Przeskalowanie i zmienne bezwymiarowe. Równania opisujące układy fizyczne zawierają wielkości wyrażone w zmiennych wymiarowych i często zawierają wiele parametrów. Przez wprowadzenie bezwymiarowych zmiennych zależnych i niezależnych możemy liczbę tych parametrów zmniejszyć, gdyż tylko pewne ich kombinacje pojawią się w uzyskanych przeskalowanych równaniach. Ponadto przez wybranie właściwych skal przy takim przekształcaniu równań możliwe jest oszacowanie względnych rozmiarów różnych składników tych równań. Procedura taka może prowadzić do uproszczenia zapisu modelu, co z kolei upraszcza zarówno analizę teoretyczna jak i symulacje numeryczne. Od strony matematycznej kroki, które należy podjąć, aby wykonać przeskalowanie są bardzo proste. Sprowadzają się one do wprowadzenia nowych zmiennych, które są najczęściej liniowymi funkcjami starych zmiennych oraz podstawienia ich do równań i wykonania odpowiednich przekształceń. Dla równań NPP w jednym wymiarze przeskalowanie zmiennych i przejście do wielkości bezwymiarowych przeprowadzamy następująco ([16], [11]). Wprowadzamy współczynniki skalujące (których wartości będą ustalone w dalszej kolejności) xs , ts , Es , cs dla długości, czasu, pola elektrycznego i koncentracji. Formalnie zakładamy tylko, że są to dowolne liczby nieujemne. Jednakże dalsza przydatność przeskalowania zależy od wyboru konkretnych wartości, które ujawniają ważne charakterystyki układu. (Na przykład wybór wartości dla ts jest związany z występowaniem charakterystycznej skali czasowej – może ich być kilka. W tym kontekście często pojawiającym się wyborem dla ts jest wartość `2 /D, która jest typowym czasem dyfuzji w obszarze o rozmiarze `). Tak więc owe bezwymiarowe zmienne są zdefiniowane następująco (2.14). x¯ = x/xs , t¯ = t/ts , c¯i (¯ x, t¯) = ci (xs x¯, ts t¯)/cs , ¯ E(¯ x, t¯) = E(xs x¯, ts t¯)/Es .. Po podstawieniu tych zależności do równań układu NPP, wykonaniu rachunków i uporządkowaniu otrzymamy następujący układ równań bezwymiarowych (2.15). ∂¯ ci ∂ J¯i =− , ∂ t¯ ∂ x¯.

(20) 17. ROZDZIAŁ 2. WSTĘP I CEL PRACY. gdzie (2.16). ci ¯ i ∂¯ ¯ i zi c¯i E, ¯ +D J¯i = −D ∂ x¯. (2.17). X ∂ E¯ =λ zi c¯i , ∂ x¯ i=1. r. gdzie λ=. cs xs F . Es εr ε0. Równanie (2.7) ma teraz postać ∂ I¯ = 0, ∂ x¯. (2.18) gdzie. r. (2.19). ¯ X ci ¯ t¯) = 1 ∂ E − ¯ i ∂¯ I¯ = I( zi D + λ ∂ t¯ ∂ x ¯ i=1. Ã r X. ! ¯ i c¯i E, ¯ zi2 D. i=1. jest przeskalowaną gęstością całkowitego prądu niezależną od położenia x¯. Dla natężenia pola elektrycznego stosuje się zazwyczaj współczynnik skalu2. jący Es = FRT , a dla czasu ts = Dd s gdzie Ds charakterystyczna wielkość xs współczynników dyfuzji. Współczynnik λ ma wyraża się wtedy następująco λ=. cs F 2 x2s . RT εr ε0. Gdy wybierzemy skalowanie długości szerokością membrany d, tj. xs = d (wtedy bezwymiarowa szerokość wynosi 1 tzn. 0 ≤ x¯ ≤ 1), to dla typowych warunków eksperymentalnych T = 298[K], εr (H2 O) = 78, 54, 10−6 6 d 6 10−3 [m], 1 6 ci 6 103 [mol/m3 ], uzyskujemy następujący zakres zmienności parametru λ: (2.20). 10−12 6 λ−1 6 10−4 .. W związku z tym interesujące jest pytanie o sens fizyczny tak pojawiającej się wielkości bezwymiarowej. Okazuje się, że wartość ta jest związana z tzw. długością Debye’a [11], której wartość w rozważanym kontekście wynosi r RT εr ε0 . (2.21) LD = cs F 2.

(21) 18. ROZDZIAŁ 2. WSTĘP I CEL PRACY. Ze wzoru (2.21) oraz z definicji λ widać, że λ = (długość Debye’a)2 /(typowa długość). Ze względu na bardzo małą wartość λ−1 równanie Gaussa (2.4) sugeruje, że przestrzeń, w której odbywa się elektrodyfuzja może być podzielona na dwa obszary ([7], [62], [56]): jeden to roztwór w głębi, gdzie zachodzi warunek elektroobojętności (2.22). r X. zi ci = 0,. i=1. oraz drugi typ to warstwa brzegowa o typowej szerokości. √. λ−1 , w której kon-. centruje się nieskompensowany ładunek, zatem gdzie warunek elektroobojętności nie jest spełniony. To oznacza, że długość Debye’a LD określa w fizycznych jednostkach maksymalny rozmiar obszaru występowania ładunku przestrzennego w pobliżu granic membrany. Powyższa interpretacja jest zgodna z teorią warstwy podwójnej Guya-Chapmana [3], [56]. Można również dokonać interpretacji LD w oparciu o elementarne rozważania. Załóżmy, że ładunek o gęstości przestrzennej q zajmuje obszar o grubości dx. Z prawa Gaussa mamy dE ≈ λqdx. Z drugiej strony czysto dyfuzyjny składnik całkowitego strumienia jest proporcjonalny do DS dcs /dx. Aby ten czysto dyfuzyjny składnik zrównoważył składnik migracyjny Ds cs E musi być spełniony warunek (2.23). Ds cs λqdx ≈ Ds cs /dx, √ skąd otrzymujemy ponownie dx ≈ λ−1 . Podsumowując można stwierdzić, że warunek elektroobojętności jest pewnym przybliżeniem, które czasami może być stosowane a czasami nie. Jeżeli √ typowy rozmiar układu (np. szerokość membrany) jest rzędu λ−1 [m], to należy się spodziewać, że taka aproksymacja nie będzie dawać dobrych rezultatów, gdyż występuje wtedy istotne odstępstwo od warunku elektroobojętności. Zagadnienie to jest szczegółowo omówione na przykład w przeglądowej pracy Bazant, Thorton, Ajdari [5], gdzie wykorzystuje się metodę perturbacji singularnych,4 a także w pracy Urtenov, Kirillova, Seidova, Nikonenko [94], gdzie są wykorzystywane metody analityczne i numeryczne do 4. W przypadku metody perturbacji rozwiązanie poszukiwane jest w postaci szeregu.

(22) 19. ROZDZIAŁ 2. WSTĘP I CEL PRACY. analizy warunku elektroobojętności. W pracy Perram, Stiles [77] autorzy analizują założenie elektroobojętności w przypadku dwóch jonów monowalentnych. Ponadto próbują odpowiedzieć na pytanie o zakres stosowalności modelu Hendersona i modelu Goldmana (przypadek dwóch jonów). W tym celu wykonują symulacje pełnego układu NPP dla elektrolitu 1 : 1 w jednowymiarowej membranie. Do analizy używają także metod perturbacji, ale żmudne rachunki wykonują w środowisku Mathematica (dla układu NPP z czasem). W szczególności wnioskują, że przybliżenie Hendersona dobrze się sprawdza dla grubszych membran, a przybliżenie Goldmana dla cieńszych (rzędu nanometrów) – ale podkreślmy to jeszcze, analiza dotyczyła tylko przypadku układów 1 : 1.. 2.2. Zakres stosowalności modelu NPP. Układ równań NPP jako model elektrodyfuzji jest modelem fenomenologicznym. Posługuje się on takimi pojęciami jak gęstość substancji oraz wykorzystuje makroskopowy opis zajwisk elektrycznych.. Jak większość mod-. eli fizycznych jest to model przybliżony. Łączy on klasyczny opis dyfuzji wynikającej z gradientu koncentracji z wpływem pola elektrycznego, które jest wyliczane w oparciu o równania elektrodynamiki ośrodków ciągłych używającej średniej gęstości ładunków dyfundujących jonów.5 W związku z tym względem potęg małego parametru, który występuje w równaniach. Czasami szereg taki jest zbieżny globalnie (tzn. dla całego obszaru, w którym są określone równania). Mówimy wtedy o perturbacjach regularnych. Niestety w wielu przypadkach (tak jest dla równań NPP) mamy sytuację, w której występują np. dwa różne szeregi dla różnych podobszarów. Mówimy wtedy o perturbacji singularnej (osobliwej). Jest to trudniejszy przypadek od poprzedniego. Na ogół w tym przypadku występują tzw. warstwy brzegowe (ang. boundary layer) i warstwy zewnętrzne (ang. outer layer ) czyli obszary, w których są określone te rozwinięcia. Ponadto pojawia się także tzw. warstwa wewnętrzna (ang. inner layer ) – pośrednia pomiędzy poprzednimi, w której wspomniane rozwinięcia należy dopasować (ang. matching asymptotic expansions). 5 Wyrażenie na strumień można otrzymać z gradientu potencjału ´ ³ oczywiście ∂G , gdzie G jest energią swobodną Gibbsa chemicznego, i-tego składnika µi = ∂N i T,p,Nj6=i. układu, ale w żadnym stopniu nie wpływa to na prawdziwość powyższych stwierdzeń. W takim ujęciu [60], [68], gdy potencjał chemiczny µi nie jest stały względem położenia, to strumień Ji dąży do usunięcia tej różnicy. W liniowej teorii procesów nieodwracalnych, gdy zaniedbujemy efekty krzyżowe, strumień ten jest proporcjonalny do gradientu ∇µi (mówimy, że gradient potencjału chemicznego jest siłą napędową (driving force) dyfuzji): Ji = − Lτi ∇µi , gdzie Li współczynnik fenomenologiczny Onsagera (np. [60], rów. (10.3.16)). Dla fazy naładowanej elektrycznie potencjał chemiczny (elektrochemiczny) dla układów termodynamicznie idealnych ma postać µi = µoi + τ ln ci + eφ, więc Ji = − Lτi ( cτi ∇ci + e∇φ), co po zastosowaniu związku (10.3.7) [60], Lcii = Di oraz.

(23) 20. ROZDZIAŁ 2. WSTĘP I CEL PRACY. należy skomentować zakres stosowalności takiego modelu. Występują tu dwa zagadnienia dotyczące ogólnych warunków stosowalności równań NernstaPlancka-Poissona. Pierwsze zagadnienie dotyczy pomijania efektów pola magnetycznego. Ponieważ równania, którymi się posługujemy zawierają pole elektryczne i prądy (obie wielkości mogą zależeć od czasu), więc z równań Maxwella wynika możliwość pojawienia się pola magnetycznego B. W przypadku posługiwania się sformułowaniem w oparciu o potencjały oznacza to, że nie wystarczy rozważać tylko potencjału elektrycznego φ, ale należy używać także potencjału wektorowego (oznaczanego zwykle przez A). Wtedy pola E i B wyrażają się następująco E = −∇φ −. ∂A , ∂t. B = rotA.. Szczegóły można znaleźć w podręczniku [36], Rozdział 10. Dlaczego te efekty można pominąć? Powszechnie akceptowana odpowiedź sprowadza się do zauważenia, że zmienność czasowa rozważanych procesów jest po prostu mała. Największa charakterystyczna prędkość, które pojawia się w modelowanym układzie związana jest z relaksacją warstwy brzegowej (podwójnej warstwy elektrycznej – w terminologii elektrochemicznej). Oszacowanie tej prędkości jest następujące: LD /ts = Ds /LD ≈ 1, 0 [m/s] (jeden metr na sekundę). Wartość ta w porównaniu z prędkością światła (c = 3, 0 · 108 [m/s]) jest jak widać bardzo mała, co uzasadnia podejście, w którym nie są uwzględniane pole magnetyczne i ogólnie efekty falowe. Drugie zagadnienie dotyczy granic stosowalności modeli ciągłych zamiast opisu opartego na mechanice statystycznej, zwłaszcza, że w naszych rozważaniach pojawiają się skale przestrzenne rzędu angstremów (z (2.21) wynika, że LD może być rzędu 10−9 [m]). Problem ten dotyczy przede wszystkim małych układów, których rozmiary są porównywalne lub mniejsze od długości Debye’a. Takimi układami są na przykład kanały jonowe w błonach komórkowych. Innym modelem zjawisk elektrodyfuzji jonów może być opis wykorzystujący dynamikę Browna (Brownian dynamics, BD) [66], [72], [72]. Bardziej szczegółowymi modelami są te, które wykorzystują symulacje metodami dynamiki molekularnej, w której ruch xj = xj (t) każdej molekuły jonowej i wodnej jest opisywany przez drugą zasadę Newtona, mj x¨j (t) = Fj , E = −∇φ daje Ji = −Di ∇ci + τe Di ci E, czyli strumień Nernsta-Plancka (bez dryftu) (2.3)..

(24) 21. ROZDZIAŁ 2. WSTĘP I CEL PRACY. (j = 1, . . . , N ), gdzie N jest liczbą molekuł w układzie, a Fj całkowitą siłą działającą na molekułę j. Jednakże nawet dla kanałów jonowych symulacje takie są obecnie bardzo kosztowne czasowo i trudne do wykonania, chociaż w pewnych sytuacjach możliwe, np. [83], [65]. Model oparty o ruchy brownowskie jest pewnym kompromisem pomiędzy w pełni fenomenologicznym modelem ciągłym NPP a bardzo dokładnym opisem dynamiki molekularnej. W modelu tym tylko trajektorie indywidualnych jonów są obliczane w oparciu o równanie Langevin’a, które dla transportu jonów ma postać ([66], [65]) mj. d2 xj dxj = −mj γj + FR (t) + zj eEj , 2 dt dt. gdzie mj – masa, γj – współczynnik tarcia, zj – liczba ładunkowa,. dxj dt. = vj –. prędkość j-tej molekuły, e – ładunek elementarny. Natomiast otoczenie (np. woda) jest traktowane jako ośrodek ciągły, którego wpływ na ruch molekuły jest reprezentowany przez średnią siłę tarcia, −mj γj vj , oraz przez siłę losową FR (t) wynikającą ze zderzeń losowych. Ostatni wyraz w powyższym rówaniu, zj eEj , jest to całkowita siła elektryczna działająca na cząstkę j pochodząca od pola innych cząstek, ładunków zewnętrznych czy indukowanych. Pole elektryczne jest wyliczane w oparciu o prawo Gaussa (gdy używamy potencjału elektrycznego zamiast wprost natężenia pola elektrycznego, to jest to równanie Poissona). Jedną z metod weryfikacji modelu NPP jest porównanie jego przewidywań różnych wielkości fizycznych (np. koncentracje, strumienie składników, potencjał elektryczny) z wartościami otrzymanymi w symulacjach dynamiki Browna. Na uwagę zasługują tutaj opublikowane w 2000r. prace grupy Moy, Corry, Kuyucak i Chung [72], [73], [15], w których autorzy podjęli zagadnienie przetestowania modelu NPP w odniesieniu do kanałów jonowych. Były to kanały cylindryczne o zmiennym promieniu (3Å-14Å), długościach rzędu 40Å oraz z ustalonymi lub bez ustalonych ładunków. Długość Debye’a w tych układach była rzędu 2Å. Generalny wniosek jest taki, że w przypadku układów, których promienie były porównywalne lub mniejsze od długości Debye’a, model NPP się załamuje i daje wyniki odbiegające od bardziej dokładnych modeli..

(25) 22. ROZDZIAŁ 2. WSTĘP I CEL PRACY. 2.3 2.3.1. Przegląd numerycznych metod w NPP. Równania NPP niezależne od czasu. Historia metod numerycznych modelu NPP zaczęła się właściwie z chwilą pojawienia się samego modelu. Już w 1890 roku Planck w swojej pracy [78] z 1890 roku opisał pewien sposób uzyskiwania wartości numerycznej potencjału dyfuzyjnego. Praca ta była przede wszystkim pierwszym znaczącym wkładem w teorię elektrodyfuzji. Przypadek tam rozważany dotyczył stacjonarnego stanu dyfundujących monowalentnych jonów (tzn. o liczbach ładunkowych +1 lub -1). Warunki brzegowe były typu Dirichleta – zadana wartość koncentracji oraz potencjału elektrycznego na obu brzegach. Ponadto wykorzystywany był warunek elektroobojętności. Przy takich założeniach udało się Planckowi wzór na potencjał dyfuzyjny, który jest wyrażony następująco (2.24). τ ∆φM = − ln ξ, e. gdzie pomocnicza zmienna ξ spełnia nieliniowe równanie algebraiczne jednej zmiennej (3.35). W równaniu tym występują parametry opisujące układ (koncentracje na brzegach oraz ruchliwości składników). Równanie to jest dość skomplikowane, ale nie nastręcza żadnych trudności przy numerycznym rozwiązywaniu. Planck opisał także graficzną metodę takiego rozwiązywania. Dlatego można jego pracę uznać za początek metod numerycznych w zagadnieniu NPP. Równanie podane przez Plancka można podać w nieco innej formie, którą zaproponował 90 lat później Morf [71] wraz z nowym, prostszym wyprowadzeniem. Ponadto wyprowadzenie to dotyczy nieco ogólniejszego przypadku, w którym są dwie grupy jonów: dodatnie jony mają ten sam ładunek z+ > 0, ujemne mają ten sam ładunek z− < 0. (Model Plancka odpowiada sytuacji z+ = +1 oraz z− = −1). W tym ujęciu równanie na potencjał dyfuzyjny jest postaci (2.25). ¯ + (∆φM ) − D ¯ − (∆φM ) c(d) D ∆φM = − ¯ ¯ − (∆φM ) ln c(0) , D+ (∆φM ) + D. ¯ ± (∆φM ) są określone następująco gdzie funkcje D  P P z+ ∆φ ¯ + (∆φM ) = ( pPDp cp (d))ez ∆φM −Pp Dp cp (0) ,  D ( cp (d))e + M − p cp (0) P P p (2.26) z ∆φ ( D n  D ¯ − (∆φM ) = nP cn (d))ez− ∆φM −Pn Dn cn (0) , c (0) ( c (d))e − M − n n. n n.

(26) ROZDZIAŁ 2. WSTĘP I CEL PRACY. gdzie symbol. 23. P. oznacza sumowanie po jonach dodatnich o liczbie ładunP kowej z+ > 0, symbol n oznacza sumowanie po jonach ujemnych o liczbie ładunkowej z− < 0, oraz P P cn (d) c(d) p cp (d) =P = Pn . c(0) p cp (0) n cn (0) p. (Równość powyższych ilorazów sum po dodatnich (cp ) i ujemnych (cn ) jonach wynikaja z przyjętych w modelu założeń: warunku elektroobojętności oraz występowania jonów tylko o ładunkach z+ lub z− ).6 Ponieważ równanie (2.25) ma postać równania na punkt stały7 ∆φM = f (∆φM ), więc naturalna metoda numeryczna, którą stosuje Morf [71] polega na prostych iteracjach ∆φM,k+1 = f (∆φM,k ) k = 0, 1, 2 . . . . Takie podejście wymaga podania wartości początkowej ∆φM,0 , która jak pisze Morf może być dowolna (ale w przedstawionych symulacjach zawsze startował od ∆φM,0 = 0). Liczba iteracji, która wystarczała do otrzymania potencjału z dokładnością do ±0, 01 mV była w zakresie od 2 do 27 (dla rozważanych układów KCL, NaCl, HCl i NaOH). W odniesieniu do metody prostych iteracji należy podnieść jednakże jedno zastrzeżenie. Otóż stwierdzenie, iż punkt startu ∆φM,0 można wybrać dowolnie nie musi być prawdziwe. Zależy to bowiem od własności funkcji f = f (∆φM ). Gdy spełnia ona warunek „zwężania”, czyli, gdy istnieje stała 0 ≤ L < 1 taka, że |f (φ1 ) − f (φ2 )| ≤ L|φ1 − φ2 |, to faktycznie jako punkt startu można wybrać dowolną wartość ∆φM,0 . Wynika to z twierdzenia Banacha o punkcie stałym [30]. Ale nie jest wiadome P Z warunku elektroobojętności mamy dla x ∈ [0, d] : zgrupowai zi ci (x) = 0, co poP niuPwzględem jonów dodatnich (z = z ) i ujemnych (z = z ) daje z i +P i − + p cp (x) + P z− n cn (x) = 0, P czyli z+ p cp (x)P = −z− n cnP (x). Podstawiając do tej równości x=0 P i x = d mamy z+ P c (d) = −z c (d), z c (0) = −z c (0), co po podzie− p p n n P− n n P + p pP leniu stronami daje p cp (d)/ p cp (0) = n cn (d)/ n cn (0). 7 Jeżeli dana jest funkcja f : X → X, to punktem stałym tej funkcji nazywamy dowolny element x ∈ X taki, że f (x) = x. Teoria punktów stałych odwzorowań (jak i ich miejsc zerowych) oraz sposoby ich wyznaczania są ważnymi działami zarówno matematyki jak i metod numerycznych. 6.

(27) ROZDZIAŁ 2. WSTĘP I CEL PRACY. 24. czy podana w pracy Morfa funkcja spełnia ten warunek. Funkcja ta jest określona przez wyrażenia z (2.25) i (2.26) i jest ona skomplikowana, tak więc sprawdzenie globalnego warunku „zwężania” może być trudne.8 Z drugiej strony postępowanie praktyczne polega na startowaniu od wartości ∆φM,0 (najlepiej jak najbliżej spodziewanej wartości potencjału) i sprawdzaniu czy iteracje się stabilizują. Jeżeli tak, to oznacza że otrzymaliśmy potencjał membranowy ∆φM . Jeżeli nie, to należy spróbować rozpocząć iteracje od innej wartości początkowej ∆φM,0 . Na koniec omawiania modelu Plancka należy jeszcze podkreślić, że oba opisane rozwiązania dotyczą jedynie wyznaczania potencjału dyfuzyjnego ∆φM = φ(d− ) − φ(0+ ). Metody te nie dają bezpośrednio informacji o profilu potencjału, φ = φ(x) oraz profilów koncentracji, ci = ci (x). W 1943 roku Goldman w pracy [35] – opisującej m. in. impedancję membran – podał rozwiązanie w układzie stacjonarnym, które nie wykorzystywało założenia elektroobojętności. W zamian za to przyjął on inne założenie – stałość pola elektrycznego w membranie (a zatem liniowość potencjału elektrycznego względem położenia). Założenie to pozwoliło na analityczne scałkowanie równań i podanie wyrażenia na różnicę potencjału elektrycznego pomiędzy granicami membrany. W latach 1949-1952 Hodgkin, Huxley i Katz wykorzystali idee sformułowane przez Goldmana do opisu propagacji potencjału membranowego w aksonie kałamarnicy9 [45]. Samo wyznaczenie potencjału w oparciu o uproszczenie Goldmana nie 8. Autor rozprawy wykonał taką analizę stosując pomocnicze podstawienia, które sprowadzają funkcję f do postaci wymiernej – dokładniej do ilorazu dwóch funkcji kwadratowych. Współczynniki tych funkcji zależą m. in. od koncentracji na brzegach i od współczynników dyfuzji. Uzyskane w ten sposób wyrażenie na L daje warunek „zwężanie”, jest on jednak skomplikowany i w praktyce trudno go stosować do ogólnych rozważań. Ale dla danych konkretnych wartości ci (d), ci (0), Di możne je podstawić i otrzymać wartość liczbową L. Dla wspomnianych obliczeń w pracy Morfa [71] stała L jest rzędu 0, 72 (czyli mniejsza od 1), co potwierdza, że iteracje powinny być zbieżne do rozwiązania. 9 Akson (neuryt) to wypustka komórki nerwowej (neuronu), ma postać długiego cienkiego włókna – może mieć długość od kilku mikrometrów do 1m oraz średnicę rzędu 1µm – przez które przewodzone są sygnały nerwowe. Sygnały te są natury elektrochemicznej i mają postać impulsu potencjału elektrycznego, który wędruje wzdłuż zewnętrznej membrany aksonu. Najważniejsze jony, które biorą udział w powstawaniu potencjału to kationy N a+ , K + , Ca2+ i aniony Cl− . Aksony są podstawowymi liniami transmisyjnymi układu nerwowego, a ich wiązki budują nerwy. W 1952r. Alan Lloyd Hodgkin i Andrew Huxley opisali model wyjaśniający jonowy mechanizm generujący inicjalizację i propagację potencjału elektrycznego w aksonie kałamarnicy. Otrzymali za to w 1963r. Nagrodę Nobla w dziedzinie fizjologii i medycyny..

(28) 25. ROZDZIAŁ 2. WSTĘP I CEL PRACY. stanowi praktycznie żadnego problemu numerycznego. W przypadku, który omawia on w swej pracy (wszystkie jony są monowalentne, tzn. +1 lub −1 jest on dany w sposób jawny. Na przykład dla trzech jonów o liczbach ładunkowych z1 = z2 = +1, z3 = −1) wzór na potencjał ma postać10 (2.27). τ D1 c1 (d) + D2 c2 (d) + D3 c3 (0) ∆φM = − ln . e D1 c1 (0) + D2 c2 (0) + D3 c3 (d). Warto tutaj nadmienić, że w ogólnym przypadku (dowolna liczba składników r, dowolne liczby ładunkowe zi dla i = 1, . . . , r) założenie o stałości natężenia pola elektrycznego prowadzi do nieliniowego równania algebraicznego na potencjał ∆φM . Autor wyprowadził takie równanie przy założeniach Goldmana (szczegóły w Dodatku), tutaj podany jest tylko rezultat r X ai − bi ezi ∆φM i=1. 1 − ezi ∆φM. = 0,. gdzie ai = zi2 Di ci (d), bi = zi2 Di ci (0). W Dodatku pokazano, że z powyższego równania wynika (2.27), gdy zi ∈ {−1, +1}. Ponadto są tam pokazane inne przypadki, które prowadzą do wzorów analitycznych na potencjał ∆φM . Tak więc mamy równanie z jedną niewiadomą postaci f (∆φM ) = 0. Do znajdowania miejsca zerowego mamy do dyspozycji cały arsenał metod numerycznych (metoda Newtona, metoda bisekcji, metoda siecznych itd.) [80], [55]. Należy jeszcze zaznaczyć, że uproszczony model Goldmana pozwala na dokładne scałkowanie stacjonarnego układu NPP, gdy są spełnione pewne szczególne warunki dotyczące koncentracji jonów w otaczającym membranę roztworze. Warunki takie podali Arndt, Bond i Roper [1] (1970): (a) wszystkie jony dodatnie mają tę samą liczbę ładunkową z+ , a wszystkie ujemne tę samą liczbę ładunkową z− ; (b) na brzegu zachodzi warunek elektroobojętności; (c) całkowita koncentracja jonów dodatnich dla x = 0 jest taka sama jak dla x = d (to samo będzie dla jonów ujemnych – na mocy (a) i (b)). Oznacza to, że w tym przypadku dysponujemy dokładnymi wyrażeniami analitycznymi na profile koncentracji. Pozwala to stosować ten przypadek do testowania metod numerycznych rozwiązujących stacjonarny układ NPP, gdy spełnione są powyższe warunki. 10. Ściśle rzecz biorąc tradycyjne równanie Goldmana ma nieco inną postać, w którym zamiast współczynników dyfuzji Di występują tzw. współczynniki przepuszczalności błony, Pi , (np. [62])..

(29) ROZDZIAŁ 2. WSTĘP I CEL PRACY. 26. Kolejnym krokiem były przede wszystkim prace Schlögla [87] (1954), w których uogólnił on wyniki Plancka na przypadek jonów o dowolnych liczbach ładunkowych. Pozostałe założenie były takie same jak u Plancka.. 2.3.2. Równania NPP zależne od czasu. Pierwsze głębsze analizy zależnych od czasu równań NPP (stan niestacjonarny) pojawiły się w 1965.. Co ciekawe, dokładnie w tym samym nu-. merze 5 czasopisma Biophysical Journal ukazały się dwie prace dotyczące tego zagadnienia. Były to: Cohen, Cooley The Numerical Solution of the Time-Dependent Nernst-Planck Equations [16] oraz Conti, Eisenman The Non-Steady-State Membrane Potential of Ion Exchangers with Fixed Sites [17]. Jednakże to praca Cohen-Cooley zasługuje tutaj na szczególne wyeksponowanie, gdyż była to pierwsza próba zastosowania na szeroką skalę metod numerycznych z wykorzystaniem komputerów do zależnych od czasu równań NPP. W publikacji tej autorzy zaproponowali – jak się wydaje – po raz pierwszy tak rozbudowaną procedurę numeryczną, której implementacja wymagała użycia programu i komputera. Algorytm numeryczny, którym się posłużyli opiera się na różnicach skończonych na stałej siatce przestrzennej i zmiennej siatce czasowej.11 Wbrew opinii spopularyzowanej przez Brumleve, Buck [11],12 procedura numeryczna całkowania równań NPP zrealizowana przez Cohena i Cooley’a nie jest schematem jawnym!13 Jednakże w istocie autorzy ci stosują pewien wariant schematu niejawnego przy obliczaniu wartości koncentracji i pola elektrycznego w kolejnym momencie tk+1 na podstawie znanych wartości w momencie tk . Szczegóły procedury są jednak opisane bezpośrednio, bez wykorzystania powszechnie używanej terminologii – stąd być może wspomniana błędna opinia. W 1975 Sandifer i Buck [85] użyli mieszanej metody do rozwiązyania 11. Węzły przestrzenne {xi } spełniają warunek xi+1 −xi = h. Wielkość h > 0 jest zadana na początku obliczeń i nie ulega modyfikacjom w trakcie wykonywania się procedury. Natomiast krok czasowego tk+1 = tk + δt może się zmieniać w zależności od wyliczanego pomocniczo błędu. Nie jest to jednak w pełni metoda adaptacyjna, gdyż nie następuje zwiększanie kroku czasowego w kolejnym etapie, gdy błąd jest mały. 12 str. 2 „In 1965, Cohen and Cooley developed a simulation procedure for the complete diffusion-migration problem. Although their procedure was explicite (non-predictive) (...)” (podkreślenie moje). 13 Wyjaśnienie różnicy pomiędzy schematem jawnym (ang. explicite) a niejawnym (ang. implicite) jest w dalszej części rozdziału..

(30) 27. ROZDZIAŁ 2. WSTĘP I CEL PRACY. układu NPP: niejawnej do obliczania pola elektrycznego, jawnej dla obliczania koncentracji. Ponieważ użyto numerycznych schematów jawnych, więc spowodowało to naturalne dla tych metod ograniczenie zredukowanego współczynnika dyfuzji przez 1/2.14 To ograniczenie staje się szczególnie uciążliwe dla układów o dużej grubości, gdyż wtedy ze względu na duże gradienty przy brzegach krok przestrzenny musi być mały, a w konsekwencji krok czasowy bardzo mały (∆t ∼ (∆x)2 , patrz przypis). Należy także dodać, że praca Cohen-Cooley była nowatorska pod jednym jeszcze względem. Dotyczy to analizy założenia elektroobojętności i zastosowania w tym celu metody perturbacji (w istocie była to perturbacja singularna) [54]. Jakkolwiek temat ten nie został w pełni rozwinięty, to jednak został on wprowadzony do literatury z zakresu chemii fizycznej. W szczególności autorzy pokazali, że w stanie stacjonarnym elektroobojętność zachodzi w przypadku, gdy koncentracje są równe na obu brzegach, ale gdy są różne, to założenie to może być słuszne we wnętrzu membrany z wyjątkiem cienkich warstw w pobliżu brzegu. W roku 1978 ukazała się praca Brumleve’a i Bucka [11], która była istotnym postępem w dziedzinie obliczeń numerycznych układu NPP i ich zastosowań. Po pierwsze rozważany problem dotyczył stanów niestacjonarnych (stan stacjonarny pojawia się również, ale jako przypadek graniczny do którego dąży układ (t → ∞)). Po drugie autorzy nie używali założenia elektroobojętności, tylko pełne równanie Gaussa (w modelu numerycznym autorzy preferowali zapis równań z natężeniem pola elektrycznego zamiast potencjału elektrycznego). Technicznie rzecz biorąc równanie Gaussa zostało zastąpione równaniem na pole elektryczne z prądem przesunięcia (co już było wykorzystane w pracy Cohena i Cooley’a [16]). Po trzecie użyto warunków brzegowych typu Neumanna, w których strumienie na brzegu zależą od kon14. Zjawisko to można prześledzić na uproszczonym modelowym problemie czystej dyfuzji ∂2c pojedynczego składnika ∂c ∂t = D ∂x2 . Obliczanie jawne przybliżonej wartości koncentracji k+1. k. k. k. k. )−c (xi ) w chwili tk+1 na podstawie chwili tk ma postać c (xi∆t = D c (xi+1 )−2c∆x(x2i )+c (xi−1 ) , ¯ k (xi+1 ) + (1 − 2D)c ¯ k (xi ) + Dc ¯ k (xi+1 ), gdzie zreco można przepisać tak: ck+1 (xi ) = Dc 2 ¯ dukowany współczynnik D = D∆t/∆x . Wyliczanie iteracyjne kolejnych wartości ck ¯ > 0, czyli D ¯ < 1/2 (np. [80], str. wg podanego wzoru nie jest rozbieżne, gdy 1 − 2D 585). Oznacza to, że dla stabilności metody musi być spełniony warunek ∆t ∼ (∆x)2 co wymusza bardo mały krok czasowy..

(31) ROZDZIAŁ 2. WSTĘP I CEL PRACY. 28. centracji składników w membranie i na zewnątrz następująco  ←  Ji (0, t) = ~kiL ciL − k iL ci (0, t), (2.28)  J (d, t) = −~k c + ← k iR ci (d, t), i iR iR ←. ←. gdzie współczynniki k iL , ~kiL , k iR , ~kiR to stałe heterogeniczne szybkości reakcji lub przenikania na brzegach, ciL , ciR to koncentracje na zewnątrz membrany (odpowiednio z lewej i prawej strony). Warunki takie są bardziej fizyczne, a co istotne są ogólniejsze niż proste warunki typu Dirichleta narzucające wartości koncentracji na brzegach (w pracy z 2006r. Perram, Stiles [77] autorzy wręcz stwierdzają, że warunki typu Dirichleta są „eksperymentem myślowym”), (2.29). ci (0, t) = ciL , ci (d, t) = ciR .. Co więcej, na warunki brzegowe Dirichleta można patrzeć jak na graniczny przypadek poprzednich warunków: gdy stałe heterogeniczne ki stają się bardzo duże (ki → ∞), to warunki (2.28) przechodzą w warunki Dirichleta. Po czwarte wreszcie, autorzy zaprezentowali efektywny algorytm numeryczny do całkowania tak ogólnego modelu i przedstawili szereg wyników swoich symulacji. Procedura numeryczna, która wykorzystali jest niejawna (implicite), opiera się o metodę różnic skończonych z pełną jednoczesną dyskretyzacją przestrzenną i czasową. Dzięki temu, że schemat jest niejawny nie ma ograniczeń na zredukowany współczynnik dyfuzji, a w konsekwencji wymogu używania bardzo małego kroku czasowego ∆t ∼ (∆x)2 (tak jak to miało miejsce w pracy Sandifer, Buck [85]). Metody zaprezentowane w pracy Brumleve’a-Bucka były później wielokrotnie wykorzystywane przez innych autorów w szczególności do interpretowania potencjałów i koncentracji w membranach jonoselektywnych (Sokalski, Lewenstam [89]) czy do badania widm impedancyjnych układów elektrochemicznych (Danielewski, Kucza, Lewenstam [26], [61]).. 2.3.3. Stosowane metody numeryczne. Wszystkie omówione powyżej algorytmy numeryczne bazowały na metodach różnicowych. Oznacza to, że procedurę obliczeniową uzyskuje się poprzez.

(32) 29. ROZDZIAŁ 2. WSTĘP I CEL PRACY. zastąpienie pochodnych odpowiednimi różnicami skończonymi. Występuje tutaj cała gama możliwości. Można stosować różny rząd aproksymacji pochodnych (w pierwszym rzędzie zależy to od liczby punktów, w których przeprowadza się wartościowanie funkcji), np. wyrażenia f 0 (x) = f 0 (x) =. 1 (f (x + h) − f (x − h)) + 2h 1 (f (x − 2h) − 8f (x − h) 12h. O(h2 ), + 8f (x + h) − f (x + 2h)) + O(h4 ),. przedstawiają aproksymacje centralne odpowiednio dwu- i cztero-punktowe dla siatki jednorodnej (stała odległość pomiędzy węzłami równa h > 0). Zauważmy, że drugi schemat jest dokładniejszy w porównaniu z pierwszym (błąd jest rzędu h4 a w pierwszym przypadku h2 ) ale wymaga więcej wartościowań funkcji f = f (x), tak więc wydłuż czas obliczeń.Więcej przykładów różnic skończonych jest w Dodatku B. Oprócz liczby punktów używanych do przybliżania pochodnych (w powyższym przykładzie były to dwa lub cztery punkty) można jeszcze modyfikować odległość pomiędzy nimi uzyskując w ogólności schematy ze stałą siatką (tzn. xi+1 − xi = h = const) lub schematy ze zmienną siatką (tzn. xi+1 − xi = hi ). Ten drugi sposób daje większą elastyczność przy aproksymacji równań, ale wymaga pewnej wstępnej wiedzy pozwalającej zagęszczać siatkę tam gdzie występuję większa zmienność koncentracji lub potencjału (czyli, gdzie występują duże gradienty). W przypadku membran jonoselektywnych takie duże gradienty mogą wystąpić przede wszystkim w pobliżu brzegów. Jak już to zasygnalizowano przy omawianiu prac [16] i [11], schematy różnicowe mogą być również podzielone na tzw. jawne i niejawne (ang. explicite oraz implicite). W przypadku metod jawnych uzyskanie wartości w kolejnym kroku obliczeń nie wymaga rozwiązywanie dodatkowego równania (lub układu równań) – wystarczy po prostu zwykłe wartościowanie (obliczenie wartości funkcji). Dla metod niejawnych uzyskanie wartości w kolejnym kroku wymaga rozwiązania dodatkowo pomocniczego równania, a zatem koszt takiej metody jest większy. Różnicę tę dobrze ilustruje najprostszy przykład – całkowanie pojedynczego równania różniczkowego zwyczajnego metodą Eulera. Rozważmy klasyczny problem początkowy dla układu równań różniczkowych zwyczajnych, w którym niewiadomymi są funkcje xi (t), określony.

(33) 30. ROZDZIAŁ 2. WSTĘP I CEL PRACY. następująco ½ (2.30). xi 0 (t) = fi (t, x1 (t), . . . , xn (t)), xi (t0 ) = xi,0 ,. (i = 1, . . . , n). gdzie fi : Rn+1 → R są danymi funkcjami. W bardziej zwartej notacji wektorowej, gdzie x = (x1 , . . . , xn ) ∈ Rn , f = (f1 , . . . , fn ), układ taki można zapisać w skrócie tak ½ (2.31). x0 (t) = f (t, x(t)), x(t0 ) = x0 .. Zadaniem praktycznie każdej metody numerycznej całkowania układu równań różniczkowych zwyczajnych jest podanie ciągu wektorów xk ∈ Rn (w przypadku pojedynczego równania jest to ciąg liczb xk ∈ R), które przybliżają rozwiązanie dokładne x = x(t) w wybranych punktach {tk }, czyli xk ' x(tk ) dla k = 0, 1, 2, . . .. W przypadku jawnej metody Eulera ciąg przybliżeń xk jest określony następująco: ½ 0 x = x0 , (2.32) xk+1 = xk + hf (tk , xk ) (dla k = 1, 2, . . .). Podstawowe znaczenie ma tutaj sposób przejścia od przybliżenia xk do xk+1 . Otóż przejście od kroku k do k + 1 odbywa się następująco: xk+1 = xk + hf (tk , xk ) – tak więc obliczanie kolejnej wartość przybliżenia polega na zwykłym wartościowaniu funkcji, f (tk , xk ).. Dlatego używane jest określenie. metoda jawna, gdyż sama postać równości (2.32) daje natychmiast kolejną wartość przybliżenia. Natomiast w przypadku metod niejawnych – np. niejawnej metody Eulera – przejście od kroku k do k + 1 wymaga więcej pracy i wiąże się z rozwiązaniem pewnego dodatkowego problemu. W przypadku wspomnianej niejawnej metody Eulera zależność ta wygląda następująco: (2.33). xk+1 = xk + hf (tk+1 , xk+1 ).. Widać teraz podstawową różnicę w porównaniu z metodą jawną (2.32): szukana wartość xk+1 występuje też po prawej stronie równości. Tak więc, aby otrzymać kolejne przybliżenie, xk+1 , musimy rozwiązać układ równań (na ogół nieliniowy) względem x o postaci (2.34). x = xk + hf (tk+1 , x)..

(34) ROZDZIAŁ 2. WSTĘP I CEL PRACY. 31. W tym przypadku uzyskanie kolejnej wartości przybliżenia xk+1 nie jest takie proste, jak w przypadku jawnej metody Eulera (2.32). Tym razem musimy rozwiązać dodatkowy problem, który dla niejawnej metody Eulera ma postać powyższego (2.34) układu algebraicznego. Układ taki może być rozwiązywany różnymi metodami. Najczęściej stosuje się procedurę prostych iteracji lub procedury bazujące na metodzie Newtona-Raphsona. Metoda ta dla układu n równań postaci F (x) = 0 (x ∈ Rn ) polega na obliczaniu kolejnych przybliżeń xj (2.35). xj+1 = xj − [DF (xj )]−1 · F (xj ) (dla j = 1, 2, . . .),. które powinny zbliżać się do rozwiązania (limj→∞ xj = x). We wzorze tym: DF (xj ) macierz pochodnej funkcji F w punkcie xj , wykładnik (−1 ) oznacza odwracanie macierzy, a kropka (·) mnożenie macierzy.15 Metoda jest bardzo szybko zbieżna, ale trzeba pamiętać, że wymaga dodatkowo obliczenia macierzy pochodnych oraz wybrania jako punktu startowego do iteracji (2.35), punktu dostatecznie bliskiego rozwiązaniu ([80], str. 281-287). W przypadku metody elementów skończonych (MES, ang. FEM, Finite Elements Method ) proces uzyskiwania algorytmu numerycznego jest odmienny od stosowanego w metodach różnicowych. Będzie to szczegółowiej opisane w rozdziale poświęconym MES-owi w problemie NPP. Tutaj podamy tylko kilka uwag ogólnych. Po pierwsze równanie różniczkowe (zwyczajne lub cząstkowe) należy zapisać w innej, ogólniejszej postaci. Postać tę wyprowadza się z pewnych tożsamości całkowych (w szczególności istotną rolę odgrywa tutaj wzór na całkowanie przez części). Postać uogólniona (zwana także postacią słabą lub wariacyjną) jest istotnie ogólniejsza: każde rozwiązanie klasyczne jest rozwiązaniem uogólnionym, ale nie na odwrót. Po drugie dyskretyzacja numeryczna w MES opiera się właśnie na uogólnionym sformułowaniu i uzyskana jest w ten sposób, że konstruuje się przestrzeń funkcyjną, która w pewnym sensie aproksymuje przestrzeń, w której znaj15. W praktyce wzór w postaci (2.35), która zawiera odwracanie macierzy DF jest rzadko stosowany – z wyjątkiem przypadku jednowymiarowego, gdy mamy pojedyncze równanie j ) F (x) = 0, (x ∈ R), wtedy ma on postać xj+1 = xj − FF0(x (xj ) (metoda stycznych). Dla układu równań lepiej jest równość (2.35) zapisać tak: DF (xj ) · (xj+1 − xj ) = −F (xj ). Jest to układ równań liniowych z macierzą współczynników DF (xj ) i prawą stroną −F (xj ). Zatem algorytm ma postać: (1) wybieramy punkt startu x0 ; (2) w pętli wykonujemy: (i) rozwiązywanie układu równań względem z: DF (xj )z = −F (xj ), (ii) podstawienie: xj+1 = xj + z..