Modelado y estudio fluidomecánico de los actuadores de plasma SDBD

Equation Chapter 1 Section 1 Proyecto Fin de Carrera Ingeniería Aeroespacial Modelado y estudio fluidomecánico de los actuadores de plasma SDBD Autor: Pablo Dugo Domínguez Tutor: Miguel Pérez-Saborid Sánchez-Pastor Dpto. Ing. Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 ii Proyecto Fin de Carrera Ingeniería Aeroespacial Modelado y estudio fluidomecánico de los actuadores de plasma SDBD Autor: Pablo Dugo Domínguez Tutor: Miguel Pérez-Saborid Sánchez-Pastor Profesor titular de Universidad Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 iv Proyecto Fin de Carrera: Modelado y estudio fluidomecánico de los actuadores de plasma SDBD Autor: Pablo Dugo Domínguez Tutor: Miguel Pérez-Saborid Sánchez-Pastor El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2019 El Secretario del Tribunal vi A mi familia y amigos Agradecimientos Hace cuatro años, aproximadamente en esta misma fecha, daba mis primeros pasos en la Escuela Técnica Superior de Ingeniería, siendo consciente del camino que estaba tomando, de los sacrificios a los que iba a tener que someterme, al esfuerzo que tenía que estar dispuesto a hacer por mis sueños, por mi futuro. Me prometí a mí mismo y a mi familia y amigos que haría todo lo posible por cumplir mis objetivos, que daría el 100%, que todo el esfuerzo, todas las horas de estudio, todas las noches en vela, toda privación… tendrían su recompensa. Dentro del tópico sobre la acción y reacción, la descripción que se hace en Interstellar, la película que motivaría mis estudios en Ingeniería Aeroespacial, de la Tercera Ley de Newton es esta: ‘’La única forma que conocen los humanos de avanzar es dejar algo atrás’’. Jamás imaginé que una frase me acompañaría durante tanto tiempo, a veces estrechando mis perspectivas, a veces ensanchando mis deseos y ambiciones. Pero, como muchas personas dicen, tras la tormenta siempre viene el Sol y la calma y durante este largo y arduo camino, siempre encuentras personas dispuestas a ayudarte, incluso sin ofrecer nada a cambio. Sin lugar a dudas, he tenido la suerte de conocer personas muy valiosas que siempre llevaré conmigo dentro, no sólo compañeros de clase sino amigos, amigos de verdad. Gracias a ellos y a ellas este camino ha sido mucho más fácil; sin ese ánimo cuando más lo necesitaba, sin esa palmada en la espalda que tanto anhelaba, sin esos mensajes de: ‘’Ánimo, tú puedes con todo’’, ‘’La suerte es sólo para quién la necesita’’, ‘’Esto no es nada para ti’’… nada habría sido posible o, al menos, no con este resultado. Y, de la misma forma que conoces a nuevas personas que antes no formaban parte de tu vida, conoces realmente a los que ya estaban en ella: a mis amigos y amigas de toda la vida, a mi familia, a mis padres, a mis hermanos…, GRACIAS, de todo corazón, por haber estado ahí siempre. Sin duda, las personas que han aportado algo positivo a mi día a día en estos últimos cuatro años, estén o no en estos momentos siendo parte de mi vida, son el regalo más valioso que me llevo. Y es que, como el viento, las personas vienen y van, pero algunas, aquellas que están destinadas a estar contigo, nunca te fallan. Gracias a mis mejores amigos Josema, Migue y Andrés, a quienes he conocido durante la carrera, gracias por esa ayuda desinteresada y por todos los grandes momentos vividos durante todo este tiempo, sé que me llevo tres grandes amigos para toda la vida; gracias a mis hermanos, Víctor y Pilar, por darme esos ánimos en los momentos más delicados, por sacar siempre mis virtudes; gracias a mis padres, Toni y Pilar, por ser los cimientos de mi vida; gracias al resto de mi familia y amigos, por siempre confiar en mí. Gracias, en especial, a mi compañero de clase, compañero de Trabajo de Fin de Grado y amigo, Migue ,por este gran equipo que hemos formado durante la realización de este proyecto, por abrirme los brazos con una sonrisa cuando salió la propuesta de hacerlo juntos, codo con codo. Gracias también a nuestro tutor, Miguel Pérez- Saborid, por toda la ayuda prestada y el tiempo dedicado, con el que hemos podido entablar una magnífica relación y seguro que en un futuro no muy lejano, una buena amistad. Por último, a la Universidad de Sevilla, al personal no docente, a todos los profesores que han aportado su granito de arena, gracias por todo. Pablo Dugo Domínguez Sevilla, 2019 viii Resumen Los beneficios aerodinámicos que suponen el control del flujo que incide, por ejemplo, sobre la superficie de un perfil aerodinámico son múltiples: retraso de la transición de régimen laminar a turbulento, retraso del punto de separación, reducción de la resistencia aerodinámica, aumento de la sustentación… La obtención de dichos resultados ha sido objeto de estudio en los últimos años por investigadores, quienes encontraron una solución sencilla, barata y ligera para realizar estrategias de control de flujo activo: el uso de actuadores de plasma. Los actuadores de plasma son dispositivos eléctricos capaces de acelerar un flujo por medio de la creación de un campo eléctrico. En concreto, los actuadores de plasma de Descarga con Barrera Dieléctrica (DBD) han sido los mejores posicionados debido a sus ventajas: una corta respuesta en el tiempo, su bajo consumo de potencia, ausencia de partes móviles (baja complejidad mecánica), bajo peso, fáciles de constuir y de adaptar a la superficie que queremos estudiar es la carta de presentación de este dispositivo. Para poner en evidencia las ventajas de los actuadores DBD, existen en el mercado múltiples códigos comerciales, pero su complejidad, tanto matemática como computacional se escapa de nuestra zona de estudio. Para obtener resultados de una forma más sencilla, sin perder la exactitud y la precisión que se requieren, se ha desarrollado un modelo matemático mediante el sistema de cómputo numérico MATLAB basado en el modelo de Suzen-Huang (SH). Dichos resultados, además, serán correctamente contrastados con publicaciones e investigaciones precedentes para corroborar la validez del estudio. La ejecución de nuestro modelo se llevará a cabo, además, para distintas configuraciones físicas del actuador, maximizando aquellos parámetros que mejoran el comportamiento de nuestro actuador DBD con el fin de abarcar al máximo el problema en cuestión. Abstract The aerodynamic benefits of controlling the Flow that affects, for example, the Surface of an airfoil are multiple: delay of transition from laminar regime to turbulent regime, delay of separation point, reduction of drag, increase of lift…The obtaining of these results has been an object of study in recent years by researchers, who have found a simple, cheap and lightweight solution to perform active flow control strategies: the use of plasma actuators. Plasma actuators are electrical devices capable of accelerating a flow through the creation of an electric field. Specifically, the Dielectric Barrier Discharge (DBD) actuators have been the best positioned due to their advantages: a short response time, low power consumption, absence of moving parts (low mechanical complexity), low weight, easy to build and adapt to the surface we want to study… is the letter of introduction of this device. In order to highlight the advantages of DBD actuatores, there are multiple comercial codes on the Marketplace, but their complexity, both mathematical and computational, escapes from our study área. To obtain results in a simpler way, without losing the accuracy and the precision that is required, a mathematical model has been developed through the numerical computation system MATLAB based on the Suzen-Huang (SH) model. In addition, these results will be contrasted with previous publications and research to corroborate the validity of the study. Besides, the execution of our model will be carried out for the different physical configurations of the actuator, maximizing those parameters that improve the behaviour of our DBD in order to cover the problema as much as posible. xvi Notación 𝑛𝑒 Densidad de electrones libres 𝑛𝑖 Densidad de iones positivos 𝑛0 Densidad de átomos o moléculas neutros 𝑇𝑒 Temperatura de los electrones o temperatura cinética media 𝑇𝑖 Temperatura de los iones 𝜆𝐷 Longitud de Debye 𝜀0 Permitividad del vacío 𝑘 Constante de Boltzmann 𝑒 Carga del electrón 𝛼 Grado de ionización 𝐷𝐵𝐷 Dielectric Barrier Discharge MEE Multiple Encapsulated Electrode 𝜌𝑐 Densidad de carga neta en el plasma 𝜌𝑐,max Densidad de carga máxima en el plasma 𝜙𝑀𝐴𝑋 Amplitud del voltaje aplicado 𝐸 󰇍 Campo eléctrico Ex Componente x del campo eléctrico E𝑧 Componente z del campo eléctrico Φ Potencial eléctrico 𝜀 Permitividad 𝜀𝑟 Permitividad relativa del medio 𝜙 Potencial debido a campos eléctricos externos 𝜑 Potencial debido a la densidad de carga en el plasma 𝐺 Distribución semi-Gaussiana 𝑓 Variación temporal del voltaje aplicado 𝜇 Parámetro de localización del máximo de la campana Gaussiana 𝜎 Parámetro de ratio de decrecimiento de la distribución Gaussiana 𝜌 Densidad del fluido 𝜇 Viscosidad dinámica del fluido 𝑉 󰇍 Velocidad del fluido 𝑉𝑥 Componente x de la velocidad del fluido 𝑉𝑧 Componente z de la velocidad del fluido 𝑝 Presión del fluido 𝐹𝑒 󰇍 󰇍 󰇍 Fuerza electrohidrodinámica 𝑅𝑒 Número de Reynolds 𝐷𝑐 Fuerza Coulombiana adimensional T Período de la frecuencia del voltaje aplicado 𝑓𝑐 Frecuencia del voltaje aplicado t Tiempo real Ψ Función de corriente u Componente x de la velocidad adimensional v Componente z de la velocidad adimensional 𝑢0 Velocidad de referencia 𝜉 Componente x adimensional 𝜂 Componente z adimensional 𝜏 Tiempo adimensional 𝐸∗ Campo eléctrico adimensional E𝜉 Componente x del campo eléctrico adimensional E𝜂 Componente z del campo eléctrico adimensional 𝜌𝑐∗ Densidad de carga adimensional 𝜙∗ Potencial debido a campos eléctricos externos adimensional 1. INTRODUCCIÓN En esta sección, primero introducimos los antecedentes que llevan al estudio de los actuadores de plasma como herramientas para el control del flujo. Posteriormente, incluimos una breve síntesis acerca de los objetivos que se desean alcanzar y la motivación que nos lleva al estudio de los actuadores de plasma. Finalmente se realiza una descripción exhaustiva de la estructura del trabajo que en estas líneas se presenta. 1.1 Antecedentes Es bien conocido en campos como la física o la química que cuando modificamos las condiciones de temperatura y presión de cualquier sustancia, pueden obtenerse distintos estados o fases llamados estados de agregación de la materia. Dichos estados están, como no podía ser de otra forma, íntimamente relacionados con las fuerzas de cohesión de las partículas que constituyen la materia. También sabemos que estos estados de agregación tienen propiedades y características dispares. De forma natural son conocidos cuatro estados o fases: sólido (por ejemplo, un cubito de hielo), líquido (el agua de los mares y océanos), gaseoso (el vapor de agua utilizado en una central nuclear) y plasmático o de plasma (aurora boreal). Figura 1. Estados de agregación de la materia. Fuente de la imagen: [20] Tal y como hemos señalado, el plasma es considerado desde hace años como el único estado cuya existencia no puede darse de forma aislada, como es el caso de los otros tres, pues debe existir ineludiblemente alguna forma de campo que lo ionice, como puede ser la temperatura o la electricidad. Los usos del plasma son más que cotidianos hoy en día, desde lámparas hasta televisores, pasando por tubos fluorescentes y otros tipos de luces. También es importante remarcar la aparición natural del estado de plasma, como puede ser en algunas capas de la atmósfera, los rayos en una tormenta o el fenómeno de la aurora boreal que frecuentemente se produce en los polos. La aplicación del plasma ha motivado al desarrollo de tecnologías como la propulsión iónica para vehículos espaciales, concebida desde hace 50 años, aunque no ha sido hasta la fecha actual cuando se ha convertido en una realidad [1]; también se ha utilizado ampliamente en los últimos años en campos como la magnetohidrodinámica (MHD) [2]. Con el tiempo, estas aplicaciones han florecido en el marco de la aerodinámica o la mecánica de fluidos, que son los temas con los que más se relaciona este trabajo. 1. Introducción 12 12 Dentro de este campo, se llevan estudiando desde hace ya varios años los denominados actuadores de plasma para el control de las superficies de las aeronaves. A diferencia de los actuadores servo-mecánicos típicos de los aviones, los actuadores de plasma son un tipo de actuadores aerodinámicos de control de flujo basados en la formación de plasma. Éstos son también conocidos como actuadores de descarga de barrera dieléctrica o DBD – Dielectric Barrier Discharge Plasma Actuator-. Estos novedosos actuadores se llevan probando eficazmente desde 2005 y los resultados arrojados de los estudios son realmente prometedores, por lo que su incorporación para la mejora de la sustentación de superficies aerodinámicas está cerca de ser una realidad [3]. También han sido probados como generadores de vórtices con el objetivo de retardar la separación del flujo aerodinámico y así mantener durante más tiempo el régimen laminar [4]. Por último, otras posibles aplicaciones de estos actuadores están relacionadas con la reducción de la resistencia aerodinámica [5] o incluso como una alternativa a los alerones o los flaps [6]. 1.2 Motivación del problema y objetivos La manipulación del flujo, por ejemplo, sobre las alas de los aviones, han generado numerosos estudios [7], pues es una solución que resulta muy beneficiosa. Dichos estudios se centran básicamente en el control de la capa límite. El control exitoso de esta región permite sacar un mayor rendimiento a la eficiencia aerodinámica, aumentando la sustentación y disminuyendo la resistencia aerodinámca, permitiendo velocidades de despegue y aterrizaje menores, pistas más cortas, mejor maniobrabilidad en entornos desfavorables y por supuesto, un menor consumo de combustible. Resulta interesante, por tanto, el estudio de mecanismos que permitan la manipulación del flujo, como pueden ser el uso de actuadores de plasma. Además de las ventajas señaladas con anterioridad, es preciso valorar a los actuadores de plasma como una alternativa viable a los dispositivos mecánicos de control, o incluso el uso junto a éstos con el objetivo de garantizar la seguridad mediante la redundancia de sistemas y que permitirá, además, la reducción de peso y el alivio de la complejidad mecánica que suponen los servos de control. De esta forma el objetivo es desarrollar un modelo de actuador DBD que sea capaz de simular su actuación, desde el punto de vista hidrodinámico, con el fin de analizar su comportamiento frente a diferentes parámetros y condiciones de contorno. Ello se realizará de una forma sencilla pero compacta desde el entorno del software matemático MATLAB, mediante un modelo simple, transparente y flexible que permita, sin grandes recursos computacionales, obtener resultados fiables y reales del comportamiento del dispositivo. Resumidamente, la finalidad de este estudio es: • Descripción de los principios y fundamentos básicos de la física que envuelve a los actuadores DBD. • Explicación de las ecuaciones que gobiernan el comportamiento de los actuadores DBD. • Presentación de la arquitectura de los actuadores DBD. • Simulación del comportamiento de un actuador DBD desde el punto de vista hidrodinámico. • Comparación de los resultados obtenidos de la simulación con investigaciones antecesoras. Es importante remarcar que, debido a la amplitud del modelado del actuador DBD en cuestión, este estudio se encuentra dividido en dos partes bien diferenciadas, aunque complementarias: el enfoque eléctrico del actuador DBD, realizado por mi compañero y amigo Miguel Ángel Domínguez Jiménez es desarrollado en [43], en el que se estudia la densidad de carga, el potencial eléctrico y las fuerzas electrohidrodinámicas que aparecen y el enfoque hidrodinámico, que es del que trata este estudio, encargado de la implementación del problema fluidomecánico y el estudio del campo de velocidades inducidas por el actuador. 1.3 Estructura del trabajo El estudio está dividido en cinco capítulos y un anexo: ▪ Capítulo 1: Introducción. En esta sección se realiza una descripción de los principios básicos del plasma y los fundamentos teóricos de la física que envuelve a los actuadores DBD, mostrando la arquitectura del dispositivo y su aplicación en el proceso de control del flujo, así como las ventajas y las limitaciones de dichos dispositivos. Se divide en las subsecciones: 1.1 Antecedentes, 1.2 Motivación del problema y objetivos, 1.3 Estructura del trabajo y 1.4 Principios básicos. ▪ Capítulo 2: Modelado del problema. En esta sección se exponen las ecuaciones que rigen el comportamiento de los actuadores DBD, centrándonos tanto en el problema eléctrico como en el problema hidrodinámico. Se divide en subsecciones: 2.1 Modelo de Suzen-Huang 2.2. Formulación del problema eléctrico, 2.3 Formulación del problema hidrodinámico • Capítulo 3: Simulación del actuador de plasma DBD. En esta sección realiza una explicación exhaustiva de los métodos numéricos utilizados para la realización del código numérico en MATLAB que simula la actuación de nuestro actuador DBD. Se divide en las siguientes subsecciones: 3.1 El método de colocación, 3.2 Resolución numérica del problema eléctrico y 3.3 Resolución numérica del problema hidrodinámico. ▪ Capítulo 4: Análisis de resultados. Se recogen y analizan en esta sección los resultados obtenidos de las diferentes simulaciones, comparándolas con investigaciones predecesoras que nos permitan, no sólo validar el estudio, sino extender su flexibilidad a futuras investigaciones. En la última parte de este capítulo se analiza el coste computacional necesario para llevar a cabo nuestras simulaciones. Se divide en las siguientes subsecciones: 4.1 Validación del código del problema hidrodinámico, 4.2 Influencia del resto de parámetros en el problema hidrodinámico, 4.3 Análisis del viento iónico producido y 4.4 Coste computacional. ▪ Capítulo 5: Conclusiones, posibles mejoras y líneas de estudio futuras. En esta sección se explican las conclusiones alcanzadas. Además, se analizará cualitativamente el estudio teniendo en cuenta los objetivos iniciales propuestos para verificar que, efectivamente, se ha abordado completamente el problema en cuestión. Posteriormente, se discuten y se presentan posibles mejoras en la metodología utilizada en el estudio que podrían enriquecerlo y hacerlo más completo. Finalmente, se efectúan una serie de propuestas de futuro obtenidas a raíz de las diferentes investigaciones llevadas a cabo, donde se establecen posibles líneas de investigación para profundizar aún más en el campo de los actuadores de plasma, generando nuevas preguntas e ideas que podrían abrir interesantes vías de trabajo. Se divide en las siguientes subsecciones: 5.1 Conslusiones y 5.2 Posibles mejoras y líneas de estudio futuras. ▪ Anexo: programas de MATLAB. Por último, en este bloque quedan recogidos los códigos numéricos elaborados en MATLAB. Se divide en las siguientes subsecciones: Anexo 1: Cálculo de la distribución de densidad de carga, Anexo 2: Cálculo de la distribución de potencial eléctrico y Anexo 3: Funciones auxiliares. 1. Introducción 14 14 1.4 Principios básicos 1.4.1 ¿Qué es el plasma? En este apartado se pretende realizar una descripción exhaustiva de lo que conocemos como plasma o estado plasmático con el fin de que el lector se adentre en el fundamento básico que compone a los actuadores de plasma, que es el asunto del que trata el trabajo. Los siguientes apartados tratarán de explicarnos cómo se descubrió, cómo funciona y en qué condiciones se presenta el estado de plasma. Antes de ellos, es preciso dar una breve definición del tema a tratar: el plasma. Se conoce al plasma en campos como la química o la física, al estado de la materia en el que la mayoría de los átomos se encuentran ionizados, existiendo una cierta cantidad de electrones que tienen la capacidad de moverse libremente sin encontrarse ligados a ninguno de estos átomos. Es decir, es un fluido formado por iones positivos y por electrones, concentrados en la misma cantidad, además de átomos neutros. Este es el motivo por el cual el plasma es un gran conductor eléctrico y responde de manera eficaz en presencia de campos electromagnéticos 1.4.1.1 Historia del plasma El concepto de plasma data aproximadamente de la tercera década del siglo XX de la mano del ingeniero metalúrgico, físico y químico estadounidense Irvin Langmuir (1881-1957). Además, unos años más tarde, en 1932, sería galardonado con el Nobel de Química. Sin embargo, las evidencias de que un nuevo tipo de energía o estado existían eran anteriores [15]: • El ámbar y su electrificación (1698): En 1698 un científico originario de Inglaterra descubrió la aparición de una chispa cuando estudiaba la electrificación del ámbar. Frotando con asiduidad dicho material, reveló una pequeña descarga eléctrica en el aire, hecho que solo podia producirse mediante la creación de una cierta cantidad de cargas eléctricas libres en el aire como para que éste se convirtiera en un gas conductor: un plasma • Benjamín Franklin, el experimento del rayo y la cometa (1752): En 1752, el científico e inventor estadounidense Benjamín Franklin (1706-1790) hizo un experimento con el que quería demostrar la naturaleza eléctrica de los rayos. Años atrás había estado trabajando en experimentos cuyo fundamento principal era la electricidad. Aquel día una tormenta azotaba en el lugar en que se encontraba Franklin, por lo que se le ocurrió volar una cometa hecha de varillas metálicas y seda para comprobar que los rayos estaban repletos de carga eléctrica. Sus evidencias no estaban tan equivocadas y cuando echó a volar la cometa ésta comenzó a atraer los rayos. Dicho experimento le permitió demostrar que los rayos, de naturaleza eléctrica, podían ser atraidos por un material metálico y que esto podia utilizarse para el uso civil: así nació el pararrayos. • Michael Faraday: descargas en arco (1830): Desde hacía ya una década, allá por 1830, el científico inglés Michael Faraday (1791-1867) había estado estudiando las descargas eléctricas. En uno de sus experimentos descubrió que en las descargas eléctricas en arco había zonas donde se producía una reducción de la presión, y que ésto traía consigo la emision de una luz tenue. • William Crookes: descargas eléctricas en gases (1880): A finales del siglo XIX, el físico inglés William Crookes (1832-1919) descubrió, mientras experimentaba con descargas eléctricas en gases que la zona del gas en que se producía la descarga se comportaba de forma diferente a un gas normal. Nombró al plasma como ‘gas radiante’. También identificó a los electrones, que en esa época eran todavía desconocidos como ‘rayos catódicos’. • Joseph J. Thompson: desviación de rayos catódicos (1897): Casi a principios del siglo XX, Thompson (1856-1940) realizó un conjunto de experimentos en los que midió la desviación de lo que en aquella época se conocía como rayos catódicos cuando se utilizaban campos electromagnéticos. Dichos experimentos le permitieron descubir a los electrones y proponer el primer modelo atómico. Años más tarde, fue galardonado con el Premio Nobel de Física (1906). • Irving Langmuir: descargas eléctricas en gases (1929): Las investigaciones de descargas eléctricas en gases iniciadas por Irving Langmuir en 1923 y junto con el conocimiento de que éstas ionizaban a los átomos del gas dieron sus frutos en 1929, cuando el término de plasma se incluía en el informe de un trabajo donde se estudiaba la nube de electrones que se veía en el interior del gas durante una descarga. Esta nube de electrones brillaba y se movía de una forma gelatinosa que a Langmuir le recordaba al plasma de la sangre. • Lev Landau: teoría estadística del plasma (1936) y los estudios predecesores: No sería hasta finales del primer tercio del siglo XX cuando el físico soviético Lev Landau (1908-1968) desarrollara una teoría estadística que describiera el comportamiento y el origen del plasma. Posterior a Landau, el científico sueco Hannes Alfvén (1908-1995) desarrolló las ecuaciones que describen el movimiento de un fluido eléctrico al ser sometidos a un campo magnético. Tras Alfvén, de nuevo Landau describió matemáticamente la interacción entre las partículas y las ondas en un plasma. Se considera que la física de plasmas modernas tiene su origen en estos trabajos. De hecho, los trabajos de Alfvén acerca del plasma le valieron para ser galardonado con el premio Nobel en 1970. De esta forma la física de plasmas ha tenido un enorme impulso en las últimas décadas y aún quedan muchos asuntos que estudiar, constituyendo una de las áreas de mayor dificultad en al física de hoy. 1.4.1.2 Definición de plasma Un plasma se trata de un gas energizado hasta el punto en que iones positivos de átomos y sus respectivos electrones liberados se encuentran confinados en el mismo espacio. Cuando un material alcanza este estado se dice que está en estado de plasma. De esta forma, el plasma no es otra cosa que un gas ionizado en el que existe un movimiento libre entre iones y electrones. Figura 2. Diferencias entre gas y plasma. Fuente de la imagen: [21] Además, encontramos las siguientes características en los plasmas: • Carecen de forma y volumen fijos. • Son menos densos que los sólidos y los líquidos. 1. Introducción 16 16 • Todos o algunos de los electrones han sido liberados de sus átomos. • Formados por iones positivos y electrones que poseen libertad de movimiento. • Poseen carga neta cero (en su totalidad) pero no a nivel de partículas individuales (nivel local). • Son altamente influenciables por campos electromagnéticos externos. • Tienen la habilidad de crear campos magnéticos dentro de sí mismos debido a la atracción electrostática entre los iones positivos y los electrones, que influyen en la dinámica de las partículas, creando un comportamiento colectivo entre iones y electrones. • Son conductores de la electricidad. Por otro lado, es importante remarcar los procesos más importantes que tienen lugar en un plasma [15]: • Ionización La ionización es el fenómeno desencadenante del plasma, en el cual se produce la formación de iones (aniones y cationes). Como consecuencia, los electrones liberados en la ionización pueden ser acelerados por campos electromagnéticos, que, a su vez, a través de sucesivos choques, pueden producir una ionización en cadena, estableciendo una corriente eléctrica. Figura 3. Proceso de ionización. Fuente de la imagen: [15] Otros procesos derivados de la ionización pueden ser, por un lado, la excitación atómica interna por impacto de electrones y, por otro lado, la desexcitación y emisión de fotones, cuyo análisis permite conocer las especies que contiene. Figura 4. Excitación y Desexcitación de átomos. Fuente de la imagen: [22] • Disociación molecular: Mediante este proceso las moléculas se rompen liberando átomos y radicales libres que reaccionan químicamente y forman a su vez nuevas especies. Es decir, un plasma es un medio muy reactivo químicamente. Figura 5. Disociación molecular. Fuente de la imagen: [23] • Reacciones con las superficies: Puesto que las partículas del plasma están cargadas con mucha energía, son capaces de arrancar partículas de las superficies cercanas y que se incorporan al plasma. Al mismo tiempo, partículas del plasma pasan a depositarse en las superficies. Aprovechando esta idea se están haciendo investigaciones para el uso del plasma en la limpieza de superficies. Figura 6. Limpieza con plasma. Fuente de la imagen: [24] 1.4.1.3 Tipos de plasma Los diferentes tipos de plasma suelen clasificarse en función de dos parámetros fundamentales [14]: • Densidad de electrones: mide la cantidad de electrones y de iones positivos que hay en un volumen de plasma. En sí es una medida del grado de ionización del plasma. • Temperatura: mide el grado de agitación térmica de las partículas. Está relacionada con la energía cinética de éstas. Debido a la temperatura distinguimos entre plasmas no térmicos o fríos, que sólo tienen una pequeña fracción de sus átomos ionizados (bajo grado de ionización) y plasmas calientes, que son creados debido al aumento de temperatura de un gas hasta que los electrones están tan energizados que pueden liberarse de sus átomos. Podemos observar en la siguiente imagen una clasificación de los diferentes tipos de plasma en función de la temperatura y la densidad de electrones. 1. Introducción 24 24 1.4.2.2 Arquitectura del actuador La configuración mecánica del actuador aerodinámico de plasma se muestra en la siguiente imagen: Figura 15. Arquitectura de un actuador DBD. Fuente de la imagen: [4] El actuador de plasma DBD [4] consiste en un conjunto de electrodos finos (en nuestro caso son dos electrodos hechos de cinta de lámina de cobre) dispuestos a lo largo de una superficie. Uno de los electrodos está expuesto al aire, mientras que otro está encapsulado en dieléctrico, en este caso Kapton. Los electrodos están además desplazados uno sobre el otro en la dirección del flujo. Entre dichos electrodos se aplica una diferencia de potencial mediante corriente alterna que varía con el tiempo. El hecho de que se utilicen voltajes de corriente alterna es debido a que una corriente continua no es capaz de atravesar el dieléctrico. Por otro lado, el uso de material dieléctrico se debe a que producen limitaciones al paso de la corriente y evitan que se produzcan arcos eléctricos. Los parámetros más importantes que intervienen en la disposición de los actuadores DBD son: • Longitud y espesor de los electrodos • Espesor de dieléctrico • Desplazamiento entre los electrodos (gap) • Tipo de dieléctrico • Voltaje y frecuencia de la fuente de alimentación La optimización de dichas variables ha sido extensamente estudiada por investigadores [10]. Por otro lado, además de los diferentes valores que pueden darse a los parámetros anteriores, multitud de configuraciones físicas pueden darse, tal y como se observa en las figuras 16, 17 y 18: Figura 16. Arquitectura de un único actuador DBD. Fuente de la imagen: [4] Figura 17. Arquitectura de un actuador DBD doble. Fuente de la imagen: [4] Figura 18. Arquitectura de multiactuadores DBD basados en cuatro actuadores DBD únicos. Fuente de la imagen: [4] Por otro lado, también es interesante destacar lo que se denomina Actuador de plasma de electrodo encapsulado múltiple (MEE) [7]: 1. Introducción 26 26 Figura 19. Arquitectura de diferentes actuadores de plasma MEE a)MEE5, b)MEE6, c)MEE7, d)MEE8. Fuente de la imagen: [7] La diferencia de este actuador respecto a los anteriores es que, en lugar de tener un único electrodo encapsulado uniforme, tiene múltiples electrodos encapsulados que se distribuyen a lo largo de la capa dieléctrica. La profundidad máxima, el ancho total y la polaridad del electrodo encapsulado permanecen iguales en todas las configuraciones. En concreto, las configuraciones MEE constan de dos electrodos encapsulados. En la siguiente tabla se pueden observar los datos de las siguientes configuraciones: Tabla 2. Parámetros de los diferentes actuadores MEE. Fuente de la imagen: [7] Por último, es importante hablar de la función que cumplen los dieléctricos en el actuador de plasma. El material dieléctrico utilizado en la fabricación de actuadores de plasma presenta un rol importante en la eficiencia del dispositivo y también en su durabilidad. Además, juega un papel fundamental en la actuación del propio actuador. Se denomina dieléctrico a un material que tiene una baja conductividad eléctrica y que se caracteriza por formar dipolos eléctricos cuando es excitado mediante la acción de un campo eléctrico. Dicha propiedad puede cuantificarse mediante lo que se denomina permitividad relativa o constante dieléctrica del material. Figura 20. Efecto de un dieléctrico. Fuente de la imagen:[29] Todos los materiales dieléctricos se pueden considerar aislantes, pero se vuelven conductores cuando se sobrepasa el campo de ruptura del dieléctrico, es decir, cuando la tensión supera la tensión máxima o rigidez dieléctrica. Además, es importante remarcar que todos los dieléctricos son aislantes, pero no todos los aislantes son dieléctricos. La importancia de la presencia del dieléctrico es incuestionable. Como hemos comentado, el dieléctrico produce una oposición al paso de la corriente, y su principal objetivo es la de evitar la aparición de arcos eléctricos. En recientes investigaciones [12] se están estudiando los materiales utilizados como dieléctrico con el objetivo de minimizar la degradación de dichos materiales cuando son utilizados en actuadores de plasma Se presenta en la siguiente tabla los principales dieléctricos utilizados en los actuadores de plasma junto con los valores, tanto de la constante dieléctrica como de la tensión máxima o rigidez dieléctrica. Tipo de Dieléctrico Constante de rigidez Rigidez dieléctrica (KV/mm) PLA (poly-lactic acid) 3.2 50-64 PIB (polyisobutylene) rubber 2.82 40 Kapton 2.5-3.5 154-303 Acetoxy silicone 3.6 20 Tabla 3. Tipos de dieléctricos y propiedades. 1. Introducción 28 28 1.4.2.3 Ventajas y limitaciones Los estudios realizados durante las últimas décadas, de los cuales sólo hemos señalado algunos, han demostrado la posibilidad de utilizar éstos como estrategias de control activo. Se resumirán en esta sección las ventajas y limitaciones que suponen el uso de estos dispositivos de plasma [8]. Ventajas: • Posibilidad de ubicar estos actuadores en diferentes configuraciones y superficies, como cuerpos aerodinámicos. • Interferencias aerodinámicas insignificantes cuando los dispositivos están inactivos. • Tiempos de actuación altos (bajo delay). • Bajo consumo energético. • Bajo peso. • Fáciles de construir. • Baja complejidad mecánica • Posibilidad de ajustar la estrategia de activación/desactivación en función de la condición dinámica del fluido en particular. Limitaciones: • Velocidad inducida por debajo de los 10 m/s, lo que obliga a utilizar varios actuadores con el objetivo de satisfacer las necesidades aeronáuticas • Baja eficiencia de conversión eléctrica en cinética. • Uso de alta tensión potencialmente peligrosa y generación de ruido electromagnético. • Eficacia en el control de flujo para Reynolds por debajo de aproximadamente 3000-5000. 2. MODELADO DEL PROBLEMA Se exponen en este capítulo las ecuaciones que gobiernan la física del actuador de plasma DBD. Se desarrollarán, por tanto, las formulaciones matemáticas que envuelven el problema eléctrico y el problema fluido. Previamente a ello es necesario poner en relevancia el modelo de Suzen-Huang, que es el que utilizaremos para realizar el modelado del problema. 2.1 Modelo de Suzen-Huang En el caso de los actuadores de plasma [10], existen modelos de plasma autoconsistentes y modelos fenomenológicos para estimar la fuerza electrostática y la física del plasma subyacente. Durante un tiempo se trabajó una, dos o múltiples especies de plasma; trabajos más recientes durante la última década se han centrado en modelos de plasma de transporte de tres especies del plasma: electrones, iones positivos y negativos. La dinámica y la evolución de las especies de plasma con ciclos de voltaje positivo y negativo se han explorado con estos modelos de plasma autosuficientes, discutiendo las observaciones con las investigaciones experimentales sobre el flujo inducido debido al movimiento de las especies de plasma en el interior de un campo eléctrico externo. El desarrollo del Modelo de Suzen-Huang (SH) en 2005, permitió resolver los modelos del plasma de una forma más ingenieril y no mediante una física del plasma detallada, con el fin de explorar diferentes aplicaciones prácticas de los actuadores DBD con menor capacidad computacional que los modelos de plasma autoconsistentes completos. Es importante remarcar que el modelo SH no tiene en cuenta la dinámica de descarga, pues maneja modelos multiespecificicos más sofisticados. Su dependencia con la distribución experimental de carga hace que este modelo sea inadecuado para muchos investigadores que desean estudiar la dinámica de descarga asociada a los actuadores DBD; sin embargo, dicho modelo no pierde validez pues fue desarrollado principalmente para estudiar aplicaciones de control de flujo. Desde entonces, han sido muchas las publicaciones realizadas en todo el mundo estudiando los diferentes aspectos del modelo SH y sus diferentes aplicaciones ingenieriles: • Suzen et al. (2005) [30] demostró el modelo como un método fenomenológico efectivo al estudiar numéricamente la reducción de la separación del flujo a lo largo de una pala de turbina a baja presión. Luego mejoró el modelo para extender su aplicación a geometrías de actuador más compicadas. • Belson et al (2012) [31] presentaron un estudio de DNS sobre control de la transición del flujo. • Sato y cols.m [32] realizaron estudios LES de alta fidelidad sobre el control del flujo en un perfil NACA 0015 con un actuador de plasma DBD. • Un estudio paramétrico detallado para investigar los impactos de los dos parámetros escalares del modelo escalar de SH, la longitud de Debye y la densidad de carga máxima, fue realizado por Skote e Ibrahim. Ibrahim et al. (2012) [33] realizaron un estudio de control de flujo con un actuador de chorro sintético de plasma lineal (L-PSJA) y encontró que su modelo de SH modificado proporcionó unos mejores valores de velocidad que la forma estándar del modelo SH de 2005. • Skote et al. (2016) [34] publicaron un estudio con su modelo de SH modificado, que utiliza una L-PSJA para controlar los flujos de estela cilíndricos, que es un problema fundamental en la aerodinámica. • El rendimiento de los actuadores de plasma también se ha probado numéricamente con el modelo SH en flujos de canal: Mahfoze et al. (2017) [35] presentan simulaciones del modelo de SH de la reducción de la resistencia aerodinámica de flujos en canales, al alinear los actuadores de plasma cerca de las paredes de los canales para proporcionar chorros orientados en sentido transversal. • T. Brauner et al. (2016) [36] compararon los resultados de sus simulaciones en el actuador DBD realizados con los modelos SH y DNS de alto orden, con otro conjunto de resultados de simulación. 2. Modelado del problema 30 30 Así, la mayoría de las investigaciones validaron el modelo SH en sus respectivas aplicaciones, imponiendo al modelo SH como una herramienta potente, de bajo coste computacional y que evita el uso de modelos de sistemas complejos, alejados de las ecuaciones fluidas simples, como puede ser la formación de filamentación durante descargas. El funcionamiento de los actuadores DBD, se puede explicar con una breve descripción de la velocidad inducida del viento iónico. La diferencia de potencial entre los dos electrodos induce un fuerte campo eléctrico entre ambos, haciendo que las partículas cargadas se muevan bajo la influencia de dicho campo eléctrico. El impulso de las partículas cargadas se transfiere a las moléculas neutras del fluido, transmitiendo el movimiento según la fuerza de Coulomb inducida. Este movimiento fluido se denomina viento iónico. La zona de separación entre ambos electrodos (gap) es la región de máximo campo eléctrico, y como consecuencia, la región donde se produce la fuerza máxima de Coulomb. La componente vertical de esta fuerza hace que el fluido se dirija hacia la pared dieléctrica. Por otro lado, la componente horizontal de la fuerza empuja el fluido a lo largo de la pared dieléctrica, aguas abajo. El efecto combinado de ambas componentes produce un efecto de succción sobre la superficie de descarga, conectando de nuevo el flujo separado a la superficie del actuador. De esta forma, el modelo SH imita una descarga de plasma inducida por el actuador DBD. Con la distribución de carga obtenida mediante dicho modelo, se calcula la fuerza electrohidrodinámica EHD por unidad de volumen, que posteriormente se incorpora como término dominante o fuente en las ecuaciones de Navier-Stokes. 2.2 Formulación del problema eléctrico Una descripción matemática del modelo de Suzen-Huang es proporcionada en las referencias [30, 33, 44, 45]. Se presenta brevemente el modelo en este apartado. El objetivo es, como hemos comentado, la obtención de la fuerza electrohidrodinámica EHD por unidad de volumen. Dicha fuerza, despreciando fuerzas magnéticas, se puede expresar como: 𝐹𝑒 󰇍 󰇍 󰇍 =𝜌𝑐𝐸 󰇍 (1) donde 𝜌𝑐 representa la densidad de carga neta y 𝐸 󰇍 el campo eléctrico. Por tanto, el cálculo de 𝜌𝑐 y 𝐸 󰇍 es el siguiente paso. Por un lado, si las variaciones temporales del campo magnético resultan despreciables (hecho ocurre normalmente en el caso de plasmas), las ecuaciones de Maxwell determinan: ∇×𝐸 󰇍 ≈0 (2) lo que implica que el campo eléctrico deriva de un potencial, que se puede expresar como: 𝐸 󰇍 =−∇Φ (3) Según la ley de Gauss: ∇·(𝜀 𝐸 󰇍 )= 𝜌𝑐 (4) que también puede expresarse en función del potencial ∇·(𝜀 ∇Φ)=−𝜌𝑐 (5) donde 𝜀 es la permitividad, que viene dada por 𝜀=𝜀𝑟𝜀0, siendo 𝜀𝑟la permisividad relativa del medio y 𝜀0 es la permitividad del vacío. Por otro lado, puesto que las partículas de gas están débilmente ionizadas, podemos asumir que el potencial Φ se puede descomponer en dos partes: una debido al campo eléctrico externo, 𝜑 y la otra es debida a un potencial debido a la densidad de carga neta en el plasma, 𝜙. Por tanto: Φ=𝜙+𝜑 (6) Además, si asumimos que la longitude de Debye es pequeña y que la carga localmente en la pared no es grande, es posible considerar que: • La distribución de partículas cargadas en el dominio se rigen por el potential en la pared debido a la carga eléctrica en la pared, 𝜑. • El potencial asociado a campos eléctricos externos, 𝜙, tiene poca influencia en la distribución de dichas partículas cargadas. Por tanto, podemos escribir dos ecuaciones separadas en términos de estos dos potenciales, uno para el campo eléctrico externo debido al voltaje aplicado entre los dos electrodos: ∇·(𝜀𝑟∇𝜙)=0 (7) y otra que representa el potencial debido a las partículas cargadas ∇·(𝜀𝑟∇𝜑)=−𝜌𝑐 𝜀0 (8) El objetivo será por tanto, resolver ambas ecuaciones. La ecuación (7) se resuelve para obtener el potencial eléctrico 𝜙, usando el voltaje aplicado en los electrodos como condiciones de contorno. Las condiciones de contorno y el dominio computacional utilizados para resolver dicha ecuación pueden observarse en la figura 21: Figura 21. Dominio y condiciones de contorno para resolver el potencial 𝝓. Fuente de la imagen: [10] En dicha imagen puede observarse que la ecuación (7) se resuelve para los valores correctos de 𝜀𝑟. En el caso del dominio fluido se ha tomado el valor del aire 𝜀𝑟=1 y en el caso del dieléctrico se ha utilizado Kapton, con 𝜀𝑟=2.7. Por otro lado, en la región de separación dieléctrico-aire, la media entre 𝜀𝑟1 𝑦 𝜀𝑟2 debe ser usada con el objetivo de conservar el campo eléctrico. Además, es necesario también considerar el potencial eléctrico debido a la carga neta en el plasma, por lo que 2. Modelado del problema 32 32 tenemos que resolver la ecuación (8). La densidad de carga dentro del plasma en cualquier punto es: 𝜌𝑐 𝜀0=𝑒(𝑛𝑖−𝑛𝑒) 𝜀0=−𝑒𝑛0 𝜀0 [exp(−𝑒𝜑 𝑘𝑇𝑖)−exp(𝑒𝜑 𝑘𝑇𝑒)] (9) Donde 𝑛𝑖 es la densidad de iones positivos, 𝑛𝑒, la densidad de electrones y 𝑛0 la densidad de plasma neutro. T es la temperatura de las especies. Las densidades de las especies se obtienen mediante las siguientes ecuaciones: 𝑛𝑒=𝑛0·exp(𝑒𝜑 𝑘𝑇𝑒) (10) 𝑛𝑖=𝑛0·exp(−𝑒𝜑 𝑘𝑇𝑖) (11) Haciendo un desarrollo en series de Taylor para la función exponencial e introduciendo la definición de longitud de Debye: 1 𝜆𝐷=𝑒2𝑛0 𝜀0[1 𝑘𝑇𝑖+1 𝑘𝑇𝑒] (12) la ecuación (8) pasa a ser ∇·(𝜀𝑟∇𝜌𝑐)=𝜌𝑐 𝜆𝐷 2 (13) La ecuación (15) se resuelve para obtener la densdad neta de carga, 𝜌𝑐, sólo en el dominio del aire usando las condiciones de contorno, tal y como se observa en la figura 22. Figura 22. Dominio y condiciones de contorno que nos permiten calcular la densidad de carga 𝝆𝒄. Fuente de la imagen: [10] Como se puede observar, un gradiente normal e igual a cero es impuesto para la densidad de carga en las paredes excepto en la región que cubre el electrodo encapsulaso y en las paredes exteriores. Sobre la superficie de dieléctrico que cubre el electrodo encapsulado, se establece una distribución de densidad de carga que varía a lo largo del electrodo y con el tiempo, tal y como se expresa en la siguiente ecuación: 𝜌𝑐,𝑤(𝑥,𝑡)=𝜌𝑐𝑚𝑎𝑥 𝐺(𝑥) 𝑓(𝑡) (14) donde 𝜌𝑐𝑚𝑎𝑥 es el máximo valor de densidad de, G(x) es una distribución semi-Gaussiana que modela la distribución del plasma a lo largo del electrodo encapsulado y 𝑓(𝑡) representa la variación temporal del voltaje aplicado sobre los electrodos. Resultados experimentales sugieren que la distribución semi-Gaussiana es: 𝐺(𝑥)=𝑒[−(𝑥−𝜇)2 2𝜎2] (15) con 𝑥≥0, donde 𝜇 es el parámetro de localización que indica el máximo en x, y 𝜎 es el parámetro de escala que determina el ratio de decrecimiento. Tenemos, por tanto, todo lo necesario para calcular la fuerza electrohidrodinámica (EHD): 𝐹𝑒 󰇍 󰇍 󰇍 =𝜌𝑐𝐸 󰇍 =−𝜌𝑐∇𝜙 (16) Es importante remarcar que la ecuación (7) puede resolverse de forma adimensional. Igualmente, la ecuación (13) puede resolverse computacionalmente obteniendo la distribución de densidad de carga adimensional. Introduciendo las siguientes variables adimensionales: 𝑥=𝜉∗𝐿 ; 𝑧=𝜂∗𝐿 ; 𝜌𝑐=𝜌𝑐∗ 𝜌𝑐𝑀𝐴𝑋 𝑓(𝑡) ; 𝜙=𝜙∗𝜙𝑀𝐴𝑋 𝑓(𝑡) (17) donde 𝜙𝑚𝑎𝑥 es la amplitud del voltaje aplicado y 𝜌𝑐𝑀𝐴𝑋 es la máxima densidad de carga, que se ha obtenido experimentalmente. Por otro lado, 𝑓(𝑡)=𝑆𝑒𝑛(2𝜋𝑓𝑏∗𝑡) es la variación temporal del voltaje, que dependerá del tipo de onda que apliquemos, con 𝑓𝑏 la frecuencia del voltaje aplicado, 𝜌𝑐∗ es la distribución de densidad de carga adimensional y 𝜙∗ es la distribución de potencial eléctrico adimensional. Por último, L es un parámetro adimensional, que coincide con el valor del gap entre los electrodos. Definiendo el operador adimensional ∇= 1𝐿∇∗ las ecuaciones (7) y (13) se transforman en: ∇∗·(𝜀𝑟∇∗𝜙∗)=0 (18) ∇∗·(𝜀𝑟∇∗𝜌𝑐∗)=𝜌𝑐∗(𝐿 𝜆𝐷)2 (19) Por tanto, los problemas que hay que resolver son: Figura 23. Dominio y condiciones de contorno para resolver el potencial 𝝓∗ adimensional 3. Simulación del actuador de plasma DBD 40 40 Llamando g(x) al polinomio interpolador de Lagrange para la derivada y haciendo un proceso análogo: 𝑔(𝑥)=∑𝐿𝑗(𝑥)·𝑓𝑗′ 𝑁 𝑗=1 (45) 𝑓𝑖′′=𝑔(𝑥𝑖)′=∑𝐿𝑗′(𝑥𝑖)·𝑓𝑗′= ∑(∑𝐿𝑘 ′(𝑥𝑗)·𝑓𝑗) 𝑁 𝑘=1 𝐿𝑗′(𝑥𝑖) 𝑁 𝑗=1 𝑁 𝑗=1 (46) Así, la función y sus derivadas pueden escribirse de una forma más sencilla en forma vectorial y matricial como: 𝑓 =[𝑓1 𝑓2 … 𝑓𝑁], 𝑓 ′=[𝑓1′ 𝑓2′ … 𝑓3′], 𝑓 ′′=[𝑓1′′ 𝑓2′′ … 𝑓3′′] (47) De forma compacta: 𝑓 ′=𝐿𝑝𝑥𝑓 𝑓 ′′=𝐿𝑝𝑥𝑓 ′=𝐿𝑝𝑥 2𝑓 (48) Siendo 𝐿𝑝𝑥=(𝐿1′(𝑥1) ⋯ 𝐿𝑁 ′(𝑥1) ⋮ ⋱ ⋮ 𝐿1′(𝑥𝑁)⋯ 𝐿𝑁 ′(𝑥𝑁)) (49) Aplicando el mismo procedimiento, pueden obtenerse valores para aproximar las derivadas de cualquier orden n de f(x) como 𝑓 𝑛=𝐿𝑝𝑥 𝑛·𝑓 (50) Por otro lado, es importante elegir correctamente los N nodos x1, x2, …, xN en el interior del dominio [a,b], pudiendo ser éstos arbitrarios. Una buena elección de dichos nodos es crucial para que la función f(x) sea una buena aproximación en el intervalo [a , b] a la función incógnita. La mejor elección de los nodos se basa en utilizar los denominados nodos de Chebyshev, descubiertos por el matemático ruso Chebyshev, quien encontró que el polinomio de grado N debe ser elegido mediante la siguiente fórmula: 𝑇𝑁(𝑥)=2𝑥𝑇𝑁−1(𝑥)−𝑇𝑁−2(𝑥) (51) con 𝑇1(𝑥)=1 y 𝑇2(𝑥)=𝑥. Así, los puntos x1, x2, …, xN deben coincidir con las raíces del polinomio 𝑇𝑁(𝑥) y cumplen: 𝑥𝑖=𝑎+𝑏 2+𝑎−𝑏 2·cos[𝜋 (𝑖−1) 𝑁−1] 𝑖=1…𝑁 (52) Resulta interesante, además, comentar algunos aspectos acerca de los nodos de Chebyshev. Por un lado, los nodos de Chebyshev son aquellos que permiten obtener la mejor aproximación mediante polinomios de la solución. Tal y como se puede observar en la siguiente figura, la aproximación obtenida es mucho mejor que mediante nodos equiespaciados Figura 26. Comparación aproximación polinómica mediante nodos equiespaciados y mediante nodos de Chebyshev. Fuente de la imagen: [39] Por otro lado, los nodos de Chebyshev no distan la misma distancia unos de otros, concentrándose en los extremos del intervalo. Esto es interesante, por ejemplo, cuando se estudia un fluido en las inmediaciones de una pared. 3.1.2 El método de colocación 2D: nodos de Chebyshev Extenderemos en esta sección lo expuesto anteriormente para el caso de problemas bidimensionales. Teniendo en cuena una función cuyas variables naturales son x y z, Φ=Φ (𝑥,𝑧), definida en un dominio con forma rectangular definido por 0≤𝑥≤𝐿 y 0≤𝑧≤𝐻 y con un total de Nx · Nz nodos, cumpliéndose { (xm, zn) : m=1, … Nx , n = 1, … Nz } y que para cualquier valor de x se puede realizar una interpolación de Lagrange para la dependencia con z, de la forma: Φ(𝑥,𝑧)=∑𝐿𝑧,𝑛(𝑧)·Φ(𝑥,𝑧𝑛) 𝑁𝑧 𝑛=1 (53) Haciendo uso de la fórmula de interpolación de lagrange para cada una de las funciones de la coordenada x, Φ(𝑥,𝑧𝑛), llegamos a: Φ(𝑥,𝑧)=∑∑[ 𝑁𝑧 𝑛=1 𝐿𝑥,𝑚(𝑥)𝐿𝑧,𝑛(𝑧)Φ(𝑥𝑚,𝑧𝑛) 𝑁𝑥 𝑚=1 ] (54) siendo Φ(𝑥𝑚,𝑧𝑛) el valor que toma la función en el nodo (xm,zn) y 𝐿𝑧,𝑛, 𝐿𝑥,𝑚 los polinomios interpolantes de Lagrange para las variables x y z, 3. Simulación del actuador de plasma DBD 42 42 𝐿𝑥,𝑚(𝑥)=∏𝑥−𝑥𝑗 𝑥𝑚−𝑥𝑗 𝐿𝑧,𝑛(𝑧)=∏𝑧−𝑧𝑗 𝑧𝑛−𝑧𝑗 𝑁𝑧 𝑗=1 𝑗≠𝑛 𝑁𝑥 𝑗=1 𝑗≠𝑚 (55) Análogamente a como se hizo en problemas unidimensionales, se calculan las derivadas parciales de la función Φ(𝑥,𝑧) como: Φ𝑥(𝑥,𝑧)≡𝜕Φ(𝑥,𝑧) 𝜕𝑥 =∑∑𝐿𝑥,𝑚 ′(𝑥)·𝐿𝑧,𝑛(𝑧)·Φ(𝑥𝑚,𝑧𝑛) 𝑁𝑧 𝑛=1 𝑁𝑥 𝑚=1 (56) Φ𝑧(𝑥,𝑧)≡𝜕Φ(𝑥,𝑧) 𝜕𝑧 =∑∑𝐿𝑥,𝑚(𝑥)·𝐿𝑧,𝑛 ′(𝑧)·Φ(𝑥𝑚,𝑧𝑛) 𝑁𝑧 𝑛=1 𝑁𝑥 𝑚=1 (57) Además, elegimos nodos de Chebyshev tanto en x como en z, por lo que los puntos en el interior de nuestro dominio serán de la forma: 𝑥𝑖=𝑥min +𝑥𝑚𝑎𝑥 2+𝑥𝑚𝑖𝑛−𝑥𝑚𝑎𝑥 2·cos[𝜋 (𝑖−1) 𝑁𝑥−1] 𝑑𝑜𝑛𝑑𝑒 𝑖=1…𝑁𝑥 (58) 𝑧𝑗=𝑧𝑚𝑖𝑛+𝑧𝑚𝑎𝑥 2+𝑧𝑚𝑖𝑛−𝑧𝑚𝑎𝑥 2·cos[𝜋 (𝑗−1) 𝑁𝑧−1] 𝑑𝑜𝑛𝑑𝑒 𝑗=1…𝑁𝑧 (59) Una vez que se realiza el mallado del dominio en el conjunto de Nx · Nz puntos (𝑥𝑖,𝑧𝑗) tal y como se muestra en la siguiente figura Figura 27. Mallado del dominio de forma rectangular. Fuente de la imagen: [39] es conveniente, al igual que se hace en el método de los elementos finitos, hacer una asignación biunívoca mediante un único índice, de la forma: (𝑥𝑖,𝑧𝑗)↔𝐼=(𝑖−1)·𝑁𝑧+𝑗 (60) Obsérvese que mediante el úso de un único índice se pueden considerar los productos entre los interpolantes de Lagrange y sus derivadas como elementos de matrices Dx y Dz definidas como: 𝐷𝑥(𝐼,𝐾)=𝐿𝑥,𝑚 ′(𝑥𝑖)𝐿𝑧,𝑛(𝑧𝑗) 𝑦 𝐷𝑧(𝐼,𝐾)=𝐿𝑥.𝑚(𝑥𝑖)𝐿𝑧,𝑛 ′(𝑧𝑗) 𝐼,𝐾=1,…,𝑁𝑡 (61) con 𝐼=(𝑖−1)·𝑁𝑧+𝑗 , 𝐾=(𝑚−1)·𝑁𝑧+𝑛 y 𝑁𝑡=𝑁𝑥·𝑁𝑧 el número total de puntos. Simplificando, las expresiones de las derivadas parciales de Φ(𝑥,𝑧): 𝜕Φ 𝜕𝑥|𝐼=∑𝐷𝑥(𝐼,𝐾)·Φ𝐾 𝑁𝑡 𝐾=1 (62) 𝜕Φ 𝜕𝑧|𝐼=∑𝐷𝑧(𝐼,𝐾)·Φ𝐾 𝑁𝑡 𝐾=1 (63) donde se han definidos los valores de la función en los nodos como Φ𝐾=Φ(𝑥𝑚,𝑧𝑛). Utilizando notación vectorial y matricial Φ 󰇍 󰇍 󰇍 ≡[Φ1 Φ2 … Φ𝑁𝑡]; Φ 󰇍 󰇍 󰇍 𝑥≡ [ Φ𝑥1 ′ Φ𝑥2 ′ … Φ𝑥𝑁𝑡 ′ ] ; Φ 󰇍 󰇍 󰇍 𝑧≡ [ Φ𝑧1 ′ Φ𝑧2 ′ … Φ𝑧𝑁𝑡 ′ ] (64) y haciendo uso de las matrices Dx y Dz : Φ𝑥 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑥·Φ 󰇍 󰇍 󰇍 Φ𝑧 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑧·Φ 󰇍 󰇍 󰇍 (65) Análogamente al problema unidimensional, las derivadas sucesivas se pueden calcular como Φ𝑥𝑥 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑥·Φ𝑥 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑥2Φ 󰇍 󰇍 󰇍 Φ𝑧𝑧 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑧·Φ𝑧 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑧2Φ 󰇍 󰇍 󰇍 Φ𝑥𝑧 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑥·𝐷𝑧·Φ 󰇍 󰇍 󰇍 =𝐷𝑧·𝐷𝑥·Φ 󰇍 󰇍 󰇍 =Φ𝑧𝑥 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 (66) 3.1.3 El Método de colocación 2D: diferencias finitas Análogamente al caso del Método de colocación 2D: nodos de Chebyshev, es posible también calcular las matrices de derivadas para el caso denominado como ‘’diferencias finitas’’, haciendo uso de los polinomios de Lagrange según se ha implementado en las funciones auxiliares que aparecen en el Anexo 3: Funciones auxiliares. Se entenderá por el término ‘’diferencias finitas’’ el considerar un entorno en cada punto de la malla que, distinto a Chebyshev, contiene un número de puntos menor que el de la malla, en nuestro caso 7, lo que permite hallar derivadas de cuarto orden con precisión de ℎ3 y matrices mucho más dispersas, lo que es especialmente conveniente cuando la malla en alguna dirección es grande. En el caso de requerir más información, acudir a [46]. 3. Simulación del actuador de plasma DBD 44 44 3.2. Resolución numérica del problema eléctrico Se expone en este apartado la resolución numérica de la distribución de densidad de carga y de la distribución de potencial eléctrico. Un desarrollo más extenso se lleva a cabo en el Trabajo de Fin de Grado de M. A. Domínguez [43]. Presentamos aquí un breve resumen. 3.2.1 Aplicación del método de colocación: Ecuación de Poisson La ecuación que controla la distribución de densidad de carga en el interior del dominio fluido y la distribución de potencial eléctrico en todo el actuador DBD es una ecuación de tipo Poisson. Previo a su resolución será necesario, por tanto, resolver una ecuación de Poisson de carácter general, hecho que se realiza a continuación. Consideremos la ecuación de Poisson: 𝜕2Φ 𝜕𝑥2+𝜕2Φ 𝜕𝑧2=𝑓(𝑥,𝑧) (67) y un dominio de longitud L y altura H (rectangular) donde el término fuente 𝑓(𝑥,𝑧) es conocido. Supongamos que las condiciones de contorno que deben cumplirse en este dominio son las siguientes: Φ(0,𝑧)=0, Φ(𝑥,0)=0, Φ(𝐿,𝑧)=0, Φ(𝑥,𝐻)=0 (68) Adimensionalizando: 𝑥=𝑥𝐻, 𝑧=𝑧𝐻, Φ =Φ Φ0 (69) siendo Φ0 un valor de referencia de Φ y sustituyendo, llegamos a la siguiente ecuación: 𝜕2Φ  𝜕𝑥2+𝜕2Φ  𝜕𝑧2=𝑓󰆹(𝑥,𝑧) (70) donde 𝑓󰆹=𝑓𝐻2 Φ0. Por otro lado, las condiciones de contorno serán ahora: Φ (0,𝑧)=0, Φ (𝑥,0)=0, Φ (𝐿𝐻,𝑧)=0, Φ (𝑥,1)=0 (71) Haciendo uso de una malla de Nt puntos siendo Nt=Nx x Nz y aplicando lo expuesto anteriormente en el apartado 3.1, se llega a: (𝐷𝑥2+𝐷𝑧2)Φ 󰇍 󰇍 󰇍 ≡𝐷𝐿·Φ 󰇍 󰇍 󰇍 =𝑓 (72) donde los vectores Φ 󰇍 󰇍 󰇍 y 𝑓 contienen los valores de las funciones Φ (𝑥,𝑧) y 𝑓󰆹(𝑥,𝑧) en los diferentes nodos de la malla I=(i-1) Nz + j. Por último, antes de resolver el sistema algebraico al que hemos llegado, se deben introducir las condiciones de contorno del problema. Para ello, se modifica la matriz del sistema 𝐷𝐿 y el vector 𝑓 , de forma que si el valor en un nodo I coincide con el de un nodo situado en algún lado de la cavidad, entonces todos los elementos de la fila I de la matriz 𝐷𝐿 se anulan excepto el correspondiente a la diagonal, que toma el valor unidad, 𝐷𝐿(𝐼,𝐼)=1 y se anula también la componente 𝑓 (𝐼) puesto que en este caso se ha impuesto que Φ  se anule en todo el contorno. 3.2.2 Resolución numérica de la distribución de densidad de carga En este apartado se va a implementar un programa de MATLAB para la resolución numérica, utilizando el método de colocación explicado anteriormente, de la ecuación (19) y condiciones de contorno que se tienen en la figura 24, que gobiernan la distribución de densidad de carga adimensional en el actuador DBD, cuyo dominio es un dominio rectangular de altura H y longitud L. En nuestro caso, la distribución de densidad de carga adimensional se obtiene en el dominio fluido, pues en el dominio del dieléctrico ésta es nula. Por este motivo, se utilizará en lo que sigue el subíndice a (más adelante se utilizará el subíndice b para el dieléctrico). La discretización del problema se lleva a cabo introduciendo un mallado mediante nodos de Chebyshev en la dirección z y mediante diferencias finitas en la dirección x: 𝜂𝑗=𝜂𝑚𝑖𝑛,𝑎+𝜂𝑚𝑎𝑥,𝑎 2+𝜂𝑚𝑖𝑛,𝑎−𝜂𝑚𝑎𝑥,𝑎 2·cos[𝜋 (𝑗−1) 𝑁𝑧−1] 𝑑𝑜𝑛𝑑𝑒 𝑗=1…𝑁𝑧 𝜉𝑖=𝜉min+𝑖−1 𝑁−1(𝜉max−𝜉min ) 𝑑𝑜𝑛𝑑𝑒 𝑖=1…𝑁 (73) donde 𝜉𝑚𝑖𝑛 ,𝜉𝑚𝑎𝑥 ,𝜂𝑚𝑖𝑛,𝑎 𝑦 𝜂𝑚𝑎𝑥,𝑎 son valores conocidos. En nuestro caso, N no coincide exactamente con 𝑁𝑥 , pues tenemos dividido el dominio en 𝜉 en tres intervalos, cada uno con un número total de puntos igual a 𝑁1 ,𝑁2 𝑦 𝑁3, respectivamente, cumpliéndose que 𝑁𝑥=𝑁1+ 𝑁2+ 𝑁3. Este espaciado en 𝜉 será común tanto para el cálculo de la densidad de carga como del potencial eléctrico y el campo de velocidades. Tal y como se explicó en el apartado del método de colocación será necesario enumerar los nodos mediante índices simples, lo que puede hacerse mediante la expresión: (𝜉𝑖,𝜂𝑗𝑎)↔𝐼𝑎=(𝑖−1)·𝑁𝑧𝑎+𝑗 (74) lo que permite expresar el producto entre polinomios interpolantes y sus derivadas como elementos de matrices, y con ello poder escribir el problema en forma vectorial. Una vez se ha definido la malla para la discretización especial, se explica cómo se aplica el método de colocación visto en el apartado 3.1 y su posterior aplicación en el apartado 3.2.1. Como sabemos, la expresión que gobierna la distribución de densidad de carga es la siguiente: ∇∗·(𝜀𝑟∇∗𝜌𝑐∗)=𝜌𝑐∗(λD 𝐿)2 (75) Ésta se puede reescribir como: 𝜀𝑟 ∇∗2ρc∗=𝜀𝑟(𝜕2𝜌𝑐∗ 𝜕𝜉2+𝜕2𝜌𝑐∗ 𝜕𝜂2)=𝜌𝑐∗(λD 𝐿)2 (76) Para discretizar la ecuación de Poisson, definimos la distribución de densidad de carga adimensional en los 𝑁𝑡𝑎=𝑁𝑥𝑎∗𝑁𝑧𝑎 nodos como: 𝜌𝑐∗(𝜉,𝜂𝑎)≡𝜌𝑐∗(𝐼)≡𝜌𝑐∗ 󰇍 󰇍 󰇍 󰇍 (77) donde 𝜌𝑐∗ 󰇍 󰇍 󰇍 󰇍 designa, en la notación usada en el Código de MATLAB, el vector de dimension 𝑁𝑡𝑎 𝑥 1, cuyas componentes I albergarán los valores de ρc(𝜉,𝜂𝑎) en el nodo I. Teniendo en cuenta lo expuesto en la aplicación del método de colocación para la ecuación de Poisson, la discretización correspondiente al operador laplaciano adimensional es: 𝐷𝐿𝑎 ∗=(𝐷𝜉2+𝐷𝜂𝑎 2) (78) 3. Simulación del actuador de plasma DBD 46 46 Por tanto, tenemos la siguiente expresión: ∇∗2ρc∗|𝐼=∑𝐷𝐿𝑎 ∗(𝐼𝑎,𝐾𝑎) ρc∗(Ka)≡(DLa 𝜌𝑐∗ 󰇍 󰇍 󰇍 󰇍 )𝐼𝑎 𝑁𝑡 𝐾𝑎=1 (79) Así, llegamos a: 𝜀𝑟DL∗𝑎 𝜌𝑐∗ 󰇍 󰇍 󰇍 󰇍 =𝜌𝑐∗ 󰇍 󰇍 󰇍 󰇍 (λD 𝐿)2 (80) que se puede reescribir como: (𝜀𝑟DLa ∗−𝐸𝑦𝑒(𝑁𝑡𝑎)(λD 𝐿)2)𝜌𝑐∗ 󰇍 󰇍 󰇍 󰇍 =0 (81) donde 𝐸𝑦𝑒(𝑁𝑡𝑎) es la matriz identidad de orden 𝑁𝑡𝑎. Por tanto, el sistema a resolver es: 𝐴 𝜌𝑐∗ 󰇍 󰇍 󰇍 󰇍 =𝑏 󰇍 (82) donde obviamente, 𝐴=(𝜀𝑟DL∗a−𝐸𝑦𝑒(𝑁𝑡𝑎)(λD 𝐿)2) y 𝑏 󰇍 =0 󰇍 . El último paso será, antes de resolver el sistema de ecuaciones algebraico que nos permitirá obtener la solución, introducir las condiciones de contorno de la figura 24 tanto para la matriz del sistema A como en el término fuente 𝑏 󰇍 . Las modificaciones que deben realizarse para poder implementar las condiciones de contorno son: • Si el índice 𝐼𝑎 coincide con el de un nodo perteneciente a una de las paredes verticales o a la pared horizontal superior del dominio, todos los elementos de la fila 𝐼𝑎 en la matriz A se anulan, excepto el componente correspondiente a la diagonal, que toma 1. También se anula la componente 𝐼𝑎 del vector 𝑏 debido a que se impone que la densidad de carga se anule en los extremos del dominio. • Si el índice 𝐼𝑎 coincide con el de un nodo perteneciente a la pared horizontal inferior, que se corresponde con la entrefase dieléctrico-fluido, se debe cumplir que la derivada de la densidad de carga adimensional respecto a la dirección normal a la pared sea nula en toda la superficie excepto en la superficie correspondiente a la proyección del electrodo embebido, donde se impone que la densidad de carga adimensional se distribuye como una distribución semi-gaussiana. Así, se impone que toda la fila 𝐼𝑎 de la matriz A sea igual que la fila 𝐼𝑎 de la matriz 𝐷𝜂𝑎 y que el término fuente en la componente 𝐼𝑎 sea igual al valor de la derivada, en este caso 0. • Por último, si 𝐼𝑎 coincide con el de un nodo perteneciente a la superficie correspondiente al electrodo embebido, se impone que la fila I de la matriz A sea nula excepto en la diagonal, donde valdrá 1. La fila 𝐼𝑎 del término fuente valdrá 𝐺(𝜉). 3.2.3 Resolución numérica de la distribución de potencial eléctrico En este apartado se va a implementar un programa de MATLAB para la resolución numérica, utilizando el método de colocación explicado anteriormente, de la ecuación (18) y condiciones de contorno que se tienen en la figura 23, que gobiernan la distribución de potencial eléctrico adimensional en el actuador DBD. A diferencia de la distribución adimensional de densidad de carga que únicamente se calculaba en el dominio fluido, el potencial eléctrico deberá calcularse tanto en la región fluida como en la dieléctrica, pues presentan diferente permitividad relativa. En nuestro caso, el dominio será un dominio rectangular que acopla la región fluida, comprendida entre 𝜉𝑚𝑖𝑛 ,𝜉𝑚𝑎𝑥 ,𝜂min,𝑎 𝑦 𝜂𝑚𝑎𝑥,𝑎 y la región dieléctrica, delimitada por 𝜉𝑚𝑖𝑛 ,𝜉𝑚𝑎𝑥 ,𝜂min,𝑏 𝑦 𝜂𝑚𝑎𝑥,𝑏 Al igual que antes, la discretización del problema se lleva a cabo introduciendo un mallado mediante nodos de Chebyshev en la dirección z y mediante diferencias finitas en la dirección x tanto en la región fluida (subíndice a) como en la región dieléctrica (subíndice b): 𝜂𝑗𝑎,𝑏=𝜂min,(𝑎,𝑏) +𝜂𝑚𝑎𝑥,(𝑎,𝑏) 2+𝜂min,(𝑎,𝑏) −𝜂𝑚𝑎𝑥,(𝑎,𝑏) 2·cos[𝜋 (𝑗−1) 𝑁𝑧−1] 𝑑𝑜𝑛𝑑𝑒 𝑗=1…𝑁𝑧𝑎,𝑏 𝜉𝑖=𝜉min+𝑖−1 𝑁−1(𝜉max−𝜉min ) 𝑑𝑜𝑛𝑑𝑒 𝑖=1…𝑁 (83) donde 𝜉𝑚𝑖𝑛 ,𝜉𝑚𝑎𝑥 ,𝜂min,(𝑎,𝑏) 𝑦 𝜂𝑚𝑎𝑥,(𝑎,𝑏) son, de nuevo, valores conocidos. Tal y como se explicó en el apartado del método de colocación será necesario enumerar los nodos mediante la relación biunívoca, que se hace mediante la expresión: (𝜉𝑖𝑎,𝑏,𝜂𝑗𝑎,𝑏)↔𝐼𝑎,𝑏=(𝑖−1)·𝑁𝑧𝑎,𝑏+𝑗 (84) lo que permite escribir el problema en forma vectorial. Una vez se ha definido la malla para la discretización especial, se explica cómo se aplica el método de colocación. La expresión que gobierna la distribución adimensional de potencial eléctrico es, como sabemos, la siguiente: ∇∗·(𝜀𝑟∇∗𝜙𝑎,𝑏 ∗)=0 (85) Dicha expresión, se puede reescribir como: ∇∗2(𝜙∗)=𝜕2𝜙𝑎,𝑏 ∗ 𝜕𝜉2+𝜕2𝜙𝑎,𝑏 ∗ 𝜕𝜂2=0 (86) Para discretizar dicha ecuación, definimos la distribución de potencial eléctrico adimensional en los 𝑁𝑡𝑎,𝑏= 𝑁𝑥∗𝑁𝑧𝑎,𝑏 nodos como: 𝜙𝑎,𝑏 ∗(𝜉,𝜂𝑎,𝑏)≡𝜙∗(𝐼𝑎,𝑏)≡𝜙𝑎,𝑏 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 (87) donde 𝜙𝑎,𝑏 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 designa, en la notación usada en el Código de MATLAB, el vector de dimension 𝑁𝑡𝑎,𝑏 𝑥 1, cuyas componentes I albergarán los valores de 𝜙𝑎,𝑏 ∗(𝜉,𝜂𝑎,𝑏). De nuevo, la discretización correspondiente al operador laplaciano adimensional es: 𝐷𝐿𝑎,𝑏 ∗=(𝐷𝜉2+𝐷𝜂𝑎,𝑏 2) (88) Por lo que ∇∗2𝜙𝑎,𝑏 ∗|𝐼= ∑ 𝐷𝐿∗(𝐼𝑎,𝑏,𝐾𝑎,𝑏) 𝜙𝑎,𝑏 ∗(Ka,b)≡(𝐷𝐿𝑎,𝑏 ∗ 𝜙𝑎,𝑏 ∗ 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 )𝐼𝑎,𝑏 𝑁𝑡𝑎,𝑏 𝐾𝑎,𝑏=1 (89) Así, llegamos a 𝐷𝐿𝑎,𝑏 ∗ 𝜙𝑎,𝑏 ∗ 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 =0 (90) Por tanto, el sistema a resolver es 3. Simulación del actuador de plasma DBD 48 48 𝐴𝑎,𝑏 𝜙𝑎,𝑏 ∗ 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 =𝑏𝑎,𝑏 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 (91) donde 𝐴𝑎,𝑏=𝐷𝐿𝑎,𝑏 ∗ y 𝑏𝑎,𝑏 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 =0 󰇍 De esta forma, la resolución del problema podrá realizarse por dos vías: • Imponer las condiciones de contorno en el dominio fluido, en el dominio del dieléctrico y en la entrefase, resolviendo ambos sistemas de ecuaciones por separado. 𝐴𝑎 𝜙𝑎∗ 󰇍 󰇍 󰇍 󰇍 󰇍 =𝑏𝑎 󰇍 󰇍 󰇍 󰇍 𝐴𝑏 𝜙𝑏∗ 󰇍 󰇍 󰇍 󰇍 󰇍 =𝑏𝑏 󰇍 󰇍 󰇍 󰇍 (92) • Imponer las condiciones de contorno en el dominio fluido, en el dominio del dieléctrico y en la entrefase, acoplando las dos regiones para resolver el problema completo con un solo sistema, mediante matrices que a su vez tengan en cuenta los aspectos asociados a la región fluida (𝐴𝑎), aquellos asociados a la región dieléctrica (𝐴𝑏) y además la interacción entre ambos medios a través de la entrefase común (𝐵𝑎 y 𝐵𝑏) ,esto es, un sistema del tipo: [𝐴𝑎𝐵𝑎 𝐵𝑏𝐴𝑏][𝜙𝑎∗ 󰇍 󰇍 󰇍 󰇍 󰇍 𝜙𝑏∗ 󰇍 󰇍 󰇍 󰇍 󰇍 ]=[𝑏𝑎 󰇍 󰇍 󰇍 󰇍 𝑏𝑏 󰇍 󰇍 󰇍 󰇍 ] (93) En nuestro caso, la última opción ha sido la desarrollada en este estudio. Por último, será necesario imponer las condiciones de contorno de la figura 23. Las modificaciones que deben realizarse para poder implementar las condiciones de contorno son: • Si 𝐼𝑎,𝑏 coincide con un nodo perteneciente a las paredes horizontales o verticales excepto en la entrefase entre ambos dominios, se debe cumplir que la derivada del potencial eléctrico adimensional respecto a la dirección normal a la pared sea nulo en toda la superficie. Para ello la fila 𝐼𝑎,𝑏 de las matrices 𝐴𝑎,𝑏 se impone que sea igual a la fila I de las matrices 𝐷𝜉𝑎,𝑏(en las paredes verticales) y 𝐷𝜉𝑎,𝑏(en las paredes horizontales), mientras que la componente 𝐼𝑎,𝑏 del término fuente 𝑏𝑎,𝑏 se impone que sea igual al valor de la derivada, en este caso 0. • Si 𝐼𝑎,𝑏 coincide con un nodo perteneciente al electrodo expuesto, se anula la fila 𝐼𝑎,𝑏 de la matriz 𝐴𝑎,𝑏 excepto en la diagonal, que vale 1. Por otro lado, la componente 𝐼𝑎,𝑏 del término fuente 𝑏𝑎,𝑏 se hace idénticamente 1. • Si 𝐼𝑎,𝑏 coincide con un nodo perteneciente al electrodo embebido, se anula la fila 𝐼𝑎,𝑏 de la matriz 𝐴𝑎,𝑏 excepto en la diagonal, que vale 1. Por otro lado, la componente 𝐼𝑎,𝑏 del término fuente 𝑏𝑎,𝑏 se hace 0. • Por ultimo, en la entrefase será necesario imponer dos condiciones: la continuidad del campo eléctrico en la dirección 𝜉 y la continuidad en la dirección 𝜂 del producto 𝜀𝑎,𝑏𝐸𝜂 siendo 𝜀𝑎,𝑏 la permitividad de cada uno de los dominios. Para ello, se hace lo siguiente: la fila 𝐼𝑎 de la matriz 𝐴𝑎 se hace igual a la fila 𝐼𝑎 de la matriz 𝐷𝜉𝑎, la fila 𝐼𝑏 de la matriz 𝐵𝑎 se hace igual a la fila 𝐼𝑏 de la matriz −𝐷𝜂𝑎y por último, se anula la componente 𝐼𝑎 del término 𝑏𝑎. Posteriormente, se hace lo siguiente: la fila 𝐼𝑏 de la matriz 𝐴𝑏 se hace igual a la fila 𝐼𝑎 de la matriz 𝜀𝑏𝐷𝜂𝑏, la fila 𝐼𝑏 de la matriz 𝐵𝑏 se hace igual a la fila 𝐼𝑏 de la matriz −𝜀𝑎𝐷𝜂𝑎y por último, se anula la componente 𝐼𝑏 del término 𝑏𝑏. 3.2.4 Programación del problema eléctrico en MATLAB La lógica empleada en la simulación del problema eléctrico con el objetivo de calcular la distribución de densidad de carga, la distribución de potencial eléctrico y la fuerza electrohidrodinámica es la siguiente: Figura 28. Lógica empleada en MATLAB para resolver el problema eléctrico Puesto que este estudio se ha centrado en los aspectos hidrodinámicos del Actuador DBD, se adjuntan en el Anexo de programas de MATLAB los programas utilizados para el cálculo de la distribución de densidad de carga, la distribución de potencial eléctrico y la fuerza electrohidrodinámica. En la sección 3.3.3 se explicará de una forma exhaustiva la lógica empleada en la simulación del problema hidrodinámico, con el fin de dejar constancia de la transparencia del Código utilizado, a diferencia de los softwares comerciales utilizados en las diferentes bibliografías aportadas. Inicio del programa Definición de parámetros Mallado del dominio Definición del sistema de ecuaciones Condiciones de contorno Resolución del sistema Representación de resultados Fin del programa 3. Simulación del actuador de plasma DBD 56 56 • Mallado del dominio: se realiza el mallado del dominio que nos permite discretizar las ecuaciones. Éste se realiza tanto mediante nodos de Chebyshev como con diferencias finitas. • Cálculo de matrices: en este bloque se realiza el cálculo de las matrices 𝐷𝑥, 𝐷𝑧 ,𝐷𝐿,𝐷𝐿2, 𝐷𝑥𝐿, 𝐷𝑧𝐿 necesarias para aplicar la discretización basada en el método de colocación. Este bloque hace uso de funciones auxiliares que se muestran en el Anexo: Programas de MATLAB • Carga de datos del problema eléctrico: en este bloque se ejecuta el problema eléctrico, de forma que se obtienen tanto la distribución de densidad de carga como el campo eléctrico. Éstos, como ya se ha comentado con anterioridad, son introducidos en las ecuaciones de Navier-Stokes para obtener el campo de velocidades. • Definición del sistema de ecuaciones y condiciones de contorno: En este bloque se definen las matrices y vectores necesarios para la resolución del sistema y se imponen las condiciones de contorno que deben cumplirse para satisfacer el problema de acuerdo con lo expuesto en el apartado 3.3.2 3. Simulación del actuador de plasma DBD 58 58 • Resolución del sistema y representación de resultados. Mediante una simulación a lo largo del tiempo se resuelve el sistema de ecuaciones descrito y se grafican los resultados que se deseen. En concreto, se calcula la distribución de la función de corriente en el dominio; posteriormente se calcula el campo de velocidades y, finalmente, se grafican a lo largo del tiempo. • Fin del programa 3. Simulación del actuador de plasma DBD 60 60 4 ANÁLISIS DE RESULTADOS El código numérico presentado en la sección 3.3.3 nos permite obtener, mediante pequeñas modificaciones, resultados que pueden ser comparados con otras investigaciones. En primer lugar, se llevará a cabo la validación del código correspondiente al problema hidrodinámico, que se realiza en la sección 4.1, justificando los resultados obtenidos con los de la literatura en [10]. Posteriormente, en la sección 4.2, se estudia el comportamiento del actuador DBD frente a la variación de parámetros como la densidad de carga máxima o la longitud de Debye, que no son estudiados en [10], pero que tienen una relevancia importante en estudios como el de [42]. En la sección 4.3, se analiza el campo de velocidades del viento iónico producido mediante una simulación temporal, de acuerdo con [18]. Por último, en la sección 4.4, se estudia el coste computacional necesario para llevar a cabo las simulaciones. 4.1 Validación del código del problema hidrodinámico 4.1.1 Geometría del actuador de plasma DBD Se presenta en la figura 30 la configuración geométrica del actuador de plasma utilizada en la validación del código fluidodinámico con la literatura en [10] Figura 30. Geometría del actuador de plasma DBD • Espesor de dieléctrico: 3 𝑚𝑚 • Longitud del electrodo expuesto: 10 𝑚𝑚 • Longitud del electrodo embebido: 20 𝑚𝑚 • Distancia de separación (gap) entre electrodos: 1 𝑚𝑚 • Dominio fluido: 90 𝑚𝑚 𝑥 20 𝑚𝑚 • Dominio dieléctrico: 90 𝑚𝑚 𝑥 3 𝑚𝑚 Análisis de resultados 62 62 4.1.2 Parámetros Los parámetros utilizados en las siguientes simulaciones son, de acuerdo con [10]: • Amplitud del voltaje AC aplicado: 𝜙𝑀𝐴𝑋=20 𝑘𝑉 • Densidad de carga máxima: 𝜌𝑐,𝑀𝐴𝑋=0.001𝐶 𝑚3 • Frecuencia de la señal AC: 𝑓𝑏=1 𝑘𝐻𝑧 • Densidad del aire: 𝜌=1.225𝑘𝑔 𝑚3 • Viscosidad dinámica del aire: 𝜇=1.85∗10−5𝑘𝑔 𝑚 𝑠 • Longitud de Debye: 𝜆𝐷=0.001 𝑚 • Longitud de referencia: 𝐿=0.001 𝑚 • Velocidad de referencia: 𝑢0=√𝜌𝑐,𝑀𝐴𝑋 𝜌 𝜙𝑀𝐴𝑋=4.04𝑚 𝑠 • Número de Reynolds: 𝑅𝑒=𝜌 𝑢0 𝐿 𝜇=267.65 • Número adimensional correspondiente a la fuerza de Coulomb: 𝐷𝑐=𝜌𝑐,𝑀𝐴𝑋 𝜙𝑀𝐴𝑋 𝜌 𝑢02=1 En los siguientes apartados se estudiarán el efecto de diferentes parámetros en la velocidad máxima inducida del viento iónico, manteniéndose el resto constantes. 4.1.3 Efecto del espesor de dieléctrico Como se demuestra en estudios experimentales como los de [19], el espesor de dieléctrico influye sustancialmente tanto en la fuerza inducida como en la velocidad inducida del viento iónico. La evolución a lo largo del tiempo de la velocidad inducida máxima del viento iónico en función del espesor de dieléctrico queda recogida en la siguiente gráfica: Figura 31. Efecto del espesor de dieléctrico en la velocidad máxima obtenido con MATLAB Este resultado es acorde con la optimización llevada a cabo en [43], donde el espesor óptimo de dieléctrico es 1 mm, lo cual hace que la fuerza media generada sea máxima. Se puede observar cómo, a excepción del caso en el que el espesor del dieléctrico es 1 mm, la gráfica se adecua bastante bien a los resultados obtenidos en [10]: Figura 32. Efecto del espesor de dieléctrico en la velocidad máxima obtenido con Oracle3D. Fuente de la imagen: [10] Esta última gráfica sugiere, además, que el espesor de dieléctrico óptimo es el de 3 mm, tal y como se puede observar en la figura 33: Figura 33. Efecto del espesor de dieléctrico en la componente x de fuerza electrohidrodinámica obtenido con Oracle3D. Fuente de la imagen: [10] Tal y como se expone en [10], la caída de esta fuerza electrohidrodinámica y también de la velocidad inducida del viento iónico para el espesor de dieléctrico de 1 mm, puede ser atribuida a la influencia del dominio del campo eléctrico: con un espesor de 1 mm, los electrodos están muy cerca, y a pesar de que se incrementa el campo eléctrico localmente, esta cercanía reduce el dominio o área de influencia de la fuerza total de Coulomb. Análisis de resultados 64 64 Este efecto, como queda expuesto en la figura 30, no somos capaces de recogerlo con nuestro código en MATLAB. Por otro lado, es interesante graficar la variación de los perfiles de velocidad con el espesor de dieléctrico. En la siguiente figura se presentan los perfiles de velocidad en función del espesor de dieléctrico en x=7 mm y t=36 ms. Figura 34. Perfiles de velocidad con variación de espesor de dieléctrico en x=7 mm y t=36 ms obtenidos con MATLAB En comparación con los obtenidos en [10], en la figura 35: Figura 35. Perfiles de velocidad con variación de espesor de dieléctrico en x=7 mm y t=36 ms obtenidos con Oracle3D. Fuente de la imagen: [10] Se puede observar que, cualitativamente, los perfiles de velocidad son similares que los obtenidos en [10]. De esta forma, a medida que disminuye el espesor de dieléctrico, el campo de velocidades (en una determinada sección x de nuestro actuador) aumenta. 4.1.4 Efecto de la separación (gap) entre los electrodos La separación entre los electrodos es otro parámetro geométrico que participa activamente en el comportamiento del actuador de plasma DBD. Considerando diferentes gaps (0 mm, 1 mm, 3mm y 5mm), se pueden obtener diferentes valores de la velocidad máxima inducida del viento iónico a lo largo del tiempo. Haciendo las modificaciones necesarias en nuestro programa de MATLAB, se muestra en la figura 36 la evolución de la velocidad máxima frente al gap entre los electrodos: Figura 36. Efecto del gap entre los electrodos en las velocidades máximas obtenido con MATLAB Comparando con [10] en la figura 37: Figura 37. Efecto del gap entre los electrodos en las velocidades máximas obtenido con Oracle3D. Fuente de la imagen: [10] En efecto, tal y como se observa en la figura 38, la fuerza componente x de la fuerza volumétrica aumenta a Análisis de resultados 72 72 Figura 49. Efecto de la longitud de Debye en la velocidad máxima inducida. Fuente de la imagen: [42] Se puede observar que los efectos en la velocidad inducida son sustancialmente mayores cuando incrementamos la longitud de Debye. De esta forma, cuando la longitud de Debye es 𝜆𝐷=0.00012 𝑚 las velocidades máximas que se alcanzan son del orden de 1.4 m/s, mientras que, cuando la longitud de Debye aumenta hasta 𝜆𝐷= 0.001 𝑚, las velocidades máximas son aproximadamente de 3.5 m/s. Además, hay un movimiento del fluido desde el electrodo expuesto hacia el electrodo embebido, debido a que la fuerza actúa hacia esa dirección, lo que provoca que el punto de máxima velocidad inducida se desplace aguas abajo. 4.2.2 Efecto de la densidad de carga máxima El siguiente parámetro que analizamos es el del valor de la densidad de carga máxima, que la denotamos como 𝜌𝑐,𝑀𝐴𝑋. Los parámetros que se toman para este estudio son, de acuerdo con [42]: • Espesor de dieléctrico: 0.127 mm • Distancia de separación (gap) entre electrodos: 0.5 mm • Longitud de Debye: 𝜆𝐷=0.17 𝑚𝑚 • Amplitud del voltaje AC aplicado: 𝜙𝑀𝐴𝑋=5 𝑘𝑉 • Longitud del electrodo expuesto: 10 𝑚𝑚 • Longitud del electrodo embebido: 20 𝑚𝑚 Graficando las velocidades máximas inducidas en función de la densidad de carga máxima y comparando con [42], se puede observar que, teniendo en cuenta que la fuerza electrohidrodinámica varía linealmente con el valor de la densidad de carga máxima, pues 𝐹𝑏=𝜌𝑐(−∇𝜙), la velocidad máxima inducida también lo hará aproximadamente; por ello, para 𝜌𝑐,𝑀𝐴𝑋=0.00012𝐶 𝑚3 , se alcanzan velocidades máximas de aproximadamente 0.08 m/s, mientras que para 𝜌𝑐,𝑀𝐴𝑋=0.001𝐶 𝑚3 , es de unos 0.5 m/s. Además, es importante hacer notar que el punto de máxima velocidad inducida permanece prácticamente estático e invariable con 𝜌𝑐,𝑀𝐴𝑋, al contrario que ocurría cuando la longitude de Debye se incrementaba, donde sí que existía un desplazamiento. Figura 50. Campo de velocidades inducido para 𝝆𝒄,𝑴𝑨𝑿=𝟎.𝟎𝟎𝟎𝟏𝟐𝑪 𝒎𝟑 obtenido con MATLAB Figura 51. Campo de velocidades inducido para 𝝆𝒄,𝑴𝑨𝑿=𝟎.𝟎𝟎𝟏𝑪 𝒎𝟑 obtenido con MATLAB Figura 52. Efecto de 𝝆𝒄,𝑴𝑨𝑿 en la velocidad máxima inducida. Fuente de la imagen: [42] Análisis de resultados 74 74 4.3 Análisis del viento iónico producido Una vez que hemos estudiado la influencia de los diferentes parámetros en el funcionamiento de nuestro actuador de plasma DBD, es interesante, por último, estudiar la evolución temporal del campo de velocidades. En concordancia con [18], estudiaremos el campo de velocidades inducidas en diferentes instantes de tiempo real de la simulación: en 0.4 ms, en 1.6 ms, en 4.8 ms, en 9.6 ms, en 19.2 ms, en 30 ms y en 40 ms. Figura 53. Campo de velocidades para t=0.4 ms obtenido con MATLAB Figura 54. Campo de velocidades para t=1.6 ms obtenido con MATLAB Figura 55. Campo de velocidades para t=4.8 ms obtenido con MATLAB Figura 56. Campo de velocidades para t=9.6 ms obtenido con MATLAB Figura 57. Campo de velocidades para t=19.2 ms obtenido con MATLAB Análisis de resultados 76 76 Figura 58. Campo de velocidades para t=30 ms obtenido con MATLAB Figura 59. Campo de velocidades para t=40 ms obtenido con MATLAB Figura 60. Campo de velocidades obtenido con Oracle3D. Fuente de la imagen: [18] Análisis de resultados 78 78 Figura 61. Campo de velocidades vectorial obtenido con Oracle3D. Fuente de la imagen: [18] De acuerdo con las anteriores gráficas, el proceso es el siguiente: un campo eléctrico en el dominio es aplicado mediante una diferencia de potencial o voltaje entre los electrodos. Este campo eléctrico actúa sobre las especies cargadas y un viento iónico es inducido en el dominio. Tras esto, un vórtice debido a este viento iónico aparece, el cual es algo característico de este tipo de descargas, como puede observarse no sólo en [18], sino también en [19] y [36], que se extiende aguas abajo. Como queda constatado, el centro del vórtice se extiende lentamente en la dirección x, al igual que ocurría con el máximo valor de la velocidad inducida al variar el valor de la longitud de debye. Como se puede observer, la dinámica del fluido ocurre en la escala de milisegundos. En 0.4 ms se observa que no se adecua a la perfección con [18], sin embargo, en ese valor de tiempo puede considerarse que el actuador aún se encuentra en el proceso de arranque. Entre 1.6 ms y 4.8 ms se muestra cómo el flujo (jet) se va extendiendo aguas abajo hasta que, en 9.6 ms, se comienza a visualizar la formación del vórtice. Una vez que se ha formado el vórtice, el centro de éste avanza a lo largo del tiempo (aproximadamente en x=12 mm en 19.2 ms, x=15 mm en 30 ms y finalmente x=20 mm en 40 ms), siendo éste rodeado unicamente tanto por fluido en reposo como por el propio viento iónico producido, pues no existe ningún otro flujo externo. 4.4 Coste computacional Resulta interesante, previo a las conclusiones, el estudio del coste computacional necesario para llevar a cabo las simulaciones del problema fluidomecánico. La simulación se lleva a cabo con un HP Pavilion con procesador Intel Core i5-5200U de dos núcleos a 2.20 GHz con 6,00 GB de RAM, siendo el tiempo total de simulación de 200 ms. Los resultados obtenidos quedan expuestos en la siguiente tabla: Número de nodos de Chebyshev en z 15 20 30 40 50 Tiempo total 29.421 s 60.311 s 136.586 s 266.142 s 397.104 s Tabla 4. Tiempos de ejecución del problema fluidomecánico Para las simulaciones recogidas en la tabla 4 se han utilizado 300 nodos en la dirección x en diferencias finitas. Se puede observar cómo, el aumento del número de nodos de Chebyshev en z supone un enorme coste computacional. Por ello, todas las simulaciones realizadas han sido llevadas a cabo con 30-40 nodos en z y entre 200 y 300 nodos en x. Conclusiones, posibles mejoras y líneas de estudio futuras 80 80 5 CONCLUSIONES, POSIBLES MEJORAS Y LÍNEAS DE ESTUDIO FUTURAS 5.1 Conclusiones Una vez llegados a este punto, es el momento de realizar un balance acerca de los objetivos y conclusiones alcanzados. El objetivo principal de este Proyecto de fin de grado ha sido poner de manifiesto que, es posible caracterizar y desarrollar un método numérico sencillo y transparente que permite simular el comportamiento de un actuador DBD con una amplia flexibilidad en cuanto al valor de los parámetros del problema y las condiciones de contorno, hacienda uso de conocimientos básicos de programación y cálculo. La implementación en la herramienta MATLAB de dicho método numérico proporciona un mejor entendimiento del problema y del funcionamiento del dispositivo, creando un entorno abierto a las modificaciones de parámetros, condiciones de contorno e incluso de la forma en que los resultados son expuestos, lejos de los códigos comerciales de tipo ‘’caja negra’’, en los que el usuario no entiende qué se está haciendo en cada paso del proceso. Para llevar a cabo lo descrito, se ha utilizado el modelo de Suzen-Huang para calcular la distribución de densidad de carga y la distribución de potencial eléctrico, para luego incorporarlas a las ecuaciones de Navier Stokes y resolver el campo de velocidades, haciendo uso del método de colocación y la definición de la función de corriente. En el capítulo 1 se han expuesto los fundamentos básicos que rodean a los actuadores de plasma, así como su funcionamiento y su arquitectura, sus ventajas y sus limitaciones. En el capítulo 2 se han desarrollado las ecuaciones que gobiernan la física de los actuadores de plasma, desde el enfoque eléctrico hasta el enfoque fluidomecánico hacienda uso del modelo de Suzen-Huang y se han expuesto las condiciones de contorno que más tarde deberán cumplirse en la simulación del dispositivo. En el capítulo 3 se ha expuesto y desarrollado el método de colocación basado tanto en polinomios de Lagrange en nodos de Chebyshev como en diferencias finitas y la discretización utilizada para resolver tanto el problema eléctrico como el problema hidrodinámico. En el capítulo 4 se ha resuelto numéricamente el problema expuesto en el capítulo 2 y se ha simulado la actuación del dispositivo teniendo en cuenta la influencia de parámetros como el espesor de dieléctrico, el gap entre los electrodos, la amplitud, la frecuencia y la forma de onda utilizada, la densidad de carga máxima o la longitud de Debye, con el fin de validar el funcionamiento de nuestro actuador con investigaciones precedentes. En la última parte de dicho capítulo se analiza el coste computacional de la simulación en función del mallado realizado, con el fin de evidenciar la accesibilidad a los resultados que tiene el Código desarrollado desde cualquier ordenador personal. De esta forma, y en base a los resultados obtenidos, puede concluirse que: • Tanto el espesor de dieléctrico como el gap entre los electrodos son parámetros que juegan un rol importantísimo en el comportamiento de nuestro actuador. Así, se determine que el óptimo para el espesor de dieléctrico era de 3 mm, mientras que no existe limitación para el gap, siendo de 1 mm el óptimo alcanzado en este Proyecto. • A pesar de que que la frecuencia del voltaje aplicado no tiene una influencia significante en el campo de velocidades inducido, la forma de onda y la amplitud del voltaje aplicado sí que protagonizan un papel imprescindible, determinando la forma de onda cuadrada y un voltaje de unos 20-24 kV como los óptimos para maximizar el campo de velocidades inducido. Sin embargo, tal y como explicábamos en la sección 4.1.7 en base a investigaciones precedentes, por sus fluctuaciones y su alto consumo de energía eléctrica, la forma de onda cuadrada es reemplazada por la forma de onda sinusoidal. • También veíamos en la sección 4.2 cómo, tanto la densidad de carga máxima como el valor de la longitud de Debye, influenciaban en gran medida el valor de las velocidades máximas inducidas del viento iónico, tomándose 𝜆𝐷=0.001 𝑚 y 𝜌𝑐,𝑀𝐴𝑋=0.001𝐶 𝑚3 como valores estándar. • Por ultimo, en la sección 4.4 ha quedado constatado cómo, a partir de 30-40 nodos de chebychev en z el coste computacional se dispara, utilizando un total de 40 para todas las simulaciones realizadas. Por tanto, podemos confirmar el cumplimiento de los objetivos iniciales propuestos, premisa fundamental al inicio del trabajo. 5.2 Posibles mejoras y líneas de estudio futuras Una vez que quedan expuestas las conclusiones alcanzadas y los objetivos iniciales propuestos han sido cumplidos y abordados, resulta de interés, de cara a futuros proyectos de investigación enumerar las posibles mejoras que podrían beneficiar, completar y mejorar este en cuestión: • Estudio del efecto del espesor de los electrodos en el campo de velocidades del actuador DBD, que aunque, teniendo en cuenta el estudio realizado en [43], en el que se constata que la influencia del espesor es ínfima en cuanto a la distribución de la densidad de carga y la distribución de potencial eléctrico, resultaría interesante confirmarlo también en el campo de velocidades inducido. • Estudiar alternativas a las condiciones de contorno estáticas impuestas (donde se ha tenido en cuenta que inicialmente el fluido está en reposo) en el problema fluidodinámico, como podría ser la aplicación de un caudal inicial al flujo. • De forma complementaria a las condiciones periódicas estudiadas en [43], que no son de interés en el problema hidrodinámico, sería relevante desarrollar la implementación y simulación de multiactuadores DBD en serie (figura 18), con el fin de estudiar la influencia que tiene cada actuador DBD único sobre el flujo total. • Estudio de la influencia que tendrían el uso de diferentes condiciones de contorno del problema eléctrico sobre el campo de velocidades inducido. • Simulación del actuador DBD sobre dominios que no sean rectangulares, por ejemplo una geometría roma como podría ser un cilindro o sobre una geometría propia de un cuerpo fuselado, como podría ser un perfil NACA, con el fin de estudiar el campo de presiones, la sustentación, la resistencia aerodinámica, el punto de desprendimiento… para extender los resultados alcanzados. • Estudio del efecto de la humedad y la temperatura en la actuación del dispositivo. • Estudio del comportamiento del actuador DBD mediante el uso de otros materiales dieléctricos. • Construcción de un actuador de plasma DBD para validar los resultados de las simulaciones mediante la medida de la fuerza electrohidrodinámica y la velocidad del viento iónico inducido. ANEXO: Programas de matlab 88 88 %%%%%%%% GRÁFICOS %%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% % figure(2);contourf(xmat, zmat, phimat,10); colorbar; % title('Distribución del potencial eléctrico \phi^*') % xlabel('x (mm)') % ylabel('z (mm)') % axis([-30 60 -3 20]) % %%%% Graficamos los electrodos % %1)Electrodo expuesto % eexp_x=[-10 0]; % eexp_y=[0 0]; % %2)Electrodo embebido % eemb_x=[1 21]; % eemb_y=[-3 -3]; % hold on % plot(eexp_x,eexp_y,'r','LineWidth',3) % plot(eemb_x,eemb_y,'b','Linewidth',3) % figure(3) % contourf(xmat, zmat, forcexmat,10);colorbar; % title('Distribución de Fe_x') % xlabel('x (mm)') % ylabel('z (mm)') % %%%% Graficamos los electrodos % %1)Electrodo expuesto % eexp_x=[-10 0]; % eexp_y=[0 0]; % %2)Electrodo embebido % eemb_x=[1 21]; % eemb_y=[-3 -3]; % plot(eexp_x,eexp_y,'r','LineWidth',3) % plot(eemb_x,eemb_y,'b','Linewidth',3) % figure(4) % contourf(xmat, zmat, forcezmat,10);colorbar; % title('Distribución de Fe_z') % xlabel('x (mm)') % ylabel('z (mm)') % %%%% Graficamos los electrodos % %1)Electrodo expuesto % eexp_x=[-10 0]; % eexp_y=[0 0]; % %2)Electrodo embebido % eemb_x=[1 21]; % eemb_y=[-3 -3]; % plot(eexp_x,eexp_y,'r','LineWidth',3) % plot(eemb_x,eemb_y,'b','Linewidth',3) % figure(5) % contourf(xmat, zmat, forcemat,15);colorbar; % title('Distribución de F (módulo)') % xlabel('x (mm)') % ylabel('z (mm)') % axis([-2 8 0 2]) % %%%% Graficamos los electrodos % %1)Electrodo expuesto % eexp_x=[-10 0]; % eexp_y=[0 0]; % %2)Electrodo embebido % eemb_x=[1 21]; % eemb_y=[-3 -3]; % hold on % plot(eexp_x,eexp_y,'r','LineWidth',3) % plot(eemb_x,eemb_y,'g','Linewidth',3) ANEXO 3: Funciones auxiliares function [lp]=dCheby( vp, Np) for m=1:Np, Denp(m)=1; for j=1:Np, if j==m, else Denp(m)=Denp(m)*(vp(m)-vp(j)); end end for k=1:Np, sumk=0; for i=1:Np, if i==m, else prod=1; for l=1:Np, if (l-i)*(l-m)==0, else prod=prod*(vp(k)-vp(l)); end end sumk=sumk+prod; end end lp(k,m)=sumk/Denp(m); end end function[df d2f d3f d4f]=Derivs(x, N) df=sparse(N,N); d2f=sparse(N,N); ANEXO: Programas de matlab 90 90 d3f=sparse(N,N); d4f=sparse(N,N); Nst=9; Nhst=(Nst-1)/2; k=1; [df]=Dmatrix(k, Nst, Nhst, N, x); k=2; [d2f]=Dmatrix(k, Nst, Nhst, N, x); k=3; [d3f]=Dmatrix(k, Nst, Nhst, N, x); k=4; [d4f]=Dmatrix(k, Nst, Nhst, N, x); function D=Dmatrix(k, Nst, Nhst, N, x) D=sparse(N,N); for i=1:N, xbar=x(i); if i <= Nhst, xst(1:Nst)=x(1:Nst); [c]=fdcoeffV(Nst,k,xbar,xst); D(i,1:Nst)=c; end if i > Nhst && i<= N-Nhst, xst(1:Nst)=x((i-Nhst):(i+Nhst)); [c]=fdcoeffV(Nst,k,xbar,xst); D(i,(i-Nhst):(i+Nhst))=c; end if i > N-Nhst, xst(1:Nst)=x((N-Nst+1):N); [c]=fdcoeffV(Nst,k,xbar,xst); D(i,(N-Nst+1):N)=c; end end function c = fdcoeffV(n,k,xbar,x) A = ones(n,n); xrow = (x(:)-xbar); % displacements as a row vector. hmax=abs(x(n)-x(1)); for i=2:n A(i,:) = ((xrow/hmax) .^(i-1)) ./ factorial(i-1); end b = zeros(n,1); % b is right hand side, b(k+1) = 1/hmax^k; % so k’th derivative term remains% Ainv=inv(A); c = (A\b)'; % solve system for coefficients REFERENCIAS [1] BREWER, G. R.; CURRIE, M. R.; KNECHTLI, R. C. Ionic and plasma propulsion for space vehicles. Proceedings of the IRE, 1961, vol. 49, no 12, p. 1789-1821. [2] JANG, Jaesung; LEE, Seung S. Theoretical and experimental study of MHD (magnetohydrodynamic) micropump. Sensors and Actuators A: Physical, 2000, vol. 80, no 1, p. 84-89. [3] MOREAU, Eric. Airflow control by non-thermal plasma actuators. Journal of physics D: applied physics, 2007, vol. 40, no 3, p. 605. [4] MOREAU, E., et al. Surface dielectric barrier discharge plasma actuators. ERCOFTAC Bulletin, 2013, vol. 94, p. 5-10 [5] JUKES, Tim, et al. Turbulent boundary-layer control for drag reduction using surface plasma. En 2nd AIAA flow control conference. 2004. p. 2216. [6] HE, Chuan; CORKE, Thomas C.; PATEL, Mehul P. Plasma flaps and slats: an application of weakly ionized plasma actuators. Journal of Aircraft, 2009, vol. 46, no 3, p. 864-873. [7] ERFANI, Rasool, et al. Development of DBD plasma actuators: the double encapsulated electrode. Acta Astronautica, 2015, vol. 109, p. 132-143. [8] NERETTI, Gabriele. Active flow control by using plasma actuators. Recent Progress in Some Aircraft Technologies, 2016, p. 57-76. [9] AKANSU, Y. E.; KARAKAYA, F.; ŞANLISOY, A. Active control of flow around NACA 0015 airfoil by using DBD plasma actuator. En EPJ web of conferences. EDP Sciences, 2013. p. 01008. [10] SETH, U. K., et al. Parametric study of a DBD plasma actuation based on the Suzen-Huang model. Journal of Electrostatics, 2018, vol. 93, p. 1-9. [11] SUNG, Y., et al. Aerodynamic modification of flow over bluff objects by plasma actuation. Experiments in Fluids, 2006, vol. 41, no 3, p. 479-486. [12] RODRIGUES, F. F.; PASCOA, J. C.; TRANCOSSI, Michele. Experimental analysis of alternative dielectric materials for DBD plasma actuators. En ASME 2018 International Mechanical Engineering Congress and Exposition. American Society of Mechanical Engineers Digital Collection, 2018. [13] MCKARTHY, Kieran J., et al. Desarrollo de diagnósticos para plasmas: Catalizador de nueva física. Revista Española de Física, 27, 25-39., 2013. [14] DONKO, Z.; HARTMANN, P.; KALMAN, G. J. Strongly coupled plasma liquids. arXiv preprint arXiv:0710.5229, 2007. [15] TANARRO, Isabel. PLASMA, El Cuarto Estado de la Materia. 2011. [16] GOLDSTON, Robert J.; RUTHERFORD, Paul Harding. Introduction to plasma physics. CRC Press, 1995. [17] ORLOV, Dmitriy M. Modelling and simulation of single dielectric barrier discharge plasma actuators. 2006. [18] U.K Seth. Application to the numerical simulation of various typical Electro-hydro-dynamic flows. PhD tesis. University of Poitiers, 2018 [19] BENARD, Nicolas; MOREAU, Eric. Electrical and mechanical characteristics of surface AC dielectric barrier discharge plasma actuators applied to airflow control. Experiments in Fluids, 2014, vol. 55, no 11, p. 1846. [20] https://www.areaciencias.com/fisica/plasma.html [21] http://www.gaiaciencia.com/2017/09/el-plasma-ese-gran-desconocido/ [22] http://saberesyciencias.com.mx/2015/01/07/como-se-produce-la-luz/ Referencias 92 92 [23] http://ocw.uc3m.es/ciencia-e-oin/quimica-de-los-materiales/Material-de-clase/tema-4.-solidos-liquidos-y- disoluciones-ii/skinless_view [24] https://www.plasmatreat.es/tecnologia_del_plasma/tecnologia_openair_plasma.html [25] https://www.nationalgeographic.es/espacio/el-sol-0 [26] https://concepto.de/estado-plasmatico/ [27] https://www.meteorologiaenred.com [28] https://www.elpais.com.uy [29] http://hyperphysics.phy-astr.gsu.edu/hbasees/electric/dielec.html [30] SUZEN, Yildirim, et al. Numerical simulations of plasma based flow control applications. En 35th AIAA Fluid Dynamics Conference and Exhibit. 2005. p. 4633. [31] BELSON, Brandt, et al. Comparison of plasma actuators in simulations and experiments for control of bypass transition. En 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition. 2012. p. 1141. [32] SATO, Makoto, et al. LES of Separated-flow Controlled by DBD Plasma Actuator around NACA 0015 over Reynolds Number Range of 10^ 4-10^ 6. En 53rd AIAA Aerospace Sciences Meeting. 2015. p. 0308. [33] IBRAHIM, I. H.; SKOTE, M. Simulations of the linear plasma synthetic jet actuator utilizing a modified Suzen-Huang model. Physics of Fluids, 2012, vol. 24, no 11, p. 113602. [34] SKOTE, Martin; IBRAHIM, Imran Halimi. Utilizing the L-PSJA for controlling cylindrical wake flow. International Journal of Numerical Methods for Heat & Fluid Flow, 2016, vol. 26, no 5, p. 1593-1616. [35] MAHFOZE, O.; LAIZET, S. Skin-friction drag reduction in a channel flow with streamwise-aligned plasma actuators. International Journal of Heat and Fluid Flow, 2017, vol. 66, p. 83-94. [36] BRAUNER, Timothy, et al. Modelling of dielectric barrier discharge plasma actuators for direct numerical simulations. En 8th AIAA Flow Control Conference. 2016. p. 3774. [37] MORALES, Carlos. Convección natural de Rayleigh-Bénard en recintos rectangulares (Proyecto Mecánica de Fluidos Avanzada). Trabajo de Fin de Grado. Universidad de Sevilla, Sevilla, 2016. [38] RUIZ, Pablo. El método de colocación para el problema de convección de Rayleigh-Bénard, Proyecto de fin de Carrera, Universidad de Sevilla, Sevilla, 2016 [39] https://www.pybonacci.org/2013/08/15/ajuste-e-interpolacion-unidimensionales-basicos-en-python-con- scipy/ [40] BENARD, N.; MOREAU, E. Role of the electric waveform supplying a dielectric barrier discharge plasma actuator. Applied Physics Letters, 2012, vol. 100, no 19, p. 193503. [41] SATO, Shintaro; OHNISHI, Naofumi. Influence of voltage waveform on electrohydrodynamic force in a dielectric-barrier-discharge plasma actuator. En 55th AIAA Aerospace Sciences Meeting. 2017. p. 1804. [42] HALIMI BIN IBRAHIM, Imran; SKOTE, Martin. Effects of the scalar parameters in the Suzen-Huang model on plasma actuator characteristics. International Journal of Numerical Methods for Heat & Fluid Flow, 2013, vol. 23, no 6, p. 1076-1103. [43] DOMIGUEZ, Miguel A. Modelado y estudio de un actuador SDBD para el control de flujo. Caracterización de la fuerza hidrodinámica generada, Trabajo de fin de grado, Universidad de Sevilla, Sevilla, 2019. [44] CORKE, Thomas C.; ENLOE, C. Lon; WILKINSON, Stephen P. Dielectric barrier discharge plasma actuators for flow control. Annual review of fluid mechanics, 2010, vol. 42, p. 505-529. [45] LATEN, Jeff B.; LEBEAU, Raymond P. Improving the Performance of a Plasma Actuator Model for DBD and Multi-Encapsulated Electrode Actuators. En 55th AIAA Aerospace Sciences Meeting. 2017. p. 1808. [46] PÉREZ-SABORID, M. Notas sobre el cálculo general y eficiente de matrices de derivadas mediante los polinomios de Lagrange y sus aplicaciones, ETSI. Sin publicar, 2019. GLOSARIO DBD: Dielectric Barrier Discharge AC: Corriente alterna SH: Suzen-Huang NLP: Término no lineal Glosario 94 94