Física de plasma espacial utilizando simulação computacional de partículas

Space plasma physics via computational particle simulation

ARTIGOS GERAIS

Física de plasma espacial utilizando simulação computacional de partículas

Space plasma physics via computational particle simulation

F.J.R. Simões Jr.I, ¹; E. Costa Jr.^II; M.V. Alves e F.R. Cardoso^II

^IDepartamento de Física, Universidade Federal de Pelotas, Campus UFPel Capão do Leão, Pelotas, RS, Brasil

^IILaboratório Associado de Plasma, INPE, São José dos Campos, SP, Brasil

RESUMO

Neste artigo convidamos os leitores a descobrir o mundo dos "experimentos computacionais" (ou simulações computacionais) para plasma espacial. O experimento computacional é semelhante ao experimento tradicional realizado em laboratório. Ele não é, e não será capaz de solucionar todos os problemas físicos que permanecem incompreendidos, mas fornece uma poderosa ferramenta para auxiliar os físicos (cientistas) a buscar soluções dos mesmos. Apresentamos neste trabalho algumas características das simulações computacionais e exemplos onde estas simulações podem ser aplicadas. Para uma melhor compreensão é recomendável que o leitor tenha conhecimentos básicos em física de plasmas e eletromagnetismo, uma vez que utilizamos estes conhecimentos para desenvolver este trabalho.

Palavras-chave: plasmas espaciais, simulação computacional, códigos PIC.

ABSTRACT

In this paper we invite the readers to discover the world of "computational experiments" (or computer simulation) of space plasmas. The computational experiment is similar to those traditionally made in laboratory. It is not purpose to solve all the physical problems that remain misunderstood, but provides a powerful tool to help physicists (scientists) in the search of their solution. We present some characteristics of computer simulations and examples where these simulations can be applied. For a better understanding of this paper it is recommended that the reader have a basic knowledge in plasma physics and electromagnetism, since we use these concepts to develop this work.

Keywords: space plasmas, Computer simulation, PIC codes.

1. Introdução

Tradicionalmente, o método científico envolve a interpretação mútua entre os experimentos de laboratório e as teorias feitas em "papel e caneta". No laboratório, os experimentos são realizados diversas vezes de forma controlada, permitindo ao cientista observar qual o comportamento do sistema e a descrição física mais adequada para explicar o fenômeno estudado. Por vários anos as teorias foram obtidas com sucesso desta forma. Nas últimas décadas, com o avanço tecnológico, um outro tipo de experimento foi introduzido para auxiliar os cientistas a modelar a forma complexa dos fenômenos que regem as leis da natureza. Este tipo de experimento é chamado de "experimento computacional" ou "simulação computacional".

O experimento computacional surgiu para diminuir a lacuna existente entre teoria e experimentos de laboratório. A idéia básica de um experimento computacional é simular de forma controlada o comportamento físico de um sistema complexo resolvendo um conjunto de equações matemáticas baseado em um modelo físico-matemático fundamental e bem aceito pela comunidade científica.

Embora similares às técnicas tradicionais usadas em laboratório nas quais os parâmetros físicos são variados de maneira controlada, os experimentos computacionais fornecem algumas vantagens: podemos obter informações detalhadas do sistema; os diagnósticos são não invasivos, ou seja, não perturbam o sistema; os efeitos físicos podem ser considerados ou não, permitindo identificar qual é o agente causador mais significante para um determinado fenômeno observado e principalmente, os experimentos podem ser reproduzidos de forma semelhante, sem influência de agentes externos. Apresentam ainda a vantagem de permitirem a investigação de fenômenos lineares, não lineares e dependentes do tempo. Em geral, nos experimentos tradicionais a obtenção dos diagnósticos perturba o sistema estudado.

Quanto maior for a lacuna entre a teoria e experimento, maior será a vantagem de se usar um experimento computacional. Em plasmas espaciais os experimentos são realizados através de observações de satélites, experimentos ativos no espaço e medidas a partir da Terra, ativas e passivas. Adquirir informações através destes experimentos é bastante onerosa e demanda muito tempo. Devido à complexidade dos fenômenos que ocorrem em plasmas espaciais e a dificuldade de se obter informações in loco, muitas vezes apenas em uma única região do espaço, os resultados das observações podem levar a uma interpretação e compreensão errônea do fenômeno. Em alguns casos, os experimentos computacionais podem reproduzir de forma global determinados fenômenos que não poderiam ser tratados analiticamente ou experimentalmente.

A natureza complexa dos problemas encontrados em física de plasma espacial tem motivado o interesse e desenvolvimento de códigos para o desenvolvimento das teorias de plasmas espaciais. Em adição, simulações numéricas também ajudam no desenvolvimento de sistemas computacionais de alta performance e previsões de comportamento de novos equipamentos e experimentos de física de plasma, como reatores de fusão termo-nuclear controlada (TOKAMAK) e outros.

A simulação computacional de plasma compreende duas grandes áreas de atuação, uma baseada na descrição cinética e outra na descrição de fluidos, Fig. 1. Num simulação baseada na descrição de fluidos resolve-se numericamente as equações considerando o plasma como um fluido magnetizado e que pode ou não estar eletricamente neutro. Um exemplo de aplicação da aproximação de fluido é as simulação magneto-hidrodinâmica (MHD) [1]. A simulação MHD considera o plasmas como um único fluido e é utilizada para tratar fenômenos de baixa frequência. Adotando os coeficientes apropriados, como conservação da massa, momento e energia, esta classe de simulação é usual para o estudo e entendimento da dinâmica macroscópica de um sistema.

