En física del plasma , el método de partículas en celdas ( PIC ) se refiere a una técnica utilizada para resolver una determinada clase de ecuaciones diferenciales parciales . En este método, las partículas individuales (o elementos de fluido) en un marco de Lagrange se rastrean en un espacio de fase continuo , mientras que los momentos de la distribución, como las densidades y las corrientes, se calculan simultáneamente en puntos de malla eulerianos (estacionarios) .
Los métodos PIC ya se utilizaban en 1955, [1] incluso antes de que estuvieran disponibles los primeros compiladores Fortran . El método ganó popularidad para la simulación de plasma a finales de los años 50 y principios de los 60 gracias a Buneman , Dawson , Hockney, Birdsall, Morse y otros. En aplicaciones de física de plasma , el método consiste en seguir las trayectorias de partículas cargadas en campos electromagnéticos (o electrostáticos) autoconsistentes calculados en una malla fija. [2]
Para muchos tipos de problemas, el método PIC clásico inventado por Buneman, Dawson, Hockney, Birdsall, Morse y otros es relativamente intuitivo y sencillo de implementar. Esto probablemente explica gran parte de su éxito, en particular para la simulación de plasma, para la cual el método generalmente incluye los siguientes procedimientos:
Los modelos que incluyen interacciones de partículas únicamente a través de los campos promedio se denominan PM (partícula-malla). Los que incluyen interacciones binarias directas se denominan PP (partícula-partícula). Los modelos con ambos tipos de interacciones se denominan PP-PM o P 3 M .
Desde los primeros días, se ha reconocido que el método PIC es susceptible a errores debido al llamado ruido de partículas discretas . [3] Este error es de naturaleza estadística y hoy en día sigue siendo menos comprendido que el de los métodos tradicionales de cuadrícula fija, como los esquemas eulerianos o semilagrangianos .
Los algoritmos PIC geométricos modernos se basan en un marco teórico muy diferente. Estos algoritmos utilizan herramientas de variedad discreta, formas diferenciales interpoladas e integradores simplécticos canónicos o no canónicos para garantizar la invariancia de calibración y la conservación de la carga, la energía-momento y, lo que es más importante, la estructura simpléctica infinitamente dimensional del sistema partícula-campo. [4] [5] Estas características deseadas se atribuyen al hecho de que los algoritmos PIC geométricos se basan en el marco teórico de campo más fundamental y están directamente vinculados a la forma perfecta, es decir, el principio variacional de la física.
En el ámbito de la investigación del plasma se investigan sistemas de diferentes especies (electrones, iones, neutros, moléculas, partículas de polvo, etc.). El conjunto de ecuaciones asociadas a los códigos PIC son, por tanto, la fuerza de Lorentz como ecuación de movimiento, resuelta en el llamado empujador o motor de partículas del código, y las ecuaciones de Maxwell que determinan los campos eléctrico y magnético , calculadas en el solucionador (de campos) .
Los sistemas reales estudiados suelen ser extremadamente grandes en términos de la cantidad de partículas que contienen. Para que las simulaciones sean eficientes o posibles, se utilizan las llamadas superpartículas . Una superpartícula (o macropartícula ) es una partícula computacional que representa muchas partículas reales; pueden ser millones de electrones o iones en el caso de una simulación de plasma o, por ejemplo, un elemento de vórtice en una simulación de fluido. Se permite reescalar el número de partículas, porque la aceleración de la fuerza de Lorentz depende solo de la relación carga-masa, por lo que una superpartícula seguirá la misma trayectoria que una partícula real.
El número de partículas reales correspondientes a una superpartícula debe elegirse de modo que se puedan recopilar suficientes estadísticas sobre el movimiento de las partículas. Si existe una diferencia significativa entre la densidad de las distintas especies en el sistema (por ejemplo, entre iones y neutros), se pueden utilizar para ellas proporciones independientes entre partículas reales y superpartículas.
Incluso con superpartículas, la cantidad de partículas simuladas suele ser muy grande (> 10 5 ), y a menudo el impulsor de partículas es la parte que más tiempo consume del PIC, ya que debe realizarse para cada partícula por separado. Por lo tanto, se requiere que el impulsor sea de alta precisión y velocidad y se dedica mucho esfuerzo a optimizar los diferentes esquemas.
Los esquemas utilizados para el transportador de partículas se pueden dividir en dos categorías: solucionadores implícitos y explícitos. Mientras que los solucionadores implícitos (por ejemplo, el esquema de Euler implícito) calculan la velocidad de la partícula a partir de los campos ya actualizados, los solucionadores explícitos utilizan solo la fuerza antigua del paso de tiempo anterior y, por lo tanto, son más simples y rápidos, pero requieren un paso de tiempo menor. En la simulación PIC se utiliza el método de salto de rana , un método explícito de segundo orden. [6] También se utiliza el algoritmo de Boris , que cancela el campo magnético en la ecuación de Newton-Lorentz. [7] [8]
Para aplicaciones de plasma, el método de salto adopta la siguiente forma:
donde el subíndice se refiere a cantidades "antiguas" del paso de tiempo anterior, a cantidades actualizadas del siguiente paso de tiempo (es decir, ), y las velocidades se calculan entre los pasos de tiempo habituales .
Las ecuaciones del esquema de Boris que se sustituyen en las ecuaciones anteriores son:
con
y .
Debido a su excelente precisión a largo plazo, el algoritmo de Boris es el estándar de facto para el avance de una partícula cargada. Se ha descubierto que la excelente precisión a largo plazo del algoritmo no relativista de Boris se debe al hecho de que conserva el volumen del espacio de fases, aunque no es simpléctico. El límite global del error de energía que normalmente se asocia con los algoritmos simplécticos sigue siendo válido para el algoritmo de Boris, lo que lo convierte en un algoritmo eficaz para la dinámica multiescala de los plasmas. También se ha demostrado [9] que se puede mejorar el impulso relativista de Boris para que preserve el volumen y tenga una solución de velocidad constante en campos E y B cruzados.
Los métodos más comúnmente utilizados para resolver las ecuaciones de Maxwell (o más generalmente, ecuaciones diferenciales parciales (EDP)) pertenecen a una de las siguientes tres categorías:
Con la FDM, el dominio continuo se reemplaza por una cuadrícula discreta de puntos, en la que se calculan los campos eléctricos y magnéticos . Las derivadas se aproximan entonces con las diferencias entre los valores de los puntos de la cuadrícula vecinos y, de este modo, las ecuaciones diferenciales parciales se convierten en ecuaciones algebraicas.
Utilizando el método de elementos finitos, el dominio continuo se divide en una malla discreta de elementos. Las ecuaciones diferenciales parciales se tratan como un problema de valores propios y, en un primer momento, se calcula una solución de prueba utilizando funciones de base que se localizan en cada elemento. La solución final se obtiene entonces mediante optimización hasta alcanzar la precisión requerida.
También los métodos espectrales, como la transformada rápida de Fourier (FFT), transforman las EDP en un problema de valores propios, pero esta vez las funciones base son de orden superior y están definidas globalmente en todo el dominio. El dominio en sí no está discretizado en este caso, permanece continuo. Nuevamente, se encuentra una solución de prueba insertando las funciones base en la ecuación de valores propios y luego se optimiza para determinar los mejores valores de los parámetros de prueba iniciales.
El nombre "partícula en celda" se origina en la forma en que las macrocantidades de plasma ( densidad numérica , densidad de corriente , etc.) se asignan a las partículas de simulación (es decir, la ponderación de partículas ). Las partículas pueden estar situadas en cualquier lugar del dominio continuo, pero las macrocantidades se calculan solo en los puntos de la malla, al igual que los campos. Para obtener las macrocantidades, se supone que las partículas tienen una "forma" determinada por la función de forma.
donde es la coordenada de la partícula y el punto de observación. Quizás la opción más fácil y utilizada para la función de forma es el llamado esquema de nube en celda (CIC), que es un esquema de ponderación de primer orden (lineal). Cualquiera que sea el esquema, la función de forma tiene que satisfacer las siguientes condiciones: [10] isotropía espacial, conservación de carga y aumento de la precisión (convergencia) para términos de orden superior.
Los campos obtenidos del solucionador de campos se determinan solo en los puntos de la cuadrícula y no se pueden usar directamente en el transportador de partículas para calcular la fuerza que actúa sobre las partículas, sino que se deben interpolar a través de la ponderación del campo :
donde el subíndice etiqueta el punto de la cuadrícula. Para garantizar que las fuerzas que actúan sobre las partículas se obtengan de forma autoconsistente, la forma de calcular las macrocantidades a partir de las posiciones de las partículas en los puntos de la cuadrícula y de interpolar los campos desde los puntos de la cuadrícula a las posiciones de las partículas también debe ser consistente, ya que ambos aparecen en las ecuaciones de Maxwell . Sobre todo, el esquema de interpolación de campos debe conservar el momento . Esto se puede lograr eligiendo el mismo esquema de ponderación para partículas y campos y asegurando la simetría espacial apropiada (es decir, sin fuerza propia y cumpliendo la ley de acción-reacción ) del solucionador de campos al mismo tiempo [10]
Como se requiere que el solucionador de campo esté libre de fuerzas propias, dentro de una celda el campo generado por una partícula debe disminuir al disminuir la distancia desde la partícula y, por lo tanto, las fuerzas entre partículas dentro de las celdas se subestiman. Esto se puede equilibrar con la ayuda de colisiones de Coulomb entre partículas cargadas. Simular la interacción para cada par de un sistema grande sería computacionalmente demasiado costoso, por lo que se han desarrollado varios métodos de Monte Carlo en su lugar. Un método ampliamente utilizado es el modelo de colisión binaria , [11] en el que las partículas se agrupan según su celda, luego estas partículas se emparejan aleatoriamente y, finalmente, los pares colisionan.
En un plasma real, pueden intervenir muchas otras reacciones, desde colisiones elásticas, como las colisiones entre partículas cargadas y neutras, hasta colisiones inelásticas, como la colisión de ionización entre electrones y partículas neutras, y reacciones químicas; cada una de ellas requiere un tratamiento independiente. La mayoría de los modelos de colisión que manejan colisiones entre partículas cargadas y neutras utilizan el esquema directo de Monte Carlo , en el que todas las partículas llevan información sobre su probabilidad de colisión, o el esquema de colisión nula [12] [13], que no analiza todas las partículas sino que utiliza la probabilidad máxima de colisión para cada especie cargada.
Como en todo método de simulación, también en PIC, el paso de tiempo y el tamaño de la cuadrícula deben elegirse adecuadamente, de modo que los fenómenos de escala de tiempo y longitud de interés se resuelvan adecuadamente en el problema. Además, el paso de tiempo y el tamaño de la cuadrícula afectan la velocidad y la precisión del código.
Para una simulación de plasma electrostático que utiliza un esquema de integración temporal explícito (por ejemplo, leapfrog, que es el más comúnmente utilizado), se deben cumplir dos condiciones importantes con respecto al tamaño de la cuadrícula y el paso de tiempo para garantizar la estabilidad de la solución:
que se puede derivar considerando las oscilaciones armónicas de un plasma unidimensional no magnetizado. La última condición es estrictamente necesaria, pero las consideraciones prácticas relacionadas con la conservación de la energía sugieren utilizar una restricción mucho más estricta donde el factor 2 se reemplaza por un número un orden de magnitud más pequeño. El uso de es típico. [10] [14] No es sorprendente que la escala de tiempo natural en el plasma esté dada por la frecuencia inversa del plasma y la escala de longitud por la longitud de Debye .
Para una simulación explícita de plasma electromagnético, el paso de tiempo también debe satisfacer la condición CFL :
donde , y es la velocidad de la luz.
Dentro de la física del plasma, la simulación PIC se ha utilizado con éxito para estudiar las interacciones láser-plasma, la aceleración de electrones y el calentamiento de iones en la ionosfera auroral , la magnetohidrodinámica , la reconexión magnética , así como el gradiente de temperatura de iones y otras microinestabilidades en tokamaks , además de descargas de vacío y plasmas polvorientos .
Los modelos híbridos pueden utilizar el método PIC para el tratamiento cinético de algunas especies, mientras que otras especies (que son maxwellianas ) se simulan con un modelo de fluido.
Las simulaciones PIC también se han aplicado fuera de la física del plasma a problemas de mecánica de sólidos y fluidos . [15] [16]
{{cite book}}
: CS1 maint: DOI inactivo a partir de mayo de 2024 ( enlace ){{cite journal}}
: CS1 maint: DOI inactivo a partir de mayo de 2024 ( enlace )