Difusão ambipolar em plasmas

De Física Computacional
Revisão de 12h05min de 5 de abril de 2021 por Leirbag (discussão | contribs)
Ir para navegação Ir para pesquisar

Equação da difusão ambipolar

A difusão é o modo como um fluido se 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

e ,

onde , , , , e , 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, é maior que , 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.

Difusao ambipolar.png [1]

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

,

onde e , sendo a frequência de colisão ambipolar e o coeficiente de difusão ambipolar, que pode ser escrito como [3].

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

.

O Método

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


No nosso caso e .

Discretizando as variáveis do problema, temos que

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

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

Omitindo todos os temos de ordem e isolando , obtemos

sendo .

Essa é a equação para resolver o problema para , mas nessecitamos ainda de uma maneira de determinar a paritr de . Para isso assumimos que a função é inicialmente estacionária e fazemos

Substituindo na equação 4 para obtemos

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 e seu erro é

Resultados e Discussão

Aplicamos o método descrito acima para simular a evolução da densidade de um plasma se difundindo em um tubo de largura . Fizemos o plasma inicialmente concentrado na região central do tubo: para e para fora dessa região. Usamos e , e criamos gifs mostrando a evlução temporal da função de densidade para diferentes valores de e . Para garantir a estabilidade do método, essas constantes devem ser tais que .

Quanto à solução nas bordas, fizemos de duas manerias: A primeira foi com condições de contorno fixas em 0 () e represenda o caso em que as bordas são um sumidouro, como se fosse um tubo aberto. A segunda foi usando condições de contorno periódicas () e representa o caso de um tubo fechado.

Evolução temporal da densidade do plasma para diferentes valores de e
Alt text
Alt text
Alt text
Alt text
Evolução temporal da densidade do plasma para diferentes valores de e
Alt text
Alt text
Alt text
Alt text


Podemos observar que é o parâmetro que domina a velocidade com que que a densidade decai, o que é esperado, uma vez que um coeficiente de difusão maior faz o plasma se difundir mais rápido. Já parece estar ligado à "suavidade" da distribuição, sendo que com frequências baixas começam a aparecer diversos picos de densidade.

Programas Utilizados

Para implementar o método computacionalmente e criar os gifs foram usados códigos em python. O código abaixo é a solução para as bordas fixas em 0:

import numpy as np
import matplotlib.pyplot as plt
import imageio

dt = 0.01
dx = 0.1
L = 10
T = 50

#constantes do plasma
nu_a = 0.1
Da = 0.5

#constantes para a eq da onda
c = np.sqrt(nu_a*Da)
h = nu_a/2
s = (c*dt/dx)**2

x = np.linspace(0, L, int(L/dx)) #array com as coordenadas espaciais
t = np.linspace(0, T, int(T/dt)) #array com as coordenadas temporais
n = np.zeros((len(t), len(x)))   #matriz com a densidade n(x,t)

#fazemos o plasma inicialmente concentrado em uma regiao
for i in range(int(len(x)/4),int(3*len(x)/4)):
    n[0,i] = 1
plt.plot(x,n[0]) #plota estado inicial da funcao
plt.title(r'$\nu_a=$'+str(nu_a)+'   $D_a=$'+str(Da))
plt.xlabel('L')
plt.ylabel('n(x)')
plt.ylim([0,1.1])
plt.xlim([0,L])
plt.text(8,0.9,'T = 0.0')
plt.savefig('n_0.png')
plt.clf()

#calculamos o próximo passo considerando dn/dt = 0 inicialmente    
for i in range(1,len(x)-1):
    n[1,i] = 0.5*(2*(1-s)*n[0,i] + s*(n[0,i+1] + n[0,i-1]))

#calculamos a posterior evolucao
for k in range(1, len(t)-1):
    for i in range(1,len(x)-1): #isso fixa os contornos em 0
        n[k+1,i] = (1/(1+h*dt))*(2*(1-s)*n[k,i] - (1-h*dt)*n[k-1,i] + s*(n[k,i+1] + n[k,i-1]))
    if k*dt - int(k*dt) == 0: #plota figuras para valores inteiros de t
       plt.plot(x,n[k])
       plt.title(r'$\nu_a=$'+str(nu_a)+'   $D_a=$'+str(Da))
       plt.xlabel('L')
       plt.ylabel('n(x)')
       plt.xlim([0,L])
       plt.ylim([0,1.1])
       plt.text(8,0.9,'T = '+str(k*dt))
       plt.savefig('n_'+str(int(k*dt))+'.png')
       plt.clf()
    