Já as simulações cinéticas podem se constituir da resolução numérica das equações cinéticas, de Vlasov ou Fokker-Plank [1], ou da resolução numérica da equação do movimento de várias partículas carregadas, sujeitas a ação de campos eletromagnéticos, tanto os aplicados externamente quanto aos gerados pelo movimento das próprias partículas. Neste último caso, o método utilizado recebe o nome de simulação por partículas. As simulações via partículas permitem a interpretação de efeitos localizados, como efeitos não lineares por exemplo, instabilidades de ondas, difusão, aquecimento e aceleração de partículas em plasmas espaciais.

Os pioneiros em experimentos computacionais via partículas foram O. Buneman e J.M. Dawson que mostraram que a utilização de métodos apropriados tornam possíveis as simulações de sistemas simples, com algumas centenas de partículas, obtendo o comportamento coletivo de plasmas reais [2]. Com o surgimento de novos algoritmos e novas tecnologias na área computacional, foi possível expandir as técnicas de simulação de sistemas unidimensionais com poucas partículas para sistemas mais complexos envolvendo milhões de partículas em duas ou três dimensões.

Outra forma de simulação são os chamados códigos híbridos, no qual diferentes tratamentos são aplicados para diferentes componentes de um dado plasma. Estes códigos utilizam as chamadas "partículas-hidrodinâmicas", sendo as equações de fluido resolvidas por métodos de partículas [2].

Na realização de um experimento computacional permite-se que o sistema evolua a partir das condições iniciais de interesse. Os resultados obtidos podem ser comparados com os resultados experimentais e devem estar de acordo com a teoria, podendo ainda prever comportamentos em novos experimentos [3].

Dentre as técnicas de simulação computacional utilizadas para investigar plasmas, uma das mais antigas é a simulação por partículas, largamente utilizada para investigar fenômenos de plasma na magnetosfera terrestre, magnetosfera de outros planetas, vento solar e regiões de choque, foreshocks, de planetas e cometas [4]. Os modelos de simulação utilizando a aproximação de partículas podem ser uni, bi ou tridimensionais, dependendo do tipo de sistema estudado e dos diagnósticos necessários para a compreensão do fenômeno, podendo ser eletrostáticos, magnetostáticos e eletromagnéticos.

2. Importância das escalas de variação temporal e espacial na escolha dos modelos numéricos utilizados

Estimar as escalas de tempo e espaço para estudar um determinado fenômeno e escolher um modelo numérico que melhor se encaixe nestas escalas é um dos aspectos fundamentais das simulações computacionais. Para resumir a idéia do problema das escalas de tempo e de espaço vamos considerar um exemplo utilizando a Fig. (2.) [4]. A partir desta figura podemos observar as diferentes regiões onde cada modelo numérico é mais adequado, de acordo com as características físicas do sistema.

Em sistemas onde as dimensões de espaço e tempo são da ordem de L > 10⁴ km e T > 10² s a aproximação MHD é a mais apropriada, embora o regime híbrido ainda seja válido para grandes escalas, ele é mais aplicável em escalas menores do que as do regime MHD. No caso em que as escalas são compreendidas entre 10² < L < 10⁴km e 1 < T < 10²s a aproximação híbrida é mais adequada. Quando as dimensões do sistema estudado abrangem pequenos períodos de plasma e dimensões que são alguns comprimentos de Debye a aproximação por partícula é mais apropriada.

Caso a análise feita indique que o melhor mátodo a seguir é o de simulação por partículas, o procedimento é relativamente simples: segue-se as trajetórias das partículas no tempo e no espaço de acordo com a equação de movimento e as equações de Maxwell. Para realizar uma simulação, devemos nos lembrar que os computadores apresentam certos limites. As equações são resolvidas apenas em um número finito de pontos das região de interesse, chamados pontos de grade, que pode ser espacial ou temporal. Assim, vários métodos foram desenvolvidos para permitir a solução numérica das equações diferencias numa região bem definida, que devem ser discretizadas [5]. Em simulação via partículas, o método de leap-frog é utilizado para avançar a posição das partículas no tempo.

3. Modelo de simulação cinética

As simulações cinéticas têm sido aplicadas com sucesso no tratamento de problemas de física básica, nos quais, a função de distribuição das partículas desvia-se consideravelmente da distribuição Maxwelliana, quando ocorre aquecimento estocástico, aprisionamento de partículas ou ressonância onda-partícula [1].

As simulações MHD são, geralmente, aplicadas a problemas de grande escala, diretamente relacionados ao comportamento de dispositivos experimentais ou, por exemplo, para simular a estrutura global e dinâmica da magnetosfera em grande escala.

Em geral, os modelos de simulação cinética seguem o diagrama mostrado na Fig. 3., onde, primeiro determina-se o tipo de simulação e as aproximações que devem ser incluídas no modelo, o tipo de sistema estudado (eletrostático ou eletromagnético), a geometria, condições iniciais e condições de contorno. Estes três últimos são muito importantes, uma vez que, serão eles que incluirão todas as condições físicas do sistema estudado, como por exemplo, a região onde o sistema se encontra e as características desta região.

