dr inż. Jerzy J. Jasielec

Model dyfuzji-reakcji

PMiKNoM seminaria - zajęcia 10

Wstęp

Dyfuzja chlorków w betonach a porowatość i wiązanie.

Model wykonany na zajęciach 8 zakładał że transport jonów chlorkowych odbywa się w jednorodnym ośrodku ciągłym. W rzeczywistości, płyta betonowa jest ośrodkiem porowatym, a dyfuzja jonów zachodzi głównie w cieczy porowej, w porach betonu. W celu uwzględnienia porowatości lewą stronę prawa Fick'a należy pomnożyć przez porowatość `\phi`:

$$ \phi \frac{\partial c}{\partial t} = D \frac{ \partial^2 c}{\partial x^2} $$

Dodatkowo zachodzi reakcja adsorpcji i część jonów zostaje związana w fazie stałej:

$$ \phi \frac{\partial c}{\partial t} = D \frac{ \partial^2 c}{\partial x^2} - r $$

Człon reakcyjny `r` oznacza jony, które zostaną zaadsorbowane w czasie `\partial t`, a zmiana stężenie chlorków związanych `{\partial C_b}/{\partial t}` opisywana jest równaniem:

$$ (1 -\phi) \rho_s \frac{\partial C_b}{\partial t} = r $$

W równaniu pojawia się gęstość fazy stałej `rho_s`, ponieważ `c` i `C_b` wyrażone są w innych jednostkach, t.j. [`kg`/`m^3`] i [`kg`/`kg`].

Adsorpcja (człon reakcyjny) może być opisana przy pomocy izoterm: liniowej, Langmuira lub Freunlicha. My użyjemy tej ostatniej. Szczegóły modelu w tym równania, rozwiązanie numeryczne oraz algorytm programu znajdują się poniżej.

Czasozależny model dyfuzji-reakcji - podstawowe wzory

Dyfuzja i wiązanie

Po uwzględnieniu informacji na temat porowatości betonu oraz wiązania chlorków w fazie stałej, prawo zachowania masy przybiera postać:

$$ \left\{ \begin{array}{l} \phi \frac{\partial c}{\partial t} = D \frac{ \partial^2 c}{\partial x^2} - r \\ (1 -\phi) \rho_s \frac{\partial C_b}{\partial t} = r \end{array} \right. $$

gdzie `c` oznacza stężenie chlorków wolnych w cieczy porowej, `C_b` oznacza stężenie chlorków związanych, `\phi` oznacza porowatość, `D` to współczynnik dyfuzji, `x` to położenie, a `t` oznacza czas.

Człon reakcyjny, `r` wyznaczany jest przy użyciu izotermy Freundlicha:

$$ r = k [c - (C_b \text{/} K_b)^{1\text{/}\eta}] $$

gdzie `k`, `K_b` oraz `\eta` są stałymi reakcji.

Gęstość fazy stałej w betonie, `\rho_s` wyliczana jest z sumy ważonej:

$$ \rho_s = \frac{\rho_c - \phi \cdot \rho_w}{1 - \phi} $$

gdzie `rho_c` oznacza gęstość całkowitą betonu, a `rho_w = 1000 {kg}/m^3` to gęstość wody.

Warunki początkowe i brzegowe

Zakładamy stałe stężenia na brzegach (warunki brzegowe Dirichleta):

$$ \left\{ \begin{array}{l} c(0,t) = c_L = const \\ c(l,t) = c_R = const \end{array} \right. $$

Na początku procesu, próbka nie zawiera chlorków:

$$ \left\{ \begin{array}{l} c(x,0) = 0 \\ C_b(x,0) = 0 \end{array} \right. $$

Całkowite stężenie chlorków

Całkowite stężenie chlorków (podana w `kg`/`kg` betonu) jest sumą chlorków wolnych i związanych:

$$ C_{t o t}(x,t) = C_{f r e e}(x,t) + C_{b o u n d}(x,t), $$

gdzie `C_{f r e e}(x,t) = \frac{\phi \cdot c(x,t)}{\rho_c}` oraz `C_{b o u n d}(x,t) = \frac{(1- \phi) \cdot \rho_s \cdot C_b(x,t)}{\rho_c}`.

Model dyfuzji-reakcji - rozwiązanie numeryczne

Rozwiązanie numeryczne powyższych równań bazuje na Metodzie Linii (MoL) oraz na metodzie jawnej Eulera.

Metoda Linii

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 \equiv c(x_i)` oraz `C_{b,i} \equiv C_b(x_i)`. Drugą pochodną stężenia po położeniu zastępujemy ilorazem różnicowym, otrzymując następujący układ równań:

$$ \left\{ \begin{array}{lll} \phi \frac{dc_i}{dt} = \frac{D}{h^2} (c_{i+1}-2c_i+c_{i-1}) - k [c_i - (C_{b,i} \text{/} K_b)^{1\text{/}\eta}] & & \text{ dla } i=1..n \\ (1 -\phi) \rho_s \frac{dC_{b,i}}{dt} = k [c_i - (C_{b,i} \text{/} K_b)^{1\text{/}\eta}] & & \text{ dla } i=0..n+1 \end{array} \right. $$

Metoda Eulera

Jeżeli `c_{i}` i `C_{b,i}` oznaczają stężenia w czasie `t_j`, a `\overline c_{i}` i `\overline C_{b,i}` oznaczają stężenia w czasie `t_{j+1}`, to wtedy:

$$ \left\{ \begin{array}{ll} \overline c_{i} = c_{i} + dt \cdot \frac{dc_{i}}{dt} \\ \overline C_{b,i} = C_{b,i} + dt \cdot \frac{dC_{b,i}}{dt} \end{array} \right. $$

Po wstawieniu wzoru na pochodną po czasie otrzymujemy:

$$ \left\{ \begin{array}{ll} \overline c_{i} = c_{i} + \alpha \left( c_{i+1}-2c_i+c_{i-1} \right) - \beta \bigl( c_i - (C_{b,i} \text{/} K_b)^{1\text{/}\eta} \bigr) \\ \overline C_{b,i} = C_{b,i} + \gamma \bigl( c_i - (C_{b,i} \text{/} K_b)^{1\text{/}\eta} \bigr) \end{array} \right. $$

gdzie `\alpha = \frac{dt \cdot D}{\phi \cdot h^2}`, `\beta = \frac{dt \cdot k}{\phi}` and `\gamma = \frac{dt \cdot k}{(1 -\phi) \rho_s}`.

Model dyfuzji-reakcji - algorytm programu

Zadeklaruj stałe globalne `c_L`, `c_R`, `l`, `D`, `t_{end}`, `n`, `m`, a także `row`, `roc`, `fi`, `k`, `Kb` oraz `eta`.
Procedurę główną rozpocznij od deklaracji zmiennych `h`, `dt`, `ros`, `\alpha`, `\beta`, `\gamma`, `i`, `j` oraz tablic `c(0 \text{ to } n+1)`, `c_\text{new}(1 \text{ to } n)` oraz `C_b(0 \text{ to } n+1)`.
Wylicz wartość kroku `h` i kroku `dt`, gęstość `ros` oraz zmienne `\alpha`, `\beta`, `\gamma`.
Zapisz warunek początkowy dla wszystkich punktów: `c(i)=0` i `C_b(i)=0` .
Zapisz warunek brzegowy: `c(0)=c_L` i `c(n+1)=c_R`.
Wykonaj pętlę zagnieżdżoną, w której dla każdego kroku czasowego (`j=1` do `m`) wykonane będą trzy pętle wewnętrzne:
- Program wyliczy nowe wartości stężeń chlorków wolnych, zapisując je w tablicy `c_\text{new}`, t.j.
  dla `i = 1` do `n` `c_{\text{new},i} = c_{i} + \alpha ( c_{i+1}-2c_i+c_{i-1} ) - \beta ( c_i - (C_{b,i} \text{/} K_b)^{1\text{/}\eta} )`.
- Program wyliczy nowe wartości stężeń chlorków związanych, zapisując w tablicy `C_b`, t.j.:
  dla `i = 0` do `n+1` `C_{b,i} = C_{b,i} + \gamma ( c_i - (C_{b,i} \text{/} K_b)^{1\text{/}\eta} )`.
- Program podmieni tablicę stężeń chlorków wolnych, tj.
  dla `i = 1` do `n` `c_i=c_{\text{new},i}`.
Wyświetl wyniki, w kolejnych kolumnach: `i \cdot h`, `C_{f r e e}(i)`, `C_{b o u n d}(i)`, `C_\{t o t}(i)`. (Patrz równanie 2.3.)

Model Dyfuzji-Reakcji - program

Model Dyfuzji-Reakcji - implementacja w VBA.
Porównanie wyników z aplikacją napisaną w JavaScript.