Stosowane metody numeryczne

  J_i(0, t) = ~k_iLc_iL−^←k_iLc_i(0, t), J_i(d, t) = −~k_iRc_iR+^←k_iRc_i(d, t),

gdzie współczynniki^←kiL, ~kiL,^←kiR, ~kiRto stałe heterogeniczne szybkości reakcji lub przenikania na brzegach, c_iL, c_iR 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) c_i(0, t) = c_iL, c_i(d, t) = c_iR.

Co więcej, na warunki brzegowe Dirichleta można patrzeć jak na graniczny przypadek poprzednich warunków: gdy stałe heterogeniczne k_i stają się bar-dzo duże (k_i → ∞), 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

(im-plicite), opiera się o metodę różnic skończonych z pełną jednoczesną

dyskre-tyzacją 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 wielokrot-nie wykorzystywane przez innych autorów w szczególności do interpretowa-nia potencjałów i koncentracji w membranach jonoselektywnych (Sokalski, Lewenstam [89]) czy do badania widm impedancyjnych układów elektro-chemicznych (Danielewski, Kucza, Lewenstam [26], [61]).

2.3.3 Stosowane metody numeryczne

Wszystkie omówione powyżej algorytmy numeryczne bazowały na metodach

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

f0(x) = 1

2h(f (x + h) − f (x − h)) + O(h2),

f0(x) = 1

12h(f (x − 2h) − 8f (x − h) + 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). Za-uważ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.

x_i+1− x_i = h = const) lub schematy ze zmienną siatką (tzn. x_i+1− x_i = h_i). 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 wys-tę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.