Após determinar corretamente o modelo numérico que melhor representa o sistema físico estudado, as simulações cinéticas via partículas utilizam as soluções das equações de Maxwell, baseadas em termos fontes auto-consistentes (densidade de plasma e correntes) [4], as quais são geradas a partir da posição e velocidade das partículas carregadas que representam o plasma submetido aos campos eletromagnéticos. Durante as simulações é necessário levar em conta os campos externos e os campos gerados pelo movimento das próprias partículas carregadas. Frequentemente, este processo acarreta no avanço das partículas num pequeno intervalo de tempo, ou passo temporal, Δt, para coletar os termos fontes que são utilizados para calcular os campos. Uma vez obtidos os novos campos, as partículas podem ser movidas novamente de forma que se obtenha os termos fontes, atualizados, e assim o processo é repetido tantas vezes quanto necessário. Este processo é ilustrado pelo ciclo dentro do quadro da Fig. 3., que é visto em mais detalhes na Fig. 4.

A simulação segue passo a passo, Δt, usando o método numérico que garanta estabilidade e precisão numérica suficientes. O passo temporal deve ser pequeno quando comparado com o período da onda estudada, ω_pΔt << 1, onde ω_p é a frequência local do plasma. O passo temporal deve ser pequeno o suficiente para que seja possível observar as variações dos fenômenos estudados.

Usualmente os campos são calculados na grade espacial, a partir da carga e da densidade de corrente. A precisão requer que o espaçamento da grade, Δx para uma simulação unidimensional, seja pequeno quando comparado com o menor comprimento de onda de interesse, kΔx << 1. A utilização de uma grade temporal e espacial inadequada pode introduzir um falso comportamento físico no sistema, minimizado através da escolha de parâmetros adequados para a simulação [2].

O ciclo se inicia em t = 0 com as condições iniciais apropriadas para a posição x_i e a velocidade v_i das partículas, que são representadas com o sub-índice i = 1,2,3,...N_total, onde N_total é o número máximo de partículas do sistema. Os campos são obtidos apenas na grade espacial, pontos discretos do espaço representados com o sub-índice j = 1,2,3,...P_total, E_j, B_j, onde P_total é o número de pontos da grade espacial. A partir da velocidade e posição das partículas calculam-se as densidades de carga e corrente na grade espacial, usando uma função de ponderação, a mesma utilizada para calcular os campos elétrico e magnético. A força, que é utilizada para mover as partículas, é obtida através de uma interpolação a partir da grade, novamente, através da função de ponderação; este ciclo segue durante toda a simulação [2,5-8].

4. Discretização espaço - temporal

Como mencionado anteriormente, vários métodos foram desenvolvidos para permitir a solução numérica de equações diferenciais que descrevem um sistema físico [5]. Um dos métodos mais utilizados para discretizar as equações diferenciais é o de diferenças finitas [5]. Uma boa aproximação para o cálculo da derivada primeira, df/dx, e derivada segunda, d²f/dx², de uma função f contínua, que varia lentamente em um intervalo Δx, pode ser obtida através das Eqs. (1) e (2) utilizando o método de diferenças finitas

As Figs. 5 e 6 apresentam o esquema geométrico para a obtenção das derivadas pelo método de diferenças finitas.

Em simulação com partículas o espaço e o tempo devem ser discretizados. A discretização espacial é introduzida por duas razões. A primeira razão surge devido a forma no qual força que atua sobre as partículas é calculada [8], porque, ao invés de calcular toda a contribuição da força Coulombiana a partir de todas as partículas, a força atuando em uma super-partícula é calculada pelas quantidades dos campos definidos nos pontos da grade próximo a ela. A segunda razão surge porque a super-partícula possui tamanho finito sobre uma certa região do espaço, de forma que uma resolução espacial menor que o tamanho da super-partícula é desnecessária e sem sentido [14]. Normalmente, o espaçamento da grade espacial varia de um a três comprimentos de Debye, este valor foi obtido através de experimentos numéricos [2,9]. A estabilidade numérica está diretamente relacionada a este fator, se for escolhido um espaçamento maior que três comprimentos de Debye surgirá instabilidade numérica [8] e assim, os resultados da simulação não terão significado físico.

Outro fator importante é o número de partículas por ponto na grade, quanto maior o número de partícula por ponto na grade, menor serão as flutuações numéricas relacionadas ao cálculo dos campos eletromagnéticos. Um estudo mais detalhado sobre instabilidade numérica pode ser encontrado nas Refs. [15] e [16].

A discretização temporal é um artifício inevitável em qualquer aproximação numérica para resolver equações diferenciais parciais. Uma questão que sempre surge é: qual valor atribuir para o passo no tempo? A resposta será: sempre o menor possível. A escolha deste valor deve ser de tal forma que exista estabilidade numérica durante a simulação. Uma condição que evita a instabilidade numérica e geralmente é utilizada em simulação computacional é a condição de Courant-Fredericks-Lewy (CFL) dada por

onde Δx é o espaçamento da grade, Δt o passo temporal e v_max a velocidade máxima que as partículas do sistema podem adquirir. Esta condição garante que em um passo temporal a distância percorrida pelas partículas com velocidade v_max não será maior que Δx.

Por exemplo, esta condição pode ser obtida através da propagação de uma onda eletromagnética no vácuo, cuja relação de dispersão é dada por ω² = k²c², onde ω é a frequência da onda, k o número de onda e c a velocidade da luz no vácuo. Supondo que esta onda eletromagnética possa ser representada por A(x,t) = A₀exp[i(k·x - ωt)] onde A(x,t) é a amplitude da onda num instante t e A₀ é a amplitude no instante t = 0, utilizando o método de diferenças finitas, centradas no espaço e no tempo a derivada de A(x,t) em relação a x será dada por

