Físico-química de los fenómenos transitorios en descargas eléctricas

Memoria para optar al título de Doctor en Física Francisco José Durán Olivencia Departamento de Electrónica y Electromagnetismo, Facultad de Física UNIVERSIDAD DE SEVILLA  Francisco José Durán Olivencia Físico-Química de los Fenómenos Transitorios en las Descargas Eléctricas. Memoria para optar al título de Doctor en Física Licenciado en Física por la Universidad de Sevilla DIRECTOR: Francisco Pontiga Romero Profesor Titular en el Dpto. de Física Aplicada II Escuela Técnica Superior de Ingeniería de Edificación Universidad de Sevilla Co-responsable de la línea de investigación Físico-Química de las Descargas Eléctricas en Gases en el grupo EHD MGC: Electrohidrodinámica y Medios Granulares Cohesivos CO-DIRECTOR: Antonio Castellanos Mata Catedrático en el Dpto. Electrónica y Electromagnetismo Facultad de Física Universidad de Sevilla Responsable principal del grupo EHD MGC: Electrohidrodinámica y Medios Granulares Cohesivos INSTITUCIÓN: Dpto. Electrónica y Electromagnetismo, Facultad de Física Dpto. de Física Aplicada II, Escuela Técnica Superior de Ingeniería de Edificación Universidad de Sevilla Av. Reina Mercedes S/N, C.P. 41012, Sevilla, España. Junio, 2015. TESIS DOCTORAL Francisco José Durán Olivencia de los fenómenos en Memoria para optar al título de Doctor en Física Tesis inscrita en el Dpto. de Electrónica y Electromagnetismo Facultad de Física, Universidad de Sevilla, España Junio, 2015 Francisco José Durán Olivencia Licenciado en Física, Universidad de Sevilla DIRECTOR: Francisco Pontiga Romero Profesor Titular Dpto. de Física Aplicada II Universidad de Sevilla CO-DIRECTOR: Antonio Castellanos Mata Catedrático Dpto. de Electrónica y Electromagnetismo Universidad de Sevilla v Dedicatoria “Los locos abren los caminos que más tarde recorren los sabios.” CARLOS DOSSI, Locura. Dedicado a mis padres, a mi herman a, a mi hermano, a mi esposa Sara, a mi hijo Adrián y a todas aquellas personas que han confiado en mí y me han ayudado a desarrollar este trabajo. xi Prólogo En 1752 Benjamin Franklin desarrolló una teoría sobre los rayos (lightning), que supuso un hito histórico en la comprensión este tipo de fenómenos. Para él, los rayos podían entenderse como un flujo de electricidad a través del aire, lo que está esencialmente en consonancia con la interpretación actual. Desde un punto de vista fenomenológico no es difícil mostrar, entonces, que la distribución de las partículas (electrones e iones) implicadas en el flujo de tal densidad de corriente puede describirse mediante ecuaciones de continuidad con términos de reacción (creación/destrucción). De hecho, una buena parte de la literatura sobre las descargas eléctricas está basada en este enfoque. Por el contrario, para una imagen microscópica de las descargas (en el contexto de gases débilmente ionizados) se acepta que la ecuación de Boltzmann (BE) contiene toda la información estadística del sistema. En este trabajo, ambos enfoques se han sometido a un análisis profundo que delimita la validez de tales planteamientos. El primer capítulo introduce el concepto de descargas eléctricas en el contexto de los plasmas, y señala (de forma cualitativa) el efecto avalancha, el mecanismo de Townsend, el concepto de streamer y la distorsión del campo eléctrico local como precursores de los diferentes modos de propagación de las descargas. En la última parte de este capítulo se presenta la clasificación tradicional de los regímenes estacionarios (DC) de las descargas eléctricas a baja presión basada en la curva característica corrientetensión. Este tipo de descargas ha copado una parte significativa de las investigaciones en el pasado, y como resultado de ellas se han desarrollado muchas aplicaciones prácticas de uso cotidiano. Sin embargo, el estudio de las descargas transitorias despierta hoy en día tanto o más interés entre los investigadores, gracias al desarrollo de nuevas técnicas experimentales y computacionales, y a sus posibles aplicaciones tecnológicas. En el escenario previamente citado, respecto de las ecuaciones macroscópicas y la versión microscópica vía BE, la primera parte del capítulo 2 se centra en poner en contexto el ámbito de validez de la BE. Una limitación que no se tiene en consideración a menudo cuando se formula la BE en el marco de las descargas es el hecho de que esta ecuación está suscrita a la teoría cinética de gases. Esto significa que el concepto de colisión tiene una carácter localizado, tal y como sugiere la intuición física. Sin embargo, cuando se interpreta las descargas como un sistema de partículas cargadas Prólogo xii que interaccionan entre sí en presencia de campos electromagnéticos, esta imagen intuitiva puede fallar. Prueba de ello es que la BE se obtiene de una jerarquía BBGKI, donde están contenidos, por defecto, todos los tipos de correlaciones entre partículas. A partir de esta jerarquía es posible establecer dos límites que permiten reducir la comprensión de un sistema de muchas partículas a partir de la función de distribución de una sola partícula. Por estas razones, no conviene aceptar ad hoc el uso de la BE. De manera que para delimitar su validez, en primer lugar se analiza la situación en el contexto del equilibrio termodinámico local (LTE). En estas condiciones surge el concepto de longitud de Debye, lo que permite discernir la distancia a partir de la cual una partícula cargada experimenta un campo medio, resultado de la acción de todas las partículas restantes, sin atender a interacciones individuales. Esto supone una enorme simplificación que reduce el cálculo del campo eléctrico por medio de una ecuación de Poisson. Para el caso en el que el sistema se encuentra fuera del equilibrio, la aproximación de campo medio no resulta tan evidente. No obstante, la condición de régimen hidrodinámico permite establecer una relación de escalas (espaciales y/o temporales) que delimitan la validez de la BE. Así pues, la validez de las ecuaciones fluidas, obtenidas a partir de la BE, está igualmente vinculada a estos mismos límites. Esta es una idea extremadamente potente para actuar con limitadores físicos en el contexto numérico. Otra de los puntos clave en el modelado de las descargas que se aborda en este capítulo es la interacción con las superficies (electrodos y/o paredes del volumen de la descarga). A este respecto, hay que advertir que este trabajo solo plantea los principales mecanismos implicados en este tipo de interacciones, cuantificándolos y de este modo estableciendo un criterio para discriminar que tipo de procesos incorporar en el modelo. Se omiten en esta parte (obviamente, por la dimensión en complejidad) todos los efectos de superficies asociados a flujos tangenciales, desorción, recombinación, etc. El capítulo 3 expone el método de los elementos finitos (FEM) como herramienta para simular las descargas. No obstante, además de la revisión e introducción de las nociones elementales que subyacen a esta técnica, en este capítulo se exploran todas las dificultades inherentes a la formulación estándar en problemas de transporte. Así, se muestra como modificar las funciones de forma para estabilizar los problemas de convección dominante (en el contexto de la discretización espacial) y qué papel juega la discretización temporal en la estabilización de los términos de reacción. Para finalizar este capítulo se comparan los resultados obtenidos mediante elementos finitos (esquema estabilizado Crank-Nicholson (CN)+streamline upwind Petrov-Galerkin (SUPG)) con los presentados vía diferencias finitas (FD), con limitador flux–corrected transport (FCT), a partir del problema de Prólogo xiii la descarga corona en el régimen de pulsos de Trichel. Por último, el capítulo 4 trata el problema de los pulsos de Trichel bajo una descripción multiespecie. Este capítulo aúna todo lo visto en los capítulos 2 y 3, integrando un planteamiento físico, químico y numérico. La definición del tamaño mínimo de la malla muestra la conexión entre el modelo físico y el numérico. En cuanto a la descripción química se aborda la construcción de un modelo químico consistente reducido de 10 especies sobre la base de la distribución de pérdidas energéticas. Además se presta una especial atención a la hora de definir el mecanismo de ionización dominante. Así, no solo se analizan las vías directas, disociativa y no disociativa, sino también cuál es el papel de la ionización por pasos en el modelo y su relación con la ionización directa. El papel de este tipo de procesos en el desarrollo de streamers justifica la importancia de cuantificar la relación entre la ionización directa y por pasos. Por último, los resultados de las simulaciones muestran la compleja e intrincada dinámica que siguen las especies durante la evolución de la descarga. xv xvii Publicaciones Aunque una parte importante del contenido de este trabajo no ha sido publicado, algunas de las ideas y figuras aparecen en la siguiente lista de publicaciones: 1. Durán-Olivencia, F. J., Pontiga, F. & Castellanos, A. Multi-species simulation of Trichel pulses in oxygen. XXI Europhys. Conf. At. Mol. Phys. Ioniz. Gases (2012). 2. Durán-Olivencia, F. J., Pontiga, F. & Castellanos, A. 3D Multispecies Simulation in Low Pressure Oxygen. XXXI International Conference on Phenomena in Ionized Gases (2013). 3. Pontiga F., Durán-Olivencia, F. J. & Castellanos, A. Plasma-chemical simulation of negative corona near the inception voltaje. 66th Annual Gaseous Electronics Conference (2013). 4. Durán-Olivencia, F. J. Modeling and diagnostics of a multi-streamer surface DBD based jet for biomedical applications. COST Action, ECOST-STSM-MP1101-090114-040765 (2014). 5. Durán-Olivencia, F. J., Pontiga, F. & Castellanos, A. Numerical Modelling of Electrical Discharges. J. Phys. Conf. Ser. 490, 012209 (2014). 6. Durán-Olivencia, F. J., Pontiga, F. & Castellanos, A. Multi-Species Simulation of Negative Corona Pulses in Oxygen. 42nd Advanced Problems in Mechanics (2014). 7. Durán-Olivencia, F. J., Pontiga, F. & Castellanos, A. Multi-species simulation of Trichel pulses in oxygen. J. Phys. D: Appl. Phys. 47, 415203 (2014). 8. Pontiga F., Durán-Olivencia, F. J. & Castellanos, A. Plasma Chemistry Simulation of Trichel Pulses in Oxygen. Tlemcen (2014) xix xxvii Contenidos Dedicatoria v Abstract vii Resumen vii Prólogo xi Publicaciones xvii Agradecimientos xxi Símbolos xxxi Acrónimos y Abreviaturas xxxix Descargas Eléctricas __________________________________________ 1 1.1 Introducción .................................. 1 1.2 El efecto avalancha .............................. 4 Mecanismo de Townsend—4 El concepto de streamer—7 1.3 Modos estacionarios de una descarga DC .............. 10 1.4 Esta Tesis: Objetivos ............................ 13 Bibliografía 17 Un Modelo Microscópico _____________________________________21 2.1 Introducción ................................. 21 2.2 Una perspectiva microscópica ...................... 22 2.3 La aproximación de acoplamiento débil ............... 26 Colisiones entre partículas en las descargas—27 Factorización de 𝛼𝛽 (2)—30 2.4 La ecuación de Boltzmann ........................ 37 2.5 La aproximación fluida .......................... 39 2.6 Efectos de superficie ............................ 42 Emisión termoiónica—43 Fotoemisión (efecto fotoeléctrico)—45 Emisión por efecto campo—45 Emisión secundaria—46 Bibliografía 51 El método: FEM _____________________________________________55 3.1 Introducción ................................. 55 Contenidos xxviii 3.2 El método de los residuos ponderados ................. 57 Equivalencia entre las ecuaciones diferenciales y la formulación débil—57 Formulación Bubnov-Galerkin (método de Galerkin)—59 Formulación matricial—59 3.3 El problema de convección-difusión .................. 62 Representación estándar—62 Representación autoadjunta—65 Streamline Upwind Petrov-Galerkin (SUPG)—66 3.4 Discretización temporal .......................... 69 El algoritmo—70 3.5 Precisión y Estabilidad ........................... 72 3.6 Un test numérico .............................. 76 3.7 Algunas conclusiones finales ....................... 83 Bibliografía 85 Pulsos de Trichel en 𝟐 ______________________________________ 89 4.1 Introducción ................................. 90 4.2 Formulación del problema ........................ 93 Construcción de un modelo químico consistente—93 Descripción multifluida: una cuestión de escalas—100 Ecuaciones fluidas. Aproximaciones—102 Efectos de superficie. Interacción descarga-electrodo—104 4.3 Descripción del modelo ......................... 106 Sistema de ecuaciones de continuidad—108 Condiciones de contorno—110 Condiciones iniciales—112 4.4 Aproximación numérica ......................... 112 4.5 Evolución de los pulsos de Trichel .................. 114 El primer pulso de Trichel—114 Pulsos de Trichel regulares—121 4.6 Conclusiones ................................ 131 Apéndice 133 Bibliografía 135 Conclusiones ______________________________________________ 141 xxix de xxxi Símbolos (1) diferencia potencial entre electrodos (2) volumen del sistema distancia entre electrodos 0 campo eléctrico externo 0 intensidad de corriente inicial debida a la ionización natural (1) primer coeficiente de Townsend (recorrido medio para que el que tiene lugar un proceso de ionización) (2)factor de peso que actúa sobre el término estabilizador en la formulación streamline upwind Petrov-Galerkin (3) (superíndice) ∗ : coeficiente de excitación (1) densidad numérica del medio (2) número total de partículas 𝛼 densidad numérica de la especie presión 𝑠𝑒 segundo coeficiente de Townsend generalizado energía a de ionización 0 trabajo de extracción de un electrón en el metal función escalón unitario 𝑐𝑎𝑡 intensidad de corriente total en el cátodo 2 coeficiente generalizado para la captura electrónica (recorrido medio para que un proceso de este tipo tenga lugar) 3 coeficiente generalizado para la recombinación electrónica (recorrido medio para que un proceso de este tipo tenga lugar) 𝑒𝑓𝑓 valor efectivo del primer coeficiente de Townsed ( 𝑒𝑓𝑓 , donde 2 3 ) 𝑡 tensión umbral para el mecanismo de ruptura o ignición de Townsend 𝐴 diferencia de potencial aplicada a los electrodos Símbolos xxxii 𝑐 intensidad de corriente en una descarga tipo corona (superíndice) 𝜈 número de especies que componen un sistema de partículas 𝑁 vector de posición para un sistema de partículas 𝑁 vector velocidad para un sistema de partículas función de distribución 𝑖 vector de posición para la partícula i - ésima 𝑖 vector velocidad para la partícula i - ésima 𝛼 (especie cargada) contribución de arrastre debida al campo eléctrico 𝑖 vector de aceleración total que sufre la partícula i velocidad de la luz en el vacío tiempo derivada sustancial o material 𝑖 carga de la partícula /especie i 𝑖 masa de la partícula i 𝑖𝑗 potencial eléctrico entre las partículas i y j ( 𝑠 ) función de distribución reducida para s partículas campo eléctrico campo magnético (superíndice) ( ext ) referido a los campos aplicados externamente a un sistema de partículas cargadas 𝛼 función de distribución reducida para una partícula de la especie 𝛼𝛽 función de distribución por pareja para las especies y . 𝛼𝛽𝛾 función de distribución para una terna de partículas de las especies , y 𝛼𝛽 función de correlación por parejas 12 distancia relativa entre las partículas y 𝑐 radio de correlación 𝐷 longitud de Debye parámetro del plasma 𝐵 constante de Boltzmann Símbolos xxxiii temperatura 𝐻 𝐻 escalas características de la relajación macroscópica : espacial y temporal, respectivamente 𝑚𝑝𝑓 𝑚𝑓𝑝 escalas características de la fase cinética : espacial (recorrido libre medio) y temporal, respectivamente campo eléctrico medio (promedio en el espacio de las fases) creado por todas las partículas cargadas campo magnético medio (promedio en el espacio de las fases) creado por todas las partículas cargadas 𝛼 𝑐 término de colisión en la ecuación de Boltzmann vector densidad de corriente macroscópica sección eficaz de colisión ángulo sólido 𝛼𝛽 velocidad relativa entre una partícula de la especie y otra de la especie 𝛼 𝑐 promedio temporal del término de colisión de Boltzmann 𝑐 escala espacial característica para una colisión binaria (1) parámetro de gas diluido (2)probabilidad de emisión secundaria por impacto (3)(subíndice) 𝑝 : factor de eficiencia para los fotoelectrones emitidos en la superficie del cátodo ángulo de dispersado en una colisión binaria 𝑅 masa reducida 0 permitividad eléctrica del vacío parámetro de impacto 𝑒𝑛 𝑒𝑒 frecuencias de colisión electrón - partícula neutra (medio) y electrón - electrón, respectivamente 𝑚 frecuencia de colisiones elásticas 𝛼 velocidad media de la especie 𝛼 tensor de presiones para la especie 𝛼 fricción neta debida al efecto de las colisiones para la especie 𝛼 flujo de energía térmico o flujo de calor 𝛼 término de reacción (ec. de continuidad de la especie Símbolos xxxiv ) 𝐸 , 𝛼 flujo de calor asociado al efecto de las colisiones para la especie (1) potencial químico (2)coeficiente de absorción (radiación) del medio (3) (subíndice ) 𝛼 : movilidad de la especie 𝐹 energía de Fermi frecuencia de la luz ( ) constante de Planck ( ) 𝐹𝑁 𝐹𝑁 constantes para la relación de Fowler - Nordheim 𝑀 electrón en el metal X ∗ energía del estado excitado que captura el electrón del metal en un proceso de neutralización de Auger defecto de energía X + − 𝑀 energía de interacción del ion con la superficie metálica a una distancia de ésta X − 𝑀 energía de interacción del átomo con la superficie metálica a una distancia de ésta coeficiente de difusión 𝑠 normal a la superficie (apuntando desde el volumen a la superficie) (subíndice) 𝑝 Referido a la pared (1) interior del dominio (2)ángulo sólido definido por un disco de carga (método de los discos) desde el punto del cátodo situado en el eje de simetría frontera del dominio región que encierra el problema D contorno del tipo Dirichlet N contorno del tipo Neumann función arbitraria utilizada para la ponderación en el método de los residuos ponderados (discretización espacial) conjunto de funciones (solución) admisibles en la formulación débil del problema Símbolos xxxv conjunto de variaciones admisibles en la formulación débil del problema 1 espacio de Sobolev de grado 1 y orden 2 2 ó 2 funciones de cuadrado integrable sobre el dominio ( ) o el contorno ( ), respectivamente 𝑚 funciones de clase m , esto es, con derivada m - ésima continua ℎ aproximación de dimensión finita del espacio funcional ℎ aproximación de dimensión finita del espacio funcional 𝑘 espacio de polinomios de grado k ℎ y ℎ funciones solución y variación/test, respectivamente, definidas en sendos espacios aproximados 𝐸 número total de elementos (discretización espacial) número total de nodos 𝑖 función de forma asociada al punto i de la malla 𝑖 valor noda l (punto i de la malla) matriz de amortiguamiento ( damping matrix ) matriz de rigidez ( stiffness matrix ) matriz de carga ( load matrix ) 𝑖 tamaño de la malla en el elemento i número de Péclet función arbitraria utilizada para la ponderación en el método de los residuos ponderados (discretización temporal) residuo matriz tangente o Jacobiana operador (numérico) asociado a la ecuación diferencial en cuestión factor de amplificación número de Courant número de difusión 𝐷 error de difusión asociado al esquema numérico Φ error de dispersión (o error en el fase) asociado al 1 Descargas Eléctricas Contenido 1.1 Introducción .................................. 1 1.2 El efecto avalancha .............................. 4 Mecanismo de Townsend—4 El concepto de streamer—7 1.3 Modos estacionarios de una descarga DC .............. 10 1.4 Esta Tesis: Objetivos ............................ 13 Bibliografía 17 1.1 Introducción Originalmente, el concepto de la descarga eléctrica aparece ligado con la descarga de un condensador en un circuito eléctrico. El condensador se introduce separando dos electrodos metálicos e incorporando entre éstos un material dieléctrico (un gas en el caso más simple). En estas condiciones, si la tensión aplicada a los electrodos es lo suficientemente elevada, se produce la ruptura dieléctrica del gas, formándose un nuevo estado en el que el gas está ionizado. Este proceso, en el que un gas inicialmente neutro (que actuaba como aislante) se convierte en conductor, se conoce como breakdown, ruptura o ignición. Cuando esto ocurre, el circuito se cierra y en consecuencia el condensador se descarga. De ahí que el término de descarga se extienda de manera natural a todos aquellos procesos en los que surgen Capítulo “ The variety and complexity of the electrical phenomena which occur when matter is present in the electric field are in marked contrast to the simplicity of the phenomena when the ether alone is involved.” J.J. THOMSON, The Discharge of Electricity Through Gases, 1898. 2 Capítulo 1 | Descargas Eléctricas este tipo de corrientes, derivadas de una etapa de ionización previa. Actualmente, las descargas eléctricas se consideran una rama más de la física de los plasmas. El término “plasma” fue acuñado por Irving Langmuir1 en 1928 para referirse a la zona cuasi neutra que aparecía detrás de las capas límites formadas en las proximidades de los electrodosa. Aunque no existe una versión definitiva de la razón que impulsó a Langmuir para adoptar esta terminología, Rogoff2 sugiere que la elección estuvo motivada por el significado griego de la palabra, en alusión al carácter moldeable de algo. Es precisamente en este contexto en el que a mediados del siglo XIX se adoptó por primera vez dicha palabra griega. Entonces, el fisiólogo checo Jan Evangelista Purkinje introdujo el término plasma para referirse al fluido remanente tras eliminar todo el material corpuscular de la sangre. El concepto de plasma en la física abarca las peculiaridades de la materia ionizada, motivo por el cual se conoce como el cuarto estado de la materia. La presencia de portadores de carga libre hace de los plasmas sistemas conductores (en términos eléctricos), con una fuerte predisposición a la interacción interna y muy susceptibles a la presencia de campos electromagnéticos. La esencia de todas las peculiaridades características de los plasmas puede reducirse a la comprensión de la dinámica de un conjunto de partículas cargadas en el seno de campos electromagnéticos. Sin embargo, el marco teórico necesario para ello manifiesta la complejidad y el carácter interdisciplinar que conlleva afrontar este tipo de fenómenos con un mínimo de rigurosidad. Aun a día de hoy, después de décadas de investigación las descargas eléctricas, muchas de las cuestiones fundamentales continúan abiertas. Desde luego, una parte importante del esfuerzo realizado históricamente se ha centrado en descifrar aquellos aspectos cualitativos que permiten explicar la taxonomía de las descargas eléctricas. Como resultado de ello, muchas de las características cruciales para entender la evolución espacio-temporal de las descargas empiezan a comprenderse en todas sus facetas. Un elemento fundamental como precursor de la ruptura es el efecto avalancha que experimentan los electrones. Este proceso se desarrolla en aquellas zonas donde existe un campo eléctrico relativamente elevado que conduce a un proceso en cadena vía ionización. Habitualmente, la avalancha suele iniciarse a partir de un pequeño número de electrones semillas que aparecen en la región entre electrodos de forma accidental, originados por la radiación a “Except near the electrodes, where there are sheats containing few electrons, the ionized gas contains ions and electrons in about equal numbers, so that the resultant space charge is very small. We shall use the name plasma to describe this region containing balanced charges of ions and electrons.” 1.1 Introducción 3 cósmica u otra fuente natural. Sin embargo, a fin de obtener unas condiciones iniciales reproducibles, los electrones semillas también pueden crearse artificialmente por medio de fotoelectrones emitidos en el cátodo al irradiar éste con una fuente ultravioleta. La ruptura o breakdown es un fenómeno que esencialmente tiene lugar a partir de un cierto umbral energético (threshold). Por este motivo, en aquellas zonas donde el campo eléctrico es elevado y los electrones alcanzan la energía umbral de ionización del gas (medio), las propiedades de éste se verán afectadas significativamente. Por otra parte, en aquellas regiones donde no se alcanza dicho umbral energético, el estado del medio no variará de manera notable, al menos en primera instancia. A posteriori, estas zonas desempeñarán un papel fundamental como sumidero electrónico junto con el propio electrodo (ánodo). El caso de los gases electronegativos es especialmente interesante dado que la carga iónica negativa acumulada y el barrido de ésta son dos de los parámetros que condicionan el desarrollo de nuevos frentes de avalancha en las proximidades del electrodo. El proceso realimentado o multiplicativo de las avalanchas desencadenado por la presencia de un campo eléctrico externo que acelera los electrones se conoce como mecanismo de Townsend. Este mecanismo de ignición permite explicar el desarrollo de las descargas principalmente a presiones bajas y moderadas. En el caso de descargas automantenidas debe existir además una realimentación entre los electrones secundarios emitidos por el cátodo, los iones positivos y el campo eléctrico. A presiones superiores (del orden de la atmosférica) es necesario introducir un efecto añadido que describa la aparición de canales de ionización conocidos como streamers, que pueden evolucionar para convertirse en leaders y llegar incluso a comunicar ambos electrodos (sparks). La distorsión del campo local creada por la propia descarga juega un papel fundamental en el desarrollo de los streamers, aunque también hay implicados otros factores en la realimentación como la ionización por pasos o la propia radiación emitida por el streamer. El concepto de streamer junto con el mecanismo de Townsend y la distorsión del campo eléctrico representan los principales procesos utilizados para explicar la evolución de las descargas. Este último puede llegar a ser determinante en el efecto de multiplicación por avalancha por medio de la carga espacial acumulada en el transcurso de la descarga. Sin embargo, la acción combinada del campo eléctrico externo junto con el campo creado por el resto de partículas cargadas constituye el principal problema, aunque no el único, cuando se requiere una descripción cuantitativa. Este último campo depende de las posiciones de las partículas cargadas en cualquier punto del espacio, por lo que puede variar rápidamente en el espacio y tiempo. Como 4 Capítulo 1 | Descargas Eléctricas resultado de ello, se obtienen un conjunto de ecuaciones de movimiento que forman un complejo sistema de ecuaciones no lineales fuertemente acopladas entre sí. En este contexto, la aproximación fluida (§2.5) se ha utilizado históricamente para dar una descripción de grano grueso que permitiese transformar las ecuaciones en un conjunto “más asequible”. No obstante, algunas de las hipótesis implícitas en el tratamiento fluido requieren cierta cautela y pueden fallar completamente en según qué zonas de la descarga (esta discusión se seguirá en más detalle en §4.2). 1.2 El efecto avalancha Mecanismo de Townsend Consideremos dos electrodos planos separados por una distancia d y entre los que se aplica una diferencia de potencial V. Si el cátodo se expone a radiación ultravioleta de frecuencia adecuada, la emisión de fotoelectrones por el metal, generará una corriente inicial, 0 , asociada con el arrastre de éstos hacia el ánodo. El desplazamiento de este primer frente electrónico puede generar nuevos electrones por medio de la ionización de las moléculas del gas (véase Fig. 1.2-1). El primer coeficiente de Townsend, , caracteriza este proceso de ionización por avalancha indicando el número de electrones producidos (vía ionización directa) por unidad de longitud, 𝑥 𝑒 𝑒 . En la literatura3–5 suele encontrarse en términos reducidos, bien como o donde y representa la presión y densidad del medio Fig. 1.2 - 1 : Esquema ilustrativo de la evolución de la avalancha electrónica en un campo uniforme donde la carga espacial no llega a distorsionar el campo externo, 0 . 1.2 El efecto avalancha 5 respectivamente. En cada proceso de ionización, además de un nuevo electrón se genera un ion positivo que migrará hacia el cátodo. Como resultado de ello, por cada electrón primario que abandona en el cátodo se recolectan 𝛼𝑑 iones positivos en dicho electrodo. Por otra parte, el segundo coeficiente generalizado de Townsend, 𝒔𝒆, describe los diferentes procesos implicados en la emisión secundaria de electrones en la superficie del cátodo, tales como: 1. Emisión secundaria debido al impacto iónico en el cátodo (proceso de neutralización de Auger6). 2. Emisión de fotoelectrones en el cátodo debido al efecto de la radiación creada en el transcurso de la descarga (fotoemisión). 3. Emisión de electrones secundarios asociados al impacto de estados excitados y metaestables sobre la pared del cátodo (proceso de desexcitación de Auger6). A menudo suele utilizarse una relación experimental3 (con una precisión en torno al 50% en superficies metálicas “limpias”) como primera aproximación de este coeficiente, 𝑠𝑒 0 0 1.2 . 1 donde representa la energía de ionización del medio, 0 es el trabajo de extracción de un electrón y 0 describe la función escalón unitario. Típicamente los valores de 𝑠𝑒 suelen estar comprendidos en un rango7 entre y . El efecto de la avalancha inicial junto con las posibles avalanchas secundarias conduce a una contribución electrónica a la corriente total en el cátodo de acuerdo con la relación: 𝑐𝑎𝑡 0 𝑐𝑎𝑡 𝑠𝑒 𝛼𝑑 1.2 . 2 donde 0 representa la contribución de los fotoelectrones primarios y el segundo término, 𝑐𝑎𝑡 𝑠𝑒 𝛼𝑑 , justifica la contribución debida a la emisión secundaria en el cátodo. La corriente electrónica en el ánodo coincide con la corriente externa del circuito, dada la ausencia de contribución iónica en dicho electrodo. En tal caso, la corriente total en el ánodo resulta, 𝑐𝑎𝑡 𝛼𝑑, lo que conduce a la conocida fórmula de Townsend: 0 𝛼𝑑 𝑠𝑒 𝛼𝑑 1.2.3 La validez de esta expresión está condicionada a que las densidades de electrones e iones positivos sean suficientemente bajas, de forma que los 6 Capítulo 1 | Descargas Eléctricas procesos de recombinación electrón-ion puedan despreciarse. En estas circunstancias, los electrodos constituyen las principales vías para la eliminación de la carga. No obstante, a medida que la intensidad de corriente aumenta, la recombinación irá ganando importancia, y puede llegar a ser un factor limitante en el desarrollo de las sucesivas avalanchas secundarias. Para incorporar el efecto de procesos como la captura, 2, o recombinación electrónica, 3 , en la relación 1.2.3 basta considerar un valor efectivo para el primer coeficiente de Townsend 𝑒𝑓𝑓 , siendo 2 3. La expresión 1.2.3 proporciona el valor de la corriente para una presión, , y campo eléctrico, , dados, siempre que sea menor que un cierto valor crítico, 𝑡. Este último queda definido por la condición 𝑠𝑒 𝛼𝑑 1.2 . 4 Puede parecer que en el proceso límite con tendiendo a 𝑡, el valor de la corriente crece indefinidamente. Sin embargo, en la práctica, esta aproximación está limitada por la estabilidad de la propia fuente de alimentación. En un caso particular, fijada la distancia entre electrodos, , queda determinado el valor del campo externo 0, que a su vez fija el valor de las cantidades y 𝑠𝑒. Siguiendo esta cadena de relaciones entre 0 (tomando 𝑠𝑒 constante), la tensión umbral 𝑡 que satisface la relación 1.2.4 es aquella para la cual 𝑡 𝑠𝑒 1.2.5 Así la tensión umbral conduce a un valor del primer coeficiente de Townsend (relacionado con la tasa de ionización) para el que se cumple que 𝑠𝑒 𝛼𝑑 . Es decir, dado que por cada electrón primario que abandona el cátodo se obtienen 𝛼𝑑 iones que impactan sobre dicho electrodo, la condición anterior (ec. 1.2.4 y/o 1.2.5 ) expresa que por cada electrón del frente primario se produce un electrón por emisión secundaria en el cátodo. Por este motivo 1.2.4 no solo refleja la condición para el mecanismo de ruptura o ignición de Townsend, sino que define la transición entre descargas automantenidas, , y no automantenidas . Tanto en modelos teóricos como numéricos, el valor suele ajustarse por medio de la relación experimental propuesta por Townsend8,3,7, − 𝐵𝑝 𝐸 1.2.6 donde y son constantes de ajuste experimentales. En este caso 1.2.6 sustituida sobre 1.2.5 conduce a una relación para la tensión umbral, 1.2 El efecto avalancha 7 𝑡 𝑡 1.2.7 donde 𝑠𝑒 Las curvas 𝑡 descritas por la relación 1.2.7 se conocen como curvas de Paschen. Los valores para constantes y pueden deducirse a partir del ajuste de 1.2.6 para el coeficiente de ionización8,3,7. En Fig. 1.2-2 se ha representado el caso particular del oxígeno a . El concepto de streamer El mecanismo de Townsend describe apropiadamente una descarga estacionaria en el rango de presiones medias y bajas ( ), dentro de un sobrevoltaje moderado. Incluyendo el efecto de la distorsión del campo eléctrico en las proximidades del cátodo es posible explicar, incluso, la transición desde la descarga oscura de Townsend hacia una descarga de tipo glow (véase §1.3). Sin embargo cuando los valores de y el sobrevoltaje son elevados, el desarrollo de la descarga es mucho más rápido que el predicho por el efecto de sucesivas avalanchas iniciadas por medio de las emisión secundaria en el cátodo. Es decir, no hay tiempo suficiente para la multiplicación por avalancha desde el cátodo, incluso cuando se considera la Fig. 1.2 - 2 : Curva de Paschen para el 2 a (el valor de las constantes se ha tomado de Lieberman y Lichtenberg7). 8 Capítulo 1 | Descargas Eléctricas fotoemisión. En este contexto surge el concepto de streamer8–10. El streamer representa un estrecho canal ionizado que aparece entre los electrodos como consecuencia de una primera avalancha muy intensa. Este es el canal que “seguirá” la carga positiva hacia el ánodo. A diferencia del mecanismo de Townsend, vinculado a la emisión secundaria en el cátodo, el concepto de streamer explica precisamente aquella región para la cual el desarrollo de la descarga es independiente del material usado en el cátodo (la aparición de streamers puede tener lugar incluso cuando se emplean electrodos de aleaciones con coeficientes de emisión secundaria del orden de, 𝑠𝑒 −8 ). En este caso los electrones secundarios que inician las sucesivas avalanchas se deben a los fotones que caen dentro de este canal ionizado. Estos fotones son emitidos por la nube de átomos excitados que resultan tanto de la primera como de las sucesivas avalanchas secundarias. La distorsión del campo eléctrico local en el frente del streamer, a diferencia de la situación en el cátodo, no llega a generar una región de ionización lo suficientemente extensa. Por este motivo, el papel de la interacción entre el frente del streamer y la radiación generada por la propia descarga es crucial. Por tanto, ahora el ciclo de realimentación positiva se construye por medio de la interacción entre la distorsión del campo eléctrico local, estados excitados y procesos del tipo ionización por pasos en los que electrones de baja energía resultantes de la ionización desempeñan un papel fundamental (el proceso de ionización por pasos puede llegar a ser entre 3 y 4 veces más rápido que la ionización directa, véase §4.2). La transición desde la avalancha al streamer ocurre cuando el campo eléctrico local se hace comparable al campo eléctrico impuesto externamente Fig. 1.2 - 3 : (a) Evolución de la descarga entre electrodos. (b) Distorsión del campo eléctrico en la avalancha. ( a ) ( b ) 1.2 El efecto avalancha 9 (véase Fig. 1.2-3). Si la distancia entre electrodos es relativamente corta, esta transición solo ocurre en las proximidades del ánodo. Este tipo de streamers que crecen en la dirección del cátodo se conocen como streamers dirigidos al cátodo (cathode-directed streamers) o streamers positivos. Cuando la separación entre electrodos y el sobrevoltaje son suficientemente grandes, el streamer puede ocurrir lejos del ánodo. Estos streamers se conocen como Fig. 1.2 - 4 : Esquema representativo de la evolución de un streamer positivo: (a) Propagación del streamer vía fotoionización. (b) Distorsión del campo eléctrico en las proximidades del frente del streamer. Fig. 1.2 - 5 : Esquema representativo de la evolución de un str eamer negativo: (a) Propagación del streamer . (b) Distorsión de campo eléctrico en las proximidades del frente del streamer. 16 Capítulo 1 | Descargas Eléctricas descarga. 7. La interacción del plasma con la superficie supone todo un reto incluso en los casos más sencillos. En la inmensa mayoría de simulaciones a nivel de volumen, el tratamiento de interacción descarga/plasma/radiación-superficie se limita a un conjunto elemental de procesos basados en la neutralización y desexcitación (una discusión detallada de los principales mecanismos considerados en el contexto de las descargas puede verse en §2.6). Conforme a las dificultades que plantea una descripción microscópica consistente, los objetivos de esta Tesis pueden resumirse en los siguientes puntos: 1. Plantear un modelo microscópico consistente que describa el régimen transitorio de la descarga corona en oxígeno (incluyendo radicales, estados excitados, estados metaestables, estados vibracionales, electrones, iones, especies neutras, etc). Este tipo de modelado ofrece una descripción no solo microscópica sino también multiespecie que revela el papel de cada una de éstas en según qué etapa de la descarga (capítulos 2 y 4). 2. Definir un esquema numérico basado en elementos finitos (FEM) para simular la descarga. Debido a las inestabilidades asociadas a los problemas de convección dominante en la formulación estándar, se utilizaran esquemas de tipo Petrov-Galerkin (§3.3), donde el espacio de soluciones no coincide con el espacio de las funciones test (o funciones peso). Así mismo, se aprovecha la discretización temporal para estabilizar el esquema en regiones de reacción dominante (§3.4- §3.6). 3. Definir una formulación consistente para los procesos superficiales en el marco de una descripción de Langmuir-Hinshelwood bajo una aproximación de campo medio. De esta forma todas las especies adsorbidas se consideran distribuidas aleatoriamente en la superficie y se consigue prescindir de la interacción con los estados adsorbentes en la superficie (§2.6 y §4.2-§4.3). Bibliografía 17 Bibliografía 1. Langmuir, I. Oscillations in Ionized Gases. Proc. Natl. Acad. Sci. U. S. A. 14, 627–637 (1928). 2. Rogoff, G. L. Special issue on applications of partially ionized plasmas. IEEE Trans. Plasma Sci. 19, 989–990 (1991). 3. Raizer, Y. P. Gas Discharges Physics. (Springer-Verlag Berlin Heidelberg, 1991). 4. Fridman, A. Plasma Chemistry. (Cambridge University Press, 2008). 5. Fridman, A. & Kennedy, L. A. Plasma Physics and Engineering. (CRC Press, 2011). 6. Hagstrum, H. D. Theory of Auger Ejection Electrons of Electrons from Metal by Ions. Phys. Rev. 96, 336–365 (1954). 7. Lieberman, M. A. & Lichtenberg, A. J. Principles of Plasma Discharges and Materials Processing. (John Wiley & Sons, 2005). 8. Meek, J. M. & Craggs, J. D. Electrical Breakdown of Gases. (Wiley, 1978). 9. Raether, H. Electron Avalanches and Breakdown in Gases. (Butterworth, 1964). 10. Loeb, L. B. Basic Processes of Gaseous Electronics. (University of California Press, 1960). 11. Davies, A. J., Davies, C. S. & Evans, C. J. Computer simulation of rapidly developing gaseous discharges. Proc. Inst. Electr. Eng. 118, 816–823 (1971). 12. Druyvesteyn, M. & Penning, F. The Mechanism of Electrical Discharges in Gases of Low Pressure. Rev. Mod. Phys. 12, 87–174 (1940). 13. Roth, J. R. Industrial Plasma Engineering: Volume 1 - Principles. 1, (CRC Press, 1995). Bibliografía 19 un CAPÍTULO 2 21 Un Modelo Microscópico Contenido 2.1 Introducción ................................. 21 2.2 Una perspectiva microscópica ...................... 22 2.3 La aproximación de acoplamiento débil ............... 26 Colisiones entre partículas en las descargas—27 Factorización de 𝛼𝛽 (2)—30 2.4 La ecuación de Boltzmann ........................ 37 2.5 La aproximación fluida .......................... 39 2.6 Efectos de superficie ............................ 42 Emisión termoiónica—43 Fotoemisión (efecto fotoeléctrico)—45 Emisión por efecto campo—45 Emisión secundaria—46 Bibliografía 51 2.1 Introducción La descripción microscópica de la evolución temporal de las descargas puede reducirse al comportamiento de un conjunto de partículas cargadas en el seno de campos electromagnéticos. Esta Tesis se encuadra dentro de la teoría clásica de transporte1 (teoría cinética de gases), donde fundamentalmente se Capítulo “ …First the variations in velocity develop majestically, then the equations of state enter on one side, the equations of the central motion on the other; ever higher surges the chaos of formulae. Suddenly, four words sound out:‘Put n=5’. The evil demon v (velocity) disappears…” LUDWIG BOLTZMANN, Vorlesungen über Gastheorie, 1898. 22 Capítulo 2 | Un Modelo Microscópico interpreta el plasma como un conjunto de partículas cargadas que interaccionan vía colisiones binarias en presencia de campos electromagnéticos. Esto significa que se adopta la ecuación de Boltzmann (BE) para describir la dinámica microscópica de este tipo de sistemas. Aunque la naturaleza de las fuerzas electromagnéticas y el largo alcance de éstas contradicen la hipótesis de interacción de corto alcance implícita en la BE, desde un punto de vista de la teoría de perturbaciones tiene sentido utilizar la BE como aproximación de primer orden. Una amplia gama de los fenómenos asociados a los plasmas pueden explicarse bajo este tipo de aproximaciones. Por el contrario todos aquellos efectos asociados a modos colectivos y con ello relacionados con el carácter de largo alcance de la interacción de Coulomb (como por ejemplo el amortiguamiento de Landau o las inestabilidades en el espacio de velocidades entre otros) quedan fuera del ámbito de este trabajo y requieren una imagen más exacta del sistema (teoría anómala del transporte2,3). El objetivo de este capítulo es introducir una imagen microscópica de las descargas eléctricas en el ámbito de la teoría cinética para obtener finalmente las ecuaciones macroscópicas que gobiernan la denominada aproximación fluida. Entender el contexto de las aproximaciones que se realizan en la descripción fluida es crucial para resolver algunos problemas numéricos que pueden derivarse de las limitaciones del modelo físico (en §4.2 retomaremos una discusión detallada acerca de cómo interactúan el modelo físico-químico y el numérico, y la necesidad de integrar ambos enfoques para una descripción coherente de las descargas). Es por esto que se presta especial atención al desarrollo e hipótesis implícitas en la descripción fluida a partir de la BE. El planteamiento que se sigue para este propósito parte del sistema de ecuaciones BBGKI (generalización de la BE, véase §2.2) para llegar propiamente a la BE. 2.2 Una perspectiva microscópica La descripción de un sistema de partículas deja de ser una cuestión de la mecánica clásica en tanto , no solo por la dimensión del sistema de ecuaciones a considerar sino también por la imposibilidad práctica de fijar las condiciones iniciales para cada una de las partículas que definen el sistema. El estado dinámico de un sistema aislado de partículas con grados de libertad ( , en el caso más simple ), entre las que se distinguen 2.2 Una perspectiva microscópica 23 especiesa, queda definido completamente por medio de sus correspondientes posiciones y velocidades 𝑁 1 α 𝑁 𝜈 𝑁 1 𝛼 𝑁 𝜈 2.2 . 1 En el marco de la mecánica clásica 5,6 la evolución del sistema queda determinada a partir de la trayectoria seguida por éste en el espacio de las fases. Cada uno de los puntos de este espacio (ec. 2.2.1 ) determina las características energéticas del sistema, y con ello la trayectoria fásica queda ligada a la evolución energética del mismo. Bajo este enfoque microscópico, la evolución de cualquier propiedad del sistema requiere la solución de las ecuaciones de movimiento de todas las partículas que lo forman. Esta situación no sólo es inverosímil dada la inexactitud con la que pueden especificarse las condiciones iniciales, sino también por la imposibilidad práctica de resolver un sistema de ecuaciones de este tipob. La física estadística7,8 salva esta dificultad a través de la función de distribución, 𝑁 𝑁 𝑁 𝑁 2.2 . 2 a partir de la cual se mide la probabilidad de encontrar el sistema en un elemento infinitesimal de volumen del espacio fásico, en torno al punto 𝑁𝑁. En este contexto, cualquier propiedad microscópica, 𝑁 𝑁 , tiene su reflejo en una magnitud macroscópica, , a partir del promedio de sobre el espacio de fases (véase Fig. 2.2-1) Por tanto en una descripción estadística (clásica), el estado del sistema queda completamente caracterizado por medio de la función de distribución, . En concreto, cuando las partículas sigan una dinámica de Hamilton5,6, la evolución de la función de distribución del sistema viene dada por una ecuación de Liouville5,9 𝑖 𝐫 𝑖 𝑖 𝐯 𝑖 𝑁 𝑖 = 1 2.2.4 donde, en el caso más simple, los términos de aceleración para un conjunto a En principio se consideran idénticas todas las partículas de una misma especie. b Conviene tener presente que un sistema de partículas libres para un mol de sustancia resultaría en un sistema de 3 𝐴 ecuaciones de movimiento (con 𝐴 el número de Avogadro). 𝑁 𝑁 𝑁 𝑁 𝑁 𝑁 2.2.3 24 Capítulo 2 | Un Modelo Microscópico de partículas cargadas pueden escribirse según la expresiónc: 𝑖 𝑖 𝑖 ext ) 𝐫 ( 𝑖 , 𝑗 ) 𝑖𝑗 2.2.5 La primera componente a la derecha de 2.2.5 refleja la contribución del campo eléctrico externo, mientras que la segunda alude a la interacción interna, donde 𝑖𝑗 representa el potencial eléctrico. Así, en la relación 2.2.5 se asume la aproximación cuasielectrostática, donde no se consideran campo magnéticos externos aplicados y se desprecian los inducidos por las corrientes internas asociadas a la descarga, ya que el orden de estas corrientes (en el ámbito de la descarga corona, donde se encuadra este trabajo) varía en el rango de a . Desafortunadamente, establecer la relación que vincula y las variables naturales del problema es una tarea extremadamente compleja, cuando no imposible. No obstante, una descripción del sistema a ese nivel detalle es innecesaria en la gran mayoría de las situaciones prácticas. En tales casos, el interés suele estar dirigido hacia un grupo específico de partículas, por lo que suele prescindirse de toda la información redundante acerca de las restantes partículas integrando sobre las coordenadas y velocidades de estas últimas. Se define así la función de distribución reducida, (𝑠) , de acuerdo con la expresiónd,7 c Se utiliza por sencillez un sistema Gaussiano de unidades. d Es posible encontrar en la literatura versiones alternativas a la definición de las funciones Fig. 2.2 - 1 : Representación esquemática que refleja la diferencia entre la descripción de la mecánica clásica y la mecánica estadística clásica para un sistema de muchas partículas: en la primera (izquierda) la función toma un valor definido. En cambio, en el segundo caso (derecha), para fijado, se considera como una función aleatoria. 2.2 Una perspectiva microscópica 25 ( 𝑠 ) 1 1 𝑠 𝑠 𝑠 𝑠 + 1 𝑠 + 1 𝑁 𝑁 2.2.6 Haciendo uso de 2.2.6 , la relación 2.2.4 conduce a la expresión: ( 𝑠 ) 𝑖 𝐫𝑖 𝑠 𝑖=1 ( 𝑠 ) 𝑠 𝑖 𝐫𝑖 𝑠+1 𝑁 𝑁 𝑖=𝑠+1 𝑠 𝑁 𝑖 = 1 𝑖 𝐯 𝑖 𝑠 + 1 𝑁 2.2.7 donde 𝑖 𝑖 𝑖. La ecuación 2.2.7 puede transformarse en una cadena de ecuaciones entrelazadas conocida como la jerarquía de ecuaciones BBGKIe,10–13 ( 𝑠 ) 𝑖 𝐫𝑖 𝑠 𝑖=1 ( 𝑠 ) 𝑖 𝐯𝑖 ( 𝑠 ) 𝑠 𝑖=1 𝑠 𝑖 = 1 𝑖 𝑖 𝐫 𝑖 𝑖 , 𝑠 + 1 𝐯 𝑖 (𝑠+1) 𝑠 + 1 𝑠 + 1 2.2.8 donde 𝑖𝑖 𝑖 ext) 𝐫𝑖 𝑠 𝑗 = 1 ( 𝑖 ≠ 𝑗 ) 𝑖𝑗 2.2.9 El conjunto no cerrado de ecuaciones 2.2.8 contiene toda la información estadística del sistema. Sin embargo, en la práctica, y en particular en la teoría cinética y teoría del transporte, la mayor parte de esta información es irrelevante. El carácter no cerrado de las ecuaciones 2.2.8 muestra la dificultad esencial de la física estadística: para determinar (2) se requiere conocer (3), y así sucesivamente. Por tanto, una solución exacta del problema requeriría la solución de las (prácticamente infinitas) ecuaciones acopladas para las correspondientes funciones de distribución reducidas. Existe una extensa bibliografía14,8,15,f donde se exponen y se discuten en detalle las distintas aproximaciones que se consideran en el marco de la física de distribución en términos de la condición de normalización. En este trabajo se adopta el criterio indicado en la referencia donde se normaliza al volumen. e Yvon, Bogolyubov, Born y Green, y Kirkwood obtuvieron el sistema de ecuaciones 2.2.8 de manera independiente, de ahí el acrónimo con el que se conocen estas ecuaciones en la literatura. f Se citan aquí solo algunas de las fuentes clásicas. 32 Capítulo 2 | Un Modelo Microscópico la distancia a partir de la cual una partícula deja de percibir la interacción individual (partícula-partícula), para interaccionar vía campo medio creado por todas las restantes partículas. Es por esto, que desempeña un papel esencial a hora de contemplar interacciones débiles y, con ello, el rango de validez de la BE. Para realizar una estimación representativa de este parámetro es posible considerar aquellas regiones dominadas por la ionización, donde las especies predominantes serán los electrones e iones positivos producidos vía ionización directa. En tal caso, en el intervalo en el que las densidades locales sean aproximadamente iguales, 𝑒 𝑖, la longitud de Debye resulta 𝐷 𝑒 2 𝐵 − 1 / 2 𝑒 − 3 1 / 2 2.3.16 Típicamente en descargas eléctricas, 𝑒4 y 𝑒12 −3, por lo que 𝐷. De acuerdo con 2.3.13 𝛼𝛽 𝛼 𝛽 𝛼𝛽 𝛼 𝛽 𝐵 12 − 𝓇 12 / 𝜆 𝐷 − 𝓻 12 / 𝜆 𝐷 12 𝐷 2.3.17 donde 𝑒 𝐷 3 2.3.18 La magnitud se conoce como parámetro del plasma (adimensional), y expresa la inversa del número de partículas (electrones) contenidos en una celda cúbica de lado igual 𝐷. Atendiendo al significado físico de 𝐷, cuando 12 𝐷, de 2.3.17 resulta 𝛼𝛽 2.3.19 Este hecho manifiesta que en el caso , cualquiera que sea el valor de 12, la función de distribución por pareja podrá factorizarse por medio de las correspondientes funciones de una sola partícula. Típicamente en descargas eléctricas, el valor suele estar comprendido en el rango −5 −3 . Este hecho constata la validez de esta aproximación en condiciones de equilibrio termodinámico local (LTE), esto es, 𝛼𝛽 𝛼 𝛽 2.3 . 20 Validez fuera del equilibrio termodinámico. Método de Dupree18 Si bien el resultado 2.3.19 pone de manifiesto la bondad de la aproximación 2.3 La aproximación de acoplamiento débil 33 de acoplamiento débil en el equilibrio termodinámico, solo podemos esperar que esta simplificación siga siendo válida en situaciones de fuera del equilibrio. No obstante, tanto para sistemas que interaccionan débilmente como para sistemas con una interacción arbitraria pero con baja densidad, puede establecerse que una aproximación del tipo 2.3.20 es válida (véase Ajiezer19 para una discusión detallada). En esencia, esta aproximación puede interpretarse en términos de las diferentes escalas de tiempo que siguen las funciones de distribución de una y muchas partículas. Este es el argumento que sigue Bogolyubov11 (conjetura de Bogolyubov) para explicar la evolución en el tiempo que sigue el sistema cuando parte de una situación fuera del equilibrio. Para ello define esta evolución de acuerdo con tres escalas características de tiempo8,14: Fase 1: Inicial. Relacionada con el periodo de tiempo, 𝑐 𝑐 , siendo 𝑐 el radio de correlación y la velocidad media de las partículas. Es claro que esta fase refleja el efecto de las correlaciones entre las partículas. Estas funciones sólo son significativas en cierta región del espacio, y con ello durante un intervalo de tiempo limitado. Por tanto, puede entenderse esta fase como una versión extendida del concepto de colisión dado en la teoría cinética de gases. Fase 2: Cinética. Asociada al tiempo medio entre colisiones, 𝑚𝑓𝑝 𝑚𝑓𝑝 , donde 𝑚𝑓𝑝 representa la distancia media entre colisiones sucesivas y la velocidad media de las partículas. Fase 3: Hidrodinámica. Comprende la escala de tiempo de la relajación macroscópica. La longitud característica de esta escala, 𝐻, mide la razón de la variación espacial de (1) en sistemas no homogéneos. A menudo suele estar relacionada con el gradiente local de alguna de las magnitudes macroscópicas tales como densidad, temperatura, etc, 𝐻 − 1 𝐻 𝐻 2.3.21 donde representa la velocidad media de las partículas. Este período se caracteriza por el hecho de que después de algunas colisiones, el sistema alcanza un equilibrio local. A partir de estos estados de cuasiequilibrio, el sistema puede describirse por medio de la función de distribución de una partícula, e incluso haciendo uso de los primeros momentos. La conjetura de Bogolyubov (o régimen hidrodinámico, como se conoce en un contexto más amplio que el de la teoría cinética de gases) responde a la intuición física por la que se entiende que en el proceso de relajación, las fases anteriormente citadas preservan dicho orden con un carácter perfectamente 34 Capítulo 2 | Un Modelo Microscópico segregado, es decir, 𝑐 𝑚𝑓𝑝 𝐻 2.3 . 22 o bien, que, 𝐷 𝑐 𝑚𝑓𝑝 𝐻 2.3 . 23 La dinámica secuencial a la que responde la intuición de Bogolyubov permite interpretar una evolución segregada de las funciones de distribución 𝛼 𝛼𝛽 𝛼𝛽𝛾 En estas condiciones tiene sentido una aproximación de primer orden (en el contexto de la aprox. de acoplamiento débil) para el problema de las colisiones binarias, donde 𝛼𝛽 𝛼 𝛽 Antes de aclarar la relación entre estas escalas de tiempo y las funciones de distribución conviene tener presente las limitaciones asociadas a este orden de escalas (ec. 2.3.22 y/o 2.3.23 ). Particularmente en el caso de las descargas en gases pueden darse situaciones donde el régimen hidrodinámico no se satisface, tales como los frentes de choques que pueden desarrollarse en la evolución de la descarga donde el gradiente de densidad o de alguna otra variable macroscópica tiene lugar en una región del espacio extremadamente pequeña (como ejemplo de este tipo de situaciones puede verse §3.6 y/o §4.5). Salvando este tipo de situaciones, la relación entre las escalas de tiempo definidas por Bogolyubov y las funciones de distribución puede interpretarse de acuerdo con el hecho de que en la práctica todas las ecuaciones cinéticas surgen de algún desarrollo perturbativo. En este sentido la aproximación de acoplamiento débil sugiere como estado de referencia aquél en el que todas las interacciones quedan “apagadas”. Es decir, se asume que el potencial de interacción entre las partículas es despreciable a cualquier distancia. En estas circunstancias, de acuerdo con la interpretación del parámetro del plasma (ec. 2.3.18 ), se tiene 𝑖𝑗 𝑖𝑗 2.3.24 A partir de 2.3.24 , se define un sistema débilmente acoplado en las siguientes condiciones 0 2.3.25 donde 0 representa una distancia finita arbitraria. Así, al imponer la condición 2.3.25 sobre los dos primeros órdenes de 2.3 La aproximación de acoplamiento débil 35 la BBGKI (ec. 2.2.8 ) se obtiene una evolución segregada para las funciones de distribución, tal y como prescribe el régimen hidrodinámico. En efecto, tomando los dos primeros órdenes de 2.2.8 , para α se tiene: 𝛼 𝛼 𝐫𝛼 𝛼 𝛼 𝛼 ext ) 𝐯𝛼 𝛼 𝛽 𝛼 𝛼 𝐫 𝛼 𝛼𝛽 𝐯 𝛼 𝛼𝛽 𝛽 𝛽 𝜈 𝛽 = 1 2.3.26 donde se ha utilizado que 𝛽 y 𝛽 𝛽 . Y, para el caso de 𝛼𝛽: 𝛼𝛽 𝛼 𝐫𝛼 𝛼𝛽 𝛽 𝐫𝛽 𝛼𝛽 𝛼 𝛼 ext) 𝐫𝛼𝛼𝛽 𝐯𝛼 𝛼𝛽 𝛽 𝛽 ext) 𝐫𝛽𝛽𝛼 𝐯𝛽 𝛼𝛽 𝛽𝛼 𝛼 𝜈 𝛽=1 𝐫𝛼𝛼𝛽 𝐯𝛼 𝛼𝛽𝛾 𝛽𝛽 𝛼 𝛽 𝛽 𝜈 𝛼 = 1 𝐫 𝛽 𝛽𝛼 𝐯 𝛽 𝛼𝛽𝛾 𝛼 𝛼 2.3.27 que, nuevamente, hace uso de la condición 𝛼 y 𝛽. Luego, la aproximación de acoplamiento débil (ec. 2.3.25 ) sobre 2.3.26 y 2.3.27 , supone que el potencial de interacción entre partículas aparece como una corrección de primer orden. Esto significa que en una aproximación de orden cero, respecto del parámetro , no se considera el efecto de interacción entre partículas, 𝛼 𝛼𝛽 2.3.28 Como resultado de ello puede despreciarse el término integral en 2.3.26 , resultando una derivada material de 𝛼 igualada a cero, es decir, 𝛼 permanece constante. Por tanto, durante la fase inicial, 𝑐, 𝛼 no varía. La siguiente fase contenida en el régimen hidrodinámico (escala cinética) contempla la variación de 𝛼 debido al efecto de la interacción entre partículas. De modo que en una corrección de primer orden sobre la fase previa se desprecian los dos últimos términos en 2.3.27 , ya que 𝛼𝛽𝛾 2 2.3.29 Esto es así, dado que las correlaciones de tres cuerpos requieren de, al menos, 36 Capítulo 2 | Un Modelo Microscópico dos procesos de interacción (sea por ejemplo: y , donde los índices se refieren a las etiquetas de las partículas implicadas, , y , respectivamente). Lo que significa que 𝛼𝛽 permanece constante durante esta fase. Por el contrario, la presencia de 𝛼𝛽 en 2.3.26 hace que en esta escala de tiempo, 𝛼 varíe. De este modo, en la escala cinética y en condiciones de régimen hidrodinámico, 𝛼𝛽 (y en general los órdenes superiores) sufre un borrado de memoria, ya que su dependencia en el tiempo se refleja exclusivamente a partir de la variación de 𝛼, ( 𝑠 ) 1 𝑠 𝑡 ≫ 𝜏 𝑐 ( 𝑠 ) 1 𝑠 ( 1 ) ′ 2.3.30 donde 𝑖 𝑖 𝑖 . Por tanto, fuera del equilibrio, la aproximación de acoplamiento débil induce la secuencia de escalas prescritas por el régimen hidrodinámico. En estas condiciones, las funciones de distribución admiten un desarrollo perturbativo basado en un desarrollo en potencias de 𝛼 𝛼 , 0 𝛼𝛽 𝛼 𝛽𝛼𝛽 𝛼𝛽,0 𝛼𝛽,1 𝛼𝛽𝛾 𝛼𝛽𝛾 , 0 𝛼𝛽𝛾 , 1 2 𝛼𝛽𝛾 , 2 2.3.31 De este modo, en una aproximación de primer orden para el problema de las colisiones binarias, el sistema 2.2.8 puede reducirse a su primer orden, ya que de acuerdo con 2.3.31 𝛼𝛽 𝛼 𝛽 Finalmente, tanto en el contexto de LTE como en el de régimen hidrodinámico con predominio de colisiones binarias es posible describir un sistema de muchas partículas (ec. 2.2.8 ) a partir de la función de distribución de una sola partícula, 𝛼 Condiciones éstas que delimitan la ámbito de la BE, ya que constituyen las bases de esta afamada ecuación. Limitaciones de una descripción asintótica El carácter lineal de la ecuación de Liouville 2.2.4 facilita el uso de desarrollos perturbativos para resolver el sistema de ecuaciones 2.2.8 . Sin embargo, en el estudio de las descargas eléctricas y, en general, de los plasmas, la evolución secular j de la soluciones puede conducir a que la representación dada por un desarrollo perturbativo sea incorrecta. Un ejemplo sencillo pueden ser aquellas propiedades macroscópicas que relajan con el tiempo de acuerdo con la expresión j Funciones que crecen sistemáticamente en el tiempo, tales como potencias del tiempo, t. 2.4 La ecuación de Boltzmann 37 − 𝜖𝑡 2 2 3 3 donde representa el parámetro de orden del desarrollo. En este caso, aun cuando tome valores pequeños que justifiquen un desarrollo limitado a los primeros órdenes, la descripción asintótica para tiempos suficientemente largos dejará de ser válida. Este hecho pone de manifiesto la importancia práctica que supone fijar ciertos límites a una descripción basada en este tipo de desarrollos, donde se trunca la descripción del sistema a unos pocos términos. En esencia, esta característica muestra la incapacidad de un desarrollo en pocos términos para describir la relajación a tiempos largos de un sistema macroscópico. 2.4 La ecuación de Boltzmann Como hemos visto la aproximación de acoplamiento débil contiene dos de las hipótesis cruciales en la BE: la condición de gas diluido, que se traduce en el predominio de las colisiones binarias; y la hipótesis de caos molecular que permite reducir el problema de colisiones binarias en la BBGKI al primer orden de ésta jerarquía de ecuaciones, factorizando la función de distribución por pareja, esto es, 𝛼𝛽 𝛼 𝛽 En estas condiciones, la expresión 2.3.26 puede reescribirse en la forma: 𝛼 𝛼 𝛼 tot ) 𝐯 𝛼 𝛼 𝛼 𝑐 2.4.1 donde se han utilizado las relaciones, 𝛼 𝐫𝛼 tot) ext) 𝛽 𝐫𝛼 𝛽 𝛼 𝛽 𝛽 𝛽 𝛽 𝜈 𝛽=1 𝛼 𝑐 𝛽 𝛼𝛽 int) 𝛼𝛽 int) 𝐯 𝛼 𝛼 𝛽 𝛽 𝛽 𝜈 𝛽 = 1 2.4.2 siendo 𝛼𝛽 int) la componente de la aceleración debida a la interacción entre partículas vía campo medio; el factor 𝛼𝛽 int) 𝛼𝛽 int) en el término de colisión, la aceleración debida a la interacción partícula-partícula; y el campo medio generado por todas las partículas. Este tipo de campos también 38 Capítulo 2 | Un Modelo Microscópico se conocen como campos autoconcordantes o autoconsistentesk. Además, estos campos internos satisfacen las ecuaciones de Maxwell (aprox. cuasielectrostática, véase §2.2) 𝛽 𝛽 𝛽 𝜈 𝛽=1 2.4.3 Por tanto la versión 2.4.1 para el primer orden de la BBGKI distingue dos tipos de términos de colisión (en primer orden de aproximación): por un lado un término de colisión asociado a la interacción vía campo medio y, por otra parte, un término de colisión relacionado con la interacción partículapartícula (colisión de núcleo duro, de acuerdo los resultados mostrados en §2.3). De manera que este último responde al prototipo de término de colisión de la BE. En la práctica, la escala temporal observable de la descarga cumple que 𝑐, y con ello el término de colisión 𝛼 𝑐 en 2.4.1 puede transformarse, finalmente, en un término de colisión de tipo Boltzmann, esto es, 𝛼 𝑐 𝛼𝛽 𝛼𝛽 𝛼𝛽 𝛼 ′ 𝛽 ′ 𝛼 𝛽 2.4.4 donde los superíndices “ ” aluden a las funciones de distribución después de la colisión, mientras que las funciones 𝛼 y 𝛽 están referidas a la situación previa a la colisión. La transición desde el término de colisión de 2.4.2 a 2.4.4 , siguiendo los argumentos de Kirkwood13,20, se deduce mediante un promedio temporal sobre la escala 𝑐. Esto significa que 2.4.1 se promedia en un intervalo de tiempo ∗, suficientemente largo en comparación con 𝑐, pero a su vez suficientemente corto con relación a 𝑚𝑓𝑝 (escala de tiempo para la variación de 𝛼) 𝛼 ∗ 𝛼 𝜏 ∗ 0 2.4.5 Como resultado de ello 2.4.1 queda como sigue 𝛼 𝛼 𝛼 tot ) 𝐯 𝛼 𝛼 𝛼 𝑐 2.4.6 donde k La notacion de autoconsistente para este tipo de campos se debe a que sus fuentes dependen de la función de distribución, 𝛽, y ésta a su vez depende del campo a través de la ecuación cinética. Por tanto, ambos deben calcularse simultáneamente de una manera consistente. 2.5 La aproximación fluida 39 𝛼 𝑐 𝛽 ∗ 𝛼𝛽 int ) 𝛼𝛽 int ) 𝐯 𝛼 𝛼 𝛽 𝛽 𝛽 𝜏 ∗ 0 𝜈 𝛽 = 1 2.4.7 La relación 2.4.7 puede transformarse en un término de colisión del tipo (Kirkwood13,20): 𝛼 𝑐 𝛼 ′ 𝛽 ′ 𝛼 𝛽 𝛽 2.4.8 en las condiciones anteriormente citadas: 1. Dominio de colisiones binarias. 2. Hipótesis de caos molecular: 𝛼𝛽 𝛼 𝛽 Finalmente, la forma clásica del término de colisión de Boltzmann (ec. 2.4.4 ) puede deducirse al integrar sobre la distancia relativa ( 𝛼𝛽 𝛼 𝛽 ) en 2.4.8 , adoptando para ello el sistema en coordenadas cilíndricas . La dirección se toma antiparalela al vector de velocidad relativa, 𝛼𝛽 𝛼 𝛽 (véase Fig. 2.4-1). 2.5 La aproximación fluida El modelo macroscópico o aproximación fluida en las descargas está basado en el método de los momentos de Grad21. Esta aproximación caracteriza el estado del fluido por medio de 13 magnitudes escalares: densidad, 𝛼 ; velocidad media, 𝛼 ; tensor de presión, 𝛼 ; y el vector flujo de calor, 𝛼 . La conexión entre el modelo microscópico derivado de la BE y las Fig. 2.4 - 1 : Esquema ilustrativo de la geometría de una colisión binaria. 40 Capítulo 2 | Un Modelo Microscópico magnitudes termodinámicas viene dada por los momentos de la función de distribución. De este modo, 𝛼 𝛼 2.5.1 define la densidad de la especie como una función que depende solo de la posición para un instante dado, . Si la integración en 2.5.1 se realizar sobre la ecuación que gobierna la dinámica de 𝛼 (ec. 2.4.6 ) en lugar de 𝛼, entonces se obtiene una expresión para la dinámica de la densidad en conexión con la dinámica microscópica, esto es, 𝛼 𝛼 𝛼 𝛼𝛽 𝛽 2.5.2 Esta relación responde al prototipo de ecuación de continuidad, donde 𝛼 𝛼 𝛼 𝛼 𝛼 2.5.3 representa la velocidad media de dicha especie y 𝛼𝛽 el efecto de creación/destrucción de partículas de la especie como resultado de la interacción con las restantes especies 𝛼𝛽 𝛼 𝑐 − 𝛽 2.5.4 Para resolver la ecuación de continuidad es necesario conocer la ecuación para la velocidad media. Para ello, se sigue el mismo proceso que para la densidad, en relación ahora con 2.5.3 , de lo que se obtiene la ecuación de conservación del momento 𝛼 𝛼 𝛼 𝛼 𝛼 tot ) 𝛼 𝛼 2.5.5 donde, 𝛼 𝛼 𝛼 𝛼 𝛼 2.5.6 representa el tensor de presiones y 𝛼 la fricción neta debido al efecto de las colisiones 𝛼 𝛼 𝛼𝛽 𝛽 𝛼 𝛼 𝑐 − 𝛽 𝛽 2.5.7 En particular, cuando la distorsión de la función de distribución respecto de la situación de equilibrio termodinámico es relativamente débil, 2.4.7 2.5 La aproximación fluida 41 puede reescribirse en términos del operador de colisión de Krook22, esto es, 𝛼 𝑐 𝛼 𝛼 , 0 2.5.8 donde 𝛼,0 representa la función de distribución en el equilibrio termodinámico. Resulta, entonces, 𝛼 𝛼 𝛼 𝛼𝛽 ∗ 𝛼 𝛽 𝛽 2.5.9 donde 𝛼𝛽 ∗ 𝛽 𝛽 𝛾 𝛾 𝛾 𝑚 2.5.10 representa la frecuencia de colisiones entre las especies y . La relación 2.5.10 está referida a las colisiones elásticas cuya sección eficaz ( 𝑚) suele ser dominante en la mayoría de las situaciones prácticas para las descargas. En cualquier caso, la ecuación 2.5.10 es fácilmente generalizable a cualquier tipo de proceso, tales como la ionización directa, captura electrónica, recombinación, excitación, ionización por pasos, etc. En general, la frecuencia de colisión asociada a un proceso de interacción entre electrones () y las moléculas del medio ( ) puede expresarse como sigue 𝑒𝑛 𝑒𝑛 − 1 𝑛 𝑒𝑛 𝑒 𝑣 𝑛 𝑒 𝑒 𝑒𝑛 3 + ∞ 0 2.5.11 donde el factor 𝑒𝑛 𝑒 𝑣 representa la tasa de reacción, siendo 𝑒𝑛 la sección eficaz de colisión electrón-molécula neutra del gas para el proceso en cuestión. Nuevamente, para resolver la ecuación del momento 2.5.5 es necesario conocer la dinámica del tensor de presiones. Para ello se sigue el mismo proceso que el utilizado para la densidad y conservación del momento, de lo que resulta: 𝛼 𝛼 𝛼 𝛼 𝛼𝛽 𝛼 𝐸 , 𝛼 2.5.12 donde 𝛼 representa el vector flujo de energía térmica o flujo de calor 𝛼 𝛼 𝛼 2 𝛼 𝛼 2.5.13 Por otra parte, el último término de 2.5.12 representa el flujo de calor asociado al efecto de las colisiones, 48 Capítulo 2 | Un Modelo Microscópico emisión secundaria de un electrón. En el primer caso (neutralización, véase Fig. 2.6-4–(a)) el electrón es capturado en uno de los estados excitados del ion positivo + 𝑀 ′ ∗ 2.6 . 12 Este mecanismo conduce a un proceso de relajación en el que se emite un fotón de energía igual a X ∗ X 2.6 . 13 donde X ∗ X + 𝑀′ . El segundo caso (emisión secundaria, véase Fig. 2.6-4–(b)) tiene lugar cuando el electrón sustraído es capturado por el estado fundamental del ion, + 𝑀 ′ 2.6 . 14 En estas circunstancias el defecto de energía 𝑀 ′ X + X 2.6 . 15 puede ser absorbido por un electrón secundario en el cátodo, y en caso de superar el trabajo de salida podrá abandonar la superficie metálica. Por este motivo la emisión secundaria vía mecanismo de Auger suele representarse Fig. 2.6 - 4 : Diagrama esquemático del proceso de neutralización de Auger que tiene lugar cuando un ion “lento” se aproxima a la superficie metálica: (a) neutralización; y (b) emisión secundaria. 2.6 Efectos de superficie 49 por la secuencia + 𝑀 ′ 𝑀 ′′ 𝑘 2.6 . 16 donde ′0, ′′ 0 y 𝑘 0 𝑘 (véase Fig. 2.6-4–(b)). De acuerdo con la Fig. 2.6-4–(b), la energía máxima con la que el electrón puede escapar del metal viene dada por la expresión 𝑘 máx. X + ′ 0 2.6 . 17 donde se hace uso de las relaciones derivadas de la conservación de la energía en 2.6.16 , esto es, X + − 𝑀 X + X − 𝑀 𝑘 2.6.18 de las que se deduce que: ′ ′′ 𝑘 0 X + ′ X + ′ X + X + − 𝑀 X − 𝑀 2.6.19 siendo X+−𝑀 y X−𝑀 las energía de interacción del ion y el átomo con la superficie metálica a una distancia , respectivamente. En la práctica la condición 2.6.17 suele utilizarse para definir la relación experimental del coeficiente 25 (ec. 1.2.1 ). Emisión secundaria estimulada por impacto de estados excitados Los procesos de desexcitación vía mecanismos de Auger pueden describirse de un modo completamente análogo al proceso descrito en 2.6.12 (véase la referencia: Hagstrum24). En estos casos, la secuencia puede reducirse a la relación ∗ 𝑀 ′′ 𝑘 2.6 . 20 Este mecanismo, a diferencia del proceso de neutralización previo, conduce a la emisión secundaria de un electrón del metal. Un balance energético similar a 2.6.18 conduce a las relaciones: X ∗ − 𝑀 X ∗ X − 𝑀 𝑘 2.6.21 de las que se deduce, ′′ 𝑘 0 X + ′ X + ′ X ∗ X ∗ − 𝑀 X − 𝑀 2.6.22 Luego, en este caso la energía cinética máxima con la que el electrón puede escapar del metal viene dada por la expresión: 50 Capítulo 2 | Un Modelo Microscópico 𝑘 máx. X ∗ ′ 0 2.6 . 23 Comparando 2.6.23 con 2.6.17 puede verse que el mecanismo de emisión secundaria por medio de la neutralización de un ion domina sobre la emisión secundaria estimulada por impacto de estados excitados. Una estimación de esta valoración puede hacerse despreciando la energía de interacción y considerando exclusivamente las energías: X+ y X∗. Teniendo en cuenta que típicamente 0, es claro que por medio de la emisión secundaria estimulada por impacto de estados excitados solo una pequeña proporción de los electrones en el metal consiguen abandonarlo. Por este motivo se considerará (véase §4.2) que el principal mecanismo para la emisión secundaria está determinado por la neutralización de los iones positivos en las proximidades del cátodo. Bibliografía 51 Bibliografía 1. Balescu, R. Transport processes in plasmas. 1 Classical transport theory. (Elsevier Science Publishers B.V., 1988). 2. Klages, R., Radons, G. & Sokolov, I. M. Anomalous transport: foundations and applications. (John Wiley & Sons, 2008). 3. Balescu, R. Aspects of anomalous transport in plasmas. (CRC Press, 2010). 4. Lebowitz, J. L. & Montroll, E. W. Nonequilibrium Phenomena: The Boltzmann Equation. (Elsevier Science Ltd, 1983). 5. Goldstein, H., Poole, C. & Safko, J. Classical Mechanics. (Addison Wesley, 2001). 6. Gantmájer, F. R. Mecánica analítica. (URRS, 2003). 7. Krall, N. A. & Trivelpiece, A. W. Principles of plasma physics. (McGraw-Hill, 1973). 8. Balescu, R. Equilibrium and nonequilibrium statistical mechanics. NASA STIRecon Tech. Rep. A 76, (1975). 9. Liboff, R. L. Kinetic theory. (Springer-Verlag, 2003). 10. Yvon, J. La théorie statistique des fluides et l’équation d’état. 203, (Hermann & cie, 1935). 11. Bogoliubov, N. N. Problems of dynamical theory in statistical physics. J PhysUSSR 10, 256 (1946). 12. Born, M. & Green, H. S. A general kinetic theory of liquids. I. The molecular distribution functions. Proc. R. Soc. Lond. Ser. Math. Phys. Sci. 188, 10–18 (1946). 13. Kirkwood, J. G. The statistical mechanical theory of transport processes I. General theory. J. Chem. Phys. 14, 180–201 (1946). 14. Montgomery, D. C. & Tidman, D. A. Plasma Kinetic Theory. (McGraw-Hill, 1964). 15. Ichimaru, S. Statistical Plasma Physics, Volume I: Basic Principles. (Westview Press, 2008). 16. Mayer, J. E. & Mayer, M. G. Statistical Mechanics of Fluids. (Wiley, 1940). 17. Landau, L. D. & Lifshitz, E. M. Statistical physics. (Butterworth Heinemann, Oxford, 1991). 18. Dupree, T. H. Dynamics of Ionized Gases. Phys. Fluids 4, 696 (1961). 19. Ajiezer, A. I. & Peletminski, S. V. Métodos de la Física Estadística. (MIR, 1981). 20. Kirkwood, J. G. The statistical mechanical theory of transport processes II. Transport in gases. J. Chem. Phys. 15, 72–76 (1947). 21. Grad, H. Note on N-dimensional hermite polynomials. Commun. Pure Appl. Math. 2, 325–330 (1949). 22. Bhatnagar, P. L., Gross, E. P. & Krook, M. A Model for Collision Processes in 52 Capítulo 2 | Un Modelo Microscópico Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems. Phys. Rev. 94, 511–525 (1954). 23. Fowler, R. H. & Nordheim, L. Electron Emission in Intense Electric Fields. Proc. R. Soc. Lond. Math. Phys. Eng. Sci. 119, 173–181 (1928). 24. Hagstrum, H. D. Theory of Auger Ejection Electrons of Electrons from Metal by Ions. Phys. Rev. 96, 336–365 (1954). 25. Raizer, Y. P. Gas Discharges Physics. (Springer-Verlag Berlin Heidelberg, 1991). Bibliografía 53 el CAPÍTULO 3 55 El método: FEM Contenido 3.1 Introducción ................................. 55 3.2 El método de los residuos ponderados ................ 57 Equivalencia entre las ecuaciones diferenciales y la formulación débil—57 Formulación Bubnov-Galerkin (método de Galerkin)—59 Formulación matricial—59 3.3 El problema de convección-difusión .................. 62 Representación estándar—62 Representación autoadjunta—65 Stream Upwind Petrov-Galerkin (SUPG)—66 3.4 Discretización temporal .......................... 69 El algoritmo—70 3.5 Precisión y Estabilidad .......................... 72 3.6 Un test numérico .............................. 76 3.7 Algunas conclusiones finales ....................... 83 Bibliografía 85 3.1 Introducción La estrategia natural cuando se abordan sistemas complejos (aquellos que no admiten una simplificación extrema como para poder reducir su comportamiento a pocos grados de libertad) consiste en interpretar el comportamiento de éstos a través de sus partes elementales (elementos). La Capítulo “ Experience is what you get when you don't get what you want.” DAN STANFORD 56 Capítulo 3 | El método: FEM esperanza de este planteamiento radica en dos aspectos fundamentales: el hecho de que el grado de complejidad del comportamiento de estos elementos sea ínfimo en comparación con el del problema global, y la posibilidad de reconstruir la dinámica del sistema como un todo a partir del ensamblaje de las partes. Históricamente, el método de los elementos finitos nace con los trabajos de Lord Rayleigh1 (1870) y W. Ritz2 (1909) en el contexto de métodos basados en principios variacionales. Posteriormente, el trabajo de B. Galerkin3 (1915) supuso una generalización del método para operadores que admitían una formulación funcional o bien que permitían una formulación débil. En 1943 Courant4 hizo una contribución significativa al método de Ritz usando funciones continuas definidas por secciones. Algunos años más tarde (1956), Clough5 desarrolló un método idéntico al usado por Courant para estudiar la rigidez y flexión de estructuras complejas. Fue entonces cuando Clough6 introdujo el concepto de elemento finito. En la actualidad existen multitud de herramientas de cálculo basadas en este método. En este trabajo se ha utilizado el software Comsol Multhiphysics para implementar las ecuaciones de la aproximación fluida en el ámbito de las descargas eléctricas. En particular, se ha hecho uso del módulo de ecuaciones matemáticas en formulación débil, con la intención de controlar todos los aspectos de estabilidad del esquema numérico. Desde luego, una discusión detallada del método de elementos finitos (FEM) va más allá de los objetivos de este trabajo, aunque merece la pena explorar las dificultades numéricas asociadas a la representación estándar en problemas de transporte. En este sentido, toda la atención se centra en el diseño de las funciones de forma y técnicas de discretización temporal. De hecho, en estos dos puntos se concentra la mayor parte del factor humano en el desarrollo del método. A lo largo de este capítulo se presenta la formulación streamline upwind Petrov–Galerkin7 (SUPG) como una alternativa para evitar las inestabilidades introducidas por la formulación estándar en problemas de convección dominante. No obstante, en problemas transitorios, la difusión artificial introducida por esta técnica puede llegar a ser excesiva sin un refinamiento adecuado de la malla. Para evitar este efecto de suavizado continuo de la solución se introduce, en el ámbito de la discretización temporal, una representación que combina los efectos de una descripción implícita y explícita ( –method 8–10). El resultado de esta derivada temporal promediada permite construir esquemas de alta resolución basados en el concepto de limitadores, originalmente desarrollado por Van Leer11,12 y Boris-Book13,14. En este escenario, se aprovechan las características de un esquema (Crank-Nicholson –CN–), ya que no producen amortiguación alguna en la 3.2 El método de los residuos ponderados 57 Formulación fuerte solución. Sin embargo, la imposibilidad de amortiguar los posibles errores numéricos conduce a inestabilidades, a menos que la evolución sea suficientemente suavea. Por este motivo, se construye un esquema híbrido, CN+SUPG, que permite filtrar la región de altas frecuencias asociadas a fuertes gradientes o discontinuidades. En la parte final del capítulo, para testear la bondad del esquema numérico se analiza el problema de la descarga corona (R. Morrow15) y se comparan los resultados con los originales, obtenidos por medio de un esquema de alta resolución (flux–corrected transport13 –FCT–) en diferencias finitas (FD). Se introduce entonces el efecto de los términos de reacción en la estabilidad del esquema, y se analiza el papel de la discretización temporal sobre la base de un análisis espectral del error16,17 (J. Von Neumann). 3.2 El método de los residuos ponderados Equivalencia entre las ecuaciones diferenciales y la formulación débil La ecuación diferencial que caracteriza un problema de convección-difusión- reacción (CDR) viene dada por la expresión: 0 𝑝 D 𝑠 𝑝 N 3.2.1 donde 𝑠 representa el vector normal a la superficie, apuntando desde el volumen (dominio) del problema hacia la superficie; mientras que define el interior del dominio en cuestión; y D N la frontera del mismo. Sea la solución de 3.2.1 , y una función arbitraria (suficientemente regular) tal que: D entonces, un planteamiento en residuos ponderado de 3.2.1 y su posterior integración por partes conducen a la expresión: 𝑡 ΩΩ 𝑝 Γ N Ω 3.2.2 a En la inmensa mayoría de las ocasiones se requerirían pasos de tiempo inasumibles en relación a la escala de tiempo de problema. 64 Capítulo 3 | El método: FEM define el número de Reynolds de la celda (en ocasiones denominado número de Péclet cuando el coeficiente de difusión no representa una viscosidad cinemática). La expresión 3.3.5 puede resolverse siguiendo el método del modo normal desarrollado por Godunov y Ryabenkii20,21 (1964) para resolver de problemas estacionarios en diferencias finitas (FD). La esencia de este método consiste en buscar soluciones de la forma: 𝑖 𝑖 3.3 . 7 Como resultado de aplicar 3.3.7 sobre 3.3.5 se obtiene una ecuación cuadrática con la forma: 2 3.3 . 8 donde se ha despreciado el efecto del término de la derecha en 3.3.5 . Las raíces de 3.3.8 resultan 1 2 3.3.9 y con ello la solución de 3.3.5 viene dada por la expresión 𝑖 1 2 𝑖 3.3.10 Esta relación manifiesta la inestabilidad inherente a la representación estándar para , dado que el valor de 𝑖 depende de la paridad nodal. Fig. 3.3 - 1 : Inestabilidad introducida por la representación estándar en problemas de convección dominante. Se muestran los resultados del problema: 𝑥𝑥𝑥 con y para diferentes . 3.3 El problema de convección-difusión 65 El efecto de estas oscilaciones en problema de convección dominante puede llegar a distorsionar por completo la solución del problema (véase Fig. 3.3-1). Representación autoadjunta Como hemos visto en la sección precedente (ec. 3.3.3 ), la formulación estándar conduce a matrices de rigidez simétricas, en problemas de difusión. Además, también se ha hecho hincapié en que, en esas condiciones, el método de Galerkin (formulación estándar) es óptimo en la norma de energía18,19. Sin embargo, cuando los términos de convección son significativos esta representación conduce a inestabilidades. Surge entonces la cuestión de si las funciones de variación/peso pueden definirse de tal forma que desaparezca la asimetría introducida por los términos de convección. Para ello es posible multiplicar la formulación débil por una función de “buen comportamiento”, , Ω 3.3.11 donde está definida de acuerdo con 3.2.6 . Resulta entonces la expresión Ω Ω Ω 3.3.12 Esta relación generará matrices simétricas cuando 3.3.13 a partir de lo cual puede definirse la función : − 2 Pe 𝑥 / ℎ 3.3 . 14 Utilizando la expresión 3.3.14 sobre 3.3.12 , y haciendo uso de los subespacios aproximados, se obtiene la ecuación para un nodo interno 2 Pe Pe −Pe 2𝑖−1 2 Pe − 2 Pe Pe −Pe 2𝑖−1 −2Pe Pe − Pe 2 𝑖 + 1 2 3.3.15 Como puede comprobarse, cuando , el nodo aguas arriba ( 𝑖−1) está ponderado por exceso en comparación con la representación estándar. Del mismo modo, el nodo aguas abajo (𝑖+1) se pondera por defecto en comparación con 3.3.5 , por lo que resulta un problema equivalente a una 66 Capítulo 3 | El método: FEM discretización de primer orden en FD para la parte convectiva. Por este motivo es de esperar que representación 3.3.15 sea más estable que su equivalente 3.3.5 . Streamline Upwind Petrov-Galerkin (SUPG) Para entender la dificultades en la formulación estándar es posible establecer una equivalencia entre 3.3.5 y un esquema centrado en FD. Por otra parte, es un hecho conocido que las dificultades que surgen en estos últimos esquemas se deben al efecto del pesado simétrico de los nodos vecinos. Luego, parece que el primer paso para estabilizar el esquema pasa por romper la simetría de las funciones de variación/peso. Para ello es posible establecer una correspondencia con la representación autoadjunta, construyendo las funciones de forma modificadas en el subespacio aproximado ℎ de acuerdo con la expresión: 𝑖 ∗ 𝑖 𝑖 𝑖 + 1 3.3.16 donde la función se define como sigue 𝑖 𝑖 𝑖 con 𝑖𝑖−1 𝑖. Utilizando ahora 3.3.16 para las funciones de peso en 3.3.1 se obtiene una ecuación nodal (para nodos internos) con la forma: 𝑖−1 𝑖 𝑖 + 1 𝑖 3.3.17 o bien, en términos del número de Péclet ( ) 𝑖 − 1 𝑖 𝑖 + 1 2 𝑖 3.3.18 Imponiendo la condición de igualdad entre 3.3.15 y 3.3.18 , se obtiene 2 Pe Pe − P e 2 3.3.19 lo que conduce a la expresión óptima para el factor de peso sobre 𝑖 ∗: óptimo 3.3.20 3.3 El problema de convección-difusión 67 Finalmente, a partir de 3.3.18 puede deducirse la estructura de las funciones de forma modificadas, dada la “similaridad” (en términos de molécula computacionalf) del término estabilizador con un término difusivo. Resulta, entonces, 𝑖 ∗ 𝑖 𝑖 3.3.21 La Fig. 3.3-2 muestra el efecto del término estabilizador 3.3.21 en problemas de convección dominante, analizando el impacto del factor de peso sobre el término estabilizador. Este factor cuantifica en qué medida la redistribución de pesos, entre los nodos aguas abajo y aguas arriba, conduce a un término convectivo evaluado en la dirección del flujo (método streamline upwind, SU ). El valor crítico para puede deducirse a partir de 3.3.18 , imponiendo la condición límite respecto de la influencia del nodo aguas abajo ( 𝑖+1 ), 3.3 . 22 de donde se obtiene f El concepto de molécula computacional se estableció, originalmente, en el contexto de FD, y alude a la discretización de los diferentes términos de la ecuación (convección, difusión, reacción, etc) por medio de una representación gráfica (tipo grafo) de los nodos implicados en la misma así como los correspondientes pesos asociados. Fig. 3.3 - 2 : Efecto del pesado upwind en problemas de convección dominante. Se muestran los resultados equivalentes al problema ilustrado en la Fig. 3.3-1: 𝑥𝑥𝑥 con y para diferentes . 68 Capítulo 3 | El método: FEM crítico 3.3.23 La otra situación límite contemplada en la Fig. 3.3-2 es el caso : upwind Pe → + ∞ 3.3.24 Estrictamente hablando, la definición del factor de peso en 3.3.24 como upwind es un abuso del lenguaje, derivado del uso de esquemas de primer orden del tipo upwind en el contexto original de los esquemas con limitadores en FD (Van Leer11,12 y Boris-Book13,14). En tales casos, la presencia de estos términos permitía preservar el carácter monótono (condición necesaria para la convergencia, Teorema de Godunov8 (1959)) de los esquemas cuando aquellos términos de mayor precisión (órdenes superiores a , en general) inducían inestabilidades. Por esta razón hemos adoptado esta filosofía en la notación, ya que de acuerdo con 3.3.17 , el factor de peso definido en 3.3.24 conduce a un esquema upwind (de primer orden) para el término convectivo. Esto significa que la precisión del término convectivo se reduce de un segundo orden (bajo la configuración original, equivalente a un Fig. 3.3 - 3 : Representación de los diferentes pesos para el término upwind (), y su impacto sobre la estabilidad del esquema. 3.4 Discretización temporal 69 término del tipo leapfrog en FD) a un primer orden (bajo la configuración upwind). Como es sabido en FD, los problemas hiperbólicos bajo esquemas implícitos donde el término convectivo es de tipo upwind (y de primer orden) son incondicionalmente estables. De ahí que se recurra a este tipo de pesado en situaciones localmente difíciles de resolver mediante esquemas de órdenes superiores. La representación de los diferentes en términos del número de Péclet (Fig. 3.3-3) resume el efecto estabilizador de este tipo de funciones de forma. De este modo, la Fig. 3.3-3 muestra una zona inestable para valores de inferiores a crítico, mientras que toda la región por encima de óptimo conduce a un exceso de difusión que se traduce en una pérdida de precisión. 3.4 Discretización temporal La representación estándar (ec. 3.2.14 ) asume un planteamiento prismático9 o de discretización parcial sobre la variable tiempo, . Esto significa que mientras las funciones de forma dan una solución aproximada en función de la discretización espacial del dominio, toda la información temporal queda contenida en los valores nodales. De este modo, 3.2.14 puede reescribirse mostrando la dependencia explícita con y ℎ 𝑖 𝑖 𝑖 𝑀 − 1 𝑖 = 1 3.4.1 Como resultado ello, la ecuación semidiscreta (ec. 3.2.18 ) se convierte en una ecuación diferencial ordinaria (ODE), siendo la variable independiente. Una de las alternativas posibles para la discretización temporal pasa por un tratamiento similar al seguido para la discretización espacial, donde se hace uso de un planteamiento en residuos ponderados. Con ello, al igual que ocurre en la parte espacial, la ecuación semidiscreta se cumple en promedio sobre cada paso de tiempo. Esta es la esencia de los métodos de paso único basados en Backward Difference (BDF) o -Generalized 9,22. En el caso de un esquema de primer orden: 𝑖 𝑖 𝑛 𝑖 𝑛 + 1 𝑖 𝑛 3.4.2 donde 𝑖 𝑛𝑖𝑛 con 𝑛, la expresión en residuos ponderados para la ecuación semidiscreta Δ 𝑡 0 Δ 𝑡 0 𝑖 3.4.3 70 Capítulo 3 | El método: FEM puede transformarse en la expresión 𝑛 + 1 𝑛 𝑛 𝑛 + 1 𝑛 3.4.4 donde Δ 𝑡 0 Δ 𝑡 0 Δ 𝑡 0 Δ 𝑡 0 3.4.5 Por tanto, 3.4.4 puede reorganizarse para obtener el sistema de ecuaciones que define los nuevos valores nodales 𝑛 + 1 − 1 𝑛 3.4 . 6 El algoritmo En la práctica, la relación 3.4.6 no puede resolverse exactamente, por lo que se recurre a un proceso iterativo basado en el residuo numérico de la formulación en cuestión. De esta manera, además de las aproximaciones inducidas por la discretización espacial y temporal (que introducen un error vinculado a la calidad de sendas aproximaciones respecto del comportamiento espacial y temporal de la solución real del problema), debe considerarse que 3.4.6 solo se satisface de manera aproximada debido a la discretización numérica, es decir, el método de resolución numérico seguido en 3.4.6 añade una tercera dimensión a la discretización del problema y con ello un error (numérico) relacionado con el condicionamiento del sistemag 3.4.6 . En este sentido se define el residuo asociado a la formulación 3.4.4 de acuerdo con la relación que sigue: 𝑛 + 1 𝑛 + 1 𝑛 3.4.7 donde, 𝑛 + 1 𝑛 + 1 3.4.8 g En el ámbito del análisis numérico, el condicionamiento de un sistema de ecuaciones mide la “sensibilidad” (en términos de convergencia) del sistema frente a pequeñas variaciones en las variables a resolver. En sentido, asumido el error numérico, se entiende que las soluciones en cada paso de tiempo son del tipo: 𝑛+1 𝑛+1, siendo 𝑛+1 el error numérico del sistema. Este tipo de errores así como la convergencia del método usado para resolver el sistema suele analizarse a partir de lo que se conoce como número de condición, cond −1 donde representa la matriz de coeficientes asociada al sistema de ecuaciones, . 3.4 Discretización temporal 71 define la matriz tangente o Jacobiana del problema. Para resolver el sistema 𝑛+1 se sigue un algoritmo basado en el método iterativo de Newton-Raphson con un predictor lineal: Paso 1. En primer lugar, se calcula una solución aproximada para el instante 𝑛+1 a partir de un predictor lineal: 0 𝑛 + 1 𝑛 𝑛 3.4.9 donde el subíndice que precede a la base indica el índice de la iteración. En cuanto a 𝑛, puede obtenerse de la ecuación semidiscreta dado que 𝑛 satisface dicha ecuación, 𝑛 − 1 𝑛 3.4 . 10 La solución del sistema de ecuaciones 3.4.10 puede llevarse a cabo tanto con resolvedores directos como iterativos, dependiendo de la naturaleza del sistema en cuestión y de los grados de libertad del mismo. Paso 2. A partir de 3.4.7 se obtiene el residuo correspondiente a la primera iteración, esto es, 0 𝑛 + 1 0 𝑛 + 1 3.4.11 Paso 3. Conocido el residuo asociado a la iteración, cuando no supere el criterio de convergencia se calcula el factor de amortiguamiento por medio de éste (ec. 3.4.8 ): 𝑛 + 1 − 1 0 𝑛 + 1 3.4.12 Paso 4. De acuerdo con 3.4.9 el vector de los nuevos valores nodales para resulta 1 𝑛 + 1 0 𝑛 + 1 𝑛 + 1 1 3.4.13 Paso 5. Finalmente se estudia el error sobre el residuo, esto es, 1 𝑛 + 1 3.4.14 a partir de cual se establece: 1 𝑛 + 1 𝑛 + 1 1 𝑛 + 1 3.4.15 donde representa una función escalón o función de Heaviside. Paso 6. Si después de un máximo de iteraciones no se alcanza solución convergente, se reduce el paso de tiempo, , a una fracción de éste de acuerdo con algún criterio (se utiliza por defecto una reducción de un orden de magnitud). 72 Capítulo 3 | El método: FEM 3.5 Precisión y Estabilidad Para entender las dificultades derivadas de las aproximaciones implícitas en el esquema numérico conviene tener presente que, numéricamente, solo es posible aspirar a resolver una ecuación diferencial equivalente8 (EDE). Es por esto que no existe ningún misterio en el hecho de que el comportamiento del sistema, en ocasiones, “se aleje de lo esperado”. Este hecho manifiesta la importancia de cuantificar y predecir el comportamiento de los errores numéricos asociados a un esquema. Para ello se utiliza un análisis espectral del error o método de von Neumann16,17. La esencia de este método reside en la linealidad del esquema numérico. En este sentido, la solución calculada a partir de 3.4.6 puede escribirse como solución aproximada exacta más el error de truncadoh (incluye errores numéricos), esto es, 𝑖 𝑛 𝑖 𝑛 𝑖 𝑛 3.5 . 1 donde 𝑖 𝑛 indica el error asociado al punto de la malla , siendo el índice de tiempo; y 𝑖 𝑛 representa el solución exacta de la ecuación diferencial original bajo las aproximaciones implícitas en la discretización espacial y temporal. De este modo, conforme a la linealidad del esquemai se tiene: 𝑖 𝑛 𝑖 𝑛 𝑖 𝑛 𝑖 𝑛 3.5 . 2 Por tanto, en estas condiciones, estudiar el comportamiento del error equivale a estudiar el comportamiento de la propia solución numérica. Así, para una función arbitraria definida sobre la malla cuya solución es 𝑖 𝑛, la descomposición correspondiente en una serie finita de Fourier viene dada por la expresión: 𝑖 𝑛 𝑙 𝑛 𝑁 𝑙 = − 𝑁 𝑗 𝑘 𝑙 𝑥 𝑖 𝑙 𝑛 𝑗 𝑘 𝑙 ⋅ 𝑖 ℎ 𝑁 𝑙 = − 𝑁 𝑙 𝑛 𝑗𝑖 𝜙 𝑙 𝑁 𝑙 = − 𝑁 3.5.3 h No debe confundirse este error de truncado con el error asociado a la aritmética del ordenador. Este error está aquí referido a todos los términos truncados en la aproximación espacial y temporal. Así por ejemplo, de acuerdo con la discretización espacial hemos visto que los términos convectivos adoptan una estructura equivalente (en FD) a un esquema centrado o bien upwind de primer orden, según sea el valor del factor de peso . Esto significa que en un desarrollo en serie de Taylor, en torno al nodo en cuestión, estaríamos truncando a partir del tercer orden, en el mejor caso. Por este motivo, la ecuación diferencial (discretizada) que resolvemos no sería la original sino una relación equivalente (EDE) que contempla el problema original y los términos truncados. i Se utiliza la letra caligráfica para referirnos al operador (numérico) asociado a la ecuación diferencial del problema en cuestión. Así, tomando como ejemplo el problema de convección-difusión 3.3.1 , equivale a la discretización espacial correspondiente, dada en 3.3.5 . 3.5 Precisión y Estabilidad 73 donde 𝑙 𝑛 representa amplitud del armónico l-ésimo y . En 3.5.3 se han utilizado las relaciones: 𝑙 0 0 máx 𝑙 𝑙 3.5.4 La condición de estabilidad de von Neumann establece que un esquema será estable en tanto que la amplitud de un armónico arbitrario no crezca en el tiempo. Para ello se define el factor de amplificación, , como sigue: 𝑛 + 1 𝑛 3.5.5 Con esta magnitud, la condición de estabilidad puede expresarse como: 𝑙 3.5 . 6 En términos del esquema 3.4.6 , esto significa aplicar la descomposición 3.5.3 sobre el soporte correspondiente j: 𝑇 𝑖+1 𝑛+1 𝑖 𝑛+1 𝑖−1 𝑛+1 𝑇 𝑖+1 𝑛+1 𝑖 𝑛+1 𝑖−1 𝑛+1 Resulta, entonces, el factor de amplificación: 𝜗 2 2 2 2 3.5.7 donde, 2 3.5.8 representan el número de Courant y de difusión, respectivamente. Finalmente, a partir del factor de amplificación pueden definirse dos magnitudes que permiten cuantificar el error asociado al esquema: j Seguimos aquí la notación adoptada por Zienkiewicz9, donde se distingue entre matrices () y vectores o matrices de una columna ( ). 80 Capítulo 3 | El método: FEM quedan bien resueltos, pueden surgir dificultades numéricas motivadas por la separación de carga. Por tanto, la evolución del campo eléctrico es un excelente indicador de las diferentes escalas de tiempo y sus efectos. En este sentido, puede apreciarse (Fig. 3.6-3) un ligero desfase respecto de los resultados obtenidos por medio de FCT en FD. Este diferencia de fase está en consonancia con la calidad de los esquemas FCT en lo que respecta al error de dispersiónb, y los resultados obtenidos para el esquema CN ( ) bajo elementos finitos (Fig. 3.5-4). En cuanto al segundo punto del test (errores de difusión), los resultados mostrados en la Fig. 3.6-4 ponen de manifiesto la importancia de controlar este tipo de errores. Como puede verse en la Fig. 3.6-2, la distorsión del campo eléctrico induce la formación de frentes de choque en las diferentes especies, y con ello el comportamiento del esquema a altas frecuencias es crucial para la estabilidad del problema numérico. En principio, en dicha figura (Fig. 3.6-2) solo se muestran las especies negativas porque el comportamiento “asimétrico” del esquema, respecto de los términos de reacción (generación o destrucción), hace que la estabilidad del problema numérico esté determinada por la evolución de los términos de captura b Los esquemas FCT no genera un error de dispersión en problemas de convección pura con velocidad constante. Fig. 3.6 - 4 : Comparación de la evolución del campo eléctrico y el término de reacción asociado a la captura electrónica ( 2 3 ). En la gráfica de la derecha, el pico que aparece en torno a , tanto para como para , no es un efecto numérico y se debe a la estructura de los coeficientes 2 y 3. Sí lo son, en cambio, las posteriores oscilacion es en el caso de . 3.6 Un test numérico 81 electrónica (esto es, términos de reacción de destrucción). Antes de explicar por qué esto es así, conviene resaltar los resultados de la Fig. 3.6-4, donde se muestra el papel de la ponderación en el tiempo en la estabilidad del esquema. Aunque los resultados mostrados en las figuras Fig. 3.5-1 y Fig. 3.5-3 parecen contradecir los presentados en la Fig. 3.6-4 (es decir, exhibía una mayor atenuación, tanto en la Fig. 3.5-1 como en la Fig. 3.5-3, y con ello se esperaría una mayor estabilidad), en realidad no ocurre así, ya que estas figuras (Fig. 3.5-1 y Fig. 3.5-3) responden al caso de un problema tipo convección-difusión (CD). De este modo, al introducir los términos de reacción, el factor de amplificación resulta (CDR): 𝜗 3.6.4 donde, 2 2 3.6.5 siendo el número de Damköhler, es la frecuencia de colisión asociada al término de reacción en cuestión, la velocidad relativa entre especies y el tamaño de la malla. Las funciones y están en correspondencia Fig. 3.6 - 5 : Error de difusión en esquema de primer orden, CN+SUPG, para ecuaciones CDR. 82 Capítulo 3 | El método: FEM con la relación 3.5.7 , y con ello representa el número de difusión, el número de Courant y En estas condiciones, en la Fig. 3.6-5 se muestran tres situaciones específicas para el número de Damköhler que definen los casos límites para la estabilidad del problema, dado un valor de . Las gráficas para exhiben, en general, una atenuación respecto del efecto real de ionización o captura. No obstante, es significativa la ligera amplificación mostrada en la región de bajas frecuencias por el esquema implícito ( ). En cualquier caso, es de esperar que para ambos esquemas muestren un buen comportamiento. Por otra parte, la gráfica para pone de manifiesto la inestabilidad de los términos de captura debido al efecto de amplificación mostrado tanto para como para . Además esta figura también revela que los esquemas implícitos ( ) son más sensibles a estas inestabilidades casi con independencia de la fase en cuestión. En cambio, los esquemas CN ( ) salvan relativamente bien este tipo de situaciones con elevados, siempre y cuando se mantenga el esquema en la zona de bajas frecuencias. Esta conclusión permite explicar el comportamiento mostrado en la Fig. 3.6-4, donde la inestabilidad de los términos de captura distorsiona significativamente la solución. Falta por demostrar, entonces, que las situaciones representadas en la Fig. 3.6-5 son realmente los casos determinantes en este problema. Para ello, basta reconsiderar la definición del número de Damköhler en los siguientes términos: 𝑔𝑎𝑠 3.6.6 donde representa la sección eficaz del proceso en cuestión (ionización o captura) y 𝑔𝑎𝑠 la densidad del gas ( 2). En las condiciones del problema, a una presión de , se tiene ( 𝑔𝑎𝑠 24 −3): 𝛼 − 19 2 ℎ = 10 − 5 m 𝜂 − 18 2 ℎ = 10 −5 m 3.6.7 Por tanto, con un tamaño de malla −5 , el término de ionización se mantiene en el límite de , por lo que se resuelve perfectamente. En cambio, los términos de captura se hacen inestables ( ) al tomar su valor máximo. Precisamente esta situación tiene lugar cuando el coeficiente de captura 3 predomina sobre 2. La región de colapso del campo eléctrico (donde ) que surge tras la distorsión del mismo se corresponde con dicha zona. Por último, como puede verse en la Fig. 3.6-6, la transición de 2 a 3 genera un pico de alta frecuencia en la zona comprendida entre el 3.7 Un test numérico 83 gradiente máximo de y la región de colapso del campo eléctrico. Esto explica que la presencia del pico en la Fig. 3.6-4 en torno a no sea un efecto numérico y sí las oscilaciones posteriores, tal y como se indicaba en dicha figura (Fig. 3.6-4). 3.7 Algunas conclusiones finales A modo de resumen, los aspectos fundamentales abordados en este capítulo giran en torno a las dificultades numéricas analizadas en las descargas eléctricas, como es el caso del problema propuesto en la sección del test numérico (§3.6). Se organizan estas características en los siguientes puntos: 1. En problemas de convección dominante, la formulación estándar en el FEM introduce inestabilidades cuando . La formulación SUPG permite salvar estas situaciones rompiendo con la simetría espacial de las funciones de forma (Fig. 3.3-2). 2. En la formulación SUPG, el peso asociado al término estabilizador () es crítico en la estabilidad del problema (Fig. 3.3-2). La relación 3.3.20 define el valor óptimo para dicho peso en caso de convección dominante (Fig. 3.3-3). 3. La contrapartida a la estabilidad asociada a la formulación SUPG viene dada por el hecho de reducir la precisión, ya que, en última instancia se transforma Fig. 3.6 - 6 : Representación de las secciones eficaces para el problema de Morrow 15 : ionización ( ; captura electrónica — 2 cuerpos— ( 2 ); captura electrónica —3 cuerpos— ( 3 ). 84 Capítulo 3 | El método: FEM la parte convectiva de una diferencia centrada a un término backward de primer orden. Por este motivo, en problemas transitorios, la difusión artificial introducida por el término upwind puede ser excesiva. 4. Para evitar el exceso de difusión artificial debido al término upwind, se plantea como alternativa, en el caso problemas transitorios, un esquema híbrido entre un CN ( ) y SUPG. Este último, actuando a modo de limitador, permite mantener la calidad del error de difusión del primero filtrando la región de altas frecuencias (Fig. 3.5-3). 5. En situaciones de reacción dominante, la ponderación temporal juega un papel crucial en lo que respecta a la sensibilidad del esquema frente a las discontinuidades introducidas por los términos de reacción (Fig. 3.6-5). 6. En general, la sensibilidad del esquema CN+SUPG no es simétrica en lo que respecta a términos de reacción fuentes o sumideros. Esta diferencia se acentúa con el aumento de (Fig. 3.6-5). Bibliografía 85 Bibliografía 1. Strutt, J. W. On the theory of resonance. Philos. Trans. R. Soc. Lond. 161, 77– 118 (1871). 2. Ritz, W. Über eine neue Methode zur Lösung gewisser Variationsprobleme der mathematischen Physik. (1909). 3. Galerkin, B. G. Series solution of some problems of elastic equilibrium of rods and plates. Vestn Inzh Tekh 19, 897–908 (1915). 4. Courant, R. & others. Variational methods for the solution of problems of equilibrium and vibrations. Bull Amer Math Soc 49, 1–23 (1943). 5. Turner, M. J. Stiffness and deflection analysis of complex structures. J. Aeronaut. Sci. Inst. Aeronaut. Sci. 23, (1956). 6. Clough, R. W. The finite element method in plane stress analysis. (1960). 7. Brooks, A. N. & Hughes, T. J. R. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Eng. 32, 199–259 (1982). 8. Donea, J. & Huerta, A. Finite Element Methods for Flow Problems. (Wiley, 2003). 9. Zienkiewicz, O. C., Taylor, R. L. & Zhu, J. Z. The Finite Element Method: Its Basis and Fundamentals, Sixth Edition. (Butterworth-Heinemann, 2005). 10. Zienkiewicz, O. C., Taylor, R. L. & Nithiarasu, P. The Finite Element Method for Fluid Dynamics, Sixth Edition. (Butterworth-Heinemann, 2005). 11. Van Leer, B. Towards the ultimate conservative difference scheme I. The quest of monotonicity. in Proceedings of the Third International Conference on Numerical Methods in Fluid Mechanics 163–168 (Springer, 1973). 12. Van Leer, B. Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme. J. Comput. Phys. 14, 361–370 (1974). 13. Boris, J. P. & Book, D. L. Flux-corrected transport. I. SHASTA, A fluid transport algorithm that works. J. Comput. Phys. 11, 38–69 (1973). 14. Boris, J. P. & Book, D. L. Solution of continuity equations by the method of flux-corrected transport. Control. Fusion 1, 85–129 (1976). 15. Morrow, R. Theory of negative corona in oxygen. Phys. Rev. A 32, 1799–1809 (1985). 16. Crank, J., Nicolson, P. & Hartree, D. R. A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Math. Proc. Camb. Philos. Soc. 43, 50 (1947). 17. Charney, J. G., FjöRtoft, R. & Neumann, J. Numerical Integration of the Barotropic Vorticity Equation. Tellus 2, 237–254 (1950). 86 Capítulo 3 | El método: FEM 18. Reddy, J. N. Applied Functional Analysis and Variational Methods in Engineering. (McGraw Hill Higher Education, 1986). 19. Brenner, S. C. & Scott, R. The Mathematical Theory of Finite Element Methods. (Springer Science & Business Media, 2008). 20. Godunov, S. K. A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics. Mat. Sb. 89, 271– 306 (1959). 21. Godunov, S. K. & Ryabenkii, V. S. Theory of Difference Schemes. (North- Holland Publ., Amsterdam, 1964). 22. Chung, J. & Hulbert, G. M. A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalized-α Method. J. Appl. Mech. 60, 371–375 (1993). 23. Durán-Olivencia, F. J., Pontiga, F. & Castellanos, A. Multi-species simulation of Trichel pulses in oxygen. XXI Europhys. Conf. At. Mol. Phys. Ioniz. Gases (2012). 24. Durán-Olivencia, F. J., Pontiga, F. & Castellanos, A. Numerical Modelling of Electrical Discharges. J. Phys. Conf. Ser. 490, 012209 (2014). 25. Durán-Olivencia, F. J., Pontiga, F. & Castellanos, A. Multi-species simulation of Trichel pulses in oxygen. J. Phys. Appl. Phys. 47, 415203 (2014). Bibliografía 87 de en CAPÍTULO 4 89 Pulsos de Trichel en Contenido 4.1 Introducción ................................. 90 4.2 Formulación del problema ........................ 93 Construcción de un modelo químico consistente—93 Descripción multifluida: una cuestión de escalas—100 Ecuaciones fluidas. Aproximaciones—102 Efectos de superficie. Interacción descarga-electrodo—104 4.3 Descripción del modelo ..........................106 Sistema de ecuaciones de continuidad—108 Condiciones de contorno—110 Condiciones iniciales—112 4.4 Aproximación numérica .........................112 4.5 Evolución de los pulsos de Trichel ...................114 El primer pulso de Trichel—114 Pulsos de Trichel regulares—121 4.6 Conclusiones ................................131 Apéndice 133 Bibliografía 135 Capítulo —“And you do Addition?” the White Queen asked. “What's one and one and one and one and one and one and one and one and one and one?” —“I don't know”, said Alice. “I lost count”. —“She can't do Addition”, the Red Queen interrupted. LEWIS CARROLL, Through the Looking Glass, 1871. 96 Capítulo 4 | Pulsos de Trichel 2 los iones negativos a tenor de los procesos considerados en la Fig. 4.2-2. El hecho de considerar el ozono en el modelo responde a dos aspectos fundamentales: el primero de ellos, como se muestra en la Fig. 4.2-4–(a), se debe a la estrecha relación que existe entre el proceso de formación de ozono Fig. 4.2 - 2 : Distribución de perdidas energéticas basada en los cálculos de Eliasson y Kogelschatz 26 . Fig. 4.2 - 3 : Modelo simplificado de los niveles de energía para los estados neutros del 2 en la región de Frank- Condon (el diagrama completo puede verse en apéndice, Krupenie 52 ). 4.2 Formulación del problema 97 y el mecanismo de disociación de 2 , siendo este último unos de los procesos fundamentales señalados anteriormente. El otro aspecto está relacionado con la absorción de radiación por parte del 3 . Aunque el modelo no contempla, por sencillez, el efecto que la radiación tiene en el desarrollo de la descarga, el ozono representa una de las principales especies implicadas en la absorción Fig. 4.2 - 4 : (a) Representación esquemática del proceso de formación de ozono en oxígeno y su relación con el mecanismo de disociación (Eliasson y Kogelschatz26 ). (b) Sección eficaz para la absorción de radiación (Eliasson y Kogelschatz26). 98 Capítulo 4 | Pulsos de Trichel 2 energética de la radiación en el rango UV (Fig. 4.2-4–(b)). Además, su poder de oxidación la convierte en una especie de interés en numerosas aplicaciones prácticas. Las consideraciones hasta aquí expuestas sirven para definir las líneas generales del modelo plasma-químico. Sin embargo, éste debe completarse con la inclusión de todos los procesos importantes que involucren a las especies consideradas. La selección de dichos procesos se ha realizado teniendo presente la importancia de sus tasas de reacción. Así mismo, se han considerado los resultados obtenidos en estudios previos de la descarga corona estacionaria19, donde se realiza una selección de aquellas reacciones Tabla 4.2 - 1 : Modelo simp lificado de la cinética química de los pulsos de Trichel en 2. Nº (ref.) Reacción Ct e . de reacción, 𝟑 − 𝟏 R0 27 2 2 R1 28 2 2 + R2 27 2 − R3 29 3 − 2 R4 29 3 2 − R5 29 2 2 − 2 R626 2 1 g 2 2 − 2 − 30 2 R7 27 2 R827 2 2 1 g R926 2 − 2 1 g 2 − 10 R10 26 3 − 3 2 − 12 R11 26 3 − 2 − 13 R12 26 − 2 2 R13 26 − 2 3 − 15 R14 26 − 2 2 − R15 26 2 − 3 3 − 2 − 10 R16 26 − 3 3 − − 10 R17 26 − 2 3 − 2 − 30 2 R18 26 2 + 2 4 + 2 − 30 3 . 2 2 R19 26 4 + 2 2 + 2 − 6 4 − 5030 / 𝑇 R20 26 4 + 2 + 3 − 10 R21 26 2 3 2 − 35 663 / 𝑇 2 R22 26 3 2 − 11 − 2300 / 𝑇 R2326 2 1 g 2 2 − 18 0 . 8 R2426 2 1 g 3 2 2 − 11 − 2840 / 𝑇 4.2 Formulación del problema 99 cuya contribución a la cinética plasma-química es superior al 1%. El resultado de este proceso de selección queda recogido en la Tabla 4.2-1 donde se muestra el conjunto simplificado de reacciones utilizadas para estudiar la cinética química de los pulsos de Trichel en 2 . Dado que algunas reacciones dependen del campo eléctrico, en la Fig. 4.2-5–(a) se han representado las secciones eficaces de los procesos por impacto electrónico Fig. 4.2 - 5 : (a) Secciones eficaces de los procesos por impacto electrónico incorporados en el modelo (LXCat 27–29 ). (b) Constantes de reacción como función del campo eléctrico para R12 y R14 (Eliasson y Kogelschatz 26 ). ( b ) ( a ) 100 Capítulo 4 | Pulsos de Trichel 2 (LXCat27–29). Por otra parte, en la Fig. 4.2-5-(b) se muestran las tasas de reacción R12 y R14 como función del campo eléctrico reducido (Eliasson y Kogelschatz26). Descripción multifluida: una cuestión de escalas Como toda aproximación, al usar la descripción fluida en el caso de las descargas deben establecerse ciertos límites fuera de los cuales el modelo no tiene sentido. Para entender el ámbito de la aproximación fluida en el contexto de los pulsos de Trichel pueden utilizarse algunos de los resultados obtenidos en el capítulo 2. Sin embargo, antes de ir desgranando las implicaciones de estos resultados, conviene rescatar una parte de lo visto en el capítulo 3 en lo que respecta a las dificultades numéricas. Las ecuaciones macroscópicas para la descarga pueden reducirse a un sistema de ecuaciones de continuidad con términos “fuentes” y “sumideros”, donde los términos de flujo se descomponen en una contribución de arrastre (debida al campo eléctrico) y otra de difusión (debida a las colisiones moleculares). Desde un punto de vista matemático esto significa reducir la dinámica de las especies a una relación del tipo convección-difusión-reacción (CDR). En el capítulo 3 se mostró, para un caso concreto, las dificultades numéricas derivadas de la naturaleza de estas ecuaciones. No obstante, desde una perspectiva numérica, todas estas dificultades son una consecuencia de la imposibilidad de implementar la ecuación diferencial ED) sin cometer errores en la discretización. Si pudiéramos resolver la ecuación con , los problemas numéricos derivados de los errores de discretización deberían desaparecer (al margen de inestabilidades inherentes a la ecuación diferencial). En efecto, la condición de consistencia implícita en el método (FEM) garantiza que para se resolvería exactamente la ED. Aunque, otra cosa muy distinta sería el sentido físico de la solución así obtenida. Sin embargo, como se muestra a continuación, la elección del tamaño de la malla (y con ello todo lo que supone en términos de errores numéricos) está determinada por la validez de la aproximación fluida. Por tanto ecuación diferencial equivalente (EDE), discretización espacial y temporal, errores numéricos y aproximación fluida deben tratarse como un todo para integrar una imagen coherente de la descarga. En la Fig. 4.2-6 se muestra una representación esquemática de las diferentes escalas espaciales, la relación entre éstas y las condiciones para una descripción fluida consistente. El concepto de longitud de Debye (véase §2.3, ec. 2.3.15 ) define la distancia a partir de la cual dos partículas dejan de estar correlacionadas vía interacción de Coulomb. Este punto define un límite (inferior) en lo que respecta a la validez de la aproximación de campo medio. A partir de éste es posible simplificar el estudio del campo eléctrico por medio 4.2 Formulación del problema 101 de una distribución de carga “media”, de tal modo que cada una de las partículas cargadas solo experimenta el campo creado por todas las restantes. Como resultado de ello se obtiene una ecuación de Poisson similar a la situación macroscópica (véase §2.3–§2.4, ec. 2.4.3 ). En estas condiciones, siempre que se mantenga el régimen hidrodinámico (véase §2.3, ec. 2.3.22 –2.3.23 ) 𝑐 𝐷 𝑚𝑓𝑝 𝐻 4.2 . 7 una escala intermedia entre 𝐷 y 𝐻 parece un buen punto para definir . De este modo, la escala hidrodinámica puede resolverse adecuadamente y por otra parte la escala asociada a 𝐷 sería imperceptible al modelo fluido. Esto resultaría en una imagen fluida coherente. Sin embargo, en ocasiones la escala hidrodinámica requiere algún refinamiento en la malla capaz de resolver distancias del orden de 𝐷 , como es el caso de los frentes de choque desarrollados en las proximidades del electrodo en §3.6. Este es un punto de contacto con el problema numérico, ya que en tal caso los nodos de un elemento están situados a una distancia de 𝐷 , y con ello el condicionamiento del sistema de ecuaciones podría verse afectado drásticamente. Por este motivo, para seguir usando la aproximación de campo medio en el interior de cada elemento es necesario establecer que el efecto de estas correlaciones es despreciable. Para aclarar esta situación es posible dividir las interacciones dentro de la esfera de Debye en dos tipos (véase §2.3): interacciones de corto Fig. 4.2 - 6 : Descripción esquemática de la aproximación fluida para estudiar la evolución de los pulsos de Trichel en el contexto del régimen hidrodinámico. 102 Capítulo 4 | Pulsos de Trichel 2 alcance (colisiones de “núcleo duro”) e interacciones de largo alcance (interacción de Coulomb). El efecto de estas últimas es el único preocupante en relación con la malla, ya que las colisiones de núcleo duro tienen una escala característica del orden del ngström. En este sentido, la ecuación 2.3.19 muestra que estos modos colectivos pueden despreciarse cuando . En el caso de las descargas, y en el marco de gases débilmente ionizados, se cumple que −4 (véase §2.3), por lo que estos modos colectivos pueden despreciarse siempre que el sistema esté, al menos, en equilibrio termodinámico local (LTE). Por otra parte, en situaciones fuera del equilibrio, la relación 2.3.9 muestra que el efecto de las correlaciones de largo alcance es despreciable, lo que indica el predomino de las colisiones binarias y la validez de la BE en el caso de las descargas. No obstante, está condición está sujeta a la restricción de régimen hidrodinámico. Por tanto, dado que ambos criterios tienen alguna limitación para definir con claridad la validez del modelo fluido con 𝐷, se tomará 𝐷 como límite inferior de la escala numérica. Este hecho es crucial a la hora de elegir el esquema numérico, a tenor de las inestabilidades que pueden introducir los términos de reacción (una discusión detallada del impacto numérico de los términos de reacción en la discretización espacial y temporal del problema puede verse en §3.6). Ecuaciones fluidas. Aproximaciones Las 13 magnitudes que definen la aproximación fluida (método de Grad, véase §2.5) para cada especie, pueden concentrarse en una ecuación CDR bajo determinadas condiciones. A continuación se enumeran cada una de estas simplificaciones: 1. Plasma no magnetizado: Típicamente los valores de corriente eléctrica en el caso de la descarga corona no superan los . En particular, la región de los pulsos de Trichel se caracteriza por corrientes del orden de . En estas condiciones puede prescindirse de los efectos magnéticos en 2.3.18 . 2. Dispersión elástica dominante: En las condiciones de la descargas corona (véase Fig. 4.2-5–(a)), la dispersión elástica constituye el principal mecanismo que contribuye a la fricción neta debido al efecto de las colisiones (véase §2.5). 3. Fuerzas macroscópicas: Los términos de aceleración, 𝛼 , en 2.5.5 pueden despreciarse dado que se supone que tanto el campo eléctrico como cualquier otra fuerza macroscópica varían lentamente en comparación con la dinámica de la descarga (régimen hidrodinámico, véase §2.3 ec. 2.3.22 – 2.3.23 ). Por otra parte el término inercial, 𝛼 𝛼 en 2.5.5 , bajo dispersión elástica dominante y en régimen 4.2 Formulación del problema 103 hidrodinámico, tenderá a ser menor (tanto más a presiones superiores) que el términos de colisión a la derecha de 2.5.5 . Esto es debido a la isotropía que introduce el efecto de las colisiones sobre 𝛼 a escala macroscópica. Como resultado de ello la ecuación para la conservación del momento puede reescribirse con la forma: 𝛼 𝛼 𝛼 𝛼 𝑚 𝛼 4.2 . 8 donde 𝛼 𝛼 𝑚 y 𝛼, representan la masa, carga, frecuencia de colisiones elásticas ( ) y densidad de la especie , respectivamente. 4. Plasma templado (warm plasma30): Esta simplificación supone que la temperatura es suficientemente baja como para despreciar los momentos asociados a la velocidad térmica en comparación con la componente de arrastre. Como consecuencia de ello solo se contempla el efecto de una presión escalar, despreciando así todos los términos no diagonales del tensor de presiones (ec. 2.5.6 ). En estas condiciones, en el límite isotérmico ( 𝛼 𝐵 𝛼 𝛼), la conservación del momento para la especie (ec. 4.2.8 ) puede reescribirse de la forma: 𝛼 𝛼 𝛼 𝛼 𝛼 𝛼 𝛼 4.2 . 9 donde 𝛼 𝛼 𝛼 𝛼 y 𝛼 𝐵 𝛼 𝛼 𝛼 representan el coeficiente de movilidad y difusión de la especie , respectivamente. Finalmente sustituyendo 4.2.9 sobre la ecuación de continuidad correspondiente (ec. 2.5.15 ) se obtiene la ecuación CDR para la especie 𝛼 𝛼 𝛼 𝛼 𝛼 4.2.10 donde 𝛼 𝛼 𝛼 , representa la contribución de arrastre en 4.2.9 . 5. Aproximación de campo medio: En el contexto del régimen hidrodinámico (véase §2.3, ec. 2.3.22 – 2.3.23 , subsección previa y Fig. 4.2-6), el acoplamiento entre el campo eléctrico y las ecuaciones fluidas de las diferentes especies puede realizarse por medio de la aproximación de campo medio (§2.4, ec. 2.4.3 ). Con ello, las diferentes especies solo perciben un campo medio creado por todas las partículas, lo que conduce a una ecuación de Poisson: 2 0 0 𝛼 𝛼 𝛼 4.2.11 siendo el potencial eléctrico y 0 la permitividad eléctrica del vacío. 104 Capítulo 4 | Pulsos de Trichel 2 Efectos de superficie. Interacción descarga-electrodo En general, la interacción que tiene lugar entre el plasma y la superficie de los electrodos es un fenómeno extremadamente complejo31–33, en el que están interrelacionados muchos de los parámetros de la descarga, así como de la naturaleza y del estado de la superficie del electrodo. Por este motivo, la interacción de la descarga con las paredes de los electrodos se simplifica sobre la base de las siguientes consideraciones: 1. En relación con las partículas cargadas, se considera que quedan atrapadas en cuanto impactan con la superficie vía proceso de neutralización de Auger (§2.6, ec. 2.6.16 ). 2. Para los estados excitados ( 21g), se asume que el impacto con la pared conduce a un proceso de desexcitación vía proceso de desexcitación de Auger (§2.6, ec. 2.6.20 ). 3. En lo que respecta a la cinética de los estados adsorbidos, como es el caso de los estados metaestables, se considera que la cinética de Langmuir-Hinshelwood33 es aplicable. Esta aproximación asume que todas las partículas adsorbidas son “acomodadas” en la superficie (en equilibrio térmico) antes de tomar parte en otra reacción. Esta simplificación se complementa con la aproximación de campo medio (en el contexto de reacciones catalíticas en superficies), por la que se supone que las diferentes especies adsorbidas se distribuyen aleatoriamente en la superficie sin que exista ningún tipo de interacción entre éstas. Como resultado de ello, pueden despreciarse todos los efectos relacionados con la dinámica de las especies adsorbidas (movilidad, difusión, desorción y recombinación, etc), asociados a la componente tangencial del flujo. Se interpreta así que la interacción con la superficie queda determinada, exclusivamente, por la componente normal del flujo (véase Fig. 4.2-7). Para implementar las condiciones de contorno se sigue una estrategia similar a la empleada por Hagelaar et al.34, donde se conecta (en el sentido de continuidad) el flujo en el volumen con el flujo en la superficie35 (véase Fig. 4.2-7). Este último puede describirse por medio de la expresión: 𝛼 𝑠 𝜕Ω 𝛼𝑡ℎ,𝛼 𝛼 𝛼 𝛼 𝑠 𝛼 𝛼 𝑠 𝑘 𝑘 𝑠 𝑘 4.2.12 donde 𝛼 𝛼 (siendo 𝛼 el coeficiente de reflexión), 𝑡ℎ,𝛼 la velocidad térmica de la especie , 𝑘 el coeficiente de emisión secundaria de la especie a partir de alguna reacción entre la especie y la superficie del electrodo, 4.2 Formulación del problema 105 𝑠 el vector normal a la superficie apuntado desde el volumen a la superficie y la función 𝛼 𝑠 𝛼 𝑠 4.2.13 Los dos primeros términos en el corchete a la derecha de 4.2.12 representan el flujo de partículas que impactan con la pared, donde la velocidad se ha descompuesto en su contribución térmica y de arrastre. Esta última refleja el flujo de partículas a través de la pared con velocidad media 𝛼 , mientras que la componente térmica establece la contribución aleatoria dentro de una región delimitada por el recorrido libre medio desde la pared. El último término del corchete a la derecha de 4.2.12 representa el flujo por difusión a través de la pared del electrodo. El factor se debe al hecho de que el proceso de difusión en la pared solo tiene lugar en el semiplano que contiene al volumen de la descarga (Fig. 4.2-7). Finalmente el término sumatorio en lado derecho de 4.2.12 está asociado con las reacciones de superficie de otras especies que conllevan emisión secundaria de la especie . De acuerdo con las simplificaciones consideradas en la subsección previa, el flujo de partículas en el volumen puede expresarse como sigue (ec. 4.2.9 ) Fig. 4.2 - 7 : Representación esquemática de los procesos de superficie. En la parte derecha el cuadrado negro indica el punto de la malla más cercano a la pared del electrodo, mientras que el círculo negro muestra el punto límite a partir del cual se activa un proceso de Auger (neutralización o desexcitación).