Difusão ambipolar em plasmas

De Física Computacional
Ir para navegação Ir para pesquisar

Equação da difusão ambipolar

Esse trabalho tem como objetivo demonstrar uma resolução numérica para o caso unidimensional da difusão ambipolar de um plasma(gás formado de elétrons e íons). 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.


IMAGEM (1)

Quando a partição é removida ha uma pertubação instantanea na condição de equilibrio. Os elétrons por ser mais leve difundem-se mais rápido que os íons devido a sua menor massa.

Essa separação de cargas gera um campo elétrico de polarização que aumenta a taxa de difusão dos íons e diminui a dos elétrons até que ocorra um equilíbrio.

Esses movimentos coletivos são caracterizados por uma frequência natural de oscilação(Oscilações de Langmuir).


IMAGEM (2):


Como mostrado por Shimony e Cahn[1], esse problema é descrito por uma equação de onda amortecida para a função de densidade Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n(\vec r, t)} :

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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 Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u^2 = \nu_a D_a} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \alpha = 1/D_a} , sendo Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \nu_a} a frequência de colisão ambipolar e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle D_a} o coeficiente de difusão ambipolar.


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

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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 [2]. Começamos com a forma mais usual de escrever a equação da onda amortecida unidimencional

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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 Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 2h = \alpha u^2 = \nu_a} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle c^2 = u^2} .

Discretizando as variáveis do problema, temos que

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle x_i = i\Delta x \qquad i = 0,1,2,...,I} ,

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle t_k = k\Delta t \qquad k = 0,1,2,...,K} .

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

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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) } ,

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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) } ,

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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 Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle O\{\Delta t^2,\Delta x^2\}} e isolando Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u_i^{k_1}} , obtemos

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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 Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle s = (c^2\Delta t^2/\Delta x^2) } .

Essa é a equação para resolver o problema para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle k \geq 1} , porém o problema envolve uma derivada de segunda ordem no tempo, o que faz com que precisemos saber os dois passos anteriores para calcular o próximo. Então necessitamos ainda de uma maneira de determinar Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u_i^1} a paritr de Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u_i^0} , para então calcular os demais passos. Para isso assumimos que a função é inicialmente estacionária e fazemos

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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 Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle k = 0} obtemos

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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 (no caso desse trabalho usamos bordas fixas em 0 e também condições periódicas de contorno), podemos calcular a evolução temporal da função de densidade. Esse método é estável para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 0 \leq s \leq 1} e seu erro é

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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\} } .

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 Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle L=10} . Fizemos o plasma inicialmente concentrado na região central do tubo: Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n(x,0)=1} para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle L/4 \leq x \leq 3L/4} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n(x,0)=0} para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle x} fora dessa região. Usamos Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta x = 0.1} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta t = 0.01} , que resulta em Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle s = 0.01\nu_aD_a} , e criamos gifs mostrando a evlução temporal da função de densidade para diferentes valores de Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \nu_i} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle D_a} . Para garantir a estabilidade do método, essas constantes devem ser tais que Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 0 \leq \nu_iD_a \leq (\Delta x / \Delta t)^2} .

Quanto à solução nas bordas, fizemos de duas manerias: A primeira foi com condições de contorno fixas em 0 (Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n(0,t) = n(L,t) = 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 (Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n_{I+1}^k = n_0^k} ) e representa o caso de um tubo fechado.

Evolução temporal da densidade do plasma para diferentes valores de Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \nu_a} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle D_a}
Alt text
Alt text
Alt text
Alt text
Evolução temporal da densidade do plasma para diferentes valores de Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \nu_a} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle D_a}
Alt text
Alt text
Alt text
Alt text


Podemos observar que Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle D_a} é 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á Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \nu_a} 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 Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n_i^1} 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. ↑ Z. Shimony and J. H. Cahn, "Time-dependent ambipolar diffusion waves", The Physics of Fluids 8, 1704 (1965)
  2. ↑ H. Najafi and F. Izadi, "Comparison of two finite-difference methods for solving the damped wave equation", viXra, 2016
Disponível em “http://fiscomp.if.ufrgs.br/index.php?title=Difusão_ambipolar_em_plasmas&oldid=4147”