dr inż. Jerzy J. Jasielec

Model Nernsta-Plancka-Poissona

PMiKNoM seminaria - zajęcia 11 i 12

Wstęp

Elektro-dyfuzja jest to ruch jonów, wynikający z dwóch procesów:

dyfuzji - procesu, którego siłą napędową jest gradient stężeń
migracji - procesu, którego siłą napędową jest gradient potencjału

Procesy elektro-dyfuzji omawiane są w wielu dziedzinach nauki, wliczając naukę o półprzewodnikach, chemię analityczną czy neuro-fiziologię. Narzędzie do opisu elektro-dyfuzji znane jest jako model Nernsta-Plancka-Poissona i może być zastosowane do opisu procesów zachodzących w skali nano, jak i w skali makro. Więcej szczegółów w:

K. Szyszkiewicz, J.J. Jasielec, M. Danielewski and R. Filipek J. Electrochem. Soc. 2017, 164, E3559-E3568.

Model NPP - podstawowe wzory

Model NPP jest problemem początkowo-brzegowym. Poniżej pokazane jest rozwiązanie dla warunków brzegowych Dirichleta zarówno dla stężeń, jak i dla potencjału.

Prawo zachowania masy

Reakcje nie są uwzględnione, a zatem prawo zachowania masy przyjmuje postać:

$$ \frac{\partial c_i}{\partial t} = - \frac{\partial J_i}{\partial x} $$

gdzie `c_i` oznacza stężenie `i`-tego jonu, a `J_i` jego strumień; `x` oznacza przestrzeń, a `t` reprezentuje czas.

Równanie Nernsta-Plancka

Równanie Nernsta-Plancka opisuje strumień jonu, uwzględniając wpływ procesów dyfuzji oraz migracji (jako suma dwóch składników):

$$ J_i= - D_i \frac{\partial c_i}{\partial x} - \frac{F}{RT} z_i D_i c_i \frac{\partial \varphi}{\partial x} $$

gdzie `D_i` oznacza współczynnik dyfuzji i-tego jonu, a `z_i` jego ładunek; `\varphi` representuje potencjał elektryczny; `F` to stała Faradaya, `R` to stała gazowa, a `T` oznacza temperaturę.

Równanie Poissona

Równanie Poissona łączy potencjał z powstającym zaburzeniem elektro-obojętności:

$$ \frac{\partial^2 \varphi}{\partial x^2} = - \frac{F}{\varepsilon} \sum_i z_i c_i $$

gdzie `\varepsilon` oznacza przenikalność elektryczną ośrodka.

Warunki początkowe i brzegowe

Zakładamy stałe początkowe stężenia wszyskich jonów:

`c_i (x,0) = c_{i,ST} = const`

Stosujemy warunki brzegowe Dirichleta dla potencjału, t.j. zakładamy stałe wartości potencjału na brzegach:

$$ \left\{ \begin{array}{l} \varphi(0,t) = \varphi_L = const \\ \varphi(l,t) = \varphi_R = const \\ \end{array} \right. $$

Warunki brzegowe Dirichleta stosujemy również dla stężeń:

$$ \left\{ \begin{array}{l} c_i(0,t) = c_{i,L} \cdot k_{in} / k_{out} = const \\ c_i(l,t) = c_{i,R} \cdot k_{in} / k_{out} = const \end{array} \right. $$

Zamiennie można zastosować warunki brzegowe Changa-Jaffego:

$$ \left\{ \begin{array}{l} J_i(0,t) = k_{in} \cdot c_{i,L} - k_{out} \cdot c_i(0,t)\\ J_i(l,t) = k_{out} \cdot c_i(l,t) - k_{in} \cdot c_{i,R} \end{array} \right. $$

gdzie `k_{i n}` oraz `k_{out}` są heterogenicznymi stałymi przejścia `i`-tego wchodzącego do membrany i wychodzącego z membrany.

Model NPP - rozwiązanie numeryczne

Rozwiązanie numeryczne powyższych równań opiera się na Metodzie Linii oraz metodzie jawnej Eulera.

Metoda Linii dla stężeń

Dzielimy rozpatrywany odcinek, na `n + 1` przedziałów, otrzymując punkty `x_0` do `x_{n+1}`. Odległość między kolejnymi punktami wynosi `h`. Definiujemy `c_{i,k} \equiv c_i(x_k)`, oraz definiujemy punkty pomocnicze ` x_{k-1/2} ` oraz ` x_{k+1/2} ` w odległości `h/2` od punktu `x_k`. W tych punktach wyznaczać będziemy strumienie (`J_{i,k-1/2}` oraz `J_{i,k+1/2}`).

Prawo zachowania masy dla `i`-tego składnika w `k`-tym punkcie jest wyrażone jako:

$$ \frac{dc_{i,k}}{dt} = - \frac{J_{i,k+1/2}-J_{i,k-1/2}}{h} $$

Równanie Nernsta-Plancka dla strumieni `J_{i,k-1/2}` oraz `J_{i,k+1/2}`, zapisujemy następująco:

$$ \left\{ \begin{array}{l} J_{i,k-1/2}= - D_i \frac{d c_{i,k-1/2}}{dx} - \frac{F}{RT} z_i D_i c_{i,k-1/2} \frac{d \varphi_{k-1/2}}{dx} \\ J_{i,k+1/2}= - D_i \frac{d c_{i,k+1/2}}{dx} - \frac{F}{RT} z_i D_i c_{i,k+1/2} \frac{d \varphi_{k+1/2}}{dx} \end{array} \right. $$

Następnie pochodne, występujące w równaniu Nernsta-Plancka, możemy zastąpić ilorazami różnicowymi otrzymując:

$$ \left\{ \begin{array}{l} J_{i,k-1/2}= - D_i \frac{c_{i,k}-c_{i,k-1}}{h} - \frac{F}{RT} z_i D_i \frac{c_{i,k}+c_{i,k-1}}{2} \frac{\varphi_k -\varphi_{k-1}}{h} \\ J_{i,k+1/2}= - D_i \frac{c_{i,k+1}-c_{i,k}}{h} - \frac{F}{RT} z_i D_i \frac{c_{i,k+1}+c_{i,k}}{2} \frac{\varphi_{k+1} -\varphi_k}{h} \end{array} \right. $$

Po podstawieniu `J_{i,k-1/2}` oraz `J_{i,k+1/2}` do prawa zachowania masy otrzymujemy:

$$ \frac{dc_{i,k}}{dt} = \frac{D_i}{h^2} \left( c_{i,k+1}-2c_{i,k}+c_{i,k-1} \right) + \frac{F z_i D_i}{2RT h^2} \left( (c_{i,k+1}+c_{i,k})(\varphi_{k+1} -\varphi_k) - (c_{i,k}+c_{i,k-1})(\varphi_k -\varphi_{k-1}) \right). $$

Metoda Eulera

Jeżeli `c_{i,k}` oznacza stężenie w czasie `t_j`, a `c\text{_}n e w_{i,k}` oznacza stężenie w czasie `t_{j+1}`, to wtedy:

$$ c\text{_}new_{i,k} = c_{i,k} + dt \cdot \frac{dc_{i,k}}{dt} $$

Po wstawieniu wzoru na pochodną po czasie otrzymujemy:

$$ c\text{_}new_{i,k} = c_{i,k} + \alpha_i \left( c_{i,k+1}-2c_{i,k}+c_{i,k-1} \right) + \beta_i \bigl( (c_{i,k+1}+c_{i,k})(\varphi_{k+1} -\varphi_k) - (c_{i,k}+c_{i,k-1})(\varphi_k -\varphi_{k-1}) \bigr) $$

gdzie `\alpha_i = \frac{dt \cdot D_i}{h^2}` oraz `\beta_i = \frac{dt \cdot F z_i D_i}{2RT h^2}`

Metoda Linii dla równania Poissona

Pochodną po położeniu zastępujemy ilorazem różnicowym otrzymując:

$$ \frac{\varphi_{k+1}-2\varphi_k+\varphi_{k-1}}{h^2} = - \frac{F}{\varepsilon} \sum_i z_i c_{i,k} $$

Po pomnożeniu obustronnie przez `h^2` otrzymujemy równanie `i`-tego elementu:

$$ \varphi_{k-1} - 2\varphi_k + \varphi_{k+1} = - \gamma \sum_i z_i c_{i,k} $$

gdzie `\gamma = \frac{h^2 F}{\varepsilon}`

Wiedząc, że ` \varphi_0 = \varphi_L = const` oraz `\varphi_{n+1} =\varphi_R =const` do powyższego równania podstawiamy `k=1` oraz `k=n`, otrzymując:

$$ -2\varphi_1 + \varphi_2 = - \varphi_L - \gamma \sum_i z_i c_{i,1} $$ $$ \varphi_{n-1} - 2\varphi_{n} = - \varphi_R - \gamma \sum_i z_i c_{i,n} $$

Powyższe trzy równania stanowią układ równań macierzy trójdiagonalnej:

$$ \left[ \begin{array}{ccccccccc|c} -2 & 1 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & - \gamma \sum z_i c_{i,1} - \varphi_L \\ 1 & -2 & 1 & 0 & \cdots & 0 & 0 & 0 & 0 & - \gamma \sum z_i c_{i,2} \\ 0 & 1 & -2 & 1 & \cdots & 0 & 0 & 0 & 0 & - \gamma \sum z_i c_{i,3} \\ 0 & 0 & 1 & -2 & \cdots & 0 & 0 & 0 & 0 & - \gamma \sum z_i c_{i,4} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots\\ 0 & 0 & 0 & 0 & \cdots & -2 & 1 & 0 & 0 & - \gamma \sum z_i c_{i,n-3} \\ 0 & 0 & 0 & 0 & \cdots & 1 & -2 & 1 & 0 & - \gamma \sum z_i c_{i,n-2} \\ 0 & 0 & 0 & 0 & \cdots & 0 & 1 & -2 & 1 & - \gamma \sum z_i c_{i,n-1} \\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 1 & -2 & - \gamma \sum z_i c_{i,n} - \varphi_R \\ \end{array} \right] $$

Model NPP - algorytm programu

Deklaracje - globalne stałe i zmienne:

Zadeklaruj stałe globalne `l`, `t_{end}`, `p`, `n`, `m`, `F`, `R`, `T` oraz `\varepsilon`.
Zadeklaruj tablice dwuwymiarowe `c(1` to `p`, `0` to `n+1)` oraz `c__\text{new}(1` to `p`, `1` to `n)`, a także jednowymiarową tablicę `varphi(0` to `n+1)`.
Zadeklaruj tablice `D(1` to `p)`, `z(1` to `p)`, `c_L(1` to `p)`, `c_{ST}(1` to `p)`, `c_R(1` to `p)`, oraz zmienne `varphi_L` i `varphi_R`.
Zadeklaruj zmienne `h`, `dt` oraz `gamma`; tablice `alpha(1` to `p)` oraz `beta(1` to `p)`; a także zmienne `i`, `j`, `k`.

Procedura obliczająca potencjał elektryczny - Sub CalculatePotential():

Zadeklaruj tablice `a(1` to `n-1)`, `b(1` to `n)`, `c c(1` to `n-1)`, `dd(1` to `n)` oraz zmienną `s`.
Przypisz wartości diagonali, t.j. dla `k=1` do `n`: `dd(k)=-2`.
Przypisz wartości naddiagonalnej i poddiagonalnej, t.j. dla `k=1` do `n-1`: `c c(k)=1` oraz `a(k)=1`.
Przypisz wartości wyrazu wolnego używając pętli (dla `k=1` do `n`):
- wyzeruj zmienną `s`, t.j. `s = 0`,
- oblicz `\sum_i z_i c_{i,k}`, t.j. dla `i=1` do `p`: `s = s + z(i) \cdot c(i,k)`,
- przypisz `b(k) = -\gamma \cdot s`.
Zmodyfikuj wartość `b` dla pierwszego i ostatniego wiersza macierzy, t.j. `b(1)=b(1) - \varphi_L` oraz `b(n)=b(n) - \varphi_R`.
Rozwiąż macierz trójdiagonalną (tak jak to było robione na zajęciach 4).
Przepisz wyniki do tablicy potencjału, t.j. dla `k=1` do `n`: `\varphi (k)=b(k)`.

Procedura główna - Sub NPP():

Przypisz wartości do tablic `D`, `z`, `c_L`, `c_{ST}`, `c_R` oraz zmiennych `varphi_L` i `varphi_R`.
Wylicz wartość kroku `h=l/{n+1}` i kroku `dt=t_{end}/m` oraz wartości `alpha(i)`, `beta(i)` oraz `gamma`.
Przypisz wartości początkowe stężeń `c(i,k)=c_{ST}(i)`. (pętla zagnieżdżona dla `k=1` do `n` oraz dla `i=1` do `p`).
Przypisz wartości brzegowe stężeń, czyli dla `i=1` do `p`: `c(i,0)=c_L(i)` oraz `c(i,n+1)=c_R(i)`.
Przypisz wartości brzegowe potencjału, t.j. `varphi(0)=varphi_L` oraz `varphi(n+1)=varphi_R`.
Wywołaj procedurę liczącą potencjał, t.j. `Call` `Calcu latePotential`.
Wykonaj pętlę zagnieżdżoną (dla `j=1` do `m`), w której dla każdego kroku czasowego:
- program wyliczy nowe wartości stężeń wszystkich jonów (dla `k=1` do `n` oraz dla `i=1` do `p`), zapisując je w tablicy `c_\text{new}` (patrz wzór 3.2.);
- program podmieni tablicę stężeń, tj. (dla `k=1` do `n` oraz dla `i=1` do `p`) `c(i,k)=c_\text{new}(i,k)`;
- program wywoła funkcję wyliczającą potencjał, t.j. `Call` `Calcu latePotential`.
Zapisz wyniki do arkusza (w kolejnych kolumnach: położenie, stężenia kolejnych jonów i potencjał).

Model NPP - program

Model Nernsta-Plancka-Poissona - implementacja w VBA.
Porównanie wyników z aplikacją NPP napisaną w JavaScript.