#criamos gifs usando os plots
images = []
for k in range(T):
    images.append(imageio.imread('n_'+str(k)+'.png'))
imageio.mimsave('difusao_ambipolar_'+str(nu_a)+'_'+str(Da)+'.gif', images, format='GIF', duration=1./10)

O código para as condições de contorno periódicas é levemente diferente, exigindo o acrécimo de duas duas linhas dentro do terceiro loop. Vale notar que aqui não nos preocupamos com o contorno no cáculo de porque a função nos contornos é 0 nos primeiros passos, independente de adotarmos ou não contornos periódicos.

import numpy as np
import matplotlib.pyplot as plt
import imageio

dt = 0.01
dx = 0.1
L = 10
T = 50

#constantes do plasma
nu_a = 0.1
Da = 0.5

#constantes para a eq da onda
c = np.sqrt(nu_a*Da)
h = nu_a/2
s = (c*dt/dx)**2

x = np.linspace(0, L, int(L/dx)) #array com as coordenadas espaciais
t = np.linspace(0, T, int(T/dt)) #array com as coordenadas temporais
n = np.zeros((len(t), len(x)))   #matriz com a densidade n(x,t)

#fazemos o plasma inicialmente concentrado em uma regiao
for i in range(int(len(x)/4),int(3*len(x)/4)):
    n[0,i] = 1
plt.plot(x,n[0]) #plota estado inicial da funcao
plt.title(r'$\nu_a=$'+str(nu_a)+'   $D_a=$'+str(Da))
plt.xlabel('L')
plt.ylabel('n(x)')
plt.ylim([0,1.1])
plt.xlim([0,L])
plt.text(8,0.9,'T = 0.0')
plt.savefig('n_0.png')
plt.clf()

#calculamos o próximo passo considerando dn/dt = 0 inicialmente    
for i in range(1,len(x)-1):
    n[1,i] = 0.5*(2*(1-s)*n[0,i] + s*(n[0,i+1] + n[0,i-1]))

#calculamos a posterior evolucao
for k in range(1, len(t)-1):
    for i in range(1,len(x)-1):
        n[k+1,i] = (1/(1+h*dt))*(2*(1-s)*n[k,i] - (1-h*dt)*n[k-1,i] + s*(n[k,i+1] + n[k,i-1]))
    #condicoes de contorno periodicas
    n[k+1,0] = (1/(1+h*dt))*(2*(1-s)*n[k,0] - (1-h*dt)*n[k-1,0] + s*(n[k,1] + n[k,-1]))
    n[k+1,-1] = (1/(1+h*dt))*(2*(1-s)*n[k,-1] - (1-h*dt)*n[k-1,-1] + s*(n[k,0] + n[k,-2]))
    if k*dt - int(k*dt) == 0: #plota figuras para valores inteiros de t
       plt.plot(x,n[k])
       plt.title(r'$\nu_a=$'+str(nu_a)+'   $D_a=$'+str(Da))
       plt.xlabel('L')
       plt.ylabel('n(x)')
       plt.xlim([0,L])
       plt.ylim([0,1.1])
       plt.text(8,0.9,'T = '+str(k*dt))
       plt.savefig('n_'+str(int(k*dt))+'.png')
       plt.clf()
    
#criamos gifs usando os plots
images = []
for k in range(T):
    images.append(imageio.imread('n_'+str(k)+'.png'))
imageio.mimsave('difusao_ambipolar_'+str(nu_a)+'_'+str(Da)+'PBC.gif', images, format='GIF', duration=1./10)


Referências

  1. http://www.enigmatic-consulting.com/semiconductor_processing/CVD_Fundamentals/plasmas/ambipolar_diffusion.html
  2. Z. Shimony and J. H. Cahn, "Time-dependent ambipolar diffusion waves", The Physics of Fluids 8, 1704 (1965)
  3. . http://uigelz.eecs.umich.edu/classes/pub/eecs517/handouts/derivation_ambipolar_diffusion_v02.pdf
  4. H. Najafi and F. Izadi, "Comparison of two finite-difference methods for solving the damped wave equation", viXra, 2016