Difusão ambipolar em plasmas

Equação da difusão ambipolar

A difusão é o modo como um fluido de dilui em um meio. Estudar as equações que governam esse fenômeno e as formas de resolvê-las é de extremo interesse para a física de fluidos e de plasmas, entre outras áreas. Aqui mostramos uma resolução numérica para o caso unidimensional da difusão ambipolar de um plasma (gás formado de elétrons e íons) envolto em um gás neutro, ou seja, o caso de um plasma se espalhando por um tubo.

Diferentemente de um gás de átomos/moléculas neutros(as), os plasmas são menos livres ao se moverem por causa das interações eletromagnéticas envolvidas no movimento das cargas, como a força de Coulomb e a força magnética. Na difusão de plasmas em um gás neutro, os coeficientes de difusão dos elétrons e dos íons são tipicamente dados por

$D_e=\frac{k_b{T}_e}{m_e{\nu }_e}$ e $D_i=\frac{h_b{T}_i}{m_i{\nu }_i}$

onde $T_e$, $T_i$, $m_e$, $m_i$, ${\nu }_e$ e ${\nu }_i$, são as temperaturas, massas e frequências de colisão dos elétrons e íons com os átomos neutros. Devido à massa do elétron ser muito menor que a massa de um íon, $D_e$ é maior que $D_i$, então quando um plasma começa a se difundir, incialmente os elétrons se espalham mais rapidamente que os íons, o que gera um campo elétrico que freia os elétron e acelera os íons. Chamamos esse processo de difusão ambipolar.

Como mostrado por Shimony e Cahn, esse problema é descrito por uma equação de onda amortecida

${\mathrm{\nabla }}^{2}n(\overrightarrow{r},t)=\frac{1}{u^2}\frac{{\partial }^{2}n(\overrightarrow{r},t)}{\partial t^2}+\alpha \frac{\partial n(\overrightarrow{r},t)}{\partial t}(1)$

onde $u^2={\nu }_a{D}_a$ e $\alpha =1/D_a$, sendo ${\nu }_a$ a frequência de colisão ambipolar e $D_a$ o coeficiente de difusão ambipolar, que pode ser escrito como $D_a=D_i(1+T_e/T_i)$ .

Como tratamos do caso unidimensional, a equação 1 torna-se

$\frac{{\partial }^{2}n(x,t)}{\partial x^2}=\frac{1}{u^2}\frac{{\partial }^{2}n(x,t)}{\partial t^2}+\alpha \frac{\partial n(x,t)}{\partial t}(2)$

O Método

A resolução numérica do problema foi baseada no artigo de Najafi e Izadi . Começamos com a forma mais usual de escrever a equação da onda amortecida unidimencional

$\frac{{\partial }^{2}n(x,t)}{\partial t^2}+2h\frac{\partial n(x,t)}{\partial t}=c^2\frac{{\partial }^{2}n(x,t)}{\partial x^2}(3)$

No nosso caso $2h=\alpha u^2={\nu }_a$ e $c^2=u^2$.

Discretizando as variáveis do problema, temos que

$x_i=i\mathrm{\Delta }xi=0,1,2,...,I$

$t_k=k\mathrm{\Delta }tk=0,1,2,...,K$

Substituindo as derivadas que aparecem na equação por diferenças finitas, obtemos

$\frac{{\partial }^{2}n}{\partial t^2}{|}_i^k=\frac{n_i^{k+1}-2n_i^k+n_i^{k-1}}{\mathrm{\Delta }{t}^2}-\frac{\mathrm{\Delta }{t}^2}{12}\frac{{\partial }^{2}n}{\partial t^4}{|}_i^k+O(\mathrm{\Delta }{t}^4)$

$\frac{\partial n}{\partial t}{|}_i^k=\frac{n_i^{k+1}-n_i^{k-1}}{2\mathrm{\Delta }t}-\frac{\mathrm{\Delta }{t}^2}{6}\frac{{\partial }^{3}n}{\partial t^3}{|}_i^k+O(\mathrm{\Delta }{t}^4)$

$\frac{{\partial }^{2}n}{\partial x^2}{|}_i^k=\frac{n_{i+1}^k-2n_i^k+n_{i-1}^k}{\mathrm{\Delta }{x}^2}-\frac{\mathrm{\Delta }{x}^2}{12}\frac{{\partial }^{2}n}{\partial x^4}{|}_i^k+O(\mathrm{\Delta }{x}^4)$

Substituindo essas relações na equação 3, obtemos

$[\frac{n_i^{k+1}-2n_i^k+n_i^{k-1}}{\mathrm{\Delta }{t}^2}-\frac{\mathrm{\Delta }{t}^2}{12}\frac{{\partial }^{4}n}{\partial t^4}{|}_i^k+O(\mathrm{\Delta }{t}^4)]+2h[\frac{n_i^{k+1}-n_i^{k-1}}{2\mathrm{\Delta }t}-\frac{\mathrm{\Delta }{t}^2}{6}\frac{{\partial }^{3}n}{\partial t^3}{|}_i^k+O(\mathrm{\Delta }{t}^4)]=c^2[\frac{n_{i+1}^k-2n_i^k+n_{i-1}^k}{\mathrm{\Delta }{x}^2}-\frac{\mathrm{\Delta }{x}^2}{12}\frac{{\partial }^{4}n}{\partial x^4}{|}_i^k+O(\mathrm{\Delta }{x}^4)]$

Omitindo todos os temos de ordem $O\mathrm{\Delta }{t}^2,\mathrm{\Delta }{x}^2$ e isolando $u_i^{k_1}$, obtemos

$u_i^{k+1}=\frac{1}{1+h\mathrm{\Delta }t}[2(1-s)n_i^k-(1-h\mathrm{\Delta }t)n_i^{k-1}+s(n_{i+1}^k+n_{i-1}^k)](4)$

sendo $s=(c^2\mathrm{\Delta }{t}^2/\mathrm{\Delta }{x}^2)$.

Essa é a equação para resolver o problema para $k\ge 1$, mas nessecitamos ainda de uma maneira de determinar $u_i^1$ a paritr de $u_i^0$. Para isso assumimos que a função é inicialmente estacionária e fazemos

$n_t(x,0)=0\Rightarrow \frac{n_i^1-n_i^{-1}}{\mathrm{\Delta }t}=0\Rightarrow n_i^1=n_i^{-1}$

Substituindo na equação 4 para $k=0$ obtemos

$n_i^1=\frac{1}{2}[2(1-s)n_i^0+s(n_{i+1}^0+n_{i-1}^0)](5)$

Com as equações 4 e 5, e tomando as devidas condições de contorno nas bordas, podemos calcular a evolução temporal da função de densidade. Esse método é estável para $0\le s\le 1$ e seu erro é

$E\{n_i^k\}=-\frac{\mathrm{\Delta }{t}^2}{12}\frac{{\partial }^{4}n}{\partial t^4}{|}_i^k-\frac{h\mathrm{\Delta }{t}^2}{3}\frac{{\partial }^{3}n}{\partial t^3}{|}_i^k+c^2\frac{\mathrm{\Delta }{x}^2}{12}\frac{{\partial }^{4}n}{\partial x^4}{|}_i^k+O\{\mathrm{\Delta }{t}^4,\mathrm{\Delta }{x}^4\}$

Resultados e Discussão

Evolução temporal da densidade do plasma para diferentes valores de ${\nu }_a$ e $D_a$