Comparando ΔA/Δx com a derivada parcial espacial ∂A/∂x podemos ver que o número de onda k pode ser substituído por K representado por

da mesma maneira a frequência ω pode ser substituída por uma frequência Ω dada por

Substituindo Ω e K na relação de dispersão para a onda eletromagnética anterior obtemos,

para o maior comprimento de onda k_max = π/Δx teremos

Se cΔt/Δx > 1, ω torna-se complexo, dando origem a instabilidade numérica, se cΔt/Δx = 1 o sistema será marginalmente estável e se cΔt/Δx < 1 o sistema será estável de forma que não será introduzida instabilidade numérica, esta é a conhecida condição CFL [9].

Um exemplo geométrico da condição CFL pode ser visto na Fig. 8. O caso (a), a curva contínua representada por Δx/Δt = v = c onde v é a velocidade de propagação da informação no sistema, e está sobre a curva tracejada que caracteriza a velocidade c, mantendo o sistema marginalmente estável. Caso (b), a curva Δx/Δt = v > c tornará o sistema numericamente instável, significando que as informações se propagarão com velocidade maior que c. Em (c), Δx/Δt = v < c o sistema torna-se estável, de forma que as informações se propagam com velocidade menor que c, em [17] pode se obter um tratamento detalhado de soluções numéricas de equações hiperbólicas.

5. Super-partículas

Podemos dizer que o método de simulação por partículas consiste em acompanhar o movimento de um grande número de partículas sob a ação de forças produzidas pelo próprio movimento das partículas, como resultado da interação entre elas, e/ou campos externamente aplicados. Em geral, a região de interesse contém um número extremamente alto de partículas. Por exemplo, no meio interplanetário (região próxima ao sol), o plasma apresenta densidade de partículas típica da ordem de 10⁴ a 10⁶ m^-3. Fazendo uma estimativa direta podemos ver que essa quantidade de partículas não é plausível na maioria dos computadores atuais pois, o número de operações aritméticas por partículas é bem elevado e assim, o tempo gasto para resolver as forças envolvidas no sistema, por passo de integração, seria muito grande, tornando as simulações impraticáveis.

Para contornar este problema foi criado o conceito de super-partícula, que é um modelo matemático que representa muitas partículas de um plasma real com tamanho finito, com sua carga distribuída sobre uma região finita do espaço [9].

As super-partículas foram inicialmente introduzidas por dois grupos de pesquisas para suprir flutuações estatísticas e colisões de pequeno alcance causadas por uma função delta [10-13,15]. Isto é, para partículas pontuais, numericamente, o potencial tende a infinito quando o raio da partícula tende a zero. Vários métodos foram desenvolvidos para determinar as densidades de cargas e correntes das super-partículas nos pontos da grade espacial. Um dos métodos de ponderação de área desenvolvido foi o ponto da grade mais próximo, Nearest Grid Point - (NPG), que pondera igualmente todas as partículas que se encontram em uma determinada distância do ponto da grade considerado [13]. Outros esquemas de ponderação das cargas para a grade são os denominados Nuvem em Célula, Cloud In Cell (CIC), ou Partícula em Célula, Particle In Cell (PIC), [11]. Estes dois métodos possuem uma pequena diferença: no método CIC, a posição da partícula determina o seu centro, e para o PIC a partícula é limitada pelas posições dos pontos da grade mais próximos, independente de sua posição na célula [2].

As super-partículas podem assumir qualquer forma, sendo que, a diferença fundamental entre elas é a maneira de acumular a carga nos pontos da grade a partir das suas posições. Entretanto, as mais utilizadas em simulação são denominadas de forma quadrada (1), forma triangular (2) e forma Gaussiana (3) [2,9,14], como apresentado na Fig. 7.

Uma vez que a super-partícula é idealizada para representar muitas partículas de um plasma real, a densidade de carga, massa e energia das super-partículas devem ser as mesmas das partículas reais [9], ou seja,

Densidade de carga → N_sQ_s = N_rQ_r;

Densidade de massa → N_sMs = N_rMr;

Densidade de energia → N_sk_BT_s = N_rk_BT_r.

Os sub-índices s e r representam a super-partícula e partícula real, respectivamente, N, Q, M, k_B e T são, densidade numérica, carga, massa, constante de Boltzmann e temperatura, respectivamente. Estas definições garantem que todos os parâmetros físicos básicos permaneçam os mesmos para as partículas reais e para as super-partículas durante a simulação. Por exemplo, razão carga-massa, frequência de plasma, frequência ciclotrônica, velocidade térmica e comprimento de Debye, , onde ω_pe é a frequência eletrônica de plasma, e v_the é a velocidade térmica dos elétrons. As propriedades físicas do plasma devem ser reproduzidas durante a simulação.

6. Código de simulação