ex-plicite oraz imex-plicite). 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 (oblicze-nie wartości funkcji). Dla metod (oblicze-niejawnych uzyska(oblicze-nie wartości w kolej-nym kroku wymaga rozwiązania dodatkowo pomocniczego równania, a za-tem 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óżnicz-kowych zwyczajnych, w którym niewiadomymi są funkcje xi(t), określony

następująco (2.30)

½

xi⁰(t) = fi(t, x1(t), . . . , xn(t)), (i = 1, . . . , n)

x_i(t₀) = x_i,0,

gdzie f_i : Rn+1 → R są danymi funkcjami. W bardziej zwartej notacji

wektorowej, gdzie x = (x₁, . . . , x_n) ∈ Rn, f = (f₁, . . . , f_n), układ taki można zapisać w skrócie tak

(2.31)

½

x0(t) = f (t, x(t)),

x(t₀) = x₀.

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

przy-bliżają rozwiązanie dokładne x = x(t) w wybranych punktach {t_k}, czyli xk ' x(t_k) dla k = 0, 1, 2, . . .. W przypadku jawnej metody Eulera ciąg przybliżeń xk jest określony następująco:

(2.32)

½

x0 = x₀,

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 (t_k, xk) – tak więc obliczanie kolejnej wartość przybliżenia polega na zwy-kłym wartościowaniu funkcji, f (t_k, 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. nie-jawnej 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 (t_k+1, xk+1).

Widać teraz podstawową różnicę w porównaniu z metodą jawną (2.32): szu-kana 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

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) x^j+1 = x^j − [DF (x^j)]⁻¹· F (x^j) (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) ozna-cza odwracanie macierzy, a kropka (·) mnożenie macierzy.15 Metoda jest bardzo szybko zbieżna, ale trzeba pamiętać, że wymaga dodatkowo oblicze-nia macierzy pochodnych oraz wybraoblicze-nia 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

odmi-enny od stosowanego w metodach różnicowych. Będzie to szczegółowiej opisane w rozdziale poświęconym MES-owi w problemie NPP. Tutaj po-damy tylko kilka uwag ogólnych. Po pierwsze równanie różniczkowe (zwycza-jne 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

znaj-15W 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 F (x) = 0, (x ∈ R), wtedy ma on postać xj+1 = xj −_F^{F (x}0(x^jj⁾) (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.

dują się rozwiązania problemu. Następnie próbuje się zapisać wyjściowe rów-nanie w tej przybliżonej przestrzeni. W praktyce obliczeniowej oznacza to, że przybliżonego rozwiązania u_h16 poszukuje się w postaci sumy względem specjalnie dobranych funkcji bazowych {ϕ_k}

(2.36) uh(x, t) =

n

X

k=1

λk(t)ϕk(x).

We wzorze tym występują jako nieznane funkcje λ_k(t). Zadaniem metody elementów skończonych jest określenie sposobu wyznaczania tych funkcji. Prowadzi to do układu równań różniczkowych zwyczajnych na funkcje {λ_k(t)}. W niniejszej pracy będą m. in. wyprowadzone takie wzory dla układu NPP w Rozdziale 4 i użyte do obliczenia koncentracji składników oraz potencjału w membranie.

Bardzo istotnym elementem w MES jest właśnie dobór funkcji bazowych

{ϕk}. Jest on zawsze związany z pewnym podziałem obszaru

przestrzen-nego, w którym określone są równania. Podział ten zwany jest triangulacją, gdyż w przypadku dwuwymiarowym jest on oparty o trójkąty (w przypadku trójwymiarowym o czworościany).

Dominującymi metodami numerycznymi wykorzystywanymi dotychczas w elektrochemii są metody oparte o różnice skończone. Świadczy o tym chociażby trzecie wydanie książki ”Digital Simulation in Electrochemistry” Dietera Britza [10], w którym nie ma żadnych praktycznie żadnych informacji o MES17 oraz przegląd dostępnej literatury naukowej.

Należy również podkreślić, że w pracach opisujących transport ładunków w półprzewodnikach – gdzie używane równania są podobne do tych z elek-trochemii – dominującymi metodami numerycznymi są także metody oparte o różnice skończone. Najczęściej bazują one na różnych wariantach pochodzą-cych z pracy Scharfetter, Gummel [86]. Oryginalna metoda Gummela-Sharfetera wykorzystuje schemat Cranka-Nicolsona do kroku czasowego.18 Szczegółowy

16Indeks h w uhma podkreślić, że jest to przybliżenie, które zależy od stopnia aproksy-macji przestrzeni rozwiązań przestrzenią rozwiązań przybliżonych. W praktyce taką przestrzeń konstruuje się w oparciu o wielościany, których średnice nie przekraczają h. (Dla jednego wymiaru wielościany to odcinki).

17Dokładniej: nie było żadnej wzmianki w dwóch wcześniejszych wydaniach książki [10], natomiast w najnowszym jest tylko krótka informacja na pół strony.

18Metoda Cranka-Nicolsona jest kolejnym przykładem metody niejawnej. Jest ona często używana do rozwiązywania zarówno równania Schrödingera jak i klasycznego

równa-opis metod różnicowych stosowanych w dziedzinie półprzewodników można znaleźć w monografii Selberherr [88].

nia dyfuzji. W przypadku jednowymiarowym oba te równania mogą być zapisane następu-jąco: ∂u

∂t = F (x, t, u,∂u ∂x,∂2u

∂x2). Jeżeli oznaczymy uk

i ma przybliżać wartość u(xi, tk), to schemat Cranka-Nicolsona ma postać

uk+1 i = uk i +¹ 2^h ¡ Fk i + Fk+1 i ¢ , gdzie Fk i = F ³ xi, tk, uk i,∂u ∂x(xi, tk),∂2u ∂x2(xi, tk) ´ .

Elektrodyfuzja w stanie stacjonarnym

3.1 Wprowadzenie

Ogólny problem NPP jest opisany ewolucyjnym układem równań różniczko-wych cząstkoróżniczko-wych (3.1)    ∂ci ∂t = −divJ_i (i = 1, . . . , r), div(εE) = ρ, Ji = −Di∇ci+ e τziDiciE.

Rozwiązania takiego układu – uzupełnionego odpowiednimi warunkami po-czątkowymi i brzegowymi są funkcjami położenia i czasu. Jednakże w pewnych sytuacjach interesują nas przypadki szczególne, np. gdy proces jest w stanie, który nie zmienia się w czasie.

Z punktu widzenia termodynamiki należy wyróżnić i precyzyjnie odróż-niać dwie klasy takich stanów (Callen [12]): stany równowagi i stany

stacjo-narne.

• Stan równowagi – proste układy mogą znajdować się w szczególnych

stanach scharakteryzowanych przez kilka parametrów makroskopowych. W reprezentacji entropijnej oznacza to, że dana jest entropia jako funkcja pozostałych parametrów ekstensywnych, U, V, N1, . . ., Nr, gdzie

U – energia wewnętrzna, V – objętość, N_i – ilość substancji i. Oznacza to, że dana jest ekstensywna funkcja

(3.2) S = S(U, V, N₁, . . . , N_r).

Ekstensywność entropii S może być matematycznie wyrażona przez własność jednorodności stopnia 1, która oznacza, że

(3.3) S(λU, λV, λN₁, . . . , λN_r) = λS(U, V, N₁, . . . , N_r), 34

dla dowolnego λ > 0.

Istnienie takich szczególnych stanów to podstawowy postulat klasycznej termodynamiki fenomenologicznej (termostatyki). W terminologii Gibbsa [12], związek (3.2) nazywamy relacją fundamentalną. Zawarta jest w nim pełna informacja termodynamiczna o układzie. W stanie równowagi wszys-tkie parametry nie zależą od czasu i są stałe w całej przestrzeni zajmowanej przez dane fazy. Od strony czysto formalnej termodynamika klasyczna jest teorią badającą związki pomiędzy takim stanami równowagi materii.

• Stan stacjonarny – to taki stan, w którym funkcje opisujące układ

nie zależą od czasu. W takich stanach czas nie występuje jako zmienna fizyczna, ale wielkości fizyczne mogą zależeć od położenia. Podsta-wowa różnica pomiędzy stanami równowagi a stanami stacjonarnymi jest taka, że w tych drugich może występować zależność od zmiennej przestrzennej. Jest to na ogół związane z występowaniem niezerowych strumieni (ale niezależnych od czasu).

Tak więc, mówiąc niezbyt precyzyjnie można stany równowagi scharak-teryzować jako takie stany stacjonarne, w których nie występują strumienie makroskopowe.