Belousov-Zhabotinsky

Belousov-Zhabotinsky Reaction

A reação de Belousov-Zhabotinsky (BZ) consiste em uma família de reações químicas oscilatórias descobertas inicialmente por Belousov, e posteriormente analisadas por Zhabotinsky. A reação consiste em 3BrO₃⁻ + 5CH₂(CO₂H)₂ + 3H⁺ → 3BrCH(CO₂H)₂ + 4CO₂ + 5H₂O + 2CH₂O₂, e demonstra um comportamento oscilatório não linear até atingir o equilíbrio químico. A interação entre a reação e a difusão dos produtos químicos no espaço resultará na auto-organização de ondas viajantes dinâmicas. Seu mecanismo original, foi descrito através de 27 espécies químicas e um total de 80 reações.

Reação de Belousov-Zhabotinsky em uma placa de Petri.

Oregonator

Oregonator é um modelo matemático utilizado para descrever de forma mais simples a dinâmica da reação BZ, desenvolvido por Field e Noyes (1974). Foi um modelo não espacial originalmente composto por três variáveis de estado, onde posteriormente, vemos que tornam-se apenas duas. O mecanismo é, inicialmente composto por cinco etapas irreversíveis, onde, A = 3BrO₃ ^-, B = 5CH₂(COOH)₂; 2HCOOH, 3BrCH(COOH)₂ (no geral, estas e demais espécies orgânicas); P = HOBr; X = HBrO₂; Y = Br^-; Z = forma oxidada do catalisador e f é o coeficiente estequiométrico. Observamos também. suas respectivas equações de taxa, onde v é a taxa da reação e k_i corresponde às constantes de taxa de reação:

A + Y

$\longrightarrow$

X + P

v1 = k1 [A][Y]

X + Y

$\longrightarrow$

2 P

v2 = k2 [X][Y]

A + X

$\longrightarrow$

2 X + 2 Z

v3 = k3 [A][X]

2 X

$\longrightarrow$

A + P

v4 = k4 [X]2

B + Z

$\longrightarrow$

${\frac {1}{2}}f$ Y

v5 = k5 [B][Z]

Para construir o modelo Oregonator, é necessário supor que as concentrações de A e B permaneçam constantes (estão associadas às concentrações iniciais dos precursores). Posteriormente, deve-se aplicar as técnicas padrão de cinética química para obter o modelo dinâmico considerando X, Y e Z como variáveis dinâmicas, assumindo que as reações químicas são elementares, ou seja, os coeficientes estequiométricos coincidem com a potência das variáveis dinâmicas. Considerando $\tau$ como o tempo, vemos as seguintes equações de velocidade:

${\frac {d[X]}{d\tau }}=k_{1}[A][Y]-k_{2}[X][Y]+k_{3}[A][X]-2k_{4}[X]^{2}$

${\frac {d[Y]}{d\tau }}=-k_{1}[A][Y]-k_{2}[X][Y]+{\frac {1}{2}}fk_{5}[B][Z]$

${\frac {d[Z]}{d\tau }}=2k_{3}[A][X]-k_{5}[B][Z]$

A análise é simplificada convertendo essas equações em uma forma adimensional:

$x\equiv {\tfrac {2k_{4}[X]}{k_{3}[Y]}}$

$y\equiv {\tfrac {k_{2}[X]}{k_{3}[A]}}$

$z\equiv {\tfrac {k_{5}k_{4}[B][Z]}{(k_{3}[A])^{2}}}$

$t\equiv k_{5}[B]\tau$

A partir de operações algébricas com as equações acima, obtemos para x, y e z, o seguinte sistema de equações diferenciais não lineares:

${\begin{cases}{\frac {dx}{dt}}={\frac {qy-xy+x(1-x)}{\epsilon }}\\{\frac {dy}{dt}}={\frac {-qy-xy+fz}{\epsilon '}}\\{\frac {dz}{dt}}=x-z\\\end{cases}}$

Onde $\epsilon \equiv {\tfrac {k_{5}[B]}{k_{3}[A]}}$, $\epsilon '\equiv {\tfrac {2k_{4}k_{5}[B]}{k_{2}k_{3}[A]}}$ e $q\equiv {\tfrac {2k_{4}k_{1}}{k_{2}k_{3}}}$, e seus respectivos valores típicos são $\epsilon \approx 10^{-2}$,$\epsilon '\approx 10^{-5}$ e $q\approx 10^{-4}$. Como o parâmetro $\epsilon '\approx 10^{-5}$ (que é obtido através de reações experimentais), é possível ver que y muda na escala do tempo de forma muito mais rápida que as demais variaveis, e então, assumimos que $y(h,t)$ é sempre determinado pelos valores instantâneos de $x$ e $z$, e assim, reescrevemos y como $y\equiv {\tfrac {fz}{q+x}}$:. Deste modo, as equações são reduzidas para :

$\epsilon {\frac {dx}{dt}}=x(1-x)+f{\frac {q-x}{q+x}}z$

${\frac {dz}{dt}}=x-z$

Dessa forma, o modelo de Oregonator mostra a forma típica de um sistema de feedback químico, ou seja, a variável x, que será reescrito como u, funciona como um ativador, enquanto a variável z, que será reescrita como v, tem o papel de inibidor. Se para as equações termos associados à difusão são adicionados, onde D_u e D_v são os coeficientes de difusão adimensionais, e $\nabla ^{2}$ é o operador Laplaciano, então, o sistema torna-se:

$\epsilon {\frac {du}{dt}}=u(1-u)+f{\frac {q-u}{q+u}}v+D_{u}\nabla ^{2}u$

${\frac {dv}{dt}}=u-v+D_{v}\nabla ^{2}v$

Implementação

Antes de discretizarmos a equação para que assim possamos utiliza-la em um código, explicaremos brevemente métodos e fórmulas utilizados para isso.

Método FTCS (Forward Time Centered Space)

De modo a resolver numericamente as equações descritas acima, serão utilizado o método FTCS (Forward Time Centered Space), que consiste em um método para resolver equações parciais através da derivada parcial de primeira ordem no tempo por uma diferença finita e progressiva e a derivada parcial de segunda ordem no espaço por uma diferença centrada, como vemos logo abaixo:

${\frac {\partial f}{\partial t}}\approx {\frac {f(h,t+dt)-f(h,dt)}{\Delta t}}={\frac {f_{i}^{n+1}-f_{i}^{n}}{\Delta t}}$

${\frac {\partial ^{2}f}{\partial h^{2}}}\approx {\frac {f(h,t-dt)-2f(h,t)+f(h,t+1)}{\Delta h^{2}}}={\frac {f_{i}^{n-1}-2f_{i}^{n}+f_{i}^{n+1}}{\Delta h^{2}}}$

Laplaciano

O Laplaciano pode tanto ser representado por $\nabla ^{2}$ quanto por $\Delta$. O laplaciano é a soma de todas as derivadas parciais simples de segunda ordem:

Seja $u:\mathbb {R} ^{n}\to \mathbb {R}$

$\Delta u=\sum _{i=1}^{n}{\frac {\partial ^{2}u}{\partial x_{i}^{2}}}$

Como o código para Belousov-Zhabotinsky terá apenas 2 dimensões o $\Delta f(h,t)$ ficará da seguinte forma:

$\nabla ^{2}f(h,t)={\frac {\partial ^{2}f}{\partial x^{2}}}+{\frac {\partial ^{2}f}{\partial y^{2}}}$

$\nabla ^{2}f(h,t)={\frac {f(x+dh,y,t)+f(x-dh,y,t)+f(x,y+dh,t)+f(x,y-dh,t)-4f(x,y,t)}{\Delta h^{2}}}$

De modo a simplificar a expressão para analises posteriores, reescreveremos desta forma, considerando uma analise 2D:

$\nabla ^{2}f(h,t)={\frac {f_{R}+f_{L}+f_{U}+f_{D}-4f_{C}}{\Delta h^{2}}}$

Aplicação dos Métodos para a Reação de Belousov-Zhabotinsky

Considerando a equação $\epsilon {\frac {\partial u}{\partial t}}=u(1-u)+f{\frac {q-u}{q+u}}v+D_{u}\nabla ^{2}u$:

$\epsilon {\frac {u_{C}^{n+1}-u_{C}^{n}}{\Delta t}}=u_{C}(1-u_{C})+f{\frac {q-u_{C}}{q+u_{C}}}v_{C}+D_{u}\left({\frac {u_{R}+u_{L}+f_{U}+u_{D}-4u_{C}}{\Delta h^{2}}}\right)$

$u_{C}^{n+1}-u_{C}^{n}=\left[u_{C}(1-u_{C})+f{\frac {q-u_{C}}{q+u_{C}}}v_{C}+D_{u}\left({\frac {u_{R}+u_{L}+u_{U}+u_{D}-4u_{C}}{\Delta h^{2}}}\right)\right]{\frac {\Delta t}{\epsilon }}$

$u_{C}^{n+1}=u_{C}+\left[u_{C}(1-u_{C})+f{\frac {q-u_{C}}{q+u_{C}}}v_{C}+D_{u}\left({\frac {u_{R}+u_{L}+u_{U}+u_{D}-4u_{C}}{\Delta h^{2}}}\right)\right]{\frac {\Delta t}{\epsilon }}$

Considerando a equação ${\frac {\partial v}{\partial t}}=u-v+D_{v}\nabla ^{2}v$:

${\frac {v_{C}^{n+1}-v_{C}^{n}}{\Delta t}}=u_{C}-v_{C}+D_{v}\left({\frac {v_{R}+v_{L}+v_{U}+v_{D}-4v_{C}}{\Delta h^{2}}}\right)$

$v_{C}^{n+1}-v_{C}^{n}=\left[u_{C}-v_{C}+D_{v}\left({\frac {v_{R}+v_{L}+v_{U}+v_{D}-4v_{C}}{\Delta h^{2}}}\right)\right]\Delta t$

$v_{C}^{n+1}=v_{C}+\left[u_{C}-v_{C}+D_{v}\left({\frac {v_{R}+v_{L}+v_{U}+v_{D}-4v_{C}}{\Delta h^{2}}}\right)\right]\Delta t$

Resultados

Analise da concentração de u com $(e,q,f)=(0.2,10^{-3},1)$ $(D_{u},D_{v})=(10^{-5},10^{-5})$
t = 0.	t = 12000.	t = 25000.
t = 38000.	t = 50000.	t = 63000.
t = 75000.	t = 88000.	t = 100000.

BZ com $(e,q,f)=(0.2,10^{-3},1)$ $(D_{u},D_{v})=(10^{-5},10^{-5})$
BZ da concentração de $u$ até t = 20k.	BZ da concentração de $v$ até t = 20k.

Programas Utilizados

Simulaçao Belousov-Zhabotinsky