O código numérico utilizado neste trabalho é o KEMPO 1D "Kyoto university's ElectroMagnetic Particle cOde" [9,14]. Este código pode ser obtido gratuitamente no endereço eletrônico (http://www.terrapub.co.jp/e-library/cspp/), com pequenas adaptações este código pode ser utilizado em um PC.

6.1. Metodologia usada para resolver as equações básicas

Como já foi dito antes, as equações básicas utilizadas são: a equação de movimento e equações de Maxwell.

A equação de movimento para uma partícula (identificada pelo subíndice i) com carga q e massa m é dada por

As equações básicas acima são escritas no sistema de unidade SI. Escrevendo F = q(E+v×B) e utilizando o método de diferenças finitas as equações acima podem ser escritas na forma

onde o índice sobrescrito representa a discretização temporal das equações.

A Fig. 9 ilustra como a evolução temporal da posição e velocidade de cada partícula, Eqs. (11) e (12), assim como a centralização no tempo do método de leap-frog. O método utilizado avança v^t e x^t para v^t+Δt e x^t+Δt. Podemos notar que v e x são defasados em t/2.

As partículas são avançadas da posição x usando a velocidade v_x. Em cada passo temporal Δt, a posição das partículas são avançadas duas vezes, cada uma delas por um passo Δt/2, conforme o modelo abaixo,

A velocidade das partículas é obtida integrando a equação de movimento (10) que através do método de diferenças finitas torna-se

Introduzindo novas variáveis v^- e v⁺ na forma

Substituindo as Eqs. (16) e (17) na Eq. (15) obtemos

Aplicando o produto escalar na Eq. (18) com (v⁺ + v^-), vamos obter

O termo do lado direito da Eq. (18) torna-se nulo na Eq. (19). Isso ocorre porque o vetor obtido no produto vetorial é perpendicular a ambos os vetores, assim, o produto escalar deste novo vetor com qualquer um dos anteriores será nulo.

Portanto, verificamos que o termo da força magnética é apenas responsável pela variação da direção de movimento das partículas, não alterando a magnitude da velocidade. A partir da Eq. (19) obtemos

A partir da Fig. 10 vemos que a Eq. (18) representa apenas uma rotação com ângulo θ dado por

A velocidade das partículas é obtida seguindo quatro passos conforme mostrado abaixo, este procedimento é chamado de método de Buneman-Boris [2]

onde E^t e B^t são os campos elétrico e magnético interpolados linearmente a partir dos valores obtidos nos pontos da grade. A quantidade (q/m)_s é a razão carga massa das espécies s, e v, E e B são quantidades vetoriais.

As componentes dos campos elétrico E ≡ (E_x, E_y, E_z) e magnético B ≡ (B_x, B_y, B_z), são obtidas somente sobre os pontos discretos de uma grade espacial, a partir da densidade de carga ρ e densidade de corrente J definidas nestes mesmos pontos da grade. Estas densidades são obtidas a partir das velocidades e posições das partículas,

onde J ≡ (J_x, J_y, J_z), c e μ₀ são o vetor densidade de corrente, velocidade da luz e permeabilidade magnética no vácuo, respectivamente. Em simulação, os valores para a permissividade ε₀ e permeabilidade μ₀ no vácuo podem ser definidos arbitrariamente, apenas com a condição que devem satisfaz a relação

Adotando um sistema uni-dimensional ao longo do eixo x, o campo elétrico E_x deve satisfazer a condição inicial dada pela equação de Poisson [14] e pela equação de Gauss

onde ρ e ε₀ são a densidade de carga e a permissividade elétrica, respectivamente. A equação de Poisson é resolvida somente como condição inicial do sistema, pois esta condição é satisfeita automaticamente se as Eqs. (26) e (27) forem resolvidas corretamente no tempo; a densidade de corrente J deve satisfazer a equação da continuidade e o campo magnético B_x, deve satisfazer à condição inicial

Esta condição garante que B_x seja constante no espaço e no tempo, uma vez que consideramos apenas o caso uni-dimensional e não temos os termos de B_x nas equações de Maxwell (26) e (27).

A densidade de corrente J e a densidade de carga ρ são obtidas a partir da equação de movimento para as partículas.

O campo elétrico é calculado integrando a equação abaixo

que para um sistema unidimensional pode ser escrita para cada uma de suas componentes como

O campo magnético é obtido integrando a equação abaixo

a equação acima pode ser escrita em suas componentes

A densidade de carga é calculada a partir das super-partículas, as quais são consideradas de forma quadrada, como mostra a Fig. 11. Uma super-partícula numa posição x_p tem uma distribuição de carga q/Δx no intervalo x_p - Δx/2 < x_p < x_p + Δx/2. Por outro lado, cada ponto da grade, em X_j, tem um território que cobre um intervalo, X_j - Δx/2 < X_j < X_j + Δx/2. Dessa forma, a carga q da super-partícula é dividida para os pontos adjacentes da grade, proporcionalmente à área compartilhada pelo pontos da grade. Numericamente, escrevemos Q(x_p - X_j)/Δx designado por ρ(X_j+1) e q(X_j+1 - x_p)/Δx designado por ρ(X_j).

A densidade de corrente é obtida a partir de equação da continuidade

que, através do método de diferenças finitas, torna-se

onde J é a componente J_x.

Semelhante ao cálculo da densidade de carga, a densidade de corrente é obtida satisfazendo a Eq. (38). Podemos supor duas situações desde que a partícula não se mova mais de um espaçamento da grade Δx em um passo temporal Δt, isto é

onde x_p é a posição da partícula no instante t. O primeiro caso é mostrado na Fig. 12, ambos x_p(t) e x_p(t + Δt) estão no mesmo ponto da grade, entre X_i e X_i+1. O segundo caso é mostrado na Fig. 13, neste caso x_p(t) e x_p(t + Δt) estão em pontos da grade diferentes.

No primeiro caso, a densidade de corrente J_i+1/2 em X_i+1/2 é obtida calculando a quantidade de carga passando no ponto X_i+1/2 em um passo temporal Δt

onde q_A e q_B são dados por

No segundo caso, o movimento das partículas contribui para a corrente nos pontos X_i+1/2 e X_i+3/2 que são dados por

Nas Figs. 12 e 13, consideramos as velocidades das partículas positivas, para velocidade negativa é necessário multiplicar o lado direito das Eqs. (40) e (42) por (-1).

6.2. Grade espacial

No código de simulação KEMPO são definidas duas grades, uma temporal e outra espacial. Para a grade espacial a posição j sobre a grade é definida por jΔx inteiro, (j = 1,2,3,...,N_x), e uma outra posição também é definida sobre a grade, de forma que o próximo ponto intermediário na grade seja dado por , (j = 1,2,3,...,N_x). Esta definição facilita a interpolação de diferenças finitas centrada no espaço para as derivadas espaciais das equações de Maxwell.

Dentro do código as componentes E_y, B_y, J_y e ρ são definidas sobre as j-ésimas posições da grade espacial, jΔx (j = 1,2,3,...,N_x), enquanto que as componentes E_z, B_z, J_x, J_z são definidas sobre as -ésimas posições da grade espacial, , (j = 1,2,3,...,N_x), sendo Δx o espaçamento entre os pontos da grade espacial. Na Fig. (14) é apresentado o esquema de como são calculadas a densidade, as componentes dos campos eletromagnéticos e a densidade de corrente para as partículas na grade espacial. Os quadrados marcam os pontos inteiros e os círculos representam os pontos intermediários da grade espacial.

6.3. Grade temporal

Da mesma forma que foi realizado para a grade espacial, a evolução temporal é realizada calculando as quantidades em pontos inteiros e meio pontos de grade temporal. Os pontos inteiros na grade temporal são determinados por t e os pontos intermediários por (n+1/2)Δt. O campo elétrico E e o campo magnético B são calculados pelo método de leap-frog sendo E calculado em passos temporais inteiros t e B calculado em passos temporais intermediários (n+1/2)Δt. No entanto, o avanço de Δt para o campo magnético é feito em dois passos de Δt/2; logo após o primeiro passo Δt/2, seu valor é utilizado para avançar a posição das partículas (calculadas em múltiplos inteiros de Δt)[9].

A posição x das partículas também é calculada nos pontos inteiros da grade temporal, enquanto a velocidade v é obtida nos pontos intermediários da grade temporal, pelo método de leap-frog. Semelhante ao que ocorre para os campos, a posição é avançada duas vezes com um meio passo temporal (Δt/2), para obter o valor intermediário e calcular J em um ponto intermediário da grade temporal. A densidade de corrente é calculada a partir da posição e da velocidade das partículas. A densidade de carga é calculada nos passos temporais inteiros e é utilizada para obter o campo elétrico. Este processo é ilustrado na Fig. (15) e se repete quantas vezes for necessário para obter os resultados da simulação.

6.4. Normalização e sistema de unidades

Para realizar as simulações não é necessário definir um sistema de unidades real, como CGS ou SI. O importante é definir as razões entre as quantidades do sistema utilizado, isto é, a razão entre o campo magnético da onda e o campo magnetostático ou a razão entre a energia cinética do plasma ambiente e a energia total do sistema, etc. Durante as simulações as quantidades físicas são normalizadas pelos parâmetros básicos do sistema tornando-se adimensionais. Entretanto, a seleção dos parâmetros básicos podem ser diferentes dependendo do modelo físico adotado. Para desenvolver um código de simulação aplicável a vários problemas de física e com um amplo sistema de variáveis considerando um conjunto de parâmetros básicos para obter as equações fundamentais. Estes parâmetros são:

1. Frequência angular (plasma, cíclotron, frequência da onda, etc.) ω_pi, Ω_c1, ω;

2. Comprimento do sistema L_x;

3. Razão carga-massa (q/m)_i;

4. Número de super partículas no sistema N_i;

onde o sub-índice i denota a i-ésima espécie de partículas. A frequência ciclotrônica é especificamente definida para a espécie 1, e está relacionada com o intensidade do campo magnético ambiente. Os valores destas quatro quantidades são dados de forma arbitrária, exceto que as razões entre as quantidades no mesmo sistema de unidade, ω_pi/Ω_c1, Δx/L_x ou (q/m)₂/(q/m)₁, sejam mantidas as mesmas das quantidades físicas reais. O número de super partículas não tem relação com a densidade numérica de partículas do plasma real e N₁ e N₂ são independentes um do outro.

As equações básicas são escritas de tal forma que serão idênticas às equações de movimento e Maxwell no sistema internacional. Os valores para a permissividade elétrica ε₀ e permeabilidade magnética μ₀ podem ser definidos arbitrariamente, apenas com a condição que devem satisfazer a relação

As quantidades físicas são calculadas seguindo as relações obtidas a partir das equações básicas. A frequência ciclotrônica para a espécie 1 e a frequência de plasma para as espécies i são dadas por

onde n_i é a densidade de partículas da espécie i que pode ser obtida por

A partir das Eqs. (44) e (45) obtemos as seguintes quantidades físicas

- Carga da partícula

- Massa da partícula

- Campo magnetostático

Podemos notar que a massa m_i e a carga q_i têm pouca importância no sistema "físico" computacional. Para realizar as simulações, os parâmetros mais importantes, fisicamente, são as densidade de massa n_im_i e a densidade de carga n_iq_i, que são dadas por

Os valores de ω_pi e (q/m)_i para as diferentes espécies devem ser levados em conta durante as simulações. Uma vez definidas as grades espacial e temporal, bem como o sistema de equações que serão utilizados na simulação juntamente com as condições de contorno adequadas, pode-se utilizar as simulações por partículas para estudar uma série de fenômenos de física dos plasmas como veremos a seguir.

7. Exemplos de simulação de partículas

Através das simulações por partículas podemos obter muita informação a respeito do comportamento físico de plasmas espaciais e de laboratório.

A seguir, apresentaremos alguns resultados obtidos com simulações por partículas, para fenômenos clássicos de física dos plasma, que ocorrem tanto em plasmas de laboratórios quanto em plasmas espaciais. Procuraremos ainda exemplificar a importância da definição dos parâmetros utilizados de forma a garantir a boa solução do problema.

No primeiro exemplo, apresentamos o resultado de simulação de uma oscilação eletromagnética em um plasma com o objetivo de visualizar a importância de garantir que a condição CFL seja satisfeita ao longo das simulações. A Tabela 1 apresenta os parâmetros utilizados na simulação, onde q/m é a razão carga massa, PCH é o ângulo entre o número de onda k e o campo B, v_pa e v_pe são as velocidades térmicas paralela e perpendicular das partículas, NX, NTIME e NP são número de pontos na grade espacial, número de passos temporal e o número de partículas de uma determinada espécie (NS = 1) envolvidas na simulação. Os outros parâmetros estão definidos no texto.

No caso (a) a condição de CFL é satisfeita (Δx > v_maxΔt) e os parâmetros utilizados são os mesmos da Tabela 1. Podemos perceber que durante a simulação a energia total do sistema (T) permanece constante. No caso (b) consideramos (Δx ≈ v_maxΔt) de forma que a condição de CFL seja satisfeita marginalmente. Neste caso, podemos perceber variação na energia total do sistema. Esta variação está relacionada com as flutuações numéricas, uma vez que não introduzimos nenhum termo fonte no sistema. Neste caso a variação da energia total é resultado das instabilidades numéricas. No caso (b) acima citado, não poderíamos confiar plenamente nos resultados das simulações, uma vez que as instabilidades numéricas poderiam modificar o comportamento físico do sistema.

Outro exemplo de simulação pode ser visto na Fig. 17, no qual mostramos as relações de dispersão para um plasma magnetizado e o espaço de fase das partículas. Neste exemplo consideramos os parâmetros da Tabela 1 com os seguintes parâmetros diferentes, NX = 256, NTIME = 1024, ω_pe = 2,0 e v_pa = v_pe = 1,0. Consideramos propagação paralela ao campo magnético externo k||B₀. Podemos observar os modos circularmente polarizado a direita (R) de alta e baixa frequência, cuja relação de dispersão é dada por

e o modo circularmente polarizado a esquerda (L), definido como

Na Fig. 18 mostramos uma simulação de instabilidade de dois feixes [1]. Os parâmetros utilizados são mostrados na Tabela 2, sendo v_D(a) e v_D(b) as velocidades de deriva para as espécies (a) e (b) respectivamente.

Podemos observar quatro quadros na Fig. 18: O primeiro quadro (1) é o espaço de fase no instante inicial da simulação. Podemos observar os dois feixes de elétrons com suas velocidades características. No quadro (2) observamos que os elétrons começam a mudar suas velocidades devido ao aprisionamento dos elétrons pelo potencial eletrostático, caracterizado pelo aumento da energia total do campo elétrico (E) (quadro superior esquerdo). Neste instante começa a formação de vórtices devido a interação dos dois feixes. O terceiro quadro (3) apresenta configuração final da simulação, no qual se observa a mistura dos dois feixes gerada pelo potencial eletrostático.

A combinação de experimentos computacionais, teoria e experimentos espaciais tem fornecido uma grande ferramenta para a compreensão de fenômenos de larga escala e processos cinéticos não lineares em plasmas espaciais.

Simulações cinéticas têm sido utilizadas para investigar fenômenos de plasma na magnetosfera terrestre, no vento solar e na região do arco de choque de planetas e cometas [4,14,16], entre outros.

Uma das maiores áreas de investigação de plasma espacial são ondas de choque não colisionais. O tipo mais comum de choque magnético são os choques rápidos, que surgem do escoamento das ondas magnetosônicas. Uma das regiões mais estudadas é o arco de choque terrestre, bem como a região de choque de outros planetas e cometas. O arco de choque é formado devido a interação do vento solar e a frente da magnetosfera em direção ao Sol. A interação da magnetosfera terrestre com o vento solar causa uma variação na velocidade do vento solar, aquecimento do plasma local e deflexão das partículas na direção do vento solar. Esta região envolve diversos fenômenos com várias escalas de tempo e espaço, e tem sido estudada cineticamente através de um grande número de técnicas, tais como, códigos implícitos [18], códigos explícitos [19], códigos por partículas e híbridos [20].

Outra região muito estudada em plasmas espaciais é a região acima do arco de choque da Terra conhecida como frente de choque. Esta região é caracterizada por partículas que escoam ou são refletidas a partir do arco de choque, retornando na direção do vento solar. A geometria desta região resulta em uma parte dominada pelos elétrons e outra dominada pelos íons que fluem em sentido contrário ao vento solar.

A frente de choque dos íons é uma região caracterizada pela presença de ondas de baixa frequência bem definidas (ULF - ultra low frequency) com períodos característicos de meio a um minuto. Acredita-se que as ondas observadas são associadas à instabilidades eletromagnéticas geradas por feixes de íons, que envolvem campos de deriva dos íons contrapropagantes alinhados ao vento solar [21]. O processo de geração de ondas envolve uma ressonância ciclotrônica do feixe de íons com uma onda eletromagnética. Neste processo os elétrons não têm efeitos significantes, sugerindo que os fenômenos relevantes podem ser estudados numericamente usando um modelo no qual os elétrons podem ser tratados passivamente mas, os efeitos cinéticos dos íons são necessários. O modelo de simulação híbrida tem se mostrado a melhor ferramenta para explicar tais processos na frente de choque [3].

Um exemplo clássico de instabilidade cinética é a excitação de ondas de Langmuir devida à instabilidade resultante da existência de uma "corcova" na cauda da função de distribuição de velocidades dos elétrons, Fig. 19. Este processo é utilizado para explicar as ondas de Langmuir observadas em diferentes regiões do plasma espacial e serve de base para a teoria das explosões de rádio solar tipo III. As simulações de partículas são utilizadas para investigar o processo de emissão de explosões de radio solar tipo III na região de reconexão magnética durante uma explosão solar solar flare [22].

8. Considerações finais

Neste trabalho apresentamos algumas características a respeito das simulações por partículas aplicadas à simulação de plasmas. Mostramos alguns critérios necessários para a realização das simulações computacionais e os passos que devem ser seguidos para a evolução temporal e espacial de um conjunto de partículas de um plasma magnetizado.

Ressaltamos o cuidado que deve ser tomado nas condições iniciais da simulação para evitar as instabilidades numéricas, que introduzem um falso comportamento físico do sistema. No caso onde não são incluídos termos fontes externos ao sistema as instabilidades numéricas são caracterizadas pela variação da energia total (T) do sistema, como pode ser observado na Fig. 16 caso (b).

Apresentamos alguns exemplos de diagnósticos que podem ser obtidos através das simulações computacionais, como as relações de dispersão de ondas que estão relacionadas com o índice de refração do meio no qual estas ondas se propagam. Outro tipo de diagnóstico muito utilizado é o espaço de fase. Através deste pode se acompanhar a trajetória das partículas e observar seu comportamento devido às interações com os campos eletromagnéticos. Outros diagnósticos podem ser obtidos através das simulações computacionais dependendo do tipo de fenômenos estudados.

As simulações computacionais não são e provavelmente nunca serão a solução para todos os problemas físicos estudados, mas representam uma poderosa ferramenta que pode auxiliar os pesquisadores a resolver diversos problemas de física que continuam sem respostas, principalmente problemas de sistemas não lineares e dependente do tempo.

Recebido em 9/4/2010; Aceito em 28/2/2011; Publicado em 30/3/2011

• [1] F.F. Chen, Introduction to Plasma Physics and Controlled Fusion (Plenum Press, New York, 1984), v. I.
• [2] C.K. Birdsall and A.B. Langdon, Plasma Physics via Computer Simulation (Institute of Physics Publishing, London, 1991).
• [3] M.V. Alves, Kinetic Simulation in Plasmas: A General View and Some Application INP-7171-PRE/3096 (1999).
• [4] D. Winske and N. Omidi, J. Geophysics Research 101, 287 (1996).
• [5] D. Potter, Computational Physics (John Wiley & Sons, 1973).
• [6] R.W. Hockney and J.W. Eastwood, Computer Simulation Using Particles (McGraw-Hill, Bristol, 1981).
• [7] T. Tajima, Computational Plasma Physics with Application to Fusion and Astrophysics (Addison-wesley Publishing Company, New York, 1989).
• [8] J.M. Dawson, Reviews of Modern Physics , 403 (1983).
• [9] H. Matsumoto and Y. Omura Computer Space Plasma Physics: Simulation Thecniques and Software (Terra scientific publishing company, Tokyo, 1993).
• [10] R.W. Hocney, Physics of Fluids 1826, (1969).
• [11] C.K. Birdsall and D. Fuss, J. of Computational Physics, , 141 (1997).
• [12] H. Okuda and C.K. Birdsall, Physics of Fluids , 2123 (1970).
• [13] R.L. Morse and C.W. Nielson, Physics of Fluids , 2418, (1969).
• [14] H. Matsumoto and T. Sato Computer Simulation of Space Plasmas (Terra Scientific Publishing Company, Tokyo, 1985).
• [15] A.B. Langdon and C.K. Birdsall, Physics of Fluids 2115, (1970).
• [16] H. Okuda and J.M. Dowson, Physics of Fluids , 408 (1973).
• [17] A.O. Fortuna, Técnicas Computacionais Para Dinâmica dos Fluidos (Editora da Universidade de São Paulo, São Paulo, 2000).
• [18] D. Krauss-Varban, Geophysical Research Letters , 3271 (1995).
• [19] P. Savoini and B. Lembege, J. of Geophysics Research-Space Physics , 6609 (1994).
• [20] D. Winske and M.M. Leroy, J. of Geophysics Research-Space Physics 7327 (1984).
• [21] S.P. Gary, J. of Geophysics Research-Space Physics , 4331 (1981).
• [22] J.I. Sakai, T. Kitamoto and S. Satio, The Astrophysical Journal , L157 (2005).

Datas de Publicação

• Publicação nesta coleção
  18 Abr 2011
• Data do Fascículo
  Mar 2011

Histórico

• Recebido
  09 Abr 2010
• Aceito
  28 Fev 2011