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

$\nabla ^{2}n({\vec {r}},t)={\frac {1}{u^{2}}}{\frac {\partial ^{2}n({\vec {r}},t)}{\partial t^{2}}}+\alpha {\frac {\partial n({\vec {r}},t)}{\partial t}}\qquad (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}}\qquad (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}}}\qquad (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\Delta x\qquad i=0,1,2,...,I$

$t_{k}=k\Delta t\qquad k=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}}{\Delta t^{2}}}-{\frac {\Delta t^{2}}{12}}{\frac {\partial ^{2}n}{\partial t^{4}}}|_{i}^{k}+O(\Delta t^{4})$

${\frac {\partial n}{\partial t}}|_{i}^{k}={\frac {n_{i}^{k+1}-n_{i}^{k-1}}{2\Delta t}}-{\frac {\Delta t^{2}}{6}}{\frac {\partial ^{3}n}{\partial t^{3}}}|_{i}^{k}+O(\Delta t^{4})$

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

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

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

Omitindo todos os temos de ordem $O{\Delta t^{2},\Delta x^{2}}$ e isolando $u_{i}^{k_{1}}$, obtemos

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

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

Essa é a equação para resolver o problema para $k\geq 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}}{\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})]\qquad (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\leq s\leq 1$ e seu erro é

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